Diff for "AJMPublic/orient/energy-calculations"

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

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
   Read rank 1
     #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
   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
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%.