# Short MD run. It should take less than 1 minutes on 4 cores.

# Input variables.
variable run    string tobermorite_s02_3 # Simulation name. 
variable rst    string tobermorite_s02_1 # Restart name. 
variable ts     equal  1            # timestep
variable tf     equal  300.         # Equilibrium temperature
variable p      equal  1.           # Equilibrium pressure
variable cutoff equal  12.          # For interactions in real space
variable skin   equal  2.           # Skin distance for neighbour list
variable cl     equal  1000         # correlation length for averaging
variable s      equal  5            # sample interval for averaging
variable prod   equal  5000         # Production steps

# PBC
boundary p p p

# Force field setup.
units           real
atom_style      full
bond_style      harmonic
angle_style     harmonic
pair_style      lj/cut/coul/long ${cutoff}
pair_modify     mix arithmetic
special_bonds lj/coul 0 1 1

# Neighbour list.
neighbor        ${skin} bin

# ----------------- Restart Section -----------------
read_restart ${rst}.last.restart

# Bond coefficients.
bond_coeff 1  553.9350 1.0 # O-H in tobermorite and water

# Angle coefficients.
angle_coeff 1   45.7696 109.47 # H-O-H
angle_coeff 2  30.0000 109.47
angle_coeff 3  340.0  138.00 # adapted from Heinz et al.(2005)
angle_coeff 4  340.0  102.70

# Non-bondeds.
pair_coeff 1  1  0.1554      3.16556  # O-O
pair_coeff 2 2 5.0298e-6   5.56170	# Ca-Ca hydroxide
pair_coeff 3 3 1.e-01      2.87199  # Ca-Ca aqueous (Koneshan, 1998)
pair_coeff 4  4  0.1554      3.16556  # Oh-Oh
pair_coeff 5  5  0.0         0.0      # Ho-Ho
pair_coeff 6 6 1.8405e-6   3.302027 # Si-Si
pair_coeff 7  7  0.1554      3.16556  # Ow-Ow (SPC/E water model)
pair_coeff 8  8  0.0         0.0      # Hw-Hw (SPC/E water model)
pair_coeff 9 9 0.1554      3.16556  # O-O
pair_coeff 10  10  1.3298e-6   4.7943   # Al-Al (VIII)
pair_coeff 11  11  1.8405e-6   3.7064   # Al-Al (IV)

# Long-range interactions.
kspace_style    pppm 1.e-5

# Extra variables
variable mol_unitary equal 20*3*3*2

# Groups for computing diffusion coefficients
group water type 7 8
group unit molecule <= ${mol_unitary}
group frac molecule > ${mol_unitary}
group g1 intersect water unit # structural water
group g2 intersect water frac # interlayer water

# Derived variables.
# thermo interval.
variable d       equal ${cl}*$s
# dump interval ("sample interval" frames).
variable dcycle  equal ${cl}*$s
variable tcouple equal ${ts}*100
variable pcouple equal ${ts}*1000

# Use the COM of the molecule to compute the MSD.
compute ac1 g1 chunk/atom molecule
compute ac2 g2 chunk/atom molecule
compute  msqdis1 all msd/chunk ac1
compute  msqdis2 all msd/chunk ac2
# average over all molecules.
variable msqdis1 equal ave(c_msqdis1[4])
variable msqdis2 equal ave(c_msqdis2[4])
# Same frequency as thermo.
fix sum1 all vector $d v_msqdis1
fix sum2 all vector $d v_msqdis2
variable diff1 equal slope(f_sum1)/6./($d*${ts}) # Angs2/fs.
variable diff2 equal slope(f_sum2)/6./($d*${ts})

# Output.
thermo          100
thermo_style custom step etotal evdwl ecoul elong ebond eangle ke pe temp press vol density
thermo_modify flush yes

# Trajectory.
#dump TRJ all dcd ${dcycle} ${run}.dcd
#dump_modify TRJ unwrap yes

# Thermalisation and relaxation, NPT ensemble.
timestep 	 ${ts}
if "$(is_defined(variable,ti)&&is_defined(variable,tf))" then &
"fix NPT all npt temp ${ti} ${tf} ${tcouple} tri $p $p ${pcouple}" &
"variable ti delete" &
else "fix NPT all npt temp ${tf} ${tf} ${tcouple} tri $p $p ${pcouple}"
run  500
unfix NPT