Optimized 1rdm

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,7 +55,6 @@
     end do
   end do
 
-  write(6,*) ' provided integrals (PQ|xx) '
 END_PROVIDER
 
 
@@ -116,7 +115,6 @@ 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
 
 
@@ -146,6 +144,5 @@ 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,7 +84,6 @@
       end do
     end do
   end do
-  write(6,*) ' transformed PQxx'
 
 END_PROVIDER
 
@@ -176,7 +175,6 @@ 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
 
@@ -267,7 +265,6 @@ 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,24 +12,32 @@ subroutine run
   implicit none
   double precision               :: energy_old, energy
   logical                        :: converged
+  integer                        :: iteration
   converged = .False.
 
   energy = 0.d0
-!  do while (.not.converged)
+  mo_label = "MCSCF"
-    N_det = 1
+  iteration = 1
-    TOUCH N_det psi_det psi_coef 
+  do while (.not.converged)
     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
+
-!  enddo
+    call write_time(6)
+    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,70 +1,19 @@
 use bitmasks
 
 BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
-  BEGIN_DOC
-  ! the first-order density matrix in the basis of the starting MOs
-  ! 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
   implicit none
-  integer                        :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
+  BEGIN_DOC
-  integer                        :: ierr
+  ! the first-order density matrix in the basis of the starting MOs.
-  integer(bit_kind)              :: det_mu(N_int,2)
+  ! matrix is state averaged.
-  integer(bit_kind)              :: det_mu_ex(N_int,2)
+  END_DOC
-  integer(bit_kind)              :: det_mu_ex1(N_int,2)
+  integer                        :: t,u
-  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)
   
-  write(6,*) ' providing density matrices D0 and P0 '
+  do u=1,n_act_orb
-  
-  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
-      ipart=list_act(t)
+      D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + &
-      do u=1,n_act_orb
+                  one_e_dm_mo_beta_average ( list_act(t), list_act(u) ) 
-        ihole=list_act(u)
+    enddo
-        ! apply E_tu
+  enddo
-        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
 
@@ -90,7 +39,9 @@ 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
    
-   write(6,*) ' providing density matrices D0 and P0 '
+  if (bavard) then
+    write(6,*) ' providing density matrix P0'
+  endif
 
    P0tuvx = 0.d0
    

src/casscf/det_manip.irp.f

@@ -31,6 +31,8 @@ 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
@@ -50,13 +52,6 @@ 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,32 +1,3 @@
-       subroutine driver_optorb
+subroutine driver_optorb
-         implicit none
+  implicit none
-         integer :: i,j
+end 
-
-         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,7 +6,6 @@ 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)]
@@ -87,9 +86,11 @@ BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
     norm_grad+=gradvec(indx)*gradvec(indx)
   end do
   norm_grad=sqrt(norm_grad)
-  write(6,*)
+  if (bavard) then
-  write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
+    write(6,*)
-  write(6,*)
+    write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
+    write(6,*)
+  endif
   
   
 END_PROVIDER
@@ -118,17 +119,11 @@ 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
@@ -176,9 +171,11 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
     norm_grad+=gradvec2(indx)*gradvec2(indx)
   end do
   norm_grad=sqrt(norm_grad)
-  write(6,*)
+  if (bavard) then
-  write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
+    write(6,*)
-  write(6,*)
+    write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
+    write(6,*)
+  endif
   
 END_PROVIDER
 

src/casscf/hessian.irp.f

@@ -14,8 +14,10 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
   character*3                    :: iexc,jexc
   real*8                         :: res
   
-  write(6,*) ' providing Hessian matrix hessmat '
+  if (bavard) then
-  write(6,*) '  nMonoEx = ',nMonoEx
+    write(6,*) ' providing Hessian matrix hessmat '
+    write(6,*) '  nMonoEx = ',nMonoEx
+  endif
   
   do indx=1,nMonoEx
     do jndx=1,nMonoEx
@@ -32,8 +34,6 @@ 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,8 +198,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
   real*8                         :: hessmat_iatb
   real*8                         :: hessmat_taub
   
-  write(6,*) ' providing Hessian matrix hessmat2 '
+  if (bavard) then
-  write(6,*) '  nMonoEx = ',nMonoEx
+    write(6,*) ' providing Hessian matrix hessmat2 '
+    write(6,*) '  nMonoEx = ',nMonoEx
+  endif
   
   indx=1
   do i=1,n_core_orb
