From 767906a051bcff97e9ca3aa02d5f63f9efd531cf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Mar 2017 12:11:54 +0100 Subject: [PATCH] Localized MOs --- plugins/Hartree_Fock/localize_mos.irp.f | 75 +++++++++++++++++++++++ src/MO_Basis/EZFIO.cfg | 8 ++- src/MO_Basis/cholesky_mo.irp.f | 80 +++++++++++++------------ src/MO_Basis/mos.irp.f | 1 + src/Nuclei/nuclei.irp.f | 27 ++++++++- 5 files changed, 152 insertions(+), 39 deletions(-) create mode 100644 plugins/Hartree_Fock/localize_mos.irp.f diff --git a/plugins/Hartree_Fock/localize_mos.irp.f b/plugins/Hartree_Fock/localize_mos.irp.f new file mode 100644 index 00000000..8a665c64 --- /dev/null +++ b/plugins/Hartree_Fock/localize_mos.irp.f @@ -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 diff --git a/src/MO_Basis/EZFIO.cfg b/src/MO_Basis/EZFIO.cfg index 5aec39e0..368b70a0 100644 --- a/src/MO_Basis/EZFIO.cfg +++ b/src/MO_Basis/EZFIO.cfg @@ -20,7 +20,13 @@ doc: MO occupation numbers interface: ezfio 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] type: character*(32) doc: Ao_md5 -interface: ezfio \ No newline at end of file +interface: ezfio diff --git a/src/MO_Basis/cholesky_mo.irp.f b/src/MO_Basis/cholesky_mo.irp.f index 97b6abd2..65184c1e 100644 --- a/src/MO_Basis/cholesky_mo.irp.f +++ b/src/MO_Basis/cholesky_mo.irp.f @@ -1,8 +1,20 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) implicit none BEGIN_DOC -! Cholesky decomposition of AO Density matrix to -! generate MOs +! Cholesky decomposition of AO Density matrix +! +! 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 integer, intent(in) :: n,m, LDC, LDP 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) double precision:: tol 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)) tol=tol_in @@ -41,40 +50,37 @@ subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) deallocate(W,work) 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 BEGIN_DOC - ! Density matrix in MO basis - END_DOC - integer :: i,j,k - 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 +! Singular value decomposition of the AO Density matrix +! +! n : Number of AOs -BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ] - implicit none - BEGIN_DOC - ! Density matrix in MO basis (virtual MOs) +! 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 - integer :: i,j,k - mo_density_matrix_virtual = 0.d0 - do k=1,mo_tot_num - do j=1,ao_num - do i=1,ao_num - mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + & - (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k) - enddo - enddo - enddo -END_PROVIDER + integer, intent(in) :: n,m, LDC, LDP + double precision, intent(in) :: P(LDP,n) + double precision, intent(out) :: C(LDC,m) + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:) + + allocate(W(LDC,n),work(2*n)) + call svd(P,LDP,C,LDC,W,size(W,1),m,n) + + deallocate(W,work) +end diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 69abf7b3..19835395 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -258,3 +258,4 @@ subroutine mix_mo_jk(j,k) enddo end + diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index c4729713..b4da5fb1 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -169,7 +169,7 @@ END_PROVIDER 'Nuclear repulsion energy') END_PROVIDER -BEGIN_PROVIDER [ character*(128), element_name, (54)] +BEGIN_PROVIDER [ character*(128), element_name, (78)] BEGIN_DOC ! Array of the name of element, sorted by nuclear charge (integer) END_DOC @@ -227,4 +227,29 @@ BEGIN_PROVIDER [ character*(128), element_name, (54)] element_name(52) = 'Te' element_name(53) = 'I' 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