src/casscf/bavard.irp.f

@@ -1,6 +1,6 @@
 ! -*- F90 -*- 
 BEGIN_PROVIDER [logical, bavard]
      bavard=.true.
-!     bavard=.false.
+     bavard=.false.
 END_PROVIDER
 

src/casscf/bielec.irp.f

@@ -55,6 +55,7 @@
     end do
   end do
 
+  write(6,*) ' provided integrals (PQ|xx) '
 END_PROVIDER
 
 
@@ -115,6 +116,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a
       end do
     end do
   end do
+  write(6,*) ' provided integrals (Px|xQ) '
 END_PROVIDER
 
 
@@ -144,5 +146,6 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
       end do
     end do
   end do
+  write(6,*) ' provided integrals (tu|xP) '
 END_PROVIDER
 

src/casscf/bielec_natorb.irp.f

@@ -84,6 +84,7 @@
       end do
     end do
   end do
+  write(6,*) ' transformed PQxx'
 
 END_PROVIDER
 
@@ -175,6 +176,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+
       end do
     end do
   end do
+  write(6,*) ' transformed PxxQ '
 
 END_PROVIDER
 
@@ -265,6 +267,7 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
       end do 
     end do  
   end do  
+  write(6,*) ' transformed tuvP '
 
 END_PROVIDER
 

src/casscf/casscf.irp.f

@@ -12,32 +12,24 @@ subroutine run
   implicit none
   double precision               :: energy_old, energy
   logical                        :: converged
-  integer                        :: iteration
   converged = .False.
 
   energy = 0.d0
-  mo_label = "MCSCF"
+!  do while (.not.converged)
-  iteration = 1
+    N_det = 1
-  do while (.not.converged)
+    TOUCH N_det psi_det psi_coef 
     call run_cipsi
 
+    write(6,*) ' total energy = ',eone+etwo+ecore
+    mo_label = "MCSCF"
+    mo_label = "Natural"
+    mo_coef(:,:) = NatOrbsFCI(:,:)
+    call save_mos
+
+    call driver_optorb
     energy_old = energy
     energy = eone+etwo+ecore
-
+    converged = dabs(energy - energy_old) < 1.d-10
-    call write_time(6)
+!  enddo
-    call write_int(6,iteration,'CAS-SCF iteration')
-    call write_double(6,energy,'CAS-SCF energy')
-    call write_double(6,energy_improvement, 'Predicted energy improvement')
-
-    converged = dabs(energy_improvement) < thresh_scf
-
-    mo_coef = NewOrbs
-    call save_mos
-    call map_deinit(mo_integrals_map)
-    N_det = 1
-    iteration += 1
-    FREE mo_integrals_map mo_two_e_integrals_in_map psi_det psi_coef
-    SOFT_TOUCH mo_coef N_det
-  enddo
 
 end

src/casscf/densities.irp.f

@@ -1,19 +1,70 @@
 use bitmasks
 
 BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
-  implicit none
   BEGIN_DOC
-  ! the first-order density matrix in the basis of the starting MOs.
+  ! the first-order density matrix in the basis of the starting MOs
-  ! matrix is state averaged.
+  ! matrices are state averaged
+  !
+  ! we use the spin-free generators of mono-excitations
+  ! E_pq destroys q and creates p
+  ! D_pq   =     <0|E_pq|0> = D_qp
+  !
   END_DOC
-  integer                        :: t,u
+  implicit none
+  integer                        :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
+  integer                        :: ierr
+  integer(bit_kind)              :: det_mu(N_int,2)
+  integer(bit_kind)              :: det_mu_ex(N_int,2)
+  integer(bit_kind)              :: det_mu_ex1(N_int,2)
+  integer(bit_kind)              :: det_mu_ex2(N_int,2)
+  real*8                         :: phase1,phase2,term
+  integer                        :: nu1,nu2
+  integer                        :: ierr1,ierr2
+  real*8                         :: cI_mu(N_states)
   
-  do u=1,n_act_orb
+  write(6,*) ' providing density matrices D0 and P0 '
+  
+  D0tu = 0.d0
+  
+  ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
+  do mu=1,n_det
+    call det_extract(det_mu,mu,N_int)
+    do istate=1,n_states
+      cI_mu(istate)=psi_coef(mu,istate)
+    end do
     do t=1,n_act_orb
