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>