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_temp_csld.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "compute.h"
24 #include "error.h"
25 #include "force.h"
26 #include "group.h"
27 #include "input.h"
28 #include "memory.h"
29 #include "modify.h"
30 #include "random_mars.h"
31 #include "update.h"
32 #include "variable.h"
33
34 #include <cmath>
35 #include <cstring>
36
37 using namespace LAMMPS_NS;
38 using namespace FixConst;
39
40 enum{NOBIAS,BIAS};
41 enum{CONSTANT,EQUAL};
42
43 /* ---------------------------------------------------------------------- */
44
FixTempCSLD(LAMMPS * lmp,int narg,char ** arg)45 FixTempCSLD::FixTempCSLD(LAMMPS *lmp, int narg, char **arg) :
46 Fix(lmp, narg, arg),
47 vhold(nullptr), tstr(nullptr), id_temp(nullptr), random(nullptr)
48 {
49 if (narg != 7) error->all(FLERR,"Illegal fix temp/csld command");
50
51 // CSLD thermostat should be applied every step
52
53 restart_global = 1;
54 nevery = 1;
55 scalar_flag = 1;
56 ecouple_flag = 1;
57 global_freq = nevery;
58 dynamic_group_allow = 1;
59 extscalar = 1;
60
61 tstr = nullptr;
62 if (utils::strmatch(arg[3],"^v_")) {
63 tstr = utils::strdup(arg[3]+2);
64 tstyle = EQUAL;
65 } else {
66 t_start = utils::numeric(FLERR,arg[3],false,lmp);
67 t_target = t_start;
68 tstyle = CONSTANT;
69 }
70
71 t_stop = utils::numeric(FLERR,arg[4],false,lmp);
72 t_period = utils::numeric(FLERR,arg[5],false,lmp);
73 int seed = utils::inumeric(FLERR,arg[6],false,lmp);
74
75 // error checks
76
77 if (t_period <= 0.0) error->all(FLERR,"Illegal fix temp/csld command");
78 if (seed <= 0) error->all(FLERR,"Illegal fix temp/csld command");
79
80 random = new RanMars(lmp,seed + comm->me);
81
82 // create a new compute temp style
83 // id = fix-ID + temp, compute group = fix group
84
85 id_temp = utils::strdup(std::string(id) + "_temp");
86 modify->add_compute(fmt::format("{} {} temp",id_temp,group->names[igroup]));
87 tflag = 1;
88
89 vhold = nullptr;
90 nmax = -1;
91 energy = 0.0;
92 }
93
94 /* ---------------------------------------------------------------------- */
95
~FixTempCSLD()96 FixTempCSLD::~FixTempCSLD()
97 {
98 delete [] tstr;
99
100 // delete temperature if fix created it
101
102 if (tflag) modify->delete_compute(id_temp);
103 delete [] id_temp;
104
105 delete random;
106 memory->destroy(vhold);
107 vhold = nullptr;
108 nmax = -1;
109 }
110
111 /* ---------------------------------------------------------------------- */
112
setmask()113 int FixTempCSLD::setmask()
114 {
115 int mask = 0;
116 mask |= END_OF_STEP;
117 return mask;
118 }
119
120 /* ---------------------------------------------------------------------- */
121
init()122 void FixTempCSLD::init()
123 {
124
125 // we cannot handle constraints via rattle or shake correctly.
126
127 int has_shake = 0;
128 for (int i = 0; i < modify->nfix; i++)
129 if ((strcmp(modify->fix[i]->style,"shake") == 0)
130 || (strcmp(modify->fix[i]->style,"rattle") == 0)) ++has_shake;
131
132 if (has_shake > 0)
133 error->all(FLERR,"Fix temp/csld is not compatible with fix rattle or fix shake");
134
135 // check variable
136
137 if (tstr) {
138 tvar = input->variable->find(tstr);
139 if (tvar < 0)
140 error->all(FLERR,"Variable name for fix temp/csld does not exist");
141 if (input->variable->equalstyle(tvar)) tstyle = EQUAL;
142 else error->all(FLERR,"Variable for fix temp/csld is invalid style");
143 }
144
145 int icompute = modify->find_compute(id_temp);
146 if (icompute < 0)
147 error->all(FLERR,"Temperature ID for fix temp/csld does not exist");
148 temperature = modify->compute[icompute];
149
150 if (modify->check_rigid_group_overlap(groupbit))
151 error->warning(FLERR,"Cannot thermostat atoms in rigid bodies");
152
153 if (temperature->tempbias) which = BIAS;
154 else which = NOBIAS;
155 }
156
157 /* ---------------------------------------------------------------------- */
158
end_of_step()159 void FixTempCSLD::end_of_step()
160 {
161
162 // set current t_target
163 // if variable temp, evaluate variable, wrap with clear/add
164
165 double delta = update->ntimestep - update->beginstep;
166
167 if (delta != 0.0) delta /= update->endstep - update->beginstep;
168 if (tstyle == CONSTANT)
169 t_target = t_start + delta * (t_stop-t_start);
170 else {
171 modify->clearstep_compute();
172 t_target = input->variable->compute_equal(tvar);
173 if (t_target < 0.0)
174 error->one(FLERR,
175 "Fix temp/csld variable returned negative temperature");
176 modify->addstep_compute(update->ntimestep + nevery);
177 }
178
179 double t_current = temperature->compute_scalar();
180 double ekin_old = t_current * 0.5 * temperature->dof * force->boltz;
181
182 // there is nothing to do, if there are no degrees of freedom
183
184 if (temperature->dof < 1) return;
185
186 double * const * const v = atom->v;
187 const int * const mask = atom->mask;
188 const int * const type = atom->type;
189 const int nlocal = atom->nlocal;
190
191 // adjust holding space, if needed and copy existing velocities
192
193 if (nmax < atom->nlocal) {
194 nmax = atom->nlocal + 1;
195 memory->destroy(vhold);
196 memory->create(vhold,nmax,3,"csld:vhold");
197 }
198
199 // The CSLD thermostat is a linear combination of old and new velocities,
200 // where the new ones are randomly chosen from a gaussian distribution.
201 // see Bussi and Parrinello, Phys. Rev. E (2007).
202
203 for (int i = 0; i < nlocal; i++) {
204 if (mask[i] & groupbit) {
205 double m;
206 if (atom->rmass_flag) m = atom->rmass[i];
207 else m = atom->mass[type[i]];
208
209 const double factor = 1.0/sqrt(m);
210 const double vx = random->gaussian() * factor;
211 vhold[i][0] = v[i][0];
212 v[i][0] = vx;
213 const double vy = random->gaussian() * factor;
214 vhold[i][1] = v[i][1];
215 v[i][1] = vy;
216 const double vz = random->gaussian() * factor;
217 vhold[i][2] = v[i][2];
218 v[i][2] = vz;
219 }
220 }
221
222 // mixing factors
223 const double c1 = exp(-update->dt/t_period);
224 const double c2 = sqrt((1.0-c1*c1)*t_target/temperature->compute_scalar());
225
226 if (which == NOBIAS) {
227 for (int i = 0; i < nlocal; i++) {
228 if (mask[i] & groupbit) {
229 v[i][0] = vhold[i][0]*c1 + v[i][0]*c2;
230 v[i][1] = vhold[i][1]*c1 + v[i][1]*c2;
231 v[i][2] = vhold[i][2]*c1 + v[i][2]*c2;
232 }
233 }
234 } else {
235 for (int i = 0; i < nlocal; i++) {
236 if (mask[i] & groupbit) {
237 temperature->remove_bias(i,vhold[i]);
238 v[i][0] = vhold[i][0]*c1 + v[i][0]*c2;
239 v[i][1] = vhold[i][1]*c1 + v[i][1]*c2;
240 v[i][2] = vhold[i][2]*c1 + v[i][2]*c2;
241 temperature->restore_bias(i,v[i]);
242 }
243 }
244 }
245
246 // tally the kinetic energy transferred between heat bath and system
247
248 t_current = temperature->compute_scalar();
249 energy += ekin_old - t_current * 0.5 * temperature->dof * force->boltz;
250 }
251
252 /* ---------------------------------------------------------------------- */
253
modify_param(int narg,char ** arg)254 int FixTempCSLD::modify_param(int narg, char **arg)
255 {
256 if (strcmp(arg[0],"temp") == 0) {
257 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
258 if (tflag) {
259 modify->delete_compute(id_temp);
260 tflag = 0;
261 }
262 delete [] id_temp;
263 id_temp = utils::strdup(arg[1]);
264
265 int icompute = modify->find_compute(id_temp);
266 if (icompute < 0)
267 error->all(FLERR,"Could not find fix_modify temperature ID");
268 temperature = modify->compute[icompute];
269
270 if (temperature->tempflag == 0)
271 error->all(FLERR,
272 "Fix_modify temperature ID does not compute temperature");
273 if (temperature->igroup != igroup && comm->me == 0)
274 error->warning(FLERR,"Group for fix_modify temp != fix group");
275 return 2;
276 }
277 return 0;
278 }
279
280 /* ---------------------------------------------------------------------- */
281
reset_target(double t_new)282 void FixTempCSLD::reset_target(double t_new)
283 {
284 t_target = t_start = t_stop = t_new;
285 }
286
287 /* ---------------------------------------------------------------------- */
288
compute_scalar()289 double FixTempCSLD::compute_scalar()
290 {
291 return energy;
292 }
293
294 /* ----------------------------------------------------------------------
295 pack entire state of Fix into one write
296 ------------------------------------------------------------------------- */
297
write_restart(FILE * fp)298 void FixTempCSLD::write_restart(FILE *fp)
299 {
300 const int PRNGSIZE = 98+2+3;
301 int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy
302 double *list = nullptr;
303 if (comm->me == 0) {
304 list = new double[nsize];
305 list[0] = energy;
306 list[1] = comm->nprocs;
307 }
308 double state[PRNGSIZE];
309 random->get_state(state);
310 MPI_Gather(state,PRNGSIZE,MPI_DOUBLE,list+2,PRNGSIZE,MPI_DOUBLE,0,world);
311
312 if (comm->me == 0) {
313 int size = nsize * sizeof(double);
314 fwrite(&size,sizeof(int),1,fp);
315 fwrite(list,sizeof(double),nsize,fp);
316 delete[] list;
317 }
318 }
319
320 /* ----------------------------------------------------------------------
321 use state info from restart file to restart the Fix
322 ------------------------------------------------------------------------- */
323
restart(char * buf)324 void FixTempCSLD::restart(char *buf)
325 {
326 double *list = (double *) buf;
327
328 energy = list[0];
329 int nprocs = (int) list[1];
330 if (nprocs != comm->nprocs) {
331 if (comm->me == 0)
332 error->warning(FLERR,"Different number of procs. Cannot restore RNG state.");
333 } else random->set_state(list+2+comm->me*103);
334 }
335
336 /* ----------------------------------------------------------------------
337 extract thermostat properties
338 ------------------------------------------------------------------------- */
339
extract(const char * str,int & dim)340 void *FixTempCSLD::extract(const char *str, int &dim)
341 {
342 dim=0;
343 if (strcmp(str,"t_target") == 0) {
344 return &t_target;
345 }
346 return nullptr;
347 }
348