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

- 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

- 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