diff --git a/mol/H2O.xyz b/mol/H2O.xyz index 00a490d..2e291f6 100644 --- a/mol/H2O.xyz +++ b/mol/H2O.xyz @@ -1,5 +1,5 @@ 3 O 0.0000 0.0000 0.0000 -H 0.9591 0.0000 0.0000 -H -0.2373 0.9293 0.0000 +H 0.7571 0.0000 0.5861 +H -0.7571 0.0000 0.5861 diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index 7e9b976..ddda574 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -59,11 +59,13 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, integer :: max_diis,n_diis double precision :: rcond double precision,allocatable :: err_diis(:,:) - double precision,allocatable :: Om_diis(:,:) + double precision,allocatable :: Phi_diis(:,:) double precision,allocatable :: err(:) - double precision,allocatable :: Om(:) + double precision,allocatable :: Phi(:) double precision :: alpha + integer ::p,q,r,s,pqrs + ! Output variables ! None @@ -75,15 +77,15 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! DIIS parameters - ! max_diis = 10 - ! n_diis = 0 - ! rcond = 0d0 + max_diis = 1 + n_diis = 0 + rcond = 1d0 - ! allocate(err_diis(nS+nOO+nVV,max_diis),Om_diis(nS+nOO+nVV,max_diis)) - ! allocate(err(nS+nOO+nVV),Om(nS+nOO+nVV)) + allocate(err_diis(2*nOrb**4,max_diis),Phi_diis(2*nOrb**4,max_diis)) + allocate(err(2*nOrb**4),Phi(2*nOrb**4)) - ! err_diis(:,:) = 0d0 - ! Om_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + Phi_diis(:,:) = 0d0 ! Start @@ -365,38 +367,64 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:) ! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:) - err_eh = maxval(abs(old_eh_Phi - eh_Phi)) - err_pp = maxval(abs(old_pp_Phi - pp_Phi)) + err_eh = maxval(abs(eh_Phi - old_eh_Phi)) + err_pp = maxval(abs(pp_Phi - old_pp_Phi)) + + call matout(nOrb**2,nOrb**2,eh_Phi - old_eh_Phi) + call matout(nOrb**2,nOrb**2,pp_Phi - old_pp_Phi) + + !--------------------! + ! DIIS extrapolation ! + !--------------------! + + pqrs = 0 + do p=1,nOrb + do q=1,nOrb + do r=1,nOrb + do s=1,nOrb + pqrs = pqrs + 1 + + err( pqrs) = eh_Phi(p,q,r,s) - old_eh_Phi(p,q,r,s) + err(nOrb**4+pqrs) = pp_Phi(p,q,r,s) - old_pp_Phi(p,q,r,s) + + Phi( pqrs) = eh_Phi(p,q,r,s) + Phi(nOrb**4+pqrs) = pp_Phi(p,q,r,s) + + end do + end do + end do + end do + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,2*nOrb**4,2*nOrb**4,n_diis,err_diis,Phi_diis,err,Phi) + print*,rcond + + end if + + pqrs = 0 + do p=1,nOrb + do q=1,nOrb + do r=1,nOrb + do s=1,nOrb + pqrs = pqrs + 1 + + eh_Phi(p,q,r,s) = Phi( pqrs) + pp_Phi(p,q,r,s) = Phi(nOrb**4+pqrs) + + end do + end do + end do + end do old_eh_Phi(:,:,:,:) = eh_Phi(:,:,:,:) old_pp_Phi(:,:,:,:) = pp_Phi(:,:,:,:) ! Free memory + deallocate(eh_Phi,pp_Phi) - !--------------------! - ! DIIS extrapolation ! - !--------------------! - - ! err( 1:nS ) = eh_Om(:) - old_eh_Om(:) - ! err( nS+1:nS+nVV ) = ee_Om(:) - old_ee_Om(:) - ! err(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) - old_hh_Om(:) - - ! Om( 1:nS ) = eh_Om(:) - ! Om( nS+1:nS+nVV ) = ee_Om(:) - ! Om(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) - - ! if(max_diis > 1) then - - ! n_diis = min(n_diis+1,max_diis) - ! call DIIS_extrapolation(rcond,nS+nOO+nVV,nS+nOO+nVV,n_diis,err_diis,Om_diis,err,Om) - - ! end if - - ! eh_Om(:) = Om( 1:nS ) - ! ee_Om(:) = Om( nS+1:nS+nVV ) - ! hh_Om(:) = Om(nVV+nS+1:nS+nVV+nOO) - write(*,*) '----------------------------------------' write(*,*) ' Two-body (kernel) convergence ' write(*,*) '----------------------------------------' @@ -404,7 +432,6 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, write(*,'(1X,A30,F10.6)')'Error for pp channel = ',err_pp write(*,*) '----------------------------------------' write(*,*) - ! Convergence criteria err_2b = max(err_eh,err_pp) diff --git a/src/Parquet/G_screened_integrals.f90 b/src/Parquet/G_screened_integrals.f90 index 594c887..9d40221 100644 --- a/src/Parquet/G_screened_integrals.f90 +++ b/src/Parquet/G_screened_integrals.f90 @@ -38,9 +38,9 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho rho(p,q,ia) = rho(p,q,ia) & + (ERI(q,j,p,b) - ERI(q,j,b,p) & - - eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X & + - 1d0*eh_Phi(q,j,b,p) + 1d0*pp_Phi(q,j,p,b)) * X & + (ERI(q,b,p,j) - ERI(q,b,j,p) & - - eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y + - 1d0*eh_Phi(q,b,j,p) + 1d0*pp_Phi(q,b,p,j)) * Y end do end do @@ -108,7 +108,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do d=c+1,nOrb-nR cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & - + eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X1(cd,ab) + + 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X1(cd,ab) end do end do @@ -117,7 +117,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do l=k+1,nO kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & - + eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y1(kl,ab) + + 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab) end do end do @@ -133,7 +133,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do d=c+1,nOrb-nR cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & - + eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X2(cd,ij) + + 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X2(cd,ij) end do end do @@ -142,7 +142,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do l=k+1,nO kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & - + eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y2(kl,ij) + + 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij) end do end do end do diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index bf66039..fe78983 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -198,8 +198,6 @@ program QuAcK call read_dipole_integrals(working_dir,nBas,dipole_int_AO) call wall_time(end_int) - call matout(nBas,nBas,dipole_int_AO(:,:,1)) - t_int = end_int - start_int write(*,*) write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 1e-integrals = ',t_int,' seconds'