Merge pull request #289 from QuantumPackage/dev-stable

Improve 5idx

scripts/qp_import_trexio.py

@@ -17,6 +17,7 @@ import numpy as np
 from functools import reduce
 from ezfio import ezfio
 from docopt import docopt
+import qp_bitmasks
 
 try:
   import trexio
@@ -453,6 +454,20 @@ def write_ezfio(trexio_filename, filename):
     else:
         print("None")
 
+    print("Determinant\t\t...\t", end=' ')
+    alpha = [ i for i in range(num_alpha) ]
+    beta  = [ i for i in range(num_beta) ]
+    if trexio.has_mo_spin(trexio_file):
+       spin = trexio.read_mo_spin(trexio_file)
+       beta  = [ i for i in range(mo_num) if spin[i] == 1 ]
+       beta  = [ beta[i] for i in range(num_beta) ]
+
+    alpha = qp_bitmasks.BitMask(alpha)
+    beta  = qp_bitmasks.BitMask(beta )
+    print(alpha)
+    print(beta)
+    print("OK")
+
 
 
 

scripts/utility/qp_bitmasks.py

@@ -22,7 +22,7 @@ def int_to_string(s):
         assert s>=0
     AssertionError
     """
-    assert type(s) in (int, long)
+    assert type(s) == int
     assert s>=0
     return '{s:0b}'.format(s=s)
 
@@ -62,7 +62,7 @@ def int_to_bitmask(s,bit_kind_size=BIT_KIND_SIZE):
     ['1111111111111111111111111111111111111111111111111111111111110110']
     >>>
     """
-    assert type(s) in (int, long)
+    assert type(s) == int
     if s < 0:
         s = s + (1 << bit_kind_size)
     return ['{s:0{width}b}'.format(s=s,width=bit_kind_size)]
@@ -104,7 +104,7 @@ class BitMask(object):
     return self._data_int[i]
 
   def __setitem__(self,i,value):
-    if type(value) in (int,long):
+    if type(value) == int :
         self._data_int[i] = value
     elif type(value) == str:
         s = string_to_bitmask(value,bit_kind_size=self.bit_kind_size)[0]

src/ao_basis/cosgtos.irp.f

@@ -6,13 +6,14 @@ BEGIN_PROVIDER [ logical, use_cosgtos  ]
 
   logical                        :: has
   PROVIDE ezfio_filename
+  use_cosgtos = .False.
   if (mpi_master) then
     call ezfio_has_ao_basis_use_cosgtos(has)
     if (has) then
 !      write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
       call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
     else
-      use_cosgtos = .False.
+      call ezfio_set_ao_basis_use_cosgtos(use_cosgtos)
     endif
   endif
   IRP_IF MPI_DEBUG

src/bi_ort_ints/bi_ort_ints.irp.f

@@ -55,6 +55,7 @@ subroutine test_5idx
  implicit none
  integer :: i,k,j,l,m,n,ipoint
  double precision :: accu, contrib,new,ref
+ double precision, external :: three_e_5_idx_exch12_bi_ort
  i = 1
  k = 1
  n = 0
@@ -64,18 +65,21 @@ subroutine test_5idx
    do j = 1, mo_num
     do l = 1, mo_num
      do m = 1, mo_num
+!      if (dabs(three_e_5_idx_direct_bi_ort(m,l,j,k,i) - three_e_5_idx_exch12_bi_ort(m,l,i,k,j)) > 1.d-10) then 
+!         stop
+!      endif
 
-      new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
+!      new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
+!      ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
+!      contrib = dabs(new - ref)
-      accu += contrib
+!      accu += contrib
-      if(contrib .gt. 1.d-10)then
+!      if(contrib .gt. 1.d-10)then
-       print*,'direct'
+!       print*,'direct'
-       print*,i,k,j,l,m
+!       print*,i,k,j,l,m
-       print*,ref,new,contrib
+!       print*,ref,new,contrib
-       stop
+!       stop
-      endif
+!      endif
-
+!
       new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
       ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
       contrib = dabs(new - ref)
