37 KiB
Diagonal hessian
The hessian of the CI energy with respects to the orbital rotation is : (C-c C-x C-l)
\begin{align*} H_{pq,rs} &= \dfrac{\partial^2 E(x)}{\partial x_{pq}^2} \\ &= \mathcal{P}_{pq} \mathcal{P}_{rs} [ \frac{1}{2} \sum_u [\delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) + \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_r^u)] -(h_p^s \gamma_r^q + h_r^q \gamma_p^s) \\ &+ \frac{1}{2} \sum_{tuv} [\delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} + v_{uv}^{st} \Gamma_{pt}^{uv}) + \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})] \\ &+ \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs} \Gamma_{pr}^{uv}) - \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts}) \end{align*}With pq a permutation operator :
\begin{align*} \mathcal{P}_{pq}= 1 - (p \leftrightarrow q) \end{align*} \begin{align*} \mathcal{P}_{pq} \mathcal{P}_{rs} &= (1 - (p \leftrightarrow q))(1 - (r \leftrightarrow s)) \\ &= 1 - (p \leftrightarrow q) - (r \leftrightarrow s) + (p \leftrightarrow q, r \leftrightarrow s) \end{align*}Where p,q,r,s,t,u,v are general spatial orbitals mo_num : the number of molecular orbitals $$h$$ : One electron integrals $$\gamma$$ : One body density matrix (state average in our case) $$v$$ : Two electron integrals $$\Gamma$$ : Two body density matrice (state average in our case)
Source : Seniority-based coupled cluster theory
- Chem. Phys. 141, 244104 (2014); https://doi.org/10.1063/1.4904384
Thomas M. Henderson, Ireneusz W. Bulik, Tamar Stein, and Gustavo E. Scuseria
Here for the diagonal of the hessian it's a little more complicated than for the hessian. It's not just compute the diagonal terms of the hessian because of the permutations.
The hessian is (p,q,r,s), so the diagonal terms are (p,q,p,q). But with the permutations : p <-> q, r <-> s, p <-> q and r <-> s, we have a diagonal term, if : p = r and q = s, => (p,q,p,q) or q = r and p = s, => (p,q,q,p)
For that reason, we will use 2D temporary arrays to store the elements. One for the terms (p,q,p,q) and an other for the terms of kind (p,q,q,p). We will also use a 1D temporary array to store the terms of the kind (p,p,p,p) due to the kronoecker delta.
Compute the diagonal hessian of energy with respects to orbital rotations By diagonal hessian we mean, diagonal elements of the hessian
Provided:
| mo_num | integer | number of MOs |
| mo_one_e_integrals(mo_num,mo_num) | double precision | mono-electronic integrals |
| one_e_dm_mo(mo_num,mo_num) | double precision | one e- density matrix (state average) |
| two_e_dm_mo(mo_num,mo_num,mo_num) | double precision | two e- density matrix (state average) |
Input:
| n | integer | mo_num*(mo_num-1)/2 |
Output:
| H(n,n) | double precision | Hessian matrix |
| h_tmpr(mo_num,mo_num,mo_num,mo_num) | double precision | Complete hessian matrix before the tranformation |
| in n by n matrix |
Internal:
| hessian(mo_num,mo_num,mo_num,mo_num) | double precision | temporary array containing the hessian before |
| the permutations | ||
| p, q, r, s | integer | indexes of the hessian elements |
| t, u, v | integer | indexes for the sums |
| pq, rs | integer | indexes for the transformation of the hessian |
| (4D -> 2D) | ||
| t1,t2,t3 | double precision | time to compute the hessian |
| t4,t5,t6 | double precision | time to compute the differ each element |
| tmp_bi_int_3(mo_num,mo_num,mo_num) | double precision | 3 indexes temporary array for the bielectronic integrals (private) |
| tmp_bi_int_3_shared(mo_num,mo_num,mo_num) | double precision | 3 indexes temporary array for the bielectronic integrals (shared) |
| tmp_2rdm_3(mo_num,mo_num,mo_num) | double precision | 3 indexes temporary array for the 2 body density matrix (private) |
| tmp_2rdm_3_shared(mo_num,mo_num,mo_num) | double precision | 3 indexes temporary array for the 2 body density matrix (shared) |
| tmp_accu(mo_num,mo_num) | double precision | temporary array (private) |
| tmp_accu_shared(mo_num,mo_num) | double precision | temporary array (shared) |
| tmp_accu_1(mo_num) | double precision | temporary array (private) |
| tmp_accu_1_shared(mo_num) | double precision | temporary array (shared) |
| tmp_h_pppp(mo_num) | double precision | matrix containing the hessien elements hessian(p,p,p,p) |
| tmp_h_pqpq(mo_num,mo_num) | double precision | matrix containing the hessien elements hessian(p,q,p,q) |
| tmp_h_pqqp(mo_num,mo_num) | double precision | matrix containing the hessien elements hessian(p,q,q,p) |
Function:
| get_two_e_integral | double precision | bi-electronic integrals |
subroutine diag_hessian_opt(n,H)!, h_tmpr)
use omp_lib
include 'constants.h'
implicit none
! Variables
! in
integer, intent(in) :: n
! out
double precision, intent(out) :: H(n)!,n), h_tmpr(mo_num,mo_num,mo_num,mo_num)
! internal
!double precision, allocatable :: hessian(:,:,:,:)!, h_tmpr(:,:,:,:)
integer :: p,q,k
integer :: r,s,t,u,v
integer :: pq,rs
integer :: istate
double precision :: t1,t2,t3,t4,t5,t6
double precision, allocatable :: tmp_bi_int_3(:,:,:),tmp_bi_int_3_shared(:,:,:)
double precision, allocatable :: tmp_2rdm_3(:,:,:),tmp_2rdm_3_shared(:,:,:)
double precision, allocatable :: tmp_accu(:,:)
double precision, allocatable :: tmp_accu_shared(:,:), tmp_accu_1_shared(:)
double precision, allocatable :: tmp_h_pppp(:), tmp_h_pqpq(:,:), tmp_h_pqqp(:,:)
! Function
double precision :: get_two_e_integral
print*,''
print*,'--- diagonal hessian---'
print*,'Use the diagonal hessian'
! Allocation of shared arrays
!allocate(hessian(mo_num,mo_num,mo_num,mo_num))!,h_tmpr(mo_num,mo_num,mo_num,mo_num))
allocate(tmp_h_pppp(mo_num),tmp_h_pqpq(mo_num,mo_num),tmp_h_pqqp(mo_num,mo_num))
allocate(tmp_2rdm_3_shared(mo_num,mo_num,mo_num))
allocate(tmp_bi_int_3_shared(mo_num,mo_num,mo_num))
allocate(tmp_accu_1_shared(mo_num),tmp_accu_shared(mo_num,mo_num))
! OMP
call omp_set_max_active_levels(1)
!$OMP PARALLEL &
!$OMP PRIVATE( &
!$OMP p,q,r,s, tmp_accu,k, &
!$OMP u,v,t, tmp_bi_int_3, tmp_2rdm_3) &
!$OMP SHARED(H, tmp_h_pppp, tmp_h_pqpq, tmp_h_pqqp, &
!$OMP mo_num,n, mo_one_e_integrals, one_e_dm_mo, &
!$OMP tmp_bi_int_3_shared, tmp_2rdm_3_shared,tmp_accu_shared, &
!$OMP tmp_accu_1_shared,two_e_dm_mo,mo_integrals_map,t1,t2,t3,t4,t5,t6) &
!$OMP DEFAULT(NONE)
! Allocation of the private arrays
allocate(tmp_2rdm_3(mo_num,mo_num,mo_num),tmp_bi_int_3(mo_num,mo_num,mo_num))
allocate(tmp_accu(mo_num,mo_num))Initialization of the arrays
!!$OMP DO
!do s = 1,mo_num
! do r = 1, mo_num
! do q = 1, mo_num
! do p = 1, mo_num
! hessian(p,q,r,s) = 0d0
! enddo
! enddo
! enddo
!enddo
!!$OMP END DO
!$OMP DO
do p = 1, mo_num
tmp_h_pppp(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqpq(p,q) = 0d0
enddo
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = 0d0
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t1)
!$OMP END MASTERLine 1, term 1
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s))) then
if (q==r) then do u = 1, mo_num
hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( & mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
- mo_one_e_integrals(s,u) * one_e_dm_mo(p,u))
enddo endif endif enddo enddo enddo enddo
With optimization :
Part 1 : p=r and q=s and q=r
hessian(p,q,r,s) -> hessian(p,p,p,p)
0.5d0 * ( & mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
- mo_one_e_integrals(s,u) * one_e_dm_mo(p,u)) = 0.5d0 * ( & mo_one_e_integrals(u,p) * one_e_dm_mo(u,p) &
- mo_one_e_integrals(p,u) * one_e_dm_mo(p,u)) = mo_one_e_integrals(u,p) * one_e_dm_mo(u,p)
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(p) = tmp_accu_1_shared(p) + mo_one_e_integrals(u,p) * one_e_dm_mo(u,p)
enddo
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
tmp_h_pppp(p) = tmp_h_pppp(p) + tmp_accu_1_shared(p)
enddo
!$OMP END DOPart 2 : q=r and p=s and q=r
hessian(p,q,r,s) -> hessian(p,q,q,p)
0.5d0 * ( & mo_one_e_integrals(u,p) * one_e_dm_mo(u,s) &
- mo_one_e_integrals(s,u) * one_e_dm_mo(p,u)) = 0.5d0 * ( & mo_one_e_integrals(u,p) * one_e_dm_mo(u,p) &
- mo_one_e_integrals(p,u) * one_e_dm_mo(p,u)) = mo_one_e_integrals(u,p) * one_e_dm_mo(u,p)
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(p) = tmp_accu_1_shared(p) + mo_one_e_integrals(u,p) * one_e_dm_mo(u,p)
enddo
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) + tmp_accu_1_shared(p)
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6= t5-t4
print*,'l1 1',t6
!$OMP END MASTERLine 1, term 2
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s))) then
if (p==s) then do u = 1, mo_num
hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( & mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
- mo_one_e_integrals(q,u) * one_e_dm_mo(r,u))
enddo endif endif enddo enddo enddo enddo
Part 1 : p=r and q=s and p=s
hessian(p,q,r,s) -> hessian(p,p,p,p)
0.5d0 * (& mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
- mo_one_e_integrals(q,u) * one_e_dm_mo(r,u)) = 0.5d0 * ( & mo_one_e_integrals(u,p) * one_e_dm_mo(u,p) &
- mo_one_e_integrals(p,u) * one_e_dm_mo(p,u)) = mo_one_e_integrals(u,p) * one_e_dm_mo(u,p)
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(p) = tmp_accu_1_shared(p) + mo_one_e_integrals(u,p) * one_e_dm_mo(u,p)
enddo
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
tmp_h_pppp(p) = tmp_h_pppp(p) + tmp_accu_1_shared(p)
enddo
!$OMP END DOPart 2 : q=r and p=s and p=s
hessian(p,q,r,s) -> hessian(p,q,q,p)
0.5d0 * (& mo_one_e_integrals(u,r) * one_e_dm_mo(u,q) &
- mo_one_e_integrals(q,u) * one_e_dm_mo(r,u)) = 0.5d0 * ( & mo_one_e_integrals(u,q) * one_e_dm_mo(u,q) &
- mo_one_e_integrals(q,u) * one_e_dm_mo(q,u)) = mo_one_e_integrals(u,q) * one_e_dm_mo(u,q)
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(q) = tmp_accu_1_shared(q) + mo_one_e_integrals(u,q) * one_e_dm_mo(u,q)
enddo
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) + tmp_accu_1_shared(q)
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6= t5-t4
print*,'l1 2',t6
!$OMP END MASTERLine 1, term 3
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s))) then
hessian(p,q,r,s) = hessian(p,q,r,s) &
- mo_one_e_integrals(s,p) * one_e_rdm_mo(r,q) &
- mo_one_e_integrals(q,r) * one_e_rdm_mo(p,s)
endif enddo enddo enddo enddo
With optimization :
Part 1 : p=r and q=s
hessian(p,q,r,s) -> hessian(p,q,p,q)
- mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
- mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
=
- mo_one_e_integrals(q,p) * one_e_dm_mo(p,q) &
- mo_one_e_integrals(q,p) * one_e_dm_mo(p,q)
=
- 2d0 mo_one_e_integrals(q,p) * one_e_dm_mo(p,q)
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqpq(p,q) = tmp_h_pqpq(p,q) &
- 2d0 * mo_one_e_integrals(q,p) * one_e_dm_mo(p,q)
enddo
enddo
!$OMP END DOPart 2 : q=r and p=s
hessian(p,q,r,s) -> hessian(p,q,p,q)
- mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
- mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
=
- mo_one_e_integrals(q,p) * one_e_dm_mo(p,q) &
- mo_one_e_integrals(q,p) * one_e_dm_mo(p,q)
=
- 2d0 mo_one_e_integrals(q,p) * one_e_dm_mo(p,q)
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) &
- 2d0 * mo_one_e_integrals(p,p) * one_e_dm_mo(q,q)
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6= t5-t4
print*,'l1 3',t6
!$OMP END MASTERLine 2, term 1
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s))) then
if (q==r) then do t = 1, mo_num do u = 1, mo_num do v = 1, mo_num
hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( & get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
- get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v))
enddo enddo enddo endif endif enddo enddo enddo enddo
With optimization :
Part 1 : p=r and q=s and q=r
hessian(p,q,r,s) -> hessian(p,p,p,p)
0.5d0 * ( & get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
- get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v)) = 0.5d0 * ( & get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,p,t) &
- get_two_e_integral(p,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v)) = get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,p,t)
Just re-order the index and use 3D temporary arrays for optimal memory accesses.
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do t = 1, mo_num
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_bi_int_3(u,v,p) = get_two_e_integral(u,v,p,t,mo_integrals_map)
enddo
enddo
enddo
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3(u,v,p) = two_e_dm_mo(u,v,p,t)
enddo
enddo
enddo
!$OMP CRITICAL
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(p) = tmp_accu_1_shared(p) &
+ tmp_bi_int_3(u,v,p) * tmp_2rdm_3(u,v,p)
enddo
enddo
enddo
!$OMP END CRITICAL
enddo
!$OMP END DO
!$OMP DO
do p =1, mo_num
tmp_h_pppp(p) = tmp_h_pppp(p) + tmp_accu_1_shared(p)
enddo
!$OMP END DOPart 2 : q=r and p=s and q=r
hessian(p,q,r,s) -> hessian(p,q,q,p)
0.5d0 * ( & get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,s,t) &
- get_two_e_integral(s,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v)) = 0.5d0 * ( & get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,p,t) &
- get_two_e_integral(p,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v)) = get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,p,t)
Just re-order the index and use 3D temporary arrays for optimal memory accesses.
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do t = 1, mo_num
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_bi_int_3(u,v,p) = get_two_e_integral(u,v,p,t,mo_integrals_map)
enddo
enddo
enddo
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3(u,v,p) = two_e_dm_mo(u,v,p,t)
enddo
enddo
enddo
!$OMP CRITICAL
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(p) = tmp_accu_1_shared(p) + &
tmp_bi_int_3(u,v,p) * tmp_2rdm_3(u,v,p)
enddo
enddo
enddo
!$OMP END CRITICAL
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) + tmp_accu_1_shared(p)
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6 = t5-t4
print*,'l2 1',t6
!$OMP END MASTERLine 2, term 2
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s))) then
if (p==s) then do t = 1, mo_num do u = 1, mo_num do v = 1, mo_num
hessian(p,q,r,s) = hessian(p,q,r,s) + 0.5d0 * ( & get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
- get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t))
enddo enddo enddo endif endif enddo enddo enddo enddo
With optimization :
Part 1 : p=r and q=s and p=s
hessian(p,q,r,s) -> hessian(p,p,p,p)
0.5d0 * ( & get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
- get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t)) = 0.5d0 * ( & get_two_e_integral(p,t,u,v,mo_integrals_map) * two_e_dm_mo(p,t,u,v) &
- get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,p,t)) = get_two_e_integral(u,v,p,t,mo_integrals_map) * two_e_dm_mo(u,v,p,t)
Just re-order the index and use 3D temporary arrays for optimal memory accesses.
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do p = 1, mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do t = 1, mo_num
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_bi_int_3(u,v,p) = get_two_e_integral(u,v,p,t,mo_integrals_map)
enddo
enddo
enddo
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3(u,v,p) = two_e_dm_mo(u,v,p,t)
enddo
enddo
enddo
!$OMP CRITICAL
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(p) = tmp_accu_1_shared(p) +&
tmp_bi_int_3(u,v,p) * tmp_2rdm_3(u,v,p)
enddo
enddo
enddo
!$OMP END CRITICAL
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
tmp_h_pppp(p) = tmp_h_pppp(p) + tmp_accu_1_shared(p)
enddo
!$OMP END DOPart 2 : q=r and p=s and p=s
hessian(p,q,r,s) -> hessian(p,q,q,p)
0.5d0 * ( & get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(r,t,u,v) &
- get_two_e_integral(u,v,r,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t)) = 0.5d0 * ( & get_two_e_integral(q,t,u,v,mo_integrals_map) * two_e_dm_mo(q,t,u,v) &
- get_two_e_integral(u,v,q,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t)) = get_two_e_integral(u,v,q,t,mo_integrals_map) * two_e_dm_mo(u,v,q,t)
Just re-order the index and use 3D temporary arrays for optimal memory accesses.
!$OMP DO
do p = 1,mo_num
tmp_accu_1_shared(p) = 0d0
enddo
!$OMP END DO
!$OMP DO
do t = 1, mo_num
do q = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_bi_int_3(u,v,q) = get_two_e_integral(u,v,q,t,mo_integrals_map)
enddo
enddo
enddo
do q = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3(u,v,q) = two_e_dm_mo(u,v,q,t)
enddo
enddo
enddo
!$OMP CRITICAL
do q = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_accu_1_shared(q) = tmp_accu_1_shared(q) +&
tmp_bi_int_3(u,v,q) * tmp_2rdm_3(u,v,q)
enddo
enddo
enddo
!$OMP END CRITICAL
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) + tmp_accu_1_shared(p)
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6 = t5-t4
print*,'l2 2',t6
!$OMP END MASTERLine 3, term 1
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)))) then
do u = 1, mo_num do v = 1, mo_num
hessian(p,q,r,s) = hessian(p,q,r,s) &
- get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
- get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v)
enddo enddo endif
enddo enddo enddo enddo
With optimization
Part 1 : p=r and q=s
hessian(p,q,r,s) -> hessian(p,q,p,q)
get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
- get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v) = get_two_e_integral(u,v,p,p,mo_integrals_map) * two_e_dm_mo(u,v,q,q) &
- get_two_e_integral(q,q,u,v,mo_integrals_map) * two_e_dm_mo(p,p,u,v) = 2d0 * get_two_e_integral(u,v,p,p,mo_integrals_map) * two_e_dm_mo(u,v,q,q)
Arrays of the kind (u,v,p,p) can be transform in 4D arrays (u,v,p). Using u,v as one variable a matrix multiplication appears. $$c_{p,q} = \sum_{uv} a_{p,uv} b_{uv,q}$$
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do q = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3_shared(u,v,q) = two_e_dm_mo(u,v,q,q)
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_bi_int_3_shared(u,v,p) = get_two_e_integral(u,v,p,p,mo_integrals_map)
enddo
enddo
enddo
!$OMP END DO
call dgemm('T','N', mo_num, mo_num, mo_num*mo_num, 1d0, tmp_bi_int_3_shared,&
mo_num*mo_num, tmp_2rdm_3_shared, mo_num*mo_num, 0d0, tmp_accu, mo_num)
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqpq(p,q) = tmp_h_pqpq(p,q) + tmp_accu(p,q) + tmp_accu(q,p)
enddo
enddo
!$OMP END DOPart 2 : q=r and p=s
hessian(p,q,r,s) -> hessian(p,q,q,p)
get_two_e_integral(u,v,p,r,mo_integrals_map) * two_e_dm_mo(u,v,q,s) &
- get_two_e_integral(q,s,u,v,mo_integrals_map) * two_e_dm_mo(p,r,u,v) = get_two_e_integral(u,v,p,q,mo_integrals_map) * two_e_dm_mo(u,v,q,p) &
- get_two_e_integral(q,p,u,v,mo_integrals_map) * two_e_dm_mo(p,q,u,v) = 2d0 * get_two_e_integral(u,v,p,q,mo_integrals_map) * two_e_dm_mo(u,v,q,p)
Just re-order the indexes and use 3D temporary arrays for optimal memory accesses.
!$OMP MASTER
call wall_time(t4)
!$OMP END MASTER
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_bi_int_3(u,v,p) = 2d0 * get_two_e_integral(u,v,q,p,mo_integrals_map)
enddo
enddo
enddo
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3(u,v,p) = two_e_dm_mo(u,v,p,q)
enddo
enddo
enddo
do p = 1, mo_num
do v = 1, mo_num
do u = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) &
+ tmp_bi_int_3(u,v,p) * tmp_2rdm_3(u,v,p)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6= t5-t4
print*,'l3 1',t6
!$OMP END MASTERLine 3, term 2
Without optimization :
do p = 1, mo_num do q = 1, mo_num do r = 1, mo_num do s = 1, mo_num
! Permutations if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) & .or. ((p==s) .and. (q==r))) then
do t = 1, mo_num do u = 1, mo_num
hessian(p,q,r,s) = hessian(p,q,r,s) &
- get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
- get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
- get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
- get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
enddo enddo
endif
enddo enddo enddo enddo
With optimization :
Part 1 : p=r and q=s
hessian(p,q,r,s) -> hessian(p,q,p,q)
- get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
- get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
- get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
- get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
=
- get_two_e_integral(q,t,p,u,mo_integrals_map) * two_e_dm_mo(p,t,q,u) &
- get_two_e_integral(t,q,p,u,mo_integrals_map) * two_e_dm_mo(t,p,q,u) &
- get_two_e_integral(q,u,p,t,mo_integrals_map) * two_e_dm_mo(p,u,q,t) &
- get_two_e_integral(q,u,t,p,mo_integrals_map) * two_e_dm_mo(p,u,t,q)
=
- 2d0 * get_two_e_integral(q,t,p,u,mo_integrals_map) * two_e_dm_mo(p,t,q,u) &
- 2d0 * get_two_e_integral(t,q,p,u,mo_integrals_map) * two_e_dm_mo(t,p,q,u)
=
- 2d0 * get_two_e_integral(q,u,p,t,mo_integrals_map) * two_e_dm_mo(q,u,p,t) &
- 2d0 * get_two_e_integral(t,q,p,u,mo_integrals_map) * two_e_dm_mo(t,p,q,u)
Just re-order the indexes and use 3D temporary arrays for optimal memory accesses.
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!----------
! Part 1.1
!----------
! - 2d0 * get_two_e_integral(q,u,p,t,mo_integrals_map) * two_e_dm_mo(q,u,p,t)
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_accu_shared(p,q) = 0d0
enddo
enddo
!$OMP END DO
!$OMP DO
do t = 1, mo_num
do p = 1, mo_num
do u = 1, mo_num
do q = 1, mo_num
tmp_bi_int_3(q,u,p) = 2d0 * get_two_e_integral(q,u,p,t,mo_integrals_map)
enddo
enddo
enddo
do p = 1, mo_num
do u = 1, mo_num
do q = 1, mo_num
tmp_2rdm_3(q,u,p) = two_e_dm_mo(q,u,p,t)
enddo
enddo
enddo
!$OMP CRITICAL
do p = 1, mo_num
do u = 1, mo_num
do q = 1, mo_num
tmp_accu_shared(p,q) = tmp_accu_shared(p,q) &
- tmp_bi_int_3(q,u,p) * tmp_2rdm_3(q,u,p)
enddo
enddo
enddo
!$OMP END CRITICAL
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqpq(p,q) = tmp_h_pqpq(p,q) + tmp_accu_shared(p,q)
enddo
enddo
!$OMP END DOJust re-order the indexes and use 3D temporary arrays for optimal memory accesses.
!--------
! Part 1.2
!--------
! - 2d0 * get_two_e_integral(t,q,p,u,mo_integrals_map) * two_e_dm_mo(t,p,q,u)
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_accu_shared(p,q) = 0d0
enddo
enddo
!$OMP END DO
!$OMP DO
do u = 1, mo_num
do p = 1, mo_num
do q = 1, mo_num
do t = 1, mo_num
tmp_bi_int_3(t,q,p) = 2d0*get_two_e_integral(t,q,p,u,mo_integrals_map)
enddo
enddo
enddo
do p= 1, mo_num
do q = 1, mo_num
do t = 1, mo_num
tmp_2rdm_3(t,q,p) = two_e_dm_mo(t,p,q,u)
enddo
enddo
enddo
!$OMP CRITICAL
do q = 1, mo_num
do p = 1, mo_num
do t = 1, mo_num
tmp_accu_shared(p,q) = tmp_accu_shared(p,q) &
- tmp_bi_int_3(t,q,p) * tmp_2rdm_3(t,q,p)
enddo
enddo
enddo
!$OMP END CRITICAL
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqpq(p,q) = tmp_h_pqpq(p,q) + tmp_accu_shared(p,q)
enddo
enddo
!$OMP END DOPart 2 : q=r and p=s
hessian(p,q,r,s) -> hessian(p,q,p,q)
- get_two_e_integral(s,t,p,u,mo_integrals_map) * two_e_dm_mo(r,t,q,u) &
- get_two_e_integral(t,s,p,u,mo_integrals_map) * two_e_dm_mo(t,r,q,u) &
- get_two_e_integral(q,u,r,t,mo_integrals_map) * two_e_dm_mo(p,u,s,t) &
- get_two_e_integral(q,u,t,r,mo_integrals_map) * two_e_dm_mo(p,u,t,s)
=
- get_two_e_integral(p,t,p,u,mo_integrals_map) * two_e_dm_mo(q,t,q,u) &
- get_two_e_integral(t,p,p,u,mo_integrals_map) * two_e_dm_mo(t,q,q,u) &
- get_two_e_integral(q,u,q,t,mo_integrals_map) * two_e_dm_mo(p,u,p,t) &
- get_two_e_integral(q,u,t,q,mo_integrals_map) * two_e_dm_mo(p,u,t,p)
=
- get_two_e_integral(p,t,p,u,mo_integrals_map) * two_e_dm_mo(q,t,q,u) &
- get_two_e_integral(q,t,q,u,mo_integrals_map) * two_e_dm_mo(p,t,p,u) &
- get_two_e_integral(t,u,p,p,mo_integrals_map) * two_e_dm_mo(t,q,q,u) &
- get_two_e_integral(t,u,q,q,mo_integrals_map) * two_e_dm_mo(t,p,p,u)
=
- get_two_e_integral(t,p,u,p,mo_integrals_map) * two_e_dm_mo(t,q,u,q) &
- get_two_e_integral(t,q,u,q,mo_integrals_map) * two_e_dm_mo(p,t,p,u) &
- get_two_e_integral(t,u,p,p,mo_integrals_map) * two_e_dm_mo(q,u,t,q) &
- get_two_e_integral(t,u,q,q,mo_integrals_map) * two_e_dm_mo(p,u,t,p)
Arrays of the kind (t,p,u,p) can be transformed in 3D arrays. By doing so and using t,u as one variable, a matrix multiplication appears : $$c_{p,q} = \sum_{tu} a_{p,tu} b_{tu,q}$$
!----------
! Part 2.1
!----------
! - get_two_e_integral(t,p,u,p,mo_integrals_map) * two_e_dm_mo(t,q,u,q) &
! - get_two_e_integral(t,q,u,q,mo_integrals_map) * two_e_dm_mo(p,t,p,u)
!$OMP DO
do q = 1, mo_num
do u = 1, mo_num
do t = 1, mo_num
tmp_2rdm_3_shared(t,u,q) = two_e_dm_mo(t,q,u,q)
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do p = 1, mo_num
do u = 1, mo_num
do t = 1, mo_num
tmp_bi_int_3_shared(t,u,p) = get_two_e_integral(t,p,u,p,mo_integrals_map)
enddo
enddo
enddo
!$OMP END DO
call dgemm('T','N', mo_num, mo_num, mo_num*mo_num, 1d0, tmp_bi_int_3_shared,&
mo_num*mo_num, tmp_2rdm_3_shared, mo_num*mo_num, 0d0, tmp_accu, mo_num)
!$OMP DO
do p = 1, mo_num
do q = 1, mo_num
tmp_h_pqqp(q,p) = tmp_h_pqqp(q,p) - tmp_accu(q,p) - tmp_accu(p,q)
enddo
enddo
!$OMP END DOArrays of the kind (t,u,p,p) can be transformed in 3D arrays. By doing so and using t,u as one variable, a matrix multiplication appears : $$c_{p,q} = \sum_{tu} a_{p,tu} b_{tu,q}$$
!--------
! Part 2.2
!--------
! - get_two_e_integral(t,u,p,p,mo_integrals_map) * two_e_dm_mo(q,u,t,q) &
! - get_two_e_integral(t,u,q,q,mo_integrals_map) * two_e_dm_mo(p,u,t,p)
!$OMP DO
do p = 1, mo_num
do u = 1, mo_num
do t = 1, mo_num
tmp_bi_int_3_shared(t,u,p) = get_two_e_integral(t,u,p,p,mo_integrals_map)
enddo
enddo
enddo
!$OMP END DO
!$OMP DO
do q = 1, mo_num
do t = 1, mo_num
do u = 1, mo_num
tmp_2rdm_3_shared(u,t,q) = two_e_dm_mo(q,u,t,q)
enddo
enddo
enddo
!$OMP END DO
call dgemm('T','N', mo_num, mo_num, mo_num*mo_num, 1d0, tmp_2rdm_3_shared,&
mo_num*mo_num, tmp_bi_int_3_shared, mo_num*mo_num, 0d0, tmp_accu, mo_num)
!$OMP DO
do q = 1, mo_num
do p = 1, mo_num
tmp_h_pqqp(p,q) = tmp_h_pqqp(p,q) - tmp_accu(p,q) - tmp_accu(q,p)
enddo
enddo
!$OMP END DO
!$OMP MASTER
CALL wall_TIME(t5)
t6= t5-t4
print*,'l3 2',t6
!$OMP END MASTER
!$OMP MASTER
CALL wall_TIME(t2)
t2 = t2 - t1
print*, 'Time to compute the hessian :', t2
!$OMP END MASTERDeallocation of private arrays
In the OMP section !
deallocate(tmp_2rdm_3,tmp_bi_int_3)
deallocate(tmp_accu)Permutations
As we mentioned before there are two permutation operator in the formula : Hessian(p,q,r,s) = P_pq P_rs […] => Hessian(p,q,r,s) = (p,q,r,s) - (q,p,r,s) - (p,q,s,r) + (q,p,s,r)
!!$OMP DO
!do p = 1, mo_num
! hessian(p,p,p,p) = hessian(p,p,p,p) + tmp_h_pppp(p)
!enddo
!!$OMP END DO
!!$OMP DO
!do q = 1, mo_num
! do p = 1, mo_num
! hessian(p,q,p,q) = hessian(p,q,p,q) + tmp_h_pqpq(p,q)
! enddo
!enddo
!!$OMP END DO
!
!!$OMP DO
!do q = 1, mo_num
! do p = 1, mo_num
! hessian(p,q,q,p) = hessian(p,q,q,p) + tmp_h_pqqp(p,q)
! enddo
!enddo
!!$OMP END DO
!!$OMP DO
!do s = 1, mo_num
! do r = 1, mo_num
! do q = 1, mo_num
! do p = 1, mo_num
! h_tmpr(p,q,r,s) = (hessian(p,q,r,s) - hessian(q,p,r,s) - hessian(p,q,s,r) + hessian(q,p,s,r))
! enddo
! enddo
! enddo
!enddo
!!$OMP END DO4D -> 2D matrix
We need a 2D matrix for the Newton method's. Since the Hessian is "antisymmetric" : $$H_{pq,rs} = -H_{rs,pq}$$ We can write it as a 2D matrix, N by N, with N = mo_num(mo_num-1)/2 with p<q and r<s
!$OMP MASTER
CALL wall_TIME(t4)
!$OMP END MASTER
!$OMP DO
do pq = 1, n
call vec_to_mat_index(pq,p,q)
do rs = 1, n
call vec_to_mat_index(rs,r,s)
!H(pq,rs) = h_tmpr(p,q,r,s)
if (pq == rs) then
k = pq
if (p == r) then
H(k) = tmp_h_pqpq(p,q) + tmp_h_pqpq(q,p) - tmp_h_pqqp(p,q) - tmp_h_pqqp(q,p)
elseif (p == s) then
H(k) = - tmp_h_pqpq(p,q) - tmp_h_pqpq(q,p) + tmp_h_pqqp(p,q) + tmp_h_pqqp(q,p)
endif
endif
enddo
enddo
!$OMP END DO
!!$OMP MASTER
!call wall_TIME(t5)
!t6 = t5-t4
!print*,'4D -> 2D :',t6
!!$OMP END MASTER
!$OMP END PARALLEL
call omp_set_max_active_levels(4)
! Display
!if (debug) then
! print*,'2D diag Hessian matrix'
! do pq = 1, n
! write(*,'(100(F10.5))') H(pq,:)
! enddo
!endifDeallocation of shared arrays, end
!deallocate(hessian)!,h_tmpr)
deallocate(tmp_h_pppp,tmp_h_pqpq,tmp_h_pqqp)
deallocate(tmp_accu_1_shared, tmp_accu_shared)
print*,'---diagonal_hessian'
end subroutine