Pyfrag old

Theoretical Chemistry – ACMM – Department of Chemistry – Faculty of Sciences – Vrije Universiteit Amsterdam

PyFrag

« Research
PyFrag Old
Research Tools »

“PyFrag – Streamlining your reaction path analysis” is available
in the Journal of Computational Chemistry.

Reference: W.J. van Zeist, C.Fonseca Guerra, F.M. Bickelhaupt; J. Comput. Chem. 29, 312-315, 2008.

Pyfrag Released!

The program PyFrag has been released and is available for download. PyFrag enables user-friendly and routine explorations and analyses of one- or multidimensional potential energy surfaces with the Amsterdam Density Functional (ADF) package. The Python source can be found in the download section. PyFrag is written by Willem-Jan van Zeist and is under development by Lando P. Wolters at the Theoretical Chemistry department at the Vrije Universiteit in Amsterdam. So far it has been mainly used for research in the subgroup of Prof. Dr. F.M. Bickelhaupt.

Why PyFrag?

Recently, the group of Prof. Dr. Bickelhaupt has focused on explaining chemical reactivity by means of the ‘Extended Activation Strain’ model. PyFrag serves as a valuable tool for this kind of research. Moreover, PyFrag facilitates an easy scan and analysis of any given (possibly multidimensional) potential energy surface. Such a PES can be constructed by PyFrag or can be taking from the output of for example Intrinisc Reaction Coordinate calculations. A short description of the activation strain model and the energy decomposition analysis with ADF is given below.

Extended Activation Strain Model & Energy Decomposition Analysis

Our group has used the Activation Strain model of chemical reactivity to gain insight into various fundamental chemical reactions such as SN(2) reactions and oxidative insertion/reductive elimination. In this model, the activation energy ∆E≠ is decomposed into the activation strain ∆E≠strain and the TS interaction ∆E≠int:

∆E≠= ∆E≠strain + ∆E≠int

The activation strain ∆E≠strain is the strain energy associated with deforming the reactants from their equilibrium geometry to the geometry they acquire in the activated complex. The TS interaction ∆E≠int is the actual interaction energy between the deformed reactants in the transition state. However, a problem is that the position of the TS along the reaction coordinate (ζ = ζTS) has a large effect on the magnitude of ∆E≠strain = ∆E≠strain(ζTS) and ∆E≠int = ∆E≠int(ζTS). It is therefore unsatisfactory to view the decomposition solely at the TS-geometry.

To tackle this problem, we have extended the Activation Strain model from a single-point analysis of the TS to an analysis along the reaction coordinate ζ:

∆E(ζ) = ∆Estrain(ζ) + ∆Eint(ζ)

The energy profile along the reaction coordinate ζ, can be for example computed using the Intrinsic Reaction Coordinate (IRC) method but any potential energy surface can apply. In the case of, for example, a reaction in which a bond is broken, the reaction coordinate can be suitably chosen as the bondlength of that particular bond.
The interaction ∆Eint(ζ) between the strained reactants is further analyzed in the conceptual framework provided by the Kohn-Sham molecular orbital model. Thus, ∆Eint is further decomposed into three physically meaningful terms using a quantitative energy decomposition scheme developed by Ziegler and Rauk.

∆Eint = ∆Velstat + ∆EPauli + ∆Eoi

The term ∆Velstat corresponds to the classical electrostatic interaction between the unperturbed charge distributions of the deformed reactants and is usually attractive. The Pauli repulsion ∆EPauli comprises the destabilizing interactions between occupied orbitals and is responsible for the steric repulsion. The orbital interaction ∆Eoi accounts for charge transfer (interaction between occupied orbitals on one moiety with unoccupied orbitals of the other, including the HOMO–LUMO interactions) and polarization (empty–occupied orbital mixing on one fragment due to the presence of another fragment).

For a list of references, see the Literature section.

Description

Pyfrag2007

