mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-01 05:33:40 +02:00
Compare commits
3 Commits
6531181316
...
2ef517488c
Author | SHA1 | Date | |
---|---|---|---|
2ef517488c | |||
a128c20afa | |||
5902f3231e |
@ -1,6 +1,6 @@
|
|||||||
! -*- F90 -*-
|
! -*- F90 -*-
|
||||||
BEGIN_PROVIDER [logical, bavard]
|
BEGIN_PROVIDER [logical, bavard]
|
||||||
bavard=.true.
|
bavard=.true.
|
||||||
bavard=.false.
|
! bavard=.false.
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -55,7 +55,6 @@
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
write(6,*) ' provided integrals (PQ|xx) '
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
@ -116,7 +115,6 @@ 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
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' provided integrals (Px|xQ) '
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
@ -146,6 +144,5 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' provided integrals (tu|xP) '
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -84,7 +84,6 @@
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed PQxx'
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -176,7 +175,6 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed PxxQ '
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -267,7 +265,6 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed tuvP '
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -12,24 +12,32 @@ subroutine run
|
|||||||
implicit none
|
implicit none
|
||||||
double precision :: energy_old, energy
|
double precision :: energy_old, energy
|
||||||
logical :: converged
|
logical :: converged
|
||||||
|
integer :: iteration
|
||||||
converged = .False.
|
converged = .False.
|
||||||
|
|
||||||
energy = 0.d0
|
energy = 0.d0
|
||||||
! do while (.not.converged)
|
mo_label = "MCSCF"
|
||||||
N_det = 1
|
iteration = 1
|
||||||
TOUCH N_det psi_det psi_coef
|
do while (.not.converged)
|
||||||
call run_cipsi
|
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_old = energy
|
||||||
energy = eone+etwo+ecore
|
energy = eone+etwo+ecore
|
||||||
converged = dabs(energy - energy_old) < 1.d-10
|
|
||||||
! enddo
|
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
|
||||||
|
|
||||||
end
|
end
|
||||||
|
@ -1,70 +1,19 @@
|
|||||||
use bitmasks
|
use bitmasks
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
|
BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
|
||||||
BEGIN_DOC
|
|
||||||
! 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
|
|
||||||
implicit none
|
implicit none
|
||||||
integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
|
BEGIN_DOC
|
||||||
integer :: ierr
|
! the first-order density matrix in the basis of the starting MOs.
|
||||||
integer(bit_kind) :: det_mu(N_int,2)
|
! matrix is state averaged.
|
||||||
integer(bit_kind) :: det_mu_ex(N_int,2)
|
END_DOC
|
||||||
integer(bit_kind) :: det_mu_ex1(N_int,2)
|
integer :: t,u
|
||||||
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)
|
|
||||||
|
|
||||||
write(6,*) ' providing density matrices D0 and P0 '
|
do u=1,n_act_orb
|
||||||
|
|
||||||
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
|
do t=1,n_act_orb
|
||||||
ipart=list_act(t)
|
D0tu(t,u) = one_e_dm_mo_alpha_average( list_act(t), list_act(u) ) + &
|
||||||
do u=1,n_act_orb
|
one_e_dm_mo_beta_average ( list_act(t), list_act(u) )
|
||||||
ihole=list_act(u)
|
enddo
|
||||||
! apply E_tu
|
enddo
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
@ -90,7 +39,9 @@ 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_ex1, det_mu_ex11, det_mu_ex12
|
||||||
integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
|
integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
|
||||||
|
|
||||||
write(6,*) ' providing density matrices D0 and P0 '
|
if (bavard) then
|
||||||
|
write(6,*) ' providing density matrix P0'
|
||||||
|
endif
|
||||||
|
|
||||||
P0tuvx = 0.d0
|
P0tuvx = 0.d0
|
||||||
|
|
||||||
|
@ -31,6 +31,8 @@ subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, &
|
|||||||
! get the number in the list
|
! get the number in the list
|
||||||
found=.false.
|
found=.false.
|
||||||
nu=0
|
nu=0
|
||||||
|
|
||||||
|
!TODO BOTTLENECK
|
||||||
do while (.not.found)
|
do while (.not.found)
|
||||||
nu+=1
|
nu+=1
|
||||||
if (nu.gt.N_det) then
|
if (nu.gt.N_det) then
|
||||||
@ -50,13 +52,6 @@ subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, &
|
|||||||
end do
|
end do
|
||||||
end if
|
end if
|
||||||
end do
|
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
|
end if
|
||||||
!
|
!
|
||||||
! we found the new string, the phase, and possibly the number in the list
|
! we found the new string, the phase, and possibly the number in the list
|
||||||
|
@ -1,32 +1,3 @@
|
|||||||
subroutine driver_optorb
|
subroutine driver_optorb
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i,j
|
end
|
||||||
|
|
||||||
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,7 +6,6 @@ BEGIN_PROVIDER [ integer, nMonoEx ]
|
|||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
|
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
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
|
BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
|
||||||
@ -87,9 +86,11 @@ BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
|
|||||||
norm_grad+=gradvec(indx)*gradvec(indx)
|
norm_grad+=gradvec(indx)*gradvec(indx)
|
||||||
end do
|
end do
|
||||||
norm_grad=sqrt(norm_grad)
|
norm_grad=sqrt(norm_grad)
|
||||||
write(6,*)
|
if (bavard) then
|
||||||
write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
|
write(6,*)
|
||||||
write(6,*)
|
write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
|
||||||
|
write(6,*)
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
@ -118,17 +119,11 @@ subroutine calc_grad_elem(ihole,ipart,res)
|
|||||||
call do_signed_mono_excitation(det_mu,det_mu_ex,nu &
|
call do_signed_mono_excitation(det_mu,det_mu_ex,nu &
|
||||||
,ihole,ipart,ispin,phase,ierr)
|
,ihole,ipart,ispin,phase,ierr)
|
||||||
if (ierr.eq.1) then
|
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 &
|
call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int &
|
||||||
,N_det,N_det,N_states,i_H_psi_array)
|
,N_det,N_det,N_states,i_H_psi_array)
|
||||||
do istate=1,N_states
|
do istate=1,N_states
|
||||||
res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
|
res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
|
||||||
end do
|
end do
|
||||||
! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res
|
|
||||||
end if
|
end if
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
@ -176,9 +171,11 @@ BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
|
|||||||
norm_grad+=gradvec2(indx)*gradvec2(indx)
|
norm_grad+=gradvec2(indx)*gradvec2(indx)
|
||||||
end do
|
end do
|
||||||
norm_grad=sqrt(norm_grad)
|
norm_grad=sqrt(norm_grad)
|
||||||
write(6,*)
|
if (bavard) then
|
||||||
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
|
write(6,*)
|
||||||
write(6,*)
|
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
|
||||||
|
write(6,*)
|
||||||
|
endif
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -14,8 +14,10 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
|
|||||||
character*3 :: iexc,jexc
|
character*3 :: iexc,jexc
|
||||||
real*8 :: res
|
real*8 :: res
|
||||||
|
|
||||||
write(6,*) ' providing Hessian matrix hessmat '
|
if (bavard) then
|
||||||
write(6,*) ' nMonoEx = ',nMonoEx
|
write(6,*) ' providing Hessian matrix hessmat '
|
||||||
|
write(6,*) ' nMonoEx = ',nMonoEx
|
||||||
|
endif
|
||||||
|
|
||||||
do indx=1,nMonoEx
|
do indx=1,nMonoEx
|
||||||
do jndx=1,nMonoEx
|
do jndx=1,nMonoEx
|
||||||
@ -32,8 +34,6 @@ BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)]
|
|||||||
jpart=excit(2,jndx)
|
jpart=excit(2,jndx)
|
||||||
jexc=excit_class(jndx)
|
jexc=excit_class(jndx)
|
||||||
call calc_hess_elem(ihole,ipart,jhole,jpart,res)
|
call calc_hess_elem(ihole,ipart,jhole,jpart,res)
|
||||||
! write(6,*) ' Hessian ',ihole,'->',ipart &
|
|
||||||
! ,' (',iexc,')',jhole,'->',jpart,' (',jexc,')',res
|
|
||||||
hessmat(indx,jndx)=res
|
hessmat(indx,jndx)=res
|
||||||
hessmat(jndx,indx)=res
|
hessmat(jndx,indx)=res
|
||||||
end do
|
end do
|
||||||
@ -198,8 +198,10 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
|
|||||||
real*8 :: hessmat_iatb
|
real*8 :: hessmat_iatb
|
||||||
real*8 :: hessmat_taub
|
real*8 :: hessmat_taub
|
||||||
|
|
||||||
write(6,*) ' providing Hessian matrix hessmat2 '
|
if (bavard) then
|
||||||
write(6,*) ' nMonoEx = ',nMonoEx
|
write(6,*) ' providing Hessian matrix hessmat2 '
|
||||||
|
write(6,*) ' nMonoEx = ',nMonoEx
|
||||||
|
endif
|
||||||
|
|
||||||
indx=1
|
indx=1
|
||||||
do i=1,n_core_orb
|
do i=1,n_core_orb
|
||||||
@ -214,7 +216,6 @@ BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)]
|
|||||||
do u=ustart,n_act_orb
|
do u=ustart,n_act_orb
|
||||||
hessmat2(indx,jndx)=hessmat_itju(i,t,j,u)
|
hessmat2(indx,jndx)=hessmat_itju(i,t,j,u)
|
||||||
hessmat2(jndx,indx)=hessmat2(indx,jndx)
|
hessmat2(jndx,indx)=hessmat2(indx,jndx)
|
||||||
! write(6,*) ' result I :',i,t,j,u,indx,jndx,hessmat(indx,jndx),hessmat2(indx,jndx)
|
|
||||||
jndx+=1
|
jndx+=1
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
@ -294,7 +295,6 @@ real*8 function hessmat_itju(i,t,j,u)
|
|||||||
integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
|
integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
|
||||||
real*8 :: term,t2
|
real*8 :: term,t2
|
||||||
|
|
||||||
! write(6,*) ' hessmat_itju ',i,t,j,u
|
|
||||||
ii=list_core(i)
|
ii=list_core(i)
|
||||||
tt=list_act(t)
|
tt=list_act(t)
|
||||||
if (i.eq.j) then
|
if (i.eq.j) then
|
||||||
@ -340,8 +340,6 @@ real*8 function hessmat_itju(i,t,j,u)
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
!!! write(6,*) ' direct diff ',i,t,j,u,term,term2
|
|
||||||
!!! term=term2
|
|
||||||
end if
|
end if
|
||||||
else
|
else
|
||||||
! it/ju
|
! it/ju
|
||||||
@ -382,7 +380,6 @@ real*8 function hessmat_itja(i,t,j,a)
|
|||||||
integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
|
integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
|
||||||
real*8 :: term
|
real*8 :: term
|
||||||
|
|
||||||
! write(6,*) ' hessmat_itja ',i,t,j,a
|
|
||||||
! it/ja
|
! it/ja
|
||||||
ii=list_core(i)
|
ii=list_core(i)
|
||||||
tt=list_act(t)
|
tt=list_act(t)
|
||||||
@ -416,7 +413,6 @@ 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
|
integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3
|
||||||
real*8 :: term
|
real*8 :: term
|
||||||
|
|
||||||
! write(6,*) ' hessmat_itua ',i,t,u,a
|
|
||||||
ii=list_core(i)
|
ii=list_core(i)
|
||||||
tt=list_act(t)
|
tt=list_act(t)
|
||||||
t3=t+n_core_orb
|
t3=t+n_core_orb
|
||||||
@ -457,7 +453,6 @@ real*8 function hessmat_iajb(i,a,j,b)
|
|||||||
implicit none
|
implicit none
|
||||||
integer :: i,a,j,b,ii,aa,jj,bb
|
integer :: i,a,j,b,ii,aa,jj,bb
|
||||||
real*8 :: term
|
real*8 :: term
|
||||||
! write(6,*) ' hessmat_iajb ',i,a,j,b
|
|
||||||
|
|
||||||
ii=list_core(i)
|
ii=list_core(i)
|
||||||
aa=list_virt(a)
|
aa=list_virt(a)
|
||||||
@ -495,7 +490,6 @@ real*8 function hessmat_iatb(i,a,t,b)
|
|||||||
integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
|
integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
|
||||||
real*8 :: term
|
real*8 :: term
|
||||||
|
|
||||||
! write(6,*) ' hessmat_iatb ',i,a,t,b
|
|
||||||
ii=list_core(i)
|
ii=list_core(i)
|
||||||
aa=list_virt(a)
|
aa=list_virt(a)
|
||||||
tt=list_act(t)
|
tt=list_act(t)
|
||||||
@ -552,7 +546,6 @@ real*8 function hessmat_taub(t,a,u,b)
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
term=t1+t2+t3
|
term=t1+t2+t3
|
||||||
! write(6,*) ' Hess taub ',t,a,t1,t2,t3
|
|
||||||
else
|
else
|
||||||
bb=list_virt(b)
|
bb=list_virt(b)
|
||||||
! ta/tb b/=a
|
! ta/tb b/=a
|
||||||
|
@ -14,10 +14,12 @@
|
|||||||
occnum(list_act(i))=occ_act(n_act_orb-i+1)
|
occnum(list_act(i))=occ_act(n_act_orb-i+1)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
write(6,*) ' occupation numbers '
|
if (bavard) then
|
||||||
do i=1,mo_num
|
write(6,*) ' occupation numbers '
|
||||||
write(6,*) i,occnum(i)
|
do i=1,mo_num
|
||||||
end do
|
write(6,*) i,occnum(i)
|
||||||
|
end do
|
||||||
|
endif
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -32,13 +34,11 @@ END_PROVIDER
|
|||||||
|
|
||||||
call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
|
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
|
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
|
integer :: nmx
|
||||||
real*8 :: xmx
|
real*8 :: xmx
|
||||||
@ -152,7 +152,6 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
|
|||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed P0tuvx '
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -198,7 +197,6 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
|
|||||||
one_ints_no(j,list_act(p))=d(p)
|
one_ints_no(j,list_act(p))=d(p)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed one_ints '
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
@ -226,148 +224,5 @@ BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
|
|||||||
NatOrbsFCI(j,list_act(p))=d(p)
|
NatOrbsFCI(j,list_act(p))=d(p)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed orbitals '
|
|
||||||
END_PROVIDER
|
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,222 +1,178 @@
|
|||||||
! -*- F90 -*-
|
|
||||||
BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
|
BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i,j
|
BEGIN_DOC
|
||||||
do i=1,nMonoEx+1
|
! Single-excitation matrix
|
||||||
do j=1,nMonoEx+1
|
END_DOC
|
||||||
SXmatrix(i,j)=0.D0
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
do i=1,nMonoEx
|
integer :: i,j
|
||||||
SXmatrix(1,i+1)=gradvec2(i)
|
|
||||||
SXmatrix(1+i,1)=gradvec2(i)
|
|
||||||
end do
|
|
||||||
|
|
||||||
do i=1,nMonoEx
|
do i=1,nMonoEx+1
|
||||||
do j=1,nMonoEx
|
do j=1,nMonoEx+1
|
||||||
SXmatrix(i+1,j+1)=hessmat2(i,j)
|
SXmatrix(i,j)=0.D0
|
||||||
SXmatrix(j+1,i+1)=hessmat2(i,j)
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
|
||||||
|
|
||||||
if (bavard) then
|
do i=1,nMonoEx
|
||||||
do i=2,nMonoEx+1
|
SXmatrix(1,i+1)=gradvec2(i)
|
||||||
write(6,*) ' diagonal of the Hessian : ',i,hessmat2(i,i)
|
SXmatrix(1+i,1)=gradvec2(i)
|
||||||
end do
|
end do
|
||||||
end if
|
|
||||||
|
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
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)]
|
BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,nMonoEx+1)]
|
||||||
&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)]
|
&BEGIN_PROVIDER [real*8, SXeigenval, (nMonoEx+1)]
|
||||||
END_PROVIDER
|
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
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, SXvector, (nMonoEx+1)]
|
BEGIN_PROVIDER [real*8, SXvector, (nMonoEx+1)]
|
||||||
&BEGIN_PROVIDER [real*8, energy_improvement]
|
&BEGIN_PROVIDER [real*8, energy_improvement]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: ierr,matz,i
|
BEGIN_DOC
|
||||||
real*8 :: c0
|
! Best eigenvector of the single-excitation matrix
|
||||||
|
END_DOC
|
||||||
|
integer :: ierr,matz,i
|
||||||
|
real*8 :: c0
|
||||||
|
|
||||||
call lapack_diag(SXeigenval,SXeigenvec,SXmatrix,nMonoEx+1,nMonoEx+1)
|
if (bavard) then
|
||||||
write(6,*) ' SXdiag : lowest 5 eigenvalues '
|
write(6,*) ' SXdiag : lowest 5 eigenvalues '
|
||||||
write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
|
write(6,*) ' 1 - ',SXeigenval(1),SXeigenvec(1,1)
|
||||||
write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2)
|
write(6,*) ' 2 - ',SXeigenval(2),SXeigenvec(1,2)
|
||||||
write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3)
|
write(6,*) ' 3 - ',SXeigenval(3),SXeigenvec(1,3)
|
||||||
write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4)
|
write(6,*) ' 4 - ',SXeigenval(4),SXeigenvec(1,4)
|
||||||
write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5)
|
write(6,*) ' 5 - ',SXeigenval(5),SXeigenvec(1,5)
|
||||||
write(6,*)
|
write(6,*)
|
||||||
write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
|
write(6,*) ' SXdiag : lowest eigenvalue = ',SXeigenval(1)
|
||||||
energy_improvement = SXeigenval(1)
|
endif
|
||||||
|
energy_improvement = SXeigenval(1)
|
||||||
|
|
||||||
integer :: best_vector
|
integer :: best_vector
|
||||||
real*8 :: best_overlap
|
real*8 :: best_overlap
|
||||||
best_overlap=0.D0
|
best_overlap=0.D0
|
||||||
do i=1,nMonoEx+1
|
do i=1,nMonoEx+1
|
||||||
if (SXeigenval(i).lt.0.D0) then
|
if (SXeigenval(i).lt.0.D0) then
|
||||||
if (abs(SXeigenvec(1,i)).gt.best_overlap) then
|
if (abs(SXeigenvec(1,i)).gt.best_overlap) then
|
||||||
best_overlap=abs(SXeigenvec(1,i))
|
best_overlap=abs(SXeigenvec(1,i))
|
||||||
best_vector=i
|
best_vector=i
|
||||||
end if
|
end if
|
||||||
end if
|
end if
|
||||||
end do
|
end do
|
||||||
|
|
||||||
write(6,*) ' SXdiag : eigenvalue for best overlap with '
|
energy_improvement = SXeigenval(best_vector)
|
||||||
write(6,*) ' previous orbitals = ',SXeigenval(best_vector)
|
|
||||||
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
|
||||||
|
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, NewOrbs, (ao_num,mo_num) ]
|
BEGIN_PROVIDER [real*8, NewOrbs, (ao_num,mo_num) ]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i,j,ialph
|
BEGIN_DOC
|
||||||
|
! Updated orbitals
|
||||||
|
END_DOC
|
||||||
|
integer :: i,j,ialph
|
||||||
|
|
||||||
! form the exponential of the Orbital rotations
|
call dgemm('N','T', ao_num,mo_num,mo_num,1.d0, &
|
||||||
call get_orbrotmat
|
NatOrbsFCI, size(NatOrbsFCI,1), &
|
||||||
! form the new orbitals
|
Umat, size(Umat,1), 0.d0, &
|
||||||
do i=1,ao_num
|
NewOrbs, size(NewOrbs,1))
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, Tpotmat, (mo_num,mo_num) ]
|
BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
|
||||||
&BEGIN_PROVIDER [real*8, Umat, (mo_num,mo_num) ]
|
implicit none
|
||||||
&BEGIN_PROVIDER [real*8, wrkline, (mo_num) ]
|
BEGIN_DOC
|
||||||
&BEGIN_PROVIDER [real*8, Tmat, (mo_num,mo_num) ]
|
! Orbital rotation matrix
|
||||||
END_PROVIDER
|
END_DOC
|
||||||
|
integer :: i,j,indx,k,iter,t,a,ii,tt,aa
|
||||||
|
logical :: converged
|
||||||
|
|
||||||
subroutine get_orbrotmat
|
real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num)
|
||||||
implicit none
|
real*8 :: Tmat(mo_num,mo_num)
|
||||||
integer :: i,j,indx,k,iter,t,a,ii,tt,aa
|
real*8 :: f
|
||||||
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
|
||||||
|
|
||||||
! the orbital rotation matrix T
|
! Form the exponential
|
||||||
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
|
|
||||||
|
|
||||||
indx=1
|
Tpotmat(:,:)=0.D0
|
||||||
do i=1,n_core_orb
|
Umat(:,:) =0.D0
|
||||||
ii=list_core(i)
|
do i=1,mo_num
|
||||||
do t=1,n_act_orb
|
Tpotmat(i,i)=1.D0
|
||||||
tt=list_act(t)
|
Umat(i,i) =1.d0
|
||||||
indx+=1
|
end do
|
||||||
Tmat(ii,tt)= SXvector(indx)
|
iter=0
|
||||||
Tmat(tt,ii)=-SXvector(indx)
|
converged=.false.
|
||||||
end do
|
do while (.not.converged)
|
||||||
end do
|
iter+=1
|
||||||
do i=1,n_core_orb
|
f = 1.d0 / dble(iter)
|
||||||
ii=list_core(i)
|
Tpotmat2(:,:) = Tpotmat(:,:) * f
|
||||||
do a=1,n_virt_orb
|
call dgemm('N','N', mo_num,mo_num,mo_num,1.d0, &
|
||||||
aa=list_virt(a)
|
Tpotmat2, size(Tpotmat2,1), &
|
||||||
indx+=1
|
Tmat, size(Tmat,1), 0.d0, &
|
||||||
Tmat(ii,aa)= SXvector(indx)
|
Tpotmat, size(Tpotmat,1))
|
||||||
Tmat(aa,ii)=-SXvector(indx)
|
Umat(:,:) = Umat(:,:) + Tpotmat(:,:)
|
||||||
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
|
|
||||||
|
|
||||||
write(6,*) ' the T matrix '
|
converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
|
||||||
do indx=1,nMonoEx
|
end do
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
@ -42,8 +42,6 @@
|
|||||||
end do
|
end do
|
||||||
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 =nuclear_repulsion
|
||||||
ecore_bis=nuclear_repulsion
|
ecore_bis=nuclear_repulsion
|
||||||
do i=1,n_core_orb
|
do i=1,n_core_orb
|
||||||
@ -98,24 +96,6 @@
|
|||||||
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,*)
|
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user