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 /* ----------------------------------------------------------------------
16    Contributing author: Steven Vandenbrande
17 ------------------------------------------------------------------------- */
18 
19 #include "angle_cross.h"
20 
21 #include <cmath>
22 #include "atom.h"
23 #include "neighbor.h"
24 #include "domain.h"
25 #include "comm.h"
26 #include "force.h"
27 #include "math_const.h"
28 #include "memory.h"
29 #include "error.h"
30 
31 
32 using namespace LAMMPS_NS;
33 using namespace MathConst;
34 
35 #define SMALL 0.001
36 
37 /* ---------------------------------------------------------------------- */
38 
AngleCross(LAMMPS * lmp)39 AngleCross::AngleCross(LAMMPS *lmp) : Angle(lmp) {}
40 
41 /* ---------------------------------------------------------------------- */
42 
~AngleCross()43 AngleCross::~AngleCross()
44 {
45   if (copymode) return;
46 
47   if (allocated) {
48     memory->destroy(setflag);
49     memory->destroy(kss);
50     memory->destroy(kbs0);
51     memory->destroy(kbs1);
52     memory->destroy(r00);
53     memory->destroy(r01);
54     memory->destroy(theta0);
55   }
56 }
57 
58 /* ---------------------------------------------------------------------- */
59 
compute(int eflag,int vflag)60 void AngleCross::compute(int eflag, int vflag)
61 {
62   int i1,i2,i3,n,type;
63   double delx1,dely1,delz1,delx2,dely2,delz2;
64   double eangle,f1[3],f3[3];
65   double dtheta;
66   double dr1,dr2,tk1,tk2,aa1,aa2,aa11,aa12,aa21,aa22;
67   double rsq1,rsq2,r1,r2,c,s,b1,b2;
68   double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22;
69 
70   eangle = 0.0;
71   ev_init(eflag,vflag);
72 
73   double **x = atom->x;
74   double **f = atom->f;
75   int **anglelist = neighbor->anglelist;
76   int nanglelist = neighbor->nanglelist;
77   int nlocal = atom->nlocal;
78   int newton_bond = force->newton_bond;
79 
80   for (n = 0; n < nanglelist; n++) {
81     i1 = anglelist[n][0];
82     i2 = anglelist[n][1];
83     i3 = anglelist[n][2];
84     type = anglelist[n][3];
85 
86     // 1st bond
87 
88     delx1 = x[i1][0] - x[i2][0];
89     dely1 = x[i1][1] - x[i2][1];
90     delz1 = x[i1][2] - x[i2][2];
91 
92     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
93     r1 = sqrt(rsq1);
94 
95     // 2nd bond
96 
97     delx2 = x[i3][0] - x[i2][0];
98     dely2 = x[i3][1] - x[i2][1];
99     delz2 = x[i3][2] - x[i2][2];
100 
101     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
102     r2 = sqrt(rsq2);
103 
104     // angle (cos and sin)
105 
106     c = delx1*delx2 + dely1*dely2 + delz1*delz2;
107     c /= r1*r2;
108 
109     if (c > 1.0) c = 1.0;
110     if (c < -1.0) c = -1.0;
111 
112     s = sqrt(1.0 - c*c);
113     if (s < SMALL) s = SMALL;
114     s = 1.0/s;
115 
116     // force & energy for bond-bond term
117     dr1 = r1 - r00[type];
118     dr2 = r2 - r01[type];
119     tk1 = kss[type] * dr1;
120     tk2 = kss[type] * dr2;
121 
122     f1[0] = -delx1*tk2/r1;
123     f1[1] = -dely1*tk2/r1;
124     f1[2] = -delz1*tk2/r1;
125 
126     f3[0] = -delx2*tk1/r2;
127     f3[1] = -dely2*tk1/r2;
128     f3[2] = -delz2*tk1/r2;
129 
130     if (eflag) eangle = kss[type]*dr1*dr2;
131 
132     // force & energy for bond-angle term
133     dtheta = acos(c) - theta0[type];
134 
135     aa1 = s * dr1 * kbs0[type];
136     aa2 = s * dr2 * kbs1[type];
137 
138     aa11 = aa1 * c / rsq1;
139     aa12 = -aa1 / (r1 * r2);
140     aa21 = aa2 * c / rsq1;
141     aa22 = -aa2 / (r1 * r2);
142 
143     vx11 = (aa11 * delx1) + (aa12 * delx2);
144     vx12 = (aa21 * delx1) + (aa22 * delx2);
145     vy11 = (aa11 * dely1) + (aa12 * dely2);
146     vy12 = (aa21 * dely1) + (aa22 * dely2);
147     vz11 = (aa11 * delz1) + (aa12 * delz2);
148     vz12 = (aa21 * delz1) + (aa22 * delz2);
149 
150     aa11 = aa1 * c / rsq2;
151     aa21 = aa2 * c / rsq2;
152 
153     vx21 = (aa11 * delx2) + (aa12 * delx1);
154     vx22 = (aa21 * delx2) + (aa22 * delx1);
155     vy21 = (aa11 * dely2) + (aa12 * dely1);
156     vy22 = (aa21 * dely2) + (aa22 * dely1);
157     vz21 = (aa11 * delz2) + (aa12 * delz1);
158     vz22 = (aa21 * delz2) + (aa22 * delz1);
159 
160     b1 = kbs0[type] * dtheta / r1;
161     b2 = kbs1[type] * dtheta / r2;
162 
163     f1[0] -= vx11 + b1*delx1 + vx12;
164     f1[1] -= vy11 + b1*dely1 + vy12;
165     f1[2] -= vz11 + b1*delz1 + vz12;
166 
167     f3[0] -= vx21 + b2*delx2 + vx22;
168     f3[1] -= vy21 + b2*dely2 + vy22;
169     f3[2] -= vz21 + b2*delz2 + vz22;
170 
171     if (eflag) eangle += kbs0[type]*dr1*dtheta + kbs1[type]*dr2*dtheta;
172 
173     // apply force to each of 3 atoms
174 
175     if (newton_bond || i1 < nlocal) {
176       f[i1][0] += f1[0];
177       f[i1][1] += f1[1];
178       f[i1][2] += f1[2];
179     }
180 
181     if (newton_bond || i2 < nlocal) {
182       f[i2][0] -= f1[0] + f3[0];
183       f[i2][1] -= f1[1] + f3[1];
184       f[i2][2] -= f1[2] + f3[2];
185     }
186 
187     if (newton_bond || i3 < nlocal) {
188       f[i3][0] += f3[0];
189       f[i3][1] += f3[1];
190       f[i3][2] += f3[2];
191     }
192 
193     if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
194                          delx1,dely1,delz1,delx2,dely2,delz2);
195   }
196 }
197 
198 /* ---------------------------------------------------------------------- */
199 
allocate()200 void AngleCross::allocate()
201 {
202   allocated = 1;
203   int n = atom->nangletypes;
204 
205   memory->create(kss,n+1,"angle:kss");
206   memory->create(kbs0,n+1,"angle:kbs0");
207   memory->create(kbs1,n+1,"angle:kbs1");
208   memory->create(r00,n+1,"angle:r00");
209   memory->create(r01,n+1,"angle:r01");
210   memory->create(theta0,n+1,"angle:theta0");
211   memory->create(setflag,n+1,"angle:setflag");
212 
213   for (int i = 1; i <= n; i++)
214     setflag[i] = 0;
215 }
216 
217 /* ----------------------------------------------------------------------
218    set coeffs
219 ------------------------------------------------------------------------- */
220 
coeff(int narg,char ** arg)221 void AngleCross::coeff(int narg, char **arg)
222 {
223   if (narg != 7) error->all(FLERR,"Incorrect args for angle coefficients");
224   if (!allocated) allocate();
225 
226   int ilo,ihi;
227   utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
228 
229   int count = 0;
230 
231     double kss_one = utils::numeric(FLERR,arg[1],false,lmp);
232     double kbs0_one = utils::numeric(FLERR,arg[2],false,lmp);
233     double kbs1_one = utils::numeric(FLERR,arg[3],false,lmp);
234     double r0_one = utils::numeric(FLERR,arg[4],false,lmp);
235     double r1_one = utils::numeric(FLERR,arg[5],false,lmp);
236     double theta0_one = utils::numeric(FLERR,arg[6],false,lmp);
237 
238     for (int i = ilo; i <= ihi; i++) {
239       kss[i] = kss_one;
240       kbs0[i] = kbs0_one;
241       kbs1[i] = kbs1_one;
242       r00[i] = r0_one;
243       r01[i] = r1_one;
244       // Convert theta0 from degrees to radians
245       theta0[i] = theta0_one*MY_PI/180.0;
246       setflag[i] = 1;
247       count++;
248     }
249 
250   if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
251 }
252 
253 /* ---------------------------------------------------------------------- */
254 
equilibrium_angle(int i)255 double AngleCross::equilibrium_angle(int i)
256 {
257   return theta0[i];
258 }
259 
260 /* ----------------------------------------------------------------------
261    proc 0 writes out coeffs to restart file
262 ------------------------------------------------------------------------- */
263 
write_restart(FILE * fp)264 void AngleCross::write_restart(FILE *fp)
265 {
266   fwrite(&kss[1],sizeof(double),atom->nangletypes,fp);
267   fwrite(&kbs0[1],sizeof(double),atom->nangletypes,fp);
268   fwrite(&kbs1[1],sizeof(double),atom->nangletypes,fp);
269   fwrite(&r00[1],sizeof(double),atom->nangletypes,fp);
270   fwrite(&r01[1],sizeof(double),atom->nangletypes,fp);
271   fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
272 }
273 
274 /* ----------------------------------------------------------------------
275    proc 0 reads coeffs from restart file, bcasts them
276 ------------------------------------------------------------------------- */
277 
read_restart(FILE * fp)278 void AngleCross::read_restart(FILE *fp)
279 {
280   allocate();
281 
282   if (comm->me == 0) {
283     utils::sfread(FLERR,&kss[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
284     utils::sfread(FLERR,&kbs0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
285     utils::sfread(FLERR,&kbs1[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
286     utils::sfread(FLERR,&r00[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
287     utils::sfread(FLERR,&r01[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
288     utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
289   }
290 
291   MPI_Bcast(&kss[1],atom->nangletypes,MPI_DOUBLE,0,world);
292   MPI_Bcast(&kbs0[1],atom->nangletypes,MPI_DOUBLE,0,world);
293   MPI_Bcast(&kbs1[1],atom->nangletypes,MPI_DOUBLE,0,world);
294   MPI_Bcast(&r00[1],atom->nangletypes,MPI_DOUBLE,0,world);
295   MPI_Bcast(&r01[1],atom->nangletypes,MPI_DOUBLE,0,world);
296   MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
297 
298   for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
299 }
300 
301 /* ----------------------------------------------------------------------
302    proc 0 writes to data file
303 ------------------------------------------------------------------------- */
304 
write_data(FILE * fp)305 void AngleCross::write_data(FILE *fp)
306 {
307   for (int i = 1; i <= atom->nangletypes; i++)
308     fprintf(fp,"%d %g %g %g %g %g %g\n",
309             i,kss[i],kbs0[i],kbs1[i],r00[i],r01[i],theta0[i]/MY_PI*180.0);
310 }
311 
312 /* ---------------------------------------------------------------------- */
313 
single(int type,int i1,int i2,int i3)314 double AngleCross::single(int type, int i1, int i2, int i3)
315 {
316   double **x = atom->x;
317 
318   double delx1 = x[i1][0] - x[i2][0];
319   double dely1 = x[i1][1] - x[i2][1];
320   double delz1 = x[i1][2] - x[i2][2];
321   domain->minimum_image(delx1,dely1,delz1);
322   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
323 
324   double delx2 = x[i3][0] - x[i2][0];
325   double dely2 = x[i3][1] - x[i2][1];
326   double delz2 = x[i3][2] - x[i2][2];
327   domain->minimum_image(delx2,dely2,delz2);
328   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
329 
330   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
331   c /= r1*r2;
332   if (c > 1.0) c = 1.0;
333   if (c < -1.0) c = -1.0;
334 
335   double s = sqrt(1.0 - c*c);
336   if (s < SMALL) s = SMALL;
337   s = 1.0/s;
338 
339   double dtheta = acos(c) - theta0[type];
340   double dr1 = r1 - r00[type];
341   double dr2 = r2 - r01[type];
342   double energy = kss[type]*dr1*dr2+kbs0[type]*dr1*dtheta + kbs1[type]*dr2*dtheta;
343   return energy;
344 }
345