matteoferla/Fragmenstein

View on GitHub
fragmenstein/openmm/openvictor.py

Summary

Maintainability
A
50 mins
Test Coverage
from ..victor import Victor
from .fritz import Fritz
from ..error import FragmensteinError
from rdkit import Chem
from rdkit.Chem import AllChem
from typing import List, Optional, Dict, Union, Callable
import os, json, time

class OpenVictor(Victor):
    uses_pyrosetta = False

    def _calculate_combination_thermo(self):
        return self._calculate_thermo()

    def _calculate_placement_thermo(self):
        return self._calculate_thermo()

    def _calculate_thermo(self):
        if self.is_covalent:
            raise NotImplementedError('OpenVictor does not support covalent ligands')
        self.journal.debug(f'{self.long_name} - Starting system setup')
        # in _calculate_*_chem this was set:
        # self.mol = self.monster.positioned_mol
        restraint_k = self.settings['mm_restraint_k'] # default 1000.
        tolerance = self.settings['mm_tolerance']  # 10 * mmu.kilocalorie_per_mole / (mmu.nano * mmu.meter)
        maxIterations = self.settings['mm_max_iterations'] # 0 is infinite
        self.fritz = Fritz(prepped_mol=self.preminimized_undummied_mol,
                           pdb_block=self.apo_pdbblock,
                           resn=self.ligand_resn,
                           restraining_atom_indices=self._get_restraining_atom_indices(),
                           restraint_k=restraint_k,
                           mobile_radius=self.settings['mm_mobile_radius'],
                           )
        self.unminimized_pdbblock = self.fritz.to_pdbblock()
        self._data: Dict = self.fritz.reanimate(tolerance=tolerance, maxIterations=maxIterations)
        while not self.quick_reanimation and self._data['binding_dG'] > 0 * self.fritz.molar_energy_unit:
            restraint_k = 2.
            self.fritz.alter_restraint(restraint_k / 3.)
            self._data: Dict = self.fritz.reanimate(tolerance=tolerance, maxIterations=maxIterations)
            self.journal.debug(f'{self.long_name} - restraints at {restraint_k}')
        self._data['restraint_k'] = restraint_k
        self._data['origins'] = self.monster.origin_from_mol()
        self.energy_score = {k: v.value_in_unit(self.fritz.molar_energy_unit)
                                    if hasattr(v, 'value_in_unit') else v for k, v in self._data.items()
                                                                 if k not in ('minimized_pdb', 'origins')}
        self.energy_score['unit'] = str(self.fritz.molar_energy_unit)
        self.minimized_pdbblock = self.fritz.to_pdbblock()
        self.minimized_mol = self.fritz.to_mol()
        self.mrmsd = self._calculate_rmsd()
        self.checkpoint()
        self.tock = time.time()
        return self.summarize()

    def _process_settings(self):
        pass

    def _get_restraining_atom_indices(self) -> List[int]:
        # Place has '_Stdev': 1.887379141862766e-15, '_Origin': '["x0395.9", "x0434.6"]', '_Max': 0.37818712299601354}
        # Combine has '_ori_i': 100, '_ori_name': 'x0395', '_x': 9.309, '_y': -5.402, '_z': 26.27
        # and occassionally '_ring_i'
        # The problem stemming from the collasing business
        # for now, all atoms are the same.
        restrainables: List[int] = []
        for atom in self.monster.positioned_mol.GetAtoms():
            if atom.GetSymbol() == 'H':
                continue
            elif atom.HasProp('_ori_name') or atom.HasProp('_Origin'):
                restrainables.append( atom.GetIdx() )
            else:
                # no restraint
                pass
        return restrainables

    def checkpoint(self, save_outputs: Optional[bool] = None):
        if save_outputs is True:
            pass
        elif self.settings.get('save_outputs', False):
            return
        self.journal.debug(f'{self.long_name} - saving data to disk at {self.work_path}')
        self.make_output_folder()
        with open(os.path.join(self.work_path, self.long_name, self.long_name + '.holo_unminimised.pdb'), 'w') as w:
            w.write(self.unminimized_pdbblock)
        with open(os.path.join(self.work_path, self.long_name, self.long_name + '.holo_minimised.pdb'), 'w') as w:
            w.write(self.minimized_pdbblock)
        for hit in self.hits:
            Chem.MolToMolFile(hit, os.path.join(self.work_path, self.long_name, hit.GetProp('_Name') + '.mol'))
        Chem.MolToMolFile(self.monster.positioned_mol,
                          os.path.join(self.work_path, self.long_name, self.long_name + '.positioned.mol'))
        Chem.MolToMolFile(self.fritz.prepped_mol,
                          os.path.join(self.work_path, self.long_name, self.long_name + '.prepped.mol'))
        Chem.MolToMolFile(self.minimized_mol,
                          os.path.join(self.work_path, self.long_name, self.long_name + '.minimised.mol'))
        with open(os.path.join(self.work_path, self.long_name, self.long_name + '.minimised.json'), 'w') as w:
            json.dump({**self.energy_score, 'origins': self._data['origins']}, w)

    def summarize(self):
        if self.error_msg:
            if self.monster is None:
                N_constrained_atoms = float('nan')
                N_unconstrained_atoms = float('nan')
            elif self.monster.positioned_mol is None:
                N_constrained_atoms = float('nan')
                N_unconstrained_atoms = float('nan')
            else:
                N_constrained_atoms = self.constrained_atoms
                N_unconstrained_atoms = self.unconstrained_heavy_atoms
            return {'name': self.long_name,
                    'smiles': self.smiles,
                    'error': self.error_msg,
                    'mode': self.merging_mode,
                    '∆∆G': float('nan'),
                    '∆G_bound': float('nan'),
                    '∆G_unbound': float('nan'),
                    'comRMSD': float('nan'),
                    'N_constrained_atoms': N_constrained_atoms,
                    'N_unconstrained_atoms': N_unconstrained_atoms,
                    'runtime': self.tock - self.tick,
                    'regarded': self.monster.matched,
                    'disregarded': self.monster.unmatched
                    }
        else:
            return {'name': self.long_name,
                    'smiles': self.smiles,
                    'error': self.error_msg,
                    'mode': self.merging_mode,
                    '∆∆G': self.energy_score.get('binding_dG', float('nan')),
                    'comRMSD': self.mrmsd.mrmsd,
                    'N_constrained_atoms': self.constrained_atoms,
                    'N_unconstrained_atoms': self.unconstrained_heavy_atoms,
                    'runtime': self.tock - self.tick,
                    'regarded': self.monster.matched,
                    'disregarded': self.monster.unmatched,
                    'origins': self._data['origins'],
                    }

    def _fix_minimized(self, *args, **kwargs):
        raise NotImplementedError('OpenVictor does not support covalent ligands')