Merge branch 'casscf' of github.com:QuantumPackage/qp2 into casscf

src/bitmask/core_inact_act_virt.irp.f

@@ -141,6 +141,10 @@ END_PROVIDER
  n_act_orb_tmp = 0
  n_virt_orb_tmp = 0
  n_del_orb_tmp = 0
+ core_bitmask = 0_bit_kind
+ inact_bitmask = 0_bit_kind
+ act_bitmask = 0_bit_kind
+ virt_bitmask = 0_bit_kind
  do i = 1, mo_num
   if(mo_class(i) == 'Core')then
    n_core_orb_tmp += 1

src/casscf/bavard.irp.f

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

src/casscf/casscf.irp.f

@@ -4,7 +4,8 @@ program casscf
 ! TODO : Put the documentation of the program here
   END_DOC
   no_vvvv_integrals = .True.
-  SOFT_TOUCH no_vvvv_integrals
+  pt2_max = 0.02
+  SOFT_TOUCH no_vvvv_integrals pt2_max
   call run
 end
 
@@ -19,8 +20,7 @@ subroutine run
   mo_label = "MCSCF"
   iteration = 1
   do while (.not.converged)
-    call run_cipsi
+    call run_stochastic_cipsi
-
     energy_old = energy
     energy = eone+etwo+ecore
 
@@ -30,14 +30,19 @@ subroutine run
     call write_double(6,energy_improvement, 'Predicted energy improvement')
 
     converged = dabs(energy_improvement) < thresh_scf
+    pt2_max = dabs(energy_improvement / pt2_relative_error)
 
     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
+    N_det = N_det/2
-    SOFT_TOUCH mo_coef N_det
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    read_wf = .True.
+    FREE mo_integrals_map mo_two_e_integrals_in_map
+    SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef
+
   enddo
 
 end

src/casscf/densities.irp.f

@@ -29,7 +29,9 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
    !
    END_DOC
    implicit none
-   integer                        :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
+   integer                        :: t,u,v,x
+   integer                        :: tt,uu,vv,xx
+   integer                        :: mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
    integer                        :: ierr
    real*8                         :: phase1,phase11,phase12,phase2,phase21,phase22
    integer                        :: nu1,nu2,nu11,nu12,nu21,nu22
@@ -44,121 +46,21 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
   endif
 
   P0tuvx= 0.d0
-   
+  do istate=1,N_states
-   ! 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)
-       do u=1,n_act_orb
-         ihole=list_act(u)
-         ! 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
-             ! and we fill P0_tvvu
-             do v=1,n_act_orb
-               P0tuvx(t,v,v,u)-=term
-             end do
-           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
-             do v=1,n_act_orb
-               P0tuvx(t,v,v,u)-=term
-             end do
-           end do
-         end if
-       end do
-     end do
-   end do
-   ! now we do the double excitation E_tu E_vx |0>
-   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 v=1,n_act_orb
-       ipart=list_act(v)
-       do x=1,n_act_orb
-         ihole=list_act(x)
-         ! apply E_vx
-         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)
-         ! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0>
-         if (ierr1.eq.1) then
-           do t=1,n_act_orb
-             jpart=list_act(t)
-             do u=1,n_act_orb
-               jhole=list_act(u)
-               call det_copy(det_mu_ex1,det_mu_ex11,N_int)
-               call det_copy(det_mu_ex1,det_mu_ex12,N_int)
-               call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11&
-                   ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12)
-               if (nu11.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)&
-                       *phase11*phase1
-                 end do
-               end if
-               if (nu12.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)&
-                       *phase12*phase1
-                 end do
-               end if
-             end do
-           end do
-         end if
-         
-         ! we apply E_tu to the second resultant determinant
-         if (ierr2.eq.1) then
-           do t=1,n_act_orb
-             jpart=list_act(t)
-             do u=1,n_act_orb
-               jhole=list_act(u)
-               call det_copy(det_mu_ex2,det_mu_ex21,N_int)
-               call det_copy(det_mu_ex2,det_mu_ex22,N_int)
-               call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21&
-                   ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22)
-               if (nu21.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)&
-                       *phase21*phase2
-                 end do
-               end if
-               if (nu22.ne.-1) then
-                 do istate=1,n_states
-                   P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)&
-                       *phase22*phase2
-                 end do
-               end if
-             end do
-           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
+    xx = list_act(x)
     do v = 1, n_act_orb
+     vv = list_act(v)
      do u = 1, n_act_orb
+      uu = list_act(u)
       do t = 1, n_act_orb
-           P0tuvx(t,u,v,x)*=0.5D0/dble(N_states)
+       tt = list_act(t)
+       P0tuvx(t,u,v,x) =                                   &
+         state_average_weight(istate) *                    &
+         ( two_rdm_alpha_beta_mo (tt,uu,vv,xx,istate) +    &
+           two_rdm_alpha_alpha_mo(tt,uu,vv,xx,istate) +    &
+           two_rdm_beta_beta_mo  (tt,uu,vv,xx,istate) )
+      enddo
      enddo 
     enddo
    enddo

