Thermostat.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Jul 31, 17:18

Thermostat.cpp
View Options

	#include "Thermostat.h"
	#include "ATC_Coupling.h"
	#include "ATC_Error.h"
	#include "PrescribedDataManager.h"
	#include "ThermalTimeIntegrator.h"
	#include "TransferOperator.h"

	using namespace std;

	namespace ATC {

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class Thermostat
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	Thermostat::Thermostat(ATC_Coupling * atc,
	const string & regulatorPrefix) :
	AtomicRegulator(atc,regulatorPrefix),
	lambdaMaxIterations_(myLambdaMaxIterations)
	{
	// nothing to do
	}

	//--------------------------------------------------------
	// modify:
	// parses and adjusts thermostat state based on
	// user input, in the style of LAMMPS user input
	//--------------------------------------------------------
	bool Thermostat::modify(int narg, char **arg)
	{
	bool foundMatch = false;

	int argIndex = 0;
	if (strcmp(arg[argIndex],"thermal")==0) {
	argIndex++;

	// thermostat type
	/*! \page man_control_thermal fix_modify AtC control thermal
	\section syntax
	fix_modify AtC control thermal <control_type> <optional_args>
	- control_type (string) = none \| rescale \| hoover \| flux\n

	fix_modify AtC control thermal rescale <frequency> \n
	- frequency (int) = time step frequency for applying velocity rescaling \n

	fix_modify AtC control thermal hoover \n

	fix_modify AtC control thermal flux <boundary_integration_type(optional)> <face_set_id(optional)>\n
	- boundary_integration_type (string) = faceset \| interpolate\n
	- face_set_id (string), optional = id of boundary face set, if not specified
	(or not possible when the atomic domain does not line up with
	mesh boundaries) defaults to an atomic-quadrature approximate
	evaulation, does not work with interpolate\n
	\section examples
	<TT> fix_modify AtC control thermal none </TT> \n
	<TT> fix_modify AtC control thermal rescale 10 </TT> \n
	<TT> fix_modify AtC control thermal hoover </TT> \n
	<TT> fix_modify AtC control thermal flux </TT> \n
	<TT> fix_modify AtC control thermal flux faceset bndy_faces </TT> \n
	\section description
	Sets the energy exchange mechansim from the finite elements to the atoms, managed through a control algorithm. Rescale computes a scale factor for each atom to match the finite element temperature. Hoover is a Gaussian least-constraint isokinetic thermostat enforces that the nodal restricted atomic temperature matches the finite element temperature. Flux is a similar mode, but rather adds energy to the atoms based on conservation of energy. Hoover and flux allows the prescription of sources or fixed temperatures on the atoms.
	\section restrictions
	only for be used with specific transfers :
	thermal (rescale, hoover, flux), two_temperature (flux) \n
	rescale not valid with time filtering activated
	\section related
	\section default
	none\n
	rescale frequency is 1\n
	flux boundary_integration_type is interpolate
	*/
	if (strcmp(arg[argIndex],"none")==0) { // restore defaults
	regulatorTarget_ = NONE;
	couplingMode_ = UNCOUPLED;
	howOften_ = 1;
	boundaryIntegrationType_ = NO_QUADRATURE;
	foundMatch = true;
	}
	else if (strcmp(arg[argIndex],"rescale")==0) {
	argIndex++;
	howOften_ = atoi(arg[argIndex]);
	if (howOften_ < 1) {
	throw ATC_Error("Bad rescaling thermostat frequency");
	}
	else {
	regulatorTarget_ = FIELD;
	couplingMode_ = UNCOUPLED;
	boundaryIntegrationType_ = NO_QUADRATURE;
	foundMatch = true;
	}
	}
	else if (strcmp(arg[argIndex],"hoover")==0) {
	regulatorTarget_ = DYNAMICS;
	couplingMode_ = FIXED;
	howOften_ = 1;
	boundaryIntegrationType_ = NO_QUADRATURE;
	foundMatch = true;
	}
	else if (strcmp(arg[argIndex],"flux")==0) {
	regulatorTarget_ = DYNAMICS;
	couplingMode_ = FLUX;
	howOften_ = 1;
	argIndex++;

	boundaryIntegrationType_ = atc_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_);
	foundMatch = true;
	}
	// set parameters for numerical matrix solutions unique to this thermostat
	/*! \page man_control_thermal_correction_max_iterations fix_modify AtC control thermal correction_max_iterations
	\section syntax
	fix_modify AtC control thermal correction_max_iterations <max_iterations>
	- max_iterations (int) = maximum number of iterations that will be used by iterative matrix solvers\n

	\section examples
	<TT> fix_modify AtC control thermal correction_max_iterations 10 </TT> \n
	\section description
	Sets the maximum number of iterations to compute the 2nd order in time correction term for lambda with the fractional step method. The method uses the same tolerance as the controller's matrix solver.
	\section restrictions
	only for use with thermal physics using the fractional step method.
	\section related
	\section default
	correction_max_iterations is 20
	*/
	else if (strcmp(arg[argIndex],"correction_max_iterations")==0) {
	argIndex++;
	lambdaMaxIterations_ = atoi(arg[argIndex]);
	if (lambdaMaxIterations_ < 1) {
	throw ATC_Error("Bad correction maximum iteration count");
	}
	foundMatch = true;
	}
	}

	if (!foundMatch)
	foundMatch = AtomicRegulator::modify(narg,arg);
	if (foundMatch)
	needReset_ = true;
	return foundMatch;
	}

	//--------------------------------------------------------
	// reset_lambda_contribution:
	// resets the thermostat generated power to a
	// prescribed value
	//--------------------------------------------------------
	void Thermostat::reset_lambda_contribution(const DENS_MAT & target)
	{
	DENS_MAN * lambdaPowerFiltered = regulator_data("LambdaPowerFiltered",1);
	*lambdaPowerFiltered = target;
	}

	//--------------------------------------------------------
	// construct_methods:
	// instantiations desired regulator method(s)

	// dependence, but in general there is also a
	// time integrator dependence. In general the
	// precedence order is:
	// time filter -> time integrator -> thermostat
	// In the future this may need to be added if
	// different types of time integrators can be
	// specified.
	//--------------------------------------------------------
	void Thermostat::construct_methods()
	{
	// get data associated with stages 1 & 2 of ATC_Method::initialize
	AtomicRegulator::construct_methods();

	if (atc_->reset_methods()) {
	// eliminate existing methods
	delete_method();

	// update time filter
	TimeIntegrator::TimeIntegrationType myIntegrationType = (atc_->time_integrator(TEMPERATURE))->time_integration_type();
	TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
	if (timeFilterManager->need_reset() ) {
	if (myIntegrationType == TimeIntegrator::GEAR) {
	timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT);
	}
	else if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
	timeFilter_ = timeFilterManager->construct(TimeFilterManager::EXPLICIT_IMPLICIT);
	}
	}

	if (timeFilterManager->filter_dynamics()) {
	switch (regulatorTarget_) {
	case NONE: {
	regulatorMethod_ = new RegulatorMethod(this);
	break;
	}
	case FIELD: { // error check, rescale and filtering not supported together
	throw ATC_Error("Cannot use rescaling thermostat with time filtering");
	break;
	}
	case DYNAMICS: {
	switch (couplingMode_) {
	case FIXED: {
	if (use_lumped_lambda_solve()) {
	throw ATC_Error("Thermostat:construct_methods - lumped lambda solve cannot be used with Hoover thermostats");
	}
	if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
	if (md_flux_nodes(TEMPERATURE)) {
	if (!md_fixed_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
	// there are fluxes but no fixed or coupled nodes
	regulatorMethod_ = new ThermostatIntegratorFluxFiltered(this,lambdaMaxIterations_);
	}
	else {
	// there are both fixed and flux nodes
	regulatorMethod_ = new ThermostatFluxFixedFiltered(this,lambdaMaxIterations_);
	}
	}
	else {
	// there are only fixed nodes
	regulatorMethod_ = new ThermostatIntegratorFixedFiltered(this,lambdaMaxIterations_);
	}
	}
	else {
	regulatorMethod_ = new ThermostatHooverVerletFiltered(this);
	}
	break;
	}
	case FLUX: {
	if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
	if (use_lumped_lambda_solve()) {
	throw ATC_Error("Thermostat:construct_methods - lumped lambda solve has been depricated for fractional step thermostats");
	}
	if (md_fixed_nodes(TEMPERATURE)) {
	if (!md_flux_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
	// there are fixed nodes but no fluxes
	regulatorMethod_ = new ThermostatIntegratorFixedFiltered(this,lambdaMaxIterations_);
	}
	else {
	// there are both fixed and flux nodes
	regulatorMethod_ = new ThermostatFluxFixedFiltered(this,lambdaMaxIterations_);
	}
	}
	else {
	// there are only flux nodes
	regulatorMethod_ = new ThermostatIntegratorFluxFiltered(this,lambdaMaxIterations_);
	}
	}
	else {
	if (use_localized_lambda()) {
	if (!((atc_->prescribed_data_manager())->no_fluxes(TEMPERATURE)) &&
	atc_->boundary_integration_type() != NO_QUADRATURE) {
	throw ATC_Error("Cannot use flux coupling with localized lambda");
	}
	}
	regulatorMethod_ = new ThermostatPowerVerletFiltered(this);
	}
	break;
	}
	default:
	throw ATC_Error("Unknown coupling mode in Thermostat::initialize");
	}
	break;
	}
	default:
	throw ATC_Error("Unknown thermostat type in Thermostat::initialize");
	}
	}
	else {
	switch (regulatorTarget_) {
	case NONE: {
	regulatorMethod_ = new RegulatorMethod(this);
	break;
	}
	case FIELD: {
	if (atc_->temperature_def()==KINETIC)
	regulatorMethod_ = new ThermostatRescale(this);
	else if (atc_->temperature_def()==TOTAL)
	regulatorMethod_ = new ThermostatRescaleMixedKePe(this);
	else
	throw ATC_Error("Unknown temperature definition");
	break;
	}
	case DYNAMICS: {
	switch (couplingMode_) {
	case FIXED: {
	if (use_lumped_lambda_solve()) {
	throw ATC_Error("Thermostat:construct_methods - lumped lambda solve cannot be used with Hoover thermostats");
	}
	if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
	if (md_flux_nodes(TEMPERATURE)) {
	if (!md_fixed_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
	// there are fluxes but no fixed or coupled nodes
	regulatorMethod_ = new ThermostatIntegratorFlux(this,lambdaMaxIterations_);
	}
	else {
	// there are both fixed and flux nodes
	regulatorMethod_ = new ThermostatFluxFixed(this,lambdaMaxIterations_);
	}
	}
	else {
	// there are only fixed nodes
	regulatorMethod_ = new ThermostatIntegratorFixed(this,lambdaMaxIterations_);
	}
	}
	else {
	regulatorMethod_ = new ThermostatHooverVerlet(this);
	}
	break;
	}
	case FLUX: {
	if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
	if (use_lumped_lambda_solve()) {
	throw ATC_Error("Thermostat:construct_methods - lumped lambda solve has been depricated for fractional step thermostats");
	}
	if (md_fixed_nodes(TEMPERATURE)) {
	if (!md_flux_nodes(TEMPERATURE) && (boundaryIntegrationType_ == NO_QUADRATURE)) {
	// there are fixed nodes but no fluxes
	regulatorMethod_ = new ThermostatIntegratorFixed(this,lambdaMaxIterations_);
	}
	else {
	// there are both fixed and flux nodes
	regulatorMethod_ = new ThermostatFluxFixed(this,lambdaMaxIterations_);
	}
	}
	else {
	// there are only flux nodes
	regulatorMethod_ = new ThermostatIntegratorFlux(this,lambdaMaxIterations_);
	}
	}
	else {
	if (use_localized_lambda()) {
	if (!((atc_->prescribed_data_manager())->no_fluxes(TEMPERATURE)) &&
	atc_->boundary_integration_type() != NO_QUADRATURE) {
	throw ATC_Error("Cannot use flux coupling with localized lambda");
	}
	}
	regulatorMethod_ = new ThermostatPowerVerlet(this);
	}
	break;
	}
	default:
	throw ATC_Error("Unknown coupling mode in Thermostat::initialize");
	}
	break;
	}
	default:
	throw ATC_Error("Unknown thermostat target in Thermostat::initialize");
	}
	}

	AtomicRegulator::reset_method();
	}
	else {
	set_all_data_to_used();
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatShapeFunction
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatShapeFunction::ThermostatShapeFunction(AtomicRegulator * thermostat,
	const string & regulatorPrefix) :
	RegulatorShapeFunction(thermostat,regulatorPrefix),
	mdMassMatrix_(atc_->set_mass_mat_md(TEMPERATURE)),
	atomVelocities_(NULL)
	{
	fieldMask_(TEMPERATURE,FLUX) = true;
	lambda_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); // data associated with stage 3 in ATC_Method::initialize
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatShapeFunction::construct_transfers()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());
	RegulatorShapeFunction::construct_transfers();

	// get atom velocity data from manager
	atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);

	// construct lambda evaluated at atom locations
	atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
	lambda_,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
	}

	//---------------------------------------------------------
	// set_weights:
	// set the diagonal weighting matrix to be the atomic
	// temperatures
	//---------------------------------------------------------
	void ThermostatShapeFunction::set_weights()
	{
	if (this->use_local_shape_functions()) {
	VelocitySquaredMapped * myWeights = new VelocitySquaredMapped(atc_,lambdaAtomMap_);
	weights_ = myWeights;
	(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
	regulatorPrefix_+"AtomVelocitySquaredMapped");
	}
	else {
	VelocitySquared * myWeights = new VelocitySquared(atc_);
	weights_ = myWeights;
	(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
	regulatorPrefix_+"AtomVelocitySquared");
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatRescale
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatRescale::ThermostatRescale(AtomicRegulator * thermostat) :
	ThermostatShapeFunction(thermostat),
	nodalTemperature_(atc_->field(TEMPERATURE)),
	atomVelocityRescalings_(NULL)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatRescale::construct_transfers()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// set up node mappings
	create_node_maps();

	// set up data for linear solver
	shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
	(atc_->interscale_manager()).add_per_atom_sparse_matrix(shapeFunctionMatrix_,
	regulatorPrefix_+"LambdaCouplingMatrixEnergy");
	linearSolverType_ = AtomicRegulator::CG_SOLVE;

	// base class transfers
	ThermostatShapeFunction::construct_transfers();

	// velocity rescaling factor
	atomVelocityRescalings_ = new AtomicVelocityRescaleFactor(atc_,atomLambdas_);
	interscaleManager.add_per_atom_quantity(atomVelocityRescalings_,
	regulatorPrefix_+"AtomVelocityRescaling");
	}

	//---------------------------------------------------------
	// set_weights:
	// set the diagonal weighting matrix to be the atomic
	// temperatures
	//---------------------------------------------------------
	void ThermostatRescale::set_weights()
	{
	weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicEnergyForTemperature");
	}

	//--------------------------------------------------------
	// apply_post_corrector:
	// apply the thermostat in the post corrector phase
	//--------------------------------------------------------
	void ThermostatRescale::apply_post_corrector(double dt)
	{
	compute_thermostat(dt);

	// application of rescaling lambda due
	apply_to_atoms(atomVelocities_);
	}

	//--------------------------------------------------------
	// compute_thermostat
	// manages the solution of the
	// thermostat equations and variables
	//--------------------------------------------------------
	void ThermostatRescale::compute_thermostat(double dt)
	{
	// compute right-hand side
	this->set_rhs(_rhs_);

	// solve equations
	solve_for_lambda(_rhs_,lambda_->set_quantity());
	}

	//--------------------------------------------------------
	// set_rhs:
	// constructs the RHS vector with the target
	// temperature
	//--------------------------------------------------------
	void ThermostatRescale::set_rhs(DENS_MAT & rhs)
	{
	rhs = mdMassMatrix_.quantity()*nodalTemperature_.quantity();
	}

	//--------------------------------------------------------
	// apply_lambda_to_atoms:
	// applies the velocity rescale with an existing lambda
	// note oldAtomicQuantity and dt are not used
	//--------------------------------------------------------
	void ThermostatRescale::apply_to_atoms(PerAtomQuantity<double> * atomVelocities)
	{
	atomVelocities = atomVelocityRescalings_->quantity();
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatRescale::output(OUTPUT_LIST & outputData)
	{
	DENS_MAT & lambda(lambda_->set_quantity());
	if ((atc_->lammps_interface())->rank_zero()) {
	outputData["LambdaEnergy"] = λ
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatRescaleMixedKePe
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatRescaleMixedKePe::ThermostatRescaleMixedKePe(AtomicRegulator * thermostat) :
	ThermostatRescale(thermostat),
	nodalAtomicFluctuatingPotentialEnergy_(NULL)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatRescaleMixedKePe::construct_transfers()
	{
	ThermostatRescale::construct_transfers();
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// get fluctuating PE at nodes
	nodalAtomicFluctuatingPotentialEnergy_ =
	interscaleManager.dense_matrix("NodalAtomicFluctuatingPotentialEnergy");
	}

	//---------------------------------------------------------
	// set_weights:
	// set the diagonal weighting matrix to be the atomic
	// temperatures
	//---------------------------------------------------------
	void ThermostatRescaleMixedKePe::set_weights()
	{
	weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicTwiceKineticEnergy");
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatRescaleMixedKePe::initialize()
	{
	ThermostatRescale::initialize();
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// multipliers for KE and PE
	AtomicEnergyForTemperature * atomEnergyForTemperature =
	static_cast<AtomicEnergyForTemperature * >(interscaleManager.per_atom_quantity("AtomicEnergyForTemperature"));
	keMultiplier_ = atomEnergyForTemperature->kinetic_energy_multiplier();
	peMultiplier_ = 2. - keMultiplier_;
	keMultiplier_ /= 2.; // account for use of 2 X KE in matrix equation
	}

	//--------------------------------------------------------
	// set_rhs:
	// accounts for potential energy contribution to
	// definition of atomic temperature
	//--------------------------------------------------------
	void ThermostatRescaleMixedKePe::set_rhs(DENS_MAT & rhs)
	{
	ThermostatRescale::set_rhs(rhs);
	rhs -= peMultiplier_*(nodalAtomicFluctuatingPotentialEnergy_->quantity());
	rhs /= keMultiplier_;
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatFsSolver
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatFsSolver::ThermostatFsSolver(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	RegulatorShapeFunction(thermostat,regulatorPrefix),
	lambdaMaxIterations_(lambdaMaxIterations),
	rhsLambdaSquared_(NULL),
	dtFactor_(1.)
	{
	fieldMask_(TEMPERATURE,FLUX) = true;
	lambda_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1); // data associated with stage 3 in ATC_Method::initialize
	}

	//--------------------------------------------------------
	// initialize
	// creates mapping from all nodes to those to which
	// the thermostat applies
	//--------------------------------------------------------
	void ThermostatFsSolver::initialize()
	{
	RegulatorShapeFunction::initialize();

	rhsMap_.resize(overlapToNodeMap_->nRows(),1);
	DENS_MAT rhsMapGlobal(nNodes_,1);
	const set<int> & applicationNodes(applicationNodes_->quantity());
	for (int i = 0; i < nNodes_; i++) {
	if (applicationNodes.find(i) != applicationNodes.end()) {
	rhsMapGlobal(i,0) = 1.;
	}
	else {
	rhsMapGlobal(i,0) = 0.;
	}
	}
	map_unique_to_overlap(rhsMapGlobal,rhsMap_);
	}

	//---------------------------------------------------------
	// set_weights:
	// set the diagonal weighting matrix to be the atomic
	// temperatures
	//---------------------------------------------------------
	void ThermostatFsSolver::set_weights()
	{
	if (this->use_local_shape_functions()) {
	VelocitySquaredMapped * myWeights = new VelocitySquaredMapped(atc_,lambdaAtomMap_);
	weights_ = myWeights;
	(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
	regulatorPrefix_+"AtomVelocitySquaredMapped");
	}
	else {
	VelocitySquared * myWeights = new VelocitySquared(atc_);
	weights_ = myWeights;
	(atc_->interscale_manager()).add_per_atom_quantity(myWeights,
	regulatorPrefix_+"AtomVelocitySquared");
	}
	}

	//--------------------------------------------------------
	// compute_lambda:
	// solves linear system for lambda, if the
	// bool is true it iterators to a non-linear solution
	//--------------------------------------------------------
	void ThermostatFsSolver::compute_lambda(const DENS_MAT & rhs,
	bool iterateSolution)
	{
	// solve linear system for lambda guess
	DENS_MAT & lambda(lambda_->set_quantity());
	solve_for_lambda(rhs,lambda);

	// iterate to solution
	if (iterateSolution) {
	iterate_lambda(rhs);
	}
	}

	//--------------------------------------------------------
	// iterate_lambda:
	// iteratively solves the equations for lambda
	// for the higher order dt corrections, assuming
	// an initial guess for lambda
	//--------------------------------------------------------
	void ThermostatFsSolver::iterate_lambda(const MATRIX & rhs)
	{
	int nNodeOverlap = overlapToNodeMap_->nRows();
	DENS_VEC _lambdaOverlap_(nNodeOverlap);
	DENS_MAT & lambda(lambda_->set_quantity());
	map_unique_to_overlap(lambda,_lambdaOverlap_);
	double factor = 0.5dtFactor_atc_->dt();

	_lambdaOld_.resize(nNodes_,1);
	_rhsOverlap_.resize(nNodeOverlap,1);
	map_unique_to_overlap(rhs,_rhsOverlap_);
	_rhsTotal_.resize(nNodeOverlap);

	// solve assuming we get initial guess for lambda
	double error(-1.);
	for (int i = 0; i < lambdaMaxIterations_; ++i) {
	_lambdaOld_ = lambda;

	// solve the system with the new rhs
	const DENS_MAT & rhsLambdaSquared(rhsLambdaSquared_->quantity());
	for (int i = 0; i < nNodeOverlap; i++) {
	if (rhsMap_(i,0) == 1.) {
	_rhsTotal_(i) = _rhsOverlap_(i,0) + factor*rhsLambdaSquared(i,0);
	}
	else {
	_rhsTotal_(i) = 0.;
	}
	}
	matrixSolver_->execute(_rhsTotal_,_lambdaOverlap_);

	// check convergence
	map_overlap_to_unique(_lambdaOverlap_,lambda);
	lambda_->force_reset();
	DENS_MAT difference = lambda-_lambdaOld_;
	error = difference.col_norm()/_lambdaOld_.col_norm();
	if (error < tolerance_)
	break;
	}

	if (error >= tolerance_) {
	stringstream message;
	message << "WARNING: Iterative solve for lambda failed to converge after " << lambdaMaxIterations_ << " iterations, final tolerance was " << error << "\n";
	ATC::LammpsInterface::instance()->print_msg(message.str());
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatGlcFs
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatGlcFs::ThermostatGlcFs(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	RegulatorMethod(thermostat,regulatorPrefix),
	lambdaSolver_(NULL),
	mdMassMatrix_(atc_->set_mass_mat_md(TEMPERATURE)),
	atomVelocities_(NULL),
	temperature_(atc_->field(TEMPERATURE)),
	timeFilter_(atomicRegulator_->time_filter()),
	nodalAtomicLambdaPower_(NULL),
	lambdaPowerFiltered_(NULL),
	atomLambdas_(NULL),
	atomThermostatForces_(NULL),
	atomMasses_(NULL),
	isFirstTimestep_(true),
	nodalAtomicEnergy_(NULL),
	atomPredictedVelocities_(NULL),
	nodalAtomicPredictedEnergy_(NULL),
	firstHalfAtomForces_(NULL)
	{
	// construct/obtain data corresponding to stage 3 of ATC_Method::initialize
	nodalAtomicLambdaPower_ = thermostat->regulator_data(regulatorPrefix_+"NodalAtomicLambdaPower",1);
	lambdaPowerFiltered_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaPowerFiltered",1);
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatGlcFs::construct_transfers()
	{
	lambdaSolver_->construct_transfers();
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// get atom velocity data from manager
	atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);

	// construct lambda evaluated at atom locations
	atomLambdas_ = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy");
	if (!atomLambdas_) {
	DENS_MAN * lambda = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1);
	atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
	lambda,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
	}

	// get data from manager
	atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
	nodalAtomicEnergy_ = interscaleManager.dense_matrix("NodalAtomicEnergy");

	// thermostat forces based on lambda and the atomic velocities
	atomThermostatForces_ = new AtomicThermostatForce(atc_,atomLambdas_);
	interscaleManager.add_per_atom_quantity(atomThermostatForces_,
	regulatorPrefix_+"AtomThermostatForce");

	// predicted temperature quantities: atom velocities, atom energies, and restricted atom energies
	atomPredictedVelocities_ = new AtcAtomQuantity<double>(atc_,atc_->nsd());
	// MAKE THINGS WORK WITH ONLY ONE PREDICTED VELOCITY, CHECK IT EXISTS
	interscaleManager.add_per_atom_quantity(atomPredictedVelocities_,
	regulatorPrefix_+"AtomicPredictedVelocities");
	AtomicEnergyForTemperature * atomPredictedEnergyForTemperature = new TwiceKineticEnergy(atc_,
	atomPredictedVelocities_);
	interscaleManager.add_per_atom_quantity(atomPredictedEnergyForTemperature,
	regulatorPrefix_+"AtomicPredictedTwiceKineticEnergy");
	nodalAtomicPredictedEnergy_ = new AtfShapeFunctionRestriction(atc_,
	atomPredictedEnergyForTemperature,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_dense_matrix(nodalAtomicPredictedEnergy_,
	regulatorPrefix_+"NodalAtomicPredictedEnergy");
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatGlcFs::initialize()
	{
	RegulatorMethod::initialize();

	TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
	if (!timeFilterManager->end_equilibrate()) {
	// we should reset lambda and lambdaForce to zero in this case
	// implies an initial condition of 0 for the filtered nodal lambda power
	// initial conditions will always be needed when using time filtering
	// however, the fractional step scheme must assume the instantaneous
	// nodal lambda power is 0 initially because all quantities are in delta form
	lambdaSolver_->initialize(); // ensures initial lambda force is zero
	lambdaSolver_->set_lambda_to_value(0.);
	*nodalAtomicLambdaPower_ = 0.; // energy change due to thermostats
	*lambdaPowerFiltered_ = 0.; // filtered energy change due to thermostats
	}
	else {
	lambdaSolver_->initialize();
	// we can grab lambda power variables using time integrator and atc transfer in cases for equilibration
	}

	// sets up time filter for cases where variables temporally filtered
	if (timeFilterManager->need_reset()) {
	// the form of this integrator implies no time filters that require history data can be used
	timeFilter_->initialize(nodalAtomicLambdaPower_->quantity());
	}

	atomThermostatForces_->quantity(); // initialize
	atomThermostatForces_->fix_quantity();
	firstHalfAtomForces_ = atomThermostatForces_; // initialize
	#ifdef OBSOLETE
	compute_rhs_map();
	#endif
	}

	//--------------------------------------------------------
	// reset_atom_materials:
	// resets the localized atom to material map
	//--------------------------------------------------------
	void ThermostatGlcFs::reset_atom_materials(const Array<int> & elementToMaterialMap,
	const MatrixDependencyManager<DenseMatrix, int> * atomElement)
	{
	lambdaSolver_->reset_atom_materials(elementToMaterialMap,
	atomElement);
	}

	//--------------------------------------------------------
	// apply_to_atoms:
	// determines what if any contributions to the
	// atomic moition is needed for
	// consistency with the thermostat
	// and computes the instantaneous induced power
	//--------------------------------------------------------
	void ThermostatGlcFs::apply_to_atoms(PerAtomQuantity<double> * atomicVelocity,
	const DENS_MAN * nodalAtomicEnergy,
	const DENS_MAT & lambdaForce,
	DENS_MAT & nodalAtomicLambdaPower,
	double dt)
	{
	// compute initial contributions to lambda power
	nodalAtomicLambdaPower = nodalAtomicEnergy->quantity();
	nodalAtomicLambdaPower *= -1.;

	// apply lambda force to atoms
	_velocityDelta_ = lambdaForce;
	_velocityDelta_ /= atomMasses_->quantity();
	_velocityDelta_ *= dt;
	(*atomicVelocity) += _velocityDelta_;

	// finalize lambda power
	nodalAtomicLambdaPower += nodalAtomicEnergy->quantity();
	}

	//--------------------------------------------------------
	// full_prediction:
	// flag to perform a full prediction calcalation
	// for lambda rather than using the old value
	//--------------------------------------------------------
	bool ThermostatGlcFs::full_prediction()
	{
	if (isFirstTimestep_ \|\| ((atc_->atom_to_element_map_type() == EULERIAN)
	&& (atc_->atom_to_element_map_frequency() > 1)
	&& (atc_->step() % atc_->atom_to_element_map_frequency() == 0 ))) {
	return true;
	}
	return false;
	}

	//--------------------------------------------------------
	// apply_predictor:
	// apply the thermostat to the atoms in the first step
	// of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatGlcFs::apply_pre_predictor(double dt)
	{
	DENS_MAT & myLambdaPowerFiltered(lambdaPowerFiltered_->set_quantity());
	DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());

	// update filtered power
	timeFilter_->apply_pre_step1(myLambdaPowerFiltered,myNodalAtomicLambdaPower,dt); // equivalent update to measure power as change in energy due to thermostat

	// apply lambda force to atoms and compute instantaneous lambda power for first half of time step
	this->apply_to_atoms(atomVelocities_,nodalAtomicEnergy_,
	firstHalfAtomForces_->quantity(),
	myNodalAtomicLambdaPower,0.5*dt);

	// update nodal variables for first half of time step
	this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy1_,0.5*dt);

	// start update of filtered lambda power
	myNodalAtomicLambdaPower = 0.; // temporary power for first part of update
	timeFilter_->apply_post_step1(myLambdaPowerFiltered,myNodalAtomicLambdaPower,dt);
	}

	//--------------------------------------------------------
	// apply_pre_corrector:
	// apply the thermostat to the atoms in the first part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatGlcFs::apply_pre_corrector(double dt)
	{
	(*atomPredictedVelocities_) = atomVelocities_->quantity();

	// do full prediction if we just redid the shape functions
	if (full_prediction()) {
	this->compute_lambda(dt);

	atomThermostatForces_->unfix_quantity(); // allow update of atomic force applied by lambda
	}

	// apply lambda force to atoms and compute instantaneous lambda power to predict second half of time step
	DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
	apply_to_atoms(atomPredictedVelocities_,
	nodalAtomicPredictedEnergy_,
	firstHalfAtomForces_->quantity(),
	myNodalAtomicLambdaPower,0.5*dt);

	if (full_prediction())
	atomThermostatForces_->fix_quantity();

	// update predicted nodal variables for second half of time step
	this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy2_,0.5*dt);
	// following manipulations performed this way for efficiency
	deltaEnergy1_ += deltaEnergy2_;
	atc_->apply_inverse_mass_matrix(deltaEnergy1_,TEMPERATURE);
	temperature_ += deltaEnergy1_;
	}

	//--------------------------------------------------------
	// apply_post_corrector:
	// apply the thermostat to the atoms in the second part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatGlcFs::apply_post_corrector(double dt)
	{
	// remove predicted power effects
	DENS_MAT & myTemperature(temperature_.set_quantity());
	atc_->apply_inverse_mass_matrix(deltaEnergy2_,TEMPERATURE);
	myTemperature -= deltaEnergy2_;

	// set up equation and update lambda
	this->compute_lambda(dt);

	// apply lambda force to atoms and compute instantaneous lambda power for second half of time step
	DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
	// allow computation of force applied by lambda using current velocities
	atomThermostatForces_->unfix_quantity();
	atomThermostatForces_->quantity();
	atomThermostatForces_->fix_quantity();
	apply_to_atoms(atomVelocities_,nodalAtomicEnergy_,
	atomThermostatForces_->quantity(),
	myNodalAtomicLambdaPower,0.5*dt);

	// finalize filtered lambda power by adding latest contribution
	timeFilter_->apply_post_step2(lambdaPowerFiltered_->set_quantity(),
	myNodalAtomicLambdaPower,dt);

	// update nodal variables for second half of time step
	this->add_to_energy(myNodalAtomicLambdaPower,deltaEnergy2_,0.5*dt);
	atc_->apply_inverse_mass_matrix(deltaEnergy2_,TEMPERATURE);
	myTemperature += deltaEnergy2_;


	isFirstTimestep_ = false;
	}

	//--------------------------------------------------------
	// compute_lambda:
	// sets up and solves linear system for lambda, if the
	// bool is true it iterators to a non-linear solution
	//--------------------------------------------------------
	void ThermostatGlcFs::compute_lambda(double dt,
	bool iterateSolution)
	{
	// set up rhs for lambda equation
	this->set_thermostat_rhs(rhs_,0.5*dt);

	// solve system
	lambdaSolver_->compute_lambda(rhs_,iterateSolution);
	#ifdef OBSOLETE
	// solve linear system for lambda guess
	DENS_MAT & lambda(lambda_->set_quantity());
	solve_for_lambda(rhs_,lambda);

	// iterate to solution
	if (iterateSolution) {
	iterate_lambda(rhs_);
	}
	#endif
	}

	//--------------------------------------------------------
	// compute_boundary_flux
	// default computation of boundary flux based on
	// finite
	//--------------------------------------------------------
	void ThermostatGlcFs::compute_boundary_flux(FIELDS & fields)
	{

	lambdaSolver_->compute_boundary_flux(fields);
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatGlcFs::output(OUTPUT_LIST & outputData)
	{
	_lambdaPowerOutput_ = nodalAtomicLambdaPower_->quantity();
	// approximate value for lambda power
	double dt = LammpsInterface::instance()->dt();
	_lambdaPowerOutput_ *= (2./dt);
	DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity());
	if ((atc_->lammps_interface())->rank_zero()) {
	outputData[regulatorPrefix_+"LambdaEnergy"] = λ
	outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_);
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatSolverFlux
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatSolverFlux::ThermostatSolverFlux(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	ThermostatFsSolver(thermostat,lambdaMaxIterations,regulatorPrefix)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatSolverFlux::construct_transfers()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// set up node mappings
	create_node_maps();

	// set up data for linear solver
	shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
	interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
	regulatorPrefix_+"LambdaCouplingMatrixEnergy");
	if (elementMask_) {
	lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
	interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
	regulatorPrefix_+"LambdaAtomMap");
	}
	if (atomicRegulator_->use_localized_lambda()) {
	linearSolverType_ = AtomicRegulator::RSL_SOLVE;
	}
	else {
	linearSolverType_ = AtomicRegulator::CG_SOLVE;
	}

	// base class transfers
	ThermostatFsSolver::construct_transfers();

	// add transfers for computation of extra RHS term accounting of O(lambda^2)
	// lambda squared followed by fractional step RHS contribution
	atomLambdas_ = (interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy"));
	if (!atomLambdas_) {
	atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
	lambda_,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
	}
	LambdaSquared * lambdaSquared = new LambdaSquared(atc_,
	interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
	weights_,
	atomLambdas_);
	interscaleManager.add_per_atom_quantity(lambdaSquared,
	regulatorPrefix_+"LambdaSquaredMapped");
	rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_);
	interscaleManager.add_dense_matrix(rhsLambdaSquared_,
	regulatorPrefix_+"RhsLambdaSquared");
	}

	//--------------------------------------------------------
	// construct_regulated_nodes:
	// constructs the set of nodes being regulated
	//--------------------------------------------------------
	void ThermostatSolverFlux::construct_regulated_nodes()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// matrix requires all entries even if localized for correct lumping
	regulatedNodes_ = interscaleManager.set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
	if (!regulatedNodes_) {
	regulatedNodes_ = new RegulatedNodes(atc_);
	interscaleManager.add_set_int(regulatedNodes_,
	regulatorPrefix_+"ThermostatRegulatedNodes");
	}

	// if localized monitor nodes with applied fluxes
	if (atomicRegulator_->use_localized_lambda()) {
	if ((atomicRegulator_->coupling_mode() == Thermostat::FLUX) && (atomicRegulator_->boundary_integration_type() != NO_QUADRATURE)) {
	// include boundary nodes
	applicationNodes_ = new FluxBoundaryNodes(atc_);

	boundaryNodes_ = new BoundaryNodes(atc_);
	interscaleManager.add_set_int(boundaryNodes_,
	regulatorPrefix_+"ThermostatBoundaryNodes");
	}
	else {
	// fluxed nodes only
	applicationNodes_ = new FluxNodes(atc_);
	}
	interscaleManager.add_set_int(applicationNodes_,
	regulatorPrefix_+"ThermostatApplicationNodes");
	}
	else {
	applicationNodes_ = regulatedNodes_;
	}

	// special set of boundary elements for boundary flux quadrature
	if ((atomicRegulator_->boundary_integration_type() == FE_INTERPOLATION)
	&& (atomicRegulator_->use_localized_lambda())) {
	elementMask_ = interscaleManager.dense_matrix_bool(regulatorPrefix_+"BoundaryElementMask");
	if (!elementMask_) {
	elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
	interscaleManager.add_dense_matrix_bool(elementMask_,
	regulatorPrefix_+"BoundaryElementMask");
	}
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatIntegratorFlux
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatIntegratorFlux::ThermostatIntegratorFlux(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	ThermostatGlcFs(thermostat,lambdaMaxIterations,regulatorPrefix),
	heatSource_(atc_->atomic_source(TEMPERATURE))
	{
	lambdaSolver_ = new ThermostatSolverFlux(thermostat,
	lambdaMaxIterations,
	regulatorPrefix);
	}

	//--------------------------------------------------------
	// add_to_temperature
	// add in contributions from lambda power and boundary
	// flux to the FE temperature
	//--------------------------------------------------------
	void ThermostatIntegratorFlux::add_to_energy(const DENS_MAT & nodalLambdaPower,
	DENS_MAT & deltaEnergy,
	double dt)
	{
	deltaEnergy.resize(nNodes_,1);
	const DENS_MAT & myBoundaryFlux(boundaryFlux_[TEMPERATURE].quantity());
	for (int i = 0; i < nNodes_; i++) {
	deltaEnergy(i,0) = nodalLambdaPower(i,0) + dt*myBoundaryFlux(i,0);
	}
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatIntegratorFlux::initialize()
	{
	ThermostatGlcFs::initialize();

	// timestep factor
	lambdaSolver_->set_timestep_factor(1.);
	}

	//--------------------------------------------------------
	// set_thermostat_rhs:
	// sets up the right-hand side including boundary
	// fluxes (coupling & prescribed), heat sources, and
	// fixed (uncoupled) nodes
	//--------------------------------------------------------
	void ThermostatIntegratorFlux::set_thermostat_rhs(DENS_MAT & rhs,
	double dt)
	{

	// only tested with flux != 0 + ess bc = 0

	// (a) for flux based :
	// form rhs : 2/3kB * W_I^-1 * \int N_I r dV
	// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
	// fluxes are set in ATC transfer
	const DENS_MAT & heatSource(heatSource_.quantity());
	#if true
	const set<int> & applicationNodes((lambdaSolver_->application_nodes())->quantity());
	rhs.resize(nNodes_,1);
	for (int i = 0; i < nNodes_; i++) {
	if (applicationNodes.find(i) != applicationNodes.end()) {
	rhs(i,0) = heatSource(i,0);
	}
	else {
	rhs(i,0) = 0.;
	}
	}
	#else
	rhs.resize(nNodes_,1);
	for (int i = 0; i < nNodes_; i++) {
	rhs(i,0) = heatSource(i,0);
	}
	#endif
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatSolverFixed
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatSolverFixed::ThermostatSolverFixed(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	ThermostatFsSolver(thermostat,lambdaMaxIterations,regulatorPrefix)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatSolverFixed::construct_transfers()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// set up node mappings
	create_node_maps();

	// determine if map is needed and set up if so
	if (this->use_local_shape_functions()) {
	lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
	interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
	regulatorPrefix_+"LambdaAtomMap");
	shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
	lambdaAtomMap_,
	nodeToOverlapMap_);
	}
	else {
	shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
	}
	interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
	regulatorPrefix_+"LambdaCouplingMatrixEnergy");
	linearSolverType_ = AtomicRegulator::CG_SOLVE;

	// base class transfers, e.g. weights
	ThermostatFsSolver::construct_transfers();

	// add transfers for computation of extra RHS term accounting of O(lambda^2)
	// lambda squared followed by fractional step RHS contribution
	atomLambdas_ = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaEnergy");
	if (!atomLambdas_) {
	atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
	lambda_,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_per_atom_quantity(atomLambdas_,regulatorPrefix_+"AtomLambdaEnergy");
	}
	if (lambdaAtomMap_) {
	LambdaSquaredMapped * lambdaSquared = new LambdaSquaredMapped(atc_,
	lambdaAtomMap_,
	interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
	weights_,
	atomLambdas_);
	interscaleManager.add_per_atom_quantity(lambdaSquared,
	regulatorPrefix_+"LambdaSquared");
	rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_);
	}
	else {
	LambdaSquared * lambdaSquared = new LambdaSquared(atc_,
	interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS),
	weights_,
	atomLambdas_);
	interscaleManager.add_per_atom_quantity(lambdaSquared,
	regulatorPrefix_+"LambdaSquaredMapped");
	rhsLambdaSquared_ = new AtfShapeFunctionRestriction(atc_,lambdaSquared,shapeFunctionMatrix_);
	}
	interscaleManager.add_dense_matrix(rhsLambdaSquared_,
	regulatorPrefix_+"RhsLambdaSquared");
	}

	//--------------------------------------------------------
	// construct_regulated_nodes:
	// constructs the set of nodes being regulated
	//--------------------------------------------------------
	void ThermostatSolverFixed::construct_regulated_nodes()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());
	regulatedNodes_ = interscaleManager.set_int(regulatorPrefix_+"ThermostatRegulatedNodes");

	if (!regulatedNodes_) {
	if (!atomicRegulator_->use_localized_lambda()) {
	regulatedNodes_ = new RegulatedNodes(atc_);
	}
	else if (atomicRegulator_->coupling_mode() == AtomicRegulator::FLUX) {
	regulatedNodes_ = new FixedNodes(atc_);
	}
	else if (atomicRegulator_->coupling_mode() == AtomicRegulator::FIXED) {
	// include boundary nodes
	regulatedNodes_ = new FixedBoundaryNodes(atc_);
	}
	else {
	throw ATC_Error("ThermostatSolverFixed::construct_regulated_nodes - couldn't determine set of regulated nodes");
	}

	interscaleManager.add_set_int(regulatedNodes_,
	regulatorPrefix_+"ThermostatRegulatedNodes");
	}

	applicationNodes_ = regulatedNodes_;

	// special set of boundary elements for defining regulated atoms
	if (atomicRegulator_->use_localized_lambda()) {
	elementMask_ = interscaleManager.dense_matrix_bool(regulatorPrefix_+"BoundaryElementMask");
	if (!elementMask_) {
	elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
	interscaleManager.add_dense_matrix_bool(elementMask_,
	regulatorPrefix_+"BoundaryElementMask");
	}
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatIntegratorFixed
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatIntegratorFixed::ThermostatIntegratorFixed(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	ThermostatGlcFs(thermostat,lambdaMaxIterations,regulatorPrefix),
	atomThermostatForcesPredVel_(NULL),
	filterCoefficient_(1.)
	{
	lambdaSolver_ = new ThermostatSolverFixed(thermostat,
	lambdaMaxIterations,
	regulatorPrefix);
	}

	//--------------------------------------------------------
	// construct_regulated_nodes:
	// constructs the set of nodes being regulated
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::construct_transfers()
	{
	ThermostatGlcFs::construct_transfers();

	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// predicted forces for halving update
	atomThermostatForcesPredVel_ = new AtomicThermostatForce(atc_,atomLambdas_,atomPredictedVelocities_);
	interscaleManager.add_per_atom_quantity(atomThermostatForcesPredVel_,
	regulatorPrefix_+"AtomThermostatForcePredictedVelocity");
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::initialize()
	{
	ThermostatGlcFs::initialize();
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// set KE multiplier
	AtomicEnergyForTemperature * atomEnergyForTemperature =
	static_cast<AtomicEnergyForTemperature * >(interscaleManager.per_atom_quantity("AtomicEnergyForTemperature"));
	keMultiplier_ = atomEnergyForTemperature->kinetic_energy_multiplier();

	// reset data to zero
	deltaFeEnergy_.reset(nNodes_,1);
	deltaNodalAtomicEnergy_.reset(nNodes_,1);

	// initialize filtered energy
	TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
	if (!timeFilterManager->end_equilibrate()) {
	nodalAtomicEnergyFiltered_ = nodalAtomicEnergy_->quantity();
	}

	// timestep factor
	lambdaSolver_->set_timestep_factor(0.5);
	}

	//--------------------------------------------------------
	// halve_force:
	// flag to halve the lambda force for improved
	// accuracy
	//--------------------------------------------------------
	bool ThermostatIntegratorFixed::halve_force()
	{
	if (isFirstTimestep_ \|\| ((atc_->atom_to_element_map_type() == EULERIAN)
	&& (atc_->atom_to_element_map_frequency() > 1)
	&& (atc_->step() % atc_->atom_to_element_map_frequency() == 1))) {
	return true;
	}
	return false;
	}

	//--------------------------------------------------------
	// initialize_delta_nodal_atomic_energy:
	// initializes storage for the variable tracking
	// the change in the nodal atomic energy
	// that has occured over the past timestep
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::initialize_delta_nodal_atomic_energy(double dt)
	{
	// initialize delta energy
	const DENS_MAT & myNodalAtomicEnergy(nodalAtomicEnergy_->quantity());
	initialNodalAtomicEnergy_ = myNodalAtomicEnergy;
	initialNodalAtomicEnergy_ *= -1.; // initially stored as negative for efficiency
	timeFilter_->apply_pre_step1(nodalAtomicEnergyFiltered_.set_quantity(),
	myNodalAtomicEnergy,dt);
	}

	//--------------------------------------------------------
	// compute_delta_nodal_atomic_energy:
	// computes the change in the nodal atomic energy
	// that has occured over the past timestep
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::compute_delta_nodal_atomic_energy(double dt)
	{
	// set delta energy based on predicted atomic velocities
	const DENS_MAT & myNodalAtomicEnergy(nodalAtomicEnergy_->quantity());
	timeFilter_->apply_post_step1(nodalAtomicEnergyFiltered_.set_quantity(),
	myNodalAtomicEnergy,dt);
	deltaNodalAtomicEnergy_ = initialNodalAtomicEnergy_;
	deltaNodalAtomicEnergy_ += myNodalAtomicEnergy;
	}

	//--------------------------------------------------------
	// compute_lambda:
	// sets up and solves linear system for lambda
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::compute_lambda(double dt,
	bool iterateSolution)
	{
	// compute predicted changes in nodal atomic energy
	compute_delta_nodal_atomic_energy(dt);

	// change in finite element energy
	deltaFeEnergy_ = initialFeEnergy_;
	deltaFeEnergy_ += (mdMassMatrix_.quantity())*(temperature_.quantity());

	ThermostatGlcFs::compute_lambda(dt,iterateSolution);
	}

	//--------------------------------------------------------
	// apply_predictor:
	// apply the thermostat to the atoms in the first step
	// of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::apply_pre_predictor(double dt)
	{
	// initialize values to be track change in finite element energy over the timestep
	initialize_delta_nodal_atomic_energy(dt);
	initialFeEnergy_ = -1.((mdMassMatrix_.quantity())(temperature_.quantity())); // initially stored as negative for efficiency

	ThermostatGlcFs::apply_pre_predictor(dt);
	}

	//--------------------------------------------------------
	// apply_pre_corrector:
	// apply the thermostat to the atoms in the first part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::apply_pre_corrector(double dt)
	{
	// do full prediction if we just redid the shape functions
	if (full_prediction()) {
	firstHalfAtomForces_ = atomThermostatForces_; // reset in case this time step needed special treatment
	_tempNodalAtomicEnergyFiltered_ = nodalAtomicEnergyFiltered_.quantity();
	}

	ThermostatGlcFs::apply_pre_corrector(dt);

	if (full_prediction()) {
	// reset temporary variables
	nodalAtomicEnergyFiltered_ = _tempNodalAtomicEnergyFiltered_;
	}

	if (halve_force()) {
	// save old velocities if we are doing halving calculation of lambda force
	// copy velocities over into temporary storage
	(*atomPredictedVelocities_) = atomVelocities_->quantity();
	}
	}

	//--------------------------------------------------------
	// apply_post_corrector:
	// apply the thermostat to the atoms in the second part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::apply_post_corrector(double dt)
	{

	bool halveForce = halve_force();

	ThermostatGlcFs::apply_post_corrector(dt);

	// update filtered energy with lambda power
	DENS_MAT & myNodalAtomicLambdaPower(nodalAtomicLambdaPower_->set_quantity());
	timeFilter_->apply_post_step2(nodalAtomicEnergyFiltered_.set_quantity(),
	myNodalAtomicLambdaPower,dt);

	if (halveForce) {
	// Halve lambda force due to fixed temperature constraints
	// 1) makes up for poor initial condition
	// 2) accounts for possibly large value of lambda when atomic shape function values change
	// from eulerian mapping after more than 1 timestep
	// avoids unstable oscillations arising from
	// thermostat having to correct for error introduced in lambda changing the
	// shape function matrices
	lambdaSolver_->scale_lambda(0.5);
	firstHalfAtomForces_ = atomThermostatForcesPredVel_;
	atomThermostatForcesPredVel_->unfix_quantity();
	}
	else {
	firstHalfAtomForces_ = atomThermostatForces_;
	}
	}

	//--------------------------------------------------------
	// add_to_temperature
	// add in contributions from lambda power and boundary
	// flux to the FE temperature
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::add_to_energy(const DENS_MAT & nodalLambdaPower,
	DENS_MAT & deltaEnergy,
	double dt)
	{
	deltaEnergy.resize(nNodes_,1);

	SetDependencyManager<int> * myRegulatedNodes =
	(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
	const set<int> & regulatedNodes(myRegulatedNodes->quantity());

	for (int i = 0; i < nNodes_; i++) {
	if (regulatedNodes.find(i) != regulatedNodes.end()) {
	deltaEnergy(i,0) = 0.;
	}
	else {
	deltaEnergy(i,0) = nodalLambdaPower(i,0);
	}
	}
	}

	//--------------------------------------------------------
	// set_thermostat_rhs:
	// sets up the right-hand side including boundary
	// fluxes (coupling & prescribed), heat sources, and
	// fixed (uncoupled) nodes
	//--------------------------------------------------------
	void ThermostatIntegratorFixed::set_thermostat_rhs(DENS_MAT & rhs,
	double dt)
	{
	// for essential bcs (fixed nodes) :
	// form rhs : (delThetaV - delTheta)/dt
	SetDependencyManager<int> * myRegulatedNodes =
	(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
	const set<int> & regulatedNodes(myRegulatedNodes->quantity());
	double factor = (1./dt)/keMultiplier_;
	rhs.resize(nNodes_,1);

	for (int i = 0; i < nNodes_; i++) {
	if (regulatedNodes.find(i) != regulatedNodes.end()) {
	rhs(i,0) = factor*(deltaNodalAtomicEnergy_(i,0) - deltaFeEnergy_(i,0));
	}
	else {
	rhs(i,0) = 0.;
	}
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatIntegratorFluxFiltered
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatIntegratorFluxFiltered::ThermostatIntegratorFluxFiltered(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	ThermostatIntegratorFlux(thermostat,lambdaMaxIterations,regulatorPrefix)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatIntegratorFluxFiltered::initialize()
	{
	ThermostatIntegratorFlux::initialize();

	TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
	if (!timeFilterManager->end_equilibrate()) {
	// always must start as zero because of filtering scheme
	heatSourceOld_.reset(nNodes_,1);
	instantHeatSource_.reset(nNodes_,1);
	timeStepSource_.reset(nNodes_,1);
	}
	}

	//--------------------------------------------------------
	// apply_post_corrector:
	// apply the thermostat to the atoms in the second part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatIntegratorFluxFiltered::apply_post_corrector(double dt)
	{
	// compute lambda
	ThermostatIntegratorFlux::apply_post_corrector(dt);

	// store data needed for filter inversion of heat flux for thermostat rhs
	instantHeatSource_ = rhs_;
	heatSourceOld_ = heatSource_.quantity();
	}

	//--------------------------------------------------------
	// add_to_temperature
	// add in contributions from lambda power and boundary
	// flux to the FE temperature
	//--------------------------------------------------------
	void ThermostatIntegratorFluxFiltered::add_to_energy(const DENS_MAT & nodalLambdaPower,
	DENS_MAT & deltaEnergy,
	double dt)
	{
	deltaEnergy.reset(nNodes_,1);
	double coef = timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
	const DENS_MAT & myBoundaryFlux(boundaryFlux_[TEMPERATURE].quantity());
	for (int i = 0; i < nNodes_; i++) {
	deltaEnergy(i,0) = coefnodalLambdaPower(i,0) + dtmyBoundaryFlux(i,0);
	}
	}

	//--------------------------------------------------------
	// set_thermostat_rhs:
	// sets up the right-hand side including boundary
	// fluxes (coupling & prescribed), heat sources, and
	// fixed (uncoupled) nodes
	//--------------------------------------------------------
	void ThermostatIntegratorFluxFiltered::set_thermostat_rhs(DENS_MAT & rhs,
	double dt)
	{

	// only tested with flux != 0 + ess bc = 0

	// (a) for flux based :
	// form rhs : 2/3kB * W_I^-1 * \int N_I r dV
	// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
	// fluxes are set in ATC transfer

	// invert heatSource_ to get unfiltered source
	// relevant coefficients from time filter

	double coefF1 = timeFilter_->filtered_coefficient_pre_s1(2.*dt);
	double coefF2 = timeFilter_->filtered_coefficient_post_s1(2.*dt);
	double coefU1 = timeFilter_->unfiltered_coefficient_pre_s1(2.*dt);
	double coefU2 = timeFilter_->unfiltered_coefficient_post_s1(2.*dt);

	const DENS_MAT & heatSource(heatSource_.quantity());
	SetDependencyManager<int> * myApplicationNodes =
	(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatApplicationNodes");
	const set<int> & applicationNodes(myApplicationNodes->quantity());
	rhs.resize(nNodes_,1);
	for (int i = 0; i < nNodes_; i++) {
	if (applicationNodes.find(i) != applicationNodes.end()) {
	rhs(i,0) = heatSource(i,0) - coefF1coefF2heatSourceOld_(i,0) - coefU1coefF2instantHeatSource_(i,0);
	rhs(i,0) /= coefU2;
	}
	else {
	rhs(i,0) = 0.;
	}
	}
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatIntegratorFluxFiltered::output(OUTPUT_LIST & outputData)
	{
	_lambdaPowerOutput_ = lambdaPowerFiltered_->quantity();
	// approximate value for lambda power
	double dt = LammpsInterface::instance()->dt();
	_lambdaPowerOutput_ *= (2./dt);
	DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity());
	if ((atc_->lammps_interface())->rank_zero()) {
	outputData[regulatorPrefix_+"Lambda"] = λ
	outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_);
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatIntegratorFixedFiltered
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatIntegratorFixedFiltered::ThermostatIntegratorFixedFiltered(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	const string & regulatorPrefix) :
	ThermostatIntegratorFixed(thermostat,lambdaMaxIterations,regulatorPrefix)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// initialize_delta_nodal_atomic_energy:
	// initializes storage for the variable tracking
	// the change in the nodal atomic energy
	// that has occured over the past timestep
	//--------------------------------------------------------


	void ThermostatIntegratorFixedFiltered::initialize_delta_nodal_atomic_energy(double dt)
	{
	// initialize delta energy
	DENS_MAT & myNodalAtomicEnergyFiltered(nodalAtomicEnergyFiltered_.set_quantity());
	initialNodalAtomicEnergy_ = myNodalAtomicEnergyFiltered;
	initialNodalAtomicEnergy_ *= -1.; // initially stored as negative for efficiency
	timeFilter_->apply_pre_step1(myNodalAtomicEnergyFiltered,
	nodalAtomicEnergy_->quantity(),dt);
	}

	//--------------------------------------------------------
	// compute_delta_nodal_atomic_energy:
	// computes the change in the nodal atomic energy
	// that has occured over the past timestep
	//--------------------------------------------------------
	void ThermostatIntegratorFixedFiltered::compute_delta_nodal_atomic_energy(double dt)
	{
	// set delta energy based on predicted atomic velocities
	DENS_MAT & myNodalAtomicEnergyFiltered(nodalAtomicEnergyFiltered_.set_quantity());
	timeFilter_->apply_post_step1(myNodalAtomicEnergyFiltered,
	nodalAtomicEnergy_->quantity(),dt);
	deltaNodalAtomicEnergy_ = initialNodalAtomicEnergy_;
	deltaNodalAtomicEnergy_ += myNodalAtomicEnergyFiltered;
	}

	//--------------------------------------------------------
	// add_to_temperature
	// add in contributions from lambda power and boundary
	// flux to the FE temperature
	//--------------------------------------------------------
	void ThermostatIntegratorFixedFiltered::add_to_energy(const DENS_MAT & nodalLambdaPower,
	DENS_MAT & deltaEnergy,
	double dt)
	{
	deltaEnergy.resize(nNodes_,1);
	SetDependencyManager<int> * myRegulatedNodes =
	(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
	const set<int> & regulatedNodes(myRegulatedNodes->quantity());
	double coef = timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
	for (int i = 0; i < nNodes_; i++) {
	if (regulatedNodes.find(i) != regulatedNodes.end()) {
	deltaEnergy(i,0) = 0.;
	}
	else {
	deltaEnergy(i,0) = coef*nodalLambdaPower(i,0);
	}
	}
	}
	//--------------------------------------------------------
	// set_thermostat_rhs:
	// sets up the right-hand side for fixed
	// (coupling & prescribed) temperature values
	//--------------------------------------------------------
	void ThermostatIntegratorFixedFiltered::set_thermostat_rhs(DENS_MAT & rhs,
	double dt)
	{
	// (b) for essential bcs (fixed nodes):
	// form rhs : (delThetaV - delTheta)/dt
	SetDependencyManager<int> * myRegulatedNodes =
	(atc_->interscale_manager()).set_int(regulatorPrefix_+"ThermostatRegulatedNodes");
	const set<int> & regulatedNodes(myRegulatedNodes->quantity());
	double factor = (1./dt)/keMultiplier_;
	factor /= timeFilter_->unfiltered_coefficient_post_s1(2.*dt);
	rhs.resize(nNodes_,1);

	for (int i = 0; i < nNodes_; i++) {
	if (regulatedNodes.find(i) != regulatedNodes.end()) {
	rhs(i,0) = factor*(deltaNodalAtomicEnergy_(i,0) - deltaFeEnergy_(i,0));
	}
	else {
	rhs(i,0) = 0.;
	}
	}
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatIntegratorFixedFiltered::output(OUTPUT_LIST & outputData)
	{
	_lambdaPowerOutput_ = lambdaPowerFiltered_->quantity();
	// approximate value for lambda power
	double dt = LammpsInterface::instance()->dt();
	_lambdaPowerOutput_ *= (2./dt);
	DENS_MAT & lambda((atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaEnergy",1))->set_quantity());
	if ((atc_->lammps_interface())->rank_zero()) {
	outputData[regulatorPrefix_+"Lambda"] = λ
	outputData[regulatorPrefix_+"NodalLambdaPower"] = &(_lambdaPowerOutput_);
	}
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatFluxFixed
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatFluxFixed::ThermostatFluxFixed(AtomicRegulator * thermostat,
	int lambdaMaxIterations,
	bool constructThermostats) :
	RegulatorMethod(thermostat),
	thermostatFlux_(NULL),
	thermostatFixed_(NULL),
	thermostatBcs_(NULL)
	{
	if (constructThermostats) {
	thermostatFlux_ = new ThermostatIntegratorFlux(thermostat,lambdaMaxIterations,regulatorPrefix_+"Flux");
	thermostatFixed_ = new ThermostatIntegratorFixed(thermostat,lambdaMaxIterations,regulatorPrefix_+"Fixed");

	// need to choose BC type based on coupling mode
	if (thermostat->coupling_mode() == AtomicRegulator::FLUX) {
	thermostatBcs_ = thermostatFlux_;
	}
	else if (thermostat->coupling_mode() == AtomicRegulator::FIXED) {
	thermostatBcs_ = thermostatFixed_;
	}
	else {
	throw ATC_Error("ThermostatFluxFixed:create_thermostats - invalid thermostat type provided");
	}
	}
	}

	//--------------------------------------------------------
	// Destructor
	//--------------------------------------------------------
	ThermostatFluxFixed::~ThermostatFluxFixed()
	{
	if (thermostatFlux_) delete thermostatFlux_;
	if (thermostatFixed_) delete thermostatFixed_;
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatFluxFixed::construct_transfers()
	{
	thermostatFlux_->construct_transfers();
	thermostatFixed_->construct_transfers();
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatFluxFixed::initialize()
	{
	thermostatFixed_->initialize();
	thermostatFlux_->initialize();
	}

	//--------------------------------------------------------
	// apply_predictor:
	// apply the thermostat to the atoms in the first step
	// of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatFluxFixed::apply_pre_predictor(double dt)
	{
	thermostatFixed_->apply_pre_predictor(dt);
	thermostatFlux_->apply_pre_predictor(dt);
	}

	//--------------------------------------------------------
	// apply_pre_corrector:
	// apply the thermostat to the atoms in the first part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatFluxFixed::apply_pre_corrector(double dt)
	{
	thermostatFlux_->apply_pre_corrector(dt);
	if (thermostatFixed_->full_prediction()) {
	atc_->set_fixed_nodes();
	}
	thermostatFixed_->apply_pre_corrector(dt);
	}

	//--------------------------------------------------------
	// apply_post_corrector:
	// apply the thermostat to the atoms in the second part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatFluxFixed::apply_post_corrector(double dt)
	{
	thermostatFlux_->apply_post_corrector(dt);
	atc_->set_fixed_nodes();
	thermostatFixed_->apply_post_corrector(dt);
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatFluxFixed::output(OUTPUT_LIST & outputData)
	{
	thermostatFlux_->output(outputData);
	thermostatFixed_->output(outputData);
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatFluxFixedFiltered
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatFluxFixedFiltered::ThermostatFluxFixedFiltered(AtomicRegulator * thermostat,
	int lambdaMaxIterations) :
	ThermostatFluxFixed(thermostat,lambdaMaxIterations,false)
	{
	thermostatFlux_ = new ThermostatIntegratorFluxFiltered(thermostat,lambdaMaxIterations,regulatorPrefix_+"Flux");
	thermostatFixed_ = new ThermostatIntegratorFixedFiltered(thermostat,lambdaMaxIterations,regulatorPrefix_+"Fixed");

	// need to choose BC type based on coupling mode
	if (thermostat->coupling_mode() == AtomicRegulator::FLUX) {
	thermostatBcs_ = thermostatFlux_;
	}
	else if (thermostat->coupling_mode() == AtomicRegulator::FIXED) {
	thermostatBcs_ = thermostatFixed_;
	}
	else {
	throw ATC_Error("ThermostatFluxFixed:create_thermostats - invalid thermostat type provided");
	}
	}
	//--------------------------------------------------------
	// Class ThermostatGlc
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	//--------------------------------------------------------
	ThermostatGlc::ThermostatGlc(AtomicRegulator * thermostat) :
	ThermostatShapeFunction(thermostat),
	timeFilter_(atomicRegulator_->time_filter()),
	lambdaPowerFiltered_(NULL),
	atomThermostatForces_(NULL),
	prescribedDataMgr_(atc_->prescribed_data_manager()),
	atomMasses_(NULL)
	{
	// consistent with stage 3 of ATC_Method::initialize
	lambdaPowerFiltered_= atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaPowerFiltered",1);
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatGlc::construct_transfers()
	{
	ThermostatShapeFunction::construct_transfers();
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// get data from manager
	atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);

	// thermostat forces based on lambda and the atomic velocities
	AtomicThermostatForce * atomThermostatForces = new AtomicThermostatForce(atc_);
	interscaleManager.add_per_atom_quantity(atomThermostatForces,
	regulatorPrefix_+"AtomThermostatForce");
	atomThermostatForces_ = atomThermostatForces;
	}

	//--------------------------------------------------------
	// apply_to_atoms:
	// determines what if any contributions to the
	// atomic moition is needed for
	// consistency with the thermostat
	//--------------------------------------------------------
	void ThermostatGlc::apply_to_atoms(PerAtomQuantity<double> * atomVelocities,
	const DENS_MAT & lambdaForce,
	double dt)
	{
	_velocityDelta_ = lambdaForce;
	_velocityDelta_ /= atomMasses_->quantity();
	_velocityDelta_ *= dt;
	(*atomVelocities) += _velocityDelta_;
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatPowerVerlet
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatPowerVerlet::ThermostatPowerVerlet(AtomicRegulator * thermostat) :
	ThermostatGlc(thermostat),
	nodalTemperatureRoc_(atc_->field_roc(TEMPERATURE)),
	heatSource_(atc_->atomic_source(TEMPERATURE)),
	nodalAtomicPower_(NULL),
	nodalAtomicLambdaPower_(NULL)
	{
	// do nothing
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatPowerVerlet::construct_transfers()
	{
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	// set up node mappings
	create_node_maps();

	// determine if mapping is needed and set up if so
	if (atomicRegulator_->use_localized_lambda()) {
	lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
	interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
	regulatorPrefix_+"LambdaAtomMap");
	}

	// set up linear solver
	if (atomicRegulator_->use_lumped_lambda_solve()) {
	shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
	linearSolverType_ = AtomicRegulator::RSL_SOLVE;
	}
	else {
	if (lambdaAtomMap_) {
	shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
	lambdaAtomMap_,
	nodeToOverlapMap_);
	}
	else {
	shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
	}
	linearSolverType_ = AtomicRegulator::CG_SOLVE;
	}
	interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
	regulatorPrefix_+"LambdaCouplingMatrixEnergy");

	// base class transfers, e.g. weights
	ThermostatGlc::construct_transfers();

	// get managed data
	nodalAtomicPower_ = interscaleManager.dense_matrix("NodalAtomicPower");

	// power induced by lambda
	DotTwiceKineticEnergy * atomicLambdaPower =
	new DotTwiceKineticEnergy(atc_,atomThermostatForces_);
	interscaleManager.add_per_atom_quantity(atomicLambdaPower,
	regulatorPrefix_+"AtomicLambdaPower");

	// restriction to nodes of power induced by lambda
	nodalAtomicLambdaPower_ = new AtfShapeFunctionRestriction(atc_,
	atomicLambdaPower,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_dense_matrix(nodalAtomicLambdaPower_,
	regulatorPrefix_+"NodalAtomicLambdaPower");
	}

	//--------------------------------------------------------
	// initialize
	// initializes all method data
	//--------------------------------------------------------
	void ThermostatPowerVerlet::initialize()
	{
	ThermostatGlc::initialize();

	// sets up time filter for cases where variables temporally filtered
	TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
	if (!timeFilterManager->end_equilibrate()) {
	_nodalAtomicLambdaPowerOut_ = 0.;
	*lambdaPowerFiltered_ = 0.;
	timeFilter_->initialize(lambdaPowerFiltered_->quantity());
	}
	}

	//--------------------------------------------------------
	// apply_predictor:
	// apply the thermostat to the atoms in the first step
	// of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatPowerVerlet::apply_pre_predictor(double dt)
	{
	atomThermostatForces_->unfix_quantity();
	compute_thermostat(0.5*dt);

	// apply lambda force to atoms
	const DENS_MAT & thermostatForces(atomThermostatForces_->quantity());
	atomThermostatForces_->fix_quantity();
	apply_to_atoms(atomVelocities_,thermostatForces,0.5*dt);
	}

	//--------------------------------------------------------
	// apply_pre_corrector:
	// apply the thermostat to the atoms in the first part
	// of the corrector step of the Verlet algorithm
	//--------------------------------------------------------
	void ThermostatPowerVerlet::apply_pre_corrector(double dt)
	{
	atomThermostatForces_->unfix_quantity();
	compute_thermostat(0.5*dt);

	// apply lambda force to atoms
	const DENS_MAT & thermostatForces(atomThermostatForces_->quantity());
	atomThermostatForces_->fix_quantity();
	apply_to_atoms(atomVelocities_,thermostatForces,0.5*dt);
	}

	//--------------------------------------------------------
	// add_to_rhs:
	// determines what if any contributions to the
	// finite element equations are needed for
	// consistency with the thermostat
	//--------------------------------------------------------
	void ThermostatPowerVerlet::add_to_rhs(FIELDS & rhs)
	{
	rhs[TEMPERATURE] += nodalAtomicLambdaPower_->quantity() + boundaryFlux_[TEMPERATURE].quantity();
	}

	//--------------------------------------------------------
	// set_thermostat_rhs:
	// sets up the right-hand side including boundary
	// fluxes (coupling & prescribed), heat sources, and
	// fixed (uncoupled) nodes
	//--------------------------------------------------------
	void ThermostatPowerVerlet::set_thermostat_rhs(DENS_MAT & rhs_nodes)
	{
	// (a) for flux based :
	// form rhs : \int N_I r dV
	// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
	// fluxes are set in ATC transfer
	rhs_nodes = heatSource_.quantity();

	// (b) for ess. bcs
	// form rhs : {sum_a (2 * N_Ia * v_ia * f_ia) - (dtheta/dt)_I}

	// replace rhs for prescribed nodes
	const DENS_MAT & myNodalAtomicPower(nodalAtomicPower_->quantity());
	const DIAG_MAT & myMdMassMatrix(mdMassMatrix_.quantity());
	const DENS_MAT & myNodalTemperatureRoc(nodalTemperatureRoc_.quantity());
	for (int i = 0; i < nNodes_; i++) {
	if (prescribedDataMgr_->is_fixed(i,TEMPERATURE,0)) {
	rhs_nodes(i,0) = 0.5(myNodalAtomicPower(i,0) - myMdMassMatrix(i,i)myNodalTemperatureRoc(i,0));
	}
	}
	}

	//--------------------------------------------------------
	// compute_thermostat:
	// sets up and solves the thermostat equations since
	// they are the same at different parts of the time
	// step
	//--------------------------------------------------------
	void ThermostatPowerVerlet::compute_thermostat(double dt)
	{
	// set up rhs
	set_thermostat_rhs(_rhs_);

	// solve linear system for lambda
	DENS_MAT & myLambda(lambda_->set_quantity());
	solve_for_lambda(_rhs_,myLambda);

	nodalAtomicLambdaPower_->unfix_quantity(); // enable computation of force applied by lambda
	timeFilter_->apply_pre_step1(lambdaPowerFiltered_->set_quantity(),
	nodalAtomicLambdaPower_->quantity(),dt);
	nodalAtomicLambdaPower_->fix_quantity();
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatPowerVerlet::output(OUTPUT_LIST & outputData)
	{
	_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity();
	DENS_MAT & lambda(lambda_->set_quantity());
	if ((atc_->lammps_interface())->rank_zero()) {
	outputData["lambda"] = λ
	outputData["nodalLambdaPower"] = &(_nodalAtomicLambdaPowerOut_);
	}
	}

	//--------------------------------------------------------
	// finish:
	// final tasks after a run
	//--------------------------------------------------------
	void ThermostatPowerVerlet::finish()
	{
	_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity();
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatHooverVerlet
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatHooverVerlet::ThermostatHooverVerlet(AtomicRegulator * thermostat) :
	ThermostatPowerVerlet(thermostat),
	lambdaHoover_(NULL),
	nodalAtomicHooverLambdaPower_(NULL)
	{
	// set up data consistent with stage 3 of ATC_Method::initialize
	lambdaHoover_ = atomicRegulator_->regulator_data(regulatorPrefix_+"LambdaHoover",1);

	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatHooverVerlet::construct_transfers()
	{
	ThermostatPowerVerlet::construct_transfers();
	InterscaleManager & interscaleManager(atc_->interscale_manager());


	FtaShapeFunctionProlongation * atomHooverLambdas = new FtaShapeFunctionProlongation(atc_,
	lambdaHoover_,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_per_atom_quantity(atomHooverLambdas,
	regulatorPrefix_+"AtomHooverLambda");
	AtomicThermostatForce * atomHooverThermostatForces = new AtomicThermostatForce(atc_,atomHooverLambdas);
	interscaleManager.add_per_atom_quantity(atomHooverThermostatForces,
	regulatorPrefix_+"AtomHooverThermostatForce");
	SummedAtomicQuantity<double> * atomTotalThermostatForces =
	new SummedAtomicQuantity<double>(atc_,atomThermostatForces_,atomHooverThermostatForces);
	interscaleManager.add_per_atom_quantity(atomTotalThermostatForces,
	regulatorPrefix_+"AtomTotalThermostatForce");
	atomThermostatForces_ = atomTotalThermostatForces;

	// transfers dependent on time integration method
	DotTwiceKineticEnergy * atomicHooverLambdaPower =
	new DotTwiceKineticEnergy(atc_,atomHooverThermostatForces);
	interscaleManager.add_per_atom_quantity(atomicHooverLambdaPower,
	regulatorPrefix_+"AtomicHooverLambdaPower");

	nodalAtomicHooverLambdaPower_ = new AtfShapeFunctionRestriction(atc_,
	atomicHooverLambdaPower,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_dense_matrix(nodalAtomicHooverLambdaPower_,
	regulatorPrefix_+"NodalAtomicHooverLambdaPower");
	}

	//--------------------------------------------------------
	// add_to_rhs:
	// determines what if any contributions to the
	// finite element equations are needed for
	// consistency with the thermostat
	//--------------------------------------------------------
	void ThermostatHooverVerlet::add_to_rhs(FIELDS & rhs)
	{
	rhs[TEMPERATURE] += _nodalAtomicLambdaPowerOut_;
	}

	//--------------------------------------------------------
	// compute_thermostat:
	// sets up and solves the thermostat equations since
	// they are the same at different parts of the time
	// step
	//--------------------------------------------------------
	void ThermostatHooverVerlet::compute_thermostat(double dt)
	{
	// apply prescribed/extrinsic sources and fixed nodes
	ThermostatPowerVerlet::compute_thermostat(0.5*dt);
	_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); // save power from lambda in power-based thermostat

	// set up Hoover rhs
	set_hoover_rhs(_rhs_);

	// solve linear system for lambda
	DENS_MAT & myLambda(lambdaHoover_->set_quantity());
	solve_for_lambda(_rhs_,myLambda);

	// compute force applied by lambda
	// compute nodal atomic power from Hoover coupling
	// only add in contribution to uncoupled nodes
	if (atomicRegulator_->use_localized_lambda())
	add_to_lambda_power(atomThermostatForces_->quantity(),0.5*dt);
	}

	//--------------------------------------------------------
	// set_hoover_rhs:
	// sets up the right-hand side for fixed value,
	// i.e. Hoover coupling
	//--------------------------------------------------------
	void ThermostatHooverVerlet::set_hoover_rhs(DENS_MAT & rhs)
	{
	// form rhs : sum_a ( N_Ia * v_ia * f_ia) - 0.5M_MD(dtheta/dt)_I
	rhs = nodalAtomicPower_->quantity();
	rhs -= mdMassMatrix_.quantity()*nodalTemperatureRoc_.quantity();
	rhs /= 2.;
	}

	//--------------------------------------------------------
	// add_to_nodal_lambda_power:
	// determines the power exerted by the Hoover
	// thermostat at each FE node
	//--------------------------------------------------------
	void ThermostatHooverVerlet::add_to_lambda_power(const DENS_MAT & myLambdaForce,
	double dt)
	{
	_myNodalLambdaPower_ = nodalAtomicHooverLambdaPower_->quantity();
	const INT_ARRAY & nodeToOverlapMap(nodeToOverlapMap_->quantity());
	for (int i = 0; i < nNodes_; ++i) {
	if (nodeToOverlapMap(i,0)==-1)
	_nodalAtomicLambdaPowerOut_(i,0) += _myNodalLambdaPower_(i,0);
	else
	_myNodalLambdaPower_(i,0) = 0.;
	}
	timeFilter_->apply_post_step1(lambdaPowerFiltered_->set_quantity(),_myNodalLambdaPower_,dt);
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatPowerVerletFiltered
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatPowerVerletFiltered::ThermostatPowerVerletFiltered(AtomicRegulator * thermostat) :
	ThermostatPowerVerlet(thermostat),
	nodalTemperature2Roc_(atc_->field_2roc(TEMPERATURE)),
	fieldsRoc_(atc_->fields_roc()),
	filterScale_((atc_->time_filter_manager())->filter_scale())
	{
	heatSourceRoc_.reset(nNodes_,1);
	fluxRoc_[TEMPERATURE].reset(nNodes_,1);
	}

	//--------------------------------------------------------
	// compute_boundary_flux
	// also sets time derivatives of boundary flux and
	// heat sources
	//--------------------------------------------------------
	void ThermostatPowerVerletFiltered::compute_boundary_flux(FIELDS & fields)
	{
	ThermostatPowerVerlet::compute_boundary_flux(fields);

	// compute boundary flux rate of change
	fluxRoc_[TEMPERATURE] = 0.;
	atc_->compute_boundary_flux(fieldMask_,
	fieldsRoc_,
	fluxRoc_,
	atomMaterialGroups_,
	shpFcnDerivs_);



	// compute extrinsic model rate of change
	(atc_->extrinsic_model_manager()).set_sources(fieldsRoc_,fluxRoc_);
	heatSourceRoc_ = fluxRoc_[TEMPERATURE].quantity();
	}

	//--------------------------------------------------------
	// add_to_rhs:
	// determines what if any contributions to the
	// finite element equations are needed for
	// consistency with the thermostat
	//--------------------------------------------------------
	void ThermostatPowerVerletFiltered::add_to_rhs(FIELDS & rhs)
	{
	rhs[TEMPERATURE] += lambdaPowerFiltered_->quantity() + boundaryFlux_[TEMPERATURE].quantity();
	}

	//--------------------------------------------------------
	// set_thermostat_rhs:
	// sets up the right-hand side including boundary
	// fluxes (coupling & prescribed), heat sources, and
	// fixed (uncoupled) nodes
	//--------------------------------------------------------
	void ThermostatPowerVerletFiltered::set_thermostat_rhs(DENS_MAT & rhs_nodes)
	{
	// (a) for flux based :
	// form rhs : \int N_I r dV
	// vs Wagner, CMAME, 2008 eq(24) RHS_I = 2/(3kB) flux_I
	// fluxes are set in ATC transfer
	rhs_nodes = heatSource_.quantity() + filterScale_*heatSourceRoc_.quantity();

	// (b) for ess. bcs
	// form rhs : {sum_a (N_Ia * v_ia * f_ia) - 0.5*(dtheta/dt)_I}
	const DENS_MAT & myNodalAtomicPower(nodalAtomicPower_->quantity());
	const DIAG_MAT & myMdMassMatrix(mdMassMatrix_.quantity());
	const DENS_MAT & myNodalTemperatureRoc(nodalTemperatureRoc_.quantity());
	const DENS_MAT & myNodalTemperature2Roc(nodalTemperature2Roc_.quantity());
	for (int i = 0; i < nNodes_; i++) {
	if (prescribedDataMgr_->is_fixed(i,TEMPERATURE,0)) {
	rhs_nodes(i,0) = 0.5(myNodalAtomicPower(i,0) - myMdMassMatrix(i,i)(myNodalTemperatureRoc(i,0)+ filterScale_*myNodalTemperature2Roc(i,0)));
	}
	}
	}

	//--------------------------------------------------------
	// output:
	// adds all relevant output to outputData
	//--------------------------------------------------------
	void ThermostatPowerVerletFiltered::output(OUTPUT_LIST & outputData)
	{
	outputData["lambda"] = &(lambda_->set_quantity());
	outputData["nodalLambdaPower"] = &(lambdaPowerFiltered_->set_quantity());
	}

	//--------------------------------------------------------
	//--------------------------------------------------------
	// Class ThermostatHooverVerletFiltered
	//--------------------------------------------------------
	//--------------------------------------------------------

	//--------------------------------------------------------
	// Constructor
	// Grab references to ATC and thermostat data
	//--------------------------------------------------------
	ThermostatHooverVerletFiltered::ThermostatHooverVerletFiltered(AtomicRegulator * thermostat) :
	ThermostatPowerVerletFiltered(thermostat),
	lambdaHoover_(NULL),
	nodalAtomicHooverLambdaPower_(NULL)
	{
	// consistent with stage 3 of ATC_Method::initialize
	lambdaHoover_ = atomicRegulator_->regulator_data("LambdaHoover",1);
	}

	//--------------------------------------------------------
	// constructor_transfers
	// instantiates or obtains all dependency managed data
	//--------------------------------------------------------
	void ThermostatHooverVerletFiltered::construct_transfers()
	{
	ThermostatPowerVerletFiltered::construct_transfers();
	InterscaleManager & interscaleManager(atc_->interscale_manager());

	FtaShapeFunctionProlongation * atomHooverLambdas = new FtaShapeFunctionProlongation(atc_,
	lambdaHoover_,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_per_atom_quantity(atomHooverLambdas,
	regulatorPrefix_+"AtomHooverLambda");
	AtomicThermostatForce * atomHooverThermostatForces = new AtomicThermostatForce(atc_,atomHooverLambdas);
	interscaleManager.add_per_atom_quantity(atomHooverThermostatForces,
	regulatorPrefix_+"AtomHooverThermostatForce");
	SummedAtomicQuantity<double> * atomTotalThermostatForces =
	new SummedAtomicQuantity<double>(atc_,atomThermostatForces_,atomHooverThermostatForces);
	interscaleManager.add_per_atom_quantity(atomTotalThermostatForces,
	regulatorPrefix_+"AtomTotalThermostatForce");
	atomThermostatForces_ = atomTotalThermostatForces;

	// transfers dependent on time integration method
	DotTwiceKineticEnergy * atomicHooverLambdaPower =
	new DotTwiceKineticEnergy(atc_,atomHooverThermostatForces);
	interscaleManager.add_per_atom_quantity(atomicHooverLambdaPower,
	regulatorPrefix_+"AtomicHooverLambdaPower");

	nodalAtomicHooverLambdaPower_ = new AtfShapeFunctionRestriction(atc_,
	atomicHooverLambdaPower,
	interscaleManager.per_atom_sparse_matrix("Interpolant"));
	interscaleManager.add_dense_matrix(nodalAtomicHooverLambdaPower_,
	regulatorPrefix_+"NodalAtomicHooverLambdaPower");
	}

	//--------------------------------------------------------
	// compute_thermostat:
	// sets up and solves the thermostat equations since
	// they are the same at different parts of the time
	// step
	//--------------------------------------------------------
	void ThermostatHooverVerletFiltered::compute_thermostat(double dt)
	{
	// apply prescribed/extrinsic sources and fixed nodes
	ThermostatPowerVerletFiltered::compute_thermostat(0.5*dt);
	_nodalAtomicLambdaPowerOut_ = nodalAtomicLambdaPower_->quantity(); // save power from lambda in power-based thermostat

	// set up Hoover rhs
	set_hoover_rhs(_rhs_);

	// solve linear system for lambda
	DENS_MAT & myLambda(lambdaHoover_->set_quantity());
	solve_for_lambda(_rhs_,myLambda);


	// compute force applied by lambda
	// compute nodal atomic power from Hoover coupling
	// only add in contribution to uncoupled nodes
	if (atomicRegulator_->use_localized_lambda())
	add_to_lambda_power(atomThermostatForces_->quantity(),0.5*dt);
	}

	//--------------------------------------------------------
	// set_hoover_rhs:
	// sets up the right-hand side for fixed value,
	// i.e. Hoover coupling
	//--------------------------------------------------------
	void ThermostatHooverVerletFiltered::set_hoover_rhs(DENS_MAT & rhs)
	{
	// form rhs : sum_a (N_Ia * v_ia * f_ia) - 0.5M_MD(dtheta/dt)_I
	rhs = nodalAtomicPower_->quantity();
	rhs -= mdMassMatrix_.quantity()(nodalTemperatureRoc_.quantity() + filterScale_nodalTemperature2Roc_.quantity());
	rhs /= 2.;
	}

	//--------------------------------------------------------
	// add_to_nodal_lambda_power:
	// determines the power exerted by the Hoover
	// thermostat at each FE node
	//--------------------------------------------------------
	void ThermostatHooverVerletFiltered::add_to_lambda_power(const DENS_MAT & myLambdaForce,
	double dt)
	{
	_myNodalLambdaPower_ = nodalAtomicHooverLambdaPower_->quantity();
	const INT_ARRAY nodeToOverlapMap(nodeToOverlapMap_->quantity());
	for (int i = 0; i < nNodes_; ++i) {
	if (nodeToOverlapMap(i,0)==-1)
	_nodalAtomicLambdaPowerOut_(i,0) += _myNodalLambdaPower_(i,0);
	else
	_myNodalLambdaPower_(i,0) = 0.;
	}
	timeFilter_->apply_post_step1(lambdaPowerFiltered_->set_quantity(),_myNodalLambdaPower_,dt);
	}

	//--------------------------------------------------------
	// add_to_rhs:
	// determines what if any contributions to the
	// finite element equations are needed for
	// consistency with the thermostat
	//--------------------------------------------------------
	void ThermostatHooverVerletFiltered::add_to_rhs(FIELDS & rhs)
	{
	rhs[TEMPERATURE] += lambdaPowerFiltered_->quantity();
	}

	};

Thermostat.cppNo OneTemporaryActions

File Metadata

Thermostat.cppView Options

Event Timeline

Thermostat.cpp
No OneTemporary
Actions

Thermostat.cpp
View Options