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