diff --git a/examples/basis.Hu b/examples/basis.Hu index 6f77b3e..0f3f20d 100644 --- a/examples/basis.Hu +++ b/examples/basis.Hu @@ -1,6 +1,6 @@ 1 1 S 1 1.00 - 1.0000000 1.0000000 +1 1.0000000 1.0000000 2 1 S 1 1.00 - 1.0000000 1.0000000 +1 1.0000000 1.0000000 diff --git a/examples/molecule.F2 b/examples/molecule.F2 index af08593..aff52a6 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.64 + F 0. 0. 2.6682933 diff --git a/examples/molecule.LiH b/examples/molecule.LiH index 1d6bb3f..a1a3c86 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. 2.994 + H 0. 0. 3.0141132 diff --git a/input/basis b/input/basis index 69b24c8..20bc6a7 100644 --- a/input/basis +++ b/input/basis @@ -1,30 +1,49 @@ -1 6 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 +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 S 1 - 1 0.3258000 1.0000000 + 1 0.0280500 1.0000000 S 1 - 1 0.1027000 1.0000000 + 1 0.0086400 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 P 1 - 1 1.4070000 1.0000000 + 1 0.0240300 1.0000000 P 1 - 1 0.3880000 1.0000000 + 1 0.0057900 1.0000000 D 1 - 1 1.0570000 1.0000000 -2 6 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 -S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 + 1 0.1239000 1.0000000 D 1 - 1 1.0570000 1.0000000 + 1 0.0725000 1.0000000 +2 5 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 +S 1 + 1 0.1220000 1.0000000 +S 1 + 1 0.0297400 1.0000000 +P 1 + 1 0.7270000 1.0000000 +P 1 + 1 0.1410000 1.0000000 + diff --git a/input/methods b/input/methods index 2aa67c1..8e05978 100644 --- a/input/methods +++ b/input/methods @@ -7,10 +7,10 @@ # CIS RPA RPAx ppRPA ADC F F F F F # GF2 GF3 - F F + T F # G0W0 evGW qsGW - T F F -# G0T0 evGT qsGT F F F +# G0T0 evGT qsGT + T F F # MCMP2 F diff --git a/input/molecule b/input/molecule index 81c624a..a1a3c86 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 2 2 2 0 0 # Znuc x y z - H 0. 0. 0. - H 0. 0. 1.4 + Li 0. 0. 0. + H 0. 0. 3.0141132 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..dfbf65a 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + Li 0.0000000000 0.0000000000 0.0000000000 + H 0.0000000000 0.0000000000 1.5950001314 diff --git a/input/options b/input/options index c47156a..5a8d22a 100644 --- a/input/options +++ b/input/options @@ -5,11 +5,11 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 F 1 # CIS/TDHF/BSE: singlet triplet - T F + T T # GF: maxSCF thresh DIIS n_diis renormalization 64 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 F F F F F 0.000 + 256 0.00001 T 5 F F T F F F F 0.000 # ACFDT: AC Kx XBS T F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/input/weight b/input/weight index 69b24c8..20bc6a7 100644 --- a/input/weight +++ b/input/weight @@ -1,30 +1,49 @@ -1 6 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 +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 S 1 - 1 0.3258000 1.0000000 + 1 0.0280500 1.0000000 S 1 - 1 0.1027000 1.0000000 + 1 0.0086400 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 P 1 - 1 1.4070000 1.0000000 + 1 0.0240300 1.0000000 P 1 - 1 0.3880000 1.0000000 + 1 0.0057900 1.0000000 D 1 - 1 1.0570000 1.0000000 -2 6 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 -S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 + 1 0.1239000 1.0000000 D 1 - 1 1.0570000 1.0000000 + 1 0.0725000 1.0000000 +2 5 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 +S 1 + 1 0.1220000 1.0000000 +S 1 + 1 0.0297400 1.0000000 +P 1 + 1 0.7270000 1.0000000 +P 1 + 1 0.1410000 1.0000000 + diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix.f90 index 3715424..2fb580c 100644 --- a/src/QuAcK/Bethe_Salpeter_A_matrix.f90 +++ b/src/QuAcK/Bethe_Salpeter_A_matrix.f90 @@ -46,4 +46,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) + 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 07d4e70..cf7573b 100644 --- a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 @@ -46,4 +46,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) end subroutine Bethe_Salpeter_B_matrix diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index d2d55c1..f11de8d 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -101,11 +101,11 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif ! Find all the roots of the QP equation if necessary - call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) +! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) ! Dump results -! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) ! Compute the RPA correlation energy diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index ca8a956..93be12e 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -70,19 +70,19 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1 do ab=1,nVV do ij=1,nOO - write(42,*) ab,ij,B(ab,ij) + 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 - write(43,*) ab,cd,C(ab,cd) + 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 - write(44,*) ij,kl,D(ij,kl) + if(abs(D(ij,kl)) > 1d-15) write(44,*) ij,kl,D(ij,kl) end do end do diff --git a/src/QuAcK/obj/.gitingore b/src/QuAcK/obj/.gitingore deleted file mode 100644 index 5761abc..0000000 --- a/src/QuAcK/obj/.gitingore +++ /dev/null @@ -1 +0,0 @@ -*.o diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index c26ea66..7063d7b 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/src/QuAcK/sort_ppRPA.f90 @@ -15,6 +15,13 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! Local variables integer :: pq,ab,ij + double precision,allocatable :: M(:,:) + double precision,allocatable :: Z1(:,:) + double precision,allocatable :: Z2(:,:) + double precision,allocatable :: S1(:,:) + double precision,allocatable :: S2(:,:) + double precision,allocatable :: O1(:,:) + double precision,allocatable :: O2(:,:) ! Output variables @@ -25,6 +32,13 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) double precision,intent(out) :: X2(nVV,nOO) double precision,intent(out) :: Y2(nOO,nOO) +! Memory allocation + + allocate(M(nOO+nVV,nOO+nVV), & + Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO), & + S1(nVV,nVV),S2(nOO,nOO), & + O1(nVV,nVV),O2(nOO,nOO)) + ! Initializatiom Omega1(:) = 0d0 @@ -35,6 +49,20 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) X2(:,:) = 0d0 Y2(:,:) = 0d0 +! Compute metric + + M(:,:) = 0d0 + + do ab=1,nVV + M(ab,ab) = 1d0 + end do + + do ij=1,nOO + M(nVV+ij,nVV+ij) = -1d0 + end do + +! Start sorting eigenvectors + ab = 0 ij = 0 @@ -44,15 +72,17 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ab = ab + 1 Omega1(ab) = Omega(pq) - X1(1:nVV,ab) = Z( 1: nVV,pq) - Y1(1:nOO,ab) = Z(nVV+1:nOO+nVV,pq) + Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq) +! X1(1:nVV,ab) = Z( 1: nVV,pq) +! Y1(1:nOO,ab) = Z(nVV+1:nOO+nVV,pq) else ij = ij + 1 Omega2(ij) = Omega(pq) - X2(1:nVV,ij) = Z( 1: nVV,pq) - Y2(1:nOO,ij) = Z(nVV+1:nOO+nVV,pq) + Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq) +! X2(1:nVV,ij) = Z( 1: nVV,pq) +! Y2(1:nOO,ij) = Z(nVV+1:nOO+nVV,pq) end if @@ -69,12 +99,29 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) call matout(nOO,1,Omega2(:)) write(*,*) - 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)) +! Orthogonalize eigenvectors + + S1 = matmul(transpose(Z1),matmul(M,Z1)) + S2 = - matmul(transpose(Z2),matmul(M,Z2)) + + call orthogonalization_matrix(1,nVV,S1,O1) + call orthogonalization_matrix(1,nOO,S2,O2) + + Z1 = matmul(Z1,O1) + Z2 = matmul(Z2,O2) + + X1(1:nVV,1:nVV) = Z1( 1: nVV,1:nVV) + Y1(1:nOO,1:nVV) = Z1(nVV+1:nOO+nVV,1:nVV) + + 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)) end subroutine sort_ppRPA