md1.f90

分子动力学程序

if (VelAcc) then
   print '(a)','#          (positions, velocities, accelerations read from file)'
else
   print '(a)','#                                (only positions read from file)'
endif

print '(a,i8,a,f7.4,a,f12.4)', &
   '# Number of steps:',Nsteps,', time step:',deltat, &
   ', total time:',Nsteps*deltat
print '(a,i6)','# Number of atoms:',N
print '(a,3f12.6,a,f15.3)', &
   '# Box size:',BoxSize,', Volume:',volume
if (ChangeRho) then
   print '(a,f12.6,a)','# Density:',density,' (changed)'
else
   print '(a,f12.6,a)','# Density:',density,' (unchanged)'
endif
if (ConstantT) then
   print '(a,f12.6)','# Constant T run with T =',TRequested
else
   print '(a)','# Free evolution run.'
endif

!
!  Now print headers of columns
!


print '(a,/,a,/,a)', '#', &
'#  Step   Temperature     Kinetic      Potential   Total Energy    Pressure',&
'# -----  ------------  ------------  ------------  ------------  ------------'

end subroutine Initial_Printout

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  TIME EVOLUTION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     ROUTINES

subroutine Evolve_Sample
!
!  This is the main subroutine, controlling the time evolution of the system.
!
use Particles
use Simulation_Control
use Statistics
implicit none
integer step
double precision :: ene_kin_aver,ene_pot_aver,ene_tot_aver,temperature,pressure
double precision :: chi



!
!  We need to have the initial temperature ready in case we are going
!  at constant T:
!
call Compute_Temperature(ene_kin_aver,temperature) 
!
!  "Velocity Verlet" integrator (see e.g. Allen and Tildesley book, p. 81).
!  Simple velocity scaling (done on velocities at the previous step)
!  applied when ConstantT is enabled.
!



time: do step=1,Nsteps
   call Refold_Positions

   pos = pos + deltat*vel + 0.5d0*(deltat**2)*acc      ! r(t+dt)
   if (ConstantT .and. (temperature > 0) ) then	  !  veloc rescale for const T
      call Compute_Temperature(ene_kin_aver,temperature) ! T(t)
      chi = sqrt( Trequested / temperature )
      vel = chi*vel + 0.5d0*deltat*acc                 ! v(t+dt/2)
   else                          !  regular constant E dynamics
      vel = vel + 0.5d0*deltat*acc                     ! v(t+dt/2)
   endif
  
 call Compute_Forces                                 ! a(t+dt),ene_pot,virial

  
   vel = vel + 0.5d0*deltat*acc                        ! v(t+dt)
   call Compute_Temperature(ene_kin_aver,temperature)  ! at t+dt, also ene_kin
   ene_pot_aver = sum( ene_pot ) / N
   !!
   !!add by myself
   print*,'ene_pot_aver=',ene_pot_aver
   !!
   !!

   ene_tot_aver = ene_kin_aver + ene_pot_aver
!
!  For the pressure calculation, see the Allen and Tildesley book, section 2.4
!
   pressure = density*temperature + virial/volume
   print '(1x,i6,5f14.6)',step,temperature, &
                          ene_kin_aver,ene_pot_aver,ene_tot_aver,pressure
!     accumulate statistics:
   temperature_sum = temperature_sum + temperature
   ene_kin_sum     = ene_kin_sum     + ene_kin_aver
   ene_pot_sum     = ene_pot_sum     + ene_pot_aver
   pressure_sum    = pressure_sum    + pressure
enddo time
end subroutine Evolve_Sample


subroutine Refold_Positions
!
!  Particles that left the box are refolded back into the box by
!  periodic boundary conditions 
!
use Particles
implicit none
where ( pos >  0.5d0 ) pos = pos - 1.d0
where ( pos < -0.5d0 ) pos = pos + 1.d0
end subroutine Refold_Positions


