Day 1 Feb-19-2024

cell2mol/c2m_driver.py

@@ -87,6 +87,7 @@
     ######################
     print(f"running cell2mol with debug={debug}")
 
+## def __init__(self, refcode: str, labels: list, pos: list, cellvec: list, cellparam: list, warning_list: list) -> None:
     cell = cell2mol(infopath, refcode, output_dir, step, debug)
     print_types = "gmol"
     save_cell(cell, print_types, output_dir, debug=debug)

cell2mol/c2m_module.py

@@ -19,9 +19,10 @@
 from typing import Tuple
 import sklearn
 from cell2mol import __file__
-from cell2mol.elementdata import ElementData
 
+from cell2mol.other import handle_error
 
+from cell2mol.elementdata import ElementData
 elemdatabase = ElementData()
 
 ##################################################################################
@@ -472,40 +473,48 @@ def assign_spin_old (cell: object, debug: int=0) -> object:
 ##################################################################################
 ################################## MAIN ##########################################
 ##################################################################################
-def cell2mol(infopath: str, refcode: str, output_dir: str, step: int=3, debug: int=1) -> object:
-
+def cell2mol(infopath: str, refcode: str, output_dir: str, step: int=3, cellpath: str=None, debug: int=1) -> object:
+    ## Step == 1: Only Cell Reconstruction
+    ## Step == 2: Only Charge Assignment (cell reconstruction must have been run before)
+    ## Step == 3: Both Cell Reconstruction and Charge Assignment (cell reconstruction must have been run before)
+
     if step == 1 or step == 3:
 
         tini = time.time()
 
-        # Read reference molecules from info file
+        # Reads reference molecules from info file, as well as labels and coordinates
         labels, pos, ref_labels, ref_fracs, cellvec, cellparam = readinfo(infopath)
 
-        # Initialize cell object
-        warning_list = []
+        # Initiates cell
         newcell = cell(refcode, labels, pos, cellvec, cellparam, warning_list)
-        if debug >= 1: print("[Refcode]", newcell.refcode)
 
+        # Loads the reference molecules and checks_missing_H
+        newcell.get_references_molecules(ref_labels, ref_fracs, debug=debug)      ## Evaluates boolean variable self.has_isolated_H. If true, indicates a problem with the cif
+        if newcell.has_isolated_H: handle_error(1)
+        newcell.check_missing_H(debug=debug)                                      ## Evaluates boolean variable self.has_missing_H. If true, indicates a problem with the cif
+        if newcell.has_missing_H:  handle_error(2)
+
         # Cell Reconstruction
-        if debug >= 1: print("===================================== step 1 : Cell reconstruction =====================================\n")
-        newcell = reconstruct(newcell, ref_labels, ref_fracs, debug=debug)
+        newcell.reconstruct(debug=debug)
+        if newcell.is_fragmented:  handle_error(3)
+
         tend = time.time()
-        if debug >= 1: print(f"\nTotal execution time for Cell Reconstruction: {tend - tini:.2f} seconds")
+        if debug >= 1: print(f"\nCell Reconstruction Finished Normally. Total execution time: {tend - tini:.2f} seconds")
 
     elif step == 2:
+        from cell2mol.readwrite import load_binary
         if debug >= 1: print("\n***Runing only Charge Assignment***")
+        if cellpath is None: cellpath = os.path.join(output_dir, "Cell_{}.gmol".format(refcode))
+        newcell = load_binary(cellpath)
         if debug >= 1: print("\nCell object loaded with pickle")
-        cellpath = os.path.join(output_dir, "Cell_{}.gmol".format(refcode))
-        newcell = load_cell_reset_charges(cellpath)
+        newcell.reset_charges() 
+        ### newcell = load_cell_reset_charges(cellpath)       --- Split into load_binary and newcell.reset_charges
+
     else:
         if debug >= 1: print("Step number is incorrect. Only values 1, 2 or 3 are accepted")
         sys.exit(1)
 
-    if not any(newcell.warning_after_reconstruction):
-        if step == 1 or step == 3:
-            if debug >= 1: print("Cell reconstruction successfully finished.\n")
-        elif step == 2:
-            if debug >= 1: print("No Warnings in loaded Cell object in cell reconstruction \n")
+    if not newcell.has_missing_H and not newcell.has_isolated_H and not newcell.is_fragmented:
 
         if step == 1:
             pass

cell2mol/cell_operations.py

