mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-01 21:53:40 +02:00
373 lines
10 KiB
Fortran
373 lines
10 KiB
Fortran
subroutine first_diag_hessian_list_opt(tmp_n,m,list,H)!, h_tmpr)
|
|
|
|
include 'constants.h'
|
|
|
|
implicit none
|
|
|
|
!===========================================================================
|
|
! Compute the diagonal hessian of energy with respects to orbital rotations
|
|
!===========================================================================
|
|
|
|
!===========
|
|
! Variables
|
|
!===========
|
|
|
|
! in
|
|
integer, intent(in) :: tmp_n, m, list(m)
|
|
! tmp_n : integer, tmp_n = m*(m-1)/2
|
|
|
|
! out
|
|
double precision, intent(out) :: H(tmp_n)!, h_tmpr(m,m,m,m)
|
|
! H : n by n double precision matrix containing the 2D hessian
|
|
|
|
! internal
|
|
double precision, allocatable :: hessian(:,:,:,:), tmp(:,:),h_tmpr(:,:,:,:)
|
|
integer :: p,q, tmp_p,tmp_q
|
|
integer :: r,s,t,u,v,tmp_r,tmp_s,tmp_t,tmp_u,tmp_v
|
|
integer :: pq,rs,tmp_pq,tmp_rs
|
|
double precision :: t1,t2,t3
|
|
! hessian : mo_num 4D double precision matrix containing the hessian before the permutations
|
|
! h_tmpr : mo_num 4D double precision matrix containing the hessian after the permutations
|
|
! p,q,r,s : integer, indexes of the 4D hessian matrix
|
|
! t,u,v : integer, indexes to compute hessian elements
|
|
! pq,rs : integer, indexes for the conversion from 4D to 2D hessian matrix
|
|
! t1,t2,t3 : double precision, t3 = t2 - t1, time to compute the hessian
|
|
|
|
! Function
|
|
double precision :: get_two_e_integral
|
|
! get_two_e_integral : double precision function, two e integrals
|
|
|
|
! Provided :
|
|
! mo_one_e_integrals : mono e- integrals
|
|
! get_two_e_integral : two e- integrals
|
|
! one_e_dm_mo_alpha, one_e_dm_mo_beta : one body density matrix
|
|
! two_e_dm_mo : two body density matrix
|
|
|
|
print*,'---first_diag_hess_list---'
|
|
|
|
!============
|
|
! Allocation
|
|
!============
|
|
|
|
allocate(hessian(m,m,m,m),tmp(tmp_n,tmp_n),h_tmpr(mo_num,mo_num,mo_num,mo_num))
|
|
|
|
!=============
|
|
! Calculation
|
|
!=============
|
|
|
|
! From Anderson et. al. (2014)
|
|
! The Journal of Chemical Physics 141, 244104 (2014); doi: 10.1063/1.4904384
|
|
|
|
! LaTeX formula :
|
|
|
|
!\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_u^r)]
|
|
!-(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_{ps}^{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*}
|
|
|
|
!================
|
|
! Initialization
|
|
!================
|
|
hessian = 0d0
|
|
|
|
CALL wall_time(t1)
|
|
|
|
!========================
|
|
! First line, first term
|
|
!========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! Permutations
|
|
if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
|
|
.or. ((p==s) .and. (q==r))) then
|
|
|
|
if (q==r) then
|
|
do u = 1, mo_num
|
|
|
|
hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_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
|
|
|
|
!=========================
|
|
! First line, second term
|
|
!=========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! Permutations
|
|
if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
|
|
.or. ((p==s) .and. (q==r))) then
|
|
|
|
if (p==s) then
|
|
do u = 1, mo_num
|
|
|
|
hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_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
|
|
|
|
!========================
|
|
! First line, third term
|
|
!========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! Permutations
|
|
if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
|
|
.or. ((p==s) .and. (q==r))) then
|
|
|
|
hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) &
|
|
- mo_one_e_integrals(s,p) * one_e_dm_mo(r,q) &
|
|
- mo_one_e_integrals(q,r) * one_e_dm_mo(p,s)
|
|
|
|
endif
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
!=========================
|
|
! Second line, first term
|
|
!=========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! Permutations
|
|
if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
|
|
.or. ((p==s) .and. (q==r))) then
|
|
|
|
if (q==r) then
|
|
do t = 1, mo_num
|
|
do u = 1, mo_num
|
|
do v = 1, mo_num
|
|
|
|
hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_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
|
|
|
|
!==========================
|
|
! Second line, second term
|
|
!==========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! Permutations
|
|
if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
|
|
.or. ((p==s) .and. (q==r))) then
|
|
|
|
if (p==s) then
|
|
do t = 1, mo_num
|
|
do u = 1, mo_num
|
|
do v = 1, mo_num
|
|
|
|
hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_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
|
|
|
|
!========================
|
|
! Third line, first term
|
|
!========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! Permutations
|
|
if (((p==r) .and. (q==s)) .or. ((q==r) .and. (p==s)) &
|
|
.or. ((p==s) .and. (q==r))) then
|
|
|
|
do u = 1, mo_num
|
|
do v = 1, mo_num
|
|
|
|
hessian(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_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
|
|
|
|
!=========================
|
|
! Third line, second term
|
|
!=========================
|
|
do tmp_p = 1, m
|
|
p = list(tmp_p)
|
|
do tmp_q = 1, m
|
|
q = list(tmp_q)
|
|
do tmp_r = 1, m
|
|
r = list(tmp_r)
|
|
do tmp_s = 1, m
|
|
s = list(tmp_s)
|
|
|
|
! 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(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_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
|
|
|
|
CALL wall_time(t2)
|
|
t2 = t2 - t1
|
|
print*, 'Time to compute the hessian :', t2
|
|
|
|
!==============
|
|
! Permutations
|
|
!==============
|
|
|
|
! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
|
|
! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
|
|
! H(pq,rs) : p<q and r<s
|
|
|
|
do tmp_r = 1, m
|
|
do tmp_s = 1, m
|
|
do tmp_q = 1, m
|
|
do tmp_p = 1, m
|
|
|
|
h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s) = hessian(tmp_p,tmp_q,tmp_r,tmp_s) - hessian(tmp_q,tmp_p,tmp_r,tmp_s) &
|
|
- hessian(tmp_p,tmp_q,tmp_s,tmp_r) + hessian(tmp_q,tmp_p,tmp_s,tmp_r)
|
|
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
!========================
|
|
! 4D matrix -> 2D matrix
|
|
!========================
|
|
|
|
! Convert the hessian mo_num * mo_num * mo_num * mo_num matrix in a
|
|
! 2D n * n matrix (n = mo_num*(mo_num-1)/2)
|
|
! H(pq,rs) : p<q and r<s
|
|
|
|
! 4D mo_num matrix to 2D n matrix
|
|
do tmp_rs = 1, tmp_n
|
|
call vec_to_mat_index(tmp_rs,tmp_r,tmp_s)
|
|
do tmp_pq = 1, tmp_n
|
|
call vec_to_mat_index(tmp_pq,tmp_p,tmp_q)
|
|
tmp(tmp_pq,tmp_rs) = h_tmpr(tmp_p,tmp_q,tmp_r,tmp_s)
|
|
enddo
|
|
enddo
|
|
|
|
do p = 1, tmp_n
|
|
H(p) = tmp(p,p)
|
|
enddo
|
|
|
|
! Display
|
|
if (debug) then
|
|
print*,'2D diag Hessian matrix'
|
|
do tmp_pq = 1, tmp_n
|
|
write(*,'(100(F10.5))') tmp(tmp_pq,:)
|
|
enddo
|
|
endif
|
|
|
|
!==============
|
|
! Deallocation
|
|
!==============
|
|
|
|
deallocate(hessian,h_tmpr,tmp)
|
|
|
|
print*,'---End first_diag_hess_list---'
|
|
|
|
end subroutine
|