Page MenuHomec4science
Files F91413339

LammpsInterface.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Nov 10, 21:14

LammpsInterface.cpp
View Options

// Header file for this class
#include "LammpsInterface.h"
// LAMMPS includes
#include "lammps.h"
#include "lmptype.h"
#include "atom.h" // x, v, f
#include "domain.h" // for basing locations on regions
#include "region.h" // region bounding box and style
#include "force.h" // boltzman constant
#include "group.h" // atom masks
#include "memory.h" // grow atom information
#include "compute.h" // computes
#include "modify.h" //
#include "neighbor.h" // neighbors
#include "neigh_list.h" // neighbor list
#include "update.h" // timestepping information
#include "pair.h" // pair potentials
#include "lattice.h" // lattice parameters
#include "comm.h" //
// ATC includes
#include "ATC_Error.h"
#include "MatrixLibrary.h"
// Other include files
#include "mpi.h"
#include <cstring>
namespace ATC {
LammpsInterface * LammpsInterface::myInstance_ = NULL;
// -----------------------------------------------------------------
// instance()
// -----------------------------------------------------------------
LammpsInterface * LammpsInterface::instance()
{
if (myInstance_ == NULL) {
myInstance_ = new LammpsInterface();
}
return myInstance_;
}
// -----------------------------------------------------------------
// constructor
// -----------------------------------------------------------------
LammpsInterface::LammpsInterface()
: lammps_(NULL),
atomPE_(NULL),
commRank_(0)
{
}
// -----------------------------------------------------------------
// general interface methods
// -----------------------------------------------------------------
MPI_Comm LammpsInterface::world() { return lammps_->world; }
// -----------------------------------------------------------------
// atom interface methods
// -----------------------------------------------------------------
int LammpsInterface::nlocal() { return lammps_->atom->nlocal; }
int LammpsInterface::nghost() { return lammps_->atom->nghost; }
double LammpsInterface::natoms() { return lammps_->atom->natoms; }
int LammpsInterface::nmax() { return lammps_->atom->nmax; }
int LammpsInterface::ntypes() { return lammps_->atom->ntypes; }
double ** LammpsInterface::xatom() { return lammps_->atom->x; }
const double ** LammpsInterface::xatom() const { return (const double**)(lammps_->atom->x); }
double ** LammpsInterface::vatom() { return lammps_->atom->v; }
double ** LammpsInterface::fatom() { return lammps_->atom->f; }
int * LammpsInterface::atom_mask() { return lammps_->atom->mask; }
int * LammpsInterface::atom_type() { return lammps_->atom->type; }
int * LammpsInterface::atom_tag() { return lammps_->atom->tag; }
double * LammpsInterface::atom_mass() { return lammps_->atom->mass; }
double LammpsInterface::atom_mass(int iType) { return lammps_->atom->mass[iType]; }
double * LammpsInterface::atom_rmass() { return lammps_->atom->rmass; }
double * LammpsInterface::atom_charge() { return lammps_->atom->q; }
// -----------------------------------------------------------------
// domain interface methods
// -----------------------------------------------------------------
int LammpsInterface::dimension() { return lammps_->domain->dimension; }
int LammpsInterface::nregion() { return lammps_->domain->nregion; }
void LammpsInterface::get_box_bounds(double & boxxlo, double & boxxhi,
double & boxylo, double & boxyhi,
double & boxzlo, double &boxzhi)
{
if (lammps_->domain->triclinic == 0) {
boxxlo = lammps_->domain->boxlo[0];
boxxhi = lammps_->domain->boxhi[0];
boxylo = lammps_->domain->boxlo[1];
boxyhi = lammps_->domain->boxhi[1];
boxzlo = lammps_->domain->boxlo[2];
boxzhi = lammps_->domain->boxhi[2];
}
else {
boxxlo = lammps_->domain->boxlo_bound[0];
boxxhi = lammps_->domain->boxhi_bound[0];
boxylo = lammps_->domain->boxlo_bound[1];
boxyhi = lammps_->domain->boxhi_bound[1];
boxzlo = lammps_->domain->boxlo_bound[2];
boxzhi = lammps_->domain->boxhi_bound[2];
}
}
int LammpsInterface::xperiodic() { return lammps_->domain->xperiodic; }
int LammpsInterface::yperiodic() { return lammps_->domain->yperiodic; }
int LammpsInterface::zperiodic() { return lammps_->domain->zperiodic; }
int LammpsInterface::nperiodic()
{
int nprd = 0;
if ( lammps_->domain->xperiodic > 0 ) { nprd++ ; }
if ( lammps_->domain->yperiodic > 0 ) { nprd++ ; }
if ( lammps_->domain->zperiodic > 0 ) { nprd++ ; }
return nprd;
}
double LammpsInterface::domain_xprd() { return lammps_->domain->xprd; }
double LammpsInterface::domain_yprd() { return lammps_->domain->yprd; }
double LammpsInterface::domain_zprd() { return lammps_->domain->zprd; }
double LammpsInterface::domain_xy() { return lammps_->domain->xy; }
double LammpsInterface::domain_xz() { return lammps_->domain->xz; }
double LammpsInterface::domain_yz() { return lammps_->domain->yz; }
int LammpsInterface::domain_triclinic() { return lammps_->domain->triclinic; }
void LammpsInterface::get_box_periodicity(int & xperiodic,
int & yperiodic,
int & zperiodic)
{
xperiodic = lammps_->domain->xperiodic;
yperiodic = lammps_->domain->yperiodic;
zperiodic = lammps_->domain->zperiodic;
}
int LammpsInterface::get_region_id(const char * regionName) {
int nregion = this->nregion();
for (int iregion = 0; iregion < nregion; iregion++) {
if (strcmp(regionName, region_name(iregion)) == 0) {
return iregion;
}
}
throw ATC_Error(0,"Region has not been defined");
return -1;
}
// -----------------------------------------------------------------
// update interface methods
// -----------------------------------------------------------------
LammpsInterface::UnitsType LammpsInterface::units_style(void)
{
if (strcmp(lammps_->update->unit_style,"lj") == 0) return LJ;
else if (strcmp(lammps_->update->unit_style,"real") == 0) return REAL;
else if (strcmp(lammps_->update->unit_style,"metal") == 0) return METAL;
else return UNKNOWN;
}
// -----------------------------------------------------------------
// lattice interface methods
// -----------------------------------------------------------------
double LammpsInterface::xlattice() { return lammps_->domain->lattice->xlattice; }
double LammpsInterface::ylattice() { return lammps_->domain->lattice->ylattice; }
double LammpsInterface::zlattice() { return lammps_->domain->lattice->zlattice; }
LammpsInterface::LatticeType LammpsInterface::lattice_style()
{
if (lammps_->domain->lattice)
return (LammpsInterface::LatticeType)lammps_->domain->lattice->style;
else
throw ATC_Error(0,"Lattice has not been defined");
}
//* retuns the number of basis vectors
int LammpsInterface::get_n_basis()
{
return lammps_->domain->lattice->nbasis;
}
//* returns the basis vectors, transformed to the box coords
void LammpsInterface::get_basis(double **basis)
{
LAMMPS_NS::Lattice *lattice = lammps_->domain->lattice;
int i,j;
double origin[3] = {0.0, 0.0, 0.0};
lattice->lattice2box(origin[0], origin[1], origin[2]);
for (i=0; i<get_n_basis(); i++)
{
memcpy(basis[i],lattice->basis[i],3*sizeof(double));
lattice->lattice2box(basis[i][0], basis[i][1], basis[i][2]);
for (j=0; j<3; j++) basis[i][j] -= origin[j];
}
}
//* gets the unit cell vectors
void LammpsInterface::get_unit_cell(double *a1, double *a2, double *a3)
{
int i, j;
double *a[3] = {a1,a2,a3};
double origin[3] = {0.0,0.0,0.0};
LAMMPS_NS::Lattice *lattice = lammps_->domain->lattice;
// transform origin
lattice->lattice2box(origin[0], origin[1], origin[2]);
// copy reference lattice vectors
memcpy(a[0], lattice->a1, 3*sizeof(double));
memcpy(a[1], lattice->a2, 3*sizeof(double));
memcpy(a[2], lattice->a3, 3*sizeof(double));
for (i=0; i<3; i++)
{
lattice->lattice2box(a[i][0], a[i][1], a[i][2]);
for (j=0; j<3; j++) a[i][j] -= origin[j];
}
}
//* gets number of atoms in a unit cell
int LammpsInterface::num_atoms_per_cell(void)
{
int naCell = 0;
LatticeType type = lattice_style();
if (type == LammpsInterface::SC) naCell = 1;
else if (type == LammpsInterface::BCC) naCell = 2;
else if (type == LammpsInterface::FCC) naCell = 4;
else if (type == LammpsInterface::DIAMOND) naCell = 8;
else if (comm_rank()==0) {
//{throw ATC_Error(0,"lattice style not currently supported by ATC");}
cout << "ATC WARNING: Cannot get number of atoms per cell from lattice\n";
naCell = 1; //HACK to enable us to keep going since this is only used to compute volume per atom
// ATC modes with a user specified atomic volume or using only volumetric quantities are fine
}
return naCell;
}
//* gets tributary volume for an atom
double LammpsInterface::volume_per_atom(void)
{
double naCell = num_atoms_per_cell();
double volPerAtom =
xlattice() * ylattice() * zlattice() / naCell;
return volPerAtom;
}
//* gets lattice basis
void LammpsInterface::get_lattice(MATRIX &N, MATRIX &B)
{
int nbasis = get_n_basis();
double **basis = new double*[nbasis];
N.reset(3,3);
B.reset(3,nbasis);
for (int i=0; i<nbasis; i++) basis[i] = column(B,i).get_ptr();
get_basis(basis);
get_unit_cell(column(N,0).get_ptr(),
column(N,1).get_ptr(),
column(N,2).get_ptr());
delete [] basis;
}
// -----------------------------------------------------------------
// force interface methods
// -----------------------------------------------------------------
double LammpsInterface::boltz() { return lammps_->force->boltz; }
double LammpsInterface::mvv2e() { return lammps_->force->mvv2e; }
double LammpsInterface::ftm2v() { return lammps_->force->ftm2v; }
double LammpsInterface::nktv2p() { return lammps_->force->nktv2p; }
double LammpsInterface::qqr2e() { return lammps_->force->qqr2e; }
double LammpsInterface::qe2f() { return lammps_->force->qe2f; }
double LammpsInterface::dielectric() { return lammps_->force->dielectric; }
double LammpsInterface::qqrd2e() { return lammps_->force->qqrd2e; }
double LammpsInterface::pair_force(int i, int j, double rsq,
double & fmag_over_rmag)
{
int itype = (lammps_->atom->type)[i];
int jtype = (lammps_->atom->type)[j];
// return value is the energy
if (rsq < (lammps_->force->pair->cutsq)[itype][jtype]) {
return lammps_->force->pair->single(i,j,itype,jtype,
rsq,1.0,1.0,fmag_over_rmag);
}
return 0.0;
}
double LammpsInterface::pair_cutoff()
{
return lammps_->force->pair->cutforce;
}
int LammpsInterface::single_enable()
{
return lammps_->force->pair->single_enable;
}
//* Boltzmann's constant in M,L,T,t units
double LammpsInterface::kBoltzmann() {
return (lammps_->force->boltz)/(lammps_->force->mvv2e);
}
//* Dulong-Petit heat capacity
double LammpsInterface::heat_capacity() {
double rhoCp = dimension()*kBoltzmann()/volume_per_atom();
return rhoCp;
}
//* reference mass density
double LammpsInterface::mass_density()
{
const int ntypes = lammps_->atom->ntypes;
const int *mass_setflag = lammps_->atom->mass_setflag;
const int *type = lammps_->atom->type;
const double *mass = lammps_->atom->mass;
const double *rmass = lammps_->atom->rmass;
// NOTE currently assumes all atoms have same mass and volume
// in the future, mass and volume will be different but density should be
// an atom indepedent quantity
if (mass) {
if (type) return mass[type[0]]/volume_per_atom();
// no type array - just use first type that has a set mass
for (int i=1; i<=ntypes; i++) {
if (mass_setflag[i]) return mass[i]/volume_per_atom();
}
// NOTE: no masses specified in input file should we warn the user of this?
return 0.0;
}
// NOTE is this valid - lammps likes to not use 0 index
if (rmass) return rmass[0]/volume_per_atom();
return 0.0;
}
// -----------------------------------------------------------------
// group interface methods
// -----------------------------------------------------------------
int LammpsInterface::ngroup() { return lammps_->group->ngroup; }
int LammpsInterface::group_bit(int iGroup) { return lammps_->group->bitmask[iGroup]; }
int LammpsInterface::find_group(const char * c) { return lammps_->group->find(c); }
int LammpsInterface::group_inverse_mask(int iGroup)
{
return lammps_->group->inversemask[iGroup];
}
char * LammpsInterface::group_name(int iGroup)
{
return lammps_->group->names[iGroup];
}
void LammpsInterface::group_bounds(int iGroup, double * b)
{
lammps_->group->bounds(iGroup, b);
}
// -----------------------------------------------------------------
// memory interface methods
// -----------------------------------------------------------------
double * LammpsInterface::create_1d_double_array(int nlo, int nhi, const char *name) {
double *array;
return lammps_->memory->create1d_offset(array, nlo, nhi, name);
}
void LammpsInterface::destroy_1d_double_array(double * d, int i) {
lammps_->memory->destroy1d_offset(d, i);
}
double ** LammpsInterface::create_2d_double_array(int n1, int n2, const char *name) {
double **array;
return lammps_->memory->create(array, n1, n2, name);
}
void LammpsInterface::destroy_2d_double_array(double **d) {
lammps_->memory->destroy(d);
}
double **LammpsInterface::grow_2d_double_array(double **array,
int n1,
int n2,
const char *name)
{
return lammps_->memory->grow(array, n1, n2, name);
}
int ** LammpsInterface::create_2d_int_array(int n1, int n2, const char *name) {
int **array;
return lammps_->memory->create(array, n1, n2, name);
}
void LammpsInterface::destroy_2d_int_array(int **i) {
lammps_->memory->destroy(i);
}
int ** LammpsInterface::grow_2d_int_array(int **array, int n1, int n2, const char *name) {
return lammps_->memory->grow(array, n1, n2, name);
}
// -----------------------------------------------------------------
// update interface methods
// -----------------------------------------------------------------
double LammpsInterface::dt() { return lammps_->update->dt; }
int LammpsInterface::ntimestep() { return lammps_->update->ntimestep; }
int LammpsInterface::nsteps() { return lammps_->update->nsteps; }
// -----------------------------------------------------------------
// neighbor list interface methods
// -----------------------------------------------------------------
void LammpsInterface::init_list(int id, LAMMPS_NS::NeighList *ptr) { list_ = ptr; }
int LammpsInterface::neighbor_list_inum() { return list_->inum; }
int * LammpsInterface::neighbor_list_numneigh() { return list_->numneigh; }
int * LammpsInterface::neighbor_list_ilist() { return list_->ilist; }
int ** LammpsInterface::neighbor_list_firstneigh() { return list_->firstneigh; }
int LammpsInterface::neighbor_ago() { return lammps_->neighbor->ago; }
// -----------------------------------------------------------------
// region interface methods
// -----------------------------------------------------------------
char * LammpsInterface::region_name(int iRegion)
{
return lammps_->domain->regions[iRegion]->id;
}
char * LammpsInterface::region_style(int iRegion)
{
return lammps_->domain->regions[iRegion]->style;
}
double LammpsInterface::region_xlo(int iRegion)
{
return lammps_->domain->regions[iRegion]->extent_xlo;
}
double LammpsInterface::region_xhi(int iRegion)
{
return lammps_->domain->regions[iRegion]->extent_xhi;
}
double LammpsInterface::region_ylo(int iRegion)
{
return lammps_->domain->regions[iRegion]->extent_ylo;
}
double LammpsInterface::region_yhi(int iRegion)
{
return lammps_->domain->regions[iRegion]->extent_yhi;
}
double LammpsInterface::region_zlo(int iRegion)
{
return lammps_->domain->regions[iRegion]->extent_zlo;
}
double LammpsInterface::region_zhi(int iRegion)
{
return lammps_->domain->regions[iRegion]->extent_zhi;
}
double LammpsInterface::region_xscale(int iRegion)
{
return lammps_->domain->regions[iRegion]->xscale;
}
double LammpsInterface::region_yscale(int iRegion)
{
return lammps_->domain->regions[iRegion]->yscale;
}
double LammpsInterface::region_zscale(int iRegion)
{
return lammps_->domain->regions[iRegion]->zscale;
}
int LammpsInterface::region_match(int iRegion, double x, double y, double z) {
return lammps_->domain->regions[iRegion]->match(x,y,z);
}
// -----------------------------------------------------------------
// compute methods
// -----------------------------------------------------------------
int LammpsInterface::find_compute(const char* tag)
{
// a clunky way to safely get rid of the const
int n = strlen(tag) + 1;
char* tag_copy = new char[n];
strcpy(tag_copy,tag);
int icompute = lammps_->modify->find_compute(tag_copy);
if (icompute < 0) {
string msg("Could not find compute ");
msg += tag;
throw ATC_Error(0,msg);
}
return icompute;
}
LAMMPS_NS::Compute* LammpsInterface::get_compute(const char* tag)
{
int icompute = find_compute(tag);
LAMMPS_NS::Compute* cmpt = lammps_->modify->compute[icompute];
return cmpt;
}
double** LammpsInterface::compute_vector_data(const char* tag)
{
LAMMPS_NS::Compute* cmpt = get_compute(tag);
if (!(cmpt->invoked_flag & INVOKED_PERATOM)) {
cmpt->compute_peratom();
cmpt->invoked_flag |= INVOKED_PERATOM;
}
return cmpt->array_atom;
}
double* LammpsInterface::compute_scalar_data(const char* tag)
{
LAMMPS_NS::Compute* cmpt = get_compute(tag);
if (!(cmpt->invoked_flag & INVOKED_PERATOM)) {
cmpt->compute_peratom();
cmpt->invoked_flag |= INVOKED_PERATOM;
}
return cmpt->vector_atom;
}
int LammpsInterface::compute_ncols(const char* tag)
{
int icompute = find_compute(tag);
int ncols = lammps_->modify->compute[icompute]->size_peratom_cols;
if (ncols == 0) ncols = 1; // oddity of lammps, used as flag for scalar
return ncols;
}
// -----------------------------------------------------------------
// compute pe/atom interface methods
// -----------------------------------------------------------------
int LammpsInterface::atomPE_create(void)
{
//char * list[3] = {"atcPE","all","pe/atom"};
char * list[4] = {"atcPE","all","pe/atom","pair"};
int icompute = lammps_->modify->find_compute(list[0]);
if (icompute < 0) {
lammps_->modify->add_compute(3,list);
icompute = lammps_->modify->find_compute(list[0]);
}
if (! atomPE_ ) {
atomPE_ = lammps_->modify->compute[icompute];
}
return icompute;
}
void LammpsInterface::atomPE_init(void)
{
if (atomPE_) {
atomPE_->init();
}
else {
throw ATC_Error(0,"no atomPE compute");
}
}
void LammpsInterface::atomPE_addstep(LAMMPS_NS::bigint step)
{
atomPE_->addstep(step);
}
double * LammpsInterface::atomPE_compute(void)
{
if (atomPE_) {
atomPE_->compute_peratom();
return atomPE_->vector_atom;
}
else {
return NULL;
}
}
/* ---------------------------------------------------------------------- */
void LammpsInterface::unwrap_coordinates(int iatom, double* xatom)
{
double **x = lammps_->atom->x;
int *mask = lammps_->atom->mask;
int *image = lammps_->atom->image;
int nlocal = lammps_->atom->nlocal;
double *h = lammps_->domain->h;
double xprd = lammps_->domain->xprd;
double yprd = lammps_->domain->yprd;
double zprd = lammps_->domain->zprd;
int xbox,ybox,zbox;
// for triclinic, need to unwrap current atom coord via h matrix
// NOTE: Using current box dimensions.
if (lammps_->domain->triclinic == 0) {
xbox = (image[iatom] & 1023) - 512;
ybox = (image[iatom] >> 10 & 1023) - 512;
zbox = (image[iatom] >> 20) - 512;
xatom[0] = x[iatom][0] + xbox*xprd;
xatom[1] = x[iatom][1] + ybox*yprd;
xatom[2] = x[iatom][2] + zbox*zprd;
} else {
xbox = (image[iatom] & 1023) - 512;
ybox = (image[iatom] >> 10 & 1023) - 512;
zbox = (image[iatom] >> 20) - 512;
xatom[0] = x[iatom][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
xatom[1] = x[iatom][1] + h[1]*ybox + h[3]*zbox;
xatom[2] = x[iatom][2] + h[2]*zbox;
}
}
// -----------------------------------------------------------------
// other methods
// -----------------------------------------------------------------
/** Return lammps pointer -- only as a last resort! */
LAMMPS_NS::LAMMPS * LammpsInterface::get_lammps_ptr() { return lammps_; }
}

Event Timeline

c4science · Help