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 /* ----------------------------------------------------------------------
15 Contributing authors: Mario Orsi & Wei Ding (QMUL), m.orsi@qmul.ac.uk
16 ------------------------------------------------------------------------- */
17
18 #include "angle_dipole.h"
19
20 #include "atom.h"
21 #include "comm.h"
22 #include "domain.h"
23 #include "error.h"
24 #include "force.h"
25 #include "math_const.h"
26 #include "memory.h"
27 #include "neighbor.h"
28
29 #include <cmath>
30
31 using namespace LAMMPS_NS;
32 using namespace MathConst;
33
34 /* ---------------------------------------------------------------------- */
35
AngleDipole(LAMMPS * lmp)36 AngleDipole::AngleDipole(LAMMPS *lmp) : Angle(lmp)
37 {
38 k = nullptr;
39 gamma0 = nullptr;
40 }
41
42 /* ---------------------------------------------------------------------- */
43
~AngleDipole()44 AngleDipole::~AngleDipole()
45 {
46 if (allocated) {
47 memory->destroy(setflag);
48 memory->destroy(k);
49 memory->destroy(gamma0);
50 }
51 }
52
53 /* ---------------------------------------------------------------------- */
54
compute(int eflag,int vflag)55 void AngleDipole::compute(int eflag, int vflag)
56 {
57 int iRef, iDip, iDummy, n, type;
58 double delx, dely, delz;
59 double eangle, tangle, fi[3], fj[3];
60 double r, cosGamma, deltaGamma, kdg, rmu;
61
62 eangle = 0.0;
63 ev_init(eflag, vflag);
64
65 double **x = atom->x; // position vector
66 double **mu = atom->mu; // point-dipole components and moment magnitude
67 double **torque = atom->torque;
68 int **anglelist = neighbor->anglelist;
69 int nanglelist = neighbor->nanglelist;
70 int nlocal = atom->nlocal;
71 int newton_bond = force->newton_bond;
72
73 double **f = atom->f;
74 double delTx, delTy, delTz;
75 double fx, fy, fz, fmod, fmod_sqrtff;
76
77 if (!newton_bond) error->all(FLERR, "'newton' flag for bonded interactions must be 'on'");
78
79 for (n = 0; n < nanglelist; n++) {
80 iDip = anglelist[n][0]; // dipole whose orientation is to be restrained
81 iRef = anglelist[n][1]; // reference atom toward which dipole will point
82 iDummy = anglelist[n][2]; // dummy atom - irrelevant to the interaction
83 type = anglelist[n][3];
84
85 delx = x[iRef][0] - x[iDip][0];
86 dely = x[iRef][1] - x[iDip][1];
87 delz = x[iRef][2] - x[iDip][2];
88
89 r = sqrt(delx * delx + dely * dely + delz * delz);
90
91 rmu = r * mu[iDip][3];
92 cosGamma = (mu[iDip][0] * delx + mu[iDip][1] * dely + mu[iDip][2] * delz) / rmu;
93 deltaGamma = cosGamma - cos(gamma0[type]);
94 kdg = k[type] * deltaGamma;
95
96 if (eflag) eangle = kdg * deltaGamma; // energy
97
98 tangle = 2.0 * kdg / rmu;
99
100 delTx = tangle * (dely * mu[iDip][2] - delz * mu[iDip][1]);
101 delTy = tangle * (delz * mu[iDip][0] - delx * mu[iDip][2]);
102 delTz = tangle * (delx * mu[iDip][1] - dely * mu[iDip][0]);
103
104 torque[iDip][0] += delTx;
105 torque[iDip][1] += delTy;
106 torque[iDip][2] += delTz;
107
108 // Force couple that counterbalances dipolar torque
109 fx = dely * delTz - delz * delTy; // direction (fi): - r x (-T)
110 fy = delz * delTx - delx * delTz;
111 fz = delx * delTy - dely * delTx;
112
113 fmod = sqrt(delTx * delTx + delTy * delTy + delTz * delTz) / r; // magnitude
114 fmod_sqrtff = fmod / sqrt(fx * fx + fy * fy + fz * fz);
115
116 fi[0] = fx * fmod_sqrtff;
117 fi[1] = fy * fmod_sqrtff;
118 fi[2] = fz * fmod_sqrtff;
119
120 fj[0] = -fi[0];
121 fj[1] = -fi[1];
122 fj[2] = -fi[2];
123
124 f[iDip][0] += fj[0];
125 f[iDip][1] += fj[1];
126 f[iDip][2] += fj[2];
127
128 f[iRef][0] += fi[0];
129 f[iRef][1] += fi[1];
130 f[iRef][2] += fi[2];
131
132 if (evflag) // virial = rij.fi = 0 (fj = -fi & fk = 0)
133 ev_tally(iRef, iDip, iDummy, nlocal, newton_bond, eangle, fj, fi, 0.0, 0.0, 0.0, 0.0, 0.0,
134 0.0);
135 }
136 }
137
138 /* ---------------------------------------------------------------------- */
139
allocate()140 void AngleDipole::allocate()
141 {
142 allocated = 1;
143 int n = atom->nangletypes;
144
145 memory->create(k, n + 1, "angle:k");
146 memory->create(gamma0, n + 1, "angle:gamma0");
147
148 memory->create(setflag, n + 1, "angle:setflag");
149 for (int i = 1; i <= n; i++) setflag[i] = 0;
150 }
151
152 /* ----------------------------------------------------------------------
153 set coeffs for one or more types
154 ------------------------------------------------------------------------- */
155
coeff(int narg,char ** arg)156 void AngleDipole::coeff(int narg, char **arg)
157 {
158 if (narg != 3) error->all(FLERR, "Incorrect args for angle coefficients");
159 if (!allocated) allocate();
160
161 int ilo, ihi;
162 utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
163
164 double k_one = utils::numeric(FLERR, arg[1], false, lmp);
165 double gamma0_one = utils::numeric(FLERR, arg[2], false, lmp);
166
167 // convert gamma0 from degrees to radians
168
169 int count = 0;
170 for (int i = ilo; i <= ihi; i++) {
171 k[i] = k_one;
172 gamma0[i] = gamma0_one / 180.0 * MY_PI;
173 setflag[i] = 1;
174 count++;
175 }
176
177 if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
178 }
179
180 /* ----------------------------------------------------------------------
181 set coeffs for one or more types
182 ------------------------------------------------------------------------- */
183
init_style()184 void AngleDipole::init_style()
185 {
186 if (!atom->mu_flag || !atom->torque_flag)
187 error->all(FLERR, "Angle style dipole requires atom attributes mu and torque");
188 }
189
190 /* ----------------------------------------------------------------------
191 used by SHAKE
192 ------------------------------------------------------------------------- */
193
equilibrium_angle(int i)194 double AngleDipole::equilibrium_angle(int i)
195 {
196 return gamma0[i];
197 }
198
199 /* ----------------------------------------------------------------------
200 proc 0 writes out coeffs to restart file
201 ------------------------------------------------------------------------- */
202
write_restart(FILE * fp)203 void AngleDipole::write_restart(FILE *fp)
204 {
205 fwrite(&k[1], sizeof(double), atom->nangletypes, fp);
206 fwrite(&gamma0[1], sizeof(double), atom->nangletypes, fp);
207 }
208
209 /* ----------------------------------------------------------------------
210 proc 0 reads coeffs from restart file, bcasts them
211 ------------------------------------------------------------------------- */
212
read_restart(FILE * fp)213 void AngleDipole::read_restart(FILE *fp)
214 {
215 allocate();
216
217 if (comm->me == 0) {
218 utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
219 utils::sfread(FLERR, &gamma0[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
220 }
221 MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world);
222 MPI_Bcast(&gamma0[1], atom->nangletypes, MPI_DOUBLE, 0, world);
223
224 for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
225 }
226
227 /* ----------------------------------------------------------------------
228 proc 0 writes to data file
229 ------------------------------------------------------------------------- */
230
write_data(FILE * fp)231 void AngleDipole::write_data(FILE *fp)
232 {
233 for (int i = 1; i <= atom->nangletypes; i++) fprintf(fp, "%d %g %g\n", i, k[i], gamma0[i]);
234 }
235
236 /* ----------------------------------------------------------------------
237 used by ComputeAngleLocal
238 ------------------------------------------------------------------------- */
239
single(int type,int iRef,int iDip,int)240 double AngleDipole::single(int type, int iRef, int iDip, int /*iDummy*/)
241 {
242 double **x = atom->x; // position vector
243 double **mu = atom->mu; // point-dipole components and moment magnitude
244
245 double delx = x[iRef][0] - x[iDip][0];
246 double dely = x[iRef][1] - x[iDip][1];
247 double delz = x[iRef][2] - x[iDip][2];
248
249 domain->minimum_image(delx, dely, delz);
250
251 double r = sqrt(delx * delx + dely * dely + delz * delz);
252 double rmu = r * mu[iDip][3];
253 double cosGamma = (mu[iDip][0] * delx + mu[iDip][1] * dely + mu[iDip][2] * delz) / rmu;
254 double deltaGamma = cosGamma - cos(gamma0[type]);
255 double kdg = k[type] * deltaGamma;
256
257 return kdg * deltaGamma; // energy
258 }
259