Source code for scm.conformers.plams.interface

from os.path import join as opj

from scm.plams import SingleJob, Results, Molecule, KFFile, AMSJob, AMSResults, Settings, Units
from scm.plams.core.basejob import Job
from typing import List
import math

from ..conformers.conformers_rotamers import Conformers

__all__ = ["ConformersJob", "ConformersResults"]


[docs]class ConformersResults(Results):
    """
    A specialized |Results| subclass for accessing the results of ConformersJob.
    Conformers are sorted by energy, from lowest to highest.
    """

    rkfname = "conformers.rkf"

[docs]    def __init__(self, *args, **kwargs):
        Results.__init__(self, *args, **kwargs)
        self._conformers = None
        self.rkf = None

[docs]    def rkfpath(self) -> str:
        """Absolute path to the 'conformers.rkf' results file"""
        return opj(self.job.path, self.rkfname)

[docs]    def get_lowest_conformer(self) -> Molecule:
        """Return the conformer with the lowest energy"""
        return self._conformers.get_molecule(0)

[docs]    def get_conformers(self) -> List[Molecule]:
        """
        Return a list containing all conformers found.
        The conformers are sorted according to their energy, the first element being the lowest energy conformer.
        """
        return self._conformers.get_conformers()

[docs]    def get_relative_energies(self, unit="au") -> List[float]:
        """
        Return the relative energies of the conformers i.e. the energy of the conformer minus the energy of the lowest conformer found.
        This list is sorted according to the energy of the conformers, the first element corresponding to the lowest energy conformer.
        So, by definition, the first element will have an energy of 0.
        """
        e = self._conformers.get_energies()
        return Units.convert(e, "kcal/mol", unit)

[docs]    def get_energies(self, unit="au") -> List[float]:
        """
        Return the energies of the conformers.
        This list is sorted according to the energy of the conformers, the first element corresponding to the lowest energy conformer.
        """
        e = self._conformers.get_all_energies()
        return Units.convert(e, "kcal/mol", unit)

[docs]    def get_lowest_energy(self, unit="au") -> float:
        """Return the energy of the lowest-energy conformer."""
        return Units.convert(self._conformers.energies[0], "kcal/mol", unit)

[docs]    def get_boltzmann_distribution(self, temperature) -> List[float]:
        """
        Return the Boltzmann distribution at a given temperature: exp^(E_i/kB*temperature) / (sum_j exp^(E_j/kB*temperature)), where E_i is the energy of conformer i.
        This list is sorted according to the energy of the conformers, the first element corresponding to the lowest energy (and highest probability) conformer.
        The temperature is in Kelvin.
        """
        if temperature <= 0:
            raise ValueError(f"temperature ({temperature}) should be a positive number.")
        kB = 3.166819e-6  # Hartree / Kelvin
        weights = [math.exp(-e / (kB * temperature)) for e in self.get_relative_energies()]
        denominator = sum(weights)
        dist = [w / denominator for w in weights]
        return dist

[docs]    def __str__(self):
        return str(self._conformers)

    def get_energy_landscape(self):
        el = AMSResults.EnergyLandscape(None)
        for energy, mol in zip(self.get_energies(), self.get_conformers()):
            state = AMSResults.EnergyLandscape.State(el, None, energy, mol, 1, False)
            el._states.append(state)
        return el

[docs]    def collect(self):
        """Collect files present in the job folder.

        Use parent method from |Results| to get a list of all files in the results folder.
        Then instantiate ``self.rkf`` to be a |KFFile| instance for the main ``conformers.rkf`` output file.
        Also instantiate ``self._conformers`` to be a Conformers instance built from it.

        This method is called automatically during the final part of the job execution and there is no need to call it manually.
        """
        Results.collect(self)

        if self.rkfname in self.files:
            rkf_path = opj(self.job.path, self.rkfname)
            self.rkf = KFFile(rkf_path)
            self._conformers = Conformers()
            self._conformers.prepare_state(Molecule(rkf_path))
            self._conformers.read(rkf_path, reorder=False, filetype="rkf")


[docs]class ConformersJob(SingleJob):

    _result_type = ConformersResults
    _command = "conformers"

[docs]    def __init__(self, name="conformers", molecule=None, **kwargs):
        Job.__init__(self, name=name, **kwargs)

        if molecule is None:
            self.molecule = None
        elif isinstance(molecule, Molecule):
            self.molecule = molecule.copy()
        else:
            raise NotImplementedError("TODO: Implement multiple molecules input.")

[docs]    def check(self):
        return True

[docs]    def get_input(self):
        sett = self.settings.copy()

        # If there are references to ConformersResults in the settings, expand them into full paths:
        def expand_results_into_paths(s):
            for k, v in s.items():
                if isinstance(v, Settings):
                    expand_results_into_paths(v)
                elif isinstance(v, ConformersResults):
                    s[k] = v.rkfpath()
                elif isinstance(v, list) and any([isinstance(x, ConformersResults) for x in v]):
                    s[k] = [x.rkfpath() if isinstance(x, ConformersResults) else x for x in v]

        expand_results_into_paths(sett)

        return AMSJob(settings=sett, molecule=self.molecule).get_input()

[docs]    def get_runscript(self):
        ret = ""
        if "preamble_lines" in self.settings.runscript:
            for line in self.settings.runscript.preamble_lines:
                ret += f"{line}\n"
        ret += 'AMS_JOBNAME="{}" AMS_RESULTSDIR=. $AMSBIN/conformers'.format(self.name)
        if "nproc" in self.settings.runscript:
            ret += " -n {}".format(self.settings.runscript.nproc)
        ret += ' <"{}"'.format(self._filename("inp"))
        if self.settings.runscript.stdout_redirect:
            ret += ' >"{}"'.format(self._filename("out"))
        ret += "\n\n"
        return ret