PROGRAM MD_2D_Diff
*******************************************************************
*     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

      REAL R(97)
      EXTERNAL RAN


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

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

      AMASS=26.981538
      DEL=1.

      NSTEP=50

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 gewuenschte Temperatur in Kelvin ein:'

      READ (*,*) TEMPI

      IF (TEMPI.EQ.0.) THEN

        WRITE (*,*) 'Temperatur muss groesser als Null sein.'
        WRITE (*,*) 'Neue Eingabe:'
        READ (*,*) TEMPI

      ENDIF


      WRITE (*,*) 'Gebe einen Anfangswert ISEED (Natuerliche Zahl):'

      READ (*,*) ISEED

      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 (*,*) 'xmgrace Trajektorie_Al_Diff.dat'
      WRITE (*,*) 'anschauen'


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

      EKIN=2.*TEMPI*8.6173431E-5

      V=SQRT(2.*EKIN/AMASS)

      X0=X0*A
      Y0=Y0*A

      ANG=360.*RAN(ISEED,IX1,IX2,IX3,R)

      WRITE (*,*) 'Anfangsrichtung: ',ANG

      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

      ISTEP=INT(N/NSTEP)

      WRITE (*,*) 'Anzahl der Schritte: ',N, ISTEP

      DO 10, I=1,ISTEP+1

         NI=NSTEP*(I-1)
         CALL VERLET (NSTEP,X0,Y0,VX,VY,XF,YF,VXF,VYF,DEL,AMASS,NI)

         EKINF=.5*AMASS*(VXF**2+VYF**2)
C        2 Freiheitsgrade => Ekin = 2 * (1/2) * k_B T                    
         TEMP=EKIN/8.6173431E-5

         X0=XF
         Y0=YF
         ANG=(180.*ATAN2(VYF,VXF)/PI)
         ANG=ANG+15.0*(2.*RAN(ISEED,IX1,IX2,IX3,R)-1.)

         EKIN=TEMPI*8.6173431E-5
         V=SQRT(2.*EKIN/AMASS)
         VX=V*COS(ANG*PI/180.)
         VY=V*SIN(ANG*PI/180.)

 10    CONTINUE

      END


      SUBROUTINE VERLET (N,X0,Y0,VX0,VY0,XF,YF,VXF,VYF,DEL,AMASS,NI)
*********************************************************************
*     VERLET Algorithmus
*********************************************************************
      
      REAL R(97)
      EXTERNAL RAN


      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,*) (NI+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

*
*
      REAL FUNCTION RAN(IDUM,IX1,IX2,IX3,R)
**********************************************************************
*     Returns a random deviate between 0.0 and 1.0. Set IDUM to
*     any negative value to initialize or reinitialize the sequence.
*     Numerical Recipes p.196
**********************************************************************
      REAL R(97)
      INTEGER  IX1,IX2,IX3,IC1,IC2,IC3,IA1,IA2,IA3,IFF
      PARAMETER (M1=259200,IA1=7141,IC1=54733,RM1=1./M1)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)

      DATA IFF /0/

      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF=1
         IX1=MOD(IC1-IDUM,M1)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX2=MOD(IX1,M2)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX3=MOD(IX1,M3)
         DO 11, J=1,97
           IX1=MOD(IA1*IX1+IC1,M1)
           IX2=MOD(IA2*IX2+IC2,M2)
           R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
 11      CONTINUE
         IDUM=1

       ENDIF

       IX1=MOD(IA1*IX1+IC1,M1)
       IX2=MOD(IA2*IX2+IC2,M2)
       IX3=MOD(IA3*IX3+IC3,M3)
       
       J=1+(97*IX3)/M3
       
*       IF (J.GT.97.OR.J.LT.1) PAUSE

       RAN=R(J)  


       R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

       RETURN
       END
*
*