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

This commit is contained in:
sjplimp
2013-07-23 15:10:07 +00:00
parent 366484eb34
commit b49a5519f2
11 changed files with 22 additions and 340 deletions

View File

@ -177,7 +177,7 @@ void PairLJSDKCoulLongGPU::init_style()
// setup force tables
if (ncoultablebits) init_tables();
if (ncoultablebits) init_tables(cut_coul,NULL);
int maxspecial=0;
if (atom->molecular)

View File

@ -219,12 +219,10 @@ void TAD::command(int narg, char **arg)
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,
"Step CPU N M Status Barrier Margin t_lo delt_lo\n"
);
"Step CPU N M Status Barrier Margin t_lo delt_lo\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,
"Step CPU N M Status Barrier Margin t_lo delt_lo\n"
);
"Step CPU N M Status Barrier Margin t_lo delt_lo\n");
}
ulogfile_lammps = universe->ulogfile;
@ -285,14 +283,9 @@ void TAD::command(int narg, char **arg)
while (update->ntimestep < update->endstep) {
dynamics();
fix_event->store_state();
quench();
event_flag = check_event();
MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);

View File

@ -65,26 +65,6 @@ void PairCGCMMCoulLong::allocate()
memory->create(cut_coulsq,n+1,n+1,"paircg:cut_coulsq");
}
/* ----------------------------------------------------------------------
free memory for tables used in pair computations
------------------------------------------------------------------------- */
void PairCGCMMCoulLong::free_tables()
{
memory->destroy(rtable);
memory->destroy(drtable);
memory->destroy(ftable);
memory->destroy(dftable);
memory->destroy(ctable);
memory->destroy(dctable);
memory->destroy(etable);
memory->destroy(detable);
memory->destroy(vtable);
memory->destroy(dvtable);
memory->destroy(ptable);
memory->destroy(dptable);
}
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::init_style()
@ -109,7 +89,7 @@ void PairCGCMMCoulLong::init_style()
// setup force tables
if (ncoultablebits) init_tables();
if (ncoultablebits) init_tables(cut_coul_global,cut_respa);
}
/* ---------------------------------------------------------------------- */
@ -128,171 +108,6 @@ double PairCGCMMCoulLong::init_one(int i, int j)
/* ---------------------------------------------------------------------- */
void PairCGCMMCoulLong::init_tables()
{
int masklo,maskhi;
double r,grij,expm2,derfc,rsw;
double qqrd2e = force->qqrd2e;
tabinnersq = tabinner*tabinner;
init_bitmap(tabinner,cut_coul_global,ncoultablebits,
masklo,maskhi,ncoulmask,ncoulshiftbits);
int ntable = 1;
for (int i = 0; i < ncoultablebits; i++) ntable *= 2;
// linear lookup tables of length N = 2^ncoultablebits
// stored value = value at lower edge of bin
// d values = delta from lower edge to upper edge of bin
if (ftable) free_tables();
memory->create(rtable,ntable,"pair:rtable");
memory->create(ftable,ntable,"pair:ftable");
memory->create(ctable,ntable,"pair:ctable");
memory->create(etable,ntable,"pair:etable");
memory->create(drtable,ntable,"pair:drtable");
memory->create(dftable,ntable,"pair:dftable");
memory->create(dctable,ntable,"pair:dctable");
memory->create(detable,ntable,"pair:detable");
if (cut_respa == NULL) {
vtable = ptable = dvtable = dptable = NULL;
} else {
memory->create(vtable,ntable,"pair:vtable");
memory->create(ptable,ntable,"pair:ptable");
memory->create(dvtable,ntable,"pair:dvtable");
memory->create(dptable,ntable,"pair:dptable");
}
union_int_float_t rsq_lookup;
union_int_float_t minrsq_lookup;
int itablemin;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
} else {
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
}
}
}
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
for (int i = 0; i < ntablem1; i++) {
drtable[i] = 1.0/(rtable[i+1] - rtable[i]);
dftable[i] = ftable[i+1] - ftable[i];
dctable[i] = ctable[i+1] - ctable[i];
detable[i] = etable[i+1] - etable[i];
}
if (cut_respa) {
for (int i = 0; i < ntablem1; i++) {
dvtable[i] = vtable[i+1] - vtable[i];
dptable[i] = ptable[i+1] - ptable[i];
}
}
// get the delta values for the last table entries
// tables are connected periodically between 0 and ntablem1
drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]);
dftable[ntablem1] = ftable[0] - ftable[ntablem1];
dctable[ntablem1] = ctable[0] - ctable[ntablem1];
detable[ntablem1] = etable[0] - etable[ntablem1];
if (cut_respa) {
dvtable[ntablem1] = vtable[0] - vtable[ntablem1];
dptable[ntablem1] = ptable[0] - ptable[ntablem1];
}
// get the correct delta values at itablemax
// smallest r is in bin itablemin
// largest r is in bin itablemax, which is itablemin-1,
// or ntablem1 if itablemin=0
// deltas at itablemax only needed if corresponding rsq < cut*cut
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq_lookup.f < cut_coulsq_global) {
rsq_lookup.f = cut_coulsq_global;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
c_tmp = qqrd2e/r;
e_tmp = qqrd2e/r * derfc;
} else {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
c_tmp = 0.0;
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
} else {
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
c_tmp = qqrd2e/r;
}
}
}
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
if (cut_respa) {
dvtable[itablemax] = v_tmp - vtable[itablemax];
dptable[itablemax] = p_tmp - ptable[itablemax];
}
}
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- *
* the real compute work is done in the PairCMMCommon::eval_XXX<>() templates
* in the common PairCG class. Through using templates we can have one

