From f7b60f03f0d914d23faebe90843e50862f98d134 Mon Sep 17 00:00:00 2001 From: m-titarenko Date: Sat, 13 Jun 2026 01:00:14 +0300 Subject: [PATCH] Fragment placer --- setup.py | 2 +- siman/core/molecule.py | 155 +++++++++++++++++++++++++++ tests/test10/acn.xyz | 8 ++ tests/test10/pf6.xyz | 9 ++ tests/test10/test_fragment_placer.py | 63 +++++++++++ 5 files changed, 236 insertions(+), 1 deletion(-) create mode 100644 tests/test10/acn.xyz create mode 100644 tests/test10/pf6.xyz create mode 100644 tests/test10/test_fragment_placer.py diff --git a/setup.py b/setup.py index f4c1f33..76cd231 100644 --- a/setup.py +++ b/setup.py @@ -18,7 +18,7 @@ packages=setuptools.find_packages(), # packages=['siman'], # data_files=[(os.path.expanduser("~"), ['siman/simanrc.py'])], - install_requires=['numpy', 'tabulate', 'pymatgen', 'pandas', 'scipy', 'six', 'matplotlib', 'ase', 'paramiko', 'adjustText', 'mp-api'], + install_requires=['numpy', 'tabulate', 'pymatgen', 'pandas', 'scipy', 'six', 'matplotlib', 'ase', 'paramiko', 'adjustText', 'mp-api', 'tblite'], classifiers=[ "Programming Language :: Python :: 3", "License :: OSI Approved :: GNU General Public License (GPL)", diff --git a/siman/core/molecule.py b/siman/core/molecule.py index 7779bc2..f338cca 100644 --- a/siman/core/molecule.py +++ b/siman/core/molecule.py @@ -8,6 +8,15 @@ from siman.small_functions import makedir from siman.header import printlog, runBash +import numpy as np +from ase.io import read +from ase.data import vdw_radii +from ase.optimize import BFGS +from scipy.spatial.distance import cdist +from scipy.spatial.transform import Rotation +from tblite.interface import Calculator +from tblite.ase import TBLite + class Molecule(Molecule_pymatgen): """Class for molecule structure representation based on pymatgen Molecule """ @@ -62,3 +71,149 @@ def jmol(self, program = 'jmol'): def get_volume(self): printlog('Molecule has no volume!', imp = 'n') return 0 + + +class FragmentPlacer: + """ + Places an ionic or molecular fragment near a host molecule + + The placement procedure consists of: + 1. Calculating of atomic charges by xTB. + 2. Selecting candidate binding sites based on atomic charges. + 3. Generating random fragment positions and orientations around the selected sites. + 4. Ranking generated poses using a scoring function. + 5. Performing an optional xTB geometry optimization of the best structure. + + Parameters: + clash_scale: float, scaling factor applied to vdw radii when detecting steric overlaps + n_sites: int, number of candidate binding sites considered + n_trials: int, number of random poses generated for site + distance_range: tuple(float, float), minimum and maximum placement distances (Å) + between the host site and the fragment center + """ + + def __init__(self, clash_scale=0.9, n_sites=5, n_trials=1000, distance_range=(2.0, 5.0)): + + self.clash_scale = clash_scale + self.n_sites = n_sites + self.n_trials = n_trials + self.distance_range = distance_range + + def place_fragment(self, host_file, fragment_file, fragment_charge, output=None, optimize=True): + """ + Place a fragment near a host molecule + + Parameters: + host_file: str, host structure file (all types, that can be read by ASE) + fragment_file: str, fragment structure file (-//-) + fragment_charge: int, charge of the fragment + output: str, file name for output structure + optimize (optional): If True, perform geometry optimization of final structure + + Returns: + ase.Atoms - final complex geometry + """ + + host = read(host_file) + fragment = read(fragment_file) + + host_q = self.get_xtb_charges(host) + frag_q = self.get_xtb_charges(fragment) + + sites = self.find_sites(host_q, fragment_charge) + + best_score = np.inf + best_pose = None + + for site in sites: + + for trial in range(self.n_trials): + + pose = self.generate_pose(host, fragment, site) + + score = self.score_pose(host, + host_q, + pose, + frag_q) + + if score < best_score: + best_score = score + best_pose = pose + + structure = host + best_pose + + if optimize: + structure = self.optimize(structure) + + structure.write(output) + + return structure + + def get_xtb_charges(self, atoms): + + calc = Calculator(method="GFN2-xTB", + numbers=atoms.numbers, + positions=atoms.positions) + res = calc.singlepoint() + + return np.array(res["charges"]) + + def find_sites(self, charges, fragment_charge): + + if fragment_charge < 0: + idx = np.argsort(-charges) + + elif fragment_charge > 0: + idx = np.argsort(charges) + + else: + idx = np.argsort(np.abs(charges))[::-1] + + return idx[:self.n_sites] + + def generate_pose(self, host, fragment, site_idx): + + frag = fragment.copy() + + coords = frag.positions - frag.get_center_of_mass() + coords = Rotation.random().apply(coords) + + direction = np.random.normal(size=3) + direction /= np.linalg.norm(direction) + + distance = np.random.uniform(*self.distance_range) + + target = host.positions[site_idx] + distance * direction + + frag.positions = coords + target + + return frag + + def score_pose(self, host, host_q, frag, frag_q): + + D = cdist(host.positions, frag.positions) + D[D < 0.5] = 0.5 + + electrostatic = np.sum(host_q[:, None] * frag_q[None, :] / D) + + r_host = np.array([vdw_radii[a.number] for a in host]) + r_frag = np.array([vdw_radii[a.number] for a in frag]) + + cutoff = self.clash_scale * (r_host[:, None] + r_frag[None, :]) + + overlap = np.clip(cutoff - D, 0.0, None) + penalty = np.sum(overlap**2) + + dmin = D.min() + + return electrostatic + 1000.0 * penalty + 0.1 * dmin + + def optimize(self, atoms, method="GFN2-xTB", fmax=0.05, steps=300): + + mol = atoms.copy() + mol.calc = TBLite(method=method) + + opt = BFGS(mol) + opt.run(fmax=fmax, steps=steps) + + return mol diff --git a/tests/test10/acn.xyz b/tests/test10/acn.xyz new file mode 100644 index 0000000..a78d823 --- /dev/null +++ b/tests/test10/acn.xyz @@ -0,0 +1,8 @@ +6 +6342 +N 1.26080 0.00000 0.00000 +C -1.36130 0.00000 0.00000 +C 0.10060 0.00000 0.00000 +H -1.75000 -0.83010 0.59740 +H -1.75010 -0.10220 -1.01750 +H -1.75000 0.93240 0.42020 diff --git a/tests/test10/pf6.xyz b/tests/test10/pf6.xyz new file mode 100644 index 0000000..0614109 --- /dev/null +++ b/tests/test10/pf6.xyz @@ -0,0 +1,9 @@ +7 +XYZ file generated by Avogadro. +F 0.00000 1.63533 0.00000 +P 0.00000 0.00000 0.00000 +F -0.00000 -1.63533 0.00000 +F 1.63533 -0.00000 0.00000 +F -1.63533 -0.00000 0.00000 +F -0.00000 0.00000 1.63532 +F 0.00000 0.00000 -1.63532 diff --git a/tests/test10/test_fragment_placer.py b/tests/test10/test_fragment_placer.py new file mode 100644 index 0000000..aee85e8 --- /dev/null +++ b/tests/test10/test_fragment_placer.py @@ -0,0 +1,63 @@ +""" +Test to check work of fragment_placer for molecules + +Author: Marina Titarenko + +To do: +Addability to check reaction sites and place ions/other molecules according to the reaction sites found +""" +from siman.core.molecule import FragmentPlacer + + +import os +import numpy as np +from ase.io import read +from ase.data import vdw_radii +from scipy.spatial.distance import cdist + + +def test(output, molecule_file): + #output - file with final geometry of complex + #molecule_file - str, file with initial geometry of host molecule + + if os.path.isfile(output): + print("CORRECT, Output file created") + else: + print("Output file missing") + return + + atoms = read(output) + molecule_atoms = read(molecule_file) + host_atoms = len(molecule_atoms) + + host = atoms[:host_atoms] + frag = atoms[host_atoms:] + + D = cdist(host.positions, frag.positions) + + r_host = np.array([vdw_radii[a.number] for a in host]) + r_frag = np.array([vdw_radii[a.number] for a in frag]) + + cutoff = 0.85 * (r_host[:, None] + r_frag[None, :]) + + if np.any(D < cutoff): + print("Atom overlaps detected") + else: + print("NO OVERLAPPING") + + +if __name__ == "__main__": + + molecule = "acn.xyz" + anion= "pf6.xyz" + fragment_charge = -1 + output = "acn_pf6_complex.xyz" + + placer = FragmentPlacer() + acn_comp = placer.place_fragment(host_file = molecule, + fragment_file = anion, + fragment_charge = -1, + output = output) + + test(output, molecule) +