Adding Libxc support

.github/mrchem-codecov.yml

@@ -10,3 +10,4 @@ dependencies:
   - ninja
   - nlohmann_json
   - xcfun
+  - libxc

.github/mrchem-gha.yml

@@ -9,3 +9,4 @@ dependencies:
   - ninja
   - nlohmann_json
   - xcfun
+  - libxc

cmake/custom/main.cmake

@@ -56,6 +56,7 @@ set(_build_type ${CMAKE_BUILD_TYPE})
 # Order IS important here!
 include(${PROJECT_SOURCE_DIR}/external/upstream/fetch_nlohmann_json.cmake)
 include(${PROJECT_SOURCE_DIR}/external/upstream/fetch_xcfun.cmake)
+include(${PROJECT_SOURCE_DIR}/external/upstream/fetch_libxc.cmake)
 include(${PROJECT_SOURCE_DIR}/external/upstream/fetch_eigen3.cmake)
 include(${PROJECT_SOURCE_DIR}/external/upstream/fetch_mrcpp.cmake)
 # reset CMAKE_BUILD_TYPE to whatever it was for MRChem

doc/users/user_inp.rst

@@ -239,7 +239,7 @@ Here we specify the exchange-correlation functional used in DFT
 
     DFT {
       spin = false                          # Use spin-polarized functionals
-      density_cutoff = 0.0                  # Cutoff to set XC potential to zero
+      density_cutoff = 1e-11                # Cutoff to set XC potential to zero
     $functionals
     <func1>     1.0                         # Functional name and coefficient
     <func2>     1.0                         # Functional name and coefficient

external/upstream/fetch_eigen3.cmake

@@ -1,4 +1,4 @@
-find_package(Eigen3 CONFIG QUIET
+find_package(Eigen3 3.4 CONFIG QUIET
   NO_CMAKE_PATH
   NO_CMAKE_PACKAGE_REGISTRY
   )

external/upstream/fetch_libxc.cmake

@@ -0,0 +1,36 @@
+find_package(Libxc QUIET CONFIG)
+
+if(TARGET Libxc::xc)
+  get_target_property(_loc Libxc::xc LOCATION)
+  message(STATUS "Libxc::xc: ${_loc} (found version ${Libxc_VERSION})")
+else()
+  message(STATUS "Suitable LibXC could not be located. Fetching and building!")
+  include(FetchContent)
+  FetchContent_Declare(libxc_sources
+    GIT_REPOSITORY https://gitlab.com/libxc/libxc
+    GIT_TAG 7.0.0
+  )
+  set(CMAKE_INSTALL_INCLUDEDIR "include/" CACHE STRING "" FORCE) # Creates Libxc subdir in install
+  set(BUILD_TESTING     OFF CACHE BOOL "Build LibXC tests" FORCE)
+  set(ENABLE_TESTS      OFF CACHE BOOL "Enable LibXC tests" FORCE)
+  set(BUILD_SHARED_LIBS ON  CACHE BOOL "Build LibXC shared libs" FORCE)
+  set(ENABLE_FORTRAN    OFF CACHE BOOL "Build LibXC Fortran bindings" FORCE)
+  set(DISABLE_FXC       ON  CACHE BOOL "Disable 2nd derivatives (Fxc)" FORCE)
+  set(DISABLE_KXC       ON  CACHE BOOL "Disable 3rd derivatives (Kxc)" FORCE)
+  set(DISABLE_LXC       ON  CACHE BOOL "Disable 4th derivatives (Lxc)" FORCE)
+
+  FetchContent_MakeAvailable(libxc_sources)
+
+  if(TARGET xc)
+    if(NOT TARGET Libxc::xc)
+      add_library(Libxc::xc ALIAS xc)
+    endif()
+
+    target_include_directories(xc
+      INTERFACE
+        $<BUILD_INTERFACE:${libxc_sources_SOURCE_DIR}/src>
+        $<BUILD_INTERFACE:${libxc_sources_BINARY_DIR}>
+        $<INSTALL_INTERFACE:include/Libxc> 
+    )
+  endif()
+endif()

python/mrchem/helpers.py

@@ -122,6 +122,10 @@ def write_scf_fock(user_dict, wf_dict, origin):
             },
         }
 
