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