diff --git a/scripts/qp_cipsi_rsh b/scripts/qp_cipsi_rsh deleted file mode 120000 index c3d4376b..00000000 --- a/scripts/qp_cipsi_rsh +++ /dev/null @@ -1 +0,0 @@ -/home/scemama/qp2/plugins/qp_plugins_lct/stable/rsdft_cipsi/qp_cipsi_rsh \ No newline at end of file diff --git a/scripts/qp_exc_energy.py b/scripts/qp_exc_energy.py index e08866e3..e3060349 100755 --- a/scripts/qp_exc_energy.py +++ b/scripts/qp_exc_energy.py @@ -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 diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 02eedf53..ef871538 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -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) ] diff --git a/src/ao_one_e_ints/ao_ortho_canonical.irp.f b/src/ao_one_e_ints/ao_ortho_canonical.irp.f index 523e49f7..e015c89e 100644 --- a/src/ao_one_e_ints/ao_ortho_canonical.irp.f +++ b/src/ao_one_e_ints/ao_ortho_canonical.irp.f @@ -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") diff --git a/src/mo_basis/mo_class.irp.f b/src/mo_basis/mo_class.irp.f index 7705e414..8017e119 100644 --- a/src/mo_basis/mo_class.irp.f +++ b/src/mo_basis/mo_class.irp.f @@ -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 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 7ebc4638..72c1328d 100644 --- a/src/mo_one_e_ints/ao_to_mo.irp.f +++ b/src/mo_one_e_ints/ao_to_mo.irp.f @@ -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 + + +