ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Cleaned up commented-out and debugging stuffs, removed irrelevant changes to lj/cut/dipole/cut, reverted unwanted changes in the PPPMGPU destructor, fixed unresolved conflicts in tinker.py, updated the userbinsize==0 case in atom.cpp and using Force::pair_match() as suggested. Internal timing stuffs need work.

Browse Source
This commit is contained in:
Trung Nguyen
2023-01-15 15:41:54 -06:00
parent c21f2faa1f
commit 67574601ed
12 changed files with 19 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
cmake/CMakeLists.txt
View File

@ -395,7 +395,6 @@ endif()
pkg_depends(ML-IAP ML-SNAP)
pkg_depends(MPIIO MPI)
pkg_depends(ATC MANYBODY)
pkg_depends(AMOEBA KSPACE)
pkg_depends(LATBOLTZ MPI)
pkg_depends(SCAFACOS MPI)
pkg_depends(AMOEBA KSPACE)

2
examples/amoeba/in.ubiquitin
View File

@ -34,8 +34,6 @@ pair_coeff * * amoeba_ubiquitin.prm amoeba_ubiquitin.key
special_bonds lj/coul 0.5 0.5 0.5 one/five yes
# setup force components this way so can dump them (AMOEBA or HIPPO also needs them for now)
fix fhal all store/state 0 fx fy fz
fix frepulse all store/state 0 fx fy fz
fix fdisp all store/state 0 fx fy fz

40
src/AMOEBA/amoeba_induce.cpp
View File

