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

bin/qp_reset

@@ -97,6 +97,8 @@ if [[ $dets -eq 1 ]] ; then
   rm --force -- ${ezfio}/determinants/psi_{det,coef}.gz
   rm --force -- ${ezfio}/determinants/n_det_qp_edit
   rm --force -- ${ezfio}/determinants/psi_{det,coef}_qp_edit.gz
+  rm --force -- ${ezfio}/tc_bi_ortho/psi_{l,r}_coef_bi_ortho.gz 
+
 fi
 
 if [[ $mos -eq 1 ]] ; then

configure

@@ -211,6 +211,7 @@ EOF
             execute << EOF
               cd "\${QP_ROOT}"/external
               wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
+              rm -rf trexio-${VERSION}
               tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
               cd trexio-${VERSION}
               ./configure --prefix=\${QP_ROOT} --without-hdf5 CFLAGS='-g'
@@ -224,6 +225,7 @@ EOF
             execute << EOF
               cd "\${QP_ROOT}"/external
               wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
+              rm -rf trexio-${VERSION}
               tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
               cd trexio-${VERSION}
               ./configure --prefix=\${QP_ROOT} CFLAGS="-g"
@@ -231,15 +233,28 @@ EOF
 EOF
     elif [[ ${PACKAGE} = qmckl ]] ; then
 
-            VERSION=0.5.3
+            VERSION=0.5.4
             execute << EOF
               cd "\${QP_ROOT}"/external
               wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz
+              rm -rf qmckl-${VERSION}
               tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
               cd qmckl-${VERSION}
               ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc CFLAGS='-g'
               make && make -j 4 check && make install
 EOF
+    elif [[ ${PACKAGE} = qmckl-intel ]] ; then
+
+            VERSION=0.5.4
+            execute << EOF
+              cd "\${QP_ROOT}"/external
+              wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz
+              rm -rf qmckl-${VERSION}
+              tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
+              cd qmckl-${VERSION}
+              ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc --with-icc --with-ifort CFLAGS='-g'
+              make && make -j 4 check && make install
+EOF
  
 
     elif [[ ${PACKAGE} = gmp ]] ; then
@@ -378,13 +393,13 @@ fi
 
 TREXIO=$(find_lib -ltrexio)
 if [[ ${TREXIO} = $(not_found) ]] ; then
-        error "TREXIO (trexio,trexio-nohdf5) is not installed. If you don't have HDF5, use trexio-nohdf5"
+        error "TREXIO (trexio | trexio-nohdf5) is not installed. If you don't have HDF5, use trexio-nohdf5"
         fail
 fi
 
 QMCKL=$(find_lib -lqmckl)
 if [[ ${QMCKL} = $(not_found) ]] ; then
-        error "QMCkl (qmckl) is not installed."
+        error "QMCkl (qmckl | qmckl-intel) is not installed."
         fail
 fi
 

scripts/module/module_handler.py

@@ -115,9 +115,7 @@ def get_l_module_descendant(d_child, l_module):
             except KeyError:
                 print("Error: ", file=sys.stderr)
                 print("`{0}` is not a submodule".format(module), file=sys.stderr)
-                print("Check the typo (spelling, case, '/', etc.) ", file=sys.stderr)
+                raise
-#                pass
-                sys.exit(1)
 
     return list(set(l))
 

src/ao_one_e_ints/kin_ao_ints.irp.f

@@ -52,7 +52,7 @@
     !$OMP DEFAULT(NONE) &
     !$OMP PRIVATE(A_center,B_center,power_A,power_B,&
     !$OMP  overlap_y, overlap_z, overlap, &
-    !$OMP  alpha, beta,i,j,c,d_a_2,d_2,deriv_tmp, &
+    !$OMP  alpha, beta, n, l, i,j,c,d_a_2,d_2,deriv_tmp, &
     !$OMP  overlap_x0,overlap_y0,overlap_z0) &
     !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
     !$OMP  ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &

src/ao_tc_eff_map/NEED

@@ -1,4 +1,4 @@
-ao_two_e_erf_ints
+ao_two_e_ints
 mo_one_e_ints
 ao_many_one_e_ints
 dft_utils_in_r

src/ao_two_e_erf_ints/EZFIO.cfg

@@ -1,13 +0,0 @@
-[io_ao_two_e_integrals_erf]
-type: Disk_access
-doc: Read/Write |AO| integrals with the long range interaction from/to disk [ Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: None
-
-[mu_erf]
-type: double precision
-doc: cutting of the interaction in the range separated model
-interface: ezfio,provider,ocaml
-default: 0.5
-ezfio_name: mu_erf
-

src/ao_two_e_erf_ints/NEED

@@ -1 +0,0 @@
-ao_two_e_ints

src/ao_two_e_erf_ints/README.rst

@@ -1,19 +0,0 @@
-======================
-ao_two_e_erf_ints
-======================
-
-Here, all two-electron integrals (:math:`erf(\mu r_{12})/r_{12}`) are computed.
-As they have 4 indices and many are zero, they are stored in a map, as defined
-in :file:`utils/map_module.f90`.
-
-The main parameter of this module is :option:`ao_two_e_erf_ints mu_erf` which is the range-separation parameter.
-
-To fetch an |AO| integral, use the
-`get_ao_two_e_integral_erf(i,j,k,l,ao_integrals_erf_map)` function.
-
-
-The conventions are:
-* For |AO| integrals : (ij|kl) = (11|22) = <ik|jl> = <12|12>
-
-
-

src/ao_two_e_ints/EZFIO.cfg

@@ -35,3 +35,15 @@ type: logical
 doc: Perform Cholesky decomposition of AO integrals
 interface: ezfio,provider,ocaml
 default: False
+
+[io_ao_two_e_integrals_erf]
+type: Disk_access
+doc: Read/Write |AO| erf integrals from/to disk [ Write | Read | None ] 
+interface: ezfio,provider,ocaml
+default: None
+
+[use_only_lr]
+type: logical
+doc: If true, use only the long range part of the two-electron integrals instead of 1/r12
+interface: ezfio, provider, ocaml
+default: False

src/ao_two_e_ints/NEED

@@ -1,3 +1,4 @@
+hamiltonian
 ao_one_e_ints
 pseudo
 bitmask

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -21,9 +21,9 @@ double precision function ao_two_e_integral(i, j, k, l)
   double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
   double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
 
-  double precision               :: ao_two_e_integral_schwartz_accel
+  double precision, external     :: ao_two_e_integral_erf
-
+  double precision, external     :: ao_two_e_integral_cosgtos
-  double precision               :: ao_two_e_integral_cosgtos
+  double precision, external     :: ao_two_e_integral_schwartz_accel
 
 
   if(use_cosgtos) then
@@ -31,9 +31,11 @@ double precision function ao_two_e_integral(i, j, k, l)
 
     ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l)
 
-  else
+  else if (use_only_lr) then
 
-    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_erf(i, j, k, l)
+
+  else 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)
 
@@ -121,8 +123,6 @@ double precision function ao_two_e_integral(i, j, k, l)
 
   endif
 
-  endif
-
 end
 
 ! ---

src/casscf_cipsi/EZFIO.cfg

@@ -73,3 +73,9 @@ type: logical
 doc: If |true|, the pt2_max value in the CIPSI iterations will automatically adapt, otherwise it is fixed at the value given in the EZFIO folder
 interface: ezfio,provider,ocaml
 default: True
