cas_ful -> cas_full

ocaml/qptypes_generator.ml

@@ -154,8 +154,8 @@ let input_ezfio = "
 
 * N_int_number : int
   determinants_n_int
-  1 : 30
-  N_int > 30
+  1 : 128
+  N_int > 128
 
 * Det_number : int
   determinants_n_det

plugins/local/basis_correction/README.rst

@@ -12,7 +12,7 @@ This basis set correction relies mainy on :
        When HF is a qualitative representation of the electron pairs (i.e. weakly correlated systems), such an approach for \mu(r) is OK. 
        See for instance JPCL, 10, 2931-2937 (2019) for typical flavours of the results. 
        Thanks to the trivial nature of such a two-body rdm, the equation (22) of J. Chem. Phys. 149, 194301 (2018) can be rewritten in a very efficient way, and therefore the limiting factor of such an approach is the AO->MO four-index transformation of the two-electron integrals.  
-    b) "mu_of_r_potential = cas_ful" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder). 
+    b) "mu_of_r_potential = cas_full" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder). 
       If the CAS is properly chosen (i.e. the CAS-like wave function qualitatively represents the wave function of the systems), then such an approach is OK for \mu(r) even in the case of strong correlation. 
 
   +) The use of DFT correlation functionals with multi-determinant reference (Ecmd). These functionals are originally defined in the RS-DFT framework (see for instance Theor.  Chem.  Acc.114,  305(2005)) and design to capture short-range correlation effects. A important quantity arising in the Ecmd is the exact on-top pair density of the system, and the main differences of approximated Ecmd relies on different approximations for the exact on-top pair density.  

plugins/local/basis_correction/pbe_on_top.irp.f

@@ -39,7 +39,7 @@
    grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
    grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
    
-   if(mu_of_r_potential == "cas_ful")then
+   if(mu_of_r_potential == "cas_full")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
     on_top =  on_top_cas_mu_r(ipoint,istate)
    else
@@ -101,7 +101,7 @@
    grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
    grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
    
-   if(mu_of_r_potential == "cas_ful")then
+   if(mu_of_r_potential == "cas_full")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
     on_top = on_top_cas_mu_r(ipoint,istate)
    else
@@ -163,7 +163,7 @@
    grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
    grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
    
-   if(mu_of_r_potential == "cas_ful")then
+   if(mu_of_r_potential == "cas_full")then
     ! You take the on-top of the CAS wave function which is computed with mu(r) 
     on_top = on_top_cas_mu_r(ipoint,istate)
    else

plugins/local/basis_correction/print_routine.irp.f

@@ -4,8 +4,8 @@ subroutine print_basis_correction
  provide mu_average_prov
  if(mu_of_r_potential.EQ."hf")then
   provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
- else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
-  provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r 
+ else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated")then
+  provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
   provide ecmd_pbe_on_top_mu_of_r ecmd_pbe_on_top_su_mu_of_r
  endif
 
@@ -25,7 +25,7 @@ subroutine print_basis_correction
  if(mu_of_r_potential.EQ."hf")then
    print*, ''
    print*,'Using a HF-like two-body density to define mu(r)'
-   print*,'This assumes that HF is a qualitative representation of the wave function ' 
+   print*,'This assumes that HF is a qualitative representation of the wave function '
    print*,'********************************************'
    print*,'Functionals more suited for weak correlation'
    print*,'********************************************'
@@ -38,10 +38,10 @@ subroutine print_basis_correction
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-UEG       , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
    enddo
 
-  else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
+  else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
    print*, ''
    print*,'Using a CAS-like two-body density to define mu(r)'
-   print*,'This assumes that the CAS is a qualitative representation of the wave function ' 
+   print*,'This assumes that the CAS is a qualitative representation of the wave function '
    print*,'********************************************'
    print*,'Functionals more suited for weak correlation'
    print*,'********************************************'
@@ -56,14 +56,14 @@ subroutine print_basis_correction
    print*,''
    print*,'********************************************'
    print*,'********************************************'
-   print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) ' 
+   print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) '
    print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
    do istate = 1, N_states
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-OT        , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
    enddo
    print*,''
    print*,'********************************************'
-   print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)' 
+   print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)'
    print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
    do istate = 1, N_states
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD SU-PBE-OT     , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)

scripts/qp_import_trexio.py

@@ -142,6 +142,7 @@ def write_ezfio(trexio_filename, filename):
     try:
         basis_type = trexio.read_basis_type(trexio_file)
 
+        print ("BASIS TYPE: ", basis_type.lower())
         if basis_type.lower() in ["gaussian", "slater"]:
             shell_num   = trexio.read_basis_shell_num(trexio_file)
             prim_num    = trexio.read_basis_prim_num(trexio_file)

