95 lines
3.0 KiB
Plaintext
95 lines
3.0 KiB
Plaintext
# Example for 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 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
|
|
|
|
# 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
|
|
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[1] f_transform[4]
|
|
thermo_modify colname f_transform[1] lambda colname f_transform[4] EPot_mixed
|
|
thermo_modify press pressure
|
|
|
|
thermo 100
|
|
run 20000
|