View File

@ -47,8 +47,6 @@ namespace LAMMPS_NS {
protected:
void allocate();
void init_tables();
void free_tables();
};
}

View File

@ -1,4 +1,4 @@
/* ----------------------------------------------------------------------
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@ -78,12 +78,6 @@ namespace LAMMPS_NS {
double cut_coul_global, cut_coulsq_global, kappa, g_ewald;
double **cut_coul, **cut_coulsq;
// tables
double tabinnersq;
double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable;
double *etable,*detable,*ptable,*dptable,*vtable,*dvtable;
int ncoulshiftbits,ncoulmask;
// r-RESPA parameters
double *cut_respa;

View File

@ -371,9 +371,9 @@ void PairLJSDKCoulLong::init_style()
error->all(FLERR,"Pair style requires a KSpace style");
g_ewald = force->kspace->g_ewald;
// setup force tables
// setup force tables (no rRESPA support yet)
if (ncoultablebits) init_tables();
if (ncoultablebits) init_tables(cut_coul,NULL);
}
/* ----------------------------------------------------------------------
@ -438,115 +438,6 @@ double PairLJSDKCoulLong::init_one(int i, int j)
return cut;
}
/* ----------------------------------------------------------------------
setup force tables used in compute routines
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::init_tables()
{
int masklo,maskhi;
double r,grij,expm2,derfc;
double qqrd2e = force->qqrd2e;
tabinnersq = tabinner*tabinner;
init_bitmap(tabinner,cut_coul,ncoultablebits,
masklo,maskhi,ncoulmask,ncoulshiftbits);
int ntable = 1;
for (int i = 0; i < ncoultablebits; i++) ntable *= 2;
// linear lookup tables of length N = 2^ncoultablebits
// stored value = value at lower edge of bin
// d values = delta from lower edge to upper edge of bin
if (ftable) free_tables();
memory->create(rtable,ntable,"pair:rtable");
memory->create(ftable,ntable,"pair:ftable");
memory->create(ctable,ntable,"pair:ctable");
memory->create(etable,ntable,"pair:etable");
memory->create(drtable,ntable,"pair:drtable");
memory->create(dftable,ntable,"pair:dftable");
memory->create(dctable,ntable,"pair:dctable");
memory->create(detable,ntable,"pair:detable");
union_int_float_t rsq_lookup;
union_int_float_t minrsq_lookup;
int itablemin;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
for (int i = 0; i < ntablem1; i++) {
drtable[i] = 1.0/(rtable[i+1] - rtable[i]);
dftable[i] = ftable[i+1] - ftable[i];
dctable[i] = ctable[i+1] - ctable[i];
detable[i] = etable[i+1] - etable[i];
}
// get the delta values for the last table entries
// tables are connected periodically between 0 and ntablem1
drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]);
dftable[ntablem1] = ftable[0] - ftable[ntablem1];
dctable[ntablem1] = ctable[0] - ctable[ntablem1];
detable[ntablem1] = etable[0] - etable[ntablem1];
// get the correct delta values at itablemax
// smallest r is in bin itablemin
// largest r is in bin itablemax, which is itablemin-1,
// or ntablem1 if itablemin=0
// deltas at itablemax only needed if corresponding rsq < cut*cut
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
c_tmp = qqrd2e/r;
e_tmp = qqrd2e/r * derfc;
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
@ -659,22 +550,6 @@ void PairLJSDKCoulLong::write_data_all(FILE *fp)
epsilon[i][j],sigma[i][j],cut_lj[i][j]);
}
/* ----------------------------------------------------------------------
free memory for tables used in pair computations
------------------------------------------------------------------------- */
void PairLJSDKCoulLong::free_tables()
{
memory->destroy(rtable);
memory->destroy(drtable);
memory->destroy(ftable);
memory->destroy(dftable);
memory->destroy(ctable);
memory->destroy(dctable);
memory->destroy(etable);
memory->destroy(detable);
}
/* ---------------------------------------------------------------------- */
double PairLJSDKCoulLong::single(int i, int j, int itype, int jtype,

View File

@ -64,14 +64,7 @@ class PairLJSDKCoulLong : public Pair {
double cut_lj_global;
double g_ewald;
double tabinnersq;
double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable;
double *etable,*detable;
int ncoulshiftbits,ncoulmask;
void allocate();
void init_tables();
void free_tables();
private:
template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval();

View File

@ -64,6 +64,18 @@ void Error::universe_one(const char *file, int line, const char *str)
MPI_Abort(universe->uworld,1);
}
/* ----------------------------------------------------------------------
called by one proc in universe
prints a warning message to the screen
------------------------------------------------------------------------- */
void Error::universe_warn(const char *file, int line, const char *str)
{
if (universe->uscreen)
fprintf(universe->uscreen,"WARNING on proc %d: %s (%s:%d)\n",
universe->me,str,file,line);
}
/* ----------------------------------------------------------------------
called by all procs in one world
close all output, screen, and log files in world

View File

@ -24,6 +24,7 @@ class Error : protected Pointers {
void universe_all(const char *, int, const char *);
void universe_one(const char *, int, const char *);
void universe_warn(const char *, int, const char *);
void all(const char *, int, const char *);
void one(const char *, int, const char *);

View File

@ -204,7 +204,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
if (logflag == 0) {
universe->ulogfile = fopen("log.lammps","w");
if (universe->ulogfile == NULL)
error->universe_one(FLERR,"Cannot open log.lammps");
error->universe_warn(FLERR,"Cannot open log.lammps for writing");
} else if (strcmp(arg[logflag],"none") == 0)
universe->ulogfile = NULL;
else {

View File

@ -85,6 +85,7 @@ Output::Output(LAMMPS *lmp) : Pointers(lmp)
dump = NULL;
restart_flag = restart_flag_single = restart_flag_double = 0;
restart_every_single = restart_every_double = 0;
last_restart = -1;
restart1 = restart2a = restart2b = NULL;
var_restart_single = var_restart_double = NULL;