Source code for scm.plams.interfaces.adfsuite.quickjobs

import numpy as np

from scm.plams.core.functions import config, delete_job, finish, init
from scm.plams.core.settings import Settings
from scm.plams.interfaces.adfsuite.ams import AMSJob
from scm.plams.interfaces.adfsuite.amsworker import AMSWorker
from scm.plams.mol.molecule import Molecule
from typing import Optional

__all__ = ["preoptimize", "refine_density", "refine_lattice", "shakemd"]


[docs]def preoptimize( molecule: Molecule, model: str = "UFF", settings: Settings = None, nproc: int = 1, maxiterations: int = 100 ): """ Returns an optimized Molecule (or list of optimized molecules) molecule: Molecule or list of Molecules Molecule to optimize model: str Shorthand for some model, e.g. 'UFF' settings: Settings Custom engine settings (overrides ``model``) nproc: int Number of processes maxiterations: int Maximum number of iterations for the geometry optimization. """ my_settings = _get_quick_settings(model, settings, nproc) single_molecule = isinstance(molecule, Molecule) input_molecules = [molecule] if single_molecule else molecule output_molecules = [] with AMSWorker(my_settings) as worker: for i, mol in enumerate(input_molecules): results = worker.GeometryOptimization( name=f"preoptimization_{i}", molecule=mol, maxiterations=maxiterations, pretendconverged=True ) output_molecules.append(results.get_main_molecule()) if single_molecule: return output_molecules[0] else: return output_molecules
[docs]def shakemd( molecule: Molecule, density: Optional[float] = None, # kg/m^3 nsteps: int = 4000, model: str = "UFF", settings: Optional[Settings] = None, nproc: int = 1, ) -> Molecule: """ Performs SHAKE MD simulation with all bonds fixed. If ``density`` is specified, a lattice deformation will also be applied. The temperature is 100 K and the time step 0.5 fs. Returns: a Molecule corresponding to the final MD step. molecule: Molecule The initial structure density: float, optional If specified, the molecule must have a 3D lattice. Target density in kg/m^3 (1000x the density in g/cm^3) nsteps: int Number of MD steps. model: str e.g. 'UFF' settings: Settings Engine settings (overrides ``model``) """ called_plams_init = _ensure_init() s = _get_quick_settings(model, settings, nproc) s.input.ams.Task = "MolecularDynamics" s.input.ams.MolecularDynamics.NSteps = nsteps s.input.ams.MolecularDynamics.Thermostat.Type = "Berendsen" s.input.ams.MolecularDynamics.Thermostat.Tau = 100 s.input.ams.MolecularDynamics.Thermostat.Temperature = 100 s.input.ams.MolecularDynamics.TimeStep = 0.5 s.input.ams.MolecularDynamics.InitialVelocities.Temperature = 100 s.input.ams.MolecularDynamics.Shake.All = "bonds * *" if density is not None: temp_mol = molecule.copy() temp_mol.set_density(density) l = temp_mol.lattice target_lattice_str = f""" {l[0][0]} {l[0][1]} {l[0][2]} {l[1][0]} {l[1][1]} {l[1][2]} {l[2][0]} {l[2][1]} {l[2][2]} """ s.input.ams.MolecularDynamics.Deformation.StartStep = max(nsteps - 3000, 1) s.input.ams.MolecularDynamics.Deformation.StopStep = nsteps s.input.ams.MolecularDynamics.Deformation.TargetLattice._1 = target_lattice_str job = AMSJob(settings=s, molecule=molecule, name="shakemd") job.run() out = job.results.get_main_molecule() delete_job(job) if called_plams_init: finish() return out
[docs]def refine_density( molecule: Molecule, density: float, step_size=50, model: str = "UFF", settings: Settings = None, nproc: int = 1, maxiterations: int = 100, ): """ Performs a series of geometry optimizations with densities approaching ``density``. This can be useful if you want to compress a system to a given density, but cannot just use apply_strain() (because apply_strain() also scales the bond lengths). This function can be useful if for example packmol does not succeed to pack molecules with the desired density. Packmol can then generate a structure with a lower density, and this function can be used to increase the density to the desired value. Returns: a Molecule with the requested density. molecule: Molecule The molecule must have a 3D lattice density: float Target density in kg/m^3 (1000x the density in g/cm^3) step_size: float Step size for the density (in kg/m^3). Set step_size to a large number to only use 1 step. model: str e.g. 'UFF' settings: Settings Engine settings (overrides ``model``) maxiterations: int maximum number of iterations for the geometry optimization. """ assert len(molecule.lattice) == 3 tolerance = 1e-3 my_settings = _get_quick_settings(model, settings, nproc) current_density = molecule.get_density() # kg/m^3 output_molecule = molecule.copy() counter = 0 with AMSWorker(my_settings) as worker: while current_density < density - tolerance or current_density > density + tolerance: counter += 1 if current_density < density - step_size: new_density = current_density + step_size elif current_density > density + step_size: new_density = current_density - step_size else: new_density = density output_molecule.set_density(new_density) results = worker.GeometryOptimization( name=f"preoptimization_{counter}", molecule=output_molecule, maxiterations=maxiterations, pretendconverged=True, ) output_molecule = results.get_main_molecule() current_density = output_molecule.get_density() # kg/m^3 return output_molecule
[docs]def refine_lattice( molecule: Molecule, lattice, n_points=None, max_strain=0.15, model: str = "UFF", settings: Settings = None, nproc: int = 1, maxiterations: int = 10, ): """ Returns a ``Molecule`` for which the lattice of the ``molecule`` is transformed to ``lattice``, by performing short geometry optimizations (each for at most ``maxiterations``) on gradually distorted lattices (linearly interpolating from the original lattice to the new lattice using ``n_points`` points). This can be useful for transforming an orthorhombic box of a liquid into a non-orthorhombic box of a liquid, where the gradual transformation of the lattice ensures that the molecules do not become too distorted. If init() has been called before calling this function, the job will be run in the current PLAMS working directory and will be deleted when the job finishes. Returns: a Molecule with the requested lattice. If the refinement fails, ``None`` is returned. molecule: Molecule The initial molecule lattice: list of list of float List with 1, 2, or 3 elements. Each element is a list of float with 3 elements each. For example, ``lattice=[[10, 0, 0],[-5, 5, 0],[0, 0, 12]]``. n_points: None or int >=2 Number of points used for the linear interpolation. If None, n_points will be chosen such that the maximum strain for any step is at most ``max_strain`` compared to the original lattice vector lengths. max_strain: float Only if ``n_points=None``, use this value to determine the maximum allowed strain from one step to the next (as a fraction of the length of the original lattice vectors). model: str e.g. 'UFF' settings: Settings Engine settings (overrides ``model``) nproc: int Number of processes used by the job maxiterations: int maximum number of iterations for the geometry optimizations. """ assert len(lattice) == len( molecule.lattice ), f"Different number of lattice vectors: len(molecule.lattice) = {len(molecule.lattice)}, len(lattice) = {len(lattice)}" assert all(len(x) == 3 for x in lattice), f"Lattice vectors must have three components. Lattice: {lattice}" assert all(len(x) == 3 for x in molecule.lattice), f"Lattice vectors must have three components. Lattice: {lattice}" assert ( len(lattice) >= 1 and len(lattice) <= 3 ), f"{len(lattice)} lattice vectors given but must be between 1 and 3. Lattice: {lattice}" def lattice2str(latt): return "\n".join(" ".join(str(j) for j in i) for i in latt) if n_points is None: lat1 = np.array(molecule.lattice) lat2 = np.array(lattice) strain = np.linalg.norm(lat2, axis=1) / np.linalg.norm(lat1, axis=1) n_points = np.ceil(np.max(np.abs(strain - 1) / max_strain)) n_points = int(max(n_points, 2)) assert n_points >= 2, f"Expected n_points >=2, but received {n_points}" called_plams_init = _ensure_init() s = _get_quick_settings(model, settings, nproc) s.input.ams.task = "PESScan" s.input.ams.PESScan.ScanCoordinate.FromLattice._1 = lattice2str(molecule.lattice) s.input.ams.PESScan.ScanCoordinate.ToLattice._1 = lattice2str(lattice) s.input.ams.PESScan.ScanCoordinate.NPoints = n_points s.input.ams.GeometryOptimization.MaxIterations = maxiterations s.input.ams.GeometryOptimization.PretendConverged = "Yes" job = AMSJob(settings=s, molecule=molecule, name="refine_density") job.run() final_molecule = None if job.ok(): results = job.results.get_pesscan_results(molecules=True) final_molecule = results["Molecules"][-1] delete_job(job) if called_plams_init: finish() return final_molecule
def _ensure_init(): if config.init: called_plams_init = False else: s = Settings() s.erase_workdir = True s.log.stdout = 0 init(config_settings=s) called_plams_init = True return called_plams_init def model_to_settings(model: str): """ Returns Settings """ settings = Settings() model = model.lower() if model == "uff": settings.input.ForceField.Type = "UFF" elif model == "gaff": settings.input.ForceField.Type = "GAFF" settings.input.ForceField.AnteChamberIntegration = "Yes" elif model == "gfnff": settings.input.GFNFF elif model == "ani-2x" or model == "ani2x": settings.input.MLPotential.Model = "ANI-2x" elif model == "gfn1xtb" or model == "gfn1-xtb": settings.input.DFTB.Model = "GFN1-xTB" elif model == "m3gnet-up-2022" or model == "m3gnetup2022": settings.input.MLPotential.Model = "M3GNet-UP-2022" else: raise ValueError("Unknown model: {}".format(model)) return settings def _get_quick_settings(model, settings, nproc): if settings is None: my_settings = model_to_settings(model) else: my_settings = settings.copy() if nproc: my_settings.runscript.nproc = nproc return my_settings