+
+[small_active_space]
+type: logical
+doc: If |true|, the pt2_max value in the CIPSI is set to 10-10 and will not change
+interface: ezfio,provider,ocaml
+default: False

src/casscf_cipsi/README.rst

@@ -3,3 +3,45 @@ casscf
 ======
 
 |CASSCF| program with the CIPSI algorithm.
+
+Example of inputs
+-----------------
+
+a) Small active space : standard CASSCF 
+---------------------------------------
+Let's do O2 (triplet) in aug-cc-pvdz with the following geometry (xyz format, Bohr units)
+3
+
+ O           0.0000000000        0.0000000000       -1.1408000000
+ O           0.0000000000        0.0000000000        1.1408000000
+
+# Create the ezfio folder 
+qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz
+
+# Start with an ROHF guess 
+qp run scf | tee ${EZFIO_FILE}.rohf.out
+
+# Get the ROHF energy for check 
+qp get hartree_fock energy # should be -149.4684509
+
+# Define the full valence active space: the two 1s are doubly occupied, the other 8 valence orbitals are active 
+# CASSCF(12e,10orb) 
+qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
+
+# Specify that you want an near exact CASSCF, i.e. the CIPSI selection will stop at pt2_max = 10^-10
+qp set casscf_cipsi small_active_space True 
+# RUN THE CASSCF 
+qp run casscf | tee ${EZFIO_FILE}.casscf.out
+# you should find around -149.7243542
+
+
+b) Large active space : Exploit the selected CI in the active space 
+-------------------------------------------------------------------
+#Let us start from the small active space calculation orbitals and add another 10 virtuals: CASSCF(12e,20orb)
+qp set_mo_class -c "[1-2]" -a "[3-20]" -v "[21-46]"
+# As this active space is larger, you unset the small_active_space feature 
+qp set casscf_cipsi small_active_space False
+# As it is a large active space, the energy convergence thereshold is set to be 0.0001
+qp run casscf | tee ${EZFIO_FILE}.casscf_large.out
+# you should find around -149.9046
+

src/casscf_cipsi/casscf.irp.f

@@ -8,17 +8,23 @@ program casscf
 !  touch no_vvvv_integrals
   n_det_max_full = 500
   touch n_det_max_full
+  if(small_active_space)then
+   pt2_relative_error = 0.00001
+  else
+   thresh_scf = 1.d-4
    pt2_relative_error = 0.04
+  endif
   touch pt2_relative_error 
-!  call run_stochastic_cipsi
   call run
 end
 
 subroutine run
   implicit none
-  double precision               :: energy_old, energy, pt2_max_before, ept2_before,delta_E
+  double precision               :: energy_old, energy, pt2_max_before,delta_E
   logical                        :: converged,state_following_casscf_cipsi_save
-  integer                        :: iteration
+  integer                        :: iteration,istate
+  double precision, allocatable :: E_PT2(:), PT2(:), Ev(:), ept2_before(:)
+  allocate(E_PT2(N_states), PT2(N_states), Ev(N_states), ept2_before(N_states))
   converged = .False.
 
   energy = 0.d0
@@ -28,13 +34,20 @@ subroutine run
   state_following_casscf = .True.
   touch state_following_casscf
   ept2_before = 0.d0
+  if(small_active_space)then
+   pt2_max = 1.d-10
+   SOFT_TOUCH pt2_max
+  else
    if(adaptive_pt2_max)then
     pt2_max = 0.005
     SOFT_TOUCH pt2_max
    endif
+  endif
   do while (.not.converged)
     print*,'pt2_max = ',pt2_max
-    call run_stochastic_cipsi
+    call run_stochastic_cipsi(Ev,PT2)
+    print*,'Ev,PT2',Ev(1),PT2(1)
+    E_PT2(1:N_states) = Ev(1:N_states) + PT2(1:N_states)
     energy_old = energy
     energy = eone+etwo+ecore
     pt2_max_before = pt2_max
@@ -42,15 +55,13 @@ subroutine run
     call write_time(6)
     call write_int(6,iteration,'CAS-SCF iteration = ')
     call write_double(6,energy,'CAS-SCF energy = ')
-    if(n_states == 1)then
+!    if(n_states == 1)then
-     double precision :: E_PT2, PT2
+!     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
-     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
+!     call ezfio_get_casscf_cipsi_energy(PT2)
-     call ezfio_get_casscf_cipsi_energy(PT2)
+     call write_double(6,E_PT2(1:N_states),'E + PT2 energy = ')
-     PT2 -= E_PT2
+     call write_double(6,PT2(1:N_states),'  PT2          = ')
-     call write_double(6,E_PT2,'E + PT2 energy = ')
-     call write_double(6,PT2,'  PT2          = ')
      call write_double(6,pt2_max,' PT2_MAX       = ')
-    endif
+!    endif
 
     print*,''
     call write_double(6,norm_grad_vec2,'Norm of gradients = ')
@@ -65,10 +76,14 @@ subroutine run
     else if (criterion_casscf == "gradients")then
      converged = norm_grad_vec2 < thresh_scf
     else if (criterion_casscf == "e_pt2")then
-     delta_E = dabs(E_PT2 - ept2_before)
+     delta_E = 0.d0
+     do istate = 1, N_states
+      delta_E += dabs(E_PT2(istate) - ept2_before(istate))
+     enddo
      converged = dabs(delta_E) < thresh_casscf
     endif
     ept2_before = E_PT2
+    if(.not.small_active_space)then
      if(adaptive_pt2_max)then
       pt2_max = dabs(energy_improvement / (pt2_relative_error))
       pt2_max = min(pt2_max, pt2_max_before)
@@ -76,6 +91,7 @@ subroutine run
        pt2_max = max(pt2_max,pt2_min_casscf)
       endif
      endif
+    endif
     print*,''
     call write_double(6,pt2_max, 'PT2_MAX for next iteration = ')
 
@@ -94,9 +110,11 @@ subroutine run
      read_wf = .True.
      call clear_mo_map
      SOFT_TOUCH mo_coef N_det psi_det psi_coef
+     if(.not.small_active_space)then
       if(adaptive_pt2_max)then
         SOFT_TOUCH pt2_max  
       endif
+     endif
      if(iteration .gt. 3)then
       state_following_casscf = state_following_casscf_cipsi_save
       soft_touch state_following_casscf
@@ -104,6 +122,25 @@ subroutine run
     endif
 
   enddo
+     integer :: i
+    print*,'Converged CASSCF '
+    print*,'--------------------------'
+    write(6,*) ' occupation numbers of orbitals '
+    do i=1,mo_num
+      write(6,*) i,occnum(i)
+    end do
+    print*,'--------------'
+!
+!     write(6,*)
+!     write(6,*) ' the diagonal of the inactive effective Fock matrix '
+!     write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
+!     write(6,*)
+  print*,'Fock MCSCF'
+  do i = 1, mo_num
+   write(*,*)i,mcscf_fock_diag_mo(i)
+!   write(*,*)mcscf_fock_alpha_mo(i,i)
+  enddo
+
 
 end
 

src/casscf_cipsi/densities.irp.f