@@ -214,7 +216,6 @@ 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
@@ -294,7 +295,6 @@ 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,8 +340,6 @@ 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
@@ -382,7 +380,6 @@ 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)
@@ -416,7 +413,6 @@ 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
@@ -457,7 +453,6 @@ 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)
@@ -495,7 +490,6 @@ 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)
@@ -552,7 +546,6 @@ 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,10 +14,12 @@
     occnum(list_act(i))=occ_act(n_act_orb-i+1)
   end do
 
-  write(6,*) ' occupation numbers '
+  if (bavard) then
-  do i=1,mo_num
+    write(6,*) ' occupation numbers '
-    write(6,*) i,occnum(i)
+    do i=1,mo_num
-  end do
+      write(6,*) i,occnum(i)
+    end do
+  endif
 
 END_PROVIDER
 
@@ -32,14 +34,12 @@ END_PROVIDER
  
  call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
  
- write(6,*) ' found occupation numbers as '
- do i=1,n_act_orb
-   write(6,*) i,occ_act(i)
- end do
- 
  if (bavard) then
-   !
+   write(6,*) ' found occupation numbers as '
-   
+   do i=1,n_act_orb
+     write(6,*) i,occ_act(i)
+   end do
+ 
    integer                        :: nmx
    real*8                         :: xmx
    do i=1,n_act_orb
@@ -152,7 +152,6 @@ 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
 
@@ -198,7 +197,6 @@ 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
 
 
@@ -226,148 +224,5 @@ 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,222 +1,178 @@
-! -*- F90 -*-
 BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
-      implicit none
+  implicit none
-      integer :: i,j
+  BEGIN_DOC
-      do i=1,nMonoEx+1
+  ! Single-excitation matrix
-       do j=1,nMonoEx+1
+  END_DOC
-        SXmatrix(i,j)=0.D0
+  
-       end do
+  integer                        :: i,j
-      end do
+  
-
+  do i=1,nMonoEx+1
-      do i=1,nMonoEx
+    do j=1,nMonoEx+1
-       SXmatrix(1,i+1)=gradvec2(i)
+      SXmatrix(i,j)=0.D0
-       SXmatrix(1+i,1)=gradvec2(i)
+    end do
-      end do
+  end do
-
+  
-      do i=1,nMonoEx
+  do i=1,nMonoEx
-       do j=1,nMonoEx
+    SXmatrix(1,i+1)=gradvec2(i)
-        SXmatrix(i+1,j+1)=hessmat2(i,j)
+    SXmatrix(1+i,1)=gradvec2(i)
-        SXmatrix(j+1,i+1)=hessmat2(i,j)
+  end do
-       end do
+  
-      end do
+  do i=1,nMonoEx
-
+    do j=1,nMonoEx
-   if (bavard) then
+      SXmatrix(i+1,j+1)=hessmat2(i,j)
-      do i=2,nMonoEx+1
+      SXmatrix(j+1,i+1)=hessmat2(i,j)
-       write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
+    end do
-      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)]
- END_PROVIDER
+  implicit none
+  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
-      integer :: ierr,matz,i
+  BEGIN_DOC
-      real*8 :: c0
+  ! Best eigenvector of the single-excitation matrix
+  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
 
-      call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
+  c0=SXeigenvec(1,best_vector)
-      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)
 
