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