Skip to content

Commit

Permalink
squash commit
Browse files Browse the repository at this point in the history
  • Loading branch information
Gabrielgerez committed Jan 18, 2023
1 parent 5bc83e4 commit 837b038
Show file tree
Hide file tree
Showing 27 changed files with 341,154 additions and 82 deletions.
17 changes: 10 additions & 7 deletions python/mrchem/CUBEparser.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,15 +41,15 @@ def parse_files(user_dict, direction=None, state=None):
found = False
for key, val in cube_guess_dict.items():
if (direction is not None):
data_type = "_".join(key.split("_")[2:]) + f"_{direction:d}"
data_type = "_".join(key.split("_")[2:]) + f"_rsp_{direction:d}"
path_list = _get_paths(Path(val), rsp=True, direction=direction)
elif (state is not None):
data_type = "_".join(key.split("_")[2:]) + f"_{state:d}"
path_list = _get_paths(Path(val), exc=True, state=state)
data_type = "_".join(key.split("_")[2:]) + f"_exc_{state:d}"
path_list = _get_paths(Path(val), exc=True, state=state)
else:
data_type = "_".join(key.split("_")[2:])
path_list = _get_paths(Path(val))

if path_list:
found = found or True
_write_cube_vectors(path_list, data_type, world_unit, pc, vector_dir)
Expand All @@ -61,7 +61,6 @@ def parse_files(user_dict, direction=None, state=None):

def _write_cube_vectors(path_list, data_type, world_unit, pc, vector_dir):
cube_list = []

if not vector_dir.is_dir():
vector_dir.mkdir()

Expand All @@ -80,13 +79,17 @@ def _write_cube_vectors(path_list, data_type, world_unit, pc, vector_dir):

def _get_paths(path, rsp=False, exc=False, direction=None, state=None):
directory = path.parent
prefix = path.name if (not rsp) else f"{path.name}_rsp_{direction:d}" if (not exc) else f"{path.name}_exc_{state:d}"
prefix = path.name
if (rsp):
prefix = f"{path.name}_rsp_{direction:d}"
elif (exc):
prefix = f"{path.name}_exc_{state:d}"


if directory.is_dir():
path_l = [file.resolve() for file in directory.glob(f"{prefix}*.cube")]
else:
path_l = []

return path_l


Expand Down
15 changes: 7 additions & 8 deletions python/mrchem/helpers.py
Original file line number Diff line number Diff line change
Expand Up @@ -344,12 +344,12 @@ def write_rsp_calc(omega, user_dict, origin):
"file_y_p": f"{file_dict['guess_y_p']}_rsp_{dir:d}",
"file_y_a": f"{file_dict['guess_y_a']}_rsp_{dir:d}",
"file_y_b": f"{file_dict['guess_y_b']}_rsp_{dir:d}",
"file_CUBE_x_p": f"{vector_dir}CUBE_x_p_{dir:d}_vector.json",
"file_CUBE_x_a": f"{vector_dir}CUBE_x_a_{dir:d}_vector.json",
"file_CUBE_x_b": f"{vector_dir}CUBE_x_b_{dir:d}_vector.json",
"file_CUBE_y_p": f"{vector_dir}CUBE_y_p_{dir:d}_vector.json",
"file_CUBE_y_a": f"{vector_dir}CUBE_y_a_{dir:d}_vector.json",
"file_CUBE_y_b": f"{vector_dir}CUBE_y_b_{dir:d}_vector.json",
"file_CUBE_x_p": f"{vector_dir}CUBE_x_p_rsp_{dir:d}_vector.json",
"file_CUBE_x_a": f"{vector_dir}CUBE_x_a_rsp_{dir:d}_vector.json",
"file_CUBE_x_b": f"{vector_dir}CUBE_x_b_rsp_{dir:d}_vector.json",
"file_CUBE_y_p": f"{vector_dir}CUBE_y_p_rsp_{dir:d}_vector.json",
"file_CUBE_y_a": f"{vector_dir}CUBE_y_a_rsp_{dir:d}_vector.json",
"file_CUBE_y_b": f"{vector_dir}CUBE_y_b_rsp_{dir:d}_vector.json",
}
if rsp_dict["write_orbitals"]:
path_orbitals = rsp_dict["path_orbitals"]
Expand Down Expand Up @@ -395,7 +395,6 @@ def write_exc_calc(user_dict, origin):

exc_calc["states"] = []
for r in range(exc_dict["n_states"]):

exc_state = {}