The PyFrag program (originally released as PyFrag2007) is a ‘wrap-around’ for the Amsterdam Density Functional (ADF) package and facilitates the extension of the fragment analysis method implemented in ADF along an entire potential energy surfaces. The purpose is to make analyses of reaction paths and other (in principle also multidimensional) potential energy surfaces easier and more transparent. PyFrag also enables a user-friendly analysis of reaction paths in terms of the Extended Activation Strain model of chemical reactivity.

PyFrag = Python + ADF

PyFrag is a program written in the highly transferable Python programming language. Running it on any given system should pose no problems giving there is are working copies of Python and the ADF package present. PyFrag requires an ADF input script as a basis and constructs the necessary ADF jobs from it by reading additional PyFrag input statements. By reading or constructing a series of geometry points a series of ADF fragment analysis calculations will be performed, the results of which are gathered in a neat text data file. This facilitates a quick scan and analysis of any (possibly multidimensional) potential energy surface. For more info on usage and features, see the documentation page.

Documentation

Main Features & Usage

How does it work?
The main idea is that you make an input script much like the one you normally use for ADF, adding statements in a separate block that are understandable by PyFrag. In here you will define the fragments and the file from which PyFrag will collect geometrical data, which the program will use to set up the required ADF calculations.

How to run the program on the command line?
pyfrag.py pyfrag_job_script {temp_dir}

If no temporary directory is given, a suitable one will be created. The working directory, where the output files of pyfrag will be generated, is the one in which the pyfrag job script is present. It is the user’s own responsibility to set all the right ADF environment variables, so that ADF can run properly. In principle, nothing needs to be changed compared to running a single ADF job.

The script will look something like this:

PYFRAG INPUT

Type = IRC (or LT or SP)

Statements for pyfrag in whatever order you like.

# This will be a comment, ignored by PyFrag

END PYFRAG INPUT

$ADFBIN/adf << eor >

GEOMETRY

SINGLE POINT # It should be single point!

END

Other ADF statements of your choosing

END INPUT # The end input statement is important so that PyFrag nows when the ADF part ends.

eor

The program will substitute or add atom coordinates and copy TAPE21 statements so that they are appropriate for the calculations. Use tcsub with option -a pyfrag to submit the pyfrag input script in the usual manner. Error messages will be printed in the .err file. In the .out file, you will find the submitted adf input scripts. Outputfiles of individual calculations are stored in a subfolder (default /fragmentfiles). After the computation, the file ‘fragment_energies.txt’ in the main directory will contain all the data. TAPE21s will, by default, be removed.

The data file: fragment_energies.txt
By default this file will contain, as columns:
• The filenames of the fragment analysis calculations, which are numbered for each step in the calculations.
• The total interaction energy (kcal/mol)
• The Pauli repulsion (kcal/mol)
• The Electrostatic interaction (kcal/mol)
• The Total Steric interaction, Pauli + Electrostatic (kcal/mol)
• The Orbital interaction (kcal/mol)

By adding additional statements (such as print hirshfeld) more columns of data will be printed.

Main Specifications

(Unless otherwise stated, every key is by default disabled)

Firstly, on a line, everything after a ‘#’ will be ignored, allowing for comments in your inputfile.

Type = IRC / LT / SP (obligatory)
Type specifies mode in which pyfrag should be run. For analysis of an IRC calculation set type = IRC, for Linear transit set type = LT. Type = SP will set pyfrag in the mode in which you supply a geometry and single point calculations are performed with substitutions of one or more variables in a supplied geometry.

fa_name = name (default = fa_filename)
This will specify the name of the fragment analysis output files. Numbers will be added to the end of the filename for the different points, followed by a .o extension.
fragment dir #dir_name (default = fragment_files)
Specify the name for the directory in which the output files of the calculations will be place.
data file #file_name (default = fragment_energies.txt)
Specify the name of the data file that pyfrag generates.

save tape21s
Setting this option will generate a directory TAPE21s within your fragment directory, where the TAPE21 files of the calculations will be placed.

print bond #atomnr1 #atomnr2 #bond_diff
Specify the bond length to be printed for each geometry step, just indicate the atom numbers as they appear in the input order of the total molecule. Specifying bond_diff as well subtracts this value from the actual bond length.