src/hartree_fock/fock_matrix_hf.irp.f

@@ -174,7 +174,6 @@ END_PROVIDER
 
  allocate (X(cholesky_ao_num))
 
-
  ! X(j) = \sum_{mn} SCF_density_matrix_ao(m,n) * cholesky_ao(m,n,j)
  call dgemm('T','N',cholesky_ao_num,1,ao_num*ao_num,1.d0,            &
      cholesky_ao, ao_num*ao_num,                                     &

src/mo_two_e_ints/mo_bi_integrals_erf.irp.f

@@ -31,37 +31,144 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_erf_in_map ]
 
   PROVIDE mo_class
 
-  real                           :: map_mb
-
   mo_two_e_integrals_erf_in_map = .True.
   if (read_mo_two_e_integrals_erf) then
     print*,'Reading the MO integrals_erf'
     call map_load_from_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
     print*, 'MO integrals_erf provided'
     return
-  else
-    PROVIDE ao_two_e_integrals_erf_in_map
   endif
 
-  !  call four_index_transform_block(ao_integrals_erf_map,mo_integrals_erf_map, &
-  !      mo_coef, size(mo_coef,1),                                      &
-  !      1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num,                    &
-  !      1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
-   call add_integrals_to_map_erf(full_ijkl_bitmask_4)
-    integer*8                      :: get_mo_erf_map_size, mo_erf_map_size
-    mo_erf_map_size = get_mo_erf_map_size()
+  PROVIDE ao_two_e_integrals_erf_in_map
 
-! print*,'Molecular integrals ERF provided:'
-! print*,' Size of MO ERF map           ', map_mb(mo_integrals_erf_map) ,'MB'
-! print*,' Number of MO ERF integrals: ',  mo_erf_map_size
-  if (write_mo_two_e_integrals_erf) then
+  print *,  ''
+  print *,  'AO -> MO ERF integrals transformation'
+  print *,  '-------------------------------------'
+  print *,  ''
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+      if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
+        call four_idx_dgemm_erf
+      else
+        call add_integrals_to_map_erf(full_ijkl_bitmask_4)
+      endif
+
+  call wall_time(wall_2)
+  call cpu_time(cpu_2)
+
+  integer*8                      :: get_mo_erf_map_size, mo_erf_map_size
+  mo_erf_map_size = get_mo_erf_map_size()
+
+  double precision, external     :: map_mb
+  print*,'Molecular integrals provided:'
+  print*,' Size of MO map           ', map_mb(mo_integrals_erf_map) ,'MB'
+  print*,' Number of MO integrals: ',  mo_erf_map_size
+  print*,' cpu  time :',cpu_2 - cpu_1, 's'
+  print*,' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
+
+  if (write_mo_two_e_integrals_erf.and.mpi_master) then
     call ezfio_set_work_empty(.False.)
     call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
-    call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals_erf("Read")
+    call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals_erf('Read')
   endif
 
 END_PROVIDER
 
+subroutine four_idx_dgemm_erf
+  implicit none
+  integer :: p,q,r,s,i,j,k,l
+  double precision, allocatable :: a1(:,:,:,:)
+  double precision, allocatable :: a2(:,:,:,:)
+
+  if (ao_num > 1289) then
+    print *,  irp_here, ': Integer overflow in ao_num**3'
+  endif
+
+  allocate (a1(ao_num,ao_num,ao_num,ao_num))
+
+  print *, 'Getting AOs'
+  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,r,s)
+  do s=1,ao_num
+    do r=1,ao_num
+      do q=1,ao_num
+        call get_ao_two_e_integrals_erf(q,r,s,ao_num,a1(1,q,r,s))
+      enddo
+    enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+
+  print *, '1st transformation'
+  ! 1st transformation
+  allocate (a2(ao_num,ao_num,ao_num,mo_num))
+  call dgemm('T','N', (ao_num*ao_num*ao_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*ao_num*ao_num))
+
+  ! 2nd transformation
+  print *, '2nd transformation'
+  deallocate (a1)
+  allocate (a1(ao_num,ao_num,mo_num,mo_num))
+  call dgemm('T','N', (ao_num*ao_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (ao_num*ao_num*mo_num))
+
+  ! 3rd transformation
+  print *, '3rd transformation'
+  deallocate (a2)
+  allocate (a2(ao_num,mo_num,mo_num,mo_num))
+  call dgemm('T','N', (ao_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*mo_num*mo_num))
+
+  ! 4th transformation
+  print *, '4th transformation'
+  deallocate (a1)
+  allocate (a1(mo_num,mo_num,mo_num,mo_num))
+  call dgemm('T','N', (mo_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (mo_num*mo_num*mo_num))
+
+  deallocate (a2)
+
+  integer :: n_integrals, size_buffer
+  integer(key_kind)  , allocatable :: buffer_i(:)
+  real(integral_kind), allocatable :: buffer_value(:)
+  size_buffer = min(ao_num*ao_num*ao_num,16000000)
+
+  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals)
+  allocate ( buffer_i(size_buffer), buffer_value(size_buffer) )
+
+  n_integrals = 0
+  !$OMP DO
+  do l=1,mo_num
+    do k=1,mo_num
+      do j=1,l
+        do i=1,k
+            if (abs(a1(i,j,k,l)) < mo_integrals_threshold) then
+              cycle
+            endif
+            n_integrals += 1
+            buffer_value(n_integrals) = a1(i,j,k,l)
+            !DIR$ FORCEINLINE
+            call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals))
+            if (n_integrals == size_buffer) then
+              call map_append(mo_integrals_erf_map, buffer_i, buffer_value, n_integrals)
+              n_integrals = 0
+            endif
+        enddo
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+
+  call map_append(mo_integrals_erf_map, buffer_i, buffer_value, n_integrals)
+
+  deallocate(buffer_i, buffer_value)
+  !$OMP END PARALLEL
+
+  deallocate (a1)
+
+  call map_sort(mo_integrals_erf_map)
+  call map_unique(mo_integrals_erf_map)
+
+end subroutine
+
+
 
  BEGIN_PROVIDER [ double precision, mo_two_e_int_erf_jj_from_ao, (mo_num,mo_num) ]
 &BEGIN_PROVIDER [ double precision, mo_two_e_int_erf_jj_exchange_from_ao, (mo_num,mo_num) ]

