LammpsInterface.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Jun 30, 10:33

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_; }

	}

LammpsInterface.cppNo OneTemporaryActions

File Metadata

LammpsInterface.cppView Options

Event Timeline

LammpsInterface.cpp
No OneTemporary
Actions

LammpsInterface.cpp
View Options