1 // clang-format off
2 /* -*- c++ -*- ----------------------------------------------------------
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 Tanmoy Sanyal, M.Scott Shell, UC Santa Barbara
18 David Rosenberger, TU Darmstadt
19 ------------------------------------------------------------------------- */
20
21 #include "pair_local_density.h"
22
23 #include "atom.h"
24 #include "citeme.h"
25 #include "comm.h"
26 #include "error.h"
27 #include "force.h"
28 #include "memory.h"
29 #include "neigh_list.h"
30 #include "neighbor.h"
31
32 #include <cstring>
33
34 using namespace LAMMPS_NS;
35
36 #define MAXLINE 1024
37
38 static const char cite_pair_local_density[] =
39 "pair_style local/density command:\n\n"
40 "@Article{Sanyal16,\n"
41 " author = {T.Sanyal and M.Scott Shell},\n"
42 " title = {Coarse-grained models using local-density potentials optimized with the relative entropy: Application to implicit solvation},\n"
43 " journal = {J.~Chem.~Phys.},\n"
44 " year = 2016,\n"
45 " DOI = doi.org/10.1063/1.4958629"
46 "}\n\n"
47 "@Article{Sanyal18,\n"
48 " author = {T.Sanyal and M.Scott Shell},\n"
49 " title = {Transferable coarse-grained models of liquid-liquid equilibrium using local density potentials optimized with the relative entropy},\n"
50 " journal = {J.~Phys.~Chem. B},\n"
51 " year = 2018,\n"
52 " DOI = doi.org/10.1021/acs.jpcb.7b12446"
53 "}\n\n";
54
55 /* ---------------------------------------------------------------------- */
56
PairLocalDensity(LAMMPS * lmp)57 PairLocalDensity::PairLocalDensity(LAMMPS *lmp) : Pair(lmp)
58 {
59 restartinfo = 0;
60 one_coeff = 1;
61 single_enable = 1;
62
63 // stuff read from tabulated file
64 nLD = 0;
65 nrho = 0;
66 rho_min = nullptr;
67 rho_max = nullptr;
68 a = nullptr;
69 b = nullptr;
70 c0 = nullptr;
71 c2 = nullptr;
72 c4 = nullptr;
73 c6 = nullptr;
74 uppercut = nullptr;
75 lowercut = nullptr;
76 uppercutsq = nullptr;
77 lowercutsq = nullptr;
78 frho = nullptr;
79 rho = nullptr;
80
81 // splined arrays
82 frho_spline = nullptr;
83
84 // per-atom arrays
85 nmax = 0;
86 fp = nullptr;
87 localrho = nullptr;
88
89 // set comm size needed by this pair
90 comm_forward = 1;
91 comm_reverse = 1;
92
93 // cite publication
94 if (lmp->citeme) lmp->citeme->add(cite_pair_local_density);
95 }
96
97 /* ----------------------------------------------------------------------
98 check if allocated, since class can be destructed when incomplete
99 ------------------------------------------------------------------------- */
100
~PairLocalDensity()101 PairLocalDensity::~PairLocalDensity()
102 {
103
104 memory->destroy(localrho);
105 memory->destroy(fp);
106
107 if (allocated) {
108 memory->destroy(setflag);
109 memory->destroy(cutsq);
110 }
111
112 memory->destroy(frho_spline);
113
114 memory->destroy(rho_min);
115 memory->destroy(rho_max);
116 memory->destroy(delta_rho);
117 memory->destroy(c0);
118 memory->destroy(c2);
119 memory->destroy(c4);
120 memory->destroy(c6);
121 memory->destroy(uppercut);
122 memory->destroy(lowercut);
123 memory->destroy(uppercutsq);
124 memory->destroy(lowercutsq);
125 memory->destroy(frho);
126 memory->destroy(rho);
127
128 memory->destroy(a);
129 memory->destroy(b);
130 }
131
132 /* ---------------------------------------------------------------------- */
133
compute(int eflag,int vflag)134 void PairLocalDensity::compute(int eflag, int vflag)
135 {
136
137 int i,j,ii,jj,m,k,inum,jnum,itype,jtype;
138 double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
139 double rsqinv, phi, uLD, dphi, evdwl,fpair;
140 double p, *coeff;
141 int *ilist,*jlist,*numneigh,**firstneigh;
142
143 ev_init(eflag,vflag);
144 phi = uLD = evdwl = fpair = rsqinv = 0.0;
145
146 /* localrho = LD at each atom
147 fp = derivative of embedding energy at each atom for each LD potential
148 uLD = embedding energy of each atom due to each LD potential*/
149
150 // grow LD and fp arrays if necessary
151 // need to be atom->nmax in length
152
153 if (atom->nmax > nmax) {
154 memory->destroy(localrho);
155 memory->destroy(fp);
156 nmax = atom->nmax;
157 memory->create(localrho, nLD, nmax, "pairLD:localrho");
158 memory->create(fp, nLD, nmax, "pairLD:fp");
159 }
160
161 double **x = atom->x;
162 double **f = atom->f;
163 int *type = atom->type;
164 int nlocal = atom->nlocal;
165 int newton_pair = force->newton_pair;
166
167 inum = list->inum;
168 ilist = list->ilist;
169 numneigh = list->numneigh;
170 firstneigh = list->firstneigh;
171
172 // zero out LD and fp
173
174 if (newton_pair) {
175 m = nlocal + atom->nghost;
176 for (k = 0; k < nLD; k++) {
177 for (i = 0; i < m; i++) {
178 localrho[k][i] = 0.0;
179 fp[k][i] = 0.0;
180 }
181 }
182 }
183 else {
184 for (k = 0; k < nLD; k++) {
185 for (i = 0; i < nlocal; i++) {
186 localrho[k][i] = 0.0;
187 fp[k][i] = 0.0;
188 }
189 }
190 }
191
192 // loop over neighs of central atoms and types of LDs
193
194 for (ii = 0; ii < inum; ii++) {
195 i = ilist[ii];
196 xtmp = x[i][0];
197 ytmp = x[i][1];
198 ztmp = x[i][2];
199 itype = type[i];
200 jlist = firstneigh[i];
201 jnum = numneigh[i];
202
203 for (jj = 0; jj < jnum; jj++) {
204 j = jlist[jj];
205 j &= NEIGHMASK;
206 jtype = type[j];
207
208 // calculate distance-squared between i,j atom-types
209
210 delx = xtmp - x[j][0];
211 dely = ytmp - x[j][1];
212 delz = ztmp - x[j][2];
213 rsq = delx*delx + dely*dely + delz*delz;
214
215 // calculating LDs based on central and neigh filters
216
217 for (k = 0; k < nLD; k++) {
218 if (rsq < lowercutsq[k]) {
219 phi = 1.0;
220 }
221 else if (rsq > uppercutsq[k]) {
222 phi = 0.0;
223 }
224 else {
225 phi = c0[k] + rsq * (c2[k] + rsq * (c4[k] + c6[k]*rsq));
226 }
227 localrho[k][i] += (phi * b[k][jtype]);
228
229 /*checking for both i,j is necessary
230 since a half neighbor list is processed.*/
231
232 if (newton_pair || j<nlocal) {
233 localrho[k][j] += (phi * b[k][itype]);
234 }
235 }
236 }
237 }
238
239 // communicate and sum LDs over all procs
240 if (newton_pair) comm->reverse_comm_pair(this);
241
242 //
243
244 for (ii = 0; ii < inum; ii++) {
245 i = ilist[ii];
246 itype = type[i];
247 uLD = 0.0;
248
249 for (k = 0; k < nLD; k++) {
250
251 /*skip over this loop if the LD potential
252 is not intendend for central atomtype <itype>*/
253 if (!(a[k][itype])) continue;
254
255 // linear extrapolation at rho_min and rho_max
256
257 if (localrho[k][i] <= rho_min[k]) {
258 coeff = frho_spline[k][0];
259 fp[k][i] = coeff[2];
260 uLD += a[k][itype] * ( coeff[6] + fp[k][i]*(localrho[k][i] - rho_min[k]) );
261 }
262 else if (localrho[k][i] >= rho_max[k]) {
263 coeff = frho_spline[k][nrho-2];
264 fp[k][i] = coeff[0] + coeff[1] + coeff[2];
265 uLD += a[k][itype] * ( (coeff[3] + coeff[4] + coeff[5] + coeff[6]) + fp[k][i]*(localrho[k][i] - rho_max[k]) );
266 }
267 else {
268 p = (localrho[k][i] - rho_min[k]) / delta_rho[k];
269 m = static_cast<int> (p);
270 m = MAX(0, MIN(m, nrho-2));
271 p -= m;
272 p = MIN(p, 1.0);
273 coeff = frho_spline[k][m];
274 fp[k][i] = (coeff[0]*p + coeff[1])*p + coeff[2];
275 uLD += a[k][itype] * (((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
276 }
277 }
278
279 if (eflag) {
280 if (eflag_global) eng_vdwl += uLD;
281 if (eflag_atom) eatom[i] += uLD;
282 }
283 }
284
285 // communicate LD and fp to all procs
286
287 comm->forward_comm_pair(this);
288
289 // compute forces on each atom
290 // loop over neighbors of my atoms
291
292 for (ii = 0; ii < inum; ii++) {
293 i = ilist[ii];
294 xtmp = x[i][0];
295 ytmp = x[i][1];
296 ztmp = x[i][2];
297 itype = type[i];
298
299 jlist = firstneigh[i];
300 jnum = numneigh[i];
301
302 for (jj = 0; jj < jnum; jj++) {
303 j = jlist[jj];
304 j &= NEIGHMASK;
305 jtype = type[j];
306
307 // calculate square of distance between i,j atoms
308
309 delx = xtmp - x[j][0];
310 dely = ytmp - x[j][1];
311 delz = ztmp - x[j][2];
312 rsq = delx*delx + dely*dely + delz*delz;
313
314 // calculate force between two atoms
315 fpair = 0.0;
316 if (rsq < cutforcesq) { // global cutoff check
317 rsqinv = 1.0/rsq;
318 for (k = 0; k < nLD; k++) {
319 if (rsq >= lowercutsq[k] && rsq < uppercutsq[k]) {
320 dphi = rsq * (2.0*c2[k] + rsq * (4.0*c4[k] + 6.0*c6[k]*rsq));
321 fpair += -(a[k][itype]*b[k][jtype]*fp[k][i] + a[k][jtype]*b[k][itype]*fp[k][j]) * dphi;
322 }
323 }
324 fpair *= rsqinv;
325
326 f[i][0] += delx*fpair;
327 f[i][1] += dely*fpair;
328 f[i][2] += delz*fpair;
329 if (newton_pair || j < nlocal) {
330 f[j][0] -= delx*fpair;
331 f[j][1] -= dely*fpair;
332 f[j][2] -= delz*fpair;
333 }
334
335 /*eng_vdwl has already been completely built,
336 so no need to add anything here*/
337
338 if (eflag) evdwl = 0.0;
339
340 if (evflag) ev_tally(i,j,nlocal,newton_pair,
341 evdwl,0.0,fpair,delx,dely,delz);
342 }
343
344 }
345 }
346
347 if (vflag_fdotr) virial_fdotr_compute();
348 }
349
350
351 /* ----------------------------------------------------------------------
352 allocate all arrays
353 ------------------------------------------------------------------------- */
354
allocate()355 void PairLocalDensity::allocate()
356 {
357 allocated = 1;
358 int n = atom->ntypes;
359
360 memory->create(cutsq,n+1,n+1,"pair:cutsq");
361
362 memory->create(setflag,n+1,n+1,"pair:setflag");
363 for (int i = 1; i <= n; i++)
364 for (int j = i; j <= n; j++)
365 setflag[i][j] = 0;
366 }
367
368 /* ----------------------------------------------------------------------
369 global settings
370 ------------------------------------------------------------------------- */
371
settings(int narg,char **)372 void PairLocalDensity::settings(int narg, char ** /* arg */)
373 {
374 if (narg > 0) error->all(FLERR,"Illegal pair_style command");
375 }
376
377 /* ----------------------------------------------------------------------
378 set coeffs for all type pairs
379 read tabulated LD input file
380 ------------------------------------------------------------------------- */
381
coeff(int narg,char ** arg)382 void PairLocalDensity::coeff(int narg, char **arg)
383 {
384 int i, j;
385 if (!allocated) allocate();
386
387 if (narg != 3) error->all(FLERR,"Incorrect args for pair coefficients");
388
389 // insure I,J args are * *
390
391 if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
392 error->all(FLERR,"Incorrect args for pair coefficients");
393
394 // parse LD file
395
396 parse_file(arg[2]);
397
398
399 // clear setflag since coeff() called once with I,J = * *
400
401 for (i = 1; i <= atom->ntypes; i++)
402 for (j = i; j <= atom->ntypes; j++)
403 setflag[i][j] = 0;
404
405 // set setflag for all i,j type pairs
406
407 int count = 0;
408 for (i = 1; i <= atom->ntypes; i++) {
409 for (j = i; j <= atom->ntypes; j++) {
410 setflag[i][j] = 1;
411 count++;
412 }
413 }
414 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
415 }
416
417 /* ----------------------------------------------------------------------
418 init specific to this pair style
419 ------------------------------------------------------------------------- */
420
init_style()421 void PairLocalDensity::init_style()
422 {
423 // spline rho and frho arrays
424 // request half neighbor list
425
426 array2spline();
427
428 // half neighbor request
429 neighbor->request(this);
430 }
431
432 /* ----------------------------------------------------------------------
433 init for one type pair i,j and corresponding j,i
434 ------------------------------------------------------------------------- */
435
init_one(int,int)436 double PairLocalDensity::init_one(int /* i */, int /* j */)
437 {
438 // single global cutoff = max of all uppercuts read in from LD file
439
440 cutmax = 0.0;
441 for (int k = 0; k < nLD; k++)
442 cutmax = MAX(cutmax,uppercut[k]);
443
444 cutforcesq = cutmax*cutmax;
445
446 return cutmax;
447 }
448
449
450 /*--------------------------------------------------------------------------
451 pair_write functionality for this pair style that gives just a snap-shot
452 of the LD potential without doing an actual MD run
453 ---------------------------------------------------------------------------*/
454
single(int,int,int itype,int jtype,double rsq,double,double,double & fforce)455 double PairLocalDensity::single(int /* i */, int /* j */, int itype, int jtype,
456 double rsq, double /* factor_coul */,
457 double /* factor_lj */, double &fforce)
458 {
459 int m, k, index;
460 double rsqinv, p, uLD;
461 double *coeff, **LD;
462 double dFdrho, phi, dphi;
463
464 uLD = dFdrho = dphi = 0.0;
465
466 memory->create(LD, nLD, 3, "pairLD:LD");
467 for (k = 0; k < nLD; k++) {
468 LD[k][1] = 0.0; // itype:- 1
469 LD[k][2] = 0.0; // jtype:- 2
470 }
471
472 rsqinv = 1.0/rsq;
473 for (k = 0; k < nLD; k++) {
474 if (rsq < lowercutsq[k]) {
475 phi = 1.0;
476 }
477 else if (rsq > uppercutsq[k]) {
478 phi = 0.0;
479 }
480 else {
481 phi = c0[k] + rsq * (c2[k] + rsq * (c4[k] + c6[k]*rsq));
482 }
483 LD[k][1] += (phi * b[k][jtype]);
484 LD[k][2] += (phi * b[k][itype]);
485 }
486
487 for (k = 0; k < nLD; k++) {
488 if (a[k][itype]) index = 1;
489 if (a[k][jtype]) index = 2;
490
491 if (LD[k][index] <= rho_min[k]) {
492 coeff = frho_spline[k][0];
493 dFdrho = coeff[2];
494 uLD += a[k][itype] * ( coeff[6] + dFdrho*(LD[k][index] - rho_min[k]) );
495 }
496 else if (LD[k][index] >= rho_max[k]) {
497 coeff = frho_spline[k][nrho-1];
498 dFdrho = coeff[0] + coeff[1] + coeff[2];
499 uLD += a[k][itype] * ( (coeff[3] + coeff[4] + coeff[5] + coeff[6]) + dFdrho*(LD[k][index] - rho_max[k]) );
500 }
501 else {
502 p = (LD[k][index] - rho_min[k]) / delta_rho[k];
503 m = static_cast<int> (p);
504 m = MAX(0, MIN(m, nrho-2));
505 p -= m;
506 p = MIN(p, 1.0);
507 coeff = frho_spline[k][m];
508 dFdrho = (coeff[0]*p + coeff[1])*p + coeff[2];
509 uLD += a[k][itype] * (((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
510 }
511
512 if (rsq < lowercutsq[k]) {
513 dphi = 0.0;
514 }
515 else if (rsq > uppercutsq[k]) {
516 dphi = 0.0;
517 }
518 else {
519 dphi = rsq * (2.0*c2[k] + rsq * (4.0*c4[k] + 6.0*c6[k]*rsq));
520 }
521 fforce += -(a[k][itype]*b[k][jtype]*dFdrho + a[k][jtype]*b[k][itype]*dFdrho) * dphi *rsqinv;
522 }
523 memory->destroy(LD);
524
525 return uLD;
526 }
527
528 /*--------------------------------------------------------------------
529 spline the array frho read in from the file to create
530 frho_spline
531 ---------------------------------------------------------------------- */
532
array2spline()533 void PairLocalDensity::array2spline() {
534 memory->destroy(frho_spline);
535 memory->create(frho_spline, nLD, nrho, 7, "pairLD:frho_spline");
536
537 for (int k = 0; k < nLD; k++)
538 interpolate_cbspl(nrho, delta_rho[k], frho[k], frho_spline[k]);
539
540 }
541
542 /* ----------------------------------------------------------------------
543 (one-dimensional) cubic spline interpolation sub-routine,
544 which determines the coeffs for a clamped cubic spline
545 given tabulated data
546 ------------------------------------------------------------------------*/
547
interpolate_cbspl(int n,double delta,double * f,double ** spline)548 void PairLocalDensity::interpolate_cbspl(int n, double delta,
549 double *f, double **spline)
550 {
551 /* inputs:
552 n number of interpolating points
553
554 f array containing function values to
555 be interpolated; f[i] is the function
556 value corresponding to x[i]
557 ('x' refers to the independent var)
558
559 delta difference in tabulated values of x
560
561 outputs: (packaged as columns of the coeff matrix)
562 coeff_b coeffs of linear terms
563 coeff_c coeffs of quadratic terms
564 coeff_d coeffs of cubic terms
565 spline matrix that collects b,c,d
566
567
568 other parameters:
569 fpa derivative of function at x=a
570 fpb derivative of function at x=b
571 */
572
573 double *dl, *dd, *du;
574 double *coeff_b, *coeff_c, *coeff_d;
575 double fpa, fpb;
576
577 int i;
578
579 coeff_b = new double [n];
580 coeff_c = new double [n];
581 coeff_d = new double [n];
582 dl = new double [n];
583 dd = new double [n];
584 du = new double [n];
585
586 // initialize values
587 for ( i = 0; i<n; i++) {
588 coeff_b[i] = coeff_c[i] = coeff_d[i] = 0.;
589 dl[i] = dd[i] = du[i] = 0.;
590 }
591
592 // set slopes at beginning and end
593 fpa = 0.;
594 fpb = 0.;
595
596 for (i = 0; i < n-1; i++) {
597 dl[i] = du[i] = delta;
598 }
599
600 dd[0] = 2.0 * delta;
601 dd[n-1] = 2.0 * delta;
602 coeff_c[0] = ( 3.0 / delta ) * ( f[1] - f[0] ) - 3.0 * fpa;
603 coeff_c[n-1] = 3.0 * fpb - ( 3.0 / delta ) * ( f[n-1] - f[n-2] );
604 for (i = 0; i < n-2; i++) {
605 dd[i+1] = 4.0 * delta;
606 coeff_c[i+1] = ( 3.0 / delta ) * ( f[i+2] - f[i+1] ) -
607 ( 3.0 / delta ) * ( f[i+1] - f[i] );
608 }
609
610 // tridiagonal solver
611 for (i = 0; i < n-1; i++) {
612 du[i] /= dd[i];
613 dd[i+1] -= dl[i]*du[i];
614 }
615
616 coeff_c[0] /= dd[0];
617 for (i = 1; i < n; i++)
618 coeff_c[i] = ( coeff_c[i] - dl[i-1] * coeff_c[i-1] ) / dd[i];
619
620 for (i = n-2; i >= 0; i--)
621 coeff_c[i] -= coeff_c[i+1] * du[i];
622
623 for (i = 0; i < n-1; i++) {
624 coeff_d[i] = ( coeff_c[i+1] - coeff_c[i] ) / ( 3.0 * delta );
625 coeff_b[i] = ( f[i+1] - f[i] ) / delta - delta * ( coeff_c[i+1] + 2.0*coeff_c[i] ) / 3.0;
626 }
627
628 // normalize
629 for (i = 0; i < n-1; i++) {
630 coeff_b[i] = coeff_b[i] * delta ;
631 coeff_c[i] = coeff_c[i] * delta*delta ;
632 coeff_d[i] = coeff_d[i] * delta*delta*delta;
633 }
634
635 //copy to coefficient matrix
636 for ( i = 0; i < n; i++) {
637 spline[i][3] = coeff_d[i];
638 spline[i][4] = coeff_c[i];
639 spline[i][5] = coeff_b[i];
640 spline[i][6] = f[i];
641 spline[i][2] = spline[i][5]/delta;
642 spline[i][1] = 2.0*spline[i][4]/delta;
643 spline[i][0] = 3.0*spline[i][3]/delta;
644 }
645
646 delete [] coeff_b;
647 delete [] coeff_c;
648 delete [] coeff_d;
649 delete [] du;
650 delete [] dd;
651 delete [] dl;
652 }
653
654 /* ----------------------------------------------------------------------
655 read potential values from tabulated LD input file
656 ------------------------------------------------------------------------- */
657
parse_file(char * filename)658 void PairLocalDensity::parse_file(char *filename) {
659
660 int k, n;
661 int me = comm->me;
662 FILE *fptr;
663 char line[MAXLINE];
664 double ratio, lc2, uc2, denom;
665
666 if (me == 0) {
667 fptr = fopen(filename, "r");
668 if (fptr == nullptr)
669 error->one(FLERR,"Cannot open Local Density potential file {}: {}",filename,utils::getsyserror());
670 }
671
672 double *ftmp; // tmp var to extract the complete 2D frho array from file
673
674 // broadcast number of LD potentials and number of (rho,frho) pairs
675 if (me == 0) {
676
677 // first 2 comment lines ignored
678 utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
679 utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
680
681 // extract number of potentials and number of (frho, rho) points
682 utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
683 sscanf(line, "%d %d", &nLD, &nrho);
684 utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
685 }
686
687 MPI_Bcast(&nLD,1,MPI_INT,0,world);
688 MPI_Bcast(&nrho,1,MPI_INT,0,world);
689
690 // setting up all arrays to be read from files and broadcasted
691 memory->create(uppercut, nLD, "pairLD:uppercut");
692 memory->create(lowercut, nLD, "pairLD:lowercut");
693 memory->create(uppercutsq, nLD, "pairLD:uppercutsq");
694 memory->create(lowercutsq, nLD, "pairLD:lowercutsq");
695 memory->create(c0, nLD, "pairLD:c0");
696 memory->create(c2, nLD, "pairLD:c2");
697 memory->create(c4, nLD, "pairLD:c4");
698 memory->create(c6, nLD, "pairLD:c6");
699 memory->create(rho_min, nLD, "pairLD:rho_min");
700 memory->create(rho_max, nLD, "pairLD:rho_max");
701 memory->create(delta_rho, nLD,"pairLD:delta_rho");
702 memory->create(ftmp, nrho*nLD, "pairLD:ftmp");
703
704 // setting up central and neighbor atom filters
705 memory->create(a, nLD, atom->ntypes+1 , "pairLD:a");
706 memory->create(b, nLD, atom->ntypes+1, "pairLD:b");
707 if (me == 0) {
708 for (n = 1; n <= atom->ntypes; n++) {
709 for (k = 0; k < nLD; k++) {
710 a[k][n] = 0;
711 b[k][n] = 0;
712 }
713 }
714 }
715
716 // read file block by block
717
718 if (me == 0) {
719 for (k = 0; k < nLD; k++) {
720
721 // parse upper and lower cut values
722 if (fgets(line,MAXLINE,fptr)==nullptr) break;
723 sscanf(line, "%lf %lf", &lowercut[k], &uppercut[k]);
724
725 // parse and broadcast central atom filter
726 utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error);
727 char *tmp = strtok(line, " /t/n/r/f");
728 while (tmp != nullptr) {
729 a[k][atoi(tmp)] = 1;
730 tmp = strtok(nullptr, " /t/n/r/f");
731 }
732
733 // parse neighbor atom filter
734 utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error);
735 tmp = strtok(line, " /t/n/r/f");
736 while (tmp != nullptr) {
737 b[k][atoi(tmp)] = 1;
738 tmp = strtok(nullptr, " /t/n/r/f");
739 }
740
741 // parse min, max and delta rho values
742 utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error);
743 sscanf(line, "%lf %lf %lf", &rho_min[k], &rho_max[k], &delta_rho[k]);
744 // recompute delta_rho from scratch for precision
745 delta_rho[k] = (rho_max[k] - rho_min[k]) / (nrho - 1);
746
747 // parse tabulated frho values from each line into temporary array
748 for (n = 0; n < nrho; n++) {
749 utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
750 sscanf(line, "%lf", &ftmp[k*nrho + n]);
751 }
752
753 // ignore blank line at the end of every block
754 utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
755
756 // set coefficients for local density indicator function
757 uc2 = uppercut[k] * uppercut[k];
758 uppercutsq[k] = uc2;
759 lc2 = lowercut[k] * lowercut[k];
760 lowercutsq[k] = lc2;
761 ratio = lc2/uc2;
762 denom = 1.0 - ratio;
763 denom = denom*denom*denom;
764 c0[k] = (1 - 3.0 * ratio) / denom;
765 c2[k] = (6.0 * ratio) / (uc2 * denom);
766 c4[k] = -(3.0 + 3.0*ratio) / (uc2*uc2 * denom);
767 c6[k] = 2.0 / (uc2*uc2*uc2 * denom);
768 }
769 }
770
771 // Broadcast all parsed arrays
772 MPI_Bcast(&lowercut[0], nLD, MPI_DOUBLE, 0, world);
773 MPI_Bcast(&uppercut[0], nLD, MPI_DOUBLE, 0, world);
774 MPI_Bcast(&lowercutsq[0], nLD, MPI_DOUBLE, 0, world);
775 MPI_Bcast(&uppercutsq[0], nLD, MPI_DOUBLE, 0, world);
776 MPI_Bcast(&c0[0], nLD, MPI_DOUBLE, 0, world);
777 MPI_Bcast(&c2[0], nLD, MPI_DOUBLE, 0, world);
778 MPI_Bcast(&c4[0], nLD, MPI_DOUBLE, 0, world);
779 MPI_Bcast(&c6[0], nLD, MPI_DOUBLE, 0, world);
780 for (k = 0; k < nLD; k++) {
781 MPI_Bcast(&a[k][1], atom->ntypes, MPI_INT, 0, world);
782 MPI_Bcast(&b[k][1], atom->ntypes, MPI_INT, 0, world);
783 }
784 MPI_Bcast(&rho_min[0], nLD, MPI_DOUBLE, 0, world);
785 MPI_Bcast(&rho_max[0], nLD, MPI_DOUBLE, 0, world);
786 MPI_Bcast(&delta_rho[0], nLD, MPI_DOUBLE, 0, world);
787 MPI_Bcast(&ftmp[0], nLD*nrho, MPI_DOUBLE, 0, world);
788
789 if (me == 0) fclose(fptr);
790
791 // set up rho and frho arrays
792 memory->create(rho, nLD, nrho, "pairLD:rho");
793 memory->create(frho, nLD, nrho, "pairLD:frho");
794
795 for (k = 0; k < nLD; k++) {
796 for (n = 0; n < nrho; n++) {
797 rho[k][n] = rho_min[k] + n*delta_rho[k];
798 frho[k][n] = ftmp[k*nrho + n];
799 }
800 }
801
802 // delete temporary array
803 memory->destroy(ftmp);
804 }
805
806 /* ----------------------------------------------------------------------
807 communication routines
808 ------------------------------------------------------------------------- */
809
pack_comm(int n,int * list,double * buf,int,int *)810 int PairLocalDensity::pack_comm(int n, int *list, double *buf,
811 int /* pbc_flag */, int * /* pbc */) {
812 int i,j,k;
813 int m;
814
815 m = 0;
816 for (i = 0; i < n; i++) {
817 j = list[i];
818 for (k = 0; k < nLD; k++) {
819 buf[m++] = fp[k][j];
820 }
821 }
822
823 return nLD;
824 }
825
826 /* ---------------------------------------------------------------------- */
827
unpack_comm(int n,int first,double * buf)828 void PairLocalDensity::unpack_comm(int n, int first, double *buf) {
829
830 int i,k,m,last;
831
832 m = 0;
833 last = first + n;
834 for (i = first; i < last; i++) {
835 for (k = 0; k < nLD; k++) {
836 fp[k][i] = buf[m++];
837 }
838 }
839 }
840
841 /* ---------------------------------------------------------------------- */
842
pack_reverse_comm(int n,int first,double * buf)843 int PairLocalDensity::pack_reverse_comm(int n, int first, double *buf) {
844
845 int i,k,m,last;
846
847 m = 0;
848 last = first + n;
849 for (i = first; i < last; i++) {
850 for (k = 0; k < nLD; k++) {
851 buf[m++] = localrho[k][i];
852 }
853 }
854 return nLD;
855 }
856
857 /* ---------------------------------------------------------------------- */
858
unpack_reverse_comm(int n,int * list,double * buf)859 void PairLocalDensity::unpack_reverse_comm(int n, int *list, double *buf) {
860
861 int i,j,k;
862 int m;
863
864 m = 0;
865 for (i = 0; i < n; i++) {
866 j = list[i];
867 for (k = 0; k < nLD; k++) {
868 localrho[k][j] += buf[m++];
869 }
870 }
871 }
872
873 /* ----------------------------------------------------------------------
874 memory usage of local atom-based arrays
875 ------------------------------------------------------------------------- */
876
memory_usage()877 double PairLocalDensity::memory_usage()
878 {
879 double bytes = (double)maxeatom * sizeof(double);
880 bytes += (double)maxvatom*6 * sizeof(double);
881 bytes += (double)2 * (nmax*nLD) * sizeof(double);
882 return bytes;
883 }
884
885