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

src/ao_one_e_ints/pot_ao_ints.irp.f

@@ -104,6 +104,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
     IF(do_pseudo) THEN
        ao_integrals_n_e += ao_pseudo_integrals
     ENDIF
+    IF(point_charges) THEN
+       ao_integrals_n_e += ao_integrals_pt_chrg
+    ENDIF
 
   endif
 

src/ao_two_e_ints/cholesky.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
  ! Number of Cholesky vectors in AO basis
  END_DOC
 
- cholesky_ao_num_guess = ao_num*ao_num / 2
+ cholesky_ao_num_guess = ao_num*ao_num 
 END_PROVIDER
 
  BEGIN_PROVIDER [ integer, cholesky_ao_num ]
@@ -44,19 +44,12 @@ END_PROVIDER
    do m=0,9
      do l=1+m,ao_num,10
        !$OMP DO SCHEDULE(dynamic)
-       do j=1,l
+       do j=1,ao_num
          do k=1,ao_num
-           do i=1,min(k,j)
+           do i=1,ao_num
              if (ao_two_e_integral_zero(i,j,k,l)) cycle
              integral = get_ao_two_e_integral(i,j,k,l, ao_integrals_map)
              ao_integrals(i,k,j,l) = integral
-             ao_integrals(k,i,j,l) = integral
-             ao_integrals(i,k,l,j) = integral
-             ao_integrals(k,i,l,j) = integral
-             ao_integrals(j,l,i,k) = integral
-             ao_integrals(j,l,k,i) = integral
-             ao_integrals(l,j,i,k) = integral
-             ao_integrals(l,j,k,i) = integral
            enddo
          enddo
        enddo

src/ccsd/EZFIO.cfg

@@ -0,0 +1,11 @@
+[energy]
+type: double precision
+doc: CCSD energy
+interface: ezfio
+
+[energy_t]
+type: double precision
+doc: CCSD(T) energy
+interface: ezfio
+
+

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -135,8 +135,11 @@ subroutine run_ccsd_space_orb
   write(*,'(A15,1pE10.2,A3)')' Conv        = ', max_r
   print*,''
 
-  call write_t1(nO,nV,t1)
+  if (write_amplitudes) then
-  call write_t2(nO,nV,t2)
+    call write_t1(nO,nV,t1)
+    call write_t2(nO,nV,t2)
+    call ezfio_set_utils_cc_io_amplitudes('Read')
+  endif
 
   ! Deallocation
   if (cc_update_method == 'diis') then
@@ -147,6 +150,7 @@ subroutine run_ccsd_space_orb
 
   ! CCSD(T)
   double precision :: e_t
+  e_t = 0.d0
 
   if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then
 
@@ -182,8 +186,7 @@ subroutine run_ccsd_space_orb
     print*,''
   endif
 
-  print*,'Reference determinant:'
+  call save_energy(uncorr_energy + energy, e_t)
-  call print_det(det,N_int)
 
   deallocate(t1,t2)
 

src/ccsd/ccsd_spin_orb_sub.irp.f

@@ -269,8 +269,11 @@ subroutine run_ccsd_spin_orb
   write(*,'(A15,1pE10.2,A3)')' Conv        = ', max_r
   print*,''
 
-  call write_t1(nO,nV,t1)
+  if (write_amplitudes) then
-  call write_t2(nO,nV,t2)
+    call write_t1(nO,nV,t1)
+    call write_t2(nO,nV,t2)
+    call ezfio_set_utils_cc_io_amplitudes('Read')
+  endif
 
   ! Deallocate
   if (cc_update_method == 'diis') then
@@ -284,8 +287,9 @@ subroutine run_ccsd_spin_orb
   deallocate(v_ovoo,v_oovo)
   deallocate(v_ovvo,v_ovov,v_oovv)
   
+  double precision :: t_corr
+  t_corr = 0.d0
   if (cc_par_t .and. elec_alpha_num  +elec_beta_num > 2) then
-    double precision :: t_corr
     print*,'CCSD(T) calculation...'
     call wall_time(ta)
     !allocate(v_vvvo(nV,nV,nV,nO))
@@ -307,8 +311,8 @@ subroutine run_ccsd_spin_orb
     write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + t_corr, ' Ha'
     print*,''
   endif
-  print*,'Reference determinant:'
+
-  call print_det(det,N_int)
+  call save_energy(uncorr_energy + energy, t_corr)
   
   deallocate(f_oo,f_ov,f_vv,f_o,f_v)
   deallocate(v_ooov,v_vvoo,t1,t2)

src/ccsd/save_energy.irp.f

