scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-06-02 11:25:23 +02:00
Code Issues Releases Wiki Activity

changed hessian in order to better select the blocks

Browse Source
This commit is contained in:
Emmanuel Giner 2021-07-02 16:02:33 +02:00
parent 1edbd0b890
commit 98c3948d6a
5 changed files with 528 additions and 313 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
devel/casscf/gradient.irp.f
View File

@ -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

486
devel/casscf/hessian.irp.f
View File

@ -1,312 +1,5 @@
use bitmasks
BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
BEGIN_DOC
! calculate the orbital hessian 2 <Psi| E_pq H E_rs |Psi>
! + <Psi| E_pq E_rs H |Psi> + <Psi| E_rs E_pq H |Psi> 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 <Psi| E_pq H E_rs |Psi>
! + <Psi| E_pq E_rs H |Psi> + <Psi| E_rs E_pq H |Psi>
! 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

310
devel/casscf/hessian_old.irp.f Normal file
View File

@ -0,0 +1,310 @@
use bitmasks
BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)]
BEGIN_DOC
! calculate the orbital hessian 2 <Psi| E_pq H E_rs |Psi>
! + <Psi| E_pq E_rs H |Psi> + <Psi| E_rs E_pq H |Psi> 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 <Psi| E_pq H E_rs |Psi>
! + <Psi| E_pq E_rs H |Psi> + <Psi| E_rs E_pq H |Psi>
! 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

6
devel/casscf/neworbs.irp.f
View File

@ -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

4
devel/casscf/swap_orb.irp.f
View File

@ -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