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 “900*18.0154*10/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 profile.mp 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*(step*0.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