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