ATC_Method.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Jul 5, 13:03

ATC_Method.cpp
View Options

	// ATC headers
	#include "ATC_Method.h"
	#include "LammpsInterface.h"
	#include "FE_Engine.h"
	#include "Array.h"
	#include "Array2D.h"
	#include "ATC_Error.h"
	#include "Function.h"
	#include "PrescribedDataManager.h"
	#include "TimeIntegrator.h"
	#include "PhysicsModel.h"
	#include "PerAtomQuantityLibrary.h"
	#include "TransferLibrary.h"
	#include "KernelFunction.h"
	#include "Utility.h"
	#include "FieldManager.h"

	#include <fstream>
	#include <sstream>
	#include <iostream>

	using ATC_Utility::sgn;
	using ATC_Utility::to_string;
	using ATC_Utility::is_dbl;

	using std::stringstream;
	using std::ifstream;
	using std::ofstream;
	using std::string;
	using std::map;
	using std::set;
	using std::vector;
	using std::pair;

	namespace ATC {

	ATC_Method::ATC_Method(string groupName, double ** & perAtomArray, LAMMPS_NS::Fix * thisFix) :
	nodalAtomicVolume_(NULL),
	needReset_(true),
	lammpsInterface_(LammpsInterface::instance()),
	interscaleManager_(this),
	timeFilterManager_(this),
	integrateInternalAtoms_(false),
	atomTimeIntegrator_(NULL),
	ghostManager_(this),
	feEngine_(NULL),
	initialized_(false),
	meshDataInitialized_(false),
	localStep_(0),
	sizeComm_(8), // 3 positions + 1 material id * 2 for output
	atomCoarseGrainingPositions_(NULL),
	atomGhostCoarseGrainingPositions_(NULL),
	atomProcGhostCoarseGrainingPositions_(NULL),
	atomReferencePositions_(NULL),
	nNodes_(0),
	nsd_(lammpsInterface_->dimension()),
	xref_(NULL),
	readXref_(false),
	needXrefProcessorGhosts_(false),
	trackDisplacement_(false),
	needsAtomToElementMap_(true),
	atomElement_(NULL),
	atomGhostElement_(NULL),
	internalElementSet_(""),
	atomMasses_(NULL),
	atomPositions_(NULL),
	atomVelocities_(NULL),
	atomForces_(NULL),
	parallelConsistency_(true),
	outputNow_(false),
	outputTime_(true),
	outputFrequency_(0),
	sampleFrequency_(0),
	sampleCounter_(0),
	peScale_(1./(lammpsInterface_->mvv2e())),
	keScale_(1.),
	scalarFlag_(0),
	vectorFlag_(0),
	sizeVector_(0),
	scalarVectorFreq_(0),
	sizePerAtomCols_(4),
	perAtomOutput_(NULL),
	perAtomArray_(perAtomArray),
	extScalar_(0),
	extVector_(0),
	extList_(NULL),
	thermoEnergyFlag_(0),
	atomVolume_(NULL),
	atomicWeightsWriteFlag_(false),
	atomicWeightsWriteFrequency_(0),
	atomWeightType_(LATTICE),
	domainDecomposition_(REPLICATED_MEMORY),
	groupbit_(0),
	groupbitGhost_(0),
	needProcGhost_(false),
	groupTag_(groupName),
	nLocal_(0),
	nLocalTotal_(0),
	nLocalGhost_(0),
	atomToElementMapType_(LAGRANGIAN),
	atomToElementMapFrequency_(0),
	regionID_(-1),
	mdMassNormalization_(false),
	kernelBased_(false),
	kernelOnTheFly_(false),
	kernelFunction_(NULL),
	bondOnTheFly_(false),
	accumulant_(NULL),
	accumulantMol_(NULL),
	accumulantMolGrad_(NULL),
	accumulantWeights_(NULL),
	accumulantInverseVolumes_(&invNodeVolumes_),
	accumulantBandwidth_(0),
	useRestart_(false),
	hasRefPE_(false),
	setRefPE_(false),
	setRefPEvalue_(false),
	refPEvalue_(0.),
	readRefPE_(false),
	nodalRefPotentialEnergy_(NULL),
	simTime_(0.0),
	stepCounter_(0)
	{
	lammpsInterface_->print_msg_once("version "+version());
	lammpsInterface_->set_fix_pointer(thisFix);
	interscaleManager_.set_lammps_data_prefix();
	grow_arrays(lammpsInterface_->nmax());
	feEngine_ = new FE_Engine(lammpsInterface_->world());


	lammpsInterface_->create_compute_pe_peratom();
	}

	ATC_Method::~ATC_Method()
	{
	lammpsInterface_->destroy_2d_double_array(xref_);
	lammpsInterface_->destroy_2d_double_array(perAtomOutput_);
	if (atomTimeIntegrator_) delete atomTimeIntegrator_;
	if (feEngine_) delete feEngine_;
	}

	//--------------------------------------------------
	// pack_fields
	// bundle all allocated field matrices into a list
	// for output needs
	//--------------------------------------------------
	void ATC_Method::pack_fields(RESTART_LIST & data)
	{
	map<FieldName,int>::const_iterator field;
	for (field = fieldSizes_.begin(); field!=fieldSizes_.end(); field++) {
	FieldName thisField = field->first;
	string fieldName = field_to_string(thisField);
	string matrixName;
	// copy all fields from ATC_Method.h
	matrixName = "fields_" + fieldName;
	data[matrixName] = & fields_[thisField].set_quantity();
	matrixName = "dot_fields_" + fieldName;
	data[matrixName] = & dot_fields_[thisField].set_quantity();
	matrixName = "ddot_fields_" + fieldName;
	data[matrixName] = & ddot_fields_[thisField].set_quantity();
	matrixName = "dddot_fields_" + fieldName;
	data[matrixName] = & dddot_fields_[thisField].set_quantity();
	matrixName = "NodalAtomicFieldsRoc_" + fieldName;
	data[matrixName] = & nodalAtomicFieldsRoc_[thisField].set_quantity();
	}
	}

	//--------------------------------------------------
	// write_restart_file
	// bundle matrices that need to be saved and call
	// fe_engine to write the file
	//--------------------------------------------------
	void ATC_Method::write_restart_data(string fileName, RESTART_LIST & data)
	{
	pack_fields(data);
	feEngine_->write_restart_file(fileName,&data);
	}

	//--------------------------------------------------
	// read_restart_file
	// bundle matrices that need to be saved and call
	// fe_engine to write the file
	//--------------------------------------------------
	void ATC_Method::read_restart_data(string fileName, RESTART_LIST & data)
	{
	pack_fields(data);
	feEngine_->read_restart_file(fileName,&data);
	}

	//--------------------------------------------------
	// Interactions with LAMMPS fix commands
	// parse input command and pass on to finite element engine
	// or physics specific transfers if necessary
	// revert to physics-specific transfer if no command matches input
	// first keyword is unique to particular class
	// base class keyword matching must apply to ALL physics
	// order: derived, base, owned objects
	//--------------------------------------------------
	bool ATC_Method::modify(int narg, char **arg)
	{
	int argIdx=0;
	bool match = false;

	// gateways to other modules e.g. extrinsic, control, mesh
	// pass off to fe engine
	if (strcmp(arg[argIdx],"mesh")==0) {
	match = feEngine_->modify(narg, arg);
	if (feEngine_->has_mesh() && !meshDataInitialized_)
	this->initialize_mesh_data();
	}
	// pass off to time filter
	else if (strcmp(arg[argIdx],"filter")==0) {
	argIdx++;
	match = timeFilterManager_.modify(narg-argIdx,&arg[argIdx]);

	// consistentInitialization_ = false;
	}
	// pass off to kernel function manager
	else if (strcmp(arg[argIdx],"kernel")==0) {
	argIdx++;

	if (kernelFunction_) {
	//delete kernelFunction_;
	//resetKernelFunction_ = true;
	}
	kernelFunction_ = KernelFunctionMgr::instance()->function(&arg[argIdx],narg-argIdx);
	if (kernelFunction_) match = true;
	else ATC_Error("no matching kernel found");
	feEngine_->set_kernel(kernelFunction_);

	accumulantMol_=&kernelAccumulantMol_; // KKM add
	accumulantMolGrad_=&kernelAccumulantMolGrad_; // KKM add
	}
	// pass off to ghost manager
	else if (strcmp(arg[argIdx],"boundary_dynamics")==0) {
	argIdx++;
	match = ghostManager_.modify(narg-argIdx,&arg[argIdx]);
	}

	// parsing handled here
	else {
	if (strcmp(arg[argIdx],"parallel_consistency")==0) {
	argIdx++;
	//if (!kernelFunction_) { throw ATC_Error("on_the_fly requires a kernel function"); }
	if (strcmp(arg[argIdx],"off")==0) parallelConsistency_ = false;
	else parallelConsistency_ = true;
	match = true;
	}
	/*! \page man_hardy_on_the_fly fix_modify AtC on_the_fly
	\section syntax
	fix_modify AtC on_the_fly <bond \| kernel> <optional on \| off> \n - bond \| kernel (keyword) = specifies on-the-fly calculation of bond or
	kernel matrix elements \n
	- on \| off (keyword) = activate or discontinue on-the-fly mode \n
	\section examples
	<TT> fix_modify AtC on_the_fly bond on </TT> \n <TT> fix_modify AtC on_the_fly kernel </TT> \n
	<TT> fix_modify AtC on_the_fly kernel off </TT> \n
	\section description
	Overrides normal mode of pre-calculating and storing bond pair-to-node a
	nd
	kernel atom-to-node matrices. If activated, will calculate elements of t
	hese
	matrices during repeated calls of field computations (i.e. "on-the-fly") and not store them for
	future use. \n on flag is optional - if omitted, on_the_fly will be activated for the s
	pecified
	matrix. Can be deactivated using off flag. \n
	\section restrictions
	Must be used with the hardy/field type of AtC fix
	( see \ref man_fix_atc )
	\section related
	\section default
	By default, on-the-fly calculation is not active (i.e. off). However, code does a memory allocation check to determine if it can store all needed bond and kernel matrix ele ments. If this allocation fails, on-the-fly is activated. \n
	*/

	else if (strcmp(arg[argIdx],"on_the_fly")==0) {
	argIdx++;
	//if (!kernelFunction_) { throw ATC_Error("on_the_fly requires a kernel function"); }
	if (strcmp(arg[argIdx],"bond")==0) {
	argIdx++;
	bondOnTheFly_ = true;
	if (narg > argIdx && strcmp(arg[argIdx],"off")==0) bondOnTheFly_ = false;
	}
	else if (strcmp(arg[argIdx],"kernel")==0) {
	argIdx++;
	kernelOnTheFly_ = true;
	if (narg > argIdx && strcmp(arg[argIdx],"off")==0) kernelOnTheFly_ = false;
	}
	else { throw ATC_Error("unsupported on_the_fly type"); }
	match = true;
	}

	/*! \page man_output fix_modify AtC output
	\section syntax
	fix_modify AtC output <filename_prefix> <frequency>
	[text \| full_text \| binary \| vector_components \| tensor_components ]
	fix_modify AtC output index [step \| time ]
	- filename_prefix (string) = prefix for data files
	- frequency (integer) = frequency of output in time-steps
	- options (keyword/s): \n
	text = creates text output of index, step and nodal variable values for unique nodes \n
	full_text = creates text output index, nodal id, step, nodal coordinates and nodal variable values for unique and image nodes \n
	binary = creates binary Ensight output \n
	vector_components = outputs vectors as scalar components \n
	tensor_components = outputs tensor as scalar components
	(use this for Paraview)\n

	\section examples
	<TT> fix_modify AtC output heatFE 100 </TT> \n
	<TT> fix_modify AtC output hardyFE 1 text tensor_components </TT> \n
	<TT> fix_modify AtC output hardyFE 10 text binary tensor_components </TT> \n
	<TT> fix_modify AtC output index step </TT> \n
	\section description
	Creates text and/or binary (Ensight, "gold" format) output of nodal/mesh data
	which is transfer/physics specific. Output indexed by step or time is possible.
	\section restrictions
	\section related
	see \ref man_fix_atc
	\section default
	no default format
	output indexed by time
	*/
	else if (strcmp(arg[argIdx],"output")==0) {
	argIdx++;
	/*! \page man_output_nodeset fix_modify AtC output nodeset
	\section syntax
	fix_modify AtC output nodeset <nodeset_name> <operation>
	- nodeset_name (string) = name of nodeset to be operated on
	- operation (keyword/s): \n
	sum = creates nodal sum over nodes in specified nodeset \n
	\section examples
	<TT> fix_modify AtC output nodeset nset1 sum </TT> \n
	\section description
	Performs operation over the nodes belonging to specified nodeset
	and outputs resulting variable values to GLOBALS file.
	\section restrictions
	\section related
	see \ref man_fix_atc
	\section default
	none
	*/
	if (strcmp(arg[argIdx],"nodeset")==0) {
	argIdx++;
	string nset = arg[argIdx++];
	if (strcmp(arg[argIdx],"sum")==0) {
	argIdx++;
	string field = arg[argIdx];
	pair < string, FieldName > id(nset,string_to_field(field));
	nsetData_[id] = NODESET_SUM;
	match = true;
	}
	else if (strcmp(arg[argIdx],"average")==0) {
	argIdx++;
	string field = arg[argIdx];
	pair < string, FieldName > id(nset,string_to_field(field));
	nsetData_[id] = NODESET_AVERAGE;
	match = true;
	}
	}

	/*! \page man_boundary_integral fix_modify AtC output boundary_integral
	\section syntax
	fix_modify AtC output boundary_integral [field] faceset [name]
	- field (string) : name of hardy field
	- name (string) : name of faceset
	\section examples
	<TT> fix_modify AtC output boundary_integral stress faceset loop1 </TT> \n
	\section description
	Calculates a surface integral of the given field dotted with the
	outward normal of the faces and puts output in the "GLOBALS" file
	\section restrictions
	Must be used with the hardy/field type of AtC fix
	( see \ref man_fix_atc )
	\section related
	\section default
	none
	*/

	/*! \page man_contour_integral fix_modify AtC output contour_integral
	\section syntax
	fix_modify AtC output contour_integral [field] faceset [name] <axis [x \| y \| z
	]>
	- field (string) : name of hardy field
	- name (string) : name of faceset
	- axis (string) : x or y or z
	\section examples
	<TT> fix_modify AtC output contour_integral stress faceset loop1 </TT> \n
	\section description
	Calculates a surface integral of the given field dotted with the
	outward normal of the faces and puts output in the "GLOBALS" file
	\section restrictions
	Must be used with the hardy/field type of AtC fix
	( see \ref man_fix_atc )
	\section related
	\section default
	none
	*/

	else if ( (strcmp(arg[argIdx],"boundary_integral")==0)
	\|\| (strcmp(arg[argIdx],"contour_integral")==0) ) {
	FacesetIntegralType type = BOUNDARY_INTEGRAL;
	if (strcmp(arg[argIdx],"contour_integral")==0)
	type = CONTOUR_INTEGRAL;
	argIdx++;
	string field(arg[argIdx++]);
	if(strcmp(arg[argIdx],"faceset")==0) {
	argIdx++;
	string name(arg[argIdx++]);
	pair <string,string> pair_name(name,field);
	fsetData_[pair_name] = type;
	match = true;
	}
	} // end "boundary_integral" \|\| "contour_integral"

	/*! \page man_output_elementset fix_modify AtC output elementset
	\section syntax
	fix_modify AtC output volume_integral <eset_name> <field>
	- set_name (string) = name of elementset to be integrated over
	- fieldname (string) = name of field to integrate
	csum = creates nodal sum over nodes in specified nodeset \n
	\section examples
	<TT> fix_modify AtC output eset1 mass_density </TT> \n
	\section description
	Performs volume integration of specified field over elementset
	and outputs resulting variable values to GLOBALS file.
	\section restrictions
	\section related
	see \ref man_fix_atc
	\section default
	none
	*/

	else if ( (strcmp(arg[argIdx],"volume_integral")==0) ) {
	argIdx++;
	string name(arg[argIdx++]);
	string field(arg[argIdx++]);
	pair <string,FieldName> pair_name(name,string_to_field(field));
	if (++argIdx < narg) { // keyword average
	esetData_[pair_name] = ELEMENTSET_AVERAGE;
	}
	else {
	esetData_[pair_name] = ELEMENTSET_TOTAL;
	}
	match = true;
	}

	else if (strcmp(arg[argIdx],"now")==0) {
	argIdx++;
	double dt = 1.0;
	if (argIdx < narg) {
	dt = atof(arg[argIdx++]);
	}
	update_time(dt);
	update_step();
	outputNow_ = true;
	this->output();
	outputNow_ = false;
	match = true;
	}
	else
	if (strcmp(arg[argIdx],"index")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"step")==0) { outputTime_ = false; }
	else { outputTime_ = true; }
	match = true;
	}
	else {
	outputPrefix_ = arg[argIdx++];
	outputFrequency_ = atoi(arg[argIdx++]);
	bool ensight_output = false, full_text_output = false;
	bool text_output = false, vect_comp = false, tensor_comp = false;
	int rank = lammpsInterface_->comm_rank();
	for (int i = argIdx; i<narg; ++i) {
	if (strcmp(arg[i],"full_text")==0) full_text_output = true;
	else if (strcmp(arg[i],"text")==0) text_output = true;
	else if (strcmp(arg[i],"binary")==0) ensight_output = true;
	else if (strcmp(arg[i],"vector_components")==0) vect_comp = true;
	else if (strcmp(arg[i],"tensor_components")==0) tensor_comp = true;
	else { throw ATC_Error(" output: unknown keyword "); }
	}
	if (outputFrequency_>0) {
	set<int> otypes;
	if (full_text_output \|\| text_output) {
	lammpsInterface_->print_msg_once("Warning : text output can create _LARGE_ files");
	}
	if (full_text_output) otypes.insert(FULL_GNUPLOT);
	if (text_output) otypes.insert(GNUPLOT);
	if (ensight_output) otypes.insert(ENSIGHT);
	if (ntracked() > 0) {
	string fstem = field_to_string(SPECIES_CONCENTRATION);
	string istem = field_to_intrinsic_name(SPECIES_CONCENTRATION);
	vector<string> tnames = tracked_names();
	vector<string> fnames;
	vector<string> inames;
	for (unsigned int i = 0; i < tnames.size(); i++) {
	fnames.push_back(fstem+tnames[i]);
	inames.push_back(istem+tnames[i]);
	}
	feEngine_->add_field_names(fstem,fnames);
	feEngine_->add_field_names(istem,inames);
	}
	feEngine_->initialize_output(rank,outputPrefix_,otypes);
	if (vect_comp)
	feEngine_->output_manager()
	->set_option(OUTPUT_VECTOR_COMPONENTS,true);
	if (tensor_comp)
	feEngine_->output_manager()
	->set_option(OUTPUT_TENSOR_COMPONENTS,true);
	}
	match = true;
	}
	}
	else if (strcmp(arg[argIdx],"write")==0) {
	argIdx++;
	FieldName thisField;
	int thisIndex;
	parse_field(arg,argIdx,thisField,thisIndex);
	string nsetName(arg[argIdx++]);
	string filename(arg[argIdx++]);
	stringstream f;
	set<int> nodeSet = (feEngine_->fe_mesh())->nodeset(nsetName);
	set<int>::const_iterator iset;
	const DENS_MAT & field =(fields_.find(thisField)->second).quantity();
	for (iset = nodeSet.begin(); iset != nodeSet.end(); iset++) {
	int inode = *iset;
	double v = field(inode,thisIndex);
	f << inode << " " << std::setprecision(17) << v << "\n";
	}
	LammpsInterface::instance()->write_file(filename,f.str());
	match = true;
	}
	// add a species for tracking
	/*! \page man_add_species fix_modify AtC add_species
	\section syntax
	fix_modify_AtC add_species <TAG> <group\|type> <ID> \n
	- <TAG> = tag for tracking a species \n
	- group\|type = LAMMPS defined group or type of atoms \n
	- <ID> = name of group or type number \n
	\section examples
	<TT> fix_modify AtC add_species gold type 1 </TT> \n
	<TT> group GOLDGROUP type 1 </TT> \n
	<TT> fix_modify AtC add_species gold group GOLDGROUP </TT>
	\section description
	Associates a tag with all atoms of a specified type or within a specified group. \n
	\section restrictions
	\section related
	\section default
	No defaults for this command.
	*/
	else if (strcmp(arg[argIdx],"add_species")==0) {
	argIdx++;
	string speciesTag = arg[argIdx];
	string tag = arg[argIdx];
	argIdx++;
	if (strcmp(arg[argIdx],"group")==0) {
	if (narg-argIdx == 2) {
	string name = arg[++argIdx];
	int id = lammpsInterface_->group_bit(name);
	groupList_.push_back(id);
	groupNames_.push_back(tag);
	}
	else {
	while (++argIdx < narg) {
	string name = arg[argIdx];
	int id = lammpsInterface_->group_bit(name);
	string tag = speciesTag+"-"+name;
	groupList_.push_back(id);
	groupNames_.push_back(tag);
	}
	}
	}
	else if (strcmp(arg[argIdx],"type")==0) {
	if (narg-argIdx == 2) {
	int id = atoi(arg[++argIdx]);
	typeList_.push_back(id);
	typeNames_.push_back(tag);
	}
	else {
	while (++argIdx < narg) {
	int id = atoi(arg[argIdx]);
	string tag = speciesTag+"_"+to_string(id);
	typeList_.push_back(id);
	typeNames_.push_back(tag);
	}
	}
	}
	else {
	throw ATC_Error("ATC_Method: add_species only handles groups or types"); }
	match = true;
	}

	// remove species from tracking

	/*! \page man_remove_species fix_modify AtC remove_species
	\section syntax
	fix_modify_AtC delete_species <TAG> \n

	- <TAG> = tag for tracking a species \n
	\section examples
	<TT> fix_modify AtC remove_species gold </TT> \n
	\section description
	Removes tag designated for tracking a specified species. \n
	\section restrictions
	\section related
	\section default
	No defaults for this command.
	*/
	else if (strcmp(arg[argIdx],"delete_species")==0) {
	argIdx++;
	string tag = arg[argIdx++];
	if (strcmp(arg[argIdx],"group")==0) {
	for (unsigned int j = 0; j < groupList_.size(); j++) {
	if (tag == groupNames_[j]) {
	groupList_.erase(groupList_.begin()+j);
	groupNames_.erase(groupNames_.begin()+j);
	break;
	}
	}
	}
	else if (strcmp(arg[argIdx],"type")==0) {
	for (unsigned int j = 0; j < typeList_.size(); j++) {
	if (tag == typeNames_[j]) {
	typeList_.erase(typeList_.begin()+j);
	typeNames_.erase(typeNames_.begin()+j);
	break;
	}
	}
	}
	else {
	throw ATC_Error("ATC_Method: delete_species only handles groups or types"); }
	match = true;

	}

	// add a molecule for tracking
	/*! \page man_add_molecule fix_modify AtC add_molecule
	\section syntax
	fix_modify_AtC add_molecule <small\|large> <TAG> <GROUP_NAME> \n

	- small\|large = can be small if molecule size < cutoff radius, must be large otherwise \n
	- <TAG> = tag for tracking a species \n
	- <GROUP_NAME> = name of group that tracking will be applied to \n
	\section examples
	<TT> group WATERGROUP type 1 2 </TT> \n
	<TT> fix_modify AtC add_molecule small water WATERGROUP </TT> \n
	\section description
	Associates a tag with all molecules corresponding to a specified group. \n
	\section restrictions
	\section related
	\section default
	No defaults for this command.
	*/
	else if (strcmp(arg[argIdx],"add_molecule")==0) {
	argIdx++;
	MolSize size;
	if (strcmp(arg[argIdx],"small")==0) {
	size = MOL_SMALL;
	//needXrefProcessorGhosts_ = true;
	needProcGhost_ = true;
	}
	else
	throw ATC_Error("ATC_CouplingMass: Bad molecule size in add_molecule");
	argIdx++;
	string moleculeTag = arg[argIdx];

	argIdx++;
	int groupBit = lammpsInterface_->group_bit(arg[argIdx]);
	moleculeIds_[moleculeTag] = pair<MolSize,int>(size,groupBit);
	match = true;
	}

	// remove molecule from tracking
	/*! \page man_remove_molecule fix_modify AtC remove_molecule
	\section syntax
	fix_modify_AtC remove_molecule <TAG> \n

	- <TAG> = tag for tracking a molecule type \n
	\section examples
	<TT> fix_modify AtC remove_molecule water </TT> \n
	\section description
	Removes tag designated for tracking a specified set of molecules. \n
	\section restrictions
	\section related
	\section default
	No defaults for this command.
	*/
	else if (strcmp(arg[argIdx],"remove_molecule")==0) {
	argIdx++;
	string moleculeTag = arg[argIdx];
	moleculeIds_.erase(moleculeTag);

	taggedDensMan_.erase(moleculeTag);
	}

	/*! \page man_boundary fix_modify AtC boundary
	\section syntax
	fix_modify AtC boundary type <atom-type-id>
	- <atom-type-id> = type id for atoms that represent a ficticious
	boundary internal to the FE mesh
	\section examples
	<TT> fix_modify AtC boundary type ghost_atoms </TT>
	\section description
	Command to define the atoms that represent the ficticious
	boundary internal to the FE mesh. For fully overlapped MD/FE
	domains with periodic boundary conditions no boundary atoms should
	be defined.
	\section restrictions
	\section default
	none
	*/
	else if (strcmp(arg[argIdx],"boundary")==0) {
	argIdx++;
	groupTagGhost_ = arg[argIdx++];
	match = true;
	}

	/*! \page man_internal_atom_integrate fix_modify AtC internal_atom_integrate
	\section syntax
	fix_modify AtC internal_atom_integrate <on \| off>
	<TT> fix_modify AtC internal_atom_integrate on </TT>
	\section description
	Has AtC perform time integration for the atoms in the group on which it operates. This does not include boundary atoms.
	\section restrictions
	AtC must be created before any fixes doing time integration. It must be on for coupling methods which impose constraints on velocities during the first verlet step, e.g. control momentum glc_velocity.
	\section default
	on for coupling methods, off for post-processors
	off
	*/
	else if (strcmp(arg[argIdx],"internal_atom_integrate")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"off")==0) {
	integrateInternalAtoms_ = false;
	match = true;
	}
	else {
	integrateInternalAtoms_ = true;
	match = true;
	}
	}

	/*! \page man_internal_element_set fix_modify AtC internal_element_set
	\section syntax
	fix_modify AtC internal_element_set <element-set-name>
	- <element-set-name> = name of element set defining internal region, or off
	\section examples
	<TT> fix_modify AtC internal_element_set myElementSet </TT>
	<TT> fix_modify AtC internal_element_set off </TT>
	\section description
	Enables AtC to base the region for internal atoms to be an element set.
	If no ghost atoms are used, all the AtC atoms must be constrained to remain
	in this element set by the user, e.g., with walls. If boundary atoms are
	used in conjunction with Eulerian atom maps
	AtC will partition all atoms of a boundary or internal type to be of type internal
	if they are in the internal region or to be of type boundary otherwise.
	\section restrictions
	If boundary atoms are used in conjunction with Eulerian atom maps, the Eulerian
	reset frequency must be an integer multiple of the Lammps reneighbor frequency
	\section related
	see \ref atom_element_map_type and \ref boundary
	\section default
	off
	*/
	else if (strcmp(arg[argIdx],"internal_element_set")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"off")==0) {
	internalElementSet_ = "";
	match = true;
	}
	else {
	internalElementSet_ = string(arg[argIdx]);
	const set<int> & elementSet((feEngine_->fe_mesh())->elementset(internalElementSet_)); // check it exists and is not trivial
	if (elementSet.size()==0) throw ATC_Error("internal_element_set - element set " + internalElementSet_ + " has no elements");
	match = true;
	}
	}

	/*! \page man_atom_weight fix_modify AtC atom_weight
	\section syntax
	fix_modify AtC atom_weight <method> <arguments>
	- <method> = \n
	value: atoms in specified group assigned constant value given \n
	lattice: volume per atom for specified lattice type (e.g. fcc) and parameter \n
	element: element volume divided among atoms within element \n
	region: volume per atom determined based on the atom count in the MD regions and their volumes. Note: meaningful only if atoms completely fill all the regions. \n
	group: volume per atom determined based on the atom count in a group and its volume\n
	read_in: list of values for atoms are read-in from specified file \n
	\section examples
	<TT> fix_modify atc atom_weight constant myatoms 11.8 </TT> \n
	<TT> fix_modify atc atom_weight lattice </TT> \n
	<TT> fix_modify atc atom_weight read-in atm_wt_file.txt </TT> \n
	\section description
	Command for assigning the value of atomic weights used for atomic integration in
	atom-continuum coupled simulations.
	\section restrictions
	Use of lattice option requires a lattice type and parameter is already specified.
	\section related
	\section default
	lattice
	*/
	else if (strcmp(arg[argIdx],"atom_weight")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"constant")==0) {
	argIdx++;
	atomWeightType_ = USER;
	int groupbit = -1;
	if (strcmp(arg[argIdx],"all")==0) {
	}
	else {
	groupbit = lammpsInterface_->group_bit(arg[argIdx]);
	}
	argIdx++;
	double value = atof(arg[argIdx]);
	Valpha_[groupbit] = value;
	match = true;
	}
	else if (strcmp(arg[argIdx],"lattice")==0) {
	atomWeightType_ = LATTICE;
	match = true;
	}
	else if (strcmp(arg[argIdx],"element")==0) {
	atomWeightType_ = ELEMENT;
	match = true;
	}
	else if (strcmp(arg[argIdx],"region")==0) {
	atomWeightType_ = REGION;
	match = true;
	}
	else if (strcmp(arg[argIdx],"group")==0) {
	atomWeightType_ = GROUP;
	match = true;
	}
	else if (strcmp(arg[argIdx],"multiscale")==0) {
	atomWeightType_ = MULTISCALE;
	match = true;
	}
	else if (strcmp(arg[argIdx],"node")==0) {
	atomWeightType_ = NODE;
	match = true;
	}
	else if (strcmp(arg[argIdx],"node_element")==0) {
	atomWeightType_ = NODE_ELEMENT;
	match = true;
	}
	else if (strcmp(arg[argIdx],"read_in")==0) {
	atomWeightType_ = READ_IN;
	argIdx++;
	atomicWeightsFile_ = arg[argIdx];
	match = true;
	}
	if (match) {
	needReset_ = true;
	}
	}

	/*! \page man_decomposition fix_modify AtC decomposition
	\section syntax
	fix_modify AtC decomposition <type>
	- <type> = \n
	replicated_memory: nodal information replicated on each processor \n
	distributed_memory: only owned nodal information on processor \n
	\section examples
	<TT> fix_modify atc decomposition distributed_memory </TT> \n
	\section description
	Command for assigning the distribution of work and memory for parallel runs.
	\section restrictions
	replicated_memory is appropriate for simulations were the number of nodes << number of atoms
	\section related
	\section default
	replicated_memory
	*/
	else if (strcmp(arg[argIdx],"decomposition")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"replicated_memory")==0) {
	domainDecomposition_ = REPLICATED_MEMORY;
	match = true;
	}
	else if (strcmp(arg[argIdx],"distributed_memory")==0) {
	domainDecomposition_ = DISTRIBUTED_MEMORY;
	match = true;
	}
	}

	/*! \page man_write_atom_weights fix_modify AtC write_atom_weights
	\section syntax
	fix_modify AtC write_atom_weights <filename> <frequency>
	- <filename> = name of file that atomic weights are written to \n
	- <frequency> = how often writes will occur \n
	\section examples
	<TT> fix_modify atc write_atom_weights atm_wt_file.txt 10 </TT> \n
	\section description
	Command for writing the values of atomic weights to a specified file.
	\section restrictions
	\section related
	\section default
	*/
	else if (strcmp(arg[argIdx],"write_atom_weights")==0) {
	argIdx++;
	atomicWeightsFile_ = arg[argIdx];
	argIdx++;
	atomicWeightsWriteFrequency_ = atoi(arg[argIdx]);
	atomicWeightsWriteFlag_ = true;
	match = true;
	}


	/*! \page man_reset_time fix_modify AtC reset_time
	\section syntax
	fix_modify AtC reset_time <value>
	\section examples
	<TT> fix_modify atc reset_time 0.0 </TT> \n
	\section description
	Resets the simulation time counter.
	\section restrictions
	\section related
	\section default
	*/
	else if (strcmp(arg[argIdx],"reset_time")==0) {
	argIdx++;
	set_time();
	if (narg > argIdx) {
	double time = atof(arg[argIdx]);
	set_time(time);
	}
	match = true;
	}

	/*! \page man_reset_time fix_modify AtC kernel_bandwidth
	\section syntax
	fix_modify AtC kernel_bandwidth <value>
	\section examples
	<TT> fix_modify atc kernel_bandwidth 8 </TT> \n
	\section description
	Sets a maximum parallel bandwidth for the kernel functions during parallel communication. If the command is not issued, the default will be to assume the bandwidth of the kernel matrix corresponds to the number of sampling locations.
	\section restrictions
	Only is used if kernel functions are being used.
	\section related
	\section default
	Number of sample locations.
	*/
	else if (strcmp(arg[argIdx],"kernel_bandwidth")==0) {
	argIdx++;
	accumulantBandwidth_ = atoi(arg[argIdx]);
	match = true;
	}

	/*! \page man_reset_atomic_reference_positions fix_modify AtC reset_atomic_reference_positions
	\section syntax
	fix_modify AtC reset_atomic_reference_positions
	\section examples
	<TT> fix_modify atc reset_atomic_reference_positions
	\section description
	Resets the atomic positions ATC uses to perform point to field operations.
	In can be used to use perfect lattice sites in ATC but a thermalized or
	deformed lattice in LAMMPS.
	\section restrictions
	\section related
	\section default
	Default is off
	*/
	else if (strcmp(arg[argIdx],"reset_atomic_reference_positions")==0) {
	argIdx++;
	xRefFile_ = arg[argIdx];
	readXref_ = true;
	match = true;
	}

	/*! \page man_set fix_modify AtC set
	\section syntax
	fix_modify AtC set reference_potential_energy <value_or_filename(optional)>
	- value (double) : optional user specified zero point for PE in native LAMMPS energy units \n
	- filename (string) : optional user specified string for file of nodal PE values to be read-in
	\section examples
	<TT> fix_modify AtC set reference_potential_energy </TT> \n
	<TT> fix_modify AtC set reference_potential_energy -0.05 </TT> \n
	<TT> fix_modify AtC set reference_potential_energy myPEvalues </TT> \n
	\section description
	Used to set various quantities for the post-processing algorithms.
	It sets the zero point for the potential energy density using
	the value provided for all nodes, or from the current
	configuration of the lattice if no value is provided, or
	values provided within the specified filename.
	\section restrictions
	Must be used with the hardy/field type of AtC fix
	( see \ref man_fix_atc )
	\section related
	\section default
	Defaults to lammps zero point i.e. isolated atoms
	*/
	else if (strcmp(arg[argIdx],"set")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"reference_potential_energy")==0) {
	argIdx++;
	setRefPE_ = true;
	if (narg > argIdx) {
	string a(arg[argIdx]);
	if (is_dbl(a)) {
	double value = atof(arg[argIdx]);
	refPEvalue_ = value;
	setRefPEvalue_ = true;
	}
	else {
	nodalRefPEfile_ = arg[argIdx];
	readRefPE_ = true;
	}
	}
	match = true;
	}
	} // end "set"



	/*! \page man_atom_element_map fix_modify AtC atom_element_map
	\section syntax
	fix_modify AtC atom_element_map <eulerian\|lagrangian> <frequency> \n
	- frequency (int) : frequency of updating atom-to-continuum maps based on the
	current configuration - only for eulerian
	\section examples
	<TT> fix_modify atc atom_element_map eulerian 100 </TT>
	\section description
	Changes frame of reference from eulerian to lagrangian and sets the
	frequency for which the map from atoms to elements is reformed and
	all the attendant data is recalculated.
	\section restrictions
	Cannot change map type after initialization.
	\section related
	\section default
	lagrangian
	*/
	else if (strcmp(arg[argIdx],"atom_element_map")==0) {
	argIdx++;
	if (strcmp(arg[argIdx],"eulerian")==0) {
	atomToElementMapType_ = EULERIAN;
	argIdx++;
	atomToElementMapFrequency_ = atoi(arg[argIdx]);
	}
	else {
	atomToElementMapType_ = LAGRANGIAN;
	atomToElementMapFrequency_ = 0;
	}
	match = true;
	needReset_ = true;
	}

	/*! \page man_read_restart fix_modify AtC read_restart
	\section syntax
	fix_modify AtC read_restart [file_name] \n
	\section examples
	<TT> fix_modify AtC read_restart ATC_state </TT> \n
	\section description
	Reads the current state of the fields from a named text-based restart
	file.
	\section restrictions
	The restart file only contains fields and their time derivatives.
	The reference positions of the atoms and the commands that initialize
	the fix are not saved e.g. an identical mesh containing the same atoms
	will have to be recreated.
	\section related
	see write_restart \ref man_write_restart
	\section default
	none
	*/
	else if (strcmp(arg[argIdx],"read_restart")==0) {
	argIdx++;
	restartFileName_ = arg[argIdx];
	useRestart_ = true;
	match = true;
	}

	/*! \page man_write_restart fix_modify AtC write_restart
	\section syntax
	fix_modify AtC write_restart [file_name] \n
	\section examples
	<TT> fix_modify AtC write_restart restart.mydata </TT> \n
	\section description
	Dumps the current state of the fields to a named text-based restart file.
	This done when the command is invoked and not repeated, unlike the
	similar lammps command.
	\section restrictions
	The restart file only contains fields and their time derivatives.
	The reference positions of the atoms and the commands that initialize
	the fix are not saved e.g. an identical mesh containing the same atoms
	will have to be recreated.
	\section related
	see read_restart \ref man_read_restart
	\section default
	none
	*/
	else if (strcmp(arg[argIdx],"write_restart")==0) {
	argIdx++;
	string restartFileName(arg[argIdx]);
	RESTART_LIST data;
	write_restart_data(restartFileName,data);
	match = true;
	}

	} // end else

	return match; // return to FixATC

	}

	//--------------------------------------------------
	// helper function for parser
	// handles : "displacement x" and "temperature" by indexing argIdx
	// for fluxes : only normal fluxes can be prescribed
	//--------------------------------------------------
	void ATC_Method::parse_field(/const/ char ** args, int & argIdx,
	FieldName & thisField)
	{
	string thisName = args[argIdx++];
	thisField = string_to_field(thisName);
	map<FieldName,int>::const_iterator iter = fieldSizes_.find(thisField);
	if (iter == fieldSizes_.end()) {
	throw ATC_Error("Bad field name: "+thisName);
	}
	}

	//--------------------------------------------------
	// helper function for parser
	// handles : "displacement x" and "temperature" by indexing argIdx
	// for fluxes : only normal fluxes can be prescribed
	//--------------------------------------------------
	void ATC_Method::parse_field(/const/ char ** args, int & argIdx,
	FieldName & thisField, int & thisIndex)
	{
	string thisName = args[argIdx++];
	if (args[argIdx] == NULL) {
	throw ATC_Error("Need to give field '"+thisName+"' more args");
	}
	thisField = string_to_field(thisName);
	map<FieldName,int>::const_iterator iter = fieldSizes_.find(thisField);
	if (iter == fieldSizes_.end()) {
	throw ATC_Error("Bad field name: "+thisName);
	}
	string thisDim = args[argIdx];
	thisIndex = 0;
	if (string_to_index(thisDim,thisIndex)) {
	if ( !( thisIndex < fieldSizes_[thisField]) ) {
	throw ATC_Error("Bad field dimension "+thisDim);
	}
	argIdx++;
	}
	}

	//-------------------------------------------------------------------
	// this sets the peratom output

	void ATC_Method::update_peratom_output()
	{
	perAtomArray_ = perAtomOutput_;
	// copy values
	for (int i = 0; i < lammpsInterface_->nlocal(); i++) {
	for (int j = 0; j < nsd_; j++) {
	perAtomOutput_[i][j] = xref_[i][j];
	}
	for (int j = nsd_; j < sizePerAtomCols_; j++) {
	perAtomOutput_[i][j] = 0;
	}
	}
	int indx = nsd_;
	if (atomVolume_->nRows() > 0) { // kernel Hardy does not compute these
	const DIAG_MAT & myAtomicWeights(atomVolume_->quantity());
	for (int i = 0; i < nLocal_; i++) {
	double wg = myAtomicWeights(i,i);
	if (wg > 0) {
	int ii = internalToAtom_(i);
	perAtomOutput_[ii][indx] = 1./wg;
	}
	}
	}
	}

	void ATC_Method::adjust_xref_pbc()
	{

	int nlocal = lammpsInterface_->nlocal();
	int xperiodic = lammpsInterface_->xperiodic();
	int yperiodic = lammpsInterface_->yperiodic();
	int zperiodic = lammpsInterface_->zperiodic();
	double **x = lammpsInterface_->xatom();
	double boxxlo,boxxhi;
	double boxylo,boxyhi;
	double boxzlo,boxzhi;

	lammpsInterface_->box_bounds(boxxlo,boxxhi,
	boxylo,boxyhi,
	boxzlo,boxzhi);
	// bool changed = false;
	for (int i = 0; i < nlocal; i++) {
	if (xperiodic) {
	if (x[i][0] < boxxlo) {
	xref_[i][0] += Xprd_;
	// changed = true;
	}
	if (x[i][0] >= boxxhi) {
	xref_[i][0] -= Xprd_;
	// changed = true;
	}
	}

	if (yperiodic) {
	if (x[i][1] < boxylo) {
	xref_[i][1] += Yprd_;
	// changed = true;
	}
	if (x[i][1] >= boxyhi) {
	xref_[i][1] -= Yprd_;
	// changed = true;
	}
	}

	if (zperiodic) {
	if (x[i][2] < boxzlo) {
	xref_[i][2] += Zprd_;
	// changed = true;
	}
	if (x[i][2] >= boxzhi) {
	xref_[i][2] -= Zprd_;
	// changed = true;
	}
	}
	}

	// propagate reset if needed
	if (atomToElementMapType_ == LAGRANGIAN) {
	if (atomCoarseGrainingPositions_) {
	atomCoarseGrainingPositions_->force_reset();
	}
	}
	else if (atomReferencePositions_) {
	atomReferencePositions_->force_reset();
	}

	}
	//-------------------------------------------------------------------
	void ATC_Method::initialize()
	{
	feEngine_->partition_mesh();
	// initialized_ is set to true by derived class initialize()
	// localStep_ is a counter within a run or minimize
	localStep_ = 0;
	// STEP 1) get basic information data from Lammps/fix
	// 1a) group ids for all internal atoms
	groupbit_ = lammpsInterface_->group_bit(groupTag_);

	// 1b) group ids for ghost atoms
	groupbitGhost_ = 0;
	if (!groupTagGhost_.empty()) {
	groupbitGhost_ = lammpsInterface_->group_bit(groupTagGhost_);
	}

	// 1c) periodicity and box bounds/lengths
	if (!initialized_) {

	lammpsInterface_->box_periodicity(periodicity[0],
	periodicity[1],
	periodicity[2]);
	lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0],
	box_bounds[0][1],box_bounds[1][1],
	box_bounds[0][2],box_bounds[1][2]);
	for (int k = 0; k < nsd_; k++) {
	box_length[k] = box_bounds[1][k] - box_bounds[0][k];
	}

	lammpsInterface_->set_reference_box();

	// get periodicity data from lammps for parallel exchange to adjust for periodicity
	Xprd_ = lammpsInterface_->domain_xprd();
	Yprd_ = lammpsInterface_->domain_yprd();
	Zprd_ = lammpsInterface_->domain_zprd();
	// box_length[0] = Xprd_;
	// box_length[1] = Yprd_;
	// box_length[2] = Zprd_;
	XY_ = lammpsInterface_->domain_xy();
	XZ_ = lammpsInterface_->domain_xz();
	YZ_ = lammpsInterface_->domain_yz();
	}

	// STEP 2 computational geometry
	// 2a) get basic information from continuum/FE
	this->set_continuum_data();

	// STEP 2b) set up data structures for computational geometry
	if (this->reset_methods()) {
	// clear memory manager
	interscaleManager_.clear_temporary_data();
	atomVolume_ = NULL;

	// reference positions and energy
	if (!initialized_) {
	double **x = lammpsInterface_->xatom();
	for (int i = 0; i < lammpsInterface_->nmax(); i++) {
	for (int j = 0; j < nsd_; j++) {
	xref_[i][j] = x[i][j];
	}
	}

	// re-write non-ghosts' xref with values from a file
	if (readXref_) {
	bool success = read_atomic_ref_positions(xRefFile_.c_str());
	if (!success)
	throw ATC_Error("Error reading atomic reference positions");
	readXref_ = false;
	}

	// ensure initial configuration is consistent with element set

	if (internalElementSet_.size() && groupbitGhost_) {
	int *mask = lammpsInterface_->atom_mask();
	int nlocal = lammpsInterface_->nlocal();
	const FE_Mesh * feMesh = feEngine_->fe_mesh();
	const set<int> & elementSet(feMesh->elementset(internalElementSet_));
	int element;
	DENS_VEC coords(nsd_);
	for (int i = 0; i < nlocal; ++i) {
	if (mask[i] & groupbit_ \|\| mask[i] & groupbitGhost_) {
	for (int j = 0; j < nsd_; j++) {
	coords(j) = xref_[i][j];
	}
	element = feMesh->map_to_element(coords);
	if (elementSet.find(element) == elementSet.end()) {
	mask[i] \|= groupbitGhost_;
	mask[i] &= ~groupbit_;
	}
	else {
	mask[i] &= ~groupbitGhost_;
	mask[i] \|= groupbit_;
	}
	}
	}
	}

	// set up maps from lammps to atc indexing
	reset_nlocal();
	}

	this->set_computational_geometry();
	}

	// 2c) basic data regarding atomic system, e.g. atom coordinates
	if (atomToElementMapType_ == EULERIAN) {
	reset_coordinates();
	}

	// STEP 3) set up variables which will be integrated in time
	this->construct_time_integration_data();

	// STEP 4) instantiate all the various specific algorithms and methods
	this->construct_methods();

	// STEP 5) construct dependency-managed data
	// 5b) all other transfer operators
	// needs to be done before every run in case options have changed or the atoms have been changed by the user

	if (this->reset_methods()) {
	// construct all the needed data structures
	this->construct_transfers();

	// allocate all space needed for lammps arrays
	interscaleManager_.grow_arrays(lammpsInterface_->nmax());
	}
	// reset all computes invoked flags and lammps data
	interscaleManager_.lammps_force_reset();

	// STEP 6) initialize data
	// 6b) size quantities which use pack_comm
	interscaleManager_.size_comm_quantities();

	// 6c) set coarse-graining functions and atomic weights
	if (!initialized_) {
	// FE_Engine allocates all required memory
	// assume initial atomic position is the reference position for now

	// \int_\Omega N_I dV : static if the mesh is
	NodeVolumes_.reset(nNodes_,nNodes_);
	invNodeVolumes_.reset(nNodes_,nNodes_);
	feEngine_->compute_lumped_mass_matrix(NodeVolumes_);
	invNodeVolumes_.set_quantity() = NodeVolumes_.inv();
	}
	atomVolume_->set_reset();

	// 6d) reference values
	this->set_reference_potential_energy();

	// 6e) atomic output for 0th step
	update_peratom_output();

	massMatInv_.reset(nNodes_,nNodes_);
	feEngine_->compute_lumped_mass_matrix(massMatInv_);
	for (int i = 0; i < nNodes_; ++i) {
	massMatInv_(i,i) = 1./massMatInv_(i,i);
	}

	// clear need for resets
	needReset_ = false;

	}
	//-------------------------------------------------------------------
	void ATC_Method::set_continuum_data()
	{
	// initialize finite element engine and get basic properties
	if (!initialized_) {
	feEngine_->initialize();
	if (nsd_!=feEngine_->nsd()) {
	throw ATC_Error("Spatial dimensions inconsistent between LAMMPS and ATC");
	}
	nNodes_ = feEngine_->num_nodes();
	}
	}

	//-------------------------------------------------------------------
	void ATC_Method::set_computational_geometry()
	{
	// set positions used for coarse-graining operators




	if (!initialized_) {
	if (atomToElementMapType_ == EULERIAN) {
	FundamentalAtomQuantity * atomPositionsAll = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,ALL);
	ClonedAtomQuantity<double> * myAtomPositions =
	new ClonedAtomQuantity<double>(this,atomPositionsAll,INTERNAL);
	atomCoarseGrainingPositions_ = myAtomPositions;
	interscaleManager_.add_per_atom_quantity(myAtomPositions,
	"AtomicCoarseGrainingPositions");

	if (trackDisplacement_) {
	XrefWrapper * myAtomReferencePositions = new XrefWrapper(this);
	atomReferencePositions_ = myAtomReferencePositions;
	interscaleManager_.add_per_atom_quantity(myAtomReferencePositions,
	"AtomicReferencePositions");
	atomReferencePositions_->set_memory_type(PERSISTENT);
	}

	if (groupbitGhost_) {
	myAtomPositions = new ClonedAtomQuantity<double>(this,atomPositionsAll,GHOST);
	atomGhostCoarseGrainingPositions_ = myAtomPositions;
	interscaleManager_.add_per_atom_quantity(myAtomPositions,
	"AtomicGhostCoarseGrainingPositions");
	}
	if(needProcGhost_){
	FundamentalAtomQuantity * atomPositionsAll = interscaleManager_.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION,PROC_GHOST);
	ClonedAtomQuantity<double> * myAtomPositions =
	new ClonedAtomQuantity<double>(this,atomPositionsAll,PROC_GHOST);
	atomProcGhostCoarseGrainingPositions_ = myAtomPositions;
	interscaleManager_.add_per_atom_quantity(myAtomPositions,
	"AtomicProcGhostCoarseGrainingPositions");
	}
	}
	else {
	XrefWrapper * myAtomPositions = new XrefWrapper(this);
	atomCoarseGrainingPositions_ = myAtomPositions;
	interscaleManager_.add_per_atom_quantity(myAtomPositions,
	"AtomicCoarseGrainingPositions");
	atomReferencePositions_ = atomCoarseGrainingPositions_;

	if (groupbitGhost_) {
	myAtomPositions = new XrefWrapper(this,GHOST);
	atomGhostCoarseGrainingPositions_ = myAtomPositions;
	interscaleManager_.add_per_atom_quantity(myAtomPositions,
	"AtomicGhostCoarseGrainingPositions");
	}
	if (needProcGhost_) {
	XrefWrapper * myAtomPositions = new XrefWrapper(this);
	atomProcGhostCoarseGrainingPositions_ = myAtomPositions;
	interscaleManager_.add_per_atom_quantity(myAtomPositions,
	"AtomicProcGhostCoarseGrainingPositions");
	}
	}
	atomCoarseGrainingPositions_->set_memory_type(PERSISTENT);
	if (atomGhostCoarseGrainingPositions_) atomGhostCoarseGrainingPositions_->set_memory_type(PERSISTENT);
	if (atomProcGhostCoarseGrainingPositions_) atomProcGhostCoarseGrainingPositions_->set_memory_type(PERSISTENT);
	}

	// Add in atom to element map if using shape functions
	if (needsAtomToElementMap_) {
	atomElement_ = new AtomToElementMap(this);
	interscaleManager_.add_per_atom_int_quantity(atomElement_,"AtomElement");
	}
	}

	//-------------------------------------------------------------------
	void ATC_Method::construct_methods()
	{

	if (this->reset_methods()) {
	if (atomTimeIntegrator_) delete atomTimeIntegrator_;
	if (integrateInternalAtoms_) {
	atomTimeIntegrator_ = new AtomTimeIntegratorType(this,INTERNAL);
	}
	else {
	atomTimeIntegrator_ = new AtomTimeIntegrator();
	}

	// set up integration schemes for ghosts
	ghostManager_.construct_methods();
	}
	}

	//-------------------------------------------------------------------
	void ATC_Method::construct_transfers()
	{
	this->construct_interpolant();

	this->construct_molecule_transfers();

	atomTimeIntegrator_->construct_transfers();
	ghostManager_.construct_transfers();



	}
	//-------------------------------------------------------------------
	PerAtomDiagonalMatrix<double> * ATC_Method::create_atom_volume()
	{
	if (atomVolume_) {
	return atomVolume_;
	}
	else {
	// set variables to compute atomic weights
	DENS_MAN * nodalVolume(NULL);
	switch (atomWeightType_) {
	case USER:
	atomVolume_ = new AtomVolumeUser(this,Valpha_);
	break;
	case LATTICE:
	atomVolume_ = new AtomVolumeLattice(this);
	break;
	case ELEMENT:
	atomVolume_ = new AtomVolumeElement(this);
	break;
	case REGION:
	atomVolume_ = new AtomVolumeRegion(this);
	break;
	case GROUP:
	atomVolume_ = new AtomVolumeGroup(this,Valpha_);
	break;
	case MULTISCALE:
	if (!shpFcn_) {
	throw ATC_Error("ATC_Method::create_atom_volume - Multiscale algorithm requires an interpolant");
	}
	nodalVolume = new NodalAtomVolume(this,shpFcn_);
	interscaleManager_.add_dense_matrix(nodalVolume,"NodalAtomVolume");
	atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_);
	break;
	case NODE:
	if (!shpFcn_) {
	throw ATC_Error("ATC_Method::create_atom_volume - Node algorithm requires an interpolant");
	}
	nodalVolume = new NodalVolume(this,shpFcn_);
	interscaleManager_.add_dense_matrix(nodalVolume,"NodalVolume");
	atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_);
	break;
	case NODE_ELEMENT:
	if (!shpFcn_) {
	throw ATC_Error("ATC_Method::create_atom_volume - Node-Element algorithm requires an interpolant");
	}
	nodalVolume = new NodalAtomVolumeElement(this,shpFcn_);
	interscaleManager_.add_dense_matrix(nodalVolume,"NodalAtomVolumeElement");
	atomVolume_ = new FtaShapeFunctionProlongationDiagonalMatrix(this,nodalVolume,shpFcn_);
	break;
	case READ_IN:
	atomVolume_ = new AtomVolumeFile(this,atomicWeightsFile_);
	break;
	}
	if (atomVolume_) {
	interscaleManager_.add_per_atom_diagonal_matrix(atomVolume_,"AtomVolume");
	}
	else {
	throw ATC_Error("ATC_Method::create_atom_volume - bad option for atom volume algorithm");
	}

	return atomVolume_;
	}
	}
	//--------------------------------------------------------
	void ATC_Method::init_integrate()
	{
	atomTimeIntegrator_->init_integrate_velocity(dt());
	ghostManager_.init_integrate_velocity(dt());
	// account for other fixes doing time integration
	interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY);

	atomTimeIntegrator_->init_integrate_position(dt());
	ghostManager_.init_integrate_position(dt());
	// account for other fixes doing time integration
	interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_POSITION);
	}
	//-------------------------------------------------------------------
	void ATC_Method::post_init_integrate()
	{
	ghostManager_.post_init_integrate();
	}
	//-------------------------------------------------------------------
	void ATC_Method::pre_exchange()
	{
	adjust_xref_pbc();
	// call interscale manager to sync atc per-atom data with lammps array ahead of parallel communication
	interscaleManager_.prepare_exchange();

	// change types based on moving from internal region to ghost region
	if ((atomToElementMapType_ == EULERIAN) && (step() % atomToElementMapFrequency_ == 0)) {
	ghostManager_.pre_exchange();
	}
	}
	//-------------------------------------------------------------------
	void ATC_Method::setup_pre_exchange()
	{
	adjust_xref_pbc();
	// call interscale manager to sync atc per-atom data with lammps array ahead of parallel communication
	interscaleManager_.prepare_exchange();
	}
	//-------------------------------------------------------------------
	void ATC_Method::pre_neighbor()
	{
	// reset quantities arising from atom exchange
	reset_nlocal();

	interscaleManager_.post_exchange();

	// forward_comm should go here
	}
	//-------------------------------------------------------------------
	void ATC_Method::min_post_force()
	{
	post_force();
	}
	//-------------------------------------------------------------------
	void ATC_Method::post_force()
	{
	// this resets allow for the possibility of other fixes modifying positions and velocities, e.g. walls, but reduces efficiency
	interscaleManager_.lammps_force_reset();
	}
	//--------------------------------------------------------
	void ATC_Method::final_integrate()
	{
	atomTimeIntegrator_->final_integrate(dt());
	ghostManager_.final_integrate(dt());
	// account for other fixes doing time integration
	interscaleManager_.fundamental_force_reset(LammpsInterface::ATOM_VELOCITY);
	}
	//-------------------------------------------------------------------
	void ATC_Method::post_final_integrate()
	{
	if (atomicWeightsWriteFlag_ && (step() % atomicWeightsWriteFrequency_ == 0)) {
	write_atomic_weights(atomicWeightsFile_,atomVolume_->quantity());
	}
	}
	//-------------------------------------------------------------------
	void ATC_Method::end_of_step()
	{
	localStep_ += 1;
	}
	//--------------------------------------------------------------
	void ATC_Method::finish()
	{
	// FE Engine
	if (feEngine_) feEngine_->finish();
	feEngine_->departition_mesh();
	}

	//--------------------------------------------------------------
	/** method to add new fields to the included list */
	//--------------------------------------------------------------
	void ATC_Method::add_fields(map<FieldName,int> & newFieldSizes)
	{
	map<FieldName,int>::const_iterator field;
	for (field = newFieldSizes.begin(); field!=newFieldSizes.end(); field++) {
	FieldName thisField = field->first;
	int thisSize = field->second;
	if (fieldSizes_.find(thisField)==fieldSizes_.end()) {
	fieldSizes_[thisField] = thisSize;
	}
	}
	}

	//-------------------------------------------------------------------
	void ATC_Method::set_reference_potential_energy(void)
	{
	if (setRefPE_) {
	if (setRefPEvalue_) {
	nodalRefPotentialEnergy_->set_quantity() = refPEvalue_;
	setRefPEvalue_ = false;
	}
	else if (readRefPE_) {
	if (LammpsInterface::instance()->rank_zero()) {
	stringstream ss;
	ss << "reading reference potential energy from " << nodalRefPEfile_;
	LammpsInterface::instance()->print_msg(ss.str());
	}
	(nodalRefPotentialEnergy_->set_quantity()).from_file(nodalRefPEfile_);
	readRefPE_ = false;
	}
	else {
	hasRefPE_ = false;
	SPAR_MAN * referenceAccumulant = interscaleManager_.sparse_matrix("ReferenceAccumulant");
	if (referenceAccumulant) {
	referenceAccumulant->set_quantity() = accumulant_->quantity();
	}
	DIAG_MAN * referenceAccumulantInverseVolumes = interscaleManager_.diagonal_matrix("ReferenceAccumulantInverseVolumes");
	if (referenceAccumulantInverseVolumes) {
	referenceAccumulantInverseVolumes->set_quantity() = accumulantInverseVolumes_->quantity();
	}
	PAQ * atomicRefPe = interscaleManager_.per_atom_quantity("AtomicReferencePotential");
	if (!atomicRefPe) {
	throw ATC_Error("ATC_Method::set_reference_potential_energy - atomic reference PE object was not created during construct_transfers");
	}
	PAQ* pe = interscaleManager_.per_atom_quantity("AtomicPotentialEnergy");
	if (!pe) {
	throw ATC_Error("ATC_Method::set_reference_potential_energy - atomic PE object was not created during construct_transfers");
	}
	atomicRefPe->set_quantity() = pe->quantity();
	atomicRefPe->fix_quantity();
	}
	setRefPE_ = false;
	hasRefPE_ = true;
	}
	}
	//-------------------------------------------------------------------


	//=================================================================
	// memory management and processor information exchange
	//=================================================================


	//-----------------------------------------------------------------
	// number of doubles
	//-----------------------------------------------------------------
	int ATC_Method::doubles_per_atom() const
	{

	int doubles = 4;
	doubles += interscaleManager_.memory_usage();
	return doubles;
	}

	//-----------------------------------------------------------------
	// memory usage of local atom-based arrays
	//-----------------------------------------------------------------
	int ATC_Method::memory_usage()
	{
	int bytes = doubles_per_atom();
	bytes = lammpsInterface_->nmax() sizeof(double);
	return bytes;
	}

	//-----------------------------------------------------------------
	// allocate local atom-based arrays
	//-----------------------------------------------------------------
	void ATC_Method::grow_arrays(int nmax)
	{
	xref_ =
	lammpsInterface_->grow_2d_double_array(xref_,nmax,3,"fix_atc:xref");

	perAtomOutput_ =
	lammpsInterface_->grow_2d_double_array(perAtomOutput_,nmax,sizePerAtomCols_,"fix_atc:perAtomOutput");
	interscaleManager_.grow_arrays(nmax);
	}

	//-----------------------------------------------------------------
	// copy values within local atom-based arrays
	//-----------------------------------------------------------------
	void ATC_Method::copy_arrays(int i, int j)
	{
	xref_[j][0] = xref_[i][0];
	xref_[j][1] = xref_[i][1];
	xref_[j][2] = xref_[i][2];

	for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) {
	perAtomOutput_[j][ii] = perAtomOutput_[i][ii];
	}
	interscaleManager_.copy_arrays(i,j);
	}

	//-----------------------------------------------------------------
	// pack values in local atom-based arrays for exchange with another proc
	//-----------------------------------------------------------------
	int ATC_Method::pack_exchange(int i, double *buf)
	{
	buf[0] = xref_[i][0];
	buf[1] = xref_[i][1];
	buf[2] = xref_[i][2];

	int j = 4;
	for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) {
	buf[j++] = perAtomOutput_[i][ii];
	}
	int interscaleSizeComm = interscaleManager_.pack_exchange(i,&buf[j]);
	return sizeComm_ + interscaleSizeComm;
	}

	//-----------------------------------------------------------------
	// unpack values in local atom-based arrays from exchange with another proc
	//-----------------------------------------------------------------
	int ATC_Method::unpack_exchange(int nlocal, double *buf)
	{
	xref_[nlocal][0] = buf[0];
	xref_[nlocal][1] = buf[1];
	xref_[nlocal][2] = buf[2];

	int j = 4;
	for (int ii = 0 ; ii < sizePerAtomCols_ ; ii++ ) {
	perAtomOutput_[nlocal][ii] = buf[j++];
	}
	int interscaleSizeComm = interscaleManager_.unpack_exchange(nlocal,&buf[j]);
	return sizeComm_ + interscaleSizeComm;
	}

	//-----------------------------------------------------------------
	// pack values in local atom-based arrays from exchange with another proc
	//-----------------------------------------------------------------
	int ATC_Method::pack_comm(int n, int list, double buf,
	int pbc_flag, int *pbc)
	{
	int i,j,m;
	double dx = 0,dy = 0,dz = 0;

	int * num_bond = lammpsInterface_->num_bond();
	int ** bond_atom = lammpsInterface_->bond_atom();

	m = 0;
	if (pbc_flag == 0) {
	for (i = 0; i < n; i++) {
	j = list[i];
	buf[m++] = xref_[j][0];
	buf[m++] = xref_[j][1];
	buf[m++] = xref_[j][2];

	if (num_bond) {
	buf[m++] = num_bond[j];
	for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) {
	buf[m++] = bond_atom[j][ii];
	}
	}
	}
	}
	else {
	if (lammpsInterface_->domain_triclinic() == 0) {
	dx = pbc[0]*Xprd_;
	dy = pbc[1]*Yprd_;
	dz = pbc[2]*Zprd_;
	}
	else {
	dx = pbc[0]Xprd_ + pbc[5]XY_ + pbc[4]*XZ_;
	dy = pbc[1]Yprd_ + pbc[3]YZ_;
	dz = pbc[2]*Zprd_;
	}
	for (i = 0; i < n; i++) {
	j = list[i];
	buf[m++] = xref_[j][0] + dx;
	buf[m++] = xref_[j][1] + dy;
	buf[m++] = xref_[j][2] + dz;

	if (num_bond) {
	buf[m++] = num_bond[j];
	for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) {
	buf[m++] = bond_atom[j][ii];
	}
	}

	}
	}
	return m; // total amount of data sent
	}

	//-----------------------------------------------------------------
	// unpack values in local atom-based arrays from exchange with another proc
	//-----------------------------------------------------------------
	void ATC_Method::unpack_comm(int n, int first, double *buf)
	{
	int i,m,last;

	int * num_bond = lammpsInterface_->num_bond();
	int ** bond_atom = lammpsInterface_->bond_atom();

	m = 0;
	last = first + n;
	for (i = first; i < last; i++) {
	xref_[i][0] = buf[m++];
	xref_[i][1] = buf[m++];
	xref_[i][2] = buf[m++];

	if (num_bond) {
	num_bond[i] = static_cast<int>(buf[m++]);
	for (int ii = 0; ii < lammpsInterface_->bond_per_atom(); ii++) {
	bond_atom[i][ii] = static_cast<int>(buf[m++]);
	}
	}

	}
	}

	//-----------------------------------------------------------------
	//
	//-----------------------------------------------------------------
	int ATC_Method::comm_forward()
	{
	int size = 3;
	if (lammpsInterface_->num_bond())
	{ size += lammpsInterface_->bond_per_atom()+1; }
	return size;
	}

	//-----------------------------------------------------------------
	//
	//-----------------------------------------------------------------
	void ATC_Method::reset_nlocal()
	{
	nLocalTotal_ = lammpsInterface_->nlocal();
	const int * mask = lammpsInterface_->atom_mask();
	nLocal_ = 0;
	nLocalGhost_ = 0;

	for (int i = 0; i < nLocalTotal_; ++i) {
	if (mask[i] & groupbit_) nLocal_++;
	if (mask[i] & groupbitGhost_) nLocalGhost_++;
	}

	// set up internal & ghost maps

	if (nLocal_>0) {
	// set map
	internalToAtom_.resize(nLocal_);
	int j = 0;
	// construct internalToAtom map
	// : internal index -> local lammps atom index
	for (int i = 0; i < nLocalTotal_; ++i) {
	if (mask[i] & groupbit_) internalToAtom_(j++) = i;
	}
	#ifdef EXTENDED_ERROR_CHECKING
	stringstream ss;
	ss << "Nlocal = " << nLocal_ << " but only found " << j << "atoms";
	if (j!=nLocal_) throw ATC_Error(ss.str());
	#endif
	// construct reverse map
	atomToInternal_.clear();
	for (int i = 0; i < nLocal_; ++i) {
	atomToInternal_[internalToAtom_(i)] = i;
	}
	}
	if (nLocalGhost_>0) {
	// set map
	ghostToAtom_.resize(nLocalGhost_);
	int j = 0;
	for (int i = 0; i < nLocalTotal_; ++i) {
	if (mask[i] & groupbitGhost_) ghostToAtom_(j++) = i;
	}
	}

	//WIP_JAT this should not be needed at all, but a memory problem with sparse matrices requires them to be reset (possibly related to note in SparseMatrix-inl.h::_delete())
	interscaleManager_.reset_nlocal();

	}

	//-------------------------------------------------------------------
	void ATC_Method::reset_coordinates()
	{
	// update coarse graining positions for internal and ghost atoms
	atomCoarseGrainingPositions_->unfix_quantity();
	atomCoarseGrainingPositions_->quantity();
	atomCoarseGrainingPositions_->fix_quantity();
	if (atomGhostCoarseGrainingPositions_) {
	atomGhostCoarseGrainingPositions_->unfix_quantity();
	atomGhostCoarseGrainingPositions_->quantity();
	atomGhostCoarseGrainingPositions_->fix_quantity();
	}
	if (atomProcGhostCoarseGrainingPositions_) {
	atomProcGhostCoarseGrainingPositions_->unfix_quantity();
	atomProcGhostCoarseGrainingPositions_->quantity();
	atomProcGhostCoarseGrainingPositions_->fix_quantity();
	}
	}

	//-----------------------------------------------------------------
	//
	//-----------------------------------------------------------------
	void ATC_Method::write_atomic_weights(const string filename, const DIAG_MAT & atomicVolumeMatrix)
	{
	int nlocal = lammpsInterface_->nlocal();
	int nlocalmax;
	LammpsInterface::instance()->int_allmax(&nlocal,&nlocalmax);
	int natoms = int(lammpsInterface_->natoms());
	ofstream out;
	const char* fname = &filename[0];

	// create tag to local id map for this processor
	map <int,int> id2tag;
	map <int,int>::const_iterator itr;
	int * atom_tag = lammpsInterface_->atom_tag();
	for (int i = 0; i < nlocal; ++i) {
	id2tag[i] = atom_tag[i];
	}

	int comm_rank = LammpsInterface::instance()->comm_rank();
	int nprocs;
	LammpsInterface::instance()->int_allmax(&comm_rank,&nprocs);
	nprocs += 1;

	if (comm_rank == 0) {
	out.open(fname);
	// print header lines
	out << "Atomic Weights for LAMMPS/atc analysis\n";
	out << " \n";
	out << natoms << " Atoms in system\n";
	out << " \n";
	// print atomic weights from proc 0
	for(int i = 0; i < nlocal; i++) {
	out << id2tag[i] << " " << atomicVolumeMatrix(i,i) << "\n";
	}
	}

	if (nprocs > 1) {
	int max_size,send_size;
	send_size = nlocal;
	LammpsInterface::instance()->int_allmax(&send_size,&max_size);

	if (comm_rank == 0) {
	int intbuf[max_size];
	double buf[max_size];
	for (int iproc = 1; iproc < nprocs; iproc++) {
	LammpsInterface::instance()->int_recv(intbuf,max_size,iproc);
	LammpsInterface::instance()->recv(buf,max_size,iproc);
	for (int i = 0; i < max_size; i++) {
	out << intbuf[i] << " " << buf[i] << "\n";
	}
	}
	} else {
	int intbuf[send_size];
	double buf[send_size];
	for (int i = 0; i < send_size; i++) {
	intbuf[i] = id2tag[i];
	buf[i] = atomicVolumeMatrix(i,i);
	}
	LammpsInterface::instance()->int_send(intbuf,send_size);
	LammpsInterface::instance()->send(buf,send_size);
	}
	}

	if (comm_rank == 0) {
	out.close();
	}
	}

	//-----------------------------------------------------------------
	//
	//-----------------------------------------------------------------
	void ATC_Method::compute_consistent_md_mass_matrix(const SPAR_MAT & shapeFunctionMatrix,
	SPAR_MAT & mdMassMatrix) const
	{

	int nCols = shapeFunctionMatrix.nCols();
	DENS_MAT massMatrixLocal(nCols,nCols);
	DENS_MAT denseMdMassMatrix(nCols,nCols);
	if (nLocal_>0)
	massMatrixLocal = shapeFunctionMatrix.transMat(shapeFunctionMatrix);

	lammpsInterface_->allsum(massMatrixLocal.ptr(),
	denseMdMassMatrix.ptr(),
	denseMdMassMatrix.size());
	mdMassMatrix.reset(denseMdMassMatrix,1.e-10);
	}

	//=================================================================
	// Interscale operators
	//=================================================================
	// in the spirit of the current design of ATC: atoms local, nodes global




	bool ATC_Method::nodal_influence(const int groupbit,
	set<int> & nset, set<int> & aset, double tol)
	{
	int nghost = nodal_influence(groupbit,nset,aset,true,tol);
	int local_nghost = nghost;
	lammpsInterface_->int_allsum(&local_nghost,&nghost);
	if (nghost == 0) {
	nodal_influence(groupbit,nset,aset,false,tol);
	}
	return (nghost > 0);
	}
	int ATC_Method::nodal_influence(const int groupbit,
	set<int> & nset, set<int> & aset, bool ghost, double tol)
	{
	Array<int> & amap = (ghost) ? ghostToAtom_ : internalToAtom_;
	int natoms = (ghost) ? nLocalGhost_ : nLocal_;
	DENS_MAT influence(nNodes_,1);
	DENS_MAT atomInfluence(natoms,1);
	const int *mask = lammpsInterface_->atom_mask();
	for (int i = 0; i < natoms; i++) {
	if (mask[amap(i)] & groupbit) {
	atomInfluence(i,0) = 1;
	aset.insert(i);
	}
	}
	// relies on shape functions
	if (ghost) {
	restrict_volumetric_quantity(atomInfluence,influence,(interscaleManager_.per_atom_sparse_matrix("InterpolantGhost"))->quantity());
	}
	else {
	restrict_volumetric_quantity(atomInfluence,influence);
	}

	DENS_MAT localInfluence = influence;
	lammpsInterface_->allsum(localInfluence.ptr(),
	influence.ptr(),
	influence.size());

	for (int i = 0; i < nNodes_; i++) {
	if (fabs(influence(i,0)) > tol) { nset.insert(i); }
	}
	return aset.size();
	}


	//--------------------------------------------------------
	void ATC_Method::restrict_volumetric_quantity(const MATRIX & atomData,
	MATRIX & nodeData,
	const SPAR_MAT & shpFcn)
	{
	// computes nodeData = NDeltaVAtomatomData where N are the shape functions
	DENS_MAT workNodeArray(nodeData.nRows(),nodeData.nCols());
	//DENS_MAT workNodeArray;


	if (atomData.nRows() > 0) { // or shpFcn_???
	workNodeArray = shpFcn.transMat(atomData);
	}
	int count = nodeData.nRows()*nodeData.nCols();
	lammpsInterface_->allsum(workNodeArray.ptr(),nodeData.ptr(),count);
	return;
	}


	//--------------------------------------------------------
	void ATC_Method::restrict_volumetric_quantity(const MATRIX & atomData,
	MATRIX & nodeData)
	{
	restrict_volumetric_quantity(atomData,nodeData,shpFcn_->quantity());
	return;
	}


	//--------------------------------------------------------
	void ATC_Method::prolong(const MATRIX & nodeData,
	MATRIX & atomData)
	{
	// computes the finite element interpolation at atoms atomData = N*nodeData
	if (nLocal_>0) {
	atomData = (shpFcn_->quantity())*nodeData;
	}
	return;
	}

	//========================================================
	// FE functions
	//========================================================

	//--------------------------------------------------------
	void ATC_Method::output()
	{
	// if (lammpsInterface_->comm_rank() == 0) {
	compute_nodeset_output();
	compute_faceset_output();
	compute_elementset_output();
	// }
	}
	//--------------------------------------------------------
	void ATC_Method::compute_nodeset_output(void)
	{
	map< pair <string, FieldName>, NodesetOperationType >::const_iterator iter;
	for (iter = nsetData_.begin(); iter != nsetData_.end();iter++){
	pair <string, FieldName> id = iter->first;
	string nsetName = id.first;
	FieldName field = id.second;
	double sum = 0.0;
	const set<int> nset = feEngine_->fe_mesh()->nodeset(nsetName);
	const DENS_MAT & thisField = (fields_.find(field)->second).quantity();
	set< int >::const_iterator itr;
	for (itr = nset.begin(); itr != nset.end();itr++){
	int node = *itr;
	sum += thisField(node,0);
	}
	string name = nsetName + "_" + field_to_string(field);
	if (iter->second == NODESET_AVERAGE) {
	sum /= nset.size();
	name = "average_"+name;
	}
	feEngine_->add_global(name, sum);
	}
	}
	//--------------------------------------------------------
	void ATC_Method::compute_faceset_output(void)
	{
	map < pair<string,string>, FacesetIntegralType >::const_iterator iter;
	DENS_MAT values;
	for (iter = fsetData_.begin(); iter != fsetData_.end(); iter++) {
	string bndyName = (iter->first).first;
	string fieldName = (iter->first).second;
	const set< PAIR > & faceSet = (feEngine_->fe_mesh())->faceset(bndyName);
	ATOMIC_DATA::iterator data_iter = filteredData_.find(fieldName);
	if (data_iter == filteredData_.end()) {
	string msg = "Specified fieldName "+fieldName+
	" not found in filteredData_ while attempting surface integration";
	throw ATC_Error(msg);
	}
	const DENS_MAT & data = ((data_iter->second).quantity());
	string stem = bndyName + "_" + fieldName + "_";
	bool tf = false;
	if (iter->second == CONTOUR_INTEGRAL) {
	stem = "contour_"+stem;
	tf = true;
	}
	feEngine_->field_surface_flux(data,faceSet,values,tf);
	for (int i = 0; i < values.nRows() ; ++i ) {
	string name = stem + to_string(i+1);
	feEngine_->add_global(name, values(i,0));
	}
	}
	}
	//--------------------------------------------------------
	void ATC_Method::compute_elementset_output(void)
	{
	map< pair <string, FieldName>, ElementsetOperationType >::const_iterator iter;
	for (iter = esetData_.begin(); iter != esetData_.end();iter++){
	pair <string, FieldName> id = iter->first;
	string esetName = id.first;
	FieldName field = id.second;
	const ESET eset = feEngine_->fe_mesh()->elementset(esetName);
	const DENS_MAT & thisField = (fields_.find(field)->second).quantity();
	DENS_VEC total = feEngine_->integrate(thisField,eset);
	string name = esetName + "_" + field_to_string(field);
	if (iter->second == ELEMENTSET_AVERAGE) {
	DENS_MAT ones(nNodes_,0); ones = 1;
	DENS_VEC V = feEngine_->integrate(ones,eset);
	total /= V[0];
	name = "average_"+name;
	}
	if (total.size() == 1) {
	feEngine_->add_global(name, total[0]);
	}
	else {
	for (int i = 0; i < total.size(); i++) {
	feEngine_->add_global(name+to_string(i), total[i]);
	}
	}
	}
	}


	//=================================================================
	//
	//=================================================================
	//--------------------------------------------------------
	bool ATC_Method::read_atomic_ref_positions(const char * filename)
	{
	int nlocal = lammpsInterface_->nlocal();
	ifstream in;
	const int lineSize = 256;
	char line[lineSize];

	// create tag to local id map for this processor
	map <int,int> tag2id;
	map <int,int>::const_iterator itr;
	int * atom_tag = lammpsInterface_->atom_tag();
	for (int i = 0; i < nlocal; ++i) {
	tag2id[atom_tag[i]] = i;
	}

	// get number of atoms
	int natoms = 0;
	int ncols = 0;
	int style = LammpsInterface::CHARGE_STYLE;
	if (LammpsInterface::instance()->rank_zero()) {
	in.open(filename);
	string msg;
	string name = filename;
	msg = "no "+name+" reference file found";
	if (! in.good()) throw ATC_Error(msg);
	in.getline(line,lineSize); // header
	in.getline(line,lineSize); // blank line
	in >> natoms;
	stringstream ss;
	ss << "found " << natoms << " atoms in reference " << filename ;
	while(in.good()) {
	in.getline(line,lineSize);
	string str(line);
	int pos = str.find("Atoms");
	if (pos > -1) {
	in.getline(line,lineSize); // blank line
	break;
	}
	}
	in.getline(line,lineSize);
	std::vector<std::string> tokens;
	ATC_Utility::command_strings(line, tokens);
	ncols = tokens.size();
	switch (ncols) {
	// atomic: id type x y z
	case 5:
	case 8:
	ss << " style:atomic";
	style = LammpsInterface::ATOMIC_STYLE;
	break;
	// charge: id type q x y z
	// molecule : id molecule-ID type x y z
	case 6:
	ss << " style:charge";
	style = LammpsInterface::CHARGE_STYLE;
	break;
	// full : id molecule-ID type q x y z
	case 7:
	ss << " style:full";
	style = LammpsInterface::FULL_STYLE;
	break;
	default:
	throw ATC_Error("cannot determine atom style, columns:"+to_string(ncols));
	break;
	}
	LammpsInterface::instance()->print_msg(ss.str());
	in.close();
	}
	LammpsInterface::instance()->int_broadcast(&natoms);

	// read atoms and assign
	if (LammpsInterface::instance()->rank_zero()) {
	in.open(filename);
	while(in.good()) {
	in.getline(line,lineSize);
	string str(line);
	int pos = str.find("Atoms");
	if (pos > -1) {
	in.getline(line,lineSize); // blank line
	break;
	}
	}
	}
	int nread = 0,type = -1, tag = -1, count = 0, mId = -1;
	double x[3]={0,0,0}, q =0;
	while (nread < natoms) {
	if (LammpsInterface::instance()->rank_zero()) {
	in.getline(line,lineSize);
	stringstream ss (line,stringstream::in \| stringstream::out);
	if (style == LammpsInterface::CHARGE_STYLE)
	ss >> tag >> type >> q >> x[0] >> x[1] >> x[2];
	else if (style == LammpsInterface::FULL_STYLE)
	ss >> tag >> mId >> type >> q >> x[0] >> x[1] >> x[2];
	else
	ss >> tag >> type >> x[0] >> x[1] >> x[2];
	nread++;
	}
	LammpsInterface::instance()->int_broadcast(&nread);
	LammpsInterface::instance()->int_broadcast(&tag);
	LammpsInterface::instance()->broadcast(x,3);
	itr = tag2id.find(tag);
	if (itr != tag2id.end()) {
	int iatom = itr->second;
	xref_[iatom][0] = x[0];
	xref_[iatom][1] = x[1];
	xref_[iatom][2] = x[2];
	count++;
	}
	}
	if (LammpsInterface::instance()->rank_zero()) {
	in.close();
	stringstream ss;
	ss << "read " << natoms << " reference positions";
	LammpsInterface::instance()->print_msg(ss.str());
	}
	if (count != nlocal)
	throw ATC_Error("reset "+to_string(count)+" atoms vs "+to_string(nlocal));


	return true;
	}

	//--------------------------------------------------------
	void ATC_Method::remap_ghost_ref_positions(void)
	{

	int nlocal = lammpsInterface_->nlocal();
	int nghost = lammpsInterface_->nghost();

	double box_bounds[2][3];
	lammpsInterface_->box_bounds(box_bounds[0][0],box_bounds[1][0],
	box_bounds[0][1],box_bounds[1][1],
	box_bounds[0][2],box_bounds[1][2]);
	double xlo = box_bounds[0][0], xhi = box_bounds[1][0];
	double ylo = box_bounds[0][1], yhi = box_bounds[1][1];
	double zlo = box_bounds[0][2], zhi = box_bounds[1][2];

	double box_length[3];
	for (int k = 0; k < 3; k++) {
	box_length[k] = box_bounds[1][k] - box_bounds[0][k];
	}
	double Lx = box_length[0], Ly = box_length[1], Lz = box_length[2];

	// create tag to local id map for this processor
	map <int,int> tag2id;
	map <int,int>::const_iterator itr;
	int * atom_tag = lammpsInterface_->atom_tag();
	for (int i = 0; i < nlocal; ++i) {
	tag2id[atom_tag[i]] = i;
	}

	// loop over ghosts
	double ** x = lammpsInterface_->xatom();
	for (int j = nlocal; j < nlocal + nghost; j++) {
	int tag = atom_tag[j];
	int i = tag2id[tag];
	//itr = tag2id.find(tag);
	//if (itr != tag2id.end())
	double* xj = x[j];
	double* Xj = xref_[j];
	//double Xj[3];
	double* Xi = xref_[i];
	// the assumption is that xref_[j] has been shuffled
	// so make an image of xref_[i] that is close to x[j]
	if (xj[0] <= xlo) Xj[0] = Xi[0] -Lx;
	if (xj[0] >= xhi) Xj[0] = Xi[0] +Lx;
	if (xj[1] <= ylo) Xj[1] = Xi[1] -Ly;
	if (xj[1] >= yhi) Xj[1] = Xi[1] +Ly;
	if (xj[2] <= zlo) Xj[2] = Xi[2] -Lz;
	if (xj[2] >= zhi) Xj[2] = Xi[2] +Lz;
	}
	}
	};

ATC_Method.cppNo OneTemporaryActions

File Metadata

ATC_Method.cppView Options

Event Timeline

ATC_Method.cpp
No OneTemporary
Actions

ATC_Method.cpp
View Options