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

gitpeterwind · 2024-12-30T10:32:22Z

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 MRChem PR.

…exTrees

robertodr · 2025-01-18T13:51:12Z

src/treebuilders/SquareCalculator.h


-    void calcNode(MWNode<D> &node_o) override {
+    void calcNode(MWNode<D, double> &node_o) {


Same comment here as well

robertodr · 2025-01-18T13:53:50Z

src/treebuilders/add.cpp

 }

-template void add<1>(double prec,
-                     FunctionTree<1> &out,
+template void add<1, double>(double prec,


yes, a macro would go a very long way to reduce the boilerplate also for these instantiations!

robertodr · 2025-01-18T13:59:05Z

src/treebuilders/apply.cpp

+            else defaultMetric[i][j] = 0.0;
+        }
+    }
+    if (metric == nullptr) {


C++17 introduced std::optional to deal with arguments that may or may not be passed. We should use that here.

robertodr · 2025-01-18T14:00:41Z

src/treebuilders/apply.cpp

+            else defaultMetric[i][j] = 0.0;
+        }
+    }
+    if (metric == nullptr) {


Same comment about std::optional applies also here.

robertodr · 2025-01-18T14:01:02Z

src/treebuilders/apply.cpp

+            else defaultMetric[i][j] = 0.0;
+        }
+    }
+    if (metric == nullptr) {


std::optional also here.

robertodr · 2025-01-18T14:01:35Z

src/treebuilders/apply.cpp

@@ -308,6 +461,46 @@ template <int D> void apply(FunctionTree<D> &out, DerivativeOperator<D> &oper, F
    print::separator(10, ' ');
 }

+template <int D> void apply(CompFunction<D> &out, DerivativeOperator<D> &oper, CompFunction<D> &inp, int dir, ComplexDouble metric[4][4]) {
+    //TODO: sums and not only each components independently


TODO still relevant?

Yes, but I wait for the implementation until we have real examples: the solution could take several forms (when to introduce the multiplicative constants and the like). (I completed the comment)

robertodr · 2025-01-18T14:01:52Z

src/treebuilders/apply.cpp

+                    apply(*out.CompD[ocomp], oper, *inp.CompD[icomp], dir);
+                    if (std::norm(metric[icomp][ocomp] - 1.0) > MachinePrec) {
+                        if(std::imag(metric[icomp][ocomp]) < MachinePrec) out.CompD[ocomp]->rescale(std::real(metric[icomp][ocomp]));
+                        else out.func_ptr->data.c1[ocomp] *= metric[icomp][ocomp]; //TODO: multiply c1 in rescale?


TODO still relevant?

I rephrased: // To consider: multiply c1 in rescale?
The "c1 constant" is somewhat new and it is not yet completely integrated into the logic. It is a "soft" multiplicative constant; soft meaning the same as the conjugate flag for a function: not explicitly applied before it is required.
So far it is used for the complex spin operators: we do not want to convert all real functions to complex just because there is a complex multiplicative factor in the operator.

robertodr · 2025-01-18T14:09:12Z

src/treebuilders/map.cpp

@@ -66,12 +66,12 @@ namespace mrcpp {
 *
 */
 template <int D>
-void map(double prec, FunctionTree<D> &out, FunctionTree<D> &inp, FMap fmap, int maxIter, bool absPrec) {
+void map(double prec, FunctionTree<D, double> &out, FunctionTree<D, double> &inp, FMap<double, double> fmap, int maxIter, bool absPrec) {


is map not relevant when one has float or std::complex scalar types?

I don't know, but since it is not used I prefer to not implement before we know what we want to use it for.

robertodr · 2025-01-18T14:16:28Z

src/treebuilders/multiply.h

+template <int D> ComplexDouble dot(FunctionTree<D, ComplexDouble> &bra,
+                            FunctionTree<D, double> &ket);
+template <int D> ComplexDouble dot(FunctionTree<D, double> &bra,
+                            FunctionTree<D, ComplexDouble> &ket);


I think it's enough to have one function, but with three template parameters (similarly could be done for all arithmetic functions, actually)

Suggested change

template <int D> ComplexDouble dot(FunctionTree<D, ComplexDouble> &bra,

FunctionTree<D, double> &ket);

template <int D> ComplexDouble dot(FunctionTree<D, double> &bra,

FunctionTree<D, ComplexDouble> &ket);

template <int D, typename T, typename U = T, typename V = decltype(std::declval<T>() * std::declval<U>())>

V dot(FunctionTree<D, T> &bra, FunctionTree<D, U> &ket);

you can see why I dislike implementing templates in .cpp files 😅 (cough cough @ilfreddy)

Note that in the C++ standard mixing float and std::complex<double> (or double and std::complex<float>) is not possible, so the above should be combined with a compile-time exclusion of such a possibility (through SFINAE or something) Another thing I have in my archives.

This is nice!

robertodr · 2025-01-18T14:19:58Z

src/trees/FunctionNode.h

-
-template <int D> double dot_scaling(const FunctionNode<D> &bra, const FunctionNode<D> &ket);
-template <int D> double dot_wavelet(const FunctionNode<D> &bra, const FunctionNode<D> &ket);
+template <int D>


again, I think these signatures can be declared only once using the "trick" above.

I did not use the trick for the dot_scaling and dot_wavelet, because with the conjugate flags it became messy.

robertodr · 2025-01-18T14:22:57Z

src/trees/FunctionTree.cpp

+/** @brief Read a previously stusing MADNESS conventions for n, l and index order.ored tree assuming text/ASCII format,
+ *   in a representation
+ * @param[in] file: File name
+ * @note This tree must have the exact same MRA the one that was saved(?)
+ */


The docstring here seems a bit mangled. Can you fix it?

robertodr · 2025-01-18T14:24:19Z

src/trees/FunctionTree.cpp

+        if (std::abs(coord[d][0] + L) > TXT_thres) std::cout<<coord[d][0]<<" "<<L<<std::endl;;
+        if (std::abs(coord[d][0] + L) > TXT_thres) NOT_IMPLEMENTED_ABORT;
+        if (std::abs(coord[d][1] - L) > TXT_thres) std::cout<<coord[d][1]<<" "<<L<<std::endl;;
+        if (std::abs(coord[d][1] - L) > TXT_thres) NOT_IMPLEMENTED_ABORT;


I don't understand why the 2nd and 4th statements

It did not seem to work when the "scaling factors" where not 1.0. At present the box must be with sizes which are powers of 2.

robertodr · 2025-01-18T14:27:54Z

src/trees/FunctionTree.cpp

                                      std::vector<int> &indices,
                                      std::vector<int> &parent_indices,
                                      std::vector<double> &scalefac,
                                      int &max_index,
-                                      MWTree<D> &refTree,
-                                      std::vector<MWNode<D> *> *refNodes) {
+                                      MWTree<D, double> &refTree,


here (and some other places) you explicitly fix the scalar type of some function inputs. Why is that?

robertodr · 2025-01-18T14:30:02Z

src/trees/FunctionTree.cpp

+    return out;
+}
+
+    /*


I think it's fine to delete commented-out code.

robertodr · 2025-01-18T14:35:20Z

src/trees/MWNode.cpp

-#else
-    for (int i = start; i < start + size; i++) { sq_norm += c[i] * c[i]; }
-#endif
+//#ifdef HAVE_BLAS


commented-out code can be removed

robertodr · 2025-01-18T14:38:12Z

src/trees/MultiResolutionAnalysis.cpp

@@ -146,7 +146,12 @@ template <int D> bool MultiResolutionAnalysis<D>::operator==(const MultiResoluti
 *  For more information about the meaning of equality for BoundingBox and ScalingBasis objets, see their respective classes.
 */
 template <int D> bool MultiResolutionAnalysis<D>::operator!=(const MultiResolutionAnalysis<D> &mra) const {
-    return !(*this == mra);
+    if (this->basis != mra.basis) std::cout<<"diff basis "<<this->basis<<std::endl <<"and  "<< mra.basis<<std::endl;


are the print statements needed here? if yes, maybe the should use the printer. Also, no need to check the condition twice, just do a non-inline if-statement.

robertodr

Can you run clang-format on all the C++ files this PR changes?
When giving a default for the scalar type (e.g. typename T = double) I think it should be enough to do it on the parent classes. This is however just to save some typing :)
I think we should have compile-time checks on whether the scalar type T makes sense or not, e.g. only one of float, double, std::complex<float>, std::complex<double>. I have some code to do this that we can reuse.
For the metric arguments of some functions in apply we may use std::optional<Eigen::Matrix<std::complex<double>, 4, 4>>.
There's a lot more boilerplate in the explicit instantiations. The obvious solution is to not put the implementation in .cpp files, but leave everything in the header. @ilfreddy (and maybe @stigrj?) feels strongly against it, so we'll have to make do with some preprocessor macros and/or some metaprogramming tricks.
I think we should allow to mix scalar types, but it can be left for a future refinement, as it involves some tricky gymnastics with template parameters.

robertodr · 2025-01-18T14:51:42Z

src/utils/math_utils.h

@@ -67,6 +67,7 @@ double matrix_norm_1(const Eigen::MatrixXd &M);
 double matrix_norm_2(const Eigen::MatrixXd &M);

 void apply_filter(double *out, double *in, const Eigen::MatrixXd &filter, int kp1, int kp1_dm1, double fac);
+void apply_filter(ComplexDouble *out, ComplexDouble *in, const Eigen::MatrixXd &filter, int kp1, int kp1_dm1, double fac);


I think the two can be unified with if constexpr

robertodr

in parallel.h/parallel.cpp I think the implementations for double and std::complex<double> could be unified, so we reduce a bit the boilerplate (and possibly extend to float)

gitpeterwind · 2025-01-28T11:04:01Z

I have now included all the changes you proposed, or left a comment.
It is a pity that I did not use the if constexpr from the start. Actually I did, but then went back because it was not compliant with the std14.
About using float I am sceptical: for small system there is no point, and for large "HPC" systems the error would be too large. Roughly a HPC calculation would take 10¹² operation at least, giving sqrt(10⁻¹²) = 10⁻⁶ relative numerical error, 6 digits error of about 7 significant digits is not useful.
About reducing the boilerplate instantiations, I do not think it is important: it is not those parts that take time to implement, and they do not reduce code readability (while macros might!).

ilfreddy · 2025-01-28T12:58:11Z

Does this include BOTH the complex number template AND the MADNESS input stuff?

gitpeterwind · 2025-01-28T12:59:51Z

yes

gitpeterwind and others added 30 commits April 25, 2024 19:21

template for native complex trees

56c96dc

update examples for native complex trees

0760089

Merge branch 'master' of github.com:MRChemSoft/mrcpp into NativeCompl…

b8c9085

…exTrees

remove Complex template from Gaussians

f271ea9

test with derivative of Complex Exponential

2627f64

Merge branch 'master' into NativeComplexTrees

365bd0c

mrchem compatibility

fa478e8

Four components class and apply

a595da6

More native Complex tree operations: add, mult, Bank, send, rescale

303a545

Merge branch 'NativeComplexTrees' into FourComponents

daaed5c

FunctionnTree deepcopy, addinplace, creator from components function

fbc96f7

Component Functions: mpi send and recv

3aa778a

Component Functions: norms

9ce1b76

CompFunction destructor

3057d0a

FourComp for share, add, multiply rescale, dot

9b24281

removed typename template and upgrades

f734374

moved ComplexFunctions.h from parallel.cpp to Bank.h

47471f8

removed the mpi:: in parallel.cpp

aa6f254

Function vector with pointers

52a6d52

explicit real and complex reduce and broadcast

0971589

all functions upgraded to component functions

6e6e881

added #ifdef MRCPP_HAS_MPI in parallel.cpp

da6e5e9

moved #ifdef MRCPP_HAS_MPI in parallel.cpp

6ffb77f

compiles together with UseComponent branch of mrchem

36ee92b

removed ComplexFunction from api

0df20bc

shared CompFunctions trees

1e33df7

dagger, serial rotate and overlap

e95e128

multiply, apply, rescale, mixed types

5d8ee63

orthogonalize for complex type

6624a73

all tests passes

1b29dea

robertodr reviewed Jan 18, 2025

View reviewed changes

src/treebuilders/apply.cpp Outdated

else defaultMetric[i][j] = 0.0;

}

}

if (metric == nullptr) {

Copy link

Contributor

robertodr Jan 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

std::optional also here.

robertodr reviewed Jan 18, 2025

View reviewed changes

src/trees/FunctionTree.cpp Outdated

return out;

}

/*

Copy link

Contributor

robertodr Jan 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's fine to delete commented-out code.

robertodr reviewed Jan 18, 2025

View reviewed changes

robertodr requested changes Jan 18, 2025

View reviewed changes

robertodr reviewed Jan 18, 2025

View reviewed changes

gitpeterwind added 4 commits January 24, 2025 12:22

clang-format

f1a5ef4

removed the defaultMetric from apply.cpp

71302e5

std17 and PR review requests

da5e41f

std17 and PR review requests

71d6e2d

std17 and PR review requests

f46ac4e

gitpeterwind mentioned this pull request Jan 28, 2025

template based native complex trees #244

Closed

ilfreddy mentioned this pull request Jan 28, 2025

Four components #245

Closed

Merge branch 'master' into text

9b45d23

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

gitpeterwind commented Dec 30, 2024

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

gitpeterwind Jan 27, 2025

robertodr Jan 18, 2025

gitpeterwind Jan 27, 2025

robertodr Jan 18, 2025

gitpeterwind Jan 27, 2025

robertodr Jan 18, 2025

gitpeterwind Jan 27, 2025

robertodr Jan 18, 2025

gitpeterwind Jan 27, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

gitpeterwind Jan 27, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr Jan 18, 2025

robertodr left a comment

robertodr Jan 18, 2025

robertodr left a comment

gitpeterwind commented Jan 28, 2025

ilfreddy commented Jan 28, 2025

gitpeterwind commented Jan 28, 2025


		void calcNode(MWNode<D> &node_o) override {
		void calcNode(MWNode<D, double> &node_o) {

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

Are you sure you want to change the base?

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

Conversation

gitpeterwind commented Dec 30, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

robertodr left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

robertodr left a comment

Choose a reason for hiding this comment

gitpeterwind commented Jan 28, 2025

ilfreddy commented Jan 28, 2025

gitpeterwind commented Jan 28, 2025