@@ -0,0 +1,13 @@
+subroutine save_energy(E,ET)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E, ET
+  call ezfio_set_ccsd_energy(E)
+  if (ET /= 0.d0) then
+    call ezfio_set_ccsd_energy_t(E+ET)
+  endif
+end
+
+

src/hartree_fock/10.hf.bats

@@ -43,12 +43,11 @@ python write_pt_charges.py ${EZFIO}
 qp set nuclei point_charges True
 qp run scf | tee ${EZFIO}.pt_charges.out
   energy="$(ezfio get hartree_fock energy)"
-good=-92.76613324421798
+good=-92.79920682236470
   eq $energy $good $thresh
 rm -rf $EZFIO
 }
 
-
 @test "H2_1" { # 1s
   run h2_1.ezfio -1.005924963288527
 }
@@ -85,6 +84,8 @@ rm -rf $EZFIO
   run hcn.ezfio -92.88717500035233
 }
 
+
+
 @test "B-B" { # 3s
   run b2_stretched.ezfio -48.9950585434279
 }

src/mo_optimization/my_providers.irp.f

@@ -1,141 +0,0 @@
-! Dimensions of MOs
-
-
-BEGIN_PROVIDER [ integer, n_mo_dim ]
-  implicit none
-  BEGIN_DOC
-  ! Number of different pairs (i,j) of MOs we can build,
-  ! with i>j
-  END_DOC
-
-  n_mo_dim = mo_num*(mo_num-1)/2
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, n_mo_dim_core ]
-  implicit none 
-  BEGIN_DOC
-  ! Number of different pairs (i,j) of core MOs we can build,
-  ! with i>j
-  END_DOC
-
-  n_mo_dim_core = dim_list_core_orb*(dim_list_core_orb-1)/2
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, n_mo_dim_act ]
-  implicit none
-  BEGIN_DOC
-  ! Number of different pairs (i,j) of active MOs we can build,
-  ! with i>j
-  END_DOC
-
-  n_mo_dim_act = dim_list_act_orb*(dim_list_act_orb-1)/2
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, n_mo_dim_inact ]
-  implicit none 
-  BEGIN_DOC
-  ! Number of different pairs (i,j) of inactive MOs we can build,
-  ! with i>j
-  END_DOC
-
-  n_mo_dim_inact = dim_list_inact_orb*(dim_list_inact_orb-1)/2
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, n_mo_dim_virt ]
-  implicit none 
-  BEGIN_DOC
-  ! Number of different pairs (i,j) of virtual MOs we can build,
-  ! with i>j
-  END_DOC
-
-  n_mo_dim_virt = dim_list_virt_orb*(dim_list_virt_orb-1)/2
-
-END_PROVIDER
-
-! Energies/criterions
-
-BEGIN_PROVIDER [ double precision, my_st_av_energy ]
-  implicit none
-  BEGIN_DOC
-  ! State average CI energy
-  END_DOC
-
-  !call update_st_av_ci_energy(my_st_av_energy)
-  call state_average_energy(my_st_av_energy)
-
-END_PROVIDER
-
-! With all the MOs
-
-BEGIN_PROVIDER [ double precision, my_gradient_opt, (n_mo_dim) ]
-&BEGIN_PROVIDER [ double precision, my_CC1_opt ]
-  implicit none
-  BEGIN_DOC
-  ! - Gradient of the energy with respect to the MO rotations, for all the MOs.
-  ! - Maximal element of the gradient in absolute value 
-  END_DOC
-
-  double precision :: norm_grad
-
-  PROVIDE mo_two_e_integrals_in_map
-
-  call gradient_opt(n_mo_dim, my_gradient_opt, my_CC1_opt, norm_grad)
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, my_hessian_opt, (n_mo_dim, n_mo_dim) ]
-  implicit none
-  BEGIN_DOC
-  ! - Gradient of the energy with respect to the MO rotations, for all the MOs.
-  ! - Maximal element of the gradient in absolute value 
-  END_DOC
-
-  double precision, allocatable :: h_f(:,:,:,:)
-
-  PROVIDE mo_two_e_integrals_in_map
-
-  allocate(h_f(mo_num, mo_num, mo_num, mo_num))
-
-  call hessian_list_opt(n_mo_dim, my_hessian_opt, h_f)
-
-END_PROVIDER
-
-! With the list of active MOs
-! Can be generalized to any mo_class by changing the list/dimension
-
-BEGIN_PROVIDER [ double precision, my_gradient_list_opt, (n_mo_dim_act) ]
-&BEGIN_PROVIDER [ double precision, my_CC2_opt ]
-  implicit none
-  BEGIN_DOC
-  ! - Gradient of the energy with respect to the MO rotations, only for the active MOs !
-  ! - Maximal element of the gradient in absolute value 
-  END_DOC
-
-  double precision :: norm_grad
-
-  PROVIDE mo_two_e_integrals_in_map !one_e_dm_mo two_e_dm_mo mo_one_e_integrals 
-
-  call gradient_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_gradient_list_opt, my_CC2_opt, norm_grad)
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, my_hessian_list_opt, (n_mo_dim_act, n_mo_dim_act) ]
-  implicit none
-  BEGIN_DOC
-  ! - Gradient of the energy with respect to the MO rotations, only for the active MOs !
-  ! - Maximal element of the gradient in absolute value 
-  END_DOC
-
-  double precision, allocatable :: h_f(:,:,:,:)
-
-  PROVIDE mo_two_e_integrals_in_map
-
-  allocate(h_f(dim_list_act_orb, dim_list_act_orb, dim_list_act_orb, dim_list_act_orb))
-
-  call hessian_list_opt(n_mo_dim_act, dim_list_act_orb, list_act, my_hessian_list_opt, h_f)
-
-END_PROVIDER

