Page MenuHomec4science
Files F107087313

domain_lammps_minimize.cc-old
No OneTemporary
Actions

File Metadata

Created
Fri, Apr 4, 10:38

domain_lammps_minimize.cc-old
View Options

/**
* @file domain_lammps_minimize.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Thu Jul 31 22:41:23 2014
*
* @brief This is the generic LAMMPS model capable of static solve
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but
* WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#define TIMER
#define NEED_LAMMPS
#define LOCAL_MODULE MOD_MD
//#define CHECK_STABILITY
/* -------------------------------------------------------------------------- */
#include "domain_lammps_minimize.hh"
#include "communicator.hh"
#include "import_lammps.hh"
#include "lammps_common.hh"
#include "lm_common.hh"
#include "min_cg.h"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
// LAMMPS include files
#include "angle.h"
#include "bond.h"
#include "dihedral.h"
#include "improper.h"
#include "integrate.h"
#include "kspace.h"
#include "pair.h"
#include "thermo.h"
#include "verlet.h"
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
/* patch for till min cg */
#define EPS_ENERGY 1.0e-8
#define BACKTRACK_SLOPE 0.5
#define ALPHA_REDUCE 0.5
enum { FAIL, MAXITER, MAXEVAL, ETOL, FTOL }; // same as in other min classes
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammpsMinimize<Dim>::init() {
if (this->getGeomConstrained().getDim() == 1)
LM_FATAL("Lammps cannot work in 1D");
/* --------------------------------------------------------------------------
*/
// extracted from minimize.cpp in LAMMPS project
/* --------------------------------------------------------------------------
*/
if (this->domain->box_exist == 0)
this->error->all("Minimize command before simulation box is defined");
this->update->etol = this->mincg_etol;
this->update->ftol = this->mincg_ftol;
this->update->nsteps = nb_step;
// this->update->nsteps = nb_step_next_event;
this->update->max_eval = this->mincg_maxeval * nb_step;
if (this->update->etol < 0.0 || this->update->ftol < 0.0)
this->error->all("Illegal minimize command");
this->update->whichflag = 2;
this->update->beginstep = this->update->firststep = this->update->ntimestep;
// this->update->endstep = this->update->laststep = this->update->firststep +
// this->update->nsteps;
this->update->endstep = this->update->laststep =
this->update->firststep + this->update->nsteps;
if (this->update->laststep < 0 || this->update->laststep > UINT_MAX)
this->error->all("Too many iterations");
LAMMPS_NS::LAMMPS::init();
this->setCurrentStep(0);
this->update->minimize->setup();
this->initDOFs();
LAMMPS_NS::MinCG *min_ptr =
static_cast<LAMMPS_NS::MinCG *>(this->update->minimize);
min_ptr->reset_vectors();
/* --------------------------------------------------------------------------
*/
// init of min_cg : extracted from MinCG::iterate at min_cg.cpp file
/* --------------------------------------------------------------------------
*/
int i, m, n; //,fail,ntimestep;
double *fatom, *gatom, *hatom;
// initialize working vectors
for (i = 0; i < this->update->minimize->nvec; i++) {
min_ptr->h[i] = min_ptr->fvec[i];
min_ptr->g[i] = min_ptr->fvec[i];
}
Geometry &geom_constrained = this->getGeomConstrained();
for (int j = 0; j < min_ptr->nvec / 3; ++j) {
Real x = lammps_main_object->atom->x0[j][0];
Real y = lammps_main_object->atom->x0[j][1];
Real z = lammps_main_object->atom->x0[j][2];
if (geom_constrained.contains(x, y, z)) {
min_ptr->g[3 * j + 0] = 0;
min_ptr->g[3 * j + 1] = 0;
min_ptr->g[3 * j + 2] = 0;
min_ptr->fvec[3 * j + 0] = 0;
min_ptr->fvec[3 * j + 1] = 0;
min_ptr->fvec[3 * j + 2] = 0;
min_ptr->h[3 * j] = 0;
min_ptr->h[3 * j + 1] = 0;
min_ptr->h[3 * j + 2] = 0;
}
}
if (min_ptr->nextra_atom)
for (m = 0; m < min_ptr->nextra_atom; m++) {
fatom = min_ptr->fextra_atom[m];
gatom = min_ptr->gextra_atom[m];
hatom = min_ptr->hextra_atom[m];
n = min_ptr->extra_nlen[m];
for (i = 0; i < n; i++)
hatom[i] = gatom[i] = fatom[i];
}
if (min_ptr->nextra_global)
for (i = 0; i < min_ptr->nextra_global; i++)
min_ptr->hextra[i] = min_ptr->gextra[i] = min_ptr->fextra[i];
gg = min_ptr->fnorm_sqr();
/* --------------------------------------------------------------------------
*/
// initialization for DOFs in constrained geometry
Real *f = lammps_main_object->atom->f[0];
for (i = 0; i < this->atom->nlocal; i++) {
min_ptr->h[3 * i + 0] = min_ptr->g[3 * i + 0] = f[3 * i + 0];
min_ptr->h[3 * i + 1] = min_ptr->g[3 * i + 1] = f[3 * i + 1];
min_ptr->h[3 * i + 2] = min_ptr->g[3 * i + 2] = f[3 * i + 2];
double x = lammps_main_object->atom->x0[i][0];
double y = lammps_main_object->atom->x0[i][1];
double z = lammps_main_object->atom->x0[i][2];
auto *geom_ptr = &this->getGeomConstrained();
if (geom_ptr && this->getGeomConstrained().contains(x, y, z)) {
lammps_main_object->atom->f[i][0] = 0;
lammps_main_object->atom->f[i][1] = 0;
lammps_main_object->atom->f[i][2] = 0;
}
}
/* --------------------------------------------------------------------------
*/
// activate the migration manager
this->ref_manager->setState(true);
#ifdef TRACE_ATOM
// for (int i = 0;
// i < this->atom->nlocal;
// ++i)
// {
// Real X[3] = {this->atom->x0[i][0],
// this->atom->x0[i][1],
// this->atom->x0[i][2] };
// SET_INTERNAL_TRACE_INDEX(X,i);
// IF_TRACED(X,"traced atom has been detected at position " <<
// internal_index);
// }
VIEW_ATOM(RefLammps<Dim>);
#endif
}
template <UInt Dim>
DomainLammpsMinimize<Dim>::DomainLammpsMinimize(LMID ID, CommGroup &group)
: LMObject(ID), DomainLammps<Dim>(ID, group) {}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammpsMinimize<Dim>::performStep1() {
VIEW_ATOM(RefLammps<Dim>);
DomainLammps<Dim>::performStep1();
VIEW_ATOM(RefLammps<Dim>);
///////////////////////////////////////////////////////
this->update->ntimestep++;
this->update->eflag_global++;
this->update->vflag_global++;
this->update->integrate->ev_set(this->update->ntimestep);
// this->modify->initial_integrate(0);
this->echangeAtomes();
LM_TOIMPLEMENT;
// this->atoms.incRelease();
///////////////////////////////////////////////////////
LAMMPS_NS::MinCG *min_ptr =
static_cast<LAMMPS_NS::MinCG *>(this->update->minimize);
int nlimit = static_cast<int>(std::min(MAXSMALLINT, (int)min_ptr->ndoftotal));
int i, m, n;
UInt fail;
UInt ntimestep;
Real beta, dot[2], dotall[2];
Real *fatom, *gatom, *hatom;
Real *fvec = lammps_main_object->atom->f[0];
ntimestep = this->update->ntimestep;
min_ptr->niter++;
// line minimization along direction h from current atom->x
min_ptr->eprevious = min_ptr->ecurrent;
min_ptr->neval = 0;
// line minimization along direction h from current atom->x
if (lm_my_proc_id == 0) {
DUMP("previous energy was: " << min_ptr->eprevious
<< ", current energy is: "
<< min_ptr->ecurrent,
DBG_MESSAGE);
}
min_ptr->eprevious = min_ptr->ecurrent;
Geometry &geom_constrained = this->getGeomConstrained();
for (int j = 0; j < min_ptr->nvec / 3; ++j) {
Real *X = lammps_main_object->atom->x0[j];
if (geom_constrained.template contains<Dim>(X)) {
min_ptr->g[3 * j + 0] = 0;
min_ptr->g[3 * j + 1] = 0;
min_ptr->g[3 * j + 2] = 0;
min_ptr->fvec[3 * j + 0] = 0;
min_ptr->fvec[3 * j + 1] = 0;
min_ptr->fvec[3 * j + 2] = 0;
min_ptr->h[3 * j] = 0;
min_ptr->h[3 * j + 1] = 0;
min_ptr->h[3 * j + 2] = 0;
}
IF_TRACED(X, "L'index est " << j << "h = ");
}
LAMMPS_NS::MinLineSearch::FnPtr tmp = min_ptr->linemin;
fail = (min_ptr->*(tmp))(min_ptr->ecurrent, min_ptr->alpha_final);
if (fail)
DUMP("linescan FAIL", DBG_MESSAGE);
// function evaluation criterion
// if (min_ptr->neval >= lammps_main_object->update->max_eval)
// LM_FATAL("min_ptr->neval " << min_ptr->neval << " ; MAXEVAL " <<
// lammps_main_object->update->max_eval);
// energy tolerance criterion
if (fabs(min_ptr->ecurrent - min_ptr->eprevious) <=
lammps_main_object->update->etol * 0.5 *
(fabs(min_ptr->ecurrent) + fabs(min_ptr->eprevious) + EPS_ENERGY))
// Real de = min_ptr->ecurrent-min_ptr->eprevious;
// Real tol = lammps_main_object->update->etol;
// Real E = .5*fabs(min_ptr->ecurrent + min_ptr->eprevious)+EPS_ENERGY;
// //check whether de is a small negative number
// if (de > -tol*E && de <= 0)
{
Real deltaEtol = min_ptr->ecurrent - min_ptr->eprevious;
DUMPBYPROC("converge in ETOL: deltaEtol = "
<< fabs(deltaEtol) << "("
<< fabs(deltaEtol) /
(0.5 * (fabs(min_ptr->ecurrent) +
fabs(min_ptr->eprevious) + EPS_ENERGY))
<< ")",
DBG_MESSAGE, 0);
current_step = nb_step_next_event - 1;
DUMP("move to next_event step" << nb_step_next_event, DBG_MESSAGE);
if (min_ptr->niter != 1) {
this->update->minimize->setup_minimal(1);
}
}
// force tolerance criterion
fvec = lammps_main_object->atom->f[0];
dot[0] = dot[1] = 0.0;
for (i = 0; i < min_ptr->nvec; i++) {
dot[0] += fvec[i] * fvec[i];
dot[1] += fvec[i] * min_ptr->g[i];
}
if (min_ptr->nextra_atom)
for (m = 0; m < min_ptr->nextra_atom; m++) {
fatom = min_ptr->fextra_atom[m];
gatom = min_ptr->gextra_atom[m];
n = min_ptr->extra_nlen[m];
for (i = 0; i < n; i++) {
dot[0] += fatom[i] * fatom[i];
dot[1] += fatom[i] * gatom[i];
}
}
MPI_Allreduce(dot, dotall, 2, MPI_DOUBLE, MPI_SUM, lammps_main_object->world);
if (min_ptr->nextra_global)
for (i = 0; i < min_ptr->nextra_global; i++) {
dotall[0] += min_ptr->fextra[i] * min_ptr->fextra[i];
dotall[1] += min_ptr->fextra[i] * min_ptr->gextra[i];
}
if (dotall[0] <
lammps_main_object->update->ftol * lammps_main_object->update->ftol) {
DUMP("converge in FTOL: " << sqrt(dotall[0]) << " < "
<< lammps_main_object->update->ftol,
DBG_MESSAGE);
current_step = nb_step_next_event - 1;
DUMP("move to next_event step" << nb_step_next_event, DBG_MESSAGE);
this->update->minimize->setup_minimal(0);
}
if (lm_my_proc_id == 0) {
DUMP("Epot = " << min_ptr->ecurrent << ", delta Epot = "
<< min_ptr->ecurrent - min_ptr->eprevious << " (tol is "
<< lammps_main_object->update->etol * 0.5 *
(fabs(min_ptr->ecurrent) + fabs(min_ptr->eprevious) +
EPS_ENERGY)
<< "), |force| = " << sqrt(dotall[0]) << " (tol is "
<< lammps_main_object->update->ftol << ")",
DBG_MESSAGE);
}
// update h from new f = -Grad(x) and old g
// beta = dotall[0]/gg would be Fletcher-Reeves CG
beta = std::max(0.0, (dotall[0] - dotall[1]) / gg);
if ((min_ptr->niter + 1) % nlimit == 0)
beta = 0.0;
gg = dotall[0];
for (i = 0; i < min_ptr->nvec; i++) {
min_ptr->g[i] = fvec[i];
min_ptr->h[i] = min_ptr->g[i] + beta * min_ptr->h[i];
}
if (min_ptr->nextra_atom)
for (m = 0; m < min_ptr->nextra_atom; m++) {
fatom = min_ptr->fextra_atom[m];
gatom = min_ptr->gextra_atom[m];
hatom = min_ptr->hextra_atom[m];
n = min_ptr->extra_nlen[m];
for (i = 0; i < n; i++) {
gatom[i] = fatom[i];
hatom[i] = gatom[i] + beta * hatom[i];
}
}
if (min_ptr->nextra_global)
for (i = 0; i < min_ptr->nextra_global; i++) {
min_ptr->gextra[i] = min_ptr->fextra[i];
min_ptr->hextra[i] = min_ptr->gextra[i] + beta * min_ptr->hextra[i];
}
// reinitialize CG if new search direction h is not downhill
dot[0] = 0.0;
for (i = 0; i < min_ptr->nvec; i++)
dot[0] += min_ptr->g[i] * min_ptr->h[i];
if (min_ptr->nextra_atom)
for (m = 0; m < min_ptr->nextra_atom; m++) {
gatom = min_ptr->gextra_atom[m];
hatom = min_ptr->hextra_atom[m];
n = min_ptr->extra_nlen[m];
for (i = 0; i < n; i++)
dot[0] += gatom[i] * hatom[i];
}
MPI_Allreduce(dot, dotall, 1, MPI_DOUBLE, MPI_SUM, this->world);
if (min_ptr->nextra_global)
for (i = 0; i < min_ptr->nextra_global; i++)
dotall[0] += min_ptr->gextra[i] * min_ptr->hextra[i];
if (dotall[0] <= 0.0) {
for (i = 0; i < min_ptr->nvec; i++)
min_ptr->h[i] = min_ptr->g[i];
if (min_ptr->nextra_atom)
for (m = 0; m < min_ptr->nextra_atom; m++) {
gatom = min_ptr->gextra_atom[m];
hatom = min_ptr->hextra_atom[m];
n = min_ptr->extra_nlen[m];
for (i = 0; i < n; i++)
hatom[i] = gatom[i];
}
if (min_ptr->nextra_global)
for (i = 0; i < min_ptr->nextra_global; i++)
min_ptr->hextra[i] = min_ptr->gextra[i];
}
// output for thermo, dump, restart files
if (lammps_main_object->output->next == ntimestep) {
lammps_main_object->timer->stamp();
lammps_main_object->output->write(ntimestep);
lammps_main_object->timer->stamp(TIME_OUTPUT);
}
lammps_main_object->output->write(ntimestep);
VIEW_ATOM(RefLammps<Dim>);
if (current_step == nb_step_next_event - 2) {
// it check 2 step earlier than given relax time, because
// stimulation_next_event is visited one step earlier than the input value,
// and
// current time step is setted one step before of nb_step_next_event
LM_FATAL("Not converged within given relaxation time\n");
}
if (current_step == nb_step_next_event - 2) {
LM_FATAL("Not converged within given relaxation time\n");
}
// if (lm_my_proc_id == 0)
// DUMP("E = " << FORMATREAL(min_ptr->ecurrent) << " FNORM = " <<
// dotall[0],DBG_MESSAGE);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammpsMinimize<Dim>::performStep3() {}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammpsMinimize<Dim>::performStep2() {
if (!this->getCommGroup().amIinGroup())
return;
LAMMPS_NS::MinCG *min_ptr =
static_cast<LAMMPS_NS::MinCG *>(this->update->minimize);
VIEW_ATOM(RefLammps<Dim>);
STARTTIMER("AtomStepForces");
this->echangeAtomes();
// //This is ugly fix to start minimization from restart file
// static bool restart = false;
// if (restart) {
min_ptr->reset_vectors();
// }
// restart = true;
UInt eflag, vflag;
eflag = 1;
vflag = 0; // virial_every;
if (LAMMPS_NS::LAMMPS::output->next_thermo == this->update->ntimestep) {
eflag = 1;
vflag = 0; // virial_thermo;
}
this->force_clear(vflag);
this->timer->stamp();
if (this->atom->molecular) {
if (this->force->bond)
this->force->bond->compute(eflag, vflag);
if (this->force->angle)
this->force->angle->compute(eflag, vflag);
if (this->force->dihedral)
this->force->dihedral->compute(eflag, vflag);
if (this->force->improper)
this->force->improper->compute(eflag, vflag);
this->timer->stamp(TIME_BOND);
}
if (this->force->pair) {
this->force->pair->compute(eflag, vflag);
this->timer->stamp(TIME_PAIR);
}
if (this->force->kspace) {
this->force->kspace->compute(eflag, vflag);
this->timer->stamp(TIME_KSPACE);
}
if (this->force->newton) {
this->comm->reverse_comm();
this->timer->stamp(TIME_COMM);
}
if (this->modify->n_post_force)
this->modify->post_force(vflag);
// LammpsSys::output->thermo->compute(1);
Geometry &geom_constrained = this->getGeomConstrained();
for (UInt j = 0; j < static_cast<UInt>(min_ptr->nvec) / 3; ++j) {
Real x = this->atom->x0[j][0];
Real y = this->atom->x0[j][1];
Real z = this->atom->x0[j][2];
if (geom_constrained.contains(x, y, z)) {
this->atom->f[j][0] = 0;
this->atom->f[j][1] = 0;
this->atom->f[j][2] = 0;
}
}
STOPTIMER("AtomStepForces");
VIEW_ATOM(RefLammps<Dim>);
}
/* -------------------------------------------------------------------------- */
template class DomainLammpsMinimize<2>;
template class DomainLammpsMinimize<3>;
__END_LIBMULTISCALE__

Event Timeline

c4science · Help