diff --git a/src/KSPACE/pair_lj_long_tip4p_long.cpp b/src/KSPACE/pair_lj_long_tip4p_long.cpp index f2825c32b5..31623fc732 100755 --- a/src/KSPACE/pair_lj_long_tip4p_long.cpp +++ b/src/KSPACE/pair_lj_long_tip4p_long.cpp @@ -484,6 +484,7 @@ void PairLJLongTIP4PLong::compute_inner() int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; + // atom->nmax > nmax will occur during setup if (atom->nmax > nmax) { nmax = atom->nmax; memory->destroy(hneigh); @@ -999,6 +1000,19 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) // initialize hneigh[2] to 0 every step int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->destroy(hneigh); + memory->create(hneigh,nmax,3,"pair:hneigh"); + memory->destroy(newsite); + memory->create(newsite,nmax,3,"pair:newsite"); + } + if (neighbor->ago == 0) { + for (i = 0; i < nall; i++) hneigh[i][0] = -1; + for (i = 0; i < nall; i++) hneigh[i][2] = 0; + } double **f = atom->f; double **x = atom->x; @@ -1236,7 +1250,6 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) cforce = forcecoul * r2inv; fvirial = (forcecoul + respa_coul) * r2inv; - double fvcf = fvirial/cforce; // if i,j are not O atoms, force is applied directly // if i or j are O atoms, force is on fictitious atom & partitioned @@ -1291,15 +1304,27 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) f[iH2][2] += fH[2]; if (vflag) { + fd[0] = delx*fvirial; + fd[1] = dely*fvirial; + fd[2] = delz*fvirial; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + domain->closest_image(x[i],x[iH1],xH1); domain->closest_image(x[i],x[iH2],xH2); - v[0] = (x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0])*fvcf; - v[1] = (x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1])*fvcf; - v[2] = (x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2])*fvcf; - v[3] = (x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1])*fvcf; - v[4] = (x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2])*fvcf; - v[5] = (x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2])*fvcf; + v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] = x[i][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] = x[i][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] = x[i][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; } vlist[n++] = i; vlist[n++] = iH1; @@ -1349,15 +1374,28 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag) f[jH2][2] += fH[2]; if (vflag) { + + fd[0] = -delx*fvirial; + fd[1] = -dely*fvirial; + fd[2] = -delz*fvirial; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + domain->closest_image(x[j],x[jH1],xH1); domain->closest_image(x[j],x[jH2],xH2); - v[0] += (x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0])*fvcf; - v[1] += (x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1])*fvcf; - v[2] += (x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2])*fvcf; - v[3] += (x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1])*fvcf; - v[4] += (x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2])*fvcf; - v[5] += (x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2])*fvcf; + v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] += x[j][0]*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] += x[j][0]*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] += x[j][1]*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; } vlist[n++] = j; vlist[n++] = jH1; diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 6d631ba891..13ee4679c2 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -70,7 +70,7 @@ PPPMDisp::PPPMDisp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) triclinic_support = 0; pppmflag = dispersionflag = 1; - accuracy_relative = atof(arg[0]); + accuracy_relative = fabs(force->numeric(FLERR,arg[0])); nfactors = 3; factors = new int[nfactors]; @@ -1137,7 +1137,7 @@ void PPPMDisp::compute(int eflag, int vflag) if (slabflag) slabcorr(eflag); if (function[0]) energy += energy_1; - if (function[1]) energy += energy_6; + if (function[1] + function[2]) energy += energy_6; // convert atoms back from lamda to box coords @@ -2838,12 +2838,15 @@ void PPPMDisp::calc_csum() MPI_Allreduce(neach,neach_all,ntypes+1,MPI_INT,MPI_SUM,world); // copmute csumij and csumi - + double d1, d2; if (function[1]){ for (i=1; i<=ntypes; i++) { for (j=1; j<=ntypes; j++) { csumi[i] += neach_all[j]*B[i]*B[j]; - csumij += neach_all[i]*neach_all[j]*B[i]*B[j]; + d1 = neach_all[i]*B[i]; + d2 = neach_all[j]*B[j]; + csumij += d1*d2; + //csumij += neach_all[i]*neach_all[j]*B[i]*B[j]; } } } else { @@ -2851,7 +2854,10 @@ void PPPMDisp::calc_csum() for (j=1; j<=ntypes; j++) { for (k=0; k<=6; k++) { csumi[i] += neach_all[j]*B[7*i + k]*B[7*(j+1)-k-1]; - csumij += neach_all[i]*neach_all[j]*B[7*i + k]*B[7*(j+1)-k-1]; + d1 = neach_all[i]*B[7*i + k]; + d2 = neach_all[j]*B[7*(j+1)-k-1]; + csumij += d1*d2; + //csumij += neach_all[i]*neach_all[j]*B[7*i + k]*B[7*(j+1)-k-1]; } } } diff --git a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp index a86c077fc7..2c34e54c4d 100644 --- a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp +++ b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp @@ -425,6 +425,30 @@ void PairLJLongTIP4PLongOMP::compute_outer(int eflag, int vflag) const int order6 = ewald_order&(1<<6); const int nall = atom->nlocal + atom->nghost; + + // reallocate hneigh_thr & newsite_thr if necessary + // initialize hneigh_thr[0] to -1 on steps when reneighboring occured + // initialize hneigh_thr[2] to 0 every step + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memory->destroy(hneigh_thr); + memory->create(hneigh_thr,nmax,"pair:hneigh_thr"); + memory->destroy(newsite_thr); + memory->create(newsite_thr,nmax,"pair:newsite_thr"); + } + + int i; + // tag entire list as completely invalid after a neighbor + // list update, since that can change the order of atoms. + if (neighbor->ago == 0) { + for (i = 0; i < nall; i++) hneigh_thr[i].a = -1; + // indicate that the coordinates for the M point need to + // be updated. this needs to be done only if neighbor list + // has been updated in compute_outer + for (i = 0; i < nall; i++) hneigh_thr[i].t = 0; + } + const int nthreads = comm->nthreads; const int inum = listouter->inum; @@ -1798,7 +1822,6 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th cforce = forcecoul * r2inv; fvirial = (forcecoul + respa_coul) * r2inv; - double fvcf = fvirial/cforce; // if i,j are not O atoms, force is applied directly // if i or j are O atoms, force is on fictitious atom & partitioned @@ -1853,15 +1876,28 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th f[iH2][2] += fH[2]; if (EVFLAG) { + + fd[0] = delx*fvirial; + fd[1] = dely*fvirial; + fd[2] = delz*fvirial; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + domain->closest_image(&x[i].x,&x[iH1].x,xH1); domain->closest_image(&x[i].x,&x[iH2].x,xH2); - v[0] = (x[i].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0])*fvcf; - v[1] = (x[i].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1])*fvcf; - v[2] = (x[i].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2])*fvcf; - v[3] = (x[i].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1])*fvcf; - v[4] = (x[i].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2])*fvcf; - v[5] = (x[i].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2])*fvcf; + v[0] = x[i].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] = x[i].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] = x[i].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] = x[i].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] = x[i].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] = x[i].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; } vlist[n++] = i; vlist[n++] = iH1; @@ -1911,15 +1947,28 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th f[jH2][2] += fH[2]; if (EVFLAG) { + + fd[0] = -delx*fvirial; + fd[1] = -dely*fvirial; + fd[2] = -delz*fvirial; + + fO[0] = fd[0]*(1 - alpha); + fO[1] = fd[1]*(1 - alpha); + fO[2] = fd[2]*(1 - alpha); + + fH[0] = 0.5 * alpha * fd[0]; + fH[1] = 0.5 * alpha * fd[1]; + fH[2] = 0.5 * alpha * fd[2]; + domain->closest_image(&x[j].x,&x[jH1].x,xH1); domain->closest_image(&x[j].x,&x[jH2].x,xH2); - v[0] += (x[j].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0])*fvcf; - v[1] += (x[j].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1])*fvcf; - v[2] += (x[j].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2])*fvcf; - v[3] += (x[j].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1])*fvcf; - v[4] += (x[j].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2])*fvcf; - v[5] += (x[j].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2])*fvcf; + v[0] += x[j].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0]; + v[1] += x[j].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1]; + v[2] += x[j].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2]; + v[3] += x[j].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1]; + v[4] += x[j].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2]; + v[5] += x[j].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2]; } vlist[n++] = j; vlist[n++] = jH1;