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_viscous.h"
16
17 #include "atom.h"
18 #include "error.h"
19 #include "respa.h"
20 #include "update.h"
21
22 #include <cstring>
23
24 using namespace LAMMPS_NS;
25 using namespace FixConst;
26
27 /* ---------------------------------------------------------------------- */
28
FixViscous(LAMMPS * lmp,int narg,char ** arg)29 FixViscous::FixViscous(LAMMPS *lmp, int narg, char **arg) :
30 Fix(lmp, narg, arg),
31 gamma(nullptr)
32 {
33 dynamic_group_allow = 1;
34
35 if (narg < 4) error->all(FLERR,"Illegal fix viscous command");
36
37 double gamma_one = utils::numeric(FLERR,arg[3],false,lmp);
38 gamma = new double[atom->ntypes+1];
39 for (int i = 1; i <= atom->ntypes; i++) gamma[i] = gamma_one;
40
41 // optional args
42
43 int iarg = 4;
44 while (iarg < narg) {
45 if (strcmp(arg[iarg],"scale") == 0) {
46 if (iarg+3 > narg) error->all(FLERR,"Illegal fix viscous command");
47 int itype = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
48 double scale = utils::numeric(FLERR,arg[iarg+2],false,lmp);
49 if (itype <= 0 || itype > atom->ntypes)
50 error->all(FLERR,"Illegal fix viscous command");
51 gamma[itype] = gamma_one * scale;
52 iarg += 3;
53 } else error->all(FLERR,"Illegal fix viscous command");
54 }
55
56 respa_level_support = 1;
57 ilevel_respa = 0;
58 }
59
60 /* ---------------------------------------------------------------------- */
61
~FixViscous()62 FixViscous::~FixViscous()
63 {
64 delete [] gamma;
65 }
66
67 /* ---------------------------------------------------------------------- */
68
setmask()69 int FixViscous::setmask()
70 {
71 int mask = 0;
72 mask |= POST_FORCE;
73 mask |= POST_FORCE_RESPA;
74 mask |= MIN_POST_FORCE;
75 return mask;
76 }
77
78 /* ---------------------------------------------------------------------- */
79
init()80 void FixViscous::init()
81 {
82 int max_respa = 0;
83
84 if (utils::strmatch(update->integrate_style,"^respa")) {
85 ilevel_respa = max_respa = ((Respa *) update->integrate)->nlevels-1;
86 if (respa_level >= 0) ilevel_respa = MIN(respa_level,max_respa);
87 }
88 }
89
90 /* ---------------------------------------------------------------------- */
91
setup(int vflag)92 void FixViscous::setup(int vflag)
93 {
94 if (utils::strmatch(update->integrate_style,"^verlet"))
95 post_force(vflag);
96 else {
97 ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
98 post_force_respa(vflag,ilevel_respa,0);
99 ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
100 }
101 }
102
103 /* ---------------------------------------------------------------------- */
104
min_setup(int vflag)105 void FixViscous::min_setup(int vflag)
106 {
107 post_force(vflag);
108 }
109
110 /* ---------------------------------------------------------------------- */
111
post_force(int)112 void FixViscous::post_force(int /*vflag*/)
113 {
114 // apply drag force to atoms in group
115 // direction is opposed to velocity vector
116 // magnitude depends on atom type
117
118 double **v = atom->v;
119 double **f = atom->f;
120 int *mask = atom->mask;
121 int *type = atom->type;
122 int nlocal = atom->nlocal;
123
124 double drag;
125
126 for (int i = 0; i < nlocal; i++)
127 if (mask[i] & groupbit) {
128 drag = gamma[type[i]];
129 f[i][0] -= drag*v[i][0];
130 f[i][1] -= drag*v[i][1];
131 f[i][2] -= drag*v[i][2];
132 }
133 }
134
135 /* ---------------------------------------------------------------------- */
136
post_force_respa(int vflag,int ilevel,int)137 void FixViscous::post_force_respa(int vflag, int ilevel, int /*iloop*/)
138 {
139 if (ilevel == ilevel_respa) post_force(vflag);
140 }
141
142 /* ---------------------------------------------------------------------- */
143
min_post_force(int vflag)144 void FixViscous::min_post_force(int vflag)
145 {
146 post_force(vflag);
147 }
148