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

Updated pppm/disp/dielectric for long-range energy calculation (eflag_global is true)

Browse Source
This commit is contained in:
Trung Nguyen
2022-12-30 16:02:26 -06:00
parent 652c237b5e
commit d5e897ccbb
4 changed files with 80 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files

136
src/DIELECTRIC/pppm_dielectric.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/DIELECTRIC/pppm_dielectric.h
View File

@ -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;

72
src/DIELECTRIC/pppm_disp_dielectric.cpp
View File

@ -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;

4
src/DIELECTRIC/pppm_disp_dielectric.h
View File

@ -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;
};