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 This software is distributed under the GNU General Public License.
8
9 See the README file in the top-level LAMMPS directory.
10 ------------------------------------------------------------------------- */
11
12 /* ----------------------------------------------------------------------
13 Contributing author: Axel Kohlmeyer (Temple U)
14 ------------------------------------------------------------------------- */
15
16 #include "pair_born_coul_long_omp.h"
17
18 #include "atom.h"
19 #include "comm.h"
20 #include "force.h"
21 #include "neigh_list.h"
22 #include "suffix.h"
23
24 #include <cmath>
25
26 #include "omp_compat.h"
27 using namespace LAMMPS_NS;
28
29 #define EWALD_F 1.12837917
30 #define EWALD_P 0.3275911
31 #define A1 0.254829592
32 #define A2 -0.284496736
33 #define A3 1.421413741
34 #define A4 -1.453152027
35 #define A5 1.061405429
36
37 /* ---------------------------------------------------------------------- */
38
PairBornCoulLongOMP(LAMMPS * lmp)39 PairBornCoulLongOMP::PairBornCoulLongOMP(LAMMPS *lmp) :
40 PairBornCoulLong(lmp), ThrOMP(lmp, THR_PAIR)
41 {
42 suffix_flag |= Suffix::OMP;
43 respa_enable = 0;
44 }
45
46 /* ---------------------------------------------------------------------- */
47
compute(int eflag,int vflag)48 void PairBornCoulLongOMP::compute(int eflag, int vflag)
49 {
50 ev_init(eflag,vflag);
51
52 const int nall = atom->nlocal + atom->nghost;
53 const int nthreads = comm->nthreads;
54 const int inum = list->inum;
55
56 #if defined(_OPENMP)
57 #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
58 #endif
59 {
60 int ifrom, ito, tid;
61
62 loop_setup_thr(ifrom, ito, tid, inum, nthreads);
63 ThrData *thr = fix->get_thr(tid);
64 thr->timer(Timer::START);
65 ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
66
67 if (evflag) {
68 if (eflag) {
69 if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
70 else eval<1,1,0>(ifrom, ito, thr);
71 } else {
72 if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
73 else eval<1,0,0>(ifrom, ito, thr);
74 }
75 } else {
76 if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
77 else eval<0,0,0>(ifrom, ito, thr);
78 }
79
80 thr->timer(Timer::PAIR);
81 reduce_thr(this, eflag, vflag, thr);
82 } // end of omp parallel region
83 }
84
85 /* ---------------------------------------------------------------------- */
86
87 template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
eval(int iifrom,int iito,ThrData * const thr)88 void PairBornCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
89 {
90 int i,j,ii,jj,jnum,itype,jtype;
91 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
92 double rsq,r2inv,r6inv,r,rexp,forcecoul,forceborn,factor_coul,factor_lj;
93 double grij,expm2,prefactor,t,erfc;
94 int *ilist,*jlist,*numneigh,**firstneigh;
95
96 evdwl = ecoul = 0.0;
97
98 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
99 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
100 const double * _noalias const q = atom->q;
101 const int * _noalias const type = atom->type;
102 int nlocal = atom->nlocal;
103 double *special_coul = force->special_coul;
104 double *special_lj = force->special_lj;
105 double qqrd2e = force->qqrd2e;
106 double fxtmp,fytmp,fztmp;
107
108 ilist = list->ilist;
109 numneigh = list->numneigh;
110 firstneigh = list->firstneigh;
111
112 // loop over neighbors of my atoms
113
114 for (ii = iifrom; ii < iito; ++ii) {
115
116 i = ilist[ii];
117 qtmp = q[i];
118 xtmp = x[i].x;
119 ytmp = x[i].y;
120 ztmp = x[i].z;
121 itype = type[i];
122 jlist = firstneigh[i];
123 jnum = numneigh[i];
124 fxtmp=fytmp=fztmp=0.0;
125
126 for (jj = 0; jj < jnum; jj++) {
127 j = jlist[jj];
128 factor_lj = special_lj[sbmask(j)];
129 factor_coul = special_coul[sbmask(j)];
130 j &= NEIGHMASK;
131
132 delx = xtmp - x[j].x;
133 dely = ytmp - x[j].y;
134 delz = ztmp - x[j].z;
135 rsq = delx*delx + dely*dely + delz*delz;
136 jtype = type[j];
137
138 if (rsq < cutsq[itype][jtype]) {
139 r2inv = 1.0/rsq;
140 r = sqrt(rsq);
141
142 if (rsq < cut_coulsq) {
143 grij = g_ewald * r;
144 expm2 = exp(-grij*grij);
145 t = 1.0 / (1.0 + EWALD_P*grij);
146 erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
147 prefactor = qqrd2e * qtmp*q[j]/r;
148 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
149 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
150 } else forcecoul = 0.0;
151
152 if (rsq < cut_ljsq[itype][jtype]) {
153 r6inv = r2inv*r2inv*r2inv;
154 rexp = exp((sigma[itype][jtype]-r)*rhoinv[itype][jtype]);
155 forceborn = born1[itype][jtype]*r*rexp - born2[itype][jtype]*r6inv
156 + born3[itype][jtype]*r2inv*r6inv;
157 } else forceborn = 0.0;
158
159 fpair = (forcecoul + factor_lj*forceborn)*r2inv;
160
161 fxtmp += delx*fpair;
162 fytmp += dely*fpair;
163 fztmp += delz*fpair;
164 if (NEWTON_PAIR || j < nlocal) {
165 f[j].x -= delx*fpair;
166 f[j].y -= dely*fpair;
167 f[j].z -= delz*fpair;
168 }
169
170 if (EFLAG) {
171 if (rsq < cut_coulsq) {
172 ecoul = prefactor*erfc;
173 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
174 } else ecoul = 0.0;
175 if (rsq < cut_ljsq[itype][jtype]) {
176 evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv
177 + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype];
178 evdwl *= factor_lj;
179 } else evdwl = 0.0;
180 }
181
182 if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
183 evdwl,ecoul,fpair,delx,dely,delz,thr);
184 }
185 }
186 f[i].x += fxtmp;
187 f[i].y += fytmp;
188 f[i].z += fztmp;
189 }
190 }
191
192 /* ---------------------------------------------------------------------- */
193
memory_usage()194 double PairBornCoulLongOMP::memory_usage()
195 {
196 double bytes = memory_usage_thr();
197 bytes += PairBornCoulLong::memory_usage();
198
199 return bytes;
200 }
201