lukasplevac/easyconfigs-it4i
1
0
Fork 0
mirror of https://code.it4i.cz/sccs/easyconfigs-it4i.git synced 2025-04-17 12:10:50 +01:00
Code Issues Packages Projects Releases Wiki Activity
easyconfigs-it4i/m/MRCC/mrcc_files/boys.f
Lukas Krupcik eb1ea94130 merge master
2024-07-25 10:27:17 +02:00

1448 lines
44 KiB
Fortran
Raw Blame History

************************************************************************
subroutine init_boys(maxit,nmatrix,xdip,ydip,zdip,usum,work)
************************************************************************
* Boys localization
************************************************************************
implicit none
integer maxit,nmatrix,i,j,k,it
real*8 norm,tol,xii,xjj,xij,yii,yjj,yij,zii,zjj,zij,sina,cosa
real*8 xdip(nmatrix,nmatrix),ydip(nmatrix,nmatrix),work(*),a,b
real*8 zdip(nmatrix,nmatrix),usum(nmatrix,nmatrix),angle,z
logical lsym
C
call boys_tol(lsym,tol)
C usum matrix is unity
call dfillzero(usum,nmatrix**2)
do i=1,nmatrix
usum(i,i)=1.d0
enddo
C Boys iteration
norm=0.d0
angle=0.d0 ! NP
if(maxit.ne.0) write(*,"(' Step Convergence')")
do it=1,maxit
norm=0.d0
c C$OMP PARALLEL DO Schedule(Dynamic)
c C$OMP& DEFAULT(PRIVATE)
c C$OMP& SHARED(nmatrix,tol,usum,xdip,ydip,zdip,work)
c C$OMP& REDUCTION(+:norm)
do j=2,nmatrix
do i=1,j-1
xjj=xdip(j,j)
yjj=ydip(j,j)
zjj=zdip(j,j)
xii=xdip(i,i)
yii=ydip(i,i)
zii=zdip(i,i)
xij=xdip(i,j)
yij=ydip(i,j)
zij=zdip(i,j)
C Calculate rotation angle
a=xij*xij+yij*yij+zij*zij-0.25d0*
$ (xii*xii+yii*yii+zii*zii+xjj*xjj+yjj*yjj+zjj*zjj-2.d0*
$ (xii*xjj+yii*yjj+zii*zjj))
b=xij*(xii-xjj)+yij*(yii-yjj)+zij*(zii-zjj)
if(a.ne.0.d0) angle=0.25d0*datan(b/a)
z=0.d0
if(a.ge.0.d0) z=dasin(1.d0)/2.d0
if(b.gt.0.d0) z=-z
angle=angle+z
if(dabs(4.d0*a)+dabs(2.d0*b).lt.0.1d0*tol) angle=0.d0
if (lsym.and.it.eq.1.and.dabs(angle).lt.1.d-4) angle=1.d-4 ! random mixing in case of symmetry
sina=dsin(angle)
cosa=dcos(angle)
norm=norm+1.d0-cosa**2
C Calculate cumulative the rotation matrix
c C$OMP CRITICAL (us)
work(1:nmatrix) =cosa*usum(i,1:nmatrix)-
$ sina*usum(j,1:nmatrix)
work(nmatrix+1:2*nmatrix)=sina*usum(i,1:nmatrix)+
$ cosa*usum(j,1:nmatrix)
usum(i,1:nmatrix)=work(1:nmatrix)
usum(j,1:nmatrix)=work(nmatrix+1:2*nmatrix)
C Transform dipole integrals, x
work(1:nmatrix) =cosa*xdip(i,1:nmatrix)-
$ sina*xdip(j,1:nmatrix)
work(nmatrix+1:2*nmatrix)=sina*xdip(i,1:nmatrix)+
$ cosa*xdip(j,1:nmatrix)
xdip(i,1:nmatrix)=work(1:nmatrix)
xdip(j,1:nmatrix)=work(nmatrix+1:2*nmatrix)
work(1:nmatrix) =cosa*xdip(1:nmatrix,i)-
$ sina*xdip(1:nmatrix,j)
work(nmatrix+1:2*nmatrix)=sina*xdip(1:nmatrix,i)+
$ cosa*xdip(1:nmatrix,j)
xdip(1:nmatrix,i)=work(1:nmatrix)
xdip(1:nmatrix,j)=work(nmatrix+1:2*nmatrix)
C Transform dipole integrals, y
work(1:nmatrix) =cosa*ydip(i,1:nmatrix)-
$ sina*ydip(j,1:nmatrix)
work(nmatrix+1:2*nmatrix)=sina*ydip(i,1:nmatrix)+
$ cosa*ydip(j,1:nmatrix)
ydip(i,1:nmatrix)=work(1:nmatrix)
ydip(j,1:nmatrix)=work(nmatrix+1:2*nmatrix)
work(1:nmatrix) =cosa*ydip(1:nmatrix,i)-
$ sina*ydip(1:nmatrix,j)
work(nmatrix+1:2*nmatrix)=sina*ydip(1:nmatrix,i)+
$ cosa*ydip(1:nmatrix,j)
ydip(1:nmatrix,i)=work(1:nmatrix)
ydip(1:nmatrix,j)=work(nmatrix+1:2*nmatrix)
C Transform dipole integrals, z
work(1:nmatrix) =cosa*zdip(i,1:nmatrix)-
$ sina*zdip(j,1:nmatrix)
work(nmatrix+1:2*nmatrix)=sina*zdip(i,1:nmatrix)+
$ cosa*zdip(j,1:nmatrix)
zdip(i,1:nmatrix)=work(1:nmatrix)
zdip(j,1:nmatrix)=work(nmatrix+1:2*nmatrix)
work(1:nmatrix) =cosa*zdip(1:nmatrix,i)-
$ sina*zdip(1:nmatrix,j)
work(nmatrix+1:2*nmatrix)=sina*zdip(1:nmatrix,i)+
$ cosa*zdip(1:nmatrix,j)
zdip(1:nmatrix,i)=work(1:nmatrix)
zdip(1:nmatrix,j)=work(nmatrix+1:2*nmatrix)
c C$OMP END CRITICAL (us)
enddo
enddo
c C$OMP END PARALLEL DO
write(*,"(i4,f17.12)") it,norm
if(norm.lt.tol) exit
enddo
C
return
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to transform a matrix (a dipole moment or a second
C moment integral) with umatrix=exp(kappa)
C***********************************************************************
C***********************************************************************
C
subroutine transform_integral(
&nmatrix,
&umatrix,
&x,
&work)
C
implicit none
C
C Size of matrices:
integer nmatrix
C
C Full umatrix=exp(kappa) matrix to transform a matrix
real*8 umatrix
C
C The full matrix to be transformate (input/output)
real*8 x
C
C Work array for dgemm subroutine
real*8 work
C
C
call dgemm('n','n',nmatrix,nmatrix,nmatrix,1.d0,umatrix,nmatrix,
& x,nmatrix,0.d0,work,nmatrix)
call dgemm('n','t',nmatrix,nmatrix,nmatrix,1.d0,work,nmatrix,
& umatrix,nmatrix,0.d0,x,nmatrix)
C
end
C
C***********************************************************************
C***********************************************************************
C Solve newton equation by (micro)iteration and transform integrals
C by (macro)iteration:
C xi2*kappa=-xi1
C The equation is solved by Jacobi method, where the preconditioning
C matrix contains the diagonal elements of the four indexed tensor
C xi2 (Hessian matrix). xi1 (gradientvmatrix) and kappa are two
C indexed antisymmetrix matrices. Their upper triangle part is
C ordered in (nmatrix-1)*nmatrix/2 long vector rowwise.
C
C Journal of Chemical Physics 134, 194104 (2011), B. Jansík et al.
C***********************************************************************
C***********************************************************************
C
subroutine transform(
&m,
&nmatrix,
&tresh_macro,
&tresh_micro,
&itol,
&xdip,
&ydip,
&zdip,
&uboys,
&xi1,
&xi2,
&xi2prec,
&kappa,
&kappaold,
&q,
&omega1,
&xim,
&kappam,
&kq,
&kx,
&ky,
&kz,
&basis_vec,
&tr_basis_vec,
&reduced_hessian,
&max_macro_it,
&max_micro_it,
&work,line_search_mem,lwork)
implicit none
C
C Power of the orbital variances
real*8 m
C
C Size of matrices
integer nmatrix
C
C Tresholds
real*8 tresh_macro,tresh_micro,itol
C
C Dipole moment and second moment matrices
real*8 xdip(nmatrix*nmatrix)
real*8 ydip(nmatrix*nmatrix)
real*8 zdip(nmatrix*nmatrix)
C
C Cumulated full usum matrix, in which the product of the
C umatrix=exp(kappa) transformation matrices is stored
real*8 uboys(nmatrix,nmatrix)
C
C xi1 matrix (right hand of the Newton equation), and xi2 matrix (left
C side of the Newton equation. They are antisymmetric matrices. Their
C upper triangle part is ordered in (nmatrix-1)*nmatrix/2 long vectors
C rowwise
real*8 xi1((nmatrix-1)*nmatrix/2)
real*8 xi2((nmatrix-1)*nmatrix/2)
C
C xi2prec preconditioning matrix (diagonal of the four indexed Hessian
C tensor). Its diagonal is ordered in an (nmatrix-1)*nmatrix/2 long
C vector rowwise
real*8 xi2prec((nmatrix-1)*nmatrix/2 )
C
C Kappa matrix (antisymmetric). Its upper triangle part is ordered in an
C (nmatrix-1)*nmatrix/2 long vector rowwise
real*8 kappa((nmatrix-1)*nmatrix/2)
C
C A temporary kappa matrix to calculate the norm of the difference of
C the kappas before and after the macroiteration
real*8 kappaold((nmatrix-1)*nmatrix/2)
C
C Work arrays for calculate_xi1 and calculate_xi2 subroutines
real*8 q(*),omega1(*),xim(*),kq(*),kx(*),ky(*),kz(*)
real*8 work(*),kappam(*)
C
real*8 dnrm2
real*8 norm,norm1
c Variables for trust-region
integer max_macro_it,max_micro_it
integer nvec,lwork
real*8 trust_radius
real*8 basis_vec((nmatrix-1)*nmatrix/2,max_micro_it)
real*8 tr_basis_vec((nmatrix-1)*nmatrix/2,max_micro_it)
real*8 reduced_hessian(max_micro_it,max_micro_it)
real*8 line_search_mem(lwork)
real*8 alpha,mu,norm1_init,grad_norm,res_init
real*8 reduction_global,reduction_local,reduction
c Variables for tolarance
logical lsym
real*8 tol
C
integer max_it1
integer max_it2
C
integer i,ii,jj,nn
C
integer info,ic,ihc,ia,isa,ieigval,ieigvec,iscr,dim
integer alg
integer itmax,dblalloc,idamax
itmax=50
dim=53
C
nn=(nmatrix-1)*nmatrix/2
C calculating tolerance
call boys_tol(lsym,tol)
C
C Initialize variables
ii=1
max_it1=max_macro_it
max_it2=max_micro_it
reduction_global=0.05d0
reduction_local=0.005d0
reduction=reduction_global
C Calculating xi1 gradient vector
call calculate_xi1(m,nmatrix,xdip,ydip,zdip,xi1,xim,q,omega1)
i=idamax(nn,xi1,1)
norm=dabs(xi1(i))
c write(*,"(' Step',' Convergence ')")
write(*,"(' Step Gradient')")
write(*,"(i4,3X,ES12.5E2)") ii,norm
alg=2
C
C Macro iteration
do while (tol.le.norm .and. ii.le.max_it1)
nvec=0
trust_radius=0.4d0
grad_norm=dnrm2(nn,xi1,1)
C
C xi2prec preconditioning matrix
call calculate_prec(m,nmatrix,xdip,ydip,zdip,xi2prec,q,omega1)
C Modify preconditioning matrix if it is singular
call modify_prec(nmatrix,xi2prec)
C Initialize kappa matrix (zero matrix) and other variables
call dfillzero(kappa,nn)
call dfillzero(xi2,nn) ! H*0=0
jj=1
C
i=idamax(nn,xi1,1)
norm1=dabs(xi1(i))
norm1_init=norm1
res_init=norm1
c if(res_init.lt.1.0d-3) reduction=reduction_local
C Microiteration
do while (reduction*res_init.le.norm1 .and. nvec.le.max_it2)
if(alg.eq.1) then ! simple Newton iteration
c call jacobi(nmatrix,kappa,xi1,xi2,xi2prec)
call gauss_seidel(nmatrix,kappa,xi1,xi2,xi2prec)
C H*kappa
call calculate_xi2(m,nmatrix,xdip,ydip,zdip,
& xi1,xi2,kappa,xim,
& q,omega1,kq,kx,ky,kz)
C Calculate residual
call dcopy(nn,xi2,1,kappaold,1)
call daxpy(nn,1.0d0,xi1,1,kappaold,1)
elseif(alg.eq.2) then !trust-region iteration
C Generating trial vectors
if(nvec.eq.0) then
basis_vec(1:nn,1)=xi1(1:nn)/grad_norm
call calculate_xi2(m,nmatrix,xdip,ydip,zdip,
& xi1,tr_basis_vec(1,1),
& basis_vec(1,1),xim,q,
& omega1,kq,kx,ky,kz)
nvec=1
if(nn.ge.2) then
basis_vec(1:nn,2)=0.0d0
i=minloc(xi2prec,1)
basis_vec(i,2)=1.0d0
call gsch(basis_vec,basis_vec(1,2),work,nn,nn,1)
call calculate_xi2(m,nmatrix,xdip,ydip,zdip,
& xi1,tr_basis_vec(1,2),
& basis_vec(1,2),xim,q,
& omega1,kq,kx,ky,kz)
nvec=2
endif
if(nn.ge.3) then
basis_vec(1:nn,3)=xi1(1:nn)+
& tr_basis_vec(1:nn,1)/grad_norm
call gsch(basis_vec,basis_vec(1,3),work,nn,nn,2)
call calculate_xi2(m,nmatrix,xdip,ydip,zdip,
& xi1,tr_basis_vec(1,3),
& basis_vec(1,3),xim,q,
& omega1,kq,kx,ky,kz)
nvec=3
endif
else
nvec=nvec+1
do i=1,nn
basis_vec(i,nvec)=kappaold(i)/(xi2prec(i)-mu)
enddo
call gsch(basis_vec,basis_vec(1,nvec),work,nn,nn,nvec-1)
call calculate_xi2(m,nmatrix,xdip,ydip,zdip,
& xi1,tr_basis_vec(1,nvec),
& basis_vec(1,nvec),xim,q,
& omega1,kq,kx,ky,kz)
endif
if(jj.eq.6.or.jj.eq.11) then
if(norm1.gt.0.1d0*res_init)trust_radius=0.5d0*trust_radius
endif
call build_red_hess(basis_vec,tr_basis_vec,
& reduced_hessian,max_micro_it,
& nvec,nn,jj.eq.1)
call ah_method(reduced_hessian,max_micro_it,xi1,
& trust_radius,nvec,kappam,alpha,mu,
& work,nvec+1,grad_norm,
& work(1+(1+nvec)**2+nvec+1),
& 2*max_micro_it**2,work(1+(1+nvec)**2))
c Calculate kappa, and H*kappa
call dgemv('n',nn,nvec,1.0d0,basis_vec,nn,kappam,1,0.0d0,
& kappa,1)
call dgemv('n',nn,nvec,1.0d0,tr_basis_vec,nn,kappam,1,
& 0.0d0,kappaold,1)
c Calculate residual (it is -1*residual!!)
call daxpy(nn,1.d0,xi1,1,kappaold,1)
call daxpy(nn,-mu,kappa,1,kappaold,1)
endif
C
i=idamax(nn,kappaold,1)
norm1=dabs(kappaold(i))
jj=jj+1
if(jj.gt.max_it2) then
write(*,"(' Solution of Newton equation failed.',/,
& ' Iteration not converged.',/,
& ' Program aborted.')")
call mrccend(1)
endif
C
enddo
if(mu.ne.0.0d0.and.lwork.gt.5*nmatrix**2+nmatrix+nn) then
call line_search(m,nmatrix,xdip,ydip,zdip,q,kappa,kx,ky,kz,
& line_search_mem,
& line_search_mem(1+nmatrix**2+nmatrix))
endif
C Calculate the exponent of kappa matrix (umatrix)
call exp_kappa(nmatrix,kappa,ky,kx,kz,kq,work)
C Calculate the product of exp(kappa) matrices
call sum_u(nmatrix,uboys,ky,work)
C Transform orbitals
call transform_integrals(nmatrix,ky,xdip,ydip,zdip,q,work)
call calculate_xi1(m,nmatrix,xdip,ydip,zdip,xi1,xim,q,omega1)
C Check convergence
i=idamax(nn,xi1,1)
norm=dabs(xi1(i))
if (ii.gt.max_it1) then
write(*,"(' Minimization of orbital variances failed.',/,
& ' Iteration not converged.',/,
& ' Program aborted.')")
endif
ii=ii+1
write(*,"(i4,3X,ES12.5E2)") ii,norm
C
enddo
end
C
************************************************************************
subroutine build_red_hess(b,btr,a,lda,m,n,linit)
************************************************************************
* Build reduced hessian matrix for trust region optimization
************************************************************************
implicit none
integer lda, n, m, i, j
double precision b(n,m),a(lda,*),btr(n,m)
double precision ddot
logical linit
if(linit) then
do i = 1, m
do j = i, m
a(i, j) = ddot(n, btr(1, i), 1, b(1, j), 1)
a(j, i) = a(i, j)
enddo
enddo
else
do i = 1, m
a(i, m) = ddot(n, btr(1, m), 1, b(1, i), 1)
a(m, i) = a(i, m)
enddo
endif
end subroutine
************************************************************************
subroutine ah_method(a, lda, grad, h, n, x, alpha, mu,
& aa, ldaa, gn, work, lwork, w)
************************************************************************
* Augemented Hessian method
************************************************************************
implicit none
integer lda, n, ldaa, lwork, phase, i
integer info
double precision a(lda, *), x(*), grad(*), h
double precision aa(n, *), w(*), alpha
double precision mu, mu1, mu2, gn
double precision tol_low, tol_up, work(*)
double precision x_norm
logical cond
double precision dnrm2
integer*4 isyev
equivalence(isyev,info)
logical lhard_case
tol_low = h*0.9d0
tol_up = h
cond = .true.
phase = 1
aa(1:n, 1:n) = a(1:n, 1:n)
call dsyev('v', 'u', n, aa, n, w, work, lwork, info)
! right hand side
do i=1,n
work(i)=-gn*aa(1,i)
enddo
! solving the newton eq.
do i=1,n
x(i)=work(i)/w(i)
enddo
if(dnrm2(n,x,1).lt.tol_up) then
cond=.false.
mu=0.0d0
call dgemv('n',n,n,1.0d0,aa,n,x,1,0.0d0,work,1)
x(1:n)=work(1:n)
return
endif
c check for hard case
if(dabs(work(1)) .lt. 1.0d-7) then
mu=w(1)
w(1:n)=aa(:,1)
aa(1:n,1:n)=a(1:n,1:n)
x(1:n)=0.0d0
x(1)=-gn
call hard_case(n,aa,w,mu,x,work,lwork,tol_low,lhard_case)
if(lhard_case) then
mu=1.0d0
return
endif
endif
c regular case, everything is fine
c recalculating overwritten variables
aa(1:n, 1:n) = a(1:n, 1:n)
call dsyev('v', 'u', n, aa, n, w, work, lwork, info)
! right hand side
do i=1,n
work(i)=-gn*aa(1,i)
enddo
mu1=w(1)
mu2=mu1-10.0d0
mu=mu2
do while(cond)
! solving the lvl shifted newton eq. using diagonalization
do i=1,n
x(i)=work(i)/(w(i)-mu)
enddo
c bisection search
x_norm = dnrm2(n, x, 1)
select case(phase)
case(1) ! searching for the right domain: [mu2, mu1]
if(x_norm .gt. tol_low .and. x_norm .lt. tol_up) then
cond = .false. ! solution found
elseif(x_norm .lt. tol_up) then ! found the right region
mu=0.5d0*(mu1+mu2)
phase = 2
else
mu1=mu
mu2=mu2-1.0d1
mu=mu2
endif
case(2) ! we have the domain, actual bisection search starts
if(x_norm .lt. tol_up .and. x_norm .gt. tol_low) then
cond = .false. ! solution found
elseif(x_norm .gt. tol_up) then
mu1=mu
mu=0.5d0*(mu1+mu2)
else
mu2=mu
mu=0.5d0*(mu1+mu2)
endif
end select
enddo
call dgemv('n',n,n,1.0d0,aa,n,x,1,0.0d0,work,1)
x(1:n)=work(1:n)
end subroutine
************************************************************************
subroutine hard_case(ndim,amat,eigvec,eigval,x,work,lwork,
& trust_radius,lhard_case)
************************************************************************
************************************************************************
implicit none
integer ndim,lwork
double precision amat(ndim,ndim),eigvec(ndim),x(ndim)
double precision work(lwork)
double precision eigval,trust_radius
logical lhard_case
integer rank,info,i
integer*4 isyev
equivalence(isyev,info)
double precision norm,tau
double precision dnrm2
do i=1,ndim
amat(i,i)=amat(i,i)-eigval
enddo
info=0
call dgelss(ndim,ndim,1,amat,ndim,x,ndim,work,1.0d-8,rank,
& work(ndim+1),lwork-ndim,info)
if(info /= 0) then
write(*,'(A)') 'Error in localization!'
write(*,'(A,I10)') 'Error code: ',info
endif
norm = dnrm2(ndim,x,1)
if(norm .lt. trust_radius) then
tau = trust_radius - norm
x = x + tau*eigvec
lhard_case=.true.
else
lhard_case=.false.
endif
end subroutine
C***********************************************************************
C***********************************************************************
C Subroutine to calculate the right side of the Newton equation (xi1
C matrix)
C***********************************************************************
C***********************************************************************
C
subroutine calculate_xi1(
&m,
&nmatrix,
&xdip,
&ydip,
&zdip,
&xi1,
&xi1m,
&q,
&omega1)
C
implicit none
C
C Power of the orbital variances
real*8 m
C
C Size of matrices
integer nmatrix
C
C Dipole moment and second moment matrices
real*8 xdip(nmatrix,*)
real*8 ydip(nmatrix,*)
real*8 zdip(nmatrix,*)
C
C xi1 matrix (right hand of the Newton equation)
real*8 xi1((nmatrix-1)*nmatrix/2)
C
C Work arrays for full xi1m,Q matrices
real*8 xi1m(*),q(nmatrix,nmatrix)
C
C Work array for omega**(m-1) vector
real*8 omega1(nmatrix)
C
integer i,k
C
C Calculate omega**(m-1) vector (sum of the orbital variances)
if(m.eq.1.0d0) then
do i=1,nmatrix
omega1(i)=1.0d0
enddo
elseif(m.eq.2.0d0) then
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)
enddo
else
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)**(m-1.d0)
enddo
endif
C
C Calculate xi1 matrix (antisymmetric). Its upper triangle part is
C ordered in an (nmatrix-1)*nmatrix/2 long vector rowwise.
C
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,K)
do k=1, nmatrix
do i=1,nmatrix
xi1m(i+(k-1)*nmatrix)=-2.0d0*(xdip(i,i)*xdip(i,k)+
& ydip(i,i)*ydip(i,k)+
& zdip(i,i)*zdip(i,k))+q(i,k)
xi1m(i+(k-1)*nmatrix)=2.0d0*m*omega1(i)*xi1m(i+(k-1)*nmatrix)
enddo
enddo
C$OMP END PARALLEL DO
C
C {...}^0
call antisymm(nmatrix,xi1m,xi1)
C
end
C***********************************************************************
C***********************************************************************
C Subroutine to calculate the left side of the Newton equation (xi2
C matrix)
C***********************************************************************
C***********************************************************************
C
subroutine calculate_xi2(
&m,
&nmatrix,
&xdip,
&ydip,
&zdip,
&xi1,
&xi2,
&kappa,
&xi2m,
&q,
&omega1,
&kq,
&kx,
&ky,
&kz)
C
implicit none
C
C Power of orbital variances
real*8 m
C
C Size of matrices
integer nmatrix
C
C Dipole moment and second moment matrices
real*8 xdip(nmatrix,*)
real*8 ydip(nmatrix,*)
real*8 zdip(nmatrix,*)
C
C xi1,xi2,kappa matrices (antisymmetric)
real*8 xi1((nmatrix-1)*nmatrix/2),xi2((nmatrix-1)*nmatrix/2)
real*8 kappa((nmatrix-1)*nmatrix/2)
C
C Work arrays for full xi1,xi2,kappa,Q matrices
real*8 xi2m(*),q(nmatrix,nmatrix)
C
C Work arrays for omega**(m+1) and omega**(m+2) vector (sum of the
C orbital variances)
real*8 omega1(nmatrix)
C
C Work arrays for the matrix products of the q,x,y,z matrices and kappa
C matrix
real*8 kq(nmatrix,nmatrix)
real*8 kx(nmatrix,nmatrix),ky(nmatrix,nmatrix)
real*8 kz(nmatrix,nmatrix)
C
integer i,k,idamax
C
C
C Calculate omega**(m-1) and omega**(m-2) vectors
if(m.eq.1.0d0) then
do i=1,nmatrix
omega1(i)=1.0d0
enddo
elseif(m.eq.2.0d0) then
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)
enddo
else
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)**(m-1.d0)
enddo
endif
C
C Calculate kq, kx, ky, kz temporary matrices
call create_full(nmatrix,kappa,xi2m)
call dsymm('r','u',nmatrix,nmatrix,1.d0,q,nmatrix,xi2m,nmatrix,
& 0.d0,kq,nmatrix)
c kq-qk=kq+(kq)^T
call symmat2(kq,nmatrix)
call dsymm('r','u',nmatrix,nmatrix,1.d0,xdip,nmatrix,xi2m,
& nmatrix,0.d0,kx,nmatrix)
call symmat2(kx,nmatrix)
call dsymm('r','u',nmatrix,nmatrix,1.d0,ydip,nmatrix,xi2m,
& nmatrix,0.d0,ky,nmatrix)
call symmat2(ky,nmatrix)
call dsymm('r','u',nmatrix,nmatrix,1.d0,zdip,nmatrix,xi2m,
& nmatrix,0.d0,kz,nmatrix)
call symmat2(kz,nmatrix)
C
C
C Calculate xi2 matrix (antisymmetric). Its upper triangle part is
C ordered in an (nmatrix-1)*nmatrix/2 long vector rowwise.
C
C -4*sum_x(dia(kx)*x)=-4*sum(dia(kx+(kx)^T)/2*x)=-2*sum(...)
C -2*sum_x(dia(x)*(kx-xk))
C kq-qk
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k,i)
do k=1, nmatrix
do i=1,nmatrix
xi2m(i+(k-1)*nmatrix)=-2.0d0*(kx(i,i)*xdip(i,k)+
& ky(i,i)*ydip(i,k)+
& kz(i,i)*zdip(i,k)+
& kx(i,k)*xdip(i,i)+
& ky(i,k)*ydip(i,i)+
& kz(i,k)*zdip(i,i))+kq(i,k)
xi2m(i+(k-1)*nmatrix)=xi2m(i+(k-1)*nmatrix)*omega1(i)*2.0d0*m
enddo
enddo
C$OMP END PARALLEL DO
if(m.ne.1.0d0) then
if(m.eq.2.0d0) then
do i=1,nmatrix
omega1(i)=kx(i,i)*xdip(i,i)+ky(i,i)*ydip(i,i)+
& kz(i,i)*zdip(i,i)-0.5d0*kq(i,i)
enddo
elseif(m.eq.3.0d0) then
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)*
& (kx(i,i)*xdip(i,i)+ky(i,i)*ydip(i,i)+
& kz(i,i)*zdip(i,i)-0.5d0*kq(i,i))
enddo
else
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)**(m-2.d0)*
& (kx(i,i)*xdip(i,i)+ky(i,i)*ydip(i,i)+
& kz(i,i)*zdip(i,i)-0.5d0*kq(i,i))
enddo
endif
omega1=m*(m-1.0d0)*omega1
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k,i)
do k=1, nmatrix
do i=1,nmatrix
! sum_x dia(x)*x
kx(i,k)=8.0d0*(xdip(i,i)*xdip(i,k)+ydip(i,i)*ydip(i,k)+
& zdip(i,i)*zdip(i,k))-4.0d0*q(i,k)
enddo
enddo
C$OMP END PARALLEL DO
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(k,i)
do k=1, nmatrix
do i=1,nmatrix
xi2m(i+(k-1)*nmatrix)=xi2m(i+(k-1)*nmatrix)+
& omega1(i)*kx(i,k)
enddo
enddo
C$OMP END PARALLEL DO
endif
C
C -1/2*k*xi1
C Create a temporary full matrix (antisymmetric)
call create_full(nmatrix,kappa,kx)
call create_full(nmatrix,xi1,ky)
call dgemm('n','n',nmatrix,nmatrix,nmatrix,-0.5d0,kx,nmatrix,
& ky,nmatrix,1.d0,xi2m,nmatrix)
C {...}^0
call antisymm(nmatrix,xi2m,xi2)
C
end
C***********************************************************************
C***********************************************************************
C Subroutine to calculate the preconditioning matrix (diagonal of
C the Hessian matrix)
C***********************************************************************
C***********************************************************************
C
subroutine calculate_prec(
&m,
&nmatrix,
&xdip,
&ydip,
&zdip,
&xi2prec,
&q,
&omega1)
C
implicit none
C
C Power of orbital variances
real*8 m
C
C Size of matrices
integer nmatrix
C
C Dipole moment and second moment matrices
real*8 xdip(nmatrix,*)
real*8 ydip(nmatrix,*)
real*8 zdip(nmatrix,*)
C
C xi2prec preconditioning matrix
real*8 xi2prec(*)
C
C Work array for full q matrix
real*8 q(nmatrix,*)
C
C Work arrays for omega**(m+1) and omega**(m+2) vector (sum of the
C orbital variances)
real*8 omega1(*)
C
integer i,k,l
C
real*8 temp1,temp2
C
C Calculate omega**(m+1) and omega**(m+2) vectors
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)**(m-1.d0)
enddo
C
C Eq.(37) -> Calculate (xi2prec)_kl,kl matrix element
do k=1,nmatrix
do l=(k+1),nmatrix
C
C 0. Calculate second sums of 4. formula
C temp2=sum_y(y_lk*y_kk)
C temp3=sum_y(y_kl*y_ll)
temp1=(xdip(l,k)*xdip(k,k)+ydip(l,k)*ydip(k,k)+
& zdip(l,k)*zdip(k,k))
temp2=(xdip(k,l)*xdip(l,l)+ydip(k,l)*ydip(l,l)+
& zdip(k,l)*zdip(l,l))
C
C 1. 2*m*(Q_kk-Q_ll)*(Omega_l^(m-1)-Omega_k^(m-1))-
C 8*m*sum_x(x_lk^2)*(Omega_k^(m-1)+Omega_l^(m-1))+
xi2prec(nmatrix*(k-1)-k*(k-1)/2+l-k)=2.d0*m*
&(q(k,k)-q(l,l))*(omega1(l)-omega1(k))-8.d0*m*
&(
& (omega1(k)+omega1(l))*
& (xdip(l,k)**2.d0+ydip(l,k)**2.d0+zdip(l,k)**2.d0)
&)+4.d0*m*
C
C 2. 4*m*Omega_k^(m-1)*sum_x(x_kk^2-x_ll*x_kk)+
C 4*m*Omega_l^(m-1)*sum_x(x_ll^2-x_ll*x_kk)+
&(
& omega1(k)*
& (xdip(k,k)**2.d0+ydip(k,k)**2.d0+zdip(k,k)**2.d0-
& xdip(l,l)*xdip(k,k)-ydip(l,l)*ydip(k,k)-zdip(l,l)*zdip(k,k))+
& omega1(l)*
& (xdip(l,l)**2.d0+ydip(l,l)**2.d0+zdip(l,l)**2.d0-
& xdip(l,l)*xdip(k,k)-ydip(l,l)*ydip(k,k)-zdip(l,l)*zdip(k,k))
&)
enddo
enddo
if(m.ne.1.0d0) then
do k=1,nmatrix
do l=(k+1),nmatrix
do i=1,nmatrix
omega1(i)=(q(i,i)-xdip(i,i)**2.d0-ydip(i,i)**2.d0-
& zdip(i,i)**2.d0)**(m-2.d0)
enddo
xi2prec(nmatrix*(k-1)-k*(k-1)/2+l-k)=
& xi2prec(nmatrix*(k-1)-k*(k-1)/2+l-k)
C 4*m*(m-1)*Q_lk^2*(Omega_k^(m-2)+Omega_l^(m-2))-
&+4.d0*m*(m-1.d0)*q(l,k)**2.d0*
&(
& omega1(k)+omega1(l)
&)-16.d0*m*(m-1.d0)*q(l,k)*
C
C 3. 16*m*(m-1)*Q_lk*
C sum_x(x_lk*(Omega_k^(m-2)*x_kk+Omega_l^(m-2)*x_ll))
&(
& xdip(l,k)*
& (omega1(k)*xdip(k,k)+omega1(l)*xdip(l,l))+
& ydip(l,k)*
& (omega1(k)*ydip(k,k)+omega1(l)*ydip(l,l))+
& zdip(l,k)*
& (omega1(k)*zdip(k,k)+omega1(l)*zdip(l,l))
&)+16.d0*m*(m-1.d0)*
C
C 4. 16*m*(m-1)*Omega_k^(m-2)*sum_x(x_lk*x_kk)*sum_y(y_lk*y_kk)+
C 16*m*(m-1)*Omega_l^(m-2)*sum_x(x_kl*x_ll)*sum_y(y_kl*y_ll)
&(
& omega1(k)*temp1**2.d0+omega1(l)*temp2**2.d0
&)
enddo
enddo
endif
C
C
end
C
C***********************************************************************
C***********************************************************************
C%%% Subroutine to calculate the exponent of a matrix
C***********************************************************************
C***********************************************************************
C
subroutine exp_kappa(
&nmatrix,
&kappa,
&umatrix,
&kappam,
&kappa2,
&work1,work2)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C Kappa matrix
real*8 kappa(*)
C
C Full umatrix=exp(kappa) matrix
real*8 umatrix(nmatrix,*)
C
C Work array for full kappa matrix
real*8 kappam(nmatrix**2),kappa2(nmatrix**2)
C
C Work array for kpower matrix and for dgemm subroutine
real*8 work1(nmatrix**2),work2(nmatrix)
real*8 small_work(10)
C
real*8 dnrm2
real*8 norm,nfactor,nn
C
integer i,j
C
C
C Create a temporary full kappa matrix (antisymmetric)
call create_full(nmatrix,kappa,kappam)
C
c Calculate kappa**2 and its eigendeco (kappa**2=V*eps*V')
call dgemm('n','n',nmatrix,nmatrix,nmatrix,1.0d0,kappam,nmatrix,
& kappam,nmatrix,0.0d0,kappa2,nmatrix)
if(nmatrix.ge.3) then
call dsyev('v','u',nmatrix,kappa2,nmatrix,work2,work1,
& nmatrix**2,i)
else
call dsyev('v','u',nmatrix,kappa2,nmatrix,work2,small_work,10,i)
endif
c b*V'*kappa
call dgemm('t','n',nmatrix,nmatrix,nmatrix,1.0d0,kappa2,nmatrix,
& kappam,nmatrix,0.0d0,umatrix,nmatrix)
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
do i=1,nmatrix
work2(i)=dsqrt(dabs(work2(i)))
if(work2(i).lt.1.0d-7) then
kappam(i)=1.0d0
else
kappam(i)=dsin(work2(i))/work2(i)
endif
enddo
C$OMP END PARALLEL DO
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
do i=1,nmatrix
do j=1,nmatrix
umatrix(j,i)=kappam(j)*umatrix(j,i)
enddo
enddo
C$OMP END PARALLEL DO
c a*V'
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
do i=1,nmatrix
if(work2(i).lt.1.0d-7) then
kappam(i)=1.0d0
else
kappam(i)=dcos(work2(i))
endif
enddo
C$OMP END PARALLEL DO
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
do i=1,nmatrix
do j=1,nmatrix
umatrix(j,i)=umatrix(j,i)+kappam(j)*kappa2(i+(j-1)*nmatrix)
enddo
enddo
C$OMP END PARALLEL DO
c V*(a*V1+b*V1*kappa)
call dcopy(nmatrix**2,umatrix,1,kappam,1)
call dgemm('n','n',nmatrix,nmatrix,nmatrix,1.0d0,kappa2,nmatrix,
& kappam,nmatrix,0.0d0,umatrix,nmatrix)
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to transform orbitals (dipole moment and second moment
C integrals) with umatrix=exp(kappa)
C***********************************************************************
C***********************************************************************
C
subroutine transform_integrals(
&nmatrix,
&umatrix,
&xdip,
&ydip,
&zdip,
&q,
&work)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C Full umatrix=exp(kappa) matrix to transform the orbitals (dipole
C moment and second oment integrals)
real*8 umatrix
C
C Dipole moment and second moment integrals (input/outputs)
real*8 xdip
real*8 ydip
real*8 zdip
real*8 q
C
C Work array for transform_integral subroutine
real*8 work
C
C
call transform_integral(nmatrix,umatrix,xdip,work)
call transform_integral(nmatrix,umatrix,ydip,work)
call transform_integral(nmatrix,umatrix,zdip,work)
call transform_integral(nmatrix,umatrix,q,work)
C
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to make a step of the Gauss-Seidel method (iteration)
C***********************************************************************
C***********************************************************************
C
subroutine gauss_seidel(
&nmatrix,
&kappa,
&xi1,
&xi2,
&xi2prec)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C kappa,xi1,xi2,xi2prec matrices
real*8 kappa(*),xi1(*),xi2(*),xi2prec(*)
real*8 sigma
C
integer i,nn
nn=(nmatrix-1)*nmatrix/2
do i=1,nn
sigma=xi2(i)-xi2prec(i)*kappa(i)
kappa(i)=(-xi1(i)-sigma)/xi2prec(i)
enddo
end subroutine
C***********************************************************************
C***********************************************************************
C Subroutine to make a step of the Jacobian method (iteration)
C***********************************************************************
C***********************************************************************
C
subroutine jacobi(
&nmatrix,
&kappa,
&xi1,
&xi2,
&xi2prec)
C
C One step of the Jacobian method (iteration)
C
C x_new = A_0^(-1)*b-A_0^(-1)*A*x_old+E*x_old
C
C A_0 -> preconditioning matrix, xi2prec (four indexed, diagonal)
C b -> -xi1
C x -> kappa
C A*x -> xi2
C E -> unit matrix
C
implicit none
C
C Size of matrices
integer nmatrix
C
C Coeffitient to control the step size of the Jacobian overrelaxation
C method. If omega=1.d0 then it is the original Jacobian method
real*8 omega
parameter (omega=1.d0)
C
C kappa,xi1,xi2,xi2prec matrices
real*8 kappa(*),xi1(*),xi2(*),xi2prec(*)
C
integer i
C
C
C Do one step of the Jacobian method (iteration)
do i=1,(nmatrix-1)*nmatrix/2
kappa(i)=kappa(i)-omega*(1.d0/xi2prec(i)*xi1(i)+
& 1.d0/xi2prec(i)*xi2(i))
enddo
C
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to calculate the cumulated full usum matrix, in which
C the product of the umatrix=exp(kappa) transformation matrices is
C stored.
C***********************************************************************
C***********************************************************************
C
subroutine sum_u(
&nmatrix,
&usum,
&umatrix,
&work)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C Cumulated full usum matrix, in which the product of the
C umatrix=exp(kappa) transformation matrices is stored
real*8 usum
C
C Full umatrix=exp(kappa) matrix to transform the dipole moment and
C second moment matrices every step of the macroiteration
real*8 umatrix
C
C Work array for dgemm subroutine
real*8 work
C
C
call dgemm('n','n',nmatrix,nmatrix,nmatrix,1.d0,umatrix,nmatrix,
& usum,nmatrix,0.d0,work,nmatrix)
call dcopy(nmatrix**2,work,1,usum,1)
C
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to create full matrix from the upper triangle part of
C an antisymmetric matrix (rowwise ordered)
C***********************************************************************
C***********************************************************************
C
subroutine create_full(
&nmatrix,
&x,
&xfull)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C The vector to be transformed
real*8 x(*)
C
C The full matrix
real*8 xfull(nmatrix,*)
C
integer i,j
C
C
C Create the full antisymmetric matrix (the diagonal is zero)
do i=1,nmatrix
do j=1,nmatrix
if (i.eq.j) then
xfull(i,j)=0.d0
elseif (i.lt.j) then
xfull(i,j)=x(nmatrix*(i-1)-i*(i-1)/2+j-i)
xfull(j,i)=-x(nmatrix*(i-1)-i*(i-1)/2+j-i)
endif
enddo
enddo
C
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to modify preconditioning matrix if it
C is singular
C***********************************************************************
C***********************************************************************
C
subroutine modify_prec(
&nmatrix,
&xi2prec)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C xi2prec preconditioning matrix
real*8 xi2prec(*)
C
C Treshold for elements of the preconditioning matrix
C%%% Szerintem globalisba kene!!!
real*8 tresh
parameter (tresh=1.d-5)
C
integer i
C
C
do i=1,(nmatrix-1)*nmatrix/2
if ((xi2prec(i).le.0.d0).and.(xi2prec(i).ge.(-1.d0*tresh))) then
xi2prec(i)=xi2prec(i)-tresh
elseif ((xi2prec(i).ge.0.d0).and.(xi2prec(i).le.tresh))then
xi2prec(i)=xi2prec(i)+tresh
endif
enddo
C
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to create full diagonal matrix from a vector or the
C diagonal elements of a full matrix
C***********************************************************************
C***********************************************************************
C
subroutine diag(
&nmatrix,
&job,
&array,
&diagonal)
C
implicit none
C
C Size of full diagonal matrix
integer nmatrix
C
C Job control: if (job.eq.0) then array is a vector, if (job.ne.0) then
C array is a full matrix
integer job
C
C Array to diagonalize
real*8 array
C
C Diagonal matrix
real*8 diagonal
C
C
call dfillzero(diagonal,nmatrix**2)
if (job.eq.0) then
call dcopy(nmatrix,array,1,diagonal,nmatrix+1)
else
call dcopy(nmatrix,array,nmatrix+1,diagonal,nmatrix+1)
endif
C
end
C
C***********************************************************************
C***********************************************************************
C Subroutine to calculate the odd component of a full matrix:
C {A}°=A-A+
C Only the upper triangle part of the antisymmetric matrix {A}° is
C stored in a vector ordered in rowwise.
C***********************************************************************
C***********************************************************************
C
subroutine antisymm(
&nmatrix,
&matrixin,
&matrixout)
C
implicit none
C
C Size of matrices
integer nmatrix
C
C Full matrix on input
real*8 matrixin(nmatrix,*)
C
C Antisymmetrized matrix. Only the upper triangle part of tha matrix is
C stored in this vector, ordered rowwise
real*8 matrixout(*)
C
integer i,j
C
C
do i=1,nmatrix-1
do j=i+1,nmatrix
matrixout(nmatrix*(i-1)-i*(i-1)/2+j-i)=matrixin(i,j)-
& matrixin(j,i)
enddo
enddo
C
end
************************************************************************
subroutine line_search(m,nmatrix,xdip,ydip,zdip,q,
&kappa,umatrix,kappam1,kappam2,temp,work)
************************************************************************
* Line search for localization
************************************************************************
implicit none
integer nmatrix
real*8 m
real*8 xdip(nmatrix,nmatrix),ydip(nmatrix,nmatrix)
real*8 zdip(nmatrix,nmatrix),q(nmatrix,nmatrix)
real*8 kappa((nmatrix-1)*nmatrix/2)
real*8 umatrix(nmatrix,nmatrix)
real*8 kappam1(nmatrix,nmatrix)
real*8 kappam2(nmatrix,nmatrix)
real*8 temp(nmatrix*nmatrix+nmatrix)
real*8 work(4*nmatrix**2+(nmatrix-1)*nmatrix/2)
integer idipx2,idipy2,idipz2,iq2,ikappa2,nn,i,j
integer ntrial
real*8 alpha0,alpha
parameter(alpha0=0.5d0,ntrial=20)
real*8 vals(ntrial)
idipx2=1
idipy2=idipx2+nmatrix**2
idipz2=idipy2+nmatrix**2
iq2=idipz2+nmatrix**2
ikappa2=iq2+nmatrix**2
nn=(nmatrix-1)*nmatrix/2
call dcopy(nmatrix**2,xdip,1,work(idipx2),1)
call dcopy(nmatrix**2,ydip,1,work(idipy2),1)
call dcopy(nmatrix**2,zdip,1,work(idipz2),1)
call dcopy(nmatrix**2,q,1,work(iq2),1)
call dcopy(nn,kappa,1,work(ikappa2),1)
call dscal(nn,alpha0,work(ikappa2),1)
do i=1,ntrial
call exp_kappa(nmatrix,work(ikappa2),umatrix,kappam1,kappam2,
& temp,temp(1+nmatrix**2))
call transform_integrals(nmatrix,umatrix,work(idipx2),
& work(idipy2),work(idipz2),work(iq2),temp)
vals(i)=0.0d0
do j=1,nmatrix
vals(i)=vals(i)+work(iq2+(j-1)*(1+nmatrix))
& -work(idipx2+(j-1)*(1+nmatrix))**2
& -work(idipy2+(j-1)*(1+nmatrix))**2
& -work(idipz2+(j-1)*(1+nmatrix))**2
enddo
if(m.ne.1.0d0) vals(i)=vals(i)**m
enddo
i=minloc(vals,1)
c write(*,'(A,f)') 'Line search: ', dble(i)*alpha0
call dcopy(nn,work(ikappa2),1,kappa,1)
call dscal(nn,dble(i),kappa,1)
c call exp_kappa(nmatrix,work(ikappa2),umatrix,kappam1,kappam2,
c & temp,temp(1+nmatrix**2))
c call transform_integrals(nmatrix,umatrix,xdip,ydip,zdip,q,temp)
C
end
************************************************************************
subroutine boys_tol(lsym,tol)
************************************************************************
* calculating tolerance
************************************************************************
double precision tol
logical lsym
integer i
character*4 scftol,localcorrsymm
character*8 pntgrp
call getkey('scftol',6,scftol,4)
read(scftol,*) i
tol=min(1.d-9,10.d0**(-i))
call getkey('localcorrsymm',13,localcorrsymm,4)
call getvar('pntgrp ',pntgrp)
lsym=.false.
if((localcorrsymm.ne.'off '.and.pntgrp.ne.'C1 ').or.
$ localcorrsymm.eq.'forc') then
tol=tol/1000.d0
lsym=.true.
endif
end subroutine
Reference in New Issue View Git Blame Copy Permalink