      program dynamics
      parameter(maxdim=999)
      real c(maxdim),xmass(maxdim),tau,
     $V(maxdim,3),XDELTA(maxdim),temp
      character*2 sym(maxdim)
      character*2 choice
      COMMON /CSEED/ ISEED
      write(*,*)' ISEED =?'
      read(*,*) ISEED 

       open(unit=1,file='xyz')
       read(1,*) nat
       write(*,*)nat
     
      do i=1,nat
       read(1,*) sym(i),c(3*i-2),c(3*i-1),c(3*i)
       write(*,*) sym(i),c(3*i-2),c(3*i-1),c(3*i)
      enddo
     
      close(1)
     
      do i=1,nat
       if(sym(i).eq.'au'.or.sym(i).eq.'Au'.or.sym(i).eq.'AU') then
        xmass(i)=359047.964
        elseif (sym(i).eq.'c'.or.sym(i).eq.'C') then
        xmass(i)=21892.7109375
        elseif (sym(i).eq.'o'.or.sym(i).eq.'O') then
        xmass(i)=29163.55664062
        elseif(sym(i).eq.'ag'.or.sym(i).eq.'Ag'.or.sym(i).eq.'AG') then
        xmass(i)=196631.6863
        elseif (sym(i).eq.'h'.or.sym(i).eq.'H') then
        xmass(i)=1837.374491
        elseif (sym(i).eq.'N'.or.sym(i).eq.'n') then
        xmass(i)=25532.68181
        elseif(sym(i).eq.'li'.or.sym(i).eq.'Li'.or.sym(i).eq.'LI') then
        xmass(i)=12653.03019
        elseif(sym(i).eq.'na'.or.sym(i).eq.'Na'.or.sym(i).eq.'NA') then
        xmass(i)=41907.77622
        elseif(sym(i).eq.'mg'.or.sym(i).eq.'Mg'.or.sym(i).eq.'MG') then
        xmass(i)=44305.40221
        elseif(sym(i).eq.'sc'.or.sym(i).eq.'Sc'.or.sym(i).eq.'SC') then
        xmass(i)=81949.5903
        elseif(sym(i).eq.'ti'.or.sym(i).eq.'Ti'.or.sym(i).eq.'TI') then
        xmass(i)=87256.36218
        elseif(sym(i).eq.'br'.or.sym(i).eq.'Br'.or.sym(i).eq.'BR') then
        xmass(i)=145656.225
        elseif (sym(i).eq.'V'.or.sym(i).eq.'v') then
        xmass(i)=92860.66728
        elseif(sym(i).eq.'nb'.or.sym(i).eq.'Nb'.or.sym(i).eq.'NB') then
        xmass(i)=169357.9289
        elseif(sym(i).eq.'ta'.or.sym(i).eq.'Ta'.or.sym(i).eq.'TA') then
        xmass(i)=329847.7738
        elseif(sym(i).eq.'si'.or.sym(i).eq.'Si'.or.sym(i).eq.'SI') then
        xmass(i)=51196.77561
        elseif(sym(i).eq.'f' .or.sym(i).eq.'F') then
        xmass(i)=34628.81288
        elseif (sym(i).eq.'s'.or.sym(i).eq.'S') then
        xmass(i)=58447.413      
        elseif(sym(i).eq.'cl'.or.sym(i).eq.'Cl'.or.sym(i).eq.'CL') then
        xmass(i)=64629.86192
        elseif(sym(i).eq.'Mo'.or.sym(i).eq.'mo'.or.sym(i).eq.'MO') then
        xmass(i)=174896.03
        elseif(sym(i).eq.'zr'.or.sym(i).eq.'Zr'.or.sym(i).eq.'ZR') then
        xmass(i)=166287.32812500
        elseif(sym(i).eq.'Q'.or.sym(i).eq.'q') then
        xmass(i)=0.00
       endif
      enddo
    
      write(*,111) 
111   format("MD with constant energy (e) or constant temperature (t) ?"
     $)
      read(*,*) choice
    
 100  format(a,1x,f24.8,1x,f16.12,1x,f16.12,1x,f16.12)
 