@@ -0,0 +1,140 @@
+#!/usr/bin/env python
+
+import numpy as np
+
+#######################################################
+def frac2cart_fromcellvec(frac_coord, cellvec):
+    """ Convert fractional coordinates to cartesian coordinates
+    
+    Args:
+        frac_coord (list): list of fractional coordinates
+        cellvec (list): list of cell vectors
+        
+    Returns:
+        cartesian (list): list of cartesian coordinates
+    """
+
+    cartesian = []
+    for idx, frac in enumerate(frac_coord):
+        xcar = (
+            frac[0] * cellvec[0][0] + frac[1] * cellvec[1][0] + frac[2] * cellvec[2][0]
+        )
+        ycar = (
+            frac[0] * cellvec[0][1] + frac[1] * cellvec[1][1] + frac[2] * cellvec[2][1]
+        )
+        zcar = (
+            frac[0] * cellvec[0][2] + frac[1] * cellvec[1][2] + frac[2] * cellvec[2][2]
+        )
+        cartesian.append([float(xcar), float(ycar), float(zcar)])
+    return cartesian
+
+#######################################################
+def frac2cart_fromparam(frac_coord, cellparam):
+    """ Convert fractional coordinates to cartesian coordinates
+    
+    Args:
+        frac_coord (list): list of fractional coordinates
+        cellparam (list): list of cell parameters
+        
+    Returns:
+        cartesian (list): list of cartesian coordinates
+    """
+
+    a = cellparam[0]
+    b = cellparam[1]
+    c = cellparam[2]
+    alpha = np.radians(cellparam[3])
+    beta = np.radians(cellparam[4])
+    gamma = np.radians(cellparam[5])
+
+    volume = (
+        a
+        * b
+        * c
+        * np.sqrt(
+            1
+            - np.cos(alpha) ** 2
+            - np.cos(beta) ** 2
+            - np.cos(gamma) ** 2
+            + 2 * np.cos(alpha) * np.cos(beta) * np.cos(gamma)
+        )
+    )
+
+    m = np.zeros((3, 3))
+    m[0][0] = a
+    m[0][1] = b * np.cos(gamma)
+    m[0][2] = c * np.cos(beta)
+    m[1][0] = 0
+    m[1][1] = b * np.sin(gamma)
+    m[1][2] = c * ((np.cos(alpha) - np.cos(beta) * np.cos(gamma)) / np.sin(gamma))
+    m[2][0] = 0
+    m[2][1] = 0
+    m[2][2] = volume / (a * b * np.sin(gamma))
+
+    cartesian = []
+    for idx, frac in enumerate(frac_coord):
+        xcar = frac[0] * m[0][0] + frac[1] * m[0][1] + frac[2] * m[0][2]
+        ycar = frac[0] * m[1][0] + frac[1] * m[1][1] + frac[2] * m[1][2]
+        zcar = frac[0] * m[2][0] + frac[1] * m[2][1] + frac[2] * m[2][2]
+        cartesian.append([float(xcar), float(ycar), float(zcar)])
+    return cartesian
+
+############################;###########################
+def translate(vector, coords, cellvec):
+    newcoord = []
+    for idx, coord in enumerate(coords):
+        newx = (
+            coord[0]
+            + vector[0] * cellvec[0][0]
+            + vector[1] * cellvec[1][0]
+            + vector[2] * cellvec[2][0]
+        )
+        newy = (
+            coord[1]
+            + vector[0] * cellvec[0][1]
+            + vector[1] * cellvec[1][1]
+            + vector[2] * cellvec[2][1]
+        )
+        newz = (
+            coord[2]
+            + vector[0] * cellvec[0][2]
+            + vector[1] * cellvec[1][2]
+            + vector[2] * cellvec[2][2]
+        )
+        newcoord.append([float(newx), float(newy), float(newz)])
+    return newcoord
+
+#######################################################
+def cart2frac(cartCoords, cellvec):
+    """ Translate coordinates by a vector
+    Args:
+        vector (list): list of vector components
+        coords (list): list of coordinates
+        cellvec (list): list of cell vectors
+    Returns:
+        newcoord (list): list of translated coordinates
+    """
+
+    newcoord = []
+    for idx, coord in enumerate(coords):
+        newx = (
+            coord[0]
+            + vector[0] * cellvec[0][0]
+            + vector[1] * cellvec[1][0]
+            + vector[2] * cellvec[2][0]
+        )
+        newy = (
+            coord[1]
+            + vector[0] * cellvec[0][1]
+            + vector[1] * cellvec[1][1]
+            + vector[2] * cellvec[2][1]
+        )
+        newz = (
+            coord[2]
+            + vector[0] * cellvec[0][2]
+            + vector[1] * cellvec[1][2]
+            + vector[2] * cellvec[2][2]
+        )
+        newcoord.append([float(newx), float(newy), float(newz)])
+
+    return newcoord