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_coul_long_soft_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
PairCoulLongSoftOMP(LAMMPS * lmp)39 PairCoulLongSoftOMP::PairCoulLongSoftOMP(LAMMPS *lmp) :
40 PairCoulLongSoft(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 PairCoulLongSoftOMP::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 PairCoulLongSoftOMP::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,ecoul,fpair;
92 double r,rsq,forcecoul,factor_coul;
93 double grij,expm2,prefactor,t,erfc;
94 double denc;
95 int *ilist,*jlist,*numneigh,**firstneigh;
96
97 ecoul = 0.0;
98
99 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
100 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
101 const double * _noalias const q = atom->q;
102 const int * _noalias const type = atom->type;
103 const int nlocal = atom->nlocal;
104 const double * _noalias const special_coul = force->special_coul;
105 const 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_coul = special_coul[sbmask(j)];
129 j &= NEIGHMASK;
130
131 delx = xtmp - x[j].x;
132 dely = ytmp - x[j].y;
133 delz = ztmp - x[j].z;
134 rsq = delx*delx + dely*dely + delz*delz;
135 jtype = type[j];
136
137 if (rsq < cut_coulsq) {
138
139 r = sqrt(rsq);
140 grij = g_ewald * r;
141 expm2 = exp(-grij*grij);
142 t = 1.0 / (1.0 + EWALD_P*grij);
143 erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
144
145 denc = sqrt(lam2[itype][jtype] + rsq);
146 prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / (denc*denc*denc);
147
148 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
149 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
150
151 fpair = forcecoul;
152
153 fxtmp += delx*fpair;
154 fytmp += dely*fpair;
155 fztmp += delz*fpair;
156 if (NEWTON_PAIR || j < nlocal) {
157 f[j].x -= delx*fpair;
158 f[j].y -= dely*fpair;
159 f[j].z -= delz*fpair;
160 }
161
162 if (EFLAG) {
163 prefactor = qqrd2e * lam1[itype][jtype] * qtmp*q[j] / denc;
164 ecoul = prefactor*erfc;
165 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
166 }
167
168 if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
169 0.0,ecoul,fpair,delx,dely,delz,thr);
170 }
171 }
172 f[i].x += fxtmp;
173 f[i].y += fytmp;
174 f[i].z += fztmp;
175 }
176 }
177
178 /* ---------------------------------------------------------------------- */
179
memory_usage()180 double PairCoulLongSoftOMP::memory_usage()
181 {
182 double bytes = memory_usage_thr();
183 bytes += PairCoulLongSoft::memory_usage();
184
185 return bytes;
186 }
187