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_mm3.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
AngleMM3(LAMMPS * lmp)39 AngleMM3::AngleMM3(LAMMPS *lmp) : Angle(lmp) {}
40
41 /* ---------------------------------------------------------------------- */
42
~AngleMM3()43 AngleMM3::~AngleMM3()
44 {
45 if (copymode) return;
46
47 if (allocated) {
48 memory->destroy(setflag);
49 memory->destroy(theta0);
50 memory->destroy(k2);
51 }
52 }
53
54 /* ---------------------------------------------------------------------- */
55
compute(int eflag,int vflag)56 void AngleMM3::compute(int eflag, int vflag)
57 {
58 int i1,i2,i3,n,type;
59 double delx1,dely1,delz1,delx2,dely2,delz2;
60 double eangle,f1[3],f3[3];
61 double dtheta,dtheta2,dtheta3,dtheta4,de_angle;
62 double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
63
64 eangle = 0.0;
65 ev_init(eflag,vflag);
66
67 double **x = atom->x;
68 double **f = atom->f;
69 int **anglelist = neighbor->anglelist;
70 int nanglelist = neighbor->nanglelist;
71 int nlocal = atom->nlocal;
72 int newton_bond = force->newton_bond;
73
74 for (n = 0; n < nanglelist; n++) {
75 i1 = anglelist[n][0];
76 i2 = anglelist[n][1];
77 i3 = anglelist[n][2];
78 type = anglelist[n][3];
79
80 // 1st bond
81
82 delx1 = x[i1][0] - x[i2][0];
83 dely1 = x[i1][1] - x[i2][1];
84 delz1 = x[i1][2] - x[i2][2];
85
86 rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
87 r1 = sqrt(rsq1);
88
89 // 2nd bond
90
91 delx2 = x[i3][0] - x[i2][0];
92 dely2 = x[i3][1] - x[i2][1];
93 delz2 = x[i3][2] - x[i2][2];
94
95 rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
96 r2 = sqrt(rsq2);
97
98 // angle (cos and sin)
99
100 c = delx1*delx2 + dely1*dely2 + delz1*delz2;
101 c /= r1*r2;
102
103 if (c > 1.0) c = 1.0;
104 if (c < -1.0) c = -1.0;
105
106 s = sqrt(1.0 - c*c);
107 if (s < SMALL) s = SMALL;
108 s = 1.0/s;
109
110 // force & energy for angle term
111
112 dtheta = acos(c) - theta0[type];
113 dtheta2 = dtheta*dtheta;
114 dtheta3 = dtheta2*dtheta;
115 dtheta4 = dtheta3*dtheta;
116 // MM3 angle term, taking into account that dtheta is expressed in rad
117 de_angle = 2.0*k2[type]*dtheta*(1.0-1.203211*dtheta+0.367674*dtheta2-0.3239159*dtheta3+0.711270*dtheta4);
118
119 a = -de_angle*s;
120 a11 = a*c / rsq1;
121 a12 = -a / (r1*r2);
122 a22 = a*c / rsq2;
123
124 f1[0] = a11*delx1 + a12*delx2;
125 f1[1] = a11*dely1 + a12*dely2;
126 f1[2] = a11*delz1 + a12*delz2;
127
128 f3[0] = a22*delx2 + a12*delx1;
129 f3[1] = a22*dely2 + a12*dely1;
130 f3[2] = a22*delz2 + a12*delz1;
131 // MM3 angle term, taking into account that dtheta is expressed in rad
132 if (eflag) eangle = k2[type]*dtheta2*(1.0-0.802141*dtheta+0.183837*dtheta2-0.131664*dtheta3+0.237090*dtheta4);
133
134 // apply force to each of 3 atoms
135
136 if (newton_bond || i1 < nlocal) {
137 f[i1][0] += f1[0];
138 f[i1][1] += f1[1];
139 f[i1][2] += f1[2];
140 }
141
142 if (newton_bond || i2 < nlocal) {
143 f[i2][0] -= f1[0] + f3[0];
144 f[i2][1] -= f1[1] + f3[1];
145 f[i2][2] -= f1[2] + f3[2];
146 }
147
148 if (newton_bond || i3 < nlocal) {
149 f[i3][0] += f3[0];
150 f[i3][1] += f3[1];
151 f[i3][2] += f3[2];
152 }
153
154 if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
155 delx1,dely1,delz1,delx2,dely2,delz2);
156 }
157 }
158
159 /* ---------------------------------------------------------------------- */
160
allocate()161 void AngleMM3::allocate()
162 {
163 allocated = 1;
164 int n = atom->nangletypes;
165
166 memory->create(setflag,n+1,"angle:setflag");
167 memory->create(k2,n+1,"angle:k2");
168 memory->create(theta0,n+1,"angle:theta0");
169 for (int i = 1; i <= n; i++)
170 setflag[i] = 0;
171 }
172
173 /* ----------------------------------------------------------------------
174 set coeffs
175 else -> Angle coeffs
176 ------------------------------------------------------------------------- */
177
coeff(int narg,char ** arg)178 void AngleMM3::coeff(int narg, char **arg)
179 {
180 if (narg != 3) error->all(FLERR,"Incorrect args for angle coefficients");
181 if (!allocated) allocate();
182
183 int ilo,ihi;
184 utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
185
186 int count = 0;
187
188 double k2_one = utils::numeric(FLERR,arg[1],false,lmp);
189 double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
190
191 // convert theta0 from degrees to radians
192
193 for (int i = ilo; i <= ihi; i++) {
194 k2[i] = k2_one;
195 theta0[i] = theta0_one/180.0 * MY_PI;
196 setflag[i] = 1;
197 count++;
198 }
199
200 if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
201
202 }
203
204 /* ---------------------------------------------------------------------- */
205
equilibrium_angle(int i)206 double AngleMM3::equilibrium_angle(int i)
207 {
208 return theta0[i];
209 }
210
211 /* ----------------------------------------------------------------------
212 proc 0 writes out coeffs to restart file
213 ------------------------------------------------------------------------- */
214
write_restart(FILE * fp)215 void AngleMM3::write_restart(FILE *fp)
216 {
217 fwrite(&k2[1],sizeof(double),atom->nangletypes,fp);
218 fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
219 }
220
221 /* ----------------------------------------------------------------------
222 proc 0 reads coeffs from restart file, bcasts them
223 ------------------------------------------------------------------------- */
224
read_restart(FILE * fp)225 void AngleMM3::read_restart(FILE *fp)
226 {
227 allocate();
228
229 if (comm->me == 0) {
230 utils::sfread(FLERR,&k2[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
231 utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
232 }
233
234 MPI_Bcast(&k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
235 MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
236
237 for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
238 }
239
240 /* ----------------------------------------------------------------------
241 proc 0 writes to data file
242 ------------------------------------------------------------------------- */
243
write_data(FILE * fp)244 void AngleMM3::write_data(FILE *fp)
245 {
246 for (int i = 1; i <= atom->nangletypes; i++)
247 fprintf(fp,"%d %g %g\n",
248 i,k2[i],theta0[i]/MY_PI*180.0);
249 }
250
251 /* ---------------------------------------------------------------------- */
252
single(int type,int i1,int i2,int i3)253 double AngleMM3::single(int type, int i1, int i2, int i3)
254 {
255 double **x = atom->x;
256
257 double delx1 = x[i1][0] - x[i2][0];
258 double dely1 = x[i1][1] - x[i2][1];
259 double delz1 = x[i1][2] - x[i2][2];
260 domain->minimum_image(delx1,dely1,delz1);
261 double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
262
263 double delx2 = x[i3][0] - x[i2][0];
264 double dely2 = x[i3][1] - x[i2][1];
265 double delz2 = x[i3][2] - x[i2][2];
266 domain->minimum_image(delx2,dely2,delz2);
267 double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
268
269 double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
270 c /= r1*r2;
271 if (c > 1.0) c = 1.0;
272 if (c < -1.0) c = -1.0;
273
274 double s = sqrt(1.0 - c*c);
275 if (s < SMALL) s = SMALL;
276 s = 1.0/s;
277
278 double dtheta = acos(c) - theta0[type];
279 double dtheta2 = dtheta*dtheta;
280 double dtheta3 = dtheta2*dtheta;
281 double dtheta4 = dtheta3*dtheta;
282
283 double energy = k2[type]*dtheta2*(1.0-0.802141*dtheta+0.183837*dtheta2-0.131664*dtheta3+0.237090*dtheta4);
284
285 return energy;
286 }
287