## Orient Energy Calculations

Main page Orient

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

1. Molecular geometries in XYZ format.
2. One or more dimer configurations at which interaction energies are to be evaluated.
3. Specification of UNITS used in various parts of the input.

Electrostatic energy

1. Distributed multipoles. The molecular geometry is included in the multipole file. The multipoles are usually defined in the laboratory axis system.

Polarization energy

1. Distributed multipoles. In laboratory axis system.
2. Distributed polarizabilities. In either the laboratory axis system or a set of local axis frames.
3. If local axes are used, these need to be specified in an Orient AXIS file.
4. A damping model. Specify damping function type and damping parameters.
5. Iterations used and convergence criterion.

Dispersion energy

1. Dispersion model. This can be part of the potential specification.
2. 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.
3. A damping model.

Short-range repulsion

1. Parameters specified in a potential file.
2. 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

1. Multipole moments and molecular geometry H2O_daTZ_DMA2_L4.mom

2. Local axis definition H2O.axes

3. Polarizabilities defined in the local axis frame H2O_daTZ_ref_wt4_L1_static.pol

4. 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

1. 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.
2. 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
#include H2O_daTZ_ref_wt4_L1_static.pol
End
Limit rank 1 for O H1 H2
End

Next 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
#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
End

Switch 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.

AJMPublic/orient/energy-calculations (last edited 2021-04-14 12:56:33 by apw109)