program_guess_type = user_guess_type
Expand Down Expand Up @@ -478,7 +477,7 @@ def write_rsp_fock(user_dict, wf_dict):
# Exchange-Correlation
if wf_dict["method_type"] in ["dft"]:
func_dict = []
for line in wf_dict["dft_funcs"].split("\n" ):
for line in wf_dict["dft_funcs"].split("\n"):
sp = line.split()
if len(sp) > 0:
func = sp[0].lower()
Expand Down
13 changes: 2 additions & 11 deletions src/chemistry/Molecule.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -74,27 +74,16 @@ Molecule::Molecule(const std::vector<std::string> &coord_str, int c, int m)
}

void Molecule::initPerturbedOrbitals(bool dynamic, int n_states) {
MSG_INFO("before reserving the size of the vectors");
std::cout << "size of state_vector to initialize " << n_states << "\n";
// this->orbitals_x.reserve(n_states);
// this->orbitals_y.reserve(n_states);
// MSG_INFO("after reserving the size of the vectors");
std::cout << "size of state_vector after reserving " << this->orbitals_x.size() << "\n";
MSG_INFO("before loop to add vector pointers to nstate vecotr ");
for (auto i = 0; i < n_states; i++) {
auto orbital_x = std::make_shared<OrbitalVector>();
std::cout << "in the loop \n";
this->orbitals_x.push_back(orbital_x);
// this->orbitals_x[i] = std::make_shared<OrbitalVector>();
std::cout << "size of state_vector to initialize " << this->orbitals_x.size() << "\n";
if (dynamic) {
auto orbital_y = std::make_shared<OrbitalVector>();
this->orbitals_y.push_back(orbital_y);
} else {
this->orbitals_y.push_back(orbital_x);
}
}
std::cout << "size of state_vector after loop " << this->orbitals_x.size() << "\n";
}

