Merge pull request #4 from eljost/graphs

Graphs & sympy patch

.gitignore

@@ -129,3 +129,11 @@ dmypy.json
 .pyre/
 
 version.py
+# Calculated graphs
+*.xml.tar.gz
+
+*.o
+*.mod
+*.png
+*.svg
+*.f90

README.md

@@ -1,13 +1,33 @@
 # sympleints
 Molecular integrals over **shells** of Gaussian basis functions using sympy.
 
-By providing the appropriate commands, `sympleints` is able to generate integral code that already includes
-the Cartesian-to-spherical transformation (`--sph`) and normalization factors (`--norm-pgto`).
+By providing the appropriate commands, `sympleints` is able to generate integral code that
+already includes the Cartesian-to-spherical transformation (`--sph`) and/or normalization
+factors (`--normalize pgto|cgto|none`).
+
+A sensible way to handle normalization is outlined in the notebook
+`ressources/cgto_normalization.ipynb`.
 
 The generated python code makes heavy use of `numpy` and by providing arguments with appropriate
 shapes, all functions can be evaluted using arrays of exponents and contraction coefficients
 (contraced gaussian functions).
 
+Presently, the library implements two approaches:
+
+1. Explicit generation of full integral expressions via explicit implementation of reccurence
+   relations, followed by common-subexpression-elimination.
+2. Generation of dependence graphs for integrals and their translation into actual code.
+
+Both approaches have their advantages and drawbacks.
+
+1. Programming the reccurence relations is cumbersome, having an explicit expression for the integral at hand is very  
+    convenient. Generation of kinetic energy integrals is very simple, as the appropriate products of overalp and intermediate  
+    kinetic energy integrals are easily multiplied using sympy. Derivatives are probably easily obtained by differentiation, but I  
+    never explored this. The resulting code is often slow(er), compared to the second approach.
+2. Formulating the reccurence relations is very easy and the resulting code can be much faster. Currently, the actual code generation  
+     is less automated compared to the first approach, but one also has greater flexibility. Right now, the second approach yields only  
+     Fortran code, but extension to other languages would be trivial, as long as there is an associated sympy-Printer.
+
 ## Installation
 ```bash
 git clone git@github.com:eljost/sympleints.git
@@ -19,9 +39,11 @@ After successful installation the `sympleints` command should be available in th
 
 ## Example usage
 
+Currently, only the first approach is exposed via the command line interface.
+
 ```bash
 $ sympleints --help
-usage: sympleints [-h] [--lmax LMAX] [--lauxmax LAUXMAX] [--write] [--out-dir OUT_DIR] [--keys KEYS [KEYS ...]] [--sph] [--norm-pgto]
+usage: sympleints [-h] [--lmax LMAX] [--lauxmax LAUXMAX] [--write] [--out-dir OUT_DIR] [--keys KEYS [KEYS ...]] [--sph] [--opt-basic] [--normalize {pgto,cgto,none}]
 
 options:
   -h, --help            show this help message and exit
@@ -32,7 +54,8 @@ options:
   --keys KEYS [KEYS ...]
                         Generate only certain expressions. Possible keys are: (gto, ovlp, dpm, dqpm, qpm, kin, coul, 2c2e, 3c2e_sph). If not given, all expressions are generated.
   --sph
