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