diff --git a/src/USER-DPD/fix_eos_table_rx.cpp b/src/USER-DPD/fix_eos_table_rx.cpp
index 9a67afd72afa350b8572d0cf0056355df4d5ff20..1dfd58c4658b0372692327562c689f0617951d20 100644
--- a/src/USER-DPD/fix_eos_table_rx.cpp
+++ b/src/USER-DPD/fix_eos_table_rx.cpp
@@ -28,6 +28,12 @@
#define MAXLINE 1024
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
using namespace LAMMPS_NS;
using namespace FixConst;
@@ -37,17 +43,18 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), ntables(0), tables(NULL),
tables2(NULL), dHf(NULL), eosSpecies(NULL)
{
- if (narg != 8) error->all(FLERR,"Illegal fix eos/table/rx command");
+ if (narg != 8 && narg != 10) error->all(FLERR,"Illegal fix eos/table/rx command");
restart_peratom = 1;
nevery = 1;
- bool rx_flag = false;
+ rx_flag = false;
+ nspecies = 1;
for (int i = 0; i < modify->nfix; i++)
- if (strncmp(modify->fix[i]->style,"rx",2) == 0) rx_flag = true;
- if (!rx_flag) error->all(FLERR,"FixEOStableRX requires a fix rx command.");
-
- nspecies = atom->nspecies_dpd;
- if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
+ if (strncmp(modify->fix[i]->style,"rx",2) == 0){
+ rx_flag = true;
+ nspecies = atom->nspecies_dpd;
+ if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
+ }
if (strcmp(arg[3],"linear") == 0) tabstyle = LINEAR;
else error->all(FLERR,"Unknown table style in fix eos/table/rx");
@@ -113,8 +120,25 @@ FixEOStableRX::FixEOStableRX(LAMMPS *lmp, int narg, char **arg) :
ntables++;
}
- // Read the Formation Enthalpies
- read_file(arg[7]);
+ // Read the Formation Enthalpies and Correction Coefficients
+ dHf = new double[nspecies];
+ energyCorr = new double[nspecies];
+ tempCorrCoeff = new double[nspecies];
+ moleculeCorrCoeff= new double[nspecies];
+ for (int ii=0; ii<nspecies; ii++){
+ dHf[ii] = 0.0;
+ energyCorr[ii] = 0.0;
+ tempCorrCoeff[ii] = 0.0;
+ moleculeCorrCoeff[ii] = 0.0;
+ }
+
+ if(rx_flag) read_file(arg[7]);
+ else dHf[0] = atof(arg[7]);
+
+ if(narg==10){
+ energyCorr[0] = atof(arg[8]);
+ tempCorrCoeff[0] = atof(arg[9]);
+ }
comm_forward = 3;
comm_reverse = 2;
@@ -136,6 +160,9 @@ FixEOStableRX::~FixEOStableRX()
delete [] dHf;
delete [] eosSpecies;
+ delete [] energyCorr;
+ delete [] tempCorrCoeff;
+ delete [] moleculeCorrCoeff;
}
/* ---------------------------------------------------------------------- */
@@ -269,10 +296,9 @@ void FixEOStableRX::end_of_step()
void FixEOStableRX::read_file(char *file)
{
- int params_per_line = 2;
- char **words = new char*[params_per_line+1];
-
- dHf = new double[nspecies];
+ int min_params_per_line = 2;
+ int max_params_per_line = 5;
+ char **words = new char*[max_params_per_line+1];
// open file on proc 0
@@ -313,7 +339,7 @@ void FixEOStableRX::read_file(char *file)
// concatenate additional lines until have params_per_line words
- while (nwords < params_per_line) {
+ while (nwords < min_params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
@@ -330,7 +356,7 @@ void FixEOStableRX::read_file(char *file)
nwords = atom->count_words(line);
}
- if (nwords != params_per_line)
+ if (nwords != min_params_per_line && nwords != max_params_per_line)
error->all(FLERR,"Incorrect format in eos table/rx potential file");
// words = ptrs to all words in line
@@ -342,8 +368,14 @@ void FixEOStableRX::read_file(char *file)
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(words[0],&atom->dname[ispecies][0]) == 0) break;
- if (ispecies < nspecies)
+ if (ispecies < nspecies){
dHf[ispecies] = atof(words[1]);
+ if(nwords > min_params_per_line+1){
+ energyCorr[ispecies] = atof(words[2]);
+ tempCorrCoeff[ispecies] = atof(words[3]);
+ moleculeCorrCoeff[ispecies] = atof(words[4]);
+ }
+ }
}
delete [] words;
@@ -545,27 +577,33 @@ void FixEOStableRX::param_extract(Table *tb, char *line)
error->one(FLERR,"Invalid keyword in fix eos/table/rx parameters");
word = strtok(NULL," \t\n\r\f");
- while (word) {
- for (ispecies = 0; ispecies < nspecies; ispecies++)
- if (strcmp(word,&atom->dname[ispecies][0]) == 0){
- eosSpecies[ncolumn] = ispecies;
- ncolumn++;
- break;
+ if(rx_flag){
+ while (word) {
+ for (ispecies = 0; ispecies < nspecies; ispecies++)
+ if (strcmp(word,&atom->dname[ispecies][0]) == 0){
+ eosSpecies[ncolumn] = ispecies;
+ ncolumn++;
+ break;
+ }
+ if (ispecies == nspecies){
+ printf("name=%s not found in species list\n",word);
+ error->one(FLERR,"Invalid keyword in fix eos/table/rx parameters");
}
- if (ispecies == nspecies){
- printf("name=%s not found in species list\n",word);
- error->one(FLERR,"Invalid keyword in fix eos/table/rx parameters");
+ word = strtok(NULL," \t\n\r\f");
}
- word = strtok(NULL," \t\n\r\f");
- }
- for (int icolumn = 0; icolumn < ncolumn; icolumn++)
- if(eosSpecies[icolumn]==-1)
- error->one(FLERR,"EOS data is missing from fix eos/table/rx tabe");
- if(ncolumn != nspecies){
- printf("ncolumns=%d nspecies=%d\n",ncolumn,nspecies);
- error->one(FLERR,"The number of columns in fix eos/table/rx does not match the number of species");
+ for (int icolumn = 0; icolumn < ncolumn; icolumn++)
+ if(eosSpecies[icolumn]==-1)
+ error->one(FLERR,"EOS data is missing from fix eos/table/rx tabe");
+ if(ncolumn != nspecies){
+ printf("ncolumns=%d nspecies=%d\n",ncolumn,nspecies);
+ error->one(FLERR,"The number of columns in fix eos/table/rx does not match the number of species");
+ }
+ } else {
+ eosSpecies[0] = 0;
+ ncolumn++;
}
+
if (tb->ninput == 0) error->one(FLERR,"fix eos/table/rx parameters did not set N");
}
@@ -653,11 +691,27 @@ double FixEOStableRX::splint(double *xa, double *ya, double *y2a, int n, double
void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
{
- int itable;
- double fraction, uTmp, nTotal;
+ int itable, nPG;
+ double fraction, uTmp, nMolecules, nTotal, nTotalPG;
+ double tolerance = 1.0e-10;
ui = 0.0;
nTotal = 0.0;
+ nTotalPG = 0.0;
+ nPG = 0;
+
+ if(rx_flag){
+ for(int ispecies=0;ispecies<nspecies;ispecies++){
+ nTotal += atom->dvector[ispecies][id];
+ if(fabs(moleculeCorrCoeff[ispecies]) > tolerance){
+ nPG++;
+ nTotalPG += atom->dvector[ispecies][id];
+ }
+ }
+ } else {
+ nTotal = 1.0;
+ }
+
for(int ispecies=0;ispecies<nspecies;ispecies++){
Table *tb = &tables[ispecies];
thetai = MAX(thetai,tb->lo);
@@ -669,9 +723,13 @@ void FixEOStableRX::energy_lookup(int id, double thetai, double &ui)
uTmp = tb->e[itable] + fraction*tb->de[itable];
uTmp += dHf[ispecies];
- // mol fraction form:
- ui += atom->dvector[ispecies][id]*uTmp;
- nTotal += atom->dvector[ispecies][id];
+ uTmp += tempCorrCoeff[ispecies]*thetai; // temperature correction
+ uTmp += energyCorr[ispecies]; // energy correction
+ if(nPG > 0) ui += moleculeCorrCoeff[ispecies]*nTotalPG/double(nPG); // molecule correction
+
+ if(rx_flag) nMolecules = atom->dvector[ispecies][id];
+ else nMolecules = 1.0;
+ ui += nMolecules*uTmp;
}
}
ui = ui - double(nTotal+1.5)*force->boltz*thetai;
@@ -690,6 +748,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
double maxit = 100;
double temp;
double delta = 0.001;
+ double tolerance = 1.0e-10;
// Store the current thetai in t1
t1 = MAX(thetai,tb->lo);
@@ -713,7 +772,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
// Apply the Secant Method
for(it=0; it<maxit; it++){
- if(fabs(f2-f1)<1e-15){
+ if(fabs(f2-f1) < MY_EPSILON){
if(isnan(f1) || isnan(f2)) error->one(FLERR,"NaN detected in secant solver.");
temp = t1;
temp = MAX(temp,tb->lo);
@@ -724,7 +783,7 @@ void FixEOStableRX::temperature_lookup(int id, double ui, double &thetai)
break;
}
temp = t2 - f2*(t2-t1)/(f2-f1);
- if(fabs(temp-t2) < 1e-6) break;
+ if(fabs(temp-t2) < tolerance) break;
f1 = f2;
t1 = t2;
t2 = temp;
diff --git a/src/USER-DPD/fix_eos_table_rx.h b/src/USER-DPD/fix_eos_table_rx.h
index 078cf1e2e1c71ea38eacee620ccc98f4a0123290..8c26d133a52e364bf34bce6c4d7a7b9009adde57 100644
--- a/src/USER-DPD/fix_eos_table_rx.h
+++ b/src/USER-DPD/fix_eos_table_rx.h
@@ -67,7 +67,7 @@ class FixEOStableRX : public Fix {
void read_file(char *);
- double *dHf;
+ double *dHf,*energyCorr,*tempCorrCoeff,*moleculeCorrCoeff;
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
@@ -76,6 +76,7 @@ class FixEOStableRX : public Fix {
int *eosSpecies;
int ncolumn;
+ bool rx_flag;
};
}
diff --git a/src/USER-DPD/pair_exp6_rx.cpp b/src/USER-DPD/pair_exp6_rx.cpp
index a1add9a4d2d60ad0a64ca7572fe641bcc66b098f..7815a81b17e853bfb88f553658c713449fd303e8 100644
--- a/src/USER-DPD/pair_exp6_rx.cpp
+++ b/src/USER-DPD/pair_exp6_rx.cpp
@@ -35,6 +35,12 @@ using namespace MathSpecial;
#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 oneFluidApproxParameter (-1)
#define isOneFluidApprox(_site) ( (_site) == oneFluidApproxParameter )
@@ -47,17 +53,17 @@ using namespace MathSpecial;
struct PairExp6ParamDataType
{
int n;
- double *epsilon1, *alpha1, *rm1, *fraction1,
- *epsilon2, *alpha2, *rm2, *fraction2,
- *epsilonOld1, *alphaOld1, *rmOld1, *fractionOld1,
- *epsilonOld2, *alphaOld2, *rmOld2, *fractionOld2;
+ double *epsilon1, *alpha1, *rm1, *mixWtSite1,
+ *epsilon2, *alpha2, *rm2, *mixWtSite2,
+ *epsilonOld1, *alphaOld1, *rmOld1, *mixWtSite1old,
+ *epsilonOld2, *alphaOld2, *rmOld2, *mixWtSite2old;
// Default constructor -- nullify everything.
PairExp6ParamDataType(void)
- : n(0), epsilon1(NULL), alpha1(NULL), rm1(NULL), fraction1(NULL),
- epsilon2(NULL), alpha2(NULL), rm2(NULL), fraction2(NULL),
- epsilonOld1(NULL), alphaOld1(NULL), rmOld1(NULL), fractionOld1(NULL),
- epsilonOld2(NULL), alphaOld2(NULL), rmOld2(NULL), fractionOld2(NULL)
+ : n(0), epsilon1(NULL), alpha1(NULL), rm1(NULL), mixWtSite1(NULL),
+ epsilon2(NULL), alpha2(NULL), rm2(NULL), mixWtSite2(NULL),
+ epsilonOld1(NULL), alphaOld1(NULL), rmOld1(NULL), mixWtSite1old(NULL),
+ epsilonOld2(NULL), alphaOld2(NULL), rmOld2(NULL), mixWtSite2old(NULL)
{}
};
@@ -71,6 +77,7 @@ PairExp6rx::PairExp6rx(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
mol2param = NULL;
+ fractionalWeighting = true;
}
/* ---------------------------------------------------------------------- */
@@ -89,6 +96,11 @@ PairExp6rx::~PairExp6rx()
memory->destroy(cutsq);
memory->destroy(cut);
}
+ if(scalingFlag == POLYNOMIAL){
+ memory->destroy(coeffAlpha);
+ memory->destroy(coeffEps);
+ memory->destroy(coeffRm);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -130,10 +142,10 @@ void PairExp6rx::compute(int eflag, int vflag)
double epsilon2_j,alpha2_j,rm2_j;
double evdwlOldEXP6_12, evdwlOldEXP6_21, fpairOldEXP6_12, fpairOldEXP6_21;
double evdwlEXP6_12, evdwlEXP6_21;
- double fractionOld1_i, fractionOld1_j;
- double fractionOld2_i, fractionOld2_j;
- double fraction1_i, fraction1_j;
- double fraction2_i, fraction2_j;
+ double mixWtSite1old_i, mixWtSite1old_j;
+ double mixWtSite2old_i, mixWtSite2old_j;
+ double mixWtSite1_i, mixWtSite1_j;
+ double mixWtSite2_i, mixWtSite2_j;
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
@@ -153,38 +165,38 @@ void PairExp6rx::compute(int eflag, int vflag)
memory->create( PairExp6ParamData.epsilon1 , np_total, "PairExp6ParamData.epsilon1");
memory->create( PairExp6ParamData.alpha1 , np_total, "PairExp6ParamData.alpha1");
memory->create( PairExp6ParamData.rm1 , np_total, "PairExp6ParamData.rm1");
- memory->create( PairExp6ParamData.fraction1 , np_total, "PairExp6ParamData.fraction1");
+ memory->create( PairExp6ParamData.mixWtSite1 , np_total, "PairExp6ParamData.mixWtSite1");
memory->create( PairExp6ParamData.epsilon2 , np_total, "PairExp6ParamData.epsilon2");
memory->create( PairExp6ParamData.alpha2 , np_total, "PairExp6ParamData.alpha2");
memory->create( PairExp6ParamData.rm2 , np_total, "PairExp6ParamData.rm2");
- memory->create( PairExp6ParamData.fraction2 , np_total, "PairExp6ParamData.fraction2");
+ memory->create( PairExp6ParamData.mixWtSite2 , np_total, "PairExp6ParamData.mixWtSite2");
memory->create( PairExp6ParamData.epsilonOld1 , np_total, "PairExp6ParamData.epsilonOld1");
memory->create( PairExp6ParamData.alphaOld1 , np_total, "PairExp6ParamData.alphaOld1");
memory->create( PairExp6ParamData.rmOld1 , np_total, "PairExp6ParamData.rmOld1");
- memory->create( PairExp6ParamData.fractionOld1 , np_total, "PairExp6ParamData.fractionOld1");
+ memory->create( PairExp6ParamData.mixWtSite1old , np_total, "PairExp6ParamData.mixWtSite1old");
memory->create( PairExp6ParamData.epsilonOld2 , np_total, "PairExp6ParamData.epsilonOld2");
memory->create( PairExp6ParamData.alphaOld2 , np_total, "PairExp6ParamData.alphaOld2");
memory->create( PairExp6ParamData.rmOld2 , np_total, "PairExp6ParamData.rmOld2");
- memory->create( PairExp6ParamData.fractionOld2 , np_total, "PairExp6ParamData.fractionOld2");
+ memory->create( PairExp6ParamData.mixWtSite2old , np_total, "PairExp6ParamData.mixWtSite2old");
for (i = 0; i < np_total; ++i)
{
- getParamsEXP6 (i, PairExp6ParamData.epsilon1[i],
+ getMixingWeights (i, PairExp6ParamData.epsilon1[i],
PairExp6ParamData.alpha1[i],
PairExp6ParamData.rm1[i],
- PairExp6ParamData.fraction1[i],
+ PairExp6ParamData.mixWtSite1[i],
PairExp6ParamData.epsilon2[i],
PairExp6ParamData.alpha2[i],
PairExp6ParamData.rm2[i],
- PairExp6ParamData.fraction2[i],
+ PairExp6ParamData.mixWtSite2[i],
PairExp6ParamData.epsilonOld1[i],
PairExp6ParamData.alphaOld1[i],
PairExp6ParamData.rmOld1[i],
- PairExp6ParamData.fractionOld1[i],
+ PairExp6ParamData.mixWtSite1old[i],
PairExp6ParamData.epsilonOld2[i],
PairExp6ParamData.alphaOld2[i],
PairExp6ParamData.rmOld2[i],
- PairExp6ParamData.fractionOld2[i]);
+ PairExp6ParamData.mixWtSite2old[i]);
}
}
@@ -208,19 +220,19 @@ void PairExp6rx::compute(int eflag, int vflag)
epsilon1_i = PairExp6ParamData.epsilon1[i];
alpha1_i = PairExp6ParamData.alpha1[i];
rm1_i = PairExp6ParamData.rm1[i];
- fraction1_i = PairExp6ParamData.fraction1[i];
+ mixWtSite1_i = PairExp6ParamData.mixWtSite1[i];
epsilon2_i = PairExp6ParamData.epsilon2[i];
alpha2_i = PairExp6ParamData.alpha2[i];
rm2_i = PairExp6ParamData.rm2[i];
- fraction2_i = PairExp6ParamData.fraction2[i];
+ mixWtSite2_i = PairExp6ParamData.mixWtSite2[i];
epsilonOld1_i = PairExp6ParamData.epsilonOld1[i];
alphaOld1_i = PairExp6ParamData.alphaOld1[i];
rmOld1_i = PairExp6ParamData.rmOld1[i];
- fractionOld1_i = PairExp6ParamData.fractionOld1[i];
+ mixWtSite1old_i = PairExp6ParamData.mixWtSite1old[i];
epsilonOld2_i = PairExp6ParamData.epsilonOld2[i];
alphaOld2_i = PairExp6ParamData.alphaOld2[i];
rmOld2_i = PairExp6ParamData.rmOld2[i];
- fractionOld2_i = PairExp6ParamData.fractionOld2[i];
+ mixWtSite2old_i = PairExp6ParamData.mixWtSite2old[i];
}
for (jj = 0; jj < jnum; jj++) {
@@ -255,19 +267,19 @@ void PairExp6rx::compute(int eflag, int vflag)
epsilon1_j = PairExp6ParamData.epsilon1[j];
alpha1_j = PairExp6ParamData.alpha1[j];
rm1_j = PairExp6ParamData.rm1[j];
- fraction1_j = PairExp6ParamData.fraction1[j];
+ mixWtSite1_j = PairExp6ParamData.mixWtSite1[j];
epsilon2_j = PairExp6ParamData.epsilon2[j];
alpha2_j = PairExp6ParamData.alpha2[j];
rm2_j = PairExp6ParamData.rm2[j];
- fraction2_j = PairExp6ParamData.fraction2[j];
+ mixWtSite2_j = PairExp6ParamData.mixWtSite2[j];
epsilonOld1_j = PairExp6ParamData.epsilonOld1[j];
alphaOld1_j = PairExp6ParamData.alphaOld1[j];
rmOld1_j = PairExp6ParamData.rmOld1[j];
- fractionOld1_j = PairExp6ParamData.fractionOld1[j];
+ mixWtSite1old_j = PairExp6ParamData.mixWtSite1old[j];
epsilonOld2_j = PairExp6ParamData.epsilonOld2[j];
alphaOld2_j = PairExp6ParamData.alphaOld2[j];
rmOld2_j = PairExp6ParamData.rmOld2[j];
- fractionOld2_j = PairExp6ParamData.fractionOld2[j];
+ mixWtSite2old_j = PairExp6ParamData.mixWtSite2old[j];
}
// A2. Apply Lorentz-Berthelot mixing rules for the i-j pair
@@ -368,9 +380,9 @@ void PairExp6rx::compute(int eflag, int vflag)
}
if (isite1 == isite2)
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12;
+ evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwlOldEXP6_12;
else
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwlOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*evdwlOldEXP6_21;
+ evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwlOldEXP6_12 + sqrt(mixWtSite2old_i*mixWtSite1old_j)*evdwlOldEXP6_21;
evdwlOld *= factor_lj;
@@ -451,8 +463,8 @@ void PairExp6rx::compute(int eflag, int vflag)
//
// Apply Mixing Rule to get the overall force for the CG pair
//
- if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12;
- else fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpairOldEXP6_12 + sqrt(fractionOld2_i*fractionOld1_j)*fpairOldEXP6_21;
+ if (isite1 == isite2) fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpairOldEXP6_12;
+ else fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpairOldEXP6_12 + sqrt(mixWtSite2old_i*mixWtSite1old_j)*fpairOldEXP6_21;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
@@ -463,8 +475,8 @@ void PairExp6rx::compute(int eflag, int vflag)
f[j][2] -= delz*fpair;
}
- if (isite1 == isite2) evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12;
- else evdwl = sqrt(fraction1_i*fraction2_j)*evdwlEXP6_12 + sqrt(fraction2_i*fraction1_j)*evdwlEXP6_21;
+ if (isite1 == isite2) evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwlEXP6_12;
+ else evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwlEXP6_12 + sqrt(mixWtSite2_i*mixWtSite1_j)*evdwlEXP6_21;
evdwl *= factor_lj;
uCGnew[i] += 0.5*evdwl;
@@ -484,19 +496,19 @@ void PairExp6rx::compute(int eflag, int vflag)
if (PairExp6ParamData.epsilon1 ) memory->destroy(PairExp6ParamData.epsilon1);
if (PairExp6ParamData.alpha1 ) memory->destroy(PairExp6ParamData.alpha1);
if (PairExp6ParamData.rm1 ) memory->destroy(PairExp6ParamData.rm1);
- if (PairExp6ParamData.fraction1 ) memory->destroy(PairExp6ParamData.fraction1);
+ if (PairExp6ParamData.mixWtSite1 ) memory->destroy(PairExp6ParamData.mixWtSite1);
if (PairExp6ParamData.epsilon2 ) memory->destroy(PairExp6ParamData.epsilon2);
if (PairExp6ParamData.alpha2 ) memory->destroy(PairExp6ParamData.alpha2);
if (PairExp6ParamData.rm2 ) memory->destroy(PairExp6ParamData.rm2);
- if (PairExp6ParamData.fraction2 ) memory->destroy(PairExp6ParamData.fraction2);
+ if (PairExp6ParamData.mixWtSite2 ) memory->destroy(PairExp6ParamData.mixWtSite2);
if (PairExp6ParamData.epsilonOld1 ) memory->destroy(PairExp6ParamData.epsilonOld1);
if (PairExp6ParamData.alphaOld1 ) memory->destroy(PairExp6ParamData.alphaOld1);
if (PairExp6ParamData.rmOld1 ) memory->destroy(PairExp6ParamData.rmOld1);
- if (PairExp6ParamData.fractionOld1) memory->destroy(PairExp6ParamData.fractionOld1);
+ if (PairExp6ParamData.mixWtSite1old) memory->destroy(PairExp6ParamData.mixWtSite1old);
if (PairExp6ParamData.epsilonOld2 ) memory->destroy(PairExp6ParamData.epsilonOld2);
if (PairExp6ParamData.alphaOld2 ) memory->destroy(PairExp6ParamData.alphaOld2);
if (PairExp6ParamData.rmOld2 ) memory->destroy(PairExp6ParamData.rmOld2);
- if (PairExp6ParamData.fractionOld2) memory->destroy(PairExp6ParamData.fractionOld2);
+ if (PairExp6ParamData.mixWtSite2old) memory->destroy(PairExp6ParamData.mixWtSite2old);
}
}
@@ -526,10 +538,20 @@ void PairExp6rx::allocate()
void PairExp6rx::settings(int narg, char **arg)
{
- if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+ if (narg < 1) error->all(FLERR,"Illegal pair_style command");
cut_global = force->numeric(FLERR,arg[0]);
+ // optional keywords
+
+ int iarg = 1;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
+ else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
+ else error->all(FLERR,"Illegal pair_style command");
+ iarg++;
+ }
+
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
@@ -547,7 +569,7 @@ void PairExp6rx::settings(int narg, char **arg)
void PairExp6rx::coeff(int narg, char **arg)
{
- if (narg < 7 || narg > 8) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg < 6 || narg > 9) error->all(FLERR,"Incorrect args for pair coefficients");
bool rx_flag = false;
for (int i = 0; i < modify->nfix; i++)
@@ -624,21 +646,36 @@ void PairExp6rx::coeff(int narg, char **arg)
params[iparam].potentialType = exp6PotentialType;
else
error->all(FLERR,"params[].potential type unknown");
-
- //printf("params[%d].name= %s ispecies= %d potential= %s potentialType= %d\n", iparam, params[iparam].name, params[iparam].ispecies, params[iparam].potential, params[iparam].potentialType);
}
}
delete[] site1;
delete[] site2;
site1 = site2 = NULL;
- fuchslinR = force->numeric(FLERR,arg[5]);
- fuchslinEpsilon = force->numeric(FLERR,arg[6]);
-
setup();
double cut_one = cut_global;
- if (narg == 8) cut_one = force->numeric(FLERR,arg[7]);
+ if (strcmp(arg[5],"exponent") == 0){
+ scalingFlag = EXPONENT;
+ exponentR = force->numeric(FLERR,arg[6]);
+ exponentEpsilon = force->numeric(FLERR,arg[7]);
+ if (narg > 9) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg == 9) cut_one = force->numeric(FLERR,arg[8]);
+ } else if (strcmp(arg[5],"polynomial") == 0){
+ scalingFlag = POLYNOMIAL;
+ memory->create(coeffAlpha,6,"pair:coeffAlpha");
+ memory->create(coeffEps,6,"pair:coeffEps");
+ memory->create(coeffRm,6,"pair:coeffRm");
+ read_file2(arg[6]);
+ if (narg > 8) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg == 8) cut_one = force->numeric(FLERR,arg[7]);
+ } else if (strcmp(arg[5],"none") == 0){
+ scalingFlag = NONE;
+ if (narg > 7) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg == 7) cut_one = force->numeric(FLERR,arg[6]);
+ } else {
+ error->all(FLERR,"Incorrect args for pair coefficients");
+ }
int count = 0;
for (int i = ilo; i <= ihi; i++) {
@@ -780,6 +817,95 @@ void PairExp6rx::read_file(char *file)
/* ---------------------------------------------------------------------- */
+void PairExp6rx::read_file2(char *file)
+{
+ int params_per_line = 7;
+ char **words = new char*[params_per_line+1];
+
+ // open file on proc 0
+
+ FILE *fp;
+ fp = NULL;
+ if (comm->me == 0) {
+ fp = fopen(file,"r");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open polynomial file %s",file);
+ error->one(FLERR,str);
+ }
+ }
+
+ // one set of params can span multiple lines
+ int n,nwords;
+ 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;
+
+ // concatenate additional lines until have params_per_line words
+
+ while (nwords < params_per_line) {
+ n = strlen(line);
+ if (comm->me == 0) {
+ ptr = fgets(&line[n],MAXLINE-n,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);
+ if ((ptr = strchr(line,'#'))) *ptr = '\0';
+ nwords = atom->count_words(line);
+ }
+
+ if (nwords != params_per_line)
+ error->all(FLERR,"Incorrect format in polynomial file");
+
+ // words = ptrs to all words in line
+
+ nwords = 0;
+ words[nwords++] = strtok(line," \t\n\r\f");
+ while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
+
+ if (strcmp(words[0],"alpha") == 0){
+ for (int ii=1; ii<params_per_line; ii++)
+ coeffAlpha[ii-1] = atof(words[ii]);
+ }
+ if (strcmp(words[0],"epsilon") == 0){
+ for (int ii=1; ii<params_per_line; ii++)
+ coeffEps[ii-1] = atof(words[ii]);
+ }
+ if (strcmp(words[0],"rm") == 0){
+ for (int ii=1; ii<params_per_line; ii++)
+ coeffRm[ii-1] = atof(words[ii]);
+ }
+ }
+
+ delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
void PairExp6rx::setup()
{
int i,j,n;
@@ -877,7 +1003,7 @@ void PairExp6rx::read_restart_settings(FILE *fp)
/* ---------------------------------------------------------------------- */
-void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm1, double &fraction1,double &epsilon2,double &alpha2,double &rm2,double &fraction2,double &epsilon1_old,double &alpha1_old,double &rm1_old, double &fraction1_old,double &epsilon2_old,double &alpha2_old,double &rm2_old,double &fraction2_old) const
+void PairExp6rx::getMixingWeights(int id,double &epsilon1,double &alpha1,double &rm1, double &mixWtSite1,double &epsilon2,double &alpha2,double &rm2,double &mixWtSite2,double &epsilon1_old,double &alpha1_old,double &rm1_old, double &mixWtSite1old,double &epsilon2_old,double &alpha2_old,double &rm2_old,double &mixWtSite2old) const
{
int iparam, jparam;
double rmi, rmj, rmij, rm3ij;
@@ -885,11 +1011,16 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
double alphai, alphaj, alphaij;
double epsilon_old, rm3_old, alpha_old;
double epsilon, rm3, alpha;
- double fractionOFA, fractionOFA_old;
- double nTotalOFA, nTotalOFA_old;
- double nTotal, nTotal_old;
double xMolei, xMolej, xMolei_old, xMolej_old;
+ double fractionOFAold, fractionOFA;
+ double fractionOld1, fraction1;
+ double fractionOld2, fraction2;
+ double nMoleculesOFAold, nMoleculesOFA;
+ double nMoleculesOld1, nMolecules1;
+ double nMoleculesOld2, nMolecules2;
+ double nTotal, nTotalOld;
+
rm3 = 0.0;
epsilon = 0.0;
alpha = 0.0;
@@ -897,32 +1028,32 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rm3_old = 0.0;
alpha_old = 0.0;
fractionOFA = 0.0;
- fractionOFA_old = 0.0;
- nTotalOFA = 0.0;
- nTotalOFA_old = 0.0;
+ fractionOFAold = 0.0;
+ nMoleculesOFA = 0.0;
+ nMoleculesOFAold = 0.0;
nTotal = 0.0;
- nTotal_old = 0.0;
+ nTotalOld = 0.0;
// Compute the total number of molecules in the old and new CG particle as well as the total number of molecules in the fluid portion of the old and new CG particle
for (int ispecies = 0; ispecies < nspecies; ispecies++){
nTotal += atom->dvector[ispecies][id];
- nTotal_old += atom->dvector[ispecies+nspecies][id];
+ nTotalOld += atom->dvector[ispecies+nspecies][id];
iparam = mol2param[ispecies];
if (iparam < 0 || params[iparam].potentialType != exp6PotentialType ) continue;
if (isOneFluidApprox(isite1) || isOneFluidApprox(isite2)) {
if (isite1 == params[iparam].ispecies || isite2 == params[iparam].ispecies) continue;
- nTotalOFA_old += atom->dvector[ispecies+nspecies][id];
- nTotalOFA += atom->dvector[ispecies][id];
+ nMoleculesOFAold += atom->dvector[ispecies+nspecies][id];
+ nMoleculesOFA += atom->dvector[ispecies][id];
}
}
- if(nTotal < 1e-8 || nTotal_old < 1e-8)
- error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
+ if(nTotal < MY_EPSILON || nTotalOld < MY_EPSILON)
+ error->all(FLERR,"The number of molecules in CG particle is less than 10*DBL_EPSILON.");
// Compute the mole fraction of molecules within the fluid portion of the particle (One Fluid Approximation)
- fractionOFA_old = nTotalOFA_old / nTotal_old;
- fractionOFA = nTotalOFA / nTotal;
+ fractionOFAold = nMoleculesOFAold / nTotalOld;
+ fractionOFA = nMoleculesOFA / nTotal;
for (int ispecies = 0; ispecies < nspecies; ispecies++) {
iparam = mol2param[ispecies];
@@ -938,8 +1069,10 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
alpha1 = params[iparam].alpha;
// Compute the mole fraction of Site1
- fraction1_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
- fraction1 = atom->dvector[ispecies][id]/nTotal;
+ nMoleculesOld1 = atom->dvector[ispecies+nspecies][id];
+ nMolecules1 = atom->dvector[ispecies][id];
+ fractionOld1 = nMoleculesOld1/nTotalOld;
+ fraction1 = nMolecules1/nTotal;
}
// If Site2 matches a pure species, then grab the parameters
@@ -952,7 +1085,9 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
alpha2 = params[iparam].alpha;
// Compute the mole fraction of Site2
- fraction2_old = atom->dvector[ispecies+nspecies][id]/nTotal_old;
+ nMoleculesOld2 = atom->dvector[ispecies+nspecies][id];
+ nMolecules2 = atom->dvector[ispecies][id];
+ fractionOld2 = atom->dvector[ispecies+nspecies][id]/nTotalOld;
fraction2 = atom->dvector[ispecies][id]/nTotal;
}
@@ -962,8 +1097,10 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rmi = params[iparam].rm;
epsiloni = params[iparam].epsilon;
alphai = params[iparam].alpha;
- xMolei = atom->dvector[ispecies][id]/nTotalOFA;
- xMolei_old = atom->dvector[ispecies+nspecies][id]/nTotalOFA_old;
+ if(nMoleculesOFA<MY_EPSILON) xMolei = 0.0;
+ else xMolei = atom->dvector[ispecies][id]/nMoleculesOFA;
+ if(nMoleculesOFAold<MY_EPSILON) xMolei_old = 0.0;
+ else xMolei_old = atom->dvector[ispecies+nspecies][id]/nMoleculesOFAold;
for (int jspecies = 0; jspecies < nspecies; jspecies++) {
jparam = mol2param[jspecies];
@@ -972,15 +1109,17 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
rmj = params[jparam].rm;
epsilonj = params[jparam].epsilon;
alphaj = params[jparam].alpha;
- xMolej = atom->dvector[jspecies][id]/nTotalOFA;
- xMolej_old = atom->dvector[jspecies+nspecies][id]/nTotalOFA_old;
+ if(nMoleculesOFA<MY_EPSILON) xMolej = 0.0;
+ else xMolej = atom->dvector[jspecies][id]/nMoleculesOFA;
+ if(nMoleculesOFAold<MY_EPSILON) xMolej_old = 0.0;
+ else xMolej_old = atom->dvector[jspecies+nspecies][id]/nMoleculesOFAold;
rmij = (rmi+rmj)/2.0;
rm3ij = rmij*rmij*rmij;
epsilonij = sqrt(epsiloni*epsilonj);
alphaij = sqrt(alphai*alphaj);
- if(fractionOFA_old > 0.0){
+ if(fractionOFAold > 0.0){
rm3_old += xMolei_old*xMolej_old*rm3ij;
epsilon_old += xMolei_old*xMolej_old*rm3ij*epsilonij;
alpha_old += xMolei_old*xMolej_old*rm3ij*epsilonij*alphaij;
@@ -996,7 +1135,7 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
if (isOneFluidApprox(isite1)){
rm1 = cbrt(rm3);
- if(rm1 < 1e-16) {
+ if(rm1 < MY_EPSILON) {
rm1 = 0.0;
epsilon1 = 0.0;
alpha1 = 0.0;
@@ -1004,11 +1143,11 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon1 = epsilon / rm3;
alpha1 = alpha / epsilon1 / rm3;
}
-
+ nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
fraction1 = fractionOFA;
rm1_old = cbrt(rm3_old);
- if(rm1_old < 1e-16) {
+ if(rm1_old < MY_EPSILON) {
rm1_old = 0.0;
epsilon1_old = 0.0;
alpha1_old = 0.0;
@@ -1016,42 +1155,21 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon1_old = epsilon_old / rm3_old;
alpha1_old = alpha_old / epsilon1_old / rm3_old;
}
- fraction1_old = fractionOFA_old;
-
- // Fuchslin-Like Exp-6 Scaling
- double powfuch = 0.0;
- if(fuchslinEpsilon < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon1 = 0.0;
- else epsilon1 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon1_old = 0.0;
- else epsilon1_old *= 1.0/powfuch;
-
- } else {
- epsilon1 *= pow(nTotalOFA,fuchslinEpsilon);
- epsilon1_old *= pow(nTotalOFA_old,fuchslinEpsilon);
- }
-
- if(fuchslinR < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinR);
- if(powfuch<1e-15) rm1 = 0.0;
- else rm1 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinR);
- if(powfuch<1e-15) rm1_old = 0.0;
- else rm1_old *= 1.0/powfuch;
-
- } else {
- rm1 *= pow(nTotalOFA,fuchslinR);
- rm1_old *= pow(nTotalOFA_old,fuchslinR);
+ nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
+ fractionOld1 = fractionOFAold;
+
+ if(scalingFlag == EXPONENT){
+ exponentScaling(nMoleculesOFA,epsilon1,rm1);
+ exponentScaling(nMoleculesOFAold,epsilon1_old,rm1_old);
+ } else if(scalingFlag == POLYNOMIAL){
+ polynomialScaling(nMoleculesOFA,alpha1,epsilon1,rm1);
+ polynomialScaling(nMoleculesOFAold,alpha1_old,epsilon1_old,rm1_old);
}
}
if (isOneFluidApprox(isite2)){
rm2 = cbrt(rm3);
- if(rm2 < 1e-16) {
+ if(rm2 < MY_EPSILON) {
rm2 = 0.0;
epsilon2 = 0.0;
alpha2 = 0.0;
@@ -1059,10 +1177,11 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon2 = epsilon / rm3;
alpha2 = alpha / epsilon2 / rm3;
}
+ nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
fraction2 = fractionOFA;
rm2_old = cbrt(rm3_old);
- if(rm2_old < 1e-16) {
+ if(rm2_old < MY_EPSILON) {
rm2_old = 0.0;
epsilon2_old = 0.0;
alpha2_old = 0.0;
@@ -1070,64 +1189,96 @@ void PairExp6rx::getParamsEXP6(int id,double &epsilon1,double &alpha1,double &rm
epsilon2_old = epsilon_old / rm3_old;
alpha2_old = alpha_old / epsilon2_old / rm3_old;
}
- fraction2_old = fractionOFA_old;
-
- // Fuchslin-Like Exp-6 Scaling
- double powfuch = 0.0;
- if(fuchslinEpsilon < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon2 = 0.0;
- else epsilon2 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinEpsilon);
- if(powfuch<1e-15) epsilon2_old = 0.0;
- else epsilon2_old *= 1.0/powfuch;
-
- } else {
- epsilon2 *= pow(nTotalOFA,fuchslinEpsilon);
- epsilon2_old *= pow(nTotalOFA_old,fuchslinEpsilon);
- }
-
- if(fuchslinR < 0.0){
- powfuch = pow(nTotalOFA,-fuchslinR);
- if(powfuch<1e-15) rm2 = 0.0;
- else rm2 *= 1.0/powfuch;
-
- powfuch = pow(nTotalOFA_old,-fuchslinR);
- if(powfuch<1e-15) rm2_old = 0.0;
- else rm2_old *= 1.0/powfuch;
-
- } else {
- rm2 *= pow(nTotalOFA,fuchslinR);
- rm2_old *= pow(nTotalOFA_old,fuchslinR);
+ nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
+ fractionOld2 = fractionOFAold;
+
+ if(scalingFlag == EXPONENT){
+ exponentScaling(nMoleculesOFA,epsilon2,rm2);
+ exponentScaling(nMoleculesOFAold,epsilon2_old,rm2_old);
+ } else if(scalingFlag == POLYNOMIAL){
+ polynomialScaling(nMoleculesOFA,alpha2,epsilon2,rm2);
+ polynomialScaling(nMoleculesOFAold,alpha2_old,epsilon2_old,rm2_old);
}
}
// Check that no fractions are less than zero
- if(fraction1 < 0.0){
- if(fraction1 < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fraction1 < 0.0 || nMolecules1 < 0.0){
+ if(fraction1 < -MY_EPSILON || nMolecules1 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
+ nMolecules1 = 0.0;
fraction1 = 0.0;
}
- if(fraction2 < 0.0){
- if(fraction2 < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fraction2 < 0.0 || nMolecules2 < 0.0){
+ if(fraction2 < -MY_EPSILON || nMolecules2 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
+ nMolecules2 = 0.0;
fraction2 = 0.0;
}
- if(fraction1_old < 0.0){
- if(fraction1_old < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fractionOld1 < 0.0 || nMoleculesOld1 < 0.0){
+ if(fractionOld1 < -MY_EPSILON || nMoleculesOld1 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
- fraction1_old = 0.0;
+ nMoleculesOld1 = 0.0;
+ fractionOld1 = 0.0;
}
- if(fraction2_old < 0.0){
- if(fraction2_old < -1.0e-10){
- error->all(FLERR,"Computed fraction less than -1.0e-10");
+ if(fractionOld2 < 0.0 || nMoleculesOld2 < 0.0){
+ if(fractionOld2 < -MY_EPSILON || nMoleculesOld2 < -MY_EPSILON){
+ error->all(FLERR,"Computed fraction less than -10*DBL_EPSILON");
}
- fraction2_old = 0.0;
+ nMoleculesOld2 = 0.0;
+ fractionOld2 = 0.0;
}
+
+ if(fractionalWeighting){
+ mixWtSite1old = fractionOld1;
+ mixWtSite1 = fraction1;
+ mixWtSite2old = fractionOld2;
+ mixWtSite2 = fraction2;
+ } else {
+ mixWtSite1old = nMoleculesOld1;
+ mixWtSite1 = nMolecules1;
+ mixWtSite2old = nMoleculesOld2;
+ mixWtSite2 = nMolecules2;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairExp6rx::exponentScaling(double phi, double &epsilon, double &rm) const
+{
+ double powfuch;
+
+ if(exponentEpsilon < 0.0){
+ powfuch = pow(phi,-exponentEpsilon);
+ if(powfuch<MY_EPSILON) epsilon = 0.0;
+ else epsilon *= 1.0/powfuch;
+ } else {
+ epsilon *= pow(phi,exponentEpsilon);
+ }
+
+ if(exponentR < 0.0){
+ powfuch = pow(phi,-exponentR);
+ if(powfuch<MY_EPSILON) rm = 0.0;
+ else rm *= 1.0/powfuch;
+ } else {
+ rm *= pow(phi,exponentR);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairExp6rx::polynomialScaling(double phi, double &alpha, double &epsilon, double &rm) const
+{
+ double phi2 = phi*phi;
+ double phi3 = phi2*phi;
+ double phi4 = phi2*phi2;
+ double phi5 = phi2*phi3;
+
+ alpha = coeffAlpha[0]*phi5 + coeffAlpha[1]*phi4 + coeffAlpha[2]*phi3 + coeffAlpha[3]*phi2 + coeffAlpha[4]*phi + coeffAlpha[5];
+ epsilon = coeffEps[0]*phi5 + coeffEps[1]*phi4 + coeffEps[2]*phi3 + coeffEps[3]*phi2 + coeffEps[4]*phi + coeffEps[5];
+ rm = coeffRm[0]*phi5 + coeffRm[1]*phi4 + coeffRm[2]*phi3 + coeffRm[3]*phi2 + coeffRm[4]*phi + coeffRm[5];
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-DPD/pair_exp6_rx.h b/src/USER-DPD/pair_exp6_rx.h
index dd9fa22a48096f5b8ce3a7b3d07aa3b2b83aa2f2..c16875c96021702ecb532aed737b688daf8d7303 100644
--- a/src/USER-DPD/pair_exp6_rx.h
+++ b/src/USER-DPD/pair_exp6_rx.h
@@ -39,6 +39,7 @@ class PairExp6rx : public Pair {
protected:
enum{LINEAR};
+ enum{NONE,EXPONENT,POLYNOMIAL};
double cut_global;
double **cut;
double **epsilon,**rm,**alpha;
@@ -59,12 +60,18 @@ class PairExp6rx : public Pair {
int nspecies;
void read_file(char *);
+ void read_file2(char *);
void setup();
int isite1, isite2;
char *site1, *site2;
- double fuchslinR, fuchslinEpsilon;
- void getParamsEXP6(int, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &) const;
+ void getMixingWeights(int, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &, double &) const;
+ double exponentR, exponentEpsilon;
+ int scalingFlag;
+ void exponentScaling(double, double &, double &) const;
+ void polynomialScaling(double, double &, double &, double &) const;
+ double *coeffAlpha, *coeffEps, *coeffRm;
+ bool fractionalWeighting;
inline double func_rin(const double &) const;
inline double expValue(const double) const;
@@ -128,7 +135,7 @@ E: Potential file has duplicate entry.
Self-explanatory
-E: The number of molecules in CG particle is less than 1e-8.
+E: The number of molecules in CG particle is less than 10*DBL_EPSILON.
Self-explanatory. Check the species concentrations have been properly set
and check the reaction kinetic solver parameters in fix rx to more for
diff --git a/src/USER-DPD/pair_multi_lucy_rx.cpp b/src/USER-DPD/pair_multi_lucy_rx.cpp
index cd107f1519e8f5918eb6ecc71ad50f4086c207d0..6e20f206d5b81a114e10e3997044bb6b60adb876 100644
--- a/src/USER-DPD/pair_multi_lucy_rx.cpp
+++ b/src/USER-DPD/pair_multi_lucy_rx.cpp
@@ -43,6 +43,12 @@ enum{NONE,RLINEAR,RSQ};
#define MAXLINE 1024
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
#define oneFluidParameter (-1)
#define isOneFluid(_site) ( (_site) == oneFluidParameter )
@@ -72,6 +78,7 @@ PairMultiLucyRX::PairMultiLucyRX(LAMMPS *lmp) : Pair(lmp),
comm_forward = 1;
comm_reverse = 1;
+ fractionalWeighting = true;
}
/* ---------------------------------------------------------------------- */
@@ -112,9 +119,9 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
int nghost = atom->nghost;
int newton_pair = force->newton_pair;
- double fractionOld1_i,fractionOld1_j;
- double fractionOld2_i,fractionOld2_j;
- double fraction1_i;
+ double mixWtSite1old_i,mixWtSite1old_j;
+ double mixWtSite2old_i,mixWtSite2old_j;
+ double mixWtSite1_i;
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
@@ -124,20 +131,20 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
int jtable;
double *rho = atom->rho;
- double *fractionOld1 = NULL;
- double *fractionOld2 = NULL;
- double *fraction1 = NULL;
- double *fraction2 = NULL;
+ double *mixWtSite1old = NULL;
+ double *mixWtSite2old = NULL;
+ double *mixWtSite1 = NULL;
+ double *mixWtSite2 = NULL;
{
const int ntotal = nlocal + nghost;
- memory->create(fractionOld1, ntotal, "PairMultiLucyRX::fractionOld1");
- memory->create(fractionOld2, ntotal, "PairMultiLucyRX::fractionOld2");
- memory->create(fraction1, ntotal, "PairMultiLucyRX::fraction1");
- memory->create(fraction2, ntotal, "PairMultiLucyRX::fraction2");
+ memory->create(mixWtSite1old, ntotal, "PairMultiLucyRX::mixWtSite1old");
+ memory->create(mixWtSite2old, ntotal, "PairMultiLucyRX::mixWtSite2old");
+ memory->create(mixWtSite1, ntotal, "PairMultiLucyRX::mixWtSite1");
+ memory->create(mixWtSite2, ntotal, "PairMultiLucyRX::mixWtSite2");
for (int i = 0; i < ntotal; ++i)
- getParams(i, fractionOld1[i], fractionOld2[i], fraction1[i], fraction2[i]);
+ getMixingWeights(i, mixWtSite1old[i], mixWtSite2old[i], mixWtSite1[i], mixWtSite2[i]);
}
inum = list->inum;
@@ -162,9 +169,9 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
double fy_i = 0.0;
double fz_i = 0.0;
- fractionOld1_i = fractionOld1[i];
- fractionOld2_i = fractionOld2[i];
- fraction1_i = fraction1[i];
+ mixWtSite1old_i = mixWtSite1old[i];
+ mixWtSite2old_i = mixWtSite2old[i];
+ mixWtSite1_i = mixWtSite1[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@@ -179,8 +186,8 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) {
fpair = 0.0;
- fractionOld1_j = fractionOld1[j];
- fractionOld2_j = fractionOld2[j];
+ mixWtSite1old_j = mixWtSite1old[j];
+ mixWtSite2old_j = mixWtSite2old[j];
tb = &tables[tabindex[itype][jtype]];
if (rho[i]*rho[i] < tb->innersq || rho[j]*rho[j] < tb->innersq){
@@ -235,8 +242,8 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
- if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
- else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
+ if (isite1 == isite2) fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpair;
+ else fpair = (sqrt(mixWtSite1old_i*mixWtSite2old_j) + sqrt(mixWtSite2old_i*mixWtSite1old_j))*fpair;
fx_i += delx*fpair;
fy_i += dely*fpair;
@@ -268,8 +275,8 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
} else error->one(FLERR,"Only LOOKUP and LINEAR table styles have been implemented for pair multi/lucy/rx");
evdwl *=(pi*cutsq[itype][itype]*cutsq[itype][itype])/84.0;
- evdwlOld = fractionOld1_i*evdwl;
- evdwl = fraction1_i*evdwl;
+ evdwlOld = mixWtSite1old_i*evdwl;
+ evdwl = mixWtSite1_i*evdwl;
uCG[i] += evdwlOld;
uCGnew[i] += evdwl;
@@ -281,10 +288,10 @@ void PairMultiLucyRX::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
- memory->destroy(fractionOld1);
- memory->destroy(fractionOld2);
- memory->destroy(fraction1);
- memory->destroy(fraction2);
+ memory->destroy(mixWtSite1old);
+ memory->destroy(mixWtSite2old);
+ memory->destroy(mixWtSite1);
+ memory->destroy(mixWtSite2);
}
/* ----------------------------------------------------------------------
@@ -311,7 +318,7 @@ void PairMultiLucyRX::allocate()
void PairMultiLucyRX::settings(int narg, char **arg)
{
- if (narg != 2) error->all(FLERR,"Illegal pair_style command");
+ if (narg < 2) error->all(FLERR,"Illegal pair_style command");
// new settings
@@ -322,6 +329,16 @@ void PairMultiLucyRX::settings(int narg, char **arg)
tablength = force->inumeric(FLERR,arg[1]);
if (tablength < 2) error->all(FLERR,"Illegal number of pair table entries");
+ // optional keywords
+
+ int iarg = 2;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
+ else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
+ else error->all(FLERR,"Illegal pair_style command");
+ iarg++;
+ }
+
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
@@ -928,9 +945,14 @@ void PairMultiLucyRX::computeLocalDensity()
/* ---------------------------------------------------------------------- */
-void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOld2, double &fraction1, double &fraction2)
+void PairMultiLucyRX::getMixingWeights(int id, double &mixWtSite1old, double &mixWtSite2old, double &mixWtSite1, double &mixWtSite2)
{
- double fractionOld, fraction;
+ double fractionOFAold, fractionOFA;
+ double fractionOld1, fraction1;
+ double fractionOld2, fraction2;
+ double nMoleculesOFAold, nMoleculesOFA;
+ double nMoleculesOld1, nMolecules1;
+ double nMoleculesOld2, nMolecules2;
double nTotal, nTotalOld;
nTotal = 0.0;
@@ -941,32 +963,56 @@ void PairMultiLucyRX::getParams(int id, double &fractionOld1, double &fractionOl
}
if (isOneFluid(isite1) == false){
- fractionOld1 = atom->dvector[isite1+nspecies][id]/nTotalOld;
- fraction1 = atom->dvector[isite1][id]/nTotal;
+ nMoleculesOld1 = atom->dvector[isite1+nspecies][id];
+ nMolecules1 = atom->dvector[isite1][id];
+ fractionOld1 = nMoleculesOld1/nTotalOld;
+ fraction1 = nMolecules1/nTotal;
}
if (isOneFluid(isite2) == false){
- fractionOld2 = atom->dvector[isite2+nspecies][id]/nTotalOld;
- fraction2 = atom->dvector[isite2][id]/nTotal;
+ nMoleculesOld2 = atom->dvector[isite2+nspecies][id];
+ nMolecules2 = atom->dvector[isite2][id];
+ fractionOld2 = nMoleculesOld2/nTotalOld;
+ fraction2 = nMolecules2/nTotal;
}
if (isOneFluid(isite1) || isOneFluid(isite2)){
- fractionOld = 0.0;
- fraction = 0.0;
+ nMoleculesOFAold = 0.0;
+ nMoleculesOFA = 0.0;
+ fractionOFAold = 0.0;
+ fractionOFA = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if (isite1 == ispecies || isite2 == ispecies) continue;
- fractionOld += atom->dvector[ispecies+nspecies][id] / nTotalOld;
- fraction += atom->dvector[ispecies][id] / nTotal;
+ nMoleculesOFAold += atom->dvector[ispecies+nspecies][id];
+ nMoleculesOFA += atom->dvector[ispecies][id];
+ fractionOFAold += atom->dvector[ispecies+nspecies][id] / nTotalOld;
+ fractionOFA += atom->dvector[ispecies][id] / nTotal;
}
if (isOneFluid(isite1)){
- fractionOld1 = fractionOld;
- fraction1 = fraction;
+ nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld1 = fractionOFAold;
+ fraction1 = fractionOFA;
}
if (isOneFluid(isite2)){
- fractionOld2 = fractionOld;
- fraction2 = fraction;
+ nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld2 = fractionOFAold;
+ fraction2 = fractionOFA;
}
}
+
+ if(fractionalWeighting){
+ mixWtSite1old = fractionOld1;
+ mixWtSite1 = fraction1;
+ mixWtSite2old = fractionOld2;
+ mixWtSite2 = fraction2;
+ } else {
+ mixWtSite1old = nMoleculesOld1;
+ mixWtSite1 = nMolecules1;
+ mixWtSite2old = nMoleculesOld2;
+ mixWtSite2 = nMolecules2;
+ }
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-DPD/pair_multi_lucy_rx.h b/src/USER-DPD/pair_multi_lucy_rx.h
index 2913716c5a5ddce63670b869ff8064cc78ea78a3..87d2ed7ae8194f9fa0f03dfd85fdc0312461458d 100644
--- a/src/USER-DPD/pair_multi_lucy_rx.h
+++ b/src/USER-DPD/pair_multi_lucy_rx.h
@@ -78,7 +78,8 @@ class PairMultiLucyRX : public Pair {
int nspecies;
char *site1, *site2;
int isite1, isite2;
- void getParams(int, double &, double &, double &, double &);
+ void getMixingWeights(int, double &, double &, double &, double &);
+ bool fractionalWeighting;
};
diff --git a/src/USER-DPD/pair_table_rx.cpp b/src/USER-DPD/pair_table_rx.cpp
index 902d0e5bb4397144599637f8c03c977b7784ba6c..8498d752e31396b79101c362568beb06659d8b24 100644
--- a/src/USER-DPD/pair_table_rx.cpp
+++ b/src/USER-DPD/pair_table_rx.cpp
@@ -35,6 +35,12 @@ enum{NONE,RLINEAR,RSQ,BMP};
#define MAXLINE 1024
+#ifdef DBL_EPSILON
+ #define MY_EPSILON (10.0*DBL_EPSILON)
+#else
+ #define MY_EPSILON (10.0*2.220446049250313e-16)
+#endif
+
#define OneFluidValue (-1)
#define isOneFluid(_site_) ( (_site_) == OneFluidValue )
@@ -44,6 +50,7 @@ PairTableRX::PairTableRX(LAMMPS *lmp) : Pair(lmp)
{
ntables = 0;
tables = NULL;
+ fractionalWeighting = true;
}
/* ---------------------------------------------------------------------- */
@@ -82,21 +89,6 @@ void PairTableRX::compute(int eflag, int vflag)
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
- double *fractionOld1, *fractionOld2;
- double *fraction1, *fraction2;
-
- {
- const int ntotal = atom->nlocal + atom->nghost;
-
- memory->create(fractionOld1, ntotal, "PairTableRx::compute::fractionOld1");
- memory->create(fractionOld2, ntotal, "PairTableRx::compute::fractionOld2");
- memory->create(fraction1, ntotal, "PairTableRx::compute::fraction1");
- memory->create(fraction2, ntotal, "PairTableRx::compute::fraction2");
-
- for (int i = 0; i < ntotal; ++i)
- getParams(i, fractionOld1[i], fractionOld2[i], fraction1[i], fraction2[i]);
- }
-
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
@@ -104,13 +96,29 @@ void PairTableRX::compute(int eflag, int vflag)
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
- double fractionOld1_i, fractionOld1_j;
- double fractionOld2_i, fractionOld2_j;
- double fraction1_i, fraction1_j;
- double fraction2_i, fraction2_j;
+ double mixWtSite1old_i, mixWtSite1old_j;
+ double mixWtSite2old_i, mixWtSite2old_j;
+ double mixWtSite1_i, mixWtSite1_j;
+ double mixWtSite2_i, mixWtSite2_j;
double *uCG = atom->uCG;
double *uCGnew = atom->uCGnew;
+ double *mixWtSite1old = NULL;
+ double *mixWtSite2old = NULL;
+ double *mixWtSite1 = NULL;
+ double *mixWtSite2 = NULL;
+
+ {
+ const int ntotal = atom->nlocal + atom->nghost;
+ memory->create(mixWtSite1old, ntotal, "PairTableRx::compute::mixWtSite1old");
+ memory->create(mixWtSite2old, ntotal, "PairTableRx::compute::mixWtSite2old");
+ memory->create(mixWtSite1, ntotal, "PairTableRx::compute::mixWtSite1");
+ memory->create(mixWtSite2, ntotal, "PairTableRx::compute::mixWtSite2");
+
+ for (int i = 0; i < ntotal; ++i)
+ getMixingWeights(i, mixWtSite1old[i], mixWtSite2old[i], mixWtSite1[i], mixWtSite2[i]);
+ }
+
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
@@ -130,10 +138,10 @@ void PairTableRX::compute(int eflag, int vflag)
double uCGnew_i = 0.0;
double fx_i = 0.0, fy_i = 0.0, fz_i = 0.0;
- fractionOld1_i = fractionOld1[i];
- fractionOld2_i = fractionOld2[i];
- fraction1_i = fraction1[i];
- fraction2_i = fraction2[i];
+ mixWtSite1old_i = mixWtSite1old[i];
+ mixWtSite2old_i = mixWtSite2old[i];
+ mixWtSite1_i = mixWtSite1[i];
+ mixWtSite2_i = mixWtSite2[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@@ -147,10 +155,10 @@ void PairTableRX::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
- fractionOld1_j = fractionOld1[j];
- fractionOld2_j = fractionOld2[j];
- fraction1_j = fraction1[j];
- fraction2_j = fraction2[j];
+ mixWtSite1old_j = mixWtSite1old[j];
+ mixWtSite2old_j = mixWtSite2old[j];
+ mixWtSite1_j = mixWtSite1[j];
+ mixWtSite2_j = mixWtSite2[j];
tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq)
@@ -186,8 +194,8 @@ void PairTableRX::compute(int eflag, int vflag)
value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value;
}
- if (isite1 == isite2) fpair = sqrt(fractionOld1_i*fractionOld2_j)*fpair;
- else fpair = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*fpair;
+ if (isite1 == isite2) fpair = sqrt(mixWtSite1old_i*mixWtSite2old_j)*fpair;
+ else fpair = (sqrt(mixWtSite1old_i*mixWtSite2old_j) + sqrt(mixWtSite2old_i*mixWtSite1old_j))*fpair;
fx_i += delx*fpair;
fy_i += dely*fpair;
@@ -208,11 +216,11 @@ void PairTableRX::compute(int eflag, int vflag)
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) *
tb->deltasq6;
if (isite1 == isite2){
- evdwlOld = sqrt(fractionOld1_i*fractionOld2_j)*evdwl;
- evdwl = sqrt(fraction1_i*fraction2_j)*evdwl;
+ evdwlOld = sqrt(mixWtSite1old_i*mixWtSite2old_j)*evdwl;
+ evdwl = sqrt(mixWtSite1_i*mixWtSite2_j)*evdwl;
} else {
- evdwlOld = (sqrt(fractionOld1_i*fractionOld2_j) + sqrt(fractionOld2_i*fractionOld1_j))*evdwl;
- evdwl = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*evdwl;
+ evdwlOld = (sqrt(mixWtSite1old_i*mixWtSite2old_j) + sqrt(mixWtSite2old_i*mixWtSite1old_j))*evdwl;
+ evdwl = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*evdwl;
}
evdwlOld *= factor_lj;
evdwl *= factor_lj;
@@ -238,10 +246,10 @@ void PairTableRX::compute(int eflag, int vflag)
}
if (vflag_fdotr) virial_fdotr_compute();
- memory->destroy(fractionOld1);
- memory->destroy(fractionOld2);
- memory->destroy(fraction1);
- memory->destroy(fraction2);
+ memory->destroy(mixWtSite1old);
+ memory->destroy(mixWtSite2old);
+ memory->destroy(mixWtSite1);
+ memory->destroy(mixWtSite2);
}
/* ----------------------------------------------------------------------
@@ -291,6 +299,8 @@ void PairTableRX::settings(int narg, char **arg)
else if (strcmp(arg[iarg],"msm") == 0) msmflag = 1;
else if (strcmp(arg[iarg],"dispersion") == 0) dispersionflag = 1;
else if (strcmp(arg[iarg],"tip4p") == 0) tip4pflag = 1;
+ else if (strcmp(arg[iarg],"fractional") == 0) fractionalWeighting = true;
+ else if (strcmp(arg[iarg],"molecular") == 0) fractionalWeighting = false;
else error->all(FLERR,"Illegal pair_style command");
iarg++;
}
@@ -1059,17 +1069,17 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
int tlm1 = tablength - 1;
Table *tb = &tables[tabindex[itype][jtype]];
- double fraction1_i, fraction1_j;
- double fraction2_i, fraction2_j;
- double fractionOld1_i, fractionOld1_j;
- double fractionOld2_i, fractionOld2_j;
+ double mixWtSite1_i, mixWtSite1_j;
+ double mixWtSite2_i, mixWtSite2_j;
+ double mixWtSite1old_i, mixWtSite1old_j;
+ double mixWtSite2old_i, mixWtSite2old_j;
fraction = 0.0;
a = 0.0;
b = 0.0;
- getParams(i,fractionOld1_i,fractionOld2_i,fraction1_i,fraction2_i);
- getParams(j,fractionOld1_j,fractionOld2_j,fraction1_j,fraction2_j);
+ getMixingWeights(i,mixWtSite1old_i,mixWtSite2old_i,mixWtSite1_i,mixWtSite2_i);
+ getMixingWeights(j,mixWtSite1old_j,mixWtSite2old_j,mixWtSite1_j,mixWtSite2_j);
if (rsq < tb->innersq) error->one(FLERR,"Pair distance < table inner cutoff");
@@ -1102,8 +1112,8 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
fforce = factor_lj * value;
}
- if (isite1 == isite2) fforce = sqrt(fraction1_i*fraction2_j)*fforce;
- else fforce = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*fforce;
+ if (isite1 == isite2) fforce = sqrt(mixWtSite1_i*mixWtSite2_j)*fforce;
+ else fforce = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*fforce;
if (tabstyle == LOOKUP)
phi = tb->e[itable];
@@ -1113,8 +1123,8 @@ double PairTableRX::single(int i, int j, int itype, int jtype, double rsq,
phi = a * tb->e[itable] + b * tb->e[itable+1] +
((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6;
- if (isite1 == isite2) phi = sqrt(fraction1_i*fraction2_j)*phi;
- else phi = (sqrt(fraction1_i*fraction2_j) + sqrt(fraction2_i*fraction1_j))*phi;
+ if (isite1 == isite2) phi = sqrt(mixWtSite1_i*mixWtSite2_j)*phi;
+ else phi = (sqrt(mixWtSite1_i*mixWtSite2_j) + sqrt(mixWtSite2_i*mixWtSite1_j))*phi;
return factor_lj*phi;
}
@@ -1141,46 +1151,74 @@ void *PairTableRX::extract(const char *str, int &dim)
/* ---------------------------------------------------------------------- */
-void PairTableRX::getParams(int id, double &fractionOld1, double &fractionOld2, double &fraction1, double &fraction2)
+void PairTableRX::getMixingWeights(int id, double &mixWtSite1old, double &mixWtSite2old, double &mixWtSite1, double &mixWtSite2)
{
- double nTotal = 0.0;
- double nTotalOld = 0.0;
+ double fractionOFAold, fractionOFA;
+ double fractionOld1, fraction1;
+ double fractionOld2, fraction2;
+ double nMoleculesOFAold, nMoleculesOFA;
+ double nMoleculesOld1, nMolecules1;
+ double nMoleculesOld2, nMolecules2;
+ double nTotal, nTotalOld;
+
+ nTotal = 0.0;
+ nTotalOld = 0.0;
for (int ispecies = 0; ispecies < nspecies; ++ispecies){
nTotal += atom->dvector[ispecies][id];
nTotalOld += atom->dvector[ispecies+nspecies][id];
}
- if(nTotal < 1e-8 || nTotalOld < 1e-8)
- error->all(FLERR,"The number of molecules in CG particle is less than 1e-8.");
+ if(nTotal < MY_EPSILON || nTotalOld < MY_EPSILON)
+ error->all(FLERR,"The number of molecules in CG particle is less than 10*DBL_EPSILON.");
if (isOneFluid(isite1) == false){
- fractionOld1 = atom->dvector[isite1+nspecies][id]/nTotalOld;
- fraction1 = atom->dvector[isite1][id]/nTotal;
+ nMoleculesOld1 = atom->dvector[isite1+nspecies][id];
+ nMolecules1 = atom->dvector[isite1][id];
+ fractionOld1 = nMoleculesOld1/nTotalOld;
+ fraction1 = nMolecules1/nTotal;
}
if (isOneFluid(isite2) == false){
- fractionOld2 = atom->dvector[isite2+nspecies][id]/nTotalOld;
- fraction2 = atom->dvector[isite2][id]/nTotal;
+ nMoleculesOld2 = atom->dvector[isite2+nspecies][id];
+ nMolecules2 = atom->dvector[isite2][id];
+ fractionOld2 = nMoleculesOld2/nTotalOld;
+ fraction2 = nMolecules2/nTotal;
}
if (isOneFluid(isite1) || isOneFluid(isite2)){
- double fractionOld = 0.0;
- double fraction = 0.0;
+ nMoleculesOFAold = 0.0;
+ nMoleculesOFA = 0.0;
+ fractionOFAold = 0.0;
+ fractionOFA = 0.0;
for (int ispecies = 0; ispecies < nspecies; ispecies++){
if (isite1 == ispecies || isite2 == ispecies) continue;
-
- fractionOld += atom->dvector[ispecies+nspecies][id]/nTotalOld;
- fraction += atom->dvector[ispecies][id]/nTotal;
+ nMoleculesOFAold += atom->dvector[ispecies+nspecies][id];
+ nMoleculesOFA += atom->dvector[ispecies][id];
+ fractionOFAold += atom->dvector[ispecies+nspecies][id]/nTotalOld;
+ fractionOFA += atom->dvector[ispecies][id]/nTotal;
}
-
if(isOneFluid(isite1)){
- fractionOld1 = fractionOld;
- fraction1 = fraction;
+ nMoleculesOld1 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules1 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld1 = fractionOFAold;
+ fraction1 = fractionOFA;
}
-
if(isOneFluid(isite2)){
- fractionOld2 = fractionOld;
- fraction2 = fraction;
+ nMoleculesOld2 = 1.0-(nTotalOld-nMoleculesOFAold);
+ nMolecules2 = 1.0-(nTotal-nMoleculesOFA);
+ fractionOld2 = fractionOFAold;
+ fraction2 = fractionOFA;
}
}
-}
+ if(fractionalWeighting){
+ mixWtSite1old = fractionOld1;
+ mixWtSite1 = fraction1;
+ mixWtSite2old = fractionOld2;
+ mixWtSite2 = fraction2;
+ } else {
+ mixWtSite1old = nMoleculesOld1;
+ mixWtSite1 = nMolecules1;
+ mixWtSite2old = nMoleculesOld2;
+ mixWtSite2 = nMolecules2;
+ }
+}
diff --git a/src/USER-DPD/pair_table_rx.h b/src/USER-DPD/pair_table_rx.h
index f04ebced207fba9ad3cb8bd1e4b041aa13c367f8..c6afe6a8d59eea99f021add838112be13250683d 100644
--- a/src/USER-DPD/pair_table_rx.h
+++ b/src/USER-DPD/pair_table_rx.h
@@ -72,7 +72,8 @@ class PairTableRX : public Pair {
int nspecies;
char *site1, *site2;
int isite1, isite2;
- void getParams(int, double &, double &, double &, double &);
+ void getMixingWeights(int, double &, double &, double &, double &);
+ bool fractionalWeighting;
};
@@ -163,7 +164,7 @@ When using pair style table with a long-range KSpace solver, the
cutoffs for all atom type pairs must all be the same, since the
long-range solver starts at that cutoff.
-E: The number of molecules in CG particle is less than 1e-8
+E: The number of molecules in CG particle is less than 10*DBL_EPSILON
Self-explanatory. Check the species concentrations have been properly set
and check the reaction kinetic solver parameters in fix rx to more for