Merge pull request #279 from MaterSim/network

Running a systematic check with vscode

pyxtal/__init__.py

@@ -1,16 +1,21 @@
 """
-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
 import itertools
 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"]

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):

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

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")

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
 

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())
 

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)
 

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):
         """

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

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:

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):

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

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:]: