add the trial function to make H-bond

doc/pyxtal.interface.ase_opt.rst

@@ -0,0 +1,7 @@
+pyxtal.interface.ase_opt module
+===============================
+
+.. automodule:: pyxtal.interface.ase_opt
+   :members:
+   :undoc-members:
+   :show-inheritance:

doc/pyxtal.interface.rst

@@ -13,7 +13,7 @@ Submodules
    :maxdepth: 4
 
    pyxtal.interface.LJ
-   pyxtal.interface.ani
+   pyxtal.interface.ase_opt
    pyxtal.interface.charmm
    pyxtal.interface.dftb
    pyxtal.interface.gulp

doc/requirements.txt

@@ -9,3 +9,4 @@ ase>=3.18.0
 sphinx_rtd_theme
 pyocse>=0.1.1
 toml
+mace

pyxtal/db.py

@@ -412,10 +412,10 @@ def make_db_from_CSD(dbname, codes):
 
 class database:
     """
-    This is a database class to process crystal data
+    This is a database class to process crystal data.
 
     Args:
-        db_name: *.db format from ase database
+        db_name: `*.db` format from ase database
     """
 
     def __init__(self, db_name):
@@ -630,7 +630,7 @@ class database_topology:
     This is a database class to process atomic crystal data
 
     Args:
-        db_name (str): *.db format from ase database
+        db_name (str): `*.db` format from ase database
         rank (int): default 0
         size (int): default 1
         ltol (float): lattice tolerance
@@ -992,6 +992,7 @@ def clean_structures_pmg(self, ids=(None, None), min_id=None, dtol=5e-2, criteri
         Clean up the db by removing the duplicate structures.
         Here we check the follow criteria same density and pymatgen matcher.
         The criteria should look like the following,
+
         {'CN': {'C': 3},
          'cutoff': 1.8,
          'MAX_energy': -8.00,

pyxtal/lego/builder.py

@@ -3,6 +3,7 @@
 building blocks (e.g., SiO2 with 4-coordined Si and 2-coordinated O)
 1. Generate all possible wp combinations
 2. For each wp combination,
+
     2.1 generate structure randomly
     2.2 optimize the geomtry
     2.3 parse the coordination

pyxtal/molecular_crystal.py

@@ -430,7 +430,7 @@ def _set_mol_wyckoffs(self, id, numMol, pyxtal_mol, valid_ori, mol_wyks):
                         pt[-1] = 0.5
 
                     ms0 = self._set_orientation(pyxtal_mol, pt, oris, wp)
-                    if ms0 is not None:
+                    if not ms0.short_dist():
                         # Check current WP against existing WP's
                         passed_wp_check = True
                         # print("Checking", ms0)
@@ -483,12 +483,10 @@ 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.no_short_dist():
-            return ms0
-        else:
-            # Maximize the separation if needed
-            if len(pyxtal_mol.mol) > 1 and ori.degrees > 0:
-               return ms0.optimize_orientation_by_dist(self.ori_attempts)
+        if len(pyxtal_mol.mol) > 1 and ori.degrees > 0:
+            if ms0.short_dist():
+                ms0.optimize_orientation_by_dist(self.ori_attempts)
+        return ms0
 
     def _check_consistency(self, site, numMol):
         """

pyxtal/molecule.py

@@ -1635,7 +1635,7 @@ def load_dict(cls, dicts):
         axis = dicts["axis"]
         return cls(matrix, degrees, axis)
 
-    def change_orientation(self, angle="random", flip=False):
+    def change_orientation(self, angle="random", flip=False, update=True):
         """
         Change the orientation of molecule by applying a rotation.
 
@@ -1652,25 +1652,38 @@ def change_orientation(self, angle="random", flip=False):
         if self.degrees >= 1:
             # Choose the axis
             if self.axis is None:
-                axis = self.random_state.random(3) - 0.5
-                self.axis = axis / np.linalg.norm(axis)
+                self.set_axis()
 
             # Parse the angle
             if angle == "random":
                 angle = (self.random_state.random() - 1) * np.pi * 2
-            self.angle = angle
+            #self.angle = angle
 
             # Update the matrix
-            r1 = Rotation.from_rotvec(self.angle * self.axis)
+            r1 = Rotation.from_rotvec(angle * self.axis)
 
             # Optionally flip the molecule
             if self.degrees == 2 and flip and self.random_state.random() > 0.5:
                 ax = self.random_state.choice(["x", "y", "z"])
                 angle0 = self.random_state.choice([90, 180, 270])
                 r2 = Rotation.from_euler(ax, angle0, degrees=True)
                 r1 = r2 * r1
-            self.r = r1 * self.r
-            self.matrix = self.r.as_matrix()
+
+            r = r1 * self.r
+            matrix = r.as_matrix()
+            if update:
+                self.r = r
+                self.matrix = matrix
+            return matrix
+        else:
+            return self.matrix
+
+
+    def set_axis(self):
+        if self.degrees == 2:
+            axis = self.random_state.random(3) - 0.5
+            self.axis = axis / np.linalg.norm(axis)
+            self.angle = 0
 
     def rotate_by_matrix(self, matrix, ignore_constraint=True):
         """