diff --git a/PyDuck.py b/PyDuck.py index c74dc22..a8e9033 100644 --- a/PyDuck.py +++ b/PyDuck.py @@ -18,6 +18,7 @@ parser.add_argument('-b', '--basis', type=str, required=True, help='Name of the parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.') parser.add_argument('-c', '--charge', type=int, default=0, help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0') parser.add_argument('--cartesian', default=False, action='store_true', help='Add this option if you want to use cartesian basis functions.') +parser.add_argument('--print_2e', default=False, action='store_true', help='Add this option if you want to print 2e-integrals.') parser.add_argument('-fc', '--frozen_core', type=bool, default=False, help='Freeze core MOs. Default is false') parser.add_argument('-m', '--multiplicity', type=int, default=1, help='Spin multiplicity. Default is 1 therefore singlet') parser.add_argument('--working_dir', type=str, default=QuAcK_dir, help='Set a working directory to run the calculation.') @@ -32,6 +33,7 @@ frozen_core=args.frozen_core multiplicity=args.multiplicity xyz=args.xyz + '.xyz' cartesian=args.cartesian +print_2e=args.print_2e working_dir=args.working_dir #Read molecule @@ -90,11 +92,10 @@ t1e = mol.intor('int1e_kin') #Kinetic energy matrix elements dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators x,y,z = dipole[0],dipole[1],dipole[2] -norb = len(ovlp) +norb = len(ovlp) # nBAS_AOs subprocess.call(['rm', working_dir + '/int/nBas.dat']) f = open(working_dir+'/int/nBas.dat','w') -f.write(str(norb)) -f.write(' ') +f.write(" {} ".format(str(norb))) f.close() @@ -122,7 +123,6 @@ write_matrix_to_file(y,norb,working_dir+'/int/y.dat') subprocess.call(['rm', working_dir + '/int/z.dat']) write_matrix_to_file(z,norb,working_dir+'/int/z.dat') -#Write two-electron integrals eri_ao = mol.intor('int2e') def write_tensor_to_file(tensor,size,file,cutoff=1e-15): @@ -132,12 +132,23 @@ def write_tensor_to_file(tensor,size,file,cutoff=1e-15): for k in range(i,size): for l in range(j,size): if abs(tensor[i][k][j][l]) > cutoff: + #f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l])) f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l])) f.write('\n') f.close() -subprocess.call(['rm', working_dir + '/int/ERI.dat']) -write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat') +# Write two-electron integrals +if print_2e: + # (formatted) + subprocess.call(['rm', working_dir + '/int/ERI.dat']) + write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat') +else: + # (binary) + subprocess.call(['rm', working_dir + '/int/ERI.bin']) + # chem -> phys notation + eri_ao = eri_ao.transpose(0, 2, 1, 3) + eri_ao.tofile('int/ERI.bin') + #Execute the QuAcK fortran program subprocess.call(QuAcK_dir+'/bin/QuAcK') diff --git a/src/AOtoMO/AOtoMO.f90 b/src/AOtoMO/AOtoMO.f90 index 2192bc7..8383273 100644 --- a/src/AOtoMO/AOtoMO.f90 +++ b/src/AOtoMO/AOtoMO.f90 @@ -1,26 +1,31 @@ -subroutine AOtoMO(nBas,C,A,B) +subroutine AOtoMO(nBas, nOrb, C, M_AOs, M_MOs) -! Perform AO to MO transformation of a matrix A for given coefficients c + ! Perform AO to MO transformation of a matrix M_AOs for given coefficients c + ! M_MOs = C.T M_AOs C implicit none -! Input variables + integer, intent(in) :: nBas, nOrb + double precision, intent(in) :: C(nBas,nOrb) + double precision, intent(in) :: M_AOs(nBas,nBas) - integer,intent(in) :: nBas - double precision,intent(in) :: C(nBas,nBas) - double precision,intent(in) :: A(nBas,nBas) + double precision, intent(out) :: M_MOs(nOrb,nOrb) -! Local variables + double precision, allocatable :: AC(:,:) - double precision,allocatable :: AC(:,:) + allocate(AC(nBas,nOrb)) -! Output variables + !AC = matmul(M_AOs, C) + !M_MOs = matmul(transpose(C), AC) - double precision,intent(out) :: B(nBas,nBas) + call dgemm("N", "N", nBas, nOrb, nBas, 1.d0, & + M_AOs(1,1), nBas, C(1,1), nBas, & + 0.d0, AC(1,1), nBas) - allocate(AC(nBas,nBas)) + call dgemm("T", "N", nOrb, nOrb, nBas, 1.d0, & + C(1,1), nBas, AC(1,1), nBas, & + 0.d0, M_MOs(1,1), nOrb) - AC = matmul(A,C) - B = matmul(transpose(C),AC) + deallocate(AC) end subroutine diff --git a/src/AOtoMO/AOtoMO_ERI_RHF.f90 b/src/AOtoMO/AOtoMO_ERI_RHF.f90 index c93e111..f9f64a2 100644 --- a/src/AOtoMO/AOtoMO_ERI_RHF.f90 +++ b/src/AOtoMO/AOtoMO_ERI_RHF.f90 @@ -1,4 +1,7 @@ -subroutine AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + +! --- + +subroutine AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm @@ -7,32 +10,51 @@ subroutine AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: c(nBas,nOrb) ! Local variables - double precision,allocatable :: scr(:,:,:,:) - integer :: mu,nu,la,si - integer :: i,j,k,l + double precision,allocatable :: a1(:,:,:,:) + double precision,allocatable :: a2(:,:,:,:) ! Output variables - double precision,intent(out) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(out) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) ! Memory allocation - allocate(scr(nBas,nBas,nBas,nBas)) + allocate(a2(nBas,nBas,nBas,nOrb)) + allocate(a1(nBas,nBas,nOrb,nOrb)) ! Four-index transform via semi-direct O(N^5) algorithm - call dgemm('T','N',nBas**3,nBas,nBas,1d0,ERI_AO,nBas,c(1,1),size(c,1),0d0,scr,nBas**3) - - call dgemm('T','N',nBas**3,nBas,nBas,1d0,scr,nBas,c(1,1),size(c,1),0d0,ERI_MO,nBas**3) + call dgemm( 'T', 'N', nBas*nBas*nBas, nOrb, nBas, 1.d0 & + , ERI_AO(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, a2(1,1,1,1), nBas*nBas*nBas) - call dgemm('T','N',nBas**3,nBas,nBas,1d0,ERI_MO,nBas,c(1,1),size(c,1),0d0,scr,nBas**3) + call dgemm( 'T', 'N', nBas*nBas*nOrb, nOrb, nBas, 1.d0 & + , a2(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, a1(1,1,1,1), nBas*nBas*nOrb) - call dgemm('T','N',nBas**3,nBas,nBas,1d0,scr,nBas,c(1,1),size(c,1),0d0,ERI_MO,nBas**3) + deallocate(a2) + allocate(a2(nBas,nOrb,nOrb,nOrb)) + + call dgemm( 'T', 'N', nBas*nOrb*nOrb, nOrb, nBas, 1.d0 & + , a1(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, a2(1,1,1,1), nBas*nOrb*nOrb) + + deallocate(a1) + + call dgemm( 'T', 'N', nOrb*nOrb*nOrb, nOrb, nBas, 1.d0 & + , a2(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, ERI_MO(1,1,1,1), nOrb*nOrb*nOrb) + + deallocate(a2) end subroutine + + + + diff --git a/src/AOtoMO/MOtoAO.f90 b/src/AOtoMO/MOtoAO.f90 index faffa92..a5ffaed 100644 --- a/src/AOtoMO/MOtoAO.f90 +++ b/src/AOtoMO/MOtoAO.f90 @@ -1,31 +1,39 @@ -subroutine MOtoAO(nBas,S,C,B,A) +subroutine MOtoAO(nBas, nOrb, S, C, M_MOs, M_AOs) -! Perform MO to AO transformation of a matrix A for a given metric S -! and coefficients c + ! Perform MO to AO transformation of a matrix M_AOs for a given metric S + ! and coefficients c + ! + ! M_AOs = S C M_MOs (S C).T implicit none -! Input variables + integer, intent(in) :: nBas, nOrb + double precision, intent(in) :: S(nBas,nBas) + double precision, intent(in) :: C(nBas,nOrb) + double precision, intent(in) :: M_MOs(nOrb,nOrb) + double precision, intent(out) :: M_AOs(nBas,nBas) - integer,intent(in) :: nBas - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: C(nBas,nBas) - double precision,intent(in) :: B(nBas,nBas) + double precision, allocatable :: SC(:,:),BSC(:,:) -! Local variables - double precision,allocatable :: SC(:,:),BSC(:,:) + allocate(SC(nBas,nOrb), BSC(nOrb,nBas)) -! Output variables + !SC = matmul(S, C) + !BSC = matmul(M_MOs, transpose(SC)) + !M_AOs = matmul(SC, BSC) - double precision,intent(out) :: A(nBas,nBas) + call dgemm("N", "N", nBas, nOrb, nBas, 1.d0, & + S(1,1), nBas, C(1,1), nBas, & + 0.d0, SC(1,1), nBas) -! Memory allocation + call dgemm("N", "T", nOrb, nBas, nOrb, 1.d0, & + M_MOs(1,1), nOrb, SC(1,1), nBas, & + 0.d0, BSC(1,1), nOrb) - allocate(SC(nBas,nBas),BSC(nBas,nBas)) + call dgemm("N", "N", nBas, nBas, nOrb, 1.d0, & + SC(1,1), nBas, BSC(1,1), nOrb, & + 0.d0, M_AOs(1,1), nBas) - SC = matmul(S,C) - BSC = matmul(B,transpose(SC)) - A = matmul(SC,BSC) + deallocate(SC, BSC) end subroutine diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index aa8923b..e4c774f 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -1,5 +1,8 @@ -subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF,thresh,max_diis,nBas,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) ! Coupled-cluster module @@ -24,18 +27,18 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) ! Local variables @@ -48,11 +51,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCD) then call wall_time(start_CC) - call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call CCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CCD = ',t_CC,' seconds' write(*,*) end if @@ -64,12 +67,12 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doDCD) then call wall_time(start_CC) - call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & + call DCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for DCD = ',t_CC,' seconds' write(*,*) end if @@ -83,11 +86,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCSD) then call wall_time(start_CC) - call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CCSD or CCSD(T)= ',t_CC,' seconds' write(*,*) end if @@ -99,11 +102,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dodrCCD) then call wall_time(start_CC) - call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call drCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for direct ring CCD = ',t_CC,' seconds' write(*,*) end if @@ -115,11 +118,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dorCCD) then call wall_time(start_CC) - call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call rCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for rCCD = ',t_CC,' seconds' write(*,*) end if @@ -131,11 +134,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(docrCCD) then call wall_time(start_CC) - call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call crCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for crossed-ring CCD = ',t_CC,' seconds' write(*,*) end if @@ -147,11 +150,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dolCCD) then call wall_time(start_CC) - call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call lCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ladder CCD = ',t_CC,' seconds' write(*,*) end if @@ -163,12 +166,13 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dopCCD) then call wall_time(start_CC) - call pCCD(dotest,maxSCF,thresh,max_diis,nBas,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) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pair CCD = ',t_CC,' seconds' write(*,*) end if diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index bcfd50b..e238950 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -1,4 +1,8 @@ -subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,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) ! pair CCD module @@ -12,15 +16,10 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + integer,intent(in) :: nBas, nOrb, nC, nO, nV, nR + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) @@ -94,18 +93,20 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH O = nO - nC V = nV - nR - N = O + V + N = O + V ! nOrb - nC - nR !------------------------------------! ! Star Loop for orbital optimization ! !------------------------------------! allocate(ERI_MO(N,N,N,N)) - allocate(c(N,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)) - c(:,:) = cHF(nC+1:nBas-nR,nC+1:nBas-nR) + do i = 1, N + c(:,i) = cHF(:,nC+i) + enddo CvgOrb = 1d0 nItOrb = 0 @@ -121,13 +122,14 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH ! Transform integrals - h = matmul(transpose(c),matmul(Hc(nC+1:nBas-nR,nC+1:nBas-nR),c)) - call AOtoMO_ERI_RHF(N,c,ERI_AO(nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR),ERI_MO) + 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)) ! Form energy denominator eO(:) = eHF(nC+1:nO) - eV(:) = eHF(nO+1:nBas-nR) + eV(:) = eHF(nO+1:nOrb-nR) do i=1,O do a=1,V @@ -701,18 +703,18 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH deallocate(Kap) write(*,*) 'e^kappa' - call matout(N,N,ExpKap) + call matout(N, N, ExpKap) write(*,*) write(*,*) 'Old orbitals' - call matout(N,N,c) + call matout(nBas, N, c) write(*,*) - c = matmul(c,ExpKap) + c = matmul(c, ExpKap) deallocate(ExpKap) write(*,*) 'Rotated orbitals' - call matout(N,N,c) + call matout(nBas, N, c) write(*,*) end do diff --git a/src/CI/RCI.f90 b/src/CI/RCI.f90 index 3762baf..3a9905e 100644 --- a/src/CI/RCI.f90 +++ b/src/CI/RCI.f90 @@ -1,5 +1,8 @@ -subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & - epsHF,EHF,cHF,S) + +! --- + +subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, & + nC, nO, nV, nR, nS, ERI, dipole_int, epsHF, EHF) ! Configuration interaction module @@ -18,18 +21,16 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n logical,intent(in) :: singlet logical,intent(in) :: triplet - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: epsHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables @@ -42,11 +43,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n if(doCIS) then call wall_time(start_CI) - call RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) + call RCIS(dotest,singlet,triplet,doCIS_D,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) call wall_time(end_CI) t_CI = end_CI - start_CI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CIS = ',t_CI,' seconds' write(*,*) end if @@ -58,11 +59,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n if(doCID) then call wall_time(start_CI) - call CID(dotest,singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF) + call CID(dotest,singlet,triplet,nOrb,nC,nO,nV,nR,ERI,epsHF,EHF) call wall_time(end_CI) t_CI = end_CI - start_CI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CID = ',t_CI,' seconds' write(*,*) end if @@ -74,11 +75,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n if(doCISD) then call wall_time(start_CI) - call CISD(dotest,singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF) + call CISD(dotest,singlet,triplet,nOrb,nC,nO,nV,nR,ERI,epsHF,EHF) call wall_time(end_CI) t_CI = end_CI - start_CI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CISD = ',t_CI,' seconds' write(*,*) end if @@ -91,11 +92,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n call wall_time(start_CI) write(*,*) ' FCI is not yet implemented! Sorry.' -! call FCI(nBas,nC,nO,nV,nR,ERI,epsHF) +! call FCI(nOrb,nC,nO,nV,nR,ERI,epsHF) call wall_time(end_CI) t_CI = end_CI - start_CI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for FCI = ',t_CI,' seconds' write(*,*) end if diff --git a/src/CI/RCIS.f90 b/src/CI/RCIS.f90 index 84489e3..aed3196 100644 --- a/src/CI/RCIS.f90 +++ b/src/CI/RCIS.f90 @@ -41,7 +41,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in ! Memory allocation - allocate(A(nS,nS),Om(nS)) + allocate(A(nS,nS), Om(nS)) ! Compute CIS matrix @@ -117,4 +117,6 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in end if + deallocate(A, Om) + end subroutine diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90 index 7cb3c40..53aeb1e 100644 --- a/src/GF/RG0F2.f90 +++ b/src/GF/RG0F2.f90 @@ -51,7 +51,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGFlin(nBas),eGF(nBas)) + allocate(SigC(nBas), Z(nBas), eGFlin(nBas), eGF(nBas)) ! Frequency-dependent second-order contribution @@ -133,4 +133,6 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, end if + deallocate(SigC, Z, eGFlin, eGF) + end subroutine diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 2fe50f2..0f1deaf 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -1,7 +1,10 @@ -subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, & - dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & - nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, & - dipole_int_AO,dipole_int,PHF,cHF,epsHF) + +! --- + +subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm, maxSCF, & + thresh, max_diis, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, linearize, & + eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, & + S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF) ! Green's function module @@ -39,7 +42,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -47,18 +50,18 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max integer,intent(in) :: nS double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: epsHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -71,12 +74,13 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doG0F2) then call wall_time(start_GF) - call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & - nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF) + call RG0F2(dotest, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, & + linearize, eta, regularize, nOrb, nC, nO, nV, nR, nS, & + ENuc, EHF, ERI_MO, dipole_int_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -89,12 +93,12 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max call wall_time(start_GF) call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & - ERI,dipole_int,epsHF) + singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EHF, & + ERI_MO,dipole_int_MO,epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -106,12 +110,14 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doqsGF2) then call wall_time(start_GF) - call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF) + call qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & + dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, & + rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, S, & + X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGF2 = ',t_GF,' seconds' write(*,*) end if @@ -123,11 +129,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doufG0F02) then call wall_time(start_GF) - call ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF) + call ufRG0F02(dotest, nOrb, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0F02 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0F02 = ',t_GF,' seconds' write(*,*) end if @@ -139,11 +145,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doG0F3) then call wall_time(start_GF) - call RG0F3(dotest,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) + call RG0F3(dotest, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds' write(*,*) end if @@ -155,11 +161,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max if(doevGF3) then call wall_time(start_GF) - call evRGF3(dotest,maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF) + call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds' write(*,*) end if diff --git a/src/GF/evRGF2.f90 b/src/GF/evRGF2.f90 index 4dc2610..c3c4637 100644 --- a/src/GF/evRGF2.f90 +++ b/src/GF/evRGF2.f90 @@ -62,7 +62,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + allocate(SigC(nBas), Z(nBas), eGF(nBas), eOld(nBas), error_diis(nBas,max_diis), e_diis(nBas,max_diis)) ! Initialization @@ -189,4 +189,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si end if + deallocate(SigC, Z, eGF, eOld, error_diis, e_diis) + end subroutine diff --git a/src/GF/print_qsRGF2.f90 b/src/GF/print_qsRGF2.f90 index 845c5a0..feb590f 100644 --- a/src/GF/print_qsRGF2.f90 +++ b/src/GF/print_qsRGF2.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF,dipole) + +! --- + +subroutine print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, c, & + SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF, dipole) ! Print one-electron energies and other stuff for qsGF2 @@ -7,17 +11,17 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGF(nBas) - double precision,intent(in) :: c(nBas) - double precision,intent(in) :: SigC(nBas,nBas) - double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGF(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: ET double precision,intent(in) :: EV double precision,intent(in) :: EJ @@ -53,7 +57,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do q=1,nBas + do q = 1, nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|' end do @@ -102,12 +106,12 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ, write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGF2 MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas, nOrb, c) write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGF2 MO energies' write(*,'(A50)') '---------------------------------------' - call matout(nBas,1,eGF) + call matout(nOrb, 1, eGF) write(*,*) end if diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 8222e56..651c876 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -1,6 +1,10 @@ -subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & + dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, & + rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, & + S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GF2 calculation @@ -29,25 +33,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas,nOrb,nC,nO,nV,nR,nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: ispin integer :: n_diis double precision :: EqsGF2 @@ -94,7 +98,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Stuff - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! TDA @@ -105,9 +109,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Memory allocation - allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), & - error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGF(nOrb)) + allocate(eOld(nOrb)) + + allocate(c(nBas,nOrb)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(error(nBas,nBas)) + + allocate(Z(nOrb)) + allocate(SigC(nOrb,nOrb)) + + allocate(SigCp(nBas,nBas)) + + allocate(error_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Initialization @@ -117,7 +139,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si Conv = 1d0 P(:,:) = PHF(:,:) eOld(:) = eHF(:) - eGF(:) = eHF(:) + eGF(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 @@ -135,25 +157,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Buid Hartree matrix - call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + call exchange_matrix_AO_basis(nBas, P, ERI_AO, K) ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute self-energy and renormalization factor if(regularize) then - call GF2_reg_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z) + call GF2_reg_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) else - call GF2_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z) + call GF2_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) end if @@ -161,7 +183,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas,S,c,SigC,SigCp) + call MOtoAO(nBas, nOrb, S, c, SigC, SigCp) ! Solve the quasi-particle equation @@ -169,28 +191,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Compute commutator and convergence criteria - error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + error = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F) ! DIIS extrapolation - n_diis = min(n_diis+1,max_diis) + n_diis = min(n_diis+1, max_diis) if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) else n_diis = 0 end if ! Diagonalize Hamiltonian in AO basis - Fp = matmul(transpose(X),matmul(F,X)) + Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGF) - c = matmul(X,cp) - SigCp = matmul(transpose(c),matmul(SigCp,c)) + call diagonalize_matrix(nOrb, cp, eGF) + c = matmul(X, cp) ! Compute new density matrix in the AO basis - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) ! Save quasiparticles energy for next cycle @@ -203,23 +224,23 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! 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 - Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) + Ex = 0.25d0*trace_matrix(nBas, matmul(P, K)) ! Correlation energy - call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec) + call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec) ! Total energy @@ -230,8 +251,9 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si ! Print results !------------------------------------------------------------------------ - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole) + call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) + call print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, & + c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole) end do !------------------------------------------------------------------------ @@ -248,19 +270,21 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis) stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis) ! Perform BSE calculation if(dophBSE) then - call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE) + call GF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, & + nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -278,7 +302,8 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si if(doppBSE) then - call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE) + call GF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & + nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 6fa0bdb..2aedac1 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -197,7 +197,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved end do do is=1,nspin - call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 index 4225d36..74d1d12 100644 --- a/src/GT/RGT.f90 +++ b/src/GT/RGT.f90 @@ -1,7 +1,11 @@ -subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis, & - doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF) + +! --- + +subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & + maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & + doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, & + nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, & + V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! T-matrix module @@ -44,7 +48,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -52,18 +56,18 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -78,11 +82,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0pp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0pp = ',t_GT,' seconds' write(*,*) end if @@ -95,11 +99,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTpp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTpp = ',t_GT,' seconds' write(*,*) end if @@ -111,13 +115,13 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do if(doqsGTpp) then call wall_time(start_GT) - call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, & - PHF,cHF,eHF) + call qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, dBSE, & + dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, & + nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTpp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTpp = ',t_GT,' seconds' write(*,*) end if @@ -129,11 +133,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do if(doufG0T0pp) then call wall_time(start_GT) - call ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0T0pp(dotest,TDA_T,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0T0pp = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0T0pp = ',t_GT,' seconds' write(*,*) end if @@ -146,11 +150,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0eh = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0eh = ',t_GT,' seconds' write(*,*) end if @@ -163,11 +167,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do call wall_time(start_GT) call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTeh = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTeh = ',t_GT,' seconds' write(*,*) end if @@ -179,13 +183,14 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do if(doqsGTeh) then call wall_time(start_GT) - call qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, & - PHF,cHF,eHF) + call qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, & + dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, & + ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & + Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTeh = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTeh = ',t_GT,' seconds' write(*,*) end if diff --git a/src/GT/print_qsRGTeh.f90 b/src/GT/print_qsRGTeh.f90 index 9f8f266..e1e4f95 100644 --- a/src/GT/print_qsRGTeh.f90 +++ b/src/GT/print_qsRGTeh.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + +! --- + +subroutine print_qsRGTeh(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, & + Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) ! Print one-electron energies and other stuff for qsGTeh @@ -7,7 +11,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -19,11 +23,11 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGT(nBas) - double precision,intent(in) :: c(nBas) - double precision,intent(in) :: SigC(nBas,nBas) - double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGT(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: EqsGT double precision,intent(in) :: dipole(ncart) @@ -58,7 +62,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ '|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas + do p=1,nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' end do @@ -109,13 +113,13 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTeh MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTeh MO energies' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eGT) + call vecout(nOrb, eGT) write(*,*) end if diff --git a/src/GT/print_qsRGTpp.f90 b/src/GT/print_qsRGTpp.f90 index c3bebfa..8eab4c5 100644 --- a/src/GT/print_qsRGTpp.f90 +++ b/src/GT/print_qsRGTpp.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + +! --- + +subroutine print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, Z, & + ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) ! Print one-electron energies and other stuff for qsGT @@ -7,7 +11,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -19,11 +23,11 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGT(nBas) - double precision,intent(in) :: c(nBas) - double precision,intent(in) :: SigC(nBas,nBas) - double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGT(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: EqsGT double precision,intent(in) :: dipole(ncart) @@ -58,7 +62,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas + do p=1,nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' end do @@ -109,13 +113,13 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTpp MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTpp MO energies' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eGT) + call vecout(nOrb, eGT) write(*,*) end if diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90 index b814001..5f6acab 100644 --- a/src/GT/qsRGTeh.f90 +++ b/src/GT/qsRGTeh.f90 @@ -1,6 +1,10 @@ -subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA, & - dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, & + dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, & + ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & + Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GTeh calculation @@ -33,31 +37,31 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables logical :: dRPA = .false. integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: ispin integer :: n_diis double precision :: ET @@ -113,7 +117,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Stuff - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! TDA for T @@ -131,9 +135,29 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Memory allocation - allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), & - rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(Aph(nS,nS), Bph(nS,nS), Om(nS), XpY(nS,nS), XmY(nS,nS)) + + allocate(eGT(nOrb)) + allocate(eOld(nOrb)) + allocate(Z(nOrb)) + + allocate(c(nBas,nOrb)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(Sig(nOrb,nOrb)) + + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(Sigp(nBas,nBas)) + allocate(err(nBas,nBas)) + + allocate(err_diis(nBas_Sq,max_diis), F_diis(nBas_Sq,max_diis)) + + allocate(rhoL(nOrb,nOrb,nS), rhoR(nOrb,nOrb,nS)) + ! Initialization @@ -169,12 +193,12 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute linear response - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) - if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) + if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -182,17 +206,17 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Compute correlation part of the self-energy - call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) + call GTeh_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) - if(regularize) call GTeh_regularization(nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) + if(regularize) call GTeh_regularization(nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) - call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) + call GTeh_self_energy(eta,nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas,S,c,Sig,Sigp) + call MOtoAO(nBas, nOrb, S, c, Sig, Sigp) ! Solve the quasi-particle equation @@ -207,7 +231,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) end if @@ -215,9 +239,8 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGT) + call diagonalize_matrix(nOrb, cp, eGT) c = matmul(X,cp) - Sigp = matmul(transpose(c),matmul(Sigp,c)) ! Compute new density matrix in the AO basis @@ -255,7 +278,8 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + call print_qsRGTeh(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, Sig, & + Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) end do !------------------------------------------------------------------------ @@ -272,13 +296,15 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis) + stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,err,err_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis) ! Perform BSE calculation diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 2ff7165..15a29dc 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -1,6 +1,9 @@ -subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, & - dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, & + dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, & + nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -31,25 +34,26 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nC,nO,nV,nR,nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: ispin integer :: iblock integer :: n_diis @@ -119,7 +123,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Stuff - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! TDA for T @@ -137,16 +141,30 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Memory allocation - allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas), & - error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGT(nOrb)) + allocate(eOld(nOrb)) + allocate(Z(nOrb)) - allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & - Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & - Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & - Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + allocate(c(nBas,nOrb)) + + allocate(Fp(nOrb,nOrb)) + allocate(cp(nOrb,nOrb)) + allocate(Sig(nOrb,nOrb)) + + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(error(nBas,nBas)) + allocate(Sigp(nBas,nBas)) + + allocate(error_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) + + allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nOrb,nOrb,nVVs)) + allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nOrb,nOrb,nOOs)) + allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nOrb,nOrb,nVVt)) + allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nOrb,nOrb,nOOt)) ! Initialization @@ -182,7 +200,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute linear response @@ -191,9 +209,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -204,9 +222,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) @@ -217,25 +235,25 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Compute correlation part of the self-energy - iblock = 3 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + iblock = 3 + call GTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - iblock = 4 - call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + iblock = 4 + call GTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) if(regularize) then - call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) - call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) + call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) + call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) end if - call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & + call GTpp_self_energy(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas,S,c,Sig,Sigp) + call MOtoAO(nBas, nOrb, S, c, Sig, Sigp) ! Solve the quasi-particle equation @@ -249,22 +267,21 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d n_diis = min(n_diis+1,max_diis) if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) else n_diis = 0 end if ! Diagonalize Hamiltonian in AO basis - Fp = matmul(transpose(X),matmul(F,X)) + Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGT) - c = matmul(X,cp) - Sigp = matmul(transpose(c),matmul(Sigp,c)) + call diagonalize_matrix(nOrb, cp, eGT) + c = matmul(X, cp) ! Compute new density matrix in the AO basis - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) ! Save quasiparticles energy for next cycle @@ -298,7 +315,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + call print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, & + eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, & + EqsGT, dipole) end do !------------------------------------------------------------------------ @@ -315,19 +334,24 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) + deallocate(Om1s, X1s, Y1s, rho1s) + deallocate(Om2s, X2s, Y2s, rho2s) + deallocate(Om1t, X1t, Y1t, rho1t) + deallocate(Om2t, X2t, Y2t, rho2t) stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis) ! Perform BSE calculation if(dophBSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) @@ -363,7 +387,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d end if - call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE) @@ -391,4 +415,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d end if + deallocate(Om1s, X1s, Y1s, rho1s) + deallocate(Om2s, X2s, Y2s, rho2s) + deallocate(Om1t, X1t, Y1t, rho1t) + deallocate(Om2t, X2t, Y2t, rho2t) + end subroutine diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90 index 5f539bb..81c1cda 100644 --- a/src/GT/qsUGTpp.f90 +++ b/src/GT/qsUGTpp.f90 @@ -280,7 +280,7 @@ subroutine qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B end do do ispin=1,nspin - call MOtoAO(nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin)) + call MOtoAO(nBas,nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin)) end do ! Solve the quasi-particle equation diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index 696039d..48408f7 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -1,7 +1,10 @@ -subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF) + +! --- + +subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF, thresh, max_diis, doACFDT, & + exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, & + linearize, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, & + S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Restricted GW module @@ -43,7 +46,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -51,18 +54,18 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -76,11 +79,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0W0 = ',t_GW,' seconds' write(*,*) end if @@ -93,11 +96,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre call wall_time(start_GW) call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGW = ',t_GW,' seconds' write(*,*) end if @@ -109,13 +112,14 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doqsGW) then call wall_time(start_GW) - call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & - dipole_int_AO,dipole_int,PHF,cHF,eHF) + call qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & + TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, & + ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, & + dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' write(*,*) end if @@ -127,13 +131,15 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doSRGqsGW) then call wall_time(start_GW) - call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, & - singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & - dipole_int_AO,dipole_int,PHF,cHF,eHF) + call SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, & + dophBSE, dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & + nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, & + ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, & + PHF, cHF, eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' write(*,*) end if @@ -145,11 +151,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doufG0W0) then call wall_time(start_GW) + ! TODO call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0W0 = ',t_GW,' seconds' write(*,*) end if @@ -161,11 +168,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre if(doufGW) then call wall_time(start_GW) + ! TODO call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufGW = ',t_GW,' seconds' write(*,*) end if diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 9fbbd8c..807991b 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -1,6 +1,10 @@ -subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA, & - dBSE,dTDA,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, & + BSE, BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nNuc, & + ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, & + X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -32,30 +36,30 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: ispin integer :: ixyz integer :: n_diis @@ -114,7 +118,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Stuff - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! TDA for W @@ -132,9 +136,32 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Memory allocation - allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS), & - Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGW(nOrb)) + allocate(eOld(nOrb)) + allocate(Z(nOrb)) + + allocate(c(nBas,nOrb)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(SigC(nOrb,nOrb)) + + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(error(nBas,nBas)) + allocate(SigCp(nBas,nBas)) + + allocate(Aph(nS,nS)) + allocate(Bph(nS,nS)) + allocate(Om(nS)) + allocate(XpY(nS,nS)) + allocate(XmY(nS,nS)) + allocate(rho(nOrb,nOrb,nS)) + + allocate(error_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Initialization @@ -174,22 +201,22 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call wall_time(tao1) - do ixyz=1,ncart - call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas, nOrb, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) call wall_time(tao2) - tao = tao + tao2 -tao1 + tao = tao + tao2 - tao1 ! Compute linear response call wall_time(tlr1) - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -203,13 +230,13 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call wall_time(tex1) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) + call GW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) call wall_time(tex2) tex=tex+tex2-tex1 call wall_time(tsrg1) - call SRG_self_energy(flow,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call SRG_self_energy(flow,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) call wall_time(tsrg2) @@ -218,7 +245,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Make correlation self-energy Hermitian and transform it back to AO basis call wall_time(tmo1) - call MOtoAO(nBas,S,c,SigC,SigCp) + call MOtoAO(nBas, nOrb, S, c, SigC, SigCp) call wall_time(tmo2) tmo = tmo + tmo2 - tmo1 ! Solve the quasi-particle equation @@ -234,18 +261,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) end if ! Diagonalize Hamiltonian in AO basis - Fp = matmul(transpose(X),matmul(F,X)) + Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGW) - c = matmul(X,cp) + call diagonalize_matrix(nOrb, cp, eGW) + c = matmul(X, cp) - call AOtoMO(nBas,c,SigCp,SigC) + call AOtoMO(nBas, nOrb, c, SigCp, SigC) ! Compute new density matrix in the AO basis @@ -283,7 +310,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) + call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, & + SigC, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGW, dipole) end do !------------------------------------------------------------------------ @@ -300,6 +328,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis) + stop end if @@ -313,17 +343,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Cumulant expansion - call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) + call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis) ! Perform BSE calculation if(BSE) then - call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & + nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) if(exchange_kernel) then @@ -357,7 +388,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) + call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, BSE, singlet, triplet, & + eta, nOrb, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/SRG_qsUGW.f90 b/src/GW/SRG_qsUGW.f90 index 0705847..6939afb 100644 --- a/src/GW/SRG_qsUGW.f90 +++ b/src/GW/SRG_qsUGW.f90 @@ -184,8 +184,8 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS !-------------------------------------------------- do ixyz=1,ncart - call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block @@ -228,7 +228,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Make correlation self-energy Hermitian and transform it back to AO basis do is=1,nspin - call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation @@ -279,7 +279,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Back-transform self-energy do is=1,nspin - call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) + call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) end do ! Compute density matrix diff --git a/src/GW/print_qsRGW.f90 b/src/GW/print_qsRGW.f90 index d09a670..9fd695c 100644 --- a/src/GW/print_qsRGW.f90 +++ b/src/GW/print_qsRGW.f90 @@ -1,4 +1,8 @@ -subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) + +! --- + +subroutine print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, & + Z, ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) ! Print useful information about qsRGW calculation @@ -7,7 +11,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -19,11 +23,11 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E double precision,intent(in) :: EcRPA double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: c(nBas) - double precision,intent(in) :: SigC(nBas,nBas) - double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: EqsGW double precision,intent(in) :: dipole(ncart) @@ -59,7 +63,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E '|','#','|','e_HF (eV)','|','Sig_GW (eV)','|','Z','|','e_GW (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas + do p=1,nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' end do @@ -110,13 +114,13 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' Restricted qsGW orbital coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' Restricted qsGW orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eGW) + call vecout(nOrb, eGW) write(*,*) end if diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 6e10f75..1756a0b 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -1,6 +1,10 @@ -subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, & - ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! --- + +subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & + TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, & + ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, & + ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -34,15 +38,15 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) @@ -50,14 +54,14 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: ispin integer :: ixyz integer :: n_diis @@ -112,7 +116,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Stuff - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! TDA for W @@ -130,10 +134,31 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Memory allocation - allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & - J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), & - Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & - err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(eGW(nOrb)) + allocate(Z(nOrb)) + + allocate(c(nBas,nOrb)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(SigC(nOrb,nOrb)) + + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(err(nBas,nBas)) + allocate(SigCp(nBas,nBas)) + + allocate(Aph(nS,nS)) + allocate(Bph(nS,nS)) + allocate(Om(nS)) + allocate(XpY(nS,nS)) + allocate(XmY(nS,nS)) + allocate(rho(nOrb,nOrb,nS)) + + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Initialization @@ -160,38 +185,38 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Build Hartree-exchange matrix - call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) - call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J) + call exchange_matrix_AO_basis(nBas, P, ERI_AO, K) ! AO to MO transformation of two-electron integrals - do ixyz=1,ncart - call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute linear response - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, eGW, ERI_MO, Aph) + if(.not.TDA_W) call phLR_B(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phLR(TDA_W, nS, Aph, Bph, EcRPA, Om, XpY, XmY) if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) ! Compute correlation part of the self-energy - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) + call GW_excitation_density(nOrb, nC, nO, nR, nS, ERI_MO, XpY, rho) - if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho) + if(regularize) call GW_regularization(nOrb, nC, nO, nV, nR, nS, eGW, Om, rho) - call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call GW_self_energy(eta, nOrb, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z) ! Make correlation self-energy Hermitian and transform it back to AO basis SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas,S,c,SigC,SigCp) + call MOtoAO(nBas, nOrb, S, c, SigC, SigCp) ! Solve the quasi-particle equation @@ -205,19 +230,19 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! 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 @@ -228,26 +253,27 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) end if ! Diagonalize Hamiltonian in AO basis - Fp = matmul(transpose(X),matmul(F,X)) + Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGW) - c = matmul(X,cp) - call AOtoMO(nBas,c,SigCp,SigC) + call diagonalize_matrix(nOrb, cp, eGW) + c = matmul(X, cp) + call AOtoMO(nBas, nOrb, c, SigCp, SigC) ! Density matrix - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) ! Print results - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole) + call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) + call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, & + ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) end do !------------------------------------------------------------------------ @@ -264,19 +290,21 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis) stop end if ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) + deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis) ! Perform BSE calculation if(dophBSE) then - call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & + nOrb, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) if(exchange_kernel) then @@ -310,7 +338,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) + call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, dophBSE, singlet, triplet, & + eta, nOrb, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -327,7 +356,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) + call GW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & + nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, eGW, EcBSE) EcBSE(2) = 3d0*EcBSE(2) diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index ea57d9e..69e3f11 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -184,8 +184,8 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE !-------------------------------------------------- do ixyz=1,ncart - call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block @@ -232,7 +232,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE end do do is=1,nspin - call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation @@ -283,7 +283,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Back-transform self-energy do is=1,nspin - call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) + call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) end do ! Compute density matrix diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index 60fd9d2..5659a04 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -123,7 +123,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Guess coefficients and density matrices - call mo_guess(nBas2,guess_type,S,H,X,C) + call mo_guess(nBas2,nBas2,guess_type,S,H,X,C) ! Construct super density matrix @@ -227,7 +227,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Level-shifting - if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,C,F) + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,C,F) ! Transform Fock matrix in orthogonal basis diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 8b1e441..1ae1b27 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,5 +1,8 @@ -subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,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) ! Perform restricted Hartree-Fock calculation @@ -16,7 +19,7 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -26,14 +29,14 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: n_diis double precision :: ET double precision :: EV @@ -56,8 +59,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Output variables double precision,intent(out) :: ERHF - double precision,intent(out) :: eHF(nBas) - double precision,intent(inout):: c(nBas,nBas) + double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: P(nBas,nBas) ! Hello world @@ -70,24 +73,37 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Useful quantities - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! Memory allocation - allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), & - Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + + allocate(err(nBas,nBas)) + allocate(F(nBas,nBas)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Guess coefficients and density matrix - call mo_guess(nBas,guess_type,S,Hc,X,c) - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + 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, & + c(1,1), nBas, c(1,1), nBas, & + 0.d0, P(1,1), nBas) ! Initialization - n_diis = 0 + n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 - rcond = 0d0 + rcond = 0d0 Conv = 1d0 nSCF = 0 @@ -110,31 +126,35 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! 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 + call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) + call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) + endif ! 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 @@ -144,25 +164,37 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN if(max_diis > 1) then - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,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) call level_shifting(level_shift,nBas,nO,S,c,F) + if(level_shift > 0d0 .and. Conv > thresh) then + call level_shifting(level_shift, nBas, nOrb, nO, S, c, F) + endif ! Diagonalize Fock matrix - Fp = matmul(transpose(X),matmul(F,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eHF) - c = matmul(X,cp) + if(nBas .eq. nOrb) then + Fp = matmul(transpose(X), matmul(F, X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nOrb, cp, eHF) + c = matmul(X, cp) + else + Fp = matmul(transpose(c), matmul(F, c)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nOrb, cp, eHF) + c = matmul(c, cp) + endif ! Density matrix - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + !P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) + call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & + c(1,1), nBas, c(1,1), nBas, & + 0.d0, P(1,1), nBas) ! Dump results @@ -185,14 +217,16 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) + stop end if ! Compute dipole moments - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_RHF(nBas,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 @@ -205,4 +239,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN end if + deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) + end subroutine diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index b0bdb30..d200139 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -1,6 +1,9 @@ -subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO, & - X,ERHF,e,c,P) + +! --- + +subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas, nOrb, nC, nO, nV, nR, S, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, & + X, ERHF, e, c, P) ! Search for RHF solutions @@ -13,7 +16,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -26,11 +29,11 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -59,8 +62,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Output variables double precision,intent(out) :: ERHF - double precision,intent(out) :: e(nBas) - double precision,intent(inout):: c(nBas,nBas) + double precision,intent(out) :: e(nOrb) + double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: P(nBas,nBas) ! Memory allocation @@ -76,7 +79,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !-------------------! nS = (nO - nC)*(nV - nR) - allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas,nBas),ExpR(nBas,nBas)) + allocate(Aph(nS,nS), Bph(nS,nS), AB(nS,nS), Om(nS)) + allocate(R(nOrb,nOrb), ExpR(nOrb,nOrb)) !------------------! ! Search algorithm ! @@ -92,12 +96,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !---------------------! call wall_time(start_HF) - call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,e,c,P) + call RHF(.false., maxSCF, thresh, max_diis, guess, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, e, c, P) call wall_time(end_HF) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' write(*,*) !----------------------------------! @@ -108,14 +112,14 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,*) write(*,*) 'AO to MO transformation... Please be patient' write(*,*) - do ixyz=1,ncart - call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds' write(*,*) !-------------------------------------------------------------! @@ -124,12 +128,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ispin = 1 - call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph) - call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph) + call phLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) AB(:,:) = Aph(:,:) + Bph(:,:) - call diagonalize_matrix(nS,AB,Om) + call diagonalize_matrix(nS, AB, Om) Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' @@ -156,6 +160,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN if(eig < 0 .or. eig > nS) then write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...' write(*,*) + deallocate(Aph, Bph, AB, Om, R, ExpR) stop end if @@ -164,15 +169,15 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN R(:,:) = 0d0 ia = 0 do i=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nOrb-nR ia = ia + 1 R(a,i) = +AB(ia,eig) R(i,a) = -AB(ia,eig) end do end do - call matrix_exponential(nBas,R,ExpR) - c = matmul(c,ExpR) + call matrix_exponential(nOrb, R, ExpR) + c = matmul(c, ExpR) else @@ -191,4 +196,6 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !---------------! end do + deallocate(Aph, Bph, AB, Om, R, ExpR) + end subroutine diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 index d7b5385..154a4f9 100644 --- a/src/HF/RHF_stability.f90 +++ b/src/HF/RHF_stability.f90 @@ -30,7 +30,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) ! Memory allocation - allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS)) + allocate(A(nS,nS), B(nS,nS), AB(nS,nS), Om(nS)) !-------------------------------------------------------------! ! Stability analysis: Real RHF -> Real RHF @@ -148,5 +148,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) end if write(*,*)'-------------------------------------------------------------' write(*,*) + + deallocate(A, B, AB, Om) end subroutine diff --git a/src/HF/RMOM.f90 b/src/HF/RMOM.f90 index a49746d..fc43b76 100644 --- a/src/HF/RMOM.f90 +++ b/src/HF/RMOM.f90 @@ -190,6 +190,6 @@ subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) EK = 0.5d0*trace_matrix(nBas,matmul(P,K)) ERHF = ET + EV + EJ + EK - call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) + call print_RHF(nBas,nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) end subroutine diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index 7d86359..5a0ff02 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -1,5 +1,8 @@ -subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,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) ! Perform restricted open-shell Hartree-Fock calculation @@ -16,7 +19,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -28,14 +31,14 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_Sq integer :: n_diis double precision :: Conv double precision :: rcond @@ -62,8 +65,8 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Output variables double precision,intent(out) :: EROHF - double precision,intent(out) :: eHF(nBas) - double precision,intent(inout):: c(nBas,nBas) + double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: Ptot(nBas,nBas) ! Hello world @@ -76,19 +79,30 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Useful stuff - nBasSq = nBas*nBas + nBas_Sq = nBas*nBas ! Memory allocation - allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), & - P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), & - err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + 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)) + + allocate(Fp(nOrb,nOrb)) + allocate(cp(nOrb,nOrb)) + + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Guess coefficients and demsity matrices - call mo_guess(nBas,guess_type,S,Hc,X,c) - do ispin=1,nspin - P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin)))) + call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) + + 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 Ptot(:,:) = P(:,:,1) + P(:,:,2) @@ -120,51 +134,51 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! 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,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 @@ -176,7 +190,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,Ftot) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,Ftot) end if @@ -185,28 +199,29 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(level_shift,nBas,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(nBas,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 - P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin)))) + 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 Ptot(:,:) = P(:,:,1) + P(:,:,2) @@ -231,6 +246,8 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) + stop end if @@ -238,7 +255,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Compute final UHF energy call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_ROHF(nBas,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 @@ -248,4 +265,6 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN end if + deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) + end subroutine diff --git a/src/HF/ROHF_fock_matrix.f90 b/src/HF/ROHF_fock_matrix.f90 index bb23662..a88c57a 100644 --- a/src/HF/ROHF_fock_matrix.f90 +++ b/src/HF/ROHF_fock_matrix.f90 @@ -1,4 +1,7 @@ -subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) + +! --- + +subroutine ROHF_fock_matrix(nBas, nOrb, nOa, nOb, S, c, FaAO, FbAO, FAO) ! Construct the ROHF Fock matrix in the AO basis ! For open shells, the ROHF Fock matrix in the MO basis reads @@ -17,12 +20,12 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nOa integer,intent(in) :: nOb double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: c(nBas,nOrb) double precision,intent(inout):: FaAO(nBas,nBas) double precision,intent(inout):: FbAO(nBas,nBas) @@ -46,7 +49,7 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) ! Memory allocation - allocate(F(nBas,nBas),Fa(nBas,nBas),Fb(nBas,nBas)) + allocate(F(nOrb,nOrb), Fa(nOrb,nOrb), Fb(nOrb,nOrb)) ! Roothan canonicalization parameters @@ -61,14 +64,14 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) ! Number of closed, open, and virtual orbitals - nC = min(nOa,nOb) + nC = min(nOa, nOb) nO = abs(nOa - nOb) nV = nBas - nC - nO ! Block-by-block Fock matrix - call AOtoMO(nBas,c,FaAO,Fa) - call AOtoMO(nBas,c,FbAO,Fb) + call AOtoMO(nBas, nOrb, c, FaAO, Fa) + call AOtoMO(nBas, nOrb, c, FbAO, Fb) F(1:nC, 1:nC ) = aC*Fa(1:nC, 1:nC ) + bC*Fb(1:nC, 1:nC ) F(1:nC, nC+1:nC+nO ) = Fb(1:nC, nC+1:nC+nO ) @@ -82,8 +85,10 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) F(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) = Fa(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) F(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) = aV*Fa(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) + bV*Fb(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) - call MOtoAO(nBas,S,c,F,FAO) - call MOtoAO(nBas,S,c,Fa,FaAO) - call MOtoAO(nBas,S,c,Fb,FbAO) + call MOtoAO(nBas, nOrb, S, c, F, FAO) + call MOtoAO(nBas, nOrb, S, c, Fa, FaAO) + call MOtoAO(nBas, nOrb, S, c, Fb, FbAO) + + deallocate(F, Fa, Fb) end subroutine diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 919cb0f..faa7de1 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -85,7 +85,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Guess coefficients and demsity matrices do ispin=1,nspin - call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin)) + call mo_guess(nBas,nBas,guess_type,S,Hc,X,c(:,:,ispin)) P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) end do @@ -186,7 +186,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) + call level_shifting(level_shift,nBas,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) end do end if diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 index 30ae0e2..dcc1025 100644 --- a/src/HF/UHF_search.f90 +++ b/src/HF/UHF_search.f90 @@ -124,8 +124,8 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Transform dipole-related integrals do ixyz=1,ncart - call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block diff --git a/src/HF/cRHF.f90 b/src/HF/cRHF.f90 index 9640a48..36cfa41 100644 --- a/src/HF/cRHF.f90 +++ b/src/HF/cRHF.f90 @@ -94,7 +94,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Guess coefficients and density matrix - call mo_guess(nBas,guess_type,S,Hc,X,c) + call mo_guess(nBas,nBas,guess_type,S,Hc,X,c) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) ! Initialization @@ -166,7 +166,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Level shift - if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F) + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,c,F) ! Diagonalize Fock matrix @@ -207,7 +207,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Compute dipole moments call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_RHF(nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) + call print_RHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) ! Testing zone diff --git a/src/HF/core_guess.f90 b/src/HF/core_guess.f90 index 7e45e8e..8d34444 100644 --- a/src/HF/core_guess.f90 +++ b/src/HF/core_guess.f90 @@ -1,4 +1,4 @@ -subroutine core_guess(nBas,Hc,X,c) +subroutine core_guess(nBas, nOrb, Hc, X, c) ! Core guess of the molecular orbitals for HF calculation @@ -6,9 +6,9 @@ subroutine core_guess(nBas,Hc,X,c) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) ! Local variables @@ -18,16 +18,19 @@ subroutine core_guess(nBas,Hc,X,c) ! Output variables - double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: c(nBas,nOrb) ! Memory allocation - allocate(cp(nBas,nBas),e(nBas)) + allocate(cp(nOrb,nOrb), e(nOrb)) ! Core guess - cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - call diagonalize_matrix(nBas,cp,e) - c(:,:) = matmul(X(:,:),cp(:,:)) + cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:))) + + call diagonalize_matrix(nOrb, cp, e) + c(:,:) = matmul(X(:,:), cp(:,:)) + + deallocate(cp, e) end subroutine diff --git a/src/HF/huckel_guess.f90 b/src/HF/huckel_guess.f90 index 658c853..7afd0a0 100644 --- a/src/HF/huckel_guess.f90 +++ b/src/HF/huckel_guess.f90 @@ -1,4 +1,4 @@ -subroutine huckel_guess(nBas,S,Hc,X,c) +subroutine huckel_guess(nBas, nOrb, S, Hc, X, c) ! Hickel guess @@ -6,10 +6,10 @@ subroutine huckel_guess(nBas,S,Hc,X,c) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) ! Local variables @@ -20,7 +20,7 @@ subroutine huckel_guess(nBas,S,Hc,X,c) ! Output variables - double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: c(nBas,nOrb) ! Memory allocation @@ -32,9 +32,9 @@ subroutine huckel_guess(nBas,S,Hc,X,c) ! GWH approximation - do mu=1,nBas + do mu = 1, nBas F(mu,mu) = Hc(mu,mu) - do nu=mu+1,nBas + do nu = mu+1, nBas F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu)) F(nu,mu) = F(mu,nu) @@ -42,6 +42,8 @@ subroutine huckel_guess(nBas,S,Hc,X,c) end do end do - call core_guess(nBas,F,X,c) + call core_guess(nBas, nOrb, F, X, c) + + deallocate(F) end subroutine diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index fef8d87..2046569 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -1,4 +1,7 @@ -subroutine mo_guess(nBas,guess_type,S,Hc,X,c) + +! --- + +subroutine mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) ! Guess of the molecular orbitals for HF calculation @@ -6,15 +9,15 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: guess_type double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) ! Output variables - double precision,intent(inout) :: c(nBas,nBas) + double precision,intent(inout) :: c(nBas,nOrb) if(guess_type == 0) then @@ -24,12 +27,12 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c) elseif(guess_type == 1) then write(*,*) 'Core guess...' - call core_guess(nBas,Hc,X,c) + call core_guess(nBas, nOrb, Hc, X, c) elseif(guess_type == 2) then write(*,*) 'Huckel guess...' - call huckel_guess(nBas,S,Hc,X,c) + call huckel_guess(nBas, nOrb, S, Hc, X, c) elseif(guess_type == 3) then diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 3a06d31..277db55 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -1,4 +1,7 @@ -subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) + +! --- + +subroutine print_RHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -7,10 +10,10 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: ENuc double precision,intent(in) :: ET double precision,intent(in) :: EV @@ -75,13 +78,13 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital coefficients ' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,cHF) + call matout(nBas, nOrb, cHF) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eHF) + call vecout(nOrb, eHF) write(*,*) end subroutine diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90 index 9d29535..0297bb4 100644 --- a/src/HF/print_ROHF.f90 +++ b/src/HF/print_ROHF.f90 @@ -1,14 +1,17 @@ -subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) + +! --- + +subroutine print_ROHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROHF, dipole) ! Print one- and two-electron energies and other stuff for RoHF calculation implicit none include 'parameters.h' - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO(nspin) - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: c(nBas,nOrb) double precision,intent(in) :: ENuc double precision,intent(in) :: ET(nspin) double precision,intent(in) :: EV(nspin) @@ -31,7 +34,7 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) do ispin=1,nspin if(nO(ispin) > 0) then HOMO(ispin) = eHF(nO(ispin)) - if(nO(ispin) < nBas) then + if(nO(ispin) < nOrb) then LUMO(ispin) = eHF(nO(ispin)+1) else LUMO(ispin) = 0d0 @@ -102,13 +105,13 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') 'ROHF orbital coefficients ' write(*,'(A50)') '-----------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' ROHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eHF) + call vecout(nOrb, eHF) write(*,*) end subroutine diff --git a/src/MP/RMP.f90 b/src/MP/RMP.f90 index 9fba986..daf3f1e 100644 --- a/src/MP/RMP.f90 +++ b/src/MP/RMP.f90 @@ -41,7 +41,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_MP) t_MP = end_MP - start_MP - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP2 = ',t_MP,' seconds' write(*,*) end if @@ -57,7 +57,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_MP) t_MP = end_MP - start_MP - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP2 = ',t_MP,' seconds' write(*,*) end if diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 6af981b..c6a77c9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -15,7 +15,7 @@ program QuAcK logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh - integer :: nNuc,nBas + integer :: nNuc, nBas, nOrb integer :: nC(nspin) integer :: nO(nspin) integer :: nV(nspin) @@ -31,6 +31,7 @@ program QuAcK double precision,allocatable :: X(:,:) double precision,allocatable :: dipole_int_AO(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) + double precision,allocatable :: Uvec(:,:), Uval(:) double precision :: start_QuAcK,end_QuAcK,t_QuAcK double precision :: start_int ,end_int ,t_int @@ -68,6 +69,10 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest + integer :: i, j, j0 + double precision :: acc_d, acc_nd + double precision, allocatable :: tmp1(:,:), tmp2(:,:) + !-------------! ! Hello World ! !-------------! @@ -120,15 +125,16 @@ program QuAcK doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA) -!------------------------------------------------! -! Read input information ! -!------------------------------------------------! -! nC = number of core orbitals ! -! nO = number of occupied orbitals ! -! nV = number of virtual orbitals (see below) ! -! nR = number of Rydberg orbitals ! -! nBas = number of basis functions (see below) ! -!------------------------------------------------! +!---------------------------------------------------! +! Read input information ! +!---------------------------------------------------! +! nC = number of core orbitals ! +! nO = number of occupied orbitals ! +! nV = number of virtual orbitals (see below) ! +! nR = number of Rydberg orbitals ! +! nBas = number of basis functions in AOs ! +! nOrb = number of basis functions in MOs ! +!---------------------------------------------------! call read_molecule(nNuc,nO,nC,nR) allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) @@ -141,7 +147,7 @@ program QuAcK ! Read basis set information from PySCF ! !---------------------------------------! - call read_basis_pyscf(nBas,nO,nV) + call read_basis_pyscf(nBas, nO, nV) !--------------------------------------! ! Read one- and two-electron integrals ! @@ -149,26 +155,84 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), & - ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart)) + allocate(S(nBas,nBas)) + allocate(T(nBas,nBas)) + allocate(V(nBas,nBas)) + allocate(Hc(nBas,nBas)) + allocate(ERI_AO(nBas,nBas,nBas,nBas)) + allocate(dipole_int_AO(nBas,nBas,ncart)) ! Read integrals call wall_time(start_int) - call read_integrals(nBas,S,T,V,Hc,ERI_AO) - call read_dipole_integrals(nBas,dipole_int_AO) + call read_integrals(nBas, S(1,1), T(1,1), V(1,1), Hc(1,1), ERI_AO(1,1,1,1)) + call read_dipole_integrals(nBas, dipole_int_AO) call wall_time(end_int) t_int = end_int - start_int write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds' write(*,*) ! Compute orthogonalization matrix - call orthogonalization_matrix(nBas,S,X) + !call orthogonalization_matrix(nBas, S, X) + + allocate(Uvec(nBas,nBas), Uval(nBas)) + + Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas) + call diagonalize_matrix(nBas, Uvec, Uval) + + nOrb = 0 + do i = 1, nBas + if(Uval(i) > 1d-6) then + Uval(i) = 1d0 / dsqrt(Uval(i)) + nOrb = nOrb + 1 + else + write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + end if + end do + + write(*,'(A38)') '--------------------------------------' + write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas + write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nOrb + write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nOrb)/dble(nBas)) + write(*,'(A38)') '--------------------------------------' + write(*,*) + + j0 = nBas - nOrb + allocate(X(nBas,nOrb)) + do j = j0+1, nBas + do i = 1, nBas + X(i,j-j0) = Uvec(i,j) * Uval(j) + enddo + enddo + + deallocate(Uvec, Uval) + + !! check if X.T S X = 1_(nOrb,nOrb) + !allocate(tmp1(nOrb,nBas), tmp2(nOrb,nOrb)) + !call dgemm("T", "N", nOrb, nBas, nBas, 1.d0, & + ! X(1,1), nBas, S(1,1), nBas, & + ! 0.d0, tmp1(1,1), nOrb) + !call dgemm("N", "N", nOrb, nOrb, nBas, 1.d0, & + ! tmp1(1,1), nOrb, X(1,1), nBas, & + ! 0.d0, tmp2(1,1), nOrb) + !acc_d = 0.d0 + !acc_nd = 0.d0 + !do i = 1, nOrb + ! !write(*,'(1000(F15.7,2X))') (tmp2(i,j), j = 1, nOrb) + ! acc_d = acc_d + tmp2(i,i) + ! do j = 1, nOrb + ! if(j == i) cycle + ! acc_nd = acc_nd + dabs(tmp2(j,i)) + ! enddo + !enddo + !print*, ' diag part: ', dabs(acc_d - dble(nOrb)) / dble(nOrb) + !print*, ' non-diag part: ', acc_nd + !deallocate(tmp1, tmp2) !---------------------! ! Choose QuAcK branch ! @@ -200,7 +264,7 @@ program QuAcK dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -216,7 +280,7 @@ program QuAcK dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -228,13 +292,13 @@ program QuAcK !--------------------------! if(doGQuAcK) & - call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & - maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & - maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & - maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & + maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & + maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & + maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) !-----------! @@ -256,7 +320,7 @@ program QuAcK call wall_time(end_QuAcK) t_QuAcK = end_QuAcK - start_QuAcK - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for QuAcK = ',t_QuAcK,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds' write(*,*) end program diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 9632c32..5725c4d 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -2,7 +2,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -29,7 +29,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh - integer,intent(in) :: nNuc,nBas + integer,intent(in) :: nNuc,nBas,nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -42,7 +42,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) @@ -109,8 +109,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Memory allocation ! !-------------------! - allocate(cHF(nBas,nBas),eHF(nBas),PHF(nBas,nBas), & - dipole_int_MO(nBas,nBas,ncart),ERI_MO(nBas,nBas,nBas,nBas)) + allocate(cHF(nBas,nOrb)) + allocate(eHF(nOrb)) + allocate(PHF(nBas,nBas)) + allocate(dipole_int_MO(nOrb,nOrb,ncart)) + allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) !---------------------! ! Hartree-Fock module ! @@ -119,12 +122,12 @@ 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,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) call wall_time(end_HF) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' write(*,*) end if @@ -132,12 +135,12 @@ 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,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) call wall_time(end_HF) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ROHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ROHF = ',t_HF,' seconds' write(*,*) end if @@ -154,18 +157,18 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Read and transform dipole-related integrals - do ixyz=1,ncart - call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas, nOrb, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do ! 4-index transform - call AOtoMO_ERI_RHF(nBas,cHF,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, cHF, ERI_AO, ERI_MO) call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds' write(*,*) !-----------------------------------! @@ -177,11 +180,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dostab) then call wall_time(start_stab) - call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO) + call RHF_stability(nOrb, nC, nO, nV, nR, nS, eHF, ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds' write(*,*) end if @@ -189,12 +192,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dosearch) then call wall_time(start_stab) - call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X,ERHF,eHF,cHF,PHF) + call RHF_search(maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas, nOrb, nC, nO, nV, nR, S, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, & + dipole_int_MO, X, ERHF, eHF, cHF, PHF) call wall_time(end_stab) t_stab = end_stab - start_stab - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds' write(*,*) end if @@ -208,11 +212,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doMP) then call wall_time(start_MP) - call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call RMP(dotest, doMP2, doMP3, reg_MP, nOrb, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, ERHF, eHF) call wall_time(end_MP) t_MP = end_MP - start_MP - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP = ',t_MP,' seconds' write(*,*) end if @@ -227,12 +231,12 @@ 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,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) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CC = ',t_CC,' seconds' write(*,*) end if @@ -246,12 +250,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doCI) then call wall_time(start_CI) - call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & - eHF,ERHF,cHF,S) + call RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, & + nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, ERHF) call wall_time(end_CI) t_CI = end_CI - start_CI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CI = ',t_CI,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CI = ',t_CI,' seconds' write(*,*) end if @@ -265,12 +269,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doRPA) then call wall_time(start_RPA) - call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S) + call RRPA(dotest, dophRPA, dophRPAx, docrRPA, doppRPA, TDA, doACFDT, exchange_kernel, singlet, triplet, & + nOrb, nC, nO, nV, nR, nS, ENuc, ERHF, ERI_MO, dipole_int_MO, eHF) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds' write(*,*) end if @@ -284,14 +288,14 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGF) then call wall_time(start_GF) - call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & - dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, & - nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & - dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + call RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm_GF, maxSCF_GF, & + thresh_GF, max_diis_GF, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, lin_GF, & + eta_GF, reg_GF, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, & + S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -305,14 +309,14 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGW) then call wall_time(start_GW) - call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, & - doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & - ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + call RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF_GW, thresh_GW, max_diis_GW, & + doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, & + lin_GW, eta_GW, reg_GW, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, & + V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GW = ',t_GW,' seconds' write(*,*) end if @@ -326,14 +330,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGT) then call wall_time(start_GT) - call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & - TDA_T,TDA,dBSE,dTDA,singlet,triplet,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + call RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & + maxSCF_GT, thresh_GT, max_diis_GT, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, & + TDA_T, TDA, dBSE, dTDA, singlet, triplet, lin_GT, eta_GT, reg_GT, nNuc, ZNuc, rNuc, ENuc, & + nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, & + dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GT = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GT = ',t_GT,' seconds' write(*,*) end if diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index e3453b0..6e45b0d 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -163,8 +163,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do ! Read and transform dipole-related integrals do ixyz=1,ncart - call AOtoMO(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block diff --git a/src/RPA/RRPA.f90 b/src/RPA/RRPA.f90 index 93bf457..ab38932 100644 --- a/src/RPA/RRPA.f90 +++ b/src/RPA/RRPA.f90 @@ -1,5 +1,5 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,cHF,S) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Random-phase approximation module @@ -29,8 +29,6 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -49,7 +47,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker call wall_time(end_RPA) t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds' write(*,*) end if @@ -65,7 +63,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker call wall_time(end_RPA) t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPAx = ',t_RPA,' seconds' write(*,*) end if @@ -81,7 +79,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker call wall_time(end_RPA) t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if @@ -97,7 +95,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker call wall_time(end_RPA) t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 index 8435b32..8b627fa 100644 --- a/src/utils/level_shifting.f90 +++ b/src/utils/level_shifting.f90 @@ -1,4 +1,4 @@ -subroutine level_shifting(level_shift,nBas,nO,S,c,F) +subroutine level_shifting(level_shift, nBas, nOrb, nO, S, c, F) ! Perform level-shifting on the Fock matrix @@ -7,10 +7,10 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F) ! Input variables double precision,intent(in) :: level_shift - integer,intent(in) :: nBas + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: c(nBas,nOrb) ! Local variables @@ -23,15 +23,17 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F) double precision,intent(inout):: F(nBas,nBas) - allocate(F_MO(nBas,nBas),Sc(nBas,nBas)) + allocate(F_MO(nOrb,nOrb), Sc(nBas,nOrb)) - F_MO(:,:) = matmul(transpose(c),matmul(F,c)) + F_MO(:,:) = matmul(transpose(c), matmul(F, c)) - do a=nO+1,nBas + do a = nO+1, nOrb F_MO(a,a) = F_MO(a,a) + level_shift end do - Sc(:,:) = matmul(S,c) - F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc))) + Sc(:,:) = matmul(S, c) + F(:,:) = matmul(Sc, matmul(F_MO, transpose(Sc))) + + deallocate(F_MO, Sc) end subroutine diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index d7e0089..af781e1 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -1,4 +1,7 @@ -subroutine orthogonalization_matrix(nBas,S,X) + +! --- + +subroutine orthogonalization_matrix(nBas, S, X) ! Compute the orthogonalization matrix X @@ -35,14 +38,32 @@ subroutine orthogonalization_matrix(nBas,S,X) if(ortho_type == 1) then + ! + ! S V = V s where + ! + ! V.T V = 1 and s > 0 (S is positive def) + ! + ! S = V s V.T + ! = V s^0.5 s^0.5 V.T + ! = V s^0.5 V.T V s^0.5 V.T + ! = S^0.5 S^0.5 + ! + ! where + ! + ! S^0.5 = V s^0.5 V.T + ! + ! X = S^(-0.5) + ! = V s^(-0.5) V.T + ! + ! write(*,*) ! write(*,*) ' Lowdin orthogonalization' ! write(*,*) Uvec = S - call diagonalize_matrix(nBas,Uvec,Uval) + call diagonalize_matrix(nBas, Uvec, Uval) - do i=1,nBas + do i = 1, nBas if(Uval(i) < thresh) then @@ -50,7 +71,7 @@ subroutine orthogonalization_matrix(nBas,S,X) end if - Uval(i) = 1d0/sqrt(Uval(i)) + Uval(i) = 1d0 / dsqrt(Uval(i)) end do @@ -63,13 +84,13 @@ subroutine orthogonalization_matrix(nBas,S,X) ! write(*,*) Uvec = S - call diagonalize_matrix(nBas,Uvec,Uval) + call diagonalize_matrix(nBas, Uvec, Uval) - do i=1,nBas + do i = 1, nBas if(Uval(i) > thresh) then - Uval(i) = 1d0/sqrt(Uval(i)) + Uval(i) = 1d0 / dsqrt(Uval(i)) else @@ -79,7 +100,7 @@ subroutine orthogonalization_matrix(nBas,S,X) end do - call AD(nBas,Uvec,Uval) + call AD(nBas, Uvec, Uval) X = Uvec elseif(ortho_type == 3) then @@ -117,3 +138,6 @@ subroutine orthogonalization_matrix(nBas,S,X) end if end subroutine + +! --- + diff --git a/src/utils/read_basis_pyscf.f90 b/src/utils/read_basis_pyscf.f90 index a677323..42dfde2 100644 --- a/src/utils/read_basis_pyscf.f90 +++ b/src/utils/read_basis_pyscf.f90 @@ -1,4 +1,4 @@ -subroutine read_basis_pyscf(nBas,nO,nV) +subroutine read_basis_pyscf(nBas_AOs, nO, nV) ! Read basis set information from PySCF @@ -14,23 +14,23 @@ subroutine read_basis_pyscf(nBas,nO,nV) ! Output variables integer,intent(out) :: nV(nspin) - integer,intent(out) :: nBas + integer,intent(out) :: nBas_AOs !------------------------------------------------------------------------ ! Primary basis set information !------------------------------------------------------------------------ open(unit=3,file='int/nBas.dat') - read(3,*) nBas + read(3, *) nBas_AOs close(unit=3) - write(*,'(A28)') '------------------' - write(*,'(A28,1X,I16)') 'Number of basis functions',nBas - write(*,'(A28)') '------------------' - write(*,*) +! write(*,'(A38)') '--------------------------------------' +! write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs +! write(*,'(A38)') '--------------------------------------' +! write(*,*) ! Number of virtual orbitals - nV(:) = nBas - nO(:) + nV(:) = nBas_AOs - nO(:) end subroutine diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index 9e215e4..ef8e1d9 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -1,4 +1,4 @@ -subroutine read_integrals(nBas,S,T,V,Hc,G) +subroutine read_integrals(nBas_AOs, S, T, V, Hc, G) ! Read one- and two-electron integrals from files @@ -7,7 +7,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs ! Local variables @@ -18,11 +18,11 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) ! Output variables - double precision,intent(out) :: S(nBas,nBas) - double precision,intent(out) :: T(nBas,nBas) - double precision,intent(out) :: V(nBas,nBas) - double precision,intent(out) :: Hc(nBas,nBas) - double precision,intent(out) :: G(nBas,nBas,nBas,nBas) + double precision,intent(out) :: S(nBas_AOs,nBas_AOs) + double precision,intent(out) :: T(nBas_AOs,nBas_AOs) + double precision,intent(out) :: V(nBas_AOs,nBas_AOs) + double precision,intent(out) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) ! Open file with integrals @@ -35,7 +35,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) open(unit=8 ,file='int/Ov.dat') open(unit=9 ,file='int/Kin.dat') open(unit=10,file='int/Nuc.dat') - open(unit=11,file='int/ERI.dat') open(unit=21,file='int/x.dat') open(unit=22,file='int/y.dat') @@ -75,31 +74,29 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) Hc(:,:) = T(:,:) + V(:,:) -! Read nuclear integrals +! Read 2e-integrals - G(:,:,:,:) = 0d0 - do - read(11,*,end=11) mu,nu,la,si,ERI +! ! formatted file +! open(unit=11, file='int/ERI.dat') +! G(:,:,:,:) = 0d0 +! do +! read(11,*,end=11) mu, nu, la, si, ERI +! ERI = lambda*ERI +! G(mu,nu,la,si) = ERI ! <12|34> +! G(la,nu,mu,si) = ERI ! <32|14> +! G(mu,si,la,nu) = ERI ! <14|32> +! G(la,si,mu,nu) = ERI ! <34|12> +! G(si,mu,nu,la) = ERI ! <41|23> +! G(nu,la,si,mu) = ERI ! <23|41> +! G(nu,mu,si,la) = ERI ! <21|43> +! G(si,la,nu,mu) = ERI ! <43|21> +! end do +! 11 close(unit=11) - ERI = lambda*ERI -! <12|34> - G(mu,nu,la,si) = ERI -! <32|14> - G(la,nu,mu,si) = ERI -! <14|32> - G(mu,si,la,nu) = ERI -! <34|12> - G(la,si,mu,nu) = ERI -! <41|23> - G(si,mu,nu,la) = ERI -! <23|41> - G(nu,la,si,mu) = ERI -! <21|43> - G(nu,mu,si,la) = ERI -! <43|21> - G(si,la,nu,mu) = ERI - end do - 11 close(unit=11) + ! binary file + open(unit=11, file='int/ERI.bin', form='unformatted', access='stream') + read(11) G + close(11) ! Print results @@ -107,24 +104,24 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Overlap integrals' write(*,'(A28)') '----------------------' - call matout(nBas,nBas,S) + call matout(nBas_AOs,nBas_AOs,S) write(*,*) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Kinetic integrals' write(*,'(A28)') '----------------------' - call matout(nBas,nBas,T) + call matout(nBas_AOs,nBas_AOs,T) write(*,*) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Nuclear integrals' write(*,'(A28)') '----------------------' - call matout(nBas,nBas,V) + call matout(nBas_AOs,nBas_AOs,V) write(*,*) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Electron repulsion integrals' write(*,'(A28)') '----------------------' - do la=1,nBas - do si=1,nBas - call matout(nBas,nBas,G(1,1,la,si)) + do la=1,nBas_AOs + do si=1,nBas_AOs + call matout(nBas_AOs, nBas_AOs, G(1,1,la,si)) end do end do write(*,*) diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index 80b8322..316a87c 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -65,7 +65,7 @@ subroutine diagonal_matrix(N,D,A) end subroutine !------------------------------------------------------------------------ -subroutine matrix_exponential(N,A,ExpA) +subroutine matrix_exponential(N, A, ExpA) ! Compute Exp(A) @@ -81,7 +81,7 @@ subroutine matrix_exponential(N,A,ExpA) ! Memory allocation - allocate(W(N,N),tau(N),t(N,N)) + allocate(W(N,N), tau(N), t(N,N)) ! Initialize @@ -89,8 +89,8 @@ subroutine matrix_exponential(N,A,ExpA) ! Diagonalize - W(:,:) = - matmul(A,A) - call diagonalize_matrix(N,W,tau) + W(:,:) = - matmul(A, A) + call diagonalize_matrix(N, W, tau) ! do i=1,N ! tau(i) = max(abs(tau(i)),1d-14) @@ -99,16 +99,18 @@ subroutine matrix_exponential(N,A,ExpA) ! Construct cos part - call diagonal_matrix(N,cos(tau),t) - t(:,:) = matmul(t,transpose(W)) - ExpA(:,:) = ExpA(:,:) + matmul(W,t) + call diagonal_matrix(N, cos(tau), t) + t(:,:) = matmul(t, transpose(W)) + ExpA(:,:) = ExpA(:,:) + matmul(W, t) ! Construct sin part - call diagonal_matrix(N,sin(tau)/tau,t) - t(:,:) = matmul(t,transpose(W)) - t(:,:) = matmul(t,A) - ExpA(:,:) = ExpA(:,:) + matmul(W,t) + call diagonal_matrix(N, sin(tau)/tau, t) + t(:,:) = matmul(t, transpose(W)) + t(:,:) = matmul(t, A) + ExpA(:,:) = ExpA(:,:) + matmul(W, t) + + deallocate(W, tau, t) end subroutine