Initial POD procedure

src/apps/molresponse/FrequencyResponse.cpp

@@ -584,7 +584,6 @@ std::pair<X_space, X_space> QuadraticResponse::compute_first_order_fock_matrix_t
         FAXB.y[i] = copy(world, transform(world, B.y[i], fax_dagger, true), true);
         FBXA.x[i] = copy(world, transform(world, A.x[i], fbx, true), true);
         FBXA.y[i] = copy(world, transform(world, A.y[i], fb_dagger, true), true);
-
     }
 
     FAXB.truncate();
@@ -830,3 +829,59 @@ X_space QuadraticResponse::compute_exchange_term(World &world, const X_space &A,
 
 
 //
+void PODResponse::compute_pod_modes(World &world) {
+    //
+    // 1. Convert the x_data into a response matrix
+    // 2. Compute the SVD of the response matrix
+    // 3. Compute the POD modes
+
+
+    // num calculations
+    auto m = x_data.size();
+    auto n = x_data[0].first.num_orbitals();
+
+
+    auto all_chi = X_space(world, m, n);
+    if (world.rank() == 0) print("m: ", m);
+    if (world.rank() == 0) print("n: ", n);
+
+    for (auto i = 0; i < m; i++) {
+
+
+        // first step is to copy x_data_i into response_matrix_i
+        auto x_data_i = x_data[i].first;
+        for (auto j = 0; j < n; j++) {
+            auto index = i * 3 + j;
+            all_chi.x[index] = copy(world, x_data_i.x[j]);
+            all_chi.y[index] = copy(world, x_data_i.y[j]);
+        }
+    }
+
+    auto A = inner(all_chi, all_chi);
+
+    // A = U * diag(s) * VT    for A real
+    // A = U * diag(s) * VH    for A complex
+
+    Tensor<double> VT(m, m);
+    Tensor<double> U(m, m);
+    Tensor<double> sigma(4);
+    svd(A, U, sigma, VT);
+
+    // create a json printing out the singular values
+    json j = {};
+    j["sigma"] = tensor_to_json(sigma);
+    j["VT"] = tensor_to_json(VT);
+    j["U"] = tensor_to_json(U);
+    j["A"] = tensor_to_json(A);
+
+    if (world.rank() == 0) print("sigma: ", sigma);
+    if (world.rank() == 0) print("VT: ", VT);
+    if (world.rank() == 0) print("U: ", U);
+    if (world.rank() == 0) print("A: ", A);
+
+
+    if (world.rank() == 0) {
+        std::ofstream o("pod.json");
+        o << std::setw(4) << j << std::endl;
+    }
+}

src/apps/molresponse/FrequencyResponse.hpp

@@ -216,6 +216,191 @@ class QuadraticResponse : public ResponseBase {
                                                     const X_space &phi0);
 };
 
