Dear LAMMPS developers,
My system contains several polymers with rigid rings, which are surrounded by water solvent. I noticed in the LAMMPS manual that there is an explanation “For hybrid systems with both rigid bodies and solvent particles, you can thermostat only the solvent particles that surround one or more rigid bodies by appropriate choice of groups in the compute and fix commands for temperature and thermostatting. The solvent interactions with the rigid bodies should then effectively thermostat the rigid body temperature as well without use of the Langevin or Nose/Hoover options associated with the fix rigid commands.”
In my system, only the rings of the polymers are treated as separate rigid bodies, while the other part of the polymers (such as the methyl groups bonded with the rigid rings) are flexible. I built two groups, “rings” and “notring”, for the group of the rings and the other part of the system. My question is, can I still only thermostat the group “notring”? (the command lines are shown as follows)
fix 2 notring nvt temp $T T {Tdamp}
fix 3 rings rigid/nve molecule
I’m worrying about that this may not be reasonable because the group “notring” contains not only the solvent but also the flexible part of the polymers that bond to the rigid rings.
Thank you so much!
Best regards,
Jibao
Dear LAMMPS developers,
My system contains several polymers with rigid rings, which are surrounded
by water solvent. I noticed in the LAMMPS manual that there is an
explanation "For hybrid systems with both rigid bodies and solvent
particles, you can thermostat only the solvent particles that surround one
or more rigid bodies by appropriate choice of groups in the compute and fix
commands for temperature and thermostatting. The solvent interactions with
the rigid bodies should then effectively thermostat the rigid body
temperature as well without use of the Langevin or Nose/Hoover options
associated with the fix rigid commands."
In my system, only the rings of the polymers are treated as separate rigid
bodies, while the other part of the polymers (such as the methyl groups
bonded with the rigid rings) are flexible. I built two groups, "rings" and
"notring", for the group of the rings and the other part of the system. My
question is, can I still only thermostat the group "notring"? (the command
lines are shown as follows)
fix 2 notring nvt temp $T T {Tdamp}
fix 3 rings rigid/nve molecule
I'm worrying about that this may not be reasonable because the group
"notring" contains not only the solvent but also the flexible part of the
polymers that bond to the rigid rings.
fix nvt doesn't care about this. it only see atoms and their forces and
velocities. also, you can now use fix rigid/nvt to thermostat the rigid
degrees of freedom in a consistent matter. if you have many (small) rigid
rings, you may also want to consider the rigid/small versions instead of
plain rigid as the latter requires non-parallel global communication and
causes additional synchronization.
axel.
Dear Axel,
Thank you for your answer.
Another question about the temperature of the polymer system: I found that the output of the ensemble average of the temperature for the group of rigid rings was about 7 K lower than the one for the group of the solvent for a run of 10 ns in nve ensemble. Although the total energy of the system was conserved during the 10 ns simulation, the difference of the temperatures could only be reduced by further decreasing the timestep.
My question is, is the difference of the temperatures a signal that the timestep was not actually conserving the energy? Do you have an explanation of why there is a difference between the temperatures of the molecules and the solvent?
Thank you!
Jibao
How do you measure the temperature of the rings?
Steve
Dear Axel,
Thank you for your answer.
Another question about the temperature of the polymer system: I found that
the output of the ensemble average of the temperature for the group of
rigid rings was about 7 K lower than the one for the group of the solvent
for a run of 10 ns in *nve* ensemble. Although the total energy of the
system was conserved during the 10 ns simulation, the difference of the
temperatures could only be reduced by further decreasing the timestep.
rigid body integrators are more sensitive toward the choice of time step
and often require a much shorter time step than a corresponding system of
individual atoms. that primarily affects the rotational degrees of freedom.
there also is the question about how well kinetic energy is transferred
between the individual degrees of freedom and thus whether your system has,
in fact, reached equilibrium (or rather the equipartitioning condition is
reached). another open question is how you assign degrees of freedom for
such a complex system, hence steve's question.
axel.
Dear Steve,
For the temperature of the rings, I created a group of particles named “rings” including all the rings, and set a compute “compute temp_rings rings temp”; and a group of particles called “solvent” including all the solvent, and set a compute “compute temp_solvent solvent temp” for the temperature of the solvent. I found that the outputs of the ensemble average of these two temperatures are different (on the order of several degrees) in the nve ensemble. The difference depends on the time step: larger timestep results in larger difference, and sufficiently small timestep eliminates the difference. But the ensemble average of the total energy was conserved in the nve run (8 ns). Maybe it actually was not conserved, which would be shown in a much longer run.
If you use compute temp on rigid bodies, you are only left with
3 degrees-of-freedom (DOF) per body. How many rings do you
have. If you have a small #, then the “tempeature” of them could
vary by several degrees depending on whether you include/exclude
3 DOF for the center of mass. See the discussion of compute_modify
extra for details. You can also look at the output of various fix rigid
commands, which give you the temperature of the rigid bodies
and account for all of this.
There should be no problem restarting with read_restart remap so
far as I know (in the current version of LAMMPS). Though normally
you should not almost never need to use remap.
Steve
If you use compute temp on rigid bodies, you are only left with
3 degrees-of-freedom (DOF) per body. How many rings do you
have. If you have a small #, then the “tempeature” of them could
vary by several degrees depending on whether you include/exclude
3 DOF for the center of mass. See the discussion of compute_modify
extra for details. You can also look at the output of various fix rigid
commands, which give you the temperature of the rigid bodies
and account for all of this.
But that should be independent of the length of the timestep, should it not?
Axel
Dear Steve and Axel,
Thanks for your answer.
I think the the “compute temp” command in the latest version of LAMMPS has already considered the DOF of the rigid bodies. I notice that there is an explanation on this in the manual for the command: “This compute subtracts out degrees-of-freedom due to fixes that constrain molecular motion, such as fix shake and fix rigid.” And also, if there is a problem on the DOF, I think that should be independent of the length of the timestep.
Actually, the difference between temperatures of solvent and solute molecules also exists in systems without rigid bodies. Maybe this is because the timestep is not actually conserving the energy, which can only emerge in a long enough nve simulation.
Thanks,
Jibao