From d81e574ff0be87bbda3d2b17ec185b0b7ad149ed Mon Sep 17 00:00:00 2001
From: qzhu2017 <alecfans@gmail.com>
Date: Fri, 4 Oct 2024 21:52:57 -0400
Subject: [PATCH] Running a systematic check with vscode

---
 pyxtal/__init__.py          |  18 ++++---
 pyxtal/database/element.py  |   2 +-
 pyxtal/db.py                |  17 +++---
 pyxtal/interface/gulp.py    |   2 +-
 pyxtal/lattice.py           |   3 +-
 pyxtal/lego/SO3.py          |  16 +++---
 pyxtal/lego/builder.py      |   5 +-
 pyxtal/molecular_crystal.py |  50 +++---------------
 pyxtal/molecule.py          |   2 +-
 pyxtal/optimize/QRS.py      |   2 +-
 pyxtal/optimize/base.py     |   4 +-
 pyxtal/optimize/common.py   |   3 +-
 pyxtal/symmetry.py          |   1 +
 pyxtal/wyckoff_site.py      | 100 ++++++++++++++++++++++++++----------
 pyxtal/wyckoff_split.py     |   2 -
 15 files changed, 120 insertions(+), 107 deletions(-)

diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py
index 58cd8678..bd51dbb9 100644
--- a/pyxtal/__init__.py
+++ b/pyxtal/__init__.py
@@ -1,5 +1,8 @@
 """
-main pyxtal module to create the pyxtal class
+PyXtal: A Python library for crystal structure generation and manipulation.
+
+This module initializes the pyxtal package, providing essential imports and
+utility functions.
 """
 
 # Standard Libraries
@@ -7,10 +10,12 @@
 import json
 from copy import deepcopy
 
+# Third-Party Libraries
 import numpy as np
 from ase import Atoms
 from pymatgen.core.structure import Molecule, Structure
 
+# PyXtal Modules
 from pyxtal.block_crystal import block_crystal
 from pyxtal.crystal import random_crystal
 from pyxtal.io import read_cif, structure_from_ext, write_cif
@@ -20,13 +25,12 @@
 from pyxtal.representation import representation, representation_atom
 from pyxtal.symmetry import Group, Wyckoff_position
 from pyxtal.tolerance import Tol_matrix
-
-# PyXtal imports #avoid *
 from pyxtal.version import __version__
 from pyxtal.viz import display_atomic, display_cluster, display_molecular
 from pyxtal.wyckoff_site import atom_site, mol_site
 from pyxtal.wyckoff_split import wyckoff_split
 
+# Uncomment the following line to set the package name
 # name = "pyxtal"
 
 
