fluctuating charge water model

Hello all:
Thanks to Ray's hard work for implementing the reax-independent QEq routines. I've been having some issues running a water box of SPC/Fq waters however. I'm attempting to follow to simulations of Rick and Stuart [J. Chem. Phys. 101, 6141 (1994)].

I always get a segmentation fault when using the qeq/point, qeq/dynamic and qeq/shielded options. The qeq/slater option runs, however there are convergence issues and very nonphysical charges/electrostatic after 1 or 2 steps. I've minimized the system without qeq first, so I'm sure the starting structure is reasonable. The QEq parameters are from the Stuart paper.

I should point out that this is only a problem if I use the long range electrostatics solvers. If I switch to a truncated coulomb pair style (say lj/charmm/coul/charmm), the system is stable. I wonder if the problem is in evaluating the self energy term in the ewald routines. For fixed charges, I believe that this is evaluated once at the start of the simulation and is taken as a constant thereafter. However, for fluctuating charges, it would appear to me that this term would have to be reevaluated every step. I'm not well versed enough with the LAMMPS source to know where to look (of course this has probably been thought of before and already implemented).

Anyway, I currently stuck. I know that this is a relatively new feature in lammps, but if someone has any experience/advice in running these sort of systems, it would be highly appreciated.

Best,

Mike

param1.qeq (71 Bytes)

in.spcfq.box (2.19 KB)

data.spcfq.box (81 KB)

Hello all:

Thanks to Ray’s hard work for implementing the reax-independent QEq routines. I’ve been having some issues running a water box of SPC/Fq waters however. I’m attempting to follow to simulations of Rick and Stuart [J. Chem. Phys. 101, 6141 (1994)].

I always get a segmentation fault when using the qeq/point, qeq/dynamic and qeq/shielded options. The qeq/slater option runs, however there are convergence issues and very nonphysical charges/electrostatic after 1 or 2 steps. I’ve minimized the system without qeq first, so I’m sure the starting structure is reasonable. The QEq parameters are from the Stuart paper.

I should point out that this is only a problem if I use the long range electrostatics solvers. If I switch to a truncated coulomb pair style (say lj/charmm/coul/charmm), the system is stable.

I wonder if the problem is in evaluating the self energy term in the ewald routines. For fixed charges, I believe that this is evaluated once at the start of the simulation and is taken as a constant thereafter. However, for fluctuating charges, it would appear to me that this term would have to be reevaluated every step. I’m not well versed enough with the LAMMPS source to know where to look (of course this has probably been thought of before and already implemented).

Anyway, I currently stuck. I know that this is a relatively new feature in lammps, but if someone has any experience/advice in running these sort of systems, it would be highly appreciated.

Best,

Mike

param1.qeq (71 Bytes)

in.spcfq.box (2.19 KB)

data.spcfq.box (81 KB)

Hi Michael,

Michael,

Thanks to Axel’s help, the segmentation fault is now fixed. Please update your source code with the attached files, thanks.

I have also made a QEq parameter file from Stuart’s paper. Somehow I had to change X parameters to get published charges.

fix_qeq_dynamic.cpp (7.21 KB)

fix_qeq_point.cpp (4.94 KB)

fix_qeq_shielded.cpp (6.74 KB)

fix_qeq.h (3.21 KB)

fix_qeq.cpp (17.6 KB)

I am not sure if this works for that paper. In that paper the charges in each molecule is constrained to be zero. I don’t think this constrain is implemented in qeq/variant.

Best,
Yi

You are correct that intermolecular charge transfer is not constrained in fix qeq/variants, but the entire system is constrained to be neutral. In Stuart’s paper, both QEq formulation are provided and fix qeq/variants just implement the more general (and more physical) case - intermolecular charge transfer allowed so that each molecule can carry non-zero charge.

Fix qeq/slater has the Stuart’s formulation, it just solves charge using the matrix inversion method, as opposed to the extended Lagragian method used by Stuart et al.

What is your reasoning that fix qeq/variants do not work for Stuart’s paper?

Ray

I agree all the interaction are included in qeq/variant. And I think it is a great work.

However, in that paper, most of the results are got without intermolecular charge transfer. They mentioned that with intermolecular charge transfer the total system would be metallic and will cause problems.

Yi

Yes, you are right and that is a known problem for QEq. However, if you have a condensed phase (solid or liquid), then the problem identified by Stuart et al in their paper (when molecules are separated in vacuo) does not exist. Therefore, it is still reasonable to use QEq for liquid water, which is what Michael is interested in.

Thanks,
Ray

Yes, I agree that it works for liquid water. In that paper, bulk water systems are also calculated without intermolecular charge transfer. So, I think it if Michael is trying to get the same bulk water properties with that paper it may failed. If he is trying to use the parameters, I think it works.

Best,
Yi

These files will be in the next patch.

Thanks, Ray

Steve

Hi Ray:
Thanks for this (and sorry for not getting back to you until now).
May I ask:
How did you come up with the chi (electronegativity) parameters? Trial and error? In the Stuart paper, they only report chi(oxygen) - chi(hydrogen) = 73.33 kcal/mol. How did you come up with the values you reported?

Thanks
–Mike

I have lost track. The values I reported where?

In ReaxFF and COMB papers, these QEq parameters are fitted to DFT/experiments so that they reproduce the ionization energies and/or Mulliken (or Bader) charges.

Ray

Hi Ray:
Yes it has been a while so no worries. I would just curious about your statement below in your Sept 14th email last year:

"Michael,

Thanks to Axel’s help, the segmentation fault is now fixed. Please update your source code with the attached files, thanks.

I have also made a QEq parameter file from Stuart’s paper. Somehow I had to change X parameters to get published charges.