+
+//Create a quadratic response class
+
+class PODResponse : public ResponseBase {
+
+    // A quadratic response class needs X_space vectors and one ground state.
+    // It will compute the 3rd order response property at the 3 given frequencies
+    // Beta(omegaA;omegaB,OmegaC)=tr(xA,vBC)+tr(muA,pBC)+tr(muA,qBC)
+    // Where xA, xB, xC are the response functions at the frequencies omegaA, omegaB, omegaC
+    // And 2nd order right hand side perturbation vector vBC(xB,xC) and
+    // pBC and qBC are the homogeneous components of the 2nd order density matrix response
+    // made entirely from first order vectors xB, xC
+    // vBC(xB,xC) = Q*(
+public:
+    PODResponse(World &world, const CalcParams &params, RHS_Generator rhs_generator)
+        : ResponseBase(world, params), generator(std::move(rhs_generator)) {
+        FunctionDefaults<3>::set_cubic_cell(-r_params.L(), r_params.L());
+        FunctionDefaults<3>::set_truncate_mode(1);
+        auto thresh = r_params.protocol()[0];
+
+        int k;
+        // Allow for imprecise conversion of threshold
+        if (thresh >= 0.9e-2) k = 4;
+        else if (thresh >= 0.9e-4)
+            k = 6;
+        else if (thresh >= 0.9e-6)
+            k = 8;
+        else if (thresh >= 0.9e-8)
+            k = 10;
+        else
+            k = 12;
+
+        // k defaults to make sense with thresh, override by providing k in
+        // input file
+        if (r_params.k() == -1) {
+            FunctionDefaults<3>::set_k(k);
+        } else {
+            FunctionDefaults<3>::set_k(r_params.k());
+        }
+
+        // Set Function Defaults
+        FunctionDefaults<3>::set_thresh(thresh);
+        FunctionDefaults<3>::set_refine(true);
+        FunctionDefaults<3>::set_initial_level(2);
+
+
+        FunctionDefaults<3>::set_autorefine(false);
+        FunctionDefaults<3>::set_apply_randomize(false);
+        FunctionDefaults<3>::set_project_randomize(false);
+        GaussianConvolution1DCache<double>::map.clear();
+
+        double safety = 0.1;
+        vtol = FunctionDefaults<3>::get_thresh() * safety;
+        shared_coulomb_operator = poperatorT(CoulombOperatorPtr(world, r_params.lo(), thresh));
+        gradop = gradient_operator<double, 3>(world);
+        potential_manager = std::make_shared<PotentialManager>(molecule, "a");
+        potential_manager->make_nuclear_potential(world);
+        // GaussianConvolution1DCache<double>::map.clear();//(TODO:molresponse-What
+        // Create the masking function
+        mask = real_function_3d(real_factory_3d(world).f(mask3).initial_level(4).norefine());
+
+        ground_density = make_ground_density(world);
+        ground_density.truncate(FunctionDefaults<3>::get_thresh());
+        // Basic print
+        if (world.rank() == 0) {
+            print("\nSolving NDIM=", 3, " with thresh", thresh, "    k", FunctionDefaults<3>::get_k(), "  dconv",
+                  std::max(thresh, r_params.dconv()), "\n");
+        }
+        if (world.rank() == 0) { print("Successfully set protocol"); }
+
+        // ground state orbitals change
+        bool redo = false;
+        // Verify ground state orbitals have correct k
+        if (FunctionDefaults<3>::get_k() != ground_orbitals[0].k()) {
+            // Re-read orbitals from the archive (assuming
+            // the archive has orbitals stored at a higher
+            if (world.rank() == 0) { print("check k: ground orbitals"); }
+            // k value than what was previously computed
+            ground_calc.read(world);
+            if (world.rank() == 0) { print("check k: read ground orbitals"); }
+            // k value than what was previously computed
+            reconstruct(world, ground_orbitals);
+            if (world.rank() == 0) { print("check k: reconstruct ground orbitals"); }
+            // Reset correct k (its set in g_params.read)
+            FunctionDefaults<3>::set_k(k);
+            // Project each ground state to correct k
+            for (auto &orbital: ground_orbitals) {
+                orbital = project(orbital, FunctionDefaults<3>::get_k(), thresh, false);
+            }
+            world.gop.fence();
+            if (world.rank() == 0) { print("check k: project ground orbitals"); }
+            // Clean up a bit
+            truncate(world, ground_orbitals);
+            if (world.rank() == 0) { print("check k: truncate ground orbitals"); }
+            // Now that ground orbitals have correct k lets make the ground density
+            // again
+            ground_density = make_ground_density(world);
+            ground_density.truncate(FunctionDefaults<3>::get_thresh());
+            if (world.rank() == 0) { print("check k: make ground density"); }
+            // Ground state orbitals changed, clear old hamiltonian
+            redo = true;
+        }
+        // Recalculate ground state hamiltonian here
+        if (redo or !hamiltonian.has_data()) {
+            if (world.rank() == 0) { print("check k: re-do hamiltonian"); }
+            auto [HAM, HAM_NO_DIAG] = ComputeHamiltonianPair(world);
+            if (world.rank() == 0) { print("check k: output hamiltonian"); }
+            // TODO this doesn't seem right...
+            hamiltonian = HAM;
+            ham_no_diag = HAM_NO_DIAG;
+        }
+
+        // If we stored the potential, check that too
+        if (r_params.store_potential()) {
+            if (FunctionDefaults<3>::get_k() != stored_potential[0][0].k()) {
+                // Project the potential into correct k
+                for (auto &potential_vector: stored_potential) {
+                    reconstruct(world, potential_vector);
+                    for (auto &vi: potential_vector) { vi = project(vi, FunctionDefaults<3>::get_k(), thresh, false); }
+                    world.gop.fence();
+                }
+            }
+            if (FunctionDefaults<3>::get_k() != stored_v_coul.k())
+                stored_v_coul = project(stored_v_coul, FunctionDefaults<3>::get_k(), thresh, false);
+            if (FunctionDefaults<3>::get_k() != stored_v_nuc.k())
+                stored_v_nuc = project(stored_v_nuc, FunctionDefaults<3>::get_k(), thresh, false);
+        }
+        // project the mask
+        if (FunctionDefaults<3>::get_k() != mask.k()) {
+            mask = project(mask, FunctionDefaults<3>::get_k(), thresh, false);
+            if (world.rank() == 0) { print("check k: project mask"); }
+        }
+        if (world.rank() == 0) { print("check k: project Chi"); }
+        // Make sure everything is done before leaving
+        world.gop.fence();
+    }
+
+
+    void set_x_data(World &world, const std::vector<double> &freqABC, const std::vector<path> &restart_files) {
+        this->frequencies = freqABC;
+        if (freqABC.size() != 3) { throw std::runtime_error("Quadratic response requires 3 freqABC"); }
+        // print ABC frequencies
+        print("ABC frequencies: ", freqABC[0], " ", freqABC[1], " ", freqABC[2]);
+
+        for (size_t i = 0; i < freqABC.size(); i++) {
+            auto omega = freqABC[i];
+
+            if (omega == 0.0) {
+                frequency_contexts[i].set_strategy(
+                        std::make_unique<static_inner_product>(), std::make_unique<J1StrategyStable>(),
+                        std::make_unique<K1StrategyStatic>(), std::make_unique<VXC1StrategyStandard>(),
+                        std::make_unique<StaticDensityStrategy>(), std::make_unique<LoadFrequencyXSpace>());
+            } else {
+                frequency_contexts[i].set_strategy(
+                        std::make_unique<full_inner_product>(), std::make_unique<J1StrategyStable>(),
+                        std::make_unique<K1StrategyFull>(), std::make_unique<VXC1StrategyStandard>(),
+                        std::make_unique<FullDensityStrategy>(), std::make_unique<LoadFrequencyXSpace>());
+            }
+
+            // print omega before load_x_space
+            print("omega: ", omega);
+            x_data[i] = frequency_contexts[i].load_x_space(world, restart_files[i], r_params, omega);
+            ::check_k(world, x_data[i].first, FunctionDefaults<3>::get_thresh(), r_params.k());
+        }
+    };
+
+    void initialize(World &world) override;
+
+    void load(World &world, const std::string &name) override;
+
+    void save(World &world, const std::string &name) override;
+    void iterate(World &world) override;
+
+    void compute_pod_modes(World &world);
+
+
+private:
+    std::vector<Context> frequency_contexts;
+    std::vector<double> frequencies;
+    std::vector<XData> x_data;
+    std::pair<X_space, X_space> setup_XBC(World &world);
+    RHS_Generator generator;
+};
+
+
 class FrequencyResponse : public ResponseBase {
 
 public:

src/apps/molresponse/testing/CMakeLists.txt

@@ -15,6 +15,7 @@ set(MY_EXECUTABLES
         test_development
         x_space_testing
         mad-quadratic
+        pod_response
         )
 
 add_mad_executable(test_schema_json "madness_catch_main.cc;qcschema_json_testing.cpp" "MADresponse;MADchem")

src/apps/molresponse/testing/pod_response.cpp

@@ -0,0 +1,80 @@
+//
+// Created by adrianhurtado on 1/1/22.
+//
+#include "ResponseExceptions.hpp"
+#include "madness/tensor/tensor_json.hpp"
+#include "madness/world/worldmem.h"
+#include "response_functions.h"
+#include "runners.hpp"
+
+#if defined(HAVE_SYS_TYPES_H) && defined(HAVE_SYS_STAT_H) && defined(HAVE_UNISTD_H)
+
+#include <sys/stat.h>
+#include <unistd.h>
+
+static inline int file_exists(const char *input_name) {
+    struct stat buffer {};
+    size_t rc = stat(input_name, &buffer);
+    return (rc == 0);
+}
+
+#endif
+
+using path = std::filesystem::path;
+
+using namespace madness;
+
+
+auto main(int argc, char *argv[]) -> int {
+
+    madness::initialize(argc, argv);
+
+    {
+        World world(SafeMPI::COMM_WORLD);
+        startup(world, argc, argv, true);
+
+        try {
+            sleep(5);
+            std::cout.precision(6);
+            if (argc != 6) {
+                std::cout << "Wrong number of inputs" << std::endl;
+                return 1;
+            }
+            const std::string molecule_name{argv[1]};
+            const std::string xc{argv[2]};
+            const std::string op{argv[3]};
+            const std::string precision{argv[4]};
+            const std::string static_calc{argv[5]};
+
+
+            auto schema = runSchema(world, xc);
+            auto m_schema = moldftSchema(world, molecule_name, xc, schema);
+            auto f_schema = frequencySchema(world, schema, m_schema, op, static_calc == "true");
+
+
+            runPODResponse(world, f_schema, precision);
+            world.gop.fence();
+        } catch (const SafeMPI::Exception &e) {
+            print(e.what());
+            error("caught an MPI exception");
+        } catch (const madness::MadnessException &e) {
+            print(e);
+            error("caught a MADNESS exception");
+        } catch (const madness::TensorException &e) {
+            print(e.what());
+            error("caught a Tensor exception");
+        } catch (const nlohmann::detail::exception &e) {
+            print(e.what());
+            error("Caught JSON exception");
+        } catch (const std::filesystem::filesystem_error &ex) {
+            std::cerr << ex.what() << "\n";
+        } catch (const std::exception &e) {
+            print(e.what());
+            error("caught an STL exception");
+        } catch (...) { error("caught unhandled exception"); }
+        // Nearly all memory will be freed at this point
+        print_stats(world);
+    }
+    finalize();
+    return 0;
+}