From 1cbb3d2fe66ada34311798b6965ca47b5d96734f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 2 Sep 2020 18:50:30 +0200 Subject: [PATCH 01/14] Replaced overlap bby pt2 matrix --- src/ao_one_e_ints/pseudopot.f90 | 1 + src/cipsi/selection.irp.f | 23 +++++++++++------------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90 index aaa456ba..6c8e5c83 100644 --- a/src/ao_one_e_ints/pseudopot.f90 +++ b/src/ao_one_e_ints/pseudopot.f90 @@ -1862,6 +1862,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) qk = dble(q) two_qkmp1 = 2.d0*(qk+mk)+1.d0 do k=0,q-1 + ! possible FPE here. To be checked s_q_k = two_qkmp1*qk*inverses(k)*s_q_k ! if (s_q_k < 1.d-32) then ! s_q_k = 0.d0 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 071d9294..c325020f 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -780,24 +780,23 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif enddo - ! Gram-Schmidt using input overlap matrix - do istate=1,N_states - do jstate=1,istate-1 - if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle - coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate) - enddo - enddo +! ! Gram-Schmidt using input overlap matrix +! do istate=1,N_states +! do jstate=1,istate-1 +! if ( (pt2_overlap(jstate,istate) == 0.d0).or.(pt2_overlap(jstate,jstate) == 0.d0) ) cycle +! coef(istate) = coef(istate) - pt2_overlap(jstate,istate)/pt2_overlap(jstate,jstate) * coef(jstate) +! enddo +! enddo do istate=1, N_states - do jstate=1,N_states - pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) - enddo - enddo - do istate=1,N_states alpha_h_psi = mat(istate, p1, p2) e_pert = coef(istate) * alpha_h_psi + do jstate=1,N_states + pt2_data % overlap(jstate,istate) += coef(jstate) * alpha_h_psi + enddo + pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi pt2_data % pt2(istate) += e_pert From 82e4299ad6e72742a7388f4e3a889c663c7a67b0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Sep 2020 18:58:07 +0200 Subject: [PATCH 02/14] Introduced nxn diagonalization for PT2 --- src/cipsi/selection.irp.f | 89 ++++++++++++++++++++++++++-------- src/utils/linear_algebra.irp.f | 2 +- 2 files changed, 70 insertions(+), 21 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index e599737c..024dd435 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -1,4 +1,3 @@ - use bitmasks BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] @@ -144,6 +143,23 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] END_PROVIDER +BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! If true, the excitation is banned in the selection. Useful with local MOs. + END_DOC + banned_excitation = .False. + integer :: i,j + double precision :: buffer(mo_num) + do j=1,mo_num + call get_mo_two_e_integrals_exch_ii(j,j,mo_num,buffer,mo_integrals_map) + buffer = dabs(buffer) + do i=1,mo_num + banned_excitation(i,j) = buffer(i) < 1.d-15 + enddo + enddo +END_PROVIDER + subroutine get_mask_phase(det1, pm, Nint) use bitmasks @@ -288,6 +304,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp + PROVIDE banned_excitation monoAdo = .true. monoBdo = .true. @@ -607,7 +624,8 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d h2 = hole_list(i2,s2) call apply_hole(pmask, s2,h2, mask, ok, N_int) - banned = .false. + banned(:,:,1) = banned_excitation(:,:) + banned(:,:,2) = banned_excitation(:,:) do j=1,mo_num bannedOrb(j, 1) = .true. bannedOrb(j, 2) = .true. @@ -674,7 +692,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d logical :: ok integer :: s1, s2, p1, p2, ib, j, istate, jstate integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert(N_states), coef(N_states), X(N_states) + double precision :: e_pert(N_states), coef(N_states) double precision :: delta_E, val, Hii, w, tmp, alpha_h_psi double precision, external :: diag_H_mat_elem_fock double precision :: E_shift @@ -769,10 +787,16 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! call occ_pattern_of_det(det,occ,N_int) ! call occ_pattern_to_dets_size(occ,n,elec_alpha_num,N_int) + e_pert = 0.d0 + coef = 0.d0 + logical :: do_diag + do_diag = .False. do istate=1,N_states delta_E = E0(istate) - Hii + E_shift alpha_h_psi = mat(istate, p1, p2) + if (alpha_h_psi == 0.d0) cycle + val = alpha_h_psi + alpha_h_psi tmp = dsqrt(delta_E * delta_E + val * val) if (delta_E < 0.d0) then @@ -784,13 +808,38 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d else coef(istate) = alpha_h_psi / delta_E endif - if (e_pert(istate) < 0.d0) then - X(istate) = -dsqrt(-e_pert(istate)) - else - X(istate) = dsqrt(e_pert(istate)) - endif enddo + do_diag = maxval(dabs(coef)) > 0.001d0 + + double precision :: eigvalues(N_states+1) + double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2) + integer :: iwork(3+5*(N_states+1)), info, k ,n + + if (do_diag) then + double precision :: pt2_matrix(N_states+1,N_states+1) + pt2_matrix(N_states+1,N_states+1) = Hii+E_shift + do istate=1,N_states + pt2_matrix(:,istate) = 0.d0 + pt2_matrix(istate,istate) = E0(istate) + pt2_matrix(istate,N_states+1) = mat(istate,p1,p2) + pt2_matrix(N_states+1,istate) = mat(istate,p1,p2) + enddo + + call DSYEVD( 'V', 'U', N_states+1, pt2_matrix, N_states+1, eigvalues, & + work, size(work), iwork, size(iwork), info ) + if (info /= 0) then + print *, 'error in '//irp_here + stop -1 + endif + pt2_matrix = dabs(pt2_matrix) + iwork = maxloc(pt2_matrix,1) + do k=1,N_states + e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k)) + enddo + endif + + ! ! Gram-Schmidt using input overlap matrix ! do istate=1,N_states ! do jstate=1,istate-1 @@ -835,25 +884,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case(5) ! Variance selection w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate) - do jstate=1,N_states - if (istate == jstate) cycle - w = w + alpha_h_psi*mat(jstate,p1,p2) * s_weight(istate,jstate) - enddo +! do jstate=1,N_states +! if (istate == jstate) cycle +! w = w + dabs(alpha_h_psi*mat(jstate,p1,p2)) * s_weight(istate,jstate) +! enddo case(6) w = w - coef(istate) * coef(istate) * s_weight(istate,istate) - do jstate=1,N_states - if (istate == jstate) cycle - w = w + coef(istate)*coef(jstate) * s_weight(istate,jstate) - enddo +! do jstate=1,N_states +! if (istate == jstate) cycle +! w = w + dabs(coef(istate)*coef(jstate)) * s_weight(istate,jstate) +! enddo case default ! Energy selection w = w + e_pert(istate) * s_weight(istate,istate) - do jstate=1,N_states - if (istate == jstate) cycle - w = w - dabs(X(istate))*X(jstate) * s_weight(istate,jstate) - enddo +! do jstate=1,N_states +! if (istate == jstate) cycle +! w = w + dabs(X(istate)*X(jstate)) * s_weight(istate,jstate) +! enddo end select end do diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index a8dea97a..cf43edd5 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -19,7 +19,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) double precision,allocatable :: A_tmp(:,:) allocate (A_tmp(LDA,n)) - A_tmp = A + A_tmp(:,:) = A(:,:) ! Find optimal size for temp arrays allocate(work(1)) From 161517e0ea4d520fd4ee9a5db7ba3aa9fb3d8d61 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Sep 2020 19:21:19 +0200 Subject: [PATCH 03/14] Using max instead of sum in pt2 --- src/cipsi/selection.irp.f | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 024dd435..27e47f03 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -810,7 +810,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif enddo - do_diag = maxval(dabs(coef)) > 0.001d0 + do_diag = sum(dabs(coef)) > 0.001d0 double precision :: eigvalues(N_states+1) double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2) @@ -883,14 +883,16 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case(5) ! Variance selection - w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate) +! w = w - alpha_h_psi * alpha_h_psi * s_weight(istate,istate) + w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) ! do jstate=1,N_states ! if (istate == jstate) cycle ! w = w + dabs(alpha_h_psi*mat(jstate,p1,p2)) * s_weight(istate,jstate) ! enddo case(6) - w = w - coef(istate) * coef(istate) * s_weight(istate,istate) +! w = w - coef(istate) * coef(istate) * s_weight(istate,istate) + w = min(w,- coef(istate) * coef(istate) * s_weight(istate,istate)) ! do jstate=1,N_states ! if (istate == jstate) cycle ! w = w + dabs(coef(istate)*coef(jstate)) * s_weight(istate,jstate) @@ -898,7 +900,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case default ! Energy selection - w = w + e_pert(istate) * s_weight(istate,istate) +! w = w + e_pert(istate) * s_weight(istate,istate) + w = min(w, e_pert(istate) * s_weight(istate,istate)) ! do jstate=1,N_states ! if (istate == jstate) cycle ! w = w + dabs(X(istate)*X(jstate)) * s_weight(istate,jstate) From e8388681812b76f2a302094a3283f70ab6b5722e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 25 Sep 2020 15:14:56 +0200 Subject: [PATCH 04/14] Fixed Schwartz and added banned-excitations --- ocaml/Units.ml | 4 +-- src/ao_two_e_ints/two_e_integrals.irp.f | 2 +- src/cipsi/selection.irp.f | 18 -------------- src/mo_two_e_ints/map_integrals.irp.f | 33 +++++++++++++++++++++++++ 4 files changed, 36 insertions(+), 21 deletions(-) diff --git a/ocaml/Units.ml b/ocaml/Units.ml index ab0e944c..7a7cb39e 100644 --- a/ocaml/Units.ml +++ b/ocaml/Units.ml @@ -4,8 +4,8 @@ type units = | Angstrom ;; -let angstrom_to_bohr = 1. /. 0.52917721092 -let bohr_to_angstrom = 0.52917721092 +let angstrom_to_bohr = 1. /. 0.52917721067121 +let bohr_to_angstrom = 0.52917721067121 ;; diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index b6e959d7..8c6b3875 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -436,7 +436,7 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ] !$OMP SCHEDULE(dynamic) do i=1,ao_num do k=1,i - ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,k,i,k)) + ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,i,k,k)) ao_two_e_integral_schwartz(k,i) = ao_two_e_integral_schwartz(i,k) enddo enddo diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 27e47f03..f1eb3425 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -143,24 +143,6 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] END_PROVIDER -BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] - implicit none - BEGIN_DOC - ! If true, the excitation is banned in the selection. Useful with local MOs. - END_DOC - banned_excitation = .False. - integer :: i,j - double precision :: buffer(mo_num) - do j=1,mo_num - call get_mo_two_e_integrals_exch_ii(j,j,mo_num,buffer,mo_integrals_map) - buffer = dabs(buffer) - do i=1,mo_num - banned_excitation(i,j) = buffer(i) < 1.d-15 - enddo - enddo -END_PROVIDER - - subroutine get_mask_phase(det1, pm, Nint) use bitmasks implicit none diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 661add2e..2221e0c9 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -99,6 +99,10 @@ double precision function get_two_e_integral(i,j,k,l,map) type(map_type), intent(inout) :: map real(integral_kind) :: tmp PROVIDE mo_two_e_integrals_in_map mo_integrals_cache + if (banned_excitation(i,k) .or. banned_excitation(j,l)) then + get_two_e_integral = 0.d0 + return + endif ii = l-mo_integrals_cache_min ii = ior(ii, k-mo_integrals_cache_min) ii = ior(ii, j-mo_integrals_cache_min) @@ -159,6 +163,11 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) ! return !DEBUG + out_val(1:sze) = 0.d0 + if (banned_excitation(j,l)) then + return + endif + ii0 = l-mo_integrals_cache_min ii0 = ior(ii0, k-mo_integrals_cache_min) ii0 = ior(ii0, j-mo_integrals_cache_min) @@ -172,6 +181,7 @@ subroutine get_mo_two_e_integrals(j,k,l,sze,out_val,map) q = q+shiftr(s*s-s,1) do i=1,sze + if (banned_excitation(i,k)) cycle ii = ior(ii0, i-mo_integrals_cache_min) if (iand(ii, -128) == 0) then ii_8 = ior( shiftl(ii0_8,7), int(i,8)-mo_integrals_cache_min_8) @@ -272,6 +282,29 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) end +BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] + implicit none + use map_module + BEGIN_DOC + ! If true, the excitation is banned in the selection. Useful with local MOs. + END_DOC + banned_excitation = .False. + integer :: i,j + integer(key_kind) :: idx + double precision :: tmp +! double precision :: buffer(mo_num) + do j=1,mo_num + do i=1,j-1 + call two_e_integrals_index(i,j,j,i,idx) + !DIR$ FORCEINLINE + call map_get(mo_integrals_map,idx,tmp) + banned_excitation(i,j) = dabs(tmp) < 1.d-15 + banned_excitation(j,i) = banned_excitation(i,j) + enddo + enddo +END_PROVIDER + + integer*8 function get_mo_map_size() implicit none From 88b798fe22d0c3d5c718686d11772a281c31308c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 30 Sep 2020 20:48:27 +0200 Subject: [PATCH 05/14] Added randomized SVD routines --- src/utils/linear_algebra.irp.f | 155 ++++++++++++++++++++++++++++++++- 1 file changed, 152 insertions(+), 3 deletions(-) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index cf43edd5..deb4f4ca 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -11,7 +11,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) integer, intent(in) :: LDA, LDU, LDVt, m, n double precision, intent(in) :: A(LDA,n) - double precision, intent(out) :: U(LDU,m) + double precision, intent(out) :: U(LDU,min(m,n)) double precision,intent(out) :: Vt(LDVt,n) double precision,intent(out) :: D(min(m,n)) double precision,allocatable :: work(:) @@ -24,14 +24,14 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) ! Find optimal size for temp arrays allocate(work(1)) lwork = -1 - call dgesvd('A','A', m, n, A_tmp, LDA, & + call dgesvd('S','S', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 lwork = max(int(work(1)), 5*MIN(M,N)) deallocate(work) allocate(work(lwork)) - call dgesvd('A','A', m, n, A_tmp, LDA, & + call dgesvd('S','S', m, n, A_tmp, LDA, & D, U, LDU, Vt, LDVt, work, lwork, info) deallocate(work,A_tmp) @@ -42,6 +42,128 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) end +subroutine eigSVD(A,LDA,U,LDU,D,Vt,LDVt,m,n) + implicit none + BEGIN_DOC +! Algorithm 3 of https://arxiv.org/pdf/1810.06860.pdf +! +! A(m,n) = U(m,n) D(n) Vt(n,n) with m>n + END_DOC + integer, intent(in) :: LDA, LDU, LDVt, m, n + double precision, intent(in) :: A(LDA,n) + double precision, intent(out) :: U(LDU,n) + double precision,intent(out) :: Vt(LDVt,n) + double precision,intent(out) :: D(n) + + integer :: i,j,k + + if (m= 0.d0) then + exit + endif + D(k) = dsqrt(-D(k)) + call dscal(n, 1.d0/D(k), V(1,k), 1) + enddo + D(k:n) = 0.d0 + k=k-1 + call dgemm('N','N',m,n,k,1.d0,A,size(A,1),V,size(V,1),0.d0,U,size(U,1)) + +end + + +subroutine randomized_svd(A,LDA,U,LDU,D,Vt,LDVt,m,n,q,r) + implicit none + include 'constants.include.F' + BEGIN_DOC +! Randomized SVD: rank r, q power iterations +! +! 1. Sample column space of A with P: Z = A.P where P is random with r+p columns. +! +! 2. Power iterations : Z <- X . (Xt.Z) +! +! 3. Z = Q.R +! +! 4. Compute SVD on projected Qt.X = U' . S. Vt +! +! 5. U = Q U' + END_DOC + + integer, intent(in) :: LDA, LDU, LDVt, m, n, q, r + double precision, intent(in) :: A(LDA,n) + double precision, intent(out) :: U(LDU,r) + double precision,intent(out) :: Vt(LDVt,r) + double precision,intent(out) :: D(r) + integer :: i, j, k + + double precision,allocatable :: Z(:,:), P(:,:), Y(:,:), UY(:,:) + double precision :: r1,r2 + allocate(P(n,r), Z(m,r)) + + ! P is a normal random matrix (n,r) + do k=1,r + do i=1,n + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + P(i,k) = r1*dcos(r2) + enddo + enddo + + ! Z(m,r) = A(m,n).P(n,r) + call dgemm('N','N',m,r,n,1.d0,A,size(A,1),P,size(P,1),0.d0,Z,size(Z,1)) + + ! Power iterations + do k=1,q + ! P(n,r) = At(n,m).Z(m,r) + call dgemm('T','N',n,r,m,1.d0,A,size(A,1),Z,size(Z,1),0.d0,P,size(P,1)) + ! Z(m,r) = A(m,n).P(n,r) + call dgemm('N','N',m,r,n,1.d0,A,size(A,1),P,size(P,1),0.d0,Z,size(Z,1)) + enddo + + deallocate(P) + + ! QR factorization of Z + call ortho_svd(Z,size(Z,1),m,r) + + allocate(Y(r,n), UY(r,r)) + ! Y(r,n) = Zt(r,m).A(m,n) + call dgemm('T','N',r,n,m,1.d0,Z,size(Z,1),A,size(A,1),0.d0,Y,size(Y,1)) + + ! SVD of Y + call svd(Y,size(Y,1),UY,size(UY,1),D,Vt,size(Vt,1),r,n) + deallocate(Y) + + ! U(m,r) = Z(m,r).UY(r,r) + call dgemm('N','N',m,r,r,1.d0,Z,size(Z,1),UY,size(UY,1),0.d0,U,size(U,1)) + deallocate(UY,Z) + +end subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n) implicit none @@ -807,6 +929,33 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m,cutoff) end +subroutine ortho_svd(A,LDA,m,n) + implicit none + BEGIN_DOC + ! Orthogonalization via fast SVD + ! + ! A : matrix to orthogonalize + ! + ! LDA : leftmost dimension of A + ! + ! m : Number of rows of A + ! + ! n : Number of columns of A + ! + END_DOC + integer, intent(in) :: m,n, LDA + double precision, intent(inout) :: A(LDA,n) + if (m < n) then + call ortho_qr(A,LDA,m,n) + endif + double precision, allocatable :: U(:,:), D(:), Vt(:,:) + allocate(U(m,n), D(n), Vt(n,n)) + call SVD(A,LDA,U,size(U,1),D,Vt,size(Vt,1),m,n) + A(1:m,1:n) = U(1:m,1:n) + deallocate(U,D, Vt) + +end + subroutine ortho_qr(A,LDA,m,n) implicit none BEGIN_DOC From 28cb4a516702fe0da6b2a7e63b0f1fe8d5cb9b10 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Oct 2020 11:35:45 +0200 Subject: [PATCH 06/14] Change default in davidson --- src/davidson/EZFIO.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/davidson/EZFIO.cfg b/src/davidson/EZFIO.cfg index f820730c..4343cc1e 100644 --- a/src/davidson/EZFIO.cfg +++ b/src/davidson/EZFIO.cfg @@ -8,7 +8,7 @@ default: 1.e-10 type: logical doc: Thresholds of Davidson's algorithm is set to E(rPT2)*threshold_davidson_from_pt2 interface: ezfio,provider,ocaml -default: true +default: false [n_states_diag] type: States_number From b1484089386b4321527a4456ff0fc1fa6d737f35 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Oct 2020 16:58:07 +0200 Subject: [PATCH 07/14] Fixed maxloc --- src/cipsi/selection.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index f1eb3425..00ac915e 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -815,7 +815,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d stop -1 endif pt2_matrix = dabs(pt2_matrix) - iwork = maxloc(pt2_matrix,1) + iwork(1:N_states) = maxloc(pt2_matrix,1) do k=1,N_states e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k)) enddo From 862a1804a4e7e3081771f8e0b9007ba750712a6f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 24 Oct 2020 14:11:04 +0200 Subject: [PATCH 08/14] Fixed Maxloc --- src/cipsi/selection.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 00ac915e..08595aea 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -815,7 +815,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d stop -1 endif pt2_matrix = dabs(pt2_matrix) - iwork(1:N_states) = maxloc(pt2_matrix,1) + iwork(1:N_states+1) = maxloc(pt2_matrix,DIM=1) do k=1,N_states e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k)) enddo From cbc33201a27d86a1d139ad6fba739a33d0decc3b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 25 Oct 2020 13:45:34 +0100 Subject: [PATCH 09/14] Test files --- src/fci/40.fci.bats | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index 28faa758..f86c7488 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -59,7 +59,7 @@ function run_stoch() { @test "HCO" { # 12.2868s qp set_file hco.ezfio - run -113.389298880564 3.e-4 100000 + run -113.390292448857 3.e-4 100000 } @test "H2O2" { # 12.9214s @@ -175,7 +175,7 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file cu_nh3_4_2plus.ezfio qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]" - run -1862.98614665139 3.e-04 100000 + run -1862.97568806589 3.e-04 100000 } @test "HCN" { # 20.3273s From e33bd5060dbf138353b1c7e6b49b25cf0353dbc6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Oct 2020 13:42:36 +0100 Subject: [PATCH 10/14] Moved CASSCF into qp_plugins_scemama --- src/casscf/50.casscf.bats | 49 --- src/casscf/EZFIO.cfg | 31 -- src/casscf/MORALITY | 1 - src/casscf/NEED | 4 - src/casscf/README.rst | 5 - src/casscf/bavard.irp.f | 6 - src/casscf/bielec.irp.f | 155 -------- src/casscf/bielec_natorb.irp.f | 369 ------------------ src/casscf/casscf.irp.f | 58 --- src/casscf/class.irp.f | 12 - src/casscf/densities.irp.f | 63 --- src/casscf/det_manip.irp.f | 125 ------ src/casscf/driver_optorb.irp.f | 3 - src/casscf/get_energy.irp.f | 51 --- src/casscf/grad_old.irp.f | 74 ---- src/casscf/gradient.irp.f | 171 --------- src/casscf/hessian.irp.f | 656 -------------------------------- src/casscf/mcscf_fock.irp.f | 80 ---- src/casscf/natorb.irp.f | 231 ----------- src/casscf/neworbs.irp.f | 221 ----------- src/casscf/reorder_orb.irp.f | 70 ---- src/casscf/save_energy.irp.f | 9 - src/casscf/superci_dm.irp.f | 207 ---------- src/casscf/swap_orb.irp.f | 132 ------- src/casscf/test_pert_2rdm.irp.f | 29 -- src/casscf/tot_en.irp.f | 101 ----- 26 files changed, 2913 deletions(-) delete mode 100644 src/casscf/50.casscf.bats delete mode 100644 src/casscf/EZFIO.cfg delete mode 100644 src/casscf/MORALITY delete mode 100644 src/casscf/NEED delete mode 100644 src/casscf/README.rst delete mode 100644 src/casscf/bavard.irp.f delete mode 100644 src/casscf/bielec.irp.f delete mode 100644 src/casscf/bielec_natorb.irp.f delete mode 100644 src/casscf/casscf.irp.f delete mode 100644 src/casscf/class.irp.f delete mode 100644 src/casscf/densities.irp.f delete mode 100644 src/casscf/det_manip.irp.f delete mode 100644 src/casscf/driver_optorb.irp.f delete mode 100644 src/casscf/get_energy.irp.f delete mode 100644 src/casscf/grad_old.irp.f delete mode 100644 src/casscf/gradient.irp.f delete mode 100644 src/casscf/hessian.irp.f delete mode 100644 src/casscf/mcscf_fock.irp.f delete mode 100644 src/casscf/natorb.irp.f delete mode 100644 src/casscf/neworbs.irp.f delete mode 100644 src/casscf/reorder_orb.irp.f delete mode 100644 src/casscf/save_energy.irp.f delete mode 100644 src/casscf/superci_dm.irp.f delete mode 100644 src/casscf/swap_orb.irp.f delete mode 100644 src/casscf/test_pert_2rdm.irp.f delete mode 100644 src/casscf/tot_en.irp.f diff --git a/src/casscf/50.casscf.bats b/src/casscf/50.casscf.bats deleted file mode 100644 index a0db725d..00000000 --- a/src/casscf/50.casscf.bats +++ /dev/null @@ -1,49 +0,0 @@ -#!/usr/bin/env bats - -source $QP_ROOT/tests/bats/common.bats.sh -source $QP_ROOT/quantum_package.rc - - -function run_stoch() { - thresh=$2 - test_exe casscf || skip - qp set perturbation do_pt2 True - qp set determinants n_det_max $3 - qp set davidson threshold_davidson 1.e-10 - qp set davidson n_states_diag 4 - qp run casscf | tee casscf.out - energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)" - eq $energy1 $1 $thresh -} - -@test "F2" { # 18.0198s - rm -rf f2_casscf - qp_create_ezfio -b aug-cc-pvdz ../input/f2.zmt -o f2_casscf - qp set_file f2_casscf - qp run scf - qp set_mo_class --core="[1-6,8-9]" --act="[7,10]" --virt="[11-46]" - run_stoch -198.773366970 1.e-4 100000 -} - -@test "N2" { # 18.0198s - rm -rf n2_casscf - qp_create_ezfio -b aug-cc-pvdz ../input/n2.xyz -o n2_casscf - qp set_file n2_casscf - qp run scf - qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]" - run_stoch -109.0961643162 1.e-4 100000 -} - -@test "N2_stretched" { - rm -rf n2_stretched_casscf - qp_create_ezfio -b aug-cc-pvdz -m 7 ../input/n2_stretched.xyz -o n2_stretched_casscf - qp set_file n2_stretched_casscf - qp run scf | tee scf.out - qp set_mo_class --core="[1-4]" --act="[5-10]" --virt="[11-46]" - qp set electrons elec_alpha_num 7 - qp set electrons elec_beta_num 7 - run_stoch -108.7860471300 1.e-4 100000 -# - -} - diff --git a/src/casscf/EZFIO.cfg b/src/casscf/EZFIO.cfg deleted file mode 100644 index 4e4d3d3a..00000000 --- a/src/casscf/EZFIO.cfg +++ /dev/null @@ -1,31 +0,0 @@ -[energy] -type: double precision -doc: Calculated Selected |FCI| energy -interface: ezfio -size: (determinants.n_states) - -[energy_pt2] -type: double precision -doc: Calculated |FCI| energy + |PT2| -interface: ezfio -size: (determinants.n_states) - -[cisd_guess] -type: logical -doc: If true, the CASSCF starts with a CISD wave function -interface: ezfio,provider,ocaml -default: True - -[state_following_casscf] -type: logical -doc: If |true|, the CASSCF will try to follow the guess CI vector and orbitals -interface: ezfio,provider,ocaml -default: False - - -[level_shift_casscf] -type: Positive_float -doc: Energy shift on the virtual MOs to improve SCF convergence -interface: ezfio,provider,ocaml -default: 0.005 - diff --git a/src/casscf/MORALITY b/src/casscf/MORALITY deleted file mode 100644 index 9701a647..00000000 --- a/src/casscf/MORALITY +++ /dev/null @@ -1 +0,0 @@ -the CASCF can be obtained if a proper guess is given to the WF part diff --git a/src/casscf/NEED b/src/casscf/NEED deleted file mode 100644 index d9da718e..00000000 --- a/src/casscf/NEED +++ /dev/null @@ -1,4 +0,0 @@ -cipsi -selectors_full -generators_cas -two_body_rdm diff --git a/src/casscf/README.rst b/src/casscf/README.rst deleted file mode 100644 index 08bfd95b..00000000 --- a/src/casscf/README.rst +++ /dev/null @@ -1,5 +0,0 @@ -====== -casscf -====== - -|CASSCF| program with the CIPSI algorithm. diff --git a/src/casscf/bavard.irp.f b/src/casscf/bavard.irp.f deleted file mode 100644 index 463c3ea4..00000000 --- a/src/casscf/bavard.irp.f +++ /dev/null @@ -1,6 +0,0 @@ -! -*- F90 -*- -BEGIN_PROVIDER [logical, bavard] -! bavard=.true. - bavard=.false. -END_PROVIDER - diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f deleted file mode 100644 index 0a44f994..00000000 --- a/src/casscf/bielec.irp.f +++ /dev/null @@ -1,155 +0,0 @@ -BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] - BEGIN_DOC - ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active - ! indices are unshifted orbital numbers - END_DOC - implicit none - integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 - real*8 :: mo_two_e_integral - - bielec_PQxx(:,:,:,:) = 0.d0 - PROVIDE mo_two_e_integrals_in_map - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,ii,j,jj,i3,j3) & - !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx, & - !$OMP n_act_orb,mo_integrals_map,list_act) - - !$OMP DO - do i=1,n_core_inact_orb - ii=list_core_inact(i) - do j=i,n_core_inact_orb - jj=list_core_inact(j) - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j),mo_integrals_map) - bielec_PQxx(:,:,j,i)=bielec_PQxx(:,:,i,j) - end do - do j=1,n_act_orb - jj=list_act(j) - j3=j+n_core_inact_orb - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i,j3),mo_integrals_map) - bielec_PQxx(:,:,j3,i)=bielec_PQxx(:,:,i,j3) - end do - end do - !$OMP END DO - - - !$OMP DO - do i=1,n_act_orb - ii=list_act(i) - i3=i+n_core_inact_orb - do j=i,n_act_orb - jj=list_act(j) - j3=j+n_core_inact_orb - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,bielec_PQxx(1,1,i3,j3),mo_integrals_map) - bielec_PQxx(:,:,j3,i3)=bielec_PQxx(:,:,i3,j3) - end do - end do - !$OMP END DO - - !$OMP END PARALLEL - -END_PROVIDER - - - -BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] - BEGIN_DOC - ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active - ! indices are unshifted orbital numbers - END_DOC - implicit none - integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 - double precision, allocatable :: integrals_array(:,:) - real*8 :: mo_two_e_integral - - PROVIDE mo_two_e_integrals_in_map - bielec_PxxQ = 0.d0 - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,ii,j,jj,i3,j3,integrals_array) & - !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PxxQ, & - !$OMP n_act_orb,mo_integrals_map,list_act) - - allocate(integrals_array(mo_num,mo_num)) - - !$OMP DO - do i=1,n_core_inact_orb - ii=list_core_inact(i) - do j=i,n_core_inact_orb - jj=list_core_inact(j) - call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) - do q=1,mo_num - do p=1,mo_num - bielec_PxxQ(p,i,j,q)=integrals_array(p,q) - bielec_PxxQ(p,j,i,q)=integrals_array(q,p) - end do - end do - end do - do j=1,n_act_orb - jj=list_act(j) - j3=j+n_core_inact_orb - call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) - do q=1,mo_num - do p=1,mo_num - bielec_PxxQ(p,i,j3,q)=integrals_array(p,q) - bielec_PxxQ(p,j3,i,q)=integrals_array(q,p) - end do - end do - end do - end do - !$OMP END DO - - - ! (ip|qj) - !$OMP DO - do i=1,n_act_orb - ii=list_act(i) - i3=i+n_core_inact_orb - do j=i,n_act_orb - jj=list_act(j) - j3=j+n_core_inact_orb - call get_mo_two_e_integrals_ij(ii,jj,mo_num,integrals_array,mo_integrals_map) - do q=1,mo_num - do p=1,mo_num - bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q) - bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p) - end do - end do - end do - end do - !$OMP END DO - - deallocate(integrals_array) - !$OMP END PARALLEL - -END_PROVIDER - - -BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] - BEGIN_DOC - ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active - ! index p runs over the whole basis, t,u,v only over the active orbitals - END_DOC - implicit none - integer :: i,j,k,p,t,u,v - double precision, external :: mo_two_e_integral - PROVIDE mo_two_e_integrals_in_map - - !$OMP PARALLEL DO DEFAULT(NONE) & - !$OMP PRIVATE(i,j,k,p,t,u,v) & - !$OMP SHARED(mo_num,n_act_orb,list_act,bielecCI) - do p=1,mo_num - do j=1,n_act_orb - u=list_act(j) - do k=1,n_act_orb - v=list_act(k) - do i=1,n_act_orb - t=list_act(i) - bielecCI(i,k,j,p) = mo_two_e_integral(t,u,v,p) - end do - end do - end do - end do - !$OMP END PARALLEL DO - -END_PROVIDER diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f deleted file mode 100644 index 9968530c..00000000 --- a/src/casscf/bielec_natorb.irp.f +++ /dev/null @@ -1,369 +0,0 @@ - BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] - BEGIN_DOC - ! integral (pq|xx) in the basis of natural MOs - ! indices are unshifted orbital numbers - END_DOC - implicit none - integer :: i,j,k,l,t,u,p,q - double precision, allocatable :: f(:,:,:), d(:,:,:) - - - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(j,k,l,p,d,f) & - !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & - !$OMP bielec_PQxx_no,bielec_PQxx,list_act,natorbsCI) - - allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), & - d(n_act_orb,mo_num,n_core_inact_act_orb)) - - !$OMP DO - do l=1,n_core_inact_act_orb - bielec_PQxx_no(:,:,:,l) = bielec_PQxx(:,:,:,l) - - do k=1,n_core_inact_act_orb - do j=1,mo_num - do p=1,n_act_orb - f(p,j,k)=bielec_PQxx_no(list_act(p),j,k,l) - end do - end do - end do - call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0, & - natorbsCI, size(natorbsCI,1), & - f, n_act_orb, & - 0.d0, & - d, n_act_orb) - do k=1,n_core_inact_act_orb - do j=1,mo_num - do p=1,n_act_orb - bielec_PQxx_no(list_act(p),j,k,l)=d(p,j,k) - end do - end do - - do j=1,mo_num - do p=1,n_act_orb - f(p,j,k)=bielec_PQxx_no(j,list_act(p),k,l) - end do - end do - end do - call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0, & - natorbsCI, n_act_orb, & - f, n_act_orb, & - 0.d0, & - d, n_act_orb) - do k=1,n_core_inact_act_orb - do p=1,n_act_orb - do j=1,mo_num - bielec_PQxx_no(j,list_act(p),k,l)=d(p,j,k) - end do - end do - end do - end do - !$OMP END DO NOWAIT - - deallocate (f,d) - - allocate (f(mo_num,mo_num,n_act_orb),d(mo_num,mo_num,n_act_orb)) - - !$OMP DO - do l=1,n_core_inact_act_orb - - do p=1,n_act_orb - do k=1,mo_num - do j=1,mo_num - f(j,k,p) = bielec_PQxx_no(j,k,n_core_inact_orb+p,l) - end do - end do - end do - call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0, & - f, mo_num*mo_num, & - natorbsCI, n_act_orb, & - 0.d0, & - d, mo_num*mo_num) - do p=1,n_act_orb - do k=1,mo_num - do j=1,mo_num - bielec_PQxx_no(j,k,n_core_inact_orb+p,l)=d(j,k,p) - end do - end do - end do - end do - !$OMP END DO NOWAIT - - !$OMP BARRIER - - !$OMP DO - do l=1,n_core_inact_act_orb - do p=1,n_act_orb - do k=1,mo_num - do j=1,mo_num - f(j,k,p) = bielec_PQxx_no(j,k,l,n_core_inact_orb+p) - end do - end do - end do - call dgemm('N','N',mo_num*mo_num,n_act_orb,n_act_orb,1.d0, & - f, mo_num*mo_num, & - natorbsCI, n_act_orb, & - 0.d0, & - d, mo_num*mo_num) - do p=1,n_act_orb - do k=1,mo_num - do j=1,mo_num - bielec_PQxx_no(j,k,l,n_core_inact_orb+p)=d(j,k,p) - end do - end do - end do - end do - !$OMP END DO - - deallocate (f,d) - !$OMP END PARALLEL - -END_PROVIDER - - - -BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] - BEGIN_DOC - ! integral (px|xq) in the basis of natural MOs - ! indices are unshifted orbital numbers - END_DOC - implicit none - integer :: i,j,k,l,t,u,p,q - double precision, allocatable :: f(:,:,:), d(:,:,:) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(j,k,l,p,d,f) & - !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & - !$OMP bielec_PxxQ_no,bielec_PxxQ,list_act,natorbsCI) - - - allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), & - d(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb)) - - !$OMP DO - do j=1,mo_num - bielec_PxxQ_no(:,:,:,j) = bielec_PxxQ(:,:,:,j) - do l=1,n_core_inact_act_orb - do k=1,n_core_inact_act_orb - do p=1,n_act_orb - f(p,k,l) = bielec_PxxQ_no(list_act(p),k,l,j) - end do - end do - end do - call dgemm('T','N',n_act_orb,n_core_inact_act_orb**2,n_act_orb,1.d0, & - natorbsCI, size(natorbsCI,1), & - f, n_act_orb, & - 0.d0, & - d, n_act_orb) - do l=1,n_core_inact_act_orb - do k=1,n_core_inact_act_orb - do p=1,n_act_orb - bielec_PxxQ_no(list_act(p),k,l,j)=d(p,k,l) - end do - end do - end do - end do - !$OMP END DO NOWAIT - - deallocate (f,d) - - allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), & - d(n_act_orb,mo_num,n_core_inact_act_orb)) - - !$OMP DO - do k=1,mo_num - do l=1,n_core_inact_act_orb - do j=1,mo_num - do p=1,n_act_orb - f(p,j,l) = bielec_PxxQ_no(j,n_core_inact_orb+p,l,k) - end do - end do - end do - call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0, & - natorbsCI, size(natorbsCI,1), & - f, n_act_orb, & - 0.d0, & - d, n_act_orb) - do l=1,n_core_inact_act_orb - do j=1,mo_num - do p=1,n_act_orb - bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p,j,l) - end do - end do - end do - end do - !$OMP END DO NOWAIT - - deallocate(f,d) - - allocate(f(mo_num,n_core_inact_act_orb,n_act_orb), & - d(mo_num,n_core_inact_act_orb,n_act_orb) ) - - !$OMP DO - do k=1,mo_num - do p=1,n_act_orb - do l=1,n_core_inact_act_orb - do j=1,mo_num - f(j,l,p) = bielec_PxxQ_no(j,l,n_core_inact_orb+p,k) - end do - end do - end do - call dgemm('N','N',mo_num*n_core_inact_act_orb,n_act_orb,n_act_orb,1.d0, & - f, mo_num*n_core_inact_act_orb, & - natorbsCI, size(natorbsCI,1), & - 0.d0, & - d, mo_num*n_core_inact_act_orb) - do p=1,n_act_orb - do l=1,n_core_inact_act_orb - do j=1,mo_num - bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(j,l,p) - end do - end do - end do - end do - !$OMP END DO NOWAIT - - !$OMP BARRIER - - !$OMP DO - do l=1,n_core_inact_act_orb - do p=1,n_act_orb - do k=1,n_core_inact_act_orb - do j=1,mo_num - f(j,k,p) = bielec_PxxQ_no(j,k,l,n_core_inact_orb+p) - end do - end do - end do - call dgemm('N','N',mo_num*n_core_inact_act_orb,n_act_orb,n_act_orb,1.d0, & - f, mo_num*n_core_inact_act_orb, & - natorbsCI, size(natorbsCI,1), & - 0.d0, & - d, mo_num*n_core_inact_act_orb) - do p=1,n_act_orb - do k=1,n_core_inact_act_orb - do j=1,mo_num - bielec_PxxQ_no(j,k,l,n_core_inact_orb+p)=d(j,k,p) - end do - end do - end do - end do - !$OMP END DO NOWAIT - deallocate(f,d) - !$OMP END PARALLEL - -END_PROVIDER - - -BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] - BEGIN_DOC - ! integrals (tu|vp) in the basis of natural MOs - ! index p runs over the whole basis, t,u,v only over the active orbitals - END_DOC - implicit none - integer :: i,j,k,l,t,u,p,q - double precision, allocatable :: f(:,:,:), d(:,:,:) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(j,k,l,p,d,f) & - !$OMP SHARED(n_core_inact_act_orb,mo_num,n_act_orb,n_core_inact_orb, & - !$OMP bielecCI_no,bielecCI,list_act,natorbsCI) - - allocate (f(n_act_orb,n_act_orb,mo_num), & - d(n_act_orb,n_act_orb,mo_num)) - - !$OMP DO - do l=1,mo_num - bielecCI_no(:,:,:,l) = bielecCI(:,:,:,l) - do k=1,n_act_orb - do j=1,n_act_orb - do p=1,n_act_orb - f(p,j,k)=bielecCI_no(p,j,k,l) - end do - end do - end do - call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0, & - natorbsCI, size(natorbsCI,1), & - f, n_act_orb, & - 0.d0, & - d, n_act_orb) - do k=1,n_act_orb - do j=1,n_act_orb - do p=1,n_act_orb - bielecCI_no(p,j,k,l)=d(p,j,k) - end do - end do - - do j=1,n_act_orb - do p=1,n_act_orb - f(p,j,k)=bielecCI_no(j,p,k,l) - end do - end do - end do - call dgemm('T','N',n_act_orb,n_act_orb*n_act_orb,n_act_orb,1.d0, & - natorbsCI, n_act_orb, & - f, n_act_orb, & - 0.d0, & - d, n_act_orb) - do k=1,n_act_orb - do p=1,n_act_orb - do j=1,n_act_orb - bielecCI_no(j,p,k,l)=d(p,j,k) - end do - end do - end do - - do p=1,n_act_orb - do k=1,n_act_orb - do j=1,n_act_orb - f(j,k,p)=bielecCI_no(j,k,p,l) - end do - end do - end do - call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0, & - f, n_act_orb*n_act_orb, & - natorbsCI, n_act_orb, & - 0.d0, & - d, n_act_orb*n_act_orb) - - do p=1,n_act_orb - do k=1,n_act_orb - do j=1,n_act_orb - bielecCI_no(j,k,p,l)=d(j,k,p) - end do - end do - end do - end do - !$OMP END DO - - !$OMP DO - do l=1,n_act_orb - do p=1,n_act_orb - do k=1,n_act_orb - do j=1,n_act_orb - f(j,k,p)=bielecCI_no(j,k,l,list_act(p)) - end do - end do - end do - call dgemm('N','N',n_act_orb*n_act_orb,n_act_orb,n_act_orb,1.d0, & - f, n_act_orb*n_act_orb, & - natorbsCI, n_act_orb, & - 0.d0, & - d, n_act_orb*n_act_orb) - - do p=1,n_act_orb - do k=1,n_act_orb - do j=1,n_act_orb - bielecCI_no(j,k,l,list_act(p))=d(j,k,p) - end do - end do - end do - end do - !$OMP END DO - - deallocate(d,f) - !$OMP END PARALLEL - - -END_PROVIDER - diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f deleted file mode 100644 index 950cfd55..00000000 --- a/src/casscf/casscf.irp.f +++ /dev/null @@ -1,58 +0,0 @@ -program casscf - implicit none - BEGIN_DOC -! TODO : Put the documentation of the program here - END_DOC - call reorder_orbitals_for_casscf - no_vvvv_integrals = .True. - touch no_vvvv_integrals - pt2_max = 0.02 - SOFT_TOUCH no_vvvv_integrals pt2_max - call run_stochastic_cipsi - call run -end - -subroutine run - implicit none - double precision :: energy_old, energy - logical :: converged,state_following_casscf_save - integer :: iteration - converged = .False. - - energy = 0.d0 - mo_label = "MCSCF" - iteration = 1 - state_following_casscf_save = state_following_casscf - state_following_casscf = .True. - touch state_following_casscf - do while (.not.converged) - call run_stochastic_cipsi - energy_old = energy - energy = eone+etwo+ecore - - call write_time(6) - call write_int(6,iteration,'CAS-SCF iteration') - call write_double(6,energy,'CAS-SCF energy') - call write_double(6,energy_improvement, 'Predicted energy improvement') - - converged = dabs(energy_improvement) < thresh_scf - pt2_max = dabs(energy_improvement / pt2_relative_error) - - mo_coef = NewOrbs - mo_occ = occnum - call save_mos - iteration += 1 - N_det = max(N_det/2 ,N_states) - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - read_wf = .True. - call clear_mo_map - SOFT_TOUCH mo_coef N_det pt2_max psi_det psi_coef - if(iteration .gt. 3)then - state_following_casscf = state_following_casscf_save - touch state_following_casscf - endif - - enddo - -end diff --git a/src/casscf/class.irp.f b/src/casscf/class.irp.f deleted file mode 100644 index 7360a661..00000000 --- a/src/casscf/class.irp.f +++ /dev/null @@ -1,12 +0,0 @@ - BEGIN_PROVIDER [ logical, do_only_1h1p ] -&BEGIN_PROVIDER [ logical, do_only_cas ] -&BEGIN_PROVIDER [ logical, do_ddci ] - implicit none - BEGIN_DOC - ! In the CAS case, all those are always false except do_only_cas - END_DOC - do_only_cas = .True. - do_only_1h1p = .False. - do_ddci = .False. -END_PROVIDER - diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f deleted file mode 100644 index d181d732..00000000 --- a/src/casscf/densities.irp.f +++ /dev/null @@ -1,63 +0,0 @@ -use bitmasks - -BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ] - implicit none - BEGIN_DOC - ! the first-order density matrix in the basis of the starting MOs. - ! matrix is state averaged. - END_DOC - integer :: t,u - - do u=1,n_act_orb - do t=1,n_act_orb - D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + & - one_e_dm_mo_beta_average ( list_act(t), list_act(u) ) - enddo - enddo - -END_PROVIDER - -BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] - BEGIN_DOC - ! The second-order density matrix in the basis of the starting MOs ONLY IN THE RANGE OF ACTIVE MOS - ! The values are state averaged - ! - ! We use the spin-free generators of mono-excitations - ! E_pq destroys q and creates p - ! D_pq = <0|E_pq|0> = D_qp - ! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0> - ! - ! P0tuvx(p,q,r,s) = chemist notation : 1/2 <0|E_pq E_rs - delta_qr E_ps|0> - END_DOC - implicit none - integer :: t,u,v,x - integer :: tt,uu,vv,xx - integer :: mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart - integer :: ierr - real*8 :: phase1,phase11,phase12,phase2,phase21,phase22 - integer :: nu1,nu2,nu11,nu12,nu21,nu22 - integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22 - real*8 :: cI_mu(N_states),term - integer(bit_kind), dimension(N_int,2) :: det_mu, det_mu_ex - integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12 - integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22 - - if (bavard) then - write(6,*) ' providing the 2 body RDM on the active part' - endif - - P0tuvx= 0.d0 - do istate=1,N_states - do x = 1, n_act_orb - do v = 1, n_act_orb - do u = 1, n_act_orb - do t = 1, n_act_orb - ! 1 1 2 2 1 2 1 2 - P0tuvx(t,u,v,x) = state_av_act_2_rdm_spin_trace_mo(t,v,u,x) - enddo - enddo - enddo - enddo - enddo - -END_PROVIDER diff --git a/src/casscf/det_manip.irp.f b/src/casscf/det_manip.irp.f deleted file mode 100644 index d8c309a4..00000000 --- a/src/casscf/det_manip.irp.f +++ /dev/null @@ -1,125 +0,0 @@ -use bitmasks - -subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, & - ispin,phase,ierr) - BEGIN_DOC - ! we create the mono-excitation, and determine, if possible, - ! the phase and the number in the list of determinants - END_DOC - implicit none - integer(bit_kind) :: key1(N_int,2),key2(N_int,2) - integer(bit_kind), allocatable :: keytmp(:,:) - integer :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin - real*8 :: phase - logical :: found - allocate(keytmp(N_int,2)) - - nu=-1 - phase=1.D0 - ierr=0 - call det_copy(key1,key2,N_int) - ! write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin - ! call print_det(key2,N_int) - call do_single_excitation(key2,ihole,ipart,ispin,ierr) - ! write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin - ! call print_det(key2,N_int) - ! write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr - if (ierr.eq.1) then - ! excitation is possible - ! get the phase - call get_single_excitation(key1,key2,exc,phase,N_int) - ! get the number in the list - found=.false. - nu=0 - - !TODO BOTTLENECK - do while (.not.found) - nu+=1 - if (nu.gt.N_det) then - ! the determinant is possible, but not in the list - found=.true. - nu=-1 - else - call det_extract(keytmp,nu,N_int) - integer :: i,ii - found=.true. - do ii=1,2 - do i=1,N_int - if (keytmp(i,ii).ne.key2(i,ii)) then - found=.false. - end if - end do - end do - end if - end do - end if - ! - ! we found the new string, the phase, and possibly the number in the list - ! -end subroutine do_signed_mono_excitation - -subroutine det_extract(key,nu,Nint) - BEGIN_DOC - ! extract a determinant from the list of determinants - END_DOC - implicit none - integer :: ispin,i,nu,Nint - integer(bit_kind) :: key(Nint,2) - do ispin=1,2 - do i=1,Nint - key(i,ispin)=psi_det(i,ispin,nu) - end do - end do -end subroutine det_extract - -subroutine det_copy(key1,key2,Nint) - use bitmasks ! you need to include the bitmasks_module.f90 features - BEGIN_DOC - ! copy a determinant from key1 to key2 - END_DOC - implicit none - integer :: ispin,i,Nint - integer(bit_kind) :: key1(Nint,2),key2(Nint,2) - do ispin=1,2 - do i=1,Nint - key2(i,ispin)=key1(i,ispin) - end do - end do -end subroutine det_copy - -subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2 & - ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr) - BEGIN_DOC - ! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q) - ! we may create two determinants as result - ! - END_DOC - implicit none - integer(bit_kind) :: key_in(N_int,2),key_out1(N_int,2) - integer(bit_kind) :: key_out2(N_int,2) - integer :: ihole,ipart,ierr,jerr,nu1,nu2 - integer :: ispin - real*8 :: phase1,phase2 - - ! write(6,*) ' applying E_',ipart,ihole,' on determinant ' - ! call print_det(key_in,N_int) - - ! spin alpha - ispin=1 - call do_signed_mono_excitation(key_in,key_out1,nu1,ihole & - ,ipart,ispin,phase1,ierr) - ! if (ierr.eq.1) then - ! write(6,*) ' 1 result is ',nu1,phase1 - ! call print_det(key_out1,N_int) - ! end if - ! spin beta - ispin=2 - call do_signed_mono_excitation(key_in,key_out2,nu2,ihole & - ,ipart,ispin,phase2,jerr) - ! if (jerr.eq.1) then - ! write(6,*) ' 2 result is ',nu2,phase2 - ! call print_det(key_out2,N_int) - ! end if - -end subroutine do_spinfree_mono_excitation - diff --git a/src/casscf/driver_optorb.irp.f b/src/casscf/driver_optorb.irp.f deleted file mode 100644 index 2e3e02dc..00000000 --- a/src/casscf/driver_optorb.irp.f +++ /dev/null @@ -1,3 +0,0 @@ -subroutine driver_optorb - implicit none -end diff --git a/src/casscf/get_energy.irp.f b/src/casscf/get_energy.irp.f deleted file mode 100644 index cfb26b59..00000000 --- a/src/casscf/get_energy.irp.f +++ /dev/null @@ -1,51 +0,0 @@ -program print_2rdm - implicit none - BEGIN_DOC - ! get the active part of the bielectronic energy on a given wave function. - ! - ! useful to test the active part of the spin trace 2 rdms - END_DOC -!no_vvvv_integrals = .True. - read_wf = .True. -!touch read_wf no_vvvv_integrals -!call routine -!call routine_bis - call print_grad -end - -subroutine print_grad - implicit none - integer :: i - do i = 1, nMonoEx - if(dabs(gradvec2(i)).gt.1.d-5)then - print*,'' - print*,i,gradvec2(i),excit(:,i) - endif - enddo -end - -subroutine routine - integer :: i,j,k,l - integer :: ii,jj,kk,ll - double precision :: accu(4),twodm,thr,act_twodm2,integral,get_two_e_integral - thr = 1.d-10 - - - accu = 0.d0 - do ll = 1, n_act_orb - l = list_act(ll) - do kk = 1, n_act_orb - k = list_act(kk) - do jj = 1, n_act_orb - j = list_act(jj) - do ii = 1, n_act_orb - i = list_act(ii) - integral = get_two_e_integral(i,j,k,l,mo_integrals_map) - accu(1) += state_av_act_2_rdm_spin_trace_mo(ii,jj,kk,ll) * integral - enddo - enddo - enddo - enddo - print*,'accu = ',accu(1) - -end diff --git a/src/casscf/grad_old.irp.f b/src/casscf/grad_old.irp.f deleted file mode 100644 index d60a60c8..00000000 --- a/src/casscf/grad_old.irp.f +++ /dev/null @@ -1,74 +0,0 @@ - -BEGIN_PROVIDER [real*8, gradvec_old, (nMonoEx)] - BEGIN_DOC - ! calculate the orbital gradient by hand, i.e. for - ! each determinant I we determine the string E_pq |I> (alpha and beta - ! separately) and generate - ! sum_I c_I is then the pq component of the orbital - ! gradient - ! E_pq = a^+_pa_q + a^+_Pa_Q - END_DOC - implicit none - integer :: ii,tt,aa,indx,ihole,ipart,istate - real*8 :: res - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - call calc_grad_elem(ihole,ipart,res) - gradvec_old(indx)=res - end do - - real*8 :: norm_grad - norm_grad=0.d0 - do indx=1,nMonoEx - norm_grad+=gradvec_old(indx)*gradvec_old(indx) - end do - norm_grad=sqrt(norm_grad) - if (bavard) then - write(6,*) - write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad - write(6,*) - endif - - -END_PROVIDER - -subroutine calc_grad_elem(ihole,ipart,res) - BEGIN_DOC - ! eq 18 of Siegbahn et al, Physica Scripta 1980 - ! we calculate 2 , q=hole, p=particle - END_DOC - implicit none - integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) - real*8 :: i_H_psi_array(N_states),phase - allocate(det_mu(N_int,2)) - allocate(det_mu_ex(N_int,2)) - - res=0.D0 - - do mu=1,n_det - ! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 - ! do the monoexcitation on it - call det_copy(det_mu,det_mu_ex,N_int) - call do_signed_mono_excitation(det_mu,det_mu_ex,nu & - ,ihole,ipart,ispin,phase,ierr) - if (ierr.eq.1) then - call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase - end do - end if - end do - end do - - ! state-averaged gradient - res*=2.D0/dble(N_states) - -end subroutine calc_grad_elem - diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f deleted file mode 100644 index e717e822..00000000 --- a/src/casscf/gradient.irp.f +++ /dev/null @@ -1,171 +0,0 @@ -use bitmasks - -BEGIN_PROVIDER [ integer, nMonoEx ] - BEGIN_DOC - ! Number of single excitations - END_DOC - implicit none - nMonoEx=n_core_inact_orb*n_act_orb+n_core_inact_orb*n_virt_orb+n_act_orb*n_virt_orb -END_PROVIDER - - BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] -&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] - BEGIN_DOC - ! a list of the orbitals involved in the excitation - END_DOC - - implicit none - integer :: i,t,a,ii,tt,aa,indx - indx=0 - do ii=1,n_core_inact_orb - i=list_core_inact(ii) - do tt=1,n_act_orb - t=list_act(tt) - indx+=1 - excit(1,indx)=i - excit(2,indx)=t - excit_class(indx)='c-a' - end do - end do - - do ii=1,n_core_inact_orb - i=list_core_inact(ii) - do aa=1,n_virt_orb - a=list_virt(aa) - indx+=1 - excit(1,indx)=i - excit(2,indx)=a - excit_class(indx)='c-v' - end do - end do - - do tt=1,n_act_orb - t=list_act(tt) - do aa=1,n_virt_orb - a=list_virt(aa) - indx+=1 - excit(1,indx)=t - excit(2,indx)=a - excit_class(indx)='a-v' - end do - end do - - if (bavard) then - write(6,*) ' Filled the table of the Monoexcitations ' - do indx=1,nMonoEx - write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' & - ,excit(2,indx),' ',excit_class(indx) - end do - end if - -END_PROVIDER - -BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] - BEGIN_DOC - ! calculate the orbital gradient from density - ! matrices and integrals; Siegbahn et al, Phys Scr 1980 - ! eqs 14 a,b,c - END_DOC - implicit none - integer :: i,t,a,indx - real*8 :: gradvec_it,gradvec_ia,gradvec_ta - real*8 :: norm_grad - - indx=0 - do i=1,n_core_inact_orb - do t=1,n_act_orb - indx+=1 - gradvec2(indx)=gradvec_it(i,t) - end do - end do - - do i=1,n_core_inact_orb - do a=1,n_virt_orb - indx+=1 - gradvec2(indx)=gradvec_ia(i,a) - end do - end do - - do t=1,n_act_orb - do a=1,n_virt_orb - indx+=1 - gradvec2(indx)=gradvec_ta(t,a) - end do - end do - - norm_grad=0.d0 - do indx=1,nMonoEx - norm_grad+=gradvec2(indx)*gradvec2(indx) - end do - norm_grad=sqrt(norm_grad) - write(6,*) - write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad - write(6,*) - -END_PROVIDER - -real*8 function gradvec_it(i,t) - BEGIN_DOC - ! the orbital gradient core/inactive -> active - ! we assume natural orbitals - END_DOC - implicit none - integer :: i,t - - integer :: ii,tt,v,vv,x,y - integer :: x3,y3 - - ii=list_core_inact(i) - tt=list_act(t) - gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii)) - gradvec_it-=occnum(tt)*Fipq(ii,tt) - do v=1,n_act_orb - vv=list_act(v) - do x=1,n_act_orb - x3=x+n_core_inact_orb - do y=1,n_act_orb - y3=y+n_core_inact_orb - gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) - end do - end do - end do - gradvec_it*=2.D0 -end function gradvec_it - -real*8 function gradvec_ia(i,a) - BEGIN_DOC - ! the orbital gradient core/inactive -> virtual - END_DOC - implicit none - integer :: i,a,ii,aa - - ii=list_core_inact(i) - aa=list_virt(a) - gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) - gradvec_ia*=2.D0 - -end function gradvec_ia - -real*8 function gradvec_ta(t,a) - BEGIN_DOC - ! the orbital gradient active -> virtual - ! we assume natural orbitals - END_DOC - implicit none - integer :: t,a,tt,aa,v,vv,x,y - - tt=list_act(t) - aa=list_virt(a) - gradvec_ta=0.D0 - gradvec_ta+=occnum(tt)*Fipq(aa,tt) - do v=1,n_act_orb - do x=1,n_act_orb - do y=1,n_act_orb - gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa) - end do - end do - end do - gradvec_ta*=2.D0 - -end function gradvec_ta - diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f deleted file mode 100644 index 52be1b76..00000000 --- a/src/casscf/hessian.irp.f +++ /dev/null @@ -1,656 +0,0 @@ -use bitmasks - -BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)] - BEGIN_DOC - ! calculate the orbital hessian 2 - ! + + by hand, - ! determinant per determinant, as for the gradient - ! - ! we assume that we have natural active orbitals - END_DOC - implicit none - integer :: indx,ihole,ipart - integer :: jndx,jhole,jpart - character*3 :: iexc,jexc - real*8 :: res - - if (bavard) then - write(6,*) ' providing Hessian matrix hessmat ' - write(6,*) ' nMonoEx = ',nMonoEx - endif - - do indx=1,nMonoEx - do jndx=1,nMonoEx - hessmat(indx,jndx)=0.D0 - end do - end do - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - iexc=excit_class(indx) - do jndx=indx,nMonoEx - jhole=excit(1,jndx) - jpart=excit(2,jndx) - jexc=excit_class(jndx) - call calc_hess_elem(ihole,ipart,jhole,jpart,res) - hessmat(indx,jndx)=res - hessmat(jndx,indx)=res - end do - end do - -END_PROVIDER - -subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) - BEGIN_DOC - ! eq 19 of Siegbahn et al, Physica Scripta 1980 - ! we calculate 2 - ! + + - ! average over all states is performed. - ! no transition between states. - END_DOC - implicit none - integer :: ihole,ipart,ispin,mu,istate - integer :: jhole,jpart,jspin - integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:) - integer(bit_kind), allocatable :: det_nu(:,:) - integer(bit_kind), allocatable :: det_mu_pq(:,:) - integer(bit_kind), allocatable :: det_mu_rs(:,:) - integer(bit_kind), allocatable :: det_nu_rs(:,:) - integer(bit_kind), allocatable :: det_mu_pqrs(:,:) - integer(bit_kind), allocatable :: det_mu_rspq(:,:) - real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 - real*8 :: i_H_j_element - allocate(det_mu(N_int,2)) - allocate(det_nu(N_int,2)) - allocate(det_mu_pq(N_int,2)) - allocate(det_mu_rs(N_int,2)) - allocate(det_nu_rs(N_int,2)) - allocate(det_mu_pqrs(N_int,2)) - allocate(det_mu_rspq(N_int,2)) - integer :: mu_pq_possible - integer :: mu_rs_possible - integer :: nu_rs_possible - integer :: mu_pqrs_possible - integer :: mu_rspq_possible - - res=0.D0 - - ! the terms <0|E E H |0> - do mu=1,n_det - ! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 - ! do the monoexcitation pq on it - call det_copy(det_mu,det_mu_pq,N_int) - call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & - ,ihole,ipart,ispin,phase,mu_pq_possible) - if (mu_pq_possible.eq.1) then - ! possible, but not necessarily in the list - ! do the second excitation - do jspin=1,2 - call det_copy(det_mu_pq,det_mu_pqrs,N_int) - call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& - ,jhole,jpart,jspin,phase2,mu_pqrs_possible) - ! excitation possible - if (mu_pqrs_possible.eq.1) then - call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 - end do - end if - ! try the de-excitation with opposite sign - call det_copy(det_mu_pq,det_mu_pqrs,N_int) - call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& - ,jpart,jhole,jspin,phase2,mu_pqrs_possible) - phase2=-phase2 - ! excitation possible - if (mu_pqrs_possible.eq.1) then - call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 - end do - end if - end do - end if - ! exchange the notion of pq and rs - ! do the monoexcitation rs on the initial determinant - call det_copy(det_mu,det_mu_rs,N_int) - call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & - ,jhole,jpart,ispin,phase2,mu_rs_possible) - if (mu_rs_possible.eq.1) then - ! do the second excitation - do jspin=1,2 - call det_copy(det_mu_rs,det_mu_rspq,N_int) - call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& - ,ihole,ipart,jspin,phase3,mu_rspq_possible) - ! excitation possible (of course, the result is outside the CAS) - if (mu_rspq_possible.eq.1) then - call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 - end do - end if - ! we may try the de-excitation, with opposite sign - call det_copy(det_mu_rs,det_mu_rspq,N_int) - call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& - ,ipart,ihole,jspin,phase3,mu_rspq_possible) - phase3=-phase3 - ! excitation possible (of course, the result is outside the CAS) - if (mu_rspq_possible.eq.1) then - call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 - end do - end if - end do - end if - ! - ! the operator E H E, we have to do a double loop over the determinants - ! we still have the determinant mu_pq and the phase in memory - if (mu_pq_possible.eq.1) then - do nu=1,N_det - call det_extract(det_nu,nu,N_int) - do jspin=1,2 - call det_copy(det_nu,det_nu_rs,N_int) - call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & - ,jhole,jpart,jspin,phase2,nu_rs_possible) - ! excitation possible ? - if (nu_rs_possible.eq.1) then - call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) - do istate=1,N_states - res+=2.D0*i_H_j_element*psi_coef(mu,istate) & - *psi_coef(nu,istate)*phase*phase2 - end do - end if - end do - end do - end if - end do - end do - - ! state-averaged Hessian - res*=1.D0/dble(N_states) - -end subroutine calc_hess_elem - -BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] - BEGIN_DOC - ! explicit hessian matrix from density matrices and integrals - ! of course, this will be used for a direct Davidson procedure later - ! we will not store the matrix in real life - ! formulas are broken down as functions for the 6 classes of matrix elements - ! - END_DOC - implicit none - integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart,indx_shift - - real*8 :: hessmat_itju - real*8 :: hessmat_itja - real*8 :: hessmat_itua - real*8 :: hessmat_iajb - real*8 :: hessmat_iatb - real*8 :: hessmat_taub - - if (bavard) then - write(6,*) ' providing Hessian matrix hessmat2 ' - write(6,*) ' nMonoEx = ',nMonoEx - endif - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(hessmat2,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & - !$OMP PRIVATE(i,indx,jndx,j,ustart,t,u,a,bstart,indx_shift) - - !$OMP DO - do i=1,n_core_inact_orb - do t=1,n_act_orb - indx = t + (i-1)*n_act_orb - jndx=indx - do j=i,n_core_inact_orb - if (i.eq.j) then - ustart=t - else - ustart=1 - end if - do u=ustart,n_act_orb - hessmat2(jndx,indx)=hessmat_itju(i,t,j,u) - jndx+=1 - end do - end do - do j=1,n_core_inact_orb - do a=1,n_virt_orb - hessmat2(jndx,indx)=hessmat_itja(i,t,j,a) - jndx+=1 - end do - end do - do u=1,n_act_orb - do a=1,n_virt_orb - hessmat2(jndx,indx)=hessmat_itua(i,t,u,a) - jndx+=1 - end do - end do - end do - end do - !$OMP END DO NOWAIT - - indx_shift = n_core_inact_orb*n_act_orb - !$OMP DO - do a=1,n_virt_orb - do i=1,n_core_inact_orb - indx = a + (i-1)*n_virt_orb + indx_shift - jndx=indx - do j=i,n_core_inact_orb - if (i.eq.j) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - hessmat2(jndx,indx)=hessmat_iajb(i,a,j,b) - jndx+=1 - end do - end do - do t=1,n_act_orb - do b=1,n_virt_orb - hessmat2(jndx,indx)=hessmat_iatb(i,a,t,b) - jndx+=1 - end do - end do - end do - end do - !$OMP END DO NOWAIT - - indx_shift += n_core_inact_orb*n_virt_orb - !$OMP DO - do a=1,n_virt_orb - do t=1,n_act_orb - indx = a + (t-1)*n_virt_orb + indx_shift - jndx=indx - do u=t,n_act_orb - if (t.eq.u) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - hessmat2(jndx,indx)=hessmat_taub(t,a,u,b) - jndx+=1 - end do - end do - end do - end do - !$OMP END DO - - !$OMP END PARALLEL - - do jndx=1,nMonoEx - do indx=1,jndx-1 - hessmat2(indx,jndx) = hessmat2(jndx,indx) - enddo - enddo - - -END_PROVIDER - -real*8 function hessmat_itju(i,t,j,u) - BEGIN_DOC - ! the orbital hessian for core/inactive -> active, core/inactive -> active - ! i, t, j, u are list indices, the corresponding orbitals are ii,tt,jj,uu - ! - ! we assume natural orbitals - END_DOC - implicit none - integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj - real*8 :: term,t2 - - ii=list_core_inact(i) - tt=list_act(t) - if (i.eq.j) then - if (t.eq.u) then - ! diagonal element - term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) & - -2.D0*(Fipq(ii,ii)+Fapq(ii,ii)) - term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i)) - term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt) & - -bielec_pqxx_no(tt,tt,i,i)) - term-=occnum(tt)*Fipq(tt,tt) - do v=1,n_act_orb - vv=list_act(v) - do x=1,n_act_orb - xx=list_act(x) - term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(vv,xx,i,i) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & - bielec_pxxq_no(vv,i,i,xx)) - do y=1,n_act_orb - term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) - end do - end do - end do - else - ! it/iu, t != u - uu=list_act(u) - term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu)) - term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & - -bielec_PQxx_no(tt,uu,i,j)) - term-=occnum(tt)*Fipq(uu,tt) - term-=(occnum(tt)+occnum(uu)) & - *(3.D0*bielec_PxxQ_no(tt,i,i,uu)-bielec_PQxx_no(uu,tt,i,i)) - do v=1,n_act_orb - vv=list_act(v) - ! term-=D0tu(u,v)*Fipq(tt,vv) ! published, but inverting t and u seems more correct - do x=1,n_act_orb - xx=list_act(x) - term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,i) & - +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & - *bielec_pxxq_no(vv,i,i,xx)) - do y=1,n_act_orb - term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(u,v,y,xx) - end do - end do - end do - end if - else - ! it/ju - jj=list_core_inact(j) - uu=list_act(u) - if (t.eq.u) then - term=occnum(tt)*Fipq(ii,jj) - term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) - else - term=0.D0 - end if - term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & - -bielec_PQxx_no(tt,uu,i,j)) - term-=(occnum(tt)+occnum(uu))* & - (4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & - -bielec_PQxx_no(uu,tt,i,j)) - do v=1,n_act_orb - vv=list_act(v) - do x=1,n_act_orb - xx=list_act(x) - term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,j) & - +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & - *bielec_pxxq_no(vv,i,j,xx)) - end do - end do - end if - - term*=2.D0 - hessmat_itju=term - -end function hessmat_itju - -real*8 function hessmat_itja(i,t,j,a) - BEGIN_DOC - ! the orbital hessian for core/inactive -> active, core/inactive -> virtual - END_DOC - implicit none - integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y - real*8 :: term - - ! it/ja - ii=list_core_inact(i) - tt=list_act(t) - jj=list_core_inact(j) - aa=list_virt(a) - term=2.D0*(4.D0*bielec_pxxq_no(aa,j,i,tt) & - -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt)) - term-=occnum(tt)*(4.D0*bielec_pxxq_no(aa,j,i,tt) & - -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt)) - if (i.eq.j) then - term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt)) - term-=0.5D0*occnum(tt)*Fipq(aa,tt) - do v=1,n_act_orb - do x=1,n_act_orb - do y=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa) - end do - end do - end do - end if - term*=2.D0 - hessmat_itja=term - -end function hessmat_itja - -real*8 function hessmat_itua(i,t,u,a) - BEGIN_DOC - ! the orbital hessian for core/inactive -> active, active -> virtual - END_DOC - implicit none - integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 - real*8 :: term - - ii=list_core_inact(i) - tt=list_act(t) - t3=t+n_core_inact_orb - uu=list_act(u) - u3=u+n_core_inact_orb - aa=list_virt(a) - if (t.eq.u) then - term=-occnum(tt)*Fipq(aa,ii) - else - term=0.D0 - end if - term-=occnum(uu)*(bielec_pqxx_no(aa,ii,t3,u3)-4.D0*bielec_pqxx_no(aa,uu,t3,i)& - +bielec_pxxq_no(aa,t3,u3,ii)) - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - do x=1,n_act_orb - integer :: x3 - xx=list_act(x) - x3=x+n_core_inact_orb - term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,ii,v3,x3) & - +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u)) & - *bielec_pqxx_no(aa,xx,v3,i)) - end do - end do - if (t.eq.u) then - term+=Fipq(aa,ii)+Fapq(aa,ii) - end if - term*=2.D0 - hessmat_itua=term - -end function hessmat_itua - -real*8 function hessmat_iajb(i,a,j,b) - BEGIN_DOC - ! the orbital hessian for core/inactive -> virtual, core/inactive -> virtual - END_DOC - implicit none - integer :: i,a,j,b,ii,aa,jj,bb - real*8 :: term - - ii=list_core_inact(i) - aa=list_virt(a) - if (i.eq.j) then - if (a.eq.b) then - ! ia/ia - term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii)) - term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,aa)-bielec_pqxx_no(aa,aa,i,i)) - else - bb=list_virt(b) - ! ia/ib - term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb)) - term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,bb)-bielec_pqxx_no(aa,bb,i,i)) - end if - else - ! ia/jb - jj=list_core_inact(j) - bb=list_virt(b) - term=2.D0*(4.D0*bielec_pxxq_no(aa,i,j,bb)-bielec_pqxx_no(aa,bb,i,j) & - -bielec_pxxq_no(aa,j,i,bb)) - if (a.eq.b) then - term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) - end if - end if - term*=2.D0 - hessmat_iajb=term - -end function hessmat_iajb - -real*8 function hessmat_iatb(i,a,t,b) - BEGIN_DOC - ! the orbital hessian for core/inactive -> virtual, active -> virtual - END_DOC - implicit none - integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 - real*8 :: term - - ii=list_core_inact(i) - aa=list_virt(a) - tt=list_act(t) - bb=list_virt(b) - t3=t+n_core_inact_orb - term=occnum(tt)*(4.D0*bielec_pxxq_no(aa,i,t3,bb)-bielec_pxxq_no(aa,t3,i,bb)& - -bielec_pqxx_no(aa,bb,i,t3)) - if (a.eq.b) then - term-=Fipq(tt,ii)+Fapq(tt,ii) - term-=0.5D0*occnum(tt)*Fipq(tt,ii) - do v=1,n_act_orb - do x=1,n_act_orb - do y=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,ii) - end do - end do - end do - end if - term*=2.D0 - hessmat_iatb=term - -end function hessmat_iatb - -real*8 function hessmat_taub(t,a,u,b) - BEGIN_DOC - ! the orbital hessian for act->virt,act->virt - END_DOC - implicit none - integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y - integer :: v3,x3 - real*8 :: term,t1,t2,t3 - - tt=list_act(t) - aa=list_virt(a) - if (t == u) then - if (a == b) then - ! ta/ta - t1=occnum(tt)*Fipq(aa,aa) - t2=0.D0 - t3=0.D0 - t1-=occnum(tt)*Fipq(tt,tt) - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & - bielec_pxxq_no(aa,x3,v3,aa)) - do y=1,n_act_orb - t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) - end do - end do - end do - term=t1+t2+t3 - else - bb=list_virt(b) - ! ta/tb b/=a - term=occnum(tt)*Fipq(aa,bb) - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & - *bielec_pxxq_no(aa,x3,v3,bb)) - end do - end do - end if - else - ! ta/ub t/=u - uu=list_act(u) - bb=list_virt(b) - term=0.D0 - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & - +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & - *bielec_pxxq_no(aa,x3,v3,bb)) - end do - end do - if (a.eq.b) then - term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) - do v=1,n_act_orb - do y=1,n_act_orb - do x=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) - term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt) - end do - end do - end do - end if - - end if - - term*=2.D0 - hessmat_taub=term - -end function hessmat_taub - -BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] - BEGIN_DOC - ! the diagonal of the Hessian, needed for the Davidson procedure - END_DOC - implicit none - integer :: i,t,a,indx,indx_shift - real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(hessdiag,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & - !$OMP PRIVATE(i,indx,t,a,indx_shift) - - !$OMP DO - do i=1,n_core_inact_orb - do t=1,n_act_orb - indx = t + (i-1)*n_act_orb - hessdiag(indx)=hessmat_itju(i,t,i,t) - end do - end do - !$OMP END DO NOWAIT - - indx_shift = n_core_inact_orb*n_act_orb - !$OMP DO - do a=1,n_virt_orb - do i=1,n_core_inact_orb - indx = a + (i-1)*n_virt_orb + indx_shift - hessdiag(indx)=hessmat_iajb(i,a,i,a) - end do - end do - !$OMP END DO NOWAIT - - indx_shift += n_core_inact_orb*n_virt_orb - !$OMP DO - do a=1,n_virt_orb - do t=1,n_act_orb - indx = a + (t-1)*n_virt_orb + indx_shift - hessdiag(indx)=hessmat_taub(t,a,t,a) - end do - end do - !$OMP END DO - !$OMP END PARALLEL - -END_PROVIDER diff --git a/src/casscf/mcscf_fock.irp.f b/src/casscf/mcscf_fock.irp.f deleted file mode 100644 index e4568405..00000000 --- a/src/casscf/mcscf_fock.irp.f +++ /dev/null @@ -1,80 +0,0 @@ -BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] - BEGIN_DOC - ! the inactive Fock matrix, in molecular orbitals - END_DOC - implicit none - integer :: p,q,k,kk,t,tt,u,uu - - do q=1,mo_num - do p=1,mo_num - Fipq(p,q)=one_ints_no(p,q) - end do - end do - - ! the inactive Fock matrix - do k=1,n_core_inact_orb - kk=list_core_inact(k) - do q=1,mo_num - do p=1,mo_num - Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q) - end do - end do - end do - - if (bavard) then - integer :: i - write(6,*) - write(6,*) ' the diagonal of the inactive effective Fock matrix ' - write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) - write(6,*) - end if - - -END_PROVIDER - - -BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] - BEGIN_DOC - ! the active active Fock matrix, in molecular orbitals - ! we create them in MOs, quite expensive - ! - ! for an implementation in AOs we need first the natural orbitals - ! for forming an active density matrix in AOs - ! - END_DOC - implicit none - integer :: p,q,k,kk,t,tt,u,uu - - Fapq = 0.d0 - - ! the active Fock matrix, D0tu is diagonal - do t=1,n_act_orb - tt=list_act(t) - do q=1,mo_num - do p=1,mo_num - Fapq(p,q)+=occnum(tt) & - *(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q)) - end do - end do - end do - - if (bavard) then - integer :: i - write(6,*) - write(6,*) ' the effective Fock matrix over MOs' - write(6,*) - - write(6,*) - write(6,*) ' the diagonal of the inactive effective Fock matrix ' - write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) - write(6,*) - write(6,*) - write(6,*) ' the diagonal of the active Fock matrix ' - write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) - write(6,*) - end if - - -END_PROVIDER - - diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f deleted file mode 100644 index 9ce90304..00000000 --- a/src/casscf/natorb.irp.f +++ /dev/null @@ -1,231 +0,0 @@ - BEGIN_PROVIDER [real*8, occnum, (mo_num)] - implicit none - BEGIN_DOC - ! MO occupation numbers - END_DOC - - integer :: i - occnum=0.D0 - do i=1,n_core_inact_orb - occnum(list_core_inact(i))=2.D0 - end do - - do i=1,n_act_orb - occnum(list_act(i))=occ_act(i) - end do - - if (bavard) then - write(6,*) ' occupation numbers ' - do i=1,mo_num - write(6,*) i,occnum(i) - end do - endif - -END_PROVIDER - - - BEGIN_PROVIDER [ real*8, natorbsCI, (n_act_orb,n_act_orb) ] -&BEGIN_PROVIDER [ real*8, occ_act, (n_act_orb) ] - implicit none - BEGIN_DOC - ! Natural orbitals of CI - END_DOC - integer :: i, j - double precision :: Vt(n_act_orb,n_act_orb) - -! call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb) - call svd(D0tu, size(D0tu,1), natorbsCI,size(natorbsCI,1), occ_act, Vt, size(Vt,1),n_act_orb,n_act_orb) - - if (bavard) then - write(6,*) ' found occupation numbers as ' - do i=1,n_act_orb - write(6,*) i,occ_act(i) - end do - - integer :: nmx - real*8 :: xmx - do i=1,n_act_orb - ! largest element of the eigenvector should be positive - xmx=0.D0 - nmx=0 - do j=1,n_act_orb - if (abs(natOrbsCI(j,i)).gt.xmx) then - nmx=j - xmx=abs(natOrbsCI(j,i)) - end if - end do - xmx=sign(1.D0,natOrbsCI(nmx,i)) - do j=1,n_act_orb - natOrbsCI(j,i)*=xmx - end do - - write(6,*) ' Eigenvector No ',i - write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb) - end do - end if - -END_PROVIDER - - -BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] - implicit none - BEGIN_DOC - ! 4-index transformation of 2part matrices - END_DOC - integer :: i,j,k,l,p,q - real*8 :: d(n_act_orb) - - ! index per index - ! first quarter - P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:) - - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(p,j,k,l)=d(p) - end do - end do - end do - end do - ! 2nd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(j,p,k,l)=d(p) - end do - end do - end do - end do - ! 3rd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(j,k,p,l)=d(p) - end do - end do - end do - end do - ! 4th quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx_no(j,k,l,p)=d(p) - end do - end do - end do - end do - -END_PROVIDER - - - -BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)] - implicit none - BEGIN_DOC - ! Transformed one-e integrals - END_DOC - integer :: i,j, p, q - real*8 :: d(n_act_orb) - one_ints_no(:,:)=mo_one_e_integrals(:,:) - - ! 1st half-trf - do j=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - one_ints_no(list_act(p),j)=d(p) - end do - end do - - ! 2nd half-trf - do j=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - do q=1,n_act_orb - d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - one_ints_no(j,list_act(p))=d(p) - end do - end do -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, NatOrbsCI_mos, (mo_num, mo_num) ] - implicit none - BEGIN_DOC - ! Rotation matrix from current MOs to the CI natural MOs - END_DOC - integer :: p,q - - NatOrbsCI_mos(:,:) = 0.d0 - - do q = 1,mo_num - NatOrbsCI_mos(q,q) = 1.d0 - enddo - - do q = 1,n_act_orb - do p = 1,n_act_orb - NatOrbsCI_mos(list_act(p),list_act(q)) = natorbsCI(p,q) - enddo - enddo -END_PROVIDER - - -BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)] - implicit none - BEGIN_DOC -! FCI natural orbitals - END_DOC - - call dgemm('N','N', ao_num,mo_num,mo_num,1.d0, & - mo_coef, size(mo_coef,1), & - NatOrbsCI_mos, size(NatOrbsCI_mos,1), 0.d0, & - NatOrbsFCI, size(NatOrbsFCI,1)) -END_PROVIDER - diff --git a/src/casscf/neworbs.irp.f b/src/casscf/neworbs.irp.f deleted file mode 100644 index 06a89318..00000000 --- a/src/casscf/neworbs.irp.f +++ /dev/null @@ -1,221 +0,0 @@ -BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] - implicit none - BEGIN_DOC - ! Single-excitation matrix - END_DOC - - integer :: i,j - - do i=1,nMonoEx+1 - do j=1,nMonoEx+1 - SXmatrix(i,j)=0.D0 - end do - end do - - do i=1,nMonoEx - SXmatrix(1,i+1)=gradvec2(i) - SXmatrix(1+i,1)=gradvec2(i) - end do - - do i=1,nMonoEx - do j=1,nMonoEx - SXmatrix(i+1,j+1)=hessmat2(i,j) - SXmatrix(j+1,i+1)=hessmat2(i,j) - end do - end do - - do i = 1, nMonoEx - SXmatrix(i+1,i+1) += level_shift_casscf - enddo - if (bavard) then - do i=2,nMonoEx - write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i) - end do - end if - - -END_PROVIDER - - BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)] -&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)] - implicit none - BEGIN_DOC - ! Eigenvectors/eigenvalues of the single-excitation matrix - END_DOC - call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1) - if (bavard) then - write(6,*) ' SXdiag : lowest 5 eigenvalues ' - write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1) - if(nmonoex.gt.0)then - write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2) - write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3) - write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4) - write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5) - endif - write(6,*) - write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1) - endif -END_PROVIDER - - BEGIN_PROVIDER [real*8, energy_improvement] - implicit none - if(state_following_casscf)then - energy_improvement = SXeigenval(best_vector_ovrlp_casscf) - else - energy_improvement = SXeigenval(1) - endif - END_PROVIDER - - - - BEGIN_PROVIDER [ integer, best_vector_ovrlp_casscf ] -&BEGIN_PROVIDER [ double precision, best_overlap_casscf ] - implicit none - integer :: i - double precision :: c0 - best_overlap_casscf = 0.D0 - best_vector_ovrlp_casscf = -1000 - do i=1,nMonoEx+1 - if (SXeigenval(i).lt.0.D0) then - if (abs(SXeigenvec(1,i)).gt.best_overlap_casscf) then - best_overlap_casscf=abs(SXeigenvec(1,i)) - best_vector_ovrlp_casscf = i - end if - end if - end do - if(best_vector_ovrlp_casscf.lt.0)then - best_vector_ovrlp_casscf = minloc(SXeigenval,nMonoEx+1) - endif - c0=SXeigenvec(1,best_vector_ovrlp_casscf) - if (bavard) then - write(6,*) ' SXdiag : eigenvalue for best overlap with ' - write(6,*) ' previous orbitals = ',SXeigenval(best_vector_ovrlp_casscf) - write(6,*) ' weight of the 1st element ',c0 - endif - END_PROVIDER - - BEGIN_PROVIDER [double precision, SXvector, (nMonoEx+1)] - implicit none - BEGIN_DOC - ! Best eigenvector of the single-excitation matrix - END_DOC - integer :: i - double precision :: c0 - c0=SXeigenvec(1,best_vector_ovrlp_casscf) - do i=1,nMonoEx+1 - SXvector(i)=SXeigenvec(i,best_vector_ovrlp_casscf)/c0 - end do - END_PROVIDER - - -BEGIN_PROVIDER [double precision, NewOrbs, (ao_num,mo_num) ] - implicit none - BEGIN_DOC - ! Updated orbitals - END_DOC - integer :: i,j,ialph - - if(state_following_casscf)then - print*,'Using the state following casscf ' - call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & - NatOrbsFCI, size(NatOrbsFCI,1), & - Umat, size(Umat,1), 0.d0, & - NewOrbs, size(NewOrbs,1)) - - level_shift_casscf *= 0.5D0 - level_shift_casscf = max(level_shift_casscf,0.002d0) - !touch level_shift_casscf - else - if(best_vector_ovrlp_casscf.ne.1.and.n_orb_swap.ne.0)then - print*,'Taking the lowest root for the CASSCF' - print*,'!!! SWAPPING MOS !!!!!!' - level_shift_casscf *= 2.D0 - level_shift_casscf = min(level_shift_casscf,0.5d0) - print*,'level_shift_casscf = ',level_shift_casscf - NewOrbs = switch_mo_coef - !mo_coef = switch_mo_coef - !soft_touch mo_coef - !call save_mos_no_occ - !stop - else - level_shift_casscf *= 0.5D0 - level_shift_casscf = max(level_shift_casscf,0.002d0) - !touch level_shift_casscf - call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, & - NatOrbsFCI, size(NatOrbsFCI,1), & - Umat, size(Umat,1), 0.d0, & - NewOrbs, size(NewOrbs,1)) - endif - endif - -END_PROVIDER - -BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ] - implicit none - BEGIN_DOC - ! Orbital rotation matrix - END_DOC - integer :: i,j,indx,k,iter,t,a,ii,tt,aa - logical :: converged - - real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num) - real*8 :: Tmat(mo_num,mo_num) - real*8 :: f - - ! the orbital rotation matrix T - Tmat(:,:)=0.D0 - indx=1 - do i=1,n_core_inact_orb - ii=list_core_inact(i) - do t=1,n_act_orb - tt=list_act(t) - indx+=1 - Tmat(ii,tt)= SXvector(indx) - Tmat(tt,ii)=-SXvector(indx) - end do - end do - do i=1,n_core_inact_orb - ii=list_core_inact(i) - do a=1,n_virt_orb - aa=list_virt(a) - indx+=1 - Tmat(ii,aa)= SXvector(indx) - Tmat(aa,ii)=-SXvector(indx) - end do - end do - do t=1,n_act_orb - tt=list_act(t) - do a=1,n_virt_orb - aa=list_virt(a) - indx+=1 - Tmat(tt,aa)= SXvector(indx) - Tmat(aa,tt)=-SXvector(indx) - end do - end do - - ! Form the exponential - - Tpotmat(:,:)=0.D0 - Umat(:,:) =0.D0 - do i=1,mo_num - Tpotmat(i,i)=1.D0 - Umat(i,i) =1.d0 - end do - iter=0 - converged=.false. - do while (.not.converged) - iter+=1 - f = 1.d0 / dble(iter) - Tpotmat2(:,:) = Tpotmat(:,:) * f - call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, & - Tpotmat2, size(Tpotmat2,1), & - Tmat, size(Tmat,1), 0.d0, & - Tpotmat, size(Tpotmat,1)) - Umat(:,:) = Umat(:,:) + Tpotmat(:,:) - - converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30) - end do -END_PROVIDER - - - diff --git a/src/casscf/reorder_orb.irp.f b/src/casscf/reorder_orb.irp.f deleted file mode 100644 index 3cb90522..00000000 --- a/src/casscf/reorder_orb.irp.f +++ /dev/null @@ -1,70 +0,0 @@ -subroutine reorder_orbitals_for_casscf - implicit none - BEGIN_DOC -! routine that reorders the orbitals of the CASSCF in terms block of core, active and virtual - END_DOC - integer :: i,j,iorb - integer, allocatable :: iorder(:),array(:) - allocate(iorder(mo_num),array(mo_num)) - do i = 1, n_core_orb - iorb = list_core(i) - array(iorb) = i - enddo - - do i = 1, n_inact_orb - iorb = list_inact(i) - array(iorb) = mo_num + i - enddo - - do i = 1, n_act_orb - iorb = list_act(i) - array(iorb) = 2 * mo_num + i - enddo - - do i = 1, n_virt_orb - iorb = list_virt(i) - array(iorb) = 3 * mo_num + i - enddo - - do i = 1, mo_num - iorder(i) = i - enddo - call isort(array,iorder,mo_num) - double precision, allocatable :: mo_coef_new(:,:) - allocate(mo_coef_new(ao_num,mo_num)) - do i = 1, mo_num - mo_coef_new(:,i) = mo_coef(:,iorder(i)) - enddo - mo_coef = mo_coef_new - touch mo_coef - - list_core_reverse = 0 - do i = 1, n_core_orb - list_core(i) = i - list_core_reverse(i) = i - mo_class(i) = "Core" - enddo - - list_inact_reverse = 0 - do i = 1, n_inact_orb - list_inact(i) = i + n_core_orb - list_inact_reverse(i+n_core_orb) = i - mo_class(i+n_core_orb) = "Inactive" - enddo - - list_act_reverse = 0 - do i = 1, n_act_orb - list_act(i) = n_core_inact_orb + i - list_act_reverse(n_core_inact_orb + i) = i - mo_class(n_core_inact_orb + i) = "Active" - enddo - - list_virt_reverse = 0 - do i = 1, n_virt_orb - list_virt(i) = n_core_inact_orb + n_act_orb + i - list_virt_reverse(n_core_inact_orb + n_act_orb + i) = i - mo_class(n_core_inact_orb + n_act_orb + i) = "Virtual" - enddo - touch list_core_reverse list_core list_inact list_inact_reverse list_act list_act_reverse list_virt list_virt_reverse - -end diff --git a/src/casscf/save_energy.irp.f b/src/casscf/save_energy.irp.f deleted file mode 100644 index 8729c5af..00000000 --- a/src/casscf/save_energy.irp.f +++ /dev/null @@ -1,9 +0,0 @@ -subroutine save_energy(E,pt2) - implicit none - BEGIN_DOC -! Saves the energy in |EZFIO|. - END_DOC - double precision, intent(in) :: E(N_states), pt2(N_states) - call ezfio_set_casscf_energy(E(1:N_states)) - call ezfio_set_casscf_energy_pt2(E(1:N_states)+pt2(1:N_states)) -end diff --git a/src/casscf/superci_dm.irp.f b/src/casscf/superci_dm.irp.f deleted file mode 100644 index ee831c35..00000000 --- a/src/casscf/superci_dm.irp.f +++ /dev/null @@ -1,207 +0,0 @@ - BEGIN_PROVIDER [double precision, super_ci_dm, (mo_num,mo_num)] - implicit none - BEGIN_DOC -! density matrix of the super CI matrix, in the basis of NATURAL ORBITALS OF THE CASCI WF -! -! This is obtained from annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 -! -! WARNING ::: in the equation B3.d there is a TYPO with a forgotten MINUS SIGN (see variable mat_tmp_dm_super_ci ) - END_DOC - super_ci_dm = 0.d0 - integer :: i,j,iorb,jorb - integer :: a,aorb,b,borb - integer :: t,torb,v,vorb,u,uorb,x,xorb - double precision :: c0,ci - c0 = SXeigenvec(1,1) - ! equation B3.a of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - ! loop over the core/inact - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(iorb,iorb) = 2.d0 ! first term of B3.a - ! loop over the core/inact - do j = 1, n_core_inact_orb - jorb = list_core_inact(j) - ! loop over the virtual - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(jorb,iorb) += -2.d0 * lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,jorb) ! second term in B3.a - enddo - do t = 1, n_act_orb - torb = list_act(t) - ! thrid term of the B3.a - super_ci_dm(jorb,iorb) += - lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(jorb,torb) * (2.d0 - occ_act(t)) - enddo - enddo - enddo - - ! equation B3.b of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do t = 1, n_act_orb - torb = list_act(t) - super_ci_dm(iorb,torb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) - super_ci_dm(torb,iorb) = c0 * lowest_super_ci_coef_mo(torb,iorb) * (2.d0 - occ_act(t)) - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(iorb,torb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - super_ci_dm(torb,iorb) += - lowest_super_ci_coef_mo(aorb,iorb) * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - enddo - enddo - enddo - - ! equation B3.c of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(aorb,iorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) - super_ci_dm(iorb,aorb) = 2.d0 * c0 * lowest_super_ci_coef_mo(aorb,iorb) - enddo - enddo - - ! equation B3.d of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do t = 1, n_act_orb - torb = list_act(t) - super_ci_dm(torb,torb) = occ_act(t) ! first term of equation B3.d - do x = 1, n_act_orb - xorb = list_act(x) - super_ci_dm(torb,torb) += - occ_act(x) * occ_act(t)* mat_tmp_dm_super_ci(x,x) ! second term involving the ONE-rdm - enddo - do u = 1, n_act_orb - uorb = list_act(u) - - ! second term of equation B3.d - do x = 1, n_act_orb - xorb = list_act(x) - do v = 1, n_act_orb - vorb = list_act(v) - super_ci_dm(torb,uorb) += 2.d0 * P0tuvx_no(v,x,t,u) * mat_tmp_dm_super_ci(v,x) ! second term involving the TWO-rdm - enddo - enddo - - ! third term of equation B3.d - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(torb,uorb) += lowest_super_ci_coef_mo(iorb,torb) * lowest_super_ci_coef_mo(iorb,uorb) * (2.d0 - occ_act(t) - occ_act(u)) - enddo - - enddo - enddo - - ! equation B3.e of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do t = 1, n_act_orb - torb = list_act(t) - do a = 1, n_virt_orb - aorb = list_virt(a) - super_ci_dm(aorb,torb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - super_ci_dm(torb,aorb) += c0 * lowest_super_ci_coef_mo(aorb,torb) * occ_act(t) - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(aorb,torb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) - super_ci_dm(torb,aorb) += lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,torb) * (2.d0 - occ_act(t)) - enddo - enddo - enddo - - ! equation B3.f of the annex B of Roos et. al. Chemical Physics 48 (1980) 157-173 - do a = 1, n_virt_orb - aorb = list_virt(a) - do b = 1, n_virt_orb - borb= list_virt(b) - - ! First term of equation B3.f - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - super_ci_dm(borb,aorb) += 2.d0 * lowest_super_ci_coef_mo(iorb,aorb) * lowest_super_ci_coef_mo(iorb,borb) - enddo - - ! Second term of equation B3.f - do t = 1, n_act_orb - torb = list_act(t) - super_ci_dm(borb,aorb) += lowest_super_ci_coef_mo(torb,aorb) * lowest_super_ci_coef_mo(torb,borb) * occ_act(t) - enddo - enddo - enddo - - END_PROVIDER - - BEGIN_PROVIDER [double precision, superci_natorb, (ao_num,mo_num) -&BEGIN_PROVIDER [double precision, superci_nat_occ, (mo_num) - implicit none - call general_mo_coef_new_as_svd_vectors_of_mo_matrix_eig(super_ci_dm,mo_num,mo_num,mo_num,NatOrbsFCI,superci_nat_occ,superci_natorb) - -END_PROVIDER - - BEGIN_PROVIDER [double precision, mat_tmp_dm_super_ci, (n_act_orb,n_act_orb)] - implicit none - BEGIN_DOC - ! computation of the term in [ ] in the equation B3.d of Roos et. al. Chemical Physics 48 (1980) 157-173 - ! - ! !!!!! WARNING !!!!!! there is a TYPO: a MINUS SIGN SHOULD APPEAR in that term - END_DOC - integer :: a,aorb,i,iorb - integer :: x,xorb,v,vorb - mat_tmp_dm_super_ci = 0.d0 - do v = 1, n_act_orb - vorb = list_act(v) - do x = 1, n_act_orb - xorb = list_act(x) - do a = 1, n_virt_orb - aorb = list_virt(a) - mat_tmp_dm_super_ci(x,v) += lowest_super_ci_coef_mo(aorb,vorb) * lowest_super_ci_coef_mo(aorb,xorb) - enddo - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - ! MARK THE MINUS SIGN HERE !!!!!!!!!!! BECAUSE OF TYPO IN THE ORIGINAL PAPER - mat_tmp_dm_super_ci(x,v) -= lowest_super_ci_coef_mo(iorb,vorb) * lowest_super_ci_coef_mo(iorb,xorb) - enddo - enddo - enddo - END_PROVIDER - - BEGIN_PROVIDER [double precision, lowest_super_ci_coef_mo, (mo_num,mo_num)] - implicit none - integer :: i,j,iorb,jorb - integer :: a, aorb,t, torb - double precision :: sqrt2 - - sqrt2 = 1.d0/dsqrt(2.d0) - do i = 1, nMonoEx - iorb = excit(1,i) - jorb = excit(2,i) - lowest_super_ci_coef_mo(iorb,jorb) = SXeigenvec(i+1,1) - lowest_super_ci_coef_mo(jorb,iorb) = SXeigenvec(i+1,1) - enddo - - ! a_{it} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do t = 1, n_act_orb - torb = list_act(t) - lowest_super_ci_coef_mo(torb,iorb) *= (2.d0 - occ_act(t))**(-0.5d0) - lowest_super_ci_coef_mo(iorb,torb) *= (2.d0 - occ_act(t))**(-0.5d0) - enddo - enddo - - ! a_{ia} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do a = 1, n_virt_orb - aorb = list_virt(a) - lowest_super_ci_coef_mo(aorb,iorb) *= sqrt2 - lowest_super_ci_coef_mo(iorb,aorb) *= sqrt2 - enddo - enddo - - ! a_{ta} of the equation B.2 of Roos et. al. Chemical Physics 48 (1980) 157-173 - do a = 1, n_virt_orb - aorb = list_virt(a) - do t = 1, n_act_orb - torb = list_act(t) - lowest_super_ci_coef_mo(torb,aorb) *= occ_act(t)**(-0.5d0) - lowest_super_ci_coef_mo(aorb,torb) *= occ_act(t)**(-0.5d0) - enddo - enddo - - END_PROVIDER - diff --git a/src/casscf/swap_orb.irp.f b/src/casscf/swap_orb.irp.f deleted file mode 100644 index 5d442157..00000000 --- a/src/casscf/swap_orb.irp.f +++ /dev/null @@ -1,132 +0,0 @@ - BEGIN_PROVIDER [double precision, SXvector_lowest, (nMonoEx)] - implicit none - integer :: i - do i=2,nMonoEx+1 - SXvector_lowest(i-1)=SXeigenvec(i,1) - enddo - END_PROVIDER - - BEGIN_PROVIDER [double precision, thresh_overlap_switch] - implicit none - thresh_overlap_switch = 0.5d0 - END_PROVIDER - - BEGIN_PROVIDER [integer, max_overlap, (nMonoEx)] -&BEGIN_PROVIDER [integer, n_max_overlap] -&BEGIN_PROVIDER [integer, dim_n_max_overlap] - implicit none - double precision, allocatable :: vec_tmp(:) - integer, allocatable :: iorder(:) - allocate(vec_tmp(nMonoEx),iorder(nMonoEx)) - integer :: i - do i = 1, nMonoEx - iorder(i) = i - vec_tmp(i) = -dabs(SXvector_lowest(i)) - enddo - call dsort(vec_tmp,iorder,nMonoEx) - n_max_overlap = 0 - do i = 1, nMonoEx - if(dabs(vec_tmp(i)).gt.thresh_overlap_switch)then - n_max_overlap += 1 - max_overlap(n_max_overlap) = iorder(i) - endif - enddo - dim_n_max_overlap = max(1,n_max_overlap) - END_PROVIDER - - BEGIN_PROVIDER [integer, orb_swap, (2,dim_n_max_overlap)] -&BEGIN_PROVIDER [integer, index_orb_swap, (dim_n_max_overlap)] -&BEGIN_PROVIDER [integer, n_orb_swap ] - implicit none - use bitmasks ! you need to include the bitmasks_module.f90 features - integer :: i,imono,iorb,jorb,j - n_orb_swap = 0 - do i = 1, n_max_overlap - imono = max_overlap(i) - iorb = excit(1,imono) - jorb = excit(2,imono) - if (excit_class(imono) == "c-a" .and.hessmat2(imono,imono).gt.0.d0)then ! core --> active rotation - n_orb_swap += 1 - orb_swap(1,n_orb_swap) = iorb ! core - orb_swap(2,n_orb_swap) = jorb ! active - index_orb_swap(n_orb_swap) = imono - else if (excit_class(imono) == "a-v" .and.hessmat2(imono,imono).gt.0.d0)then ! active --> virtual rotation - n_orb_swap += 1 - orb_swap(1,n_orb_swap) = jorb ! virtual - orb_swap(2,n_orb_swap) = iorb ! active - index_orb_swap(n_orb_swap) = imono - endif - enddo - - integer,allocatable :: orb_swap_tmp(:,:) - allocate(orb_swap_tmp(2,dim_n_max_overlap)) - do i = 1, n_orb_swap - orb_swap_tmp(1,i) = orb_swap(1,i) - orb_swap_tmp(2,i) = orb_swap(2,i) - enddo - - integer(bit_kind), allocatable :: det_i(:),det_j(:) - allocate(det_i(N_int),det_j(N_int)) - logical, allocatable :: good_orb_rot(:) - allocate(good_orb_rot(n_orb_swap)) - integer, allocatable :: index_orb_swap_tmp(:) - allocate(index_orb_swap_tmp(dim_n_max_overlap)) - index_orb_swap_tmp = index_orb_swap - good_orb_rot = .True. - integer :: icount,k - do i = 1, n_orb_swap - if(.not.good_orb_rot(i))cycle - det_i = 0_bit_kind - call set_bit_to_integer(orb_swap(1,i),det_i,N_int) - call set_bit_to_integer(orb_swap(2,i),det_i,N_int) - do j = i+1, n_orb_swap - det_j = 0_bit_kind - call set_bit_to_integer(orb_swap(1,j),det_j,N_int) - call set_bit_to_integer(orb_swap(2,j),det_j,N_int) - icount = 0 - do k = 1, N_int - icount += popcnt(ior(det_i(k),det_j(k))) - enddo - if (icount.ne.4)then - good_orb_rot(i) = .False. - good_orb_rot(j) = .False. - exit - endif - enddo - enddo - icount = n_orb_swap - n_orb_swap = 0 - do i = 1, icount - if(good_orb_rot(i))then - n_orb_swap += 1 - index_orb_swap(n_orb_swap) = index_orb_swap_tmp(i) - orb_swap(1,n_orb_swap) = orb_swap_tmp(1,i) - orb_swap(2,n_orb_swap) = orb_swap_tmp(2,i) - endif - enddo - - if(n_orb_swap.gt.0)then - print*,'n_orb_swap = ',n_orb_swap - endif - do i = 1, n_orb_swap - print*,'imono = ',index_orb_swap(i) - print*,orb_swap(1,i),'-->',orb_swap(2,i) - enddo - END_PROVIDER - - BEGIN_PROVIDER [double precision, switch_mo_coef, (ao_num,mo_num)] - implicit none - integer :: i,j,iorb,jorb - switch_mo_coef = NatOrbsFCI - do i = 1, n_orb_swap - iorb = orb_swap(1,i) - jorb = orb_swap(2,i) - do j = 1, ao_num - switch_mo_coef(j,jorb) = NatOrbsFCI(j,iorb) - enddo - do j = 1, ao_num - switch_mo_coef(j,iorb) = NatOrbsFCI(j,jorb) - enddo - enddo - - END_PROVIDER diff --git a/src/casscf/test_pert_2rdm.irp.f b/src/casscf/test_pert_2rdm.irp.f deleted file mode 100644 index 7c40de0f..00000000 --- a/src/casscf/test_pert_2rdm.irp.f +++ /dev/null @@ -1,29 +0,0 @@ -program test_pert_2rdm - implicit none - read_wf = .True. - touch read_wf -!call get_pert_2rdm - integer :: i,j,k,l,ii,jj,kk,ll - double precision :: accu , get_two_e_integral, integral - accu = 0.d0 - print*,'n_orb_pert_rdm = ',n_orb_pert_rdm - do ii = 1, n_orb_pert_rdm - i = list_orb_pert_rdm(ii) - do jj = 1, n_orb_pert_rdm - j = list_orb_pert_rdm(jj) - do kk = 1, n_orb_pert_rdm - k= list_orb_pert_rdm(kk) - do ll = 1, n_orb_pert_rdm - l = list_orb_pert_rdm(ll) - integral = get_two_e_integral(i,j,k,l,mo_integrals_map) -! if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then -! print*,i,j,k,l -! print*,pert_2rdm_provider(ii,jj,kk,ll) * integral,pert_2rdm_provider(ii,jj,kk,ll), pert_2rdm_provider(ii,jj,kk,ll), integral -! endif - accu += pert_2rdm_provider(ii,jj,kk,ll) * integral - enddo - enddo - enddo - enddo - print*,'accu = ',accu -end diff --git a/src/casscf/tot_en.irp.f b/src/casscf/tot_en.irp.f deleted file mode 100644 index 1d70e087..00000000 --- a/src/casscf/tot_en.irp.f +++ /dev/null @@ -1,101 +0,0 @@ - BEGIN_PROVIDER [real*8, etwo] -&BEGIN_PROVIDER [real*8, eone] -&BEGIN_PROVIDER [real*8, eone_bis] -&BEGIN_PROVIDER [real*8, etwo_bis] -&BEGIN_PROVIDER [real*8, etwo_ter] -&BEGIN_PROVIDER [real*8, ecore] -&BEGIN_PROVIDER [real*8, ecore_bis] - implicit none - integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3 - real*8 :: e_one_all,e_two_all - e_one_all=0.D0 - e_two_all=0.D0 - do i=1,n_core_inact_orb - ii=list_core_inact(i) - e_one_all+=2.D0*mo_one_e_integrals(ii,ii) - do j=1,n_core_inact_orb - jj=list_core_inact(j) - e_two_all+=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i) - end do - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_inact_orb - do u=1,n_act_orb - uu=list_act(u) - u3=u+n_core_inact_orb - e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) & - -bielec_PQxx(tt,ii,i,u3)) - end do - end do - end do - do t=1,n_act_orb - tt=list_act(t) - do u=1,n_act_orb - uu=list_act(u) - e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu) - do v=1,n_act_orb - v3=v+n_core_inact_orb - do x=1,n_act_orb - x3=x+n_core_inact_orb - e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxx(tt,uu,v3,x3) - end do - end do - end do - end do - ecore =nuclear_repulsion - ecore_bis=nuclear_repulsion - do i=1,n_core_inact_orb - ii=list_core_inact(i) - ecore +=2.D0*mo_one_e_integrals(ii,ii) - ecore_bis+=2.D0*mo_one_e_integrals(ii,ii) - do j=1,n_core_inact_orb - jj=list_core_inact(j) - ecore +=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i) - ecore_bis+=2.D0*bielec_PxxQ(ii,i,j,jj)-bielec_PxxQ(ii,j,j,ii) - end do - end do - eone =0.D0 - eone_bis=0.D0 - etwo =0.D0 - etwo_bis=0.D0 - etwo_ter=0.D0 - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_inact_orb - do u=1,n_act_orb - uu=list_act(u) - u3=u+n_core_inact_orb - eone +=D0tu(t,u)*mo_one_e_integrals(tt,uu) - eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu) - do i=1,n_core_inact_orb - ii=list_core_inact(i) - eone +=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) & - -bielec_PQxx(tt,ii,i,u3)) - eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQ(tt,u3,i,ii) & - -bielec_PxxQ(tt,i,i,uu)) - end do - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_inact_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_inact_orb - real*8 :: h1,h2,h3 - h1=bielec_PQxx(tt,uu,v3,x3) - h2=bielec_PxxQ(tt,u3,v3,xx) - h3=bielecCI(t,u,v,xx) - etwo +=P0tuvx(t,u,v,x)*h1 - etwo_bis+=P0tuvx(t,u,v,x)*h2 - etwo_ter+=P0tuvx(t,u,v,x)*h3 - if ((h1.ne.h2).or.(h1.ne.h3)) then - write(6,9901) t,u,v,x,h1,h2,h3 - 9901 format('aie: ',4I4,3E20.12) - end if - end do - end do - end do - end do - -END_PROVIDER - - From 827a3976fca97e46dc15aa7ecf378da66703914d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Oct 2020 16:43:17 +0100 Subject: [PATCH 11/14] Fix travis --- src/fci/40.fci.bats | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index f86c7488..89effe4f 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -59,7 +59,7 @@ function run_stoch() { @test "HCO" { # 12.2868s qp set_file hco.ezfio - run -113.390292448857 3.e-4 100000 + run -113.389297812482 6.e-4 100000 } @test "H2O2" { # 12.9214s From 9f49da3936e8b83e60d6de8f0d7d236f0678f3b6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 8 Nov 2020 16:15:32 +0100 Subject: [PATCH 12/14] Fixed rPT2 broken in prev commit --- src/cipsi/selection.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 6c1adfde..fd434c1e 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -835,7 +835,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d alpha_h_psi = mat(istate, p1, p2) do jstate=1,N_states - pt2_data % overlap(jstate,istate) += coef(jstate) * alpha_h_psi + pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) enddo pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi From 750c227e23433e37ef27c2f231a18972e2c4a8dd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 8 Nov 2020 16:52:39 +0100 Subject: [PATCH 13/14] Introduced pt2_min_parallel_tasks --- src/cipsi/pt2_stoch_routines.irp.f | 7 +++++-- src/perturbation/EZFIO.cfg | 6 ++++++ 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 31f27e1d..0966039b 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -17,9 +17,12 @@ END_PROVIDER pt2_F(:) = int(sqrt(float(pt2_n_tasks_max))) do i=1,pt2_n_0(1+pt2_N_teeth/4) - pt2_F(i) = pt2_n_tasks_max + pt2_F(i) = pt2_n_tasks_max*pt2_min_parallel_tasks enddo - do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/4), N_det_generators + do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/4), pt2_n_0(pt2_N_teeth-pt2_N_teeth/10) + pt2_F(i) = pt2_min_parallel_tasks + enddo + do i=1+pt2_n_0(pt2_N_teeth-pt2_N_teeth/10), N_det_generators pt2_F(i) = 1 enddo diff --git a/src/perturbation/EZFIO.cfg b/src/perturbation/EZFIO.cfg index 168c3b0d..f4bd8e65 100644 --- a/src/perturbation/EZFIO.cfg +++ b/src/perturbation/EZFIO.cfg @@ -16,6 +16,12 @@ doc: The selection process stops when the largest variance (for all the state) i interface: ezfio,provider,ocaml default: 0.0 +[pt2_min_parallel_tasks] +type: integer +doc: Minimum number of tasks in PT2 calculation +interface: ezfio,provider,ocaml +default: 1 + [pt2_relative_error] type: Normalized_float doc: Stop stochastic |PT2| when the relative error is smaller than `pT2_relative_error` From 96f26a3516cb8a74332bf836094afa6247b2b637 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 8 Nov 2020 17:11:27 +0100 Subject: [PATCH 14/14] Fixed threshold_from_pt2 --- src/davidson/diagonalization_hs2_dressed.irp.f | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index aa748628..bd339367 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -181,7 +181,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ endif overlap = 0.d0 - PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse threshold_davidson_pt2 + PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse threshold_davidson_pt2 threshold_davidson_from_pt2 call write_time(6) write(6,'(A)') '' @@ -600,7 +600,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ if ((itertot>1).and.(iter == 1)) then - !don't print + !don't print continue else write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:3,1:N_st) @@ -608,9 +608,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ! Check convergence if (iter > 1) then + if (threshold_davidson_from_pt2) then converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson_pt2 - endif - + else + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + endif do k=1,N_st if (residual_norm(k) > 1.e8) then