domain_lammps.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Nov 15, 10:25

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 "lm_common.hh"
	#include "lammps_common.hh"
	#include "domain_lammps.hh"
	#include "reference_manager_lammps.hh"
	#include "import_lammps.hh"
	#include "filter_manager.hh"
	#include <iomanip>
	#include <iostream>
	#include <iterator>
	#include "trace_atom.hh"
	/* -------------------------------------------------------------------------- */
	//LAMMPS includes
	#include "integrate.h"
	#include "pair.h"
	#include "bond.h"
	#include "angle.h"
	#include "dihedral.h"
	#include "improper.h"
	#include "kspace.h"
	#include "change_box.h"
	#include "displace_box.h"
	#include "verlet.h"
	#include "variable.h"
	#include "fix_minimize.h"
	#include "min_cg.h"
	/* -------------------------------------------------------------------------- */
	#include "compute_lammps.hh"
	#include "geometry_manager.hh"
	#include "communicator.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) {
	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) {
	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(){
	Cube & c = GeometryManager::getManager().getGeometry(this->newGeomBox)->getBoundingBox();
	this->setDomainDimensions(c.getXmin().getValues(),c.getXmax().getValues());
	}

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



	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());

	std::map<UInt,GeomID>::iterator it = geom_by_type.begin();
	std::map<UInt,GeomID>::iterator end = geom_by_type.end();

	while (it != end){
	Geometry * g = GeometryManager::getManager().getObject(it->second);
	UInt itype = it->first;
	ImportLammpsInterface::setGeom(g,itype);
	++it;
	}

	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");
	}
	/* -------------------------------------------------------------------------- */
	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
	}
	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

	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(DomainID ID,CommGroup GID):
	DomainAtomic<ContainerLammps<Dim>,Dim>(ID,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 Real * xmin,const Real * 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) {
	timer->stamp();
	comm->forward_comm();
	timer->stamp(TIME_COMM);
	}
	else {
	//if starts from restart file, this part has to be fulfilled.
	this->restart_start=false;
	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);
	}
	if(this->when_change_box == current_time){
	//if one change box size, this part has to be fulfilled.
	this->restart_start=false;
	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);

	}
	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]));
	}
	if (&this->getGeomConstrained()){
	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;
	}

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

	/* 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 ocnditions.
	*/
	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);

	/* LMKEYWORD 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