@@ -15,6 +15,35 @@ BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
     enddo
   enddo
   
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, D0tu_alpha_ao, (ao_num, ao_num)]
+&BEGIN_PROVIDER [double precision, D0tu_beta_ao, (ao_num, ao_num)]
+ implicit none
+ integer :: i,ii,j,u,t,uu,tt
+ double precision, allocatable :: D0_tmp_alpha(:,:),D0_tmp_beta(:,:)
+ allocate(D0_tmp_alpha(mo_num, mo_num),D0_tmp_beta(mo_num, mo_num))
+ D0_tmp_beta = 0.d0
+ D0_tmp_alpha = 0.d0
+ do i = 1, n_core_inact_orb
+  ii = list_core_inact(i)
+  D0_tmp_alpha(ii,ii) = 1.d0
+  D0_tmp_beta(ii,ii) = 1.d0
+ enddo
+ print*,'Diagonal elements of the 1RDM in the active space'
+ do u=1,n_act_orb
+   uu = list_act(u)
+   print*,uu,one_e_dm_mo_alpha_average(uu,uu),one_e_dm_mo_beta_average(uu,uu)
+   do t=1,n_act_orb
+    tt = list_act(t)
+    D0_tmp_alpha(tt,uu) = one_e_dm_mo_alpha_average(tt,uu)
+    D0_tmp_beta(tt,uu) = one_e_dm_mo_beta_average(tt,uu)
+   enddo
+ enddo
+
+ call mo_to_ao_no_overlap(D0_tmp_alpha,mo_num,D0tu_alpha_ao,ao_num)
+ call mo_to_ao_no_overlap(D0_tmp_beta,mo_num,D0tu_beta_ao,ao_num)
+
 END_PROVIDER 
 
 BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]

src/casscf_cipsi/mcscf_fock.irp.f

@@ -77,4 +77,119 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
    
 END_PROVIDER
  
+ BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_ao, (ao_num, ao_num)] 
+&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_ao, (ao_num, ao_num)] 
+ implicit none
+ BEGIN_DOC
+  ! mcscf_fock_alpha_ao are set to usual Fock like operator but computed with the MCSCF densities on the AO basis 
+ END_DOC
+ SCF_density_matrix_ao_alpha = D0tu_alpha_ao
+ SCF_density_matrix_ao_beta = D0tu_beta_ao
+ soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta 
+ mcscf_fock_beta_ao = fock_matrix_ao_beta
+ mcscf_fock_alpha_ao = fock_matrix_ao_alpha
+END_PROVIDER 
 
