Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev

bin/qp_convert_output_to_ezfio

@@ -335,7 +335,7 @@ def write_ezfio(res, filename):
 def get_full_path(file_path):
     file_path = os.path.expanduser(file_path)
     file_path = os.path.expandvars(file_path)
-    file_path = os.path.abspath(file_path)
+#    file_path = os.path.abspath(file_path)
     return file_path
 
 

src/ao_one_e_ints/pot_ao_ints.irp.f

@@ -3,6 +3,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
   !  Nucleus-electron interaction, in the |AO| basis set.
   !
   !  :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
+  !
+  !  These integrals also contain the pseudopotential integrals.
   END_DOC
   implicit none
   double precision               :: alpha, beta, gama, delta
@@ -75,12 +77,12 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
 
     !$OMP END DO
     !$OMP END PARALLEL
-  endif
-
     IF (DO_PSEUDO) THEN
        ao_integrals_n_e += ao_pseudo_integrals
     ENDIF
 
+  endif
+
 
   if (write_ao_integrals_n_e) then
     call ezfio_set_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)

src/ao_one_e_ints/screening.irp.f

@@ -3,14 +3,11 @@ logical function ao_one_e_integral_zero(i,k)
   integer, intent(in) :: i,k
 
   ao_one_e_integral_zero = .False.
-  if (.not.(read_ao_one_e_integrals.or.is_periodic)) then
+  if (.not.((io_ao_integrals_overlap/='None').or.is_periodic)) then
     if (ao_overlap_abs(i,k) < ao_integrals_threshold) then
         ao_one_e_integral_zero = .True.
         return
     endif
   endif
-  if (ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then
-      ao_one_e_integral_zero = .True.
-  endif
 end
 

src/ao_two_e_ints/screening.irp.f

@@ -8,8 +8,8 @@ logical function ao_two_e_integral_zero(i,j,k,l)
         ao_two_e_integral_zero = .True.
         return
     endif
-  endif
     if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k)  < ao_integrals_threshold) then
         ao_two_e_integral_zero = .True.
     endif
+  endif
 end

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -18,8 +18,7 @@ double precision function ao_two_e_integral(i,j,k,l)
 
    if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
      ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
-     return
+   else
-   endif
 
     dim1 = n_pt_max_integrals
 
@@ -101,6 +100,7 @@ double precision function ao_two_e_integral(i,j,k,l)
 
     endif
 
+  endif
 end
 
 double precision function ao_two_e_integral_schwartz_accel(i,j,k,l)
@@ -343,8 +343,6 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
   integer                        :: kk, m, j1, i1, lmax
   character*(64)                 :: fmt
 
-  integral = ao_two_e_integral(1,1,1,1)
-
   double precision               :: map_mb
   PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals
   if (read_ao_two_e_integrals) then
@@ -352,14 +350,18 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
     call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
     print*, 'AO integrals provided'
     ao_two_e_integrals_in_map = .True.
-    return
+  else
-  endif
 
     print*, 'Providing the AO integrals'
     call wall_time(wall_0)
     call wall_time(wall_1)
     call cpu_time(cpu_1)
 
+    if (.True.) then
+      ! Avoid openMP
+      integral = ao_two_e_integral(1,1,1,1)
+    endif
+
     integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
     call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals')
 
@@ -414,6 +416,8 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
       call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
     endif
 
+  endif
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num)  ]

src/becke_numerical_grid/EZFIO.cfg

@@ -14,3 +14,22 @@ type: double precision
 doc: threshold on the weight of a given grid point
 interface: ezfio,provider,ocaml
 default: 1.e-20
+
+[my_grid_becke]
+type: logical
+doc: if True, the number of angular and radial grid points are read from EZFIO
+interface: ezfio,provider,ocaml
+default: False
+
+[my_n_pt_r_grid]
+type: integer
+doc: Number of radial grid points given from input
+interface: ezfio,provider,ocaml
+default: 300
+
+[my_n_pt_a_grid]
+type: integer
+doc: Number of angular grid points given from input. Warning, this number cannot be any integer. See file list_angular_grid
+interface: ezfio,provider,ocaml
+default: 1202
+

src/becke_numerical_grid/grid_becke.irp.f

@@ -8,6 +8,7 @@
  !
  ! These numbers are automatically set by setting the grid_type_sgn parameter
  END_DOC
+if(.not.my_grid_becke)then
  select case (grid_type_sgn)
     case(0)
      n_points_radial_grid = 23
