PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-09 21:03:58 +02:00
Code Releases Activity
94662d3da0
QuantumPackage/src/ccsd/org/ccsd_space_orb.org
ydamour fadbddc869 add ccsd
2023-03-13 14:08:32 +01:00

2122 lines
53 KiB
Org Mode
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

* ccsd with spatial orbitals
Scuseria, Gustavo E.; Janssen, Curtis L.; Schaefer, Henry
F. (1988). An efficient reformulation of the closed-shell coupled
cluster single and double excitation (CCSD) equations. The Journal of
Chemical Physics, 89(12), 7382–. doi:10.1063/1.455269
* Code
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb.irp.f
program ccsd
implicit none
read_wf = .True.
touch read_wf
call run_ccsd_space_orb
end
#+end_src
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine run_ccsd_space_orb
implicit none
integer :: i,j,k,l,a,b,c,d,tmp_a,tmp_b,tmp_c,tmp_d
integer :: u,v,gam,beta,tmp_gam,tmp_beta
integer :: nb_iter
double precision :: get_two_e_integral
double precision :: uncorr_energy,energy, max_elem, max_r, max_r1, max_r2,ta,tb
logical :: not_converged
double precision, allocatable :: t2(:,:,:,:), r2(:,:,:,:), tau(:,:,:,:)
double precision, allocatable :: t1(:,:), r1(:,:)
double precision, allocatable :: H_oo(:,:), H_vv(:,:), H_vo(:,:)
double precision, allocatable :: all_err(:,:), all_t(:,:)
integer, allocatable :: list_occ(:), list_vir(:)
integer(bit_kind) :: det(N_int,2)
integer :: nO, nV, nOa, nOb, nVa, nVb, n_spin(4)
PROVIDE mo_two_e_integrals_in_map
det = psi_det(:,:,cc_ref)
print*,'Reference determinant:'
call print_det(det,N_int)
! Extract number of occ/vir alpha/beta spin orbitals
!call extract_n_spin(det,n_spin)
nOa = cc_nOa !n_spin(1)
nOb = cc_nOb !n_spin(2)
nVa = cc_nVa !n_spin(3)
nVb = cc_nVb !n_spin(4)
! Check that the reference is a closed shell determinant
if (cc_ref_is_open_shell) then
call abort
endif
! Number of occ/vir spatial orb
nO = nOa
nV = nVa
allocate(list_occ(nO),list_vir(nV))
list_occ = cc_list_occ
list_vir = cc_list_vir
! Debug
!call extract_list_orb_space(det,nO,nV,list_occ,list_vir)
!print*,'occ',list_occ
!print*,'vir',list_vir
allocate(t2(nO,nO,nV,nV), r2(nO,nO,nV,nV))
allocate(tau(nO,nO,nV,nV))
allocate(t1(nO,nV), r1(nO,nV))
allocate(H_oo(nO,nO), H_vv(nV,nV), H_vo(nV,nO))
if (cc_update_method == 'diis') then
allocate(all_err(nO*nV+nO*nO*nV*nV,cc_diis_depth), all_t(nO*nV+nO*nO*nV*nV,cc_diis_depth))
all_err = 0d0
all_t = 0d0
endif
if (elec_alpha_num /= elec_beta_num) then
print*, 'Only for closed shell systems'
print*, 'elec_alpha_num=',elec_alpha_num
print*, 'elec_beta_num =',elec_beta_num
print*, 'abort'
call abort
endif
! Init
call guess_t1(nO,nV,cc_space_f_o,cc_space_f_v,cc_space_f_ov,t1)
call guess_t2(nO,nV,cc_space_f_o,cc_space_f_v,cc_space_v_oovv,t2)
call update_tau_space(nO,nV,t1,t2,tau)
!print*,'hf_energy', hf_energy
call det_energy(det,uncorr_energy)
print*,'Det energy', uncorr_energy
call ccsd_energy_space(nO,nV,tau,t1,energy)
print*,'Guess energy', uncorr_energy+energy, energy
nb_iter = 0
not_converged = .True.
max_r1 = 0d0
max_r2 = 0d0
write(*,'(A77)') ' -----------------------------------------------------------------------------'
write(*,'(A77)') ' | It. | E(CCSD) (Ha) | Correlation (Ha) | Conv. T1 | Conv. T2 |'
write(*,'(A77)') ' -----------------------------------------------------------------------------'
call wall_time(ta)
do while (not_converged)
call compute_H_oo(nO,nV,t1,t2,tau,H_oo)
call compute_H_vv(nO,nV,t1,t2,tau,H_vv)
call compute_H_vo(nO,nV,t1,t2,H_vo)
! Residue
call compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
call compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
max_r = max(max_r1,max_r2)
! Update
if (cc_update_method == 'diis') then
!call update_t_ccsd(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
!call update_t_ccsd_diis(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2)
call update_t_ccsd_diis_v3(nO,nV,nb_iter,cc_space_f_o,cc_space_f_v,r1,r2,t1,t2,all_err,all_t)
! Standard update as T = T - Delta
elseif (cc_update_method == 'none') then
call update_t1(nO,nV,cc_space_f_o,cc_space_f_v,r1,t1)
call update_t2(nO,nV,cc_space_f_o,cc_space_f_v,r2,t2)
else
print*,'Unkonw cc_method_method: '//cc_update_method
endif
call update_tau_space(nO,nV,t1,t2,tau)
! Energy
call ccsd_energy_space(nO,nV,tau,t1,energy)
write(*,'(A3,I6,A3,F18.12,A3,F16.12,A3,1pE10.2,A3,1pE10.2,A2)') ' | ',nb_iter,' | ', uncorr_energy+energy,' | ', energy,' | ', max_r1,' | ', max_r2,' |'
nb_iter = nb_iter + 1
if (max_r < cc_thresh_conv .or. nb_iter > cc_max_iter) then
not_converged = .False.
endif
enddo
write(*,'(A77)') ' -----------------------------------------------------------------------------'
call wall_time(tb)
print*,'Time: ',tb-ta, ' s'
print*,''
if (max_r < cc_thresh_conv) then
write(*,'(A30,I6,A11)') ' Successful convergence after ', nb_iter, ' iterations'
else
write(*,'(A26,I6,A11)') ' Failed convergence after ', nb_iter, ' iterations'
endif
print*,''
write(*,'(A15,F18.12,A3)') ' E(CCSD) = ', uncorr_energy+energy, ' Ha'
write(*,'(A15,F18.12,A3)') ' Correlation = ', energy, ' Ha'
write(*,'(A15,1pE10.2,A3)')' Conv = ', max_r
print*,''
call write_t1(nO,nV,t1)
call write_t2(nO,nV,t2)
! Deallocation
if (cc_update_method == 'diis') then
deallocate(all_err,all_t)
endif
deallocate(H_vv,H_oo,H_vo,r1,r2,tau)
! CCSD(T)
double precision :: e_t
if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then
! Dumb way
!call wall_time(ta)
!call ccsd_par_t_space(nO,nV,t1,t2,e_t)
!call wall_time(tb)
!print*,'Time: ',tb-ta, ' s'
!print*,''
!write(*,'(A15,F18.12,A3)') ' E(CCSD(T)) = ', uncorr_energy + energy + e_t, ' Ha'
!write(*,'(A15,F18.12,A3)') ' E(T) = ', e_t, ' Ha'
!write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + e_t, ' Ha'
!print*,''
! New
print*,'Computing (T) correction...'
call wall_time(ta)
call ccsd_par_t_space_v2(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v &
,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t)
call wall_time(tb)
print*,'Time: ',tb-ta, ' s'
print*,''
write(*,'(A15,F18.12,A3)') ' E(CCSD(T)) = ', uncorr_energy + energy + e_t, ' Ha'
write(*,'(A15,F18.12,A3)') ' E(T) = ', e_t, ' Ha'
write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + e_t, ' Ha'
print*,''
endif
print*,'Reference determinant:'
call print_det(det,N_int)
deallocate(t1,t2)
end
#+END_SRC
* Energy
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine ccsd_energy_space(nO,nV,tau,t1,energy)
implicit none
integer, intent(in) :: nO, nV
double precision, intent(in) :: tau(nO,nO,nV,nV)
double precision, intent(in) :: t1(nO,nV)
double precision, intent(out) :: energy
! internal
integer :: i,j,a,b
double precision :: e
energy = 0d0
!$omp parallel &
!$omp shared(nO,nV,energy,tau,t1,&
!$omp cc_space_f_vo,cc_space_w_oovv) &
!$omp private(i,j,a,b,e) &
!$omp default(none)
e = 0d0
!$omp do
do i = 1, nO
do a = 1, nV
e = e + 2d0 * cc_space_f_vo(a,i) * t1(i,a)
enddo
enddo
!$omp end do nowait
!$omp do
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
e = e + tau(i,j,a,b) * cc_space_w_oovv(i,j,a,b)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp critical
energy = energy + e
!$omp end critical
!$omp end parallel
end
#+END_SRC
* T
** Tau
#+begin_src f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine update_tau_space(nO,nV,t1,t2,tau)
implicit none
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV)
! out
double precision, intent(out) :: tau(nO,nO,nV,nV)
! internal
integer :: i,j,a,b
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,tau,t2,t1) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
tau(i,j,a,b) = t2(i,j,a,b) + t1(i,a) * t1(j,b)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end
#+end_src
* Residual equations
** R1
*** R1
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
implicit none
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV), tau(nO,nO,nV,nV)
double precision, intent(in) :: H_oo(nO,nO), H_vv(nV,nV), H_vo(nV,nO)
! out
double precision, intent(out) :: r1(nO,nV), max_r1
! internal
integer :: u,i,j,beta,a,b
!$omp parallel &
!$omp shared(nO,nV,r1,cc_space_f_ov) &
!$omp private(u,beta) &
!$omp default(none)
!$omp do
do beta = 1, nV
do u = 1, nO
r1(u,beta) = cc_space_f_ov(u,beta)
enddo
enddo
!$omp end do
!$omp end parallel
! r1(u,beta) = r1(u,beta) - 2d0 * cc_space_f_vo(a,i) * t1(i,beta) * t1(u,a)
! cc_space_f_vo(a,i) * t1(i,beta) -> X1(nV,nV), O(nV*nV*nO)
! X1(a,beta) * t1(u,a) -> O(nO*nV*nV)
! cc_space_f_vo(a,i) * t1(u,a) -> X1(nO,nO), O(nO*nO*nV)
! X1(i,u) * t1(i,beta) -> O(nO*nO*nV)
!do beta = 1, nV
! do u = 1, nO
! do i = 1, nO
! do a = 1, nV
! r1(u,beta) = r1(u,beta) - 2d0 * cc_space_f_vo(a,i) * t1(i,beta) * t1(u,a)
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_oo(:,:)
allocate(X_oo(nO,nO))
call dgemm('N','N', nO, nO, nV, &
-2d0, t1 , size(t1,1), &
cc_space_f_vo, size(cc_space_f_vo,1), &
0d0, X_oo , size(X_oo,1))
call dgemm('T','N', nO, nV, nO, &
1d0, X_oo, size(X_oo,2), &
t1 , size(t1,1), &
1d0, r1 , size(r1,1))
deallocate(X_oo)
! r1(u,beta) = r1(u,beta) + H_vv(a,beta) * t1(u,a)
!do beta = 1, nV
! do u = 1, nO
! do a = 1, nV
! r1(u,beta) = r1(u,beta) + H_vv(a,beta) * t1(u,a)
! enddo
! enddo
!enddo
call dgemm('N','N', nO, nV, nV, &
1d0, t1 , size(t1,1), &
H_vv, size(H_vv,1), &
1d0, r1 , size(r1,1))
! r1(u,beta) = r1(u,beta) - H_oo(u,i) * t1(i,beta)
!do beta = 1, nV
! do u = 1, nO
! do i = 1, nO
! r1(u,beta) = r1(u,beta) - H_oo(u,i) * t1(i,beta)
! enddo
! enddo
!enddo
call dgemm('N','N', nO, nV, nO, &
-1d0, H_oo, size(H_oo,1), &
t1 , size(t1,1), &
1d0, r1, size(r1,1))
!r1(u,beta) = r1(u,beta) + H_vo(a,i) * (2d0 * t2(i,u,a,beta) - t2(u,i,a,beta) + t1(u,a) * t1(i,beta))
! <=>
! r1(u,beta) = r1(u,beta) + H_vo(a,i) * X(a,i,u,beta)
!do beta = 1, nV
! do u = 1, nO
! do i = 1, nO
! do a = 1, nV
! r1(u,beta) = r1(u,beta) + H_vo(a,i) * &
! (2d0 * t2(i,u,a,beta) - t2(u,i,a,beta) + t1(u,a) * t1(i,beta))
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_voov(:,:,:,:)
allocate(X_voov(nV, nO, nO, nV))
!$omp parallel &
!$omp shared(nO,nV,X_voov,t2,t1) &
!$omp private(u,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
do beta = 1, nV
do u = 1, nO
do i = 1, nO
do a = 1, nV
X_voov(a,i,u,beta) = 2d0 * t2(i,u,a,beta) - t2(u,i,a,beta) + t1(u,a) * t1(i,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemv('T', nV*nO, nO*nV, &
1d0, X_voov, size(X_voov,1) * size(X_voov,2), &
H_vo , 1, &
1d0, r1 , 1)
deallocate(X_voov)
! r1(u,beta) = r1(u,beta) + (2d0 * cc_space_v_voov(a,u,i,beta) - cc_space_v_ovov(u,a,i,beta)) * t1(i,a)
! <=>
! r1(u,beta) = r1(u,beta) + X(i,a,u,beta)
!do beta = 1, nV
! do u = 1, nO
! do i = 1, nO
! do a = 1, nV
! r1(u,beta) = r1(u,beta) + (2d0 * cc_space_v_voov(a,u,i,beta) - cc_space_v_ovov(u,a,i,beta)) * t1(i,a)
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_ovov(:,:,:,:)
allocate(X_ovov(nO, nV, nO, nV))
!$omp parallel &
!$omp shared(nO,nV,cc_space_v_ovov,cc_space_v_voov,X_ovov) &
!$omp private(u,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
do beta = 1, nV
do u = 1, nO
do a = 1, nv
do i = 1, nO
X_ovov(i,a,u,beta) = 2d0 * cc_space_v_voov(a,u,i,beta) - cc_space_v_ovov(u,a,i,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemv('T', nO*nV, nO*nV, &
1d0, X_ovov, size(X_ovov,1) * size(X_ovov,2), &
t1 , 1, &
1d0, r1 , 1)
deallocate(X_ovov)
! r1(u,beta) = r1(u,beta) + (2d0 * cc_space_v_vvov(a,b,i,beta) - cc_space_v_vvov(b,a,i,beta)) * tau(i,u,a,b)
! r1(u,beta) = r1(u,beta) + W(a,b,i,beta) * T(u,a,b,i)
!do beta = 1, nV
! do u = 1, nO
! do i = 1, nO
! do a = 1, nV
! do b = 1, nV
! r1(u,beta) = r1(u,beta) + (2d0 * cc_space_v_vvov(a,b,i,beta) - cc_space_v_vvov(b,a,i,beta)) * tau(i,u,a,b)
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: W_vvov(:,:,:,:), T_vvoo(:,:,:,:)
allocate(W_vvov(nV,nV,nO,nV), T_vvoo(nV,nV,nO,nO))
!$omp parallel &
!$omp shared(nO,nV,cc_space_v_vvov,W_vvov,T_vvoo,tau) &
!$omp private(b,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
do beta = 1, nV
do i = 1, nO
do b = 1, nV
do a = 1, nV
W_vvov(a,b,i,beta) = 2d0 * cc_space_v_vvov(a,b,i,beta) - cc_space_v_vvov(b,a,i,beta)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do i = 1, nO
do b = 1, nV
do a = 1, nV
do u = 1, nO
T_vvoo(a,b,i,u) = tau(i,u,a,b)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp end parallel
call dgemm('T','N',nO,nV,nO*nV*nV, &
1d0, T_vvoo, size(T_vvoo,1) * size(T_vvoo,2) * size(T_vvoo,3), &
W_vvov, size(W_vvov,1) * size(W_vvov,2) * size(W_vvov,3), &
1d0, r1 , size(r1,1))
deallocate(W_vvov,T_vvoo)
! r1(u,beta) = r1(u,beta) - (2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)) * tau(i,j,a,beta)
! r1(u,beta) = r1(u,beta) - W(i,j,a,u) * tau(i,j,a,beta)
!do beta = 1, nV
! do u = 1, nO
! do i = 1, nO
! do j = 1, nO
! do a = 1, nV
! r1(u,beta) = r1(u,beta) - (2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)) * tau(i,j,a,beta)
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: W_oovo(:,:,:,:)
allocate(W_oovo(nO,nO,nV,nO))
!$omp parallel &
!$omp shared(nO,nV,cc_space_v_vooo,W_oovo) &
!$omp private(u,a,i,j) &
!$omp default(none)
!$omp do collapse(3)
do u = 1, nO
do a = 1, nV
do j = 1, nO
do i = 1, nO
W_oovo(i,j,a,u) = 2d0 * cc_space_v_vooo(a,u,i,j) - cc_space_v_vooo(a,u,j,i)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('T','N', nO, nV, nO*nO*nV, &
-1d0, W_oovo, size(W_oovo,1) * size(W_oovo,2) * size(W_oovo,3), &
tau , size(tau,1) * size(tau,2) * size(tau,3), &
1d0, r1 , size(r1,1))
deallocate(W_oovo)
max_r1 = 0d0
do a = 1, nV
do i = 1, nO
if (dabs(r1(i,a)) > max_r1) then
max_r1 = dabs(r1(i,a))
endif
enddo
enddo
! Change the sign for consistency with the code in spin orbitals
!$omp parallel &
!$omp shared(nO,nV,r1) &
!$omp private(a,i) &
!$omp default(none)
!$omp do
do a = 1, nV
do i = 1, nO
r1(i,a) = -r1(i,a)
enddo
enddo
!$omp end do
!$omp end parallel
end
#+end_src
*** Intermediates
**** H_oo
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_H_oo(nO,nV,t1,t2,tau,H_oo)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(in) :: tau(nO, nO, nV, nV)
double precision, intent(out) :: H_oo(nO, nO)
integer :: a,tmp_a,k,b,l,c,d,tmp_c,tmp_d,i,j,u
!H_oo = 0d0
!do i = 1, nO
! do u = 1, nO
! H_oo(u,i) = cc_space_f_oo(u,i)
! do j = 1, nO
! do a = 1, nV
! do b = 1, nV
! !H_oo(u,i) = H_oo(u,i) + (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(b,a,i,j)) * tau(u,j,a,b)
! !H_oo(u,i) = H_oo(u,i) + cc_space_w_vvoo(a,b,i,j) * tau(u,j,a,b)
! H_oo(u,i) = H_oo(u,i) + cc_space_w_oovv(i,j,a,b) * tau(u,j,a,b)
! enddo
! enddo
! enddo
!
! enddo
!enddo
! H_oo(u,i) = cc_space_f_oo(u,i)
!$omp parallel &
!$omp shared(nO,H_oo,cc_space_f_oo) &
!$omp private(i,u) &
!$omp default(none)
!$omp do
do i = 1, nO
do u = 1, nO
H_oo(u,i) = cc_space_f_oo(u,i)
enddo
enddo
!$omp end do
!$omp end parallel
! H_oo(u,i) += cc_space_w_oovv(i,j,a,b) * tau(u,j,a,b)
! H_oo(u,i) += tau(u,j,a,b) * cc_space_w_oovv(i,j,a,b)
call dgemm('N','T', nO, nO, nO*nV*nV, &
1d0, tau , size(tau,1), &
cc_space_w_oovv, size(cc_space_w_oovv,1), &
1d0, H_oo , size(H_oo,1))
end
#+END_SRC
**** H_vv
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_H_vv(nO,nV,t1,t2,tau,H_vv)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(in) :: tau(nO, nO, nV, nV)
double precision, intent(out) :: H_vv(nV, nV)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u, beta
!H_vv = 0d0
!do beta = 1, nV
! do a = 1, nV
! H_vv(a,beta) = cc_space_f_vv(a,beta)
! do j = 1, nO
! do i = 1, nO
! do b = 1, nV
! !H_vv(a,beta) = H_vv(a,beta) - (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(a,b,j,i)) * tau(i,j,beta,b)
! H_vv(a,beta) = H_vv(a,beta) - cc_space_w_vvoo(a,b,i,j) * tau(i,j,beta,b)
! enddo
! enddo
! enddo
!
! enddo
!enddo
double precision, allocatable :: tmp_tau(:,:,:,:)
allocate(tmp_tau(nV,nO,nO,nV))
! H_vv(a,beta) = cc_space_f_vv(a,beta)
!$omp parallel &
!$omp shared(nV,nO,H_vv,cc_space_f_vv,tmp_tau,tau) &
!$omp private(a,beta,i,j,b) &
!$omp default(none)
!$omp do
do beta = 1, nV
do a = 1, nV
H_vv(a,beta) = cc_space_f_vv(a,beta)
enddo
enddo
!$omp end do nowait
! H_vv(a,beta) = H_vv(a,beta) - cc_space_w_vvoo(a,b,i,j) * tau(i,j,beta,b)
! H_vv(a,beta) = H_vv(a,beta) - cc_space_w_vvoo(a,b,i,j) * tmp_tau(b,i,j,beta)
!$omp do collapse(3)
do beta = 1, nV
do j = 1, nO
do i = 1, nO
do b = 1, nV
tmp_tau(b,i,j,beta) = tau(i,j,beta,b)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nV,nV,nO*nO*nV, &
-1d0, cc_space_w_vvoo, size(cc_space_w_vvoo,1), &
tmp_tau , size(tmp_tau,1) * size(tmp_tau,2) * size(tmp_tau,3), &
1d0, H_vv , size(H_vv,1))
deallocate(tmp_tau)
end
#+END_SRC
**** H_vo
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_H_vo(nO,nV,t1,t2,H_vo)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: H_vo(nV, nO)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u, beta
!H_vo = 0d0
!do i = 1, nO
! do a = 1, nV
! H_vo(a,i) = cc_space_f_vo(a,i)
! do j = 1, nO
! do b = 1, nV
! !H_vo(a,i) = H_vo(a,i) + (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(b,a,i,j)) * t1(j,b)
! H_vo(a,i) = H_vo(a,i) + cc_space_w_vvoo(a,b,i,j) * t1(j,b)
! enddo
! enddo
!
! enddo
!enddo
double precision, allocatable :: w(:,:,:,:)
allocate(w(nV,nO,nO,nV))
!$omp parallel &
!$omp shared(nV,nO,H_vo,cc_space_f_vo,w,cc_space_w_vvoo,t1) &
!$omp private(a,beta,i,j,b) &
!$omp default(none)
!$omp do
do i = 1, nO
do a = 1, nV
H_vo(a,i) = cc_space_f_vo(a,i)
enddo
enddo
!$omp end do nowait
! H_vo(a,i) = H_vo(a,i) + cc_space_w_vvoo(a,b,i,j) * t1(j,b)
! H_vo(a,i) = H_vo(a,i) + w(a,i,j,b) * t1(j,b)
!$omp do collapse(3)
do b = 1, nV
do j = 1, nO
do i = 1, nO
do a = 1, nV
w(a,i,j,b) = cc_space_w_vvoo(a,b,i,j)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemv('N',nV*nO, nO*nV, &
1d0, w , size(w,1) * size(w,2), &
t1 , 1, &
1d0, H_vo, 1)
deallocate(w)
end
#+END_SRC
** R2
*** R2
#+begin_src f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
implicit none
! in
integer, intent(in) :: nO, nV
double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV), tau(nO,nO,nV,nV)
double precision, intent(in) :: H_oo(nO,nO), H_vv(nV,nV), H_vo(nV,nO)
! out
double precision, intent(out) :: r2(nO,nO,nV,nV), max_r2
! internal
double precision, allocatable :: g_occ(:,:), g_vir(:,:), J1(:,:,:,:), K1(:,:,:,:)
double precision, allocatable :: A1(:,:,:,:), B1(:,:,:,:)
integer :: u,v,i,j,beta,gam,a,b
allocate(g_occ(nO,nO), g_vir(nV,nV))
allocate(J1(nO,nV,nV,nO), K1(nO,nV,nO,nV))
allocate(A1(nO,nO,nO,nO))
call compute_g_occ(nO,nV,t1,t2,H_oo,g_occ)
call compute_g_vir(nO,nV,t1,t2,H_vv,g_vir)
call compute_A1(nO,nV,t1,t2,tau,A1)
call compute_J1(nO,nV,t1,t2,cc_space_v_ovvo,cc_space_v_ovoo, &
cc_space_v_vvvo,cc_space_v_vvoo,J1)
call compute_K1(nO,nV,t1,t2,cc_space_v_ovoo,cc_space_v_vvoo, &
cc_space_v_ovov,cc_space_v_vvov,K1)
! Residual
!r2 = 0d0
!$omp parallel &
!$omp shared(nO,nV,r2,cc_space_v_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = cc_space_v_oovv(u,v,beta,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do j = 1, nO
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! + A1(u,v,i,j) * tau(i,j,beta,gam)
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
call dgemm('N','N',nO*nO,nV*nV,nO*nO, &
1d0, A1, size(A1,1) * size(A1,2), &
tau, size(tau,1) * size(tau,2), &
1d0, r2, size(r2,1) * size(r2,2))
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! do b = 1, nv
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! + B1(a,b,beta,gam) * tau(u,v,a,b)
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
allocate(B1(nV,nV,nV,nV))
call compute_B1(nO,nV,t1,t2,B1)
call dgemm('N','N',nO*nO,nV*nV,nV*nV, &
1d0, tau, size(tau,1) * size(tau,2), &
B1 , size(B1,1) * size(B1,2), &
1d0, r2, size(r2,1) * size(r2,2))
deallocate(B1)
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! + g_vir(a,beta) * t2(u,v,a,gam) &
! + g_vir(a,gam) * t2(v,u,a,beta) ! P
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_oovv(:,:,:,:),Y_oovv(:,:,:,:)
allocate(X_oovv(nO,nO,nV,nV),Y_oovv(nO,nO,nV,nV))
!$omp parallel &
!$omp shared(nO,nV,t2,X_oovv) &
!$omp private(u,v,gam,a) &
!$omp default(none)
!$omp do collapse(3)
do a = 1, nV
do gam = 1, nV
do v = 1, nO
do u = 1, nO
X_oovv(u,v,gam,a) = t2(u,v,gam,a)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO*nO*nV,nV,nV, &
1d0, X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3), &
g_vir, size(g_vir,1), &
0d0, Y_oovv, size(Y_oovv,1) * size(Y_oovv,2) * size(Y_oovv,3))
!$omp parallel &
!$omp shared(nO,nV,r2,Y_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) + Y_oovv(u,v,beta,gam) + Y_oovv(v,u,gam,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! - g_occ(u,i) * t2(i,v,beta,gam) &
! - g_occ(v,i) * t2(i,u,gam,beta) ! P
! enddo
! enddo
! enddo
! enddo
!enddo
call dgemm('N','N',nO,nO*nV*nV,nO, &
1d0, g_occ , size(g_occ,1), &
t2 , size(t2,1), &
0d0, X_oovv, size(X_oovv,1))
!$omp parallel &
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) - X_oovv(u,v,beta,gam) - X_oovv(v,u,gam,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_oovv)
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! + cc_space_v_ovvv(u,a,beta,gam) * t1(v,a) &
! + cc_space_v_ovvv(v,a,gam,beta) * t1(u,a) ! P
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_vovv(:,:,:,:)
allocate(X_vovv(nV,nO,nV,nV))
!$omp parallel &
!$omp shared(nO,nV,X_vovv,cc_space_v_ovvv) &
!$omp private(u,a,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do u = 1, nO
do a = 1, nV
X_vovv(a,u,beta,gam) = cc_space_v_ovvv(u,a,beta,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO,nO*nV*nV,nV, &
1d0, t1 , size(t1,1), &
X_vovv, size(X_vovv,1), &
0d0, Y_oovv, size(Y_oovv,1))
!$omp parallel &
!$omp shared(nO,nV,r2,Y_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) + Y_oovv(v,u,beta,gam) + Y_oovv(u,v,gam,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! - cc_space_v_ovov(u,a,i,gam) * t1(i,beta) * t1(v,a) &
! - cc_space_v_ovov(v,a,i,beta) * t1(i,gam) * t1(u,a) ! P
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_vovo(:,:,:,:), Y_vovv(:,:,:,:)
allocate(X_vovo(nV,nO,nV,nO), Y_vovv(nV,nO,nV,nV),X_oovv(nO,nO,nV,nV))
!$omp parallel &
!$omp shared(nO,nV,X_vovo,cc_space_v_ovov) &
!$omp private(u,v,gam,i) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do gam = 1, nV
do u = 1, nO
do a = 1, nV
X_vovo(a,u,gam,i) = cc_space_v_ovov(u,a,i,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nV*nO*nV,nV,nO, &
1d0, X_vovo, size(X_vovo,1) * size(X_vovo,2) * size(X_vovo,3), &
t1 , size(t1,1), &
0d0, Y_vovv, size(Y_vovv,1) * size(Y_vovv,2) * size(Y_vovv,3))
call dgemm('N','N',nO,nO*nV*nV,nV, &
1d0, t1, size(t1,1), &
Y_vovv, size(Y_vovv,1), &
0d0, X_oovv, size(X_oovv,1))
!$omp parallel &
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) - X_oovv(v,u,gam,beta) - X_oovv(u,v,beta,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_vovo,Y_vovv)
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! - cc_space_v_oovo(u,v,beta,i) * t1(i,gam) &
! - cc_space_v_oovo(v,u,gam,i) * t1(i,beta) ! P
! enddo
! enddo
! enddo
! enddo
!enddo
call dgemm('N','N',nO*nO*nV,nV,nO, &
1d0, cc_space_v_oovo, size(cc_space_v_oovo,1) * size(cc_space_v_oovo,2) * size(cc_space_v_oovo,3), &
t1 , size(t1,1), &
0d0, X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3))
!$omp parallel &
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) - X_oovv(u,v,beta,gam) - X_oovv(v,u,gam,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do i = 1, nO
! do a = 1, nV
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! - cc_space_v_ovvo(u,a,beta,i) * t1(v,a) * t1(i,gam) &
! - cc_space_v_ovvo(v,a,gam,i) * t1(u,a) * t1(i,beta) ! P
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: Y_oovo(:,:,:,:)
allocate(X_vovo(nV,nO,nV,nO), Y_oovo(nO,nO,nV,nO))
!$omp parallel &
!$omp shared(nO,nV,X_vovo,cc_space_v_ovvo) &
!$omp private(a,v,gam,i) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do gam = 1, nV
do v = 1, nO
do a = 1, nV
X_vovo(a,v,gam,i) = cc_space_v_ovvo(v,a,gam,i)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO,nO*nV*nO,nV, &
1d0, t1, size(t1,1), &
X_vovo, size(X_vovo,1), &
0d0, Y_oovo, size(Y_oovo,1))
call dgemm('N','N',nO*nO*nV, nV, nO, &
1d0, Y_oovo, size(Y_oovo,1) * size(Y_oovo,2) * size(Y_oovo,3), &
t1 , size(t1,1), &
0d0, X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3))
!$omp parallel &
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) - X_oovv(u,v,gam,beta) - X_oovv(v,u,beta,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_vovo,Y_oovo)
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! + 0.5d0 * (2d0 * J1(u,a,beta,i) - K1(u,a,i,beta)) * &
! (2d0 * t2(i,v,a,gam) - t2(i,v,gam,a)) &
! + 0.5d0 * (2d0 * J1(v,a,gam,i) - K1(v,a,i,gam)) * &
! (2d0 * t2(i,u,a,beta) - t2(i,u,beta,a)) ! P
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_ovvo(:,:,:,:), Y_voov(:,:,:,:), Z_ovov(:,:,:,:)
allocate(X_ovvo(nO,nV,nV,nO), Y_voov(nV,nO,nO,nV),Z_ovov(nO,nV,nO,nV))
!$omp parallel &
!$omp shared(nO,nV,X_ovvo,Y_voov,K1,J1,t2) &
!$omp private(u,v,gam,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do a = 1, nV
do beta = 1, nV
do u = 1, nO
X_ovvo(u,beta,a,i) = 0.5d0 * (2d0 * J1(u,a,beta,i) - K1(u,a,i,beta))
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do gam = 1, nV
do v = 1, nO
do i = 1, nO
do a = 1, nV
Y_voov(a,i,v,gam) = 2d0 * t2(i,v,a,gam) - t2(i,v,gam,a)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N', nO*nV,nO*nV,nV*nO, &
1d0, X_ovvo, size(X_ovvo,1) * size(X_ovvo,2), &
Y_voov, size(Y_voov,1) * size(Y_voov,2), &
0d0, Z_ovov, size(Z_ovov,1) * size(Z_ovov,2))
!$omp parallel &
!$omp shared(nO,nV,r2,Z_ovov) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) + Z_ovov(u,beta,v,gam) + Z_ovov(v,gam,u,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_ovvo,Y_voov)
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! - 0.5d0 * K1(u,a,i,beta) * t2(i,v,gam,a) &
! - 0.5d0 * K1(v,a,i,gam) * t2(i,u,beta,a) !P
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_ovov(:,:,:,:),Y_ovov(:,:,:,:)
allocate(X_ovov(nO,nV,nO,nV),Y_ovov(nO,nV,nO,nV))
!$omp parallel &
!$omp shared(nO,nV,r2,K1,X_ovov,Y_ovov,t2) &
!$omp private(u,a,i,beta,gam) &
!$omp default(none)
!$omp do collapse(3)
do beta = 1, nV
do u = 1, nO
do a = 1, nV
do i = 1, nO
X_ovov(i,a,u,beta) = 0.5d0 * K1(u,a,i,beta)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do gam = 1, nV
do v = 1, nO
do a = 1, nV
do i = 1, nO
Y_ovov(i,a,v,gam) = t2(i,v,gam,a)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('T','N',nO*nV,nO*nV,nO*nV, &
1d0, X_ovov, size(X_ovov,1) * size(X_ovov,2), &
Y_ovov, size(Y_ovov,1) * size(Y_ovov,2), &
0d0, Z_ovov, size(Y_ovov,1) * size(Y_ovov,2))
!$omp parallel &
!$omp shared(nO,nV,r2,Z_ovov) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) - Z_ovov(u,beta,v,gam) - Z_ovov(v,gam,u,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
!do gam = 1, nV
! do beta = 1, nV
! do v = 1, nO
! do u = 1, nO
! do a = 1, nV
! do i = 1, nO
! r2(u,v,beta,gam) = r2(u,v,beta,gam) &
! - K1(u,a,i,gam) * t2(i,v,beta,a) &
! - K1(v,a,i,beta) * t2(i,u,gam,a) ! P
! enddo
! enddo
! enddo
! enddo
! enddo
!enddo
!$omp parallel &
!$omp shared(nO,nV,K1,X_ovov,Z_ovov,t2) &
!$omp private(u,v,gam,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
do a = 1, nV
do i = 1, nO
do gam = 1, nV
do u = 1, nO
X_ovov(u,gam,i,a) = K1(u,a,i,gam)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do beta = 1, nV
do v = 1, nO
do a = 1, nV
do i = 1, nO
Z_ovov(i,a,v,beta) = t2(i,v,beta,a)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO*nV,nO*nV,nO*nV, &
1d0, X_ovov, size(X_ovov,1) * size(X_ovov,2), &
Y_ovov, size(Y_ovov,1) * size(Y_ovov,2), &
0d0, Z_ovov, size(Y_ovov,1) * size(Y_ovov,2))
!$omp parallel &
!$omp shared(nO,nV,r2,Z_ovov) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
do u = 1, nO
r2(u,v,beta,gam) = r2(u,v,beta,gam) - Z_ovov(u,gam,v,beta) - Z_ovov(v,beta,u,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_ovov,Y_ovov,Z_ovov)
! Change the sign for consistency with the code in spin orbitals
!$omp parallel &
!$omp shared(nO,nV,r2) &
!$omp private(i,j,a,b) &
!$omp default(none)
!$omp do collapse(3)
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
r2(i,j,a,b) = -r2(i,j,a,b)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
max_r2 = 0d0
do b = 1, nV
do a = 1, nV
do j = 1, nO
do i = 1, nO
if (dabs(r2(i,j,a,b)) > max_r2) then
max_r2 = dabs(r2(i,j,a,b))
endif
enddo
enddo
enddo
enddo
deallocate(g_occ,g_vir,J1,K1,A1)
end
#+end_src
*** Intermediates
**** A1
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_A1(nO,nV,t1,t2,tau,A1)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(in) :: tau(nO, nO, nV, nV)
double precision, intent(out) :: A1(nO, nO, nO, nO)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta
!A1 = 0d0
!do j = 1, nO
! do i = 1, nO
! do v = 1, nO
! do u = 1, nO
! A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
! do a = 1, nV
! A1(u,v,i,j) = A1(u,v,i,j) &
! + cc_space_v_ovoo(u,a,i,j) * t1(v,a) &
! + cc_space_v_vooo(a,v,i,j) * t1(u,a)
!
! do b = 1, nV
! A1(u,v,i,j) = A1(u,v,i,j) + cc_space_v_vvoo(a,b,i,j) * tau(u,v,a,b)
! enddo
! enddo
!
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_vooo(:,:,:,:), Y_oooo(:,:,:,:)
allocate(X_vooo(nV,nO,nO,nO), Y_oooo(nO,nO,nO,nO))
! A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
!$omp parallel &
!$omp shared(nO,nV,A1,cc_space_v_oooo,cc_space_v_ovoo,X_vooo) &
!$omp private(u,v,i,j) &
!$omp default(none)
!$omp do collapse(3)
do j = 1, nO
do i = 1, nO
do v = 1, nO
do u = 1, nO
A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
enddo
enddo
enddo
enddo
!$omp end do nowait
! A1(u,v,i,j) += cc_space_v_ovoo(u,a,i,j) * t1(v,a) &
!$omp do collapse(3)
do j = 1, nO
do i = 1, nO
do u = 1, nO
do a = 1, nV
X_vooo(a,u,i,j) = cc_space_v_ovoo(u,a,i,j)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N', nO, nO*nO*nO, nV, &
1d0, t1 , size(t1,1), &
X_vooo, size(X_vooo,1), &
0d0, Y_oooo, size(Y_oooo,1))
!$omp parallel &
!$omp shared(nO,nV,A1,Y_oooo) &
!$omp private(u,v,i,j) &
!$omp default(none)
!$omp do collapse(3)
do j = 1, nO
do i = 1, nO
do v = 1, nO
do u = 1, nO
A1(u,v,i,j) = A1(u,v,i,j) + Y_oooo(v,u,i,j)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_vooo,Y_oooo)
! A1(u,v,i,j) += cc_space_v_vooo(a,v,i,j) * t1(u,a)
call dgemm('N','N', nO, nO*nO*nO, nV, &
1d0, t1 , size(t1,1), &
cc_space_v_vooo, size(cc_space_v_vooo,1), &
1d0, A1 , size(A1,1))
! A1(u,v,i,j) += cc_space_v_vvoo(a,b,i,j) * tau(u,v,a,b)
call dgemm('N','N', nO*nO, nO*nO, nV*nV, &
1d0, tau , size(tau,1) * size(tau,2), &
cc_space_v_vvoo, size(cc_space_v_vvoo,1) * size(cc_space_v_vvoo,2), &
1d0, A1 , size(A1,1) * size(A1,2))
end
#+END_SRC
**** B1
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_B1(nO,nV,t1,t2,B1)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: B1(nV, nV, nV, nV)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
!B1 = 0d0
!do gam = 1, nV
! do beta = 1, nV
! do b = 1, nV
! do a = 1, nV
! B1(a,b,beta,gam) = cc_space_v_vvvv(a,b,beta,gam)
! do i = 1, nO
! B1(a,b,beta,gam) = B1(a,b,beta,gam) &
! - cc_space_v_vvvo(a,b,beta,i) * t1(i,gam) &
! - cc_space_v_vvov(a,b,i,gam) * t1(i,beta)
! enddo
!
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_vvvo(:,:,:,:), Y_vvvv(:,:,:,:)
allocate(X_vvvo(nV,nV,nV,nO), Y_vvvv(nV,nV,nV,nV))
! B1(a,b,beta,gam) = cc_space_v_vvvv(a,b,beta,gam)
!$omp parallel &
!$omp shared(nO,nV,B1,cc_space_v_vvvv,cc_space_v_vvov,X_vvvo) &
!$omp private(a,b,beta,gam) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do b = 1, nV
do a = 1, nV
B1(a,b,beta,gam) = cc_space_v_vvvv(a,b,beta,gam)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do i = 1, nO
do gam = 1, nV
do b = 1, nV
do a = 1, nV
X_vvvo(a,b,gam,i) = cc_space_v_vvov(a,b,i,gam)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
! B1(a,b,beta,gam) -= cc_space_v_vvvo(a,b,beta,i) * t1(i,gam) &
call dgemm('N','N', nV*nV*nV, nV, nO, &
-1d0, cc_space_v_vvvo, size(cc_space_v_vvvo,1) * size(cc_space_v_vvvo,2) * size(cc_space_v_vvvo,3), &
t1 , size(t1,1), &
1d0, B1 , size(B1,1) * size(B1,2) * size(B1,3))
! B1(a,b,beta,gam) -= cc_space_v_vvov(a,b,i,gam) * t1(i,beta)
call dgemm('N','N', nV*nV*nV, nV, nO, &
-1d0, X_vvvo, size(X_vvvo,1) * size(X_vvvo,2) * size(X_vvvo,3), &
t1 , size(t1,1), &
0d0, Y_vvvv, size(Y_vvvv,1) * size(Y_vvvv,2) * size(Y_vvvv,3))
!$omp parallel &
!$omp shared(nV,B1,Y_vvvv) &
!$omp private(a,b,beta,gam) &
!$omp default(none)
!$omp do collapse(3)
do gam = 1, nV
do beta = 1, nV
do b = 1, nV
do a = 1, nV
B1(a,b,beta,gam) = B1(a,b,beta,gam) + Y_vvvv(a,b,gam,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_vvvo,Y_vvvv)
end
#+END_SRC
**** g_occ
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_g_occ(nO,nV,t1,t2,H_oo,g_occ)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV), H_oo(nO, nO)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: g_occ(nO, nO)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
!g_occ = 0d0
!do i = 1, nO
! do u = 1, nO
! g_occ(u,i) = H_oo(u,i)
!
! do a = 1, nV
! g_occ(u,i) = g_occ(u,i) + cc_space_f_vo(a,i) * t1(u,a)
!
! do j = 1, nO
! g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a)
! enddo
!
! enddo
! enddo
!enddo
call dgemm('N','N',nO,nO,nV, &
1d0, t1, size(t1,1), &
cc_space_f_vo, size(cc_space_f_vo,1), &
0d0, g_occ, size(g_occ,1))
!$omp parallel &
!$omp shared(nO,nV,g_occ,H_oo, cc_space_v_ovoo,t1) &
!$omp private(i,j,a,u) &
!$omp default(none)
!$omp do
do i = 1, nO
do u = 1, nO
g_occ(u,i) = g_occ(u,i) + H_oo(u,i)
enddo
enddo
!$omp end do
!$omp do collapse(1)
do i = 1, nO
do j = 1, nO
do a = 1, nV
do u = 1, nO
g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
end
#+END_SRC
**** g_vir
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_g_vir(nO,nV,t1,t2,H_vv,g_vir)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV), H_vv(nV, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: g_vir(nV, nV)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
!g_vir = 0d0
!do beta = 1, nV
! do a = 1, nV
! g_vir(a,beta) = H_vv(a,beta)
!
! do i = 1, nO
! g_vir(a,beta) = g_vir(a,beta) - cc_space_f_vo(a,i) * t1(i,beta)
!
! do b = 1, nV
! g_vir(a,beta) = g_vir(a,beta) + (2d0 * cc_space_v_vvvo(a,b,beta,i) - cc_space_v_vvvo(b,a,beta,i)) * t1(i,b)
! enddo
!
! enddo
! enddo
!enddo
call dgemm('N','N',nV,nV,nO, &
-1d0, cc_space_f_vo , size(cc_space_f_vo,1), &
t1 , size(t1,1), &
0d0, g_vir, size(g_vir,1))
!$omp parallel &
!$omp shared(nO,nV,g_vir,H_vv, cc_space_v_vvvo,t1) &
!$omp private(i,b,a,beta) &
!$omp default(none)
!$omp do
do beta = 1, nV
do a = 1, nV
g_vir(a,beta) = g_vir(a,beta) + H_vv(a,beta)
enddo
enddo
!$omp end do
!$omp do collapse(1)
do beta = 1, nV
do i = 1, nO
do b = 1, nV
do a = 1, nV
g_vir(a,beta) = g_vir(a,beta) + (2d0 * cc_space_v_vvvo(a,b,beta,i) - cc_space_v_vvvo(b,a,beta,i)) * t1(i,b)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
end
#+END_SRC
**** J1
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(in) :: v_ovvo(nO,nV,nV,nO), v_ovoo(nO,nV,nO,nO)
double precision, intent(in) :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO)
double precision, intent(out) :: J1(nO, nV, nV, nO)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
!J1 = 0d0
!do i = 1, nO
! do beta = 1, nV
! do a = 1, nV
! do u = 1, nO
! J1(u,a,beta,i) = cc_space_v_ovvo(u,a,beta,i)
! do j = 1, nO
! J1(u,a,beta,i) = J1(u,a,beta,i) &
! - cc_space_v_ovoo(u,a,j,i) * t1(j,beta)
! enddo
! do b = 1, nV
! J1(u,a,beta,i) = J1(u,a,beta,i) &
! + cc_space_v_vvvo(b,a,beta,i) * t1(u,b)
! enddo
! do j = 1, nO
! do b = 1, nV
! J1(u,a,beta,i) = J1(u,a,beta,i) &
! - cc_space_v_vvoo(a,b,i,j) * (0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)) &
! + 0.5d0 * (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(b,a,i,j)) * t2(u,j,beta,b)
! enddo
! enddo
!
! enddo
! enddo
! enddo
!enddo
double precision, allocatable :: X_ovoo(:,:,:,:), Y_ovov(:,:,:,:)
allocate(X_ovoo(nO,nV,nO,nO),Y_ovov(nO,nV,nO,nV))
!$omp parallel &
!$omp shared(nO,nV,J1,v_ovvo,v_ovoo,X_ovoo) &
!$omp private(i,j,a,u,beta) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do beta = 1, nV
do a = 1, nV
do u = 1, nO
J1(u,a,beta,i) = v_ovvo(u,a,beta,i)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do j = 1, nO
do i = 1, nO
do a = 1, nV
do u = 1, nO
X_ovoo(u,a,i,j) = v_ovoo(u,a,j,i)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO*nV*nO,nV,nO, &
-1d0, X_ovoo, size(X_ovoo,1) * size(X_ovoo,2) * size(X_ovoo,3), &
t1 , size(t1,1), &
0d0, Y_ovov, size(Y_ovov,1) * size(Y_ovov,2) * size(Y_ovov,3))
!$omp parallel &
!$omp shared(nO,nV,J1,Y_ovov) &
!$omp private(i,beta,a,u) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do beta = 1, nV
do a = 1, nV
do u = 1, nO
J1(u,a,beta,i) = J1(u,a,beta,i) + Y_ovov(u,a,i,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_ovoo)
! v_vvvo(b,a,beta,i) * t1(u,b)
call dgemm('N','N',nO,nV*nV*nO,nV, &
1d0, t1 , size(t1,1), &
v_vvvo, size(v_vvvo,1), &
1d0, J1 , size(J1,1))
!- cc_space_v_vvoo(a,b,i,j) * (0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)) &
double precision, allocatable :: X_voov(:,:,:,:), Z_ovvo(:,:,:,:)
allocate(X_voov(nV,nO,nO,nV), Z_ovvo(nO,nV,nV,nO))
!$omp parallel &
!$omp shared(nO,nV,t2,t1,Y_ovov,X_voov,v_vvoo) &
!$omp private(i,beta,a,u,b,j) &
!$omp default(none)
!$omp do collapse(3)
do b = 1, nV
do j = 1, nO
do beta = 1, nV
do u = 1, nO
Y_ovov(u,beta,j,b) = 0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do b = 1, nV
do j = 1, nO
do i = 1, nO
do a = 1, nV
X_voov(a,i,j,b) = v_vvoo(a,b,i,j)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','T',nO*nV,nV*nO,nO*nV, &
-1d0, Y_ovov, size(Y_ovov,1) * size(Y_ovov,2), &
X_voov, size(X_voov,1) * size(X_voov,2), &
0d0, Z_ovvo, size(Z_ovvo,1) * size(Z_ovvo,2))
deallocate(X_voov)
double precision, allocatable :: X_ovvo(:,:,:,:), Y_vovo(:,:,:,:)
allocate(X_ovvo(nO,nV,nV,nO),Y_vovo(nV,nO,nV,nO))
!$omp parallel &
!$omp shared(nO,nV,J1,Z_ovvo,t2,Y_vovo,v_vvoo,X_ovvo) &
!$omp private(i,beta,a,u,j,b) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do beta = 1, nV
do a = 1, nV
do u = 1, nO
J1(u,a,beta,i) = J1(u,a,beta,i) + Z_ovvo(u,beta,a,i)
enddo
enddo
enddo
enddo
!$omp end do nowait
!+ 0.5d0 * (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(b,a,i,j)) * t2(u,j,beta,b)
!$omp do collapse(3)
do j = 1, nO
do b = 1, nV
do i = 1, nO
do a = 1, nV
Y_vovo(a,i,b,j) = 0.5d0 * (2d0 * v_vvoo(a,b,i,j) - v_vvoo(b,a,i,j))
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do j = 1, nO
do b = 1, nV
do beta = 1, nV
do u = 1, nO
X_ovvo(u,beta,b,j) = t2(u,j,beta,b)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','T',nO*nV,nV*nO,nV*nO, &
1d0, X_ovvo, size(X_ovvo,1) * size(X_ovvo,2), &
Y_vovo, size(Y_vovo,1) * size(Y_vovo,2), &
0d0, Z_ovvo, size(Z_ovvo,1) * size(Z_ovvo,2))
!$omp parallel &
!$omp shared(nO,nV,J1,Z_ovvo) &
!$omp private(i,beta,a,u) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
do beta = 1, nV
do a = 1, nV
do u = 1, nO
J1(u,a,beta,i) = J1(u,a,beta,i) + Z_ovvo(u,beta,a,i)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_ovvo,Z_ovvo,Y_ovov)
end
#+END_SRC
**** K1
#+BEGIN_SRC f90 :comments org :tangle ccsd_space_orb_sub.irp.f
subroutine compute_K1(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(in) :: v_vvoo(nV,nV,nO,nO), v_ovov(nO,nV,nO,nV)
double precision, intent(in) :: v_vvov(nV,nV,nO,nV), v_ovoo(nO,nV,nO,nO)
double precision, intent(out) :: K1(nO, nV, nO, nV)
double precision, allocatable :: X(:,:,:,:), Y(:,:,:,:), Z(:,:,:,:)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
!K1 = 0d0
!do beta = 1, nV
! do i = 1, nO
! do a = 1, nV
! do u = 1, nO
! K1(u,a,i,beta) = cc_space_v_ovov(u,a,i,beta)
! do j = 1, nO
! K1(u,a,i,beta) = K1(u,a,i,beta) &
! - cc_space_v_ovoo(u,a,i,j) * t1(j,beta)
! enddo
! do b = 1, nV
! K1(u,a,i,beta) = K1(u,a,i,beta) &
! + cc_space_v_vvov(b,a,i,beta) * t1(u,b)
! enddo
! do j = 1, nO
! do b = 1, nV
! K1(u,a,i,beta) = K1(u,a,i,beta) &
! - cc_space_v_vvoo(b,a,i,j) * (0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta))
! enddo
! enddo
!
! enddo
! enddo
! enddo
!enddo
allocate(X(nV,nO,nV,nO),Y(nO,nV,nV,nO),Z(nO,nV,nV,nO))
!$omp parallel &
!$omp shared(nO,nV,K1,X,Y,v_vvoo,v_ovov,t1,t2) &
!$omp private(i,beta,a,u,j,b) &
!$omp default(none)
!$omp do collapse(3)
do beta = 1, nV
do i = 1, nO
do a = 1, nV
do u = 1, nO
K1(u,a,i,beta) = v_ovov(u,a,i,beta)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do i = 1, nO
do a = 1, nV
do j = 1, nO
do b = 1, nV
X(b,j,a,i) = - v_vvoo(b,a,i,j)
enddo
enddo
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do j = 1, nO
do b = 1, nV
do beta = 1, nV
do u = 1, nO
Y(u,beta,b,j) = 0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO*nV*nO,nV,nO, &
-1d0, v_ovoo, size(v_ovoo,1) * size(v_ovoo,2) * size(v_ovoo,3), &
t1 , size(t1,1), &
1d0, K1 , size(K1,1) * size(K1,2) * size(K1,3))
call dgemm('N','N',nO,nV*nO*nV,nV, &
1d0, t1 , size(t1,1), &
v_vvov, size(v_vvov,1), &
1d0, K1 , size(K1,1))
! Y(u,beta,b,j) * X(b,j,a,i) = Z(u,beta,a,i)
call dgemm('N','N',nV*nO,nO*nV,nV*nO, &
1d0, Y, size(Y,1) * size(Y,2), &
X, size(X,1) * size(X,2), &
0d0, Z, size(Z,1) * size(Z,2))
!$omp parallel &
!$omp shared(nO,nV,K1,Z) &
!$omp private(i,beta,a,u) &
!$omp default(none)
!$omp do collapse(3)
do beta = 1, nV
do i = 1, nO
do a = 1, nV
do u = 1, nO
K1(u,a,i,beta) = K1(u,a,i,beta) + Z(u,beta,a,i)
enddo
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X,Y,Z)
end
#+END_SRC
View Git Blame Copy Permalink