Merge branch 'master' into ParallelNCFlex

Bring up to date with master

README.rst

@@ -12,7 +12,7 @@ useful routines for:
 - Elastic properties
 
 In addition to domain-specific routines, it also implements a set of
-general-purpose, low-level utilies:
+general-purpose, low-level utilities:
 
 - Efficient neighbour lists
 - Atomic strain

docs/api.rst

@@ -41,7 +41,7 @@ Analysis tools
    matscipy.neighbours
    matscipy.numerical
    matscipy.rings
-   matscipy.spatial_correlation_functions
+   matscipy.spatial_correlation_function
 
 
 Interatomic potentials
@@ -52,12 +52,15 @@ Interatomic potentials
    :recursive:
 
    matscipy.calculators.calculator
-   matscipy.calculators.pair_potential.PairPotential
-   matscipy.calculators.polydisperse.Polydisperse
-   matscipy.calculators.mcfm.MultiClusterForceMixingPotential
-   matscipy.calculators.eam.EAM
-   matscipy.calculators.ewald.Ewald
-   matscipy.calculators.manybody.Manybody
+   matscipy.calculators.committee
+   matscipy.calculators.eam
+   matscipy.calculators.ewald
+   matscipy.calculators.fitting
+   matscipy.calculators.manybody
+   matscipy.calculators.mcfm
+   matscipy.calculators.pair_potential
+   matscipy.calculators.polydisperse
+   matscipy.calculators.supercell_calculator
 
 
 Utility functions

matscipy/calculators/committee/committee.py

@@ -48,10 +48,14 @@ def __init__(self, committee, atoms=None):
     def calculate(self, atoms=None, properties=['energy', 'forces', 'stress'], system_changes=all_changes):
         """Calculates committee (mean) values and variances."""
 
+        # ACEHAL compatibility: Ignore ACEHAL explicit requests for committee_energy etc
+        properties = [p for p in properties if "committee" not in p]
+
         logger.info(f'Calculating properties {properties} with committee')
         super().calculate(atoms, properties, system_changes)
 
         property_committee = {k_i: [] for k_i in properties}
+        self.results_extra = {}
 
         for cm_i in self.committee.members:
             cm_i.calculator.calculate(atoms=atoms, properties=properties, system_changes=system_changes)
