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 author: Axel Kohlmeyer (Temple U)
17 ------------------------------------------------------------------------- */
18
19 #include "fix_qeq_comb_omp.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 "pair_comb.h"
29 #include "respa.h"
30 #include "update.h"
31
32 #include <cmath>
33
34 using namespace LAMMPS_NS;
35 using namespace FixConst;
36
37 /* ---------------------------------------------------------------------- */
38
FixQEQCombOMP(LAMMPS * lmp,int narg,char ** arg)39 FixQEQCombOMP::FixQEQCombOMP(LAMMPS *lmp, int narg, char **arg) :
40 FixQEQComb(lmp, narg, arg)
41 {
42 if (narg < 5) error->all(FLERR,"Illegal fix qeq/comb/omp command");
43 }
44
45 /* ---------------------------------------------------------------------- */
46
init()47 void FixQEQCombOMP::init()
48 {
49 if (!atom->q_flag)
50 error->all(FLERR,"Fix qeq/comb/omp requires atom attribute q");
51
52 if (nullptr != force->pair_match("comb3",0))
53 error->all(FLERR,"No support for comb3 currently available in OPENMP");
54
55 comb = (PairComb *) force->pair_match("comb/omp",1);
56 if (comb == nullptr)
57 comb = (PairComb *) force->pair_match("comb",1);
58 if (comb == nullptr)
59 error->all(FLERR,"Must use pair_style comb or "
60 "comb/omp with fix qeq/comb/omp");
61
62 if (utils::strmatch(update->integrate_style,"^respa")) {
63 ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
64 if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
65 }
66
67 ngroup = group->count(igroup);
68 if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms");
69 }
70
71 /* ---------------------------------------------------------------------- */
72
post_force(int)73 void FixQEQCombOMP::post_force(int /* vflag */)
74 {
75 int i,ii,iloop,loopmax,inum,*ilist;
76 double heatpq,qmass,dtq,dtq2;
77 double enegchkall,enegmaxall;
78
79 if (update->ntimestep % nevery) return;
80
81 // reallocate work arrays if necessary
82 // qf = charge force
83 // q1 = charge displacement
84 // q2 = tmp storage of charge force for next iteration
85
86 if (atom->nmax > nmax) {
87 memory->destroy(qf);
88 memory->destroy(q1);
89 memory->destroy(q2);
90 nmax = atom->nmax;
91 memory->create(qf,nmax,"qeq:qf");
92 memory->create(q1,nmax,"qeq:q1");
93 memory->create(q2,nmax,"qeq:q2");
94 vector_atom = qf;
95 }
96
97 // more loops for first-time charge equilibrium
98
99 iloop = 0;
100 if (firstflag) loopmax = 500;
101 else loopmax = 200;
102
103 // charge-equilibration loop
104
105 if (me == 0 && fp)
106 fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
107 update->ntimestep);
108
109 heatpq = 0.05;
110 qmass = 0.016;
111 dtq = 0.01;
112 dtq2 = 0.5*dtq*dtq/qmass;
113
114 double enegchk = 0.0;
115 double enegtot = 0.0;
116 double enegmax = 0.0;
117
118 double *q = atom->q;
119 int *mask = atom->mask;
120
121 inum = comb->list->inum;
122 ilist = comb->list->ilist;
123
124 for (ii = 0; ii < inum; ii++) {
125 i = ilist[ii];
126 q1[i] = q2[i] = qf[i] = 0.0;
127 }
128
129 for (iloop = 0; iloop < loopmax; iloop ++) {
130 for (ii = 0; ii < inum; ii++) {
131 i = ilist[ii];
132 if (mask[i] & groupbit) {
133 q1[i] += qf[i]*dtq2 - heatpq*q1[i];
134 q[i] += q1[i];
135 }
136 }
137 comm->forward_comm_fix(this);
138
139 if (comb) enegtot = comb->yasu_char(qf,igroup);
140 enegtot /= ngroup;
141 enegchk = enegmax = 0.0;
142
143 for (ii = 0; ii < inum ; ii++) {
144 i = ilist[ii];
145 if (mask[i] & groupbit) {
146 q2[i] = enegtot-qf[i];
147 enegmax = MAX(enegmax,fabs(q2[i]));
148 enegchk += fabs(q2[i]);
149 qf[i] = q2[i];
150 }
151 }
152
153 MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
154 enegchk = enegchkall/ngroup;
155 MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
156 enegmax = enegmaxall;
157
158 if (enegchk <= precision && enegmax <= 100.0*precision) break;
159
160 if (me == 0 && fp)
161 fprintf(fp," iteration: %d, enegtot %.6g, "
162 "enegmax %.6g, fq deviation: %.6g\n",
163 iloop,enegtot,enegmax,enegchk);
164
165 for (ii = 0; ii < inum; ii++) {
166 i = ilist[ii];
167 if (mask[i] & groupbit)
168 q1[i] += qf[i]*dtq2 - heatpq*q1[i];
169 }
170 }
171
172 if (me == 0 && fp) {
173 if (iloop == loopmax)
174 fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
175 else
176 fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
177 iloop,enegchk);
178 }
179 }
180