ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2011-12-08 17:56:37 +00:00
parent 9692aa6aaf
commit 16a2afc642
15 changed files with 111 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
src/KSPACE/pppm.cpp
View File

@ -22,6 +22,7 @@
#include "stdlib.h"
#include "math.h"
#include "pppm.h"
#include "math_const.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
@ -35,8 +36,6 @@
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -1050,7 +1049,8 @@ double PPPM::rms(double h, double prd, bigint natoms,
compute difference in real-space and kspace RMS precision
------------------------------------------------------------------------- */
double PPPM::diffpr(double h_x, double h_y, double h_z, double q2, double **acons)
double PPPM::diffpr(double h_x, double h_y, double h_z, double q2,
double **acons)
{
double lprx,lpry,lprz,kspace_prec,real_prec;
double xprd = domain->xprd;
@ -1738,6 +1738,7 @@ void PPPM::fieldforce()
}
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;

15
src/KSPACE/pppm_tip4p.cpp
View File

@ -40,6 +40,18 @@ using namespace LAMMPS_NS;
PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp, int narg, char **arg) :
PPPM(lmp, narg, arg) {}
/* ---------------------------------------------------------------------- */
void PPPMTIP4P::init()
{
// TIP4P PPPM requires newton on, b/c it computes forces on ghost atoms
if (force->newton == 0)
error->all(FLERR,"Kspace style pppm/tip4p requires newton on");
PPPM::init();
}
/* ----------------------------------------------------------------------
find center grid pt for each of my particles
check that full stencil for the particle will fit in my 3d brick
@ -206,13 +218,14 @@ void PPPMTIP4P::fieldforce()
}
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
if (type[i] != typeO) {
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
f[i][2] += qfactor*ekz;
} else {
} else {
fx = qfactor * ekx;
fy = qfactor * eky;
fz = qfactor * ekz;

1
src/KSPACE/pppm_tip4p.h
View File

@ -28,6 +28,7 @@ class PPPMTIP4P : public PPPM {
public:
PPPMTIP4P(class LAMMPS *, int, char **);
virtual ~PPPMTIP4P () {};
void init();
protected:
virtual void particle_map();

1
src/MANYBODY/pair_airebo.h
View File

@ -163,6 +163,7 @@ class PairAIREBO : public Pair {
};
/* kronecker delta function returning a double */
inline double kronecker(const int a, const int b) const {
return (a == b) ? 1.0 : 0.0;
};

3
src/REPLICA/neb.cpp
View File

@ -215,6 +215,7 @@ void NEB::run()
// break induced if converged
// damped dynamic min styles insure all replicas converge together
timer->init();
timer->barrier_start(TIME_LOOP);
while (update->minimize->niter < n1steps) {
@ -284,7 +285,9 @@ void NEB::run()
// break induced if converged
// damped dynamic min styles insure all replicas converge together
timer->init();
timer->barrier_start(TIME_LOOP);
while (update->minimize->niter < n2steps) {
update->minimize->run(nevery);
print_status();

5
src/REPLICA/prd.cpp
View File

@ -261,8 +261,11 @@ void PRD::command(int narg, char **arg)
quench();
ncoincident = 0;
share_event(0,0);
timer->init();
timer->barrier_start(TIME_LOOP);
time_start = timer->array[TIME_LOOP];
log_event();
// do full init/setup since are starting all replicas after event
@ -343,6 +346,7 @@ void PRD::command(int narg, char **arg)
update->integrate->setup();
timer->barrier_start(TIME_LOOP);
if (t_corr > 0) replicate(ireplica);
if (temp_flag == 0) {
if (ireplica == universe->iworld)
@ -350,6 +354,7 @@ void PRD::command(int narg, char **arg)
MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[ireplica],
universe->uworld);
}
timer->barrier_stop(TIME_LOOP);
time_comm += timer->array[TIME_LOOP];

9
src/REPLICA/tad.cpp
View File

@ -250,8 +250,11 @@ void TAD::command(int narg, char **arg)
// This should work with if uncommented, but does not
// if (universe->iworld == 0) {
fix_event->store_state();
quench();
timer->init();
timer->barrier_start(TIME_LOOP);
time_start = timer->array[TIME_LOOP];
fix_event->store_event_tad(update->ntimestep);
@ -731,7 +734,7 @@ void TAD::perform_neb(int ievent)
memory->destroy(buf_init);
memory->destroy(buf_final);
// Run NEB
// run NEB
int beginstep_hold = update->beginstep;
int endstep_hold = update->endstep;
@ -763,7 +766,7 @@ void TAD::perform_neb(int ievent)
universe->uscreen = uscreen_lammps;
}
// Extract barrier energy from NEB
// extract barrier energy from NEB
if (universe->iworld == 0)
fix_event_list[ievent]->ebarrier = neb->ebf;
@ -786,7 +789,7 @@ void TAD::perform_neb(int ievent)
delete [] args;
// Clean up
// clean up
modify->delete_fix("neb");
delete neb;

