Maxwellian distribution of velocities with Shake algorithm present

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.

Any thoughts?

• In the RATTLE method for ordinary MD (not hybrid MC), I think there is a velocity correction step that eliminates atomic velocities along constrained bonds. Is this kind of velocity correction step present in your HMC...? If not, it seems natural that the measured kinetic energy becomes greater than 1/2 kT * Ndof (= 3*np - ncons). – roygvib Oct 1 '15 at 18:09
• The measured KE during the course of the simulation does fluctuate around 1/2kT *NDOF(3*np-ncons). Otherwise you would never be able to equilibrate a simulation at the correct temperature set point with the shake algorithm present. Now that I have equilibrated using MD with shake present, I am trying to run umbrella sampling using HMC, and thus I am using this thermalize routine. It seeks strange to me that I would have to implement the rattle algorithm now, given that I didn't use it during the equilibration. What do you think? – user3225087 Oct 1 '15 at 21:16
• Hmm... Because I have never programmed HMC by myself, I cannot go further into details (for example, I don't know how velocities are usually treated in HMC). But if your program is already generating the correct KE, yes, RATTLE seems not necessary... To get more attention from people familiar with HMC, it may be better to attach more tags like "Monte Carlo", "constraint", "equilibration", etc :) – roygvib Oct 1 '15 at 21:37
• Ok I will add those tags! HMC isn't anything special! Its just as follows: (1) initialize momenta according to boltzmann distribution (2) run short MD trajectory (50-100 steps) (3) accept or reject trajectory based off change in hamiltonian, delta_H = delta_potential + delta_KE. So if you know how to do it in standard MD, thats what I need! – user3225087 Oct 2 '15 at 1:01
• If the velocities are handled right during the MD integration steps (i.e., velocities are kept orthogonal to the constraint surface due to chemical bonds etc), how about doing "re-scaling" all the velocities uniformly after a few steps in the NVE evolution? For example, vel(:)%x = vel(:)%x * sqrt( target_temp / get_measured_temp( vel ) ) at istep = 2 and the same for y and z. This way, it may be possible to eliminate the normal components of the velocities that were generated by the "thermalize" routine if the integrator in the code automatically eliminates such components. – roygvib Oct 3 '15 at 12:12