LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-19 19:52:15 +02:00
Code Issues Releases Wiki Activity

SCF works with non-square matrices

Browse Source
This commit is contained in:
Anthony Scemama 2015-11-27 11:42:46 +01:00
parent 88c47b1f73
commit ad4dc469cc
2 changed files with 42 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
plugins/Hartree_Fock/Fock_matrix.irp.f
View File

@ -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), &

4
plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
View File

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