Diffusion lammps (efficient_neuronet) Browse

ATC_Method.cpp
6eca661e4a54efficient_neuronet

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

lammps/lib/atc/ATC_Method.cpp6eca661e4a54efficient_neuronet

ATC_Method.cppView Options

lammps/lib/atc/ATC_Method.cpp
6eca661e4a54efficient_neuronet

ATC_Method.cpp
View Options