From d156263f540ebfcfebeddefe7eb0a6c9f78c1121 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julien=20Gu=C3=A9nol=C3=A9?= Date: Fri, 11 Aug 2017 13:39:13 +0200 Subject: [PATCH] Add relax_box function. Adapted documentation. --- doc/src/min_modify.txt | 14 ++++++- doc/src/min_style.txt | 2 + src/min.cpp | 18 +++++++++ src/min.h | 3 ++ src/min_adaptglok.cpp | 84 ++++++++++++++++++++++++++++++++++++++++-- src/min_adaptglok.h | 4 ++ 6 files changed, 119 insertions(+), 6 deletions(-) diff --git a/doc/src/min_modify.txt b/doc/src/min_modify.txt index b564d7696b..a97bf8eba5 100644 --- a/doc/src/min_modify.txt +++ b/doc/src/min_modify.txt @@ -28,7 +28,10 @@ keyword = {dmax} or {line} or {delaystep} or {dt_grow} or {dt_shrink} {tmin} value = {float} for adaptglok {integrator} value = {eulerimplicit} or {eulerexplicit} or {verlet} or {leapfrog} for adaptglok - {halfstepback} value = {yes} or {no} for adaptglok :pre + {halfstepback} value = {yes} or {no} for adaptglok :pre + {relaxbox} value = {no} or {iso} or {axial} for adaptglok + {relaxbox_modulus} value = {float} for adaptglok + {relaxbox_rate} value = {float} for adaptglok :ule [Examples:] @@ -84,6 +87,12 @@ to adapt them for particular cases. In particular for debugging purpose, it is possible to switch off the the backstep when downhill, by using and {halfstepback}. +The style {adaptglok} allows for a basic box relaxation to a stress-free +state by choosing {relaxbox} {iso} or {axial}. One should give an +approximated value for the bulk modulus {relaxbox_modulus} of the system. +The default value of the relaxation rate {relaxbox_rate} should be +optimum in most situations. + [Restrictions:] none [Related commands:] @@ -94,4 +103,5 @@ the backstep when downhill, by using and {halfstepback}. The option defaults are dmax = 0.1, line = quadratic, integrator = eulerimplicit, delaystep = 20, dt_grow = 1.1, dt_shrink = 0.5, alpha0 = 0.25, alpha_shrink = 0.99, -tmax = 2.0, tmin = 0.02, adaptstep = yes and halfstepback = yes. +tmax = 2.0, tmin = 0.02, adaptstep = yes, halfstepback = yes, +relaxbox = no, relaxbox_modulus = 1000000 and relaxbox_rate = 0.33. diff --git a/doc/src/min_style.txt b/doc/src/min_style.txt index e2a1c1fb73..b93e8f99e7 100644 --- a/doc/src/min_style.txt +++ b/doc/src/min_style.txt @@ -66,6 +66,8 @@ fit the published and unpublished optimizations of the method described in "(Bitzek)"_#Bitzek. It includes a temporization of the variable timestep, a backstep when downhill and different integration schemes. +Note that a basic scheme for having a stress-free simulation box is +implemented within that style. Either the {quickmin} and {fire} styles are useful in the context of nudged elastic band (NEB) calculations via the "neb"_neb.html command. diff --git a/src/min.cpp b/src/min.cpp index a650464eba..df92fbf153 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -63,7 +63,10 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp) tmax = 2.0; tmin = 0.02; integrator = 0; + relaxbox_modulus = 1000000; + relaxbox_rate = 0.33; halfstepback_flag = 1; + relaxbox_flag = 0; elist_global = elist_atom = NULL; vlist_global = vlist_atom = NULL; @@ -687,6 +690,21 @@ void Min::modify_params(int narg, char **arg) else if (strcmp(arg[iarg+1],"no") == 0) halfstepback_flag = 0; else error->all(FLERR,"Illegal min_modify command"); iarg += 2; + } else if (strcmp(arg[iarg],"relaxbox") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); + if (strcmp(arg[iarg+1],"no") == 0) relaxbox_flag = 0; + else if (strcmp(arg[iarg+1],"iso") == 0) relaxbox_flag = 1; + else if (strcmp(arg[iarg+1],"axial") == 0) relaxbox_flag = 2; + else error->all(FLERR,"Illegal min_modify command"); + iarg += 2; + } else if (strcmp(arg[iarg],"relaxbox_modulus") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); + relaxbox_modulus = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"relaxbox_rate") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); + relaxbox_rate = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; } else if (strcmp(arg[iarg],"integrator") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); if (strcmp(arg[iarg+1],"eulerimplicit") == 0) integrator = 0; diff --git a/src/min.h b/src/min.h index 53db7afd83..1a9e4be8b5 100644 --- a/src/min.h +++ b/src/min.h @@ -63,7 +63,10 @@ class Min : protected Pointers { double alpha0,alpha_shrink; // mixing velocities+forces coefficient (adaptglok) double tmax,tmin; // timestep max, min (adaptglok) int integrator; // choose the style of time integrator (adaptglok) + double relaxbox_modulus; // Bulk modulus used for box relax (for adaptglok) + double relaxbox_rate; // for box relaxation to 0 pressure (for adaptglok) int halfstepback_flag; // 1: half step backward when v.f <= 0.0 (adaptglok) + int relaxbox_flag; // 1: box relaxation iso; 2: axial int nelist_global,nelist_atom; // # of PE,virial computes to check int nvlist_global,nvlist_atom; diff --git a/src/min_adaptglok.cpp b/src/min_adaptglok.cpp index 199dee743c..08ea84e510 100644 --- a/src/min_adaptglok.cpp +++ b/src/min_adaptglok.cpp @@ -20,6 +20,10 @@ #include "output.h" #include "timer.h" #include "error.h" +#include "variable.h" +#include "modify.h" +#include "compute.h" +#include "domain.h" using namespace LAMMPS_NS; @@ -42,6 +46,13 @@ void MinAdaptGlok::init() dtmin = tmin * dt; alpha = alpha0; last_negative = ntimestep_fire = update->ntimestep; + + if (relaxbox_flag){ + int icompute = modify->find_compute("thermo_temp"); + temperature = modify->compute[icompute]; + icompute = modify->find_compute("thermo_press"); + pressure = modify->compute[icompute]; + } } /* ---------------------------------------------------------------------- */ @@ -69,6 +80,68 @@ void MinAdaptGlok::reset_vectors() if (nvec) fvec = atom->f[0]; } +/* ---------------------------------------------------------------------- + save current box state for converting atoms to lamda coords +------------------------------------------------------------------------- */ + +void MinAdaptGlok::save_box_state() +{ + boxlo[0] = domain->boxlo[0]; + boxlo[1] = domain->boxlo[1]; + boxlo[2] = domain->boxlo[2]; + + for (int i = 0; i < 6; i++) + h_inv[i] = domain->h_inv[i]; +} + +/* ---------------------------------------------------------------------- + deform the simulation box and remap the particles +------------------------------------------------------------------------- */ + +void MinAdaptGlok::relax_box() +{ + int i,n; + + // rescale simulation box from linesearch starting point + // scale atom coords for all atoms or only for fix group atoms + + double **x = atom->x; + double epsilon; + double *current_pressure_v; + int *mask = atom->mask; + n = atom->nlocal + atom->nghost; + save_box_state(); + + // convert pertinent atoms and rigid bodies to lamda coords + + domain->x2lamda(n); + + // ensure the virial is tallied + + update->vflag_global = update->ntimestep; + + // compute pressure and change simulation box + + pressure->compute_scalar(); + pressure->compute_vector(); + epsilon = pressure->scalar / relaxbox_modulus; + current_pressure_v = pressure->vector; + for (int i = 0; i < 3; i++) { + if (relaxbox_flag == 2) epsilon = current_pressure_v[i] / relaxbox_modulus; + domain->boxhi[i] += domain->boxhi[i] * epsilon * relaxbox_rate;; + } + + // reset global and local box to new size/shape + + domain->set_global_box(); + domain->set_local_box(); + + // convert atoms and back to box coords + + domain->lamda2x(n); + save_box_state(); +} + /* ---------------------------------------------------------------------- */ int MinAdaptGlok::iterate(int maxiter) @@ -123,7 +196,11 @@ int MinAdaptGlok::iterate(int maxiter) ntimestep = ++update->ntimestep; niter++; - // vdotfall = v dot f + // Relax the simulation box + + if (relaxbox_flag) relax_box(); + + // vdotfall = v dot f // Euler || Leap Frog integration @@ -132,6 +209,7 @@ int MinAdaptGlok::iterate(int maxiter) double **f = atom->f; } int nlocal = atom->nlocal; + double **x = atom->x; vdotf = 0.0; for (int i = 0; i < nlocal; i++) @@ -243,9 +321,7 @@ int MinAdaptGlok::iterate(int maxiter) MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld); } - double **x = atom->x; - - // Semi-implicit Euler integration + // Semi-implicit Euler integration if (integrator == 0) { diff --git a/src/min_adaptglok.h b/src/min_adaptglok.h index 8e1bf33062..8ae8f0a947 100644 --- a/src/min_adaptglok.h +++ b/src/min_adaptglok.h @@ -31,12 +31,16 @@ class MinAdaptGlok : public Min { void init(); void setup_style(); void reset_vectors(); + void save_box_state(); + void relax_box(); int iterate(int); private: double dt,dtmax,dtmin,dtinit; double alpha; bigint last_negative,ntimestep_fire; + class Compute *temperature,*pressure; + double boxlo[3],h_inv[6]; }; }