I am trying to run short NVE trajectories of rigid water molecules with the SHAKE algorithm. I need to create velocities from the Boltzmann distribution at the start of the run. I implement “fix shake” to remove degrees of freedom before creating the velocities via the velocity command at the desired temperature. However, after 1 time step of NVE integration, the temperature (and total energy) drop. We think this is caused by SHAKE adding forces to constrain bond lengths and angles in the first time step. For the remainder of the run, the energy remains constant. Is there a way to set the velocities for rigid molecules from the Boltzmann distribution without this energy drop in the first time step? Excerpt of my input code and a log file is below.
So long as fix shake is defined before you create
velocitities, the temperature rescaling performed
by the velocity command should include the removal
of the SHAKE degrees-of-freedom, and this shouldn't
happen. You are printing out a temperature that
includes the removal of SHAKE dof?
Scratch my previous comment. I was thinking
of another issue which leads to an incorrect
temp on step 0. Can you verify that
you get the same behavior using fix rigid
instead of fix shake? I think you're correct
that the initial components of velocity
internal to the rigid body are being removed
at the first time step. I suppose we need
an option on the velocity command to allow
for init of rigid bodies. Do you want just
the translational temperature or also the
rotational temperature?
You could always work around this by
equilibrating with fix nvt or fix rigid/nvt
or fix langevin for a short time before
running pure NVE.
I ran my code with fix rigid instead of shake and I do NOT get the drop in energy in the first time step like I did using shake. First, I have to set the velocities of all the atoms from the Boltzmann distribution and then apply fix rigid. I think (correct me if I’m wrong) that fix rigid will determine the center of mass motion (translation and rotation) based on the velocities of the water atoms. From then onward the NVE integration conserves energy. Below is the code I used and the output logfile.
I am doing hybrid MC/MD so I cannot do a short equilibration before collecting data. Fix rigid appears to allow for this type of simulation but because of the energy drop in the first time step, fix shake cannot be used for this type of method at the moment.
Also, it seems like fix rigid is a form of the settle algorithm; except settle steps forward and then corrects whereas fix rigid, in a way, corrects and then steps forward. Do you know if this is correct?
Thanks for the help,
Mark
group water molecule > 0
velocity water create 300.0 83792 mom yes rot yes dist gaussian
So long as fix shake is defined before you create
velocitities, the temperature rescaling performed
by the velocity command should include the removal
of the SHAKE degrees-of-freedom, and this shouldn't
happen. You are printing out a temperature that
includes the removal of SHAKE dof?
steve,
it looks like exactly this is not happening correctly.
if you initialize a system of water molecules held rigid with fix
shake to 450K you get 300K afterwards.
If you define the fix shake before the velocity command,
then you will get 450 on timestep 0, but not if you
flip the order of the 2 commands (since then the velocity
command doesn't know about the SHAKE degrees of freedom).
The issue that Mark is asking about is that when you then
equllibrate with NVE the system will not stay near 450.
So I looked into this a bit more for fix shake and fix rigid.
Your descriptions of what they do is basically correct, and I think
that explains why their behavior is different. Neither is
doing what you want. Though fix rigid seems better,
it is not generating body translation and rotation from the Boltz distribution.
I think the "right" solution is to initialize the velocity of a set
of rigid bodies to a desired T by setting the COM velocity
of each body to a value sampled from the Boltzmann dist
for T, and likewise the angular velocity of the body to an
omega sampled from the the Boltz for T so that all 6 DOF
for the body have been populated correctly. And then
setting the individual velocity of the atoms in the body
to the values that correspond to Vcm and omega for that body. Then
NVE integration either via fix rigid or SHAKE + fix nve should
simply conserve energy.
Such an initialization option could possibly be added to
the velocity command, but it would be non-trivial. This is
b/c there is no simple way for the velocity command to know
about the rigid bodies defined explicitly by fix rigid (their properties
are often not fully defined at the point in time a velocity command
would be used) or implicitly by fix shake.
Normally I would suggest you just equlibrate the system via fix rigid/nvt
or fix nvt or fix langevin (the latter two work fine with SHAKE). This
will rapidly produce the Boltmann distribution in Vcom and omega
that you want.
But in your case, since you don't want a thermostat, I think the
simplest solution (for me, not for you), is to define the velocity
of each atom in the data file. This means you have to pre-compute
them like I outlined above.
Mark,
Steve comment below is correct. It all burns down to the energy
equipartition theorem where each
degree of freedom accounts for 1/2KbT. In a rigid body you have 6 of
those all related to kinetic energy (translational/rotational) terms
and nothing to potential. However, finding the individual atomic
velocities by reverse engineering from the COM is in general an
underdefined problem because your system is likely to have more
unknowns variables (3N velocities) that know variables (6COM
velocities).
What you need to do is to find a matrix that yields 6 eigenvectors
corresponding to infinitesimal translations and rotations in R3. After
normalization, those eigenvectors will give you 6 "orthonormal"
directions in R3, 3 of which will correspond to pure COM translations
and the other 3 to pure rotations about the COM. As the velocities are
nothing but just time derivatives of the positions, the aforementioned
eigenvectors when chosen for the velocities will render the right
directions to be used when generating atomic velocities from the
Boltzman distribution that truly satisfy the rigid body nature of your
system. All left to be done is to assign a multiplicative factor (6 of
them) to each NORMALIZED eigenvector and then draw the values of the
multiplicative factor from the corresponding Boltzman dist at the
desired temperature.
One way to construct the matrix that I am talking about can be found
in: http://dx.doi.org/10.1063/1.454172
It is the R matrix referred to in appendix D. Just keep in mind that
in the paper the authors employ mass weighted Cartesian coordinates
thus you have to make sure you work with the eigenvectors in real
coordinates because LAMMPS only deals with pure Cartesian coordinates.
Carlos
PS: There might be an easier way to achieve your goal.
Carlos, Steve, Axel,
Thank you very much for looking into this issue and helping pointing me in the right direction to accomplish what I want with LAMMPS. I really appreciate it.
-Mark