PROGRAM MD_2D
*******************************************************************
*     Programm zur Durchfuehrung von Molekulardynamiksimuationen
*     mit dem Verlet Algorithmus. Propagation eines Metallteilchens
*     mit der Masse von Al (26,981538 amu) in einem periodischen
*     zweidimensionalen Potential
*     Zeitschritt 1fs, Gitterkonstante A= 3 Angstrom
******************************************************************

      COMMON /DATA/ A,PI,POT0

      OPEN (10,FILE="Trajektorie_Al.dat",
     &            FORM="FORMATTED",STATUS="UNKNOWN")
      OPEN (25,FILE="Trajektorie_Al_all.dat",
     &      FORM="FORMATTED",STATUS="UNKNOWN")

      POT0=1.
      A=3.
      PI=ACOS(-1.)

      AMASS=26.981538
      DEL=1.

C   X0 x-Koordinate des Anfangspunktes in der Einheitszelle in Einheiten von a
C   Y0 y-Koordinate des Anfangspunktes in der Einheitszelle in Einheiten von a
C   EKIN  kinetische Energie am Anfang in eV
C   ANG Winkel der Anfangsrichtung
C   AMASS Masse in amu des Teilchens
C   DEL  Zeitschritt in fs
C   N Anzahl der Zeitschritte

      Write (*,*) "Molekulardynamiksimulationen"

CAG      WRITE (*,*) "Gebe die Anfangskoordinaten X0 und Y0 in der"
CAG      WRITE (*,*) "Einheitszelle an (0<X0,Y0<1):"
CAG      READ (*,*) X0,Y0

      X0=0.5
      y0=0.5

      WRITE (*,*) 'Gebe die kinetische Energie in eV an (0<Ekin<1):'

      READ (*,*) EKIN

      WRITE (*,*) 'Gebe die Anfangsrichtung an (0<Ang<360):'

      READ (*,*) ANG

      WRITE (*,*) 'Gebe die Anzahl der Zeitschritte an (0< N):'

      READ (*,*) N

      WRITE (*,*) 'Nach Vollendung des Programmes kan man sich die'
      WRITE (*,*) 'berechnete Trajektorie mit'
      WRITE (*,*) 'Matlab (Windows)'
      WRITE (*,*) 'anschauen'


      DEL=DEL*1.E-15
      AMASS=AMASS*1.03643533E-28

      V=SQRT(2.*EKIN/AMASS)

      X0=X0*A
      Y0=Y0*A

      VX=V*COS(ANG*PI/180.)
      VY=V*SIN(ANG*PI/180.)

      VPOT=POT(X0,Y0)
      EKIN=.5*AMASS*(VX**2+VY**2)
      ETOT=EKIN+VPOT

      WRITE (10,*) X0,Y0
      WRITE (25,*) 0,X0,Y0,EKIN,VPOT,ETOT


      CALL VERLET (N,X0,Y0,VX,VY,XF,YF,VXF,VYF,DEL,AMASS)

      END


      SUBROUTINE VERLET (N,X0,Y0,VX0,VY0,XF,YF,VXF,VYF,DEL,AMASS)
*********************************************************************
*     VERLET Algorithmus
*********************************************************************
      


      CALL ABLEITUNG(X0,Y0,FX,FY)
 
      X1=X0+DEL*VX0+(.5*FX*DEL**2)/AMASS
      Y1=Y0+DEL*VY0+(.5*FY*DEL**2)/AMASS

      DO 100, I=1,N
 
         VPOT=POT(X1,Y1)
         CALL ABLEITUNG(X1,Y1,FX,FY)

         X2=2.*X1-X0+FX*DEL**2/AMASS
         Y2=2.*Y1-Y0+FY*DEL**2/AMASS
         VX1=(X2-X0)/(2.*DEL)
         VY1=(Y2-Y0)/(2.*DEL)

         EKIN=.5*AMASS*(VX1**2+VY1**2)
         ETOT=EKIN+VPOT
     
         WRITE (10,*) X1,Y1
         WRITE (25,*) I,X1,Y1,EKIN,VPOT,ETOT

         X0=X1
         X1=X2
         Y0=Y1
         Y1=Y2

 100  CONTINUE


      XF=X2
      YF=Y2
      VXF=VX1         
      VYF=VY1         

      END

  
      SUBROUTINE ABLEITUNG(X,Y,FX,FY)
****************************************************************
*     Kraft = -Gradient des Potentials
***************************************************************
      COMMON /DATA/ A,PI,POT0

      ARG=2.*PI/A

      POT=.25*POT0*(2.+COS(ARG*X)+COS(ARG*Y))

     
      FX=.25*POT0*ARG*SIN(ARG*X)
      FY=.25*POT0*ARG*SIN(ARG*Y)
      
      END
*
*
      REAL FUNCTION POT(X,Y)
  
      COMMON /DATA/ A,PI,POT0
    
      ARGX=2.*PI*X/A
      ARGY=2.*PI*Y/A

      POT=.25*POT0*(2.+COS(ARGX)+COS(ARGY))
 
      END