Diffusion LibMultiScale (features/new_parralelism) Browse

new_parralelism

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__

libmultiscale/src/md/lammps/domain_lammps.ccb05d5f2aba78features/new_parralelism

domain_lammps.ccView Options

libmultiscale/src/md/lammps/domain_lammps.cc
b05d5f2aba78features/new_parralelism

domain_lammps.cc
View Options