Merge branch 'relativistic_projector' into 'master'

Relativistic projectors

See merge request npneq/inq!1172

Author: xavierandrade

external_libs/pseudopod

@@ -16,6 +16,7 @@
 #include <hamiltonian/projector.hpp>
 #include <hamiltonian/projector_all.hpp>
 #include <hamiltonian/projector_fourier.hpp>
+#include <hamiltonian/relativistic_projector.hpp>
 #include <hamiltonian/scalar_potential.hpp>
 #include <input/environment.hpp>
 #include <operations/transform.hpp>
@@ -48,6 +49,7 @@ class ks_hamiltonian {
 	bool non_local_in_fourier_;
 	std::unordered_map<std::string, projector_fourier> projectors_fourier_map_;
 	std::vector<std::unordered_map<std::string, projector_fourier>::iterator> projectors_fourier_;
+	std::list<relativistic_projector> projectors_rel_;
 	states::ks_states states_;
 
 #ifdef ENABLE_CUDA
@@ -72,15 +74,20 @@ class ks_hamiltonian {
 
 		std::list<projector> projectors;
 
-		projectors_fourier_map_.clear();			
+		projectors_fourier_map_.clear();
 
 		for(int iatom = 0; iatom < ions.size(); iatom++){
-			if(non_local_in_fourier_){
-				auto insert = projectors_fourier_map_.emplace(ions.symbol(iatom), projector_fourier(basis, pot.pseudo_for_element(ions.species(iatom))));
+			auto && ps = pot.pseudo_for_element(ions.species(iatom));
+
+			if(ps.has_total_angular_momentum()){
+				projectors_rel_.emplace_back(basis, pot.double_grid(), ps, ions.positions()[iatom], iatom);
+				if(projectors_rel_.back().empty()) projectors_rel_.pop_back();
+			} else if(non_local_in_fourier_){
+				auto insert = projectors_fourier_map_.emplace(ions.symbol(iatom), projector_fourier(basis, ps));
 				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.species(iatom)), ions.positions()[iatom], iatom);
-				if(projectors.back().empty()) projectors.pop_back(); 
+				projectors.emplace_back(basis, pot.double_grid(), ps, ions.positions()[iatom], iatom);
+				if(projectors.back().empty()) projectors.pop_back();
 			}
 		}
 
@@ -136,6 +143,8 @@ class ks_hamiltonian {
 			vnlphi.fill(0.0);
 
 			projectors_all_.apply(proj, vnlphi, phi.kpoint() + uniform_vector_potential_);
+
+			for(auto & pr : projectors_rel_) pr.apply(phi, vnlphi, phi.kpoint() + uniform_vector_potential_);
 
 			return vnlphi;
 		}
@@ -156,7 +165,9 @@ class ks_hamiltonian {
 			return en;
 
 		} else {
-			return projectors_all_.energy(phi, phi.kpoint() + uniform_vector_potential_, occupations, reduce_states);
+			auto en = projectors_all_.energy(phi, phi.kpoint() + uniform_vector_potential_, occupations, reduce_states);
+			for(auto & pr : projectors_rel_) en += pr.energy(phi, occupations, phi.kpoint() + uniform_vector_potential_);
+			return en;
 		}
 
 	}
