Example of Use¶
These commands were used to generate the output data in LAMMPS Simulation Data of Alchemical Processes.
Import Packages¶
[1]:
import numpy as np
import generate_alchemical_lammps_inputs as gali
Lennard-Jones Dimer¶
A dimer made of two Lennard-Jones beads is solvated at T*=0.75. Here we generate a section of the LAMMPS input file to compute the solvation free energy of the dimer using MBAR. Note that this output will also contain the information to calculate the solvation free energy with thermodynamic integration also.
[5]:
# Simple input where lambda pair changes
output = gali.generate_mbar_input(
[1.0, 0.0], # The parameter will change from 1 to 0
-0.2, # The parameter will be reduced by 0.2 between each step, so there will be 5 windows
[
[ # This list represents an adaptation, more than one may be listed (e.g., change charge and nonbonded at the same time although that isn't recommended)
"pair", # May be "pair" or "atom", where charges can only be changed as per atom properties
(
"lj/cut/soft", # Pair style as defined in lammps
"lambda", # Define the changing parameter to be "lambda", an option for this pair style in the LAMMPS documentation (https://docs.lammps.org/fix_adapt_fep.html)
1, # Solvent bead type (or range)
2 # Sulute bead type (or range)
)
]
],
0.75, # Temperature
# output_file="your/output/here", # Optional file to write to
n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)
n_run_prod_steps=5000000, # Duration of production run for each step (i.e., window)
)
for line in output:
print(line[:-2])
# Variables and System Condition
variable TK equal 0.7
variable freq equal 100
variable runtime_equil equal 100000
variable runtime_prod equal 500000
variable pinst equal pres
variable tinst equal tem
variable vinst equal vo
variable pe equal p
variable deltacdm equal 0.01 # delta used in central different method for derivative in T
variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T
variable nblocks equal
variable lambdas vector [1.0,0.8,0.6,0.4,0.2,0.0
thermo ${freq
# Set-up Loo
variable runid loop 1 ${nblocks} pa
label runloop
# Adjust param for the box and equilibrat
variable param equal v_lambdas[v_runid
if "${runid} == 1" then
"jump SELF skipequil
variable ind equal v_runid-
variable param0 equal v_lambdas[v_ind
variable paramramp equal ramp(v_param0,v_param
fix ADAPT all adapt/fep 1
pair lj/cut/soft lambda 1 2 v_paramramp
thermo_style custom step v_paramramp temp press pe evdwl enthalp
run ${runtime_equil} # Run Ram
unfix ADAP
fix ADAPT all adapt/fep ${freq}
pair lj/cut/soft lambda 1 2 v_param
thermo_style custom step v_param temp press pe evdwl enthalp
run ${runtime_equil} # Run Equi
label skipequil
write_data files/lj-cut-soft-lambda_${param}.dat
# Initialize compute
thermo_style custom step v_param temp press pe evdwl enthalp
compute FEPdb all fep ${TK}
pair lj/cut/soft lambda 1 2 v_deltacdm2
volume ye
compute FEPdf all fep ${TK}
pair lj/cut/soft lambda 1 2 v_deltacdm
volume ye
fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_tinst v_pinst v_vinst v_pe
c_FEPdb[1] c_FEPdf[1] file files/ti_lj-cut-soft-lambda_${param}.txt
variable delta0 equal v_lambdas[1]-v_para
compute FEP000 all fep ${TK}
pair lj/cut/soft lambda 1 2 v_delta0
volume ye
variable param000 equal v_param+v_delta
fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_pe
c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}.txt
variable delta1 equal v_lambdas[2]-v_para
compute FEP001 all fep ${TK}
pair lj/cut/soft lambda 1 2 v_delta1
volume ye
variable param001 equal v_param+v_delta
fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_pe
c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}.txt
variable delta2 equal v_lambdas[3]-v_para
compute FEP002 all fep ${TK}
pair lj/cut/soft lambda 1 2 v_delta2
volume ye
variable param002 equal v_param+v_delta
fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_pe
c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}.txt
variable delta3 equal v_lambdas[4]-v_para
compute FEP003 all fep ${TK}
pair lj/cut/soft lambda 1 2 v_delta3
volume ye
variable param003 equal v_param+v_delta
fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_pe
c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}.txt
variable delta4 equal v_lambdas[5]-v_para
compute FEP004 all fep ${TK}
pair lj/cut/soft lambda 1 2 v_delta4
volume ye
variable param004 equal v_param+v_delta
fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_pe
c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}.txt
variable delta5 equal v_lambdas[6]-v_para
compute FEP005 all fep ${TK}
pair lj/cut/soft lambda 1 2 v_delta5
volume ye
variable param005 equal v_param+v_delta
fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_pe
c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}.txt
# dump TRAJ all custom ${freq} files/dump_lj-cut-soft-lambda_${param}.lammpstrj id mol type element xu yu z
run ${runtime_prod}
uncompute FEPd
uncompute FEPd
if "${runid} != 1" then
"unfix ADAPT
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
unfix FEPou
# undump TRA
next runi
jump SELF runloop
thermo_style custom v_param temp press pe evdwl enthalp
write_data final.data nocoef
Solvated Benzene¶
Step 1: Scale Charges of Benzene to Zero¶
The package LAMMPS cannot remove the Coulombic interactions between some atoms (e.g., benzene and water) while allowing other to remain (e.g., H in benzene with C in benzene). We must then scale the charges of the entire molecule to zero.
[6]:
# Input where charges change
output = gali.generate_mbar_input(
[1.0, 0.0], # Scale charges from full strength to none
-0.2, # Decrease in steps of 0.2 (5 windows)
[
["atom", ("charge", 3, -0.115000)], # Change the atom property, charge, of atom type 3, with a full strength value of -0.115000
["atom", ("charge", 4, 0.115000)], # Change another atom property...
],
300, # Temperature in K
# output_file="your/output/here", # Optional file to write to
n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)
n_run_prod_steps=5000000, # Duration of production run for each step (i.e., window)
)
for line in output:
print(line[:-2])
# Variables and System Condition
variable TK equal 30
variable freq equal 100
variable runtime_equil equal 100000
variable runtime_prod equal 500000
variable pinst equal pres
variable tinst equal tem
variable vinst equal vo
variable pe equal p
variable deltacdm equal 0.01 # delta used in central different method for derivative in T
variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T
variable nblocks equal
variable lambdas vector [1.0,0.8,0.6,0.4,0.2,0.0
variable deltacdmA0 equal -0.115*v_deltacd
variable deltacdm2A0 equal 0.115*v_deltacd
variable deltacdmA1 equal 0.115*v_deltacd
variable deltacdm2A1 equal -0.115*v_deltacd
thermo ${freq
# Set-up Loo
variable runid loop 1 ${nblocks} pa
label runloop
# Adjust param for the box and equilibrat
variable param equal v_lambdas[v_runid
if "${runid} == 1" then
"jump SELF skipequil
variable ind equal v_runid-
variable param0 equal v_lambdas[v_ind
variable paramramp equal ramp(v_param0,v_param
variable paramrampA0 equal -0.115*ramp(v_param0,v_param
variable paramA0 equal -0.115*v_para
variable paramrampA1 equal 0.115*ramp(v_param0,v_param
variable paramA1 equal 0.115*v_para
fix ADAPT all adapt/fep 1
atom charge 3 v_paramrampA0
atom charge 4 v_paramrampA1
thermo_style custom step v_paramramp v_paramrampA0 v_paramrampA1 temp press pe evdwl enthalp
run ${runtime_equil} # Run Ram
unfix ADAP
fix ADAPT all adapt/fep ${freq}
atom charge 3 v_paramA0
atom charge 4 v_paramA1
thermo_style custom step v_param v_paramA0 v_paramA1 temp press pe evdwl enthalp
run ${runtime_equil} # Run Equi
label skipequil
write_data files/charge_${param}.dat
# Initialize compute
thermo_style custom step v_param temp press pe evdwl enthalp
compute FEPdb all fep ${TK}
atom charge 3 v_deltacdm2A0
atom charge 4 v_deltacdm2A1
volume ye
compute FEPdf all fep ${TK}
atom charge 3 v_deltacdmA0
atom charge 4 v_deltacdmA1
volume ye
fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_tinst v_pinst v_vinst v_pe
c_FEPdb[1] c_FEPdf[1] file files/ti_charge_${param}.txt
variable delta0 equal v_lambdas[1]-v_para
variable deltaA00 equal -0.115*v_delta
variable deltaA01 equal 0.115*v_delta
compute FEP000 all fep ${TK}
atom charge 3 v_deltaA00
atom charge 4 v_deltaA01
volume ye
variable param000 equal v_param+v_delta
fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_pe
c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_charge_${param}_${param000}.txt
variable delta1 equal v_lambdas[2]-v_para
variable deltaA10 equal -0.115*v_delta
variable deltaA11 equal 0.115*v_delta
compute FEP001 all fep ${TK}
atom charge 3 v_deltaA10
atom charge 4 v_deltaA11
volume ye
variable param001 equal v_param+v_delta
fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_pe
c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_charge_${param}_${param001}.txt
variable delta2 equal v_lambdas[3]-v_para
variable deltaA20 equal -0.115*v_delta
variable deltaA21 equal 0.115*v_delta
compute FEP002 all fep ${TK}
atom charge 3 v_deltaA20
atom charge 4 v_deltaA21
volume ye
variable param002 equal v_param+v_delta
fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_pe
c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_charge_${param}_${param002}.txt
variable delta3 equal v_lambdas[4]-v_para
variable deltaA30 equal -0.115*v_delta
variable deltaA31 equal 0.115*v_delta
compute FEP003 all fep ${TK}
atom charge 3 v_deltaA30
atom charge 4 v_deltaA31
volume ye
variable param003 equal v_param+v_delta
fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_pe
c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_charge_${param}_${param003}.txt
variable delta4 equal v_lambdas[5]-v_para
variable deltaA40 equal -0.115*v_delta
variable deltaA41 equal 0.115*v_delta
compute FEP004 all fep ${TK}
atom charge 3 v_deltaA40
atom charge 4 v_deltaA41
volume ye
variable param004 equal v_param+v_delta
fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_pe
c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_charge_${param}_${param004}.txt
variable delta5 equal v_lambdas[6]-v_para
variable deltaA50 equal -0.115*v_delta
variable deltaA51 equal 0.115*v_delta
compute FEP005 all fep ${TK}
atom charge 3 v_deltaA50
atom charge 4 v_deltaA51
volume ye
variable param005 equal v_param+v_delta
fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_pe
c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_charge_${param}_${param005}.txt
# dump TRAJ all custom ${freq} files/dump_charge_${param}.lammpstrj id mol type element xu yu z
run ${runtime_prod}
uncompute FEPd
uncompute FEPd
if "${runid} != 1" then
"unfix ADAPT
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
unfix FEPou
# undump TRA
next runi
jump SELF runloop
thermo_style custom v_param temp press pe evdwl enthalp
write_data final.data nocoef
Step 2a: Remove Non-bonded Interaction between Benzene and Water without Trajectory¶
Because saving the trajectories for a free energy calculation can be quite large, this is the preferred method. However, if later you determine that you need an additional value of lambda, you’ll need to rerun the entire calculation.
[7]:
# Input where pair changes after charges change
output = gali.generate_mbar_input(
[1.0, 0.0],
-0.0625, # Remove within 16 steps
[
[
"pair", # change a pair style (in the previous step we changed an atom property)
(
"lj/cut/soft",
"lambda",
"1*2", # Range of atom types for water (in format required by LAMMPS)
"3*4" # Range of atom types for the solute, benzene (in format required by LAMMPS)
)
]
],
300,
parameter_array=np.array([1,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0]), # Overwrite parameter_change with custom lambda spacing
fix_adapt_changes2 = [ # Keep the change in charges done in the previous step, so defined again here
["atom", ("charge", 3, -0.115000)],
["atom", ("charge", 4, 0.115000)],
],
parameter2_value=0, # Define the final scaling factor for the parameter change in step 1
# output_file="your/output/here", # Optional file to write to
n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)
n_run_prod_steps=5000000, # Duration of production run for each step (i.e., window)
)
for line in output:
print(line[:-2])
# Variables and System Condition
variable TK equal 30
variable freq equal 100
variable runtime_equil equal 100000
variable runtime_prod equal 500000
variable pinst equal pres
variable tinst equal tem
variable vinst equal vo
variable pe equal p
variable deltacdm equal 0.01 # delta used in central different method for derivative in T
variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T
variable nblocks equal 1
variable lambdas vector [1.0,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0
# Set Previous Chang
variable param2 equal
variable param2A0 equal -0.115*v_param
variable param2A1 equal 0.115*v_param
fix ADAPT2 all adapt/fep 1
atom charge 3 v_param2A0
atom charge 4 v_param2A1
thermo ${freq
# Set-up Loo
variable runid loop 1 ${nblocks} pa
label runloop
# Adjust param for the box and equilibrat
variable param equal v_lambdas[v_runid
if "${runid} == 1" then
"jump SELF skipequil
variable ind equal v_runid-
variable param0 equal v_lambdas[v_ind
variable paramramp equal ramp(v_param0,v_param
fix ADAPT all adapt/fep 1
pair lj/cut/soft lambda 1*2 3*4 v_paramramp
thermo_style custom step v_paramramp temp press pe evdwl enthalp
run ${runtime_equil} # Run Ram
unfix ADAP
fix ADAPT all adapt/fep ${freq}
pair lj/cut/soft lambda 1*2 3*4 v_param
thermo_style custom step v_param temp press pe evdwl enthalp
run ${runtime_equil} # Run Equi
label skipequil
write_data files/lj-cut-soft-lambda_${param}.dat
# Initialize compute
thermo_style custom step v_param temp press pe evdwl enthalp
compute FEPdb all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_deltacdm2
volume ye
compute FEPdf all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_deltacdm
volume ye
fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_param2 v_tinst v_pinst v_vinst v_pe
c_FEPdb[1] c_FEPdf[1] file files/ti_lj-cut-soft-lambda_${param}_charge--0.0_0.txt
variable delta0 equal v_lambdas[1]-v_para
compute FEP000 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta0
volume ye
variable param000 equal v_param+v_delta
fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_param2 v_pe
c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}_charge--0.0.txt
variable delta1 equal v_lambdas[2]-v_para
compute FEP001 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta1
volume ye
variable param001 equal v_param+v_delta
fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_param2 v_pe
c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}_charge--0.0.txt
variable delta2 equal v_lambdas[3]-v_para
compute FEP002 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta2
volume ye
variable param002 equal v_param+v_delta
fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_param2 v_pe
c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}_charge--0.0.txt
variable delta3 equal v_lambdas[4]-v_para
compute FEP003 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta3
volume ye
variable param003 equal v_param+v_delta
fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_param2 v_pe
c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}_charge--0.0.txt
variable delta4 equal v_lambdas[5]-v_para
compute FEP004 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta4
volume ye
variable param004 equal v_param+v_delta
fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_param2 v_pe
c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}_charge--0.0.txt
variable delta5 equal v_lambdas[6]-v_para
compute FEP005 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta5
volume ye
variable param005 equal v_param+v_delta
fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_param2 v_pe
c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}_charge--0.0.txt
variable delta6 equal v_lambdas[7]-v_para
compute FEP006 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta6
volume ye
variable param006 equal v_param+v_delta
fix FEPout006 all ave/time ${freq} 1 ${freq} v_param v_param006 v_param2 v_pe
c_FEP006[1] c_FEP006[2] c_FEP006[3] file files/mbar_lj-cut-soft-lambda_${param}_${param006}_charge--0.0.txt
variable delta7 equal v_lambdas[8]-v_para
compute FEP007 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta7
volume ye
variable param007 equal v_param+v_delta
fix FEPout007 all ave/time ${freq} 1 ${freq} v_param v_param007 v_param2 v_pe
c_FEP007[1] c_FEP007[2] c_FEP007[3] file files/mbar_lj-cut-soft-lambda_${param}_${param007}_charge--0.0.txt
variable delta8 equal v_lambdas[9]-v_para
compute FEP008 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta8
volume ye
variable param008 equal v_param+v_delta
fix FEPout008 all ave/time ${freq} 1 ${freq} v_param v_param008 v_param2 v_pe
c_FEP008[1] c_FEP008[2] c_FEP008[3] file files/mbar_lj-cut-soft-lambda_${param}_${param008}_charge--0.0.txt
variable delta9 equal v_lambdas[10]-v_para
compute FEP009 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta9
volume ye
variable param009 equal v_param+v_delta
fix FEPout009 all ave/time ${freq} 1 ${freq} v_param v_param009 v_param2 v_pe
c_FEP009[1] c_FEP009[2] c_FEP009[3] file files/mbar_lj-cut-soft-lambda_${param}_${param009}_charge--0.0.txt
variable delta10 equal v_lambdas[11]-v_para
compute FEP010 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta10
volume ye
variable param010 equal v_param+v_delta1
fix FEPout010 all ave/time ${freq} 1 ${freq} v_param v_param010 v_param2 v_pe
c_FEP010[1] c_FEP010[2] c_FEP010[3] file files/mbar_lj-cut-soft-lambda_${param}_${param010}_charge--0.0.txt
variable delta11 equal v_lambdas[12]-v_para
compute FEP011 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta11
volume ye
variable param011 equal v_param+v_delta1
fix FEPout011 all ave/time ${freq} 1 ${freq} v_param v_param011 v_param2 v_pe
c_FEP011[1] c_FEP011[2] c_FEP011[3] file files/mbar_lj-cut-soft-lambda_${param}_${param011}_charge--0.0.txt
variable delta12 equal v_lambdas[13]-v_para
compute FEP012 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta12
volume ye
variable param012 equal v_param+v_delta1
fix FEPout012 all ave/time ${freq} 1 ${freq} v_param v_param012 v_param2 v_pe
c_FEP012[1] c_FEP012[2] c_FEP012[3] file files/mbar_lj-cut-soft-lambda_${param}_${param012}_charge--0.0.txt
variable delta13 equal v_lambdas[14]-v_para
compute FEP013 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta13
volume ye
variable param013 equal v_param+v_delta1
fix FEPout013 all ave/time ${freq} 1 ${freq} v_param v_param013 v_param2 v_pe
c_FEP013[1] c_FEP013[2] c_FEP013[3] file files/mbar_lj-cut-soft-lambda_${param}_${param013}_charge--0.0.txt
variable delta14 equal v_lambdas[15]-v_para
compute FEP014 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta14
volume ye
variable param014 equal v_param+v_delta1
fix FEPout014 all ave/time ${freq} 1 ${freq} v_param v_param014 v_param2 v_pe
c_FEP014[1] c_FEP014[2] c_FEP014[3] file files/mbar_lj-cut-soft-lambda_${param}_${param014}_charge--0.0.txt
variable delta15 equal v_lambdas[16]-v_para
compute FEP015 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta15
volume ye
variable param015 equal v_param+v_delta1
fix FEPout015 all ave/time ${freq} 1 ${freq} v_param v_param015 v_param2 v_pe
c_FEP015[1] c_FEP015[2] c_FEP015[3] file files/mbar_lj-cut-soft-lambda_${param}_${param015}_charge--0.0.txt
# dump TRAJ all custom ${freq} files/dump_lj-cut-soft-lambda_${param}_charge--0.0_0.lammpstrj id mol type element xu yu z
run ${runtime_prod}
uncompute FEPd
uncompute FEPd
if "${runid} != 1" then
"unfix ADAPT
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
unfix FEPou
# undump TRA
next runi
jump SELF runloop
thermo_style custom v_param temp press pe evdwl enthalp
write_data final.data nocoef
Step 2b: Remove Non-bonded Interaction between Benzene and Water with Trajectory¶
Generate the input file to write the trajectories needed to calculate the free energy of solvation. In the next step, we generate the input to calculate the values needed for MBAR FOR ONE WINDOW using rerun and the trajectories. This might be desirable if you don’t know the number of lambda windows or spacing that your system requires.
[8]:
# Input where pair changes after charges change
ouptut = gali.generate_traj_input(
[1.0, 0.0],
-0.0625,
[
["pair", ("lj/cut/soft", "lambda", "1*2", "3*4")]
],
300,
parameter_array=np.array([1,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0]), # Overwrite parameter_change with custom lambda spacing
fix_adapt_changes2 = [
["atom", ("charge", 3, -0.115000)],
["atom", ("charge", 4, 0.115000)],
],
parameter2_value=0,
# output_file="your/output/here", # Optional file to write to
n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)
n_run_prod_steps=5000000, # Duration of production run for each step (i.e., window)
)
for line in output:
print(line[:-2])
# Variables and System Condition
variable TK equal 30
variable freq equal 100
variable runtime_equil equal 100000
variable runtime_prod equal 500000
variable pinst equal pres
variable tinst equal tem
variable vinst equal vo
variable pe equal p
variable deltacdm equal 0.01 # delta used in central different method for derivative in T
variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T
variable nblocks equal 1
variable lambdas vector [1.0,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0
# Set Previous Chang
variable param2 equal
variable param2A0 equal -0.115*v_param
variable param2A1 equal 0.115*v_param
fix ADAPT2 all adapt/fep 1
atom charge 3 v_param2A0
atom charge 4 v_param2A1
thermo ${freq
# Set-up Loo
variable runid loop 1 ${nblocks} pa
label runloop
# Adjust param for the box and equilibrat
variable param equal v_lambdas[v_runid
if "${runid} == 1" then
"jump SELF skipequil
variable ind equal v_runid-
variable param0 equal v_lambdas[v_ind
variable paramramp equal ramp(v_param0,v_param
fix ADAPT all adapt/fep 1
pair lj/cut/soft lambda 1*2 3*4 v_paramramp
thermo_style custom step v_paramramp temp press pe evdwl enthalp
run ${runtime_equil} # Run Ram
unfix ADAP
fix ADAPT all adapt/fep ${freq}
pair lj/cut/soft lambda 1*2 3*4 v_param
thermo_style custom step v_param temp press pe evdwl enthalp
run ${runtime_equil} # Run Equi
label skipequil
write_data files/lj-cut-soft-lambda_${param}.dat
# Initialize compute
thermo_style custom step v_param temp press pe evdwl enthalp
compute FEPdb all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_deltacdm2
volume ye
compute FEPdf all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_deltacdm
volume ye
fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_param2 v_tinst v_pinst v_vinst v_pe
c_FEPdb[1] c_FEPdf[1] file files/ti_lj-cut-soft-lambda_${param}_charge--0.0_0.txt
variable delta0 equal v_lambdas[1]-v_para
compute FEP000 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta0
volume ye
variable param000 equal v_param+v_delta
fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_param2 v_pe
c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}_charge--0.0.txt
variable delta1 equal v_lambdas[2]-v_para
compute FEP001 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta1
volume ye
variable param001 equal v_param+v_delta
fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_param2 v_pe
c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}_charge--0.0.txt
variable delta2 equal v_lambdas[3]-v_para
compute FEP002 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta2
volume ye
variable param002 equal v_param+v_delta
fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_param2 v_pe
c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}_charge--0.0.txt
variable delta3 equal v_lambdas[4]-v_para
compute FEP003 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta3
volume ye
variable param003 equal v_param+v_delta
fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_param2 v_pe
c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}_charge--0.0.txt
variable delta4 equal v_lambdas[5]-v_para
compute FEP004 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta4
volume ye
variable param004 equal v_param+v_delta
fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_param2 v_pe
c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}_charge--0.0.txt
variable delta5 equal v_lambdas[6]-v_para
compute FEP005 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta5
volume ye
variable param005 equal v_param+v_delta
fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_param2 v_pe
c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}_charge--0.0.txt
variable delta6 equal v_lambdas[7]-v_para
compute FEP006 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta6
volume ye
variable param006 equal v_param+v_delta
fix FEPout006 all ave/time ${freq} 1 ${freq} v_param v_param006 v_param2 v_pe
c_FEP006[1] c_FEP006[2] c_FEP006[3] file files/mbar_lj-cut-soft-lambda_${param}_${param006}_charge--0.0.txt
variable delta7 equal v_lambdas[8]-v_para
compute FEP007 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta7
volume ye
variable param007 equal v_param+v_delta
fix FEPout007 all ave/time ${freq} 1 ${freq} v_param v_param007 v_param2 v_pe
c_FEP007[1] c_FEP007[2] c_FEP007[3] file files/mbar_lj-cut-soft-lambda_${param}_${param007}_charge--0.0.txt
variable delta8 equal v_lambdas[9]-v_para
compute FEP008 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta8
volume ye
variable param008 equal v_param+v_delta
fix FEPout008 all ave/time ${freq} 1 ${freq} v_param v_param008 v_param2 v_pe
c_FEP008[1] c_FEP008[2] c_FEP008[3] file files/mbar_lj-cut-soft-lambda_${param}_${param008}_charge--0.0.txt
variable delta9 equal v_lambdas[10]-v_para
compute FEP009 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta9
volume ye
variable param009 equal v_param+v_delta
fix FEPout009 all ave/time ${freq} 1 ${freq} v_param v_param009 v_param2 v_pe
c_FEP009[1] c_FEP009[2] c_FEP009[3] file files/mbar_lj-cut-soft-lambda_${param}_${param009}_charge--0.0.txt
variable delta10 equal v_lambdas[11]-v_para
compute FEP010 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta10
volume ye
variable param010 equal v_param+v_delta1
fix FEPout010 all ave/time ${freq} 1 ${freq} v_param v_param010 v_param2 v_pe
c_FEP010[1] c_FEP010[2] c_FEP010[3] file files/mbar_lj-cut-soft-lambda_${param}_${param010}_charge--0.0.txt
variable delta11 equal v_lambdas[12]-v_para
compute FEP011 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta11
volume ye
variable param011 equal v_param+v_delta1
fix FEPout011 all ave/time ${freq} 1 ${freq} v_param v_param011 v_param2 v_pe
c_FEP011[1] c_FEP011[2] c_FEP011[3] file files/mbar_lj-cut-soft-lambda_${param}_${param011}_charge--0.0.txt
variable delta12 equal v_lambdas[13]-v_para
compute FEP012 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta12
volume ye
variable param012 equal v_param+v_delta1
fix FEPout012 all ave/time ${freq} 1 ${freq} v_param v_param012 v_param2 v_pe
c_FEP012[1] c_FEP012[2] c_FEP012[3] file files/mbar_lj-cut-soft-lambda_${param}_${param012}_charge--0.0.txt
variable delta13 equal v_lambdas[14]-v_para
compute FEP013 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta13
volume ye
variable param013 equal v_param+v_delta1
fix FEPout013 all ave/time ${freq} 1 ${freq} v_param v_param013 v_param2 v_pe
c_FEP013[1] c_FEP013[2] c_FEP013[3] file files/mbar_lj-cut-soft-lambda_${param}_${param013}_charge--0.0.txt
variable delta14 equal v_lambdas[15]-v_para
compute FEP014 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta14
volume ye
variable param014 equal v_param+v_delta1
fix FEPout014 all ave/time ${freq} 1 ${freq} v_param v_param014 v_param2 v_pe
c_FEP014[1] c_FEP014[2] c_FEP014[3] file files/mbar_lj-cut-soft-lambda_${param}_${param014}_charge--0.0.txt
variable delta15 equal v_lambdas[16]-v_para
compute FEP015 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta15
volume ye
variable param015 equal v_param+v_delta1
fix FEPout015 all ave/time ${freq} 1 ${freq} v_param v_param015 v_param2 v_pe
c_FEP015[1] c_FEP015[2] c_FEP015[3] file files/mbar_lj-cut-soft-lambda_${param}_${param015}_charge--0.0.txt
# dump TRAJ all custom ${freq} files/dump_lj-cut-soft-lambda_${param}_charge--0.0_0.lammpstrj id mol type element xu yu z
run ${runtime_prod}
uncompute FEPd
uncompute FEPd
if "${runid} != 1" then
"unfix ADAPT
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
uncompute FEP01
unfix FEPout01
unfix FEPou
# undump TRA
next runi
jump SELF runloop
thermo_style custom v_param temp press pe evdwl enthalp
write_data final.data nocoef
[9]:
output = gali.generate_rerun_mbar(
0.5, # Value of lambda to make a rerun script for
[1.0, 0.0],
-0.0625,
[
["pair", ("lj/cut/soft", "lambda", "1*2", "3*4")]
],
300,
parameter_array=np.array([1,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0]),
fix_adapt_changes2 = [
["atom", ("charge", 3, -0.115000)],
["atom", ("charge", 4, 0.115000)],
],
parameter2_value=0,
# output_file="your/output/here", # Optional file to write to
)
for line in output:
print(line[:-2])
# Variables and System Condition
variable TK equal 30
variable freq equal 100
variable pinst equal pres
variable tinst equal tem
variable vinst equal vo
variable pe equal p
variable param equal 0.
variable lambdas vector [1.0,0.95,0.9,0.85,0.8,0.75,0.7,0.65,0.6,0.5,0.4,0.3,0.2,0.1,0.05,0.0
thermo ${freq
# Set Previous Chang
variable param2 equal
variable param2A0 equal -0.115*v_param
variable param2A1 equal 0.115*v_param
fix ADAPT2 all adapt/fep 1
atom charge 3 v_param2A0
atom charge 4 v_param2A1
# Initialize compute
fix ADAPT all adapt/fep ${freq}
pair lj/cut/soft lambda 1*2 3*4 v_param
variable delta0 equal 0.
compute FEP000 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta0
volume ye
variable param000 equal v_param+v_delta
fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_param2 v_pe
c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_lj-cut-soft-lambda_${param}_${param000}_charge--0.0.txt
variable delta1 equal 0.4
compute FEP001 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta1
volume ye
variable param001 equal v_param+v_delta
fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_param2 v_pe
c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_lj-cut-soft-lambda_${param}_${param001}_charge--0.0.txt
variable delta2 equal 0.
compute FEP002 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta2
volume ye
variable param002 equal v_param+v_delta
fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_param2 v_pe
c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_lj-cut-soft-lambda_${param}_${param002}_charge--0.0.txt
variable delta3 equal 0.3
compute FEP003 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta3
volume ye
variable param003 equal v_param+v_delta
fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_param2 v_pe
c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_lj-cut-soft-lambda_${param}_${param003}_charge--0.0.txt
variable delta4 equal 0.
compute FEP004 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta4
volume ye
variable param004 equal v_param+v_delta
fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_param2 v_pe
c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_lj-cut-soft-lambda_${param}_${param004}_charge--0.0.txt
variable delta5 equal 0.2
compute FEP005 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta5
volume ye
variable param005 equal v_param+v_delta
fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_param2 v_pe
c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_lj-cut-soft-lambda_${param}_${param005}_charge--0.0.txt
variable delta6 equal 0.
compute FEP006 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta6
volume ye
variable param006 equal v_param+v_delta
fix FEPout006 all ave/time ${freq} 1 ${freq} v_param v_param006 v_param2 v_pe
c_FEP006[1] c_FEP006[2] c_FEP006[3] file files/mbar_lj-cut-soft-lambda_${param}_${param006}_charge--0.0.txt
variable delta7 equal 0.1
compute FEP007 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta7
volume ye
variable param007 equal v_param+v_delta
fix FEPout007 all ave/time ${freq} 1 ${freq} v_param v_param007 v_param2 v_pe
c_FEP007[1] c_FEP007[2] c_FEP007[3] file files/mbar_lj-cut-soft-lambda_${param}_${param007}_charge--0.0.txt
variable delta8 equal 0.
compute FEP008 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta8
volume ye
variable param008 equal v_param+v_delta
fix FEPout008 all ave/time ${freq} 1 ${freq} v_param v_param008 v_param2 v_pe
c_FEP008[1] c_FEP008[2] c_FEP008[3] file files/mbar_lj-cut-soft-lambda_${param}_${param008}_charge--0.0.txt
variable delta9 equal 0.
compute FEP009 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta9
volume ye
variable param009 equal v_param+v_delta
fix FEPout009 all ave/time ${freq} 1 ${freq} v_param v_param009 v_param2 v_pe
c_FEP009[1] c_FEP009[2] c_FEP009[3] file files/mbar_lj-cut-soft-lambda_${param}_${param009}_charge--0.0.txt
variable delta10 equal -0.
compute FEP010 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta10
volume ye
variable param010 equal v_param+v_delta1
fix FEPout010 all ave/time ${freq} 1 ${freq} v_param v_param010 v_param2 v_pe
c_FEP010[1] c_FEP010[2] c_FEP010[3] file files/mbar_lj-cut-soft-lambda_${param}_${param010}_charge--0.0.txt
variable delta11 equal -0.
compute FEP011 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta11
volume ye
variable param011 equal v_param+v_delta1
fix FEPout011 all ave/time ${freq} 1 ${freq} v_param v_param011 v_param2 v_pe
c_FEP011[1] c_FEP011[2] c_FEP011[3] file files/mbar_lj-cut-soft-lambda_${param}_${param011}_charge--0.0.txt
variable delta12 equal -0.
compute FEP012 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta12
volume ye
variable param012 equal v_param+v_delta1
fix FEPout012 all ave/time ${freq} 1 ${freq} v_param v_param012 v_param2 v_pe
c_FEP012[1] c_FEP012[2] c_FEP012[3] file files/mbar_lj-cut-soft-lambda_${param}_${param012}_charge--0.0.txt
variable delta13 equal -0.
compute FEP013 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta13
volume ye
variable param013 equal v_param+v_delta1
fix FEPout013 all ave/time ${freq} 1 ${freq} v_param v_param013 v_param2 v_pe
c_FEP013[1] c_FEP013[2] c_FEP013[3] file files/mbar_lj-cut-soft-lambda_${param}_${param013}_charge--0.0.txt
variable delta14 equal -0.4
compute FEP014 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta14
volume ye
variable param014 equal v_param+v_delta1
fix FEPout014 all ave/time ${freq} 1 ${freq} v_param v_param014 v_param2 v_pe
c_FEP014[1] c_FEP014[2] c_FEP014[3] file files/mbar_lj-cut-soft-lambda_${param}_${param014}_charge--0.0.txt
variable delta15 equal -0.
compute FEP015 all fep ${TK}
pair lj/cut/soft lambda 1*2 3*4 v_delta15
volume ye
variable param015 equal v_param+v_delta1
fix FEPout015 all ave/time ${freq} 1 ${freq} v_param v_param015 v_param2 v_pe
c_FEP015[1] c_FEP015[2] c_FEP015[3] file files/mbar_lj-cut-soft-lambda_${param}_${param015}_charge--0.0.txt
rerun files/dump_lj-cut-soft-lambda_${param}_charge--0.0_0.lammpstrj every ${freq} dump x y z
Step 3: Restore the Charges of Benzene Atoms in Vacuum¶
Because the previous two steps were missing the intramolecular charge interactions of benzene, and the end of the second system resulted in benzene being completely decoupled from water, we must now determine the free energy of returning the charges to complete the cycle.
[10]:
# Input where charges change
output = gali.generate_mbar_input(
[0.0, 1.0], # Scale the charges from zero to in full force
0.2,
[
["atom", ("charge", 1, -0.115000)],
["atom", ("charge", 2, 0.115000)],
],
300,
# output_file="your/output/here", # Optional file to write to
n_run_equil_steps=1000000, # Duration of equilibration run for each step (i.e., window)
n_run_prod_steps=5000000, # Duration of production run for each step (i.e., window)
)
for line in output:
print(line[:-2])
# Variables and System Condition
variable TK equal 30
variable freq equal 100
variable runtime_equil equal 100000
variable runtime_prod equal 500000
variable pinst equal pres
variable tinst equal tem
variable vinst equal vo
variable pe equal p
variable deltacdm equal 0.01 # delta used in central different method for derivative in T
variable deltacdm2 equal -0.01 # delta used in central different method for derivative in T
variable nblocks equal
variable lambdas vector [0.0,0.2,0.4,0.6,0.8,1.0
variable deltacdmA0 equal -0.115*v_deltacd
variable deltacdm2A0 equal 0.115*v_deltacd
variable deltacdmA1 equal 0.115*v_deltacd
variable deltacdm2A1 equal -0.115*v_deltacd
thermo ${freq
# Set-up Loo
variable runid loop 1 ${nblocks} pa
label runloop
# Adjust param for the box and equilibrat
variable param equal v_lambdas[v_runid
if "${runid} == 1" then
"jump SELF skipequil
variable ind equal v_runid-
variable param0 equal v_lambdas[v_ind
variable paramramp equal ramp(v_param0,v_param
variable paramrampA0 equal -0.115*ramp(v_param0,v_param
variable paramA0 equal -0.115*v_para
variable paramrampA1 equal 0.115*ramp(v_param0,v_param
variable paramA1 equal 0.115*v_para
fix ADAPT all adapt/fep 1
atom charge 1 v_paramrampA0
atom charge 2 v_paramrampA1
thermo_style custom step v_paramramp v_paramrampA0 v_paramrampA1 temp press pe evdwl enthalp
run ${runtime_equil} # Run Ram
unfix ADAP
fix ADAPT all adapt/fep ${freq}
atom charge 1 v_paramA0
atom charge 2 v_paramA1
thermo_style custom step v_param v_paramA0 v_paramA1 temp press pe evdwl enthalp
run ${runtime_equil} # Run Equi
label skipequil
write_data files/charge_${param}.dat
# Initialize compute
thermo_style custom step v_param temp press pe evdwl enthalp
compute FEPdb all fep ${TK}
atom charge 1 v_deltacdm2A0
atom charge 2 v_deltacdm2A1
volume ye
compute FEPdf all fep ${TK}
atom charge 1 v_deltacdmA0
atom charge 2 v_deltacdmA1
volume ye
fix FEPout all ave/time ${freq} 1 ${freq} v_param v_deltacdm v_tinst v_pinst v_vinst v_pe
c_FEPdb[1] c_FEPdf[1] file files/ti_charge_${param}.txt
variable delta0 equal v_lambdas[1]-v_para
variable deltaA00 equal -0.115*v_delta
variable deltaA01 equal 0.115*v_delta
compute FEP000 all fep ${TK}
atom charge 1 v_deltaA00
atom charge 2 v_deltaA01
volume ye
variable param000 equal v_param+v_delta
fix FEPout000 all ave/time ${freq} 1 ${freq} v_param v_param000 v_pe
c_FEP000[1] c_FEP000[2] c_FEP000[3] file files/mbar_charge_${param}_${param000}.txt
variable delta1 equal v_lambdas[2]-v_para
variable deltaA10 equal -0.115*v_delta
variable deltaA11 equal 0.115*v_delta
compute FEP001 all fep ${TK}
atom charge 1 v_deltaA10
atom charge 2 v_deltaA11
volume ye
variable param001 equal v_param+v_delta
fix FEPout001 all ave/time ${freq} 1 ${freq} v_param v_param001 v_pe
c_FEP001[1] c_FEP001[2] c_FEP001[3] file files/mbar_charge_${param}_${param001}.txt
variable delta2 equal v_lambdas[3]-v_para
variable deltaA20 equal -0.115*v_delta
variable deltaA21 equal 0.115*v_delta
compute FEP002 all fep ${TK}
atom charge 1 v_deltaA20
atom charge 2 v_deltaA21
volume ye
variable param002 equal v_param+v_delta
fix FEPout002 all ave/time ${freq} 1 ${freq} v_param v_param002 v_pe
c_FEP002[1] c_FEP002[2] c_FEP002[3] file files/mbar_charge_${param}_${param002}.txt
variable delta3 equal v_lambdas[4]-v_para
variable deltaA30 equal -0.115*v_delta
variable deltaA31 equal 0.115*v_delta
compute FEP003 all fep ${TK}
atom charge 1 v_deltaA30
atom charge 2 v_deltaA31
volume ye
variable param003 equal v_param+v_delta
fix FEPout003 all ave/time ${freq} 1 ${freq} v_param v_param003 v_pe
c_FEP003[1] c_FEP003[2] c_FEP003[3] file files/mbar_charge_${param}_${param003}.txt
variable delta4 equal v_lambdas[5]-v_para
variable deltaA40 equal -0.115*v_delta
variable deltaA41 equal 0.115*v_delta
compute FEP004 all fep ${TK}
atom charge 1 v_deltaA40
atom charge 2 v_deltaA41
volume ye
variable param004 equal v_param+v_delta
fix FEPout004 all ave/time ${freq} 1 ${freq} v_param v_param004 v_pe
c_FEP004[1] c_FEP004[2] c_FEP004[3] file files/mbar_charge_${param}_${param004}.txt
variable delta5 equal v_lambdas[6]-v_para
variable deltaA50 equal -0.115*v_delta
variable deltaA51 equal 0.115*v_delta
compute FEP005 all fep ${TK}
atom charge 1 v_deltaA50
atom charge 2 v_deltaA51
volume ye
variable param005 equal v_param+v_delta
fix FEPout005 all ave/time ${freq} 1 ${freq} v_param v_param005 v_pe
c_FEP005[1] c_FEP005[2] c_FEP005[3] file files/mbar_charge_${param}_${param005}.txt
# dump TRAJ all custom ${freq} files/dump_charge_${param}.lammpstrj id mol type element xu yu z
run ${runtime_prod}
uncompute FEPd
uncompute FEPd
if "${runid} != 1" then
"unfix ADAPT
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
uncompute FEP00
unfix FEPout00
unfix FEPou
# undump TRA
next runi
jump SELF runloop
thermo_style custom v_param temp press pe evdwl enthalp
write_data final.data nocoef