1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 /* ----------------------------------------------------------------------
15 Contributing author: Mike Brown (SNL)
16 ------------------------------------------------------------------------- */
17
18 #include "pair_resquared.h"
19
20 #include "atom.h"
21 #include "atom_vec_ellipsoid.h"
22 #include "comm.h"
23 #include "error.h"
24 #include "force.h"
25 #include "math_extra.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "neighbor.h"
29
30 #include <cmath>
31
32 using namespace LAMMPS_NS;
33
34 /* ---------------------------------------------------------------------- */
35
PairRESquared(LAMMPS * lmp)36 PairRESquared::PairRESquared(LAMMPS *lmp) :
37 Pair(lmp), cr60(pow(60.0, 1.0 / 3.0)), b_alpha(45.0 / 56.0)
38 {
39 single_enable = 0;
40
41 cr60 = pow(60.0, 1.0 / 3.0);
42 b_alpha = 45.0 / 56.0;
43 solv_f_a = 3.0 / (16.0 * atan(1.0) * -36.0);
44 solv_f_r = 3.0 / (16.0 * atan(1.0) * 2025.0);
45 }
46
47 /* ----------------------------------------------------------------------
48 free all arrays
49 ------------------------------------------------------------------------- */
50
~PairRESquared()51 PairRESquared::~PairRESquared()
52 {
53 if (allocated) {
54 memory->destroy(setflag);
55 memory->destroy(cutsq);
56
57 memory->destroy(form);
58 memory->destroy(epsilon);
59 memory->destroy(sigma);
60 memory->destroy(shape1);
61 memory->destroy(shape2);
62 memory->destroy(well);
63 memory->destroy(cut);
64 memory->destroy(lj1);
65 memory->destroy(lj2);
66 memory->destroy(lj3);
67 memory->destroy(lj4);
68 memory->destroy(offset);
69 delete[] lshape;
70 delete[] setwell;
71 }
72 }
73
74 /* ---------------------------------------------------------------------- */
75
compute(int eflag,int vflag)76 void PairRESquared::compute(int eflag, int vflag)
77 {
78 int i, j, ii, jj, inum, jnum, itype, jtype;
79 double evdwl, one_eng, rsq, r2inv, r6inv, forcelj, factor_lj;
80 double fforce[3], ttor[3], rtor[3], r12[3];
81 int *ilist, *jlist, *numneigh, **firstneigh;
82 RE2Vars wi, wj;
83
84 evdwl = 0.0;
85 ev_init(eflag, vflag);
86
87 double **x = atom->x;
88 double **f = atom->f;
89 double **tor = atom->torque;
90 int *type = atom->type;
91 int nlocal = atom->nlocal;
92 double *special_lj = force->special_lj;
93 int newton_pair = force->newton_pair;
94
95 inum = list->inum;
96 ilist = list->ilist;
97 numneigh = list->numneigh;
98 firstneigh = list->firstneigh;
99
100 // loop over neighbors of my atoms
101
102 for (ii = 0; ii < inum; ii++) {
103 i = ilist[ii];
104 itype = type[i];
105
106 // not a LJ sphere
107
108 if (lshape[itype] != 0.0) precompute_i(i, wi);
109
110 jlist = firstneigh[i];
111 jnum = numneigh[i];
112
113 for (jj = 0; jj < jnum; jj++) {
114 j = jlist[jj];
115 factor_lj = special_lj[sbmask(j)];
116 j &= NEIGHMASK;
117
118 // r12 = center to center vector
119
120 r12[0] = x[j][0] - x[i][0];
121 r12[1] = x[j][1] - x[i][1];
122 r12[2] = x[j][2] - x[i][2];
123 rsq = MathExtra::dot3(r12, r12);
124 jtype = type[j];
125
126 // compute if less than cutoff
127
128 if (rsq < cutsq[itype][jtype]) {
129 fforce[0] = fforce[1] = fforce[2] = 0.0;
130
131 switch (form[itype][jtype]) {
132
133 case SPHERE_SPHERE:
134 r2inv = 1.0 / rsq;
135 r6inv = r2inv * r2inv * r2inv;
136 forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
137 forcelj *= -r2inv;
138 if (eflag) {
139 one_eng = r6inv * (r6inv * lj3[itype][jtype] - lj4[itype][jtype]);
140 one_eng -= offset[itype][jtype];
141 }
142 fforce[0] = r12[0] * forcelj;
143 fforce[1] = r12[1] * forcelj;
144 fforce[2] = r12[2] * forcelj;
145 break;
146
147 case SPHERE_ELLIPSE:
148 precompute_i(j, wj);
149 if (newton_pair || j < nlocal) {
150 one_eng = resquared_lj(j, i, wj, r12, rsq, fforce, rtor, true);
151 tor[j][0] += rtor[0] * factor_lj;
152 tor[j][1] += rtor[1] * factor_lj;
153 tor[j][2] += rtor[2] * factor_lj;
154 } else
155 one_eng = resquared_lj(j, i, wj, r12, rsq, fforce, rtor, false);
156 break;
157
158 case ELLIPSE_SPHERE:
159 one_eng = resquared_lj(i, j, wi, r12, rsq, fforce, ttor, true);
160 tor[i][0] += ttor[0] * factor_lj;
161 tor[i][1] += ttor[1] * factor_lj;
162 tor[i][2] += ttor[2] * factor_lj;
163 break;
164
165 default:
166 precompute_i(j, wj);
167 one_eng = resquared_analytic(i, j, wi, wj, r12, rsq, fforce, ttor, rtor);
168 tor[i][0] += ttor[0] * factor_lj;
169 tor[i][1] += ttor[1] * factor_lj;
170 tor[i][2] += ttor[2] * factor_lj;
171 if (newton_pair || j < nlocal) {
172 tor[j][0] += rtor[0] * factor_lj;
173 tor[j][1] += rtor[1] * factor_lj;
174 tor[j][2] += rtor[2] * factor_lj;
175 }
176 break;
177 }
178
179 fforce[0] *= factor_lj;
180 fforce[1] *= factor_lj;
181 fforce[2] *= factor_lj;
182 f[i][0] += fforce[0];
183 f[i][1] += fforce[1];
184 f[i][2] += fforce[2];
185
186 if (newton_pair || j < nlocal) {
187 f[j][0] -= fforce[0];
188 f[j][1] -= fforce[1];
189 f[j][2] -= fforce[2];
190 }
191
192 if (eflag) evdwl = factor_lj * one_eng;
193
194 if (evflag)
195 ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fforce[0], fforce[1], fforce[2],
196 -r12[0], -r12[1], -r12[2]);
197 }
198 }
199 }
200
201 if (vflag_fdotr) virial_fdotr_compute();
202 }
203
204 /* ----------------------------------------------------------------------
205 allocate all arrays
206 ------------------------------------------------------------------------- */
207
allocate()208 void PairRESquared::allocate()
209 {
210 allocated = 1;
211 const int n = atom->ntypes + 1;
212
213 memory->create(setflag, n, n, "pair:setflag");
214 for (int i = 1; i < n; i++)
215 for (int j = i; j < n; j++) setflag[i][j] = 0;
216
217 memory->create(cutsq, n, n, "pair:cutsq");
218
219 memory->create(form, n, n, "pair:form");
220 memory->create(epsilon, n, n, "pair:epsilon");
221 memory->create(sigma, n, n, "pair:sigma");
222 memory->create(shape1, n, 3, "pair:shape1");
223 memory->create(shape2, n, 3, "pair:shape2");
224 memory->create(well, n, 3, "pair:well");
225 memory->create(cut, n, n, "pair:cut");
226 memory->create(lj1, n, n, "pair:lj1");
227 memory->create(lj2, n, n, "pair:lj2");
228 memory->create(lj3, n, n, "pair:lj3");
229 memory->create(lj4, n, n, "pair:lj4");
230 memory->create(offset, n, n, "pair:offset");
231
232 lshape = new double[n];
233 setwell = new int[n];
234 for (int i = 1; i < n; i++) setwell[i] = 0;
235 }
236
237 /* ----------------------------------------------------------------------
238 global settings
239 ------------------------------------------------------------------------- */
240
settings(int narg,char ** arg)241 void PairRESquared::settings(int narg, char **arg)
242 {
243 if (narg != 1) error->all(FLERR, "Illegal pair_style command");
244
245 cut_global = utils::numeric(FLERR, arg[0], false, lmp);
246
247 // reset cutoffs that have been explicitly set
248
249 if (allocated) {
250 int i, j;
251 for (i = 1; i <= atom->ntypes; i++)
252 for (j = i; j <= atom->ntypes; j++)
253 if (setflag[i][j]) cut[i][j] = cut_global;
254 }
255 }
256
257 /* ----------------------------------------------------------------------
258 set coeffs for one or more type pairs
259 ------------------------------------------------------------------------- */
260
coeff(int narg,char ** arg)261 void PairRESquared::coeff(int narg, char **arg)
262 {
263 if (narg < 10 || narg > 11) error->all(FLERR, "Incorrect args for pair coefficients");
264 if (!allocated) allocate();
265
266 int ilo, ihi, jlo, jhi;
267 utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
268 utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
269
270 double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
271 double sigma_one = utils::numeric(FLERR, arg[3], false, lmp);
272 double eia_one = utils::numeric(FLERR, arg[4], false, lmp);
273 double eib_one = utils::numeric(FLERR, arg[5], false, lmp);
274 double eic_one = utils::numeric(FLERR, arg[6], false, lmp);
275 double eja_one = utils::numeric(FLERR, arg[7], false, lmp);
276 double ejb_one = utils::numeric(FLERR, arg[8], false, lmp);
277 double ejc_one = utils::numeric(FLERR, arg[9], false, lmp);
278
279 double cut_one = cut_global;
280 if (narg == 11) cut_one = utils::numeric(FLERR, arg[10], false, lmp);
281
282 int count = 0;
283 for (int i = ilo; i <= ihi; i++) {
284 for (int j = MAX(jlo, i); j <= jhi; j++) {
285 epsilon[i][j] = epsilon_one;
286 sigma[i][j] = sigma_one;
287 cut[i][j] = cut_one;
288 if (eia_one != 0.0 || eib_one != 0.0 || eic_one != 0.0) {
289 well[i][0] = eia_one;
290 well[i][1] = eib_one;
291 well[i][2] = eic_one;
292 if (eia_one == 1.0 && eib_one == 1.0 && eic_one == 1.0)
293 setwell[i] = 2;
294 else
295 setwell[i] = 1;
296 }
297 if (eja_one != 0.0 || ejb_one != 0.0 || ejc_one != 0.0) {
298 well[j][0] = eja_one;
299 well[j][1] = ejb_one;
300 well[j][2] = ejc_one;
301 if (eja_one == 1.0 && ejb_one == 1.0 && ejc_one == 1.0)
302 setwell[j] = 2;
303 else
304 setwell[j] = 1;
305 }
306 setflag[i][j] = 1;
307 count++;
308 }
309 }
310
311 if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
312 }
313
314 /* ----------------------------------------------------------------------
315 init specific to this pair style
316 ------------------------------------------------------------------------- */
317
init_style()318 void PairRESquared::init_style()
319 {
320 avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
321 if (!avec) error->all(FLERR, "Pair resquared requires atom style ellipsoid");
322
323 neighbor->request(this, instance_me);
324
325 // per-type shape precalculations
326 // require that atom shapes are identical within each type
327
328 for (int i = 1; i <= atom->ntypes; i++) {
329 if (!atom->shape_consistency(i, shape1[i][0], shape1[i][1], shape1[i][2]))
330 error->all(FLERR, "Pair resquared requires atoms with same type have same shape");
331 if (setwell[i]) {
332 shape2[i][0] = shape1[i][0] * shape1[i][0];
333 shape2[i][1] = shape1[i][1] * shape1[i][1];
334 shape2[i][2] = shape1[i][2] * shape1[i][2];
335 lshape[i] = shape1[i][0] * shape1[i][1] * shape1[i][2];
336 }
337 }
338 }
339
340 /* ----------------------------------------------------------------------
341 init for one type pair i,j and corresponding j,i
342 ------------------------------------------------------------------------- */
343
init_one(int i,int j)344 double PairRESquared::init_one(int i, int j)
345 {
346 if (setwell[i] == 0 || setwell[j] == 0)
347 error->all(FLERR, "Pair resquared epsilon a,b,c coeffs are not all set");
348
349 int ishape = 0;
350 if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0) ishape = 1;
351 int jshape = 0;
352 if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0) jshape = 1;
353
354 if (ishape == 0 && jshape == 0) {
355 form[i][j] = SPHERE_SPHERE;
356 form[j][i] = SPHERE_SPHERE;
357 } else if (ishape == 0) {
358 form[i][j] = SPHERE_ELLIPSE;
359 form[j][i] = ELLIPSE_SPHERE;
360 } else if (jshape == 0) {
361 form[i][j] = ELLIPSE_SPHERE;
362 form[j][i] = SPHERE_ELLIPSE;
363 } else {
364 form[i][j] = ELLIPSE_ELLIPSE;
365 form[j][i] = ELLIPSE_ELLIPSE;
366 }
367
368 // allow mixing only for LJ spheres
369
370 if (setflag[i][j] == 0) {
371 if (setflag[j][i] == 0) {
372 if (ishape == 0 && jshape == 0) {
373 epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], sigma[i][i], sigma[j][j]);
374 sigma[i][j] = mix_distance(sigma[i][i], sigma[j][j]);
375 cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
376 } else
377 error->all(FLERR, "Pair resquared epsilon and sigma coeffs are not all set");
378 }
379 epsilon[i][j] = epsilon[j][i];
380 sigma[i][j] = sigma[j][i];
381 cut[i][j] = cut[j][i];
382 }
383
384 lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
385 lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
386 lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
387 lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
388
389 if (offset_flag && (cut[i][j] > 0.0)) {
390 double ratio = sigma[i][j] / cut[i][j];
391 offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio, 12.0) - pow(ratio, 6.0));
392 } else
393 offset[i][j] = 0.0;
394
395 epsilon[j][i] = epsilon[i][j];
396 sigma[j][i] = sigma[i][j];
397 lj1[j][i] = lj1[i][j];
398 lj2[j][i] = lj2[i][j];
399 lj3[j][i] = lj3[i][j];
400 lj4[j][i] = lj4[i][j];
401 offset[j][i] = offset[i][j];
402
403 return cut[i][j];
404 }
405
406 /* ----------------------------------------------------------------------
407 proc 0 writes to restart file
408 ------------------------------------------------------------------------- */
409
write_restart(FILE * fp)410 void PairRESquared::write_restart(FILE *fp)
411 {
412 write_restart_settings(fp);
413
414 int i, j;
415 for (i = 1; i <= atom->ntypes; i++) {
416 fwrite(&setwell[i], sizeof(int), 1, fp);
417 if (setwell[i]) fwrite(&well[i][0], sizeof(double), 3, fp);
418 for (j = i; j <= atom->ntypes; j++) {
419 fwrite(&setflag[i][j], sizeof(int), 1, fp);
420 if (setflag[i][j]) {
421 fwrite(&epsilon[i][j], sizeof(double), 1, fp);
422 fwrite(&sigma[i][j], sizeof(double), 1, fp);
423 fwrite(&cut[i][j], sizeof(double), 1, fp);
424 }
425 }
426 }
427 }
428
429 /* ----------------------------------------------------------------------
430 proc 0 reads from restart file, bcasts
431 ------------------------------------------------------------------------- */
432
read_restart(FILE * fp)433 void PairRESquared::read_restart(FILE *fp)
434 {
435 read_restart_settings(fp);
436 allocate();
437
438 int i, j;
439 int me = comm->me;
440 for (i = 1; i <= atom->ntypes; i++) {
441 if (me == 0) utils::sfread(FLERR, &setwell[i], sizeof(int), 1, fp, nullptr, error);
442 MPI_Bcast(&setwell[i], 1, MPI_INT, 0, world);
443 if (setwell[i]) {
444 if (me == 0) utils::sfread(FLERR, &well[i][0], sizeof(double), 3, fp, nullptr, error);
445 MPI_Bcast(&well[i][0], 3, MPI_DOUBLE, 0, world);
446 }
447 for (j = i; j <= atom->ntypes; j++) {
448 if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
449 MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
450 if (setflag[i][j]) {
451 if (me == 0) {
452 utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
453 utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error);
454 utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
455 }
456 MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
457 MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
458 MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
459 }
460 }
461 }
462 }
463
464 /* ----------------------------------------------------------------------
465 proc 0 writes to restart file
466 ------------------------------------------------------------------------- */
467
write_restart_settings(FILE * fp)468 void PairRESquared::write_restart_settings(FILE *fp)
469 {
470 fwrite(&cut_global, sizeof(double), 1, fp);
471 fwrite(&mix_flag, sizeof(int), 1, fp);
472 }
473
474 /* ----------------------------------------------------------------------
475 proc 0 reads from restart file, bcasts
476 ------------------------------------------------------------------------- */
477
read_restart_settings(FILE * fp)478 void PairRESquared::read_restart_settings(FILE *fp)
479 {
480 int me = comm->me;
481 if (me == 0) {
482 utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
483 utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
484 }
485 MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
486 MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
487 }
488
489 /* ----------------------------------------------------------------------
490 Precompute per-particle temporaries for RE-squared calculation
491 ------------------------------------------------------------------------- */
492
precompute_i(const int i,RE2Vars & ws)493 void PairRESquared::precompute_i(const int i, RE2Vars &ws)
494 {
495 double aTs[3][3]; // A1'*S1^2
496 int *ellipsoid = atom->ellipsoid;
497 AtomVecEllipsoid::Bonus *bonus = avec->bonus;
498 MathExtra::quat_to_mat_trans(bonus[ellipsoid[i]].quat, ws.A);
499 MathExtra::transpose_diag3(ws.A, well[atom->type[i]], ws.aTe);
500 MathExtra::transpose_diag3(ws.A, shape2[atom->type[i]], aTs);
501 MathExtra::diag_times3(shape2[atom->type[i]], ws.A, ws.sa);
502 MathExtra::times3(aTs, ws.A, ws.gamma);
503 MathExtra::rotation_generator_x(ws.A, ws.lA[0]);
504 MathExtra::rotation_generator_y(ws.A, ws.lA[1]);
505 MathExtra::rotation_generator_z(ws.A, ws.lA[2]);
506 for (int m = 0; m < 3; m++) {
507 MathExtra::times3(aTs, ws.lA[m], ws.lAtwo[m]);
508 MathExtra::transpose_times3(ws.lA[m], ws.sa, ws.lAsa[m]);
509 MathExtra::plus3(ws.lAsa[m], ws.lAtwo[m], ws.lAsa[m]);
510 }
511 }
512
513 /* ----------------------------------------------------------------------
514 Compute the derivative of the determinant of m, using m and the
515 derivative of m (m2)
516 ------------------------------------------------------------------------- */
517
518 // clang-format off
det_prime(const double m[3][3],const double m2[3][3])519 double PairRESquared::det_prime(const double m[3][3], const double m2[3][3])
520 {
521 double ans;
522 ans = m2[0][0]*m[1][1]*m[2][2] - m2[0][0]*m[1][2]*m[2][1] -
523 m[1][0]*m2[0][1]*m[2][2] + m[1][0]*m2[0][2]*m[2][1] +
524 m[2][0]*m2[0][1]*m[1][2] - m[2][0]*m2[0][2]*m[1][1] +
525 m[0][0]*m2[1][1]*m[2][2] - m[0][0]*m2[1][2]*m[2][1] -
526 m2[1][0]*m[0][1]*m[2][2] + m2[1][0]*m[0][2]*m[2][1] +
527 m[2][0]*m[0][1]*m2[1][2] - m[2][0]*m[0][2]*m2[1][1] +
528 m[0][0]*m[1][1]*m2[2][2] - m[0][0]*m[1][2]*m2[2][1] -
529 m[1][0]*m[0][1]*m2[2][2] + m[1][0]*m[0][2]*m2[2][1] +
530 m2[2][0]*m[0][1]*m[1][2] - m2[2][0]*m[0][2]*m[1][1];
531 return ans;
532 }
533 // clang-format on
534
535 /* ----------------------------------------------------------------------
536 Compute the energy, force, torque for a pair (INTEGRATED-INTEGRATED)
537 ------------------------------------------------------------------------- */
538
resquared_analytic(const int i,const int j,const RE2Vars & wi,const RE2Vars & wj,const double * r,const double rsq,double * fforce,double * ttor,double * rtor)539 double PairRESquared::resquared_analytic(const int i, const int j, const RE2Vars &wi,
540 const RE2Vars &wj, const double *r, const double rsq,
541 double *fforce, double *ttor, double *rtor)
542 {
543 int *type = atom->type;
544
545 // pair computations for energy, force, torque
546
547 double z1[3], z2[3]; // A1*rhat # don't need to store
548 double v1[3], v2[3]; // inv(S1^2)*z1 # don't need to store
549 double sigma1, sigma2; // 1/sqrt(z1'*v1)
550 double sigma1p2, sigma2p2; // sigma1^2
551 double rnorm; // L2 norm of r
552 double rhat[3]; // r/rnorm
553 double s[3]; // inv(gamma1+gamma2)*rhat
554 double sigma12; // 1/sqrt(0.5*s'*rhat)
555 double H12[3][3]; // gamma1/sigma1+gamma2/sigma2
556 double dH; // det(H12)
557 double lambda; // dS1/sigma1p2+dS2/sigma2p2
558 double nu; // sqrt(dH/(sigma1+sigma2))
559 double w[3]; // inv(A1'*E1*A1+A2'*E2*A2)*rhat
560 double h12; // rnorm-sigma12;
561 double eta; // lambda/nu
562 double chi; // 2*rhat'*w
563 double sprod; // dS1*dS2
564 double sigh; // sigma/h12
565 double tprod; // eta*chi*sigh
566 double Ua, Ur; // attractive/repulsive parts of potential
567
568 // pair computations for force, torque
569
570 double sec; // sigma*eta*chi
571 double sigma1p3, sigma2p3; // sigma1^3
572 double vsigma1[3], vsigma2[3]; // sigma1^3*v1;
573 double sigma12p3; // sigma12^3
574 double gsigma1[3][3], gsigma2[3][3]; // -gamma1/sigma1^2
575 double tsig1sig2; // eta/(2*(sigma1+sigma2))
576 double tdH; // eta/(2*dH)
577 double teta1, teta2; // 2*eta/lambda*dS1/sigma1p3
578 double fourw[3]; // 4*w;
579 double spr[3]; // 0.5*sigma12^3*s
580 double hsec; // h12+[3,b_alpha]*sec
581 double dspu; // 1/h12 - 1/hsec + temp
582 double pbsu; // 3*sigma/hsec
583 double dspr; // 7/h12-1/hsec+temp
584 double pbsr; // b_alpha*sigma/hsec;
585 double u[3]; // (-rhat(i)*rhat+eye(:,i))/rnorm
586 double u1[3], u2[3]; // A1*u
587 double dsigma1, dsigma2; // u1'*vsigma1 (force) p'*vsigma1 (tor)
588 double dH12[3][3]; // dsigma1*gsigma1 + dsigma2*gsigma2
589 double ddH; // derivative of det(H12)
590 double deta, dchi, dh12; // derivatives of eta,chi,h12
591 double dUr, dUa; // derivatives of Ua,Ur
592
593 // pair computations for torque
594
595 double fwae[3]; // -fourw'*aTe
596 double p[3]; // lA*rhat
597
598 rnorm = sqrt(rsq);
599 rhat[0] = r[0] / rnorm;
600 rhat[1] = r[1] / rnorm;
601 rhat[2] = r[2] / rnorm;
602
603 // energy
604
605 double temp[3][3];
606 MathExtra::plus3(wi.gamma, wj.gamma, temp);
607 int ierror = MathExtra::mldivide3(temp, rhat, s);
608 if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3");
609
610 sigma12 = 1.0 / sqrt(0.5 * MathExtra::dot3(s, rhat));
611 MathExtra::matvec(wi.A, rhat, z1);
612 MathExtra::matvec(wj.A, rhat, z2);
613 v1[0] = z1[0] / shape2[type[i]][0];
614 v1[1] = z1[1] / shape2[type[i]][1];
615 v1[2] = z1[2] / shape2[type[i]][2];
616 v2[0] = z2[0] / shape2[type[j]][0];
617 v2[1] = z2[1] / shape2[type[j]][1];
618 v2[2] = z2[2] / shape2[type[j]][2];
619 sigma1 = 1.0 / sqrt(MathExtra::dot3(z1, v1));
620 sigma2 = 1.0 / sqrt(MathExtra::dot3(z2, v2));
621 H12[0][0] = wi.gamma[0][0] / sigma1 + wj.gamma[0][0] / sigma2;
622 H12[0][1] = wi.gamma[0][1] / sigma1 + wj.gamma[0][1] / sigma2;
623 H12[0][2] = wi.gamma[0][2] / sigma1 + wj.gamma[0][2] / sigma2;
624 H12[1][0] = wi.gamma[1][0] / sigma1 + wj.gamma[1][0] / sigma2;
625 H12[1][1] = wi.gamma[1][1] / sigma1 + wj.gamma[1][1] / sigma2;
626 H12[1][2] = wi.gamma[1][2] / sigma1 + wj.gamma[1][2] / sigma2;
627 H12[2][0] = wi.gamma[2][0] / sigma1 + wj.gamma[2][0] / sigma2;
628 H12[2][1] = wi.gamma[2][1] / sigma1 + wj.gamma[2][1] / sigma2;
629 H12[2][2] = wi.gamma[2][2] / sigma1 + wj.gamma[2][2] / sigma2;
630 dH = MathExtra::det3(H12);
631 sigma1p2 = sigma1 * sigma1;
632 sigma2p2 = sigma2 * sigma2;
633 lambda = lshape[type[i]] / sigma1p2 + lshape[type[j]] / sigma2p2;
634 nu = sqrt(dH / (sigma1 + sigma2));
635 MathExtra::times3(wi.aTe, wi.A, temp);
636 double temp2[3][3];
637 MathExtra::times3(wj.aTe, wj.A, temp2);
638 MathExtra::plus3(temp, temp2, temp);
639 ierror = MathExtra::mldivide3(temp, rhat, w);
640 if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3");
641
642 h12 = rnorm - sigma12;
643 eta = lambda / nu;
644 chi = 2.0 * MathExtra::dot3(rhat, w);
645 sprod = lshape[type[i]] * lshape[type[j]];
646 sigh = sigma[type[i]][type[j]] / h12;
647 tprod = eta * chi * sigh;
648
649 double stemp = h12 / 2.0;
650 Ua = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
651 (shape1[type[j]][0] + stemp) * (shape1[type[j]][1] + stemp) * (shape1[type[j]][2] + stemp);
652 Ua = (1.0 + 3.0 * tprod) * sprod / Ua;
653 Ua = epsilon[type[i]][type[j]] * Ua / -36.0;
654
655 stemp = h12 / cr60;
656 Ur = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
657 (shape1[type[j]][0] + stemp) * (shape1[type[j]][1] + stemp) * (shape1[type[j]][2] + stemp);
658 Ur = (1.0 + b_alpha * tprod) * sprod / Ur;
659 Ur = epsilon[type[i]][type[j]] * Ur * pow(sigh, 6.0) / 2025.0;
660
661 // force
662
663 sec = sigma[type[i]][type[j]] * eta * chi;
664 sigma12p3 = pow(sigma12, 3.0);
665 sigma1p3 = sigma1p2 * sigma1;
666 sigma2p3 = sigma2p2 * sigma2;
667 vsigma1[0] = -sigma1p3 * v1[0];
668 vsigma1[1] = -sigma1p3 * v1[1];
669 vsigma1[2] = -sigma1p3 * v1[2];
670 vsigma2[0] = -sigma2p3 * v2[0];
671 vsigma2[1] = -sigma2p3 * v2[1];
672 vsigma2[2] = -sigma2p3 * v2[2];
673 gsigma1[0][0] = -wi.gamma[0][0] / sigma1p2;
674 gsigma1[0][1] = -wi.gamma[0][1] / sigma1p2;
675 gsigma1[0][2] = -wi.gamma[0][2] / sigma1p2;
676 gsigma1[1][0] = -wi.gamma[1][0] / sigma1p2;
677 gsigma1[1][1] = -wi.gamma[1][1] / sigma1p2;
678 gsigma1[1][2] = -wi.gamma[1][2] / sigma1p2;
679 gsigma1[2][0] = -wi.gamma[2][0] / sigma1p2;
680 gsigma1[2][1] = -wi.gamma[2][1] / sigma1p2;
681 gsigma1[2][2] = -wi.gamma[2][2] / sigma1p2;
682 gsigma2[0][0] = -wj.gamma[0][0] / sigma2p2;
683 gsigma2[0][1] = -wj.gamma[0][1] / sigma2p2;
684 gsigma2[0][2] = -wj.gamma[0][2] / sigma2p2;
685 gsigma2[1][0] = -wj.gamma[1][0] / sigma2p2;
686 gsigma2[1][1] = -wj.gamma[1][1] / sigma2p2;
687 gsigma2[1][2] = -wj.gamma[1][2] / sigma2p2;
688 gsigma2[2][0] = -wj.gamma[2][0] / sigma2p2;
689 gsigma2[2][1] = -wj.gamma[2][1] / sigma2p2;
690 gsigma2[2][2] = -wj.gamma[2][2] / sigma2p2;
691 tsig1sig2 = eta / (2.0 * (sigma1 + sigma2));
692 tdH = eta / (2.0 * dH);
693 teta1 = 2.0 * eta / lambda;
694 teta2 = teta1 * lshape[type[j]] / sigma2p3;
695 teta1 = teta1 * lshape[type[i]] / sigma1p3;
696 fourw[0] = 4.0 * w[0];
697 fourw[1] = 4.0 * w[1];
698 fourw[2] = 4.0 * w[2];
699 spr[0] = 0.5 * sigma12p3 * s[0];
700 spr[1] = 0.5 * sigma12p3 * s[1];
701 spr[2] = 0.5 * sigma12p3 * s[2];
702
703 stemp = 1.0 / (shape1[type[i]][0] * 2.0 + h12) + 1.0 / (shape1[type[i]][1] * 2.0 + h12) +
704 1.0 / (shape1[type[i]][2] * 2.0 + h12) + 1.0 / (shape1[type[j]][0] * 2.0 + h12) +
705 1.0 / (shape1[type[j]][1] * 2.0 + h12) + 1.0 / (shape1[type[j]][2] * 2.0 + h12);
706 hsec = h12 + 3.0 * sec;
707 dspu = 1.0 / h12 - 1.0 / hsec + stemp;
708 pbsu = 3.0 * sigma[type[i]][type[j]] / hsec;
709
710 stemp = 1.0 / (shape1[type[i]][0] * cr60 + h12) + 1.0 / (shape1[type[i]][1] * cr60 + h12) +
711 1.0 / (shape1[type[i]][2] * cr60 + h12) + 1.0 / (shape1[type[j]][0] * cr60 + h12) +
712 1.0 / (shape1[type[j]][1] * cr60 + h12) + 1.0 / (shape1[type[j]][2] * cr60 + h12);
713 hsec = h12 + b_alpha * sec;
714 dspr = 7.0 / h12 - 1.0 / hsec + stemp;
715 pbsr = b_alpha * sigma[type[i]][type[j]] / hsec;
716
717 for (int m = 0; m < 3; m++) {
718 u[0] = -rhat[m] * rhat[0];
719 u[1] = -rhat[m] * rhat[1];
720 u[2] = -rhat[m] * rhat[2];
721 u[m] += 1.0;
722 u[0] /= rnorm;
723 u[1] /= rnorm;
724 u[2] /= rnorm;
725 MathExtra::matvec(wi.A, u, u1);
726 MathExtra::matvec(wj.A, u, u2);
727 dsigma1 = MathExtra::dot3(u1, vsigma1);
728 dsigma2 = MathExtra::dot3(u2, vsigma2);
729 dH12[0][0] = dsigma1 * gsigma1[0][0] + dsigma2 * gsigma2[0][0];
730 dH12[0][1] = dsigma1 * gsigma1[0][1] + dsigma2 * gsigma2[0][1];
731 dH12[0][2] = dsigma1 * gsigma1[0][2] + dsigma2 * gsigma2[0][2];
732 dH12[1][0] = dsigma1 * gsigma1[1][0] + dsigma2 * gsigma2[1][0];
733 dH12[1][1] = dsigma1 * gsigma1[1][1] + dsigma2 * gsigma2[1][1];
734 dH12[1][2] = dsigma1 * gsigma1[1][2] + dsigma2 * gsigma2[1][2];
735 dH12[2][0] = dsigma1 * gsigma1[2][0] + dsigma2 * gsigma2[2][0];
736 dH12[2][1] = dsigma1 * gsigma1[2][1] + dsigma2 * gsigma2[2][1];
737 dH12[2][2] = dsigma1 * gsigma1[2][2] + dsigma2 * gsigma2[2][2];
738 ddH = det_prime(H12, dH12);
739 deta = (dsigma1 + dsigma2) * tsig1sig2;
740 deta -= ddH * tdH;
741 deta -= dsigma1 * teta1 + dsigma2 * teta2;
742 dchi = MathExtra::dot3(u, fourw);
743 dh12 = rhat[m] + MathExtra::dot3(u, spr);
744 dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu;
745 dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr;
746 fforce[m] = dUr * Ur + dUa * Ua;
747 }
748
749 // torque on i
750
751 MathExtra::vecmat(fourw, wi.aTe, fwae);
752
753 for (int i = 0; i < 3; i++) {
754 MathExtra::matvec(wi.lA[i], rhat, p);
755 dsigma1 = MathExtra::dot3(p, vsigma1);
756 dH12[0][0] = wi.lAsa[i][0][0] / sigma1 + dsigma1 * gsigma1[0][0];
757 dH12[0][1] = wi.lAsa[i][0][1] / sigma1 + dsigma1 * gsigma1[0][1];
758 dH12[0][2] = wi.lAsa[i][0][2] / sigma1 + dsigma1 * gsigma1[0][2];
759 dH12[1][0] = wi.lAsa[i][1][0] / sigma1 + dsigma1 * gsigma1[1][0];
760 dH12[1][1] = wi.lAsa[i][1][1] / sigma1 + dsigma1 * gsigma1[1][1];
761 dH12[1][2] = wi.lAsa[i][1][2] / sigma1 + dsigma1 * gsigma1[1][2];
762 dH12[2][0] = wi.lAsa[i][2][0] / sigma1 + dsigma1 * gsigma1[2][0];
763 dH12[2][1] = wi.lAsa[i][2][1] / sigma1 + dsigma1 * gsigma1[2][1];
764 dH12[2][2] = wi.lAsa[i][2][2] / sigma1 + dsigma1 * gsigma1[2][2];
765 ddH = det_prime(H12, dH12);
766 deta = tsig1sig2 * dsigma1 - tdH * ddH;
767 deta -= teta1 * dsigma1;
768 double tempv[3];
769 MathExtra::matvec(wi.lA[i], w, tempv);
770 dchi = -MathExtra::dot3(fwae, tempv);
771 MathExtra::matvec(wi.lAtwo[i], spr, tempv);
772 dh12 = -MathExtra::dot3(s, tempv);
773
774 dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu;
775 dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr;
776 ttor[i] = -(dUa * Ua + dUr * Ur);
777 }
778
779 // torque on j
780
781 if (!(force->newton_pair || j < atom->nlocal)) return Ua + Ur;
782
783 MathExtra::vecmat(fourw, wj.aTe, fwae);
784
785 for (int i = 0; i < 3; i++) {
786 MathExtra::matvec(wj.lA[i], rhat, p);
787 dsigma2 = MathExtra::dot3(p, vsigma2);
788 dH12[0][0] = wj.lAsa[i][0][0] / sigma2 + dsigma2 * gsigma2[0][0];
789 dH12[0][1] = wj.lAsa[i][0][1] / sigma2 + dsigma2 * gsigma2[0][1];
790 dH12[0][2] = wj.lAsa[i][0][2] / sigma2 + dsigma2 * gsigma2[0][2];
791 dH12[1][0] = wj.lAsa[i][1][0] / sigma2 + dsigma2 * gsigma2[1][0];
792 dH12[1][1] = wj.lAsa[i][1][1] / sigma2 + dsigma2 * gsigma2[1][1];
793 dH12[1][2] = wj.lAsa[i][1][2] / sigma2 + dsigma2 * gsigma2[1][2];
794 dH12[2][0] = wj.lAsa[i][2][0] / sigma2 + dsigma2 * gsigma2[2][0];
795 dH12[2][1] = wj.lAsa[i][2][1] / sigma2 + dsigma2 * gsigma2[2][1];
796 dH12[2][2] = wj.lAsa[i][2][2] / sigma2 + dsigma2 * gsigma2[2][2];
797 ddH = det_prime(H12, dH12);
798 deta = tsig1sig2 * dsigma2 - tdH * ddH;
799 deta -= teta2 * dsigma2;
800 double tempv[3];
801 MathExtra::matvec(wj.lA[i], w, tempv);
802 dchi = -MathExtra::dot3(fwae, tempv);
803 MathExtra::matvec(wj.lAtwo[i], spr, tempv);
804 dh12 = -MathExtra::dot3(s, tempv);
805
806 dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu;
807 dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr;
808 rtor[i] = -(dUa * Ua + dUr * Ur);
809 }
810
811 return Ua + Ur;
812 }
813
814 /* ----------------------------------------------------------------------
815 Compute the energy, force, torque for a pair (INTEGRATED-LJ)
816 ------------------------------------------------------------------------- */
817
resquared_lj(const int i,const int j,const RE2Vars & wi,const double * r,const double rsq,double * fforce,double * ttor,bool calc_torque)818 double PairRESquared::resquared_lj(const int i, const int j, const RE2Vars &wi, const double *r,
819 const double rsq, double *fforce, double *ttor, bool calc_torque)
820 {
821 int *type = atom->type;
822
823 // pair computations for energy, force, torque
824
825 double rnorm; // L2 norm of r
826 double rhat[3]; // r/rnorm
827 double s[3]; // inv(gamma1)*rhat
828 double sigma12; // 1/sqrt(0.5*s'*rhat)
829 double w[3]; // inv(A1'*E1*A1+I)*rhat
830 double h12; // rnorm-sigma12;
831 double chi; // 2*rhat'*w
832 double sigh; // sigma/h12
833 double tprod; // chi*sigh
834 double Ua, Ur; // attractive/repulsive parts of potential
835
836 // pair computations for force, torque
837
838 double sec; // sigma*chi
839 double sigma12p3; // sigma12^3
840 double fourw[3]; // 4*w;
841 double spr[3]; // 0.5*sigma12^3*s
842 double hsec; // h12+[3,b_alpha]*sec
843 double dspu; // 1/h12 - 1/hsec + temp
844 double pbsu; // 3*sigma/hsec
845 double dspr; // 7/h12-1/hsec+temp
846 double pbsr; // b_alpha*sigma/hsec;
847 double u[3]; // (-rhat(i)*rhat+eye(:,i))/rnorm
848 double dchi, dh12; // derivatives of chi,h12
849 double dUr, dUa; // derivatives of Ua,Ur
850 double h12p3; // h12^3
851
852 // pair computations for torque
853
854 double fwae[3]; // -fourw'*aTe
855 double p[3]; // lA*rhat
856
857 // distance of closest approach correction
858
859 double aTs[3][3]; // A1'*S1^2
860 double gamma[3][3]; // A1'*S1^2*A
861 double lAtwo[3][3][3]; // A1'*S1^2*wi.lA
862 double scorrect[3];
863 double half_sigma = sigma[type[i]][type[j]] / 2.0;
864 scorrect[0] = shape1[type[i]][0] + half_sigma;
865 scorrect[1] = shape1[type[i]][1] + half_sigma;
866 scorrect[2] = shape1[type[i]][2] + half_sigma;
867 scorrect[0] = scorrect[0] * scorrect[0] / 2.0;
868 scorrect[1] = scorrect[1] * scorrect[1] / 2.0;
869 scorrect[2] = scorrect[2] * scorrect[2] / 2.0;
870 MathExtra::transpose_diag3(wi.A, scorrect, aTs);
871 MathExtra::times3(aTs, wi.A, gamma);
872 for (int ii = 0; ii < 3; ii++) MathExtra::times3(aTs, wi.lA[ii], lAtwo[ii]);
873
874 rnorm = sqrt(rsq);
875 rhat[0] = r[0] / rnorm;
876 rhat[1] = r[1] / rnorm;
877 rhat[2] = r[2] / rnorm;
878
879 // energy
880
881 int ierror = MathExtra::mldivide3(gamma, rhat, s);
882 if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3");
883
884 sigma12 = 1.0 / sqrt(0.5 * MathExtra::dot3(s, rhat));
885 double temp[3][3];
886 MathExtra::times3(wi.aTe, wi.A, temp);
887 temp[0][0] += 1.0;
888 temp[1][1] += 1.0;
889 temp[2][2] += 1.0;
890 ierror = MathExtra::mldivide3(temp, rhat, w);
891 if (ierror) error->all(FLERR, "Bad matrix inversion in mldivide3");
892
893 h12 = rnorm - sigma12;
894 chi = 2.0 * MathExtra::dot3(rhat, w);
895 sigh = sigma[type[i]][type[j]] / h12;
896 tprod = chi * sigh;
897
898 h12p3 = pow(h12, 3.0);
899 double sigmap3 = pow(sigma[type[i]][type[j]], 3.0);
900 double stemp = h12 / 2.0;
901 Ua = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
902 h12p3 / 8.0;
903 Ua = (1.0 + 3.0 * tprod) * lshape[type[i]] / Ua;
904 Ua = epsilon[type[i]][type[j]] * Ua * sigmap3 * solv_f_a;
905
906 stemp = h12 / cr60;
907 Ur = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
908 h12p3 / 60.0;
909 Ur = (1.0 + b_alpha * tprod) * lshape[type[i]] / Ur;
910 Ur = epsilon[type[i]][type[j]] * Ur * sigmap3 * pow(sigh, 6.0) * solv_f_r;
911
912 // force
913
914 sec = sigma[type[i]][type[j]] * chi;
915 sigma12p3 = pow(sigma12, 3.0);
916 fourw[0] = 4.0 * w[0];
917 fourw[1] = 4.0 * w[1];
918 fourw[2] = 4.0 * w[2];
919 spr[0] = 0.5 * sigma12p3 * s[0];
920 spr[1] = 0.5 * sigma12p3 * s[1];
921 spr[2] = 0.5 * sigma12p3 * s[2];
922
923 stemp = 1.0 / (shape1[type[i]][0] * 2.0 + h12) + 1.0 / (shape1[type[i]][1] * 2.0 + h12) +
924 1.0 / (shape1[type[i]][2] * 2.0 + h12) + 3.0 / h12;
925 hsec = h12 + 3.0 * sec;
926 dspu = 1.0 / h12 - 1.0 / hsec + stemp;
927 pbsu = 3.0 * sigma[type[i]][type[j]] / hsec;
928
929 stemp = 1.0 / (shape1[type[i]][0] * cr60 + h12) + 1.0 / (shape1[type[i]][1] * cr60 + h12) +
930 1.0 / (shape1[type[i]][2] * cr60 + h12) + 3.0 / h12;
931 hsec = h12 + b_alpha * sec;
932 dspr = 7.0 / h12 - 1.0 / hsec + stemp;
933 pbsr = b_alpha * sigma[type[i]][type[j]] / hsec;
934
935 for (int m = 0; m < 3; m++) {
936 u[0] = -rhat[m] * rhat[0];
937 u[1] = -rhat[m] * rhat[1];
938 u[2] = -rhat[m] * rhat[2];
939 u[m] += 1.0;
940 u[0] /= rnorm;
941 u[1] /= rnorm;
942 u[2] /= rnorm;
943 dchi = MathExtra::dot3(u, fourw);
944 dh12 = rhat[m] + MathExtra::dot3(u, spr);
945 dUa = pbsu * dchi - dh12 * dspu;
946 dUr = pbsr * dchi - dh12 * dspr;
947 fforce[m] = dUr * Ur + dUa * Ua;
948 }
949
950 // torque on i
951
952 if (calc_torque) {
953 MathExtra::vecmat(fourw, wi.aTe, fwae);
954
955 for (int m = 0; m < 3; m++) {
956 MathExtra::matvec(wi.lA[m], rhat, p);
957 double tempv[3];
958 MathExtra::matvec(wi.lA[m], w, tempv);
959 dchi = -MathExtra::dot3(fwae, tempv);
960 MathExtra::matvec(lAtwo[m], spr, tempv);
961 dh12 = -MathExtra::dot3(s, tempv);
962
963 dUa = pbsu * dchi - dh12 * dspu;
964 dUr = pbsr * dchi - dh12 * dspr;
965 ttor[m] = -(dUa * Ua + dUr * Ur);
966 }
967 }
968
969 return Ua + Ur;
970 }
971