print angle #atomnr1 #atomnr2 #atomnr3 #angle_diff
Specify the angle between atoms 1, 2 and 3 to be printed for each geometry step, just indicate the atom numbers as they appear in the input order of the total molecule. Specifying angle_diff as well subtracts this value from the actual angle.

print plane angle #atomnr1 #atomnr2 #atomnr3 #atomnr4 #atomnr5 #angle_diff
Specify the angle between the planes designated by the set of atoms 1, 2 and 3 and the set of atomes 3, 4 and 5 to be printed for each geometry step, just indicate the atom numbers as they appear in the input order of the total molecule. Specifying angle_diff as well subtracts this value from the actual angle.

print average distance
Prints the average change in cartesian coordinates after each step. This is done accumatively so after each step the average change in distance is increased.

print irrep oi
The value of the orbital interaction energy will be printed per available irrep.

print pauli decomp
Print de decomposition of the pauli energy term into its kinetic and electrostatic parts.

print strain frag1/frag2 (default disabled)
Specifies if you want the strain energy to be printed for the fragments. Specify the equilibrium energies (in kcal/mol) for the fragments. This value will then simply be subtracted from the energy of the fragment in question. The program will print the individual strain values for each fragment, plus the total strain and total energy.
Example:
print strain frag1 -301.01
print strain frag2 -19.02

print hirshfeld
This will print the hirshfeld charges for the two fragments. The order of the fragments as used internally by ADF may differ from what you would expect. Note also that hirshfeld charges as computed in a fragment analysis differ from those obtained in a ‘normal’ calculation.

Print VDD #atomnrs
Prints the VDD charges on the atoms with numbers as given by #atomnrs. Note also that VDD charges as computed in a fragment analysis differ from those obtained in a ‘normal’ calculation.
Example: print VDD 1 2

print dipole
Prints the dipole moment (debye).

Print SFO overlap #irrep #orb1|#orb2
This statement will print the overlap between orbital numbers #orb1 and #orb2 in irrep #irrep as they appear in the fragment analysis calculation. The numbers of the orbtials inserted for #orb have to include the core orbitals!
Examples:
print SFO overlap AA 17|18
print SFO overlap AAA 8|8
But also:
print SFO overlap HOMO|LUMO
print SFO overlap HOMO-1|LUMO+1
This will print the overlap between the HOMO’s of fragment one (frag1) and the LUMOs of fragment two (orbitals as found on the fragment calculations). If the orbitals found in this way differ in symmetry, an overlap value of zero is returned.

Print SFO population #irrep #orb (default disabled)
This statement will print the population of orbital number #orb in irrep #irrep. The numbers of the orbtitals inserted for #orb have to include the core orbitals!
Examples:
print SFO population AA 17
print SFO population AAA 8
But also:
Print SFO population frag2 HOMO
print SFO population frag1 LUMO

Print Orb Energy frag# #irrep #orb
This statement will print the energy of a fragment orbtial number #orb in irrep #irrep. The numbers of the orbtitals inserted for #orb have to include the core orbitals!
Examples:
print Orb Energy frag1 AA 17 —> Orbital numbers as on the fragment calculations!
print Orb Energy frag2 AAA 8
But also:
Print Orb Energy frag2 HOMO
print Orb Energy frag1 LUMO

Starting point #nr (default = 0)
If, for any reason, the calculation crashed at a certain cycle, you can specify a new, arbitrary, starting point. The number entered corresponds to the extension number on the fragment analysis output files that are generated. The program will continue to write data into the existing fragment_energies.txt file, simply appending new lines at the end of it.

points range #nrs, odd/even
With this you may more closely specify which points are used by pyfrag.
Example:
points range odd: Only use odd numbers (similar: even)
points range 1, 4-10: Use numbers 1 and 4 through 10.

densf frag# orb (feature in development)
Will perform densf and cntrs calculations for desired orbitals. At the moment this will simply generate a contour file of the specified orbitals, position in the x-y plane.
Example:
densf frag1 HOMO
densf frag2 LUMO
densf fa HOMO

