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 #include "pair_lj_cut_coul_debye.h"
16 
17 #include <cmath>
18 #include "atom.h"
19 #include "neigh_list.h"
20 #include "force.h"
21 #include "comm.h"
22 #include "error.h"
23 
24 
25 using namespace LAMMPS_NS;
26 
27 /* ---------------------------------------------------------------------- */
28 
PairLJCutCoulDebye(LAMMPS * lmp)29 PairLJCutCoulDebye::PairLJCutCoulDebye(LAMMPS *lmp) : PairLJCutCoulCut(lmp) {}
30 
31 /* ---------------------------------------------------------------------- */
32 
compute(int eflag,int vflag)33 void PairLJCutCoulDebye::compute(int eflag, int vflag)
34 {
35   int i,j,ii,jj,inum,jnum,itype,jtype;
36   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
37   double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
38   double r,rinv,screening;
39   int *ilist,*jlist,*numneigh,**firstneigh;
40 
41   evdwl = ecoul = 0.0;
42   ev_init(eflag,vflag);
43 
44   double **x = atom->x;
45   double **f = atom->f;
46   double *q = atom->q;
47   int *type = atom->type;
48   int nlocal = atom->nlocal;
49   double *special_coul = force->special_coul;
50   double *special_lj = force->special_lj;
51   int newton_pair = force->newton_pair;
52   double qqrd2e = force->qqrd2e;
53 
54   inum = list->inum;
55   ilist = list->ilist;
56   numneigh = list->numneigh;
57   firstneigh = list->firstneigh;
58 
59   // loop over neighbors of my atoms
60 
61   for (ii = 0; ii < inum; ii++) {
62     i = ilist[ii];
63     qtmp = q[i];
64     xtmp = x[i][0];
65     ytmp = x[i][1];
66     ztmp = x[i][2];
67     itype = type[i];
68     jlist = firstneigh[i];
69     jnum = numneigh[i];
70 
71     for (jj = 0; jj < jnum; jj++) {
72       j = jlist[jj];
73       factor_lj = special_lj[sbmask(j)];
74       factor_coul = special_coul[sbmask(j)];
75       j &= NEIGHMASK;
76 
77       delx = xtmp - x[j][0];
78       dely = ytmp - x[j][1];
79       delz = ztmp - x[j][2];
80       rsq = delx*delx + dely*dely + delz*delz;
81       jtype = type[j];
82 
83       if (rsq < cutsq[itype][jtype]) {
84         r2inv = 1.0/rsq;
85 
86         if (rsq < cut_coulsq[itype][jtype]) {
87           r = sqrt(rsq);
88           rinv = 1.0/r;
89           screening = exp(-kappa*r);
90           forcecoul = qqrd2e * qtmp*q[j] * screening * (kappa + rinv);
91         } else forcecoul = 0.0;
92 
93         if (rsq < cut_ljsq[itype][jtype]) {
94           r6inv = r2inv*r2inv*r2inv;
95           forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
96         } else forcelj = 0.0;
97 
98         fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
99 
100         f[i][0] += delx*fpair;
101         f[i][1] += dely*fpair;
102         f[i][2] += delz*fpair;
103         if (newton_pair || j < nlocal) {
104           f[j][0] -= delx*fpair;
105           f[j][1] -= dely*fpair;
106           f[j][2] -= delz*fpair;
107         }
108 
109         if (eflag) {
110           if (rsq < cut_coulsq[itype][jtype])
111             ecoul = factor_coul * qqrd2e * qtmp*q[j] * rinv * screening;
112           else ecoul = 0.0;
113           if (rsq < cut_ljsq[itype][jtype]) {
114             evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
115               offset[itype][jtype];
116             evdwl *= factor_lj;
117           } else evdwl = 0.0;
118         }
119 
120         if (evflag) ev_tally(i,j,nlocal,newton_pair,
121                              evdwl,ecoul,fpair,delx,dely,delz);
122       }
123     }
124   }
125 
126   if (vflag_fdotr) virial_fdotr_compute();
127 }
128 
129 /* ----------------------------------------------------------------------
130    global settings
131 ------------------------------------------------------------------------- */
132 
settings(int narg,char ** arg)133 void PairLJCutCoulDebye::settings(int narg, char **arg)
134 {
135   if (narg < 2 || narg > 3) error->all(FLERR,"Illegal pair_style command");
136 
137   kappa = utils::numeric(FLERR,arg[0],false,lmp);
138   cut_lj_global = utils::numeric(FLERR,arg[1],false,lmp);
139   if (narg == 2) cut_coul_global = cut_lj_global;
140   else cut_coul_global = utils::numeric(FLERR,arg[2],false,lmp);
141 
142   // reset cutoffs that were previously set from data file
143 
144   if (allocated) {
145     int i,j;
146     for (i = 1; i <= atom->ntypes; i++)
147       for (j = i+1; j <= atom->ntypes; j++)
148         if (setflag[i][j] == 1) {
149           cut_lj[i][j] = cut_lj_global;
150           cut_coul[i][j] = cut_coul_global;
151         }
152   }
153 }
154 
155 /* ----------------------------------------------------------------------
156    proc 0 writes to restart file
157 ------------------------------------------------------------------------- */
158 
write_restart_settings(FILE * fp)159 void PairLJCutCoulDebye::write_restart_settings(FILE *fp)
160 {
161   fwrite(&cut_lj_global,sizeof(double),1,fp);
162   fwrite(&cut_coul_global,sizeof(double),1,fp);
163   fwrite(&kappa,sizeof(double),1,fp);
164   fwrite(&offset_flag,sizeof(int),1,fp);
165   fwrite(&mix_flag,sizeof(int),1,fp);
166 }
167 
168 /* ----------------------------------------------------------------------
169    proc 0 reads from restart file, bcasts
170 ------------------------------------------------------------------------- */
171 
read_restart_settings(FILE * fp)172 void PairLJCutCoulDebye::read_restart_settings(FILE *fp)
173 {
174   if (comm->me == 0) {
175     utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error);
176     utils::sfread(FLERR,&cut_coul_global,sizeof(double),1,fp,nullptr,error);
177     utils::sfread(FLERR,&kappa,sizeof(double),1,fp,nullptr,error);
178     utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
179     utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
180   }
181   MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
182   MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
183   MPI_Bcast(&kappa,1,MPI_DOUBLE,0,world);
184   MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
185   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
186 }
187 
188 /* ---------------------------------------------------------------------- */
189 
single(int i,int j,int itype,int jtype,double rsq,double factor_coul,double factor_lj,double & fforce)190 double PairLJCutCoulDebye::single(int i, int j, int itype, int jtype,
191                                   double rsq,
192                                   double factor_coul, double factor_lj,
193                                   double &fforce)
194 {
195   double r2inv,r6inv,r,rinv,screening,forcecoul,forcelj,phicoul,philj;
196 
197   r2inv = 1.0/rsq;
198   if (rsq < cut_coulsq[itype][jtype]) {
199     r = sqrt(rsq);
200     rinv = 1.0/r;
201     screening = exp(-kappa*r);
202     forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] *
203       screening * (kappa + rinv);
204   } else forcecoul = 0.0;
205   if (rsq < cut_ljsq[itype][jtype]) {
206     r6inv = r2inv*r2inv*r2inv;
207     forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
208   } else forcelj = 0.0;
209   fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
210 
211   double eng = 0.0;
212   if (rsq < cut_coulsq[itype][jtype]) {
213     phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * screening;
214     eng += factor_coul*phicoul;
215   }
216   if (rsq < cut_ljsq[itype][jtype]) {
217     philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
218       offset[itype][jtype];
219     eng += factor_lj*philj;
220   }
221 
222   return eng;
223 }
224