Diffusion LibMultiScale (features/cadd) Browse

cadd

domain_lammps.cc
View Options

			/**
			* @file domain_lammps.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 model wrapping LAMMPS
			*
			* @section LICENSE
			*
			* Copyright INRIA and CEA
			*
			* The LibMultiScale is a C++ parallel framework for the multiscale
			* coupling methods dedicated to material simulations. This framework
			* provides an API which makes it possible to program coupled simulations
			* and integration of already existing codes.
			*
			* This Project was initiated in a collaboration between INRIA Futurs Bordeaux
			* within ScAlApplix team and CEA/DPTA Ile de France.
			* The project is now continued at the Ecole Polytechnique Fédérale de Lausanne
			* within the LSMS/ENAC laboratory.
			*
			* This software is governed by the CeCILL-C license under French law and
			* abiding by the rules of distribution of free software. You can use,
			* modify and/ or redistribute the software under the terms of the CeCILL-C
			* license as circulated by CEA, CNRS and INRIA at the following URL
			* "http://www.cecill.info".
			*
			* As a counterpart to the access to the source code and rights to copy,
			* modify and redistribute granted by the license, users are provided only
			* with a limited warranty and the software's author, the holder of the
			* economic rights, and the successive licensors have only limited
			* liability.
			*
			* In this respect, the user's attention is drawn to the risks associated
			* with loading, using, modifying and/or developing or reproducing the
			* software by the user in light of its specific status of free software,
			* that may mean that it is complicated to manipulate, and that also
			* therefore means that it is reserved for developers and experienced
			* professionals having in-depth computer knowledge. Users are therefore
			* encouraged to load and test the software's suitability as regards their
			* requirements in conditions enabling the security of their systems and/or
			* data to be ensured and, more generally, to use and operate it in the
			* same conditions as regards security.
			*
			* The fact that you are presently reading this means that you have had
			* knowledge of the CeCILL-C license and that you accept its terms.
			*
			*/

			/* -------------------------------------------------------------------------- */
			//#define CHECK_STABILITY
			/* -------------------------------------------------------------------------- */
			#include "domain_lammps.hh"
			#include "factory_multiscale.hh"
			#include "import_lammps.hh"
			#include "lammps_common.hh"
			#include "lm_common.hh"
			#include "reference_manager_lammps.hh"
			#include "trace_atom.hh"
			/* -------------------------------------------------------------------------- */
			// LAMMPS includes
			#include <angle.h>
			#include <bond.h>
			#include <change_box.h>
			#include <dihedral.h>
			#include <displace_box.h>
			#include <fix_minimize.h>
			#include <improper.h>
			#include <integrate.h>
			#include <kspace.h>
			#include <min_cg.h>
			#include <pair.h>
			#include <variable.h>
			#include <verlet.h>
			/* -------------------------------------------------------------------------- */
			#include <functional>
			#include <iomanip>
			#include <iostream>
			#include <iterator>
			/* -------------------------------------------------------------------------- */
			#include "compute_lammps.hh"
			#include "geometry_manager.hh"
			#include "lm_communicator.hh"
			#include "lm_parser.hh"
			/* -------------------------------------------------------------------------- */
			LAMMPS_NS::LAMMPS *lammps_main_object = NULL;
			/* -------------------------------------------------------------------------- */

			__BEGIN_LIBMULTISCALE__

			template <UInt Dim> void DomainLammps<Dim>::init() {
			if (!this->getCommGroup().amIinGroup())
			return;

			this->initModel();
			this->update->whichflag = 1;
			LAMMPS_NS::LAMMPS::init();
			this->update->integrate->setup();
			this->initDOFs();
			this->ref_manager->setState(true);
			this->atoms.getBoundingBox() = this->getDomainBoundingBox();
			}

			/* -------------------------------------------------------------------------- */
			template <UInt Dim> ArrayView DomainLammps<Dim>::gradient() {
			return getField(_force);
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> ArrayView DomainLammps<Dim>::primal() {
			return getField(_position);
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim>
			auto DomainLammps<Dim>::primalTimeDerivative() -> ArrayView {
			return getField(_velocity);
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> auto DomainLammps<Dim>::mass() -> ArrayView {
			return getField(_mass);
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> auto DomainLammps<Dim>::acceleration() -> ArrayView {
			return getField(_acceleration);
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim>
			auto DomainLammps<Dim>::getField(FieldType field_type) -> ArrayView {
			switch (field_type) {
			case _position0:
			return ArrayView(&lammps_main_object->atom->x0[0][0],
			lammps_main_object->atom->nlocal, 3);
			case _velocity:
			return ArrayView(&lammps_main_object->atom->v[0][0],
			lammps_main_object->atom->nlocal, 3);
			case _position:
			return ArrayView(&lammps_main_object->atom->x[0][0],
			lammps_main_object->atom->nlocal, 3);
			case _force:
			return ArrayView(&lammps_main_object->atom->f[0][0],
			lammps_main_object->atom->nlocal, 3);
			case _mass:
			return ArrayView(&lammps_main_object->atom->mass[0],
			lammps_main_object->atom->ntypes, 1);
			case _acceleration:
			return ArrayView(this->acceleration_field.data(),
			lammps_main_object->atom->nlocal, 3);
			default:
			LM_FATAL("Not accessible field: " << field_type);
			}
			}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::loadLammpsConfigFile() {

			// forward LM variables to lammps

			std::map<std::string, double> &vars = Parser::getAlgebraicVariables();
			std::map<std::string, double>::iterator it = vars.begin();
			std::map<std::string, double>::iterator end = vars.end();
			while (it != end) {
			std::string varname = it->first;
			std::stringstream varval;
			varval << std::setprecision(16) << std::scientific << it->second;
			LAMMPS_NS::Variable *lmpvar =
			static_cast<LAMMPS_NS::LAMMPS *>(this)->input->variable;

			char cvarname = const_cast<char >(varname.c_str());

			char cvarval[30];
			varval.str().copy(cvarval, 29);
			cvarval[varval.str().length()] = 0;

			if (lmpvar->find(cvarname) >= 0)
			LM_FATAL("cannot export LibMultiscale variable "
			<< varname << " to LAMMPS since variable already exists");

			DUMP(cvarname << " = " << cvarval, DBG_INFO);
			lmpvar->add_equal(cvarname, cvarval);

			++it;
			}

			std::string file_to_read = lammpsfile;
			DUMP("lammpsfile"
			<< " = " << lammpsfile,
			DBG_INFO);

			if (create_header_flag && this->getCommGroup().getMyRank() == 0) {
			std::string unit_name;
			if (code_unit_system == metal_unit_system)
			unit_name = "metal";
			else if (code_unit_system == atomic_unit_system)
			unit_name = "real";
			else if (code_unit_system == real_unit_system)
			unit_name = "si";
			else
			LM_FATAL("unknown units");

			file_to_read = "./generated_lammps.config";
			std::ofstream file_lammps(file_to_read.c_str());

			file_lammps << "units " << unit_name << std::endl;
			file_lammps << "dimension " << Dim << std::endl;
			file_lammps << "atom_style atomic" << std::endl;
			file_lammps << "boundary ";
			for (UInt i = 0; i < 3; ++i) {
			file_lammps << boundaries[i] << " ";
			}
			file_lammps << std::endl;
			file_lammps << "lattice " << lattice << " " << std::setprecision(15)
			<< static_cast<Real>(lattice_size);

			file_lammps << " orient x ";
			for (UInt i = 0; i < 3; ++i)
			file_lammps << lattice_orient_x[i] << " ";
			if (Dim > 1) {
			file_lammps << "orient y ";
			for (UInt i = 0; i < 3; ++i)
			file_lammps << lattice_orient_y[i] << " ";
			}
			if (Dim == 3) {
			file_lammps << "orient z ";
			for (UInt i = 0; i < 3; ++i)
			file_lammps << lattice_orient_z[i] << " ";
			}

			file_lammps << "spacing ";
			for (UInt i = 0; i < 3; ++i)
			file_lammps << lattice_spacing[i] << " ";

			file_lammps << "origin ";
			for (UInt i = 0; i < Dim; ++i)
			file_lammps << Real(lattice_origin[i]) << " ";
			for (UInt i = Dim; i < 3; ++i)
			file_lammps << 0.0 << " ";
			file_lammps << std::endl;

			file_lammps << "region box block ";
			for (UInt i = 0; i < 6; ++i)
			file_lammps << replica[i] << " ";
			file_lammps << std::endl;
			file_lammps << "create_box 1 box" << std::endl;
			file_lammps << "create_atoms 1 box" << std::endl;
			file_lammps << std::endl << std::endl;

			if (triclinic == 1 && floorf(tilt[1] * 100 + 0.5) / 100 != 0.0) {
			// double value;
			std::string direction;
			Real tilting = floorf(tilt[1] * 100 + 0.5) / 100;
			if (tilt[0] == 0.0) {
			// value = tilt[1]/lattice_size/lattice_spacing[0];
			direction = "xy";
			} else if (tilt[0] == 1.0) {
			// value = tilt[1]/lattice_size/lattice_spacing[0];
			direction = "xz";
			} else if (tilt[0] == 2.0) {
			// value = tilt[1]/lattice_size/lattice_spacing[0];
			direction = "yz";
			} else
			LM_FATAL("tilt must be happend in xy(0), xz(1) or yz(2) axis!!!");

			file_lammps << "change_box triclinic" << std::endl;
			file_lammps << "displace_box all " << direction << " final " << tilting
			<< " remap none units box" << std::endl;
			// file_lammps << "displace_box all "<<direction<<" final "<<value<<"
			// remap none"<< std::endl;
			}

			std::ifstream in(lammpsfile.c_str());
			std::istreambuf_iterator<char> src(in.rdbuf());
			std::istreambuf_iterator<char> end;
			std::ostream_iterator<char> dest(file_lammps);
			std::copy(src, end, dest);
			}
			DUMP("opening lammps config file " << file_to_read, DBG_INFO_STARTUP);
			infile = fopen(file_to_read.c_str(), "rb");
			if (infile == NULL)
			LM_FATAL("cannot open " << file_to_read << " as a lammps config file !!!");
			input->file();
			DUMP("lammps file loaded", DBG_INFO);

			// register lammps computes to libmultiscale system
			UInt ncompute = this->modify->ncompute;

			for (UInt i = 0; i < ncompute; ++i) {
			auto _ptr = std::make_shared<ComputeLammps<LAMMPS_NS::Compute>>(
			(this->modify->compute[i]), this);
			_ptr->setCommGroup(this->getCommGroup());
			FilterManager::getManager().addObject(_ptr);
			}

			// register lammps computes to libmultiscale system
			UInt nfix = this->modify->nfix;

			for (UInt i = 0; i < nfix; ++i) {
			auto _ptr = std::make_shared<ComputeLammps<LAMMPS_NS::Fix>>(
			(this->modify->fix[i]), this);
			_ptr->setCommGroup(this->getCommGroup());
			FilterManager::getManager().addObject(_ptr);
			}
			}

			/* -------------------------------------------------------------------------- */
			template <UInt Dim> void DomainLammps<Dim>::changeBox() {
			auto &geom_manager = GeometryManager::getManager();
			auto c = geom_manager.getGeometry(this->newGeomBox)->getBoundingBox();
			this->setDomainBoundingBox(c);
			}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::initModel() {

			lammps_main_object = this;

			ref_manager = std::make_shared<ReferenceManagerLammps<Dim>>(this->atoms);
			this->atoms.setRefManager(ref_manager);
			ref_manager->setCommGroup(this->getCommGroup());

			ImportLammpsInterface::createImportLammps(ref_manager);

			if (this->getGeom().getDim() == 1)
			LM_FATAL("LAMMPS cannot work in 1D");

			ImportLammpsInterface::setGeom(&this->getGeom());

			// for (auto &&_type : geom_by_type) {
			// Geometry &g = GeometryManager::getManager().getObject(_type.second);
			// UInt itype = _type.first;
			// ImportLammpsInterface::setGeom(&g, itype);
			// }

			loadLammpsConfigFile();

			this->atoms.init();

			update->whichflag = 0;
			update->firststep = 0;
			update->laststep = nb_step;
			update->ntimestep = 0;
			update->beginstep = 0;
			update->endstep = update->laststep;

			if (flag_units) {
			force->ftm2v = code_unit_system.ft_m2v;
			// force->nktv2p = 1;
			force->mvv2e = code_unit_system.mvv2e;
			}

			if (this->domain->dimension != Dim)
			LM_FATAL("mismatch dimension: check lammps config file");

			if (std::string(this->atom->atom_style) != "atomic")
			LM_FATAL("can only work with atomic style");

			this->acceleration_field.resize(atom->nlocal, 3);
			this->acceleration_field.setZero();
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> void DomainLammps<Dim>::initDOFs() {
			{
			UInt nall;
			if (force->newton)
			nall = atom->nlocal + atom->nghost;
			else
			nall = atom->nlocal;

			Real **x = atom->x;
			Real **x0 = atom->x0;

			for (UInt i = 0; i < nall; i++) {
			x0[i][0] = x[i][0];
			x0[i][1] = x[i][1];
			x0[i][2] = x[i][2];
			}
			}
			if (this->shouldRestart()) {
			this->restart_start = true;
			this->readXMLFile(this->restartfile);

			if (this->newGeomBox != invalidGeom &&
			current_step >= this->when_change_box) {
			this->changeBox();
			this->enforceCompatibility(); // perform migration after change box to
			// keep periodicity
			this->when_change_box = current_step + nb_step; // to avoid repetitions
			} else if (neighbor->decide()) {
			this->enforceCompatibility();
			this->resetBox();
			}
			this->updateGradient();

			LAMMPS_NS::LAMMPS::output->setup(false);
			}
			// conversion of masses from g/mol to Kg
			if (flag_units) {
			UnitsConverter uc;
			uc.setOutUnits(code_unit_system);
			if (strcmp(update->unit_style, "real") == 0)
			uc.setInUnits(atomic_unit_system);
			else if (strcmp(update->unit_style, "metal") == 0)
			uc.setInUnits(metal_unit_system);
			else
			LM_FATAL("not known lammps unit system:"
			<< " things need to be done to integrate");

			uc.computeConversions();

			for (int i = 1; i <= atom->ntypes; i++)
			atom->mass[i] *= uc.getConversion<Mass>();
			}

			#ifdef TRACE_ATOM
			for (int i = 0; i < atom->nlocal; ++i) {
			Real X[3] = {atom->x0[i][0], atom->x0[i][1], 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> Real DomainLammps<Dim>::potentialEnergy() {
			return update->minimize->pe_compute->compute_scalar();
			}

			/* -------------------------------------------------------------------------- */
			extern "C" {
			static UInt fargc = 1;
			static char prog[20] = "lammps";
			static char *___tmp = &prog[0];
			static char **fargv = &___tmp;
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> DomainLammps<Dim>::~DomainLammps() {}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim>
			DomainLammps<Dim>::DomainLammps(const LMID &ID, CommGroup &group)
			: LMObject(ID), DomainAtomic<ContainerLammps<Dim>, Dim>(ID, group),
			LAMMPS(fargc, fargv, group.getMPIComm()), flag_units(true), triclinic(0) {

			newGeomBox = invalidGeom;
			when_change_box = 0;
			lattice = "";
			lattice_size = 0.;
			for (UInt i = 0; i < 6; ++i)
			replica[i] = 0;
			create_header_flag = false;
			restart_start = false;

			for (UInt i = 0; i < 3; ++i) {
			lattice_orient_x[i] = 0;
			lattice_orient_y[i] = 0;
			lattice_orient_z[i] = 0;
			lattice_origin[i] = 0.;
			}
			lattice_orient_x[0] = 1;
			lattice_orient_y[1] = 1;
			lattice_orient_z[2] = 1;

			lattice_spacing[0] = 1.0;
			lattice_spacing[1] = 1.0;
			lattice_spacing[2] = 1.0;
			}

			/* -------------------------------------------------------------------------- */
			template <UInt Dim> Cube DomainLammps<Dim>::getSubDomainBoundingBox() {
			Cube c;
			c.setXmin(VectorView<3>(domain->sublo));
			c.setXmax(VectorView<3>(domain->subhi));
			return c;
			}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim>
			void DomainLammps<Dim>::setSubDomainBoundingBox(const Cube &cube) {

			VectorView<3>(domain->sublo) = cube.getXmin();
			VectorView<3>(domain->subhi) = cube.getXmax();

			if (neighbor->decide())
			this->resetBox();
			}

			/* -------------------------------------------------------------------------- */
			template <UInt Dim> Cube DomainLammps<Dim>::getDomainBoundingBox() {
			Cube c;
			c.setXmin(VectorView<3>(domain->boxlo));
			c.setXmax(VectorView<3>(domain->boxhi));
			return c;
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim>
			void DomainLammps<Dim>::setDomainBoundingBox(const Cube &cube) {
			VectorView<3>(domain->boxlo) = cube.getXmin();
			VectorView<3>(domain->boxhi) = cube.getXmax();

			if (neighbor->decide())
			this->resetBox();
			}
			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::resetBox() {
			atom->setup();
			if (domain->triclinic)
			domain->x2lamda(atom->nlocal);
			domain->pbc();
			domain->reset_box();
			comm->setup();
			if (neighbor->style)
			neighbor->setup_bins();
			comm->exchange();
			if (atom->sortfreq > 0)
			atom->sort();
			comm->borders();
			if (domain->triclinic)
			domain->lamda2x(atom->nlocal + atom->nghost);
			neighbor->build();
			}
			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::enforceCompatibility() {
			if (this->getCommGroup().amIinGroup()) {
			this->echangeAtomes();
			this->acceleration_field.resize(atom->nlocal, 3);
			}
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> void DomainLammps<Dim>::echangeAtomes() {
			VIEW_ATOM(RefLammps<Dim>);
			UInt nflag = neighbor->decide();

			// if (nflag == 0 && this->restart_start==false && this->update_pad == false)
			// {
			if (nflag == 0 && this->restart_start == false) {
			timer->stamp();
			comm->forward_comm();
			timer->stamp(TIME_COMM);
			} else {
			if (modify->n_pre_exchange)
			modify->pre_exchange();
			if (domain->triclinic)
			domain->x2lamda(atom->nlocal);

			domain->pbc();
			if (domain->box_change) {
			domain->reset_box();
			comm->setup();
			if (neighbor->style)
			neighbor->setup_bins();
			}
			timer->stamp();

			comm->exchange();
			comm->borders();
			if (domain->triclinic)
			domain->lamda2x(atom->nlocal + atom->nghost);
			timer->stamp(TIME_COMM);
			if (modify->n_pre_neighbor)
			modify->pre_neighbor();
			neighbor->build();
			timer->stamp(TIME_NEIGHBOR);
			this->restart_start = false;
			}

			if (this->when_change_box == current_time) {
			// if one change box size, this part has to be fulfilled.
			if (modify->n_pre_exchange)
			modify->pre_exchange();
			if (domain->triclinic)
			domain->x2lamda(atom->nlocal);
			domain->pbc();
			if (domain->box_change) {
			domain->reset_box();
			comm->setup();
			if (neighbor->style)
			neighbor->setup_bins();
			}
			timer->stamp();
			comm->exchange();
			comm->borders();
			if (domain->triclinic)
			domain->lamda2x(atom->nlocal + atom->nghost);
			timer->stamp(TIME_COMM);
			if (modify->n_pre_neighbor)
			modify->pre_neighbor();
			neighbor->build();
			timer->stamp(TIME_NEIGHBOR);
			this->restart_start = false;
			}
			VIEW_ATOM(RefLammps<Dim>);
			}

			/* -------------------------------------------------------------------------- */
			template <UInt Dim> bool DomainLammps<Dim>::isPeriodic(UInt dir) {
			switch (dir) {
			case 0:
			return this->domain->xperiodic;
			break;
			case 1:
			return this->domain->yperiodic;
			break;
			case 2:
			return this->domain->zperiodic;
			break;
			default:
			LM_FATAL("called on invalid-dimension");
			return false;
			}
			}

			/* -------------------------------------------------------------------------- */
			template <UInt Dim> void DomainLammps<Dim>::force_clear(UInt) {
			UInt i;

			// clear global force array
			// nall includes ghosts only if either newton flag is set

			UInt nall;
			if (force->newton)
			nall = atom->nlocal + atom->nghost;
			else
			nall = atom->nlocal;

			Real **f = atom->f;
			for (i = 0; i < nall; i++) {
			f[i][0] = 0.0;
			f[i][1] = 0.0;
			f[i][2] = 0.0;
			}

			// if (integ.torqueflag) {
			// Real **torque = atom->torque;
			// for (i = 0; i < nall; i++) {
			// torque[i][0] = 0.0;
			// torque[i][1] = 0.0;
			// torque[i][2] = 0.0;
			// }
			// }

			// if (granflag) {
			// Real **phia = atom->phia;
			// for (i = 0; i < nall; i++) {
			// phia[i][0] = 0.0;
			// phia[i][1] = 0.0;
			// phia[i][2] = 0.0;
			// }
			// }

			// clear pair force array if using it this timestep to compute virial

			// if (vflag == 2 && pairflag) {
			// if (atom->nmax > maxpair) {
			// maxpair = atom->nmax;
			// memory->destroy_2d_Real_array(f_pair);
			// f_pair = memory->create_2d_Real_array(maxpair,3,"verlet:f_pair");
			// }
			// for (i = 0; i < nall; i++) {
			// f_pair[i][0] = 0.0;
			// f_pair[i][1] = 0.0;
			// f_pair[i][2] = 0.0;
			// }
			// }
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> void DomainLammps<Dim>::setCurrentStep(UInt ts) {
			current_step = ts;

			// content of the Update::reset_timestep function
			update->eflag_global = update->vflag_global = ts;
			for (int i = 0; i < modify->ncompute; i++) {
			modify->compute[i]->invoked_scalar = -1;
			modify->compute[i]->invoked_vector = -1;
			modify->compute[i]->invoked_array = -1;
			modify->compute[i]->invoked_peratom = -1;
			modify->compute[i]->invoked_local = -1;
			// modify->compute[i]->eflag_global = -1;
			}
			update->ntimestep = ts;

			if (update->laststep < ts)
			update->laststep += ts;
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> void DomainLammps<Dim>::setTimeStep(Real ts) {
			update->dt = ts;
			for (int i = 0; i < modify->nfix; ++i)
			modify->fix[i]->init();
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> UInt DomainLammps<Dim>::getCurrentStep() {
			return update->ntimestep;
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> Real DomainLammps<Dim>::getTimeStep() { return update->dt; }
			/* -------------------------------------------------------------------------- */

			template <UInt Dim>
			int DomainLammps<Dim>::getOrientValue(UInt axis, UInt comp) {
			if (axis == 0) {
			return this->lattice_orient_x[comp];
			} else if (axis == 1) {
			return this->lattice_orient_y[comp];
			} else if (axis == 2) {
			return this->lattice_orient_z[comp];
			} else {
			LM_FATAL("wrong axis value: " << axis);
			}
			}
			/* -------------------------------------------------------------------------- */
			template <UInt Dim> Real DomainLammps<Dim>::getLatticeConstant() {
			return this->lattice_size;
			}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::newmarkPredictor() {
			if (this->newGeomBox != invalidGeom &&
			current_step >= this->when_change_box) {
			this->changeBox();
			}

			VIEW_ATOM(RefLammps<Dim>);
			STARTTIMER("AtomStepPosition");
			this->update->ntimestep++;
			this->update->integrate->ev_set(this->update->ntimestep);
			this->modify->initial_integrate(0);
			this->echangeAtomes();
			STOPTIMER("AtomStepPosition");
			VIEW_ATOM(RefLammps<Dim>);
			this->getOutput()->incRelease();
			}

			template <UInt Dim> void DomainLammps<Dim>::newmarkCorrector() {
			VIEW_ATOM(RefLammps<Dim>);
			DUMP("AtomStepVelocities x " << this->atom->x, DBG_INFO);
			STARTTIMER("AtomStepVelocities");
			this->modify->final_integrate();
			if (this->modify->n_end_of_step)
			this->modify->end_of_step();

			if (LAMMPS_NS::LAMMPS::output->next == this->update->ntimestep) {
			this->timer->stamp();
			LAMMPS_NS::LAMMPS::output->write(this->update->ntimestep);
			this->timer->stamp(TIME_OUTPUT);
			}
			STOPTIMER("AtomStepVelocities");
			VIEW_ATOM(RefLammps<Dim>);
			this->getOutput()->incRelease();
			}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::updateAcceleration() {
			if (!this->getCommGroup().amIinGroup())
			return;

			double *rmass = atom->rmass;
			double *mass = atom->mass;
			int *type = atom->type;
			auto f = this->gradient();

			std::function<Real(UInt)> compute_mass_inverse;

			if (rmass)
			compute_mass_inverse = [&](UInt i) {
			return code_unit_system.ft_m2v / rmass[i];
			};
			else
			compute_mass_inverse = [&](UInt i) {
			return code_unit_system.ft_m2v / mass[type[i]];
			};

			for (UInt i = 0; i < UInt(lammps_main_object->atom->nlocal); ++i) {
			Real mass_inverse = compute_mass_inverse(i);
			this->acceleration_field(i, 0) = mass_inverse * f(i, 0);
			this->acceleration_field(i, 1) = mass_inverse * f(i, 1);
			this->acceleration_field(i, 2) = mass_inverse * f(i, 2);
			}
			}

			/* -------------------------------------------------------------------------- */

			template <UInt Dim> void DomainLammps<Dim>::updateGradient() {
			if (!this->getCommGroup().amIinGroup())
			return;
			VIEW_ATOM(RefLammps<Dim>);
			STARTTIMER("AtomStepForces");
			UInt eflag, vflag;
			eflag = this->update->integrate->eflag;
			vflag = this->update->integrate->vflag;

			this->force_clear(vflag);
			if (this->modify->n_pre_force)
			this->modify->pre_force(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);

			STOPTIMER("AtomStepForces");
			VIEW_ATOM(RefLammps<Dim>);

			this->getOutput()->incRelease();
			}

			/* -------------------------------------------------------------------------- */
			/* LMDESC LAMMPS_BASE
			This plugin make the interface with lammps
			*/
			/* LMHERITANCE domain_md */

			/* LM EXAMPLE
			Section MultiScale AtomsUnits\\
			...\\
			GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\
			MODEL LAMMPS md\\
			...\\
			endSection\\ \\

			Section LAMMPS:md AtomsUnits \\
			LAMMPS_FILE lammps.in \\
			DOMAIN_GEOMETRY 1 \\
			endSection
			*/

			template <UInt Dim> void DomainLammps<Dim>::declareParams() {

			DomainAtomic<ContainerLammps<Dim>, Dim>::declareParams();

			/*LMKEYWORD LAMMPS_FILE
			This specifies the configuration file to be used for lammps.
			If the create header keyword is used then only the potential
			definition should be into that file. Otherwise, a complete
			lammps configuration file can be specified.
			*/

			this->parseKeyword("LAMMPS_FILE", lammpsfile);

			/* LMKEYWORD CREATE_HEADER
			A header is automatically generated based on the information
			passed by keywords LATTICE, LATTICE\_SIZE, BOUNDARY and REPLICA
			*/
			this->parseTag("CREATE_HEADER", create_header_flag, false);

			/* LMKEYWORD LATTICE
			Specifies the lattice to be used to generate the lammps header.
			Admissible values are lammps keywords
			(see \url{http://lammps.sandia.gov/doc/lattice.html})
			*/
			this->parseKeyword("LATTICE", lattice, "");

			/* LMKEYWORD LATTICE_SIZE
			Lattice length parameter to generate the primitive lattice
			*/

			this->parseKeyword("LATTICE_SIZE", lattice_size, 0.);

			/* LMKEYWORD LATTICE_ORIGIN
			Lattice length parameter to generate the primitive lattice
			*/

			this->parseVectorKeyword("LATTICE_ORIGIN", 3, lattice_origin,
			VEC_DEFAULTS(0., 0., 0.));

			/* LMKEYWORD LATTICE_ORIENTX
			Lattice length parameter to generate the primitive lattice
			*/

			this->parseVectorKeyword("LATTICE_ORIENTX", 3, lattice_orient_x,
			VEC_DEFAULTS(1, 0, 0));

			/* LMKEYWORD LATTICE_ORIENTY
			Lattice length parameter to generate the primitive lattice
			*/

			this->parseVectorKeyword("LATTICE_ORIENTY", 3, lattice_orient_y,
			VEC_DEFAULTS(0, 1, 0));

			/* LMKEYWORD LATTICE_ORIENTZ
			Lattice length parameter to generate the primitive lattice
			*/

			this->parseVectorKeyword("LATTICE_ORIENTZ", 3, lattice_orient_z,
			VEC_DEFAULTS(0, 0, 1));

			/* LMKEYWORD SPACING
			The minimal normalized lattice sizes
			*/

			this->parseVectorKeyword("SPACING", 3, lattice_spacing,
			VEC_DEFAULTS(1., 1., 1.));

			/* LMKEYWORD BOUNDARY
			Sequence of lammps code for boundary.
			(see \url{http://lammps.sandia.gov/doc/boundary.html})
			*/
			this->parseVectorKeyword("BOUNDARY", 3, boundaries, VEC_DEFAULTS("", "", ""));

			/* LMKEYWORD REPLICA
			For the creation of atoms a region, specified as the number of
			replicas of the lattice must be provided, before calling the

			create atoms command.
			*/
			this->parseVectorKeyword("REPLICA", 6, replica,
			VEC_DEFAULTS(0, 0, 0, 0, 0, 0));

			/* LMKEYWORD CHANGE_DOMAIN_BOX
			This allows to specify live reshaping of the bounding box.
			Useful to enforce stresses states or dislocation displacement
			field with periodic boundary conditions.
			*/
			this->parseKeyword("CHANGE_DOMAIN_BOX", newGeomBox, invalidGeom);

			/* LMKEYWORD CHANGE_DOMAIN_BOX_TIME
			Gives a timestep where the reshaping provided by CHANGE\_DOMAIN\_BOX
			should occur.
			*/
			this->parseKeyword("CHANGE_DOMAIN_BOX_TIME", when_change_box, 0u);

			/* LMKEYWORD TRICLINIC
			This changes the orthogonal simulation box to triclinic box
			*/
			this->parseKeyword("TRICLINIC", triclinic, 0);

			/* LMKEYWORD TILT
			This tils the simulation box by given amount value
			0: xy-direction; 1: xz-direction; 2: yz-direction
			(TILT 0 0.25: tilt simulation box in xy-direction by 0.25)
			*/
			this->parseVectorKeyword("TILT", 2, tilt, VEC_DEFAULTS(0., 0.));

			/* LMKEYWORD UNITS_CONVERSION
			This is for debug purpose: force libmultiscale not to touch the
			units conversion factors of lammps.
			*/
			this->parseTag("UNITS_CONVERSION", flag_units, true);
			}
			/* -------------------------------------------------------------------------- */

			template class DomainLammps<2>;
			template class DomainLammps<3>;
			__END_LIBMULTISCALE__

libmultiscale/src/md/lammps/domain_lammps.cc2fb34b74795efeatures/cadd

domain_lammps.ccView Options

libmultiscale/src/md/lammps/domain_lammps.cc
2fb34b74795efeatures/cadd

domain_lammps.cc
View Options