I am simulating water with a flexible SPC/E-like model which has an intermolecular LJ potential for O-O, O-H, and H-H. The potential specifies an intermolecular potential "U_intermolecular = U_lj + U_couloumb", and an intramolecular potential "U_intramolecular = 1/2 K_theta (theta - theta_0)^2 + 1/2 K_r (r - r_0)^2". With this in mind, should I be using the "special_bonds" command to prevent 'double-counting' the lj and coulomb pair energy. Additionally, isn't this also the case if we have a long-range TIP4P, SPC/E, or TIP3P potential for water?
Thanks in advance,
My input file (with special_bonds) follows:
units real
dimension 3
boundary p p p
atom_style full
read_data initial.water.dat
group ox type 2
group hy type 1
set group ox charge -0.8476
set group hy charge 0.4238
### Flexible SPC/E Potential Parameters ###
### Zhang et al., Fluid Phase Equilibria, 262 (2007) 210-216 ###
pair_style lj/cut/coul/long 10.0
pair_coeff 1 1 0.6287 3.1169
pair_coeff 1 2 0.142723 2.04845
pair_coeff 2 2 0.0324 0.98
bond_style harmonic
bond_coeff 1 1480 0.9611
angle_style harmonic
angle_coeff 1 353 109.4712
kspace_style pppm 1.0e-5 #final npt relaxation
special_bonds lj/coul 0.0 0.0 0.0
velocity all create 298.0 1234546 dist gaussian
#fix 2 all npt temp 298.0 298.0 100.0 iso 1.0 1.0 1000.0
#Compute
compute Oke ox ke
compute totalke all ke
compute totalpe all pe
#Variables
variable totale equal (c_totalke+c_totalpe)/count(ox) #total energy per H2O molecule
variable keperO equal c_Oke/count(ox) #average total energy of oxygen atoms
#Computes
#Compute averages and save to file
fix 3 all ave/time 1 1000 1000 c_thermo_temp c_thermo_press c_totalke c_totalpe v_totale v_keperO file ave.spc.dat
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
thermo 100
thermo_style custom step temp etotal press
run
timestep 0.5
run 10