search

MRChemSoft / mrchem

MultiResolution Chemistry
GNU Lesser General Public License v3.0
27 stars 21 forks source link

Input parsing cleaning and PhysConsts from QCelemental #404

Closed Brakjen closed 2 years ago

Brakjen commented 2 years ago

Summary of changes

On the Python side

    - electron g factor        : CODATA: -2.00231930436256
                                 MRChem:  2.0023193043618

    - fine structure constant  : CODATA: 7.2973525693e-3
                                 MRChem: 7.2973525664

    - au -> wavenumbers        : CODATA: 219474.6313632
                                 MRChem: 219471.5125976648

    - au -> debye              : CODATA: 2.5417464739297717
                               : MRChem: 2.54174623105

On the c++ side

ToDo

stigrj commented 2 years ago

I believe this file python/mrchem/input_parser/plumbing/init.py should not be committed. It was the one causing problems before (reason still unknown). Are you now able to run with file in source?

Brakjen commented 2 years ago

Ah, that one slipped through. Very annoying. Why has this file started causing problems?

Another thing. The ANG_2_BOHR conversion factor is present both in helpers.py and CUBEparser.py. The cube parser cannot import from helpers due to circular import errors. A simple fix solution could be to collect the constants we use during input parsing in its own file, and import from there. Perhaps ideally python and c++ should access the same set of constants (this is perhaps what you mentioned in #171)

Gabrielgerez commented 2 years ago

That li_solv has always had weird behavior, I'll check it out

stigrj commented 2 years ago

@robertodr do you know why this file suddenly started to appear?

I think one nice way to handle the constants is to let Python define them and feed them into the C++ program through the JSON input. Then we also have complete control and can tweak the parameters without recompiling, by intercepting the input parser and edit the JSON file before execution.

robertodr commented 2 years ago

@Andersmb Which version of pyparsing is installed locally? parselglossy cannot yet deal with v3 of that library, I think that's what you're seeing. I opened a PR sometime ago to fix that https://github.com/dev-cafe/parselglossy/pull/109 Didn't get around to finishing it.

robertodr commented 2 years ago

The quick solution is to ignore any change in that plumbing when committing.

Brakjen commented 2 years ago

The quick solution is to ignore any change in that plumbing when committing.

I have version 3.0.7 installed...

robertodr commented 2 years ago

The quick solution is to ignore any change in that plumbing when committing.

I have version 3.0.7 installed...

Mystery understood then.

Brakjen commented 2 years ago

@Gabrielgerez I discovered why the li_solv was failing. The Molecule.cavity_coords section had a duplicated sphere

"cavity_coords": [
        {
          "center": [
            0.0,
            0.0,
            0.0
          ],
          "radius": 4.0
        },
        {
          "center": [
            0.0,
            0.0,
            0.0
          ],
          "radius": 4.0
        }
      ]

Just deleting the extra sphere lead to a passing test.

codecov[bot] commented 2 years ago

Codecov Report

Merging #404 (d135280) into master (c1959c6) will decrease coverage by 0.09%. The diff coverage is 52.11%.

@@            Coverage Diff             @@
##           master     #404      +/-   ##
==========================================
- Coverage   68.14%   68.04%   -0.10%     
==========================================
  Files         180      182       +2     
  Lines       15170    15208      +38     
==========================================
+ Hits        10337    10349      +12     
- Misses       4833     4859      +26
Impacted Files Coverage Δ
src/properties/Magnetizability.h 0.00% <0.00%> (ø)
src/qmoperators/one_electron/H_BM_dia.h 0.00% <0.00%> (ø)
src/qmoperators/one_electron/H_MB_dia.h 0.00% <0.00%> (ø)
src/qmoperators/one_electron/H_M_pso.h 0.00% <0.00%> (ø)
src/utils/print_utils.cpp 76.02% <0.00%> (-12.07%) :arrow_down:
src/chemistry/PhysicalConstants.cpp 50.00% <50.00%> (ø)
src/scf_solver/HelmholtzVector.cpp 92.98% <50.00%> (ø)
src/chemistry/PhysicalConstants.h 71.42% <71.42%> (ø)
src/analyticfunctions/HarmonicOscillatorFunction.h 100.00% <100.00%> (ø)
src/analyticfunctions/HydrogenFunction.cpp 61.34% <100.00%> (ø)
... and 11 more

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update c1959c6...d135280. Read the comment docs.

Brakjen commented 2 years ago

@stigrj Many unit tests will now fail because they are not being fed the constants generated by the input parser. How best solve this?

stigrj commented 2 years ago

Hmm yes, we probably need some defaults on the C++ side. Did you get around the linking issue?

robertodr commented 2 years ago

As for the singleton, I would use the pattern given in Alexandrescu's "Modern C++ Design" which uses static local references and avoids spilling around pointers and checking for NULLs.

PhysicalConstants.h

namespace mrchem {
class PhysicalConstants final {
  public:
     static PhysicalConstants& getInstance(const nlhomann::json & constants = /* */);
     static auto get(const std::string& key) const {
         return constants_[key];
     }
     PhysicalConstants() = default;
   ~PhysicalConstants() = default;
   PhysicalConstants(const PhysicalConstants&) = delete; 
   PhysicalConstants& operator=(const PhysicalConstants&) = delete;
   PhysicalConstants(PhysicalConstants&&) = delete; 
   PhysicalConstants& operator=(PhysicalConstants&&) = delete;
  private:
     PhysicalConstants(const nlhomann::json & constants) : constants_{constants} {}
     nlohmann::json constants_;
};
}

PhysicalConstants.cpp

namespace mrchem {
PhysicalConstants& PhysicalConstants::getInstance(const nlhomann::json & constants)
{
    static PhysicalConstants obj(constants);
    return obj;
}
}

Edit

Anders: I managed to edit your reply instead of quoting it 🤦 Some text may be missing....

Brakjen commented 2 years ago

@robertodr I will update singleton design.

The tests should pass now, but I am not sure if it is the best design. Now some constants are explicitly set in the singleton constructor. I added a method to MRChemConstants.py for printing a predefined subset of constants, to help avoid errors (you just need to copy/paste the output into the c++ file).

@stigrj Yes, the errors arose when I stored the instance in a variable. Accessing the constants without making a variable works:

double pi = mrchem::PhysicalConstants::getInstance()->get("pi");

Edit

The tests should pass when I manage to add qcelemental as a python dependency. Should I just update the Pipfile? Like so:

[[source]]
url = "https://pypi.python.org/simple"
verify_ssl = true
name = "pypi"

[packages]
qcelemental = ">=0.24.0"

[dev-packages]
Pygments = "*"
recommonmark = "*"
Sphinx = ">=2.0"
sphinx_rtd_theme = "*"
sphinxcontrib-bibtex = "*"
breathe = "*"
pyyaml = "*"
yapf = "*"
parselglossy = ">=0.7"
Brakjen commented 2 years ago

@robertodr I am not able to get this to work. After sorting out a few typos, I end up with Undefined symbols for architecture x86_64: errors during compilation, ala:

Undefined symbols for architecture x86_64:
  "mrchem::PhysicalConstants::constants_", referenced from:
      mrchem::driver::build_fock_operator(nlohmann::basic_json<std::__1::map, std::__1::vector, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, bool, long long, unsigned long long, double, std::__1::allocator, nlohmann::adl_serializer> const&, mrchem::Molecule&, mrchem::FockBuilder&, int) in driver.cpp.o
      mrchem::H_M_pso::H_M_pso(mrchem::MomentumOperator, mrchem::NuclearGradientOperator) in driver.cpp.o
      mrchem::H_MB_dia::H_MB_dia(mrchem::PositionOperator, mrchem::NuclearGradientOperator) in driver.cpp.o
      mrchem::H_BM_dia::H_BM_dia(mrchem::PositionOperator, mrchem::NuclearGradientOperator) in driver.cpp.o
      mrchem::AngularFunction::calcConstant() const in HydrogenFunction.cpp.o
      mrchem::chemistry::compute_nuclear_density(double, mrchem::Nuclei const&, double) in chemistry_utils.cpp.o
      mrchem::SCFEnergy::print(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&) const in Molecule.cpp.o
      ...
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

I currently have in PhysicalConstants.h

#pragma once
#include <nlohmann/json.hpp>

namespace mrchem {
class PhysicalConstants final {
  public:
     static PhysicalConstants& getInstance(const nlohmann::json &constants);
     static double get(const std::string& key) {
         return constants_[key];
     }
     PhysicalConstants() = default;
   ~PhysicalConstants() = default;
   PhysicalConstants(const PhysicalConstants&) = delete; 
   PhysicalConstants& operator=(const PhysicalConstants&) = delete;
   PhysicalConstants(PhysicalConstants&&) = delete; 
   PhysicalConstants& operator=(PhysicalConstants&&) = delete;

  private:
     PhysicalConstants(const nlohmann::json & constants) { constants_ = constants; }
     static nlohmann::json constants_;
};
}

and in PhysicalConstants.cpp

#include "PhysicalConstants.h"
#include <nlohmann/json.hpp>

namespace mrchem {
PhysicalConstants& PhysicalConstants::getInstance(const nlohmann::json & constants) {
    static PhysicalConstants obj(constants);
    return obj;
}
}

I think it is means that the private static member constants_ is declared but not defined anywhere. Perhaps I broke the class when trying to get rid of the compiler errors after copy/pasting...

stigrj commented 2 years ago

I guess this is not in the latest commit you pushed, because that one work fine for me?

Two comments:

Brakjen commented 2 years ago

@stigrj I did not push the version I does not work. The version here does not use Roberto's singleton suggestion.

I updated the python side so that only the constants we need are added to the program input. The "constants" section now looks like this:

"constants": {
      "angstrom2bohrs": 1.8897261246257702,
      "atomic_unit_of_bohr_magneton": 0.5000000000000764,
      "atomic_unit_of_nuclear_magneton": 0.0002723085107443953,
      "c_au": 137.035999084,
      "dipmom_au2debye": 2.5417464739297717,
      "electron_g_factor": -2.00231930436256,
      "fine_structure_constant": 0.0072973525693,
      "hartree2ev": 27.211386245988,
      "hartree2kJmol": 2625.4996394798254,
      "hartree2kcalmol": 627.5094740630558,
      "hartree2simagnetizability": 78.9451185,
      "hartree2wavenumbers": 219474.6313632,
      "pi": 3.141592653589793,
      "pi_sqrt": 1.7724538509055159
    }

I see no mention in the ORCA manual that changing the light speed in a relativistic calculation affects any other aspect of that code. But I am not sure what makes sense.

Brakjen commented 2 years ago

Also, how do I add qcelemental as a python dependency? As of now, the source of everything needed on the python side (that is not part of the standard library) is explicitly included, to avoid having the user install packages themselves? Do we do the same with qcelemental?

stigrj commented 2 years ago

To fix the tests on GitHub Actions I think you just have to add qcelemental to .github/mrchem-gha.yml. For the general user I'm not sure how we want to do this, perhaps it's not too much to ask of the user to install a simple python package? @robertodr ?

robertodr commented 2 years ago

It's not too onerous to install a Python package. Consider though that we removed the dependency on runtest to make development work run smoother. runtest was used to only run tests and we could rewrite that functionality in pure Python. qcelemental might be used to do more than just supply units of measure so it is a slightly different use case. If there's a way to make the dependency optional, you'd make a few developers happier 😄

stigrj commented 2 years ago

We could use qcelemental in the same way we use parselglossy, e.i. only by developers once in a while to update a JSON file with the parameters we need. Or, it could be used as part of the parselglossy workflow of generating the input template, so that default values are fetched directly from qcelemental when we run parselglossy?

robertodr commented 2 years ago

We could use qcelemental in the same way we use parselglossy, e.i. only by developers once in a while to update a JSON file with the parameters we need. Or, it could be used as part of the parselglossy workflow of generating the input template, so that default values are fetched directly from qcelemental when we run parselglossy?

sounds reasonable, yes. As long as you don't dump all constants and conversion factors the JSON shouldn't be too large.

Brakjen commented 2 years ago

@robertodr

I have updated the singleton class, but I still get problems during linking:

[ 76%] Linking CXX executable ../bin/mrchem-tests
Undefined symbols for architecture x86_64:
  "mrchem::PhysicalConstants::constants_", referenced from:
      mrchem::PhysicalConstants::PhysicalConstants(nlohmann::basic_json<std::__1::map, std::__1::vector, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, bool, long long, unsigned long long, double, std::__1::allocator, nlohmann::adl_serializer> const&) in PhysicalConstants.cpp.o
Undefined symbols for architecture x86_64:
  "mrchem::PhysicalConstants::constants_", referenced from:
      mrchem::PhysicalConstants::PhysicalConstants(nlohmann::basic_json<std::__1::map, std::__1::vector, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, bool, long long, unsigned long long, double, std::__1::allocator, nlohmann::adl_serializer> const&) in PhysicalConstants.cpp.o
ld: symbol(s) not found for architecture x86_64
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make[2]: *** [bin/mrchem.x] Error 1
make[2]: *** [bin/mrchem-pilot.x] Error 1
make[1]: *** [src/CMakeFiles/mrchem.x.dir/all] Error 2
make[1]: *** Waiting for unfinished jobs....
make[1]: *** [pilot/CMakeFiles/mrchem-pilot.x.dir/all] Error 2
Undefined symbols for architecture x86_64:
  "mrchem::PhysicalConstants::constants_", referenced from:
      mrchem::PhysicalConstants::PhysicalConstants(nlohmann::basic_json<std::__1::map, std::__1::vector, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, bool, long long, unsigned long long, double, std::__1::allocator, nlohmann::adl_serializer> const&) in PhysicalConstants.cpp.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make[2]: *** [bin/mrchem-tests] Error 1
make[1]: *** [tests/CMakeFiles/mrchem-tests.dir/all] Error 2
make: *** [all] Error 2

Do you understand why the linker complains? The updated code has been pushed.

robertodr commented 2 years ago

Not at a glance. I need to run this locally to figure it out.

Brakjen commented 2 years ago

Defining/"implementing" constants_ in PhysicalConstants.cpp seemed to do the trick:

nlohmann::json PhysicalConstants::constants_ = nlohmann::json()
Brakjen commented 2 years ago

The user can now overwrite the constant defaults in the user input, e.g.

Constants {
    pi = 3
}

The keywords to the new input section are added automatically by python/mrchem/update_input_parser.py, which reads the template and builds the Constants section based on the constants defined in MRChemPhysConstants. It also runs parselglossy and syncs the user ref file.

Brakjen commented 2 years ago

Is it normal behavior that the total energy is identical for one calculation with normal pi and one with pi = 1000000.0?

stigrj commented 2 years ago

That sounds strange. Maybe it's not used? There's also a pi in MRCPP, so maybe that is the one that is actually used?

stigrj commented 2 years ago

Come to think of it, maybe pi should be fixed and provided by MRCPP? Then we leave only physical constants to MRChem

Brakjen commented 2 years ago

HelmholtVector::apply uses pi, but this method is not currently being called during a normal SCF (we do a H() instead). pi is also referenced in Cavity.cpp and SCRF.cpp, so it might be relevant there. And some unit tests fail when drastically changing the value of pi :)

stigrj commented 2 years ago

I see MRCPP's pi is used once in tests/qmfunctions/orbital_vector.cpp. Can you just replace all the pi's with mrcpp::pi and remove it from MRChem's constants?

Brakjen commented 2 years ago

Btw, I also added a PhysicalConstants::Print(int plevel) method (being called right after initialization in mrchem.cpp). I set the print threshold to zero, but I can increase it to 1.... The printout looks like this

===========================================================================
                 Physical Constants (truncated precision)
---------------------------------------------------------------------------
angstrom2bohrs                        =  1.88972612e+00
dipmom_au2debye                       =  2.54174647e+00
electron_g_factor                     = -2.00231930e+00
fine_structure_constant               =  7.29735257e-03
hartree2ev                            =  2.72113862e+01
hartree2kcalmol                       =  6.27509474e+02
hartree2kjmol                         =  2.62549964e+03
hartree2simagnetizability             =  7.89451185e+01
hartree2wavenumbers                   =  2.19474631e+05
light_speed                           =  1.37035999e+02
===========================================================================
Brakjen commented 2 years ago

I think all the suggested edits have been made.