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 #include "pair_coul_exclude.h"
16
17 #include "atom.h"
18 #include "comm.h"
19 #include "error.h"
20 #include "force.h"
21 #include "memory.h"
22 #include "neigh_list.h"
23 #include "neighbor.h"
24
25 #include <cmath>
26 #include <cstring>
27
28 using namespace LAMMPS_NS;
29
30 /* ---------------------------------------------------------------------- */
31
PairCoulExclude(LAMMPS * lmp)32 PairCoulExclude::PairCoulExclude(LAMMPS *lmp) : Pair(lmp) {
33 writedata = 1;
34 }
35
36 /* ---------------------------------------------------------------------- */
37
~PairCoulExclude()38 PairCoulExclude::~PairCoulExclude()
39 {
40 if (copymode) return;
41
42 if (allocated) {
43 memory->destroy(setflag);
44 memory->destroy(cutsq);
45 }
46 }
47
48 /* ---------------------------------------------------------------------- */
49
compute(int eflag,int vflag)50 void PairCoulExclude::compute(int eflag, int vflag)
51 {
52 int i,j,ii,jj,inum,jnum;
53 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
54 double rsq,r2inv,rinv,forcecoul,factor_coul;
55 int *ilist,*jlist,*numneigh,**firstneigh;
56
57 ecoul = 0.0;
58 ev_init(eflag,vflag);
59
60 double **x = atom->x;
61 double **f = atom->f;
62 double *q = atom->q;
63 int nlocal = atom->nlocal;
64 double *special_coul = force->special_coul;
65 int newton_pair = force->newton_pair;
66 double qqrd2e = force->qqrd2e;
67
68 inum = list->inum;
69 ilist = list->ilist;
70 numneigh = list->numneigh;
71 firstneigh = list->firstneigh;
72
73 // loop over neighbors of my atoms
74
75 for (ii = 0; ii < inum; ii++) {
76 i = ilist[ii];
77 qtmp = q[i];
78 xtmp = x[i][0];
79 ytmp = x[i][1];
80 ztmp = x[i][2];
81 jlist = firstneigh[i];
82 jnum = numneigh[i];
83
84 for (jj = 0; jj < jnum; jj++) {
85 j = jlist[jj];
86
87 // skip over all non-excluded pairs and subtract excluded coulomb
88
89 if (sbmask(j) == 0) continue;
90 factor_coul = special_coul[sbmask(j)] - 1.0;
91
92 j &= NEIGHMASK;
93
94 delx = xtmp - x[j][0];
95 dely = ytmp - x[j][1];
96 delz = ztmp - x[j][2];
97 rsq = delx*delx + dely*dely + delz*delz;
98
99 r2inv = 1.0/rsq;
100 rinv = sqrt(r2inv);
101 forcecoul = qqrd2e * qtmp*q[j]*rinv;
102 fpair = factor_coul*forcecoul * r2inv;
103
104 f[i][0] += delx*fpair;
105 f[i][1] += dely*fpair;
106 f[i][2] += delz*fpair;
107 if (newton_pair || j < nlocal) {
108 f[j][0] -= delx*fpair;
109 f[j][1] -= dely*fpair;
110 f[j][2] -= delz*fpair;
111 }
112
113 if (eflag)
114 ecoul = factor_coul * qqrd2e * qtmp*q[j]*rinv;
115
116 if (evflag) ev_tally(i,j,nlocal,newton_pair,0.0,ecoul,fpair,delx,dely,delz);
117 }
118 }
119
120 if (vflag_fdotr) virial_fdotr_compute();
121 }
122
123 /* ----------------------------------------------------------------------
124 allocate all arrays
125 ------------------------------------------------------------------------- */
126
allocate()127 void PairCoulExclude::allocate()
128 {
129 allocated = 1;
130 int n = atom->ntypes;
131
132 memory->create(setflag,n+1,n+1,"pair:setflag");
133 memory->create(cutsq,n+1,n+1,"pair:setflag");
134 for (int i = 1; i <= n; i++)
135 for (int j = i; j <= n; j++)
136 setflag[i][j] = 0;
137 }
138
139 /* ----------------------------------------------------------------------
140 global settings
141 ------------------------------------------------------------------------- */
142
settings(int narg,char ** arg)143 void PairCoulExclude::settings(int narg, char **arg)
144 {
145 if (narg != 1) error->all(FLERR,"Illegal pair_style command");
146
147 cut_global = utils::numeric(FLERR,arg[0],false,lmp);
148 }
149
150 /* ----------------------------------------------------------------------
151 set coeffs for one or more type pairs
152 ------------------------------------------------------------------------- */
153
coeff(int narg,char ** arg)154 void PairCoulExclude::coeff(int narg, char **arg)
155 {
156 if (narg != 2)
157 error->all(FLERR,"Incorrect args for pair coefficients");
158 if (!allocated) allocate();
159
160 int ilo,ihi,jlo,jhi;
161 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
162 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
163
164 int count = 0;
165 for (int i = ilo; i <= ihi; i++) {
166 for (int j = MAX(jlo,i); j <= jhi; j++) {
167 setflag[i][j] = 1;
168 count++;
169 }
170 }
171
172 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
173 }
174
175
176 /* ----------------------------------------------------------------------
177 init specific to this pair style
178 ------------------------------------------------------------------------- */
179
init_style()180 void PairCoulExclude::init_style()
181 {
182 if (!atom->q_flag)
183 error->all(FLERR,"Pair style coul/exclude requires atom attribute q");
184
185 neighbor->request(this,instance_me);
186 }
187
188 /* ----------------------------------------------------------------------
189 init for one type pair i,j and corresponding j,i
190 ------------------------------------------------------------------------- */
191
init_one(int,int)192 double PairCoulExclude::init_one(int /*i*/, int /*j*/)
193 {
194 return cut_global;
195 }
196
197 /* ----------------------------------------------------------------------
198 proc 0 writes to restart file
199 ------------------------------------------------------------------------- */
200
write_restart(FILE * fp)201 void PairCoulExclude::write_restart(FILE *fp)
202 {
203 write_restart_settings(fp);
204
205 int i,j;
206 for (i = 1; i <= atom->ntypes; i++) {
207 for (j = i; j <= atom->ntypes; j++) {
208 fwrite(&setflag[i][j],sizeof(int),1,fp);
209 }
210 }
211 }
212
213 /* ----------------------------------------------------------------------
214 proc 0 reads from restart file, bcasts
215 ------------------------------------------------------------------------- */
216
read_restart(FILE * fp)217 void PairCoulExclude::read_restart(FILE *fp)
218 {
219 read_restart_settings(fp);
220 allocate();
221
222 int i,j;
223 int me = comm->me;
224 for (i = 1; i <= atom->ntypes; i++) {
225 for (j = i; j <= atom->ntypes; j++) {
226 if (me == 0) {
227 utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
228 }
229 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
230 }
231 }
232 }
233
234 /* ----------------------------------------------------------------------
235 proc 0 writes to restart file
236 ------------------------------------------------------------------------- */
237
write_restart_settings(FILE * fp)238 void PairCoulExclude::write_restart_settings(FILE *fp)
239 {
240 fwrite(&cut_global,sizeof(double),1,fp);
241 fwrite(&offset_flag,sizeof(int),1,fp);
242 fwrite(&mix_flag,sizeof(int),1,fp);
243 }
244
245 /* ----------------------------------------------------------------------
246 proc 0 reads from restart file, bcasts
247 ------------------------------------------------------------------------- */
248
read_restart_settings(FILE * fp)249 void PairCoulExclude::read_restart_settings(FILE *fp)
250 {
251 if (comm->me == 0) {
252 utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
253 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
254 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
255 }
256 MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
257 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
258 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
259 }
260
261 /* ----------------------------------------------------------------------
262 proc 0 writes to data file
263 ------------------------------------------------------------------------- */
264
write_data(FILE * fp)265 void PairCoulExclude::write_data(FILE *fp)
266 {
267 for (int i = 1; i <= atom->ntypes; i++)
268 fprintf(fp,"%d\n",i);
269 }
270
271 /* ----------------------------------------------------------------------
272 proc 0 writes all pairs to data file
273 ------------------------------------------------------------------------- */
274
write_data_all(FILE * fp)275 void PairCoulExclude::write_data_all(FILE *fp)
276 {
277 for (int i = 1; i <= atom->ntypes; i++)
278 for (int j = i; j <= atom->ntypes; j++)
279 fprintf(fp,"%d %d\n",i,j);
280 }
281
282 /* ---------------------------------------------------------------------- */
283
single(int i,int j,int,int,double rsq,double factor_coul,double,double & fforce)284 double PairCoulExclude::single(int i, int j, int /*itype*/, int /*jtype*/,
285 double rsq, double factor_coul, double /*factor_lj*/,
286 double &fforce)
287 {
288 double r2inv,rinv,forcecoul,phicoul;
289
290 factor_coul -= 1.0; // we want to subtract the excluded coulomb
291 r2inv = 1.0/rsq;
292 rinv = sqrt(r2inv);
293 forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*rinv;
294 fforce = factor_coul*forcecoul * r2inv;
295
296 phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*rinv;
297 return factor_coul*phicoul;
298 }
299
300 /* ---------------------------------------------------------------------- */
301
extract(const char * str,int & dim)302 void *PairCoulExclude::extract(const char *str, int &dim)
303 {
304 dim = 0;
305 if (strcmp(str,"cut_coul") == 0) return (void *) &cut_global;
306 return nullptr;
307 }
308