diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index f851325a3a..f0e9a86ef0 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -41,6 +41,7 @@ using namespace MathConst; #define MAXLINE 1024 #define DELTA 4 #define PGDELTA 1 +#define MAXNEIGH 24 /* ---------------------------------------------------------------------- */ @@ -132,26 +133,15 @@ void PairComb::compute(int eflag, int vflag) double yaself; double potal,fac11,fac11e; double vionij,fvionij,sr1,sr2,sr3,Eov,Fov; - int sht_jnum, *sht_jlist; + int sht_jnum, *sht_jlist, nj; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = vflag_atom = 0; // Build short range neighbor list - // int every=neighbor->every; - // int ntimestep=update->ntimestep; - // if(ntimestep <= 1 || (ntimestep % every == 0)) - Short_neigh(); - // grow coordination array if necessary - if (atom->nmax > nmax) { - memory->destroy(NCo); - memory->destroy(bbij); - nmax = atom->nmax; - memory->create(NCo,nmax,"pair:NCo"); - memory->create(bbij,nmax,nmax,"pair:bbij"); - } + Short_neigh(); double **x = atom->x; double **f = atom->f; @@ -183,6 +173,7 @@ void PairComb::compute(int eflag, int vflag) ztmp = x[i][2]; iq = q[i]; NCo[i] = 0; + nj = 0; iparam_i = elem2param[itype][itype][itype]; // self energy, only on i atom @@ -279,8 +270,7 @@ void PairComb::compute(int eflag, int vflag) if (cor_flag) { for (jj = 0; jj < sht_jnum; jj++) { - j = sht_jlist[jj]; - j &= NEIGHMASK; + j = sht_jlist[jj]; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -301,7 +291,6 @@ void PairComb::compute(int eflag, int vflag) for (jj = 0; jj < sht_jnum; jj++) { j = sht_jlist[jj]; - j &= NEIGHMASK; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -316,6 +305,7 @@ void PairComb::compute(int eflag, int vflag) rsq1 = vec3_dot(delr1,delr1); if (rsq1 > params[iparam_ij].cutsq) continue; + nj ++; // accumulate bondorder zeta for each i-j interaction via loop over k @@ -325,7 +315,6 @@ void PairComb::compute(int eflag, int vflag) for (kk = 0; kk < sht_jnum; kk++) { k = sht_jlist[kk]; if (j == k) continue; - k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -345,7 +334,7 @@ void PairComb::compute(int eflag, int vflag) if (cuo_flag1 && cuo_flag2) cuo_flag = 1; else cuo_flag = 0; - force_zeta(¶ms[iparam_ij],eflag,i,j,rsq1,zeta_ij, + force_zeta(¶ms[iparam_ij],eflag,i,nj,rsq1,zeta_ij, iq,jq,fpair,prefactor,evdwl); // over-coordination correction for HfO2 @@ -370,7 +359,6 @@ void PairComb::compute(int eflag, int vflag) for (kk = 0; kk < sht_jnum; kk++) { k = sht_jlist[kk]; if (j == k) continue; - k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -551,7 +539,6 @@ void PairComb::init_style() int irequest = neighbor->request(this); neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->full = 1; - neighbor->requests[irequest]->ghost = 1; pgsize = neighbor->pgsize; oneatom = neighbor->oneatom; @@ -822,6 +809,8 @@ void PairComb::setup() params[m].lcut = params[m].coulcut; params[m].lcutsq = params[m].lcut*params[m].lcut; + + params[m].gamma = 1.0; // for the change in pair_comb.h } // set cutmax to max of all params @@ -831,7 +820,7 @@ void PairComb::setup() for (m = 0; m < nparams; m++) { if (params[m].cut > cutmax) cutmax = params[m].cut; if (params[m].lcut > cutmax) cutmax = params[m].lcut; - if (params[m].cutsq > cutmin) cutmin = params[m].cutsq+1.0; + if (params[m].cutsq > cutmin) cutmin = params[m].cutsq+0.2; if (params[m].hfocor > 0.0001) cor_flag = 1; } } @@ -958,7 +947,7 @@ double PairComb::elp(Param *param, double rsqij, double rsqik, comtt += param->aconf *(4.0-(rmu-c123)*(rmu-c123)); } - return 1.0 * fck * comtt; + return 0.5 * fck * comtt; } return 0.0; @@ -1024,10 +1013,10 @@ void PairComb::flp(Param *param, double rsqij, double rsqik, com4k = fcj * fck_d * comtt; com5 = fcj * fck * comtt_d; - ffj1 =-1.0*(com5/(rij*rik)); - ffj2 = 1.0*(com5*rmu/rsqij); + ffj1 =-0.5*(com5/(rij*rik)); + ffj2 = 0.5*(com5*rmu/rsqij); ffk1 = ffj1; - ffk2 = 1.0*(-com4k/rik+com5*rmu/rsqik); + ffk2 = 0.5*(-com4k/rik+com5*rmu/rsqik); } else { ffj1 = 0.0; ffj2 = 0.0; @@ -1361,15 +1350,7 @@ void PairComb::sm_table() nntypes = int((n+1)*n/2); // interaction types ncoul = int((rc-drin)/dra)+1; -/* - memory->destroy(intype); - memory->destroy(fafb); - memory->destroy(dfafb); - memory->destroy(ddfafb); - memory->destroy(phin); - memory->destroy(dphin); - memory->destroy(erpaw); -*/ + // allocate arrays memory->create(intype,n,n,"pair:intype"); @@ -1380,7 +1361,7 @@ void PairComb::sm_table() memory->create(dphin,ncoul,nntypes,"pair:dphin"); memory->create(erpaw,25000,2,"pair:erpaw"); memory->create(NCo,nmax,"pair:NCo"); - memory->create(bbij,nmax,nmax,"pair:bbij"); + memory->create(bbij,nmax,MAXNEIGH,"pair:bbij"); memory->create(sht_num,nmax,"pair:sht_num"); sht_first = (int **) memory->smalloc(nmax*sizeof(int *), "pair:sht_first"); @@ -1636,8 +1617,8 @@ double PairComb::yasu_char(double *qf_fix, int &igroup) int *ilist,*jlist,*numneigh,**firstneigh; double iq,jq,fqi,fqj,fqij,fqjj; double potal,fac11,fac11e,sr1,sr2,sr3; - int mr1,mr2,mr3,inty; - int sht_jnum, *sht_jlist; + int mr1,mr2,mr3,inty,nj; + double **x = atom->x; double *q = atom->q; @@ -1674,6 +1655,7 @@ double PairComb::yasu_char(double *qf_fix, int &igroup) for (ii = 0; ii < inum; ii ++) { i = ilist[ii]; itag = tag[i]; + nj = 0; if (mask[i] & groupbit) { itype = map[type[i]]; xtmp = x[i][0]; @@ -1690,8 +1672,6 @@ double PairComb::yasu_char(double *qf_fix, int &igroup) jlist = firstneigh[i]; jnum = numneigh[i]; - sht_jlist = sht_first[i]; - sht_jnum = sht_num[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -1756,10 +1736,11 @@ double PairComb::yasu_char(double *qf_fix, int &igroup) iparam_ij = elem2param[itype][jtype][jtype]; if (rsq1 > params[iparam_ij].cutsq) continue; + nj ++; // charge force in Aij and Bij - qfo_short(¶ms[iparam_ij],i,j,rsq1,iq,jq,fqij,fqjj); + qfo_short(¶ms[iparam_ij],i,nj,rsq1,iq,jq,fqij,fqjj); fqi += fqij; qf[j] += fqjj; } qf[i] += fqi; @@ -2050,7 +2031,7 @@ double PairComb::memory_usage() double bytes = maxeatom * sizeof(double); bytes += maxvatom*6 * sizeof(double); bytes += nmax * sizeof(int); - bytes += nmax * nmax * sizeof(int); + bytes += MAXNEIGH * nmax * sizeof(int); return bytes; } /* ---------------------------------------------------------------------- */ @@ -2069,12 +2050,13 @@ void PairComb::Short_neigh() int ntype = atom->ntypes; if (atom->nmax > nmax) { - nmax = int(1.0 * atom->nmax); - memory->sfree(sht_num); memory->sfree(sht_first); - memory->create(sht_num,nmax,"pair:sht_num"); + nmax = atom->nmax; sht_first = (int **) memory->smalloc(nmax*sizeof(int *), - "pair:sht_first"); + "pair:sht_first"); + memory->grow(sht_num,nmax,"pair:sht_num"); + memory->grow(NCo,nmax,"pair:NCo"); + memory->grow(bbij,nmax,MAXNEIGH,"pair:bbij"); } inum = list->inum; @@ -2108,7 +2090,6 @@ void PairComb::Short_neigh() j = jlist[jj]; j &= NEIGHMASK; jtype = type[j]; - iparam_ij = elem2param[itype][jtype][jtype]; delrj[0] = xtmp - x[j][0]; delrj[1] = ytmp - x[j][1];