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: Julien Guénolé, CNRS and
17 Erik Bitzek, FAU Erlangen-Nuernberg
18 ------------------------------------------------------------------------- */
19
20 #include "min_fire.h"
21
22 #include "atom.h"
23 #include "comm.h"
24 #include "error.h"
25 #include "force.h"
26 #include "output.h"
27 #include "timer.h"
28 #include "universe.h"
29 #include "update.h"
30
31 #include <cmath>
32
33 using namespace LAMMPS_NS;
34
35 // EPS_ENERGY = minimum normalization for energy tolerance
36
37 #define EPS_ENERGY 1.0e-8
38
39 /* ---------------------------------------------------------------------- */
40
MinFire(LAMMPS * lmp)41 MinFire::MinFire(LAMMPS *lmp) : Min(lmp) {}
42
43 /* ---------------------------------------------------------------------- */
44
init()45 void MinFire::init()
46 {
47 Min::init();
48
49 // simple parameters validation
50
51 if (tmax < tmin) error->all(FLERR,"tmax has to be larger than tmin");
52 if (dtgrow < 1.0) error->all(FLERR,"dtgrow has to be larger than 1.0");
53 if (dtshrink > 1.0) error->all(FLERR,"dtshrink has to be smaller than 1.0");
54
55 dt = update->dt;
56 dtmax = tmax * dt;
57 dtmin = tmin * dt;
58 alpha = alpha0;
59 last_negative = ntimestep_start = update->ntimestep;
60 vdotf_negatif = 0;
61 }
62
63 /* ---------------------------------------------------------------------- */
64
setup_style()65 void MinFire::setup_style()
66 {
67 double **v = atom->v;
68 int nlocal = atom->nlocal;
69
70 // print the parameters used within fire into the log
71
72 const char *s1[] = {"eulerimplicit","verlet","leapfrog","eulerexplicit"};
73 const char *s2[] = {"no","yes"};
74
75 if (comm->me == 0 && logfile) {
76 fprintf(logfile," Parameters for fire: \n"
77 " dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin "
78 " integrator halfstepback \n"
79 " %4g %9i %6g %8g %6g %11g %4g %4g %13s %12s \n",
80 dmax, delaystep, dtgrow, dtshrink, alpha0, alphashrink, tmax, tmin,
81 s1[integrator], s2[halfstepback_flag]);
82 }
83
84 // initialize the velocities
85
86 for (int i = 0; i < nlocal; i++)
87 v[i][0] = v[i][1] = v[i][2] = 0.0;
88 flagv0 = 1;
89 }
90
91 /* ----------------------------------------------------------------------
92 set current vector lengths and pointers
93 called after atoms have migrated
94 ------------------------------------------------------------------------- */
95
reset_vectors()96 void MinFire::reset_vectors()
97 {
98 // atomic dof
99
100 nvec = 3 * atom->nlocal;
101 if (nvec) xvec = atom->x[0];
102 if (nvec) fvec = atom->f[0];
103 }
104
105 /* ---------------------------------------------------------------------- */
106
iterate(int maxiter)107 int MinFire::iterate(int maxiter)
108 {
109 bigint ntimestep;
110 double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall;
111 double scale1,scale2;
112 double dtvone,dtv,dtf,dtfm;
113 int flag,flagall;
114
115 alpha_final = 0.0;
116
117 // Leap Frog integration initialization
118
119 if (integrator == 2) {
120
121 double **f = atom->f;
122 double **v = atom->v;
123 double *rmass = atom->rmass;
124 double *mass = atom->mass;
125 int *type = atom->type;
126 int nlocal = atom->nlocal;
127
128 energy_force(0);
129 neval++;
130
131 dtf = -0.5 * dt * force->ftm2v;
132
133 if (rmass) {
134 for (int i = 0; i < nlocal; i++) {
135 dtfm = dtf / rmass[i];
136 v[i][0] = dtfm * f[i][0];
137 v[i][1] = dtfm * f[i][1];
138 v[i][2] = dtfm * f[i][2];
139 }
140 } else {
141 for (int i = 0; i < nlocal; i++) {
142 dtfm = dtf / mass[type[i]];
143 v[i][0] = dtfm * f[i][0];
144 v[i][1] = dtfm * f[i][1];
145 v[i][2] = dtfm * f[i][2];
146 }
147 }
148 }
149
150 for (int iter = 0; iter < maxiter; iter++) {
151
152 if (timer->check_timeout(niter))
153 return TIMEOUT;
154
155 ntimestep = ++update->ntimestep;
156 niter++;
157
158 // pointers
159
160 int nlocal = atom->nlocal;
161 double **v = atom->v;
162 double **f = atom->f;
163 double **x = atom->x;
164 double *rmass = atom->rmass;
165 double *mass = atom->mass;
166 int *type = atom->type;
167
168 // vdotfall = v dot f
169
170 vdotf = 0.0;
171 for (int i = 0; i < nlocal; i++)
172 vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2];
173 MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world);
174
175 // sum vdotf over replicas, if necessary
176 // this communicator would be invalid for multiprocess replicas
177
178 if (update->multireplica == 1) {
179 vdotf = vdotfall;
180 MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
181 }
182
183 // if (v dot f) > 0:
184 // v = (1-alpha) v + alpha |v| Fhat
185 // |v| = length of v, Fhat = unit f
186 // Only: (1-alpha) and alpha |v| Fhat is calculated here
187 // the modificatin of v is made wihtin the integration, after v update
188 // if more than delaystep since v dot f was negative:
189 // increase timestep, update global timestep and decrease alpha
190
191 if (vdotfall > 0.0) {
192 vdotv = 0.0;
193 vdotf_negatif = 0;
194 for (int i = 0; i < nlocal; i++)
195 vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
196 MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
197
198 // sum vdotv over replicas, if necessary
199 // this communicator would be invalid for multiprocess replicas
200
201 if (update->multireplica == 1) {
202 vdotv = vdotvall;
203 MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
204 }
205
206 fdotf = 0.0;
207 for (int i = 0; i < nlocal; i++)
208 fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
209 MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world);
210
211 // sum fdotf over replicas, if necessary
212 // this communicator would be invalid for multiprocess replicas
213
214 if (update->multireplica == 1) {
215 fdotf = fdotfall;
216 MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
217 }
218
219 scale1 = 1.0 - alpha;
220 if (fdotfall <= 1e-20) scale2 = 0.0;
221 else scale2 = alpha * sqrt(vdotvall/fdotfall);
222
223 if (ntimestep - last_negative > delaystep) {
224 dt = MIN(dt*dtgrow,dtmax);
225 update->dt = dt;
226 alpha *= alphashrink;
227 }
228
229 // else (v dot f) <= 0
230 // if more than delaystep since starting the relaxation:
231 // reset alpha
232 // if dt*dtshrink > dtmin:
233 // decrease timestep
234 // update global timestep (for thermo output)
235 // half step back within the dynamics: x(t) = x(t-0.5*dt)
236 // reset velocities: v = 0
237
238 } else {
239 last_negative = ntimestep;
240 int delayflag = 1;
241 if (ntimestep - ntimestep_start < delaystep && delaystep_start_flag)
242 delayflag = 0;
243 if (delayflag) {
244 alpha = alpha0;
245 if (dt*dtshrink >= dtmin) {
246 dt *= dtshrink;
247 update->dt = dt;
248 }
249 }
250
251 // stopping criterion while stuck in a local bassin of the PES
252
253 vdotf_negatif++;
254 if (max_vdotf_negatif > 0 && vdotf_negatif > max_vdotf_negatif)
255 return MAXVDOTF;
256
257 // inertia correction
258
259 if (halfstepback_flag) {
260 for (int i = 0; i < nlocal; i++) {
261 x[i][0] -= 0.5 * dt * v[i][0];
262 x[i][1] -= 0.5 * dt * v[i][1];
263 x[i][2] -= 0.5 * dt * v[i][2];
264 }
265 }
266
267 for (int i = 0; i < nlocal; i++)
268 v[i][0] = v[i][1] = v[i][2] = 0.0;
269 flagv0 = 1;
270 }
271
272 // evaluates velocties to estimate wether dtv has to be limited
273 // required when v have been reset
274
275 if (flagv0) {
276 dtf = dt * force->ftm2v;
277 energy_force(0);
278 neval++;
279
280 if (rmass) {
281 for (int i = 0; i < nlocal; i++) {
282 dtfm = dtf / rmass[i];
283 v[i][0] = dtfm * f[i][0];
284 v[i][1] = dtfm * f[i][1];
285 v[i][2] = dtfm * f[i][2];
286 }
287 } else {
288 for (int i = 0; i < nlocal; i++) {
289 dtfm = dtf / mass[type[i]];
290 v[i][0] = dtfm * f[i][0];
291 v[i][1] = dtfm * f[i][1];
292 v[i][2] = dtfm * f[i][2];
293 }
294 }
295 }
296
297 // limit timestep so no particle moves further than dmax
298
299 dtvone = dt;
300
301 for (int i = 0; i < nlocal; i++) {
302 vmax = MAX(fabs(v[i][0]),fabs(v[i][1]));
303 vmax = MAX(vmax,fabs(v[i][2]));
304 if (dtvone*vmax > dmax) dtvone = dmax/vmax;
305 }
306
307 MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
308
309 // reset velocities when necessary
310
311 if (flagv0) {
312 for (int i = 0; i < nlocal; i++)
313 v[i][0] = v[i][1] = v[i][2] = 0.0;
314 }
315
316 // min dtv over replicas, if necessary
317 // this communicator would be invalid for multiprocess replicas
318
319 if (update->multireplica == 1) {
320 dtvone = dtv;
321 MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
322 }
323
324 // Dynamic integration scheme:
325 // 0: semi-implicit Euler
326 // 1: velocity Verlet
327 // 2: leapfrog (initial half step before the iteration loop)
328 // 3: explicit Euler
329
330 // Semi-implicit Euler OR Leap Frog integration
331
332 if (integrator == 0 || integrator == 2) {
333
334 dtf = dtv * force->ftm2v;
335
336 if (rmass) {
337 for (int i = 0; i < nlocal; i++) {
338 dtfm = dtf / rmass[i];
339 v[i][0] += dtfm * f[i][0];
340 v[i][1] += dtfm * f[i][1];
341 v[i][2] += dtfm * f[i][2];
342 if (vdotfall > 0.0) {
343 v[i][0] = scale1*v[i][0] + scale2*f[i][0];
344 v[i][1] = scale1*v[i][1] + scale2*f[i][1];
345 v[i][2] = scale1*v[i][2] + scale2*f[i][2];
346 }
347 x[i][0] += dtv * v[i][0];
348 x[i][1] += dtv * v[i][1];
349 x[i][2] += dtv * v[i][2];
350 }
351 } else {
352 for (int i = 0; i < nlocal; i++) {
353 dtfm = dtf / mass[type[i]];
354 v[i][0] += dtfm * f[i][0];
355 v[i][1] += dtfm * f[i][1];
356 v[i][2] += dtfm * f[i][2];
357 if (vdotfall > 0.0) {
358 v[i][0] = scale1*v[i][0] + scale2*f[i][0];
359 v[i][1] = scale1*v[i][1] + scale2*f[i][1];
360 v[i][2] = scale1*v[i][2] + scale2*f[i][2];
361 }
362 x[i][0] += dtv * v[i][0];
363 x[i][1] += dtv * v[i][1];
364 x[i][2] += dtv * v[i][2];
365 }
366 }
367
368 eprevious = ecurrent;
369 ecurrent = energy_force(0);
370 neval++;
371
372 // Velocity Verlet integration
373
374 } else if (integrator == 1) {
375
376 dtf = 0.5 * dtv * force->ftm2v;
377
378 if (rmass) {
379 for (int i = 0; i < nlocal; i++) {
380 dtfm = dtf / rmass[i];
381 v[i][0] += dtfm * f[i][0];
382 v[i][1] += dtfm * f[i][1];
383 v[i][2] += dtfm * f[i][2];
384 if (vdotfall > 0.0) {
385 v[i][0] = scale1*v[i][0] + scale2*f[i][0];
386 v[i][1] = scale1*v[i][1] + scale2*f[i][1];
387 v[i][2] = scale1*v[i][2] + scale2*f[i][2];
388 }
389 x[i][0] += dtv * v[i][0];
390 x[i][1] += dtv * v[i][1];
391 x[i][2] += dtv * v[i][2];
392 }
393 } else {
394 for (int i = 0; i < nlocal; i++) {
395 dtfm = dtf / mass[type[i]];
396 v[i][0] += dtfm * f[i][0];
397 v[i][1] += dtfm * f[i][1];
398 v[i][2] += dtfm * f[i][2];
399 if (vdotfall > 0.0) {
400 v[i][0] = scale1*v[i][0] + scale2*f[i][0];
401 v[i][1] = scale1*v[i][1] + scale2*f[i][1];
402 v[i][2] = scale1*v[i][2] + scale2*f[i][2];
403 }
404 x[i][0] += dtv * v[i][0];
405 x[i][1] += dtv * v[i][1];
406 x[i][2] += dtv * v[i][2];
407 }
408 }
409
410 eprevious = ecurrent;
411 ecurrent = energy_force(0);
412 neval++;
413
414 if (rmass) {
415 for (int i = 0; i < nlocal; i++) {
416 dtfm = dtf / rmass[i];
417 v[i][0] += dtfm * f[i][0];
418 v[i][1] += dtfm * f[i][1];
419 v[i][2] += dtfm * f[i][2];
420 }
421 } else {
422 for (int i = 0; i < nlocal; i++) {
423 dtfm = dtf / mass[type[i]];
424 v[i][0] += dtfm * f[i][0];
425 v[i][1] += dtfm * f[i][1];
426 v[i][2] += dtfm * f[i][2];
427 }
428 }
429
430 // Standard Euler integration
431
432 } else if (integrator == 3) {
433
434 dtf = dtv * force->ftm2v;
435
436 if (rmass) {
437 for (int i = 0; i < nlocal; i++) {
438 dtfm = dtf / rmass[i];
439 if (vdotfall > 0.0) {
440 v[i][0] = scale1*v[i][0] + scale2*f[i][0];
441 v[i][1] = scale1*v[i][1] + scale2*f[i][1];
442 v[i][2] = scale1*v[i][2] + scale2*f[i][2];
443 }
444 x[i][0] += dtv * v[i][0];
445 x[i][1] += dtv * v[i][1];
446 x[i][2] += dtv * v[i][2];
447 v[i][0] += dtfm * f[i][0];
448 v[i][1] += dtfm * f[i][1];
449 v[i][2] += dtfm * f[i][2];
450 }
451 } else {
452 for (int i = 0; i < nlocal; i++) {
453 dtfm = dtf / mass[type[i]];
454 if (vdotfall > 0.0) {
455 v[i][0] = scale1*v[i][0] + scale2*f[i][0];
456 v[i][1] = scale1*v[i][1] + scale2*f[i][1];
457 v[i][2] = scale1*v[i][2] + scale2*f[i][2];
458 }
459 x[i][0] += dtv * v[i][0];
460 x[i][1] += dtv * v[i][1];
461 x[i][2] += dtv * v[i][2];
462 v[i][0] += dtfm * f[i][0];
463 v[i][1] += dtfm * f[i][1];
464 v[i][2] += dtfm * f[i][2];
465 }
466 }
467
468 eprevious = ecurrent;
469 ecurrent = energy_force(0);
470 neval++;
471 }
472
473 // velocities have been evaluated
474
475 flagv0 = 0;
476
477 // energy tolerance criterion
478 // only check after delaystep elapsed since velocties reset to 0
479 // sync across replicas if running multi-replica minimization
480 // reset the timestep to the initial value
481
482 if (update->etol > 0.0 && ntimestep-last_negative > delaystep) {
483 if (update->multireplica == 0) {
484 if (fabs(ecurrent-eprevious) <
485 update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
486 return ETOL;
487 } else {
488 if (fabs(ecurrent-eprevious) <
489 update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
490 flag = 0;
491 else flag = 1;
492 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
493 if (flagall == 0)
494 return ETOL;
495 }
496 }
497
498 // force tolerance criterion
499 // sync across replicas if running multi-replica minimization
500 // reset the timestep to the initial value
501
502 fdotf = 0.0;
503 if (update->ftol > 0.0) {
504 if (normstyle == MAX) fdotf = fnorm_max(); // max force norm
505 else if (normstyle == INF) fdotf = fnorm_inf(); // inf force norm
506 else if (normstyle == TWO) fdotf = fnorm_sqr(); // Euclidean force 2-norm
507 else error->all(FLERR,"Illegal min_modify command");
508 if (update->multireplica == 0) {
509 if (fdotf < update->ftol*update->ftol) return FTOL;
510 } else {
511 if (fdotf < update->ftol*update->ftol) flag = 0;
512 else flag = 1;
513 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
514 if (flagall == 0) return FTOL;
515 }
516 }
517
518 // output for thermo, dump, restart files
519
520 if (output->next == ntimestep) {
521 timer->stamp();
522 output->write(ntimestep);
523 timer->stamp(Timer::OUTPUT);
524 }
525 }
526
527 return MAXITER;
528 }
529