From f249fd6c2d0fa3cc0de229c0f6ee853868ac6eac Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 14 Oct 2024 09:16:38 -0500 Subject: [PATCH 01/40] inclusion magnetic field section: if has_magnetic_field call zeeman_coupling --- src/hamiltonian/self_consistency.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 6ad80382..52ffe7f5 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,12 @@ 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()) { + std::cout << "MAGNETIC FIELD ACTIVE" << std::endl; + zeeman_coupling zc_(spin_density.set_size(), pert_.uniform_magnetic_field(time)); + } } From 1836257333b214b9071ff91b8a50ce6d3040c299 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 14 Oct 2024 09:18:05 -0500 Subject: [PATCH 02/40] correction in the non collinear calculation of nvxc: the factor 2 is needed for the product with the spin density --- src/hamiltonian/xc_term.hpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/hamiltonian/xc_term.hpp b/src/hamiltonian/xc_term.hpp index 3a8e9074..90d56c2f 100644 --- a/src/hamiltonian/xc_term.hpp +++ b/src/hamiltonian/xc_term.hpp @@ -102,8 +102,16 @@ 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(), + [v = begin(vxc.matrix())] GPU_LAMBDA (auto is, auto ip){ + if (is >= 2){ + v[ip][is] = 2.0*v[ip][is]; + } + }); + } nvxc_ += operations::integral_product_sum(spin_density, vxc); return nvxc_; } From b2bd4c4459b4686f7870de0f7e81a37ee282c837 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 14 Oct 2024 09:18:58 -0500 Subject: [PATCH 03/40] zeeman coupling class initialization - only constructor and basic test --- src/hamiltonian/zeeman_coupling.hpp | 44 +++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 src/hamiltonian/zeeman_coupling.hpp diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp new file mode 100644 index 00000000..e1155481 --- /dev/null +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -0,0 +1,44 @@ +/* -*- indent-tabs-mode: t -*- */ + +#ifndef ZEEMAN_COUPL_HPP +#define ZEEMAN_COUPL_HPP + +namespace inq { +namespace hamiltonian { + +class zeeman_coupling { + +private: + + vector3 B_; + int spin_components_; + +public: + + zeeman_coupling(int const spin_components, vector3 const & B): + spin_components_(spin_components), + B_(B) + { + assert(spin_components_ > 1); + std::cout << "SPIN COMPONENTS : " << spin_components_ << std::endl; + std::cout << B_ << std::endl; + } + +}; + +} +} +#endif + +#ifdef INQ_HAMILTONIAN_ZEEMAN_COUPL_UNIT_TEST +#undef INQ_HAMILTONIAN_ZEEMAN_COUPL_UNIT_TEST + +TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { + + parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; + + SECTION("Spin polarized vz calculation") { + auto par = input::parallelization(comm); + } +} +#endif \ No newline at end of file From 25ff28235ee13bde9b07c91a95fe9ca86d2b6e42 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 14 Oct 2024 09:20:06 -0500 Subject: [PATCH 04/40] simple homogeneous magnetic field class analogous to simple_electric_field --- src/perturbations/magnetic_field.hpp | 52 ++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 src/perturbations/magnetic_field.hpp diff --git a/src/perturbations/magnetic_field.hpp b/src/perturbations/magnetic_field.hpp new file mode 100644 index 00000000..dd63e29d --- /dev/null +++ b/src/perturbations/magnetic_field.hpp @@ -0,0 +1,52 @@ +/* -*- indent-tabs-mode: t -*- */ + +#ifndef INQ__PERTURBATIONS__MAGNETIC_FIELD +#define INQ__PERTURBATIONS__MAGNETIC_FIELD + +#include + +namespace inq { +namespace perturbations { + +class magnetic_field : public none { + + vector3 uniform_magnetic_field_; + +public: + + magnetic_field(vector3 const & field): + uniform_magnetic_field_(field) + { + } + + auto has_uniform_magnetic_field() const { + return true; + } + + auto uniform_magnetic_field(double /*time*/) const { + return uniform_magnetic_field_; + } + + template + friend OStream & operator<<(OStream & out, magnetic_field const & self){ + return out; + } + +}; + +} +} +#endif + +#ifdef INQ__PERTURBATIONS__MAGNETIC_FIELD_UNIT_TEST +#undef INQ__PERTURBATIONS__MAGNETIC_FIELD_UNIT_TEST + +TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { + + perturbations::magnetic_field bfield{{0.0, 0.0, 1.0}}; + + CHECK(bfield.has_magnetic_field()); + CHECK(bfield.uniform_magnetic_field(/*time*/ 0.0) == vector3{0.0, 0.0, 1.0}); + CHECK(bfield.uniform_magnetic_field(/*time*/ 1000.0) == vector3{0.0, 0.0, 1.0}); +} +#endif \ No newline at end of file From 08885b3e037accb802d3058eb2699c0e5952ffa9 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 14 Oct 2024 09:20:45 -0500 Subject: [PATCH 05/40] inclusion of a simple magnetic field test for hydrogen atom --- tests/ext_mag_field.cpp | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 tests/ext_mag_field.cpp diff --git a/tests/ext_mag_field.cpp b/tests/ext_mag_field.cpp new file mode 100644 index 00000000..e1be08d2 --- /dev/null +++ b/tests/ext_mag_field.cpp @@ -0,0 +1,19 @@ +/* -*- indent-tabs-mode: t -*- */ + +#include +#include + +using namespace inq; +using namespace inq::magnitude; + +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}); + + auto electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); + ground_state::initial_guess(ions, electrons); + perturbations::magnetic_field bfield{{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(1000).mixing(0.1), bfield); +} \ No newline at end of file From d35c1cf84662abcd4008cca9767440d03bcc0a3e Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 04:35:50 -0500 Subject: [PATCH 06/40] added external B field declaration passed to pert_.magnetic_field and then to zeeman coupling --- src/hamiltonian/self_consistency.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 52ffe7f5..a5ffafac 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -160,7 +160,10 @@ class self_consistency { if (pert_.has_magnetic_field()) { std::cout << "MAGNETIC FIELD ACTIVE" << std::endl; - zeeman_coupling zc_(spin_density.set_size(), pert_.uniform_magnetic_field(time)); + basis::field> B(spin_density.basis()); + B.fill(vector3 {0.0, 0.0, 0.0}); + pert_.magnetic_field(time, B); + zeeman_coupling zc_(spin_density.set_size()); } } From bb0fa4d96b629ed9cbd1c449fe48eb046ec612e6 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 04:38:11 -0500 Subject: [PATCH 07/40] modified zeeman_coupling - only private member is num_spin_components - the B field is not stored as a member --- src/hamiltonian/zeeman_coupling.hpp | 30 ++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index e1155481..3094a7ef 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -1,7 +1,7 @@ /* -*- indent-tabs-mode: t -*- */ -#ifndef ZEEMAN_COUPL_HPP -#define ZEEMAN_COUPL_HPP +#ifndef ZEEMAN_COUPLING_HPP +#define ZEEMAN_COUPLING_HPP namespace inq { namespace hamiltonian { @@ -10,18 +10,15 @@ class zeeman_coupling { private: - vector3 B_; int spin_components_; public: - zeeman_coupling(int const spin_components, vector3 const & B): - spin_components_(spin_components), - B_(B) + zeeman_coupling(int const spin_components): + spin_components_(spin_components) { assert(spin_components_ > 1); - std::cout << "SPIN COMPONENTS : " << spin_components_ << std::endl; - std::cout << B_ << std::endl; + std::cout << "SPIN COMPONENTS : " << spin_components_ << std::endl; } }; @@ -30,15 +27,26 @@ class zeeman_coupling { } #endif -#ifdef INQ_HAMILTONIAN_ZEEMAN_COUPL_UNIT_TEST -#undef INQ_HAMILTONIAN_ZEEMAN_COUPL_UNIT_TEST +#ifdef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST +#undef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST + +#include TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { + using namespace inq; + using namespace inq::magnitude; + parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; - SECTION("Spin polarized vz calculation") { + 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 B{{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(1000).mixing(0.1), B); } } #endif \ No newline at end of file From f5319d755c8badd6aa5146f4e45ebef2020b2af7 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 04:39:17 -0500 Subject: [PATCH 08/40] removed older version of magnetic field implementation --- src/perturbations/magnetic_field.hpp | 52 ---------------------------- 1 file changed, 52 deletions(-) delete mode 100644 src/perturbations/magnetic_field.hpp diff --git a/src/perturbations/magnetic_field.hpp b/src/perturbations/magnetic_field.hpp deleted file mode 100644 index dd63e29d..00000000 --- a/src/perturbations/magnetic_field.hpp +++ /dev/null @@ -1,52 +0,0 @@ -/* -*- indent-tabs-mode: t -*- */ - -#ifndef INQ__PERTURBATIONS__MAGNETIC_FIELD -#define INQ__PERTURBATIONS__MAGNETIC_FIELD - -#include - -namespace inq { -namespace perturbations { - -class magnetic_field : public none { - - vector3 uniform_magnetic_field_; - -public: - - magnetic_field(vector3 const & field): - uniform_magnetic_field_(field) - { - } - - auto has_uniform_magnetic_field() const { - return true; - } - - auto uniform_magnetic_field(double /*time*/) const { - return uniform_magnetic_field_; - } - - template - friend OStream & operator<<(OStream & out, magnetic_field const & self){ - return out; - } - -}; - -} -} -#endif - -#ifdef INQ__PERTURBATIONS__MAGNETIC_FIELD_UNIT_TEST -#undef INQ__PERTURBATIONS__MAGNETIC_FIELD_UNIT_TEST - -TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { - - perturbations::magnetic_field bfield{{0.0, 0.0, 1.0}}; - - CHECK(bfield.has_magnetic_field()); - CHECK(bfield.uniform_magnetic_field(/*time*/ 0.0) == vector3{0.0, 0.0, 1.0}); - CHECK(bfield.uniform_magnetic_field(/*time*/ 1000.0) == vector3{0.0, 0.0, 1.0}); -} -#endif \ No newline at end of file From bc76a580804746e666b96f3fbb679f0c71afdc42 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 04:40:48 -0500 Subject: [PATCH 09/40] updated ext_mag_field test - run it on spin polarized O2 under external B field --- tests/ext_mag_field.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/tests/ext_mag_field.cpp b/tests/ext_mag_field.cpp index e1be08d2..5045d799 100644 --- a/tests/ext_mag_field.cpp +++ b/tests/ext_mag_field.cpp @@ -1,19 +1,27 @@ /* -*- indent-tabs-mode: t -*- */ #include -#include +#include using namespace inq; using namespace inq::magnitude; +inq::utils::match match(3.0e-4); + int main (int argc, char ** argv){ auto env = input::environment{}; + auto d = 1.21_A; auto ions = systems::ions(systems::cell::cubic(10.0_b)); - ions.insert("H", {0.0_b, 0.0_b, 0.0_b}); + ions.insert("O", {0.0_b, 0.0_b, d/2}); + ions.insert("O", {0.0_b, 0.0_b,-d/2}); auto electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); ground_state::initial_guess(ions, electrons); - perturbations::magnetic_field bfield{{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(1000).mixing(0.1), bfield); + perturbations::magnetic B{{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(1000).mixing(0.1), B); + auto mag = observables::total_magnetization(electrons.spin_density()); + + match.check("magnetization direction", mag/mag.length(), {0.0, 0.0, 1.0}); + return match.fail(); } \ No newline at end of file From 4624ebd1fb8f8cc083d54b9edd2d863054685183 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 04:42:51 -0500 Subject: [PATCH 10/40] new version of magnetic field implementation: only local magnetic 3vector is stored - MagneticField used to update the external magnetic field on real space grid --- src/perturbations/magnetic.hpp | 66 ++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100644 src/perturbations/magnetic.hpp diff --git a/src/perturbations/magnetic.hpp b/src/perturbations/magnetic.hpp new file mode 100644 index 00000000..a8ccee85 --- /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 & B): + magnetic_vector_(B) + { + } + + auto has_magnetic_field() const { + return true; + } + + template + void magnetic_field(const double time, MagneticField & magnetic) const { + gpu::run(magnetic.basis().local_size(), + [b = begin(magnetic.linear()), mv = magnetic_vector_] GPU_LAMBDA (auto ip){ + b[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 B{{0.0, 0.0, 1.0}}; + + basis::real_space bas(systems::cell::cubic(5.0_b), /*spacing*/ 0.1, comm); + basis::field> Bfield(bas); + Bfield.fill(vector3 {0.0, 0.0, 0.0}); + + CHECK(B.has_magnetic_field()); + B.magnetic_field(/*time*/ 0.0, Bfield); + CHECK(Bfield.linear()[0] == vector3{0.0, 0.0, 1.0}); + CHECK(Bfield.linear()[1] == vector3{0.0, 0.0, 1.0}); + + B.magnetic_field(/*time*/ 1000.0, Bfield); + CHECK(Bfield.linear()[0] == vector3{0.0, 0.0, 2.0}); +} +#endif \ No newline at end of file From 6a25a6cf315bc8b1393317c01a550cfda29c805e Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Thu, 17 Oct 2024 10:09:50 -0500 Subject: [PATCH 11/40] call zeeman coupling calculation --- src/hamiltonian/self_consistency.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index a5ffafac..6f7ddf95 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -164,6 +164,9 @@ class self_consistency { B.fill(vector3 {0.0, 0.0, 0.0}); pert_.magnetic_field(time, B); zeeman_coupling zc_(spin_density.set_size()); + auto nvz = 0.0; + zc_(spin_density, B, hamiltonian.scalar_potential_, nvz); + std::cout << " nvz -> " << nvz << std::endl; } } From b73d03ec2271285b463ac1e9b142f3b294622cab Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Thu, 17 Oct 2024 10:10:22 -0500 Subject: [PATCH 12/40] zeeman coupling () implementation --- src/hamiltonian/zeeman_coupling.hpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 3094a7ef..3ffb35e3 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -3,6 +3,8 @@ #ifndef ZEEMAN_COUPLING_HPP #define ZEEMAN_COUPLING_HPP +#include + namespace inq { namespace hamiltonian { @@ -18,7 +20,18 @@ class zeeman_coupling { spin_components_(spin_components) { assert(spin_components_ > 1); - std::cout << "SPIN COMPONENTS : " << spin_components_ << std::endl; + std::cout << "SPIN COMPONENTS : " << spin_components_ << std::endl; + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void operator()(SpinDensityType const & spin_density, MagneticField const & B, VKSType & vks, double & nvz) const { + + basis::field_set vz(vks.skeleton()); + vz.fill(0.0); + + assert(vz.set_size() == spin_components_); } }; From cd20b440a5cb25664693752e7e8ba8c19f26c8af Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 21:22:48 -0500 Subject: [PATCH 13/40] ext. magnetic field for O2 --- tests/ext_mag_field.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/ext_mag_field.cpp b/tests/ext_mag_field.cpp index 5045d799..c56f6523 100644 --- a/tests/ext_mag_field.cpp +++ b/tests/ext_mag_field.cpp @@ -16,11 +16,12 @@ int main (int argc, char ** argv){ ions.insert("O", {0.0_b, 0.0_b, d/2}); ions.insert("O", {0.0_b, 0.0_b,-d/2}); - auto electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(2).spin_polarized()); + auto electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(10).spin_polarized()); ground_state::initial_guess(ions, electrons); - perturbations::magnetic B{{0.0, 0.0, 1.0}}; + perturbations::magnetic B{{0.0, 0.0, 0.4}}; auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.1), B); auto mag = observables::total_magnetization(electrons.spin_density()); + std::cout << mag << std::endl; match.check("magnetization direction", mag/mag.length(), {0.0, 0.0, 1.0}); return match.fail(); From 82e4a18365af2866d3de9bbc8b2176b9ba2296c7 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Thu, 17 Oct 2024 21:23:48 -0500 Subject: [PATCH 14/40] added: compute_vz, compute_nvz, process_potential --- src/hamiltonian/zeeman_coupling.hpp | 49 +++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 3ffb35e3..10fb8b19 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -32,8 +32,57 @@ class zeeman_coupling { vz.fill(0.0); assert(vz.set_size() == spin_components_); + + compute_vz(B, vz); + + process_potential(vz, vks); + + nvz += compute_nvz(spin_density, vz); + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void compute_vz(MagneticField const & B, VZType & vz) const { + + if (vz.set_size() == 4) { + } + else { + assert(vz.set_size() == 2); + gpu::run(vz.basis().local_size(), + [v = begin(vz.matrix()), b = begin(B.linear())] GPU_LAMBDA (auto ip) { + v[ip][0] +=-b[ip][2]; + v[ip][1] += b[ip][2]; + }); + } } + //////////////////////////////////////////////////////////////////////////////////////////// + + template + double compute_nvz(SpinDensityType const & spin_density, VZType & vz) const { + + auto nvz_ = 0.0; + if (spin_density.set_size() == 4) { + gpu::run(spin_density.local_set_size(), spin_density.basis().local_size(), + [v = begin(vz.matrix())] GPU_LAMBDA (auto is, auto ip) { + if (is >= 2) v[ip][is] = 2.0*v[ip][is]; + }); + } + nvz_ += operations::integral_product_sum(spin_density, vz); + return nvz_; + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void process_potential(VZType const & vz, VKSType & vks) const { + + gpu::run(vz.local_set_size(), vz.basis().local_size(), + [v = begin(vz.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip) { + vk[ip][is] += v[ip][is]; + }); + } }; } From 5599894a7cd2b6c73cf2e9b5de699108ca15f445 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 21 Oct 2024 09:42:40 -0500 Subject: [PATCH 15/40] included zeeman energy internal variable and output method --- src/hamiltonian/energy.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/hamiltonian/energy.hpp b/src/hamiltonian/energy.hpp index 643aaf98..62e8d80f 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 nvz_ = 0.0; #ifdef ENABLE_CUDA public: @@ -149,6 +150,14 @@ namespace hamiltonian { nvxc_ = val; } + auto & nvz() const { + return nvz_; + } + + void nvz(double const & val) { + nvz_ = 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 + "/nvz", nvz_, 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 + "/nvz", en.nvz_, 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, " nvz = %20.12f Ha\n", self.nvz()); tfm::format(out, "\n"); return out; From 66e2e4849b82c4d37f17af04dec26f49db44237e Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 21 Oct 2024 09:43:24 -0500 Subject: [PATCH 16/40] save zeeman energy --- src/hamiltonian/self_consistency.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 6f7ddf95..92a1c01a 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -166,7 +166,8 @@ class self_consistency { zeeman_coupling zc_(spin_density.set_size()); auto nvz = 0.0; zc_(spin_density, B, hamiltonian.scalar_potential_, nvz); - std::cout << " nvz -> " << nvz << std::endl; + energy.nvz(nvz); + std::cout << " nvz -> " << energy.nvz() << std::endl; } } From 788fd9b93293e75252acd722be42948a6ec53a9b Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 21 Oct 2024 09:45:39 -0500 Subject: [PATCH 17/40] added to internal methods to compute the zeeman energy from the wave functions - the test units use these to compare the computed zeeman energy during SCF iteration and this post process calculation --- src/hamiltonian/zeeman_coupling.hpp | 99 +++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 10fb8b19..a4024244 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -46,6 +46,13 @@ class zeeman_coupling { void compute_vz(MagneticField const & B, VZType & vz) const { if (vz.set_size() == 4) { + gpu::run(vz.basis().local_size(), + [v = begin(vz.matrix()), b = begin(B.linear())] GPU_LAMBDA (auto ip) { + v[ip][0] +=-b[ip][2]; + v[ip][1] += b[ip][2]; + v[ip][2] +=-b[ip][0]; + v[ip][3] +=-b[ip][1]; + }); } else { assert(vz.set_size() == 2); @@ -83,6 +90,54 @@ class zeeman_coupling { vk[ip][is] += v[ip][is]; }); } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void compute_psi_vz_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VZType const & vz, RFType & rfield) { + + assert(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); + assert(std::get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); + + if (vz.set_size() == 2){ + gpu::run(phi.local_set_size(), phi.basis().local_size(), + [ph = begin(phi.matrix()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { + rf[ip] += occ[ist]*v[ip][spi]*norm(ph[ip][ist]); + }); + } + else { + assert(vz.set_size() == 4); + gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), + [ph = begin(phi.spinor_array()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { + auto offdiag = v[ip][2] + complex{0.0, 1.0}*v[ip][3]; + auto cross = 2.0*occ[ist]*real(offdiag*ph[ip][1][ist]*conj(ph[ip][0][ist])); + rf[ip] += occ[ist]*v[ip][0]*norm(ph[ip][0][ist]); + rf[ip] += occ[ist]*v[ip][1]*norm(ph[ip][1][ist]); + rf[ip] += cross; + }); + } + } + + //////////////////////////////////////////////////////////////////////////////////////////// + + template + void eval_psi_vz_psi(CommType & comm, SpinDensityType const & spin_density, MagneticField const & B, occupations_array_type const & occupations, kpin_type const & kpin, double & nvz) { + + basis::field_set vz(spin_density.skeleton()); + vz.fill(0.0); + compute_vz(B, vz); + + basis::field rfield(vz.basis()); + rfield.fill(0.0); + int iphi = 0; + for (auto & phi : kpin) { + compute_psi_vz_psi_ofr(occupations[iphi], phi, vz, rfield); + iphi++; + } + + rfield.all_reduce(comm); + nvz += operations::integral(rfield); + } }; } @@ -93,11 +148,13 @@ class zeeman_coupling { #undef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST #include +#include TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { using namespace inq; using namespace inq::magnitude; + using Catch::Approx; parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; @@ -109,6 +166,48 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { ground_state::initial_guess(ions, electrons); perturbations::magnetic B{{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(1000).mixing(0.1), B); + auto mag = observables::total_magnetization(electrons.spin_density()); + CHECK(mag[0]/mag.length() == 0.0); + CHECK(mag[1]/mag.length() == 0.0); + CHECK(mag[2]/mag.length() ==-1.0); + auto nvz = result.energy.nvz(); + Approx target = Approx(nvz).epsilon(1.e-10); + + hamiltonian::zeeman_coupling zc_(electrons.states().num_density_components()); + basis::field> Bfield(electrons.spin_density().basis()); + Bfield.fill(vector3 {0.0, 0.0, 0.0}); + B.magnetic_field(0.0, Bfield); + auto nvz2 = 0.0; + zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), nvz2); + CHECK(nvz2 == 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 B{{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(10000).mixing(0.01), B); + auto mag = observables::total_magnetization(electrons.spin_density()); + auto mx = mag[0]/mag.length(); + auto my = mag[1]/mag.length(); + auto mz = mag[2]/mag.length(); + CHECK(abs(mx) < 1.e-7); + CHECK(abs(my) < 1.e-7); + CHECK(abs(mz + 1.0) < 1.e-7); + + auto nvz = result.energy.nvz(); + Approx target = Approx(nvz).epsilon(1.e-10); + hamiltonian::zeeman_coupling zc_(electrons.states().num_density_components()); + basis::field> Bfield(electrons.spin_density().basis()); + Bfield.fill(vector3 {0.0, 0.0, 0.0}); + B.magnetic_field(0.0, Bfield); + auto nvz2 = 0.0; + zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), nvz2); + CHECK(nvz2 == target); } } #endif \ No newline at end of file From 057a8c4a4bc1a822008e8a4d53ad85c16cdf297e Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Mon, 21 Oct 2024 09:46:44 -0500 Subject: [PATCH 18/40] included the non collinear part of the test and the collinear with opposite B field --- tests/ext_mag_field.cpp | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/tests/ext_mag_field.cpp b/tests/ext_mag_field.cpp index c56f6523..235d63cc 100644 --- a/tests/ext_mag_field.cpp +++ b/tests/ext_mag_field.cpp @@ -22,7 +22,23 @@ int main (int argc, char ** argv){ auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.1), B); auto mag = observables::total_magnetization(electrons.spin_density()); std::cout << mag << std::endl; - match.check("magnetization direction", mag/mag.length(), {0.0, 0.0, 1.0}); + + perturbations::magnetic B2{{0.0, 0.0, -0.4}}; + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.1), B2); + auto mag2 = observables::total_magnetization(electrons.spin_density()); + std::cout << mag2 << std::endl; + match.check("magnetization direction", mag2/mag2.length(), {0.0, 0.0, -1.0}); + + perturbations::magnetic B3{{0.001, 0.0, -0.4}}; + electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(10).spin_non_collinear()); + ground_state::initial_guess(ions, electrons); + result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.01), B3); + auto mag3 = observables::total_magnetization(electrons.spin_density()); + std::cout << mag3 << std::endl; + vector3 e_v = {0.001, 0.0, -0.4}; + e_v = e_v/sqrt(norm(e_v)); + match.check("magnetization direction", mag3/mag3.length(), e_v); + return match.fail(); } \ No newline at end of file From 509011cb6f39bf7907ecfe9c5180652b206493cd Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 23 Oct 2024 07:52:31 -0500 Subject: [PATCH 19/40] changed nvz -> zeeman_ener_ --- src/hamiltonian/energy.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/hamiltonian/energy.hpp b/src/hamiltonian/energy.hpp index 62e8d80f..80f5944b 100644 --- a/src/hamiltonian/energy.hpp +++ b/src/hamiltonian/energy.hpp @@ -27,7 +27,7 @@ namespace hamiltonian { double xc_ = 0.0; double nvxc_ = 0.0; double exact_exchange_ = 0.0; - double nvz_ = 0.0; + double zeeman_ener_ = 0.0; #ifdef ENABLE_CUDA public: @@ -150,12 +150,12 @@ namespace hamiltonian { nvxc_ = val; } - auto & nvz() const { - return nvz_; + auto & zeeman_energy() const { + return zeeman_ener_; } - void nvz(double const & val) { - nvz_ = val; + void zeeman_energy(double const & val) { + zeeman_ener_ = val; } auto & exact_exchange() const { @@ -195,7 +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 + "/nvz", nvz_, error_message); + utils::save_value(comm, dirname + "/zeeman_energy", zeeman_ener_, error_message); } static auto load(std::string const & dirname) { @@ -211,7 +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 + "/nvz", en.nvz_, error_message); + utils::load_value(dirname + "/zeeman_energy", en.zeeman_ener_, error_message); return en; } @@ -230,7 +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, " nvz = %20.12f Ha\n", self.nvz()); + tfm::format(out, " zeeman-energy = %20.12f Ha\n", self.zeeman_energy()); tfm::format(out, "\n"); return out; From 8ea850279f9f79ced69a80c6d9958e05df8ff5c9 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 23 Oct 2024 07:55:12 -0500 Subject: [PATCH 20/40] changed nvz with zeeman_ener variable + removed cout statements --- src/hamiltonian/self_consistency.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 92a1c01a..5ef33cd2 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -159,15 +159,13 @@ class self_consistency { // THE MAGNETIC FIELD if (pert_.has_magnetic_field()) { - std::cout << "MAGNETIC FIELD ACTIVE" << std::endl; basis::field> B(spin_density.basis()); B.fill(vector3 {0.0, 0.0, 0.0}); pert_.magnetic_field(time, B); zeeman_coupling zc_(spin_density.set_size()); - auto nvz = 0.0; - zc_(spin_density, B, hamiltonian.scalar_potential_, nvz); - energy.nvz(nvz); - std::cout << " nvz -> " << energy.nvz() << std::endl; + auto zeeman_ener = 0.0; + zc_(spin_density, B, hamiltonian.scalar_potential_, zeeman_ener); + energy.zeeman_energy(zeeman_ener); } } From ce293a59f9cd994ed23b7dc04a53cb7b938bf5db Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 23 Oct 2024 07:56:44 -0500 Subject: [PATCH 21/40] replaced nvz with zeeman_ener variable + cout removed --- src/hamiltonian/zeeman_coupling.hpp | 41 ++++++++++++++--------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index a4024244..9433412f 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -20,13 +20,12 @@ class zeeman_coupling { spin_components_(spin_components) { assert(spin_components_ > 1); - std::cout << "SPIN COMPONENTS : " << spin_components_ << std::endl; } //////////////////////////////////////////////////////////////////////////////////////////// template - void operator()(SpinDensityType const & spin_density, MagneticField const & B, VKSType & vks, double & nvz) const { + void operator()(SpinDensityType const & spin_density, MagneticField const & B, VKSType & vks, double & zeeman_ener) const { basis::field_set vz(vks.skeleton()); vz.fill(0.0); @@ -37,7 +36,7 @@ class zeeman_coupling { process_potential(vz, vks); - nvz += compute_nvz(spin_density, vz); + zeeman_ener += compute_zeeman_energy(spin_density, vz); } //////////////////////////////////////////////////////////////////////////////////////////// @@ -67,17 +66,17 @@ class zeeman_coupling { //////////////////////////////////////////////////////////////////////////////////////////// template - double compute_nvz(SpinDensityType const & spin_density, VZType & vz) const { + double compute_zeeman_energy(SpinDensityType const & spin_density, VZType & vz) const { - auto nvz_ = 0.0; + auto zeeman_ener_ = 0.0; if (spin_density.set_size() == 4) { gpu::run(spin_density.local_set_size(), spin_density.basis().local_size(), [v = begin(vz.matrix())] GPU_LAMBDA (auto is, auto ip) { if (is >= 2) v[ip][is] = 2.0*v[ip][is]; }); } - nvz_ += operations::integral_product_sum(spin_density, vz); - return nvz_; + zeeman_ener_ += operations::integral_product_sum(spin_density, vz); + return zeeman_ener_; } //////////////////////////////////////////////////////////////////////////////////////////// @@ -121,7 +120,7 @@ class zeeman_coupling { //////////////////////////////////////////////////////////////////////////////////////////// template - void eval_psi_vz_psi(CommType & comm, SpinDensityType const & spin_density, MagneticField const & B, occupations_array_type const & occupations, kpin_type const & kpin, double & nvz) { + void eval_psi_vz_psi(CommType & comm, SpinDensityType const & spin_density, MagneticField const & B, occupations_array_type const & occupations, kpin_type const & kpin, double & zeeman_ener) { basis::field_set vz(spin_density.skeleton()); vz.fill(0.0); @@ -136,7 +135,7 @@ class zeeman_coupling { } rfield.all_reduce(comm); - nvz += operations::integral(rfield); + zeeman_ener += operations::integral(rfield); } }; @@ -165,21 +164,21 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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 B{{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(1000).mixing(0.1), B); + 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), B); auto mag = observables::total_magnetization(electrons.spin_density()); CHECK(mag[0]/mag.length() == 0.0); CHECK(mag[1]/mag.length() == 0.0); CHECK(mag[2]/mag.length() ==-1.0); - auto nvz = result.energy.nvz(); - Approx target = Approx(nvz).epsilon(1.e-10); + auto zeeman_ener = result.energy.zeeman_energy(); + Approx target = Approx(zeeman_ener).epsilon(1.e-10); hamiltonian::zeeman_coupling zc_(electrons.states().num_density_components()); basis::field> Bfield(electrons.spin_density().basis()); Bfield.fill(vector3 {0.0, 0.0, 0.0}); B.magnetic_field(0.0, Bfield); - auto nvz2 = 0.0; - zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), nvz2); - CHECK(nvz2 == target); + auto zeeman_ener2 = 0.0; + zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target); } SECTION("Spin non collinear zeeman calculation") { @@ -190,7 +189,7 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { ground_state::initial_guess(ions, electrons); perturbations::magnetic B{{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(10000).mixing(0.01), B); + 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), B); auto mag = observables::total_magnetization(electrons.spin_density()); auto mx = mag[0]/mag.length(); auto my = mag[1]/mag.length(); @@ -199,15 +198,15 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { CHECK(abs(my) < 1.e-7); CHECK(abs(mz + 1.0) < 1.e-7); - auto nvz = result.energy.nvz(); - Approx target = Approx(nvz).epsilon(1.e-10); + auto zeeman_ener = result.energy.zeeman_energy(); + Approx target = Approx(zeeman_ener).epsilon(1.e-10); hamiltonian::zeeman_coupling zc_(electrons.states().num_density_components()); basis::field> Bfield(electrons.spin_density().basis()); Bfield.fill(vector3 {0.0, 0.0, 0.0}); B.magnetic_field(0.0, Bfield); - auto nvz2 = 0.0; - zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), nvz2); - CHECK(nvz2 == target); + auto zeeman_ener2 = 0.0; + zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target); } } #endif \ No newline at end of file From 82da8027f3b2098c4488f239a8ea4daee2377aaf Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 23 Oct 2024 08:16:46 -0500 Subject: [PATCH 22/40] removed file because it takes too long to run the tests --- tests/ext_mag_field.cpp | 44 ----------------------------------------- 1 file changed, 44 deletions(-) delete mode 100644 tests/ext_mag_field.cpp diff --git a/tests/ext_mag_field.cpp b/tests/ext_mag_field.cpp deleted file mode 100644 index 235d63cc..00000000 --- a/tests/ext_mag_field.cpp +++ /dev/null @@ -1,44 +0,0 @@ -/* -*- indent-tabs-mode: t -*- */ - -#include -#include - -using namespace inq; -using namespace inq::magnitude; - -inq::utils::match match(3.0e-4); - -int main (int argc, char ** argv){ - auto env = input::environment{}; - - auto d = 1.21_A; - auto 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}); - - auto electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(10).spin_polarized()); - ground_state::initial_guess(ions, electrons); - perturbations::magnetic B{{0.0, 0.0, 0.4}}; - auto result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.1), B); - auto mag = observables::total_magnetization(electrons.spin_density()); - std::cout << mag << std::endl; - match.check("magnetization direction", mag/mag.length(), {0.0, 0.0, 1.0}); - - perturbations::magnetic B2{{0.0, 0.0, -0.4}}; - result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.1), B2); - auto mag2 = observables::total_magnetization(electrons.spin_density()); - std::cout << mag2 << std::endl; - match.check("magnetization direction", mag2/mag2.length(), {0.0, 0.0, -1.0}); - - perturbations::magnetic B3{{0.001, 0.0, -0.4}}; - electrons = systems::electrons(env.par(), ions, options::electrons{}.cutoff(30.0_Ha).extra_states(10).spin_non_collinear()); - ground_state::initial_guess(ions, electrons); - result = ground_state::calculate(ions, electrons, options::theory{}.lda(), inq::options::ground_state{}.steepest_descent().energy_tolerance(1.e-8_Ha).max_steps(1000).mixing(0.01), B3); - auto mag3 = observables::total_magnetization(electrons.spin_density()); - std::cout << mag3 << std::endl; - vector3 e_v = {0.001, 0.0, -0.4}; - e_v = e_v/sqrt(norm(e_v)); - match.check("magnetization direction", mag3/mag3.length(), e_v); - - return match.fail(); -} \ No newline at end of file From 44568391b4f62b6064c704a9e4d339974ad3bfc7 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Wed, 13 Nov 2024 10:22:13 -0600 Subject: [PATCH 23/40] replaced B with uniform_magnetic --- src/hamiltonian/self_consistency.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/hamiltonian/self_consistency.hpp b/src/hamiltonian/self_consistency.hpp index 5ef33cd2..c7dbc824 100644 --- a/src/hamiltonian/self_consistency.hpp +++ b/src/hamiltonian/self_consistency.hpp @@ -159,12 +159,12 @@ class self_consistency { // THE MAGNETIC FIELD if (pert_.has_magnetic_field()) { - basis::field> B(spin_density.basis()); - B.fill(vector3 {0.0, 0.0, 0.0}); - pert_.magnetic_field(time, B); + 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, B, hamiltonian.scalar_potential_, zeeman_ener); + zc_(spin_density, uniform_magnetic, hamiltonian.scalar_potential_, zeeman_ener); energy.zeeman_energy(zeeman_ener); } From da4459f255b5bc212cfe2da0bd244bbde37e2a16 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Wed, 13 Nov 2024 10:25:06 -0600 Subject: [PATCH 24/40] moved out of the class the test functions - replaced B with magnetic_uniform and Bfield with mag_field - removed process_potential - modified the if condition inside compute_vz --- src/hamiltonian/zeeman_coupling.hpp | 169 +++++++++++++--------------- 1 file changed, 77 insertions(+), 92 deletions(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 9433412f..e1006d98 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -24,42 +24,40 @@ class zeeman_coupling { //////////////////////////////////////////////////////////////////////////////////////////// - template - void operator()(SpinDensityType const & spin_density, MagneticField const & B, VKSType & vks, double & zeeman_ener) const { + template + void operator()(SpinDensityType const & spin_density, basis::field> const & magnetic_field, VKSType & vks, double & zeeman_ener) const { basis::field_set vz(vks.skeleton()); vz.fill(0.0); assert(vz.set_size() == spin_components_); - compute_vz(B, vz); + compute_vz(magnetic_field, vz); - process_potential(vz, vks); + gpu::run(vz.local_set_size(), vz.basis().local_size(), + [v = begin(vz.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip) { + vk[ip][is] += v[ip][is]; + }); zeeman_ener += compute_zeeman_energy(spin_density, vz); } //////////////////////////////////////////////////////////////////////////////////////////// - template - void compute_vz(MagneticField const & B, VZType & vz) const { + template + void compute_vz(basis::field> const & magnetic_field, VZType & vz) const { + gpu::run(vz.basis().local_size(), + [v = begin(vz.matrix()), magnetic_ = begin(magnetic_field.linear())] GPU_LAMBDA (auto ip) { + v[ip][0] +=-magnetic_[ip][2]; + v[ip][1] += magnetic_[ip][2]; + }); if (vz.set_size() == 4) { - gpu::run(vz.basis().local_size(), - [v = begin(vz.matrix()), b = begin(B.linear())] GPU_LAMBDA (auto ip) { - v[ip][0] +=-b[ip][2]; - v[ip][1] += b[ip][2]; - v[ip][2] +=-b[ip][0]; - v[ip][3] +=-b[ip][1]; - }); - } - else { - assert(vz.set_size() == 2); - gpu::run(vz.basis().local_size(), - [v = begin(vz.matrix()), b = begin(B.linear())] GPU_LAMBDA (auto ip) { - v[ip][0] +=-b[ip][2]; - v[ip][1] += b[ip][2]; - }); + gpu::run(vz.basis().local_size(), + [v = begin(vz.matrix()), magnetic_ = begin(magnetic_field.linear())] GPU_LAMBDA (auto ip) { + v[ip][2] +=-magnetic_[ip][0]; + v[ip][3] +=-magnetic_[ip][1]; + }); } } @@ -78,76 +76,65 @@ class zeeman_coupling { zeeman_ener_ += operations::integral_product_sum(spin_density, vz); return zeeman_ener_; } +}; - //////////////////////////////////////////////////////////////////////////////////////////// - - template - void process_potential(VZType const & vz, VKSType & vks) const { +} +} +#endif - gpu::run(vz.local_set_size(), vz.basis().local_size(), - [v = begin(vz.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip) { - vk[ip][is] += v[ip][is]; - }); - } +#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 & vz, RFType & rfield) { +template +void compute_psi_vz_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VZType const & vz, RFType & rfield) { - assert(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); - assert(std::get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); + assert(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); + assert(std::get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); - if (vz.set_size() == 2){ - gpu::run(phi.local_set_size(), phi.basis().local_size(), - [ph = begin(phi.matrix()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { - rf[ip] += occ[ist]*v[ip][spi]*norm(ph[ip][ist]); - }); - } - else { - assert(vz.set_size() == 4); - gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), - [ph = begin(phi.spinor_array()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { - auto offdiag = v[ip][2] + complex{0.0, 1.0}*v[ip][3]; - auto cross = 2.0*occ[ist]*real(offdiag*ph[ip][1][ist]*conj(ph[ip][0][ist])); - rf[ip] += occ[ist]*v[ip][0]*norm(ph[ip][0][ist]); - rf[ip] += occ[ist]*v[ip][1]*norm(ph[ip][1][ist]); - rf[ip] += cross; - }); - } + if (vz.set_size() == 2){ + gpu::run(phi.local_set_size(), phi.basis().local_size(), + [ph = begin(phi.matrix()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { + rf[ip] += occ[ist]*v[ip][spi]*norm(ph[ip][ist]); + }); } + else { + assert(vz.set_size() == 4); + gpu::run(phi.local_spinor_set_size(), phi.basis().local_size(), + [ph = begin(phi.spinor_array()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { + auto offdiag = v[ip][2] + complex{0.0, 1.0}*v[ip][3]; + auto cross = 2.0*occ[ist]*real(offdiag*ph[ip][1][ist]*conj(ph[ip][0][ist])); + rf[ip] += occ[ist]*v[ip][0]*norm(ph[ip][0][ist]); + rf[ip] += occ[ist]*v[ip][1]*norm(ph[ip][1][ist]); + rf[ip] += cross; + }); + } +} - //////////////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////////////// - template - void eval_psi_vz_psi(CommType & comm, SpinDensityType const & spin_density, MagneticField const & B, occupations_array_type const & occupations, kpin_type const & kpin, double & zeeman_ener) { +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 vz(spin_density.skeleton()); - vz.fill(0.0); - compute_vz(B, vz); - - basis::field rfield(vz.basis()); - rfield.fill(0.0); - int iphi = 0; - for (auto & phi : kpin) { - compute_psi_vz_psi_ofr(occupations[iphi], phi, vz, rfield); - iphi++; - } - - rfield.all_reduce(comm); - zeeman_ener += operations::integral(rfield); + basis::field_set vz(spin_density.skeleton()); + vz.fill(0.0); + hamiltonian::zeeman_coupling zc_(spin_density.set_size()); + zc_.compute_vz(magnetic_field, vz); + + basis::field rfield(vz.basis()); + rfield.fill(0.0); + int iphi = 0; + for (auto & phi : kpin) { + compute_psi_vz_psi_ofr(occupations[iphi], phi, vz, rfield); + iphi++; } -}; + rfield.all_reduce(comm); + zeeman_ener += operations::integral(rfield); } -} -#endif - -#ifdef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST -#undef INQ_HAMILTONIAN_ZEEMAN_COUPLING_UNIT_TEST - -#include -#include TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { @@ -163,8 +150,8 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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 B{{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), B); + 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(mag[0]/mag.length() == 0.0); CHECK(mag[1]/mag.length() == 0.0); @@ -172,12 +159,11 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { auto zeeman_ener = result.energy.zeeman_energy(); Approx target = Approx(zeeman_ener).epsilon(1.e-10); - hamiltonian::zeeman_coupling zc_(electrons.states().num_density_components()); - basis::field> Bfield(electrons.spin_density().basis()); - Bfield.fill(vector3 {0.0, 0.0, 0.0}); - B.magnetic_field(0.0, Bfield); + 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; - zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), zeeman_ener2); + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); CHECK(zeeman_ener2 == target); } @@ -187,9 +173,9 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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 B{{0.0, 0.0, -1.0}}; + 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), B); + 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()); auto mx = mag[0]/mag.length(); auto my = mag[1]/mag.length(); @@ -200,12 +186,11 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { auto zeeman_ener = result.energy.zeeman_energy(); Approx target = Approx(zeeman_ener).epsilon(1.e-10); - hamiltonian::zeeman_coupling zc_(electrons.states().num_density_components()); - basis::field> Bfield(electrons.spin_density().basis()); - Bfield.fill(vector3 {0.0, 0.0, 0.0}); - B.magnetic_field(0.0, Bfield); + 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; - zc_.eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), Bfield, electrons.occupations(), electrons.kpin(), zeeman_ener2); + eval_psi_vz_psi(electrons.kpin_states_comm(), electrons.spin_density(), mag_field, electrons.occupations(), electrons.kpin(), zeeman_ener2); CHECK(zeeman_ener2 == target); } } From df1037564311ab6bfb3d1379a16451f2f29dc028 Mon Sep 17 00:00:00 2001 From: jacopos86 Date: Wed, 13 Nov 2024 10:26:09 -0600 Subject: [PATCH 25/40] changed B -> uniform_magnetic and Bfield -> mag_field --- src/perturbations/magnetic.hpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/perturbations/magnetic.hpp b/src/perturbations/magnetic.hpp index a8ccee85..3159a6e4 100644 --- a/src/perturbations/magnetic.hpp +++ b/src/perturbations/magnetic.hpp @@ -14,8 +14,8 @@ class magnetic : public none { public: - magnetic(vector3 const & B): - magnetic_vector_(B) + magnetic(vector3 const & value): + magnetic_vector_(value) { } @@ -26,8 +26,8 @@ class magnetic : public none { template void magnetic_field(const double time, MagneticField & magnetic) const { gpu::run(magnetic.basis().local_size(), - [b = begin(magnetic.linear()), mv = magnetic_vector_] GPU_LAMBDA (auto ip){ - b[ip] += mv; + [magnetic_ = begin(magnetic.linear()), mv = magnetic_vector_] GPU_LAMBDA (auto ip){ + magnetic_[ip] += mv; }); } @@ -49,18 +49,18 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { parallel::communicator comm{boost::mpi3::environment::get_world_instance()}; - perturbations::magnetic B{{0.0, 0.0, 1.0}}; + 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> Bfield(bas); - Bfield.fill(vector3 {0.0, 0.0, 0.0}); + basis::field> mag_field(bas); + mag_field.fill(vector3 {0.0, 0.0, 0.0}); - CHECK(B.has_magnetic_field()); - B.magnetic_field(/*time*/ 0.0, Bfield); - CHECK(Bfield.linear()[0] == vector3{0.0, 0.0, 1.0}); - CHECK(Bfield.linear()[1] == vector3{0.0, 0.0, 1.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}); - B.magnetic_field(/*time*/ 1000.0, Bfield); - CHECK(Bfield.linear()[0] == vector3{0.0, 0.0, 2.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 From 80a6a6ace2bf10bfc99bf1e0b8bc4577e332566e Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 07:12:20 -0600 Subject: [PATCH 26/40] replaced integral_sum with integral_partial_sum inside normalize --- src/observables/density.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/observables/density.hpp b/src/observables/density.hpp index d2b819e2..af3e7c6e 100644 --- a/src/observables/density.hpp +++ b/src/observables/density.hpp @@ -102,7 +102,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(), From e649559f60fdd8252f15ed8cb62d8ae28b69f854 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 07:13:29 -0600 Subject: [PATCH 27/40] defined integral_partial_sum to sum the density field components only up to a certain index --- src/operations/integral.hpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) 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; From a82dff7ad7eb11687387d10ced9c4ee8c0c0c072 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 08:04:08 -0600 Subject: [PATCH 28/40] corrected magnetization along y sign --- src/observables/magnetization.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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; From 709604d309284528ac61b701d1f06f2be66da551 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 08:05:45 -0600 Subject: [PATCH 29/40] added rotate_total_magnetization function - this is used to initialize the magnetization along a certain direction in the initial_guess function --- src/observables/density.hpp | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/src/observables/density.hpp b/src/observables/density.hpp index af3e7c6e..bc1ed6ea 100644 --- a/src/observables/density.hpp +++ b/src/observables/density.hpp @@ -15,6 +15,7 @@ #include #include #include +#include namespace inq { namespace observables { @@ -133,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; + }); + } +} + } } } From a17d69cba1ade039b229306e40ed62fa5e1fa21b Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 08:09:39 -0600 Subject: [PATCH 30/40] added the magnet_dir optional variable in case we initialize the magnetization along a given direction + series of tests to see if the magnetization direction and total charge is correct after initialization --- src/ground_state/initial_guess.hpp | 123 ++++++++++++++++++++++++++++- 1 file changed, 122 insertions(+), 1 deletion(-) 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 From c26ac9bdb16410f530160a9f34112251ef52ea31 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 09:43:46 -0600 Subject: [PATCH 31/40] compute_psi_vxc_psi_ofr and eval_psi_vxc_psi moved out of the class because used only in tests - added functionals() to get the list of functionals that is used inside eval_psi_vxc_psi - in compute_vxc I changed the name of few internal variables inside the GPU call - inside compute_nvxc I corrected the sign of the y component --- src/hamiltonian/xc_term.hpp | 187 +++++++++++++++++++----------------- 1 file changed, 97 insertions(+), 90 deletions(-) diff --git a/src/hamiltonian/xc_term.hpp b/src/hamiltonian/xc_term.hpp index 90d56c2f..b9974688 100644 --- a/src/hamiltonian/xc_term.hpp +++ b/src/hamiltonian/xc_term.hpp @@ -106,9 +106,12 @@ class xc_term { auto nvxc_ = 0.0; if (spin_density.set_size() == 4) { gpu::run(spin_density.local_set_size(), spin_density.basis().local_size(), - [v = begin(vxc.matrix())] GPU_LAMBDA (auto is, auto ip){ - if (is >= 2){ - v[ip][is] = 2.0*v[ip][is]; + [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]; } }); } @@ -157,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; } }); } @@ -185,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){ @@ -322,6 +252,12 @@ class xc_term { //////////////////////////////////////////////////////////////////////////////////////////// + auto & functionals() const { + return functionals_; + } + + //////////////////////////////////////////////////////////////////////////////////////////// + }; } } @@ -332,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(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); + assert(std::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){ @@ -510,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); } @@ -531,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); } @@ -552,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); } From 81493a5b88657d848c2b0088a40ffa1b1965fac7 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 11:12:33 -0600 Subject: [PATCH 32/40] changed variable name vz with zeeman_pot -> this is the zeeman potential added to vks --- src/hamiltonian/zeeman_coupling.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index e1006d98..45169587 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -27,19 +27,19 @@ class zeeman_coupling { template void operator()(SpinDensityType const & spin_density, basis::field> const & magnetic_field, VKSType & vks, double & zeeman_ener) const { - basis::field_set vz(vks.skeleton()); - vz.fill(0.0); + basis::field_set zeeman_pot(vks.skeleton()); + zeeman_pot.fill(0.0); - assert(vz.set_size() == spin_components_); + assert(zeeman_pot.set_size() == spin_components_); - compute_vz(magnetic_field, vz); + compute_vz(magnetic_field, zeeman_pot); - gpu::run(vz.local_set_size(), vz.basis().local_size(), - [v = begin(vz.matrix()), vk = begin(vks.matrix())] GPU_LAMBDA (auto is, auto ip) { - vk[ip][is] += v[ip][is]; + 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, vz); + zeeman_ener += compute_zeeman_energy(spin_density, zeeman_pot); } //////////////////////////////////////////////////////////////////////////////////////////// From 2faef75553650ec73bd5c5d30ebff4e55da8c093 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:22:23 -0600 Subject: [PATCH 33/40] replaced vz everywhere with zeeman_pot - vz used inside the GPU call for zeeman_pot matrix - replaced compute_vz with compute_zeeman_pot - inside compute_zeeman_energy correction of the sign for the y component - modified the CHECK inside unit tests to use margin --- src/hamiltonian/zeeman_coupling.hpp | 80 +++++++++++++++-------------- 1 file changed, 41 insertions(+), 39 deletions(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 45169587..3aeb1211 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -32,7 +32,7 @@ class zeeman_coupling { assert(zeeman_pot.set_size() == spin_components_); - compute_vz(magnetic_field, zeeman_pot); + 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) { @@ -45,18 +45,18 @@ class zeeman_coupling { //////////////////////////////////////////////////////////////////////////////////////////// template - void compute_vz(basis::field> const & magnetic_field, VZType & vz) const { + void compute_zeeman_potential(basis::field> const & magnetic_field, VZType & zeeman_pot) const { - gpu::run(vz.basis().local_size(), - [v = begin(vz.matrix()), magnetic_ = begin(magnetic_field.linear())] GPU_LAMBDA (auto ip) { - v[ip][0] +=-magnetic_[ip][2]; - v[ip][1] += magnetic_[ip][2]; + 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 (vz.set_size() == 4) { - gpu::run(vz.basis().local_size(), - [v = begin(vz.matrix()), magnetic_ = begin(magnetic_field.linear())] GPU_LAMBDA (auto ip) { - v[ip][2] +=-magnetic_[ip][0]; - v[ip][3] +=-magnetic_[ip][1]; + 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]; }); } } @@ -64,16 +64,21 @@ class zeeman_coupling { //////////////////////////////////////////////////////////////////////////////////////////// template - double compute_zeeman_energy(SpinDensityType const & spin_density, VZType & vz) const { + 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(), - [v = begin(vz.matrix())] GPU_LAMBDA (auto is, auto ip) { - if (is >= 2) v[ip][is] = 2.0*v[ip][is]; + [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, vz); + zeeman_ener_ += operations::integral_product_sum(spin_density, zeeman_pot); return zeeman_ener_; } }; @@ -90,25 +95,25 @@ class zeeman_coupling { using namespace inq; template -void compute_psi_vz_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VZType const & vz, RFType & rfield) { +void compute_psi_vz_psi_ofr(occupations_array_type const & occupations, field_set_type const & phi, VZType const & zeeman_pot, RFType & rfield) { assert(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); assert(std::get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); - if (vz.set_size() == 2){ + if (zeeman_pot.set_size() == 2){ gpu::run(phi.local_set_size(), phi.basis().local_size(), - [ph = begin(phi.matrix()), rf = begin(rfield.linear()), v = begin(vz.matrix()), occ = begin(occupations), spi = phi.spin_index()] GPU_LAMBDA (auto ist, auto ip) { - rf[ip] += occ[ist]*v[ip][spi]*norm(ph[ip][ist]); + [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(vz.set_size() == 4); + 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()), v = begin(vz.matrix()), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto ip) { - auto offdiag = v[ip][2] + complex{0.0, 1.0}*v[ip][3]; - auto cross = 2.0*occ[ist]*real(offdiag*ph[ip][1][ist]*conj(ph[ip][0][ist])); - rf[ip] += occ[ist]*v[ip][0]*norm(ph[ip][0][ist]); - rf[ip] += occ[ist]*v[ip][1]*norm(ph[ip][1][ist]); + [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; }); } @@ -119,16 +124,16 @@ void compute_psi_vz_psi_ofr(occupations_array_type const & occupations, field_se 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 vz(spin_density.skeleton()); - vz.fill(0.0); + basis::field_set zeeman_pot(spin_density.skeleton()); + zeeman_pot.fill(0.0); hamiltonian::zeeman_coupling zc_(spin_density.set_size()); - zc_.compute_vz(magnetic_field, vz); + zc_.compute_zeeman_potential(magnetic_field, zeeman_pot); - basis::field rfield(vz.basis()); + 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, vz, rfield); + compute_psi_vz_psi_ofr(occupations[iphi], phi, zeeman_pot, rfield); iphi++; } @@ -140,6 +145,7 @@ 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()}; @@ -153,9 +159,9 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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(mag[0]/mag.length() == 0.0); - CHECK(mag[1]/mag.length() == 0.0); - CHECK(mag[2]/mag.length() ==-1.0); + 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); @@ -177,12 +183,8 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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()); - auto mx = mag[0]/mag.length(); - auto my = mag[1]/mag.length(); - auto mz = mag[2]/mag.length(); - CHECK(abs(mx) < 1.e-7); - CHECK(abs(my) < 1.e-7); - CHECK(abs(mz + 1.0) < 1.e-7); + 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); From e07f2bf9ea93f838086d94a916fabd1778709bc5 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:33:38 -0600 Subject: [PATCH 34/40] included test with magnetic field oriented along x-y -> check both final magnetic orientation and zeeman energy --- src/hamiltonian/zeeman_coupling.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 3aeb1211..a1ce7270 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -194,6 +194,23 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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); } } #endif \ No newline at end of file From 9a4d97015b8345cfcc1a879ba200b083b0df0283 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:36:38 -0600 Subject: [PATCH 35/40] added test with external magnetic field (1,-1,0) --- src/hamiltonian/zeeman_coupling.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index a1ce7270..642d7069 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -211,6 +211,23 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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(), bfield, electrons.occupations(), electrons.kpin(), zeeman_ener2); + CHECK(zeeman_ener2 == target3); } } #endif \ No newline at end of file From 553162db0030b1d8069e1702befd8ad67d3315f8 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:39:09 -0600 Subject: [PATCH 36/40] added test case with external magnetic field along (1,1,1) --- src/hamiltonian/zeeman_coupling.hpp | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 642d7069..48af77db 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -226,8 +226,25 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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(), bfield, electrons.occupations(), electrons.kpin(), zeeman_ener2); + 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); } } #endif \ No newline at end of file From 19f1cd2783554b5bb6295292dede076ab6b91506 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:41:01 -0600 Subject: [PATCH 37/40] added test with external magnetic field along (0,-1,1) --- src/hamiltonian/zeeman_coupling.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 48af77db..8fab2d90 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -245,6 +245,23 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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); } } #endif \ No newline at end of file From 1d45e9b1bf29d30d8b09e19ec4b96f4ea16a6892 Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:43:12 -0600 Subject: [PATCH 38/40] added test with external magnetic field (1,-2,1.5) --- src/hamiltonian/zeeman_coupling.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 8fab2d90..b476fdea 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -262,6 +262,23 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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); } } #endif \ No newline at end of file From 3ea6669c0d270cf49ef46c6fb20be20b4e49a5bf Mon Sep 17 00:00:00 2001 From: Jacopo Simoni Date: Wed, 20 Nov 2024 13:44:47 -0600 Subject: [PATCH 39/40] added test with external magnetic field (4,-2,1) --- src/hamiltonian/zeeman_coupling.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index b476fdea..3f6bc0ec 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -279,6 +279,23 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { 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 \ No newline at end of file From d70ccc22b50296f32a6dffac065345ada4d6d6b6 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Wed, 11 Dec 2024 10:49:08 -0800 Subject: [PATCH 40/40] Replaced std::get by get. --- src/hamiltonian/xc_term.hpp | 4 ++-- src/hamiltonian/zeeman_coupling.hpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/hamiltonian/xc_term.hpp b/src/hamiltonian/xc_term.hpp index b9974688..d2176b61 100644 --- a/src/hamiltonian/xc_term.hpp +++ b/src/hamiltonian/xc_term.hpp @@ -273,8 +273,8 @@ 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(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); - assert(std::get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); + 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(), diff --git a/src/hamiltonian/zeeman_coupling.hpp b/src/hamiltonian/zeeman_coupling.hpp index 3f6bc0ec..aaf5fad6 100644 --- a/src/hamiltonian/zeeman_coupling.hpp +++ b/src/hamiltonian/zeeman_coupling.hpp @@ -97,8 +97,8 @@ 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(std::get<1>(sizes(phi.spinor_array())) == phi.spinor_dim()); - assert(std::get<2>(sizes(phi.spinor_array())) == phi.local_spinor_set_size()); + 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(), @@ -298,4 +298,4 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) { CHECK(zeeman_ener2 == target7); } } -#endif \ No newline at end of file +#endif