9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-04 02:48:24 +01:00

Cleaned neworbs

This commit is contained in:
Anthony Scemama 2019-06-25 23:10:19 +02:00
parent 6531181316
commit 5902f3231e
2 changed files with 155 additions and 209 deletions

View File

@ -21,10 +21,6 @@ subroutine run
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

View File

@ -1,222 +1,172 @@
! -*- F90 -*-
BEGIN_PROVIDER [real*8, SXmatrix, (nMonoEx+1,nMonoEx+1)]
implicit none
integer :: i,j
do i=1,nMonoEx+1
do j=1,nMonoEx+1
SXmatrix(i,j)=0.D0
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
implicit none
BEGIN_DOC
! Single-excitation matrix
END_DOC
integer :: i,j
do i=1,nMonoEx+1
do j=1,nMonoEx+1
SXmatrix(i,j)=0.D0
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
END_PROVIDER
BEGIN_PROVIDER [real*8, SXeigenvec, (nMonoEx+1,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, energy_improvement]
implicit none
integer :: ierr,matz,i
real*8 :: c0
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
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
implicit none
BEGIN_DOC
! Best eigenvector of the single-excitation matrix
END_DOC
integer :: ierr,matz,i
real*8 :: c0
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
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
integer :: i,j,ialph
! 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
implicit none
BEGIN_DOC
! Updated orbitals
END_DOC
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))
END_PROVIDER
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
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
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
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
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
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
real*8 :: Tpotmat (mo_num,mo_num), Tpotmat2 (mo_num,mo_num)
real*8 :: Tmat(mo_num,mo_num)
real*8 :: f
! 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
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(:,:)
converged = ( sum(abs(Tpotmat(:,:))) < 1.d-6).or.(iter>30)
end do
END_PROVIDER