Skip to content

Commit

Permalink
Merge branch 'non_collinear_XC' into 'master'
Browse files Browse the repository at this point in the history
main modifications needed to implement non collinear XC functional

Closes #152 and #157

See merge request npneq/inq!1082
  • Loading branch information
xavierandrade committed Sep 28, 2024
2 parents 0b3d750 + a003539 commit 3eee860
Show file tree
Hide file tree
Showing 6 changed files with 317 additions and 35 deletions.
8 changes: 4 additions & 4 deletions src/hamiltonian/scalar_potential.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ void scalar_potential_add(basis::field_set<basis::real_space, PotentialType> con
gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(),
[pot = begin(potential.matrix()), vph = begin(vphi.spinor_array()), ph = begin(phi.spinor_array()), shift] GPU_LAMBDA (auto ist, auto ip){
auto offdiag = pot[ip][2] + complex{0.0, 1.0}*pot[ip][3];
vph[ip][0][ist] += (pot[ip][0] + shift)*ph[ip][0][ist] + offdiag*ph[ip][1][ist];
vph[ip][1][ist] += (pot[ip][1] + shift)*ph[ip][1][ist] + conj(offdiag)*ph[ip][0][ist];
vph[ip][0][ist] += (pot[ip][0] + shift)*ph[ip][0][ist] + conj(offdiag)*ph[ip][1][ist];
vph[ip][1][ist] += (pot[ip][1] + shift)*ph[ip][1][ist] + offdiag*ph[ip][0][ist];
});
}

Expand Down Expand Up @@ -68,8 +68,8 @@ states::orbital_set<basis::real_space, complex> scalar_potential(basis::field_se
gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(),
[pot = begin(potential.matrix()), vph = begin(vphi.spinor_array()), ph = begin(phi.spinor_array()), shift] GPU_LAMBDA (auto ist, auto ip){
auto offdiag = pot[ip][2] + complex{0.0, 1.0}*pot[ip][3];
vph[ip][0][ist] = (pot[ip][0] + shift)*ph[ip][0][ist] + offdiag*ph[ip][1][ist];
vph[ip][1][ist] = (pot[ip][1] + shift)*ph[ip][1][ist] + conj(offdiag)*ph[ip][0][ist];
vph[ip][0][ist] = (pot[ip][0] + shift)*ph[ip][0][ist] + conj(offdiag)*ph[ip][1][ist];
vph[ip][1][ist] = (pot[ip][1] + shift)*ph[ip][1][ist] + offdiag*ph[ip][0][ist];
});

}
Expand Down
4 changes: 4 additions & 0 deletions src/hamiltonian/xc_functional.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,10 @@ namespace hamiltonian {
}
return refs;
}

auto & nspin() const {
return nspin_;
}

private:

Expand Down
149 changes: 133 additions & 16 deletions src/hamiltonian/xc_term.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,16 @@
#include <basis/field.hpp>
#include <solvers/poisson.hpp>
#include <observables/density.hpp>
#include <observables/magnetization.hpp>
#include <operations/add.hpp>
#include <operations/integral.hpp>
#include <options/theory.hpp>
#include <hamiltonian/xc_functional.hpp>
#include <hamiltonian/atomic_potential.hpp>
#include <perturbations/none.hpp>
#include <utils/profiling.hpp>
#include <systems/electrons.hpp>
#include <cmath>

