Merge pull request #52 from BoothGroup/cleanup

Maintenance bits

ebcc/brueckner.py

@@ -8,6 +8,7 @@
 from ebcc import NullLogger
 from ebcc import numpy as np
 from ebcc import util
+from ebcc.damping import DIIS
 from ebcc.precision import types
 
 
@@ -29,12 +30,15 @@ class Options:
     diis_space : int, optional
         Number of amplitudes to use in DIIS extrapolation. Default value is
         `12`.
+    damping : float, optional
+        Damping factor for DIIS extrapolation. Default value is `0.0`.
     """
 
     e_tol: float = 1e-8
     t_tol: float = 1e-8
     max_iter: int = 20
     diis_space: int = 12
+    damping: float = 0.0
 
 
 class BruecknerREBCC:
@@ -79,6 +83,7 @@ def __init__(self, cc, log=None, options=None, **kwargs):
         cc.log.info(" > t_tol:  %s", options.t_tol)
         cc.log.info(" > max_iter:  %s", options.max_iter)
         cc.log.info(" > diis_space:  %s", options.diis_space)
+        cc.log.info(" > damping:  %s", options.damping)
         cc.log.debug("")
 
     def get_rotation_matrix(self, u_tot=None, diis=None, t1=None):
@@ -276,8 +281,9 @@ def kernel(self):
                 self.cc.kernel()
 
         # Set up DIIS:
-        diis = lib.diis.DIIS()
+        diis = DIIS()
         diis.space = self.options.diis_space
+        diis.damping = self.options.damping
 
         # Initialise coefficients:
         mo_coeff_new = np.array(self.cc.mo_coeff, copy=True, dtype=types[float])
@@ -287,9 +293,17 @@ def kernel(self):
 
         self.cc.log.output("Solving for Brueckner orbitals.")
         self.cc.log.info(
-            "%4s %16s %10s %16s %16s", "Iter", "Energy (corr.)", "Converged", "Δ(Energy)", "|T1|"
+            "%4s %16s %18s %8s %13s %13s",
+            "Iter",
+            "Energy (corr.)",
+            "Energy (tot.)",
+            "Conv.",
+            "Δ(Energy)",
+            "|T1|",
+        )
+        self.log.info(
+            "%4d %16.10f %18.10f %8s", 0, self.cc.e_corr, self.cc.e_tot, self.cc.converged
         )
-        self.cc.log.info("%4d %16.10f %10s", 0, self.cc.e_corr, self.cc.converged)
 
         converged = False
         for niter in range(1, self.options.max_iter + 1):
@@ -322,8 +336,14 @@ def kernel(self):
             de = abs(e_prev - self.cc.e_tot)
             dt = self.get_t1_norm()
 
-            self.cc.log.info(
-                "%4d %16.10f %10s %16.5g %16.5g", niter, self.cc.e_corr, self.cc.converged, de, dt
+            self.log.info(
+                "%4d %16.10f %18.10f %8s %13.3e %13.3e",
+                niter,
+                self.cc.e_corr,
+                self.cc.e_tot,
+                self.cc.converged,
+                de,
+                dt,
             )
 
             # Check for convergence:

ebcc/damping.py

@@ -0,0 +1,29 @@
+"""Damping and DIIS control."""
+
+from pyscf.lib import diis
+
+
+class DIIS(diis.DIIS):
+    """Direct inversion in the iterative subspace."""
+
+    def __init__(self, space=6, min_space=1, damping=0.0):
+        super().__init__(incore=True)
+        self.verbose = 0
+        self.space = space
+        self.min_space = min_space
+        self.damping = damping
+
+    def update(self, x, xerr=None):
+        """Extrapolate a vector."""
+
+        # Extrapolate the vector
+        x = super().update(x, xerr=xerr)
+
+        # Apply damping
+        if self.damping:
+            nd = self.get_num_vec()
+            if nd > 1:
+                xprev = self.get_vec(self.get_num_vec() - 1)
+                x = (1.0 - self.damping) * x + self.damping * xprev
+
+        return x

ebcc/rebcc.py

@@ -10,6 +10,7 @@
 from ebcc.ansatz import Ansatz
 from ebcc.brueckner import BruecknerREBCC
 from ebcc.cderis import RCDERIs
+from ebcc.damping import DIIS
 from ebcc.dump import Dump
 from ebcc.eris import RERIs
 from ebcc.fock import RFock
@@ -46,13 +47,16 @@ class Options:
     diis_space : int, optional
         Number of amplitudes to use in DIIS extrapolation. Default value is
         `12`.
+    damping : float, optional
+        Damping factor for DIIS extrapolation. Default value is `0.0`.
     """
 
     shift: bool = True
     e_tol: float = 1e-8
     t_tol: float = 1e-8
     max_iter: int = 200
     diis_space: int = 12
