rename TDHF and co

src/QuAcK/ACDFT.f90

@@ -1,127 +0,0 @@
-subroutine ACDFT(scaled_screening,dRPA,TDA,BSE,singlet_manifold,triplet_manifold, & 
-                 nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho)
-
-! Compute the correlation energy via the adiabatic connection dissipation fluctuation theorem
-
-  implicit none
-  include 'parameters.h'
-  include 'quadrature.h'
-
-! Input variables
-
-  logical,intent(in)            :: scaled_screening
-  logical,intent(in)            :: dRPA
-  logical,intent(in)            :: TDA
-  logical,intent(in)            :: BSE
-  logical,intent(in)            :: singlet_manifold
-  logical,intent(in)            :: triplet_manifold
-
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-
-  double precision              :: Omega(nS,nspin)
-  double precision              :: XpY(nS,nS,nspin)
-  double precision              :: XmY(nS,nS,nspin)
-  double precision              :: rho(nBas,nBas,nS,nspin)
-
-! Local variables
-
-  integer                       :: ispin
-  logical                       :: adiabatic_connection
-  integer                       :: iAC
-  double precision              :: lambda
-  double precision,allocatable  :: Ec(:,:)
-  double precision,allocatable  :: EcAC(:,:)
-
-! Memory allocation
-
-  allocate(Ec(nAC,nspin),EcAC(nAC,nspin))
-
-  if(singlet_manifold) then
-
-    ispin = 1
-    Ec(:,ispin)   = 0d0
-    EcAC(:,ispin) = 0d0
-
-    write(*,*) '--------------'
-    write(*,*) 'Singlet states'
-    write(*,*) '--------------'
-    write(*,*) 
-
-    write(*,*) '-----------------------------------------------------------------------------------'
-    write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(V x P_lambda)'
-    write(*,*) '-----------------------------------------------------------------------------------'
-
-    do iAC=1,nAC
- 
-      lambda = rAC(iAC)
-
-      if(scaled_screening) then
-
-        call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
-                             rho(:,:,:,ispin),Ec(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
-        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
-
-      end if
-  
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
-                       rho(:,:,:,ispin),Ec(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
-
-      call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
-
-      write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,Ec(iAC,ispin),EcAC(iAC,ispin)
-
-    end do
-
-    write(*,*) '-----------------------------------------------------------------------------------'
-    write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
-    write(*,*) '-----------------------------------------------------------------------------------'
-    write(*,*)
-
-  end if
- 
-  if(triplet_manifold) then
-
-    ispin = 2
-    Ec(:,ispin)   = 0d0
-    EcAC(:,ispin) = 0d0
-
-    write(*,*) '--------------'
-    write(*,*) 'Triplet states'
-    write(*,*) '--------------'
-    write(*,*) 
-
-    write(*,*) '-----------------------------------------------------------------------------------'
-    write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(V x P_lambda)'
-    write(*,*) '-----------------------------------------------------------------------------------'
-
-    do iAC=1,nAC
- 
-      lambda = rAC(iAC)
-
-      if(scaled_screening) then
-
-        call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
-                             rho(:,:,:,ispin),Ec(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
-        call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
-
-      end if  
-
-      call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
-                           rho(:,:,:,ispin),Ec(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
-
-      call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
-
-      write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,Ec(iAC,ispin),EcAC(iAC,ispin)
-
-    end do
-
-    write(*,*) '-----------------------------------------------------------------------------------'
-    write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin))
-    write(*,*) '-----------------------------------------------------------------------------------'
-    write(*,*)
-
-  end if
-
-end subroutine ACDFT

src/QuAcK/Ec_AC.f90

@@ -1,92 +0,0 @@
-subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC)
-
-! Compute the correlation energy via the adiabatic connection formula
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  integer,intent(in)            :: ispin
-  logical,intent(in)            :: dRPA
-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: XpY(nS,nS)
-  double precision,intent(in)   :: XmY(nS,nS)
-
-! Local variables
-
-  integer                       :: i,j,a,b
-  integer                       :: ia,jb,kc
-  double precision              :: delta_spin
-  double precision              :: delta_dRPA
-  double precision,allocatable  :: P(:,:)
-  double precision,allocatable  :: Ap(:,:)
-  double precision,allocatable  :: Bp(:,:)
-  double precision,allocatable  :: X(:,:)
-  double precision,allocatable  :: Y(:,:)
-  double precision,external     :: trace_matrix
-
-! Output variables
-
-  double precision,intent(out)  :: EcAC
-
-! Singlet or triplet manifold?
-
-  delta_spin = 0d0
-  if(ispin == 1) delta_spin = +1d0
-  if(ispin == 2) delta_spin = -1d0
-
-! Direct RPA
-
-  delta_dRPA = 0d0
-  if(dRPA) delta_dRPA = 1d0
-
-! Memory allocation
-
-  allocate(P(nS,nS),Ap(nS,nS),Bp(nS,nS),X(nS,nS),Y(nS,nS))
-
-! Compute P = (X+Y)(X+Y) - 1
-
-  P(:,:) = matmul(transpose(XpY),XpY)
-
-  do ia=1,nS
-    P(ia,ia) = P(ia,ia) - 1d0
-  enddo
-
-! Compute Aiajb = (ia|bj) and Biajb = (ia|jb)
-
-  ia = 0
-  do i=nC+1,nO
-    do a=nO+1,nBas-nR
-      ia = ia + 1
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
-
-            Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) 
-            Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,b,a) 
-
-        enddo
-      enddo
-    enddo
-  enddo
-
-! Compute Tr(A x P)
-
-! EcAC = trace_matrix(nS,matmul(Ap,P))
-
-! print*,'EcAC =',EcAC
-
-  X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:))
-  Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:))
-
-  EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) &
-       + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) &
-       - trace_matrix(nS,Ap)
-
-! print*,'EcAC =',EcAC
-
-end subroutine Ec_AC
-

