added the possibility to diagonalize the SXmatrix with Davidson

devel/casscf/EZFIO.cfg

@@ -29,6 +29,12 @@ doc: If |true|, only the DIAGONAL part of the hessian is retained for the CASSCF
 interface: ezfio,provider,ocaml
 default: False
 
+[hess_cv_cv]
+type: logical                                                                                                                                
+doc: If |true|, the core-virtual - core-virtual part of the hessian is computed 
+interface: ezfio,provider,ocaml
+default: True 
+
 
 [level_shift_casscf]
 type: Positive_float

devel/casscf/NEED

@@ -2,3 +2,4 @@ cipsi
 selectors_full
 generators_cas
 two_body_rdm
+dav_general_mat

devel/casscf/casscf.irp.f

@@ -6,8 +6,10 @@ program casscf
   call reorder_orbitals_for_casscf
   no_vvvv_integrals = .True.
   touch no_vvvv_integrals
-  pt2_max = 0.02
+  pt2_max = 0.005
   SOFT_TOUCH pt2_max
+  n_det_max_full = 500
+  touch n_det_max_full
   call run_stochastic_cipsi
   call run
 end
@@ -31,10 +33,13 @@ subroutine run
     energy = eone+etwo+ecore
 
     call write_time(6)
-    call write_int(6,iteration,'CAS-SCF iteration')
+    call write_int(6,iteration,'CAS-SCF iteration = ')
-    call write_double(6,energy,'CAS-SCF energy')
+    call write_double(6,energy,'CAS-SCF energy = ')
-    call write_double(6,norm_grad_vec2,'Norm of gradients')
+    call write_double(6,norm_grad_vec2,'Norm of gradients = ')
-    call write_double(6,energy_improvement, 'Predicted energy improvement')
+    call write_double(6,norm_grad_vec2_tab(1), '     Core-active  gradients = ')
+    call write_double(6,norm_grad_vec2_tab(2), '     Core-virtual gradients = ')
+    call write_double(6,norm_grad_vec2_tab(3), '   Active-virtual gradients = ')
+    call write_double(6,energy_improvement, 'Predicted energy improvement = ')
 
     converged = dabs(energy_improvement) < thresh_scf
     pt2_max = dabs(energy_improvement / pt2_relative_error)

devel/casscf/dav_sx_mat.irp.f

@@ -0,0 +1,44 @@
+
+
+subroutine davidson_diag_sx_mat(N_st, u_in, energies)
+ implicit none
+ integer, intent(in) :: N_st
+ double precision, intent(out) :: u_in(nMonoEx+1,n_states_diag), energies(N_st)
+ integer :: i,j,N_st_tmp, dim_in, sze, N_st_diag_in
+ integer, allocatable :: list_guess(:)
+ double precision, allocatable :: H_jj(:)
+ logical    :: converged                                                                                                
+ N_st_diag_in = n_states_diag
+ provide SXmatrix
+ sze = nMonoEx+1
+ dim_in = sze
+ allocate(H_jj(sze), list_guess(sze))
+ H_jj(1) = 0.d0
+ N_st_tmp = 1
+ list_guess(1) = 1
+ do j = 2, nMonoEx+1
+  H_jj(j) = SXmatrix(j,j)
+  if(H_jj(j).lt.0.d0)then 
+   list_guess(N_st_tmp) = j
+   N_st_tmp += 1
+  endif
+ enddo
+ if(N_st_tmp .ne. N_st)then
+  print*,'Pb in davidson_diag_sx_mat'
+  print*,'N_st_tmp .ne. N_st'
+  print*,N_st_tmp, N_st
+  stop
+ endif
+ print*,'Number of possibly interesting states = ',N_st
+ print*,'Corresponding diagonal elements of the SX matrix '
+ u_in = 0.d0
+ do i = 1, N_st
+  j = list_guess(i) 
+  print*,'i,j',i,j
+  print*,'SX(i,i) = ',H_jj(j)
+  u_in(j,i) = 1.d0
+ enddo
+ call davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,SXmatrix)
+ print*,'energies = ',energies
+
+end

devel/casscf/gradient.irp.f

@@ -95,6 +95,7 @@ END_PROVIDER
 
  BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
 &BEGIN_PROVIDER [real*8, norm_grad_vec2]
+&BEGIN_PROVIDER [real*8, norm_grad_vec2_tab, (3)]
   BEGIN_DOC
   ! calculate the orbital gradient <Psi| H E_pq |Psi> from density
   ! matrices and integrals; Siegbahn et al, Phys Scr 1980
@@ -105,10 +106,12 @@ END_PROVIDER
   real*8                         :: gradvec_it,gradvec_ia,gradvec_ta
   
   indx=0
+  norm_grad_vec2_tab = 0.d0
   do i=1,n_core_inact_orb
     do t=1,n_act_orb
       indx+=1
       gradvec2(indx)=gradvec_it(i,t)
+      norm_grad_vec2_tab(1) += gradvec2(indx)*gradvec2(indx)
     end do
   end do
   
@@ -116,6 +119,7 @@ END_PROVIDER
     do a=1,n_virt_orb
       indx+=1
       gradvec2(indx)=gradvec_ia(i,a)
+      norm_grad_vec2_tab(2) += gradvec2(indx)*gradvec2(indx)
     end do
   end do
   
@@ -123,6 +127,7 @@ END_PROVIDER
     do a=1,n_virt_orb
       indx+=1
       gradvec2(indx)=gradvec_ta(t,a)
+      norm_grad_vec2_tab(3) += gradvec2(indx)*gradvec2(indx)
     end do
   end do
   
