From 6e6a8ac82ac3294ec255dfaca08cfece8bc5e5d5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Feb 2016 10:56:39 +0100 Subject: [PATCH 1/2] Corrected bug in spherical guess MOs --- src/Determinants/connected_to_ref.irp.f | 5 +- src/MO_Basis/ao_ortho_canonical.irp.f | 150 ++++++++++++++++++------ src/Utils/LinearAlgebra.irp.f | 1 - 3 files changed, 116 insertions(+), 40 deletions(-) diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index dc7698b5..b93b18b6 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -95,12 +95,13 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) exit endif enddo - i += 1 - if (i > N_det) then + if (i >= N_det) then return endif + i += 1 + !DIR$ FORCEINLINE do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or. & diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 62a72149..1916bfe4 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -1,3 +1,87 @@ + BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num_align,ao_num)] +&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ] + implicit none + BEGIN_DOC +! matrix of the coefficients of the mos generated by the +! orthonormalization by the S^{-1/2} canonical transformation of the aos +! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital + END_DOC + integer :: i + integer, external :: ao_power_index + integer :: ibegin,j,k + ao_cart_to_sphe_coef(:,:) = 0.d0 + ! Assume order provided by ao_power_index + i = 1 + ao_cart_to_sphe_num = 0 + do while (i <= ao_num) + select case ( ao_l(i) ) + case (0) + ao_cart_to_sphe_coef(i,i) = 1.d0 + i += 1 + ao_cart_to_sphe_num += 1 + BEGIN_TEMPLATE + case ($SHELL) + if (ao_power(i,1) == $SHELL) then + do k=1,size(cart_to_sphe_$SHELL,2) + do j=1,size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k) + enddo + enddo + i += size(cart_to_sphe_$SHELL,1) + ao_cart_to_sphe_num += size(cart_to_sphe_$SHELL,2) + endif + SUBST [ SHELL ] + 1;; + 2;; + 3;; + 4;; + 5;; + 6;; + 7;; + 8;; + 9;; + END_TEMPLATE + case default + stop 'Error in ao_cart_to_sphe' + end select + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_cart_to_sphe_num,ao_cart_to_sphe_num) ] + implicit none + BEGIN_DOC + ! AO overlap matrix in the spherical basis set + END_DOC + double precision, allocatable :: S(:,:) + allocate (S(ao_cart_to_sphe_num,ao_num)) + + call dgemm('T','N',ao_cart_to_sphe_num,ao_num,ao_num, 1.d0, & + ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), & + ao_overlap,size(ao_overlap,1), 0.d0, & + S, size(S,1)) + + call dgemm('N','N',ao_cart_to_sphe_num,ao_cart_to_sphe_num,ao_num, 1.d0, & + S, size(S,1), & + ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), 0.d0, & + ao_cart_to_sphe_overlap,size(ao_cart_to_sphe_overlap,1)) + + deallocate(S) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_num) ] + implicit none + BEGIN_DOC + ! AO_cart_to_sphe_coef^(-1) + END_DOC + + call get_pseudo_inverse(ao_cart_to_sphe_coef,ao_num,ao_cart_to_sphe_num, & + ao_cart_to_sphe_inv, size(ao_cart_to_sphe_coef,1)) +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ] implicit none @@ -8,47 +92,39 @@ END_DOC integer :: i ao_ortho_canonical_coef(:,:) = 0.d0 + do i=1,ao_num + ao_ortho_canonical_coef(i,i) = 1.d0 + enddo + if (ao_cartesian) then - do i=1,ao_num - ao_ortho_canonical_coef(i,i) = 1.d0 - enddo + + ao_ortho_canonical_num = ao_num + call ortho_canonical(ao_overlap,size(ao_overlap,1), & + ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & + ao_ortho_canonical_num) + + else - integer, external :: ao_power_index - integer :: ibegin,j,k - ! Assume order provided by ao_power_index - i = 1 - do while (i <= ao_num) - select case ( ao_l(i) ) - case (0) - ao_ortho_canonical_coef(i,i) = 1.d0 - i += 1 - BEGIN_TEMPLATE - case ($SHELL) - if (ao_power(i,1) == $SHELL) then - do k=1,size(cart_to_sphe_$SHELL,2) - do j=1,size(cart_to_sphe_$SHELL,1) - ao_ortho_canonical_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k) - enddo - enddo - i += size(cart_to_sphe_$SHELL,1) - endif - SUBST [ SHELL ] - 1;; - 2;; - 3;; - 4;; - 5;; - 6;; - 7;; - 8;; - 9;; - END_TEMPLATE - case default - stop 'Error in sphe_to_cart' - end select + + double precision, allocatable :: S(:,:) + + allocate(S(ao_cart_to_sphe_num,ao_cart_to_sphe_num)) + S = 0.d0 + do i=1,ao_cart_to_sphe_num + S(i,i) = 1.d0 enddo + + ao_ortho_canonical_num = ao_cart_to_sphe_num + call ortho_canonical (ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & + ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num) + + call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, & + ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & + S, size(S,1), & + 0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1)) + + deallocate(S) endif - call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num) END_PROVIDER BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index bcfd43b8..dfa5d982 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -88,7 +88,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) endif enddo do i=m+1,n - print *, D(i) D(i) = 0.d0 enddo From ceb4ef8d08b027084870e56f1ee45bfb70248798 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 4 Feb 2016 12:43:24 +0100 Subject: [PATCH 2/2] Updated tests --- tests/bats/qp.bats | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index ab686d71..71604759 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -100,11 +100,11 @@ function run_FCI() { } @test "SCF H2O cc-pVDZ" { - run_HF h2o.ezfio -76.0270218692377 + run_HF h2o.ezfio -76.0230273511649 } @test "FCI H2O cc-pVDZ" { - run_FCI h2o.ezfio 10000 -76.2382562245107 -76.2433933430971 + run_FCI h2o.ezfio 10000 -76.2315960087713 -76.237428352199 } @test "CAS_SD H2O cc-pVDZ" { @@ -113,10 +113,10 @@ function run_FCI() { ezfio set_file $INPUT ezfio set perturbation do_pt2_end False ezfio set determinants n_det_max 1000 - qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-25]" + qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-23]" qp_run cas_sd_selected $INPUT energy="$(ezfio get cas_sd energy)" - eq $energy -76.2219816183745 1.E-6 + eq $energy -76.2085511296300 1.E-6 } @test "MRCC H2O cc-pVDZ" { @@ -128,7 +128,8 @@ function run_FCI() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2313853324571 1.E-3 + eq $energy -76.2165731870755 1.E-3 + } @@ -138,11 +139,11 @@ function run_FCI() { } @test "SCF H2O VDZ pseudo" { - run_HF h2o_pseudo.ezfio -16.9483703904632 + run_HF h2o_pseudo.ezfio -16.9457263818675 } @test "FCI H2O VDZ pseudo" { - run_FCI h2o_pseudo.ezfio 2000 -17.1550015533975 -17.1645044211228 + run_FCI h2o_pseudo.ezfio 2000 -17.1476897854369 -17.1598005211929 } #=== Convert