domain_lammps.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, May 21, 13:18

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 = LAMMPS_NS::LAMMPS::input->variable;

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

	char cvarval[30];
	varval.str().copy(cvarval, 29);
	cvarval[varval.str().length()] = 0;
	kk 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->acquireContext(*this);
	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->acquireContext(*this);
	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->acquireContext(*this);

	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__

domain_lammps.ccNo OneTemporaryActions

File Metadata

domain_lammps.ccView Options

Event Timeline

domain_lammps.cc
No OneTemporary
Actions

domain_lammps.cc
View Options