/*---------------------------------------------------------------------- pmd.h is an include file for a parallel MD program, pmd.c. ----------------------------------------------------------------------*/ #include #include #include "mpi.h" /* Constants------------------------------------------------------------ NMAX = Manimum # of atoms per processor NEMAX = Maximum # of augmented (= resident + copied) atoms NDBUF = Size of a double-precision buffer, dbuf > 6*(# of boundary atoms for each neighbor) NBMAX = Maximum # of copied boundary atoms per neighbor. NCLMAX = Maximum # of cells per processor. RCUT = Potential cut-off length MOVED_OUT: Signifies a moved-out resident atom in function atom_move. EMPTY: Signifies the end of a linked list. ----------------------------------------------------------------------*/ #define NMAX 100000 #define NEMAX 200000 #define NDBUF 300000 #define NBMAX 100000 #define NCLMAX 100000 #define RCUT 2.5 #define MOVED_OUT -1.0e10 #define EMPTY -1 /* Constants for the random number generator */ #define D2P31M 2147483647.0 #define DMUL 16807.0 /* Variables------------------------------------------------------------ al[0|1|2] = Box length per processor in the x|y|z direction. n = # of resident atoms in this processor. nb = # of copied boundary atoms from neighbors. nglob = Total # of atoms summed over processors. r[NEMAX][3]: r[i][0|1|2] is the x|y|z coordinate of atom i (including the copied atoms). rv[NEMAX][3]: rv[i][0|1|2] is the x|y|z velocity of atom i (including the copied atoms). ra[NMAX][3]: ra[i][0|1|2] is the x|y|z acceleration on atom i. dbuf[NDBUF]: Buffer for sending double-precision data dbufr[NDBUF]: receiving vproc[0|1|2] = # of processors in the x|y|z direction. nproc = # of processors = vproc[0]*vproc[1]*vproc[2]. sid = Sequential processor ID. vid[3] = Vector processor ID; sid = vid[0]*vproc[1]*vproc[2] + vid[1]*vproc[2] + vid[2]. NN[6]: NN[ku] is the node ID of the neighbor specified by a neighbor. index, ku. The neighbor index is defined as: ku = 0: xlow (West ); 1: xhigh (East ); 2: ylow (South); 3: yhigh (North); 4: zlow (Down ); 5: zhigh (Up ). sv[6][3]: sv[ku][] is the shift vector to the ku-th neighbor. myparity[0|1|2] = Parity of vector processor ID in the x|y|z direction. lsb[6][NBMAX]: lsb[ku][0] is the total # of boundary atoms to be sent to neighbor ku; lsb[ku][k] is the atom ID, used in r, of the k-th atom to be sent. status: Returned by MPI message-passing routines. cpu: Elapsed wall-clock time in seconds. comt: Communication time in seconds. lc[3]: lc[0|1|2] is the # of cells in the x|y|z direction. rc[3]: rc[0|1|2] is the length of a cell in the x|y|z direction. lscl[NEMAX]: Linked cell lists. head[NCLMAX]: Headers for the linked cell lists. kinEnergy = Kinetic energy. potEnergy = Potential energy. totEnergy = Total energy. temperature = Current temperature. stepCount = Current time step. ----------------------------------------------------------------------*/ double al[3]; int n,nb,nglob; double r[NEMAX][3],rv[NEMAX][3],ra[NMAX][3]; double dbuf[NDBUF],dbufr[NDBUF]; int vproc[3] = {1,1,2}, nproc = 2; int sid,vid[3],nn[6],myparity[3]; double sv[6][3]; int lsb[6][NBMAX]; MPI_Status status; double cpu,comt; int head[NCLMAX],lscl[NEMAX],lc[3]; double rc[3]; double kinEnergy,potEnergy,totEnergy,temperature; int stepCount; double DeltaTH; /* Half the time step */ double Uc, Duc; /* Potential cut-off parameters */ /* Input data----------------------------------------------------------- Control data: pmd.in. ----------------------------------------------------------------------*/ int InitUcell[3]; /* Number of unit cells per processor */ double Density; /* Number density of atoms (in reduced unit) */ double InitTemp; /* Starting temperature (in reduced unit) */ double DeltaT; /* Size of a time step (in reduced unit) */ int StepLimit; /* Number of time steps to be simulated */ int StepAvg; /* Reporting interval for statistical data */ /* Functions & function prototypes------------------------------------*/ double SignR(double v,double x) {if (x > 0) return v; else return -v;} double Dmod(double a, double b) { int n; n = (int) (a/b); return (a - b*n); } double RandR(double *seed) { *seed = Dmod(*seed*DMUL,D2P31M); return (*seed/D2P31M); } void RandVec3(double *p, double *seed) { double x,y,s = 2.0; while (s > 1.0) { x = 2.0*RandR(seed) - 1.0; y = 2.0*RandR(seed) - 1.0; s = x*x + y*y; } p[2] = 1.0 - 2.0*s; s = 2.0*sqrt(1.0 - s); p[0] = s*x; p[1] = s*y; } void init_params(); void set_topology(); void init_conf(); void single_step(); void half_kick(); void atom_copy(); void compute_accel(); void eval_props(); void atom_move(); int bbd(double* ri, int ku); int bmv(double* ri, int ku); /*--------------------------------------------------------------------*/