#!/usr/bin/env amspython from scm.plams import * import numpy as np def main(): """ Three calculations: 1. PESScan providing initial guess for the TS search 2. TS Search + frequencies 3. Reactants (single point) + frequencies The first time you run this script, set the previous_pesscan_job_dir, previous_ts_job_dir, and previous_reactants_job_dir to None. That will run the calculations. If you have already run the calculations, set the variables to the directories containing the previous jobs to just load results from them For example previous_pesscan_job_dir = 'plams_workdir/pes_scan' previous_ts_job_dir = 'plams_workdir/ts_search' previous_reactants_job_dir = 'plams_workdir/reactants_frequencies' """ init() previous_pesscan_job_dir = None previous_ts_job_dir = None previous_reactants_job_dir = None if previous_ts_job_dir: ts_job = AMSJob.load_external(previous_ts_job_dir) else: if previous_pesscan_job_dir: pesscan_job = AMSJob.load_external(previous_pesscan_job_dir) else: pesscan_job = run_pesscan_job('ethene-metallocenium.xyz', charge=+1.0) ts_guess_molecule = get_highest_energy_molecule(pesscan_job) ts_job = run_ts_job(ts_guess_molecule) print_ts_results(ts_job) gibbs_energy_ts = get_gibbs_energy(ts_job.results, unit='kcal/mol') if previous_reactants_job_dir: reactants_job = AMSJob.load_external(previous_reactants_job_dir) else: reactants_job = run_reactants_job('ethene-metallocenium.xyz', charge=+1.0) gibbs_energy_reactants = get_gibbs_energy(reactants_job.results, unit='kcal/mol') print_final_results(gibbs_energy_ts, gibbs_energy_reactants, T=298.15) finish() def get_dftb_settings(): # Engine settings: GFN-1xTB with implicit Tolene solvation settings = Settings() settings.Model = 'GFN1-xTB' settings.Solvation.Solvent = 'Toluene' return settings def run_reactants_job(xyz_file, charge=0.0): molecule = Molecule(xyz_file) molecule.properties.charge = charge settings = Settings() settings.input.dftb = get_dftb_settings() settings.input.ams.Task = 'SinglePoint' settings.input.ams.Properties.NormalModes = 'Yes' job = AMSJob(molecule=molecule, settings=settings, name='reactants_frequencies') job.run(watch=True) return job def run_pesscan_job(xyz_file, charge=0.0): # Generate a molecule object from xyz coordinates molecule = Molecule(xyz_file) molecule.properties.charge = charge settings = Settings() # Engine settings: GFN-1xTB with implicit Tolene solvation settings.input.dftb = get_dftb_settings() # AMS settings for PES scan (see online tutorial for choice of scan coordinates) settings.input.ams.Task = 'PESScan' settings.input.ams.PESScan.ScanCoordinate.nPoints = '20' settings.input.ams.PESScan.ScanCoordinate.Distance = ['1 28 3.205 2.3', '29 25 3.295 1.5'] # Loosened convergence criteria settings.input.ams.GeometryOptimization.Convergence.Energy = '5.0e-5 [Hartree]' settings.input.ams.GeometryOptimization.Convergence.Gradients = '5.0e-3' settings.input.ams.GeometryOptimization.Convergence.Step = '5.0e-3 [Angstrom]' # setting up the AMS job with the molecule and settings object job = AMSJob(molecule=molecule, settings=settings, name='pes_scan') # run the calulation job.run(watch=True) return job def get_highest_energy_molecule(pesscan_job : AMSJob): res = pesscan_job.results.get_pesscan_results() # res['PES'] contains all the energies from the PESScan # see the AMSResults PLAMS documentation for more information about the other keys in res highest_energy = np.max(res['PES']) highest_energy_index = np.argmax(res['PES']) # res['Molecules'] contains all the converged geometries (as PLAMS Molecules) highest_energy_geo = res['Molecules'][highest_energy_index] print("Highest energy in the PESScan job ({:.6f} hartree) was found for converged geometry #{}".format(highest_energy, highest_energy_index+1)) return highest_energy_geo def run_ts_job(molecule : Molecule): settings = Settings() # Engine settings: GFN-1xTB with implicit Tolene solvation settings.input.dftb = get_dftb_settings() # AMS settings for TS search settings.input.ams.Task = 'TransitionStateSearch' settings.input.ams.Properties.NormalModes = 'Yes' settings.input.ams.GeometryOptimization.InitialHessian.Type = 'Calculate' # setting up the AMS job with the molecule and settings object job = AMSJob(molecule=molecule, settings=settings, name='ts_search') # run the calulation result = job.run(watch=True) return job def get_gibbs_energy(result: AMSResults, unit='hartree'): gibbs_energy = result.readrkf('Thermodynamics', 'Gibbs free Energy', file='engine') gibbs_energy *= Units.convert(1, 'hartree', unit) return gibbs_energy def print_ts_results(ts_job): result = ts_job.results energy = result.get_energy(unit='kcal/mol') gibbs_energy = get_gibbs_energy(result, unit='kcal/mol') frequencies = result.get_frequencies(unit='cm^-1') optimized_molecule = result.get_main_molecule() # read the classification of the PES point from file dftb.rkf character = result.readrkf('AMSResults', 'PESPointCharacter', file='engine') # print results print('== Results ==') print('Transition state geometry:') print(optimized_molecule) print('Character : {}'.format(character)) print('Energy : {:.3f} kcal/mol'.format(energy)) print('Gibbs Energy: {:.3f} kcal/mol'.format(gibbs_energy)) for freq in frequencies: if freq < 0: print('Imaginary frequency : {:.3f} cm^-1'.format(freq)) else: break print('End of TS job summary') print('------------------------') def print_final_results(gibbs_energy_ts : float, gibbs_energy_reactants : float, T=298.15): barrier = gibbs_energy_ts - gibbs_energy_reactants # kcal/ mol R = 1.9872 # cal / mol*K k_B = Units.constants['Boltzmann'] # J / (mol*K) h = 6.626e-34 # J * s rate_constant = (k_B*T/h) * np.exp(-barrier*1000/(R*T)) print("Results summary") print("Reactants Gibbs Energy: {:.3f} kcal/mol".format(gibbs_energy_reactants)) print("TS Gibbs Energy: {:.3f} kcal/mol".format(gibbs_energy_ts)) print("Free energy barrier: {:.3f} kcal/mol".format(barrier)) print("Rate constant: {:.3f} s⁻¹ at {} K".format(rate_constant, T)) if __name__ == '__main__': main()