print relative energies #nr
Makes a new data file with all energies and other values scaled to start from reference point nr.

EXTRA frag1/frag2/fa
ADF input statements here
END EXTRA frag1/frag2/fa
(default disabled)
You might want additional input for different parts of the calculation. Here frag1 and frag2 indicate the single point calculations for fragment 1 and 2, respectively. Extra fa inserts statements for the fragment analysis calculation. These statements will be useful when the charges on the fragments differ, in which case the following statements arise:

EXTRA frag1

CHARGE 0

END EXTRA frag1

EXTRA frag2

CHARGE 1

END EXTRA frag2

EXTRA fa

CHARGE 1

END EXTRA fa

Important is that, in this case, you shouldn’t have any charge statements in the ADF input file following the pyfrag input specifications which would result in double charge statements in the ADF scripts.

IRC and LT Specifications

The input specifications for analyzing IRC or Linear Transit output files are identical. You only need to specify one of them so the program reads the files in the correct way.

output file = #filename (obligatory)
output file 2 = #filename (optional)
Specify the filename of the LT or IRC output file you wish to analyze, both tape 21 and normal output files can be given. A TAPE21 file should have a .t21 extension! Please make sure of this. You can specify two output files for an IRC t21 (backward and forward). TAPE21 files containing only a backward or forward path should be indicated as such:
Only forward: filename = fw.t21
Only backward: filename = bw.t21

frag1/frag2 = name1/name2 (obligatory)

#atomnrs #element names
end frag1/frag2
These statements allows you to define the fragments, and the names that should be used for them, in the analysis. The use of the frag1 and frag2 statements is obligatory. For each of the fragments you supply a list of the numbers for the atoms as they exist on the LT or IRC file, followed by the name of the element. The program will check if they match with the order as in the IRC or LT calculation file. It will print a statement in the error log if it thinks the order is wrong.

frag1 = Methane_Substrate

1 C

2 H

4 H

5 H

6 H

end frag1

frag2 = PdPH3_Catalyst

3 Pd

7 P

8 H

9 H

10 H

end frag2

SP Specifications

The single point option allows you to directly run a series of single point fragment calculations. To do this, you can chose variables to be substituted into the geometry blocks given to pyfrag. If you have fixed fragments and only the total geometry of the molecule changes, you can supply the program with existing TAPE21s:

give tape21 frag1 = t21_name.t21
give tape21 frag2 = t21_name2.t21

Specify the fa_matrix as shown below, make sure you have defined a variable (BONDLENGTH by default).
In the case of an internal matrix, the matrix block should also incorporate the Geovar section. In case of internal coordinates use the ‘ATOMS INTERNAL’ header for the atoms block, and not ‘ATOMS ZMATRIX’, as this will confuse the program.
If you don’t have TAPE21s the program will calculate the fragments as it finds them in the atoms block according to the ‘f=frag1’ and ‘f=frag2’ statements. Note that this cannot be done if an internal matrix is supplied, since PyFrag will be unable to extract the fragment data from it.

matrix

ATOMS

O 0.0 0.0 0.0 f=frag1

H 1.0 0.0 0.0 f=frag2

H 0.0 BONDLENGTH 0.0 f=frag2

END

end matrix

The next step is of course to define how the variable BONDLENGTH (or other variables) should be substituted in each subsequent calculation. The easiest way is to simply put:

GeomVar = 1.0, 1.2, 1.5, 1.8, 2.0

This will by default then substitute these values for the variable BONDLENGTH wherever it is present (either fragments or fa_matrix).
Otherwise you can create a series of values by using:

GeomVar steps = 1.0, 2.0, 11

Where the values are, respectively, the starting bond length, the last bond length and the total number of calculations (steps) to perform. This example of course then produces the bond lengths 1.0, 1.1, 1.2 … 2.0.

In order to introduce more than one variable, simply add a variable name, for example:

GeomVar ANGLE = 10, 12, 14, 16, 18, 20
GeomVar steps BONDL1 = 1.0, 2.0, 6

