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: Andres Jaramillo-Botero
17 ------------------------------------------------------------------------- */
18
19 #include <cmath>
20
21 #include "fix_langevin_eff.h"
22 #include "atom.h"
23 #include "update.h"
24 #include "modify.h"
25 #include "compute.h"
26 #include "domain.h"
27 #include "input.h"
28 #include "variable.h"
29 #include "random_mars.h"
30 #include "memory.h"
31 #include "error.h"
32 #include "group.h"
33
34 using namespace LAMMPS_NS;
35 using namespace FixConst;
36
37 enum{NOBIAS,BIAS};
38 enum{CONSTANT,EQUAL,ATOM};
39
40 #define SINERTIA 0.4 // moment of inertia prefactor for sphere
41 #define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
42
43 /* ---------------------------------------------------------------------- */
44
FixLangevinEff(LAMMPS * lmp,int narg,char ** arg)45 FixLangevinEff::FixLangevinEff(LAMMPS *lmp, int narg, char **arg) :
46 FixLangevin(lmp, narg, arg)
47 {
48 erforcelangevin = nullptr;
49 }
50
51 /* ---------------------------------------------------------------------- */
52
~FixLangevinEff()53 FixLangevinEff::~FixLangevinEff()
54 {
55 memory->destroy(erforcelangevin);
56 }
57
58 /* ---------------------------------------------------------------------- */
59
post_force(int)60 void FixLangevinEff::post_force(int /*vflag*/)
61 {
62 if (tallyflag) post_force_tally();
63 else post_force_no_tally();
64 }
65
66 /* ---------------------------------------------------------------------- */
67
post_force_no_tally()68 void FixLangevinEff::post_force_no_tally()
69 {
70 double gamma1,gamma2,t_target;
71
72 double **v = atom->v;
73 double **f = atom->f;
74 double *ervel = atom->ervel;
75 double *erforce = atom->erforce;
76 int *spin = atom->spin;
77 int *type = atom->type;
78 int *mask = atom->mask;
79 int nlocal = atom->nlocal;
80 double mefactor = domain->dimension/4.0;
81 double sqrtmefactor = sqrt(mefactor);
82
83 double delta = update->ntimestep - update->beginstep;
84 delta /= update->endstep - update->beginstep;
85
86 // set current t_target and t_sqrt
87 // if variable temp, evaluate variable, wrap with clear/add
88 // reallocate tforce array if necessary
89
90 if (tstyle == CONSTANT) {
91 t_target = t_start + delta * (t_stop-t_start);
92 tsqrt = sqrt(t_target);
93 } else {
94 modify->clearstep_compute();
95 if (tstyle == EQUAL) {
96 t_target = input->variable->compute_equal(tvar);
97 if (t_target < 0.0)
98 error->one(FLERR,"Fix langevin/eff variable returned negative temperature");
99 tsqrt = sqrt(t_target);
100 } else {
101 if (atom->nmax > maxatom2) {
102 maxatom2 = atom->nmax;
103 memory->destroy(tforce);
104 memory->create(tforce,maxatom2,"langevin/eff:tforce");
105 }
106 input->variable->compute_atom(tvar,igroup,tforce,1,0);
107 for (int i = 0; i < nlocal; i++)
108 if (mask[i] & groupbit)
109 if (tforce[i] < 0.0)
110 error->one(FLERR,
111 "Fix langevin/eff variable returned negative temperature");
112 }
113 modify->addstep_compute(update->ntimestep + 1);
114 }
115
116 // apply damping and thermostat to atoms in group
117 // for BIAS:
118 // calculate temperature since some computes require temp
119 // computed on current nlocal atoms to remove bias
120 // test v = 0 since some computes mask non-participating atoms via v = 0
121 // and added force has extra term not multiplied by v = 0
122 // for ZEROFLAG:
123 // sum random force over all atoms in group
124 // subtract sum/particles from each atom in group
125
126 double fran[4],fsum[4],fsumall[4];
127 fsum[0] = fsum[1] = fsum[2] = fsum[3] = 0.0;
128
129 int particles = group->count(igroup);
130 if (zeroflag) {
131 if (particles == 0)
132 error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons");
133 }
134
135 // find number of electrons in group
136 int dof,fix_dof;
137 dof = domain->dimension * particles;
138 fix_dof = 0;
139 for (int i = 0; i < modify->nfix; i++)
140 fix_dof += modify->fix[i]->dof(igroup);
141
142 // extra_dof = domain->dimension
143 dof -= domain->dimension + fix_dof;
144
145 int one = 0;
146 for (int i = 0; i < nlocal; i++)
147 if (mask[i] & groupbit) {
148 if (abs(spin[i])==1) one++;
149 }
150 int nelectrons, dofelectrons, dofnuclei;
151 MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
152 dofelectrons = domain->dimension*nelectrons;
153 dofnuclei = dof-dofelectrons;
154
155 // thermal partitioning factor between nuclei and electrons
156 // extra dof from electron size
157 double gfactor3=(double) (dof+nelectrons)/dofnuclei;
158
159 if (tbiasflag == NOBIAS) {
160 for (int i = 0; i < nlocal; i++) {
161 if (mask[i] & groupbit) {
162 if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
163 gamma1 = gfactor1[type[i]] * gfactor3;
164 gamma2 = gfactor2[type[i]] * tsqrt;
165 fran[0] = gamma2*(random->uniform()-0.5);
166 fran[1] = gamma2*(random->uniform()-0.5);
167 fran[2] = gamma2*(random->uniform()-0.5);
168 f[i][0] += gamma1*v[i][0] + fran[0];
169 f[i][1] += gamma1*v[i][1] + fran[1];
170 f[i][2] += gamma1*v[i][2] + fran[2];
171 fsum[0] += fran[0];
172 fsum[1] += fran[1];
173 fsum[2] += fran[2];
174 if (abs(spin[i])==1) {
175 fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5);
176 erforce[i] += mefactor*gamma1*ervel[i]+fran[3];
177 fsum[3] += fran[3];
178 }
179 }
180 }
181 } else if (tbiasflag == BIAS) {
182 temperature->compute_scalar();
183 for (int i = 0; i < nlocal; i++) {
184 if (mask[i] & groupbit) {
185 if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
186 gamma1 = gfactor1[type[i]] * gfactor3;
187 gamma2 = gfactor2[type[i]] * tsqrt;
188 temperature->remove_bias(i,v[i]);
189 fran[0] = gamma2*(random->uniform()-0.5);
190 fran[1] = gamma2*(random->uniform()-0.5);
191 fran[2] = gamma2*(random->uniform()-0.5);
192 if (v[i][0] != 0.0)
193 f[i][0] += gamma1*v[i][0] + fran[0];
194 if (v[i][1] != 0.0)
195 f[i][1] += gamma1*v[i][1] + fran[1];
196 if (v[i][2] != 0.0)
197 f[i][2] += gamma1*v[i][2] + fran[2];
198 fsum[0] += fran[0];
199 fsum[1] += fran[1];
200 fsum[2] += fran[2];
201 if (abs(spin[i])==1) {
202 fran[3] = sqrtmefactor*gamma2*(random->uniform()-0.5);
203 if (ervel[i] != 0.0) erforce[i] += mefactor*gamma1*ervel[i]+fran[3];
204 fsum[3] += fran[3];
205 }
206 temperature->restore_bias(i,v[i]);
207 }
208 }
209 }
210
211 // set total force to zero
212
213 if (zeroflag) {
214 MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
215 fsumall[0] /= particles;
216 fsumall[1] /= particles;
217 fsumall[2] /= particles;
218 fsumall[3] /= nelectrons;
219 for (int i = 0; i < nlocal; i++) {
220 if (mask[i] & groupbit) {
221 f[i][0] -= fsumall[0];
222 f[i][1] -= fsumall[1];
223 f[i][2] -= fsumall[2];
224 if (abs(spin[i])==1) erforce[i] -= fsumall[3];
225 }
226 }
227 }
228 }
229
230 /* ---------------------------------------------------------------------- */
231
post_force_tally()232 void FixLangevinEff::post_force_tally()
233 {
234 double gamma1,gamma2,t_target;
235
236 // reallocate flangevin and erforcelangevin if necessary
237
238 if (atom->nmax > maxatom1) {
239 memory->destroy(flangevin);
240 memory->destroy(erforcelangevin);
241 maxatom1 = atom->nmax;
242 memory->create(flangevin,maxatom1,3,"langevin/eff:flangevin");
243 memory->create(erforcelangevin,maxatom1,"langevin/eff:erforcelangevin");
244 }
245
246 double **v = atom->v;
247 double **f = atom->f;
248 double *erforce = atom->erforce;
249 double *ervel = atom->ervel;
250 int *spin = atom->spin;
251 double mefactor = domain->dimension/4.0;
252 double sqrtmefactor = sqrt(mefactor);
253
254 int *type = atom->type;
255 int *mask = atom->mask;
256 int nlocal = atom->nlocal;
257
258 double delta = update->ntimestep - update->beginstep;
259 delta /= update->endstep - update->beginstep;
260
261 // set current t_target and t_sqrt
262 // if variable temp, evaluate variable, wrap with clear/add
263 // reallocate tforce array if necessary
264
265 if (tstyle == CONSTANT) {
266 t_target = t_start + delta * (t_stop-t_start);
267 tsqrt = sqrt(t_target);
268 } else {
269 modify->clearstep_compute();
270 if (tstyle == EQUAL) {
271 t_target = input->variable->compute_equal(tvar);
272 if (t_target < 0.0)
273 error->one(FLERR,"Fix langevin/eff variable returned negative temperature");
274 tsqrt = sqrt(t_target);
275 } else {
276 if (atom->nmax > maxatom2) {
277 maxatom2 = atom->nmax;
278 memory->destroy(tforce);
279 memory->create(tforce,maxatom2,"langevin/eff:tforce");
280 }
281 input->variable->compute_atom(tvar,igroup,tforce,1,0);
282 for (int i = 0; i < nlocal; i++)
283 if (mask[i] & groupbit)
284 if (tforce[i] < 0.0)
285 error->one(FLERR,
286 "Fix langevin/eff variable returned negative temperature");
287 }
288 modify->addstep_compute(update->ntimestep + 1);
289 }
290
291 // apply damping and thermostat to appropriate atoms
292 // for BIAS:
293 // calculate temperature since some computes require temp
294 // computed on current nlocal atoms to remove bias
295 // test v = 0 since some computes mask non-participating atoms via v = 0
296 // and added force has extra term not multiplied by v = 0
297
298 int particles = group->count(igroup);
299 if (zeroflag) {
300 if (particles == 0)
301 error->all(FLERR,"Cannot zero Langevin force of 0 atoms/electrons");
302 }
303
304 // find number of electrons in group
305 int dof,fix_dof;
306 dof = domain->dimension * particles;
307 fix_dof = 0;
308 for (int i = 0; i < modify->nfix; i++)
309 fix_dof += modify->fix[i]->dof(igroup);
310
311 // extra_dof = domain->dimension
312 dof -= domain->dimension + fix_dof;
313
314 int one = 0;
315 for (int i = 0; i < nlocal; i++)
316 if (mask[i] & groupbit) {
317 if (abs(spin[i])==1) one++;
318 }
319 int nelectrons, dofelectrons, dofnuclei;
320 MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world);
321 dofelectrons = domain->dimension*nelectrons;
322 dofnuclei = dof-dofelectrons;
323
324 // thermal partitioning factor between nuclei and electrons
325 // extra dof from electron size
326 double gfactor3=(double) (dof+nelectrons)/dofnuclei;
327
328 if (tbiasflag == NOBIAS) {
329 for (int i = 0; i < nlocal; i++) {
330 if (mask[i] & groupbit) {
331 if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
332 gamma1 = gfactor1[type[i]] * gfactor3;
333 gamma2 = gfactor2[type[i]] * tsqrt;
334 flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
335 flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
336 flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
337 f[i][0] += flangevin[i][0];
338 f[i][1] += flangevin[i][1];
339 f[i][2] += flangevin[i][2];
340 if (abs(spin[i])==1) {
341 erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5);
342 erforce[i] += erforcelangevin[i];
343 }
344 }
345 }
346 } else if (tbiasflag == BIAS) {
347 temperature->compute_scalar();
348 for (int i = 0; i < nlocal; i++) {
349 if (mask[i] & groupbit) {
350 if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
351 gamma1 = gfactor1[type[i]] * gfactor3;
352 gamma2 = gfactor2[type[i]] * tsqrt;
353 temperature->remove_bias(i,v[i]);
354 flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
355 flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
356 flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
357 if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
358 else flangevin[i][0] = 0.0;
359 if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
360 else flangevin[i][1] = 0.0;
361 if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
362 else flangevin[i][2] = 0.0;
363 if (abs(spin[i])==1) {
364 erforcelangevin[i] = mefactor*gamma1*ervel[i]+sqrtmefactor*gamma2*(random->uniform()-0.5);
365 if (ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
366 else erforcelangevin[i] = 0.0;
367 }
368 temperature->restore_bias(i,v[i]);
369 }
370 }
371 }
372 }
373
374 /* ----------------------------------------------------------------------
375 tally energy transfer to thermal reservoir
376 ------------------------------------------------------------------------- */
377
end_of_step()378 void FixLangevinEff::end_of_step()
379 {
380 if (!tallyflag) return;
381
382 double **v = atom->v;
383 int *mask = atom->mask;
384 int nlocal = atom->nlocal;
385 int *spin = atom->spin;
386
387 energy_onestep = 0.0;
388
389 for (int i = 0; i < nlocal; i++)
390 if (mask[i] & groupbit) {
391 energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
392 flangevin[i][2]*v[i][2];
393 if (abs(spin[i])==1) energy_onestep += erforcelangevin[i];
394 }
395 energy += energy_onestep*update->dt;
396 }
397
398 /* ---------------------------------------------------------------------- */
399
compute_scalar()400 double FixLangevinEff::compute_scalar()
401 {
402 if (!tallyflag || flangevin == nullptr || erforcelangevin == nullptr) return 0.0;
403
404 // capture the very first energy transfer to thermal reservoir
405
406 double **v = atom->v;
407 int *mask = atom->mask;
408 int nlocal = atom->nlocal;
409 int *spin = atom->spin;
410
411 if (update->ntimestep == update->beginstep) {
412 energy_onestep = 0.0;
413 for (int i = 0; i < nlocal; i++)
414 if (mask[i] & groupbit) {
415 energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
416 flangevin[i][2]*v[i][2];
417 if (abs(spin[i])==1) energy_onestep += erforcelangevin[i];
418 }
419 energy = 0.5*energy_onestep*update->dt;
420 }
421
422 double energy_me = energy - 0.5*energy_onestep*update->dt;
423
424 double energy_all;
425 MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
426 return -energy_all;
427 }
428