Update mean-field scripts

docs/environment.yml

@@ -6,6 +6,7 @@ dependencies:
 - nbsphinx
 - sphinx-gallery
 - sphinx-rtd-theme
+- versioneer
 - ase =3.22.1
 - coveralls
 - coverage

pyiron_contrib/__init__.py

@@ -49,6 +49,7 @@ def __init__(self, *args, **kwargs):
         "GleMD": "pyiron_contrib.atomistics.ipi.ipi_jobs",
         "PigletMD": "pyiron_contrib.atomistics.ipi.ipi_jobs",
         "QuasiHarmonicApproximation": "pyiron_contrib.atomistics.atomistics.master.qha",
+        "TemperatureRemapping": "pyiron_contrib.atomistics.tra.temperature_remapping",
     }
 )
 

pyiron_contrib/atomistics/mean_field/core/bond_analysis.py

@@ -154,9 +154,9 @@ def _get_nn(structure, cutoff_radius):
         neighbors = structure.get_neighbors(
             num_neighbors=None, cutoff_radius=cutoff_radius
         )
-        nn_distances = np.around(neighbors.distances, decimals=5)
+        nn_distances = np.around(neighbors.distances, decimals=8)
         _, _, n_bonds_per_shell = np.unique(
-            np.around(nn_distances[0], decimals=5),
+            np.around(nn_distances[0], decimals=8),
             return_inverse=True,
             return_counts=True,
         )
@@ -174,7 +174,7 @@ def _set_n_shells(nn_distances, n_shells=None):
             n_shells
         """
         if n_shells is None:
-            nn_dists = np.around(nn_distances[0], decimals=5)
+            nn_dists = np.around(nn_distances[0], decimals=8)
             n_shells = int(np.unique(nn_dists, return_index=True)[1][-1] + 1)
         return n_shells
 
@@ -198,46 +198,22 @@ def _populate_shells(n_bonds_per_shell, vectors, indices):
             sums += n
         return vectors_per_shell
 
-    @staticmethod
-    def _rotation_matrix_from_vectors(vec_1, vec_2):
-        """
-        Find the rotation matrix that aligns vec_1 to vec_2.
-        Args:
-            vec_1: A 3d "source" vector
-            vec_2: A 3d "destination" vector
-
-        Returns:
-            A transformation matrix (3x3) which when applied to vec1, aligns it with vec2.
-        """
-        a, b = (vec_1 / np.linalg.norm(vec_1)).reshape(3), (
-            vec_2 / np.linalg.norm(vec_2)
-        ).reshape(3)
-        v = np.cross(a, b)
-        c = np.dot(a, b)
-        if np.any(v):  # if not all zeros then
-            s = np.linalg.norm(v)
-            kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
-            return np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s**2))
-        elif not np.any(v) and c < 0.0:
-            return np.eye(3) * -1  # for opposite directions
-        else:
-            return np.eye(3)  # for identical directions
-
     def _get_irreducible_bond_vector_per_shell(self, nn, n_bonds_per_shell):
         """
         Get one irreducible bond vector per nearest neighbor shell. If symmetries are known, this irreducible bond can
             be used to generate all other bonds for the corresponding shell.
         """
         # arrange bond vectors according to their shells
         bond_vectors_per_shell = self._populate_shells(
-            vectors=np.around(nn.vecs[0], decimals=5),
+            vectors=np.around(nn.vecs[0], decimals=8),
             indices=nn.indices[0],
             n_bonds_per_shell=n_bonds_per_shell,
         )
         # save 1 irreducible bond vector per shell
         per_shell_irreducible_bond_vectors = np.array(
             [b[0] for b in bond_vectors_per_shell]
         )
+
         # get all the rotations from spglib
         data = self.input.structure.get_symmetry()
         # account for only 0 translation rotations
@@ -247,23 +223,21 @@ def _get_irreducible_bond_vector_per_shell(self, nn, n_bonds_per_shell):
         # and sort them doing the following:
         per_shell_0K_rotations = []
         all_nn_bond_vectors_list = []
-        for b in enumerate(per_shell_irreducible_bond_vectors):
+        for b in per_shell_irreducible_bond_vectors:
             # get 'unique bonds' from the spglib data
-            unique_bonds = np.unique(np.dot(b[1], all_rotations), axis=0)
+            unique_bonds = np.unique(np.dot(b, all_rotations), axis=0)
             args_list = []
             for bond in unique_bonds:
                 # collect the arguments of the rotations that that give the unique bonds
                 args = []
-                for i, bd in enumerate(np.dot(b[1], all_rotations)):
-                    if np.array_equal(
-                        np.round(bd, decimals=3), np.round(bond, decimals=3)
-                    ):
+                for i, bd in enumerate(np.dot(b, all_rotations)):
+                    if np.array_equal(np.round(bd, decimals=8), np.round(bond, decimals=8)):
                         args.append(i)
                 args_list.append(args[0])
             # sort the arguments and append the rotations to a list...
             per_shell_0K_rotations.append(all_rotations[np.sort(args_list)])
             # and the unique bonds also, to another list...
-            all_nn_bond_vectors_list.append(np.dot(b[1], per_shell_0K_rotations[-1]))
+            all_nn_bond_vectors_list.append(np.dot(b, per_shell_0K_rotations[-1]))
         # and clump the unique bonds into a single array
         all_nn_bond_vectors_list = np.array(
             [
@@ -306,18 +280,13 @@ def _get_transformations(self, n_bonds_per_shell, all_nn_bond_vectors):
                 # second bond is normal to the first (transverse1). If multiple normal bonds, select 1
                 try:
                     b2 = bonds[
-                        np.argwhere(np.round(bonds @ b1, decimals=5) == 0.0).flatten()[
+                        np.argwhere(np.round(bonds@b1, decimals=8) == 0.0).flatten()[
                             0
                         ]
                     ]
                 except IndexError:
-                    # if in case a normal bond is not found for the shell, traverse through the previous shells as well
-                    b2 = all_bonds[
-                        np.argwhere(
-                            np.round(all_bonds @ b1, decimals=5) == 0.0
-                        ).flatten()[0]
-                    ]
-                # third bond is then  normal to both the first and second bonds (transverse2)
+                    # if in case a normal bond is not found for the shell, traverse through the other input shells as well
+                    b2 = all_bonds[np.argwhere(np.round(all_bonds@b1, decimals=8) == 0.).flatten()[0]]
                 b3 = np.cross(b1, b2)
                 if b1.dot(np.cross(b2, b3)) < 0.0:  # if this condition is not met
                     b2, b3 = b3, b2  # reverse b2 and b3
@@ -348,7 +317,7 @@ def _get_bond_indexed_neighbor_list(
         Returns:
             per_shell_bond_indexed_neighbor_list
         """