-      D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + &
+      ipart=list_act(t)
-                  one_e_dm_mo_beta_average ( list_act(t), list_act(u) ) 
+      do u=1,n_act_orb
-    enddo
+        ihole=list_act(u)
-  enddo
+        ! apply E_tu
+        call det_copy(det_mu,det_mu_ex1,N_int)
+        call det_copy(det_mu,det_mu_ex2,N_int)
+        call do_spinfree_mono_excitation(det_mu,det_mu_ex1           &
+            ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
+        ! det_mu_ex1 is in the list
+        if (nu1.ne.-1) then
+          do istate=1,n_states
+            term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
+            D0tu(t,u)+=term
+          end do
+        end if
+        ! det_mu_ex2 is in the list
+        if (nu2.ne.-1) then
+          do istate=1,n_states
+            term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
+            D0tu(t,u)+=term
+          end do
+        end if
+      end do
+    end do
+  end do
+  
+  ! we average by just dividing by the number of states
+  do x=1,n_act_orb
+    do v=1,n_act_orb
+      D0tu(v,x)*=1.0D0/dble(N_states)
+    end do
+  end do
   
 END_PROVIDER
 
@@ -39,9 +90,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
    integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12
    integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
    
-  if (bavard) then
+   write(6,*) ' providing density matrices D0 and P0 '
-    write(6,*) ' providing density matrix P0'
-  endif
 
    P0tuvx = 0.d0
    

src/casscf/det_manip.irp.f

@@ -31,8 +31,6 @@ subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart,       &
     ! get the number in the list
     found=.false.
     nu=0
-
-    !TODO BOTTLENECK
     do while (.not.found)
       nu+=1
       if (nu.gt.N_det) then
@@ -52,6 +50,13 @@ subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart,       &
         end do
       end if
     end do
+    !         if (found) then
+    !          if (nu.eq.-1) then
+    !           write(6,*) ' image not found in the list, thus nu = ',nu
+    !          else
+    !           write(6,*) ' found in the list as No ',nu,' phase = ',phase
+    !          end if
+    !         end if
   end if
   !
   ! we found the new string, the phase, and possibly the number in the list

src/casscf/driver_optorb.irp.f

@@ -1,3 +1,32 @@
-subroutine driver_optorb
+       subroutine driver_optorb
-  implicit none
+         implicit none
-end 
+         integer :: i,j
+
+         write(6,*)
+!        write(6,*) ' <0|H|0>     (qp)      = ',psi_energy_with_nucl_rep(1)
+         write(6,*) ' energy improvement    = ',energy_improvement
+!        write(6,*) ' new energy            = ',psi_energy_with_nucl_rep(1)+energy_improvement
+         write(6,*)
+
+         write(6,*)
+         write(6,*) '   creating new orbitals '
+         do i=1,mo_num
+          write(6,*) ' Orbital No ',i
+          write(6,'(5F14.6)') (NewOrbs(j,i),j=1,mo_num)
+          write(6,*)
+         end do
+
+         mo_label = "Natural"
+         do i=1,mo_num
+          do j=1,ao_num
+           mo_coef(j,i)=NewOrbs(j,i)
+          end do
+         end do
+         call save_mos
+         call map_deinit(mo_integrals_map)
+         FREE mo_integrals_map mo_coef mo_two_e_integrals_in_map
+ 
+         write(6,*)
+         write(6,*) '   ... all done '
+
+       end 

src/casscf/gradient.irp.f

@@ -6,6 +6,7 @@ BEGIN_PROVIDER [ integer, nMonoEx ]
   END_DOC
   implicit none
   nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
+  write(6,*) ' nMonoEx = ',nMonoEx
 END_PROVIDER
 
  BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
@@ -86,11 +87,9 @@ BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
     norm_grad+=gradvec(indx)*gradvec(indx)
   end do
   norm_grad=sqrt(norm_grad)
-  if (bavard) then
+  write(6,*)
-    write(6,*)
+  write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
-    write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
+  write(6,*)
-    write(6,*)
-  endif
   
   
 END_PROVIDER
@@ -119,11 +118,17 @@ subroutine calc_grad_elem(ihole,ipart,res)
       call do_signed_mono_excitation(det_mu,det_mu_ex,nu             &
           ,ihole,ipart,ispin,phase,ierr)
       if (ierr.eq.1) then
+        !       write(6,*)
+        !       write(6,*) ' mu = ',mu
+        !       call print_det(det_mu,N_int)
+        !         write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin
+        !         call print_det(det_mu_ex,N_int)
         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
+        !         write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res
       end if
     end do
   end do
@@ -171,11 +176,9 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
     norm_grad+=gradvec2(indx)*gradvec2(indx)
   end do
   norm_grad=sqrt(norm_grad)
-  if (bavard) then
+  write(6,*)
-    write(6,*)
+  write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
-    write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
+  write(6,*)
-    write(6,*)
-  endif
   
 END_PROVIDER
 

src/casscf/hessian.irp.f

@@ -14,10 +14,8 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
   character*3                    :: iexc,jexc
   real*8                         :: res
   
-  if (bavard) then
+  write(6,*) ' providing Hessian matrix hessmat '
-    write(6,*) ' providing Hessian matrix hessmat '
+  write(6,*) '  nMonoEx = ',nMonoEx
-    write(6,*) '  nMonoEx = ',nMonoEx
-  endif
   
   do indx=1,nMonoEx
     do jndx=1,nMonoEx
@@ -34,6 +32,8 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
       jpart=excit(2,jndx)
       jexc=excit_class(jndx)
       call calc_hess_elem(ihole,ipart,jhole,jpart,res)
+      !      write(6,*) ' Hessian ',ihole,'->',ipart                 &
+      !           ,' (',iexc,')',jhole,'->',jpart,' (',jexc,')',res
       hessmat(indx,jndx)=res
       hessmat(jndx,indx)=res
     end do
@@ -198,10 +198,8 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
   real*8                         :: hessmat_iatb
   real*8                         :: hessmat_taub
   
-  if (bavard) then
+  write(6,*) ' providing Hessian matrix hessmat2 '
-    write(6,*) ' providing Hessian matrix hessmat2 '
+  write(6,*) '  nMonoEx = ',nMonoEx
-    write(6,*) '  nMonoEx = ',nMonoEx
-  endif
   
   indx=1
   do i=1,n_core_orb
@@ -216,6 +214,7 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
         do u=ustart,n_act_orb
           hessmat2(indx,jndx)=hessmat_itju(i,t,j,u)
           hessmat2(jndx,indx)=hessmat2(indx,jndx)
+          !          write(6,*) ' result I :',i,t,j,u,indx,jndx,hessmat(indx,jndx),hessmat2(indx,jndx)
           jndx+=1
         end do
       end do
@@ -295,6 +294,7 @@ real*8 function hessmat_itju(i,t,j,u)
   integer                        :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
   real*8                         :: term,t2
   
+  !      write(6,*) ' hessmat_itju ',i,t,j,u
   ii=list_core(i)
   tt=list_act(t)
   if (i.eq.j) then
@@ -340,6 +340,8 @@ real*8 function hessmat_itju(i,t,j,u)
           end do
         end do
       end do
+      !!!         write(6,*) ' direct diff ',i,t,j,u,term,term2
+      !!!         term=term2
     end if
   else
     ! it/ju
@@ -380,6 +382,7 @@ real*8 function hessmat_itja(i,t,j,a)
   integer                        :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
   real*8                         :: term
   
+  !      write(6,*) ' hessmat_itja ',i,t,j,a
   ! it/ja
   ii=list_core(i)
   tt=list_act(t)
@@ -413,6 +416,7 @@ real*8 function hessmat_itua(i,t,u,a)
   integer                        :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3
   real*8                         :: term
   
+  !      write(6,*) ' hessmat_itua ',i,t,u,a
   ii=list_core(i)
   tt=list_act(t)
   t3=t+n_core_orb
@@ -453,6 +457,7 @@ real*8 function hessmat_iajb(i,a,j,b)
   implicit none
   integer                        :: i,a,j,b,ii,aa,jj,bb
   real*8                         :: term
+  !      write(6,*) ' hessmat_iajb ',i,a,j,b
   
   ii=list_core(i)
   aa=list_virt(a)
@@ -490,6 +495,7 @@ real*8 function hessmat_iatb(i,a,t,b)
   integer                        :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
   real*8                         :: term
   
+  !      write(6,*) ' hessmat_iatb ',i,a,t,b
   ii=list_core(i)
   aa=list_virt(a)
   tt=list_act(t)
@@ -546,6 +552,7 @@ real*8 function hessmat_taub(t,a,u,b)
         end do
       end do
       term=t1+t2+t3
+      !        write(6,*) ' Hess taub ',t,a,t1,t2,t3
     else
       bb=list_virt(b)
       ! ta/tb b/=a

src/casscf/natorb.irp.f

@@ -14,12 +14,10 @@
     occnum(list_act(i))=occ_act(n_act_orb-i+1)
   end do
 
-  if (bavard) then
+  write(6,*) ' occupation numbers '
-    write(6,*) ' occupation numbers '
+  do i=1,mo_num
-    do i=1,mo_num
+    write(6,*) i,occnum(i)
-      write(6,*) i,occnum(i)
+  end do
-    end do
-  endif
 
 END_PROVIDER
 
@@ -34,12 +32,14 @@ END_PROVIDER
  
  call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
  
- if (bavard) then
+ write(6,*) ' found occupation numbers as '
-   write(6,*) ' found occupation numbers as '
+ do i=1,n_act_orb
-   do i=1,n_act_orb
+   write(6,*) i,occ_act(i)
-     write(6,*) i,occ_act(i)
+ end do
-   end do
  
+ if (bavard) then
+   !
+   
    integer                        :: nmx
    real*8                         :: xmx
    do i=1,n_act_orb
@@ -152,6 +152,7 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
       end do
     end do
   end do
+  write(6,*) ' transformed P0tuvx '
 
 END_PROVIDER
 
@@ -197,6 +198,7 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
       one_ints_no(j,list_act(p))=d(p)
     end do
   end do
+  write(6,*) ' transformed one_ints '
 END_PROVIDER
 
 
@@ -224,5 +226,148 @@ BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
       NatOrbsFCI(j,list_act(p))=d(p)
     end do
   end do
+  write(6,*) ' transformed orbitals '
 END_PROVIDER
 
+
+
+
+
+
+subroutine trf_to_natorb()
+  implicit none
+  BEGIN_DOC
+  ! save the diagonal somewhere, in inverse order
+  ! 4-index-transform the 2-particle density matrix over active orbitals
+  ! correct the bielectronic integrals
+  ! correct the monoelectronic integrals
+  ! put integrals on file, as well orbitals, and the density matrices
+  !
+  END_DOC
+  integer                        :: i,j,k,l,t,u,p,q,pp
+  real*8                         :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb)
+  
+  ! we recalculate total energies
+  write(6,*)
+  write(6,*)  ' recalculating energies after the transformation '
+  write(6,*)
+  write(6,*)
+  real*8                         :: e_one_all
+  real*8                         :: e_two_all
+  integer                        :: ii
+  integer                        :: jj
+  integer                        :: t3
+  integer                        :: tt
+  integer                        :: u3
+  integer                        :: uu
+  integer                        :: v
+  integer                        :: v3
+  integer                        :: vv
+  integer                        :: x
+  integer                        :: x3
+  integer                        :: xx
+  
+  e_one_all=0.D0
+  e_two_all=0.D0
+  do i=1,n_core_orb
+    ii=list_core(i)
+    e_one_all+=2.D0*one_ints_no(ii,ii)
+    do j=1,n_core_orb
+      jj=list_core(j)
+      e_two_all+=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
+    end do
+    do t=1,n_act_orb
+      tt=list_act(t)
+      t3=t+n_core_orb
+      e_two_all += occnum(list_act(t)) *                            &
+          (2.d0*bielec_PQxx_no(tt,tt,i,i) - bielec_PQxx_no(tt,ii,i,t3))
+    end do
+  end do
+
+    
+
+  do t=1,n_act_orb
+    tt=list_act(t)
+    e_one_all += occnum(list_act(t))*one_ints_no(tt,tt)
+    do u=1,n_act_orb
+      uu=list_act(u)
+      do v=1,n_act_orb
+        v3=v+n_core_orb
+        do x=1,n_act_orb
+          x3=x+n_core_orb
+          e_two_all  +=P0tuvx_no(t,u,v,x)*bielec_PQxx_no(tt,uu,v3,x3)
+        end do
+      end do
+    end do
+  end do
+  write(6,*) ' e_one_all = ',e_one_all
+  write(6,*) ' e_two_all = ',e_two_all
+  ecore    =nuclear_repulsion
+  ecore_bis=nuclear_repulsion
+  do i=1,n_core_orb
+    ii=list_core(i)
+    ecore    +=2.D0*one_ints_no(ii,ii)
+    ecore_bis+=2.D0*one_ints_no(ii,ii)
+    do j=1,n_core_orb
+      jj=list_core(j)
+      ecore    +=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
+      ecore_bis+=2.D0*bielec_PxxQ_no(ii,i,j,jj)-bielec_PxxQ_no(ii,j,j,ii)
+    end do
+  end do
+  eone    =0.D0
+  eone_bis=0.D0
+  etwo    =0.D0
+  etwo_bis=0.D0
+  etwo_ter=0.D0
+  do t=1,n_act_orb
+    tt=list_act(t)
+    t3=t+n_core_orb
+    eone     += occnum(list_act(t))*one_ints_no(tt,tt)
+    eone_bis += occnum(list_act(t))*one_ints_no(tt,tt)
+    do i=1,n_core_orb
+      ii=list_core(i)
+      eone     += occnum(list_act(t)) * &
+        (2.D0*bielec_PQxx_no(tt,tt,i,i ) - bielec_PQxx_no(tt,ii,i,t3))
+      eone_bis += occnum(list_act(t)) * &
+        (2.D0*bielec_PxxQ_no(tt,t3,i,ii) - bielec_PxxQ_no(tt,i ,i,tt))
+    end do
+    do u=1,n_act_orb
+      uu=list_act(u)
+      u3=u+n_core_orb
+      do v=1,n_act_orb
+        vv=list_act(v)
+        v3=v+n_core_orb
+        do x=1,n_act_orb
+          xx=list_act(x)
+          x3=x+n_core_orb
+          real*8                         :: h1,h2,h3
+          h1=bielec_PQxx_no(tt,uu,v3,x3)
+          h2=bielec_PxxQ_no(tt,u3,v3,xx)
+          h3=bielecCI_no(t,u,v,xx)
+          etwo    +=P0tuvx_no(t,u,v,x)*h1
+          etwo_bis+=P0tuvx_no(t,u,v,x)*h2
+          etwo_ter+=P0tuvx_no(t,u,v,x)*h3
+          if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
+            write(6,9901) t,u,v,x,h1,h2,h3
+            9901 format('aie: ',4I4,3E20.12)
+          end if
+        end do
+      end do
+    end do
+  end do
+  
+  write(6,*) ' energy contributions '
+  write(6,*) '     core energy       = ',ecore,' using PQxx integrals '
+  write(6,*) '     core energy (bis) = ',ecore,' using PxxQ integrals '
+  write(6,*) '     1el  energy       = ',eone ,' using PQxx integrals '
+  write(6,*) '     1el  energy (bis) = ',eone ,' using PxxQ integrals '
+  write(6,*) '     2el  energy       = ',etwo    ,' using PQxx integrals '
+  write(6,*) '     2el  energy (bis) = ',etwo_bis,' using PxxQ integrals '
+  write(6,*) '     2el  energy (ter) = ',etwo_ter,' using tuvP integrals '
+  write(6,*) ' ----------------------------------------- '
+  write(6,*) '     sum of all        = ',eone+etwo+ecore
+  write(6,*)
+  SOFT_TOUCH ecore ecore_bis eone eone_bis etwo etwo_bis etwo_ter 
+  
+end subroutine trf_to_natorb
+

