domain_lammps.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, May 22, 11:58

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"
	#include <iomanip>
	#include <iostream>
	#include <iterator>
	/* -------------------------------------------------------------------------- */
	// 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 "communicator.hh"
	#include "compute_lammps.hh"
	#include "geometry_manager.hh"
	#include "lm_parser.hh"
	/* -------------------------------------------------------------------------- */
	LAMMPS_NS::LAMMPS *lammps_main_object = NULL;
	/* -------------------------------------------------------------------------- */

	__BEGIN_LIBMULTISCALE__

	template <UInt Dim> void DomainLammps<Dim>::init() { this->initModel(); }

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

	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;

	if (create_header_flag &&
	Communicator::getCommunicator().groupRank(lm_my_proc_id,
	this->getGroupID()) == 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 << 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) {
	LM_TOIMPLEMENT;

	// ComputeLammps<LAMMPS_NS::Compute> *_ptr =
	// new ComputeLammps<LAMMPS_NS::Compute>(*(this->modify->compute[i]),
	// *this);
	// _ptr->setRelease(INITIAL_MODEL_RELEASE);
	// _ptr->setCommGroup(this->getGroupID());
	// FilterManager::getManager().addObject(_ptr);
	}

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

	for (UInt i = 0; i < nfix; ++i) {
	LM_TOIMPLEMENT;
	// ComputeLammps<LAMMPS_NS::Fix> *_ptr =
	// new ComputeLammps<LAMMPS_NS::Fix>((this->modify->fix[i]), this);
	// _ptr->setRelease(INITIAL_MODEL_RELEASE);
	// _ptr->setCommGroup(this->getGroupID());
	// FilterManager::getManager().addObject(_ptr);
	}
	}

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

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

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

	template <UInt Dim> void DomainLammps<Dim>::initModel() {
	lammps_main_object = this;

	ReferenceManagerLammps<Dim> *refMgr =
	new ReferenceManagerLammps<Dim>(this->atoms);
	this->atoms.setRefManager(refMgr);
	refMgr->setMPIComm(
	Communicator::getCommunicator().getMpiGroup(this->getGroupID()));

	ImportLammpsInterface::createImportLammps(*refMgr);

	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->update_pad = false;
	}
	/* -------------------------------------------------------------------------- */
	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->performMigration(); // perform migration after change box to keep
	// periodicity
	this->when_change_box = current_step + nb_step; // to avoid repetitions
	} else if (neighbor->decide()) {
	this->performMigration();
	this->resetBox();
	}
	this->performStep2();

	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

	LM_TOIMPLEMENT;
	// this->atoms.setRelease(INITIAL_MODEL_RELEASE);
	}
	/* -------------------------------------------------------------------------- */

	template <UInt Dim> Real DomainLammps<Dim>::getEcin() { return 0.0; }

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

	template <UInt Dim> Real DomainLammps<Dim>::getEpot() {
	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;
	//#include "../../common/openmp.hh"
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> DomainLammps<Dim>::~DomainLammps() {
	if (this->atoms.hasRefManager())
	delete (&this->atoms.getRefManager());
	}

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

	template <UInt Dim>
	DomainLammps<Dim>::DomainLammps(LMID ID, CommGroup GID)
	: LMObject(ID), DomainAtomic<ContainerLammps<Dim>, Dim>(GID),
	LAMMPS(fargc, fargv, Communicator::getCommunicator().getMpiGroup(GID)),
	flag_units(true), triclinic(0) {
	// external_work = 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>
	void DomainLammps<Dim>::getSubDomainDimensions(Real xmin, Real xmax) {
	for (UInt i = 0; i < 3; ++i) {
	xmin[i] = domain->sublo[i];
	xmax[i] = domain->subhi[i];
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void DomainLammps<Dim>::getDomainDimensions(Real xmin, Real xmax) {

	for (UInt i = 0; i < 3; ++i) {
	xmin[i] = domain->boxlo[i];
	xmax[i] = domain->boxhi[i];
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void DomainLammps<Dim>::setDomainDimensions(const Vector<Dim> &xmin,
	const Vector<Dim> &xmax) {
	for (UInt i = 0; i < 3; ++i) {
	domain->boxlo[i] = xmin[i];
	domain->boxhi[i] = xmax[i];
	}
	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>::performMigration() {
	if (Communicator::getCommunicator().amIinGroup(this->getGroupID())) {
	echangeAtomes();
	}
	}
	/* -------------------------------------------------------------------------- */
	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;
	this->update_pad = 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->update_pad = false;
	this->restart_start = false;
	}
	VIEW_ATOM(RefLammps<Dim>);
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> Real DomainLammps<Dim>::getFdotDir() {
	// Real *f = atom->f[0];
	Real fdotdirme = 0.0;
	UInt n = update->minimize->nvec;
	for (UInt i = 0; i < n; i++) {
	// fdotdirme += f[i]*update->minimize->h[i];
	}
	Real fdotdirall;
	MPI_Allreduce(&fdotdirme, &fdotdirall, 1, MPI_DOUBLE, MPI_SUM,
	MPI_COMM_WORLD);
	return fdotdirall;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> Real DomainLammps<Dim>::getFMax() {
	Real fmax = 0;
	Real fmaxall;
	Real *f = atom->f[0];
	UInt n = update->minimize->nvec;
	for (UInt i = 0; i < n; i++)
	fmax = std::max(fabs(f[i]), fmax);
	MPI_Allreduce(&fmax, &fmaxall, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
	return fmaxall;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> Real DomainLammps<Dim>::getDirMax() {
	Real fmax = 0;
	Real fmaxall;
	// UInt n = update->minimize->nvec;
	// for (UInt i = 0; i < n; i++) fmax = max(fabs(update->minimize->h[i]),fmax);
	MPI_Allreduce(&fmax, &fmaxall, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
	return fmaxall;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> Real DomainLammps<Dim>::getFNorm2() {
	Real fnorm2all = 0.0;
	Real *f = atom->f[0];
	Real fnorm2me = 0.0;
	UInt n = update->minimize->nvec;
	for (UInt i = 0; i < n; i++)
	fnorm2me += f[i] * f[i];
	MPI_Allreduce(&fnorm2me, &fnorm2all, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
	return fnorm2all;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> Real DomainLammps<Dim>::getFdotOldF() {
	Real fdotoldf;
	// Real *f = atom->f[0];
	Real fdotoldfme = 0.0;
	// UInt n = update->minimize->nvec;
	// for (UInt i = 0; i < n; i++) fdotoldfme += f[i]*update->minimize->g[i];
	MPI_Allreduce(&fdotoldfme, &fdotoldf, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
	return fdotoldf;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void DomainLammps<Dim>::updateDirection(Real beta) {
	for (int i = 0; i < update->minimize->nvec; i++) {
	// update->minimize->h[i] = update->minimize->g[i] +
	// beta*update->minimize->h[i];
	// max_h = MAX(max_h,fabs(min_ptr->h[i]));
	}

	auto *geom_ptr = &this->getGeomConstrained();
	if (geom_ptr) {
	for (int j = 0; j < update->minimize->nvec / 3; ++j) {
	Real x = atom->x0[j][0];
	Real y = atom->x0[j][1];
	Real z = atom->x0[j][2];
	if (this->getGeomConstrained().contains(x, y, z)) {
	// update->minimize->h[3*j] = 0;
	// update->minimize->h[3*j+1] = 0;
	// update->minimize->h[3*j+2] = 0;
	}
	// double norm =
	// min_ptr->h[3j]min_ptr->h[3*j]
	// + min_ptr->h[3j+1]min_ptr->h[3*j+1]
	// + min_ptr->h[3j+2]min_ptr->h[3*j+2];
	// norm = sqrt(norm);
	// double factor = exp (-max_h/100/norm );
	// if (norm < max_h/100) {
	// min_ptr->h[3j] = factor;
	// min_ptr->h[3j+1] = factor;
	// min_ptr->h[3j+2] = factor;
	// }
	// if (norm < max_h/10) {
	// min_ptr->h[3*j] = 0;
	// min_ptr->h[3*j+1] = 0;
	// min_ptr->h[3*j+2] = 0;
	// }
	}
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void DomainLammps<Dim>::saveForceVector() {
	// Real *f = atom->f[0];
	for (int i = 0; i < update->minimize->nvec; i++) {
	// update->minimize->g[i] = f[i];
	}
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim>
	void DomainLammps<Dim>::displaceTowardsDirection(Real delta) {
	// Real * x = atom->x[0];
	// for (int i = 0; i < update->minimize->nvec; i++) {
	// // std::cerr << "x[" << i << "/" << update->minimize->ndof << "](" <<
	// x[i] << ") += " << delta << "*" << update->minimize->h[i] << std::endl;
	// x[i] += delta*update->minimize->h[i];
	// }
	// double local_external_work = 0;
	// double * v = lammps_main_object->atom->v[0];
	// for (int i = 0 ; i < min_ptr->ndof ; ++i){
	// local_external_work += deltamin_ptr->h[i]v[i];
	// }
	// double global_external_work = 0;
	// MPI_Allreduce(&local_external_work,&global_external_work,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
	// external_work += global_external_work;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> Real DomainLammps<Dim>::getExternalWork() {
	// return external_work;
	return 0;
	}
	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void DomainLammps<Dim>::setExternalWork(Real work) {
	// external_work = work;
	}
	/* -------------------------------------------------------------------------- */
	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 non-dimension");
	return false;
	}
	}

	/* -------------------------------------------------------------------------- */
	template <UInt Dim> void DomainLammps<Dim>::force_clear(UInt vflag) {
	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> void DomainLammps<Dim>::updatePAD() {

	// if (modify->n_pre_exchange)
	modify->pre_exchange();

	if (domain->triclinic)
	domain->x2lamda(atom->nlocal);

	domain->pbc();

	domain->reset_box();
	comm->setup();
	neighbor->setup_bins();
	// 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);

	// if (modify->n_pre_neighbor)
	modify->pre_neighbor();

	neighbor->build();

	// this->performStep2();

	this->update_pad = false;
	this->restart_start = false;
	}
	/* -------------------------------------------------------------------------- */
	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;
	}
	/* -------------------------------------------------------------------------- */
	/* 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);

	/* L2MKEYWORD GEOM_BY_TYPE
	Assign a geometry for the atom type provided
	*/

	// this->parseKeyword("GEOM_BY_TYPE", geom_by_type);

	/* LMKEYWORD MINCG_ETOL
	obsolete
	*/
	this->parseKeyword("MINCG_ETOL", mincg_etol, 1e-8);

	/* LMKEYWORD MINCG_FTOL
	osolete
	*/

	this->parseKeyword("MINCG_FTOL", mincg_ftol, 1e-8);

	/* LMKEYWORD MINCG_MAXEVAL
	obsolete
	*/

	this->parseKeyword("MINCG_MAXEVAL", mincg_maxeval, 1000000u);

	// // /* LMK EYWORD MINCG_DMAX
	// // obsolete
	// // */

	// this->parseKeyword("MINCG_DMAX",mincg_dmax,.1);
	}
	/* -------------------------------------------------------------------------- */

	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