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: Paolo Raiteri (Curtin University)
17 ------------------------------------------------------------------------- */
18
19 #include "improper_distance.h"
20
21 #include <cmath>
22 #include "atom.h"
23 #include "comm.h"
24 #include "neighbor.h"
25 #include "domain.h"
26 #include "force.h"
27 #include "memory.h"
28 #include "error.h"
29
30
31 using namespace LAMMPS_NS;
32
33 #define TOLERANCE 0.05
34 #define SMALL 0.001
35
36 /* ---------------------------------------------------------------------- */
37
ImproperDistance(LAMMPS * lmp)38 ImproperDistance::ImproperDistance(LAMMPS *lmp) : Improper(lmp) {}
39
40 /* ---------------------------------------------------------------------- */
41
~ImproperDistance()42 ImproperDistance::~ImproperDistance()
43 {
44 if (allocated) {
45 memory->destroy(setflag);
46 memory->destroy(k);
47 memory->destroy(chi);
48 }
49 }
50
51 /* ---------------------------------------------------------------------- */
52
compute(int eflag,int vflag)53 void ImproperDistance::compute(int eflag, int vflag)
54 {
55 int i1,i2,i3,i4,n,type;
56 double xab, yab, zab; // bond 1-2
57 double xac, yac, zac; // bond 1-3
58 double xad, yad, zad; // bond 1-4
59 double xbc, ybc, zbc; // bond 2-3
60 double xbd, ybd, zbd; // bond 2-4
61 double xna, yna, zna, rna; // normal
62 double da;
63
64 double eimproper,f1[3],f2[3],f3[3],f4[3];
65 // double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2;
66 // double s12,c,s,domega,a,a11,a22,a33,a12,a13,a23;
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 // 1 is the central atom
88 // 2-3-4 are ment to be equivalent
89 // I need the bonds between 2-3 and 2-4 to get the plane normal
90 // Then I need the bond 1-2 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 xna = ybc*zbd - zbc*ybd;
123 yna = -(xbc*zbd - zbc*xbd);
124 zna = xbc*ybd - ybc*xbd;
125 rna = 1.0 / sqrt(xna*xna+yna*yna+zna*zna);
126 xna *= rna;
127 yna *= rna;
128 zna *= rna;
129
130 da = xna*xab + yna*yab + zna*zab;
131
132 domega = k[type]*da*da + chi[type]*da*da*da*da;
133 //printf("%3i %3i %3i %3i %10.5f %10.5f \n",i1,i2,i3,i4,da,domega);
134 a = 2.0* (k[type]*da + 2.0*chi[type]*da*da*da);
135
136 if (eflag) eimproper = domega;
137
138 f1[0] = a*( xna);
139 f1[1] = a*( yna);
140 f1[2] = a*( zna);
141
142 f2[0] = a*( -xna -yab*(zbd-zbc)*rna +zab*(ybd-ybc)*rna -da*( -yna*(zbd-zbc) + zna*(ybd-ybc) )*rna);
143 f2[1] = a*( +xab*(zbd-zbc)*rna -yna +zab*(xbc-xbd)*rna -da*( +xna*(zbd-zbc) + zna*(xbc-xbd) )*rna);
144 f2[2] = a*( -xab*(ybd-ybc)*rna -yab*(xbc-xbd)*rna -zna -da*( +xna*(ybc-ybd) - yna*(xbc-xbd) )*rna);
145
146 f3[0] = a*( ( yab*zbd -zab*ybd ) *rna +da*( -yna*zbd +zna*ybd )*rna);
147 f3[1] = a*( ( -xab*zbd +zab*xbd ) *rna +da*( +xna*zbd -zna*xbd )*rna);
148 f3[2] = a*( ( +xab*ybd -yab*xbd ) *rna +da*( -xna*ybd +yna*xbd )*rna);
149
150 f4[0] = a*( ( -yab*zbc +zab*ybc ) *rna -da*( -yna*zbc +zna*ybc )*rna);
151 f4[1] = a*( ( +xab*zbc -zab*xbc ) *rna -da*( +xna*zbc -zna*xbc )*rna);
152 f4[2] = a*( ( -xab*ybc +yab*xbc ) *rna -da*( -xna*ybc +yna*xbc )*rna);
153 //printf("%10.5f %10.5f %10.5f \n",f1[0],f1[1],f1[2]);
154 //printf("%10.5f %10.5f %10.5f \n",f2[0],f2[1],f2[2]);
155 //printf("%10.5f %10.5f %10.5f \n",f3[0],f3[1],f3[2]);
156 //printf("%10.5f %10.5f %10.5f \n",f4[0],f4[1],f4[2]);
157
158 // apply force to each of 4 atoms
159
160 if (newton_bond || i1 < nlocal) {
161 f[i1][0] += f1[0];
162 f[i1][1] += f1[1];
163 f[i1][2] += f1[2];
164 }
165
166 if (newton_bond || i2 < nlocal) {
167 f[i2][0] += f2[0];
168 f[i2][1] += f2[1];
169 f[i2][2] += f2[2];
170 }
171
172 if (newton_bond || i3 < nlocal) {
173 f[i3][0] += f3[0];
174 f[i3][1] += f3[1];
175 f[i3][2] += f3[2];
176 }
177
178 if (newton_bond || i4 < nlocal) {
179 f[i4][0] += f4[0];
180 f[i4][1] += f4[1];
181 f[i4][2] += f4[2];
182 }
183
184 if (evflag)
185 ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f2,f3,f4,
186 xab,yab,zab,xac,yac,zac,xad-xac,yad-yac,zad-zac);
187 }
188 }
189
190 /* ---------------------------------------------------------------------- */
191
allocate()192 void ImproperDistance::allocate()
193 {
194 allocated = 1;
195 int n = atom->nimpropertypes;
196
197 memory->create(k,n+1,"improper:k");
198 memory->create(chi,n+1,"improper:chi");
199
200 memory->create(setflag,n+1,"improper:setflag");
201 for (int i = 1; i <= n; i++) setflag[i] = 0;
202 }
203
204 /* ----------------------------------------------------------------------
205 set coeffs for one type
206 ------------------------------------------------------------------------- */
207
coeff(int narg,char ** arg)208 void ImproperDistance::coeff(int narg, char **arg)
209 {
210 // if (which > 0) return;
211 if (narg != 3) error->all(FLERR,"Incorrect args for improper coefficients");
212 if (!allocated) allocate();
213
214 int ilo,ihi;
215 utils::bounds(FLERR,arg[0],1,atom->nimpropertypes,ilo,ihi,error);
216
217 double k_one = utils::numeric(FLERR,arg[1],false,lmp);
218 double chi_one = utils::numeric(FLERR,arg[2],false,lmp);
219
220 // convert chi from degrees to radians
221
222 int count = 0;
223 for (int i = ilo; i <= ihi; i++) {
224 k[i] = k_one;
225 chi[i] = chi_one;
226 setflag[i] = 1;
227 count++;
228 }
229
230 if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients");
231 }
232
233 /* ----------------------------------------------------------------------
234 proc 0 writes out coeffs to restart file
235 ------------------------------------------------------------------------- */
236
write_restart(FILE * fp)237 void ImproperDistance::write_restart(FILE *fp)
238 {
239 fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp);
240 fwrite(&chi[1],sizeof(double),atom->nimpropertypes,fp);
241 }
242
243 /* ----------------------------------------------------------------------
244 proc 0 reads coeffs from restart file, bcasts them
245 ------------------------------------------------------------------------- */
246
read_restart(FILE * fp)247 void ImproperDistance::read_restart(FILE *fp)
248 {
249 allocate();
250
251 if (comm->me == 0) {
252 utils::sfread(FLERR,&k[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
253 utils::sfread(FLERR,&chi[1],sizeof(double),atom->nimpropertypes,fp,nullptr,error);
254 }
255 MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
256 MPI_Bcast(&chi[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
257
258 for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;
259 }
260
261 /* ----------------------------------------------------------------------
262 proc 0 writes to data file
263 ------------------------------------------------------------------------- */
264
write_data(FILE * fp)265 void ImproperDistance::write_data(FILE *fp)
266 {
267 for (int i = 1; i <= atom->nimpropertypes; i++)
268 fprintf(fp,"%d %g %g\n",i,k[i],chi[i]);
269 }
270