src/casscf/neworbs.irp.f

@@ -1,178 +1,222 @@
+! -*- F90 -*-
 BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
-  implicit none
+      implicit none
-  BEGIN_DOC
+      integer :: i,j
-  ! Single-excitation matrix
+      do i=1,nMonoEx+1
-  END_DOC
+       do j=1,nMonoEx+1
-  
+        SXmatrix(i,j)=0.D0
-  integer                        :: i,j
+       end do
-  
+      end do
-  do i=1,nMonoEx+1
+
-    do j=1,nMonoEx+1
+      do i=1,nMonoEx
-      SXmatrix(i,j)=0.D0
+       SXmatrix(1,i+1)=gradvec2(i)
-    end do
+       SXmatrix(1+i,1)=gradvec2(i)
-  end do
+      end do
-  
+
-  do i=1,nMonoEx
+      do i=1,nMonoEx
-    SXmatrix(1,i+1)=gradvec2(i)
+       do j=1,nMonoEx
-    SXmatrix(1+i,1)=gradvec2(i)
+        SXmatrix(i+1,j+1)=hessmat2(i,j)
-  end do
+        SXmatrix(j+1,i+1)=hessmat2(i,j)
-  
+       end do
-  do i=1,nMonoEx
+      end do
-    do j=1,nMonoEx
+
-      SXmatrix(i+1,j+1)=hessmat2(i,j)
+   if (bavard) then
-      SXmatrix(j+1,i+1)=hessmat2(i,j)
+      do i=2,nMonoEx+1
-    end do
+       write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
-  end do
+      end do
-  
+   end if
-  if (bavard) then
+      
-    do i=2,nMonoEx+1
+
-      write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
-    end do
-  end if
-  
-  
 END_PROVIDER
 
  BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)]
 &BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)]