-  --norm-pgto
+  --opt-basic           Turn on basic optimizations in CSE.
+  --normalize {pgto,cgto,none}
 ```
 
 To get a quick impression some overlap integrals up to p-functions are generated via:
@@ -62,11 +85,11 @@ this, please see the module `pysisyphus.wavefunction.shells` in the
     3. Quadratic moment (quadrupole moment) integrals (order 2)
     4. Extension to higher multipoles would be trivial
   2. Kinetic energy integrals
-  3. Nuclear attraction integrals
+  3. Nuclear attraction integrals  (avaiable via first and second approach)
 
 ### 2-electron integrals
-  1. 2-center-2-electron integrals
-  2. 3-center-2-electron-integrals
+  1. 2-center-2-electron integrals (avaiable via first and second approach)
+  2. 3-center-2-electron-integrals (avaiable via first and second approach)
 
 Together, these two types of 2-electron integrals allow the implementation of density fitting.
 
@@ -76,7 +99,8 @@ of `sympleints`. Suitable code for the Boys-function is found in the
 
 ## Advantages
 The resulting Python code is surely not optimal, but very convenient to have and easy to use. Besides an implementation
-of the Boys-function, the resulting code only depends numpy. Understanding the integral generator is still possible.
+of the Boys-function, the resulting code only depends numpy. Understanding the integral generator is still possible for a
+simple minded python programmer as myself.
 
 Calculation of the overlap matrix with 1300 spherical basis functions (Tris(bipyridine)ruthenium(II))
 takes 9.5 s, calculation of the quadratic moment integrals requires 20.2 s. Maybe, these timings make the previous sentence more suitable for its placement in the Limiations section ;)

setup.cfg

@@ -23,13 +23,15 @@ install_requires =
     colorama
     fprettify
     jinja2
+    matplotlib
     networkx
     numpy
     pytest
-    sympy==1.10.1  # Currently fails with 1.11 because of GH issue #24236
+    sympy
 
 
 [options.entry_points]
 console_scripts =
     sympleints = sympleints.main:run_cli
+    sympleints-graph = sympleints.graphs.main:run
     symplebench = sympleints.benchmark:run

sympleints/FortranRenderer.py

@@ -1,6 +1,7 @@
 import os
 import re
 import tempfile
+import warnings
 
 from sympy.codegen.ast import Assignment
 from sympy.printing.fortran import FCodePrinter
@@ -66,6 +67,20 @@ def _print_IndexedSlice(self, expr):
         return f"{expr.org_name}({start}:{stop})"
 
 
+def get_fortran_print_func(**print_settings):
+    # This allows using the 'boys' function without producing an error
+    _print_settings = {
+        "allow_unknown_functions": True,
+        # Without disabling contract some expressions will raise ValueError.
+        "contract": False,
+        "standard": 2008,
+        "source_format": "free",
+    }
+    _print_settings.update(print_settings)
+    print_func = FCodePrinterMod(_print_settings).doprint
+    return print_func
+
+
 def make_fortran_comment(comment_str):
     return "".join([f"! {line}\n" for line in comment_str.strip().split("\n")])
 
@@ -103,15 +118,10 @@ def get_argument_declaration(self, functions, contracted=False):
     def render_function(
         self, functions, repls, reduced, shape, shape_iter, args, name, doc_str=""
     ):
-        # This allows using the 'boys' function without producing an error
-        print_settings = {
-            "allow_unknown_functions": True,
-            # Without disabling contract some expressions will raise ValueError.
-            "contract": False,
-            "standard": 2008,
-            "source_format": "free",
-        }
-        print_func = FCodePrinterMod(print_settings).doprint
+        if functions.primitive:
+            warnings.warn("primitive=True is currently ignored by FortranRenderer!")
+
+        print_func = get_fortran_print_func()
         assignments = [Assignment(lhs, rhs) for lhs, rhs in repls]
         repl_lines = [print_func(as_) for as_ in assignments]
         results = [print_func(red) for red in reduced]

sympleints/Functions.py

@@ -23,6 +23,7 @@ class Functions:
     full_name: Optional["str"] = None
     l_aux_max: Optional[int] = None
     spherical: bool = False
+    primitive: bool = False
 
     def __post_init__(self):
         assert self.l_max >= 0

sympleints/PythonRenderer.py

@@ -32,7 +32,8 @@ def render_function(
         # Here, we expect the orbital exponents and the contraction coefficients
         # to be 2d/3d/... numpy arrays. Then we can utilize array broadcasting
         # to evalute the integrals over products of primitive basis functions.
-        result_lines = [f"numpy.sum({line})" for line in result_lines]
+        if not functions.primitive:
+            result_lines = [f"numpy.sum({line})" for line in result_lines]
         # Drop ncomponents for simple integrals, as the python code can deal with
         # contracted GTOs via array broadcasting.
         if functions.ncomponents == 1:

sympleints/__init__.py

@@ -1,3 +1,11 @@
+import warnings
+
+from sympleints.patch_sympy import patch_assignment_printing
+
+patch_assignment_printing()
+warnings.warn("Monkeypatched NumpyPrinter._print_Assignment!")
+
+
 from sympleints.helpers import (
     canonical_order,
     get_center,

sympleints/config.py

@@ -1,3 +1,9 @@
+from pathlib import Path
+
+
 L_MAX = 4  # Up to G by default
 L_AUX_MAX = 6  # Up to I by default
 PREC = 16
+# Base name for arrays in generated code
+ARR_BASE_NAME = "work"
+CACHE_DIR = Path(".sympleints")

sympleints/defs/multipole.py

@@ -1,8 +1,9 @@
 import functools
 
-from sympy import pi, sqrt
+from sympy import pi, sqrt, S
 
 from sympleints import shell_iter
+from sympleints.sym_solid_harmonics import mOrder, Rlm_polys_up_to_l_iter
 from sympleints.defs import TwoCenter1d, RecurStrategy, Strategy
 
 
@@ -54,3 +55,84 @@ def gen_diag_quadrupole_shell(La_tot, Lb_tot, a, b, A, B, R=(0.0, 0.0, 0.0)):
             exprs.append(gen_multipole_3d(La, Lb, a, b, A, B, Le, R))
     lmns = lmns + lmns + lmns
     return exprs, lmns
+
+
+def gen_multipole_sph_shell(La_tot, Lb_tot, a, b, A, B, Le_tot=2):
+    """Integrals up to spherical quadrupoles.
+
+    The minus is NOT taken into account here, but must be taken into account
+    in the code, that utilizes these integrals.
+
+    Different approaches are possible:
+
+    Take one origin argument and calculate all multipoles w.r.t this origin,
+    or use a different origin for all primitive pairs, and shift later.
+    """
+    # Here, the multipole origin is the natural center of the Gaussian overlap
+    # distribution and depends on the orbital exponents.
+    R = (a * A + b * B) / (a + b)
+
+    """Step 1:
+
+    Generate the necessary regular solid harmonic operators. We do this
+    by generating the appropriate regular solid harmonic polynomials.
+    They can be converted to sympy.Poly-objects in the basis of x, y and z.
+
+    From these Poly-objects the required coefficients to combine the original
+    multipole integrals are easily extracted."""
+
+    # m runs from -l to +l. This is not Stone's ordering, which is 0, +1, -1,
+    # +2, -2, ...!
+    polys, lms = zip(*Rlm_polys_up_to_l_iter(Le_tot, order=mOrder.NATURAL))
+    poly_terms = [poly.terms() for poly in polys]
+
+    # Gather unique multipole-operators from all terms (Le values)
+    unique_Le_operators = set()
+    L_tot = La_tot + Lb_tot
+    for pterms in poly_terms:
+        for Le, _ in pterms:
+            unique_Le_operators.add(Le)
+
+    Le_exprs = dict()
+
+    lmns = list(shell_iter((La_tot, Lb_tot)))
+    # Generate the required basic multipole integrals and store the expressions in
+    # a dictionary, with the operator-tuple as key.
+    for Le in unique_Le_operators:
+        for La, Lb in lmns:
+            Le_expr = gen_multipole_3d(La, Lb, a, b, A, B, Le, R)
+            Le_exprs.setdefault(Le, list()).append(Le_expr)
+
+    # Combine the basic multipole integrals into final integrals over spherical
+    # multipole operators.
+    exprs = list()
+    all_lmns = list()
+
+    # Loop over all operators that we wan't to generate.
+    for pterms in poly_terms:
+        # Every spherical multipole integral is the sum of multiple basic multipole integrals.
+        # Depending on the angular momenta of the involved basis functions some integrals
+        # may be 0, so we initialize all expressions with 0.
+        tmp_exprs = [S.Zero] * len(lmns)
+
+        # Don't generate unneccesary expressions, as the product of two primitive
+        # Gaussians with angular momenta La and Lb gives rise to multipoles up to order
+        # La + Lb at most.
+        # So, 2 s-orbitals produce a charge, a s- and a p-orbital give rise to a charge
+        # and a dipole moment etc.
+        Le0, _ = pterms[0]
+        Le0 = sum(Le0)
+        if Le0 > L_tot:
+            exprs.extend(tmp_exprs)
+            all_lmns += lmns
+            continue
+
+        # Build up the final expressions. Every operator can comprise multiple basic
+        # integrals. Here, we loop over all terms in the polynomial and add the respective
+        # basic integrals, multiplied by the appropriate coefficient.
+        for Le, coeff in pterms:
+            for i, expr in enumerate(Le_exprs[Le]):
+                tmp_exprs[i] += coeff * expr
+        exprs.extend(tmp_exprs)
+        all_lmns += lmns
+    return exprs, all_lmns