/** @brief Return number of electrons */
Expand Down Expand Up @@ -225,6 +214,7 @@ void Molecule::printGeometry() const {
void Molecule::printEnergies(const std::string &txt) const {
energy.print(txt);
epsilon.print(txt);
if (txt == "final") omega.print(txt);
}

/** @brief Pretty output of molecular properties
Expand Down Expand Up @@ -258,6 +248,7 @@ nlohmann::json Molecule::json() const {

json_out["scf_energy"] = energy.json();
json_out["orbital_energies"] = epsilon.json();
json_out["excitation_energies"] = omega.json();
if (not dipole.empty()) json_out["dipole_moment"] = {};
if (not quadrupole.empty()) json_out["quadrupole_moment"] = {};
if (not polarizability.empty()) json_out["polarizability"] = {};
Expand Down
5 changes: 5 additions & 0 deletions src/chemistry/Molecule.h
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@
#include "properties/Magnetizability.h"
#include "properties/NMRShielding.h"
#include "properties/OrbitalEnergies.h"
#include "properties/ExcitationEnergies.h"
#include "properties/Polarizability.h"
#include "properties/QuadrupoleMoment.h"
#include "properties/SCFEnergy.h"
Expand Down Expand Up @@ -101,6 +102,8 @@ class Molecule final {
auto getOrbitalsX_p(int state = 0) const { return this->orbitals_x[state]; }
auto getOrbitalsY_p(int state = 0) const { return this->orbitals_y[state]; }
auto getCavity_p() const { return this->cavity; }
auto getStatesX() const { return this->orbitals_x; }
auto getStatesY() const { return this->orbitals_y; }

nlohmann::json json() const;
void printGeometry() const;
Expand All @@ -112,6 +115,7 @@ class Molecule final {

SCFEnergy &getSCFEnergy() { return this->energy; }
OrbitalEnergies &getOrbitalEnergies() { return this->epsilon; }
ExcitationEnergies &getExcitationEnergies() { return this->omega; }
DipoleMoment &getDipoleMoment(const std::string &id) { return this->dipole.at(id); }
QuadrupoleMoment &getQuadrupoleMoment(const std::string &id) { return this->quadrupole.at(id); }
Polarizability &getPolarizability(const std::string &id) { return this->polarizability.at(id); }
Expand Down Expand Up @@ -140,6 +144,7 @@ class Molecule final {
// Properties
SCFEnergy energy{};
OrbitalEnergies epsilon{};
ExcitationEnergies omega{};
PropertyMap<DipoleMoment> dipole{};
PropertyMap<QuadrupoleMoment> quadrupole{};
PropertyMap<Polarizability> polarizability{};
Expand Down
42 changes: 28 additions & 14 deletions src/driver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -700,7 +700,7 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) {
*/
json driver::rsp::run(const json &json_rsp, Molecule &mol) {
auto run_exc = json_rsp.contains("states");
if (not run_exc) {
if (!run_exc) {
print_utils::headline(0, "Computing Linear Response Wavefunction");
} else {
print_utils::headline(0, "Optimizing Transition Orbitals and Frequencies");
Expand Down Expand Up @@ -743,13 +743,13 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) {
auto run_tda = false;
auto iter_var = 0;
RankOneOperator<3> h_1;
if (not run_exc) {
if ( !run_exc) {
omega = json_rsp["frequency"];
dynamic = json_rsp["dynamic"];
iter_var = json_rsp["components"].size();
mol.initPerturbedOrbitals(dynamic, 1);
const auto &json_pert = json_rsp["perturbation"];
auto h_1 = driver::get_operator<3>(json_pert["operator"], json_pert);
h_1 = driver::get_operator<3>(json_pert["operator"], json_pert);
json_out["perturbation"] = json_pert["operator"];
json_out["frequency"] = omega;
json_out["components"] = {};
Expand All @@ -761,19 +761,19 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) {
} else {
MSG_ERROR("Response calculation badly initialized. No components or excited states")
}
FockBuilder F_1;
const auto &json_fock_1 = json_rsp["fock_operator"];
driver::build_fock_operator(json_fock_1, mol, F_1, 1);

const auto &json_fock = json_rsp["fock_operator"];
for (auto i = 0; i < iter_var; i++) {
FockBuilder F_1;
if (run_exc) {
driver::build_fock_operator(json_fock, mol, F_1, 1, i, run_tda);
} else {
driver::build_fock_operator(json_fock, mol, F_1, 1);
F_1.clear();
driver::build_fock_operator(json_fock_1, mol, F_1, 1, i, run_tda);
}

json comp_out = {};
const auto &json_comp = not run_exc ? json_rsp["components"][i] : json_rsp["states"][i];
if (not run_exc) { F_1.perturbation() = h_1[i]; }
const auto &json_comp = !run_exc ? json_rsp["components"][i] : json_rsp["states"][i];
if (!run_exc) { F_1.perturbation() = h_1[i]; }

///////////////////////////////////////////////////////////
/////////////// Setting Up Initial Guess //////////////
Expand Down Expand Up @@ -814,9 +814,24 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) {
solver.setOrthPrec(orth_prec);
comp_out["exc_solver"] = solver.optimize(mol, F_0, F_1, i);
json_out["success"] = comp_out["exc_solver"]["converged"];
//TODO implement write orbitals for the excited states and property calculations.
if (json_out["success"]) {
if (json_comp.contains("write_orbitals")) rsp::write_orbitals(json_comp["write_orbitals"], mol, dynamic);
if (json_rsp.contains("properties")) rsp::calc_properties(json_rsp["properties"], mol, 2, json_rsp["frequency"], i);
Timer t_lap;
t_lap.start();
mrcpp::print::header(2, "Excitation energy");
double w_au = comp_out["exc_solver"]["frequency"];
auto w_cm = PhysicalConstants::get("hartree2wavenumbers") * w_au;
auto l_nm = (1.0e7 / w_cm);

mrcpp::print::header(0, "Exciation energy (state " + std::to_string(i+1) + ")");
print_utils::scalar(0, "Wavelength", l_nm, "(nm)");
print_utils::scalar(0, "Frequency", w_au, "(au)");
print_utils::scalar(0, " ", w_cm, "(cm-1)");
mrcpp::print::separator(0, '-');
mrcpp::print::separator(0, '=', 2);
mrcpp::print::footer(2, t_lap, 2);

}
mol.getOrbitalsX(i).clear(); // Clear orbital vector
mol.getOrbitalsY(i).clear(); // Clear orbital vector
Expand Down Expand Up @@ -853,17 +868,15 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) {
///////////////////////////////////////////////////////////
//////////// Compute Response Properties //////////////
///////////////////////////////////////////////////////////

if (json_out["success"]) {
if (json_comp.contains("write_orbitals")) rsp::write_orbitals(json_comp["write_orbitals"], mol, dynamic);
if (json_rsp.contains("properties")) rsp::calc_properties(json_rsp["properties"], mol, 2, json_rsp["frequency"], i);
if (json_rsp.contains("properties")) rsp::calc_properties(json_rsp["properties"], mol, i, json_rsp["frequency"]);
}
mol.getOrbitalsX(0).clear(); // Clear orbital vector
mol.getOrbitalsY(0).clear(); // Clear orbital vector
json_out["components"].push_back(comp_out);
}
}

F_0.clear();
mpi::barrier(mpi::comm_orb);
if (run_exc) {
Expand Down Expand Up @@ -1041,6 +1054,7 @@ void driver::rsp::calc_properties(const json &json_prop, Molecule &mol, int dir,
mrcpp::print::footer(2, t_lap, 2);
if (plevel == 1) mrcpp::print::time(1, "NMR shielding (para)", t_lap);
}


if (json_prop.contains("hyperpolarizability")) MSG_ERROR("Hyperpolarizability not implemented");
if (json_prop.contains("geometric_derivative")) MSG_ERROR("Geometric derivative not implemented");
Expand Down
90 changes: 90 additions & 0 deletions src/properties/ExcitationEnergies.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
/*
* MRChem, a numerical real-space code for molecular electronic structure
* calculations within the self-consistent field (SCF) approximations of quantum
* chemistry (Hartree-Fock and Density Functional Theory).
* Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
*
* This file is part of MRChem.
*
* MRChem is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* MRChem is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with MRChem. If not, see <https://www.gnu.org/licenses/>.
*
* For information on the complete list of contributors to MRChem, see:
* <https://mrchem.readthedocs.io/>
*/

#pragma once

#include <nlohmann/json.hpp>

#include "mrchem.h"
#include "utils/print_utils.h"

/** @class ExcitationEnergies
*
* @brief Simple POD container to hold the Excitation energies
*
*/

namespace mrchem {

// clang-format off
class ExcitationEnergies final {
public:
std::vector<double> &getOmega() { return this->omega; }

const std::vector<double> &getOmega() const { return this->omega; }

void print(const std::string &id) const {
auto pprec = 2 * mrcpp::Printer::getPrecision();
auto w0 = mrcpp::Printer::getWidth() - 1;
auto w1 = 5;
auto w2 = 2 * w0 / 9;
auto w3 = w0 - 3 * w1 - 3 * w2;

std::stringstream o_head;
o_head << std::setw(w1) << "r";
o_head << std::string(w1 + w1 + w3 - 1, ' ') << ':';
o_head << std::setw(3 * w2) << "Omega";

// could use this to print both the starting guess and final
// converged frequencies, could help to check the quality of the guess
mrcpp::print::header(0, "Excitation Energies (" + id + ")");
println(0, o_head.str());
mrcpp::print::separator(0, '-');

for (int i = 0; i < this->omega.size(); i++) {
std::stringstream o_txt;
o_txt << std::setw(w1 - 1) << i+1;
print_utils::scalar(0, o_txt.str(), this->omega[i], "(au)", pprec);
}
mrcpp::print::separator(0, '=', 2);
}

nlohmann::json json() const {
const std::vector<double> &ome = getOmega();
nlohmann::json json_out;
for (auto i = 0; i < ome.size(); i++) {
std::stringstream state_key;
state_key << "state-" << i+1;
json_out[state_key.str()] = ome[i];
}
return json_out;
}

private:
std::vector<double> omega;
};
// clang-format on

} // namespace mrchem
1 change: 1 addition & 0 deletions src/properties/properties_fwd.h
Original file line number Diff line number Diff line change
Expand Up @@ -33,5 +33,6 @@ class Magnetizability;
class NMRShielding;
class Polarizability;
class OrbitalEnergies;
class ExcitationEnergies;

} // namespace mrchem
25 changes: 25 additions & 0 deletions src/qmfunctions/OrbitalIterator.cpp
Original file line number Diff line number Diff line change
@@ -1,3 +1,28 @@
/*
* MRChem, a numerical real-space code for molecular electronic structure
* calculations within the self-consistent field (SCF) approximations of quantum
* chemistry (Hartree-Fock and Density Functional Theory).
* Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
*
* This file is part of MRChem.
*
* MRChem is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* MRChem is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with MRChem. If not, see <https://www.gnu.org/licenses/>.
*
* For information on the complete list of contributors to MRChem, see:
* <https://mrchem.readthedocs.io/>
*/

/*
* Mrchem, a numerical real-space code for molecular electronic structure
* calculations within the self-consistent field (SCF) approximations of quantum
Expand Down
Loading

0 comments on commit 837b038

Please sign in to comment.