Generating multiple conformers for a compound¶

import os
import matplotlib.pyplot as plt
from scm.conformers import ConformersJob
from scm.plams import from_smiles, init, finish, Settings, delete_job, KFFile, AMSJob, Molecule, CRSJob, Units

class COSKFConformers:
    '''
    Args:
        mol_identifier     (str)   : A smiles string or file (.xyz, .pdb, .mol, or mol2) with which to input a molecule
    Keyword Args:
        generator          (str)   : The method used to generate conformers.  Can be one of ['RDKit','CREST']
        initial_conformers (int)   : The number of initial conformers to sample
        coskf_name         (str)   : the filename base for the .coskf files this script produces.  multiple conformers will be saved with the filename base appended with a ``_n`` where n is the index of the conformer.
        max_energy_range   (float) : the largest acceptable energy (relative to the lowest energy conformer) for an adf calculation to be performed.  The energy is compared at the dftb level.
        max_num_confs      (int)   : the maximum number of conformers to produce .coskf files for
        print_output       (bool)  : Whether to print more detailed information about the progress of the script.
        dftb_solvent       (str)   : A solvent to use with GFN-xTB (DFTB).  Can be one of ['Acetone','Acetonitrile','CHCl3','CS2','DMSO','Ether','H2O','Methanol','THF','Toluene']
        coskf_dir          (str)   : The directory where produced .coskf files should be saved.  Defaults to the <plams run directory>/conformers
    '''
    def __init__(self,
                 mol_identifier : str,
                 generator          = 'RDKit',
                 initial_conformers = 500,
                 coskf_name         = None,
                 max_energy_range   = 2.0,
                 max_num_confs      = None,
                 print_output       = False,
                 dftb_solvent       = None,
                 coskf_dir          = None
                 ):

        self.print_output       = print_output
        self.coskf_name         = coskf_name
        self.generator          = generator
        self.initial_conformers = initial_conformers
        self.dftb_solvent       = dftb_solvent
        self.coskf_map          = {}
        self.coskf_dir          = coskf_dir

        self.mol                = self.process_molecule(mol_identifier)
        if self.mol is None:
            print("There was a problem reading the molecule: ", mol_identifier)
            return

        self.confs0 = self.get_conformers()
        self.confs1 = self.conformers_xtb()

        if max_num_confs and len(self.confs1.get_relative_energies())>max_num_confs:
            num_max_energy   = self.confs1.get_relative_energies()[max_num_confs-1]
            num_max_energy   = Units.convert(num_max_energy,'Hartree','kcal/mol') + 0.001
            max_energy_range = min(max_energy_range,num_max_energy) if max_energy_range else num_max_energy

        self.confs2 = self.conformers_adf(max_energy_range=max_energy_range)
        replay = self.conformers_cosmo()
        self.make_coskfs(replay)

        # Remove the results from the replay job because they can be large
        delete_job(replay.job)


    def process_molecule(self, mol_identifier):
        '''
        This function tries to process any molecular identifier passed to the COSKFConformers constructor.
        This first checks for a filename ending with .xyz, .pdb, .mol, or mol2.  If this is not found, the
        molecule is assumed to be input as a SMILES string.
        '''
        if any( mol_identifier.endswith(x) for x in ['.xyz','.pdb','.mol','mol2'] ):
            try:
                return Molecule(mol_identifier)
            except:
                return None
        try:
            return from_smiles(mol_identifier)
        except:
            return None

    def get_conformers(self):
        '''
        Optimize conformers at the UFF level
        '''
        sett = Settings()
        if self.generator == 'RDKit':
            sett.input.AMS.Generator.RDKit
            sett.input.AMS.Generator.RDKit.InitialNConformers = self.initial_conformers

        conformers_uff = ConformersJob(name="conformers_uff", molecule=self.mol).run()
        if self.print_output:
            print("Conformers at the UFF level:")
            print(conformers_uff)
        return conformers_uff

    def conformers_xtb(self):
        '''
        Optimize conformers with DFTB
        '''
        sett = Settings()
        sett.input.AMS.Task = "Optimize"
        sett.input.AMS.InputConformersSet = self.confs0
        sett.input.DFTB
        if self.dftb_solvent:
            sett.input.DFTB.Solvation.Solvent = self.dftb_solvent
        conformers_xtb = ConformersJob(name="conformers_xtb", settings=sett).run()
        if self.print_output:
            print(f"Conformers at the GFN1-xTB level with solvent {self.dftb_solvent}:")
            print(conformers_xtb)
        return conformers_xtb

    def conformers_adf(self,max_energy_range=None):
        '''
        Optimize conformers with ADF and the special settings for CRS
        '''
        sett = Settings()
        sett.input.AMS.Task = "Optimize"
        sett.input.AMS.GeometryOptimization.UseAMSWorker = "False"
        sett.input.AMS.InputConformersSet = self.confs1
        self._add_default_crs_settings(sett)
        if max_energy_range:
            sett.input.AMS.InputMaxEnergy = f'{max_energy_range} [kcal/mol]'
        conformers_adf = ConformersJob(name="conformers", settings=sett).run()
        if self.print_output:
            print("Conformers at the DFT level:")
            print(conformers_adf)
        return conformers_adf

    def conformers_cosmo(self):
        '''
        Replay job on the conformer set to generate all the ADF engine output files with the COSMO sections.
        These are then converted to .coskf files afterwards.
        '''
        sett = Settings()
        self._add_default_crs_settings(sett)
        sett.input.AMS.Task = "Replay"
        sett.input.AMS.Replay.File = self.confs2["conformers.rkf"]
        sett.input.AMS.Replay.StoreAllResultFiles = "True"
        sett.input.ADF.Symmetry = "NOSYM"
        self._add_default_crs_solvation(sett)
        replay = AMSJob(name="replay", settings=sett).run()
        return replay

    def make_coskfs(self,replay):
        '''
        Copy the COSMO sections from all the files back to the folder with the conformers
        '''
        base_name = self.coskf_name if self.coskf_name is not None else "conformer"
        if self.coskf_dir is None:
            self.coskf_dir = self.confs2.job.path
        for i, E in enumerate(self.confs2.get_energies("Ha")):
            if f"Frame{i+1}.rkf" in replay:
                cosmo_section = replay.read_rkf_section("COSMO", f"Frame{i+1}")
                cosmo_section["Gas Phase Bond Energy"] = E
                name = f"{base_name}_{i}.coskf"
                coskf = KFFile(os.path.join(self.coskf_dir, name), autosave=False)
                for key, val in cosmo_section.items():
                    coskf.write("COSMO", key, val)
                coskf.save()
                self.coskf_map[i] = name

    def _add_default_crs_settings(self,sett):
        '''
        This function adds the default settings for optimizing geometries in cosmo-rs/-sac.
        '''
        sett.input.ADF.Basis.Type        = "TZP"
        sett.input.ADF.Basis.Core        = "Small"
        sett.input.ADF.BeckeGrid.Quality = "Good"
        sett.input.ADF.XC.GGA            = "BP86"

    def _add_default_crs_solvation(self,sett):
        '''
        This function adds the solvation block required for generating sigma profiles.
        Note that in this section, only the radii of the most common elements have been added.  If your compound
        contains an element not listed below, you'll have to add its symbol and radius in the same format as the others.
        '''
        sett.input.ADF.Solvation = {
            "Surf": "Delley",
            "Solv": "Name=CRS Cav0=0.0 Cav1=0.0",
            "Charged": "Method=Conj Corr",
            "C-Mat": "Exact",
            "SCF": "Var All",
            "Radii": {
                "H": 1.30,
                "C": 2.00,
                "N": 1.83,
                "O": 1.72,
                "F": 1.72,
                "Si": 2.48,
                "P": 2.13,
                "S": 2.16,
                "Cl": 2.05,
                "Br": 2.16,
                "I": 2.32,
            },
        }

