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=.false.
+!     bavard=.true.
+     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
-    SOFT_TOUCH mo_coef N_det
+    N_det = N_det/2
+    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
@@ -43,125 +45,25 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
     write(6,*) ' providing density matrix P0'
   endif
 
-   P0tuvx = 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)
-       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
-     do v=1,n_act_orb
-       do u=1,n_act_orb
-         do t=1,n_act_orb
-           P0tuvx(t,u,v,x)*=0.5D0/dble(N_states)
-         end do
-       end do
-     end do
-   end do
+  P0tuvx= 0.d0
+  do istate=1,N_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
+       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
+  enddo
 
 END_PROVIDER

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,mo_num
-    do j=1,mo_num
-      do k=1,mo_num
-        do l=1,mo_num
-         twodm = coussin_peter_two_rdm_mo(i,j,k,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/README.rst

@@ -2,5 +2,7 @@
 two_body_rdm
 ============
 
-Contains the two rdms (aa,bb,ab) stored as plain arrays 
+Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as
+maps, with pysicists notation, consistent with the two-electron integrals in the
+MO basis.
 

src/two_body_rdm/all_2rdm_routines.irp.f

@@ -1,5 +1,4 @@
-
- subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze)
+subroutine all_two_rdm_dm_nstates_openmp(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze)
    use bitmasks
    implicit none
    BEGIN_DOC
@@ -36,10 +35,10 @@
      call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
    enddo
    
- end
+end
  
  
- subroutine all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
+subroutine all_two_rdm_dm_nstates_openmp_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
    use bitmasks
    implicit none
    BEGIN_DOC
@@ -69,9 +68,9 @@
        case default
        call all_two_rdm_dm_nstates_openmp_work_N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
    end select
- end
+end
  
-BEGIN_TEMPLATE
+ BEGIN_TEMPLATE
  
 subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
    use bitmasks
@@ -121,8 +120,8 @@ subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,b
    !!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
        !    !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,          &
        !    !$OMP          psi_bilinear_matrix_columns,              &
- !    !$OMP          psi_det_alpha_unique, psi_det_beta_unique,      &
- !    !$OMP          n_det_alpha_unique, n_det_beta_unique, N_int,   &
+       !    !$OMP          psi_det_alpha_unique, psi_det_beta_unique,&
+       !    !$OMP          n_det_alpha_unique, n_det_beta_unique, N_int,&
        !    !$OMP          psi_bilinear_matrix_transp_rows,          &
        !    !$OMP          psi_bilinear_matrix_transp_columns,       &
        !    !$OMP          psi_bilinear_matrix_transp_order, N_st,   &
@@ -131,7 +130,7 @@ subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,b
        !    !$OMP          psi_bilinear_matrix_transp_rows_loc,      &
        !    !$OMP          istart, iend, istep, irp_here, v_t, s_t,  &
        !    !$OMP          ishift, idx0, u_t, maxab)                 &
- !    !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,     &
+       !    !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,&
        !    !$OMP          lcol, lrow, l_a, l_b,                     &
        !    !$OMP          buffer, doubles, n_doubles,               &
        !    !$OMP          tmp_det2, idx, l, kcol_prev,              &
@@ -222,9 +221,9 @@ subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,b
      enddo
      
    enddo
-! !$OMP END DO
+   ! !$OMP END DO
    
-! !$OMP DO SCHEDULE(dynamic,64)
+   ! !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
      
      
@@ -431,13 +430,13 @@ subroutine all_two_rdm_dm_nstates_openmp_work_$N_int(big_array_aa,big_array_bb,b
    
 end
  
-SUBST [ N_int ]
+ SUBST [ N_int ]
  
-1;;
-2;;
-3;;
-4;;
-N_int;;
+ 1;;
+ 2;;
+ 3;;
+ 4;;
+ N_int;;
  
-END_TEMPLATE
+ END_TEMPLATE
  

src/two_body_rdm/two_rdm.irp.f

@@ -1,25 +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
- 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,:) = 0.5d0 * (two_rdm_alpha_beta_mo(i,j,k,l,:) + two_rdm_alpha_beta_mo(i,j,k,l,:)) & 
-                                         +          two_rdm_alpha_alpha_mo(i,j,k,l,:)                                  &
-                                         +          two_rdm_beta_beta_mo(i,j,k,l,:)         
-    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)]
@@ -45,7 +23,7 @@
   call wall_time(cpu_1)
   print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0)
   
- END_PROVIDER 
+END_PROVIDER
 
 
  BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
@@ -80,5 +58,5 @@
   call wall_time(cpu_1)
   print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0)
   
- END_PROVIDER 
+END_PROVIDER