Geometry optimization of water¶
If you’re a first-time PLAMS user, also check out the Getting started PLAMS tutorial!
Note
To execute this PLAMS script:
$AMSBIN/plams water_optimization.py
This example shows how to perform a geometry optimization of a water molecule and compute the vibrational normal modes using GFN1-xTB.
If you do not have a DFTB license, remove the line with DFTB settings
and instead set settings.input.ForceField.Type = 'UFF'
Initial imports¶
These two lines are not needed if you run PLAMS using the
$AMSBIN/plams program. They are only needed if you use
$AMSBIN/amspython.
from scm.plams import *
init()
PLAMS working folder: plams_workdir
Initial structure¶
# You could also load the geometry from an xyz file:
# molecule = Molecule('path/my_molecule.xyz')
# or generate a molecule from SMILES:
# molecule = from_smiles('O')
molecule = Molecule()
molecule.add_atom(Atom(symbol='O', coords=(0,0,0)))
molecule.add_atom(Atom(symbol='H', coords=(1,0,0)))
molecule.add_atom(Atom(symbol='H', coords=(0,1,0)))
try: plot_molecule(molecule) # plot molecule in a Jupyter Notebook in AMS2023+
except NameError: pass
Calculation settings¶
The calculation settings are stored in a Settings object, which is a
type of nested dictionary.
settings = Settings()
settings.input.ams.Task = 'GeometryOptimization'
settings.input.ams.Properties.NormalModes = 'Yes'
settings.input.DFTB.Model = 'GFN1-xTB'
#settings.input.ForceField.Type = 'UFF' # set this instead of DFTB if you do not have a DFTB license. You will then not be able to extract the HOMO and LUMO energies.
Create an AMSJob¶
job = AMSJob(molecule=molecule, settings=settings, name='water_optimization')
You can check the input to AMS by calling the get_input() method:
print("-- input to the job --")
print(job.get_input())
print("-- end of input --")
-- input to the job --
Properties
NormalModes Yes
End
Task GeometryOptimization
system
Atoms
O 0.0000000000 0.0000000000 0.0000000000
H 1.0000000000 0.0000000000 0.0000000000
H 0.0000000000 1.0000000000 0.0000000000
End
End
Engine DFTB
Model GFN1-xTB
EndEngine
-- end of input --
Run the job¶
job.run();
[23.01|18:55:53] JOB water_optimization STARTED
[23.01|18:55:53] JOB water_optimization RUNNING
[23.01|18:55:55] JOB water_optimization FINISHED
[23.01|18:55:55] JOB water_optimization SUCCESSFUL
Main results files: ams.rkf and dftb.rkf¶
The paths to the main binary results files ams.rkf and dftb.rkf
can be retrieved as follows:
print(job.results.rkfpath(file='ams'))
print(job.results.rkfpath(file='engine'))
plams_workdir/water_optimization/ams.rkf
plams_workdir/water_optimization/dftb.rkf
Optimized coordinates¶
optimized_molecule = job.results.get_main_molecule()
print("Optimized coordinates")
print("---------------------")
print(optimized_molecule)
print("---------------------")
Optimized coordinates
---------------------
Atoms:
1 O 0.066921 0.066921 0.000000
2 H 1.012042 -0.078963 0.000000
3 H -0.078963 1.012042 0.000000
---------------------
try: plot_molecule(optimized_molecule) # plot molecule in a Jupyter Notebook in AMS2023+
except NameError: pass
Optimized bond lengths and angle¶
Unlike python lists, where the index of the first element is 0, the index of the first atom in the molecule object is 1.
bond_length = optimized_molecule[1].distance_to(optimized_molecule[2])
print('O-H bond length: {:.3f} angstrom'.format(bond_length))
O-H bond length: 0.956 angstrom
bond_angle = optimized_molecule[1].angle(optimized_molecule[2], optimized_molecule[3])
print('Bond angle : {:.1f} degrees'.format(Units.convert(bond_angle, 'rad', 'degree')))
Bond angle : 107.5 degrees
Calculation timing¶
timings = job.results.get_timings()
print("Timings")
print("-------")
for key, value in timings.items():
print(f'{key:<20s}: {value:.3f} seconds')
print("-------")
Timings
-------
elapsed : 0.953 seconds
system : 0.045 seconds
cpu : 0.712 seconds
-------
Energy¶
energy = job.results.get_energy(unit='kcal/mol')
print('Energy : {:.3f} kcal/mol'.format(energy))
Energy : -3618.400 kcal/mol
Vibrational frequencies¶
frequencies = job.results.get_frequencies(unit='cm^-1')
print("Frequencies")
print("-----------")
for freq in frequencies:
print(f'{freq:.3f} cm^-1')
print("-----------")
Frequencies
-----------
1427.924 cm^-1
3674.507 cm^-1
3785.960 cm^-1
-----------
Dipole moment¶
import numpy as np
try:
dipole_moment = np.linalg.norm(np.array(job.results.get_dipolemoment()))
dipole_moment *= Units.convert(1.0, 'au', 'debye')
print('Dipole moment: {:.3f} debye'.format(dipole_moment))
except KeyError:
print("Couldn't extract the dipole moment")
Dipole moment: 1.830 debye
HOMO, LUMO, and HOMO-LUMO gap¶
Note: The methods for extracting HOMO, LUMO, and HOMO-LUMO gap only exist in AMS2023 and later.
try:
homo = job.results.get_homo_energies(unit='eV')[0]
lumo = job.results.get_lumo_energies(unit='eV')[0]
homo_lumo_gap = job.results.get_smallest_homo_lumo_gap(unit='eV')
print('HOMO : {:.3f} eV'.format(homo))
print('LUMO : {:.3f} eV'.format(lumo))
print('HOMO-LUMO gap : {:.3f} eV'.format(homo_lumo_gap))
except AttributeError:
print("Methods to extract HOMO and LUMO require AMS2023 or later")
except KeyError:
print("Couldn't extract the HOMO and LUMO.")
HOMO : -13.593 eV
LUMO : -4.206 eV
HOMO-LUMO gap : 9.387 eV
Read results directly from binary .rkf files¶
You can also read results directly from the binary .rkf files. Use the “expert mode” of the KFbrowser program that comes with AMS to find out which section and variable to read.
Below, we show how to extract the AMSResults%Energy variable from
the dftb.rkf file. This is the same number that was extracted previously
using the job.results.get_energy() method.
energy = job.results.readrkf('AMSResults', 'Energy', file='engine')
print(f"Energy from the engine .rkf file (in hartree): {energy}")
Energy from the engine .rkf file (in hartree): -5.766288141072482
Finish PLAMS¶
The finish() method is called automatically if you run the script
with $AMSBIN/plams. You should only call it if you use
$AMSBIN/amspython to run the script.
finish()
[23.01|18:55:55] PLAMS run finished. Goodbye