I am writing a code to perform hybrid monte carlo molecular dynamics. To do this, I need to have a code to initialize the velocities of all particles according to a maxwell distribution. The code is as follows
subroutine thermalize implicit none double precision :: MassInv double precision :: RATIO integer :: i !--- temp is temperature set point value STDtemp = sqrt(temp) do i = 1 ,np MassInv = 1.0d0/(mass(position(i)%type)) MassInv = sqrt(MassInv*kboltz*kcal2amu) !---Gaussian() is a functions that returns gaussian random variable n(0,1) v(i)%x = Gaussian()*STDtemp*MassInv*ratio v(i)%y = Gaussian()*STDtemp*MassInv*ratio v(i)%z = Gaussian()*STDtemp*MassInv*ratio enddo end subroutine thermalize
To measure the temperature
T = 2.0*KE/(kboltz*(3*np - NCONS)
where KE is the kinetic energy, np is the number of atoms, kbotlz is boltzmanns constant. I was trying to run some tests with temp=380, and found that my thermalize function was giving me initial kinetic energies that translated to temperatures around T=550 K. When I turned shake off, and hence ncons = 0, the temperature produced by thermalize was 379. So I think I need to modify the thermalize command to account for the constrained DOF's
To try to account for the shake, the ratio of the KE with shake and without shake present for a fixed temperature is
1/2kT*(3np-ncons) = 1/2kT*(3np)
dividing and rearranging
KE_shake/KE_noshake = (3np-ncons)/np
Thus, I am thinking about multiplying by the sqrt of that ratio.
v(i)%x = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np) v(i)%y = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np) v(i)%z = Gaussian()*STDtemp*MassInv*ratio*sqrt((3np-ncons)/np)
This gives me the correct temperature (~380). However, when I thermalize my particles, and run an NVE trajectory of a previously equilibrated fluid (temperature should fluctuate around 380), the temperature drops to about 300. I have no idea if this is statistically mechanically correct.