7
src/REPLICA/temper.cpp
View File

@ -68,7 +68,8 @@ void Temper::command(int narg, char **arg)
error->all(FLERR,"Must have more than one processor partition to temper");
if (domain->box_exist == 0)
error->all(FLERR,"Temper command before simulation box is defined");
if (narg != 6 && narg != 7) error->universe_all(FLERR,"Illegal temper command");
if (narg != 6 && narg != 7)
error->universe_all(FLERR,"Illegal temper command");
int nsteps = atoi(arg[0]);
nevery = atoi(arg[1]);
@ -87,7 +88,8 @@ void Temper::command(int narg, char **arg)
// swap frequency must evenly divide total # of timesteps
if (nevery == 0) error->universe_all(FLERR,"Invalid frequency in temper command");
if (nevery == 0)
error->universe_all(FLERR,"Invalid frequency in temper command");
nswaps = nsteps/nevery;
if (nswaps*nevery != nsteps)
error->universe_all(FLERR,"Non integer # of swaps in temper command");
@ -203,6 +205,7 @@ void Temper::command(int narg, char **arg)
print_status();
}
timer->init();
timer->barrier_start(TIME_LOOP);
for (int iswap = 0; iswap < nswaps; iswap++) {

40
src/comm.cpp
View File

@ -50,6 +50,7 @@ using namespace LAMMPS_NS;
#define BIG 1.0e20
enum{SINGLE,MULTI};
enum{NONE,MULTIPLE};
/* ----------------------------------------------------------------------
setup MPI and allocate buffer space
@ -70,6 +71,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
cutghostuser = 0.0;
ghost_velocity = 0;
recv_from_partition = send_to_partition = -1;
other_partition_style = NONE;
// use of OpenMP threads
// query OpenMP for number of threads/process set by user at run-time
@ -159,11 +161,23 @@ void Comm::set_processors(int narg, char **arg)
int irecv = atoi(arg[iarg+2]);
if (isend < 1 || isend > universe->nworlds ||
irecv < 1 || irecv > universe->nworlds || isend == irecv)
error->all(FLERR,"Invalid partition in processors part command");
if (isend-1 == universe->iworld) send_to_partition = irecv-1;
if (irecv-1 == universe->iworld) recv_from_partition = isend-1;
error->all(FLERR,"Invalid partitions in processors part command");
if (isend-1 == universe->iworld) {
if (send_to_partition >= 0)
error->all(FLERR,
"Sending partition in processors part command "
"is already a sender");
send_to_partition = irecv-1;
}
if (irecv-1 == universe->iworld) {
if (recv_from_partition >= 0)
error->all(FLERR,
"Receiving partition in processors part command "
"is already a receiver");
recv_from_partition = isend-1;
}
if (strcmp(arg[iarg+3],"multiple") == 0) {
if (universe->iworld == irecv-1) other_partition_style = 0;
if (universe->iworld == irecv-1) other_partition_style = MULTIPLE;
} else error->all(FLERR,"Illegal processors command");
iarg += 4;
} else error->all(FLERR,"Illegal processors command");
@ -1229,8 +1243,9 @@ void Comm::procs2box()
while (ipx <= nprocs) {
valid = 1;
if (user_procgrid[0] && ipx != user_procgrid[0]) valid = 0;
if (other_procgrid[0]) {
if (other_partition_style == 0 && other_procgrid[0] % ipx) valid = 0;
if (other_partition_style) {
if (other_partition_style == MULTIPLE && other_procgrid[0] % ipx)
valid = 0;
}
if (nprocs % ipx) valid = 0;
if (!valid) {
@ -1242,8 +1257,9 @@ void Comm::procs2box()
while (ipy <= nprocs/ipx) {
valid = 1;
if (user_procgrid[1] && ipy != user_procgrid[1]) valid = 0;
if (other_procgrid[1]) {
if (other_partition_style == 0 && other_procgrid[1] % ipy) valid = 0;
if (other_partition_style) {
if (other_partition_style == MULTIPLE && other_procgrid[1] % ipy)
valid = 0;
}
if ((nprocs/ipx) % ipy) valid = 0;
if (!valid) {
@ -1254,8 +1270,9 @@ void Comm::procs2box()
ipz = nprocs/ipx/ipy;
valid = 1;
if (user_procgrid[2] && ipz != user_procgrid[2]) valid = 0;
if (other_procgrid[2]) {
if (other_partition_style == 0 && other_procgrid[2] % ipz) valid = 0;
if (other_partition_style) {
if (other_partition_style == MULTIPLE && other_procgrid[2] % ipz)
valid = 0;
}
if (domain->dimension == 2 && ipz != 1) valid = 0;
if (!valid) {
@ -1276,7 +1293,8 @@ void Comm::procs2box()
ipx++;
}
if (procgrid[0] == 0) error->one(FLERR,"Could not layout grid of processors");
if (procgrid[0] == 0)
error->all(FLERR,"Could not layout grid of processors",1);
}
/* ----------------------------------------------------------------------

11
src/error.cpp
View File

@ -27,6 +27,7 @@ Error::Error(LAMMPS *lmp) : Pointers(lmp) {}
/* ----------------------------------------------------------------------
called by all procs in universe
close all output, screen, and log files in world and universe
no abort, so insure all procs in universe call, else will hang
------------------------------------------------------------------------- */
void Error::universe_all(const char *file, int line, const char *str)
@ -53,6 +54,7 @@ void Error::universe_all(const char *file, int line, const char *str)
/* ----------------------------------------------------------------------
called by one proc in universe
forces abort of entire universe if any proc in universe calls
------------------------------------------------------------------------- */
void Error::universe_one(const char *file, int line, const char *str)
@ -66,9 +68,13 @@ void Error::universe_one(const char *file, int line, const char *str)
/* ----------------------------------------------------------------------
called by all procs in one world
close all output, screen, and log files in world
insure all procs in world call, else will hang
if abort = 0 (default):
if only one world in universe calls, universe will hang
if abort = 1: force abort of entire universe if any world in universe calls
------------------------------------------------------------------------- */
void Error::all(const char *file, int line, const char *str)
void Error::all(const char *file, int line, const char *str, int abort)
{
MPI_Barrier(world);
@ -84,6 +90,7 @@ void Error::all(const char *file, int line, const char *str)
if (screen && screen != stdout) fclose(screen);
if (logfile) fclose(logfile);
if (abort) MPI_Abort(world,1);
MPI_Finalize();
exit(1);
}
@ -92,6 +99,7 @@ void Error::all(const char *file, int line, const char *str)
called by one proc in world
write to world screen only if non-NULL on this proc
always write to universe screen
forces abort of entire world (and universe) if any proc in world calls
------------------------------------------------------------------------- */
void Error::one(const char *file, int line, const char *str)
@ -132,6 +140,7 @@ void Error::message(const char *file, int line, char *str, int logflag)
/* ----------------------------------------------------------------------
called by all procs in one world
close all output, screen, and log files in world
no abort, so insure all procs in world call, else will hang
------------------------------------------------------------------------- */
void Error::done()

2
src/error.h
View File

@ -25,7 +25,7 @@ class Error : protected Pointers {
void universe_all(const char *, int, const char *);
void universe_one(const char *, int, const char *);
void all(const char *, int, const char *);
void all(const char *, int, const char *, int = 0);
void one(const char *, int, const char *);
void warning(const char *, int, const char *, int = 1);
void message(const char *, int, char *, int = 1);

2
src/lammps.cpp
View File

@ -482,6 +482,7 @@ void LAMMPS::post_create()
/* ----------------------------------------------------------------------
initialize top-level classes
do not initialize Timer class, other classes like Run() do that explicitly
------------------------------------------------------------------------- */
void LAMMPS::init()
@ -500,7 +501,6 @@ void LAMMPS::init()
neighbor->init(); // neighbor must come after force, modify
comm->init(); // comm must come after force, modify, neighbor, atom
output->init(); // output must come after domain, force, modify
timer->init();
}
/* ----------------------------------------------------------------------

1
src/minimize.cpp
View File

@ -53,6 +53,7 @@ void Minimize::command(int narg, char **arg)
lmp->init();
update->minimize->setup();
timer->init();
timer->barrier_start(TIME_LOOP);
update->minimize->run(update->nsteps);
timer->barrier_stop(TIME_LOOP);

12
src/run.cpp
View File

@ -168,11 +168,9 @@ void Run::command(int narg, char **arg)
if (preflag || update->first_update == 0) {
lmp->init();
update->integrate->setup();
} else {
timer->init();
output->setup(0);
}
} else output->setup(0);
timer->init();
timer->barrier_start(TIME_LOOP);
update->integrate->run(nsteps);
timer->barrier_stop(TIME_LOOP);
@ -208,11 +206,9 @@ void Run::command(int narg, char **arg)
if (preflag || iter == 0) {
lmp->init();
update->integrate->setup();
} else {
timer->init();
output->setup(0);
}
} else output->setup(0);
timer->init();
timer->barrier_start(TIME_LOOP);
update->integrate->run(nsteps);
timer->barrier_stop(TIME_LOOP);

32
src/verlet_split.cpp
View File

@ -215,11 +215,14 @@ void VerletSplit::run(int n)
{
int nflag,ntimestep,sortflag;
// reset LOOP timer due to imbalance in setup() on 2 partitions
// setup initial Rspace <-> Kspace comm params
// sync both partitions before start timer
MPI_Barrier(universe->uworld);
timer->array[TIME_LOOP] = MPI_Wtime();
timer->init();
timer->barrier_start(TIME_LOOP);
// setup initial Rspace <-> Kspace comm params
rk_setup();
// flags for timestepping iterations
@ -315,6 +318,13 @@ void VerletSplit::run(int n)
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
// TIP4P PPPM puts forces on ghost atoms, so must reverse_comm()
if (tip4p_flag && force->newton) {
comm->reverse_comm();
timer->stamp(TIME_COMM);
}
}
// comm and sum Kspace forces back to Rspace procs
@ -384,8 +394,10 @@ void VerletSplit::rk_setup()
MPI_Gatherv(atom->q,n,MPI_DOUBLE,atom->q,qsize,qdisp,MPI_DOUBLE,0,block);
// for TIP4P also need to send type and tag
// KSpace proc needs to acquire ghost atoms and map all its atoms
// for TIP4P also need to send atom type and tag
// KSpace procs need to acquire ghost atoms and map all their atoms
// map_clear() call is in lieu of comm->exchange() which performs map_clear
// borders() call acquires ghost atoms and maps them
if (tip4p_flag) {
MPI_Gatherv(atom->type,n,MPI_INT,atom->type,qsize,qdisp,MPI_INT,0,block);
@ -394,7 +406,7 @@ void VerletSplit::rk_setup()
if (triclinic) domain->x2lamda(atom->nlocal);
if (domain->box_change) comm->setup();
timer->stamp();
atom->map_clear(); // in lieu of comm->exchange()
atom->map_clear();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
@ -442,6 +454,14 @@ void VerletSplit::r2k_comm()
force->kspace->setup();
}
}
// for TIP4P, Kspace partition needs to update its ghost atoms
if (tip4p_flag && !master) {
timer->stamp();
comm->forward_comm();
timer->stamp(TIME_COMM);
}
}
/* ----------------------------------------------------------------------