mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-01 05:33:40 +02:00
Compare commits
No commits in common. "2ef517488c9038b641a4f3c95ca01cb2d38b7181" and "6531181316c131f30d54a3653d17b597c0a43f3b" have entirely different histories.
2ef517488c
...
6531181316
@ -1,6 +1,6 @@
|
||||
! -*- F90 -*-
|
||||
BEGIN_PROVIDER [logical, bavard]
|
||||
bavard=.true.
|
||||
! bavard=.false.
|
||||
bavard=.false.
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -55,6 +55,7 @@
|
||||
end do
|
||||
end do
|
||||
|
||||
write(6,*) ' provided integrals (PQ|xx) '
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -115,6 +116,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_a
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' provided integrals (Px|xQ) '
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -144,5 +146,6 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' provided integrals (tu|xP) '
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -84,6 +84,7 @@
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' transformed PQxx'
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -175,6 +176,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' transformed PxxQ '
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -265,6 +267,7 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' transformed tuvP '
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -12,32 +12,24 @@ subroutine run
|
||||
implicit none
|
||||
double precision :: energy_old, energy
|
||||
logical :: converged
|
||||
integer :: iteration
|
||||
converged = .False.
|
||||
|
||||
energy = 0.d0
|
||||
mo_label = "MCSCF"
|
||||
iteration = 1
|
||||
do while (.not.converged)
|
||||
! do while (.not.converged)
|
||||
N_det = 1
|
||||
TOUCH N_det psi_det psi_coef
|
||||
call run_cipsi
|
||||
|
||||
write(6,*) ' total energy = ',eone+etwo+ecore
|
||||
mo_label = "MCSCF"
|
||||
mo_label = "Natural"
|
||||
mo_coef(:,:) = NatOrbsFCI(:,:)
|
||||
call save_mos
|
||||
|
||||
call driver_optorb
|
||||
energy_old = energy
|
||||
energy = eone+etwo+ecore
|
||||
|
||||
call write_time(6)
|
||||
call write_int(6,iteration,'CAS-SCF iteration')
|
||||
call write_double(6,energy,'CAS-SCF energy')
|
||||
call write_double(6,energy_improvement, 'Predicted energy improvement')
|
||||
|
||||
converged = dabs(energy_improvement) < thresh_scf
|
||||
|
||||
mo_coef = NewOrbs
|
||||
call save_mos
|
||||
call map_deinit(mo_integrals_map)
|
||||
N_det = 1
|
||||
iteration += 1
|
||||
FREE mo_integrals_map mo_two_e_integrals_in_map psi_det psi_coef
|
||||
SOFT_TOUCH mo_coef N_det
|
||||
enddo
|
||||
converged = dabs(energy - energy_old) < 1.d-10
|
||||
! enddo
|
||||
|
||||
end
|
||||
|
@ -1,19 +1,70 @@
|
||||
use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! the first-order density matrix in the basis of the starting MOs.
|
||||
! matrix is state averaged.
|
||||
! the first-order density matrix in the basis of the starting MOs
|
||||
! matrices are state averaged
|
||||
!
|
||||
! we use the spin-free generators of mono-excitations
|
||||
! E_pq destroys q and creates p
|
||||
! D_pq = <0|E_pq|0> = D_qp
|
||||
!
|
||||
END_DOC
|
||||
integer :: t,u
|
||||
implicit none
|
||||
integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
|
||||
integer :: ierr
|
||||
integer(bit_kind) :: det_mu(N_int,2)
|
||||
integer(bit_kind) :: det_mu_ex(N_int,2)
|
||||
integer(bit_kind) :: det_mu_ex1(N_int,2)
|
||||
integer(bit_kind) :: det_mu_ex2(N_int,2)
|
||||
real*8 :: phase1,phase2,term
|
||||
integer :: nu1,nu2
|
||||
integer :: ierr1,ierr2
|
||||
real*8 :: cI_mu(N_states)
|
||||
|
||||
do u=1,n_act_orb
|
||||
write(6,*) ' providing density matrices D0 and P0 '
|
||||
|
||||
D0tu = 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
|
||||
D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + &
|
||||
one_e_dm_mo_beta_average ( list_act(t), list_act(u) )
|
||||
enddo
|
||||
enddo
|
||||
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
|
||||
D0tu(t,u)+=term
|
||||
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
|
||||
D0tu(t,u)+=term
|
||||
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
|
||||
D0tu(v,x)*=1.0D0/dble(N_states)
|
||||
end do
|
||||
end do
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -39,9 +90,7 @@ BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
|
||||
integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12
|
||||
integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' providing density matrix P0'
|
||||
endif
|
||||
write(6,*) ' providing density matrices D0 and P0 '
|
||||
|
||||
P0tuvx = 0.d0
|
||||
|
||||
|
@ -31,8 +31,6 @@ subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, &
|
||||
! get the number in the list
|
||||
found=.false.
|
||||
nu=0
|
||||
|
||||
!TODO BOTTLENECK
|
||||
do while (.not.found)
|
||||
nu+=1
|
||||
if (nu.gt.N_det) then
|
||||
@ -52,6 +50,13 @@ subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, &
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
! if (found) then
|
||||
! if (nu.eq.-1) then
|
||||
! write(6,*) ' image not found in the list, thus nu = ',nu
|
||||
! else
|
||||
! write(6,*) ' found in the list as No ',nu,' phase = ',phase
|
||||
! end if
|
||||
! end if
|
||||
end if
|
||||
!
|
||||
! we found the new string, the phase, and possibly the number in the list
|
||||
|
@ -1,3 +1,32 @@
|
||||
subroutine driver_optorb
|
||||
implicit none
|
||||
end
|
||||
subroutine driver_optorb
|
||||
implicit none
|
||||
integer :: i,j
|
||||
|
||||
write(6,*)
|
||||
! write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1)
|
||||
write(6,*) ' energy improvement = ',energy_improvement
|
||||
! write(6,*) ' new energy = ',psi_energy_with_nucl_rep(1)+energy_improvement
|
||||
write(6,*)
|
||||
|
||||
write(6,*)
|
||||
write(6,*) ' creating new orbitals '
|
||||
do i=1,mo_num
|
||||
write(6,*) ' Orbital No ',i
|
||||
write(6,'(5F14.6)') (NewOrbs(j,i),j=1,mo_num)
|
||||
write(6,*)
|
||||
end do
|
||||
|
||||
mo_label = "Natural"
|
||||
do i=1,mo_num
|
||||
do j=1,ao_num
|
||||
mo_coef(j,i)=NewOrbs(j,i)
|
||||
end do
|
||||
end do
|
||||
call save_mos
|
||||
call map_deinit(mo_integrals_map)
|
||||
FREE mo_integrals_map mo_coef mo_two_e_integrals_in_map
|
||||
|
||||
write(6,*)
|
||||
write(6,*) ' ... all done '
|
||||
|
||||
end
|
||||
|
@ -6,6 +6,7 @@ BEGIN_PROVIDER [ integer, nMonoEx ]
|
||||
END_DOC
|
||||
implicit none
|
||||
nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
|
||||
write(6,*) ' nMonoEx = ',nMonoEx
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
|
||||
@ -86,11 +87,9 @@ BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
|
||||
norm_grad+=gradvec(indx)*gradvec(indx)
|
||||
end do
|
||||
norm_grad=sqrt(norm_grad)
|
||||
if (bavard) then
|
||||
write(6,*)
|
||||
write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
|
||||
write(6,*)
|
||||
endif
|
||||
write(6,*)
|
||||
write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
|
||||
write(6,*)
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
@ -119,11 +118,17 @@ subroutine calc_grad_elem(ihole,ipart,res)
|
||||
call do_signed_mono_excitation(det_mu,det_mu_ex,nu &
|
||||
,ihole,ipart,ispin,phase,ierr)
|
||||
if (ierr.eq.1) then
|
||||
! write(6,*)
|
||||
! write(6,*) ' mu = ',mu
|
||||
! call print_det(det_mu,N_int)
|
||||
! write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin
|
||||
! call print_det(det_mu_ex,N_int)
|
||||
call i_H_psi(det_mu_ex,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
|
||||
end do
|
||||
! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
@ -171,11 +176,9 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
|
||||
norm_grad+=gradvec2(indx)*gradvec2(indx)
|
||||
end do
|
||||
norm_grad=sqrt(norm_grad)
|
||||
if (bavard) then
|
||||
write(6,*)
|
||||
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
|
||||
write(6,*)
|
||||
endif
|
||||
write(6,*)
|
||||
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
|
||||
write(6,*)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -14,10 +14,8 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
|
||||
character*3 :: iexc,jexc
|
||||
real*8 :: res
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' providing Hessian matrix hessmat '
|
||||
write(6,*) ' nMonoEx = ',nMonoEx
|
||||
endif
|
||||
write(6,*) ' providing Hessian matrix hessmat '
|
||||
write(6,*) ' nMonoEx = ',nMonoEx
|
||||
|
||||
do indx=1,nMonoEx
|
||||
do jndx=1,nMonoEx
|
||||
@ -34,6 +32,8 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
|
||||
jpart=excit(2,jndx)
|
||||
jexc=excit_class(jndx)
|
||||
call calc_hess_elem(ihole,ipart,jhole,jpart,res)
|
||||
! write(6,*) ' Hessian ',ihole,'->',ipart &
|
||||
! ,' (',iexc,')',jhole,'->',jpart,' (',jexc,')',res
|
||||
hessmat(indx,jndx)=res
|
||||
hessmat(jndx,indx)=res
|
||||
end do
|
||||
@ -198,10 +198,8 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
|
||||
real*8 :: hessmat_iatb
|
||||
real*8 :: hessmat_taub
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' providing Hessian matrix hessmat2 '
|
||||
write(6,*) ' nMonoEx = ',nMonoEx
|
||||
endif
|
||||
write(6,*) ' providing Hessian matrix hessmat2 '
|
||||
write(6,*) ' nMonoEx = ',nMonoEx
|
||||
|
||||
indx=1
|
||||
do i=1,n_core_orb
|
||||
@ -216,6 +214,7 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
|
||||
do u=ustart,n_act_orb
|
||||
hessmat2(indx,jndx)=hessmat_itju(i,t,j,u)
|
||||
hessmat2(jndx,indx)=hessmat2(indx,jndx)
|
||||
! write(6,*) ' result I :',i,t,j,u,indx,jndx,hessmat(indx,jndx),hessmat2(indx,jndx)
|
||||
jndx+=1
|
||||
end do
|
||||
end do
|
||||
@ -295,6 +294,7 @@ real*8 function hessmat_itju(i,t,j,u)
|
||||
integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
|
||||
real*8 :: term,t2
|
||||
|
||||
! write(6,*) ' hessmat_itju ',i,t,j,u
|
||||
ii=list_core(i)
|
||||
tt=list_act(t)
|
||||
if (i.eq.j) then
|
||||
@ -340,6 +340,8 @@ real*8 function hessmat_itju(i,t,j,u)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!!! write(6,*) ' direct diff ',i,t,j,u,term,term2
|
||||
!!! term=term2
|
||||
end if
|
||||
else
|
||||
! it/ju
|
||||
@ -380,6 +382,7 @@ real*8 function hessmat_itja(i,t,j,a)
|
||||
integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
|
||||
real*8 :: term
|
||||
|
||||
! write(6,*) ' hessmat_itja ',i,t,j,a
|
||||
! it/ja
|
||||
ii=list_core(i)
|
||||
tt=list_act(t)
|
||||
@ -413,6 +416,7 @@ real*8 function hessmat_itua(i,t,u,a)
|
||||
integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3
|
||||
real*8 :: term
|
||||
|
||||
! write(6,*) ' hessmat_itua ',i,t,u,a
|
||||
ii=list_core(i)
|
||||
tt=list_act(t)
|
||||
t3=t+n_core_orb
|
||||
@ -453,6 +457,7 @@ real*8 function hessmat_iajb(i,a,j,b)
|
||||
implicit none
|
||||
integer :: i,a,j,b,ii,aa,jj,bb
|
||||
real*8 :: term
|
||||
! write(6,*) ' hessmat_iajb ',i,a,j,b
|
||||
|
||||
ii=list_core(i)
|
||||
aa=list_virt(a)
|
||||
@ -490,6 +495,7 @@ real*8 function hessmat_iatb(i,a,t,b)
|
||||
integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
|
||||
real*8 :: term
|
||||
|
||||
! write(6,*) ' hessmat_iatb ',i,a,t,b
|
||||
ii=list_core(i)
|
||||
aa=list_virt(a)
|
||||
tt=list_act(t)
|
||||
@ -546,6 +552,7 @@ real*8 function hessmat_taub(t,a,u,b)
|
||||
end do
|
||||
end do
|
||||
term=t1+t2+t3
|
||||
! write(6,*) ' Hess taub ',t,a,t1,t2,t3
|
||||
else
|
||||
bb=list_virt(b)
|
||||
! ta/tb b/=a
|
||||
|
@ -14,12 +14,10 @@
|
||||
occnum(list_act(i))=occ_act(n_act_orb-i+1)
|
||||
end do
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' occupation numbers '
|
||||
do i=1,mo_num
|
||||
write(6,*) i,occnum(i)
|
||||
end do
|
||||
endif
|
||||
write(6,*) ' occupation numbers '
|
||||
do i=1,mo_num
|
||||
write(6,*) i,occnum(i)
|
||||
end do
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -34,11 +32,13 @@ END_PROVIDER
|
||||
|
||||
call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
|
||||
|
||||
write(6,*) ' found occupation numbers as '
|
||||
do i=1,n_act_orb
|
||||
write(6,*) i,occ_act(i)
|
||||
end do
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' found occupation numbers as '
|
||||
do i=1,n_act_orb
|
||||
write(6,*) i,occ_act(i)
|
||||
end do
|
||||
!
|
||||
|
||||
integer :: nmx
|
||||
real*8 :: xmx
|
||||
@ -152,6 +152,7 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' transformed P0tuvx '
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -197,6 +198,7 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
|
||||
one_ints_no(j,list_act(p))=d(p)
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' transformed one_ints '
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -224,5 +226,148 @@ BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
|
||||
NatOrbsFCI(j,list_act(p))=d(p)
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' transformed orbitals '
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine trf_to_natorb()
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! save the diagonal somewhere, in inverse order
|
||||
! 4-index-transform the 2-particle density matrix over active orbitals
|
||||
! correct the bielectronic integrals
|
||||
! correct the monoelectronic integrals
|
||||
! put integrals on file, as well orbitals, and the density matrices
|
||||
!
|
||||
END_DOC
|
||||
integer :: i,j,k,l,t,u,p,q,pp
|
||||
real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb)
|
||||
|
||||
! we recalculate total energies
|
||||
write(6,*)
|
||||
write(6,*) ' recalculating energies after the transformation '
|
||||
write(6,*)
|
||||
write(6,*)
|
||||
real*8 :: e_one_all
|
||||
real*8 :: e_two_all
|
||||
integer :: ii
|
||||
integer :: jj
|
||||
integer :: t3
|
||||
integer :: tt
|
||||
integer :: u3
|
||||
integer :: uu
|
||||
integer :: v
|
||||
integer :: v3
|
||||
integer :: vv
|
||||
integer :: x
|
||||
integer :: x3
|
||||
integer :: xx
|
||||
|
||||
e_one_all=0.D0
|
||||
e_two_all=0.D0
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
e_one_all+=2.D0*one_ints_no(ii,ii)
|
||||
do j=1,n_core_orb
|
||||
jj=list_core(j)
|
||||
e_two_all+=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
|
||||
end do
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
t3=t+n_core_orb
|
||||
e_two_all += occnum(list_act(t)) * &
|
||||
(2.d0*bielec_PQxx_no(tt,tt,i,i) - bielec_PQxx_no(tt,ii,i,t3))
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
e_one_all += occnum(list_act(t))*one_ints_no(tt,tt)
|
||||
do u=1,n_act_orb
|
||||
uu=list_act(u)
|
||||
do v=1,n_act_orb
|
||||
v3=v+n_core_orb
|
||||
do x=1,n_act_orb
|
||||
x3=x+n_core_orb
|
||||
e_two_all +=P0tuvx_no(t,u,v,x)*bielec_PQxx_no(tt,uu,v3,x3)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' e_one_all = ',e_one_all
|
||||
write(6,*) ' e_two_all = ',e_two_all
|
||||
ecore =nuclear_repulsion
|
||||
ecore_bis=nuclear_repulsion
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
ecore +=2.D0*one_ints_no(ii,ii)
|
||||
ecore_bis+=2.D0*one_ints_no(ii,ii)
|
||||
do j=1,n_core_orb
|
||||
jj=list_core(j)
|
||||
ecore +=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
|
||||
ecore_bis+=2.D0*bielec_PxxQ_no(ii,i,j,jj)-bielec_PxxQ_no(ii,j,j,ii)
|
||||
end do
|
||||
end do
|
||||
eone =0.D0
|
||||
eone_bis=0.D0
|
||||
etwo =0.D0
|
||||
etwo_bis=0.D0
|
||||
etwo_ter=0.D0
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
t3=t+n_core_orb
|
||||
eone += occnum(list_act(t))*one_ints_no(tt,tt)
|
||||
eone_bis += occnum(list_act(t))*one_ints_no(tt,tt)
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
eone += occnum(list_act(t)) * &
|
||||
(2.D0*bielec_PQxx_no(tt,tt,i,i ) - bielec_PQxx_no(tt,ii,i,t3))
|
||||
eone_bis += occnum(list_act(t)) * &
|
||||
(2.D0*bielec_PxxQ_no(tt,t3,i,ii) - bielec_PxxQ_no(tt,i ,i,tt))
|
||||
end do
|
||||
do u=1,n_act_orb
|
||||
uu=list_act(u)
|
||||
u3=u+n_core_orb
|
||||
do v=1,n_act_orb
|
||||
vv=list_act(v)
|
||||
v3=v+n_core_orb
|
||||
do x=1,n_act_orb
|
||||
xx=list_act(x)
|
||||
x3=x+n_core_orb
|
||||
real*8 :: h1,h2,h3
|
||||
h1=bielec_PQxx_no(tt,uu,v3,x3)
|
||||
h2=bielec_PxxQ_no(tt,u3,v3,xx)
|
||||
h3=bielecCI_no(t,u,v,xx)
|
||||
etwo +=P0tuvx_no(t,u,v,x)*h1
|
||||
etwo_bis+=P0tuvx_no(t,u,v,x)*h2
|
||||
etwo_ter+=P0tuvx_no(t,u,v,x)*h3
|
||||
if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
|
||||
write(6,9901) t,u,v,x,h1,h2,h3
|
||||
9901 format('aie: ',4I4,3E20.12)
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
write(6,*) ' energy contributions '
|
||||
write(6,*) ' core energy = ',ecore,' using PQxx integrals '
|
||||
write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals '
|
||||
write(6,*) ' 1el energy = ',eone ,' using PQxx integrals '
|
||||
write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals '
|
||||
write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals '
|
||||
write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals '
|
||||
write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals '
|
||||
write(6,*) ' ----------------------------------------- '
|
||||
write(6,*) ' sum of all = ',eone+etwo+ecore
|
||||
write(6,*)
|
||||
SOFT_TOUCH ecore ecore_bis eone eone_bis etwo etwo_bis etwo_ter
|
||||
|
||||
end subroutine trf_to_natorb
|
||||
|
||||
|
@ -1,178 +1,222 @@
|
||||
! -*- F90 -*-
|
||||
BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Single-excitation matrix
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j
|
||||
do i=1,nMonoEx+1
|
||||
do j=1,nMonoEx+1
|
||||
SXmatrix(i,j)=0.D0
|
||||
end do
|
||||
end do
|
||||
|
||||
integer :: i,j
|
||||
do i=1,nMonoEx
|
||||
SXmatrix(1,i+1)=gradvec2(i)
|
||||
SXmatrix(1+i,1)=gradvec2(i)
|
||||
end do
|
||||
|
||||
do i=1,nMonoEx+1
|
||||
do j=1,nMonoEx+1
|
||||
SXmatrix(i,j)=0.D0
|
||||
end do
|
||||
end do
|
||||
do i=1,nMonoEx
|
||||
do j=1,nMonoEx
|
||||
SXmatrix(i+1,j+1)=hessmat2(i,j)
|
||||
SXmatrix(j+1,i+1)=hessmat2(i,j)
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=1,nMonoEx
|
||||
SXmatrix(1,i+1)=gradvec2(i)
|
||||
SXmatrix(1+i,1)=gradvec2(i)
|
||||
end do
|
||||
|
||||
do i=1,nMonoEx
|
||||
do j=1,nMonoEx
|
||||
SXmatrix(i+1,j+1)=hessmat2(i,j)
|
||||
SXmatrix(j+1,i+1)=hessmat2(i,j)
|
||||
end do
|
||||
end do
|
||||
|
||||
if (bavard) then
|
||||
do i=2,nMonoEx+1
|
||||
write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
|
||||
end do
|
||||
end if
|
||||
if (bavard) then
|
||||
do i=2,nMonoEx+1
|
||||
write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
|
||||
end do
|
||||
end if
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)]
|
||||
&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Eigenvectors/eigenvalues of the single-excitation matrix
|
||||
END_DOC
|
||||
call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
|
||||
END_PROVIDER
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [real*8, SXvector, (nMonoEx+1)]
|
||||
&BEGIN_PROVIDER [real*8, energy_improvement]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Best eigenvector of the single-excitation matrix
|
||||
END_DOC
|
||||
integer :: ierr,matz,i
|
||||
real*8 :: c0
|
||||
implicit none
|
||||
integer :: ierr,matz,i
|
||||
real*8 :: c0
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' SXdiag : lowest 5 eigenvalues '
|
||||
write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
|
||||
write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2)
|
||||
write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3)
|
||||
write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4)
|
||||
write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5)
|
||||
write(6,*)
|
||||
write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
|
||||
endif
|
||||
energy_improvement = SXeigenval(1)
|
||||
call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
|
||||
write(6,*) ' SXdiag : lowest 5 eigenvalues '
|
||||
write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
|
||||
write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2)
|
||||
write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3)
|
||||
write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4)
|
||||
write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5)
|
||||
write(6,*)
|
||||
write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
|
||||
energy_improvement = SXeigenval(1)
|
||||
|
||||
integer :: best_vector
|
||||
real*8 :: best_overlap
|
||||
best_overlap=0.D0
|
||||
do i=1,nMonoEx+1
|
||||
if (SXeigenval(i).lt.0.D0) then
|
||||
if (abs(SXeigenvec(1,i)).gt.best_overlap) then
|
||||
best_overlap=abs(SXeigenvec(1,i))
|
||||
best_vector=i
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
integer :: best_vector
|
||||
real*8 :: best_overlap
|
||||
best_overlap=0.D0
|
||||
do i=1,nMonoEx+1
|
||||
if (SXeigenval(i).lt.0.D0) then
|
||||
if (abs(SXeigenvec(1,i)).gt.best_overlap) then
|
||||
best_overlap=abs(SXeigenvec(1,i))
|
||||
best_vector=i
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
energy_improvement = SXeigenval(best_vector)
|
||||
|
||||
if (bavard) then
|
||||
write(6,*) ' SXdiag : eigenvalue for best overlap with '
|
||||
write(6,*) ' previous orbitals = ',SXeigenval(best_vector)
|
||||
write(6,*) ' weight of the 1st element ',c0
|
||||
endif
|
||||
|
||||
c0=SXeigenvec(1,best_vector)
|
||||
|
||||
do i=1,nMonoEx+1
|
||||
SXvector(i)=SXeigenvec(i,best_vector)/c0
|
||||
end do
|
||||
write(6,*) ' SXdiag : eigenvalue for best overlap with '
|
||||
write(6,*) ' previous orbitals = ',SXeigenval(best_vector)
|
||||
energy_improvement = SXeigenval(best_vector)
|
||||
|
||||
c0=SXeigenvec(1,best_vector)
|
||||
write(6,*) ' weight of the 1st element ',c0
|
||||
do i=1,nMonoEx+1
|
||||
SXvector(i)=SXeigenvec(i,best_vector)/c0
|
||||
! write(6,*) ' component No ',i,' : ',SXvector(i)
|
||||
end do
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [real*8, NewOrbs, (ao_num,mo_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Updated orbitals
|
||||
END_DOC
|
||||
integer :: i,j,ialph
|
||||
implicit none
|
||||
integer :: i,j,ialph
|
||||
|
||||
call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, &
|
||||
NatOrbsFCI, size(NatOrbsFCI,1), &
|
||||
Umat, size(Umat,1), 0.d0, &
|
||||
NewOrbs, size(NewOrbs,1))
|
||||
! form the exponential of the Orbital rotations
|
||||
call get_orbrotmat
|
||||
! form the new orbitals
|
||||
do i=1,ao_num
|
||||
do j=1,mo_num
|
||||
NewOrbs(i,j)=0.D0
|
||||
end do
|
||||
end do
|
||||
|
||||
do ialph=1,ao_num
|
||||
do i=1,mo_num
|
||||
wrkline(i)=mo_coef(ialph,i)
|
||||
end do
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
NewOrbs(ialph,i)+=Umat(i,j)*wrkline(j)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Orbital rotation matrix
|
||||
END_DOC
|
||||
integer :: i,j,indx,k,iter,t,a,ii,tt,aa
|
||||
logical :: converged
|
||||
BEGIN_PROVIDER [real*8, Tpotmat, (mo_num,mo_num) ]
|
||||
&BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
|
||||
&BEGIN_PROVIDER [real*8, wrkline, (mo_num) ]
|
||||
&BEGIN_PROVIDER [real*8, Tmat, (mo_num,mo_num) ]
|
||||
END_PROVIDER
|
||||
|
||||
real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num)
|
||||
real*8 :: Tmat(mo_num,mo_num)
|
||||
real*8 :: f
|
||||
subroutine get_orbrotmat
|
||||
implicit none
|
||||
integer :: i,j,indx,k,iter,t,a,ii,tt,aa
|
||||
real*8 :: sum
|
||||
logical :: converged
|
||||
|
||||
! the orbital rotation matrix T
|
||||
Tmat(:,:)=0.D0
|
||||
indx=1
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
indx+=1
|
||||
Tmat(ii,tt)= SXvector(indx)
|
||||
Tmat(tt,ii)=-SXvector(indx)
|
||||
end do
|
||||
end do
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
do a=1,n_virt_orb
|
||||
aa=list_virt(a)
|
||||
indx+=1
|
||||
Tmat(ii,aa)= SXvector(indx)
|
||||
Tmat(aa,ii)=-SXvector(indx)
|
||||
end do
|
||||
end do
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
do a=1,n_virt_orb
|
||||
aa=list_virt(a)
|
||||
indx+=1
|
||||
Tmat(tt,aa)= SXvector(indx)
|
||||
Tmat(aa,tt)=-SXvector(indx)
|
||||
end do
|
||||
end do
|
||||
|
||||
! Form the exponential
|
||||
! the orbital rotation matrix T
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
Tmat(i,j)=0.D0
|
||||
Umat(i,j)=0.D0
|
||||
Tpotmat(i,j)=0.D0
|
||||
end do
|
||||
Tpotmat(i,i)=1.D0
|
||||
end do
|
||||
|
||||
Tpotmat(:,:)=0.D0
|
||||
Umat(:,:) =0.D0
|
||||
do i=1,mo_num
|
||||
Tpotmat(i,i)=1.D0
|
||||
Umat(i,i) =1.d0
|
||||
end do
|
||||
iter=0
|
||||
converged=.false.
|
||||
do while (.not.converged)
|
||||
iter+=1
|
||||
f = 1.d0 / dble(iter)
|
||||
Tpotmat2(:,:) = Tpotmat(:,:) * f
|
||||
call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, &
|
||||
Tpotmat2, size(Tpotmat2,1), &
|
||||
Tmat, size(Tmat,1), 0.d0, &
|
||||
Tpotmat, size(Tpotmat,1))
|
||||
Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
|
||||
indx=1
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
indx+=1
|
||||
Tmat(ii,tt)= SXvector(indx)
|
||||
Tmat(tt,ii)=-SXvector(indx)
|
||||
end do
|
||||
end do
|
||||
do i=1,n_core_orb
|
||||
ii=list_core(i)
|
||||
do a=1,n_virt_orb
|
||||
aa=list_virt(a)
|
||||
indx+=1
|
||||
Tmat(ii,aa)= SXvector(indx)
|
||||
Tmat(aa,ii)=-SXvector(indx)
|
||||
end do
|
||||
end do
|
||||
do t=1,n_act_orb
|
||||
tt=list_act(t)
|
||||
do a=1,n_virt_orb
|
||||
aa=list_virt(a)
|
||||
indx+=1
|
||||
Tmat(tt,aa)= SXvector(indx)
|
||||
Tmat(aa,tt)=-SXvector(indx)
|
||||
end do
|
||||
end do
|
||||
|
||||
converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
|
||||
end do
|
||||
write(6,*) ' the T matrix '
|
||||
do indx=1,nMonoEx
|
||||
i=excit(1,indx)
|
||||
j=excit(2,indx)
|
||||
! if (abs(Tmat(i,j)).gt.1.D0) then
|
||||
! write(6,*) ' setting matrix element ',i,j,' of ',Tmat(i,j),' to ' &
|
||||
! , sign(1.D0,Tmat(i,j))
|
||||
! Tmat(i,j)=sign(1.D0,Tmat(i,j))
|
||||
! Tmat(j,i)=-Tmat(i,j)
|
||||
! end if
|
||||
if (abs(Tmat(i,j)).gt.1.D-9) write(6,9901) i,j,excit_class(indx),Tmat(i,j)
|
||||
9901 format(' ',i4,' -> ',i4,' (',A3,') : ',E14.6)
|
||||
end do
|
||||
|
||||
write(6,*)
|
||||
write(6,*) ' forming the matrix exponential '
|
||||
write(6,*)
|
||||
! form the exponential
|
||||
iter=0
|
||||
converged=.false.
|
||||
do while (.not.converged)
|
||||
iter+=1
|
||||
! add the next term
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
Umat(i,j)+=Tpotmat(i,j)
|
||||
end do
|
||||
end do
|
||||
! next power of T, we multiply Tpotmat with Tmat/iter
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
wrkline(j)=Tpotmat(i,j)/dble(iter)
|
||||
Tpotmat(i,j)=0.D0
|
||||
end do
|
||||
do j=1,mo_num
|
||||
do k=1,mo_num
|
||||
Tpotmat(i,j)+=wrkline(k)*Tmat(k,j)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
! Convergence test
|
||||
sum=0.D0
|
||||
do i=1,mo_num
|
||||
do j=1,mo_num
|
||||
sum+=abs(Tpotmat(i,j))
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' Iteration No ',iter,' Sum = ',sum
|
||||
if (sum.lt.1.D-6) then
|
||||
converged=.true.
|
||||
end if
|
||||
if (iter.ge.NItExpMax) then
|
||||
stop ' no convergence '
|
||||
end if
|
||||
end do
|
||||
write(6,*)
|
||||
write(6,*) ' Converged ! '
|
||||
write(6,*)
|
||||
|
||||
end subroutine get_orbrotmat
|
||||
|
||||
BEGIN_PROVIDER [integer, NItExpMax]
|
||||
NItExpMax=100
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
@ -42,6 +42,8 @@
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
write(6,*) ' e_one_all = ',e_one_all
|
||||
write(6,*) ' e_two_all = ',e_two_all
|
||||
ecore =nuclear_repulsion
|
||||
ecore_bis=nuclear_repulsion
|
||||
do i=1,n_core_orb
|
||||
@ -96,6 +98,24 @@
|
||||
end do
|
||||
end do
|
||||
|
||||
write(6,*) ' energy contributions '
|
||||
write(6,*) ' core energy = ',ecore,' using PQxx integrals '
|
||||
write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals '
|
||||
write(6,*) ' 1el energy = ',eone ,' using PQxx integrals '
|
||||
write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals '
|
||||
write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals '
|
||||
write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals '
|
||||
write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals '
|
||||
write(6,*) ' ----------------------------------------- '
|
||||
write(6,*) ' sum of all = ',eone+etwo+ecore
|
||||
write(6,*)
|
||||
write(6,*) ' nuclear (qp) = ',nuclear_repulsion
|
||||
write(6,*) ' core energy (qp) = ',core_energy
|
||||
write(6,*) ' 1el energy (qp) = ',psi_energy_h_core(1)
|
||||
write(6,*) ' 2el energy (qp) = ',psi_energy_two_e(1)
|
||||
write(6,*) ' nuc + 1 + 2 (qp) = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1)
|
||||
write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user