Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/md/lammps/domain_lammps.cc
2afe98e0d677master

Last Change
Raw File

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 <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 <timer.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 = nullptr;
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <UInt Dim> void DomainLammps<Dim>::init() {
DomainAtomic<ContainerLammps<Dim>, Dim>::init();
if (!this->getCommGroup().amIinGroup())
return;
this->initModel();
this->update->whichflag = 1;
LAMMPS_NS::LAMMPS::init();
this->update->integrate->setup(0);
this->initDOFs();
this->ref_manager->setState(true);
this->atoms.getBoundingBox() = this->getDomainBoundingBox();
this->updateAcceleration();
}
/* -------------------------------------------------------------------------- */
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:
if (not is_mass_per_particle)
return ArrayView(&lammps_main_object->atom->mass[0],
lammps_main_object->atom->ntypes, 1);
return ArrayView(&lammps_main_object->atom->rmass[0],
lammps_main_object->atom->natoms, 1);
case _acceleration:
return ArrayView(this->acceleration_field.data(),
lammps_main_object->atom->nlocal, 3);
case _angular_velocity:
if (not lammps_main_object->atom->omega) {
LM_FATAL("created atoms do not have angle degrees of freedom");
}
return ArrayView(&lammps_main_object->atom->omega[0][0],
lammps_main_object->atom->nlocal, 3);
case _angular_acceleration:
return ArrayView(this->angular_acceleration_field.data(),
lammps_main_object->atom->nlocal, 3);
case _torque:
if (not lammps_main_object->atom->torque) {
LM_FATAL("created atoms do not have angle degrees of freedom");
}
return ArrayView(&lammps_main_object->atom->torque[0][0],
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
auto &vars = Parser::getAlgebraicVariables();
for (auto &&[varname, val] : vars) {
std::stringstream varval;
varval << std::setprecision(16) << std::scientific << val;
LAMMPS_NS::Variable *lmpvar = LAMMPS_NS::LAMMPS::input->variable;
char cvarname[256];
varname.copy(cvarname, 255);
cvarname[varname.length()] = 0;
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);
}
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 == real_unit_system)
unit_name = "real";
else if (code_unit_system == si_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");
this->acceleration_field.resize(atom->nlocal, 3);
this->acceleration_field.setZero();
if (std::string(this->atom->atom_style) != "atomic") {
this->is_mass_per_particle = true;
this->angular_acceleration_field.resize(atom->nlocal, 3);
this->angular_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(real_unit_system);
else if (strcmp(update->unit_style, "metal") == 0)
uc.setInUnits(metal_unit_system);
else if (strcmp(update->unit_style, "si") == 0)
uc.setInUnits(si_unit_system);
else
LM_FATAL("not known lammps unit system:"
<< " things need to be done to integrate");
uc.computeConversions();
if (not is_mass_per_particle) {
for (int i = 1; i <= atom->ntypes; i++)
atom->mass[i] *= uc.getConversion<Mass>();
} else {
for (int i = 0; i <= atom->natoms; i++)
atom->rmass[i] *= uc.getConversion<Mass>();
}
}
#ifdef LIBMULTISCALE_TRACE_ATOM
for (int i = 0; i < atom->nlocal; ++i) {
Vector<3> X = {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() const {
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() const {
Cube c;
if (!this->getCommGroup().amIinGroup())
return 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(false);
}
/* -------------------------------------------------------------------------- */
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(LAMMPS_NS::Timer::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(LAMMPS_NS::Timer::COMM);
if (modify->n_pre_neighbor)
modify->pre_neighbor();
neighbor->build(0);
timer->stamp(LAMMPS_NS::Timer::NEIGH);
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(LAMMPS_NS::Timer::COMM);
if (modify->n_pre_neighbor)
modify->pre_neighbor();
neighbor->build(false);
timer->stamp(LAMMPS_NS::Timer::NEIGH);
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 (lammps_main_object->atom->torque) {
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->changeRelease();
}
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(LAMMPS_NS::Timer::OUTPUT);
}
STOPTIMER("AtomStepVelocities");
VIEW_ATOM(RefLammps<Dim>);
this->changeRelease();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void DomainLammps<Dim>::updateAcceleration() {
if (!this->getCommGroup().amIinGroup())
return;
double *rmass = atom->rmass;
double *mass = atom->mass;
double *radius = atom->radius;
double inertia = 0.4;
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);
}
if (std::string(this->atom->atom_style) != "atomic") {
auto torque = this->getField(_torque);
for (UInt i = 0; i < UInt(lammps_main_object->atom->nlocal); ++i) {
Real mass_inverse = code_unit_system.ft_m2v / inertia /
(radius[i] * radius[i] * rmass[i]);
this->angular_acceleration_field(i, 0) = mass_inverse * torque(i, 0);
this->angular_acceleration_field(i, 1) = mass_inverse * torque(i, 1);
this->angular_acceleration_field(i, 2) = mass_inverse * torque(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);
timer->stamp(LAMMPS_NS::Timer::MODIFY);
}
if (this->force->pair && this->force->pair->compute_flag) {
this->force->pair->compute(eflag, vflag);
this->timer->stamp(LAMMPS_NS::Timer::PAIR);
}
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(LAMMPS_NS::Timer::BOND);
}
if (this->force->kspace) {
this->force->kspace->compute(eflag, vflag);
this->timer->stamp(LAMMPS_NS::Timer::KSPACE);
}
if (modify->n_pre_reverse) {
modify->pre_reverse(eflag, vflag);
timer->stamp(LAMMPS_NS::Timer::MODIFY);
}
if (this->force->newton) {
this->comm->reverse_comm();
this->timer->stamp(LAMMPS_NS::Timer::COMM);
}
if (this->modify->n_post_force)
this->modify->post_force(vflag);
STOPTIMER("AtomStepForces");
VIEW_ATOM(RefLammps<Dim>);
this->changeRelease();
}
/* -------------------------------------------------------------------------- */
/* LMDESC LAMMPS
This plugin make the interface with lammps
*/
/* LMHERITANCE domain_md */
/* LM EXAMPLE
Section MultiScale RealUnits\\
...\\
GEOMETRY myGeom CUBE BBOX -1 1 -1 1 -1 1\\
MODEL LAMMPS md\\
...\\
endSection\\ \\
Section LAMMPS:md RealUnits \\
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__
c4science · Help