+
+ BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_mo, (mo_num, mo_num)] 
+&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_mo, (mo_num, mo_num)] 
+ implicit none
+ BEGIN_DOC
+  ! Mo_mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities on the MO basis 
+ END_DOC
+
+ call ao_to_mo(mcscf_fock_alpha_ao,ao_num,mcscf_fock_alpha_mo,mo_num)
+ call ao_to_mo(mcscf_fock_beta_ao,ao_num,mcscf_fock_beta_mo,mo_num)
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision, mcscf_fock_mo, (mo_num,mo_num) ]
+&BEGIN_PROVIDER [ double precision, mcscf_fock_diag_mo, (mo_num)]
+   implicit none
+   BEGIN_DOC
+   ! MCSF Fock matrix on the MO basis.
+   ! For open shells, the ROHF Fock Matrix is ::
+   !
+   !       |  Rcc  |  F^b  |  Fcv  |
+   !       |-----------------------|
+   !       |  F^b  |  Roo  |  F^a  |
+   !       |-----------------------|
+   !       |  Fcv  |  F^a  |  Rvv  |
+   !
+   ! C: Core, O: Open, V: Virtual 
+   ! 
+   ! Rcc = Acc Fcc^a + Bcc Fcc^b
+   ! Roo = Aoo Foo^a + Boo Foo^b
+   ! Rvv = Avv Fvv^a + Bvv Fvv^b
+   ! Fcv = (F^a + F^b)/2
+   ! 
+   ! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO)
+   ! A,B: Coupling parameters
+   !
+   ! J. Chem. Phys. 133, 141102 (2010), https://doi.org/10.1063/1.3503173
+   ! Coupling parameters from J. Chem. Phys. 125, 204110 (2006); https://doi.org/10.1063/1.2393223.
+   !       cc   oo   vv
+   !  A  -0.5  0.5  1.5
+   !  B   1.5  0.5 -0.5
+   ! 
+   END_DOC
+   integer                        :: i,j,n
+   if (elec_alpha_num == elec_beta_num) then
+     mcscf_fock_mo = mcscf_fock_alpha_mo
+   else
+     ! Core
+     do j = 1, elec_beta_num
+       ! Core
+       do i = 1, elec_beta_num
+         mcscf_fock_mo(i,j) = - 0.5d0 * mcscf_fock_alpha_mo(i,j) &
+                               + 1.5d0 * mcscf_fock_beta_mo(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         mcscf_fock_mo(i,j) = mcscf_fock_beta_mo(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         mcscf_fock_mo(i,j) =   0.5d0 * mcscf_fock_alpha_mo(i,j) &
+                               + 0.5d0 * mcscf_fock_beta_mo(i,j)
+       enddo
+     enddo
+     ! Open
+     do j = elec_beta_num+1, elec_alpha_num
+       ! Core
+       do i = 1, elec_beta_num
+         mcscf_fock_mo(i,j) = mcscf_fock_beta_mo(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         mcscf_fock_mo(i,j) =   0.5d0 * mcscf_fock_alpha_mo(i,j) &
+                               + 0.5d0 * mcscf_fock_beta_mo(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         mcscf_fock_mo(i,j) = mcscf_fock_alpha_mo(i,j)
+       enddo
+     enddo
+     ! Virtual
+     do j = elec_alpha_num+1, mo_num
+       ! Core
+       do i = 1, elec_beta_num
+         mcscf_fock_mo(i,j) =   0.5d0 * mcscf_fock_alpha_mo(i,j) &
+                               + 0.5d0 * mcscf_fock_beta_mo(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         mcscf_fock_mo(i,j) = mcscf_fock_alpha_mo(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         mcscf_fock_mo(i,j) =   1.5d0 * mcscf_fock_alpha_mo(i,j) &
+                               - 0.5d0 * mcscf_fock_beta_mo(i,j)
+       enddo
+     enddo
+   endif
+
+ do i = 1, mo_num
+  mcscf_fock_diag_mo(i) = mcscf_fock_mo(i,i)
+ enddo
+END_PROVIDER 

src/cipsi/stochastic_cipsi.irp.f

@@ -1,10 +1,11 @@
-subroutine run_stochastic_cipsi
+subroutine run_stochastic_cipsi(Ev,PT2) 
   use selection_types
   implicit none
   BEGIN_DOC
 ! Selected Full Configuration Interaction with Stochastic selection and PT2.
   END_DOC
   integer                        :: i,j,k
+  double precision, intent(out)  :: Ev(N_states), PT2(N_states) 
   double precision, allocatable  :: zeros(:)
   integer                        :: to_select
   type(pt2_type)                 :: pt2_data, pt2_data_err
@@ -79,12 +80,14 @@ subroutine run_stochastic_cipsi
     to_select = max(N_states_diag, to_select)
 
 
+    Ev(1:N_states) = psi_energy_with_nucl_rep(1:N_states)
     call pt2_dealloc(pt2_data)
     call pt2_dealloc(pt2_data_err)
     call pt2_alloc(pt2_data, N_states)
     call pt2_alloc(pt2_data_err, N_states)
     call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error,to_select) ! Stochastic PT2 and selection
 
+    PT2(1:N_states) = pt2_data % pt2(1:N_states)
     correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref)  /     &
                     (psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref)
     correlation_energy_ratio = min(1.d0,correlation_energy_ratio)

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -286,7 +286,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
       ! Small
       h(N_st_diag*itermax,N_st_diag*itermax),                        &
-      h_p(N_st_diag*itermax,N_st_diag*itermax),                      &
+!      h_p(N_st_diag*itermax,N_st_diag*itermax),                      &
       y(N_st_diag*itermax,N_st_diag*itermax),                        &
       s_(N_st_diag*itermax,N_st_diag*itermax),                       &
       s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
@@ -340,7 +340,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       exit
     endif
 
-    do iter=1,itermax-1
+    iter = 0
+    do while (iter < itermax-1)
+      iter += 1
+!    do iter=1,itermax-1
 
       shift  = N_st_diag*(iter-1)
       shift2 = N_st_diag*iter
@@ -430,30 +433,30 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
       call dgemm('T','N', shift2, shift2, sze,                       &
           1.d0, U, size(U,1), W, size(W,1),                          &
-          0.d0, h, size(h_p,1))
+          0.d0, h, size(h,1))
       call dgemm('T','N', shift2, shift2, sze,                       &
           1.d0, U, size(U,1), U, size(U,1),                          &
           0.d0, s_tmp, size(s_tmp,1))
 
-      ! Penalty method
+!      ! Penalty method
-      ! --------------
+!      ! --------------
-
+!
-      if (s2_eig) then
+!      if (s2_eig) then
-        h_p = s_
+!        h_p = s_
-        do k=1,shift2
+!        do k=1,shift2
-          h_p(k,k) = h_p(k,k) - expected_s2
+!          h_p(k,k) = h_p(k,k) - expected_s2
-        enddo
+!        enddo
-        if (only_expected_s2) then
+!        if (only_expected_s2) then
-          alpha = 0.1d0
+!          alpha = 0.1d0
-          h_p = h + alpha*h_p
+!          h_p = h + alpha*h_p
-        else
+!        else
-          alpha = 0.0001d0
+!          alpha = 0.0001d0
-          h_p = h + alpha*h_p
+!          h_p = h + alpha*h_p
-        endif
+!        endif
-      else
+!      else
-        h_p = h
+!        h_p = h
-        alpha = 0.d0
+!        alpha = 0.d0
-      endif
+!      endif
 
       ! Diagonalize h_p
       ! ---------------
@@ -473,8 +476,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        call dsygv(1,'V','U',shift2,y,size(y,1), &
           s_tmp,size(s_tmp,1), lambda, work,lwork,info)
        deallocate(work)
-       if (info /= 0) then
+       if (info > 0) then
-         stop 'DSYGV Diagonalization failed'
+         ! Numerical errors propagate. We need to reduce the number of iterations
+         itermax = iter-1
+         exit
        endif
 
       ! Compute Energy for each eigenvector

src/determinants/density_matrix.irp.f

@@ -493,3 +493,101 @@ subroutine get_occupation_from_dets(istate,occupation)
   enddo
 end
 
+BEGIN_PROVIDER [double precision, difference_dm, (mo_num, mo_num, N_states)]
+ implicit none
+ BEGIN_DOC
+! difference_dm(i,j,istate) = dm(i,j,1) - dm(i,j,istate)
+ END_DOC
+ integer :: istate
+ do istate = 1, N_states
+   difference_dm(:,:,istate) =  one_e_dm_mo_alpha(:,:,1) + one_e_dm_mo_beta(:,:,1) & 
+                             - (one_e_dm_mo_alpha(:,:,istate) + one_e_dm_mo_beta(:,:,istate))
+ enddo
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, difference_dm_eigvect, (mo_num, mo_num, N_states) ]
+&BEGIN_PROVIDER [double precision, difference_dm_eigval, (mo_num,  N_states) ]
+ implicit none
+ BEGIN_DOC
+! eigenvalues and eigevenctors of the difference_dm
+ END_DOC
+ integer :: istate,i
+ do istate = 2, N_states
+  call lapack_diag(difference_dm_eigval(1,istate),difference_dm_eigvect(1,1,istate)& 
+                  ,difference_dm(1,1,istate),mo_num,mo_num)
+  print*,'Eigenvalues of difference_dm for state ',istate
+  do i = 1, mo_num
+   print*,i,difference_dm_eigval(i,istate)
+  enddo
+ enddo
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ integer         , n_attachment,  (N_states)]
+&BEGIN_PROVIDER [ integer         , n_dettachment, (N_states)]
+&BEGIN_PROVIDER [ integer         , list_attachment,  (mo_num,N_states)]
+&BEGIN_PROVIDER [ integer         , list_dettachment, (mo_num,N_states)]
+ implicit none
+ integer :: i,istate
+ integer :: list_attachment_tmp(mo_num)
+ n_attachment = 0
+ n_dettachment = 0
+ do istate = 2, N_states
+  do i = 1, mo_num
+   if(difference_dm_eigval(i,istate).lt.0.d0)then ! dettachment_orbitals
+    n_dettachment(istate) += 1
+    list_dettachment(n_dettachment(istate),istate) = i ! they are already sorted 
+   else
+    n_attachment(istate) += 1
+    list_attachment_tmp(n_attachment(istate)) = i ! they are not sorted 
+   endif
+  enddo
+  ! sorting the attachment
+  do i = 0, n_attachment(istate) - 1
+   list_attachment(i+1,istate) = list_attachment_tmp(n_attachment(istate) - i)
+  enddo
+ enddo
+ 
+END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision,  attachment_numbers_sorted,  (mo_num, N_states)]
+&BEGIN_PROVIDER [ double precision,  dettachment_numbers_sorted, (mo_num, N_states)]
+ implicit none
+ integer :: i,istate
+ do istate = 2, N_states
+  print*,'dettachment'
+  do i = 1, n_dettachment(istate)
+   dettachment_numbers_sorted(i,istate) = difference_dm_eigval(list_dettachment(i,istate),istate)
+   print*,i,list_dettachment(i,istate),dettachment_numbers_sorted(i,istate)
+  enddo
+  print*,'attachment'
+  do i = 1, n_attachment(istate)
+   attachment_numbers_sorted(i,istate) = difference_dm_eigval(list_attachment(i,istate),istate)
+   print*,i,list_attachment(i,istate),attachment_numbers_sorted(i,istate)
+  enddo
+ enddo
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [ double precision,  attachment_orbitals, (ao_num, mo_num, N_states)]
+&BEGIN_PROVIDER [ double precision, dettachment_orbitals, (ao_num, mo_num, N_states)]
+ implicit none
+ integer :: i,j,k,istate
+ attachment_orbitals = 0.d0
+ dettachment_orbitals = 0.d0
+ do istate = 2, N_states
+   do i = 1, n_dettachment(istate)
+    do j = 1, mo_num
+     do k = 1, ao_num
+      dettachment_orbitals(k,list_dettachment(i,istate),istate) += mo_coef(k,j) * difference_dm_eigvect(j,list_dettachment(i,istate),istate) 
+     enddo
+    enddo
+   enddo
+   do i = 1, n_attachment(istate)
+    do j = 1, mo_num
+     do k = 1, ao_num
+      attachment_orbitals(k,i,istate) += mo_coef(k,j) * difference_dm_eigvect(j,list_attachment(i,istate),istate) 
+     enddo
+    enddo
+   enddo
+ enddo
+
+END_PROVIDER 

