From b8f17c77b6e7cf79caf534b603d3503c7305bf20 Mon Sep 17 00:00:00 2001 From: linpz Date: Thu, 6 May 2021 21:25:42 +0800 Subject: [PATCH 1/3] 1. refactor class H_XC_pw --- ABACUS.develop/source/src_pw/H_XC_pw.cpp | 40 +++++++++---------- ABACUS.develop/source/src_pw/H_XC_pw.h | 7 ++-- ABACUS.develop/source/src_pw/forces.cpp | 6 ++- ABACUS.develop/source/src_pw/potential.cpp | 10 ++--- .../source/src_pw/stress_func_cc.cpp | 6 ++- 5 files changed, 35 insertions(+), 34 deletions(-) diff --git a/ABACUS.develop/source/src_pw/H_XC_pw.cpp b/ABACUS.develop/source/src_pw/H_XC_pw.cpp index 04a06ce90f..2475ae17b0 100644 --- a/ABACUS.develop/source/src_pw/H_XC_pw.cpp +++ b/ABACUS.develop/source/src_pw/H_XC_pw.cpp @@ -5,21 +5,21 @@ double H_XC_pw::etxc; double H_XC_pw::vtxc; -void H_XC_pw::v_xc +std::tuple H_XC_pw::v_xc ( const int &nrxx, // number of real-space grid const int &ncxyz, // total number of charge grid const double &omega, // volume of cell - double **rho_in, - double *rho_core, // core charge density - matrix &v) + const double*const*const rho_in, + const double*const rho_core) // core charge density { TITLE("H_XC_pw","v_xc"); timer::tick("H_XC_pw","v_xc"); //Exchange-Correlation potential Vxc(r) from n(r) - etxc = 0.0; - vtxc = 0.0; + double et_xc = 0.0; + double vt_xc = 0.0; + matrix v(NSPIN, nrxx); // the square of the e charge // in Rydeberg unit, so * 2.0. @@ -51,9 +51,9 @@ void H_XC_pw::v_xc XC_Functional::xc(arhox, ex, ec, vx[0], vc[0]); v(0,ir) = e2 * (vx[0] + vc[0]); // consider the total charge density - etxc += e2 * (ex + ec) * rhox; + et_xc += e2 * (ex + ec) * rhox; // only consider rho_in - vtxc += v(0, ir) * rho_in[0][ir]; + vt_xc += v(0, ir) * rho_in[0][ir]; } // endif } //enddo } @@ -97,9 +97,9 @@ void H_XC_pw::v_xc v(is, ir) = e2 * (vx[is] + vc[is]); } - etxc += e2 * (ex + ec) * rhox; + et_xc += e2 * (ex + ec) * rhox; - vtxc += v(0, ir) * rho_in[0][ir] + v(1, ir) * rho_in[1][ir]; + vt_xc += v(0, ir) * rho_in[0][ir] + v(1, ir) * rho_in[1][ir]; } } @@ -129,10 +129,10 @@ void H_XC_pw::v_xc XC_Functional::xc_spin( arhox, zeta, ex, ec, vx[0], vx[1], vc[0], vc[1] ); - etxc += e2 * ( ex + ec ) * rhox; + et_xc += e2 * ( ex + ec ) * rhox; v(0, ir) = e2*( 0.5 * ( vx[0] + vc[0] + vx[1] + vc[1] ) ); - vtxc += v(0,ir) * rho_in[0][ir]; + vt_xc += v(0,ir) * rho_in[0][ir]; double vs = 0.5 * ( vx[0] + vc[0] - vx[1] - vc[1] ); if ( amag > vanishing_charge ) @@ -140,7 +140,7 @@ void H_XC_pw::v_xc for(int ipol = 1;ipol< 4;ipol++) { v(ipol, ir) = e2 * vs * rho_in[ipol][ir] / amag; - vtxc += v(ipol,ir) * rho_in[ipol][ir]; + vt_xc += v(ipol,ir) * rho_in[ipol][ir]; }//end do }//end if }//end if @@ -150,16 +150,16 @@ void H_XC_pw::v_xc // add gradient corrections (if any) // mohan modify 2009-12-15 - GGA_PW::gradcorr(etxc, vtxc, v); + GGA_PW::gradcorr(et_xc, vt_xc, v); - // parallel code : collect vtxc,etxc + // parallel code : collect vt_xc,et_xc // mohan add 2008-06-01 - Parallel_Reduce::reduce_double_pool( etxc ); - Parallel_Reduce::reduce_double_pool( vtxc ); + Parallel_Reduce::reduce_double_pool( et_xc ); + Parallel_Reduce::reduce_double_pool( vt_xc ); - etxc *= omega / ncxyz; - vtxc *= omega / ncxyz; + et_xc *= omega / ncxyz; + vt_xc *= omega / ncxyz; timer::tick("H_XC_pw","v_xc"); - return; + return std::make_tuple(et_xc,vt_xc,v); } diff --git a/ABACUS.develop/source/src_pw/H_XC_pw.h b/ABACUS.develop/source/src_pw/H_XC_pw.h index a8f6dc0d00..bef5b032d7 100644 --- a/ABACUS.develop/source/src_pw/H_XC_pw.h +++ b/ABACUS.develop/source/src_pw/H_XC_pw.h @@ -25,13 +25,12 @@ class H_XC_pw static double vtxc; // compute the exchange-correlation energy - static void v_xc( + static std::tuple v_xc( const int &nrxx, // number of real-space grid const int &ncxyz, // total number of charge grid const double &omega, // volume of cell - double** rho_in, - double *rho_core, // core charge density - matrix &v); + const double*const*const rho_in, + const double*const rho_core); // core charge density }; diff --git a/ABACUS.develop/source/src_pw/forces.cpp b/ABACUS.develop/source/src_pw/forces.cpp index 2b8d9e46e3..d4c47c09d5 100644 --- a/ABACUS.develop/source/src_pw/forces.cpp +++ b/ABACUS.develop/source/src_pw/forces.cpp @@ -551,8 +551,10 @@ void Forces::cal_force_ew(matrix& forceion) void Forces::cal_force_cc(matrix& forcecc) { // recalculate the exchange-correlation potential. - matrix vxc(NSPIN, pw.nrxx); - H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core, vxc); + const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core); + H_XC_pw::etxc = std::get<0>(etxc_vtxc_v); // may delete? + H_XC_pw::vtxc = std::get<1>(etxc_vtxc_v); // may delete? + const matrix vxc = std::get<2>(etxc_vtxc_v); complex * psiv = new complex [pw.nrxx]; ZEROS(psiv, pw.nrxx); diff --git a/ABACUS.develop/source/src_pw/potential.cpp b/ABACUS.develop/source/src_pw/potential.cpp index 770d110756..a5635ee115 100644 --- a/ABACUS.develop/source/src_pw/potential.cpp +++ b/ABACUS.develop/source/src_pw/potential.cpp @@ -292,15 +292,13 @@ void Potential::v_of_rho //---------------------------------------------------------- #ifdef USE_LIBXC - 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); + const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core); #endif + 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); //---------------------------------------------------------- // calculate the Hartree potential diff --git a/ABACUS.develop/source/src_pw/stress_func_cc.cpp b/ABACUS.develop/source/src_pw/stress_func_cc.cpp index 0a09c7f13d..b8e8c359fb 100644 --- a/ABACUS.develop/source/src_pw/stress_func_cc.cpp +++ b/ABACUS.develop/source/src_pw/stress_func_cc.cpp @@ -34,8 +34,10 @@ void Stress_Func::stress_cc(matrix& sigma, const bool is_pw) } //recalculate the exchange-correlation potential - matrix vxc(NSPIN, pw.nrxx); - H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core, vxc); + const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core); + H_XC_pw::etxc = std::get<0>(etxc_vtxc_v); // may delete? + H_XC_pw::vtxc = std::get<1>(etxc_vtxc_v); // may delete? + const matrix vxc = std::get<2>(etxc_vtxc_v); complex * psic = new complex [pw.nrxx]; From 7593056846d136ae6d7175640efeca9e2aad31ab Mon Sep 17 00:00:00 2001 From: linpz Date: Thu, 6 May 2021 21:25:42 +0800 Subject: [PATCH 2/3] 1. refactor class H_XC_pw --- ABACUS.develop/source/src_pw/H_XC_pw.cpp | 40 +++++++++---------- ABACUS.develop/source/src_pw/H_XC_pw.h | 7 ++-- ABACUS.develop/source/src_pw/forces.cpp | 6 ++- ABACUS.develop/source/src_pw/potential.cpp | 10 ++--- .../source/src_pw/stress_func_cc.cpp | 6 ++- 5 files changed, 35 insertions(+), 34 deletions(-) diff --git a/ABACUS.develop/source/src_pw/H_XC_pw.cpp b/ABACUS.develop/source/src_pw/H_XC_pw.cpp index 04a06ce90f..2475ae17b0 100644 --- a/ABACUS.develop/source/src_pw/H_XC_pw.cpp +++ b/ABACUS.develop/source/src_pw/H_XC_pw.cpp @@ -5,21 +5,21 @@ double H_XC_pw::etxc; double H_XC_pw::vtxc; -void H_XC_pw::v_xc +std::tuple H_XC_pw::v_xc ( const int &nrxx, // number of real-space grid const int &ncxyz, // total number of charge grid const double &omega, // volume of cell - double **rho_in, - double *rho_core, // core charge density - matrix &v) + const double*const*const rho_in, + const double*const rho_core) // core charge density { TITLE("H_XC_pw","v_xc"); timer::tick("H_XC_pw","v_xc"); //Exchange-Correlation potential Vxc(r) from n(r) - etxc = 0.0; - vtxc = 0.0; + double et_xc = 0.0; + double vt_xc = 0.0; + matrix v(NSPIN, nrxx); // the square of the e charge // in Rydeberg unit, so * 2.0. @@ -51,9 +51,9 @@ void H_XC_pw::v_xc XC_Functional::xc(arhox, ex, ec, vx[0], vc[0]); v(0,ir) = e2 * (vx[0] + vc[0]); // consider the total charge density - etxc += e2 * (ex + ec) * rhox; + et_xc += e2 * (ex + ec) * rhox; // only consider rho_in - vtxc += v(0, ir) * rho_in[0][ir]; + vt_xc += v(0, ir) * rho_in[0][ir]; } // endif } //enddo } @@ -97,9 +97,9 @@ void H_XC_pw::v_xc v(is, ir) = e2 * (vx[is] + vc[is]); } - etxc += e2 * (ex + ec) * rhox; + et_xc += e2 * (ex + ec) * rhox; - vtxc += v(0, ir) * rho_in[0][ir] + v(1, ir) * rho_in[1][ir]; + vt_xc += v(0, ir) * rho_in[0][ir] + v(1, ir) * rho_in[1][ir]; } } @@ -129,10 +129,10 @@ void H_XC_pw::v_xc XC_Functional::xc_spin( arhox, zeta, ex, ec, vx[0], vx[1], vc[0], vc[1] ); - etxc += e2 * ( ex + ec ) * rhox; + et_xc += e2 * ( ex + ec ) * rhox; v(0, ir) = e2*( 0.5 * ( vx[0] + vc[0] + vx[1] + vc[1] ) ); - vtxc += v(0,ir) * rho_in[0][ir]; + vt_xc += v(0,ir) * rho_in[0][ir]; double vs = 0.5 * ( vx[0] + vc[0] - vx[1] - vc[1] ); if ( amag > vanishing_charge ) @@ -140,7 +140,7 @@ void H_XC_pw::v_xc for(int ipol = 1;ipol< 4;ipol++) { v(ipol, ir) = e2 * vs * rho_in[ipol][ir] / amag; - vtxc += v(ipol,ir) * rho_in[ipol][ir]; + vt_xc += v(ipol,ir) * rho_in[ipol][ir]; }//end do }//end if }//end if @@ -150,16 +150,16 @@ void H_XC_pw::v_xc // add gradient corrections (if any) // mohan modify 2009-12-15 - GGA_PW::gradcorr(etxc, vtxc, v); + GGA_PW::gradcorr(et_xc, vt_xc, v); - // parallel code : collect vtxc,etxc + // parallel code : collect vt_xc,et_xc // mohan add 2008-06-01 - Parallel_Reduce::reduce_double_pool( etxc ); - Parallel_Reduce::reduce_double_pool( vtxc ); + Parallel_Reduce::reduce_double_pool( et_xc ); + Parallel_Reduce::reduce_double_pool( vt_xc ); - etxc *= omega / ncxyz; - vtxc *= omega / ncxyz; + et_xc *= omega / ncxyz; + vt_xc *= omega / ncxyz; timer::tick("H_XC_pw","v_xc"); - return; + return std::make_tuple(et_xc,vt_xc,v); } diff --git a/ABACUS.develop/source/src_pw/H_XC_pw.h b/ABACUS.develop/source/src_pw/H_XC_pw.h index a8f6dc0d00..bef5b032d7 100644 --- a/ABACUS.develop/source/src_pw/H_XC_pw.h +++ b/ABACUS.develop/source/src_pw/H_XC_pw.h @@ -25,13 +25,12 @@ class H_XC_pw static double vtxc; // compute the exchange-correlation energy - static void v_xc( + static std::tuple v_xc( const int &nrxx, // number of real-space grid const int &ncxyz, // total number of charge grid const double &omega, // volume of cell - double** rho_in, - double *rho_core, // core charge density - matrix &v); + const double*const*const rho_in, + const double*const rho_core); // core charge density }; diff --git a/ABACUS.develop/source/src_pw/forces.cpp b/ABACUS.develop/source/src_pw/forces.cpp index 2b8d9e46e3..d4c47c09d5 100644 --- a/ABACUS.develop/source/src_pw/forces.cpp +++ b/ABACUS.develop/source/src_pw/forces.cpp @@ -551,8 +551,10 @@ void Forces::cal_force_ew(matrix& forceion) void Forces::cal_force_cc(matrix& forcecc) { // recalculate the exchange-correlation potential. - matrix vxc(NSPIN, pw.nrxx); - H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core, vxc); + const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core); + H_XC_pw::etxc = std::get<0>(etxc_vtxc_v); // may delete? + H_XC_pw::vtxc = std::get<1>(etxc_vtxc_v); // may delete? + const matrix vxc = std::get<2>(etxc_vtxc_v); complex * psiv = new complex [pw.nrxx]; ZEROS(psiv, pw.nrxx); diff --git a/ABACUS.develop/source/src_pw/potential.cpp b/ABACUS.develop/source/src_pw/potential.cpp index 770d110756..a5635ee115 100644 --- a/ABACUS.develop/source/src_pw/potential.cpp +++ b/ABACUS.develop/source/src_pw/potential.cpp @@ -292,15 +292,13 @@ void Potential::v_of_rho //---------------------------------------------------------- #ifdef USE_LIBXC - 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); + const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core); #endif + 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); //---------------------------------------------------------- // calculate the Hartree potential diff --git a/ABACUS.develop/source/src_pw/stress_func_cc.cpp b/ABACUS.develop/source/src_pw/stress_func_cc.cpp index 0a09c7f13d..b8e8c359fb 100644 --- a/ABACUS.develop/source/src_pw/stress_func_cc.cpp +++ b/ABACUS.develop/source/src_pw/stress_func_cc.cpp @@ -34,8 +34,10 @@ void Stress_Func::stress_cc(matrix& sigma, const bool is_pw) } //recalculate the exchange-correlation potential - matrix vxc(NSPIN, pw.nrxx); - H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core, vxc); + const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core); + H_XC_pw::etxc = std::get<0>(etxc_vtxc_v); // may delete? + H_XC_pw::vtxc = std::get<1>(etxc_vtxc_v); // may delete? + const matrix vxc = std::get<2>(etxc_vtxc_v); complex * psic = new complex [pw.nrxx]; From 8432362b3ab5657367ba67bf0d0cdd7e5bfcc0e6 Mon Sep 17 00:00:00 2001 From: linpz Date: Thu, 6 May 2021 23:17:14 +0800 Subject: [PATCH 3/3] 1. refactor Potential::v_of_rho() --- ABACUS.develop/source/src_lcao/ELEC_scf.cpp | 10 +++++----- ABACUS.develop/source/src_pw/efield.cpp | 4 ++-- ABACUS.develop/source/src_pw/efield.h | 4 ++-- ABACUS.develop/source/src_pw/electrons.cpp | 4 ++-- ABACUS.develop/source/src_pw/potential.cpp | 22 ++++++++++----------- ABACUS.develop/source/src_pw/potential.h | 4 +++- ABACUS.develop/source/src_pw/sto_elec.cpp | 4 ++-- 7 files changed, 26 insertions(+), 26 deletions(-) diff --git a/ABACUS.develop/source/src_lcao/ELEC_scf.cpp b/ABACUS.develop/source/src_lcao/ELEC_scf.cpp index b65d190f66..7cea8f06c5 100644 --- a/ABACUS.develop/source/src_lcao/ELEC_scf.cpp +++ b/ABACUS.develop/source/src_lcao/ELEC_scf.cpp @@ -189,7 +189,7 @@ void ELEC_scf::scf(const int &istep) // so be careful here, make sure // rho1 and rho2 are the same rho. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); en.delta_escf(); if (vext == 0) { @@ -417,7 +417,7 @@ void ELEC_scf::scf(const int &istep) if(!conv_elec) { // option 1 - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); en.delta_escf(); // option 2 @@ -426,16 +426,16 @@ void ELEC_scf::scf(const int &istep) // use real E_tot functional. //------------------------------ /* - pot.v_of_rho(CHR.rho_save, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho_save, CHR.rho); en.calculate_etot(); en.print_etot(conv_elec, istep, iter, dr2, 0.0, ETHR, avg_iter,0); - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); en.delta_escf(); */ } else { - pot.v_of_rho(CHR.rho, pot.vnew); + pot.vnew = pot.v_of_rho(CHR.rho, CHR.rho_core); //(used later for scf correction to the forces ) pot.vnew -= pot.vr; en.descf = 0.0; diff --git a/ABACUS.develop/source/src_pw/efield.cpp b/ABACUS.develop/source/src_pw/efield.cpp index be6140706e..1ed3ecbeb6 100644 --- a/ABACUS.develop/source/src_pw/efield.cpp +++ b/ABACUS.develop/source/src_pw/efield.cpp @@ -13,7 +13,7 @@ double Efield::etotefield; // energy correction due to the field double Efield::bvec[3]; double Efield::bmod; -void Efield::add_efield(double* rho, double* v_in) +void Efield::add_efield(const double*const rho, double* v_in) { TITLE("Efield","add_efield"); @@ -207,7 +207,7 @@ void Efield::add_efield(double* rho, double* v_in) } void Efield::compute_el_dip(const double &emaxpos, const double &eopreg, - const int &edir, const double *rho, double &e_dipole)const + const int &edir, const double*const rho, double &e_dipole)const { TITLE("Efield","compute_el_dip"); diff --git a/ABACUS.develop/source/src_pw/efield.h b/ABACUS.develop/source/src_pw/efield.h index 2bbb5433c1..98a6c9a396 100644 --- a/ABACUS.develop/source/src_pw/efield.h +++ b/ABACUS.develop/source/src_pw/efield.h @@ -10,7 +10,7 @@ class Efield Efield(); ~Efield(); - void add_efield(double* rho, double* v_in); + void add_efield(const double*const rho, double* v_in); // efield. static int edir; // direction of the field @@ -26,7 +26,7 @@ class Efield private: void compute_el_dip(const double& emaxpos, const double &eopreg, - const int &edir, const double* rho, double& e_dipole)const; + const int &edir, const double*const rho, double& e_dipole)const; void compute_ion_dip(const double& emaxpos, const double& eoperg, const int &edir, double& ion_dipole, const double& bmod, const double* bg)const; diff --git a/ABACUS.develop/source/src_pw/electrons.cpp b/ABACUS.develop/source/src_pw/electrons.cpp index bccf567c75..34e7d8a7a6 100644 --- a/ABACUS.develop/source/src_pw/electrons.cpp +++ b/ABACUS.develop/source/src_pw/electrons.cpp @@ -286,7 +286,7 @@ void Electrons::self_consistent(const int &istep) if (!conv_elec) { // not converged yet, calculate new potential from mixed charge density - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); // because = are calculated after sum // band, using output charge density. @@ -308,7 +308,7 @@ void Electrons::self_consistent(const int &istep) // mohan fix bug 2012-06-05, // the new potential V(PL)+V(H)+V(xc) - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); //cout<<"Exc = "<v_of_rho(CHR.rho, vr); + this->vr = this->v_of_rho(CHR.rho, CHR.rho_core); //---------------------------------------------------------- // (4) total potentials @@ -276,17 +276,15 @@ void Potential::set_local_pot( // The XC potential is computed in real space, while the // Hartree potential is computed in reciprocal space. //========================================================== -void Potential::v_of_rho -( - double **rho_in, - matrix &v_in -) +matrix Potential::v_of_rho( + const double*const*const rho_in, + const double * const rho_core_in) { TITLE("Potential","v_of_rho"); - v_in.zero_out(); - timer::tick("Potential","v_of_rho",'E'); + matrix v(NSPIN,pw.nrxx); + //---------------------------------------------------------- // calculate the exchange-correlation potential //---------------------------------------------------------- @@ -298,12 +296,12 @@ void Potential::v_of_rho #endif 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); + v += std::get<2>(etxc_vtxc_v); //---------------------------------------------------------- // calculate the Hartree potential //---------------------------------------------------------- - v_in += H_Hartree_pw::v_hartree(ucell, pw, NSPIN, rho_in); + v += H_Hartree_pw::v_hartree(ucell, pw, NSPIN, rho_in); // mohan add 2011-06-20 if(EFIELD && DIPOLE) @@ -311,11 +309,11 @@ void Potential::v_of_rho Efield EFID; for (int is = 0;is < NSPIN;is++) { - EFID.add_efield(rho_in[is], &v_in.c[is*pw.nrxx]); + EFID.add_efield(rho_in[is], &v.c[is*pw.nrxx]); } } timer::tick("Potential","v_of_rho",'E'); - return; + return v; } //end subroutine v_of_rho diff --git a/ABACUS.develop/source/src_pw/potential.h b/ABACUS.develop/source/src_pw/potential.h index 3580c2a668..f1ac36974c 100644 --- a/ABACUS.develop/source/src_pw/potential.h +++ b/ABACUS.develop/source/src_pw/potential.h @@ -41,7 +41,9 @@ class Potential ComplexMatrix &sf // structure factors ); - void v_of_rho(double** rho_in, matrix &v); + matrix v_of_rho( + const double*const*const rho_in, + const double * const rho_core_in); void set_vr_eff(void); diff --git a/ABACUS.develop/source/src_pw/sto_elec.cpp b/ABACUS.develop/source/src_pw/sto_elec.cpp index 8cd94ffd0e..c5bb04ea8c 100644 --- a/ABACUS.develop/source/src_pw/sto_elec.cpp +++ b/ABACUS.develop/source/src_pw/sto_elec.cpp @@ -215,7 +215,7 @@ void Stochastic_Elec::scf_stochastic(const int &istep) if (!conv_elec) { // not converged yet, calculate new potential from mixed charge density - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); // because = are calculated after sum // band, using output charge density. @@ -235,7 +235,7 @@ void Stochastic_Elec::scf_stochastic(const int &istep) } // the new potential V(PL)+V(H)+V(xc) - pot.v_of_rho(CHR.rho, pot.vr); + pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core); //( vnew used later for scf correction to the forces ) pot.vnew = pot.vr - pot.vnew;