diff --git a/src/apps/molresponse/CMakeLists.txt b/src/apps/molresponse/CMakeLists.txt index 12149df3054..a159f2195c5 100644 --- a/src/apps/molresponse/CMakeLists.txt +++ b/src/apps/molresponse/CMakeLists.txt @@ -16,6 +16,7 @@ set(MOLRESPONSE_SOURCES property.cc global_functions.cc timer.cc + Plot_VTK.cc basic_operators.cc response_parameters.cpp x_space.cc @@ -28,6 +29,7 @@ set(MOLRESPONSE_HEADERS x_space.h global_functions.h timer.h + Plot_VTK.h basic_operators.h response_parameters.h ) @@ -79,6 +81,7 @@ set(RALL_SOURCES property.cc global_functions.cc timer.cc + Plot_VTK.cc basic_operators.cc response_parameters.cpp x_space.cc @@ -92,6 +95,7 @@ set(RALL_HEADERS global_functions.h timer.h basic_operators.h + Plot_VTK.h response_parameters.h x_space.h ) diff --git a/src/apps/molresponse/Plot_VTK.cc b/src/apps/molresponse/Plot_VTK.cc index f04eb76d76d..9a909782f6a 100644 --- a/src/apps/molresponse/Plot_VTK.cc +++ b/src/apps/molresponse/Plot_VTK.cc @@ -12,15 +12,22 @@ */ #include "Plot_VTK.h" - #include - -#include -//#include #include #include -#include "../chem/molecule.h" +#if defined(__has_include) +#if __has_include() +#define MADCHEM_HAS_STD_FILESYSTEM +// is not reliably usable on Linux with gcc < 9 +#if defined(__GNUC__) +#if __GNUC__ >= 7 && __GNUC__ < 9 +#undef MADCHEM_HAS_STD_FILESYSTEM +#endif +#endif +#if defined(MADCHEM_HAS_STD_FILESYSTEM) + +#include namespace madness { void write_molecules_to_file(const Molecule &molecule, const std::string &geo_file) { @@ -36,11 +43,8 @@ namespace madness { // Write the data size_t Natoms = molecule.natom(); for (size_t i = 0; i < Natoms; i++) { - std::fprintf(f, - "%5s %16.12f %16.12f %16.12f\n", - atomic_number_to_symbol(molecule.get_atomic_number(i)).c_str(), - coords[i][0], - coords[i][1], + std::fprintf(f, "%5s %16.12f %16.12f %16.12f\n", + atomic_number_to_symbol(molecule.get_atomic_number(i)).c_str(), coords[i][0], coords[i][1], coords[i][2]); } @@ -48,7 +52,9 @@ namespace madness { fclose(f); } - void do_response_orbital_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, const vector_real_function_3d &ground_orbs, const response_matrix &responseMatrix) { + void do_response_orbital_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, + const vector_real_function_3d &ground_orbs, + const response_matrix &responseMatrix) { // Stuff needed to plot // Vector box_lo{-L / 4, -L / 4, -L / 4}; @@ -81,15 +87,7 @@ namespace madness { auto orb_file = vtk_dir + "/phi0_" + std::to_string(orb_num) + ".vts"; filename = orb_file.c_str(); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(phi0_i, - "ground_orbtial", - world, - filename, - box_lo, - box_hi, - points, - true, - false); + plotvtk_data(phi0_i, "ground_orbtial", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); orb_num++; }); @@ -104,14 +102,7 @@ namespace madness { auto orb_file = vtk_dir + "/" + "x" + std::to_string(state_number) + std::to_string(orb_num) + ".vts"; filename = orb_file.c_str(); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(xi, - "response_x_orbitals", - world, - filename, - box_lo, - box_hi, - points, - true, + plotvtk_data(xi, "response_x_orbitals", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); orb_num++; @@ -121,14 +112,7 @@ namespace madness { auto orb_file = vtk_dir + "/" + "y" + std::to_string(state_number) + std::to_string(orb_num) + ".vts"; filename = orb_file.c_str(); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(yi, - "response_y_orbitals", - world, - filename, - box_lo, - box_hi, - points, - true, + plotvtk_data(yi, "response_y_orbitals", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); orb_num++; @@ -136,7 +120,9 @@ namespace madness { state_number++; }); } - void do_response_density_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, const real_function_3d &ground_density, const vector_real_function_3d &response_density) { + void do_response_density_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, + const real_function_3d &ground_density, + const vector_real_function_3d &response_density) { // Stuff needed to plot // Vector box_lo{-L / 4, -L / 4, -L / 4}; @@ -156,15 +142,7 @@ namespace madness { auto density_file = vtk_dir + "/" + "rho0.vts"; filename = density_file.c_str(); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(ground_density, - "ground_density", - world, - filename, - box_lo, - box_hi, - points, - true, - false); + plotvtk_data(ground_density, "ground_density", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); //***********************************ground density plot @@ -174,27 +152,13 @@ namespace madness { filename = density_file.c_str(); auto field_name = "r_density_" + std::to_string(state_number); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(rho_i, - field_name.c_str(), - world, - filename, - box_lo, - box_hi, - points, - true, - false); + plotvtk_data(rho_i, field_name.c_str(), world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); state_number++; }); } - void do_vtk_plots(World &world, - int npt_plot, - double L, - int lowest_orbital, - int highest_orbital, - const Molecule &molecule, - std::vector densities, - const std::string &name) { + void do_vtk_plots(World &world, int npt_plot, double L, int lowest_orbital, int highest_orbital, + const Molecule &molecule, std::vector densities, const std::string &name) { // Stuff needed to plot // Vector box_lo{-L, -L, -L}; @@ -239,15 +203,7 @@ namespace madness { // npts // binary // plot refile - plotvtk_data(densities[i], - "density", - world, - filename, - box_lo, - box_hi, - points, - true, - false); + plotvtk_data(densities[i], "density", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); } std::string b; @@ -256,24 +212,11 @@ namespace madness { b = vtk_dir + "/" + "total-electrondensity.vts"; filename = b.c_str(); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(rho, - "total-electrondensity", - world, - filename, - box_lo, - box_hi, - points, - true, - false); + plotvtk_data(rho, "total-electrondensity", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); } - void do_vtk_plots(World &world, - int npt_plot, - double L, - Molecule molecule, - real_function_3d &rho_0, - std::vector &rho_omega, - std::vector &ground_orbitals, + void do_vtk_plots(World &world, int npt_plot, double L, Molecule molecule, real_function_3d &rho_0, + std::vector &rho_omega, std::vector &ground_orbitals, X_space &Chi) { std::string vtk_dir = "vtk_plots"; std::filesystem::create_directories(vtk_dir); @@ -297,11 +240,8 @@ namespace madness { // Write the data size_t Natoms = molecule.natom(); for (size_t i = 0; i < Natoms; i++) { - std::fprintf(f, - "%5s %16.12f %16.12f %16.12f\n", - atomic_number_to_symbol(molecule.get_atomic_number(i)).c_str(), - coords[i][0], - coords[i][1], + std::fprintf(f, "%5s %16.12f %16.12f %16.12f\n", + atomic_number_to_symbol(molecule.get_atomic_number(i)).c_str(), coords[i][0], coords[i][1], coords[i][2]); } // Clean up @@ -310,15 +250,7 @@ namespace madness { filename = rho0_file.c_str(); plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(rho_0, - "ground_density", - world, - filename, - box_lo, - box_hi, - points, - true, - false); + plotvtk_data(rho_0, "ground_density", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); // ground orbitals std::string g_orb_file = "phi_"; @@ -328,14 +260,7 @@ namespace madness { filename = fname.c_str(); // VTK plotting stuff plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(ground_orbitals[i], - "ground_orbital", - world, - filename, - box_lo, - box_hi, - points, - true, + plotvtk_data(ground_orbitals[i], "ground_orbital", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); } @@ -343,19 +268,11 @@ namespace madness { std::string x_orb_file = "x_"; for (size_t i = 0; i < Chi.num_states(); ++i) { for (size_t j = 0; j < Chi.num_orbitals(); ++j) { - fname = vtk_dir + "/" + x_orb_file + std::to_string(i) + "_" + - std::to_string(j) + ".vts"; + fname = vtk_dir + "/" + x_orb_file + std::to_string(i) + "_" + std::to_string(j) + ".vts"; filename = fname.c_str(); // VTK plotting stuff plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(Chi.x[i][j], - "x_orbitals", - world, - filename, - box_lo, - box_hi, - points, - true, + plotvtk_data(Chi.x[i][j], "x_orbitals", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); } @@ -363,19 +280,11 @@ namespace madness { std::string y_orb_file = "y_"; for (size_t i = 0; i < Chi.num_states(); ++i) { for (size_t j = 0; j < Chi.num_orbitals(); ++j) { - fname = vtk_dir + "/" + y_orb_file + std::to_string(i) + "_" + - std::to_string(j) + ".vts"; + fname = vtk_dir + "/" + y_orb_file + std::to_string(i) + "_" + std::to_string(j) + ".vts"; filename = fname.c_str(); // VTK plotting stuff plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(Chi.y[i][j], - "y_orbitals", - world, - filename, - box_lo, - box_hi, - points, - true, + plotvtk_data(Chi.y[i][j], "y_orbitals", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); } @@ -387,14 +296,7 @@ namespace madness { filename = fname.c_str(); // VTK plotting stuff plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true); - plotvtk_data(rho_omega[i], - "transition_density", - world, - filename, - box_lo, - box_hi, - points, - true, + plotvtk_data(rho_omega[i], "transition_density", world, filename, box_lo, box_hi, points, true, false); plotvtk_end<3>(world, filename, true); } @@ -405,4 +307,11 @@ namespace madness { // // } + }// namespace madness + +#endif +#endif +#endif + +#include "../chem/molecule.h" diff --git a/src/apps/molresponse/Plot_VTK.h b/src/apps/molresponse/Plot_VTK.h index fe2fdba2be6..81fd50824a2 100644 --- a/src/apps/molresponse/Plot_VTK.h +++ b/src/apps/molresponse/Plot_VTK.h @@ -22,31 +22,41 @@ #include "../chem/molecule.h" #include "x_space.h" + +#if defined(__has_include) +#if __has_include() +#define MADCHEM_HAS_STD_FILESYSTEM +// is not reliably usable on Linux with gcc < 9 +#if defined(__GNUC__) +#if __GNUC__ >= 7 && __GNUC__ < 9 +#undef MADCHEM_HAS_STD_FILESYSTEM +#endif +#endif + +#if defined(MADCHEM_HAS_STD_FILESYSTEM) + +#include namespace madness { void write_molecules_to_file(const Molecule &molecule, const std::string &geo_file); - void do_response_orbital_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, const vector_real_function_3d &ground_orbs, const response_matrix &responseMatrix); - void do_response_density_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, const real_function_3d &ground_density, const vector_real_function_3d &response_density); - - void do_vtk_plots(World &world, - int npt_plot, - double L, - int lowest_orbital, - int highest_orbital, - const Molecule &molecule, - std::vector densities, - const std::string &name); - - void do_vtk_plots(World &world, - int npt_plot, - double L, - Molecule molecule, - real_function_3d &rho_0, - std::vector &rho_omega, - std::vector &ground_orbitals, + void do_response_orbital_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, + const vector_real_function_3d &ground_orbs, + const response_matrix &responseMatrix); + void do_response_density_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule, + const real_function_3d &ground_density, + const vector_real_function_3d &response_density); + + void do_vtk_plots(World &world, int npt_plot, double L, int lowest_orbital, int highest_orbital, + const Molecule &molecule, std::vector densities, const std::string &name); + + void do_vtk_plots(World &world, int npt_plot, double L, Molecule molecule, real_function_3d &rho_0, + std::vector &rho_omega, std::vector &ground_orbitals, X_space &Chi); + }// namespace madness +#endif +#endif +#endif #endif// SRC_APPS_molresponse_PLOT_VTK_H_ -// Dueces diff --git a/src/apps/molresponse/ResponseBase.cpp b/src/apps/molresponse/ResponseBase.cpp index 057260bece5..b70c7078da9 100644 --- a/src/apps/molresponse/ResponseBase.cpp +++ b/src/apps/molresponse/ResponseBase.cpp @@ -17,7 +17,69 @@ #if defined(MADCHEM_HAS_STD_FILESYSTEM) #include +#include + +void ResponseBase::plotResponseOrbitals(World &world, size_t iteration, const response_space &x_response, + const response_space &y_response, ResponseParameters const &responseParameters, + GroundStateCalculation const &g_params) { + std::filesystem::create_directories("plots/densities"); + std::filesystem::create_directory("plots/orbitals"); + + // TESTING + // get transition density + // num orbitals + size_t n = x_response[0].size(); + size_t m = x_response.size(); + + real_function_3d rho0 = dot(world, ground_orbitals, ground_orbitals); + std::vector rho1 = transition_density(world, ground_orbitals, x_response, y_response); + std::string dir("xyz"); + // for plot_name size + size_t buffSize = 500; + char plot_name[buffSize]; + double Lp = std::min(responseParameters.L(), 24.0); + // Doing line plots along each axis + for (int d = 0; d < 3; d++) { + // print ground_state + plotCoords plt(d, Lp); + // plot ground density + if (iteration == 1) { + snprintf(plot_name, buffSize, "plots/densities/rho0_%c_0.plt", dir[d]); + plot_line(plot_name, 5001, plt.lo, plt.hi, rho0); + } + // plot ground orbitals and transition densities in xyz directions + snprintf(plot_name, buffSize, "plots/orbitals/phi0_%c_0.plt", dir[d]); + plot_line(plot_name, 5001, plt.lo, plt.hi, ground_orbitals); + snprintf(plot_name, buffSize, "plots/densities/rho1_%c.plt", dir[d]); + plot_line(plot_name, 5001, plt.lo, plt.hi, rho1); + + for (int b = 0; b < static_cast(m); b++) { + // plot x and y transition orbitals + snprintf(plot_name, buffSize, "plots/orbitals/phix_%c_%d.plt", dir[d], static_cast(b)); + plot_line(plot_name, 5001, plt.lo, plt.hi, x_response[b]); + snprintf(plot_name, buffSize, "plots/orbitals/phiy_%c_%d.plt", dir[d], static_cast(b)); + plot_line(plot_name, 5001, plt.lo, plt.hi, y_response[b]); + } + } + world.gop.fence(); + // END TESTING +} + +void PlotGroundDensityVTK(World &world, const ResponseBase &calc) { + auto [ground_calc, molecule, r_params] = calc.get_parameter(); + auto ground_orbitals = calc.get_orbitals(); + + if (r_params.plot_initial()) { + if (world.rank() == 0) print("\n Plotting ground state densities.\n"); + if (r_params.plot_l() > 0.0) + do_vtk_plots(world, int(r_params.plot_pts()), r_params.plot_l(), 0, int(r_params.num_orbitals()), molecule, + square(world, ground_orbitals), "ground"); + else + do_vtk_plots(world, int(r_params.plot_pts()), r_params.L() / 2.0, 0, int(r_params.num_orbitals()), molecule, + square(world, ground_orbitals), "ground"); + } +} #endif #endif #endif @@ -28,18 +90,15 @@ /// \param world /// \param params ResponseBase::ResponseBase(World &world, const CalcParams ¶ms) - : r_params(params.response_parameters), molecule(params.molecule), - ground_calc(params.ground_calculation), - ground_orbitals(ground_calc.orbitals()), - ground_energies(ground_calc.get_energies()), + : r_params(params.response_parameters), molecule(params.molecule), ground_calc(params.ground_calculation), + ground_orbitals(ground_calc.orbitals()), ground_energies(ground_calc.get_energies()), Chi(world, r_params.num_states(), r_params.num_orbitals()) { // Broadcast to all the other nodes world.gop.broadcast_serializable(r_params, 0); world.gop.broadcast_serializable(ground_energies, 0); world.gop.broadcast_serializable(molecule, 0); - xcf.initialize(r_params.xc(), !r_params.spinrestricted(), world, - r_params.print_level() >= 3); + xcf.initialize(r_params.xc(), !r_params.spinrestricted(), world, r_params.print_level() >= 3); r_params.to_json(j_molresponse); // Set the Box Size and Truncate Mode @@ -68,15 +127,12 @@ void ResponseBase::check_k(World &world, double thresh, int k) { 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"); - } + 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); + orbital = project(orbital, FunctionDefaults<3>::get_k(), thresh, false); } world.gop.fence(); if (world.rank() == 0) { print("check k: project ground orbitals"); } @@ -107,19 +163,14 @@ void ResponseBase::check_k(World &world, double thresh, int 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); - } + 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); + 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); + stored_v_nuc = project(stored_v_nuc, FunctionDefaults<3>::get_k(), thresh, false); } // project the mask if (FunctionDefaults<3>::get_k() != mask.k()) { @@ -141,8 +192,7 @@ void ResponseBase::check_k(World &world, double thresh, int k) { // /// \param world /// \return -auto ResponseBase::ComputeHamiltonianPair(World &world) const - -> std::pair, Tensor> { +auto ResponseBase::ComputeHamiltonianPair(World &world) const -> std::pair, Tensor> { // Basic output if (r_params.print_level() >= 1) molresponse::start_timer(world); auto phi = ground_orbitals; @@ -170,9 +220,7 @@ auto ResponseBase::ComputeHamiltonianPair(World &world) const // Note: No negative as the formula above // has one as well, so they cancel Tensor T = - 1.0 / 2.0 * - (matrix_inner(world, fx, fx) + matrix_inner(world, fy, fy) + - matrix_inner(world, fz, fz)); + 1.0 / 2.0 * (matrix_inner(world, fx, fx) + matrix_inner(world, fy, fy) + matrix_inner(world, fz, fz)); // Construct phiVphi // v_nuc first @@ -215,23 +263,20 @@ auto ResponseBase::ComputeHamiltonianPair(World &world) const for (const auto &phi_i: phi) { /// Multiplies a function against a vector of functions using sparsity of a /// and v[i] --- q[i] = a * v[i] - auto psif = mul_sparse(world, phi_i, phi, - FunctionDefaults<3>::get_thresh()); + auto psif = mul_sparse(world, phi_i, phi, FunctionDefaults<3>::get_thresh()); truncate(world, psif); psif = apply(world, *op, psif); truncate(world, psif); // Save the potential here if we are saving it if (r_params.store_potential()) { stored_potential.push_back(psif); } - psif = mul_sparse(world, phi_i, psif, - FunctionDefaults<3>::get_thresh()); + psif = mul_sparse(world, phi_i, psif, FunctionDefaults<3>::get_thresh()); gaxpy(world, 1.0, Kphi, 1.0, psif); } // Only use the exchange above if HF: Tensor phiVphi; if (r_params.xc() == "hf") { // Construct phiVphi - phiVphi = matrix_inner(world, phi, v_phi0) - - matrix_inner(world, phi, Kphi); + phiVphi = matrix_inner(world, phi, v_phi0) - matrix_inner(world, phi, Kphi); } else {// DFT XCOperator xcop = make_xc_operator(world); @@ -250,14 +295,10 @@ auto ResponseBase::ComputeHamiltonianPair(World &world) const auto new_hamiltonian = T + phiVphi; for (int64_t i = 0; i < new_hamiltonian.dim(0); i++) { - for (int64_t j = i + 1; j < new_hamiltonian.dim(1); j++) { - new_hamiltonian(j, i) = new_hamiltonian(i, j); - } + for (int64_t j = i + 1; j < new_hamiltonian.dim(1); j++) { new_hamiltonian(j, i) = new_hamiltonian(i, j); } } double traceOfHamiltonian(0); - for (int64_t i = 0; i < new_hamiltonian.dim(0); i++) { - traceOfHamiltonian += new_hamiltonian(i, i); - } + for (int64_t i = 0; i < new_hamiltonian.dim(0); i++) { traceOfHamiltonian += new_hamiltonian(i, i); } if (world.rank() == 0) { print("Trace of Hamiltonian"); print(traceOfHamiltonian); @@ -266,12 +307,10 @@ auto ResponseBase::ComputeHamiltonianPair(World &world) const // (T+phiVphi) - Lambda * eye // Copy new_hamiltonian and zero the diagonal auto new_hamiltonian_no_diag = copy(new_hamiltonian); - for (size_t i = 0; i < num_orbitals; i++) - new_hamiltonian_no_diag(long(i), long(i)) = 0.0; + for (size_t i = 0; i < num_orbitals; i++) new_hamiltonian_no_diag(long(i), long(i)) = 0.0; // End timer - if (r_params.print_level() >= 1) - molresponse::end_timer(world, " Create grnd ham:"); + if (r_params.print_level() >= 1) molresponse::end_timer(world, " Create grnd ham:"); return {new_hamiltonian, new_hamiltonian_no_diag}; } @@ -301,14 +340,11 @@ auto ResponseBase::Coulomb(World &world) const -> real_function_3d { // TODO Create apply_operator(f) generalized function in place of coulomb -auto ResponseBase::make_xc_operator(World &world) const - -> XCOperator { +auto ResponseBase::make_xc_operator(World &world) const -> XCOperator { return {world, r_params.xc(), false, ground_density, ground_density}; } -auto ResponseBase::update_density(World &world, const X_space &chi, - const vecfuncT &old_density) const - -> vecfuncT { +auto ResponseBase::update_density(World &world, const X_space &chi, const vecfuncT &old_density) const -> vecfuncT { auto density = copy(world, old_density); auto calc_type = r_params.calc_type(); auto thresh = FunctionDefaults<3>::get_thresh(); @@ -338,8 +374,7 @@ auto ResponseBase::update_density(World &world, const X_space &chi, return density; } -auto ResponseBase::make_density(World &world, const X_space &chi) const - -> vecfuncT { +auto ResponseBase::make_density(World &world, const X_space &chi) const -> vecfuncT { auto density = vector_real_function_3d(chi.num_states()); auto calc_type = r_params.calc_type(); auto thresh = FunctionDefaults<3>::get_thresh(); @@ -377,28 +412,21 @@ void ResponseBase::load_balance_chi(World &world) { real_function_3d v_nuclear; v_nuclear = potential_manager->vnuclear(); for (auto &b: Chi.active) { - for (auto &xij: Chi.x[b]) { - lb.add_tree(xij, lbcost(1.0, 8.0), false); - } + for (auto &xij: Chi.x[b]) { lb.add_tree(xij, lbcost(1.0, 8.0), false); } } if (r_params.omega() != 0) { for (auto &b: Chi.active) { - for (auto &yij: Chi.y[b]) { - lb.add_tree(yij, lbcost(1.0, 8.0), false); - } + for (auto &yij: Chi.y[b]) { lb.add_tree(yij, lbcost(1.0, 8.0), false); } } } world.gop.fence(); - FunctionDefaults<3>::redistribute( - world, - lb.load_balance( - r_params.loadbalparts()));// 6.0 needs retuning after + FunctionDefaults<3>::redistribute(world, + lb.load_balance(r_params.loadbalparts()));// 6.0 needs retuning after world.gop.fence(); molresponse::end_timer(world, "Load balancing"); } -auto ResponseBase::make_bsh_operators_response(World &world, double &shift, - const double omega) const +auto ResponseBase::make_bsh_operators_response(World &world, double &shift, const double omega) const -> std::vector { if (r_params.print_level() >= 1) molresponse::start_timer(world); @@ -410,90 +438,61 @@ auto ResponseBase::make_bsh_operators_response(World &world, double &shift, int p = 0; std::for_each(ops.begin(), ops.end(), [&](auto &operator_p) { double mu = sqrt(-2.0 * (ground_energies(p++) + omega + shift)); - operator_p = - poperatorT(BSHOperatorPtr3D(world, mu, r_params.lo(), tol)); + operator_p = poperatorT(BSHOperatorPtr3D(world, mu, r_params.lo(), tol)); }); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "make bsh operators response"); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "make bsh operators response"); } return ops; // End timer } -auto ResponseBase::compute_gamma(World &world, const gamma_orbitals &density, - const XCOperator &xc) const +auto ResponseBase::compute_gamma(World &world, const gamma_orbitals &density, const XCOperator &xc) const -> X_space { - std::shared_ptr>> old_pmap = - FunctionDefaults<3>::get_pmap(); + std::shared_ptr>> old_pmap = FunctionDefaults<3>::get_pmap(); auto c_xc = xcf.hf_exchange_coefficient(); - auto [chi_alpha, phi0, rho1] = - orbital_load_balance(world, density, r_params.loadbalparts()); + auto [chi_alpha, phi0, rho1] = orbital_load_balance(world, density, r_params.loadbalparts()); QProjector projector(world, phi0); - X_space W = X_space::zero_functions(world, chi_alpha.num_states(), - chi_alpha.num_orbitals()); + X_space W = X_space::zero_functions(world, chi_alpha.num_states(), chi_alpha.num_orbitals()); auto apply_projector = [&](auto &xi) { return projector(xi); }; // apply the exchange kernel to rho if necessary if (r_params.print_level() >= 1) { molresponse::start_timer(world); } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_make_density"); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma_make_density"); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } - auto J = response_context.compute_j1(world, chi_alpha, rho1, phi0, - shared_coulomb_operator); - inner_to_json(world, "j1", response_context.inner(chi_alpha, J), - iter_function_data); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); - } + auto J = response_context.compute_j1(world, chi_alpha, rho1, phi0, shared_coulomb_operator); + inner_to_json(world, "j1", response_context.inner(chi_alpha, J), iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); } // Create Coulomb potential on ground_orbitals if (r_params.print_level() >= 1) { molresponse::start_timer(world); } auto K = response_context.compute_k1(world, chi_alpha, phi0); - inner_to_json(world, "k1", response_context.inner(chi_alpha, K), - iter_function_data); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); - } - if (r_params.print_level() >= 20) { - print_inner(world, "old xKx", chi_alpha, K); - } + inner_to_json(world, "k1", response_context.inner(chi_alpha, K), iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); } + if (r_params.print_level() >= 20) { print_inner(world, "old xKx", chi_alpha, K); } molresponse::start_timer(world); - X_space gamma = X_space::zero_functions(world, chi_alpha.num_states(), - chi_alpha.num_orbitals()); + X_space gamma = X_space::zero_functions(world, chi_alpha.num_states(), chi_alpha.num_orbitals()); gamma = 2 * J - c_xc * K; if (c_xc != 1.0) { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } W = response_context.compute_VXC1(world, chi_alpha, rho1, phi0, xc); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "XC[omega]", "XC[omega]", - iter_timing); - } - inner_to_json(world, "v1_xc", response_context.inner(chi_alpha, W), - iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "XC[omega]", "XC[omega]", iter_timing); } + inner_to_json(world, "v1_xc", response_context.inner(chi_alpha, W), iter_function_data); gamma += (1.0 - c_xc) * W; } //gamma.truncate(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma add", "gamma_truncate_add", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma add", "gamma_truncate_add", iter_timing); } // project out ground state if (r_params.print_level() >= 1) { molresponse::start_timer(world); } gamma = oop_apply(gamma, apply_projector); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_project", "gamma_project", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma_project", "gamma_project", iter_timing); } if (r_params.print_level() >= 20) { molresponse::start_timer(world); print_inner(world, "xJx", chi_alpha, J); @@ -518,21 +517,15 @@ auto ResponseBase::compute_gamma(World &world, const gamma_orbitals &density, } -auto ResponseBase::compute_theta_X(World &world, const X_space &chi, - const vector_real_function_3d &rho1, - const XCOperator &xc, - const std::string &calc_type) const - -> X_space { +auto ResponseBase::compute_theta_X(World &world, const X_space &chi, const vector_real_function_3d &rho1, + const XCOperator &xc, const std::string &calc_type) const -> X_space { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } bool compute_Y = calc_type == "full"; X_space Theta_X = X_space(world, chi.num_states(), chi.num_orbitals()); if (r_params.print_level() >= 1) { molresponse::start_timer(world); } X_space V0X = compute_V0X(world, chi, xc, compute_Y); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "compute_V0X", "compute_V0X", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "compute_V0X", "compute_V0X", iter_timing); } if (r_params.print_level() >= 20) { print_inner(world, "xV0x", chi, V0X); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } X_space E0X(world, chi.num_states(), chi.num_orbitals()); @@ -544,16 +537,10 @@ auto ResponseBase::compute_theta_X(World &world, const X_space &chi, E0X.y = E0X.x.copy(); E0X.x.truncate_rf(); } - inner_to_json(world, "E0", response_context.inner(chi, E0X), - iter_function_data); - if (r_params.print_level() >= 20) { - print_inner(world, "xE0x", chi, E0X); - } - } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "compute_E0X", "compute_E0X", - iter_timing); + inner_to_json(world, "E0", response_context.inner(chi, E0X), iter_function_data); + if (r_params.print_level() >= 20) { print_inner(world, "xE0x", chi, E0X); } } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "compute_E0X", "compute_E0X", iter_timing); } X_space gamma; if (r_params.print_level() >= 1) { molresponse::start_timer(world); } if (calc_type != "tda") { @@ -563,46 +550,31 @@ auto ResponseBase::compute_theta_X(World &world, const X_space &chi, } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_compute", "gamma_compute", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma_compute", "gamma_compute", iter_timing); } - inner_to_json(world, "gamma_x", response_context.inner(chi, gamma), - iter_function_data); + inner_to_json(world, "gamma_x", response_context.inner(chi, gamma), iter_function_data); if (r_params.print_level() >= 1) { molresponse::start_timer(world); } Theta_X = (V0X - E0X) + gamma; - inner_to_json(world, "theta_x", response_context.inner(chi, Theta_X), - iter_function_data); + inner_to_json(world, "theta_x", response_context.inner(chi, Theta_X), iter_function_data); world.gop.fence(); // Theta_X.truncate(); if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "compute_ThetaX_add", - "compute_ThetaX_add", iter_timing); - } - if (r_params.print_level() >= 20) { - print_inner(world, "xThetax", chi, Theta_X); - } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "compute_ThetaX", "compute_ThetaX", - iter_timing); + molresponse::end_timer(world, "compute_ThetaX_add", "compute_ThetaX_add", iter_timing); } + if (r_params.print_level() >= 20) { print_inner(world, "xThetax", chi, Theta_X); } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "compute_ThetaX", "compute_ThetaX", iter_timing); } return Theta_X; } -auto ResponseBase::compute_gamma_full(World &world, - const gamma_orbitals &density, - const XCOperator &xc) const - -> X_space { - std::shared_ptr>> old_pmap = - FunctionDefaults<3>::get_pmap(); +auto ResponseBase::compute_gamma_full(World &world, const gamma_orbitals &density, + const XCOperator &xc) const -> X_space { + std::shared_ptr>> old_pmap = FunctionDefaults<3>::get_pmap(); - auto [chi_alpha, phi0, rho1] = - orbital_load_balance(world, density, r_params.loadbalparts()); + auto [chi_alpha, phi0, rho1] = orbital_load_balance(world, density, r_params.loadbalparts()); QProjector projector(world, phi0); size_t num_states = chi_alpha.num_states(); @@ -617,46 +589,32 @@ auto ResponseBase::compute_gamma_full(World &world, X_space W = X_space::zero_functions(world, num_states, num_orbitals); if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_zero_functions", - "gamma_zero_functions", iter_timing); + molresponse::end_timer(world, "gamma_zero_functions", "gamma_zero_functions", iter_timing); } auto apply_projector = [&](auto &xi) { return projector(xi); }; // apply the exchange kernel to rho if necessary if (r_params.print_level() >= 1) { molresponse::start_timer(world); } auto rho_b = make_density(world, chi_alpha); - auto J = response_context.compute_j1(world, chi_alpha, rho_b, phi0, - shared_coulomb_operator); + auto J = response_context.compute_j1(world, chi_alpha, rho_b, phi0, shared_coulomb_operator); world.gop.fence(); - inner_to_json(world, "j1", response_context.inner(chi_alpha, J), - iter_function_data); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); - } + inner_to_json(world, "j1", response_context.inner(chi_alpha, J), iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); } // Create Coulomb potential on ground_orbitals if (xcf.hf_exchange_coefficient() != 1.0) { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } W = response_context.compute_VXC1(world, chi_alpha, rho_b, phi0, xc); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "XC[omega]", "XC[omega]", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "XC[omega]", "XC[omega]", iter_timing); } } - inner_to_json(world, "v1_xc", response_context.inner(chi_alpha, W), - iter_function_data); + inner_to_json(world, "v1_xc", response_context.inner(chi_alpha, W), iter_function_data); if (r_params.print_level() >= 1) { molresponse::start_timer(world); } auto K = response_context.compute_k1(world, chi_alpha, phi0); - inner_to_json(world, "k1", response_context.inner(chi_alpha, K), - iter_function_data); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); - } - if (r_params.print_level() >= 20) { - print_inner(world, "old xKx", chi_alpha, K); - } + inner_to_json(world, "k1", response_context.inner(chi_alpha, K), iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); } + if (r_params.print_level() >= 20) { print_inner(world, "old xKx", chi_alpha, K); } molresponse::start_timer(world); X_space gamma(world, num_states, num_orbitals); auto c_xc = xcf.hf_exchange_coefficient(); @@ -667,19 +625,13 @@ auto ResponseBase::compute_gamma_full(World &world, if (world.rank() == 0) { print("gamma: += W"); } } //gamma.truncate(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma add", "gamma_truncate_add", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma add", "gamma_truncate_add", iter_timing); } // project out ground state if (r_params.print_level() >= 1) { molresponse::start_timer(world); } gamma = oop_apply(gamma, apply_projector); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_project", "gamma_project", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma_project", "gamma_project", iter_timing); } if (r_params.print_level() >= 20) { molresponse::start_timer(world); print_inner(world, "xJx", chi_alpha, J); @@ -704,14 +656,11 @@ auto ResponseBase::compute_gamma_full(World &world, } -auto ResponseBase::compute_gamma_static(World &world, - const gamma_orbitals &gammaOrbitals, - const XCOperator &xc) const - -> X_space { +auto ResponseBase::compute_gamma_static(World &world, const gamma_orbitals &gammaOrbitals, + const XCOperator &xc) const -> X_space { auto old_pmap = FunctionDefaults<3>::get_pmap(); - auto [xy, phi0, rho1] = - orbital_load_balance(world, gammaOrbitals, r_params.loadbalparts()); + auto [xy, phi0, rho1] = orbital_load_balance(world, gammaOrbitals, r_params.loadbalparts()); QProjector projector(world, phi0); size_t num_states = xy.num_states(); size_t num_orbitals = xy.num_orbitals(); @@ -727,30 +676,23 @@ auto ResponseBase::compute_gamma_static(World &world, X_space KX = X_space::zero_functions(world, num_states, num_orbitals); X_space KY = X_space::zero_functions(world, num_states, num_orbitals); if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_zero_functions", - "gamma_zero_functions", iter_timing); + molresponse::end_timer(world, "gamma_zero_functions", "gamma_zero_functions", iter_timing); } // apply the exchange kernel to rho if necessary if (r_params.print_level() >= 1) { molresponse::start_timer(world); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } auto rho = make_density(world, xy); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "compute density J[omega]"); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "compute density J[omega]"); } int b = 0; for (const auto &rho_b: rho) { auto temp_J = apply(*shared_coulomb_operator, rho_b); J.x[b++] = mul(world, temp_J, phi0); } - std::transform(J.x.begin(), J.x.end(), J.x.begin(), - [&](auto &jxi) { return projector(jxi); }); + std::transform(J.x.begin(), J.x.end(), J.x.begin(), [&](auto &jxi) { return projector(jxi); }); J.y = J.x.copy(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); - } - inner_to_json(world, "j1", response_context.inner(xy, J), - iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); } + inner_to_json(world, "j1", response_context.inner(xy, J), iter_function_data); if (xcf.hf_exchange_coefficient() != 1.0) { auto compute_wx = [&, &phi0 = phi0](auto rho_alpha) { auto xc_rho = xc.apply_xc_kernel(rho_alpha); @@ -758,47 +700,30 @@ auto ResponseBase::compute_gamma_static(World &world, }; if (r_params.print_level() >= 1) { molresponse::start_timer(world); } std::transform(rho.begin(), rho.end(), W.x.begin(), compute_wx); - std::transform(W.x.begin(), W.x.end(), W.x.begin(), - [&](auto &wxi) { return projector(wxi); }); + std::transform(W.x.begin(), W.x.end(), W.x.begin(), [&](auto &wxi) { return projector(wxi); }); W.y = W.x.copy(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "XC[omega]", "XC[omega]", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "XC[omega]", "XC[omega]", iter_timing); } } - inner_to_json(world, "v1_xc", response_context.inner(xy, W), - iter_function_data); + inner_to_json(world, "v1_xc", response_context.inner(xy, W), iter_function_data); if (r_params.print_level() >= 1) { molresponse::start_timer(world); } K = response_exchange_multiworld(phi0, xy, false); - std::transform(K.x.begin(), K.x.end(), K.x.begin(), - [&](auto &kxi) { return projector(kxi); }); - inner_to_json(world, "k1", response_context.inner(xy, K), - iter_function_data); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); - } - if (r_params.print_level() >= 20) { - print_inner(world, "new static KX", xy, K); - } + std::transform(K.x.begin(), K.x.end(), K.x.begin(), [&](auto &kxi) { return projector(kxi); }); + inner_to_json(world, "k1", response_context.inner(xy, K), iter_function_data); + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); } + if (r_params.print_level() >= 20) { print_inner(world, "new static KX", xy, K); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } gamma = 2 * J - K * xcf.hf_exchange_coefficient() + W; - inner_to_json(world, "gamma_x_beforeQ", response_context.inner(xy, gamma), - iter_function_data); + inner_to_json(world, "gamma_x_beforeQ", response_context.inner(xy, gamma), iter_function_data); if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_truncate_add", - "gamma_truncate_add", iter_timing); + molresponse::end_timer(world, "gamma_truncate_add", "gamma_truncate_add", iter_timing); } // project out ground state if (r_params.print_level() >= 1) { molresponse::start_timer(world); } //for (size_t i = 0; i < num_states; i++) { gamma.X[i] = projector(gamma.X[i]); } - std::transform(gamma.x.begin(), gamma.x.end(), gamma.x.begin(), - [&](auto &gxi) { return projector(gxi); }); + std::transform(gamma.x.begin(), gamma.x.end(), gamma.x.begin(), [&](auto &gxi) { return projector(gxi); }); gamma.y = gamma.x.copy(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_project", "gamma_project", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma_project", "gamma_project", iter_timing); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } J.clear(); @@ -812,22 +737,17 @@ auto ResponseBase::compute_gamma_static(World &world, } if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_clear_functions", - "gamma_clear_functions", iter_timing); + molresponse::end_timer(world, "gamma_clear_functions", "gamma_clear_functions", iter_timing); } // Done return gamma; // Get sizes } -auto ResponseBase::compute_gamma_tda(World &world, - const gamma_orbitals &density, - const XCOperator &xc) const +auto ResponseBase::compute_gamma_tda(World &world, const gamma_orbitals &density, const XCOperator &xc) const -> X_space { - auto [d_alpha, phi0, rho1] = - orbital_load_balance(world, density, r_params.loadbalparts()); - std::shared_ptr>> oldpmap = - FunctionDefaults<3>::get_pmap(); + auto [d_alpha, phi0, rho1] = orbital_load_balance(world, density, r_params.loadbalparts()); + std::shared_ptr>> oldpmap = FunctionDefaults<3>::get_pmap(); size_t num_states = d_alpha.num_states(); size_t num_orbitals = d_alpha.num_orbitals(); @@ -840,8 +760,7 @@ auto ResponseBase::compute_gamma_tda(World &world, response_space k1_x(world, num_states, num_orbitals); response_space W(world, num_states, num_orbitals); if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_zero_functions", - "gamma_zero_functions", iter_timing); + molresponse::end_timer(world, "gamma_zero_functions", "gamma_zero_functions", iter_timing); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } @@ -855,9 +774,7 @@ auto ResponseBase::compute_gamma_tda(World &world, }; std::transform(rho.begin(), rho.end(), J.begin(), compute_jx); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "J[omega]", "J[omega]", iter_timing); } // Create Coulomb potential on ground_orbitals if (xcf.hf_exchange_coefficient() != 1.0) { @@ -871,10 +788,7 @@ auto ResponseBase::compute_gamma_tda(World &world, std::transform(rho.begin(), rho.end(), W.begin(), compute_wx); W = W.copy(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "XC[omega]", "XC[omega]", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "XC[omega]", "XC[omega]", iter_timing); } } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } @@ -885,9 +799,7 @@ auto ResponseBase::compute_gamma_tda(World &world, k1_x[b] = newK(x, phi0, phi0); } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "K[omega]", "K[omega]", iter_timing); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } k1_x.truncate_rf(); @@ -896,8 +808,7 @@ auto ResponseBase::compute_gamma_tda(World &world, gamma.x = (J * 2) - k1_x * xcf.hf_exchange_coefficient() + W; if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_truncate_add", - "gamma_truncate_add", iter_timing); + molresponse::end_timer(world, "gamma_truncate_add", "gamma_truncate_add", iter_timing); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } @@ -906,10 +817,7 @@ auto ResponseBase::compute_gamma_tda(World &world, gamma.x[i] = projector(gamma.x[i]); truncate(world, gamma.x[i]); } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_project", "gamma_project", - iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "gamma_project", "gamma_project", iter_timing); } if (r_params.print_level() >= 20) { print("------------------------ Gamma Functions Norms " @@ -932,18 +840,15 @@ auto ResponseBase::compute_gamma_tda(World &world, } if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "gamma_clear_functions", - "gamma_clear_functions", iter_timing); + molresponse::end_timer(world, "gamma_clear_functions", "gamma_clear_functions", iter_timing); } // Done world.gop.fence(); return gamma; } -auto ResponseBase::compute_lambda_X(World &world, const X_space &chi, - XCOperator &xc, - const std::string &calc_type) const - -> X_space { +auto ResponseBase::compute_lambda_X(World &world, const X_space &chi, XCOperator &xc, + const std::string &calc_type) const -> X_space { // compute bool compute_Y = calc_type == "full"; @@ -969,14 +874,11 @@ auto ResponseBase::compute_lambda_X(World &world, const X_space &chi, X_space gamma; // compute if (calc_type == "full") { - gamma = compute_gamma_full( - world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); + gamma = compute_gamma_full(world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); } else if (calc_type == "static") { - gamma = compute_gamma_static( - world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); + gamma = compute_gamma_static(world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); } else { - gamma = compute_gamma_tda( - world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); + gamma = compute_gamma_tda(world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); } if (r_params.print_level() >= 20) { auto gamma_mx = inner(Chi_truncated, gamma); @@ -999,9 +901,8 @@ auto ResponseBase::compute_lambda_X(World &world, const X_space &chi, return Lambda_X; } -auto ResponseBase::compute_response_potentials( - World &world, const X_space &chi, XCOperator &xc, - const std::string &calc_type) const +auto ResponseBase::compute_response_potentials(World &world, const X_space &chi, XCOperator &xc, + const std::string &calc_type) const -> std::tuple { // compute bool compute_Y = calc_type == "full"; @@ -1036,19 +937,15 @@ auto ResponseBase::compute_response_potentials( X_space gamma; // compute if (calc_type == "full") { - gamma = compute_gamma_full( - world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); + gamma = compute_gamma_full(world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); } else if (calc_type == "static") { - gamma = compute_gamma_static( - world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); + gamma = compute_gamma_static(world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); } else { - gamma = compute_gamma_tda( - world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); + gamma = compute_gamma_tda(world, {chi, ground_orbitals, vector_real_function_3d{}}, xc); } - X_space Lambda_X( - world, m, - n);// = X_space(world, chi.num_states(), chi.num_orbitals()); + X_space Lambda_X(world, m, + n);// = X_space(world, chi.num_states(), chi.num_orbitals()); Lambda_X = (T0X + V0X - E0X) + gamma; @@ -1061,9 +958,8 @@ auto ResponseBase::compute_response_potentials( // J0=J[ground_density] // K0=K[ground_density]f // EXC0=W[ground_density] -auto ResponseBase::compute_V0X(World &world, const X_space &X, - const XCOperator &xc, - bool compute_Y) const -> X_space { +auto ResponseBase::compute_V0X(World &world, const X_space &X, const XCOperator &xc, bool compute_Y) const + -> X_space { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } size_t m = X.num_states(); @@ -1084,9 +980,7 @@ auto ResponseBase::compute_V0X(World &world, const X_space &X, } else {// Already pre-computed v_nuc = stored_v_nuc; } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "V0_nuc", "V0_nuc", iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "V0_nuc", "V0_nuc", iter_timing); } // Coulomb Potential J0*f if (r_params.print_level() >= 1) { molresponse::start_timer(world); } @@ -1098,21 +992,15 @@ auto ResponseBase::compute_V0X(World &world, const X_space &X, } else {// Already pre-computed v_j0 = stored_v_coul; } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "J[0]", "J[0]", iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "J[0]", "J[0]", iter_timing); } if (xcf.hf_exchange_coefficient() != 1.0) { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } v_xc = xc.make_xc_potential(); - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "XC[0]", "XC[0]", iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "XC[0]", "XC[0]", iter_timing); } } else { // make a zero function - v_xc = Function( - FunctionFactory(world).fence(false).initial_level( - 1)); + v_xc = Function(FunctionFactory(world).fence(false).initial_level(1)); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } K0 = ground_exchange_multiworld(ground_orbitals, X, compute_Y); @@ -1127,18 +1015,12 @@ auto ResponseBase::compute_V0X(World &world, const X_space &X, if (r_params.print_level() >= 15) { - inner_to_json(world, "v0_nuc", response_context.inner(X, v_nuc * X), - iter_function_data); - inner_to_json(world, "j0", response_context.inner(X, v_j0 * X), - iter_function_data); - inner_to_json(world, "k0", response_context.inner(X, K0), - iter_function_data); - inner_to_json(world, "v0_xc", response_context.inner(X, v_xc * X), - iter_function_data); - } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "K[0]", "K[0]", iter_timing); + inner_to_json(world, "v0_nuc", response_context.inner(X, v_nuc * X), iter_function_data); + inner_to_json(world, "j0", response_context.inner(X, v_j0 * X), iter_function_data); + inner_to_json(world, "k0", response_context.inner(X, K0), iter_function_data); + inner_to_json(world, "v0_xc", response_context.inner(X, v_xc * X), iter_function_data); } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "K[0]", "K[0]", iter_timing); } if (r_params.print_level() >= 1) { molresponse::start_timer(world); } @@ -1150,26 +1032,21 @@ auto ResponseBase::compute_V0X(World &world, const X_space &X, V0 = X * v0; V0 += -c_xc * K0; V0.truncate(); - inner_to_json(world, "v0", response_context.inner(X, V0), - iter_function_data); + inner_to_json(world, "v0", response_context.inner(X, V0), iter_function_data); } else { V0 = X.copy(); V0.x = v0 * V0.x; V0.x += -c_xc * K0.x; V0.x.truncate_rf(); V0.y = V0.x.copy(); - inner_to_json(world, "v0", response_context.inner(X, V0), - iter_function_data); + inner_to_json(world, "v0", response_context.inner(X, V0), iter_function_data); } if (r_params.print_level() >= 20) { print_inner(world, "xV0x", X, V0); } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "V0_add", "V0_add", iter_timing); - } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "V0_add", "V0_add", iter_timing); } return V0; } -auto ResponseBase::compute_TX(World &world, const X_space &x, - bool compute_Y) const -> X_space { +auto ResponseBase::compute_TX(World &world, const X_space &x, bool compute_Y) const -> X_space { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } X_space T0X = X_space(world, x.num_states(), x.num_orbitals()); @@ -1219,9 +1096,8 @@ auto ResponseBase::compute_TX(World &world, const X_space &x, return T0X; } // Returns the ground state fock operator applied to functions f -auto ResponseBase::compute_F0X(World &world, const X_space &X, - const XCOperator &xc, - bool compute_Y) const -> X_space { +auto ResponseBase::compute_F0X(World &world, const X_space &X, const XCOperator &xc, bool compute_Y) const + -> X_space { // Debugging output molresponse::start_timer(world); @@ -1271,11 +1147,8 @@ auto ResponseBase::compute_F0X(World &world, const X_space &X, return F0X; } -auto ResponseBase::update_residual(World &world, const X_space &chi, - const X_space &g_chi, - const std::string &calc_type, - const Tensor &old_residuals, - const X_space &xres_old) -> residuals { +auto ResponseBase::update_residual(World &world, const X_space &chi, const X_space &g_chi, const std::string &calc_type, + const Tensor &old_residuals, const X_space &xres_old) -> residuals { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } size_t m = chi.x.size(); size_t n = chi.x.size_orbitals(); @@ -1288,28 +1161,20 @@ auto ResponseBase::update_residual(World &world, const X_space &chi, res = chi - g_chi; auto rx = to_response_matrix(res); auto gx = to_response_matrix(g_chi); - for (const auto &b: chi.active) { - residual_norms(b) = norm2(world, rx[b]) / norm2(world, gx[b]); - } + for (const auto &b: chi.active) { residual_norms(b) = norm2(world, rx[b]) / norm2(world, gx[b]); } } else { res.x = chi.x - g_chi.x; - for (const auto &b: chi.active) { - residual_norms(b) = - norm2(world, res.x[b]) / norm2(world, g_chi.x[b]); - } + for (const auto &b: chi.active) { residual_norms(b) = norm2(world, res.x[b]) / norm2(world, g_chi.x[b]); } } if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "compute_bsh_residual", - "compute_bsh_residual", iter_timing); + molresponse::end_timer(world, "compute_bsh_residual", "compute_bsh_residual", iter_timing); } // Next calculate 2-norm of these vectors of differences return {res, residual_norms}; } -auto ResponseBase::kain_x_space_update(World &world, const X_space &chi, - const X_space &residual_chi, - response_solver &kain_x_space) - -> X_space { +auto ResponseBase::kain_x_space_update(World &world, const X_space &chi, const X_space &residual_chi, + response_solver &kain_x_space) -> X_space { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } size_t m = chi.num_states(); size_t n = chi.num_orbitals(); @@ -1326,22 +1191,14 @@ auto ResponseBase::kain_x_space_update(World &world, const X_space &chi, std::copy(temp.begin() + n, temp.end(), kain_update.y[i].begin()); }; } else { - for (const auto &i: Chi.active) { - kain_update.x[i] = - kain_x_space[i].update(chi.x[i], residual_chi.x[i]); - } - } - if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "kain_x_update", "kain_x_update", - iter_timing); + for (const auto &i: Chi.active) { kain_update.x[i] = kain_x_space[i].update(chi.x[i], residual_chi.x[i]); } } + if (r_params.print_level() >= 1) { molresponse::end_timer(world, "kain_x_update", "kain_x_update", iter_timing); } return kain_update; } -void ResponseBase::x_space_step_restriction(World &world, - const X_space &old_Chi, - X_space &temp, bool restrict_y, +void ResponseBase::x_space_step_restriction(World &world, const X_space &old_Chi, X_space &temp, bool restrict_y, const double &max_bsh_rotation) { size_t m = old_Chi.num_states(); size_t n = old_Chi.num_orbitals(); @@ -1361,21 +1218,17 @@ void ResponseBase::x_space_step_restriction(World &world, auto norm_xb = norm2(world, m_old[b]); auto max_step = max_bsh_rotation;//norm;//* norm_xb; if (world.rank() == 0) { - print("---------------- step restriction :", b, - " ------------------"); + print("---------------- step restriction :", b, " ------------------"); if (world.rank() == 0) { print("X[b]: ", norm_xb); } if (world.rank() == 0) { print("deltaX[b]: ", step_size); } - if (world.rank() == 0) { - print("max_step = max_rotation*norm_X: ", max_step); - } + if (world.rank() == 0) { print("max_step = max_rotation*norm_X: ", max_step); } } if (step_size > max_step && step_size < 25) { // and if the step size is less thant 10% the vector norm double s = .80 * max_step / step_size; if (world.rank() == 0) { if (r_params.print_level() > 1) - print(" restricting step for response-state: ", b, - " step size", s); + print(" restricting step for response-state: ", b, " step size", s); } gaxpy(world, s, m_new[b], (1.0 - s), m_old[b], false); } @@ -1387,107 +1240,29 @@ void ResponseBase::x_space_step_restriction(World &world, auto norm_xb = norm2(world, old_Chi.x[b]); auto max_step = max_bsh_rotation;//norm;//* norm_xb; if (world.rank() == 0) { - print("---------------- step restriction :", b, - " ------------------"); + print("---------------- step restriction :", b, " ------------------"); if (world.rank() == 0) { print("X[b]: ", norm_xb); } if (world.rank() == 0) { print("deltaX[b]: ", step_size); } - if (world.rank() == 0) { - print("max_step = max_rotation*norm_X: ", max_step); - } + if (world.rank() == 0) { print("max_step = max_rotation*norm_X: ", max_step); } } - if (step_size > max_step && - step_size < 10000 * FunctionDefaults<3>::get_thresh()) { + if (step_size > max_step && step_size < 10000 * FunctionDefaults<3>::get_thresh()) { // and if the step size is less thant 10% the vector norm double s = .80 * max_step / step_size; if (world.rank() == 0) { if (r_params.print_level() > 1) - print(" restricting step for response-state: ", b, - " step size", s); + print(" restricting step for response-state: ", b, " step size", s); } gaxpy(world, s, temp.x[b], (1.0 - s), old_Chi.x[b], false); } } } - if (world.rank() == 0) { - print("----------------End Step Restriction -----------------"); - } + if (world.rank() == 0) { print("----------------End Step Restriction -----------------"); } if (r_params.print_level() >= 1) { - molresponse::end_timer(world, "x_space_restriction", - "x_space_restriction", iter_timing); + molresponse::end_timer(world, "x_space_restriction", "x_space_restriction", iter_timing); } } -//void ResponseBase::plotResponseOrbitals( -// World &world, size_t iteration, const response_space &x_response, -// const response_space &y_response, -// ResponseParameters const &responseParameters, -// GroundStateCalculation const &g_params) { -// std::filesystem::create_directories("plots/densities"); -// std::filesystem::create_directory("plots/orbitals"); -// -// // TESTING -// // get transition density -// // num orbitals -// size_t n = x_response[0].size(); -// size_t m = x_response.size(); -// -// real_function_3d rho0 = dot(world, ground_orbitals, ground_orbitals); -// std::vector rho1 = -// transition_density(world, ground_orbitals, x_response, y_response); -// std::string dir("xyz"); -// // for plot_name size -// size_t buffSize = 500; -// char plot_name[buffSize]; -// double Lp = std::min(responseParameters.L(), 24.0); -// // Doing line plots along each axis -// for (int d = 0; d < 3; d++) { -// // print ground_state -// plotCoords plt(d, Lp); -// // plot ground density -// if (iteration == 1) { -// snprintf(plot_name, buffSize, "plots/densities/rho0_%c_0.plt", -// dir[d]); -// plot_line(plot_name, 5001, plt.lo, plt.hi, rho0); -// } -// // plot ground orbitals and transition densities in xyz directions -// snprintf(plot_name, buffSize, "plots/orbitals/phi0_%c_0.plt", dir[d]); -// plot_line(plot_name, 5001, plt.lo, plt.hi, ground_orbitals); -// snprintf(plot_name, buffSize, "plots/densities/rho1_%c.plt", dir[d]); -// plot_line(plot_name, 5001, plt.lo, plt.hi, rho1); -// -// for (int b = 0; b < static_cast(m); b++) { -// // plot x and y transition orbitals -// snprintf(plot_name, buffSize, "plots/orbitals/phix_%c_%d.plt", -// dir[d], static_cast(b)); -// plot_line(plot_name, 5001, plt.lo, plt.hi, x_response[b]); -// snprintf(plot_name, buffSize, "plots/orbitals/phiy_%c_%d.plt", -// dir[d], static_cast(b)); -// plot_line(plot_name, 5001, plt.lo, plt.hi, y_response[b]); -// } -// } -// world.gop.fence(); -// // END TESTING -//} - - -//void PlotGroundDensityVTK(World &world, const ResponseBase &calc) { -// auto [ground_calc, molecule, r_params] = calc.get_parameter(); -// auto ground_orbitals = calc.get_orbitals(); - -// if (r_params.plot_initial()) { -// if (world.rank() == 0) print("\n Plotting ground state densities.\n"); -// if (r_params.plot_l() > 0.0) -// do_vtk_plots(world, int(r_params.plot_pts()), r_params.plot_l(), 0, -// int(r_params.num_orbitals()), molecule, -// square(world, ground_orbitals), "ground"); -// else -// do_vtk_plots(world, int(r_params.plot_pts()), r_params.L() / 2.0, 0, -// int(r_params.num_orbitals()), molecule, -// square(world, ground_orbitals), "ground"); -// } -//} - /// Push back empty json onto "protocol_data" field /// Writes protocol to that new field /// Sets up the iteration data json with an empty {} @@ -1501,10 +1276,8 @@ void protocol_to_json(json &j, const double proto) { j["protocol_data"][proto_index]["iter_data"] = {}; } -void ResponseBase::function_data_to_json(json &j_mol_in, size_t iter, - const Tensor &x_norms, - const Tensor &x_abs_norms, - const Tensor &rho_norms, +void ResponseBase::function_data_to_json(json &j_mol_in, size_t iter, const Tensor &x_norms, + const Tensor &x_abs_norms, const Tensor &rho_norms, const Tensor &rho_abs_norms) { json j = {}; @@ -1540,9 +1313,7 @@ void ResponseBase::solve(World &world) { // protocol if (first_protocol) { if (r_params.restart()) { - if (world.rank() == 0) { - print(" Restarting from file:", r_params.restart_file()); - } + if (world.rank() == 0) { print(" Restarting from file:", r_params.restart_file()); } load(world, r_params.restart_file()); first_protocol = false; } else { @@ -1550,15 +1321,11 @@ void ResponseBase::solve(World &world) { if (world.rank() == 0) { print("Successfully initialized "); } } check_k(world, iter_thresh, FunctionDefaults<3>::get_k()); - if (world.rank() == 0) { - print("Successfully check K first initialization "); - } + if (world.rank() == 0) { print("Successfully check K first initialization "); } first_protocol = false; } else { check_k(world, iter_thresh, FunctionDefaults<3>::get_k()); - if (world.rank() == 0) { - print("Successfully check K not first initialization "); - } + if (world.rank() == 0) { print("Successfully check K not first initialization "); } } protocol_to_json(j_molresponse, iter_thresh); // Now actually ready to iterate... @@ -1568,20 +1335,19 @@ void ResponseBase::solve(World &world) { converged_to_json(j_molresponse); if (r_params.plot()) { auto r_matrix = to_response_matrix(Chi); -// do_response_orbital_vtk_plots(world, r_params.plot_pts(), r_params.L(), -// molecule, ground_orbitals, r_matrix); + // do_response_orbital_vtk_plots(world, r_params.plot_pts(), r_params.L(), + // molecule, ground_orbitals, r_matrix); auto response_densities = make_density(world, Chi); -// do_response_density_vtk_plots(world, r_params.plot_pts(), r_params.L(), -// molecule, ground_density, -// response_densities); + // do_response_density_vtk_plots(world, r_params.plot_pts(), r_params.L(), + // molecule, ground_density, + // response_densities); } // Plot the response function if desired } -void check_k(World &world, X_space &Chi, - double thresh = FunctionDefaults<3>::get_thresh(), +void check_k(World &world, X_space &Chi, double thresh = FunctionDefaults<3>::get_thresh(), int k = FunctionDefaults<3>::get_k()) { if (0 != Chi.x.size()) { if (FunctionDefaults<3>::get_k() != Chi.x[0].at(0).k()) { @@ -1589,18 +1355,12 @@ void check_k(World &world, X_space &Chi, for (auto &xi: Chi.x) { reconstruct(world, xi); - for (auto &xij: xi) { - xij = project(xij, FunctionDefaults<3>::get_k(), thresh, - false); - } + for (auto &xij: xi) { xij = project(xij, FunctionDefaults<3>::get_k(), thresh, false); } world.gop.fence(); } for (auto &yi: Chi.y) { reconstruct(world, yi); - for (auto &yij: yi) { - yij = project(yij, FunctionDefaults<3>::get_k(), thresh, - false); - } + for (auto &yij: yi) { yij = project(yij, FunctionDefaults<3>::get_k(), thresh, false); } world.gop.fence(); } Chi.truncate(); @@ -1616,8 +1376,7 @@ void check_k(World &world, X_space &Chi, /// \param f /// \param magnitude /// \return -auto add_randomness(World &world, const response_space &f, double magnitude) - -> response_space { +auto add_randomness(World &world, const response_space &f, double magnitude) -> response_space { // Copy input functions response_space f_copy = f.copy(); @@ -1679,20 +1438,15 @@ void normalize(World &world, X_space &Chi) { } } -auto solid_harmonics(World &world, int n) - -> std::map, real_function_3d> { +auto solid_harmonics(World &world, int n) -> std::map, real_function_3d> { // Container to return std::map, real_function_3d> result; // Create the basic x, y, z, constant and zero - real_function_3d x = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{1, 0, 0}))); - real_function_3d y = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{0, 1, 0}))); - real_function_3d z = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{0, 0, 1}))); - real_function_3d c = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{0, 0, 0}))); + real_function_3d x = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{1, 0, 0}))); + real_function_3d y = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{0, 1, 0}))); + real_function_3d z = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{0, 0, 1}))); + real_function_3d c = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{0, 0, 0}))); real_function_3d zero = real_factory_3d(world); // Add in first few, since they're simple @@ -1707,12 +1461,10 @@ auto solid_harmonics(World &world, int n) // Calculate ends of this row first result[std::vector{l + 1, l + 1}] = sqrt(pow(2, kronecker(l, 0) * (2 * l) / (2 * l + 1))) * - (x * result[std::vector{l, l}] - - (1 - kronecker(l, 0) * y * result[std::vector{l, -l}])); + (x * result[std::vector{l, l}] - (1 - kronecker(l, 0) * y * result[std::vector{l, -l}])); result[std::vector{l + 1, -l - 1}] = sqrt(pow(2, kronecker(l, 0) * (2 * l) / (2 * l + 1))) * - (y * result[std::vector{l, l}] + - (1 - kronecker(l, 0) * x * result[std::vector{l, -l}])); + (y * result[std::vector{l, l}] + (1 - kronecker(l, 0) * x * result[std::vector{l, -l}])); // Formula below calls for some functions that don't exist. // Need zeroes where that would occur @@ -1725,8 +1477,7 @@ auto solid_harmonics(World &world, int n) result[std::vector{l + 1, m}] = 1.0 / std::sqrt((l + m + 1) * (l - m + 1)) * ((2 * l + 1) * z * result[std::vector{l, m}] - - sqrt((l + m) * (l - m)) * (x * x + y * y + z * z) * - result[std::vector{l - 1, m}]); + sqrt((l + m) * (l - m)) * (x * x + y * y + z * z) * result[std::vector{l - 1, m}]); } } @@ -1745,9 +1496,8 @@ auto solid_harmonics(World &world, int n) } -vector_real_function_3d -transition_density(World &world, const vector_real_function_3d &orbitals, - const response_space &x, const response_space &y) { +vector_real_function_3d transition_density(World &world, const vector_real_function_3d &orbitals, + const response_space &x, const response_space &y) { // Get sizes // Check sizes and then run the algorithm //size_t m = x.size(); @@ -1759,13 +1509,11 @@ transition_density(World &world, const vector_real_function_3d &orbitals, //xx.truncate_rf(); // yy.truncate_rf(); //truncate(world, phi0); - std::vector densities = - zero_functions(world, x.size(), true); + std::vector densities = zero_functions(world, x.size(), true); // Return container - auto compute_density = [&world, &orbitals](const auto &x_alpha, - const auto &y_alpha) { + auto compute_density = [&world, &orbitals](const auto &x_alpha, const auto &y_alpha) { /* for (const auto &xij: x_alpha) { print("xij !!", xij.max_depth(), " ", (void *) xij.get_impl().get()); @@ -1785,8 +1533,7 @@ transition_density(World &world, const vector_real_function_3d &orbitals, print("phi_i !!", phi_i.max_depth(), " ", (void *) phi_i.get_impl().get()); } */ - std::transform(x.begin(), x.end(), y.begin(), densities.begin(), - compute_density); + std::transform(x.begin(), x.end(), y.begin(), densities.begin(), compute_density); world.gop.fence(); //truncate(world, densities, FunctionDefaults<3>::get_thresh(), true); return densities; @@ -1801,9 +1548,7 @@ transition_density(World &world, const vector_real_function_3d &orbitals, * @param load_balance * @return */ -auto ResponseBase::orbital_load_balance(World &world, - const gamma_orbitals &gammaOrbitals, - const double load_balance) +auto ResponseBase::orbital_load_balance(World &world, const gamma_orbitals &gammaOrbitals, const double load_balance) -> gamma_orbitals { auto X = std::get<0>(gammaOrbitals); auto psi0 = std::get<1>(gammaOrbitals); @@ -1816,22 +1561,14 @@ auto ResponseBase::orbital_load_balance(World &world, molresponse::start_timer(world); LoadBalanceDeux<3> lb(world); - for (const auto &phi0_i: psi0) { - lb.add_tree(phi0_i, lbcost(1.0, 8.0), false); - } + for (const auto &phi0_i: psi0) { lb.add_tree(phi0_i, lbcost(1.0, 8.0), false); } - for (const auto &rho1_i: rho1) { - lb.add_tree(rho1_i, lbcost(1.0, 8.0), false); - } + for (const auto &rho1_i: rho1) { lb.add_tree(rho1_i, lbcost(1.0, 8.0), false); } for (const auto &xi: X.x) { - for (const auto &xij: xi) { - lb.add_tree(xij, lbcost(1.0, 8.0), false); - } + for (const auto &xij: xi) { lb.add_tree(xij, lbcost(1.0, 8.0), false); } } for (const auto &yi: X.y) { - for (const auto &yij: yi) { - lb.add_tree(yij, lbcost(1.0, 8.0), false); - } + for (const auto &yij: yi) { lb.add_tree(yij, lbcost(1.0, 8.0), false); } } world.gop.fence(); // newpamap is the new pmap just based on the orbitals @@ -1854,8 +1591,7 @@ auto ResponseBase::orbital_load_balance(World &world, } } -void ResponseBase::analyze_vectors(World &world, const vecfuncT &x, - const std::string &response_state) { +void ResponseBase::analyze_vectors(World &world, const vecfuncT &x, const std::string &response_state) { molresponse::start_timer(world); AtomicBasisSet sto3g("sto-3g"); vecfuncT ao = project_ao_basis(world, sto3g); @@ -1870,14 +1606,10 @@ void ResponseBase::analyze_vectors(World &world, const vecfuncT &x, rsq = inner(world, x, mul_sparse(world, frsq, x, vtol)); for (int axis = 0; axis < 3; ++axis) { // x y z - functionT fdip = - factoryT(world) - .functor(functorT(new madness::DipoleFunctor(axis))) - .initial_level(4); + functionT fdip = factoryT(world).functor(functorT(new madness::DipoleFunctor(axis))).initial_level(4); dip(axis, _) = inner(world, x, mul_sparse(world, fdip, x, vtol)); // - ^2-^2-^2 - for (int i = 0; i < nmo1; ++i) - rsq(i) -= dip(axis, i) * dip(axis, i); + for (int i = 0; i < nmo1; ++i) rsq(i) -= dip(axis, i) * dip(axis, i); } } molresponse::end_timer(world, "Analyze vectors"); @@ -1892,26 +1624,18 @@ void ResponseBase::analyze_vectors(World &world, const vecfuncT &x, printf("ncoeff=%.2e:", (double) ncoeffi); - printf("center=(%.2f,%.2f,%.2f) : radius=%.2f\n", dip(0, i), - dip(1, i), dip(2, i), sqrt(rsq(i))); + printf("center=(%.2f,%.2f,%.2f) : radius=%.2f\n", dip(0, i), dip(1, i), dip(2, i), sqrt(rsq(i))); sto3g.print_anal(molecule, C(i, _)); printf("total number of coefficients = %.8e\n\n", double(ncoeff)); } } } -auto ResponseBase::project_ao_basis_only(World &world, - const AtomicBasisSet &aobasis, - const Molecule &mol) -> vecfuncT { +auto ResponseBase::project_ao_basis_only(World &world, const AtomicBasisSet &aobasis, const Molecule &mol) -> vecfuncT { vecfuncT ao = vecfuncT(aobasis.nbf(mol)); for (int i = 0; i < aobasis.nbf(mol); ++i) { - functorT aofunc(new madchem::AtomicBasisFunctor( - aobasis.get_atomic_basis_function(mol, i))); - ao[i] = factoryT(world) - .functor(aofunc) - .truncate_on_project() - .nofence() - .truncate_mode(1); + functorT aofunc(new madchem::AtomicBasisFunctor(aobasis.get_atomic_basis_function(mol, i))); + ao[i] = factoryT(world).functor(aofunc).truncate_on_project().nofence().truncate_mode(1); } world.gop.fence(); truncate(world, ao); @@ -1919,8 +1643,7 @@ auto ResponseBase::project_ao_basis_only(World &world, return ao; } -auto ResponseBase::project_ao_basis(World &world, const AtomicBasisSet &aobasis) - -> vecfuncT { +auto ResponseBase::project_ao_basis(World &world, const AtomicBasisSet &aobasis) -> vecfuncT { // Make at_to_bf, at_nbf ... map from atom to first bf on atom, and nbf/atom std::vector at_to_bf, at_nbf; aobasis.atoms_to_bfn(molecule, at_to_bf, at_nbf); @@ -1953,8 +1676,7 @@ void ResponseBase::output_json() { void ResponseBase::converged_to_json(json &j) { j["converged"] = all_done; } -void ResponseBase::print_inner(World &world, const std::string &name, - const X_space &left, const X_space &right) { +void ResponseBase::print_inner(World &world, const std::string &name, const X_space &left, const X_space &right) { auto m_val = inner(left, right); world.gop.fence(); if (world.rank() == 0) { @@ -1963,10 +1685,8 @@ void ResponseBase::print_inner(World &world, const std::string &name, } } - -auto transition_densityTDA(World &world, - const vector_real_function_3d &orbitals, - const response_space &x) -> vector_real_function_3d { +auto transition_densityTDA(World &world, const vector_real_function_3d &orbitals, const response_space &x) + -> vector_real_function_3d { // Get sizes size_t m = x.size(); @@ -1976,8 +1696,7 @@ auto transition_densityTDA(World &world, xx.truncate_rf(); truncate(world, phi0); - std::vector densities = - zero_functions(world, m); + std::vector densities = zero_functions(world, m); // dot xi with phi0 auto f = [&world, &phi0](auto xi) { return dot(world, xi, phi0); }; @@ -1992,20 +1711,16 @@ auto transition_densityTDA(World &world, // Transforms the given matrix of functions according to the give // transformation matrix. Used to update orbitals / potential -response_space transform(World &world, const response_space &f, - const Tensor &U) { +response_space transform(World &world, const response_space &f, const Tensor &U) { // Return container response_space result; // Go element by element for (unsigned int i = 0; i < f.size(); i++) { // Temp for the result of one row - std::vector temp = - zero_functions_compressed(world, f[0].size()); + std::vector temp = zero_functions_compressed(world, f[0].size()); - for (unsigned int j = 0; j < f.size(); j++) { - gaxpy(world, 1.0, temp, U(j, i), f[j]); - } + for (unsigned int j = 0; j < f.size(); j++) { gaxpy(world, 1.0, temp, U(j, i), f[j]); } // Add to temp to result result.push_back(temp); @@ -2017,8 +1732,7 @@ response_space transform(World &world, const response_space &f, return result; } -auto transform(World &world, const X_space &x, const Tensor &U) - -> X_space { +auto transform(World &world, const X_space &x, const Tensor &U) -> X_space { // Return container X_space result(world, x.num_states(), x.num_orbitals()); @@ -2028,8 +1742,7 @@ auto transform(World &world, const X_space &x, const Tensor &U) return result; } -auto expectation(World &world, const response_space &A, const response_space &B) - -> Tensor { +auto expectation(World &world, const response_space &A, const response_space &B) -> Tensor { // Get sizes MADNESS_ASSERT(!A[0].empty()); MADNESS_ASSERT(A[0].size() == B[0].size()); @@ -2054,9 +1767,7 @@ auto expectation(World &world, const response_space &A, const response_space &B) */ // Run over dimension two // each vector in orbital has dim_1 response functoins associated - for (size_t p = 0; p < dim_2; p++) { - result += matrix_inner(world, A_t[p], B_t[p]); - } + for (size_t p = 0; p < dim_2; p++) { result += matrix_inner(world, A_t[p], B_t[p]); } // Done return result; @@ -2076,17 +1787,14 @@ void print_norms(World &world, const response_space &f) { if (world.rank() == 0) print(norms); } -response_space select_functions(World &world, response_space f, - Tensor &energies, size_t k, +response_space select_functions(World &world, response_space f, Tensor &energies, size_t k, size_t print_level) { // Container for result response_space answer; // Debugging output if (print_level >= 1) { - if (world.rank() == 0) - print("\n Selecting the", k, - "lowest excitation energy components.\n"); + if (world.rank() == 0) print("\n Selecting the", k, "lowest excitation energy components.\n"); } // Get rid of extra functions and save @@ -2101,8 +1809,7 @@ response_space select_functions(World &world, response_space f, // Basic output if (print_level >= 1) { - if (world.rank() == 0) - print(" The selected components have excitation energies:"); + if (world.rank() == 0) print(" The selected components have excitation energies:"); if (world.rank() == 0) print(energies); } @@ -2208,12 +1915,9 @@ auto make_xyz_functions(World &world) -> vector_real_function_3d { // Container to return // Create the basic x, y, z, constant and zero - real_function_3d x = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{1, 0, 0}))); - real_function_3d y = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{0, 1, 0}))); - real_function_3d z = real_factory_3d(world).functor( - real_functor_3d(new MomentFunctor(std::vector{0, 0, 1}))); + real_function_3d x = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{1, 0, 0}))); + real_function_3d y = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{0, 1, 0}))); + real_function_3d z = real_factory_3d(world).functor(real_functor_3d(new MomentFunctor(std::vector{0, 0, 1}))); std::vector funcs = {x, y, z}; return funcs; @@ -2252,8 +1956,7 @@ response_timing::response_timing() : iter(0) { wall_time_data.insert({"E0X", std::vector(0)}); wall_time_data.insert({"E0mDX", std::vector(0)}); wall_time_data.insert({"subspace_reduce", std::vector(0)}); - wall_time_data.insert( - {"diagonalize_response_matrix", std::vector(0)}); + wall_time_data.insert({"diagonalize_response_matrix", std::vector(0)}); cpu_time_data.insert({"iter_total", std::vector(0)}); cpu_time_data.insert({"update", std::vector(0)}); @@ -2285,8 +1988,7 @@ response_timing::response_timing() : iter(0) { cpu_time_data.insert({"E0X", std::vector(0)}); cpu_time_data.insert({"E0mDX", std::vector(0)}); cpu_time_data.insert({"subspace_reduce", std::vector(0)}); - cpu_time_data.insert( - {"diagonalize_response_matrix", std::vector(0)}); + cpu_time_data.insert({"diagonalize_response_matrix", std::vector(0)}); } void response_timing::print_data() { @@ -2301,23 +2003,17 @@ void response_timing::print_data() { * values.second=cpu_time * @param values */ -void response_timing::add_data( - std::map> values) { +void response_timing::add_data(std::map> values) { // print("ADDING DATA"); iter++; - std::for_each( - wall_time_data.begin(), wall_time_data.end(), [&values](auto &v) { - // print(v.first, " : ", values[v.first]); - v.second.push_back( - values[v.first] - .first);// .first to get first value of pair wall_time - }); + std::for_each(wall_time_data.begin(), wall_time_data.end(), [&values](auto &v) { + // print(v.first, " : ", values[v.first]); + v.second.push_back(values[v.first].first);// .first to get first value of pair wall_time + }); std::for_each(cpu_time_data.begin(), cpu_time_data.end(), [&values](auto &v) { //print(v.first, " : ", values[v.first]); - v.second.push_back( - values[v.first] - .second);// .first to get first value of pair wall_time + v.second.push_back(values[v.first].second);// .first to get first value of pair wall_time }); } @@ -2327,13 +2023,9 @@ void response_timing::to_json(json &j) { j["time_data"] = json(); j["time_data"]["iterations"] = iter; j["time_data"]["wall_time"] = json(); - for (const auto &e: wall_time_data) { - j["time_data"]["wall_time"][e.first] = e.second; - } + for (const auto &e: wall_time_data) { j["time_data"]["wall_time"][e.first] = e.second; } j["time_data"]["cpu_time"] = json(); - for (const auto &e: cpu_time_data) { - j["time_data"]["cpu_time"][e.first] = e.second; - } + for (const auto &e: cpu_time_data) { j["time_data"]["cpu_time"][e.first] = e.second; } } /** * add the pair of s wall_time and cpu_time to the time_data and wall_data maps @@ -2346,13 +2038,10 @@ void response_data::add_data(std::map> values) { // print("ADDING DATA"); iter++; std::for_each(function_data.begin(), function_data.end(), [&values](auto &v) { - v.second.push_back( - values[v.first]);// .first to get first value of pair wall_time + v.second.push_back(values[v.first]);// .first to get first value of pair wall_time }); } -void response_data::add_convergence_targets(double p_thresh, - double p_density_target, - double p_bsh_target) { +void response_data::add_convergence_targets(double p_thresh, double p_density_target, double p_bsh_target) { this->thresh.push_back(p_thresh); this->density_target.push_back(p_density_target); this->bsh_target.push_back(p_bsh_target); @@ -2418,13 +2107,11 @@ response_data::response_data() : iter(0) { function_data.insert({"d", std::vector>(0)}); function_data.insert({"density_norms", std::vector>(0)}); function_data.insert({"r_d", std::vector>(0)}); - function_data.insert( - {"x_relative_residuals", std::vector>(0)}); + function_data.insert({"x_relative_residuals", std::vector>(0)}); } -void inner_to_json(World &world, const std::string &name, - const Tensor &m_val, +void inner_to_json(World &world, const std::string &name, const Tensor &m_val, std::map> &data) { data[name] = m_val; } diff --git a/src/apps/molresponse/ResponseBase.hpp b/src/apps/molresponse/ResponseBase.hpp index 7fa9cc524e7..632edd94a6e 100644 --- a/src/apps/molresponse/ResponseBase.hpp +++ b/src/apps/molresponse/ResponseBase.hpp @@ -26,18 +26,14 @@ using namespace madness; class ComputeDensityStrategy { public: virtual ~ComputeDensityStrategy() = default; - virtual vector_real_function_3d - compute_density(World &world, const X_space &x, - const vector_real_function_3d &phi0, - const vector_real_function_3d &rho1, bool update) const = 0; + virtual vector_real_function_3d compute_density(World &world, const X_space &x, const vector_real_function_3d &phi0, + const vector_real_function_3d &rho1, bool update) const = 0; }; class StaticDensityStrategy : public ComputeDensityStrategy { public: - vector_real_function_3d compute_density(World &world, const X_space &x, - const vector_real_function_3d &phi0, - const vector_real_function_3d &rho1, - bool update) const override { + vector_real_function_3d compute_density(World &world, const X_space &x, const vector_real_function_3d &phi0, + const vector_real_function_3d &rho1, bool update) const override { vector_real_function_3d rho_new; if (update) { @@ -59,10 +55,8 @@ class StaticDensityStrategy : public ComputeDensityStrategy { }; class FullDensityStrategy : public ComputeDensityStrategy { public: - vector_real_function_3d compute_density(World &world, const X_space &x, - const vector_real_function_3d &phi0, - const vector_real_function_3d &rho1, - bool update) const override { + vector_real_function_3d compute_density(World &world, const X_space &x, const vector_real_function_3d &phi0, + const vector_real_function_3d &rho1, bool update) const override { vector_real_function_3d rho_new; if (update) { rho_new = copy(world, rho1); @@ -89,22 +83,17 @@ class FullDensityStrategy : public ComputeDensityStrategy { class VXC1Strategy { public: virtual ~VXC1Strategy() = default; - virtual X_space compute_VXC1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const XCOperator &xc) const = 0; + virtual X_space compute_VXC1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const XCOperator &xc) const = 0; }; class VXC1StrategyStandard : public VXC1Strategy { public: - X_space compute_VXC1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const XCOperator &xc) const override { + X_space compute_VXC1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const XCOperator &xc) const override { - X_space W = X_space::zero_functions(world, x.num_states(), - x.num_orbitals()); + X_space W = X_space::zero_functions(world, x.num_states(), x.num_orbitals()); auto compute_wx = [&, &phi0 = phi0](auto rho_alpha) { auto xc_rho = xc.apply_xc_kernel(rho_alpha); @@ -120,21 +109,16 @@ class VXC1StrategyStandard : public VXC1Strategy { class J1Strategy { public: virtual ~J1Strategy() = default; - virtual X_space compute_J1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const poperatorT &coulomb_ops) const = 0; + virtual X_space compute_J1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const poperatorT &coulomb_ops) const = 0; }; class J1StrategyFull : public J1Strategy { public: - X_space compute_J1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const poperatorT &coulomb_ops) const override { + X_space compute_J1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const poperatorT &coulomb_ops) const override { - X_space J = X_space::zero_functions(world, x.num_states(), - x.num_orbitals()); + X_space J = X_space::zero_functions(world, x.num_states(), x.num_orbitals()); vector_real_function_3d temp_J(3); for (const auto &b: x.active) { temp_J[b] = apply(*coulomb_ops, rho1[b]); @@ -148,13 +132,10 @@ class J1StrategyFull : public J1Strategy { class J1StrategyStable : public J1Strategy { public: - X_space compute_J1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const poperatorT &coulomb_ops) const override { + X_space compute_J1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const poperatorT &coulomb_ops) const override { - X_space J = X_space::zero_functions(world, x.num_states(), - x.num_orbitals()); + X_space J = X_space::zero_functions(world, x.num_states(), x.num_orbitals()); if (world.rank() == 0) { print("J1StrategyStable"); } vector_real_function_3d temp_J(3); for (const auto &b: x.active) { @@ -174,8 +155,7 @@ class J1StrategyStable : public J1Strategy { class K1Strategy { public: virtual ~K1Strategy() = default; - virtual X_space compute_K1(World &world, const X_space &x, - const vector_real_function_3d &phi0) const = 0; + virtual X_space compute_K1(World &world, const X_space &x, const vector_real_function_3d &phi0) const = 0; static auto make_k(const vecfuncT &ket, const vecfuncT &bra) { auto &world = ket[0].world(); @@ -189,11 +169,9 @@ class K1Strategy { class K1StrategyFull : public K1Strategy { public: - X_space compute_K1(World &world, const X_space &x, - const vector_real_function_3d &phi0) const override { + X_space compute_K1(World &world, const X_space &x, const vector_real_function_3d &phi0) const override { - auto K = X_space::zero_functions(world, x.num_states(), - x.num_orbitals()); + auto K = X_space::zero_functions(world, x.num_states(), x.num_orbitals()); vector_real_function_3d k1x, k1y, k2x, k2y; @@ -201,14 +179,10 @@ class K1StrategyFull : public K1Strategy { vector_real_function_3d yb; if (world.rank() == 0) { print("K1StrategyFull"); } - auto k1x_temp = - create_response_matrix(x.num_states(), x.num_orbitals()); - auto k1y_temp = - create_response_matrix(x.num_states(), x.num_orbitals()); - auto k2x_temp = - create_response_matrix(x.num_states(), x.num_orbitals()); - auto k2y_temp = - create_response_matrix(x.num_states(), x.num_orbitals()); + auto k1x_temp = create_response_matrix(x.num_states(), x.num_orbitals()); + auto k1y_temp = create_response_matrix(x.num_states(), x.num_orbitals()); + auto k2x_temp = create_response_matrix(x.num_states(), x.num_orbitals()); + auto k2y_temp = create_response_matrix(x.num_states(), x.num_orbitals()); for (const auto &b: x.active) { @@ -234,10 +208,8 @@ class K1StrategyFull : public K1Strategy { class K1StrategyStatic : public K1Strategy { public: - X_space compute_K1(World &world, const X_space &x, - const vector_real_function_3d &phi0) const override { - X_space K = X_space::zero_functions(world, x.num_states(), - x.num_orbitals()); + X_space compute_K1(World &world, const X_space &x, const vector_real_function_3d &phi0) const override { + X_space K = X_space::zero_functions(world, x.num_states(), x.num_orbitals()); vector_real_function_3d k1x, k1y, k2x, k2y; const double lo = 1e-10; @@ -261,8 +233,7 @@ class inner_strategy { public: virtual ~inner_strategy() = default; - [[nodiscard]] virtual Tensor - compute_inner(const X_space &x, const X_space &y) const = 0; + [[nodiscard]] virtual Tensor compute_inner(const X_space &x, const X_space &y) const = 0; }; class Context { @@ -275,23 +246,16 @@ class Context { std::unique_ptr density_strategy_; public: - explicit Context( - std::unique_ptr &&innerStrategy = {}, - std::unique_ptr &&j1Strategy = {}, - std::unique_ptr &&k1Strategy = {}, - std::unique_ptr &&vxc1trategy = {}, - std::unique_ptr &&densityStrategy = {}) - : inner_strategy_(std::move(innerStrategy)), - j1_strategy_(std::move(j1Strategy)), - k1_strategy_(std::move(k1Strategy)), - vxc1_strategy_(std::move(vxc1trategy)), + explicit Context(std::unique_ptr &&innerStrategy = {}, + std::unique_ptr &&j1Strategy = {}, std::unique_ptr &&k1Strategy = {}, + std::unique_ptr &&vxc1trategy = {}, + std::unique_ptr &&densityStrategy = {}) + : inner_strategy_(std::move(innerStrategy)), j1_strategy_(std::move(j1Strategy)), + k1_strategy_(std::move(k1Strategy)), vxc1_strategy_(std::move(vxc1trategy)), density_strategy_(std::move(densityStrategy)) {} - void - set_strategy(std::unique_ptr &&strategy, - std::unique_ptr &&j1Strategy, - std::unique_ptr &&K1Strategy, - std::unique_ptr &&vxc1Strategy, - std::unique_ptr &&densityStrategy) { + void set_strategy(std::unique_ptr &&strategy, std::unique_ptr &&j1Strategy, + std::unique_ptr &&K1Strategy, std::unique_ptr &&vxc1Strategy, + std::unique_ptr &&densityStrategy) { inner_strategy_ = std::move(strategy); j1_strategy_ = std::move(j1Strategy); k1_strategy_ = std::move(K1Strategy); @@ -299,65 +263,50 @@ class Context { density_strategy_ = std::move(densityStrategy); } - [[nodiscard]] Tensor inner(const X_space &x, - const X_space &y) const { + [[nodiscard]] Tensor inner(const X_space &x, const X_space &y) const { if (inner_strategy_) { return inner_strategy_->compute_inner(x, y); } else { - throw madness::MadnessException("Inner product Strategy isn't set", - "Need to set a strategy", 2, 455, + throw madness::MadnessException("Inner product Strategy isn't set", "Need to set a strategy", 2, 455, "inner", "ResponseBase.hpp"); } } - X_space compute_j1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const poperatorT &coulomb_ops) const { + X_space compute_j1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const poperatorT &coulomb_ops) const { if (j1_strategy_) { return j1_strategy_->compute_J1(world, x, rho1, phi0, coulomb_ops); } else { - throw madness::MadnessException("Compute J1 Strategy isn't set", - "Need to set a strategy", 2, 455, - "inner", "ResponseBase.hpp"); + throw madness::MadnessException("Compute J1 Strategy isn't set", "Need to set a strategy", 2, 455, "inner", + "ResponseBase.hpp"); } } - X_space compute_k1(World &world, const X_space &x, - const vector_real_function_3d &phi0) const { + X_space compute_k1(World &world, const X_space &x, const vector_real_function_3d &phi0) const { if (k1_strategy_) { return k1_strategy_->compute_K1(world, x, phi0); } else { - throw madness::MadnessException("Compute K1 Strategy isn't set", - "Need to set a strategy", 2, 455, - "inner", "ResponseBase.hpp"); + throw madness::MadnessException("Compute K1 Strategy isn't set", "Need to set a strategy", 2, 455, "inner", + "ResponseBase.hpp"); } } - X_space compute_VXC1(World &world, const X_space &x, - const vector_real_function_3d &rho1, - const vector_real_function_3d &phi0, - const XCOperator &xc) const { + X_space compute_VXC1(World &world, const X_space &x, const vector_real_function_3d &rho1, + const vector_real_function_3d &phi0, const XCOperator &xc) const { if (vxc1_strategy_) { return vxc1_strategy_->compute_VXC1(world, x, rho1, phi0, xc); } else { - throw madness::MadnessException("Compute VXC1 Strategy isn't set", - "Need to set a strategy", 2, 455, + throw madness::MadnessException("Compute VXC1 Strategy isn't set", "Need to set a strategy", 2, 455, "inner", "ResponseBase.hpp"); } } - vector_real_function_3d compute_density(World &world, const X_space &x, - const vector_real_function_3d &phi0, - const vector_real_function_3d &rho1, - bool update) const { + vector_real_function_3d compute_density(World &world, const X_space &x, const vector_real_function_3d &phi0, + const vector_real_function_3d &rho1, bool update) const { if (density_strategy_) { - return density_strategy_->compute_density(world, x, phi0, rho1, - update); + return density_strategy_->compute_density(world, x, phi0, rho1, update); } else { - throw madness::MadnessException( - "Compute Density Strategy isn't set", - "Need to set a strategy", 2, 455, "inner", - "ResponseBase.hpp"); + throw madness::MadnessException("Compute Density Strategy isn't set", "Need to set a strategy", 2, 455, + "inner", "ResponseBase.hpp"); } } }; @@ -365,25 +314,19 @@ class Context { class full_inner_product : public inner_strategy { public: - [[nodiscard]] Tensor - compute_inner(const X_space &x, const X_space &y) const override { + [[nodiscard]] Tensor compute_inner(const X_space &x, const X_space &y) const override { return inner(x, y); } }; class static_inner_product : public inner_strategy { public: - [[nodiscard]] Tensor - compute_inner(const X_space &x, const X_space &y) const override { + [[nodiscard]] Tensor compute_inner(const X_space &x, const X_space &y) const override { return response_space_inner(x.x, y.x); } }; -typedef std::vector> - response_solver; -typedef std::vector< - XNonlinearSolver> - response_function_solver; +typedef std::vector> response_solver; +typedef std::vector> response_function_solver; class response_timing { std::map> wall_time_data; @@ -412,8 +355,7 @@ class response_data { void to_json(json &j); void add_data(std::map> values); - void add_convergence_targets(double p_thresh, double p_density_target, - double p_bsh_target); + void add_convergence_targets(double p_thresh, double p_density_target, double p_bsh_target); }; class ResponseTester; @@ -424,8 +366,7 @@ struct residuals { }; -using gamma_orbitals = - std::tuple; +using gamma_orbitals = std::tuple; class ResponseBase { public: @@ -440,13 +381,9 @@ class ResponseBase { virtual void iterate(World &world) = 0; //virtual void iterate(); - auto get_parameter() const -> CalcParams { - return {ground_calc, molecule, r_params}; - } + auto get_parameter() const -> CalcParams { return {ground_calc, molecule, r_params}; } - auto get_orbitals() const -> vector_real_function_3d { - return ground_orbitals; - } + auto get_orbitals() const -> vector_real_function_3d { return ground_orbitals; } auto get_chi() const -> X_space { return Chi.copy(); }; @@ -507,14 +444,11 @@ class ResponseBase { std::shared_ptr potential_manager; // shared pointers to Operators - poperatorT - shared_coulomb_operator;// shared pointer to seperated convolution operator + poperatorT shared_coulomb_operator;// shared pointer to seperated convolution operator std::vector> gradop; // Stored functions - mutable real_function_3d - stored_v_nuc;// Stored nuclear potential from ground state - mutable real_function_3d - stored_v_coul;// Stored coulomb potential from ground state + mutable real_function_3d stored_v_nuc; // Stored nuclear potential from ground state + mutable real_function_3d stored_v_coul;// Stored coulomb potential from ground state // Ground state orbitals and energies vector_real_function_3d ground_orbitals{}; @@ -533,8 +467,7 @@ class ResponseBase { functionT ground_density;// ground state density - mutable response_space - stored_potential;// The ground state potential, stored only + mutable response_space stored_potential;// The ground state potential, stored only // if store_potential is true (default is double vtol{}; @@ -580,22 +513,19 @@ class ResponseBase { double safety = 0.1; vtol = FunctionDefaults<3>::get_thresh() * safety; - shared_coulomb_operator = - poperatorT(CoulombOperatorPtr(world, r_params.lo(), thresh)); + shared_coulomb_operator = poperatorT(CoulombOperatorPtr(world, r_params.lo(), thresh)); gradop = gradient_operator(world); potential_manager = std::make_shared(molecule, "a"); potential_manager->make_nuclear_potential(world); // GaussianConvolution1DCache::map.clear();//(TODO:molresponse-What // Create the masking function - mask = real_function_3d( - real_factory_3d(world).f(mask3).initial_level(4).norefine()); + 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", + print("\nSolving NDIM=", 3, " with thresh", thresh, " k", FunctionDefaults<3>::get_k(), " dconv", std::max(thresh, r_params.dconv()), "\n"); } } @@ -604,8 +534,7 @@ class ResponseBase { auto make_ground_density(World &world) const -> functionT; - auto ComputeHamiltonianPair(World &world) const - -> std::pair, Tensor>; + auto ComputeHamiltonianPair(World &world) const -> std::pair, Tensor>; auto Coulomb(World &world) const -> real_function_3d; @@ -619,82 +548,71 @@ class ResponseBase { void load_balance_chi(World &world); - auto make_bsh_operators_response(World &world, double &shift, - const double omega) const - -> vector; + auto make_bsh_operators_response(World &world, double &shift, const double omega) const -> vector; - auto kain_x_space_update(World &world, const X_space &chi, - const X_space &residual_chi, + auto kain_x_space_update(World &world, const X_space &chi, const X_space &residual_chi, response_solver &kain_x_space) -> X_space; - void x_space_step_restriction(World &world, const X_space &old_Chi, - X_space &temp, bool restrict_y, + void x_space_step_restriction(World &world, const X_space &old_Chi, X_space &temp, bool restrict_y, const double &max_bsh_rotation); +#if defined(__has_include) +#if __has_include() +#define MADCHEM_HAS_STD_FILESYSTEM +// is not reliably usable on Linux with gcc < 9 +#if defined(__GNUC__) +#if __GNUC__ >= 7 && __GNUC__ < 9 +#undef MADCHEM_HAS_STD_FILESYSTEM +#endif +#endif +#if defined(MADCHEM_HAS_STD_FILESYSTEM) + void PlotGroundDensityVTK(World &world, const ResponseBase &calc); + void plotResponseOrbitals(World &world, size_t iteration, const response_space &x_response, + const response_space &y_response, const ResponseParameters &responseParameters, + const GroundStateCalculation &g_params); + +#endif +#endif +#endif + + + static auto orbital_load_balance(World &world, const gamma_orbitals &, double load_balance) -> gamma_orbitals; + + auto compute_gamma_tda(World &world, const gamma_orbitals &density, const XCOperator &xc) const + -> X_space; -// void plotResponseOrbitals(World &world, size_t iteration, -// const response_space &x_response, -// const response_space &y_response, -// const ResponseParameters &responseParameters, -// const GroundStateCalculation &g_params); - - static auto orbital_load_balance(World &world, const gamma_orbitals &, - double load_balance) -> gamma_orbitals; - - auto compute_gamma_tda(World &world, const gamma_orbitals &density, - const XCOperator &xc) const -> X_space; - - auto compute_gamma_static(World &world, const gamma_orbitals &, - const XCOperator &xc) const -> X_space; + auto compute_gamma_static(World &world, const gamma_orbitals &, const XCOperator &xc) const -> X_space; - auto compute_gamma_full(World &world, const gamma_orbitals &, - const XCOperator &xc) const -> X_space; - auto compute_gamma(World &world, const gamma_orbitals &, - const XCOperator &xc) const -> X_space; + auto compute_gamma_full(World &world, const gamma_orbitals &, const XCOperator &xc) const -> X_space; + auto compute_gamma(World &world, const gamma_orbitals &, const XCOperator &xc) const -> X_space; - auto compute_V0X(World &world, const X_space &X, - const XCOperator &xc, bool compute_Y) const - -> X_space; + auto compute_V0X(World &world, const X_space &X, const XCOperator &xc, bool compute_Y) const -> X_space; - auto compute_lambda_X(World &world, const X_space &chi, - XCOperator &xc, + auto compute_lambda_X(World &world, const X_space &chi, XCOperator &xc, const std::string &calc_type) const -> X_space; - auto compute_theta_X(World &world, const X_space &chi, - const vector_real_function_3d &rho1, - const XCOperator &xc, - const std::string &calc_type) const -> X_space; + auto compute_theta_X(World &world, const X_space &chi, const vector_real_function_3d &rho1, + const XCOperator &xc, const std::string &calc_type) const -> X_space; - auto compute_F0X(World &world, const X_space &X, - const XCOperator &xc, bool compute_Y) const - -> X_space; + auto compute_F0X(World &world, const X_space &X, const XCOperator &xc, bool compute_Y) const -> X_space; - void analyze_vectors(World &world, const vecfuncT &x, - const std::string &response_state); + void analyze_vectors(World &world, const vecfuncT &x, const std::string &response_state); - auto project_ao_basis(World &world, const AtomicBasisSet &aobasis) - -> vecfuncT; + auto project_ao_basis(World &world, const AtomicBasisSet &aobasis) -> vecfuncT; - static auto project_ao_basis_only(World &world, - const AtomicBasisSet &aobasis, - const Molecule &mol) -> vecfuncT; + static auto project_ao_basis_only(World &world, const AtomicBasisSet &aobasis, const Molecule &mol) -> vecfuncT; void converged_to_json(json &j); - auto update_residual(World &world, const X_space &chi, const X_space &g_chi, - const std::string &calc_type, - const Tensor &old_residuals, - const X_space &xres_old) -> residuals; + auto update_residual(World &world, const X_space &chi, const X_space &g_chi, const std::string &calc_type, + const Tensor &old_residuals, const X_space &xres_old) -> residuals; - auto compute_response_potentials(World &world, const X_space &chi, - XCOperator &xc, - const std::string &calc_type) const - -> std::tuple; + auto compute_response_potentials(World &world, const X_space &chi, XCOperator &xc, + const std::string &calc_type) const -> std::tuple; // compute exchange |i> - auto exchangeHF(const vecfuncT &ket, const vecfuncT &bra, - const vecfuncT &vf) const -> vecfuncT { + auto exchangeHF(const vecfuncT &ket, const vecfuncT &bra, const vecfuncT &vf) const -> vecfuncT { World &world = ket[0].world(); auto n = bra.size(); auto nf = ket.size(); @@ -714,8 +632,7 @@ class ResponseBase { // i-j sym for (int i = 0; i < n; ++i) { // for each |i> - vecfuncT psi_f = - mul_sparse(world, bra[i], vf, mul_tol, true);/// was vtol + vecfuncT psi_f = mul_sparse(world, bra[i], vf, mul_tol, true);/// was vtol truncate(world, psi_f, tol, true); // apply to vector of products .. ... psi_f = apply(world, *shared_coulomb_operator, psi_f); @@ -730,17 +647,13 @@ class ResponseBase { return Kf; } - static void print_inner(World &world, const std::string &name, - const X_space &left, const X_space &right); + static void print_inner(World &world, const std::string &name, const X_space &left, const X_space &right); - void function_data_to_json(json &j_mol_in, size_t iter, - const Tensor &x_norms, - const Tensor &x_abs_norms, - const Tensor &rho_norms, + void function_data_to_json(json &j_mol_in, size_t iter, const Tensor &x_norms, + const Tensor &x_abs_norms, const Tensor &rho_norms, const Tensor &rho_res_norms); X_space compute_TX(World &world, const X_space &X, bool compute_Y) const; - vecfuncT update_density(World &world, const X_space &chi, - const vecfuncT &old_density) const; + vecfuncT update_density(World &world, const X_space &chi, const vecfuncT &old_density) const; }; @@ -755,8 +668,7 @@ class ResponseBase { void check_k(World &world, X_space &Chi, double thresh, int k); -auto add_randomness(World &world, const response_space &f, double magnitude) - -> response_space; +auto add_randomness(World &world, const response_space &f, double magnitude) -> response_space; void normalize(World &world, response_space &f); @@ -768,8 +680,7 @@ static auto kronecker(size_t l, size_t n) -> double { return 0.0; } -auto solid_harmonics(World &world, int n) - -> std::map, real_function_3d>; +auto solid_harmonics(World &world, int n) -> std::map, real_function_3d>; /*** * @brief Prints the norms of the functions of a response space @@ -787,8 +698,8 @@ auto make_xyz_functions(World &world) -> vector_real_function_3d; // Selects from a list of functions and energies the k functions with the // lowest energy -auto select_functions(World &world, response_space f, Tensor &energies, - size_t k, size_t print_level) -> response_space; +auto select_functions(World &world, response_space f, Tensor &energies, size_t k, size_t print_level) + -> response_space; // Sorts the given tensor of eigenvalues and // response functions @@ -808,26 +719,20 @@ auto gram_schmidt(World &world, const response_space &f) -> response_space; /// \param x /// \param y /// \return -auto transition_density(World &world, const vector_real_function_3d &orbitals, - const response_space &x, const response_space &y) - -> vector_real_function_3d; +auto transition_density(World &world, const vector_real_function_3d &orbitals, const response_space &x, + const response_space &y) -> vector_real_function_3d; -auto transition_densityTDA(World &world, - const vector_real_function_3d &orbitals, - const response_space &x) -> vector_real_function_3d; +auto transition_densityTDA(World &world, const vector_real_function_3d &orbitals, const response_space &x) + -> vector_real_function_3d; -auto transform(World &world, const response_space &f, const Tensor &U) - -> response_space; +auto transform(World &world, const response_space &f, const Tensor &U) -> response_space; -auto transform(World &world, const X_space &x, const Tensor &U) - -> X_space; +auto transform(World &world, const X_space &x, const Tensor &U) -> X_space; // result(i,j) = inner(a[i],b[j]).sum() -auto expectation(World &world, const response_space &A, const response_space &B) - -> Tensor; +auto expectation(World &world, const response_space &A, const response_space &B) -> Tensor; -void inner_to_json(World &world, const std::string &name, - const Tensor &m_val, +void inner_to_json(World &world, const std::string &name, const Tensor &m_val, std::map> &data); class ResponseTester { @@ -838,21 +743,18 @@ class ResponseTester { p->check_k(world, thresh, FunctionDefaults<3>::get_k()); } - static X_space compute_gamma_full(World &world, ResponseBase *p, - double thresh) { + static X_space compute_gamma_full(World &world, ResponseBase *p, double thresh) { XCOperator xc = p->make_xc_operator(world); return X_space{}; } X_space compute_lambda_X(World &world, ResponseBase *p, double thresh) { XCOperator xc = p->make_xc_operator(world); - X_space gamma = - p->compute_lambda_X(world, p->Chi, xc, p->r_params.calc_type()); + X_space gamma = p->compute_lambda_X(world, p->Chi, xc, p->r_params.calc_type()); return gamma; } - std::pair compute_VFOX(World &world, ResponseBase *p, - bool compute_y) { + std::pair compute_VFOX(World &world, ResponseBase *p, bool compute_y) { XCOperator xc = p->make_xc_operator(world); X_space V = p->compute_V0X(world, p->Chi, xc, compute_y); X_space F = p->compute_F0X(world, p->Chi, xc, compute_y);