git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15253 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2016-07-01 23:31:52 +00:00
parent 2a1b47172d
commit a329a724b2
7 changed files with 167 additions and 93 deletions

View File

@ -77,7 +77,6 @@ void AtomVecDPD::grow(int n)
uCG = memory->grow(atom->uCG,nmax,"atom:uCG");
uCGnew = memory->grow(atom->uCGnew,nmax,"atom:uCGnew");
duChem = memory->grow(atom->duChem,nmax,"atom:duChem");
ssaAIR = memory->grow(atom->ssaAIR,nmax,"atom:ssaAIR");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -101,7 +100,6 @@ void AtomVecDPD::grow_reset()
uCG = atom->uCG;
uCGnew = atom->uCGnew;
duChem = atom->duChem;
ssaAIR = atom->ssaAIR;
}
/* ----------------------------------------------------------------------
@ -126,7 +124,6 @@ void AtomVecDPD::copy(int i, int j, int delflag)
uChem[j] = uChem[i];
uCG[j] = uCG[i];
uCGnew[j] = uCGnew[i];
ssaAIR[j] = ssaAIR[i];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -516,41 +513,6 @@ int AtomVecDPD::pack_border_hybrid(int n, int *list, double *buf)
return m;
}
/* ----------------------------------------------------------------------
convert atom coords into the ssa active interaction region number
------------------------------------------------------------------------- */
int AtomVecDPD::coord2ssaAIR(double *x)
{
int ix, iy, iz;
ix = iy = iz = 0;
if (x[2] < domain->sublo[2]) iz = -1;
if (x[2] > domain->subhi[2]) iz = 1;
if (x[1] < domain->sublo[1]) iy = -1;
if (x[1] > domain->subhi[1]) iy = 1;
if (x[0] < domain->sublo[0]) ix = -1;
if (x[0] > domain->subhi[0]) ix = 1;
if(iz < 0){
return 0;
} else if(iz == 0){
if( iy<0 ) return 0; // bottom left/middle/right
if( (iy==0) && (ix<0) ) return 0; // left atoms
if( (iy==0) && (ix==0) ) return 1; // Locally owned atoms
if( (iy==0) && (ix>0) ) return 3; // Right atoms
if( (iy>0) && (ix==0) ) return 2; // Top-middle atoms
if( (iy>0) && (ix!=0) ) return 4; // Top-right and top-left atoms
} else { // iz > 0
if((ix==0) && (iy==0)) return 5; // Back atoms
if((ix==0) && (iy!=0)) return 6; // Top-back and bottom-back atoms
if((ix!=0) && (iy==0)) return 7; // Left-back and right-back atoms
if((ix!=0) && (iy!=0)) return 8; // Back corner atoms
}
return 0;
}
/* ---------------------------------------------------------------------- */
void AtomVecDPD::unpack_border(int n, int first, double *buf)
@ -573,7 +535,6 @@ void AtomVecDPD::unpack_border(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
ssaAIR[i] = coord2ssaAIR(x[i]);
}
if (atom->nextra_border)
@ -607,7 +568,6 @@ void AtomVecDPD::unpack_border_vel(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
ssaAIR[i] = coord2ssaAIR(x[i]);
}
if (atom->nextra_border)
@ -648,7 +608,6 @@ int AtomVecDPD::unpack_border_hybrid(int n, int first, double *buf)
uChem[i] = buf[m++];
uCG[i] = buf[m++];
uCGnew[i] = buf[m++];
ssaAIR[i] = coord2ssaAIR(x[i]);
}
return m;
}
@ -710,7 +669,6 @@ int AtomVecDPD::unpack_exchange(double *buf)
uChem[nlocal] = buf[m++];
uCG[nlocal] = buf[m++];
uCGnew[nlocal] = buf[m++];
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
@ -801,7 +759,6 @@ int AtomVecDPD::unpack_restart(double *buf)
uCond[nlocal] = buf[m++];
uMech[nlocal] = buf[m++];
uChem[nlocal] = buf[m++];
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
double **extra = atom->extra;
if (atom->nextra_store) {
@ -842,7 +799,6 @@ void AtomVecDPD::create_atom(int itype, double *coord)
uCG[nlocal] = 0.0;
uCGnew[nlocal] = 0.0;
duChem[nlocal] = 0.0;
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
atom->nlocal++;
}
@ -883,7 +839,6 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
uChem[nlocal] = 0.0;
uCG[nlocal] = 0.0;
uCGnew[nlocal] = 0.0;
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
atom->nlocal++;
}
@ -896,7 +851,6 @@ void AtomVecDPD::data_atom(double *coord, tagint imagetmp, char **values)
int AtomVecDPD::data_atom_hybrid(int nlocal, char **values)
{
dpdTheta[nlocal] = atof(values[0]);
ssaAIR[nlocal] = 1; /* coord2ssaAIR(x[nlocal]) */
return 1;
}
@ -978,7 +932,6 @@ bigint AtomVecDPD::memory_usage()
if (atom->memcheck("uCG")) bytes += memory->usage(uCG,nmax);
if (atom->memcheck("uCGnew")) bytes += memory->usage(uCGnew,nmax);
if (atom->memcheck("duChem")) bytes += memory->usage(duChem,nmax);
if (atom->memcheck("ssaAIR")) bytes += memory->usage(ssaAIR,nmax);
return bytes;
}

