diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 9b8dba78..7b243bb4 100644 --- a/src/casscf/densities.irp.f +++ b/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 diff --git a/src/casscf/test_two_rdm.irp.f b/src/casscf/test_two_rdm.irp.f index 562d15a6..9abe0aa0 100644 --- a/src/casscf/test_two_rdm.irp.f +++ b/src/casscf/test_two_rdm.irp.f @@ -8,11 +8,11 @@ program print_two_rdm 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) + 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)) if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then print*,'' print*,'sum' diff --git a/src/two_body_rdm/two_rdm.irp.f b/src/two_body_rdm/two_rdm.irp.f index 1c299bba..a75a92cc 100644 --- a/src/two_body_rdm/two_rdm.irp.f +++ b/src/two_body_rdm/two_rdm.irp.f @@ -1,23 +1,27 @@ - - BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] +BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num)] 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 + integer :: i,j,k,l, istate + coussin_peter_two_rdm_mo = 0.d0 + 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) = & + state_average_weight(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 +END_PROVIDER BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]