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 Based on code by Paolo Raiteri (Curtin U) and Giovanni Bussi (SISSA)
18 ------------------------------------------------------------------------- */
19
20 #include "fix_temp_csvr.h"
21
22 #include "atom.h"
23 #include "comm.h"
24 #include "compute.h"
25 #include "error.h"
26 #include "force.h"
27 #include "group.h"
28 #include "input.h"
29 #include "modify.h"
30 #include "random_mars.h"
31 #include "update.h"
32 #include "variable.h"
33
34 #include <cstring>
35 #include <cmath>
36
37 using namespace LAMMPS_NS;
38 using namespace FixConst;
39
40 enum{NOBIAS,BIAS};
41 enum{CONSTANT,EQUAL};
42
43 static constexpr int PRNGSIZE = 98+2+3;
44 /* ---------------------------------------------------------------------- */
45
FixTempCSVR(LAMMPS * lmp,int narg,char ** arg)46 FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
47 Fix(lmp, narg, arg),
48 tstr(nullptr), id_temp(nullptr), random(nullptr)
49 {
50 if (narg != 7) error->all(FLERR,"Illegal fix temp/csvr command");
51
52 // CSVR thermostat should be applied every step
53
54 restart_global = 1;
55 nevery = 1;
56 scalar_flag = 1;
57 ecouple_flag = 1;
58 global_freq = nevery;
59 dynamic_group_allow = 1;
60 extscalar = 1;
61
62 tstr = nullptr;
63 if (utils::strmatch(arg[3],"^v_")) {
64 tstr = utils::strdup(arg[3]+2);
65 tstyle = EQUAL;
66 } else {
67 t_start = utils::numeric(FLERR,arg[3],false,lmp);
68 t_target = t_start;
69 tstyle = CONSTANT;
70 }
71
72 t_stop = utils::numeric(FLERR,arg[4],false,lmp);
73 t_period = utils::numeric(FLERR,arg[5],false,lmp);
74 int seed = utils::inumeric(FLERR,arg[6],false,lmp);
75
76 // error checks
77
78 if (t_period <= 0.0) error->all(FLERR,"Illegal fix temp/csvr command");
79 if (seed <= 0) error->all(FLERR,"Illegal fix temp/csvr command");
80
81 random = new RanMars(lmp,seed + comm->me);
82
83 // create a new compute temp style
84 // id = fix-ID + temp, compute group = fix group
85
86 id_temp = utils::strdup(std::string(id) + "_temp");
87 modify->add_compute(fmt::format("{} {} temp",id_temp,group->names[igroup]));
88 tflag = 1;
89
90 nmax = -1;
91 energy = 0.0;
92 }
93
94 /* ---------------------------------------------------------------------- */
95
~FixTempCSVR()96 FixTempCSVR::~FixTempCSVR()
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 nmax = -1;
107 }
108
109 /* ---------------------------------------------------------------------- */
110
setmask()111 int FixTempCSVR::setmask()
112 {
113 int mask = 0;
114 mask |= END_OF_STEP;
115 return mask;
116 }
117
118 /* ---------------------------------------------------------------------- */
119
init()120 void FixTempCSVR::init()
121 {
122 // check variable
123
124 if (tstr) {
125 tvar = input->variable->find(tstr);
126 if (tvar < 0)
127 error->all(FLERR,"Variable name for fix temp/csvr does not exist");
128 if (input->variable->equalstyle(tvar)) tstyle = EQUAL;
129 else error->all(FLERR,"Variable for fix temp/csvr is invalid style");
130 }
131
132 int icompute = modify->find_compute(id_temp);
133 if (icompute < 0)
134 error->all(FLERR,"Temperature ID for fix temp/csvr does not exist");
135 temperature = modify->compute[icompute];
136
137 if (temperature->tempbias) which = BIAS;
138 else which = NOBIAS;
139 }
140
141 /* ---------------------------------------------------------------------- */
142
end_of_step()143 void FixTempCSVR::end_of_step()
144 {
145 // set current t_target
146 // if variable temp, evaluate variable, wrap with clear/add
147
148 double delta = update->ntimestep - update->beginstep;
149
150 if (delta != 0.0) delta /= update->endstep - update->beginstep;
151 if (tstyle == CONSTANT)
152 t_target = t_start + delta * (t_stop-t_start);
153 else {
154 modify->clearstep_compute();
155 t_target = input->variable->compute_equal(tvar);
156 if (t_target < 0.0)
157 error->one(FLERR,
158 "Fix temp/csvr variable returned negative temperature");
159 modify->addstep_compute(update->ntimestep + nevery);
160 }
161
162 const double t_current = temperature->compute_scalar();
163 const double efactor = 0.5 * temperature->dof * force->boltz;
164 const double ekin_old = t_current * efactor;
165 const double ekin_new = t_target * efactor;
166
167 // there is nothing to do, if there are no degrees of freedom
168
169 if (temperature->dof < 1) return;
170
171 // compute velocity scaling factor on root node and broadcast
172
173 double lamda;
174 if (comm->me == 0) {
175 lamda = resamplekin(ekin_old, ekin_new);
176 }
177 MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world);
178
179 double * const * const v = atom->v;
180 const int * const mask = atom->mask;
181 const int nlocal = atom->nlocal;
182
183 if (which == NOBIAS) {
184 for (int i = 0; i < nlocal; i++) {
185 if (mask[i] & groupbit) {
186 v[i][0] *= lamda;
187 v[i][1] *= lamda;
188 v[i][2] *= lamda;
189 }
190 }
191 } else {
192 for (int i = 0; i < nlocal; i++) {
193 if (mask[i] & groupbit) {
194 temperature->remove_bias(i,v[i]);
195 v[i][0] *= lamda;
196 v[i][1] *= lamda;
197 v[i][2] *= lamda;
198 temperature->restore_bias(i,v[i]);
199 }
200 }
201 }
202
203 // tally the kinetic energy transferred between heat bath and system
204
205 energy += ekin_old * (1.0 - lamda*lamda);
206 }
207
208 /* ---------------------------------------------------------------------- */
209
modify_param(int narg,char ** arg)210 int FixTempCSVR::modify_param(int narg, char **arg)
211 {
212 if (strcmp(arg[0],"temp") == 0) {
213 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
214 if (tflag) {
215 modify->delete_compute(id_temp);
216 tflag = 0;
217 }
218 delete [] id_temp;
219 id_temp = utils::strdup(arg[1]);
220
221 int icompute = modify->find_compute(id_temp);
222 if (icompute < 0)
223 error->all(FLERR,"Could not find fix_modify temperature ID");
224 temperature = modify->compute[icompute];
225
226 if (temperature->tempflag == 0)
227 error->all(FLERR,
228 "Fix_modify temperature ID does not compute temperature");
229 if (temperature->igroup != igroup && comm->me == 0)
230 error->warning(FLERR,"Group for fix_modify temp != fix group");
231 return 2;
232 }
233 return 0;
234 }
235
236 /* ---------------------------------------------------------------------- */
237
gamdev(const int ia)238 double FixTempCSVR::gamdev(const int ia)
239 {
240 int j;
241 double am,e,s,v1,v2,x,y;
242
243 if (ia < 1) return 0.0;
244 if (ia < 6) {
245 x=1.0;
246 for (j=1; j<=ia; j++)
247 x *= random->uniform();
248
249 // make certain, that -log() doesn't overflow.
250 if (x < 2.2250759805e-308)
251 x = 708.4;
252 else
253 x = -log(x);
254 } else {
255 restart:
256 do {
257 do {
258 do {
259 v1 = random->uniform();
260 v2 = 2.0*random->uniform() - 1.0;
261 } while (v1*v1 + v2*v2 > 1.0);
262
263 y=v2/v1;
264 am=ia-1;
265 s=sqrt(2.0*am+1.0);
266 x=s*y+am;
267 } while (x <= 0.0);
268
269 if (am*log(x/am)-s*y < -700 || v1<0.00001) {
270 goto restart;
271 }
272
273 e=(1.0+y*y)*exp(am*log(x/am)-s*y);
274 } while (random->uniform() > e);
275 }
276 return x;
277 }
278
279 /* -------------------------------------------------------------------
280 returns the sum of n independent gaussian noises squared
281 (i.e. equivalent to summing the square of the return values of nn
282 calls to gasdev)
283 ---------------------------------------------------------------------- */
284
sumnoises(int nn)285 double FixTempCSVR::sumnoises(int nn) {
286 if (nn == 0) {
287 return 0.0;
288 } else if (nn == 1) {
289 const double rr = random->gaussian();
290 return rr*rr;
291 } else if (nn % 2 == 0) {
292 return 2.0 * gamdev(nn / 2);
293 } else {
294 const double rr = random->gaussian();
295 return 2.0 * gamdev((nn-1) / 2) + rr*rr;
296 }
297 }
298
299 /* -------------------------------------------------------------------
300 returns the scaling factor for velocities to thermalize
301 the system so it samples the canonical ensemble
302 ---------------------------------------------------------------------- */
303
resamplekin(double ekin_old,double ekin_new)304 double FixTempCSVR::resamplekin(double ekin_old, double ekin_new) {
305 const double tdof = temperature->dof;
306 const double c1 = exp(-update->dt/t_period);
307 const double c2 = (1.0-c1)*ekin_new/ekin_old/tdof;
308 const double r1 = random->gaussian();
309 const double r2 = sumnoises(tdof - 1);
310
311 const double scale = c1 + c2*(r1*r1+r2) + 2.0*r1*sqrt(c1*c2);
312 return sqrt(scale);
313 }
314
315 /* ---------------------------------------------------------------------- */
316
reset_target(double t_new)317 void FixTempCSVR::reset_target(double t_new)
318 {
319 t_target = t_start = t_stop = t_new;
320 }
321
322 /* ---------------------------------------------------------------------- */
323
compute_scalar()324 double FixTempCSVR::compute_scalar()
325 {
326 return energy;
327 }
328
329 /* ----------------------------------------------------------------------
330 pack entire state of Fix into one write
331 ------------------------------------------------------------------------- */
332
write_restart(FILE * fp)333 void FixTempCSVR::write_restart(FILE *fp)
334 {
335 int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy
336 double *list = nullptr;
337 if (comm->me == 0) {
338 list = new double[nsize];
339 list[0] = energy;
340 list[1] = comm->nprocs;
341 }
342 double state[PRNGSIZE];
343 random->get_state(state);
344 MPI_Gather(state,PRNGSIZE,MPI_DOUBLE,list+2,PRNGSIZE,MPI_DOUBLE,0,world);
345
346 if (comm->me == 0) {
347 int size = nsize * sizeof(double);
348 fwrite(&size,sizeof(int),1,fp);
349 fwrite(list,sizeof(double),nsize,fp);
350 delete[] list;
351 }
352 }
353
354 /* ----------------------------------------------------------------------
355 use state info from restart file to restart the Fix
356 ------------------------------------------------------------------------- */
357
restart(char * buf)358 void FixTempCSVR::restart(char *buf)
359 {
360 double *list = (double *) buf;
361
362 energy = list[0];
363 int nprocs = (int) list[1];
364 if (nprocs != comm->nprocs) {
365 if (comm->me == 0)
366 error->warning(FLERR,"Different number of procs. Cannot restore RNG state.");
367 } else random->set_state(list+2+comm->me*PRNGSIZE);
368 }
369
370 /* ----------------------------------------------------------------------
371 extract thermostat properties
372 ------------------------------------------------------------------------- */
373
extract(const char * str,int & dim)374 void *FixTempCSVR::extract(const char *str, int &dim)
375 {
376 dim=0;
377 if (strcmp(str,"t_target") == 0) {
378 return &t_target;
379 }
380 return nullptr;
381 }
382