From 16bb8a1439fba129139061be0e7ebaf313c21bb5 Mon Sep 17 00:00:00 2001 From: Mingjian Wen Date: Wed, 17 Apr 2019 16:58:18 -0500 Subject: [PATCH] Clean up comments --- src/USER-MISC/pair_drip.cpp | 174 ++++++++++++------------------------ src/USER-MISC/pair_drip.h | 82 +++++++++-------- 2 files changed, 101 insertions(+), 155 deletions(-) diff --git a/src/USER-MISC/pair_drip.cpp b/src/USER-MISC/pair_drip.cpp index 39c48eb56a..11721f048e 100644 --- a/src/USER-MISC/pair_drip.cpp +++ b/src/USER-MISC/pair_drip.cpp @@ -41,49 +41,41 @@ using namespace LAMMPS_NS; #define MAXLINE 1024 #define DELTA 4 -#define PGDELTA 1 #define HALF 0.5 /* ---------------------------------------------------------------------- */ PairDRIP::PairDRIP(LAMMPS *lmp) : Pair(lmp) { - // initialize element to parameter maps single_enable = 0; - nelements = 0; - elements = NULL; - nparams = maxparam = 0; + restartinfo = 0; + params = NULL; + nearest3neigh = NULL; + elements = NULL; elem2param = NULL; map = NULL; + nelements = 0; cutmax = 0.0; - -//j nmax = 0; -//j maxlocal = 0; -//j ipage = NULL; -//j pgsize = oneatom = 0; - - nearest3neigh = NULL; } /* ---------------------------------------------------------------------- */ PairDRIP::~PairDRIP() { -// delete [] ipage; - if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); + delete [] map; } - if (elements) + if (elements != NULL) { for (int i = 0; i < nelements; i++) delete [] elements[i]; - delete [] elements; + delete [] elements; + elements = NULL; + } memory->destroy(params); memory->destroy(elem2param); - if (allocated) delete [] map; - memory->destroy(nearest3neigh); } @@ -99,7 +91,6 @@ void PairDRIP::init_style() error->all(FLERR,"Pair style drip requires atom attribute molecule"); // need a full neighbor list, including neighbors of ghosts - int irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->full = 1; @@ -181,7 +172,6 @@ void PairDRIP::coeff(int narg, char **arg) } } - read_file(arg[2]); int count = 0; @@ -216,8 +206,8 @@ void PairDRIP::read_file(char *filename) int params_per_line = 14; char **words = new char*[params_per_line+1]; memory->sfree(params); - params = NULL; - nparams = maxparam = 0; + int nparams = 0; + int maxparam = 0; // open file on proc 0 @@ -326,9 +316,7 @@ void PairDRIP::read_file(char *filename) // set max cutoff if(params[nparams].rcut > cutmax) cutmax = params[nparams].rcut; - nparams++; - //if(nparams >= pow(atom->ntypes,3)) break; } memory->destroy(elem2param); @@ -353,18 +341,15 @@ void PairDRIP::read_file(char *filename) void PairDRIP::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype,k,l,kk,ll; tagint itag,jtag; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,fpair1,fpair2, r, rsq; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,r,rsq; int *ilist,*jlist,*numneigh,**firstneigh; double ni[DIM]; double dni_dri[DIM][DIM], dni_drnb1[DIM][DIM]; double dni_drnb2[DIM][DIM], dni_drnb3[DIM][DIM]; - - evdwl = 0.0; ev_init(eflag,vflag); double **x = atom->x; @@ -381,7 +366,6 @@ void PairDRIP::compute(int eflag, int vflag) find_nearest3neigh(); - // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; itag = tag[i]; @@ -392,8 +376,7 @@ void PairDRIP::compute(int eflag, int vflag) jlist = firstneigh[i]; jnum = numneigh[i]; - - // normal and its derivatives w.r.t. atom i and its 3 nearest neighs + // normal and its derivatives w.r.t. atom i and its 3 nearest neighbors calc_normal(i, ni, dni_dri,dni_drnb1, dni_drnb2, dni_drnb3); double fi[DIM] = {0., 0., 0.}; @@ -412,7 +395,6 @@ void PairDRIP::compute(int eflag, int vflag) Param& p = params[iparam_ij]; double rcutsq = p.rcutsq; - // only include the interation between different layers if (rsq < rcutsq && atom->molecule[i] != atom->molecule[j]) { @@ -425,13 +407,14 @@ void PairDRIP::compute(int eflag, int vflag) ni, dni_dri, dni_drnb1, dni_drnb2, dni_drnb3, fi, fj); if (eflag) evdwl = HALF * (phi_repul + phi_attr); + else evdwl = 0.0; if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,0,0,0,0); f[j][0] += fj[0]; f[j][1] += fj[1]; f[j][2] += fj[2]; - // *2 since v_tally has a 0.5 coeff + // multiply 2 since v_tally has a 0.5 coeff fj[0] *= 2; fj[1] *= 2; fj[2] *= 2; if (vflag_atom) v_tally(j, fj, x[j]); @@ -442,7 +425,7 @@ void PairDRIP::compute(int eflag, int vflag) f[i][1] += fi[1]; f[i][2] += fi[2]; - // *2 since v_tally has a 0.5 coeff + // multiply 2 since v_tally has a 0.5 coeff fi[0] *= 2; fi[1] *= 2; fi[2] *= 2; if (vflag_atom) v_tally(i, fi, x[i]); @@ -452,45 +435,15 @@ void PairDRIP::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); - - -printf("@@@ evflags in compute\n"); -printf("@@@ eflag_either=%d\n", eflag_either); -printf("@@@ eflag_global=%d\n", eflag_global); -printf("@@@ eflag_atom=%d\n", eflag_atom); -printf("@@@ vflag_either=%d\n", vflag_either); -printf("@@@ vflag_global=%d\n", vflag_global); -printf("@@@ vflag_atom=%d\n", vflag_atom); -printf("@@@ vflag_fdotr=%d\n", vflag_fdotr); - - -printf("@@@@@@@@@@@@@@@@@@@@@@@ virial\n"); -printf("%f, %f, %f, %f, %f, %f\n", virial[0], virial[1], virial[2], virial[3], virial[4], virial[5]); -printf("@@@@@@@@@@@@@@@@@@@@@@@ virial fdotr\n"); -virial[0]= virial[1]= virial[2]= virial[3]= virial[4]= virial[5]=0.; -virial_fdotr_compute(); -printf("%f, %f, %f, %f, %f, %f\n", virial[0], virial[1], virial[2], virial[3], virial[4], virial[5]); - -printf("@@@@@@@@@@@@@@@@@@@@@@@ virial from atom virial\n"); - - int allnum = list->inum + list->gnum; - double v[6] = {0., 0., 0., 0., 0., 0.}; - for (int kk=0; kkf; double **x = atom->x; - // params double C0 = p.C0; double C2 = p.C2; double C4 = p.C4; @@ -544,8 +498,6 @@ double PairDRIP::calc_repulsive(int const i, int const j, int nbj2 = nearest3neigh[j][1]; int nbj3 = nearest3neigh[j][2]; - double r = sqrt(rsq); - double fnbi1[DIM]; double fnbi2[DIM]; double fnbi3[DIM]; @@ -566,9 +518,9 @@ double PairDRIP::calc_repulsive(int const i, int const j, V3 drhosqij_drnb2; V3 drhosqij_drnb3; + double r = sqrt(rsq); - // derivative of rhosq w.r.t coordinates of atoms i, j, and the nearests 3 - // neighs of i + // derivative of rhosq w.r.t. atoms i j and the nearests 3 neighs of i get_drhosqij(rvec, ni, dni_dri, dni_drnb1, dni_drnb2, dni_drnb3, drhosqij_dri, drhosqij_drj, drhosqij_drnb1, drhosqij_drnb2, drhosqij_drnb3); @@ -588,10 +540,11 @@ double PairDRIP::calc_repulsive(int const i, int const j, double dtp; double tp = tap(r, cutoff, dtp); - /* exponential part */ + // exponential part double V1 = exp(-lambda * (r - z0)); double dV1 = -V1 * lambda; + // total energy double phi = tp * V1 * V2; for (int k = 0; k < DIM; k++) { @@ -601,16 +554,15 @@ double PairDRIP::calc_repulsive(int const i, int const j, fi[k] += tmp; fj[k] -= tmp; - // the following incldue the transverse decay part tdij and the dihedral part gij + // contributions from the transverse decay part tdij and the dihedral part gij + // derivative of V2 contribute to atoms i, j fi[k] -= HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_dri[k] + dgij_dri[k]); fj[k] -= HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drj[k] + dgij_drj[k]); - // derivative of V2 contribute to nearest 3 neighs of atom i fnbi1[k] =- HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drnb1[k] + dgij_drk1[k]); fnbi2[k] =- HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drnb2[k] + dgij_drk2[k]); fnbi3[k] =- HALF * tp * V1 * ((dtdij + dgij_drhosq) * drhosqij_drnb3[k] + dgij_drk3[k]); - // derivative of V2 contribute to nearest 3 neighs of atom j fnbj1[k] =- HALF * tp * V1 * dgij_drl1[k]; fnbj2[k] =- HALF * tp * V1 * dgij_drl2[k]; @@ -627,8 +579,7 @@ double PairDRIP::calc_repulsive(int const i, int const j, } if (vflag_atom) { - - // *2 since v_tally has a 0.5 coeff + // multiply since v_tally has a 0.5 coeff for (int k = 0; k < DIM; k++) { fnbi1[k]*=2; fnbi2[k]*=2; @@ -649,16 +600,13 @@ double PairDRIP::calc_repulsive(int const i, int const j, } - /* ---------------------------------------------------------------------- */ void PairDRIP::find_nearest3neigh() { - int i,j,ii,jj,n,allnum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *ilist,*jlist,*numneigh,**firstneigh; - //int *neighptr; double **x = atom->x; int *type = atom->type; @@ -683,12 +631,11 @@ void PairDRIP::find_nearest3neigh() jlist = firstneigh[i]; jnum = numneigh[i]; - // init nb1 to be the 1st nearest neigh, nb3 the 3rd nearest int nb1 = -1; int nb2 = -1; int nb3 = -1; - double nb1_rsq = 1.1e10; + double nb1_rsq = 1.0e10 + 1; double nb2_rsq = 2.0e10; double nb3_rsq = 3.0e10; @@ -741,13 +688,11 @@ void PairDRIP::find_nearest3neigh() } // loop over ii } - /* ---------------------------------------------------------------------- */ void PairDRIP::calc_normal(int const i, double * const normal, V3 *const dn_dri, V3 *const dn_drk1, V3 *const dn_drk2, V3 *const dn_drk3) { - int k1 = nearest3neigh[i][0]; int k2 = nearest3neigh[i][1]; int k3 = nearest3neigh[i][2]; @@ -764,9 +709,9 @@ void PairDRIP::calc_normal(int const i, double * const normal, deriv_cross(x[k1], x[k2], x[k3], normal, dn_drk1, dn_drk2, dn_drk3); } - /* ---------------------------------------------------------------------- */ -void PairDRIP::get_drhosqij( double const* rij, double const* ni, + +void PairDRIP::get_drhosqij(double const* rij, double const* ni, V3 const* dni_dri, V3 const* dni_drn1, V3 const* dni_drn2, V3 const* dni_drn3, double* const drhosq_dri, double* const drhosq_drj, @@ -795,12 +740,10 @@ void PairDRIP::get_drhosqij( double const* rij, double const* ni, } } +/* ---------------------------------------------------------------------- + derivartive of transverse decay function f(rho) w.r.t. rho +------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- */ - - -// derivartive of transverse decay function f(rho) w.r.t rho double PairDRIP::td(double C0, double C2, double C4, double delta, double const* const rvec, double r, const double* const n, @@ -810,7 +753,8 @@ double PairDRIP::td(double C0, double C2, double C4, double delta, rho_sq = r * r - n_dot_r * n_dot_r; - if (rho_sq < 0) { // in case n is [0, 0, 1] and rho_sq is negative due to numerical error + // in case n is [0, 0, 1] and rho_sq is negative due to numerical error + if (rho_sq < 0) { rho_sq = 0; } @@ -822,9 +766,10 @@ double PairDRIP::td(double C0, double C2, double C4, double delta, return td; } +/* ---------------------------------------------------------------------- + derivartive of dihedral angle func gij w.r.t rho, and atom positions +------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- */ -// derivartive of dihedral angle func gij w.r.t rho, and atom positions double PairDRIP::dihedral( const int i, const int j, Param& p, double const rhosq, double& d_drhosq, double* const d_dri, double* const d_drj, @@ -935,9 +880,10 @@ double PairDRIP::dihedral( return dihe; } +/* ---------------------------------------------------------------------- + compute cos(omega_kijl) and the derivateives +------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- */ -// compute cos(omega_kijl) and the derivateives double PairDRIP::deriv_cos_omega( double const* rk, double const* ri, double const* rj, double const* rl, double* const dcos_drk, double* const dcos_dri, double* const dcos_drj, double* const dcos_drl) @@ -954,10 +900,12 @@ double PairDRIP::deriv_cos_omega( double const* rk, double const* ri, double deijl_drl[DIM][DIM]; - // ejik and derivatives (Note the dejik_dri ... returned are actually the transpose) + // ejik and derivatives + // Note the returned dejik_dri ... are actually the transpose deriv_cross(ri, rj, rk, ejik, dejik_dri, dejik_drj, dejik_drk); - // flip sign because deriv_cross computes rij cross rik, here we need rji cross rik + // flip sign + // deriv_cross computes rij cross rik, here we need rji cross rik for (int m = 0; m < DIM; m++) { ejik[m] = -ejik[m]; for (int n = 0; n < DIM; n++) { @@ -1002,11 +950,8 @@ double PairDRIP::deriv_cos_omega( double const* rk, double const* ri, return cos_omega; } - - - /* ---------------------------------------------------------------------- */ -// tap cutoff function + double PairDRIP::tap(double r, double cutoff, double& dtap) { double t; @@ -1027,9 +972,8 @@ double PairDRIP::tap(double r, double cutoff, double& dtap) return t; } - /* ---------------------------------------------------------------------- */ -// tap rho + double PairDRIP::tap_rho(double rhosq, double cut_rhosq, double& drhosq) { double roc_sq; @@ -1047,11 +991,12 @@ double PairDRIP::tap_rho(double rhosq, double cut_rhosq, double& drhosq) } -/* ---------------------------------------------------------------------- */ -// Compute the normalized cross product of two vector rkl, rkm, and the -// derivates w.r.t rk, rl, rm. -// NOTE, the dcross_drk, dcross_drl, and dcross_drm is actually the transpose -// of the actual one. +/* ---------------------------------------------------------------------- + Compute the normalized cross product of two vector rkl, rkm, and the + derivates w.r.t rk, rl, rm. + NOTE, the returned dcross_drk, dcross_drl, and dcross_drm are actually the + transpose. +------------------------------------------------------------------------- */ void PairDRIP::deriv_cross( double const* rk, double const* rl, double const* rm, double* const cross, V3 *const dcross_drk, @@ -1134,6 +1079,3 @@ void PairDRIP::deriv_cross( double const* rk, double const* rl, double const* rm } } - - - diff --git a/src/USER-MISC/pair_drip.h b/src/USER-MISC/pair_drip.h index 13f3391355..c4a130d226 100644 --- a/src/USER-MISC/pair_drip.h +++ b/src/USER-MISC/pair_drip.h @@ -11,6 +11,17 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing author: Mingjian Wen (University of Minnesota) + e-mail: wenxx151@umn.edu + based on "pair_style kolmogorov/crespi/full" by Wengen Ouyang + + This implements the DRIP model as described in + M. Wen, S. Carr, S. Fang, E. Kaxiras, and E. B. Tadmor, Phys. Rev. B, 98, + 235404 (2018). +------------------------------------------------------------------------- */ + + #ifdef PAIR_CLASS PairStyle(drip,PairDRIP) @@ -46,70 +57,67 @@ class PairDRIP : public Pair { int ielement,jelement; double C0,C2,C4,C,delta,lambda,A,z0,B,eta,rhocut,rcut; double rhocutsq, rcutsq; - double delta2inv,z06; }; Param *params; // parameter set for I-J interactions + int ** nearest3neigh; // nearest 3 neighbors of atoms char **elements; // names of unique elements int **elem2param; // mapping from element pairs to parameters int *map; // mapping from atom types to elements int nelements; // # of unique elements - int nparams; // # of stored parameter sets - int maxparam; // max # of parameter sets double cutmax; // max cutoff for all species - int ** nearest3neigh; // nearest 3 neighbors of atoms void read_file(char * ); void allocate(); // DRIP specific functions - double calc_attractive(int const i, int const j, Param& p, - double const rsq, double const * rvec, double * const fi, double * const fj); + double calc_attractive(int const, int const, Param&, double const, + double const *, double * const, double * const); - double calc_repulsive(int const i, int const j, - Param& p, double const rsq, double const * rvec, - double const * ni, V3 const * dni_dri, - V3 const * dni_drnb1, V3 const * dni_drnb2, V3 const * dni_drnb3, - double * const fi, double * const fj); + double calc_repulsive(int const , int const , + Param& , double const , double const * , + double const * , V3 const * , + V3 const * , V3 const * , V3 const * , + double * const , double * const ); void find_nearest3neigh(); - void calc_normal(int const i, double * const normal, - V3 *const dn_dri, V3 *const dn_drk1, V3 *const dn_drk2, V3 *const dn_drk3); + void calc_normal(int const , double * const , + V3 *const , V3 *const , V3 *const , V3 *const ); -void get_drhosqij( double const* rij, double const* ni, - V3 const* dni_dri, V3 const* dni_drn1, - V3 const* dni_drn2, V3 const* dni_drn3, - double* const drhosq_dri, double* const drhosq_drj, - double* const drhosq_drn1, double* const drhosq_drn2, - double* const drhosq_drn3); +void get_drhosqij( double const* , double const* , + V3 const* , V3 const* , + V3 const* , V3 const* , + double* const , double* const , + double* const , double* const , + double* const ); - double td(double C0, double C2, double C4, double delta, - double const* const rvec, double r, - const double* const n, - double& rho_sq, double& dtd); + double td(double , double , double , double , + double const* const , double , + const double* const , + double& , double& ); double dihedral( - const int i, const int j, Param& p, double const rhosq, double& d_drhosq, - double* const d_dri, double* const d_drj, - double* const d_drk1, double* const d_drk2, double* const d_drk3, - double* const d_drl1, double* const d_drl2, double* const d_drl3); + const int , const int , Param& , double const , double& , + double* const , double* const , + double* const , double* const , double* const , + double* const , double* const , double* const ); - double deriv_cos_omega( double const* rk, double const* ri, - double const* rj, double const* rl, double* const dcos_drk, - double* const dcos_dri, double* const dcos_drj, double* const dcos_drl); + double deriv_cos_omega( double const* , double const* , + double const* , double const* , double* const , + double* const , double* const , double* const ); - double tap(double r, double cutoff, double& dtap); + double tap(double , double , double& ); - double tap_rho(double rhosq, double cut_rhosq, double& drhosq); + double tap_rho(double , double , double& ); -void deriv_cross( double const* rk, double const* rl, double const* rm, - double* const cross, V3 *const dcross_drk, - V3 *const dcross_drl, V3 *const dcross_drm); +void deriv_cross( double const* , double const* , double const* , + double* const , V3 *const , + V3 *const , V3 *const ); // inline functions inline double dot(double const* x, double const* y) { @@ -123,7 +131,6 @@ void deriv_cross( double const* rk, double const* rl, double const* rm, } } - }; } @@ -153,7 +160,4 @@ E: No enough neighbors to construct normal Cannot find three neighbors within cutoff of the target atom. Check the configuration. - - - */