@@ -130,6 +135,9 @@ END_PROVIDER
   do indx=1,nMonoEx
     norm_grad_vec2+=gradvec2(indx)*gradvec2(indx)
   end do
+  do i = 1, 3
+   norm_grad_vec2_tab(i) = dsqrt(norm_grad_vec2_tab(i))
+  enddo
   norm_grad_vec2=sqrt(norm_grad_vec2)
   if(bavard)then
    write(6,*)

devel/casscf/hessian.irp.f

@@ -453,34 +453,36 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
   !$OMP END DO NOWAIT
   !$OMP END PARALLEL
 
+
+ if(hess_cv_cv)then
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP SHARED(hessmat,n_c_v_prov,list_idx_c_v,n_core_inact_orb,n_virt_orb,mat_idx_c_v) &
  !$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx)
-
+  !$OMP DO
- !$OMP DO
+!!!!! < Core-VIRTUAL | H |Core-VIRTUAL >
-!!!! < Core-VIRTUAL | H |Core-VIRTUAL >
- ! Core-VIRTUAL excitations 
- do indx_tmp = 1, n_c_v_prov
-  indx = list_idx_c_v(1,indx_tmp)
-  i    = list_idx_c_v(2,indx_tmp)
-  a    = list_idx_c_v(3,indx_tmp)
   ! Core-VIRTUAL excitations 
-  do j = 1, n_core_inact_orb
+  do indx_tmp = 1, n_c_v_prov
-   if (i.eq.j) then
+   indx = list_idx_c_v(1,indx_tmp)
-     bstart=a
+   i    = list_idx_c_v(2,indx_tmp)
-   else
+   a    = list_idx_c_v(3,indx_tmp)
-     bstart=1
+   ! Core-VIRTUAL excitations 
-   end if
+   do j = 1, n_core_inact_orb
-   do b=bstart,n_virt_orb
+    if (i.eq.j) then
-    jndx = mat_idx_c_v(j,b)
+      bstart=a
-    hessmat(jndx,indx) = hessmat_iajb(i,a,j,b)
+    else
-    hessmat(indx,jndx) = hessmat(jndx,indx)
+      bstart=1
+    end if
+    do b=bstart,n_virt_orb
+     jndx = mat_idx_c_v(j,b)
+     hessmat(jndx,indx) = hessmat_iajb(i,a,j,b)
+     hessmat(indx,jndx) = hessmat(jndx,indx)
+    enddo
    enddo
   enddo
- enddo
  
   !$OMP END DO NOWAIT
   !$OMP END PARALLEL
+ endif
 
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP SHARED(hessmat,n_c_v_prov,n_a_v_prov,list_idx_c_v,list_idx_a_v) &

devel/casscf/neworbs.irp.f

@@ -1,4 +1,5 @@
-BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
+ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
+&BEGIN_PROVIDER [integer, n_guess_sx_mat ]
   implicit none
   BEGIN_DOC
   ! Single-excitation matrix
@@ -32,13 +33,18 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
   do i = 1, nMonoEx
    SXmatrix(i+1,i+1) += level_shift_casscf
   enddo
+  n_guess_sx_mat = 1
+  do i = 1, nMonoEx
+   if(SXmatrix(i+1,i+1).lt.0.d0 )then
+    n_guess_sx_mat += 1
+   endif
+  enddo
   if (bavard) then
     do i=2,nMonoEx
       write(6,*) ' diagonal of the Hessian : ',i,hessmat(i,i)
     end do
   end if
   
-  
 END_PROVIDER
 
  BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)]
@@ -47,19 +53,40 @@ END_PROVIDER
   BEGIN_DOC
   ! Eigenvectors/eigenvalues of the single-excitation matrix
   END_DOC
-  call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
+  if(nMonoEx+1.gt.n_det_max_full)then
+!   if(bavard)then
+    print*,'Using the Davidson algorithm to diagonalize the SXmatrix'
+!   endif
+   double precision, allocatable :: u_in(:,:),energies(:)
+   allocate(u_in(nMonoEx+1,n_states_diag),energies(n_guess_sx_mat))
+   call davidson_diag_sx_mat(n_guess_sx_mat, u_in, energies)
+   integer :: i,j
+   SXeigenvec = 0.d0
+   SXeigenval = 0.d0
+   do i = 1, n_guess_sx_mat
+    SXeigenval(i) = energies(i)
+    do j = 1, nMonoEx+1
+     SXeigenvec(j,i) = u_in(j,i)
+    enddo
+   enddo
+  else
+!   if(bavard)then
+    print*,'Diagonalize the SXmatrix with Jacobi'
+!   endif
+   call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
+  endif
   if (bavard) then
-    write(6,*) ' SXdiag : lowest 5 eigenvalues '
+    write(6,*) ' SXdiag : lowest eigenvalues '
     write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
-    if(nmonoex.gt.0)then
+    if(n_guess_sx_mat.gt.0)then
      write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2)
      write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3)
      write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4)
      write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5)
     endif
+  endif
     write(6,*)
     write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
-  endif
 END_PROVIDER
 
  BEGIN_PROVIDER [real*8, energy_improvement]
@@ -82,8 +109,8 @@ END_PROVIDER
   best_vector_ovrlp_casscf = -1000
   do i=1,nMonoEx+1
     if (SXeigenval(i).lt.0.D0) then
-      if (abs(SXeigenvec(1,i)).gt.best_overlap_casscf) then
+      if (dabs(SXeigenvec(1,i)).gt.best_overlap_casscf) then
-        best_overlap_casscf=abs(SXeigenvec(1,i))
+        best_overlap_casscf=dabs(SXeigenvec(1,i))
         best_vector_ovrlp_casscf = i
       end if
     end if