Merge branch 'force_reorganization' into 'master'

Create the forces_stress object

See merge request npneq/inq!1138

src/ground_state/calculator.hpp

@@ -18,14 +18,14 @@
 #include <hamiltonian/ks_hamiltonian.hpp>
 #include <hamiltonian/self_consistency.hpp>
 #include <hamiltonian/energy.hpp>
-#include <hamiltonian/forces.hpp>
 #include <basis/field_set.hpp>
 #include <operations/randomize.hpp>
 #include <operations/overlap.hpp>
 #include <operations/orthogonalize.hpp>
 #include <operations/preconditioner.hpp>
 #include <operations/integral.hpp>
 #include <observables/density.hpp>
+#include <observables/forces_stress.hpp>
 #include <parallel/gather.hpp>
 #include <mixers/linear.hpp>
 #include <mixers/broyden.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 = hamiltonian::calculate_forces(ions_, electrons, ham_);
+			res.forces = observables::forces_stress{ions_, electrons, ham_}.forces;
 		}
 
 		if(solver_.calc_forces() and electrons.states().spinor_dim() == 2) {

src/observables/forces_stress.hpp

@@ -0,0 +1,131 @@
+/* -*- indent-tabs-mode: t -*- */
+
+#ifndef INQ__OBSERVABLES__FORCES_STRESS
+#define INQ__OBSERVABLES__FORCES_STRESS
+
+// 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
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include <observables/density.hpp>
+#include <operations/gradient.hpp>
+#include <solvers/poisson.hpp>
+#include <systems/ions.hpp>
+#include <systems/electrons.hpp>
+#include <utils/raw_pointer_cast.hpp>
+
+namespace inq {
+namespace observables {
+
+template <typename LongRangeType, typename ShortRangeType, typename GDensityType>
+struct loc_pot {
+
+	LongRangeType v1;
+	ShortRangeType v2;
+	GDensityType gdensityp;
+
+	GPU_FUNCTION auto operator()(long ip) const {
+		return (v1[ip] + v2[ip])*gdensityp[ip];
+	}
+};
+
+struct forces_stress {
+	gpu::array<vector3<double>, 1> forces;
+	gpu::array<double, 2>          stress;
+
+	forces_stress() = default;
+
+	template <typename HamiltonianType>
+	forces_stress(systems::ions const & ions, systems::electrons const & electrons, HamiltonianType const & ham):
+		forces(ions.size()),
+		stress({3, 3}, 0.0)
+	{
+		calculate(ions, electrons, ham);
+	}
+
+private:
+
+	template <typename HamiltonianType>
+	void calculate(const systems::ions & ions, systems::electrons const & electrons, HamiltonianType const & ham){
+
+		CALI_CXX_MARK_FUNCTION;
+
+		basis::field<basis::real_space, vector3<double, covariant>> gdensity(electrons.density_basis());
+		gdensity.fill(vector3<double, covariant>{0.0, 0.0, 0.0});
+
+		gpu::array<vector3<double>, 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++;
+		}
+
+		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<vector3<double>, 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<decltype(begin(ionic_long_range.linear())), decltype(begin(ionic_short_range.linear())), decltype(begin(gdensity.linear()))>
+																	 {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<double *>(raw_pointer_cast(forces_local.data_elements())), 3*forces_local.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
+	}
+
+};
+
+}
+}
+#endif
+
+#ifdef INQ_OBSERVABLES_FORCES_STRESS_UNIT_TEST
+#undef INQ_OBSERVABLES_FORCES_STRESS_UNIT_TEST
+
+#include <catch2/catch_all.hpp>
+#include <basis/real_space.hpp>
+
+TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){
+
+	using namespace inq;
+	using namespace Catch::literals;
+
+}
+#endif
+

src/real_time/propagate.hpp

@@ -10,9 +10,9 @@
 
 #include <systems/ions.hpp>
 #include <hamiltonian/self_consistency.hpp>
-#include <hamiltonian/forces.hpp>
 #include <operations/overlap_diagonal.hpp>
 #include <observables/dipole.hpp>
+#include <observables/forces_stress.hpp>
 #include <options/real_time.hpp>
 #include <perturbations/none.hpp>
 #include <ionic/propagator.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(hamiltonian::calculate_forces(ions, electrons, ham)){};
-		if(ion_propagator.needs_force()) forces = hamiltonian::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<double, covariant>{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 = hamiltonian::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);
 

src/real_time/viewables.hpp

@@ -10,10 +10,10 @@
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
 #include <systems/ions.hpp>
-#include <hamiltonian/forces.hpp>
 #include <operations/overlap_diagonal.hpp>
-#include <observables/dipole.hpp>
 #include <observables/current.hpp>
+#include <observables/dipole.hpp>
+#include <observables/forces_stress.hpp>
 #include <perturbations/none.hpp>
 #include <systems/electrons.hpp>
 #include <real_time/crank_nicolson.hpp>
@@ -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::forces_stress{ions_, electrons_, ham_}.forces;
 		return forces_;
 	}