move points_2_plane to utils;

update `Aligning band edged to reference potentials` section in tutorial to use package function rather than internal code

macrodensity/plotting.py

@@ -25,8 +25,8 @@
     gradient_magnitude
 )
 from macrodensity.io import read_gulp_potential, read_vasp_density
-from macrodensity.tools import create_plotting_mesh, points_2_plane
-from macrodensity.utils import matrix_2_abc, vector_2_abscissa, numbers_2_grid
+from macrodensity.tools import create_plotting_mesh
+from macrodensity.utils import matrix_2_abc, vector_2_abscissa, numbers_2_grid, points_2_plane
 
 MODULE_DIR = os.path.dirname(os.path.abspath(__file__))
 plt.style.use(f"{MODULE_DIR}/macrodensity.mplstyle")

macrodensity/tools.py

@@ -51,12 +51,9 @@ def bulk_interstitial_alignment(
     """ 
 
     ## GETTING POTENTIAL
-    vasp_pot, NGX, NGY, NGZ, lattice = read_vasp_density(locpot,quiet=True)
-    vector_a,vector_b,vector_c,av,bv,cv = matrix_2_abc(lattice)
-    resolution_x = vector_a/NGX
-    resolution_y = vector_b/NGY
-    resolution_z = vector_c/NGZ
-    grid_pot, electrons = density_2_grid(vasp_pot,NGX,NGY,NGZ,Format="VASP")
+    vasp_pot, NGX, NGY, NGZ, lattice = read_vasp_density(locpot, quiet=True)
+    vector_a, vector_b, vector_c, av, bv, cv = matrix_2_abc(lattice)
+    grid_pot, electrons = density_2_grid(vasp_pot, NGX, NGY, NGZ, Format="VASP")
 
     ## GETTING BAND EDGES
     band_extrema = get_band_extrema(outcar)
@@ -67,31 +64,33 @@ def bulk_interstitial_alignment(
     interstitial_potentials = []
     interstitial_variances = []
     for interstice in interstices:
-        locpot_extract = volume_average(origin=interstice,cube=cube_size,grid=grid_pot,nx=NGX,ny=NGY,nz=NGZ)
+        locpot_extract = volume_average(
+            origin=interstice, cube=cube_size, grid=grid_pot, nx=NGX, ny=NGY, nz=NGZ
+        )
         interstitial_potentials.append(locpot_extract[0])
         interstitial_variances.append(locpot_extract[1])
 
     ## CALCULATING ALIGNED BAND ENERGIES
     sum_interstitial_potential = 0
     for ele in interstitial_potentials:
         sum_interstitial_potential += ele
-    average_interstitial_potential = sum_interstitial_potential/len(interstitial_potentials)
-    VB_aligned = round(VB_eigenvalue - average_interstitial_potential,2)
-    CB_aligned = round(CB_eigenvalue - average_interstitial_potential,2)
+    average_interstitial_potential = sum_interstitial_potential / len(interstitial_potentials)
+    VB_aligned = round(VB_eigenvalue - average_interstitial_potential, 2)
+    CB_aligned = round(CB_eigenvalue - average_interstitial_potential, 2)
 
     ## PRINTING
-    if print_output == True:
+    if print_output:
         print("Reading band edges from file: "+str(outcar))
         print("Reading potential from file: "+str(locpot))
         print("Interstital variances: "+str(interstitial_variances))
         print("VB_aligned      CB_aligned")
         print("--------------------------------")
-        print(VB_aligned,"         ",CB_aligned)
+        print(VB_aligned, "         ", CB_aligned)
 
-    return VB_aligned, CB_aligned, interstitial_variances
+    return (VB_aligned, CB_aligned, interstitial_variances)
 
 
-def subs_potentials(A: np.ndarray,B: np.ndarray,tol: float) -> np.ndarray:
+def subs_potentials(A: np.ndarray, B: np.ndarray, tol: float) -> np.ndarray:
     """
     Subtract potentials between two datasets based on a tolerance value.
 
@@ -201,7 +200,7 @@ def match_resolution(A: np.ndarray, B: np.ndarray) -> tuple:
     return A_new, B_new
 
 
-def spline_generate(A: np.ndarray,new_res_factor: int) -> np.ndarray:
+def spline_generate(A: np.ndarray, new_res_factor: int) -> np.ndarray:
     """
     Generate a new dataset with higher resolution using cubic spline interpolation.
 
@@ -232,7 +231,7 @@ def spline_generate(A: np.ndarray,new_res_factor: int) -> np.ndarray:
     return B
 
 
-def matched_spline_generate(A: np.ndarray,B: np.ndarray, V_A: np.ndarray, V_B: np.ndarray) -> tuple:
+def matched_spline_generate(A: np.ndarray, B: np.ndarray, V_A: np.ndarray, V_B: np.ndarray) -> tuple:
     """
     Generate matched datasets with the same resolution using cubic spline interpolation.
 
@@ -293,7 +292,7 @@ def matched_spline_generate(A: np.ndarray,B: np.ndarray, V_A: np.ndarray, V_B: n
     return TD_A, TD_B
 
 
-def scissors_shift(potential: np.ndarray,delta: float) -> np.ndarray:
+def scissors_shift(potential: np.ndarray, delta: float) -> np.ndarray:
     """
     Shift the potential values by a constant amount.
 
@@ -320,7 +319,7 @@ def scissors_shift(potential: np.ndarray,delta: float) -> np.ndarray:
     return shifted_potential
 
 
-def extend_potential(potential: np.ndarray,extension: float,vector: list) -> np.ndarray:
+def extend_potential(potential: np.ndarray, extension: float, vector: list) -> np.ndarray:
     """
     Extend a dataset by duplicating potential values along a specified vector direction.
 
@@ -393,7 +392,7 @@ def sort_potential(potential: np.ndarray) -> np.ndarray:
     return sorted_potential
 
 
-def diff_potentials(potential_a: np.ndarray, potential_b: np.ndarray,start: float,end: float,tol: float=0.04) -> np.ndarray:
+def diff_potentials(potential_a: np.ndarray, potential_b: np.ndarray,start: float, end: float, tol: float=0.04) -> np.ndarray:
     """
     Subtract potential values between two datasets within a specified range.
 
@@ -482,7 +481,7 @@ def translate_grid(potential: np.ndarray, translation: float, periodic: bool=Fal
     return sorted_potential_trans
 
 
-def create_plotting_mesh(NGX: int,NGY: int,NGZ: int,pc: np.ndarray,grad: np.ndarray) -> np.ndarray:
+def create_plotting_mesh(NGX: int, NGY: int, NGZ: int, pc: np.ndarray, grad: np.ndarray) -> np.ndarray:
     """
     Creates a plotting mesh based on the given grid data and plane coefficients.
 
@@ -532,42 +531,7 @@ def create_plotting_mesh(NGX: int,NGY: int,NGZ: int,pc: np.ndarray,grad: np.ndar
     return plane
 
 
-def points_2_plane(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> np.ndarray:
-    """
-    Calculates the plane coefficients from three points in space.
-
-    Parameters:
-        a (numpy.ndarray): First point with shape (3,).
-
-        b (numpy.ndarray): Second point with shape (3,).
-
-        c (numpy.ndarray): Third point with shape (3,).
-
-    Returns:
-        numpy.ndarray: An array containing the plane coefficients with shape (4,).
-
-    Example:
-        >>> # Sample points in space
-        >>> a = np.array([0, 0, 0])
-        >>> b = np.array([1, 0, 0])
-        >>> c = np.array([0, 1, 0])
-        >>> # Calculate plane coefficients
-        >>> plane_coefficients = points_2_plane(a, b, c)
-        >>> print(plane_coefficients)
-    """
-    coefficients = np.zeros(shape=(4))
-
-    ca = c - a
-    ba = b - a
-    normal = np.cross(ba,ca)
-    d = -np.dot(normal,a)
-    A, B, C = normal[0], normal[1], normal[2]
-    D = d
-    coefficients = np.array([A, B, C, D])
-    return coefficients
-
-
-def get_third_coordinate(plane_coeff: np.ndarray,NGX: int,NGY: int) -> list:
+def get_third_coordinate(plane_coeff: np.ndarray, NGX: int, NGY: int) -> list:
     """
     Computes the third coordinate of the plane for given plane coefficients.
 

macrodensity/utils.py

@@ -186,4 +186,39 @@ def numbers_2_grid(a: tuple, NGX: int, NGY: int, NGZ: int) -> np.ndarray:
     a_grid[1] = round(float(a[1])*NGY)
     a_grid[2] = round(float(a[2])*NGZ)
 
-    return a_grid
+    return a_grid
+
+
+def points_2_plane(a: np.ndarray, b: np.ndarray, c: np.ndarray) -> np.ndarray:
+    """
+    Calculates the plane coefficients from three points in space.
+
+    Parameters:
+        a (numpy.ndarray): First point with shape (3,).
+
+        b (numpy.ndarray): Second point with shape (3,).
+
+        c (numpy.ndarray): Third point with shape (3,).
+
+    Returns:
+        numpy.ndarray: An array containing the plane coefficients with shape (4,).
+
+    Example:
+        >>> # Sample points in space
+        >>> a = np.array([0, 0, 0])
+        >>> b = np.array([1, 0, 0])
+        >>> c = np.array([0, 1, 0])
+        >>> # Calculate plane coefficients
+        >>> plane_coefficients = points_2_plane(a, b, c)
+        >>> print(plane_coefficients)
+    """
+    coefficients = np.zeros(shape=(4))
+
+    ca = c - a
+    ba = b - a
+    normal = np.cross(ba,ca)
+    d = -np.dot(normal,a)
+    A, B, C = normal[0], normal[1], normal[2]
+    D = d
+    coefficients = np.array([A, B, C, D])
+    return coefficients

tests/unit_tests.py

@@ -243,8 +243,13 @@ def test_bulk_interstitial_alignment(self):
                     __name__, path_join('../tests', 'LOCPOT.test'))
         Outcar = pkg_resources.resource_filename(
                     __name__, path_join('../tests', 'OUTCAR.test'))
-        out = md.bulk_interstitial_alignment(interstices=([0.5,0.5,0.5],[0.25,0.25,0.25]),outcar=Outcar,locpot=Locpot,cube_size=[2,2,2])
-        self.assertEqual(out,(-3.24, -1.72, [1.8665165271901357e-05, 6.277207757909537e-06]))
+        out = md.bulk_interstitial_alignment(
+            interstices=([0.5,0.5,0.5],[0.25,0.25,0.25]),
+            outcar=Outcar,
+            locpot=Locpot,
+            cube_size=[2,2,2]
+        )
+        self.assertEqual(out, (-3.24, -1.72, [1.8665165271901357e-05, 6.277207757909537e-06]))
 
     def test_moving_cube(self):
         '''Tests the moving_cube function'''

tutorials/MD_workbook.ipynb

@@ -458,122 +458,52 @@
    "source": [
     "## Aligning band edged to reference potentials\n",
     "\n",
-    "To align the band edges of a material to the reference potential within its interstitial regions, we can use the function ``.\n",
-    " \n",
-    "Input parameters include the positions of interstitial spaces, the `VASP OUTCAR` and `LOCPOT` filenames, and the size of a cube defined by `LOCPOT` FFT mesh points. The function's output consists of the aligned valence band, aligned conduction band, and variances within the interstitial regions."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Input Section"
+    "To align the band edges of a material to the reference potential within its interstitial regions, we can use the function `macrodensity.tools.bulk_interstitial_alignment`. For this, we need to specify the positions of interstitial spaces, the `VASP OUTCAR` and `LOCPOT` filenames, and the size of a cube defined by `LOCPOT` FFT mesh points. The function's output consists of the aligned valence band, aligned conduction band, and variances within the interstitial regions."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "interstices = ([0.5,0.5,0.5], [0.25,0.25,0.25])\n",
-    "outcar = 'OUTCAR.test'\n",
-    "locpot = 'LOCPOT.test'\n",
-    "cube_size = [2,2,2]"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Getting potential values and band edges "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
+      "Band edges before alignment: 2.8952 eV; 4.411 eV\n",
+      "\n",
+      "Aligning the band edges...\n",
       "Reading header information...\n",
-      "Reading 3D data using Pandas...\n"
+      "Reading 3D data using Pandas...\n",
+      "\n",
+      "Band edges after alignment: -3.24 eV; -1.72 eV\n"
      ]
     }
    ],
    "source": [
-    "# GETTING POTENTIAL\n",
-    "vasp_pot, NGX, NGY, NGZ, Lattice = md.read_vasp_density(locpot, quiet=True)\n",
-    "vector_a, vector_b, vector_c, av, bv, cv = md.matrix_2_abc(Lattice)\n",
-    "resolution_x = vector_a/NGX\n",
-    "resolution_y = vector_b/NGY\n",
-    "resolution_z = vector_c/NGZ\n",
-    "grid_pot, electrons = md.density_2_grid(vasp_pot, NGX, NGY, NGZ)\n",
+    "import macrodensity\n",
     "\n",
-    "# GETTING BAND EDGES\n",
-    "band_extrema = md.get_band_extrema(outcar)\n",
-    "VB_eigenvalue = band_extrema[0]\n",
-    "CB_eigenvalue = band_extrema[1]"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Calculating reference states and the aligned band energies "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 16,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Reading band edges from file: OUTCAR.test\n",
-      "Reading potential from file: LOCPOT.test\n",
-      "Interstital variances: [1.8665165271901357e-05, 6.277207757909537e-06]\n",
-      "VB_aligned (eV)      CB_aligned (eV)\n",
-      "--------------------------------\n",
-      "-3.24               -1.72\n"
-     ]
-    }
-   ],
-   "source": [
-    "## CALCULATING REFERENCE STATE\n",
-    "interstitial_potentials = []\n",
-    "interstitial_variances = []\n",
-    "for interstice in interstices:\n",
-    "    locpot_extract = md.volume_average(\n",
-    "        origin=interstice,\n",
-    "        cube=cube_size,\n",
-    "        grid=grid_pot,\n",
-    "        nx=NGX,\n",
-    "        ny=NGY,\n",
-    "        nz=NGZ\n",
-    "    )\n",
-    "    interstitial_potentials.append(locpot_extract[0])\n",
-    "    interstitial_variances.append(locpot_extract[1])\n",
+    "interstices = ([0.5,0.5,0.5], [0.25,0.25,0.25])\n",
+    "outcar = 'OUTCAR.test'\n",
+    "locpot = 'LOCPOT.test'\n",
+    "cube_size = [2,2,2]\n",
     "\n",
-    "## CALCULATING ALIGNED BAND ENERGIES\n",
-    "sum_interstitial_potential = 0\n",
-    "for ele in interstitial_potentials:\n",
-    "    sum_interstitial_potential += ele\n",
-    "average_interstitial_potential = sum_interstitial_potential/len(interstitial_potentials)\n",
-    "VB_aligned = round(VB_eigenvalue - average_interstitial_potential,2)\n",
-    "CB_aligned = round(CB_eigenvalue - average_interstitial_potential,2)\n",
+    "# Before alignment\n",
     "\n",
-    "## PRINTING\n",
-    "print(\"Reading band edges from file: \" + str(outcar))\n",
-    "print(\"Reading potential from file: \" + str(locpot))\n",
-    "print(\"Interstital variances: \" + str(interstitial_variances))\n",
-    "print(\"VB_aligned (eV)      CB_aligned (eV)\")\n",
-    "print(\"--------------------------------\")\n",
-    "print(VB_aligned, \"             \", CB_aligned)\n"
+    "# GETTING BAND EDGES\n",
+    "VB_eigenvalue, CB_eigenvalue = macrodensity.io.get_band_extrema(outcar)\n",
+    "print(f\"Band edges before alignment: {VB_eigenvalue} eV; {CB_eigenvalue} eV\")\n",
+    "\n",
+    "# Aligning the band edges\n",
+    "print(\"\\nAligning the band edges...\")\n",
+    "VB_eigenvalue, CB_eigenvalue, interstitial_variance = macrodensity.tools.bulk_interstitial_alignment(\n",
+    "    interstices=interstices,\n",
+    "    outcar=outcar,\n",
+    "    locpot=locpot,\n",
+    "    cube_size=cube_size,\n",
+    "    print_output=False,\n",
+    ")\n",
+    "print(f\"\\nBand edges after alignment: {VB_eigenvalue} eV; {CB_eigenvalue} eV\")"
    ]
   }
  ],