src/mo_two_e_ints/cholesky.irp.f

@@ -27,6 +27,8 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_ao_num, mo_num,
  double precision, allocatable :: buffer(:,:)
 
  print *, 'AO->MO Transformation of Cholesky vectors  .'
+
+ call set_multiple_levels_omp(.False.)
  !$OMP PARALLEL PRIVATE(i,j,k,buffer)
  allocate(buffer(mo_num,mo_num))
  !$OMP DO SCHEDULE(static)

src/nuclei/nuclei.irp.f

@@ -206,7 +206,12 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
      enddo
      nuclear_repulsion *= 0.5d0
      if(point_charges)then
-      nuclear_repulsion += pt_chrg_nuclei_interaction + pt_chrg_interaction
+      print*,'bear nuclear repulsion = ',nuclear_repulsion 
+      print*,'adding the interaction between the nuclein and the point charges'
+      print*,'to the usual nuclear repulsion '
+      nuclear_repulsion += pt_chrg_nuclei_interaction 
+      print*,'new nuclear repulsion =  ',nuclear_repulsion 
+      print*,'WARNING: we do not add the interaction between the point charges themselves'
      endif
    end if
 

src/nuclei/point_charges.irp.f

@@ -205,5 +205,8 @@ BEGIN_PROVIDER [ double precision, pt_chrg_nuclei_interaction]
  enddo
  print*,'Interaction between point charges and nuclei'
  print*,'pt_chrg_nuclei_interaction = ',pt_chrg_nuclei_interaction
+ if(point_charges)then
+  provide pt_chrg_interaction
+ endif
 END_PROVIDER 
 

src/utils/linear_algebra.irp.f

@@ -1831,7 +1831,7 @@ double precision, intent(in)     :: tol
 
 integer, dimension(:), allocatable          :: piv
 double precision, dimension(:), allocatable :: work
-character, parameter :: uplo = "U"
+character, parameter :: uplo = 'L'
 integer :: LDA
 integer :: info
 integer :: k, l, rank0
@@ -1848,14 +1848,14 @@ if (rank > rank0) then
 end if
 
 do k = 1, ndim
-  A(k+1:ndim, k) = 0.00D+0
+  A(k,k+1:ndim) = 0.00D+0
 end do
 ! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
 ! to do the swapping in-place
 U(:,:) = 0.00D+0
 do k = 1, ndim
   l = piv(k)
-  U(l, 1:rank) = A(1:rank, k)
+  U(l, 1:rank) = A(k,1:rank)
 end do
 
 end subroutine pivoted_cholesky

src/utils_cc/EZFIO.cfg

@@ -46,17 +46,11 @@ doc: Guess used to initialize the T2 amplitudes. none -> 0, MP -> perturbation t
 interface: ezfio,ocaml,provider
 default: MP
 
-[cc_write_t1]
+[io_amplitudes]
-type: logical
+type: Disk_access
-doc: If true, it will write on disk the T1 amplitudes at the end of the calculation.
+doc: Read/Write |CCSD| amplitudes from/to disk [ Write | Read | None ]
-interface: ezfio,ocaml,provider
+interface: ezfio,provider,ocaml
-default: False
+default: None
-
-[cc_write_t2]
-type: logical
-doc: If true, it will write on disk the T2 amplitudes at the end of the calculation.
-interface: ezfio,ocaml,provider
-default: False
 
 [cc_par_t]
 type: logical

src/utils_cc/guess_t.irp.f

@@ -91,16 +91,17 @@ subroutine write_t1(nO,nV,t1)
   double precision, intent(in) :: t1(nO, nV)
 
   ! internal
-  integer                      :: i,a
+  integer                      :: i,a, iunit
+  integer, external :: getunitandopen
 
