CamCASP/Programming/2

CamCASP => Programming => The DF-INTEGRAL module

Outline

This module consists of modules to calculate integrals needed for density-fitting, the density-fitting module, and modules to calculate integrals *using* primitive integrals in the auxiliary basis set and the density-fitting solution.

Through module df_int_operations and a lot of procedure overloading, a number of logical and computational operations are hidden from the higher-level energy subroutines. For example, to calculate the generalized OOOOAAAB-type 4-index Coulomb integral needed for the first-order exchange repulsion, subroutine calculate_e1exchS2 makes a call to open_gen_2eint('OOOO','AAAB',...) which does the following:

• DF: Performs the density-fitting for monomer A and dimer density-fitting. The integrals needed for the density-fitting will be automatically calculated.
• Integrals: The generalized integrals need various overlap, nuclear and 2-index Coulomb integrals. These will be identified and calculated.
• Assembly: All the components will be assembled into the required integrals.

These modules also need to decide on what needs re-calculation when molecules rotate. If the molecular basis alters, then everything needs to be re-calculated, but if all that happens is a rotation, then some molecular quantities need not be calculated again. In particular, the density-fitting solution remains invariant with respect to rotations. But rotations mix symmetry components of a shell. For example, (px,py,pz) would be mixed into (py,pz,-px) under a rotation by 90 degrees about the z-axis. This means that you cannot directly combine the density-fitting solution with primitive integrals calculated with the molecule in another orientation. In the next section we will describe what should be done.

Molecular rotations

First of all, as far as the dimer is concerned, a rotated molecule implies the dimer has changed. So *everything* needs to be recalculated. So let's look at monomer rotation only.

To see what's going on, consider the integral: <math>(aa'|bb')</math>. Using density-fitting (ignore integral robustness issues for now) we expand each pair of molecular orbitals as

<math>

 |aa') = \sum_k D^A_{aa',k} |k)

</math>

and

<math>

 |bb') = \sum_l D^B_{bb',l} |l)

</math>

Therefore

<math>

 (aa'|bb') = D^A_{aa',k} (k|l) D^B_{bb',l}.

</math>

Now, if the monomers are rotated by rotation matrices <math>R^A</math> and <math>R^B</math> the new 4-index integral is

<math>(\hat{R}(aa')|\hat{R'}(bb'))</math>.

Now, a rotation does two things: (1) rotates the *centers* of the basis functions, and (2) rotates the symmetry components of each shell. So,

<math>

 \hat{R}|aa') = \sum_k D^A_{aa',k} \hat{R}|k) = \sum_k D^A_{aa',k} {R^{A'A}_{k'k}}^{-1} |k'),

</math>

where <math>R^{AA'}_{kk'}</math> is a *block-diagonal* rotation matrix that does the rotation of shell components (it is block-diagonal as it does not mix components between shells) and the functions <math>|k')</math> are basis functions with their centers rotated. Notice that we have used the inverse of the rotation matrix.(explanation needed? later)

Therefore, our integral can be written as