subroutine Compute_Forces
!
!  Compute forces on atoms from the positions, using the Lennard-Jones 
!  potential.
!  Note double nested loop, giving O(N^2) time: this is a SLOW ROUTINE,
!  unsuitable for large systems.
!
use Particles
use Potential
implicit none
double precision, dimension(DIM) :: Sij,Rij
double precision :: Rsqij,phi,dphi
double precision :: rm2,rm6,rm12
integer :: i,j
!
!  Reset to zero potential energies, forces, virial term
!
ene_pot = 0.d0
acc = 0.d0
virial = 0.d0
!
!  Loop over all pairs of particles
!
do i = 1,N-1                                      ! looping an all pairs
   do j = i+1,N
      Sij = pos(:,i) - pos(:,j)                   ! distance vector between i j
 
	  
	  where ( abs(Sij) > 0.5d0 )                  ! (in box scaled units)
         Sij = Sij - sign(1.d0,Sij)               ! periodic boundary conditions
      end where                                   ! applied where needed.
      Rij = BoxSize*Sij                           ! go to real space units
      Rsqij = dot_product(Rij,Rij)                ! compute square distance
      if ( Rsqij < Rcutoff**2 ) then              ! particles are interacting?
                                                  ! compute Lennard-Jones potntl
         rm2 = 1.d0/Rsqij                         !  1/r^2
         rm6 = rm2**3                             !  1/r^6
         rm12 = rm6**2                            !  1/r^12
         phi  = 4.d0 * ( rm12 - rm6 ) - phicutoff !  4[1/r^12 - 1/r^6] - phi(Rc)
                                      !  The following is dphi = -(1/r)(dV/dr)
         dphi = 24.d0*rm2*( 2.d0*rm12 - rm6 )     !  24[2/r^14 - 1/r^8]
         ene_pot(i) = ene_pot(i) + 0.5d0*phi      ! accumulate energy
         ene_pot(j) = ene_pot(j) + 0.5d0*phi      ! (i and j share it)
         virial = virial - dphi*Rsqij             ! accumul. virial=sum r(dV/dr)
         acc(:,i) = acc(:,i) + dphi*Sij           ! accumulate forces
         acc(:,j) = acc(:,j) - dphi*Sij           !    (Fji = -Fij)
      endif

   enddo
   	!!
   !!add by myself
   print*,'Rsq',i,j,'=',Rsqij
   print*,'Rcutoff**2=',Rcutoff**2
   !!
   !!

enddo
virial = - virial/DIM                             ! definition of virial term
end subroutine Compute_Forces


subroutine Compute_Temperature(ene_kin_aver,temperature)
!
!  Starting from the velocities currently stored in vel, it updates
!  the kinetic energy array ene_kin, and computes and returns the
!  averaged (on particles) kinetic energy ene_kin_aver and the
!  instantaneous temperature.
!
use Particles
implicit none
double precision :: ene_kin_aver, temperature
double precision, dimension(DIM) :: real_vel
integer :: i
do i=1,N
   real_vel = BoxSize*vel(:,i)                        ! real space velocity of i
   ene_kin(i) = 0.5d0*dot_product(real_vel,real_vel)  ! kin en of each atom 
enddo
ene_kin_aver = sum( ene_kin ) / N
temperature = 2.d0*ene_kin_aver/DIM
end subroutine Compute_Temperature


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  TERMINATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    ROUTINES

subroutine Terminate
!
!  Termination procedure (called once at the end, after time evolution
!  and before program termination)
!
use Particles
implicit none
!
!  Print a line with averages, etc
!
call Print_Statistics
!
!  Write the final sample (positions, velocities, accelerations) on file
!
call Write_Sample
!
!  Deallocate all the dynamical arrays to clean memory ('ecological' practice)
!
deallocate( pos )
deallocate( vel )
deallocate( acc )
deallocate( ene_pot )
deallocate( ene_kin )
end subroutine Terminate


subroutine Print_Statistics
!
!  Print the mean value, averaged during the run, of the statistical 
!  quantities which have been accumulated
!  
use Simulation_Control
use Statistics
implicit none
if ( Nsteps <= 0 ) return               ! protection against '0 steps run'
print '(a,/,a,5f14.6)','#','# Means', &
                     temperature_sum / Nsteps , &
                     ene_kin_sum / Nsteps , &
                     ene_pot_sum / Nsteps , &
                     ( ene_kin_sum + ene_pot_sum ) / Nsteps , &
                     pressure_sum / Nsteps         
end subroutine Print_Statistics


subroutine Write_Sample
!
!  Writes the final sample to file unit 2
!
use Particles
use Simulation_Control
implicit none
double precision, dimension(DIM) :: PosAtomReal,VelAtomReal,AccAtomReal
integer :: i,lis

lis = len_trim(SampOut)
open(unit=2,file=SampOut(1:lis),status='unknown', &
     action='write',err=700)
VelAcc = .TRUE.  ! we are going to write velocities and accelerations in SampOut
!
!  The '%' signals to ignore this line to the BallRoom program producing
!  an image from the simulation
!
write(2,'(A1,L2,I7,3E23.15)',err=900) '%',VelAcc,N,BoxSize
!
!  Multiply coordinates (which are scaled by the box size) by BoxSize in
!  order to have them in real, unscaled units, then write them in the
!  file (one line per atom).
!
do i=1,N
   PosAtomReal = pos(:,i)*BoxSize
   write(2,'(1X,3E23.15)',err=900) PosAtomReal
enddo
!
!  Do the same for velocities and accelerations.
!
do i=1,N
   VelAtomReal = vel(:,i)*BoxSize
   write(2,'(A1,3E23.15)',err=900) '%',VelAtomReal
enddo
do i=1,N
   AccAtomReal = acc(:,i)*BoxSize
   write(2,'(A1,3E23.15)',err=900) '%',AccAtomReal
enddo
!
!  file written successfully if we get to this point
!
close(unit=2)
return
!
!  handling of various kinds of errors
!
700 continue
   print*,'Write_Sample: WARNING: cannot open ',SampOut,' for write'
   return
900 continue
   print*,'Write_Sample: WARNING: error when writing ',SampOut
   close(unit=2)
   return
end subroutine Write_Sample