fix_rx.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Jun 22, 23:20

fix_rx.cpp
View Options

	/* ----------------------------------------------------------------------
	LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
	http://lammps.sandia.gov, Sandia National Laboratories
	Steve Plimpton, sjplimp@sandia.gov

	Copyright (2003) Sandia Corporation. Under the terms of Contract
	DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
	certain rights in this software. This software is distributed under
	the GNU General Public License.

	See the README file in the top-level LAMMPS directory.
	------------------------------------------------------------------------- */

	#include <stdio.h>
	#include <stdlib.h>
	#include <string.h>
	#include <math.h>
	#include "fix_rx.h"
	#include "atom.h"
	#include "error.h"
	#include "group.h"
	#include "modify.h"
	#include "force.h"
	#include "memory.h"
	#include "comm.h"
	#include "update.h"
	#include "domain.h"
	#include "neighbor.h"
	#include "neigh_list.h"
	#include "neigh_request.h"
	#include "math_special.h"
	#include "pair_dpd_fdt_energy.h"

	#include <float.h> // DBL_EPSILON
	#include <vector> // std::vector<>
	#include <algorithm> // std::max
	#include <cmath> // std::fmod

	using namespace LAMMPS_NS;
	using namespace FixConst;
	using namespace MathSpecial;

	enum{NONE,HARMONIC};
	enum{LUCY};

	#define MAXLINE 1024
	#define DELTA 4

	#ifdef DBL_EPSILON
	#define MY_EPSILON (10.0*DBL_EPSILON)
	#else
	#define MY_EPSILON (10.0*2.220446049250313e-16)
	#endif

	#define SparseKinetics_enableIntegralReactions (true)
	#define SparseKinetics_invalidIndex (-1)

	namespace /* anonymous */
	{

	typedef double TimerType;
	TimerType getTimeStamp(void) { return MPI_Wtime(); }
	double getElapsedTime( const TimerType &t0, const TimerType &t1) { return t1-t0; }

	} // end namespace

	/* ---------------------------------------------------------------------- */

	FixRX::FixRX(LAMMPS lmp, int narg, char *arg) :
	Fix(lmp, narg, arg), mol2param(NULL), nreactions(0),
	params(NULL), Arr(NULL), nArr(NULL), Ea(NULL), tempExp(NULL),
	stoich(NULL), stoichReactants(NULL), stoichProducts(NULL), kR(NULL),
	pairDPDE(NULL), dpdThetaLocal(NULL), sumWeights(NULL), sparseKinetics_nu(NULL),
	sparseKinetics_nuk(NULL), sparseKinetics_inu(NULL), sparseKinetics_isIntegralReaction(NULL),
	kineticsFile(NULL), id_fix_species(NULL),
	id_fix_species_old(NULL), fix_species(NULL), fix_species_old(NULL)
	{
	if (narg < 7 \|\| narg > 12) error->all(FLERR,"Illegal fix rx command");
	nevery = 1;

	nreactions = maxparam = 0;
	params = NULL;
	mol2param = NULL;
	pairDPDE = NULL;
	id_fix_species = NULL;
	id_fix_species_old = NULL;

	const int Verbosity = 1;

	// Keep track of the argument list.
	int iarg = 3;

	// Read the kinetic file in arg[3].
	kineticsFile = arg[iarg++];

	// Determine the local temperature averaging method in arg[4].
	wtFlag = 0;
	localTempFlag = NONE;

	{
	char *word = arg[iarg++];
	if (strcmp(word,"none") == 0){
	wtFlag = 0;
	localTempFlag = NONE;
	}
	else if (strcmp(word,"lucy") == 0){
	wtFlag = LUCY;
	localTempFlag = HARMONIC;
	}
	else
	error->all(FLERR,"Illegal fix rx local temperature weighting technique");
	}

	// Select either sparse and dense matrix
	// representations of the stoichiometric matrix.
	useSparseKinetics = true;
	{
	char *word = arg[iarg++];

	if (strcmp(word,"sparse") == 0)
	useSparseKinetics = true;
	else if (strcmp(word,"dense") == 0)
	useSparseKinetics = false;
	else {
	std::string errmsg = "Illegal command " + std::string(word)
	+ " expected \"sparse\" or \"dense\"\n";
	error->all(FLERR, errmsg.c_str());
	}

	if (comm->me == 0 and Verbosity > 1){
	std::string msg = "FixRX: matrix format is ";
	if (useSparseKinetics)
	msg += std::string("sparse");
	else
	msg += std::string("dense");

	error->message(FLERR, msg.c_str());
	}
	}

	// Determine the ODE solver/stepper strategy in arg[6].
	odeIntegrationFlag = ODE_LAMMPS_RK4;

	{
	char *word = arg[iarg++];
	if (strcmp(word,"lammps_rk4") == 0 \|\| strcmp(word,"rk4") == 0)
	odeIntegrationFlag = ODE_LAMMPS_RK4;
	else if (strcmp(word,"lammps_rkf45") == 0 \|\| strcmp(word,"rkf45") == 0)
	odeIntegrationFlag = ODE_LAMMPS_RKF45;
	else {
	std::string errmsg = "Illegal ODE integration type: " + std::string(word);
	error->all(FLERR, errmsg.c_str());
	}
	}

	/// Set the default ODE parameters here. Modify with arg[].
	/// 'minSteps' has a different meaning for RK4 and RKF45.
	/// RK4: This is the # of steps that will be taken with h = dt_dpd / minSteps;
	/// RKF45: This sets as h0 = dt_dpd / minSteps. If minSteps == 0, RKF45 will
	/// estimate h0 internally. h will be adjusted as needed on subsequent steps.
	minSteps = 1;
	maxIters = 100;
	relTol = 1.0e-6;
	absTol = 1.0e-8;

	diagnosticFrequency = 0;
	for (int i = 0; i < numDiagnosticCounters; ++i){
	diagnosticCounter[i] = 0;
	diagnosticCounterPerODE[i] = NULL;
	}

	if (odeIntegrationFlag == ODE_LAMMPS_RK4 && narg==8){
	char *word = arg[iarg++];
	minSteps = atoi( word );

	if (comm->me == 0 and Verbosity > 1){
	char msg[128];
	sprintf(msg, "FixRX: RK4 numSteps= %d", minSteps);
	error->message(FLERR, msg);
	}
	}
	else if (odeIntegrationFlag == ODE_LAMMPS_RK4 && narg>8){
	error->all(FLERR,"Illegal fix rx command. Too many arguments for RK4 solver.");
	}
	else if (odeIntegrationFlag == ODE_LAMMPS_RKF45){
	// Must have four options.
	if (narg < 11)
	error->all(FLERR,"Illegal fix rx command. Too few arguments for RKF45 solver.");

	minSteps = atoi( arg[iarg++] );
	maxIters = atoi( arg[iarg++] );
	relTol = strtod( arg[iarg++], NULL);
	absTol = strtod( arg[iarg++], NULL);

	if (iarg < narg)
	diagnosticFrequency = atoi( arg[iarg++] );

	// maxIters must be at least minSteps.
	maxIters = std::max( minSteps, maxIters );

	if (comm->me == 0 and Verbosity > 1){
	//printf("FixRX: RKF45 minSteps= %d maxIters= %d absTol= %e relTol= %e\n", minSteps, maxIters, absTol, relTol);
	char msg[128];
	sprintf(msg, "FixRX: RKF45 minSteps= %d maxIters= %d relTol= %.1e absTol= %.1e diagnosticFrequency= %d", minSteps, maxIters, relTol, absTol, diagnosticFrequency);
	error->message(FLERR, msg);
	}
	}

	// Initialize/Create the sparse matrix database.
	sparseKinetics_nu = NULL;
	sparseKinetics_nuk = NULL;
	sparseKinetics_inu = NULL;
	sparseKinetics_isIntegralReaction = NULL;
	sparseKinetics_maxReactants = 0;
	sparseKinetics_maxProducts = 0;
	sparseKinetics_maxSpecies = 0;
	}

	/* ---------------------------------------------------------------------- */

	FixRX::~FixRX()
	{
	printf("Inside FixRX::~FixRX copymode= %d\n", copymode);
	if (copymode) return;

	// De-Allocate memory to prevent memory leak
	for (int ii = 0; ii < nreactions; ii++){
	delete [] stoich[ii];
	delete [] stoichReactants[ii];
	delete [] stoichProducts[ii];
	}
	delete [] Arr;
	delete [] nArr;
	delete [] Ea;
	delete [] tempExp;
	delete [] stoich;
	delete [] stoichReactants;
	delete [] stoichProducts;
	delete [] kR;
	delete [] id_fix_species;
	delete [] id_fix_species_old;

	if (useSparseKinetics){
	memory->destroy( sparseKinetics_nu );
	memory->destroy( sparseKinetics_nuk );
	memory->destroy( sparseKinetics_inu );
	memory->destroy( sparseKinetics_isIntegralReaction );
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::post_constructor()
	{
	int maxspecies = 1000;
	int nUniqueSpecies = 0;
	bool match;

	for (int i = 0; i < modify->nfix; i++)
	if (strncmp(modify->fix[i]->style,"property/atom",13) == 0)
	error->all(FLERR,"fix rx cannot be combined with fix property/atom");

	char *tmpspecies = new char[maxspecies];
	for(int jj=0; jj < maxspecies; jj++)
	tmpspecies[jj] = NULL;

	// open file on proc 0

	FILE *fp;
	fp = NULL;
	if (comm->me == 0) {
	fp = force->open_potential(kineticsFile);
	if (fp == NULL) {
	char str[128];
	sprintf(str,"Cannot open rx file %s",kineticsFile);
	error->one(FLERR,str);
	}
	}

	// Assign species names to tmpspecies array and determine the number of unique species

	int n,nwords;
	char line[MAXLINE],*ptr;
	int eof = 0;
	char * word;

	while (1) {
	if (comm->me == 0) {
	ptr = fgets(line,MAXLINE,fp);
	if (ptr == NULL) {
	eof = 1;
	fclose(fp);
	} else n = strlen(line) + 1;
	}
	MPI_Bcast(&eof,1,MPI_INT,0,world);
	if (eof) break;
	MPI_Bcast(&n,1,MPI_INT,0,world);
	MPI_Bcast(line,n,MPI_CHAR,0,world);

	// strip comment, skip line if blank

	if ((ptr = strchr(line,'#'))) *ptr = '\0';
	nwords = atom->count_words(line);
	if (nwords == 0) continue;

	// words = ptrs to all words in line

	nwords = 0;
	word = strtok(line," \t\n\r\f");
	while (word != NULL){
	word = strtok(NULL, " \t\n\r\f");
	match=false;
	for(int jj=0;jj<nUniqueSpecies;jj++){
	if(strcmp(word,tmpspecies[jj])==0){
	match=true;
	break;
	}
	}
	if(!match){
	if(nUniqueSpecies+1>=maxspecies)
	error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
	tmpspecies[nUniqueSpecies] = new char[strlen(word)+1];
	strcpy(tmpspecies[nUniqueSpecies],word);
	nUniqueSpecies++;
	}
	word = strtok(NULL, " \t\n\r\f");
	if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0) break;
	word = strtok(NULL, " \t\n\r\f");
	}
	}
	atom->nspecies_dpd = nUniqueSpecies;
	nspecies = atom->nspecies_dpd;

	// new id = fix-ID + FIX_STORE_ATTRIBUTE
	// new fix group = group for this fix

	id_fix_species = NULL;
	id_fix_species_old = NULL;

	n = strlen(id) + strlen("_SPECIES") + 1;
	id_fix_species = new char[n];
	n = strlen(id) + strlen("_SPECIES_OLD") + 1;
	id_fix_species_old = new char[n];

	strcpy(id_fix_species,id);
	strcat(id_fix_species,"_SPECIES");
	strcpy(id_fix_species_old,id);
	strcat(id_fix_species_old,"_SPECIES_OLD");

	char *newarg = new char[nspecies+5];
	char *newarg2 = new char[nspecies+5];
	newarg[0] = id_fix_species;
	newarg[1] = group->names[igroup];
	newarg[2] = (char *) "property/atom";
	newarg2[0] = id_fix_species_old;
	newarg2[1] = group->names[igroup];
	newarg2[2] = (char *) "property/atom";
	for(int ii=0; ii<nspecies; ii++){
	char str1[2+strlen(tmpspecies[ii])+1];
	char str2[2+strlen(tmpspecies[ii])+4];
	strcpy(str1,"d_");
	strcpy(str2,"d_");
	strncat(str1,tmpspecies[ii],strlen(tmpspecies[ii]));
	strncat(str2,tmpspecies[ii],strlen(tmpspecies[ii]));
	strncat(str2,"Old",3);
	newarg[ii+3] = new char[strlen(str1)+1];
	newarg2[ii+3] = new char[strlen(str2)+1];
	strcpy(newarg[ii+3],str1);
	strcpy(newarg2[ii+3],str2);
	}
	newarg[nspecies+3] = (char *) "ghost";
	newarg[nspecies+4] = (char *) "yes";
	newarg2[nspecies+3] = (char *) "ghost";
	newarg2[nspecies+4] = (char *) "yes";

	modify->add_fix(nspecies+5,newarg,1);
	fix_species = (FixPropertyAtom *) modify->fix[modify->nfix-1];
	restartFlag = modify->fix[modify->nfix-1]->restart_reset;

	modify->add_fix(nspecies+5,newarg2,1);
	fix_species_old = (FixPropertyAtom *) modify->fix[modify->nfix-1];

	if(nspecies==0) error->all(FLERR,"There are no rx species specified.");

	for(int jj=0;jj<nspecies;jj++) {
	delete[] tmpspecies[jj];
	delete[] newarg[jj+3];
	delete[] newarg2[jj+3];
	}

	delete[] newarg;
	delete[] newarg2;
	delete[] tmpspecies;

	read_file( kineticsFile );

	if (useSparseKinetics)
	this->initSparse();

	// set comm size needed by this Pair
	comm_forward = nspecies*2;
	comm_reverse = 2;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::initSparse()
	{
	const int Verbosity = 1;

	if (comm->me == 0 and Verbosity > 1){
	for (int k = 0; k < nspecies; ++k)
	printf("atom->dname[%d]= %s\n", k, atom->dname[k]);

	printf("stoich[][]\n");
	for (int i = 0; i < nreactions; ++i){
	int nreac_i = 0, nprod_i = 0;
	printf("%d: ", i);
	for (int k = 0; k < nspecies; ++k){
	printf(" %g", stoich[i][k]);
	if (stoich[i][k] < 0.0) nreac_i++;
	else if (stoich[i][k] > 0.0) nprod_i++;
	}
	printf(" : %d %d\n", nreac_i, nprod_i);
	}

	printf("stoichReactants[][]\n");
	for (int i = 0; i < nreactions; ++i){
	int nreac_i = 0;
	printf("%d: ", i);
	for (int k = 0; k < nspecies; ++k){
	printf(" %g", stoichReactants[i][k]);
	if (stoichReactants[i][k] > 0.0) nreac_i++;
	}
	printf(" : %d\n", nreac_i);
	}

	printf("stoichProducts[][]\n");
	for (int i = 0; i < nreactions; ++i){
	int nprod_i = 0;
	printf("%d: ", i);
	for (int k = 0; k < nspecies; ++k){
	printf(" %g", stoichProducts[i][k]);
	if (stoichProducts[i][k] > 0.0) nprod_i++;
	}
	printf(" : %d\n", nprod_i);
	}
	} // if (Verbose)

	// 1) Measure the sparsity of stoich[][]
	int nzeros = 0;
	int mxprod = 0;
	int mxreac = 0;
	int mxspec = 0;
	int nIntegral = 0;
	for (int i = 0; i < nreactions; ++i){
	int nreac_i = 0, nprod_i = 0;
	std::string pstr, rstr;
	bool allAreIntegral = true;
	for (int k = 0; k < nspecies; ++k){
	if (stoichReactants[i][k] == 0 and stoichProducts[i][k] == 0)
	nzeros++;

	if (stoichReactants[i][k] > 0.0){
	allAreIntegral &= (std::fmod( stoichReactants[i][k], 1.0 ) == 0.0);

	nreac_i++;
	if (rstr.length() > 0)
	rstr += " + ";

	char digit[6];
	sprintf(digit, "%4.1f ", stoichReactants[i][k]); rstr += digit;
	rstr += atom->dname[k];
	}
	if (stoichProducts[i][k] > 0.0){
	allAreIntegral &= (std::fmod( stoichProducts[i][k], 1.0 ) == 0.0);

	nprod_i++;
	if (pstr.length() > 0)
	pstr += " + ";

	char digit[6];
	sprintf(digit, "%4.1f ", stoichProducts[i][k]); pstr += digit;

	pstr += atom->dname[k];
	}
	}
	if (comm->me == 0 and Verbosity > 1)
	printf("rx%3d: %d %d %d ... %s %s %s\n", i, nreac_i, nprod_i, allAreIntegral, rstr.c_str(), /reversible[i]/ (false) ? "<=>" : "=", pstr.c_str());

	mxreac = std::max( mxreac, nreac_i );
	mxprod = std::max( mxprod, nprod_i );
	mxspec = std::max( mxspec, nreac_i + nprod_i );
	if (allAreIntegral) nIntegral++;
	}

	if (comm->me == 0 and Verbosity > 1){
	char msg[256];
	sprintf(msg, "FixRX: Sparsity of Stoichiometric Matrix= %.1f%% non-zeros= %d nspecies= %d nreactions= %d maxReactants= %d maxProducts= %d maxSpecies= %d integralReactions= %d", 100(double(nzeros) / (nspecies nreactions)), nzeros, nspecies, nreactions, mxreac, mxprod, (mxreac + mxprod), SparseKinetics_enableIntegralReactions);
	error->message(FLERR, msg);
	}

	// Allocate the sparse matrix data.
	{
	sparseKinetics_maxSpecies = (mxreac + mxprod);
	sparseKinetics_maxReactants = mxreac;
	sparseKinetics_maxProducts = mxprod;

	memory->create( sparseKinetics_nu , nreactions, sparseKinetics_maxSpecies, "sparseKinetics_nu");
	memory->create( sparseKinetics_nuk, nreactions, sparseKinetics_maxSpecies, "sparseKinetics_nuk");

	for (int i = 0; i < nreactions; ++i)
	for (int k = 0; k < sparseKinetics_maxSpecies; ++k){
	sparseKinetics_nu [i][k] = 0.0;
	sparseKinetics_nuk[i][k] = SparseKinetics_invalidIndex; // Initialize with an invalid index.
	}

	if (SparseKinetics_enableIntegralReactions){
	memory->create( sparseKinetics_inu, nreactions, sparseKinetics_maxSpecies, "sparseKinetics_inu");
	memory->create( sparseKinetics_isIntegralReaction, nreactions, "sparseKinetics_isIntegralReaction");

	for (int i = 0; i < nreactions; ++i){
	sparseKinetics_isIntegralReaction[i] = false;
	for (int k = 0; k < sparseKinetics_maxSpecies; ++k)
	sparseKinetics_inu[i][k] = 0;
	}
	}
	}

	// Measure the distribution of the # of moles for the ::fastpowi function.
	std::vector<int> nu_bin(10);

	for (int i = 0; i < nreactions; ++i){
	int nreac_i = 0, nprod_i = 0;
	bool isIntegral_i = true;
	for (int k = 0; k < nspecies; ++k){
	if (stoichReactants[i][k] > 0.0){
	const int idx = nreac_i;
	sparseKinetics_nu [i][idx] = stoichReactants[i][k];
	sparseKinetics_nuk[i][idx] = k;

	isIntegral_i &= (std::fmod( stoichReactants[i][k], 1.0 ) == 0.0);
	if (SparseKinetics_enableIntegralReactions){
	sparseKinetics_inu[i][idx] = (int)sparseKinetics_nu[i][idx];
	if (isIntegral_i){
	if (sparseKinetics_inu[i][idx] >= nu_bin.size())
	nu_bin.resize( sparseKinetics_inu[i][idx] );

	nu_bin[ sparseKinetics_inu[i][idx] ] ++;
	}
	}

	nreac_i++;
	}
	if (stoichProducts[i][k] > 0.0){
	const int idx = sparseKinetics_maxReactants + nprod_i;
	sparseKinetics_nu [i][idx] = stoichProducts[i][k];
	sparseKinetics_nuk[i][idx] = k;

	isIntegral_i &= (std::fmod( sparseKinetics_nu[i][idx], 1.0 ) == 0.0);
	if (SparseKinetics_enableIntegralReactions){
	sparseKinetics_inu[i][idx] = (int) sparseKinetics_nu[i][idx];
	if (isIntegral_i){
	if (sparseKinetics_inu[i][idx] >= nu_bin.size())
	nu_bin.resize( sparseKinetics_inu[i][idx] );

	nu_bin[ sparseKinetics_inu[i][idx] ] ++;
	}
	}

	nprod_i++;
	}
	}

	if (SparseKinetics_enableIntegralReactions)
	sparseKinetics_isIntegralReaction[i] = isIntegral_i;
	}

	if (comm->me == 0 and Verbosity > 1){
	for (int i = 1; i < nu_bin.size(); ++i)
	if (nu_bin[i] > 0)
	printf("nu_bin[%d] = %d\n", i, nu_bin[i]);

	for (int i = 0; i < nreactions; ++i){
	std::string pstr, rstr;

	for (int kk = 0; kk < sparseKinetics_maxReactants; kk++){
	const int k = sparseKinetics_nuk[i][kk];
	if (k != SparseKinetics_invalidIndex){
	if (rstr.length() > 0)
	rstr += " + ";

	char digit[6];
	if (SparseKinetics_enableIntegralReactions and sparseKinetics_isIntegralReaction[i])
	sprintf(digit,"%d ", sparseKinetics_inu[i][kk]);
	else
	sprintf(digit,"%4.1f ", sparseKinetics_nu[i][kk]);
	rstr += digit;
	rstr += atom->dname[k];
	}
	}

	for (int kk = sparseKinetics_maxReactants; kk < sparseKinetics_maxSpecies; kk++){
	const int k = sparseKinetics_nuk[i][kk];
	if (k != SparseKinetics_invalidIndex){
	if (pstr.length() > 0)
	pstr += " + ";

	char digit[6];
	if (SparseKinetics_enableIntegralReactions and sparseKinetics_isIntegralReaction[i])
	sprintf(digit,"%d ", sparseKinetics_inu[i][kk]);
	else
	sprintf(digit,"%4.1f ", sparseKinetics_nu[i][kk]);
	pstr += digit;
	pstr += atom->dname[k];
	}
	}
	if (comm->me == 0 and Verbosity > 1)
	printf("rx%3d: %s %s %s\n", i, rstr.c_str(), /reversible[i]/ (false) ? "<=>" : "=", pstr.c_str());
	}
	// end for nreactions
	}
	// end if Verbose
	}

	/* ---------------------------------------------------------------------- */

	int FixRX::setmask()
	{
	int mask = 0;
	mask \|= PRE_FORCE;
	return mask;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::init()
	{
	pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy",1);
	if (pairDPDE == NULL)
	pairDPDE = (PairDPDfdtEnergy *) force->pair_match("dpd/fdt/energy/kk",1);

	if (pairDPDE == NULL)
	error->all(FLERR,"Must use pair_style dpd/fdt/energy with fix rx");

	bool eos_flag = false;
	for (int i = 0; i < modify->nfix; i++)
	if (strcmp(modify->fix[i]->style,"eos/table/rx") == 0) eos_flag = true;
	if(!eos_flag) error->all(FLERR,"fix rx requires fix eos/table/rx to be specified");

	// need a half neighbor list
	// built whenever re-neighboring occurs

	int irequest = neighbor->request(this,instance_me);
	neighbor->requests[irequest]->pair = 0;
	neighbor->requests[irequest]->fix = 1;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::init_list(int, class NeighList* ptr)
	{
	this->list = ptr;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::setup_pre_force(int vflag)
	{
	int nlocal = atom->nlocal;
	int nghost = atom->nghost;
	int *mask = atom->mask;
	int newton_pair = force->newton_pair;
	double tmp;
	int ii;

	if(restartFlag){
	restartFlag = 0;
	}
	else
	{
	int ode_counter[4] = {0};

	UserRHSData userData;
	userData.kFor = new double[nreactions];
	userData.rxnRateLaw = new double[nreactions];

	double rwork = new double[8nspecies];

	if(localTempFlag){
	int count = nlocal + (newton_pair ? nghost : 0);
	dpdThetaLocal = new double[count];
	memset(dpdThetaLocal, 0, sizeof(double)*count);
	computeLocalTemperature();
	}

	for (int id = 0; id < nlocal; id++)
	for (int ispecies=0; ispecies<nspecies; ispecies++){
	tmp = atom->dvector[ispecies][id];
	atom->dvector[ispecies+nspecies][id] = tmp;
	}

	for (int i = 0; i < nlocal; i++)
	if (mask[i] & groupbit){

	// Set the reaction rate constants to zero: no reactions occur at step 0
	for(int irxn=0;irxn<nreactions;irxn++)
	userData.kFor[irxn] = 0.0;

	if (odeIntegrationFlag == ODE_LAMMPS_RK4)
	rk4(i, rwork, &userData);
	else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
	rkf45(i, rwork, &userData, ode_counter);
	}

	// Communicate the updated momenta and velocities to all nodes
	comm->forward_comm_fix(this);
	if(localTempFlag) delete [] dpdThetaLocal;

	delete [] userData.kFor;
	delete [] userData.rxnRateLaw;
	delete [] rwork;
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::pre_force(int vflag)
	{
	TimerType timer_start = getTimeStamp();

	int nlocal = atom->nlocal;
	int nghost = atom->nghost;
	int *mask = atom->mask;
	double *dpdTheta = atom->dpdTheta;
	int newton_pair = force->newton_pair;

	if(localTempFlag){
	int count = nlocal + (newton_pair ? nghost : 0);
	dpdThetaLocal = new double[count];
	memset(dpdThetaLocal, 0, sizeof(double)*count);
	computeLocalTemperature();
	}

	TimerType timer_localTemperature = getTimeStamp();

	// Zero the counters for the ODE solvers.
	int nSteps = 0;
	int nIters = 0;
	int nFuncs = 0;
	int nFails = 0;

	if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency == 1)
	{
	memory->create( diagnosticCounterPerODE[StepSum], nlocal, "FixRX::diagnosticCounterPerODE");
	memory->create( diagnosticCounterPerODE[FuncSum], nlocal, "FixRX::diagnosticCounterPerODE");
	}

	#pragma omp parallel \
	reduction(+: nSteps, nIters, nFuncs, nFails )
	{
	double rwork = new double[8nspecies];

	UserRHSData userData;
	userData.kFor = new double[nreactions];
	userData.rxnRateLaw = new double[nreactions];

	int ode_counter[4] = { 0 };

	#pragma omp for schedule(runtime)
	for (int i = 0; i < nlocal; i++)
	{
	if (mask[i] & groupbit)
	{
	double theta;
	if (localTempFlag)
	theta = dpdThetaLocal[i];
	else
	theta = dpdTheta[i];

	//Compute the reaction rate constants
	for (int irxn = 0; irxn < nreactions; irxn++)
	userData.kFor[irxn] = Arr[irxn]pow(theta,nArr[irxn])exp(-Ea[irxn]/force->boltz/theta);

	if (odeIntegrationFlag == ODE_LAMMPS_RK4)
	rk4(i, rwork, &userData);
	else if (odeIntegrationFlag == ODE_LAMMPS_RKF45)
	rkf45(i, rwork, &userData, ode_counter);
	}
	}

	nSteps += ode_counter[0];
	nIters += ode_counter[1];
	nFuncs += ode_counter[2];
	nFails += ode_counter[3];

	delete [] rwork;
	delete [] userData.kFor;
	delete [] userData.rxnRateLaw;

	} // end parallel region

	TimerType timer_ODE = getTimeStamp();

	// Communicate the updated momenta and velocities to all nodes
	comm->forward_comm_fix(this);
	if(localTempFlag) delete [] dpdThetaLocal;

	TimerType timer_stop = getTimeStamp();

	double time_ODE = getElapsedTime(timer_localTemperature, timer_ODE);

	printf("me= %d total= %g temp= %g ode= %g comm= %g nlocal= %d nfc= %d %d\n", comm->me,
	getElapsedTime(timer_start, timer_stop),
	getElapsedTime(timer_start, timer_localTemperature),
	getElapsedTime(timer_localTemperature, timer_ODE),
	getElapsedTime(timer_ODE, timer_stop), nlocal, nFuncs, nSteps);

	// Warn the user if a failure was detected in the ODE solver.
	if (nFails > 0){
	char sbuf[128];
	sprintf(sbuf,"in FixRX::pre_force, ODE solver failed for %d atoms.", nFails);
	error->warning(FLERR, sbuf);
	}

	// Compute and report ODE diagnostics, if requested.
	if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency != 0){
	// Update the counters.
	diagnosticCounter[StepSum] += nSteps;
	diagnosticCounter[FuncSum] += nFuncs;
	diagnosticCounter[TimeSum] += time_ODE;
	diagnosticCounter[AtomSum] += nlocal;
	diagnosticCounter[numDiagnosticCounters-1] ++;

	if ( (diagnosticFrequency > 0 &&
	((update->ntimestep - update->firststep) % diagnosticFrequency) == 0) \|\|
	(diagnosticFrequency < 0 && update->ntimestep == update->laststep) )
	this->odeDiagnostics();

	for (int i = 0; i < numDiagnosticCounters; ++i)
	if (diagnosticCounterPerODE[i])
	memory->destroy( diagnosticCounterPerODE[i] );
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::read_file(char *file)
	{
	nreactions = 0;

	// open file on proc 0

	FILE *fp;
	fp = NULL;
	if (comm->me == 0) {
	fp = force->open_potential(file);
	if (fp == NULL) {
	char str[128];
	sprintf(str,"Cannot open rx file %s",file);
	error->one(FLERR,str);
	}
	}

	// Count the number of reactions from kinetics file

	int n,nwords,ispecies;
	char line[MAXLINE],*ptr;
	int eof = 0;

	while (1) {
	if (comm->me == 0) {
	ptr = fgets(line,MAXLINE,fp);
	if (ptr == NULL) {
	eof = 1;
	fclose(fp);
	} else n = strlen(line) + 1;
	}
	MPI_Bcast(&eof,1,MPI_INT,0,world);
	if (eof) break;
	MPI_Bcast(&n,1,MPI_INT,0,world);
	MPI_Bcast(line,n,MPI_CHAR,0,world);

	// strip comment, skip line if blank

	if ((ptr = strchr(line,'#'))) *ptr = '\0';
	nwords = atom->count_words(line);
	if (nwords == 0) continue;

	nreactions++;
	}

	// open file on proc 0
	if (comm->me == 0) fp = force->open_potential(file);

	// read each reaction from kinetics file
	eof=0;
	char * word;
	double tmpStoich;
	double sign;

	Arr = new double[nreactions];
	nArr = new double[nreactions];
	Ea = new double[nreactions];
	tempExp = new double[nreactions];
	stoich = new double*[nreactions];
	stoichReactants = new double*[nreactions];
	stoichProducts = new double*[nreactions];
	for (int ii=0;ii<nreactions;ii++){
	stoich[ii] = new double[nspecies];
	stoichReactants[ii] = new double[nspecies];
	stoichProducts[ii] = new double[nspecies];
	}
	kR = new double[nreactions];
	for (int ii=0;ii<nreactions;ii++){
	for (int jj=0;jj<nspecies;jj++){
	stoich[ii][jj] = 0.0;
	stoichReactants[ii][jj] = 0.0;
	stoichProducts[ii][jj] = 0.0;
	}
	}

	nreactions=0;
	sign = -1.0;
	while (1) {
	if (comm->me == 0) {
	ptr = fgets(line,MAXLINE,fp);
	if (ptr == NULL) {
	eof = 1;
	fclose(fp);
	} else n = strlen(line) + 1;
	}
	MPI_Bcast(&eof,1,MPI_INT,0,world);
	if (eof) break;
	MPI_Bcast(&n,1,MPI_INT,0,world);
	MPI_Bcast(line,n,MPI_CHAR,0,world);

	// strip comment, skip line if blank

	if ((ptr = strchr(line,'#'))) *ptr = '\0';
	nwords = atom->count_words(line);
	if (nwords == 0) continue;

	// words = ptrs to all words in line

	nwords = 0;
	word = strtok(line," \t\n\r\f");
	while (word != NULL){
	tmpStoich = atof(word);
	word = strtok(NULL, " \t\n\r\f");
	for (ispecies = 0; ispecies < nspecies; ispecies++){
	if (strcmp(word,&atom->dname[ispecies][0]) == 0){
	stoich[nreactions][ispecies] += sign*tmpStoich;
	if(sign<0.0)
	stoichReactants[nreactions][ispecies] += tmpStoich;
	else stoichProducts[nreactions][ispecies] += tmpStoich;
	break;
	}
	}
	if(ispecies==nspecies){
	if (comm->me) {
	fprintf(stderr,"%s mol fraction is not found in data file\n",word);
	fprintf(stderr,"nspecies=%d ispecies=%d\n",nspecies,ispecies);
	}
	error->all(FLERR,"Illegal fix rx command");
	}
	word = strtok(NULL, " \t\n\r\f");
	if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
	if(strcmp(word,"=") == 0) sign = 1.0;
	if(strcmp(word,"+") != 0 && strcmp(word,"=") != 0){
	if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
	Arr[nreactions] = atof(word);
	word = strtok(NULL, " \t\n\r\f");
	if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
	nArr[nreactions] = atof(word);
	word = strtok(NULL, " \t\n\r\f");
	if(word==NULL) error->all(FLERR,"Missing parameters in reaction kinetic equation");
	Ea[nreactions] = atof(word);
	sign = -1.0;
	break;
	}
	word = strtok(NULL, " \t\n\r\f");
	}
	nreactions++;
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::setupParams()
	{
	int i,j,n;

	// set mol2param for all combinations
	// must be a single exact match to lines read from file

	memory->destroy(mol2param);
	memory->create(mol2param,nspecies,"pair:mol2param");

	for (i = 0; i < nspecies; i++) {
	n = -1;
	for (j = 0; j < nreactions; j++) {
	if (i == params[j].ispecies) {
	if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
	n = j;
	}
	}
	mol2param[i] = n;
	}
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::rk4(int id, double rwork, void v_params)
	{
	double *k1 = rwork;
	double *k2 = k1 + nspecies;
	double *k3 = k2 + nspecies;
	double *k4 = k3 + nspecies;
	double *y = k4 + nspecies;
	double *yp = y + nspecies;

	const int numSteps = minSteps;

	const double h = update->dt / double(numSteps);

	// Update ConcOld
	for (int ispecies = 0; ispecies < nspecies; ispecies++)
	{
	const double tmp = atom->dvector[ispecies][id];
	atom->dvector[ispecies+nspecies][id] = tmp;
	y[ispecies] = tmp;
	}

	// Run the requested steps with h.
	for (int step = 0; step < numSteps; step++)
	{
	// k1
	rhs(0.0,y,k1,v_params);

	// k2
	for (int ispecies = 0; ispecies < nspecies; ispecies++)
	yp[ispecies] = y[ispecies] + 0.5hk1[ispecies];

	rhs(0.0,yp,k2,v_params);

	// k3
	for (int ispecies = 0; ispecies < nspecies; ispecies++)
	yp[ispecies] = y[ispecies] + 0.5hk2[ispecies];

	rhs(0.0,yp,k3,v_params);

	// k4
	for (int ispecies = 0; ispecies < nspecies; ispecies++)
	yp[ispecies] = y[ispecies] + h*k3[ispecies];

	rhs(0.0,yp,k4,v_params);

	for (int ispecies = 0; ispecies < nspecies; ispecies++)
	y[ispecies] += h*(k1[ispecies]/6.0 + k2[ispecies]/3.0 + k3[ispecies]/3.0 + k4[ispecies]/6.0);

	} // end for (int step...

	// Store the solution back in atom->dvector.
	for (int ispecies = 0; ispecies < nspecies; ispecies++){
	if(y[ispecies] < -MY_EPSILON)
	error->one(FLERR,"Computed concentration in RK4 solver is < -10*DBL_EPSILON");
	else if(y[ispecies] < MY_EPSILON)
	y[ispecies] = 0.0;
	atom->dvector[ispecies][id] = y[ispecies];
	}
	}

	/* ---------------------------------------------------------------------- */

	// f1 = dt*f(t,x)
	// f2 = dtf(t+ c20dt,x + c21*f1)
	// f3 = dtf(t+ c30dt,x + c31f1 + c32f2)
	// f4 = dtf(t+ c40dt,x + c41f1 + c42f2 + c43*f3)
	// f5 = dtf(t+dt,x + c51f1 + c52f2 + c53f3 + c54*f4)
	// f6 = dtf(t+ c60dt,x + c61f1 + c62f2 + c63f3 + c64f4 + c65*f5)
	//
	// fifth-order runge-kutta integration
	// x5 = x + b1f1 + b3f3 + b4f4 + b5f5 + b6*f6
	// fourth-order runge-kutta integration
	// x = x + a1f1 + a3f3 + a4f4 + a5f5

	void FixRX::rkf45_step (const int neq, const double h, double y[], double y_out[], double rwk[], void* v_param)
	{
	const double c21=0.25;
	const double c31=0.09375;
	const double c32=0.28125;
	const double c41=0.87938097405553;
	const double c42=-3.2771961766045;
	const double c43=3.3208921256258;
	const double c51=2.0324074074074;
	const double c52=-8.0;
	const double c53=7.1734892787524;
	const double c54=-0.20589668615984;
	const double c61=-0.2962962962963;
	const double c62=2.0;
	const double c63=-1.3816764132554;
	const double c64=0.45297270955166;
	const double c65=-0.275;
	const double a1=0.11574074074074;
	const double a3=0.54892787524366;
	const double a4=0.5353313840156;
	const double a5=-0.2;
	const double b1=0.11851851851852;
	const double b3=0.51898635477583;
	const double b4=0.50613149034201;
	const double b5=-0.18;
	const double b6=0.036363636363636;

	// local dependent variables (5 total)
	double* f1 = &rwk[ 0];
	double* f2 = &rwk[ neq];
	double* f3 = &rwk[2*neq];
	double* f4 = &rwk[3*neq];
	double* f5 = &rwk[4*neq];
	double* f6 = &rwk[5*neq];

	// scratch for the intermediate solution.
	//double* ytmp = &rwk[6*neq];
	double* ytmp = y_out;

	// 1)
	rhs (0.0, y, f1, v_param);

	for (int k = 0; k < neq; k++){
	f1[k] *= h;
	ytmp[k] = y[k] + c21 * f1[k];
	}

	// 2)
	rhs(0.0, ytmp, f2, v_param);

	for (int k = 0; k < neq; k++){
	f2[k] *= h;
	ytmp[k] = y[k] + c31 * f1[k] + c32 * f2[k];
	}

	// 3)
	rhs(0.0, ytmp, f3, v_param);

	for (int k = 0; k < neq; k++) {
	f3[k] *= h;
	ytmp[k] = y[k] + c41 * f1[k] + c42 * f2[k] + c43 * f3[k];
	}

	// 4)
	rhs(0.0, ytmp, f4, v_param);

	for (int k = 0; k < neq; k++) {
	f4[k] *= h;
	ytmp[k] = y[k] + c51 * f1[k] + c52 * f2[k] + c53 * f3[k] + c54 * f4[k];
	}

	// 5)
	rhs(0.0, ytmp, f5, v_param);

	for (int k = 0; k < neq; k++) {
	f5[k] *= h;
	ytmp[k] = y[k] + c61f1[k] + c62f2[k] + c63f3[k] + c64f4[k] + c65*f5[k];
	}

	// 6)
	rhs(0.0, ytmp, f6, v_param);

	for (int k = 0; k < neq; k++)
	{
	//const double f6 = h * ydot[k];
	f6[k] *= h;

	// 5th-order solution.
	const double r5 = b1f1[k] + b3f3[k] + b4f4[k] + b5f5[k] + b6*f6[k];

	// 4th-order solution.
	const double r4 = a1f1[k] + a3f3[k] + a4f4[k] + a5f5[k];

	// Truncation error: difference between 4th and 5th-order solutions.
	rwk[k] = fabs(r5 - r4);

	// Update solution.
	//y_out[k] = y[k] + r5; // Local extrapolation
	y_out[k] = y[k] + r4;
	}

	return;
	}

	int FixRX::rkf45_h0 (const int neq, const double t, const double t_stop,
	const double hmin, const double hmax,
	double& h0, double y[], double rwk[], void* v_params)
	{
	// Set lower and upper bounds on h0, and take geometric mean as first trial value.
	// Exit with this value if the bounds cross each other.

	// Adjust upper bound based on ydot ...
	double hg = sqrt(hmin*hmax);

	//if (hmax < hmin)
	//{
	// h0 = hg;
	// return;
	//}

	// Start iteration to find solution to ... {WRMS norm of (h0^2 y'' / 2)} = 1

	double *ydot = rwk;
	double *y1 = ydot + neq;
	double *ydot1 = y1 + neq;

	const int max_iters = 10;
	bool hnew_is_ok = false;
	double hnew = hg;
	int iter = 0;

	// compute ydot at t=t0
	rhs (t, y, ydot, v_params);

	while(1)
	{
	// Estimate y'' with finite-difference ...

	for (int k = 0; k < neq; k++)
	y1[k] = y[k] + hg * ydot[k];

	// compute y' at t1
	rhs (t + hg, y1, ydot1, v_params);

	// Compute WRMS norm of y''
	double yddnrm = 0.0;
	for (int k = 0; k < neq; k++){
	double ydd = (ydot1[k] - ydot[k]) / hg;
	double wterr = ydd / (relTol * fabs( y[k] ) + absTol);
	yddnrm += wterr * wterr;
	}

	yddnrm = sqrt( yddnrm / double(neq) );

	//std::cout << "iter " << _iter << " hg " << hg << " y'' " << yddnrm << std::endl;
	//std::cout << "ydot " << ydot[neq-1] << std::endl;

	// should we accept this?
	if (hnew_is_ok \|\| iter == max_iters){
	hnew = hg;
	if (iter == max_iters)
	fprintf(stderr, "ERROR_HIN_MAX_ITERS\n");
	break;
	}

	// Get the new value of h ...
	hnew = (yddnrmhmaxhmax > 2.0) ? sqrt(2.0 / yddnrm) : sqrt(hg * hmax);

	// test the stopping conditions.
	double hrat = hnew / hg;

	// Accept this value ... the bias factor should bring it within range.
	if ( (hrat > 0.5) && (hrat < 2.0) )
	hnew_is_ok = true;

	// If y'' is still bad after a few iterations, just accept h and give up.
	if ( (iter > 1) && hrat > 2.0 ) {
	hnew = hg;
	hnew_is_ok = true;
	}

	//printf("iter=%d, yddnrw=%e, hnew=%e, hmin=%e, hmax=%e\n", iter, yddnrm, hnew, hmin, hmax);

	hg = hnew;
	iter ++;
	}

	// bound and bias estimate
	h0 = hnew * 0.5;
	h0 = fmax(h0, hmin);
	h0 = fmin(h0, hmax);
	//printf("h0=%e, hmin=%e, hmax=%e\n", h0, hmin, hmax);

	return (iter + 1);
	}

	void FixRX::odeDiagnostics(void)
	{
	TimerType timer_start = getTimeStamp();

	// Compute:
	// 1) Average # of ODE integrator steps and RHS evaluations per atom globally.
	// 2) RMS # of ...
	// 3) Average # of ODE steps and RHS evaluations per MPI task.
	// 4) RMS # of ODE steps and RHS evaluations per MPI task.
	// 5) MAX # of ODE steps and RHS evaluations per MPI task.
	//
	// ... 1,2 are for ODE control diagnostics.
	// ... 3-5 are for load balancing diagnostics.
	//
	// To do this, we'll need to
	// a) Allreduce (sum) the sum of nSteps / nFuncs. Dividing by atom->natoms
	// gives the avg # of steps/funcs per atom globally.
	// b) Reduce (sum) to root the sum of squares of the differences.
	// i) Sum_i (steps_i - avg_steps_global)^2
	// ii) Sum_i (funcs_i - avg_funcs_global)^2
	// iii) (avg_steps_local - avg_steps_global)^2
	// iv) (avg_funcs_local - avg_funcs_global)^2

	const int numCounters = numDiagnosticCounters-1;

	// # of time-steps for averaging.
	const int nTimes = this->diagnosticCounter[numDiagnosticCounters-1];

	// # of ODE's per time-step (on average).
	//const int nODEs = this->diagnosticCounter[AtomSum] / nTimes;

	// Sum up the sums from each task.
	double sums[numCounters];
	double my_vals[numCounters];
	double max_per_proc[numCounters];
	double min_per_proc[numCounters];

	if(1)
	{
	static bool firstStep = true;

	static TimerType oldTimeStamp (-1);

	TimerType now = getTimeStamp();

	// Query the fix database and look for rx_weight for the balance fix.
	int type_flag = -1;
	int rx_weight_index = atom->find_custom( "rx_weight", /0:int, 1:float/ type_flag );

	// Compute the average # of neighbors.
	double averageNumNeighbors = 0;
	{
	const int inum = pairDPDE->list->inum;
	const int* ilist = pairDPDE->list->ilist;
	const int* numneigh = pairDPDE->list->numneigh;

	for (int ii = 0; ii < inum; ++ii)
	{
	const int i = ilist[ii];
	averageNumNeighbors += numneigh[i];
	}

	averageNumNeighbors /= inum;
	}

	printf("me= %d nst= %g nfc= %g time= %g nlocal= %g lmpnst= %g weight_idx= %d 1st= %d aveNeigh= %g\n", comm->me, this->diagnosticCounter[0], this->diagnosticCounter[1], this->diagnosticCounter[2], this->diagnosticCounter[3], this->diagnosticCounter[4], rx_weight_index, firstStep, averageNumNeighbors);

	if (rx_weight_index != -1 && !firstStep && 0)
	{
	double *rx_weight = atom->dvector[rx_weight_index];

	const int nlocal = atom->nlocal;
	const int *mask = atom->mask;

	if (odeIntegrationFlag == ODE_LAMMPS_RKF45 && diagnosticFrequency == 1)
	{
	const double total_time = getElapsedTime( oldTimeStamp, now );
	const double fixrx_time = this->diagnosticCounter[TimeSum];
	const double time_ratio = fixrx_time / total_time;

	double tsum = 0.0;
	double tmin = 100000, tmax = 0;
	for (int i = 0; i < nlocal; ++i)
	if (mask[i] & groupbit)
	{
	double nfunc_ratio = double( diagnosticCounterPerODE[FuncSum][i] ) / diagnosticCounter[FuncSum];
	rx_weight[i] = nfunc_ratio * fixrx_time + (total_time - fixrx_time) / nlocal;
	tmin = fmin( tmin, rx_weight[i] );
	tmax = fmax( tmax, rx_weight[i] );
	tsum += rx_weight[i];
	//rx_weight[i] = (double) diagnosticCounterPerODE[FuncSum][i];
	}

	printf("me= %d total= %g fixrx= %g ratio= %g tsum= %g %g %g %g\n", comm->me, total_time, fixrx_time, time_ratio, tsum, (total_time - fixrx_time) / nlocal, tmin, tmax);
	}
	else
	{
	error->warning(FLERR, "Dynamic load balancing enabled but per-atom weights not available.");

	for (int i = 0; i < nlocal; ++i)
	if (mask[i] & groupbit)
	rx_weight[i] = 1.0;
	}
	}

	firstStep = false;
	oldTimeStamp = now;
	}

	// Compute counters per dpd time-step.
	for (int i = 0; i < numCounters; ++i){
	my_vals[i] = this->diagnosticCounter[i] / nTimes;
	//printf("my sum[%d] = %f %d\n", i, my_vals[i], comm->me);
	}

	MPI_Allreduce (my_vals, sums, numCounters, MPI_DOUBLE, MPI_SUM, world);

	MPI_Reduce (my_vals, max_per_proc, numCounters, MPI_DOUBLE, MPI_MAX, 0, world);
	MPI_Reduce (my_vals, min_per_proc, numCounters, MPI_DOUBLE, MPI_MIN, 0, world);

	const double nODEs = sums[numCounters-1];

	double avg_per_atom[numCounters], avg_per_proc[numCounters];

	// Averages per-ODE and per-proc per time-step.
	for (int i = 0; i < numCounters; ++i){
	avg_per_atom[i] = sums[i] / nODEs;
	avg_per_proc[i] = sums[i] / comm->nprocs;
	}

	// Sum up the differences from each task.
	double sum_sq[2*numCounters];
	double my_sum_sq[2*numCounters];
	for (int i = 0; i < numCounters; ++i){
	double diff_i = my_vals[i] - avg_per_proc[i];
	my_sum_sq[i] = diff_i * diff_i;
	}

	double max_per_ODE[numCounters], min_per_ODE[numCounters];

	// Process the per-ODE RMS of the # of steps/funcs
	if (diagnosticFrequency == 1){
	double my_max[numCounters], my_min[numCounters];

	const int nlocal = atom->nlocal;
	const int *mask = atom->mask;

	for (int i = 0; i < numCounters; ++i){
	my_sum_sq[i+numCounters] = 0;
	my_max[i] = 0;
	my_min[i] = DBL_MAX;

	if (diagnosticCounterPerODE[i] != NULL){
	for (int j = 0; j < nlocal; ++j)
	if (mask[j] & groupbit){
	double diff = double(diagnosticCounterPerODE[i][j]) - avg_per_atom[i];
	my_sum_sq[i+numCounters] += diff*diff;

	my_max[i] = std::max( my_max[i], (double)diagnosticCounterPerODE[i][j] );
	my_min[i] = std::min( my_min[i], (double)diagnosticCounterPerODE[i][j] );
	}
	}
	}

	MPI_Reduce (my_sum_sq, sum_sq, 2*numCounters, MPI_DOUBLE, MPI_SUM, 0, world);

	MPI_Reduce (my_max, max_per_ODE, numCounters, MPI_DOUBLE, MPI_MAX, 0, world);
	MPI_Reduce (my_min, min_per_ODE, numCounters, MPI_DOUBLE, MPI_MIN, 0, world);
	}
	else
	MPI_Reduce (my_sum_sq, sum_sq, numCounters, MPI_DOUBLE, MPI_SUM, 0, world);

	TimerType timer_stop = getTimeStamp();
	double time_local = getElapsedTime( timer_start, timer_stop );

	if (comm->me == 0){
	char smesg[128];

	#define print_mesg(smesg) {\
	if (screen) fprintf(screen,"%s\n", smesg); \
	if (logfile) fprintf(logfile,"%s\n", smesg); }

	sprintf(smesg, "FixRX::ODE Diagnostics: # of iters \|# of rhs evals\| run-time (sec) \| # atoms");
	print_mesg(smesg);

	sprintf(smesg, " AVG per ODE : %-12.5g \| %-12.5g \| %-12.5g", avg_per_atom[0], avg_per_atom[1], avg_per_atom[2]);
	print_mesg(smesg);

	// only valid for single time-step!
	if (diagnosticFrequency == 1){
	double rms_per_ODE[numCounters];
	for (int i = 0; i < numCounters; ++i)
	rms_per_ODE[i] = sqrt( sum_sq[i+numCounters] / nODEs );

	sprintf(smesg, " RMS per ODE : %-12.5g \| %-12.5g ", rms_per_ODE[0], rms_per_ODE[1]);
	print_mesg(smesg);

	sprintf(smesg, " MAX per ODE : %-12.5g \| %-12.5g ", max_per_ODE[0], max_per_ODE[1]);
	print_mesg(smesg);

	sprintf(smesg, " MIN per ODE : %-12.5g \| %-12.5g ", min_per_ODE[0], min_per_ODE[1]);
	print_mesg(smesg);
	}

	sprintf(smesg, " AVG per Proc : %-12.5g \| %-12.5g \| %-12.5g \| %-12.5g", avg_per_proc[StepSum], avg_per_proc[FuncSum], avg_per_proc[TimeSum], avg_per_proc[AtomSum]);
	print_mesg(smesg);

	if (comm->nprocs > 1){
	double rms_per_proc[numCounters];
	for (int i = 0; i < numCounters; ++i)
	rms_per_proc[i] = sqrt( sum_sq[i] / comm->nprocs );

	sprintf(smesg, " RMS per Proc : %-12.5g \| %-12.5g \| %-12.5g \| %-12.5g", rms_per_proc[0], rms_per_proc[1], rms_per_proc[2], rms_per_proc[AtomSum]);
	print_mesg(smesg);

	sprintf(smesg, " MAX per Proc : %-12.5g \| %-12.5g \| %-12.5g \| %-12.5g", max_per_proc[0], max_per_proc[1], max_per_proc[2], max_per_proc[AtomSum]);
	print_mesg(smesg);

	sprintf(smesg, " MIN per Proc : %-12.5g \| %-12.5g \| %-12.5g \| %-12.5g", min_per_proc[0], min_per_proc[1], min_per_proc[2], min_per_proc[AtomSum]);
	print_mesg(smesg);
	}

	sprintf(smesg, " AVG'd over %d time-steps", nTimes);
	print_mesg(smesg);
	sprintf(smesg, " AVG'ing took %g sec", time_local);
	print_mesg(smesg);

	#undef print_mesg

	}

	// Reset the counters.
	for (int i = 0; i < numDiagnosticCounters; ++i)
	diagnosticCounter[i] = 0;

	return;
	}

	void FixRX::rkf45(int id, double rwork, void v_param, int ode_counter[])
	{
	// Rounding coefficient.
	const double uround = DBL_EPSILON;

	// Adaption limit (shrink or grow)
	const double adaption_limit = 4.0;

	//double y = new double[8nspecies + nreactions];
	double *y = rwork;

	const int neq = nspecies;

	// Update ConcOld and initialize the ODE solution vector y[].
	for (int ispecies = 0; ispecies < nspecies; ispecies++){
	const double tmp = atom->dvector[ispecies][id];
	atom->dvector[ispecies+nspecies][id] = tmp;
	y[ispecies] = tmp;
	}

	// Integration length.
	const double t_stop = update->dt; // DPD time-step.

	// Safety factor on the adaption. very specific but not necessary .. 0.9 is common.
	const double hsafe = 0.840896415;

	// Time rounding factor.
	const double tround = t_stop * uround;

	// Counters for diagnostics.
	int nst = 0; // # of steps (accepted)
	int nit = 0; // # of iterations total
	int nfe = 0; // # of RHS evaluations

	// Min/Max step-size limits.
	const double h_min = 100.0 * tround;
	const double h_max = (minSteps > 0) ? t_stop / double(minSteps) : t_stop;

	// Set the initial step-size. 0 forces an internal estimate ... stable Euler step size.
	double h = (minSteps > 0) ? t_stop / double(minSteps) : 0.0;

	double t = 0.0;

	if (h < h_min){
	//fprintf(stderr,"hin not implemented yet\n");
	//exit(-1);
	nfe = rkf45_h0 (neq, t, t_stop, h_min, h_max, h, y, y + neq, v_param);
	}

	//printf("t= %e t_stop= %e h= %e\n", t, t_stop, h);

	// Integrate until we reach the end time.
	while (fabs(t - t_stop) > tround){
	double *yout = y + neq;
	double *eout = yout + neq;

	// Take a trial step.
	rkf45_step (neq, h, y, yout, eout, v_param);

	// Estimate the solution error.
	// ... weighted 2-norm of the error.
	double err2 = 0.0;
	for (int k = 0; k < neq; k++){
	const double wterr = eout[k] / (relTol * fabs( y[k] ) + absTol);
	err2 += wterr * wterr;
	}

	double err = fmax( uround, sqrt( err2 / double(nspecies) ));

	// Accept the solution?
	if (err <= 1.0 \|\| h <= h_min){
	t += h;
	nst++;

	for (int k = 0; k < neq; k++)
	y[k] = yout[k];
	}

	// Adjust h for the next step.
	double hfac = hsafe * sqrt( sqrt( 1.0 / err ) );

	// Limit the adaption.
	hfac = fmax( hfac, 1.0 / adaption_limit );
	hfac = fmin( hfac, adaption_limit );

	// Apply the adaption factor...
	h *= hfac;

	// Limit h.
	h = fmin( h, h_max );
	h = fmax( h, h_min );

	// Stretch h if we're within 5% ... and we didn't just fail.
	if (err <= 1.0 && (t + 1.05*h) > t_stop)
	h = t_stop - t;

	// And don't overshoot the end.
	if (t + h > t_stop)
	h = t_stop - t;

	nit++;
	nfe += 6;

	if (maxIters && nit > maxIters){
	//fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop);
	//nFails ++;
	ode_counter[3] ++;
	break;
	// We should set an error here so that the solution is not used!
	}

	} // end while

	ode_counter[0] += nst;
	ode_counter[1] += nit;
	ode_counter[2] += nfe;

	//if (diagnosticFrequency == 1 && diagnosticCounterPerODE[StepSum] != NULL)
	if (diagnosticCounterPerODE[StepSum] != NULL){
	diagnosticCounterPerODE[StepSum][id] = nst;
	diagnosticCounterPerODE[FuncSum][id] = nfe;
	}
	//printf("id= %d nst= %d nit= %d\n", id, nst, nit);

	// Store the solution back in atom->dvector.
	for (int ispecies = 0; ispecies < nspecies; ispecies++){
	if(y[ispecies] < -1.0e-10)
	error->one(FLERR,"Computed concentration in RKF45 solver is < -1.0e-10");
	else if(y[ispecies] < MY_EPSILON)
	y[ispecies] = 0.0;
	atom->dvector[ispecies][id] = y[ispecies];
	}
	}

	/* ---------------------------------------------------------------------- */

	int FixRX::rhs(double t, const double y, double dydt, void *params)
	{
	// Use the sparse format instead.
	if (useSparseKinetics)
	return this->rhs_sparse( t, y, dydt, params);
	else
	return this->rhs_dense ( t, y, dydt, params);
	}

	/* ---------------------------------------------------------------------- */

	int FixRX::rhs_dense(double t, const double y, double dydt, void *params)
	{
	UserRHSData userData = (UserRHSData ) params;

	double *rxnRateLaw = userData->rxnRateLaw;
	double *kFor = userData->kFor;

	const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;
	const int nspecies = atom->nspecies_dpd;

	for(int ispecies=0; ispecies<nspecies; ispecies++)
	dydt[ispecies] = 0.0;

	// Construct the reaction rate laws
	for(int jrxn=0; jrxn<nreactions; jrxn++){
	double rxnRateLawForward = kFor[jrxn];

	for(int ispecies=0; ispecies<nspecies; ispecies++){
	const double concentration = y[ispecies]/VDPD;
	rxnRateLawForward *= pow(concentration,stoichReactants[jrxn][ispecies]);
	}
	rxnRateLaw[jrxn] = rxnRateLawForward;
	}

	// Construct the reaction rates for each species
	for(int ispecies=0; ispecies<nspecies; ispecies++)
	for(int jrxn=0; jrxn<nreactions; jrxn++)
	dydt[ispecies] += stoich[jrxn][ispecies]VDPDrxnRateLaw[jrxn];

	return 0;
	}

	/* ---------------------------------------------------------------------- */

	int FixRX::rhs_sparse(double t, const double y, double dydt, void *v_params) const
	{
	UserRHSData userData = (UserRHSData ) v_params;

	const double VDPD = domain->xprd * domain->yprd * domain->zprd / atom->natoms;

	#define kFor (userData->kFor)
	#define kRev (NULL)
	#define rxnRateLaw (userData->rxnRateLaw)
	#define conc (dydt)
	#define maxReactants (this->sparseKinetics_maxReactants)
	#define maxSpecies (this->sparseKinetics_maxSpecies)
	#define nuk (this->sparseKinetics_nuk)
	#define nu (this->sparseKinetics_nu)
	#define inu (this->sparseKinetics_inu)
	#define isIntegral(idx) (SparseKinetics_enableIntegralReactions \
	&& this->sparseKinetics_isIntegralReaction[idx])

	for (int k = 0; k < nspecies; ++k)
	conc[k] = y[k] / VDPD;

	// Construct the reaction rate laws
	for (int i = 0; i < nreactions; ++i)
	{
	double rxnRateLawForward;
	if (isIntegral(i)){
	rxnRateLawForward = kFor[i] * powint( conc[ nuk[i][0] ], inu[i][0]);
	for (int kk = 1; kk < maxReactants; ++kk){
	const int k = nuk[i][kk];
	if (k == SparseKinetics_invalidIndex) break;
	//if (k != SparseKinetics_invalidIndex)
	rxnRateLawForward *= powint( conc[k], inu[i][kk] );
	}
	} else {
	rxnRateLawForward = kFor[i] * pow( conc[ nuk[i][0] ], nu[i][0]);
	for (int kk = 1; kk < maxReactants; ++kk){
	const int k = nuk[i][kk];
	if (k == SparseKinetics_invalidIndex) break;
	//if (k != SparseKinetics_invalidIndex)
	rxnRateLawForward *= pow( conc[k], nu[i][kk] );
	}
	}

	rxnRateLaw[i] = rxnRateLawForward;
	}

	// Construct the reaction rates for each species from the
	// Stoichiometric matrix and ROP vector.
	for (int k = 0; k < nspecies; ++k)
	dydt[k] = 0.0;

	for (int i = 0; i < nreactions; ++i){
	// Reactants ...
	dydt[ nuk[i][0] ] -= nu[i][0] * rxnRateLaw[i];
	for (int kk = 1; kk < maxReactants; ++kk){
	const int k = nuk[i][kk];
	if (k == SparseKinetics_invalidIndex) break;
	//if (k != SparseKinetics_invalidIndex)
	dydt[k] -= nu[i][kk] * rxnRateLaw[i];
	}

	// Products ...
	dydt[ nuk[i][maxReactants] ] += nu[i][maxReactants] * rxnRateLaw[i];
	for (int kk = maxReactants+1; kk < maxSpecies; ++kk){
	const int k = nuk[i][kk];
	if (k == SparseKinetics_invalidIndex) break;
	//if (k != SparseKinetics_invalidIndex)
	dydt[k] += nu[i][kk] * rxnRateLaw[i];
	}
	}

	// Add in the volume factor to convert to the proper units.
	for (int k = 0; k < nspecies; ++k)
	dydt[k] *= VDPD;

	#undef kFor
	#undef kRev
	#undef rxnRateLaw
	#undef conc
	#undef maxReactants
	#undef maxSpecies
	#undef nuk
	#undef nu
	#undef inu
	#undef isIntegral
	//#undef invalidIndex

	return 0;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::computeLocalTemperature()
	{
	int i,j,ii,jj,inum,jnum,itype,jtype;
	double xtmp,ytmp,ztmp,delx,dely,delz;
	double rsq;
	int ilist,jlist,numneigh,*firstneigh;

	double **x = atom->x;
	int *type = atom->type;
	int nlocal = atom->nlocal;
	int nghost = atom->nghost;
	int newton_pair = force->newton_pair;

	// local temperature variables
	double wij=0.0;
	double *dpdTheta = atom->dpdTheta;

	// Initialize the local temperature weight array
	int sumWeightsCt = nlocal + (newton_pair ? nghost : 0);
	sumWeights = new double[sumWeightsCt];
	memset(sumWeights, 0, sizeof(double)*sumWeightsCt);

	inum = list->inum;
	ilist = list->ilist;
	numneigh = list->numneigh;
	firstneigh = list->firstneigh;

	// loop over neighbors of my atoms
	for (ii = 0; ii < inum; ii++) {
	i = ilist[ii];
	xtmp = x[i][0];
	ytmp = x[i][1];
	ztmp = x[i][2];
	itype = type[i];
	jlist = firstneigh[i];
	jnum = numneigh[i];

	for (jj = 0; jj < jnum; jj++) {
	j = jlist[jj];
	j &= NEIGHMASK;
	jtype = type[j];

	delx = xtmp - x[j][0];
	dely = ytmp - x[j][1];
	delz = ztmp - x[j][2];
	rsq = delxdelx + delydely + delz*delz;

	if (rsq < pairDPDE->cutsq[itype][jtype]) {
	double rcut = sqrt(pairDPDE->cutsq[itype][jtype]);
	double rij = sqrt(rsq);
	double ratio = rij/rcut;

	// Lucy's Weight Function
	if(wtFlag==LUCY){
	wij = (1.0+3.0ratio) (1.0-ratio)(1.0-ratio)(1.0-ratio);
	dpdThetaLocal[i] += wij/dpdTheta[j];
	if (newton_pair \|\| j < nlocal)
	dpdThetaLocal[j] += wij/dpdTheta[i];
	}

	sumWeights[i] += wij;
	if (newton_pair \|\| j < nlocal)
	sumWeights[j] += wij;
	}
	}
	}
	if (newton_pair) comm->reverse_comm_fix(this);

	// self-interaction for local temperature
	for (i = 0; i < nlocal; i++){

	// Lucy Weight Function
	if(wtFlag==LUCY){
	wij = 1.0;
	dpdThetaLocal[i] += wij / dpdTheta[i];
	}
	sumWeights[i] += wij;

	// Normalized local temperature
	dpdThetaLocal[i] = dpdThetaLocal[i] / sumWeights[i];

	if(localTempFlag == HARMONIC)
	dpdThetaLocal[i] = 1.0 / dpdThetaLocal[i];

	}

	delete [] sumWeights;
	}

	/* ---------------------------------------------------------------------- */

	int FixRX::pack_forward_comm(int n, int list, double buf, int pbc_flag, int *pbc)
	{
	int ii,jj,m;
	double tmp;

	m = 0;
	for (ii = 0; ii < n; ii++) {
	jj = list[ii];
	for(int ispecies=0;ispecies<nspecies;ispecies++){
	tmp = atom->dvector[ispecies][jj];
	buf[m++] = tmp;
	tmp = atom->dvector[ispecies+nspecies][jj];
	buf[m++] = tmp;
	}
	}
	return m;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::unpack_forward_comm(int n, int first, double *buf)
	{
	int ii,m,last;
	double tmp;

	m = 0;
	last = first + n ;
	for (ii = first; ii < last; ii++){
	for(int ispecies=0;ispecies<nspecies;ispecies++){
	tmp = buf[m++];
	atom->dvector[ispecies][ii] = tmp;
	tmp = buf[m++];
	atom->dvector[ispecies+nspecies][ii] = tmp;
	}
	}
	}

	/* ---------------------------------------------------------------------- */

	int FixRX::pack_reverse_comm(int n, int first, double *buf)
	{
	int i,m,last;

	m = 0;
	last = first + n;
	for (i = first; i < last; i++) {
	buf[m++] = dpdThetaLocal[i];
	buf[m++] = sumWeights[i];
	}
	return m;
	}

	/* ---------------------------------------------------------------------- */

	void FixRX::unpack_reverse_comm(int n, int list, double buf)
	{
	int i,j,m;

	m = 0;
	for (i = 0; i < n; i++) {
	j = list[i];

	dpdThetaLocal[j] += buf[m++];
	sumWeights[j] += buf[m++];
	}
	}

fix_rx.cppNo OneTemporaryActions

File Metadata

fix_rx.cppView Options

Event Timeline

fix_rx.cpp
No OneTemporary
Actions

fix_rx.cpp
View Options