src/mu_of_r/EZFIO.cfg

@@ -6,7 +6,7 @@ size: (becke_numerical_grid.n_points_final_grid,determinants.n_states)
 
 [mu_of_r_potential]
 type: character*(32)
-doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated | pure_act]
+doc: type of potential for the mu(r) interaction: can be [ hf| cas_full | cas_truncated | pure_act]
 interface: ezfio, provider, ocaml
 default: hf
 

src/mu_of_r/mu_of_r_conditions.irp.f

@@ -26,7 +26,7 @@
   do ipoint = 1, n_points_final_grid
    if(mu_of_r_potential.EQ."hf")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_hf(ipoint)
-   else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
+   else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_psi_cas(ipoint,istate)
    else 
     print*,'you requested the following mu_of_r_potential'

src/mu_of_r/test_proj_op.irp.f

@@ -9,7 +9,7 @@ program projected_operators
   ! orbitals coming from core 
   no_core_density = .True.
   touch no_core_density
-  mu_of_r_potential = "cas_ful"
+  mu_of_r_potential = "cas_full"
   touch mu_of_r_potential 
   print*,'Using Valence Only functions'
 !  call test_f_HF_valence_ab

src/two_body_rdm/state_av_full_orb_2_rdm.irp.f

@@ -8,7 +8,7 @@
 !
 !                                     = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
+! WHERE ALL ORBITALS (i,j,k,l) BELONG TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
 !
@@ -149,7 +149,7 @@
 !
 !                                     = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,alpha} a_{l,alpha} a_{k,alpha} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
+! WHERE ALL ORBITALS (i,j,k,l) BELONG TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
 !
@@ -262,7 +262,7 @@
 !
 !                                     = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,beta} a^{\dagger}_{j,beta} a_{l,beta} a_{k,beta} |Psi_{istate}>
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
+! WHERE ALL ORBITALS (i,j,k,l) BELONG TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
 !
@@ -376,7 +376,7 @@
 !                                     = \sum_{istate} w(istate) * \sum_{sigma,sigma'} <Psi_{istate}| a^{\dagger}_{i,sigma} a^{\dagger'}_{j,sigma} a_{l,sigma'} a_{k,sigma} |Psi_{istate}>
 !
 !
-! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
+! WHERE ALL ORBITALS (i,j,k,l) BELONG TO ALL OCCUPIED ORBITALS : core, inactive and active
 !
 ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
 !
@@ -619,3 +619,4 @@
   !$OMP END PARALLEL
 
  END_PROVIDER
+

src/utils/constants.include.F

@@ -1,6 +1,6 @@
 integer, parameter :: max_dim = 511
 integer, parameter :: SIMD_vector = 32
-integer, parameter :: N_int_max = 32
+integer, parameter :: N_int_max = 128
 
 double precision, parameter :: pi =  dacos(-1.d0)
 double precision, parameter :: inv_pi =  1.d0/dacos(-1.d0)