src/casscf/neworbs.irp.f

@@ -25,7 +25,7 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
   end do
   
   if (bavard) then
-    do i=2,nMonoEx+1
+    do i=2,nMonoEx
       write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
     end do
   end if
@@ -77,14 +77,14 @@ END_PROVIDER
   
   energy_improvement = SXeigenval(best_vector)
   
+  c0=SXeigenvec(1,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)
-
   do i=1,nMonoEx+1
     SXvector(i)=SXeigenvec(i,best_vector)/c0
   end do

src/casscf/test_two_rdm.irp.f

@@ -1,30 +0,0 @@
-program print_two_rdm
-  implicit none
-  integer                        :: i,j,k,l
-  read_wf = .True.
-  TOUCH read_wf
-
-  double precision, parameter    :: thr = 1.d-15
-
-  double precision :: accu,twodm
-  accu = 0.d0
-  do i=1,n_act_orb
-    do j=1,n_act_orb
-      do k=1,n_act_orb
-        do l=1,n_act_orb
-         twodm = coussin_peter_two_rdm_mo(list_act(i),list_act(j),list_act(k),list_act(l),1)
-            if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then
-             print*,''
-             print*,'sum'
-             write(*,'(3X,4(I2,X),3(F16.13,X))'),  i, j, k, l, twodm,P0tuvx(i,j,k,l),dabs(twodm - P0tuvx(i,j,k,l))
-             print*,''
-            endif
-            accu += dabs(twodm - P0tuvx(i,j,k,l))
-        enddo
-      enddo
-    enddo
-  enddo
-  print*,'accu = ',accu
-  print*,'<accu> ',accu / dble(mo_num**4)
-
-end

src/cipsi/cipsi.irp.f

@@ -13,6 +13,7 @@ subroutine run_cipsi
   rss = memory_of_double(N_states)*4.d0
   call check_mem(rss,irp_here)
 
+  N_iter = 1
   allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))
 
   double precision               :: hf_energy_ref

src/cipsi/pt2_stoch_routines.irp.f

@@ -135,7 +135,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
     PROVIDE psi_occ_pattern_hii det_to_occ_pattern
   endif
 
-  if (N_det < max(4,N_states)) then
+  if (N_det <= max(4,N_states)) then
     pt2=0.d0
     variance=0.d0
     norm=0.d0
@@ -719,6 +719,15 @@ END_PROVIDER
 
   double precision :: rss
   double precision, external :: memory_of_double, memory_of_int
+  if (N_det_generators == 1) then
+    pt2_w = 1.d0
+    pt2_cw = 1.d0
+    pt2_W_T = 1.d0
+    pt2_u_0 = 1.d0
+    pt2_n_0 = 1
+    return
+  endif
+
   rss = memory_of_double(2*N_det_generators+1)
   call check_mem(rss,irp_here)
 
@@ -754,7 +763,7 @@ END_PROVIDER
     end if
     pt2_n_0(1) += 1
     if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then
-      stop "teeth building failed"
+      print *, "teeth building failed"
     end if
   end do
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

src/cipsi/stochastic_cipsi.irp.f

@@ -12,6 +12,7 @@ subroutine run_stochastic_cipsi
   double precision, external :: memory_of_double
   PROVIDE H_apply_buffer_allocated N_generators_bitmask
 
+  N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 

src/generators_cas/generators.irp.f

@@ -55,6 +55,7 @@ END_PROVIDER
       nongen(inongen) = i
     endif
   enddo
+  ASSERT (m == N_det_generators)
 
   psi_det_sorted_gen(:,:,:N_det_generators) = psi_det_generators(:,:,:N_det_generators)
   psi_coef_sorted_gen(:N_det_generators, :) = psi_coef_generators(:N_det_generators, :)

src/two_body_rdm/two_rdm.irp.f

@@ -1,27 +1,3 @@
-BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
-  implicit none
-  BEGIN_DOC
-  !  coussin_peter_two_rdm_mo(i,j,k,l) = the two rdm that peter wants for his CASSCF
-  END_DOC
-  integer                        :: i,j,k,l, istate
-  do istate = 1,N_states
-    do l = 1, mo_num
-      do k = 1, mo_num
-        do j = 1, mo_num
-          do i = 1, mo_num
-            coussin_peter_two_rdm_mo  (i,j,k,l,istate) =             &
-                two_rdm_alpha_beta_mo (i,j,k,l,istate) +             &
-                two_rdm_alpha_alpha_mo(i,j,k,l,istate) +             &
-                two_rdm_beta_beta_mo  (i,j,k,l,istate)
-          enddo
-        enddo
-      enddo
-    enddo
-  enddo
-  
-END_PROVIDER
-
-
  BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
 &BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
 &BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]