###############################################################################
#
#
#  This input script is a modified version of the example script spce_ehex.lmp
#  which is part of the supplementary (open access) material of the paper
# 
#  P. Wirnsberger, D. Frenkel and C. Dellago, 
#  "An enhanced version of the heat exchange algorithm with excellent energy 
#  conservation properties", J. Chem. Phys. 143, 124104 (2015).
#
#  The full article is available on arXiv: http://arxiv.org/pdf/1507.07081.
#
#
#  Description: 
#  ------------
# 
#  This file is a LAMMPS input script for carrying out a NEMD simulation of
#  SPC/E water using the eHEX/a algorithm. The run produces two files:
#  "out.Tspce_ehex" contains the temperature profile and "out.Espce_ehex" the time
#  evolution of the total energy.
#
###############################################################################

log                 log.spce_ehex

# energy flux into the reservoir
 variable F         	equal 0.075

# timestep              
 variable dt             equal 3.0  

# simulation time for the production run (1 ns)
 variable tprod      	equal 1000000

# total number of timesteps
 variable Nprod      	equal floor(${tprod}/${dt})

# parameters for the SPC/E model
 variable epsOO       	equal 0.15535
 variable sigOO       	equal 3.166
 variable theta       	equal 109.47

# long-range and short-range cutoffs, respectively
 variable cutC        	equal (xhi-xlo)/2.
 variable cutLJ  	equal 11

# specification of units, spatial dimensions, boundary conditions and atom-style
 units        real
 dimension    3
 boundary     p p p
 atom_style   full
 read_data    "data.spce"

# group atoms to molecules
 group    O type 2
 group    H type 1
 group    water type 1 2

# define the pair style with long-range Coulomb interaction
# and short-range LJ interaction
 pair_style    	lj/cut/coul/long ${cutLJ} ${cutC}
 pair_coeff    	2 2 ${epsOO} ${sigOO}
 pair_coeff    	1 2 0 0
 pair_coeff    	1 1 0 0

# use Ewald summation with a precision of 1.e-5
 kspace_style   ewald 1.e-5

# use harmonic bonds between sites of a molecules
# NOTE: this will not have any effects as we use RATTLE to keep the bonds fixed,
#       but it is recommended.
 bond_style    	harmonic              
 angle_style   	harmonic              
 bond_coeff    	1 1000.00 1.000
 angle_coeff   	1 100.0 ${theta}

# increase neigbor skin because of the large timestep
 neighbor  4.5 bin

# use standard correction terms for the truncated tail of the LJ potential
 pair_modify tail yes

 variable    Nsamp    equal 10
 variable    Nrepeat  equal floor(${Nprod}/${Nsamp})
 variable    Nevery   equal ${Nsamp}*${Nrepeat}
 
# compute the centre of mass velocity of the box (vcmx, vcmy, vcmz)
 variable vcmx equal "vcm(all,x)"
 variable vcmy equal "vcm(all,y)"
 variable vcmz equal "vcm(all,z)"
 variable vcm2 equal v_vcmx*v_vcmx+v_vcmy*v_vcmy+v_vcmz*v_vcmz

# compute temperature, pressure, potential energy, kinetic energy and total energy
 compute   cT  all temp
 compute   cP  all pressure thermo_temp
 compute   cPe all pe
 compute   cKe all ke
 variable  vE  equal c_cKe+c_cPe

# specify the reservoir extents
 variable Lz          equal zhi-zlo
 variable delta       equal 4 
 variable dz          equal ${Lz}/60
 variable zlo_Thi     equal -${Lz}/4.-${delta}/2.
 variable zhi_Thi     equal ${zlo_Thi}+${delta}
 variable zlo_Tlo     equal ${Lz}/4.-${delta}/2.
 variable zhi_Tlo     equal ${zlo_Tlo}+${delta}

# create regions of low and high temperature and apply thermostats
 region   Thi_region     block INF INF INF INF ${zlo_Thi} ${zhi_Thi}
 region   Tlo_region     block INF INF INF INF ${zlo_Tlo} ${zhi_Tlo}

# compute temperature of reservoirs using 3 degrees of freedom for every atom 
 compute  cTlo   water temp/region Tlo_region
 compute  cThi   water temp/region Thi_region

# rescale temperature to correct for the constraint bonds (6 instead of 9 degrees of freedom per molecule)
 variable  Tlo_act    equal c_cTlo/2*3
 variable  Thi_act    equal c_cThi/2*3

# thermostat the reservoirs using the eHEX algorithm
# NOTE: add the keyword "hex" at the end of each of the two following lines
#       if you want to use the HEX algorithm.

 fix fHi all ehex 1  ${F}   region Thi_region com constrain
 fix fLo all ehex 1 -${F}   region Tlo_region com constrain

# use velocity Verlet integration
 fix fNVE all nve

# calculate the (kinetic) temperature from the kinetic
# energy per atom
# kB is Boltzmann's constant
# NOTE: For simplicity, we do not subtract the centre of mass
#       velocity of the individual slabs in this example script. 
#       However, we did take this into account in the publication. 
#       (The differences are negligible for our setup.)

 variable  kB  equal  0.001987204
 compute   ke  water  ke/atom
 variable  T   atom   c_ke/${kB}

# use RATTLE with a precision of 1.e-10
 fix    fRattle all rattle 1e-10 400 0 b 1 a 1

# output the timestep, temperatures (average, cold reservoir, hot reservoir), energies (kinetic, potential and total),
# pressure and squared com velocity every 100 timesteps
 reset_timestep  0
 timestep        ${dt}

 thermo_style    custom step temp v_Tlo_act v_Thi_act ke pe etotal press v_vcm2
 thermo_modify   flush yes
 thermo          1000

 compute         cchT  all  chunk/atom bin/1d z lower ${dz}
 fix             fchT  all  ave/chunk  ${Nsamp} ${Nrepeat} ${Nevery} cchT v_T file out.Tspce_ehex

 fix             fE    all  ave/time 10 500 5000 v_vE file out.Espce_ehex
 run             ${Nprod}