Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

external/ezfio

@@ -1 +1 @@
-Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c
+Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93

external/irpf90

@@ -1 +1 @@
-Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271
+Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102

src/casscf_tc_bi/grad_dm.irp.f

@@ -50,11 +50,13 @@
     end do
   end do
   
+! print*,'DM  grad'
   do t=1,n_act_orb
     tt=list_act(t)
     do a=1,n_virt_orb
       aa=list_virt(a)
       indx = mat_idx_a_v(t,a) 
+!      print*,indx,t,a
       call gradvec_tc_ta(tt,aa,res_l, res_r)
       do fff = 0,3
        gradvec_tc_l(fff,indx)=res_l(fff)
@@ -79,7 +81,28 @@ subroutine gradvec_tc_ia(i,a,res_l, res_r)
  res_r = 0.d0
  res_l(1) = -2 * mo_bi_ortho_tc_one_e(a,i)
  res_r(1) = -2 * mo_bi_ortho_tc_one_e(i,a)
- 
+ integer :: j,t,r,jj,tt,rr
+ do jj = 1, n_core_inact_orb
+  j = list_core_inact(jj)
+  res_r(2) += -2.d0 * ( 2.d0 * mo_bi_ortho_tc_two_e(j,i,j,a) - mo_bi_ortho_tc_two_e(i,j,j,a))
+  res_l(2) -= -2.d0 * ( 2.d0 * mo_bi_ortho_tc_two_e(j,a,j,i) - mo_bi_ortho_tc_two_e(j,a,i,j))
+ enddo
+ do tt = 1, n_act_orb
+  t = list_act(tt)
+  do rr = 1, n_act_orb
+   r = list_act(rr)
+   res_r(2) += -0.5d0 * (                                                                                    &
+   tc_transition_matrix_mo(r,t,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(r,i,t,a) - mo_bi_ortho_tc_two_e(i,r,t,a))  & 
+  +tc_transition_matrix_mo(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(t,i,r,a) - mo_bi_ortho_tc_two_e(i,t,r,a))  &
+   )
+   res_l(2) -= -0.5d0 * (                                                                                    &
+   tc_transition_matrix_mo(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(t,a,r,i) - mo_bi_ortho_tc_two_e(t,a,i,r))  &
+  +tc_transition_matrix_mo(r,t,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(r,a,t,i) - mo_bi_ortho_tc_two_e(r,a,i,t))  & 
+   )
+  enddo
+ enddo
+ res_r(0) = res_r(1) + res_r(2) + res_r(3)
+ res_l(0) = res_l(1) + res_l(2) + res_l(3)
 end
 
 subroutine gradvec_tc_it(i,t,res_l, res_r)
@@ -93,19 +116,51 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
  END_DOC
  integer, intent(in) :: i,t
  double precision, intent(out) :: res_l(0:3),res_r(0:3)
- integer :: rr,r,ss,s,m,mm
+ integer :: rr,r,j,jj,u,uu,v,vv
- double precision :: dm
  res_r = 0.d0
  res_l = 0.d0
  res_r(1)  += -2.d0 * mo_bi_ortho_tc_one_e(i,t) 
- res_l(1)  +=  2.D0 * mo_bi_ortho_tc_one_e(t,i) 
+ res_l(1)  -= -2.D0 * mo_bi_ortho_tc_one_e(t,i) 
 
  do rr = 1, n_act_orb
   r = list_act(rr)
   res_r(1) +=  mo_bi_ortho_tc_one_e(i,r) * tc_transition_matrix_mo(t,r,1,1)   
-  res_l(1) += -mo_bi_ortho_tc_one_e(r,i) * tc_transition_matrix_mo(r,t,1,1)
+  res_l(1) -=  mo_bi_ortho_tc_one_e(r,i) * tc_transition_matrix_mo(r,t,1,1)
  enddo
 
