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