<math>

 (\hat{R}(aa')|\hat{R'}(bb')) = D^A_{aa',k} {R^{A'A}_{k'k}}^{-1} (k'|l') {R^{B'B}_{l'l}}^{-1} D^B_{bb',l},

</math>

So we can now do the following:

• Calculate the matrix <math>S^{A'B'}_{k',l'}=(k'|l')</math>, where the primes indicate that the molecules are rotated.
• Now perform the block-rotations on <math>S^{A'B'}</math> and transform monomer A' back to A and B' back to B

<math>

 \tilde{S}^{AB}_{k,l} = {R^{A'A}_{k'k}}^{-1} S^{A'B'}_{k',l'} {R^{B'B}_{l'l}}^{-1}

</math> or <math>

 \tilde{S}^{AB}_{k,l} = R^{AA'}_{kk'} S^{A'B'}_{k',l'} R^{BB'}_{ll'}.

</math>

• We now construct the 4-index integral

<math>

 (aa'|bb') = D^A_{aa',k} \tilde{S}^{AB}_{k,l} D^B_{bb',l}.

</math>

So, our algorithm for calculating rotated integrals is

1. Perform the monomer density-fitting for A and B. This is done once in the molecular reference frame in which the molecular orbitals and Hessians have been obtained. This solution is recalculated only when the molecular geometry or basis alters.
2. The rotation (and translation) parameters for A and B are read in and the following steps performed:
   1. Rotation of the sites of molecules A and B. MOs are NOT rotated!!! After this step, we have the primed (i.e., rotated) molecules A' and B'.
   2. Calculate the integrals <math>S^{A'B'}_{k',l'}=(k'|l')</math>.
   3. Now calculate the block-diagonal rotation matrices <math>R^{AA'}</math> and <math>R^{BB'}</math>.
   4. Perform rotations of symmetries: <math>\tilde{S}^{AB}_{k,l} = R^{AA'}_{kk'} S^{A'B'}_{k',l'} R^{BB'}_{ll'}</math>. I.e., <math>{\tilde{S}}^{AB} = {R}^{AA'} {S}^{A'B'} ({R}^{BB'})^T </math>·
   5. Finallly: <math>(aa'|bb') = D^A_{aa',k} {\tilde{S}}^{AB}_{k,l} D^B_{bb',l} </math>.
3. Repeat the previous step as often as is required.

Implementation

So how is this done in practice? Let's look at the steps involved in the calculation of the density overlap integrals. These are defined as <math>

 (\rho^A|\rho^B) = \int \rho^A(r) \rho^B(r) dr.

</math>

About the code

The decisions to recalculate and/or rotate are made in subroutine check_if_integrals_are_done in df_integral_utilities.F90.

• If df_redo_df_on_rotation = .FALSE.:
  • (moA,moA) : no recalculation : The integrals involving just MOs of either A or B, but not both, do not need to be recalculated.
  • (moA,moA,auxA): no recalculation or rotation : The integrals involving MOs of either A or B, but not both, and additionally, auxiliary basis functions of the same monomer do not need to be recalculated or rotated. Note that, if they are recalculated, then they will need to have their symmetry components rotated to be consistent with the DF solution, but this would result in the very same integrals once again. So this is a waste!
  • (auxA,auxA): no recalculation or rotation : Integrals involving auxiliary basis functions of either A or B (not both) do not need recalculation or rotation.
  • (moA,moB): recalculation : Integrals involving MOs of both A and B will need to be recalculated.
  • (auxA,auxB): recalculation and rotation : Integrals involving auxiliary basis functions of both A and B (but not AB) will need to be recalculated and rotated. This category also includes integrals like NUCA,__XB that uses the sites of one monomer and the auxiliary basis functions of the other.
  • (auxAB,auxA): recalculation and partial rotation : Integrals involving auxiliary basis functions of the dimer AB will always need to be recalculated, but the rotation is necessary only if one of the indices involves auxiliary basis functions of either A or B. So <math>S^{A,AB}</math> needs recalculation and needs to have the symmetry components of the left index rotated.
• If df_redo_df_on_rotation = .TRUE.:
  • (moA,moA): no recalculation.
  • (moA,moA,auxA): recalculation, but no rotation.
  • (auxA,auxA): recalculation, but no rotation.
  • (moA,moB): recalculation.
  • (auxA,auxB): recalculation, but no rotation.
  • (auxAB,auxA): recalculation, but no rotation.

The advantage of using df_redo_df_on_rotation = .FALSE. is clear. A lot of recalculations are avoided. The biggest saving is the DF is not re-done for each monomer geometry.

Integral types and rotation

Every integral type has the logical fields ForceRotate and Rotate. The former is set to FALSE by default. It can be set to TRUE on an integral-by-integral basis using the BEGIN DF-INTS...END block. The latter, I%Rotate, is set for every integral when they are defined. The above rules are used to decide the value of this field.

Integrals will have their symmetry components rotated IF:

1. I%ForceRotate = .TRUE., or
2. I%Rotate = .TRUE. && REDO_DF_ON_ROTATION = .FALSE..

Integrals that should be rotated

This is a list of integrals that should be rotated if REDO_DF_ON_ROTATION = .FALSE.:

1. OVRL,XAXB
2. KE__,XAXB
3. NUCA,__XB
4. NUCB,__XA
5. PROB,__XA
6. PROB,__XB
7. S___,_A_B
8. SC__,_A_B
9. S___,A_AB *left-rotation only
10. S___,B_AB *left-rotation only

Special Cases

Some integrals depend on parameters like the probe charge (PROB), regularization parameter (NUCA,NUCB,PROB) and DF constraints (SC__). So if these parameters change, these integrals will be re-calculated. Now, some of these, for example, NUCA,__XA, have %Rotate=.false., and so are not rotated when calculated. This can lead to an error. Consider the follow sequence of steps:

1. REDO_DF_ON_ROTATION=.false.
2. DF performed.
3. NUCA,__XA calculated. It is consistent with the DF solution.
4. Monomer A rotated.
5. NUCA,__XA not rotated because %Rotate=.false. Still consistent with DF.
6. reg_eta changed.
7. NUCA,__XA is now re-calculated using the rotated monomer A, but symmetry blocks are not rotated because %Rotate=.false.
8. It is now inconsistent with DF solution.

It is unlikely that this will happen as this is not normal usage of the code, but it could and it would be awful to track down. So, I need a work-around.

Effects integrals:

1. NUCA,__XA
2. NUCB,__XB
3. SC__,_A_A
4. SC__,_B_B

One solution is to use the %ForceRotate field to force rotation when this happens. This field will need to be reset after rotation is performed or else a rotation will be performed on a subsequent call even if nothing has changed.

A cleaner solution might be to set %Rotate=.true. for these integrals and re-calculate them when the molecules rotate. Then the RECALCULATE-ROTATE operations will always be performed on these integrals. This will be a waste of CPU most of the time.

NEEDS RESOLVING