<> == Orient Energy Calculations == Main page [[AJMPublic/orient|Orient]] See also '''Orient''' manual found at [[http://www-stone.ch.cam.ac.uk/programs.html|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. 1. One or more dimer configurations at which interaction energies are to be evaluated. 1. 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. 1. Distributed polarizabilities. In either the laboratory axis system or a set of local axis frames. 1. If local axes are used, these need to be specified in an Orient AXIS file. 1. A damping model. Specify damping function type and damping parameters. 1. Iterations used and convergence criterion. Dispersion energy 1. Dispersion model. This can be part of the potential specification. 1. 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. 1. A damping model. Short-range repulsion 1. Parameters specified in a potential file. 1. 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 [[AJMPublic/orient/examples/water2-orient-energy-calculation|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 [[AJMPublic/orient/examples/water-moments|H2O_daTZ_DMA2_L4.mom]] 1. Local axis definition [[AJMPublic/orient/examples/water-localaxes|H2O.axes]] 1. Polarizabilities defined in the local axis frame [[AJMPublic/orient/examples/water-WSML1-pol|H2O_daTZ_ref_wt4_L1_static.pol]] 1. Dimer geometry file [[AJMPublic/orient/examples/water-globalminconfig|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 '''TYPE'''s 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. 1. 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 [[AJMPublic/orient/examples/water-moments|H2O_daTZ_DMA2_L4.mom]] which is '''#include'''d 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 [[AJMPublic/orient/examples/water-localaxes|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 [[AJMPublic/orient/examples/water-WSML1-pol|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 ''. 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 [[ajmorientexampleswater-globalminconfig|globalminconfig.geom]] which is '''#include'''d. 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 [[AJMPublic/orient/examples/water2-orient-energy-calculation|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 [[AJMPublic/orient/examples/water2-orient-energy-calculation-output|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%. {{{#!wiki note Note See Stone's book, chapter 10 in the second edition for more on many-body effects and the polarization.}}}