Problems with fix/qEq

Anthony_DeFilippo · July 30, 2015, 4:37pm

Hello,
I am responding to an old thread, but I include it below, as it appears
relevant to my case.

Late last year I developed what seemed to be a workable set of qeq
parameters for the qeq/point model with the 23-Nov-2014 LAMMPS version.
When I just recently went back to re-run some things with an updated LAMMPS
version, the partial charges returned by the model were rather different.
Some detective work brought me to this email thread and I found that the
January 30th patch (http://lammps.sandia.gov/patches/patch.30Jan15) was
what updated the qeq_point code so that now the qeq/point model is 'more
correct,' as evidenced by the fact that the results of the
examples/qeq/in.qeq.buck problem now are the same when using qeq/dynamic
and qeq/point.

My question is this: Since my calculations and model tuning were
unfortunately made with an outdated version of the code, is there perhaps a
logical conversion from parameters developed using the old (pre 30 January
2015) qeq/point model to work with the updated LAMMPS? As a simplified
example of my question, for example, is there a logical modification that
could be made to the parameters in examples/qeq/param.qeq2 such that
running the script examples/qeq/in.qeq.buck with fix qeq/point selected
would produce the same results as were produced by the original parameters
using a pre-30Jan2015 version of LAMMPS? It seems that dividing all qeq
parameters by 2 and running the new code makes the results closer to what
they were before, but they are not exactly the same, so perhaps that is
just a coincidence. I am pasting the results from the before and after code
as well as the qeq parameters below the message. (I of course don't expect
any trial-and-error on anyone's part, just hoping for some insight into
what actually changed about the code)

It appears that the relevant changes had to do with how and under what
conditions LAMMPS updates certain values in the "H matrix," and there is a
factor of 2 difference in one of the quantities between the two codes,
though I have trouble interpreting the actual implications of such
changes. If my parameters cannot be converted to work with the new code, my
question then becomes, would it make any sense to continue running an old
version of the code, or is the old qeq/point model simply unreliable?

Thanks, and thanks to all for all of your work on LAMMPS.

Sincerely,
Anthony DeFilippo

Output from examples/qeq/in.qeq.buck using qeq/point with 28 Jul 2015
version of code and original parameters
Step PotEng q1 q2 qtot
       0 -3432.1801 0.85228288 -0.42614144 -2.8421709e-13
       1 -3432.1845 0.85228293 -0.42614147 -2.8421709e-13
       2 -3432.2576 0.85228264 -0.42614132 -5.6843419e-14
       3 -3432.5072 0.85229828 -0.42614914 3.9790393e-13
       4 -3435.0114 0.85263921 -0.42631961 2.2737368e-13
       5 -3439.3347 0.8532404 -0.4266202 -2.2737368e-13
       6 -3444.3866 0.8539392 -0.4269696 1.7053026e-13
       7 -3448.1162 0.85442995 -0.42721497 1.1368684e-13
       8 -3449.5089 0.85456235 -0.42728117 3.4106051e-13
       9 -3450.6008 0.85463914 -0.42731957 1.7053026e-13

Output from examples/qeq/in.qeq.buck using qeq/point with 23 Nov 2014
version of code
Step PotEng q1 q2 qtot
       0 -4203.3633 0.96050687 -0.48025343 -1.7053026e-13
       1 -4203.3665 0.96050718 -0.48025359 1.1368684e-13
       2 -4203.4316 0.96050724 -0.48025362 2.8421709e-13
       3 -4203.9046 0.96055215 -0.48027607 -5.6843419e-13
       4 -4208.8342 0.9611769 -0.48058845 -4.5474735e-13
       5 -4218.2838 0.96239138 -0.48119569 -2.8421709e-13
       6 -4228.2874 0.96366808 -0.48183404 -2.2737368e-13
       7 -4235.6522 0.96457429 -0.48228714 -4.5474735e-13
       8 -4238.2575 0.96485121 -0.4824256 3.9790393e-13
       9 -4239.3066 0.96492333 -0.48246167 -1.1368684e-13
      10 -4242.2467 0.96522386 -0.48261193 -1.7053026e-13

original parameters param.qeq2
1 0.00000 7.25028 0.01 0.772871 0.000000
2 11.26882 15.37920 0.01 0.243072 0.000000

Output from examples/qeq/in.qeq.buck using qeq/point with 28 Jul 2015
version of code and all parameters in param.qeq2 divided by 2:

Step PotEng q1 q2 qtot
       0 -4349.5182 0.97955931 -0.48977965 -2.8421709e-13
       1 -4349.5206 0.97955946 -0.48977973 6.8212103e-13
       2 -4349.5783 0.97955878 -0.48977939 1.7053026e-13
       3 -4350.026 0.97960074 -0.48980037 0
       4 -4356.9253 0.98046878 -0.49023439 -2.8421709e-13
       5 -4369.953 0.98211939 -0.49105969 2.2737368e-13
       6 -4383.7305 0.98385791 -0.49192895 6.8212103e-13
       7 -4395.4175 0.98531242 -0.49265621 -6.2527761e-13
       8 -4398.0506 0.98559377 -0.49279688 -5.6843419e-14
       9 -4400.2942 0.98581556 -0.49290778 3.4106051e-13
      10 -4403.0957 0.9861013 -0.49305065 -2.8421709e-13

Updated QEQ parameters with qeq/point relevant parameters divided by 2:
1 0.00000 3.62514 0.01 0.772871 0.000000
2 5.634410 7.689600 0.01 0.243072 0.000000

Ray_Shan · August 9, 2015, 5:35pm

Hi Anthony,

Thanks for the report. The 1/30 version did not make changes to qeq, perhaps it is the 1/19 version? But the 1/19 patch would only affect if you have fix npt or fix box/relax that changes box dimensions. Can you send me your full input deck so I can take a look at it? If it is work in progress for piblication, you can send me your qeq parameters privately.

Thanks,
Ray

Anthony_DeFilippo · August 10, 2015, 4:49pm

Hi Ray, thanks for offering to look at this. I will send you some things privately when I get a chance to compile them into a better-commented format, and hopefully we can come up with an answer that can be generalized to benefit the LAMMPS community.

It very well could have been the 1/19 fix that made the change, as I am not familiar enough with the syntax of changelogs such as this to say that this is certainly when the change occurred. It just looks to me like the http://lammps.sandia.gov/patches/patch.30Jan15 link shows a change to the file fix_qeq_point.cpp:

diff -Naur lammps-21Jan15/src/QEQ/fix_qeq_point.cpp lammps-30Jan15/src/QEQ/fix_qeq_point.cpp
— lammps-21Jan15/src/QEQ/fix_qeq_point.cpp 2015-01-19 16:49:20.000000000 -0700
+++ lammps-30Jan15/src/QEQ/fix_qeq_point.cpp 2015-01-30 15:22:31.573389000 -0700
…

if( flag ) {

if (r_sqr <= cutoff_sq) {
H.jlist[m_fill] = j;
r = sqrt(r_sqr);

H.val[m_fill] = 1.0/r;

H.val[m_fill] = 0.5/r;
m_fill++;

Either way, I tested the stable LAMMPS releases surrounding January 2014 (lammps-9Dec14 and lammps-10Feb2015) and they return different results when running the in.qeq.buck example with qeq/point selected. That example case runs nve, so perhaps the difference in tutorial results indicates that there was also some other change to qeq that affects more than just npt cases. The case I will send you is run in npt, so perhaps there have been two relevant changes that I will have to deal with.

Thanks again,
Anthony