+    # Exchange-Correlation library
+    if wf_dict["method_type"] in ["dft"]:
+        fock_dict["xc_library"] = user_dict["DFT"]["xc_library"],
+
     # External electric field
     if len(user_dict["ExternalFields"]["electric_field"]) > 0:
         fock_dict["external_operator"] = {
@@ -249,6 +253,8 @@ def write_scf_guess(user_dict, wf_dict):
         "file_CUBE_p": f"{vector_dir}CUBE_p_vector.json",
         "file_CUBE_a": f"{vector_dir}CUBE_a_vector.json",
         "file_CUBE_b": f"{vector_dir}CUBE_b_vector.json",
+        "xc_library": user_dict["DFT"]["xc_library"],
+        "cutoff": user_dict["DFT"]["density_cutoff"],
     }
     return guess_dict
 
@@ -279,6 +285,7 @@ def write_scf_solver(user_dict, wf_dict):
         "energy_thrs": scf_dict["energy_thrs"],
         "orbital_thrs": scf_dict["orbital_thrs"],
         "helmholtz_prec": user_dict["Precisions"]["helmholtz_prec"],
+        "xc_library": user_dict["DFT"]["xc_library"],
         "deltascf_method": scf_dict["deltascf_method"],
     }
 
@@ -466,6 +473,12 @@ def write_rsp_fock(user_dict, wf_dict):
             },
         }
 
+    # Exchange-Correlation library
+    if wf_dict["method_type"] in ["dft"]:
+        fock_dict["xc_library"] = {
+            "xc_library": user_dict["DFT"]["xc_library"],
+        }
+
     # Reaction
     if user_dict["WaveFunction"]["environment"].lower() != "none":
         fock_dict["reaction_operator"] = _reaction_operator_handler(user_dict, rsp=True)
@@ -496,6 +509,8 @@ def write_rsp_solver(user_dict, wf_dict, d):
         "property_thrs": user_dict["Response"]["property_thrs"],
         "helmholtz_prec": user_dict["Precisions"]["helmholtz_prec"],
         "orth_prec": 1.0e-14,
+        "xc_library": user_dict["DFT"]["xc_library"],
+        "cutoff": user_dict["DFT"]["density_cutoff"],
     }
     return solver_dict
 

python/mrchem/input_parser/api.py

@@ -338,15 +338,22 @@ def stencil() -> JSONDict:
                                             'name': 'azora_potential_path',
                                             'type': 'str'}],
                         'name': 'ZORA'},
