From 13b170fc2f8ab7dd394e2a83659a8ac4764eb86e Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 20 Apr 2020 12:28:19 +0200 Subject: [PATCH] star implementation of CISD --- input/basis | 44 +----- input/methods | 10 +- input/molecule | 4 +- input/molecule.xyz | 2 +- input/weight | 44 +----- src/QuAcK/CISD.f90 | 204 +++++++++++++++++++++++++ src/QuAcK/QuAcK.f90 | 61 ++++++-- src/QuAcK/linear_response_A_matrix.f90 | 2 +- src/QuAcK/ppRPA.f90 | 3 - src/QuAcK/read_methods.f90 | 25 +-- 10 files changed, 292 insertions(+), 107 deletions(-) create mode 100644 src/QuAcK/CISD.f90 diff --git a/input/basis b/input/basis index 6f3d2a9..6796e3b 100644 --- a/input/basis +++ b/input/basis @@ -1,39 +1,9 @@ -1 10 -S 8 - 1 24350.0000000 0.0005020 - 2 3650.0000000 0.0038810 - 3 829.6000000 0.0199970 - 4 234.0000000 0.0784180 - 5 75.6100000 0.2296760 - 6 26.7300000 0.4327220 - 7 9.9270000 0.3506420 - 8 1.1020000 -0.0076450 -S 8 - 1 24350.0000000 -0.0001180 - 2 3650.0000000 -0.0009150 - 3 829.6000000 -0.0047370 - 4 234.0000000 -0.0192330 - 5 75.6100000 -0.0603690 - 6 26.7300000 -0.1425080 - 7 9.9270000 -0.1777100 - 8 1.1020000 0.6058360 +1 3 +S 3 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 2.8360000 1.0000000 -S 1 - 1 0.3782000 1.0000000 -P 3 - 1 54.7000000 0.0171510 - 2 12.4300000 0.1076560 - 3 3.6790000 0.3216810 + 1 0.2976000 1.0000000 P 1 - 1 1.1430000 1.0000000 -P 1 - 1 0.3300000 1.0000000 -D 1 - 1 4.0140000 1.0000000 -D 1 - 1 1.0960000 1.0000000 -F 1 - 1 2.5440000 1.0000000 - - + 1 1.2750000 1.0000000 diff --git a/input/methods b/input/methods index bd52455..ff1d101 100644 --- a/input/methods +++ b/input/methods @@ -6,13 +6,15 @@ F F F # drCCD rCCD lCCD pCCD F F F F -# CIS RPA RPAx ppRPA ADC - F F F F F +# CIS CID CISD + F F T +# RPA RPAx ppRPA + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - T F F + F F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index edeba31..c78e87e 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 5 5 0 0 + 1 1 1 0 0 # Znuc x y z - Ne 0.0 0.0 0.0 + He 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 1c70680..797b5fc 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,3 @@ 1 - Ne 0.0000000000 0.0000000000 0.0000000000 + He 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index 6f3d2a9..6796e3b 100644 --- a/input/weight +++ b/input/weight @@ -1,39 +1,9 @@ -1 10 -S 8 - 1 24350.0000000 0.0005020 - 2 3650.0000000 0.0038810 - 3 829.6000000 0.0199970 - 4 234.0000000 0.0784180 - 5 75.6100000 0.2296760 - 6 26.7300000 0.4327220 - 7 9.9270000 0.3506420 - 8 1.1020000 -0.0076450 -S 8 - 1 24350.0000000 -0.0001180 - 2 3650.0000000 -0.0009150 - 3 829.6000000 -0.0047370 - 4 234.0000000 -0.0192330 - 5 75.6100000 -0.0603690 - 6 26.7300000 -0.1425080 - 7 9.9270000 -0.1777100 - 8 1.1020000 0.6058360 +1 3 +S 3 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 2.8360000 1.0000000 -S 1 - 1 0.3782000 1.0000000 -P 3 - 1 54.7000000 0.0171510 - 2 12.4300000 0.1076560 - 3 3.6790000 0.3216810 + 1 0.2976000 1.0000000 P 1 - 1 1.1430000 1.0000000 -P 1 - 1 0.3300000 1.0000000 -D 1 - 1 4.0140000 1.0000000 -D 1 - 1 1.0960000 1.0000000 -F 1 - 1 2.5440000 1.0000000 - - + 1 1.2750000 1.0000000 diff --git a/src/QuAcK/CISD.f90 b/src/QuAcK/CISD.f90 new file mode 100644 index 0000000..611d940 --- /dev/null +++ b/src/QuAcK/CISD.f90 @@ -0,0 +1,204 @@ +subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF) + +! Perform configuration interaction with singles and doubles + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + logical :: dump_trans = .false. + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb,iajb,kcld + integer :: ishift,jshift + integer :: ispin + integer :: nS + integer :: nD + integer :: nSD + double precision,external :: Kronecker_delta + double precision,allocatable :: H(:,:),Omega(:) + +! Hello world + + write(*,*) + write(*,*)'******************************************************' + write(*,*)'| Configuration Interaction with Singles and Doubles |' + write(*,*)'******************************************************' + write(*,*) + +! Compute CIS matrix + + if(singlet_manifold) then + + ispin = 1 + + ! Dimensions + + nS = (nO - nC)*(nV - nR) + nD = (nO - nC)*(nO - nC + 1)/2*(nV - nR)*(nV - nR + 1)/2 + nSD = 1 + nS + nD + + print*,'nS = ',nS + print*,'nD = ',nD + print*,'nSD = ',nSD + + ! Memory allocation + + allocate(H(nSD,nSD),Omega(nSD)) + + ! 0D block + + ishift = 0 + jshift = 1 + nS + + iajb = 0 + + do i=nC+1,nO + do a=1,nV-nR + do j=i,nO + do b=a,nV-nR + + iajb = iajb + 1 + H(ishift+1,jshift+iajb) = ERI(i,j,nO+a,nO+b) + + end do + end do + end do + end do + + ! SS block + + ishift = 1 + jshift = 1 + + ia = 0 + jb = 0 + + do i=nC+1,nO + do a=1,nV-nR + + ia = ia + 1 + + do j=nC+1,nO + do b=1,nV-nR + + jb = jb + 1 + + H(ishift+ia,jshift+jb) & + = Kronecker_delta(i,j)*Kronecker_delta(a,b)*(eHF(nO+a) - eHF(i)) & + + ERI(nO+a,j,i,nO+b) - ERI(nO+a,j,nO+b,i) + + end do + end do + + end do + end do + + ! SD block + + ishift = 1 + jshift = 1 + nS + + ia = 0 + kcld = 0 + + do i=nC+1,nO + do a=1,nV-nR + + ia = ia + 1 + + do k=nC+1,nO + do c=1,nV-nR + do l=k,nO + do d=c,nV-nR + + kcld = kcld + 1 + + H(ishift+ia,jshift+kcld) & + = Kronecker_delta(i,k)*(ERI(nO+a,l,nO+c,nO+d) - ERI(nO+a,l,nO+d,nO+c)) & + - Kronecker_delta(i,l)*(ERI(nO+a,k,nO+c,nO+d) - ERI(nO+a,k,nO+d,nO+c)) & + - Kronecker_delta(a,c)*(ERI(k,l,i,nO+d) - ERI(k,l,nO+d,i)) & + + Kronecker_delta(a,d)*(ERI(k,l,i,nO+c) - ERI(k,l,nO+c,i)) + + end do + end do + end do + end do + + end do + end do + + ! DD block + + ishift = 1 + nS + jshift = 1 + nS + + iajb = 0 + kcld = 0 + + do i=nC+1,nO + do a=1,nV-nR + do j=i,nO + do b=a,nV-nR + + iajb = iajb + 1 + + do k=nC+1,nO + do c=1,nV-nR + do l=k,nO + do d=c,nV-nR + + kcld = kcld + 1 + +! H(ishift+iajb,jshift+kcld) & +! = Kronecker_delta(i,k)*(ERI(a,l,c,d) - ERI(a,l,d,c)) & +! - Kronecker_delta(i,l)*(ERI(a,k,c,d) - ERI(a,k,d,c)) & +! - Kronecker_delta(a,c)*(ERI(k,l,i,d) - ERI(k,l,d,i)) & +! + Kronecker_delta(a,d)*(ERI(k,l,i,c) - ERI(k,l,c,i)) + + end do + end do + end do + end do + + end do + end do + end do + end do + + call diagonalize_matrix(nSD,H,Omega) + call print_excitation('CISD ',ispin,nS,Omega) + + if(dump_trans) then + print*,'Singlet CISD transition vectors' + call matout(nSD,nSD,H) + write(*,*) + endif + + endif + +! if(triplet_manifold) then + +! ispin = 2 +! +! call diagonalize_matrix(nSD,H,Omega) +! call print_excitation('CISD ',ispin,nSD,Omega) + +! if(dump_trans) then +! print*,'Triplet CIS transition vectors' +! call matout(nSD,nSD,H) +! write(*,*) +! endif + +! endif + +end subroutine CISD diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 7c8c831..ee79980 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -8,8 +8,9 @@ program QuAcK logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD - logical :: doCIS,doRPA,doRPAx - logical :: doppRPA,doADC + logical :: doCIS,doCID,doCISD + logical :: doRPA,doRPAx,doppRPA + logical :: doADC logical :: doG0F2,doevGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW logical :: doG0T0,doevGT,doqsGT @@ -61,6 +62,8 @@ program QuAcK double precision :: start_CCD ,end_CCD ,t_CCD double precision :: start_CCSD ,end_CCSD ,t_CCSD double precision :: start_CIS ,end_CIS ,t_CIS + double precision :: start_CID ,end_CID ,t_CID + double precision :: start_CISD ,end_CISD ,t_CISD double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_RPAx ,end_RPAx ,t_RPAx double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA @@ -128,8 +131,8 @@ program QuAcK doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doRPA,doRPAx, & - doppRPA,doADC, & + doCIS,doCID,doCISD, & + doRPA,doRPAx,doppRPA, & doG0F2,doevGF2,doG0F3,doevGF3, & doG0W0,doevGW,doqsGW, & doG0T0,doevGT,doqsGT, & @@ -491,6 +494,38 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Compute CID excitations +!------------------------------------------------------------------------ + + if(doCID) then + + call cpu_time(start_CID) +! call CID(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF) + call cpu_time(end_CID) + + t_CID = end_CID - start_CID + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CID,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute CISD excitations +!------------------------------------------------------------------------ + + if(doCISD) then + + call cpu_time(start_CISD) + call CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF) + call cpu_time(end_CISD) + + t_CISD = end_CISD - start_CISD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CISD,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute (direct) RPA excitations !------------------------------------------------------------------------ @@ -546,18 +581,18 @@ program QuAcK ! Compute ADC excitations !------------------------------------------------------------------------ - if(doADC) then +! if(doADC) then - call cpu_time(start_ADC) - call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & - nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO) - call cpu_time(end_ADC) +! call cpu_time(start_ADC) +! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & +! nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO) +! call cpu_time(end_ADC) - t_ADC = end_ADC - start_ADC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ADC = ',t_ADC,' seconds' - write(*,*) +! t_ADC = end_ADC - start_ADC +! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ADC = ',t_ADC,' seconds' +! write(*,*) - end if +! end if !------------------------------------------------------------------------ ! Compute G0F2 electronic binding energies diff --git a/src/QuAcK/linear_response_A_matrix.f90 b/src/QuAcK/linear_response_A_matrix.f90 index 1a1d391..f2c20e0 100644 --- a/src/QuAcK/linear_response_A_matrix.f90 +++ b/src/QuAcK/linear_response_A_matrix.f90 @@ -16,7 +16,7 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, ! Local variables double precision :: delta_spin,delta_dRPA - double precision :: Kronecker_delta + double precision,external :: Kronecker_delta integer :: i,j,a,b,ia,jb diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index 7c10627..3c45cc5 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -56,9 +56,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER nOO = nO*(nO+1)/2 nVV = nV*(nV+1)/2 -! print*,'nOO = ',nOO -! print*,'nVV = ',nVV - ! Memory allocation allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index ae8ede1..ec9afc1 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -2,8 +2,8 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doRPA,doRPAx, & - doppRPA,doADC, & + doCIS,doCID,doCISD, & + doRPA,doRPAx,doppRPA, & doG0F2,doevGF2,doG0F3,doevGF3, & doG0W0,doevGW,doqsGW, & doG0T0,doevGT,doqsGT, & @@ -19,7 +19,8 @@ subroutine read_methods(doRHF,doUHF,doMOM, & logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doCCD,doCCSD,doCCSDT logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD - logical,intent(out) :: doCIS,doRPA,doRPAx,doppRPA,doADC + logical,intent(out) :: doCIS,doCID,doCISD + logical,intent(out) :: doRPA,doRPAx,doppRPA logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0T0,doevGT,doqsGT @@ -53,10 +54,12 @@ subroutine read_methods(doRHF,doUHF,doMOM, & do_pCCD = .false. doCIS = .false. + doCID = .false. + doCISD = .false. + doRPA = .false. doRPAx = .false. doppRPA = .false. - doADC = .false. doG0F2 = .false. doevGF2 = .false. @@ -108,12 +111,16 @@ subroutine read_methods(doRHF,doUHF,doMOM, & ! Read excited state methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 + read(1,*) answer1,answer2,answer3 if(answer1 == 'T') doCIS = .true. - if(answer2 == 'T') doRPA = .true. - if(answer3 == 'T') doRPAx = .true. - if(answer4 == 'T') doppRPA = .true. - if(answer5 == 'T') doADC = .true. + if(answer2 == 'T') doCID = .true. + if(answer3 == 'T') doCISD = .true. + + read(1,*) + read(1,*) answer1,answer2,answer3 + if(answer1 == 'T') doRPA = .true. + if(answer2 == 'T') doRPAx = .true. + if(answer3 == 'T') doppRPA = .true. ! Read Green function methods