Davidson N_states_diag incrementing

configure

@@ -487,9 +487,12 @@ def create_ninja_and_rc(l_installed):
         'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")),
         'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
         'function qp_prepend_export () {',
-        '  #Prepend path $2:${!1}. Add the semicolon only if ${!1} is defined',
-        '  eval "value_1=\"\${$1}\""',
-        '  echo ${value_1:+${2}:${value_1}}',
+        'eval "value_1="\${$1}""',
+        'if [[ -z $value_1 ]] ; then',
+        '    echo "${1}=${2}:"',
+        'else',
+        '    echo "${1}=${2}:${value_1}"',
+        'fi',
         '}',
         'export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")',
         'export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)',

scripts/ezfio_interface/qp_convert_output_to_ezfio.py

@@ -28,7 +28,8 @@ except:
     sys.exit(1)
 else:
     sys.path = [ QP_ROOT + "/install/EZFIO/Python", 
-                 QP_ROOT + "/resultsFile", 
+                 QP_ROOT + "/install/resultsFile", 
+                 QP_ROOT + "/install", 
                  QP_ROOT + "/scripts"] + sys.path
 
 # ~#~#~#~#~#~ #

scripts/qp_set_frozen_core.py

@@ -12,14 +12,8 @@ if filename == '-q': filename = sys.argv[2]
 
 ezfio.set_filename(filename)
 
-if ezfio.pseudo_do_pseudo:
-  if '-q' in sys.argv:
-    print 0
-    sys.exit(0)
-  print "I will not set frozen core with pseudopotentials."
-  sys.exit(0)
-
 nb = 0
+if not ezfio.pseudo_do_pseudo:
   for charge in ezfio.nuclei_nucl_charge:
     if charge < 5:
         pass

src/Davidson/diagonalization_hs2_dressed.irp.f

@@ -65,12 +65,27 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
     enddo
   endif
 
-  call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
+  integer :: N_st_diag_local
+  double precision, allocatable :: energies_local(:), s2_out_local(:), u_in_local(:,:)
+  logical :: converged
+  converged = .False.
+  call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state,converged)
+  N_st_diag_local = N_st_diag
+  do while (.not.converged)
+    N_st_diag_local += N_states
+    allocate (energies_local(N_st_diag_local), s2_out_local(N_st_diag_local), u_in_local(sze,N_st_diag_local))
+    u_in_local(1:sze,1:N_st_diag) = u_in(1:sze,1:N_st_diag) 
+    call davidson_diag_hjj_sjj(dets_in,u_in_local,H_jj,s2_out_local,energies_local,dim_in,sze,N_st,N_st_diag_local,Nint,dressing_state,converged)
+    energies(1:N_st_diag) = energies_local(1:N_st_diag)
+    s2_out(1:N_st_diag) = s2_out_local(1:N_st_diag)
+    u_in(1:sze,1:N_st_diag) = u_in_local(1:sze,1:N_st_diag)
+    deallocate (energies_local, s2_out_local, u_in_local)
+  enddo
   deallocate (H_jj)
 end
 
 
-subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
+subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state,converged)
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -105,13 +120,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   
   integer                        :: iter
   integer                        :: i,j,k,l,m
-  logical                        :: converged
+  logical, intent(inout)         :: converged
   
   double precision, external     :: u_dot_v, u_dot_u
   
   integer                        :: k_pairs, kl
   
-  integer                        :: iter2
+  integer                        :: iter2, itertot
   double precision, allocatable  :: W(:,:),  U(:,:), S(:,:), overlap(:,:)
   double precision, allocatable  :: y(:,:), h(:,:), lambda(:), s2(:)
   double precision, allocatable  :: c(:), s_(:,:), s_tmp(:,:)
@@ -133,6 +148,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   endif
   
   itermax = max(3,min(davidson_sze_max, sze/N_st_diag))
+  itertot = 0
   
   PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse
   
@@ -220,6 +236,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
   
   do while (.not.converged)
+    itertot = itertot+1
+    if (itertot == 5) then
+      exit
+    endif
 
     do k=1,N_st_diag
       do i=1,sze

src/MO_Basis/ao_ortho_canonical.irp.f

@@ -2,9 +2,7 @@
 &BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ]
   implicit none
   BEGIN_DOC
-! matrix of the coefficients of the mos generated by the 
-! orthonormalization by the S^{-1/2} canonical transformation of the aos
-! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
+! ao_cart_to_sphe coefficients of the current basis set
   END_DOC
   integer :: i
   integer, external              :: ao_power_index
@@ -117,7 +115,6 @@ END_PROVIDER
       ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), &
       ao_ortho_canonical_num)
 
-
   else
 
     double precision, allocatable :: S(:,:)
@@ -137,6 +134,13 @@ END_PROVIDER
       S, size(S,1), &
       0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1))
 
+!integer :: j
+!do i=1,ao_num
+! do j=1,ao_num
+!  print *,  i,j, ao_ortho_canonical_coef(i,j)
+! enddo
+!enddo
+!stop
     deallocate(S)
   endif
 END_PROVIDER

src/MO_Basis/utils.irp.f

@@ -112,6 +112,7 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
   double precision, intent(in)   :: matrix(lda,n)
   
   integer :: i,j
+  double precision :: accu
   double precision, allocatable  :: mo_coef_new(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A
   
@@ -137,12 +138,16 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
   write (6,'(A)') 'Eigenvalues'
   write (6,'(A)')  '-----------'
   write (6,'(A)') ''
-  write (6,'(A)')  '======== ================'
+  write (6,'(A)')  '======== ================ ================'
+  write (6,'(A)')  '   MO       Eigenvalue       Cumulative   '
+  write (6,'(A)')  '======== ================ ================'
 
+  accu = 0.d0
   do i=1,m
-    write (6,'(I8,1X,F16.10)')  i,D(i)
+    accu = accu + D(i)
+    write (6,'(I8,1X,F16.10,1X,F16.10)')  i,D(i), accu
   enddo
-  write (6,'(A)')  '======== ================'
+  write (6,'(A)')  '======== ================ ================'
   write (6,'(A)')  ''
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))