PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-27 18:54:53 +02:00
Code Releases Activity

added src/ao_two_e_ints/screening.irp.f

Browse Source
This commit is contained in:
eginer 2025-04-25 14:09:39 +02:00
parent 8ce2be608c
commit 96cfcb2498
10 changed files with 344 additions and 383 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/ao_cart_basis/EZFIO.cfg
View File

@ -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|

2
src/ao_cart_two_e_ints/cholesky.irp.f
View File

@ -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

150
src/ao_two_e_ints/ao_bi_integrals.irp.f
View File

@ -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

154
src/ao_two_e_ints/map_integrals.irp.f
View File

@ -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 ]

31
src/ao_two_e_ints/schwartz.irp.f Normal file
View File

@ -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

16
src/ao_two_e_ints/screening.irp.f Normal file
View File

@ -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

152
src/hartree_fock/fock_matrix_hf.irp.f
View File

@ -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

64
src/mo_basis/mos.irp.f
View File

@ -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

78
src/mo_one_e_ints/ao_to_mo.irp.f
View File

@ -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

74
src/mo_one_e_ints/pot_mo_pseudo_ints.irp.f
View File

@ -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