-        nn_vecs = np.around(nn.vecs, decimals=5)
+        nn_vecs = np.around(nn.vecs, decimals=8)
         nn_indices = nn.indices
         sums = 0
         per_shell_bond_indexed_neighbor_list = []
@@ -358,7 +327,7 @@ def _get_bond_indexed_neighbor_list(
             bond_vectors_list = all_nn_bond_vectors[sums : sums + n].copy()
             sums += n
             # initialize the M_ij matrix with shape [atoms x unique_bonds]
-            x_0 = np.around(structure.positions, decimals=5)  # zero Kelvin positions
+            x_0 = np.around(structure.positions, decimals=8)  # zero Kelvin positions
             M_ij = np.zeros([len(x_0), len(nn_vecs_per_shell[0])]).astype(int)
             # populate the M_ij matrix
             for i, per_atom_nn in enumerate(nn_vecs_per_shell):
@@ -742,47 +711,30 @@ def _get_corr_column_value(axis):
         else:
             raise ValueError("choose between long, t1, and t2")
 
-    def _get_correlations(self, shell, bond_x, bond_y, axis_x, axis_y, n_bins):
-        per_shell_r_t1_t2_bond_vectors = self.get_r_t1_t2_bond_vectors()
-        all_bonds = per_shell_r_t1_t2_bond_vectors[shell]
-        n_bonds = len(all_bonds)
-        if (bond_x >= n_bonds) or (bond_y >= n_bonds):
-            raise ValueError(
-                "there are only {} bonds in shell {}".format(n_bonds, shell)
-            )
+    def _get_correlations(self, shell_x, shell_y, bond_x, bond_y, axis_x, axis_y, n_bins):
+        per_shell_l_t1_t2_bond_vectors = self.get_l_t1_t2_bond_vectors()
+        all_bonds_x = per_shell_l_t1_t2_bond_vectors[shell_x]
+        all_bonds_y = per_shell_l_t1_t2_bond_vectors[shell_y]
+        if (bond_x >= len(all_bonds_x)):
+            raise ValueError("there are only {} bonds in shell {}".format(len(all_bonds_x), shell_x))
+        elif (bond_y >= len(all_bonds_y)):
+            raise ValueError("there are only {} bonds in shell {}".format(len(all_bonds_y), shell_y))
         axis_0 = self._get_corr_column_value(axis_x)
         axis_1 = self._get_corr_column_value(axis_y)
-        return self._get_rho_corr_and_uncorr(
-            x_data=all_bonds[bond_y, :, axis_0],
-            y_data=all_bonds[bond_x, :, axis_1],
-            n_bins=n_bins,
-        )
+        return self._get_rho_corr_and_uncorr(x_data=all_bonds_x[bond_x, :, axis_0], y_data=all_bonds_y[bond_y, :, axis_1],
+                                             n_bins=n_bins)
 
-    def get_mutual_information(
-        self, shell=0, bond_x=0, bond_y=0, axis_x="long", axis_y="long", n_bins=101
-    ):
-        rho_corr, rho_uncorr, _, _ = self._get_correlations(
-            shell=shell,
-            bond_x=bond_x,
-            bond_y=bond_y,
-            axis_x=axis_x,
-            axis_y=axis_y,
-            n_bins=n_bins,
-        )
-        sel = (rho_uncorr > 0.0) * (rho_corr > 0.0)
-        return np.sum(rho_corr[sel] * np.log(rho_corr[sel] / rho_uncorr[sel]))
+    def get_mutual_information(self, shell_x=0, shell_y=0, bond_x=0, bond_y=0, axis_x='long', axis_y='long', n_bins=101):
+        rho_corr, rho_uncorr, _, _ = self._get_correlations(shell_x=shell_x, shell_y=shell_y, bond_x=bond_x, bond_y=bond_y, axis_x=axis_x,
+                                                            axis_y=axis_y, n_bins=n_bins)
+        sel = (rho_uncorr>0.0)*(rho_corr>0.0)
+        return np.sum(rho_corr[sel]*np.log(rho_corr[sel]/rho_uncorr[sel]))
+
+    def plot_correlations(self, shell_x=0, shell_y=0, bond_x=0, bond_y=0, axis_x='long', axis_y='long', n_bins=101):
+        rho_corr, rho_uncorr, x_edges, y_edges = self._get_correlations(shell_x=shell_x, shell_y=shell_y,
+                                                                        bond_x=bond_x, bond_y=bond_y, 
+                                                                        axis_x=axis_x, axis_y=axis_y, n_bins=n_bins)
 
-    def plot_correlations(
-        self, shell=0, bond_x=0, bond_y=0, axis_x="long", axis_y="long", n_bins=101
-    ):
-        rho_corr, rho_uncorr, x_edges, y_edges = self._get_correlations(
-            shell=shell,
-            bond_x=bond_x,
-            bond_y=bond_y,
-            axis_x=axis_x,
-            axis_y=axis_y,
-            n_bins=n_bins,
-        )
         rho_diff = rho_corr - rho_uncorr
         cm = LinearSegmentedColormap.from_list(
             "my_spec", [myc["b"], (1, 1, 1), myc["o"]], N=n_bins
@@ -796,8 +748,8 @@ def plot_correlations(
         plt.imshow(rho_diff[::-1, :], cmap=cm, extent=plims)
         plt.xticks(np.around(np.linspace(plims[0], plims[1], 5), decimals=2))
         plt.yticks(np.around(np.linspace(plims[2], plims[3], 5), decimals=2))
-        plt.xlabel(axis_x + "_" + str(bond_x))
-        plt.ylabel(axis_y + "_" + str(bond_y))
+        plt.xlabel('shell_' + str(shell_x) + '_' + axis_x + '_' + str(bond_x))
+        plt.ylabel('shell_' + str(shell_y) + '_' + axis_y + '_' + str(bond_y))
         plt.show()
 
 

pyiron_contrib/atomistics/mean_field/core/hessian.py

@@ -0,0 +1,165 @@
+# coding: utf-8
+# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
+# Distributed under the terms of "New BSD License", see the LICENSE file.
+
+import os
+import numpy as np
+from phonopy.structure.grid_points import get_qpoints
+
+class GenerateHessian():
+    """
+    Class to generate the Hessian blocks for a reference atom in a bulk fcc or bcc crystal.
+    Used to generate covariance matrices alpha_{rpl} for the mean-field model.
+
+    Args:
+        project (str): Name of the pyiron project for storing the jobs.
+        ref_job (pyiron Job): Reference job containing an appropriate structure and an interatomic potential.
+        potential (potential function): A valid interatomic potential.
+        ref_atom (int): The id of the atom which will be displaced. Defaults to 0.
+        delta_x (float): Displacement of the reference atom from its equilibrium position. Defaults to 0.01 Angstroms.
+        qpoints (float): Qpoints along each direction (the mesh will be qpoints x qpoints x qpoints). Defaults to 10.
+        cutoff_radius (float): The cutoff radius upto which the nearest neighbors are considered. Defaults to None, consider all atoms as neighbors.
+    """
+    def __init__(self, project, ref_job, potential, structure=None, ref_atom=0, delta_x=0.01, qpoints=10, cutoff_radius=None):
+        self.project = project
+        self.ref_job = ref_job
+        self.potential = potential
+        self.structure = structure
+        self.ref_atom = ref_atom
+        self.delta_x = delta_x
+        self.qpoints = qpoints
+        self.cutoff_radius = cutoff_radius
+
+        self._jobs = None
+        if self.structure is None:
+            self.structure = self.ref_job.structure.copy()
+
+    def _run_job(self, job_name, i, j):
+        """
+        Helper class to run the displacement job.
+        """
+        job = self.ref_job.copy_template(project=self.project, new_job_name=job_name)
+        if self.potential is not None:
+            job.potential = self.potential
+        dx = np.zeros_like(job.structure.positions)
+        if i%2 == 0:
+            dx[self.ref_atom][j] = self.delta_x
+        else:
+            dx[self.ref_atom][j] = -self.delta_x
+        job.structure.positions += dx
+        job.calc_static()
+        job.run()
+
+    def run_jobs(self, delete_existing_jobs=False):
+        """
+        Run the displacement jobs.
+
+        Parameters:
+            delete_existing_jobs (boolean): Delete the existing displacement jobs.
+        """
+        job_list = self.project.job_table().job.to_list()
+        for i in range(6):
+            # make displacements for the ref atom along x, y, z in the + and - directions
+            # arragned as [x+, x-, y+, y-, z+, z-]
+            job_name = 'disp_job_' + str(i)
+            if job_name not in job_list:
+                    self._run_job(job_name=job_name, i=i, j=int((i-(i%2))/2))
+            else:
+                job = self.project.inspect(job_name)
+                if job.status in ['aborted'] or delete_existing_jobs:
+                    self.project.remove_job(job_name)
+                    self._run_job(job_name=job_name, i=i, j=int((i-(i%2))/2))
+
+    def load_jobs(self):
+        """
+        Load the displacement jobs.
+        """
+        self._jobs = [self.project.inspect('disp_job_' + str(i)) for i in range(6)]
+
+    def get_forces(self):
+        """
+        Get the forces on all the atoms from each displacement job.
+
+        Returns:
+            forces (np.ndarray): 6 x forces on all the atoms.
+        """
+        self.load_jobs()
+        return np.array([job['output/generic/forces'][-1] for job in self._jobs])
+
+    def get_hessian_direct(self):
+        """
+        Get the Hessian blocks for the reference atom in direct space.
+
+        Returns:
+            hessian (np.ndarray): n_atoms x 3 x 3 array of hessians.
+        """
+
+        forces = self.get_forces()
+        forces_xyz = np.array([(forces[0]-forces[1]).T, 
+                               (forces[2]-forces[3]).T, 
+                               (forces[4]-forces[5]).T])/2.
+        hessian_direct = (-forces_xyz/self.delta_x).transpose(2, 0, 1)  # shape (n_atoms, 3, 3)
+        return hessian_direct
+
+    def get_qpoint_vectors(self):
+        """
+        Get qpoint vectors and their weights (without rotational symmetry).
+
+        Returns:
+            qpoint_vectors (np.ndarray): Reduced qpoints**3 x 3 vectors.
+            qpoint_weights (np.array): Weights of the reduced qpoint vectors.
+        """
+        structure = self.structure.copy()
+        sym = structure.get_symmetry()
+        reciprocal_cell = sym.get_primitive_cell().cell.reciprocal()
+        # n = np.arange(self.qpoints**3)
+        # ix = np.array([n//self.qpoints**2, (n//self.qpoints)%self.qpoints, n%self.qpoints]).T.astype(float)
+        # qpoint_vectors = (2.0*np.pi/float(self.qpoints))*(ix+0.5)@reciprocal_cell
+        # qpoint_vectors -= qpoint_vectors.mean(0)
+        # qpoint_weights = np.ones(len(qpoint_vectors))
+        mesh = [self.qpoints, self.qpoints, self.qpoints]
+        grid, qpoint_weights = get_qpoints(mesh_numbers=mesh, reciprocal_lattice=reciprocal_cell)
+        qpoint_vectors = (2.0*np.pi*grid)@reciprocal_cell
+        return qpoint_vectors, qpoint_weights
+
+    def get_hessian_reciprocal(self, structure=None, hessian_direct=None, rewrite=False):
+        """
+        Get the Hessian blocks for the reference atom in reciprocal space, along with the qpoint vectors and their weights.
+
+        Parameters:
+            rewrite (boolean): Whether or not to rewrite the reciprocal hessians, qpoint vectors, and qpoint weights. Defaults to False.
+
+        Returns:
+            hessian_reciprocal (np.ndarray): Reduced qpoints**3 x 3 x 3 Hessians.
+            qpoint_vectors (np.ndarray): Reduced qpoints**3 x 3 vectors.
+            qpoint_weights (np.array): Weights of the reduced qpoint vectors.
+        """
+        if not os.path.exists(os.path.join(self.project.path, 'resources')):
+            os.mkdir(os.path.join(self.project.path, 'resources'))
+        exists = [os.path.exists(os.path.join(self.project.path, 'resources', file)) for file in ['hessian_reciprocal.npy', 'qpoint_vectors.npy', 'qpoint_weights.npy']]
+        if np.all(exists) and (not rewrite):
+            return (np.load(os.path.join(self.project.path, 'resources', 'hessian_reciprocal.npy')),
+                    np.load(os.path.join(self.project.path, 'resources', 'qpoint_vectors.npy')),
+                    np.load(os.path.join(self.project.path, 'resources', 'qpoint_weights.npy'))
+                    )
+        else:
+            if hessian_direct is None:
+                hessian_direct = self.get_hessian_direct()
+            qpoint_vectors, qpoint_weights = self.get_qpoint_vectors()
+            structure = self.structure.copy()
+            X = structure.positions.copy()
+            if self.cutoff_radius is not None:
+                select = structure.get_neighborhood(positions=X[self.ref_atom], cutoff_radius=self.cutoff_radius, num_neighbors=None).indices
+                if len(select)>structure.get_number_of_atoms():
+                    select = np.ones(structure.get_number_of_atoms(), dtype=bool)
+            else:
+                sq_trace = np.einsum('ijk,ikj->i', hessian_direct, hessian_direct)
+                select = sq_trace > 1e-3
+            dX = structure.find_mic(X-X[self.ref_atom])
+            q_dX = qpoint_vectors@dX[select].T
+            hessian_reciprocal = np.einsum('il,ijk->ljk', np.exp(1j*q_dX.T), hessian_direct[select])
+
+            np.save(os.path.join(self.project.path, 'resources', 'hessian_reciprocal.npy'), hessian_reciprocal)
+            np.save(os.path.join(self.project.path, 'resources', 'qpoint_vectors.npy'), qpoint_vectors)
+            np.save(os.path.join(self.project.path, 'resources', 'qpoint_weights.npy'), qpoint_weights)
+            return hessian_reciprocal, qpoint_vectors, qpoint_weights