Hi all, I have two questions and I am asking both of them in one single post. I am trying to compute solvation free energy of methanol in water using TI. I have come up with this input file for my system
# -- Initialise system --
units real
atom_style full
# -- read data --
read_data system.data
# -- set charges --
set type 1 charge -0.8476
set type 2 charge 0.4238
set type 3 charge -0.683
set type 4 charge 0.418
set type 5 charge 0.04
set type 6 charge 0.145
# -- set variables for TI --
variable lj equal 0.3
variable lj1 equal 1
variable dlj equal 1
variable ks equal 1
variable ks1 equal 0
variable dks equal 1
# -- non bonded interactions --
pair_style lj/cut/coul/long/soft 1 0.5 0.0 10.0 10.0
pair_coeff 1 1 0.1554 3.16557 1.0
pair_coeff 2 2 0.0 0.0 1.0
pair_coeff 3 3 0.17 3.12 1.0
pair_coeff 4 4 0.0 0.0 1.0
pair_coeff 5 5 0.03 2.5 1.0
pair_coeff 6 6 0.066 3.5 1.0
# mixed rules
pair_coeff 1 2 0 0 1.0
pair_coeff 1 3 0.162536 3.1427 v_lj
pair_coeff 1 4 0 0 v_lj
pair_coeff 1 5 0.0682788 2.81317 v_lj
pair_coeff 1 6 0.101274 3.32859 v_lj
pair_coeff 2 3 0 0 v_lj
pair_coeff 2 4 0 0 v_lj
pair_coeff 2 5 0 0 v_lj
pair_coeff 2 6 0 0 v_lj
pair_coeff 3 4 0 0 1.0
pair_coeff 3 5 0.0714143 2.79285 1.0
pair_coeff 3 6 0.105925 3.30454 1.0
pair_coeff 4 5 0 0 1.0
pair_coeff 4 6 0 0 1.0
pair_coeff 5 6 0.0444972 2.95804 1.0
# -- kspace --
kspace_style pppm 0.0001
# -- special bonds --
special_bonds lj/coul 0.0 0.0 0.5
# -- bonded interactions --
bond_style harmonic
bond_coeff 1 553.0 0.945
bond_coeff 2 320.0 1.41
bond_coeff 3 340.0 1.09
bond_coeff 4 600.0 1.0
# -- angular potential --
angle_style harmonic
angle_coeff 1 55.0 108.5
angle_coeff 2 33.0 107.8
angle_coeff 3 35.0 109.5
angle_coeff 4 75.0 109.47
# -- dihedral potential ( only methanol ) --
dihedral_style opls
dihedral_coeff 1 0.0 0.0 0.352 0.0
# -- write ff data --
write_data ff.data pair ij
# -- water shake --
group spc type 1:2
fix myshake spc shake 0.0001 10 100 b 4 a 4
# -- TI --
compute TI all ti lj/cut/coul/long 3*6 v_lj1 v_dlj kspace * v_ks1 v_dks
# -- simulation protocol --
timestep 0.5
minimize 0.0 0.0 1000 10000
velocity all create 300.0 5463576
fix relax1 all npt temp 300.0 300.0 50.0 iso 1.0 1.0 100.0
run 20000
unfix relax1
fix relax2 all npt temp 300.0 300.0 50.0 iso 1.0 1.0 100.0
thermo_style custom step temp press density c_TI
thermo 10
log cg.log
run 20000
write_data prod.data nocoeff
Input structure file:
system.data (88.5 KB)
I read, in a paper to which I am comparing my results to, that soft-core potential was used in their computation. Excerpt from the paper: A linear dependence of the electrostatic interactions with the coupling parameter was imposed. For all four solutes, the soft-core function of Beuler et al. Thus, I made this input script.
I run this script and as one might expect I get an error:
ERROR: Incorrect args for pair coefficients (src/src/FEP/pair_lj_cut_coul_long_soft.cpp:604)
Last command: pair_coeff 2 2 0.0 0.0 1.0
My first question is:
1 ) How to fix this 
The second question pertains to TI.
2) Since the soft core potential already have lambda which scales the interaction, how exactly should I write the compute ti command? I have attempted to write something, does that make sense? If not, how should one go about doing it?
Thank you so much.
Best regards,
Ved