@@ -86,51 +90,52 @@ subroutine test_5idx
        print*,ref,new,contrib
        stop
       endif
+
+!
+!      new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
+!      ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
+!      contrib = dabs(new - ref)
+!      accu += contrib
+!      if(contrib .gt. 1.d-10)then
+!       print*,'cycle1'
+!       print*,i,k,j,l,m
+!       print*,ref,new,contrib
+!       stop
+!      endif
+!
+!      new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
+!      ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
+!      contrib = dabs(new - ref)
+!      accu += contrib
+!      if(contrib .gt. 1.d-10)then
+!       print*,'cycle2'
+!       print*,i,k,j,l,m
+!       print*,ref,new,contrib
+!       stop
+!      endif
+!
+!      new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
+!      ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
+!      contrib = dabs(new - ref)
+!      accu += contrib
+!      if(contrib .gt. 1.d-10)then
+!       print*,'exch23'
+!       print*,i,k,j,l,m
+!       print*,ref,new,contrib
+!       stop
+!      endif
+!
+!      new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
+!      ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
+!      contrib = dabs(new - ref)
+!      accu += contrib
+!      if(contrib .gt. 1.d-10)then
+!       print*,'exch13'
+!       print*,i,k,j,l,m
+!       print*,ref,new,contrib
+!       stop
+!      endif
 !
-      new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
-      accu += contrib
-      if(contrib .gt. 1.d-10)then
-       print*,'cycle1'
-       print*,i,k,j,l,m
-       print*,ref,new,contrib
-       stop
-      endif
-
-      new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
-      accu += contrib
-      if(contrib .gt. 1.d-10)then
-       print*,'cycle2'
-       print*,i,k,j,l,m
-       print*,ref,new,contrib
-       stop
-      endif
-
-      new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
-      accu += contrib
-      if(contrib .gt. 1.d-10)then
-       print*,'exch23'
-       print*,i,k,j,l,m
-       print*,ref,new,contrib
-       stop
-      endif
-
-      new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
-      ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
-      contrib = dabs(new - ref)
-      accu += contrib
-      if(contrib .gt. 1.d-10)then
-       print*,'exch13'
-       print*,i,k,j,l,m
-       print*,ref,new,contrib
-       stop
-      endif
-
      enddo
     enddo
    enddo

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -1,7 +1,11 @@
 ! ---
+double precision function three_e_5_idx_exch12_bi_ort(m,l,i,k,j) result(integral)
+  implicit none
+  integer, intent(in) :: m,l,j,k,i
+  integral = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
+end
 
  BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
-&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
 &BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
 &BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
 &BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
@@ -14,6 +18,8 @@
   ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
   !
   ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
+  !
+
   END_DOC
 
   implicit none
@@ -74,7 +80,6 @@
       do j = 1, mo_num
         do l = 1, mo_num
             three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j)
-            three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,k,j) - tmp_mat(k,j,l,i)
         enddo
       enddo
     enddo
@@ -128,7 +133,6 @@
       do j = 1, mo_num
         do l = 1, mo_num
             three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i)
-            three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
         enddo
       enddo
     enddo

src/tc_bi_ortho/31.tc_bi_ortho.bats

@@ -4,12 +4,16 @@ source $QP_ROOT/tests/bats/common.bats.sh
 source $QP_ROOT/quantum_package.rc
 
 
+function get_e() {
+  grep "eigval_right_tc_bi_orth" $1 | cut -d '=' -f 2 | xargs
+}
+
 function run_Ne() {
   qp set_file Ne_tc_scf
   qp run cisd
   qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out
   eref=-128.77020441279302
-  energy="$(grep "eigval_right_tc_bi_orth =" Ne_tc_scf.cisd_tc_bi_ortho.out)"
+  energy=$(get_e Ne_tc_scf.cisd_tc_bi_ortho.out)
   eq $energy $eref 1e-6
 }
 
@@ -24,7 +28,7 @@ function run_C() {
   qp run cisd
   qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out
   eref=-37.757536149952514
-  energy="$(grep "eigval_right_tc_bi_orth =" C_tc_scf.cisd_tc_bi_ortho.out)"
+  energy=$(get_e C_tc_scf.cisd_tc_bi_ortho.out)
   eq $energy $eref 1e-6
 }
 
