mirror of
https://code.it4i.cz/sccs/easyconfigs-it4i.git
synced 2025-04-16 19:50:50 +01:00
2074 lines
76 KiB
Fortran
2074 lines
76 KiB
Fortran
************************************************************************
|
|
subroutine pml(nmatrix,nbasis,natoms,natrange,s,c,v,work,lsym)
|
|
************************************************************************
|
|
* Pipek-Mezey localization
|
|
************************************************************************
|
|
implicit none
|
|
integer, parameter :: maxit = 500
|
|
integer nmatrix,nbasis,natoms,iatoms,natrange(2,natoms),n1,n2
|
|
real*8 norm,s(nbasis,nbasis),c(nbasis,nmatrix),v(nbasis,nmatrix)
|
|
real*8 sina,cosa,a,b,angle,qii,qij,ci,cj,vi,vj,denom,work,tol
|
|
integer i,j,mu,it
|
|
logical lsym,ldorot
|
|
character*4 scftol
|
|
c write(6,*) (s(i,i),i=1,nbasis)
|
|
C
|
|
c write(6,*) 'Beg'
|
|
c do i=1,nmatrix
|
|
c write(6,"(10000f10.6)") (c(i,j),j=1,nbasis)
|
|
c enddo
|
|
C Read variables
|
|
call getvar('natrange ',natrange)
|
|
call getkey('scftol',6,scftol,4)
|
|
read(scftol,*) i
|
|
tol=min(1.d-12,10.d0**(-i))
|
|
C Read S matrix
|
|
c rewind(ifile)
|
|
c call roeint(work,work,s,ifile,nbasis)
|
|
C V = S * C
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,c,nbasis,0.d0,v,
|
|
$nbasis)
|
|
C PM iteration
|
|
c call timer
|
|
write(*,"(' Step Convergence')")
|
|
do it=1,maxit
|
|
ldorot=lsym.and.it.eq.1
|
|
norm=0.d0
|
|
C$OMP PARALLEL DO Schedule(Dynamic)
|
|
C$OMP& DEFAULT(PRIVATE)
|
|
C$OMP& SHARED(natoms,natrange,c,v,ldorot,nbasis,nmatrix,tol)
|
|
C$OMP& REDUCTION(+:norm)
|
|
do i=2,nmatrix
|
|
do j=1,i-1
|
|
C Calculate Mulliken charges and parameters A and B
|
|
a=0.d0
|
|
b=0.d0
|
|
do iatoms=1,natoms
|
|
n1=natrange(1,iatoms)+1
|
|
n2=natrange(2,iatoms)
|
|
qii=dot_product(c(n1:n2,i),v(n1:n2,i))-
|
|
$ dot_product(c(n1:n2,j),v(n1:n2,j))
|
|
qij=dot_product(c(n1:n2,i),v(n1:n2,j))+
|
|
$ dot_product(c(n1:n2,j),v(n1:n2,i))
|
|
qij=0.5d0*qij
|
|
a=a+qij**2-0.25d0*qii**2
|
|
b=b+qij*qii
|
|
enddo
|
|
C Calculate rotation angle
|
|
denom=a*a+b*b
|
|
if(denom.gt.0.1d0*tol) then
|
|
angle=sign(0.25d0*dacos(-a/dsqrt(denom)),b)
|
|
if(ldorot.and.dabs(angle).lt.1.d-4) angle=0.1d0
|
|
else
|
|
angle=0.d0
|
|
endif
|
|
C 2x2 rotation
|
|
sina=dsin(angle)
|
|
cosa=dcos(angle)
|
|
C$OMP CRITICAL (cv)
|
|
do mu=1,nbasis
|
|
ci= cosa*c(mu,i)+sina*c(mu,j)
|
|
cj=-sina*c(mu,i)+cosa*c(mu,j)
|
|
c(mu,i)=ci
|
|
c(mu,j)=cj
|
|
vi= cosa*v(mu,i)+sina*v(mu,j)
|
|
vj=-sina*v(mu,i)+cosa*v(mu,j)
|
|
v(mu,i)=vi
|
|
v(mu,j)=vj
|
|
enddo
|
|
C$OMP END CRITICAL (cv)
|
|
norm=norm+1.d0-cosa**2
|
|
enddo
|
|
enddo
|
|
C$OMP END PARALLEL DO
|
|
c do j=1,nmatrix
|
|
c write(6,"(i3,1000f10.6)") j,(c(i,j),i=1,nbasis)
|
|
c enddo
|
|
write(*,"(i4,f17.12)") it,norm
|
|
if(norm.lt.tol) exit
|
|
enddo
|
|
c call timer
|
|
C
|
|
c write(6,*) 'End'
|
|
c do i=1,nbasis
|
|
c write(6,"(10000f10.6)") (c(i,j),j=1,nmatrix)
|
|
c enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine pmll(nmatrix,nbasis,natoms,natrange,s,c,v,work,lsym,
|
|
$ifile,step,devparr,iout)
|
|
************************************************************************
|
|
* Pipek-Mezey localization with Lowdin charges
|
|
************************************************************************
|
|
implicit none
|
|
integer, parameter :: maxit = 500
|
|
integer nmatrix,nbasis,natoms,iatoms,natrange(2,natoms),iout
|
|
real*8 norm,s(nbasis,nbasis),c(nbasis,nmatrix),v(nbasis,nmatrix)
|
|
real*8 sina,cosa,a,b,angle,qii,qjj,qij,ci,cj,vi,vj,denom,work,tol
|
|
real*8 devparr
|
|
integer i,j,mu,it,ifile,n1,n2,step
|
|
logical lsym,ldorot
|
|
character*4 scftol
|
|
C Read variables
|
|
call getvar('natrange ',natrange)
|
|
call getkey('scftol',6,scftol,4)
|
|
read(scftol,*) i
|
|
tol=min(max(10.d0**(-step),1.d-12),10.d0**(-i),devparr)
|
|
if(devparr.lt.10.d0**(-dble(i)/2.d0)) tol=min(tol,1.d-12)
|
|
C Read S^1/2
|
|
open(ifile,file='SROOT',form='UNFORMATTED')
|
|
call roeint(work,work,s,ifile,nbasis)
|
|
C V = S^1/2 * C
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,c,nbasis,0.d0,v,
|
|
$nbasis)
|
|
C PM iteration
|
|
c call timer
|
|
write(iout,"(' Step Convergence')")
|
|
do it=1,maxit
|
|
ldorot=lsym.and.it.eq.1
|
|
norm=0.d0
|
|
C$OMP PARALLEL DO Schedule(Dynamic)
|
|
C$OMP& DEFAULT(PRIVATE)
|
|
C$OMP& SHARED(natoms,natrange,v,ldorot,nbasis,nmatrix,tol)
|
|
C$OMP& REDUCTION(+:norm)
|
|
do i=2,nmatrix
|
|
do j=1,i-1
|
|
C Calculate Mulliken charges and parameters A and B
|
|
a=0.d0
|
|
b=0.d0
|
|
do iatoms=1,natoms
|
|
n1=natrange(1,iatoms)+1
|
|
n2=natrange(2,iatoms)
|
|
C Exponent=2
|
|
qii=dot_product(v(n1:n2,i),v(n1:n2,i))-
|
|
$ dot_product(v(n1:n2,j),v(n1:n2,j))
|
|
qij=dot_product(v(n1:n2,i),v(n1:n2,j))
|
|
a=a+qij*qij-0.25d0*qii*qii
|
|
b=b+qij*qii
|
|
C Exponent=4
|
|
c qii=dot_product(v(n1:n2,i),v(n1:n2,i))
|
|
c qjj=dot_product(v(n1:n2,j),v(n1:n2,j))
|
|
c qij=dot_product(v(n1:n2,i),v(n1:n2,j))
|
|
c a=a-qii**4-qjj**4+6.d0*(qii**2+qjj**2)*qij**2+
|
|
c $qii**3*qjj+qii*qjj**3
|
|
c b=b+4.d0*qij*(qii**3-qjj**3)
|
|
enddo
|
|
C Calculate rotation angle
|
|
denom=a*a+b*b
|
|
if(denom.gt.0.1d0*tol) then
|
|
angle=sign(0.25d0*dacos(-a/dsqrt(denom)),b)
|
|
if(ldorot.and.dabs(angle).lt.1.d-4) angle=0.1d0
|
|
else
|
|
angle=0.d0
|
|
endif
|
|
C 2x2 rotation
|
|
sina=dsin(angle)
|
|
cosa=dcos(angle)
|
|
C$OMP CRITICAL (cpmll)
|
|
do mu=1,nbasis
|
|
vi= cosa*v(mu,i)+sina*v(mu,j)
|
|
vj=-sina*v(mu,i)+cosa*v(mu,j)
|
|
v(mu,i)=vi
|
|
v(mu,j)=vj
|
|
enddo
|
|
C$OMP END CRITICAL (cpmll)
|
|
norm=norm+1.d0-cosa**2
|
|
enddo
|
|
enddo
|
|
C$OMP END PARALLEL DO
|
|
write(iout,"(i4,f17.12)") it,norm
|
|
if(norm.lt.tol) exit
|
|
enddo
|
|
C Read S^-1/2
|
|
call roeint(work,work,s,ifile,nbasis)
|
|
close(ifile)
|
|
C C = S^1/2 * V
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,v,nbasis,0.d0,c,
|
|
$nbasis)
|
|
c call timer
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine iao(nb,mb,no,s1,s1is,s2is,s12,c,ct,w,o,ot,p12,d)
|
|
************************************************************************
|
|
* Calculate IAO atomic charges
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,no,mb,i,j
|
|
real*8 s1(nb,nb),s1is(nb,nb),s2is(mb,mb),s12(nb,mb),c(nb,no),w(*)
|
|
real*8 ct(nb,no),o(nb,nb),ot(nb,nb),p12(nb,mb),d(nb,nb)
|
|
C
|
|
call dgemm('t','n',mb,no,nb,1.d0,s12,nb,c,nb,0.d0,w,mb)
|
|
call dsymm('l','u',mb,no,1.d0,s2is,mb,w,mb,0.d0,w(mb*nb+1),mb)
|
|
call dsymm('l','u',mb,no,1.d0,s2is,mb,w(mb*nb+1),mb,0.d0,w,mb)
|
|
call dgemm('n','n',nb,no,mb,1.d0,s12,nb,w,mb,0.d0,w(nb*nb+1),nb)
|
|
call dsymm('l','u',nb,no,1.d0,s1is,nb,w(nb*nb+1),nb,0.d0,w,nb)
|
|
call dsymm('l','u',nb,no,1.d0,s1is,nb,w,nb,0.d0,ct,nb)
|
|
call orth(nb,no,s1,ct,w)
|
|
call dsyrk('u','n',nb,no,1.d0,c,nb,0.d0,w,nb)
|
|
call dsymm('l','u',nb,nb,1.d0,w,nb,s1,nb,0.d0,o,nb)
|
|
call dsyrk('u','n',nb,no,1.d0,ct,nb,0.d0,w,nb)
|
|
call dsymm('l','u',nb,nb,1.d0,w,nb,s1,nb,0.d0,ot,nb)
|
|
call dsymm('l','u',nb,mb,1.d0,s1is,nb,s12,nb,0.d0,w,nb)
|
|
call dsymm('l','u',nb,mb,1.d0,s1is,nb,w,nb,0.d0,p12,nb)
|
|
call dgemm('n','n',nb,nb,nb,1.d0,o,nb,ot,nb,0.d0,w,nb)
|
|
call dscal(nb*nb,-1.d0,o ,1)
|
|
call dscal(nb*nb,-1.d0,ot,1)
|
|
do i=1,nb
|
|
o(i,i) =1.d0+o (i,i)
|
|
ot(i,i)=1.d0+ot(i,i)
|
|
enddo
|
|
call dgemm('n','n',nb,nb,nb,1.d0,o,nb,ot,nb,1.d0,w,nb)
|
|
call dgemm('n','n',nb,mb,nb,1.d0,w,nb,p12,nb,0.d0,s12,nb)
|
|
call orth(nb,mb,s1,s12,w)
|
|
call dsymm('l','u',nb,mb,1.d0,s1,nb,s12,nb,0.d0,w,nb)
|
|
call dsymm('l','u',nb,mb,1.d0,d ,nb,w ,nb,0.d0,w(nb*nb+1),nb)
|
|
call dsymm('l','u',nb,mb,1.d0,s1,nb,w(nb*nb+1),nb,0.d0,s2is,nb)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine ibo(nmatrix,nbasis,natoms,natrange,s,c,work,lsym,am,
|
|
$mnbasis)
|
|
************************************************************************
|
|
* Intrinsic bond orbitals (IBOs)
|
|
************************************************************************
|
|
implicit none
|
|
integer nmatrix,nbasis,natoms,iatoms,natrange(2,natoms,4),mnbasis
|
|
real*8 norm,s(nbasis,nbasis),c(mnbasis,nmatrix)
|
|
real*8 ciao(mnbasis,nmatrix)
|
|
real*8 sina,cosa,a,b,angle,qii,qjj,qij,ci,cj,vi,vj,denom,work,tol
|
|
real*8 am(nbasis,mnbasis)
|
|
integer i,j,mu,it
|
|
logical lsym,ldorot
|
|
character*4 scftol
|
|
c write(6,*) 'Beg'
|
|
c do i=1,mnbasis
|
|
c write(6,"(10000f10.5)") (c(i,j),j=1,nmatrix)
|
|
c enddo
|
|
C Read variables
|
|
call getvar('natrange ',natrange)
|
|
call getkey('scftol',6,scftol,4)
|
|
read(scftol,*) i
|
|
tol=min(1.d-10,10.d0**(-i))
|
|
C Transform MO coeffs to IAO basis
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,c,nbasis,0.d0,
|
|
$work,nbasis)
|
|
call dgemm('t','n',mnbasis,nmatrix,nbasis,1.d0,am,nbasis,work,
|
|
$nbasis,0.d0,c,mnbasis)
|
|
c do i=1,mnbasis
|
|
c write(6,"(10000f9.5)") (c(i,j),j=1,nmatrix)
|
|
c enddo
|
|
C PM iteration
|
|
write(*,"(' Step Convergence')")
|
|
do it=1,300
|
|
ldorot=lsym.and.it.eq.1
|
|
norm=0.d0
|
|
C$OMP PARALLEL DO Schedule(Dynamic)
|
|
C$OMP& DEFAULT(PRIVATE)
|
|
C$OMP& SHARED(natoms,natrange,c,ldorot,mnbasis,nmatrix,tol)
|
|
C$OMP& REDUCTION(+:norm)
|
|
do i=2,nmatrix
|
|
do j=1,i-1
|
|
C Calculate Mulliken charges and parameters A and B
|
|
a=0.d0
|
|
b=0.d0
|
|
do iatoms=1,natoms
|
|
C Exponent=2
|
|
c qii=0.d0
|
|
c qij=0.d0
|
|
c do mu=natrange(1,iatoms,4)+1,natrange(2,iatoms,4)
|
|
c qii=qii+c(mu,i)*c(mu,i)-c(mu,j)*c(mu,j)
|
|
c qij=qij+c(mu,i)*c(mu,j)
|
|
c enddo
|
|
c a=a+qij**2-0.25d0*qii**2
|
|
c b=b+qij*qii
|
|
C Exponent=4
|
|
qii=0.d0
|
|
qjj=0.d0
|
|
qij=0.d0
|
|
do mu=natrange(1,iatoms,4)+1,natrange(2,iatoms,4)
|
|
qii=qii+c(mu,i)*c(mu,i)
|
|
qjj=qjj+c(mu,j)*c(mu,j)
|
|
qij=qij+c(mu,i)*c(mu,j)
|
|
enddo
|
|
a=a-qii**4-qjj**4+6.d0*(qii**2+qjj**2)*qij**2+
|
|
$qii**3*qjj+qii*qjj**3
|
|
b=b+4.d0*qij*(qii**3-qjj**3)
|
|
enddo
|
|
C Calculate rotation angle
|
|
denom=a*a+b*b
|
|
if(denom.gt.0.1d0*tol) then
|
|
angle=sign(0.25d0*dacos(-a/dsqrt(denom)),b)
|
|
if(ldorot.and.dabs(angle).lt.1.d-4) angle=0.1d0
|
|
else
|
|
angle=0.d0
|
|
endif
|
|
C 2x2 rotation
|
|
sina=dsin(angle)
|
|
cosa=dcos(angle)
|
|
C$OMP CRITICAL (cibo)
|
|
do mu=1,mnbasis
|
|
ci= cosa*c(mu,i)+sina*c(mu,j)
|
|
cj=-sina*c(mu,i)+cosa*c(mu,j)
|
|
c(mu,i)=ci
|
|
c(mu,j)=cj
|
|
enddo
|
|
C$OMP END CRITICAL (cibo)
|
|
norm=norm+1.d0-cosa**2
|
|
enddo
|
|
enddo
|
|
C$OMP END PARALLEL DO
|
|
c do j=1,nmatrix
|
|
c write(6,"(i3,1000f10.6)") j,(c(i,j),i=1,nbasis)
|
|
c enddo
|
|
write(*,"(i4,f17.12)") it,norm
|
|
if(norm.lt.tol) exit
|
|
enddo
|
|
C
|
|
c do i=1,nmatrix
|
|
c do iatoms=1,natoms
|
|
c qii=0.d0
|
|
c do mu=natrange(1,iatoms,4)+1,natrange(2,iatoms,4)
|
|
c qii=qii+c(mu,i)*c(mu,i)
|
|
c enddo
|
|
c if(qii.gt.0.02d0) write(*,"(2i4,f8.3)") i,iatoms,2.d0*qii
|
|
c enddo
|
|
c enddo
|
|
c do i=1,mnbasis
|
|
c write(6,"(10000f9.5)") (c(i,j),j=1,nmatrix)
|
|
c enddo
|
|
C Transform MO coeffs back to AO basis
|
|
call dgemm('n','n',nbasis,nmatrix,mnbasis,1.d0,am,nbasis,c,
|
|
$mnbasis,0.d0,work,nbasis)
|
|
call dcopy(nbasis*nmatrix,work,1,c,1)
|
|
C
|
|
c write(6,*) 'End'
|
|
c do i=1,nbasis
|
|
c write(6,"(10000f9.5)") (c(i,j),j=1,nmatrix)
|
|
c enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine pmlg(nmatrix,nbasis,natoms,natrange,s,c,v,work,lsym,nc)
|
|
************************************************************************
|
|
* Pipek-Mezey localization
|
|
************************************************************************
|
|
implicit none
|
|
integer nmatrix,nbasis,natoms,iatoms,natrange(2,natoms),nc
|
|
real*8 norm,s(nbasis,nbasis),c(nbasis,nmatrix),v(nbasis,nmatrix)
|
|
real*8 sina,cosa,a,b,angle,qii,qij,ci,cj,vi,vj,denom,work,tol
|
|
integer i,j,mu,it
|
|
logical lsym,ldorot
|
|
character*4 scftol
|
|
c write(6,*) (s(i,i),i=1,nbasis)
|
|
C
|
|
c write(6,*) 'Beg'
|
|
c do i=1,nmatrix
|
|
c write(6,"(10000f10.6)") (c(i,j),j=1,nbasis)
|
|
c enddo
|
|
C Read variables
|
|
call getvar('natrange ',natrange)
|
|
call getkey('scftol',6,scftol,4)
|
|
read(scftol,*) i
|
|
tol=min(1.d-12,10.d0**(-i))
|
|
C Read S matrix
|
|
c rewind(ifile)
|
|
c call roeint(work,work,s,ifile,nbasis)
|
|
C V = S * C
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,c,nbasis,0.d0,v,
|
|
$nbasis)
|
|
C PM iteration
|
|
write(*,"(' Step Convergence')")
|
|
do it=1,300
|
|
ldorot=lsym.and.it.eq.1
|
|
norm=0.d0
|
|
do i=nc+1,nmatrix
|
|
do j=1,i-1
|
|
C Calculate Mulliken charges and parameters A and B
|
|
a=0.d0
|
|
b=0.d0
|
|
do iatoms=8,natoms
|
|
qii=0.d0
|
|
qij=0.d0
|
|
do mu=natrange(1,iatoms)+1,natrange(2,iatoms)
|
|
qii=qii+c(mu,i)*v(mu,i)-c(mu,j)*v(mu,j)
|
|
qij=qij+c(mu,i)*v(mu,j)+c(mu,j)*v(mu,i)
|
|
enddo
|
|
qij=0.5d0*qij
|
|
a=a+qij**2-0.25d0*qii**2
|
|
b=b+qij*qii
|
|
enddo
|
|
C Calculate rotation angle
|
|
denom=a*a+b*b
|
|
if(denom.gt.0.1d0*tol) then
|
|
angle=sign(0.25d0*dacos(-a/dsqrt(denom)),b)
|
|
if(ldorot.and.dabs(angle).lt.1.d-4) angle=0.1d0
|
|
else
|
|
angle=0.d0
|
|
endif
|
|
C 2x2 rotation
|
|
sina=dsin(angle)
|
|
cosa=dcos(angle)
|
|
do mu=1,nbasis
|
|
ci= cosa*c(mu,i)+sina*c(mu,j)
|
|
cj=-sina*c(mu,i)+cosa*c(mu,j)
|
|
vi= cosa*v(mu,i)+sina*v(mu,j)
|
|
vj=-sina*v(mu,i)+cosa*v(mu,j)
|
|
c(mu,i)=ci
|
|
c(mu,j)=cj
|
|
v(mu,i)=vi
|
|
v(mu,j)=vj
|
|
enddo
|
|
norm=norm+1.d0-cosa**2
|
|
enddo
|
|
enddo
|
|
c do j=1,nmatrix
|
|
c write(6,"(i3,1000f10.6)") j,(c(i,j),i=1,nbasis)
|
|
c enddo
|
|
write(*,"(i4,f17.12)") it,norm
|
|
if(norm.lt.tol) exit
|
|
enddo
|
|
C
|
|
c write(6,*) 'End'
|
|
c do i=1,nmatrix
|
|
c write(6,"(10000f10.6)") (c(i,j),j=1,nbasis)
|
|
c enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine pmlgrad(nmatrix,nbasis,natoms,natrange,s,c,v,work,lsym)
|
|
************************************************************************
|
|
* Pipek-Mezey localization with gradient
|
|
************************************************************************
|
|
implicit none
|
|
integer nmatrix,nbasis,natoms,iatoms,natrange(2,natoms)
|
|
real*8 norm,s(nbasis,nbasis),c(nbasis,nmatrix),v(nbasis,nmatrix)
|
|
real*8 sina,cosa,b,angle,qii,qij,ci,cj,vi,vj,denom,work,tol
|
|
real*8 q(nmatrix,nmatrix),gr(nmatrix,nmatrix),dnrm2
|
|
real*8 u(nmatrix,nmatrix),hs(nmatrix,nmatrix)
|
|
integer i,j,mu,it
|
|
logical lsym,ldorot
|
|
character*4 scftol
|
|
C Read variables
|
|
call getvar('natrange ',natrange)
|
|
call getkey('scftol',6,scftol,4)
|
|
read(scftol,*) i
|
|
tol=min(1.d-12,10.d0**(-i))
|
|
C V = S * C
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,c,nbasis,0.d0,v,
|
|
$nbasis)
|
|
C PM iteration
|
|
write(*,"(' Step Convergence')")
|
|
do it=1,300
|
|
ldorot=lsym.and.it.eq.1
|
|
gr=0.d0
|
|
hs=0.d0
|
|
do iatoms=1,natoms
|
|
C Calculate Q_ij^A
|
|
call dsyr2k('l','t',nmatrix,natrange(2,iatoms)-
|
|
$natrange(1,iatoms),0.5d0,c(natrange(1,iatoms)+1,1),nbasis,
|
|
$v(natrange(1,iatoms)+1,1),nbasis,0.d0,q,nmatrix)
|
|
C Calculate gradient and diagonal Hessian
|
|
do i=1,nmatrix
|
|
do j=1,i-1
|
|
gr(i,j)=gr(i,j)+(q(i,i)-q(j,j))*q(i,j)
|
|
hs(i,j)=hs(i,j)+(q(i,i)-q(j,j))*q(j,j)-
|
|
$ (q(i,i)-q(j,j))*q(i,i)+2.d0*q(i,j)**2
|
|
enddo
|
|
enddo
|
|
enddo
|
|
do i=1,nmatrix
|
|
do j=1,i-1
|
|
if(dabs(hs(i,j)).gt.0.1d0*tol) then
|
|
gr(i,j)=gr(i,j)/hs(i,j)
|
|
else
|
|
gr(i,j)=-0.01d0*gr(i,j)
|
|
endif
|
|
if(dabs(gr(i,j)).gt.1.d0) gr(i,j)=sign(0.1d0,gr(i,j))
|
|
gr(j,i)=-gr(i,j)
|
|
enddo
|
|
enddo
|
|
norm=dnrm2(nmatrix**2,gr,1)
|
|
if(ldorot) then
|
|
do i=1,nmatrix
|
|
do j=1,i-1
|
|
if(dabs(gr(i,j)).lt.tol) then
|
|
gr(i,j)=gr(i,j)+0.1d0
|
|
gr(j,i)=-gr(i,j)
|
|
endif
|
|
enddo
|
|
enddo
|
|
endif
|
|
call exp_k(nmatrix,tol,u,gr,work)
|
|
call dgemm('n','n',nbasis,nmatrix,nmatrix,1.d0,c,nbasis,
|
|
$u,nmatrix,0.d0,work,nbasis)
|
|
call dcopy(nbasis*nmatrix,work,1,c,1)
|
|
write(*,"(i4,f17.12)") it,norm
|
|
if(norm.lt.tol) exit
|
|
call dgemm('n','n',nbasis,nmatrix,nmatrix,1.d0,v,nbasis,
|
|
$u,nmatrix,0.d0,work,nbasis)
|
|
call dcopy(nbasis*nmatrix,work,1,v,1)
|
|
enddo
|
|
C
|
|
write(6,*) 'End'
|
|
do i=1,nbasis
|
|
write(6,"(10000f10.6)") (c(i,j),j=1,nmatrix)
|
|
enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
C***********************************************************************
|
|
subroutine exp_k(nmatrix,itol,umatrix,kappa,work)
|
|
C
|
|
implicit none
|
|
C
|
|
C Size of matrices
|
|
integer nmatrix
|
|
C
|
|
C Treshold for the Taylor-series
|
|
real*8 itol
|
|
C
|
|
C Length of the Taylor-series of exp(kappa)
|
|
integer nseries
|
|
C
|
|
C Full umatrix=exp(kappa) matrix
|
|
real*8 umatrix(nmatrix,*)
|
|
C
|
|
C Work array for full kappa matrix
|
|
real*8 kappa
|
|
C
|
|
C Work array for kpower matrix and for dgemm subroutine
|
|
real*8 work(*)
|
|
C
|
|
real*8 dnrm2
|
|
real*8 norm,nfactor,nn
|
|
C
|
|
integer i
|
|
C
|
|
C
|
|
C Initialize umatrix=exp(kappa) (unit matrices)
|
|
call dfillzero(umatrix,nmatrix**2)
|
|
do i=1,nmatrix
|
|
umatrix(i,i)=1.d0
|
|
enddo
|
|
C
|
|
C work(1)=kpower, it contains the powers od the full kappa matrix
|
|
C kpower=kappa, umatrix=1+kappa
|
|
call dcopy(nmatrix**2,kappa,1,work(1),1)
|
|
call daxpy(nmatrix**2,1.d0,work(1),1,umatrix,1)
|
|
C
|
|
nseries=32
|
|
norm=2.d0*itol
|
|
nfactor=1.d0
|
|
nn=2.d0
|
|
C
|
|
C Loop over nseries
|
|
do while (nn.le.nseries+1 .and. itol.le.norm)
|
|
nfactor=nfactor/nn
|
|
C
|
|
C kpower=kpower*kappa, umatrix=1/n!*kpower+umatrix
|
|
call dgemm('n','n',nmatrix,nmatrix,nmatrix,1.d0,work(1),nmatrix,
|
|
& kappa,nmatrix,0.d0,work(nmatrix**2+1),nmatrix)
|
|
call dcopy(nmatrix**2,work(nmatrix**2+1),1,work(1),1)
|
|
call daxpy(nmatrix**2,nfactor,work(1),1,umatrix,1)
|
|
C
|
|
norm=dnrm2(nmatrix**2,work(1),1)
|
|
norm=norm*nfactor
|
|
C
|
|
nn=nn+1.d0
|
|
enddo
|
|
C
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine pmlvec(nmatrix,nbasis,natoms,natrange,s,c,v,work,lsym)
|
|
************************************************************************
|
|
* Pipek-Mezey localization
|
|
************************************************************************
|
|
implicit none
|
|
integer nmatrix,nbasis,natoms,iatoms,natrange(2,natoms),k
|
|
real*8 norm,s(nbasis,nbasis),c(nbasis,nmatrix),v(nbasis,nmatrix)
|
|
real*8 sina,cosa,a,b,angle,qii,qij,ci,cj,vi,vj,denom,work,tol
|
|
real*8 q(nmatrix,nmatrix),am(nmatrix,nmatrix),bm(nmatrix,nmatrix)
|
|
real*8 u(nmatrix,nmatrix)
|
|
integer i,j,mu,it
|
|
logical lsym,ldorot
|
|
character*4 scftol
|
|
c write(6,*) (s(i,i),i=1,nbasis)
|
|
C
|
|
c write(6,*) 'Beg'
|
|
c do i=1,nmatrix
|
|
c write(6,"(10000f10.6)") (c(i,j),j=1,nbasis)
|
|
c enddo
|
|
C Read variables
|
|
call getvar('natrange ',natrange)
|
|
call getkey('scftol',6,scftol,4)
|
|
read(scftol,*) i
|
|
tol=min(1.d-12,10.d0**(-i))
|
|
C Read S matrix
|
|
c rewind(ifile)
|
|
c call roeint(work,work,s,ifile,nbasis)
|
|
C V = S * C
|
|
call dsymm('l','l',nbasis,nmatrix,1.d0,s,nbasis,c,nbasis,0.d0,v,
|
|
$nbasis)
|
|
C PM iteration
|
|
write(*,"(' Step Convergence')")
|
|
do it=1,300
|
|
ldorot=lsym.and.it.eq.1
|
|
norm=0.d0
|
|
am=0.d0
|
|
bm=0.d0
|
|
u=0.d0
|
|
C Calculate Mulliken charges and parameters A and B
|
|
do iatoms=1,natoms
|
|
call dsyr2k('l','t',nmatrix,natrange(2,iatoms)-
|
|
$natrange(1,iatoms),0.5d0,c(natrange(1,iatoms)+1,1),nbasis,
|
|
$v(natrange(1,iatoms)+1,1),nbasis,0.d0,q,nmatrix)
|
|
do i=2,nmatrix
|
|
do j=1,i-1
|
|
qii=q(i,i)-q(j,j)
|
|
qij=q(i,j)
|
|
am(i,j)=am(i,j)+qij**2-0.25d0*qii**2
|
|
bm(i,j)=bm(i,j)+qij*qii
|
|
enddo
|
|
enddo
|
|
enddo
|
|
do i=2,nmatrix
|
|
do j=1,i-1
|
|
C Calculate rotation angle
|
|
a=am(i,j)
|
|
b=bm(i,j)
|
|
denom=a*a+b*b
|
|
if(denom.gt.0.1d0*tol) then
|
|
angle=sign(0.25d0*dacos(-a/dsqrt(denom)),b)
|
|
if(ldorot.and.dabs(angle).lt.1.d-4) angle=0.1d0
|
|
else
|
|
angle=0.d0
|
|
endif
|
|
C 2x2 rotation
|
|
sina=dsin(angle)
|
|
cosa=dcos(angle)
|
|
u=0.d0
|
|
do k=1,nmatrix
|
|
u(k,k)=1.d0
|
|
enddo
|
|
u(i,i)=cosa
|
|
u(i,j)=-sina
|
|
u(j,i)=sina
|
|
u(j,j)=cosa
|
|
call dgemm('n','n',nbasis,nmatrix,nmatrix,1.d0,c,nbasis,
|
|
$u,nmatrix,0.d0,work,nbasis)
|
|
call dcopy(nbasis*nmatrix,work,1,c,1)
|
|
call dgemm('n','n',nbasis,nmatrix,nmatrix,1.d0,v,nbasis,
|
|
$u,nmatrix,0.d0,work,nbasis)
|
|
call dcopy(nbasis*nmatrix,work,1,v,1)
|
|
norm=norm+1.d0-cosa**2
|
|
enddo
|
|
enddo
|
|
c do j=1,nmatrix
|
|
c write(6,"(i3,1000f10.6)") j,(c(i,j),i=1,nbasis)
|
|
c enddo
|
|
write(*,"(i4,f17.12)") it,norm
|
|
if(norm.lt.tol) exit
|
|
enddo
|
|
C
|
|
write(6,*) 'End'
|
|
do i=1,nbasis
|
|
write(6,"(10000f10.6)") (c(i,j),j=1,nmatrix)
|
|
enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
C***********************************************************************
|
|
C***********************************************************************
|
|
C Subroutine to calculate initial guess by Cholesky decomposing the
|
|
C density matrix of the occupied orbitals
|
|
C***********************************************************************
|
|
C***********************************************************************
|
|
C
|
|
subroutine init_cholesky(
|
|
&nmatrix,
|
|
&nbasis,
|
|
&mocoef_chol,
|
|
&density,
|
|
&workp,
|
|
&work,
|
|
&jpvt)
|
|
C
|
|
implicit none
|
|
C
|
|
C Size of matrices, number of basis functions
|
|
integer nmatrix,nbasis
|
|
C
|
|
C MO coefficients from the Cholesky decomposition to check
|
|
C%%% localization -> majd nem kell!
|
|
real*8 mocoef_chol(nmatrix,*)
|
|
C
|
|
C Density matrix
|
|
real*8 density(nbasis,*)
|
|
C
|
|
C Work array for permutate_chol subroutine
|
|
real*8 workp
|
|
C
|
|
C Work arrays for dchdc subroutine
|
|
real*8 work
|
|
C
|
|
integer jpvt
|
|
integer info
|
|
C
|
|
integer i,j
|
|
C
|
|
C
|
|
C Cholesky decomposition of the density matrix
|
|
call ifillzero(jpvt,nbasis)
|
|
call dchdc(density,nbasis,nbasis,work,jpvt,1,info)
|
|
do i=1,nbasis
|
|
do j=i,nbasis
|
|
density(j,i)=density(i,j)
|
|
if (i.ne.j) density(i,j)=0.d0
|
|
enddo
|
|
enddo
|
|
call permutate_chol(nmatrix,nbasis,jpvt,density,workp)
|
|
do i=1,nmatrix
|
|
do j=1,nbasis
|
|
mocoef_chol(i,j)=density(j,i)
|
|
enddo
|
|
enddo
|
|
C
|
|
end
|
|
C
|
|
C***********************************************************************
|
|
C***********************************************************************
|
|
C Subroutine to permutate Cholesky components
|
|
C***********************************************************************
|
|
C***********************************************************************
|
|
C
|
|
subroutine permutate_chol(
|
|
&nmatrix,
|
|
&nbasis,
|
|
&jpvt,
|
|
&low_triang,
|
|
&work)
|
|
C
|
|
implicit none
|
|
C
|
|
C Size of matrices, number of basis functions
|
|
integer nmatrix,nbasis
|
|
C
|
|
C Array of permutation indexes
|
|
integer jpvt(*)
|
|
C
|
|
C Matrix to permutate (input/output)
|
|
real*8 low_triang(nbasis,*)
|
|
C
|
|
C Work array
|
|
real*8 work(nbasis,*)
|
|
C
|
|
integer i,j
|
|
C
|
|
C
|
|
do i=1,nbasis
|
|
do j=1,nmatrix
|
|
work(jpvt(i),j)=low_triang(i,j)
|
|
enddo
|
|
enddo
|
|
call dcopy(nmatrix*nbasis,work,1,low_triang,1)
|
|
C
|
|
end
|
|
************************************************************************
|
|
subroutine set_density(
|
|
&nmatrix,
|
|
&nbasis,
|
|
&density,
|
|
&mocoef_cmo)
|
|
************************************************************************
|
|
C%%% Subroutine to create density matrix from canonical molecular
|
|
C coefficients -> veglegesen elvielg suruseget kapunk rogton!
|
|
************************************************************************
|
|
implicit none
|
|
C
|
|
C Size of matrices, number of basis functions
|
|
integer nmatrix,nbasis
|
|
C
|
|
C Density matrix
|
|
real*8 density
|
|
C
|
|
C MO coefficients (canonical)
|
|
real*8 mocoef_cmo
|
|
C
|
|
C
|
|
C Create density matrix
|
|
call dgemm('t','n',nbasis,nbasis,nmatrix,1.d0,mocoef_cmo,nmatrix,
|
|
& mocoef_cmo,nmatrix,0.d0,density,nbasis)
|
|
C
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine spade(nbasis,nMOs_total,nMOs_active,MOs,
|
|
$r8heap_nbfXX2,scr,eig,Sroot,natoms,minpfile,natrange,verbosity,
|
|
$ncorenew,iheap1a,iheap1b,lembed,job)
|
|
! $ncorenew,lembed,job)
|
|
************************************************************************
|
|
* 1) Truncate orbitals on the active atoms
|
|
* 2) Build a projector using these orbitals
|
|
* 3) Decompose the projector, and use the eigenvectors for orbital rotation
|
|
************************************************************************
|
|
use error_handler
|
|
implicit none
|
|
integer, intent(in) :: nbasis
|
|
integer, intent(in) :: nMOs_total
|
|
integer, intent(in) :: natoms
|
|
integer, intent(in) :: minpfile
|
|
integer, intent(in) :: natrange(2,natoms)
|
|
integer, intent(in) :: verbosity
|
|
!
|
|
integer, intent(in), target :: iheap1a(natoms)
|
|
integer, intent(in), target :: iheap1b(natoms)
|
|
!
|
|
integer, intent(out) :: ncorenew
|
|
integer, intent(out) :: nMOs_active
|
|
double precision, intent(in) :: Sroot(nbasis,nbasis)
|
|
double precision, intent(inout) :: MOs(nbasis,nMOs_total)
|
|
double precision, intent(in), target :: r8heap_nbfXX2
|
|
$ (nbasis,nMOs_total)
|
|
double precision, intent(inout) :: eig(nMOs_total+1)
|
|
double precision, intent(inout) :: scr(*)
|
|
character(len=1 ), intent(in) :: job
|
|
logical, intent(in) :: lembed
|
|
! pointers
|
|
double precision, pointer :: scr_MOs(:,:) => null()
|
|
integer, pointer :: ncorb(:) => null()
|
|
integer, pointer :: active_atoms(:) => null()
|
|
! local
|
|
character(len=4) :: symm
|
|
integer :: i,j,k
|
|
integer :: MO_specifier
|
|
integer :: istat
|
|
scr_MOs => r8heap_nbfXX2
|
|
ncorb => iheap1a
|
|
active_atoms => iheap1b
|
|
call getvar('ncorb ',ncorb)
|
|
call getkey('symm',4,symm,4)
|
|
if( symm .ne. '0 ' ) call no_implementation_error
|
|
$( " The SPADE localization cannot handle symmetric systems.")
|
|
call read_embedded_subsystem
|
|
call transform_MOs_to_orthogonal_basis
|
|
call truncate_orbitals_to_embedded_subsystem
|
|
call construct_projector
|
|
call decompose_projector
|
|
call determine_subsystem_boundary
|
|
call rotate_MOs
|
|
call print_results
|
|
scr_MOs => null()
|
|
|
|
contains
|
|
|
|
subroutine read_embedded_subsystem
|
|
character(len=32) :: dft
|
|
character(len=16) :: cscr16 !HB
|
|
MO_specifier = 0
|
|
open(minpfile,file='MINP',status='OLD',action='READ',
|
|
$ form='FORMATTED',iostat=istat,position='REWIND')
|
|
if(lembed) then !HB
|
|
call embedat
|
|
$ (natoms,minpfile,active_atoms,6,'embed ',5)
|
|
read(minpfile,*,iostat=istat) dft
|
|
if( istat .ne. 0 ) call io_error
|
|
$("Failed to read the 3rd line of 'embed' during spade"//
|
|
$"localization","pml.f")
|
|
if(trim(dft).eq.'user'.or.trim(dft).eq.'userd') then
|
|
read(minpfile,*,iostat=istat) j
|
|
if( istat .ne. 0 ) call io_error
|
|
$("Failed to read the 4th line of 'embed' during spade"//
|
|
$"localization","pml.f")
|
|
do k=1,j
|
|
read(minpfile,*,iostat=istat)
|
|
if( istat .ne. 0 ) call io_error
|
|
$("Failed to read after the 4th line of 'embed' during spade"//
|
|
$"localization","pml.f")
|
|
enddo
|
|
if(trim(dft).eq.'userd') then
|
|
read(minpfile,*,iostat=istat) j
|
|
do k=1,j
|
|
read(minpfile,*,iostat=istat)
|
|
if( istat .ne. 0 ) call io_error
|
|
$("Failed to read after the 4th line of 'embed' during spade"//
|
|
$"localization","pml.f")
|
|
enddo
|
|
endif
|
|
endif
|
|
read(minpfile,*,iostat=istat) MO_specifier
|
|
if( istat .ne. 0 ) call io_error
|
|
$("Failed to read MO-specifier line of 'embed' during spade"//
|
|
$"localization","pml.f")
|
|
if(job.eq.'v')
|
|
$ read(minpfile,*,iostat=istat) MO_specifier
|
|
else !HB
|
|
if(job.eq.'c') call getkeym('orblocc',7,cscr16,16)
|
|
if(job.eq.'o') call getkeym('orbloco',7,cscr16,16)
|
|
if(job.eq.'v') call getkeym('orblocv',7,cscr16,16)
|
|
read(minpfile,'(a)',iostat=istat) cscr16
|
|
if( istat .ne. 0 ) call io_error
|
|
$("Failed to read MO-specifier line of 'embed' during spade"//
|
|
$"localization","pml.f")
|
|
if( len_trim(cscr16) .eq. 0 )
|
|
$ call illegal_value_error
|
|
$(" The line following the orbloc keyword is empty.",
|
|
$"Empty line.","Atom IDs in the usual format.")
|
|
backspace(minpfile)
|
|
call readlinelist
|
|
$( natoms , minpfile , active_atoms , 6 , 'atoms ')
|
|
endif !HB
|
|
close(minpfile)
|
|
end subroutine read_embedded_subsystem
|
|
!
|
|
subroutine transform_MOs_to_orthogonal_basis
|
|
call dsymm
|
|
$('l' , 'l' , nbasis , nMOs_total , 1.d0 , Sroot , nbasis , MOs ,
|
|
$ nbasis , 0.d0 , scr_MOs , nbasis )
|
|
end subroutine transform_MOs_to_orthogonal_basis
|
|
!
|
|
subroutine truncate_orbitals_to_embedded_subsystem
|
|
integer :: iatom
|
|
ncorenew=0
|
|
do iatom=1,natoms
|
|
if( active_atoms( iatom ) .eq. 1 ) then
|
|
ncorenew = ncorenew + ncorb( iatom )
|
|
else
|
|
scr_MOs( natrange(1,iatom)+1 : natrange(2,iatom) ,
|
|
$ 1 : nMOs_total )
|
|
$ = 0.d0
|
|
endif
|
|
enddo
|
|
end subroutine truncate_orbitals_to_embedded_subsystem
|
|
!
|
|
subroutine construct_projector
|
|
call dsyrk
|
|
$( 'u' , 't' , nMOs_total , nbasis , -1.d0 , scr_MOs ,
|
|
$ nbasis, 0.d0 , scr,nMOs_total )
|
|
end subroutine construct_projector
|
|
!
|
|
subroutine decompose_projector
|
|
call dsyev
|
|
$( 'V' , 'U' , nMOs_total , scr , nMOs_total , eig ,
|
|
$ scr(nMOs_total*nMOs_total+1) , 3*nMOs_total*nMOs_total , i )
|
|
end subroutine decompose_projector
|
|
!
|
|
subroutine determine_subsystem_boundary
|
|
double precision :: tmp
|
|
eig=dsqrt(dabs(eig))
|
|
c eig=-eig
|
|
if(MO_specifier .eq. 0 .and. nMOs_total .ge. 2 ) then
|
|
eig( nMOs_total + 1 ) = 0.d0
|
|
nMOs_active = 1
|
|
tmp = eig(1) - eig(2)
|
|
do i = 2 , nMOs_total
|
|
if( eig(1) .lt. 1d-6 ) then
|
|
nMOs_active = 0
|
|
exit
|
|
endif
|
|
if( eig(i) - eig(i+1) .ge. tmp ) then
|
|
nMOs_active = i
|
|
tmp = eig(i) - eig(i+1)
|
|
endif
|
|
enddo
|
|
else if( MO_specifier .eq. 0 .and. nMOs_total .lt. 2 ) then
|
|
nMOs_active = 0
|
|
if( nMOs_total .eq. 1 .and. eig(1) .gt. 0.5d0 )
|
|
$ nMOs_active = 1
|
|
else if( MO_specifier .gt. 0 .and.
|
|
$ MO_specifier .le. nMOs_total) then
|
|
nMOs_active = MO_specifier
|
|
else
|
|
call out_of_range_error(" The MO number specifier"
|
|
$//" (3rd or 4th line after the keyword 'embed' or 'corembed') "
|
|
$//"is out of range.", nMOs_total , .true., 0 , .true. )
|
|
endif
|
|
end subroutine determine_subsystem_boundary
|
|
!
|
|
subroutine rotate_MOs
|
|
call dgemm
|
|
$( 'n' , 'n' , nbasis , nMOs_total , nMOs_total , 1.d0 ,
|
|
$ MOs,nbasis , scr , nMOs_total , 0.d0 , scr_MOs , nbasis )
|
|
call dcopy(nbasis*nMOs_total , scr_MOs , 1 , MOs , 1 )
|
|
end subroutine rotate_MOs
|
|
!
|
|
subroutine print_results
|
|
if( verbosity .ge. 3 ) then
|
|
write(6,*)
|
|
write(6,'(a)') !HB
|
|
$' Square root of natural orbital occupation numbers:'
|
|
write(6,"(7f10.6)") eig(1:nMOs_total)
|
|
write(6,*)
|
|
write(6,'(a,i6)')
|
|
$' Number of SPADE MOs that are localized in the subsystem: ',
|
|
$ nMOs_active
|
|
write(6,'(a)') ' (based on the change of the '//
|
|
$'occupation numbers)' !HB
|
|
endif
|
|
end subroutine print_results
|
|
end subroutine spade
|
|
C
|
|
************************************************************************
|
|
subroutine fragloc1(nb,dens,scr,eig,si,ss,c,densa,natoms,minpfile,
|
|
$natrange,no,ne,fac,nbnew,verbosity,densb,lenv,ncorenew,sm)
|
|
************************************************************************
|
|
* Create orbitals from density matrix for a fragment
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ne,i,natoms,minpfile,ind(natoms),natrange(2,natoms)
|
|
integer iatoms,jatoms,no,nf,nbnew,verbosity,ncorb(natoms),ncorenew
|
|
real*8 scr(*),dens(nb,nb),si(nb,nb),ss(nb,nb),fac,densb(nb,nb)
|
|
real*8 c(nb,no),eig(nb),densa(nb,nb),ddot,sm(nb,nb)
|
|
logical lenv
|
|
C Read atoms of embedded subsystem
|
|
open(minpfile,file='MINP')
|
|
call embedat(natoms,minpfile,ind,6,'embed ',5)
|
|
read(minpfile,*)
|
|
read(minpfile,*) ne
|
|
close(minpfile)
|
|
call getvar('ncorb ',ncorb)
|
|
C Construct density and orbitals for embedded subsystem
|
|
densa=0.d0
|
|
nbnew=0
|
|
ncorenew=0
|
|
if(.false.) then
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.1) then
|
|
nbnew=nbnew+natrange(2,iatoms)-natrange(1,iatoms)
|
|
ncorenew=ncorenew+ncorb(iatoms)
|
|
endif
|
|
do jatoms=1,natoms
|
|
if(ind(iatoms).eq.1.or.ind(jatoms).eq.1) then
|
|
densa(natrange(1,iatoms)+1:natrange(2,iatoms),
|
|
$ natrange(1,jatoms)+1:natrange(2,jatoms))=
|
|
$ dens (natrange(1,iatoms)+1:natrange(2,iatoms),
|
|
$ natrange(1,jatoms)+1:natrange(2,jatoms))
|
|
endif
|
|
enddo
|
|
enddo
|
|
else
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.1) then
|
|
nbnew=nbnew+natrange(2,iatoms)-natrange(1,iatoms)
|
|
ncorenew=ncorenew+ncorb(iatoms)
|
|
do jatoms=1,natoms
|
|
if(ind(jatoms).eq.1) then
|
|
densa(natrange(1,iatoms)+1:natrange(2,iatoms),
|
|
$ natrange(1,jatoms)+1:natrange(2,jatoms))=
|
|
$ dens (natrange(1,iatoms)+1:natrange(2,iatoms),
|
|
$ natrange(1,jatoms)+1:natrange(2,jatoms))
|
|
endif
|
|
enddo
|
|
endif
|
|
enddo
|
|
endif
|
|
c write(6,*) 'dens'
|
|
c write(6,"(13f10.6)") dens
|
|
c write(6,*) 'densa1'
|
|
c write(6,"(13f10.6)") densa
|
|
C Transform density matrix to orthogonal basis
|
|
call dsymm('l','l',nb,nb,1.d0,densa,nb,ss,nb,0.d0,scr,nb)
|
|
call dsymm('l','l',nb,nb,-1.d0,ss,nb,scr,nb,0.d0,densa,nb)
|
|
c write(6,*) 'densa2'
|
|
c write(6,"(13f10.6)") densa
|
|
C Diagonalize density matrix to get NOs
|
|
call dsyev('V','U',nb,densa,nb,eig,scr,3*nb*nb,i)
|
|
c write(6,*) 'densa3'
|
|
c write(6,"(13f10.6)") densa
|
|
if(ne.eq.0) then
|
|
ne=1
|
|
do while(dabs(eig(ne)).gt.fac*0.75d0.and.ne.le.no)
|
|
ne=ne+1
|
|
enddo
|
|
ne=ne-1
|
|
endif
|
|
if(verbosity.ge.3) then
|
|
write(6,*)
|
|
write(6,*) 'Natural orbital occupation numbers:'
|
|
write(6,"(7f10.6)") -eig
|
|
endif
|
|
c call dsymm('l','u',nb,ne,1.d0,si,nb,densa,nb,0.d0,c,nb)
|
|
C Project NOs onto the occupied space
|
|
call dsymm('l','u',nb,ne,1.d0,ss,nb,densa,nb,0.d0,scr,nb)
|
|
call dsymm('l','l',nb,ne,0.5d0,dens,nb,scr,nb,0.d0,c,nb)
|
|
C Orthogonalize projected NOs
|
|
call dsymm('l','u',nb,ne,1.d0,ss,nb,c,nb,0.d0,scr(nb*nb+1),nb)
|
|
call dsyrk('u','t',ne,nb,1.d0,scr(nb*nb+1),nb,0.d0,scr,ne)
|
|
call filllo(scr,ne)
|
|
call invsqrt(scr,ne,scr(ne*ne+1),6,i,0.d0,1)
|
|
if(i.ne.0) call mrccend(1)
|
|
call dsymm('r','u',nb,ne,1.d0,scr,ne,c,nb,0.d0,scr(ne*ne+1),nb)
|
|
call dcopy(nb*ne,scr(ne*ne+1),1,c,1)
|
|
c call orth(nb,ne,s,c,scr)
|
|
C Calculate active-space density
|
|
call dsyrk('u','n',nb,ne,fac,c,nb,0.d0,densa,nb)
|
|
call filllo(densa,nb)
|
|
if(.not.lenv) return
|
|
C Construct density for environment
|
|
call dcopy(nb*nb,dens,1,densb,1)
|
|
call daxpy(nb*nb,-1.d0,densa,1,densb,1)
|
|
C Construct orbitals for environment
|
|
call dsymm('l','l',nb,nb,1.d0,densb,nb,ss,nb,0.d0,scr,nb)
|
|
call dsymm('l','l',nb,nb,-1.d0,ss,nb,scr,nb,0.d0,densb,nb)
|
|
call dsyev('V','U',nb,densb,nb,eig,scr,3*nb*nb,i)
|
|
call dsymm('l','u',nb,no-ne,1.d0,si,nb,densb,nb,0.d0,c(1,ne+1),nb)
|
|
c write(6,*) 'eig2',-sum(eig(1:nb))
|
|
c write(6,"(13f10.6)") -eig
|
|
c call dsyrk('u','n',nb,no,fac,c,nb,0.d0,dens,nb)
|
|
c call filllo(dens,nb)
|
|
c write(6,*) 'dens'
|
|
c write(6,"(13f10.6)") dens
|
|
c write(6,*) 'c',ddot(nb*nb,densa,1,s,1),ddot(nb*nb,densb,1,s,1)
|
|
c do i=1,nb
|
|
c write(6,"(1000f10.6)") c(i,1:no)
|
|
c enddo
|
|
C Pipek-Mezey localization
|
|
c open(unit=999,file='OEINT',status='old',form='unformatted')
|
|
c call roeint(scr(nb*nb+1),scr(nb*nb+1),scr,999,nb)
|
|
c close(999)
|
|
c call pml(no,nb,natoms,natrange,scr,c,scr(nb*nb+1),scr(2*nb*nb+1),
|
|
c $.false.)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine fragloc(nb,dens,scr,eig,si,ss,ca,cb,natoms,minpfile,
|
|
$natrange,no,ne,fac,nbnew,verbosity,densb,lenv,ncorenew,sm,mb,pm,
|
|
$sd,scrfile1)
|
|
************************************************************************
|
|
* Create orbitals from natural orbitals
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ne,nevold,i,natoms,ind(natoms),natrange(2,natoms,4),mb
|
|
integer iatoms,jatoms,katoms,no,verbosity,ncorb(natoms),ncorenew
|
|
integer iat(natoms),nn,minpfile,nea1,nea2,neb1,neb2,nea,neb,nbnew
|
|
integer scrfile1,nnn
|
|
real*8 scr(*),dens(nb,nb),si(nb,nb),ss(nb,nb),fac,densb(nb,nb),eps
|
|
real*8 ca(nb,no),eig(nb),cb(nb,nb),ddot,sm(nb,nb),pm(nb,mb)
|
|
real*8 sd(mb,mb)
|
|
logical lenv
|
|
parameter(eps=0.05d0)
|
|
C Project density onto the minimal basis, transform to orthogonal basis
|
|
open(scrfile1,file='S12MAT',form='UNFORMATTED')
|
|
call roeint(scr,scr,sm ,scrfile1,mb)
|
|
call rtdmx (scr,scr,densb,scrfile1,nb,mb)
|
|
close(scrfile1)
|
|
call dsymm('r','l',nb,mb,1.d0,sm,mb,densb,nb,0.d0,pm,nb)
|
|
call dsymm('l','l',nb,mb,1.d0,dens,nb,pm,nb,0.d0,scr,nb)
|
|
call dgemm('t','n',mb,mb,nb,-1.d0,pm,nb,scr,nb,0.d0,sd,mb)
|
|
C Read atoms of embedded subsystem
|
|
open(minpfile,file='MINP')
|
|
call embedat(natoms,minpfile,ind,6,'embed ',5)
|
|
close(minpfile)
|
|
call getvar('ncorb ',ncorb)
|
|
C Transform density matrix to orthogonal basis
|
|
call dsymm('l','l',nb,nb,1.d0,dens,nb,ss,nb,0.d0,scr,nb)
|
|
call dsymm('l','l',nb,nb,-1.d0,ss,nb,scr,nb,0.d0,sm,nb)
|
|
C Loop over atoms
|
|
nn=0
|
|
ne=0
|
|
nea=0
|
|
neb=0
|
|
c goto 4567
|
|
nea2=0
|
|
neb2=0
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.1) then
|
|
nbnew=nbnew+natrange(2,iatoms,1)-natrange(1,iatoms,1)
|
|
ncorenew=ncorenew+ncorb(iatoms)
|
|
iat(1)=iatoms
|
|
call noloc1(nb,dens,scr,eig,ss,ca,densb,1,natrange,no,fac,
|
|
$verbosity,ne,iat,sm,sd,mb,natrange(1,1,4),nnn,eps)
|
|
endif
|
|
enddo
|
|
nea1=ne
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.0) then
|
|
iat(1)=iatoms
|
|
call noloc1(nb,dens,scr,eig,ss,ca,densb,1,natrange,no,fac,
|
|
$verbosity,ne,iat,sm,sd,mb,natrange(1,1,4),nnn,eps)
|
|
endif
|
|
enddo
|
|
neb1=ne-nea1
|
|
nea=nea1
|
|
neb=neb1
|
|
call noloc2(nb,scr,ss,ca,densb,fac,dens,ne,sm,mb,pm,sd)
|
|
nevold=ne
|
|
call dcopy(nb*neb1,ca(1,nea1+1),1,cb,1)
|
|
if(ne.eq.no) goto 1234
|
|
C Loop over atom pairs
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.1) then
|
|
do jatoms=1,natoms
|
|
if((ind(jatoms).eq.1.and.jatoms.lt.iatoms).or.
|
|
$ ind(jatoms).eq.0) then
|
|
iat(1)=iatoms
|
|
iat(2)=jatoms
|
|
call noloc1(nb,dens,scr,eig,ss,ca,densb,2,natrange,no,fac,
|
|
$verbosity,ne,iat,sm,sd,mb,natrange(1,1,4),nnn,eps)
|
|
endif
|
|
enddo
|
|
endif
|
|
enddo
|
|
nea2=ne-nevold
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.0) then
|
|
do jatoms=1,iatoms-1
|
|
if(ind(jatoms).eq.0) then
|
|
iat(1)=iatoms
|
|
iat(2)=jatoms
|
|
call noloc1(nb,dens,scr,eig,ss,ca,densb,2,natrange,no,fac,
|
|
$verbosity,ne,iat,sm,sd,mb,natrange(1,1,4),nnn,eps)
|
|
endif
|
|
enddo
|
|
endif
|
|
enddo
|
|
neb2=ne-nevold-nea2
|
|
call noloc2(nb,scr,ss,ca(1,nevold+1),densb,fac,dens,ne-nevold,sm,
|
|
$mb,pm,sd)
|
|
call dcopy(nb*neb2,ca(1,nevold+nea2+1),1,cb(1,neb1+1),1)
|
|
nevold=ne
|
|
nea=nea1+nea2
|
|
neb=neb1+neb2
|
|
do i=1,nea2
|
|
ca(1:nb,nea1+i)=ca(1:nb,nea1+neb1+i)
|
|
enddo
|
|
if(ne.eq.no) goto 1234
|
|
C Construct density and remaining orbitals for the environment
|
|
c4567 continue
|
|
i=0
|
|
do iatoms=1,natoms
|
|
if(ind(iatoms).eq.0) then
|
|
i=i+1
|
|
iat(i)=iatoms
|
|
endif
|
|
enddo
|
|
ne=neb
|
|
call noloc1(nb,dens,scr,eig,ss,cb,densb,i,natrange,no,fac,
|
|
$verbosity,ne,iat,sm,sd,mb,natrange(1,1,4),nnn,eps)
|
|
ne=ne-neb
|
|
call noloc2(nb,scr,ss,cb(1,neb+1),densb,fac,dens,ne,sm,mb,pm,sd)
|
|
neb=neb+ne
|
|
C Construct orbitals for embedded subsystem
|
|
if(nea+neb.lt.no) then
|
|
nn=no-nea-neb
|
|
call dsyev('V','U',nb,sm,nb,eig,scr(nb*nb+1),3*nb*nb,i)
|
|
call dsymm('l','u',nb,nn,1.d0,si,nb,sm,nb,0.d0,ca(1,nea+1),nb)
|
|
endif
|
|
1234 continue
|
|
ne=nea+nn
|
|
call dcopy(nb*neb,cb,1,ca(1,ne+1),1)
|
|
C Reconstruct density matrix
|
|
call dsyrk('u','n',nb,no,fac,ca,nb,0.d0,dens,nb)
|
|
call filllo(dens,nb)
|
|
C Pipek-Mezey localization
|
|
c open(unit=999,file='OEINT',status='old',form='unformatted')
|
|
c call roeint(scr(nb*nb+1),scr(nb*nb+1),scr,999,nb)
|
|
c close(999)
|
|
c call pml(no,nb,natoms,natrange,scr,c,scr(nb*nb+1),scr(2*nb*nb+1),
|
|
c $.false.)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine noloc(nb,dens,scr,eig,si,ss,c,densa,natoms,natrange,no,
|
|
$fac,verbosity,densb,sm,mb,pm,sd,scrfile1,atpair)
|
|
************************************************************************
|
|
* Create orbitals from natural orbitals
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ne,nevold,i,natoms,ind(natoms),natrange(2,natoms,4),mb
|
|
integer iatoms,jatoms,katoms,no,verbosity,nn,scrfile1
|
|
real*8 scr(*),dens(nb,nb),si(nb,nb),ss(nb,nb),fac,densb(nb,nb),eps
|
|
real*8 c(nb,no),eig(nb),densa(nb,nb),ddot,sm(nb,nb),pm(nb,mb)
|
|
real*8 sd(mb,mb)
|
|
logical atpair(natoms,natoms),atl(natoms)
|
|
parameter(eps=0.05d0)
|
|
C Project density onto the minimal basis, transform to orthogonal basis
|
|
open(scrfile1,file='S12MAT',form='UNFORMATTED')
|
|
call roeint(scr,scr,sm ,scrfile1,mb)
|
|
call rtdmx (scr,scr,densa,scrfile1,nb,mb)
|
|
close(scrfile1)
|
|
call dsymm('r','l',nb,mb,1.d0,sm,mb,densa,nb,0.d0,pm,nb)
|
|
call dsymm('l','l',nb,mb,1.d0,dens,nb,pm,nb,0.d0,scr,nb)
|
|
call dgemm('t','n',mb,mb,nb,-1.d0,pm,nb,scr,nb,0.d0,sd,mb)
|
|
C Transform density matrix to orthogonal basis
|
|
call dsymm('l','l',nb,nb,1.d0,dens,nb,ss,nb,0.d0,scr,nb)
|
|
call dsymm('l','l',nb,nb,-1.d0,ss,nb,scr,nb,0.d0,sm,nb)
|
|
C
|
|
call dcopy(nb*nb,dens,1,densb,1)
|
|
ne=0
|
|
nevold=0
|
|
C Loop over atoms
|
|
do iatoms=1,natoms
|
|
ind(1)=iatoms
|
|
call noloc1(nb,densb,scr,eig,ss,c,densa,1,natrange,no,fac,
|
|
$verbosity,ne,ind,sm,sd,mb,natrange(1,1,4),nn,eps)
|
|
enddo
|
|
call noloc2(nb,scr,ss,c(1,nevold+1),densa,fac,densb,ne-nevold,sm,
|
|
$mb,pm,sd)
|
|
nevold=ne
|
|
if(ne.eq.no) goto 1234
|
|
C Loop over atom pairs
|
|
c atl=.true.
|
|
c do iatoms=1,natoms
|
|
c if(atl(iatoms)) then
|
|
c do jatoms=1,iatoms-1
|
|
c if(atl(jatoms)) then
|
|
c ind(1)=iatoms
|
|
c ind(2)=jatoms
|
|
c call noloc1(nb,densb,scr,eig,ss,c,densa,2,natrange,no,fac,
|
|
c $verbosity,ne,ind,sm,sd,mb,natrange(1,1,4),nn,0.25d0)
|
|
c if(nn.eq.1) ne=ne-1
|
|
c if(nn.gt.1) then
|
|
c atl(iatoms)=.false.
|
|
c atl(jatoms)=.false.
|
|
c exit
|
|
c endif
|
|
c endif
|
|
c enddo
|
|
c endif
|
|
c enddo
|
|
c call noloc2(nb,scr,ss,c(1,nevold+1),densa,fac,densb,ne-nevold,sm,
|
|
c $mb,pm,sd)
|
|
c nevold=ne
|
|
c if(ne.eq.no) goto 1234
|
|
C Loop over atom pairs
|
|
do iatoms=1,natoms
|
|
do jatoms=1,iatoms-1
|
|
ind(1)=iatoms
|
|
ind(2)=jatoms
|
|
call noloc1(nb,densb,scr,eig,ss,c,densa,2,natrange,no,fac,
|
|
$verbosity,ne,ind,sm,sd,mb,natrange(1,1,4),nn,eps)
|
|
enddo
|
|
enddo
|
|
call noloc2(nb,scr,ss,c(1,nevold+1),densa,fac,densb,ne-nevold,sm,
|
|
$mb,pm,sd)
|
|
nevold=ne
|
|
if(ne.eq.no) goto 1234
|
|
C Loop over atom triplets
|
|
c atpair=.true.
|
|
c do iatoms=1,natoms
|
|
c do jatoms=1,iatoms-1
|
|
c if(atpair(iatoms,jatoms)) then
|
|
c do katoms=1,jatoms-1
|
|
c if(atpair(jatoms,katoms).and.atpair(iatoms,katoms)) then
|
|
c ind(1)=iatoms
|
|
c ind(2)=jatoms
|
|
c ind(3)=katoms
|
|
c call noloc1(nb,densb,scr,eig,ss,c,densa,3,natrange,no,
|
|
c $fac,verbosity,ne,ind,sm,sd,mb,natrange(1,1,4),nn,eps)
|
|
c if(nn.gt.0) then
|
|
c atpair(iatoms,jatoms)=.false.
|
|
c atpair(jatoms,katoms)=.false.
|
|
c atpair(iatoms,katoms)=.false.
|
|
c exit
|
|
c endif
|
|
c endif
|
|
c enddo
|
|
c endif
|
|
c enddo
|
|
c enddo
|
|
c call noloc2(nb,scr,ss,c(1,nevold+1),densa,fac,densb,ne-nevold,sm,
|
|
c $mb,pm,sd)
|
|
c nevold=ne
|
|
C Construct orbitals for remaining density
|
|
if(no.ne.ne) then
|
|
nn=no-ne
|
|
write(6,*) 'Orbitals not localized:',nn
|
|
call dsymm('l','l',nb,nb,1.d0,densb,nb,ss,nb,0.d0,scr(nb*nb+1),nb)
|
|
call dsymm('l','l',nb,nb,-1.d0,ss,nb,scr(nb*nb+1),nb,0.d0,scr,nb)
|
|
call dsyev('V','U',nb,scr,nb,eig,scr(nb*nb+1),3*nb*nb,i)
|
|
c write(6,*) ne,no,nn
|
|
c write(6,*) 'eig'
|
|
c write(6,*) eig
|
|
call dsymm('l','u',nb,nn,1.d0,si,nb,scr,nb,0.d0,c(1,ne+1),nb)
|
|
ne=no
|
|
endif
|
|
c szemet
|
|
1234 continue
|
|
c call dsyrk('u','n',nb,ne,fac,c,nb,0.d0,densa,nb)
|
|
c call filllo(densa,nb)
|
|
c write(6,*) 'dens'
|
|
c write(6,"(7f10.6)") dens
|
|
c write(6,*) 'densa'
|
|
c write(6,"(7f10.6)") densa
|
|
c write(6,*) 'densb'
|
|
c write(6,"(7f10.6)") densb
|
|
c write(6,*) 'c'
|
|
c write(6,"(7f10.6)") c
|
|
c call dsymm('l','u',nb,ne,1.d0,ss,nb,c,nb,0.d0,scr(nb*nb+1),nb)
|
|
c call dsymm('l','u',nb,ne,1.d0,ss,nb,scr(nb*nb+1),nb,0.d0,scr,nb)
|
|
c call dsymm('l','l',nb,ne,0.5d0,dens,nb,scr,nb,0.d0,c,nb)
|
|
c write(6,*) 'c2'
|
|
c write(6,"(7f10.6)") c
|
|
c call dsymm('l','u',nb,ne,1.d0,ss,nb,c,nb,0.d0,scr(nb*nb+1),nb)
|
|
c call dsyrk('u','t',ne,nb,1.d0,scr(nb*nb+1),nb,0.d0,scr,ne)
|
|
c call filllo(scr,ne)
|
|
c write(6,*) 'scr'
|
|
c write(6,"(5f10.6)") scr(1:ne*ne)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine noloc11(nb,dens,scr,eig,ss,c,densa,natoms,natrange,no,
|
|
$fac,verbosity,ne,ind,sm)
|
|
************************************************************************
|
|
* Create orbitals from natural orbitals
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ne,neold,i,ind(*),natrange(2,*),nn,na,na2
|
|
integer iatoms,jatoms,iiatoms,jjatoms,no,verbosity,natoms
|
|
real*8 scr(*),dens(nb,nb),ss(nb,nb),fac,eps
|
|
real*8 c(nb,*),eig(nb),densa(nb*nb),sm(nb,nb)
|
|
parameter(eps=0.1d0)
|
|
C Calculate sqrt S for embedded subsystem
|
|
call extmat(nb,na,ind,sm,scr,natoms,natrange)
|
|
na2=na*na
|
|
call invsqrt(scr,na,scr(na2+1),6,i,0.d0,2)
|
|
if(i.ne.0) call mrccend(1)
|
|
C Construct density for embedded subsystem
|
|
call extmat(nb,na,ind,dens,densa,natoms,natrange)
|
|
C Transform density matrix to orthogonal basis
|
|
call dsymm('l','l',na,na,1.d0,densa,na,scr,na,0.d0,scr(na2+1),na)
|
|
call dsymm('l','l',na,na,-1.d0,scr,na,scr(na2+1),na,0.d0,densa,na)
|
|
C Diagonalize density matrix to get NOs
|
|
call dsyev('V','U',na,densa,na,eig,scr,3*na*na,i)
|
|
neold=ne
|
|
nn=0
|
|
do while(ne+1.le.no.and.dabs(eig(nn+1)).gt.fac*(1.d0-eps))
|
|
ne=ne+1
|
|
nn=nn+1
|
|
enddo
|
|
if(verbosity.ge.3) then
|
|
write(6,*)
|
|
write(6,"(a7,1000i4)") ' Atoms:',(ind(i),i=1,natoms)
|
|
write(6,"(a24,i4)") ' Number of NOs selected:',nn
|
|
write(6,*) 'Natural orbital occupation numbers:'
|
|
write(6,"(7f10.6)") -eig
|
|
endif
|
|
C Calculate inverse sqrt S for embedded subsystem
|
|
call extmat(nb,na,ind,sm,scr,natoms,natrange)
|
|
call invsqrt(scr,na,scr(na2+1),6,i,0.d0,1)
|
|
if(i.ne.0) call mrccend(1)
|
|
C Transform orbitals to the AO basis
|
|
call dsymm('l','l',na,nn,1.d0,scr,na,densa,na,0.d0,scr(na2+1),na)
|
|
c call dcopy(na*nn,densa,1,scr(na2+1),1)
|
|
call extmo(nb,na,nn,ind,c(1,neold+1),scr(na2+1),natoms,natrange)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine noloc1(nb,dens,scr,eig,ss,c,densa,natoms,natrange,no,
|
|
$fac,verbosity,ne,ind,sm,sd,mb,matrange,nn,eps)
|
|
************************************************************************
|
|
* Create orbitals from natural orbitals
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ne,neold,i,ind(*),natrange(2,natoms),nn,na,na2,mb
|
|
integer iatoms,jatoms,iiatoms,jjatoms,no,verbosity,natoms
|
|
integer matrange(2,natoms)
|
|
real*8 scr(*),dens(nb,nb),ss(nb,nb),fac,eps
|
|
real*8 c(nb,*),eig(nb),densa(nb,nb),sm(nb,nb),sd(mb,mb)
|
|
C Construct density for embedded subsystem in minimal basis
|
|
call extmat(mb,na,ind,sd,densa,natoms,matrange)
|
|
C Diagonalize density matrix to get NOs in minimal basis
|
|
call dsyev('V','U',na,densa,na,eig,scr,3*na*na,i)
|
|
neold=ne
|
|
nn=0
|
|
do while(ne+1.le.no.and.dabs(eig(nn+1)).gt.fac*(1.d0-eps))
|
|
ne=ne+1
|
|
nn=nn+1
|
|
enddo
|
|
if(verbosity.ge.3) then
|
|
write(6,*)
|
|
write(6,"(a7,1000i4)") ' Atoms:',(ind(i),i=1,natoms)
|
|
write(6,"(a24,i4)") ' Number of NOs selected:',nn
|
|
write(6,*) 'Natural orbital occupation numbers (minimal basis):'
|
|
write(6,"(7f10.6)") -eig(1:na)
|
|
endif
|
|
C Construct density for embedded subsystem
|
|
call extmat(nb,na,ind,sm,densa,natoms,natrange)
|
|
C Diagonalize density matrix to get NOs
|
|
call dsyev('V','U',na,densa,na,eig,scr,3*na*na,i)
|
|
c neold=ne
|
|
c nn=0
|
|
c do while(ne+1.le.no.and.dabs(eig(nn+1)).gt.fac*(1.d0-eps))
|
|
c ne=ne+1
|
|
c nn=nn+1
|
|
c enddo
|
|
if(verbosity.ge.3) then
|
|
c write(6,*)
|
|
c write(6,"(a7,1000i4)") ' Atoms:',(ind(i),i=1,natoms)
|
|
c write(6,"(a24,i4)") ' Number of NOs selected:',nn
|
|
write(6,*) 'Natural orbital occupation numbers (AO basis):'
|
|
write(6,"(7f10.6)") -eig(1:na)
|
|
endif
|
|
C Expand MO coefficients
|
|
call extmo(nb,na,nn,ind,c(1,neold+1),densa,natoms,natrange)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine noloc2(nb,scr,ss,c,densa,fac,densb,ne,sm,mb,pm,sd)
|
|
************************************************************************
|
|
* Create orbitals from natural orbitals
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ne,i,mb
|
|
real*8 scr(*),ss(nb,nb),fac,densb(nb,nb),c(nb,*),densa(nb,nb),ddot
|
|
real*8 sm(nb,nb),pm(nb,mb),sd(mb,mb)
|
|
C Return if there are no orbitals
|
|
if(ne.eq.0) return
|
|
C Project NOs onto the occupied space - previous localized orbitals
|
|
call dsymm('l','u',nb,ne,1.d0,ss,nb,c,nb,0.d0,scr,nb)
|
|
call dsymm('l','l',nb,ne,1.d0/fac,densb,nb,scr,nb,0.d0,c,nb)
|
|
C Orthogonalize projected NOs
|
|
call dsymm('l','u',nb,ne,1.d0,ss,nb,c,nb,0.d0,scr(nb*nb+1),nb)
|
|
call dsyrk('u','t',ne,nb,1.d0,scr(nb*nb+1),nb,0.d0,scr,ne)
|
|
call filllo(scr,ne)
|
|
call invsqrt(scr,ne,scr(ne*ne+1),6,i,0.d0,1)
|
|
if(i.ne.0) call mrccend(1)
|
|
call dsymm('r','u',nb,ne,1.d0,scr,ne,c,nb,0.d0,scr(ne*ne+1),nb)
|
|
call dcopy(nb*ne,scr(ne*ne+1),1,c,1)
|
|
C Calculate active-space density
|
|
call dsyrk('u','n',nb,ne,fac,c,nb,0.d0,densa,nb)
|
|
call filllo(densa,nb)
|
|
C Construct density for environment
|
|
call daxpy(nb*nb,-1.d0,densa,1,densb,1)
|
|
c szemet
|
|
c call dsymm('l','u',nb,nb,1.d0,ss,nb,ss,nb,0.d0,scr,nb)
|
|
c write(6,*) ddot(nb*nb,densa,1,scr,1),ddot(nb*nb,densb,1,scr,1)
|
|
C Transform density matrix to orthogonal basis
|
|
call dsymm('l','l',nb,nb,1.d0,densb,nb,ss,nb,0.d0,densa,nb)
|
|
call dsymm('l','l',nb,nb,-1.d0,ss,nb,densa,nb,0.d0,sm,nb)
|
|
C Project density onto the minimal basis, transform to orthogonal basis
|
|
call dsymm('l','l',nb,mb,1.d0,densb,nb,pm,nb,0.d0,scr,nb)
|
|
call dgemm('t','n',mb,mb,nb,-1.d0,pm,nb,scr,nb,0.d0,sd,mb)
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine extmat(nb,ns,ind,mb,ms,natoms,natrange)
|
|
************************************************************************
|
|
* Extract a block of a matrix
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,ns,i,j,ind(*),natrange(2,*),nn
|
|
integer iatoms,jatoms,iiatoms,jjatoms,natoms
|
|
real*8 mb(nb,nb),ms(*)
|
|
C
|
|
nn=0
|
|
ns=0
|
|
do jjatoms=1,natoms
|
|
jatoms=ind(jjatoms)
|
|
do j=natrange(1,jatoms)+1,natrange(2,jatoms)
|
|
ns=ns+1
|
|
do iiatoms=1,natoms
|
|
iatoms=ind(iiatoms)
|
|
do i=natrange(1,iatoms)+1,natrange(2,iatoms)
|
|
nn=nn+1
|
|
ms(nn)=mb(i,j)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine extmo(nb,na,nn,ind,cb,cs,natoms,natrange)
|
|
************************************************************************
|
|
* Extract MO coefficients
|
|
************************************************************************
|
|
implicit none
|
|
integer nb,na,ns,i,j,ind(*),natrange(2,*),nn,ii
|
|
integer iatoms,jatoms,iiatoms,jjatoms,natoms
|
|
real*8 cb(nb,nn),cs(na,nn)
|
|
C
|
|
cb=0.d0
|
|
ii=0
|
|
do iiatoms=1,natoms
|
|
iatoms=ind(iiatoms)
|
|
do i=natrange(1,iatoms)+1,natrange(2,iatoms)
|
|
ii=ii+1
|
|
cb(i,1:nn)=cs(ii,1:nn)
|
|
enddo
|
|
enddo
|
|
C
|
|
return
|
|
end
|
|
C
|
|
************************************************************************
|
|
subroutine pao_loc(iintln,scrfile1,oeintfile,minpfile,mocoeffile,
|
|
$iout,lembed,verbosity,natrange,nbasis,nvirt,nocc,ncore,natoms,
|
|
$nfroz,nvfroz,pao_subsys_tol,
|
|
$cmat,rheap2a,rheap2b,rheap2c,rheap2d,rheap1,iheap1)
|
|
************************************************************************
|
|
* Construct projected atomic orbitals using the active atoms of the
|
|
* embed/orblocv keyword
|
|
************************************************************************
|
|
use error_handler
|
|
implicit none
|
|
logical, intent(in) :: lembed
|
|
integer, intent(in) :: minpfile
|
|
integer, intent(in) :: oeintfile
|
|
integer, intent(in) :: scrfile1
|
|
integer, intent(in) :: mocoeffile
|
|
integer, intent(in) :: iout
|
|
integer, intent(in) :: natoms
|
|
integer, intent(in) :: nbasis
|
|
integer, intent(in) :: ncore
|
|
integer, intent(in) :: nocc
|
|
integer, intent(in) :: nvirt
|
|
integer, intent(in) :: nfroz
|
|
integer, intent(in) :: verbosity
|
|
integer, intent(in) :: natrange(2,natoms)
|
|
integer, intent(in) :: iintln
|
|
double precision, intent(in) :: pao_subsys_tol
|
|
!
|
|
integer, intent(inout) :: nvfroz
|
|
double precision, intent(inout) :: cmat(nbasis,nvirt)
|
|
!
|
|
integer, intent(in), target :: iheap1
|
|
$(nbasis)
|
|
double precision, intent(in), target :: rheap1
|
|
$(nbasis)
|
|
double precision, intent(in), target :: rheap2a
|
|
$(nbasis,nbasis)
|
|
double precision, intent(in), target :: rheap2b
|
|
$(nbasis,nbasis)
|
|
double precision, intent(in), target :: rheap2c
|
|
$(nbasis,nbasis)
|
|
double precision, intent(in), target :: rheap2d
|
|
$(nbasis,nbasis)
|
|
integer :: i,mu,iatom
|
|
integer :: nbasis_embedatoms
|
|
integer :: nvirt_embed
|
|
integer :: npao, npao_nonr
|
|
integer :: istat
|
|
character(len=16) :: cscr16
|
|
character(len=128) :: line
|
|
double precision :: temp
|
|
double precision :: ddot,dsqrt
|
|
integer, dimension(:), pointer :: actatom =>null()
|
|
double precision, dimension(:), pointer :: PAOs_eig =>null()
|
|
double precision, dimension(:,:), pointer :: MO_coeff =>null()
|
|
double precision, dimension(:,:), pointer :: read_scr =>null()
|
|
double precision, dimension(:,:), pointer :: work =>null()
|
|
double precision, dimension(:,:), pointer :: PAOs =>null()
|
|
double precision, dimension(:,:), pointer :: PAOs_nonr =>null()
|
|
double precision, dimension(:,:), pointer :: PAOs_can =>null()
|
|
double precision, dimension(:,:), pointer :: AO_overlap =>null()
|
|
double precision, dimension(:,:), pointer :: PAO_overlap =>null()
|
|
double precision, dimension(:,:), pointer :: SxPAOs =>null()
|
|
double precision, dimension(:,:), pointer :: PAO_to_AO =>null()
|
|
double precision, dimension(:,:), pointer :: AO_fock =>null()
|
|
double precision, dimension(:,:), pointer :: PAO_fock =>null()
|
|
double precision, dimension(:,:), pointer :: nonrPAO_to_canPAO
|
|
$ =>null()
|
|
logical :: debug = .false.
|
|
read_scr => rheap2a
|
|
read_scr = 0.0d0
|
|
read_scr => rheap2b
|
|
read_scr = 0.0d0
|
|
read_scr => rheap2c
|
|
read_scr = 0.0d0
|
|
read_scr => rheap2d
|
|
read_scr = 0.0d0
|
|
read_scr => null()
|
|
write(iout,'(a,e12.4)')
|
|
$' Tolerance of non-redundant subsystem PAOs: ',pao_subsys_tol
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' ****** Building non-redundant Projected Atomic Orbitals ****** '
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Constructing Projected Atomic Orbitals... '
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Reading active atom list... '
|
|
! Read the IDs of the active atoms under the embed/orblocv keywords
|
|
open(unit=minpfile,file='MINP',form='FORMATTED',status='OLD',
|
|
$ action='READ',iostat=istat)
|
|
if(istat.ne.0) call io_error
|
|
$(" Failed to open the MINP file" ,"pml.f")
|
|
actatom => iheap1
|
|
actatom=0
|
|
if( lembed ) then
|
|
call getkeym('embed',5,cscr16,8)
|
|
call readlinelist(natoms,minpfile,actatom,iout,'atoms ')
|
|
nbasis_embedatoms = 0
|
|
do iatom=1,natoms
|
|
mu = natrange(1,iatom)+1
|
|
do while( actatom(iatom) .eq. 1 .and.
|
|
$ mu .le. natrange(2,iatom) )
|
|
nbasis_embedatoms = nbasis_embedatoms + 1
|
|
mu = mu + 1
|
|
enddo
|
|
enddo
|
|
endif
|
|
cscr16=''
|
|
call getkeym('orblocv',7,cscr16,16)
|
|
read(minpfile,'(a)') line
|
|
if( len_trim(line) .eq. 0 .and. ( .not. lembed ) )then
|
|
write(iout,"(a)") " 'embed=off' and the list of atoms is"//
|
|
$" not specified under the keyword 'orblocv'"
|
|
call mrccend(1)
|
|
endif
|
|
if( len_trim(line) .eq. 0 .and. lembed ) then
|
|
call getkeym('embed',5,cscr16,8)
|
|
else
|
|
backspace(minpfile)
|
|
endif
|
|
if( verbosity.ge.3 ) then
|
|
read(minpfile,'(a)') line
|
|
write(iout,'(a)')
|
|
$' The following atoms will be used for PAO construction:'
|
|
write(iout,'(a,a)') ' ',trim(adjustl(line))
|
|
backspace(minpfile)
|
|
endif
|
|
call readlinelist(natoms,minpfile,actatom,iout,'atoms ')
|
|
close(minpfile)
|
|
read_scr => rheap2a
|
|
AO_overlap => rheap2b
|
|
MO_coeff => rheap2c
|
|
work => rheap2d
|
|
rewind(mocoeffile) ! mocoeffile is already opened in orbloc
|
|
call readmo
|
|
$( read_scr , read_scr , MO_coeff , mocoeffile ,nbasis , nbasis )
|
|
! Read the AO overlap matrix (the oeintfile is already open in the start of orbloc.f)
|
|
call roeint(read_scr , read_scr , AO_overlap , oeintfile , nbasis)
|
|
!---------------------- MO_coeff : C_i ----------!
|
|
!---------------------- work : C_i C_i^T == R ----------!
|
|
call dsyrk
|
|
$('u','n',nbasis,ncore+nocc,1.d0,MO_coeff,nbasis,0.d0,work,nbasis)
|
|
call filllo(work,nbasis)
|
|
read_scr => null()
|
|
PAOs => rheap2a
|
|
!AO_overlap => rheap2b
|
|
!MO_coeff => rheap2c
|
|
!work => rheap2d
|
|
!--------------------- work : R -------------------!
|
|
!--------------------- AO_overlap : S -------------------!
|
|
!--------------------- PAOs : - R S -------------------!
|
|
call dsymm('l','l', nbasis , nbasis , -1.d0 , work ,
|
|
$ nbasis , AO_overlap , nbasis , 0.d0 , PAOs , nbasis )
|
|
! 1 - RS
|
|
do mu=1,nbasis
|
|
PAOs(mu,mu)=PAOs(mu,mu)+1.d0
|
|
enddo
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Renormalizing and building the overlap matrix (1/2)... '
|
|
!PAOs => rheap2a
|
|
!AO_overlap => rheap2b
|
|
MO_coeff => null()
|
|
SxPAOs => rheap2c
|
|
!work => rheap2d
|
|
!---------------------- AO_overlap : S --------------!
|
|
!---------------------- PAOs : 1 - RS --------------!
|
|
!---------------------- SxPAOs : S ( 1 - RS ) --------------!
|
|
call dsymm('l','l',nbasis,nbasis,1.d0,AO_overlap,nbasis,PAOs,
|
|
$nbasis, 0.d0, SxPAOs , nbasis)
|
|
npao=0
|
|
if( debug ) write(iout,'(a,6x,a)')' PAO ID','[PAO^T x S x PAO]_ii'
|
|
do iatom=1,natoms
|
|
mu = natrange(1,iatom)+1
|
|
do while( actatom(iatom) .eq. 1 .and. mu .le. natrange(2,iatom))
|
|
!do while( mu .le. natrange(2,iatom))
|
|
! ( 1 - RS )^T S ( 1 - RS ) == N
|
|
temp=ddot(nbasis,PAOs(1,mu),1,SxPAOs(1,mu),1)
|
|
if(debug) write(iout,"(i4,12x,e12.4)") mu,temp
|
|
if( dabs(temp) .gt. pao_subsys_tol ) then ! kept PAO list: PAOs with sufficent norm after projection
|
|
npao=npao+1
|
|
temp=1.d0/dsqrt(temp)
|
|
! 1/N^(1/2) ( 1 - RS )
|
|
call dscal(nbasis,temp,PAOs(1,mu),1)
|
|
! 1/N^(1/2) S ( 1 - RS )
|
|
call dscal(nbasis,temp,SxPAOs(1,mu),1)
|
|
if( mu .ne. npao ) then
|
|
call dcopy(nbasis,PAOs(1,mu),1,PAOs(1,npao),1)
|
|
call dcopy(nbasis,SxPAOs(1,mu),1,SxPAOs(1,npao),1)
|
|
endif
|
|
endif
|
|
mu = mu + 1
|
|
enddo
|
|
enddo
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Renormalizing and building the overlap matrix (2/2)... '
|
|
actatom => null()
|
|
PAOs_eig => rheap1
|
|
!PAOs => rheap2a
|
|
AO_overlap => null()
|
|
PAO_overlap => rheap2b
|
|
PAO_to_AO => rheap2b
|
|
!SxPAOs => rheap2c
|
|
!work => rheap2d
|
|
! -------------------- PAOs : ( 1 - RS ) ( nbasis x npao )
|
|
! -------------------- SxPAO : S ( 1 - RS ) ( nbasis x npao )
|
|
! -------------------- PAO_overlap: ( 1 - RS )^T S ( 1 - RS ) == Q (npao x npao)
|
|
call dgemm('t','n', npao , npao , nbasis , 1.d0 , PAOs , nbasis ,
|
|
$ SxPAOs , nbasis , 0.d0 , PAO_overlap , nbasis)
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Diagonalizing the PAO overlap matrix... '
|
|
! Q C = E C
|
|
call dsyev('V','U', npao , PAO_overlap , nbasis , PAOs_eig ,
|
|
$ work , nbasis*nbasis , istat )
|
|
if( istat.ne.0 ) then
|
|
write(iout,'(a)' )
|
|
$'Error during PAO overlap diagonalization ( pao_loc @ pml.f )'
|
|
call mrccend(1)
|
|
endif
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Renormalizing and keeping non-redundant PAOs... '
|
|
if( debug ) write(iout,"(a)")
|
|
$" [PAO overlap eigenvalues]"
|
|
npao_nonr=0
|
|
do i=1,npao
|
|
if(debug) write(iout,"(i4,4x,e12.4)") i,PAOs_eig(i)
|
|
if( abs( PAOs_eig(i) ) .gt. pao_subsys_tol ) then
|
|
npao_nonr = npao_nonr + 1
|
|
if( i .ne. npao_nonr ) then
|
|
call dcopy
|
|
$ (npao,PAO_to_AO(1,i),1,PAO_to_AO(1,npao_nonr),1)
|
|
call dscal
|
|
$(npao,1.d0/dsqrt(PAOs_eig(i)),PAO_to_AO(1,npao_nonr),1)
|
|
endif
|
|
endif
|
|
enddo
|
|
!PAOs => rheap2a
|
|
PAO_overlap => null()
|
|
!PAO_to_AO => rheap2b
|
|
SxPAOs => null()
|
|
PAOs_nonr => rheap2c
|
|
!work => rheap2d
|
|
call dgemm( 'n' , 'n', nbasis , npao_nonr , npao , 1.d0,
|
|
$ PAOs , nbasis , PAO_to_AO , nbasis , 0.d0 , PAOs_nonr , nbasis )
|
|
c do i=1,npao_nonr
|
|
c call dcopy(nbasis,PAOs_nonr(1,i),1,cmat(1,i),1)
|
|
c enddo
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Canonization of non-redundant PAOs... '
|
|
open(unit=scrfile1,file='FOCK',form='UNFORMATTED',status='OLD',
|
|
$ action='READ',iostat=istat)
|
|
if(istat.ne.0) call io_error
|
|
$(" Failed to open the FOCK file","pml.f")
|
|
PAOs => null()
|
|
read_scr => rheap2a
|
|
PAO_to_AO => null()
|
|
AO_fock => rheap2b
|
|
!PAOs_nonr => rheap2c
|
|
!work => rheap2d
|
|
call roeint(read_scr,read_scr,AO_fock,scrfile1,nbasis)
|
|
close(scrfile1)
|
|
read_scr => null()
|
|
PAO_fock => rheap2a
|
|
nonrPAO_to_canPAO => rheap2a
|
|
!AO_fock => rheap2b
|
|
!PAOs_nonr => rheap2c
|
|
!work => rheap2d
|
|
!---------------------- AO_fock : F -------------!
|
|
!---------------------- PAOs_nonr : ( 1 - RS)' -------------!
|
|
!---------------------- work : F ( 1 - RS )' -------------!
|
|
call dsymm('l','l', nbasis , npao_nonr , 1.d0 , AO_fock , nbasis ,
|
|
$ PAOs_nonr, nbasis , 0.d0 , work , nbasis)
|
|
!---------------------- PAOs_nonr : ( 1 - RS )' ------!
|
|
!---------------------- work : F ( 1 - RS )' ------!
|
|
!---------------------- PAO_fock : ( 1 - RS)' F ( 1 - RS )' == Z ------!
|
|
call dgemm('t','n', npao_nonr , npao_nonr , nbasis , 1.d0 ,
|
|
$ PAOs_nonr , nbasis , work , nbasis , 0.d0 , PAO_fock , nbasis )
|
|
! Z C = E C
|
|
call dsyev('V','U', npao_nonr , PAO_fock , nbasis , PAOs_eig ,
|
|
$ work , nbasis*nbasis ,istat )
|
|
if( istat.ne.0 ) then
|
|
write(iout,'(a)' )
|
|
$'Error during FOCK diagonalization ( pao_loc @ pml.f )'
|
|
call mrccend(1)
|
|
endif
|
|
PAO_fock => null()
|
|
!nonrPAO_to_canPAO => rheap2a
|
|
AO_fock => null()
|
|
PAOs_can => rheap2b
|
|
!PAOs_nonr => rheap2c
|
|
!work => rheap2d
|
|
call dgemm('n','n', nbasis , npao_nonr , npao_nonr , 1.d0 ,
|
|
$ PAOs_nonr , nbasis , nonrPAO_to_canPAO , nbasis , 0.d0 ,
|
|
$ PAOs_can , nbasis )
|
|
nvfroz = nvirt - npao_nonr
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' Updating molecular orbitals... '
|
|
do i=1,npao_nonr
|
|
call dcopy(nbasis,PAOs_can(1,i),1,cmat(1,nvfroz+i),1)
|
|
enddo
|
|
call ishell('cp VARS VARS.old')
|
|
open(unit=scrfile1,file='VARS.new',form='UNFORMATTED',
|
|
$ position='ASIS',status='NEW',action='WRITE',iostat=istat)
|
|
if(istat.eq.0) write(scrfile1) 'nvfroz ',iintln,nvfroz
|
|
close(scrfile1)
|
|
call ishell('cat VARS.new VARS.old > VARS')
|
|
if(istat.eq.0) open(unit=scrfile1,file='VARS.new',iostat=istat)
|
|
if(istat.eq.0) close(scrfile1,status='delete',iostat=istat)
|
|
if(istat.eq.0) open(unit=scrfile1,file='VARS.old',iostat=istat)
|
|
if(istat.eq.0) close(scrfile1,status='delete',iostat=istat)
|
|
if( istat .ne. 0 ) call io_error
|
|
$(' Error during the modification of the VARS file',
|
|
$ 'pml.f ')
|
|
do i=1,nvfroz
|
|
cmat(1:nbasis,i) = 0.0d0
|
|
enddo
|
|
if( verbosity .gt. 3 ) write(iout,'(a)')
|
|
$' ************************************************************** '
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
$' The construction of non-redundant projected atomic orbitals
|
|
$ is finished.'
|
|
if( verbosity .ge. 3 ) write(iout,'(a)')
|
|
call print_pao_stats( nbasis , ncore , nocc , nvirt , npao ,
|
|
$ npao_nonr , lembed , iout , nfroz , nbasis_embedatoms )
|
|
PAOs_eig => null()
|
|
nonrPAO_to_canPAO => null()
|
|
PAOs_can => null()
|
|
PAOs_nonr => null()
|
|
work => null()
|
|
end subroutine
|
|
|
|
subroutine edit_molden
|
|
$(read_channel,write_channel,nbasis,ncore,nocc,nvfroz)
|
|
use error_handler
|
|
implicit none
|
|
integer, intent(in) :: read_channel
|
|
integer, intent(in) :: write_channel
|
|
integer, intent(in) :: nbasis
|
|
integer, intent(in) :: ncore
|
|
integer, intent(in) :: nocc
|
|
integer, intent(in) :: nvfroz
|
|
integer :: istat
|
|
integer :: imo,mu
|
|
character(len=128) :: read_buffer
|
|
read_buffer=''
|
|
open(unit=read_channel,file='MOLDEN',status='OLD',action='READ',
|
|
$iostat=istat)
|
|
if( istat .ne. 0 ) call io_error
|
|
$(" Cannot open MOLDEN file","pml.f")
|
|
open(unit=write_channel,file='MOLDEN_new',status='NEW',
|
|
$action='WRITE', iostat=istat)
|
|
if( istat .ne. 0 ) call io_error
|
|
$(" Cannot open MOLDEN_new file","pml.f")
|
|
do while (read_buffer(1:4).ne.'[MO]')
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
enddo
|
|
do imo=1,ncore+nocc
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
do mu=1,nbasis
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
enddo
|
|
enddo
|
|
do imo=ncore+nocc+1,ncore+nocc+nvfroz
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,"(' Ene=',f14.4)") 0.0d0
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
do mu=1,nbasis
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
enddo
|
|
enddo
|
|
do imo=ncore+nocc+nvfroz+1,nbasis
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
do mu=1,nbasis
|
|
read(read_channel,'(a)') read_buffer
|
|
write(write_channel,'(a)') read_buffer
|
|
enddo
|
|
|
|
enddo
|
|
close(read_channel)
|
|
close(write_channel)
|
|
call ishell('mv MOLDEN_new MOLDEN')
|
|
end subroutine edit_molden
|
|
|
|
subroutine print_pao_stats( nbasis , ncore , nocc , nvirt , npao ,
|
|
$ npao_nonr , lembed , iout , nfroz , nbasis_embedatoms )
|
|
implicit none
|
|
integer, intent(in) :: nbasis
|
|
integer, intent(in) :: ncore
|
|
integer, intent(in) :: nocc
|
|
integer, intent(in) :: nvirt
|
|
integer, intent(in) :: npao
|
|
integer, intent(in) :: npao_nonr
|
|
integer, intent(in) :: iout
|
|
integer, intent(in) :: nfroz
|
|
integer, intent(in) :: nbasis_embedatoms
|
|
integer :: nvirt_embed
|
|
logical, intent(in) :: lembed
|
|
write(iout,'(a)')
|
|
$' ----------------- Full-molecule-related statistics ----------- '
|
|
write(iout,'(a,i9)') ' # of AOs :',
|
|
$nbasis
|
|
write(iout,'(a,i9)') ' # of occupied orbitals :',
|
|
$ncore+nocc
|
|
write(iout,'(a,i9)') ' # of virtual orbitals :',
|
|
$nvirt
|
|
write(iout,'(a,i9)') ' # of constructed PAOs :',
|
|
$npao
|
|
write(iout,'(a,i9)') ' # of non-redundant PAOs :',
|
|
$npao_nonr
|
|
write(iout," (' Compression of virtual MOs : '
|
|
$, f4.1,' %')")
|
|
$ ABS( DBLE( npao_nonr - nvirt ) ) / DBLE( nvirt ) * 100d0
|
|
write(iout,'(a)')
|
|
if( lembed .and. nfroz .eq. 0 ) then
|
|
write(iout,'(a)') ' Warning!
|
|
$ An embedding calculation is requested without frozen orbitals.'
|
|
write(iout,'(a)')
|
|
$' Subspace-related statistics will not be printed.'
|
|
endif
|
|
if( lembed .and. nfroz .ne. 0 ) then
|
|
write(iout,'(a)')
|
|
$' ----------------- Subspace-related statistics ---------------- '
|
|
write(iout,'(a,i9)') ' # of AOs on selected atoms :',
|
|
$nbasis_embedatoms
|
|
write(iout,'(a,i9)') ' # of occupied orbitals :',
|
|
$ncore + nocc - nfroz
|
|
nvirt_embed = nbasis_embedatoms - ( ncore + nocc - nfroz )
|
|
write(iout,'(a,i9)') ' Approx. # of virtual orbitals :',
|
|
$ nvirt_embed
|
|
write(iout," (' Virtual subspace expansion : '
|
|
$, f5.1,' %')") ABS( DBLE( npao_nonr - nvirt_embed ) )
|
|
$/ DBLE( nvirt_embed ) * 100.0d0
|
|
endif
|
|
end subroutine print_pao_stats
|
|
|
|
|