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 Variant of the harmonic angle potential for use with the
19 lj/sdk potential for coarse grained MD simulations.
20 ------------------------------------------------------------------------- */
21
22 #include "angle_sdk.h"
23
24 #include <cmath>
25 #include "atom.h"
26 #include "neighbor.h"
27 #include "pair.h"
28 #include "domain.h"
29 #include "comm.h"
30 #include "force.h"
31 #include "math_const.h"
32 #include "memory.h"
33 #include "error.h"
34
35
36 #include "lj_sdk_common.h"
37
38 using namespace LAMMPS_NS;
39 using namespace MathConst;
40 using namespace LJSDKParms;
41
42 #define SMALL 0.001
43
44 /* ---------------------------------------------------------------------- */
45
AngleSDK(LAMMPS * lmp)46 AngleSDK::AngleSDK(LAMMPS *lmp) : Angle(lmp) { repflag = 0;}
47
48 /* ---------------------------------------------------------------------- */
49
~AngleSDK()50 AngleSDK::~AngleSDK()
51 {
52 if (allocated) {
53 memory->destroy(setflag);
54 memory->destroy(k);
55 memory->destroy(theta0);
56 memory->destroy(repscale);
57
58 allocated = 0;
59 }
60 }
61
62 /* ---------------------------------------------------------------------- */
63
compute(int eflag,int vflag)64 void AngleSDK::compute(int eflag, int vflag)
65 {
66 int i1,i2,i3,n,type;
67 double delx1,dely1,delz1,delx2,dely2,delz2,delx3,dely3,delz3;
68 double eangle,f1[3],f3[3],e13,f13;
69 double dtheta,tk;
70 double rsq1,rsq2,rsq3,r1,r2,c,s,a,a11,a12,a22;
71
72 eangle = 0.0;
73 ev_init(eflag,vflag);
74
75 double **x = atom->x;
76 double **f = atom->f;
77 int **anglelist = neighbor->anglelist;
78 int nanglelist = neighbor->nanglelist;
79 int nlocal = atom->nlocal;
80 int newton_bond = force->newton_bond;
81
82 for (n = 0; n < nanglelist; n++) {
83 i1 = anglelist[n][0];
84 i2 = anglelist[n][1];
85 i3 = anglelist[n][2];
86 type = anglelist[n][3];
87
88 // 1st bond
89
90 delx1 = x[i1][0] - x[i2][0];
91 dely1 = x[i1][1] - x[i2][1];
92 delz1 = x[i1][2] - x[i2][2];
93
94 rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
95 r1 = sqrt(rsq1);
96
97 // 2nd bond
98
99 delx2 = x[i3][0] - x[i2][0];
100 dely2 = x[i3][1] - x[i2][1];
101 delz2 = x[i3][2] - x[i2][2];
102
103 rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
104 r2 = sqrt(rsq2);
105
106 // angle (cos and sin)
107
108 c = delx1*delx2 + dely1*dely2 + delz1*delz2;
109 c /= r1*r2;
110
111 if (c > 1.0) c = 1.0;
112 if (c < -1.0) c = -1.0;
113
114 s = sqrt(1.0 - c*c);
115 if (s < SMALL) s = SMALL;
116 s = 1.0/s;
117
118 // 1-3 LJ interaction.
119 // we only want to use the repulsive part,
120 // and it can be scaled (or off).
121 // so this has to be done here and not in the
122 // general non-bonded code.
123
124 f13 = e13 = delx3 = dely3 = delz3 = 0.0;
125
126 if (repflag) {
127
128 delx3 = x[i1][0] - x[i3][0];
129 dely3 = x[i1][1] - x[i3][1];
130 delz3 = x[i1][2] - x[i3][2];
131 rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3;
132
133 const int type1 = atom->type[i1];
134 const int type3 = atom->type[i3];
135
136 f13=0.0;
137 e13=0.0;
138
139 if (rsq3 < rminsq[type1][type3]) {
140 const int ljt = lj_type[type1][type3];
141 const double r2inv = 1.0/rsq3;
142
143 if (ljt == LJ12_4) {
144 const double r4inv=r2inv*r2inv;
145
146 f13 = r4inv*(lj1[type1][type3]*r4inv*r4inv - lj2[type1][type3]);
147 if (eflag) e13 = r4inv*(lj3[type1][type3]*r4inv*r4inv - lj4[type1][type3]);
148
149 } else if (ljt == LJ9_6) {
150 const double r3inv = r2inv*sqrt(r2inv);
151 const double r6inv = r3inv*r3inv;
152
153 f13 = r6inv*(lj1[type1][type3]*r3inv - lj2[type1][type3]);
154 if (eflag) e13 = r6inv*(lj3[type1][type3]*r3inv - lj4[type1][type3]);
155
156 } else if (ljt == LJ12_6) {
157 const double r6inv = r2inv*r2inv*r2inv;
158
159 f13 = r6inv*(lj1[type1][type3]*r6inv - lj2[type1][type3]);
160 if (eflag) e13 = r6inv*(lj3[type1][type3]*r6inv - lj4[type1][type3]);
161 }
162
163 // make sure energy is 0.0 at the cutoff.
164 if (eflag) e13 -= emin[type1][type3];
165
166 f13 *= r2inv;
167 }
168 }
169
170 // force & energy
171
172 dtheta = acos(c) - theta0[type];
173 tk = k[type] * dtheta;
174
175 if (eflag) eangle = tk*dtheta;
176
177 a = -2.0 * tk * s;
178 a11 = a*c / rsq1;
179 a12 = -a / (r1*r2);
180 a22 = a*c / rsq2;
181
182 f1[0] = a11*delx1 + a12*delx2;
183 f1[1] = a11*dely1 + a12*dely2;
184 f1[2] = a11*delz1 + a12*delz2;
185 f3[0] = a22*delx2 + a12*delx1;
186 f3[1] = a22*dely2 + a12*dely1;
187 f3[2] = a22*delz2 + a12*delz1;
188
189 // apply force to each of the 3 atoms
190
191 if (newton_bond || i1 < nlocal) {
192 f[i1][0] += f1[0] + f13*delx3;
193 f[i1][1] += f1[1] + f13*dely3;
194 f[i1][2] += f1[2] + f13*delz3;
195 }
196
197 if (newton_bond || i2 < nlocal) {
198 f[i2][0] -= f1[0] + f3[0];
199 f[i2][1] -= f1[1] + f3[1];
200 f[i2][2] -= f1[2] + f3[2];
201 }
202
203 if (newton_bond || i3 < nlocal) {
204 f[i3][0] += f3[0] - f13*delx3;
205 f[i3][1] += f3[1] - f13*dely3;
206 f[i3][2] += f3[2] - f13*delz3;
207 }
208
209 if (evflag) {
210 ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
211 delx1,dely1,delz1,delx2,dely2,delz2);
212 if (repflag)
213 ev_tally13(i1,i3,nlocal,newton_bond,e13,f13,delx3,dely3,delz3);
214 }
215 }
216 }
217
218 /* ---------------------------------------------------------------------- */
219
allocate()220 void AngleSDK::allocate()
221 {
222 allocated = 1;
223 int n = atom->nangletypes;
224
225 memory->create(k,n+1,"angle:k");
226 memory->create(theta0,n+1,"angle:theta0");
227 memory->create(repscale,n+1,"angle:repscale");
228
229 memory->create(setflag,n+1,"angle:setflag");
230 for (int i = 1; i <= n; i++) setflag[i] = 0;
231 }
232
233 /* ----------------------------------------------------------------------
234 set coeffs for one or more types
235 ------------------------------------------------------------------------- */
236
coeff(int narg,char ** arg)237 void AngleSDK::coeff(int narg, char **arg)
238 {
239 if ((narg < 3) || (narg > 6))
240 error->all(FLERR,"Incorrect args for angle coefficients");
241
242 if (!allocated) allocate();
243
244 int ilo,ihi;
245 utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
246
247 double k_one = utils::numeric(FLERR,arg[1],false,lmp);
248 double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
249 double repscale_one=1.0;
250
251 // backward compatibility with old cg/cmm style input:
252 // this had <lj_type> <epsilon> <sigma>
253 // if epsilon is set to 0.0 we accept it as repscale 0.0
254 // otherwise assume repscale 1.0, since we were using
255 // epsilon to turn repulsion on or off.
256 if (narg == 6) {
257 repscale_one = utils::numeric(FLERR,arg[4],false,lmp);
258 if (repscale_one > 0.0) repscale_one = 1.0;
259 } else if (narg == 4) repscale_one = utils::numeric(FLERR,arg[3],false,lmp);
260 else if (narg == 3) repscale_one = 1.0;
261 else error->all(FLERR,"Incorrect args for angle coefficients");
262
263 // convert theta0 from degrees to radians and store coefficients
264
265 int count = 0;
266 for (int i = ilo; i <= ihi; i++) {
267 k[i] = k_one;
268 theta0[i] = theta0_one/180.0 * MY_PI;
269 repscale[i] = repscale_one;
270 setflag[i] = 1;
271 count++;
272 }
273
274 if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
275 }
276
277 /* ----------------------------------------------------------------------
278 error check and initialize all values needed for force computation
279 ------------------------------------------------------------------------- */
280
init_style()281 void AngleSDK::init_style()
282 {
283
284 // make sure we use an SDK pair_style and that we need the 1-3 repulsion
285
286 repflag = 0;
287 for (int i = 1; i <= atom->nangletypes; i++)
288 if (repscale[i] > 0.0) repflag = 1;
289
290 // set up pointers to access SDK LJ parameters for 1-3 interactions
291
292 if (repflag) {
293 int itmp;
294 if (force->pair == nullptr)
295 error->all(FLERR,"Angle style SDK requires use of a compatible with Pair style");
296
297 lj1 = (double **) force->pair->extract("lj1",itmp);
298 lj2 = (double **) force->pair->extract("lj2",itmp);
299 lj3 = (double **) force->pair->extract("lj3",itmp);
300 lj4 = (double **) force->pair->extract("lj4",itmp);
301 lj_type = (int **) force->pair->extract("lj_type",itmp);
302 rminsq = (double **) force->pair->extract("rminsq",itmp);
303 emin = (double **) force->pair->extract("emin",itmp);
304
305 if (!lj1 || !lj2 || !lj3 || !lj4 || !lj_type || !rminsq || !emin)
306 error->all(FLERR,"Angle style SDK is incompatible with Pair style");
307 }
308 }
309
310 /* ---------------------------------------------------------------------- */
311
equilibrium_angle(int i)312 double AngleSDK::equilibrium_angle(int i)
313 {
314 return theta0[i];
315 }
316
317 /* ----------------------------------------------------------------------
318 proc 0 writes out coeffs to restart file
319 ------------------------------------------------------------------------- */
320
write_restart(FILE * fp)321 void AngleSDK::write_restart(FILE *fp)
322 {
323 fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
324 fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
325 fwrite(&repscale[1],sizeof(double),atom->nangletypes,fp);
326 }
327
328 /* ----------------------------------------------------------------------
329 proc 0 reads coeffs from restart file, bcasts them
330 ------------------------------------------------------------------------- */
331
read_restart(FILE * fp)332 void AngleSDK::read_restart(FILE *fp)
333 {
334 allocate();
335
336 if (comm->me == 0) {
337 utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
338 utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
339 utils::sfread(FLERR,&repscale[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
340 }
341 MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
342 MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
343 MPI_Bcast(&repscale[1],atom->nangletypes,MPI_DOUBLE,0,world);
344
345 for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
346 }
347
348 /* ----------------------------------------------------------------------
349 proc 0 writes to data file
350 ------------------------------------------------------------------------- */
351
write_data(FILE * fp)352 void AngleSDK::write_data(FILE *fp)
353 {
354 for (int i = 1; i <= atom->nangletypes; i++)
355 fprintf(fp,"%d %g %g\n",i,k[i],theta0[i]/MY_PI*180.0);
356 }
357
358 /* ---------------------------------------------------------------------- */
359
ev_tally13(int i,int j,int nlocal,int newton_bond,double evdwl,double fpair,double delx,double dely,double delz)360 void AngleSDK::ev_tally13(int i, int j, int nlocal, int newton_bond,
361 double evdwl, double fpair,
362 double delx, double dely, double delz)
363 {
364 double v[6];
365
366 if (eflag_either) {
367 if (eflag_global) {
368 if (newton_bond) {
369 energy += evdwl;
370 } else {
371 if (i < nlocal)
372 energy += 0.5*evdwl;
373 if (j < nlocal)
374 energy += 0.5*evdwl;
375 }
376 }
377 if (eflag_atom) {
378 if (newton_bond || i < nlocal) eatom[i] += 0.5*evdwl;
379 if (newton_bond || j < nlocal) eatom[j] += 0.5*evdwl;
380 }
381 }
382
383 if (vflag_either) {
384 v[0] = delx*delx*fpair;
385 v[1] = dely*dely*fpair;
386 v[2] = delz*delz*fpair;
387 v[3] = delx*dely*fpair;
388 v[4] = delx*delz*fpair;
389 v[5] = dely*delz*fpair;
390
391 if (vflag_global) {
392 if (newton_bond) {
393 virial[0] += v[0];
394 virial[1] += v[1];
395 virial[2] += v[2];
396 virial[3] += v[3];
397 virial[4] += v[4];
398 virial[5] += v[5];
399 } else {
400 if (i < nlocal) {
401 virial[0] += 0.5*v[0];
402 virial[1] += 0.5*v[1];
403 virial[2] += 0.5*v[2];
404 virial[3] += 0.5*v[3];
405 virial[4] += 0.5*v[4];
406 virial[5] += 0.5*v[5];
407 }
408 if (j < nlocal) {
409 virial[0] += 0.5*v[0];
410 virial[1] += 0.5*v[1];
411 virial[2] += 0.5*v[2];
412 virial[3] += 0.5*v[3];
413 virial[4] += 0.5*v[4];
414 virial[5] += 0.5*v[5];
415 }
416 }
417 }
418
419 if (vflag_atom) {
420 if (newton_bond || i < nlocal) {
421 vatom[i][0] += 0.5*v[0];
422 vatom[i][1] += 0.5*v[1];
423 vatom[i][2] += 0.5*v[2];
424 vatom[i][3] += 0.5*v[3];
425 vatom[i][4] += 0.5*v[4];
426 vatom[i][5] += 0.5*v[5];
427 }
428 if (newton_bond || j < nlocal) {
429 vatom[j][0] += 0.5*v[0];
430 vatom[j][1] += 0.5*v[1];
431 vatom[j][2] += 0.5*v[2];
432 vatom[j][3] += 0.5*v[3];
433 vatom[j][4] += 0.5*v[4];
434 vatom[j][5] += 0.5*v[5];
435 }
436 }
437 }
438 }
439
440 /* ---------------------------------------------------------------------- */
441
single(int type,int i1,int i2,int i3)442 double AngleSDK::single(int type, int i1, int i2, int i3)
443 {
444 double **x = atom->x;
445
446 double delx1 = x[i1][0] - x[i2][0];
447 double dely1 = x[i1][1] - x[i2][1];
448 double delz1 = x[i1][2] - x[i2][2];
449 domain->minimum_image(delx1,dely1,delz1);
450 double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
451
452 double delx2 = x[i3][0] - x[i2][0];
453 double dely2 = x[i3][1] - x[i2][1];
454 double delz2 = x[i3][2] - x[i2][2];
455 domain->minimum_image(delx2,dely2,delz2);
456 double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
457
458 double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
459 c /= r1*r2;
460 if (c > 1.0) c = 1.0;
461 if (c < -1.0) c = -1.0;
462
463 double e13=0.0;
464 if (repflag) {
465
466 // 1-3 LJ interaction.
467 double delx3 = x[i1][0] - x[i3][0];
468 double dely3 = x[i1][1] - x[i3][1];
469 double delz3 = x[i1][2] - x[i3][2];
470 domain->minimum_image(delx3,dely3,delz3);
471
472 const int type1 = atom->type[i1];
473 const int type3 = atom->type[i3];
474
475 const double rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3;
476
477 if (rsq3 < rminsq[type1][type3]) {
478 const int ljt = lj_type[type1][type3];
479 const double r2inv = 1.0/rsq3;
480
481 if (ljt == LJ12_4) {
482 const double r4inv=r2inv*r2inv;
483
484 e13 = r4inv*(lj3[type1][type3]*r4inv*r4inv - lj4[type1][type3]);
485
486 } else if (ljt == LJ9_6) {
487 const double r3inv = r2inv*sqrt(r2inv);
488 const double r6inv = r3inv*r3inv;
489
490 e13 = r6inv*(lj3[type1][type3]*r3inv - lj4[type1][type3]);
491
492 } else if (ljt == LJ12_6) {
493 const double r6inv = r2inv*r2inv*r2inv;
494
495 e13 = r6inv*(lj3[type1][type3]*r6inv - lj4[type1][type3]);
496 }
497
498 // make sure energy is 0.0 at the cutoff.
499 e13 -= emin[type1][type3];
500 }
501 }
502
503 double dtheta = acos(c) - theta0[type];
504 double tk = k[type] * dtheta;
505 return tk*dtheta + e13;
506 }
507