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: Axel Kohlmeyer (Temple U), akohlmey at gmail.com
17 ------------------------------------------------------------------------- */
18
19 #include "angle_cosine_delta.h"
20 #include <cmath>
21 #include "atom.h"
22 #include "neighbor.h"
23 #include "domain.h"
24 #include "force.h"
25
26 using namespace LAMMPS_NS;
27
28 #define SMALL 0.001
29
30 /* ---------------------------------------------------------------------- */
31
AngleCosineDelta(LAMMPS * lmp)32 AngleCosineDelta::AngleCosineDelta(LAMMPS *lmp) : AngleCosineSquared(lmp) {}
33
34 /* ---------------------------------------------------------------------- */
35
compute(int eflag,int vflag)36 void AngleCosineDelta::compute(int eflag, int vflag)
37 {
38 int i1,i2,i3,n,type;
39 double delx1,dely1,delz1,delx2,dely2,delz2,theta,dtheta,dcostheta,tk;
40 double eangle,f1[3],f3[3];
41 double rsq1,rsq2,r1,r2,c,a,cot,a11,a12,a22,b11,b12,b22,c0,s0,s;
42
43 eangle = 0.0;
44 ev_init(eflag,vflag);
45
46 double **x = atom->x;
47 double **f = atom->f;
48 int **anglelist = neighbor->anglelist;
49 int nanglelist = neighbor->nanglelist;
50 int nlocal = atom->nlocal;
51 int newton_bond = force->newton_bond;
52
53 for (n = 0; n < nanglelist; n++) {
54 i1 = anglelist[n][0];
55 i2 = anglelist[n][1];
56 i3 = anglelist[n][2];
57 type = anglelist[n][3];
58
59 // 1st bond
60
61 delx1 = x[i1][0] - x[i2][0];
62 dely1 = x[i1][1] - x[i2][1];
63 delz1 = x[i1][2] - x[i2][2];
64
65 rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
66 r1 = sqrt(rsq1);
67
68 // 2nd bond
69
70 delx2 = x[i3][0] - x[i2][0];
71 dely2 = x[i3][1] - x[i2][1];
72 delz2 = x[i3][2] - x[i2][2];
73
74 rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
75 r2 = sqrt(rsq2);
76
77 // angle (cos and sin)
78
79 c = delx1*delx2 + dely1*dely2 + delz1*delz2;
80 c /= r1*r2;
81
82 if (c > 1.0) c = 1.0;
83 if (c < -1.0) c = -1.0;
84
85 theta = acos(c);
86
87 s = sqrt(1.0 - c*c);
88 if (s < SMALL) s = SMALL;
89 s = 1.0/s;
90
91 cot = c/s;
92
93 // force & energy
94
95 dtheta = theta - theta0[type];
96 dcostheta = cos(dtheta);
97 tk = k[type] * (1.0-dcostheta);
98
99 if (eflag) eangle = tk;
100
101 a = -k[type];
102
103 // expand dtheta for cos and sin contribution to force
104
105 a11 = a*c / rsq1;
106 a12 = -a / (r1*r2);
107 a22 = a*c / rsq2;
108
109 b11 = -a*c*cot / rsq1;
110 b12 = a*cot / (r1*r2);
111 b22 = -a*c*cot / rsq2;
112
113 c0 = cos(theta0[type]);
114 s0 = sin(theta0[type]);
115
116 f1[0] = (a11*delx1 + a12*delx2)*c0 + (b11*delx1 + b12*delx2)*s0;
117 f1[1] = (a11*dely1 + a12*dely2)*c0 + (b11*dely1 + b12*dely2)*s0;
118 f1[2] = (a11*delz1 + a12*delz2)*c0 + (b11*delz1 + b12*delz2)*s0;
119 f3[0] = (a22*delx2 + a12*delx1)*c0 + (b22*delx2 + b12*delx1)*s0;
120 f3[1] = (a22*dely2 + a12*dely1)*c0 + (b22*dely2 + b12*dely1)*s0;
121 f3[2] = (a22*delz2 + a12*delz1)*c0 + (b22*delz2 + b12*delz1)*s0;
122
123 // apply force to each of 3 atoms
124
125 if (newton_bond || i1 < nlocal) {
126 f[i1][0] += f1[0];
127 f[i1][1] += f1[1];
128 f[i1][2] += f1[2];
129 }
130
131 if (newton_bond || i2 < nlocal) {
132 f[i2][0] -= f1[0] + f3[0];
133 f[i2][1] -= f1[1] + f3[1];
134 f[i2][2] -= f1[2] + f3[2];
135 }
136
137 if (newton_bond || i3 < nlocal) {
138 f[i3][0] += f3[0];
139 f[i3][1] += f3[1];
140 f[i3][2] += f3[2];
141 }
142
143 if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
144 delx1,dely1,delz1,delx2,dely2,delz2);
145 }
146 }
147
148 /* ---------------------------------------------------------------------- */
149
single(int type,int i1,int i2,int i3)150 double AngleCosineDelta::single(int type, int i1, int i2, int i3)
151 {
152 double **x = atom->x;
153
154 double delx1 = x[i1][0] - x[i2][0];
155 double dely1 = x[i1][1] - x[i2][1];
156 double delz1 = x[i1][2] - x[i2][2];
157 domain->minimum_image(delx1,dely1,delz1);
158 double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
159
160 double delx2 = x[i3][0] - x[i2][0];
161 double dely2 = x[i3][1] - x[i2][1];
162 double delz2 = x[i3][2] - x[i2][2];
163 domain->minimum_image(delx2,dely2,delz2);
164 double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
165
166 double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
167 c /= r1*r2;
168 if (c > 1.0) c = 1.0;
169 if (c < -1.0) c = -1.0;
170
171 double theta = acos(c);
172 double dtheta = theta - theta0[type];
173 double dcostheta = cos(dtheta);
174 double tk = k[type] * (1.0-dcostheta);
175 return tk;
176 }
177