mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-24 21:41:46 +02:00
Localized MOs
This commit is contained in:
parent
8703a06c4f
commit
767906a051
75
plugins/Hartree_Fock/localize_mos.irp.f
Normal file
75
plugins/Hartree_Fock/localize_mos.irp.f
Normal file
@ -0,0 +1,75 @@
|
|||||||
|
program localize_mos
|
||||||
|
implicit none
|
||||||
|
integer :: rank, i,j,k
|
||||||
|
double precision, allocatable :: W(:,:)
|
||||||
|
double precision :: f, f_incr
|
||||||
|
|
||||||
|
allocate (W(ao_num,ao_num))
|
||||||
|
|
||||||
|
W = 0.d0
|
||||||
|
do k=1,elec_beta_num
|
||||||
|
do j=1,ao_num
|
||||||
|
do i=1,ao_num
|
||||||
|
W(i,j) = W(i,j) + mo_coef(i,k) * mo_coef(j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! call svd_mo(ao_num,elec_beta_num,W, size(W,1), &
|
||||||
|
! mo_coef(1,1),size(mo_coef,1))
|
||||||
|
call cholesky_mo(ao_num,elec_beta_num,W, size(W,1), &
|
||||||
|
mo_coef(1,1),size(mo_coef,1),1.d-6,rank)
|
||||||
|
print *, rank
|
||||||
|
|
||||||
|
if (elec_alpha_num>elec_alpha_num) then
|
||||||
|
W = 0.d0
|
||||||
|
do k=elec_beta_num+1,elec_alpha_num
|
||||||
|
do j=1,ao_num
|
||||||
|
do i=1,ao_num
|
||||||
|
W(i,j) = W(i,j) + mo_coef(i,k) * mo_coef(j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! call svd_mo(ao_num,elec_alpha_num-elec_beta_num,W, size(W,1), &
|
||||||
|
! mo_coef(1,1),size(mo_coef,1))
|
||||||
|
call cholesky_mo(ao_num,elec_alpha_num-elec_beta_num,W, size(W,1), &
|
||||||
|
mo_coef(1,elec_beta_num+1),size(mo_coef,1),1.d-6,rank)
|
||||||
|
print *, rank
|
||||||
|
endif
|
||||||
|
|
||||||
|
W = 0.d0
|
||||||
|
do k=elec_alpha_num+1,mo_tot_num
|
||||||
|
do j=1,ao_num
|
||||||
|
do i=1,ao_num
|
||||||
|
W(i,j) = W(i,j) + mo_coef(i,k) * mo_coef(j,k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! call svd_mo(ao_num,mo_tot_num-elec_alpha_num,W, size(W,1), &
|
||||||
|
! mo_coef(1,1),size(mo_coef,1))
|
||||||
|
call cholesky_mo(ao_num,mo_tot_num-elec_alpha_num,W, size(W,1), &
|
||||||
|
mo_coef(1,elec_alpha_num+1),size(mo_coef,1),1.d-6,rank)
|
||||||
|
print *, rank
|
||||||
|
mo_label = "Localized"
|
||||||
|
|
||||||
|
TOUCH mo_coef
|
||||||
|
|
||||||
|
W(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num)
|
||||||
|
integer :: iorder(mo_tot_num)
|
||||||
|
double precision :: s(mo_tot_num), swap(ao_num)
|
||||||
|
do k=1,mo_tot_num
|
||||||
|
iorder(k) = k
|
||||||
|
s(k) = Fock_matrix_diag_mo(k)
|
||||||
|
enddo
|
||||||
|
call dsort(s(1),iorder(1),elec_beta_num)
|
||||||
|
call dsort(s(elec_beta_num+1),iorder(elec_beta_num+1),elec_alpha_num-elec_beta_num)
|
||||||
|
call dsort(s(elec_alpha_num+1),iorder(elec_alpha_num+1),mo_tot_num-elec_alpha_num)
|
||||||
|
do k=1,mo_tot_num
|
||||||
|
mo_coef(1:ao_num,k) = W(1:ao_num,iorder(k))
|
||||||
|
print *, k, s(k)
|
||||||
|
enddo
|
||||||
|
call save_mos
|
||||||
|
|
||||||
|
end
|
@ -20,6 +20,12 @@ doc: MO occupation numbers
|
|||||||
interface: ezfio
|
interface: ezfio
|
||||||
size: (mo_basis.mo_tot_num)
|
size: (mo_basis.mo_tot_num)
|
||||||
|
|
||||||
|
[mo_class]
|
||||||
|
type: character*(32)
|
||||||
|
doc: c: core, i: inactive, a: active, v: virtual, d: deleted
|
||||||
|
interface: ezfio, provider
|
||||||
|
size: (mo_basis.mo_tot_num)
|
||||||
|
|
||||||
[ao_md5]
|
[ao_md5]
|
||||||
type: character*(32)
|
type: character*(32)
|
||||||
doc: Ao_md5
|
doc: Ao_md5
|
||||||
|
@ -1,8 +1,20 @@
|
|||||||
subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
|
subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Cholesky decomposition of AO Density matrix to
|
! Cholesky decomposition of AO Density matrix
|
||||||
! generate MOs
|
!
|
||||||
|
! n : Number of AOs
|
||||||
|
|
||||||
|
! m : Number of MOs
|
||||||
|
!
|
||||||
|
! P(LDP,n) : Density matrix in AO basis
|
||||||
|
!
|
||||||
|
! C(LDC,m) : MOs
|
||||||
|
!
|
||||||
|
! tol_in : tolerance
|
||||||
|
!
|
||||||
|
! rank : Nomber of local MOs (output)
|
||||||
|
!
|
||||||
END_DOC
|
END_DOC
|
||||||
integer, intent(in) :: n,m, LDC, LDP
|
integer, intent(in) :: n,m, LDC, LDP
|
||||||
double precision, intent(in) :: P(LDP,n)
|
double precision, intent(in) :: P(LDP,n)
|
||||||
@ -15,9 +27,6 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
|
|||||||
integer :: ipiv(n)
|
integer :: ipiv(n)
|
||||||
double precision:: tol
|
double precision:: tol
|
||||||
double precision, allocatable :: W(:,:), work(:)
|
double precision, allocatable :: W(:,:), work(:)
|
||||||
!DEC$ ATTRIBUTES ALIGN: 32 :: W
|
|
||||||
!DEC$ ATTRIBUTES ALIGN: 32 :: work
|
|
||||||
!DEC$ ATTRIBUTES ALIGN: 32 :: ipiv
|
|
||||||
|
|
||||||
allocate(W(LDC,n),work(2*n))
|
allocate(W(LDC,n),work(2*n))
|
||||||
tol=tol_in
|
tol=tol_in
|
||||||
@ -41,40 +50,37 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
|
|||||||
deallocate(W,work)
|
deallocate(W,work)
|
||||||
end
|
end
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ]
|
subroutine svd_mo(n,m,P,LDP,C,LDC)
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Density matrix in MO basis
|
! Singular value decomposition of the AO Density matrix
|
||||||
END_DOC
|
!
|
||||||
integer :: i,j,k
|
! n : Number of AOs
|
||||||
mo_density_matrix = 0.d0
|
|
||||||
do k=1,mo_tot_num
|
|
||||||
if (mo_occ(k) == 0.d0) then
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
do j=1,ao_num
|
|
||||||
do i=1,ao_num
|
|
||||||
mo_density_matrix(i,j) = mo_density_matrix(i,j) + &
|
|
||||||
mo_occ(k) * mo_coef(i,k) * mo_coef(j,k)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ]
|
! m : Number of MOs
|
||||||
implicit none
|
!
|
||||||
BEGIN_DOC
|
! P(LDP,n) : Density matrix in AO basis
|
||||||
! Density matrix in MO basis (virtual MOs)
|
!
|
||||||
|
! C(LDC,m) : MOs
|
||||||
|
!
|
||||||
|
! tol_in : tolerance
|
||||||
|
!
|
||||||
|
! rank : Nomber of local MOs (output)
|
||||||
|
!
|
||||||
END_DOC
|
END_DOC
|
||||||
integer :: i,j,k
|
integer, intent(in) :: n,m, LDC, LDP
|
||||||
mo_density_matrix_virtual = 0.d0
|
double precision, intent(in) :: P(LDP,n)
|
||||||
do k=1,mo_tot_num
|
double precision, intent(out) :: C(LDC,m)
|
||||||
do j=1,ao_num
|
|
||||||
do i=1,ao_num
|
integer :: info
|
||||||
mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + &
|
integer :: i,k
|
||||||
(2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k)
|
integer :: ipiv(n)
|
||||||
enddo
|
double precision:: tol
|
||||||
enddo
|
double precision, allocatable :: W(:,:), work(:)
|
||||||
enddo
|
|
||||||
END_PROVIDER
|
allocate(W(LDC,n),work(2*n))
|
||||||
|
call svd(P,LDP,C,LDC,W,size(W,1),m,n)
|
||||||
|
|
||||||
|
deallocate(W,work)
|
||||||
|
end
|
||||||
|
|
||||||
|
@ -258,3 +258,4 @@ subroutine mix_mo_jk(j,k)
|
|||||||
enddo
|
enddo
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -169,7 +169,7 @@ END_PROVIDER
|
|||||||
'Nuclear repulsion energy')
|
'Nuclear repulsion energy')
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ character*(128), element_name, (54)]
|
BEGIN_PROVIDER [ character*(128), element_name, (78)]
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Array of the name of element, sorted by nuclear charge (integer)
|
! Array of the name of element, sorted by nuclear charge (integer)
|
||||||
END_DOC
|
END_DOC
|
||||||
@ -227,4 +227,29 @@ BEGIN_PROVIDER [ character*(128), element_name, (54)]
|
|||||||
element_name(52) = 'Te'
|
element_name(52) = 'Te'
|
||||||
element_name(53) = 'I'
|
element_name(53) = 'I'
|
||||||
element_name(54) = 'Xe'
|
element_name(54) = 'Xe'
|
||||||
|
element_name(55) = 'Cs'
|
||||||
|
element_name(56) = 'Ba'
|
||||||
|
element_name(57) = 'La'
|
||||||
|
element_name(58) = 'Ce'
|
||||||
|
element_name(59) = 'Pr'
|
||||||
|
element_name(60) = 'Nd'
|
||||||
|
element_name(61) = 'Pm'
|
||||||
|
element_name(62) = 'Sm'
|
||||||
|
element_name(63) = 'Eu'
|
||||||
|
element_name(64) = 'Gd'
|
||||||
|
element_name(65) = 'Tb'
|
||||||
|
element_name(66) = 'Dy'
|
||||||
|
element_name(67) = 'Ho'
|
||||||
|
element_name(68) = 'Er'
|
||||||
|
element_name(69) = 'Tm'
|
||||||
|
element_name(70) = 'Yb'
|
||||||
|
element_name(71) = 'Lu'
|
||||||
|
element_name(72) = 'Hf'
|
||||||
|
element_name(73) = 'Ta'
|
||||||
|
element_name(74) = 'W'
|
||||||
|
element_name(75) = 'Re'
|
||||||
|
element_name(76) = 'Os'
|
||||||
|
element_name(77) = 'Ir'
|
||||||
|
element_name(78) = 'Pt'
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
Loading…
Reference in New Issue
Block a user