10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-03-31 14:01:36 +02:00

Add a few functions for spherical/cartesian arrays

This commit is contained in:
Anthony Scemama 2025-03-21 23:13:55 +01:00
parent 66ce7bba22
commit 976c5aeb5e
6 changed files with 183 additions and 44 deletions

View File

@ -1 +0,0 @@
/home/scemama/qp2/plugins/qp_plugins_lct/stable/rsdft_cipsi/qp_cipsi_rsh

View File

@ -1,6 +1,9 @@
#!/usr/bin/env python3
# Computes the error on the excitation energy of a CIPSI run.
# see "QUESTDB: a database of highly-accurate excitation energies for the electronic structure community"
# doi:10.1002/wcms.1517 or arXiv:2011.14675
def student(p,df):
import scipy

View File

@ -28,11 +28,15 @@ BEGIN_PROVIDER [ integer, ao_sphe_num ]
! Number of spherical AOs
END_DOC
integer :: n, i
ao_sphe_num=0
do i=1,shell_num
n = shell_ang_mom(i)
ao_sphe_num += 2*n+1
enddo
if (ao_cartesian) then
ao_sphe_num = ao_num
else
ao_sphe_num=0
do i=1,shell_num
n = shell_ang_mom(i)
ao_sphe_num += 2*n+1
enddo
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, ao_sphe_shell, (ao_sphe_num) ]

View File

@ -14,45 +14,55 @@
prev = 0
ao_cart_to_sphe_coef(:,:) = 0.d0
ao_cart_to_sphe_normalization(:) = 1.d0
! Assume order provided by ao_power_index
i = 1
ao_sphe_count = 0
do while (i <= ao_num)
select case ( ao_l(i) )
case (0)
ao_sphe_count += 1
ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0
ao_cart_to_sphe_normalization(i) = 1.d0
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k)
if (ao_cartesian) then
! Identity matrix
do i=1,ao_sphe_num
ao_cart_to_sphe_coef(i,i) = 1.d0
enddo
else
! Assume order provided by ao_power_index
i = 1
ao_sphe_count = 0
do while (i <= ao_num)
select case ( ao_l(i) )
case (0)
ao_sphe_count += 1
ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0
ao_cart_to_sphe_normalization(i) = 1.d0
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k)
enddo
enddo
enddo
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j)
enddo
i += size(cart_to_sphe_$SHELL,1)
ao_sphe_count += size(cart_to_sphe_$SHELL,2)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in ao_cart_to_sphe : angular momentum too high'
end select
enddo
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j)
enddo
i += size(cart_to_sphe_$SHELL,1)
ao_sphe_count += size(cart_to_sphe_$SHELL,2)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in ao_cart_to_sphe : angular momentum too high'
end select
enddo
endif
if (ao_sphe_count /= ao_sphe_num) then
call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num")

View File

@ -31,3 +31,42 @@ BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ]
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ integer, mo_symmetry , (mo_num) ]
implicit none
BEGIN_DOC
! MOs with the same integer belong to the same irrep.
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
if (size(mo_symmetry) == 0) return
call ezfio_has_mo_basis_mo_symmetry(has)
if (has) then
! write(6,'(A)') '.. >>>>> [ IO READ: mo_symmetry ] <<<<< ..'
call ezfio_get_mo_basis_mo_symmetry(mo_symmetry)
else
mo_symmetry(:) = 1
call ezfio_set_mo_basis_mo_symmetry(mo_symmetry)
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( mo_symmetry, (mo_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read mo_symmetry with MPI'
endif
IRP_ENDIF
! call write_time(6)
END_PROVIDER

View File

@ -73,3 +73,87 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_from_mo, (ao_num, ao_num)]
END_DOC
call mo_to_ao(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_num)
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