1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13 /* ----------------------------------------------------------------------
14 Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
15 ------------------------------------------------------------------------- */
16
17 #include "bond_oxdna_fene.h"
18
19 #include "atom.h"
20 #include "comm.h"
21 #include "error.h"
22 #include "force.h"
23 #include "memory.h"
24 #include "neighbor.h"
25 #include "update.h"
26
27 #include "atom_vec_ellipsoid.h"
28 #include "math_extra.h"
29
30 #include <cmath>
31
32 using namespace LAMMPS_NS;
33
34 /* ---------------------------------------------------------------------- */
35
~BondOxdnaFene()36 BondOxdnaFene::~BondOxdnaFene()
37 {
38 if (allocated) {
39 memory->destroy(setflag);
40 memory->destroy(k);
41 memory->destroy(Delta);
42 memory->destroy(r0);
43 }
44 }
45
46 /* ----------------------------------------------------------------------
47 compute vector COM-sugar-phosphate backbone interaction site in oxDNA
48 ------------------------------------------------------------------------- */
compute_interaction_sites(double e1[3],double[3],double[3],double r[3]) const49 void BondOxdnaFene::compute_interaction_sites(double e1[3], double /*e2*/[3], double /*e3*/[3],
50 double r[3]) const
51 {
52 constexpr double d_cs = -0.4;
53
54 r[0] = d_cs * e1[0];
55 r[1] = d_cs * e1[1];
56 r[2] = d_cs * e1[2];
57 }
58
59 /* ----------------------------------------------------------------------
60 tally energy and virial into global and per-atom accumulators
61 ------------------------------------------------------------------------- */
62
ev_tally_xyz(int i,int j,int nlocal,int newton_bond,double ebond,double fx,double fy,double fz,double delx,double dely,double delz)63 void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, double ebond, double fx,
64 double fy, double fz, double delx, double dely, double delz)
65 {
66 double ebondhalf, v[6];
67
68 if (eflag_either) {
69 if (eflag_global) {
70 if (newton_bond)
71 energy += ebond;
72 else {
73 ebondhalf = 0.5 * ebond;
74 if (i < nlocal) energy += ebondhalf;
75 if (j < nlocal) energy += ebondhalf;
76 }
77 }
78 if (eflag_atom) {
79 ebondhalf = 0.5 * ebond;
80 if (newton_bond || i < nlocal) eatom[i] += ebondhalf;
81 if (newton_bond || j < nlocal) eatom[j] += ebondhalf;
82 }
83 }
84
85 if (vflag_either) {
86 v[0] = delx * fx;
87 v[1] = dely * fy;
88 v[2] = delz * fz;
89 v[3] = delx * fy;
90 v[4] = delx * fz;
91 v[5] = dely * fz;
92
93 if (vflag_global) {
94 if (newton_bond) {
95 virial[0] += v[0];
96 virial[1] += v[1];
97 virial[2] += v[2];
98 virial[3] += v[3];
99 virial[4] += v[4];
100 virial[5] += v[5];
101 } else {
102 if (i < nlocal) {
103 virial[0] += 0.5 * v[0];
104 virial[1] += 0.5 * v[1];
105 virial[2] += 0.5 * v[2];
106 virial[3] += 0.5 * v[3];
107 virial[4] += 0.5 * v[4];
108 virial[5] += 0.5 * v[5];
109 }
110 if (j < nlocal) {
111 virial[0] += 0.5 * v[0];
112 virial[1] += 0.5 * v[1];
113 virial[2] += 0.5 * v[2];
114 virial[3] += 0.5 * v[3];
115 virial[4] += 0.5 * v[4];
116 virial[5] += 0.5 * v[5];
117 }
118 }
119 }
120
121 if (vflag_atom) {
122 if (newton_bond || i < nlocal) {
123 vatom[i][0] += 0.5 * v[0];
124 vatom[i][1] += 0.5 * v[1];
125 vatom[i][2] += 0.5 * v[2];
126 vatom[i][3] += 0.5 * v[3];
127 vatom[i][4] += 0.5 * v[4];
128 vatom[i][5] += 0.5 * v[5];
129 }
130 if (newton_bond || j < nlocal) {
131 vatom[j][0] += 0.5 * v[0];
132 vatom[j][1] += 0.5 * v[1];
133 vatom[j][2] += 0.5 * v[2];
134 vatom[j][3] += 0.5 * v[3];
135 vatom[j][4] += 0.5 * v[4];
136 vatom[j][5] += 0.5 * v[5];
137 }
138 }
139 }
140 }
141
142 /* ----------------------------------------------------------------------
143 compute function for oxDNA FENE-bond interaction
144 s=sugar-phosphate backbone site, b=base site, st=stacking site
145 ------------------------------------------------------------------------- */
compute(int eflag,int vflag)146 void BondOxdnaFene::compute(int eflag, int vflag)
147 {
148 int a, b, in, type;
149 double delf[3], delta[3], deltb[3]; // force, torque increment;;
150 double delr[3], ebond, fbond;
151 double rsq, Deltasq, rlogarg;
152 double r, rr0, rr0sq;
153 // vectors COM-backbone site in lab frame
154 double ra_cs[3], rb_cs[3];
155
156 double *qa, ax[3], ay[3], az[3];
157 double *qb, bx[3], by[3], bz[3];
158
159 double **x = atom->x;
160 double **f = atom->f;
161 double **torque = atom->torque;
162
163 AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
164 AtomVecEllipsoid::Bonus *bonus = avec->bonus;
165 int *ellipsoid = atom->ellipsoid;
166
167 int **bondlist = neighbor->bondlist;
168 int nbondlist = neighbor->nbondlist;
169 int nlocal = atom->nlocal;
170 int newton_bond = force->newton_bond;
171
172 ebond = 0.0;
173 ev_init(eflag, vflag);
174
175 // loop over FENE bonds
176
177 for (in = 0; in < nbondlist; in++) {
178
179 a = bondlist[in][1];
180 b = bondlist[in][0];
181 type = bondlist[in][2];
182
183 qa = bonus[ellipsoid[a]].quat;
184 MathExtra::q_to_exyz(qa, ax, ay, az);
185 qb = bonus[ellipsoid[b]].quat;
186 MathExtra::q_to_exyz(qb, bx, by, bz);
187
188 // vector COM-backbone site a and b
189 compute_interaction_sites(ax, ay, az, ra_cs);
190 compute_interaction_sites(bx, by, bz, rb_cs);
191
192 // vector backbone site b to a
193 delr[0] = x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0];
194 delr[1] = x[a][1] + ra_cs[1] - x[b][1] - rb_cs[1];
195 delr[2] = x[a][2] + ra_cs[2] - x[b][2] - rb_cs[2];
196 rsq = delr[0] * delr[0] + delr[1] * delr[1] + delr[2] * delr[2];
197 r = sqrt(rsq);
198
199 rr0 = r - r0[type];
200 rr0sq = rr0 * rr0;
201 Deltasq = Delta[type] * Delta[type];
202 rlogarg = 1.0 - rr0sq / Deltasq;
203
204 // if r -> Delta, then rlogarg < 0.0 which is an error
205 // issue a warning and reset rlogarg = epsilon
206 // if r > 2*Delta something serious is wrong, abort
207
208 if (rlogarg < 0.1) {
209 error->warning(FLERR, "FENE bond too long: {} {} {} {}", update->ntimestep, atom->tag[a],
210 atom->tag[b], r);
211 rlogarg = 0.1;
212 }
213
214 fbond = -k[type] * rr0 / rlogarg / Deltasq / r;
215 delf[0] = delr[0] * fbond;
216 delf[1] = delr[1] * fbond;
217 delf[2] = delr[2] * fbond;
218
219 // energy
220
221 if (eflag) { ebond = -0.5 * k[type] * log(rlogarg); }
222
223 // apply force and torque to each of 2 atoms
224
225 if (newton_bond || a < nlocal) {
226
227 f[a][0] += delf[0];
228 f[a][1] += delf[1];
229 f[a][2] += delf[2];
230
231 MathExtra::cross3(ra_cs, delf, delta);
232
233 torque[a][0] += delta[0];
234 torque[a][1] += delta[1];
235 torque[a][2] += delta[2];
236 }
237
238 if (newton_bond || b < nlocal) {
239
240 f[b][0] -= delf[0];
241 f[b][1] -= delf[1];
242 f[b][2] -= delf[2];
243
244 MathExtra::cross3(rb_cs, delf, deltb);
245
246 torque[b][0] -= deltb[0];
247 torque[b][1] -= deltb[1];
248 torque[b][2] -= deltb[2];
249 }
250
251 // increment energy and virial
252 // NOTE: The virial is calculated on the 'molecular' basis.
253 // (see G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345-392 (1986))
254
255 if (evflag)
256 ev_tally_xyz(a, b, nlocal, newton_bond, ebond, delf[0], delf[1], delf[2], x[a][0] - x[b][0],
257 x[a][1] - x[b][1], x[a][2] - x[b][2]);
258 }
259 }
260
261 /* ---------------------------------------------------------------------- */
262
allocate()263 void BondOxdnaFene::allocate()
264 {
265 allocated = 1;
266 int n = atom->nbondtypes;
267
268 memory->create(k, n + 1, "bond:k");
269 memory->create(Delta, n + 1, "bond:Delta");
270 memory->create(r0, n + 1, "bond:r0");
271 memory->create(setflag, n + 1, "bond:setflag");
272
273 for (int i = 1; i <= n; i++) setflag[i] = 0;
274 }
275
276 /* ----------------------------------------------------------------------
277 set coeffs for one type
278 ------------------------------------------------------------------------- */
279
coeff(int narg,char ** arg)280 void BondOxdnaFene::coeff(int narg, char **arg)
281 {
282 if (narg != 4) error->all(FLERR, "Incorrect args for bond coefficients in oxdna/fene");
283 if (!allocated) allocate();
284
285 int ilo, ihi;
286 utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
287
288 double k_one = utils::numeric(FLERR, arg[1], false, lmp);
289 double Delta_one = utils::numeric(FLERR, arg[2], false, lmp);
290 double r0_one = utils::numeric(FLERR, arg[3], false, lmp);
291
292 int count = 0;
293
294 for (int i = ilo; i <= ihi; i++) {
295 k[i] = k_one;
296 Delta[i] = Delta_one;
297 r0[i] = r0_one;
298 setflag[i] = 1;
299 count++;
300 }
301
302 if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients in oxdna/fene");
303 }
304
305 /* ----------------------------------------------------------------------
306 set special_bond settings and check if valid
307 ------------------------------------------------------------------------- */
308
init_style()309 void BondOxdnaFene::init_style()
310 {
311 if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0)
312 error->all(
313 FLERR,
314 "Must use 'special_bonds lj 0 1 1' with bond style oxdna/fene, oxdna2/fene or oxrna2/fene");
315 }
316
317 /* ---------------------------------------------------------------------- */
318
equilibrium_distance(int i)319 double BondOxdnaFene::equilibrium_distance(int i)
320 {
321 return r0[i];
322 }
323
324 /* ----------------------------------------------------------------------
325 proc 0 writes to restart file
326 ------------------------------------------------------------------------- */
327
write_restart(FILE * fp)328 void BondOxdnaFene::write_restart(FILE *fp)
329 {
330 fwrite(&k[1], sizeof(double), atom->nbondtypes, fp);
331 fwrite(&Delta[1], sizeof(double), atom->nbondtypes, fp);
332 fwrite(&r0[1], sizeof(double), atom->nbondtypes, fp);
333 }
334
335 /* ----------------------------------------------------------------------
336 proc 0 reads from restart file, bcasts
337 ------------------------------------------------------------------------- */
338
read_restart(FILE * fp)339 void BondOxdnaFene::read_restart(FILE *fp)
340 {
341 allocate();
342
343 if (comm->me == 0) {
344 utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
345 utils::sfread(FLERR, &Delta[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
346 utils::sfread(FLERR, &r0[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
347 }
348 MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
349 MPI_Bcast(&Delta[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
350 MPI_Bcast(&r0[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
351
352 for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
353 }
354
355 /* ----------------------------------------------------------------------
356 proc 0 writes to data file
357 ------------------------------------------------------------------------- */
358
write_data(FILE * fp)359 void BondOxdnaFene::write_data(FILE *fp)
360 {
361 for (int i = 1; i <= atom->nbondtypes; i++)
362 fprintf(fp, "%d %g %g %g\n", i, k[i], r0[i], Delta[i]);
363 }
364
365 /* ---------------------------------------------------------------------- */
366
single(int type,double rsq,int,int,double & fforce)367 double BondOxdnaFene::single(int type, double rsq, int /*i*/, int /*j*/, double &fforce)
368 {
369 double r = sqrt(rsq);
370 double rr0 = r - r0[type];
371 double rr0sq = rr0 * rr0;
372 double Deltasq = Delta[type] * Delta[type];
373 double rlogarg = 1.0 - rr0sq / Deltasq;
374
375 // if r -> Delta, then rlogarg < 0.0 which is an error
376 // issue a warning and reset rlogarg = epsilon
377 // if r > 2*Delta something serious is wrong, abort
378
379 if (rlogarg < 0.1) {
380 error->warning(FLERR, "FENE bond too long: {} {:.8}", update->ntimestep, sqrt(rsq));
381 rlogarg = 0.1;
382 }
383
384 double eng = -0.5 * k[type] * log(rlogarg);
385 fforce = -k[type] * rr0 / rlogarg / Deltasq / r;
386
387 return eng;
388 }
389