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

Made dfield0c work to compute uind and uinp correctly; need to make sure they are correct for polar_real()

Browse Source
This commit is contained in:
Trung Nguyen
2021-09-08 16:43:33 -05:00
parent 1c5d235f12
commit 8c5a116d30
3 changed files with 117 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
lib/gpu/lal_base_amoeba.cpp
View File

@ -528,14 +528,14 @@ int** BaseAmoebaT::compute_udirect2b(const int ago, const int inum_full, const i
// copy field and fieldp from device to host (_fieldp store both arrays, one after another)
_fieldp.update_host(_max_fieldp_size*8,false);
/*
printf("GPU lib: _fieldp size = %d: max fieldp size = %d\n",
this->_fieldp.cols(), _max_fieldp_size);
for (int i = 0; i < 10; i++) {
numtyp4* p = (numtyp4*)(&this->_fieldp[4*i]);
printf("i = %d; field = %f %f %f\n", i, p->x, p->y, p->z);
}
*/
return firstneigh; //nbor->host_jlist.begin()-host_start;
}

63
src/AMOEBA/amoeba_induce.cpp
View File

@ -25,6 +25,7 @@
#include "my_page.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "error.h"
using namespace LAMMPS_NS;
@ -103,14 +104,21 @@ void PairAmoeba::induce()
memory->create(usump,nlocal,3,"ameoba/induce:usump");
// get the electrostatic field due to permanent multipoles
dfield0c(field,fieldp);
// reverse comm to sum field,fieldp from ghost atoms to owned atoms
crstyle = FIELD;
comm->reverse_comm_pair(this);
/*
printf("CPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < nlocal; 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]);
}
*/
// DEBUG statements
/*
@ -135,7 +143,14 @@ void PairAmoeba::induce()
}
}
}
/*
printf("CPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < 10; i++) {
printf("i = %d: udir = %f %f %f; udirp = %f %f %f\n",
i, udir[i][0], udir[i][1], udir[i][2],
udirp[i][0], udirp[i][1], udirp[i][2]);
}
*/
// DEBUG statements
/*
@ -250,12 +265,30 @@ void PairAmoeba::induce()
cfstyle = INDUCE;
comm->forward_comm_pair(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;
comm->reverse_comm_pair(this);
/*
if (comm->me == 0) {
printf("CPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < nlocal; 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]);
}
}
*/
// DEBUG statements
/*
@ -342,6 +375,16 @@ void PairAmoeba::induce()
crstyle = FIELD;
comm->reverse_comm_pair(this);
/*
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]);
}
}
*/
// DEBUG statements
@ -537,6 +580,7 @@ void PairAmoeba::induce()
error->warning(FLERR,"AMOEBA induced dipoles did not converge");
}
// DEBUG output to dump file
if (uind_flag)
@ -553,6 +597,15 @@ void PairAmoeba::induce()
memory->destroy(usum);
memory->destroy(usump);
if (comm->me == 0) {
printf("CPU: iter = %d\n", iter);
for (i = 0; i < 20; i++) {
printf("i = %d: uind = %f %f %f; uinp = %f %f %f\n",
i, uind[i][0], uind[i][1], uind[i][2],
uinp[i][0], uinp[i][1], uinp[i][2]);
}
}
// update the lists of previous induced dipole values
// shift previous m values up to m+1, add new values at m = 0
// only when preconditioner is used

67
src/GPU/pair_amoeba_gpu.cpp
View File

@ -260,7 +260,7 @@ void PairAmoebaGPU::init_style()
}
}
// select the cutoff (off2) for neighbor list builds (the polar term for now)
// select the squared cutoff (off2) for neighbor list builds (the polar term for now)
// NOTE: induce and polar terms are using the same flags here
if (use_ewald) choose(POLAR_LONG);
@ -365,13 +365,10 @@ void PairAmoebaGPU::induce()
// get the electrostatic field due to permanent multipoles
dfield0c(field,fieldp);
// reverse comm to sum field,fieldp from ghost atoms to owned atoms
/*
crstyle = FIELD;
comm->reverse_comm_pair(this);
*/
// set induced dipoles to polarizability times direct field
for (i = 0; i < nlocal; i++) {
@ -385,7 +382,14 @@ void PairAmoebaGPU::induce()
}
}
}
/*
printf("GPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < 10; i++) {
printf("i = %d: udir = %f %f %f; udirp = %f %f %f\n",
i, udir[i][0], udir[i][1], udir[i][2],
udirp[i][0], udirp[i][1], 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 ??
@ -489,13 +493,30 @@ void PairAmoebaGPU::induce()
cfstyle = INDUCE;
comm->forward_comm_pair(this);
/*
if (comm->me == 0) {
printf("GPU: 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;
comm->reverse_comm_pair(this);
/*
if (comm->me == 0) {
printf("GPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < nlocal; 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]);
}
}
*/
// set initial conjugate gradient residual and conjugate vector
for (i = 0; i < nlocal; i++) {
@ -554,7 +575,16 @@ void PairAmoebaGPU::induce()
crstyle = FIELD;
comm->reverse_comm_pair(this);
/*
if (comm->me == 0) {
printf("GPU: 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];
@ -697,6 +727,15 @@ void PairAmoebaGPU::induce()
memory->destroy(usum);
memory->destroy(usump);
if (comm->me == 0) {
printf("GPU: iter = %d\n", iter);
for (i = 0; i < 20; i++) {
printf("i = %d: uind = %f %f %f; uinp = %f %f %f\n",
i, uind[i][0], uind[i][1], uind[i][2],
uinp[i][0], uinp[i][1], uinp[i][2]);
}
}
// update the lists of previous induced dipole values
// shift previous m values up to m+1, add new values at m = 0
// only when preconditioner is used
@ -758,7 +797,7 @@ void PairAmoebaGPU::dfield0c(double **field, double **fieldp)
// get the real space portion of the permanent field
if (rspace_flag) udirect2b(field,fieldp);
// get the self-energy portion of the permanent field
term = (4.0/3.0) * aewald*aewald*aewald / MY_PIS;
@ -768,6 +807,14 @@ void PairAmoebaGPU::dfield0c(double **field, double **fieldp)
fieldp[i][j] += term*rpole[i][j+1];
}
}
/*
printf("GPU: cutghost = %f\n", comm->cutghost[0]);
for (i = 0; i < nlocal; 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]);
}
*/
}
/* ----------------------------------------------------------------------