C     Input for MD with constant energy
 
      if(choice.eq.'e'.or.choice.eq.'E') then
         
c     write(*,*) 'Enter time-step in fs'
c     read(*,*) timestep
     
c     write(*,*) 'Enter number of steps'
c     read(*,*) nstep
     
      write(*,*) 'Enter initial kinetic energy in eV'
      read(*,*) E

      open(unit=1,file='dynamics.in')

      write(1,*) 'E'
      write(1,*) nat
c     write(1,*) timestep
c     write(1,*) nstep
      write(1,*) E
       do i=1,nat
        write(1,100) sym(i),xmass(i),1.88973*c(3*i-2),
     $1.88973*c(3*i-1),1.88973*c(3*i)
        enddo
      
      
C     Calculation of the initial velocities  

      energy=E/27.2114
      
      DO I=1,NAT
      XDELTA(I)=SQRT(2.0*energy/(NAT*xmass(I)))
      ENDDO

      DO I=1,NAT
      CALL GRNF(X,Y)
      V(I,1)=XDELTA(I)*X
      CALL GRNF(X,Y)
      V(I,2)=XDELTA(I)*X
      CALL GRNF(X,Y)
      V(I,3)=XDELTA(I)*X
      ENDDO
      DO I=1,NAT

      WRITE(1,*) V(I,1),V(I,2),V(I,3)
      ENDDO

C     Input for MD with constant temperature
     
      elseif(choice.eq.'t'.or.choice.eq.'T') then
         
c     write(*,*) 'Enter time-step in fs'
c     read(*,*) timestep
      
c     write(*,*) 'Enter number of steps'
c     read(*,*) nstep
      
      write(*,*) 'Enter temperature in K for MD'
      read(*,*) temp

c     write(*,112)  
c12   format("Parameter tau defining strength of coupling with heat bath
c    $")
c     read(*,*) tau

      open(unit=1,file='dynamics.in')

      write(1,*) 'T'
      write(1,*) nat
c     write(1,*) timestep
c     write(1,*) nstep
      write(1,*) temp
c     write(1,*) tau
       do i=1,nat
        write(1,100) sym(i),xmass(i),1.88973*c(3*i-2),
     $1.88973*c(3*i-1),1.88973*c(3*i)
       enddo
      
C     Calculation of the initial velocities  

C     E=(3N-6)*K_boltz*T/2
      energy=(temp*((3*nat)-6))/(23209*27.2114)
      DO I=1,NAT
      XDELTA(I)=SQRT(2.0*energy/(NAT*xmass(I)))
      ENDDO

      DO I=1,NAT
      CALL GRNF(X,Y)
      V(I,1)=XDELTA(I)*X
      CALL GRNF(X,Y)
      V(I,2)=XDELTA(I)*X
      CALL GRNF(X,Y)
      V(I,3)=XDELTA(I)*X
      ENDDO
      DO I=1,NAT

      WRITE(1,*) V(I,1),V(I,2),V(I,3)
      ENDDO


   
       
      endif
      
      end 

C*********************************************
C*    LINEAR CONGRUENTIAL RANDOM NUMBER      *
C*               GENERATOR                   *
C*********************************************
      FUNCTION RANF()
      DATA IA/16807/,IC/2147483647/,IQ/127773/,IR/2836/
      COMMON /CSEED/ ISEED
      IH = ISEED/IQ
            IL = MOD(ISEED,IQ)
            IT = IA*IL-IR*IH
            IF(IT.GT.0) THEN
                 ISEED = IT
            ELSE
                 ISEED = IC+IT
            END IF
            RANF = ISEED/FLOAT(IC)
      RETURN
      END
C********************************************
C*    GAUSSIAN RANDOM NUMBER GENERATOR      *
C********************************************
      SUBROUTINE GRNF(X,Y)
      COMMON /CSEED/ ISEED
            PI = 4.0*ATAN(1.0)
            R1 = -ALOG(1.0-RANF())
            R2 = 2.0*PI*RANF()
            R1 = SQRT(2.0*R1)
            X  = R1*COS(R2)
            Y  = R1*SIN(R2)
      RETURN
      END
 
      





