How to use lammps to swap the molecular center-of-mass?

Dear lammps-users
There is a problem has been bothering me for a long time.I want to use Muller-Plathe algorithm (RNEMD) to calculate thermal conductivity of water.
I know fix thermal/conductivity command can swap the different atoms now, but I don’t know how to use this command to deal with the molecular like water.From literature, we know that we should swap the molecular center-of-mass velocities,rather than atomic velocities. But I really don’t know how to use lammps to swap the molecular center-of-mass velocities. Below is my in.water,I know it is wrong,because it swaps the atoms of water,rather than center-of-mass velocities ,so the energy of system is not conserved.
Looking forward to your reply!
Jack zhuang

variable den_nvt equal “90018.015410/6.022/vol”
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable T equal 298
variable V equal vol
variable dt equal 0.5


units real
dimension 3
boundary p p p
atom_style full
read_data data.water
pair_style lj/charmm/coul/long 8.0 9.0 9.0
pair_coeff * 2 0.0000 0.0000
pair_coeff 1 1 0.15535 3.166
kspace_style pppm 0.0001
bond_style harmonic
angle_style harmonic
bond_coeff 1 554.1348 1.000
angle_coeff 1 45.7696 109.47
timestep ${dt}

min_style cg
min_modify dmax 0.2
minimize 1.0e-12 1.0e-12 5000 5000

------------- equilibration and thermalization ----------------

reset_timestep 0
velocity all create $T 102486 mom yes rot yes dist gaussian

#fix NPT all npt temp $T $T 50.0 iso 1.0 1.0 1000.0
#thermo_style custom step temp etotal vol v_den_nvt
#thermo 10000
#run 200000
#unfix NPT

fix TEMP all temp/berendsen $T $T 50
fix 1 all nve
thermo_style custom step temp etotal vol

thermo 10000
run 400000

unfix TEMP

-------------- flux calculation ---------------

compute ke all ke/atom
variable temp atom c_ke/1.5/{kB}*{kCal2J}
#fix 1 all nve
fix 2 all ave/spatial 1 1000 1000 z lower 0.05 v_temp &
file units reduced
fix 3 all thermal/conductivity 10 z 20

variable tdiff equal f_2[11][3]-f_2[1][3]
thermo_style custom step temp epair etotal f_3 v_tdiff
thermo 10000
run 1000000

-------------- thermal conductivity calculation ---------------

--------------reset fix thermal/conductivity to zero energy accumulation---------------

fix 3 all thermal/conductivity 10 z 20
fix ave all ave/time 1 1 1000 v_tdiff ave running

variable k equal (f_3/(2*(step0.5-700000)lxly+1.0e-10))${kCal2J}*1.0e25/(v_tdiff/(lz/2))

thermo_style custom step temp epair etotal f_3 v_tdiff f_ave v_k
thermo 10000
run 5000000

Last time I checked the RNEMD method in lammps was designed for atomic systems not molecular ones.

This could have been for historical reasons (Lammps started as a code for strictly atomic systems) or because the developer of the package only worked with this kind of systems and not molecular ones. Therefore, you won’t see lammps swapping velocities of the COMs. Said this one can still use the method on molecular systems but extra care should be taken on how the swapping is done and how the timestep chosen. Systems with harmonic bonds tend to show poor energy conservation when swapping velocities as done within RNEMD. There are a few discussions in the archives about this.