@@ -25,6 +26,10 @@ select case (grid_type_sgn)
       write(*,*) '!!! Quadrature grid not available !!!'
       stop
   end select
+else
+ n_points_radial_grid = my_n_pt_r_grid
+ n_points_integration_angular = my_n_pt_a_grid
+endif
 END_PROVIDER
 
 BEGIN_PROVIDER [integer, n_points_grid_per_atom]

src/becke_numerical_grid/list_angular_grid

@@ -0,0 +1,32 @@
+0006
+0014
+0026
+0038
+0050
+0074
+0086
+0110
+0146
+0170
+0194
+0230
+0266
+0302
+0350
+0434
+0590
+0770
+0974
+1202
+1454
+1730
+2030
+2354
+2702
+3074
+3470
+3890
+4334
+4802
+5294
+5810

src/cas_based_on_top/two_body_dens_rout.irp.f

@@ -0,0 +1,152 @@
+
+subroutine give_n2_ii_val_ab(r1,r2,two_bod_dens)
+ implicit none
+ BEGIN_DOC
+! contribution from purely inactive orbitals to n2_{\Psi^B}(r_1,r_2) for a CAS wave function 
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: two_bod_dens
+ integer :: i,j,m,n,i_m,i_n
+ integer :: i_i,i_j
+ double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ ! You get all orbitals in r_1 and r_2
+ allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
+ call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r2,mos_array_r2) 
+ ! You extract the inactive orbitals 
+ allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
+ do i_m = 1, n_inact_orb
+  mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m))
+ enddo
+ do i_m = 1, n_inact_orb
+  mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m))
+ enddo
+
+ two_bod_dens = 0.d0
+ ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
+ do m = 1, n_inact_orb ! electron 1 
+ ! You browse all OCCUPIED BETA  electrons in the \mathcal{A} space
+  do n = 1, n_inact_orb ! electron 2 
+   ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2)
+   two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) 
+  enddo
+ enddo
+end
+
+
+subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate)
+ BEGIN_DOC
+! contribution from inactive and active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function 
+ END_DOC
+ implicit none
+ integer, intent(in) :: istate
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: two_bod_dens
+ integer :: i,orb_i,a,orb_a,n,m,b
+ double precision :: rho
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
+ double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
+
+ two_bod_dens = 0.d0
+ ! You get all orbitals in r_1 and r_2
+ allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
+ call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r2,mos_array_r2) 
+
+ ! You extract the inactive orbitals 
+ allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
+ do i = 1, n_inact_orb
+  mos_array_inact_r1(i) = mos_array_r1(list_inact(i))
+ enddo
+ do i= 1, n_inact_orb
+  mos_array_inact_r2(i) = mos_array_r2(list_inact(i))
+ enddo
+
+ ! You extract the active orbitals 
+ allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) )
+ do i= 1, n_act_orb
+  mos_array_act_r1(i) = mos_array_r1(list_act(i))
+ enddo
+ do i= 1, n_act_orb
+  mos_array_act_r2(i) = mos_array_r2(list_act(i))
+ enddo
+
+ ! Contracted density : intermediate quantity
+ two_bod_dens = 0.d0
+ do a = 1, n_act_orb 
+  do i = 1, n_inact_orb
+   do b = 1, n_act_orb 
+    rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) 
+    two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho
+   enddo
+  enddo
+ enddo
+end
+
+
+subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate)
+ BEGIN_DOC
+! contribution from purely active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function 
+ END_DOC
+ implicit none
+ integer, intent(in) :: istate
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: two_bod_dens
+ integer :: i,orb_i,a,orb_a,n,m,b,c,d
+ double precision :: rho
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
+
+ two_bod_dens = 0.d0
+ ! You get all orbitals in r_1 and r_2
+ allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
+ call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r2,mos_array_r2) 
+
+ ! You extract the active orbitals 
+ allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) )
+ do i= 1, n_act_orb
+  mos_array_act_r1(i) = mos_array_r1(list_act(i))
+ enddo
+ do i= 1, n_act_orb
+  mos_array_act_r2(i) = mos_array_r2(list_act(i))
+ enddo
+
+ ! Contracted density : intermediate quantity
+ two_bod_dens = 0.d0
+ do a = 1, n_act_orb  ! 1
+  do b = 1, n_act_orb  ! 2
+   do c = 1, n_act_orb  ! 1
+    do d = 1, n_act_orb  ! 2
+     rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
+     two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) 
+    enddo
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine give_n2_cas(r1,r2,istate,n2_psi)
+ implicit none
+ BEGIN_DOC
+! returns mu(r), f_psi, n2_psi for a general cas wave function
+ END_DOC
+ integer, intent(in) :: istate
+ double precision, intent(in)  :: r1(3),r2(3)
+ double precision, intent(out) :: n2_psi
+ double precision :: two_bod_dens_ii
+ double precision :: two_bod_dens_ia
+ double precision :: two_bod_dens_aa
+ ! inactive-inactive part of n2_psi(r1,r2)
+ call give_n2_ii_val_ab(r1,r2,two_bod_dens_ii)
+ ! inactive-active part of n2_psi(r1,r2)
+ call give_n2_ia_val_ab(r1,r2,two_bod_dens_ia,istate)
+ ! active-active part of n2_psi(r1,r2)
+ call give_n2_aa_val_ab(r1,r2,two_bod_dens_aa,istate)
+
+ n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
+ 
+end

