Working on s8 rep of 2e-integrals

PyDuck.py

@@ -24,8 +24,10 @@
 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 @@
 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 @@
     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()
+
 
 
 

input/hpc_flags

@@ -0,0 +1,4 @@
+# if True (T), switch to HPC mode
+  F
+# if True (T), use GPU
+  F

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 
+
+! ---
+