From d5e897ccbbc629bbcf9cd99afd2e8b5710d1a5a1 Mon Sep 17 00:00:00 2001 From: Trung Nguyen Date: Fri, 30 Dec 2022 16:02:26 -0600 Subject: [PATCH] Updated pppm/disp/dielectric for long-range energy calculation (eflag_global is true) --- src/DIELECTRIC/pppm_dielectric.cpp | 136 +----------------------- src/DIELECTRIC/pppm_dielectric.h | 1 - src/DIELECTRIC/pppm_disp_dielectric.cpp | 72 +++++++++++++ src/DIELECTRIC/pppm_disp_dielectric.h | 4 + 4 files changed, 80 insertions(+), 133 deletions(-) diff --git a/src/DIELECTRIC/pppm_dielectric.cpp b/src/DIELECTRIC/pppm_dielectric.cpp index 6b274e773d..1c1ad97815 100644 --- a/src/DIELECTRIC/pppm_dielectric.cpp +++ b/src/DIELECTRIC/pppm_dielectric.cpp @@ -188,6 +188,7 @@ void PPPMDielectric::compute(int eflag, int vflag) energy = energy_before_poisson; // switch to unscaled charges to find charge density + use_qscaled = false; // redo the charge density @@ -201,7 +202,9 @@ void PPPMDielectric::compute(int eflag, int vflag) brick2fft(); - poisson(); // poisson computes elong with unscaled charges + // compute electrostatic energy with the unscaled charges and average epsilon + + poisson(); double energy_all; MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); @@ -383,137 +386,6 @@ void PPPMDielectric::make_rho() } } -/* ---------------------------------------------------------------------- - FFT-based Poisson solver for ik -------------------------------------------------------------------------- */ - -void PPPMDielectric::poisson_ik() -{ - int i,j,k,n; - double eng; - - // transform charge density (r -> k) - - n = 0; - for (i = 0; i < nfft; i++) { - work1[n++] = density_fft[i]; - work1[n++] = ZEROF; - } - - fft1->compute(work1,work1,FFT3d::FORWARD); - - // global energy and virial contribution - - double scaleinv = 1.0/(nx_pppm*ny_pppm*nz_pppm); - double s2 = scaleinv*scaleinv; - - if (eflag_global || vflag_global) { - if (vflag_global) { - n = 0; - for (i = 0; i < nfft; i++) { - eng = s2 * greensfn[i] * (work1[n]*work1[n] + work1[n+1]*work1[n+1]); - for (j = 0; j < 6; j++) virial[j] += eng*vg[i][j]; - if (eflag_global) energy += eng; - n += 2; - } - } else { - n = 0; - for (i = 0; i < nfft; i++) { - energy += - s2 * greensfn[i] * (work1[n]*work1[n] + work1[n+1]*work1[n+1]); - n += 2; - } - } - } - - // scale by 1/total-grid-pts to get rho(k) - // multiply by Green's function to get V(k) - - n = 0; - for (i = 0; i < nfft; i++) { - work1[n++] *= scaleinv * greensfn[i]; - work1[n++] *= scaleinv * greensfn[i]; - } - - // extra FFTs for per-atom energy/virial - - if (evflag_atom) poisson_peratom(); - - // triclinic system - - if (triclinic) { - poisson_ik_triclinic(); - return; - } - - // compute gradients of V(r) in each of 3 dims by transforming ik*V(k) - // FFT leaves data in 3d brick decomposition - // copy it into inner portion of vdx,vdy,vdz arrays - - // x direction gradient - - 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++) { - work2[n] = -fkx[i]*work1[n+1]; - work2[n+1] = fkx[i]*work1[n]; - n += 2; - } - - fft2->compute(work2,work2,FFT3d::BACKWARD); - - n = 0; - for (k = nzlo_in; k <= nzhi_in; k++) - for (j = nylo_in; j <= nyhi_in; j++) - for (i = nxlo_in; i <= nxhi_in; i++) { - vdx_brick[k][j][i] = work2[n]; - n += 2; - } - - // y direction gradient - - 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++) { - work2[n] = -fky[j]*work1[n+1]; - work2[n+1] = fky[j]*work1[n]; - n += 2; - } - - fft2->compute(work2,work2,FFT3d::BACKWARD); - - n = 0; - for (k = nzlo_in; k <= nzhi_in; k++) - for (j = nylo_in; j <= nyhi_in; j++) - for (i = nxlo_in; i <= nxhi_in; i++) { - vdy_brick[k][j][i] = work2[n]; - n += 2; - } - - // z direction gradient - - 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++) { - work2[n] = -fkz[k]*work1[n+1]; - work2[n+1] = fkz[k]*work1[n]; - n += 2; - } - - fft2->compute(work2,work2,FFT3d::BACKWARD); - - n = 0; - for (k = nzlo_in; k <= nzhi_in; k++) - for (j = nylo_in; j <= nyhi_in; j++) - for (i = nxlo_in; i <= nxhi_in; i++) { - vdz_brick[k][j][i] = work2[n]; - n += 2; - } -} - /* ---------------------------------------------------------------------- interpolate from grid to get electric field & force on my particles for ik ------------------------------------------------------------------------- */ diff --git a/src/DIELECTRIC/pppm_dielectric.h b/src/DIELECTRIC/pppm_dielectric.h index 4221723b1c..c266488561 100644 --- a/src/DIELECTRIC/pppm_dielectric.h +++ b/src/DIELECTRIC/pppm_dielectric.h @@ -39,7 +39,6 @@ class PPPMDielectric : public PPPM { void make_rho() override; void fieldforce_ik() override; void fieldforce_ad() override; - void poisson_ik() override; void qsum_qsq(int warning_flag = 1) override; class AtomVecDielectric *avec; diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp index e7cf2f410e..14faee8a75 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.cpp +++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp @@ -416,11 +416,56 @@ void PPPMDispDielectric::compute(int eflag, int vflag) natoms_original = atom->natoms; } + // recompute the average epsilon of all the atoms + + compute_ave_epsilon(); + // sum energy across procs and add in volume-dependent term const double qscale = force->qqrd2e * scale; if (eflag_global) { + + if (function[0]) { + + // switch to unscaled charges to find charge density + + use_qscaled = false; + + make_rho_c(); + + gc->reverse_comm(Grid3d::KSPACE,this,REVERSE_RHO,1,sizeof(FFT_SCALAR), + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + + brick2fft(nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in, + density_brick,density_fft,work1,remap); + + // compute electrostatic energy with the unscaled charges and average epsilon + + if (differentiation_flag == 1) { + poisson_ad(work1,work2,density_fft,fft1,fft2, + nx_pppm,ny_pppm,nz_pppm,nfft, + nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft, + nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in, + energy_1,greensfn, + virial_1,vg,vg2, + u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick); + + } else { + poisson_ik(work1,work2,density_fft,fft1,fft2, + nx_pppm,ny_pppm,nz_pppm,nfft, + nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft, + nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in, + energy_1,greensfn, + fkx,fky,fkz,fkx2,fky2,fkz2, + vdx_brick,vdy_brick,vdz_brick,virial_1,vg,vg2, + u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick); + + gc->forward_comm(Grid3d::KSPACE,this,FORWARD_IK,3,sizeof(FFT_SCALAR), + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + } + } + double energy_all; MPI_Allreduce(&energy_1,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy_1 = energy_all; @@ -435,6 +480,10 @@ void PPPMDispDielectric::compute(int eflag, int vflag) energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij + 1.0/12.0*pow(g_ewald_6,6)*csum; energy_1 *= qscale; + + // revert to qscaled charges (for force in the next time step) + + use_qscaled = true; } // sum virial across procs @@ -495,6 +544,28 @@ void PPPMDispDielectric::compute(int eflag, int vflag) if (triclinic) domain->lamda2x(atom->nlocal); } +/* ---------------------------------------------------------------------- + compute the average dielectric constant of all the atoms + NOTE: for dielectric use cases +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::compute_ave_epsilon() +{ + const double * const epsilon = atom->epsilon; + const int nlocal = atom->nlocal; + double epsilon_local(0.0); + +#if defined(_OPENMP) +#pragma omp parallel for default(shared) reduction(+:epsilon_local) +#endif + for (int i = 0; i < nlocal; i++) { + epsilon_local += epsilon[i]; + } + + MPI_Allreduce(&epsilon_local,&epsilon_ave,1,MPI_DOUBLE,MPI_SUM,world); + epsilon_ave /= (double)atom->natoms; +} + /* ---------------------------------------------------------------------- compute qsum,qsqsum,q2 and give error/warning if not charge neutral called initially, when particle count changes, when charges are changed @@ -564,6 +635,7 @@ void PPPMDispDielectric::make_rho_c() // (mx,my,mz) = global coords of moving stencil pt double *q = atom->q_scaled; + if (use_qscaled == false) q = atom->q; double **x = atom->x; int nlocal = atom->nlocal; diff --git a/src/DIELECTRIC/pppm_disp_dielectric.h b/src/DIELECTRIC/pppm_disp_dielectric.h index c93fb6a65a..7412762818 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.h +++ b/src/DIELECTRIC/pppm_disp_dielectric.h @@ -44,6 +44,10 @@ class PPPMDispDielectric : public PPPMDisp { void qsum_qsq(int warning_flag = 1) override; class AtomVecDielectric *avec; + bool use_qscaled; + + void compute_ave_epsilon(); + double epsilon_ave; int mu_flag; };