From b7af468f114980df42d76d7a15f60d47c2454f00 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 7 Dec 2024 02:20:05 +0100 Subject: [PATCH] Working on s8 rep of 2e-integrals --- PyDuck.py | 67 ++-- input/hardware | 2 - input/hpc_flags | 4 + src/AOtoMO/Hartree_matrix_AO_basis.f90 | 103 +++++ src/HF/RHF_hpc.f90 | 247 ++++++++++++ src/QuAcK/GQuAcK.f90 | 34 +- src/QuAcK/QuAcK.f90 | 60 ++- src/QuAcK/RQuAcK.f90 | 41 +- src/QuAcK/RQuAcK_hpc.f90 | 353 ++++++++++++++++++ src/QuAcK/UQuAcK.f90 | 18 +- src/QuAcK/read_hardware.f90 | 45 --- src/QuAcK/read_hpc_flags.f90 | 39 ++ ...ad_integrals.f90 => read_1e_integrals.f90} | 68 +--- src/utils/read_2e_integrals.f90 | 110 ++++++ 14 files changed, 1019 insertions(+), 172 deletions(-) delete mode 100644 input/hardware create mode 100644 input/hpc_flags create mode 100644 src/HF/RHF_hpc.f90 create mode 100644 src/QuAcK/RQuAcK_hpc.f90 delete mode 100644 src/QuAcK/read_hardware.f90 create mode 100644 src/QuAcK/read_hpc_flags.f90 rename src/utils/{read_integrals.f90 => read_1e_integrals.f90} (59%) create mode 100644 src/utils/read_2e_integrals.f90 diff --git a/PyDuck.py b/PyDuck.py index 6318204..45113de 100644 --- a/PyDuck.py +++ b/PyDuck.py @@ -24,8 +24,10 @@ 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('--print_2e', default=True, action='store_true', help='If True, print 2e-integrals to disk.') +parser.add_argument('--formatted_2e', default=False, action='store_true', help='Add this option if you want to print formatted 2e-integrals.') parser.add_argument('--mmap_2e', default=False, action='store_true', help='If True, avoid using DRAM when generating 2e-integrals.') +parser.add_argument('--aosym_2e', default=False, action='store_true', help='If True, use 8-fold symmetry 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.') @@ -41,7 +43,9 @@ multiplicity=args.multiplicity xyz=args.xyz + '.xyz' cartesian=args.cartesian print_2e=args.print_2e +formatted_2e=args.formatted_2e mmap_2e=args.mmap_2e +aosym_2e=args.aosym_2e working_dir=args.working_dir #Read molecule @@ -63,6 +67,7 @@ mol = gto.M( basis = input_basis, charge = charge, spin = multiplicity - 1 +# symmetry = True # Enable symmetry ) #Fix the unit for the lengths @@ -144,34 +149,46 @@ def write_tensor_to_file(tensor,size,file_name,cutoff=1e-15): f.write('\n') f.close() -# Write two-electron integrals to HD -ti_2e = time.time() if print_2e: - # (formatted) - output_file_path = working_dir + '/int/ERI.dat' - subprocess.call(['rm', '-f', output_file_path]) - eri_ao = mol.intor('int2e') - write_tensor_to_file(eri_ao, norb, output_file_path) -else: - # (binary) - output_file_path = working_dir + '/int/ERI.bin' - subprocess.call(['rm', '-f', output_file_path]) - if(mmap_2e): - # avoid using DRAM - eri_shape = (norb, norb, norb, norb) - eri_mmap = np.memmap(output_file_path, dtype='float64', mode='w+', shape=eri_shape) - mol.intor('int2e', out=eri_mmap) - for i in range(norb): - eri_mmap[i, :, :, :] = eri_mmap[i, :, :, :].transpose(1, 0, 2) - eri_mmap.flush() - del eri_mmap - else: - eri_ao = mol.intor('int2e').transpose(0, 2, 1, 3) # chem -> phys + # Write two-electron integrals to HD + ti_2e = time.time() + + if formatted_2e: + output_file_path = working_dir + '/int/ERI.dat' + subprocess.call(['rm', '-f', output_file_path]) + eri_ao = mol.intor('int2e') + write_tensor_to_file(eri_ao, norb, output_file_path) + + if aosym_2e: + output_file_path = working_dir + '/int/ERI_chem.bin' + subprocess.call(['rm', '-f', output_file_path]) + eri_ao = mol.intor('int2e', aosym='s8') + print(eri_ao.shape) f = open(output_file_path, 'w') eri_ao.tofile(output_file_path) f.close() -te_2e = time.time() -print("Wall time for writing 2e-integrals (physicist notation) to disk: {:.3f} seconds".format(te_2e - ti_2e)) + else: + output_file_path = working_dir + '/int/ERI.bin' + subprocess.call(['rm', '-f', output_file_path]) + if(mmap_2e): + # avoid using DRAM + eri_shape = (norb, norb, norb, norb) + eri_mmap = np.memmap(output_file_path, dtype='float64', mode='w+', shape=eri_shape) + mol.intor('int2e', out=eri_mmap) + for i in range(norb): + eri_mmap[i, :, :, :] = eri_mmap[i, :, :, :].transpose(1, 0, 2) + eri_mmap.flush() + del eri_mmap + else: + eri_ao = mol.intor('int2e').transpose(0, 2, 1, 3) # chem -> phys + f = open(output_file_path, 'w') + eri_ao.tofile(output_file_path) + f.close() + + te_2e = time.time() + print("Wall time for writing 2e-integrals to disk: {:.3f} seconds".format(te_2e - ti_2e)) + sys.stdout.flush() + diff --git a/input/hardware b/input/hardware deleted file mode 100644 index 99c395a..0000000 --- a/input/hardware +++ /dev/null @@ -1,2 +0,0 @@ -# if True (T), use GPU - F diff --git a/input/hpc_flags b/input/hpc_flags new file mode 100644 index 0000000..120f367 --- /dev/null +++ b/input/hpc_flags @@ -0,0 +1,4 @@ +# if True (T), switch to HPC mode + F +# if True (T), use GPU + F diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 03000a0..ed819c6 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -32,3 +32,106 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H) end do end subroutine + +! --- + +subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H) + + implicit none + + integer, intent(in) :: nBas + integer*8, intent(in) :: ERI_size + double precision, intent(in) :: P(nBas,nBas) + double precision, intent(in) :: ERI_chem(ERI_size) + double precision, intent(out) :: H(nBas,nBas) + + integer :: mu, nu, la, si + integer*8 :: munu0, munu + integer*8 :: sila0, sila + integer*8 :: munulasi0, munulasi + + integer*8, external :: Yoshimine_ind + + do nu = 1, nBas + do mu = 1, nBas + H(mu,nu) = 0.d0 + do si = 1, nBas + do la = 1, nBas + munulasi = Yoshimine_ind(mu, nu, la, si) + H(mu,nu) = H(mu,nu) + P(la,si) * ERI_chem(munulasi) + enddo + enddo + enddo + enddo + + +! do nu = 1, nBas +! munu0 = (nu * (nu + 1)) / 2 +! +! do mu = 1, nu +! munu = munu0 + mu +! munulasi0 = (munu * (munu + 1)) / 2 +! +! H(mu,nu) = 0.d0 +! +! do si = 1, nu +! sila0 = (si * (si + 1)) / 2 +! +! do la = 1, si +! sila = sila0 + la +! +! if(nu == si .and. mu < la) cycle +! +! munulasi = munulasi0 + sila +! +! H(mu,nu) = H(mu,nu) + 4.d0 * P(la,si) * ERI_chem(munulasi) +! enddo +! enddo +! enddo +! enddo +! +! +! do nu = 1, nBas +! do mu = nu+1, nBas +! H(mu,nu) = H(nu,mu) +! enddo +! enddo + + return +end subroutine + +! --- + +integer*8 function Yoshimine_ind(a, b, c, d) + + implicit none + + integer, intent(in) :: a, b, c, d + + integer*8 :: ab, cd, abcd + + if(a > b) then + ab = (a * (a - 1)) / 2 + b + else + ab = (b * (b - 1)) / 2 + a + endif + + if(c > d) then + cd = (c * (c - 1)) / 2 + d + else + cd = (d * (d - 1)) / 2 + c + endif + + if(ab > cd) then + abcd = (ab * (ab - 1)) / 2 + cd + else + abcd = (cd * (cd - 1)) / 2 + ab + endif + + Yoshimine_ind = abcd + + return +end + +! --- + diff --git a/src/HF/RHF_hpc.f90 b/src/HF/RHF_hpc.f90 new file mode 100644 index 0000000..ef6ea47 --- /dev/null +++ b/src/HF/RHF_hpc.f90 @@ -0,0 +1,247 @@ +subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,dipole_int,X,ERHF,eHF,c,P,F) + +! Perform restricted Hartree-Fock calculation + + implicit none + include 'parameters.h' + +! Input variables + + character(len=256),intent(in) :: working_dir + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + integer,intent(in) :: guess_type + double precision,intent(in) :: thresh + double precision,intent(in) :: level_shift + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nO + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + 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,nOrb) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ii, jj + integer :: nSCF + integer :: nBas_Sq + integer :: n_diis + integer*8 :: ERI_size + double precision :: diff, diff_loc + double precision :: ET + double precision :: EV + double precision :: EJ + double precision :: EK + double precision :: dipole(ncart) + double precision :: Conv + double precision :: rcond + double precision,external :: trace_matrix + double precision,allocatable :: err(:,:) + double precision,allocatable :: err_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: Fp(:,:) + double precision,allocatable :: ERI_chem(:) + double precision,allocatable :: ERI_phys(:,:,:,:), J_deb(:,:) + + +! Output variables + + double precision,intent(out) :: ERHF + double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: c(nBas,nOrb) + double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: F(nBas,nBas) + +! Hello world + + write(*,*) + write(*,*)'****************************************' + write(*,*)'* Restricted HF Calculation (HPC mode) *' + write(*,*)'****************************************' + write(*,*) + +! Useful quantities + + nBas_Sq = nBas*nBas + +! Memory allocation + + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + + allocate(err(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,nOrb,guess_type,S,Hc,X,c) + + call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & + c(1,1), nBas, c(1,1), nBas, & + 0.d0, P(1,1), nBas) + + + ERI_size = (nBas * (nBas + 1)) / 2 + ERI_size = (ERI_size * (ERI_size + 1)) / 2 + allocate(ERI_chem(ERI_size)) + call read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem) + call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J) + + allocate(J_deb(nBas,nBas)) + allocate(ERI_phys(nBas,nBas,nBas,nBas)) + call read_2e_integrals(working_dir, nBas, ERI_phys) + call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb) + + print*, maxval(dabs(J - J_deb)) + diff = 0.d0 + do ii = 1, nBas + do jj = 1, nBas + diff_loc = dabs(J(jj,ii) - J_deb(jj,ii)) + if(diff_loc .gt. 1d-13) then + print*, 'error on: ', jj, ii + print*, J(jj,ii), J_deb(jj,ii) + stop + endif + diff = diff + diff_loc + enddo + enddo + print*, 'total diff = ', diff + + stop + +! Initialization + + n_diis = 0 + F_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + rcond = 0d0 + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + + write(*,*) + write(*,*)'-----------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(RHF)','|','EJ(RHF)','|','EK(RHF)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Build Fock matrix + call Hartree_matrix_AO_basis(nBas,P,ERI_phys,J) + call exchange_matrix_AO_basis(nBas,P,ERI_phys,K) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + + ! Check convergence + 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)) + + ! Potential energy + EV = trace_matrix(nBas, matmul(P, V)) + + ! Hartree energy + EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) + + ! Exchange energy + EK = 0.25d0*trace_matrix(nBas, matmul(P, K)) + + ! Total energy + ERHF = ET + EV + EJ + EK + + ! DIIS extrapolation + if(max_diis > 1) then + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) + endif + + ! Level shift + 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(nOrb,cp,eHF) + c = matmul(X,cp) + + ! Density matrix + 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 + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') & + '|',nSCF,'|',ERHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|' + + end do + write(*,*)'-----------------------------------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + deallocate(J,K,err,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,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,ERHF,dipole) + +! Testing zone + + if(dotest) then + + call dump_test_value('R','RHF energy',ERHF) + call dump_test_value('R','RHF HOMO energy',eHF(nO)) + call dump_test_value('R','RHF LUMO energy',eHF(nO+1)) + call dump_test_value('R','RHF dipole moment',norm2(dipole)) + + end if + + deallocate(J,K,err,cp,Fp,err_diis,F_diis) + +end subroutine diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 22441fd..b82f134 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -1,17 +1,19 @@ -subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, & - 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, & - 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_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & +subroutine GQuAcK(working_dir,dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, & + nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_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_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) implicit none include 'parameters.h' + character(len=256),intent(in) :: working_dir + logical,intent(in) :: dotest logical,intent(in) :: doGHF @@ -41,7 +43,6 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix @@ -86,9 +87,11 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do double precision :: start_GW ,end_GW ,t_GW double precision :: start_GT ,end_GT ,t_GT + double precision :: start_int, end_int, t_int double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:),FHF(:,:) double precision :: EGHF double precision,allocatable :: dipole_int_MO(:,:,:) + double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: ERI_tmp(:,:,:,:) double precision,allocatable :: Ca(:,:),Cb(:,:) @@ -112,6 +115,17 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do allocate(cHF(nBas2,nBas2),eHF(nBas2),PHF(nBas2,nBas2),FHF(nBas2,nBas2), & dipole_int_MO(nBas2,nBas2,ncart),ERI_MO(nBas2,nBas2,nBas2,nBas2)) + + allocate(ERI_AO(nBas,nBas,nBas,nBas)) + call wall_time(start_int) + call read_2e_integrals(working_dir,nBas,ERI_AO) + call wall_time(end_int) + t_int = end_int - start_int + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds' + write(*,*) + + !---------------------! ! Hartree-Fock module ! !---------------------! diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 9072904..8fd0214 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -32,7 +32,6 @@ program QuAcK double precision,allocatable :: Hc(:,:) double precision,allocatable :: X(:,:),X_tmp(:,:) double precision,allocatable :: dipole_int_AO(:,:,:) - double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: Uvec(:,:), Uval(:) double precision :: start_QuAcK,end_QuAcK,t_QuAcK @@ -44,6 +43,7 @@ program QuAcK logical :: reg_MP + logical :: switch_hpc logical :: use_gpu integer :: maxSCF_CC,max_diis_CC @@ -140,7 +140,7 @@ program QuAcK ! Hardware ! !------------------! - call read_hardware(working_dir,use_gpu) + call read_hpc_flags(working_dir,switch_hpc,use_gpu) !------------------------------------! ! Read input information ! @@ -176,20 +176,19 @@ program QuAcK 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(working_dir,nBas,S,T,V,Hc,ERI_AO) + call read_1e_integrals(working_dir,nBas,S,T,V,Hc) call read_dipole_integrals(working_dir,nBas,dipole_int_AO) call wall_time(end_int) t_int = end_int - start_int write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 1e-integrals = ',t_int,' seconds' write(*,*) ! Compute orthogonalization matrix @@ -225,29 +224,44 @@ program QuAcK ! Restricted QuAcK branch ! !-------------------------! - if(doRQuAcK) & - call RQuAcK(use_gpu,doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - 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, & - doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - 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, & - TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) + if(doRQuAcK) then + + if(switch_hpc) then + call RQuAcK_hpc(working_dir,use_gpu,doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + 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, & + doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + S,T,V,Hc,X,dipole_int_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, & + TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) + else + call RQuAcK(working_dir,use_gpu,doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + 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, & + doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + S,T,V,Hc,X,dipole_int_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, & + TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) + endif + endif !---------------------------! ! Unrestricted QuAcK branch ! !---------------------------! if(doUQuAcK) & - call UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + call UQuAcK(working_dir,doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & 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, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & 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, & + S,T,V,Hc,X,dipole_int_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, & TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & @@ -257,10 +271,10 @@ program QuAcK ! Generalized QuAcK branch ! !--------------------------! if(doGQuAcK) & - call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + call GQuAcK(working_dir,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,doG0T0pp,doevGTpp,doqsGTpp, & - nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & + nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_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, & @@ -289,4 +303,10 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds' write(*,*) + deallocate(S) + deallocate(T) + deallocate(V) + deallocate(Hc) + deallocate(dipole_int_AO) + end program diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 4a03be8..16b8159 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -1,12 +1,12 @@ -subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - 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, & - doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - 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, & - TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & +subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + 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, & + doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + S,T,V,Hc,X,dipole_int_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, & + TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) ! Restricted branch of QuAcK @@ -14,6 +14,8 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD, implicit none include 'parameters.h' + character(len=256),intent(in) :: working_dir + logical,intent(in) :: use_gpu logical,intent(in) :: dotest @@ -46,7 +48,6 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD, double precision,intent(in) :: Hc(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) integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix @@ -94,12 +95,14 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD, double precision :: start_GW ,end_GW ,t_GW double precision :: start_GT ,end_GT ,t_GT + double precision :: start_int, end_int, t_int double precision,allocatable :: eHF(:) double precision,allocatable :: cHF(:,:) double precision,allocatable :: PHF(:,:) double precision,allocatable :: FHF(:,:) double precision :: ERHF double precision,allocatable :: dipole_int_MO(:,:,:) + double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz integer :: nS @@ -121,6 +124,15 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD, allocate(dipole_int_MO(nOrb,nOrb,ncart)) allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) + allocate(ERI_AO(nBas,nBas,nBas,nBas)) + call wall_time(start_int) + call read_2e_integrals(working_dir,nBas,ERI_AO) + call wall_time(end_int) + t_int = end_int - start_int + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds' + write(*,*) + !---------------------! ! Hartree-Fock module ! !---------------------! @@ -350,4 +362,13 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD, end if + + deallocate(eHF) + deallocate(cHF) + deallocate(PHF) + deallocate(FHF) + deallocate(dipole_int_MO) + deallocate(ERI_MO) + deallocate(ERI_AO) + end subroutine diff --git a/src/QuAcK/RQuAcK_hpc.f90 b/src/QuAcK/RQuAcK_hpc.f90 new file mode 100644 index 0000000..9d13a57 --- /dev/null +++ b/src/QuAcK/RQuAcK_hpc.f90 @@ -0,0 +1,353 @@ +subroutine RQuAcK_hpc(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + 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, & + doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + S,T,V,Hc,X,dipole_int_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, & + TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) + +! Restricted branch of QuAcK + + implicit none + include 'parameters.h' + + character(len=256),intent(in) :: working_dir + + logical,intent(in) :: use_gpu + + logical,intent(in) :: dotest + + logical,intent(in) :: doRHF,doROHF + logical,intent(in) :: dostab + logical,intent(in) :: dosearch + logical,intent(in) :: doMP2,doMP3 + logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD + logical,intent(in) :: doCIS,doCIS_D,doCID,doCISD,doFCI + logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA + logical,intent(in) :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3 + logical,intent(in) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW + logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp + logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh + + integer,intent(in) :: nNuc,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) :: ZNuc(nNuc),rNuc(nNuc,ncart) + + 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,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + + integer,intent(in) :: maxSCF_HF,max_diis_HF + double precision,intent(in) :: thresh_HF,level_shift,mix + integer,intent(in) :: guess_type + + logical,intent(in) :: reg_MP + + integer,intent(in) :: maxSCF_CC,max_diis_CC + double precision,intent(in) :: thresh_CC + + logical,intent(in) :: singlet + logical,intent(in) :: triplet + logical,intent(in) :: TDA + + integer,intent(in) :: maxSCF_GF,max_diis_GF,renorm_GF + double precision,intent(in) :: thresh_GF + logical,intent(in) :: lin_GF,reg_GF + double precision,intent(in) :: eta_GF + + integer,intent(in) :: maxSCF_GW,max_diis_GW + double precision,intent(in) :: thresh_GW + logical,intent(in) :: TDA_W,lin_GW,reg_GW + double precision,intent(in) :: eta_GW + + integer,intent(in) :: maxSCF_GT,max_diis_GT + double precision,intent(in) :: thresh_GT + logical,intent(in) :: TDA_T,lin_GT,reg_GT + double precision,intent(in) :: eta_GT + + logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA + logical,intent(in) :: doACFDT,exchange_kernel,doXBS + +! Local variables + + logical :: doMP,doCC,doCI,doRPA,doGF,doGW,doGT + + double precision :: start_HF ,end_HF ,t_HF + double precision :: start_stab ,end_stab ,t_stab + double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO + double precision :: start_MP ,end_MP ,t_MP + double precision :: start_CC ,end_CC ,t_CC + double precision :: start_CI ,end_CI ,t_CI + double precision :: start_RPA ,end_RPA ,t_RPA + double precision :: start_GF ,end_GF ,t_GF + double precision :: start_GW ,end_GW ,t_GW + double precision :: start_GT ,end_GT ,t_GT + + double precision :: start_int, end_int, t_int + double precision,allocatable :: eHF(:) + double precision,allocatable :: cHF(:,:) + double precision,allocatable :: PHF(:,:) + double precision,allocatable :: FHF(:,:) + double precision :: ERHF + integer :: ixyz + integer :: nS + + write(*,*) + write(*,*) '*****************************************' + write(*,*) '* Restricted Branch of QuAcK (HPC mode) *' + write(*,*) '*****************************************' + write(*,*) + +!-------------------! +! Memory allocation ! +!-------------------! + + allocate(eHF(nOrb)) + allocate(cHF(nBas,nOrb)) + allocate(PHF(nBas,nBas)) + allocate(FHF(nBas,nBas)) + +!---------------------! +! Hartree-Fock module ! +!---------------------! + + if(doRHF) then + + call wall_time(start_HF) + call RHF_hpc(working_dir,dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF) + call wall_time(end_HF) + + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' + write(*,*) + + end if + +! if(doROHF) then +! +! call wall_time(start_HF) +! call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & +! nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF) +! call wall_time(end_HF) +! +! t_HF = end_HF - start_HF +! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ROHF = ',t_HF,' seconds' +! write(*,*) +! +! end if +! +!!----------------------------------! +!! AO to MO integral transformation ! +!!----------------------------------! +! +! call wall_time(start_AOtoMO) +! +! write(*,*) +! write(*,*) 'AO to MO transformation... Please be patient' +! write(*,*) +! +! ! Read and transform dipole-related integrals +! +! 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,nOrb,cHF,ERI_AO,ERI_MO) +! +! call wall_time(end_AOtoMO) +! +! t_AOtoMO = end_AOtoMO - start_AOtoMO +! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds' +! write(*,*) +! +!!-----------------------------------! +!! Stability analysis of HF solution ! +!!-----------------------------------! +! +! nS = (nO - nC)*(nV - nR) +! +! if(dostab) then +! +! call wall_time(start_stab) +! 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 wall time for stability analysis = ',t_stab,' seconds' +! write(*,*) +! +! end if +! +! 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,nOrb,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X, & +! ERHF,eHF,cHF,PHF,FHF) +! call wall_time(end_stab) +! +! t_stab = end_stab - start_stab +! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds' +! write(*,*) +! +! end if +! +!!-----------------------! +!! Moller-Plesset module ! +!!-----------------------! +! +! doMP = doMP2 .or. doMP3 +! +! if(doMP) then +! +! call wall_time(start_MP) +! call RMP(dotest,doMP2,doMP3,reg_MP,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 wall time for MP = ',t_MP,' seconds' +! write(*,*) +! +! end if +! +!!------------------------! +!! Coupled-cluster module ! +!!------------------------! +! +! doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT .or. & +! dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD +! +! if(doCC) then +! +! call wall_time(start_CC) +! call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & +! maxSCF_CC,thresh_CC,max_diis_CC,nBas,nOrb,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO, & +! ENuc,ERHF,eHF,cHF,PHF,FHF) +! call wall_time(end_CC) +! +! t_CC = end_CC - start_CC +! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CC = ',t_CC,' seconds' +! write(*,*) +! +! end if +! +!!----------------------------------! +!! Configuration interaction module ! +!!----------------------------------! +! +! doCI = doCIS .or. doCID .or. doCISD .or. doFCI +! +! if(doCI) then +! +! call wall_time(start_CI) +! 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 wall time for CI = ',t_CI,' seconds' +! write(*,*) +! +! end if +! +!!-----------------------------------! +!! Random-phase approximation module ! +!!-----------------------------------! +! +! doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA +! +! if(doRPA) then +! +! call wall_time(start_RPA) +! call RRPA(use_gpu,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 wall time for RPA = ',t_RPA,' seconds' +! write(*,*) +! +! end if +! +!!-------------------------! +!! Green's function module ! +!!-------------------------! +! +! doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 +! +! 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,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 wall time for GF2 = ',t_GF,' seconds' +! write(*,*) +! +! end if +! +!!-----------! +!! GW module ! +!!-----------! +! +! doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW +! +! if(doGW) then +! +! call wall_time(start_GW) +! call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,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 wall time for GW = ',t_GW,' seconds' +! write(*,*) +! +! end if +! +!!-----------------! +!! T-matrix module ! +!!-----------------! +! +! doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp .or. doufG0T0pp .or. doG0T0eh .or. doevGTeh .or. doqsGTeh +! +! 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,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 wall time for GT = ',t_GT,' seconds' +! write(*,*) +! +! end if + + + return +end subroutine diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index f813a72..0ff1759 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -1,9 +1,9 @@ -subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & +subroutine UQuAcK(working_dir,dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & 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, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & 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, & + S,T,V,Hc,X,dipole_int_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, & TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & @@ -12,6 +12,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do implicit none include 'parameters.h' + character(len=256),intent(in) :: working_dir + logical,intent(in) :: dotest logical,intent(in) :: doUHF @@ -42,7 +44,6 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix @@ -90,10 +91,12 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do double precision :: start_GW ,end_GW ,t_GW double precision :: start_GT ,end_GT ,t_GT + double precision :: start_int, end_int, t_int double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:),FHF(:,:,:) double precision :: EUHF double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:) double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:) + double precision,allocatable :: ERI_AO(:,:,:,:) integer :: ixyz integer :: nS(nspin) @@ -112,6 +115,15 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas), & ERI_bbbb(nBas,nBas,nBas,nBas)) + allocate(ERI_AO(nBas,nBas,nBas,nBas)) + call wall_time(start_int) + call read_2e_integrals(working_dir,nBas,ERI_AO) + call wall_time(end_int) + t_int = end_int - start_int + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds' + write(*,*) + !---------------------! ! Hartree-Fock module ! !---------------------! diff --git a/src/QuAcK/read_hardware.f90 b/src/QuAcK/read_hardware.f90 deleted file mode 100644 index c014d44..0000000 --- a/src/QuAcK/read_hardware.f90 +++ /dev/null @@ -1,45 +0,0 @@ -subroutine read_hardware(working_dir,use_gpu) - -! Read desired methods - - implicit none - -! Input variables - - character(len=256),intent(in) :: working_dir - -! Output variables - - logical,intent(out) :: use_gpu - -! Local variables - - character(len=1) :: ans - integer :: ios - character(len=256) :: file_path - -! Open file with method specification - - file_path = trim(working_dir) // '/input/hardware' - open(unit=1, file=file_path, status='old', action='read', iostat=ios) - - if(ios /= 0) then - - use_gpu = .False. - - else - - read(1,*) - read(1,*) ans - if(ans == 'T') then - use_gpu = .true. - else - use_gpu = .False. - endif - - endif - - ! Close file with options - close(unit=1) - -end subroutine diff --git a/src/QuAcK/read_hpc_flags.f90 b/src/QuAcK/read_hpc_flags.f90 new file mode 100644 index 0000000..4921754 --- /dev/null +++ b/src/QuAcK/read_hpc_flags.f90 @@ -0,0 +1,39 @@ +subroutine read_hpc_flags(working_dir, switch_hpc, use_gpu) + + implicit none + + character(len=256), intent(in) :: working_dir + + logical, intent(out) :: switch_hpc + logical, intent(out) :: use_gpu + + character(len=1) :: ans + integer :: ios + character(len=256) :: file_path + + file_path = trim(working_dir) // '/input/hpc_flags' + open(unit=1, file=file_path, status='old', action='read', iostat=ios) + + if(ios /= 0) then + + switch_hpc = .False. + use_gpu = .False. + + else + + switch_hpc = .False. + use_gpu = .False. + + read(1,*) + read(1,*) ans + if(ans == 'T') switch_hpc = .true. + + read(1,*) + read(1,*) ans + if(ans == 'T') use_gpu = .true. + + endif + + close(unit=1) + +end subroutine diff --git a/src/utils/read_integrals.f90 b/src/utils/read_1e_integrals.f90 similarity index 59% rename from src/utils/read_integrals.f90 rename to src/utils/read_1e_integrals.f90 index 91be5fe..eecdc05 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_1e_integrals.f90 @@ -1,6 +1,6 @@ -subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) +subroutine read_1e_integrals(working_dir,nBas_AOs,S,T,V,Hc) -! Read one- and two-electron integrals from files +! Read one-electron integrals from files implicit none include 'parameters.h' @@ -13,9 +13,8 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) ! Local variables logical :: debug - integer :: mu,nu,la,si - double precision :: Ov,Kin,Nuc,ERI - double precision :: lambda + integer :: mu,nu + double precision :: Ov,Kin,Nuc ! Output variables @@ -23,26 +22,21 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) 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) - integer :: status, ios + integer :: ios character(len=256) :: file_path ! Open file with integrals debug = .false. - lambda = 1d0 - - print*, 'Scaling integrals by ',lambda - ! --- ! Read overlap integrals file_path = trim(working_dir) // '/int/Ov.dat' - open(unit=8, file=file_path, status='old', action='read', iostat=status) - if(status /= 0) then + open(unit=8, file=file_path, status='old', action='read', iostat=ios) + if(ios /= 0) then print *, "Error opening file: ", file_path stop else @@ -60,8 +54,8 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) ! Read kinetic integrals file_path = trim(working_dir) // '/int/Kin.dat' - open(unit=9, file=file_path, status='old', action='read', iostat=status) - if(status /= 0) then + open(unit=9, file=file_path, status='old', action='read', iostat=ios) + if(ios /= 0) then print *, "Error opening file: ", file_path stop else @@ -79,8 +73,8 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) ! Read nuclear integrals file_path = trim(working_dir) // '/int/Nuc.dat' - open(unit=10, file=file_path, status='old', action='read', iostat=status) - if(status /= 0) then + open(unit=10, file=file_path, status='old', action='read', iostat=ios) + if(ios /= 0) then print *, "Error opening file: ", file_path stop else @@ -99,37 +93,6 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) ! Define core Hamiltonian Hc(:,:) = T(:,:) + V(:,:) -! Read 2e-integrals - -! ! 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) - - ! binary file - file_path = trim(working_dir) // '/int/ERI.bin' - open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=status) - if(status /= 0) then - print *, "Error opening file: ", file_path - stop - else - read(11) G - endif - close(unit=11) - - ! Print results if(debug) then @@ -148,15 +111,6 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) write(*,'(A28)') '----------------------' call matout(nBas_AOs,nBas_AOs,V) write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Electron repulsion integrals' - write(*,'(A28)') '----------------------' - 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(*,*) end if end subroutine diff --git a/src/utils/read_2e_integrals.f90 b/src/utils/read_2e_integrals.f90 new file mode 100644 index 0000000..94404a4 --- /dev/null +++ b/src/utils/read_2e_integrals.f90 @@ -0,0 +1,110 @@ +subroutine read_2e_integrals(working_dir,nBas_AOs,G) + +! Read two-electron integrals from files + + implicit none + +! Input variables + + integer,intent(in) :: nBas_AOs + character(len=256),intent(in) :: working_dir + +! Local variables + + logical :: debug + integer :: mu,nu,la,si + double precision :: ERI + double precision :: lambda + +! Output variables + + double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + + integer :: ios + character(len=256) :: file_path + +! Open file with integrals + + debug = .false. + + lambda = 1d0 + + print*, 'Scaling integrals by ',lambda + + +! Read 2e-integrals + +! ! 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) + + ! binary file + file_path = trim(working_dir) // '/int/ERI.bin' + open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=ios) + if(ios /= 0) then + print *, "Error opening file: ", file_path + stop + else + read(11) G + endif + close(unit=11) + G = G * lambda + + + +! Print results + if(debug) then + write(*,'(A28)') '----------------------' + write(*,'(A28)') 'Electron repulsion integrals' + write(*,'(A28)') '----------------------' + 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(*,*) + end if + +end subroutine + +! --- + +subroutine read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem) + + implicit none + + character(len=256), intent(in) :: working_dir + integer*8, intent(in) :: ERI_size + double precision, intent(out) :: ERI_chem(ERI_size) + + integer :: ios + character(len=256) :: file_path + + file_path = trim(working_dir) // '/int/ERI_chem.bin' + open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=ios) + if(ios /= 0) then + print *, "Error opening file: ", file_path + stop + else + read(11) ERI_chem + endif + close(unit=11) + + return +end subroutine + +! --- +