diff --git a/src/USER-MISC/pair_drip.cpp b/src/USER-MISC/pair_drip.cpp index 235766f43e..4d0a941a67 100644 --- a/src/USER-MISC/pair_drip.cpp +++ b/src/USER-MISC/pair_drip.cpp @@ -12,12 +12,13 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Wengen Ouyang (Tel Aviv University) - e-mail: w.g.ouyang at gmail dot com - based on previous versions by Jaap Kroes + Contributing author: Mingjian Wen (University of Minnesota) + e-mail: wenxx151@umn.edu + based on "pair_style kolmogorov/crespi/full" by Wengen Ouyang - This is a complete version of the potential described in - [Kolmogorov & Crespi, Phys. Rev. B 71, 235415 (2005)] + 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). ------------------------------------------------------------------------- */ #include @@ -47,18 +48,19 @@ using namespace LAMMPS_NS; PairDRIP::PairDRIP(LAMMPS *lmp) : Pair(lmp) { // initialize element to parameter maps + single_enable = 0; nelements = 0; elements = NULL; nparams = maxparam = 0; params = NULL; elem2param = NULL; - cutKCsq = NULL; + cutDRIPsq = NULL; map = NULL; nmax = 0; maxlocal = 0; - KC_numneigh = NULL; - KC_firstneigh = NULL; + DRIP_numneigh = NULL; + DRIP_firstneigh = NULL; ipage = NULL; pgsize = oneatom = 0; @@ -78,8 +80,8 @@ PairDRIP::PairDRIP(LAMMPS *lmp) : Pair(lmp) PairDRIP::~PairDRIP() { - memory->destroy(KC_numneigh); - memory->sfree(KC_firstneigh); + memory->destroy(DRIP_numneigh); + memory->sfree(DRIP_firstneigh); delete [] ipage; memory->destroy(normal); memory->destroy(dnormal); @@ -97,7 +99,7 @@ PairDRIP::~PairDRIP() delete [] elements; memory->destroy(params); memory->destroy(elem2param); - memory->destroy(cutKCsq); + memory->destroy(cutDRIPsq); if (allocated) delete [] map; } @@ -112,7 +114,7 @@ void PairDRIP::compute(int eflag, int vflag) double rsq,r,rhosq1,rhosq2,exp0,exp1,exp2,r2inv,r6inv,r8inv,Tap,dTap,Vkc; double frho1,frho2,sumC1,sumC2,sumC11,sumC22,sumCff,fsum,rdsq1,rdsq2; int *ilist,*jlist,*numneigh,**firstneigh; - int *KC_neighs_i,*KC_neighs_j; + int *DRIP_neighs_i,*DRIP_neighs_j; evdwl = 0.0; ev_init(eflag,vflag); @@ -139,7 +141,7 @@ void PairDRIP::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; // Build full neighbor list - KC_neigh(); + DRIP_neigh(); // Calculate the normals calc_normal(); @@ -255,9 +257,9 @@ void PairDRIP::compute(int eflag, int vflag) f[j][2] -= fkcz + fprod2[2]*Tap; // calculate the forces acted on the neighbors of atom i from atom j - KC_neighs_i = KC_firstneigh[i]; - for (kk = 0; kk < KC_numneigh[i]; kk++) { - k = KC_neighs_i[kk]; + DRIP_neighs_i = DRIP_firstneigh[i]; + for (kk = 0; kk < DRIP_numneigh[i]; kk++) { + k = DRIP_neighs_i[kk]; if (k == i) continue; // derivatives of the product of rij and ni respect to rk, k=0,1,2, where atom k is the neighbors of atom i dprodnorm1[0] = dnormal[0][0][kk][i]*delx + dnormal[1][0][kk][i]*dely + dnormal[2][0][kk][i]*delz; @@ -276,9 +278,9 @@ void PairDRIP::compute(int eflag, int vflag) } // calculate the forces acted on the neighbors of atom j from atom i - KC_neighs_j = KC_firstneigh[j]; - for (ll = 0; ll < KC_numneigh[j]; ll++) { - l = KC_neighs_j[ll]; + DRIP_neighs_j = DRIP_firstneigh[j]; + for (ll = 0; ll < DRIP_numneigh[j]; ll++) { + l = DRIP_neighs_j[ll]; if (l == j) continue; // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz; @@ -367,8 +369,8 @@ void PairDRIP::calc_normal() } cont = 0; - jlist = KC_firstneigh[i]; - jnum = KC_numneigh[i]; + jlist = DRIP_firstneigh[i]; + jnum = DRIP_numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; @@ -488,7 +490,6 @@ void PairDRIP::calc_normal() } } } -//############################################################################################## else if(cont == 3) { // for the atoms at the edge who has only two neighbor atoms @@ -636,9 +637,8 @@ void PairDRIP::calc_normal() else { error->one(FLERR,"There are too many neighbors for calculating normals"); } - -//############################################################################################## } + } /* ---------------------------------------------------------------------- @@ -659,7 +659,7 @@ void PairDRIP::init_style() neighbor->requests[irequest]->full = 1; neighbor->requests[irequest]->ghost = 1; - // local KC neighbor list + // local DRIP neighbor list // create pages if first time or if neighbor pgsize/oneatom has changed int create = 0; @@ -684,7 +684,7 @@ void PairDRIP::init_style() create neighbor list from main neighbor list for calculating the normals ------------------------------------------------------------------------- */ -void PairDRIP::KC_neigh() +void PairDRIP::DRIP_neigh() { int i,j,ii,jj,n,allnum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; @@ -696,10 +696,10 @@ void PairDRIP::KC_neigh() if (atom->nmax > maxlocal) { maxlocal = atom->nmax; - memory->destroy(KC_numneigh); - memory->sfree(KC_firstneigh); - memory->create(KC_numneigh,maxlocal,"DRIP:numneigh"); - KC_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *), + memory->destroy(DRIP_numneigh); + memory->sfree(DRIP_firstneigh); + memory->create(DRIP_numneigh,maxlocal,"DRIP:numneigh"); + DRIP_firstneigh = (int **) memory->smalloc(maxlocal*sizeof(int *), "DRIP:firstneigh"); } @@ -708,7 +708,7 @@ void PairDRIP::KC_neigh() numneigh = list->numneigh; firstneigh = list->firstneigh; - // store all KC neighs of owned and ghost atoms + // store all DRIP neighs of owned and ghost atoms // scan full neighbor list of I ipage->reset(); @@ -735,13 +735,13 @@ void PairDRIP::KC_neigh() delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; - if (rsq != 0 && rsq < cutKCsq[itype][jtype] && atom->molecule[i] == atom->molecule[j]) { + if (rsq != 0 && rsq < cutDRIPsq[itype][jtype] && atom->molecule[i] == atom->molecule[j]) { neighptr[n++] = j; } } - KC_firstneigh[i] = neighptr; - KC_numneigh[i] = n; + DRIP_firstneigh[i] = neighptr; + DRIP_numneigh[i] = n; if (n > 3) error->one(FLERR,"There are too many neighbors for some atoms, please check your configuration"); ipage->vgot(n); if (ipage->status()) @@ -875,7 +875,7 @@ double PairDRIP::init_one(int i, int j) } /* ---------------------------------------------------------------------- - read Kolmogorov-Crespi potential file + read DRIP file ------------------------------------------------------------------------- */ void PairDRIP::read_file(char *filename) @@ -893,7 +893,7 @@ void PairDRIP::read_file(char *filename) fp = force->open_potential(filename); if (fp == NULL) { char str[128]; - snprintf(str,128,"Cannot open KC potential file %s",filename); + snprintf(str,128,"Cannot open DRIP potential file %s",filename); error->one(FLERR,str); } } @@ -944,7 +944,7 @@ void PairDRIP::read_file(char *filename) } if (nwords != params_per_line) - error->all(FLERR,"Insufficient format in KC potential file"); + error->all(FLERR,"Insufficient format in DRIP potential file"); // words = ptrs to all words in line @@ -1001,9 +1001,9 @@ void PairDRIP::read_file(char *filename) //if(nparams >= pow(atom->ntypes,3)) break; } memory->destroy(elem2param); - memory->destroy(cutKCsq); + memory->destroy(cutDRIPsq); memory->create(elem2param,nelements,nelements,"pair:elem2param"); - memory->create(cutKCsq,nelements,nelements,"pair:cutKCsq"); + memory->create(cutDRIPsq,nelements,nelements,"pair:cutDRIPsq"); for (i = 0; i < nelements; i++) { for (j = 0; j < nelements; j++) { n = -1; @@ -1015,7 +1015,7 @@ void PairDRIP::read_file(char *filename) } if (n < 0) error->all(FLERR,"Potential file is missing an entry"); elem2param[i][j] = n; - cutKCsq[i][j] = params[n].rcut*params[n].rcut; + cutDRIPsq[i][j] = params[n].rcut*params[n].rcut; } } delete [] words; @@ -1023,40 +1023,6 @@ void PairDRIP::read_file(char *filename) /* ---------------------------------------------------------------------- */ -double PairDRIP::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, - double /*factor_coul*/, double factor_lj, - double &fforce) -{ - double r,r2inv,r6inv,r8inv,forcelj,philj; - double Tap,dTap,Vkc,fpair; - - int iparam_ij = elem2param[map[itype]][map[jtype]]; - Param& p = params[iparam_ij]; - - r = sqrt(rsq); - // turn on/off taper function - if (tap_flag) { - Tap = calc_Tap(r,sqrt(cutsq[itype][jtype])); - dTap = calc_dTap(r,sqrt(cutsq[itype][jtype])); - } else {Tap = 1.0; dTap = 0.0;} - - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; - r8inv = r2inv*r6inv; - - Vkc = -p.A*p.z06*r6inv; - // derivatives - fpair = -6.0*p.A*p.z06*r8inv; - forcelj = fpair; - fforce = factor_lj*(forcelj*Tap - Vkc*dTap/r); - - if (tap_flag) philj = Vkc*Tap; - else philj = Vkc - offset[itype][jtype]; - return factor_lj*philj; -} - -/* ---------------------------------------------------------------------- */ - int PairDRIP::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { diff --git a/src/USER-MISC/pair_drip.h b/src/USER-MISC/pair_drip.h index faa05bce4d..10750985c6 100644 --- a/src/USER-MISC/pair_drip.h +++ b/src/USER-MISC/pair_drip.h @@ -39,7 +39,6 @@ class PairDRIP : public Pair { void calc_normal(); int pack_forward_comm(int, int *, double *, int, int *); void unpack_forward_comm(int, int, double *); - double single(int, int, int, int, double, double, double, double &); protected: int me; @@ -47,8 +46,8 @@ class PairDRIP : public Pair { int pgsize; // size of neighbor page int oneatom; // max # of neighbors for one atom MyPage *ipage; // neighbor list pages - int *KC_numneigh; // # of pair neighbors for each atom - int **KC_firstneigh; // ptr to 1st neighbor of each atom + int *DRIP_numneigh; // # of pair neighbors for each atom + int **DRIP_firstneigh; // ptr to 1st neighbor of each atom int tap_flag; // flag to turn on/off taper function @@ -69,7 +68,7 @@ class PairDRIP : public Pair { double cut_global; double cut_normal; double **cut; - double **cutKCsq; + double **cutDRIPsq; double **offset; double **normal; double ***dnormdri; @@ -77,7 +76,7 @@ class PairDRIP : public Pair { void read_file( char * ); void allocate(); - void KC_neigh(); + void DRIP_neigh(); /* ----Calculate the long-range cutoff term */