View File

@ -60,7 +60,6 @@ class AtomVecDPD : public AtomVec {
bigint memory_usage();
double *uCond,*uMech,*uChem,*uCG,*uCGnew,*rho,*dpdTheta;
double *duChem;
int *ssaAIR; // Shardlow Splitting Algorithm Active Interaction Region number
protected:
tagint *tag;
@ -68,7 +67,6 @@ class AtomVecDPD : public AtomVec {
imageint *image;
double **x,**v,**f;
int coord2ssaAIR(double *); // map atom coord to an AIR number
};
}

View File

@ -27,6 +27,7 @@
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_special.h"
#include "pair_dpd_fdt_energy.h"
#include <float.h> // DBL_EPSILON
@ -36,6 +37,7 @@
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathSpecial;
enum{NONE,HARMONIC};
enum{LUCY};
@ -53,18 +55,6 @@ typedef double TimerType;
TimerType getTimeStamp(void) { return MPI_Wtime(); }
double getElapsedTime( const TimerType &t0, const TimerType &t1) { return t1-t0; }
// Fast (non-IEEE) x^p function where x is a double and p is a positive, integral value.
template <typename intType>
inline double fastpowi( const double x, const intType p )
{
if (p == 1) return x;
else if (p == 2) return x*x;
else if (p == 3) return x*x*x;
else if (p == 0) return 1.0;
else
return pow(x, (double)p);
}
} // end namespace
/* ---------------------------------------------------------------------- */
@ -80,6 +70,8 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
params = NULL;
mol2param = NULL;
pairDPDE = NULL;
id_fix_species = NULL;
id_fix_species_old = NULL;
const int Verbosity = 1;
@ -130,7 +122,7 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
else
msg += std::string("dense");
error->message(__FILE__, __LINE__, msg.c_str());
error->message(FLERR, msg.c_str());
}
}
@ -172,7 +164,7 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
if (comm->me == 0 and Verbosity > 1){
char msg[128];
sprintf(msg, "FixRX: RK4 numSteps= %d", minSteps);
error->message(__FILE__, __LINE__, msg);
error->message(FLERR, msg);
}
}
else if (odeIntegrationFlag == ODE_LAMMPS_RK4 && narg>8){
@ -198,7 +190,7 @@ FixRX::FixRX(LAMMPS *lmp, int narg, char **arg) :
//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(__FILE__, __LINE__, msg);
error->message(FLERR, msg);
}
}
@ -222,11 +214,17 @@ FixRX::~FixRX()
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 );
@ -244,10 +242,10 @@ void FixRX::post_constructor()
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");
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];
char **tmpspecies = new char*[maxspecies];
for(int jj=0; jj < maxspecies; jj++)
tmpspecies[jj] = NULL;
@ -256,7 +254,7 @@ void FixRX::post_constructor()
FILE *fp;
fp = NULL;
if (comm->me == 0) {
fp = fopen(kineticsFile,"r");
fp = force->open_potential(kineticsFile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",kineticsFile);
@ -266,7 +264,7 @@ void FixRX::post_constructor()
// Assign species names to tmpspecies array and determine the number of unique species
int n,nwords,ispecies;
int n,nwords;
char line[MAXLINE],*ptr;
int eof = 0;
char * word;
@ -306,7 +304,7 @@ void FixRX::post_constructor()
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)];
tmpspecies[nUniqueSpecies] = new char[strlen(word)+1];
strcpy(tmpspecies[nUniqueSpecies],word);
nUniqueSpecies++;
}
@ -350,8 +348,8 @@ void FixRX::post_constructor()
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)];
newarg2[ii+3] = new char[strlen(str2)];
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);
}
@ -368,11 +366,15 @@ void FixRX::post_constructor()
if(nspecies==0) error->all(FLERR,"There are no rx species specified.");
for(int jj=0;jj<maxspecies;jj++) delete tmpspecies[jj];
delete tmpspecies;
for(int jj=0;jj<nspecies;jj++) {
delete[] tmpspecies[jj];
delete[] newarg[jj+3];
delete[] newarg2[jj+3];
}
delete newarg;
delete newarg2;
delete[] newarg;
delete[] newarg2;
delete[] tmpspecies;
read_file( kineticsFile );
@ -479,7 +481,7 @@ int FixRX::initSparse()
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(__FILE__, __LINE__, msg);
error->message(FLERR, msg);
}
// Allocate the sparse matrix data.
@ -757,7 +759,7 @@ void FixRX::pre_force(int vflag)
if (nFails > 0){
char sbuf[128];
sprintf(sbuf,"in FixRX::pre_force, ODE solver failed for %d atoms.", nFails);
error->warning(__FILE__, __LINE__, sbuf);
error->warning(FLERR, sbuf);
}
// Compute and report ODE diagnostics, if requested.
@ -791,7 +793,7 @@ void FixRX::read_file(char *file)
FILE *fp;
fp = NULL;
if (comm->me == 0) {
fp = fopen(file,"r");
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open rx file %s",file);
@ -828,7 +830,7 @@ void FixRX::read_file(char *file)
}
// open file on proc 0
if (comm->me == 0) fp = fopen(file,"r");
if (comm->me == 0) fp = force->open_potential(file);
// read each reaction from kinetics file
eof=0;
@ -895,8 +897,10 @@ void FixRX::read_file(char *file)
}
}
if(ispecies==nspecies){
printf("%s mol fraction is not found in data file\n",word);
printf("nspecies=%d ispecies=%d\n",nspecies,ispecies);
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");
@ -1612,12 +1616,12 @@ int FixRX::rhs_sparse(double t, const double *y, double *dydt, void *v_params) c
{
double rxnRateLawForward;
if (isIntegral(i)){
rxnRateLawForward = kFor[i] * ::fastpowi( conc[ nuk[i][0] ], inu[i][0]);
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 *= ::fastpowi( conc[k], inu[i][kk] );
rxnRateLawForward *= powint( conc[k], inu[i][kk] );
}
} else {
rxnRateLawForward = kFor[i] * pow( conc[ nuk[i][0] ], nu[i][0]);
@ -1692,7 +1696,7 @@ void FixRX::computeLocalTemperature()
int newton_pair = force->newton_pair;
// local temperature variables
double wij, fr, fr4;
double wij=0.0, fr, fr4;
double *dpdTheta = atom->dpdTheta;
// Initialize the local density and local temperature arrays
@ -1742,7 +1746,7 @@ void FixRX::computeLocalTemperature()
if(wtFlag==LUCY){
wij = (1.0+3.0*ratio) * (1.0-ratio)*(1.0-ratio)*(1.0-ratio);
dpdThetaLocal[i] += wij/dpdTheta[j];
if (newton_pair || j < nlocal)
if (newton_pair || j < nlocal)
dpdThetaLocal[j] += wij/dpdTheta[i];
}

View File

@ -105,8 +105,29 @@ FixShardlow::FixShardlow(LAMMPS *lmp, int narg, char **arg) :
if(pairDPD == NULL && pairDPDE == NULL)
error->all(FLERR,"Must use pair_style dpd/fdt or dpd/fdt/energy with fix shardlow");
if(!(atom->dpd_flag))
error->all(FLERR,"Must use atom_style dpd with fix shardlow");
// Setup the ssaAIR array
atom->ssaAIR = NULL;
grow_arrays(atom->nmax);
for (int i = 0; i < atom->nlocal; i++) {
atom->ssaAIR[i] = 1; /* coord2ssaAIR(x[i]) */
}
// Setup callbacks for maintaining atom->ssaAIR[]
atom->add_callback(0); // grow (aka exchange)
atom->add_callback(1); // restart
atom->add_callback(2); // border
}
/* ---------------------------------------------------------------------- */
FixShardlow::~FixShardlow()
{
atom->delete_callback(id, 0);
atom->delete_callback(id, 1);
atom->delete_callback(id, 2);
memory->destroy(atom->ssaAIR);
}
/* ---------------------------------------------------------------------- */
@ -506,3 +527,84 @@ void FixShardlow::unpack_reverse_comm(int n, int *list, double *buf)
}
}
}
/* ----------------------------------------------------------------------
convert atom coords into the ssa active interaction region number
------------------------------------------------------------------------- */
int FixShardlow::coord2ssaAIR(double *x)
{
int ix, iy, iz;
ix = iy = iz = 0;
if (x[2] < domain->sublo[2]) iz = -1;
if (x[2] >= domain->subhi[2]) iz = 1;
if (x[1] < domain->sublo[1]) iy = -1;
if (x[1] >= domain->subhi[1]) iy = 1;
if (x[0] < domain->sublo[0]) ix = -1;
if (x[0] >= domain->subhi[0]) ix = 1;
if(iz < 0){
return 0;
} else if(iz == 0){
if( iy<0 ) return 0; // bottom left/middle/right
if( (iy==0) && (ix<0) ) return 0; // left atoms
if( (iy==0) && (ix==0) ) return 1; // Locally owned atoms
if( (iy==0) && (ix>0) ) return 3; // Right atoms
if( (iy>0) && (ix==0) ) return 2; // Top-middle atoms
if( (iy>0) && (ix!=0) ) return 4; // Top-right and top-left atoms
} else { // iz > 0
if((ix==0) && (iy==0)) return 5; // Back atoms
if((ix==0) && (iy!=0)) return 6; // Top-back and bottom-back atoms
if((ix!=0) && (iy==0)) return 7; // Left-back and right-back atoms
if((ix!=0) && (iy!=0)) return 8; // Back corner atoms
}
return 0;
}
/* ---------------------------------------------------------------------- */
void FixShardlow::grow_arrays(int nmax)
{
memory->grow(atom->ssaAIR,nmax,"fix_shardlow:ssaAIR");
}
void FixShardlow::copy_arrays(int i, int j, int delflag)
{
atom->ssaAIR[j] = atom->ssaAIR[i];
}
void FixShardlow::set_arrays(int i)
{
atom->ssaAIR[i] = 1; /* coord2ssaAIR(x[i]) */
}
int FixShardlow::unpack_border(int n, int first, double *buf)
{
int i,last = first + n;
for (i = first; i < last; i++) {
atom->ssaAIR[i] = coord2ssaAIR(atom->x[i]);
}
return 0;
}
int FixShardlow::unpack_exchange(int i, double *buf)
{
atom->ssaAIR[i] = 1; /* coord2ssaAIR(x[i]) */
return 0;
}
void FixShardlow::unpack_restart(int i, int nth)
{
atom->ssaAIR[i] = 1; /* coord2ssaAIR(x[i]) */
}
double FixShardlow::memory_usage()
{
double bytes = 0.0;
bytes += memory->usage(atom->ssaAIR,atom->nmax);
bytes += sizeof(double)*3*atom->nghost; // v_t0[]
return bytes;
}