@@ -373,8 +377,8 @@ def from_random(
                     self.numMols = struc.numMols
                     self.molecules = struc.molecules
                     self.mol_sites = struc.mol_sites
-                    self.standard_setting = struc.mol_sites[0].wp.is_standard_setting(
-                    )
+                    wp = self.mol_sites[0].wp
+                    self.standard_setting = wp.is_standard_setting()
                 else:
                     self.numIons = struc.numIons
                     self.species = struc.species
@@ -387,6 +391,8 @@ def from_random(
             except:
                 pass
 
+            return struc
+
     def from_seed(
         self,
         seed,
@@ -3580,7 +3586,7 @@ def check_validity(self, criteria, verbose=False):
         if "Dimension" in criteria:
             try:
                 # TODO: unclear what is the criteria_cutoff
-                dim1 = self.get_dimensionality(criteria_cutoff)
+                dim1 = self.get_dimensionality(criteria['cutoff'])
             except:
                 dim1 = 3
             dim2 = criteria["Dimension"]
diff --git a/pyxtal/database/element.py b/pyxtal/database/element.py
index a1122ceb..503ba175 100644
--- a/pyxtal/database/element.py
+++ b/pyxtal/database/element.py
@@ -281,7 +281,7 @@ def __init__(self, input_value):
             self.covalent_radius = self.elements_list[pos][5]
             self.vdw_radius = self.elements_list[pos][6]
             self.metallic_radius = self.elements_list[pos][7]
-            if pos <= len(self.sf):
+            if pos < len(self.sf):
                 self.scatter = self.sf[pos]
 
     def get_all(self, pos):
diff --git a/pyxtal/db.py b/pyxtal/db.py
index c291f066..29b78722 100644
--- a/pyxtal/db.py
+++ b/pyxtal/db.py
@@ -148,7 +148,7 @@ def dftb_opt_single(id, xtal, skf_dir, steps, symmetrize, criteria, kresol=0.05)
             eng /= len(s)
 
         status = process_xtal(id, xtal, eng, criteria)
-        print(xtal.get_xtal_string(header=header, dicts=dicts))
+        print(xtal.get_xtal_string())
 
         return xtal, eng, status
     else:
@@ -614,14 +614,15 @@ def compute(self, row, work_dir, skf_dir):
             # not label information, run antechamber
             atom = self.db.get_atoms(id=row.id)
             if "gulp_info" not in row.data:
-                pmg, c_info, g_info = get_parameters(row, atom)
-                row.data = {"charmm_info": c_info, "gulp_info": g_info}
+                # pmg, c_info, g_info = get_parameters(row, atom)
+                # row.data = {"charmm_info": c_info, "gulp_info": g_info}
+                pass
             else:
                 pmg = ase2pymatgen(atom)
 
-            data = compute(row, pmg, work_dir, skf_dir)
-            self.db.update(row.id, data=data)
-            print("updated the data for", row.csd_code)
+            #data = compute(row, pmg, work_dir, skf_dir)
+            #self.db.update(row.id, data=data)
+            #print("updated the data for", row.csd_code)
 
 
 class database_topology:
@@ -1176,6 +1177,7 @@ def update_row_energy(
         use_relaxed=None,
         cmd=None,
         calc_folder=None,
+        skf_dir=None,
     ):
         """
         Update the row energy in the database for a given calculator.
@@ -1191,8 +1193,9 @@ def update_row_energy(
             ff_lib (str): Force field to use for GULP ('reaxff' by default).
             steps (int): Number of optimization steps for DFTB (default is 250).
             use_relaxed (str, optional): Use relaxed structures (e.g. 'ff_relaxed')
-            cmd (str, optional): Command for VASP calculations.
+            cmd (str, optional): Command for VASP calculations
             calc_folder (str, optional): calc_folder for GULP/VASP calculations
+            skf_dir (str, optional): Directory for DFTB potential files
 
         Functionality:
             Using the selected calculator, it updates the energy rows of the
diff --git a/pyxtal/interface/gulp.py b/pyxtal/interface/gulp.py
index 061789ea..a088d3dd 100644
--- a/pyxtal/interface/gulp.py
+++ b/pyxtal/interface/gulp.py
@@ -583,7 +583,7 @@ def write(self):
 
             # bond type
             if self.bond_type:
-                for bond in mol.molTopol.bonds:
+                for bond in site0.mol.molTopol.bonds:
                     for i in range(len(site.wp)):
                         count = i * len(coords)
                         f.write(f"connect {bond.atoms[0].id + count:4d} {bond.atoms[1].id + count:4d}\n")
diff --git a/pyxtal/lattice.py b/pyxtal/lattice.py
index 90ad5fb1..12b25ff8 100644
--- a/pyxtal/lattice.py
+++ b/pyxtal/lattice.py
@@ -4,9 +4,8 @@
 import numpy as np
 from numpy.random import Generator
 
-from pyxtal.constants import deg, ltype_keywords, rad
-
 # PyXtal imports
+from pyxtal.constants import deg, ltype_keywords, rad
 from pyxtal.msg import VolumeError
 from pyxtal.operations import angle, create_matrix
 
diff --git a/pyxtal/lego/SO3.py b/pyxtal/lego/SO3.py
index 7588e68b..9f29b64a 100644
--- a/pyxtal/lego/SO3.py
+++ b/pyxtal/lego/SO3.py
@@ -2,6 +2,8 @@
 import numpy as np
 from scipy.special import sph_harm, spherical_in
 from ase import Atoms
+from ase.neighborlist import NeighborList
+
 
 class SO3:
     '''
@@ -18,7 +20,7 @@ class SO3:
     '''
 
     def __init__(self, nmax=3, lmax=3, rcut=3.5, alpha=2.0,
-                 weight_on=False, neighborlist='ase'):
+                 weight_on=False):
         # populate attributes
         self.nmax = nmax
         self.lmax = lmax
@@ -27,7 +29,6 @@ def __init__(self, nmax=3, lmax=3, rcut=3.5, alpha=2.0,
         self._type = "SO3"
         self.cutoff_function = 'cosine'
         self.weight_on = weight_on
-        self.neighborcalc = neighborlist
         self.ncoefs = self.nmax*(self.nmax+1)//2*(self.lmax+1)
         self.tril_indices = np.tril_indices(self.nmax, k=0)
         self.ls = np.arange(self.lmax+1)
@@ -336,14 +337,9 @@ def build_neighbor_list(self, atom_ids=None):
             atom_ids = range(len(atoms))
 
         cutoffs = [self.rcut/2]*len(atoms)
-        if self.neighborcalc == 'ase':
-            from ase.neighborlist import NeighborList
-            nl = NeighborList(cutoffs, self_interaction=False, bothways=True, skin=0.0)
-            nl.update(atoms)
-        else:
-            from neighborlist import NeighborList
-            nl = NeighborList(cutoffs, self_interaction=False, bothways=True, skin=0.0)
-            nl.update(atoms, atom_ids)
+        nl = NeighborList(cutoffs, self_interaction=False, bothways=True, skin=0.0)
+        nl.update(atoms)
+
         #print(atoms, atom_ids)
         #print(atoms.get_scaled_positions())
 
diff --git a/pyxtal/lego/builder.py b/pyxtal/lego/builder.py
index 16073250..c28af612 100644
--- a/pyxtal/lego/builder.py
+++ b/pyxtal/lego/builder.py
@@ -514,7 +514,7 @@ def create_trajectory(dumpfile, trjfile, modes=['Init', 'Iter'], dim=3):
 
                 # QZ: to fix
                 xtal = pyxtal()
-                xtal = from_1d_rep(x, wps, numIons, l_type, elements, dim)
+                xtal.from_1d_rep(x, wps, numIons, l_type, elements, dim)
 
                 struc = xtal.to_ase()
                 struc.info = {'time': line_number-line_struc, 'fun': sim}
@@ -617,9 +617,8 @@ def set_descriptor_calculator(self, dtype='SO3', mykwargs={}):
                       'nmax': 2,
                       'rcut': 2.2,
                       'alpha': 1.5,
-                      # 'derivative': False,
                       'weight_on': True,
-                      'neighborlist': 'ase',
+                      #'neighborlist': 'ase',
                       }
             kwargs.update(mykwargs)
 
diff --git a/pyxtal/molecular_crystal.py b/pyxtal/molecular_crystal.py
index 5270a1df..62071ae3 100644
--- a/pyxtal/molecular_crystal.py
+++ b/pyxtal/molecular_crystal.py
@@ -4,7 +4,6 @@
 
 # Standard Libraries
 from copy import deepcopy
-
 import numpy as np
 
 from pyxtal.lattice import Lattice
@@ -454,6 +453,7 @@ def _set_mol_wyckoffs(self, id, numMol, pyxtal_mol, valid_ori, mol_wyks):
                                 return mol_sites_tmp
         return None
 
+
     def _set_orientation(self, pyxtal_mol, pt, oris, wp):
         """
         Generate valid orientations for a given molecule in a Wyckoff position.
@@ -465,19 +465,16 @@ def _set_orientation(self, pyxtal_mol, pt, oris, wp):
         - Using the bisection method is refine the orientation.
 
         Args:
-        pyxtal_mol: The pyxtal_molecule object representing the molecule.
-        pt: Position of the molecule.
-        oris: List of potential orientations.
-        wp: Wyckoff position object representing the symmetry of the site.
+            pyxtal_mol: The pyxtal_molecule object representing the molecule.
+            pt: Position of the molecule.
+            oris: List of potential orientations.
+            wp: Wyckoff position object representing the symmetry of the site.
 
         Returns:
-        ms0: A valid `mol_site` object if an acceptable orientation is found.
+            ms0: A valid `mol_site` object if an acceptable orientation is found.
              returns `None` if no valid orientation is found within the attempts.
         """
 
-        # Increment the number of attempts to generate a valid orientation
-        self.numattempts += 1
-
         # NOTE removing this copy causes tests to fail -> state not managed well
         ori = self.random_state.choice(oris).copy()
         ori.change_orientation(flip=True)
@@ -486,43 +483,12 @@ def _set_orientation(self, pyxtal_mol, pt, oris, wp):
         ms0 = mol_site(pyxtal_mol, pt, ori, wp, self.lattice)
 
         # Check if the current orientation results in valid distances
-        if ms0.short_dist():
+        if ms0.no_short_dist():
             return ms0
         else:
             # Maximize the separation if needed
             if len(pyxtal_mol.mol) > 1 and ori.degrees > 0:
-                # Define the distance function for bisection method
-                def fun_dist(angle, ori, mo, pt):
-                    # ori0 = ori.copy()
-                    ori.change_orientation(angle)
-                    ms0 = mol_site(mo, pt, ori, wp, self.lattice)
-                    return ms0.get_min_dist()
-
-                # Set initial bounds for the angle
-                angle_lo = ori.angle
-                angle_hi = angle_lo + np.pi
-                fun_lo = fun_dist(angle_lo, ori, pyxtal_mol, pt)
-                fun_hi = fun_dist(angle_hi, ori, pyxtal_mol, pt)
-
-                fun = fun_hi # Set the initial value for the function
-
-                # Refine the orientation using a bisection method
-                for _it in range(self.ori_attempts):
-                    self.numattempts += 1
-
-                    # Return as soon as a good orientation is found
-                    if (fun > 0.8) & (ms0.short_dist()):
-                        return ms0
-
-                    # Compute the midpoint angle for bisection
-                    angle = (angle_lo + angle_hi) / 2
-                    fun = fun_dist(angle, ori, pyxtal_mol, pt)
-
-                    # Update based on the function value at the midpoint
-                    if fun_lo > fun_hi:
-                        angle_hi, fun_hi = angle, fun
-                    else:
-                        angle_lo, fun_lo = angle, fun
+               return ms0.optimize_orientation_by_dist(self.ori_attempts)
 
     def _check_consistency(self, site, numMol):
         """
diff --git a/pyxtal/molecule.py b/pyxtal/molecule.py
index 51b67ce1..e4fd7511 100644
--- a/pyxtal/molecule.py
+++ b/pyxtal/molecule.py
@@ -1657,7 +1657,7 @@ def change_orientation(self, angle="random", flip=False):
 
             # Parse the angle
             if angle == "random":
-                angle = self.random_state.random() * np.pi * 2
+                angle = (self.random_state.random() - 1) * np.pi * 2
             self.angle = angle
 
             # Update the matrix
diff --git a/pyxtal/optimize/QRS.py b/pyxtal/optimize/QRS.py
index d1fc64f3..4a92b9b5 100644
--- a/pyxtal/optimize/QRS.py
+++ b/pyxtal/optimize/QRS.py
@@ -304,7 +304,7 @@ def _run(self, pool=None):
                     if self.early_termination(success_rate):
                         quit = True
 
-                elif ref_pxrd is not None:
+                elif self.ref_pxrd is not None:
                     self.count_pxrd_match(cur_xtals, matches)
 
             if self.use_mpi:
diff --git a/pyxtal/optimize/base.py b/pyxtal/optimize/base.py
index d66e483c..c3639759 100644
--- a/pyxtal/optimize/base.py
+++ b/pyxtal/optimize/base.py
@@ -928,7 +928,7 @@ def check_ref(self, reps=None, reference=None, filename="pyxtal.cif"):
                     pmg_s2 = representation(
                         rep[:-1], self.smiles).to_pyxtal().to_pymatgen()
                     pmg_s2.remove_species("H")
-                    if abs(eng1 - eng2) < 1e-2 and sm.StructureMatcher().fit(pmg_s1, pmg_s2):
+                    if abs(eng1 - eng2) < 1e-2 and self.matcher().fit(pmg_s1, pmg_s2):
                         new = False
                         break
                 if new:
@@ -982,7 +982,7 @@ def local_optimization(self, xtals, qrs=False, pool=None):
         elif self.ncpu == 1:
             return self.local_optimization_serial(xtals, qrs)
         else:
-            print(f"Local optimization by multi-threads {ncpu}")
+            print(f"Local optimization by multi-threads {self.ncpu}")
             return self.local_optimization_mproc(xtals, self.ncpu, qrs=qrs, pool=pool)
 
     def local_optimization_serial(self, xtals, qrs=False):
diff --git a/pyxtal/optimize/common.py b/pyxtal/optimize/common.py
index d5f4c040..522fab56 100644
--- a/pyxtal/optimize/common.py
+++ b/pyxtal/optimize/common.py
@@ -7,6 +7,7 @@
 from time import time
 
 import numpy as np
+from random import choice
 from ase import units
 
 from pyxtal import pyxtal
@@ -699,7 +700,7 @@ def load_reference_from_db(db_name, code=None):
     for rep in reps:
         rep = representation.from_string(rep, [smile])
         xtal1 = rep.to_pyxtal()
-        check_stable(xtal1, c_info, w_dir, skip_ani=True, optimizer=optimizer)
+        check_stable_structure(xtal1, c_info, w_dir, skip_ani=True, optimizer=optimizer)
 """
  81 11.38  6.48 11.24  96.9 1 0 0.23 0.43 0.03  -44.6   25.0   34.4  -76.6   -5.2  171.5 0 -70594.48
  81 11.38  6.48 11.24  96.9 1 0 0.23 0.43 0.03  -44.6   25.0   34.4  -76.6   -5.2  171.5 0 -70594.48
diff --git a/pyxtal/symmetry.py b/pyxtal/symmetry.py
index 2572f4ee..758117ad 100644
--- a/pyxtal/symmetry.py
+++ b/pyxtal/symmetry.py
@@ -1643,6 +1643,7 @@ def path_to_subgroup(self, H):
         paths = self.search_subgroup_paths(H)
         if len(paths) > 0:
             path = paths[0]
+            sg0 = path[0]
             pg0 = get_point_group(path[0])
             #pg0 = Group(path[0], quick=True)
             for p in path[1:]:
diff --git a/pyxtal/wyckoff_site.py b/pyxtal/wyckoff_site.py
index 40a7b423..cb44eb72 100644
--- a/pyxtal/wyckoff_site.py
+++ b/pyxtal/wyckoff_site.py
@@ -161,18 +161,6 @@ def encode(self):
         # print([self.specie, self.wp.index] + list(xyz))
         return [self.specie, self.wp.index, *list(xyz)]
 
-    def swap_axis(self, swap_id, shift=np.zeros(3)):
-        """
-        sometimes space groups like Pmm2 allows one to swap the a,b axes
-        to get an alternative representation
-        """
-        self.position += shift
-        self.position = self.position[swap_id]
-        self.position -= np.floor(self.position)
-        self.wp, _ = self.wp.swap_axis(swap_id)
-        self.site_symm = site_symm(
-            self.wp.symmetry[0], self.wp.number, dim=self.wp.dim)
-        self.update()
 
     def shift_by_swap(self, swap_id):
         """
@@ -198,8 +186,7 @@ def equivalent_set(self, tran, indices):
         self.position -= np.floor(self.position)
         self.wp = self.wp.equivalent_set(
             indices[self.wp.index])  # update the wp index
-        self.site_symm = site_symm(
-            self.wp.symmetry_m[0], self.wp.number, dim=self.wp.dim)
+        self.site_symm = self.wp.site_symm  # update the site symmetry
         self.update()
 
     def update(self, pos=None, reset_wp=False):
@@ -419,6 +406,43 @@ def update_orientation(self, angles):
         self.orientation.r = R.from_euler("zxy", angles, degrees=True)
         self.orientation.matrix = self.orientation.r.as_matrix()
 
+    def optimize_orientation_by_dist(self, ori_attempts):
+        """
+        Optimize the orientation according to the shortest distance
+        """
+        # Set initial fun value and angle bounds
+        ang_lo = self.orientation.angle
+        ang_hi = ang_lo + np.pi
+        fun_lo = self.get_min_dist(ang_lo)
+        fun_hi = self.get_min_dist(ang_hi)
+        fun = fun_hi
+
+        # Refine the orientation using a bisection method
+        for _it in range(ori_attempts):
+
+            # Return as soon as a good orientation is found
+            if (fun > 0.8) & self.no_short_dist():
+                return self
+
+            # Compute the midpoint angle for bisection
+            ang = (ang_lo + ang_hi) / 2
+            fun = self.get_min_dist(ang)
+
+            # Update based on the function value at the midpoint
+            if fun_lo > fun_hi:
+                ang_hi, fun_hi = ang, fun
+            else:
+                ang_lo, fun_lo = ang, fun
+
+            print("optimize_orientation_by_dist", _it, ang, fun)
+
+        return None
+
+    #def optimize_orientation_by_energy(self, ori_attempts):
+
+    #    V = -k*(d(D-A) - r)^2 + [sum_(cutoff*d**2)] when d<2
+
+
     def update_lattice(self, lattice):
         # QZ: Symmetrize the angle to the compatible orientation first
         self.lattice = lattice
@@ -959,7 +983,7 @@ def get_dists_WP(self, ignore=False, idx=None):
 
         return self.get_distances(coord1, coord2, ignore=ignore)
 
-    def get_min_dist(self):
+    def get_min_dist(self, angle=None):
         """
         Compute the minimum interatomic distance within the WP.
 
@@ -967,6 +991,9 @@ def get_min_dist(self):
             minimum distance
         """
         # Self image
+        if angle is not None:
+            self.orientation.change_orientation(angle)
+
         ds, _ = self.get_dists_auto()
         min_dist = np.min(ds)
 
@@ -981,7 +1008,7 @@ def get_min_dist(self):
                     min_dist = np.min(ds)
             return min_dist
 
-    def short_dist(self):
+    def no_short_dist(self):
         """
         Check if the atoms are too close within the WP.
 
@@ -1048,31 +1075,33 @@ def short_dist_with_wp2(self, wp2, tm=Tol_matrix(prototype="molecular")):
 
     def get_neighbors_auto(self, factor=1.1, max_d=4.0, ignore_E=True, detail=False, etol=-5e-2):
         """
-        Find neighboring molecules around the central molecule within a given distance threshold.
+        Find neighboring molecules within a given distance threshold.
 
         The function identifies neighboring molecules within PBC and computes the shortest
         distances and (optionally) interaction energies. it returns detailed information
         about neighboring molecule pairs, distances, and energies.
 
         Args:
-            factor (float, optional): Scaling factor for distance tolerances (default is 1.1).
-            max_d (float, optional): Maximum allowed intermolecular distance for neighbors (default is 4.0 Å).
-            ignore_E (bool, optional): If `True`, skips energy calculations (default is `True`).
-            detail (bool, optional): If `True`, returns detailed energy, molecular pairs, and distances
-                                     instead of just the shortest distances.
-            etol (float, optional): Energy tolerance for filtering pairs in detailed mode (default is -5e-2).
+            factor (float, optional): Scaling factor for distance (default is 1.1).
+            max_d (float, optional): Maximum distance for neighbors (default is 4.0 Å).
+            ignore_E (bool, optional): Skips energy calculations (default is `True`).
+            detail (bool, optional): Returns detailed eng, molecular pairs and
+                                     distances instead of only shortest distances.
+            etol (float, optional): Energy tolerance for filtering pairs in detailed
+                                    mode (default is -5e-2).
 
         Returns:
             If `detail == True`:
                 engs (list): List of interaction energies for valid molecular pairs.
-                pairs (list): List of tuples containing neighboring molecules and their relative positions.
+                pairs (list): List of tuples with neighbor molecules and positions.
                 dists (list): List of distances between neighboring molecular pairs.
             If `detail == False`:
-                min_ds (list): List of shortest distances between the central molecule and neighbors.
-                neighs (list): List of neighboring molecular coordinates (with PBC applied).
+                min_ds (list): List of shortest distances between central and neighbors.
+                neighs (list): List of neighboring molecular coordinates with PBC.
                 Ps (list): List of Wyckoff position multiplicities or translations.
-                engs (list): List of interaction energies, or `None` if energy calculation is skipped.
+                engs (list): List of energies, if energy calculation is not skipped.
         """
+
         # Compute the mol_center in Cartesian coordinate
         position = self.position - np.floor(self.position)
         mol_center = np.dot(position, self.lattice.matrix)
@@ -1083,7 +1112,7 @@ def get_neighbors_auto(self, factor=1.1, max_d=4.0, ignore_E=True, detail=False,
         # Get fractional coordinates for the central molecule
         coord1, _ = self._get_coords_and_species(first=True, unitcell=True)
 
-        # Initialize tolerance matrix for intermolecular distances (based on van der Waals radii)
+        # Initialize tolerance matrix for intermolecular distances
         tm = Tol_matrix(prototype="vdW", factor=factor)
         tols_matrix = self.molecule.get_tols_matrix(tm=tm)
 
@@ -1299,3 +1328,18 @@ def to_atom_site(self, specie=1):
         dicts["hn"] = self.wp.hall_number
         dicts["index"] = self.wp.index
         return atom_site.load_dict(dicts)
+
+
+if __name__ == "__main__":
+
+    from pyxtal.symmetry import Wyckoff_position as WP
+
+    wp = WP.from_group_and_letter(225, 'a')
+    coordinate = [0.25, 0.25, 0.25]
+    specie = 6
+    atom_site_instance = atom_site(wp=wp,
+                                   coordinate=coordinate,
+                                   specie=specie)
+
+    # Print the created instance
+    print(atom_site_instance)
diff --git a/pyxtal/wyckoff_split.py b/pyxtal/wyckoff_split.py
index 4f406b37..1f8c22e5 100644
--- a/pyxtal/wyckoff_split.py
+++ b/pyxtal/wyckoff_split.py
@@ -3,11 +3,9 @@
 """
 
 from copy import deepcopy
-
 import numpy as np
 from numpy.random import Generator
 from pymatgen.core.operations import SymmOp
-
 import pyxtal.symmetry as sym