From 7b55cfad054b2e8b7037d00c903d1a448498d28c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 16 Jun 2017 15:35:52 +0200 Subject: [PATCH] Cusp dressing almost OK --- plugins/Hartree_Fock/EZFIO.cfg | 4 +- plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f | 8 + .../LinearSystem.irp.f | 1 + .../SCF_dressed.irp.f | 13 +- .../integrals.irp.f | 550 +++++++++++++----- src/MOGuess/mo_ortho_lowdin.irp.f | 5 +- src/MO_Basis/ao_ortho_canonical.irp.f | 9 + src/MO_Basis/mos.irp.f | 71 +++ src/Utils/LinearAlgebra.irp.f | 6 +- 9 files changed, 529 insertions(+), 138 deletions(-) diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index 35879330..8db83f92 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -26,13 +26,13 @@ default: 1.e-12 type: Strictly_positive_int doc: Maximum number of SCF iterations interface: ezfio,provider,ocaml -default: 128 +default: 500 [level_shift] type: Positive_float doc: Energy shift on the virtual MOs to improve SCF convergence interface: ezfio,provider,ocaml -default: 0.0 +default: 0.3 [scf_algorithm] type: character*(32) diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f index 988528f3..18201c96 100644 --- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f +++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f @@ -71,7 +71,12 @@ END_DOC Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0 Fock_matrix_AO_beta = Fock_matrix_AO*0.5d0 + level_shift_save = level_shift + if (max_error_DIIS > 0.1d0) then + level_shift += 0.5d0 + endif TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta + level_shift = level_shift_save endif @@ -97,6 +102,9 @@ END_DOC double precision :: level_shift_save level_shift_save = level_shift + if (max_error_DIIS > 0.5d0) then + level_shift += 0.5d0 + endif do while (Delta_Energy_SCF > 0.d0) level_shift = level_shift + 0.1d0 MO_coef = eigenvectors_Fock_matrix_MO diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f index c7bd816a..490e344c 100644 --- a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f @@ -11,6 +11,7 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ] cusp_A(A,A) = slater_expo(A)/nucl_charge(A) * slater_value_at_nucl(A,A) do B=1,nucl_num cusp_A(A,B) -= slater_value_at_nucl(B,A) + ! Projector do mu=1,mo_tot_num cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A) enddo diff --git a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f index 1d8ea0df..a8c1351a 100644 --- a/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f @@ -76,14 +76,21 @@ subroutine run double precision :: EHF integer :: i_it, i, j, k - EHF = HF_energy - mo_label = "CuspDressed" + + print *, HF_energy + call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_C) -! Choose SCF algorithm + do i=1,ao_num + print *, mo_coef(i,1), cusp_corrected_mos(i,1) + enddo + mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num) + TOUCH mo_coef + EHF = HF_energy + print *, HF_energy ! call Roothaan_Hall_SCF end diff --git a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f index 08b2eb56..1111055b 100644 --- a/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f +++ b/plugins/Hartree_Fock_SlaterDressed/integrals.irp.f @@ -24,44 +24,93 @@ subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result) ! calculate the length AB between the two centers and other usful quantities - AB = (cGau(1)-cSla(1))**2d0 + (cGau(2)-cSla(2))**2d0 + (cGau(3)-cSla(3))**2d0 - AB = sqrt(AB) + AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2 + AB = dsqrt(AB) AxBx = (cGau(1)-cSla(1))/2d0 AyBy = (cGau(2)-cSla(2))/2d0 AzBz = (cGau(3)-cSla(3))/2d0 + ds = 0.d0 ! intermediate variables - t = expSla*sqrt(0.25d0/expGau) - u = sqrt(expGau)*AB + t = expSla*dsqrt(0.25d0/expGau) + u = dsqrt(expGau)*AB + double precision :: d, et2 if(AB > 0d0) then ! (s|s) - ss = (t+u)*erfc(t+u)*exp(2d0*t*(t+u)) - (t-u)*erfc(t-u)*exp(2d0*t*(t-u)) + ss = 0.d0 + + d = derfc(t+u) + if (dabs(d) > 1.d-30) then + ss = (t+u)*d*dexp(2d0*t*(t+u)) + endif + + d = derfc(t-u) + if (dabs(d) > 1.d-30) then + ss -= (t-u)*d*dexp(2d0*t*(t-u)) + endif ! (p|s) - ps = (exp(t**2d0-u**2d0)*(-4d0*t+sqrt(pi)*(exp((t-u)**2d0)*(1d0+2d0*t*(t-u))*erfc(t-u) & - + exp((t+u)**2d0)*(1d0+2d0*t*(t+u))*erfc(t+u))))/sqrt(pi) + ps = 0.d0 + if (t*t-u*u > 300.d0) then + et2 = huge(1.0) + else + et2 = dexp(t*t-u*u) + endif + if (et2 /= 0.d0) then + d = derfc(t-u) + if (d /= 0.d0) then + ps += dexp((t-u)**2)*(1d0+2d0*t*(t-u))*d + endif + d = derfc(t+u) + if (d /= 0.d0) then + ps += dexp((t+u)**2)*(1d0+2d0*t*(t+u))*d + endif + ps *= dsqrt(pi) + ps -= 4d0*t + ps *= et2/dsqrt(pi) + endif ! (d|s) - ds = 4d0*exp(2d0*t*(t-u))*t*(-((1d0+t**2d0-t*u)*erfc(t-u))+exp(4d0*t*u)*(1d0+t*(t+u))*erfc(t+u)) +! ds = 4d0*dexp(2d0*t*(t-u))*t*(-((1d0+t**2-t*u)*derfc(t-u))+dexp(4d0*t*u)*(1d0+t*(t+u))*derfc(t+u)) + ds = 0.d0 + d = derfc(t+u) + if (d /= 0.d0) then + ds = dexp(4d0*t*u)*(1d0+t*(t+u))*d + endif + d = derfc(t-u) + if (d /= 0.d0) then + ds -= (1d0+t*t-t*u)*d + endif + + if ( dabs(ds) > 1.d-100) then + ds *= 4d0*dexp(2d0*t*(t-u))*t + endif ! backward scaling ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0 - ps = ps/u**2d0-ss/u**3d0 + ps = ps/u**2-ss/u**3d0 ss = ss/u else ! concentric case - ss = 2d0*exp(t**2d0)*((-2d0*t)/sqrt(pi)+exp(t**2d0)*(1d0+2d0*t**2d0)*erfc(t)) - ps = (8d0*exp(t**2d0)*t*(-2d0*(1d0+t**2d0)+exp(t**2d0)*sqrt(pi)*t*(3d0+2d0*t**2d0)*erfc(t)))/(3d0*sqrt(pi)) + d = derfc(t) + if (d /= 0.d0) then + et2 = dexp(t*t) + ss = 2d0*et2*((-2d0*t)/dsqrt(pi)+et2*(1d0+2d0*t*t)*d) + ps = (8d0*et2*t*(-2d0*(1d0+t*t)+et2*dsqrt(pi)*t*(3d0+2d0*t*t)*d))/(3d0*dsqrt(pi)) + else + ss = 0.d0 + ps = 0.d0 + endif endif - k = t**3d0*exp(-t**2d0)*4d0*pi/expSla**(3d0/2d0) + k = t**3d0*dexp(-t*t)*4d0*pi/expSla**(3d0/2d0) ! (s|s) ss = k*ss @@ -76,9 +125,134 @@ subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result) ! (d|s) ds = k*ds - dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2d0*ds - dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2d0*ds - dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2d0*ds + dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2*ds + dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2*ds + dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2*ds + + dxys = AxBx*AyBy*ds + dxzs = AxBx*AzBz*ds + dyzs = AyBy*AzBz*ds + + select case (sum(aGau)) + case (0) + result = ss + + case (1) + if (aGau(1) == 1) then + result = pxs + else if (aGau(2) == 1) then + result = pys + else if (aGau(3) == 1) then + result = pzs + endif + + case (2) + if (aGau(1) == 2) then + result = dxxs + else if (aGau(2) == 2) then + result = dyys + else if (aGau(3) == 2) then + result = dzzs + else if (aGau(1)+aGau(2) == 2) then + result = dxys + else if (aGau(1)+aGau(3) == 2) then + result = dxzs + else if (aGau(2)+aGau(3) == 2) then + result = dyzs + endif + + case default + stop 'GauSlaOverlap not implemented' + + end select + +end +!***************************************************************************** + +!***************************************************************************** +subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla,result) + + implicit none + + BEGIN_DOC + ! Compute the kinetic energy integral between a Gaussian function + ! with arbitrary angular momemtum and a s-type Slater function + END_DOC + +! Input variables + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + double precision,intent(out) :: result + +! Final value of the integrals + double precision :: ss,ps,ds + double precision :: pxs,pys,pzs + double precision :: dxxs,dyys,dzzs,dxys,dxzs,dyzs + + double precision :: pi,E,AB,AxBx,AyBy,AzBz,t,u,k + + pi = 4d0*atan(1d0) + +! calculate the length AB between the two centers + + AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2 + AB = dsqrt(AB) + + AxBx = (cGau(1)-cSla(1))/2d0 + AyBy = (cGau(2)-cSla(2))/2d0 + AzBz = (cGau(3)-cSla(3))/2d0 + +! intermediate variables + + t = expSla*dsqrt(0.25d0/expGau) + u = dsqrt(expGau)*AB + + if(AB > 0d0) then + +! (s|s) + ss = (1d0+t*(t-u))*derfc(t-u)*dexp(2d0*t*(t-u)) - (1d0+t*(t+u))*derfc(t+u)*dexp(2d0*t*(t+u)) + +! (p|s) + ps = (dexp(t**2-2d0*t*u-u**2)*(4d0*dexp(2d0*t*u)*(1d0+t**2) & + + dsqrt(pi)*t*(-(dexp(t**2+u**2)*(3d0+2d0*t*(t-u))*derfc(t-u)) & + - dexp(2d0*t*u+(t+u)**2)*(3d0+2d0*t*(t+u))*derfc(t+u))))/dsqrt(pi) + +! (d|s) + ds = (-8d0*dexp(t**2-u**2)*u+4d0*dexp(2d0*t*(t-u))*dsqrt(pi)*t**2*((2d0+t**2-t*u)*derfc(t-u) & + - dexp(4d0*t*u)*(2d0+t*(t+u))*derfc(t+u)))/dsqrt(pi) + +! backward scaling + ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0 + ps = ps/u**2-ss/u**3d0 + ss = ss/u + + else + +! concentric case + ss = (4d0*dexp(t**2)*(1d0+t**2))/dsqrt(pi)-2d0*dexp(2d0*t**2)*t*(3d0+2d0*t**2)*derfc(t) + ps = (8d0*dexp(t**2)*(-1d0+4d0*t**2+2d0*t**4d0-dexp(t**2)*dsqrt(pi)*t**3d0*(5d0+2d0*t**2)*derfc(t)))/(3d0*dsqrt(pi)) + + endif + + k = expSla*dsqrt(expGau)*t**3d0*dexp(-t*t)*4d0*pi/expSla**(3d0/2d0) + +! (s|s) + ss = k*ss + +! (p|s) + ps = k*ps + + pxs = AxBx*ps + pys = AyBy*ps + pzs = AzBz*ps + +! (d|s) + ds = k*ds + + dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2*ds + dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2*ds + dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2*ds dxys = AxBx*AyBy*ds dxzs = AxBx*AzBz*ds @@ -121,108 +295,9 @@ end !***************************************************************************** -!***************************************************************************** -subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla) - - implicit none - - BEGIN_DOC - ! Compute the kinetic energy integral between a Gaussian function - ! with arbitrary angular momemtum and a s-type Slater function - END_DOC - -! Input variables - double precision,intent(in) :: expGau,expSla - double precision,intent(in) :: cGau(3),cSla(3) - integer,intent(in) :: aGau(3) - -! Final value of the integrals - double precision :: ss,ps,ds - double precision :: pxs,pys,pzs - double precision :: dxxs,dyys,dzzs,dxys,dxzs,dyzs - - double precision :: pi,E,AB,AxBx,AyBy,AzBz,t,u,k - - pi = 4d0*atan(1d0) - -! calculate the length AB between the two centers - - AB = (cGau(1)-cSla(1))**2d0 + (cGau(2)-cSla(2))**2d0 + (cGau(3)-cSla(3))**2d0 - AB = sqrt(AB) - - AxBx = (cGau(1)-cSla(1))/2d0 - AyBy = (cGau(2)-cSla(2))/2d0 - AzBz = (cGau(3)-cSla(3))/2d0 - -! intermediate variables - - t = expSla*sqrt(0.25d0/expGau) - u = sqrt(expGau)*AB - - if(AB > 0d0) then - -! (s|s) - ss = (1d0+t*(t-u))*erfc(t-u)*exp(2d0*t*(t-u)) - (1d0+t*(t+u))*erfc(t+u)*exp(2d0*t*(t+u)) - -! (p|s) - ps = (exp(t**2d0-2d0*t*u-u**2d0)*(4d0*exp(2d0*t*u)*(1d0+t**2d0) & - + sqrt(pi)*t*(-(exp(t**2d0+u**2d0)*(3d0+2d0*t*(t-u))*erfc(t-u)) & - - exp(2d0*t*u+(t+u)**2d0)*(3d0+2d0*t*(t+u))*erfc(t+u))))/sqrt(pi) - -! (d|s) - ds = (-8d0*exp(t**2d0-u**2d0)*u+4d0*exp(2d0*t*(t-u))*sqrt(pi)*t**2d0*((2d0+t**2d0-t*u)*erfc(t-u) & - - exp(4d0*t*u)*(2d0+t*(t+u))*erfc(t+u)))/sqrt(pi) - -! backward scaling - ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0 - ps = ps/u**2d0-ss/u**3d0 - ss = ss/u - - else - -! concentric case - ss = (4d0*exp(t**2d0)*(1d0+t**2d0))/sqrt(pi)-2d0*exp(2d0*t**2d0)*t*(3d0+2d0*t**2d0)*erfc(t) - ps = (8d0*exp(t**2d0)*(-1d0+4d0*t**2d0+2d0*t**4d0-exp(t**2)*sqrt(pi)*t**3d0*(5d0+2d0*t**2d0)*erfc(t)))/(3d0*sqrt(pi)) - - endif - - k = expSla*sqrt(expGau)*t**3d0*exp(-t**2d0)*4d0*pi/expSla**(3d0/2d0) - -! (s|s) - ss = k*ss - -! (p|s) - ps = k*ps - - pxs = AxBx*ps - pys = AyBy*ps - pzs = AzBz*ps - -! (d|s) - ds = k*ds - - dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2d0*ds - dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2d0*ds - dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2d0*ds - - dxys = AxBx*AyBy*ds - dxzs = AxBx*AzBz*ds - dyzs = AyBy*AzBz*ds - -! Print result - write(*,'(A12,F16.10)') & - '(s|T|s) = ',ss - write(*,'(A12,F16.10,3X,A12,F16.10,3X,A12,F16.10)') & - '(px|T|s) = ',pxs,'(py|T|s) = ',pys,'(pz|T|s) = ',pzs - write(*,'(A12,F16.10,3X,A12,F16.10,3X,A12,F16.10,3X,A12,F16.10,3X,A12,F16.10,3X,A12,F16.10)') & - '(dx2|T|s) = ',dxxs,'(dy2|T|s) = ',dyys,'(dz2|T|s) = ',dzzs,'(dxy|T|s) = ',dxys,'(dxz|T|s) = ',dxzs,'(dyz|T|s) = ',dyzs - -end -!***************************************************************************** - !***************************************************************************** -subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc) +subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result) implicit none @@ -237,6 +312,7 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc) integer,intent(in) :: aGau(3) double precision,intent(in) :: cNuc(3) double precision,intent(in) :: ZNuc + double precision,intent(out) :: result ! Final value of the overlap integral double precision :: ss,ps,ds,fs @@ -249,29 +325,31 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc) ! calculate the length AB between the two centers - AB = (cGau(1)-cSla(1))**2d0 + (cGau(2)-cSla(2))**2d0 + (cGau(3)-cSla(3))**2d0 - AB = sqrt(AB) + AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2 + AB = dsqrt(AB) ! intermediate variables - x = sqrt(expSla**2d0/(4d0*expGau)) - y = sqrt(expGau)*AB + x = dsqrt(expSla**2/(4d0*expGau)) + y = dsqrt(expGau)*AB if(AB > 0d0) then - ss = (1d0+x*(x+y))*erfc(x+y)*exp(2d0*x*(x+y)) - (1d0+x*(x-y))*erfc(x-y)*exp(2d0*x*(x-y)) + ss = (1d0+x*(x+y))*derfc(x+y)*dexp(2d0*x*(x+y)) - (1d0+x*(x-y))*derfc(x-y)*dexp(2d0*x*(x-y)) ss = ss/y else - ss = (4d0*E**x**2d0*(1d0+x**2d0))/sqrt(Pi)-2d0*E**(2d0*x**2d0)*x*(3d0+2d0*x**2d0)*Erfc(x) + ss = (4d0*E**x**2*(1d0+x**2))/dsqrt(Pi)-2d0*E**(2d0*x**2)*x*(3d0+2d0*x**2)*dErfc(x) endif - k = expSla*sqrt(expGau)*x**3d0*exp(-x**2)*4d0*pi/expSla**(3d0/2d0) + k = expSla*dsqrt(expGau)*x**3d0*dexp(-x*x)*4d0*pi/expSla**(3d0/2d0) ss = k*ss ! Print result write(*,*) ss + result = 0.d0 end !***************************************************************************** + double precision function BoysF0(t) implicit none double precision, intent(in) :: t @@ -280,7 +358,7 @@ double precision function BoysF0(t) pi = 4d0*atan(1d0) if(t > 0d0) then - BoysF0 = 0.5d0*sqrt(pi/t)*erf(sqrt(t)) + BoysF0 = 0.5d0*dsqrt(pi/t)*derf(dsqrt(t)) else BoysF0 = 1d0 endif @@ -288,8 +366,81 @@ double precision function BoysF0(t) end !***************************************************************************** +!TODO +subroutine GauSlaOverlap_write(expGau,cGau,aGau,expSla,cSla,result,iunit) + implicit none + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + integer,intent(in) :: iunit + double precision,intent(out) :: result + write(iunit, *) & + 'SDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{', expSla,', {',cSla(1),',',cSla(2),',',cSla(3),'} } ],' + result = 0.d0 +end -BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ] +subroutine GauSlaOverlap_read(expGau,cGau,aGau,expSla,cSla,result,iunit) + implicit none + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + integer,intent(in) :: iunit + double precision,intent(out) :: result + read(iunit, *) result +end + +subroutine GauSlaKinetic_write(expGau,cGau,aGau,expSla,cSla,result,iunit) + implicit none + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + integer,intent(in) :: iunit + double precision,intent(out) :: result + write(iunit, *) & + 'TDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{', expSla,',{',cSla(1),',',cSla(2),',',cSla(3),'} } ],' + result = 0.d0 +end + +subroutine GauSlaKinetic_read(expGau,cGau,aGau,expSla,cSla,result,iunit) + implicit none + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + integer,intent(in) :: iunit + double precision,intent(out) :: result + read(iunit, *) result +end + +subroutine GauSlaNuclear_write(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result,iunit) + implicit none + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + double precision,intent(in) :: cNuc(3) + double precision,intent(in) :: ZNuc + integer,intent(in) :: iunit + double precision,intent(out) :: result + write(iunit, *) & + 'VDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{ ', expSla,',{',cSla(1),',',cSla(2),',',cSla(3),'} }, {', ZNuc, ',{', cNuc(1),',', cNuc(2),',', cNuc(3), '} } ],' + result = 0.d0 +end + +subroutine GauSlaNuclear_read(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result,iunit) + implicit none + double precision,intent(in) :: expGau,expSla + double precision,intent(in) :: cGau(3),cSla(3) + integer,intent(in) :: aGau(3) + double precision,intent(in) :: cNuc(3) + double precision,intent(in) :: ZNuc + integer,intent(in) :: iunit + double precision,intent(out) :: result + read(iunit, *) result +end +!TODO + +BEGIN_TEMPLATE + +BEGIN_PROVIDER [ double precision, GauSla$X_matrix, (ao_num, nucl_num) ] implicit none BEGIN_DOC ! overlap matrix @@ -300,6 +451,22 @@ BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ] double precision :: expSla, res, expGau integer :: aGau(3) +!TODO + logical :: read + integer :: iunit + integer :: getunitandopen + + inquire(FILE=trim(ezfio_filename)//'/work/GauSla$X.dat',EXIST=read) + if (read) then + print *, 'READ $X' + iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSla$X.dat','r') + else + print *, 'WRITE $X' + iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSla$X.inp','w') + write(iunit,*) '{' + endif +!TODO + do k=1,nucl_num cSla(1:3) = nucl_coord_transp(1:3,k) expSla = slater_expo(k) @@ -307,28 +474,151 @@ BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ] do i=1,ao_num cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i)) aGau(1:3) = ao_power(i,1:3) - GauSlaOverlap_matrix(i,k) = 0.d0 + GauSla$X_matrix(i,k) = 0.d0 do j=1,ao_prim_num(i) expGau = ao_expo_ordered_transp(j,i) - call GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,res) - GauSlaOverlap_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res +! call GauSla$X(expGau,cGau,aGau,expSla,cSla,res) + if (read) then + call GauSla$X_read(expGau,cGau,aGau,expSla,cSla,res,iunit) + else + call GauSla$X_write(expGau,cGau,aGau,expSla,cSla,res,iunit) + endif + GauSla$X_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res enddo enddo enddo + if (.not.read) then + write(iunit,*) '0.}' + endif + close(iunit) END_PROVIDER -BEGIN_PROVIDER [ double precision, MOSlaOverlap_matrix, (mo_tot_num, nucl_num) ] +BEGIN_PROVIDER [ double precision, MOSla$X_matrix, (mo_tot_num, nucl_num) ] implicit none BEGIN_DOC ! END_DOC call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, & mo_coef_transp, size(mo_coef_transp,1), & - GauSlaOverlap_matrix, size(GauSlaOverlap_matrix,1), & - 0.d0, MOSlaOverlap_matrix, size(MOSlaOverlap_matrix,1)) + GauSla$X_matrix, size(GauSla$X_matrix,1), & + 0.d0, MOSla$X_matrix, size(MOSla$X_matrix,1)) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ] + implicit none + 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), & + GauSla$X_matrix, size(GauSla$X_matrix,1), & + 0.d0, AO_orthoSla$X_matrix, size(AO_orthoSla$X_matrix,1)) + +END_PROVIDER + +SUBST [ X ] + +Overlap ;; +Kinetic ;; + +END_TEMPLATE + +BEGIN_PROVIDER [ double precision, GauSlaNuclear_matrix, (ao_num, nucl_num) ] + implicit none + BEGIN_DOC + ! overlap matrix + END_DOC + integer :: i,j,k,A + double precision :: cGau(3) + double precision :: cSla(3) + double precision :: expSla, res, expGau, Znuc, cNuc(3) + integer :: aGau(3) + +!TODO + logical :: read + integer :: iunit + integer :: getunitandopen + + inquire(FILE=trim(ezfio_filename)//'/work/GauSlaNuclear.dat',EXIST=read) + if (read) then + print *, 'READ Nuclear' + iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSlaNuclear.dat','r') + else + print *, 'WRITE Nuclear' + iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSlaNuclear.inp','w') + write(iunit,*)'{' + endif +!TODO + + do k=1,nucl_num + cSla(1:3) = nucl_coord_transp(1:3,k) + expSla = slater_expo(k) + + do i=1,ao_num + cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i)) + aGau(1:3) = ao_power(i,1:3) + GauSlaNuclear_matrix(i,k) = 0.d0 + + do j=1,ao_prim_num(i) + expGau = ao_expo_ordered_transp(j,i) + do A=1,nucl_num + cNuc(1:3) = nucl_coord_transp(1:3,A) + Znuc = nucl_charge(A) +! call GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res) + if (read) then + call GauSlaNuclear_read(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res,iunit) + else + call GauSlaNuclear_write(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res,iunit) + endif + GauSlaNuclear_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res + enddo + enddo + + enddo + + enddo + if (.not.read) then + write(iunit,*) '0.}' + endif + close(iunit) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, GauSlaH_matrix, (ao_num, nucl_num) ] + implicit none + BEGIN_DOC + ! Core hamiltonian in AO basis + END_DOC + GauSlaH_matrix(1:ao_num,1:nucl_num) = & + GauSlaKinetic_matrix(1:ao_num,1:nucl_num) + & + GauSlaNuclear_matrix(1:ao_num,1:nucl_num) + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, MOSlaNuclear_matrix, (mo_tot_num, nucl_num) ] + implicit none + BEGIN_DOC +! + END_DOC + call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, & + mo_coef_transp, size(mo_coef_transp,1), & + GauSlaNuclear_matrix, size(GauSlaNuclear_matrix,1), & + 0.d0, MOSlaNuclear_matrix, size(MOSlaNuclear_matrix,1)) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, AO_orthoSlaH_matrix, (ao_num, nucl_num) ] + implicit none + 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)) END_PROVIDER diff --git a/src/MOGuess/mo_ortho_lowdin.irp.f b/src/MOGuess/mo_ortho_lowdin.irp.f index 90672b5e..519e4f0d 100644 --- a/src/MOGuess/mo_ortho_lowdin.irp.f +++ b/src/MOGuess/mo_ortho_lowdin.irp.f @@ -6,7 +6,9 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)] ! ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital END_DOC integer :: i,j,k,l - double precision :: tmp_matrix(ao_num_align,ao_num),accu + double precision :: accu + double precision, allocatable :: tmp_matrix(:,:) + allocate (tmp_matrix(ao_num_align,ao_num)) tmp_matrix(:,:) = 0.d0 do j=1, ao_num tmp_matrix(j,j) = 1.d0 @@ -17,6 +19,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)] ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j) enddo enddo + deallocate(tmp_matrix) END_PROVIDER BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)] diff --git a/src/MO_Basis/ao_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f index 48341129..b0eabfbd 100644 --- a/src/MO_Basis/ao_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -82,6 +82,15 @@ END_PROVIDER +BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num,ao_num)] + implicit none + BEGIN_DOC +! ao_ortho_canonical_coef^(-1) + END_DOC + call get_pseudo_inverse(ao_ortho_canonical_coef,ao_num,ao_num, & + ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef,1)) +END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ] implicit none diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 80b45efa..a7416c6b 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -68,6 +68,18 @@ END_PROVIDER endif END_PROVIDER +BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! MO coefficients in orthogonalized AO basis + END_DOC + call dgemm('T','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)) + +END_PROVIDER + BEGIN_PROVIDER [ character*(64), mo_label ] implicit none BEGIN_DOC @@ -257,3 +269,62 @@ subroutine mix_mo_jk(j,k) end +subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) + implicit none + BEGIN_DOC + ! Transform A from the AO basis to the orthogonal AO basis + END_DOC + integer, intent(in) :: LDA_ao,LDA + double precision, intent(in) :: A_ao(LDA_ao,*) + double precision, intent(out) :: A(LDA,*) + double precision, allocatable :: T(:,:) + + allocate ( T(ao_num_align,ao_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + 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','N', ao_num, ao_num, ao_num, 1.d0, & + T, ao_num_align, & + ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),& + 0.d0, A, LDA) + + deallocate(T) +end + +subroutine mo_to_ao_ortho_cano(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the AO orthogonal basis to the AO basis + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + double precision, intent(in) :: A_mo(LDA_mo) + double precision, intent(out) :: A_ao(LDA_ao) + double precision, allocatable :: T(:,:), SC(:,:) + + allocate ( SC(ao_num_align,mo_tot_num) ) + allocate ( T(mo_tot_num_align,ao_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, ao_overlap,size(ao_overlap,1), & + ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), & + 0.d0, SC, ao_num_align) + + call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, & + 1.d0, A_mo,LDA_mo, & + SC, size(SC,1), & + 0.d0, T, mo_tot_num_align) + + call dgemm('N','N', ao_num, ao_num, mo_tot_num, & + 1.d0, SC,size(SC,1), & + T, mo_tot_num_align, & + 0.d0, A_ao, LDA_ao) + + deallocate(T,SC) +end + diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index f9d0094e..51b7b080 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -72,6 +72,8 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) double precision, allocatable :: S_half(:,:) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j +call ortho_lowdin(overlap,LDA,N,C,LDC,m) +return if (n < 2) then return @@ -297,12 +299,12 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA) allocate(work(lwork)) call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info) if (info /= 0) then - print *, info, ': SVD failed' + print *, info, ':: SVD failed' stop 1 endif do i=1,n - if (abs(D(i)) > 1.d-6) then + if (dabs(D(i)) > 1.d-6) then D(i) = 1.d0/D(i) else D(i) = 0.d0