Contents
Orient Energy Calculations
Main page Orient
See also Orient manual found at Stone Programs
Calculations of interaction energies using the Orient program require an input file that is best created using the Cluster program distributed as part of the CamCASP suite. This is recommended as Cluster will generally get things like units and structure right; leaving you with simpler edits only. But this article makes no assumption and describes the main functions and usage of Orient for interaction energy calculations.
Prerequisites
All calculations require
- Molecular geometries in XYZ format.
- One or more dimer configurations at which interaction energies are to be evaluated.
- Specification of UNITS used in various parts of the input.
Electrostatic energy
- Distributed multipoles. The molecular geometry is included in the multipole file. The multipoles are usually defined in the laboratory axis system.
Polarization energy
- Distributed multipoles. In laboratory axis system.
- Distributed polarizabilities. In either the laboratory axis system or a set of local axis frames.
- If local axes are used, these need to be specified in an Orient AXIS file.
- A damping model. Specify damping function type and damping parameters.
- Iterations used and convergence criterion.
Dispersion energy
- Dispersion model. This can be part of the potential specification.
- An anisotropic dispersion model can be defined either with respect to the laboratory frame or local axis frames. If the latter, the local axes need to be specified in an Orient AXIS file.
- A damping model.
Short-range repulsion
- Parameters specified in a potential file.
- If anisotropic, an AXIS file may be required if local axes are used.
Example water dimer
Here we calculate interaction energies for the water dimer. The entire file is located here. Save this file as H2O-dimer.ornt. The various parts of this file will be explained below. This calculation uses data in a few other files. These are
Multipole moments and molecular geometry H2O_daTZ_DMA2_L4.mom
Local axis definition H2O.axes
Polarizabilities defined in the local axis frame H2O_daTZ_ref_wt4_L1_static.pol
Dimer geometry file globalminconfig.geom
You will the H2O-dimer.ornt file and the four data files listed above all in one directory to run this job. First read the explanations below. Details for how this job is run are presented towards the end.
Preamble and Parameters
The Orient commmands begin with a definition of the UNITS and Parameters. UNITS can be specified and altered at any point in the file. We have specified more sites than are necessary. You usually need $N_A + N_B + 2$ sites, where $N_A$ and $N_B$ are the number of sites in molecules A and B. All our sites are polarizable. The other parameters are necessary for large molecules (not the case here).
UNITS BOHR
Parameters
Sites 20 polarizable 20
S-functions 50000
Alphas 50000
Parameter-sets 50000
Pairs 100000
End
Types
Next we define the TYPEs of atoms and their nuclear charges. The atom types are set in the molecule geometry block or in the multipole moments file.
What are Types? We have site names/labels, say, H1, and site types site H1 could be of type H as could site H2. Types allow us to associate multiple site labels with the same type of interaction. So when we define the potential it will typically be in terms of site types and not site labels. As there will typically be fewer site types than site labels, this allows us to make more compact definitions of the potential.
However, note the following
- Sites of the same type must be symmetry equivalent! For an anisotropic potential this would mean that we need local axes defined so as to observe all symmetries. This can be tricky for a large molecule.
- The molecular multipole moments and polarizabilities are specified for site labels and not site types!
In this example we have only two site types.
Types
O Z 8
H Z 1
End
Variables
Orient allows you to define variables which take on values that need to be defined at some point. This is useful for using the code to calculate energies for many dimer/cluster configurations.
In this example we define the dimer geometries using ANGLE-AXIS coordinates. These are specified as variables which will take on values for each dimer geometry. The variables have units B is Bohr and D is degrees. Angstroms is set with A.
Variables Rx 0.0 B Ry 0.0 B Rz 0.0 B alpha 0.0 D Nx 0.0 Ny 0.0 Nz 1.0 Index 0 End
Molecule definition
Now we come to the molecule specification. The first molecule H2O_A is located at the origin. I've commented out its coordinates as the are specified in the multipole moment file H2O_daTZ_DMA2_L4.mom which is #included into the Orient command file. The molecular site types are set in the moments file. Notice that we have defined the hydrogen atoms to have different types H1 and H2. This is done even though these atoms are symmetry-equivalent. The reason for this is that our polarizability file and potential specification uses the site types H1 and H2. So if we had given both hydrogen atoms Type H, Orient would be unable to find polarizabilities or potential parameters for the hydrogen atoms.
Next we have the molecular local axis system read in from H2O.axes. The multipole moments are specified before the local axes are loaded as they were calculated in the global axis system. They will be transformed into the local axis system. See the Orient output for details of the transformation and check if it all makes sense. Notice that the local axes observe the $C_{2v}$ symmetry of water, and consequently, the two hydrogen atoms can have the same type.
The final part of the molecule specification is the Polarizability definition. In this example we read localized polarizabilities of maximum rank 1 from file H2O_daTZ_ref_wt4_L1_static.pol. These have been calculated in the local axis system so they are read in after the local axes have been defined. The Limit command can be used to restrict the ranks on specified sites. Here it is redundant. The polarizabilities are defined for site labels and not site types.
Molecule H2O_A at 0.0 0.0 0.0
! Units BOHR
! O 0.00000000 0.00000000 0.00000000 Type O
! H1 -1.45365196 0.00000000 -1.12168732 Type H1
! H2 1.45365196 0.00000000 -1.12168732 Type H2
#include H2O_daTZ_DMA2_L4.mom
End
Edit H2O_A
#include H2O.axes
Bonds Auto
End
Polarizabilities for H2O_A
! Assumed that the pols are in the local-axes
! So they are read after axes are defined
Read rank 1
#include H2O_daTZ_ref_wt4_L1_static.pol
End
Limit rank 1 for O H1 H2
EndNext we have the definition of the second molecule (there could be more than two. Orient is not restricted to the number of molecules, though the Parameters would need to be consistent with the number of molecules.). The only difference here is the first line. We use the Variables defined above to specify the location and orientation of this molecule. Orient will insert values for these variables and perform the necessary translations/rotations. Variable names should be consistent with the names in the Variables section.
Molecule H2O_B at Rx Ry Rz ROTATED BY alpha ABOUT Nx Ny Nz
! Units BOHR
! O 0.00000000 0.00000000 0.00000000 Type O
! H1 -1.45365196 0.00000000 -1.12168732 Type H1
! H2 1.45365196 0.00000000 -1.12168732 Type H2
#include H2O_daTZ_DMA2_L4.mom
End
Edit H2O_B
#include H2O.axes
Bonds Auto
End
Polarizabilities for H2O_B
! Assumed that the pols are in the local-axes
! So they are read after axes are defined
Read rank 1
#include H2O_daTZ_ref_wt4_L1_static.pol
End
Limit rank 1 for O H1 H2
End
Potetnial definition
Potentials will normally be defined in atomic units. So we use the UNITS command to set atomic units. Then we define the potential using the Pairs section. The parameters for each pair of sites is defined. In this example each pair of site-types has a different induction damping parameter. Notice that the potential is defined for site-types and not site labels! This is important, and it is important that you have site types defined consistently (with symmetry (i.e., the axis file), moments and polarizabilities).
In this example the various parts of the potential (short-range repulsion, dispersion and damping parameter are defined separately. They could be combined.
Units Hartree Bohr
Pairs
! First we define the short-range repulsion parameters
! '''*This is not meant to be a very useful potential!
! It has been taken from an early stage in the fitting process
! and includes only first-order short-range effects. I.e., no charge-transfer
H H rho alpha
00 00 0 3.671346 1.765547
00 10 1 -0.316753
10 00 1 -0.316753
END
O H rho alpha
00 00 0 4.966920 1.922795
00 10 1 -0.264044
10 00 1 -0.016025
20 00 2 0.032885
END
O O rho alpha
00 00 0 5.621948 2.154523
00 10 1 -0.052563
10 00 1 -0.052563
00 20 2 -0.019649
20 00 2 -0.019649
END
! Now define the damping factors
! Dispersion uses a common damping for all site pairs, but the
! Polarization/Induction is damped according to site-type pairs.
Dispersion damping factor 1.93
H O
Induction damping factor 1.61
End
O O
Induction damping factor 1.09
End
H H
Induction damping factor 1.80
End
! The dispersion model (isotropic L1 from PBE0/AC daTZ properties)
O O C6
00 00 0 22.10714
End
H O C6
00 00 0 4.736933
End
H H C6
00 00 0 1.024315
End
EndSwitch units back to the units we need for the output energies. And set some control options for the polarization calculations.
Switch Induce On tells Orient to calculate induced moments. But Iterate Off says do not iterate them. This is what we need to calculate ${\rm Pol}^{(2)}$, that is, the second-order polarization. If we needed the infinite-order polarization we'd set Iterate On.
The Options section allows us to set various options. In this case we define convergence options for the Induction calculation. Iterations -1 sets no iterations. It is redundant as we have already set Iterate Off. If we need the infinite-order polarization (induction is a misnomer here) then we may need to set something like Iterations 100 or an appropriate value. Convergence set the convergence threshold. Iterations will stop when this threshold is reached. So the value of Iterations is a maximum that may not be reached.
Units Bohr kJ/mol Switch Induce On Iterate Off Options Induction Iterations -1 Convergence 1e-12 End
Energy Calculation
Finally we tell Orient to calculate energies and put them in a table using the Fortran format string f12.6. The energies to be calculated are defined using Print <options>. In this case we calculate the Exchange-Repulsion (er), Electrostatic energy (es), Polarization (termed Induction in Orient, hence called 'ind'), Dispersion (disp) and the total.
The Variables line is followed by a new line containing the list of variables to be read from the geometry file globalminconfig.geom which is #included.
And the file ends with a Finish.
Energy Table Format f12.6 Print er es ind disp total Variables Index Rx Ry Rz alpha Nx Ny Nz #include globalminconfig.geom End Finish
Output
Orient is run from the command line
$ orient < INPUT-FILE >& OUTPUT-FILE &
Here the input file contains all the commands above. The entire file is located here.
Put the commands in a file called, say, H2O-dimer.ornt and run the job as follows
$ orient < H2O-dimer.ornt >& H2O-dimer.out &
The entire output is in H2O-dimer.out. Here are some snippets
Notice how Orient has combined all the parts of the potential. You must check this to make sure the parameters are as desired. Multipole moments and polarizabilities are printed separately earlier in the output.
Pair potentials Default values Pre-exponential factor 1.0000E-03 Dispersion damping type 1 Parameters 1.9300 No induction damping Types O and O Pre-exponential factor 1.0000E-03 No dispersion damping Induction damping type 1 Parameters 1.0900 Index t1 t2 j rho alpha C6 1 00 00 0 5.621948 2.154523 22.11 2 00 10 1 -0.052563 0.000000 5 10 00 1 -0.052563 0.000000 8 00 20 2 -0.019649 0.000000 22 20 00 2 -0.019649 0.000000 Types O and H Pre-exponential factor 1.0000E-03 No dispersion damping Induction damping type 1 Parameters 1.6100 Index t1 t2 j rho alpha C6 1 00 00 0 4.966920 1.922795 4.737 2 00 10 1 -0.264044 0.000000 5 10 00 1 -0.016025 0.000000 22 20 00 2 0.032885 0.000000 Types H and O Pre-exponential factor 1.0000E-03 No dispersion damping Induction damping type 1 Parameters 1.6100 Index t1 t2 j rho alpha C6 1 00 00 0 4.966920 1.922795 4.737 2 00 10 1 -0.016025 0.000000 5 10 00 1 -0.264044 0.000000 8 00 20 2 0.032885 0.000000 Types H and H Pre-exponential factor 1.0000E-03 No dispersion damping Induction damping type 1 Parameters 1.8000 Index t1 t2 j rho alpha C6 1 00 00 0 3.671346 1.765547 1.024 2 00 10 1 -0.316753 0.000000 5 10 00 1 -0.316753 0.000000
And here are the energies. The current version of Orient gets the header labels wrong. Index should be first (it appears as a real number when it should be an integer) but it appears after the geometry parameters. This is followed by the energies. The total energies are not correct for the water dimer, but this is because the potential used here is just and example and is definitely not meant for work!
Rx Ry Rz alpha Nx Ny Nz Index es er ind disp total 900.00000 -3.87666 0.00000 -2.35391 133.12505 -0.63721 0.63721 -0.43351 -55.833023 229.678533 -8.658114 -40.595142 124.592254 901.00000 -4.03818 0.00000 -2.45199 133.12505 -0.63721 0.63721 -0.43351 -46.901033 159.374179 -6.939524 -28.312383 77.221239 902.00000 -4.83352 0.00000 -2.93492 133.12505 -0.63721 0.63721 -0.43351 -22.482482 26.398108 -2.376306 -6.587206 -5.047885 903.00000 -5.65346 0.00000 -3.43278 133.12505 -0.63721 0.63721 -0.43351 -12.208447 4.148708 -0.854375 -2.091778 -11.005892 904.00000 -6.46109 0.00000 -3.92318 133.12505 -0.63721 0.63721 -0.43351 -7.384505 0.672581 -0.346162 -0.834572 -7.892658 905.00000 -8.07636 0.00000 -4.90397 133.12505 -0.63721 0.63721 -0.43351 -3.290422 0.017849 -0.075598 -0.193682 -3.541853 906.00000 -9.69164 0.00000 -5.88477 133.12505 -0.63721 0.63721 -0.43351 -1.744248 0.000480 -0.022166 -0.061243 -1.827178 907.00000 -4.36124 0.00000 -2.64814 133.12505 -0.63721 0.63721 -0.43351 -34.046794 76.756850 -4.463543 -14.819180 23.427332 908.00000 -4.52276 0.00000 -2.74622 133.12505 -0.63721 0.63721 -0.43351 -29.354722 53.275979 -3.588536 -11.053625 9.279096 909.00000 -4.68429 0.00000 -2.84430 133.12505 -0.63721 0.63721 -0.43351 -25.481864 36.982355 -2.892406 -8.387054 0.221031 910.00000 -5.00735 0.00000 -3.04046 133.12505 -0.63721 0.63721 -0.43351 -19.549184 17.826988 -1.897049 -5.044556 -8.663802 911.00000 -5.16887 0.00000 -3.13854 133.12505 -0.63721 0.63721 -0.43351 -17.260144 12.379462 -1.544724 -3.986585 -10.411992 912.00000 -5.33040 0.00000 -3.23662 133.12505 -0.63721 0.63721 -0.43351 -15.312118 8.597690 -1.262793 -3.184926 -11.162147 913.00000 -3.55360 0.00000 -2.15775 133.12505 -0.63721 0.63721 -0.43351 -82.157311 477.065797 -13.407482 -91.668841 289.832164 914.00000 -3.71513 0.00000 -2.25583 133.12505 -0.63721 0.63721 -0.43351 -67.246970 331.013054 -10.789332 -59.956555 193.020196 Finished at 144518 on 9 May 2013
Polarization with iterations
In the above example we did not use iterations in the polarization model. Let's do so with the following code
Switch Induce On Iterate On Options Induction Iterations 100 Convergence 1e-12 End
This time the output is
Rx Ry Rz alpha Nx Ny Nz Index es er ind disp total 900.00000 -3.87666 0.00000 -2.35391 133.12505 -0.63721 0.63721 -0.43351 -55.833023 229.678533 -12.056577 -40.595142 121.193791 901.00000 -4.03818 0.00000 -2.45199 133.12505 -0.63721 0.63721 -0.43351 -46.901033 159.374179 -9.262320 -28.312383 74.898444 902.00000 -4.83352 0.00000 -2.93492 133.12505 -0.63721 0.63721 -0.43351 -22.482482 26.398108 -2.774507 -6.587206 -5.446086 903.00000 -5.65346 0.00000 -3.43278 133.12505 -0.63721 0.63721 -0.43351 -12.208447 4.148708 -0.933893 -2.091778 -11.085410 904.00000 -6.46109 0.00000 -3.92318 133.12505 -0.63721 0.63721 -0.43351 -7.384505 0.672581 -0.365950 -0.834572 -7.912446 905.00000 -8.07636 0.00000 -4.90397 133.12505 -0.63721 0.63721 -0.43351 -3.290422 0.017849 -0.077571 -0.193682 -3.543826 906.00000 -9.69164 0.00000 -5.88477 133.12505 -0.63721 0.63721 -0.43351 -1.744248 0.000480 -0.022479 -0.061243 -1.827490 907.00000 -4.36124 0.00000 -2.64814 133.12505 -0.63721 0.63721 -0.43351 -34.046794 76.756850 -5.572385 -14.819180 22.318491 908.00000 -4.52276 0.00000 -2.74622 133.12505 -0.63721 0.63721 -0.43351 -29.354722 53.275979 -4.363627 -11.053625 8.504005 909.00000 -4.68429 0.00000 -2.84430 133.12505 -0.63721 0.63721 -0.43351 -25.481864 36.982355 -3.438610 -8.387054 -0.325173 910.00000 -5.00735 0.00000 -3.04046 133.12505 -0.63721 0.63721 -0.43351 -19.549184 17.826988 -2.175090 -5.044556 -8.941843 911.00000 -5.16887 0.00000 -3.13854 133.12505 -0.63721 0.63721 -0.43351 -17.260144 12.379462 -1.745580 -3.986585 -10.612848 912.00000 -5.33040 0.00000 -3.23662 133.12505 -0.63721 0.63721 -0.43351 -15.312118 8.597690 -1.409088 -3.184926 -11.308442 913.00000 -3.55360 0.00000 -2.15775 133.12505 -0.63721 0.63721 -0.43351 -82.157311 477.065797 -20.815636 -91.668841 282.424009 914.00000 -3.71513 0.00000 -2.25583 133.12505 -0.63721 0.63721 -0.43351 -67.246970 331.013054 -15.793385 -59.956555 188.016143
Notice how the polarization (induction) energies have gotten larger in magnitude. This is now the infinite-order polarization ${\rm Pol}^{(2-\infty)}$. Iterations may not seem to result in a significant change in the polarization energy, but this is only the dimer. The effect of iterations in a cluster of water molecules can be as much as 20-30%.
Note
See Stone's book, chapter 10 in the second edition for more on many-body effects and the polarization.