diff --git a/doc/schema_input.json b/doc/schema_input.json new file mode 100644 index 000000000..41554e5c6 --- /dev/null +++ b/doc/schema_input.json @@ -0,0 +1,261 @@ +"input": { + "molecule": { # Section for Molecule specification + "charge": int, # Total molecular charge + "multiplicity": int # Total spin multiplicity + "gauge_origin": array[float], # Gauge origin coordinate + "coords": array[ # Array of atoms + { # (one entry per atom) + "atom": string, # Atomic symbol + "xyz": array[float] # Nuclear Cartesian coordinate + } + ], + }, + "mpi": { # Section for MPI specification + "bank_size": int, # Number of MPI ranks in memory bank + "numerically_exact": bool, # Guarantee MPI invariant results + "shared_memory_size": int # Size (MB) of MPI shared memory blocks + }, + "mra": { # Section for MultiResolution Analysis + "basis_type": string, # Basis type (interpolating/legendre) + "basis_order": int, # Polynomial order of basis + "max_scale": int, # Maximum level of refinement + "min_scale": int, # Minimum level of refinement (root scale) + "boxes": array[int], # Number of root boxes + "corner": array[int], # Translation of first root box + "scaling_factor": array[float] # Coordinate scaling factor + }, + "printer": { # Section for printed output + "file_name": string, # Name of output file + "print_input": bool, # Dump user input in output file + "print_level": int, # Amount of printed output + "print_mpi": bool, # Use separate output file for each MPI + "print_prec": int, # Number of digits for printed output + "print_width": int # Line width of printed output + }, + "scf_calculation": { # Section for SCF specification + "fock_operator": { # Contributions to Fock operator + "kinetic_operator": { # Add Kinetic operator to Fock + "derivative": string # Type of derivative operator + }, + "nuclear_operator": { # Add Nuclear operator to Fock + "proj_prec": float, # Projection prec for potential + "smooth_prec": float, # Smoothing parameter for potential + "shared_memory": bool # Use shared memory for potential + }, + "coulomb_operator": { # Add Coulomb operator to Fock + "poisson_prec": float, # Build prec for Poisson operator + "shared_memory": bool # Use shared memory for potential + }, + "exchange_operator": { # Add Exchange operator to Fock + "poisson_prec": float, # Build prec for Poisson operator + "screen": bool # Use screening in Exchange operator + }, + "xc_operator": { # Add XC operator to Fock + "grid_prec": float, # Projection prec for density + "shared_memory": bool, # Use shared memory for potential + "xc_functional": { # XC functional specification + "spin": bool, # Use spin separated functional + "gamma": bool, # Use gamma density variable + "cutoff": float, # Cutoff value for small densities + "log_grad": bool, # Use log parametrization of gradient + "derivative": string, # Type of derivative operator + "functionals": array[ # Array of density functionals + { + "coef": float, # Numerical coefficient + "name": string # Functional name + } + ] + } + }, + "external_operator": { # Add external field operator to Fock + "electric_field": array[float], # Electric field vector + "r_O": array[float] # Gauge orgigin for electric field + } + }, + "initial_guess": { # Initial guess specification + "type": string, # Type of initial guess + "prec": float, # Precision for initial guess + "zeta": int, # Zeta quality for AO basis + "method": string, # Name of method for initial energy + "localize": bool, # Use localized orbitals + "restricted": bool, # Use spin restricted orbitals + "file_chk": string, # Path to checkpoint file + "file_basis": string, # Path to GTO basis file + "file_gto_a": string, # Path to GTO MO file (alpha) + "file_gto_b": string, # Path to GTO MO file (beta) + "file_gto_p": string, # Path to GTO MO file (paired) + "file_phi_a": string, # Path to MW orbital file (alpha) + "file_phi_b": string, # Path to MW orbital file (beta) + "file_phi_p": string # Path to MW orbital file (paired) + }, + "scf_solver": { # SCF solver specification + "kain": int, # Length of KAIN history + "max_iter": int, # Maximum number of iterations + "method": string, # Name of electronic structure method + "rotation": int, # Iterations between localize/diagonalize + "localize": bool, # Use localized orbitals + "checkpoint": bool, # Save checkpoint file + "file_chk": string, # Name of checkpoint file + "start_prec": float, # Start precision for solver + "final_prec": float, # Final precision for solver + "helmholtz_prec": float, # Precision for Helmholtz operators + "orbital_thrs": float, # Convergence threshold orbitals + "energy_thrs":float, # Convergence threshold energy + }, + "properties": { # Collection of ground state properties to compute + "dipole_moment": { # Collection of dipole moments + id (string): { # Unique identifier: 'dip-${number}' + "precision": float, # Operator precision + "operator": string, # Operator used for property + "r_O": array[float] # Operator gauge origin + } + }, + "quadrupole_moment": { # Collection of quadrupole moments + id (string): { # Unique identifier: 'quad-${number}' + "precision": float, # Operator precision + "operator": string, # Operator used for property + "r_O": array[float] # Operator gauge origin + } + } + } + }, + "rsp_calculations": { # Collection of response calculations + id (string): { # Response identifier, e.g. 'ext_el-${frequency}' + "dynamic": bool, # Use dynamic response solver + "frequency": float, # Perturbation frequency + "perturbation": { # Perturbation operator + "operator": string # Operator used in response calculation + }, + "components": array[ # Array of perturbation components + { # (one per Cartesian direction) + "initial_guess": { # Initial guess specification + "type": string, # Type of initial guess + "prec": float, # Precision for initial guess + "file_chk_x": string, # Path to checkpoint file for X + "file_chk_y": string, # Path to checkpoint file for Y + "file_x_a": string, # Path to MW file for X (alpha) + "file_x_b": string, # Path to MW file for X (beta) + "file_x_p": string, # Path to MW file for X (paired) + "file_y_a": string, # Path to MW file for Y (alpha) + "file_y_b": string, # Path to MW file for Y (beta) + "file_y_p": string # Path to MW file for Y (paired) + }, + "rsp_solver": { # Response solver specification + "kain": int, # Length of KAIN history + "max_iter": int, # Maximum number of iterations + "method": string, # Name of electronic structure method + "checkpoint": bool, # Save checkpoint file + "file_chk_x": string, # Name of X checkpoint file + "file_chk_y": string, # Name of Y checkpoint file + "orth_prec": float, # Precision for orthogonalization + "start_prec": float, # Start precision for solver + "final_prec": float, # Final precision for solver + "helmholtz_prec": float, # Precision for Helmholtz operators + "orbital_thrs": float, # Convergence threshold orbitals + "property_thrs": float # Convergence threshold property + } + } + ], + "properties": { # Collection of response properties to compute + "polarizability": { # Collection of polarizabilities + id (string): { # Unique identifier, default 'pol-${frequency}' + "precision": float, # Operator precision + "operator": string, # Operator used for property + "r_O": array[float] # Operator gauge origin + } + }, + "magnetizability": { # Collection of magnetizabilities + id (string): { # Unique identifier: 'mag-${frequency}' + "frequency": float, # Perturbation frequency + "precision": float, # Operator precision + "dia_operator": string, # Operator used for diamagnetic property + "para_operator": string, # Operator used for paramagnetic property + "derivative": string, # Operator derivative type + "r_O": array[float] # Operator gauge origin + } + }, + "nmr_shielding": { # Collection of NMR shieldings + id (string): { # Unique identifier: 'nmr-${nuc_idx}${atomic_symbol}' + "precision": float, # Operator precision + "dia_operator": string, # Operator used for diamagnetic property + "para_operator": string, # Operator used for paramagnetic property + "derivative": string, # Operator derivative type + "smoothing": float, # Operator smoothing parameter + "r_O": array[float], # Operator gauge origin + "r_K": array[float] # Nuclear coordinate + } + } + }, + "fock_operator": { # Contributions to perturbed Fock operator + "coulomb_operator": { # Add Coulomb operator to Fock + "poisson_prec": float, # Build prec for Poisson operator + "shared_memory": bool # Use shared memory for potential + }, + "exchange_operator": { # Add Exchange operator to Fock + "poisson_prec": float, # Build prec for Poisson operator + "screen": bool # Use screening in Exchange operator + }, + "xc_operator": { # Add XC operator to Fock + "grid_prec": float, # Projection prec for density + "shared_memory": bool, # Use shared memory for potential + "xc_functional": { # XC functional specification + "spin": bool, # Use spin separated functional + "gamma": bool, # Use gamma density variable + "cutoff": float, # Cutoff value for small densities + "log_grad": bool, # Use log parametrization of gradient + "derivative": string, # Type of derivative operator + "functionals": array[ # Array of density functionals + { + "coef": float, # Numerical coefficient + "name": string # Functional name + } + ] + } + } + }, + "unperturbed": { # Section for unperturbed part of response + "prec": float, # Precision used for unperturbed system + "localize": bool, # Use localized unperturbed orbitals + "fock_operator": { # Contributions to unperturbed Fock operator + "kinetic_operator": { # Add Kinetic operator to Fock + "derivative": string # Type of derivative operator + }, + "nuclear_operator": { # Add Nuclear operator to Fock + "proj_prec": float, # Projection prec for potential + "smooth_prec": float, # Smoothing parameter for potential + "shared_memory": bool # Use shared memory for potential + }, + "coulomb_operator": { # Add Coulomb operator to Fock + "poisson_prec": float, # Build prec for Poisson operator + "shared_memory": bool # Use shared memory for potential + }, + "exchange_operator": { # Add Exchange operator to Fock + "poisson_prec": float, # Build prec for Poisson operator + "screen": bool # Use screening in Exchange operator + }, + "xc_operator": { # Add XC operator to Fock + "grid_prec": float, # Projection prec for density + "shared_memory": bool, # Use shared memory for potential + "xc_functional": { # XC functional specification + "spin": bool, # Use spin separated functional + "gamma": bool, # Use gamma density variable + "cutoff": float, # Cutoff value for small densities + "log_grad": bool, # Use log parametrization of gradient + "derivative": string, # Type of derivative operator + "functionals": array[ # Array of density functionals + { + "coef": float, # Numerical coefficient + "name": string # Functional name + } + ] + } + }, + "external_operator": { # Add external field operator to Fock + "electric_field": array[float], # Electric field vector + "r_O": array[float] # Gauge orgigin for electric field + } + } + } + } + } +} diff --git a/doc/schema_output.json b/doc/schema_output.json new file mode 100644 index 000000000..da4536e2b --- /dev/null +++ b/doc/schema_output.json @@ -0,0 +1,144 @@ +"output": { + "properties": { # Collection of final properties + "charge": int, # Total molecular charge + "multiplicity": int, # Total spin multiplicity + "center_of_mass": array[float], # Center of mass coordinate + "geometry": array[ # Array of atoms + { # (one entry per atom) + "symbol": string, # Atomic symbol + "xyz": array[float] # Cartesian coordinate + } + ], + "orbital_energies": { # Collection of orbital energies + "spin": array[string], # Array of spins ('p', 'a' or 'b') + "energy": array[float], # Array of energies + "occupation": array[int], # Array of orbital occupations + "sum_occupied": float # \sum_i occupation[i]*energy[i] + }, + "scf_energy": { # Collection of energy contributions + "E_kin": float, # Kinetic energy + "E_nn": float, # Classical nuclear-nuclear interaction + "E_en": float, # Classical electron-nuclear interaction + "E_ee": float, # Classical electron-electron interaction + "E_next": float, # Classical nuclear-external field interaction + "E_eext": float, # Classical electron-external field interaction + "E_x": float, # Hartree-Fock exact exchange energy + "E_xc": float, # DFT exchange-correlation energy + "E_el": float, # Sum of electronic contributions + "E_nuc": float, # Sum of nuclear contributions + "E_tot": float # Sum of all contributions + }, + "dipole_moment": { # Collection of electric dipole moments + id (string): { # Unique identifier, default 'dip-1' + "r_O": array[float], # Gauge origin vector + "vector": array[float], # Total dipole vector + "vector_el": array[float], # Electronic dipole vector + "vector_nuc": array[float], # Nuclear dipole vector + "magnitude": float # Magnitude of total vector + } + }, + "quadrupole_moment": { # Collection of electric quadrupole moments + id (string): { # Unique identifier, default 'quad-1' + "r_O": array[float], # Gauge origin vector + "tensor": array[float], # Total quadrupole tensor + "tensor_el": array[float], # Electronic quadrupole tensor + "tensor_nuc": array[float] # Nuclear quadrupole tensor + } + }, + "polarizability": { # Collection of polarizabilities + id (string): { # Unique identifier, default 'pol-${frequency}' + "frequency": float, # Perturbation frequency + "r_O": array[float], # Gauge origin vector + "tensor": array[float], # Full polarizability tensor + "isotropic_average": float # Diagonal average + } + }, + "magnetizability": { # Collection of magnetizability + id (string): { # Unique identifier, default 'mag-${frequency}' + "frequency": float, # Perturbation frequency + "r_O": array[float], # Gauge origin vector + "tensor": array[float], # Full magnetizability tensor + "tensor_dia": array[float], # Diamagnetic tensor + "tensor_para": array[float], # Paramagnetic tensor + "isotropic_average": float # Diagonal average + } + }, + "nmr_shielding": { # Collection of NMR shielding tensors + id (string): { # Unique identifier, default 'nmr-${index}+${symbol}' + "r_O": array[float], # Gauge origin vector + "r_K": array[float], # Nuclear coordinate vector + "tensor": array[float], # Full NMR shielding tensor + "tensor_dia": array[float], # Diamagnetic tensor + "tensor_para": array[float], # Paramagnetic tensor + "diagonalized_tensor": array[float], # Diagonalized tensor used for (an)isotropy + "isotropic_average": float, # Diagonal average + "anisotropy": float # Anisotropy of tensor + } + } + }, + "scf_calculation": { # Ground state SCF calculation + "success": bool, # SCF finished successfully + "initial_energy": { # Energy computed from initial orbitals + "E_kin": float, # Kinetic energy + "E_nn": float, # Classical nuclear-nuclear interaction + "E_en": float, # Classical electron-nuclear interaction + "E_ee": float, # Classical electron-electron interaction + "E_next": float, # Classical nuclear-external field interaction + "E_eext": float, # Classical electron-external field interaction + "E_x": float, # Hartree-Fock exact exchange energy + "E_xc": float, # DFT exchange-correlation energy + "E_el": float, # Sum of electronic contributions + "E_nuc": float, # Sum of nuclear contributions + "E_tot": float # Sum of all contributions + }, + "scf_solver": { # Details from SCF optimization + "converged": bool, # Optimization converged + "wall_time": float, # Wall time (sec) for SCF optimization + "cycles": array[ # Array of SCF cycles + { # (one entry per cycle) + "energy_total": float, # Current total energy + "energy_update": float, # Current energy update + "mo_residual": float, # Current orbital residual + "wall_time": float, # Wall time (sec) for SCF cycle + "energy_terms": { # Energy contributions + "E_kin": float, # Kinetic energy + "E_nn": float, # Classical nuclear-nuclear interaction + "E_en": float, # Classical electron-nuclear interaction + "E_ee": float, # Classical electron-electron interaction + "E_next": float, # Classical nuclear-external field interaction + "E_eext": float, # Classical electron-external field interaction + "E_x": float, # Hartree-Fock exact exchange energy + "E_xc": float, # DFT exchange-correlation energy + "E_el": float, # Sum of electronic contributions + "E_nuc": float, # Sum of nuclear contributions + "E_tot": float # Sum of all contributions + } + } + ] + } + }, + "rsp_calculations": { # Collection of response calculations + id (string): { # Response identifier, e.g. 'ext_el-${frequency}' + "success": bool, # Response finished successfully + "frequency": float, # Frequency of perturbation + "perturbation": string, # Name of perturbation operator + "components": array[ # Array of operator components + { # (one entry per Cartesian direction) + "rsp_solver": { # Details from response optimization + "wall_time": float, # Wall time (sec) for response calculation + "converged": bool, # Optimization converged + "cycles": array[ # Array of response cycles + { # (one entry per cycle) + "symmetric_property": float, # Property computed from perturbation operator + "property_update": float, # Current symmetric property update + "mo_residual": float, # Current orbital residual + "wall_time": float # Wall time (sec) for response cycle + } + ] + } + } + ] + } + } +} + diff --git a/doc/users/program_json.rst b/doc/users/program_json.rst new file mode 100644 index 000000000..58127d7cf --- /dev/null +++ b/doc/users/program_json.rst @@ -0,0 +1,16 @@ +----------------------------- +The program input/output file +----------------------------- + + +Input schema +------------ + +.. literalinclude:: schema_input.json + :language: JSON + +Output schema +------------- + +.. literalinclude:: schema_output.json + :language: JSON diff --git a/doc/users/user_out.rst b/doc/users/user_out.rst new file mode 100644 index 000000000..7027e754c --- /dev/null +++ b/doc/users/user_out.rst @@ -0,0 +1,4 @@ +---------------- +User output file +---------------- + diff --git a/external/upstream/fetch_nlohmann_json.cmake b/external/upstream/fetch_nlohmann_json.cmake index d7ab8e54c..8c9ef9e78 100644 --- a/external/upstream/fetch_nlohmann_json.cmake +++ b/external/upstream/fetch_nlohmann_json.cmake @@ -1,4 +1,4 @@ -find_package(nlohmann_json 3.5 CONFIG QUIET) +find_package(nlohmann_json 3.6 CONFIG QUIET) if(TARGET nlohmann_json::nlohmann_json) get_target_property( @@ -13,7 +13,7 @@ else() FetchContent_Declare(nlohmann_json_sources QUIET URL - https://github.com/nlohmann/json/archive/v3.5.0.tar.gz + https://github.com/nlohmann/json/archive/v3.6.1.tar.gz ) FetchContent_GetProperties(nlohmann_json_sources) diff --git a/pilot/mrchem.cpp b/pilot/mrchem.cpp index f96431d68..cc5fd90a1 100644 --- a/pilot/mrchem.cpp +++ b/pilot/mrchem.cpp @@ -16,7 +16,7 @@ using namespace mrchem; int main(int argc, char **argv) { mpi::initialize(); - const auto json_input = mrenv::fetch_input(argc, argv); + const auto json_input = mrenv::fetch_json(argc, argv); mrenv::initialize(json_input); Timer timer; diff --git a/pilot/mrchem.in b/pilot/mrchem.in index 22f5692ff..2f2006701 100644 --- a/pilot/mrchem.in +++ b/pilot/mrchem.in @@ -61,7 +61,8 @@ def main(): # now that all keywords have sensible values, # we can translate user input into program input - program_dict = translate_input(user_dict) + program_dict = {} + program_dict["input"] = translate_input(user_dict) inp_name, ext_ext = os.path.splitext(inp_file_cmd) xfile = inp_name + '.json' @@ -70,7 +71,7 @@ def main(): if not dryrun: cmd = executable + ' ' + xfile - if program_dict["printer"]["print_input"]: + if program_dict["input"]["printer"]["print_input"]: subprocess.call('cat ' + inp_file_cmd, shell=True) subprocess.call(cmd, shell=True) diff --git a/src/chemistry/Molecule.cpp b/src/chemistry/Molecule.cpp index 4c3070c39..2920bd65c 100644 --- a/src/chemistry/Molecule.cpp +++ b/src/chemistry/Molecule.cpp @@ -72,10 +72,6 @@ Molecule::Molecule(const std::vector &coord_str, int c, int m) readCoordinateString(coord_str); } -void Molecule::initNuclearProperties(int nNucs) { - for (auto k = 0; k < nNucs; k++) nmr.push_back(nullptr); -} - void Molecule::initPerturbedOrbitals(bool dynamic) { if (dynamic) { this->orbitals_x = std::make_shared(); @@ -86,54 +82,6 @@ void Molecule::initPerturbedOrbitals(bool dynamic) { } } -/** @brief Return property SCFEnergy */ -SCFEnergy &Molecule::getSCFEnergy() { - if (energy == nullptr) energy = std::make_unique(); - return *energy; -} - -/** @brief Return property DipoleMoment */ -DipoleMoment &Molecule::getDipoleMoment() { - if (dipole == nullptr) dipole = std::make_unique(); - return *dipole; -} - -/** @brief Return property QuadrupoleMoment */ -QuadrupoleMoment &Molecule::getQuadrupoleMoment() { - if (quadrupole == nullptr) quadrupole = std::make_unique(); - return *quadrupole; -} - -/** @brief Return property Magnetizability */ -Magnetizability &Molecule::getMagnetizability() { - if (magnetizability == nullptr) magnetizability = std::make_unique(); - return *magnetizability; -} - -/** @brief Return property NMRShielding */ -NMRShielding &Molecule::getNMRShielding(int k) { - if (nmr.size() == 0) initNuclearProperties(getNNuclei()); - if (nmr[k] == nullptr) nmr[k] = std::make_unique(k, nuclei[k]); - return *nmr[k]; -} - -/** @brief Return property Polarizability */ -Polarizability &Molecule::getPolarizability(double omega) { - auto idx = -1; - for (auto i = 0; i < polarizability.size(); i++) { - auto omega_i = polarizability[i]->getFrequency(); - if (std::abs(omega_i - omega) < mrcpp::MachineZero) { - idx = i; - break; - } - } - if (idx < 0) { - polarizability.push_back(std::make_unique(omega)); - idx = polarizability.size() - 1; - } - return *polarizability[idx]; -} - /** @brief Return number of electrons */ int Molecule::getNElectrons() const { auto totZ = 0; @@ -253,31 +201,60 @@ void Molecule::printGeometry() const { o_sym << std::setw(w2) << nuc.getElement().getSymbol(); print_utils::coord(0, o_sym.str(), nuc.getCoord()); } + mrcpp::print::separator(0, '-'); print_utils::coord(0, "Center of mass", calcCenterOfMass()); print_utils::coord(0, "Gauge origin", getGaugeOrigin()); mrcpp::print::separator(0, '=', 2); } +void Molecule::printEnergies(const std::string &txt) const { + energy.print(txt); + epsilon.print(txt); +} + /** @brief Pretty output of molecular properties * * Only properties that have been initialized will be printed. */ void Molecule::printProperties() const { - const auto &Phi = getOrbitals(); - const auto &F_mat = getFockMatrix(); - orbital::print_eigenvalues(Phi, F_mat); + // For each std::map entry (first = id string, second = property) + // this will print the property with its id appearing in the header + for (const auto &dip : dipole) dip.second.print(dip.first); + for (const auto &qua : quadrupole) qua.second.print(qua.first); + for (const auto &pol : polarizability) pol.second.print(pol.first); + for (const auto &mag : magnetizability) mag.second.print(mag.first); + for (const auto &nmr : nmr_shielding) nmr.second.print(nmr.first); +} - if (this->energy != nullptr) this->energy->print(); - if (this->dipole != nullptr) this->dipole->print(); - if (this->quadrupole != nullptr) this->quadrupole->print(); - for (auto &pol : this->polarizability) { - if (pol != nullptr) pol->print(); - } - if (this->magnetizability != nullptr) this->magnetizability->print(); - for (auto &nmr_k : this->nmr) { - if (nmr_k != nullptr) nmr_k->print(); +nlohmann::json Molecule::json() const { + nlohmann::json json_out; + json_out["geometry"] = {}; + for (auto i = 0; i < getNNuclei(); i++) { + const auto &nuc = getNuclei()[i]; + nlohmann::json json_atom; + json_atom["symbol"] = nuc.getElement().getSymbol(); + json_atom["xyz"] = nuc.getCoord(); + json_out["geometry"].push_back(json_atom); } + json_out["charge"] = getCharge(); + json_out["multiplicity"] = getMultiplicity(); + json_out["center_of_mass"] = calcCenterOfMass(); + + json_out["scf_energy"] = energy.json(); + json_out["orbital_energies"] = epsilon.json(); + if (not dipole.empty()) json_out["dipole_moment"] = {}; + if (not quadrupole.empty()) json_out["quadrupole_moment"] = {}; + if (not polarizability.empty()) json_out["polarizability"] = {}; + if (not magnetizability.empty()) json_out["magnetizability"] = {}; + if (not nmr_shielding.empty()) json_out["nmr_shielding"] = {}; + for (const auto &dip : dipole) json_out["dipole_moment"][dip.first] = dip.second.json(); + for (const auto &qua : quadrupole) json_out["quadrupole_moment"][qua.first] = qua.second.json(); + for (const auto &pol : polarizability) json_out["polarizability"][pol.first] = pol.second.json(); + for (const auto &mag : magnetizability) json_out["magnetizability"][mag.first] = mag.second.json(); + for (const auto &nmr : nmr_shielding) json_out["nmr_shielding"][nmr.first] = nmr.second.json(); + + return json_out; } } // namespace mrchem diff --git a/src/chemistry/Molecule.h b/src/chemistry/Molecule.h index d831fca9e..f06df7cb5 100644 --- a/src/chemistry/Molecule.h +++ b/src/chemistry/Molecule.h @@ -38,12 +38,14 @@ #include #include -#include "MRCPP/Printer" +#include +#include #include "Nucleus.h" #include "properties/DipoleMoment.h" #include "properties/Magnetizability.h" #include "properties/NMRShielding.h" +#include "properties/OrbitalEnergies.h" #include "properties/Polarizability.h" #include "properties/QuadrupoleMoment.h" #include "properties/SCFEnergy.h" @@ -60,6 +62,8 @@ namespace mrchem { +template using PropertyMap = std::map; + class Molecule final { public: explicit Molecule(int c = 0, int m = 1); @@ -75,7 +79,7 @@ class Molecule final { void setCharge(int c) { this->charge = c; } void setMultiplicity(int m) { this->multiplicity = m; } - void setGaugeOrigin(mrcpp::Coord<3> &o) { this->origin = o; } + void setGaugeOrigin(const mrcpp::Coord<3> &o) { this->origin = o; } mrcpp::Coord<3> getGaugeOrigin() const { return this->origin; } mrcpp::Coord<3> calcCenterOfMass() const; @@ -97,17 +101,26 @@ class Molecule final { auto getOrbitalsX_p() const { return this->orbitals_x; } auto getOrbitalsY_p() const { return this->orbitals_y; } + nlohmann::json json() const; void printGeometry() const; + void printEnergies(const std::string &txt) const; void printProperties() const; void initPerturbedOrbitals(bool dynamic); - SCFEnergy &getSCFEnergy(); - DipoleMoment &getDipoleMoment(); - QuadrupoleMoment &getQuadrupoleMoment(); - Magnetizability &getMagnetizability(); - NMRShielding &getNMRShielding(int k); - Polarizability &getPolarizability(double omega); + SCFEnergy &getSCFEnergy() { return this->energy; } + OrbitalEnergies &getOrbitalEnergies() { return this->epsilon; } + 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); } + Magnetizability &getMagnetizability(const std::string &id) { return this->magnetizability.at(id); } + NMRShielding &getNMRShielding(const std::string &id) { return this->nmr_shielding.at(id); } + + PropertyMap &getDipoleMoments() { return this->dipole; } + PropertyMap &getQuadrupoleMoments() { return this->quadrupole; } + PropertyMap &getPolarizabilities() { return this->polarizability; } + PropertyMap &getMagnetizabilities() { return this->magnetizability; } + PropertyMap &getNMRShieldings() { return this->nmr_shielding; } protected: int charge{0}; @@ -121,14 +134,14 @@ class Molecule final { std::shared_ptr orbitals_y{nullptr}; // Properties - std::unique_ptr energy{nullptr}; - std::unique_ptr dipole{nullptr}; - std::unique_ptr quadrupole{nullptr}; - std::unique_ptr magnetizability{nullptr}; - std::vector> polarizability{}; - std::vector> nmr{}; - - void initNuclearProperties(int nNucs); + SCFEnergy energy{}; + OrbitalEnergies epsilon{}; + PropertyMap dipole{}; + PropertyMap quadrupole{}; + PropertyMap polarizability{}; + PropertyMap magnetizability{}; + PropertyMap nmr_shielding{}; + void readCoordinateFile(const std::string &file); void readCoordinateString(const std::vector &coord_str); }; diff --git a/src/driver.cpp b/src/driver.cpp index b80953e36..31faa6620 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -91,9 +91,10 @@ namespace mrchem { namespace driver { DerivativeOperator_p get_derivative(const std::string &name); -template RankOneTensorOperator get_operator(const json &json_oper); -template RankTwoTensorOperator get_operator(const json &json_oper); +template RankOneTensorOperator get_operator(const std::string &name, const json &json_oper); +template RankTwoTensorOperator get_operator(const std::string &name, const json &json_oper); void build_fock_operator(const json &input, Molecule &mol, FockOperator &F, int order); +void init_properties(const json &json_prop, Molecule &mol); namespace scf { bool guess_orbitals(const json &input, Molecule &mol); @@ -118,37 +119,86 @@ void calc_properties(const json &input, Molecule &mol, int dir, double omega); void driver::init_molecule(const json &json_mol, Molecule &mol) { print_utils::headline(0, "Initializing Molecule"); - auto charge = json_mol["charge"].get(); - auto multiplicity = json_mol["multiplicity"].get(); - auto gauge_origin = json_mol["gauge_origin"].get>(); + auto charge = json_mol["charge"]; + auto multiplicity = json_mol["multiplicity"]; + auto gauge_origin = json_mol["gauge_origin"]; mol.setCharge(charge); mol.setMultiplicity(multiplicity); mol.setGaugeOrigin(gauge_origin); - Nuclei &nuclei = mol.getNuclei(); - for (const auto &coord : json_mol["coords"].get()) { - auto atom = coord["atom"].get(); - auto xyz = coord["xyz"].get>(); + auto &nuclei = mol.getNuclei(); + for (const auto &coord : json_mol["coords"]) { + auto atom = coord["atom"]; + auto xyz = coord["xyz"]; nuclei.push_back(atom, xyz); } mol.printGeometry(); } +void driver::init_properties(const json &json_prop, Molecule &mol) { + if (json_prop.contains("dipole_moment")) { + for (const auto &item : json_prop["dipole_moment"].items()) { + const auto &id = item.key(); + const auto &r_O = item.value()["r_O"]; + auto &dip_map = mol.getDipoleMoments(); + if (not dip_map.count(id)) dip_map.insert({id, DipoleMoment(r_O)}); + } + } + if (json_prop.contains("quadrupole_moment")) { + for (const auto &item : json_prop["quadrupole_moment"].items()) { + const auto &id = item.key(); + const auto &r_O = item.value()["r_O"]; + auto &quad_map = mol.getQuadrupoleMoments(); + if (not quad_map.count(id)) quad_map.insert({id, QuadrupoleMoment(r_O)}); + } + } + if (json_prop.contains("polarizability")) { + for (const auto &item : json_prop["polarizability"].items()) { + const auto &id = item.key(); + const auto &r_O = item.value()["r_O"]; + const auto &omega = item.value()["frequency"]; + auto &pol_map = mol.getPolarizabilities(); + if (not pol_map.count(id)) pol_map.insert({id, Polarizability(omega, r_O)}); + } + } + if (json_prop.contains("magnetizability")) { + for (const auto &item : json_prop["magnetizability"].items()) { + const auto &id = item.key(); + const auto &r_O = item.value()["r_O"]; + const auto &omega = item.value()["frequency"]; + auto &mag_map = mol.getMagnetizabilities(); + if (not mag_map.count(id)) mag_map.insert({id, Magnetizability(omega, r_O)}); + } + } + if (json_prop.contains("nmr_shielding")) { + for (const auto &item : json_prop["nmr_shielding"].items()) { + const auto &id = item.key(); + const auto &r_O = item.value()["r_O"]; + const auto &r_K = item.value()["r_K"]; + auto &nmr_map = mol.getNMRShieldings(); + if (not nmr_map.count(id)) nmr_map.insert({id, NMRShielding(r_K, r_O)}); + } + } +} + /** @brief Run ground-state SCF calculation * * This function will update the ground state orbitals and the Fock * matrix of the molecule based on the chosen electronic structure * method. The resulting orbitals will be either diagonalized or - * localized at exit. Returns true if the calculation converges. + * localized at exit. Returns a JSON record of the calculation. * * After convergence the requested ground-state properties are computed. * * This function expects the "scf_calculation" subsection of the input. */ -bool driver::scf::run(const json &json_scf, Molecule &mol) { +json driver::scf::run(const json &json_scf, Molecule &mol) { print_utils::headline(0, "Computing Ground State Wavefunction"); + json json_out = {{"success", true}}; + + if (json_scf.contains("properties")) driver::init_properties(json_scf["properties"], mol); /////////////////////////////////////////////////////////// //////////////// Building Fock Operator /////////////// @@ -163,10 +213,12 @@ bool driver::scf::run(const json &json_scf, Molecule &mol) { /////////////////////////////////////////////////////////// const auto &json_guess = json_scf["initial_guess"]; - if (driver::scf::guess_orbitals(json_guess, mol)) { - driver::scf::guess_energy(json_guess, mol, F); + if (scf::guess_orbitals(json_guess, mol)) { + scf::guess_energy(json_guess, mol, F); + json_out["initial_energy"] = mol.getSCFEnergy().json(); } else { - return false; + json_out["success"] = false; + return json_out; } /////////////////////////////////////////////////////////// @@ -174,21 +226,19 @@ bool driver::scf::run(const json &json_scf, Molecule &mol) { /////////////////////////////////////////////////////////// // Run GroundStateSolver if present in input JSON - auto success = true; - auto scf_solver = json_scf.find("scf_solver"); - if (scf_solver != json_scf.end()) { - auto kain = (*scf_solver)["kain"]; - auto method = (*scf_solver)["method"]; - auto max_iter = (*scf_solver)["max_iter"]; - auto rotation = (*scf_solver)["rotation"]; - auto localize = (*scf_solver)["localize"]; - auto file_chk = (*scf_solver)["file_chk"]; - auto checkpoint = (*scf_solver)["checkpoint"]; - auto start_prec = (*scf_solver)["start_prec"]; - auto final_prec = (*scf_solver)["final_prec"]; - auto energy_thrs = (*scf_solver)["energy_thrs"]; - auto orbital_thrs = (*scf_solver)["orbital_thrs"]; - auto helmholtz_prec = (*scf_solver)["helmholtz_prec"]; + if (json_scf.contains("scf_solver")) { + auto kain = json_scf["scf_solver"]["kain"]; + auto method = json_scf["scf_solver"]["method"]; + auto max_iter = json_scf["scf_solver"]["max_iter"]; + auto rotation = json_scf["scf_solver"]["rotation"]; + auto localize = json_scf["scf_solver"]["localize"]; + auto file_chk = json_scf["scf_solver"]["file_chk"]; + auto checkpoint = json_scf["scf_solver"]["checkpoint"]; + auto start_prec = json_scf["scf_solver"]["start_prec"]; + auto final_prec = json_scf["scf_solver"]["final_prec"]; + auto energy_thrs = json_scf["scf_solver"]["energy_thrs"]; + auto orbital_thrs = json_scf["scf_solver"]["orbital_thrs"]; + auto helmholtz_prec = json_scf["scf_solver"]["helmholtz_prec"]; GroundStateSolver solver; solver.setHistory(kain); @@ -201,25 +251,21 @@ bool driver::scf::run(const json &json_scf, Molecule &mol) { solver.setHelmholtzPrec(helmholtz_prec); solver.setOrbitalPrec(start_prec, final_prec); solver.setThreshold(orbital_thrs, energy_thrs); - success = solver.optimize(mol, F); + json_out["scf_solver"] = solver.optimize(mol, F); + json_out["success"] = json_out["scf_solver"]["converged"]; } /////////////////////////////////////////////////////////// ////////// Computing Ground State Properties ////////// /////////////////////////////////////////////////////////// - if (success) { - auto json_orbs = json_scf.find("write_orbitals"); - if (json_orbs != json_scf.end()) driver::scf::write_orbitals(*json_orbs, mol); - - auto json_prop = json_scf.find("properties"); - if (json_prop != json_scf.end()) driver::scf::calc_properties(*json_prop, mol); - - auto json_plot = json_scf.find("cube_plot"); - if (json_plot != json_scf.end()) driver::scf::plot_quantities(*json_plot, mol); + if (json_out["success"]) { + if (json_scf.contains("write_orbitals")) scf::write_orbitals(json_scf["write_orbitals"], mol); + if (json_scf.contains("properties")) scf::calc_properties(json_scf["properties"], mol); + if (json_scf.contains("cube_plot")) scf::plot_quantities(json_scf["cube_plot"], mol); } - return success; + return json_out; } /** @brief Run initial guess calculation for the orbitals @@ -301,9 +347,9 @@ bool driver::scf::guess_energy(const json &json_guess, Molecule &mol, FockOperat print_utils::text(0, "Localization ", (localize) ? "On" : "Off"); mrcpp::print::separator(0, '~', 2); - Timer timer; + Timer t_scf; auto plevel = Printer::getPrintLevel(); - if (plevel == 1) mrcpp::print::header(1, "Calculating Molecular Energy"); + if (plevel == 1) mrcpp::print::header(1, "Computing molecular energy"); auto &Phi = mol.getOrbitals(); auto &nucs = mol.getNuclei(); @@ -315,10 +361,19 @@ bool driver::scf::guess_energy(const json &json_guess, Molecule &mol, FockOperat F_mat = F(Phi, Phi); mol.getSCFEnergy() = F.trace(Phi, nucs); F.clear(); - if (plevel == 1) mrcpp::print::footer(1, timer, 2); + if (not localize) orbital::diagonalize(prec, Phi, F_mat); - mol.getSCFEnergy().print(); + if (plevel == 1) mrcpp::print::footer(1, t_scf, 2); + + Timer t_eps; + mrcpp::print::header(1, "Computing orbital energies"); + OrbitalEnergies &eps = mol.getOrbitalEnergies(); + eps.getOccupation() = orbital::get_occupations(Phi); + eps.getEpsilon() = orbital::calc_eigenvalues(Phi, F_mat); + eps.getSpin() = orbital::get_spins(Phi); + mrcpp::print::footer(1, t_eps, 2); + mol.printEnergies("initial"); return true; } @@ -338,81 +393,87 @@ void driver::scf::write_orbitals(const json &json_orbs, Molecule &mol) { void driver::scf::calc_properties(const json &json_prop, Molecule &mol) { Timer t_tot, t_lap; auto plevel = Printer::getPrintLevel(); - if (plevel == 1) mrcpp::print::header(1, "Computing Ground State Properties"); + if (plevel == 1) mrcpp::print::header(1, "Computing ground state properties"); - auto &nuclei = mol.getNuclei(); auto &Phi = mol.getOrbitals(); + auto &F_mat = mol.getFockMatrix(); + auto &nuclei = mol.getNuclei(); - auto json_dip = json_prop.find("dipole_moment"); - if (json_dip != json_prop.end()) { + if (json_prop.contains("dipole_moment")) { t_lap.start(); mrcpp::print::header(2, "Computing dipole moment"); - auto prec = (*json_dip)["precision"].get(); - DipoleMoment &mu = mol.getDipoleMoment(); - - auto h = driver::get_operator<3>(*json_dip); - h.setup(prec); - mu.getNuclear() = -h.trace(nuclei).real(); - mu.getElectronic() = h.trace(Phi).real(); - h.clear(); + for (const auto &item : json_prop["dipole_moment"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["operator"]; + auto h = driver::get_operator<3>(oper_name, item.value()); + h.setup(prec); + DipoleMoment &mu = mol.getDipoleMoment(id); + mu.getNuclear() = -h.trace(nuclei).real(); + mu.getElectronic() = h.trace(Phi).real(); + h.clear(); + } mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Dipole moment", t_lap); } - auto json_quad = json_prop.find("quadrupole_moment"); - if (json_quad != json_prop.end()) { + if (json_prop.contains("quadrupole_moment")) { t_lap.start(); mrcpp::print::header(2, "Computing quadrupole moment"); - auto prec = (*json_quad)["precision"].get(); - QuadrupoleMoment &Q = mol.getQuadrupoleMoment(); - - auto h = driver::get_operator<3, 3>(*json_quad); - h.setup(prec); - Q.getNuclear() = -h.trace(nuclei).real(); - Q.getElectronic() = h.trace(Phi).real(); - h.clear(); + for (const auto &item : json_prop["quadrupole_moment"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["operator"]; + auto h = driver::get_operator<3, 3>(oper_name, item.value()); + h.setup(prec); + QuadrupoleMoment &Q = mol.getQuadrupoleMoment(id); + Q.getNuclear() = -h.trace(nuclei).real(); + Q.getElectronic() = h.trace(Phi).real(); + h.clear(); + } mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Quadrupole moment", t_lap); } - auto json_hyperpolarizability = json_prop.find("hyperpolarizability"); - if (json_hyperpolarizability != json_prop.end()) MSG_ERROR("Hyperpolarizability not implemented"); - - auto json_gradient = json_prop.find("nuclear_gradient"); - if (json_gradient != json_prop.end()) MSG_ERROR("Nuclear gradient not implemented"); - - auto json_mag = json_prop.find("magnetizability"); - if (json_mag != json_prop.end()) { + if (json_prop.contains("magnetizability")) { t_lap.start(); mrcpp::print::header(2, "Computing magnetizability (dia)"); - auto prec = (*json_mag)["precision"].get(); - Magnetizability &xi = mol.getMagnetizability(); - - auto h = driver::get_operator<3, 3>(*json_mag); - h.setup(prec); - xi.getDiamagnetic() = -h.trace(Phi).real(); - h.clear(); + for (const auto &item : json_prop["magnetizability"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["dia_operator"]; + auto h = driver::get_operator<3, 3>(oper_name, item.value()); + h.setup(prec); + Magnetizability &xi = mol.getMagnetizability(id); + xi.getDiamagnetic() = -h.trace(Phi).real(); + h.clear(); + } mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Magnetizability (dia)", t_lap); } - auto json_nmr = json_prop.find("nmr_shielding"); - if (json_nmr != json_prop.end()) { + if (json_prop.contains("nmr_shielding")) { t_lap.start(); mrcpp::print::header(2, "Computing NMR shielding (dia)"); - for (const auto &json_nuc : *json_nmr) { - auto k = json_nuc["nucleus_k"].get(); - auto prec = json_nuc["precision"].get(); - NMRShielding &sigma_k = mol.getNMRShielding(k); - - auto h = driver::get_operator<3, 3>(json_nuc); + for (const auto &item : json_prop["nmr_shielding"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["dia_operator"]; + auto h = driver::get_operator<3, 3>(oper_name, item.value()); h.setup(prec); - sigma_k.getDiamagnetic() = h.trace(Phi).real(); + NMRShielding &sigma = mol.getNMRShielding(id); + sigma.getDiamagnetic() = h.trace(Phi).real(); h.clear(); } mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "NMR shielding (dia)", t_lap); } + + if (json_prop.contains("hyperpolarizability")) MSG_ERROR("Hyperpolarizability not implemented"); + if (json_prop.contains("geometric_derivative")) MSG_ERROR("Geometric derivative not implemented"); + if (json_prop.contains("hyperfine_coupling")) MSG_ERROR("Hyperfine coupling not implemented"); + if (json_prop.contains("spin_spin_coupling")) MSG_ERROR("Spin-spin coupling not implemented"); + if (plevel == 1) mrcpp::print::footer(1, t_tot, 2); } @@ -508,24 +569,31 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) { * based on the chosen electronic structure method and perturbation * operator. Each response calculation corresponds to one particular * perturbation operator (could be a vector operator with several - * components). Returns true if the calculation converges. + * components). Returns a JSON record of the calculation. * * After convergence the requested linear response properties are computed. * * This function expects a single subsection entry in the "rsp_calculations" * vector of the input. */ -bool driver::rsp::run(const json &json_rsp, Molecule &mol) { +json driver::rsp::run(const json &json_rsp, Molecule &mol) { print_utils::headline(0, "Computing Linear Response Wavefunction"); + json json_out = {{"success", true}}; + + if (json_rsp.contains("properties")) driver::init_properties(json_rsp["properties"], mol); /////////////////////////////////////////////////////////// ///////////// Preparing Unperturbed System //////////// /////////////////////////////////////////////////////////// + Timer t_unpert; + auto plevel = Printer::getPrintLevel(); + if (plevel == 1) mrcpp::print::header(1, "Preparing unperturbed system"); + const auto &json_unpert = json_rsp["unperturbed"]; const auto &unpert_fock = json_unpert["fock_operator"]; auto unpert_loc = json_unpert["localize"]; - auto unpert_prec = json_unpert["prec"]; + auto unpert_prec = json_unpert["precision"]; auto &Phi = mol.getOrbitals(); auto &F_mat = mol.getFockMatrix(); @@ -539,6 +607,9 @@ bool driver::rsp::run(const json &json_rsp, Molecule &mol) { FockOperator F_0; driver::build_fock_operator(unpert_fock, mol, F_0, 0); F_0.setup(unpert_prec); + if (plevel == 1) mrcpp::print::footer(1, t_unpert, 2); + + if (json_rsp.contains("properties")) scf::calc_properties(json_rsp["properties"], mol); /////////////////////////////////////////////////////////// ////////////// Preparing Perturbed System ///////////// @@ -553,10 +624,12 @@ bool driver::rsp::run(const json &json_rsp, Molecule &mol) { driver::build_fock_operator(json_fock_1, mol, F_1, 1); const auto &json_pert = json_rsp["perturbation"]; - auto h_1 = driver::get_operator<3>(json_pert); - - auto success = true; + auto h_1 = driver::get_operator<3>(json_pert["operator"], json_pert); + json_out["perturbation"] = json_pert["operator"]; + json_out["frequency"] = omega; + json_out["components"] = {}; for (auto d = 0; d < 3; d++) { + json comp_out = {}; const auto &json_comp = json_rsp["components"][d]; F_1.perturbation() = h_1[d]; @@ -565,26 +638,25 @@ bool driver::rsp::run(const json &json_rsp, Molecule &mol) { /////////////////////////////////////////////////////////// const auto &json_guess = json_comp["initial_guess"]; - success = driver::rsp::guess_orbitals(json_guess, mol); + json_out["success"] = rsp::guess_orbitals(json_guess, mol); /////////////////////////////////////////////////////////// ///////////// Optimizing Perturbed Orbitals //////////// /////////////////////////////////////////////////////////// - auto rsp_solver = json_comp.find("rsp_solver"); - if (rsp_solver != json_comp.end()) { - auto kain = (*rsp_solver)["kain"]; - auto method = (*rsp_solver)["method"]; - auto max_iter = (*rsp_solver)["max_iter"]; - auto file_chk_x = (*rsp_solver)["file_chk_x"]; - auto file_chk_y = (*rsp_solver)["file_chk_y"]; - auto checkpoint = (*rsp_solver)["checkpoint"]; - auto orth_prec = (*rsp_solver)["orth_prec"]; - auto start_prec = (*rsp_solver)["start_prec"]; - auto final_prec = (*rsp_solver)["final_prec"]; - auto orbital_thrs = (*rsp_solver)["orbital_thrs"]; - auto property_thrs = (*rsp_solver)["property_thrs"]; - auto helmholtz_prec = (*rsp_solver)["helmholtz_prec"]; + if (json_comp.contains("rsp_solver")) { + auto kain = json_comp["rsp_solver"]["kain"]; + auto method = json_comp["rsp_solver"]["method"]; + auto max_iter = json_comp["rsp_solver"]["max_iter"]; + auto file_chk_x = json_comp["rsp_solver"]["file_chk_x"]; + auto file_chk_y = json_comp["rsp_solver"]["file_chk_y"]; + auto checkpoint = json_comp["rsp_solver"]["checkpoint"]; + auto orth_prec = json_comp["rsp_solver"]["orth_prec"]; + auto start_prec = json_comp["rsp_solver"]["start_prec"]; + auto final_prec = json_comp["rsp_solver"]["final_prec"]; + auto orbital_thrs = json_comp["rsp_solver"]["orbital_thrs"]; + auto property_thrs = json_comp["rsp_solver"]["property_thrs"]; + auto helmholtz_prec = json_comp["rsp_solver"]["helmholtz_prec"]; LinearResponseSolver solver(dynamic); solver.setHistory(kain); @@ -597,28 +669,27 @@ bool driver::rsp::run(const json &json_rsp, Molecule &mol) { solver.setThreshold(orbital_thrs, property_thrs); solver.setOrthPrec(orth_prec); - success = solver.optimize(omega, mol, F_0, F_1); + comp_out["rsp_solver"] = solver.optimize(omega, mol, F_0, F_1); + json_out["success"] = comp_out["rsp_solver"]["converged"]; } /////////////////////////////////////////////////////////// //////////// Compute Response Properties ////////////// /////////////////////////////////////////////////////////// - if (success) { - auto json_orbs = json_comp.find("write_orbitals"); - if (json_orbs != json_comp.end()) driver::rsp::write_orbitals(*json_orbs, mol, dynamic); - - auto json_prop = json_rsp.find("properties"); - if (json_prop != json_rsp.end()) driver::rsp::calc_properties(*json_prop, mol, d, omega); + 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, d, omega); } mol.getOrbitalsX().clear(); // Clear orbital vector mol.getOrbitalsY().clear(); // Clear orbital vector + json_out["components"].push_back(comp_out); } F_0.clear(); mol.getOrbitalsX_p().reset(); // Release shared_ptr mol.getOrbitalsY_p().reset(); // Release shared_ptr - return success; + return json_out; } /** @brief Run initial guess calculation for the response orbitals @@ -632,8 +703,8 @@ bool driver::rsp::run(const json &json_rsp, Molecule &mol) { * This function expects the "initial_guess" subsection of the input. */ bool driver::rsp::guess_orbitals(const json &json_guess, Molecule &mol) { - auto prec = json_guess["prec"]; auto type = json_guess["type"]; + auto prec = json_guess["precision"]; auto mw_xp = json_guess["file_x_p"]; auto mw_xa = json_guess["file_x_a"]; auto mw_xb = json_guess["file_x_b"]; @@ -705,59 +776,68 @@ void driver::rsp::write_orbitals(const json &json_orbs, Molecule &mol, bool dyna void driver::rsp::calc_properties(const json &json_prop, Molecule &mol, int dir, double omega) { Timer t_tot, t_lap; auto plevel = Printer::getPrintLevel(); - if (plevel == 1) mrcpp::print::header(1, "Computing Linear Response Properties"); + if (plevel == 1) mrcpp::print::header(1, "Computing linear response properties"); auto &Phi = mol.getOrbitals(); auto &X = mol.getOrbitalsX(); auto &Y = mol.getOrbitalsY(); - auto json_pol = json_prop.find("polarizability"); - if (json_pol != json_prop.end()) { + if (json_prop.contains("polarizability")) { t_lap.start(); mrcpp::print::header(2, "Computing polarizability"); - auto prec = (*json_pol)["precision"].get(); - Polarizability &alpha = mol.getPolarizability(omega); - - auto h = driver::get_operator<3>(*json_pol); - h.setup(prec); - alpha.getTensor().row(dir) = -h.trace(Phi, X, Y).real(); - h.clear(); + for (const auto &item : json_prop["polarizability"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["operator"]; + auto h = driver::get_operator<3>(oper_name, item.value()); + h.setup(prec); + Polarizability &alpha = mol.getPolarizability(id); + alpha.getTensor().row(dir) = -h.trace(Phi, X, Y).real(); + h.clear(); + } mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Polarizability", t_lap); } - auto json_mag = json_prop.find("magnetizability"); - if (json_mag != json_prop.end()) { + if (json_prop.contains("magnetizability")) { t_lap.start(); mrcpp::print::header(2, "Computing magnetizability (para)"); - auto prec = (*json_mag)["precision"].get(); - Magnetizability &xi = mol.getMagnetizability(); - - auto h = driver::get_operator<3>(*json_mag); - h.setup(prec); - xi.getParamagnetic().row(dir) = -h.trace(Phi, X, Y).real(); - h.clear(); + for (const auto &item : json_prop["magnetizability"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["para_operator"]; + auto h = driver::get_operator<3>(oper_name, item.value()); + h.setup(prec); + Magnetizability &xi = mol.getMagnetizability(id); + xi.getParamagnetic().row(dir) = -h.trace(Phi, X, Y).real(); + h.clear(); + } mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Magnetizability (para)", t_lap); } - auto json_nmr = json_prop.find("nmr_shielding"); - if (json_nmr != json_prop.end()) { + if (json_prop.contains("nmr_shielding")) { t_lap.start(); mrcpp::print::header(2, "Computing NMR shielding (para)"); - for (const auto &json_nuc : *json_nmr) { - auto k = json_nuc["nucleus_k"].get(); - auto prec = json_nuc["precision"].get(); - NMRShielding &sigma_k = mol.getNMRShielding(k); - - auto h = driver::get_operator<3>(json_nuc); + for (const auto &item : json_prop["nmr_shielding"].items()) { + const auto &id = item.key(); + const auto &prec = item.value()["precision"]; + const auto &oper_name = item.value()["para_operator"]; + auto h = driver::get_operator<3>(oper_name, item.value()); h.setup(prec); - sigma_k.getParamagnetic().row(dir) = -h.trace(Phi, X, Y).real(); + NMRShielding &sigma = mol.getNMRShielding(id); + sigma.getParamagnetic().row(dir) = -h.trace(Phi, X, Y).real(); h.clear(); } 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("geometry_derivative")) MSG_ERROR("Geometry derivative not implemented"); + if (json_prop.contains("hyperfine_coupling")) MSG_ERROR("Hyperfine coupling not implemented"); + if (json_prop.contains("spin_spin_coupling")) MSG_ERROR("Spin-spin coupling not implemented"); + if (plevel == 1) mrcpp::print::footer(1, t_tot, 2); } @@ -776,9 +856,8 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera /////////////////////////////////////////////////////////// ////////////////// Kinetic Operator /////////////////// /////////////////////////////////////////////////////////// - auto json_kinetic = json_fock.find("kinetic_operator"); - if (json_kinetic != json_fock.end()) { - auto kin_diff = (*json_kinetic)["derivative"].get(); + if (json_fock.contains("kinetic_operator")) { + auto kin_diff = json_fock["kinetic_operator"]["derivative"]; auto D_p = driver::get_derivative(kin_diff); auto T_p = std::make_shared(D_p); F.getKineticOperator() = T_p; @@ -786,21 +865,19 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera /////////////////////////////////////////////////////////// ////////////////// Nuclear Operator /////////////////// /////////////////////////////////////////////////////////// - auto json_nuclear = json_fock.find("nuclear_operator"); - if (json_nuclear != json_fock.end()) { - auto proj_prec = (*json_nuclear)["proj_prec"].get(); - auto smooth_prec = (*json_nuclear)["smooth_prec"].get(); - auto shared_memory = (*json_nuclear)["shared_memory"].get(); + if (json_fock.contains("nuclear_operator")) { + auto proj_prec = json_fock["nuclear_operator"]["proj_prec"]; + auto smooth_prec = json_fock["nuclear_operator"]["smooth_prec"]; + auto shared_memory = json_fock["nuclear_operator"]["shared_memory"]; auto V_p = std::make_shared(nuclei, proj_prec, smooth_prec, shared_memory); F.getNuclearOperator() = V_p; } /////////////////////////////////////////////////////////// ////////////////// Coulomb Operator /////////////////// /////////////////////////////////////////////////////////// - auto json_coulomb = json_fock.find("coulomb_operator"); - if (json_coulomb != json_fock.end()) { - auto poisson_prec = (*json_coulomb)["poisson_prec"].get(); - auto shared_memory = (*json_coulomb)["shared_memory"].get(); + if (json_fock.contains("coulomb_operator")) { + auto poisson_prec = json_fock["coulomb_operator"]["poisson_prec"]; + auto shared_memory = json_fock["coulomb_operator"]["shared_memory"]; auto P_p = std::make_shared(*MRA, poisson_prec); if (order == 0) { auto J_p = std::make_shared(P_p, Phi_p, shared_memory); @@ -816,17 +893,16 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera //////////////////// XC Operator ////////////////////// /////////////////////////////////////////////////////////// double exx = 1.0; - auto json_xc = json_fock.find("xc_operator"); - if (json_xc != json_fock.end()) { - auto grid_prec = (*json_xc)["grid_prec"].get(); - auto shared_memory = (*json_xc)["shared_memory"].get(); - auto json_xcfunc = (*json_xc)["xc_functional"].get(); - auto xc_spin = json_xcfunc["spin"].get(); - auto xc_gamma = json_xcfunc["gamma"].get(); - auto xc_log_grad = json_xcfunc["log_grad"].get(); - auto xc_cutoff = json_xcfunc["cutoff"].get(); - auto xc_diff = json_xcfunc["derivative"].get(); - auto xc_funcs = json_xcfunc["functionals"].get(); + if (json_fock.contains("xc_operator")) { + auto grid_prec = json_fock["xc_operator"]["grid_prec"]; + auto shared_memory = json_fock["xc_operator"]["shared_memory"]; + auto json_xcfunc = json_fock["xc_operator"]["xc_functional"]; + auto xc_spin = json_xcfunc["spin"]; + auto xc_gamma = json_xcfunc["gamma"]; + auto xc_log_grad = json_xcfunc["log_grad"]; + auto xc_cutoff = json_xcfunc["cutoff"]; + auto xc_diff = json_xcfunc["derivative"]; + auto xc_funcs = json_xcfunc["functionals"]; auto xc_order = order + 1; mrdft::Factory xc_factory(*MRA); @@ -836,8 +912,8 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera xc_factory.setLogGradient(xc_log_grad); xc_factory.setDensityCutoff(xc_cutoff); for (const auto &f : xc_funcs) { - auto name = f["name"].get(); - auto coef = f["coef"].get(); + auto name = f["name"]; + auto coef = f["coef"]; xc_factory.setFunctional(name, coef); } auto mrdft_p = xc_factory.build(); @@ -856,10 +932,9 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera /////////////////////////////////////////////////////////// ///////////////// Exchange Operator /////////////////// /////////////////////////////////////////////////////////// - auto json_exchange = json_fock.find("exchange_operator"); - if (json_exchange != json_fock.end() and exx > mrcpp::MachineZero) { - auto poisson_prec = (*json_exchange)["poisson_prec"].get(); - auto screen_prec = (*json_exchange)["screen"].get(); + if (json_fock.contains("exchange_operator") and exx > mrcpp::MachineZero) { + auto poisson_prec = json_fock["exchange_operator"]["poisson_prec"]; + auto screen_prec = json_fock["exchange_operator"]["screen"]; auto P_p = std::make_shared(*MRA, poisson_prec); if (order == 0) { auto K_p = std::make_shared(P_p, Phi_p, screen_prec); @@ -872,10 +947,9 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera /////////////////////////////////////////////////////////// ///////////////// External Operator /////////////////// /////////////////////////////////////////////////////////// - auto json_external = json_fock.find("external_operator"); - if (json_external != json_fock.end()) { - auto field = (*json_external)["electric_field"].get>(); - auto r_O = (*json_external)["r_O"].get>(); + if (json_fock.contains("external_operator")) { + auto field = json_fock["external_operator"]["electric_field"]; + auto r_O = json_fock["external_operator"]["r_O"]; auto V_ext = std::make_shared(field, r_O); F.getExtOperator() = V_ext; } @@ -883,45 +957,33 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockOpera } /** @brief Construct perturbation operator based on input keyword */ -template RankOneTensorOperator driver::get_operator(const json &json_oper) { +template RankOneTensorOperator driver::get_operator(const std::string &name, const json &json_inp) { RankOneTensorOperator h; - auto oper = json_oper["operator"].get(); - if (oper == "h_e_dip") { - auto r_O = json_oper["r_O"].get>(); - h = H_E_dip(r_O); - } else if (oper == "h_b_dip") { - auto r_O = json_oper["r_O"].get>(); - auto pert_diff = json_oper["derivative"].get(); - auto D = driver::get_derivative(pert_diff); - h = H_B_dip(D, r_O); - } else if (oper == "h_m_pso") { - auto r_K = json_oper["r_K"].get>(); - auto smoothing = json_oper["smoothing"].get(); - auto pert_diff = json_oper["derivative"].get(); - auto D = driver::get_derivative(pert_diff); - h = H_M_pso(D, r_K, smoothing); + if (name == "h_e_dip") { + h = H_E_dip(json_inp["r_O"]); + } else if (name == "h_b_dip") { + auto D = driver::get_derivative(json_inp["derivative"]); + h = H_B_dip(D, json_inp["r_O"]); + } else if (name == "h_m_pso") { + auto D = driver::get_derivative(json_inp["derivative"]); + h = H_M_pso(D, json_inp["r_K"], json_inp["smoothing"]); } else { - MSG_ERROR("Invalid operator: " << oper); + MSG_ERROR("Invalid operator: " << name); } return h; } -template RankTwoTensorOperator driver::get_operator(const json &json_oper) { +template +RankTwoTensorOperator driver::get_operator(const std::string &name, const json &json_inp) { RankTwoTensorOperator h; - auto oper = json_oper["operator"].get(); - if (oper == "h_e_quad") { - auto r_O = json_oper["r_O"].get>(); - h = H_E_quad(r_O); - } else if (oper == "h_bb_dia") { - auto r_O = json_oper["r_O"].get>(); - h = H_BB_dia(r_O); - } else if (oper == "h_bm_dia") { - auto r_O = json_oper["r_O"].get>(); - auto r_K = json_oper["r_K"].get>(); - auto smoothing = json_oper["smoothing"].get(); - h = H_BM_dia(r_O, r_K, smoothing); + if (name == "h_e_quad") { + h = H_E_quad(json_inp["r_O"]); + } else if (name == "h_bb_dia") { + h = H_BB_dia(json_inp["r_O"]); + } else if (name == "h_bm_dia") { + h = H_BM_dia(json_inp["r_O"], json_inp["r_K"], json_inp["smoothing"]); } else { - MSG_ERROR("Invalid operator: " << oper); + MSG_ERROR("Invalid operator: " << name); } return h; } @@ -943,10 +1005,12 @@ DerivativeOperator_p driver::get_derivative(const std::string &name) { return D; } -void driver::print_properties(const Molecule &mol) { +json driver::print_properties(const Molecule &mol) { print_utils::headline(0, "Printing Molecular Properties"); mol.printGeometry(); + mol.printEnergies("final"); mol.printProperties(); + return mol.json(); } } // namespace mrchem diff --git a/src/driver.h b/src/driver.h index f20c3f1c5..de6967990 100644 --- a/src/driver.h +++ b/src/driver.h @@ -32,15 +32,14 @@ class Molecule; namespace driver { void init_molecule(const nlohmann::json &input, Molecule &mol); -void print_properties(const Molecule &mol); +nlohmann::json print_properties(const Molecule &mol); namespace scf { -bool run(const nlohmann::json &input, Molecule &mol); -} // namespace scf - +nlohmann::json run(const nlohmann::json &input, Molecule &mol); +} namespace rsp { -bool run(const nlohmann::json &input, Molecule &mol); -} // namespace rsp +nlohmann::json run(const nlohmann::json &input, Molecule &mol); +} } // namespace driver } // namespace mrchem diff --git a/src/input/mrchem.in b/src/input/mrchem.in index 2d250d1e2..ce5191f85 100644 --- a/src/input/mrchem.in +++ b/src/input/mrchem.in @@ -71,16 +71,17 @@ def main(): # now that all keywords have sensible values, # we can translate user input into program input - program_dict = translate_input(user_dict) + program_dict = {} + program_dict["input"] = translate_input(user_dict) inp_name, ext_ext = os.path.splitext(inp_file_cmd) - xfile = inp_name + '_parsed.json' + xfile = inp_name + '.json' with open(xfile, 'w') as fd: fd.write(json.dumps(program_dict, indent=2)) if not dryrun: cmd = executable + ' ' + xfile - if program_dict["printer"]["print_input"]: + if program_dict["input"]["printer"]["print_input"]: subprocess.call('cat ' + inp_file_cmd, shell=True) subprocess.call(cmd, shell=True) @@ -228,8 +229,8 @@ def write_mra(user_dict, mol_dict): max_scale = 30 + min_scale mra_dict = { - "basis_type": user_dict["Basis"]["type"][0].lower(), - "order": order, + "basis_type": user_dict["Basis"]["type"].lower(), + "basis_order": order, "boxes": [2, 2, 2], "corner": [-1, -1, -1], "min_scale": min_scale, @@ -484,38 +485,19 @@ def write_scf_solver(user_dict, method_name): def write_scf_properties(user_dict, mol_dict): prop_dict = {} if user_dict["Properties"]["dipole_moment"]: - prop_dict["dipole_moment"] = { - "operator": "h_e_dip", + prop_dict["dipole_moment"] = {} + prop_dict["dipole_moment"]['dip-1'] = { "precision": user_dict["world_prec"], + "operator": "h_e_dip", "r_O": mol_dict["gauge_origin"] } if user_dict["Properties"]["quadrupole_moment"]: - prop_dict["quadrupole_moment"] = { - "operator": "h_e_quad", - "precision": user_dict["world_prec"], - "r_O": mol_dict["gauge_origin"] - } - if user_dict["Properties"]["magnetizability"]: - prop_dict["magnetizability"] = { - "operator": "h_bb_dia", + prop_dict["quadrupole_moment"] = {} + prop_dict["quadrupole_moment"]['quad-1'] = { "precision": user_dict["world_prec"], + "operator": "h_e_quad", "r_O": mol_dict["gauge_origin"] } - if user_dict["Properties"]["nmr_shielding"]: - nuc_vec = user_dict["NMRShielding"]["nucleus_k"] - prop_dict["nmr_shielding"] = [] - all_nucs = (nuc_vec[0] < 0) - nuclei = mol_dict["coords"] - for k in range(len(nuclei)): - if all_nucs or k in nuc_vec: - prop_dict["nmr_shielding"].append({ - "nucleus_k": k, - "operator": "h_bm_dia", - "precision": user_dict["world_prec"], - "smoothing": user_dict["world_prec"], - "r_O": mol_dict["gauge_origin"], - "r_K": nuclei[k]["xyz"] - }) return prop_dict @@ -533,61 +515,73 @@ def write_scf_plot(user_dict): ############################################################ def write_rsp_calculations(user_dict, mol_dict): - rsp_dict = [] + rsp_dict = {} run_pol = user_dict["Properties"]["polarizability"] - pol_freq = user_dict["Polarizability"]["frequency"] run_mag = user_dict["Properties"]["magnetizability"] run_nmr = user_dict["Properties"]["nmr_shielding"] nuc_spec = user_dict["NMRShielding"]["nuclear_specific"] if run_pol: - for omega in pol_freq: + for omega in user_dict["Polarizability"]["frequency"]: + freq_key = f'{omega:6f}' + pol_key = 'pol-' + freq_key rsp_calc = write_rsp_calc(omega, user_dict, mol_dict) rsp_calc["perturbation"] = { "operator": "h_e_dip", "r_O": mol_dict["gauge_origin"] } - prop_dict = {} - prop_dict["polarizability"] = { - "operator": "h_e_dip", + rsp_calc["properties"] = {} + rsp_calc["properties"]["polarizability"] = {} + rsp_calc["properties"]["polarizability"][pol_key] = { + "frequency": omega, "precision": user_dict["world_prec"], + "operator": "h_e_dip", "r_O": mol_dict["gauge_origin"] } - rsp_calc["properties"] = prop_dict - rsp_dict.append(rsp_calc) + rsp_key = "ext_el-" + freq_key + rsp_dict[rsp_key] = rsp_calc if run_mag or (run_nmr and not nuc_spec): + omega = 0.0 # only static magnetic response + freq_key = f'{omega:6f}' + mag_key = 'mag-' + freq_key rsp_calc = write_rsp_calc(0.0, user_dict, mol_dict) rsp_calc["perturbation"] = { "operator": "h_b_dip", "derivative": user_dict["Derivatives"]["h_b_dip"], "r_O": mol_dict["gauge_origin"] } - prop_dict = {} + rsp_calc["properties"] = {} if run_mag: - prop_dict["magnetizability"] = { - "operator": "h_b_dip", + rsp_calc["properties"]["magnetizability"] = {} + rsp_calc["properties"]["magnetizability"][mag_key] = { + "frequency": omega, "precision": user_dict["world_prec"], + "dia_operator": "h_bb_dia", + "para_operator": "h_b_dip", "derivative": user_dict["Derivatives"]["h_b_dip"], "r_O": mol_dict["gauge_origin"] } if (run_nmr and not nuc_spec): - prop_dict["nmr_shielding"] = [] + rsp_calc["properties"]["nmr_shielding"] = {} nuc_vec = user_dict["NMRShielding"]["nucleus_k"] all_nucs = (nuc_vec[0] < 0) nuclei = mol_dict["coords"] for k in range(len(nuclei)): if all_nucs or k in nuc_vec: - prop_dict["nmr_shielding"].append({ - "nucleus_k": k, - "operator": "h_m_pso", + atom_key = str(k) + nuclei[k]["atom"] + nmr_key = 'nmr-' + atom_key + rsp_calc["properties"]["nmr_shielding"][nmr_key] = { "precision": user_dict["world_prec"], + "dia_operator": "h_bm_dia", + "para_operator": "h_m_pso", "smoothing": user_dict["world_prec"], - "derivative": user_dict["Derivatives"]["h_b_dip"], + "derivative": user_dict["Derivatives"]["h_m_pso"], + "r_O": mol_dict["gauge_origin"], "r_K": nuclei[k]["xyz"] - }) - rsp_calc["properties"] = prop_dict - rsp_dict.append(rsp_calc) + } + rsp_key = "ext_mag-" + freq_key + rsp_dict[rsp_key] = rsp_calc if run_nmr and nuc_spec: nuc_vec = user_dict["NMRShielding"]["nucleus_k"] @@ -595,6 +589,8 @@ def write_rsp_calculations(user_dict, mol_dict): nuclei = mol_dict["coords"] for k in range(len(nuclei)): if (all_nucs) or (k in nuc_vec): + atom_key = str(k) + nuclei[k]["atom"] + nmr_key = 'nmr-' + atom_key rsp_calc = write_rsp_calc(0.0, user_dict, mol_dict) rsp_calc["perturbation"] = { "operator": "h_m_pso", @@ -602,17 +598,19 @@ def write_rsp_calculations(user_dict, mol_dict): "derivative": user_dict["Derivatives"]["h_m_pso"], "r_K": nuclei[k]["xyz"] } - prop_dict = {} - prop_dict["nmr_shielding"] = [] - prop_dict["nmr_shielding"].append({ - "nucleus_k": k, - "operator": "h_b_dip", + rsp_calc["properties"] = {} + rsp_calc["properties"]["nmr_shielding"] = {} + rsp_calc["properties"]["nmr_shielding"][nmr_key] = { "precision": user_dict["world_prec"], + "dia_operator": "h_bm_dia", + "para_operator": "h_b_dip", + "smoothing": user_dict["world_prec"], "derivative": user_dict["Derivatives"]["h_b_dip"], - "r_O": mol_dict["gauge_origin"] - }) - rsp_calc["properties"] = prop_dict - rsp_dict.append(rsp_calc) + "r_O": mol_dict["gauge_origin"], + "r_K": nuclei[k]["xyz"] + } + rsp_key = "nuc_mag-" + atom_key + rsp_dict[rsp_key] = rsp_calc return rsp_dict @@ -628,7 +626,7 @@ def write_rsp_calc(omega, user_dict, mol_dict): rsp_calc["dynamic"] = (omega > 1.0e-12) rsp_calc["fock_operator"] = write_rsp_fock(user_dict, mol_dict, wf_method, dft_funcs) rsp_calc["unperturbed"] = { - "prec": user_dict["world_prec"], + "precision": user_dict["world_prec"], "localize": rsp_dict["localize"], "fock_operator": write_scf_fock(user_dict, mol_dict, wf_method, dft_funcs) } @@ -637,14 +635,10 @@ def write_rsp_calc(omega, user_dict, mol_dict): for d in [0, 1, 2]: rsp_comp = {} rsp_comp["initial_guess"] = { - "prec": rsp_dict["guess_prec"], "type": rsp_dict["guess_type"].lower(), + "precision": rsp_dict["guess_prec"], "file_chk_x": rsp_dict["path_checkpoint"] + "/X_rsp_" + str(d), "file_chk_y": rsp_dict["path_checkpoint"] + "/Y_rsp_" + str(d), - "file_basis": file_dict["guess_basis"], - "file_gto_p": file_dict["guess_gto_p"], - "file_gto_a": file_dict["guess_gto_a"], - "file_gto_b": file_dict["guess_gto_b"], "file_x_p": file_dict["guess_x_p"] + "_rsp_" + str(d), "file_x_a": file_dict["guess_x_a"] + "_rsp_" + str(d), "file_x_b": file_dict["guess_x_b"] + "_rsp_" + str(d), diff --git a/src/input/template.yml b/src/input/template.yml index 41a953803..eabe1b217 100644 --- a/src/input/template.yml +++ b/src/input/template.yml @@ -98,10 +98,9 @@ sections: - name: print_level type: int default: 0 - predicates: - - 'value >= 0' docstring: | - Level of detail in the written output. Level 0 for production calculations. + Level of detail in the written output. Level 0 for production calculations, + negative level for complete silence. - name: print_input type: bool default: false @@ -190,7 +189,7 @@ sections: type: str default: 'interpolating' predicates: - - "value[0].lower() in ['i', 'l']" + - "value.lower() in ['interpolating', 'legendre']" docstring: | Polynomial type of multiwavelet basis. - name: Derivatives @@ -329,7 +328,7 @@ sections: keywords: - name: dipole_moment type: bool - default: false + default: true docstring: | Compute dipole moment. - name: quadrupole_moment @@ -640,7 +639,7 @@ sections: Length of KAIN iterative history. - name: localize type: bool - default: false + default: "user['SCF']['localize']" docstring: | Use canonical or localized unperturbed orbitals. - name: orbital_thrs diff --git a/src/mrchem.cpp b/src/mrchem.cpp index 695134e31..37d9861c4 100644 --- a/src/mrchem.cpp +++ b/src/mrchem.cpp @@ -1,15 +1,29 @@ -/** \mainpage The MRChem main program +/* + * 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) 2020 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. * - * \author Stig Rune Jensen + * This file is part of MRChem. * - * \version 1.0 + * 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. * - * \par Copyright: - * GPLv4 + * 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 . + * + * For information on the complete list of contributors to MRChem, see: + * */ -#include "MRCPP/Timer" -#include "parallel.h" + +#include #include "driver.h" #include "mrchem.h" @@ -26,22 +40,28 @@ using Timer = mrcpp::Timer; using namespace mrchem; int main(int argc, char **argv) { - const auto json_input = mrenv::fetch_input(argc, argv); + const auto json_inp = mrenv::fetch_json(argc, argv); - mrenv::initialize(json_input); - const auto &json_mol = json_input["molecule"]; - const auto &json_scf = json_input["scf_calculation"]; - const auto &json_rsps = json_input["rsp_calculations"]; + mrenv::initialize(json_inp); + const auto &mol_inp = json_inp["molecule"]; + const auto &scf_inp = json_inp["scf_calculation"]; + const auto &rsp_inp = json_inp["rsp_calculations"]; Timer timer; Molecule mol; - driver::init_molecule(json_mol, mol); - if (driver::scf::run(json_scf, mol)) { - for (const auto &json_rsp : json_rsps) driver::rsp::run(json_rsp, mol); + driver::init_molecule(mol_inp, mol); + auto scf_out = driver::scf::run(scf_inp, mol); + json rsp_out = {}; + if (scf_out["success"]) { + for (auto &i : rsp_inp.items()) rsp_out[i.key()] = driver::rsp::run(i.value(), mol); } - driver::print_properties(mol); mpi::barrier(mpi::comm_orb); + json json_out; + json_out["scf_calculation"] = scf_out; + json_out["rsp_calculations"] = rsp_out; + json_out["properties"] = driver::print_properties(mol); mrenv::finalize(timer.elapsed()); + mrenv::dump_json(json_inp, json_out); mpi::finalize(); return EXIT_SUCCESS; diff --git a/src/mrchem.h b/src/mrchem.h index 9e1d2fa72..bd9dd65d4 100644 --- a/src/mrchem.h +++ b/src/mrchem.h @@ -1,10 +1,26 @@ /* - * \date Feb 22, 2018 - * \author Stig Rune Jensen \n - * Hylleraas Centre for Quantum Molecular Sciences \n - * UiT - The Arctic University of Norway + * 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) 2020 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. * - * \breif Global objects and variables + * 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 . + * + * For information on the complete list of contributors to MRChem, see: + * */ #pragma once diff --git a/src/mrenum.h b/src/mrenum.h index ee83607c1..196d01980 100644 --- a/src/mrenum.h +++ b/src/mrenum.h @@ -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) 2020 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + /* * \date Feb 22, 2018 * \author Luca Frediani \n diff --git a/src/mrenv.cpp b/src/mrenv.cpp index 23cf30f3b..66f842cf9 100644 --- a/src/mrenv.cpp +++ b/src/mrenv.cpp @@ -1,9 +1,32 @@ -#include "MRCPP/Printer" - +/* + * 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) 2020 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#include +#include #include -#include "XCFun/xcfun.h" - #include "mrchem.h" #include "mrenv.h" #include "parallel.h" @@ -20,10 +43,9 @@ void init_printer(const json &json_print); void init_mra(const json &json_mra); void init_mpi(const json &json_mpi); void print_header(); -void print_footer(double wt); } // namespace mrenv -json mrenv::fetch_input(int argc, char **argv) { +json mrenv::fetch_json(int argc, char **argv) { const char *infile = nullptr; if (argc == 1) { infile = "STDIN"; @@ -39,29 +61,29 @@ json mrenv::fetch_input(int argc, char **argv) { ifs >> input; ifs.close(); - return input; + return input["input"]; } -void mrenv::initialize(const json &input) { - auto json_print = input.find("printer"); - auto json_mra = input.find("mra"); - auto json_mpi = input.find("mpi"); +void mrenv::initialize(const json &json_inp) { + auto json_print = json_inp.find("printer"); + auto json_mra = json_inp.find("mra"); + auto json_mpi = json_inp.find("mpi"); - if (json_mra == input.end()) MSG_ABORT("Missing MRA input!"); + if (json_mra == json_inp.end()) MSG_ABORT("Missing MRA input!"); - if (json_mra != input.end()) mrenv::init_mra(*json_mra); - if (json_mpi != input.end()) mrenv::init_mpi(*json_mpi); - if (json_print != input.end()) mrenv::init_printer(*json_print); + if (json_mra != json_inp.end()) mrenv::init_mra(*json_mra); + if (json_mpi != json_inp.end()) mrenv::init_mpi(*json_mpi); + if (json_print != json_inp.end()) mrenv::init_printer(*json_print); mrenv::print_header(); } void mrenv::init_printer(const json &json_print) { // Initialize printing - auto print_level = json_print["print_level"].get(); - auto print_prec = json_print["print_prec"].get(); - auto print_width = json_print["print_width"].get(); - auto print_mpi = json_print["print_mpi"].get(); + auto print_level = json_print["print_level"]; + auto print_prec = json_print["print_prec"]; + auto print_width = json_print["print_width"]; + auto print_mpi = json_print["print_mpi"]; auto fname = json_print["file_name"].get(); if (print_mpi) { Printer::init(print_level, mpi::world_rank, mpi::world_size, fname.c_str()); @@ -74,16 +96,16 @@ void mrenv::init_printer(const json &json_print) { void mrenv::init_mra(const json &json_mra) { // Initialize world box - auto min_scale = json_mra["min_scale"].get(); - auto max_scale = json_mra["max_scale"].get(); - auto corner = json_mra["corner"].get>(); - auto boxes = json_mra["boxes"].get>(); - auto sfac = json_mra["scaling_factor"].get>(); + int min_scale = json_mra["min_scale"]; + int max_scale = json_mra["max_scale"]; + auto corner = json_mra["corner"]; + auto boxes = json_mra["boxes"]; + auto sfac = json_mra["scaling_factor"]; mrcpp::BoundingBox<3> world(min_scale, corner, boxes, sfac); // Initialize scaling basis - auto order = json_mra["order"].get(); - auto btype = json_mra["basis_type"].get(); + auto order = json_mra["basis_order"]; + auto btype = json_mra["basis_type"]; auto max_depth = max_scale - min_scale; if (min_scale < mrcpp::MinScale) MSG_ABORT("Root scale too large"); @@ -91,10 +113,10 @@ void mrenv::init_mra(const json &json_mra) { if (max_depth > mrcpp::MaxDepth) MSG_ABORT("Max depth too large"); // Initialize global MRA - if (btype == "i") { + if (btype == "interpolating") { mrcpp::InterpolatingBasis basis(order); MRA = new mrcpp::MultiResolutionAnalysis<3>(world, basis, max_depth); - } else if (btype == "l") { + } else if (btype == "legendre") { mrcpp::LegendreBasis basis(order); MRA = new mrcpp::MultiResolutionAnalysis<3>(world, basis, max_depth); } else { @@ -103,9 +125,9 @@ void mrenv::init_mra(const json &json_mra) { } void mrenv::init_mpi(const json &json_mpi) { - mpi::numerically_exact = json_mpi["numerically_exact"].get(); - mpi::shared_memory_size = json_mpi["shared_memory_size"].get(); - mpi::bank_size = json_mpi["bank_size"].get(); + mpi::numerically_exact = json_mpi["numerically_exact"]; + mpi::shared_memory_size = json_mpi["shared_memory_size"]; + mpi::bank_size = json_mpi["bank_size"]; mpi::initialize(); // NB: must be after bank_size and init_mra but before init_printer and print_header } @@ -189,6 +211,7 @@ void mrenv::finalize(double wt) { auto hr = static_cast(wt / 3600.0); auto min = static_cast(std::fmod(wt, 3600.0) / 60.0); auto sec = static_cast(std::fmod(wt, 60.0)); + std::stringstream o_time; o_time << "Wall time : " << std::setw(2) << hr << "h" << std::setw(3) << min << "m" << std::setw(3) << sec << "s"; @@ -205,4 +228,17 @@ void mrenv::finalize(double wt) { mrcpp::print::separator(0, ' '); } +void mrenv::dump_json(const json &json_inp, const json &json_out) { + json json_tot; + json_tot["input"] = json_inp; + json_tot["output"] = json_out; + + if (mpi::grand_master()) { + std::ofstream ofs; + ofs.open("mrchem.json", std::ios::out); + ofs << json_tot.dump(2) << std::endl; + ofs.close(); + } +} + } // namespace mrchem diff --git a/src/mrenv.h b/src/mrenv.h index 7c48f8551..a77baf6e5 100644 --- a/src/mrenv.h +++ b/src/mrenv.h @@ -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) 2020 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + #pragma once #include @@ -5,10 +30,11 @@ namespace mrchem { namespace mrenv { -nlohmann::json fetch_input(int argc, char **argv); -void initialize(const nlohmann::json &input); +nlohmann::json fetch_json(int argc, char **argv); +void dump_json(const nlohmann::json &json_inp, const nlohmann::json &json_out); + +void initialize(const nlohmann::json &json_inp); void finalize(double wt); -void init_mpi(const nlohmann::json &json_mpi); } // namespace mrenv } // namespace mrchem diff --git a/src/parallel.cpp b/src/parallel.cpp index 12aacec38..8ea6602f8 100644 --- a/src/parallel.cpp +++ b/src/parallel.cpp @@ -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) 2020 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + #include #include "parallel.h" diff --git a/src/parallel.h b/src/parallel.h index 73939ebbd..95f7acbd8 100644 --- a/src/parallel.h +++ b/src/parallel.h @@ -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) 2020 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + #pragma once #include "MRCPP/Parallel" diff --git a/src/properties/DipoleMoment.h b/src/properties/DipoleMoment.h index 5d3e54c94..3002ea672 100644 --- a/src/properties/DipoleMoment.h +++ b/src/properties/DipoleMoment.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "mrchem.h" #include "utils/math_utils.h" @@ -35,13 +37,18 @@ namespace mrchem { // clang-format off class DipoleMoment final { public: + explicit DipoleMoment(const mrcpp::Coord<3> &o = {}) : r_O(o) {} + + void setOrigin(const mrcpp::Coord<3> &o) { this->r_O = o; } + const mrcpp::Coord<3> &getOrigin() const { return this->r_O; } + DoubleVector getTensor() const { return getNuclear() + getElectronic(); } DoubleVector &getNuclear() { return this->nuc_tensor; } DoubleVector &getElectronic() { return this->el_tensor; } const DoubleVector &getNuclear() const { return this->nuc_tensor; } const DoubleVector &getElectronic() const { return this->el_tensor; } - void print() const { + void print(const std::string &id) const { auto el_au = getElectronic().norm(); auto nuc_au = getNuclear().norm(); auto tot_au = getTensor().norm(); @@ -50,7 +57,9 @@ class DipoleMoment final { auto nuc_db = nuc_au * PHYSCONST::Debye; auto tot_db = tot_au * PHYSCONST::Debye; - mrcpp::print::header(0, "Dipole Moment"); + mrcpp::print::header(0, "Dipole Moment (" + id + ")"); + print_utils::coord(0, "r_O", getOrigin()); + mrcpp::print::separator(0, '-'); print_utils::vector(0, "Electronic vector", getElectronic()); print_utils::scalar(0, "Magnitude", el_au, "(au)"); print_utils::scalar(0, " ", el_db, "(Debye)"); @@ -65,7 +74,18 @@ class DipoleMoment final { mrcpp::print::separator(0, '=', 2); } + nlohmann::json json() const { + return { + {"r_O", getOrigin()}, + {"vector_nuc", print_utils::eigen_to_vector(getNuclear(), 1.0e-12)}, + {"vector_el", print_utils::eigen_to_vector(getElectronic(), 1.0e-12)}, + {"vector", print_utils::eigen_to_vector(getTensor(), 1.0e-12)}, + {"magnitude", getTensor().norm() } + }; + } + private: + mrcpp::Coord<3> r_O; DoubleVector nuc_tensor{math_utils::init_nan(3)}; DoubleVector el_tensor{math_utils::init_nan(3)}; }; diff --git a/src/properties/Magnetizability.h b/src/properties/Magnetizability.h index 3d9364fdf..7cd74a9b8 100644 --- a/src/properties/Magnetizability.h +++ b/src/properties/Magnetizability.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "mrchem.h" #include "utils/math_utils.h" @@ -35,13 +37,26 @@ namespace mrchem { // clang-format off class Magnetizability final { public: + explicit Magnetizability(double omega = 0.0, const mrcpp::Coord<3> &o = {}) : frequency(omega), r_O(o) {} + + void setFrequency(double omega) { this->frequency = omega; } + void setOrigin(const mrcpp::Coord<3> &o) { this->r_O = o; } + + double getFrequency() const { return this->frequency; } + const mrcpp::Coord<3> &getOrigin() const { return this->r_O; } + DoubleMatrix getTensor() const { return getDiamagnetic() + getParamagnetic(); } DoubleMatrix &getDiamagnetic() { return this->dia_tensor; } DoubleMatrix &getParamagnetic() { return this->para_tensor; } const DoubleMatrix &getDiamagnetic() const { return this->dia_tensor; } const DoubleMatrix &getParamagnetic() const { return this->para_tensor; } - void print() const { + void print(const std::string &id) const { + auto w_au = getFrequency(); + auto w_cm = PHYSCONST::cm_m1 * w_au; + auto dynamic = (w_au > mrcpp::MachineZero); + auto l_nm = (dynamic) ? (1.0e7 / w_cm) : 0.0; + auto iso_au_d = getDiamagnetic().trace() / 3.0; auto iso_au_p = getParamagnetic().trace() / 3.0; auto iso_au_t = iso_au_d + iso_au_p; @@ -51,7 +66,12 @@ class Magnetizability final { auto iso_si_d = iso_au_d * PHYSCONST::JT_m2; auto iso_si_p = iso_au_p * PHYSCONST::JT_m2; - mrcpp::print::header(0, "Magnetizability"); + mrcpp::print::header(0, "Magnetizability (" + id + ")"); + if (dynamic) print_utils::scalar(0, "Wavelength", l_nm, "(nm)"); + print_utils::scalar(0, "Frequency", w_au, "(au)"); + print_utils::scalar(0, " ", w_cm, "(cm-1)"); + print_utils::coord(0, "r_O", getOrigin()); + mrcpp::print::separator(0, '-'); print_utils::matrix(0, "Diamagnetic", getDiamagnetic()); print_utils::scalar(0, "Isotropic average", iso_au_d, "(au)"); print_utils::scalar(0, " ", iso_si_d, "(SI)"); @@ -66,7 +86,20 @@ class Magnetizability final { mrcpp::print::separator(0, '=', 2); } + nlohmann::json json() const { + return { + {"r_O", getOrigin()}, + {"frequency", getFrequency()}, + {"tensor_dia", print_utils::eigen_to_vector(getDiamagnetic(), 1.0e-12)}, + {"tensor_para", print_utils::eigen_to_vector(getParamagnetic(), 1.0e-12)}, + {"tensor", print_utils::eigen_to_vector(getTensor(), 1.0e-12)}, + {"isotropic_average", getTensor().trace() / 3.0 } + }; + } + private: + double frequency; + mrcpp::Coord<3> r_O; DoubleMatrix dia_tensor{math_utils::init_nan(3,3)}; DoubleMatrix para_tensor{math_utils::init_nan(3,3)}; }; diff --git a/src/properties/NMRShielding.h b/src/properties/NMRShielding.h index f4de8f1d6..c61fca0b0 100644 --- a/src/properties/NMRShielding.h +++ b/src/properties/NMRShielding.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "mrchem.h" #include "utils/math_utils.h" @@ -35,10 +37,13 @@ namespace mrchem { // clang-format off class NMRShielding final { public: - NMRShielding(int k, const Nucleus &n) : K(k), nuc(n) {} + explicit NMRShielding(const mrcpp::Coord<3> &k = {}, const mrcpp::Coord<3> &o = {}) : r_K(k), r_O(o) {} + + void setCoordK(const mrcpp::Coord<3> &k) { this->r_K = k; } + void setOrigin(const mrcpp::Coord<3> &o) { this->r_O = o; } - int getK() const { return this->K; } - const Nucleus &getNucleus() const { return this->nuc; } + const mrcpp::Coord<3> &getCoordK() const { return this->r_K; } + const mrcpp::Coord<3> &getOrigin() const { return this->r_O; } DoubleMatrix getTensor() const { return getDiamagnetic() + getParamagnetic(); } DoubleMatrix &getDiamagnetic() { return this->dia_tensor; } @@ -46,28 +51,54 @@ class NMRShielding final { const DoubleMatrix &getDiamagnetic() const { return this->dia_tensor; } const DoubleMatrix &getParamagnetic() const { return this->para_tensor; } - void print() const { - auto iso_ppm_d = getDiamagnetic().trace() / 3.0; - auto iso_ppm_p = getParamagnetic().trace() / 3.0; - auto iso_ppm_t = iso_ppm_d + iso_ppm_p; + void print(const std::string &id) const { + auto sigma = getTensor(); + Eigen::SelfAdjointEigenSolver es; + es.compute(sigma); - mrcpp::print::header(0, "NMR shielding"); - print_utils::scalar(0, "Nucleus K", getK(), getNucleus().getElement().getSymbol(), 0); + auto diag = es.eigenvalues(); + auto iso_ppm = diag.sum() / 3.0; + auto ani_ppm = diag(2) - (diag(0) + diag(1)) / 2.0; + + mrcpp::print::header(0, "NMR shielding (" + id + ")"); + print_utils::coord(0, "r_O", getOrigin()); + print_utils::coord(0, "r_K", getCoordK()); mrcpp::print::separator(0, '-'); print_utils::matrix(0, "Diamagnetic", getDiamagnetic()); - print_utils::scalar(0, "Isotropic average", iso_ppm_d, "(ppm)"); mrcpp::print::separator(0, '-'); - print_utils::matrix(0, "Paramagnetic", getParamagnetic(), -1); - print_utils::scalar(0, "Isotropic average", iso_ppm_p, "(ppm)"); + print_utils::matrix(0, "Paramagnetic", getParamagnetic()); mrcpp::print::separator(0, '-'); print_utils::matrix(0, "Total tensor", getTensor()); - print_utils::scalar(0, "Isotropic average", iso_ppm_t, "(ppm)"); + mrcpp::print::separator(0, '-'); + print_utils::vector(0, "Diagonalized tensor", diag); + print_utils::scalar(0, "Isotropic average", iso_ppm, "(ppm)"); + print_utils::scalar(0, "Anisotropy", ani_ppm, "(ppm)"); mrcpp::print::separator(0, '=', 2); } + nlohmann::json json() const { + auto sigma = getTensor(); + Eigen::SelfAdjointEigenSolver es; + es.compute(sigma); + + auto diag = es.eigenvalues(); + auto iso_ppm = diag.sum() / 3.0; + auto ani_ppm = diag(2) - (diag(0) + diag(1)) / 2.0; + return { + {"r_K", getCoordK()}, + {"r_O", getOrigin()}, + {"tensor_dia", print_utils::eigen_to_vector(getDiamagnetic(), 1.0e-12)}, + {"tensor_para", print_utils::eigen_to_vector(getParamagnetic(), 1.0e-12)}, + {"tensor", print_utils::eigen_to_vector(sigma, 1.0e-12)}, + {"diagonalized_tensor", print_utils::eigen_to_vector(diag, 1.0e-12)}, + {"isotropic_average", iso_ppm}, + {"anisotropy", ani_ppm} + }; + } + private: - const int K; - const Nucleus nuc; + mrcpp::Coord<3> r_K; + mrcpp::Coord<3> r_O; DoubleMatrix dia_tensor{math_utils::init_nan(3,3)}; DoubleMatrix para_tensor{math_utils::init_nan(3,3)}; }; diff --git a/src/properties/OrbitalEnergies.h b/src/properties/OrbitalEnergies.h new file mode 100644 index 000000000..a5f3395f9 --- /dev/null +++ b/src/properties/OrbitalEnergies.h @@ -0,0 +1,111 @@ +/* + * 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) 2020 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 . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#pragma once + +#include + +#include "mrchem.h" +#include "utils/print_utils.h" + +/** @class OrbitalEnergies + * + * @brief Simple POD container to hold the orbital eigenvalues + * + */ + +namespace mrchem { + +// clang-format off +class OrbitalEnergies final { +public: + IntVector &getSpin() { return this->spin; } + IntVector &getOccupation() { return this->occupation; } + DoubleVector &getEpsilon() { return this->epsilon; } + + const IntVector &getSpin() const { return this->spin; } + const IntVector &getOccupation() const { return this->occupation; } + const DoubleVector &getEpsilon() const { return this->epsilon; } + + 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) << "n"; + o_head << std::setw(w1) << "Occ"; + o_head << std::setw(w1) << "Spin"; + o_head << std::string(w3 - 1, ' ') << ':'; + o_head << std::setw(3 * w2) << "Epsilon"; + + mrcpp::print::header(0, "Orbital Energies (" + id + ")"); + println(0, o_head.str()); + mrcpp::print::separator(0, '-'); + + for (int i = 0; i < this->epsilon.size(); i++) { + auto sp = 'u'; + if (this->spin(i) == SPIN::Paired) sp = 'p'; + if (this->spin(i) == SPIN::Alpha) sp = 'a'; + if (this->spin(i) == SPIN::Beta) sp = 'b'; + std::stringstream o_txt; + o_txt << std::setw(w1 - 1) << i; + o_txt << std::setw(w1) << this->occupation(i); + o_txt << std::setw(w1) << sp; + print_utils::scalar(0, o_txt.str(), this->epsilon(i), "(au)", pprec); + } + auto sum_occupied = this->occupation.cast().dot(this->epsilon); + mrcpp::print::separator(0, '-'); + print_utils::scalar(0, "Sum occupied", sum_occupied, "(au)", pprec); + mrcpp::print::separator(0, '=', 2); + } + + nlohmann::json json() const { + DoubleVector occ = getOccupation().cast(); + const DoubleVector &eps = getEpsilon(); + std::vector spn; + for (auto i = 0; i < spin.size(); i++) { + if (this->spin(i) == SPIN::Paired) spn.push_back("p"); + if (this->spin(i) == SPIN::Alpha) spn.push_back("a"); + if (this->spin(i) == SPIN::Beta) spn.push_back("b"); + } + return { + {"spin", spn}, + {"occupation", print_utils::eigen_to_vector(occ, 1.0e-12)}, + {"energy", print_utils::eigen_to_vector(eps, 1.0e-12)}, + {"sum_occupied", occ.dot(eps)} + }; + } + +private: + IntVector spin; + IntVector occupation; + DoubleVector epsilon; +}; +// clang-format on + +} // namespace mrchem diff --git a/src/properties/Polarizability.h b/src/properties/Polarizability.h index 6c9c196c0..88e472128 100644 --- a/src/properties/Polarizability.h +++ b/src/properties/Polarizability.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "mrchem.h" #include "utils/math_utils.h" @@ -35,32 +37,47 @@ namespace mrchem { // clang-format off class Polarizability final { public: - explicit Polarizability(double w = 0.0) : frequency(w) {} + explicit Polarizability(double omega = 0.0, const mrcpp::Coord<3> &o = {}) : frequency(omega), r_O(o) {} + + void setFrequency(double omega) { this->frequency = omega; } + void setOrigin(const mrcpp::Coord<3> &o) { this->r_O = o; } double getFrequency() const { return this->frequency; } - double getWavelength() const { return 0.0; } + const mrcpp::Coord<3> &getOrigin() const { return this->r_O; } + DoubleMatrix &getTensor() { return this->tensor; } const DoubleMatrix &getTensor() const { return this->tensor; } - void print() const { + void print(const std::string &id) const { auto w_au = getFrequency(); auto w_cm = PHYSCONST::cm_m1 * w_au; auto dynamic = (w_au > mrcpp::MachineZero); auto l_nm = (dynamic) ? (1.0e7 / w_cm) : 0.0; auto iso_au = getTensor().trace() / 3.0; - mrcpp::print::header(0, "Polarizability"); + mrcpp::print::header(0, "Polarizability (" + id + ")"); if (dynamic) print_utils::scalar(0, "Wavelength", l_nm, "(nm)"); print_utils::scalar(0, "Frequency", w_au, "(au)"); print_utils::scalar(0, " ", w_cm, "(cm-1)"); + print_utils::coord(0, "r_O", getOrigin()); mrcpp::print::separator(0, '-'); print_utils::matrix(0, "Total tensor", getTensor()); print_utils::scalar(0, "Isotropic average", iso_au, "(au)"); mrcpp::print::separator(0, '=', 2); } + nlohmann::json json() const { + return { + {"r_O", getOrigin()}, + {"frequency", getFrequency()}, + {"tensor", print_utils::eigen_to_vector(getTensor(), 1.0e-12)}, + {"isotropic_average", getTensor().trace() / 3.0 } + }; + } + private: double frequency; + mrcpp::Coord<3> r_O; DoubleMatrix tensor{math_utils::init_nan(3,3)}; }; // clang-format on diff --git a/src/properties/QuadrupoleMoment.h b/src/properties/QuadrupoleMoment.h index ba7972865..b5677ce68 100644 --- a/src/properties/QuadrupoleMoment.h +++ b/src/properties/QuadrupoleMoment.h @@ -35,14 +35,21 @@ namespace mrchem { // clang-format off class QuadrupoleMoment final { public: + explicit QuadrupoleMoment(const mrcpp::Coord<3> &o = {}) : r_O(o) {} + + void setOrigin(const mrcpp::Coord<3> &o) { this->r_O = o; } + const mrcpp::Coord<3> &getOrigin() const { return this->r_O; } + DoubleMatrix getTensor() const { return getNuclear() + getElectronic(); } DoubleMatrix &getNuclear() { return this->nuc_tensor; } DoubleMatrix &getElectronic() { return this->el_tensor; } const DoubleMatrix &getNuclear() const { return this->nuc_tensor; } const DoubleMatrix &getElectronic() const { return this->el_tensor; } - void print() const { - mrcpp::print::header(0, "Quadrupole Moment"); + void print(const std::string &id) const { + mrcpp::print::header(0, "Quadrupole Moment (" + id + ")"); + print_utils::coord(0, "r_O", getOrigin()); + mrcpp::print::separator(0, '-'); print_utils::matrix(0, "Electronic tensor", getElectronic()); mrcpp::print::separator(0, '-'); print_utils::matrix(0, "Nuclear tensor", getNuclear()); @@ -51,7 +58,17 @@ class QuadrupoleMoment final { mrcpp::print::separator(0, '=', 2); } + nlohmann::json json() const { + return { + {"r_O", getOrigin()}, + {"tensor_nuc", print_utils::eigen_to_vector(getNuclear(), 1.0e-12)}, + {"tensor_el", print_utils::eigen_to_vector(getElectronic(), 1.0e-12)}, + {"tensor", print_utils::eigen_to_vector(getTensor(), 1.0e-12)} + }; + } + private: + mrcpp::Coord<3> r_O; DoubleMatrix nuc_tensor{math_utils::init_nan(3)}; DoubleMatrix el_tensor{math_utils::init_nan(3)}; }; diff --git a/src/properties/SCFEnergy.h b/src/properties/SCFEnergy.h index 1fbf81ecf..286053cfc 100644 --- a/src/properties/SCFEnergy.h +++ b/src/properties/SCFEnergy.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "mrchem.h" #include "utils/print_utils.h" @@ -39,10 +41,10 @@ namespace mrchem { // clang-format off class SCFEnergy final { public: - SCFEnergy() = default; - SCFEnergy(double kin, double nn, double en, - double ee, double x, double xc, - double next, double eext) : + explicit SCFEnergy(double kin = 0.0, double nn = 0.0, + double en = 0.0, double ee = 0.0, + double x = 0.0, double xc = 0.0, + double next = 0.0, double eext = 0.0) : E_kin(kin), E_nn(nn), E_en(en), E_ee(ee), E_x(x), E_xc(xc), E_next(next), E_eext(eext) { E_nuc = E_nn + E_next; @@ -54,25 +56,30 @@ class SCFEnergy final { double getElectronicEnergy() const { return this->E_el; } double getKineticEnergy() const { return this->E_kin; } + double getNuclearNuclearEnergy() const { return this->E_nn; } double getElectronNuclearEnergy() const { return this->E_en; } double getElectronElectronEnergy() const { return this->E_ee; } + double getElectronExternalEnergy() const { return this->E_eext; } + double getNuclearExternalEnergy() const { return this->E_next; } double getExchangeCorrelationEnergy() const { return this->E_xc; } double getExchangeEnergy() const { return this->E_x; } - void print() const { + void print(const std::string &id) const { auto E_au = E_nuc + E_el; auto E_eV = E_au * PHYSCONST::eV; auto E_kJ = E_au * PHYSCONST::kJ; auto E_kcal = E_au * PHYSCONST::kcal; auto pprec = 2 * mrcpp::Printer::getPrecision(); - mrcpp::print::header(0, "Molecular Energy"); + mrcpp::print::header(0, "Molecular Energy (" + id + ")"); print_utils::scalar(0, "Kinetic energy ", E_kin, "(au)", pprec, false); print_utils::scalar(0, "E-N energy ", E_en, "(au)", pprec, false); print_utils::scalar(0, "Coulomb energy ", E_ee, "(au)", pprec, false); print_utils::scalar(0, "Exchange energy ", E_x, "(au)", pprec, false); print_utils::scalar(0, "X-C energy ", E_xc, "(au)", pprec, false); print_utils::scalar(0, "Ext. field (el) ", E_eext, "(au)", pprec, false); + mrcpp::print::separator(0, '-'); + print_utils::scalar(0, "N-N energy ", E_nn, "(au)", pprec, false); print_utils::scalar(0, "Ext. field (nuc) ", E_next, "(au)", pprec, false); mrcpp::print::separator(0, '-'); print_utils::scalar(0, "Electronic energy", E_el, "(au)", pprec, false); @@ -85,6 +92,22 @@ class SCFEnergy final { mrcpp::print::separator(0, '=', 2); } + nlohmann::json json() const { + return { + {"E_kin", E_kin}, + {"E_nn", E_nn}, + {"E_en", E_en}, + {"E_ee", E_ee}, + {"E_next", E_next}, + {"E_eext", E_eext}, + {"E_x", E_x}, + {"E_xc", E_xc}, + {"E_el", E_el}, + {"E_nuc", E_nuc}, + {"E_tot", E_nuc + E_el} + }; + } + private: double E_nuc{0.0}; double E_el{0.0}; diff --git a/src/properties/properties_fwd.h b/src/properties/properties_fwd.h index 252cfc6ab..295453482 100644 --- a/src/properties/properties_fwd.h +++ b/src/properties/properties_fwd.h @@ -32,5 +32,6 @@ class DipoleMoment; class Magnetizability; class NMRShielding; class Polarizability; +class OrbitalEnergies; } // namespace mrchem diff --git a/src/qmfunctions/Orbital.cpp b/src/qmfunctions/Orbital.cpp index 630ffce1d..eb95938a2 100644 --- a/src/qmfunctions/Orbital.cpp +++ b/src/qmfunctions/Orbital.cpp @@ -44,7 +44,7 @@ Orbital::Orbital() /** @brief Constructor * * @param spin: electron spin (SPIN::Alpha/Beta/Paired) - * @param occ: occupancy + * @param occ: occupation * @param rank: MPI ownership (-1 means all MPI ranks) * * Initializes the QMFunction with NULL pointers for both real and imaginary part. @@ -88,7 +88,7 @@ Orbital &Orbital::operator=(const Orbital &orb) { /** @brief Parameter copy * - * Returns a new orbital with the same spin, occupancy and rank_id as *this orbital. + * Returns a new orbital with the same spin, occupation and rank_id as *this orbital. */ Orbital Orbital::paramCopy() const { return Orbital(this->spin(), this->occ(), this->rankID()); diff --git a/src/qmfunctions/orbital_utils.cpp b/src/qmfunctions/orbital_utils.cpp index 1095423e1..6c363d292 100644 --- a/src/qmfunctions/orbital_utils.cpp +++ b/src/qmfunctions/orbital_utils.cpp @@ -23,9 +23,9 @@ * */ -#include "MRCPP/Printer" -#include "MRCPP/Timer" -#include "MRCPP/utils/details.h" +#include +#include +#include #include "parallel.h" #include "utils/RRMaximizer.h" @@ -103,15 +103,15 @@ ComplexDouble orbital::node_norm_dot(Orbital bra, Orbital ket, bool exact) { return qmfunction::node_norm_dot(bra, ket, exact); } -/** @brief Compare spin and occupancy of two orbitals +/** @brief Compare spin and occupation of two orbitals * * Returns true if orbital parameters are the same. * */ bool orbital::compare(const Orbital &phi_a, const Orbital &phi_b) { bool comp = true; - if (compare_occ(phi_a, phi_b) < 0) { - MSG_WARN("Different occupancy"); + if (compare_occupation(phi_a, phi_b) < 0) { + MSG_WARN("Different occupation"); comp = false; } if (compare_spin(phi_a, phi_b) < 0) { @@ -121,12 +121,12 @@ bool orbital::compare(const Orbital &phi_a, const Orbital &phi_b) { return comp; } -/** @brief Compare occupancy of two orbitals +/** @brief Compare occupation of two orbitals * - * Returns the common occupancy if they match, -1 if they differ. + * Returns the common occupation if they match, -1 if they differ. * */ -int orbital::compare_occ(const Orbital &phi_a, const Orbital &phi_b) { +int orbital::compare_occupation(const Orbital &phi_a, const Orbital &phi_b) { int comp = -1; if (phi_a.occ() == phi_b.occ()) comp = phi_a.occ(); return comp; @@ -143,7 +143,7 @@ int orbital::compare_spin(const Orbital &phi_a, const Orbital &phi_b) { return comp; } -/** @brief Compare spin and occupancy of two orbital vector +/** @brief Compare spin and occupation of two orbital vector * * Returns true if orbital parameters are the same, orbital ordering * NOT taken into account. @@ -870,20 +870,20 @@ void orbital::set_spins(OrbitalVector &Phi, const IntVector &spins) { for (int i = 0; i < Phi.size(); i++) Phi[i].setSpin(spins(i)); } -/** @brief Returns a vector containing the orbital occupancies */ -IntVector orbital::get_occupancies(const OrbitalVector &Phi) { +/** @brief Returns a vector containing the orbital occupations */ +IntVector orbital::get_occupations(const OrbitalVector &Phi) { int nOrbs = Phi.size(); IntVector occ = IntVector::Zero(nOrbs); for (int i = 0; i < nOrbs; i++) occ(i) = Phi[i].occ(); return occ; } -/** @brief Assigns spin to each orbital +/** @brief Assigns occupation to each orbital * * Length of input vector must match the number of orbitals in the set. * */ -void orbital::set_occupancies(OrbitalVector &Phi, const IntVector &occ) { +void orbital::set_occupations(OrbitalVector &Phi, const IntVector &occ) { if (Phi.size() != occ.size()) MSG_ERROR("Size mismatch"); for (int i = 0; i < Phi.size(); i++) Phi[i].setOcc(occ(i)); } @@ -995,60 +995,36 @@ void orbital::print(const OrbitalVector &Phi) { mrcpp::print::separator(0, '=', 2); } -void orbital::print_eigenvalues(const OrbitalVector &Phi, const ComplexMatrix &F_mat) { - if (Phi.size() == 0) return; +DoubleVector orbital::calc_eigenvalues(const OrbitalVector &Phi, const ComplexMatrix &F_mat) { if (F_mat.cols() != Phi.size()) MSG_ABORT("Invalid Fock matrix"); + if (not orbital::orbital_vector_is_sane(Phi)) MSG_ABORT("Insane orbital vector"); - // First compute eigenvalues without rotating the orbitals DoubleVector epsilon = DoubleVector::Zero(Phi.size()); int np = orbital::size_paired(Phi); int na = orbital::size_alpha(Phi); int nb = orbital::size_beta(Phi); if (np > 0) { + Timer timer; Eigen::SelfAdjointEigenSolver es(np); es.compute(F_mat.block(0, 0, np, np)); epsilon.segment(0, np) = es.eigenvalues(); + mrcpp::print::time(1, "Diagonalize Fock matrix", timer); } if (na > 0) { + Timer timer; Eigen::SelfAdjointEigenSolver es(na); es.compute(F_mat.block(np, np, na, na)); epsilon.segment(np, na) = es.eigenvalues(); + mrcpp::print::time(1, "Diagonalize Fock matrix (alpha)", timer); } if (nb > 0) { + Timer timer; Eigen::SelfAdjointEigenSolver es(nb); es.compute(F_mat.block(np + na, np + na, nb, nb)); epsilon.segment(np + na, nb) = es.eigenvalues(); + mrcpp::print::time(1, "Diagonalize Fock matrix (beta)", timer); } - - auto pprec = Printer::getPrecision(); - auto w0 = 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) << "n"; - o_head << std::setw(w1) << "Occ"; - o_head << std::setw(w1) << "Spin"; - o_head << std::string(w3 - 1, ' ') << ':'; - o_head << std::setw(3 * w2) << "Epsilon"; - - mrcpp::print::header(0, "Orbital Energies"); - println(0, o_head.str()); - mrcpp::print::separator(0, '-'); - - auto sum_eps = 0.0; - for (int i = 0; i < epsilon.size(); i++) { - std::stringstream o_txt; - o_txt << std::setw(w1 - 1) << i; - o_txt << std::setw(w1) << Phi[i].occ(); - o_txt << std::setw(w1) << Phi[i].printSpin(); - print_utils::scalar(0, o_txt.str(), epsilon(i), "(au)", 2 * pprec); - sum_eps += Phi[i].occ() * epsilon(i); - } - mrcpp::print::separator(0, '-'); - print_utils::scalar(0, "Sum occupied", sum_eps, "(au)", 2 * pprec); - mrcpp::print::separator(0, '=', 2); + return epsilon; } /** @brief Prints statistics about the size of orbitals in an OrbitalVector diff --git a/src/qmfunctions/orbital_utils.h b/src/qmfunctions/orbital_utils.h index 95b69e420..87dbea0c4 100644 --- a/src/qmfunctions/orbital_utils.h +++ b/src/qmfunctions/orbital_utils.h @@ -33,8 +33,8 @@ namespace orbital { bool compare(const OrbitalVector &Phi_a, const OrbitalVector &Phi_b); bool compare(const Orbital &phi_a, const Orbital &phi_b); -int compare_occ(const Orbital &phi_a, const Orbital &phi_b); int compare_spin(const Orbital &phi_a, const Orbital &phi_b); +int compare_occupation(const Orbital &phi_a, const Orbital &phi_b); ComplexDouble dot(Orbital bra, Orbital ket); ComplexVector dot(OrbitalVector &Bra, OrbitalVector &Ket); @@ -83,16 +83,16 @@ int get_size_nodes(const OrbitalVector &Phi); bool orbital_vector_is_sane(const OrbitalVector &Phi); void set_spins(OrbitalVector &Phi, const IntVector &spins); -void set_occupancies(OrbitalVector &Phi, const IntVector &occ); +void set_occupations(OrbitalVector &Phi, const IntVector &occ); IntVector get_spins(const OrbitalVector &Phi); -IntVector get_occupancies(const OrbitalVector &Phi); +IntVector get_occupations(const OrbitalVector &Phi); DoubleVector get_norms(const OrbitalVector &Phi); DoubleVector get_squared_norms(const OrbitalVector &Phi); +DoubleVector calc_eigenvalues(const OrbitalVector &Phi, const ComplexMatrix &F_mat); ComplexVector get_integrals(const OrbitalVector &Phi); void print(const OrbitalVector &Phi); -void print_eigenvalues(const OrbitalVector &Phi, const ComplexMatrix &F_mat); int print_size_nodes(const OrbitalVector &Phi, const std::string &txt = "", bool all = true, int plevel = 0); } // namespace orbital diff --git a/src/qmoperators/RankZeroTensorOperator.cpp b/src/qmoperators/RankZeroTensorOperator.cpp index 4f3b0bbb7..685bf9efa 100644 --- a/src/qmoperators/RankZeroTensorOperator.cpp +++ b/src/qmoperators/RankZeroTensorOperator.cpp @@ -329,7 +329,7 @@ ComplexDouble RankZeroTensorOperator::trace(OrbitalVector &Phi) { Timer t1; RankZeroTensorOperator &O = *this; OrbitalVector OPhi = O(Phi); - ComplexVector eta = orbital::get_occupancies(Phi).cast(); + ComplexVector eta = orbital::get_occupations(Phi).cast(); ComplexVector phi_vec = orbital::dot(Phi, OPhi); std::stringstream o_name; @@ -372,7 +372,7 @@ ComplexDouble RankZeroTensorOperator::trace(OrbitalVector &Phi, OrbitalVector &X o_name << "Trace " << O.name() << "(rho_1)"; mrcpp::print::tree(2, o_name.str(), std::max(x_nodes, y_nodes), std::max(x_size, y_size), t1.elapsed()); - ComplexVector eta = orbital::get_occupancies(Phi).cast(); + ComplexVector eta = orbital::get_occupations(Phi).cast(); return eta.dot(x_vec + y_vec); } diff --git a/src/qmoperators/two_electron/FockOperator.cpp b/src/qmoperators/two_electron/FockOperator.cpp index ad173f606..8b4c13824 100644 --- a/src/qmoperators/two_electron/FockOperator.cpp +++ b/src/qmoperators/two_electron/FockOperator.cpp @@ -128,7 +128,7 @@ void FockOperator::rotate(const ComplexMatrix &U) { SCFEnergy FockOperator::trace(OrbitalVector &Phi, const Nuclei &nucs) { Timer t_tot; auto plevel = Printer::getPrintLevel(); - mrcpp::print::header(2, "Calculating molecular energy"); + mrcpp::print::header(2, "Computing molecular energy"); double E_kin = 0.0; // Kinetic energy double E_nn = 0.0; // Nuclear repulsion @@ -151,7 +151,7 @@ SCFEnergy FockOperator::trace(OrbitalVector &Phi, const Nuclei &nucs) { if (this->xc != nullptr) E_xc = this->xc->getEnergy(); if (this->ext != nullptr) E_eext = this->ext->trace(Phi).real(); mrcpp::print::footer(2, t_tot, 2); - if (plevel == 1) mrcpp::print::time(1, "Calculating molecular energy", t_tot); + if (plevel == 1) mrcpp::print::time(1, "Computing molecular energy", t_tot); return SCFEnergy{E_kin, E_nn, E_en, E_ee, E_x, E_xc, E_next, E_eext}; } @@ -159,7 +159,7 @@ SCFEnergy FockOperator::trace(OrbitalVector &Phi, const Nuclei &nucs) { ComplexMatrix FockOperator::operator()(OrbitalVector &bra, OrbitalVector &ket) { Timer t_tot; auto plevel = Printer::getPrintLevel(); - mrcpp::print::header(2, "Calculating Fock matrix"); + mrcpp::print::header(2, "Computing Fock matrix"); auto t = this->getKineticOperator(); auto v = this->potential(); @@ -171,7 +171,7 @@ ComplexMatrix FockOperator::operator()(OrbitalVector &bra, OrbitalVector &ket) { if (v.size() > 0) V += v(bra, ket); mrcpp::print::footer(2, t_tot, 2); - if (plevel == 1) mrcpp::print::time(1, "Calculating Fock matrix", t_tot); + if (plevel == 1) mrcpp::print::time(1, "Computing Fock matrix", t_tot); return T + V; } diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp index 44b4f1a28..3bb1a8ce9 100644 --- a/src/scf_solver/GroundStateSolver.cpp +++ b/src/scf_solver/GroundStateSolver.cpp @@ -37,6 +37,7 @@ using mrcpp::Printer; using mrcpp::Timer; +using nlohmann::json; namespace mrchem { @@ -187,9 +188,8 @@ void GroundStateSolver::reset() { /** @brief Run orbital optimization * + * @param mol: Molecule to optimize * @param F: FockOperator defining the SCF equations - * @param Phi_n: Orbitals to optimize - * @param F_mat: Fock matrix * * Optimize orbitals until convergence thresholds are met. This algorithm computes * the Fock matrix explicitly using the kinetic energy operator, and uses a KAIN @@ -210,10 +210,11 @@ void GroundStateSolver::reset() { * 10) Setup Fock operator * 11) Compute Fock matrix * - * Post SCF: diagonalize/localize orbitals */ -bool GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { +json GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { printParameters("Optimize ground state orbitals"); + Timer t_tot; + json json_out; KAIN kain(this->history); SCFEnergy &E_n = mol.getSCFEnergy(); @@ -237,7 +238,9 @@ bool GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { int nIter = 0; bool converged = false; + json_out["cycles"] = {}; while (nIter++ < this->maxIter or this->maxIter < 0) { + json json_cycle; std::stringstream o_header; o_header << "SCF cycle " << nIter; mrcpp::print::header(1, o_header.str(), 0, '#'); @@ -283,6 +286,7 @@ bool GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { errors = orbital::get_norms(dPhi_n); err_o = errors.maxCoeff(); err_t = errors.norm(); + json_cycle["mo_residual"] = err_t; // Update orbitals Phi_n = orbital::add(1.0, Phi_n, 1.0, dPhi_n); @@ -302,6 +306,10 @@ bool GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { auto err_p = calcPropertyError(); converged = checkConvergence(err_o, err_p); + json_cycle["energy_terms"] = E_n.json(); + json_cycle["energy_total"] = E_n.getTotalEnergy(); + json_cycle["energy_update"] = err_p; + // Rotate orbitals if (needLocalization(nIter, converged)) { ComplexMatrix U_mat = orbital::localize(orb_prec, Phi_n, F_mat); @@ -324,9 +332,12 @@ bool GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { mrcpp::print::separator(2, '=', 2); printProperty(); printMemory(); + t_scf.stop(); + json_cycle["wall_time"] = t_scf.elapsed(); mrcpp::print::footer(1, t_scf, 2, '#'); mrcpp::print::separator(2, ' ', 2); + json_out["cycles"].push_back(json_cycle); if (converged) break; } @@ -335,7 +346,17 @@ bool GroundStateSolver::optimize(Molecule &mol, FockOperator &F) { printConvergence(converged, "Total energy"); reset(); - return converged; + Timer t_eps; + mrcpp::print::header(1, "Computing orbital energies"); + OrbitalEnergies &eps = mol.getOrbitalEnergies(); + eps.getOccupation() = orbital::get_occupations(Phi_n); + eps.getEpsilon() = orbital::calc_eigenvalues(Phi_n, F_mat); + eps.getSpin() = orbital::get_spins(Phi_n); + mrcpp::print::footer(1, t_eps, 2); + + json_out["wall_time"] = t_tot.elapsed(); + json_out["converged"] = converged; + return json_out; } /** @brief Test if orbitals needs localization diff --git a/src/scf_solver/GroundStateSolver.h b/src/scf_solver/GroundStateSolver.h index 8d41230ed..199a76c3f 100644 --- a/src/scf_solver/GroundStateSolver.h +++ b/src/scf_solver/GroundStateSolver.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "SCFSolver.h" #include "properties/SCFEnergy.h" @@ -57,7 +59,7 @@ class GroundStateSolver : public SCFSolver { void setLocalize(bool loc) { this->localize = loc; } void setCheckpointFile(const std::string &file) { this->chkFile = file; } - bool optimize(Molecule &mol, FockOperator &F); + nlohmann::json optimize(Molecule &mol, FockOperator &F); protected: int rotation{0}; ///< Number of iterations between localization/diagonalization diff --git a/src/scf_solver/LinearResponseSolver.cpp b/src/scf_solver/LinearResponseSolver.cpp index f22693ff9..732da796c 100644 --- a/src/scf_solver/LinearResponseSolver.cpp +++ b/src/scf_solver/LinearResponseSolver.cpp @@ -38,6 +38,7 @@ using mrcpp::Printer; using mrcpp::Timer; +using nlohmann::json; namespace mrchem { @@ -59,8 +60,10 @@ namespace mrchem { * 5) Check for convergence * */ -bool LinearResponseSolver::optimize(double omega, Molecule &mol, FockOperator &F_0, FockOperator &F_1) { +json LinearResponseSolver::optimize(double omega, Molecule &mol, FockOperator &F_0, FockOperator &F_1) { printParameters(omega, F_1.perturbation().name()); + Timer t_tot; + json json_out; // Setup KAIN accelerators KAIN kain_x(this->history); @@ -98,7 +101,9 @@ bool LinearResponseSolver::optimize(double omega, Molecule &mol, FockOperator &F int nIter = 0; bool converged = false; + json_out["cycles"] = {}; while (nIter++ < this->maxIter or this->maxIter < 0) { + json json_cycle; std::stringstream o_header; o_header << "SCF cycle " << nIter; mrcpp::print::header(1, o_header.str(), 0, '#'); @@ -217,7 +222,6 @@ bool LinearResponseSolver::optimize(double omega, Molecule &mol, FockOperator &F mrcpp::print::header(2, "Computing symmetric property"); t_lap.start(); double prop = F_1.perturbation().trace(Phi_0, X_n, Y_n).real(); - this->property.push_back(prop); mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Computing symmetric property", t_lap); @@ -225,11 +229,18 @@ bool LinearResponseSolver::optimize(double omega, Molecule &mol, FockOperator &F F_1.clear(); // Compute errors - auto err_p = std::abs(getUpdate(this->property, nIter, false)); err_o = std::max(errors_x.maxCoeff(), errors_y.maxCoeff()); err_t = std::sqrt(errors_x.dot(errors_x) + errors_y.dot(errors_y)); - converged = checkConvergence(err_o, err_p); + json_cycle["mo_residual"] = err_t; + + // Collect convergence data this->error.push_back(err_t); + this->property.push_back(prop); + double err_p = getUpdate(this->property, nIter + 1, true); + converged = checkConvergence(err_o, err_p); + + json_cycle["symmetric_property"] = prop; + json_cycle["property_update"] = err_p; // Finalize SCF cycle if (plevel < 1) printConvergenceRow(nIter); @@ -240,16 +251,21 @@ bool LinearResponseSolver::optimize(double omega, Molecule &mol, FockOperator &F mrcpp::print::separator(2, '=', 2); printProperty(); printMemory(); + t_scf.stop(); + json_cycle["wall_time"] = t_scf.elapsed(); mrcpp::print::footer(1, t_scf, 2, '#'); mrcpp::print::separator(2, ' ', 2); + json_out["cycles"].push_back(json_cycle); if (converged) break; } printConvergence(converged, "Symmetric property"); reset(); - return converged; + json_out["wall_time"] = t_tot.elapsed(); + json_out["converged"] = converged; + return json_out; } /** @brief Pretty printing of the computed property with update */ diff --git a/src/scf_solver/LinearResponseSolver.h b/src/scf_solver/LinearResponseSolver.h index adb924b3a..47eb2b30e 100644 --- a/src/scf_solver/LinearResponseSolver.h +++ b/src/scf_solver/LinearResponseSolver.h @@ -25,6 +25,8 @@ #pragma once +#include + #include "SCFSolver.h" /** @class LinearResponseSolver @@ -42,7 +44,7 @@ class LinearResponseSolver final : public SCFSolver { : dynamic(dyn) {} ~LinearResponseSolver() override = default; - bool optimize(double omega, Molecule &mol, FockOperator &F_0, FockOperator &F_1); + nlohmann::json optimize(double omega, Molecule &mol, FockOperator &F_0, FockOperator &F_1); void setOrthPrec(double prec) { this->orth_prec = prec; } void setCheckpointFile(const std::string &file_x, const std::string &file_y) { this->chkFileX = file_x; diff --git a/src/scf_solver/SCFSolver.cpp b/src/scf_solver/SCFSolver.cpp index 54f337070..a73fc6e77 100644 --- a/src/scf_solver/SCFSolver.cpp +++ b/src/scf_solver/SCFSolver.cpp @@ -114,8 +114,8 @@ double SCFSolver::getHelmholtzPrec() { bool SCFSolver::checkConvergence(double err_o, double err_p) const { bool conv_o = false; bool conv_p = false; - if (err_o < this->orbThrs or this->orbThrs < 0.0) conv_o = true; - if (err_p < this->propThrs or this->propThrs < 0.0) conv_p = true; + if (std::abs(err_o) < this->orbThrs or this->orbThrs < 0.0) conv_o = true; + if (std::abs(err_p) < this->propThrs or this->propThrs < 0.0) conv_p = true; return (conv_o and conv_p); } diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 950d7029a..45d6e74df 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -38,6 +38,28 @@ using mrcpp::Timer; namespace mrchem { +std::vector print_utils::eigen_to_vector(const DoubleVector &inp, double thrs) { + std::vector out; + for (int i = 0; i < inp.size(); i++) { + auto out_i = inp(i); + if (std::abs(out_i) < thrs) out_i = 0.0; + out.push_back(out_i); + } + return out; +} + +std::vector print_utils::eigen_to_vector(const DoubleMatrix &inp, double thrs) { + std::vector out; + for (int i = 0; i < inp.rows(); i++) { + for (int j = 0; j < inp.cols(); j++) { + auto out_ij = inp(i, j); + if (std::abs(out_ij) < thrs) out_ij = 0.0; + out.push_back(out_ij); + } + } + return out; +} + std::string print_utils::dbl_to_str(double d, int p, bool sci) { std::stringstream o_dbl; if (sci) { diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h index d5e132e83..6ec06160b 100644 --- a/src/utils/print_utils.h +++ b/src/utils/print_utils.h @@ -39,6 +39,8 @@ void vector(int level, const std::string &txt, const DoubleVector &val, int p = void matrix(int level, const std::string &txt, const DoubleMatrix &val, int p = -1, bool s = false); void qmfunction(int level, const std::string &txt, const QMFunction &func, mrcpp::Timer &t); std::string dbl_to_str(double d, int p, bool sci); +std::vector eigen_to_vector(const DoubleVector &inp, double thrs); +std::vector eigen_to_vector(const DoubleMatrix &inpm, double thrs); } // namespace print_utils } // namespace mrchem diff --git a/tests/h2_pol_lda/reference/mrchem.stdout b/tests/h2_pol_lda/reference/mrchem.stdout index 6e1559cac..df9a00669 100644 --- a/tests/h2_pol_lda/reference/mrchem.stdout +++ b/tests/h2_pol_lda/reference/mrchem.stdout @@ -167,6 +167,18 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s ====================================================================== +====================================================================== + Orbital Energies +---------------------------------------------------------------------- + n Occ Spin : Epsilon +---------------------------------------------------------------------- + 0 1 a : (au) -0.345283880497 + 1 1 b : (au) -0.345283880497 +---------------------------------------------------------------------- + Sum occupied : (au) -0.690567760993 +====================================================================== + + ********************************************************************** *** *** diff --git a/tests/h2_scf_hf/reference/mrchem.stdout b/tests/h2_scf_hf/reference/mrchem.stdout index f76debc61..9f1ef627d 100644 --- a/tests/h2_scf_hf/reference/mrchem.stdout +++ b/tests/h2_scf_hf/reference/mrchem.stdout @@ -167,6 +167,17 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s ====================================================================== +====================================================================== + Orbital Energies +---------------------------------------------------------------------- + n Occ Spin : Epsilon +---------------------------------------------------------------------- + 0 2 p : (au) -0.611656312944 +---------------------------------------------------------------------- + Sum occupied : (au) -1.223312625889 +====================================================================== + + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Calculation : Optimize ground state orbitals Method : Hartree-Fock diff --git a/tests/h2o_energy_blyp/reference/mrchem.stdout b/tests/h2o_energy_blyp/reference/mrchem.stdout index c9abeb6e3..0eddf7f78 100644 --- a/tests/h2o_energy_blyp/reference/mrchem.stdout +++ b/tests/h2o_energy_blyp/reference/mrchem.stdout @@ -161,6 +161,20 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s ====================================================================== +====================================================================== + Orbital Energies +---------------------------------------------------------------------- + n Occ Spin : Epsilon +---------------------------------------------------------------------- + 0 2 p : (au) -18.782961402002 + 1 2 p : (au) -0.892698452791 + 2 2 p : (au) -0.447021541210 + 3 2 p : (au) -0.303255086859 + 4 2 p : (au) -0.222581060353 +---------------------------------------------------------------------- + Sum occupied : (au) -41.297035086430 +====================================================================== + ********************************************************************** *** *** diff --git a/tests/h_el_field/reference/mrchem.stdout b/tests/h_el_field/reference/mrchem.stdout index f8b07349c..33c11eae0 100644 --- a/tests/h_el_field/reference/mrchem.stdout +++ b/tests/h_el_field/reference/mrchem.stdout @@ -165,6 +165,17 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s ====================================================================== +====================================================================== + Orbital Energies +---------------------------------------------------------------------- + n Occ Spin : Epsilon +---------------------------------------------------------------------- + 0 1 a : (au) -0.491935627808 +---------------------------------------------------------------------- + Sum occupied : (au) -0.491935627808 +====================================================================== + + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Calculation : Optimize ground state orbitals Method : Hartree-Fock diff --git a/tests/li_pol_lda/reference/mrchem.stdout b/tests/li_pol_lda/reference/mrchem.stdout index 8b0634abb..a6b0e37e5 100644 --- a/tests/li_pol_lda/reference/mrchem.stdout +++ b/tests/li_pol_lda/reference/mrchem.stdout @@ -160,6 +160,19 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s ====================================================================== +====================================================================== + Orbital Energies +---------------------------------------------------------------------- + n Occ Spin : Epsilon +---------------------------------------------------------------------- + 0 1 a : (au) -1.814574765685 + 1 1 a : (au) -0.114315988999 + 2 1 b : (au) -1.805755227476 +---------------------------------------------------------------------- + Sum occupied : (au) -3.734645982160 +====================================================================== + + ********************************************************************** *** *** diff --git a/tests/li_scf_pbe0/reference/mrchem.stdout b/tests/li_scf_pbe0/reference/mrchem.stdout index b5307a8c9..4cc7f4328 100644 --- a/tests/li_scf_pbe0/reference/mrchem.stdout +++ b/tests/li_scf_pbe0/reference/mrchem.stdout @@ -155,6 +155,20 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s ====================================================================== +====================================================================== + Orbital Energies +---------------------------------------------------------------------- + n Occ Spin : Epsilon +---------------------------------------------------------------------- + 0 1 a : (au) -2.026248624388 + 1 1 a : (au) -0.139732410810 + 2 1 b : (au) -2.015998440289 +---------------------------------------------------------------------- + Sum occupied : (au) -4.181979475487 +====================================================================== + + + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Calculation : Optimize molecular orbitals Method : DFT