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