-  if (cc_write_t1) then
+  if (write_amplitudes) then
-    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T1','w')
     do a = 1, nV
       do i = 1, nO
-         write(11,'(F20.12)') t1(i,a)
+         write(iunit,'(F20.12)') t1(i,a)
       enddo
     enddo
-    close(11)
+    close(iunit)
   endif
   
 end
@@ -120,20 +121,21 @@ subroutine write_t2(nO,nV,t2)
   double precision, intent(in) :: t2(nO, nO, nV, nV)
 
   ! internal
-  integer                      :: i,j,a,b
+  integer                      :: i,j,a,b, iunit
+  integer, external :: getunitandopen
 
-  if (cc_write_t2) then
+  if (write_amplitudes) then
-    open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T2','w')
     do b = 1, nV
       do a = 1, nV
         do j = 1, nO
           do i = 1, nO
-             write(11,'(F20.12)') t2(i,j,a,b)
+             write(iunit,'(F20.12)') t2(i,j,a,b)
           enddo
         enddo
       enddo
     enddo
-    close(11)
+    close(iunit)
   endif
   
 end
@@ -153,23 +155,19 @@ subroutine read_t1(nO,nV,t1)
   double precision, intent(out) :: t1(nO, nV)
 
   ! internal
-  integer                       :: i,a
+  integer                       :: i,a, iunit
   logical                       :: ok
+  integer, external :: getunitandopen
 
-  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
+  if (read_amplitudes) then
-  if (.not. ok) then
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T1','r')
-     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
+    do a = 1, nV
-     print*, 'Do a first calculation with cc_write_t1 = True'
+      do i = 1, nO
-     print*, 'and cc_guess_t1 /= read before setting cc_guess_t1 = read'
+         read(iunit,'(F20.12)') t1(i,a)
-     call abort
+      enddo
-  endif
-  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T1')
-  do a = 1, nV
-    do i = 1, nO
-       read(11,'(F20.12)') t1(i,a)
     enddo
-  enddo
+    close(iunit)
-  close(11)
+  endif
   
 end
 
@@ -188,26 +186,23 @@ subroutine read_t2(nO,nV,t2)
   double precision, intent(out) :: t2(nO, nO, nV, nV)
 
   ! internal
-  integer                       :: i,j,a,b
+  integer                       :: i,j,a,b, iunit
   logical                       :: ok
 
-  inquire(file=trim(ezfio_filename)//'/cc_utils/T1', exist=ok)
+  integer, external :: getunitandopen
-  if (.not. ok) then
+
-     print*, 'There is no file'// trim(ezfio_filename)//'/cc_utils/T1'
+  if (read_amplitudes) then
-     print*, 'Do a first calculation with cc_write_t2 = True'
+    iunit = getUnitAndOpen(trim(ezfio_filename)//'/work/T2','r')
-     print*, 'and cc_guess_t2 /= read before setting cc_guess_t2 = read'
+    do b = 1, nV
-     call abort
+      do a = 1, nV
-  endif
+        do j = 1, nO
-  open(unit=11, file=trim(ezfio_filename)//'/cc_utils/T2')
+          do i = 1, nO
-  do b = 1, nV
+             read(iunit,'(F20.12)') t2(i,j,a,b)
-    do a = 1, nV
+          enddo
-      do j = 1, nO
-        do i = 1, nO
-           read(11,'(F20.12)') t2(i,j,a,b)
         enddo
       enddo
     enddo
-  enddo
+    close(iunit)
-  close(11)
+  endif
   
 end

src/utils_cc/org/phase.org

@@ -137,6 +137,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N
   do j = 1, 2
     k = 1
     do i = 1, n1(j)
+       if (k > n_anni(j)) exit
        if (l1(i,j) /= list_anni(k,j)) cycle
        pos_anni(k,j) = i
        k = k + 1
@@ -147,6 +148,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N
   do j = 1, 2
     k = 1
     do i = 1, n2(j)
+       if (k > n_crea(j)) exit
        if (l2(i,j) /= list_crea(k,j)) cycle
        pos_crea(k,j) = i
        k = k + 1

src/utils_cc/phase.irp.f

@@ -96,6 +96,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N
   do j = 1, 2
     k = 1
     do i = 1, n1(j)
+       if (k > n_anni(j)) exit
        if (l1(i,j) /= list_anni(k,j)) cycle
        pos_anni(k,j) = i
        k = k + 1
@@ -106,6 +107,7 @@ subroutine get_excitation_general(det1,det2,degree,n,list_anni,list_crea,phase,N
   do j = 1, 2
     k = 1
     do i = 1, n2(j)
+       if (k > n_crea(j)) exit
        if (l2(i,j) /= list_crea(k,j)) cycle
        pos_crea(k,j) = i
        k = k + 1