Diffusion LibMultiScale (features/coupler_component) Browse

coupler_compon…

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__

libmultiscale/src/md/lammps/domain_lammps_minimize.cc-old301d6de84779features/coupler_compon…

domain_lammps_minimize.cc-oldView Options

libmultiscale/src/md/lammps/domain_lammps_minimize.cc-old
301d6de84779features/coupler_compon…

domain_lammps_minimize.cc-old
View Options