-  implicit none
+ END_PROVIDER
-  BEGIN_DOC
-  ! Eigenvectors/eigenvalues of the single-excitation matrix
-  END_DOC
-  call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
-END_PROVIDER
 
  BEGIN_PROVIDER [real*8, SXvector, (nMonoEx+1)]
 &BEGIN_PROVIDER [real*8, energy_improvement]
-  implicit none
+      implicit none
-  BEGIN_DOC
+      integer :: ierr,matz,i
-  ! Best eigenvector of the single-excitation matrix
+      real*8 :: c0
-  END_DOC
-  integer                        :: ierr,matz,i
-  real*8                         :: c0
-  
-  if (bavard) then
-    write(6,*) ' SXdiag : lowest 5 eigenvalues '
-    write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
-    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)
-    write(6,*)
-    write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
-  endif
-  energy_improvement = SXeigenval(1)
-  
-  integer                        :: best_vector
-  real*8                         :: best_overlap
-  best_overlap=0.D0
-  do i=1,nMonoEx+1
-    if (SXeigenval(i).lt.0.D0) then
-      if (abs(SXeigenvec(1,i)).gt.best_overlap) then
-        best_overlap=abs(SXeigenvec(1,i))
-        best_vector=i
-      end if
-    end if
-  end do
-  
-  energy_improvement = SXeigenval(best_vector)
-  
-  if (bavard) then
-    write(6,*) ' SXdiag : eigenvalue for best overlap with '
-    write(6,*) '  previous orbitals = ',SXeigenval(best_vector)
-    write(6,*) ' weight of the 1st element ',c0
-  endif
 