@@ -180,6 +191,7 @@ class ks_hamiltonian {
 		hamiltonian::scalar_potential_add(scalar_potential_, phi.spin_index(), 0.5*phi.basis().cell().metric().norm(phi.kpoint() + uniform_vector_potential_), phi, hphi);
 		exchange_(phi, hphi);
 
+		for(auto & pr : projectors_rel_) pr.apply(phi, hphi, phi.kpoint() + uniform_vector_potential_);
 		projectors_all_.apply(proj, hphi, phi.kpoint() + uniform_vector_potential_);
 
 		return hphi;
@@ -198,7 +210,8 @@ class ks_hamiltonian {
 		auto hphi_rs = hamiltonian::scalar_potential(scalar_potential_, phi.spin_index(), 0.5*phi.basis().cell().metric().norm(phi.kpoint() + uniform_vector_potential_), phi_rs);
 
 		exchange_(phi_rs, hphi_rs);
- 
+
+		for(auto & pr : projectors_rel_) pr.apply(phi_rs, hphi_rs, phi.kpoint() + uniform_vector_potential_);
 		projectors_all_.apply(proj, hphi_rs, phi.kpoint() + uniform_vector_potential_);
 
 		auto hphi = operations::transform::to_fourier(hphi_rs);

share/unit_tests_data/He_fr.upf.gz

@@ -90,16 +90,6 @@ class projector {
 		return nproj_ == 0 or sphere_.size() == 0;
 	}
 
-	template <typename OcType, typename PhiType, typename GPhiType>
-	struct force_term {
-		OcType oc;
-		PhiType phi;
-		GPhiType gphi;
-		constexpr auto operator()(int ist, int ip) const {
-			return -2.0*oc[ist]*real(phi[ip][ist]*conj(gphi[ip][ist]));
-		}
-	};
-	
 	int num_projectors() const {
 		return nproj_;
 	}

share/unit_tests_data/Xe_fr.UPF.gz

@@ -0,0 +1,275 @@
+/* -*- indent-tabs-mode: t -*- */
+
+#ifndef INQ__HAMILTONIAN__RELATIVISTIC_PROJECTOR
+#define INQ__HAMILTONIAN__RELATIVISTIC_PROJECTOR
+
+// Copyright (C) 2019-2023 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 <pseudopod/math/sharmonic.hpp>
+
+#include <gpu/array.hpp>
+#include <math/vector3.hpp>
+#include <basis/double_grid.hpp>
+#include <basis/real_space.hpp>
+#include <basis/spherical_grid.hpp>
+#include <hamiltonian/atomic_potential.hpp>
+#include <utils/profiling.hpp>
+#include <utils/raw_pointer_cast.hpp>
+
+namespace inq {
+namespace hamiltonian {
+
+class relativistic_projector {
+
+	basis::spherical_grid sphere_;
+	int nproj_;
+	gpu::array<double, 3> beta_;
+	gpu::array<double, 1> kb_coeff_;
+	int iatom_;
+
+public: // for CUDA
+	
+	void build(basis::real_space const & basis, basis::double_grid const & double_grid, atomic_potential::pseudopotential_type const & ps) {
+
+		CALI_CXX_MARK_SCOPE("relativistic_projector::build");
+
+		assert(ps.has_total_angular_momentum());
+		
+		nproj_ = 0.0;
+		for(int iproj = 0; iproj < ps.num_projectors_l(); iproj++){
+			nproj_ += ps.projector_2j(iproj) + 1;
+		}
+		
+		beta_.reextent({nproj_, sphere_.size(), 2});
+		kb_coeff_.reextent(nproj_);
+		
+		int iproj_lm = 0;
+		for(int iproj = 0; iproj < ps.num_projectors_l(); iproj++){
+				
+			auto ll = ps.projector_l(iproj);
+			int const jj = ps.projector_2j(iproj);
+
+			assert(jj%2 == 1);
+			
+			//			std::cout << "LL = " << ll << " JJ = " << jj/2.0 << std::endl;
+
+			auto sgn = 1.0;
+			if(jj == 2*ll - 1.0) sgn = -1.0;
+
+			for(auto mj = -jj; mj <= jj; mj += 2){
+
+				auto den = sqrt(jj - sgn + 1);
+				auto cc0 = sgn*sqrt(jj/2.0 - 0.5*sgn + sgn*mj/2.0 + 0.5)/den;
+				auto cc1 =     sqrt(jj/2.0 - 0.5*sgn - sgn*mj/2.0 + 0.5)/den;
+				auto ml0 = (mj - 1)/2;
+				auto ml1 = (mj + 1)/2;
+
+				//				std::cout << cc0 << '\t' << cc1 << '\t' << ml0 << '\t' << ml1 << std::endl;												
+				
+				gpu::run(sphere_.size(),
+								 [bet = begin(beta_),
+									spline = ps.projector(iproj).function(),
+									sph = sphere_.ref(), cc0, cc1, ml0, ml1, ll, iproj_lm,
+									metric = basis.cell().metric()] GPU_LAMBDA (auto ipoint) {
+
+									 if(abs(ml0) <= ll) {
+										 bet[iproj_lm][ipoint][0] = cc0*spline(sph.distance(ipoint))*pseudo::math::sharmonic(ll, ml0, metric.to_cartesian(sph.point_pos(ipoint)));
+									 } else {
+										 bet[iproj_lm][ipoint][0] = 0.0;
+									 }
+									 if(abs(ml1) <= ll) {									 
+										 bet[iproj_lm][ipoint][1] = cc1*spline(sph.distance(ipoint))*pseudo::math::sharmonic(ll, ml1, metric.to_cartesian(sph.point_pos(ipoint)));
+									 } else {
+										 bet[iproj_lm][ipoint][1] = 0.0;
+									 }
+								 });
+
+				kb_coeff_[iproj_lm] = ps.kb_coeff(iproj);
+				
+				iproj_lm++;
+			}
+			
+		}
+	}
+	
+public:
+	relativistic_projector(const basis::real_space & basis, basis::double_grid const & double_grid, atomic_potential::pseudopotential_type const & ps, vector3<double> atom_position, int iatom):
+		sphere_(basis, atom_position, ps.projector_radius()),
+		iatom_(iatom){
+
+		build(basis, double_grid, ps);
+	}
+
+	relativistic_projector(relativistic_projector const &) = delete;		
+
+	auto empty() const {
+		return nproj_ == 0;
+	}
+	
+	auto locally_empty() const {
+		return nproj_ == 0 or sphere_.size() == 0;
+	}
+
+	int num_projectors() const {
+		return nproj_;
+	}
+		
+	auto kb_coeff(int iproj){
+		return kb_coeff_[iproj];
+	}
+
+	auto iatom() const {
+		return iatom_;
+	}
+	
+	auto & sphere() const {
+		return sphere_;
+	}
+
+	auto & beta() const {
+		return beta_;
+	}
+
+	template <typename KpointType>
+	gpu::array<complex, 3> project(states::orbital_set<basis::real_space, complex> const & phi, KpointType const & kpoint) const {
+
+		gpu::array<complex, 3> sphere_phi({sphere_.size(), phi.local_spinor_set_size(), 2});
+
+		gpu::run(phi.local_spinor_set_size(), sphere_.size(),
+						 [gr = begin(phi.spinor_hypercubic()), sph = sphere_.ref(), sgr = begin(sphere_phi)] GPU_LAMBDA (auto ist, auto ipoint){
+							 auto point = sph.grid_point(ipoint);
+							 sgr[ipoint][ist][0] = gr[point[0]][point[1]][point[2]][0][ist];
+							 sgr[ipoint][ist][1] = gr[point[0]][point[1]][point[2]][1][ist];							 
+						 });
+
+		gpu::array<complex, 3> projections({nproj_, phi.local_spinor_set_size(), 2});
+
+		gpu::run(phi.local_spinor_set_size(), nproj_,
+						 [proj = begin(projections), sgr = begin(sphere_phi), bet = begin(beta_), np = sphere_.size(), vol = phi.basis().volume_element()] GPU_LAMBDA (auto ist, auto iproj) {
+							 proj[iproj][ist][0] = 0.0;
+							 proj[iproj][ist][1] = 0.0;
+							 for(int ip = 0; ip < np; ip++) {
+								 proj[iproj][ist][0] += bet[iproj][ip][0]*sgr[ip][ist][0];
+								 proj[iproj][ist][1] += bet[iproj][ip][1]*sgr[ip][ist][1];
+							 }
+							 proj[iproj][ist][0] *= vol;
+							 proj[iproj][ist][1] *= vol;
+						 });
+
+		//missing parallel reduction of projections
+
+		return projections;
+	}
+
+	template <typename KpointType>
+	void apply(states::orbital_set<basis::real_space, complex> const & phi, states::orbital_set<basis::real_space, complex> & vnlphi, KpointType const & kpoint) const {
+
+		auto projections = project(phi, kpoint);
+
+		gpu::run(phi.local_spinor_set_size(), nproj_,
+						 [proj = begin(projections), coe = begin(kb_coeff_)] GPU_LAMBDA (auto ist, auto iproj) {
+
+							 //							 std::cout << "PROJ " <<  iproj << '\t' << proj[iproj][ist][0] << '\t' << proj[iproj][ist][1] << '\t' <<  coe[iproj] << std::endl;
+							 
+							 proj[iproj][ist][0] *= coe[iproj];
+							 proj[iproj][ist][1] *= coe[iproj];
+
+						 });
+
+		gpu::run(phi.local_spinor_set_size(), sphere_.size(),
+						 [gr = begin(vnlphi.spinor_hypercubic()), sph = sphere_.ref(), nproj = nproj_, bet = begin(beta_), proj = begin(projections)] GPU_LAMBDA (auto ist, auto ip){
+							 auto point = sph.grid_point(ip);
+							 for(int iproj = 0; iproj < nproj; iproj++) {
+								 gr[point[0]][point[1]][point[2]][0][ist] += conj(bet[iproj][ip][0])*(proj[iproj][ist][0] + proj[iproj][ist][1]);
+								 gr[point[0]][point[1]][point[2]][1][ist] += conj(bet[iproj][ip][1])*(proj[iproj][ist][0] + proj[iproj][ist][1]);
+							 }
+						 });
+	}
+
+	template <typename Occupations, typename KpointType>
+	double energy(states::orbital_set<basis::real_space, complex> const & phi, Occupations const & occupations, KpointType const & kpoint) const {
+		auto projections = project(phi, kpoint);
+		
+		return gpu::run(gpu::reduce(phi.local_spinor_set_size()), gpu::reduce(nproj_), 0.0,
+						 [proj = begin(projections), coe = begin(kb_coeff_), occ = begin(occupations)] GPU_LAMBDA (auto ist, auto iproj) {
+							 auto pp = proj[iproj][ist][0] + proj[iproj][ist][1];
+							 return real(conj(pp)*pp)*coe[iproj]*occ[ist];
+						 });
+	}
+	
+
+	friend class projector_all;
+	
+    
+};
+  
+}
+}
+#endif
+
+#ifdef INQ_HAMILTONIAN_RELATIVISTIC_PROJECTOR_UNIT_TEST
+#undef INQ_HAMILTONIAN_RELATIVISTIC_PROJECTOR_UNIT_TEST
+
+#include <config/path.hpp>
+
+#include <catch2/catch_all.hpp>
+
+TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {
+  
+	using namespace inq;
+	using namespace inq::magnitude;
+	using namespace Catch::literals;
+	
+	pseudo::math::erf_range_separation const sep(0.625);
+
+	parallel::communicator comm{boost::mpi3::environment::get_world_instance()};
+
+	basis::real_space rs(systems::cell::cubic(10.0_b), /*spacing = */ 0.49672941, comm);
+	basis::double_grid dg(false);
+
+	SECTION("He") {
+	
+		hamiltonian::atomic_potential::pseudopotential_type ps(config::path::unit_tests_data() + "He_fr.upf.gz", sep, rs.gcutoff());
+		
+		hamiltonian::relativistic_projector proj(rs, dg, ps, vector3<double>(0.0, 0.0, 0.0), 77);
+		
+		CHECK(proj.num_projectors() == 10);
+		/*
+			if(not proj.empty()){
+			CHECK(proj.kb_coeff(0) ==  7.494508815_a);
+			CHECK(proj.kb_coeff(1) ==  0.6363049519_a);
+			CHECK(proj.kb_coeff(2) == -4.2939052122_a);
+			CHECK(proj.kb_coeff(3) == -4.2939052122_a);
+			CHECK(proj.kb_coeff(4) == -4.2939052122_a);
+			CHECK(proj.kb_coeff(5) == -1.0069878791_a);
+			CHECK(proj.kb_coeff(6) == -1.0069878791_a);
+			CHECK(proj.kb_coeff(7) == -1.0069878791_a);
+			}
+			
+			CHECK(proj.iatom() == 77);
+		*/
+
+	}
+
+
+	SECTION("Xe") {
+			
+		hamiltonian::atomic_potential::pseudopotential_type ps(config::path::unit_tests_data() + "Xe_fr.UPF.gz", sep, rs.gcutoff());
+		
+		hamiltonian::relativistic_projector proj(rs, dg, ps, vector3<double>(0.0, 0.0, 0.0), 77);
+		
+		CHECK(proj.num_projectors() == 16);
+
+	}
+
+
+	
+}
+
+
+#endif
+