diff --git a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f index 40d75fc4..7c30e175 100644 --- a/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f +++ b/plugins/FOBOCI/diag_fock_inactiv_virt.irp.f @@ -29,7 +29,7 @@ subroutine diag_inactive_virt_and_update_mos label = "Canonical" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.) soft_touch mo_coef @@ -76,7 +76,7 @@ subroutine diag_inactive_virt_new_and_update_mos label = "Canonical" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.) soft_touch mo_coef diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f index 8e2f2e53..6d359696 100644 --- a/plugins/FOBOCI/new_approach.irp.f +++ b/plugins/FOBOCI/new_approach.irp.f @@ -655,7 +655,7 @@ subroutine new_approach integer :: sign sign = -1 - call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign) + call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign,.true.) soft_touch mo_coef call save_mos diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 7d194a54..05f07b22 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -449,7 +449,7 @@ subroutine save_osoci_natural_mos label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.) !if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Write" )then ! disk_access_ao_integrals = "Read" ! touch disk_access_ao_integrals @@ -584,7 +584,7 @@ subroutine set_osoci_natural_mos enddo label = "Natural" - call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1) + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.) soft_touch mo_coef deallocate(tmp,occ) diff --git a/plugins/Hartree_Fock/DIIS.irp.f b/plugins/Hartree_Fock/DIIS.irp.f index 614a9173..6ebb5879 100644 --- a/plugins/Hartree_Fock/DIIS.irp.f +++ b/plugins/Hartree_Fock/DIIS.irp.f @@ -94,7 +94,7 @@ END_PROVIDER do i=1,AO_num do j=1,AO_num - Xt(i,j) = X_Matrix_AO(j,i) + Xt(i,j) = S_half_inv(j,i) enddo enddo @@ -103,7 +103,7 @@ END_PROVIDER call dgemm('N','N',AO_num,AO_num,AO_num, & 1.d0, & Fock_matrix_AO,size(Fock_matrix_AO,1), & - X_Matrix_AO,size(X_Matrix_AO,1), & + S_half_inv,size(S_half_inv,1), & 0.d0, & eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1)) @@ -130,67 +130,10 @@ END_PROVIDER call dgemm('N','N',AO_num,AO_num,AO_num, & 1.d0, & - X_matrix_AO,size(X_matrix_AO,1), & + S_half_inv,size(S_half_inv,1), & scratch,size(scratch,1), & 0.d0, & eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1)) END_PROVIDER -BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num,AO_num) ] - - BEGIN_DOC -! Matrix X = S^{-1/2} obtained by SVD - END_DOC - - implicit none - - integer :: num_linear_dependencies - integer :: LDA, LDC - double precision, allocatable :: U(:,:),Vt(:,:), D(:) - integer :: info, i, j, k - - LDA = size(AO_overlap,1) - LDC = size(X_matrix_AO,1) - - allocate( & - U(LDC,AO_num), & - Vt(LDA,AO_num), & - D(AO_num)) - - call svd( & - AO_overlap,LDA, & - U,LDC, & - D, & - Vt,LDA, & - AO_num,AO_num) - - num_linear_dependencies = 0 - do i=1,AO_num - print*,D(i) - if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then - D(i) = 0.d0 - num_linear_dependencies += 1 - else - ASSERT (D(i) > 0.d0) - D(i) = 1.d0/sqrt(D(i)) - endif - do j=1,AO_num - X_matrix_AO(j,i) = 0.d0 - enddo - enddo - write(*,*) 'linear dependencies',num_linear_dependencies -! stop - - do k=1,AO_num - if(D(k) /= 0.d0) then - do j=1,AO_num - do i=1,AO_num - X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j) - enddo - enddo - endif - enddo - - -END_PROVIDER diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 087933c8..0764c83f 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -263,17 +263,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num) BEGIN_DOC ! Fock matrix on the MO basis END_DOC - double precision, allocatable :: T(:,:) - allocate ( T(ao_num,mo_tot_num) ) - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & - mo_coef, size(mo_coef,1), & - 0.d0, T, size(T,1)) - call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & - 1.d0, mo_coef,size(mo_coef,1), & - T, size(T,1), & - 0.d0, Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1)) - deallocate(T) + call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & + Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1)) END_PROVIDER @@ -282,17 +273,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num) BEGIN_DOC ! Fock matrix on the MO basis END_DOC - double precision, allocatable :: T(:,:) - allocate ( T(ao_num,mo_tot_num) ) - call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & - mo_coef, size(mo_coef,1), & - 0.d0, T, size(T,1)) - call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & - 1.d0, mo_coef,size(mo_coef,1), & - T, size(T,1), & - 0.d0, Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1)) - deallocate(T) + call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & + Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1)) END_PROVIDER BEGIN_PROVIDER [ double precision, HF_energy ] @@ -330,97 +312,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] enddo enddo else - double precision, allocatable :: T(:,:), M(:,:) - integer :: ierr - ! F_ao = S C F_mo C^t S - allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr) - if (ierr /=0 ) then - print *, irp_here, ' : allocation failed' - endif - -! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) -! -> M(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num) -! -> T(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - Fock_matrix_mo, size(Fock_matrix_mo,1), & - 0.d0, & - T, size(T,1)) - -! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) -! -> M(ao_num,ao_num) - call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) -! -> Fock_matrix_ao(ao_num,ao_num) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - Fock_matrix_ao, size(Fock_matrix_ao,1)) - - - deallocate(T) + call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), & + Fock_matrix_ao,size(Fock_matrix_ao,1)) endif END_PROVIDER -subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) - implicit none - integer, intent(in) :: LDFMO ! size(FMO,1) - integer, intent(in) :: LDFAO ! size(FAO,1) - double precision, intent(in) :: FMO(LDFMO,*) - double precision, intent(out) :: FAO(LDFAO,*) - - double precision, allocatable :: T(:,:), M(:,:) - integer :: ierr - ! F_ao = S C F_mo C^t S - allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr) - if (ierr /=0 ) then - print *, irp_here, ' : allocation failed' - endif - -! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) -! -> M(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num) -! -> T(ao_num,mo_tot_num) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - FMO, size(FMO,1), & - 0.d0, & - T, size(T,1)) - -! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) -! -> M(ao_num,ao_num) - call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - -! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) -! -> Fock_matrix_ao(ao_num,ao_num) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - FAO, size(FAO,1)) - deallocate(T,M) -end diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f index bf90d68d..241721ae 100644 --- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f +++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f @@ -83,7 +83,6 @@ END_DOC ! Calculate error vectors max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO)) -! max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_AO)) ! SCF energy @@ -134,7 +133,7 @@ END_DOC write(output_hartree_fock,*) if(.not.no_oa_or_av_opt)then - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) endif call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy') diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f index 5336e8e0..3d71d826 100644 --- a/plugins/Hartree_Fock/SCF.irp.f +++ b/plugins/Hartree_Fock/SCF.irp.f @@ -23,7 +23,7 @@ subroutine create_guess mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef mo_label = 'Guess' - call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label) + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.) SOFT_TOUCH mo_coef mo_label else if (mo_guess_type == "Huckel") then call huckel_guess diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index 2be14ed3..20a8abd7 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -119,7 +119,7 @@ subroutine damping_SCF write(output_hartree_fock,*) if(.not.no_oa_or_av_opt)then - call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.) endif call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f index 9cb84507..ad8c16f7 100644 --- a/plugins/Hartree_Fock/huckel.irp.f +++ b/plugins/Hartree_Fock/huckel.irp.f @@ -11,21 +11,24 @@ subroutine huckel_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label,1) + size(mo_mono_elec_integral,2),label,1,.false.) TOUCH mo_coef c = 0.5d0 * 1.75d0 do j=1,ao_num do i=1,ao_num - Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + & - ao_mono_elec_integral_diag(j)) + Fock_matrix_ao_alpha(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + & + ao_mono_elec_integral_diag(j)) + Fock_matrix_ao_beta (i,j) = Fock_matrix_ao_alpha(i,j) enddo - Fock_matrix_ao(j,j) = Fock_matrix_ao_alpha(j,j) + Fock_matrix_ao_alpha(j,j) = ao_mono_elec_integral(j,j) + ao_bi_elec_integral_alpha(j,j) + Fock_matrix_ao_beta (j,j) = Fock_matrix_ao_alpha(j,j) enddo - TOUCH Fock_matrix_ao + TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta mo_coef = eigenvectors_fock_matrix_mo SOFT_TOUCH mo_coef call save_mos + print *, 'E=', HF_energy end diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index cd2c8e7d..1d3b24e6 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -1,7 +1,7 @@ - BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num,ao_num) ] + BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num_align,ao_num) ] implicit none BEGIN_DOC ! Overlap between atomic basis functions: @@ -69,7 +69,7 @@ END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ] +BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ] implicit none BEGIN_DOC ! Overlap between absolute value of atomic basis functions: @@ -123,11 +123,68 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ] !$OMP END PARALLEL DO END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_overlap_inv,(ao_num,ao_num) ] +BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ] implicit none BEGIN_DOC ! S^-1 END_DOC - call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,ao_overlap_inv,size(ao_overlap_inv,1)) + call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1)) END_PROVIDER +BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] + + BEGIN_DOC +! Matrix X = S^{-1/2} obtained by SVD + END_DOC + + implicit none + + integer :: num_linear_dependencies + integer :: LDA, LDC + double precision, allocatable :: U(:,:),Vt(:,:), D(:) + integer :: info, i, j, k + + LDA = size(AO_overlap,1) + LDC = size(S_half_inv,1) + + allocate( & + U(LDC,AO_num), & + Vt(LDA,AO_num), & + D(AO_num)) + + call svd( & + AO_overlap,LDA, & + U,LDC, & + D, & + Vt,LDA, & + AO_num,AO_num) + + num_linear_dependencies = 0 + do i=1,AO_num + print*,D(i) + if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0/sqrt(D(i)) + endif + do j=1,AO_num + S_half_inv(j,i) = 0.d0 + enddo + enddo + write(*,*) 'linear dependencies',num_linear_dependencies +! stop + + do k=1,AO_num + if(D(k) /= 0.d0) then + do j=1,AO_num + do i=1,AO_num + S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + endif + enddo + + +END_PROVIDER diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index e8c2e7b1..2d3da314 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -337,8 +337,8 @@ end END_DOC mo_integrals_cache_min_8 = max(1_8,elec_alpha_num - 63_8) mo_integrals_cache_max_8 = min(int(mo_tot_num,8),mo_integrals_cache_min+127_8) - mo_integrals_cache_min = mo_integrals_cache_min_8 - mo_integrals_cache_max = mo_integrals_cache_max_8 + mo_integrals_cache_min = max(1,elec_alpha_num - 63) + mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+127) END_PROVIDER @@ -349,17 +349,22 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*12 END_DOC PROVIDE mo_bielec_integrals_in_map integer*8 :: i,j,k,l + integer*8 :: i4,j4,k4,l4 integer*8 :: ii integer(key_kind) :: idx real(integral_kind) :: integral FREE ao_integrals_cache - !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) + !$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral) do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8 do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8 do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8 do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8 + i4 = int(i,4) + j4 = int(j,4) + k4 = int(k,4) + l4 = int(l,4) !DIR$ FORCEINLINE - call bielec_integrals_index(i,j,k,l,idx) + call bielec_integrals_index(i4,j4,k4,l4,idx) !DIR$ FORCEINLINE call map_get(mo_integrals_map,idx,integral) ii = l-mo_integrals_cache_min_8 diff --git a/src/Integrals_Monoelec/ao_mono_ints.irp.f b/src/Integrals_Monoelec/ao_mono_ints.irp.f index 87d03ac4..8cfd25cf 100644 --- a/src/Integrals_Monoelec/ao_mono_ints.irp.f +++ b/src/Integrals_Monoelec/ao_mono_ints.irp.f @@ -1,4 +1,4 @@ - BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)] &BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_diag,(ao_num)] implicit none integer :: i,j,n,l diff --git a/src/Integrals_Monoelec/kin_ao_ints.irp.f b/src/Integrals_Monoelec/kin_ao_ints.irp.f index 8c8859fe..da4de4d7 100644 --- a/src/Integrals_Monoelec/kin_ao_ints.irp.f +++ b/src/Integrals_Monoelec/kin_ao_ints.irp.f @@ -1,6 +1,6 @@ - BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num,ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num,ao_num) ] + BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num_align,ao_num) ] implicit none integer :: i,j,n,l double precision :: f @@ -119,7 +119,7 @@ END_PROVIDER -BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num,ao_num)] +BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)] implicit none BEGIN_DOC ! array of the priminitve basis kinetic integrals @@ -134,11 +134,14 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num,ao_num)] else !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j) & - !$OMP SHARED(ao_num, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z) + !$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z) do j = 1, ao_num do i = 1, ao_num ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) ) enddo + do i = ao_num +1,ao_num_align + ao_kinetic_integral(i,j) = 0.d0 + enddo enddo !$OMP END PARALLEL DO endif diff --git a/src/Integrals_Monoelec/kin_mo_ints.irp.f b/src/Integrals_Monoelec/kin_mo_ints.irp.f index 262e4805..2301c23d 100644 --- a/src/Integrals_Monoelec/kin_mo_ints.irp.f +++ b/src/Integrals_Monoelec/kin_mo_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num,mo_tot_num)] +BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num_align,mo_tot_num)] implicit none BEGIN_DOC ! Kinetic energy integrals in the MO basis diff --git a/src/Integrals_Monoelec/mo_mono_ints.irp.f b/src/Integrals_Monoelec/mo_mono_ints.irp.f index 891ed3d5..0d912852 100644 --- a/src/Integrals_Monoelec/mo_mono_ints.irp.f +++ b/src/Integrals_Monoelec/mo_mono_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num,mo_tot_num)] +BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_tot_num)] implicit none integer :: i,j,n,l BEGIN_DOC diff --git a/src/Integrals_Monoelec/pot_ao_ints.irp.f b/src/Integrals_Monoelec/pot_ao_ints.irp.f index bc1884cd..7116d2c7 100644 --- a/src/Integrals_Monoelec/pot_ao_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num,ao_num)] +BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)] BEGIN_DOC ! interaction nuclear electron END_DOC @@ -80,7 +80,7 @@ BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num,ao_num)] END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num,ao_num,nucl_num)] + BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num_align,ao_num,nucl_num)] BEGIN_DOC ! ao_nucl_elec_integral_per_atom(i,j,k) = - ! where Rk is the geometry of the kth atom diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index 2a1eaf67..22cceab9 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num,ao_num)] +BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)] implicit none BEGIN_DOC ! Pseudo-potential integrals @@ -29,7 +29,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num,ao_num)] END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num,ao_num)] +BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)] implicit none BEGIN_DOC ! Local pseudo-potential @@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num,ao_num)] END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)] implicit none BEGIN_DOC ! Local pseudo-potential diff --git a/src/Integrals_Monoelec/pot_mo_ints.irp.f b/src/Integrals_Monoelec/pot_mo_ints.irp.f index 5810b4f3..7c7e306f 100644 --- a/src/Integrals_Monoelec/pot_mo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_mo_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num,mo_tot_num)] +BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)] implicit none BEGIN_DOC ! interaction nuclear electron on the MO basis @@ -25,7 +25,7 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num,mo_tot_num) END_PROVIDER -BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num,mo_tot_num,nucl_num)] +BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num_align,mo_tot_num,nucl_num)] implicit none BEGIN_DOC ! mo_nucl_elec_integral_per_atom(i,j,k) = - diff --git a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f index 47e6e277..f2fee5f4 100644 --- a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num,mo_tot_num)] +BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_num)] implicit none BEGIN_DOC ! interaction nuclear electron on the MO basis diff --git a/src/Integrals_Monoelec/spread_dipole_ao.irp.f b/src/Integrals_Monoelec/spread_dipole_ao.irp.f index 2ff1494f..c510d58b 100644 --- a/src/Integrals_Monoelec/spread_dipole_ao.irp.f +++ b/src/Integrals_Monoelec/spread_dipole_ao.irp.f @@ -1,6 +1,6 @@ - BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num_align,ao_num)] BEGIN_DOC ! array of the integrals of AO_i * x^2 AO_j ! array of the integrals of AO_i * y^2 AO_j @@ -69,9 +69,9 @@ - BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num_align,ao_num)] BEGIN_DOC ! array of the integrals of AO_i * x AO_j ! array of the integrals of AO_i * y AO_j @@ -139,9 +139,9 @@ !$OMP END PARALLEL DO END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num,ao_num)] - &BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num,ao_num)] + BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num_align,ao_num)] BEGIN_DOC ! array of the integrals of AO_i * d/dx AO_j ! array of the integrals of AO_i * d/dy AO_j diff --git a/src/Integrals_Monoelec/spread_dipole_mo.irp.f b/src/Integrals_Monoelec/spread_dipole_mo.irp.f index 9e21ec21..aa5ef8aa 100644 --- a/src/Integrals_Monoelec/spread_dipole_mo.irp.f +++ b/src/Integrals_Monoelec/spread_dipole_mo.irp.f @@ -1,6 +1,6 @@ - BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num,mo_tot_num)] -&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num,mo_tot_num)] -&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num,mo_tot_num)] + BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num_align,mo_tot_num)] BEGIN_DOC ! array of the integrals of MO_i * x MO_j ! array of the integrals of MO_i * y MO_j @@ -29,9 +29,9 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num,mo_tot_num)] -&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num,mo_tot_num)] -&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num,mo_tot_num)] + BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num_align,mo_tot_num)] BEGIN_DOC ! array of the integrals of MO_i * x^2 MO_j ! array of the integrals of MO_i * y^2 MO_j diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index d3e2eef9..d771feec 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -5,11 +5,5 @@ program H_CORE_guess END_DOC implicit none character*(64) :: label - label = "Guess" - call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & - size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label,1) - print *, 'save mos' - call save_mos - + call h_core_guess end diff --git a/src/MOGuess/h_core_guess_routine.irp.f b/src/MOGuess/h_core_guess_routine.irp.f index 23899160..5b4ede91 100644 --- a/src/MOGuess/h_core_guess_routine.irp.f +++ b/src/MOGuess/h_core_guess_routine.irp.f @@ -7,8 +7,7 @@ subroutine hcore_guess label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & - size(mo_mono_elec_integral,2),label,1) - print *, 'save mos' + size(mo_mono_elec_integral,2),label,1,.false.) call save_mos SOFT_TOUCH mo_coef mo_label end diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 7bcd7bd2..c2cbc5f2 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -26,7 +26,7 @@ BEGIN_PROVIDER [ integer, mo_tot_num_align ] END_PROVIDER -BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ] + BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ] implicit none BEGIN_DOC ! Molecular orbital coefficients on AO basis set @@ -35,7 +35,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ] END_DOC integer :: i, j double precision, allocatable :: buffer(:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer logical :: exists PROVIDE ezfio_filename @@ -50,9 +49,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ] do j=1,ao_num mo_coef(j,i) = buffer(j,i) enddo - do j=ao_num+1,ao_num_align - mo_coef(j,i) = 0.d0 - enddo enddo deallocate(buffer) else @@ -61,9 +57,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ] do j=1,ao_num mo_coef(j,i) = ao_ortho_canonical_coef(j,i) enddo - do j=ao_num+1,ao_num_align - mo_coef(j,i) = 0.d0 - enddo enddo endif END_PROVIDER @@ -125,7 +118,7 @@ BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num_align, mo_tot_num) ] END_DOC call dgemm('N','N', ao_num, mo_tot_num, ao_num, & - 1.d0, ao_overlap,size(ao_overlap,1), & + 1.d0, ao_overlap,size(ao_overlap,1), & mo_coef, size(mo_coef,1), & 0.d0, S_mo_coef, size(S_mo_coef,1)) @@ -159,24 +152,25 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) BEGIN_DOC ! Transform A from the AO basis to the MO basis ! - ! (S.C)t.A_ao.S.C + ! Ct.A_ao.C END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_ao(LDA_ao,ao_num) double precision, intent(out) :: A_mo(LDA_mo,mo_tot_num) double precision, allocatable :: T(:,:) - allocate ( T(mo_tot_num,ao_num) ) + allocate ( T(ao_num_align,mo_tot_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - call dgemm('T','N', mo_tot_num, ao_num, ao_num, & - 1.d0, S_mo_coef, size(S_mo_coef,1), & - A_ao, size(A_ao,1), & - 0.d0, T, size(T,1)) + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, A_ao,LDA_ao, & + mo_coef, size(mo_coef,1), & + 0.d0, T, ao_num_align) - call dgemm('N','N', mo_tot_num, mo_tot_num, ao_num, & - 1.d0, T,size(T,1), & - S_mo_coef, size(S_mo_coef,1), & - 0.d0, A_mo, size(A_mo,1)) + call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & + 1.d0, mo_coef,size(mo_coef,1), & + T, ao_num_align, & + 0.d0, A_mo, LDA_mo) deallocate(T) end @@ -186,7 +180,7 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) BEGIN_DOC ! Transform A from the MO basis to the AO basis ! - ! C.A_mo.Ct + ! (S.C).A_mo.(S.C)t END_DOC integer, intent(in) :: LDA_ao,LDA_mo double precision, intent(in) :: A_mo(LDA_mo,mo_tot_num) @@ -194,15 +188,14 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao) double precision, allocatable :: T(:,:) allocate ( T(mo_tot_num_align,ao_num) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & - 1.d0, A_mo,LDA_mo, & - mo_coef, size(mo_coef,1), & + 1.d0, A_mo,size(A_mo,1), & + S_mo_coef, size(S_mo_coef,1), & 0.d0, T, size(T,1)) call dgemm('N','N', ao_num, ao_num, mo_tot_num, & - 1.d0, mo_coef,size(mo_coef,1), & + 1.d0, S_mo_coef, size(S_mo_coef,1), & T, size(T,1), & 0.d0, A_ao, size(A_ao,1)) @@ -257,18 +250,17 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) double precision, allocatable :: T(:,:) 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('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) + ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),& + A_ao,size(A_ao,1), & + 0.d0, T, size(T,1)) - call dgemm('N','N', ao_num, ao_num, ao_num, 1.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_inv,size(ao_ortho_canonical_coef_inv,1),& - 0.d0, A, LDA) + 0.d0, A, size(A,1)) deallocate(T) end diff --git a/src/MO_Basis/utils.irp.f b/src/MO_Basis/utils.irp.f index 750e3420..4221fce4 100644 --- a/src/MO_Basis/utils.irp.f +++ b/src/MO_Basis/utils.irp.f @@ -44,11 +44,12 @@ subroutine save_mos_truncated(n) end -subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) +subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output) implicit none integer,intent(in) :: n,m, sign character*(64), intent(in) :: label double precision, intent(in) :: matrix(n,m) + logical, intent(in) :: output integer :: i,j double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:) @@ -76,22 +77,26 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign) mo_coef_new = mo_coef call lapack_diag(eigvalues,R,A,n,m) - write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**' - write (output_mo_basis,'(A)') '' - write (output_mo_basis,'(A)') 'Eigenvalues' - write (output_mo_basis,'(A)') '-----------' - write (output_mo_basis,'(A)') '' - write (output_mo_basis,'(A)') '======== ================' + if (output) then + write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**' + write (output_mo_basis,'(A)') '' + write (output_mo_basis,'(A)') 'Eigenvalues' + write (output_mo_basis,'(A)') '-----------' + write (output_mo_basis,'(A)') '' + write (output_mo_basis,'(A)') '======== ================' + endif if (sign == -1) then do i=1,m eigvalues(i) = -eigvalues(i) enddo endif - do i=1,m - write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i) - enddo - write (output_mo_basis,'(A)') '======== ================' - write (output_mo_basis,'(A)') '' + if (output) then + do i=1,m + write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i) + enddo + write (output_mo_basis,'(A)') '======== ================' + write (output_mo_basis,'(A)') '' + endif call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1)) deallocate(A,mo_coef_new,R,eigvalues)