-  c0=SXeigenvec(1,best_vector)
+      call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
+      write(6,*) ' SXdiag : lowest 5 eigenvalues '
+      write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
+      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)
+      write(6,*) 
+      write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
+      energy_improvement = SXeigenval(1)
 
-  do i=1,nMonoEx+1
+integer :: best_vector
-    SXvector(i)=SXeigenvec(i,best_vector)/c0
+real*8 :: best_overlap
-  end do
+      best_overlap=0.D0
+      do i=1,nMonoEx+1
+       if (SXeigenval(i).lt.0.D0) then
+        if (abs(SXeigenvec(1,i)).gt.best_overlap) then
+         best_overlap=abs(SXeigenvec(1,i))
+         best_vector=i
+        end if
+       end if
+      end do
+
+      write(6,*) ' SXdiag : eigenvalue for best overlap with ' 
+      write(6,*) '  previous orbitals = ',SXeigenval(best_vector)
+      energy_improvement = SXeigenval(best_vector)
+      
+      c0=SXeigenvec(1,best_vector)
+      write(6,*) ' weight of the 1st element ',c0
+      do i=1,nMonoEx+1
+       SXvector(i)=SXeigenvec(i,best_vector)/c0
+!      write(6,*) ' component No ',i,' : ',SXvector(i)
+      end do
 