@ -86,17 +86,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++) {
@ -213,16 +202,7 @@ void PairAmoeba::induce()
cfstyle = INDUCE;
comm->forward_comm(this);
/*
if (comm->me == 0) {
printf("CPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < 20; i++) {
printf("i = %d: uind = %f %f %f; udirp = %f %f %f\n",
i, uind[i][0], uind[i][1], uind[i][2],
uinp[i][0], uinp[i][1], uinp[i][2]);
}
}
*/
ufield0c(field,fieldp);
crstyle = FIELD;
@ -285,18 +265,6 @@ void PairAmoeba::induce()
crstyle = FIELD;
comm->reverse_comm(this);
//error->all(FLERR,"STOP");
/*
if (comm->me == 0) {
printf("CPU: iter = %d\n", iter);
for (i = 0; i < 10; i++) {
printf("i = %d: field = %f %f %f; fieldp = %f %f %f\n",
i, 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];
@ -413,8 +381,6 @@ void PairAmoeba::induce()
}
}
// if (comm->me == 0) printf("CG iteration count = %d\n",iter);
// terminate the calculation if dipoles failed to converge
// NOTE: could make this an error
@ -1033,9 +999,7 @@ void PairAmoeba::umutual2b(double **field, double **fieldp)
j = jlist[jj];
uindj = uind[j];
uinpj = uinp[j];
//if (i==0 && j == 10)
// printf("i = %d: j = %d: tdipdip %f %f %f %f %f %f\n",
// i, j,tdipdip[0],tdipdip[1],tdipdip[2],tdipdip[3],tdipdip[4],tdipdip[5]);
fid[0] = tdipdip[0]*uindj[0] + tdipdip[1]*uindj[1] + tdipdip[2]*uindj[2];
fid[1] = tdipdip[1]*uindj[0] + tdipdip[3]*uindj[1] + tdipdip[4]*uindj[2];
fid[2] = tdipdip[2]*uindj[0] + tdipdip[4]*uindj[1] + tdipdip[5]*uindj[2];

8
src/AMOEBA/amoeba_kspace.cpp
View File

@ -68,8 +68,6 @@ void PairAmoeba::moduli()
int maxfft = MAX(nfft1,nfft2);
maxfft = MAX(maxfft,nfft3);
//double *array = new double[bsorder];
//double *bsarray = new double[maxfft];
if (maxfft > _nfft_max) {
memory->destroy(_moduli_bsarray);
_nfft_max = maxfft;
@ -79,7 +77,6 @@ void PairAmoeba::moduli()
// compute and load the moduli values
double x = 0.0;
//bspline(x,bsorder,array);
bspline(x,bsorder,_moduli_array);
for (i = 0; i < maxfft; i++) _moduli_bsarray[i] = 0.0;
@ -88,11 +85,6 @@ void PairAmoeba::moduli()
dftmod(bsmod1,_moduli_bsarray,nfft1,bsorder);
dftmod(bsmod2,_moduli_bsarray,nfft2,bsorder);
dftmod(bsmod3,_moduli_bsarray,nfft3,bsorder);
// perform deallocation of local arrays
//delete[] array;
//delete[] bsarray;
}
/* ----------------------------------------------------------------------

6
src/AMOEBA/fix_amoeba_bitorsion.cpp
View File

@ -194,10 +194,8 @@ void FixAmoebaBiTorsion::init()
// error check that PairAmoeba or PairHiippo exist
pair = nullptr;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("amoeba/gpu",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) pair = force->pair_match("hippo/gpu",1,0);
pair = force->pair_match("^amoeba",0,0);
if (!pair) pair = force->pair_match("^hippo",0,0);
if (!pair)
error->all(FLERR,"Cannot use fix amoeba/bitorsion w/out pair amoeba/hippo");

7
src/AMOEBA/improper_amoeba.cpp
View File

@ -285,10 +285,9 @@ void ImproperAmoeba::init_style()
// check if PairAmoeba disabled improper terms
Pair *pair = nullptr;
pair = force->pair_match("amoeba",1,0);
if (!pair) pair = force->pair_match("amoeba/gpu",1,0);
if (!pair) pair = force->pair_match("hippo",1,0);
if (!pair) pair = force->pair_match("hippo/gpu",1,0);
pair = force->pair_match("^amoeba",0,0);
if (!pair) pair = force->pair_match("^hippo",0,0);
if (!pair) error->all(FLERR,"Improper amoeba could not find pair amoeba/hippo");
int tmp;

7
src/AMOEBA/pair_amoeba.cpp
View File

@ -1055,13 +1055,6 @@ void PairAmoeba::init_style()
// request standard neighbor list
// int irequest = neighbor->request(this,instance_me);
// for DEBUGGING with GPU
//neighbor->requests[irequest]->half = 0;
//neighbor->requests[irequest]->full = 1;
neighbor->add_request(this);
}

19
src/DIPOLE/pair_lj_cut_dipole_cut.cpp
View File

@ -90,8 +90,6 @@ void PairLJCutDipoleCut::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
int maxsize = 10;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
@ -104,13 +102,6 @@ void PairLJCutDipoleCut::compute(int eflag, int vflag)
jlist = firstneigh[i];
jnum = numneigh[i];
double scale_dipole = 1.0;
if (jnum > maxsize) {
scale_dipole = maxsize; //1.0/(double)maxsize;
} else {
scale_dipole = jnum; //1.0/(double)jnum;
}
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
@ -216,7 +207,7 @@ void PairLJCutDipoleCut::compute(int eflag, int vflag)
// total force
fq = scale_dipole*factor_coul*qqrd2e;
fq = factor_coul*qqrd2e;
fx = fq*forcecoulx + delx*forcelj;
fy = fq*forcecouly + dely*forcelj;
fz = fq*forcecoulz + delz*forcelj;
@ -230,7 +221,7 @@ void PairLJCutDipoleCut::compute(int eflag, int vflag)
torque[i][1] += fq*tiycoul;
torque[i][2] += fq*tizcoul;
if (newton_pair) {
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
@ -371,13 +362,7 @@ void PairLJCutDipoleCut::init_style()
if (!atom->q_flag || !atom->mu_flag || !atom->torque_flag)
error->all(FLERR,"Pair dipole/cut requires atom attributes q, mu, torque");
<<<<<<< HEAD
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
=======
neighbor->add_request(this);
>>>>>>> amoeba
}
/* ----------------------------------------------------------------------

4
src/Depend.sh
View File

@ -45,6 +45,10 @@ depend () {
# add one if statement per parent package
# add one depend() call per child package that depends on that parent
if (test $1 = "AMOEBA") then
depend GPU
fi
if (test $1 = "ASPHERE") then
depend GPU
depend OPENMP

2
src/GPU/pppm_gpu.cpp
View File

@ -102,6 +102,8 @@ PPPMGPU::PPPMGPU(LAMMPS *lmp) : PPPM(lmp)
PPPMGPU::~PPPMGPU()
{
PPPM_GPU_API(clear)(poisson_time);
destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out);
destroy_3d_offset(vd_brick,nzlo_out,nylo_out);
}
/* ----------------------------------------------------------------------

9
src/atom.cpp
View File

@ -2359,15 +2359,12 @@ void Atom::setup_sort_bins()
#ifdef LMP_GPU
if (userbinsize == 0.0) {
int ifix = modify->find_fix("package_gpu");
if (ifix >= 0) {
auto ifix = dynamic_cast<FixGPU *>(modify->get_fix_by_id("package_gpu"));
if (ifix) {
const double subx = domain->subhi[0] - domain->sublo[0];
const double suby = domain->subhi[1] - domain->sublo[1];
const double subz = domain->subhi[2] - domain->sublo[2];
FixGPU *fix = static_cast<FixGPU *>(modify->fix[ifix]);
binsize = fix->binsize(subx, suby, subz, atom->nlocal,
0.5 * neighbor->cutneighmax);
binsize = ifix->binsize(subx, suby, subz, atom->nlocal, 0.5 * neighbor->cutneighmax);
}
}
#endif

44
tools/tinker/tinker2lmp.py
View File

@ -227,11 +227,7 @@ class XYZfile(object):
print(i+1,label[i],x[i],y[i],z[i],type[i], end=' ', file=fp)
for j in bonds[i]: print(j, end=' ', file=fp)
print(file=fp)
<<<<<<< HEAD
=======
>>>>>>> develop
fp.close()
# triplet of atoms in an angle = atom 1,2,3
@ -1098,16 +1094,6 @@ for i,one in enumerate(alist):
elif len(params[3]) == 2:
nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
<<<<<<< HEAD
if nbonds != 3:
print("Center angle atom has wrong bond count")
print(" angle atom IDs:",atom1,atom2,atom3)
print(" angle atom classes:",c1,c2,c3)
print(" Tinker FF file param options:",len(params[3]))
print(" Nbonds and hydrogen count:",nbonds,hcount)
#sys.exit() NOTE: allow this for now
=======
#if nbonds != 3:
#print("Center angle atom has wrong bond count")
@ -1117,33 +1103,12 @@ for i,one in enumerate(alist):
#print(" Nbonds and hydrogen count:",nbonds,hcount)
# NOTE: allow this for now
#sys.exit()
>>>>>>> develop
if hcount == 0: which = 1
elif hcount == 1:
which = 2
m += 1
<<<<<<< HEAD
print("3-bond angle")
print(" angle atom IDs:",atom1,atom2,atom3)
print(" angle atom classes:",c1,c2,c3)
print(" Tinker FF file param options:",len(params[3]))
print(" Nbonds and hydrogen count:",nbonds,hcount)
print(" which:",which,m)
elif len(params[3]) == 3:
nbonds,hcount = xyz.angle_hbond_count(atom1,atom2,atom3,lmptype,lmpmass)
if nbonds != 4:
print("Center angle atom has wrong bond count")
print(" angle atom IDs:",atom1,atom2,atom3)
print(" angle atom classes:",c1,c2,c3)
print(" Tinker FF file param options:",len(params[3]))
print(" Nbonds and hydrogen count:",nbonds,hcount)
#sys.exit() NOTE: allow this for now
=======
#print("3-bond angle")
#print(" angle atom IDs:",atom1,atom2,atom3)
#print(" angle atom classes:",c1,c2,c3)
@ -1163,7 +1128,6 @@ for i,one in enumerate(alist):
# NOTE: allow this for now
#sys.exit()
>>>>>>> develop
if hcount == 0: which = 1
elif hcount == 1:
which = 2
@ -1207,12 +1171,8 @@ for itype in range(len(aparams)):
elif (c3,c2,c1) in badict:
n1,n2,r1,r2 = badict[(c3,c2,c1)]
else:
<<<<<<< HEAD
print("Bond-stretch angle triplet not found: %d %d %d" % (c1,c2,c3))
=======
# NOTE: just for debugging
#print("Bond-stretch angle triplet not found: %d %d %d" % (c1,c2,c3))
>>>>>>> develop
n1,n2,r1,r2 = 4*[0.0]
baparams.append((n1,n2,r1,r2))
@ -1670,11 +1630,7 @@ print("Natoms =",natoms)
print("Ntypes =",ntypes)
print("Tinker XYZ types =",len(tink2lmp))
print("Tinker PRM types =",prm.ntypes)
<<<<<<< HEAD
#print "Tinker groups =",ngroups
=======
#print("Tinker groups =",ngroups)
>>>>>>> develop
print("Nmol =",nmol)
print("Nbonds =",nbonds)
print("Nangles =",nangles)