From 5270c954c2dbf2aeebe66cd6e0091034db1d936c Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Sat, 28 Sep 2024 01:09:30 -0700 Subject: [PATCH 1/3] Move forces to observables:: --- src/ground_state/calculator.hpp | 4 ++-- src/{hamiltonian => observables}/forces.hpp | 10 +++++----- src/real_time/propagate.hpp | 8 ++++---- src/real_time/viewables.hpp | 6 +++--- 4 files changed, 14 insertions(+), 14 deletions(-) rename src/{hamiltonian => observables}/forces.hpp (95%) diff --git a/src/ground_state/calculator.hpp b/src/ground_state/calculator.hpp index ddb54fd4..2f6d02c5 100644 --- a/src/ground_state/calculator.hpp +++ b/src/ground_state/calculator.hpp @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -26,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -246,7 +246,7 @@ class calculator { auto normres = res.energy.calculate(ham_, electrons); if(solver_.calc_forces() and electrons.states().spinor_dim() == 1) { - res.forces = hamiltonian::calculate_forces(ions_, electrons, ham_); + res.forces = observables::calculate_forces(ions_, electrons, ham_); } if(solver_.calc_forces() and electrons.states().spinor_dim() == 2) { diff --git a/src/hamiltonian/forces.hpp b/src/observables/forces.hpp similarity index 95% rename from src/hamiltonian/forces.hpp rename to src/observables/forces.hpp index de6b35a9..dafa1f37 100644 --- a/src/hamiltonian/forces.hpp +++ b/src/observables/forces.hpp @@ -1,7 +1,7 @@ /* -*- indent-tabs-mode: t -*- */ -#ifndef INQ__HAMILTONIAN__FORCES -#define INQ__HAMILTONIAN__FORCES +#ifndef INQ__OBSERVABLES__FORCES +#define INQ__OBSERVABLES__FORCES // Copyright (C) 2019-2023 Lawrence Livermore National Security, LLC., Xavier Andrade, Alfredo A. Correa // @@ -17,7 +17,7 @@ #include namespace inq { -namespace hamiltonian { +namespace observables { template struct loc_pot { @@ -99,8 +99,8 @@ gpu::array, 1> calculate_forces(const systems::ions & ions, syst } #endif -#ifdef INQ_HAMILTONIAN_FORCES_UNIT_TEST -#undef INQ_HAMILTONIAN_FORCES_UNIT_TEST +#ifdef INQ_OBSERVABLES_FORCES_UNIT_TEST +#undef INQ_OBSERVABLES_FORCES_UNIT_TEST #include #include diff --git a/src/real_time/propagate.hpp b/src/real_time/propagate.hpp index a1898f58..19419517 100644 --- a/src/real_time/propagate.hpp +++ b/src/real_time/propagate.hpp @@ -10,9 +10,9 @@ #include #include -#include #include #include +#include #include #include #include @@ -72,8 +72,8 @@ void propagate(systems::ions & ions, systems::electrons & electrons, ProcessFunc energy.calculate(ham, electrons); energy.ion(ionic::interaction_energy(ions.cell(), ions, electrons.atomic_pot())); - auto forces = decltype(hamiltonian::calculate_forces(ions, electrons, ham)){}; - if(ion_propagator.needs_force()) forces = hamiltonian::calculate_forces(ions, electrons, ham); + auto forces = decltype(observables::calculate_forces(ions, electrons, ham)){}; + if(ion_propagator.needs_force()) forces = observables::calculate_forces(ions, electrons, ham); auto current = vector3{0.0, 0.0, 0.0}; if(sc.has_induced_vector_potential()) current = observables::current(ions, electrons, ham); @@ -98,7 +98,7 @@ void propagate(systems::ions & ions, systems::electrons & electrons, ProcessFunc energy.calculate(ham, electrons); - if(ion_propagator.needs_force()) forces = hamiltonian::calculate_forces(ions, electrons, ham); + if(ion_propagator.needs_force()) forces = observables::calculate_forces(ions, electrons, ham); //propagate ionic velocities to t + dt ion_propagator.propagate_velocities(dt, ions, forces); diff --git a/src/real_time/viewables.hpp b/src/real_time/viewables.hpp index 1027bf64..3b794216 100644 --- a/src/real_time/viewables.hpp +++ b/src/real_time/viewables.hpp @@ -10,10 +10,10 @@ // file, You can obtain one at https://mozilla.org/MPL/2.0/. #include -#include #include -#include #include +#include +#include #include #include #include @@ -73,7 +73,7 @@ class viewables { } auto forces() { - if(forces_.size() == 0) forces_ = hamiltonian::calculate_forces(ions_, electrons_, ham_); + if(forces_.size() == 0) forces_ = observables::calculate_forces(ions_, electrons_, ham_); return forces_; } From f64c51f8f42da73bf364c409ce7147bc62c747eb Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Sat, 28 Sep 2024 19:08:23 -0700 Subject: [PATCH 2/3] Define a structure for the forces and stress, do the calculation in the constructor. --- src/ground_state/calculator.hpp | 2 +- src/observables/forces.hpp | 119 ++++++++++++++++++-------------- src/real_time/propagate.hpp | 8 +-- src/real_time/viewables.hpp | 2 +- 4 files changed, 74 insertions(+), 57 deletions(-) diff --git a/src/ground_state/calculator.hpp b/src/ground_state/calculator.hpp index 2f6d02c5..b4e40aec 100644 --- a/src/ground_state/calculator.hpp +++ b/src/ground_state/calculator.hpp @@ -246,7 +246,7 @@ class calculator { auto normres = res.energy.calculate(ham_, electrons); if(solver_.calc_forces() and electrons.states().spinor_dim() == 1) { - res.forces = observables::calculate_forces(ions_, electrons, ham_); + res.forces = observables::forces_stress{ions_, electrons, ham_}.forces; } if(solver_.calc_forces() and electrons.states().spinor_dim() == 2) { diff --git a/src/observables/forces.hpp b/src/observables/forces.hpp index dafa1f37..56da5775 100644 --- a/src/observables/forces.hpp +++ b/src/observables/forces.hpp @@ -31,69 +31,86 @@ struct loc_pot { } }; -template -gpu::array, 1> calculate_forces(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){ +struct forces_stress { + gpu::array, 1> forces; + gpu::array stress; - CALI_CXX_MARK_FUNCTION; + forces_stress() = default; - basis::field> gdensity(electrons.density_basis()); - gdensity.fill(vector3{0.0, 0.0, 0.0}); - - gpu::array, 1> forces_non_local(ions.size(), {0.0, 0.0, 0.0}); - - auto iphi = 0; - for(auto & phi : electrons.kpin()){ - - auto gphi = operations::gradient(phi, /* factor = */ 1.0, /*shift = */ phi.kpoint()); - observables::density::calculate_gradient_add(electrons.occupations()[iphi], phi, gphi, gdensity); - - ham.projectors_all().force(phi, gphi, ions.cell().metric(), electrons.occupations()[iphi], ham.uniform_vector_potential(), forces_non_local); - - iphi++; + template + forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham): + stress({3, 3}, 0.0) + { + forces = calculate_forces(ions, electrons, ham); } - gdensity.all_reduce(electrons.kpin_states_comm()); +private: - if(electrons.full_comm().size() > 1){ - CALI_CXX_MARK_SCOPE("forces_nonlocal::reduce"); - electrons.full_comm().all_reduce_n(raw_pointer_cast(forces_non_local.data_elements()), forces_non_local.size(), std::plus<>{}); - } - - auto ionic_forces = ionic::interaction_forces(ions.cell(), ions, electrons.atomic_pot()); - - gpu::array, 1> forces_local(ions.size(), {0.0, 0.0, 0.0}); - - { CALI_CXX_MARK_SCOPE("forces_local"); + template + gpu::array, 1> calculate_forces(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){ - solvers::poisson poisson_solver; + CALI_CXX_MARK_FUNCTION; - //the force from the local potential - for(int iatom = 0; iatom < ions.size(); iatom++){ - auto ionic_long_range = poisson_solver(electrons.atomic_pot().ionic_density(electrons.states_comm(), electrons.density_basis(), ions, iatom)); - auto ionic_short_range = electrons.atomic_pot().local_potential(electrons.states_comm(), electrons.density_basis(), ions, iatom); - - auto force_cov = -gpu::run(gpu::reduce(electrons.density_basis().local_size()), - loc_pot - {begin(ionic_long_range.linear()), begin(ionic_short_range.linear()), begin(gdensity.linear())}); + basis::field> gdensity(electrons.density_basis()); + gdensity.fill(vector3{0.0, 0.0, 0.0}); + + gpu::array, 1> forces_non_local(ions.size(), {0.0, 0.0, 0.0}); + + auto iphi = 0; + for(auto & phi : electrons.kpin()){ - forces_local[iatom] = electrons.density_basis().volume_element()*ions.cell().metric().to_cartesian(force_cov); + auto gphi = operations::gradient(phi, /* factor = */ 1.0, /*shift = */ phi.kpoint()); + observables::density::calculate_gradient_add(electrons.occupations()[iphi], phi, gphi, gdensity); + + ham.projectors_all().force(phi, gphi, ions.cell().metric(), electrons.occupations()[iphi], ham.uniform_vector_potential(), forces_non_local); + + iphi++; } - - if(electrons.density_basis().comm().size() > 1){ - CALI_CXX_MARK_SCOPE("forces_local::reduce"); - electrons.density_basis().comm().all_reduce_n(reinterpret_cast(raw_pointer_cast(forces_local.data_elements())), 3*forces_local.size()); + + gdensity.all_reduce(electrons.kpin_states_comm()); + + if(electrons.full_comm().size() > 1){ + CALI_CXX_MARK_SCOPE("forces_nonlocal::reduce"); + electrons.full_comm().all_reduce_n(raw_pointer_cast(forces_non_local.data_elements()), forces_non_local.size(), std::plus<>{}); + } + + auto ionic_forces = ionic::interaction_forces(ions.cell(), ions, electrons.atomic_pot()); + + gpu::array, 1> forces_local(ions.size(), {0.0, 0.0, 0.0}); + + { CALI_CXX_MARK_SCOPE("forces_local"); + + solvers::poisson poisson_solver; + + //the force from the local potential + for(int iatom = 0; iatom < ions.size(); iatom++){ + auto ionic_long_range = poisson_solver(electrons.atomic_pot().ionic_density(electrons.states_comm(), electrons.density_basis(), ions, iatom)); + auto ionic_short_range = electrons.atomic_pot().local_potential(electrons.states_comm(), electrons.density_basis(), ions, iatom); + + auto force_cov = -gpu::run(gpu::reduce(electrons.density_basis().local_size()), + loc_pot + {begin(ionic_long_range.linear()), begin(ionic_short_range.linear()), begin(gdensity.linear())}); + + forces_local[iatom] = electrons.density_basis().volume_element()*ions.cell().metric().to_cartesian(force_cov); + } + + if(electrons.density_basis().comm().size() > 1){ + CALI_CXX_MARK_SCOPE("forces_local::reduce"); + electrons.density_basis().comm().all_reduce_n(reinterpret_cast(raw_pointer_cast(forces_local.data_elements())), 3*forces_local.size()); + } } + + gpu::array, 1> forces(ions.size()); + for(int iatom = 0; iatom < ions.size(); iatom++){ + forces[iatom] = ionic_forces[iatom] + forces_local[iatom] + forces_non_local[iatom]; + } + + // MISSING: the non-linear core correction term + + return forces; } - - gpu::array, 1> forces(ions.size()); - for(int iatom = 0; iatom < ions.size(); iatom++){ - forces[iatom] = ionic_forces[iatom] + forces_local[iatom] + forces_non_local[iatom]; - } - // MISSING: the non-linear core correction term - - return forces; -} +}; } } diff --git a/src/real_time/propagate.hpp b/src/real_time/propagate.hpp index 19419517..8010506a 100644 --- a/src/real_time/propagate.hpp +++ b/src/real_time/propagate.hpp @@ -72,8 +72,8 @@ void propagate(systems::ions & ions, systems::electrons & electrons, ProcessFunc energy.calculate(ham, electrons); energy.ion(ionic::interaction_energy(ions.cell(), ions, electrons.atomic_pot())); - auto forces = decltype(observables::calculate_forces(ions, electrons, ham)){}; - if(ion_propagator.needs_force()) forces = observables::calculate_forces(ions, electrons, ham); + auto forces = decltype(observables::forces_stress{ions, electrons, ham}.forces){}; + if(ion_propagator.needs_force()) forces = observables::forces_stress{ions, electrons, ham}.forces; auto current = vector3{0.0, 0.0, 0.0}; if(sc.has_induced_vector_potential()) current = observables::current(ions, electrons, ham); @@ -98,8 +98,8 @@ void propagate(systems::ions & ions, systems::electrons & electrons, ProcessFunc energy.calculate(ham, electrons); - if(ion_propagator.needs_force()) forces = observables::calculate_forces(ions, electrons, ham); - + if(ion_propagator.needs_force()) forces = observables::forces_stress{ions, electrons, ham}.forces; + //propagate ionic velocities to t + dt ion_propagator.propagate_velocities(dt, ions, forces); diff --git a/src/real_time/viewables.hpp b/src/real_time/viewables.hpp index 3b794216..ef6fe314 100644 --- a/src/real_time/viewables.hpp +++ b/src/real_time/viewables.hpp @@ -73,7 +73,7 @@ class viewables { } auto forces() { - if(forces_.size() == 0) forces_ = observables::calculate_forces(ions_, electrons_, ham_); + if(forces_.size() == 0) forces_ = observables::forces_stress{ions_, electrons_, ham_}.forces; return forces_; } From 410ab84b8879eb1ef4d14cde7ad67cbae15abb17 Mon Sep 17 00:00:00 2001 From: Xavier Andrade Date: Sun, 29 Sep 2024 22:29:40 -0700 Subject: [PATCH 3/3] Renamed forces.hpp to forces_stress.hpp. --- src/ground_state/calculator.hpp | 2 +- .../{forces.hpp => forces_stress.hpp} | 21 +++++++++---------- src/real_time/propagate.hpp | 2 +- src/real_time/viewables.hpp | 2 +- 4 files changed, 13 insertions(+), 14 deletions(-) rename src/observables/{forces.hpp => forces_stress.hpp} (89%) diff --git a/src/ground_state/calculator.hpp b/src/ground_state/calculator.hpp index b4e40aec..a308ce11 100644 --- a/src/ground_state/calculator.hpp +++ b/src/ground_state/calculator.hpp @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/src/observables/forces.hpp b/src/observables/forces_stress.hpp similarity index 89% rename from src/observables/forces.hpp rename to src/observables/forces_stress.hpp index 56da5775..f7a7f3dc 100644 --- a/src/observables/forces.hpp +++ b/src/observables/forces_stress.hpp @@ -1,9 +1,9 @@ /* -*- indent-tabs-mode: t -*- */ -#ifndef INQ__OBSERVABLES__FORCES -#define INQ__OBSERVABLES__FORCES +#ifndef INQ__OBSERVABLES__FORCES_STRESS +#define INQ__OBSERVABLES__FORCES_STRESS -// Copyright (C) 2019-2023 Lawrence Livermore National Security, LLC., Xavier Andrade, Alfredo A. Correa +// Copyright (C) 2019-2024 Lawrence Livermore National Security, LLC., Xavier Andrade, Alfredo A. Correa // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -39,16 +39,17 @@ struct forces_stress { template forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham): + forces(ions.size()), stress({3, 3}, 0.0) { - forces = calculate_forces(ions, electrons, ham); + calculate(ions, electrons, ham); } private: template - gpu::array, 1> calculate_forces(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){ - + void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){ + CALI_CXX_MARK_FUNCTION; basis::field> gdensity(electrons.density_basis()); @@ -100,14 +101,12 @@ struct forces_stress { } } - gpu::array, 1> forces(ions.size()); + for(int iatom = 0; iatom < ions.size(); iatom++){ forces[iatom] = ionic_forces[iatom] + forces_local[iatom] + forces_non_local[iatom]; } // MISSING: the non-linear core correction term - - return forces; } }; @@ -116,8 +115,8 @@ struct forces_stress { } #endif -#ifdef INQ_OBSERVABLES_FORCES_UNIT_TEST -#undef INQ_OBSERVABLES_FORCES_UNIT_TEST +#ifdef INQ_OBSERVABLES_FORCES_STRESS_UNIT_TEST +#undef INQ_OBSERVABLES_FORCES_STRESS_UNIT_TEST #include #include diff --git a/src/real_time/propagate.hpp b/src/real_time/propagate.hpp index 8010506a..638a801f 100644 --- a/src/real_time/propagate.hpp +++ b/src/real_time/propagate.hpp @@ -12,7 +12,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/src/real_time/viewables.hpp b/src/real_time/viewables.hpp index ef6fe314..e20b5d0e 100644 --- a/src/real_time/viewables.hpp +++ b/src/real_time/viewables.hpp @@ -13,7 +13,7 @@ #include #include #include -#include +#include #include #include #include