From fda6a452764344189ad077764a1b91bcb6635af7 Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 2 Sep 2024 22:18:53 +0200 Subject: [PATCH] passing Fock matrix and others in R branch --- src/CC/RCC.f90 | 18 ++++++------ src/CC/pCCD.f90 | 69 ++++++++++++++++++++++++++++++-------------- src/HF/RHF.f90 | 59 ++++++++++++++++++------------------- src/HF/ROHF.f90 | 65 ++++++++++++++++++++--------------------- src/QuAcK/RQuAcK.f90 | 23 +++++++++------ 5 files changed, 129 insertions(+), 105 deletions(-) diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index e4c774f..f03fcdd 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -1,8 +1,6 @@ - -! --- - -subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, & - maxSCF, thresh, max_diis, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF) +subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & + maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc, & + ERHF,eHF,cHF,PHF,FHF) ! Coupled-cluster module @@ -27,7 +25,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -36,6 +35,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: FHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) @@ -166,9 +167,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(dopCCD) then call wall_time(start_CC) - call pCCD(dotest, maxSCF, thresh, max_diis, nBas, nOrb, & - nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF) - + call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & + Hc,ERI_AO,ENuc,ERHF,eHF,cHF,PHF,FHF) call wall_time(end_CC) t_CC = end_CC - start_CC diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index e238950..e700295 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -1,8 +1,5 @@ - -! --- - -subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & - nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF) +subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, & + Hc,ERI_AO,ENuc,ERHF,eHF,cHF,PHF,FHF) ! pair CCD module @@ -16,15 +13,23 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas, nOrb, nC, nO, nV, nR + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: ENuc,ERHF double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: FHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) ! Local variables + integer :: mu,nu integer :: p,q,r,s,t,u,w integer :: pq,rs integer :: i,j,a,b @@ -35,6 +40,7 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & double precision :: CvgOrb double precision :: ECC double precision :: EcCC + double precision :: dECC double precision,allocatable :: eO(:) double precision,allocatable :: eV(:) @@ -93,18 +99,18 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & O = nO - nC V = nV - nR - N = O + V ! nOrb - nC - nR + N = O + V !------------------------------------! ! Star Loop for orbital optimization ! !------------------------------------! allocate(ERI_MO(N,N,N,N)) - allocate(c(nBas,N), h(N,N)) - allocate(eO(O), eV(V), delta_OV(O,V)) - allocate(OOOO(O,O), OOVV(O,V), OVOV(O,V), OVVO(O,V), VVVV(V,V)) + allocate(c(nBas,N),h(N,N)) + allocate(eO(O),eV(V),delta_OV(O,V)) + allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V)) - do i = 1, N + do i=1,N c(:,i) = cHF(:,nC+i) enddo @@ -116,20 +122,34 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & write(*,*)'| Orbital Optimization for pCCD |' write(*,*)'----------------------------------------------------' - do while(CvgOrb > thresh .and. nItOrb < 1) + do while(CvgOrb > thresh .and. nItOrb < maxIt) nItOrb = nItOrb + 1 ! Transform integrals - h = matmul(transpose(c), matmul(Hc, c)) + h = matmul(transpose(c),matmul(Hc,c)) - call AOtoMO_ERI_RHF(nBas, N, c(1,1), ERI_AO(1,1,1,1), ERI_MO(1,1,1,1)) + call AOtoMO_ERI_RHF(nBas,N,c(1,1),ERI_AO(1,1,1,1),ERI_MO(1,1,1,1)) ! Form energy denominator - eO(:) = eHF(nC+1:nO) - eV(:) = eHF(nO+1:nOrb-nR) + eO(:) = 0d0 + eV(:) = 0d0 + + do mu=1,nBas + do nu=1,nBas + + do i=1,O + eO(i) = eO(i) + c(mu,i)*FHF(mu,nu)*c(nu,i) + end do + + do a=1,V + eV(a) = eV(a) + c(mu,O+a)*FHF(mu,nu)*c(nu,O+a) + end do + + end do + end do do i=1,O do a=1,V @@ -170,6 +190,7 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & nItAmp = 0 ECC = ERHF EcCC = 0d0 + dECC = ECC n_diis = 0 t2(:,:) = 0d0 @@ -573,9 +594,13 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & ! Check convergence of orbital optimization CvgOrb = maxval(abs(grad)) - write(*,*) ' Iteration',nItOrb,'for pCCD orbital optimization' - write(*,*) ' Convergence of orbital gradient = ',CvgOrb + write(*,'(A10,I4,A30)') ' Iteration',nItOrb,'for pCCD orbital optimization' + write(*,*) '----------------------------------------------------------' + write(*,'(A40,F16.10,A3)') ' Convergence of orbital gradient = ',CvgOrb,' au' + write(*,'(A40,F16.10,A3)') ' Energy difference = ',ECC-dECC,' au' write(*,*) + + dECC = ECC !-------------------------! ! Compute orbital Hessian ! @@ -703,18 +728,18 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, & deallocate(Kap) write(*,*) 'e^kappa' - call matout(N, N, ExpKap) + call matout(N,N,ExpKap) write(*,*) write(*,*) 'Old orbitals' - call matout(nBas, N, c) + call matout(nBas,N,c) write(*,*) - c = matmul(c, ExpKap) + c = matmul(c,ExpKap) deallocate(ExpKap) write(*,*) 'Rotated orbitals' - call matout(nBas, N, c) + call matout(nBas,N,c) write(*,*) end do diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 1ae1b27..e66de73 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,8 +1,5 @@ - -! --- - -subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & - nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P) +subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,eHF,c,P,F) ! Perform restricted Hartree-Fock calculation @@ -19,7 +16,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nO integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -53,7 +51,6 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) double precision,allocatable :: cp(:,:) - double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) ! Output variables @@ -62,6 +59,7 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, double precision,intent(out) :: eHF(nOrb) double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: F(nBas,nBas) ! Hello world @@ -81,7 +79,6 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, allocate(K(nBas,nBas)) allocate(err(nBas,nBas)) - allocate(F(nBas,nBas)) allocate(cp(nOrb,nOrb)) allocate(Fp(nOrb,nOrb)) @@ -91,7 +88,7 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Guess coefficients and density matrix - call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) + call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c) !P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & @@ -105,8 +102,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, err_diis(:,:) = 0d0 rcond = 0d0 - Conv = 1d0 - nSCF = 0 + Conv = 1d0 + nSCF = 0 !------------------------------------------------------------------------ ! Main SCF loop @@ -126,8 +123,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Build Fock matrix - call Hartree_matrix_AO_basis(nBas, P, ERI, J) - call exchange_matrix_AO_basis(nBas, P, ERI, K) + call Hartree_matrix_AO_basis(nBas,P,ERI,J) + call exchange_matrix_AO_basis(nBas,P,ERI,K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) if(nBas .ne. nOrb) then @@ -137,24 +134,24 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Check convergence - err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F) + err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) if(nSCF > 1) Conv = maxval(abs(err)) ! Kinetic energy - ET = trace_matrix(nBas, matmul(P, T)) + ET = trace_matrix(nBas,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas, matmul(P, V)) + EV = trace_matrix(nBas,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) ! Exchange energy - EK = 0.25d0*trace_matrix(nBas, matmul(P, K)) + EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) ! Total energy @@ -164,29 +161,29 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, if(max_diis > 1) then - n_diis = min(n_diis+1, max_diis) - call DIIS_extrapolation(rcond, nBas_Sq, nBas_Sq, n_diis, err_diis, F_diis, err, F) + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) end if ! Level shift if(level_shift > 0d0 .and. Conv > thresh) then - call level_shifting(level_shift, nBas, nOrb, nO, S, c, F) + call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) endif ! Diagonalize Fock matrix if(nBas .eq. nOrb) then - Fp = matmul(transpose(X), matmul(F, X)) + Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eHF) - c = matmul(X, cp) + call diagonalize_matrix(nOrb,cp,eHF) + c = matmul(X,cp) else - Fp = matmul(transpose(c), matmul(F, c)) + Fp = matmul(transpose(c),matmul(F,c)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eHF) - c = matmul(c, cp) + call diagonalize_matrix(nOrb,cp,eHF) + c = matmul(c,cp) endif ! Density matrix @@ -217,7 +214,7 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) + deallocate(J,K,err,cp,Fp,err_diis,F_diis) stop @@ -225,8 +222,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Compute dipole moments - call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int, dipole) - call print_RHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole) + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_RHF(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,ERHF,dipole) ! Testing zone @@ -239,6 +236,6 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, end if - deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) + deallocate(J,K,err,cp,Fp,err_diis,F_diis) end subroutine diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index 5a0ff02..ed73af9 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -1,8 +1,5 @@ - -! --- - -subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, & - nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, EROHF, eHF, c, Ptot) +subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot,Ftot) ! Perform restricted open-shell Hartree-Fock calculation @@ -19,7 +16,8 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh - integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nOrb + integer,intent(in) :: nBas integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -52,7 +50,6 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, double precision,allocatable :: J(:,:,:) double precision,allocatable :: F(:,:,:) double precision,allocatable :: Fp(:,:) - double precision,allocatable :: Ftot(:,:) double precision,allocatable :: P(:,:,:) double precision,allocatable :: K(:,:,:) double precision,allocatable :: err(:,:) @@ -68,6 +65,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, double precision,intent(out) :: eHF(nOrb) double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: Ptot(nBas,nBas) + double precision,intent(out) :: Ftot(nBas,nBas) ! Hello world @@ -86,7 +84,6 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, allocate(J(nBas,nBas,nspin)) allocate(K(nBas,nBas,nspin)) allocate(F(nBas,nBas,nspin)) - allocate(Ftot(nBas,nBas)) allocate(P(nBas,nBas,nspin)) allocate(err(nBas,nBas)) @@ -98,9 +95,9 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Guess coefficients and demsity matrices - call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) + call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c) - do ispin = 1, nspin + do ispin = 1,nspin !P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin)))) call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas) end do @@ -134,51 +131,51 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Build Hartree repulsion - do ispin = 1, nspin - call Hartree_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin)) + do ispin=1,nspin + call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) end do ! Compute exchange potential - do ispin = 1, nspin - call exchange_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin)) + do ispin=1,nspin + call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin)) end do ! Build Fock operator - do ispin = 1, nspin + do ispin=1,nspin F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) end do - call ROHF_fock_matrix(nBas, nOrb, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot) + call ROHF_fock_matrix(nBas,nOrb,nO(1),nO(2),S,c,F(:,:,1),F(:,:,2),Ftot) ! Check convergence - err(:,:) = matmul(Ftot, matmul(Ptot, S)) - matmul(matmul(S, Ptot), Ftot) + err(:,:) = matmul(Ftot,matmul(Ptot,S)) - matmul(matmul(S,Ptot),Ftot) if(nSCF > 1) Conv = maxval(abs(err(:,:))) ! Kinetic energy - do ispin = 1, nspin - ET(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), T(:,:))) + do ispin=1,nspin + ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) end do ! Potential energy - do ispin = 1, nspin - EV(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), V(:,:))) + do ispin=1,nspin + EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) end do ! Hartree energy - EJ(1) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,1), J(:,:,1))) - EJ(2) = trace_matrix(nBas, matmul(P(:,:,1), J(:,:,2))) - EJ(3) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,2), J(:,:,2))) + EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) + EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) + EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) ! Exchange energy - do ispin = 1, nspin - EK(ispin) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,ispin), K(:,:,ispin))) + do ispin = 1,nspin + EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) end do ! Total energy @@ -199,27 +196,27 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(level_shift, nBas, nOrb, maxval(nO), S, c, Ftot) + call level_shifting(level_shift,nBas,nOrb,maxval(nO),S,c,Ftot) end do end if ! Transform Fock matrix in orthogonal basis - Fp(:,:) = matmul(transpose(X(:,:)), matmul(Ftot(:,:), X(:,:))) + Fp(:,:) = matmul(transpose(X(:,:)),matmul(Ftot(:,:),X(:,:))) ! Diagonalize Fock matrix to get eigenvectors and eigenvalues cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb, cp, eHF) + call diagonalize_matrix(nOrb,cp,eHF) ! Back-transform eigenvectors in non-orthogonal basis - c(:,:) = matmul(X(:,:), cp(:,:)) + c(:,:) = matmul(X(:,:),cp(:,:)) ! Compute density matrix - do ispin = 1, nspin + do ispin=1,nspin !P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin)))) call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas) end do @@ -246,7 +243,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) + deallocate(J,K,F,P,err,Fp,cp,err_diis,F_diis) stop @@ -255,7 +252,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Compute final UHF energy call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_ROHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, EROHF, dipole) + call print_ROHF(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole) ! Print test values @@ -265,6 +262,6 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, end if - deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) + deallocate(J,K,F,P,err,Fp,cp,err_diis,F_diis) end subroutine diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 2457011..3a19de4 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -92,7 +92,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d double precision :: start_GW ,end_GW ,t_GW double precision :: start_GT ,end_GT ,t_GT - double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:) + double precision,allocatable :: eHF(:) + double precision,allocatable :: cHF(:,:) + double precision,allocatable :: PHF(:,:) + double precision,allocatable :: FHF(:,:) double precision :: ERHF double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) @@ -109,9 +112,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Memory allocation ! !-------------------! - allocate(cHF(nBas,nOrb)) allocate(eHF(nOrb)) + allocate(cHF(nBas,nOrb)) allocate(PHF(nBas,nBas)) + allocate(FHF(nBas,nBas)) allocate(dipole_int_MO(nOrb,nOrb,ncart)) allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) @@ -122,8 +126,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doRHF) then call wall_time(start_HF) - call RHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & - nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF) + call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -135,8 +139,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doROHF) then call wall_time(start_HF) - call ROHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, & - nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF) + call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -163,7 +167,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! 4-index transform - call AOtoMO_ERI_RHF(nBas, nOrb, cHF, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas,nOrb,cHF,ERI_AO,ERI_MO) call wall_time(end_AOtoMO) @@ -231,8 +235,9 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doCC) then call wall_time(start_CC) - call RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, & - maxSCF_CC, thresh_CC, max_diis_CC, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF) + call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & + maxSCF_CC,thresh_CC,max_diis_CC,nBas,nOrb,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO, & + ENuc,ERHF,eHF,cHF,PHF,FHF) call wall_time(end_CC) t_CC = end_CC - start_CC