Commit 3faa57a4 authored by Tim Mattox's avatar Tim Mattox
Browse files

USER-DPD: Several updates to *_rx files: 

1) Added MY_EPSILON to handle machine precision checks
2) Removed error checks for DPD-RX; enabled use with DPD-E simulations
3) Expanded the EOS functional form to allow corrections
   in the thermo file or on the command line
4) Updated naming convention from fraction to mixWtSite*
5) Changed the name of getParams() method to getMixingWeights()
6) getMixingWeights() now handles fractional and molecular weighting
7) Added optional argument (fractional or molecular) to pair_style command
8) Added argument to specify the exp6 parameter scaling method
   NOTE: Requires additional arguments in the pair coefficients,
   thus command line areguments are NOT backward-compatible.
parent fa435fb5

src/USER-DPD/fix_eos_table_rx.cpp

+99 −40
Original line number Diff line number Diff line
@@ -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.");



    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,6 +577,7 @@ 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");



  if(rx_flag){

    while (word) {

      for (ispecies = 0; ispecies < nspecies; ispecies++)

        if (strcmp(word,&atom->dname[ispecies][0]) == 0){
@@ -566,6 +599,11 @@ void FixEOStableRX::param_extract(Table *tb, char *line) 
      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;

src/USER-DPD/fix_eos_table_rx.h

+2 −1
Original line number Diff line number Diff line
@@ -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;

  };

}

src/USER-DPD/pair_exp6_rx.cpp

+300 −149

File changed.

Preview size limit exceeded, changes collapsed.

src/USER-DPD/pair_exp6_rx.h

+10 −3
Original line number Diff line number Diff line
@@ -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

src/USER-DPD/pair_multi_lucy_rx.cpp

+86 −40
Original line number Diff line number Diff line
@@ -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;

  }

}



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

CRA Git | Maintained and supported by SUSTech CRA and CCSE