Save or read orbitals in text format, portable and readable by MADNESS group #498

gitpeterwind · 2024-12-30T10:33:59Z

Orbitals can be saved and read in text format (ascii) files. Each end node is saved with info about depth and translation.
Use keyword write_orbitals_txt = T to save in this format.

Also introduced keywords bank_per_node and use_omp_num_threads after feedback from Moritz and Roberto. (see documentation)

NB: this PR builds on top of the Four Components PR and is to be used with the corresponding MRCPP PR. It is merged with the latest master branch.

Close MRChemSoft#419 We now use qcelemental to fetch Bondi radii, and use alpha=1.2 for hydrogen --------- Co-authored-by: Roberto Di Remigio Eikås <[email protected]>

) * cavity function used directly instead of projected first * clear d_V * Update SCRF.cpp define ranks in accelerator, so that all mpi ranks use kain * cavity function used directly instead of projected first * clear d_V * Update SCRF.cpp define ranks in accelerator, so that all mpi ranks use kain * manage to pull down to 12.2 and 30.0 gb * reduced memory usage to ca. 3 gb from 8 gb * Remove print statement for memory * Small changes to desctructor and clear * Use latest MRCPP * Use default initialization in header file * Do some cleaning after feedback * Remove ´optimizer´ option from inputs * Fix tests --------- Co-authored-by: Roberto Di Remigio Eikås <[email protected]> Co-authored-by: Gabrielgerez <[email protected]>

…emSoft#473) * Convert default radii to bohr when input units are not angstrom * Update template.yml (formatting and documentation) * Refine fix for unit conversion of cavity params, add tests

* Refactor the Permittivity function for ease of future implementation work. * Add ShiftFunction to CMakeLists * Add virtual and override where necessary * Remove evalf definition from shiftFunction * Rename ShiftFunction to StepFunction * Update src/environment/Permittivity.cpp Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Update src/environment/StepFunction.cpp Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Update src/environment/StepFunction.cpp Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Update src/environment/StepFunction.cpp Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Update src/environment/StepFunction.h Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Update src/environment/StepFunction.h Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Pass cavity as a shared_ptr to save memory. * Fix rebase error * Remove use for `flipFunction` in Permittivity * Update src/environment/Permittivity.h Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Update src/environment/StepFunction.cpp Co-authored-by: Roberto Di Remigio Eikås <[email protected]> * Add print_header to detail namespace * Update src/environment/StepFunction.h Co-authored-by: Roberto Di Remigio Eikås <[email protected]> --------- Co-authored-by: Roberto Di Remigio Eikås <[email protected]>

* Use C++17 * Use std::tie to unpack tuple return value * Remove some cruft from Cavity and Permittivity * SCRF deals in densities only, not orbitals * Refactor Reaction{Potential,Operator} following Coulomb example In preparation for the computation of the response terms, we make ReactionPotential abstract and add subclass ReactionPotentialD1 to deal with one single set of orbitals. * Modify input to accept static and dynamic permittivities * Handle nonequilibrium solvation more gracefully * Sign change was necessary * Update input parser * Small fixes to include order * Rename section for permittivity outside * Add tests

in json form

robertodr · 2025-02-25T13:46:17Z

src/qmfunctions/Orbital.h

+#define spin() func_ptr->data.n1[0]
+#define occ() func_ptr->data.d1[0]


I don't understand the need for these two?

I didn't manage to use directly a function for the Orbital class, because if Phi is an OrbitalVector, Phi[i] returns not an Orbital, but a CompFunction. (OrbitalVector is class derived from a vector of CompFunction)
Phi[i].occ() = 0 for example should work. (I even tried to redefine [] so that it returns an Orbital)

src/qmfunctions/density_utils.cpp

robertodr · 2025-02-25T14:01:53Z

src/qmfunctions/orbital_utils.h

@@ -53,7 +49,7 @@ OrbitalVector param_copy(const OrbitalVector &Phi);
 OrbitalVector adjoin(OrbitalVector &Phi_a, OrbitalVector &Phi_b);
 OrbitalVector disjoin(OrbitalVector &Phi, int spin);

