# Example for an alchemical transformation of two water molecules into a hydronium and hydroxyl ion
# WARNING: This input is intended for demonstrating the method only,
# the force field parameters are mostly made up and NOT suitable for production simulations.

# set up different names for two partitions
variable name world twowater twoions

units real
atom_style full
atom_modify map array
region box block -5 5 -5 5 -5 5
boundary p p p
create_box 2 box bond/types 2 angle/types 2 &
    extra/bond/per/atom 3 extra/angle/per/atom 3 extra/special/per/atom 3

mass 1 15.9994
mass 2 1.008

pair_style lj/cut/coul/cut 10.0
pair_coeff 1 1 0.1553 3.166
pair_coeff 1 2 0.0    1.0
pair_coeff 2 2 0.0    1.0

bond_style harmonic
bond_coeff * 1000.0 1.0

angle_style harmonic
angle_coeff * 100.0 109.47

molecule water h2o.mol

# create the two molecules we want to transform ...
create_atoms 0 single -2.0 0.0 0.0 mol water 453624
create_atoms 0 single  2.0 0.0 0.0 mol water 767353

# ... and put them in a group
group transform id 1:6

# now fill the rest of the box with more water
create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33

# change topology and settings for the two states
# we cannot simply create a different topology directly or
# load a different data file because the order and position
# of all atoms must be maintained across both replica

# we first have to remove all topology data in the transform group
delete_bonds transform bond 1
delete_bonds transform angle 1 remove

# then generate different topologies for the two partitions. select by name.
if "${name} == twowater" then &
    "create_bonds single/bond 2 1 2" &
    "create_bonds single/bond 2 1 3" &
    "create_bonds single/bond 2 4 5" &
    "create_bonds single/bond 2 4 6" &
    "create_bonds single/angle 2 2 1 3" &
    "create_bonds single/angle 2 5 4 6" &
else &
    "create_bonds single/bond 2 1 2" &
    "create_bonds single/bond 2 3 4" &
    "create_bonds single/bond 2 4 5" &
    "create_bonds single/bond 2 4 6" &
    "create_bonds single/angle 2 3 4 5" &
    "create_bonds single/angle 2 5 4 6" &
    "create_bonds single/angle 2 3 4 6" &
    "set atom 1 charge -1.1354" &
    "set atom 2 charge  0.1354" &
    "set atom 3 charge 0.56775" &
    "set atom 4 charge -0.70305" &
    "set atom 5*6 charge 0.56775"

velocity all create 300.0 5463576
timestep 0.2

# define ramp variable to combine the two different partitions
if "${name} == twowater" then             &
   "variable ramp equal ramp(1.0,0.0)"    &
else                                      &
   "variable ramp equal ramp(0.0,1.0)"

# since the trajectory and forces are kept identical through fix alchemy,
# we can do fix npt simulations, but we must use the "mixed" pressure

fix integrate all npt temp 300 300 1.0 iso 1.0 1.0 10.0
fix transform all alchemy v_ramp
compute pressure all pressure/alchemy transform
fix_modify integrate press pressure

# only need to output a dump file from one partition
# if "${name} == twowater" then &
#   "dump 1 all atom 100 ${name}.lammpstrj" &
#   "dump_modify 1 sort id"

thermo_style custom step temp press etotal density pe ke f_transform f_transform[3]
thermo_modify colname f_transform lambda colname f_transform[3] EPot_mixed
thermo_modify press pressure

thermo 100
run 20000