From e07b46c77158dffa400ec5d1b9c2e556c365a667 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 17 May 2022 14:40:09 -0600 Subject: [PATCH] refactoring edits --- examples/amoeba/Compare.sh | 56 ++++----- src/AMOEBA/amoeba_induce.cpp | 215 +---------------------------------- src/AMOEBA/amoeba_polar.cpp | 17 --- src/AMOEBA/amoeba_utils.cpp | 56 +++++---- src/AMOEBA/pair_amoeba.cpp | 151 ++++++------------------ src/AMOEBA/pair_amoeba.h | 2 +- 6 files changed, 89 insertions(+), 408 deletions(-) diff --git a/examples/amoeba/Compare.sh b/examples/amoeba/Compare.sh index 9c57bb8d92..8eaebcee30 100644 --- a/examples/amoeba/Compare.sh +++ b/examples/amoeba/Compare.sh @@ -2,51 +2,51 @@ # dimer -kdiff log.water_dimer.amoeba.1 log.water_dimer.amoeba.1.test -kdiff dump.water_dimer.amoeba.1 dump.water_dimer.amoeba.1.test +kdiff.py log.water_dimer.amoeba.1 log.water_dimer.amoeba.1.test +kdiff.py dump.water_dimer.amoeba.1 dump.water_dimer.amoeba.1.test -kdiff log.water_dimer.amoeba.4 log.water_dimer.amoeba.4.test -kdiff dump.water_dimer.amoeba.4 dump.water_dimer.amoeba.4.test +kdiff.py log.water_dimer.amoeba.4 log.water_dimer.amoeba.4.test +kdiff.py dump.water_dimer.amoeba.4 dump.water_dimer.amoeba.4.test -kdiff log.water_dimer.hippo.1 log.water_dimer.hippo.1.test -kdiff dump.water_dimer.hippo.1 dump.water_dimer.hippo.1.test +kdiff.py log.water_dimer.hippo.1 log.water_dimer.hippo.1.test +kdiff.py dump.water_dimer.hippo.1 dump.water_dimer.hippo.1.test -kdiff log.water_dimer.hippo.4 log.water_dimer.hippo.4.test -kdiff dump.water_dimer.hippo.4 dump.water_dimer.hippo.4.test +kdiff.py log.water_dimer.hippo.4 log.water_dimer.hippo.4.test +kdiff.py dump.water_dimer.hippo.4 dump.water_dimer.hippo.4.test # hexamer -kdiff log.water_hexamer.amoeba.1 log.water_hexamer.amoeba.1.test -kdiff dump.water_hexamer.amoeba.1 dump.water_hexamer.amoeba.1.test +kdiff.py log.water_hexamer.amoeba.1 log.water_hexamer.amoeba.1.test +kdiff.py dump.water_hexamer.amoeba.1 dump.water_hexamer.amoeba.1.test -kdiff log.water_hexamer.amoeba.4 log.water_hexamer.amoeba.4.test -kdiff dump.water_hexamer.amoeba.4 dump.water_hexamer.amoeba.4.test +kdiff.py log.water_hexamer.amoeba.4 log.water_hexamer.amoeba.4.test +kdiff.py dump.water_hexamer.amoeba.4 dump.water_hexamer.amoeba.4.test -kdiff log.water_hexamer.hippo.1 log.water_hexamer.hippo.1.test -kdiff dump.water_hexamer.hippo.1 dump.water_hexamer.hippo.1.test +kdiff.py log.water_hexamer.hippo.1 log.water_hexamer.hippo.1.test +kdiff.py dump.water_hexamer.hippo.1 dump.water_hexamer.hippo.1.test -kdiff log.water_hexamer.hippo.4 log.water_hexamer.hippo.4.test -kdiff dump.water_hexamer.hippo.4 dump.water_dimer.hippo.4.test +kdiff.py log.water_hexamer.hippo.4 log.water_hexamer.hippo.4.test +kdiff.py dump.water_hexamer.hippo.4 dump.water_hexamer.hippo.4.test # water box -kdiff log.water_box.amoeba.1 log.water_box.amoeba.1.test -kdiff dump.water_box.amoeba.1 dump.water_box.amoeba.1.test +kdiff.py log.water_box.amoeba.1 log.water_box.amoeba.1.test +kdiff.py dump.water_box.amoeba.1 dump.water_box.amoeba.1.test -kdiff log.water_box.amoeba.32 log.water_box.amoeba.32.test -kdiff dump.water_box.amoeba.32 dump.water_box.amoeba.32.test +kdiff.py log.water_box.amoeba.32 log.water_box.amoeba.32.test +kdiff.py dump.water_box.amoeba.32 dump.water_box.amoeba.32.test -kdiff log.water_box.hippo.1 log.water_box.hippo.1.test -kdiff dump.water_box.hippo.1 dump.water_box.hippo.1.test +kdiff.py log.water_box.hippo.1 log.water_box.hippo.1.test +kdiff.py dump.water_box.hippo.1 dump.water_box.hippo.1.test -kdiff log.water_box.hippo.32 log.water_box.hippo.32.test -kdiff dump.water_box.hippo.32 dump.water_box.hippo.32.test +kdiff.py log.water_box.hippo.32 log.water_box.hippo.32.test +kdiff.py dump.water_box.hippo.32 dump.water_box.hippo.32.test # ubiquitin -kdiff log.ubi.1 log.ubi.1.test -kdiff dump.ubi.1 dump.ubi.1.test +kdiff.py log.ubi.1 log.ubi.1.test +kdiff.py dump.ubi.1 dump.ubi.1.test -kdiff log.ubi.32 log.ubi.32.test -kdiff dump.ubi.32 dump.ubi.32.test +kdiff.py log.ubi.32 log.ubi.32.test +kdiff.py dump.ubi.32 dump.ubi.32.test diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp index cbd0a2c480..c23e6c4516 100644 --- a/src/AMOEBA/amoeba_induce.cpp +++ b/src/AMOEBA/amoeba_induce.cpp @@ -40,8 +40,6 @@ enum{GORDON1,GORDON2}; #define DEBYE 4.80321 // conversion factor from q-Angs (real units) to Debye -#define UIND_DEBUG 0 - /* ---------------------------------------------------------------------- induce = induced dipole moments via pre-conditioned CG solver adapted from Tinker induce0a() routine @@ -113,17 +111,6 @@ void PairAmoeba::induce() crstyle = FIELD; comm->reverse_comm(this); - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("AAA FIELD atom %d: field %g %g %g: fieldp %g %g %g\n", - atom->tag[i], - field[i][0],field[i][1],field[i][2], - fieldp[i][0],fieldp[i][1],fieldp[i][2]); - */ - // set induced dipoles to polarizability times direct field for (i = 0; i < nlocal; i++) { @@ -138,17 +125,6 @@ void PairAmoeba::induce() } } - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("AAA UDIR atom %d: udir %g %g %g: udirp %g %g %g\n", - atom->tag[i], - DEBYE*udir[i][0],DEBYE*udir[i][1],DEBYE*udir[i][2], - DEBYE*udirp[i][0],DEBYE*udirp[i][1],DEBYE*udirp[i][2]); - */ - // get induced dipoles via the OPT extrapolation method // NOTE: any way to rewrite these loops to avoid allocating // uopt,uoptp with a optorder+1 dimension, just optorder ?? @@ -258,20 +234,6 @@ void PairAmoeba::induce() crstyle = FIELD; comm->reverse_comm(this); - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("UFIELD atom %d: uind %g %g %g uinp %g %g %g " - "field %g %g %g: fieldp %g %g %g\n", - atom->tag[i], - uind[i][0],uind[i][1],uind[i][2], - uinp[i][0],uinp[i][1],uinp[i][2], - field[i][0],field[i][1],field[i][2], - fieldp[i][0],fieldp[i][1],fieldp[i][2]); - */ - // set initial conjugate gradient residual and conjugate vector for (i = 0; i < nlocal; i++) { @@ -307,20 +269,6 @@ void PairAmoeba::induce() } } - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("CONJ atom %d: rsd %g %g %g rsdp %g %g %g " - "conj %g %g %g: conjp %g %g %g\n", - atom->tag[i], - rsd[i][0],rsd[i][1],rsd[i][2], - rsdp[i][0],rsdp[i][1],rsdp[i][2], - conj[i][0],conj[i][1],conj[i][2], - conjp[i][0],conjp[i][1],conjp[i][2]); - */ - // conjugate gradient iteration of the mutual induced dipoles while (!done) { @@ -345,19 +293,6 @@ void PairAmoeba::induce() crstyle = FIELD; comm->reverse_comm(this); - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("POST-COMM FIELD atom %d: field %g %g %g: fieldp %g %g %g\n", - atom->tag[i], - uind[i][0],uind[i][1],uind[i][2], - uinp[i][0],uinp[i][1],uinp[i][2], - field[i][0],field[i][1],field[i][2], - fieldp[i][0],fieldp[i][1],fieldp[i][2]); - */ - for (i = 0; i < nlocal; i++) { for (j = 0; j < 3; j++) { uind[i][j] = vec[i][j]; @@ -405,17 +340,6 @@ void PairAmoeba::induce() } } - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("POST-MPI UIND atom %d: uind %g %g %g: uinp %g %g %g\n", - atom->tag[i], - uind[i][0],uind[i][1],uind[i][2], - uinp[i][0],uinp[i][1],uinp[i][2]); - */ - if (pcgprec) { cfstyle = RSD; comm->forward_comm(this); @@ -424,20 +348,6 @@ void PairAmoeba::induce() comm->reverse_comm(this); } - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("POST-PRECOND atom %d: rsd %g %g %g: rsdp %g %g %g " - "zrsd %g %g %g: zrsdp %g %g %g\n", - atom->tag[i], - rsd[i][0],rsd[i][1],rsd[i][2], - rsdp[i][0],rsdp[i][1],rsdp[i][2], - zrsd[i][0],zrsd[i][1],zrsd[i][2], - zrsdp[i][0],zrsdp[i][1],zrsdp[i][2]); - */ - b = 0.0; bp = 0.0; @@ -483,32 +393,12 @@ void PairAmoeba::induce() eps = MAX(epsd,epsp); eps = DEBYE * sqrt(eps/atom->natoms); - /* - if (debug) { - if (comm->me == 0 && screen) { - fprintf(screen,"SCF induced dipole moments: " - "iter %d, RMS residual (Debye) %g\n",iter,eps); - } - } - */ - if (eps < poleps) done = true; if (eps > epsold) done = true; if (iter >= politer) done = true; // apply a "peek" iteration to the mutual induced dipoles - // DEBUG statements - - /* - printf("DONE %d\n",done); - for (i = 0; i < nlocal; i++) - printf("PRE-DONE UIND atom %d: uind %g %g %g: uinp %g %g %g\n", - atom->tag[i], - uind[i][0],uind[i][1],uind[i][2], - uinp[i][0],uinp[i][1],uinp[i][2]); - */ - if (done) { for (i = 0; i < nlocal; i++) { term = pcgpeek * poli[i]; @@ -518,17 +408,6 @@ void PairAmoeba::induce() } } } - - // DEBUG statements - - /* - for (i = 0; i < nlocal; i++) - if (atom->tag[i] == 1) - printf("POST-DONE UIND atom %d: uind %g %g %g: uinp %g %g %g\n", - atom->tag[i], - DEBYE*uind[i][0],DEBYE*uind[i][1],DEBYE*uind[i][2], - DEBYE*uinp[i][0],DEBYE*uinp[i][1],DEBYE*uinp[i][2]); - */ } // terminate the calculation if dipoles failed to converge @@ -539,11 +418,6 @@ void PairAmoeba::induce() error->warning(FLERR,"AMOEBA induced dipoles did not converge"); } - // DEBUG output to dump file - - if (UIND_DEBUG) - dump6(fp_uind,"id uindx uindy uindz uinpx uinpy uinpz",DEBYE,uind,uinp); - // deallocation of arrays memory->destroy(poli); @@ -731,14 +605,6 @@ void PairAmoeba::ufield0c(double **field, double **fieldp) field[i][j] += term*uind[i][j]; fieldp[i][j] += term*uinp[i][j]; } - - /* - // DEBUG - - printf("UMUTUAL2B SELF i %d term %g aewald %g uind %g %g %g field %g %g %g\n", - atom->tag[i],term,aewald, - uind[i][0],uind[i][1],uind[i][2],field[i][0],field[i][1],field[i][2]); - */ } } @@ -939,18 +805,6 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp, pclist[3] = rr5*yr*yr - rr3; pclist[4] = rr5*yr*zr; pclist[5] = rr5*zr*zr - rr3; - - // DEBUG - - /* - printf("PCLIST: ij %d %d: pc1-6 %g %g %g %g %g %g\n", - atom->tag[i],atom->tag[j], - pclist[0],pclist[1],pclist[2],pclist[3],pclist[4],pclist[5]); - printf(" AMOEBA: scale35 %g %g pdamp %g thole %g " - "pgamma %g facUscale %g damp %g\n", - scale3,scale5,pdamp[jtype],thole[jtype],pgamma,factor_uscale,damp); - */ - pclist += 6; } } @@ -1054,11 +908,6 @@ void PairAmoeba::umutual1(double **field, double **fieldp) } } - /* - printf("fuind %g %g %g \n",fuind[0][0],fuind[0][1],fuind[0][2]); - printf("fuinp %g %g %g \n",fuinp[0][0],fuinp[0][1],fuinp[0][2]); - */ - // gridpre = my portion of 4d grid in brick decomp w/ ghost values double ****gridpre = (double ****) ic_kspace->zero(); @@ -1106,16 +955,6 @@ void PairAmoeba::umutual1(double **field, double **fieldp) fphi_uind(gridpost,fdip_phi1,fdip_phi2,fdip_sum_phi); - // printf ("fdip_phi1_uind %g %g %g %g %g %g %g %g %g %g\n", - // fdip_phi1[0][0], - // fdip_phi1[0][1],fdip_phi1[0][2],fdip_phi1[0][3],fdip_phi1[0][4],fdip_phi1[0][5], - // fdip_phi1[0][6],fdip_phi1[0][7],fdip_phi1[0][8],fdip_phi1[0][9],fdip_phi1[0][10]); - - // printf ("fdip_phi2_uind %g %g %g %g %g %g %g %g %g %g\n", - // fdip_phi2[0][0], - // fdip_phi2[0][1],fdip_phi2[0][2],fdip_phi2[0][3],fdip_phi2[0][4],fdip_phi2[0][5], - // fdip_phi2[0][6],fdip_phi2[0][7],fdip_phi2[0][8],fdip_phi2[0][9],fdip_phi2[0][10]); - // store fractional reciprocal potentials for OPT method if (poltyp == OPT) { @@ -1217,21 +1056,6 @@ void PairAmoeba::umutual2b(double **field, double **fieldp) fkp[1] = tdipdip[1]*uinpi[0] + tdipdip[3]*uinpi[1] + tdipdip[4]*uinpi[2]; fkp[2] = tdipdip[2]*uinpi[0] + tdipdip[4]*uinpi[1] + tdipdip[5]*uinpi[2]; - // DEBUG - - /* - if (atom->tag[i] == 1 || atom->tag[j] == 1) { - printf ("TDIPDIP ij %d %d: tdd %g %g %g %g %g %g\n", - atom->tag[i],atom->tag[j], - tdipdip[0],tdipdip[1],tdipdip[2], - tdipdip[3],tdipdip[4],tdipdip[5]); - printf ("FIDFKD ij %d %d: fid %g %g %g fkd %g %g %g\n", - atom->tag[i],atom->tag[j], - fid[0],fid[1],fid[2], - fkd[0],fkd[1],fkd[2]); - } - */ - tdipdip += 6; // increment the field at each site due to this interaction @@ -1375,11 +1199,6 @@ void PairAmoeba::udirect1(double **field) fphi_mpole(gridpost,fphi); - // printf ("fphi %g %g %g %g %g %g %g %g %g %g %g %g\n", - // fphi[1][1],fphi[1][2],fphi[1][3],fphi[1][4],fphi[1][5], - // fphi[1][6],fphi[1][7],fphi[1][8],fphi[1][9],fphi[1][10], - // fphi[1][11],fphi[1][12]); - // convert the field from fractional to Cartesian fphi_to_cphi(fphi,cphi); @@ -1450,7 +1269,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp) numneigh = list->numneigh; firstneigh = list->firstneigh; - // NOTE: doesn't this have a problem if aewald is tiny ?? + // NOTE: does this have a problem if aewald is tiny ?? aesq2 = 2.0 * aewald * aewald; aesq2n = 0.0; @@ -1496,18 +1315,6 @@ void PairAmoeba::udirect2b(double **field, double **fieldp) vali = pval[i]; } - // DEBUG - - /* - if (atom->tag[i] == 4) { - printf("Atom 4 is in group %d\n",igroup); - printf("Atoms in same group:"); - for (int ig = 0; ig < atom->nlocal; ig++) - if (amgroup[ig] == igroup) printf(" %d",atom->tag[ig]); - printf("\n"); - } - */ - // evaluate all sites within the cutoff distance for (jj = 0; jj < jnum; jj++) { @@ -1546,16 +1353,6 @@ void PairAmoeba::udirect2b(double **field, double **fieldp) } } - // DEBUG - - /* - if (atom->tag[i] == 4 || atom->tag[j] == 4) { - printf("PAIR ij %d %d ij group %d %d wpdu scale %g %g %g %g\n", - atom->tag[i],atom->tag[j],igroup,jgroup, - factor_wscale,factor_pscale,factor_dscale,factor_uscale); - } - */ - r = sqrt(r2); rr1 = 1.0 / r; rr2 = rr1 * rr1; @@ -1749,16 +1546,6 @@ void PairAmoeba::udirect2b(double **field, double **fieldp) rr3ik = bn[1] - (1.0-scalek*dmpik[2])*rr3; rr5ik = bn[2] - (1.0-scalek*dmpik[4])*rr5; - // DEBUG - - /* - if (atom->tag[i] == 1 || atom->tag[j] == 1) - printf ("DAMPMUT ij %d %d: bn12 %g %g dmpik-01234 %g %g %g %g %g\n", - atom->tag[i],atom->tag[j], - bn[1],bn[2], - dmpik[0],dmpik[1],dmpik[2],dmpik[3],dmpik[4]); - */ - neighptr[n++] = j; tdipdip[ndip++] = -rr3ik + rr5ik*xr*xr; tdipdip[ndip++] = rr5ik*xr*yr; diff --git a/src/AMOEBA/amoeba_polar.cpp b/src/AMOEBA/amoeba_polar.cpp index d067fd5dc0..302121ebc6 100644 --- a/src/AMOEBA/amoeba_polar.cpp +++ b/src/AMOEBA/amoeba_polar.cpp @@ -23,10 +23,6 @@ #include "math_const.h" #include "memory.h" -// DEBUG - -#include "error.h" - using namespace LAMMPS_NS; using namespace MathConst; @@ -1436,8 +1432,6 @@ void PairAmoeba::polar_kspace() } } - //PRINT fuind - // gridpre2 = my portion of 4d grid in brick decomp w/ ghost values double ****gridpre2 = (double ****) pc_kspace->zero(); @@ -1880,17 +1874,6 @@ void PairAmoeba::polar_kspace() double ***gridpre = (double ***) p_kspace->zero(); - // DEBUG - - double psum = 0.0; - for (k = p_kspace->nzlo_out; k <= p_kspace->nzhi_out; k++) { - for (j = p_kspace->nylo_out; j <= p_kspace->nyhi_out; j++) { - for (i = p_kspace->nxlo_out; i <= p_kspace->nxhi_out; i++) { - psum += gridpre[k][j][i]; - } - } - } - // map atoms to grid grid_mpole(fmp,gridpre); diff --git a/src/AMOEBA/amoeba_utils.cpp b/src/AMOEBA/amoeba_utils.cpp index c52b8bbcac..059587a093 100644 --- a/src/AMOEBA/amoeba_utils.cpp +++ b/src/AMOEBA/amoeba_utils.cpp @@ -30,7 +30,7 @@ enum{NOFRAME,ZONLY,ZTHENX,BISECTOR,ZBISECT,THREEFOLD}; /* ---------------------------------------------------------------------- kmpole performs one-time assignment of - xaxis,yaxis,zaxis multipole neighbors to each owned atom + xyzaxis multipole neighbors to each owned atom any of the values can be 0 if not used yaxis can later be negative due to chkpole() also sets polaxe and pole[13] multipole for each owned atom @@ -52,9 +52,7 @@ void PairAmoeba::kmpole() amtype = atom->ivector[index_amtype]; int *polaxe = atom->ivector[index_polaxe]; - double *xaxis = atom->dvector[index_xaxis]; - double *yaxis = atom->dvector[index_yaxis]; - double *zaxis = atom->dvector[index_zaxis]; + double **xyzaxis = atom->darray[index_xyzaxis]; double **pole = fixpole->astore; int **nspecial = atom->nspecial; @@ -125,9 +123,9 @@ void PairAmoeba::kmpole() if (ktype == xtype) { if (ytype == 0 && !flag) { flag = 1; - zaxis[i] = ubuf(jneigh).d; - xaxis[i] = ubuf(kneigh).d; - yaxis[i] = 0.0; + xyzaxis[i][2] = ubuf(jneigh).d; + xyzaxis[i][0] = ubuf(kneigh).d; + xyzaxis[i][1] = 0.0; polaxe[i] = mpaxis[itype][iframe]; for (j = 0; j < 13; j++) pole[i][j] = fpole[itype][iframe][j]; @@ -141,9 +139,9 @@ void PairAmoeba::kmpole() mtype = amtype[m]; if (mtype == ytype && !flag) { flag = 1; - zaxis[i] = ubuf(jneigh).d; - xaxis[i] = ubuf(kneigh).d; - yaxis[i] = ubuf(mneigh).d; + xyzaxis[i][2] = ubuf(jneigh).d; + xyzaxis[i][0] = ubuf(kneigh).d; + xyzaxis[i][1] = ubuf(mneigh).d; polaxe[i] = mpaxis[itype][iframe]; for (j = 0; j < 13; j++) pole[i][j] = fpole[itype][iframe][j]; @@ -183,9 +181,9 @@ void PairAmoeba::kmpole() if (ktype == xtype) { if (ytype == 0 && !flag) { flag = 1; - zaxis[i] = ubuf(jneigh).d; - xaxis[i] = ubuf(kneigh).d; - yaxis[i] = 0.0; + xyzaxis[i][2] = ubuf(jneigh).d; + xyzaxis[i][0] = ubuf(kneigh).d; + xyzaxis[i][1] = 0.0; polaxe[i] = mpaxis[itype][iframe]; for (j = 0; j < 13; j++) pole[i][j] = fpole[itype][iframe][j]; @@ -203,9 +201,9 @@ void PairAmoeba::kmpole() if (!path) continue; if (mtype == ytype && !flag) { flag = 1; - zaxis[i] = ubuf(jneigh).d; - xaxis[i] = ubuf(kneigh).d; - yaxis[i] = ubuf(mneigh).d; + xyzaxis[i][2] = ubuf(jneigh).d; + xyzaxis[i][0] = ubuf(kneigh).d; + xyzaxis[i][1] = ubuf(mneigh).d; polaxe[i] = mpaxis[itype][iframe]; for (j = 0; j < 13; j++) pole[i][j] = fpole[itype][iframe][j]; @@ -234,8 +232,8 @@ void PairAmoeba::kmpole() if (jtype == ztype) { if (xtype == 0 && !flag) { flag = 1; - zaxis[i] = ubuf(jneigh).d; - xaxis[i] = yaxis[i] = 0.0; + xyzaxis[i][2] = ubuf(jneigh).d; + xyzaxis[i][0] = xyzaxis[i][1] = 0.0; polaxe[i] = mpaxis[itype][iframe]; for (j = 0; j < 13; j++) pole[i][j] = fpole[itype][iframe][j]; @@ -254,7 +252,7 @@ void PairAmoeba::kmpole() ztype = zpole[itype][iframe]; if (ztype == 0 && !flag) { flag = 1; - zaxis[i] = xaxis[i] = yaxis[i] = 0.0; + xyzaxis[i][2] = xyzaxis[i][0] = xyzaxis[i][2] = 0.0; polaxe[i] = mpaxis[itype][iframe]; for (j = 0; j < 13; j++) pole[i][j] = fpole[itype][iframe][j]; @@ -297,8 +295,8 @@ void PairAmoeba::chkpole(int i) double **pole = fixpole->astore; int *polaxe = atom->ivector[index_polaxe]; - double *yaxis = atom->dvector[index_yaxis]; - tagint yaxisID = (tagint) ubuf(yaxis[i]).i; + double **xyzaxis = atom->darray[index_xyzaxis]; + tagint yaxisID = (tagint) ubuf(xyzaxis[i][1]).i; // test for chirality inversion // if not, return @@ -334,7 +332,7 @@ void PairAmoeba::chkpole(int i) // flip sign in permanent yaxis, not yaxis2local if ((yaxisID < 0 && vol > 0.0) || (yaxisID > 0 && vol < 0.0)) { - yaxis[i] = -ubuf(yaxisID).d; + xyzaxis[i][1] = -ubuf(yaxisID).d; pole[i][2] = -pole[i][2]; pole[i][5] = -pole[i][5]; pole[i][7] = -pole[i][7]; @@ -716,21 +714,19 @@ void PairAmoeba::find_multipole_neighbors() tagint xaxisID,yaxisID,zaxisID; // grab current pts for xaxis,yaxis,zaxis - // xyz axis[i] = atom IDs that atom I uses for its multipole orientation + // xyzaxis[i] = atom IDs that atom I uses for its multipole orientation // can be zero if not used, in which case set local index to self // yaxis can be negative, in which case use absolute value - double *xaxis = atom->dvector[index_xaxis]; - double *yaxis = atom->dvector[index_yaxis]; - double *zaxis = atom->dvector[index_zaxis]; + double **xyzaxis = atom->darray[index_xyzaxis]; int nlocal = atom->nlocal; int nmissing = 0; for (int i = 0; i < nlocal; i++) { - xaxisID = (tagint) ubuf(xaxis[i]).i; - yaxisID = (tagint) ubuf(yaxis[i]).i; - zaxisID = (tagint) ubuf(zaxis[i]).i; + xaxisID = (tagint) ubuf(xyzaxis[i][0]).i; + yaxisID = (tagint) ubuf(xyzaxis[i][1]).i; + zaxisID = (tagint) ubuf(xyzaxis[i][2]).i; if (xaxisID) { index = atom->map(xaxisID); @@ -742,7 +738,7 @@ void PairAmoeba::find_multipole_neighbors() } else xaxis2local[i] = i; if (yaxisID) { - if (yaxis[i] < 0) yaxisID = -yaxisID; + if (xyzaxis[i][1] < 0) yaxisID = -yaxisID; index = atom->map(yaxisID); if (index == -1) nmissing++; else { diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 0c9f6f8b6e..8bae12955a 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -316,7 +316,7 @@ void PairAmoeba::compute(int eflag, int vflag) // if reneighboring step: // augment neighbor list to include 1-5 neighbor flags - // re-create xyz axis2local and red2local + // re-create red2local and xyz axis2local // re-create induce neighbor list if (neighbor->ago == 0) { @@ -710,15 +710,9 @@ void PairAmoeba::init_style() index_redID = atom->find_custom("redID",flag,cols); if (index_redID < 0 || !flag || cols) error->all(FLERR,"Pair amoeba redID is not defined"); - index_xaxis = atom->find_custom("xaxis",flag,cols); - if (index_xaxis < 0 || !flag || cols) - error->all(FLERR,"Pair amoeba xaxis is not defined"); - index_yaxis = atom->find_custom("yaxis",flag,cols); - if (index_yaxis < 0 || !flag || cols) - error->all(FLERR,"Pair amoeba yaxis is not defined"); - index_zaxis = atom->find_custom("zaxis",flag,cols); - if (index_zaxis < 0 || !flag || cols) - error->all(FLERR,"Pair amoeba zaxis is not defined"); + index_xyzaxis = atom->find_custom("xyzaxis",flag,cols); + if (index_xyzaxis < 0 || !flag || cols == 0) + error->all(FLERR,"Pair amoeba xyzaxis is not defined"); index_polaxe = atom->find_custom("polaxe",flag,cols); if (index_polaxe < 0 || flag || cols) @@ -1995,11 +1989,8 @@ void PairAmoeba::mix() double TWOSIX = pow(2.0,1.0/6.0); for (i = 1; i <= n_amclass; i++) { - - // printf("MIX i %d nclass %d eps %g sigma %g\n", - // i,n_amclass,vdwl_eps[i],vdwl_sigma[i]); - for (j = i; j <= n_amclass; j++) { + ei = vdwl_eps[i]; ej = vdwl_eps[j]; ri = vdwl_sigma[i]; @@ -2099,21 +2090,43 @@ void *PairAmoeba::extract(const char *str, int &dim) /* ---------------------------------------------------------------------- grow local vectors and arrays if necessary - keep them atom->nmax in length - NOTE: some of these do not need to grow unless nlocal > atom->nmax - these are ones that never store ghost values - could realloc them separately - e.g. thetai,igrid,fopt + some need space for ghost atoms, most do not ------------------------------------------------------------------------- */ void PairAmoeba::grow_local() { + // only reallocate if nlocal > nmax + + if (atom->nlocal > nmax) { + } + + + + // forward/reverse comm, so ghost values + uind; + uinp; + uopt; + uoptp; + rsd; + rsdp; + xred; // AMOEBA + + memory->destroy(rsd); + memory->destroy(rsdp); + memory->create(rsd,nmax,3,"amoeba:rsd"); + memory->create(rsdp,nmax,3,"amoeba:rsdp"); + + + memory->destroy(rpole); + memory->create(rpole,nmax,13,"amoeba:rpole"); + + + // free vectors and arrays memory->destroy(xaxis2local); memory->destroy(yaxis2local); memory->destroy(zaxis2local); - memory->destroy(rpole); memory->destroy(tq); if (amoeba) { @@ -2134,8 +2147,6 @@ void PairAmoeba::grow_local() memory->destroy(fieldp); memory->destroy(ufld); memory->destroy(dufld); - memory->destroy(rsd); - memory->destroy(rsdp); memory->destroy(zrsd); memory->destroy(zrsdp); @@ -2164,7 +2175,6 @@ void PairAmoeba::grow_local() memory->create(xaxis2local,nmax,"amoeba:xaxis2local"); memory->create(yaxis2local,nmax,"amoeba:yaxis2local"); memory->create(zaxis2local,nmax,"amoeba:zaxis2local"); - memory->create(rpole,nmax,13,"amoeba:rpole"); memory->create(tq,nmax,3,"amoeba:tq"); if (amoeba) { @@ -2185,8 +2195,6 @@ void PairAmoeba::grow_local() memory->create(fieldp,nmax,3,"amoeba:fieldp"); memory->create(ufld,nmax,3,"amoeba:ufld"); memory->create(dufld,nmax,6,"amoeba:dufld"); - memory->create(rsd,nmax,3,"amoeba:rsd"); - memory->create(rsdp,nmax,3,"amoeba:rsdp"); memory->create(zrsd,nmax,3,"amoeba:zrsd"); memory->create(zrsdp,nmax,3,"amoeba:zrsdp"); @@ -2212,96 +2220,3 @@ void PairAmoeba::grow_local() } } -// ---------------------------------------------------------------------- -// debug output methods -// ---------------------------------------------------------------------- - -/* ---------------------------------------------------------------------- - dump ID + 6 values from 2 (N,3) per-atom arrays - only proc 0 can write - file is already open -------------------------------------------------------------------------- */ - -void PairAmoeba::dump6(FILE *fp, const char *columns, double scale, - double **a, double **b) -{ - int i,j,m; - MPI_Status status; - MPI_Request request; - - // setup - - int size_one = 7; - int nlocal = atom->nlocal; - - char boundstr[9]; // encoding of boundary flags - domain->boundary_string(boundstr); - - int maxlocal; - MPI_Allreduce(&nlocal,&maxlocal,1,MPI_INT,MPI_MAX,world); - - double *buf; - memory->create(buf,maxlocal*size_one,"amoeba:buf"); - - // pack my data - - tagint *tag = atom->tag; - - m = 0; - for (i = 0; i < nlocal; i++) { - buf[m++] = tag[i]; - buf[m++] = scale*a[i][0]; - buf[m++] = scale*a[i][1]; - buf[m++] = scale*a[i][2]; - buf[m++] = scale*b[i][0]; - buf[m++] = scale*b[i][1]; - buf[m++] = scale*b[i][2]; - } - - // write file - - if (me == 0) { - fprintf(fp,"ITEM: TIMESTEP\n"); - fprintf(fp,BIGINT_FORMAT "\n",update->ntimestep); - fprintf(fp,"ITEM: NUMBER OF ATOMS\n"); - fprintf(fp,BIGINT_FORMAT "\n",atom->natoms); - fprintf(fp,"ITEM: BOX BOUNDS %s\n",boundstr); - fprintf(fp,"%-1.16e %-1.16e\n",domain->boxlo[0],domain->boxhi[0]); - fprintf(fp,"%-1.16e %-1.16e\n",domain->boxlo[1],domain->boxhi[1]); - fprintf(fp,"%-1.16e %-1.16e\n",domain->boxlo[2],domain->boxhi[2]); - fprintf(fp,"ITEM: ATOMS %s\n",columns); - } - - int nlines; - double tmp; - - if (me == 0) { - for (int iproc = 0; iproc < nprocs; iproc++) { - if (iproc) { - MPI_Irecv(buf,maxlocal*size_one,MPI_DOUBLE,iproc,0,world,&request); - MPI_Send(&tmp,0,MPI_INT,iproc,0,world); - MPI_Wait(&request,&status); - MPI_Get_count(&status,MPI_DOUBLE,&nlines); - nlines /= size_one; - } else nlines = nlocal; - - m = 0; - for (i = 0; i < nlines; i++) { - for (j = 0; j < size_one; j++) { - if (j == 0) fprintf(fp,"%d",static_cast (buf[m])); - else fprintf(fp," %g",buf[m]); - m++; - } - fprintf(fp,"\n"); - } - } - - } else { - MPI_Recv(&tmp,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE); - MPI_Rsend(buf,nlocal*size_one,MPI_DOUBLE,0,0,world); - } - - // clean up - - memory->destroy(buf); -} diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h index 0057219b98..5de382e2d9 100644 --- a/src/AMOEBA/pair_amoeba.h +++ b/src/AMOEBA/pair_amoeba.h @@ -159,7 +159,7 @@ class PairAmoeba : public Pair { // static per-atom properties, must persist as atoms migrate int index_amtype,index_amgroup,index_redID; - int index_xaxis,index_yaxis,index_zaxis,index_polaxe,index_pval; + int index_xyzaxis,index_polaxe,index_pval; int *amtype; // AMOEBA type, 1 to N_amtype int *amgroup; // AMOEBA polarization group, 1 to Ngroup