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:
17 James Fischer, High Performance Technologies, Inc.
18 Charles Cornwell, High Performance Technologies, Inc.
19 David Richie, Stone Ridge Technology
20 Vincent Natoli, Stone Ridge Technology
21 ------------------------------------------------------------------------- */
22
23 #include "pair_eam_opt.h"
24 #include <cmath>
25
26 #include "atom.h"
27 #include "comm.h"
28 #include "force.h"
29 #include "neigh_list.h"
30 #include "memory.h"
31 #include "update.h"
32
33 using namespace LAMMPS_NS;
34
35 /* ---------------------------------------------------------------------- */
36
PairEAMOpt(LAMMPS * lmp)37 PairEAMOpt::PairEAMOpt(LAMMPS *lmp) : PairEAM(lmp) {}
38
39 /* ---------------------------------------------------------------------- */
40
compute(int eflag,int vflag)41 void PairEAMOpt::compute(int eflag, int vflag)
42 {
43 ev_init(eflag,vflag);
44
45 if (evflag) {
46 if (eflag) {
47 if (force->newton_pair) return eval<1,1,1>();
48 else return eval<1,1,0>();
49 } else {
50 if (force->newton_pair) return eval<1,0,1>();
51 else return eval<1,0,0>();
52 }
53 } else {
54 if (force->newton_pair) return eval<0,0,1>();
55 else return eval<0,0,0>();
56 }
57 }
58
59 /* ---------------------------------------------------------------------- */
60
61 template < int EVFLAG, int EFLAG, int NEWTON_PAIR >
eval()62 void PairEAMOpt::eval()
63 {
64 typedef struct { double x,y,z; } vec3_t;
65
66 typedef struct {
67 double rhor0i,rhor1i,rhor2i,rhor3i;
68 double rhor0j,rhor1j,rhor2j,rhor3j;
69 } fast_alpha_t;
70
71 typedef struct {
72 double rhor4i,rhor5i,rhor6i;
73 double rhor4j,rhor5j,rhor6j;
74 double z2r0,z2r1,z2r2,z2r3,z2r4,z2r5,z2r6;
75 double _pad[3];
76 } fast_gamma_t;
77
78 int i,j,ii,jj,inum,jnum,itype,jtype;
79 double evdwl = 0.0;
80 double* _noalias coeff;
81
82 // grow energy array if necessary
83
84 if (atom->nmax > nmax) {
85 memory->destroy(rho);
86 memory->destroy(fp);
87 memory->destroy(numforce);
88 nmax = atom->nmax;
89 memory->create(rho,nmax,"pair:rho");
90 memory->create(fp,nmax,"pair:fp");
91 memory->create(numforce,nmax,"pair:numforce");
92 }
93
94 double** _noalias x = atom->x;
95 double** _noalias f = atom->f;
96 int* _noalias type = atom->type;
97 int nlocal = atom->nlocal;
98
99 vec3_t* _noalias xx = (vec3_t*)x[0];
100 vec3_t* _noalias ff = (vec3_t*)f[0];
101
102 double tmp_cutforcesq = cutforcesq;
103 double tmp_rdr = rdr;
104 int nr2 = nr-2;
105 int nr1 = nr-1;
106
107 inum = list->inum;
108 int* _noalias ilist = list->ilist;
109 int** _noalias firstneigh = list->firstneigh;
110 int* _noalias numneigh = list->numneigh;
111
112 int ntypes = atom->ntypes;
113 int ntypes2 = ntypes*ntypes;
114
115 fast_alpha_t* _noalias fast_alpha =
116 (fast_alpha_t*) malloc((size_t)ntypes2*(nr+1)*sizeof(fast_alpha_t));
117 for (i = 0; i < ntypes; i++) for (j = 0; j < ntypes; j++) {
118 fast_alpha_t* _noalias tab = &fast_alpha[i*ntypes*nr+j*nr];
119 if (type2rhor[i+1][j+1] >= 0) {
120 for (int m = 1; m <= nr; m++) {
121 tab[m].rhor0i = rhor_spline[type2rhor[i+1][j+1]][m][6];
122 tab[m].rhor1i = rhor_spline[type2rhor[i+1][j+1]][m][5];
123 tab[m].rhor2i = rhor_spline[type2rhor[i+1][j+1]][m][4];
124 tab[m].rhor3i = rhor_spline[type2rhor[i+1][j+1]][m][3];
125 }
126 }
127 if (type2rhor[j+1][i+1] >= 0) {
128 for (int m = 1; m <= nr; m++) {
129 tab[m].rhor0j = rhor_spline[type2rhor[j+1][i+1]][m][6];
130 tab[m].rhor1j = rhor_spline[type2rhor[j+1][i+1]][m][5];
131 tab[m].rhor2j = rhor_spline[type2rhor[j+1][i+1]][m][4];
132 tab[m].rhor3j = rhor_spline[type2rhor[j+1][i+1]][m][3];
133 }
134 }
135 }
136 fast_alpha_t* _noalias tabeight = fast_alpha;
137
138 fast_gamma_t* _noalias fast_gamma =
139 (fast_gamma_t*) malloc((size_t)ntypes2*(nr+1)*sizeof(fast_gamma_t));
140 for (i = 0; i < ntypes; i++) for (j = 0; j < ntypes; j++) {
141 fast_gamma_t* _noalias tab = &fast_gamma[i*ntypes*nr+j*nr];
142 if (type2rhor[i+1][j+1] >= 0) {
143 for (int m = 1; m <= nr; m++) {
144 tab[m].rhor4i = rhor_spline[type2rhor[i+1][j+1]][m][2];
145 tab[m].rhor5i = rhor_spline[type2rhor[i+1][j+1]][m][1];
146 tab[m].rhor6i = rhor_spline[type2rhor[i+1][j+1]][m][0];
147 }
148 }
149 if (type2rhor[j+1][i+1] >= 0) {
150 for (int m = 1; m <= nr; m++) {
151 tab[m].rhor4j = rhor_spline[type2rhor[j+1][i+1]][m][2];
152 tab[m].rhor5j = rhor_spline[type2rhor[j+1][i+1]][m][1];
153 tab[m].rhor6j = rhor_spline[type2rhor[j+1][i+1]][m][0];
154 tab[m].z2r6 = z2r_spline[type2z2r[i+1][j+1]][m][0];
155 }
156 }
157 if (type2z2r[i+1][j+1] >= 0) {
158 for (int m = 1; m <= nr; m++) {
159 tab[m].z2r0 = z2r_spline[type2z2r[i+1][j+1]][m][6];
160 tab[m].z2r1 = z2r_spline[type2z2r[i+1][j+1]][m][5];
161 tab[m].z2r2 = z2r_spline[type2z2r[i+1][j+1]][m][4];
162 tab[m].z2r3 = z2r_spline[type2z2r[i+1][j+1]][m][3];
163 tab[m].z2r4 = z2r_spline[type2z2r[i+1][j+1]][m][2];
164 tab[m].z2r5 = z2r_spline[type2z2r[i+1][j+1]][m][1];
165 tab[m].z2r6 = z2r_spline[type2z2r[i+1][j+1]][m][0];
166 }
167 }
168 }
169 fast_gamma_t* _noalias tabss = fast_gamma;
170
171 // zero out density
172
173 if (NEWTON_PAIR) {
174 int m = nlocal + atom->nghost;
175 for (i = 0; i < m; i++) rho[i] = 0.0;
176 } else for (i = 0; i < nlocal; i++) rho[i] = 0.0;
177
178 // rho = density at each atom
179 // loop over neighbors of my atoms
180
181 for (ii = 0; ii < inum; ii++) {
182 i = ilist[ii];
183 double xtmp = xx[i].x;
184 double ytmp = xx[i].y;
185 double ztmp = xx[i].z;
186 itype = type[i] - 1;
187 int* _noalias jlist = firstneigh[i];
188 jnum = numneigh[i];
189
190 double tmprho = rho[i];
191 fast_alpha_t* _noalias tabeighti = &tabeight[itype*ntypes*nr];
192
193 for (jj = 0; jj < jnum; jj++) {
194 j = jlist[jj];
195 j &= NEIGHMASK;
196
197 double delx = xtmp - xx[j].x;
198 double dely = ytmp - xx[j].y;
199 double delz = ztmp - xx[j].z;
200 double rsq = delx*delx + dely*dely + delz*delz;
201
202 if (rsq < tmp_cutforcesq) {
203 jtype = type[j] - 1;
204
205 double p = sqrt(rsq)*tmp_rdr;
206 if ((int)p <= nr2) {
207 int m = (int)p + 1;
208 p -= (double)((int)p);
209 fast_alpha_t& a = tabeighti[jtype*nr+m];
210 tmprho += ((a.rhor3j*p+a.rhor2j)*p+a.rhor1j)*p+a.rhor0j;
211 if (NEWTON_PAIR || j < nlocal) {
212 rho[j] += ((a.rhor3i*p+a.rhor2i)*p+a.rhor1i)*p+a.rhor0i;
213 }
214 } else {
215 fast_alpha_t& a = tabeighti[jtype*nr+nr1];
216 tmprho += a.rhor3j+a.rhor2j+a.rhor1j+a.rhor0j;
217 if (NEWTON_PAIR || j < nlocal) {
218 rho[j] += a.rhor3i+a.rhor2i+a.rhor1i+a.rhor0i;
219 }
220 }
221 }
222 }
223 rho[i] = tmprho;
224 }
225
226 // communicate and sum densities
227
228 if (NEWTON_PAIR) comm->reverse_comm_pair(this);
229
230 // fp = derivative of embedding energy at each atom
231 // phi = embedding energy at each atom
232 // if rho > rhomax (e.g. due to close approach of two atoms),
233 // will exceed table, so add linear term to conserve energy
234
235 for (ii = 0; ii < inum; ii++) {
236 i = ilist[ii];
237 double p = rho[i]*rdrho + 1.0;
238 int m = static_cast<int> (p);
239 m = MAX(1,MIN(m,nrho-1));
240 p -= m;
241 p = MIN(p,1.0);
242 coeff = frho_spline[type2frho[type[i]]][m];
243 fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2];
244 if (EFLAG) {
245 double phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
246 if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
247 phi *= scale[type[i]][type[i]];
248 if (eflag_global) eng_vdwl += phi;
249 if (eflag_atom) eatom[i] += phi;
250 }
251 }
252
253 // communicate derivative of embedding function
254
255 comm->forward_comm_pair(this);
256 embedstep = update->ntimestep;
257
258 // compute forces on each atom
259 // loop over neighbors of my atoms
260
261 for (ii = 0; ii < inum; ii++) {
262 i = ilist[ii];
263 double xtmp = xx[i].x;
264 double ytmp = xx[i].y;
265 double ztmp = xx[i].z;
266 int itype1 = type[i] - 1;
267 int* _noalias jlist = firstneigh[i];
268 jnum = numneigh[i];
269
270 double tmpfx = 0.0;
271 double tmpfy = 0.0;
272 double tmpfz = 0.0;
273
274 fast_gamma_t* _noalias tabssi = &tabss[itype1*ntypes*nr];
275 double* _noalias scale_i = scale[itype1+1]+1;
276 numforce[i] = 0;
277
278 for (jj = 0; jj < jnum; jj++) {
279 j = jlist[jj];
280 j &= NEIGHMASK;
281
282 double delx = xtmp - xx[j].x;
283 double dely = ytmp - xx[j].y;
284 double delz = ztmp - xx[j].z;
285 double rsq = delx*delx + dely*dely + delz*delz;
286
287 if (rsq < tmp_cutforcesq) {
288 ++numforce[i];
289 jtype = type[j] - 1;
290 double r = sqrt(rsq);
291 double rhoip,rhojp,z2,z2p;
292 double p = r*tmp_rdr;
293 if ((int)p <= nr2) {
294 int m = (int) p + 1;
295 m = MIN(m,nr-1);
296 p -= (double)((int) p);
297 p = MIN(p,1.0);
298
299 fast_gamma_t& a = tabssi[jtype*nr+m];
300 rhoip = (a.rhor6i*p + a.rhor5i)*p + a.rhor4i;
301 rhojp = (a.rhor6j*p + a.rhor5j)*p + a.rhor4j;
302 z2 = ((a.z2r3*p + a.z2r2)*p + a.z2r1)*p + a.z2r0;
303 z2p = (a.z2r6*p + a.z2r5)*p + a.z2r4;
304
305 } else {
306
307 fast_gamma_t& a = tabssi[jtype*nr+nr1];
308 rhoip = a.rhor6i + a.rhor5i + a.rhor4i;
309 rhojp = a.rhor6j + a.rhor5j + a.rhor4j;
310 z2 = a.z2r3 + a.z2r2 + a.z2r1 + a.z2r0;
311 z2p = a.z2r6 + a.z2r5 + a.z2r4;
312 }
313
314 // rhoip = derivative of (density at atom j due to atom i)
315 // rhojp = derivative of (density at atom i due to atom j)
316 // phi = pair potential energy
317 // phip = phi'
318 // z2 = phi * r
319 // z2p = (phi * r)' = (phi' r) + phi
320 // psip needs both fp[i] and fp[j] terms since r_ij appears in two
321 // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
322 // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
323 // scale factor can be applied by thermodynamic integration
324
325 double recip = 1.0/r;
326 double phi = z2*recip;
327 double phip = z2p*recip - phi*recip;
328 double psip = fp[i]*rhojp + fp[j]*rhoip + phip;
329 double fpair = -scale_i[jtype]*psip*recip;
330
331 tmpfx += delx*fpair;
332 tmpfy += dely*fpair;
333 tmpfz += delz*fpair;
334 if (NEWTON_PAIR || j < nlocal) {
335 ff[j].x -= delx*fpair;
336 ff[j].y -= dely*fpair;
337 ff[j].z -= delz*fpair;
338 }
339
340 if (EFLAG) evdwl = scale_i[jtype]*phi;
341
342 if (EVFLAG) ev_tally(i,j,nlocal,NEWTON_PAIR,
343 evdwl,0.0,fpair,delx,dely,delz);
344 }
345 }
346
347 ff[i].x += tmpfx;
348 ff[i].y += tmpfy;
349 ff[i].z += tmpfz;
350 }
351
352 free(fast_alpha); fast_alpha = 0;
353 free(fast_gamma); fast_gamma = 0;
354
355 if (vflag_fdotr) virial_fdotr_compute();
356 }
357