Dear all,
I’m doing MD simulation for a fluid system which is composed of a rigid diatomic molecule and solvent atoms.
The dynamics of the diatomic molecule is controlled by “fix rigid single” command, whereas that of the solvent atoms is by “fix nve” command.
The problem I am having is that the total linear momentum of the system is not conserved when the velocities of all atoms in the system are scaled by “velocity scale” command.
I have confirmed that the total linear momentum remains to be zero until the velocity scaling is performed.
I think that if the velocity scaling were performed correctly, the total linear momentum would be still zero (because zero*factor is still zero).
What I am suspecting is that the velocity scaling is not performed correctly on each atom in a rigid molecule.
I have found the following comment in the Lammps manual, which seems to be related to my problem:
The aggregate properties of each rigid body are calculated one time at the start of the first simulation run after this fix is specified. The properties include the position and velocity of the center-of-mass of the body, its moments of inertia, and its angular momentum. This is done using the properties of the constituent atoms of the body at that point in time (or see the infile keyword option). Thereafter, changing properties of individual atoms in the body will have no effect on a rigid body’s dynamics, unless they effect the pair_style interactions that individual particles are part of.
I have gone around this problem by writing a restart file of the diatomic molecule just before the velocity scaling.
I modified the velocity of center of mass value in the restart file in order to reflect the effect of velocity scaling, then I redefined “fix rigid single” by this new restart file.
It worked, but this requires additional modification of restart files every time velocity scaling is performed, which is inconvenient and rather ad hoc.
Could anyone suggest me for a better idea to solve this problem?
Thank you very much in advance.
Regards,
Bongsik Choi
P.S. 1. I’m using lammps-7Dec15.
P.S. 2. The following is a part of the script which I’m currently making use of.
fix 1 tracer rigid single infile dss_pn500_mr03_de60_te100_en00_init.rigid
fix 2 fluid nve
fix 2d all enforce2d
run 100000
label loopa
variable a loop 9
velocity all scale 1.0
unfix 1
fix 1 tracer rigid single infile temp$a.restart.rigid
run 100000
next a
jump in.dss_mr031_rescale loopa