diff --git a/examples/molecule.CO b/examples/molecule.CO index 827311b..4254094 100644 --- a/examples/molecule.CO +++ b/examples/molecule.CO @@ -2,4 +2,4 @@ 2 7 7 0 0 # Znuc x y z C 0. 0. 0. - O 0. 0. 2.106 + O 0. 0. 2.134 diff --git a/examples/molecule.F2 b/examples/molecule.F2 index aff52a6..981a397 100644 --- a/examples/molecule.F2 +++ b/examples/molecule.F2 @@ -2,4 +2,4 @@ 2 9 9 0 0 # Znuc x y z F 0. 0. 0. - F 0. 0. 2.6682933 + F 0. 0. 2.622798396593 diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 81c624a..779d849 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 1.4 + H 0. 0. 1.399 diff --git a/examples/molecule.LiH b/examples/molecule.LiH index a1a3c86..635c08d 100644 --- a/examples/molecule.LiH +++ b/examples/molecule.LiH @@ -2,4 +2,4 @@ 2 2 2 0 0 # Znuc x y z Li 0. 0. 0. - H 0. 0. 3.0141132 + H 0. 0. 3.017 diff --git a/input/basis b/input/basis index 20bc6a7..339fa33 100644 --- a/input/basis +++ b/input/basis @@ -1,49 +1,72 @@ -1 9 -S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 -S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 +1 15 +S 9 + 1 6601.0000000 0.0001170 + 2 989.7000000 0.0009110 + 3 225.7000000 0.0047280 + 4 64.2900000 0.0191970 + 5 21.1800000 0.0630470 + 6 7.7240000 0.1632080 + 7 3.0030000 0.3148270 + 8 1.2120000 0.3939360 + 9 0.4930000 0.1969180 +S 9 + 1 6601.0000000 -0.0000180 + 2 989.7000000 -0.0001420 + 3 225.7000000 -0.0007410 + 4 64.2900000 -0.0030200 + 5 21.1800000 -0.0101230 + 6 7.7240000 -0.0270940 + 7 3.0030000 -0.0573590 + 8 1.2120000 -0.0938950 + 9 0.4930000 -0.1210910 S 1 - 1 0.0280500 1.0000000 + 1 0.0951500 1.0000000 S 1 - 1 0.0086400 1.0000000 + 1 0.0479100 1.0000000 +S 1 + 1 0.0222000 1.0000000 P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 + 1 6.2500000 0.0033880 + 2 1.3700000 0.0193160 + 3 0.3672000 0.0791040 P 1 - 1 0.0240300 1.0000000 + 1 0.1192000 1.0000000 P 1 - 1 0.0057900 1.0000000 + 1 0.0447400 1.0000000 +P 1 + 1 0.0179500 1.0000000 D 1 - 1 0.1239000 1.0000000 + 1 0.3440000 1.0000000 D 1 - 1 0.0725000 1.0000000 -2 5 + 1 0.1530000 1.0000000 +D 1 + 1 0.0680000 1.0000000 +F 1 + 1 0.2460000 1.0000000 +F 1 + 1 0.1292000 1.0000000 +G 1 + 1 0.2380000 1.0000000 +2 10 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 - + 1 0.8380000 1.0000000 +P 1 + 1 0.2920000 1.0000000 +D 1 + 1 2.0620000 1.0000000 +D 1 + 1 0.6620000 1.0000000 +F 1 + 1 1.3970000 1.0000000 diff --git a/input/methods b/input/methods index 8e05978..2aa67c1 100644 --- a/input/methods +++ b/input/methods @@ -7,10 +7,10 @@ # CIS RPA RPAx ppRPA ADC F F F F F # GF2 GF3 - T F + F F # G0W0 evGW qsGW - F F F -# G0T0 evGT qsGT T F F +# G0T0 evGT qsGT + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index a1a3c86..635c08d 100644 --- a/input/molecule +++ b/input/molecule @@ -2,4 +2,4 @@ 2 2 2 0 0 # Znuc x y z Li 0. 0. 0. - H 0. 0. 3.0141132 + H 0. 0. 3.017 diff --git a/input/molecule.xyz b/input/molecule.xyz index dfbf65a..1982c68 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 Li 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 1.5950001314 + H 0.0000000000 0.0000000000 1.5965277602 diff --git a/input/options b/input/options index 5a8d22a..ecc7c5f 100644 --- a/input/options +++ b/input/options @@ -7,10 +7,10 @@ # CIS/TDHF/BSE: singlet triplet T T # GF: maxSCF thresh DIIS n_diis renormalization - 64 0.00001 T 5 3 + 1 0.00001 T 5 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta 256 0.00001 T 5 F F T F F F F 0.000 # ACFDT: AC Kx XBS - T F T + T F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/input/weight b/input/weight index 20bc6a7..339fa33 100644 --- a/input/weight +++ b/input/weight @@ -1,49 +1,72 @@ -1 9 -S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 -S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 +1 15 +S 9 + 1 6601.0000000 0.0001170 + 2 989.7000000 0.0009110 + 3 225.7000000 0.0047280 + 4 64.2900000 0.0191970 + 5 21.1800000 0.0630470 + 6 7.7240000 0.1632080 + 7 3.0030000 0.3148270 + 8 1.2120000 0.3939360 + 9 0.4930000 0.1969180 +S 9 + 1 6601.0000000 -0.0000180 + 2 989.7000000 -0.0001420 + 3 225.7000000 -0.0007410 + 4 64.2900000 -0.0030200 + 5 21.1800000 -0.0101230 + 6 7.7240000 -0.0270940 + 7 3.0030000 -0.0573590 + 8 1.2120000 -0.0938950 + 9 0.4930000 -0.1210910 S 1 - 1 0.0280500 1.0000000 + 1 0.0951500 1.0000000 S 1 - 1 0.0086400 1.0000000 + 1 0.0479100 1.0000000 +S 1 + 1 0.0222000 1.0000000 P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 + 1 6.2500000 0.0033880 + 2 1.3700000 0.0193160 + 3 0.3672000 0.0791040 P 1 - 1 0.0240300 1.0000000 + 1 0.1192000 1.0000000 P 1 - 1 0.0057900 1.0000000 + 1 0.0447400 1.0000000 +P 1 + 1 0.0179500 1.0000000 D 1 - 1 0.1239000 1.0000000 + 1 0.3440000 1.0000000 D 1 - 1 0.0725000 1.0000000 -2 5 + 1 0.1530000 1.0000000 +D 1 + 1 0.0680000 1.0000000 +F 1 + 1 0.2460000 1.0000000 +F 1 + 1 0.1292000 1.0000000 +G 1 + 1 0.2380000 1.0000000 +2 10 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 82.6400000 0.0020060 + 2 12.4100000 0.0153430 + 3 2.8240000 0.0755790 S 1 - 1 0.1220000 1.0000000 + 1 0.7977000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.2581000 1.0000000 +S 1 + 1 0.0898900 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 2.2920000 1.0000000 P 1 - 1 0.1410000 1.0000000 - + 1 0.8380000 1.0000000 +P 1 + 1 0.2920000 1.0000000 +D 1 + 1 2.0620000 1.0000000 +D 1 + 1 0.6620000 1.0000000 +F 1 + 1 1.3970000 1.0000000 diff --git a/src/QuAcK/ACFDT.f90 b/src/QuAcK/ACFDT.f90 index 60e021b..1e457bc 100644 --- a/src/QuAcK/ACFDT.f90 +++ b/src/QuAcK/ACFDT.f90 @@ -1,5 +1,5 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC) + nBas,nC,nO,nV,nR,nS,ERI,eHF,e,Omega,XpY,XmY,rho,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem @@ -19,6 +19,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -76,16 +77,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man if(doXBS) then - call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & + call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, & rho(:,:,:,ispin),EcAC(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,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & + ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -121,7 +123,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man if(doXBS) then - call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & + call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, & rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) @@ -130,7 +132,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & + ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix.f90 index 2fb580c..bc0b9d7 100644 --- a/src/QuAcK/Bethe_Salpeter_A_matrix.f90 +++ b/src/QuAcK/Bethe_Salpeter_A_matrix.f90 @@ -39,6 +39,7 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps enddo +! A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,a,j,b) + 4d0*lambda*chi A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,a,j,b) + 4d0*lambda*chi enddo @@ -46,7 +47,7 @@ subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, enddo enddo - print*,'BSE A' - call matout(nS,nS,A_lr) +! print*,'BSE A' +! call matout(nS,nS,A_lr) end subroutine Bethe_Salpeter_A_matrix diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 index cf7573b..3ecefb6 100644 --- a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 @@ -39,6 +39,7 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps enddo +! B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,a,b,j) + 4d0*lambda*chi B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,a,b,j) + 4d0*lambda*chi enddo @@ -46,7 +47,7 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, enddo enddo +! print*,'BSE B' +! call matout(nS,nS,B_lr) - print*,'BSE B' - call matout(nS,nS,B_lr) end subroutine Bethe_Salpeter_B_matrix diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 5c5a9db..159f3a4 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -1,4 +1,4 @@ -subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF,eG0T0) +subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) ! Perform G0W0 calculation with a T-matrix self-energy (G0T0) @@ -7,9 +7,6 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc, ! Input variables - logical,intent(in) :: BSE - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR @@ -36,9 +33,9 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc, double precision,allocatable :: SigT(:) double precision,allocatable :: Z(:) -! Output variables + double precision,allocatable :: eG0T0(:) - double precision,intent(out) :: eG0T0(nBas) +! Output variables ! Hello world @@ -64,7 +61,7 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc, Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & rho1t(nBas,nO-nC,nVVt),rho2t(nBas,nV-nR,nOOt), & - SigT(nBas),Z(nBas)) + SigT(nBas),Z(nBas),eG0T0(nBas)) !---------------------------------------------- ! Singlet manifold diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index f11de8d..75d35ff 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -166,7 +166,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif end if call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC) + nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC) if(exchange_kernel) then diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 4b8c0ee..f4eaa92 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -597,8 +597,7 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) -! call G0T0(BSE,singlet_manifold,triplet_manifold,eta, & -! nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0) + ! call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF) call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_G0T0) diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 77e7a1a..f1e28fa 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -257,7 +257,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE end if call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC) + nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/linear_response_A_matrix.f90 b/src/QuAcK/linear_response_A_matrix.f90 index 1a1d391..79272a0 100644 --- a/src/QuAcK/linear_response_A_matrix.f90 +++ b/src/QuAcK/linear_response_A_matrix.f90 @@ -47,8 +47,10 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, jb = jb + 1 A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - + (1d0 + delta_spin)*lambda*ERI(i,b,a,j) & - - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) +! + (1d0 + delta_spin)*lambda*ERI(i,b,a,j) & +! - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) + + (1d0 + delta_spin)*lambda*ERI(i,j,a,b) & + - (1d0 - delta_dRPA)*lambda*ERI(i,a,j,b) enddo enddo diff --git a/src/QuAcK/linear_response_B_matrix.f90 b/src/QuAcK/linear_response_B_matrix.f90 index 5226059..335ab1f 100644 --- a/src/QuAcK/linear_response_B_matrix.f90 +++ b/src/QuAcK/linear_response_B_matrix.f90 @@ -44,8 +44,10 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_ do b=nO+1,nBas-nR jb = jb + 1 - B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) & - - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a) +! B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) & +! - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a) + B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,b,a,j) & + - (1d0 - delta_dRPA)*lambda*ERI(i,a,b,j) enddo enddo diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 93be12e..8321cda 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -17,6 +17,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1 ! Local variables integer :: ab,cd,ij,kl + integer :: p,q,r,s double precision :: trace_matrix double precision,allocatable :: B(:,:) double precision,allocatable :: C(:,:) @@ -64,31 +65,46 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1 M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO)) - open(unit=42,file='B.dat') - open(unit=43,file='C.dat') - open(unit=44,file='D.dat') + open(unit=42,file='B.dat') + open(unit=43,file='C.dat') + open(unit=44,file='D.dat') + open(unit=45,file='ERI.dat') + open(unit=46,file='eps.dat') + + do ab=1,nVV + do ij=1,nOO + if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij) + end do + end do + + do ab=1,nVV + do cd=1,nVV + if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd) + end do + end do + + do ij=1,nOO + do kl=1,nOO + if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl) + end do + end do + + do p=1,nBas + write(46,*) p,e(p) + do q=1,nBas + do r=1,nBas + do s=1,nBas + if(abs(ERI(p,q,r,s)) > 1d-15) write(45,*) p,q,r,s,ERI(p,q,r,s) + end do + end do + end do + end do - do ab=1,nVV - do ij=1,nOO - if(abs(B(ab,ij)) > 1d-15) write(42,*) ab,ij,B(ab,ij) - end do - end do - - do ab=1,nVV - do cd=1,nVV - if(abs(C(ab,cd)) > 1d-15) write(43,*) ab,cd,C(ab,cd) - end do - end do - - do ij=1,nOO - do kl=1,nOO - if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl) - end do - end do - - close(42) - close(43) - close(44) + close(42) + close(43) + close(44) + close(45) + close(46) ! print*, 'pp-RPA matrix' diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 69fb5b5..07c54dd 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -288,7 +288,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, end if call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho,EcAC) + nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index 7063d7b..90ba9c0 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/src/QuAcK/sort_ppRPA.f90 @@ -101,7 +101,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! Orthogonalize eigenvectors - S1 = matmul(transpose(Z1),matmul(M,Z1)) + S1 = + matmul(transpose(Z1),matmul(M,Z1)) S2 = - matmul(transpose(Z2),matmul(M,Z2)) call orthogonalization_matrix(1,nVV,S1,O1) @@ -116,12 +116,17 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) X2(1:nVV,1:nOO) = Z2( 1: nVV,1:nOO) Y2(1:nOO,1:nOO) = Z2(nVV+1:nOO+nVV,1:nOO) - write(*,*) 'X1t.X1 - Y1t.Y1' - call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1)) - write(*,*) 'X2t.X2 - Y2t.Y2' - call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2)) - write(*,*) 'X1t.X2 - Y1t.Y2' - call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2)) +! write(*,*) 'Z1t.M.Z1' +! call matout(nVV,nVV,matmul(matmul(transpose(Z1),M),Z1)) +! write(*,*) 'Z2t.M.Z2' +! call matout(nOO,nOO,matmul(matmul(transpose(Z2),M),Z2)) + +! write(*,*) 'X1t.X1 - Y1t.Y1' +! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1)) +! write(*,*) 'X2t.X2 - Y2t.Y2' +! call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2)) +! write(*,*) 'X1t.X2 - Y1t.Y2' +! call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2)) end subroutine sort_ppRPA