From 5f357ee68b6511ab6f408bd6916cc6a676b332d3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 26 Jun 2017 18:11:49 +0200 Subject: [PATCH] Working on Slater dressing --- plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg | 8 +- .../LinearSystem.irp.f | 2 - .../SCF_dressed.irp.f | 6 +- .../Hartree_Fock_SlaterDressed/dressing.irp.f | 140 +++++++++--------- .../integrals.irp.f | 11 +- .../Hartree_Fock_SlaterDressed/slater.irp.f | 5 +- src/MO_Basis/ao_ortho_canonical.irp.f | 6 +- src/MO_Basis/mos.irp.f | 6 +- 8 files changed, 101 insertions(+), 83 deletions(-) diff --git a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg index 6a7d5fdc..df55c554 100644 --- a/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg +++ b/plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg @@ -14,7 +14,13 @@ interface: ezfio, provider type: double precision doc: Orthogonal AO basis size: (ao_basis.ao_num,ao_basis.ao_num) -interface: ezfio, provider +interface: ezfio + +[ao_orthoSlaOverlap] +type: double precision +doc: Orthogonal AO basis +size: (ao_basis.ao_num,nuclei.nucl_num) +interface: ezfio diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f index 4ccd09b0..61a5a49c 100644 --- a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f @@ -14,7 +14,6 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ] ! Projector do mu=1,mo_tot_num cusp_A(A,B) += AO_orthoSlaOverlap_matrix(mu,B) * ao_ortho_value_at_nucl(mu,A) -! cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A) enddo enddo enddo @@ -60,6 +59,5 @@ BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ] stop 'dgetrs failed' endif - END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f index d0998a61..8befdb91 100644 --- a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f @@ -86,8 +86,10 @@ subroutine run enddo mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num) SOFT_TOUCH mo_coef slater_coef -! call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef) -! call save_mos + call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef) + call ezfio_set_Hartree_Fock_SlaterDressed_projector(ao_ortho_canonical_coef(1:ao_num,1:ao_num)) + call ezfio_set_Hartree_Fock_SlaterDressed_ao_orthoSlaOverlap(AO_orthoSlaOverlap_matrix) + call save_mos print *, 'ci' print *, mo_coef(1:ao_num,1) print *, 'cAi' diff --git a/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f index f13f4a13..f42a0b0f 100644 --- a/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/dressing.irp.f @@ -91,74 +91,82 @@ BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num_align,mo_tot_num) cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num) slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num) return + + else + + + do idx_dressing=1,mo_tot_num + + if (idx_dressing>1) then + TOUCH idx_dressing + endif + + do j=1,mo_tot_num + do i=1,mo_tot_num + F(i,j) = Fock_matrix_mo(i,j) + enddo + enddo + + do j=1,mo_tot_num + do i=1,ao_num + M(i,j) = mo_coef(i,j) + enddo + enddo + + integer :: it + do it=1,128 + + ! print *, 'X', ao_ortho_canonical_coef(1:ao_num,1:ao_num) + ! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) + ! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num) + ! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num) + ! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) + ! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num) + ! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) + ! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) + ! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num) + ! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) + ! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num) + ! print *, 'Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) + do j=1,mo_tot_num + do i=1,mo_tot_num + Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j) + enddo + enddo + do i=1,mo_tot_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo + ! print *, 'Dressed Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num) + double precision :: conv + conv = 0.d0 + do j=1,mo_tot_num + do i=1,mo_tot_num + if (i==j) cycle + conv = max(conv,Fock_matrix_mo(i,j)) + enddo + enddo + TOUCH Fock_matrix_mo Fock_matrix_diag_mo + + mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) + TOUCH mo_coef + !print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) + !print *, '-----' + print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv) + if (conv < 1.d-5) exit + !stop + + enddo + cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing) + slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing) + enddo + + idx_dressing = 1 + mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num) + soft_TOUCH mo_coef idx_dressing slater_coef + endif - do idx_dressing=1,mo_tot_num - - if (idx_dressing>1) then - TOUCH idx_dressing - endif - - do j=1,mo_tot_num - do i=1,mo_tot_num - F(i,j) = Fock_matrix_mo(i,j) - enddo - enddo - - do j=1,mo_tot_num - do i=1,ao_num - M(i,j) = mo_coef(i,j) - enddo - enddo - - integer :: it - do it=1,128 - -! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) -! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num) -! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num) -! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) -! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) -! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num) -! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) -! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num) -! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num) -! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num) -! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num) - do j=1,mo_tot_num - do i=1,mo_tot_num - Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j) - enddo - enddo - do i=1,mo_tot_num - Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) - enddo - double precision :: conv - conv = 0.d0 - do j=1,mo_tot_num - do i=1,mo_tot_num - if (i==j) cycle - conv = max(conv,Fock_matrix_mo(i,j)) - enddo - enddo - TOUCH Fock_matrix_mo Fock_matrix_diag_mo - - mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) - TOUCH mo_coef -! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num) -! print *, '-----' - print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv) - if (conv < 1.d-5) exit - - enddo - cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing) - slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing) - enddo - - idx_dressing = 1 - mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num) - soft_TOUCH mo_coef idx_dressing slater_coef - END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f index fa4530a9..680f3e1f 100644 --- a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f @@ -344,7 +344,7 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result) ss = k*ss ! Print result - write(*,*) ss +! write(*,*) ss result = 0.d0 end @@ -520,6 +520,7 @@ BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ] END_PROVIDER + SUBST [ X ] Overlap ;; @@ -616,9 +617,9 @@ BEGIN_PROVIDER [ double precision, AO_orthoSlaH_matrix, (ao_num, nucl_num) ] BEGIN_DOC ! END_DOC - call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, & - ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & - GauSlaH_matrix, size(GauSlaH_matrix,1), & - 0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1)) + call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, & + ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & + GauSlaH_matrix, size(GauSlaH_matrix,1), & + 0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1)) END_PROVIDER diff --git a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f index 8a11b9b1..27e4f10d 100644 --- a/plugins/Hartree_Fock_SlaterDressed/slater.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/slater.irp.f @@ -8,7 +8,10 @@ BEGIN_PROVIDER [ double precision, slater_expo, (nucl_num) ] if (exists) then slater_expo(1:nucl_num) = slater_expo_ezfio(1:nucl_num) else - slater_expo(1:nucl_num) = nucl_charge(1:nucl_num) + integer :: i + do i=1,nucl_num + slater_expo(i) = nucl_charge(i) + enddo call ezfio_set_Hartree_Fock_SlaterDressed_slater_expo_ezfio(slater_expo) endif END_PROVIDER diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 2184ce4a..d77374a8 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -106,9 +106,9 @@ END_PROVIDER ao_ortho_canonical_coef(i,i) = 1.d0 enddo -!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) -!ao_ortho_canonical_num=ao_num -!return +call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) +ao_ortho_canonical_num=ao_num +return if (ao_cartesian) then diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 6906874b..8962dd00 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -75,7 +75,7 @@ BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_ ! ! C^(-1).C_mo END_DOC - call dgemm('T','N',ao_num,mo_tot_num,ao_num,1.d0, & + call dgemm('N','N',ao_num,mo_tot_num,ao_num,1.d0, & ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),& mo_coef, size(mo_coef,1), 0.d0, & mo_coef_in_ao_ortho_basis, size(mo_coef_in_ao_ortho_basis,1)) @@ -290,13 +290,13 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) allocate ( T(ao_num_align,ao_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - call dgemm('N','N', ao_num, ao_num, ao_num, & + call dgemm('T','N', ao_num, ao_num, ao_num, & 1.d0, & ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1), & A_ao,LDA_ao, & 0.d0, T, ao_num_align) - call dgemm('N','T', ao_num, ao_num, ao_num, 1.d0, & + call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, & T, size(T,1), & ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),& 0.d0, A, LDA)