-void save_orbitals(OrbitalVector &Phi, const std::string &file, int spin = -1);
+void save_orbitals(OrbitalVector &Phi, const std::string &file, int spin = -1, int text_format = 0);


why not a bool?

Yes, could have used a bool, but I though that it might be useful to have an int to, in the future, distinguish between formats (MADNESS/MRChem for example)

robertodr · 2025-02-25T14:02:30Z

src/qmfunctions/orbital_utils.h

 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);
 int print_size_nodes(const OrbitalVector &Phi, const std::string &txt = "", bool all = true, int plevel = 0);
-void saveOrbital(const std::string &file, Orbital& orb);
+void saveOrbital(const std::string &file, const Orbital& orb, int text_format = 0);


why not bool for text_format?

robertodr · 2025-02-25T14:03:03Z

src/qmfunctions/density_utils.h

@@ -36,6 +36,7 @@ void allreduce_density(double prec, Density &rho_tot, Density &rho_loc);
 void compute(double prec, Density &rho, mrcpp::GaussExp<3> &dens_exp);
 void compute(double prec, Density &rho, OrbitalVector &Phi, DensityType spin);
 void compute(double prec, Density &rho, OrbitalVector &Phi, OrbitalVector &X, OrbitalVector &Y, DensityType spin);
+    //double compute_occupation(const Orbital &phi, DensityType dens_spin);


can the commented-out code be deleted?

robertodr · 2025-02-25T14:10:22Z

src/qmfunctions/qmfunction_utils.cpp_notused

should this file be deleted?

robertodr · 2025-02-25T14:10:43Z

src/qmfunctions/qmfunction_utils.h_notused

maybe this can be deleted?

robertodr · 2025-02-25T14:16:58Z

src/qmoperators/QMPotential.cpp22aug

is this needed?

robertodr · 2025-02-25T14:20:35Z

src/qmoperators/one_electron/PositionOperator.h

-        mrcpp::cplxfunc::project(*r_x, f_x, NUMBER::Real, -1.0);
-        mrcpp::cplxfunc::project(*r_y, f_y, NUMBER::Real, -1.0);
-        mrcpp::cplxfunc::project(*r_z, f_z, NUMBER::Real, -1.0);
+        // the cast is required, because otherwise the compiler is not smart enough to understand that it returns a double


ok, wow. that's one ugly line of code the compiler forces you to write 😸 Is it because project is templated on the scalar type now? Could you post the compiler error?

/home/peter/mrchem/src/qmoperators/one_electron/PositionOperator.h:54:23: error: call of overloaded ‘project(std::__shared_ptr_access<mrchem::QMPotential, __gnu_cxx::_S_atomic, false, false>::element_type&, mrchem::PositionOperator::PositionOperator(mrcpp::Coord<3>&)::<lambda(mrcpp::Coord<3>&)>&, double)’ is ambiguous
54 | mrcpp::project(*r_x, f_x, -1.0);
| ~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~

robertodr · 2025-02-25T14:22:30Z

src/qmoperators/two_electron/ExchangePotential.cpp

