cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c utils.f c RGF c modular version c Revision 2, Jul 5, 2001 c Revision 2.2 for TASC: Jul 23, 2001 c c Michael Hirsch / Wolfgang Quapp c Uni. Leipzig subroutine vec_dcopy(arra1,arra2,ndim) c Copies ARRA1 to ARRA2 real*8 arra1,arra2 integer i,ndim dimension arra1(ndim),arra2(ndim) if(ndim .le. 0) return do 10 i = 1,ndim arra2(i) = arra1(i) 10 continue return end subroutine vec_dmult_constant(arra1,ndim,zahl,arra2) c ARRA2 = ZAHL * ARRA1 real*8 arra1,arra2,zahl integer ndim,i dimension arra1(ndim),arra2(ndim) if(ndim .le. 0)return do i = 1,ndim arra2(i) = arra1(i)*zahl end do return end subroutine mat_dcopy(arra1,arra2,ndim,mdim) c Copies matrix ARRA1 to matrix ARRA2 real*8 arra1,arra2 integer i,j,ndim,mdim dimension arra1(ndim,mdim),arra2(ndim,mdim) do 10,i = 1,ndim do 20,j = 1,mdim arra2(i,j) = arra1(i,j) 20 continue 10 continue return end subroutine vec_dinit(arra1,ndim,wert) c Set all elements of ARRA1 to the double value of WERT real*8 arra1,wert integer ndim,i dimension arra1(ndim) if(ndim .le. 0) return do i=1,ndim arra1(i)=wert end do return end subroutine vec_iinit(iarray,ndim,iwert) c Set all elements of ARRA1 to the integer value of WERT integer iarray,iwert,ndim,i dimension iarray(ndim) if(ndim .le. 0) return do i=1,ndim iarray(i)=iwert end do return end subroutine vec_oinit(oarray,ndim,owert) c Set all elements of ARRA1 to the logical value of WERT integer ndim,i logical oarray,owert dimension oarray(ndim) if(ndim .le. 0) return do i=1,ndim oarray(i)=owert end do return end subroutine mat_diag(d,n,v) c create a NxN diagonal matrix D with entries V real*8 d,v integer n,i dimension d(n,n) if(n .le. 0) return call vec_dinit(d,n*n,0.0d0) do i=1,n d(i,i)=v end do return end subroutine vec_dmult_add(arra1,arra2,ndim,zahl,arra4) c For each Komponent: ARRA4 = ARRA1 + ZAHL*ARRA2 real*8 arra1,arra2,arra4,zahl integer ndim,i dimension arra1(ndim),arra2(ndim),arra4(ndim) if(ndim .le. 0)return do 145 i = 1,ndim arra4(i) = arra1(i)+arra2(i)*zahl 145 continue return end subroutine tri_mat(dm1,nr1,nc1,dm2,nc2,dm3,nc3,dres) c dm1[nr1,nc1]*dm2[nc1,nc2]*dm3[nc2,nc3]=dres[nr1,nc3] real*8 dm1,dm2,dm3,dres,r1 integer nr1,nc1,nc2,nc3 dimension dm1(nr1,nc1),dm2(nc1,nc2),dm3(nc2,nc3), + r1(nr1,nc2),dres(nr1,nc3) call matmult(dm1,nr1,dm2,nc2,r1,nc1) call matmult(r1,nr1,dm3,nc3,dres,nc2) return end subroutine vec_mat_centry(vec,mat,npos,ndim,ncol,zkind) c Copies VEC 'i'n or 'o'ut the NPOS column of matrix MAT real*8 vec,mat integer npos,ndim,ncol,i character zkind dimension vec(ndim),mat(ndim,ncol) if(npos.gt.ncol) return if (zkind.eq.'i') then do i=1,ndim mat(i,npos)=vec(i) end do endif if (zkind.eq.'o') then do i=1,ndim vec(i)=mat(i,npos) enddo endif return end subroutine vec_mat_rentry(vec,dmat,npos,ndim,nrow,zkind) c Copies VEC 'i'n or 'o'ut the NPOS row of matrix MAT real*8 vec,dmat integer npos,ndim,nrow,i character zkind dimension vec(ndim),dmat(nrow,ndim) if(npos.gt.nrow) return if (zkind.eq.'i') then do 10,i=1,ndim 10 dmat(npos,i)=vec(i) endif if (zkind.eq.'o') then do 20,i=1,ndim 20 vec(i)=dmat(npos,i) endif return end real*8 function vec_ddot(arra1,arra2,ndim) c cartesian scalar product of ARRA1 and ARRA2 real*8 arra1,arra2 integer ndim,i dimension arra1(ndim),arra2(ndim) vec_ddot = 0 if(ndim .le. 0)return do 70 i = 1,ndim vec_ddot = vec_ddot + arra1(i) * arra2(i) 70 continue return end real*8 function adet(mat,sdet,dim) c calculate the Determinant of MAT real*8 mat,amat,det,fact,sdet,sz integer dim,itest,sp,nexp,k,dreieck,ll dimension mat(dim,dim),amat(dim,dim),sp(dim-1),fact(dim) sz=1.0d0/16.0d0 call mat_dcopy(mat,amat,dim,dim) itest=dreieck(amat,sp,dim) do 5,k=1,dim-1 5 continue if (itest.eq.-1) then adet=0.0d0 sdet=0.0d0 else nexp=0 det=1 sdet=1.0d0 do 10,k=1,dim det=det*amat(k,k) fact(k)=amat(k,k) if (k.lt.dim) then if (k.lt.sp(k)) det=-det endif 20 continue if (dabs(det).lt.sz) then det=det*16.0d0 nexp=nexp-1 goto 20 endif 30 continue if (dabs(det).ge.1.0d0) then det=det*sz nexp=nexp+1 goto 30 endif if (k.eq.dim-1) then sdet=det ll=nexp endif 10 continue if (nexp.lt.0) then nexp=-nexp do 40,k=1,nexp det=det*sz 40 continue elseif (nexp.gt.0) then do 50,k=1,nexp det=det*16.0d0 50 continue endif if (ll.lt.0) then ll=-ll do 410,k=1,ll sdet=sdet*sz 410 continue elseif (ll.gt.0) then do 510,k=1,ll sdet=sdet*16.0d0 510 continue endif adet=det endif return end subroutine recorrect(corr,cold,gmat,ocorr,nn) c manipulation of the corrector for a better convergence integer nn real*8 corr,cold,dmnorm,gmat,vec_mdot logical olen,odir,ocorr,odone dimension corr(nn),cold(nn),gmat(nn,nn) ocorr=.false. if (dmnorm(corr,gmat,nn).gt.dmnorm(cold,gmat,nn)) then olen=.true. else olen=.false. endif if (vec_mdot(corr,gmat,cold,nn).lt.0.0d0) then odir=.true. else odir=.false. endif odone=.false. if (odir.and.olen.and.ocorr) then call vec_dmult_constant(corr,nn,0.6d0,corr) odone=.true. endif ocorr=odone return end integer function dreieck(mat,sp,dim) c Dreiecksfaktorisierung of MAT ; dim(SP) = DIM-1 real*8 mat,z integer sp,dim,it,k,i,s,j dimension mat(dim,dim),sp(dim-1) it=0 do 10,k=1,dim-1 z=0.0d0 do 20,i=k,dim if (dabs(mat(i,k)).gt.z) then z=dabs(mat(i,k)) s=i endif 20 continue if (z.eq.0.0d0) then dreieck=-1 return endif sp(k)=s if (k.lt.s) then do 30,j=1,dim z=mat(k,j) mat(k,j)=mat(s,j) mat(s,j)=z 30 continue endif do 40,i=k+1,dim mat(i,k)=mat(i,k)/mat(k,k) 40 continue do 50,j=k+1,dim do 60,i=k+1,dim mat(i,j)=mat(i,j)-mat(i,k)*mat(k,j) 60 continue 50 continue 10 continue if (mat(dim,dim).eq.0.0d0) it=-1 dreieck=it return end subroutine matadd(dm1,fact,dm2,out,n,m) c OUT(N,M) = DM1(N,M) + FACT * DM1(N,M) real*8 dm1,fact,dm2,out integer m,n dimension dm1(n,m),dm2(n,m),out(n,m) call vec_dmult_add(dm1,dm2,n*m,fact,out) return end subroutine matmult(m1,d1,m2,d2,mout,dim) c matrix multiplication: m1[d1,dim]*m2[dim,d2]=mout[d1,d2] real*8 m1,m2,mout,swap integer d1,d2,dim,i,j,k dimension m1(d1,dim),m2(dim,d2),mout(d1,d2),swap(d1,d2) call vec_dinit(swap,d1*d2,0.0d0) do 10,i=1,d1 do 20,j=1,d2 do 30,k=1,dim swap(i,j)=swap(i,j)+m1(i,k)*m2(k,j) 30 continue 20 continue 10 continue call vec_dcopy(swap,mout,d1*d2) return end real*8 function det(matein,nn) c calculation of the determinant of MATEIN integer i,nn real*8 matein,mataus,w dimension matein(nn,nn),mataus(nn,nn),w(nn) call householder(matein,mataus,w,nn,nn) det=1.0d0 do 10,i=1,nn if (dabs(mataus(i,i)).lt.1.0d-25) then det=0.0d0 goto 20 endif if (dabs(det) .lt.1.0d-25) then det=0.0d0 goto 20 endif 10 det=det*mataus(i,i) if (iand(nn,1).eq.0) * det=-det 20 continue return end subroutine ortcompl(vec,compl,dmet,ndim) c calculate an Orthonormalcomplement COMPL to vector VEC real*8 basis,vec,compl,swap,det,dmgs,dmet integer i,j,ndim cc integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene cc common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene dimension vec(ndim),compl(ndim-1,ndim),basis(ndim,ndim), *swap(ndim),dmet(ndim,ndim) call vec_dinit(basis,ndim*ndim,0.0d0) do 10,i=1,ndim basis(i,i)=1.0d0 10 continue do 20,i=1,ndim if (vec(i).ne.0.0d0) then do 30,j=1,ndim swap(j)=basis(j,1) basis(j,1)=vec(j) if (i.ne.1) basis(j,i)=swap(j) 30 continue goto 40 endif 20 continue 40 continue det=dmgs(basis,dmet,ndim) call vec_dinit(compl,(ndim-1)*ndim,0.0d0) call mat_trans(basis,ndim,ndim,basis) do 160,i=2,ndim call vec_mrowcopy(i,basis,ndim,i-1,compl,ndim-1,ndim) 160 continue return end subroutine mat_trans(matin,nrow,ncol,matout) c transpose MATOUT of matrix MATIN real*8 matin,matout,swap integer nrow,ncol,i,j dimension matin(nrow,ncol),matout(ncol,nrow),swap(ncol,nrow) do 10,i=1,nrow do 20,j=1,ncol swap(j,i)=matin(i,j) 20 continue 10 continue call mat_dcopy(swap,matout,ncol,nrow) return end real*8 function dmgs(a,dmet,dim) c modified Gram-Schmidt algorithm return det(a) real*8 det,rr,a,q,r,vec_mmdot,dmet,dzero integer dim,i,k dimension a(dim,dim),q(dim,dim),r(dim,dim),dmet(dim,dim) data dzero/1.0d-55/ call vec_dinit(q,dim*dim,0.0d0) call vec_dinit(r,dim*dim,0.0d0) det=1.0d0 call mat_dcopy(a,q,dim,dim) do 10,k=1,dim if (k.eq.1) goto 21 do 20,i=1,k-1 r(i,k)=vec_mmdot(i,q,dim,k,q,dim,dim,dmet) call vec_madd(k,q,dim,1.0d0,i,q,dim,-r(i,k),k,q,dim,dim) 20 continue 21 continue rr=vec_mmdot(k,q,dim,k,q,dim,dim,dmet) r(k,k)=dsqrt(rr) det=det*r(k,k) if(rr.lt.dzero) goto 33 call vec_madd(k,q,dim,1.0d0/r(k,k),k,q,dim,0.0d0,k,q,dim,dim) 10 continue call mat_dcopy(q,a,dim,dim) dmgs=det return 33 CONTINUE call cerr(-23) return end subroutine vec_madd c (nr1,mat1,d1,fact1,nr2,mat2,d2,fact2,nrx,matx,dx,dim) c matx(nrx)=fact1*mat1(nr1)+fact2*mat2(nr2) c d1,d2,dx - columns c dim - rows real*8 mat1,mat2,matx,fact1,fact2,vec integer nr1,nr2,nrx,d1,d2,dx,dim,i dimension mat1(dim,d1),mat2(dim,d2),matx(dim,dx),vec(dim,1) call vec_dinit(vec,dim,0.0d0) do 10,i=1,dim vec(i,1)=fact1*mat1(i,nr1)+fact2*mat2(i,nr2) 10 continue call vec_mat_centry(vec,matx,nrx,dim,dx,'i') return end real*8 function vec_mmdot &(nri,matin,din,nra,matout,dout,dim,dmet) c multiply the (nri) col of (matin) and the (nra) col of c (matout). (dim) is the number of the rows with metric dmet real*8 matin,matout,prod,v1,v2,dmet,vec_mdot integer nri,din,nra,dout,dim dimension matin(dim,din),matout(dim,dout),dmet(dim,dim), + v1(dim),v2(dim) if ((nri.gt.din).or.(nra.gt.dout)) call cerr(-100) call vec_mat_centry(v1,matin,nri,din,dim,'o') call vec_mat_centry(v2,matout,nra,dout,dim,'o') prod=vec_mdot(v1,dmet,v2,dim) vec_mmdot=prod return end subroutine vec_mrowcopy(nri,matin,din,nra,matout,dout,dim) c copies the (nri) row of (matin) with (din) row and (dim) col c to the (nra) row of (matout) with (dout) row and (dim) col real*8 matin,matout integer nri,din,nra,dout,dim,i dimension matin(din,dim),matout(dout,dim) if ((nri.gt.din).or.(nra.gt.dout)) call cerr(17) do 10,i=1,dim matout(nra,i)=matin(nri,i) 10 continue return end subroutine householder(matein,mataus,w,n,sp) c householder transformation integer k,l,m,n,sp real*8 alpha,rho,sk,s,w,matein,mataus dimension matein(n,sp),mataus(n,sp),w(n) call vec_dcopy(matein,mataus,n*sp) do 10,k=1,n-1 s=0.0d0 do 20,l=k,n s=s+mataus(l,k)*mataus(l,k) 20 continue if (mataus(k,k).lt.0) then alpha=dsqrt(s) else alpha=-dsqrt(s) endif rho=dsqrt(s-mataus(k,k)*mataus(k,k)+ * (mataus(k,k)-alpha)*(mataus(k,k)-alpha)) if(rho .lt. 1.0d-63) then rho= 1.0d-63 C w(k)= 0.0d0 C goto 22 endif w(k)=(mataus(k,k)-alpha)/rho 22 continue do 30,l=k+1,n w(l)=mataus(l,k)/rho 30 continue do 40,l=k,sp sk=0.0d0 do 50,m=k,n sk=sk+w(m)*mataus(m,l) 50 continue do 60,m=k,n mataus(m,l)=mataus(m,l)-2*sk*w(m) 60 continue 40 continue 10 continue return end subroutine closys c closes internal files of mod. RGF: "store.rgf" "param.rgf" "messages.rgf" integer iochk integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene close(isto,iostat=iochk) close(ipar,iostat=iochk) close(imsg,iostat=iochk) return end subroutine print_switch(io) c Output of the switches integer io logical o_flat, o_free, o_msg, o_no_up, c o_const_sgn, o_even, o_recorr, o_nn character*40 zkey common/switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn call keytable(zkey,14) if(o_flat)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,15) if(o_free)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,16) if(o_msg)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,17) if(o_no_up)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,18) if(o_const_sgn)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,19) if(o_even)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,20) if(o_recorr)then write(io,10)zkey else write(io,20)zkey endif call keytable(zkey,21) if(o_nn)then write(io,10)zkey else write(io,20)zkey endif return 10 format(1x,A25,': On') 20 format(1x,A25,': Off') end subroutine toang(vec1,vec2,nn) c (bohr,radian)->(angstroem,degree) for internal coordinates real*8 vec1,vec2,torad,btoa integer i,nr,nn dimension vec1(nn),vec2(nn) data btoa /0.52917706d0/ torad=dacos(-1.0d0)/180.0d0 nr=(nn+6)/3-1 do i=1,nn if (i.le.nr) then vec2(i)=vec1(i)*btoa else vec2(i)=vec1(i)/torad endif end do return end subroutine message(zstr,lstr) c writes the message ZSTR with length LSTR into imsg ("message.rgf") character*40 zstr integer lstr integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene logical o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even,o_recorr, c o_nn common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene if (o_msg) write(imsg,zstr(1:lstr)) return end subroutine err(iochk,zsource,lenz) c errorhandling integer iochk,lenz integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene character*40 zsource logical o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene if (iochk.ne.0) then if (o_msg) then write(imsg,10)zsource(1:lenz) else write(ires,10)zsource(1:lenz) endif endif call cerr(iochk) return 10 format(1x,'Error-Source: ',A) end subroutine cerr(iochk) c errorhandling character*40 zkey integer iochk integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene logical o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene c errorhandling: Number of error - IOCHK, outputfile - IOFILE if ((.not.o_msg).and.(iochk.ne.0)) call exit (1) select case (iochk) case (0) return case (1:40) write(imsg,10)iochk case (-1) write(imsg,20) case (-29:-10) call errtable(zkey,iochk) write(imsg,30)zkey case (-40:-30) call errtable(zkey,iochk) write(imsg,50)zkey case default write(imsg,40)iochk end select write(imsg,100) call closys call closext call exit(1) return 10 format(1x,'IO Error with number ',i3,' occured') 20 format(1x,'Unexpected End-Of_File occoured') 30 format(1x,'IO Error occured. Source: ',A40) 40 format(1x,'strange IO Error occured. number:',i4) 50 format(1x,'Internal numeric Error occured:',/,A40) 100 format(1x,'*** Exit program with status 1***') end subroutine openext c check and open files for output, point, gradinet, hessian, b-matrix, c and energy integer iochk integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene character*40 z_store,z_para,z_msg,z_FoR,z_FoP,z_FoG,z_FoH,z_FoB, c z_iostat,z_FoE logical o_flat,opened,o_msg,o_free,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common /switch/o_flat,o_free,o_msg,o_no_up,o_const_sgn, c o_even,o_recorr,o_nn common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene common/iofiles/z_store,z_para,z_msg,z_FoR,z_FoP, c z_FoG,z_FoH,z_FoB,z_FoE opened=.false. if (o_msg) then inquire(file=z_msg,opened=opened) if (.not.opened) c open(imsg,file=z_msg,status='unknown',access='append') opened=.true. endif open (ires,file=z_FoR,status='unknown',access='append', c err=10,iostat=iochk) if (opened) write(imsg,*)z_FoR//z_iostat(iochk) open (ipnt,file=z_FoP,err=10,status='old',iostat=iochk) if (opened) write(imsg,*)z_FoP//z_iostat(iochk) open (igrd,file=z_FoG,err=10,status='old',iostat=iochk) if (opened) write(imsg,*)z_FoG//z_iostat(iochk) open (ihes,file=z_FoH,err=10,status='old',iostat=iochk) if (opened) write(imsg,*)z_FoH//z_iostat(iochk) if (.not.o_flat) then open (iwil,file=z_FoB,err=10,status='old',iostat=iochk) if (opened) write(imsg,*)z_FoB//z_iostat(iochk) endif open (iene,file=z_FoE,err=10,status='old',iostat=iochk) if (opened) write(imsg,*)z_FoE//z_iostat(iochk) return 10 continue call cerr(-12) return end subroutine closext c close files for output, point, gradinet, hessian, b-matrix integer iochk integer isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene common/iosys/isto,ipar,imsg,ires,ipnt,igrd,ihes,iwil,iene close(ires,iostat=iochk) close(ipnt,iostat=iochk) close(igrd,iostat=iochk) close(ihes,iostat=iochk) close(iwil,iostat=iochk) close(iene,iostat=iochk) return end real*8 function vec_mdot(vec1,metric,vec2,dim) c metric scalar product of VEC1 and VEC2 with respect to METRIC integer dim real*8 vec1,metric,vec2,cov2,vec_ddot dimension vec1(dim),vec2(dim),cov2(dim),metric(dim,dim) call matmult(metric,dim,vec2,1,cov2,dim) vec_mdot=vec_ddot(vec1,cov2,dim) return end real*8 function diag_inv(diag,dinv,nn) c inverse DINV of a diagonalmatrix DIAG real*8 diag,dinv,det integer nn,i dimension diag(nn,nn),dinv(nn,nn) call vec_dinit(dinv,nn*nn,0.0d0) det=1.0d0 do i=1,nn if (diag(i,i).eq.(0.0d0)) call cerr(-30) det=det*diag(i,i) dinv(i,i)=1.0d0/diag(i,i) enddo diag_inv=det return end logical function omat_inv(in,ndim,out) c calculates the inverse matrix of IN real*8 in,out integer ndim,i,j,k dimension in(ndim,ndim),out(ndim,ndim) call mat_dcopy(in,out,ndim,ndim) do 10,k=1,ndim if (out(k,k).eq.0.0d0) then omat_inv=.false. return endif out(k,k)=1.0d0/out(k,k) do 20,i=1,ndim if (i.ne.k) out(i,k)=-(out(i,k)*out(k,k)) 20 continue do 30,j=1,ndim if (j.ne.k) then do 40,i=1,ndim if (i.ne.k) out(i,j)=out(i,j)+out(i,k)*out(k,j) 40 continue out(k,j)=out(k,j)*out(k,k) endif 30 continue 10 continue omat_inv=.true. return end subroutine getmet(bmat,bpinv,gmat,ginv,nn) c calculation of GMAT and GINV real*8 bmat,bpinv,gmat,ginv integer nn logical omat_inv dimension bmat(nn,nn+6),gmat(nn,nn),ginv(nn,nn),bpinv(nn,nn+6) c build metric, its inverse, the pseudo inverse of B-matrix from B-matrix call vec_dinit(gmat,nn*nn,0.0d0) call blas_dgemm('n','t',nn,nn,nn+6,1.0d0,bmat,nn, c bmat,nn,1.0d0,gmat,nn) if (.not.omat_inv(gmat,nn,ginv)) call err(-32,'getmet',6) call vec_dinit(bpinv,nn*(nn+6),0.0d0) call blas_dgemm('n','n',nn,nn+6,nn,1.0d0 c ,ginv,nn,bmat,nn,1.0d0,bpinv,nn) return end subroutine blas_dgemm c matrix multiplikation * (mata,matb,nrow,ncol,mdi,zahl1,arra1, * ndima,arra2,ndimb,zahl2,arra3,ndimc) character*1 mata,matb integer*4 ncol,nrow,ndima,ndimb,ndimc,mdi real*8 zahl1,zahl2,arra1 real*8 arra2,arra3 real*8 zwischen integer*4 i,j,k dimension arra1(*),arra2(*),arra3(*) if (mdi.eq.0) then do 1,i=1,nrow do 2,j=1,ncol arra3(i+(j-1)*ndimc)=zahl2*arra3(i+(j-1)*ndimc) 2 continue 1 continue return endif if (mata.eq.'n'.and.matb.eq.'n') then do 3,i=1,nrow do 4,j=1,ncol zwischen = 0.0 do 5,k=1,mdi zwischen = arra1(i+(k-1)*ndima) * c arra2(k+(j-1)*ndimb) + zwischen 5 continue arra3(i+(j-1)*ndimc) = zahl1*zwischen + c zahl2*arra3(i+(j-1)*ndimc) 4 continue 3 continue return endif if (mata.eq.'n'.and.matb.eq.'t') then do 6,i=1,nrow do 7,j=1,ncol zwischen = 0.0 do 8,k=1,mdi zwischen = arra1(i+(k-1)*ndima) * c arra2(j+(k-1)*ndimb) + zwischen 8 continue arra3(i+(j-1)*ndimc) = zahl1*zwischen + c zahl2*arra3(i+(j-1)*ndimc) 7 continue 6 continue return endif if (mata.eq.'t'.and.matb.eq.'n') then do 9,i=1,nrow do 10,j=1,ncol zwischen = 0.0 do 11,k=1,mdi zwischen = arra1(k+(i-1)*ndima) * c arra2(k+(j-1)*ndimb) + zwischen 11 continue arra3(i+(j-1)*ndimc) = zahl1*zwischen + c zahl2*arra3(i+(j-1)*ndimc) 10 continue 9 continue return endif if (mata.eq.'t'.and.matb.eq.'t') then do 12,i=1,nrow do 13,j=1,ncol zwischen = 0.0 do 14,k=1,mdi zwischen = arra1(k+(i-1)*ndima) * c arra2(j+(k-1)*ndimb) + zwischen 14 continue arra3(i+(j-1)*ndimc) = zahl1*zwischen + c zahl2*arra3(i+(j-1)*ndimc) 13 continue 12 continue return endif print*,'blas_dgemm failed' stop return end subroutine buildk(dk,req,tang,rgr,nn) c build the K-Matrix real*8 dk,req,tang,rgr integer nn,i,j dimension dk(nn,nn+1),req(nn-1,nn),tang(nn),rgr(nn-1) do 10,i=1,nn-1 do 11,j=1,nn dk(i,j)=req(i,j) 11 continue 10 continue do 20,j=1,nn dk(nn,j)=tang(j) 20 continue do 30,i=1,nn-1 dk(i,nn+1)=-rgr(i) 30 continue dk(nn,nn+1)=0.0d0 return end subroutine nka(ginv,hi,eival,eivec,nn) c normal coordinate analysis real*8 schol,ginv,hi,tmp1,tmp2,eival,eivec integer nn dimension ginv(nn,nn),schol(nn,nn),tmp1(nn,nn),tmp2(nn,nn) dimension eival(nn),eivec(nn,nn) call vec_dcopy(ginv,schol,nn*nn) call cholesky(schol,nn) call vec_dinit(tmp1,nn*nn,0.0d0) call blas_dgemm $ ('n','n',nn,nn,nn,1.0d0,hi, $ nn,schol,nn,1.0d0,tmp1,nn) call vec_dinit(tmp2,nn*nn,0.0d0) call blas_dgemm $ ('t','n',nn,nn,nn,1.0d0,schol, $ nn,tmp1,nn,1.0d0,tmp2,nn) call eigen(tmp2,nn,nn,eival,eivec,tmp2,0) call matmult(schol,nn,eivec,nn,eivec,nn) return end integer function i_dfp(einmat,q,p,nn,upmat) c DFP-Update of Hessian, see c Heidrich,Kliesch,Quapp: Properties of Potential Energy Surfaces... c S.51 Eqn.(21)/(22) c c einmat Hessian to update (H_0) c q difference of gradient (g_1 - g_0) c p difference of point (x_1 - x_0) c upmat updated Hessian (H_1) real*8 einmat,q,p,upmat,qp,unit,b1,b2,b3, 1 vec_ddot,eps,sw integer nn dimension unit(nn,nn),b3(nn,nn), + einmat(nn,nn),upmat(nn,nn), + b1(nn,nn),b2(nn,nn),sw(nn,nn), + q(nn), p(nn) data eps/1.0d-30/ call vec_dinit(upmat,nn*nn,0.0d0) call mat_diag(unit,nn,1.0d0) qp=vec_ddot(q,p,nn) if(dabs(qp).lt.eps) call cerr(-31) call matmult(q,nn,p,nn,sw,1) call matadd(unit,(-1.0d0)/qp,sw,b1,nn,nn) call matmult(p,nn,q,nn,sw,1) call matadd(unit,(-1.0d0)/qp,sw,b2,nn,nn) call matmult(q,nn,q,nn,b3,1) call tri_mat(b1,nn,nn,einmat,nn,b2,nn,sw) call matadd(sw,1.0d0/qp,b3,upmat,nn,nn) i_dfp=0 2 continue return end real*8 function newton(hi,fi,step,gmat,nn) c calculates STEP = - HI^(-1) * FI and returns the steplength real*8 hi,fi,step,gmat,gl,vec_mdot integer i,nn dimension hi(nn,nn),fi(nn),step(nn),gmat(nn,nn),gl(nn,nn+1) call mat_dcopy(hi,gl,nn,nn) do 10,i=1,nn gl(i,nn+1)=-fi(i) 10 continue call linsolve(gl,step,nn) newton=dsqrt(vec_mdot(step,gmat,step,nn)) return end subroutine sig(ev,iv,nn) c count the number IV of positive, negative and zero eigenvalues EV real*8 ev integer iv,nn,i dimension ev(nn),iv(3) do i=1,3 iv(i)=0 enddo do i=1,nn if (ev(i).gt.(0.0d0)) iv(1)=iv(1)+1 if (ev(i).lt.(0.0d0)) iv(2)=iv(2)+1 if (ev(i).eq.(0.0d0)) iv(3)=iv(3)+1 enddo return end subroutine keytable(zkey,nr) c Table of KEYNAMES for the parameterfile character*40 zkey,ztab integer nr dimension ztab(30) ztab(1) ='Direction_of_Search' ztab(2) ='Dimension' ztab(3) ='Job_Name' ztab(4) ='File_of_Points' ztab(5) ='File_of_Gradient' ztab(6) ='File_of_Hessian' ztab(7) ='File_of_Results' ztab(8) ='File_of_B_matrix' ztab(9) ='Stopping_Criterion' ztab(10)='Predictor_Corrector_Criterion' ztab(11)='Length_of_Predictorsteps' ztab(12)='Min_Number_of_Steps' ztab(13)='Max_Number_of_Steps' ztab(14)='Flat_Space' ztab(15)='Free_First_Point' ztab(16)='Messages' ztab(17)='No_Update' ztab(18)='Const_Sign_Update' ztab(19)='Search_Even_Index' ztab(20)='Recorrect' ztab(21)='nul_switch' ztab(22)='nul' ztab(23)='nul' ztab(24)='nul' ztab(25)='nul' ztab(26)='nul' ztab(27)='nul' ztab(28)='nul' ztab(29)='File_of_Energy' ztab(30)='nul' zkey=ztab(nr) return end subroutine errtable(zkey,nr) c table for errormessages character*40 zkey,ztab integer nr,i dimension ztab(40) 10 continue c 1 0 0 0 0 ztab(1)='unspecific' ztab(2)='unspecific' ztab(3)='unspecific' ztab(4)='unspecific' ztab(5)='unspecific' ztab(6)='unspecific' ztab(7)='unspecific' ztab(8)='unspecific' ztab(9)='unspecific' ztab(10)='Direction_of_Search' ztab(11)='Number_of_Atoms' ztab(12)='Job_Name' ztab(13)='File_of_Points' ztab(14)='File_of_Gradient' ztab(15)='File_of_Hessian' ztab(16)='File_of_Results' ztab(17)='File_of_B_matrix' ztab(18)='Stopping_Criterion' ztab(19)='Predictor_Corrector_Criterion' ztab(20)='Length_of_Predictorsteps' ztab(21)='Min_Number_of_Steps' ztab(22)='Max_Number_of_Steps' ztab(23)='Gram-Schmidt algorithm' ztab(24)='unspecific' ztab(25)='unspecific' ztab(26)='unspecific' ztab(27)='unspecific' ztab(28)='unspecific' ztab(29)='File_of_Energy' ztab(30)='Div. by zero while inverting diag-matrix' ztab(31)='DFP update failed' ztab(32)='matrix inverting fails' ztab(33)='file not exists' ztab(34)='unspecific' ztab(35)='unspecific' ztab(36)='unspecific' ztab(37)='unspecific' ztab(38)='unspecific' ztab(39)='unspecific' ztab(40)='unspecific' if (nr.eq.0) then i=10 do while ((zkey.ne.ztab(i)).and.(i.le.22)) i=i+1 end do nr=-i return endif if ((nr.lt.0).and.(nr.ge.-40)) then zkey=ztab(-nr) return endif zkey='not known' return end subroutine loesung(matein,x,n) c intern linear equation routine real*8 matein,x integer i,j,n,k dimension matein(n,n+1),x(n) do 10,i=1,n 10 x(i)=matein(i,n+1) x(n)=x(n)/matein(n,n) do 20,i=1,n-1 k=n-i do 30,j=k+1,n 30 x(k)=x(k)-matein(k,j)*x(j) x(k)=x(k)/matein(k,k) 20 continue return end subroutine cholesky(g,ndim) c cholesky decomposition real*8 g,d integer ndim,i,j,k,l dimension g(ndim,ndim) i=ndim-1 do 1,j=1,i k=j+1 do 2,l=k,ndim g(j,l)=0.0d0 2 continue 1 continue do 3,i=1,ndim if (g(i,i).le.0.0d0) print*,'ERROR CHOL' g(i,i)=dsqrt(g(i,i)) j=i+1 if (j.gt.ndim) goto 3 d=g(i,i) do 4,k=j,ndim g(k,i)=g(k,i)/d 4 continue do 5,k=j,ndim do 6,l=k,ndim g(l,k)=g(l,k)-g(l,i)*g(k,i) 6 continue 5 continue 3 continue return end subroutine eigen (a,m,n,d,vec,e,iff) c eigenvalues and -vectors of A real*8 a,d,vec,e,eps,tol,h,g,s,f,b,p,r,c integer m,iff,n,i,j,ni,ii,l,k,j1 dimension a(m,m), d(m), vec(m,m) dimension e(m) eps = 0.278d-16 tol = 0.211758d-21 if (n.eq.1) then d(1) = a(1,1) vec(1,1) = 1.0d0 else do 30 i = 1 , n c householder's reduction c simulation of loop do 150 i=n,2,(-1) do 20 j = 1 , i vec(i,j) = a(i,j) 20 continue 30 continue do 120 ni = 2 , n ii = n + 2 - ni do 110 i = ii , ii l = i - 2 h = 0.0d0 g = vec(i,i-1) if (l.gt.0) then do 40 k = 1 , l h = h + vec(i,k)**2 40 continue s = h + g*g if (s.lt.tol) then h = 0.0d0 else if (h.gt.0) then l = l + 1 f = g g = dsqrt(s) if (f.gt.0) then g = -g end if h = s - f*g vec(i,i-1) = f - g f = 0.0d0 do 70 j = 1 , l vec(j,i) = vec(i,j)/h s = 0.0d0 do 50 k = 1 , j s = s + vec(j,k)*vec(i,k) 50 continue j1 = j + 1 if (j1.le.l) then do 60 k = j1 , l s = s + vec(k,j)*vec(i,k) 60 continue end if e(j) = s/h f = f + s*vec(j,i) 70 continue f = f/(h+h) do 80 j = 1 , l e(j) = e(j) - f*vec(i,j) 80 continue do 100 j = 1 , l f = vec(i,j) s = e(j) do 90 k = 1 , j vec(j,k) = vec(j,k) - f*e(k) - vec(i,k)*s 90 continue 100 continue end if end if c accumulation of transformation matrices d(i) = h e(i-1) = g 110 continue 120 continue d(1) = vec(1,1) vec(1,1) = 1.0d0 do 170 i = 2 , n l = i - 1 if (d(i).gt.0) then do 150 j = 1 , l s = 0.0d0 do 130 k = 1 , l s = s + vec(i,k)*vec(k,j) 130 continue do 140 k = 1 , l vec(k,j) = vec(k,j) - s*vec(k,i) 140 continue 150 continue end if d(i) = vec(i,i) vec(i,i) = 1.0d0 do 160 j = 1 , l c diagonalization of the tridiagonal matrix vec(i,j) = 0.0d0 vec(j,i) = 0.0d0 160 continue 170 continue b = 0.0d0 f = 0.0d0 e(n) = 0.0d0 do 250 l = 1 , n c test for splitting h = eps*(dabs(d(l))+dabs(e(l))) if (h.gt.b) b = h do 180 j = l , n c test for convergence if (dabs(e(j)).le.b) go to 190 180 continue 190 if (j.eq.l) then d(l) = d(l) + f else 200 p = (d(l+1)-d(l))*0.5d0/e(l) r = dsqrt(p*p+1.0d0) if (p.lt.0) then p = p - r else p = p + r end if h = d(l) - e(l)/p do 210 i = l , n c qr transformation d(i) = d(i) - h 210 continue f = f + h p = d(j) c simulation of loop do 330 i=j-1,l,(-1) c = 1.0d0 s = 0.0d0 j1 = j - 1 do 240 ni = l , j1 ii = l + j1 - ni do 230 i = ii , ii c protection against underflow of exponents g = c*e(i) h = c*p if (dabs(p).lt.dabs(e(i))) then c = p/e(i) r = dsqrt(c*c+1.0d0) e(i+1) = s*e(i)*r s = 1.0d0/r c = c/r else c = e(i)/p r = dsqrt(c*c+1.0d0) e(i+1) = s*p*r s = c/r c = 1.0d0/r end if p = c*d(i) - s*g d(i+1) = h + s*(c*g+s*d(i)) do 220 k = 1 , n h = vec(k,i+1) vec(k,i+1) = vec(k,i)*s + h*c vec(k,i) = vec(k,i)*c - h*s 220 continue 230 continue 240 continue e(l) = s*p c convergence d(l) = c*p if (dabs(e(l)).gt.b) go to 200 c ordering of eigenvalues d(l) = d(l) + f end if 250 continue if (iff.lt.1) then ni = n - 1 do 280 i = 1 , ni k = i p = d(i) j1 = i + 1 do 260 j = j1 , n if (d(j).lt.p) then k = j p = d(j) end if 260 continue if (k.ne.i) then d(k) = d(i) d(i) = p do 270 j = 1 , n p = vec(j,i) vec(j,i) = vec(j,k) vec(j,k) = p 270 continue end if c special treatment of case n = 1 280 continue end if end if return end subroutine linsolve(matein,x,n) c solves MATEIN * X = 0 real*8 matein,dre,x,w integer n dimension matein(n,n+1),dre(n,n+1),x(n),w(n) call householder(matein,dre,w,n,n+1) call loesung(dre,x,n) return end subroutine qr(gex,q,n) C QR-Zerlegung der Jacobimatrix von ge (gex) integer i,j,k,n real*8 s1,s2,s,gex,q,r,qq,rr dimension gex(n,n+1),q(n+1,n+1),r(n+1,n), * qq(n+1,n+1),rr(n+1,n) do 10,i=1,n do 20,j=1,n+1 r(j,i)=gex(i,j) 20 continue 10 continue do 30,i=1,n+1 do 40,j=1,n+1 q(i,j)=0.0d0 40 continue q(i,i)=1.0d0 30 continue do 50,i=1,n do 60,j=i+1,n+1 s1=r(i,i) s2=r(j,i) s=dsqrt(s1*s1+s2*s2) if (s.gt.0.0d0) then s1=s1/s s2=s2/s else s1=0.0d0 s2=0.0d0 endif call vec_dcopy(r,rr,(n+1)*n) do 70,k=1,n rr(i,k)=s1*r(i,k)+s2*r(j,k) rr(j,k)=s1*r(j,k)-s2*r(i,k) 70 continue call vec_dcopy(rr,r,(n+1)*n) call vec_dcopy(q,qq,(n+1)*(n+1)) do 80,k=1,n+1 qq(i,k)=s1*q(i,k)+s2*q(j,k) qq(j,k)=s1*q(j,k)-s2*q(i,k) 80 continue call vec_dcopy(qq,q,(n+1)*(n+1)) 60 continue 50 continue do 90,i=1,n+1 do 100,j=1,n+1 q(i,j)=qq(j,i) 100 continue 90 continue return end subroutine strcp(z1,z2,nn) c copies Character*NN Z1 to Z2 character*1 z1,z2 integer nn,i dimension z1(nn),z2(nn) do 10,i=1,nn 10 z2(i)=z1(i) return end subroutine recreate(ifile,zfile) c create an empty file integer iochk,ifile character*40 zfile open(ifile,file=zfile,iostat=iochk,status='unknown') close(ifile,iostat=iochk,status='delete') return end subroutine getvec(iofile,nkey,vec,nn) c reads vector VEC following the key with number NKEY from IOFILE logical oready integer nkey,i,iochk,nn,iofile character*40 zkey,zscr real*8 vec dimension vec(nn) oready=.false. rewind(iofile) call keytable(zscr,nkey) do while (.not.oready) read(iofile,*,end=10,iostat=iochk)zkey if (zkey.eq.zscr) then i=1 do while (i.le.nn) read(iofile,*,iostat=iochk)vec(i) call cerr(iochk) i=i+1 end do oready=.true. endif end do 10 continue if (.not.oready) call cerr(-10) return end subroutine gettang(curveq,gmat,tangent,nn) c This routine calculate the tangent on a curve c whereas the tangent is unique determined by: c c H*t=0 c |t|=1 c det[H,t]>0 c real*8 tangent,gmat,curveq,q,ss1,det,zero,dmnorm integer nn,i dimension tangent(nn),curveq(nn-1,nn),q(nn*nn),ss1(nn,nn) dimension gmat(nn,nn) data zero /0.0d0/ call qr(curveq,q,nn-1) call vec_dcopy(q(nn*(nn-1)+1),tangent(1),nn) call vec_dmult_constant c (tangent,nn,1/dmnorm(tangent,gmat,nn),tangent) call mat_dcopy(curveq,ss1,nn-1,nn) do 10,i=1,nn 10 ss1(nn,i)=tangent(i) if (det(ss1,nn).lt.zero) then call vec_dmult_constant(tangent,nn,-1.0d0,tangent) endif return end subroutine printv(io,val,zval,m) c writes a table ZVAL, VAL with M rows in IO integer m,io,i real*8 val character*20 zval dimension val(m),zval(m) write(io,20) write(io,30) write(io,20) do i=1,m write(io,10)zval(i),val(i) end do write(io,20) return 10 format (1x,A20,1x,f8.5) 20 format (1x,29('-')) 30 format (12x,'scalars') end real*8 function snorm(adj,grad,ginv,gmat,nn) c S norm real*8 adj,grad,ginv,dmnorm,ag,gmat integer nn dimension adj(nn,nn),grad(nn),ginv(nn,nn),ag(nn),gmat(nn,nn) call matmult(adj,nn,grad,1,ag,nn) snorm=dmnorm(ag,gmat,nn)/dmnorm(grad,ginv,nn) return end real*8 function bnorm(adj,grad,gmat,ag,nn) c Branin norm real*8 adj,grad,gmat,ag,dmnorm integer nn dimension adj(nn,nn),grad(nn),gmat(nn,nn),ag(nn) call matmult(adj,nn,grad,1,ag,nn) bnorm=dmnorm(ag,gmat,nn) call vec_dmult_constant(ag,nn,-1.0d0,ag) return end real*8 function ray(adj,grad,ginv,nn) c rayleigh coef. real*8 adj,grad,ginv,ag,vec_ddot,vec_mdot integer nn dimension adj(nn,nn),grad(nn),ginv(nn,nn),ag(nn) call matmult(adj,nn,grad,1,ag,nn) ray=vec_ddot(grad,ag,nn)/vec_mdot(grad,ginv,grad,nn) return end subroutine mat_adj(matin,ndim,matout) c adjoint matrix real*8 matin,matout,matscr,det integer ndim,i,j,k,l,ii,jj dimension matin(ndim,ndim),matout(ndim,ndim), + matscr(ndim-1,ndim-1) call vec_dinit(matscr,(ndim-1)*(ndim-1),0.0d0) do 10, i=1 , ndim do 20, j=1 , ndim do 30, ii=1 , ndim-1 do 40, jj=1 , ndim-1 if (ii.lt.i) then k=ii else k=ii+1 endif if (jj.lt.j) then l=jj else l=jj+1 endif matscr(ii,jj)=matin(k,l) c if(dabs(matscr(ii,jj)).lt.1.0d-30) then c matscr(ii,jj)=0.0d0 c endif 40 continue 30 continue matout(i,j)=det(matscr,ndim-1) if(dabs(matout(i,j)).lt.1.0d-30) goto 20 matout(i,j)= (-1.0d0)**(i+j)*matout(i,j) 20 continue 10 continue return end real*8 function cnorm(vec,nn) c cartesian NORM of VEC real*8 vec,vec_ddot integer nn dimension vec(nn) cnorm=dsqrt(vec_ddot(vec,vec,nn)) return end real*8 function dmnorm(vec,gm,nn) c metric norm of VEC with resp. to metric GM real*8 vec,vec_mdot,gm integer nn dimension vec(nn),gm(nn,nn) dmnorm=dsqrt(vec_mdot(vec,gm,vec,nn)) return end subroutine help(io,index) integer io,index c 5=stdin, 6=stdout, 0=stderror character*120 zpara,zhelp,zver,zret character*20 zrgf dimension zpara(5),zhelp(3),zver(5),zret(5) zpara(1)= c'Format of param.rgf:\n (Components of searchdir separated '// c'by space or newline)\n ¶m\n job_id=(number)' zpara(2)= c'nn=(dimension)\n z_FoP=(filename)\n z_FoG=(filename)\n z_FoH'// c'=(filename)\n z_FoB=(filename)\n z_res=(filename)' zpara(3)= c 'o_flat={T|F}\n o_free={T|F}\n o_msg={T|F}\n o_no_up={T|F}\n'// c ' o_const_sgn={T|F}\n o_even={T|F}\n o_recorr={T|F}' zpara(4)= c 'o_nn={T|F}\n crit_stop=(Fortran double)\n crit_pc=(Fortran'// c ' double)\n pstepl=(Fortran double)' zpara(5)= c 'itmin=(stepnumber)\n itmax=(stepnumber)\n /' zhelp(1)='Usage: (programname) [options]\n Options:\n '// c '-ver\t\tprint version information\n '// c '-para\t\tprint the format of »param.rgf«' zhelp(2)='-ret\t\tprint a desciption of return values\n '// c '-exact\t\tforces "exact HESSIAN"\n '// c '-update\tforces "update HESSIAN"' zhelp(3)='-reinit\treread PSTEPL,CRIT_STOP,CRIT_PC,ITMIN,ITMAX'// & '\n\t\tfrom the parameterfile' zver(1)='RGF (»Reduced Gradient Following«) is a realisation'// c '\n of the "following a reduced gradient" - method\n' zver(2)= c 'Reference:\tW.Quapp, M.Hirsch, D.Heidrich:\n '// c ' \t\tTheor Chem Acc (1998) 100:285-299' zver(3)= c '\t\tW.Quapp, M.Hirsch, O.Imig, D.Heidrich:\n \t\tJ Comput '// c 'Chem 19 (1998) 1087-1100\n' zver(4)= c 'First release: Apr 19, 2000\n Revision 1: Nov 11, 2000\n'// c ' Revision 2: Jul 05, 2001\n' zver(5)='Contact to the authors and programmers:\n '// c 'E-mail:\thirsch@mathematik.uni-leipzig.de\n '// c '\t\tquapp@mathematik.uni-leipzig.de' zret(1)= c 'Legend of return values\n 0\tafter printing any help message'// c ', like this\n 1\tparam.rgf not exists' zret(2)= c '\tOR error while reading parameters\n 4\tPREDICTOR in step '// c 'executed;\n\tnext hessian: UPDATE' zret(3)= c '5\tPREDICTOR in step executed;\n\tnext hessian: EXACT\n '// c '6\tCORRECTOR in step executed;\n\tnext hessian: UPDATE' zret(4)= c '7\tCORRECTOR in step executed;\n\tnext hessian: EXACT \n'// c '10\tmax. Stepnumber reached' zret(5)= c '\tOR stopping criterion satisfied' call getarg (0,zrgf) if (index.eq.0) write(io,*) zhelp if (index.eq.1) write(io,*) zpara if (index.eq.2) write(io,*) zver if (index.eq.3) write(io,*) zret return end character*40 function z_iostat(iostat) integer iostat select case (iostat) case(0) z_iostat='successfull io' case(100) z_iostat='error in format' case(101) z_iostat='illegal unit number' case(102) z_iostat='formatted io not allowed' case(103) z_iostat='unformatted io not allowed' case(104) z_iostat='direct io not allowed' case(105) z_iostat='sequential io not allowed' case(106) z_iostat='cant backspace file' case(107) z_iostat='null file name' case(108) z_iostat='cant stat file' case(109) z_iostat='unit not connected' case(110) z_iostat='off end of record' case(111) z_iostat='truncation failed in endfile' case(112) z_iostat='incomprehensible list input' case(113) z_iostat='out of free space' case(114) z_iostat='unit not connected' case(115) z_iostat='read unexpected character' case(116) z_iostat='bad logical input field' case(117) z_iostat='bad variable type' case(118) z_iostat='bad namelist name' case(119) z_iostat='variable not in namelist' case(120) z_iostat='no end record' case(121) z_iostat='variable count incorrect' case(122) z_iostat='subscript for scalar variable' case(123) z_iostat='invalid array section' case(124) z_iostat='substring out of bounds' case(125) z_iostat='subscript out of bounds' case(126) z_iostat='cant read file' case(127) z_iostat='cant write file' case(128) z_iostat='\"new\" file exists' case(129) z_iostat='cant append to file' case(130) z_iostat='non-positive record number' case(131) z_iostat='I/O started while already doing I/O' case default z_iostat='unknown IO status' end select return end subroutine handleargs(n,ostat) integer n,i,istat,ind logical ostat character*10 zarg character*40 zall dimension ostat(8) data zall /'-ver-para-ret-exact-update-reinit'/ do i=3,8 ostat(i)=.false. enddo do i=1,n call getarg(i,zarg) istat=len_trim(zarg) ind=index(zall,zarg(1:istat)) select case (ind) case(1) call help(6,2) call exit(0) case(5) call help(6,1) call exit(0) case(10) call help(6,3) call exit(0) case(14) ostat(3)=.true. case(20) ostat(4)=.true. case(27) ostat(5)=.true. case default call help(6,0) call exit(0) end select enddo return end c Modification since c Thu Jun 19 13:01:45 CEST 2001 subroutine func_prop(p1,p2,p3,iprop,zprop) real*8 p1,p2,p3 integer iprop,ibset character*20 zprop character*6 zmon,zext,zchg iprop=0 zmon='decr ' zext='mono ' zchg='const ' if (p3.gt.p2) then iprop=ibset(iprop,0) zmon='incr ' endif if ((p2-p1)*(p2-p3).ge.0.0d0) then iprop=ibset(iprop,1) zext='extr ' endif if (p2*p3.lt.0.0d0) then iprop=ibset(iprop,2) zchg='change' endif zprop=zmon//','//zext//','//zchg return end subroutine reread(zfile) integer job_id,nn,itmin,itmax,is character*40 z_FoP,z_FoG,z_FoH,z_FoB,z_FoR,zfile,z_FoE real*8 crit_stop,crit_pc,pstepl logical o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even, c o_recorr,o_nn,ostat dimension ostat(8) common /para/crit_stop,crit_pc,pstepl,itmin,itmax,ostat namelist /param/job_id,nn,z_FoP,z_FoG,z_FoH,z_FoB,z_FoR, c o_flat,o_free,o_msg,o_no_up,o_const_sgn,o_even,o_recorr,o_nn, c crit_stop,crit_pc,pstepl,itmin,itmax,z_FoE inquire(file=zfile,exist=ostat(1),opened=ostat(2),number=is) if (.not.ostat(1)) return if (ostat(2)) close(is) open(33,file=zfile,status='old') read(33,param,iostat=is) if (is.ne.0.and.is.ne.119) then write(*,*)'wrong format of parameter file' call help(6,1) call exit(1) endif close(33) return end C SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL c c see http://csep1.phy.ornl.gov/CSEP/MO/NODE11B.html or c T Schlick, in Rev.Comput.Chem.,Vol.III (1992) p.1-72 c c SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL SLSLSLSLSLSL c subroutine search_line(o_even,gra,yy,en,fL2,po,xx,lam1,Stepl, c SeDir,nn) implicit real*8 (a-h,p-z), integer(i-n) real*8 lam1,gra,en,po,Stepl,SeDir,xx,yy,aL,bL,fL,fL2,gL,gL2, c salt,vec_ddot,lambda integer nn logical o_even dimension gra(nn),po(nn),SeDir(nn),xx(nn),yy(nn) c c lam1 first EigenValue c xx line search point c yy gradient there c fL2 energy in test point salt=Stepl lambda=lam1 fL= en gL= vec_ddot(gra,SeDir,nn) call vec_dmult_add(po,SeDir,nn,Stepl,xx) gL2= vec_ddot(yy,SeDir,nn) 2 continue if(o_even) then If( lambda .gt. 0.0d0 ) then aL= (-2.0d0*(fL2-fL)+(gL2+ gL)*Stepl)/Stepl**3 bL= ( 3.0d0*(fL2-fL)-(gL2+2.0d0*gL)*Stepl)/Stepl**2 If((aL.ne. 0.0d0).and.(bL**2-3.0d0*aL*gL.gt.0.0d0) ) then Stepl= -(bL +Dsqrt(bL**2-3.0d0*aL*gL))/(3.0d0*al) else Stepl= -gL/(2.0d0*bl) endif else If(fL2.gt. fL) then lambda=1.0d0 goto 2 else If(fL2.lt. fL) Stepl=Stepl*1.111111d0 endif endif c else If(lambda.lt. 0.0d0) then aL=(-2.0d0*(fL2-fL)+(gL2+ gL)*Stepl)/Stepl**3 bL=( 3.0d0*(fL2-fL)-(gL2+2.0d0*gL)*Stepl)/Stepl**2 If((aL.ne. 0.0d0) .and. (bL**2-3.0d0*aL*gL.gt.0.0d0)) then Stepl= -(bL +Dsqrt(bL**2-3.0d0*aL*gL))/(3.0d0*al) else Stepl= -gL/(2.0d0*bl) endif else cccccc If(fL2.gt.fL) Stepl=Stepl c uphill in minimum bowl: LineSearch is still open !! cccccc endif endif c if(Stepl.gt.2.222d0*salt) Stepl=2.222d0*salt if(Stepl.lt.0.0d0) Stepl=salt/2.0d0 return end