In this way multiple variables may be substituted for each calculation run. The program will replace the variables wherever the are found (fragment blocks or fa_matrix). In the case of two variables, a grid of all possible combinations of the two values will be generated (compare this to the creation of a 2-dimensional PES). For more than 2 variables, please ensure that the total number of values is the same for all variables.

Downloads

Latest version - June 15, 2016 (bugfix for new version of ADF)

PyFrag2007.zip

Optimized for ADF 2016 and older. New version of PyFrag2017 can be found at here .

Please note that this version of program is no longer under development. For anyone who are intereted, please refere to PyFrag2017

Literature

For more information about PyFrag, visit http://www.few.vu.nl/~bickel/page-2/pyfrag.html

1. de Jong, G. T.; Bickelhaupt, F. M. ChemPhysChem 2007, Accepted.
2. Diefenbach, A.; de Jong, G. T.; Bickelhaupt, F. M. Mol Phys 2005, 103, 995-998.
3. Bickelhaupt, F. M. J Comp Chem 1999, 20, 114.
4. Baerends, E. J.; Autschbach, J. A.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; de Boeij, P. L.; Boerrigter, P. M.; Cavallo, L.; Chong, D. P.; Deng, L.; Dickson, R. M.; Ellis, D. E.; Fan, L.; Fischer, T. H.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Groeneveld, J. A.; Gritsenko, O. V.; Grüning, M.; Harris, F. E.; van den Hoek, P.; Jacob, C. R.; Jacobsen, H.; Jensen, L.; Kessel, G. v.; Kootstra, F.; van Lenthe, E.; McCormack, D. A.; Michalak, A.; Neugebauer, J.; Osinga, V. P.; Patchkovskii, S.; Philipsen, P. H. T.; Post, D.; Pye, C. C.; Ravenek, W.; Ros, P.; Schipper, P. R. T.; Schreckenbach, G.; Snijders, J. G.; Solà, M.; Swart, M.; Swerhone, D.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, R.; Wang, F.; Wesolowski, E.; van Wezenbeek, E.; Wiesenekker, G.; Wolff, S. K.; Woo, T. K.; Yakovlev, A. L.; Ziegler, T.; ADF 2006.01, SCM: Amsterdam, 2006.
5. te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; Van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. J Comput Chem 2001, 22, 931-967.
6. Bickelhaupt, F. M.; Baerends, E. J.; Nibbering, N. M. M. Chem Eur J 1996, 2, 196-207.
7. de Jong, G. T.; Visser, R.; Bickelhaupt, F. M. J Organomet Chem 2006, 691, 4341-4349.
8. Diefenbach, A.; de Jong, G. T.; Bickelhaupt, F. M. J Chem Theory Comput 2005, 1, 286.
9. Diefenbach, A.; Bickelhaupt, F. M. J Chem Phys 2001, 115, 4030.
10. Bickelhaupt, F. M.; Baerends, E. J. In Reviews in Computational Chemistry; Lipkowitz, K. B.; Boyd, D. B., Eds.; Wiley-VCH: New York, 2000.
11. Baerends, E. J.; Gritsenko, O. V. J Phys Chem A 1997, 101, 5383-5403.
12. Bickelhaupt, F. M.; Nibbering, N. M. M.; Van Wezenbeek, E. M.; Baerends, E. J. J Phys Chem 1992, 96, 4864-4873.
13. Ziegler, T.; Rauk, A. Theor Chim Acta 1977, 46, 1-10.
14. Ziegler, T.; Rauk, A. Inorg Chem 1979, 18, 1558-1565.
15. Ziegler, T.; Rauk, A. Inorg Chem 1979, 18, 1755-1759.
16. van Bochove, M. A.; Swart, M.; Bickelhaupt, F. M. J Am Chem Soc 2006, 128, 10738-10744.
17. W.J. van Zeist, A.H. Koers, L.P. Wolters, F.M. Bickelhaupt, J. Chem. Theory. Comput. 2008, 4, 920-928.

• Web of Science• NIST Chemistry Webbook
• Library
• VUNET

StatCounter - Free Web Tracker and Counter