src/QuAcK/G0W0.f90

@@ -173,7 +173,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
 
       end if
 
-      call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & 
+      call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & 
                  nBas,nC,nO,nV,nR,nS,ERI,eG0W0,Omega,XpY,XmY,rho)
 
     end if

src/QuAcK/QuAcK.f90

@@ -8,7 +8,7 @@ program QuAcK
   logical                       :: doMP2,doMP3,doMP2F12
   logical                       :: doCCD,doCCSD,doCCSDT
   logical                       :: do_ring_CCD,do_ladder_CCD
-  logical                       :: doCIS,doRPA,doTDHF
+  logical                       :: doCIS,doRPA,doRPAx
   logical                       :: doppRPA,doADC
   logical                       :: doGF2,doGF3
   logical                       :: doG0W0,doevGW,doqsGW
@@ -47,7 +47,7 @@ program QuAcK
   double precision              :: start_CCSD   ,end_CCSD     ,t_CCSD
   double precision              :: start_CIS    ,end_CIS      ,t_CIS
   double precision              :: start_RPA    ,end_RPA      ,t_RPA 
-  double precision              :: start_TDHF   ,end_TDHF     ,t_TDHF
+  double precision              :: start_RPAx   ,end_RPAx     ,t_RPAx
   double precision              :: start_ppRPA  ,end_ppRPA    ,t_ppRPA
   double precision              :: start_ADC    ,end_ADC      ,t_ADC
   double precision              :: start_GF2    ,end_GF2      ,t_GF2
@@ -113,7 +113,7 @@ program QuAcK
                     doMP2,doMP3,doMP2F12,      &
                     doCCD,doCCSD,doCCSDT,      &
                     do_ring_CCD,do_ladder_CCD, &
-                    doCIS,doRPA,doTDHF,        & 
+                    doCIS,doRPA,doRPAx,        & 
                     doppRPA,doADC,             &
                     doGF2,doGF3,               &
                     doG0W0,doevGW,doqsGW,      &
@@ -444,17 +444,17 @@ program QuAcK
   end if
 
 !------------------------------------------------------------------------
-! Compute TDHF excitations
+! Compute RPAx excitations
 !------------------------------------------------------------------------
 
-  if(doTDHF) then
+  if(doRPAx) then
 
-    call cpu_time(start_TDHF)
+    call cpu_time(start_RPAx)
-    call TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
+    call RPAx(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
-    call cpu_time(end_TDHF)
+    call cpu_time(end_RPAx)
 
-    t_TDHF = end_TDHF - start_TDHF
+    t_RPAx = end_RPAx - start_RPAx
-    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for TDHF = ',t_TDHF,' seconds'
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPAx,' seconds'
     write(*,*)
 
   end if

src/QuAcK/RPA.f90

@@ -93,7 +93,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
     write(*,*) '------------------------------------------------------'
     write(*,*) 
  
-    call ACDFT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, &
+    call ACFDT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, &
                nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho)
 
   end if

src/QuAcK/TDHF.f90

@@ -1,101 +0,0 @@
-subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
-
-! Perform random phase approximation calculation
-
-  implicit none
-  include 'parameters.h'
-  include 'quadrature.h'
-
-! Input variables
-
-  logical,intent(in)            :: singlet_manifold
-  logical,intent(in)            :: triplet_manifold
-  integer,intent(in)            :: nBas
-  integer,intent(in)            :: nC
-  integer,intent(in)            :: nO
-  integer,intent(in)            :: nV
-  integer,intent(in)            :: nR
-  integer,intent(in)            :: nS
-  double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
-  double precision,intent(in)   :: e(nBas)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-
-! Local variables
-
-  integer                       :: ispin
-  double precision,allocatable  :: Omega(:,:)
-  double precision,allocatable  :: XpY(:,:,:)
-  double precision,allocatable  :: XmY(:,:,:)
-
-  double precision              :: rho
-  double precision              :: EcRPA(nspin)
-
-  logical                       :: adiabatic_connection
-
-! Hello world
-
-  write(*,*)
-  write(*,*)'************************************************'
-  write(*,*)'|  Time-dependent Hartree-Fock calculation     |'
-  write(*,*)'************************************************'
-  write(*,*)
-
-! Initialization
-
-  EcRPA(:) = 0d0
-
-! Memory allocation
-
-  allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
-
-! Singlet manifold
-
-  if(singlet_manifold) then 
-
-    ispin = 1
-
-    call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
-                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
-    call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
-
-  endif
-
-! Triplet manifold 
-
-  if(triplet_manifold) then 
-
-    ispin = 2
-
-    call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
-                         EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
-    call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
-
-  endif
-
-  write(*,*)
-  write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (singlet) =',EcRPA(1)
-  write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (triplet) =',EcRPA(2)
-  write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy           =',EcRPA(1) + EcRPA(2)
-  write(*,'(2X,A40,F15.6)') 'Tr@TDHF total energy                 =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
-  write(*,*)'-------------------------------------------------------------------------------'
-  write(*,*)
-
-! Compute the correlation energy via the adiabatic connection 
-
-  adiabatic_connection = .true.
-
-  if(adiabatic_connection) then
-
-    write(*,*) '-------------------------------------------------------'
-    write(*,*) 'Adiabatic connection version of TDHF correlation energy'
-    write(*,*) '-------------------------------------------------------'
-    write(*,*)
-
-    call ACDFT(.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold, &
-               nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho)
-
-  end if
-
-end subroutine TDHF