+ do jj = 1, n_core_inact_orb
+  j = list_core_inact(jj)
+  res_r(2) += -2.d0 * (2d0 * mo_bi_ortho_tc_two_e(i,j,t,j) - mo_bi_ortho_tc_two_e(j,i,t,j)) 
+  res_l(2) -= -2.d0 * (2d0 * mo_bi_ortho_tc_two_e(t,j,i,j) - mo_bi_ortho_tc_two_e(t,j,j,i)) 
+  do rr = 1, n_act_orb
+   r = list_act(rr)
+   res_r(2) += tc_transition_matrix_mo(t,r,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(i,j,r,j) - mo_bi_ortho_tc_two_e(i,j,j,r)) 
+   res_l(2) -= tc_transition_matrix_mo(r,t,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(r,j,i,j) - mo_bi_ortho_tc_two_e(j,r,j,i)) 
+  enddo
+ enddo
+ do rr = 1, n_act_orb
+  r = list_act(rr)
+  do uu = 1, n_act_orb
+   u = list_act(uu)
+   res_r(2) += -0.5d0 * (  & 
+    tc_transition_matrix_mo(u,r,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(u,i,r,t) - mo_bi_ortho_tc_two_e(u,i,t,r))  &
+  + tc_transition_matrix_mo(r,u,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(i,r,t,u) - mo_bi_ortho_tc_two_e(i,r,u,t))  &
+    )
+   res_l(2) -= -0.5d0 * (  & 
+    tc_transition_matrix_mo(r,u,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(r,t,u,i) - mo_bi_ortho_tc_two_e(t,r,u,i))  &
+  + tc_transition_matrix_mo(u,r,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(t,u,i,r) - mo_bi_ortho_tc_two_e(u,t,i,r))  &
+    )
+   do vv = 1, n_act_orb
+    v = list_act(vv)
+    res_r(2) +=  0.5d0 * (  & 
+    mo_bi_ortho_tc_two_e(i,r,v,u) * tc_two_rdm(t,r,v,u) + mo_bi_ortho_tc_two_e(r,i,v,u) * tc_two_rdm(r,t,v,u) )
+    res_l(2) -=  0.5d0 * (  & 
+    mo_bi_ortho_tc_two_e(v,u,i,r) * tc_two_rdm(v,u,t,r) + mo_bi_ortho_tc_two_e(v,u,r,i) * tc_two_rdm(v,u,r,t) )
+   enddo
+  enddo
+ enddo
+ res_r(0) = res_r(1) + res_r(2) + res_r(3)
+ res_l(0) = res_l(1) + res_l(2) + res_l(3)
 end
 
 subroutine gradvec_tc_ta(t,a,res_l, res_r)
@@ -119,14 +174,90 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r)
  END_DOC
  integer, intent(in) :: t,a
  double precision, intent(out) :: res_l(0:3),res_r(0:3)
- integer :: rr,r,m
+ integer :: rr,r,j,jj,u,uu,v,vv
- double precision :: dm
+ double precision :: res_r_inact_test, res_r_act_test
+ double precision :: res_l_inact_test, res_l_act_test
  res_r = 0.d0
  res_l = 0.d0
  do rr = 1, n_act_orb
   r = list_act(rr) 
   res_l(1) +=  mo_bi_ortho_tc_one_e(a,r) * tc_transition_matrix_mo(t,r,1,1)
-  res_r(1) += -mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1)
+  res_r(1) -=  mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1)
  enddo
 
