From ad4dc469cc156a870f258eb0d9dffb5490e31f9d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 27 Nov 2015 11:42:46 +0100 Subject: [PATCH] SCF works with non-square matrices --- plugins/Hartree_Fock/Fock_matrix.irp.f | 47 ++++++++++++++++--- .../Hartree_Fock/HF_density_matrix_ao.irp.f | 4 +- 2 files changed, 42 insertions(+), 9 deletions(-) diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 82540de3..12ee276b 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -344,30 +344,47 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] enddo enddo else - double precision, allocatable :: T(:,:), M(:,:) + double precision, allocatable :: T(:,:), M(:,:) + integer :: ierr ! F_ao = S C F_mo C^t S - allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + allocate (T(ao_num_align,ao_num),M(ao_num_align,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)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + +! 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) endif END_PROVIDER @@ -380,23 +397,39 @@ subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) 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_align,ao_num),M(ao_num_align,ao_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + allocate (T(ao_num_align,ao_num),M(ao_num_align,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)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + +! 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), & diff --git a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f index e8585f59..597b88c7 100644 --- a/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/plugins/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_ ! S^-1 x Alpha density matrix in the AO basis x S^-1 END_DOC - call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + call dgemm('N','T',ao_num,mo_tot_num,elec_alpha_num,1.d0, & mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), 0.d0, & HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) @@ -17,7 +17,7 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_ ! S^-1 Beta density matrix in the AO basis x S^-1 END_DOC - call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + call dgemm('N','T',ao_num,mo_tot_num,elec_beta_num,1.d0, & mo_coef, size(mo_coef,1), & mo_coef, size(mo_coef,1), 0.d0, & HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1))