@@ -63,6 +67,12 @@ def calculate(self, atoms=None, properties=['energy', 'forces', 'stress'], syste
             self.results[p_i] = np.mean(property_committee[p_i], axis=0)
             self.results[f'{p_i}_uncertainty'] = np.sqrt(np.var(property_committee[p_i], ddof=1, axis=0))
 
+            # Compatibility with ACEHAL Bias Calculator
+            # https://github.com/ACEsuit/ACEHAL/blob/main/ACEHAL/bias_calc.py
+            self.results_extra[f'{p_i}_committee'] = property_committee[p_i]
+            self.results_extra[f'err_{p_i}'] = self.results[f'{p_i}_uncertainty']
+            self.results_extra[f'err_{p_i}_MAE'] = np.average([np.abs(comm - self.results[p_i]) for comm in property_committee[p_i]])
+
             if self.committee.is_calibrated_for(p_i):
                 self.results[f'{p_i}_uncertainty'] = self.committee.scale_uncertainty(self.results[f'{p_i}_uncertainty'], p_i)
             else:

matscipy/calculators/fitting.py

@@ -40,20 +40,21 @@
 from ase.units import GPa,J,m
 
 import scipy
-scipy_v = scipy.__version__
-
-if int(scipy_v.split('.')[1]) <= 14 : 
-    from scipy.optimize import minimize, leastsq, anneal, brute
-else :
-    # scipy.optimize.anneal decprecated from version 0.14.0, documentation advise to use scipy.optimize.basinhopping instead
-    from scipy.optimize import minimize, leastsq, brute
-    from scipy.signal import argrelextrema
-
+
+from scipy.optimize import minimize, leastsq, brute
+from scipy.signal import argrelextrema
+
+try:
+    from scipy.optimize import anneal
+except ImportError:
+    # FIXME! scipy.optimize.anneal decprecated from version 0.14.0, documentation advise to use scipy.optimize.basinhopping instead
+    pass
+
 
 try:
     from openopt import GLP
     have_openopt = True
-except:
+except ImportError:
     have_openopt = False
 
 ###

matscipy/elasticity.py

@@ -725,7 +725,7 @@ def generate_strained_configs(at0, symmetry='triclinic', N_steps=5, delta=1e-2):
 # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2, optimizer=None,
-                          verbose=True, graphics=False, logfile=None, **kwargs):
+                          verbose=True, graphics=False, logfile=None, stress_err=None, grad_scale=100.0*units.GPa, intercept_scale=5.0*units.GPa, **kwargs):
     """
     Compute elastic constants by linear regression of stress vs. strain
 
@@ -760,6 +760,12 @@ def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2, optimi
     logfile : bool
         Log file to write optimizer output to. Default None (i.e. suppress
         output).
+    stress_err : float
+        Assumed error on each stress measurement, in eV/A**2
+    grad_scale : float
+        Prior assumption of scale of elastic constants in eV/A**2 (default 100GPa)
+    intercept_scale : float
+        Prior assumption of scale of intercept in eV/A**2 (default 5GPa)
     **kwargs : dict
         Additional arguments to pass to `optimizer.run()` method e.g. `fmax`.
 
@@ -793,13 +799,54 @@ def fit_elastic_constants(a, symmetry='triclinic', N_steps=5, delta=1e-2, optimi
     if graphics:
         import matplotlib.pyplot as plt
 
-    def do_fit(index1, index2, stress, strain, patt):
+    def do_fit(index1, index2, stress, strain, patt, stress_err, grad_scale, intercept_scale):
+        def blr_linregress(x, y, y_sig, sigma_m, sigma_c):
+            '''
+            Use analytical blr to solve y = m * x + c given x, y, and assumed errors
+
+            y_sig : float | array
+                Assumed variance of each observation.
+            sigma_m : float
+                Prior variance of the gradient
+            sigma_c : float
+                Prior variance of the intercept
+            '''
+            # Design Matrix
+            X = np.array([x, np.ones(x.shape)]).T
+
+            # Prior Precision Matrix
+            Gamma = np.diag([1/sigma_m**2, 1/sigma_c**2])
+
+            # Likelihood Precision
+            if type(y_sig == float):
+                Lambda = np.eye(y.shape[0]) / y_sig**2
+            else:
+                Lambda = np.diag(1/y_sig ** 2)
+
+            # Solve for posterior
+            posterior_variance= X.T @ Lambda @ X + Gamma
+
+            # 2x2, so direct inverse not unstable 
+            posterior_precision = np.linalg.inv(posterior_variance)
+
+            grad, intercept = posterior_precision @ X.T @ Lambda @ y
+
+            stderr = posterior_precision[0, 0] # Only care about gradient variance
+
+            r = np.corrcoef(x, y)[0, 1] # Correlation Coefficient
+
+            return grad, intercept, r, stderr
+
         if verbose:
             print('Fitting C_%d%d' % (index1+1, index2+1))
             print('Strain %r' % strain[:,index2])
             print('Stress %r GPa' % (stress[:,index1]/units.GPa))
 
-        if scipy_stats is not None:
+
+        if stress_err is not None:
+            cijFitted, intercept, r, stderr = blr_linregress(strain[:,index2],stress[:,index1], stress_err, grad_scale, intercept_scale)
+            tt = None
+        elif scipy_stats is not None:
             cijFitted, intercept,r,tt,stderr = scipy_stats.linregress(strain[:,index2],stress[:,index1])
         else:
             cijFitted, intercept = np.polyfit(strain[:,index2],stress[:,index1], 1)
@@ -808,7 +855,7 @@ def do_fit(index1, index2, stress, strain, patt):
         if verbose:
             # print info about the fit
             print('Cij (gradient) / GPa    :    ',cijFitted/units.GPa)
-            if scipy_stats is not None:
+            if scipy_stats is not None or stress_err is not None:
                 print('Error in Cij / GPa      :    ', stderr/units.GPa)
                 if abs(r) > 0.9:
                     print('Correlation coefficient :    ',r)
@@ -899,7 +946,8 @@ def do_fit(index1, index2, stress, strain, patt):
             fitted, err = do_fit(index1, index2,
                                  stress[pattern_index,:,:],
                                  strain[pattern_index,:,:],
-                                 pattern_index)
+                                 pattern_index, 
+                                 stress_err, grad_scale, intercept_scale)
 
             index = abs(Cij_map_sym[(index1, index2)])
 

paper/paper.md

@@ -45,6 +45,9 @@ authors:
   - name: Thomas Reichenbach
     orcid: 0000-0001-7477-6248
     affiliation: 5
+  - name: Thomas Rocke
+    orcid: 0000-0002-4612-9112
+    affiliation: 3
   - name: Lakshmi Shenoy
     orcid: 0000-0001-5760-3345
     affiliation: 3
@@ -106,7 +109,7 @@ We point out that other generic multi-scale coupling packages exist. Examples of
 
 Within materials science, the package has different application domains:
 
-- **Elasticity.** Solids respond to small external loads through a reversible elastic response. The strength of the response is characterized by the elastic moduli. `matscipy.elasticity` implements functions for computing elastic moduli from small deformation that consider potential symmetries of the underlying atomic system, in particular for crystals. `matscipy` also implements analytic calculation of elastic moduli for some interatomic potentials, described in more detail below. The computation of elastic moduli is a prerequisite for multi-scale modelling of materials, as they are the most basic parameters of continuum material models. `matscipy` was used to study finite-pressure elastic constants and structural stability in crystals [@Griesser2023crystal] and glasses [@Griesser2023glass].
+- **Elasticity.** Solids respond to small external loads through a reversible elastic response. The strength of the response is characterized by the elastic moduli. `matscipy.elasticity` implements functions for computing elastic moduli from small deformation that consider potential symmetries of the underlying atomic system, in particular for crystals. The implementation also includes estimates of uncertainty on elastic moduli - either from a least-squares error, or from a Bayesian treatment if stress uncertainty is supplied. `matscipy` also implements analytic calculation of elastic moduli for some interatomic potentials, described in more detail below. The computation of elastic moduli is a prerequisite for multi-scale modelling of materials, as they are the most basic parameters of continuum material models. `matscipy` was used to study finite-pressure elastic constants and structural stability in crystals [@Griesser2023crystal] and glasses [@Griesser2023glass].
 
 - **Plasticity.** For large loads, solids can respond with irreversible deformation. One form of irreversibility is plasticity, that is carried by extended defects, the dislocations, in crystals. The module `matscipy.dislocation` implements tools for studying structure and movement of dislocations. Construction and analysis of model atomic systems is implemented for compact and dissociated screw, as well as edge dislocations in cubic crystals. The implementation supports ideal straight as well as kinked dislocations. Some of the dislocation functionality requires the `atomman` and/or `OVITO` packages as optional dependencies [@AtomMan;@Stukowski2009]. The module was employed in a study of interaction of impurities with screw dislocations in tungsten [@Grigorev2020;@Grigorev2023]. 
 

tests/test_fit_elastic_constants.py

@@ -203,6 +203,12 @@ def testtriclinic_relaxed(self):
                                              verbose=False, graphics=False)
             self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.2)
 
+        def testtriclinic_relaxed_assumed_err(self):
+            C, C_err = fit_elastic_constants(self.at0, 'triclinic',
+                                             optimizer=FIRE, fmax=self.fmax,
+                                             verbose=False, graphics=False, stress_err=0.05/self.at0.get_volume())
+            self.assertArrayAlmostEqual(C/units.GPa, self.C_ref_relaxed, tol=0.2)
+
 if __name__ == '__main__':
     unittest.main()