Questions on fix qeq/point code structure

Hello, I had some questions on the net charge constraint and code structure of fix qeq.

  1. Is the net charge constraint set to be zero, or to the same net charge seen in the first timestep? I attempted to figure this out by printing the net charge (via compute reduce) at every step. It is fluctuating, but at a very low value, so I am not sure if it is constrained to zero or the initial value.

Step c_qNet
0 1.058615e-15
1 -1.1878519e-15
2 -7.775898e-16
3 -1.477117e-15
4 -9.1940344e-17
5 -4.4235449e-17
6 1.8782718e-15
7 1.9081958e-17
8 2.0543463e-15
9 -3.6385825e-16
10 2.2846308e-15

  1. I am also trying to discern the structure of the code better, to understand how the constraints are applied.

I looked at the matrix structure in the original qeq paper (Rappe, Goddard 1991) and slide 8 of a workshop talk ( https://lammps.sandia.gov/workshops/Aug15/PDF/talk_Shan.pdf ). Slide 8 gives a layout of the matrix structure and where the next charge constraint is located. However, when looking at init_matvec() in fix_qeq_point.cpp, there appear to be some intermediate variables, and I am not entirely sure of their meaning or why the matrix solution steps are structured the way they are.

In init_matvec(), there are the variables b, t, b_s, b_t, t_hist, s_hist, which I am not sure of the meaning of. Also, how does the net charge constraint show up the matrix structure?

for( ii = 0; ii < inum; ++ii ) {
i = ilist[ii];
if (atom->mask[i] & groupbit) {
Hdia_inv[i] = 1. / eta[ atom->type[i] ];
b_s[i] = -( chi[atom->type[i]] + chizj[i] );
b_t[i] = -1.0;
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] );
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
}
}

In pre_force(), I am also not sure why the CG solver is run in 2 steps on variables b_s, s, b_t, and t:

matvecs = CG(b_s, s); // CG on s - parallel
matvecs += CG(b_t, t); // CG on t - parallel

Any pointers are greatly appreciated.

Thank you,
Anne

Hello, I had some questions on the net charge constraint and code structure of fix qeq.

  1. Is the net charge constraint set to be zero, or to the same net charge seen in the first timestep? I attempted to figure this out by printing the net charge (via compute reduce) at every step. It is fluctuating, but at a very low value, so I am not sure if it is constrained to zero or the initial value.

these numbers are effectively zero. 1e-15 is close to the limit of precision of double precision floating point numbers. if you are summing up floating point numbers of different magnitude that should add up to an integer number, you will almost never get exactly that number.
please have a look at the following explanations of the mysteries and pitfalls of floating point math to understand why:
https://floating-point-gui.de/
http://blog.reverberate.org/2014/09/what-every-computer-programmer-should.html
http://blog.reverberate.org/2016/02/06/floating-point-demystified-part2.html

Ray, who wrote the QEQ package, may have some comments on your second question.

axel.

Thank you Axel,

I also have an additional question for Ray:

  1. Regarding the H matrix in the code: Is the H matrix in the code the same as the C matrix in CD=-D (from the paper and slides), and what is the purpose of having a separate Hdia_inv vector?