-                    {   'keywords': [   {   'default': 0.0,
+                    {   'keywords': [   {   'default': 1e-11,
                                             'name': 'density_cutoff',
                                             'type': 'float'},
                                         {   'default': ' ',
                                             'name': 'functionals',
                                             'type': 'str'},
                                         {   'default': "not(user['WaveFunction']['restricted'])",
                                             'name': 'spin',
-                                            'type': 'bool'}],
+                                            'type': 'bool'},
+                                        # {   'default': True,
+                                        #     'name': 'libxc',
+                                        #     'type': 'bool'}
+                                        {   'default': 'xcfun',
+                                            'name': 'xc_library',
+                                            'type': 'str'}
+                                        ],
                         'name': 'DFT'},
                     {   'keywords': [   {   'default': True,
                                             'name': 'dipole_moment',

python/template.yml

@@ -368,6 +368,11 @@ sections:
         default: 0.0
         docstring: |
           Hard cutoff for passing density values to XCFun.
+      - name: xc_library
+        type: str
+        default: xcfun
+        docstring: |
+          Runs Libxc or XCFun
       - name: functionals
         type: str
         default: ' '

setup

@@ -47,6 +47,7 @@ def gen_cmake_command(options, arguments):
     command.append('-DENABLE_MPI={0}'.format(arguments['--mpi']))
     command.append('-DENABLE_OPENMP={0}'.format(arguments['--omp']))
     command.append('-DCMAKE_BUILD_TYPE={0}'.format(arguments['--type']))
+    command.append('-DCMAKE_POLICY_VERSION_MINIMUM=3.5')
     command.append('-G"{0}"'.format(arguments['--generator']))
     if arguments['--cmake-options'] != "''":
         command.append(arguments['--cmake-options'])

src/CMakeLists.txt

@@ -51,6 +51,7 @@ target_link_libraries(mrchem
   PRIVATE
     Eigen3::Eigen
   PUBLIC
+    Libxc::xc
     XCFun::xcfun
     MRCPP::mrcpp
     nlohmann_json::nlohmann_json

src/driver.cpp

@@ -26,6 +26,7 @@
 #include <MRCPP/MWOperators>
 #include <MRCPP/Printer>
 #include <MRCPP/Timer>
+#include <mrdft/MRDFT.h> // libxc debug
 
 #include "driver.h"
 #include <filesystem>
@@ -248,6 +249,15 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
     json json_out = {{"success", true}};
     if (json_scf.contains("properties")) driver::init_properties(json_scf["properties"], mol);
 
+    ///////////////////////////////////////////////////////////
+    //////////////////   Setting XC Library  //////////////////
+    ///////////////////////////////////////////////////////////
+    std::string xc_lib;
+
+    if (json_scf["fock_operator"].contains("xc_library")) {
+        xc_lib = json_scf["fock_operator"]["xc_library"][0].get<std::string>();
+    } else {xc_lib = "xcfun";}
+
     ///////////////////////////////////////////////////////////
     ////////////////   Building Fock Operator   ///////////////
     ///////////////////////////////////////////////////////////
@@ -312,6 +322,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
         solver.setRotation(rotation);
         solver.setLocalize(localize);
         solver.setMethodName(method);
+        solver.setLibxc((xc_lib == "libxc") ? true : false);
         solver.setRelativityName(relativity);
         solver.setEnvironmentName(environment);
         solver.setExternalFieldName(external_field);
@@ -507,14 +518,18 @@ bool driver::scf::guess_energy(const json &json_guess, Molecule &mol, FockBuilde
     auto external_field = json_guess["external_field"];
     auto localize = json_guess["localize"];
     auto rotate = json_guess["rotate"];
+    std::string xc_lib = json_guess["xc_library"].get<std::string>();
+    auto cutoff = json_guess["cutoff"];
 
     mrcpp::print::separator(0, '~');
     print_utils::text(0, "Calculation    ", "Compute initial energy");
     print_utils::text(0, "Method         ", method);
+    print_utils::text(0, "XC Library     ", (xc_lib == "libxc") ? "LibXC" : "XCFun");
     print_utils::text(0, "Relativity     ", relativity);
     print_utils::text(0, "Environment    ", environment);
     print_utils::text(0, "External fields", external_field);
     print_utils::text(0, "Precision      ", print_utils::dbl_to_str(prec, 5, true));
+    print_utils::text(0, "Density cutoff ", print_utils::dbl_to_str(cutoff, 5, true));
     print_utils::text(0, "Localization   ", (localize) ? "On" : "Off");
     mrcpp::print::separator(0, '~', 2);
 
@@ -1343,9 +1358,21 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
         auto xc_cutoff = json_xcfunc["cutoff"];
         auto xc_funcs = json_xcfunc["functionals"];
         auto xc_order = order + 1;
+        // TODO: Look over and input parser so this is not necessary
+        std::string xc_lib;
+        if (json_fock.contains("xc_library")) {
+            if(json_fock["xc_library"].is_array()){
+                xc_lib = json_fock["xc_library"][0].get<std::string>();
+            }else{
+                xc_lib = json_fock["xc_library"]["xc_library"].get<std::string>();
+            }
+        }else{
+            xc_lib = "xcfun";
+        }
 
         mrdft::Factory xc_factory(*MRA);
         xc_factory.setSpin(xc_spin);
+        xc_factory.setLibxc((xc_lib == "libxc") ? true : false);
         xc_factory.setOrder(xc_order);
         xc_factory.setDensityCutoff(xc_cutoff);
         for (const auto &f : xc_funcs) {
@@ -1354,6 +1381,9 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
             xc_factory.setFunctional(name, coef);
         }
         auto mrdft_p = xc_factory.build();
+
+        mrdft_p->functional().print_functional_references();
+
         exx = mrdft_p->functional().amountEXX();
 
         if (order == 0) {

src/initial_guess/sad.cpp

@@ -72,6 +72,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     print_utils::text(0, "Screening   ", print_utils::dbl_to_str(screen, 5, true) + " StdDev");
     print_utils::text(0, "Restricted  ", (restricted) ? "True" : "False");
     print_utils::text(0, "Functional  ", "LDA (SVWN5)");
+    print_utils::text(0, "XC Library  ", (mrdft::Factory::libxc) ? "LibXC" : "XCFun");
     print_utils::text(0, "AO basis    ", "Hydrogenic orbitals");
     print_utils::text(0, "Zeta quality", std::to_string(zeta));
     mrcpp::print::separator(0, '~', 2);
@@ -90,8 +91,8 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     MomentumOperator p(D_p);
     NuclearOperator V_nuc(nucs, prec);
     CoulombOperator J(P_p);
-    XCOperator XC(mrdft_p);
-    RankZeroOperator V = V_nuc + J + XC;
+    XCOperator XC_(mrdft_p);
+    RankZeroOperator V = V_nuc + J + XC_;
 
     auto plevel = Printer::getPrintLevel();
     if (plevel == 1) mrcpp::print::header(1, "SAD Initial Guess");
@@ -103,7 +104,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     initial_guess::sad::project_atomic_densities(prec, rho_j, nucs, screen);
 
     // Compute XC density
-    Density &rho_xc = XC.getDensity(DensityType::Total);
+    Density &rho_xc = XC_.getDensity(DensityType::Total);
     mrcpp::deep_copy(rho_xc, rho_j);
     if (plevel == 1) mrcpp::print::time(1, "Projecting GTO density", t_lap);
 
@@ -151,6 +152,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     print_utils::text(0, "Screening   ", print_utils::dbl_to_str(screen, 5, true) + " StdDev");
     print_utils::text(0, "Restricted  ", (restricted) ? "True" : "False");
     print_utils::text(0, "Functional  ", "LDA (SVWN5)");
+    print_utils::text(0, "XC Library  ", (mrdft::Factory::libxc) ? "LibXC" : "XCFun");
     print_utils::text(0, "AO basis    ", "3-21G");
     mrcpp::print::separator(0, '~', 2);
 
@@ -167,8 +169,8 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     MomentumOperator p(D_p);
     NuclearOperator V_nuc(nucs, prec);
     CoulombOperator J(P_p);
-    XCOperator XC(mrdft_p);
-    RankZeroOperator V = V_nuc + J + XC;
+    XCOperator XC_(mrdft_p);
+    RankZeroOperator V = V_nuc + J + XC_;
 
     auto plevel = Printer::getPrintLevel();
     if (plevel == 1) mrcpp::print::header(1, "SAD Initial Guess");
@@ -180,7 +182,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c
     initial_guess::sad::project_atomic_densities(prec, rho_j, nucs, screen);
 
     // Compute XC density
-    Density &rho_xc = XC.getDensity(DensityType::Total);
+    Density &rho_xc = XC_.getDensity(DensityType::Total);
     mrcpp::deep_copy(rho_xc, rho_j);
     if (plevel == 1) mrcpp::print::time(1, "Projecting GTO density", t_lap);