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 improper_umbrella.cpp Tod A Pascal (Caltech) ]
18 ------------------------------------------------------------------------- */
19 
20 #include "improper_fourier.h"
21 
22 #include "atom.h"
23 #include "comm.h"
24 #include "error.h"
25 #include "force.h"
26 #include "memory.h"
27 #include "neighbor.h"
28 
29 #include <cmath>
30 
31 using namespace LAMMPS_NS;
32 
33 #define TOLERANCE 0.05
34 #define SMALL     0.001
35 
36 /* ---------------------------------------------------------------------- */
37 
ImproperFourier(LAMMPS * lmp)38 ImproperFourier::ImproperFourier(LAMMPS *lmp) : Improper(lmp) {}
39 
40 /* ---------------------------------------------------------------------- */
41 
~ImproperFourier()42 ImproperFourier::~ImproperFourier()
43 {
44   if (allocated) {
45     memory->destroy(setflag);
46     memory->destroy(k);
47     memory->destroy(C0);
48     memory->destroy(C1);
49     memory->destroy(C2);
50     memory->destroy(all);
51   }
52 }
53 
54 /* ---------------------------------------------------------------------- */
55 
compute(int eflag,int vflag)56 void ImproperFourier::compute(int eflag, int vflag)
57 {
58   int i1,i2,i3,i4,n,type;
59   double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z;
60 
61   ev_init(eflag,vflag);
62 
63   double **x = atom->x;
64   int **improperlist = neighbor->improperlist;
65   int nimproperlist = neighbor->nimproperlist;
66 
67   for (n = 0; n < nimproperlist; n++) {
68     i1 = improperlist[n][0];
69     i2 = improperlist[n][1];
70     i3 = improperlist[n][2];
71     i4 = improperlist[n][3];
72     type = improperlist[n][4];
73 
74     // 1st bond
75 
76     vb1x = x[i2][0] - x[i1][0];
77     vb1y = x[i2][1] - x[i1][1];
78     vb1z = x[i2][2] - x[i1][2];
79 
80     // 2nd bond
81 
82     vb2x = x[i3][0] - x[i1][0];
83     vb2y = x[i3][1] - x[i1][1];
84     vb2z = x[i3][2] - x[i1][2];
85 
86     // 3rd bond
87 
88     vb3x = x[i4][0] - x[i1][0];
89     vb3y = x[i4][1] - x[i1][1];
90     vb3z = x[i4][2] - x[i1][2];
91 
92     addone(i1,i2,i3,i4, type,evflag,eflag,
93            vb1x, vb1y, vb1z,
94            vb2x, vb2y, vb2z,
95            vb3x, vb3y, vb3z);
96     if (all[type]) {
97       addone(i1,i4,i2,i3, type,evflag,eflag,
98              vb3x, vb3y, vb3z,
99              vb1x, vb1y, vb1z,
100              vb2x, vb2y, vb2z);
101       addone(i1,i3,i4,i2, type,evflag,eflag,
102              vb2x, vb2y, vb2z,
103              vb3x, vb3y, vb3z,
104              vb1x, vb1y, vb1z);
105     }
106   }
107 }
108 
addone(const int & i1,const int & i2,const int & i3,const int & i4,const int & type,const int & evflag,const int & eflag,const double & vb1x,const double & vb1y,const double & vb1z,const double & vb2x,const double & vb2y,const double & vb2z,const double & vb3x,const double & vb3y,const double & vb3z)109 void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int &i4,
110             const int &type,const int &evflag,const int &eflag,
111             const double &vb1x, const double &vb1y, const double &vb1z,
112             const double &vb2x, const double &vb2y, const double &vb2z,
113             const double &vb3x, const double &vb3y, const double &vb3z)
114 {
115   double eimproper,f1[3],f2[3],f3[3],f4[3];
116   double c,c2,a,s,projhfg,dhax,dhay,dhaz,dahx,dahy,dahz,cotphi;
117   double ax,ay,az,ra2,rh2,ra,rh,rar,rhr,arx,ary,arz,hrx,hry,hrz;
118 
119   double **f = atom->f;
120   int nlocal = atom->nlocal;
121   int newton_bond = force->newton_bond;
122 
123   eimproper = 0.0;
124 
125   // c0 calculation
126   // A = vb1 X vb2 is perpendicular to IJK plane
127 
128   ax = vb1y*vb2z-vb1z*vb2y;
129   ay = vb1z*vb2x-vb1x*vb2z;
130   az = vb1x*vb2y-vb1y*vb2x;
131   ra2 = ax*ax+ay*ay+az*az;
132   rh2 = vb3x*vb3x+vb3y*vb3y+vb3z*vb3z;
133   ra = sqrt(ra2);
134   rh = sqrt(rh2);
135   if (ra < SMALL) ra = SMALL;
136   if (rh < SMALL) rh = SMALL;
137 
138   rar = 1/ra;
139   rhr = 1/rh;
140   arx = ax*rar;
141   ary = ay*rar;
142   arz = az*rar;
143   hrx = vb3x*rhr;
144   hry = vb3y*rhr;
145   hrz = vb3z*rhr;
146 
147   c = arx*hrx+ary*hry+arz*hrz;
148 
149   // error check
150 
151   if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
152     problem(FLERR, i1, i2, i3, i4);
153 
154   if (c > 1.0) c = 1.0;
155   if (c < -1.0) c = -1.0;
156 
157   s = sqrt(1.0 - c*c);
158   if (s < SMALL) s = SMALL;
159   cotphi = c/s;
160 
161   projhfg = (vb3x*vb1x+vb3y*vb1y+vb3z*vb1z) /
162     sqrt(vb1x*vb1x+vb1y*vb1y+vb1z*vb1z);
163   projhfg += (vb3x*vb2x+vb3y*vb2y+vb3z*vb2z) /
164     sqrt(vb2x*vb2x+vb2y*vb2y+vb2z*vb2z);
165   if (projhfg > 0.0) {
166     s *= -1.0;
167     cotphi *= -1.0;
168   }
169 
170   //  force and energy
171   //  E = k ( C0 + C1 cos(w) + C2 cos(w)
172 
173   c2 = 2.0*s*s-1.0;
174   if (eflag) eimproper =  k[type]*(C0[type]+C1[type]*s+C2[type]*c2);
175 
176   // dhax = diffrence between H and A in X direction, etc
177 
178   a = k[type]*(C1[type]+4.0*C2[type]*s)*cotphi;
179   dhax = hrx-c*arx;
180   dhay = hry-c*ary;
181   dhaz = hrz-c*arz;
182 
183   dahx = arx-c*hrx;
184   dahy = ary-c*hry;
185   dahz = arz-c*hrz;
186 
187   f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
188   f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
189   f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
190 
191   f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
192   f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
193   f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
194 
195   f4[0] = dahx*rhr*a;
196   f4[1] = dahy*rhr*a;
197   f4[2] = dahz*rhr*a;
198 
199   f1[0] = -(f2[0] + f3[0] + f4[0]);
200   f1[1] = -(f2[1] + f3[1] + f4[1]);
201   f1[2] = -(f2[2] + f3[2] + f4[2]);
202 
203   // apply force to each of 4 atoms
204 
205   if (newton_bond || i1 < nlocal) {
206     f[i1][0] += f1[0];
207     f[i1][1] += f1[1];
208     f[i1][2] += f1[2];
209   }
210 
211   if (newton_bond || i2 < nlocal) {
212     f[i2][0] += f3[0];
213     f[i2][1] += f3[1];
214     f[i2][2] += f3[2];
215   }
216 
217   if (newton_bond || i3 < nlocal) {
218     f[i3][0] += f2[0];
219     f[i3][1] += f2[1];
220     f[i3][2] += f2[2];
221   }
222 
223   if (newton_bond || i4 < nlocal) {
224     f[i4][0] += f4[0];
225     f[i4][1] += f4[1];
226     f[i4][2] += f4[2];
227   }
228 
229   if (evflag)
230       ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,
231                -vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z);
232 }
233 
234 /* ---------------------------------------------------------------------- */
235 
allocate()236 void ImproperFourier::allocate()
237 {
238   allocated = 1;
239   int n = atom->nimpropertypes;
240 
241   memory->create(k,n+1,"improper:k");
242   memory->create(C0,n+1,"improper:C0");
243   memory->create(C1,n+1,"improper:C1");
244   memory->create(C2,n+1,"improper:C2");
245   memory->create(all,n+1,"improper:C2");
246 
247   memory->create(setflag,n+1,"improper:setflag");
248   for (int i = 1; i <= n; i++) setflag[i] = 0;
249 }
250 
251 /* ----------------------------------------------------------------------
252    set coeffs for one type
253 ------------------------------------------------------------------------- */
254 
coeff(int narg,char ** arg)255 void ImproperFourier::coeff(int narg, char **arg)
256 {
257 
258   if (narg != 5 && narg != 6) error->all(FLERR,"Incorrect args for improper coefficients");
259 
260   if (!allocated) allocate();
261 
262   int ilo,ihi;
263   utils::bounds(FLERR,arg[0],1,atom->nimpropertypes,ilo,ihi,error);
264 
265   double k_one = utils::numeric(FLERR,arg[1],false,lmp);
266   double C0_one = utils::numeric(FLERR,arg[2],false,lmp);
267   double C1_one = utils::numeric(FLERR,arg[3],false,lmp);
268   double C2_one = utils::numeric(FLERR,arg[4],false,lmp);
269   int all_one = 1;
270   if (narg == 6) all_one = utils::inumeric(FLERR,arg[5],false,lmp);
271 
272   // convert w0 from degrees to radians
273 
274   int count = 0;
275   for (int i = ilo; i <= ihi; i++) {
276     k[i] = k_one;
277     C0[i] = C0_one;
278     C1[i] = C1_one;
279     C2[i] = C2_one;
280     all[i] = all_one;
281     setflag[i] = 1;
282     count++;
283   }
284 
285   if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients");
286 }
287 
288 /* ----------------------------------------------------------------------
289    proc 0 writes out coeffs to restart file
290 ------------------------------------------------------------------------- */
291 
write_restart(FILE * fp)292 void ImproperFourier::write_restart(FILE *fp)
293 {
294   fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp);
295   fwrite(&C0[1],sizeof(double),atom->nimpropertypes,fp);
296   fwrite(&C1[1],sizeof(double),atom->nimpropertypes,fp);
297   fwrite(&C2[1],sizeof(double),atom->nimpropertypes,fp);
298   fwrite(&all[1],sizeof(int),atom->nimpropertypes,fp);
299 }
300 
301 /* ----------------------------------------------------------------------
302    proc 0 reads coeffs from restart file, bcasts them
303 ------------------------------------------------------------------------- */
304 
read_restart(FILE * fp)305 void ImproperFourier::read_restart(FILE *fp)
306 {
307   allocate();
308 
309   if (comm->me == 0) {
310     utils::sfread(FLERR,&k[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
311     utils::sfread(FLERR,&C0[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
312     utils::sfread(FLERR,&C1[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
313     utils::sfread(FLERR,&C2[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
314     utils::sfread(FLERR,&all[1],sizeof(int),atom->nimpropertypes,fp,nullptr,error);
315   }
316   MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
317   MPI_Bcast(&C0[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
318   MPI_Bcast(&C1[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
319   MPI_Bcast(&C2[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
320   MPI_Bcast(&all[1],atom->nimpropertypes,MPI_INT,0,world);
321 
322   for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;
323 }
324 
325 /* ----------------------------------------------------------------------
326    proc 0 writes to data file
327 ------------------------------------------------------------------------- */
328 
write_data(FILE * fp)329 void ImproperFourier::write_data(FILE *fp)
330 {
331   for (int i = 1; i <= atom->nimpropertypes; i++)
332     fprintf(fp,"%d %g %g %g %g %d\n",i,k[i],C0[i],C1[i],C2[i],all[i]);
333 }
334