src/determinants/dipole_moments.irp.f

@@ -26,10 +26,10 @@
   enddo
  enddo
  
-! print*,'electron part for z_dipole = ',z_dipole_moment
+ print*,'electron part for z_dipole = ',z_dipole_moment
-! print*,'electron part for y_dipole = ',y_dipole_moment
+ print*,'electron part for y_dipole = ',y_dipole_moment
-! print*,'electron part for x_dipole = ',x_dipole_moment
+ print*,'electron part for x_dipole = ',x_dipole_moment
-! 
+ 
  nuclei_part_z = 0.d0
  nuclei_part_y = 0.d0
  nuclei_part_x = 0.d0
@@ -38,10 +38,10 @@
   nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) 
   nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) 
  enddo
-! print*,'nuclei   part for z_dipole = ',nuclei_part_z
+ print*,'nuclei   part for z_dipole = ',nuclei_part_z
-! print*,'nuclei   part for y_dipole = ',nuclei_part_y
+ print*,'nuclei   part for y_dipole = ',nuclei_part_y
-! print*,'nuclei   part for x_dipole = ',nuclei_part_x
+ print*,'nuclei   part for x_dipole = ',nuclei_part_x
-!
+
  do istate = 1, N_states
   z_dipole_moment(istate) += nuclei_part_z
   y_dipole_moment(istate) += nuclei_part_y

src/dft_one_e/NEED

@@ -4,6 +4,4 @@ mo_one_e_ints
 mo_two_e_ints
 ao_one_e_ints
 ao_two_e_ints
-mo_two_e_erf_ints
-ao_two_e_erf_ints
 mu_of_r

src/dummy/NEED

@@ -1,6 +1,5 @@
 ao_basis
 ao_one_e_ints
-ao_two_e_erf_ints
 ao_two_e_ints
 aux_quantities
 becke_numerical_grid
@@ -24,13 +23,13 @@ functionals
 generators_cas
 generators_full
 hartree_fock
+hamiltonian
 iterations
 kohn_sham
 kohn_sham_rs
 mo_basis
 mo_guess
 mo_one_e_ints
-mo_two_e_erf_ints
 mo_two_e_ints
 mpi
 nuclei

src/fci/fci.irp.f

@@ -41,8 +41,10 @@ program fci
 
     write(json_unit,json_array_open_fmt) 'fci'
 
+    double precision, allocatable :: Ev(:),PT2(:)
+    allocate(Ev(N_states), PT2(N_states))
     if (do_pt2) then
-      call run_stochastic_cipsi
+      call run_stochastic_cipsi(Ev,PT2)
     else
       call run_cipsi
     endif

src/hamiltonian/EZFIO.cfg

@@ -0,0 +1,8 @@
+[mu_erf]
+type: double precision
+doc: cutting of the interaction in the range separated model
+interface: ezfio,provider,ocaml
+default: 0.5
+ezfio_name: mu_erf
+
+

src/hamiltonian/README.rst

@@ -0,0 +1,5 @@
+===========
+hamiltonian
+===========
+
+Parameters of the Hamiltonian.

src/iterations/print_extrapolation.irp.f

@@ -37,7 +37,7 @@ subroutine print_extrapolated_energy
     write(*,*)  'minimum PT2 ', 'Extrapolated energy ', '  Excitation (a.u)  ', '  Excitation (eV)  '
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
     do k=2,N_iter_p
