ADF Fragment Analysis Recipe¶
Use the PLAMS ADF fragment recipe to run two fragment calculations plus a full-system calculation and extract energy decomposition terms.
Initialization¶
from scm.plams import Settings, Molecule, AMSJob, Units, view, plot_image_grid
from scm.plams.recipes.adffragment import ADFFragmentJob
Define the molecules¶
For convenience we define here two molecules, normally you would read them from xyz files
def get_molecule(input_string):
job = AMSJob.from_input(input_string)
return job.molecule[""]
mol1 = get_molecule(
"""
System
Atoms
C -0.75086900 1.37782400 -2.43303700
C -0.05392100 2.51281000 -2.41769100
H -1.78964800 1.33942600 -2.09651100
H -0.30849400 0.43896500 -2.76734700
H -0.49177100 3.45043100 -2.06789100
H 0.98633900 2.54913500 -2.74329400
End
End
"""
)
mol2 = get_molecule(
"""
System
Atoms
C 0.14667300 -0.21503500 0.40053800
C 1.45297400 -0.07836900 0.12424400
C 2.23119700 1.15868100 0.12912100
C 1.78331500 2.39701500 0.38779700
H -0.48348000 0.63110600 0.67664100
H -0.33261900 -1.19332100 0.35411600
H 2.01546300 -0.97840100 -0.14506700
H 3.29046200 1.03872500 -0.12139700
H 2.45728900 3.25301000 0.35150400
H 0.74193400 2.60120700 0.64028800
End
End
"""
)
# Molecules could also be read from xyz files
# mol1 = Molecule("ethene.xyz")
# mol2 = Molecule("butadiene.xyz")
plot_image_grid(
{"1": view(mol1, guess_bonds=True), "2": view(mol2, guess_bonds=True)}, show_labels=False
);
Setup and run the job¶
common = Settings() # common settings for all 3 jobs
common.input.ams.Task = "SinglePoint"
common.input.adf.basis.type = "DZP"
common.input.adf.xc.gga = "PBE"
common.input.adf.symmetry = "NOSYM"
full = Settings() # additional settings for full system calculation
full.input.adf.etsnocv # empty block
full.input.adf.print = "etslowdin"
job = ADFFragmentJob(fragment1=mol1, fragment2=mol2, settings=common, full_settings=full)
results = job.run()
[08.04|17:13:19] JOB plamsjob STARTED
[08.04|17:13:19] JOB plamsjob RUNNING
[08.04|17:13:19] JOB plamsjob/frag1 STARTED
[08.04|17:13:19] JOB plamsjob/frag1 RUNNING
[08.04|17:13:22] JOB plamsjob/frag1 FINISHED
[08.04|17:13:22] JOB plamsjob/frag1 SUCCESSFUL
[08.04|17:13:22] JOB plamsjob/frag2 STARTED
[08.04|17:13:22] JOB plamsjob/frag2 RUNNING
[08.04|17:13:26] JOB plamsjob/frag2 FINISHED
[08.04|17:13:26] JOB plamsjob/frag2 SUCCESSFUL
[08.04|17:13:26] JOB plamsjob/full STARTED
[08.04|17:13:26] JOB plamsjob/full RUNNING
[08.04|17:13:38] JOB plamsjob/full FINISHED
[08.04|17:13:38] JOB plamsjob/full SUCCESSFUL
[08.04|17:13:38] JOB plamsjob FINISHED
[08.04|17:13:38] JOB plamsjob SUCCESSFUL
Energy Decomposition Results¶
Energy terms contributing to the bond energy (with respect to the fragments):
import pandas as pd
bond_energy = results.get_energy(unit="eV")
decom = results.get_energy_decomposition(unit="eV")
df = pd.DataFrame(
[
{"Term": t, "Energy [eV]": e}
for t, e in list(decom.items()) + [("Total Bond Energy", bond_energy)]
]
)
df
Term | Energy [eV] | |
|---|---|---|
0 | Electrostatic Energy | -0.159357 |
1 | Kinetic Energy | -0.295645 |
2 | Elstat Interaction | 0.748751 |
3 | XC Energy | -0.356182 |
4 | E_int | -1.697991 |
5 | E_Pauli | 0.187116 |
6 | E_elstat | -99.865788 |
7 | E_orb | -2.454467 |
8 | E_1 | -31.722875 |
9 | E_2 | -56.338062 |
10 | Total Bond Energy | -0.062434 |
See also¶
Python Script¶
#!/usr/bin/env python
# coding: utf-8
# ## Initialization
from scm.plams import Settings, Molecule, AMSJob, Units, view, plot_image_grid
from scm.plams.recipes.adffragment import ADFFragmentJob
# ## Define the molecules
# For convenience we define here two molecules, normally you would read them from xyz files
def get_molecule(input_string):
job = AMSJob.from_input(input_string)
return job.molecule[""]
mol1 = get_molecule(
"""
System
Atoms
C -0.75086900 1.37782400 -2.43303700
C -0.05392100 2.51281000 -2.41769100
H -1.78964800 1.33942600 -2.09651100
H -0.30849400 0.43896500 -2.76734700
H -0.49177100 3.45043100 -2.06789100
H 0.98633900 2.54913500 -2.74329400
End
End
"""
)
mol2 = get_molecule(
"""
System
Atoms
C 0.14667300 -0.21503500 0.40053800
C 1.45297400 -0.07836900 0.12424400
C 2.23119700 1.15868100 0.12912100
C 1.78331500 2.39701500 0.38779700
H -0.48348000 0.63110600 0.67664100
H -0.33261900 -1.19332100 0.35411600
H 2.01546300 -0.97840100 -0.14506700
H 3.29046200 1.03872500 -0.12139700
H 2.45728900 3.25301000 0.35150400
H 0.74193400 2.60120700 0.64028800
End
End
"""
)
# Molecules could also be read from xyz files
# mol1 = Molecule("ethene.xyz")
# mol2 = Molecule("butadiene.xyz")
plot_image_grid(
{"1": view(mol1, guess_bonds=True), "2": view(mol2, guess_bonds=True)}, show_labels=False, save_path="picture1.png"
)
# ## Setup and run the job
common = Settings() # common settings for all 3 jobs
common.input.ams.Task = "SinglePoint"
common.input.adf.basis.type = "DZP"
common.input.adf.xc.gga = "PBE"
common.input.adf.symmetry = "NOSYM"
full = Settings() # additional settings for full system calculation
full.input.adf.etsnocv # empty block
full.input.adf.print = "etslowdin"
job = ADFFragmentJob(fragment1=mol1, fragment2=mol2, settings=common, full_settings=full)
results = job.run()
# ## Energy Decomposition Results
# Energy terms contributing to the bond energy (with respect to the fragments):
import pandas as pd
bond_energy = results.get_energy(unit="eV")
decom = results.get_energy_decomposition(unit="eV")
df = pd.DataFrame(
[{"Term": t, "Energy [eV]": e} for t, e in list(decom.items()) + [("Total Bond Energy", bond_energy)]]
)
print(df)