LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Issue with fix qeq/point for SPC/E water simulation
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Issue with fix qeq/point for SPC/E water simulation


From: Ray Shan <rshan@...1795...>
Date: Tue, 22 Aug 2017 20:07:56 +0000

You can’t have a cutoff that is only 3 angstroms; at least 10 angstroms or the same as the Coulomb cutoff is reasonable.

Ray

On 8/21/17, 2:40 PM, "Baig abdullah Al muhit via lammps-users" <lammps-users@lists.sourceforge.net> wrote:

Hello everyone,
I am simulating SPC/E water model. I assigned charges with fix qeq method instead of set command. I observed several issues I am not understanding.
1. I used fix qeq/point to assign charge to 1 molecule of water without invoking any pair potential for run 0. The charge on O and H were -0.2 and 0.1, respectively for this simulation. And, after NVE dynamics run with shake, the potential energy was -13 Kcal/mol. The parameter files are from Goddard's (1991) http://pubs.acs.org/doi/pdf/10.1021/j100161a070 paper,

param file:
1 8.741 13.364 0 0 0
2 4.528 13.89 0 0 0

2. I used fix qeq/point to assign charge on O and H. But, this time I used chi and eta parameters from ffield.reax.CHO file. I multiplied eta by 2. The charge on O and H were -0.125 and 0.0625, respectively, and the potential energy was -5 Kcal/mol.

param file (ffield.reax.CHO):
1 8.5 17.9978 0 0 0
2 5.32 14.87 0 0 0

3. I used set command to assign charges O=-0.8476 and H=0.4238. After NVE dynamics run with shake, the potential energy was -211 Kcal/mol.

According to Goddard's paper, charge on H atom is 0.35, even then it's well below SPC/E model's charge. Am I applying fix qeq incorrectly? Potential energy is very low in method (1) and (2) compared to method 3. If we take, SPC/E model's charge as reference point, the first two methods are clearing giving me wrong answers.

Here's my input file:
##
units         real
atom_style     full
boundary    p p p
newton         on
dimension    3

read_data     wat.data

fix        1 all qeq/point 1 3 1.0e-6 100 param.qeq1
atom_modify    sort 0 10
group        ox type 1
group        hyd type 2
#set        type 1 charge -0.8476
#set        type 2 charge 0.4238
variable    no equal count(ox)
variable    nh equal count(hyd)

compute        qO ox property/atom q
compute        oxcharge ox reduce sum c_qO
variable    qO equal c_oxcharge/v_no

compute        qH hyd property/atom q
compute        hydcharge hyd reduce sum c_qH
variable    qH equal c_hydcharge/v_nh
thermo        100
thermo_style     custom step temp v_qO v_qH pe
run        0
unfix        1

pair_style    lj/cut/coul/long 10 10
bond_style    harmonic
angle_style    harmonic
dihedral_style  none
improper_style  none

pair_coeff     1 1 0.1553 3.166
pair_coeff     2 2 0 0

bond_coeff     1 540.6336 1.0
angle_coeff    1 50 109.47
pair_modify    mix arithmetic tail no
kspace_style    pppm 1e-6
special_bonds     lj/coul 1 1 1
neighbor    3.0 bin
neigh_modify     delay 0 every 1 one 10000 check yes
    compute    wattemp all temp
    velocity     all create 1 31216 temp wattemp dist gaussian
    minimize     1.0e-20 1.0e-20 500000 1000000
    min_style     cg
    min_modify     dmax 0.1
fix        2 all nve
fix        3 all shake 1e-4 20 0 b 1 a 1
fix        4 all temp/berendsen 1 1 100
thermo        100
thermo_style     custom step temp v_qO v_qH pe
run        10000

Regards,
Baig