-  
 END_PROVIDER
 
 
 BEGIN_PROVIDER [real*8, NewOrbs, (ao_num,mo_num) ]
-  implicit none
+      implicit none
-  BEGIN_DOC
+      integer :: i,j,ialph
-  ! Updated orbitals
+
-  END_DOC
+! form the exponential of the Orbital rotations
-  integer                        :: i,j,ialph
+      call get_orbrotmat
-  
+! form the new orbitals
-  call dgemm('N','T', ao_num,mo_num,mo_num,1.d0,                     &
+      do i=1,ao_num
-      NatOrbsFCI, size(NatOrbsFCI,1),                                &
+       do j=1,mo_num
-      Umat, size(Umat,1), 0.d0,                                      &
+        NewOrbs(i,j)=0.D0
-      NewOrbs, size(NewOrbs,1))
+       end do
-  
+      end do
+
+      do ialph=1,ao_num
+       do i=1,mo_num
+        wrkline(i)=mo_coef(ialph,i)
+       end do
+       do i=1,mo_num
+        do j=1,mo_num
+         NewOrbs(ialph,i)+=Umat(i,j)*wrkline(j)
+        end do
+       end do
+      end do
+
 END_PROVIDER
 
-BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
+ BEGIN_PROVIDER [real*8, Tpotmat, (mo_num,mo_num) ]
-  implicit none
+&BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
-  BEGIN_DOC
+&BEGIN_PROVIDER [real*8, wrkline, (mo_num) ]
-  ! Orbital rotation matrix
+&BEGIN_PROVIDER [real*8, Tmat, (mo_num,mo_num) ]
-  END_DOC
+END_PROVIDER
-  integer                        :: i,j,indx,k,iter,t,a,ii,tt,aa
+
-  logical                        :: converged
+      subroutine get_orbrotmat
-  
+      implicit none
-  real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num) 
+      integer :: i,j,indx,k,iter,t,a,ii,tt,aa
-  real*8 :: Tmat(mo_num,mo_num) 
+      real*8 :: sum
-  real*8 :: f
+      logical :: converged
-  
+
-  ! the orbital rotation matrix T
+
-  Tmat(:,:)=0.D0
+! the orbital rotation matrix T
-  indx=1
+      do i=1,mo_num
-  do i=1,n_core_orb
+       do j=1,mo_num
-    ii=list_core(i)
+        Tmat(i,j)=0.D0
-    do t=1,n_act_orb
+        Umat(i,j)=0.D0
-      tt=list_act(t)
+        Tpotmat(i,j)=0.D0
-      indx+=1
+       end do
-      Tmat(ii,tt)= SXvector(indx)
+       Tpotmat(i,i)=1.D0
-      Tmat(tt,ii)=-SXvector(indx)
+      end do
-    end do
+
-  end do
+      indx=1
-  do i=1,n_core_orb
+      do i=1,n_core_orb
-    ii=list_core(i)
+       ii=list_core(i)
-    do a=1,n_virt_orb
+       do t=1,n_act_orb
-      aa=list_virt(a)
+        tt=list_act(t)
-      indx+=1
+        indx+=1
-      Tmat(ii,aa)= SXvector(indx)
+        Tmat(ii,tt)= SXvector(indx)
-      Tmat(aa,ii)=-SXvector(indx)
+        Tmat(tt,ii)=-SXvector(indx)
-    end do
+       end do
-  end do
+      end do
-  do t=1,n_act_orb
+      do i=1,n_core_orb
-    tt=list_act(t)
+       ii=list_core(i)
-    do a=1,n_virt_orb
+       do a=1,n_virt_orb
-      aa=list_virt(a)
+        aa=list_virt(a)
-      indx+=1
+        indx+=1
-      Tmat(tt,aa)= SXvector(indx)
+        Tmat(ii,aa)= SXvector(indx)
-      Tmat(aa,tt)=-SXvector(indx)
+        Tmat(aa,ii)=-SXvector(indx)
-    end do
+       end do
-  end do
+      end do
-  
+      do t=1,n_act_orb
-  ! Form the exponential
+       tt=list_act(t)
-
+       do a=1,n_virt_orb
-  Tpotmat(:,:)=0.D0
+        aa=list_virt(a)
-  Umat(:,:)   =0.D0
+        indx+=1
-  do i=1,mo_num
+        Tmat(tt,aa)= SXvector(indx)
-    Tpotmat(i,i)=1.D0
+        Tmat(aa,tt)=-SXvector(indx)
-    Umat(i,i)   =1.d0
+       end do
-  end do
+      end do
-  iter=0
+
-  converged=.false.
+      write(6,*) ' the T matrix '
-  do while (.not.converged)
+      do indx=1,nMonoEx
-    iter+=1
+       i=excit(1,indx)
-    f = 1.d0 / dble(iter)
+       j=excit(2,indx)
-    Tpotmat2(:,:) = Tpotmat(:,:) * f
+!       if (abs(Tmat(i,j)).gt.1.D0) then
-    call dgemm('N','N', mo_num,mo_num,mo_num,1.d0,                   &
+!        write(6,*) ' setting matrix element ',i,j,' of ',Tmat(i,j),' to ' &
-        Tpotmat2, size(Tpotmat2,1),                                  &
+!             , sign(1.D0,Tmat(i,j))
-        Tmat, size(Tmat,1), 0.d0,                                    &
+!        Tmat(i,j)=sign(1.D0,Tmat(i,j))
-        Tpotmat, size(Tpotmat,1))
+!        Tmat(j,i)=-Tmat(i,j)
-    Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
+!       end if
-    
+        if (abs(Tmat(i,j)).gt.1.D-9) write(6,9901) i,j,excit_class(indx),Tmat(i,j)
-    converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
+  9901 format('   ',i4,' -> ',i4,' (',A3,') : ',E14.6)
-  end do
+      end do
+
+      write(6,*) 
+      write(6,*) ' forming the matrix exponential '
+      write(6,*) 
+! form the exponential
+      iter=0
+      converged=.false.
+      do while (.not.converged)
+       iter+=1
+! add the next term
+       do i=1,mo_num
+        do j=1,mo_num
+         Umat(i,j)+=Tpotmat(i,j)
+        end do
+       end do
+! next power of T, we multiply Tpotmat with Tmat/iter
+       do i=1,mo_num
+        do j=1,mo_num
+         wrkline(j)=Tpotmat(i,j)/dble(iter)
+         Tpotmat(i,j)=0.D0
+        end do
+        do j=1,mo_num
+         do k=1,mo_num
+          Tpotmat(i,j)+=wrkline(k)*Tmat(k,j)
+         end do
+        end do
+       end do
+! Convergence test
+       sum=0.D0
+       do i=1,mo_num
+        do j=1,mo_num
+         sum+=abs(Tpotmat(i,j))
+        end do
+       end do
+       write(6,*) ' Iteration No ',iter,' Sum = ',sum
+       if (sum.lt.1.D-6) then
+        converged=.true.
+       end if
+       if (iter.ge.NItExpMax) then
+        stop ' no convergence '
+       end if
+      end do
+      write(6,*)
+      write(6,*) ' Converged ! '
+      write(6,*)
+
+      end subroutine get_orbrotmat
+
+BEGIN_PROVIDER [integer, NItExpMax]
+   NItExpMax=100
 END_PROVIDER
