From 23d7d83091d21807e45ff395e458eda4e890d963 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 10 Nov 2021 22:41:40 +0100 Subject: [PATCH] working on ACFDT for T-matrix --- input/methods | 4 +- input/options | 2 +- mol/phenol.xyz | 15 +++++ src/MBPT/G0T0.f90 | 9 ++- src/MBPT/evGT.f90 | 4 +- src/MBPT/qsGT.f90 | 3 +- src/QuAcK/QuAcK.f90 | 2 +- src/RPA/ACFDT_Tmatrix.f90 | 61 ++++++++++---------- src/RPA/RPAx.f90 | 2 +- src/RPA/crRPA.f90 | 32 +++++------ src/RPA/ppRPA.f90 | 118 +++++++++++++++++++++++--------------- 11 files changed, 151 insertions(+), 101 deletions(-) create mode 100644 mol/phenol.xyz diff --git a/input/methods b/input/methods index 6331d53..1e67d49 100644 --- a/input/methods +++ b/input/methods @@ -5,11 +5,11 @@ # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD - T T T T + F F F F # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - T T T T + F F F T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW diff --git a/input/options b/input/options index a4093ba..bbbbbcb 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 256 0.00001 T 5 T 0.00367493 F F F F F # ACFDT: AC Kx XBS - F F T + T T T # BSE: BSE dBSE dTDA evDyn T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/mol/phenol.xyz b/mol/phenol.xyz new file mode 100644 index 0000000..f7871b3 --- /dev/null +++ b/mol/phenol.xyz @@ -0,0 +1,15 @@ +13 +Phenol; structure form HCP92 revised structure; l +C 4.555420 5.661760 4.489060 +C 4.584960 2.843420 4.498910 +C 3.378760 3.548270 4.498890 +C 3.359200 4.943960 4.500070 +C 5.758650 4.948670 4.502800 +C 5.789840 3.550210 4.500330 +H 4.622340 1.760920 4.495600 +H 3.631910 7.321310 4.523190 +H 2.474610 2.963630 4.498170 +H 2.384290 5.419790 4.496670 +H 6.695040 5.492670 4.498160 +H 6.727200 3.021210 4.497160 +O 4.537700 7.024010 4.500450 diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90 index ae946c0..52a080d 100644 --- a/src/MBPT/G0T0.f90 +++ b/src/MBPT/G0T0.f90 @@ -218,6 +218,11 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Free memory + + deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, & + Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t) + ! Compute the BSE correlation energy via the adiabatic connection if(doACFDT) then @@ -234,8 +239,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0T0,EcAC) + call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI_MO,eHF,eG0T0,EcAC) if(exchange_kernel) then diff --git a/src/MBPT/evGT.f90 b/src/MBPT/evGT.f90 index 2309e94..54405e6 100644 --- a/src/MBPT/evGT.f90 +++ b/src/MBPT/evGT.f90 @@ -281,8 +281,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC) + call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI_MO,eGT,eGT,EcAC) if(exchange_kernel) then diff --git a/src/MBPT/qsGT.f90 b/src/MBPT/qsGT.f90 index 1eeff46..1d8e32e 100644 --- a/src/MBPT/qsGT.f90 +++ b/src/MBPT/qsGT.f90 @@ -381,7 +381,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC) + call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI_MO,eGT,eGT,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 4646fb4..7af842b 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -824,7 +824,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_RPA) - call ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + call ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index c41d3d7..256abb9 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -1,5 +1,4 @@ -subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, & +subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & ERI,eT,eGT,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix @@ -27,32 +26,10 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple integer,intent(in) :: nR integer,intent(in) :: nS - integer,intent(in) :: nOOs - integer,intent(in) :: nOOt - integer,intent(in) :: nVVs - integer,intent(in) :: nVVt - double precision,intent(in) :: eT(nBas) double precision,intent(in) :: eGT(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega1s(nVVs) - double precision,intent(in) :: X1s(nVVs,nVVs) - double precision,intent(in) :: Y1s(nOOs,nVVs) - double precision,intent(in) :: Omega2s(nOOs) - double precision,intent(in) :: X2s(nVVs,nOOs) - double precision,intent(in) :: Y2s(nOOs,nOOs) - double precision,intent(in) :: rho1s(nBas,nBas,nVVs) - double precision,intent(in) :: rho2s(nBas,nBas,nOOs) - double precision,intent(in) :: Omega1t(nVVt) - double precision,intent(in) :: X1t(nVVt,nVVt) - double precision,intent(in) :: Y1t(nOOt,nVVt) - double precision,intent(in) :: Omega2t(nOOt) - double precision,intent(in) :: X2t(nVVt,nOOt) - double precision,intent(in) :: Y2t(nOOt,nOOt) - double precision,intent(in) :: rho1t(nBas,nBas,nVVt) - double precision,intent(in) :: rho2t(nBas,nBas,nOOt) - ! Local variables integer :: ispin @@ -62,6 +39,9 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision :: lambda double precision,allocatable :: Ec(:,:) + integer :: nOOs,nOOt + integer :: nVVs,nVVt + double precision :: EcRPA(nspin) double precision,allocatable :: TA(:,:) double precision,allocatable :: TB(:,:) @@ -69,14 +49,37 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) + ! Output variables double precision,intent(out) :: EcAC(nspin) +! Useful quantities + + nOOs = nO*nO + nVVs = nV*nV + + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 + ! Memory allocation - allocate(Ec(nAC,nspin)) + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(Ec(nAC,nspin)) ! Antisymmetrized kernel version @@ -185,8 +188,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 1 iblock = 3 -! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & -! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & + Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) @@ -196,8 +199,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple isp_T = 2 iblock = 4 -! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & -! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & + Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) diff --git a/src/RPA/RPAx.f90 b/src/RPA/RPAx.f90 index a7be80f..f47351d 100644 --- a/src/RPA/RPAx.f90 +++ b/src/RPA/RPAx.f90 @@ -13,8 +13,8 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet - double precision,intent(in) :: eta logical,intent(in) :: triplet + double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRPA.f90 index 38cb004..7645c6c 100644 --- a/src/RPA/crRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -106,25 +106,25 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ! Compute the correlation energy via the adiabatic connection -! if(doACFDT) then + if(doACFDT) then -! write(*,*) '-------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of crRPA correlation energy' -! write(*,*) '-------------------------------------------------------' -! write(*,*) + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of crRPA correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) -! call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) + call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1) -! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if end subroutine crRPA diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index a40e297..9066b10 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) +subroutine ppRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Perform pp-RPA calculation @@ -8,8 +8,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Input variables logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet logical,intent(in) :: triplet + double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -23,16 +26,18 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Local variables integer :: ispin - integer :: nOO - integer :: nVV - double precision,allocatable :: Omega1(:,:) - double precision,allocatable :: X1(:,:,:) - double precision,allocatable :: Y1(:,:,:) - double precision,allocatable :: Omega2(:,:) - double precision,allocatable :: X2(:,:,:) - double precision,allocatable :: Y2(:,:,:) + integer :: nS + integer :: nOOs,nOOt + integer :: nVVs,nVVt + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision :: Ec_ppRPA(nspin) + double precision :: EcAC(nspin) ! Hello world @@ -45,6 +50,25 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Initialization Ec_ppRPA(:) = 0d0 + EcAC(:) = 0d0 + +! Useful quantities + + nS = nO*nV + + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 + + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 + + ! Memory allocation + + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) + + allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) ! Singlet manifold @@ -52,25 +76,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ispin = 1 - ! Useful quantities + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, & + Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin)) - nOO = nO*(nO+1)/2 - nVV = nV*(nV+1)/2 - - ! Memory allocation - - 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,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & - Ec_ppRPA(ispin)) - - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) - - deallocate(Omega1,X1,Y1,Omega2,X2,Y2) + call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s) + call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s) endif @@ -80,26 +90,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ispin = 2 - ! Useful quantities + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, & + Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin)) - nOO = nO*(nO-1)/2 - nVV = nV*(nV-1)/2 - - ! Memory allocation - - 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,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & - Ec_ppRPA(ispin)) - - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) - - deallocate(Omega1,X1,Y1,Omega2,X2,Y2) + call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t) + call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t) endif @@ -112,4 +107,35 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Compute the correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '---------------------------------------------------------' + write(*,*) 'Adiabatic connection version of pp-RPA correlation energy' + write(*,*) '---------------------------------------------------------' + write(*,*) + + call ACFDT_Tmatrix(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI,e,e,EcAC) + + if(exchange_kernel) then + + EcAC(1) = 0.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(1) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end subroutine ppRPA