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 authors: (in addition to authors of original fix ttm)
17 Sergey Starikov (Joint Institute for High Temperatures of RAS)
18 Vasily Pisarev (Joint Institute for High Temperatures of RAS)
19 ------------------------------------------------------------------------- */
20
21 #include "fix_ttm_mod.h"
22
23 #include "atom.h"
24 #include "citeme.h"
25 #include "comm.h"
26 #include "domain.h"
27 #include "error.h"
28 #include "force.h"
29 #include "math_const.h"
30 #include "memory.h"
31 #include "random_mars.h"
32 #include "respa.h"
33 #include "potential_file_reader.h"
34 #include "tokenizer.h"
35 #include "update.h"
36
37 #include <cmath>
38 #include <cstring>
39
40 using namespace LAMMPS_NS;
41 using namespace FixConst;
42 using namespace MathConst;
43
44 // OFFSET avoids outside-of-box atoms being rounded to grid pts incorrectly
45 // SHIFT = 0.0 assigns atoms to lower-left grid pt
46 // SHIFT = 0.5 assigns atoms to nearest grid pt
47 // use SHIFT = 0.0 for now since it allows fix ave/chunk
48 // to spatially average consistent with the TTM grid
49
50 static const char cite_fix_ttm_mod[] =
51 "fix ttm/mod command:\n\n"
52 "@article{Pisarev2014,\n"
53 "author = {Pisarev, V. V. and Starikov, S. V.},\n"
54 "title = {{Atomistic simulation of ion track formation in UO2.}},\n"
55 "journal = {J.~Phys.:~Condens.~Matter},\n"
56 "volume = {26},\n"
57 "number = {47},\n"
58 "pages = {475401},\n"
59 "year = {2014}\n"
60 "}\n\n"
61 "@article{Norman2013,\n"
62 "author = {Norman, G. E. and Starikov, S. V. and Stegailov, V. V. and Saitov, I. M. and Zhilyaev, P. A.},\n"
63 "title = {{Atomistic Modeling of Warm Dense Matter in the Two-Temperature State}},\n"
64 "journal = {Contrib.~Plasm.~Phys.},\n"
65 "number = {2},\n"
66 "volume = {53},\n"
67 "pages = {129--139},\n"
68 "year = {2013}\n"
69 "}\n\n";
70
71 static constexpr int OFFSET = 16384;
72 static constexpr double SHIFT = 0.0;
73
74 /* ---------------------------------------------------------------------- */
75
FixTTMMod(LAMMPS * lmp,int narg,char ** arg)76 FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) :
77 Fix(lmp, narg, arg),
78 random(nullptr), nsum(nullptr), nsum_all(nullptr),
79 gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), flangevin(nullptr),
80 T_electron(nullptr), T_electron_old(nullptr), sum_vsq(nullptr), sum_mass_vsq(nullptr),
81 sum_vsq_all(nullptr), sum_mass_vsq_all(nullptr), net_energy_transfer(nullptr),
82 net_energy_transfer_all(nullptr)
83 {
84 if (lmp->citeme) lmp->citeme->add(cite_fix_ttm_mod);
85
86 if (narg < 8) error->all(FLERR,"Illegal fix ttm/mod command");
87
88 vector_flag = 1;
89 size_vector = 2;
90 global_freq = 1;
91 extvector = 1;
92 nevery = 1;
93 restart_peratom = 1;
94 restart_global = 1;
95
96 seed = utils::inumeric(FLERR,arg[3],false,lmp);
97
98 nxgrid = utils::inumeric(FLERR,arg[5],false,lmp);
99 nygrid = utils::inumeric(FLERR,arg[6],false,lmp);
100 nzgrid = utils::inumeric(FLERR,arg[7],false,lmp);
101
102 double tinit = 0.0;
103 infile = outfile = nullptr;
104
105 int iarg = 8;
106 while (iarg < narg) {
107 if (strcmp(arg[iarg],"set") == 0) {
108 if (iarg+2 > narg) error->all(FLERR,"Illegal fix ttm/mod command");
109 tinit = utils::numeric(FLERR,arg[iarg+1],false,lmp);
110 if (tinit <= 0.0)
111 error->all(FLERR,"Fix ttm/mod initial temperature must be > 0.0");
112 iarg += 2;
113 } else if (strcmp(arg[iarg],"infile") == 0) {
114 if (iarg+2 > narg) error->all(FLERR,"Illegal fix ttm/mod command");
115 infile = utils::strdup(arg[iarg+1]);
116 iarg += 2;
117 } else if (strcmp(arg[iarg],"outfile") == 0) {
118 if (iarg+3 > narg) error->all(FLERR,"Illegal fix ttm/mod command");
119 outevery = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
120 outfile = utils::strdup(arg[iarg+2]);
121 iarg += 3;
122 } else error->all(FLERR,"Illegal fix ttm/mod command");
123 }
124
125 // error check
126
127 if (seed <= 0)
128 error->all(FLERR,"Invalid random number seed in fix ttm/mod command");
129 if (nxgrid <= 0 || nygrid <= 0 || nzgrid <= 0)
130 error->all(FLERR,"Fix ttm/mod grid sizes must be > 0");
131
132 // check for allowed maximum number of total grid points
133
134 bigint total_ngrid = (bigint) nxgrid * nygrid * nzgrid;
135 if (total_ngrid > MAXSMALLINT)
136 error->all(FLERR,"Too many grid points in fix ttm/mod");
137 ngridtotal = total_ngrid;
138
139 // t_surface is determined by electronic temperature (not constant)
140
141 read_parameters(arg[4]);
142
143 t_surface_l = surface_l;
144 mult_factor = intensity;
145 duration = 0.0;
146 v_0_sq = v_0*v_0;
147 surface_double = double(t_surface_l)*(domain->xprd/nxgrid);
148 if ((C_limit+esheat_0) < 0.0)
149 error->all(FLERR,"Fix ttm/mod electronic_specific_heat must be >= 0.0");
150 if (electronic_density <= 0.0)
151 error->all(FLERR,"Fix ttm/mod electronic_density must be > 0.0");
152 if (gamma_p < 0.0) error->all(FLERR,"Fix ttm/mod gamma_p must be >= 0.0");
153 if (gamma_s < 0.0) error->all(FLERR,"Fix ttm/mod gamma_s must be >= 0.0");
154 if (v_0 < 0.0) error->all(FLERR,"Fix ttm/mod v_0 must be >= 0.0");
155 if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm/mod ionic_density must be > 0.0");
156 if (surface_l < 0) error->all(FLERR,"Surface coordinates must be >= 0");
157 if (surface_l >= surface_r) error->all(FLERR, "Left surface coordinate must be less than right surface coordinate");
158
159 // initialize Marsaglia RNG with processor-unique seed
160
161 random = new RanMars(lmp,seed + comm->me);
162
163 // allocate per-type arrays for force prefactors
164
165 gfactor1 = new double[atom->ntypes+1];
166 gfactor2 = new double[atom->ntypes+1];
167
168 // allocate 3d grid variables
169
170 memory->create(nsum,nxgrid,nygrid,nzgrid,"ttm/mod:nsum");
171 memory->create(nsum_all,nxgrid,nygrid,nzgrid,"ttm/mod:nsum_all");
172 memory->create(sum_vsq,nxgrid,nygrid,nzgrid,"ttm/mod:sum_vsq");
173 memory->create(sum_mass_vsq,nxgrid,nygrid,nzgrid,"ttm/mod:sum_mass_vsq");
174 memory->create(sum_vsq_all,nxgrid,nygrid,nzgrid,"ttm/mod:sum_vsq_all");
175 memory->create(sum_mass_vsq_all,nxgrid,nygrid,nzgrid,
176 "ttm/mod:sum_mass_vsq_all");
177 memory->create(T_electron_old,nxgrid,nygrid,nzgrid,"ttm/mod:T_electron_old");
178 memory->create(T_electron_first,nxgrid,nygrid,nzgrid,"ttm/mod:T_electron_first");
179 memory->create(T_electron,nxgrid,nygrid,nzgrid,"ttm/mod:T_electron");
180 memory->create(net_energy_transfer,nxgrid,nygrid,nzgrid,
181 "ttm/mod:net_energy_transfer");
182 memory->create(net_energy_transfer_all,nxgrid,nygrid,nzgrid,
183 "ttm/mod:net_energy_transfer_all");
184 flangevin = nullptr;
185 grow_arrays(atom->nmax);
186
187 // grid OFFSET to perform
188 // SHIFT to map atom to nearest or lower-left grid point
189
190 shift = OFFSET + SHIFT;
191
192 // zero out the flangevin array
193
194 for (int i = 0; i < atom->nmax; i++) {
195 flangevin[i][0] = 0.0;
196 flangevin[i][1] = 0.0;
197 flangevin[i][2] = 0.0;
198 }
199
200 atom->add_callback(Atom::GROW);
201 atom->add_callback(Atom::RESTART);
202
203 // initialize electron temperatures on grid
204
205 int ix,iy,iz;
206 for (ix = 0; ix < nxgrid; ix++)
207 for (iy = 0; iy < nygrid; iy++)
208 for (iz = 0; iz < nzgrid; iz++)
209 T_electron[ix][iy][iz] = tinit;
210
211 // if specified, read initial electron temperatures from file
212
213 if (infile) read_electron_temperatures(infile);
214 }
215
216 /* ---------------------------------------------------------------------- */
217
~FixTTMMod()218 FixTTMMod::~FixTTMMod()
219 {
220 delete random;
221 delete[] gfactor1;
222 delete[] gfactor2;
223
224 memory->destroy(nsum);
225 memory->destroy(nsum_all);
226 memory->destroy(sum_vsq);
227 memory->destroy(sum_mass_vsq);
228 memory->destroy(sum_vsq_all);
229 memory->destroy(sum_mass_vsq_all);
230 memory->destroy(T_electron_first);
231 memory->destroy(T_electron_old);
232 memory->destroy(T_electron);
233 memory->destroy(flangevin);
234 memory->destroy(net_energy_transfer);
235 memory->destroy(net_energy_transfer_all);
236 }
237
238 /* ---------------------------------------------------------------------- */
239
setmask()240 int FixTTMMod::setmask()
241 {
242 int mask = 0;
243 mask |= POST_FORCE;
244 mask |= POST_FORCE_RESPA;
245 mask |= END_OF_STEP;
246 return mask;
247 }
248
249 /* ---------------------------------------------------------------------- */
250
init()251 void FixTTMMod::init()
252 {
253 if (domain->dimension == 2)
254 error->all(FLERR,"Cannot use fix ttm/mod with 2d simulation");
255 if (domain->nonperiodic != 0)
256 error->all(FLERR,"Cannot use non-periodic boundares with fix ttm/mod");
257 if (domain->triclinic)
258 error->all(FLERR,"Cannot use fix ttm/mod with triclinic box");
259
260 // set force prefactors
261
262 for (int i = 1; i <= atom->ntypes; i++) {
263 gfactor1[i] = - gamma_p / force->ftm2v;
264 gfactor2[i] =
265 sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
266 }
267
268 for (int ix = 0; ix < nxgrid; ix++)
269 for (int iy = 0; iy < nygrid; iy++)
270 for (int iz = 0; iz < nzgrid; iz++)
271 net_energy_transfer_all[ix][iy][iz] = 0;
272
273 if (utils::strmatch(update->integrate_style,"^respa"))
274 nlevels_respa = ((Respa *) update->integrate)->nlevels;
275 }
276
277 /* ---------------------------------------------------------------------- */
278
setup(int vflag)279 void FixTTMMod::setup(int vflag)
280 {
281 if (utils::strmatch(update->integrate_style,"^verlet")) {
282 post_force_setup(vflag);
283 } else {
284 ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
285 post_force_respa_setup(vflag,nlevels_respa-1,0);
286 ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
287 }
288 }
289
290 /* ---------------------------------------------------------------------- */
291
post_force(int)292 void FixTTMMod::post_force(int /*vflag*/)
293 {
294 double **x = atom->x;
295 double **v = atom->v;
296 double **f = atom->f;
297 int *type = atom->type;
298 int *mask = atom->mask;
299 int nlocal = atom->nlocal;
300
301 double dx = domain->xprd/nxgrid;
302 double dy = domain->yprd/nygrid;
303 double dz = domain->zprd/nzgrid;
304 double gamma1,gamma2;
305
306 // apply damping and thermostat to all atoms in fix group
307
308 for (int i = 0; i < nlocal; i++) {
309 if (mask[i] & groupbit) {
310 double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
311 double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
312 double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
313 int ix = static_cast<int>(xscale*nxgrid + shift) - OFFSET;
314 int iy = static_cast<int>(yscale*nygrid + shift) - OFFSET;
315 int iz = static_cast<int>(zscale*nzgrid + shift) - OFFSET;
316 while (ix > nxgrid-1) ix -= nxgrid;
317 while (iy > nygrid-1) iy -= nygrid;
318 while (iz > nzgrid-1) iz -= nzgrid;
319 while (ix < 0) ix += nxgrid;
320 while (iy < 0) iy += nygrid;
321 while (iz < 0) iz += nzgrid;
322
323 if (T_electron[ix][iy][iz] < 0)
324 error->all(FLERR,"Electronic temperature dropped below zero");
325
326 double tsqrt = sqrt(T_electron[ix][iy][iz]);
327
328 gamma1 = gfactor1[type[i]];
329 double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
330 if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
331 gamma2 = gfactor2[type[i]] * tsqrt;
332 if (ix >= surface_l) {
333 if (ix < surface_r) {
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 double x_surf = dx*double(surface_l)+dx;
338 double x_at = x[i][0] - domain->boxlo[0];
339 int right_x = ix + 1;
340 int right_y = iy + 1;
341 int right_z = iz + 1;
342 if (right_x == nxgrid) right_x = 0;
343 if (right_y == nygrid) right_y = 0;
344 if (right_z == nzgrid) right_z = 0;
345 int left_x = ix - 1;
346 int left_y = iy - 1;
347 int left_z = iz - 1;
348 if (left_x == -1) left_x = nxgrid - 1;
349 if (left_y == -1) left_y = nygrid - 1;
350 if (left_z == -1) left_z = nzgrid - 1;
351 double T_i = T_electron[ix][iy][iz];
352 double T_ir = T_electron[right_x][iy][iz];
353 double T_iu = T_electron[ix][right_y][iz];
354 double T_if = T_electron[ix][iy][right_z];
355 double C_i = el_properties(T_electron[ix][iy][iz]).el_heat_capacity;
356 double C_ir = el_properties(T_electron[right_x][iy][iz]).el_heat_capacity;
357 double C_iu = el_properties(T_electron[ix][right_y][iz]).el_heat_capacity;
358 double C_if = el_properties(T_electron[ix][iy][right_z]).el_heat_capacity;
359 double diff_x = (x_at - x_surf)*(x_at - x_surf);
360 diff_x = pow(diff_x,0.5);
361 double len_factor = diff_x/(diff_x+free_path);
362 if (movsur == 1) {
363 if (x_at >= x_surf) {
364 flangevin[i][0] -= pres_factor/ionic_density*((C_ir*T_ir*free_path/(diff_x+free_path)/(diff_x+free_path)) +
365 (len_factor/dx)*(C_ir*T_ir-C_i*T_i));
366 flangevin[i][1] -= pres_factor/ionic_density/dy*(C_iu*T_iu-C_i*T_i);
367 flangevin[i][2] -= pres_factor/ionic_density/dz*(C_if*T_if-C_i*T_i);
368 }
369 } else {
370 flangevin[i][0] -= pres_factor/ionic_density/dx*(C_ir*T_ir-C_i*T_i);
371 flangevin[i][1] -= pres_factor/ionic_density/dy*(C_iu*T_iu-C_i*T_i);
372 flangevin[i][2] -= pres_factor/ionic_density/dz*(C_if*T_if-C_i*T_i);
373 }
374 f[i][0] += flangevin[i][0];
375 f[i][1] += flangevin[i][1];
376 f[i][2] += flangevin[i][2];
377 }
378 }
379 if (movsur == 1) {
380 if (ix < surface_l) {
381 t_surface_l = ix;
382 }
383 }
384 }
385 }
386 MPI_Allreduce(&t_surface_l,&surface_l,1,MPI_INT,MPI_MIN,world);
387 }
388
389 /* ---------------------------------------------------------------------- */
390
post_force_setup(int)391 void FixTTMMod::post_force_setup(int /*vflag*/)
392 {
393 double **f = atom->f;
394 int *mask = atom->mask;
395 int nlocal = atom->nlocal;
396
397 // apply langevin forces that have been stored from previous run
398
399 for (int i = 0; i < nlocal; i++) {
400 if (mask[i] & groupbit) {
401 f[i][0] += flangevin[i][0];
402 f[i][1] += flangevin[i][1];
403 f[i][2] += flangevin[i][2];
404 }
405 }
406 }
407
408 /* ---------------------------------------------------------------------- */
409
post_force_respa(int vflag,int ilevel,int)410 void FixTTMMod::post_force_respa(int vflag, int ilevel, int /*iloop*/)
411 {
412 if (ilevel == nlevels_respa-1) post_force(vflag);
413 }
414
415 /* ---------------------------------------------------------------------- */
416
post_force_respa_setup(int vflag,int ilevel,int)417 void FixTTMMod::post_force_respa_setup(int vflag, int ilevel, int /*iloop*/)
418 {
419 if (ilevel == nlevels_respa-1) post_force_setup(vflag);
420 }
421
422 /* ---------------------------------------------------------------------- */
423
reset_dt()424 void FixTTMMod::reset_dt()
425 {
426 for (int i = 1; i <= atom->ntypes; i++)
427 gfactor2[i] =
428 sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
429 }
430
431 /* ----------------------------------------------------------------------
432 read in ttm/mod parameters from a user-specified file
433 only called by proc 0
434 ------------------------------------------------------------------------- */
435
read_parameters(const std::string & filename)436 void FixTTMMod::read_parameters(const std::string &filename)
437 {
438 try {
439 PotentialFileReader reader(lmp, filename, "ttm/mod parameter");
440
441 // C0 (metal)
442
443 reader.next_line();
444 esheat_0 = reader.next_values(1).next_double();
445
446 // C1 (metal*10^3)
447
448 reader.next_line();
449 esheat_1 = reader.next_values(1).next_double();
450
451 // C2 (metal*10^6)
452
453 reader.next_line();
454 esheat_2 = reader.next_values(1).next_double();
455
456 // C3 (metal*10^9)
457
458 reader.next_line();
459 esheat_3 = reader.next_values(1).next_double();
460
461 // C4 (metal*10^12)
462
463 reader.next_line();
464 esheat_4 = reader.next_values(1).next_double();
465
466 // C_limit
467
468 reader.next_line();
469 C_limit = reader.next_values(1).next_double();
470
471 // Temperature damping factor
472
473 reader.next_line();
474 T_damp = reader.next_values(1).next_double();
475
476 // rho_e
477
478 reader.next_line();
479 electronic_density = reader.next_values(1).next_double();
480
481 // thermal_diffusion
482
483 reader.next_line();
484 el_th_diff = reader.next_values(1).next_double();
485
486 // gamma_p
487
488 reader.next_line();
489 gamma_p = reader.next_values(1).next_double();
490
491 // gamma_s
492
493 reader.next_line();
494 gamma_s = reader.next_values(1).next_double();
495
496 // v0
497
498 reader.next_line();
499 v_0 = reader.next_values(1).next_double();
500
501 // average intensity of pulse (source of energy) (metal units)
502
503 reader.next_line();
504 intensity = reader.next_values(1).next_double();
505
506 // coordinate of 1st surface in x-direction (in box units) - constant
507
508 reader.next_line();
509 surface_l = reader.next_values(1).next_int();
510
511 // coordinate of 2nd surface in x-direction (in box units) - constant
512
513 reader.next_line();
514 surface_r = reader.next_values(1).next_int();
515
516 // skin_layer = intensity is reduced (I=I0*exp[-x/skin_layer])
517
518 reader.next_line();
519 skin_layer = reader.next_values(1).next_int();
520
521 // width of pulse (picoseconds)
522
523 reader.next_line();
524 width = reader.next_values(1).next_double();
525
526 // factor of electronic pressure (PF) Pe = PF*Ce*Te
527
528 reader.next_line();
529 pres_factor = reader.next_values(1).next_double();
530
531 // effective free path of electrons (angstrom)
532
533 reader.next_line();
534 free_path = reader.next_values(1).next_double();
535
536 // ionic density (ions*angstrom^{-3})
537
538 reader.next_line();
539 ionic_density = reader.next_values(1).next_double();
540
541 // if movsur = 0: surface is frozen
542
543 reader.next_line();
544 movsur = reader.next_values(1).next_int();
545
546 // electron_temperature_min
547
548 reader.next_line();
549 electron_temperature_min = reader.next_values(1).next_double();
550 } catch (std::exception &e) {
551 error->one(FLERR,e.what());
552 }
553 }
554
555 /* ----------------------------------------------------------------------
556 read in initial electron temperatures from a user-specified file
557 only read by proc 0, grid values are Bcast to other procs
558 ------------------------------------------------------------------------- */
559
read_electron_temperatures(const std::string & filename)560 void FixTTMMod::read_electron_temperatures(const std::string &filename)
561 {
562 if (comm->me == 0) {
563
564 int ***T_initial_set;
565 memory->create(T_initial_set,nxgrid,nygrid,nzgrid,"ttm/mod:T_initial_set");
566 memset(&T_initial_set[0][0][0],0,ngridtotal*sizeof(int));
567
568 // read initial electron temperature values from file
569 bigint nread = 0;
570
571 try {
572 PotentialFileReader reader(lmp, filename, "electron temperature grid");
573
574 while (nread < ngridtotal) {
575 // reader will skip over comment-only lines
576 auto values = reader.next_values(4);
577 ++nread;
578
579 int ix = values.next_int();
580 int iy = values.next_int();
581 int iz = values.next_int();
582 double T_tmp = values.next_double();
583
584 // check correctness of input data
585
586 if ((ix < 0) || (ix >= nxgrid) || (iy < 0) || (iy >= nygrid) || (iz < 0) || (iz >= nzgrid))
587 throw parser_error("Fix ttm invalid grid index in fix ttm/mod grid file");
588
589 if (T_tmp < 0.0)
590 throw parser_error("Fix ttm electron temperatures must be > 0.0");
591
592 T_electron[iz][iy][ix] = T_tmp;
593 T_initial_set[iz][iy][ix] = 1;
594 }
595 } catch (std::exception &e) {
596 error->one(FLERR, e.what());
597 }
598
599 // check completeness of input data
600
601 for (int iz = 0; iz < nzgrid; iz++)
602 for (int iy = 0; iy < nygrid; iy++)
603 for (int ix = 0; ix < nxgrid; ix++)
604 if (T_initial_set[iz][iy][ix] == 0)
605 error->all(FLERR,"Fix ttm/mod infile did not set all temperatures");
606
607 memory->destroy(T_initial_set);
608 }
609
610 MPI_Bcast(&T_electron[0][0][0],ngridtotal,MPI_DOUBLE,0,world);
611 }
612
613 /* ----------------------------------------------------------------------
614 write out current electron temperatures to user-specified file
615 only written by proc 0
616 ------------------------------------------------------------------------- */
617
write_electron_temperatures(const std::string & filename)618 void FixTTMMod::write_electron_temperatures(const std::string &filename)
619 {
620 if (comm->me) return;
621
622 FILE *fp = fopen(filename.c_str(),"w");
623 if (!fp) error->one(FLERR,"Fix ttm/mod could not open output file {}: {}",
624 filename, utils::getsyserror());
625 fmt::print(fp,"# DATE: {} UNITS: {} COMMENT: Electron temperature "
626 "{}x{}x{} grid at step {}. Created by fix {}\n", utils::current_date(),
627 update->unit_style, nxgrid, nygrid, nzgrid, update->ntimestep, style);
628
629 int ix,iy,iz;
630
631 for (ix = 0; ix < nxgrid; ix++)
632 for (iy = 0; iy < nygrid; iy++)
633 for (iz = 0; iz < nzgrid; iz++) {
634 if (movsur == 1 && T_electron[ix][iy][iz] == 0.0)
635 T_electron[ix][iy][iz] = electron_temperature_min;
636 fprintf(fp,"%d %d %d %20.16g\n",ix,iy,iz,T_electron[ix][iy][iz]);
637 }
638
639 fclose(fp);
640 }
641
642 /* ---------------------------------------------------------------------- */
643
el_properties(double T_e)644 el_heat_capacity_thermal_conductivity FixTTMMod::el_properties(double T_e)
645 {
646 el_heat_capacity_thermal_conductivity properties;
647 double T_temp = T_e/1000.0, T_reduced = T_damp*T_temp;
648 double T2 = T_temp*T_temp;
649 double T3 = T2*T_temp;
650 double T4 = T3*T_temp;
651 double poly = esheat_0 + esheat_1*T_temp + esheat_2*T2 + esheat_3*T3 + esheat_4*T4;
652 properties.el_heat_capacity = electronic_density*(poly*exp(-T_reduced*T_reduced) + C_limit); // heat capacity
653 properties.el_thermal_conductivity = el_th_diff*properties.el_heat_capacity; // thermal conductivity
654 return properties;
655 }
el_sp_heat_integral(double T_e)656 double FixTTMMod::el_sp_heat_integral(double T_e)
657 {
658 double T_temp = T_e/1000.0, T_reduced = T_damp*T_temp;
659 if (T_damp != 0)
660 return electronic_density*(MY_PIS*(3*esheat_4/pow(T_damp,5)+2*esheat_2/pow(T_damp,3)+4*esheat_0/T_damp)*erf(T_reduced)+
661 4*esheat_3/pow(T_damp,4)+4*esheat_1/T_damp/T_damp-
662 ((6*esheat_4*T_temp+4*esheat_3)/pow(T_damp,4)+
663 (4*esheat_1+4*esheat_4*pow(T_temp,3)+4*esheat_3*T_temp*T_temp+4*esheat_2*T_temp)/T_damp/T_damp)*exp(-T_reduced*T_reduced))*125.0+electronic_density*C_limit*T_e;
664 else
665 return electronic_density*((esheat_0 + C_limit)*T_e + esheat_1*T_temp*T_e/2.0 + esheat_2*T_temp*T_temp*T_e/3.0 + esheat_3*pow(T_temp,3)*T_e/4.0 + esheat_4*pow(T_temp,4)*T_e/5.0);
666 }
667
668 /* ---------------------------------------------------------------------- */
669
end_of_step()670 void FixTTMMod::end_of_step()
671 {
672 double **x = atom->x;
673 double **v = atom->v;
674 int *mask = atom->mask;
675 int nlocal = atom->nlocal;
676
677 if (movsur == 1) {
678 for (int ix = 0; ix < nxgrid; ix++)
679 for (int iy = 0; iy < nygrid; iy++)
680 for (int iz = 0; iz < nzgrid; iz++) {
681 double TTT = T_electron[ix][iy][iz];
682 if (TTT > 0) {
683 if (ix < t_surface_l)
684 t_surface_l = ix;
685 }
686 }
687 }
688 for (int ix = 0; ix < nxgrid; ix++)
689 for (int iy = 0; iy < nygrid; iy++)
690 for (int iz = 0; iz < nzgrid; iz++)
691 net_energy_transfer[ix][iy][iz] = 0;
692
693 for (int i = 0; i < nlocal; i++)
694 if (mask[i] & groupbit) {
695 double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
696 double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
697 double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
698 int ix = static_cast<int>(xscale*nxgrid + shift) - OFFSET;
699 int iy = static_cast<int>(yscale*nygrid + shift) - OFFSET;
700 int iz = static_cast<int>(zscale*nzgrid + shift) - OFFSET;
701 while (ix > nxgrid-1) ix -= nxgrid;
702 while (iy > nygrid-1) iy -= nygrid;
703 while (iz > nzgrid-1) iz -= nzgrid;
704 while (ix < 0) ix += nxgrid;
705 while (iy < 0) iy += nygrid;
706 while (iz < 0) iz += nzgrid;
707 if (ix >= t_surface_l) {
708 if (ix < surface_r)
709 net_energy_transfer[ix][iy][iz] +=
710 (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
711 flangevin[i][2]*v[i][2]);
712 }
713 }
714
715 MPI_Allreduce(&net_energy_transfer[0][0][0],
716 &net_energy_transfer_all[0][0][0],
717 ngridtotal,MPI_DOUBLE,MPI_SUM,world);
718
719 double dx = domain->xprd/nxgrid;
720 double dy = domain->yprd/nygrid;
721 double dz = domain->zprd/nzgrid;
722 double del_vol = dx*dy*dz;
723 double el_specific_heat = 0.0;
724 double el_thermal_conductivity = el_properties(electron_temperature_min).el_thermal_conductivity;
725 for (int ix = 0; ix < nxgrid; ix++)
726 for (int iy = 0; iy < nygrid; iy++)
727 for (int iz = 0; iz < nzgrid; iz++)
728 {
729 if (el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity > el_thermal_conductivity)
730 el_thermal_conductivity = el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity;
731 if (el_specific_heat > 0.0)
732 {
733 if ((T_electron[ix][iy][iz] > 0.0) && (el_properties(T_electron[ix][iy][iz]).el_heat_capacity < el_specific_heat))
734 el_specific_heat = el_properties(T_electron[ix][iy][iz]).el_heat_capacity;
735 }
736 else if (T_electron[ix][iy][iz] > 0.0) el_specific_heat = el_properties(T_electron[ix][iy][iz]).el_heat_capacity;
737 }
738 // num_inner_timesteps = # of inner steps (thermal solves)
739 // required this MD step to maintain a stable explicit solve
740
741 int num_inner_timesteps = 1;
742 double inner_dt = update->dt;
743 double stability_criterion = 0.0;
744
745 for (int ix = 0; ix < nxgrid; ix++)
746 for (int iy = 0; iy < nygrid; iy++)
747 for (int iz = 0; iz < nzgrid; iz++)
748 T_electron_first[ix][iy][iz] =
749 T_electron[ix][iy][iz];
750 do {
751 for (int ix = 0; ix < nxgrid; ix++)
752 for (int iy = 0; iy < nygrid; iy++)
753 for (int iz = 0; iz < nzgrid; iz++)
754 T_electron[ix][iy][iz] =
755 T_electron_first[ix][iy][iz];
756
757 stability_criterion = 1.0 -
758 2.0*inner_dt/el_specific_heat *
759 (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
760 if (stability_criterion < 0.0) {
761 inner_dt = 0.25*el_specific_heat /
762 (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
763 }
764 num_inner_timesteps = static_cast<unsigned int>(update->dt/inner_dt) + 1;
765 inner_dt = update->dt/double(num_inner_timesteps);
766 if (num_inner_timesteps > 1000000)
767 error->warning(FLERR,"Too many inner timesteps in fix ttm/mod");
768 for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
769 ith_inner_timestep++) {
770 for (int ix = 0; ix < nxgrid; ix++)
771 for (int iy = 0; iy < nygrid; iy++)
772 for (int iz = 0; iz < nzgrid; iz++)
773 T_electron_old[ix][iy][iz] =
774 T_electron[ix][iy][iz];
775 // compute new electron T profile
776 duration = duration + inner_dt;
777 for (int ix = 0; ix < nxgrid; ix++)
778 for (int iy = 0; iy < nygrid; iy++)
779 for (int iz = 0; iz < nzgrid; iz++) {
780 int right_x = ix + 1;
781 int right_y = iy + 1;
782 int right_z = iz + 1;
783 if (right_x == nxgrid) right_x = 0;
784 if (right_y == nygrid) right_y = 0;
785 if (right_z == nzgrid) right_z = 0;
786 int left_x = ix - 1;
787 int left_y = iy - 1;
788 int left_z = iz - 1;
789 if (left_x == -1) left_x = nxgrid - 1;
790 if (left_y == -1) left_y = nygrid - 1;
791 if (left_z == -1) left_z = nzgrid - 1;
792 double skin_layer_d = double(skin_layer);
793 double ix_d = double(ix);
794 double surface_d = double(t_surface_l);
795 mult_factor = 0.0;
796 if (duration < width) {
797 if (ix >= t_surface_l) mult_factor = (intensity/(dx*skin_layer_d))*exp((-1.0)*(ix_d - surface_d)/skin_layer_d);
798 }
799 if (ix < t_surface_l) net_energy_transfer_all[ix][iy][iz] = 0.0;
800 double cr_vac = 1;
801 if (T_electron_old[ix][iy][iz] == 0) cr_vac = 0;
802 double cr_v_l_x = 1;
803 if (T_electron_old[left_x][iy][iz] == 0) cr_v_l_x = 0;
804 double cr_v_r_x = 1;
805 if (T_electron_old[right_x][iy][iz] == 0) cr_v_r_x = 0;
806 double cr_v_l_y = 1;
807 if (T_electron_old[ix][left_y][iz] == 0) cr_v_l_y = 0;
808 double cr_v_r_y = 1;
809 if (T_electron_old[ix][right_y][iz] == 0) cr_v_r_y = 0;
810 double cr_v_l_z = 1;
811 if (T_electron_old[ix][iy][left_z] == 0) cr_v_l_z = 0;
812 double cr_v_r_z = 1;
813 if (T_electron_old[ix][iy][right_z] == 0) cr_v_r_z = 0;
814 if (cr_vac != 0) {
815 T_electron[ix][iy][iz] =
816 T_electron_old[ix][iy][iz] +
817 inner_dt/el_properties(T_electron_old[ix][iy][iz]).el_heat_capacity *
818 ((cr_v_r_x*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[right_x][iy][iz]/2.0).el_thermal_conductivity*
819 (T_electron_old[right_x][iy][iz]-T_electron_old[ix][iy][iz])/dx -
820 cr_v_l_x*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[left_x][iy][iz]/2.0).el_thermal_conductivity*
821 (T_electron_old[ix][iy][iz]-T_electron_old[left_x][iy][iz])/dx)/dx +
822 (cr_v_r_y*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][right_y][iz]/2.0).el_thermal_conductivity*
823 (T_electron_old[ix][right_y][iz]-T_electron_old[ix][iy][iz])/dy -
824 cr_v_l_y*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][left_y][iz]/2.0).el_thermal_conductivity*
825 (T_electron_old[ix][iy][iz]-T_electron_old[ix][left_y][iz])/dy)/dy +
826 (cr_v_r_z*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][iy][right_z]/2.0).el_thermal_conductivity*
827 (T_electron_old[ix][iy][right_z]-T_electron_old[ix][iy][iz])/dz -
828 cr_v_l_z*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][iy][left_z]/2.0).el_thermal_conductivity*
829 (T_electron_old[ix][iy][iz]-T_electron_old[ix][iy][left_z])/dz)/dz);
830 T_electron[ix][iy][iz]+=inner_dt/el_properties(T_electron[ix][iy][iz]).el_heat_capacity*
831 (mult_factor -
832 net_energy_transfer_all[ix][iy][iz]/del_vol);
833 }
834 else T_electron[ix][iy][iz] =
835 T_electron_old[ix][iy][iz];
836 if ((T_electron[ix][iy][iz] > 0.0) && (T_electron[ix][iy][iz] < electron_temperature_min))
837 T_electron[ix][iy][iz] = T_electron[ix][iy][iz] + 0.5*(electron_temperature_min - T_electron[ix][iy][iz]);
838
839 if (el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity > el_thermal_conductivity)
840 el_thermal_conductivity = el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity;
841 if ((T_electron[ix][iy][iz] > 0.0) && (el_properties(T_electron[ix][iy][iz]).el_heat_capacity < el_specific_heat))
842 el_specific_heat = el_properties(T_electron[ix][iy][iz]).el_heat_capacity;
843 }
844 }
845 stability_criterion = 1.0 -
846 2.0*inner_dt/el_specific_heat *
847 (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
848
849 } while (stability_criterion < 0.0);
850
851 // output of grid electron temperatures to file
852
853 if (outfile && (update->ntimestep % outevery == 0))
854 write_electron_temperatures(fmt::format("{}.{}", outfile, update->ntimestep));
855 }
856
857 /* ----------------------------------------------------------------------
858 memory usage of 3d grid
859 ------------------------------------------------------------------------- */
860
memory_usage()861 double FixTTMMod::memory_usage()
862 {
863 double bytes = 0.0;
864 bytes += (double)5*ngridtotal * sizeof(int);
865 bytes += (double)14*ngridtotal * sizeof(double);
866 return bytes;
867 }
868
869 /* ---------------------------------------------------------------------- */
870
grow_arrays(int ngrow)871 void FixTTMMod::grow_arrays(int ngrow)
872 {
873 memory->grow(flangevin,ngrow,3,"ttm/mod:flangevin");
874 }
875
876 /* ----------------------------------------------------------------------
877 return the energy of the electronic subsystem or the net_energy transfer
878 between the subsystems
879 ------------------------------------------------------------------------- */
880
compute_vector(int n)881 double FixTTMMod::compute_vector(int n)
882 {
883 double e_energy = 0.0;
884 double transfer_energy = 0.0;
885
886 double dx = domain->xprd/nxgrid;
887 double dy = domain->yprd/nygrid;
888 double dz = domain->zprd/nzgrid;
889 double del_vol = dx*dy*dz;
890
891 for (int ix = 0; ix < nxgrid; ix++)
892 for (int iy = 0; iy < nygrid; iy++)
893 for (int iz = 0; iz < nzgrid; iz++) {
894 e_energy += el_sp_heat_integral(T_electron[ix][iy][iz])*del_vol;
895 transfer_energy +=
896 net_energy_transfer_all[ix][iy][iz]*update->dt;
897 }
898
899 if (n == 0) return e_energy;
900 if (n == 1) return transfer_energy;
901 return 0.0;
902 }
903
904 /* ----------------------------------------------------------------------
905 pack entire state of Fix into one write
906 ------------------------------------------------------------------------- */
907
write_restart(FILE * fp)908 void FixTTMMod::write_restart(FILE *fp)
909 {
910 double *rlist;
911 memory->create(rlist,nxgrid*nygrid*nzgrid+4,"ttm/mod:rlist");
912
913 int n = 0;
914 rlist[n++] = nxgrid;
915 rlist[n++] = nygrid;
916 rlist[n++] = nzgrid;
917 rlist[n++] = seed;
918
919 for (int ix = 0; ix < nxgrid; ix++)
920 for (int iy = 0; iy < nygrid; iy++)
921 for (int iz = 0; iz < nzgrid; iz++)
922 rlist[n++] = T_electron[ix][iy][iz];
923
924 if (comm->me == 0) {
925 int size = n * sizeof(double);
926 fwrite(&size,sizeof(int),1,fp);
927 fwrite(rlist,sizeof(double),n,fp);
928 }
929
930 memory->destroy(rlist);
931 }
932
933 /* ----------------------------------------------------------------------
934 use state info from restart file to restart the Fix
935 ------------------------------------------------------------------------- */
936
restart(char * buf)937 void FixTTMMod::restart(char *buf)
938 {
939 int n = 0;
940 double *rlist = (double *) buf;
941
942 // check that restart grid size is same as current grid size
943
944 int nxgrid_old = static_cast<int> (rlist[n++]);
945 int nygrid_old = static_cast<int> (rlist[n++]);
946 int nzgrid_old = static_cast<int> (rlist[n++]);
947
948 if (nxgrid_old != nxgrid || nygrid_old != nygrid || nzgrid_old != nzgrid)
949 error->all(FLERR,"Must restart fix ttm with same grid size");
950
951 // change RN seed from initial seed, to avoid same Langevin factors
952 // just increment by 1, since for RanMars that is a new RN stream
953
954 seed = static_cast<int> (rlist[n++]) + 1;
955 delete random;
956 random = new RanMars(lmp,seed+comm->me);
957
958 // restore global frid values
959
960 for (int ix = 0; ix < nxgrid; ix++)
961 for (int iy = 0; iy < nygrid; iy++)
962 for (int iz = 0; iz < nzgrid; iz++)
963 T_electron[ix][iy][iz] = rlist[n++];
964 }
965
966 /* ----------------------------------------------------------------------
967 pack values in local atom-based arrays for restart file
968 ------------------------------------------------------------------------- */
969
pack_restart(int i,double * buf)970 int FixTTMMod::pack_restart(int i, double *buf)
971 {
972 // pack buf[0] this way because other fixes unpack it
973
974 buf[0] = 4;
975 buf[1] = flangevin[i][0];
976 buf[2] = flangevin[i][1];
977 buf[3] = flangevin[i][2];
978 return 4;
979 }
980
981 /* ----------------------------------------------------------------------
982 unpack values from atom->extra array to restart the fix
983 ------------------------------------------------------------------------- */
984
unpack_restart(int nlocal,int nth)985 void FixTTMMod::unpack_restart(int nlocal, int nth)
986 {
987 double **extra = atom->extra;
988
989 // skip to Nth set of extra values
990 // unpack the Nth first values this way because other fixes pack them
991
992 int m = 0;
993 for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
994 m++;
995
996 flangevin[nlocal][0] = extra[nlocal][m++];
997 flangevin[nlocal][1] = extra[nlocal][m++];
998 flangevin[nlocal][2] = extra[nlocal][m++];
999 }
1000
1001 /* ----------------------------------------------------------------------
1002 maxsize of any atom's restart data
1003 ------------------------------------------------------------------------- */
1004
maxsize_restart()1005 int FixTTMMod::maxsize_restart()
1006 {
1007 return 4;
1008 }
1009
1010 /* ----------------------------------------------------------------------
1011 size of atom nlocal's restart data
1012 ------------------------------------------------------------------------- */
1013
size_restart(int)1014 int FixTTMMod::size_restart(int /*nlocal*/)
1015 {
1016 return 4;
1017 }
1018