-      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), &
+      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter_p+1-k), extrapolated_energy(k,i), &
           extrapolated_energy(k,i) - extrapolated_energy(k,1), &
           (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
     enddo

src/mo_optimization/cipsi_orb_opt.irp.f

@@ -11,11 +11,13 @@ subroutine run_optimization
   implicit none
 
   double precision :: e_cipsi, e_opt, delta_e
+  double precision, allocatable :: Ev(:),PT2(:)
   integer :: nb_iter,i
   logical :: not_converged
   character (len=100) :: filename
 
   PROVIDE psi_det psi_coef mo_two_e_integrals_in_map ao_pseudo_integrals
+  allocate(Ev(N_states),PT2(N_states))
 
   not_converged = .True.
   nb_iter = 0
@@ -38,7 +40,7 @@ subroutine run_optimization
       print*,'' 
       print*,'********** cipsi step **********'
       ! cispi calculation
-      call run_stochastic_cipsi
+      call run_stochastic_cipsi(Ev,PT2)
 
       ! State average energy after the cipsi step
       call state_average_energy(e_cipsi)

src/mo_two_e_erf_ints/EZFIO.cfg

@@ -1,6 +0,0 @@
-[io_mo_two_e_integrals_erf]
-type: Disk_access
-doc: Read/Write MO integrals with the long range interaction from/to disk [    Write | Read | None ]
-interface: ezfio,provider,ocaml
-default: None
-

src/mo_two_e_erf_ints/NEED

@@ -1,3 +0,0 @@
-ao_two_e_erf_ints
-mo_two_e_ints
-mo_basis

src/mo_two_e_erf_ints/README.rst

@@ -1,20 +0,0 @@
-======================
-mo_two_e_erf_ints
-======================
-
-Here, all two-electron integrals (:math:`erf({\mu}_{erf} * r_{12})/r_{12}`) are computed.
-As they have 4 indices and many are zero, they are stored in a map, as defined
-in :file:`Utils/map_module.f90`.
-
-The range separation parameter :math:`{\mu}_{erf}` is the variable :option:`ao_two_e_erf_ints mu_erf`.
-
-To fetch an |MO| integral, use
-`get_mo_two_e_integral_erf(i,j,k,l,mo_integrals_map_erf)`
-
-The conventions are:
-
-* For |MO| integrals : <ij|kl> = <12|12>
-
-Be aware that it might not be the same conventions for |MO| and |AO| integrals.
-
-

src/mo_two_e_ints/EZFIO.cfg

@@ -17,3 +17,10 @@ doc: If `True`, computes all integrals except for the integrals having 3 or 4 vi
 interface: ezfio,provider,ocaml
 default: false
 
+[io_mo_two_e_integrals_erf]
+type: Disk_access
+doc: Read/Write MO integrals with the long range interaction from/to disk [    Write | Read | None ]
+interface: ezfio,provider,ocaml
+default: None
+
+

src/non_hermit_dav/biorthog.irp.f

@@ -270,7 +270,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   integer,          intent(out) :: n_real_eigv
   double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
 
-  integer                       :: i, j
+  integer                       :: i, j,k 
   integer                       :: n_good
   double precision              :: thr, thr_cut, thr_diag, thr_norm
   double precision              :: accu_d, accu_nd
@@ -278,6 +278,8 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   integer,          allocatable :: list_good(:), iorder(:)
   double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
   double precision, allocatable :: S(:,:)
+  double precision, allocatable  :: phi_1_tilde(:),phi_2_tilde(:),chi_1_tilde(:),chi_2_tilde(:)
+  allocate(phi_1_tilde(n),phi_2_tilde(n),chi_1_tilde(n),chi_2_tilde(n))
 
 
   ! -------------------------------------------------------------------------------------
@@ -301,11 +303,78 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   call lapack_diag_non_sym(n, A, WR, WI, VL, VR)
   !call lapack_diag_non_sym_new(n, A, WR, WI, VL, VR)
 
-  !print *, ' '
+
-  !print *, ' eigenvalues'
+
-  !do i = 1, n
+  print *, ' '
-  !  write(*, '(1000(F16.10,X))') WR(i), WI(i)
+  print *, ' eigenvalues'
-  !enddo
+  i = 1
+  do while(i .le. n)
+    write(*, '(I3,X,1000(F16.10,X))')i, WR(i), WI(i)
+   if(.false.)then
+    if(WI(i).ne.0.d0)then
+     print*,'*****************'
+     print*,'WARNING ! IMAGINARY EIGENVALUES !!!'
+     write(*, '(1000(F16.10,X))') WR(i), WI(i+1)
+     ! phi = VR(:,i), psi = VR(:,i+1), |Phi_i> = phi + j psi , |Phi_i+1> = phi - j psi 
+     ! chi = VL(:,i), xhi = VL(:,i+1), |Chi_i> = chi + j xhi , |Chi_i+1> = chi - j xhi 
+     !                                   
+     accu_chi_phi = 0.d0
+     accu_xhi_psi = 0.d0
+     accu_chi_psi = 0.d0
+     accu_xhi_phi = 0.d0
+     double precision :: accu_chi_phi, accu_xhi_psi, accu_chi_psi, accu_xhi_phi
+     double precision :: mat_ovlp(2,2),eigval_tmp(2),eigvec(2,2),mat_ovlp_orig(2,2)
+     do j = 1, n 
+      accu_chi_phi += VL(j,i)   * VR(j,i)
+      accu_xhi_psi += VL(j,i+1) * VR(j,i+1)
+      accu_chi_psi += VL(j,i)   * VR(j,i+1)
+      accu_xhi_phi += VL(j,i+1) * VR(j,i)
+     enddo
+     mat_ovlp_orig(1,1) = accu_chi_phi
+     mat_ovlp_orig(2,1) = accu_xhi_phi
+     mat_ovlp_orig(1,2) = accu_chi_psi
+     mat_ovlp_orig(2,2) = accu_xhi_psi
+     print*,'old overlap matrix '
+     write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,1)
+     write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,2)
+
+
+     mat_ovlp(1,1) = accu_xhi_phi
+     mat_ovlp(2,1) = accu_chi_phi
+     mat_ovlp(1,2) = accu_xhi_psi
+     mat_ovlp(2,2) = accu_chi_psi
+    !print*,'accu_chi_phi = ',accu_chi_phi
+    !print*,'accu_xhi_psi = ',accu_xhi_psi
+    !print*,'accu_chi_psi = ',accu_chi_psi
+    !print*,'accu_xhi_phi = ',accu_xhi_phi
+     print*,'new overlap matrix '
+     write(*,'(100(F16.10,X))')mat_ovlp(1:2,1)
+     write(*,'(100(F16.10,X))')mat_ovlp(1:2,2)
+     call lapack_diag(eigval_tmp,eigvec,mat_ovlp,2,2)
+     print*,'eigval_tmp(1) = ',eigval_tmp(1)
+     print*,'eigvec(1) = ',eigvec(1:2,1)
+     print*,'eigval_tmp(2) = ',eigval_tmp(2)
+     print*,'eigvec(2) = ',eigvec(1:2,2)
+     print*,'*****************'
+     phi_1_tilde = 0.d0
+     phi_2_tilde = 0.d0
+     chi_1_tilde = 0.d0
+     chi_2_tilde = 0.d0
+     do j = 1, n
+      phi_1_tilde(j) += VR(j,i) * eigvec(1,1) + VR(j,i+1) * eigvec(2,1)
+      phi_2_tilde(j) += VR(j,i) * eigvec(1,2) + VR(j,i+1) * eigvec(2,2)
+      chi_1_tilde(j) += VL(j,i+1) * eigvec(1,1) + VL(j,i) * eigvec(2,1)
+      chi_2_tilde(j) += VL(j,i+1) * eigvec(1,2) + VL(j,i) * eigvec(2,2)
+     enddo
+     VR(1:n,i)   = phi_1_tilde(1:n)
+     VR(1:n,i+1) = phi_2_tilde(1:n)
+!     Vl(1:n,i)   = -chi_1_tilde(1:n)
+!     Vl(1:n,i+1) = chi_2_tilde(1:n)
+     i+=1 
+    endif
+   endif
+    i+=1
+  enddo
   !print *, ' right eigenvect bef' 
   !do i = 1, n
   !  write(*, '(1000(F16.10,X))') VR(:,i)
@@ -331,7 +400,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
   !thr    = 100d0
   thr    = Im_thresh_tcscf
   do i = 1, n
-    !print*, 'Re(i) + Im(i)', WR(i), WI(i)
+    print*, 'Re(i) + Im(i)', WR(i), WI(i)
     if(dabs(WI(i)) .lt. thr) then
       n_good += 1
     else
@@ -405,7 +474,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
 
   if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then
 
-    !print *, ' lapack vectors are normalized and bi-orthogonalized'
+    print *, ' lapack vectors are normalized and bi-orthogonalized'
     deallocate(S)
     return
 
@@ -422,13 +491,14 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
 
   else
 
-    !print *, ' lapack vectors are not normalized neither bi-orthogonalized'
+    print *, ' lapack vectors are not normalized neither bi-orthogonalized'
 
     ! ---
 
 !   call impose_orthog_degen_eigvec(n, eigval, reigvec)
 !   call impose_orthog_degen_eigvec(n, eigval, leigvec)
 
+    call reorder_degen_eigvec(n, eigval, leigvec, reigvec)
     call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)
 
 

src/non_hermit_dav/lapack_diag_non_hermit.irp.f

@@ -1857,7 +1857,7 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
   integer                       :: i, j
   double precision, allocatable :: SS(:,:)
 
-  !print *, ' check bi-orthogonality'
+  print *, ' check bi-orthogonality'
 
   ! ---
 
@@ -1865,10 +1865,10 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
             , Vl, size(Vl, 1), Vr, size(Vr, 1) &
             , 0.d0, S, size(S, 1) )
 
-  !print *, ' overlap matrix:'
+  print *, ' overlap matrix:'
-  !do i = 1, m
+  do i = 1, m
-  !  write(*,'(1000(F16.10,X))') S(i,:)
+    write(*,'(1000(F16.10,X))') S(i,:)
-  !enddo
+  enddo
 
   accu_d  = 0.d0
   accu_nd = 0.d0
@@ -1883,8 +1883,8 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
   enddo
   accu_nd = dsqrt(accu_nd) / dble(m)
 
-  !print *, ' accu_nd = ', accu_nd
+  print *, ' accu_nd = ', accu_nd
-  !print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
+  print *, ' accu_d  = ', dabs(accu_d-dble(m))/dble(m)
 
   ! ---
 
@@ -1944,6 +1944,96 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S)
 end subroutine check_orthog
 
 ! ---
