Merge branch 'species_set' into 'master'

Add a species_set class with all the species objects

See merge request npneq/inq!1098

src/hamiltonian/atomic_potential.hpp

@@ -45,20 +45,17 @@ namespace hamiltonian {
 	private:
 
 		pseudo::math::erf_range_separation sep_;
-		int natoms_;
-		double nelectrons_;
 		pseudo::set pseudo_set_;
 		std::unordered_map<std::string, pseudopotential_type> pseudopotential_list_;
 		bool has_nlcc_;
 		basis::double_grid double_grid_;
 		bool fourier_pseudo_;
 
 	public:
-		
-		template <class atom_array>
-		atomic_potential(const int natoms, const atom_array & atom_list, double gcutoff, options::electrons const & conf = {}):
+
+		template <class SpeciesList>
+		atomic_potential(SpeciesList const & species_list, double gcutoff, options::electrons const & conf = {}):
 			sep_(0.625), //this is the default from octopus
-			natoms_(natoms),
 			pseudo_set_(conf.pseudopotentials_value()),
 			double_grid_(conf.double_grid_value()),
 			fourier_pseudo_(conf.fourier_pseudo_value())
@@ -69,29 +66,20 @@ namespace hamiltonian {
 			gcutoff *= double_grid_.spacing_factor(); 
 
 			has_nlcc_ = false;
-			nelectrons_ = 0.0;
 
-			for(int iatom = 0; iatom < natoms; iatom++){
-				if(!pseudo_set_.has(atom_list[iatom])) throw std::runtime_error("inq error: pseudopotential for element " + atom_list[iatom].symbol() + " not found.");
+			for(auto const & species : species_list) {
+				if(!pseudo_set_.has(species)) throw std::runtime_error("inq error: pseudopotential for element " + species.symbol() + " not found.");
 
-				auto map_ref = pseudopotential_list_.find(atom_list[iatom].symbol());
+				if(pseudopotential_list_.find(species.symbol()) != pseudopotential_list_.end()) throw std::runtime_error("INQ Error: duplicated species");
 
-				if(map_ref == pseudopotential_list_.end()){
-
-					auto file_path = pseudo_set_.file_path(atom_list[iatom]);
-					if(atom_list[iatom].has_file()) file_path = atom_list[iatom].file_path();
+				auto file_path = pseudo_set_.file_path(species);
+				if(species.has_file()) file_path = species.file_path();
 
-					//sorry for this, emplace has a super ugly syntax
-					auto insert = pseudopotential_list_.emplace(std::piecewise_construct, std::make_tuple(atom_list[iatom].symbol()), std::make_tuple(file_path, sep_, gcutoff, atom_list[iatom].filter_pseudo()));
-					map_ref = insert.first;
-
-				}
-
-				auto & pseudo = map_ref->second;
+				//sorry for this, emplace has a super ugly syntax
+				auto insert = pseudopotential_list_.emplace(std::piecewise_construct, std::make_tuple(species.symbol()), std::make_tuple(file_path, sep_, gcutoff, species.filter_pseudo()));
+				auto & pseudo = insert.first->second;
 
-				nelectrons_ += pseudo.valence_charge();
 				has_nlcc_ = has_nlcc_ or pseudo.has_nlcc_density();
-
 			}
 
 		}
@@ -100,12 +88,19 @@ namespace hamiltonian {
 			return pseudopotential_list_.size();
 		}
 
-		const double & num_electrons() const {
-			return nelectrons_;
+		template <class SymbolsList>
+		auto num_electrons(SymbolsList const symbols_list) {
+			double nel = 0.0;
+			for(auto symbol : symbols_list) {
+				auto map_ref = pseudopotential_list_.find(symbol);
+				assert(map_ref != pseudopotential_list_.end());
+				auto & pseudo = map_ref->second;
+				nel += pseudo.valence_charge();
+			}
+			return nel;
 		}
-
-		template <class element_type>
-		const pseudopotential_type & pseudo_for_element(const element_type & el) const {
+
+		auto & pseudo_for_element(input::species const & el) const {
 			return pseudopotential_list_.at(el.symbol());
 		}
 
@@ -116,7 +111,7 @@ namespace hamiltonian {
 
 			basis::field<basis_type, double> potential(basis);
 
-			parallel::partition part(natoms_, comm);
+			parallel::partition part(ions.size(), comm);
 
 			potential.fill(0.0);
 
@@ -126,7 +121,7 @@ namespace hamiltonian {
 
 				auto atom_position = ions.positions()[iatom];
 
-				auto & ps = pseudo_for_element(ions.atoms()[iatom]);
+				auto & ps = pseudo_for_element(ions.species(iatom));
 				basis::spherical_grid sphere(basis, atom_position, ps.short_range_potential_radius());
 
 				if(not double_grid_.enabled()){
@@ -167,7 +162,7 @@ namespace hamiltonian {
 
 			CALI_CXX_MARK_FUNCTION;
 
-			parallel::partition part(natoms_, comm);
+			parallel::partition part(ions.size(), comm);
 
 			basis::field<basis_type, double> density(basis);
 
@@ -179,7 +174,7 @@ namespace hamiltonian {
 
 				auto atom_position = ions.positions()[iatom];
 
-				auto & ps = pseudo_for_element(ions.atoms()[iatom]);
+				auto & ps = pseudo_for_element(ions.species(iatom));
 				basis::spherical_grid sphere(basis, atom_position, sep_.long_range_density_radius());
 
 				//OPTIMIZATION: this should be done in parallel for atoms too
@@ -204,7 +199,7 @@ namespace hamiltonian {
 
 			CALI_CXX_MARK_FUNCTION;
 
-			parallel::partition part(natoms_, comm);
+			parallel::partition part(ions.size(), comm);
 
 			auto nspin = states.num_density_components();
 			basis::field_set<basis_type, double> density(basis, nspin);
@@ -218,7 +213,7 @@ namespace hamiltonian {
 
 				auto atom_position = ions.positions()[iatom];
 
-				auto & ps = pseudo_for_element(ions.atoms()[iatom]);
+				auto & ps = pseudo_for_element(ions.species(iatom));
 
 				if(ps.has_electronic_density()){
 
@@ -267,7 +262,7 @@ namespace hamiltonian {
 
 			CALI_CXX_MARK_FUNCTION;
 
-			parallel::partition part(natoms_, comm);
+			parallel::partition part(ions.size(), comm);
 
 			basis::field<basis_type, double> density(basis);
 
@@ -277,7 +272,7 @@ namespace hamiltonian {
 
 				auto atom_position = ions.positions()[iatom];
 
-				auto & ps = pseudo_for_element(ions.atoms()[iatom]);
+				auto & ps = pseudo_for_element(ions.species(iatom));
 
 				if(not ps.has_nlcc_density()) continue;
 
@@ -304,7 +299,6 @@ namespace hamiltonian {
 		void info(output_stream & out) const {
 			out << "ATOMIC POTENTIAL:" << std::endl;
 			out << "	Number of species		= " << num_species() << std::endl;
-			out << "	Number of electrons = " << num_electrons() << std::endl;
 			out << std::endl;
 		}
 
@@ -347,38 +341,35 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){
 	SECTION("Non-existing element"){
 		std::vector<species> el_list({"P", "X"});
 
-		CHECK_THROWS(hamiltonian::atomic_potential(el_list.size(), el_list, gcut));
+		CHECK_THROWS(hamiltonian::atomic_potential(el_list, gcut));
 	}
 
 	SECTION("Duplicated element"){
 		std::vector<species> el_list({"N", "N"});
-
-		hamiltonian::atomic_potential pot(el_list.size(), el_list.begin(), gcut);
-
-		CHECK(pot.num_species() == 1);
-		CHECK(pot.num_electrons() == 10.0_a);
-
+		CHECK_THROWS(hamiltonian::atomic_potential(el_list, gcut));
 	}
 
 	SECTION("Empty list"){
 		std::vector<species> el_list;
 
-		hamiltonian::atomic_potential pot(el_list.size(), el_list, gcut);
+		hamiltonian::atomic_potential pot(el_list, gcut);
 
 		CHECK(pot.num_species() == 0);
-		CHECK(pot.num_electrons() == 0.0_a);
 
 		CHECK(not pot.has_nlcc());
 
 	}
 
 	SECTION("CNOH"){
-		species el_list[] = {"C", "N", "O", "H"};
+		ionic::species_set el_list;
+		el_list.insert("C");
+		el_list.insert("N");
+		el_list.insert("O");
+		el_list.insert("H");
 
-		hamiltonian::atomic_potential pot(4, el_list, gcut);
+		hamiltonian::atomic_potential pot(el_list, gcut);
 
 		CHECK(pot.num_species() == 4);
-		CHECK(pot.num_electrons() == 16.0_a);
 	}
 
 	SECTION("Construct from a geometry"){
@@ -387,10 +378,10 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){
 
 		basis::real_space rs(ions.cell(), /*spacing = */ 0.49672941, comm);
 
-		hamiltonian::atomic_potential pot(ions.size(), ions.atoms(), rs.gcutoff());
+		hamiltonian::atomic_potential pot(ions.species_list(), rs.gcutoff());
 
 		CHECK(pot.num_species() == 2);
-		CHECK(pot.num_electrons() == 30.0_a);
+		CHECK(pot.num_electrons(ions.symbols()) == 30.0_a);
 
 		rs.info(std::cout);
 

src/hamiltonian/ks_hamiltonian.hpp

@@ -62,10 +62,10 @@ class ks_hamiltonian {
 
 		for(int iatom = 0; iatom < ions.size(); iatom++){
 			if(non_local_in_fourier_){
-				auto insert = projectors_fourier_map_.emplace(ions.atoms()[iatom].symbol(), projector_fourier(basis, pot.pseudo_for_element(ions.atoms()[iatom])));
+				auto insert = projectors_fourier_map_.emplace(ions.symbol(iatom), projector_fourier(basis, pot.pseudo_for_element(ions.species(iatom))));
 				insert.first->second.add_coord(basis.cell().metric().to_contravariant(ions.positions()[iatom]));
 			} else {
-				projectors.emplace_back(basis, pot.double_grid(), pot.pseudo_for_element(ions.atoms()[iatom]), ions.positions()[iatom], iatom);
+				projectors.emplace_back(basis, pot.double_grid(), pot.pseudo_for_element(ions.species(iatom)), ions.positions()[iatom], iatom);
 				if(projectors.back().empty()) projectors.pop_back(); 
 			}
 		}
@@ -263,7 +263,7 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG){
 		CHECK(rs.volume_element() == 0.125_a);
 	}
 
-	hamiltonian::atomic_potential pot(ions.size(), ions.atoms(), rs.gcutoff());
+	hamiltonian::atomic_potential pot(ions.species_list(), rs.gcutoff());
 
 	states::ks_states st(states::spin_config::UNPOLARIZED, 11.0);
 

src/ionic/brillouin.hpp

@@ -45,7 +45,7 @@ class brillouin {
 		positions[2] = 0.0;		
 
 		for(int iatom = 0; iatom < ions.size(); iatom++){
-			types[iatom] = ions.atoms()[iatom].atomic_number();
+			types[iatom] = ions.species(iatom).atomic_number();
 			auto pos = ions.cell().metric().to_contravariant(ions.cell().position_in_cell(ions.positions()[iatom]));
 			positions[3*iatom + 0] = pos[0];
 			positions[3*iatom + 1] = pos[1];

src/ionic/interaction.hpp

@@ -326,7 +326,7 @@ auto interaction_energy(const cell_type & cell, const geometry_type & geo, const
 	boost::multi::array<vector3<double>, 1> forces(geo.size());
 	boost::multi::array<double, 1> charges(geo.size());
 
-	for(int ii = 0; ii < geo.size(); ii++) charges[ii] = atomic_pot.pseudo_for_element(geo.atoms()[ii]).valence_charge();
+	for(int ii = 0; ii < geo.size(); ii++) charges[ii] = atomic_pot.pseudo_for_element(geo.species(ii)).valence_charge();
 
 	interaction_energy(geo.size(), cell, charges, geo.positions(), atomic_pot.range_separation(), energy, forces);
 
@@ -344,7 +344,7 @@ auto interaction_forces(const cell_type & cell, const geometry_type & geo, const
 	boost::multi::array<vector3<double>, 1> forces(geo.size());
 	boost::multi::array<double, 1> charges(geo.size());
 
-	for(int ii = 0; ii < geo.size(); ii++) charges[ii] = atomic_pot.pseudo_for_element(geo.atoms()[ii]).valence_charge();
+	for(int ii = 0; ii < geo.size(); ii++) charges[ii] = atomic_pot.pseudo_for_element(geo.species(ii)).valence_charge();
 
 	interaction_energy(geo.size(), cell, charges, geo.positions(), atomic_pot.range_separation(), energy, forces);
 

src/ionic/propagator.hpp

@@ -72,7 +72,7 @@ struct molecular_dynamics{
 	template <typename TypeIons, typename TypeForces>
 	auto acceleration(TypeIons& ions, TypeForces forces) const {
 
-		for(int iatom = 0; iatom < ions.size(); iatom++) forces[iatom] /= ions.atoms()[iatom].mass();
+		for(int iatom = 0; iatom < ions.size(); iatom++) forces[iatom] /= ions.species(iatom).mass();
 		return forces;
 
 	}

src/ionic/species_set.hpp

@@ -0,0 +1,82 @@
+/* -*- indent-tabs-mode: t -*- */
+
+#ifndef INQ__IONIC__SPECIES_SET
+#define INQ__IONIC__SPECIES_SET
+
+// 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 <input/species.hpp>
+
+namespace inq {
+namespace ionic {
+
+class species_set {
+
+	using container_type = 	std::unordered_map<std::string, input::species>;
+
+	container_type list_;
+
+public:
+
+	auto size() const {
+		return (long) list_.size();
+	}
+
+	auto & list() const {
+		return list_;
+	}
+
+	auto insert(input::species const & sp) {
+		list_.insert_or_assign(sp.symbol(), sp);
+	}
+
+	auto & operator[](std::string const & symbol) const{
+		return list_.at(symbol);
+	}
+
+	struct const_iterator {
+
+		container_type::const_iterator base_iter;
+
+		auto operator!=(const_iterator const & other) const {
+			return base_iter != other.base_iter;
+		}
+
+		auto operator++() {
+			return const_iterator{++base_iter};
+		}
+
+		auto operator*() const {
+			return base_iter->second;
+		}
+
+	};
+
+	auto begin() const {
+		return const_iterator{list_.begin()};
+	}
+
+	auto end() const {
+		return const_iterator{list_.end()};
+	}
+
+};
+
+}
+}
+#endif
+
+#ifdef INQ_IONIC_SPECIES_SET_UNIT_TEST
+#undef INQ_IONIC_SPECIES_SET_UNIT_TEST
+
+#include <catch2/catch_all.hpp>
+
+TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {
+
+}
+
+#endif

src/observables/dipole.hpp

@@ -47,7 +47,7 @@ vector3<double> dipole(systems::ions const & ions, const hamiltonian::atomic_pot
 	vector3<double> dip = {0.0, 0.0, 0.0};
 
 	for(int iatom = 0; iatom < ions.size(); iatom++){
-		auto zval = atomic_pot.pseudo_for_element(ions.atoms()[iatom]).valence_charge();
+		auto zval = atomic_pot.pseudo_for_element(ions.species(iatom)).valence_charge();
 		dip += proton_charge*zval*ions.positions()[iatom];
 	}
 

src/systems/electrons.hpp

@@ -135,8 +135,8 @@ class electrons {
 	template <typename KptsType = input::kpoints::list>	
 	electrons(input::parallelization const & dist, const inq::systems::ions & ions, const options::electrons & conf = {}, KptsType const & kpts = input::kpoints::gamma()):
 		brillouin_zone_(ions, kpts),
-		atomic_pot_(ions.size(), ions.atoms(), basis::real_space::gcutoff(ions.cell(), conf.spacing_value()), conf),
-		states_(conf.spin_val(), atomic_pot_.num_electrons() + conf.extra_electrons_val(), conf.extra_states_val(), conf.temperature_val(), kpts.size()),
+		atomic_pot_(ions.species_list(), basis::real_space::gcutoff(ions.cell(), conf.spacing_value()), conf),
+		states_(conf.spin_val(), atomic_pot_.num_electrons(ions.symbols()) + conf.extra_electrons_val(), conf.extra_states_val(), conf.temperature_val(), kpts.size()),
 		full_comm_(dist.cart_comm(conf.num_spin_components_val(), brillouin_zone_.size(), states_.num_states())),
 		kpin_comm_(kpin_subcomm(full_comm_)),
 		kpin_states_comm_(kpin_states_subcomm(full_comm_)),
@@ -186,8 +186,8 @@ class electrons {
 		eigenvalues_.reextent({static_cast<boost::multi::size_t>(kpin_.size()), max_local_spinor_set_size_});
 		occupations_.reextent({static_cast<boost::multi::size_t>(kpin_.size()), max_local_spinor_set_size_});
 
-		if(atomic_pot_.num_electrons() + conf.extra_electrons_val() == 0) throw std::runtime_error("inq error: the system does not have any electrons");
-		if(atomic_pot_.num_electrons() + conf.extra_electrons_val() < 0) throw std::runtime_error("inq error: the system has a negative number of electrons");		
+		if(states_.num_electrons() == 0.0) throw std::runtime_error("inq error: the system does not have any electrons");
+		if(states_.num_electrons() <  0.0) throw std::runtime_error("inq error: the system has a negative number of electrons");
 
 		print(ions);