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
- 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.
- The rotation (and translation) parameters for A and B are read in and the following steps performed:
- 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'.
- Calculate the integrals <math>S^{A'B'}_{k',l'}=(k'|l')</math>.
- Now calculate the block-diagonal rotation matrices <math>R^{AA'}</math> and <math>R^{BB'}</math>.
- 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>·
- Finallly: <math>(aa'|bb') = D^A_{aa',k} {\tilde{S}}^{AB}_{k,l} D^B_{bb',l} </math>.
- 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.