state average works

TODO

@@ -24,6 +24,8 @@
   * Example  : Simple Hartree-Fock program from scratch
   * Examples : subroutine example_module
 
+# enleverle psi_det_size for all complicated stuffs with dimension of psi_coef
+
 # Config file for Cray
 
 # Documentation de /etc

src/casscf/casscf.irp.f

@@ -41,7 +41,7 @@ subroutine run
     mo_occ  = occnum  
     call save_mos
     iteration += 1
-    N_det = N_det/2 
+    N_det = max(N_det/2 ,N_states)
     psi_det = psi_det_sorted
     psi_coef = psi_coef_sorted
     read_wf = .True.

src/casscf/get_energy.irp.f

@@ -15,7 +15,13 @@ end
 
 subroutine print_grad
  implicit none
- 
+ integer :: i 
+ do i = 1, nMonoEx
+  if(dabs(gradvec2(i)).gt.1.d-5)then
+   print*,''
+   print*,i,gradvec2(i),excit(:,i)
+  endif
+ enddo
 end
 
 subroutine routine_bis

src/casscf/grad_old.irp.f

@@ -0,0 +1,74 @@
+
+BEGIN_PROVIDER [real*8, gradvec_old, (nMonoEx)]
+  BEGIN_DOC
+  ! calculate the orbital gradient <Psi| H E_pq |Psi> by hand, i.e. for
+  ! each determinant I we determine the string E_pq |I> (alpha and beta
+  ! separately) and generate <Psi|H E_pq |I>
+  ! sum_I c_I <Psi|H E_pq |I> is then the pq component of the orbital
+  ! gradient
+  ! E_pq = a^+_pa_q + a^+_Pa_Q
+  END_DOC
+  implicit none
+  integer                        :: ii,tt,aa,indx,ihole,ipart,istate
+  real*8                         :: res
+  
+  do indx=1,nMonoEx
+    ihole=excit(1,indx)
+    ipart=excit(2,indx)
+    call calc_grad_elem(ihole,ipart,res)
+    gradvec_old(indx)=res
+  end do
+  
+  real*8                         :: norm_grad
+  norm_grad=0.d0
+  do indx=1,nMonoEx
+    norm_grad+=gradvec_old(indx)*gradvec_old(indx)
+  end do
+  norm_grad=sqrt(norm_grad)
+  if (bavard) then
+    write(6,*)
+    write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
+    write(6,*)
+  endif
+  
+  
+END_PROVIDER
+
+subroutine calc_grad_elem(ihole,ipart,res)
+  BEGIN_DOC
+  ! eq 18 of Siegbahn et al, Physica Scripta 1980
+  ! we calculate 2 <Psi| H E_pq | Psi>, q=hole, p=particle
+  END_DOC
+  implicit none
+  integer                        :: ihole,ipart,mu,iii,ispin,ierr,nu,istate
+  real*8                         :: res
+  integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
+  real*8                         :: i_H_psi_array(N_states),phase
+  allocate(det_mu(N_int,2))
+  allocate(det_mu_ex(N_int,2))
+  
+  res=0.D0
+  
+  do mu=1,n_det
+    ! get the string of the determinant
+    call det_extract(det_mu,mu,N_int)
+    do ispin=1,2
+      ! do the monoexcitation on it
+      call det_copy(det_mu,det_mu_ex,N_int)
+      call do_signed_mono_excitation(det_mu,det_mu_ex,nu             &
+          ,ihole,ipart,ispin,phase,ierr)
+      if (ierr.eq.1) then
+        call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int                &
+            ,N_det,N_det,N_states,i_H_psi_array)
+        do istate=1,N_states
+          res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
+        end do
+      end if
+    end do
+  end do
+  
+  ! state-averaged gradient
+  res*=2.D0/dble(N_states)
+  
+end subroutine calc_grad_elem
+

src/casscf/gradient.irp.f

@@ -60,79 +60,6 @@ END_PROVIDER
   
 END_PROVIDER
 
-BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
-  BEGIN_DOC
-  ! calculate the orbital gradient <Psi| H E_pq |Psi> by hand, i.e. for
-  ! each determinant I we determine the string E_pq |I> (alpha and beta
-  ! separately) and generate <Psi|H E_pq |I>
-  ! sum_I c_I <Psi|H E_pq |I> is then the pq component of the orbital
-  ! gradient
-  ! E_pq = a^+_pa_q + a^+_Pa_Q
-  END_DOC
-  implicit none
-  integer                        :: ii,tt,aa,indx,ihole,ipart,istate
-  real*8                         :: res
-  
-  do indx=1,nMonoEx
-    ihole=excit(1,indx)
-    ipart=excit(2,indx)
-    call calc_grad_elem(ihole,ipart,res)
-    gradvec(indx)=res
-  end do
-  
-  real*8                         :: norm_grad
-  norm_grad=0.d0
-  do indx=1,nMonoEx
-    norm_grad+=gradvec(indx)*gradvec(indx)
-  end do
-  norm_grad=sqrt(norm_grad)
-  if (bavard) then
-    write(6,*)
-    write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
-    write(6,*)
-  endif
-  
-  
-END_PROVIDER
-
-subroutine calc_grad_elem(ihole,ipart,res)
-  BEGIN_DOC
-  ! eq 18 of Siegbahn et al, Physica Scripta 1980
-  ! we calculate 2 <Psi| H E_pq | Psi>, q=hole, p=particle
-  END_DOC
-  implicit none
-  integer                        :: ihole,ipart,mu,iii,ispin,ierr,nu,istate
-  real*8                         :: res
-  integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
-  real*8                         :: i_H_psi_array(N_states),phase
-  allocate(det_mu(N_int,2))
-  allocate(det_mu_ex(N_int,2))
-  
-  res=0.D0
-  
-  do mu=1,n_det
-    ! get the string of the determinant
-    call det_extract(det_mu,mu,N_int)
-    do ispin=1,2
-      ! do the monoexcitation on it
-      call det_copy(det_mu,det_mu_ex,N_int)
-      call do_signed_mono_excitation(det_mu,det_mu_ex,nu             &
-          ,ihole,ipart,ispin,phase,ierr)
-      if (ierr.eq.1) then
-        call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int                &
-            ,N_det,N_det,N_states,i_H_psi_array)
-        do istate=1,N_states
-          res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
-        end do
-      end if
-    end do
-  end do
-  
-  ! state-averaged gradient
-  res*=2.D0/dble(N_states)
-  
-end subroutine calc_grad_elem
-
 BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
   BEGIN_DOC
   ! calculate the orbital gradient <Psi| H E_pq |Psi> from density

src/casscf/swap_orb.irp.f

@@ -22,13 +22,11 @@
  do i = 1, nMonoEx
   iorder(i)  = i
   vec_tmp(i) = -dabs(SXvector_lowest(i))
- !print*,'vec_tmp(i) = ',i,vec_tmp(i)
  enddo
  call dsort(vec_tmp,iorder,nMonoEx)
  n_max_overlap = 0
  do i = 1, nMonoEx
   if(dabs(vec_tmp(i)).gt.thresh_overlap_switch)then
- ! print*,vec_tmp(i),iorder(i)
    n_max_overlap += 1
    max_overlap(n_max_overlap) = iorder(i)
   endif
@@ -107,7 +105,9 @@
   endif
  enddo
 
- print*,'n_orb_swap = ',n_orb_swap
+ if(n_orb_swap.gt.0)then
+  print*,'n_orb_swap = ',n_orb_swap
+ endif
  do i = 1, n_orb_swap
   print*,'imono = ',index_orb_swap(i)
   print*,orb_swap(1,i),'-->',orb_swap(2,i)

src/determinants/EZFIO.cfg

@@ -38,7 +38,7 @@ default: True
 type: integer
 doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi))
 interface: ezfio,provider,ocaml
-default: 1
+default: 2
 
 [weight_selection]
 type: integer

src/mo_basis/mos.irp.f

@@ -91,7 +91,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ]
       enddo
     enddo
   endif
-
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ]

src/two_body_rdm/orb_range_2_rdm.irp.f

@@ -79,7 +79,9 @@
  double precision :: wall_0,wall_1
  call wall_time(wall_0)
  print*,'providing the  state average TWO-RDM ...'
- call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
+ print*,'psi_det_size = ',psi_det_size
+ print*,'N_det        = ',N_det
+ call orb_range_two_rdm_state_av(state_av_act_two_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,list_act_reverse,state_weights,ispin,psi_coef,N_states,size(psi_coef,1))
 
  call wall_time(wall_1)
  print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0

src/two_body_rdm/orb_range_routines.irp.f

@@ -147,6 +147,7 @@ subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,l
    do i=1,maxab
      idx0(i) = i
    enddo
+
    
    ! Prepare the array of all alpha single excitations
    ! -------------------------------------------------