diff --git a/pyxtal/__init__.py b/pyxtal/__init__.py index 41e0ce45..57d8e4c5 100644 --- a/pyxtal/__init__.py +++ b/pyxtal/__init__.py @@ -300,7 +300,7 @@ def from_random( pass - def from_seed(self, seed, molecules=None, tol=1e-4, a_tol=5.0, relax_H=False, ignore_HH=True, backend='pymatgen'): + def from_seed(self, seed, molecules=None, tol=1e-4, a_tol=5.0, ignore_HH=True, add_H=False, backend='pymatgen'): """ Load the seed structure from Pymatgen/ASE/POSCAR/CIFs Internally they will be handled by Pymatgen @@ -309,8 +309,8 @@ def from_seed(self, seed, molecules=None, tol=1e-4, a_tol=5.0, relax_H=False, ig seed: cif/poscar file or a Pymatgen Structure object molecules: a list of reference molecule (xyz file or Pyxtal molecule) tol: scale factor for covalent bond distance - relax_H: whether or not relax the position for hydrogen atoms in structure ignore_HH: whether or not ignore the short H-H distance when checking molecule + add_H: whether or not add H atoms """ @@ -319,7 +319,7 @@ def from_seed(self, seed, molecules=None, tol=1e-4, a_tol=5.0, relax_H=False, ig for mol in molecules: pmols.append(pyxtal_molecule(mol, fix=True)) #.mol #QZ: the default choice will not work for molecular H2, which should be rare! - struc = structure_from_ext(seed, pmols, relax_H=relax_H, ignore_HH=ignore_HH) + struc = structure_from_ext(seed, pmols, ignore_HH=ignore_HH, add_H=add_H) self.mol_sites = struc.make_mol_sites() self.group = Group(struc.wyc.number) self.lattice = struc.lattice diff --git a/pyxtal/io.py b/pyxtal/io.py index e2120bc8..2b5134c8 100644 --- a/pyxtal/io.py +++ b/pyxtal/io.py @@ -199,7 +199,7 @@ def read_cif(filename): class structure_from_ext(): - def __init__(self, struc, ref_mols, tol=0.2, relax_H=False, ignore_HH=False): + def __init__(self, struc, ref_mols, tol=0.2, ignore_HH=False, add_H=False): """ extract the mol_site information from the give cif file @@ -209,9 +209,8 @@ def __init__(self, struc, ref_mols, tol=0.2, relax_H=False, ignore_HH=False): struc: cif/poscar file or a Pymatgen Structure object ref_mols: a list of reference molecule (xyz file or Pyxtal molecule) tol: scale factor for covalent bond distance - relax_H: whether or not relax the position for hydrogen atoms in structure ignore_HH: whether or not ignore the short H-H distance when checking molecule - + add_H: whether or not add the H atoms """ for ref_mol in ref_mols: @@ -231,10 +230,13 @@ def __init__(self, struc, ref_mols, tol=0.2, relax_H=False, ignore_HH=False): print(type(struc)) raise NameError("input structure cannot be intepretted") + # reset the hydrogen position + if add_H: pmg_struc.remove_species('H') + self.ref_mols = ref_mols self.tol = tol self.diag = False - self.relax_H = relax_H + self.add_H = add_H sym_struc, number = get_symmetrized_pmg(pmg_struc) group = Group(number) @@ -243,7 +245,7 @@ def __init__(self, struc, ref_mols, tol=0.2, relax_H=False, ignore_HH=False): self.perm = [0,1,2] molecules = search_molecules_in_crystal(sym_struc, self.tol, ignore_HH=ignore_HH) - if self.relax_H: molecules = self.addh(molecules) + self.pmg_struc = sym_struc self.lattice = Lattice.from_matrix(sym_struc.lattice.matrix, ltype=group.lattice_type) self.resort(molecules) @@ -281,34 +283,36 @@ def resort(self, molecules): self.wps = [] ids_done = [] for j, mol2 in enumerate(self.ref_mols): + if self.add_H: + # create p_mol + p_mol = mol2.copy() + mol2.mol.remove_species("H") + for i, id in enumerate(ids): mol1 = molecules[id] if id not in ids_done and len(mol2.mol) == len(mol1): match, mapping = compare_mol_connectivity(mol2.mol, mol1) if match: - #print(self.numMols) - self.numMols[j] += mults[i] - #print(j, mults[i], self.numMols) - # rearrange the order - # create p_mol - p_mol = mol2.copy() + self.numMols[j] += mults[i] if len(mol1) > 1: + # rearrange the order order = [mapping[at] for at in range(len(mol1))] + xyz = mol1.cart_coords[order] + # add hydrogen positions here + if self.add_H: xyz = self.add_Hydrogens(mol2.smile, xyz) + #print(xyz) frac = np.dot(xyz, inv_lat) xyz = np.dot(frac, new_lat) center = p_mol.get_center(xyz) - #print(xyz-center) p_mol.reset_positions(xyz-center) position = np.dot(center, np.linalg.inv(new_lat)) else: position = np.dot(mol1.cart_coords[0], np.linalg.inv(new_lat)) position -= np.floor(position) - - #print(position) - #print(lat) + #print(position); print(lat) #print(p_mol.mol.cart_coords[:10] + np.dot(position, new_lat)) - # print(len(self.pmg_struc), len(self.molecule), len(self.wyc)) + # check if molecule is on the special wyckoff position if mults[i] < len(self.wyc): #Transform it to the conventional representation @@ -332,17 +336,34 @@ def resort(self, molecules): msg += molecules[id].to('xyz') raise ReadSeedError(msg) - def addh(self, molecules): + def add_Hydrogens(self, smile, xyz): """ - add hydrogen + add hydrogen for pymtagen molecule """ - from pymatgen.io.babel import BabelMolAdaptor - for mol in molecules: - ad = BabelMolAdaptor(mol) - ad.add_hydrogen() - ad.localopt() - mol = ad.pymatgen_mol - return molecules + from rdkit import Chem + from rdkit.Chem import AllChem + from rdkit.Geometry import Point3D + + #print(xyz) + m1 = Chem.MolFromSmiles(smile) + m2 = Chem.AddHs(m1) + AllChem.EmbedMolecule(m2,randomSeed=0xf00d) + m2 = Chem.RemoveHs(m2) + conf = m2.GetConformer(0) + + for i in range(conf.GetNumAtoms()): + x, y, z = xyz[i] + conf.SetAtomPosition(i, Point3D(x,y,z)) + #Chem.MolToPDBFile(m2, 'test.pdb') + #conf = m2.GetConformer(0); print(conf.GetPositions()) + + m1 = Chem.AddHs(m1) + AllChem.EmbedMolecule(m1) + AllChem.UFFOptimizeMolecule(m1) + m3 = AllChem.ConstrainedEmbed(m1,m2) + conf = m3.GetConformer(0) #; print(conf.GetPositions()) + + return conf.GetPositions() def make_mol_sites(self): """ diff --git a/pyxtal/representation.py b/pyxtal/representation.py index f74c01fb..e6f817f2 100644 --- a/pyxtal/representation.py +++ b/pyxtal/representation.py @@ -63,6 +63,10 @@ def from_string(cls, inputs, smiles, composition=None): n_cell = 6 elif g<=74: n_cell = 5 + elif g<=194: + n_cell = 4 + else: + n_cell = 3 #cubic cell = [g, diag] + inputs[2:n_cell] x = [cell] @@ -120,6 +124,8 @@ def to_pyxtal(self, smiles=None, composition=None): a, b, c, alpha, beta, gamma = v[2], v[3], v[4], 90, v[5], 90 elif ltype == 'orthorhombic': a, b, c, alpha, beta, gamma = v[2], v[3], v[4], 90, 90, 90 + elif ltype == 'tetragonal': + a, b, c, alpha, beta, gamma = v[2], v[2], v[3], 90, 90, 90 struc.lattice = Lattice.from_para(a, b, c, alpha, beta, gamma, ltype=ltype) # sites