View File

@ -27,12 +27,23 @@ namespace LAMMPS_NS {
class FixShardlow : public Fix {
public:
FixShardlow(class LAMMPS *, int, char **);
virtual ~FixShardlow() {}
~FixShardlow();
int setmask();
virtual void init_list(int,class NeighList *);
virtual void setup(int);
virtual void initial_integrate(int);
void grow_arrays(int);
void copy_arrays(int, int, int);
void set_arrays(int);
// int pack_border(int, int *, double *);
int unpack_border(int, int, double *);
int unpack_exchange(int, double *);
void unpack_restart(int, int);
double memory_usage();
protected:
int pack_reverse_comm(int, int, double *);
void unpack_reverse_comm(int, int *, double *);
@ -46,6 +57,8 @@ class FixShardlow : public Fix {
private:
class NeighList *list;
int coord2ssaAIR(double *); // map atom coord to an AIR number
};
}

View File

@ -179,8 +179,6 @@ void PairDPDfdt::settings(int narg, char **arg)
temperature = force->numeric(FLERR,arg[0]);
cut_global = force->numeric(FLERR,arg[1]);
seed = force->inumeric(FLERR,arg[2]);
if (atom->dpd_flag != 1)
error->all(FLERR,"pair_style dpd/fdt requires atom_style dpd");
// initialize Marsaglia RNG with processor-unique seed

View File

@ -75,13 +75,16 @@ PairExp6rx::PairExp6rx(LAMMPS *lmp) : Pair(lmp)
PairExp6rx::~PairExp6rx()
{
for (int i=0; i < nparams; ++i) {
delete[] params[i].name;
delete[] params[i].potential;
}
memory->destroy(params);
memory->destroy(mol2param);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
}
}
@ -631,6 +634,9 @@ void PairExp6rx::coeff(int narg, char **arg)
//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]);