-  
 
 

src/casscf/tot_en.irp.f

@@ -42,6 +42,8 @@
        end do
      end do
    end do
+   write(6,*) ' e_one_all = ',e_one_all
+   write(6,*) ' e_two_all = ',e_two_all
    ecore    =nuclear_repulsion
    ecore_bis=nuclear_repulsion
    do i=1,n_core_orb
@@ -96,6 +98,24 @@
      end do
    end do
    
+   write(6,*) ' energy contributions '
+   write(6,*) '     core energy       = ',ecore,' using PQxx integrals '
+   write(6,*) '     core energy (bis) = ',ecore,' using PxxQ integrals '
+   write(6,*) '     1el  energy       = ',eone ,' using PQxx integrals '
+   write(6,*) '     1el  energy (bis) = ',eone ,' using PxxQ integrals '
+   write(6,*) '     2el  energy       = ',etwo    ,' using PQxx integrals '
+   write(6,*) '     2el  energy (bis) = ',etwo_bis,' using PxxQ integrals '
+   write(6,*) '     2el  energy (ter) = ',etwo_ter,' using tuvP integrals '
+   write(6,*) ' ----------------------------------------- '
+   write(6,*) '     sum of all        = ',eone+etwo+ecore
+   write(6,*)
+   write(6,*) ' nuclear     (qp)      = ',nuclear_repulsion
+   write(6,*) ' core energy (qp)      = ',core_energy
+   write(6,*) ' 1el  energy (qp)      = ',psi_energy_h_core(1)
+   write(6,*) ' 2el  energy (qp)      = ',psi_energy_two_e(1)
+   write(6,*) ' nuc + 1 + 2 (qp)      = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1)
+   write(6,*) ' <0|H|0>     (qp)      = ',psi_energy_with_nucl_rep(1)
+   
 END_PROVIDER