Page MenuHomec4science
Files F91864903

domain_lammps_minimize.cc
No OneTemporary
Actions

File Metadata

Created
Fri, Nov 15, 06:46

domain_lammps_minimize.cc
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 "lm_common.hh"
#include "lammps_common.hh"
#include "domain_lammps_minimize.hh"
#include "import_lammps.hh"
#include "trace_atom.hh"
#include "min_cg.h"
#include "communicator.hh"
/* -------------------------------------------------------------------------- */
// LAMMPS include files
#include "integrate.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.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];
if (&this->getGeomConstrained() && 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->atoms.getRefManager().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(DomainID ID,CommGroup GID):
DomainLammps<Dim>(ID, GID) {
}
/* -------------------------------------------------------------------------- */
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();
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 (!Communicator::getCommunicator().amIinGroup(this->getGroupID())) 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