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 "fix_nve_tri.h"
16 #include "math_extra.h"
17 #include "atom.h"
18 #include "atom_vec_tri.h"
19 #include "domain.h"
20 #include "error.h"
21
22 using namespace LAMMPS_NS;
23 using namespace FixConst;
24
25 /* ---------------------------------------------------------------------- */
26
FixNVETri(LAMMPS * lmp,int narg,char ** arg)27 FixNVETri::FixNVETri(LAMMPS *lmp, int narg, char **arg) :
28 FixNVE(lmp, narg, arg)
29 {
30 if (narg != 3) error->all(FLERR,"Illegal fix nve/tri command");
31
32 time_integrate = 1;
33 }
34
35 /* ---------------------------------------------------------------------- */
36
setmask()37 int FixNVETri::setmask()
38 {
39 int mask = 0;
40 mask |= INITIAL_INTEGRATE;
41 mask |= FINAL_INTEGRATE;
42 mask |= INITIAL_INTEGRATE_RESPA;
43 mask |= FINAL_INTEGRATE_RESPA;
44 return mask;
45 }
46
47 /* ---------------------------------------------------------------------- */
48
init()49 void FixNVETri::init()
50 {
51 // error checks
52
53 avec = (AtomVecTri *) atom->style_match("tri");
54 if (!avec) error->all(FLERR,"Fix nve/tri requires atom style tri");
55
56 if (domain->dimension != 3)
57 error->all(FLERR,"Fix nve/tri can only be used for 3d simulations");
58
59 // check that all particles are triangles
60 // no point particles allowed
61
62 int *tri = atom->tri;
63 int *mask = atom->mask;
64 int nlocal = atom->nlocal;
65
66 for (int i = 0; i < nlocal; i++)
67 if (mask[i] & groupbit) {
68 if (tri[i] < 0) error->one(FLERR,"Fix nve/tri requires tri particles");
69 }
70
71 FixNVE::init();
72 }
73
74 /* ---------------------------------------------------------------------- */
75
initial_integrate(int)76 void FixNVETri::initial_integrate(int /*vflag*/)
77 {
78 double dtfm;
79 double omega[3];
80
81 AtomVecTri::Bonus *bonus = avec->bonus;
82 int *tri = atom->tri;
83 double **x = atom->x;
84 double **v = atom->v;
85 double **f = atom->f;
86 double **angmom = atom->angmom;
87 double **torque = atom->torque;
88 double *rmass = atom->rmass;
89 int *mask = atom->mask;
90 int nlocal = atom->nlocal;
91 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
92
93 // set timestep here since dt may have changed or come via rRESPA
94
95 dtq = 0.5 * dtv;
96
97 for (int i = 0; i < nlocal; i++)
98 if (mask[i] & groupbit) {
99 dtfm = dtf / rmass[i];
100 v[i][0] += dtfm * f[i][0];
101 v[i][1] += dtfm * f[i][1];
102 v[i][2] += dtfm * f[i][2];
103 x[i][0] += dtv * v[i][0];
104 x[i][1] += dtv * v[i][1];
105 x[i][2] += dtv * v[i][2];
106
107 // update angular momentum by 1/2 step
108
109 angmom[i][0] += dtf * torque[i][0];
110 angmom[i][1] += dtf * torque[i][1];
111 angmom[i][2] += dtf * torque[i][2];
112
113 // compute omega at 1/2 step from angmom at 1/2 step and current q
114 // update quaternion a full step via Richardson iteration
115 // returns new normalized quaternion
116
117 MathExtra::mq_to_omega(angmom[i],bonus[tri[i]].quat,
118 bonus[tri[i]].inertia,omega);
119 MathExtra::richardson(bonus[tri[i]].quat,angmom[i],omega,
120 bonus[tri[i]].inertia,dtq);
121 }
122 }
123
124 /* ---------------------------------------------------------------------- */
125
final_integrate()126 void FixNVETri::final_integrate()
127 {
128 double dtfm;
129
130 double **v = atom->v;
131 double **f = atom->f;
132 double **angmom = atom->angmom;
133 double **torque = atom->torque;
134 double *rmass = atom->rmass;
135 int *mask = atom->mask;
136 int nlocal = atom->nlocal;
137 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
138
139 // update v,omega for all particles
140 // d_omega/dt = torque / inertia
141
142 for (int i = 0; i < nlocal; i++)
143 if (mask[i] & groupbit) {
144 dtfm = dtf / rmass[i];
145 v[i][0] += dtfm * f[i][0];
146 v[i][1] += dtfm * f[i][1];
147 v[i][2] += dtfm * f[i][2];
148
149 angmom[i][0] += dtf * torque[i][0];
150 angmom[i][1] += dtf * torque[i][1];
151 angmom[i][2] += dtf * torque[i][2];
152 }
153 }
154