From 5a3c2245135cc0dbd606d94d4e726eda03e954a3 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 4 Oct 2020 09:36:03 +0200 Subject: [PATCH] dipole RHF OK --- examples/molecule.CO | 4 +-- include/parameters.h | 1 + input/methods | 2 +- src/HF/RHF.f90 | 22 +++++++++--- src/HF/dipole_moment.f90 | 53 ++++++++++++++++++++++++++++ src/HF/print_RHF.f90 | 33 +++++++++++++---- src/QuAcK/QuAcK.f90 | 14 ++++---- src/RPA/print_transition_vectors.f90 | 2 +- 8 files changed, 108 insertions(+), 23 deletions(-) create mode 100644 src/HF/dipole_moment.f90 diff --git a/examples/molecule.CO b/examples/molecule.CO index 60902b6..f65ce7a 100644 --- a/examples/molecule.CO +++ b/examples/molecule.CO @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd 2 7 7 0 0 # Znuc x y z - C 0. 0. -1.24942055 - O 0. 0. 0.89266692 + C 0. 0. 0.0 + O 0. 0. 2.132 diff --git a/include/parameters.h b/include/parameters.h index 279fe1e..d1d4fcb 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -15,6 +15,7 @@ double precision,parameter :: HaToeV = 27.21138602d0 double precision,parameter :: pmtoau = 0.0188973d0 double precision,parameter :: BoToAn = 0.529177249d0 + double precision,parameter :: auToD = 2.5415802529d0 double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0) double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) diff --git a/input/methods b/input/methods index 3a1c5a4..9a5f0de 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - F T F + T F F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 575324b..7770314 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,18 +1,28 @@ -subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,e,c,P) +subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,Enuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P) ! Perform restricted Hartree-Fock calculation implicit none + include 'parameters.h' ! Input variables integer,intent(in) :: maxSCF,max_diis,guess_type double precision,intent(in) :: thresh - integer,intent(in) :: nBas,nO + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - double precision,intent(in) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),X(nBas,nBas) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -33,6 +43,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) double precision,allocatable :: ON(:) + double precision :: dipole(ncart) ! Output variables @@ -181,6 +192,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) ERHF = ET + EV + EJ + EK - call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF) + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) end subroutine RHF diff --git a/src/HF/dipole_moment.f90 b/src/HF/dipole_moment.f90 new file mode 100644 index 0000000..49b11f8 --- /dev/null +++ b/src/HF/dipole_moment.f90 @@ -0,0 +1,53 @@ +subroutine dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) + +! Compute density matrix based on the occupation numbers + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nNuc + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ixyz + integer :: iNuc + integer :: mu,nu + +! Output variables + + double precision,intent(out) :: dipole(ncart) + +! Initialization + + dipole(:) = 0d0 + +! Loop over cartesian components + + do ixyz=1,ncart + + ! Nuclear part + + do iNuc=1,nNuc + + dipole(ixyz) = dipole(ixyz) + ZNuc(iNuc)*rNuc(iNuc,ixyz) + + end do + + ! Electronic part + + do mu=1,nBas + do nu=1,nBas + dipole(ixyz) = dipole(ixyz) - P(mu,nu)*dipole_int(mu,nu,ixyz) + enddo + enddo + + enddo + +end subroutine dipole_moment diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 0cea950..56bd1af 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -1,14 +1,29 @@ -subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF) +subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) ! Print one-electron energies and other stuff for G0W0 implicit none include 'parameters.h' - integer,intent(in) :: nBas,nO - double precision,intent(in) :: eHF(nBas),cHF(nBas,nBas),ENuc,ET,EV,EJ,EK,ERHF +! Input variables - integer :: HOMO,LUMO + integer,intent(in) :: nBas + integer,intent(in) :: nO + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ET + double precision,intent(in) :: EV + double precision,intent(in) :: EJ + double precision,intent(in) :: EK + double precision,intent(in) :: ERHF + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: ixyz + integer :: HOMO + integer :: LUMO double precision :: Gap ! HOMO and LUMO @@ -36,9 +51,13 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF) write(*,'(A32,1X,F16.10)') ' Nuclear repulsion = ',ENuc write(*,'(A32,1X,F16.10)') ' Hartree-Fock energy = ',ERHF + ENuc write(*,'(A50)') '-----------------------------------------' - write(*,'(A36,F13.6)') ' HF HOMO energy (eV) = ',eHF(HOMO)*HaToeV - write(*,'(A36,F13.6)') ' HF LUMO energy (eV) = ',eHF(LUMO)*HaToeV - write(*,'(A36,F13.6)') ' HF HOMO-LUMO gap (eV) = ',Gap*HaToeV + write(*,'(A36,F13.6)') ' HF HOMO energy (eV) = ',eHF(HOMO)*HaToeV + write(*,'(A36,F13.6)') ' HF LUMO energy (eV) = ',eHF(LUMO)*HaToeV + write(*,'(A36,F13.6)') ' HF HOMO-LUMO gap (eV) = ',Gap*HaToeV + write(*,'(A50)') '-----------------------------------------' + write(*,'(A36)') ' Dipole moment (Debye) ' + write(*,'(4A10)') 'X','Y','Z','Tot.' + write(*,'(4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD write(*,'(A50)') '-----------------------------------------' write(*,*) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ed41f7a..c84671d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -235,7 +235,8 @@ program QuAcK ! Memory allocation for one- and two-electron integrals allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), & - S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas)) + S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), & + dipole_int(nBas,nBas,ncart)) ! Read integrals @@ -249,6 +250,7 @@ program QuAcK call system('./GoQCaml') call read_integrals(nBas,S,T,V,Hc,ERI_AO) + call read_dipole_integrals(nBas,dipole_int) end if @@ -270,7 +272,8 @@ program QuAcK if(doRHF) then call cpu_time(start_HF) - call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,ERHF,eHF,cHF,PHF) + call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int,X,ERHF,eHF,cHF,PHF) call cpu_time(end_HF) t_HF = end_HF - start_HF @@ -343,9 +346,8 @@ program QuAcK ! Read and transform dipole-related integrals allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart)) - - call read_dipole_integrals(nBas,dipole_int_aa) - call read_dipole_integrals(nBas,dipole_int_bb) + dipole_int_aa(:,:,:) = dipole_int(:,:,:) + dipole_int_bb(:,:,:) = dipole_int(:,:,:) do ixyz=1,ncart call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) @@ -402,8 +404,6 @@ program QuAcK ! Read and transform dipole-related integrals - allocate(dipole_int(nBas,nBas,ncart)) - call read_dipole_integrals(nBas,dipole_int) do ixyz=1,ncart call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz)) end do diff --git a/src/RPA/print_transition_vectors.f90 b/src/RPA/print_transition_vectors.f90 index b7f6e13..ef29043 100644 --- a/src/RPA/print_transition_vectors.f90 +++ b/src/RPA/print_transition_vectors.f90 @@ -108,7 +108,7 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, end do write(*,*) - print*,2d0*sum(X(:)**2 + Y(:)**2) + print*,' = ',2d0*sum(X(:)**2 + Y(:)**2) end do