ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Change fire to fire/old and fire2 to fire. Implement normstyle in fire. Update author affiliation.

Browse Source
This commit is contained in:
Julien Guénolé
2020-01-28 14:51:23 +01:00
parent ea24ec8d6a
commit 197ba62cd9
5 changed files with 566 additions and 564 deletions
Download Patch File Download Diff File Expand all files Collapse all files

339
src/min_fire.cpp
View File

@ -11,16 +11,26 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "min_fire.h"
#include <mpi.h>
/* ----------------------------------------------------------------------
Contributing authors: Julien Guénolé, CNRS and
Erik Bitzek, FAU Erlangen-Nuernberg
------------------------------------------------------------------------- */
#include <cmath>
#include "min_fire.h"
#include "universe.h"
#include "atom.h"
#include "error.h"
#include "force.h"
#include "update.h"
#include "output.h"
#include "timer.h"
#include "error.h"
#include "variable.h"
#include "modify.h"
#include "compute.h"
#include "domain.h"
#include "neighbor.h"
#include "comm.h"
using namespace LAMMPS_NS;
@ -28,13 +38,6 @@ using namespace LAMMPS_NS;
#define EPS_ENERGY 1.0e-8
#define DELAYSTEP 5
#define DT_GROW 1.1
#define DT_SHRINK 0.5
#define ALPHA0 0.1
#define ALPHA_SHRINK 0.99
#define TMAX 10.0
/* ---------------------------------------------------------------------- */
MinFire::MinFire(LAMMPS *lmp) : Min(lmp) {}
@ -45,10 +48,18 @@ void MinFire::init()
{
Min::init();
// simple parameters validation
if (tmax < tmin) error->all(FLERR,"tmax has to be larger than tmin");
if (dtgrow < 1.0) error->all(FLERR,"dtgrow has to be larger than 1.0");
if (dtshrink > 1.0) error->all(FLERR,"dtshrink has to be smaller than 1.0");
dt = update->dt;
dtmax = TMAX * dt;
alpha = ALPHA0;
last_negative = update->ntimestep;
dtmax = tmax * dt;
dtmin = tmin * dt;
alpha = alpha0;
last_negative = ntimestep_start = update->ntimestep;
vdotf_negatif = 0;
}
/* ---------------------------------------------------------------------- */
@ -58,6 +69,22 @@ void MinFire::setup_style()
double **v = atom->v;
int nlocal = atom->nlocal;
// print the parameters used within fire into the log
const char *s1[] = {"eulerimplicit","verlet","leapfrog","eulerexplicit"};
const char *s2[] = {"no","yes"};
if (comm->me == 0 && logfile) {
fprintf(logfile," Parameters for fire: \n"
" dmax delaystep dtgrow dtshrink alpha0 alphashrink tmax tmin "
" integrator halfstepback relaxbox relaxbox_mod relaxbox_rate ptol \n"
" %4g %9i %6g %8g %6g %11g %4g %4g %13s %12s \n",
dmax, delaystep, dtgrow, dtshrink, alpha0, alphashrink, tmax, tmin,
s1[integrator], s2[halfstepback_flag]);
}
// initialize the velocities
for (int i = 0; i < nlocal; i++)
v[i][0] = v[i][1] = v[i][2] = 0.0;
}
@ -88,6 +115,39 @@ int MinFire::iterate(int maxiter)
alpha_final = 0.0;
// Leap Frog integration initialization
if (integrator == 2) {
double **f = atom->f;
double **v = atom->v;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
energy_force(0);
neval++;
dtf = -0.5 * dt * force->ftm2v;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / rmass[i];
v[i][0] = dtfm * f[i][0];
v[i][1] = dtfm * f[i][1];
v[i][2] = dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / mass[type[i]];
v[i][0] = dtfm * f[i][0];
v[i][1] = dtfm * f[i][1];
v[i][2] = dtfm * f[i][2];
}
}
}
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
@ -96,11 +156,17 @@ int MinFire::iterate(int maxiter)
ntimestep = ++update->ntimestep;
niter++;
// vdotfall = v dot f
// pointers
int nlocal = atom->nlocal;
double **v = atom->v;
double **f = atom->f;
int nlocal = atom->nlocal;
double **x = atom->x;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
// vdotfall = v dot f
vdotf = 0.0;
for (int i = 0; i < nlocal; i++)
@ -118,11 +184,14 @@ int MinFire::iterate(int maxiter)
// if (v dot f) > 0:
// v = (1-alpha) v + alpha |v| Fhat
// |v| = length of v, Fhat = unit f
// if more than DELAYSTEP since v dot f was negative:
// increase timestep and decrease alpha
// Only: (1-alpha) and alpha |v| Fhat is calculated here
// the modificatin of v is made wihtin the integration, after v update
// if more than delaystep since v dot f was negative:
// increase timestep, update global timestep and decrease alpha
if (vdotfall > 0.0) {
vdotv = 0.0;
vdotf_negatif = 0;
for (int i = 0; i < nlocal; i++)
vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world);
@ -149,36 +218,59 @@ int MinFire::iterate(int maxiter)
}
scale1 = 1.0 - alpha;
if (fdotfall == 0.0) scale2 = 0.0;
if (fdotfall <= 1e-20) scale2 = 0.0;
else scale2 = alpha * sqrt(vdotvall/fdotfall);
for (int i = 0; i < nlocal; i++) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
if (ntimestep - last_negative > delaystep) {
dt = MIN(dt*dtgrow,dtmax);
update->dt = dt;
alpha *= alphashrink;
}
if (ntimestep - last_negative > DELAYSTEP) {
dt = MIN(dt*DT_GROW,dtmax);
alpha *= ALPHA_SHRINK;
}
// else (v dot f) <= 0:
// decrease timestep, reset alpha, set v = 0
// else (v dot f) <= 0
// if more than delaystep since starting the relaxation:
// reset alpha
// if dt*dtshrink > dtmin:
// decrease timestep
// update global timestep (for thermo output)
// half step back within the dynamics: x(t) = x(t-0.5*dt)
// reset velocities: v = 0
} else {
last_negative = ntimestep;
dt *= DT_SHRINK;
alpha = ALPHA0;
int delayflag = 1;
if (ntimestep - ntimestep_start < delaystep && delaystep_start_flag)
delayflag = 0;
if (delayflag) {
alpha = alpha0;
if (dt*dtshrink >= dtmin) {
dt *= dtshrink;
update->dt = dt;
}
}
// stopping criterion while stuck in a local bassin of the PES
vdotf_negatif++;
if (max_vdotf_negatif > 0 && vdotf_negatif > max_vdotf_negatif)
return MAXVDOTF;
// inertia correction
if (halfstepback_flag) {
for (int i = 0; i < nlocal; i++) {
x[i][0] -= 0.5 * dtv * v[i][0];
x[i][1] -= 0.5 * dtv * v[i][1];
x[i][2] -= 0.5 * dtv * v[i][2];
}
}
for (int i = 0; i < nlocal; i++)
v[i][0] = v[i][1] = v[i][2] = 0.0;
}
// limit timestep so no particle moves further than dmax
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
dtvone = dt;
for (int i = 0; i < nlocal; i++) {
@ -186,6 +278,7 @@ int MinFire::iterate(int maxiter)
vmax = MAX(vmax,fabs(v[i][2]));
if (dtvone*vmax > dmax) dtvone = dmax/vmax;
}
MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world);
// min dtv over replicas, if necessary
@ -196,43 +289,161 @@ int MinFire::iterate(int maxiter)
MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld);
}
dtf = dtv * force->ftm2v;
// Dynamic integration scheme:
// 0: semi-implicit Euler
// 1: velocity Verlet
// 2: leapfrog (initial half step before the iteration loop)
// 3: explicit Euler
// Euler integration step
// Semi-implicit Euler OR Leap Frog integration
double **x = atom->x;
if (integrator == 0 || integrator == 2) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / rmass[i];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
dtf = dtv * force->ftm2v;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vdotfall > 0.0) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
}
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vdotfall > 0.0) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
}
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / mass[type[i]];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
eprevious = ecurrent;
ecurrent = energy_force(0);
neval++;
// Velocity Verlet integration
} else if (integrator == 1) {
dtf = 0.5 * dtv * force->ftm2v;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vdotfall > 0.0) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
}
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
if (vdotfall > 0.0) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
}
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
eprevious = ecurrent;
ecurrent = energy_force(0);
neval++;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
// Standard Euler integration
} else if (integrator == 3) {
dtf = dtv * force->ftm2v;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / rmass[i];
if (vdotfall > 0.0) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
}
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++) {
dtfm = dtf / mass[type[i]];
if (vdotfall > 0.0) {
v[i][0] = scale1*v[i][0] + scale2*f[i][0];
v[i][1] = scale1*v[i][1] + scale2*f[i][1];
v[i][2] = scale1*v[i][2] + scale2*f[i][2];
}
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
eprevious = ecurrent;
ecurrent = energy_force(0);
neval++;
}
eprevious = ecurrent;
ecurrent = energy_force(0);
neval++;
// energy tolerance criterion
// only check after DELAYSTEP elapsed since velocties reset to 0
// only check after delaystep elapsed since velocties reset to 0
// sync across replicas if running multi-replica minimization
// reset the timestep to the initial value
if (update->etol > 0.0 && ntimestep-last_negative > DELAYSTEP) {
if (update->etol > 0.0 && ntimestep-last_negative > delaystep) {
if (update->multireplica == 0) {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
@ -243,12 +454,14 @@ int MinFire::iterate(int maxiter)
flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return ETOL;
if (flagall == 0)
return ETOL;
}
}
// force tolerance criterion
// sync across replicas if running multi-replica minimization
// reset the timestep to the initial value
fdotf = 0.0;
if (update->ftol > 0.0) {