+subroutine reorder_degen_eigvec(n, e0, L0, R0)
+
+  implicit none
+
+  integer,          intent(in)    :: n
+  double precision, intent(in)    :: e0(n)
+  double precision, intent(inout) :: L0(n,n), R0(n,n)
+
+  logical                         :: complex_root
+  integer                         :: i, j, k, m
+  double precision                :: ei, ej, de, de_thr
+  double precision                :: accu_d, accu_nd
+  integer,          allocatable   :: deg_num(:)
+  double precision, allocatable   :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
+
+  ! ---
+
+  allocate( deg_num(n) )
+  do i = 1, n
+    deg_num(i) = 1
+  enddo
+
+  de_thr = thr_degen_tc
+
+  do i = 1, n-1
+    ei = e0(i)
+
+    ! already considered in degen vectors
+    if(deg_num(i).eq.0) cycle
+
+    do j = i+1, n
+      ej = e0(j)
+      de = dabs(ei - ej)
+
+      if(de .lt. de_thr) then
+        deg_num(i) = deg_num(i) + 1 
+        deg_num(j) = 0
+      endif
+
+    enddo
+  enddo
+  
+  do i = 1, n
+    if(deg_num(i) .gt. 1) then
+      print *, ' degen on', i, deg_num(i), e0(i)
+    endif
+  enddo
+
+  ! ---
+
+  do i = 1, n
+    m = deg_num(i)
+
+    if(m .gt. 1) then
+  
+      allocate(L(n,m))
+      allocate(R(n,m),S(m,m))
+
+      do j = 1, m
+        L(1:n,j) = L0(1:n,i+j-1)
+        R(1:n,j) = R0(1:n,i+j-1)
+      enddo
+
+      call dgemm( 'T', 'N', m, m, n, 1.d0      &
+                , L, size(L, 1), R, size(R, 1) &
+                , 0.d0, S, size(S, 1) )
+      print*,'Overlap matrix '
+      accu_nd = 0.D0
+      do j = 1, m
+       write(*,'(100(F16.10,X))')S(1:m,j)
+       do k = 1, m
+        if(j==k)cycle
+        accu_nd += dabs(S(j,k))
+       enddo
+      enddo
+      print*,'accu_nd = ',accu_nd
+!      if(accu_nd .gt.1.d-10)then
+!        stop
+!      endif
+      do j = 1, m
+        L0(1:n,i+j-1) = L(1:n,j)
+        R0(1:n,i+j-1) = R(1:n,j)
+      enddo
+
+      deallocate(L, R,S)
+
+    endif
+  enddo
+
+end subroutine reorder_degen_eigvec 
 
 subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
@@ -1987,11 +2077,11 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
     enddo
   enddo
   
-  !do i = 1, n
+  do i = 1, n
-  !  if(deg_num(i) .gt. 1) then
+    if(deg_num(i) .gt. 1) then
-  !    print *, ' degen on', i, deg_num(i), e0(i)
+      print *, ' degen on', i, deg_num(i), e0(i)
-  !  endif
+    endif
-  !enddo
+  enddo
 
   ! ---
 
@@ -2010,7 +2100,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
       ! ---
 
-      call impose_orthog_svd(n, m, L)
+!      call impose_orthog_svd(n, m, L)
       call impose_orthog_svd(n, m, R)
       !call impose_orthog_GramSchmidt(n, m, L)
       !call impose_orthog_GramSchmidt(n, m, R)
@@ -2031,6 +2121,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
       !deallocate(S, S_inv_half)
 
       call impose_biorthog_svd(n, m, L, R)
+!     call impose_biorthog_inverse(n, m, L, R)
 
       !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
 
@@ -2045,6 +2136,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
 
     endif
   enddo
+!  call impose_biorthog_inverse(n, n, L0, R0)
 
 end subroutine impose_biorthog_degen_eigvec 
 
@@ -2420,10 +2512,10 @@ subroutine impose_biorthog_svd(n, m, L, R)
             , L, size(L, 1), R, size(R, 1) &
             , 0.d0, S, size(S, 1) )
 
-  !print *, ' overlap bef SVD: '
+  print *, ' overlap bef SVD: '
-  !do i = 1, m
+  do i = 1, m
-  !  write(*, '(1000(F16.10,X))') S(i,:)
+    write(*, '(1000(F16.10,X))') S(i,:)
-  !enddo
+  enddo
 
   ! ---
  
@@ -2495,10 +2587,10 @@ subroutine impose_biorthog_svd(n, m, L, R)
             , L, size(L, 1), R, size(R, 1) &
             , 0.d0, S, size(S, 1) )
 
-  !print *, ' overlap aft SVD: '
+  print *, ' overlap aft SVD: '
-  !do i = 1, m
+  do i = 1, m
-  !  write(*, '(1000(F16.10,X))') S(i,:)
+    write(*, '(1000(F16.10,X))') S(i,:)
-  !enddo
+  enddo
 
   deallocate(S)
 
@@ -2506,6 +2598,50 @@ subroutine impose_biorthog_svd(n, m, L, R)
 
 end subroutine impose_biorthog_svd
 
+subroutine impose_biorthog_inverse(n, m, L, R)
+
+  implicit none
+
+  integer,          intent(in)    :: n, m
+  double precision, intent(inout) :: L(n,m)
+  double precision, intent(in)    :: R(n,m)
+  double precision, allocatable :: Lt(:,:),S(:,:)
+  integer :: i,j
+  allocate(Lt(m,n))
+  allocate(S(m,m))
+
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , L, size(L, 1), R, size(R, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' overlap bef SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  call get_pseudo_inverse(R,n,n,m,Lt,m,1.d-6)
+  do i = 1, m
+   do j = 1, n
+    L(j,i) = Lt(i,j)
+   enddo
+  enddo
+  ! ---
+
+  call dgemm( 'T', 'N', m, m, n, 1.d0      &
+            , L, size(L, 1), R, size(R, 1) &
+            , 0.d0, S, size(S, 1) )
+
+  print *, ' overlap aft SVD: '
+  do i = 1, m
+    write(*, '(1000(F16.10,X))') S(i,:)
+  enddo
+
+  deallocate(S,Lt)
+
+
+end subroutine impose_biorthog_svd
+
+
 ! ---
 
 subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr)

src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f

@@ -17,11 +17,25 @@ program tc_natorb_bi_ortho
   my_n_pt_a_grid = tc_grid1_a
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
+  if(j1b_type .ge. 100) then
+    my_extra_grid_becke  = .True.
+    PROVIDE tc_grid2_a tc_grid2_r
+    my_n_pt_r_extra_grid = tc_grid2_r
+    my_n_pt_a_extra_grid = tc_grid2_a
+    touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
+
+    call write_int(6, my_n_pt_r_extra_grid, 'radial  internal grid over')
+    call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
+  endif
+
+
+
   read_wf = .True.
   touch read_wf
 
   call print_energy_and_mos()
   call save_tc_natorb()
+  call print_angles_tc()
   !call minimize_tc_orb_angles()
 
 end
@@ -35,9 +49,12 @@ subroutine save_tc_natorb()
   print*,'Saving the natorbs '
 
   provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao
+  mo_l_coef = natorb_tc_leigvec_ao 
+  mo_r_coef = natorb_tc_reigvec_ao 
+  touch mo_l_coef mo_r_coef 
 