+ res_r_inact_test = 0.d0
+ res_l_inact_test = 0.d0
+ do jj = 1, n_core_inact_orb
+  j = list_core_inact(jj)
+  do rr = 1, n_act_orb
+   r = list_act(rr)
+   res_r_inact_test += -tc_transition_matrix_mo(r,t,1,1) * & 
+   (2.d0 * mo_bi_ortho_tc_two_e(r,j,a,j) - mo_bi_ortho_tc_two_e(r,j,j,a))
+   res_l_inact_test += -tc_transition_matrix_mo(t,r,1,1) * & 
+   (2.d0 * mo_bi_ortho_tc_two_e(a,j,r,j) - mo_bi_ortho_tc_two_e(j,a,r,j))
+  enddo
+ enddo
+ res_r_act_test = 0.d0
+ res_l_act_test = 0.d0
+ do rr = 1, n_act_orb
+  r = list_act(rr)
+  do vv = 1, n_act_orb
+   v = list_act(vv)
+   do uu = 1, n_act_orb
+    u = list_act(uu)
+    res_r_act_test += - (mo_bi_ortho_tc_two_e(v,r,u,a) * tc_two_rdm(r,v,t,u) & 
+                        +mo_bi_ortho_tc_two_e(v,r,a,u) * tc_two_rdm(r,v,u,t))
+    res_l_act_test += - (mo_bi_ortho_tc_two_e(u,a,v,r) * tc_two_rdm(t,u,r,v) & 
+                        +mo_bi_ortho_tc_two_e(a,u,v,r) * tc_two_rdm(u,t,r,v))
+   enddo
+  enddo
+ enddo
+ res_r_act_test *= 0.5d0
+ res_l_act_test *= 0.5d0
+ res_r(2) = res_r_inact_test + res_r_act_test
+ res_l(2) = res_l_inact_test + res_l_act_test
+
+ integer :: m,x,y
+ double precision :: res_r_inact, res_r_act
+ if(.False.)then
+   ! test quantities
+   res_r_inact = 0.d0
+   res_r_act   = 0.d0
+   do m = 1, mo_num
+    do x = 1, mo_num
+     do jj = 1, n_core_inact_orb
+      j = list_core_inact(jj) 
+      res_r_inact += 0.5d0 * mo_bi_ortho_tc_two_e(t,j,m,x) * tc_two_rdm(a,j,m,x) & 
+                  -0.5d0 * mo_bi_ortho_tc_two_e(m,j,a,x) * tc_two_rdm(m,j,t,x) & 
+                  +0.5d0 * mo_bi_ortho_tc_two_e(j,t,m,x) * tc_two_rdm(j,a,m,x) & 
+                  -0.5d0 * mo_bi_ortho_tc_two_e(x,j,m,a) * tc_two_rdm(x,j,m,t)
+     enddo
+     do rr = 1, n_act_orb
+      r = list_act(rr) 
+      res_r_act   += 0.5d0 * mo_bi_ortho_tc_two_e(t,r,m,x) * tc_two_rdm(a,r,m,x) & 
+                  -0.5d0 * mo_bi_ortho_tc_two_e(m,r,a,x) * tc_two_rdm(m,r,t,x) & 
+                  +0.5d0 * mo_bi_ortho_tc_two_e(r,t,m,x) * tc_two_rdm(r,a,m,x) & 
+                  -0.5d0 * mo_bi_ortho_tc_two_e(x,r,m,a) * tc_two_rdm(x,r,m,t)
+     enddo
+    enddo
+   enddo
+   if(dabs(res_r_inact).gt.1.d-12)then
+    if(dabs(res_r_inact_test - res_r_inact).gt.1.d-10)then
+     print*,'inact'
+     print*,'t,a',t,a
+     print*,res_r_inact_test , res_r_inact, dabs(res_r_inact_test - res_r_inact)
+    endif
+   endif
+   if(dabs(res_r_act).gt.1.d-12)then
+    if(dabs(res_r_act_test - res_r_act).gt.1.d-10)then
+     print*,'act'
+     print*,'t,a',t,a
+     print*,res_r_act_test , res_r_act, dabs(res_r_act_test - res_r_act)
+    endif
+   endif
+ endif
+
+ res_r(0) = res_r(1) + res_r(2) + res_r(3)
+ res_l(0) = res_l(1) + res_l(2) + res_l(3)
+ 
 end

src/casscf_tc_bi/grad_old.irp.f

@@ -39,11 +39,13 @@
   enddo
  enddo
 
+! print*,'old grad'
  do tt = 1, n_act_orb
   ihole = list_act(tt)
   do aa = 1, n_virt_orb
    ipart = list_virt(aa)
    indx = mat_idx_a_v(tt,aa) 
+!   print*,indx,tt,aa
    call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
    do ll = 0, 3
     gradvec_detail_left_old (ll,indx)=res_l(ll)

src/mo_optimization/EZFIO.cfg

@@ -2,7 +2,7 @@
 type: character*(32)
 doc: Define the kind of hessian for the orbital optimization full : full hessian, diag : diagonal hessian, none : no hessian
 interface: ezfio,provider,ocaml
-default: full
+default: diag
 
 [n_det_max_opt]
 type: integer
