tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refactoring edits

Browse Source
This commit is contained in:
Steve Plimpton
2022-05-17 14:40:09 -06:00
parent 510e78d4d3
commit e07b46c771
6 changed files with 89 additions and 408 deletions
Download Patch File Download Diff File Expand all files Collapse all files

56
examples/amoeba/Compare.sh
View File

@ -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

215
src/AMOEBA/amoeba_induce.cpp
View File

@ -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;

17
src/AMOEBA/amoeba_polar.cpp
View File

@ -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);

56
src/AMOEBA/amoeba_utils.cpp
View File

@ -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 {

151
src/AMOEBA/pair_amoeba.cpp
View File

@ -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<tagint> (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);
}

2
src/AMOEBA/pair_amoeba.h
View File

@ -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