-  call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao)
+  call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
-  call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao)
+  call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
   call save_ref_determinant_nstates_1()
   call ezfio_set_determinants_read_wf(.False.)
 

src/tc_bi_ortho/tc_bi_ortho_prop.irp.f

@@ -34,4 +34,19 @@ subroutine test
  do i= 1, 3
   print*,tc_bi_ortho_dipole(i,1)
  enddo
+ integer, allocatable           :: occ(:,:)
+ integer                        :: n_occ_ab(2)
+ allocate(occ(N_int*bit_kind_size,2))
+ call bitstring_to_list_ab(ref_bitmask, occ, n_occ_ab, N_int)
+ integer :: ispin,j,jorb
+ double precision :: accu
+ accu = 0.d0
+ do ispin=1, 2
+  do i = 1, n_occ_ab(ispin)
+   jorb = occ(i,ispin)
+   accu += mo_bi_orth_bipole_z(jorb,jorb)
+  enddo
+ enddo
+ print*,'accu = ',accu
+ 
 end

src/tc_bi_ortho/tc_natorb.irp.f

@@ -29,9 +29,22 @@
     write(*, '(100(F16.10,X))') -dm_tmp(:,i)
   enddo
 
+  print *, ' Transition density matrix AO'
+  do i = 1, ao_num
+    write(*, '(100(F16.10,X))') tc_transition_matrix_ao(:,i,1,1)
+  enddo
+  stop
+
   thr_d   = 1.d-6
   thr_nd  = 1.d-6
   thr_deg = 1.d-3
+  do i = 1, mo_num
+   do j = 1, mo_num
+    if(dabs(dm_tmp(j,i)).lt.thr_d)then
+     dm_tmp(j,i) = 0.d0
+    endif
+   enddo
+  enddo
 !  if(n_core_orb.ne.0)then
 !   call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & 
 !                                    , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)

src/tc_bi_ortho/tc_prop.irp.f

@@ -90,6 +90,7 @@
    enddo
   enddo
  enddo
+ print*,'tc_bi_ortho_dipole(3) elec = ',tc_bi_ortho_dipole(3,1)
 
  nuclei_part = 0.d0
  do m = 1, 3

src/tc_scf/routines_rotates.irp.f

@@ -402,6 +402,7 @@ subroutine print_energy_and_mos(good_angles)
   print *, ' TC energy              = ', TC_HF_energy
   print *, ' TC SCF energy gradient = ', grad_non_hermit
   print *, ' Max angle Left/right   = ', max_angle_left_right
+  call print_angles_tc()
 
   if(max_angle_left_right .lt. thresh_lr_angle) then
     print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'

src/tools/NEED

@@ -1,5 +1,4 @@
 fci
-mo_two_e_erf_ints
 aux_quantities
 hartree_fock
 two_body_rdm

src/tools/attachement_orb.irp.f

@@ -0,0 +1,168 @@
+program molden_detachment_attachment
+ implicit none
+ read_wf=.True.
+ touch read_wf
+ call molden_attachment
+end
+
+subroutine molden_attachment
+  implicit none
+  BEGIN_DOC
+  ! Produces a Molden file
+  END_DOC
+  character*(128)                :: output
+  integer                        :: i_unit_output,getUnitAndOpen
+  integer                        :: i,j,k,l
+  double precision, parameter :: a0 = 0.529177249d0
+
+  PROVIDE ezfio_filename
+
+  output=trim(ezfio_filename)//'.attachement.mol'
+  print*,'output = ',trim(output)
+
+  i_unit_output = getUnitAndOpen(output,'w')
+
+  write(i_unit_output,'(A)') '[Molden Format]'
+
+  write(i_unit_output,'(A)') '[Atoms] Angs'
+  do i = 1, nucl_num
+    write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))')                   &
+        trim(element_name(int(nucl_charge(i)))),                     &
+        i,                                                           &
+        int(nucl_charge(i)),                                         &
+        nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
+  enddo
+
+  write(i_unit_output,'(A)') '[GTO]'
+
+  character*(1)                  :: character_shell
+  integer                        :: i_shell,i_prim,i_ao
+  integer                        :: iorder(ao_num)
+  integer                        :: nsort(ao_num)
+
+  i_shell = 0
+  i_prim = 0
+  do i=1,nucl_num
+    write(i_unit_output,*) i, 0
+    do j=1,nucl_num_shell_aos(i)
+      i_shell +=1
+      i_ao = nucl_list_shell_aos(i,j)
+      character_shell = trim(ao_l_char(i_ao))
+      write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
+      do k = 1, ao_prim_num(i_ao)
+        i_prim +=1
+        write(i_unit_output,'(ES20.10,2X,ES20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
+      enddo
+      l = i_ao
+      do while ( ao_l(l) == ao_l(i_ao) )
+        nsort(l) = i*10000 + j*100
+        l += 1
+        if (l > ao_num) exit
+      enddo
+    enddo
+    write(i_unit_output,*)''
+  enddo
+
+
+  do i=1,ao_num
+    iorder(i) = i
+    ! p
+    if      ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 3
+    ! d
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    ! f
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 10
+    ! g
+    else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 1
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 2
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
+      nsort(i) += 3
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 4
+    else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 5
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 6
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 7
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 8
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
+      nsort(i) += 9
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
+      nsort(i) += 10
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 11
+    else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 12
+    else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 13
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
+      nsort(i) += 14
+    else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
+      nsort(i) += 15
+    endif
+  enddo
+
+  call isort(nsort,iorder,ao_num)
+  write(i_unit_output,'(A)') '[MO]'
+  integer :: istate
+  istate = 2
+   do i=1,n_dettachment(istate)
+     write (i_unit_output,*) 'Sym= 1'
+     write (i_unit_output,*) 'Ene=', dettachment_numbers_sorted(i,istate)
+     write (i_unit_output,*) 'Spin= Alpha'
+     write (i_unit_output,*) 'Occup=', dettachment_numbers_sorted(i,istate)
+     do j=1,ao_num
+       write(i_unit_output, '(I6,2X,ES20.10)') j, dettachment_orbitals(iorder(j),i,istate)
+     enddo
+   enddo
+   do i=1,n_attachment(istate)
+     write (i_unit_output,*) 'Sym= 1'
+     write (i_unit_output,*) 'Ene=', attachment_numbers_sorted(i,istate)
+     write (i_unit_output,*) 'Spin= Alpha'
+     write (i_unit_output,*) 'Occup=', attachment_numbers_sorted(i,istate)
+     do j=1,ao_num
+       write(i_unit_output, '(I6,2X,ES20.10)') j, attachment_orbitals(iorder(j),i,istate)
+     enddo
+   enddo
+  close(i_unit_output)
+end
+

src/two_body_rdm/state_av_act_2rdm.irp.f

@@ -123,7 +123,7 @@
  state_av_act_2_rdm_spin_trace_mo = state_av_act_2_rdm_ab_mo & 
                                   + state_av_act_2_rdm_aa_mo & 
                                   + state_av_act_2_rdm_bb_mo   
-
+!
 ! call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
  call wall_time(wall_2)

src/two_body_rdm/test_2_rdm.irp.f

@@ -2,7 +2,7 @@ program test_2_rdm
  implicit none
  read_wf = .True.
  touch read_wf 
-! call routine_active_only
+ call routine_active_only
  call routine_full_mos
 end