init()

coskf_name = "acetic_acid"
conformers = COSKFConformers(
    'CC(=O)O',
    print_output=True,
    coskf_name = coskf_name,
    max_energy_range=None,
    max_num_confs=100,
    dftb_solvent=None,
    coskf_dir='conformer_coskfs')

# here we'll do a simple cosmo-rs calculation with the conformers produced above
# this will be a binary mixture calculation with water (available in the cosmo-rs compound database)
settings = Settings()
settings.input.property._h = 'BINMIXCOEF'

# absolute paths are usually easier for input
path_to_confs = os.path.abspath(conformers.coskf_dir)

num_compounds = 2
compounds     = [Settings() for i in range(num_compounds)]

compounds[0].name = "acetic_acid"
form = [Settings() for i in range(len(conformers.coskf_map))]
for i,v in conformers.coskf_map.items():
    form[i]._h = os.path.join( path_to_confs, v )
compounds[0].form = form

compounds[1].name = "water"
compounds[1]._h   = os.path.join( os.getcwd(), "Water.coskf" )

settings.input.temperature = 298.15
settings.input.compound    = compounds
crs_job                    = CRSJob(settings=settings)
out                        = crs_job.run()
res                        = out.get_results()
comp_acetic                = out.get_multispecies_dist()[0]

mf1 = res['molar fraction'][0]
fig, axs = plt.subplots(2)

axs[0].plot( mf1, res['gamma'][0], label = '$\gamma_1$ (acetic)')
axs[0].plot( mf1, res['gamma'][1], label = '$\gamma_2$ (water)')

for struct, val in comp_acetic.items():
    axs[1].plot( mf1, val, label = os.path.basename(struct))

plt.setp(axs[0],ylabel='Activity coefficients')
plt.setp(axs[1],ylabel='Distribution ')
axs[0].legend()
axs[0].grid()
axs[1].legend()
axs[1].grid()
plt.xlabel('$x_1$ (acetic acid)')
plt.show()

finish()