Skip to content

Commit

Permalink
Merge branch 'develop' into dpmd
Browse files Browse the repository at this point in the history
  • Loading branch information
@dyzheng
dyzheng authored Oct 19, 2023
2 parents aca0d2a + 1d7e79f commit 0ef8ccf
Show file tree
Hide file tree
Showing 7 changed files with 120 additions and 33 deletions.
2 changes: 1 addition & 1 deletion deps/libpaw_interface
Submodule libpaw_interface updated 3 files
+0 −1 test/0_libpaw_mod.f90
+1 −1 test/get_dij.f90
+5 −3 test/prepare_libpaw.f90
6 changes: 6 additions & 0 deletions source/module_cell/module_paw/paw_cell.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,12 +219,16 @@ class Paw_Cell
int get_libpaw_ixc() {return ixc;}
int get_libpaw_xclevel() {return xclevel;}
int get_nspin() {return nspden;}

double get_epawdc() {return epawdc * 2.0;}

int get_nfft() {return nfft;}
int get_nrxx() {return nx*ny*num_z[GlobalV::RANK_IN_POOL];}
int get_val(const int it) {return paw_element_list[it].get_zval();}
int get_zat(const int it) {return paw_element_list[it].get_zat();}

double calculate_ecore();

private:
// Info to be passed to libpaw_interface:
// 1. ecut, ecutpaw : kinetic energy cutoff of the planewave basis set
Expand All @@ -242,6 +246,7 @@ class Paw_Cell

double ecut, ecutpaw;
std::vector<double> rprimd, gprimd, gmet;
std::vector<double> epsatm;
double ucvol;
std::vector<int> ngfft, ngfftdg;
int nfft;
Expand All @@ -251,6 +256,7 @@ class Paw_Cell
char* filename_list;
int xclevel, ixc;
int nspden, nsppol;
double epawdc;

// Part IV. Calling Fortran subroutines from libpaw_interface
public:
Expand Down
26 changes: 20 additions & 6 deletions source/module_cell/module_paw/paw_cell_libpaw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,7 @@ void Paw_Cell::set_libpaw_atom(const int natom_in, const int ntypat_in, const in

typat.resize(natom);
xred.resize(3 * natom);
epsatm.resize(ntypat);

for(int iat = 0; iat < natom; iat ++)
{
Expand Down Expand Up @@ -250,7 +251,7 @@ extern "C"
{
void prepare_libpaw_(double&,double&,double*,double*,double*,double&,int*,int*,
// ecut ecutpaw gmet rprimd gprimd ucvol ngfft ngfftdg
int&,int&,int*,double*,int&,int&,char*,int&,int&);
int&,int&,int*,double*,int&,int&,char*,int&,int&,double*);
// natom ntypat typat xred ixc xclevel filename_list nspden nsppol

void get_vloc_ncoret_(int*, int&,int&, int&, double*,double*,double*,double&,double*,double*,double*);
Expand All @@ -262,8 +263,8 @@ extern "C"
void get_nhat_(int&, int&, double*, int*, int&, int&, double*, double*, double&, double*, double*);
// natom,ntypat,xred, ngfft,nfft,nspden,gprimd, rprimd, ucvol, nhat, nhatgr

void calculate_dij_(int&, int&, int&, int&, int&, int&, double*, double&, double*, double*, double*);
// natom,ntypat,ixc, xclevel,nfft, nspden,xred, ucvol, gprimd, vks, vxc
void calculate_dij_(int&, int&, int&, int&, int&, int&, double*, double&, double*, double*, double*, double&);
// natom,ntypat,ixc, xclevel,nfft, nspden,xred, ucvol, gprimd, vks, vxc, epawdc

void get_dij_(int&, int&, int&, double*);
// iatom,size_dij,nspden,dij
Expand All @@ -278,7 +279,7 @@ void Paw_Cell::prepare_paw()
{
prepare_libpaw_(ecut, ecutpaw, gmet.data(), rprimd.data(), gprimd.data(), ucvol,
ngfft.data(), ngfftdg.data(), natom, ntypat, typat.data(), xred.data(),
ixc, xclevel, filename_list, nspden, nsppol);
ixc, xclevel, filename_list, nspden, nsppol, epsatm.data());
}

void Paw_Cell::get_vloc_ncoret(double* vloc, double* ncoret)
Expand Down Expand Up @@ -459,7 +460,7 @@ void Paw_Cell::calculate_dij(double* vks, double* vxc)

if(GlobalV::RANK_IN_POOL == 0)
{
calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree);
calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree,epawdc);
}

