1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/
5 Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
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_nvt_sllod_omp.h"
20
21 #include "atom.h"
22 #include "compute.h"
23 #include "domain.h"
24 #include "error.h"
25 #include "fix.h"
26 #include "fix_deform.h"
27 #include "group.h"
28 #include "math_extra.h"
29 #include "modify.h"
30
31 #include <cstring>
32
33 #include "omp_compat.h"
34
35 using namespace LAMMPS_NS;
36 using namespace FixConst;
37
38 typedef struct { double x,y,z; } dbl3_t;
39
40 /* ---------------------------------------------------------------------- */
41
FixNVTSllodOMP(LAMMPS * lmp,int narg,char ** arg)42 FixNVTSllodOMP::FixNVTSllodOMP(LAMMPS *lmp, int narg, char **arg) :
43 FixNHOMP(lmp, narg, arg)
44 {
45 if (!tstat_flag)
46 error->all(FLERR,"Temperature control must be used with fix nvt/sllod");
47 if (pstat_flag)
48 error->all(FLERR,"Pressure control can not be used with fix nvt/sllod");
49
50 // default values
51
52 if (mtchain_default_flag) mtchain = 1;
53
54
55 // create a new compute temp style
56 // id = fix-ID + temp
57
58 id_temp = utils::strdup(std::string(id) + "_temp");
59 modify->add_compute(fmt::format("{} {} temp/deform",
60 id_temp,group->names[igroup]));
61 tcomputeflag = 1;
62 }
63
64 /* ---------------------------------------------------------------------- */
65
init()66 void FixNVTSllodOMP::init()
67 {
68 FixNHOMP::init();
69
70 if (!temperature->tempbias)
71 error->all(FLERR,"Temperature for fix nvt/sllod/omp does not have a bias");
72
73 nondeformbias = 0;
74 if (strcmp(temperature->style,"temp/deform") != 0) nondeformbias = 1;
75
76 // check fix deform remap settings
77
78 int i;
79 for (i = 0; i < modify->nfix; i++)
80 if (utils::strmatch(modify->fix[i]->style,"^deform")) {
81 if (((FixDeform *) modify->fix[i])->remapflag != Domain::V_REMAP)
82 error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix "
83 "deform remap option");
84 break;
85 }
86 if (i == modify->nfix)
87 error->all(FLERR,"Using fix nvt/sllod/omp with no fix deform defined");
88 }
89
90 /* ----------------------------------------------------------------------
91 perform half-step scaling of velocities
92 -----------------------------------------------------------------------*/
93
nh_v_temp()94 void FixNVTSllodOMP::nh_v_temp()
95 {
96 // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo
97 // thermostat thermal velocity only
98 // vdelu = SLLOD correction = Hrate*Hinv*vthermal
99 // for non temp/deform BIAS:
100 // calculate temperature since some computes require temp
101 // computed on current nlocal atoms to remove bias
102
103 dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
104 const int * _noalias const mask = atom->mask;
105 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
106
107 if (nondeformbias) temperature->compute_scalar();
108
109 double h_two[6];
110 MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two);
111
112 #if defined(_OPENMP)
113 #pragma omp parallel for LMP_DEFAULT_NONE LMP_SHARED(h_two) schedule(static)
114 #endif
115 for (int i = 0; i < nlocal; i++) {
116 double vdelu0,vdelu1,vdelu2,buf[3];
117 if (mask[i] & groupbit) {
118 vdelu0 = h_two[0]*v[i].x + h_two[5]*v[i].y + h_two[4]*v[i].z;
119 vdelu1 = h_two[1]*v[i].y + h_two[3]*v[i].z;
120 vdelu2 = h_two[2]*v[i].z;
121 temperature->remove_bias_thr(i,&v[i].x,buf);
122 v[i].x = v[i].x*factor_eta - dthalf*vdelu0;
123 v[i].y = v[i].y*factor_eta - dthalf*vdelu1;
124 v[i].z = v[i].z*factor_eta - dthalf*vdelu2;
125 temperature->restore_bias_thr(i,&v[i].x,buf);
126 }
127 }
128 }
129