diff --git a/src/ground_state/initial_guess.hpp b/src/ground_state/initial_guess.hpp index 64ee9df6..fbcf4bca 100644 --- a/src/ground_state/initial_guess.hpp +++ b/src/ground_state/initial_guess.hpp @@ -20,7 +20,7 @@ namespace inq { namespace ground_state { -void initial_guess(const systems::ions & ions, systems::electrons & electrons){ +void initial_guess(const systems::ions & ions, systems::electrons & electrons, std::optional> const & magnet_dir = {}){ int iphi = 0; for(auto & phi : electrons.kpin()) { @@ -42,6 +42,11 @@ void initial_guess(const systems::ions & ions, systems::electrons & electrons){ assert(fabs(operations::integral_sum(electrons.spin_density())) > 1e-16); observables::density::normalize(electrons.spin_density(), electrons.states().num_electrons()); + if (magnet_dir) { + assert(electrons.spin_density().set_size() > 1); + auto magnet_dir_ = {magnet_dir.value()[0], magnet_dir.value()[1], magnet_dir.value()[2]}; + observables::density::rotate_total_magnetization(electrons.spin_density(), magnet_dir_); + } } } @@ -55,8 +60,124 @@ void initial_guess(const systems::ions & ions, systems::electrons & electrons){ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { using namespace inq; + using namespace inq::magnitude; using namespace Catch::literals; using Catch::Approx; + + parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; + + SECTION("Spin unpolarized initialization") { + auto par = input::parallelization(comm); + auto ions = systems::ions(systems::cell::cubic(10.0_b)); + ions.insert("H", {0.0_b, 0.0_b, 0.0_b}); + auto electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_unpolarized()); + ground_state::initial_guess(ions, electrons); + CHECK(Approx(operations::integral_sum(electrons.spin_density())).epsilon(1.e-10) == 1.0); + } + + SECTION("Spin polarized initialization") { + auto par = input::parallelization(comm); + auto ions = systems::ions(systems::cell::cubic(10.0_b)); + ions.insert("H", {0.0_b, 0.0_b, 0.0_b}); + auto electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); + ground_state::initial_guess(ions, electrons); + CHECK(Approx(operations::integral_sum(electrons.spin_density())).epsilon(1.e-10) == 1.0); + + vector3 mag_dir = {0.0, 0.0, 1.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); + ground_state::initial_guess(ions, electrons, mag_dir); + auto mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral_sum(electrons.spin_density())).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0); + CHECK(Approx(sqrt(mag[0]*mag[0]+mag[1]*mag[1])/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + + mag_dir = {0.0, 0.0, -1.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); + ground_state::initial_guess(ions, electrons, mag_dir); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral_sum(electrons.spin_density())).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == -1.0); + CHECK(Approx(sqrt(mag[0]*mag[0]+mag[1]*mag[1])/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + } + + SECTION("Spin non collinear initialization") { + auto par = input::parallelization(comm); + auto ions = systems::ions(systems::cell::cubic(10.0_b)); + ions.insert("H", {0.0_b, 0.0_b, 0.0_b}); + auto electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons); + CHECK(Approx(operations::integral_sum(electrons.spin_density())).epsilon(1.e-10) == 1.0); + + vector3 mag_dir = {0.0, 0.0, 1.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + auto ch_density = observables::density::total(electrons.spin_density()); + auto mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0); + + mag_dir = {0.0, 0.0, -1.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + ch_density = observables::density::total(electrons.spin_density()); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == -1.0); + + mag_dir = {1.0, 1.0, 0.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + ch_density = observables::density::total(electrons.spin_density()); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + + mag_dir = {-1.0, 1.0, 0.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + ch_density = observables::density::total(electrons.spin_density()); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == -1.0/sqrt(2.0)); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + + mag_dir = {1.0, -1.0, 0.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + ch_density = observables::density::total(electrons.spin_density()); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == -1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + + mag_dir = {-1.0, -1.0, 0.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + ch_density = observables::density::total(electrons.spin_density()); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == -1.0/sqrt(2.0)); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == -1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 0.0); + + mag_dir = {1.0, 1.0, 1.0}; + electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons, mag_dir); + ch_density = observables::density::total(electrons.spin_density()); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(operations::integral(ch_density)).epsilon(1.e-10) == 1.0); + CHECK(Approx(mag[0]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(3.0)); + CHECK(Approx(mag[1]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(3.0)); + CHECK(Approx(mag[2]/sqrt(norm(mag))).epsilon(1.e-10) == 1.0/sqrt(3.0)); + } } #endif diff --git a/src/hamiltonian/energy.hpp b/src/hamiltonian/energy.hpp index 643aaf98..80f5944b 100644 --- a/src/hamiltonian/energy.hpp +++ b/src/hamiltonian/energy.hpp @@ -27,6 +27,7 @@ namespace hamiltonian { double xc_ = 0.0; double nvxc_ = 0.0; double exact_exchange_ = 0.0; + double zeeman_ener_ = 0.0; #ifdef ENABLE_CUDA public: @@ -149,6 +150,14 @@ namespace hamiltonian { nvxc_ = val; } + auto & zeeman_energy() const { + return zeeman_ener_; + } + + void zeeman_energy(double const & val) { + zeeman_ener_ = val; + } + auto & exact_exchange() const { return exact_exchange_; } @@ -186,6 +195,7 @@ namespace hamiltonian { utils::save_value(comm, dirname + "/xc", xc_, error_message); utils::save_value(comm, dirname + "/nvxc", nvxc_, error_message); utils::save_value(comm, dirname + "/exact_exchange", exact_exchange_, error_message); + utils::save_value(comm, dirname + "/zeeman_energy", zeeman_ener_, error_message); } static auto load(std::string const & dirname) { @@ -201,6 +211,7 @@ namespace hamiltonian { utils::load_value(dirname + "/xc", en.xc_, error_message); utils::load_value(dirname + "/nvxc", en.nvxc_, error_message); utils::load_value(dirname + "/exact_exchange", en.exact_exchange_, error_message); + utils::load_value(dirname + "/zeeman_energy", en.zeeman_ener_, error_message); return en; } @@ -219,6 +230,7 @@ namespace hamiltonian { tfm::format(out, " nvxc = %20.12f Ha\n", self.nvxc()); tfm::format(out, " exact-exchange = %20.12f Ha\n", self.exact_exchange()); tfm::format(out, " ion = %20.12f Ha\n", self.ion()); + tfm::format(out, " zeeman-energy = %20.12f Ha\n", self.zeeman_energy()); tfm::format(out, "\n"); return out; diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 6ad80382..c7dbc824 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -155,8 +156,17 @@ class self_consistency { hamiltonian.uniform_vector_potential_ += induced_vector_potential_; } - //the magnetic perturbation is not implemented yet - assert(not pert_.has_magnetic_field()); + // THE MAGNETIC FIELD + + if (pert_.has_magnetic_field()) { + basis::field> uniform_magnetic(spin_density.basis()); + uniform_magnetic.fill(vector3 {0.0, 0.0, 0.0}); + pert_.magnetic_field(time, uniform_magnetic); + zeeman_coupling zc_(spin_density.set_size()); + auto zeeman_ener = 0.0; + zc_(spin_density, uniform_magnetic, hamiltonian.scalar_potential_, zeeman_ener); + energy.zeeman_energy(zeeman_ener); + } } diff --git a/src/hamiltonian/xc_term.hpp b/src/hamiltonian/xc_term.hpp index 3a8e9074..d2176b61 100644 --- a/src/hamiltonian/xc_term.hpp +++ b/src/hamiltonian/xc_term.hpp @@ -102,8 +102,19 @@ class xc_term { template double compute_nvxc(SpinDensityType const & spin_density, VXC const & vxc) const { - + auto nvxc_ = 0.0; + if (spin_density.set_size() == 4) { + gpu::run(spin_density.local_set_size(), spin_density.basis().local_size(), + [vx = begin(vxc.matrix())] GPU_LAMBDA (auto is, auto ip){ + if (is == 2){ + vx[ip][is] = 2.0*vx[ip][is]; + } + else if (is == 3){ + vx[ip][is] = -2.0*vx[ip][is]; + } + }); + } nvxc_ += operations::integral_product_sum(spin_density, vxc); return nvxc_; } @@ -149,20 +160,20 @@ class xc_term { 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 b0 = 0.5*(vxi[ip][0] - vxi[ip][1]); + auto v0 = 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]; + vxf[ip][0] += v0 + b0*e_mag[2]; + vxf[ip][1] += v0 - b0*e_mag[2]; + vxf[ip][2] += b0*e_mag[0]; + vxf[ip][3] += b0*e_mag[1]; } else { - vxf[ip][0] += v; - vxf[ip][1] += v; + vxf[ip][0] += v0; + vxf[ip][1] += v0; } }); } @@ -177,79 +188,6 @@ class xc_term { //////////////////////////////////////////////////////////////////////////////////////////// - template - void eval_psi_vxc_psi(CommType & comm, CoreDensityType const & core_density, SpinDensityType const & spin_density, occupations_array_type const & occupations, kpin_type const & kpin, double & nvx, double & ex) { - - if (not any_true_functional()) { - nvx += 0.0; - ex += 0.0; - } - else { - auto full_density = process_density(spin_density, core_density); - double efunc = 0.0; - basis::field_set vxc(spin_density.skeleton()); - vxc.fill(0.0); - - basis::field_set vfunc(full_density.skeleton()); - auto density_gradient = std::optional{}; - 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); - compute_vxc(spin_density, vfunc, vxc); - ex += efunc; - } - - basis::field rfield(vxc.basis()); - rfield.fill(0.0); - int iphi = 0; - for (auto & phi : kpin) { - compute_psi_vxc_psi_ofr(occupations[iphi], phi, vxc, rfield); - iphi++; - } - - rfield.all_reduce(comm); - nvx += operations::integral(rfield); - } - } - - //////////////////////////////////////////////////////////////////////////////////////////// - - template - void compute_psi_vxc_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VxcType const & vxc, basis::field & rfield) { - - assert(get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); - assert(get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); - - if (vxc.set_size() == 1) { - gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), - [ph = begin(phi.matrix()), rf = begin(rfield.linear()), vx=begin(vxc.matrix()), occ=begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { - rf[ip] += occ[ist] * vx[ip][0] * norm(ph[ip][ist]); - }); - } - else if (vxc.set_size() == 2) { - gpu::run(phi.local_set_size(), phi.basis().local_size(), - [ph = begin(phi.matrix()), rf = begin(rfield.linear()), vx = begin(vxc.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { - rf[ip] += occ[ist] * vx[ip][spi] * norm(ph[ip][ist]); - }); - } - else { - assert(vxc.set_size() == 4); - gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), - [ph = begin(phi.spinor_array()), rf = begin(rfield.linear()), vx = begin(vxc.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { - auto offdiag = vx[ip][2] + complex{0.0, 1.0}*vx[ip][3]; - auto cross = 2.0*occ[ist]*real(offdiag*ph[ip][1][ist]*conj(ph[ip][0][ist])); - rf[ip] += occ[ist]*vx[ip][0]*norm(ph[ip][0][ist]); - rf[ip] += occ[ist]*vx[ip][1]*norm(ph[ip][1][ist]); - rf[ip] += cross; - }); - } - } - - //////////////////////////////////////////////////////////////////////////////////////////// - template static void evaluate_functional(hamiltonian::xc_functional const & functional, DensityType const & density, DensityGradientType const & density_gradient, double & efunctional, basis::field_set & vfunctional){ @@ -314,6 +252,12 @@ class xc_term { //////////////////////////////////////////////////////////////////////////////////////////// + auto & functionals() const { + return functionals_; + } + + //////////////////////////////////////////////////////////////////////////////////////////// + }; } } @@ -324,6 +268,80 @@ class xc_term { #include #include +using namespace inq; + +template +void compute_psi_vxc_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VxcType const & vxc, basis::field & rfield) { + + assert(get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); + assert(get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); + + if (vxc.set_size() == 1) { + gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), + [ph = begin(phi.matrix()), rf = begin(rfield.linear()), vx=begin(vxc.matrix()), occ=begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { + rf[ip] += occ[ist] * vx[ip][0] * norm(ph[ip][ist]); + }); + } + else if (vxc.set_size() == 2) { + gpu::run(phi.local_set_size(), phi.basis().local_size(), + [ph = begin(phi.matrix()), rf = begin(rfield.linear()), vx = begin(vxc.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { + rf[ip] += occ[ist] * vx[ip][spi] * norm(ph[ip][ist]); + }); + } + else { + assert(vxc.set_size() == 4); + gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), + [ph = begin(phi.spinor_array()), rf = begin(rfield.linear()), vx = begin(vxc.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { + auto offdiag = vx[ip][2] + complex{0.0, 1.0}*vx[ip][3]; + auto cross = 2.0*occ[ist]*real(offdiag*conj(ph[ip][1][ist])*ph[ip][0][ist]); + rf[ip] += occ[ist]*vx[ip][0]*norm(ph[ip][0][ist]); + rf[ip] += occ[ist]*vx[ip][1]*norm(ph[ip][1][ist]); + rf[ip] += cross; + }); + } +} + +//////////////////////////////////////////////////////////////////////////////////////////// + +template +void eval_psi_vxc_psi(CommType & comm, options::theory interaction, CoreDensityType const & core_density, SpinDensityType const & spin_density, occupations_array_type const & occupations, kpin_type const & kpin, double & nvx, double & ex) { + + hamiltonian::xc_term xc_(interaction, spin_density.set_size()); + + if (not xc_.any_true_functional()) { + nvx += 0.0; + ex += 0.0; + } + else { + auto full_density = xc_.process_density(spin_density, core_density); + double efunc = 0.0; + basis::field_set vxc(spin_density.skeleton()); + vxc.fill(0.0); + + basis::field_set vfunc(full_density.skeleton()); + auto density_gradient = std::optional{}; + if (xc_.any_requires_gradient()) density_gradient.emplace(operations::gradient(full_density)); + + for (auto & func : xc_.functionals()){ + if (not func.true_functional()) continue; + + xc_.evaluate_functional(func, full_density, density_gradient, efunc, vfunc); + xc_.compute_vxc(spin_density, vfunc, vxc); + ex += efunc; + } + + basis::field rfield(vxc.basis()); + rfield.fill(0.0); + int iphi = 0; + for (auto & phi : kpin) { + compute_psi_vxc_psi_ofr(occupations[iphi], phi, vxc, rfield); + iphi++; + } + + rfield.all_reduce(comm); + nvx += operations::integral(rfield); + } +} TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){ @@ -502,11 +520,10 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){ Approx target = Approx(nvxc).epsilon(1.e-10); Approx target2= Approx(exc).epsilon(1.e-10); - hamiltonian::xc_term xc_(options::theory{}.lda(), electrons.spin_density().set_size()); auto core_density_ = electrons.atomic_pot().nlcc_density(electrons.states_comm(), electrons.spin_density().basis(), ions); auto nvxc2 = 0.0; auto exc2 = 0.0; - xc_.eval_psi_vxc_psi(electrons.kpin_states_comm(), core_density_, electrons.spin_density(), electrons.occupations(), electrons.kpin(), nvxc2, exc2); + eval_psi_vxc_psi(electrons.kpin_states_comm(), options::theory{}.lda(), core_density_, electrons.spin_density(), electrons.occupations(), electrons.kpin(), nvxc2, exc2); CHECK(nvxc2 == target); CHECK(exc2 == target2); } @@ -523,11 +540,10 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){ Approx target = Approx(nvxc).epsilon(1.e-10); Approx target2 = Approx(exc).epsilon(1.e-10); - hamiltonian::xc_term xc_(options::theory{}.lda(), electrons.spin_density().set_size()); auto core_density_ = electrons.atomic_pot().nlcc_density(electrons.states_comm(), electrons.spin_density().basis(), ions); auto nvxc2 = 0.0; auto exc2 = 0.0; - xc_.eval_psi_vxc_psi(electrons.kpin_states_comm(), core_density_, electrons.spin_density(), electrons.occupations(), electrons.kpin(), nvxc2, exc2); + eval_psi_vxc_psi(electrons.kpin_states_comm(), options::theory{}.lda(), core_density_, electrons.spin_density(), electrons.occupations(), electrons.kpin(), nvxc2, exc2); CHECK(nvxc2 == target); CHECK(exc2 == target2); } @@ -544,11 +560,10 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){ Approx target = Approx(nvxc).epsilon(1.e-10); Approx target2 = Approx(exc).epsilon(1.e-10); - hamiltonian::xc_term xc_(options::theory{}.lda(), electrons.spin_density().set_size()); auto core_density_ = electrons.atomic_pot().nlcc_density(electrons.states_comm(), electrons.spin_density().basis(), ions); auto nvxc2 = 0.0; auto exc2 = 0.0; - xc_.eval_psi_vxc_psi(electrons.kpin_states_comm(), core_density_, electrons.spin_density(), electrons.occupations(), electrons.kpin(), nvxc2, exc2); + eval_psi_vxc_psi(electrons.kpin_states_comm(), options::theory{}.lda(), core_density_, electrons.spin_density(), electrons.occupations(), electrons.kpin(), nvxc2, exc2); CHECK(nvxc2 == target); CHECK(exc2 == target2); } diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp new file mode 100644 index 00000000..aaf5fad6 --- /dev/null +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -0,0 +1,301 @@ +/* -*- indent-tabs-mode: t -*- */ + +#ifndef ZEEMAN_COUPLING_HPP +#define ZEEMAN_COUPLING_HPP + +#include + +namespace inq { +namespace hamiltonian { + +class zeeman_coupling { + +private: + + int spin_components_; + +public: + + zeeman_coupling(int const spin_components): + spin_components_(spin_components) + { + assert(spin_components_ > 1); + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void operator()(SpinDensityType const & spin_density, basis::field> const & magnetic_field, VKSType & vks, double & zeeman_ener) const { + + basis::field_set zeeman_pot(vks.skeleton()); + zeeman_pot.fill(0.0); + + assert(zeeman_pot.set_size() == spin_components_); + + compute_zeeman_potential(magnetic_field, zeeman_pot); + + gpu::run(zeeman_pot.local_set_size(), zeeman_pot.basis().local_size(), + [vz = begin(zeeman_pot.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip) { + vk[ip][is] += vz[ip][is]; + }); + + zeeman_ener += compute_zeeman_energy(spin_density, zeeman_pot); + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void compute_zeeman_potential(basis::field> const & magnetic_field, VZType & zeeman_pot) const { + + gpu::run(zeeman_pot.basis().local_size(), + [vz = begin(zeeman_pot.matrix()), magnetic_ = begin(magnetic_field.linear())] GPU_LAMBDA (auto ip) { + vz[ip][0] +=-magnetic_[ip][2]; + vz[ip][1] += magnetic_[ip][2]; + }); + if (zeeman_pot.set_size() == 4) { + gpu::run(zeeman_pot.basis().local_size(), + [vz = begin(zeeman_pot.matrix()), magnetic_ = begin(magnetic_field.linear())] GPU_LAMBDA (auto ip) { + vz[ip][2] +=-magnetic_[ip][0]; + vz[ip][3] +=-magnetic_[ip][1]; + }); + } + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + double compute_zeeman_energy(SpinDensityType const & spin_density, VZType & zeeman_pot) const { + + auto zeeman_ener_ = 0.0; + if (spin_density.set_size() == 4) { + gpu::run(spin_density.local_set_size(), spin_density.basis().local_size(), + [vz = begin(zeeman_pot.matrix())] GPU_LAMBDA (auto is, auto ip) { + if (is == 2) { + vz[ip][is] = 2.0*vz[ip][is]; + } + else if (is == 3) { + vz[ip][is] = -2.0*vz[ip][is]; + } + }); + } + zeeman_ener_ += operations::integral_product_sum(spin_density, zeeman_pot); + return zeeman_ener_; + } +}; + +} +} +#endif + +#ifdef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST +#undef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST + +#include +#include +using namespace inq; + +template +void compute_psi_vz_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VZType const & zeeman_pot, RFType & rfield) { + + assert(get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); + assert(get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); + + if (zeeman_pot.set_size() == 2){ + gpu::run(phi.local_set_size(), phi.basis().local_size(), + [ph = begin(phi.matrix()), rf = begin(rfield.linear()), vz = begin(zeeman_pot.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { + rf[ip] += occ[ist]*vz[ip][spi]*norm(ph[ip][ist]); + }); + } + else { + assert(zeeman_pot.set_size() == 4); + gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), + [ph = begin(phi.spinor_array()), rf = begin(rfield.linear()), vz = begin(zeeman_pot.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { + auto offdiag = vz[ip][2] + complex{0.0, 1.0}*vz[ip][3]; + auto cross = 2.0*occ[ist]*real(offdiag*conj(ph[ip][1][ist])*ph[ip][0][ist]); + rf[ip] += occ[ist]*vz[ip][0]*norm(ph[ip][0][ist]); + rf[ip] += occ[ist]*vz[ip][1]*norm(ph[ip][1][ist]); + rf[ip] += cross; + }); + } +} + +//////////////////////////////////////////////////////////////////////////////////////////// + +template +void eval_psi_vz_psi(CommType & comm, SpinDensityType const & spin_density, MagneticField const & magnetic_field, occupations_array_type const & occupations, kpin_type const & kpin, double & zeeman_ener) { + + basis::field_set zeeman_pot(spin_density.skeleton()); + zeeman_pot.fill(0.0); + hamiltonian::zeeman_coupling zc_(spin_density.set_size()); + zc_.compute_zeeman_potential(magnetic_field, zeeman_pot); + + basis::field rfield(zeeman_pot.basis()); + rfield.fill(0.0); + int iphi = 0; + for (auto & phi : kpin) { + compute_psi_vz_psi_ofr(occupations[iphi], phi, zeeman_pot, rfield); + iphi++; + } + + rfield.all_reduce(comm); + zeeman_ener += operations::integral(rfield); +} + +TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { + + using namespace inq; + using namespace inq::magnitude; + using namespace Catch::literals; + using Catch::Approx; + + parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; + + SECTION("Spin polarized zeeman calculation") { + auto par = input::parallelization(comm); + auto ions = systems::ions(systems::cell::cubic(10.0_b)); + ions.insert("H", {0.0_b, 0.0_b, 0.0_b}); + auto electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); + ground_state::initial_guess(ions, electrons); + perturbations::magnetic magnetic_uniform{{0.0, 0.0, -1.0}}; + auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform); + auto mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 0.0); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) == 0.0); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) ==-1.0); + auto zeeman_ener = result.energy.zeeman_energy(); + Approx target = Approx(zeeman_ener).epsilon(1.e-10); + + basis::field> mag_field(electrons.spin_density().basis()); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform.magnetic_field(0.0, mag_field); + auto zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target); + } + + SECTION("Spin non collinear zeeman calculation") { + auto par = input::parallelization(comm); + auto ions = systems::ions(systems::cell::cubic(10.0_b)); + ions.insert("H", {0.0_b, 0.0_b, 0.0_b}); + auto electrons = systems::electrons(par, ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_non_collinear()); + ground_state::initial_guess(ions, electrons); + perturbations::magnetic magnetic_uniform{{0.0, 0.0, -1.0}}; + + auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform); + auto mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(sqrt(mag[0]*mag[0]+mag[1]*mag[1])/mag.length()).margin(1.e-7) == 0.0); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == -1.0); + + auto zeeman_ener = result.energy.zeeman_energy(); + Approx target = Approx(zeeman_ener).epsilon(1.e-10); + basis::field> mag_field(electrons.spin_density().basis()); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform.magnetic_field(0.0, mag_field); + auto zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target); + + vector3 bvec = {1.0, 1.0, 0.0}; + bvec = bvec / bvec.length(); + perturbations::magnetic magnetic_uniform2{bvec}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform2); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == 0.0); + + zeeman_ener = result.energy.zeeman_energy(); + Approx target2 = Approx(zeeman_ener).epsilon(1.e-10); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform2.magnetic_field(0.0, mag_field); + zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target2); + + bvec = {1.0, -1.0, 0.0}; + bvec = bvec / bvec.length(); + perturbations::magnetic magnetic_uniform3{bvec}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform3); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 1.0/sqrt(2.0)); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) ==-1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == 0.0); + + zeeman_ener = result.energy.zeeman_energy(); + Approx target3 = Approx(zeeman_ener).epsilon(1.e-10); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform3.magnetic_field(0.0, mag_field); + zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target3); + + bvec = {1.0, 1.0, 1.0}; + bvec = bvec / bvec.length(); + perturbations::magnetic magnetic_uniform4{bvec}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform4); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 1.0/sqrt(3.0)); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) == 1.0/sqrt(3.0)); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == 1.0/sqrt(3.0)); + + zeeman_ener = result.energy.zeeman_energy(); + Approx target4 = Approx(zeeman_ener).epsilon(1.e-10); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform4.magnetic_field(0.0, mag_field); + zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target4); + + bvec = {0.0, -1.0, 1.0}; + bvec = bvec / bvec.length(); + perturbations::magnetic magnetic_uniform5{bvec}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform5); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 0.0); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) ==-1.0/sqrt(2.0)); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == 1.0/sqrt(2.0)); + + zeeman_ener = result.energy.zeeman_energy(); + Approx target5 = Approx(zeeman_ener).epsilon(1.e-10); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform5.magnetic_field(0.0, mag_field); + zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target5); + + bvec = {1.0, -2.0, 1.5}; + bvec = bvec / bvec.length(); + perturbations::magnetic magnetic_uniform6{bvec}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform6); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 1.0/sqrt(1.0+4.0+9.0/4)); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) ==-2.0/sqrt(1.0+4.0+9.0/4)); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == 1.5/sqrt(1.0+4.0+9.0/4)); + + zeeman_ener = result.energy.zeeman_energy(); + Approx target6 = Approx(zeeman_ener).epsilon(1.e-10); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform6.magnetic_field(0.0, mag_field); + zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target6); + + bvec = {4.0, -2.0, 1.0}; + bvec = bvec / bvec.length(); + perturbations::magnetic magnetic_uniform7{bvec}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(200).mixing(0.1), magnetic_uniform7); + mag = observables::total_magnetization(electrons.spin_density()); + CHECK(Approx(mag[0]/mag.length()).margin(1.e-7) == 4.0/sqrt(16.0+4.0+1.0)); + CHECK(Approx(mag[1]/mag.length()).margin(1.e-7) ==-2.0/sqrt(16.0+4.0+1.0)); + CHECK(Approx(mag[2]/mag.length()).margin(1.e-7) == 1.0/sqrt(16.0+4.0+1.0)); + + zeeman_ener = result.energy.zeeman_energy(); + Approx target7 = Approx(zeeman_ener).epsilon(1.e-10); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + magnetic_uniform7.magnetic_field(0.0, mag_field); + zeeman_ener2 = 0.0; + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target7); + } +} +#endif diff --git a/src/observables/density.hpp b/src/observables/density.hpp index d2b819e2..bc1ed6ea 100644 --- a/src/observables/density.hpp +++ b/src/observables/density.hpp @@ -15,6 +15,7 @@ #include #include #include +#include namespace inq { namespace observables { @@ -102,7 +103,8 @@ void normalize(FieldType & density, const double & total_charge){ CALI_CXX_MARK_FUNCTION; - auto qq = operations::integral_sum(density); + auto max_index = std::min(2, density.set_size()); + auto qq = operations::integral_partial_sum(density, max_index); assert(fabs(qq) > 1e-16); gpu::run(density.local_set_size(), density.basis().local_size(), @@ -132,6 +134,39 @@ basis::field total(basis::field_set +void rotate_total_magnetization(FieldType & density, vector3 const & magnet_dir) { + + CALI_CXX_MARK_FUNCTION + + vector3 e_v = magnet_dir/sqrt(norm(magnet_dir)); + + if (density.set_size() == 2){ + gpu::run(density.basis().local_size(), + [den = begin(density.matrix()), mv = e_v[2]] GPU_LAMBDA (auto ip){ + auto n0 = den[ip][0] + den[ip][1]; + auto m0 = den[ip][0] - den[ip][1]; + den[ip][0] = 0.5*(n0 + m0*mv); + den[ip][1] = 0.5*(n0 - m0*mv); + }); + } + else { + assert(density.set_size() == 4); + gpu::run(density.basis().local_size(), + [den = begin(density.matrix()), mv = e_v] GPU_LAMBDA (auto ip){ + auto mag = observables::local_magnetization(den[ip], 4); + auto m0 = sqrt(norm(mag)); + auto n0 = den[ip][0] + den[ip][1]; + den[ip][0] = 0.5*(n0 + m0*mv[2]); + den[ip][1] = 0.5*(n0 - m0*mv[2]); + den[ip][2] = m0*mv[0]/2.0; + den[ip][3] = -m0*mv[1]/2.0; + }); + } +} + } } } diff --git a/src/observables/magnetization.hpp b/src/observables/magnetization.hpp index 8fa93df4..841cf30a 100644 --- a/src/observables/magnetization.hpp +++ b/src/observables/magnetization.hpp @@ -21,7 +21,7 @@ GPU_FUNCTION auto local_magnetization(Density const & spin_density, int const & if(components == 4){ mag_density[0] = 2.0*spin_density[2]; - mag_density[1] = 2.0*spin_density[3]; + mag_density[1] =-2.0*spin_density[3]; } else { mag_density[0] = 0.0; mag_density[1] = 0.0; diff --git a/src/operations/integral.hpp b/src/operations/integral.hpp index fca58560..0259e851 100644 --- a/src/operations/integral.hpp +++ b/src/operations/integral.hpp @@ -38,6 +38,24 @@ auto integral_sum(basis::field_set const & phi){ return integral_value; } +template +ElementType integral_partial_sum(basis::field_set const & phi, int const & max_index){ + CALI_CXX_MARK_FUNCTION; + + assert(phi.local_set_size() >= max_index); + basis::field rphi(phi.basis()); + rphi.fill(0.0); + gpu::run(phi.basis().local_size(), + [ph = begin(phi.matrix()), rph = begin(rphi.linear()), mi = max_index] GPU_LAMBDA (auto ip){ + for (auto i=0; i 1) phi.full_comm().all_reduce_in_place_n(&integral_value, 1, std::plus<>{}); + return integral_value; +} + template double integral_abs(basis::field const & phi){ CALI_CXX_MARK_FUNCTION; diff --git a/src/perturbations/magnetic.hpp b/src/perturbations/magnetic.hpp new file mode 100644 index 00000000..3159a6e4 --- /dev/null +++ b/src/perturbations/magnetic.hpp @@ -0,0 +1,66 @@ +/* -*- indent-tabs-mode: t -*- */ + +#ifndef INQ__PERTURBATIONS__MAGNET +#define INQ__PERTURBATIONS__MAGNET + +#include + +namespace inq { +namespace perturbations { + +class magnetic : public none { + + vector3 magnetic_vector_; + +public: + + magnetic(vector3 const & value): + magnetic_vector_(value) + { + } + + auto has_magnetic_field() const { + return true; + } + + template + void magnetic_field(const double time, MagneticField & magnetic) const { + gpu::run(magnetic.basis().local_size(), + [magnetic_ = begin(magnetic.linear()), mv = magnetic_vector_] GPU_LAMBDA (auto ip){ + magnetic_[ip] += mv; + }); + } + + template + friend OStream & operator<<(OStream & out, magnetic const & self){ + return out; + } + +}; + +} +} +#endif + +#ifdef INQ_PERTURBATIONS_MAGNETIC_UNIT_TEST +#undef INQ_PERTURBATIONS_MAGNETIC_UNIT_TEST + +TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { + + parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; + + perturbations::magnetic uniform_magnetic{{0.0, 0.0, 1.0}}; + + basis::real_space bas(systems::cell::cubic(5.0_b), /*spacing*/ 0.1, comm); + basis::field> mag_field(bas); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); + + CHECK(uniform_magnetic.has_magnetic_field()); + uniform_magnetic.magnetic_field(/*time*/ 0.0, mag_field); + CHECK(mag_field.linear()[0] == vector3{0.0, 0.0, 1.0}); + CHECK(mag_field.linear()[1] == vector3{0.0, 0.0, 1.0}); + + uniform_magnetic.magnetic_field(/*time*/ 1000.0, mag_field); + CHECK(mag_field.linear()[0] == vector3{0.0, 0.0, 2.0}); +} +#endif \ No newline at end of file