src/QuAcK/evGW.f90

@@ -279,7 +279,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
 
       end if
 
-      call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
+      call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
                  nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho)
 
     end if

src/QuAcK/ppRPA.f90

@@ -21,7 +21,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
 
 ! Local variables
 
-  logical                       :: BSE
   integer                       :: ispin
   integer                       :: nOO
   integer                       :: nVV
@@ -46,10 +45,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
 
   Ec_ppRPA(:) = 0d0
 
-! Switch off Bethe-Salpeter equation for TDHF
-
-  BSE = .false. 
-
 ! Singlet manifold
 
   if(singlet_manifold) then 
@@ -66,9 +61,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
     allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & 
              Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
 
-    call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,    & 
+    call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & 
-                            Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & 
+                            Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),  & 
-                            Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & 
+                            Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),  & 
                             Ec_ppRPA(ispin))
 
     call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
@@ -95,9 +90,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
              Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
 
 
-    call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,    &
+    call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
-                            Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & 
+                            Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin),  & 
-                            Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & 
+                            Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin),  & 
                             Ec_ppRPA(ispin))
 
     call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
@@ -109,10 +104,10 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
-  write(*,'(2X,A40,F15.6)') 'pp-RPA   correlation energy (singlet) =',Ec_ppRPA(1)
+  write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy (singlet) =',Ec_ppRPA(1)
-  write(*,'(2X,A40,F15.6)') 'pp-RPA   correlation energy (triplet) =',3d0*Ec_ppRPA(2)
+  write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy (triplet) =',3d0*Ec_ppRPA(2)
-  write(*,'(2X,A40,F15.6)') 'pp-RPA   correlation energy           =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
+  write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy           =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
-  write(*,'(2X,A40,F15.6)') 'pp-RPA   total energy                 =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
+  write(*,'(2X,A40,F15.6)') 'Tr@ppRPA total energy                 =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
   write(*,*)'-------------------------------------------------------------------------------'
   write(*,*)
 

src/QuAcK/qsGW.f90

@@ -311,7 +311,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
 
       end if
 
-      call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
+      call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
                  nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho)
 
     end if

src/QuAcK/read_methods.f90

@@ -2,7 +2,7 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
                         doMP2,doMP3,doMP2F12,      & 
                         doCCD,doCCSD,doCCSDT,      & 
                         do_ring_CCD,do_ladder_CCD, &
-                        doCIS,doRPA,doTDHF,        & 
+                        doCIS,doRPA,doRPAx,        & 
                         doppRPA,doADC,             & 
                         doGF2,doGF3,               & 
                         doG0W0,doevGW,doqsGW,      & 
@@ -19,7 +19,7 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
   logical,intent(out)           :: doMP2,doMP3,doMP2F12
   logical,intent(out)           :: doCCD,doCCSD,doCCSDT
   logical,intent(out)           :: do_ring_CCD,do_ladder_CCD
-  logical,intent(out)           :: doCIS,doRPA,doTDHF,doppRPA,doADC
+  logical,intent(out)           :: doCIS,doRPA,doRPAx,doppRPA,doADC
   logical,intent(out)           :: doGF2,doGF3  
   logical,intent(out)           :: doG0W0,doevGW,doqsGW
   logical,intent(out)           :: doG0T0,doevGT,doqsGT
@@ -52,7 +52,7 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
 
   doCIS   = .false.
   doRPA   = .false.
-  doTDHF  = .false.
+  doRPAx  = .false.
   doppRPA = .false.
   doADC   = .false.
 
@@ -101,7 +101,7 @@ subroutine read_methods(doRHF,doUHF,doMOM,         &
   read(1,*) answer1,answer2,answer3,answer4,answer5
   if(answer1 == 'T') doCIS   = .true.
   if(answer2 == 'T') doRPA   = .true.
-  if(answer3 == 'T') doTDHF  = .true.
+  if(answer3 == 'T') doRPAx  = .true.
   if(answer4 == 'T') doppRPA = .true.
   if(answer5 == 'T') doADC   = .true.