Velocity and temperature for hybrid simulation

Hello there,

I am running a hybrid system with some atoms are in rigid bodies (fix shake), some aren’t.

For the velocity, I tried

velocity all create 300.0 12345
run 0 post no
velocity all scale 300.0

----------- Thermo output
Step      c_tempAll   c_tempNonrigid   c_tempRigid    PotEng      TotEng         Press
        0   349.80671      244.05641      350.35627     -242552.9      -167061.93      6795.83      
        1   349.79778      244.09576      350.34709     -242550.97     -167061.93      5547.2402    


velocity atomsNonrigid create 300.0 12345
velocity atomsRigid create 300.0 12345
run 0 post no
velocity atomsRigid scale 300.0

----------- Thermo output
Step      c_tempAll   c_tempNonrigid   c_tempRigid    PotEng      TotEng         Press   
        0   349.89941      300            350.1613      -242552.9      -167041.92      6795.7762  
        1   349.8905       300.04225      350.15212     -242550.98     -167041.93      5547.4171    

For both tests, temperatures are calculated as

compute 	tempAll all temp
compute 	tempNonrigid atomsNonrigid temp
compute 	tempRigid atomsRigid temp

I wonder which is valid and reasonable for a hybrid system, guessing is the second one that only scale velocities of the rigid bodies.



Please note that the velocity command only sets initial values. The velocity command will factor in the degrees of freedom that are removed by fix shake, but since the output of fix shake is not reliable for the first couple of steps for technical reasons, you are likely creating more problems than solving with your command sequence:

velocity atomsRigid create 300.0 12345
run 0 post no
velocity atomsRigid scale 300.0

What would be important to know is how you are time integrating the system and if you apply a thermostat to the entire system or not.
For example, if you use fix nvt (like in the peptide example) then you have just a single Nose-Hoover thermostat for the entire system (this is usually desired since thermalization is intended to be very weakly coupled). But in this scenario the actual temperature for both the rigid and non-rigid parts depends on how well those two parts can exchange kinetic energy. Most likely, it will take a significant amount of simulation time until this happens.

The alternative would be to use fix nve plus fix langevin where you thermalize individual atoms through random forces balanced by friction. This kind of dissipative thermostat can distribute kinetic energy much more efficiently, but also has much more impact on the trajectory. People often start with this setup until the expected equipartitioning is achieved and then switch to a Nose-Hoover thermostat.

Please see the graphs below for a demonstration. I have modified the “peptide” example from the LAMMPS source distribution, so that fix shake is only applied to the water molecules and then defined separate computes for temperature of all atoms, only the water atoms, and only the peptide atoms. The timestep was reduced to 0.5fs and the length of the run extended to 200000 steps.
Due to no longer having shake constraints on the hydrogen atom bonds in the peptide, the peptide atoms are not in equilibrium at the beginning and equipartitioning of the kinetic energy is not given.

As you can see from the right graph, the peptide atom temperature is lower than the target, but over the course of the simulation they slowly pick up the kinetic energy from surrounding water molecules. Since there are so many more water atoms than peptide atoms, the water temperature is only a little bit affected and fix nvt keeps the total temperature on target, but does nothing to assist the exchange of kinetic energy of the two subsystems.

I’ve then replaced fix nvt with fix nve and fix langevin with the same relaxation time as used in fix nvt. And you can see from the middle graph that the dissipative nature of fix langevin achieves equipartitioning of the temperature between the two subsystems within the first 1000 time steps.

The third graph is a mix where the first 100000 timesteps use Langevin for thermostatting and then second 100000 steps fix nvt. It is quite obvious that running fix nvt will inherit the equipartition of temperatures and not break it again.

Bottom line, getting your system in a proper equilibrium including equipartitioning to kinetic energy, may take some time. Using a dissipative thermostat can help, but one also has to factor in how much it will affect the dynamics and in most cases you want to run production with a (weak!) Nose-Hoover thermostatting. Please keep in mind that how easy the system will balance out itself depends on the geometry and details of the setup and vary by large margins.

Hi Axel,

I quite agree with your points about the thermostats.

Yes I am using nvt idea, fix nvt for the nonrigid atoms and a separated rigid/nvt for the rigid bodies (water). The timestep for the equilibration in my system is 0.5 fs and total 20000 steps, and then increases to 5fs timestep and run another 80000 steps.

After I set the initial temperature at 300 K, both temperatures for the non-rigid and rigid drop, and then increase again. This is quite strange and is it something due to the “fix shake”?



This is contradicting your original post where you wrote that you were keeping molecules rigid with fix shake. It is completely pointless to apply fix shake to atoms that are time integrated with fix rigid/nvt or vice versa. For molecules kept rigid with fix shake, you use a single fix nvt for the whole system, since with fix shake you do regular atomic position updates, that are then corrected by fix shake to fulfill the defined constraint.

5fs timestep is far too large for the system you describe. If you have water molecules (or any bonds to hydrogen atoms) kept rigid with fix shake, you may get a way with 2fs, but with fix rigid on water, you are likely stuck with 0.5fs or even less to have a proper stable and energy conserving time integration.

Not at all, what happens after the initial velocity assignment depends on how high the potential energy was of the system before. If it had been subject to a thorough minimization, it may be in a state representing a 0K configuration. In that case the temperature will drop to about half.

But what can happen in your case, if I understand your description correctly, could be that you are removing the rigid constraints twice for no good reason and that would be consistent with a temperature of about 450K for a 300K thermostat.

At any rate, it is difficult to speculate over your input and what is going wrong without actually seeing it and being able to reproduce it. Just have a look at the rhodo benchmark input deck or the peptide example input deck. Both use a rigid water potential and fix shake and start from an equilibrated structure with corresponding velocities.