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 "bond_fene_expand.h"
16 
17 #include "atom.h"
18 #include "comm.h"
19 #include "error.h"
20 #include "force.h"
21 #include "math_const.h"
22 #include "memory.h"
23 #include "neighbor.h"
24 #include "update.h"
25 
26 #include <cmath>
27 
28 using namespace LAMMPS_NS;
29 using MathConst::MY_CUBEROOT2;
30 
31 /* ---------------------------------------------------------------------- */
32 
~BondFENEExpand()33 BondFENEExpand::~BondFENEExpand()
34 {
35   if (allocated) {
36     memory->destroy(setflag);
37     memory->destroy(k);
38     memory->destroy(r0);
39     memory->destroy(epsilon);
40     memory->destroy(sigma);
41     memory->destroy(shift);
42   }
43 }
44 
45 /* ---------------------------------------------------------------------- */
46 
compute(int eflag,int vflag)47 void BondFENEExpand::compute(int eflag, int vflag)
48 {
49   int i1,i2,n,type;
50   double delx,dely,delz,ebond,fbond;
51   double rsq,r0sq,rlogarg,sr2,sr6;
52   double r,rshift,rshiftsq;
53 
54   ebond = sr6 = 0.0;
55   ev_init(eflag,vflag);
56 
57   double **x = atom->x;
58   double **f = atom->f;
59   int **bondlist = neighbor->bondlist;
60   int nbondlist = neighbor->nbondlist;
61   int nlocal = atom->nlocal;
62   int newton_bond = force->newton_bond;
63 
64   for (n = 0; n < nbondlist; n++) {
65     i1 = bondlist[n][0];
66     i2 = bondlist[n][1];
67     type = bondlist[n][2];
68 
69     delx = x[i1][0] - x[i2][0];
70     dely = x[i1][1] - x[i2][1];
71     delz = x[i1][2] - x[i2][2];
72 
73     // force from log term
74 
75     rsq = delx*delx + dely*dely + delz*delz;
76     r = sqrt(rsq);
77     rshift = r - shift[type];
78     rshiftsq = rshift*rshift;
79     r0sq = r0[type] * r0[type];
80     rlogarg = 1.0 - rshiftsq/r0sq;
81 
82     // if r -> r0, then rlogarg < 0.0 which is an error
83     // issue a warning and reset rlogarg = epsilon
84     // if r > 2*r0 something serious is wrong, abort
85 
86     if (rlogarg < 0.1) {
87       error->warning(FLERR,"FENE bond too long: {} {} {} {:.8}",
88                      update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq));
89       if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
90       rlogarg = 0.1;
91     }
92 
93     fbond = -k[type]*rshift/rlogarg/r;
94 
95     // force from LJ term
96 
97     if (rshiftsq < MY_CUBEROOT2*sigma[type]*sigma[type]) {
98       sr2 = sigma[type]*sigma[type]/rshiftsq;
99       sr6 = sr2*sr2*sr2;
100       fbond += 48.0*epsilon[type]*sr6*(sr6-0.5)/rshift/r;
101     }
102 
103     // energy
104 
105     if (eflag) {
106       ebond = -0.5 * k[type]*r0sq*log(rlogarg);
107       if (rshiftsq < MY_CUBEROOT2*sigma[type]*sigma[type])
108         ebond += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
109     }
110 
111     // apply force to each of 2 atoms
112 
113     if (newton_bond || i1 < nlocal) {
114       f[i1][0] += delx*fbond;
115       f[i1][1] += dely*fbond;
116       f[i1][2] += delz*fbond;
117     }
118 
119     if (newton_bond || i2 < nlocal) {
120       f[i2][0] -= delx*fbond;
121       f[i2][1] -= dely*fbond;
122       f[i2][2] -= delz*fbond;
123     }
124 
125     if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
126   }
127 }
128 
129 /* ---------------------------------------------------------------------- */
130 
allocate()131 void BondFENEExpand::allocate()
132 {
133   allocated = 1;
134   int n = atom->nbondtypes;
135 
136   memory->create(k,n+1,"bond:k");
137   memory->create(r0,n+1,"bond:r0");
138   memory->create(epsilon,n+1,"bond:epsilon");
139   memory->create(sigma,n+1,"bond:sigma");
140   memory->create(shift,n+1,"bond:shift");
141   memory->create(setflag,n+1,"bond:setflag");
142   for (int i = 1; i <= n; i++) setflag[i] = 0;
143 }
144 
145 /* ----------------------------------------------------------------------
146    set coeffs for one type
147 ------------------------------------------------------------------------- */
148 
coeff(int narg,char ** arg)149 void BondFENEExpand::coeff(int narg, char **arg)
150 {
151   if (narg != 6) error->all(FLERR,"Incorrect args for bond coefficients");
152   if (!allocated) allocate();
153 
154   int ilo,ihi;
155   utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error);
156 
157   double k_one = utils::numeric(FLERR,arg[1],false,lmp);
158   double r0_one = utils::numeric(FLERR,arg[2],false,lmp);
159   double epsilon_one = utils::numeric(FLERR,arg[3],false,lmp);
160   double sigma_one = utils::numeric(FLERR,arg[4],false,lmp);
161   double shift_one = utils::numeric(FLERR,arg[5],false,lmp);
162 
163   int count = 0;
164   for (int i = ilo; i <= ihi; i++) {
165     k[i] = k_one;
166     r0[i] = r0_one;
167     epsilon[i] = epsilon_one;
168     sigma[i] = sigma_one;
169     shift[i] = shift_one;
170     setflag[i] = 1;
171     count++;
172   }
173 
174   if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
175 }
176 
177 /* ----------------------------------------------------------------------
178    check if special_bond settings are valid
179 ------------------------------------------------------------------------- */
180 
init_style()181 void BondFENEExpand::init_style()
182 {
183   // special bonds should be 0 1 1
184 
185   if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 ||
186       force->special_lj[3] != 1.0) {
187     if (comm->me == 0)
188       error->warning(FLERR,"Use special bonds = 0,1,1 with bond style fene/expand");
189   }
190 }
191 
192 /* ---------------------------------------------------------------------- */
193 
equilibrium_distance(int i)194 double BondFENEExpand::equilibrium_distance(int i)
195 {
196   return 0.97*sigma[i] + shift[i];
197 }
198 
199 /* ----------------------------------------------------------------------
200    proc 0 writes to restart file
201 ------------------------------------------------------------------------- */
202 
write_restart(FILE * fp)203 void BondFENEExpand::write_restart(FILE *fp)
204 {
205   fwrite(&k[1],sizeof(double),atom->nbondtypes,fp);
206   fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);
207   fwrite(&epsilon[1],sizeof(double),atom->nbondtypes,fp);
208   fwrite(&sigma[1],sizeof(double),atom->nbondtypes,fp);
209   fwrite(&shift[1],sizeof(double),atom->nbondtypes,fp);
210 }
211 
212 /* ----------------------------------------------------------------------
213    proc 0 reads from restart file, bcasts
214 ------------------------------------------------------------------------- */
215 
read_restart(FILE * fp)216 void BondFENEExpand::read_restart(FILE *fp)
217 {
218   allocate();
219 
220   if (comm->me == 0) {
221     utils::sfread(FLERR,&k[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
222     utils::sfread(FLERR,&r0[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
223     utils::sfread(FLERR,&epsilon[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
224     utils::sfread(FLERR,&sigma[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
225     utils::sfread(FLERR,&shift[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
226   }
227   MPI_Bcast(&k[1],atom->nbondtypes,MPI_DOUBLE,0,world);
228   MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world);
229   MPI_Bcast(&epsilon[1],atom->nbondtypes,MPI_DOUBLE,0,world);
230   MPI_Bcast(&sigma[1],atom->nbondtypes,MPI_DOUBLE,0,world);
231   MPI_Bcast(&shift[1],atom->nbondtypes,MPI_DOUBLE,0,world);
232 
233   for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
234 }
235 
236 /* ----------------------------------------------------------------------
237    proc 0 writes to data file
238 ------------------------------------------------------------------------- */
239 
write_data(FILE * fp)240 void BondFENEExpand::write_data(FILE *fp)
241 {
242   for (int i = 1; i <= atom->nbondtypes; i++)
243     fprintf(fp,"%d %g %g %g %g %g\n",i,k[i],r0[i],epsilon[i],sigma[i],shift[i]);
244 }
245 
246 /* ---------------------------------------------------------------------- */
247 
single(int type,double rsq,int,int,double & fforce)248 double BondFENEExpand::single(int type, double rsq, int /*i*/, int /*j*/,
249                         double &fforce)
250 {
251   double r = sqrt(rsq);
252   double rshift = r - shift[type];
253   double rshiftsq = rshift*rshift;
254   double r0sq = r0[type] * r0[type];
255   double rlogarg = 1.0 - rshiftsq/r0sq;
256 
257   // if r -> r0, then rlogarg < 0.0 which is an error
258   // issue a warning and reset rlogarg = epsilon
259   // if r > 2*r0 something serious is wrong, abort
260 
261   if (rlogarg < 0.1) {
262     error->warning(FLERR,"FENE bond too long: {} {:.8}",
263                      update->ntimestep,sqrt(rsq));
264     if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
265     rlogarg = 0.1;
266   }
267 
268   double eng = -0.5 * k[type]*r0sq*log(rlogarg);
269   fforce = -k[type]*rshift/rlogarg/r;
270   if (rshiftsq < MY_CUBEROOT2*sigma[type]*sigma[type]) {
271     double sr2,sr6;
272     sr2 = sigma[type]*sigma[type]/rshiftsq;
273     sr6 = sr2*sr2*sr2;
274     eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
275     fforce += 48.0*epsilon[type]*sr6*(sr6-0.5)/rshift/r;
276   }
277 
278   return eng;
279 }
280