namespace inq {
namespace hamiltonian {
Expand All @@ -30,8 +33,8 @@ class xc_term {
public:

xc_term(options::theory interaction, int const spin_components){
functionals_.emplace_back(int(interaction.exchange()), spin_components);
functionals_.emplace_back(int(interaction.correlation()), spin_components);
functionals_.emplace_back(int(interaction.exchange()), std::min(spin_components, 2));
functionals_.emplace_back(int(interaction.correlation()), std::min(spin_components, 2));
}

////////////////////////////////////////////////////////////////////////////////////////////
Expand All @@ -56,17 +59,16 @@ class xc_term {
SpinDensityType full_density(spin_density.basis(), std::min(2, spin_density.set_size()));

if(spin_density.set_size() == 4) {
//for spinors convert the density to 2 components
gpu::run(spin_density.basis().local_size(),
[spi = begin(spin_density.matrix()), ful = begin(full_density.matrix()), cor = begin(core_density.linear())] GPU_LAMBDA (auto ip){
auto dtot = spi[ip][0] + spi[ip][1];
auto dd = spi[ip][0] - spi[ip][1];
auto dpol = sqrt(dd*dd + 4.0*(spi[ip][2]*spi[ip][2] + spi[ip][3]*spi[ip][3]));
auto mag = observables::local_magnetization(spi[ip], 4);
auto dpol = mag.length();
ful[ip][0] = 0.5*(dtot + dpol);
ful[ip][1] = 0.5*(dtot - dpol);
for(int ispin = 0; ispin < 2; ispin++){
if(ful[ip][ispin] < 0.0) ful[ip][ispin] = 0.0;
ful[ip][ispin] += 0.5*cor[ip];
ful[ip][ispin] += cor[ip]/2;
}
});
} else {
Expand All @@ -87,15 +89,81 @@ class xc_term {

////////////////////////////////////////////////////////////////////////////////////////////

template <typename VXC, typename VKS>
void process_potential(VXC const & vxc, VKS & vks) const {
template <typename SpinDensityType, typename VXC, typename VKS>
void process_potential(SpinDensityType const & spin_density, VXC const & vxc, VKS & vks) const {

gpu::run(vxc.local_set_size(), vxc.basis().local_size(),
if (spin_density.set_size() == 4) {
gpu::run(vxc.basis().local_size(),
[spi = begin(spin_density.matrix()), vx = begin(vxc.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto ip){
auto bxc = 0.5*(vx[ip][0] - vx[ip][1]);
auto vxc = 0.5*(vx[ip][0] + vx[ip][1]);
auto mag = observables::local_magnetization(spi[ip], 4);
auto dpol = mag.length();
if (fabs(dpol) > 1.e-7) {
auto e_mag = mag/dpol;
vk[ip][0] += vxc + bxc*e_mag[2];
vk[ip][1] += vxc - bxc*e_mag[2];
vk[ip][2] += bxc*e_mag[0];
vk[ip][3] += bxc*e_mag[1];
}
else {
vk[ip][0] += vxc;
vk[ip][1] += vxc;
}
});
}
else {
assert(spin_density.set_size() == 1 or spin_density.set_size() == 2);
gpu::run(vxc.local_set_size(), vxc.basis().local_size(),
[vx = begin(vxc.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip){
vk[ip][is] += vx[ip][is];
vk[ip][is] += vx[ip][is];
});
}

}


///////////////////////////////////////////////////////////////////////////////////////////

template <typename SpinDensityType, typename VXC>
double compute_nvxc(SpinDensityType const & spin_density, VXC const & vfunc) const {

auto nvxc_ = 0.0;
if (spin_density.set_size() == 4) {
basis::field_set<basis::real_space, double> vxc(spin_density.skeleton());
vxc.fill(0.0);
gpu::run(vfunc.basis().local_size(),
[spi = begin(spin_density.matrix()), vxi = begin(vfunc.matrix()), vxf = begin(vxc.matrix())] GPU_LAMBDA (auto ip){
auto b = 0.5*(vxi[ip][0] - vxi[ip][1]);
auto v = 0.5*(vxi[ip][0] + vxi[ip][1]);
auto mag = observables::local_magnetization(spi[ip], 4);
auto dpol = mag.length();
if (fabs(dpol) > 1.e-7) {
auto e_mag = mag/dpol;
// Vxc = [vxc+bxc^z, bxc^-; bxc^+, vxc-bxc^z]
vxf[ip][0] = v + b*e_mag[2];
vxf[ip][1] = v - b*e_mag[2];
vxf[ip][2] = 2.0*b*e_mag[0];
vxf[ip][3] = 2.0*b*e_mag[1];
}
else {
vxf[ip][0] = v;
vxf[ip][1] = v;
}
});
basis::field<basis::real_space, double> rfield(spin_density.basis());
rfield.fill(0.0);
gpu::run(spin_density.local_set_size(), spin_density.basis().local_size(),
[spi = begin(spin_density.matrix()), vx = begin(vxc.matrix()), rf = begin(rfield.linear())] GPU_LAMBDA (auto is, auto ip){
rf[ip] += vx[ip][is] * spi[ip][is];
});
nvxc_ += operations::integral(rfield);
}
else {
nvxc_ = operations::integral_product_sum(spin_density, vfunc);
}
return nvxc_;
}

////////////////////////////////////////////////////////////////////////////////////////////

template <typename SpinDensityType, typename CoreDensityType, typename VKSType>
Expand All @@ -108,8 +176,8 @@ class xc_term {
auto full_density = process_density(spin_density, core_density);

double efunc = 0.0;
basis::field_set<basis::real_space, double> vfunc(spin_density.skeleton());


basis::field_set<basis::real_space, double> vfunc(full_density.skeleton());
auto density_gradient = std::optional<decltype(operations::gradient(full_density))>{};
if(any_requires_gradient()) density_gradient.emplace(operations::gradient(full_density));

Expand All @@ -118,8 +186,56 @@ class xc_term {

evaluate_functional(func, full_density, density_gradient, efunc, vfunc);
exc += efunc;
process_potential(vfunc, vks);
nvxc += operations::integral_product_sum(spin_density, vfunc); //the core correction does not go here
process_potential(spin_density, vfunc, vks);

nvxc += compute_nvxc(spin_density, vfunc);
}
}

////////////////////////////////////////////////////////////////////////////////////////////

template <typename SpinDensityType, typename CoreDensityType, typename VxcType>
void compute_vxc(SpinDensityType const & spin_density, CoreDensityType const & core_density, VxcType & vxc) {

if(not any_true_functional()) return;
auto full_density = process_density(spin_density, core_density);
double efunc = 0.0;
basis::field_set<basis::real_space, double> vfunc(full_density.skeleton());
vfunc.fill(0.0);
vxc.fill(0.0);
auto density_gradient = std::optional<decltype(operations::gradient(full_density))>{};
if(any_requires_gradient()) density_gradient.emplace(operations::gradient(full_density));

for(auto & func : functionals_) {
if(not func.true_functional()) continue;
evaluate_functional(func, full_density, density_gradient, efunc, vfunc);
if (spin_density.set_size() == 4) {
gpu::run(vfunc.basis().local_size(),
[spi = begin(spin_density.matrix()), vxi = begin(vfunc.matrix()), vxf = begin(vxc.matrix())] GPU_LAMBDA (auto ip){
auto b = 0.5*(vxi[ip][0] - vxi[ip][1]);
auto v = 0.5*(vxi[ip][0] + vxi[ip][1]);
auto mag = observables::local_magnetization(spi[ip], 4);
auto dpol = mag.length();
if (fabs(dpol) > 1.e-7) {
auto e_mag = mag / dpol;
vxf[ip][0] += v + b*e_mag[2];
vxf[ip][1] += v - b*e_mag[2];
vxf[ip][2] += b*e_mag[0];
vxf[ip][3] += b*e_mag[1];
}
else {
vxf[ip][0] += v;
vxf[ip][1] += v;
}
});
}
else {
assert(spin_density.set_size() == 1 or spin_density.set_size() == 2);
gpu::run(vfunc.local_set_size(), vfunc.basis().local_size(),
[vxi = begin(vfunc.matrix()), vxf = begin(vxc.matrix())] GPU_LAMBDA (auto is, auto ip){
vxf[ip][is] += vxi[ip][is];
});
}
}
}

Expand All @@ -131,9 +247,10 @@ class xc_term {
CALI_CXX_MARK_FUNCTION;

auto edens = basis::field<basis::real_space, double>(density.basis());

assert(functional.nspin() == density.set_size());

if(functional.family() == XC_FAMILY_LDA){

xc_lda_exc_vxc(functional.libxc_func_ptr(), density.basis().local_size(), raw_pointer_cast(density.matrix().data_elements()),
raw_pointer_cast(edens.linear().data_elements()), raw_pointer_cast(vfunctional.matrix().data_elements()));
gpu::sync();
Expand Down
37 changes: 22 additions & 15 deletions src/observables/magnetization.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,27 @@
namespace inq {
namespace observables {

template <class Density>
GPU_FUNCTION auto local_magnetization(Density const & spin_density, int const & components) {
vector3<double> mag_density;

if(components == 4){
mag_density[0] = 2.0*spin_density[2];
mag_density[1] = 2.0*spin_density[3];
} else {
mag_density[0] = 0.0;
mag_density[1] = 0.0;
}

if(components >= 2){
mag_density[2] = spin_density[0] - spin_density[1];
} else {
mag_density[2] = 0.0;
}

return mag_density;
}

basis::field<basis::real_space, vector3<double>> magnetization(basis::field_set<basis::real_space, double> const & spin_density){

// The formula comes from here: https://gitlab.com/npneq/inq/-/wikis/Magnetization
Expand All @@ -24,21 +45,7 @@ basis::field<basis::real_space, vector3<double>> magnetization(basis::field_set<

gpu::run(magnet.basis().local_size(),
[mag = begin(magnet.linear()), den = begin(spin_density.matrix()), components = spin_density.set_size()] GPU_LAMBDA (auto ip){

if(components == 4){
mag[ip][0] = 2.0*den[ip][2];
mag[ip][1] = 2.0*den[ip][3];
} else {
mag[ip][0] = 0.0;
mag[ip][1] = 0.0;
}

if(components >= 2){
mag[ip][2] = den[ip][0] - den[ip][1];
} else {
mag[ip][2] = 0.0;
}

mag[ip] = local_magnetization(den[ip], components);
});

return magnet;
Expand Down
41 changes: 41 additions & 0 deletions tests/non_collinear_gs.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
/* -*- indent-tabs-mode: t -*- */

#include <inq/inq.hpp>
#include <stdio.h>

using namespace inq;
using namespace inq::magnitude;

inq::utils::match match(3.0e-4);

template <class EnvType>
auto compute_GS(EnvType const & env, systems::ions const & ions) {

auto electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized());
ground_state::initial_guess(ions, electrons);
auto result_col = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000));
auto mag = observables::total_magnetization(electrons.spin_density());

auto electrons_nc = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear());
ground_state::initial_guess(ions, electrons_nc);
auto result_nc = ground_state::calculate(ions, electrons_nc, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.1));
auto mag_nc = observables::total_magnetization(electrons_nc.spin_density());

match.check("total energy", result_col.energy.total(), result_nc.energy.total());
match.check("total magnetization", mag.length(), mag_nc.length());
}


int main (int argc, char ** argv) {
auto env = input::environment{};

auto ions = systems::ions(systems::cell::cubic(10.0_b));
ions.insert("H", {0.0_b, 0.0_b, 0.0_b});
compute_GS(env, ions);

auto d = 1.21_A;
ions = systems::ions(systems::cell::cubic(10.0_b));
ions.insert("O", {0.0_b, 0.0_b, d/2});
ions.insert("O", {0.0_b, 0.0_b,-d/2});
compute_GS(env, ions);
}
Loading

0 comments on commit 3eee860

Please sign in to comment.