From a329a724b2a5bf013905028ed1fc158f16541156 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 1 Jul 2016 23:31:52 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15253 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-DPD/atom_vec_dpd.cpp | 47 --------------- src/USER-DPD/atom_vec_dpd.h | 2 - src/USER-DPD/fix_rx.cpp | 80 +++++++++++++------------ src/USER-DPD/fix_shardlow.cpp | 106 +++++++++++++++++++++++++++++++++- src/USER-DPD/fix_shardlow.h | 15 ++++- src/USER-DPD/pair_dpd_fdt.cpp | 2 - src/USER-DPD/pair_exp6_rx.cpp | 8 ++- 7 files changed, 167 insertions(+), 93 deletions(-) diff --git a/src/USER-DPD/atom_vec_dpd.cpp b/src/USER-DPD/atom_vec_dpd.cpp index 586ddca02b..2bd460475d 100644 --- a/src/USER-DPD/atom_vec_dpd.cpp +++ b/src/USER-DPD/atom_vec_dpd.cpp @@ -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; } diff --git a/src/USER-DPD/atom_vec_dpd.h b/src/USER-DPD/atom_vec_dpd.h index 077d18485f..234d2ccce7 100644 --- a/src/USER-DPD/atom_vec_dpd.h +++ b/src/USER-DPD/atom_vec_dpd.h @@ -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 }; } diff --git a/src/USER-DPD/fix_rx.cpp b/src/USER-DPD/fix_rx.cpp index 33d6c1c988..5639646156 100644 --- a/src/USER-DPD/fix_rx.cpp +++ b/src/USER-DPD/fix_rx.cpp @@ -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 // 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 -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;jjme == 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]; } diff --git a/src/USER-DPD/fix_shardlow.cpp b/src/USER-DPD/fix_shardlow.cpp index 7248df2e67..1b58c2f0bd 100644 --- a/src/USER-DPD/fix_shardlow.cpp +++ b/src/USER-DPD/fix_shardlow.cpp @@ -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; +} + diff --git a/src/USER-DPD/fix_shardlow.h b/src/USER-DPD/fix_shardlow.h index 6421737b3d..4afce632ce 100644 --- a/src/USER-DPD/fix_shardlow.h +++ b/src/USER-DPD/fix_shardlow.h @@ -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 + }; } diff --git a/src/USER-DPD/pair_dpd_fdt.cpp b/src/USER-DPD/pair_dpd_fdt.cpp index 21c5e420fa..4ad4582f0f 100644 --- a/src/USER-DPD/pair_dpd_fdt.cpp +++ b/src/USER-DPD/pair_dpd_fdt.cpp @@ -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 diff --git a/src/USER-DPD/pair_exp6_rx.cpp b/src/USER-DPD/pair_exp6_rx.cpp index db3885b762..3e3e1b4e27 100644 --- a/src/USER-DPD/pair_exp6_rx.cpp +++ b/src/USER-DPD/pair_exp6_rx.cpp @@ -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]);