From 18bc94935cc9b6f36fe992d7ccfd9f34c3935131 Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 18 Jun 2020 16:48:34 +0200 Subject: [PATCH] JSON output (#290) * Add Eigen -> std::vector conversion * Allow negative print_level for complete silence * Add json property output * Add OrbitalEnergies property * Revert mrchem_parsed.json intermediate Now mrchem.json contains both the parsed input and the json output. It's still possible to give user input as JSON, but file name should still be mrchem.inp for this. * Add missing license headers * Add json SCF output * Compute dipole moment by default * Add json RSP output * Remove GTO rsp initial guess * Add schema for program input * Make Response.localize default to SCF.localize * Only MPI grand master dumps JSON * Improve JSON output * Add schema for program output * Occupancy -> occupation * Polarizability and NMR json list -> dict * Flush very small property values to zero * Tweak NMR output * Update to nlohmann::json-3.6.1 * Make use of json.contains() * Cleanup JSON output * Fix schemas * Remove orbital_energies input keyword * Update test refs * Require full basis name (legendre/interpolating) * Make default constructors for properties * Cleanup Molecule's handling of properties * Geometry derivative -> geometric derivative * Move eigen_to_vector to print_utils * Define PropertyMap * Move calculation of dia properties to RSP part * Fix printing --- doc/schema_input.json | 261 ++++++++++ doc/schema_output.json | 144 +++++ doc/users/program_json.rst | 16 + doc/users/user_out.rst | 4 + external/upstream/fetch_nlohmann_json.cmake | 4 +- pilot/mrchem.cpp | 2 +- pilot/mrchem.in | 5 +- src/chemistry/Molecule.cpp | 105 ++-- src/chemistry/Molecule.h | 45 +- src/driver.cpp | 490 ++++++++++-------- src/driver.h | 11 +- src/input/mrchem.in | 122 +++-- src/input/template.yml | 11 +- src/mrchem.cpp | 52 +- src/mrchem.h | 26 +- src/mrenum.h | 25 + src/mrenv.cpp | 98 ++-- src/mrenv.h | 32 +- src/parallel.cpp | 25 + src/parallel.h | 25 + src/properties/DipoleMoment.h | 24 +- src/properties/Magnetizability.h | 37 +- src/properties/NMRShielding.h | 61 ++- src/properties/OrbitalEnergies.h | 111 ++++ src/properties/Polarizability.h | 25 +- src/properties/QuadrupoleMoment.h | 21 +- src/properties/SCFEnergy.h | 35 +- src/properties/properties_fwd.h | 1 + src/qmfunctions/Orbital.cpp | 4 +- src/qmfunctions/orbital_utils.cpp | 70 +-- src/qmfunctions/orbital_utils.h | 8 +- src/qmoperators/RankZeroTensorOperator.cpp | 4 +- src/qmoperators/two_electron/FockOperator.cpp | 8 +- src/scf_solver/GroundStateSolver.cpp | 31 +- src/scf_solver/GroundStateSolver.h | 4 +- src/scf_solver/LinearResponseSolver.cpp | 26 +- src/scf_solver/LinearResponseSolver.h | 4 +- src/scf_solver/SCFSolver.cpp | 4 +- src/utils/print_utils.cpp | 22 + src/utils/print_utils.h | 2 + tests/h2_pol_lda/reference/mrchem.stdout | 12 + tests/h2_scf_hf/reference/mrchem.stdout | 11 + tests/h2o_energy_blyp/reference/mrchem.stdout | 14 + tests/h_el_field/reference/mrchem.stdout | 11 + tests/li_pol_lda/reference/mrchem.stdout | 13 + tests/li_scf_pbe0/reference/mrchem.stdout | 14 + 46 files changed, 1547 insertions(+), 533 deletions(-) create mode 100644 doc/schema_input.json create mode 100644 doc/schema_output.json create mode 100644 doc/users/program_json.rst create mode 100644 doc/users/user_out.rst create mode 100644 src/properties/OrbitalEnergies.h 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