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)
17 ------------------------------------------------------------------------- */
18
19 #include "fix_nh_asphere_omp.h"
20
21 #include "atom.h"
22 #include "atom_vec_ellipsoid.h"
23 #include "compute.h"
24 #include "error.h"
25 #include "math_extra.h"
26
27 #include "omp_compat.h"
28 using namespace LAMMPS_NS;
29 using namespace FixConst;
30
31 enum{NOBIAS,BIAS};
32
33 typedef struct { double x,y,z; } dbl3_t;
34
35 /* ---------------------------------------------------------------------- */
36
FixNHAsphereOMP(LAMMPS * lmp,int narg,char ** arg)37 FixNHAsphereOMP::FixNHAsphereOMP(LAMMPS *lmp, int narg, char **arg) :
38 FixNHOMP(lmp, narg, arg)
39 {
40 }
41
42 /* ---------------------------------------------------------------------- */
43
init()44 void FixNHAsphereOMP::init()
45 {
46 avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
47 if (!avec)
48 error->all(FLERR,"Compute nvt/nph/npt asphere requires atom style ellipsoid");
49
50 // check that all particles are finite-size
51 // no point particles allowed, spherical is OK
52
53 int *ellipsoid = atom->ellipsoid;
54 int *mask = atom->mask;
55 int nlocal = atom->nlocal;
56
57 for (int i = 0; i < nlocal; i++)
58 if (mask[i] & groupbit)
59 if (ellipsoid[i] < 0)
60 error->one(FLERR,"Fix nvt/nph/npt asphere requires extended particles");
61
62 FixNHOMP::init();
63 }
64
65 /* ----------------------------------------------------------------------
66 perform half-step update of angular momentum and COM velocity
67 -----------------------------------------------------------------------*/
68
nve_v()69 void FixNHAsphereOMP::nve_v()
70 {
71 dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
72 dbl3_t * _noalias const angmom = (dbl3_t *) atom->angmom[0];
73 const dbl3_t * _noalias const f = (dbl3_t *) atom->f[0];
74 const dbl3_t * _noalias const torque = (dbl3_t *) atom->torque[0];
75 const double * _noalias const rmass = atom->rmass;
76 const int * _noalias const mask = atom->mask;
77 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
78
79 // standard nve_v velocity update. for efficiency the loop is
80 // merged with FixNHOMP instead of calling it for the COM update.
81
82 #if defined(_OPENMP)
83 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
84 #endif
85 for (int i = 0; i < nlocal; i++) {
86 if (mask[i] & groupbit) {
87 const double dtfm = dtf / rmass[i];
88 v[i].x += dtfm*f[i].x;
89 v[i].y += dtfm*f[i].y;
90 v[i].z += dtfm*f[i].z;
91 angmom[i].x += dtf*torque[i].x;
92 angmom[i].y += dtf*torque[i].y;
93 angmom[i].z += dtf*torque[i].z;
94 }
95 }
96 }
97
98 /* ----------------------------------------------------------------------
99 perform full-step update of position and orientation
100 -----------------------------------------------------------------------*/
101
nve_x()102 void FixNHAsphereOMP::nve_x()
103 {
104 dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
105 const dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
106 dbl3_t * _noalias const angmom = (dbl3_t *) atom->angmom[0];
107 const double * _noalias const rmass = atom->rmass;
108 const int * _noalias const mask = atom->mask;
109 AtomVecEllipsoid::Bonus * _noalias const bonus = avec->bonus;
110 const int * _noalias const ellipsoid = atom->ellipsoid;
111 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
112
113 // set timestep here since dt may have changed or come via rRESPA
114
115 dtq = 0.5 * dtv;
116
117 // standard nve_x position update and
118 // update quaternion a full step via Richardson iteration
119 // returns new normalized quaternion
120 // principal moments of inertia
121
122 #if defined(_OPENMP)
123 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
124 #endif
125 for (int i = 0; i < nlocal; i++)
126 if (mask[i] & groupbit) {
127 double omega[3], inertia[3];
128
129 x[i].x += dtv * v[i].x;
130 x[i].y += dtv * v[i].y;
131 x[i].z += dtv * v[i].z;
132
133 // principal moments of inertia
134
135 const double * const shape = bonus[ellipsoid[i]].shape;
136 double * const quat = bonus[ellipsoid[i]].quat;
137
138 inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
139 inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
140 inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
141
142 // compute omega at 1/2 step from angmom at 1/2 step and current q
143 // update quaternion a full step via Richardson iteration
144 // returns new normalized quaternion
145
146 MathExtra::mq_to_omega(&angmom[i].x,quat,inertia,omega);
147 MathExtra::richardson(quat,&angmom[i].x,omega,inertia,dtq);
148 }
149 }
150
151 /* ----------------------------------------------------------------------
152 perform half-step temperature scaling of angular momentum
153 -----------------------------------------------------------------------*/
154
nh_v_temp()155 void FixNHAsphereOMP::nh_v_temp()
156 {
157 dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
158 dbl3_t * _noalias const angmom = (dbl3_t *) atom->angmom[0];
159 const int * _noalias const mask = atom->mask;
160 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
161
162 if (which == NOBIAS) {
163 #if defined(_OPENMP)
164 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
165 #endif
166 for (int i = 0; i < nlocal; i++) {
167 if (mask[i] & groupbit) {
168 v[i].x *= factor_eta;
169 v[i].y *= factor_eta;
170 v[i].z *= factor_eta;
171 angmom[i].x *= factor_eta;
172 angmom[i].y *= factor_eta;
173 angmom[i].z *= factor_eta;
174 }
175 }
176 } else if (which == BIAS) {
177 #if defined(_OPENMP)
178 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
179 #endif
180 for (int i = 0; i < nlocal; i++) {
181 double buf[3];
182 if (mask[i] & groupbit) {
183 temperature->remove_bias_thr(i,&v[i].x,buf);
184 v[i].x *= factor_eta;
185 v[i].y *= factor_eta;
186 v[i].z *= factor_eta;
187 temperature->restore_bias_thr(i,&v[i].x,buf);
188 angmom[i].x *= factor_eta;
189 angmom[i].y *= factor_eta;
190 angmom[i].z *= factor_eta;
191 }
192 }
193 }
194 }
195