From 7eb7e2134afe1bb26075c010be0cfd95d6d4ee60 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 19 Jun 2017 20:38:28 +0200 Subject: [PATCH 1/3] Working on dressing --- ocaml/qp_run.ml | 3 +- plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg | 6 + .../LinearSystem.irp.f | 38 ++-- .../SCF_dressed.irp.f | 19 +- .../Hartree_Fock_SlaterDressed/at_nucl.irp.f | 20 ++- .../Hartree_Fock_SlaterDressed/dressing.irp.f | 164 ++++++++++++++++++ .../integrals.irp.f | 4 +- .../ao_bielec_integrals_in_map_slave.irp.f | 2 - src/MO_Basis/ao_ortho_canonical.irp.f | 17 +- src/MO_Basis/mos.irp.f | 44 ++--- src/Utils/LinearAlgebra.irp.f | 40 ++++- src/ZMQ/utils.irp.f | 4 +- 12 files changed, 282 insertions(+), 79 deletions(-) create mode 100644 plugins/Hartree_Fock_SlaterDressed/dressing.irp.f diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml index 17d31cea..5a656d2d 100644 --- a/ocaml/qp_run.ml +++ b/ocaml/qp_run.ml @@ -134,7 +134,8 @@ let run slave exe ezfio_file = let duration = Time.diff (Time.now()) time_start |> Core.Span.to_string in Printf.printf "Wall time : %s\n\n" duration; - exit exit_code + if (exit_code <> 0) then + exit exit_code let spec = let open Command.Spec in diff --git a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg index 8e4b2847..6a7d5fdc 100644 --- a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg +++ b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg @@ -10,5 +10,11 @@ doc: Exponents of the additional Slater functions size: (nuclei.nucl_num,mo_basis.mo_tot_num) interface: ezfio, provider +[projector] +type: double precision +doc: Orthogonal AO basis +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio, provider + diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f index 490e344c..4ccd09b0 100644 --- a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f @@ -13,7 +13,8 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ] cusp_A(A,B) -= slater_value_at_nucl(B,A) ! Projector do mu=1,mo_tot_num - cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A) + cusp_A(A,B) += AO_orthoSlaOverlap_matrix(mu,B) * ao_ortho_value_at_nucl(mu,A) +! cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A) enddo enddo enddo @@ -36,18 +37,29 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ] - implicit none - BEGIN_DOC - ! Equations to solve : A.C = B - END_DOC - - double precision, allocatable :: AF(:,:) - integer :: info - allocate ( AF(nucl_num,nucl_num) ) - - call get_pseudo_inverse(cusp_A,nucl_num,nucl_num,AF,size(AF,1)) - call dgemm('N','N',nucl_num,mo_tot_num,nucl_num,1.d0, & - AF,size(AF,1), cusp_B, size(cusp_B,1), 0.d0, cusp_C, size(cusp_C,1)) + implicit none + BEGIN_DOC + ! Equations to solve : A.C = B + END_DOC + + integer :: info + integer :: ipiv(nucl_num) + double precision, allocatable :: AF(:,:) + allocate ( AF(nucl_num,nucl_num) ) + + cusp_C(1:nucl_num,1:mo_tot_num) = cusp_B(1:nucl_num,1:mo_tot_num) + AF(1:nucl_num,1:nucl_num) = cusp_A(1:nucl_num,1:nucl_num) + call dgetrf(nucl_num,nucl_num,AF,size(AF,1),ipiv,info) + if (info /= 0) then + print *, info + stop 'dgetrf failed' + endif + call dgetrs('N',nucl_num,mo_tot_num,AF,size(AF,1),ipiv,cusp_C,size(cusp_C,1),info) + if (info /= 0) then + print *, info + stop 'dgetrs failed' + endif + END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f index a8c1351a..d0998a61 100644 --- a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f @@ -76,21 +76,26 @@ subroutine run double precision :: EHF integer :: i_it, i, j, k - mo_label = "CuspDressed" + mo_label = 'None' - print *, HF_energy - - call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_C) +! print *, HF_energy do i=1,ao_num print *, mo_coef(i,1), cusp_corrected_mos(i,1) enddo mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num) - TOUCH mo_coef + SOFT_TOUCH mo_coef slater_coef +! call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef) +! call save_mos + print *, 'ci' + print *, mo_coef(1:ao_num,1) + print *, 'cAi' + print *, slater_coef - EHF = HF_energy - print *, HF_energy + +! EHF = HF_energy +! print *, HF_energy ! call Roothaan_Hall_SCF end diff --git a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f index 910c5a16..afb85b64 100644 --- a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f @@ -28,14 +28,26 @@ BEGIN_PROVIDER [ double precision , ao_value_at_nucl, (ao_num,nucl_num) ] enddo END_PROVIDER +BEGIN_PROVIDER [ double precision, ao_ortho_value_at_nucl, (ao_num,nucl_num) ] + implicit none + BEGIN_DOC +! Values of the molecular orbitals at the nucleus + END_DOC + + call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, & + ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & + ao_value_at_nucl, size(ao_value_at_nucl,1), & + 0.d0, ao_ortho_value_at_nucl,size(ao_ortho_value_at_nucl,1)) +END_PROVIDER + BEGIN_PROVIDER [ double precision, mo_value_at_nucl, (mo_tot_num,nucl_num) ] implicit none BEGIN_DOC ! Values of the molecular orbitals at the nucleus END_DOC - call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, & - mo_coef_transp, size(mo_coef_transp,1), & + call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0, & + mo_coef, size(mo_coef,1), & ao_value_at_nucl, size(ao_value_at_nucl,1), & 0.d0, mo_value_at_nucl, size(mo_value_at_nucl,1)) END_PROVIDER @@ -54,10 +66,6 @@ BEGIN_PROVIDER [ double precision , slater_value_at_nucl, (nucl_num,nucl_num) ] x = nucl_coord(i,1) - nucl_coord(k,1) y = nucl_coord(i,2) - nucl_coord(k,2) z = nucl_coord(i,3) - nucl_coord(k,3) - -! expo = slater_expo(i)*slater_expo(i)*((x*x) + (y*y) + (z*z)) -! if (expo > 160.d0) cycle -! expo = dsqrt(expo) expo = slater_expo(i) * dsqrt((x*x) + (y*y) + (z*z)) slater_value_at_nucl(i,k) = dexp(-expo) * slater_normalization(i) enddo diff --git a/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f new file mode 100644 index 00000000..f13f4a13 --- /dev/null +++ b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f @@ -0,0 +1,164 @@ +BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral_dressing, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! Dressing of the core hamiltonian in the orthogonal AO basis set + END_DOC + + integer :: i,j,k + integer :: mu, nu, lambda, A + double precision :: tmp + + ao_ortho_mono_elec_integral_dressing = 0.d0 + i = idx_dressing + do mu=1,ao_num + if (dabs(mo_coef_in_ao_ortho_basis(mu,i)) > 1.d-5) then + do A=1,nucl_num + tmp = 0.d0 + do nu=1,ao_num + tmp += AO_orthoSlaOverlap_matrix(nu,A) * ao_ortho_mono_elec_integral(mu,nu) + enddo + ao_ortho_mono_elec_integral_dressing(mu,mu) += cusp_C(A,i) * (AO_orthoSlaH_matrix(mu,A) - tmp) + enddo + ao_ortho_mono_elec_integral_dressing(mu,mu) *= 1.d0/mo_coef_in_ao_ortho_basis(mu,i) + endif + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! h core in orthogonal AO basis + END_DOC + double precision, allocatable :: T(:,:) + allocate(T(ao_num,ao_num)) + call dgemm('T','N',ao_num,ao_num,ao_num,1.d0, & + ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & + ao_mono_elec_integral, size(ao_mono_elec_integral,1), & + 0.d0, T, size(T,1)) + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + T, size(T,1), & + ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & + 0.d0, ao_ortho_mono_elec_integral, size(ao_ortho_mono_elec_integral,1)) + deallocate(T) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_dressing, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! Dressing of the core hamiltonian in the AO basis set + END_DOC + call ao_ortho_cano_to_ao(ao_ortho_mono_elec_integral_dressing,size(ao_ortho_mono_elec_integral_dressing,1),& + ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1)) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_dressing, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Dressing of the core hamiltonian in the MO basis set + END_DOC + + call ao_to_mo(ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1),& + mo_mono_elec_integral_dressing,size(mo_mono_elec_integral_dressing,1)) +END_PROVIDER + + +BEGIN_PROVIDER [ integer, idx_dressing ] + implicit none + BEGIN_DOC + ! Index of the MO which is being dressed + END_DOC + idx_dressing = 1 +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Dressing core hamiltonian in the AO basis set + END_DOC + integer :: i,j + double precision, allocatable :: F(:,:), M(:,:) + allocate(F(mo_tot_num_align,mo_tot_num),M(ao_num,mo_tot_num)) + + logical :: oneshot + +! oneshot = .True. + oneshot = .False. + + if (oneshot) then + cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num) + slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num) + return + endif + + do idx_dressing=1,mo_tot_num + + if (idx_dressing>1) then + TOUCH idx_dressing + endif + + do j=1,mo_tot_num + do i=1,mo_tot_num + F(i,j) = Fock_matrix_mo(i,j) + enddo + enddo + + do j=1,mo_tot_num + do i=1,ao_num + M(i,j) = mo_coef(i,j) + enddo + enddo + + integer :: it + do it=1,128 + +! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) +! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num) +! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num) +! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) +! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) +! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num) +! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) +! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) +! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num) +! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) +! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num) + do j=1,mo_tot_num + do i=1,mo_tot_num + Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j) + enddo + enddo + do i=1,mo_tot_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo + double precision :: conv + conv = 0.d0 + do j=1,mo_tot_num + do i=1,mo_tot_num + if (i==j) cycle + conv = max(conv,Fock_matrix_mo(i,j)) + enddo + enddo + TOUCH Fock_matrix_mo Fock_matrix_diag_mo + + mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) + TOUCH mo_coef +! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) +! print *, '-----' + print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv) + if (conv < 1.d-5) exit + + enddo + cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing) + slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing) + enddo + + idx_dressing = 1 + mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num) + soft_TOUCH mo_coef idx_dressing slater_coef + +END_PROVIDER + + diff --git a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f index 1111055b..fa4530a9 100644 --- a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f @@ -502,8 +502,8 @@ BEGIN_PROVIDER [ double precision, MOSla$X_matrix, (mo_tot_num, nucl_num) ] BEGIN_DOC ! END_DOC - call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, & - mo_coef_transp, size(mo_coef_transp,1), & + call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0, & + mo_coef, size(mo_coef,1), & GauSla$X_matrix, size(GauSla$X_matrix,1), & 0.d0, MOSla$X_matrix, size(MOSla$X_matrix,1)) END_PROVIDER diff --git a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f index 644f89c2..f24e4974 100644 --- a/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f +++ b/src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f @@ -56,7 +56,6 @@ subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, stop 'error' endif -! Activate is zmq_socket_push is a REQ IRP_IF ZMQ_PUSH IRP_ELSE integer :: idummy @@ -189,7 +188,6 @@ subroutine ao_bielec_integrals_in_map_collector rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0) -! Activate if zmq_socket_pull is a REP IRP_IF ZMQ_PUSH IRP_ELSE rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0) diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index b0eabfbd..2184ce4a 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -76,19 +76,20 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_ ! 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)) + call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),& + ao_num,ao_cart_to_sphe_num, & + ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1)) END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num,ao_num)] +BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num_align,ao_num)] implicit none BEGIN_DOC ! ao_ortho_canonical_coef^(-1) END_DOC - call get_pseudo_inverse(ao_ortho_canonical_coef,ao_num,ao_num, & - ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef,1)) + call get_inverse(ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),& + ao_num, ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1)) END_PROVIDER BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] @@ -100,11 +101,15 @@ END_PROVIDER ! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital END_DOC integer :: i - ao_ortho_canonical_coef(:,:) = 0.d0 + ao_ortho_canonical_coef = 0.d0 do i=1,ao_num ao_ortho_canonical_coef(i,i) = 1.d0 enddo +!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) +!ao_ortho_canonical_num=ao_num +!return + if (ao_cartesian) then ao_ortho_canonical_num = ao_num diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 260af392..6906874b 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -72,6 +72,8 @@ BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_ implicit none BEGIN_DOC ! MO coefficients in orthogonalized AO basis + ! + ! C^(-1).C_mo END_DOC call dgemm('T','N',ao_num,mo_tot_num,ao_num,1.d0, & ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),& @@ -155,6 +157,8 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) implicit none BEGIN_DOC ! Transform A from the AO basis to the MO basis + ! + ! C.A_ao.Ct END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_ao(LDA_ao) @@ -181,6 +185,8 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) implicit none BEGIN_DOC ! Transform A from the MO basis to the AO basis + ! + ! (S.C).A_mo.(S.C)t END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_mo(LDA_mo) @@ -273,6 +279,8 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) implicit none BEGIN_DOC ! Transform A from the AO basis to the orthogonal AO basis + ! + ! C^(-1).A_ao.Ct^(-1) END_DOC integer, intent(in) :: LDA_ao,LDA double precision, intent(in) :: A_ao(LDA_ao,*) @@ -282,13 +290,13 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) allocate ( T(ao_num_align,ao_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - call dgemm('T','N', ao_num, ao_num, ao_num, & + call dgemm('N','N', ao_num, ao_num, ao_num, & 1.d0, & ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1), & A_ao,LDA_ao, & 0.d0, T, ao_num_align) - call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, & + call dgemm('N','T', ao_num, ao_num, ao_num, 1.d0, & T, size(T,1), & ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),& 0.d0, A, LDA) @@ -296,35 +304,3 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) deallocate(T) end -subroutine mo_to_ao_ortho_cano(A_mo,LDA_mo,A_ao,LDA_ao) - implicit none - BEGIN_DOC - ! Transform A from the AO orthogonal basis to the AO basis - END_DOC - integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_mo(LDA_mo) - double precision, intent(out) :: A_ao(LDA_ao) - double precision, allocatable :: T(:,:), SC(:,:) - - allocate ( SC(ao_num_align,mo_tot_num) ) - allocate ( T(mo_tot_num_align,ao_num) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, ao_overlap,size(ao_overlap,1), & - ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & - 0.d0, SC, ao_num_align) - - call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & - 1.d0, A_mo,LDA_mo, & - SC, size(SC,1), & - 0.d0, T, mo_tot_num_align) - - call dgemm('N','N', ao_num, ao_num, mo_tot_num, & - 1.d0, SC,size(SC,1), & - T, mo_tot_num_align, & - 0.d0, A_ao, LDA_ao) - - deallocate(T,SC) -end - diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 92bcdf53..f09e98ff 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -72,8 +72,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) double precision, allocatable :: S_half(:,:) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j -!call ortho_lowdin(overlap,LDA,N,C,LDC,m) -!return if (n < 2) then return @@ -270,14 +268,42 @@ end -subroutine get_pseudo_inverse(A,m,n,C,LDA) +subroutine get_inverse(A,LDA,m,C,LDC) + implicit none + BEGIN_DOC + ! Returns the inverse of the square matrix A + END_DOC + integer, intent(in) :: m, LDA, LDC + double precision, intent(in) :: A(LDA,m) + double precision, intent(out) :: C(LDC,m) + + integer :: info,lwork + integer, allocatable :: ipiv(:) + double precision,allocatable :: work(:) + allocate (ipiv(ao_num), work(ao_num*ao_num)) + lwork = size(work) + C(1:m,1:m) = A(1:m,1:m) + call dgetrf(m,m,C,size(C,1),ipiv,info) + if (info /= 0) then + print *, info + stop 'error in inverse (dgetrf)' + endif + call dgetri(m,C,size(C,1),ipiv,work,lwork,info) + if (info /= 0) then + print *, info + stop 'error in inverse (dgetri)' + endif + deallocate(ipiv,work) +end + +subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC) implicit none BEGIN_DOC ! Find C = A^-1 END_DOC - integer, intent(in) :: m,n, LDA + integer, intent(in) :: m,n, LDA, LDC double precision, intent(in) :: A(LDA,n) - double precision, intent(out) :: C(n,m) + double precision, intent(out) :: C(LDC,m) double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:) integer :: info, lwork @@ -304,7 +330,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA) endif do i=1,n - if (dabs(D(i)) > 1.d-6) then + if (D(i)/D(1) > 1.d-10) then D(i) = 1.d0/D(i) else D(i) = 0.d0 @@ -315,7 +341,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA) do i=1,m do j=1,n do k=1,n - C(j,i) += U(i,k) * D(k) * Vt(k,j) + C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j) enddo enddo enddo diff --git a/src/ZMQ/utils.irp.f b/src/ZMQ/utils.irp.f index 4b0cd0c5..e2651073 100644 --- a/src/ZMQ/utils.irp.f +++ b/src/ZMQ/utils.irp.f @@ -228,12 +228,14 @@ function new_zmq_pull_socket() character*(8), external :: zmq_port integer(ZMQ_PTR) :: new_zmq_pull_socket + print *, 'coucou' call omp_set_lock(zmq_lock) + print *, 'pouet' if (zmq_context == 0_ZMQ_PTR) then stop 'zmq_context is uninitialized' endif IRP_IF ZMQ_PUSH - new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL) + new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_PULL) IRP_ELSE new_zmq_pull_socket = f77_zmq_socket(zmq_context, ZMQ_REP) IRP_ENDIF From 3756b0a26aedf831a58ce65d8906dea4a635d738 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Jun 2017 15:18:19 +0200 Subject: [PATCH 2/3] Added mo_energies as diff of expectation values of ionized wfs --- src/Determinants/mo_energy_expval.irp.f | 159 ++++++++++++++++++++++++ src/ZMQ/utils.irp.f | 2 - 2 files changed, 159 insertions(+), 2 deletions(-) create mode 100644 src/Determinants/mo_energy_expval.irp.f diff --git a/src/Determinants/mo_energy_expval.irp.f b/src/Determinants/mo_energy_expval.irp.f new file mode 100644 index 00000000..ea62a441 --- /dev/null +++ b/src/Determinants/mo_energy_expval.irp.f @@ -0,0 +1,159 @@ +BEGIN_PROVIDER [ double precision, mo_energy_expval, (N_states,mo_tot_num,2,2)] + use bitmasks + implicit none + BEGIN_DOC + ! Third index is spin. + ! Fourth index is 1:creation, 2:annihilation + END_DOC + integer :: i,j,k + integer :: ispin, istate + integer :: hp + double precision :: norm_out(N_states) + + integer, parameter :: hole_particle(2) = (/ -1, 1 /) + double precision :: energies(n_states) + + integer(bit_kind), allocatable :: psi_in_out(:,:,:) + double precision, allocatable :: psi_in_out_coef(:,:) + double precision :: E0(N_states), norm + double precision, parameter :: t=1.d-3 + + allocate (psi_in_out(N_int,2,N_det),psi_in_out_coef(N_det,N_states)) + mo_energy_expval = 0.d0 + + psi_in_out_coef(1:N_det,1:N_states) = psi_coef(1:N_det,1:N_states) + psi_in_out(1:N_int,1:2,1:N_det) = psi_det(1:N_int,1:2,1:N_det) + + ! Truncate the wave function + do istate=1,N_states + norm = 0.d0 + do k=1,N_det + if (dabs(psi_in_out_coef(k,istate)) < t) then + psi_in_out_coef(k,istate) = 0.d0 + endif + norm = norm + psi_in_out_coef(k,istate)*psi_in_out_coef(k,istate) + enddo + ASSERT (norm > 0.d0) + norm = 1.d0/dsqrt(norm) + psi_in_out_coef(1:N_det,istate) = psi_in_out_coef(1:N_det,istate) * norm + call au0_h_au0(E0,psi_in_out,psi_in_out_coef,N_det,size(psi_in_out_coef,1)) + enddo + + + do hp=1,2 + do ispin=1,2 + do i=1,mo_tot_num + psi_in_out_coef(1:N_det,1:N_states) = psi_coef(1:N_det,1:N_states) + psi_in_out(1:N_int,1:2,1:N_det) = psi_det(1:N_int,1:2,1:N_det) + call apply_exc_to_psi(i,hole_particle(hp),ispin, & + norm_out,psi_in_out,psi_in_out_coef, N_det,N_det,N_det,N_states) + + ! Truncate the wave function + do istate=1,N_states + norm = 0.d0 + do k=1,N_det + if (dabs(psi_in_out_coef(k,istate)) < t) then + psi_in_out_coef(k,istate) = 0.d0 + endif + norm = norm + psi_in_out_coef(k,istate)*psi_in_out_coef(k,istate) + enddo + if (norm == 0.d0) then + cycle + endif + norm = 1.d0/dsqrt(norm) + psi_in_out_coef(1:N_det,istate) = psi_in_out_coef(1:N_det,istate) * norm + enddo + call au0_h_au0(energies,psi_in_out,psi_in_out_coef,N_det,size(psi_in_out_coef,1)) + mo_energy_expval(1:N_states,i,ispin,hp) = energies(1:N_states) - E0(1:N_states) + print *, i, ispin, real(energies(1)), real(E0(1)) + enddo + enddo + + enddo + mo_energy_expval(1:N_states,1:mo_tot_num,1:2,1) = -mo_energy_expval(1:N_states,1:mo_tot_num,1:2,1) + +END_PROVIDER + + +subroutine au0_h_au0(energies,psi_in,psi_in_coef,ndet,dim_psi_coef) + use bitmasks + implicit none + integer, intent(in) :: ndet,dim_psi_coef + integer(bit_kind), intent(in) :: psi_in(N_int,2,Ndet) + double precision, intent(in) :: psi_in_coef(dim_psi_coef,N_states) + double precision, intent(out) :: energies(N_states) + + integer :: i,j, istate + double precision :: hij,accu + double precision, allocatable :: psi_coef_tmp(:) + + energies(1:N_states) = 0.d0 + do i = 1, Ndet + if(sum(dabs(psi_in_coef(i,1:N_states)))==0.d0) then + cycle + endif + call diag_H_mat_elem_au0_h_au0(psi_in(1,1,i),N_int,hij) + do istate=1,N_states + energies(istate) += psi_in_coef(i,istate) * psi_in_coef(i,istate) * hij + enddo + do j = i+1, Ndet + if(sum(dabs(psi_in_coef(j,1:N_states)))==0.d0) then + cycle + endif + call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij) + hij = hij+hij + do istate=1,N_states + energies(istate) = energies(istate) + psi_in_coef(i,istate) * psi_in_coef(j,istate) * hij + enddo + enddo + enddo +end + +subroutine diag_H_mat_elem_au0_h_au0(det_in,Nint,hii) + use bitmasks + implicit none + BEGIN_DOC + ! Computes for any determinant i. Used for wave functions with an additional electron. + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + double precision, intent(out) :: hii + + integer :: i, j, iorb, jorb + integer :: occ(Nint*bit_kind_size,2) + integer :: elec_num_tab_local(2) + + hii = 0.d0 + call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), Nint) + call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), Nint) + + ! alpha - alpha + do i = 1, elec_num_tab_local(1) + iorb = occ(i,1) + hii += mo_mono_elec_integral(iorb,iorb) + do j = i+1, elec_num_tab_local(1) + jorb = occ(j,1) + hii += mo_bielec_integral_jj_anti(jorb,iorb) + enddo + enddo + + ! beta - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + hii += mo_mono_elec_integral(iorb,iorb) + do j = i+1, elec_num_tab_local(2) + jorb = occ(j,2) + hii += mo_bielec_integral_jj_anti(jorb,iorb) + enddo + enddo + + ! alpha - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + do j = 1, elec_num_tab_local(1) + jorb = occ(j,1) + hii += mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + +end diff --git a/src/ZMQ/utils.irp.f b/src/ZMQ/utils.irp.f index e2651073..7b26ed57 100644 --- a/src/ZMQ/utils.irp.f +++ b/src/ZMQ/utils.irp.f @@ -228,9 +228,7 @@ function new_zmq_pull_socket() character*(8), external :: zmq_port integer(ZMQ_PTR) :: new_zmq_pull_socket - print *, 'coucou' call omp_set_lock(zmq_lock) - print *, 'pouet' if (zmq_context == 0_ZMQ_PTR) then stop 'zmq_context is uninitialized' endif From 5f357ee68b6511ab6f408bd6916cc6a676b332d3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Jun 2017 18:11:49 +0200 Subject: [PATCH 3/3] Working on Slater dressing --- plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg | 8 +- .../LinearSystem.irp.f | 2 - .../SCF_dressed.irp.f | 6 +- .../Hartree_Fock_SlaterDressed/dressing.irp.f | 140 +++++++++--------- .../integrals.irp.f | 11 +- .../Hartree_Fock_SlaterDressed/slater.irp.f | 5 +- src/MO_Basis/ao_ortho_canonical.irp.f | 6 +- src/MO_Basis/mos.irp.f | 6 +- 8 files changed, 101 insertions(+), 83 deletions(-) diff --git a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg index 6a7d5fdc..df55c554 100644 --- a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg +++ b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg @@ -14,7 +14,13 @@ interface: ezfio, provider type: double precision doc: Orthogonal AO basis size: (ao_basis.ao_num,ao_basis.ao_num) -interface: ezfio, provider +interface: ezfio + +[ao_orthoSlaOverlap] +type: double precision +doc: Orthogonal AO basis +size: (ao_basis.ao_num,nuclei.nucl_num) +interface: ezfio diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f index 4ccd09b0..61a5a49c 100644 --- a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f @@ -14,7 +14,6 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ] ! Projector do mu=1,mo_tot_num cusp_A(A,B) += AO_orthoSlaOverlap_matrix(mu,B) * ao_ortho_value_at_nucl(mu,A) -! cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A) enddo enddo enddo @@ -60,6 +59,5 @@ BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ] stop 'dgetrs failed' endif - END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f index d0998a61..8befdb91 100644 --- a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f @@ -86,8 +86,10 @@ subroutine run enddo mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num) SOFT_TOUCH mo_coef slater_coef -! call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef) -! call save_mos + call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef) + call ezfio_set_Hartree_Fock_SlaterDressed_projector(ao_ortho_canonical_coef(1:ao_num,1:ao_num)) + call ezfio_set_Hartree_Fock_SlaterDressed_ao_orthoSlaOverlap(AO_orthoSlaOverlap_matrix) + call save_mos print *, 'ci' print *, mo_coef(1:ao_num,1) print *, 'cAi' diff --git a/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f index f13f4a13..f42a0b0f 100644 --- a/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f @@ -91,74 +91,82 @@ BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num_align,mo_tot_num) cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num) slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num) return + + else + + + do idx_dressing=1,mo_tot_num + + if (idx_dressing>1) then + TOUCH idx_dressing + endif + + do j=1,mo_tot_num + do i=1,mo_tot_num + F(i,j) = Fock_matrix_mo(i,j) + enddo + enddo + + do j=1,mo_tot_num + do i=1,ao_num + M(i,j) = mo_coef(i,j) + enddo + enddo + + integer :: it + do it=1,128 + + ! print *, 'X', ao_ortho_canonical_coef(1:ao_num,1:ao_num) + ! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) + ! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num) + ! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num) + ! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) + ! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num) + ! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) + ! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) + ! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num) + ! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) + ! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num) + ! print *, 'Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) + do j=1,mo_tot_num + do i=1,mo_tot_num + Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j) + enddo + enddo + do i=1,mo_tot_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo + ! print *, 'Dressed Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) + double precision :: conv + conv = 0.d0 + do j=1,mo_tot_num + do i=1,mo_tot_num + if (i==j) cycle + conv = max(conv,Fock_matrix_mo(i,j)) + enddo + enddo + TOUCH Fock_matrix_mo Fock_matrix_diag_mo + + mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) + TOUCH mo_coef + !print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) + !print *, '-----' + print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv) + if (conv < 1.d-5) exit + !stop + + enddo + cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing) + slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing) + enddo + + idx_dressing = 1 + mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num) + soft_TOUCH mo_coef idx_dressing slater_coef + endif - do idx_dressing=1,mo_tot_num - - if (idx_dressing>1) then - TOUCH idx_dressing - endif - - do j=1,mo_tot_num - do i=1,mo_tot_num - F(i,j) = Fock_matrix_mo(i,j) - enddo - enddo - - do j=1,mo_tot_num - do i=1,ao_num - M(i,j) = mo_coef(i,j) - enddo - enddo - - integer :: it - do it=1,128 - -! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) -! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num) -! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num) -! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) -! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) -! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num) -! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) -! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) -! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num) -! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) -! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num) - do j=1,mo_tot_num - do i=1,mo_tot_num - Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j) - enddo - enddo - do i=1,mo_tot_num - Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) - enddo - double precision :: conv - conv = 0.d0 - do j=1,mo_tot_num - do i=1,mo_tot_num - if (i==j) cycle - conv = max(conv,Fock_matrix_mo(i,j)) - enddo - enddo - TOUCH Fock_matrix_mo Fock_matrix_diag_mo - - mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) - TOUCH mo_coef -! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) -! print *, '-----' - print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv) - if (conv < 1.d-5) exit - - enddo - cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing) - slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing) - enddo - - idx_dressing = 1 - mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num) - soft_TOUCH mo_coef idx_dressing slater_coef - END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f index fa4530a9..680f3e1f 100644 --- a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f @@ -344,7 +344,7 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result) ss = k*ss ! Print result - write(*,*) ss +! write(*,*) ss result = 0.d0 end @@ -520,6 +520,7 @@ BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ] END_PROVIDER + SUBST [ X ] Overlap ;; @@ -616,9 +617,9 @@ BEGIN_PROVIDER [ double precision, AO_orthoSlaH_matrix, (ao_num, nucl_num) ] BEGIN_DOC ! END_DOC - call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, & - ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & - GauSlaH_matrix, size(GauSlaH_matrix,1), & - 0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1)) + call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, & + ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & + GauSlaH_matrix, size(GauSlaH_matrix,1), & + 0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1)) END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f index 8a11b9b1..27e4f10d 100644 --- a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f @@ -8,7 +8,10 @@ BEGIN_PROVIDER [ double precision, slater_expo, (nucl_num) ] if (exists) then slater_expo(1:nucl_num) = slater_expo_ezfio(1:nucl_num) else - slater_expo(1:nucl_num) = nucl_charge(1:nucl_num) + integer :: i + do i=1,nucl_num + slater_expo(i) = nucl_charge(i) + enddo call ezfio_set_Hartree_Fock_SlaterDressed_slater_expo_ezfio(slater_expo) endif END_PROVIDER diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 2184ce4a..d77374a8 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -106,9 +106,9 @@ END_PROVIDER ao_ortho_canonical_coef(i,i) = 1.d0 enddo -!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) -!ao_ortho_canonical_num=ao_num -!return +call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) +ao_ortho_canonical_num=ao_num +return if (ao_cartesian) then diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 6906874b..8962dd00 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -75,7 +75,7 @@ BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_ ! ! C^(-1).C_mo END_DOC - call dgemm('T','N',ao_num,mo_tot_num,ao_num,1.d0, & + call dgemm('N','N',ao_num,mo_tot_num,ao_num,1.d0, & ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),& mo_coef, size(mo_coef,1), 0.d0, & mo_coef_in_ao_ortho_basis, size(mo_coef_in_ao_ortho_basis,1)) @@ -290,13 +290,13 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) allocate ( T(ao_num_align,ao_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - call dgemm('N','N', ao_num, ao_num, ao_num, & + call dgemm('T','N', ao_num, ao_num, ao_num, & 1.d0, & ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1), & A_ao,LDA_ao, & 0.d0, T, ao_num_align) - call dgemm('N','T', ao_num, ao_num, ao_num, 1.d0, & + call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, & T, size(T,1), & ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),& 0.d0, A, LDA)