dipole RHF OK

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
+  C     0.           0.         0.0
-  O     0.           0.         0.89266692
+  O     0.           0.         2.132

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)

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) 

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)   :: S(nBas,nBas)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas),X(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

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

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
+! Input variables
-  double precision,intent(in)        :: eHF(nBas),cHF(nBas,nBas),ENuc,ET,EV,EJ,EK,ERHF
 
-  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 HOMO      energy (eV) = ',eHF(HOMO)*HaToeV
-  write(*,'(A36,F13.6)')  ' HF LUMO      energy (eV) = ',eHF(LUMO)*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-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(*,*)
 

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))
-
+      dipole_int_aa(:,:,:) = dipole_int(:,:,:)
-      call read_dipole_integrals(nBas,dipole_int_aa)
+      dipole_int_bb(:,:,:) = dipole_int(:,:,:)
-      call read_dipole_integrals(nBas,dipole_int_bb)
       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 

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*,'<S**2> = ',2d0*sum(X(:)**2 + Y(:)**2)
 
   end do