@@ -38,7 +42,7 @@ function run_O() {
   qp run cisd
   qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out
   eref=-74.908518517716161
-  energy="$(grep "eigval_right_tc_bi_orth =" O_tc_scf.cisd_tc_bi_ortho.out)"
+  energy=$(get_e O_tc_scf.cisd_tc_bi_ortho.out)
   eq $energy $eref 1e-6
 }
 

src/tc_bi_ortho/slater_tc_3e_slow.irp.f

@@ -242,7 +242,8 @@ subroutine double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
       do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1)
        mm = occ(m,ispin)
        direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
-       exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+!       exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+       exchange_int = three_e_5_idx_direct_bi_ort(mm,p2,h1,p1,h2)
        hthree += direct_int - exchange_int
       enddo
       do m = 1, Ne(s1) ! pure contribution from s1

src/tc_bi_ortho/slater_tc_opt.irp.f

@@ -13,8 +13,7 @@ subroutine provide_all_three_ints_bi_ortho
    PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort
   endif
   if(.not.double_normal_ord.and.three_e_5_idx_term)then
-   PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort
+   PROVIDE three_e_5_idx_direct_bi_ort 
-   PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
   elseif (double_normal_ord .and. (.not. three_e_5_idx_term))then
    PROVIDE normal_two_body_bi_orth
   endif

src/tc_bi_ortho/slater_tc_opt_double.irp.f

@@ -136,11 +136,13 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
      ! exchange between (h1,p1) and (h2,p2)
      ipart=occ_particle(i,ispin)
      direct_int  = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
-     exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
+!     exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
+     exchange_int = three_e_5_idx_direct_bi_ort(ipart,p2,h1,p1,h2)
      hthree += direct_int - exchange_int
      ihole=occ_hole(i,ispin)
      direct_int  = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
-     exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
+!     exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
+     exchange_int = three_e_5_idx_direct_bi_ort(ihole,p2,h1,p1,h2)
      hthree -= direct_int - exchange_int
     enddo
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -163,11 +165,13 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
     ! exchange between (h1,p1) and (h2,p2)
     ipart=occ_particle(i,ispin)
     direct_int  = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
-    exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
+!    exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_direct_bi_ort(ipart,p2,h1,p1,h2)
     hthree += direct_int - exchange_int
     ihole=occ_hole(i,ispin)
     direct_int  = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
-    exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
+!    exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_direct_bi_ort(ihole,p2,h1,p1,h2)
     hthree -= direct_int - exchange_int
    enddo
   else                         ! (h1,p1) == alpha/(h2,p2) == beta
@@ -342,7 +346,8 @@ subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib)
    m = occ(mm,2)
    direct_int   = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
    ! exchange between (h1,p1) and (h2,p2)
-   exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+!   exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+   exchange_int = three_e_5_idx_direct_bi_ort(mm,p2,h1,p1,h2)
    contrib += direct_int - exchange_int
  enddo
 end
@@ -415,7 +420,8 @@ subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
    m = occ(mm,1)
    direct_int   = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
    ! exchange between (h1,p1) and (h2,p2)
-   exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+!   exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
+   exchange_int = three_e_5_idx_direct_bi_ort(mm,p2,h1,p1,h2)
    contrib += direct_int - exchange_int
  enddo
 end

src/tc_bi_ortho/symmetrized_3_e_int.irp.f

@@ -107,5 +107,6 @@ double precision function three_e_double_parrallel_spin(m,l,j,k,i)
   three_e_double_parrallel_spin = three_e_5_idx_direct_bi_ort(m,l,j,k,i)  ! direct
   three_e_double_parrallel_spin += three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) & ! two cyclic permutations 
   - three_e_5_idx_exch23_bi_ort(m,l,j,k,i) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i)  & ! two first exchange 
-  - three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange 
+!  - three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange 
+  - three_e_5_idx_direct_bi_ort(m,l,i,k,j) ! last exchange 
 end