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, heavily based on the
17 improper_distance code by Paolo Raiteri (Curtin University)
18 ------------------------------------------------------------------------- */
19
20 #include "improper_distharm.h"
21
22 #include <cmath>
23 #include "atom.h"
24 #include "comm.h"
25 #include "neighbor.h"
26 #include "domain.h"
27 #include "force.h"
28 #include "memory.h"
29 #include "error.h"
30
31
32 using namespace LAMMPS_NS;
33
34 #define TOLERANCE 0.05
35 #define SMALL 0.001
36
37 /* ---------------------------------------------------------------------- */
38
ImproperDistHarm(LAMMPS * lmp)39 ImproperDistHarm::ImproperDistHarm(LAMMPS *lmp) : Improper(lmp) {}
40
41 /* ---------------------------------------------------------------------- */
42
~ImproperDistHarm()43 ImproperDistHarm::~ImproperDistHarm()
44 {
45 if (allocated) {
46 memory->destroy(setflag);
47 memory->destroy(k);
48 memory->destroy(chi);
49 }
50 }
51
52 /* ---------------------------------------------------------------------- */
53
compute(int eflag,int vflag)54 void ImproperDistHarm::compute(int eflag, int vflag)
55 {
56 int i1,i2,i3,i4,n,type;
57 double xab, yab, zab; // bond 1-2
58 double xac, yac, zac; // bond 1-3
59 double xad, yad, zad; // bond 1-4
60 double xbc, ybc, zbc; // bond 2-3
61 double xbd, ybd, zbd; // bond 2-4
62 double xcd, ycd, zcd; // bond 3-4
63 double xna, yna, zna, rna; // normal
64 double da;
65
66 double eimproper,f1[3],f2[3],f3[3],f4[3];
67 double domega,a;
68
69 eimproper = 0.0;
70 ev_init(eflag,vflag);
71
72 double **x = atom->x;
73 double **f = atom->f;
74 int **improperlist = neighbor->improperlist;
75 int nimproperlist = neighbor->nimproperlist;
76 int nlocal = atom->nlocal;
77 int newton_bond = force->newton_bond;
78
79 for (n = 0; n < nimproperlist; n++) {
80 i1 = improperlist[n][0];
81 i2 = improperlist[n][1];
82 i3 = improperlist[n][2];
83 i4 = improperlist[n][3];
84 type = improperlist[n][4];
85
86 // geometry of 4-body
87 // 4 is the central atom
88 // 1-2-3 are ment to be equivalent
89 // I need the bonds between 2-3 and 3-4 to get the plane normal
90 // Then I need the bond 1-4 to project it onto the normal to the plane
91
92 // bond 1->2
93 xab = x[i2][0] - x[i1][0];
94 yab = x[i2][1] - x[i1][1];
95 zab = x[i2][2] - x[i1][2];
96 domain->minimum_image(xab,yab,zab);
97
98 // bond 1->3
99 xac = x[i3][0] - x[i1][0];
100 yac = x[i3][1] - x[i1][1];
101 zac = x[i3][2] - x[i1][2];
102 domain->minimum_image(xac,yac,zac);
103
104 // bond 1->4
105 xad = x[i4][0] - x[i1][0];
106 yad = x[i4][1] - x[i1][1];
107 zad = x[i4][2] - x[i1][2];
108 domain->minimum_image(xad,yad,zad);
109
110 // bond 2-3
111 xbc = x[i3][0] - x[i2][0];
112 ybc = x[i3][1] - x[i2][1];
113 zbc = x[i3][2] - x[i2][2];
114 domain->minimum_image(xbc,ybc,zbc);
115
116 // bond 2-4
117 xbd = x[i4][0] - x[i2][0];
118 ybd = x[i4][1] - x[i2][1];
119 zbd = x[i4][2] - x[i2][2];
120 domain->minimum_image(xbd,ybd,zbd);
121
122 // bond 3-4
123 xcd = x[i4][0] - x[i3][0];
124 ycd = x[i4][1] - x[i3][1];
125 zcd = x[i4][2] - x[i3][2];
126 domain->minimum_image(xcd,ycd,zcd);
127
128 xna = ybc*zcd - zbc*ycd;
129 yna = -(xbc*zcd - zbc*xcd);
130 zna = xbc*ycd - ybc*xcd;
131 rna = 1.0 / sqrt(xna*xna+yna*yna+zna*zna);
132 xna *= rna;
133 yna *= rna;
134 zna *= rna;
135
136 da = -(xna*xad + yna*yad + zna*zad);
137
138
139 domega = k[type]*(da - chi[type])*(da - chi[type]);
140 a = 2.0* k[type]*(da - chi[type]);
141
142 if (eflag) eimproper = domega;
143
144 f1[0] = a*( -xna);
145 f1[1] = a*( -yna);
146 f1[2] = a*( -zna);
147 f4[0] = a*( xna);
148 f4[1] = a*( yna);
149 f4[2] = a*( zna);
150
151 f2[0] = a*( yad*zcd - zad*ycd )*rna + a*da*rna*( yna*zcd - zna*ycd);
152 f2[1] = a*( zad*xcd - xad*zcd )*rna + a*da*rna*( zna*xcd - xna*zcd);
153 f2[2] = a*( xad*ycd - yad*xcd )*rna + a*da*rna*( xna*ycd - yna*xcd);
154
155 f3[0] = - a*( yad*zcd - zad*ycd )*rna - a*da*rna*( yna*zcd - zna*ycd);
156 f3[1] = - a*( zad*xcd - xad*zcd )*rna - a*da*rna*( zna*xcd - xna*zcd);
157 f3[2] = - a*( xad*ycd - yad*xcd )*rna - a*da*rna*( xna*ycd - yna*xcd);
158
159 f3[0] += -a*( yad*zbc - zad*ybc )*rna - a*da*rna*( yna*zbc - zna*ybc);
160 f3[1] += -a*( zad*xbc - xad*zbc )*rna - a*da*rna*( zna*xbc - xna*zbc);
161 f3[2] += -a*( xad*ybc - yad*xbc )*rna - a*da*rna*( xna*ybc - yna*xbc);
162 f4[0] += a*( yad*zbc - zad*ybc )*rna + a*da*rna*( yna*zbc - zna*ybc);
163 f4[1] += a*( zad*xbc - xad*zbc )*rna + a*da*rna*( zna*xbc - xna*zbc);
164 f4[2] += a*( xad*ybc - yad*xbc )*rna + a*da*rna*( xna*ybc - yna*xbc);
165
166 // apply force to each of 4 atoms
167
168 if (newton_bond || i1 < nlocal) {
169 f[i1][0] += f1[0];
170 f[i1][1] += f1[1];
171 f[i1][2] += f1[2];
172 }
173
174 if (newton_bond || i2 < nlocal) {
175 f[i2][0] += f2[0];
176 f[i2][1] += f2[1];
177 f[i2][2] += f2[2];
178 }
179
180 if (newton_bond || i3 < nlocal) {
181 f[i3][0] += f3[0];
182 f[i3][1] += f3[1];
183 f[i3][2] += f3[2];
184 }
185
186 if (newton_bond || i4 < nlocal) {
187 f[i4][0] += f4[0];
188 f[i4][1] += f4[1];
189 f[i4][2] += f4[2];
190 }
191
192 if (evflag)
193 ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f2,f3,f4,
194 xab,yab,zab,xac,yac,zac,xad-xac,yad-yac,zad-zac);
195 }
196 }
197
198 /* ---------------------------------------------------------------------- */
199
allocate()200 void ImproperDistHarm::allocate()
201 {
202 allocated = 1;
203 int n = atom->nimpropertypes;
204
205 memory->create(k,n+1,"improper:k");
206 memory->create(chi,n+1,"improper:chi");
207
208 memory->create(setflag,n+1,"improper:setflag");
209 for (int i = 1; i <= n; i++) setflag[i] = 0;
210 }
211
212 /* ----------------------------------------------------------------------
213 set coeffs for one type
214 ------------------------------------------------------------------------- */
215
coeff(int narg,char ** arg)216 void ImproperDistHarm::coeff(int narg, char **arg)
217 {
218 // if (which > 0) return;
219 if (narg != 3) error->all(FLERR,"Incorrect args for improper coefficients");
220 if (!allocated) allocate();
221
222 int ilo,ihi;
223 utils::bounds(FLERR,arg[0],1,atom->nimpropertypes,ilo,ihi,error);
224
225 double k_one = utils::numeric(FLERR,arg[1],false,lmp);
226 double chi_one = utils::numeric(FLERR,arg[2],false,lmp);
227
228 // convert chi from degrees to radians
229
230 int count = 0;
231 for (int i = ilo; i <= ihi; i++) {
232 k[i] = k_one;
233 //chi[i] = chi_one/180.0 * PI;
234 chi[i] = chi_one;
235 setflag[i] = 1;
236 count++;
237 }
238
239 if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients");
240 }
241
242 /* ----------------------------------------------------------------------
243 proc 0 writes out coeffs to restart file
244 ------------------------------------------------------------------------- */
245
write_restart(FILE * fp)246 void ImproperDistHarm::write_restart(FILE *fp)
247 {
248 fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp);
249 fwrite(&chi[1],sizeof(double),atom->nimpropertypes,fp);
250 }
251
252 /* ----------------------------------------------------------------------
253 proc 0 reads coeffs from restart file, bcasts them
254 ------------------------------------------------------------------------- */
255
read_restart(FILE * fp)256 void ImproperDistHarm::read_restart(FILE *fp)
257 {
258 allocate();
259
260 if (comm->me == 0) {
261 utils::sfread(FLERR,&k[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
262 utils::sfread(FLERR,&chi[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
263 }
264 MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
265 MPI_Bcast(&chi[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
266
267 for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;
268 }
269