Skip to content

1. refactor H_Hartree_pw::v_hartree() #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mohanchen merged 3 commits into from
May 6, 2021
Merged

1. refactor H_Hartree_pw::v_hartree() #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f53310c
1. refactor H_Hartree_pw::v_hartree()
PeizeLin May 6, 2021
fa068ed
Merge branch 'potential'
PeizeLin May 6, 2021
be8772f
1. refactor Potential_Libxc
PeizeLin May 6, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 23 additions & 55 deletions ABACUS.develop/source/src_pw/H_Hartree_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,41 +5,29 @@ double H_Hartree_pw::hartree_energy=0.0;
//--------------------------------------------------------------------
// Transform charge density to hartree potential.
//--------------------------------------------------------------------
void H_Hartree_pw::v_hartree(
matrix H_Hartree_pw::v_hartree(
const UnitCell &cell,
PW_Basis &pwb,
const Use_FFT &ufft,
const int &nspin,
matrix &v,
double** rho)
const int &nspin,
const double*const*const rho)
{
TITLE("H_Hartree_pw","v_hartree");
timer::tick("H_Hartree_pw","v_hartree");

complex<double> *Porter = ufft.porter;

// Hartree potential VH(r) from n(r)
ZEROS( Porter, pwb.nrxx );
int nspin0 = 1;
if(nspin==2) nspin0 = nspin;
vector<complex<double>> Porter(pwb.nrxx);
const int nspin0 = (nspin==2) ? 2 : 1;
for(int is=0; is<nspin0; is++)
{
for (int ir=0; ir<pwb.nrxx; ir++)
{
Porter[ir] += complex<double>( rho[is][ir], 0.0 );
}
}

//=============================
// bring rho (aux) to G space
//=============================
pwb.FFT_chg.FFT3D(Porter, -1);
pwb.FFT_chg.FFT3D(Porter.data(), -1);

//double charge;
//if (pwb.gstart == 1)
//{
// charge = cell.omega * Porter[pwb.ig2fftc[0]].real();
//}
//OUT(ofs_running, "v_h charge", charge);

//=======================================================
Expand All @@ -48,9 +36,7 @@ void H_Hartree_pw::v_hartree(

double ehart = 0.0;

complex<double> *vh_g = new complex<double>[pwb.ngmc];
ZEROS(vh_g, pwb.ngmc);

vector<complex<double>> vh_g(pwb.ngmc);
for (int ig = pwb.gstart; ig<pwb.ngmc; ig++)
{
const int j = pwb.ig2fftc[ig];
Expand All @@ -65,52 +51,37 @@ void H_Hartree_pw::v_hartree(

Parallel_Reduce::reduce_double_pool( ehart );
ehart *= 0.5 * cell.omega;

//cout << " ehart=" << ehart << endl;

H_Hartree_pw::hartree_energy = ehart;



ZEROS(Porter, pwb.nrxx);

std::fill( Porter.begin(), Porter.end(), complex<double>(0.0,0.0) );
for (int ig = 0;ig < pwb.ngmc;ig++)
{
Porter[pwb.ig2fftc[ig]] = vh_g[ig];
}

//==========================================
//transform hartree potential to real space
//==========================================
pwb.FFT_chg.FFT3D(Porter, 1);
pwb.FFT_chg.FFT3D(Porter.data(), 1);

//==========================================
//Add hartree potential to the xc potential
//==========================================

if(nspin==4)
{
for (int ir = 0;ir < pwb.nrxx;ir++)
{
v(0, ir) += Porter[ir].real();
}
}
else
{
for (int is = 0;is < nspin;is++)
{
for (int ir = 0;ir < pwb.nrxx;ir++)
{
v(is, ir) += Porter[ir].real();
}
}
}

matrix v(nspin, pwb.nrxx);
if(nspin==4)
{
for (int ir = 0;ir < pwb.nrxx;ir++)
v(0, ir) = Porter[ir].real();
}
else
{
for (int is = 0;is < nspin;is++)
for (int ir = 0;ir < pwb.nrxx;ir++)
v(is, ir) = Porter[ir].real();
}

//-----------------------------------------------------------
// we need to add this out_potential funciton back
// in near future, 2021-02-25
//-----------------------------------------------------------

//-------------------------------------------
// output the Hartree potential into a file.
//-------------------------------------------
Expand All @@ -135,8 +106,5 @@ void H_Hartree_pw::v_hartree(
*/

timer::tick("H_Hartree_pw","v_hartree");

delete[] vh_g;

return;
return v;
} // end subroutine v_h
8 changes: 3 additions & 5 deletions ABACUS.develop/source/src_pw/H_Hartree_pw.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,11 @@ class H_Hartree_pw
static double hartree_energy;

// compute the Hartree energy
static void v_hartree(
static matrix v_hartree(
const UnitCell &cell,
PW_Basis &pwb,
const Use_FFT &ufft,
const int &nspin,
matrix &v,
double** rho);
const int &nspin,
const double*const*const rho);

private:

Expand Down
9 changes: 7 additions & 2 deletions ABACUS.develop/source/src_pw/potential.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -292,15 +292,20 @@ void Potential::v_of_rho
//----------------------------------------------------------

#ifdef USE_LIBXC
Potential_Libxc::v_xc(rho_in, H_XC_pw::etxc, H_XC_pw::vtxc, v_in);
H_XC_pw::etxc = 0.0;
H_XC_pw::vtxc = 0.0;
const auto etxc_vtxc_v = Potential_Libxc::v_xc(rho_in, CHR.rho_core);
H_XC_pw::etxc += std::get<0>(etxc_vtxc_v);
H_XC_pw::vtxc += std::get<1>(etxc_vtxc_v);
v_in += std::get<2>(etxc_vtxc_v);
#else
H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core, v_in);
#endif

//----------------------------------------------------------
// calculate the Hartree potential
//----------------------------------------------------------
H_Hartree_pw::v_hartree(ucell, pw, UFFT, NSPIN, v_in, rho_in);
v_in += H_Hartree_pw::v_hartree(ucell, pw, NSPIN, rho_in);

// mohan add 2011-06-20
if(EFIELD && DIPOLE)
Expand Down
Loading