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