forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8094 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
95cf2cb86f
commit
9eeb694432
|
|
@ -14,6 +14,7 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
Contributing authors: Roy Pollock (LLNL), Paul Crozier (SNL)
|
||||
per-atom energy/virial added by German Samolyuk (ORNL), Stan Moore (BYU)
|
||||
group/group energy/force added by Stan Moore (BYU)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "mpi.h"
|
||||
|
|
@ -42,6 +43,8 @@ Ewald::Ewald(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
|
|||
{
|
||||
if (narg != 1) error->all(FLERR,"Illegal kspace_style ewald command");
|
||||
|
||||
group_group_enable = 1;
|
||||
|
||||
accuracy_relative = atof(arg[0]);
|
||||
|
||||
kmax = 0;
|
||||
|
|
@ -64,6 +67,7 @@ Ewald::Ewald(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
|
|||
Ewald::~Ewald()
|
||||
{
|
||||
deallocate();
|
||||
deallocate_groups();
|
||||
memory->destroy(ek);
|
||||
memory->destroy3d_offset(cs,-kmax_created);
|
||||
memory->destroy3d_offset(sn,-kmax_created);
|
||||
|
|
@ -261,6 +265,7 @@ void Ewald::setup()
|
|||
if (kmax > kmax_old) {
|
||||
deallocate();
|
||||
allocate();
|
||||
group_allocate_flag = 0;
|
||||
|
||||
memory->destroy(ek);
|
||||
memory->destroy3d_offset(cs,-kmax_created);
|
||||
|
|
@ -966,3 +971,199 @@ double Ewald::memory_usage()
|
|||
bytes += 2 * (2*kmax+1)*3*nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
group-group interactions
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute the Ewald total long-range force and energy for groups A and B
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Ewald::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
|
||||
{
|
||||
if (slabflag)
|
||||
error->all(FLERR,"Cannot (yet) use Kspace slab correction "
|
||||
"with compute group/group");
|
||||
|
||||
int i,k;
|
||||
|
||||
if (!group_allocate_flag) {
|
||||
allocate_groups();
|
||||
group_allocate_flag = 1;
|
||||
}
|
||||
|
||||
e2group = 0; //energy
|
||||
f2group[0] = 0; //force in x-direction
|
||||
f2group[1] = 0; //force in y-direction
|
||||
f2group[2] = 0; //force in z-direction
|
||||
|
||||
// partial and total structure factors for groups A and B
|
||||
|
||||
for (k = 0; k < kcount; k++) {
|
||||
|
||||
// group A
|
||||
|
||||
sfacrl_A[k] = 0;
|
||||
sfacim_A[k] = 0;
|
||||
sfacrl_A_all[k] = 0;
|
||||
sfacim_A_all[k] = 0;
|
||||
|
||||
// group B
|
||||
|
||||
sfacrl_B[k] = 0;
|
||||
sfacim_B[k] = 0;
|
||||
sfacrl_B_all[k] = 0;
|
||||
sfacim_B_all[k] = 0;
|
||||
}
|
||||
|
||||
double *q = atom->q;
|
||||
int nlocal = atom->nlocal;
|
||||
int *mask = atom->mask;
|
||||
|
||||
int kx,ky,kz;
|
||||
double cypz,sypz,exprl,expim;
|
||||
|
||||
// partial structure factors for groups A and B on each processor
|
||||
|
||||
for (k = 0; k < kcount; k++) {
|
||||
kx = kxvecs[k];
|
||||
ky = kyvecs[k];
|
||||
kz = kzvecs[k];
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
|
||||
if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B))
|
||||
if (BA_flag) continue;
|
||||
|
||||
if ((mask[i] & groupbit_A) || (mask[i] & groupbit_B)) {
|
||||
|
||||
cypz = cs[ky][1][i]*cs[kz][2][i] - sn[ky][1][i]*sn[kz][2][i];
|
||||
sypz = sn[ky][1][i]*cs[kz][2][i] + cs[ky][1][i]*sn[kz][2][i];
|
||||
exprl = cs[kx][0][i]*cypz - sn[kx][0][i]*sypz;
|
||||
expim = sn[kx][0][i]*cypz + cs[kx][0][i]*sypz;
|
||||
|
||||
// group A
|
||||
|
||||
if (mask[i] & groupbit_A) {
|
||||
sfacrl_A[k] += q[i]*exprl;
|
||||
sfacim_A[k] += q[i]*expim;
|
||||
}
|
||||
|
||||
// group B
|
||||
|
||||
if (mask[i] & groupbit_B) {
|
||||
sfacrl_B[k] += q[i]*exprl;
|
||||
sfacim_B[k] += q[i]*expim;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// total structure factor by summing over procs
|
||||
|
||||
MPI_Allreduce(sfacrl_A,sfacrl_A_all,kcount,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(sfacim_A,sfacim_A_all,kcount,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
MPI_Allreduce(sfacrl_B,sfacrl_B_all,kcount,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(sfacim_B,sfacim_B_all,kcount,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
double partial_group;
|
||||
|
||||
|
||||
// total group A <--> group B energy
|
||||
|
||||
for (k = 0; k < kcount; k++) {
|
||||
partial_group = sfacrl_A_all[k]*sfacrl_B_all[k] +
|
||||
sfacim_A_all[k]*sfacim_B_all[k];
|
||||
e2group += ug[k]*partial_group;
|
||||
}
|
||||
|
||||
// total charge of groups A & B, needed for self-energy correction
|
||||
|
||||
double qsum_group = 0.0, qsqsum_group = 0.0;
|
||||
|
||||
for (int i = 0; i < atom->nlocal; i++) {
|
||||
if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B)) {
|
||||
if (BA_flag) continue;
|
||||
qsqsum_group += q[i]*q[i];
|
||||
}
|
||||
|
||||
//if ((mask[i] & groupbit_A) || (mask[i] & groupbit_B))
|
||||
// qsum_group += q[i];
|
||||
}
|
||||
|
||||
double tmp;
|
||||
MPI_Allreduce(&qsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
qsum_group = tmp;
|
||||
|
||||
//MPI_Allreduce(&qsqsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
//qsqsum_group = tmp;
|
||||
|
||||
// self-energy correction
|
||||
|
||||
e2group -= g_ewald*qsqsum_group/MY_PIS;
|
||||
|
||||
// should we include the boundary condition or not?
|
||||
// e2group += MY_PI2*qsum_group*qsum_group / (g_ewald*g_ewald*volume);
|
||||
|
||||
e2group *= qscale;
|
||||
|
||||
|
||||
// total group A <--> group B force
|
||||
|
||||
for (k = 0; k < kcount; k++) {
|
||||
partial_group = sfacim_A_all[k]*sfacrl_B_all[k] -
|
||||
sfacrl_A_all[k]*sfacim_B_all[k];
|
||||
f2group[0] += eg[k][0]*partial_group;
|
||||
f2group[1] += eg[k][1]*partial_group;
|
||||
f2group[2] += eg[k][2]*partial_group;
|
||||
}
|
||||
|
||||
f2group[0] *= qscale;
|
||||
f2group[1] *= qscale;
|
||||
f2group[2] *= qscale;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate group-group memory that depends on # of K-vectors
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Ewald::allocate_groups()
|
||||
{
|
||||
// group A
|
||||
|
||||
sfacrl_A = new double[kmax3d];
|
||||
sfacim_A = new double[kmax3d];
|
||||
sfacrl_A_all = new double[kmax3d];
|
||||
sfacim_A_all = new double[kmax3d];
|
||||
|
||||
// group B
|
||||
|
||||
sfacrl_B = new double[kmax3d];
|
||||
sfacim_B = new double[kmax3d];
|
||||
sfacrl_B_all = new double[kmax3d];
|
||||
sfacim_B_all = new double[kmax3d];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
deallocate group-group memory that depends on # of K-vectors
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void Ewald::deallocate_groups()
|
||||
{
|
||||
// group A
|
||||
|
||||
delete [] sfacrl_A;
|
||||
delete [] sfacim_A;
|
||||
delete [] sfacrl_A_all;
|
||||
delete [] sfacim_A_all;
|
||||
|
||||
// group B
|
||||
|
||||
delete [] sfacrl_B;
|
||||
delete [] sfacim_B;
|
||||
delete [] sfacrl_B_all;
|
||||
delete [] sfacim_B_all;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -33,6 +33,8 @@ class Ewald : public KSpace {
|
|||
virtual void compute(int, int);
|
||||
double memory_usage();
|
||||
|
||||
void compute_group_group(int, int, int);
|
||||
|
||||
protected:
|
||||
int kxmax,kymax,kzmax;
|
||||
int kcount,kmax,kmax3d,kmax_created;
|
||||
|
|
@ -47,12 +49,23 @@ class Ewald : public KSpace {
|
|||
double *sfacrl,*sfacim,*sfacrl_all,*sfacim_all;
|
||||
double ***cs,***sn;
|
||||
|
||||
// group-group interactions
|
||||
|
||||
int group_allocate_flag;
|
||||
double *sfacrl_A,*sfacim_A,*sfacrl_A_all,*sfacim_A_all;
|
||||
double *sfacrl_B,*sfacim_B,*sfacrl_B_all,*sfacim_B_all;
|
||||
|
||||
double rms(int, double, bigint, double);
|
||||
virtual void eik_dot_r();
|
||||
void coeffs();
|
||||
virtual void allocate();
|
||||
void deallocate();
|
||||
void slabcorr();
|
||||
|
||||
// group-group interactions
|
||||
|
||||
void allocate_groups();
|
||||
void deallocate_groups();
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@
|
|||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing authors: Roy Pollock (LLNL), Paul Crozier (SNL)
|
||||
per-atom energy/virial added by Stan Moore (BYU)
|
||||
per-atom energy/virial, group/group energy/force added by Stan Moore (BYU)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "lmptype.h"
|
||||
|
|
@ -60,6 +60,8 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
|
|||
{
|
||||
if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command");
|
||||
|
||||
group_group_enable = 1;
|
||||
|
||||
accuracy_relative = atof(arg[0]);
|
||||
|
||||
nfactors = 3;
|
||||
|
|
@ -81,6 +83,9 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
|
|||
fkx = fky = fkz = NULL;
|
||||
buf1 = buf2 = buf3 = buf4 = NULL;
|
||||
|
||||
density_A_brick = density_B_brick = NULL;
|
||||
density_A_fft = density_B_fft = NULL;
|
||||
|
||||
gf_b = NULL;
|
||||
rho1d = rho_coeff = NULL;
|
||||
|
||||
|
|
@ -100,6 +105,7 @@ PPPM::~PPPM()
|
|||
delete [] factors;
|
||||
deallocate();
|
||||
deallocate_peratom();
|
||||
deallocate_groups();
|
||||
memory->destroy(part2grid);
|
||||
}
|
||||
|
||||
|
|
@ -142,6 +148,8 @@ void PPPM::init()
|
|||
deallocate();
|
||||
deallocate_peratom();
|
||||
peratom_allocate_flag = 0;
|
||||
deallocate_groups();
|
||||
group_allocate_flag = 0;
|
||||
|
||||
// extract short-range Coulombic cutoff from pair style
|
||||
|
||||
|
|
@ -2531,5 +2539,331 @@ double PPPM::memory_usage()
|
|||
bytes += 2 * nbuf_peratom * sizeof(FFT_SCALAR);
|
||||
}
|
||||
|
||||
if (group_allocate_flag) {
|
||||
bytes += 2 * nbrick * sizeof(FFT_SCALAR);
|
||||
bytes += 2 * nfft_both * sizeof(FFT_SCALAR);;
|
||||
}
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
group-group interactions
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute the PPPM total long-range force and energy for groups A and B
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
|
||||
{
|
||||
if (slabflag)
|
||||
error->all(FLERR,"Cannot (yet) use K-space slab "
|
||||
"correction with compute group/group");
|
||||
|
||||
int i,j;
|
||||
|
||||
if (!group_allocate_flag) {
|
||||
allocate_groups();
|
||||
group_allocate_flag = 1;
|
||||
}
|
||||
|
||||
e2group = 0; //energy
|
||||
f2group[0] = 0; //force in x-direction
|
||||
f2group[1] = 0; //force in y-direction
|
||||
f2group[2] = 0; //force in z-direction
|
||||
|
||||
double *q = atom->q;
|
||||
int nlocal = atom->nlocal;
|
||||
int *mask = atom->mask;
|
||||
|
||||
// convert atoms from box to lamda coords
|
||||
|
||||
if (triclinic == 0) boxlo = domain->boxlo;
|
||||
else {
|
||||
boxlo = domain->boxlo_lamda;
|
||||
domain->x2lamda(atom->nlocal);
|
||||
}
|
||||
|
||||
// map my particle charge onto my local 3d density grid
|
||||
|
||||
make_rho_groups(groupbit_A,groupbit_B,BA_flag);
|
||||
|
||||
// all procs communicate density values from their ghost cells
|
||||
// to fully sum contribution in their 3d bricks
|
||||
// remap from 3d decomposition to FFT decomposition
|
||||
|
||||
// temporarily store and switch pointers so we can
|
||||
// use brick2fft() for groups A and B (without
|
||||
// writing an additional function)
|
||||
|
||||
FFT_SCALAR ***density_brick_real = density_brick;
|
||||
FFT_SCALAR *density_fft_real = density_fft;
|
||||
|
||||
// group A
|
||||
|
||||
density_brick = density_A_brick;
|
||||
density_fft = density_A_fft;
|
||||
|
||||
brick2fft();
|
||||
|
||||
// group B
|
||||
|
||||
density_brick = density_B_brick;
|
||||
density_fft = density_B_fft;
|
||||
|
||||
brick2fft();
|
||||
|
||||
// switch back pointers
|
||||
|
||||
density_brick = density_brick_real;
|
||||
density_fft = density_fft_real;
|
||||
|
||||
// compute potential gradient on my FFT grid and
|
||||
// portion of group-group energy/force on this proc's FFT grid
|
||||
|
||||
poisson_groups(BA_flag);
|
||||
|
||||
const double qscale = force->qqrd2e * scale;
|
||||
|
||||
// total group A <--> group B energy
|
||||
|
||||
double e2group_all;
|
||||
MPI_Allreduce(&e2group,&e2group_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
e2group = e2group_all;
|
||||
|
||||
e2group *= 0.5*volume;
|
||||
|
||||
// total charge of groups A & B, needed for self-energy correction
|
||||
|
||||
double qsum_group = 0.0, qsqsum_group = 0.0;
|
||||
|
||||
for (int i = 0; i < atom->nlocal; i++) {
|
||||
if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B)) {
|
||||
if (BA_flag) continue;
|
||||
qsqsum_group += q[i]*q[i];
|
||||
}
|
||||
}
|
||||
|
||||
double tmp;
|
||||
MPI_Allreduce(&qsum_group,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
qsum_group = tmp;
|
||||
|
||||
// self-energy correction
|
||||
|
||||
e2group -= g_ewald*qsqsum_group/MY_PIS;
|
||||
|
||||
// should we include this boundary condition or not?
|
||||
// e2group += MY_PI2*qsum_group*qsum_group / (g_ewald*g_ewald*volume);
|
||||
|
||||
e2group *= qscale;
|
||||
|
||||
// total group A <--> group B force
|
||||
|
||||
double f2group_all[3];
|
||||
MPI_Allreduce(f2group,f2group_all,3,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
for (i = 0; i < 3; i++) f2group[i] = qscale*volume*f2group_all[i];
|
||||
|
||||
// convert atoms back from lamda to box coords
|
||||
|
||||
if (triclinic) domain->lamda2x(atom->nlocal);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate group-group memory that depends on # of K-vectors and order
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPM::allocate_groups()
|
||||
{
|
||||
memory->create3d_offset(density_A_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
|
||||
nxlo_out,nxhi_out,"pppm:density_A_brick");
|
||||
memory->create3d_offset(density_B_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out,
|
||||
nxlo_out,nxhi_out,"pppm:density_B_brick");
|
||||
memory->create(density_A_fft,nfft_both,"pppm:density_A_fft");
|
||||
memory->create(density_B_fft,nfft_both,"pppm:density_B_fft");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
deallocate group-group memory that depends on # of K-vectors and order
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPM::deallocate_groups()
|
||||
{
|
||||
memory->destroy3d_offset(density_A_brick,nzlo_out,nylo_out,nxlo_out);
|
||||
memory->destroy3d_offset(density_B_brick,nzlo_out,nylo_out,nxlo_out);
|
||||
memory->destroy(density_A_fft);
|
||||
memory->destroy(density_B_fft);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
create discretized "density" on section of global grid due to my particles
|
||||
density(x,y,z) = charge "density" at grid points of my 3d brick
|
||||
(nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
|
||||
in global grid for group-group interactions
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPM::make_rho_groups(int groupbit_A, int groupbit_B, int BA_flag)
|
||||
{
|
||||
int l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
||||
|
||||
// clear 3d density arrays
|
||||
|
||||
memset(&(density_A_brick[nzlo_out][nylo_out][nxlo_out]),0,
|
||||
ngrid*sizeof(FFT_SCALAR));
|
||||
|
||||
memset(&(density_B_brick[nzlo_out][nylo_out][nxlo_out]),0,
|
||||
ngrid*sizeof(FFT_SCALAR));
|
||||
|
||||
// loop over my charges, add their contribution to nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
|
||||
double *q = atom->q;
|
||||
double **x = atom->x;
|
||||
int nlocal = atom->nlocal;
|
||||
int *mask = atom->mask;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
|
||||
if ((mask[i] & groupbit_A) && (mask[i] & groupbit_B))
|
||||
if (BA_flag) continue;
|
||||
|
||||
if ((mask[i] & groupbit_A) || (mask[i] & groupbit_B)) {
|
||||
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz);
|
||||
|
||||
z0 = delvolinv * q[i];
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
y0 = z0*rho1d[2][n];
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
x0 = y0*rho1d[1][m];
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
|
||||
// group A
|
||||
|
||||
if (mask[i] & groupbit_A)
|
||||
density_A_brick[mz][my][mx] += x0*rho1d[0][l];
|
||||
|
||||
// group B
|
||||
|
||||
if (mask[i] & groupbit_B)
|
||||
density_B_brick[mz][my][mx] += x0*rho1d[0][l];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
FFT-based Poisson solver for group-group interactions
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPM::poisson_groups(int BA_flag)
|
||||
{
|
||||
int i,j,k,n;
|
||||
double eng;
|
||||
|
||||
// reuse memory (already declared)
|
||||
|
||||
FFT_SCALAR *work_A = work1;
|
||||
FFT_SCALAR *work_B = work2;
|
||||
|
||||
// transform charge density (r -> k)
|
||||
|
||||
// group A
|
||||
|
||||
n = 0;
|
||||
for (i = 0; i < nfft; i++) {
|
||||
work_A[n++] = density_A_fft[i];
|
||||
work_A[n++] = ZEROF;
|
||||
}
|
||||
|
||||
fft1->compute(work_A,work_A,1);
|
||||
|
||||
// group B
|
||||
|
||||
n = 0;
|
||||
for (i = 0; i < nfft; i++) {
|
||||
work_B[n++] = density_B_fft[i];
|
||||
work_B[n++] = ZEROF;
|
||||
}
|
||||
|
||||
fft1->compute(work_B,work_B,1);
|
||||
|
||||
// group-group energy and force contribution,
|
||||
// keep everything in reciprocal space so
|
||||
// no inverse FFTs needed
|
||||
|
||||
double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm);
|
||||
double s2 = scaleinv*scaleinv;
|
||||
|
||||
// energy
|
||||
|
||||
n = 0;
|
||||
for (i = 0; i < nfft; i++) {
|
||||
e2group += s2 * greensfn[i] *
|
||||
(work_A[n]*work_B[n] + work_A[n+1]*work_B[n+1]);
|
||||
n += 2;
|
||||
}
|
||||
|
||||
if (BA_flag) return;
|
||||
|
||||
|
||||
// multiply by Green's function and s2
|
||||
// (only for work_A so it is not squared below)
|
||||
|
||||
n = 0;
|
||||
for (i = 0; i < nfft; i++) {
|
||||
work_A[n++] *= s2 * greensfn[i];
|
||||
work_A[n++] *= s2 * greensfn[i];
|
||||
}
|
||||
|
||||
double partial_group;
|
||||
|
||||
// force, x direction
|
||||
|
||||
n = 0;
|
||||
for (k = nzlo_fft; k <= nzhi_fft; k++)
|
||||
for (j = nylo_fft; j <= nyhi_fft; j++)
|
||||
for (i = nxlo_fft; i <= nxhi_fft; i++) {
|
||||
partial_group = work_A[n+1]*work_B[n] - work_A[n]*work_B[n+1];
|
||||
f2group[0] += fkx[i] * partial_group;
|
||||
n += 2;
|
||||
}
|
||||
|
||||
// force, y direction
|
||||
|
||||
n = 0;
|
||||
for (k = nzlo_fft; k <= nzhi_fft; k++)
|
||||
for (j = nylo_fft; j <= nyhi_fft; j++)
|
||||
for (i = nxlo_fft; i <= nxhi_fft; i++) {
|
||||
partial_group = work_A[n+1]*work_B[n] - work_A[n]*work_B[n+1];
|
||||
f2group[1] += fky[j] * partial_group;
|
||||
n += 2;
|
||||
}
|
||||
|
||||
// force, z direction
|
||||
|
||||
n = 0;
|
||||
for (k = nzlo_fft; k <= nzhi_fft; k++)
|
||||
for (j = nylo_fft; j <= nyhi_fft; j++)
|
||||
for (i = nxlo_fft; i <= nxhi_fft; i++) {
|
||||
partial_group = work_A[n+1]*work_B[n] - work_A[n]*work_B[n+1];
|
||||
f2group[2] += fkz[k] * partial_group;
|
||||
n += 2;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -45,6 +45,8 @@ class PPPM : public KSpace {
|
|||
virtual void timing(int, double &, double &);
|
||||
virtual double memory_usage();
|
||||
|
||||
virtual void compute_group_group(int, int, int);
|
||||
|
||||
protected:
|
||||
int me,nprocs;
|
||||
int nfactors;
|
||||
|
|
@ -79,6 +81,13 @@ class PPPM : public KSpace {
|
|||
double *gf_b;
|
||||
FFT_SCALAR **rho1d,**rho_coeff;
|
||||
|
||||
// group-group interactions
|
||||
|
||||
int group_allocate_flag;
|
||||
FFT_SCALAR ***density_A_brick,***density_B_brick;
|
||||
FFT_SCALAR *density_A_fft,*density_B_fft;
|
||||
|
||||
|
||||
class FFT3d *fft1,*fft2;
|
||||
class Remap *remap;
|
||||
|
||||
|
|
@ -117,6 +126,13 @@ class PPPM : public KSpace {
|
|||
void compute_rho_coeff();
|
||||
void slabcorr();
|
||||
|
||||
// group-group interactions
|
||||
|
||||
virtual void allocate_groups();
|
||||
virtual void deallocate_groups();
|
||||
virtual void make_rho_groups(int, int, int);
|
||||
virtual void poisson_groups(int);
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
denominator for Hockney-Eastwood Green's function
|
||||
of x,y,z = sin(kx*deltax/2), etc
|
||||
|
|
|
|||
|
|
@ -13,6 +13,7 @@
|
|||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Naveen Michaud-Agrawal (Johns Hopkins U)
|
||||
K-space terms added by Stan Moore (BYU)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "mpi.h"
|
||||
|
|
@ -26,6 +27,7 @@
|
|||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
#include "group.h"
|
||||
#include "kspace.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
|
@ -35,7 +37,7 @@ using namespace LAMMPS_NS;
|
|||
ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg != 4) error->all(FLERR,"Illegal compute group/group command");
|
||||
if (narg < 4) error->all(FLERR,"Illegal compute group/group command");
|
||||
|
||||
scalar_flag = vector_flag = 1;
|
||||
size_vector = 3;
|
||||
|
|
@ -51,6 +53,28 @@ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) :
|
|||
error->all(FLERR,"Compute group/group group ID does not exist");
|
||||
jgroupbit = group->bitmask[jgroup];
|
||||
|
||||
pairflag = 1;
|
||||
kspaceflag = 0;
|
||||
|
||||
int iarg = 4;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"pair") == 0) {
|
||||
if (iarg+2 > narg)
|
||||
error->all(FLERR,"Illegal compute group/group command");
|
||||
if (strcmp(arg[iarg+1],"yes") == 0) pairflag = 1;
|
||||
else if (strcmp(arg[iarg+1],"no") == 0) pairflag = 0;
|
||||
else error->all(FLERR,"Illegal compute group/group command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"kspace") == 0) {
|
||||
if (iarg+2 > narg)
|
||||
error->all(FLERR,"Illegal compute group/group command");
|
||||
if (strcmp(arg[iarg+1],"yes") == 0) pairflag = 1;
|
||||
else if (strcmp(arg[iarg+1],"no") == 0) pairflag = 0;
|
||||
else error->all(FLERR,"Illegal compute group/group command");
|
||||
iarg += 2;
|
||||
} else error->all(FLERR,"Illegal compute group/group command");
|
||||
}
|
||||
|
||||
vector = new double[3];
|
||||
}
|
||||
|
||||
|
|
@ -66,17 +90,28 @@ ComputeGroupGroup::~ComputeGroupGroup()
|
|||
|
||||
void ComputeGroupGroup::init()
|
||||
{
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"No pair style defined for compute group/group");
|
||||
|
||||
// if non-hybrid, then error if single_enable = 0
|
||||
// if hybrid, let hybrid determine if sub-style sets single_enable = 0
|
||||
|
||||
if (pairflag && force->pair == NULL)
|
||||
error->all(FLERR,"No pair style defined for compute group/group");
|
||||
if (force->pair_match("hybrid",0) == NULL && force->pair->single_enable == 0)
|
||||
error->all(FLERR,"Pair style does not support compute group/group");
|
||||
|
||||
// error if Kspace style does not compute group/group interactions
|
||||
|
||||
if (kspaceflag && force->kspace == NULL)
|
||||
error->all(FLERR,"No Kspace style defined for compute group/group");
|
||||
if (kspaceflag && force->kspace->group_group_enable == 0)
|
||||
error->all(FLERR,"Kspace style does not support compute group/group");
|
||||
|
||||
if (pairflag) {
|
||||
pair = force->pair;
|
||||
cutsq = force->pair->cutsq;
|
||||
} else pair = NULL;
|
||||
|
||||
if (kspaceflag) kspace = force->kspace;
|
||||
else kspace = NULL;
|
||||
|
||||
// recheck that group 2 has not been deleted
|
||||
|
||||
|
|
@ -87,10 +122,12 @@ void ComputeGroupGroup::init()
|
|||
|
||||
// need an occasional half neighbor list
|
||||
|
||||
if (pairflag) {
|
||||
int irequest = neighbor->request((void *) this);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -106,7 +143,12 @@ double ComputeGroupGroup::compute_scalar()
|
|||
{
|
||||
invoked_scalar = invoked_vector = update->ntimestep;
|
||||
|
||||
interact();
|
||||
scalar = 0.0;
|
||||
vector[0] = vector[1] = vector[2] = 0.0;
|
||||
|
||||
if (pairflag) pair_contribution();
|
||||
if (kspaceflag) kspace_contribution();
|
||||
|
||||
return scalar;
|
||||
}
|
||||
|
||||
|
|
@ -116,12 +158,16 @@ void ComputeGroupGroup::compute_vector()
|
|||
{
|
||||
invoked_scalar = invoked_vector = update->ntimestep;
|
||||
|
||||
interact();
|
||||
scalar = 0.0;
|
||||
vector[0] = vector[1] = vector[2] = 0.0;
|
||||
|
||||
if (pairflag) pair_contribution();
|
||||
if (kspaceflag) kspace_contribution();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeGroupGroup::interact()
|
||||
void ComputeGroupGroup::pair_contribution()
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz;
|
||||
|
|
@ -148,7 +194,7 @@ void ComputeGroupGroup::interact()
|
|||
// loop over neighbors of my atoms
|
||||
// skip if I,J are not in 2 groups
|
||||
|
||||
double one[4],all[4];
|
||||
double one[4];
|
||||
one[0] = one[1] = one[2] = one[3] = 0.0;
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
|
|
@ -212,7 +258,28 @@ void ComputeGroupGroup::interact()
|
|||
}
|
||||
}
|
||||
|
||||
double all[4];
|
||||
MPI_Allreduce(one,all,4,MPI_DOUBLE,MPI_SUM,world);
|
||||
scalar = all[0];
|
||||
vector[0] = all[1]; vector[1] = all[2]; vector[2] = all[3];
|
||||
scalar += all[0];
|
||||
vector[0] += all[1]; vector[1] += all[2]; vector[2] += all[3];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeGroupGroup::kspace_contribution()
|
||||
{
|
||||
double *vector_kspace = force->kspace->f2group;
|
||||
|
||||
force->kspace->compute_group_group(groupbit,jgroupbit,0);
|
||||
scalar += force->kspace->e2group;
|
||||
vector[0] += vector_kspace[0];
|
||||
vector[1] += vector_kspace[1];
|
||||
vector[2] += vector_kspace[2];
|
||||
|
||||
// compute extra B <--> A Kspace interaction so energy matches
|
||||
// real-space style of compute group-group
|
||||
// add extra Kspace term to energy
|
||||
|
||||
force->kspace->compute_group_group(groupbit,jgroupbit,1);
|
||||
scalar += force->kspace->e2group;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -37,10 +37,14 @@ class ComputeGroupGroup : public Compute {
|
|||
char *group2;
|
||||
int jgroup,jgroupbit,othergroupbit;
|
||||
double **cutsq;
|
||||
int pairflag,kspaceflag;
|
||||
class Pair *pair;
|
||||
class NeighList *list;
|
||||
class KSpace *kspace;
|
||||
|
||||
void interact();
|
||||
void pair_contribution();
|
||||
void kspace_contribution();
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -31,6 +31,8 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
|
|||
virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0;
|
||||
|
||||
compute_flag = 1;
|
||||
group_group_enable = 0;
|
||||
|
||||
order = 5;
|
||||
gridflag = 0;
|
||||
gewaldflag = 0;
|
||||
|
|
|
|||
|
|
@ -25,9 +25,12 @@ class KSpace : protected Pointers {
|
|||
double energy; // accumulated energy
|
||||
double virial[6]; // accumlated virial
|
||||
double *eatom,**vatom; // accumulated per-atom energy/virial
|
||||
double e2group; // accumulated group-group energy
|
||||
double f2group[3]; // accumulated group-group force
|
||||
|
||||
double g_ewald;
|
||||
int nx_pppm,ny_pppm,nz_pppm;
|
||||
int group_group_enable; // 1 if style supports group/group calculation
|
||||
|
||||
int compute_flag; // 0 if skip compute()
|
||||
|
||||
|
|
@ -42,6 +45,7 @@ class KSpace : protected Pointers {
|
|||
virtual void init() = 0;
|
||||
virtual void setup() = 0;
|
||||
virtual void compute(int, int) = 0;
|
||||
virtual void compute_group_group(int, int, int) {};
|
||||
virtual void timing(int, double &, double &) {}
|
||||
virtual double memory_usage() {return 0.0;}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue