1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing authors: Ray Shan (Sandia, tnshan@sandia.gov)
17 ------------------------------------------------------------------------- */
18
19 #include "fix_qeq_comb.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "error.h"
24 #include "force.h"
25 #include "group.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "respa.h"
29 #include "update.h"
30
31 #include <cmath>
32 #include <cstring>
33
34 #include "pair_comb.h"
35 #include "pair_comb3.h"
36
37 using namespace LAMMPS_NS;
38 using namespace FixConst;
39
40 /* ---------------------------------------------------------------------- */
41
FixQEQComb(LAMMPS * lmp,int narg,char ** arg)42 FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
43 fp(nullptr), comb(nullptr), comb3(nullptr), qf(nullptr), q1(nullptr), q2(nullptr)
44 {
45 if (narg < 5) error->all(FLERR,"Illegal fix qeq/comb command");
46
47 peratom_flag = 1;
48 size_peratom_cols = 0;
49 peratom_freq = 1;
50 respa_level_support = 1;
51 ilevel_respa = 0;
52
53 nevery = utils::inumeric(FLERR,arg[3],false,lmp);
54 precision = utils::numeric(FLERR,arg[4],false,lmp);
55
56 if (nevery <= 0 || precision <= 0.0)
57 error->all(FLERR,"Illegal fix qeq/comb command");
58
59 MPI_Comm_rank(world,&me);
60
61 // optional args
62
63 int iarg = 5;
64 while (iarg < narg) {
65 if (strcmp(arg[iarg],"file") == 0) {
66 if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/comb command");
67 if (me == 0) {
68 fp = fopen(arg[iarg+1],"w");
69 if (fp == nullptr)
70 error->one(FLERR,std::string("Cannot open fix qeq/comb file ")
71 + arg[iarg+1]);
72 }
73 iarg += 2;
74 } else error->all(FLERR,"Illegal fix qeq/comb command");
75 }
76
77 nmax = atom->nmax;
78 memory->create(qf,nmax,"qeq:qf");
79 memory->create(q1,nmax,"qeq:q1");
80 memory->create(q2,nmax,"qeq:q2");
81 vector_atom = qf;
82
83 // zero the vector since dump may access it on timestep 0
84 // zero the vector since a variable may access it before first run
85
86 int nlocal = atom->nlocal;
87 for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
88
89 comm_forward = 1;
90 }
91
92 /* ---------------------------------------------------------------------- */
93
~FixQEQComb()94 FixQEQComb::~FixQEQComb()
95 {
96 if (me == 0 && fp) fclose(fp);
97 memory->destroy(qf);
98 memory->destroy(q1);
99 memory->destroy(q2);
100 }
101
102 /* ---------------------------------------------------------------------- */
103
setmask()104 int FixQEQComb::setmask()
105 {
106 int mask = 0;
107 mask |= POST_FORCE;
108 mask |= POST_FORCE_RESPA;
109 mask |= MIN_POST_FORCE;
110 return mask;
111 }
112
113 /* ---------------------------------------------------------------------- */
114
init()115 void FixQEQComb::init()
116 {
117 if (!atom->q_flag)
118 error->all(FLERR,"Fix qeq/comb requires atom attribute q");
119
120 comb3 = (PairComb3 *) force->pair_match("^comb3",0);
121 if (!comb3) comb = (PairComb *) force->pair_match("^comb",0);
122
123 if (comb == nullptr && comb3 == nullptr)
124 error->all(FLERR,"Must use pair_style comb or comb3 with fix qeq/comb");
125
126 if (utils::strmatch(update->integrate_style,"^respa")) {
127 ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
128 if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
129 }
130
131 ngroup = group->count(igroup);
132 if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
133 }
134
135 /* ---------------------------------------------------------------------- */
136
setup(int vflag)137 void FixQEQComb::setup(int vflag)
138 {
139 firstflag = 1;
140 if (utils::strmatch(update->integrate_style,"^verlet"))
141 post_force(vflag);
142 else {
143 ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
144 post_force_respa(vflag,ilevel_respa,0);
145 ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
146 }
147 firstflag = 0;
148 }
149
150 /* ---------------------------------------------------------------------- */
151
min_post_force(int vflag)152 void FixQEQComb::min_post_force(int vflag)
153 {
154 post_force(vflag);
155 }
156
157 /* ---------------------------------------------------------------------- */
158
post_force(int)159 void FixQEQComb::post_force(int /*vflag*/)
160 {
161 int i,ii,iloop,loopmax,inum,*ilist;
162 double heatpq,qmass,dtq,dtq2;
163 double enegchkall,enegmaxall;
164
165 if (update->ntimestep % nevery) return;
166
167 // reallocate work arrays if necessary
168 // qf = charge force
169 // q1 = charge displacement
170 // q2 = tmp storage of charge force for next iteration
171
172 if (atom->nmax > nmax) {
173 memory->destroy(qf);
174 memory->destroy(q1);
175 memory->destroy(q2);
176 nmax = atom->nmax;
177 memory->create(qf,nmax,"qeq:qf");
178 memory->create(q1,nmax,"qeq:q1");
179 memory->create(q2,nmax,"qeq:q2");
180 vector_atom = qf;
181 }
182
183 // more loops for first-time charge equilibrium
184
185 iloop = 0;
186 if (firstflag) loopmax = 200;
187 else loopmax = 100;
188
189 // charge-equilibration loop
190
191 if (me == 0 && fp)
192 fmt::print(fp,"Charge equilibration on step {}\n", update->ntimestep);
193
194 heatpq = 0.05;
195 qmass = 0.016;
196 dtq = 0.01;
197 dtq2 = 0.5*dtq*dtq/qmass;
198
199 double enegchk = 0.0;
200 double enegtot = 0.0;
201 double enegmax = 0.0;
202
203 double *q = atom->q;
204 int *mask = atom->mask;
205
206 if (comb) {
207 inum = comb->list->inum;
208 ilist = comb->list->ilist;
209 } else if (comb3) {
210 inum = comb3->list->inum;
211 ilist = comb3->list->ilist;
212 } else {
213 inum = 0;
214 ilist = nullptr;
215 }
216
217 for (ii = 0; ii < inum; ii++) {
218 i = ilist[ii];
219 q1[i] = q2[i] = qf[i] = 0.0;
220 }
221
222 for (iloop = 0; iloop < loopmax; iloop ++) {
223 for (ii = 0; ii < inum; ii++) {
224 i = ilist[ii];
225 if (mask[i] & groupbit) {
226 q1[i] += qf[i]*dtq2 - heatpq*q1[i];
227 q[i] += q1[i];
228 }
229 }
230
231 comm->forward_comm_fix(this);
232 enegtot = 0.0;
233 if (comb) enegtot = comb->yasu_char(qf,igroup);
234 else if (comb3) enegtot = comb3->combqeq(qf,igroup);
235
236 enegtot /= ngroup;
237 enegchk = enegmax = 0.0;
238
239 for (ii = 0; ii < inum ; ii++) {
240 i = ilist[ii];
241 if (mask[i] & groupbit) {
242 q2[i] = enegtot-qf[i];
243 enegmax = MAX(enegmax,fabs(q2[i]));
244 enegchk += fabs(q2[i]);
245 qf[i] = q2[i];
246 }
247 }
248
249 MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
250 enegchk = enegchkall/ngroup;
251 MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
252 enegmax = enegmaxall;
253
254 if (enegchk <= precision && enegmax <= 100.0*precision) break;
255
256 if (me == 0 && fp)
257 fmt::print(fp," iteration: {}, enegtot {:.6g}, "
258 "enegmax {:.6g}, fq deviation: {:.6g}\n",
259 iloop,enegtot,enegmax,enegchk);
260
261 for (ii = 0; ii < inum; ii++) {
262 i = ilist[ii];
263 if (mask[i] & groupbit)
264 q1[i] += qf[i]*dtq2 - heatpq*q1[i];
265 }
266 }
267
268 if (me == 0 && fp) {
269 if (iloop == loopmax)
270 fmt::print(fp,"Charges did not converge in {} iterations\n",iloop);
271 else
272 fmt::print(fp,"Charges converged in {} iterations to {:.10f} tolerance\n",
273 iloop,enegchk);
274 }
275 }
276
277 /* ---------------------------------------------------------------------- */
278
post_force_respa(int vflag,int ilevel,int)279 void FixQEQComb::post_force_respa(int vflag, int ilevel, int /*iloop*/)
280 {
281 if (ilevel == ilevel_respa) post_force(vflag);
282 }
283
284 /* ----------------------------------------------------------------------
285 memory usage of local atom-based arrays
286 ------------------------------------------------------------------------- */
287
memory_usage()288 double FixQEQComb::memory_usage()
289 {
290 double bytes = (double)atom->nmax*3 * sizeof(double);
291 return bytes;
292 }
293 /* ---------------------------------------------------------------------- */
294
pack_forward_comm(int n,int * list,double * buf,int,int *)295 int FixQEQComb::pack_forward_comm(int n, int *list, double *buf,
296 int /*pbc_flag*/, int * /*pbc*/)
297 {
298 int i,j,m;
299
300 m = 0;
301 for (i = 0; i < n; i ++) {
302 j = list[i];
303 buf[m++] = atom->q[j];
304 }
305 return m;
306 }
307
308 /* ---------------------------------------------------------------------- */
309
unpack_forward_comm(int n,int first,double * buf)310 void FixQEQComb::unpack_forward_comm(int n, int first, double *buf)
311 {
312 int i,m,last;
313
314 m = 0;
315 last = first + n ;
316 for (i = first; i < last; i++) atom->q[i] = buf[m++];
317 }
318