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: Loukas D. Peristeras (Scienomics SARL)
17    [ based on angle_cosine_squared.cpp Naveen Michaud-Agrawal (Johns Hopkins U)]
18 ------------------------------------------------------------------------- */
19 
20 #include "angle_fourier.h"
21 
22 #include <cmath>
23 #include "atom.h"
24 #include "neighbor.h"
25 #include "domain.h"
26 #include "comm.h"
27 #include "force.h"
28 #include "math_const.h"
29 #include "memory.h"
30 #include "error.h"
31 
32 
33 using namespace LAMMPS_NS;
34 using namespace MathConst;
35 
36 #define SMALL 0.001
37 
38 /* ---------------------------------------------------------------------- */
39 
AngleFourier(LAMMPS * lmp)40 AngleFourier::AngleFourier(LAMMPS *lmp) : Angle(lmp)
41 {
42   k = nullptr;
43   C0 = nullptr;
44   C1 = nullptr;
45   C2 = nullptr;
46 }
47 
48 /* ---------------------------------------------------------------------- */
49 
~AngleFourier()50 AngleFourier::~AngleFourier()
51 {
52   if (allocated) {
53     memory->destroy(setflag);
54     memory->destroy(k);
55     memory->destroy(C0);
56     memory->destroy(C1);
57     memory->destroy(C2);
58   }
59 }
60 
61 /* ---------------------------------------------------------------------- */
62 
compute(int eflag,int vflag)63 void AngleFourier::compute(int eflag, int vflag)
64 {
65   int i1,i2,i3,n,type;
66   double delx1,dely1,delz1,delx2,dely2,delz2;
67   double eangle,f1[3],f3[3];
68   double term;
69   double rsq1,rsq2,r1,r2,c,c2,a,a11,a12,a22;
70 
71   eangle = 0.0;
72   ev_init(eflag,vflag);
73 
74   double **x = atom->x;
75   double **f = atom->f;
76   int **anglelist = neighbor->anglelist;
77   int nanglelist = neighbor->nanglelist;
78   int nlocal = atom->nlocal;
79   int newton_bond = force->newton_bond;
80 
81   for (n = 0; n < nanglelist; n++) {
82     i1 = anglelist[n][0];
83     i2 = anglelist[n][1];
84     i3 = anglelist[n][2];
85     type = anglelist[n][3];
86 
87     // 1st bond
88 
89     delx1 = x[i1][0] - x[i2][0];
90     dely1 = x[i1][1] - x[i2][1];
91     delz1 = x[i1][2] - x[i2][2];
92 
93     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
94     r1 = sqrt(rsq1);
95 
96     // 2nd bond
97 
98     delx2 = x[i3][0] - x[i2][0];
99     dely2 = x[i3][1] - x[i2][1];
100     delz2 = x[i3][2] - x[i2][2];
101 
102     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
103     r2 = sqrt(rsq2);
104 
105     // angle (cos and sin)
106 
107     c = delx1*delx2 + dely1*dely2 + delz1*delz2;
108     c /= r1*r2;
109 
110     if (c > 1.0) c = 1.0;
111     if (c < -1.0) c = -1.0;
112 
113     // force & energy
114 
115     c2 = 2.0*c*c-1.0;
116     term = k[type]*(C0[type]+C1[type]*c+C2[type]*c2);
117 
118     if (eflag) eangle = term;
119 
120     a = k[type]*(C1[type]+4.0*C2[type]*c);
121     a11 = a*c / rsq1;
122     a12 = -a / (r1*r2);
123     a22 = a*c / rsq2;
124 
125     f1[0] = a11*delx1 + a12*delx2;
126     f1[1] = a11*dely1 + a12*dely2;
127     f1[2] = a11*delz1 + a12*delz2;
128     f3[0] = a22*delx2 + a12*delx1;
129     f3[1] = a22*dely2 + a12*dely1;
130     f3[2] = a22*delz2 + a12*delz1;
131 
132     // apply force to each of 3 atoms
133 
134     if (newton_bond || i1 < nlocal) {
135       f[i1][0] += f1[0];
136       f[i1][1] += f1[1];
137       f[i1][2] += f1[2];
138     }
139 
140     if (newton_bond || i2 < nlocal) {
141       f[i2][0] -= f1[0] + f3[0];
142       f[i2][1] -= f1[1] + f3[1];
143       f[i2][2] -= f1[2] + f3[2];
144     }
145 
146     if (newton_bond || i3 < nlocal) {
147       f[i3][0] += f3[0];
148       f[i3][1] += f3[1];
149       f[i3][2] += f3[2];
150     }
151 
152     if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
153                          delx1,dely1,delz1,delx2,dely2,delz2);
154   }
155 }
156 
157 /* ---------------------------------------------------------------------- */
158 
allocate()159 void AngleFourier::allocate()
160 {
161   allocated = 1;
162   int n = atom->nangletypes;
163 
164   memory->create(k,n+1,"angle:k");
165   memory->create(C0,n+1,"angle:C0");
166   memory->create(C1,n+1,"angle:C1");
167   memory->create(C2,n+1,"angle:C2");
168 
169   memory->create(setflag,n+1,"angle:setflag");
170   for (int i = 1; i <= n; i++) setflag[i] = 0;
171 }
172 
173 /* ----------------------------------------------------------------------
174    set coeffs for one or more types
175 ------------------------------------------------------------------------- */
176 
coeff(int narg,char ** arg)177 void AngleFourier::coeff(int narg, char **arg)
178 {
179   if (narg != 5) error->all(FLERR,"Incorrect args for angle coefficients");
180   if (!allocated) allocate();
181 
182   int ilo,ihi;
183   utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
184 
185   double k_one =  utils::numeric(FLERR,arg[1],false,lmp);
186   double C0_one = utils::numeric(FLERR,arg[2],false,lmp);
187   double C1_one = utils::numeric(FLERR,arg[3],false,lmp);
188   double C2_one = utils::numeric(FLERR,arg[4],false,lmp);
189 
190   int count = 0;
191   for (int i = ilo; i <= ihi; i++) {
192     k[i] = k_one;
193     C0[i] = C0_one;
194     C1[i] = C1_one;
195     C2[i] = C2_one;
196     setflag[i] = 1;
197     count++;
198   }
199 
200   if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
201 }
202 
203 /* ---------------------------------------------------------------------- */
204 
equilibrium_angle(int i)205 double AngleFourier::equilibrium_angle(int i)
206 {
207   double ret=MY_PI;
208   if (C2[i] != 0.0) {
209     ret = (C1[i]/4.0/C2[i]);
210     if (fabs(ret) <= 1.0) ret = acos(-ret);
211   }
212   return ret;
213 }
214 
215 /* ----------------------------------------------------------------------
216    proc 0 writes out coeffs to restart file
217 ------------------------------------------------------------------------- */
218 
write_restart(FILE * fp)219 void AngleFourier::write_restart(FILE *fp)
220 {
221   fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
222   fwrite(&C0[1],sizeof(double),atom->nangletypes,fp);
223   fwrite(&C1[1],sizeof(double),atom->nangletypes,fp);
224   fwrite(&C2[1],sizeof(double),atom->nangletypes,fp);
225 }
226 
227 /* ----------------------------------------------------------------------
228    proc 0 reads coeffs from restart file, bcasts them
229 ------------------------------------------------------------------------- */
230 
read_restart(FILE * fp)231 void AngleFourier::read_restart(FILE *fp)
232 {
233   allocate();
234 
235   if (comm->me == 0) {
236     utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
237     utils::sfread(FLERR,&C0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
238     utils::sfread(FLERR,&C1[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
239     utils::sfread(FLERR,&C2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
240   }
241   MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
242   MPI_Bcast(&C0[1],atom->nangletypes,MPI_DOUBLE,0,world);
243   MPI_Bcast(&C1[1],atom->nangletypes,MPI_DOUBLE,0,world);
244   MPI_Bcast(&C2[1],atom->nangletypes,MPI_DOUBLE,0,world);
245 
246   for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
247 }
248 
249 /* ----------------------------------------------------------------------
250    proc 0 writes to data file
251 ------------------------------------------------------------------------- */
252 
write_data(FILE * fp)253 void AngleFourier::write_data(FILE *fp)
254 {
255   for (int i = 1; i <= atom->nangletypes; i++)
256     fprintf(fp,"%d %g %g %g %g\n",i,k[i],C0[i],C1[i],C2[i]);
257 }
258 
259 /* ---------------------------------------------------------------------- */
260 
single(int type,int i1,int i2,int i3)261 double AngleFourier::single(int type, int i1, int i2, int i3)
262 {
263   double **x = atom->x;
264 
265   double delx1 = x[i1][0] - x[i2][0];
266   double dely1 = x[i1][1] - x[i2][1];
267   double delz1 = x[i1][2] - x[i2][2];
268   domain->minimum_image(delx1,dely1,delz1);
269   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
270 
271   double delx2 = x[i3][0] - x[i2][0];
272   double dely2 = x[i3][1] - x[i2][1];
273   double delz2 = x[i3][2] - x[i2][2];
274   domain->minimum_image(delx2,dely2,delz2);
275   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
276 
277   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
278   c /= r1*r2;
279   if (c > 1.0) c = 1.0;
280   if (c < -1.0) c = -1.0;
281   double c2 = 2.0*c*c-1.0;
282 
283   double eng = k[type]*(C0[type]+C1[type]*c+C2[type]*c2);
284   return eng;
285 }
286