LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

Hartree-Fock with less MOs than AOs OK

Browse Source
This commit is contained in:
Anthony Scemama 2015-11-27 19:14:36 +01:00
parent e1e2a1d3e4
commit dd78dd4324
2 changed files with 89 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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,mo_tot_num,elec_alpha_num,1.d0, &
call dgemm('N','T',ao_num,ao_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,mo_tot_num,elec_beta_num,1.d0, &
call dgemm('N','T',ao_num,ao_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))

129
plugins/Hartree_Fock/diagonalize_fock.irp.f
View File

@ -10,58 +10,103 @@
integer, allocatable :: iwork(:)
double precision, allocatable :: work(:), F(:,:), S(:,:)
allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
do j=1,ao_num
do i=1,ao_num
S(i,j) = ao_overlap(i,j)
F(i,j) = Fock_matrix_ao(i,j)
enddo
enddo
n = ao_num
lwork = 1+6*n + 2*n*n
liwork = 3 + 5*n
allocate(work(lwork), iwork(liwork) )
if (mo_tot_num == ao_num) then
! Solve H.C = E.S.C in AO basis set
lwork = -1
liwork = -1
allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
do j=1,ao_num
do i=1,ao_num
S(i,j) = ao_overlap(i,j)
F(i,j) = Fock_matrix_ao(i,j)
enddo
enddo
call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
! call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),&
! diagonal_Fock_matrix_mo, work, lwork, info)
n = ao_num
lwork = 1+6*n + 2*n*n
liwork = 3 + 5*n
allocate(work(lwork), iwork(liwork) )
lwork = -1
liwork = -1
call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
lwork = int(work(1))
liwork = iwork(1)
deallocate(work,iwork)
allocate(work(lwork), iwork(liwork) )
! deallocate(work)
! allocate(work(lwork))
if (info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
lwork = int(work(1))
liwork = iwork(1)
deallocate(work,iwork)
allocate(work(lwork), iwork(liwork) )
call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
call dsygvd(1,'v','u',ao_num,F,size(F,1),S,size(S,1),&
diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
! call dsygv(1, 'v', 'u',ao_num,F,size(F,1),S,size(S,1),&
! diagonal_Fock_matrix_mo, work, lwork, info)
if (info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
do j=1,mo_tot_num
do i=1,ao_num
eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
enddo
enddo
if (info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
do j=1,mo_tot_num
do i=1,ao_num
eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
enddo
enddo
deallocate(work, iwork, F, S)
else
! Solve H.C = E.C in MO basis set
allocate( F(mo_tot_num_align,mo_tot_num) )
do j=1,mo_tot_num
do i=1,mo_tot_num
F(i,j) = Fock_matrix_mo(i,j)
enddo
enddo
n = mo_tot_num
lwork = 1+6*n + 2*n*n
liwork = 3 + 5*n
allocate(work(lwork), iwork(liwork) )
lwork = -1
liwork = -1
call dsyevd( 'V', 'U', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
lwork = int(work(1))
liwork = iwork(1)
deallocate(work,iwork)
allocate(work(lwork), iwork(liwork) )
call dsyevd( 'V', 'U', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' failed : ', info
stop 1
endif
call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, &
mo_coef, size(mo_coef,1), F, size(F,1), &
0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
deallocate(work, iwork, F)
endif
deallocate(work, iwork, F, S)
END_PROVIDER
BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]