From 0647a9db5f7173189cf38451a7032436b21f0d57 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 4 Feb 2025 14:49:39 +0100 Subject: [PATCH] Starting optimization of open-shell ccsd --- src/ccsd/ccsd_space_orb.irp.f | 5 +- src/ccsd/ccsd_spin_orb.irp.f | 2 - src/ccsd/ccsd_spin_orb_sub.irp.f | 444 ++++++++++++++--------------- src/utils_cc/mo_integrals_cc.irp.f | 16 ++ 4 files changed, 229 insertions(+), 238 deletions(-) diff --git a/src/ccsd/ccsd_space_orb.irp.f b/src/ccsd/ccsd_space_orb.irp.f index 53028ec0..91f703a0 100644 --- a/src/ccsd/ccsd_space_orb.irp.f +++ b/src/ccsd/ccsd_space_orb.irp.f @@ -1,9 +1,10 @@ -! Code - program ccsd implicit none + BEGIN_DOC + ! Closed-shell CCSD + END_DOC read_wf = .True. touch read_wf diff --git a/src/ccsd/ccsd_spin_orb.irp.f b/src/ccsd/ccsd_spin_orb.irp.f index 6f2de11c..04344fbb 100644 --- a/src/ccsd/ccsd_spin_orb.irp.f +++ b/src/ccsd/ccsd_spin_orb.irp.f @@ -1,5 +1,3 @@ -! Prog - program ccsd implicit none diff --git a/src/ccsd/ccsd_spin_orb_sub.irp.f b/src/ccsd/ccsd_spin_orb_sub.irp.f index 09d6a0fe..fe202ebf 100644 --- a/src/ccsd/ccsd_spin_orb_sub.irp.f +++ b/src/ccsd/ccsd_spin_orb_sub.irp.f @@ -11,9 +11,9 @@ subroutine run_ccsd_spin_orb double precision, allocatable :: t1(:,:), t2(:,:,:,:), tau(:,:,:,:), tau_t(:,:,:,:) double precision, allocatable :: r1(:,:), r2(:,:,:,:) double precision, allocatable :: cF_oo(:,:), cF_ov(:,:), cF_vv(:,:) - double precision, allocatable :: cW_oooo(:,:,:,:), cW_ovvo(:,:,:,:), cW_vvvv(:,:,:,:) - - double precision, allocatable :: f_oo(:,:), f_ov(:,:), f_vv(:,:), f_o(:), f_v(:) + double precision, allocatable :: cW_oooo(:,:,:,:), cW_ovvo(:,:,:,:) !, cW_vvvv(:,:,:,:) + + double precision, allocatable :: f_o(:), f_v(:) double precision, allocatable :: v_oooo(:,:,:,:), v_vooo(:,:,:,:), v_ovoo(:,:,:,:) double precision, allocatable :: v_oovo(:,:,:,:), v_ooov(:,:,:,:), v_vvoo(:,:,:,:) double precision, allocatable :: v_vovo(:,:,:,:), v_voov(:,:,:,:), v_ovvo(:,:,:,:) @@ -24,8 +24,7 @@ subroutine run_ccsd_spin_orb double precision, allocatable :: all_err(:,:), all_t(:,:) logical :: not_converged - integer, allocatable :: list_occ(:,:), list_vir(:,:) - integer :: nO,nV,nOa,nOb,nVa,nVb,nO_m,nV_m,nO_S(2),nV_S(2),n_spin(4) + integer :: nOa,nOb,nVa,nVb,nO_m,nV_m,nO_S(2),nV_S(2),n_spin(4) integer :: nb_iter, i,j,a,b double precision :: uncorr_energy, energy, max_r, max_r1, max_r2, cc, ta, tb,ti,tf,tbi,tfi integer(bit_kind) :: det(N_int,2) @@ -33,7 +32,7 @@ subroutine run_ccsd_spin_orb 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) @@ -41,107 +40,83 @@ subroutine run_ccsd_spin_orb nVa = cc_nVa !n_spin(3) nVb = cc_nVb !n_spin(4) - ! Total number of occ/vir spin orb - nO = cc_nOab !nOa + nOb - nV = cc_nVab !nVa + nVb - ! Debug - !print*,nO,nV - ! Number of occ/vir spin orb per spin nO_S = cc_nO_S !(/nOa,nOb/) nV_S = cc_nV_S !(/nVa,nVb/) ! Debug !print*,nO_S,nV_S - ! Maximal number of occ/vir + ! Maximal number of occ/vir nO_m = cc_nO_m !max(nOa, nOb) nV_m = cc_nV_m !max(nVa, nVb) ! Debug !print*,nO_m,nV_m - - allocate(list_occ(nO_m,2), list_vir(nV_m,2)) - list_occ = cc_list_occ_spin - list_vir = cc_list_vir_spin - ! Debug - !call extract_list_orb_spin(det,nO_m,nV_m,list_occ,list_vir) - !print*,list_occ(:,1) - !print*,list_occ(:,2) - !print*,list_vir(:,1) - !print*,list_vir(:,2) ! Allocation - allocate(t1(nO,nV), t2(nO,nO,nV,nV), tau(nO,nO,nV,nV), tau_t(nO,nO,nV,nV)) - allocate(r1(nO,nV), r2(nO,nO,nV,nV)) - allocate(cF_oo(nO,nO), cF_ov(nO,nV), cF_vv(nV,nV)) - allocate(cW_oooo(nO,nO,nO,nO), cW_ovvo(nO,nV,nV,nO))!, cW_vvvv(nV,nV,nV,nV)) - allocate(v_oooo(nO,nO,nO,nO)) - !allocate(v_vooo(nV,nO,nO,nO)) - allocate(v_ovoo(nO,nV,nO,nO)) - allocate(v_oovo(nO,nO,nV,nO)) - allocate(v_ooov(nO,nO,nO,nV)) - allocate(v_vvoo(nV,nV,nO,nO)) - !allocate(v_vovo(nV,nO,nV,nO)) - !allocate(v_voov(nV,nO,nO,nV)) - allocate(v_ovvo(nO,nV,nV,nO)) - allocate(v_ovov(nO,nV,nO,nV)) - allocate(v_oovv(nO,nO,nV,nV)) - !allocate(v_vvvo(nV,nV,nV,nO)) - !allocate(v_vvov(nV,nV,nO,nV)) - !allocate(v_vovv(nV,nO,nV,nV)) - !allocate(v_ovvv(nO,nV,nV,nV)) - !allocate(v_vvvv(nV,nV,nV,nV)) - allocate(f_o(nO), f_v(nV)) - allocate(f_oo(nO, nO)) - allocate(f_ov(nO, nV)) - allocate(f_vv(nV, nV)) - + allocate(t1(cc_nOab,cc_nVab), t2(cc_nOab,cc_nOab,cc_nVab,cc_nVab), tau(cc_nOab,cc_nOab,cc_nVab,cc_nVab), tau_t(cc_nOab,cc_nOab,cc_nVab,cc_nVab)) + allocate(r1(cc_nOab,cc_nVab), r2(cc_nOab,cc_nOab,cc_nVab,cc_nVab)) + allocate(cF_oo(cc_nOab,cc_nOab), cF_ov(cc_nOab,cc_nVab), cF_vv(cc_nVab,cc_nVab)) + allocate(cW_oooo(cc_nOab,cc_nOab,cc_nOab,cc_nOab), cW_ovvo(cc_nOab,cc_nVab,cc_nVab,cc_nOab))!, cW_vvvv(cc_nVab,cc_nVab,cc_nVab,cc_nVab)) + allocate(v_oooo(cc_nOab,cc_nOab,cc_nOab,cc_nOab)) + !allocate(v_vooo(cc_nVab,cc_nOab,cc_nOab,cc_nOab)) + allocate(v_ovoo(cc_nOab,cc_nVab,cc_nOab,cc_nOab)) + allocate(v_oovo(cc_nOab,cc_nOab,cc_nVab,cc_nOab)) + allocate(v_ooov(cc_nOab,cc_nOab,cc_nOab,cc_nVab)) + allocate(v_vvoo(cc_nVab,cc_nVab,cc_nOab,cc_nOab)) + !allocate(v_vovo(cc_nVab,cc_nOab,cc_nVab,cc_nOab)) + !allocate(v_voov(cc_nVab,cc_nOab,cc_nOab,cc_nVab)) + allocate(v_ovvo(cc_nOab,cc_nVab,cc_nVab,cc_nOab)) + allocate(v_ovov(cc_nOab,cc_nVab,cc_nOab,cc_nVab)) + allocate(v_oovv(cc_nOab,cc_nOab,cc_nVab,cc_nVab)) + !allocate(v_vvvo(cc_nVab,cc_nVab,cc_nVab,cc_nOab)) + !allocate(v_vvov(cc_nVab,cc_nVab,cc_nOab,cc_nVab)) + !allocate(v_vovv(cc_nVab,cc_nOab,cc_nVab,cc_nVab)) + !allocate(v_ovvv(cc_nOab,cc_nVab,cc_nVab,cc_nVab)) + !allocate(v_vvvv(cc_nVab,cc_nVab,cc_nVab,cc_nVab)) + allocate(f_o(cc_nOab), f_v(cc_nVab)) + ! Allocation for the diis 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)) + allocate(all_err(cc_nOab*cc_nVab+cc_nOab*cc_nOab*cc_nVab*cc_nVab,cc_diis_depth), all_t(cc_nOab*cc_nVab+cc_nOab*cc_nOab*cc_nVab*cc_nVab,cc_diis_depth)) all_err = 0d0 all_t = 0d0 endif - ! Fock elements - call gen_f_spin(det, nO_m,nO_m, nO_S,nO_S, list_occ,list_occ, nO,nO, f_oo) - call gen_f_spin(det, nO_m,nV_m, nO_S,nV_S, list_occ,list_vir, nO,nV, f_ov) - call gen_f_spin(det, nV_m,nV_m, nV_S,nV_S, list_vir,list_vir, nV,nV, f_vv) - ! Diag elements - do i = 1, nO - f_o(i) = f_oo(i,i) + do i = 1, cc_nOab + f_o(i) = cc_spin_f_oo(i,i) enddo - do i = 1, nV - f_v(i) = f_vv(i,i) + do i = 1, cc_nVab + f_v(i) = cc_spin_f_vv(i,i) enddo ! Bi electronic integrals from list call wall_time(ti) ! OOOO - call gen_v_spin(nO_m,nO_m,nO_m,nO_m, nO_S,nO_S,nO_S,nO_S, list_occ,list_occ,list_occ,list_occ, nO,nO,nO,nO, v_oooo) + call gen_v_spin(nO_m,nO_m,nO_m,nO_m, nO_S,nO_S,nO_S,nO_S, cc_list_occ_spin,cc_list_occ_spin,cc_list_occ_spin,cc_list_occ_spin, cc_nOab,cc_nOab,cc_nOab,cc_nOab, v_oooo) ! OOO V - !call gen_v_spin(nV_m,nO_m,nO_m,nO_m, nV_S,nO_S,nO_S,nO_S, list_vir,list_occ,list_occ,list_occ, nV,nO,nO,nO, v_vooo) - call gen_v_spin(nO_m,nV_m,nO_m,nO_m, nO_S,nV_S,nO_S,nO_S, list_occ,list_vir,list_occ,list_occ, nO,nV,nO,nO, v_ovoo) - call gen_v_spin(nO_m,nO_m,nV_m,nO_m, nO_S,nO_S,nV_S,nO_S, list_occ,list_occ,list_vir,list_occ, nO,nO,nV,nO, v_oovo) - call gen_v_spin(nO_m,nO_m,nO_m,nV_m, nO_S,nO_S,nO_S,nV_S, list_occ,list_occ,list_occ,list_vir, nO,nO,nO,nV, v_ooov) + !call gen_v_spin(nV_m,nO_m,nO_m,nO_m, nV_S,nO_S,nO_S,nO_S, cc_list_vir_spin,cc_list_occ_spin,cc_list_occ_spin,cc_list_occ_spin, cc_nVab,cc_nOab,cc_nOab,cc_nOab, v_vooo) + call gen_v_spin(nO_m,nV_m,nO_m,nO_m, nO_S,nV_S,nO_S,nO_S, cc_list_occ_spin,cc_list_vir_spin,cc_list_occ_spin,cc_list_occ_spin, cc_nOab,cc_nVab,cc_nOab,cc_nOab, v_ovoo) + call gen_v_spin(nO_m,nO_m,nV_m,nO_m, nO_S,nO_S,nV_S,nO_S, cc_list_occ_spin,cc_list_occ_spin,cc_list_vir_spin,cc_list_occ_spin, cc_nOab,cc_nOab,cc_nVab,cc_nOab, v_oovo) + call gen_v_spin(nO_m,nO_m,nO_m,nV_m, nO_S,nO_S,nO_S,nV_S, cc_list_occ_spin,cc_list_occ_spin,cc_list_occ_spin,cc_list_vir_spin, cc_nOab,cc_nOab,cc_nOab,cc_nVab, v_ooov) ! OO VV - call gen_v_spin(nV_m,nV_m,nO_m,nO_m, nV_S,nV_S,nO_S,nO_S, list_vir,list_vir,list_occ,list_occ, nV,nV,nO,nO, v_vvoo) - !call gen_v_spin(nV_m,nO_m,nV_m,nO_m, nV_S,nO_S,nV_S,nO_S, list_vir,list_occ,list_vir,list_occ, nV,nO,nV,nO, v_vovo) - !call gen_v_spin(nV_m,nO_m,nO_m,nV_m, nV_S,nO_S,nO_S,nV_S, list_vir,list_occ,list_occ,list_vir, nV,nO,nO,nV, v_voov) - call gen_v_spin(nO_m,nV_m,nV_m,nO_m, nO_S,nV_S,nV_S,nO_S, list_occ,list_vir,list_vir,list_occ, nO,nV,nV,nO, v_ovvo) - call gen_v_spin(nO_m,nV_m,nO_m,nV_m, nO_S,nV_S,nO_S,nV_S, list_occ,list_vir,list_occ,list_vir, nO,nV,nO,nV, v_ovov) - call gen_v_spin(nO_m,nO_m,nV_m,nV_m, nO_S,nO_S,nV_S,nV_S, list_occ,list_occ,list_vir,list_vir, nO,nO,nV,nV, v_oovv) + call gen_v_spin(nV_m,nV_m,nO_m,nO_m, nV_S,nV_S,nO_S,nO_S, cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin,cc_list_occ_spin, cc_nVab,cc_nVab,cc_nOab,cc_nOab, v_vvoo) + !call gen_v_spin(nV_m,nO_m,nV_m,nO_m, nV_S,nO_S,nV_S,nO_S, cc_list_vir_spin,cc_list_occ_spin,cc_list_vir_spin,cc_list_occ_spin, cc_nVab,cc_nOab,cc_nVab,cc_nOab, v_vovo) + !call gen_v_spin(nV_m,nO_m,nO_m,nV_m, nV_S,nO_S,nO_S,nV_S, cc_list_vir_spin,cc_list_occ_spin,cc_list_occ_spin,cc_list_vir_spin, cc_nVab,cc_nOab,cc_nOab,cc_nVab, v_voov) + call gen_v_spin(nO_m,nV_m,nV_m,nO_m, nO_S,nV_S,nV_S,nO_S, cc_list_occ_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, cc_nOab,cc_nVab,cc_nVab,cc_nOab, v_ovvo) + call gen_v_spin(nO_m,nV_m,nO_m,nV_m, nO_S,nV_S,nO_S,nV_S, cc_list_occ_spin,cc_list_vir_spin,cc_list_occ_spin,cc_list_vir_spin, cc_nOab,cc_nVab,cc_nOab,cc_nVab, v_ovov) + call gen_v_spin(nO_m,nO_m,nV_m,nV_m, nO_S,nO_S,nV_S,nV_S, cc_list_occ_spin,cc_list_occ_spin,cc_list_vir_spin,cc_list_vir_spin, cc_nOab,cc_nOab,cc_nVab,cc_nVab, v_oovv) ! O VVV - !call gen_v_spin(nV_m,nV_m,nV_m,nO_m, nV_S,nV_S,nV_S,nO_S, list_vir,list_vir,list_vir,list_occ, nV,nV,nV,nO, v_vvvo) - !call gen_v_spin(nV_m,nV_m,nO_m,nV_m, nV_S,nV_S,nO_S,nV_S, list_vir,list_vir,list_occ,list_vir, nV,nV,nO,nV, v_vvov) - !call gen_v_spin(nV_m,nO_m,nV_m,nV_m, nV_S,nO_S,nV_S,nV_S, list_vir,list_occ,list_vir,list_vir, nV,nO,nV,nV, v_vovv) - !call gen_v_spin(nO_m,nV_m,nV_m,nV_m, nO_S,nV_S,nV_S,nV_S, list_occ,list_vir,list_vir,list_vir, nO,nV,nV,nV, v_ovvv) + !call gen_v_spin(nV_m,nV_m,nV_m,nO_m, nV_S,nV_S,nV_S,nO_S, cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, cc_nVab,cc_nVab,cc_nVab,cc_nOab, v_vvvo) + !call gen_v_spin(nV_m,nV_m,nO_m,nV_m, nV_S,nV_S,nO_S,nV_S, cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin,cc_list_vir_spin, cc_nVab,cc_nVab,cc_nOab,cc_nVab, v_vvov) + !call gen_v_spin(nV_m,nO_m,nV_m,nV_m, nV_S,nO_S,nV_S,nV_S, cc_list_vir_spin,cc_list_occ_spin,cc_list_vir_spin,cc_list_vir_spin, cc_nVab,cc_nOab,cc_nVab,cc_nVab, v_vovv) + !call gen_v_spin(nO_m,nV_m,nV_m,nV_m, nO_S,nV_S,nV_S,nV_S, cc_list_occ_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin, cc_nOab,cc_nVab,cc_nVab,cc_nVab, v_ovvv) ! VVVV - !call gen_v_spin(nV_m,nV_m,nV_m,nV_m, nV_S,nV_S,nV_S,nV_S, list_vir,list_vir,list_vir,list_vir, nV,nV,nV,nV, v_vvvv) + !call gen_v_spin(nV_m,nV_m,nV_m,nV_m, nV_S,nV_S,nV_S,nV_S, cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin, cc_nVab,cc_nVab,cc_nVab,cc_nVab, v_vvvv) call wall_time(tf) if (cc_dev) then print*,'Load bi elec int:',tf-ti,'s' @@ -149,11 +124,11 @@ subroutine run_ccsd_spin_orb ! Init of T t1 = 0d0 - call guess_t1(nO,nV,f_o,f_v,f_ov,t1) - call guess_t2(nO,nV,f_o,f_v,v_oovv,t2) - call compute_tau_spin(nO,nV,t1,t2,tau) - call compute_tau_t_spin(nO,nV,t1,t2,tau_t) - + call guess_t1(cc_nOab,cc_nVab,f_o,f_v,cc_spin_f_ov,t1) + call guess_t2(cc_nOab,cc_nVab,f_o,f_v,v_oovv,t2) + call compute_tau_spin(cc_nOab,cc_nVab,t1,t2,tau) + call compute_tau_t_spin(cc_nOab,cc_nVab,t1,t2,tau_t) + ! Loop init nb_iter = 0 not_converged = .True. @@ -164,9 +139,9 @@ subroutine run_ccsd_spin_orb call det_energy(det,uncorr_energy) print*,'Det energy', uncorr_energy - call ccsd_energy_spin(nO,nV,t1,t2,F_ov,v_oovv,energy) + call ccsd_energy_spin(cc_nOab,cc_nVab,t1,t2,cc_spin_F_ov,v_oovv,energy) print*,'guess energy', uncorr_energy+energy, energy - + write(*,'(A77)') ' -----------------------------------------------------------------------------' write(*,'(A77)') ' | It. | E(CCSD) (Ha) | Correlation (Ha) | Conv. T1 | Conv. T2 |' write(*,'(A77)') ' -----------------------------------------------------------------------------' @@ -179,18 +154,18 @@ subroutine run_ccsd_spin_orb ! Intermediates call wall_time(tbi) call wall_time(ti) - call compute_cF_oo(nO,nV,t1,tau_t,F_oo,F_ov,v_ooov,v_oovv,cF_oo) - call compute_cF_ov(nO,nV,t1,F_ov,v_oovv,cF_ov) - call compute_cF_vv(nO,nV,t1,tau_t,F_ov,F_vv,v_oovv,cF_vv) + call compute_cF_oo(cc_nOab,cc_nVab,t1,tau_t,cc_spin_F_oo,cc_spin_F_ov,v_ooov,v_oovv,cF_oo) + call compute_cF_ov(cc_nOab,cc_nVab,t1,cc_spin_F_ov,v_oovv,cF_ov) + call compute_cF_vv(cc_nOab,cc_nVab,t1,tau_t,cc_spin_F_ov,cc_spin_F_vv,v_oovv,cF_vv) call wall_time(tf) if (cc_dev) then print*,'Compute cFs:',tf-ti,'s' endif - + call wall_time(ti) - call compute_cW_oooo(nO,nV,t1,t2,tau,v_oooo,v_ooov,v_oovv,cW_oooo) - call compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) - !call compute_cW_vvvv(nO,nV,t1,t2,tau,v_vvvv,v_vovv,v_oovv,cW_vvvv) + call compute_cW_oooo(cc_nOab,cc_nVab,t1,t2,tau,v_oooo,v_ooov,v_oovv,cW_oooo) + call compute_cW_ovvo(cc_nOab,cc_nVab,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) + !call compute_cW_vvvv(cc_nOab,cc_nVab,t1,t2,tau,v_vvvv,v_vovv,v_oovv,cW_vvvv) call wall_time(tf) if (cc_dev) then print*,'Compute cFs:',tf-ti,'s' @@ -198,13 +173,13 @@ subroutine run_ccsd_spin_orb ! Residuals call wall_time(ti) - call compute_r1_spin(nO,nV,t1,t2,f_o,f_v,F_ov,cF_oo,cF_ov,cF_vv,v_oovo,v_ovov,r1) + call compute_r1_spin(cc_nOab,cc_nVab,t1,t2,f_o,f_v,cc_spin_F_ov,cF_oo,cF_ov,cF_vv,v_oovo,v_ovov,r1) call wall_time(tf) if (cc_dev) then print*,'Compute r1:',tf-ti,'s' endif call wall_time(ti) - call compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ovvo,v_ovoo,v_oovv,v_ovvo,r2) + call compute_r2_spin(cc_nOab,cc_nVab,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ovvo,v_ovoo,v_oovv,v_ovvo,r2) call wall_time(tf) if (cc_dev) then print*,'Compute r2:',tf-ti,'s' @@ -218,29 +193,29 @@ subroutine run_ccsd_spin_orb call wall_time(ti) ! 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,f_o,f_v,r1,r2,t1,t2,all_err,all_t) + !call update_t_ccsd(cc_nOab,cc_nVab,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2) + !call update_t_ccsd_diis(cc_nOab,cc_nVab,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err1,all_err2,all_t1,all_t2) + call update_t_ccsd_diis_v3(cc_nOab,cc_nVab,nb_iter,f_o,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,f_o,f_v,r1,t1) - call update_t2(nO,nV,f_o,f_v,r2,t2) + call update_t1(cc_nOab,cc_nVab,f_o,f_v,r1,t1) + call update_t2(cc_nOab,cc_nVab,f_o,f_v,r2,t2) else print*,'Unkonw cc_method_method: '//cc_update_method endif - call compute_tau_spin(nO,nV,t1,t2,tau) - call compute_tau_t_spin(nO,nV,t1,t2,tau_t) + call compute_tau_spin(cc_nOab,cc_nVab,t1,t2,tau) + call compute_tau_t_spin(cc_nOab,cc_nVab,t1,t2,tau_t) call wall_time(tf) if (cc_dev) then print*,'Update:',tf-ti,'s' endif ! Print - call ccsd_energy_spin(nO,nV,t1,t2,F_ov,v_oovv,energy) + call ccsd_energy_spin(cc_nOab,cc_nVab,t1,t2,cc_spin_F_ov,v_oovv,energy) call wall_time(tfi) - + write(*,'(A3,I6,A3,F18.12,A3,F16.12,A3,ES10.2,A3,ES10.2,A2)') ' | ',nb_iter,' | ', & uncorr_energy+energy,' | ', energy,' | ', max_r1,' | ', max_r2,' |' if (cc_dev) then @@ -270,8 +245,8 @@ subroutine run_ccsd_spin_orb print*,'' if (write_amplitudes) then - call write_t1(nO,nV,t1) - call write_t2(nO,nV,t2) + call write_t1(cc_nOab,cc_nVab,t1) + call write_t2(cc_nOab,cc_nVab,t2) call ezfio_set_utils_cc_io_amplitudes('Read') endif @@ -286,20 +261,20 @@ subroutine run_ccsd_spin_orb deallocate(v_oooo) deallocate(v_ovoo,v_oovo) deallocate(v_ovvo,v_ovov,v_oovv) - + double precision :: t_corr t_corr = 0.d0 if (cc_par_t .and. elec_alpha_num +elec_beta_num > 2) then print*,'CCSD(T) calculation...' call wall_time(ta) - !allocate(v_vvvo(nV,nV,nV,nO)) + !allocate(v_vvvo(cc_nVab,cc_nVab,cc_nVab,cc_nOab)) !call gen_v_spin(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, & ! cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, & ! cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, & - ! nV,nV,nV,nO, v_vvvo) + ! cc_nVab,cc_nVab,cc_nVab,cc_nOab, v_vvvo) - !call ccsd_par_t_spin(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,v_vvvo,t_corr) - call ccsd_par_t_spin_v2(nO,nV,t1,t2,f_o,f_v,f_ov,v_ooov,v_vvoo,t_corr) + !call ccsd_par_t_spin(cc_nOab,cc_nVab,t1,t2,f_o,f_v,cc_spin_f_ov,v_ooov,v_vvoo,v_vvvo,t_corr) + call ccsd_par_t_spin_v2(cc_nOab,cc_nVab,t1,t2,f_o,f_v,cc_spin_f_ov,v_ooov,v_vvoo,t_corr) !print*,'Working on it...' !call abort call wall_time(tb) @@ -313,12 +288,12 @@ subroutine run_ccsd_spin_orb endif call save_energy(uncorr_energy + energy, t_corr) - - deallocate(f_oo,f_ov,f_vv,f_o,f_v) + + deallocate(f_o,f_v) deallocate(v_ooov,v_vvoo,t1,t2) !deallocate(v_ovvv,v_vvvo,v_vovv) !deallocate(v_vvvv) - + end ! Energy @@ -354,7 +329,7 @@ subroutine ccsd_energy_spin(nO,nV,t1,t2,Fov,v_oovv,energy) do j=1,nO do a=1,nV do b=1,nV - energy = energy & + energy = energy & + 0.5d0 * v_oovv(i,j,a,b) * t1(i,a) * t1(j,b) & + 0.25d0 * v_oovv(i,j,a,b) * t2(i,j,a,b) end do @@ -375,7 +350,7 @@ subroutine compute_tau_spin(nO,nV,t1,t2,tau) double precision,intent(in) :: t2(nO,nO,nV,nV) double precision,intent(out) :: tau(nO,nO,nV,nV) - + integer :: i,j,k,l integer :: a,b,c,d @@ -463,7 +438,7 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov !$OMP v_ovov,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,f,m,n) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(1) do a=1,nV do i=1,nO @@ -494,7 +469,7 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov 1d0, t1 , size(t1,1), & cF_vv, size(cF_vv,1), & 1d0, r1 , size(r1,1)) - + !do a=1,nV ! do i=1,nO ! do m=1,nO @@ -531,7 +506,7 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov !$OMP SHARED(r1,t1,t2,X_vovf,v_ovvf,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,f,m,n) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) !do f = 1, nV @@ -546,28 +521,28 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov !enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('N','T', nO, nV, nO*nV, & -0.5d0, t2(1,1,1,f), size(t2,1), & X_vovf, size(X_vovf,1), & 1d0 , r1 , size(r1,1)) enddo - + !call dgemm('N','T', nO, nV, nO*nV*nV, & ! -0.5d0, t2 , size(t2,1), & ! X_vovv, size(X_vovv,1), & ! 1d0 , r1 , size(r1,1)) - + deallocate(X_vovf) !deallocate(X_vovv) allocate(X_oovv(nO,nO,nV,nV)) - + !$OMP PARALLEL & !$OMP SHARED(r1,t1,t2,X_oovv, & !$OMP f_o,f_v,v_oovo,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,f,m,n) & !$OMP DEFAULT(NONE) - + !do a=1,nV ! do i=1,nO ! do e=1,nV @@ -579,7 +554,7 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov ! end do ! end do !end do - + !$OMP DO collapse(3) do a = 1, nV do e = 1, nV @@ -592,12 +567,12 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('T','N', nO, nV, nO*nO*nV, & -0.5d0, v_oovo, size(v_oovo,1) * size(v_oovo,2) * size(v_oovo,3), & X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3), & 1d0 , r1 , size(r1,1)) - + !$OMP PARALLEL & !$OMP SHARED(r1,t1,X_oovv,f_o,f_v,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,f,m,n) & @@ -610,7 +585,7 @@ subroutine compute_r1_spin(nO,nV,t1,t2,f_o,f_v,Fov,cF_oo,cF_ov,cF_vv,v_oovo,v_ov enddo !$OMP END DO !$OMP END PARALLEL - + deallocate(X_oovv) end @@ -684,7 +659,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ !$OMP SHARED(r2,v_oovv,X_oovv,nO,nV) & !$OMP PRIVATE(i,j,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -697,7 +672,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ end do !$OMP END DO !$OMP END PARALLEL - + !deallocate(X_oovv) !do b=1,nV @@ -726,25 +701,25 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ 0.5d0, t2 , size(t2,1) * size(t2,2) * size(t2,3), & A_vv , size(A_vv,1), & 0d0 , X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3)) - + !$OMP PARALLEL & !$OMP SHARED(r2,v_oovv,X_oovv,nO,nV) & !$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) - X_oovv(i,j,a,b) + X_oovv(i,j,b,a) + r2(i,j,a,b) = r2(i,j,a,b) - X_oovv(i,j,a,b) + X_oovv(i,j,b,a) end do end do end do end do !$OMP END DO !$OMP END PARALLEL - + deallocate(A_vv)!,X_oovv) !do b=1,nV @@ -766,7 +741,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ !$OMP SHARED(t2,v_oovv,X_oovv,nO,nV) & !$OMP PRIVATE(i,m,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -789,13 +764,13 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ !$OMP SHARED(r2,v_oovv,Y_oovv,nO,nV) & !$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) - Y_oovv(j,i,a,b) + Y_oovv(i,j,a,b) + r2(i,j,a,b) = r2(i,j,a,b) - Y_oovv(j,i,a,b) + Y_oovv(i,j,a,b) end do end do end do @@ -821,17 +796,17 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ ! end do !end do allocate(A_oo(nO,nO),B_oovv(nO,nO,nV,nV))!,X_oovv(nO,nO,nV,nV)) - + call dgemm('N','T', nO, nO, nV, & 1d0, t1 , size(t1,1), & cF_ov, size(cF_ov,1), & 0d0, A_oo , size(A_oo,1)) - + !$OMP PARALLEL & !$OMP SHARED(t2,B_oovv,nO,nV) & !$OMP PRIVATE(i,m,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b = 1, nV do a = 1, nV @@ -844,17 +819,17 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('N','N', nO, nO*nV*nV, nO, & 0.5d0, A_oo, size(A_oo,1), & B_oovv, size(B_oovv,1), & 0d0 , X_oovv, size(X_oovv,1)) - + !$OMP PARALLEL & !$OMP SHARED(r2,X_oovv,nO,nV) & !$OMP PRIVATE(i,j,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -888,7 +863,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ 0.5d0, cW_oooo, size(cW_oooo,1) * size(cW_oooo,2), & tau , size(tau,1) * size(tau,2), & 1d0 , r2 , size(r2,1) * size(r2,2)) - + !do b=1,nV ! do a=1,nV ! do j=1,nO @@ -908,6 +883,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ ! 0.5d0, tau , size(tau,1) * size(tau,2), & ! cW_vvvv, size(cW_vvvv,1) * size(cW_vvvv,2), & ! 1d0 , r2 , size(r2,1) * size(r2,2)) + double precision :: ti,tf call wall_time(ti) call use_cW_vvvf(nO,nV,t1,t2,tau,v_oovv,r2) @@ -915,7 +891,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ if (cc_dev) then print*,'cW_vvvv:',tf-ti,'s' endif - + !do b=1,nV ! do a=1,nV ! do j=1,nO @@ -923,7 +899,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ ! do e=1,nV ! do m=1,nO - ! r2(i,j,a,b) = r2(i,j,a,b) & + ! r2(i,j,a,b) = r2(i,j,a,b) & ! + t2(i,m,a,e)*cW_ovvo(m,b,e,j) & ! - t2(j,m,a,e)*cW_ovvo(m,b,e,i) & ! - t2(i,m,b,e)*cW_ovvo(m,a,e,j) & @@ -944,7 +920,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ !$OMP SHARED(t2,A_ovov,B_ovvo,cW_ovvo,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do a = 1, nV do i = 1, nO @@ -961,24 +937,24 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ do b = 1, nV do e = 1, nV do m = 1, nO - B_ovvo(m,e,b,j) = cW_ovvo(m,b,e,j) + B_ovvo(m,e,b,j) = cW_ovvo(m,b,e,j) enddo enddo enddo enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('T','N', nO*nV, nV*nO, nO*nV, & 1d0, A_ovov, size(A_ovov,1) * size(A_ovov,2), & B_ovvo, size(B_ovvo,1) * size(B_ovvo,2), & 0d0, X_ovvo, size(X_ovvo,1) * size(X_ovvo,2)) - + !$OMP PARALLEL & !$OMP SHARED(r2,X_ovvo,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b = 1, nV do a = 1, nV @@ -992,15 +968,15 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ enddo !$OMP END DO !$OMP END PARALLEL - + deallocate(A_ovov,B_ovvo,X_ovvo) allocate(A_vvoo(nV,nV,nO,nO), B_ovoo(nO,nV,nO,nO), C_ovov(nO,nV,nO,nV)) - + !$OMP PARALLEL & !$OMP SHARED(A_vvoo,v_ovvo,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do m = 1, nO do j = 1, nO @@ -1013,22 +989,22 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('N','N', nO, nV*nO*nO, nV, & 1d0, t1 , size(t1,1), & A_vvoo, size(A_vvoo,1), & 0d0, B_ovoo, size(B_ovoo,1)) - + call dgemm('N','N', nO*nV*nO, nV, nO, & 1d0, B_ovoo, size(B_ovoo,1) * size(B_ovoo,2) * size(B_ovoo,3), & t1 , size(t1,1), & 0d0, C_ovov, size(C_ovov,1) * size(C_ovov,2) * size(C_ovov,3)) - + !$OMP PARALLEL & !$OMP SHARED(r2,C_ovov,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -1042,9 +1018,9 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ end do !$OMP END DO !$OMP END PARALLEL - + deallocate(A_vvoo, B_ovoo, C_ovov) - + !do b=1,nV ! do a=1,nV ! do j=1,nO @@ -1065,12 +1041,12 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ call gen_v_spin_3idx_i_kl(cc_nV_m,cc_nV_m,cc_nV_m,cc_nO_m, b, cc_nV_S,cc_nV_S,cc_nV_S,cc_nO_S, & cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_occ_spin, & nV,nV,nO, v_vbvo) - + !$OMP PARALLEL & !$OMP SHARED(b,A_vbov,v_vbvo,nO,nV) & !$OMP PRIVATE(i,j,a,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(2) do e = 1, nV do j = 1, nO @@ -1093,12 +1069,12 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ ! 1d0, A_vvov, size(A_vvov,1) * size(A_vvov,2) * size(A_vvov,3), & ! t1 , size(t1,1), & ! 0d0, X_vvoo, size(X_vvoo,1) * size(X_vvoo,2) * size(X_vvoo,3)) - + !$OMP PARALLEL & !$OMP SHARED(b,r2,X_vboo,nO,nV) & !$OMP PRIVATE(i,j,a,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(2) !do b = 1, nV do a = 1, nV @@ -1113,7 +1089,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ !$OMP END DO !$OMP END PARALLEL enddo - + !deallocate(A_vvov)!,X_vvoo) deallocate(A_vbov, X_vboo, v_vbvo) allocate(X_vvoo(nV,nV,nO,nO)) @@ -1132,7 +1108,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ ! end do !end do !allocate(X_vvoo(nV,nV,nO,nO)) - + call dgemm('T','N', nV, nV*nO*nO, nO, & 1d0, t1 , size(t1,1), & v_ovoo, size(v_ovoo,1), & @@ -1142,7 +1118,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ !$OMP SHARED(r2,X_vvoo,f_o,f_v,t2,nO,nV) & !$OMP PRIVATE(i,j,a,b,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -1154,7 +1130,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ end do end do !$OMP END DO - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -1167,7 +1143,7 @@ subroutine compute_r2_spin(nO,nV,t1,t2,tau,f_o,f_v,cF_oo,cF_ov,cF_vv,cW_oooo,cW_ end do !$OMP END DO !$OMP END PARALLEL - + deallocate(X_vvoo) end @@ -1182,16 +1158,16 @@ subroutine use_cF_oo(nO,nV,t1,t2,tau_t,F_oo,F_ov,v_ooov,v_oovv,r1,r2) double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV), tau_t(nO,nO,nV,nV) double precision, intent(in) :: F_oo(nO,nV), F_ov(nO,nV) double precision, intent(in) :: v_ooov(nO,nO,nO,nV), v_oovv(nO,nO,nV,nV) - + double precision, intent(inout) :: r1(nO,nV), r2(nO,nO,nV,nV) - + double precision, allocatable :: cF_oo(:,:), X_oovv(:,:,:,:),Y_oovv(:,:,:,:) integer :: i,j,m,a,b allocate(cF_oo(nO,nO)) - + call compute_cF_oo(nO,nV,t1,tau_t,F_oo,F_ov,v_ooov,v_oovv,cF_oo) - + !do a=1,nV ! do i=1,nO ! do m=1,nO @@ -1218,13 +1194,13 @@ subroutine use_cF_oo(nO,nV,t1,t2,tau_t,F_oo,F_ov,v_ooov,v_oovv,r1,r2) ! end do ! end do !end do - + allocate(Y_oovv(nO,nO,nV,nV),X_oovv(nO,nO,nV,nV)) !$OMP PARALLEL & !$OMP SHARED(t2,v_oovv,X_oovv,nO,nV) & !$OMP PRIVATE(i,m,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -1247,20 +1223,20 @@ subroutine use_cF_oo(nO,nV,t1,t2,tau_t,F_oo,F_ov,v_ooov,v_oovv,r1,r2) !$OMP SHARED(r2,v_oovv,Y_oovv,nO,nV) & !$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) - Y_oovv(j,i,a,b) + Y_oovv(i,j,a,b) + r2(i,j,a,b) = r2(i,j,a,b) - Y_oovv(j,i,a,b) + Y_oovv(i,j,a,b) end do end do end do end do !$OMP END DO !$OMP END PARALLEL - + deallocate(cF_oo,X_oovv,Y_oovv) end @@ -1274,7 +1250,7 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) integer, intent(in) :: nO,nV double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV) double precision, intent(in) :: F_ov(nO,nV), v_oovv(nO,nO,nV,nV) - + double precision, intent(inout) :: r1(nO,nV), r2(nO,nO,nV,nV) double precision, allocatable :: cF_ov(:,:), A_oo(:,:), A_vv(:,:) @@ -1282,14 +1258,14 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) integer :: i,j,a,b,e,m allocate(cF_ov(nO,nV)) - + call compute_cF_ov(nO,nV,t1,F_ov,v_oovv,cF_ov) !$OMP PARALLEL & !$OMP SHARED(r1,t2,cF_ov,nO,nV) & !$OMP PRIVATE(i,a,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(1) do a=1,nV do i=1,nO @@ -1334,22 +1310,22 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) !$OMP SHARED(nO,nV,r2,X_oovv) & !$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) - X_oovv(i,j,a,b) + X_oovv(i,j,b,a) + r2(i,j,a,b) = r2(i,j,a,b) - X_oovv(i,j,a,b) + X_oovv(i,j,b,a) end do end do end do end do !$OMP END DO !$OMP END PARALLEL - + deallocate(A_vv) - + !do b=1,nV ! do a=1,nV ! do j=1,nO @@ -1367,17 +1343,17 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) ! end do !end do allocate(A_oo(nO,nO),B_oovv(nO,nO,nV,nV))!,X_oovv(nO,nO,nV,nV)) - + call dgemm('N','T', nO, nO, nV, & 1d0, t1 , size(t1,1), & cF_ov, size(cF_ov,1), & 0d0, A_oo , size(A_oo,1)) - + !$OMP PARALLEL & !$OMP SHARED(t2,B_oovv,nO,nV) & !$OMP PRIVATE(i,m,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b = 1, nV do a = 1, nV @@ -1390,7 +1366,7 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('N','N', nO, nO*nV*nV, nO, & 0.5d0, A_oo, size(A_oo,1), & B_oovv, size(B_oovv,1), & @@ -1400,7 +1376,7 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) !$OMP SHARED(r2,X_oovv,nO,nV) & !$OMP PRIVATE(i,j,a,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b=1,nV do a=1,nV @@ -1413,9 +1389,9 @@ subroutine use_cF_ov(nO,nV,t1,t2,F_ov,v_oovv,r1,r2) end do !$OMP END DO !$OMP END PARALLEL - + deallocate(cF_ov,A_oo,B_oovv,X_oovv) - + end ! Use cF_vv @@ -1426,18 +1402,18 @@ subroutine use_cF_vv(nO,nV,t1,t2,r1,r2) integer, intent(in) :: nO,nV double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV) - + double precision, intent(inout) :: r1(nO,nV), r2(nO,nO,nV,nV) double precision, allocatable :: cF_vv(:,:) integer :: i,j,a,b,e,m allocate(cF_vv(nV,nV)) - + !call compute_cF_vv(nO,nV,t1,tau_t,F_ov,F_vv,v_oovv,v_ovvv,cF_vv) deallocate(cF_vv) - + end ! Use cW_vvvd @@ -1450,7 +1426,7 @@ subroutine use_cW_vvvf(nO,nV,t1,t2,tau,v_oovv,r2) double precision, intent(in) :: t1(nO,nV), t2(nO,nO,nV,nV), tau(nO,nO,nV,nV) double precision, intent(in) :: v_oovv(nO,nO,nV,nV) !double precision, intent(in) :: v_vovv(nV,nO,nV,nV) - + double precision, intent(inout) :: r2(nO,nO,nV,nV) double precision, allocatable :: cW_vvvf(:,:,:), v_vvvf(:,:,:), tau_f(:,:,:), v_vovf(:,:,:) @@ -1460,7 +1436,7 @@ subroutine use_cW_vvvf(nO,nV,t1,t2,tau,v_oovv,r2) allocate(cW_vvvf(nV,nV,nV),v_vvvf(nV,nV,nV),tau_f(nO,nO,nV),v_vovf(nV,nO,nV)) !PROVIDE cc_nVab - + !do b=1,nV ! do a=1,nV ! do j=1,nO @@ -1476,14 +1452,14 @@ subroutine use_cW_vvvf(nO,nV,t1,t2,tau,v_oovv,r2) ! end do ! end do !end do - + do f = 1, nV call wall_time(ti) !$OMP PARALLEL & !$OMP SHARED(tau,tau_f,f,nO,nV) & !$OMP PRIVATE(i,j,e) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(2) do e = 1, nV do j = 1, nO @@ -1515,7 +1491,7 @@ subroutine use_cW_vvvf(nO,nV,t1,t2,tau,v_oovv,r2) if (cc_dev .and. f == 1) then print*,'vovf', tf-ti endif - + call wall_time(ti) call compute_cW_vvvf(nO,nV,t1,t2,tau,f,v_vvvf,v_vovf,v_oovv,cW_vvvf) call wall_time(tf) @@ -1535,7 +1511,7 @@ subroutine use_cW_vvvf(nO,nV,t1,t2,tau,v_oovv,r2) enddo deallocate(cW_vvvf,v_vvvf,v_vovf) - + end ! cF_oo @@ -1562,7 +1538,7 @@ subroutine compute_cF_oo(nO,nV,t1,tau_t,Foo,Fov,v_ooov,v_oovv,cF_oo) !$OMP SHARED(cF_oo,Foo,t1,v_ooov,nO,nV) & !$OMP PRIVATE(i,m,n,e) & !$OMP DEFAULT(NONE) - + !do i=1,nO ! do m=1,nO ! cF_oo(m,i) = (1d0 - Kronecker_delta(m,i))*Foo(m,i) @@ -1580,7 +1556,7 @@ subroutine compute_cF_oo(nO,nV,t1,tau_t,Foo,Fov,v_ooov,v_oovv,cF_oo) cF_oo(i,i) = 0d0 end do !$OMP END DO - + do e=1,nV do n=1,nO !$OMP DO collapse(1) @@ -1620,8 +1596,8 @@ subroutine compute_cF_oo(nO,nV,t1,tau_t,Foo,Fov,v_ooov,v_oovv,cF_oo) call dgemm('N','T', nO, nO, nO*nV*nV, & 0.5d0, v_oovv, size(v_oovv,1), & tau_t , size(tau_t,1), & - 1d0 , cF_oo , size(cF_oo,1)) - + 1d0 , cF_oo , size(cF_oo,1)) + end ! cF_ov @@ -1643,7 +1619,7 @@ subroutine compute_cF_ov(nO,nV,t1,Fov,v_oovv,cF_ov) !$OMP SHARED(cF_ov,Fov,t1,v_oovv,nO,nV) & !$OMP PRIVATE(i,a,m,n,e,f) & !$OMP DEFAULT(NONE) - + !cF_ov = Fov !$OMP DO collapse(1) @@ -1659,7 +1635,7 @@ subroutine compute_cF_ov(nO,nV,t1,Fov,v_oovv,cF_ov) end do !$OMP END DO !$OMP END PARALLEL - + end ! cF_vv @@ -1677,7 +1653,7 @@ subroutine compute_cF_vv(nO,nV,t1,tau_t,Fov,Fvv,v_oovv,cF_vv) !double precision,intent(in) :: v_ovvv(nO,nV,nV,nV) double precision,intent(out) :: cF_vv(nV,nV) - + double precision, allocatable :: v_ovfv(:,:,:),X_ovfv(:,:,:) integer :: i,j,m,n integer :: a,b,e,f @@ -1699,7 +1675,7 @@ subroutine compute_cF_vv(nO,nV,t1,tau_t,Fov,Fvv,v_oovv,cF_vv) enddo !$OMP END DO !$OMP END PARALLEL - + !do e=1,nV ! do a=1,nV ! do m=1,nO @@ -1711,7 +1687,7 @@ subroutine compute_cF_vv(nO,nV,t1,tau_t,Fov,Fvv,v_oovv,cF_vv) -0.5d0, t1 , size(t1,1), & Fov , size(Fov,1), & 1d0 , cF_vv, size(cF_vv,1)) - + !do e=1,nV ! do a=1,nV ! do m=1,nO @@ -1791,7 +1767,7 @@ subroutine compute_cW_oooo(nO,nV,t1,t2,tau,v_oooo,v_ooov,v_oovv,cW_oooo) integer :: a,b,e,f double precision, allocatable :: X_oooo(:,:,:,:) - ! oooo block + ! oooo block !cW_oooo = v_oooo @@ -1809,7 +1785,7 @@ subroutine compute_cW_oooo(nO,nV,t1,t2,tau,v_oooo,v_ooov,v_oovv,cW_oooo) ! end do !end do allocate(X_oooo(nO,nO,nO,nO)) - + call dgemm('N','T', nO*nO*nO, nO, nV, & 1d0, v_ooov, size(v_ooov,1) * size(v_ooov,2) * size(v_ooov,3), & t1 , size(t1,1), & @@ -1830,14 +1806,14 @@ subroutine compute_cW_oooo(nO,nV,t1,t2,tau,v_oooo,v_ooov,v_oovv,cW_oooo) end do !$OMP END DO !$OMP END PARALLEL - + deallocate(X_oooo) - + !do m=1,nO ! do n=1,nO ! do i=1,nO ! do j=1,nO - ! + ! ! do e=1,nV ! do f=1,nV ! cW_oooo(m,n,i,j) = cW_oooo(m,n,i,j) + 0.25d0*tau(i,j,e,f)*v_oovv(m,n,e,f) @@ -1853,7 +1829,7 @@ subroutine compute_cW_oooo(nO,nV,t1,t2,tau,v_oooo,v_ooov,v_oovv,cW_oooo) 0.25d0, v_oovv , size(v_oovv,1) * size(v_oovv,2), & tau , size(tau,1) * size(tau,2), & 1.d0 , cW_oooo, size(cW_oooo,1) * size(cW_oooo,2)) - + end ! cW_ovvo @@ -1913,7 +1889,7 @@ subroutine compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) call gen_v_spin_3idx_ij_l(cc_nO_m,cc_nV_m,cc_nV_m,cc_nV_m, e, cc_nO_S,cc_nV_S,cc_nV_S,cc_nV_S, & cc_list_occ_spin,cc_list_vir_spin,cc_list_vir_spin,cc_list_vir_spin, & nO,nV,nV, v_ovev) - + call dgemm('N','T', nO*nV, nO, nV, & 1.d0, v_ovev , size(v_ovev,1) * size(v_ovev,2), & t1 , size(t1,1), & @@ -1950,14 +1926,14 @@ subroutine compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) ! end do ! end do !end do - + allocate(A_oovo(nO,nO,nV,nO), B_vovo(nV,nO,nV,nO)) - + !$OMP PARALLEL & !$OMP SHARED(A_oovo,v_oovo,nO,nV) & !$OMP PRIVATE(j,e,m,n) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do j=1,nO do e=1,nV @@ -1970,17 +1946,17 @@ subroutine compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) end do !$OMP END DO !$OMP END PARALLEL - + call dgemm('T','N', nV, nO*nV*nO, nO, & 1d0, t1 , size(t1,1), & A_oovo, size(A_oovo,1), & 0d0, B_vovo, size(B_vovo,1)) - + !$OMP PARALLEL & !$OMP SHARED(cW_ovvo,B_vovo,nO,nV) & !$OMP PRIVATE(j,e,m,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do j=1,nO do e=1,nV @@ -2015,7 +1991,7 @@ subroutine compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) !$OMP SHARED(nO,nV,A_voov,B_voov,v_oovv,t2,t1) & !$OMP PRIVATE(f,n,m,e,j,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do b = 1, nV do j = 1, nO @@ -2039,19 +2015,19 @@ subroutine compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) enddo !$OMP END DO !$OMP END PARALLEL - + call dgemm('T','N', nO*nV, nV*nO, nV*nO, & 1d0, A_voov, size(A_voov,1) * size(A_voov,2), & B_voov, size(B_voov,1) * size(B_voov,2), & 0d0, C_ovov, size(C_ovov,1) * size(C_ovov,2)) - + deallocate(A_voov,B_voov) !$OMP PARALLEL & !$OMP SHARED(cW_ovvo,C_ovov,nO,nV) & !$OMP PRIVATE(j,e,m,b) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(3) do j = 1, nO do e = 1, nV @@ -2064,7 +2040,7 @@ subroutine compute_cW_ovvo(nO,nV,t1,t2,tau,v_ovvo,v_oovo,v_oovv,cW_ovvo) enddo !$OMP END DO !$OMP END PARALLEL - + deallocate(C_ovov) end @@ -2072,7 +2048,7 @@ end ! cW_vvvv subroutine compute_cW_vvvv(nO,nV,t1,t2,tau,v_vvvv,v_vovv,v_oovv,cW_vvvv) - + implicit none integer,intent(in) :: nO,nV @@ -2154,14 +2130,14 @@ subroutine compute_cW_vvvv(nO,nV,t1,t2,tau,v_vvvv,v_vovv,v_oovv,cW_vvvv) end do !$OMP END DO !$OMP END PARALLEL - + deallocate(A_ovvv,B_vvvv) !do a=1,nV ! do b=1,nV ! do e=1,nV ! do f=1,nV - ! + ! ! do m=1,nO ! do n=1,nO ! cW_vvvv(a,b,e,f) = cW_vvvv(a,b,e,f) + 0.25d0*tau(m,n,a,b)*v_oovv(m,n,e,f) @@ -2182,7 +2158,7 @@ end ! cW_vvvf subroutine compute_cW_vvvf(nO,nV,t1,t2,tau,f,v_vvvf,v_vovf,v_oovv,cW_vvvf) - + implicit none integer,intent(in) :: nO,nV,f @@ -2207,7 +2183,7 @@ subroutine compute_cW_vvvf(nO,nV,t1,t2,tau,f,v_vvvf,v_vovf,v_oovv,cW_vvvf) !$OMP SHARED(nO,nV,cW_vvvf,A_ovvf,v_vovf,v_vvvf,f) & !$OMP PRIVATE(a,b,c,d,e,m) & !$OMP DEFAULT(NONE) - + !$OMP DO collapse(2) do c = 1, nV do b = 1, nV @@ -2248,7 +2224,7 @@ subroutine compute_cW_vvvf(nO,nV,t1,t2,tau,f,v_vvvf,v_vovf,v_oovv,cW_vvvf) 1d0, t1 , size(t1,1), & A_ovvf, size(A_ovvf,1), & 0d0, B_vvvf, size(B_vvvf,1)) - + !$OMP PARALLEL & !$OMP SHARED(nO,nV,cW_vvvf,B_vvvf,v_oovf,v_oovv,f) & !$OMP PRIVATE(a,b,c,d,e,m,n) & @@ -2264,14 +2240,14 @@ subroutine compute_cW_vvvf(nO,nV,t1,t2,tau,f,v_vvvf,v_vovf,v_oovv,cW_vvvf) end do end do !$OMP END DO NOWAIT - + !deallocate(A_ovvf,B_vvvf) !do a=1,nV ! do b=1,nV ! do e=1,nV ! do f=1,nV - ! + ! ! do m=1,nO ! do n=1,nO ! cW_vvvv(a,b,e,f) = cW_vvvv(a,b,e,f) + 0.25d0*tau(m,n,a,b)*v_oovv(m,n,e,f) @@ -2292,13 +2268,13 @@ subroutine compute_cW_vvvf(nO,nV,t1,t2,tau,f,v_vvvf,v_vovf,v_oovv,cW_vvvf) enddo enddo !$OMP END DO - !$OMP END PARALLEL - + !$OMP END PARALLEL + call dgemm('T','N', nV*nV, nV, nO*nO, & 0.25d0, tau , size(tau,1) * size(tau,2), & v_oovf , size(v_oovf,1) * size(v_oovf,2), & 1.d0 , cW_vvvf, size(cW_vvvf,1) * size(cW_vvvf,2)) - + deallocate(v_oovf) deallocate(A_ovvf,B_vvvf) diff --git a/src/utils_cc/mo_integrals_cc.irp.f b/src/utils_cc/mo_integrals_cc.irp.f index 6f21c316..eebc84ca 100644 --- a/src/utils_cc/mo_integrals_cc.irp.f +++ b/src/utils_cc/mo_integrals_cc.irp.f @@ -1006,6 +1006,22 @@ BEGIN_PROVIDER [double precision, cc_space_f_v, (cc_nVa)] END_PROVIDER + +BEGIN_PROVIDER [ double precision, cc_spin_f_oo, (cc_nOab, cc_nOab)] + implicit none + call gen_f_spin(psi_det(1,1,cc_ref), cc_nO_m, cc_nO_m, cc_nO_S, cc_nO_S, cc_list_occ_spin, cc_list_occ_spin, cc_nOab, cc_nOab, cc_spin_f_oo) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, cc_spin_f_ov, (cc_nOab, cc_nVab)] + implicit none + call gen_f_spin(psi_det(1,1,cc_ref), cc_nO_m, cc_nV_m, cc_nO_S, cc_nV_S, cc_list_occ_spin, cc_list_vir_spin, cc_nOab, cc_nVab, cc_spin_f_ov) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, cc_spin_f_vv, (cc_nVab, cc_nVab)] + implicit none + call gen_f_spin(psi_det(1,1,cc_ref), cc_nV_m, cc_nV_m, cc_nV_S, cc_nV_S, cc_list_vir_spin, cc_list_vir_spin, cc_nVab, cc_nVab, cc_spin_f_vv) +END_PROVIDER + ! Shift subroutine shift_idx_spin(s,n_S,shift)