@@ -14,7 +14,7 @@ default: 200000
 type: integer
 doc: Maximal number of iterations for the orbital optimization
 interface: ezfio,provider,ocaml
-default: 20
+default: 10
 
 [thresh_opt_max_elem_grad]
 type: double precision

src/mo_optimization/cipsi_orb_opt.irp.f

@@ -15,7 +15,7 @@ subroutine run_optimization
   logical :: not_converged
   character (len=100) :: filename
 
-  PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
+  PROVIDE psi_det psi_coef mo_two_e_integrals_in_map ao_pseudo_integrals
 
   not_converged = .True.
   nb_iter = 0

src/mo_optimization/state_average_energy.irp.f

@@ -40,16 +40,23 @@ subroutine state_average_energy(energy)
   double precision :: mono_e, bi_e
   integer :: i,j,k,l
 
+  energy = nuclear_repulsion
   ! mono electronic part
+  !$OMP PARALLEL DEFAULT(NONE) PRIVATE(i,j,k,l,mono_e, bi_e) &
+  !$OMP SHARED(mo_num, mo_integrals_map, two_e_dm_mo, one_e_dm_mo, energy, &
+  !$OMP       mo_one_e_integrals)
   mono_e = 0d0
+  !$OMP DO
   do j = 1, mo_num
     do i = 1, mo_num
       mono_e =  mono_e + mo_one_e_integrals(i,j) * one_e_dm_mo(i,j)
     enddo
   enddo
+  !$OMP END DO NOWAIT
 
   ! bi electronic part
   bi_e = 0d0
+  !$OMP DO
   do l = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -59,9 +66,13 @@ subroutine state_average_energy(energy)
       enddo
     enddo
   enddo
+  !$OMP END DO
 
   ! State average energy
-  energy = mono_e + 0.5d0 * bi_e + nuclear_repulsion
+  !$OMP CRITICAL
+  energy = energy + mono_e + 0.5d0 * bi_e
+  !$OMP END CRITICAL
+  !$OMP END PARALLEL
 
   ! Check
   !call print_energy_components

src/tc_bi_ortho/two_rdm_naive.irp.f

@@ -1,7 +1,7 @@
-BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num, mo_num)]
  implicit none
  BEGIN_DOC
- ! tc_two_rdm(p,s,q,r) = <Phi| a^dager_p
+ ! tc_two_rdm_chemist(p,s,q,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = CHEMIST NOTATION
  END_DOC
  integer :: i,j,istate,m,mm,nn
  integer                        :: exc(0:2,2,2)
@@ -13,7 +13,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
  other_spin(1) = 2
  other_spin(2) = 1
  allocate(occ(N_int*bit_kind_size,2))
- tc_two_rdm = 0.d0
+ tc_two_rdm_chemist = 0.d0
 
  do i = 1, N_det ! psi_left 
   do j = 1, N_det ! psi_right 
@@ -28,7 +28,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
      contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
     enddo
     if(degree == 2)then
-     call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+     call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
     else if(degree==1)then
      ! occupation of the determinant psi_det(j)
      call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 
@@ -39,7 +39,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
       m = occ(mm,s2)
       h2 = m 
       p2 = m 
-      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
      enddo
      ! run over the electrons of same spin than the excitation
      s2 = s1
@@ -47,7 +47,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
       m = occ(mm,s2)
       h2 = m 
       p2 = m 
-      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
      enddo
     endif
    else if(degree == 0)then
@@ -68,7 +68,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
        m = occ(mm,s2)
        h2 = m 
        p2 = m 
-       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
       enddo
       ! run over the couple of alpha-alpha electrons 
       s2 = s1
@@ -77,7 +77,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
        h2 = m 
        p2 = m 
        if(h2.le.h1)cycle
-       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
       enddo
     enddo
     s1 = 2
@@ -91,7 +91,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
        h2 = m 
        p2 = m 
        if(h2.le.h1)cycle
-       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
       enddo
     enddo
    endif
@@ -122,3 +122,22 @@ subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib)
  endif
 
 end
