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

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)
@@ -83,7 +85,7 @@ subroutine gradvec_tc_ia(i,a,res_l, res_r)
  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))
+  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)
@@ -93,12 +95,14 @@ subroutine gradvec_tc_ia(i,a,res_l, res_r)
    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 * (                                                                                    &
+   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(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(r,a,t,i) - mo_bi_ortho_tc_two_e(r,a,i,t))  &
+  +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)
@@ -116,20 +120,22 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
  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_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
@@ -140,14 +146,21 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
     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)
@@ -161,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)