+    mrcpp::FunctionTreeVector<3, double> phi_opt_vec_real;
+    mrcpp::FunctionTreeVector<3, ComplexDouble> phi_opt_vec_cplx;
+    if (RealOrbitals) {
+        if (phi_k.isreal() and phi_k.getNNodes() > 0) phi_opt_vec_real.push_back(std::make_tuple(1.0, phi_k.CompD[0]));


it should be possible to do this instead. With emplace_back the compiler should avoid temporaries:

Suggested change

if (phi_k.isreal() and phi_k.getNNodes() > 0) phi_opt_vec_real.push_back(std::make_tuple(1.0, phi_k.CompD[0]));

if (phi_k.isreal() and phi_k.getNNodes() > 0) phi_opt_vec_real.emplace_back(1.0, phi_k.CompD[0]);

robertodr · 2025-02-25T14:35:17Z

src/tensor/RankZeroOperator.cpp22aug

is this needed?

robertodr · 2025-02-25T14:35:49Z

src/tensor/tt

this seems to be a leftover?

robertodr · 2025-02-25T14:38:59Z

tests/qmfunctions/orbital.cpp

@@ -33,12 +33,13 @@ using namespace mrchem;

 namespace orbital_tests {

-auto f = [](const mrcpp::Coord<3> &r) -> double {
+std::function<double(const mrcpp::Coord<3> &r)> f = [](const mrcpp::Coord<3> &r) -> double {


Going from auto to explicit signatures for a lambda function is a step backwards. If you can't use auto then there's likely a problem in the code that should be fixed.

The error with auto is:
/home/peter/mrchem/tests/qmfunctions/orbital.cpp:122:27: error: call of overloaded ‘project(mrchem::Orbital&, orbital_tests::<lambda(mrcpp::Coord<3>&)>&, const double&)’ is ambiguous
122 | mrcpp::project(phi_2, g, prec);

i.e. the compiler is not smart enough.

And I do not agree about auto: when coding it is frustrating to not know immediately which type a variable is, and auto should be avoided as much as possible.

robertodr · 2025-02-25T14:40:22Z

tests/qmfunctions/qmfunction.cpp

        }
+        // SECTION("unitary rotation") {


why is this test commented out?

It does not make sense anymore in the new context: we do not have function with a real and imaginary part anymore, but real and complex functions.
Removed the uncommented part.

robertodr

Commented-out code can be deleted
The de-auto-ization when defining the lambdas in the tests points to some issue in the code.

gitpeterwind · 2025-02-26T09:28:14Z

All the comments were either implemented or commented (or both)

ilfreddy · 2025-02-27T08:42:00Z

tests/qmoperators/xc_hessian_pbe.cpp

The numerical differences here seem a bit large. We need to understand why it is the case.

There was a difference in the setting for the multiplication in the density

tests/qmoperators/xc_operator_blyp.cpp

ilfreddy · 2025-02-27T08:44:08Z

tests/qmoperators/xc_operator_blyp.cpp

The numerical differences seem a bit too large for comfort here as well. It might have to do with gradients (LDA was OK)

There was a difference in the setting for the multiplication in the density

Co-authored-by: Luca Frediani <[email protected]>

Gabrielgerez and others added 30 commits July 11, 2024 09:31

Fix default vdW radii (MRChemSoft#463)

328bb23

Close MRChemSoft#419 We now use qcelemental to fetch Bondi radii, and use alpha=1.2 for hydrogen --------- Co-authored-by: Roberto Di Remigio Eikås <[email protected]>

add readthedocs.yml, same as for MRCPP

6a5a348

Manually add requirements to .readthedocs.yml

096725a

attempt 2 at setting up requirements

d5becc8

Fix typo in .readthedocs.yml

bb93e8b

Fix a second typo in .readthedocs.yml

6c36768

Fix incoming build error in rtd

7405a5f

Update .readthedocs.yml

b4c5fa4

Pin sphinx

8cfc889

Update betzy_example.job

7f73667

Convert default radii to bohr when input units are not angstrom (MRCh…

2c78b6c

…emSoft#473) * Convert default radii to bohr when input units are not angstrom * Update template.yml (formatting and documentation) * Refine fix for unit conversion of cavity params, add tests

Bugfix: critical sign error in magnetic properties

9849caa

Add test for magnetic properties

f17aa5a

Validate cavity only if doing PCM

9dcd227

Tell user to specify radii explicitly when not known

5366f2b

add sqnm code

cc9733c

add readme of vcsqnm package

3092b7f

paritally finished structure of interface

341cdf9

saves

0315e90

improve print statement

ceea49e

sketch of geopt working

f24ff5c

remove unnecessary prints, adjust geopt params.

a1e3f5f

read in old previous solution as initial guess

d017dac

add parameters, mpi barrier and return final result

00eeb6f

in json form

update intrface to geompetry optimizer

b99815b

added documentation and more output.

acefaff

added mpi todo.

edabca3