From 98c3948d6a0dfeae44a9285e3969c6c69cda4ffe Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 2 Jul 2021 16:02:33 +0200 Subject: [PATCH] changed hessian in order to better select the blocks --- devel/casscf/gradient.irp.f | 35 ++- devel/casscf/hessian.irp.f | 486 ++++++++++++--------------------- devel/casscf/hessian_old.irp.f | 310 +++++++++++++++++++++ devel/casscf/neworbs.irp.f | 6 +- devel/casscf/swap_orb.irp.f | 4 +- 5 files changed, 528 insertions(+), 313 deletions(-) create mode 100644 devel/casscf/hessian_old.irp.f diff --git a/devel/casscf/gradient.irp.f b/devel/casscf/gradient.irp.f index 6795aa4..309371b 100644 --- a/devel/casscf/gradient.irp.f +++ b/devel/casscf/gradient.irp.f @@ -8,15 +8,31 @@ BEGIN_PROVIDER [ integer, nMonoEx ] nMonoEx=n_core_inact_orb*n_act_orb+n_core_inact_orb*n_virt_orb+n_act_orb*n_virt_orb END_PROVIDER + BEGIN_PROVIDER [integer, n_c_a_prov] +&BEGIN_PROVIDER [integer, n_c_v_prov] +&BEGIN_PROVIDER [integer, n_a_v_prov] + implicit none + n_c_a_prov = n_core_inact_orb * n_act_orb + n_c_v_prov = n_core_inact_orb * n_virt_orb + n_a_v_prov = n_act_orb * n_virt_orb + END_PROVIDER + BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] &BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] +&BEGIN_PROVIDER [integer, list_idx_c_a, (3,n_c_a_prov) ] +&BEGIN_PROVIDER [integer, list_idx_c_v, (3,n_c_v_prov) ] +&BEGIN_PROVIDER [integer, list_idx_a_v, (3,n_a_v_prov) ] +&BEGIN_PROVIDER [integer, mat_idx_c_a, (n_core_inact_orb,n_act_orb) +&BEGIN_PROVIDER [integer, mat_idx_c_v, (n_core_inact_orb,n_virt_orb) +&BEGIN_PROVIDER [integer, mat_idx_a_v, (n_act_orb,n_virt_orb) BEGIN_DOC ! a list of the orbitals involved in the excitation END_DOC implicit none - integer :: i,t,a,ii,tt,aa,indx + integer :: i,t,a,ii,tt,aa,indx,indx_tmp indx=0 + indx_tmp = 0 do ii=1,n_core_inact_orb i=list_core_inact(ii) do tt=1,n_act_orb @@ -25,9 +41,15 @@ END_PROVIDER excit(1,indx)=i excit(2,indx)=t excit_class(indx)='c-a' + indx_tmp += 1 + list_idx_c_a(1,indx_tmp) = indx + list_idx_c_a(2,indx_tmp) = ii + list_idx_c_a(3,indx_tmp) = tt + mat_idx_c_a(ii,tt) = indx end do end do + indx_tmp = 0 do ii=1,n_core_inact_orb i=list_core_inact(ii) do aa=1,n_virt_orb @@ -36,9 +58,15 @@ END_PROVIDER excit(1,indx)=i excit(2,indx)=a excit_class(indx)='c-v' + indx_tmp += 1 + list_idx_c_v(1,indx_tmp) = indx + list_idx_c_v(2,indx_tmp) = ii + list_idx_c_v(3,indx_tmp) = aa + mat_idx_c_v(ii,aa) = indx end do end do + indx_tmp = 0 do tt=1,n_act_orb t=list_act(tt) do aa=1,n_virt_orb @@ -47,6 +75,11 @@ END_PROVIDER excit(1,indx)=t excit(2,indx)=a excit_class(indx)='a-v' + indx_tmp += 1 + list_idx_a_v(1,indx_tmp) = indx + list_idx_a_v(2,indx_tmp) = tt + list_idx_a_v(3,indx_tmp) = aa + mat_idx_a_v(tt,aa) = indx end do end do diff --git a/devel/casscf/hessian.irp.f b/devel/casscf/hessian.irp.f index c63b903..61864ef 100644 --- a/devel/casscf/hessian.irp.f +++ b/devel/casscf/hessian.irp.f @@ -1,312 +1,5 @@ use bitmasks -BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)] - BEGIN_DOC - ! calculate the orbital hessian 2 - ! + + by hand, - ! determinant per determinant, as for the gradient - ! - ! we assume that we have natural active orbitals - END_DOC - implicit none - integer :: indx,ihole,ipart - integer :: jndx,jhole,jpart - character*3 :: iexc,jexc - real*8 :: res - - if (bavard) then - write(6,*) ' providing Hessian matrix hessmat ' - write(6,*) ' nMonoEx = ',nMonoEx - endif - - do indx=1,nMonoEx - do jndx=1,nMonoEx - hessmat(indx,jndx)=0.D0 - end do - end do - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - iexc=excit_class(indx) - do jndx=indx,nMonoEx - jhole=excit(1,jndx) - jpart=excit(2,jndx) - jexc=excit_class(jndx) - call calc_hess_elem(ihole,ipart,jhole,jpart,res) - hessmat(indx,jndx)=res - hessmat(jndx,indx)=res - end do - end do - -END_PROVIDER - -subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) - BEGIN_DOC - ! eq 19 of Siegbahn et al, Physica Scripta 1980 - ! we calculate 2 - ! + + - ! average over all states is performed. - ! no transition between states. - END_DOC - implicit none - integer :: ihole,ipart,ispin,mu,istate - integer :: jhole,jpart,jspin - integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:) - integer(bit_kind), allocatable :: det_nu(:,:) - integer(bit_kind), allocatable :: det_mu_pq(:,:) - integer(bit_kind), allocatable :: det_mu_rs(:,:) - integer(bit_kind), allocatable :: det_nu_rs(:,:) - integer(bit_kind), allocatable :: det_mu_pqrs(:,:) - integer(bit_kind), allocatable :: det_mu_rspq(:,:) - real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 - real*8 :: i_H_j_element - allocate(det_mu(N_int,2)) - allocate(det_nu(N_int,2)) - allocate(det_mu_pq(N_int,2)) - allocate(det_mu_rs(N_int,2)) - allocate(det_nu_rs(N_int,2)) - allocate(det_mu_pqrs(N_int,2)) - allocate(det_mu_rspq(N_int,2)) - integer :: mu_pq_possible - integer :: mu_rs_possible - integer :: nu_rs_possible - integer :: mu_pqrs_possible - integer :: mu_rspq_possible - - res=0.D0 - - ! the terms <0|E E H |0> - do mu=1,n_det - ! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 - ! do the monoexcitation pq on it - call det_copy(det_mu,det_mu_pq,N_int) - call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & - ,ihole,ipart,ispin,phase,mu_pq_possible) - if (mu_pq_possible.eq.1) then - ! possible, but not necessarily in the list - ! do the second excitation - do jspin=1,2 - call det_copy(det_mu_pq,det_mu_pqrs,N_int) - call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& - ,jhole,jpart,jspin,phase2,mu_pqrs_possible) - ! excitation possible - if (mu_pqrs_possible.eq.1) then - call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 - end do - end if - ! try the de-excitation with opposite sign - call det_copy(det_mu_pq,det_mu_pqrs,N_int) - call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& - ,jpart,jhole,jspin,phase2,mu_pqrs_possible) - phase2=-phase2 - ! excitation possible - if (mu_pqrs_possible.eq.1) then - call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 - end do - end if - end do - end if - ! exchange the notion of pq and rs - ! do the monoexcitation rs on the initial determinant - call det_copy(det_mu,det_mu_rs,N_int) - call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & - ,jhole,jpart,ispin,phase2,mu_rs_possible) - if (mu_rs_possible.eq.1) then - ! do the second excitation - do jspin=1,2 - call det_copy(det_mu_rs,det_mu_rspq,N_int) - call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& - ,ihole,ipart,jspin,phase3,mu_rspq_possible) - ! excitation possible (of course, the result is outside the CAS) - if (mu_rspq_possible.eq.1) then - call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 - end do - end if - ! we may try the de-excitation, with opposite sign - call det_copy(det_mu_rs,det_mu_rspq,N_int) - call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& - ,ipart,ihole,jspin,phase3,mu_rspq_possible) - phase3=-phase3 - ! excitation possible (of course, the result is outside the CAS) - if (mu_rspq_possible.eq.1) then - call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 - end do - end if - end do - end if - ! - ! the operator E H E, we have to do a double loop over the determinants - ! we still have the determinant mu_pq and the phase in memory - if (mu_pq_possible.eq.1) then - do nu=1,N_det - call det_extract(det_nu,nu,N_int) - do jspin=1,2 - call det_copy(det_nu,det_nu_rs,N_int) - call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & - ,jhole,jpart,jspin,phase2,nu_rs_possible) - ! excitation possible ? - if (nu_rs_possible.eq.1) then - call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) - do istate=1,N_states - res+=2.D0*i_H_j_element*psi_coef(mu,istate) & - *psi_coef(nu,istate)*phase*phase2 - end do - end if - end do - end do - end if - end do - end do - - ! state-averaged Hessian - res*=1.D0/dble(N_states) - -end subroutine calc_hess_elem - -BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] - BEGIN_DOC - ! explicit hessian matrix from density matrices and integrals - ! of course, this will be used for a direct Davidson procedure later - ! we will not store the matrix in real life - ! formulas are broken down as functions for the 6 classes of matrix elements - ! - END_DOC - implicit none - integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart,indx_shift - - real*8 :: hessmat_itju - real*8 :: hessmat_itja - real*8 :: hessmat_itua - real*8 :: hessmat_iajb - real*8 :: hessmat_iatb - real*8 :: hessmat_taub - - if (bavard) then - write(6,*) ' providing Hessian matrix hessmat2 ' - write(6,*) ' nMonoEx = ',nMonoEx - endif - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(hessmat2,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & - !$OMP PRIVATE(i,indx,jndx,j,ustart,t,u,a,bstart,indx_shift) - - !$OMP DO - ! (DOUBLY OCCUPIED ---> ACT ) - do i=1,n_core_inact_orb - do t=1,n_act_orb - indx = t + (i-1)*n_act_orb - jndx=indx - ! (DOUBLY OCCUPIED ---> ACT ) - do j=i,n_core_inact_orb - if (i.eq.j) then - ustart=t - else - ustart=1 - end if - do u=ustart,n_act_orb - hessmat2(jndx,indx)=hessmat_itju(i,t,j,u) - jndx+=1 - end do - end do - ! (DOUBLY OCCUPIED ---> VIRTUAL) - do j=1,n_core_inact_orb - do a=1,n_virt_orb - hessmat2(jndx,indx)=hessmat_itja(i,t,j,a) - jndx+=1 - end do - end do - ! (ACTIVE ---> VIRTUAL) - do u=1,n_act_orb - do a=1,n_virt_orb - hessmat2(jndx,indx)=hessmat_itua(i,t,u,a) - jndx+=1 - end do - end do - end do - end do - !$OMP END DO NOWAIT - - indx_shift = n_core_inact_orb*n_act_orb - !$OMP DO - ! (DOUBLY OCCUPIED ---> VIRTUAL) - do a=1,n_virt_orb - do i=1,n_core_inact_orb - indx = a + (i-1)*n_virt_orb + indx_shift - jndx=indx - ! (DOUBLY OCCUPIED ---> VIRTUAL) - do j=i,n_core_inact_orb - if (i.eq.j) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - hessmat2(jndx,indx)=hessmat_iajb(i,a,j,b) - jndx+=1 - end do - end do - ! (ACT ---> VIRTUAL) - do t=1,n_act_orb - do b=1,n_virt_orb - hessmat2(jndx,indx)=hessmat_iatb(i,a,t,b) - jndx+=1 - end do - end do - end do - end do - !$OMP END DO NOWAIT - - indx_shift += n_core_inact_orb*n_virt_orb - !$OMP DO - ! (ACT ---> VIRTUAL) - do a=1,n_virt_orb - do t=1,n_act_orb - indx = a + (t-1)*n_virt_orb + indx_shift - jndx=indx - ! (ACT ---> VIRTUAL) - do u=t,n_act_orb - if (t.eq.u) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - hessmat2(jndx,indx)=hessmat_taub(t,a,u,b) - jndx+=1 - end do - end do - end do - end do - !$OMP END DO - - !$OMP END PARALLEL - - do jndx=1,nMonoEx - do indx=1,jndx-1 - hessmat2(indx,jndx) = hessmat2(jndx,indx) - enddo - enddo - - -END_PROVIDER - real*8 function hessmat_itju(i,t,j,u) BEGIN_DOC ! the orbital hessian for core/inactive -> active, core/inactive -> active @@ -663,3 +356,182 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] !$OMP END PARALLEL END_PROVIDER + + +BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] + implicit none + integer :: i,j,t,u,a,b + integer :: indx,indx_tmp, jndx, jndx_tmp + integer :: ustart,bstart + real*8 :: hessmat_itju + real*8 :: hessmat_itja + real*8 :: hessmat_itua + real*8 :: hessmat_iajb + real*8 :: hessmat_iatb + real*8 :: hessmat_taub + ! c-a c-v a-v + ! c-a | X X X + ! c-v | X X + ! a-v | X + + provide mo_two_e_integrals_in_map + + hessmat = 0.d0 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_a_prov,list_idx_c_a,n_core_inact_orb,n_act_orb,mat_idx_c_a) & + !$OMP PRIVATE(indx_tmp,indx,i,t,j,u,ustart,jndx) + + !$OMP DO +!!!! < Core-active| H |Core-active > + ! Core-active excitations + do indx_tmp = 1, n_c_a_prov + indx = list_idx_c_a(1,indx_tmp) + i = list_idx_c_a(2,indx_tmp) + t = list_idx_c_a(3,indx_tmp) + ! Core-active excitations + do j = 1, n_core_inact_orb + if (i.eq.j) then + ustart=t + else + ustart=1 + end if + do u=ustart,n_act_orb + jndx = mat_idx_c_a(j,u) + hessmat(jndx,indx) = hessmat_itju(i,t,j,u) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_a_prov,n_c_v_prov,list_idx_c_a,list_idx_c_v) & + !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,t,j,a,jndx) + + !$OMP DO +!!!! < Core-active| H |Core-VIRTUAL > + ! Core-active excitations + do indx_tmp = 1, n_c_a_prov + indx = list_idx_c_a(1,indx_tmp) + i = list_idx_c_a(2,indx_tmp) + t = list_idx_c_a(3,indx_tmp) + ! Core-VIRTUAL excitations + do jndx_tmp = 1, n_c_v_prov + jndx = list_idx_c_v(1,jndx_tmp) + j = list_idx_c_v(2,jndx_tmp) + a = list_idx_c_v(3,jndx_tmp) + hessmat(jndx,indx) = hessmat_itja(i,t,j,a) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_a_prov,n_a_v_prov,list_idx_c_a,list_idx_a_v) & + !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,t,u,a,jndx) + + !$OMP DO +!!!! < Core-active| H |ACTIVE-VIRTUAL > + ! Core-active excitations + do indx_tmp = 1, n_c_a_prov + indx = list_idx_c_a(1,indx_tmp) + i = list_idx_c_a(2,indx_tmp) + t = list_idx_c_a(3,indx_tmp) + ! ACTIVE-VIRTUAL excitations + do jndx_tmp = 1, n_a_v_prov + jndx = list_idx_a_v(1,jndx_tmp) + u = list_idx_a_v(2,jndx_tmp) + a = list_idx_a_v(3,jndx_tmp) + hessmat(jndx,indx) = hessmat_itua(i,t,u,a) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_v_prov,list_idx_c_v,n_core_inact_orb,n_virt_orb,mat_idx_c_v) & + !$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx) + + !$OMP DO +!!!! < Core-VIRTUAL | H |Core-VIRTUAL > + ! Core-VIRTUAL excitations + do indx_tmp = 1, n_c_v_prov + indx = list_idx_c_v(1,indx_tmp) + i = list_idx_c_v(2,indx_tmp) + a = list_idx_c_v(3,indx_tmp) + ! Core-VIRTUAL excitations + do j = 1, n_core_inact_orb + if (i.eq.j) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + jndx = mat_idx_c_v(j,b) + hessmat(jndx,indx) = hessmat_iajb(i,a,j,b) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + enddo + + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_c_v_prov,n_a_v_prov,list_idx_c_v,list_idx_a_v) & + !$OMP PRIVATE(indx_tmp,jndx_tmp,indx,i,a,t,b,jndx) + + !$OMP DO +!!!! < Core-VIRTUAL | H |Active-VIRTUAL > + ! Core-VIRTUAL excitations + do indx_tmp = 1, n_c_v_prov + indx = list_idx_c_v(1,indx_tmp) + i = list_idx_c_v(2,indx_tmp) + a = list_idx_c_v(3,indx_tmp) + ! Active-VIRTUAL excitations + do jndx_tmp = 1, n_a_v_prov + jndx = list_idx_a_v(1,jndx_tmp) + t = list_idx_a_v(2,jndx_tmp) + b = list_idx_a_v(3,jndx_tmp) + hessmat(jndx,indx) = hessmat_iatb(i,a,t,b) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat,n_a_v_prov,list_idx_a_v,n_act_orb,n_virt_orb,mat_idx_a_v) & + !$OMP PRIVATE(indx_tmp,indx,t,a,u,b,bstart,jndx) + + !$OMP DO +!!!! < Active-VIRTUAL | H |Active-VIRTUAL > + ! Active-VIRTUAL excitations + do indx_tmp = 1, n_a_v_prov + indx = list_idx_a_v(1,indx_tmp) + t = list_idx_a_v(2,indx_tmp) + a = list_idx_a_v(3,indx_tmp) + ! Active-VIRTUAL excitations + do u=t,n_act_orb + if (t.eq.u) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + jndx = mat_idx_a_v(u,b) + hessmat(jndx,indx) = hessmat_taub(t,a,u,b) + hessmat(indx,jndx) = hessmat(jndx,indx) + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + +END_PROVIDER diff --git a/devel/casscf/hessian_old.irp.f b/devel/casscf/hessian_old.irp.f new file mode 100644 index 0000000..d17f1f0 --- /dev/null +++ b/devel/casscf/hessian_old.irp.f @@ -0,0 +1,310 @@ + +use bitmasks +BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)] + BEGIN_DOC + ! calculate the orbital hessian 2 + ! + + by hand, + ! determinant per determinant, as for the gradient + ! + ! we assume that we have natural active orbitals + END_DOC + implicit none + integer :: indx,ihole,ipart + integer :: jndx,jhole,jpart + character*3 :: iexc,jexc + real*8 :: res + + if (bavard) then + write(6,*) ' providing Hessian matrix hessmat_old ' + write(6,*) ' nMonoEx = ',nMonoEx + endif + + do indx=1,nMonoEx + do jndx=1,nMonoEx + hessmat_old(indx,jndx)=0.D0 + end do + end do + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + iexc=excit_class(indx) + do jndx=indx,nMonoEx + jhole=excit(1,jndx) + jpart=excit(2,jndx) + jexc=excit_class(jndx) + call calc_hess_elem(ihole,ipart,jhole,jpart,res) + hessmat_old(indx,jndx)=res + hessmat_old(jndx,indx)=res + end do + end do + +END_PROVIDER + +subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) + BEGIN_DOC + ! eq 19 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 + ! + + + ! average over all states is performed. + ! no transition between states. + END_DOC + implicit none + integer :: ihole,ipart,ispin,mu,istate + integer :: jhole,jpart,jspin + integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:) + integer(bit_kind), allocatable :: det_nu(:,:) + integer(bit_kind), allocatable :: det_mu_pq(:,:) + integer(bit_kind), allocatable :: det_mu_rs(:,:) + integer(bit_kind), allocatable :: det_nu_rs(:,:) + integer(bit_kind), allocatable :: det_mu_pqrs(:,:) + integer(bit_kind), allocatable :: det_mu_rspq(:,:) + real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 + real*8 :: i_H_j_element + allocate(det_mu(N_int,2)) + allocate(det_nu(N_int,2)) + allocate(det_mu_pq(N_int,2)) + allocate(det_mu_rs(N_int,2)) + allocate(det_nu_rs(N_int,2)) + allocate(det_mu_pqrs(N_int,2)) + allocate(det_mu_rspq(N_int,2)) + integer :: mu_pq_possible + integer :: mu_rs_possible + integer :: nu_rs_possible + integer :: mu_pqrs_possible + integer :: mu_rspq_possible + + res=0.D0 + + ! the terms <0|E E H |0> + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation pq on it + call det_copy(det_mu,det_mu_pq,N_int) + call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & + ,ihole,ipart,ispin,phase,mu_pq_possible) + if (mu_pq_possible.eq.1) then + ! possible, but not necessarily in the list + ! do the second excitation + do jspin=1,2 + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& + ,jhole,jpart,jspin,phase2,mu_pqrs_possible) + ! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 + end do + end if + ! try the de-excitation with opposite sign + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& + ,jpart,jhole,jspin,phase2,mu_pqrs_possible) + phase2=-phase2 + ! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 + end do + end if + end do + end if + ! exchange the notion of pq and rs + ! do the monoexcitation rs on the initial determinant + call det_copy(det_mu,det_mu_rs,N_int) + call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & + ,jhole,jpart,ispin,phase2,mu_rs_possible) + if (mu_rs_possible.eq.1) then + ! do the second excitation + do jspin=1,2 + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& + ,ihole,ipart,jspin,phase3,mu_rspq_possible) + ! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 + end do + end if + ! we may try the de-excitation, with opposite sign + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& + ,ipart,ihole,jspin,phase3,mu_rspq_possible) + phase3=-phase3 + ! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 + end do + end if + end do + end if + ! + ! the operator E H E, we have to do a double loop over the determinants + ! we still have the determinant mu_pq and the phase in memory + if (mu_pq_possible.eq.1) then + do nu=1,N_det + call det_extract(det_nu,nu,N_int) + do jspin=1,2 + call det_copy(det_nu,det_nu_rs,N_int) + call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & + ,jhole,jpart,jspin,phase2,nu_rs_possible) + ! excitation possible ? + if (nu_rs_possible.eq.1) then + call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) + do istate=1,N_states + res+=2.D0*i_H_j_element*psi_coef(mu,istate) & + *psi_coef(nu,istate)*phase*phase2 + end do + end if + end do + end do + end if + end do + end do + + ! state-averaged Hessian + res*=1.D0/dble(N_states) + +end subroutine calc_hess_elem + +BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] + BEGIN_DOC + ! explicit hessian matrix from density matrices and integrals + ! of course, this will be used for a direct Davidson procedure later + ! we will not store the matrix in real life + ! formulas are broken down as functions for the 6 classes of matrix elements + ! + END_DOC + implicit none + integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart,indx_shift + + real*8 :: hessmat_itju + real*8 :: hessmat_itja + real*8 :: hessmat_itua + real*8 :: hessmat_iajb + real*8 :: hessmat_iatb + real*8 :: hessmat_taub + + if (bavard) then + write(6,*) ' providing Hessian matrix hessmat_peter ' + write(6,*) ' nMonoEx = ',nMonoEx + endif + provide mo_two_e_integrals_in_map + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(hessmat_peter,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & + !$OMP PRIVATE(i,indx,jndx,j,ustart,t,u,a,bstart,indx_shift) + + !$OMP DO + ! (DOUBLY OCCUPIED ---> ACT ) + do i=1,n_core_inact_orb + do t=1,n_act_orb + indx = t + (i-1)*n_act_orb + jndx=indx + ! (DOUBLY OCCUPIED ---> ACT ) + do j=i,n_core_inact_orb + if (i.eq.j) then + ustart=t + else + ustart=1 + end if + do u=ustart,n_act_orb + hessmat_peter(jndx,indx)=hessmat_itju(i,t,j,u) + jndx+=1 + end do + end do + ! (DOUBLY OCCUPIED ---> VIRTUAL) + do j=1,n_core_inact_orb + do a=1,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_itja(i,t,j,a) + jndx+=1 + end do + end do + ! (ACTIVE ---> VIRTUAL) + do u=1,n_act_orb + do a=1,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_itua(i,t,u,a) + jndx+=1 + end do + end do + end do + end do + !$OMP END DO NOWAIT + + indx_shift = n_core_inact_orb*n_act_orb + !$OMP DO + ! (DOUBLY OCCUPIED ---> VIRTUAL) + do a=1,n_virt_orb + do i=1,n_core_inact_orb + indx = a + (i-1)*n_virt_orb + indx_shift + jndx=indx + ! (DOUBLY OCCUPIED ---> VIRTUAL) + do j=i,n_core_inact_orb + if (i.eq.j) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_iajb(i,a,j,b) + jndx+=1 + end do + end do + ! (ACT ---> VIRTUAL) + do t=1,n_act_orb + do b=1,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_iatb(i,a,t,b) + jndx+=1 + end do + end do + end do + end do + !$OMP END DO NOWAIT + + indx_shift += n_core_inact_orb*n_virt_orb + !$OMP DO + ! (ACT ---> VIRTUAL) + do a=1,n_virt_orb + do t=1,n_act_orb + indx = a + (t-1)*n_virt_orb + indx_shift + jndx=indx + ! (ACT ---> VIRTUAL) + do u=t,n_act_orb + if (t.eq.u) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat_peter(jndx,indx)=hessmat_taub(t,a,u,b) + jndx+=1 + end do + end do + end do + end do + !$OMP END DO + + !$OMP END PARALLEL + + do jndx=1,nMonoEx + do indx=1,jndx-1 + hessmat_peter(indx,jndx) = hessmat_peter(jndx,indx) + enddo + enddo + + +END_PROVIDER + diff --git a/devel/casscf/neworbs.irp.f b/devel/casscf/neworbs.irp.f index 745cd59..3970395 100644 --- a/devel/casscf/neworbs.irp.f +++ b/devel/casscf/neworbs.irp.f @@ -23,8 +23,8 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] else do i=1,nMonoEx do j=1,nMonoEx - SXmatrix(i+1,j+1)=hessmat2(i,j) - SXmatrix(j+1,i+1)=hessmat2(i,j) + SXmatrix(i+1,j+1)=hessmat(i,j) + SXmatrix(j+1,i+1)=hessmat(i,j) end do end do endif @@ -34,7 +34,7 @@ BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)] enddo if (bavard) then do i=2,nMonoEx - write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i) + write(6,*) ' diagonal of the Hessian : ',i,hessmat(i,i) end do end if diff --git a/devel/casscf/swap_orb.irp.f b/devel/casscf/swap_orb.irp.f index 5d44215..49af207 100644 --- a/devel/casscf/swap_orb.irp.f +++ b/devel/casscf/swap_orb.irp.f @@ -45,12 +45,12 @@ imono = max_overlap(i) iorb = excit(1,imono) jorb = excit(2,imono) - if (excit_class(imono) == "c-a" .and.hessmat2(imono,imono).gt.0.d0)then ! core --> active rotation + if (excit_class(imono) == "c-a" .and.hessmat(imono,imono).gt.0.d0)then ! core --> active rotation n_orb_swap += 1 orb_swap(1,n_orb_swap) = iorb ! core orb_swap(2,n_orb_swap) = jorb ! active index_orb_swap(n_orb_swap) = imono - else if (excit_class(imono) == "a-v" .and.hessmat2(imono,imono).gt.0.d0)then ! active --> virtual rotation + else if (excit_class(imono) == "a-v" .and.hessmat(imono,imono).gt.0.d0)then ! active --> virtual rotation n_orb_swap += 1 orb_swap(1,n_orb_swap) = jorb ! virtual orb_swap(2,n_orb_swap) = iorb ! active