+    damping: float = 0.0
 
 
 class REBCC(EBCC):
@@ -326,6 +330,7 @@ def __init__(
         self.log.info(" > t_tol:  %s", self.options.t_tol)
         self.log.info(" > max_iter:  %s", self.options.max_iter)
         self.log.info(" > diis_space:  %s", self.options.diis_space)
+        self.log.info(" > damping:  %s", self.options.damping)
         self.log.debug("")
         self.log.info("Ansatz: %s", self.ansatz)
         self.log.debug("")
@@ -364,13 +369,21 @@ def kernel(self, eris=None):
         e_cc = e_init = self.energy(amplitudes=amplitudes, eris=eris)
 
         self.log.output("Solving for excitation amplitudes.")
-        self.log.info("%4s %16s %16s %16s", "Iter", "Energy (corr.)", "Δ(Energy)", "Δ(Amplitudes)")
-        self.log.info("%4d %16.10f", 0, e_init)
+        self.log.info(
+            "%4s %16s %18s %13s %13s",
+            "Iter",
+            "Energy (corr.)",
+            "Energy (tot.)",
+            "Δ(Energy)",
+            "Δ(Ampl.)",
+        )
+        self.log.info("%4d %16.10f %18.10f", 0, e_init, e_init + self.mf.e_tot)
 
         if not self.ansatz.is_one_shot:
             # Set up DIIS:
-            diis = lib.diis.DIIS()
+            diis = DIIS()
             diis.space = self.options.diis_space
+            diis.damping = self.options.damping
 
             converged = False
             for niter in range(1, self.options.max_iter + 1):
@@ -388,7 +401,9 @@ def kernel(self, eris=None):
                 e_cc = self.energy(amplitudes=amplitudes, eris=eris)
                 de = abs(e_prev - e_cc)
 
-                self.log.info("%4d %16.10f %16.5g %16.5g", niter, e_cc, de, dt)
+                self.log.info(
+                    "%4d %16.10f %18.10f %13.3e %13.3e", niter, e_cc, e_cc + self.mf.e_tot, de, dt
+                )
 
                 # Check for convergence:
                 converged = de < self.options.e_tol and dt < self.options.t_tol
@@ -463,11 +478,12 @@ def solve_lambda(self, amplitudes=None, eris=None):
             lambdas = self.lambdas
 
         # Set up DIIS:
-        diis = lib.diis.DIIS()
+        diis = DIIS()
         diis.space = self.options.diis_space
+        diis.damping = self.options.damping
 
         self.log.output("Solving for de-excitation (lambda) amplitudes.")
-        self.log.info("%4s %16s", "Iter", "Δ(Amplitudes)")
+        self.log.info("%4s %13s", "Iter", "Δ(Ampl.)")
 
         converged = False
         for niter in range(1, self.options.max_iter + 1):
@@ -485,7 +501,7 @@ def solve_lambda(self, amplitudes=None, eris=None):
             lambdas = self.vector_to_lambdas(vector)
             dl = np.linalg.norm(vector - self.lambdas_to_vector(lambdas_prev), ord=np.inf)
 
-            self.log.info("%4d %16.5g", niter, dl)
+            self.log.info("%4d %13.3e", niter, dl)
 
             # Check for convergence:
             converged = dl < self.options.t_tol