Disturbing equipartition in box simulation

[bvollebregt]

Dear,

I am currently trying to perform so-called ‘Adiabatic relaxation simulations’ in which the rotational relaxation process of a (in my case) di-atomic Hydrogen molecule can be studied.

The main procedure of such a simulation is as follows:

1. Build a periodic cell of L x L x L dimensions and fill with N molecules
2. Initialize velocities by sampling from the equilibrium distribution corresponding to a specified equilibrium temperature T_inf.
3. Equilibrate the system shortly (to make sure that the velocities are distributed correctly and energies are equipartitioned)
4. Disturb the balance between translational and rotational energy
5. Let the system relax again

I have been able to perform steps 1-4 and I verified that the system is in equilibrium and the energies are equipartitioned. I compute the different energies using the velocities of the atoms. I compute translational energy using the average of the velocities of two atoms of the same molecule. Rotational velocities and energy is computed by decomposing the velocity vectors in normal and (2) tangential components respectively.

Now my aim is to re-scale the CoM velocity and angular velocities separately such that the equilibrium is disturbed. To my knowledge, LAMMPS works with atomic velocities rather than with specific angular velocities.

My main question is: Does LAMMPS have features/commands to assign/scale translational (CoM) and rotational velocities separately?

If not, I suppose the one and only method to do such thing is by basically reversing my post-processing procedure, where I went from atomic velocities to CoM and angular velocities. This will probably lead to a complex system of two equations.

[akohlmey]

Here is a possible strategy (untested) for doing what you want in LAMMPS:

• use compute chunk/atom to define chunks for your molecules
• use compute chunk/vcm to get the center of mass velocities of your molecules
• use compute chunk/spread/atom (3x) to “spread” those values to all atoms
• define atom style variables that subtract the compute/spread/atom values from vx, vy, vz (the atomic velocities). after this step your atomic velocities are only representing rotations
• use the velocity command with the scale option to scale those rotations
• define another set of atom style variables that take the previously subtracted velocities and re-add them to the scaled rotational velocities with the desired scaling factor
• use the velocity command with the set option to set the atom velocities to the values from the last set of atom style variables.

The alternative (and possibly more efficient) option would be to write a custom fix style in C++ to implement those operations. You can look at fix rigid/small to see how LAMMPS would compute the translational and rotational velocities for objects defined by molecule IDs. And then look at fix temp/berendsen to see how (and when in the MD timestep) to change velocities of atoms.

[akohlmey]

FWIW, it is probably possible to streamline the commands a little bit and include the effect of the velocity scale command directly into the atom style variable computations and the just use the velocity set command with the final result.

[bvollebregt]

Great thanks Axel, I almost got your proposed method working and I can already tell that it actually works! Nevertheless, I have one question left concerning the last step of the method; assigning the scaled velocities back to the atoms.

What I actually want is

1. Let system equilibrate (NVT)
2. Use the last velocities of equilibration and scale these
3. Assign velocities to atoms
4. Equilibrate again (NVT)

What I am currently doing is

1. Let system equilibrate (NVT)
2. Each time-step I store variables with scaled velocities
3. Scaled velocities are checked in a dump file

Note that it is no big deal that I am computing the scaled velocities for all time-steps during equilibration, however when I want to assign the scaled velocities of the last time-step using
velocity H2 set v_vxScaled v_vyScaled v_vzScaled,

I obtain the following error:
Compute used in variable between runs is not current

I suppose there are two ways to fix this:

1. Only scale the velocities of the last-time step, and then assign these
2. Extract only the last time-step part of the scaled velocities

Is this right?

FYI the following part of the code all works, except for the velocity set part

---

fix NVT all nvt temp 300 300 1.0

compute cc1 H2 chunk/atom molecule
compute myChunkVCOM H2 vcm/chunk cc1

compute comchunk all chunk/spread/atom cc1 c_myChunkVCOM[*]

variable vxRot atom (vx-c_comchunk[1])
variable vyRot atom (vy-c_comchunk[2])
variable vzRot atom (vz-c_comchunk[3])

variable vxTr atom c_comchunk[1]
variable vyTr atom c_comchunk[2]
variable vzTr atom c_comchunk[3]

variable vxRotscaled atom v_vxRot * sqrt(1.5)
variable vyRotscaled atom v_vyRot * sqrt(1.5)
variable vzRotscaled atom v_vzRot * sqrt(1.5)

variable vxTrscaled atom v_vxTr * sqrt(2)
variable vyTrscaled atom v_vyTr * sqrt(2)
variable vzTrscaled atom v_vzTr * sqrt(2)

variable vxScaled atom v_vxTrscaled+v_vxRotscaled
variable vyScaled atom v_vyTrscaled+v_vyRotscaled
variable vzScaled atom v_vzTrscaled+v_vzRotscaled

dump 1 H2 custom 1000 vRotScaled.txt id mol type v_vxRotscaled v_vyRotscaled v_vzRotscaled
dump_modify 1 sort id

dump 2 H2 custom 1000 vTrScaled.txt id mol type v_vxTrscaled v_vyTrscaled v_vzTrscaled
dump_modify 2 sort id

dump 3 H2 custom 1000 vScaled.txt id mol type v_vxScaled v_vyScaled v_vzScaled
dump_modify 3 sort id

dump 4 H2 custom 1000 chunkVcom.txt id mol type c_comchunk[*]
dump_modify 4 sort id

dump 5 H2 custom 1000 velocities.txt id mol x y z vx vy vz
dump_modify 5 sort id

run 10000

unfix NVT
velocity H2 set v_vxScaled v_vyScaled v_vzScaled

---

[akohlmey]

@bvollebregt

I suppose there are two ways to fix this:

1. Only scale the velocities of the last-time step, and then assign these
2. Extract only the last time-step part of the scaled velocities

Is this right?

No. The problem is a principal problem due to some design decisions in LAMMPS.
Computes may no have an allocated storage to access the data except during a run or they may access data that needs to be accumulated during force computation. This will both work only when evaluated during a run since LAMMPS gets informed about this before the data is used. So when using the velocity set command, you get this error. However, this limitation does not apply to the averaging fixes, which are invoked during a run but keep their storage around until the next invocation. Thus a workaround is to use fix ave/atom similar to this:

variable eqsteps index 10000

fix 0 all ave/atom 1 1 ${eqsteps} v_vxScaled v_vyScaled v_vzScaled
variable vx atom f_0[1]
variable vy atom f_0[2]
variable vz atom f_0[3]

run ${eqsteps}

velocity H2 set v_vx v_vy v_vz

[bvollebregt]

@akohlmey

variable vx atom f_0[1]
variable vy atom f_0[2]
variable vz atom f_0[3]

Awesome, that indeed works and makes sense. What would be the best way to check/verify whether the set vx/vy/vz are equal to vxscaled/vyscaled/vzscaled. One way is by comparing the dump file with a file generated with write_data data.test after the velocity set command, however I suppose there is a more efficient way of directly storing/printing the "set"velocities

[akohlmey]

@bvollebregt

What would be the best way to check/verify whether the set vx/vy/vz are equal to vxscaled/vyscaled/vzscaled.

You don’t really need the original dump file, you can generate a custom one after the velocity set command with write_dump.