diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 806e3a549a..eb6924be90 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -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; +} diff --git a/src/KSPACE/ewald.h b/src/KSPACE/ewald.h index 360c0385bc..e8cbdecf15 100644 --- a/src/KSPACE/ewald.h +++ b/src/KSPACE/ewald.h @@ -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(); }; } diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 94e71f29cb..4aced20578 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -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; + } +} diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index 852f303d5d..356cda7735 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -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 diff --git a/src/compute_group_group.cpp b/src/compute_group_group.cpp index 9ee8075fad..4cbdeb6997 100644 --- a/src/compute_group_group.cpp +++ b/src/compute_group_group.cpp @@ -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"); - pair = force->pair; - cutsq = force->pair->cutsq; + // 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 - int irequest = neighbor->request((void *) this); - neighbor->requests[irequest]->pair = 0; - neighbor->requests[irequest]->compute = 1; - neighbor->requests[irequest]->occasional = 1; + 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; } diff --git a/src/compute_group_group.h b/src/compute_group_group.h index cddcdda298..1611dc82b7 100644 --- a/src/compute_group_group.h +++ b/src/compute_group_group.h @@ -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(); }; } diff --git a/src/kspace.cpp b/src/kspace.cpp index a6f02b7ced..f6ec055520 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -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; diff --git a/src/kspace.h b/src/kspace.h index b9ffb5a5cf..851cad7dd3 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -25,10 +25,13 @@ 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() KSpace(class LAMMPS *, int, char **); @@ -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;}