diff --git a/input/basis b/input/basis index f19a2d0..1ea2746 100644 --- a/input/basis +++ b/input/basis @@ -1,30 +1,30 @@ 1 6 S 8 - 1 17880.0000000 0.0007380 - 2 2683.0000000 0.0056770 - 3 611.5000000 0.0288830 - 4 173.5000000 0.1085400 - 5 56.6400000 0.2909070 - 6 20.4200000 0.4483240 - 7 7.8100000 0.2580260 - 8 1.6530000 0.0150630 + 1 2940.0000000 0.0006800 + 2 441.2000000 0.0052360 + 3 100.5000000 0.0266060 + 4 28.4300000 0.0999930 + 5 9.1690000 0.2697020 + 6 3.1960000 0.4514690 + 7 1.1590000 0.2950740 + 8 0.1811000 0.0125870 S 8 - 1 17880.0000000 -0.0001720 - 2 2683.0000000 -0.0013570 - 3 611.5000000 -0.0067370 - 4 173.5000000 -0.0276630 - 5 56.6400000 -0.0762080 - 6 20.4200000 -0.1752270 - 7 7.8100000 -0.1070380 - 8 1.6530000 0.5670500 + 1 2940.0000000 -0.0001230 + 2 441.2000000 -0.0009660 + 3 100.5000000 -0.0048310 + 4 28.4300000 -0.0193140 + 5 9.1690000 -0.0532800 + 6 3.1960000 -0.1207230 + 7 1.1590000 -0.1334350 + 8 0.1811000 0.5307670 S 1 - 1 0.4869000 1.0000000 + 1 0.0589000 1.0000000 P 3 - 1 28.3900000 0.0460870 - 2 6.2700000 0.2401810 - 3 1.6950000 0.5087440 + 1 3.6190000 0.0291110 + 2 0.7110000 0.1693650 + 3 0.1951000 0.5134580 P 1 - 1 0.4317000 1.0000000 + 1 0.0601800 1.0000000 D 1 - 1 2.2020000 1.0000000 + 1 0.2380000 1.0000000 diff --git a/input/methods b/input/methods index 7895c92..7e8729e 100644 --- a/input/methods +++ b/input/methods @@ -3,7 +3,7 @@ # MP2 MP3 MP2-F12 T F F # CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD - F F F F F F T + F F F T F T T # CIS RPA RPAx ppRPA ADC F F F F F # G0F2 evGF2 G0F3 evGF3 diff --git a/input/molecule b/input/molecule index edeba31..6a6f6d1 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 5 5 0 0 + 1 2 2 0 0 # Znuc x y z - Ne 0.0 0.0 0.0 + Be 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 1c70680..8023e37 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,3 @@ 1 - Ne 0.0000000000 0.0000000000 0.0000000000 + Be 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index 2206914..7881bfa 100644 --- a/input/options +++ b/input/options @@ -3,7 +3,7 @@ # MP: # CC: maxSCF thresh DIIS n_diis - 64 0.0000001 F 1 + 64 0.0000001 T 5 # CIS/TDHF/BSE: singlet triplet T T # GF: maxSCF thresh DIIS n_diis lin renorm diff --git a/input/weight b/input/weight index f19a2d0..1ea2746 100644 --- a/input/weight +++ b/input/weight @@ -1,30 +1,30 @@ 1 6 S 8 - 1 17880.0000000 0.0007380 - 2 2683.0000000 0.0056770 - 3 611.5000000 0.0288830 - 4 173.5000000 0.1085400 - 5 56.6400000 0.2909070 - 6 20.4200000 0.4483240 - 7 7.8100000 0.2580260 - 8 1.6530000 0.0150630 + 1 2940.0000000 0.0006800 + 2 441.2000000 0.0052360 + 3 100.5000000 0.0266060 + 4 28.4300000 0.0999930 + 5 9.1690000 0.2697020 + 6 3.1960000 0.4514690 + 7 1.1590000 0.2950740 + 8 0.1811000 0.0125870 S 8 - 1 17880.0000000 -0.0001720 - 2 2683.0000000 -0.0013570 - 3 611.5000000 -0.0067370 - 4 173.5000000 -0.0276630 - 5 56.6400000 -0.0762080 - 6 20.4200000 -0.1752270 - 7 7.8100000 -0.1070380 - 8 1.6530000 0.5670500 + 1 2940.0000000 -0.0001230 + 2 441.2000000 -0.0009660 + 3 100.5000000 -0.0048310 + 4 28.4300000 -0.0193140 + 5 9.1690000 -0.0532800 + 6 3.1960000 -0.1207230 + 7 1.1590000 -0.1334350 + 8 0.1811000 0.5307670 S 1 - 1 0.4869000 1.0000000 + 1 0.0589000 1.0000000 P 3 - 1 28.3900000 0.0460870 - 2 6.2700000 0.2401810 - 3 1.6950000 0.5087440 + 1 3.6190000 0.0291110 + 2 0.7110000 0.1693650 + 3 0.1951000 0.5134580 P 1 - 1 0.4317000 1.0000000 + 1 0.0601800 1.0000000 D 1 - 1 2.2020000 1.0000000 + 1 0.2380000 1.0000000 diff --git a/src/QuAcK/CCD.f90 b/src/QuAcK/CCD.f90 index 4f3c3d0..20ba3a0 100644 --- a/src/QuAcK/CCD.f90 +++ b/src/QuAcK/CCD.f90 @@ -48,6 +48,11 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + ! Hello world write(*,*) @@ -110,6 +115,10 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) EcMP4 = 0d0 +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + ! Initialization allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) @@ -118,6 +127,10 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) Conv = 1d0 nSCF = 0 + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -169,6 +182,15 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) ECCD = ERHF + EcCCD + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' diff --git a/src/QuAcK/CCSD.f90 b/src/QuAcK/CCSD.f90 index e005062..3240b44 100644 --- a/src/QuAcK/CCSD.f90 +++ b/src/QuAcK/CCSD.f90 @@ -63,6 +63,14 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) double precision,allocatable :: t2(:,:,:,:) double precision,allocatable :: tau(:,:,:,:) + integer :: n_diis + double precision :: rcond1 + double precision :: rcond2 + double precision,allocatable :: err1_diis(:,:) + double precision,allocatable :: err2_diis(:,:) + double precision,allocatable :: t1_diis(:,:) + double precision,allocatable :: t2_diis(:,:) + ! Hello world write(*,*) @@ -145,9 +153,20 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) aoooo(nO,nO,nO,nO),bvvvv(nV,nV,nV,nV),hovvo(nO,nV,nV,nO), & r1(nO,nV),r2(nO,nO,nV,nV)) +! Memory allocation for DIIS + + allocate(err1_diis(nO*nV ,max_diis),t1_diis(nO*nV ,max_diis), & + err2_diis(nO*nO*nV*nV,max_diis),t2_diis(nO*nO*nV*nV,max_diis)) + Conv = 1d0 nSCF = 0 + n_diis = 0 + t1_diis(:,:) = 0d0 + t2_diis(:,:) = 0d0 + err1_diis(:,:) = 0d0 + err2_diis(:,:) = 0d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -200,6 +219,16 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) ECCSD = ERHF + EcCCSD + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond1,nO*nV ,nO*nV ,n_diis,err1_diis,t1_diis,-r1/delta_OV ,t1) + call DIIS_extrapolation(rcond2,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,err2_diis,t2_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(min(abs(rcond1),abs(rcond2)) < 1d-15) n_diis = 0 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCSD+ENuc,'|',EcCCSD,'|',Conv,'|' diff --git a/src/QuAcK/drCCD.f90 b/src/QuAcK/drCCD.f90 index df63a2d..b3cbca0 100644 --- a/src/QuAcK/drCCD.f90 +++ b/src/QuAcK/drCCD.f90 @@ -47,6 +47,12 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) + + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + ! Hello world write(*,*) @@ -103,6 +109,10 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) EcMP4 = 0d0 +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + ! Initialization allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) @@ -111,6 +121,10 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) Conv = 1d0 nSCF = 0 + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -152,6 +166,15 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) ECCD = ERHF + EcCCD + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' diff --git a/src/QuAcK/lCCD.f90 b/src/QuAcK/lCCD.f90 index 7ca2c20..c92c5e7 100644 --- a/src/QuAcK/lCCD.f90 +++ b/src/QuAcK/lCCD.f90 @@ -48,6 +48,11 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + ! Hello world write(*,*) @@ -110,6 +115,10 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) EcMP4 = 0d0 +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + ! Initialization allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) @@ -118,6 +127,10 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) Conv = 1d0 nSCF = 0 + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -159,6 +172,15 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) ECCD = ERHF + EcCCD + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' diff --git a/src/QuAcK/pCCD.f90 b/src/QuAcK/pCCD.f90 index 7c8d5bf..49e44f7 100644 --- a/src/QuAcK/pCCD.f90 +++ b/src/QuAcK/pCCD.f90 @@ -36,8 +36,11 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) double precision,allocatable :: r(:,:) double precision,allocatable :: t(:,:) - double precision,external :: trace_matrix - + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + ! Hello world write(*,*) @@ -86,6 +89,10 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) t(:,:) = - OOVV(:,:)/delta_OOVV(:,:) +! Memory allocation for DIIS + + allocate(error_diis(nO*nV,max_diis),t_diis(nO*nV,max_diis)) + ! Initialization allocate(r(nO,nV),y(nO,nO)) @@ -93,6 +100,10 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) Conv = 1d0 nSCF = 0 + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -121,7 +132,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) end do end do - ! Compute residual + ! Compute residual do i=nC+1,nO do a=1,nV-nR @@ -140,7 +151,6 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) end do end do - ! Check convergence Conv = maxval(abs(r(:,:))) @@ -158,10 +168,19 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) end do end do -! Dump results + ! Dump results ECCD = ERHF + EcCCD + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,error_diis,t_diis,-r/delta_OOVV,t) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' diff --git a/src/QuAcK/rCCD.f90 b/src/QuAcK/rCCD.f90 index e53dd52..5331be2 100644 --- a/src/QuAcK/rCCD.f90 +++ b/src/QuAcK/rCCD.f90 @@ -48,6 +48,11 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) + integer :: n_diis + double precision :: rcond + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: t_diis(:,:) + ! Hello world write(*,*) @@ -110,6 +115,11 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) EcMP4 = 0d0 + +! Memory allocation for DIIS + + allocate(error_diis(nO*nO*nV*nV,max_diis),t_diis(nO*nO*nV*nV,max_diis)) + ! Initialization allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) @@ -118,6 +128,10 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) Conv = 1d0 nSCF = 0 + n_diis = 0 + t_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ @@ -159,6 +173,15 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) ECCD = ERHF + EcCCD + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nO*nV*nV,nO*nO*nV*nV,n_diis,error_diis,t_diis,-r2/delta_OOVV,t2) + + ! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'