Skip to content

Commit

Permalink
Unified RKS/UKS unittest integrator
Browse files Browse the repository at this point in the history
  • Loading branch information
@elambros
elambros committed Aug 9, 2023
1 parent bfe96ae commit 47a8256
Showing 1 changed file with 68 additions and 151 deletions.
219 changes: 68 additions & 151 deletions tests/xc_integrator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@

using namespace GauXC;


void test_xc_integrator( ExecutionSpace ex, const RuntimeEnvironment& rt,
std::string reference_file,
ExchCXX::Functional func_key,
Expand All @@ -35,25 +36,46 @@ void test_xc_integrator( ExecutionSpace ex, const RuntimeEnvironment& rt,
using matrix_type = Eigen::MatrixXd;
Molecule mol;
BasisSet<double> basis;
matrix_type P, VXC_ref, K_ref;
matrix_type P, Pz, VXC_ref, VXCz_ref, K_ref;
double EXC_ref;
std::vector<double> EXC_GRAD_ref;
bool has_k = false, has_exc_grad = false;
bool has_k = false, has_exc_grad = false, uks = false;
{
read_hdf5_record( mol, reference_file, "/MOLECULE" );
read_hdf5_record( basis, reference_file, "/BASIS" );

HighFive::File file( reference_file, HighFive::File::ReadOnly );
auto dset = file.getDataSet("/DENSITY");

std::string den="/DENSITY";
std::string den2="/DENSITY_Z";
std::string vxc="/VXC";
std::string vxc2="VXC_Z";

if (file.exist("/DENSITY_Z")) {
den="/DENSITY_SCALAR";
vxc="/VXC_SCALAR";
uks=true;
}

auto dset = file.getDataSet(den);

auto dims = dset.getDimensions();
P = matrix_type( dims[0], dims[1] );
VXC_ref = matrix_type( dims[0], dims[1] );
P = matrix_type( dims[0], dims[1] );
VXC_ref = matrix_type( dims[0], dims[1] );
Pz = matrix_type( dims[0], dims[1] );
VXCz_ref = matrix_type( dims[0], dims[1] );

dset.read( P.data() );
dset = file.getDataSet("/VXC");
dset = file.getDataSet(vxc);
dset.read( VXC_ref.data() );

if (uks) {
dset = file.getDataSet(den2);
dset.read( Pz.data() );
dset = file.getDataSet(vxc2);
dset.read( VXCz_ref.data() );
}

dset = file.getDataSet("/EXC");
dset.read( &EXC_ref );

Expand Down Expand Up @@ -91,40 +113,59 @@ void test_xc_integrator( ExecutionSpace ex, const RuntimeEnvironment& rt,
mw.modify_weights(lb);

// Construct XC Functional
functional_type func( ExchCXX::Backend::builtin, func_key, ExchCXX::Spin::Unpolarized );
ExchCXX::Spin Spin=ExchCXX::Spin::Unpolarized; if (uks) { Spin=ExchCXX::Spin::Polarized; }
functional_type func( ExchCXX::Backend::builtin, func_key, Spin );

// Construct XCIntegrator
XCIntegratorFactory<matrix_type> integrator_factory( ex, "Replicated",
integrator_kernel, lwd_kernel, reduction_kernel );
auto integrator = integrator_factory.get_instance( func, lb );

// Integrate Density
if( check_integrate_den ) {
if( check_integrate_den and not uks ) {
auto N_EL_ref = std::accumulate( mol.begin(), mol.end(), 0ul,
[](const auto& a, const auto &b) { return a + b.Z.get(); });
auto N_EL = integrator.integrate_den( P );
CHECK( N_EL == Approx(N_EL_ref).epsilon(1e-6) );
}

// Integrate EXC/VXC
auto [ EXC, VXC ] = integrator.eval_exc_vxc( P );

// Check EXC/VXC
auto VXC_diff_nrm = ( VXC - VXC_ref ).norm();
CHECK( EXC == Approx( EXC_ref ) );
CHECK( VXC_diff_nrm / basis.nbf() < 1e-10 );
if ( not uks ) {
auto [ EXC, VXC ] = integrator.eval_exc_vxc( P );

// Check EXC/VXC
auto VXC_diff_nrm = ( VXC - VXC_ref ).norm();
CHECK( EXC == Approx( EXC_ref ) );
CHECK( VXC_diff_nrm / basis.nbf() < 1e-10 );
// Check if the integrator propagates state correctly
{
auto [ EXC1, VXC1 ] = integrator.eval_exc_vxc( P );
CHECK( EXC1 == Approx( EXC_ref ) );
auto VXC1_diff_nrm = ( VXC1 - VXC_ref ).norm();
CHECK( VXC1_diff_nrm / basis.nbf() < 1e-10 );
}
} else {
auto [ EXC, VXC, VXCz ] = integrator.eval_exc_vxc( P, Pz );

// Check EXC/VXC
auto VXC_diff_nrm = ( VXC - VXC_ref ).norm();
auto VXCz_diff_nrm = ( VXCz - VXCz_ref ).norm();
CHECK( EXC == Approx( EXC_ref ) );
CHECK( VXC_diff_nrm / basis.nbf() < 1e-10 );
CHECK( VXCz_diff_nrm / basis.nbf() < 1e-10 );
// Check if the integrator propagates state correctly
{
auto [ EXC1, VXC1, VXCz1 ] = integrator.eval_exc_vxc( P, Pz );
CHECK( EXC1 == Approx( EXC_ref ) );
auto VXC1_diff_nrm = ( VXC1 - VXC_ref ).norm();
auto VXCz1_diff_nrm = ( VXCz1 - VXCz_ref ).norm();
CHECK( VXC1_diff_nrm / basis.nbf() < 1e-10 );
CHECK( VXC1_diff_nrm / basis.nbf() < 1e-10 );
}

// Check if the integrator propagates state correctly
{
auto [ EXC1, VXC1 ] = integrator.eval_exc_vxc( P );
CHECK( EXC1 == Approx( EXC_ref ) );
auto VXC1_diff_nrm = ( VXC1 - VXC_ref ).norm();
CHECK( VXC1_diff_nrm / basis.nbf() < 1e-10 );
}


// Check EXC Grad
if( check_grad and has_exc_grad ) {
if( check_grad and has_exc_grad and not uks) {
auto EXC_GRAD = integrator.eval_exc_grad( P );
using map_type = Eigen::Map<Eigen::MatrixXd>;
map_type EXC_GRAD_ref_map( EXC_GRAD_ref.data(), mol.size(), 3 );
Expand All @@ -134,131 +175,15 @@ void test_xc_integrator( ExecutionSpace ex, const RuntimeEnvironment& rt,
}

// Check K
if( has_k and check_k ) {
if( has_k and check_k and not uks ) {
auto K = integrator.eval_exx( P );
CHECK((K - K.transpose()).norm() < std::numeric_limits<double>::epsilon()); // Symmetric
CHECK( (K - K_ref).norm() / basis.nbf() < 1e-7 );
}

}

void test_xc_integrator_uks( ExecutionSpace ex, const RuntimeEnvironment& rt,
std::string reference_file,
ExchCXX::Functional func_key,
PruningScheme pruning_scheme,
size_t quad_pad_value,
bool check_grad,
bool check_integrate_den,
bool check_k,
std::string integrator_kernel = "Default",
std::string reduction_kernel = "Default",
std::string lwd_kernel = "Default") {

// Read the reference file
using matrix_type = Eigen::MatrixXd;
Molecule mol;
BasisSet<double> basis;
matrix_type Pscalar, Pz, VXCscalar_ref, VXCz_ref, K_ref;
double EXC_ref;
std::vector<double> EXC_GRAD_ref;
bool has_k = false, has_exc_grad = false;
{
read_hdf5_record( mol, reference_file, "/MOLECULE" );
read_hdf5_record( basis, reference_file, "/BASIS" );

HighFive::File file( reference_file, HighFive::File::ReadOnly );
auto dset = file.getDataSet("/DENSITY_SCALAR");

auto dims = dset.getDimensions();
Pscalar = matrix_type( dims[0], dims[1] );
VXCscalar_ref = matrix_type( dims[0], dims[1] );
Pz = matrix_type( dims[0], dims[1] );
VXCz_ref = matrix_type( dims[0], dims[1] );

dset.read( Pscalar.data() );

dset = file.getDataSet("/DENSITY_Z");
dset.read( Pz.data() );

dset = file.getDataSet("/VXC_SCALAR");
dset.read( VXCscalar_ref.data() );

dset = file.getDataSet("/VXC_Z");
dset.read( VXCz_ref.data() );

dset = file.getDataSet("/EXC");
dset.read( &EXC_ref );

has_exc_grad = file.exist("/EXC_GRAD");
if( has_exc_grad ) {

EXC_GRAD_ref.resize( 3*mol.size() );
dset = file.getDataSet("/EXC_GRAD");
dset.read( EXC_GRAD_ref.data() );
}

has_k = file.exist("/K");
if(has_k) {
K_ref = matrix_type(dims[0], dims[1]);
dset = file.getDataSet("/K");
dset.read( K_ref.data() );
}
}

for( auto& sh : basis )
sh.set_shell_tolerance( std::numeric_limits<double>::epsilon() );

auto mg = MolGridFactory::create_default_molgrid(mol, pruning_scheme,
BatchSize(512), RadialQuad::MurrayHandyLaming, AtomicGridSizeDefault::UltraFineGrid);

// Construct Load Balancer
LoadBalancerFactory lb_factory(ExecutionSpace::Host, "Default");
auto lb = lb_factory.get_instance(rt, mol, mg, basis, quad_pad_value);

// Construct Weights Module
MolecularWeightsFactory mw_factory( ex, "Default", MolecularWeightsSettings{} );
auto mw = mw_factory.get_instance();

// Apply partition weights
mw.modify_weights(lb);

// Construct XC Functional
functional_type func( ExchCXX::Backend::builtin, func_key, ExchCXX::Spin::Polarized );

// Construct XCIntegrator
XCIntegratorFactory<matrix_type> integrator_factory( ex, "Replicated",
integrator_kernel, lwd_kernel, reduction_kernel );
auto integrator = integrator_factory.get_instance( func, lb );

// Integrate Density
if( check_integrate_den ) {
GAUXC_GENERIC_EXCEPTION("INTEGRATE_DEN NYI FOR UKS");
}

// Integrate EXC/VXC
auto [ EXC, VXCscalar, VXCz ] = integrator.eval_exc_vxc( Pscalar, Pz );

// Check EXC/VXC
auto VXCscalar_diff_nrm = ( VXCscalar - VXCscalar_ref ).norm();
auto VXCz_diff_nrm = ( VXCz - VXCz_ref ).norm();
CHECK( EXC == Approx( EXC_ref ) );
CHECK( VXCscalar_diff_nrm / basis.nbf() < 1e-10 );
CHECK( VXCz_diff_nrm / basis.nbf() < 1e-10 );
// Check if the integrator propagates state correctly
{
auto [ EXC1, VXCscalar1, VXCz1 ] = integrator.eval_exc_vxc( Pscalar, Pz );
CHECK( EXC1 == Approx( EXC_ref ) );
auto VXCscalar1_diff_nrm = ( VXCscalar1 - VXCscalar_ref ).norm();
auto VXCz1_diff_nrm = ( VXCz1 - VXCz_ref ).norm();
CHECK( VXCscalar1_diff_nrm / basis.nbf() < 1e-10 );
CHECK( VXCz1_diff_nrm / basis.nbf() < 1e-10 );
}


}


void test_integrator(std::string reference_file, ExchCXX::Functional func, PruningScheme pruning_scheme, bool uks=false) {
void test_integrator(std::string reference_file, ExchCXX::Functional func, PruningScheme pruning_scheme) {

#ifdef GAUXC_ENABLE_DEVICE
auto rt = DeviceRuntimeEnvironment(GAUXC_MPI_CODE(MPI_COMM_WORLD,) 0.9);
Expand All @@ -267,18 +192,10 @@ void test_integrator(std::string reference_file, ExchCXX::Functional func, Pruni
#endif

#ifdef GAUXC_ENABLE_HOST

if (uks) {
SECTION( "Host" ) {
test_xc_integrator_uks( ExecutionSpace::Host, rt, reference_file, func,
pruning_scheme, 1, true, false, true );
}
} else {
SECTION( "Host" ) {
test_xc_integrator( ExecutionSpace::Host, rt, reference_file, func,
pruning_scheme, 1, true, true, true );
}
}
#endif

#ifdef GAUXC_ENABLE_DEVICE
Expand Down Expand Up @@ -359,13 +276,13 @@ TEST_CASE( "XC Integrator", "[xc-integrator]" ) {
//UKS LDA Test
SECTION( "Li / SVWN5 / sto-3g" ) {
test_integrator(GAUXC_REF_DATA_PATH "/li_svwn5_sto3g_uks.bin",
ExchCXX::Functional::SVWN5, PruningScheme::Unpruned, true );
ExchCXX::Functional::SVWN5, PruningScheme::Unpruned );
}

//UKS GGA Test
SECTION( "Li / BLYP / sto-3g" ) {
test_integrator(GAUXC_REF_DATA_PATH "/li_blyp_sto3g_uks.bin",
ExchCXX::Functional::BLYP, PruningScheme::Unpruned, true );
ExchCXX::Functional::BLYP, PruningScheme::Unpruned );
}


Expand Down

0 comments on commit 47a8256

Please sign in to comment.