diff --git a/src/ao_cart_basis/EZFIO.cfg b/src/ao_cart_basis/EZFIO.cfg index 6a1fb099..dea3fc48 100644 --- a/src/ao_cart_basis/EZFIO.cfg +++ b/src/ao_cart_basis/EZFIO.cfg @@ -44,6 +44,12 @@ doc: Exponents for each primitive of each |AO| size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_prim_num_max) interface: ezfio, provider +[use_cgtos] +type: logical +doc: If true, use cgtos for AO integrals +interface: ezfio +default: False + [ao_cart_expo_im] type: double precision doc: imag part for Exponents for each primitive of each cGTOs |AO| diff --git a/src/ao_cart_two_e_ints/cholesky.irp.f b/src/ao_cart_two_e_ints/cholesky.irp.f index eb0dc0a9..8f11ad9b 100644 --- a/src/ao_cart_two_e_ints/cholesky.irp.f +++ b/src/ao_cart_two_e_ints/cholesky.irp.f @@ -457,7 +457,7 @@ END_PROVIDER write(iunit) rank write(iunit) cholesky_ao_cart$_erf close(iunit) - call ezfio_set_ao_cart_two_e_ints_io_ao_cart$_erf_cholesky('Read') + call ezfio_set_two_e_ints_keywords_io_ao_cart$_erf_cholesky('Read') endif endif diff --git a/src/ao_two_e_ints/ao_bi_integrals.irp.f b/src/ao_two_e_ints/ao_bi_integrals.irp.f index 296f6ef9..e69de29b 100644 --- a/src/ao_two_e_ints/ao_bi_integrals.irp.f +++ b/src/ao_two_e_ints/ao_bi_integrals.irp.f @@ -1,150 +0,0 @@ -subroutine ao_two_e_integrals_index(i,j,k,l,i1) - use map_module - implicit none - BEGIN_DOC - ! Computes an unique index for i,j,k,l integrals - END_DOC - integer, intent(in) :: i,j,k,l - integer(key_kind), intent(out) :: i1 - integer(key_kind) :: p,q,r,s,i2 - p = min(i,k) - r = max(i,k) - p = p+shiftr(r*r-r,1) - q = min(j,l) - s = max(j,l) - q = q+shiftr(s*s-s,1) - i1 = min(p,q) - i2 = max(p,q) - i1 = i1+shiftr(i2*i2-i2,1) -end - -BEGIN_TEMPLATE - -BEGIN_PROVIDER [ logical, ao_two_e_integrals$_erf_in_map ] - use map_module - implicit none - BEGIN_DOC - ! If True, the map of AO two-electron integrals$_erf is provided - END_DOC - integer(bit_kind) :: mask_ijkl(N_int,4) - integer(bit_kind) :: mask_ijk(N_int,3) - double precision :: cpu_1, cpu_2, wall_1, wall_2 - - - ao_two_e_integrals$_erf_in_map = .True. - if (read_ao_two_e_integrals$_erf) then - print*,'Reading the AO integrals$_erf' - call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints$_erf',ao_integrals$_erf_map) - print*, 'AO integrals$_erf provided' - return - endif - - call wall_time(wall_1) - call cpu_time(cpu_1) - - if (do_ao_cholesky) then - PROVIDE cholesky_ao_transp - else - call add_integrals$_erf_to_map_cholesky_ao - integer*8 :: get_ao_map_size, ao_map_size - ao_map_size = get_ao_map_size() - double precision, external :: map_mb - print*,'Atomic integrals$_erf provided:' - print*,' Size of AO map ', map_mb(ao_integrals$_erf_map) ,'MB' - print*,' Number of AO integrals$_erf: ', ao_map_size - endif - - call wall_time(wall_2) - call cpu_time(cpu_2) - - print*,' cpu time :',cpu_2 - cpu_1, 's' - print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' - - if (write_ao_two_e_integrals$_erf.and.mpi_master) then - call ezfio_set_work_empty(.False.) - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints$_erf',ao_integrals$_erf_map) - call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals$_erf('Read') - endif - -END_PROVIDER - -subroutine add_integrals$_erf_to_map_cholesky_ao - use bitmasks - implicit none - - BEGIN_DOC - ! Adds integrals$_erf to the AO map using Cholesky vectors - END_DOC - - integer :: i,j,k,l,m - integer :: size_buffer, n_integrals$_erf - size_buffer = min(ao_num*ao_num*ao_num,16000000) - - double precision, allocatable :: Vtmp(:,:,:) - integer(key_kind) , allocatable :: buffer_i(:) - real(integral_kind), allocatable :: buffer_value(:) - - PROVIDE cholesky_ao_transp - call set_multiple_levels_omp(.False.) - - !$OMP PARALLEL DEFAULT(SHARED) & - !$OMP PRIVATE(i,j,k,l,n_integrals$_erf,buffer_value, buffer_i, Vtmp) - allocate (buffer_i(size_buffer), buffer_value(size_buffer)) - allocate (Vtmp(ao_num,ao_num,ao_num)) - n_integrals$_erf = 0 - - !$OMP DO SCHEDULE(dynamic) - do l=1,ao_num - call dgemm('T','N',ao_num*ao_num,ao_num,cholesky_ao_num,1.d0, & - cholesky_ao_transp, cholesky_ao_num, & - cholesky_ao_transp(1,1,l), cholesky_ao_num, 0.d0, & - Vtmp, ao_num*ao_num) - - do k=1,l - do j=1,ao_num - do i=1,j - if (dabs(Vtmp(i,j,k)) > ao_integrals_threshold) then - n_integrals$_erf = n_integrals$_erf + 1 - buffer_value(n_integrals$_erf) = Vtmp(i,j,k) - !DIR$ FORCEINLINE - call ao_two_e_integrals$_erf_index(i,k,j,l,buffer_i(n_integrals$_erf)) - if (n_integrals$_erf == size_buffer) then - call map_append(ao_integrals$_erf_map, buffer_i, buffer_value, n_integrals$_erf) - n_integrals$_erf = 0 - endif - endif - enddo - enddo - enddo - enddo - !$OMP END DO NOWAIT - - if (n_integrals$_erf > 0) then - call map_append(ao_integrals$_erf_map, buffer_i, buffer_value, n_integrals$_erf) - endif - deallocate(buffer_i, buffer_value, Vtmp) - !$OMP BARRIER - !$OMP END PARALLEL - - call map_sort(ao_integrals$_erf_map) - call map_unique(ao_integrals$_erf_map) - -end - -subroutine clear_ao$_erf_map - implicit none - BEGIN_DOC - ! Frees the meaory of the AO map - END_DOC - call map_deinit(ao_integrals$_erf_map) - FREE ao_integrals$_erf_map - FREE ao_two_e_integrals$_erf_in_map -end - -SUBST [ _erf ] - -;; -_erf;; -_cgtos;; - -END_TEMPLATE diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index e1e9352e..9153fa96 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -1,9 +1,79 @@ use map_module +subroutine ao_two_e_integrals_index(i,j,k,l,i1) + use map_module + implicit none + BEGIN_DOC + ! Computes an unique index for i,j,k,l integrals + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind), intent(out) :: i1 + integer(key_kind) :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+shiftr(r*r-r,1) + q = min(j,l) + s = max(j,l) + q = q+shiftr(s*s-s,1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+shiftr(i2*i2-i2,1) +end + + !! AO Map !! ====== BEGIN_TEMPLATE +BEGIN_PROVIDER [ logical, ao_two_e_integrals$_erf_in_map ] + use map_module + implicit none + BEGIN_DOC + ! If True, the map of AO two-electron integrals$_erf is provided + END_DOC + integer(bit_kind) :: mask_ijkl(N_int,4) + integer(bit_kind) :: mask_ijk(N_int,3) + double precision :: cpu_1, cpu_2, wall_1, wall_2 + + + ao_two_e_integrals$_erf_in_map = .True. + if (read_ao_two_e_integrals$_erf) then + print*,'Reading the AO integrals$_erf' + call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints$_erf',ao_integrals$_erf_map) + print*, 'AO integrals$_erf provided' + return + endif + + call wall_time(wall_1) + call cpu_time(cpu_1) + + if (do_ao_cholesky) then + PROVIDE cholesky_ao_transp + else + call add_integrals$_erf_to_map_cholesky_ao + integer*8 :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + double precision, external :: map_mb + print*,'Atomic integrals$_erf provided:' + print*,' Size of AO map ', map_mb(ao_integrals$_erf_map) ,'MB' + print*,' Number of AO integrals$_erf: ', ao_map_size + endif + + call wall_time(wall_2) + call cpu_time(cpu_2) + + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + + if (write_ao_two_e_integrals$_erf.and.mpi_master) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints$_erf',ao_integrals$_erf_map) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals$_erf('Read') + endif + +END_PROVIDER + + BEGIN_PROVIDER [ type(map_type), ao_integrals$_erf_map ] implicit none BEGIN_DOC @@ -279,16 +349,6 @@ function get_ao$_erf_map_size() get_ao$_erf_map_size = ao_integrals$_erf_map % n_elements end -subroutine clear_ao$_erf_map - implicit none - BEGIN_DOC - ! Frees the memory of the AO map - END_DOC - call map_deinit(ao_integrals$_erf_map) - FREE ao_integrals$_erf_map -end - - subroutine insert_into_ao_integrals$_erf_map(n_integrals,buffer_i, buffer_values) use map_module implicit none @@ -379,6 +439,80 @@ integer function load_ao_integrals$_erf(filename) end +subroutine add_integrals$_erf_to_map_cholesky_ao + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals$_erf to the AO map using Cholesky vectors + END_DOC + + integer :: i,j,k,l,m + integer :: size_buffer, n_integrals$_erf + size_buffer = min(ao_num*ao_num*ao_num,16000000) + + double precision, allocatable :: Vtmp(:,:,:) + integer(key_kind) , allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_value(:) + + PROVIDE cholesky_ao_transp + call set_multiple_levels_omp(.False.) + + !$OMP PARALLEL DEFAULT(SHARED) & + !$OMP PRIVATE(i,j,k,l,n_integrals$_erf,buffer_value, buffer_i, Vtmp) + allocate (buffer_i(size_buffer), buffer_value(size_buffer)) + allocate (Vtmp(ao_num,ao_num,ao_num)) + n_integrals$_erf = 0 + + !$OMP DO SCHEDULE(dynamic) + do l=1,ao_num + call dgemm('T','N',ao_num*ao_num,ao_num,cholesky_ao_num,1.d0, & + cholesky_ao_transp, cholesky_ao_num, & + cholesky_ao_transp(1,1,l), cholesky_ao_num, 0.d0, & + Vtmp, ao_num*ao_num) + + do k=1,l + do j=1,ao_num + do i=1,j + if (dabs(Vtmp(i,j,k)) > ao_integrals_threshold) then + n_integrals$_erf = n_integrals$_erf + 1 + buffer_value(n_integrals$_erf) = Vtmp(i,j,k) + !DIR$ FORCEINLINE + call ao_two_e_integrals_index(i,k,j,l,buffer_i(n_integrals$_erf)) + if (n_integrals$_erf == size_buffer) then + call map_append(ao_integrals$_erf_map, buffer_i, buffer_value, n_integrals$_erf) + n_integrals$_erf = 0 + endif + endif + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + if (n_integrals$_erf > 0) then + call map_append(ao_integrals$_erf_map, buffer_i, buffer_value, n_integrals$_erf) + endif + deallocate(buffer_i, buffer_value, Vtmp) + !$OMP BARRIER + !$OMP END PARALLEL + + call map_sort(ao_integrals$_erf_map) + call map_unique(ao_integrals$_erf_map) + +end + + + +subroutine clear_ao$_erf_map + implicit none + BEGIN_DOC + ! Frees the meaory of the AO map + END_DOC + call map_deinit(ao_integrals$_erf_map) + FREE ao_integrals$_erf_map + FREE ao_two_e_integrals$_erf_in_map +end SUBST [ _erf ] diff --git a/src/ao_two_e_ints/schwartz.irp.f b/src/ao_two_e_ints/schwartz.irp.f new file mode 100644 index 00000000..01ca82aa --- /dev/null +++ b/src/ao_two_e_ints/schwartz.irp.f @@ -0,0 +1,31 @@ +BEGIN_TEMPLATE + +BEGIN_PROVIDER [ double precision , ao_two_e_integral$_erf_schwartz, (ao_num, ao_num) ] + implicit none + integer :: i,j + double precision :: get_ao_two_e_integral$_erf + double precision :: get_ao$_erf_integ_chol + if(do_ao_cholesky)then + do i = 1, ao_num + do j = 1, ao_num + ao_two_e_integral$_erf_schwartz(j,i) = get_ao$_erf_integ_chol(i,j,i,i) + enddo + enddo + else + PROVIDE ao_two_e_integrals$_erf_in_map + do i = 1, ao_num + do j = 1, ao_num + ao_two_e_integral$_erf_schwartz(j,i) = get_ao_two_e_integral$_erf(i, i, j, j) + enddo + enddo + endif + +END_PROVIDER + +SUBST [ _erf ] + +;; +_erf;; +_cgtos;; + +END_TEMPLATE diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f new file mode 100644 index 00000000..9c74b485 --- /dev/null +++ b/src/ao_two_e_ints/screening.irp.f @@ -0,0 +1,16 @@ +logical function ao_two_e_integral_zero(i,j,k,l) + implicit none + integer, intent(in) :: i,j,k,l + + ao_two_e_integral_zero = .False. +! if (.not.(read_ao_two_e_integrals.or.is_periodic.or.use_cgtos)) then + if (.not.(is_periodic.or.use_cgtos)) then + if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then + ao_two_e_integral_zero = .True. + return + endif + if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then + ao_two_e_integral_zero = .True. + endif + endif +end diff --git a/src/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f index 6d917322..371f2cfd 100644 --- a/src/hartree_fock/fock_matrix_hf.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -25,81 +25,83 @@ ao_two_e_integral_beta = 0.d0 if (do_direct_integrals) then - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0,& - !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2,& - !$OMP local_threshold) & - !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& - !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz,& - !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta) - - allocate(keys(1), values(1)) - allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & - ao_two_e_integral_beta_tmp(ao_num,ao_num)) - ao_two_e_integral_alpha_tmp = 0.d0 - ao_two_e_integral_beta_tmp = 0.d0 - - q = ao_num*ao_num*ao_num*ao_num - !$OMP DO SCHEDULE(static,64) - do p=1_8,q - call two_e_integrals_index_reverse(kk,ii,ll,jj,p) - if ( (kk(1)>ao_num).or. & - (ii(1)>ao_num).or. & - (jj(1)>ao_num).or. & - (ll(1)>ao_num) ) then - cycle - endif - k = kk(1) - i = ii(1) - l = ll(1) - j = jj(1) - - logical, external :: ao_two_e_integral_zero - if (ao_two_e_integral_zero(i,k,j,l)) then - cycle - endif - local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) - if (local_threshold < ao_integrals_threshold) then - cycle - endif - i0 = i - j0 = j - k0 = k - l0 = l - values(1) = 0.d0 - local_threshold = ao_integrals_threshold/local_threshold - do k2=1,8 - if (kk(k2)==0) then - cycle - endif - i = ii(k2) - j = jj(k2) - k = kk(k2) - l = ll(k2) - c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) - c1 = SCF_density_matrix_ao_alpha(k,i) - c2 = SCF_density_matrix_ao_beta(k,i) - if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then - cycle - endif - if (values(1) == 0.d0) then - values(1) = ao_two_e_integral(k0,l0,i0,j0) - endif - integral = c0 * values(1) - ao_two_e_integral_alpha_tmp(i,j) += integral - ao_two_e_integral_beta_tmp (i,j) += integral - integral = values(1) - ao_two_e_integral_alpha_tmp(l,j) -= c1 * integral - ao_two_e_integral_beta_tmp (l,j) -= c2 * integral - enddo - enddo - !$OMP END DO NOWAIT - !$OMP CRITICAL - ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp - ao_two_e_integral_beta += ao_two_e_integral_beta_tmp - !$OMP END CRITICAL - deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) - !$OMP END PARALLEL + print*,'not done yet !!' + stop +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0,& +! !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2,& +! !$OMP local_threshold) & +! !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,& +! !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz,& +! !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta) +! +! allocate(keys(1), values(1)) +! allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & +! ao_two_e_integral_beta_tmp(ao_num,ao_num)) +! ao_two_e_integral_alpha_tmp = 0.d0 +! ao_two_e_integral_beta_tmp = 0.d0 +! +! q = ao_num*ao_num*ao_num*ao_num +! !$OMP DO SCHEDULE(static,64) +! do p=1_8,q +! call two_e_integrals_index_reverse(kk,ii,ll,jj,p) +! if ( (kk(1)>ao_num).or. & +! (ii(1)>ao_num).or. & +! (jj(1)>ao_num).or. & +! (ll(1)>ao_num) ) then +! cycle +! endif +! k = kk(1) +! i = ii(1) +! l = ll(1) +! j = jj(1) +! +! logical, external :: ao_two_e_integral_zero +! if (ao_two_e_integral_zero(i,k,j,l)) then +! cycle +! endif +! local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j) +! if (local_threshold < ao_integrals_threshold) then +! cycle +! endif +! i0 = i +! j0 = j +! k0 = k +! l0 = l +! values(1) = 0.d0 +! local_threshold = ao_integrals_threshold/local_threshold +! do k2=1,8 +! if (kk(k2)==0) then +! cycle +! endif +! i = ii(k2) +! j = jj(k2) +! k = kk(k2) +! l = ll(k2) +! c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l) +! c1 = SCF_density_matrix_ao_alpha(k,i) +! c2 = SCF_density_matrix_ao_beta(k,i) +! if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then +! cycle +! endif +! if (values(1) == 0.d0) then +! values(1) = ao_two_e_integral(k0,l0,i0,j0) +! endif +! integral = c0 * values(1) +! ao_two_e_integral_alpha_tmp(i,j) += integral +! ao_two_e_integral_beta_tmp (i,j) += integral +! integral = values(1) +! ao_two_e_integral_alpha_tmp(l,j) -= c1 * integral +! ao_two_e_integral_beta_tmp (l,j) -= c2 * integral +! enddo +! enddo +! !$OMP END DO NOWAIT +! !$OMP CRITICAL +! ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp +! ao_two_e_integral_beta += ao_two_e_integral_beta_tmp +! !$OMP END CRITICAL +! deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) +! !$OMP END PARALLEL else PROVIDE ao_two_e_integrals_in_map diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 4da8ac9b..7a6f1b07 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -347,36 +347,36 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) end -BEGIN_PROVIDER [ double precision, mo_sphe_coef, (ao_sphe_num, mo_num) ] - implicit none - BEGIN_DOC - ! MO coefficients in the basis of spherical harmonics AOs. - END_DOC - double precision, allocatable, dimension (:,:) :: C0, S, tmp - allocate(C0(ao_sphe_num,mo_num), S(mo_num,mo_num), tmp(mo_num,ao_sphe_num)) - - call dgemm('T','N',ao_sphe_num,mo_num,ao_num, 1.d0, & - ao_cart_to_sphe_coef,ao_num, & - mo_coef,size(mo_coef,1), 0.d0, & - C0, size(C0,1)) - - ! C0^T S S^0 - call dgemm('T','N',mo_num,ao_sphe_num,ao_sphe_num, 1.d0, & - C0,size(C0,1), & - ao_sphe_overlap,size(ao_sphe_overlap,1), 0.d0, & - tmp, size(tmp,1)) - - call dgemm('N','N',mo_num,mo_num,ao_sphe_num, 1.d0, & - tmp, size(tmp,1), & - C0,size(C0,1), 0.d0, & - S, size(S,1)) - - integer :: m - m = ao_sphe_num - mo_sphe_coef = C0 - call ortho_lowdin(S,size(S,1),mo_num,mo_sphe_coef,ao_sphe_num,m,1.d-10) - - - deallocate (tmp, S, C0) -END_PROVIDER +!BEGIN_PROVIDER [ double precision, mo_sphe_coef, (ao_sphe_num, mo_num) ] +! implicit none +! BEGIN_DOC +! ! MO coefficients in the basis of spherical harmonics AOs. +! END_DOC +! double precision, allocatable, dimension (:,:) :: C0, S, tmp +! allocate(C0(ao_sphe_num,mo_num), S(mo_num,mo_num), tmp(mo_num,ao_sphe_num)) +! +! call dgemm('T','N',ao_sphe_num,mo_num,ao_num, 1.d0, & +! ao_cart_to_sphe_coef,ao_num, & +! mo_coef,size(mo_coef,1), 0.d0, & +! C0, size(C0,1)) +! +! ! C0^T S S^0 +! call dgemm('T','N',mo_num,ao_sphe_num,ao_sphe_num, 1.d0, & +! C0,size(C0,1), & +! ao_sphe_overlap,size(ao_sphe_overlap,1), 0.d0, & +! tmp, size(tmp,1)) +! +! call dgemm('N','N',mo_num,mo_num,ao_sphe_num, 1.d0, & +! tmp, size(tmp,1), & +! C0,size(C0,1), 0.d0, & +! S, size(S,1)) +! +! integer :: m +! m = ao_sphe_num +! mo_sphe_coef = C0 +! call ortho_lowdin(S,size(S,1),mo_num,mo_sphe_coef,ao_sphe_num,m,1.d-10) +! +! +! deallocate (tmp, S, C0) +!END_PROVIDER diff --git a/src/mo_one_e_ints/ao_to_mo.irp.f b/src/mo_one_e_ints/ao_to_mo.irp.f index 72c1328d..17c55818 100644 --- a/src/mo_one_e_ints/ao_to_mo.irp.f +++ b/src/mo_one_e_ints/ao_to_mo.irp.f @@ -79,81 +79,3 @@ END_PROVIDER ! --- -subroutine mo_to_ao_sphe(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)^\dagger$ - END_DOC - integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_mo(LDA_mo,mo_num) - double precision, intent(out) :: A_ao(LDA_ao,ao_sphe_num) - double precision, allocatable :: T(:,:) - - allocate ( T(mo_num,ao_sphe_num) ) - - call dgemm('N','T', mo_num, ao_sphe_num, mo_num, & - 1.d0, A_mo,size(A_mo,1), & - S_mo_sphe_coef, size(S_mo_sphe_coef,1), & - 0.d0, T, size(T,1)) - - call dgemm('N','N', ao_sphe_num, ao_sphe_num, mo_num, & - 1.d0, S_mo_sphe_coef, size(S_mo_sphe_coef,1), & - T, size(T,1), & - 0.d0, A_ao, size(A_ao,1)) - - deallocate(T) -end - -subroutine mo_to_ao_sphe_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) - implicit none - BEGIN_DOC - ! $C.A_{mo}.C^\dagger$ - END_DOC - integer, intent(in) :: LDA_ao,LDA_mo - double precision, intent(in) :: A_mo(LDA_mo,mo_num) - double precision, intent(out) :: A_ao(LDA_ao,ao_sphe_num) - double precision, allocatable :: T(:,:) - - allocate ( T(mo_num,ao_sphe_num) ) - - call dgemm('N','T', mo_num, ao_sphe_num, mo_num, & - 1.d0, A_mo,size(A_mo,1), & - mo_sphe_coef, size(mo_sphe_coef,1), & - 0.d0, T, size(T,1)) - - call dgemm('N','N', ao_sphe_num, ao_sphe_num, mo_num, & - 1.d0, mo_sphe_coef, size(mo_sphe_coef,1), & - T, size(T,1), & - 0.d0, A_ao, size(A_ao,1)) - - deallocate(T) -end - -BEGIN_PROVIDER [ double precision, S_mo_sphe_coef, (ao_sphe_num, mo_num) ] - implicit none - BEGIN_DOC - ! Product S.C where S is the overlap matrix in the AO basis and C the mo_sphe_coef matrix. - END_DOC - - call dgemm('N','N', ao_sphe_num, mo_num, ao_sphe_num, & - 1.d0, ao_overlap,size(ao_overlap,1), & - mo_sphe_coef, size(mo_coef,1), & - 0.d0, S_mo_coef, size(S_mo_coef,1)) - -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals_from_mo, (ao_sphe_num, ao_sphe_num)] - implicit none - BEGIN_DOC -! Integrals of the one e hamiltonian obtained from the integrals on the MO basis -! -! WARNING : this is equal to ao_one_e_integrals only if the AO and MO basis have the same number of functions - END_DOC - call mo_to_ao_sphe(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_sphe_num) -END_PROVIDER - - - diff --git a/src/mo_one_e_ints/pot_mo_pseudo_ints.irp.f b/src/mo_one_e_ints/pot_mo_pseudo_ints.irp.f index f4c1012d..660e76f5 100644 --- a/src/mo_one_e_ints/pot_mo_pseudo_ints.irp.f +++ b/src/mo_one_e_ints/pot_mo_pseudo_ints.irp.f @@ -25,44 +25,44 @@ BEGIN_PROVIDER [double precision, mo_pseudo_integrals, (mo_num,mo_num)] END_PROVIDER - -BEGIN_PROVIDER [double precision, mo_pseudo_integrals_local, (mo_num,mo_num)] - implicit none - BEGIN_DOC - ! Pseudopotential integrals in |MO| basis - END_DOC - - if (do_pseudo) then - call ao_to_mo( & - ao_pseudo_integrals_local, & - size(ao_pseudo_integrals_local,1), & - mo_pseudo_integrals_local, & - size(mo_pseudo_integrals_local,1) & - ) - else - mo_pseudo_integrals_local = 0.d0 - endif - -END_PROVIDER +! +!BEGIN_PROVIDER [double precision, mo_pseudo_integrals_local, (mo_num,mo_num)] +! implicit none +! BEGIN_DOC +! ! Pseudopotential integrals in |MO| basis +! END_DOC +! +! if (do_pseudo) then +! call ao_to_mo( & +! ao_pseudo_integrals_local, & +! size(ao_pseudo_integrals_local,1), & +! mo_pseudo_integrals_local, & +! size(mo_pseudo_integrals_local,1) & +! ) +! else +! mo_pseudo_integrals_local = 0.d0 +! endif +! +!END_PROVIDER -BEGIN_PROVIDER [double precision, mo_pseudo_integrals_non_local, (mo_num,mo_num)] - implicit none - BEGIN_DOC - ! Pseudopotential integrals in |MO| basis - END_DOC - - if (do_pseudo) then - call ao_to_mo( & - ao_pseudo_integrals_non_local, & - size(ao_pseudo_integrals_non_local,1), & - mo_pseudo_integrals_non_local, & - size(mo_pseudo_integrals_non_local,1) & - ) - else - mo_pseudo_integrals_non_local = 0.d0 - endif - -END_PROVIDER +!BEGIN_PROVIDER [double precision, mo_pseudo_integrals_non_local, (mo_num,mo_num)] +! implicit none +! BEGIN_DOC +! ! Pseudopotential integrals in |MO| basis +! END_DOC +! +! if (do_pseudo) then +! call ao_to_mo( & +! ao_pseudo_integrals_non_local, & +! size(ao_pseudo_integrals_non_local,1), & +! mo_pseudo_integrals_non_local, & +! size(mo_pseudo_integrals_non_local,1) & +! ) +! else +! mo_pseudo_integrals_non_local = 0.d0 +! endif +! +!END_PROVIDER