-integer :: best_vector
+  do i=1,nMonoEx+1
-real*8 :: best_overlap
+    SXvector(i)=SXeigenvec(i,best_vector)/c0
-      best_overlap=0.D0
+  end do
-      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
-      integer :: i,j,ialph
+  BEGIN_DOC
-
+  ! Updated orbitals
-! form the exponential of the Orbital rotations
+  END_DOC
-      call get_orbrotmat
+  integer                        :: i,j,ialph
-! form the new orbitals
+  
-      do i=1,ao_num
+  call dgemm('N','T', ao_num,mo_num,mo_num,1.d0,                     &
-       do j=1,mo_num
+      NatOrbsFCI, size(NatOrbsFCI,1),                                &
-        NewOrbs(i,j)=0.D0
+      Umat, size(Umat,1), 0.d0,                                      &
-       end do
+      NewOrbs, size(NewOrbs,1))
-      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, Tpotmat, (mo_num,mo_num) ]
+BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
-&BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
+  implicit none
-&BEGIN_PROVIDER [real*8, wrkline, (mo_num) ]
+  BEGIN_DOC
-&BEGIN_PROVIDER [real*8, Tmat, (mo_num,mo_num) ]
+  ! Orbital rotation matrix
-END_PROVIDER
+  END_DOC
-
+  integer                        :: i,j,indx,k,iter,t,a,ii,tt,aa
-      subroutine get_orbrotmat
+  logical                        :: converged
-      implicit none
+  
-      integer :: i,j,indx,k,iter,t,a,ii,tt,aa
+  real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num) 
-      real*8 :: sum
+  real*8 :: Tmat(mo_num,mo_num) 
-      logical :: converged
+  real*8 :: f
-
+  
-
+  ! the orbital rotation matrix T
-! the orbital rotation matrix T
+  Tmat(:,:)=0.D0
-      do i=1,mo_num
+  indx=1
-       do j=1,mo_num
+  do i=1,n_core_orb
-        Tmat(i,j)=0.D0
+    ii=list_core(i)
-        Umat(i,j)=0.D0
+    do t=1,n_act_orb
-        Tpotmat(i,j)=0.D0
+      tt=list_act(t)
-       end do
+      indx+=1
-       Tpotmat(i,i)=1.D0
+      Tmat(ii,tt)= SXvector(indx)
-      end do
+      Tmat(tt,ii)=-SXvector(indx)
-
+    end do
-      indx=1
+  end do
-      do i=1,n_core_orb
+  do i=1,n_core_orb
-       ii=list_core(i)
+    ii=list_core(i)
-       do t=1,n_act_orb
+    do a=1,n_virt_orb
-        tt=list_act(t)
+      aa=list_virt(a)
-        indx+=1
+      indx+=1
-        Tmat(ii,tt)= SXvector(indx)
+      Tmat(ii,aa)= SXvector(indx)
-        Tmat(tt,ii)=-SXvector(indx)
+      Tmat(aa,ii)=-SXvector(indx)
-       end do
+    end do
-      end do
+  end do
-      do i=1,n_core_orb
+  do t=1,n_act_orb
-       ii=list_core(i)
+    tt=list_act(t)
-       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(ii,aa)= SXvector(indx)
+      Tmat(tt,aa)= SXvector(indx)
-        Tmat(aa,ii)=-SXvector(indx)
+      Tmat(aa,tt)=-SXvector(indx)
-       end do
+    end do
-      end do
+  end do
-      do t=1,n_act_orb
+  
-       tt=list_act(t)
+  ! Form the exponential
-       do a=1,n_virt_orb
+
-        aa=list_virt(a)
+  Tpotmat(:,:)=0.D0
-        indx+=1
+  Umat(:,:)   =0.D0
-        Tmat(tt,aa)= SXvector(indx)
+  do i=1,mo_num
-        Tmat(aa,tt)=-SXvector(indx)
+    Tpotmat(i,i)=1.D0
-       end do
+    Umat(i,i)   =1.d0
-      end do
+  end do
-
+  iter=0
-      write(6,*) ' the T matrix '
+  converged=.false.
-      do indx=1,nMonoEx
+  do while (.not.converged)
-       i=excit(1,indx)
+    iter+=1
-       j=excit(2,indx)
+    f = 1.d0 / dble(iter)
-!       if (abs(Tmat(i,j)).gt.1.D0) then
+    Tpotmat2(:,:) = Tpotmat(:,:) * f
-!        write(6,*) ' setting matrix element ',i,j,' of ',Tmat(i,j),' to ' &
+    call dgemm('N','N', mo_num,mo_num,mo_num,1.d0,                   &
-!             , sign(1.D0,Tmat(i,j))
+        Tpotmat2, size(Tpotmat2,1),                                  &
-!        Tmat(i,j)=sign(1.D0,Tmat(i,j))
+        Tmat, size(Tmat,1), 0.d0,                                    &
-!        Tmat(j,i)=-Tmat(i,j)
+        Tpotmat, size(Tpotmat,1))
-!       end if
+    Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
-        if (abs(Tmat(i,j)).gt.1.D-9) write(6,9901) i,j,excit_class(indx),Tmat(i,j)
+    
-  9901 format('   ',i4,' -> ',i4,' (',A3,') : ',E14.6)
+    converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
-      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,8 +42,6 @@
        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
@@ -98,24 +96,6 @@
      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