Hi Everyone, I have used a geometric pair style to compute Gibbs free energy through the free energy perturbation method. My syntax to compute fep is showing an error. I am struggling to compute fep for large systems. I am attaching my input script. could you please check what is creating the error?
Atom 2, 3, 4, 12, 13, 14, 15, 16, 17, 18, 19, 29, 30 is a molecule. And I want to calculate free energy of molecule in the system.
units real
atom_style full
dimension 3
boundary p p p
pair_style lj/cut/coul/long 13.0
pair_modify tail yes
pair_modify mix geometric
bond_style harmonic
angle_style harmonic
dihedral_style opls
read_data freentf2.data
set type 1 charge 0.35 #C
set type 2 charge 0.0083 #CAA
set type 3 charge -0.1981 #CBB
set type 4 charge 0.0085 #CCC
set type 5 charge -0.17 #Ca
set type 6 charge -0.35 #Cm
set type 7 charge -0.09 #Cr
set type 8 charge -0.12 #Cs
set type 9 charge -0.24 #ct
set type 10 charge -0.24 #cw
set type 11 charge -0.16 #F
set type 12 charge 0.4122 #Hff
set type 13 charge 0.0668 #Hgg
set type 14 charge 0.0668 #Hii
set type 15 charge 0.1216 #Hjj
set type 16 charge 0.1216 #Hkk
set type 17 charge 0.4120 #Hll
set type 18 charge 0.0669 #Hmm
set type 19 charge 0.0669 #Hnn
set type 20 charge 0.18 #Ha
set type 21 charge 0.18 #Hm
set type 22 charge 0.21 #Hr
set type 23 charge 0.06 #Hs
set type 24 charge 0.08 #Ht
set type 25 charge 0.27 #Hw
set type 26 charge -0.66 #N
set type 27 charge 0.22 # Na
set type 28 charge -0.53 #O
set type 29 charge -0.5766 #OOd
set type 30 charge -0.5769 # OOe
set type 31 charge 1.02 #S
variable TK equal 300.0
variable PBAR equal 1.0
kspace_style ewald 1.0e-4
min_style sd
velocity all create ${TK} 4567712 mom yes rot yes dist gaussian
minimize 1.0e-12 1.0e-12 100 1000
timestep 1.0
neighbor 2.0 bin
neigh_modify delay 5 check yes
thermo_style multi
thermo 50
fix TPSTAT all npt temp {TK} {TK} 10 iso {PBAR} {PBAR} 100
run 1000
reset_timestep 0
variable lambda equal ramp(0.0,1.0)
variable q2 equal 0.0083v_lambda
variable q3 equal -0.1981v_lambda
variable q4 equal 0.0085v_lambda
variable q12 equal 0.4122v_lambda
variable q13 equal 0.0668v_lambda
variable q14 equal 0.0668v_lambda
variable q15 equal 0.1216v_lambda
variable q16 equal 0.1216v_lambda
variable q17 equal 0.4120v_lambda
variable q18 equal 0.0669v_lambda
variable q19 equal 0.0669v_lambda
variable q29 equal -0.5766v_lambda
variable q30 equal -0.5769*v_lambda
fix ADAPT all adapt/fep 100 &
pair lj/cut/coul/long/soft lambda 12 34 v_lambda &
atom charge 2 v_q1 &
atom charge 3 v_q2 &
atom charge 4 v_q2 &
atom charge 12 v_q2 &
atom charge 13 v_q2 &
atom charge 14 v_q2 &
atom charge 15 v_q2 &
atom charge 16 v_q2 &
atom charge 17 v_q2 &
atom charge 18 v_q2 &
atom charge 19 v_q2 &
atom charge 29 v_q2 &
atom charge 30 v_q2 &
after yes
fix PRINT all print 100 “adapt lambda = {lambda} q2 = {q2} q3 = {q3} q4 = {q4} q12 = {q12} q13 = {q13} q14 = {q14} q15 = {q15} q16 = {q16} q17 = {q17} q18 = {q18} q19 = {q19} q29 = {q29} q30 = {q30}”
variable dlambda equal 0.05
variable dq2 equal 0.0083v_lambda
variable dq3 equal -0.1981v_lambda
variable dq4 equal 0.0085v_lambda
variable dq12 equal 0.4122v_lambda
variable dq13 equal 0.0668v_lambda
variable dq14 equal 0.0668v_lambda
variable dq15 equal 0.1216v_lambda
variable dq16 equal 0.1216v_lambda
variable dq17 equal 0.4120v_lambda
variable dq18 equal 0.0669v_lambda
variable dq19 equal 0.0669v_lambda
variable dq29 equal -0.5766v_lambda
variable dq30 equal -0.5769*v_lambda
compute fep all fep ${TK} &
pair lj/cut/coul/long/soft lambda 12 34 v_dlambda &
atom charge 2 v_dq2 &
atom charge 3 v_dq1 &
atom charge 4 v_dq1 &
atom charge 12 v_dq1 &
atom charge 13 v_dq1 &
atom charge 14 v_dq1 &
atom charge 15 v_dq2 &
atom charge 16 v_dq1 &
atom charge 17 v_dq1 &
atom charge 18 v_dq1 &
atom charge 19 v_dq1 &
atom charge 29 v_dq1 &
atom charge 30 v_dq1
fix fep all ave/time 2 40 100 c_fep[1] c_fep[2] file fep01.lmp
dump TRAJ all custom 20 dump.lammpstrj id mol type element x y z xu yu zu
run 2000
write_restart restart.*.lmp
write_data data.*.lmp