diff --git a/examples/molecule.HCl b/examples/molecule.HCl index 7336f5b..88a5912 100644 --- a/examples/molecule.HCl +++ b/examples/molecule.HCl @@ -2,4 +2,4 @@ 2 9 9 0 0 # Znuc x y z H 0. 0. 0. - Cl 0. 0. 2.42 + Cl 0. 0. 2.41 diff --git a/examples/molecule.LiH b/examples/molecule.LiH index aad52fc..d4f626f 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.01 + H 0. 0. 3 diff --git a/examples/molecule.N2 b/examples/molecule.N2 index 93b448f..6e92c28 100644 --- a/examples/molecule.N2 +++ b/examples/molecule.N2 @@ -2,4 +2,4 @@ 2 7 7 0 0 # Znuc x y z N 0. 0. 0. - N 0. 0. 3.4 + N 0. 0. 2.07 diff --git a/input/basis b/input/basis index 80356bb..bb77ea0 100644 --- a/input/basis +++ b/input/basis @@ -1,64 +1,88 @@ 1 6 -S 9 -1 9.046000E+03 7.000000E-04 -2 1.357000E+03 5.389000E-03 -3 3.093000E+02 2.740600E-02 -4 8.773000E+01 1.032070E-01 -5 2.856000E+01 2.787230E-01 -6 1.021000E+01 4.485400E-01 -7 3.838000E+00 2.782380E-01 -8 7.466000E-01 1.544000E-02 -9 2.248000E-01 -2.864000E-03 -S 9 -1 9.046000E+03 -1.530000E-04 -2 1.357000E+03 -1.208000E-03 -3 3.093000E+02 -5.992000E-03 -4 8.773000E+01 -2.454400E-02 -5 2.856000E+01 -6.745900E-02 -6 1.021000E+01 -1.580780E-01 -7 3.838000E+00 -1.218310E-01 -8 7.466000E-01 5.490030E-01 -9 2.248000E-01 5.788150E-01 +S 3 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 -1 2.248000E-01 1.000000E+00 -P 4 -1 1.355000E+01 3.991900E-02 -2 2.917000E+00 2.171690E-01 -3 7.973000E-01 5.103190E-01 -4 2.185000E-01 4.622140E-01 -P 1 -1 2.185000E-01 1.000000E+00 -D 1 -1 8.170000E-01 1.0000000 -2 6 -S 9 -1 9.046000E+03 7.000000E-04 -2 1.357000E+03 5.389000E-03 -3 3.093000E+02 2.740600E-02 -4 8.773000E+01 1.032070E-01 -5 2.856000E+01 2.787230E-01 -6 1.021000E+01 4.485400E-01 -7 3.838000E+00 2.782380E-01 -8 7.466000E-01 1.544000E-02 -9 2.248000E-01 -2.864000E-03 -S 9 -1 9.046000E+03 -1.530000E-04 -2 1.357000E+03 -1.208000E-03 -3 3.093000E+02 -5.992000E-03 -4 8.773000E+01 -2.454400E-02 -5 2.856000E+01 -6.745900E-02 -6 1.021000E+01 -1.580780E-01 -7 3.838000E+00 -1.218310E-01 -8 7.466000E-01 5.490030E-01 -9 2.248000E-01 5.788150E-01 + 1 0.3258000 1.0000000 S 1 -1 2.248000E-01 1.000000E+00 -P 4 -1 1.355000E+01 3.991900E-02 -2 2.917000E+00 2.171690E-01 -3 7.973000E-01 5.103190E-01 -4 2.185000E-01 4.622140E-01 + 1 0.1027000 1.0000000 P 1 -1 2.185000E-01 1.000000E+00 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 D 1 -1 8.170000E-01 1.0000000 + 1 1.0570000 1.0000000 +2 12 +S 13 + 1 456100.0000000 0.492970E-04 + 2 68330.0000000 0.383029E-03 + 3 15550.0000000 0.200854E-02 + 4 4405.0000000 0.838558E-02 + 5 1439.0000000 0.294703E-01 + 6 520.4000000 0.878325E-01 + 7 203.1000000 0.211473E+00 + 8 83.9600000 0.365364E+00 + 9 36.2000000 0.340884E+00 + 10 15.8300000 0.102133E+00 + 11 6.3340000 0.311675E-02 + 12 2.6940000 0.105751E-02 + 13 0.4313000 0.156136E-03 +S 13 + 1 456100.0000000 -0.138304E-04 + 2 68330.0000000 -0.107279E-03 + 3 15550.0000000 -0.565083E-03 + 4 4405.0000000 -0.236135E-02 + 5 1439.0000000 -0.845886E-02 + 6 520.4000000 -0.259638E-01 + 7 203.1000000 -0.686362E-01 + 8 83.9600000 -0.141874E+00 + 9 36.2000000 -0.199319E+00 + 10 15.8300000 -0.195662E-01 + 11 6.3340000 0.499741E+00 + 12 2.6940000 0.563736E+00 + 13 0.4313000 -0.835091E-02 +S 13 + 1 456100.0000000 0.418546E-05 + 2 68330.0000000 0.324395E-04 + 3 15550.0000000 0.171105E-03 + 4 4405.0000000 0.714176E-03 + 5 1439.0000000 0.256705E-02 + 6 520.4000000 0.788552E-02 + 7 203.1000000 0.210867E-01 + 8 83.9600000 0.442264E-01 + 9 36.2000000 0.651670E-01 + 10 15.8300000 0.603012E-02 + 11 6.3340000 -0.206495E+00 + 12 2.6940000 -0.405871E+00 + 13 0.4313000 0.725661E+00 +S 1 + 1 0.9768000 1.0000000 +S 1 + 1 0.1625000 1.0000000 +P 7 + 1 663.3000000 0.240448E-02 + 2 156.8000000 0.192148E-01 + 3 49.9800000 0.885097E-01 + 4 18.4200000 0.256020E+00 + 5 7.2400000 0.436927E+00 + 6 2.9220000 0.350334E+00 + 7 0.3818000 -0.458423E-02 +P 7 + 1 663.3000000 -0.652145E-03 + 2 156.8000000 -0.519445E-02 + 3 49.9800000 -0.246938E-01 + 4 18.4200000 -0.728167E-01 + 5 7.2400000 -0.134030E+00 + 6 2.9220000 -0.947742E-01 + 7 0.3818000 0.564667E+00 +P 1 + 1 1.0220000 1.0000000 +P 1 + 1 0.1301000 1.0000000 +D 1 + 1 1.0460000 1.0000000 +D 1 + 1 0.3440000 1.0000000 +F 1 + 1 0.7060000 1.0000000 diff --git a/input/molecule b/input/molecule index 93b448f..88a5912 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 7 7 0 0 + 2 9 9 0 0 # Znuc x y z - N 0. 0. 0. - N 0. 0. 3.4 + H 0. 0. 0. + Cl 0. 0. 2.41 diff --git a/input/molecule.xyz b/input/molecule.xyz index b3ef6c5..64d6254 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - N 0.0000000000 0.0000000000 0.0000000000 - N 0.0000000000 0.0000000000 1.7992026466 + H 0.0000000000 0.0000000000 0.0000000000 + Cl 0.0000000000 0.0000000000 1.2753171701 diff --git a/input/weight b/input/weight index 80356bb..bb77ea0 100644 --- a/input/weight +++ b/input/weight @@ -1,64 +1,88 @@ 1 6 -S 9 -1 9.046000E+03 7.000000E-04 -2 1.357000E+03 5.389000E-03 -3 3.093000E+02 2.740600E-02 -4 8.773000E+01 1.032070E-01 -5 2.856000E+01 2.787230E-01 -6 1.021000E+01 4.485400E-01 -7 3.838000E+00 2.782380E-01 -8 7.466000E-01 1.544000E-02 -9 2.248000E-01 -2.864000E-03 -S 9 -1 9.046000E+03 -1.530000E-04 -2 1.357000E+03 -1.208000E-03 -3 3.093000E+02 -5.992000E-03 -4 8.773000E+01 -2.454400E-02 -5 2.856000E+01 -6.745900E-02 -6 1.021000E+01 -1.580780E-01 -7 3.838000E+00 -1.218310E-01 -8 7.466000E-01 5.490030E-01 -9 2.248000E-01 5.788150E-01 +S 3 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 -1 2.248000E-01 1.000000E+00 -P 4 -1 1.355000E+01 3.991900E-02 -2 2.917000E+00 2.171690E-01 -3 7.973000E-01 5.103190E-01 -4 2.185000E-01 4.622140E-01 -P 1 -1 2.185000E-01 1.000000E+00 -D 1 -1 8.170000E-01 1.0000000 -2 6 -S 9 -1 9.046000E+03 7.000000E-04 -2 1.357000E+03 5.389000E-03 -3 3.093000E+02 2.740600E-02 -4 8.773000E+01 1.032070E-01 -5 2.856000E+01 2.787230E-01 -6 1.021000E+01 4.485400E-01 -7 3.838000E+00 2.782380E-01 -8 7.466000E-01 1.544000E-02 -9 2.248000E-01 -2.864000E-03 -S 9 -1 9.046000E+03 -1.530000E-04 -2 1.357000E+03 -1.208000E-03 -3 3.093000E+02 -5.992000E-03 -4 8.773000E+01 -2.454400E-02 -5 2.856000E+01 -6.745900E-02 -6 1.021000E+01 -1.580780E-01 -7 3.838000E+00 -1.218310E-01 -8 7.466000E-01 5.490030E-01 -9 2.248000E-01 5.788150E-01 + 1 0.3258000 1.0000000 S 1 -1 2.248000E-01 1.000000E+00 -P 4 -1 1.355000E+01 3.991900E-02 -2 2.917000E+00 2.171690E-01 -3 7.973000E-01 5.103190E-01 -4 2.185000E-01 4.622140E-01 + 1 0.1027000 1.0000000 P 1 -1 2.185000E-01 1.000000E+00 + 1 1.4070000 1.0000000 +P 1 + 1 0.3880000 1.0000000 D 1 -1 8.170000E-01 1.0000000 + 1 1.0570000 1.0000000 +2 12 +S 13 + 1 456100.0000000 0.492970E-04 + 2 68330.0000000 0.383029E-03 + 3 15550.0000000 0.200854E-02 + 4 4405.0000000 0.838558E-02 + 5 1439.0000000 0.294703E-01 + 6 520.4000000 0.878325E-01 + 7 203.1000000 0.211473E+00 + 8 83.9600000 0.365364E+00 + 9 36.2000000 0.340884E+00 + 10 15.8300000 0.102133E+00 + 11 6.3340000 0.311675E-02 + 12 2.6940000 0.105751E-02 + 13 0.4313000 0.156136E-03 +S 13 + 1 456100.0000000 -0.138304E-04 + 2 68330.0000000 -0.107279E-03 + 3 15550.0000000 -0.565083E-03 + 4 4405.0000000 -0.236135E-02 + 5 1439.0000000 -0.845886E-02 + 6 520.4000000 -0.259638E-01 + 7 203.1000000 -0.686362E-01 + 8 83.9600000 -0.141874E+00 + 9 36.2000000 -0.199319E+00 + 10 15.8300000 -0.195662E-01 + 11 6.3340000 0.499741E+00 + 12 2.6940000 0.563736E+00 + 13 0.4313000 -0.835091E-02 +S 13 + 1 456100.0000000 0.418546E-05 + 2 68330.0000000 0.324395E-04 + 3 15550.0000000 0.171105E-03 + 4 4405.0000000 0.714176E-03 + 5 1439.0000000 0.256705E-02 + 6 520.4000000 0.788552E-02 + 7 203.1000000 0.210867E-01 + 8 83.9600000 0.442264E-01 + 9 36.2000000 0.651670E-01 + 10 15.8300000 0.603012E-02 + 11 6.3340000 -0.206495E+00 + 12 2.6940000 -0.405871E+00 + 13 0.4313000 0.725661E+00 +S 1 + 1 0.9768000 1.0000000 +S 1 + 1 0.1625000 1.0000000 +P 7 + 1 663.3000000 0.240448E-02 + 2 156.8000000 0.192148E-01 + 3 49.9800000 0.885097E-01 + 4 18.4200000 0.256020E+00 + 5 7.2400000 0.436927E+00 + 6 2.9220000 0.350334E+00 + 7 0.3818000 -0.458423E-02 +P 7 + 1 663.3000000 -0.652145E-03 + 2 156.8000000 -0.519445E-02 + 3 49.9800000 -0.246938E-01 + 4 18.4200000 -0.728167E-01 + 5 7.2400000 -0.134030E+00 + 6 2.9220000 -0.947742E-01 + 7 0.3818000 0.564667E+00 +P 1 + 1 1.0220000 1.0000000 +P 1 + 1 0.1301000 1.0000000 +D 1 + 1 1.0460000 1.0000000 +D 1 + 1 0.3440000 1.0000000 +F 1 + 1 0.7060000 1.0000000 diff --git a/scan_LiF.sh b/scan_LiF.sh index bea7182..13226bc 100755 --- a/scan_LiF.sh +++ b/scan_LiF.sh @@ -1,9 +1,9 @@ #! /bin/bash MOL="LiF" -BASIS="cc-pvdz" -R_START=3.001 -R_END=3.003 +BASIS="cc-pvtz" +R_START=3.980 +R_END=2.990 DR=0.001 for R in $(seq $R_START $DR $R_END) diff --git a/scan_LiH.sh b/scan_LiH.sh index 8026108..4eed6e1 100755 --- a/scan_LiH.sh +++ b/scan_LiH.sh @@ -2,8 +2,8 @@ MOL="LiH" BASIS="cc-pvqz" -R_START=3.010 -R_END=3.020 +R_START=2.990 +R_END=3.000 DR=0.001 for R in $(seq $R_START $DR $R_END) diff --git a/scan_N2.sh b/scan_N2.sh index 9b867a7..c62b21e 100755 --- a/scan_N2.sh +++ b/scan_N2.sh @@ -1,10 +1,10 @@ #! /bin/bash MOL="N2" -BASIS="cc-pvdz" -R_START=1.5 -R_END=3.5 -DR=0.1 +BASIS="cc-pvtz" +R_START=2.066 +R_END=2.080 +DR=0.001 for R in $(seq $R_START $DR $R_END) do diff --git a/src/QuAcK/ACFDT.f90 b/src/QuAcK/ACFDT.f90 index 6a31b16..60e021b 100644 --- a/src/QuAcK/ACFDT.f90 +++ b/src/QuAcK/ACFDT.f90 @@ -52,13 +52,14 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man end if + EcAC(:) = 0d0 + Ec(:,:) = 0d0 + ! Singlet manifold if(singlet_manifold) then ispin = 1 - EcAC(ispin) = 0d0 - Ec(:,ispin) = 0d0 write(*,*) '--------------' write(*,*) 'Singlet states' @@ -104,8 +105,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man if(triplet_manifold) then ispin = 2 - EcAC(ispin) = 0d0 - Ec(:,ispin) = 0d0 write(*,*) '--------------' write(*,*) 'Triplet states' @@ -120,20 +119,20 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man lambda = rAC(iAC) -! if(doXBS) then + if(doXBS) then -! call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & -! rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) -! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,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 + 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)) + 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) + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) end do diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 new file mode 100644 index 0000000..9d5eb3a --- /dev/null +++ b/src/QuAcK/RPAx.f90 @@ -0,0 +1,126 @@ +subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) + +! Perform random phase approximation calculation with exchange (aka TDHF) + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: singlet_manifold + double precision,intent(in) :: eta + 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 :: EcRPAx(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'***********************************************************' + write(*,*)'| Random phase approximation calculation with exchange |' + write(*,*)'***********************************************************' + write(*,*) + +! Initialization + + EcRPAx(:) = 0d0 + EcAC(:) = 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.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) + + endif + +! Triplet manifold + + if(triplet_manifold) then + + ispin = 2 + + call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) + + endif + + if(exchange_kernel) then + + EcRPAx(1) = 0.5d0*EcRPAx(1) + EcRPAx(2) = 1.5d0*EcRPAx(2) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy (singlet) =',EcRPAx(1) + write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy (triplet) =',EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of RPAx correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) + + call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, & + nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC) + + if(exchange_kernel) then + + EcAC(1) = 0.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) + + end if + + 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(*,*) + + end if + +end subroutine RPAx