Skip to content

Commit

Permalink
Merge branch 'ext_B_field' into 'master'
Browse files Browse the repository at this point in the history 
External magnetic field and Zeeman coupling

See merge request npneq/inq!1156
  • Loading branch information
@xavierandrade
xavierandrade committed Dec 11, 2024
2 parents 972817a + d70ccc2 commit 0bce1cd
Show file tree
Hide file tree
Showing 9 changed files with 671 additions and 93 deletions.
123 changes: 122 additions & 1 deletion src/ground_state/initial_guess.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<vector3<double>> const & magnet_dir = {}){

int iphi = 0;
for(auto & phi : electrons.kpin()) {
Expand All @@ -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_);
}

}
}
Expand All @@ -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

12 changes: 12 additions & 0 deletions src/hamiltonian/energy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ namespace hamiltonian {
double xc_ = 0.0;
double nvxc_ = 0.0;
double exact_exchange_ = 0.0;
double zeeman_ener_ = 0.0;

#ifdef ENABLE_CUDA
public:
Expand Down Expand Up @@ -149,6 +150,14 @@ namespace hamiltonian {
nvxc_ = val;
}

auto & zeeman_energy() const {
return zeeman_ener_;
}

void zeeman_energy(double const & val) {
zeeman_ener_ = val;
}

auto & exact_exchange() const {
return exact_exchange_;
}
Expand Down Expand Up @@ -186,6 +195,7 @@ namespace hamiltonian {
utils::save_value(comm, dirname + "/xc", xc_, error_message);
utils::save_value(comm, dirname + "/nvxc", nvxc_, error_message);
utils::save_value(comm, dirname + "/exact_exchange", exact_exchange_, error_message);
utils::save_value(comm, dirname + "/zeeman_energy", zeeman_ener_, error_message);
}

static auto load(std::string const & dirname) {
Expand All @@ -201,6 +211,7 @@ namespace hamiltonian {
utils::load_value(dirname + "/xc", en.xc_, error_message);
utils::load_value(dirname + "/nvxc", en.nvxc_, error_message);
utils::load_value(dirname + "/exact_exchange", en.exact_exchange_, error_message);
utils::load_value(dirname + "/zeeman_energy", en.zeeman_ener_, error_message);

return en;
}
Expand All @@ -219,6 +230,7 @@ namespace hamiltonian {
tfm::format(out, " nvxc = %20.12f Ha\n", self.nvxc());
tfm::format(out, " exact-exchange = %20.12f Ha\n", self.exact_exchange());
tfm::format(out, " ion = %20.12f Ha\n", self.ion());
tfm::format(out, " zeeman-energy = %20.12f Ha\n", self.zeeman_energy());
tfm::format(out, "\n");

return out;
Expand Down
14 changes: 12 additions & 2 deletions src/hamiltonian/self_consistency.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include <operations/integral.hpp>
#include <hamiltonian/xc_term.hpp>
#include <hamiltonian/atomic_potential.hpp>
#include <hamiltonian/zeeman_coupling.hpp>
#include <options/theory.hpp>
#include <perturbations/none.hpp>
#include <solvers/velocity_verlet.hpp>
Expand Down Expand Up @@ -155,8 +156,17 @@ class self_consistency {
hamiltonian.uniform_vector_potential_ += induced_vector_potential_;
}

//the magnetic perturbation is not implemented yet
assert(not pert_.has_magnetic_field());
// THE MAGNETIC FIELD

if (pert_.has_magnetic_field()) {
basis::field<basis::real_space, vector3<double>> uniform_magnetic(spin_density.basis());
uniform_magnetic.fill(vector3 {0.0, 0.0, 0.0});
pert_.magnetic_field(time, uniform_magnetic);
zeeman_coupling zc_(spin_density.set_size());
auto zeeman_ener = 0.0;
zc_(spin_density, uniform_magnetic, hamiltonian.scalar_potential_, zeeman_ener);
energy.zeeman_energy(zeeman_ener);
}

}

Expand Down
Loading

0 comments on commit 0bce1cd

Please sign in to comment.