src/casscf/casscf.irp.f

@@ -5,6 +5,7 @@ program casscf
   END_DOC
   call reorder_orbitals_for_casscf
   no_vvvv_integrals = .True.
+  touch no_vvvv_integrals
   pt2_max = 0.02
   SOFT_TOUCH no_vvvv_integrals pt2_max
   call run_stochastic_cipsi

src/davidson/davidson_parallel.irp.f

@@ -438,6 +438,11 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
   ipos=1
   do imin=1,N_det,tasksize
     imax = min(N_det,imin-1+tasksize)
+    if (imin==1) then
+      istep = 2
+    else
+      istep = 1
+    endif
     do ishift=0,istep-1
       write(task(ipos:ipos+50),'(4(I11,1X),1X,1A)') imin, imax, ishift, istep, '|'
       ipos = ipos+50

src/ezfio_files/ezfio.irp.f

@@ -1,4 +1,4 @@
-BEGIN_PROVIDER [ character*(128), ezfio_filename ]
+BEGIN_PROVIDER [ character*(1024), ezfio_filename ]
   implicit none
   BEGIN_DOC
   ! Name of EZFIO file. It is obtained from the QPACKAGE_INPUT environment
@@ -34,7 +34,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ]
   ! Adjust out-of-memory killer flag such that the current process will be
   ! killed first by the OOM killer, allowing compute nodes to survive
   integer :: getpid
-  character*(64) :: command, pidc
+  character*(1024) :: command, pidc
   write(pidc,*) getpid()
   write(command,*) 'echo 15 > /proc//'//trim(adjustl(pidc))//'/oom_adj'
   call system(command)
@@ -43,7 +43,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ character*(128), ezfio_work_dir ]
+BEGIN_PROVIDER [ character*(1024), ezfio_work_dir ]
  implicit none
  BEGIN_DOC
  ! EZFIO/work/

src/ezfio_files/get_unit_and_open.irp.f

@@ -17,7 +17,7 @@ integer function getUnitAndOpen(f,mode)
   END_DOC
 
   character*(*)      :: f
-  character*(128)    :: new_f
+  character*(256)    :: new_f
   integer            :: iunit
   logical            :: is_open, exists
   character          :: mode

src/ezfio_files/qp_stop.irp.f

@@ -1,5 +1,5 @@
- BEGIN_PROVIDER [ character*(128), qp_stop_filename ]
+ BEGIN_PROVIDER [ character*(256), qp_stop_filename ]
-&BEGIN_PROVIDER [ character*(128), qp_kill_filename ]
+&BEGIN_PROVIDER [ character*(256), qp_kill_filename ]
 &BEGIN_PROVIDER [ integer, qp_stop_variable ]
  implicit none
  BEGIN_DOC

src/mo_two_e_ints/EZFIO.cfg

@@ -11,3 +11,9 @@ interface: ezfio,provider,ocaml
 default: 1.e-15
 ezfio_name: threshold_mo
 
+[no_vvvv_integrals]
+type: logical
+doc: If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
+interface: ezfio,provider,ocaml
+default: false
+

src/mo_two_e_ints/four_idx_novvvv.irp.f

@@ -1,11 +1,11 @@
-BEGIN_PROVIDER [ logical, no_vvvv_integrals  ]
+!BEGIN_PROVIDER [ logical, no_vvvv_integrals  ]
-  implicit none
+!  implicit none
-  BEGIN_DOC
+!  BEGIN_DOC
 ! If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
