Permalink
Browse files

made xc response work again

  • Loading branch information...
fbischoff committed Aug 30, 2016
1 parent 9bfc507 commit 75b26fd84ab37eab493a9e5d2d1e312b3fc62ab7
Showing with 167 additions and 131 deletions.
  1. +1 −0 .gitignore
  2. +3 −3 src/apps/chem/SCF.cc
  3. +33 −63 src/apps/chem/SCFOperators.cc
  4. +69 −22 src/apps/chem/xcfunctional.h
  5. +61 −43 src/apps/chem/xcfunctional_libxc.cc
View
@@ -52,3 +52,4 @@ cmake/cmake-*
cmake/bin
cmake/doc
cmake/share
Makefile-prog.in
View
@@ -3424,17 +3424,17 @@ namespace madness {
world.gop.fence(); // NECESSARY
vf[XCfunctional::enum_saa] = copy( delrho[0] * delrho[0]
vf[XCfunctional::enum_saa] = copy( delrho[0] * delrho[0]
+ delrho[1] * delrho[1]
+ delrho[2] * delrho[2]); // sigma_aa
if (!param.spin_restricted && param.nbeta != 0)
{
vf[XCfunctional::enum_sab] = copy( delrho[0] * delrho[3]
vf[XCfunctional::enum_sab] = copy( delrho[0] * delrho[3]
+ delrho[1] * delrho[4]
+ delrho[2] * delrho[5]); // sigma_ab
vf[XCfunctional::enum_sbb] = copy( delrho[3] * delrho[3]
vf[XCfunctional::enum_sbb] = copy( delrho[3] * delrho[3]
+ delrho[4] * delrho[4]
+ delrho[5] * delrho[5]); // sigma_bb
}
@@ -659,7 +659,7 @@ real_function_3d XCOperator::make_xc_potential() const {
/// + 2\frac{\partial^2 f_{xc}}{\partial \rho_\alpha\sigma_{\alpha\alpha}}
/// \left(\vec\nabla \rho_a\cdot \vec \nabla\tilde\rho\right)
/// \f]
// the second partial derivatives that need to be multiplied with the density gradients
/// the second partial derivatives that need to be multiplied with the density gradients
/// \f[
/// second_{semilocal} = -\vec\nabla\cdot\left((\vec\nabla\rho)
/// \left[2\frac{\partial^2 f_{xc}}{\partial\rho_\alpha\partial\sigma_{\alpha\alpha}}\tilde\rho
@@ -686,48 +686,30 @@ real_function_3d XCOperator::apply_xc_kernel(const real_function_3d& dens_pt) co
// compute the various terms from the xc kernel
// compute the local terms: second_{local}
// compute the local terms: second_{local} * rho
real_function_3d result=multiop_values<double, xc_kernel_apply, 3>
(xc_kernel_apply(*xc, ispin, XCfunctional::kernel_second_local), xc_args);
// save(result,"local_apply");
if (xc->is_gga()) {
const real_function_3d& arho=xc_args[XCfunctional::enum_rhoa];
vecfuncT zeta(3);
zeta[0]=xc_args[XCfunctional::enum_zetaa_x];
zeta[1]=xc_args[XCfunctional::enum_zetaa_y];
zeta[2]=xc_args[XCfunctional::enum_zetaa_z];
// compute the semilocal terms, second partial derivatives
real_function_3d semilocal2a=multiop_values<double, xc_kernel_apply, 3>
(xc_kernel_apply(*xc, ispin, XCfunctional::kernel_second_semilocal), xc_args);
save(semilocal2a,"semilocal2a");
const real_function_3d& rho=xc_args[XCfunctional::enum_rhoa];
double tol=xc->get_ggatol();
real_function_3d semilocal2=binary_op(semilocal2a,rho,binary_munging(tol));
// real_function_3d semilocal2=semilocal2a;
save(semilocal2,"semilocal2_ddel");
for (int idim=0; idim<3; ++idim) {
real_function_3d ddens = xc_args[XCfunctional::enum_drhoa_x+idim];
real_function_3d ddel = 2.0*semilocal2*ddens;
save(ddel,"ddel"+stringify(idim));
Derivative<double,3> D = free_space_derivative<double,3>(world, idim);
real_function_3d vxc2=D(ddel);
save(vxc2,"vxc2"+stringify(idim));
result-=vxc2;
}
real_function_3d semilocal2=binary_op(semilocal2a,arho,binary_munging(tol));
result-=2.0* div(zeta * (arho*semilocal2)); // factor 2 missing in arho
// compute the semilocal terms, first partial derivative
real_function_3d semilocal1a=multiop_values<double, xc_kernel_apply, 3>
(xc_kernel_apply(*xc, ispin, XCfunctional::kernel_first_semilocal), xc_args);
real_function_3d semilocal1=binary_op(semilocal1a,rho,binary_munging(tol));
// real_function_3d semilocal1=semilocal1a;
// save(semilocal1a,"semilocal1a");
for (int idim=0; idim<3; ++idim) {
real_function_3d ddel = semilocal1*ddens_pt[idim];
Derivative<double,3> D = free_space_derivative<double,3>(world, idim);
real_function_3d vxc2=D(ddel);
result-=vxc2;
}
real_function_3d result1=binary_op(result,rho,binary_munging(tol));
result=result1;
real_function_3d semilocal1=binary_op(semilocal1a,arho,binary_munging(tol));
result-=div(semilocal1*ddens_pt);
}
truncate(world,xc_args);
@@ -758,7 +740,6 @@ vecfuncT XCOperator::prep_xc_args(const real_function_3d& arho,
xcargs[XCfunctional::enum_zetaa_x]=grada[0];
xcargs[XCfunctional::enum_zetaa_y]=grada[1];
xcargs[XCfunctional::enum_zetaa_z]=grada[2];
// save(xcargs[XCfunctional::enum_chi_aa],"chi");
if (have_beta) {
real_function_3d logdensb=unary_op(brho,logme());
@@ -771,10 +752,6 @@ vecfuncT XCOperator::prep_xc_args(const real_function_3d& arho,
xcargs[XCfunctional::enum_chi_bb]=chib;
xcargs[XCfunctional::enum_chi_ab]=chiab;
}
// // this is needed for proper munging of sigma_ab
// xcargs[XCfunctional::enum_sigtot]= xcargs[XCfunctional::enum_saa]
// +2.0*xcargs[XCfunctional::enum_sab]+xcargs[XCfunctional::enum_sbb];
}
world.gop.fence();
@@ -789,46 +766,39 @@ void XCOperator::prep_xc_args_response(const real_function_3d& dens_pt,
const bool have_beta=(xc->is_spin_polarized()) and (nbeta>0);
World& world=dens_pt.world();
const real_function_3d& arho=xc_args[XCfunctional::enum_rhoa];
const real_function_3d& brho=xc_args[XCfunctional::enum_rhob];
ddens_pt=vecfuncT(3);
// assign the perturbed density (spin-free ??)
// assign the perturbed density (spin-free)
xc_args[XCfunctional::enum_rho_pt]=dens_pt;
world.gop.fence();
// assign the reduced density gradients with the perturbed density for GGA
// \sigma_pt = 2.0 * \nabla \rho_\alpha \cdot \nabla\tilde\rho
// \sigma_pt_a = \nabla \rho_\alpha \cdot \nabla\tilde\rho
// \sigma_pt_b = \nabla \rho_\beta \cdot \nabla\tilde\rho
//
// using the logarithmic derivatives for rho only we get (alpha and RHF)
// \sigma_pt = 2.0 * \rho_\alpha (\nabla\zeta_\alpha \cdot \nabla\tilde\rho)
// \sigma_pt_a = \rho_\alpha (\nabla\zeta_\alpha \cdot \nabla\tilde\rho)
// we save the functions without multiplying the ground state density rho
if (xc->is_gga()) {
std::vector< std::shared_ptr<Derivative<double,3> > > gradop =
gradient_operator<double,3>(world);
arho.reconstruct(false);
brho.reconstruct(false);
dens_pt.reconstruct(false);
world.gop.fence();
ddens_pt=grad(dens_pt); // spin free
std::vector<real_function_3d> ddensa(3),ddensb(3);
for (std::size_t i=0; i<3; ++i) {
ddensa[i]=(*gradop[i])(arho, false);
if (have_beta) ddensb[i]=(*gradop[i])(brho, false);
ddens_pt[i]=(*gradop[i])(dens_pt, false);
}
std::vector<real_function_3d> zeta(3);
zeta[0]=xc_args[XCfunctional::enum_zetaa_x];
zeta[1]=xc_args[XCfunctional::enum_zetaa_y];
zeta[2]=xc_args[XCfunctional::enum_zetaa_z];
xc_args[XCfunctional::enum_sigma_pta_div_rho]=dot(world,zeta,ddens_pt); // sigma_a
// for RHF add factor 2 on rho; will be done in xcfunctional_libxc::make_libxc_args
// \sigma_pt = 2 * rho_a * sigma_pta_div_rho
world.gop.fence();
// autorefine before squaring
for (std::size_t i=0; i<3; ++i) {
ddensa[i].refine(false);
if (have_beta) ddensb[i].refine(false);
ddens_pt[i].refine(false);
if (have_beta) {
zeta[0]=xc_args[XCfunctional::enum_zetab_x];
zeta[1]=xc_args[XCfunctional::enum_zetab_y];
zeta[2]=xc_args[XCfunctional::enum_zetab_z];
xc_args[XCfunctional::enum_sigma_ptb_div_rho]= dot(world,zeta,ddens_pt); // sigma_b
}
world.gop.fence();
xc_args[XCfunctional::enum_sigma_pta]=dot(world,ddensa,ddens_pt); // sigma_a
if (have_beta) xc_args[XCfunctional::enum_sigma_ptb]= dot(world,ddensb,ddens_pt);
// save(xc_args[XCfunctional::enum_sigma_pta],"sigma_pt");
world.gop.fence();
}
world.gop.fence();
truncate(world,xc_args,extra_truncation);
@@ -43,11 +43,11 @@ struct xc_lda_potential {
class XCfunctional {
public:
/// the ordering of the intermediates is fixed, but the code can handle
/// The ordering of the intermediates is fixed, but the code can handle
/// non-initialized functions, so if e.g. no GGA is requested, all the
/// corresponding vector components may be left empty.
///
/// note the additional quantity zeta, which is defined as
/// Note the additional quantities \f$ \zeta \f$ and \f$ \chi \f$, which are defined as
/// \f$
/// \rho = exp(\zeta)
/// \f$
@@ -57,33 +57,48 @@ class XCfunctional {
/// \f$
/// The reduced gradients \sigma may then be expressed as
/// \f$
/// \sigma = |\nabla\rho|^2 = |\rho|^2 |\nabla\zeta|^2 = |\rho|^2 chi
/// \sigma = |\nabla\rho|^2 = |\rho|^2 |\nabla\zeta|^2 = |\rho|^2 \chi
/// \f$
///
/// For conversion to spin-restricted calculations add a factor of 2 for all
/// quantities with 1 spin index and a factor 4 for all quantities with 2
/// spin indices. Note the perturbed density is spin-independent, and the
/// perturbed sigma depends only on 1 spin index. The perturbed sigma is
/// reconstructed from the density and the rest
/// \f[
/// \sigma_{pt,\alpha} = \nabla \rho_\alpha \cdot \nabla \rho_{pt}
/// = \rho \left(\zeta_\alpha \cdot \nabla\rho_{pt}\right)
/// \f]
/// and for RHF
/// \f[
/// \sigma_{pt} = 2.0 * \nabla \rho_\alpha \cdot \nabla \rho_{pt}
/// = 2.0 * \rho \left(\zeta_\alpha \cdot \nabla\rho_{pt}\right)
/// \f]
enum xc_arg {
enum_rhoa=0, ///< alpha density \f$ \rho_\alpha \f$
enum_rhob=1, ///< \f$ \rho_\beta \f$
enum_rhob=1, ///< beta density \f$ \rho_\beta \f$
enum_rho_pt=2, ///< perturbed density (CPHF, TDKS) \f$ \rho_{pt} \f$
enum_drhoa_x=4, ///< \f$ \partial/{\partial x} \rho_{\alpha} \f$
enum_drhoa_y=5, ///< \f$ \partial/{\partial y} \rho_{\alpha} \f$
enum_drhoa_z=6, ///< \f$ \partial/{\partial z} \rho_{\alpha} \f$
enum_drhob_x=7, ///< \f$ \partial/{\partial x} \rho_{\beta} \f$
enum_drhob_y=8, ///< \f$ \partial/{\partial y} \rho_{\beta} \f$
enum_drhob_z=9, ///< \f$ \partial/{\partial z} \rho_{\beta} \f$
// enum_drhoa_x=4, ///< \f$ \partial/{\partial x} \rho_{\alpha} \f$
// enum_drhoa_y=5, ///< \f$ \partial/{\partial y} \rho_{\alpha} \f$
// enum_drhoa_z=6, ///< \f$ \partial/{\partial z} \rho_{\alpha} \f$
// enum_drhob_x=7, ///< \f$ \partial/{\partial x} \rho_{\beta} \f$
// enum_drhob_y=8, ///< \f$ \partial/{\partial y} \rho_{\beta} \f$
// enum_drhob_z=9, ///< \f$ \partial/{\partial z} \rho_{\beta} \f$
enum_saa=10, ///< \f$ \sigma_{aa} = \nabla \rho_{\alpha}.\nabla \rho_{\alpha} \f$
enum_sab=11, ///< \f$ \sigma_{ab} = \nabla \rho_{\alpha}.\nabla \rho_{\beta} \f$
enum_sbb=12, ///< \f$ \sigma_{bb} = \nabla \rho_{\beta}.\nabla \rho_{\beta} \f$
enum_sigtot=13, ///< \f$ \sigma = \nabla \rho.\nabla \rho \f$
enum_sigma_pta=14, ///< \f$ \nabla\rho_{\alpha}.\nabla\rho_{pt} \f$
enum_sigma_ptb=15, ///< \f$ \nabla\rho_{\beta}.\nabla\rho_{pt} \f$
enum_zetaa_x=16, ///< \f$ \partial/{\partial x} ln(rho_a) \f$
enum_zetaa_y=17, ///< \f$ \partial/{\partial y} ln(rho_a) \f$
enum_zetaa_z=18, ///< \f$ \partial/{\partial z} ln(rho_a) \f$
enum_zetab_x=19, ///< \f$ \partial/{\partial x} ln(rho_b) \f$
enum_zetab_y=20, ///< \f$ \partial/{\partial y} ln(rho_b) \f$
enum_zetab_z=21, ///< \f$ \partial/{\partial z} ln(rho_b) \f$
enum_chi_aa=22, ///< \f$ \nabla \zeta{\alpha}.\nabla \zeta{\alpha} \f$
enum_chi_ab=23, ///< \f$ \nabla \zeta{\alpha}.\nabla \zeta{\beta} \f$
enum_chi_bb=24 ///< \f$ \nabla \zeta{\beta}.\nabla \zeta{\beta} \f$
enum_sigma_pta_div_rho=14, ///< \f$ \zeta_{\alpha}.\nabla\rho_{pt} \f$
enum_sigma_ptb_div_rho=15, ///< \f$ \zeta_{\beta}.\nabla\rho_{pt} \f$
enum_zetaa_x=16, ///< \f$ \zeta_{a,x}=\partial/{\partial x} \ln(\rho_a) \f$
enum_zetaa_y=17, ///< \f$ \zeta_{a,y}=\partial/{\partial y} \ln(\rho_a) \f$
enum_zetaa_z=18, ///< \f$ \zeta_{a,z}=\partial/{\partial z} \ln(\rho_a) \f$
enum_zetab_x=19, ///< \f$ \zeta_{b,x} = \partial/{\partial x} \ln(\rho_b) \f$
enum_zetab_y=20, ///< \f$ \zeta_{b,y} = \partial/{\partial y} \ln(\rho_b) \f$
enum_zetab_z=21, ///< \f$ \zeta_{b,z} = \partial/{\partial z} \ln(\rho_b) \f$
enum_chi_aa=22, ///< \f$ \chi_{aa} = \nabla \zeta_{\alpha}.\nabla \zeta_{\alpha} \f$
enum_chi_ab=23, ///< \f$ \chi_{ab} = \nabla \zeta_{\alpha}.\nabla \zeta_{\beta} \f$
enum_chi_bb=24 ///< \f$ \chi_{bb} = \nabla \zeta_{\beta}.\nabla \zeta_{\beta} \f$
};
const static int number_xc_args=25; ///< max number of intermediates
@@ -120,10 +135,31 @@ class XCfunctional {
#ifdef MADNESS_HAS_LIBXC
std::vector< std::pair<xc_func_type*,double> > funcs;
/// convert the raw density (gradient) data to be used by the xc operators
/// involves mainly munging of the densities
/// response densities and density gradients are munged based on the
/// value of the ground state density, since they may become negative
/// and may also be much more diffuse.
/// dimensions of the output tensors are for spin-restricted and unrestricted
/// (with np the number of grid points in the box):
/// rho(np) or rho(2*np)
/// sigma(np) sigma(3*np)
/// rho_pt(np)
/// sigma_pt(2*np)
/// @param[in] t input density (gradients)
/// @param[out] rho ground state (spin) density, properly munged
/// @param[out] sigma ground state (spin) density gradients, properly munged
/// @param[out] rho_pt response density, properly munged (no spin)
/// @param[out] sigma_pt response (spin) density gradients, properly munged
/// @param[in] need_response flag if rho_pt and sigma_pt need to be calculated
void make_libxc_args(const std::vector< madness::Tensor<double> >& t,
madness::Tensor<double>& rho,
madness::Tensor<double>& sigma,
const munging_type& munging) const;
madness::Tensor<double>& rho_pt,
madness::Tensor<double>& sigma_pt,
const bool need_response) const;
void make_libxc_args_old(const std::vector< madness::Tensor<double> >& t,
madness::Tensor<double>& rho,
madness::Tensor<double>& sigma,
@@ -216,6 +252,17 @@ class XCfunctional {
return rho;
}
/// munge rho if refrho is small
/// special case for perturbed densities, which might be negative and diffuse.
/// Munge rho (e.g. the perturbed density) if the reference density refrho
/// e.g. the ground state density is small. Only where the reference density
/// is large enough DFT is numerically well-defined.
double binary_munge(double rho, double refrho) const {
if (refrho<rhotol) rho=rhomin;
return rho;
}
/// similar to the Laura's ratio thresholding, but might be more robust
void munge_xc_kernel(double& rho, double& sigma) const {
if (sigma<0.0) sigma=sigmin;
Oops, something went wrong.

0 comments on commit 75b26fd

Please sign in to comment.