geometry

Hessians, Dipole Moment Derivatives
and Harmonic Vibrational Frequencies

Harmonic force constants are computed via finite differences from analytic gradients in terms of internal coordinates. The procedure is semi-automatic since inputs for displaced geometries have to be created.
The dipole moment derivatives are calculated as finite differences also. In both cases a complete input set for gradient calculation has to be present in the main directory. In case dipole moment derivatives are computed, input files for dipole moment calculation have to be present.

Input generation

Before starting with the specific menu point in colinp for Harmonic force constants and dipole moment derivatives, a complete gradient input set has to be specified using the previous colinp menu points. You may choose MCSCF or CI/AQCC gradient.

Symmetry considerations: displaced molecular geometries created from non-totally symmetric coordinates will have lower symmetry than the original reference point of the optimized structure. Thus, new inputs matching the actual symmetry would have to be created for each case individually for optimal performance. The most convenient way is to ignore symmetry at all and create inputs in C₁ symmetry. For complete symmetry treatment see below.

Experience showed that calculations without using symmetry can give problems in the MCSCF for the construction of starting orbitals (e.g. from a SCF calculation). Moreover, one has to take care that the energy and the wavefunction show a continuous behavior on antisymmetric geometrical distortions (see reference symmetries). Therefore, the following scheme for the combination of geometry optimization and force constant calculation is recommended:

Hessian of the energy

optimize the geometry using a) symmetry, and in case of CI calculations b) interacting space restriction and c) one reference symmetry. This is the fastest alternative.
use symmetry, but no interacting space and all possible reference symmetries in the next step. The mocoef file from step 1) can be used. This calculation will give identical results to a calculation without symmetry and interacting space restriction. Optimize the molecule again starting from the geometry of step 1)
calculation without symmetry. Use the program transmo to transform the mocoef file of step 2) to a mocoef file without symmetry. Perform one more geometry optimization step to check that results are identical to step 2).
copy the input files from the main directory of step 3 into the main directory of the force constant calculation

Hessian of the nonadiabatic coupling

Prepare the input for a calculation of the nonadiabtic coupling terms in C₁ symmetry at the desired geometry.

Then you are ready to proceed. The following colinp menu point has to be entered:


     1) Job control for single point or gradient calculation
     2) Generate int. coordinates for potential energy curve
     3) Potential energy curve for one int. coordinate
-->  4) Vibrational frequencies and force constants
     5) Exit

------------------------------------------------------------

     submenu 4.1: freq.calculation options

     1) frequency calculation
     2) frequency and dipole moment derivative calculation

------------------------------------------------------------

     submenu 4.1.1/4.1.2: Individual int.coord. displacement adjustment:

   System with N internal coordinates detected!

    Int.coord. :   -disp [-0.001 ]  +disp [0.001  ]
    Int.coord. :   -disp [-0.001 ]  +disp [0.001  ]
      :                                  :
      :                                  :
      :                                  :
    Int.coord. :   -disp [-0.001 ]  +disp [0.001  ]


     calculate negative displacements only [n]

The last menu point generates a new directory structure:



     DISPLACEMENT/REFPOINT
     DISPLACEMENT/CALC.ca.d-A
     DISPLACEMENT/CALC.ca.dA
      :          
      :           
      :            
     DISPLACEMENT/CALC.cz.d-Z
     DISPLACEMENT/CALC.cz.dZ

In the DISPLACEMENT/REFPOINT directory is a copy of the unchanged input files for the reference geometry. The notation for the DISPLACEMENT/CALC.ca.d-A directories is: displacement -A in internal coordinate a. In every directory there will be now a new geometry file corresponding to the desired displacement. Moreover, the control file DISPLACEMENT/displfl is generated. Here the informations about points to be calculated are listed. It is also possible to calculate negative displacements only (reduces the computational effort by cca. 50%). In this case all positive displacements are omitted, and the numerical differentiation is performed relative to the reference point.

If only a subset of the forceconstant matrix is to be calculated, respective displacements have to be deactivated. This is done by adding the key word fixc at the end of the corresponding line in the file DISPLACEMENT/displfl.

Next steps are performed using the script: disp.pl. After calling the script disp.pl you will see on the screen:


        menu 1: COLUMBUS calculations based on the displfl file

     1) copy input files for force constant calculation
     2) copy input files for potential curve calculation
     3) perform force constant calculations (no CI restart)
     4) perform force constant calculations (using CI restart)
     5) perform potential curve calculations
     6) quit

The first point will copy all input files from the main directory into the subdirectories DISPLACEMENT/CALC.ca.dA (as listed in the file DISPLACEMENT/displfl). Thus, in these subdirectories complete inputs for individual gradient calculations for displaced geometries are located, which will be executed, as explained below. Instead of copying these same input files (usually for C₁ symmetry) into the subdirectories, specific input files for the calculation can be generated individually, in order to take into account symmetry and to achieve maximum performance. An input for geometry optimization has to be created with only one geometry optimization iteration. In the CI case the interacting space option should be disabled and all reference symmetries specified .

Executing the calculation

The calculation is started using the point 2 in the disp.pl script. The calculation is started in background with the nohup option and the points are calculated one by one. After one point is finished, the directories WORK and RESTART are removed. In this way the whole calculation requires only as much disk space as a single calculation.

Output analysis

After finishing the calculation using the script disp.pl we get in all displacement directories the corresponding gradients (nonadiabatic coupling terms) in internal coordinates (if required also dipole moment values). These data are collected automatically by the script forceconst.pl, which performs also the numerical differentiation. This program generates the file LISTINGS/forceconstls, which contains the basic informations about the performed calculations (calculation type, energy, convergence, number of performed iterations) and the calculated Hessian(s). The script generated also a file hessian in case of the Hessian of the energy and the files hessian_st1, hessian_st2 and hessian_nad in case of the Hessian of nonadiabatic couplings. These files are located in the main directory and may be used for geometry optimizations.

It is recommended to check the following informations listed in the LISTINGS/forceconstls file:

convergence of all the performed calculations
number of MCSCF iterations at every point (if the number of interations is much higher for some points than for the rest, then the character of the wave-function might have changed in this points and the ansatz for the wavefunction should be reconsidered.
number of CI/AQCC iterations. Especially when root following is used a sudden increase of the number of iterations may be the result of a calculation, where e.g. the root following did not work correctly. In this case the calculation has to be reconsidered (check the reference space, NVCIMN and NVCIMX values, ...).
energy differences relative to the reference point for all calculations. In case of shallow potential surfaces the size of the displacement has to be reconsidered (increased).
non-symmetricity in the force constant matrix. The harmonic force constant matrix should be symmetric. If large non-symmetric terms occurs, reconsider the displacement values or in the worst case the whole calculation.