-  END_DOC
+!  END_DOC
-
+!
-  no_vvvv_integrals = .False.
+!  no_vvvv_integrals = .False.
-END_PROVIDER
+!END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ]
  implicit none
@@ -56,6 +56,8 @@ subroutine four_idx_novvvv
   BEGIN_DOC
   ! Retransform MO integrals for next CAS-SCF step
   END_DOC
+  print*,'Using partial transformation'
+  print*,'It will not transform all integrals with at least 3 indices within the virtuals'
   integer                        :: i,j,k,l,n_integrals
   double precision, allocatable  :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:)
   double precision, external     :: get_ao_two_e_integral

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -189,7 +189,6 @@ subroutine add_integrals_to_map(mask_ijkl)
       two_e_tmp_2 = 0.d0
       do j1 = 1,ao_num
         call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
-        ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
       enddo
       do j1 = 1,ao_num
         kmax = 0
@@ -747,7 +746,6 @@ subroutine add_integrals_to_map_no_exit_34(mask_ijkl)
       two_e_tmp_2 = 0.d0
       do j1 = 1,ao_num
         call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
-        ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1))
       enddo
       do j1 = 1,ao_num
         kmax = 0

src/nuclei/nuclei.irp.f

@@ -211,7 +211,7 @@ END_PROVIDER
    END_DOC
    integer                        :: iunit, i
    integer, external              :: getUnitAndOpen
-   character*(128)                :: filename
+   character*(1024)               :: filename
    if (mpi_master) then
      call getenv('QP_ROOT',filename)
      filename = trim(filename)//'/data/list_element.txt'

src/two_body_rdm/act_2_rdm.irp.f

@@ -2,6 +2,8 @@
  BEGIN_PROVIDER [double precision, act_2_rdm_ab_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
  implicit none
  BEGIN_DOC
+!                             12 12
+!                 1 2 1 2 == <ij|kl>
 ! act_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons 
 ! 
 ! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>

src/utils/integration.irp.f

@@ -75,7 +75,6 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   P_new(0,1) = 0.d0
   P_new(0,2) = 0.d0
   P_new(0,3) = 0.d0
-
   !DIR$ FORCEINLINE
   call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
   if (fact_k < thresh) then

src/utils/shank.irp.f

@@ -0,0 +1,54 @@
+double precision function shank_general(array,n,nmax)
+ implicit none                                                                                                                                                         
+ integer, intent(in) :: n,nmax
+ double precision, intent(in) :: array(0:nmax) ! array of the partial sums
+ integer :: ntmp,i
+ double precision :: sum(0:nmax),shank1(0:nmax)
+ if(n.lt.3)then
+  print*,'You asked to Shank a sum but the order is smaller than 3 ...'
+  print*,'n = ',n
+  print*,'stopping ....'
+  stop
+ endif
+ ntmp = n
+ sum = array
+ i = 0
+ do while(ntmp.ge.2)
+  i += 1
+!  print*,'i = ',i
+  call shank(sum,ntmp,nmax,shank1)
+  ntmp = ntmp - 2
+  sum = shank1
+  shank_general = shank1(ntmp)
+ enddo
+end
+
+
+subroutine shank(array,n,nmax,shank1)
+ implicit none
+ integer, intent(in) :: n,nmax
+ double precision, intent(in)  :: array(0:nmax)
+ double precision, intent(out) :: shank1(0:nmax)
+ integer :: i,j
+ double precision :: shank_function
+ do i = 1, n-1
+  shank1(i-1) = shank_function(array,i,nmax)
+ enddo
+end
+
+double precision function shank_function(array,i,n)
+ implicit none
+ integer, intent(in) :: i,n
+ double precision, intent(in) :: array(0:n)
+ double precision :: b_n, b_n1
+ b_n = array(i) - array(i-1)
+ b_n1 = array(i+1) - array(i)
+ if(dabs(b_n1-b_n).lt.1.d-12)then
+  shank_function = array(i+1)
+ else
+  shank_function = array(i+1) - b_n1*b_n1/(b_n1-b_n)
+ endif
+
+end
+
+

src/zmq/utils.irp.f

@@ -585,7 +585,7 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
     stop 'Wrong end of job'
   endif
 
-  do i=1200,1,-1
+  do i=360,1,-1
     rc = f77_zmq_send(zmq_to_qp_run_socket, 'end_job '//trim(zmq_state),8+len(trim(zmq_state)),0)
     rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 512, 0)
     if (trim(message(1:13)) == 'error waiting') then