if(GlobalV::RANK_IN_POOL == 0)
Expand Down Expand Up @@ -489,7 +490,7 @@ void Paw_Cell::calculate_dij(double* vks, double* vxc)
}
}
}
calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree);
calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree,epawdc);
delete[] vks_hartree;
delete[] vxc_hartree;
#endif
Expand Down Expand Up @@ -645,4 +646,17 @@ void Paw_Cell::set_sij()
delete[] sij_libpaw;
delete[] sij;
}
}

double Paw_Cell::calculate_ecore()
{
double charge = 0.0;
double esum = 0.0;
for(int ia = 0; ia < natom; ia ++)
{
const int it = atom_type[ia];
esum += epsatm[it];
charge += paw_element_list[it].get_zval();
}
return esum * charge / ucvol * 2.0;
}
3 changes: 2 additions & 1 deletion source/module_cell/module_paw/paw_sphbes.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,7 @@ void Paw_Element::spherical_bessel_function(const int l, const double xx,
// search for case("r=a*(exp(d*i)-1)") in m_pawpsp.F90, for example
double Paw_Element::simpson_integration(std::vector<double> & f) const
{
int ir=nr;
int ir=nr-1;
while (std::abs(f[ir]) < 1e-20)
{
ir=ir-1;
Expand Down Expand Up @@ -275,6 +275,7 @@ void Paw_Element::prepare_simpson_integration(const double r_for_intg, int & mes
meshsz = ir;
}

if(meshsz == nr) meshsz -= 1;
//std::cout << "meshsz : " << meshsz << std::endl;

//get simp_fact
Expand Down
66 changes: 48 additions & 18 deletions source/module_elecstate/elecstate_energy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,10 @@
#include "elecstate_getters.h"
#include "module_base/global_variable.h"
#include "module_base/parallel_reduce.h"
#ifdef USE_PAW
#include "module_hamilt_general/module_xc/xc_functional.h"
#include "module_hamilt_pw/hamilt_pwdft/global.h"
#endif

namespace elecstate
{
Expand Down Expand Up @@ -98,41 +102,67 @@ double ElecState::cal_delta_eband() const
const double* v_eff = this->pot->get_effective_v(0);
const double* v_fixed = this->pot->get_fixed_v();
const double* v_ofk = nullptr;
if (get_xc_func_type() == 3 || get_xc_func_type() == 5)

#ifdef USE_PAW
if(GlobalV::use_paw)
{
v_ofk = this->pot->get_effective_vofk(0);
ModuleBase::matrix v_xc;
const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
= XC_Functional::v_xc(this->charge->nrxx, this->charge, &GlobalC::ucell);
v_xc = std::get<2>(etxc_vtxc_v);

for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
{
deband_aux -= this->charge->rho[0][ir] * v_xc(0,ir);
}
if (GlobalV::NSPIN == 2)
{
for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
{
deband_aux -= this->charge->rho[1][ir] * v_xc(1,ir);
}
}
}
#endif

for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
if(!GlobalV::use_paw)
{
deband_aux -= this->charge->rho[0][ir] * (v_eff[ir] - v_fixed[ir]);
if (get_xc_func_type() == 3 || get_xc_func_type() == 5)
{
deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir];
v_ofk = this->pot->get_effective_vofk(0);
}
}

if (GlobalV::NSPIN == 2)
{
v_eff = this->pot->get_effective_v(1);
v_ofk = this->pot->get_effective_vofk(1);
for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
{
deband_aux -= this->charge->rho[1][ir] * (v_eff[ir] - v_fixed[ir]);
deband_aux -= this->charge->rho[0][ir] * (v_eff[ir] - v_fixed[ir]);
if (get_xc_func_type() == 3 || get_xc_func_type() == 5)
{
deband_aux -= this->charge->kin_r[1][ir] * v_ofk[ir];
deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir];
}
}
}
else if (GlobalV::NSPIN == 4)
{
for (int is = 1; is < 4; is++)

if (GlobalV::NSPIN == 2)
{
v_eff = this->pot->get_effective_v(is);
v_eff = this->pot->get_effective_v(1);
v_ofk = this->pot->get_effective_vofk(1);
for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
{
deband_aux -= this->charge->rho[is][ir] * v_eff[ir];
deband_aux -= this->charge->rho[1][ir] * (v_eff[ir] - v_fixed[ir]);
if (get_xc_func_type() == 3 || get_xc_func_type() == 5)
{
deband_aux -= this->charge->kin_r[1][ir] * v_ofk[ir];
}
}
}
else if (GlobalV::NSPIN == 4)
{
for (int is = 1; is < 4; is++)
{
v_eff = this->pot->get_effective_v(is);
for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
{
deband_aux -= this->charge->rho[is][ir] * v_eff[ir];
}
}
}
}
Expand Down
43 changes: 40 additions & 3 deletions source/module_elecstate/fp_energy.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,8 @@
#include "fp_energy.h"
#include "module_base/global_variable.h"
#ifdef USE_PAW
#include "module_cell/module_paw/paw_cell.h"
#endif

#include <iomanip>
#include <iostream>
Expand All @@ -10,16 +14,49 @@ namespace elecstate
/// @brief calculate etot
double fenergy::calculate_etot()
{
etot = eband + deband + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield + gatefield
+ evdw + esol_el + esol_cav + edftu + edeepks_scf;
if(GlobalV::use_paw)
{
etot = eband + deband + etxc + ewald_energy - hartree_energy + demet + descf + exx + efield + gatefield
+ evdw + esol_el + esol_cav + edftu + edeepks_scf;
}
else
{
etot = eband + deband + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield + gatefield
+ evdw + esol_el + esol_cav + edftu + edeepks_scf;
}

#ifdef USE_PAW
if(GlobalV::use_paw)
{
double ecore = GlobalC::paw_cell.calculate_ecore();
double epawdc = GlobalC::paw_cell.get_epawdc();
etot += ( ecore + epawdc );
}
#endif
return etot;
}

/// @brief calculate etot_harris
double fenergy::calculate_harris()
{
etot_harris = eband + deband_harris + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield
if(GlobalV::use_paw)
{
etot_harris = eband + deband_harris + etxc + ewald_energy - hartree_energy + demet + descf + exx + efield
+ gatefield + evdw + esol_el + esol_cav + edftu + edeepks_scf;
}
else
{
etot_harris = eband + deband_harris + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield
+ gatefield + evdw + esol_el + esol_cav + edftu + edeepks_scf;
}
#ifdef USE_PAW
if(GlobalV::use_paw)
{
double ecore = GlobalC::paw_cell.calculate_ecore();
double epawdc = GlobalC::paw_cell.get_epawdc();
etot_harris += ( ecore + epawdc );
}
#endif
return etot_harris;
}

Expand Down
7 changes: 3 additions & 4 deletions tests/integrate/101_PW_15_paw/result.ref
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
etotref -261.3678632375260
etotperatomref -130.6839316188
etotref -200.9340475810059
etotperatomref -100.4670237905
pointgroupref T_d
spacegroupref O_h
nksibzref 1
totaltimeref 29.2697
1.07
totaltimeref 1.17

0 comments on commit 0ef8ccf

Please sign in to comment.