+
+
+BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! tc_two_rdm(p,q,s,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = PHYSICIST NOTATION
+ END_DOC
+ integer :: p,q,r,s
+ do r = 1, mo_num
+  do q = 1, mo_num
+   do s = 1, mo_num
+    do p = 1, mo_num
+    tc_two_rdm(p,q,s,r) = tc_two_rdm_chemist(p,s,q,r) 
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER 

src/two_body_rdm/state_av_act_2rdm.irp.f

@@ -17,12 +17,12 @@
  state_weights = state_average_weight
  integer :: ispin
  ! condition for alpha/beta spin
- print*,''
+! print*,''
- print*,''
+! print*,''
- print*,''
+! print*,''
- print*,'providint state_av_act_2_rdm_ab_mo '
+! print*,'Providing state_av_act_2_rdm_ab_mo '
  ispin = 3 
- print*,'ispin = ',ispin
+! print*,'ispin = ',ispin
  state_av_act_2_rdm_ab_mo = 0.d0
  call wall_time(wall_1)
  double precision :: wall_1, wall_2

src/two_body_rdm/state_av_full_orb_2_rdm.irp.f

@@ -16,7 +16,15 @@
 !
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
  END_DOC
+ PROVIDE n_core_orb list_core
  state_av_full_occ_2_rdm_ab_mo = 0.d0
+ !$OMP PARALLEL PRIVATE(i,j,k,l,iorb,jorb,korb,lorb) &
+ !$OMP DEFAULT(NONE) SHARED(n_act_orb, n_inact_orb, n_core_orb, &
+ !$OMP list_core, list_act, list_inact, no_core_density, &
+ !$OMP one_e_dm_mo_alpha_average, one_e_dm_mo_beta_average, &
+ !$OMP state_av_act_2_rdm_ab_mo, state_av_full_occ_2_rdm_ab_mo)
+
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -32,8 +40,10 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
    !! BETA ACTIVE - ALPHA inactive
    !!
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -45,9 +55,11 @@
      enddo
     enddo
    enddo
+ !$OMP END DO
 
    !! ALPHA ACTIVE - BETA inactive
    !!
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -59,9 +71,11 @@
      enddo
     enddo
    enddo
+ !$OMP END DO
 
    !! ALPHA INACTIVE - BETA INACTIVE
    !!
+ !$OMP DO
     do j = 1, n_inact_orb
      jorb = list_inact(j)
      do k = 1, n_inact_orb
@@ -70,6 +84,7 @@
       state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
      enddo
     enddo
+ !$OMP END DO
 
 !!!!!!!!!!!!
 !!!!!!!!!!!! if "no_core_density" then you don't put the core part
@@ -77,6 +92,7 @@
    if (.not.no_core_density)then
     !! BETA ACTIVE - ALPHA CORE
     !!
+ !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -88,9 +104,11 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
 
     !! ALPHA ACTIVE - BETA CORE
     !!
+ !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -102,9 +120,11 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
 
    !! ALPHA CORE - BETA CORE
    !!
+ !$OMP DO
     do j = 1, n_core_orb
      jorb = list_core(j)
      do k = 1, n_core_orb
@@ -113,8 +133,10 @@
       state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
      enddo
     enddo
+ !$OMP END DO
    endif
 
+ !$OMP END PARALLEL
  END_PROVIDER
 
 
@@ -136,8 +158,15 @@
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
  END_DOC
 
+ PROVIDE n_core_orb list_core
+ !$OMP PARALLEL PRIVATE(i,j,k,l,iorb,jorb,korb,lorb) &
+ !$OMP DEFAULT(NONE) SHARED(n_act_orb, n_inact_orb, n_core_orb, &
+ !$OMP list_core, list_act, list_inact, no_core_density, &
+ !$OMP one_e_dm_mo_alpha_average, one_e_dm_mo_beta_average, &
+ !$OMP state_av_act_2_rdm_aa_mo, state_av_full_occ_2_rdm_aa_mo)
    !! PURE ACTIVE PART ALPHA-ALPHA
    !!
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -152,8 +181,10 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
    !! ALPHA ACTIVE - ALPHA inactive
    !!
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -169,8 +200,10 @@
      enddo
     enddo
    enddo
+ !$OMP END DO
 
    !! ALPHA INACTIVE - ALPHA INACTIVE
+ !$OMP DO
    do j = 1, n_inact_orb
     jorb = list_inact(j)
     do k = 1, n_inact_orb
@@ -179,12 +212,14 @@
      state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -=  1.0d0
     enddo
    enddo
+ !$OMP END DO
 
 !!!!!!!!!!
 !!!!!!!!!! if "no_core_density" then you don't put the core part
 !!!!!!!!!! CAN BE USED
    if (.not.no_core_density)then
     !! ALPHA ACTIVE - ALPHA CORE
+ !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -200,8 +235,10 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
     !! ALPHA CORE - ALPHA CORE
 
+ !$OMP DO
     do j = 1, n_core_orb
      jorb = list_core(j)
      do k = 1, n_core_orb
@@ -210,8 +247,10 @@
       state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -=  1.0d0
      enddo
     enddo
+ !$OMP END DO
    endif
 
+ !$OMP END PARALLEL
  END_PROVIDER
 
  BEGIN_PROVIDER [double precision, state_av_full_occ_2_rdm_bb_mo, (n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb,n_core_inact_act_orb)]
@@ -232,8 +271,15 @@
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
  END_DOC
 
+ PROVIDE n_core_orb list_core
+ !$OMP PARALLEL PRIVATE(i,j,k,l,iorb,jorb,korb,lorb) &
+ !$OMP DEFAULT(NONE) SHARED(n_act_orb, n_inact_orb, n_core_orb, &
+ !$OMP list_core, list_act, list_inact, no_core_density, &
+ !$OMP one_e_dm_mo_alpha_average, one_e_dm_mo_beta_average, &
+ !$OMP state_av_act_2_rdm_bb_mo, state_av_full_occ_2_rdm_bb_mo)
    !! PURE ACTIVE PART beta-beta
    !!
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -248,8 +294,10 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
    !! beta ACTIVE - beta inactive
    !!
+ !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -265,8 +313,10 @@
      enddo
     enddo
    enddo
+ !$OMP END DO
 
    !! beta INACTIVE - beta INACTIVE
+ !$OMP DO
    do j = 1, n_inact_orb
     jorb = list_inact(j)
     do k = 1, n_inact_orb
@@ -275,12 +325,14 @@
      state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -=  1.0d0
     enddo
    enddo
+ !$OMP END DO
 
 !!!!!!!!!!!!
 !!!!!!!!!!!! if "no_core_density" then you don't put the core part
 !!!!!!!!!!!! CAN BE USED
    if (.not.no_core_density)then
     !! beta ACTIVE - beta CORE
+ !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -296,8 +348,10 @@
       enddo
      enddo
     enddo
+ !$OMP END DO
     !! beta CORE - beta CORE
 
+ !$OMP DO
     do j = 1, n_core_orb
      jorb = list_core(j)
      do k = 1, n_core_orb
@@ -306,7 +360,9 @@
       state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -=  1.0d0
      enddo
     enddo
+ !$OMP END DO
    endif
+ !$OMP END PARALLEL
 
  END_PROVIDER
 
@@ -329,9 +385,17 @@
 !  !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
  END_DOC
 
+ PROVIDE n_core_orb list_core
+
+ !$OMP PARALLEL PRIVATE(i,j,k,l,iorb,jorb,korb,lorb) &
+ !$OMP DEFAULT(NONE) SHARED(n_act_orb, n_inact_orb, n_core_orb, &
+ !$OMP list_core, list_act, list_inact, no_core_density, &
+ !$OMP one_e_dm_mo_alpha_average, one_e_dm_mo_beta_average, &
+ !$OMP state_av_act_2_rdm_spin_trace_mo, state_av_full_occ_2_rdm_spin_trace_mo)
    !!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!
    !! PURE ACTIVE PART SPIN-TRACE
+   !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -346,11 +410,13 @@
       enddo
      enddo
     enddo
+   !$OMP END DO
 
    !!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!
    !!!!! BETA-BETA !!!!!
    !! beta ACTIVE - beta inactive
+   !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -366,7 +432,9 @@
      enddo
     enddo
    enddo
+   !$OMP END DO
    !! beta INACTIVE - beta INACTIVE
+   !$OMP DO
    do j = 1, n_inact_orb
     jorb = list_inact(j)
     do k = 1, n_inact_orb
@@ -375,8 +443,10 @@
      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0
     enddo
    enddo
+   !$OMP END DO
    if (.not.no_core_density)then
     !! beta ACTIVE - beta CORE
+    !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -392,7 +462,9 @@
       enddo
      enddo
     enddo
+    !$OMP END DO
     !! beta CORE - beta CORE
+    !$OMP DO
     do j = 1, n_core_orb
      jorb = list_core(j)
      do k = 1, n_core_orb
@@ -401,12 +473,14 @@
       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0
      enddo
     enddo
+    !$OMP END DO
    endif
 
    !!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!
    !!!!! ALPHA-ALPHA !!!!!
    !! ALPHA ACTIVE - ALPHA inactive
+   !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -422,7 +496,9 @@
      enddo
     enddo
    enddo
+   !$OMP END DO
    !! ALPHA INACTIVE - ALPHA INACTIVE
+   !$OMP DO
    do j = 1, n_inact_orb
     jorb = list_inact(j)
     do k = 1, n_inact_orb
@@ -431,8 +507,10 @@
      state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0
     enddo
    enddo
+   !$OMP END DO
    if (.not.no_core_density)then
     !! ALPHA ACTIVE - ALPHA CORE
+    !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -448,7 +526,9 @@
       enddo
      enddo
     enddo
+    !$OMP END DO
     !! ALPHA CORE - ALPHA CORE
+    !$OMP DO
     do j = 1, n_core_orb
      jorb = list_core(j)
      do k = 1, n_core_orb
@@ -457,11 +537,13 @@
       state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -=  1.0d0
      enddo
     enddo
+    !$OMP END DO
    endif
 
    !!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!
    !!!!! ALPHA-BETA + BETA-ALPHA !!!!!
+   !$OMP DO
    do i = 1, n_act_orb
     iorb = list_act(i)
     do j = 1, n_act_orb
@@ -481,7 +563,9 @@
      enddo
     enddo
    enddo
+   !$OMP END DO
    !! ALPHA INACTIVE - BETA INACTIVE
+    !$OMP DO
     do j = 1, n_inact_orb
      jorb = list_inact(j)
      do k = 1, n_inact_orb
@@ -491,11 +575,13 @@
       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0d0
      enddo
     enddo
+    !$OMP END DO
 
 !!!!!!!!!!!!
 !!!!!!!!!!!! if "no_core_density" then you don't put the core part
 !!!!!!!!!!!! CAN BE USED
    if (.not.no_core_density)then
+    !$OMP DO
     do i = 1, n_act_orb
      iorb = list_act(i)
      do j = 1, n_act_orb
@@ -515,7 +601,9 @@
       enddo
      enddo
     enddo
+    !$OMP END DO
    !! ALPHA CORE - BETA CORE
+    !$OMP DO
     do j = 1, n_core_orb
      jorb = list_core(j)
      do k = 1, n_core_orb
@@ -525,7 +613,9 @@
       state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0D0
      enddo
     enddo
+    !$OMP END DO
 
    endif
+  !$OMP END PARALLEL
 
  END_PROVIDER

src/two_body_rdm/two_e_dm_mo.irp.f

@@ -16,6 +16,9 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)]
    two_e_dm_mo = 0.d0
    integer                        :: i,j,k,l,iorb,jorb,korb,lorb,istate
 
+   !$OMP PARALLEL DO PRIVATE(i,j,k,l,iorb,jorb,korb,lorb) &
+   !$OMP DEFAULT(NONE) SHARED(n_core_inact_act_orb, list_core_inact_act, &
+   !$OMP  two_e_dm_mo, state_av_full_occ_2_rdm_spin_trace_mo)
    do l=1,n_core_inact_act_orb
     lorb = list_core_inact_act(l)
     do k=1,n_core_inact_act_orb
@@ -29,7 +32,7 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)]
      enddo
     enddo
    enddo
-   two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) 
+   !$OMP END PARALLEL DO
 
  END_PROVIDER