subroutine dfp(inp,iout,nfile,nsave,maxit,tolx,tolg,tolrec,qua, 1 rest,iest,lint,n,*) c c this subroutine performs minimization of a certain objective function, c f, by the davidon-fletcher-powell method. diagonal and upper-diagonal c elements of the matrix v are stored in the vector v. linesearch is c included in the subroutine dirmin. the variables are converted into c their working form (vector x), and vice-versa by the subroutine convrt c calculation terminates normally if both tests on delta(x) and delta(g) c are satisfied. maximum number of cycles is maxit. c c ************** written by a.liwo, institute of chemistry c ************** university of gdansk c ************** ul.sobieskiego 18 c ************** 80952 gdansk, poland c c ************** first r-32 version (fortran iv): august 1986 c ************** ibm pc adaptation (fortran 77): may 1988 c ************** unix version (fortran 77): november 1990 c implicit real*8 (a-h,o-z) include 'DIMENSIONS.ZSCOPT' parameter (maxdim=max_paropt) parameter (maxvdim=maxdim*(maxdim+1)/2) logical fail,init,qua,rest,lint common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim), 1 vy(maxdim),g(maxdim) common /dirit/ nprint print *,"DFP n",n call convrt(n,x,.true.) nvdim=n*(n+1)/2 isave=nsave kkk=1 if (.not.rest) goto 20 print *,'reading restart information from the restart file' rewind nfile nstep=nvdim/100 nstep1=nstep if (nstep*100.lt.nvdim) nstep1=nstep+1 read(nfile) kkk,f,dftau0 do 222 i=1,n 222 read (nfile) g(i),y(i),x(i),s(i) do 111 i=1,nstep1 ns=(i-1)*100+1 ne=i*100 if (ne.gt.nvdim) ne=nvdim 111 read (nfile) (v(j),j=ns,ne) endfile nfile print *,'restart information has been read in' if (kkk.gt.1) kk1=kkk goto 21 20 f=fun(n,x,nfun) call grad(n,iest,g,x) do 201 i=1,n 201 y(i)=g(i) 21 n1=n+1 nfun=0 do 2 it=1,maxit if (kkk.gt.1) goto 17 kk1=2 ind=0 do 3 i=1,n i1=i+1 ind=ind+1 v(ind)=1. if (i1.gt.n) goto 3 do 4 j=i1,n ind=ind+1 4 v(ind)=0 3 continue do 5 i=1,n 5 d(i)=-g(i) dftau0=dirder(n,g,d) call dirmin(iout,iest,nfun,n,f,dftau0,d,y,x,s,fail,init,qua) if (nprint.gt.2) write (iout,1040) (i,x(i),i=1,n) if (fail) return 1 if (.not.init) goto 1 if (nprint.gt.-1) write (iout,1100) it,nfun if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n) 1100 format (//1x,' * * * * * * * * no point remarkably lower in ene', 1 'rgy can be found.'/1x,'calculation terminated in cycle',i3, 2 ' after',i6,' function & gradient evaluations') call convrt(n,x,.false.) * nprint=2 call csave(lint) return 1 if (n.eq.1) return xnorm=0.0d0 gnorm=0.0d0 xmax=0.0d0 gmax=0.0d0 do 163 i=1,n dxi=dabs(s(i)) gi=dabs(y(i)) if (dxi.gt.xmax) xmax=dxi if (gi.gt.gmax) gmax=gi xnorm=xnorm+dxi*dxi 163 gnorm=gnorm+gi*gi xnorm=dsqrt(xnorm/n) gnorm=dsqrt(gnorm/n) kk=1 if(nprint.gt.0)write(iout,1000)it,kk,xmax,xnorm,gmax,gnorm,f,nfun if (nprint.gt.2) write (iout,1070) (i,y(i),i=1,n) 1070 format (/5x,' * * * * gradients'/100(1x,i5,5x,1pe12.4/)) if (nfun.lt.isave) goto 23 call convrt (n,x,.false.) call psave(n,nvdim,nfun,nfile,kk,it,f,dftau0) call csave(lint) isave=isave+nsave 23 if (xmax.gt.tolx.or.gmax.gt.tolg) goto 17 if (nprint.gt.-1) write (iout,1010) it,nfun if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n) return 17 do 6 kk=kk1,n1 do 7 k=1,n gk=y(k) y(k)=gk-g(k) 7 g(k)=gk sy=0. do 8 k=1,n 8 sy=sy+s(k)*y(k) if (dabs(sy).gt.tolrec) goto 162 if (nprint.gt.-1) write (iout,1030) 1030 format (//' * * * delta(x) is perpendicular to delta(g). ', 1 'the program will try gradient direction * * *') goto 2 162 do 9 i=1,n vyi=0. do 10 j=1,i ind=(2*n-j)*(j-1)/2+i 10 vyi=vyi+v(ind)*y(j) i1=i+1 if (i1.gt.n) goto 9 do 11 j=i1,n ind=(2*n-i)*(i-1)/2+j 11 vyi=vyi+v(ind)*y(j) 9 vy(i)=vyi yvy=0. do 12 i=1,n 12 yvy=yvy+y(i)*vy(i) ind=0 do 13 i=1,n si=s(i)/sy vyi=vy(i)/yvy do 13 j=i,n ind=ind+1 v(ind)=v(ind)+si*s(j)-vyi*vy(j) 13 continue if (nprint.gt.4) call matout(nvdim,n) do 18 i=1,n if (v((2*n-i)*(i-1)/2+i).gt.0.) goto 18 if (nprint.gt.-1) write (iout,1080) i goto 2 18 continue 1080 format (/1x,i3,'-th diagonal element of the matrix v is negative.' 1 ,' restart with gradient direction follows.') do 14 i=1,n di=0. do 15 j=1,i ind=(2*n-j)*(j-1)/2+i 15 di=di-v(ind)*g(j) i1=i+1 if (i1.gt.n) goto 14 do 16 j=i1,n ind=(2*n-i)*(i-1)/2+j 16 di=di-v(ind)*g(j) 14 d(i)=di dftau0=dirder(n,g,d) do 202 i=1,n 202 y(i)=g(i) call dirmin(iout,iest,nfun,n,f,dftau0,d,y,x,s,fail,init,qua) if (fail) return 1 if (.not.init) goto 19 if (nprint.gt.-1) write (iout,1110) 1110 format (/1x,'energy from dirmin differs but slightly from the' 1 ,' initial value. restart with gradient direction follows.') goto 2 19 xnorm=0.0d0 gnorm=0.0d0 xmax=0.0d0 gmax=0.0d0 do 161 i=1,n dxi=dabs(s(i)) gi=dabs(y(i)) if (xmax.lt.dxi) xmax=dxi if (gmax.lt.gi) gmax=gi xnorm=xnorm+dxi*dxi 161 gnorm=gnorm+gi*gi xnorm=dsqrt(xnorm/n) gnorm=dsqrt(gnorm/n) if (nprint.gt.0) write (iout,1000) it,kk,xmax,xnorm,gmax,gnorm,f, 1 nfun if (nprint.gt.2) write (iout,1040) (i,x(i),i=1,n) 1040 format (/1x,' * * * * parameters:'/100(1x,i5,5x,1pe14.6/)) if (nprint.gt.2) write (iout,1070) (i,y(i),i=1,n) 1000 format (/1x,'iteration:',i4,' direction:',i4,' max. element of ', 1 'delta(x): ',1pe14.6,' rms of delta(x): ',1pe14.6/' maximum ', 2 'element of gradient: ',1pe14.6,' rms of gradient: ',1pe14.6, 3 ' energy: ',1pe14.6,' fun. eval:',i6) c if (nfun.lt.isave) goto 22 call convrt(n,x,.false.) call psave(n,nvdim,nfun,nfile,kk,it,f,dftau0) call csave(lint) isave=isave+nsave 22 if (xmax.gt.tolx.or.gmax.gt.tolg) goto 6 if (nprint.gt.-1) write (iout,1010) it,nfun if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n) if (nprint.gt.-1) write (iout,1111) 1111 format (//11x,'matrix v') if (nprint.gt.-1) call matout(nvdim,n) 1010 format (//1x,'* * * * * * convergence has been achieved in ', 1 'cycle',i4,' after',i6,' function & gradient evaluations') call convrt(n,x,.false.) * nprint=2 call csave(lint) return 6 continue 2 kkk=1 c if this point is reached, convergence hasn't been achieved. if (nprint.gt.-1) write (iout,1020) maxit,nfun if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n) if (nprint.gt.-1) write (iout,1111) if (nprint.gt.-1) call matout(nvdim,n) call convrt(n,x,.false.) * nprint=2 call csave(lint) 1020 format (//1x,' * * * * * * * convergence has not been achieved', 1 ' * * * * * * * '/1x,'number of cycles done:',i3,' number of ', 2 'function & gradient evaluations:',i6) return end ******************************************************************************* subroutine grad(n,iest,g,x) implicit real*8 (a-h,o-z) include 'DIMENSIONS.ZSCOPT' dimension g(max_paropt),x(max_paropt) integer uiparm,nf,iest external fdum double precision obg(max_paropt) call targetgrad(n,x,nf,g,uiparm,obg,fdum) return end ******************************************************************************* subroutine dirmin(iout,iest,nfun,n,f,dftau0,d,g,x,s,fail,init,qua) implicit real*8 (a-h,o-z) include 'DIMENSIONS.ZSCOPT' parameter (maxdim=max_paropt) logical fail,init,qua,weiter common /dirit/ nprint,maxexp,maxlin,maxd,expand,tolqu,dermin,amov, 1 xmin dimension x(n),s(n),d(n),g(n),z(maxdim),h(maxdim) * write (2,*) 'dirmin has been called' fail=.false. init=.true. weiter=.false. amx=dabs(d(1)) do 22 i=2,n pom=dabs(d(i)) if (pom.gt.amx) amx=pom 22 continue tau=dmin1(1.0d0,amov/amx) do 1 i=1,n s(i)=x(i) h(i)=g(i) 1 z(i)=x(i) * write (2,'(2hh:/(i3,1pe15.5))' ) (i,h(i),i=1,n) if (nprint.gt.1) write (iout,1020) dftau0 do 2 ilin=1,maxlin * write (2,'(2hh:/(i3,1pe15.5))' ) (i,h(i),i=1,n) if (weiter) goto 10 do 3 id=1,maxd do 4 iexp=1,maxexp do 5 i=1,n 5 x(i)=s(i)+tau*d(i) f1=fun(n,x,nfun) call grad(n,iest,g,x) dftau=dirder(n,g,d) if (nprint.gt.1) write (iout,1000) ilin,tau,iexp,f1,f,dftau 1000 format (1x,'dirmin iter: ',i4,' tau=',1pe12.4, 1 ' contr. step:',i4/' actual f value: ',1pe14.6,' old value: ', 2 1pe14.6,' derivative: ',1pe14.6) if (nprint.gt.4) write (iout,1080) (i,x(i),i=1,n) if (nprint.gt.5) write (iout,1070) (i,g(i),i=1,n) if (dabs(f1-f).gt.tolqu.and.amx*tau.ge.xmin) goto 21 if (f1.lt.f) goto 14 do 24 i=1,n x(i)=s(i) g(i)=h(i) 24 s(i)=x(i)-z(i) return 21 if (f1.lt.f.or.dftau.gt.0.) goto 6 4 tau=tau/expand if (nprint.gt.-1) write(iout,1010) f do 25 i=1,n g(i)=h(i) 25 x(i)=s(i) if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n) fail=.true. return 1010 format (//' * * * * * dirmin failed. least function value is: ', 1 1pe14.6,' * * * * *') 1080 format (/' * * * * parameters:'/100(1x,i5,5x,1pe14.6/)) 1070 format (/' * * * * gradients:'/100(1x,i5,5x,1pe14.6/)) 1020 format (/1x,'initial value of directional derivative: ',1pe14.6) 6 if (dftau.gt.0.) goto 10 init=.false. if (dabs(dftau/amx).le.dermin) goto 14 dftau0=dftau do 9 i=1,n h(i)=g(i) 9 s(i)=x(i) tau=tau*expand f=f1 3 continue if (nprint.gt.-1) write (iout,1030) 1030 format (//1x,' * * * * * error in dirmin. positive directional ', 1 'derivative cannot be found * * * * *') if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n) fail=.true. return 10 term=f1-f-dftau0*tau if (qua.or.dabs(dftau0/amx).lt.0.01.and.dabs(dftau/amx).lt.0.01) 1 goto 30 difd=tau*(dftau-dftau0) tt=tau*tau ttt=tt*tau aa=(tau*difd-2.*term)/ttt if (aa.eq.0.) goto 30 bb=(3.*term-difd)/tt delta=4.*bb*bb-12.*aa*dftau0 if (delta.lt.0.) goto 30 delta=dsqrt(delta) tau1=(-2.*bb-delta)/(6.*aa) tau2=(-2.*bb+delta)/(6.*aa) if (nprint.gt.1) write (iout,1113) ilin,tau1,tau2 1113 format (/1x,'dirmin iter.:',i3,' cubic interpolation', 1 ' tau1=',1pe12.4,' tau2= ',1pe12.4) if (tau1.gt.0.0d0.and.tau1.le.tau) goto 16 if (tau2.le.0.0d0.or.tau2.gt.tau) goto 30 tau1=tau2 goto 16 30 tau1=-dftau0*tau*tau/(2.*term) if (nprint.gt.1) write (iout,1040) ilin,tau1 if (tau1.lt.tau.and.tau1.gt.0.) goto 16 if (dabs(amx*tau1).lt.xmin) goto 145 weiter=.false. tau=tau/expand if (f1.gt.f) goto 2 do 17 i=1,n s(i)=x(i) h(i)=g(i) 17 d(i)=-d(i) f=f1 dftau0=-dftau goto 2 1040 format (/1x,'dirmin iteration:',i4,' tau after quadr. interpo', 1 'lation: ',1pe12.4) 16 do 11 i=1,n 11 x(i)=s(i)+tau1*d(i) if (nprint.gt.4) write (iout,1080) (i,x(i),i=1,n) if (nprint.gt.5) write (iout,1070) (i,g(i),i=1,n) f2=fun(n,x,nfun) call grad(n,iest,g,x) dftau1=dirder(n,g,d) dftn=dabs(dftau1/amx) if (nprint.gt.1) write (iout,1112) f2,dftau1 1112 format (1x,'function value: ',1pe14.6,' derivative: ',1pe14.6) if (f2.lt.f.or.dftau1.gt.0.) goto 20 if (dabs(f2-f).gt.tolqu.and.dabs(amx*tau1).gt.xmin.and.dftn.ge. 1 dermin) goto 27 do 28 i=1,n x(i)=s(i) g(i)=h(i) 28 s(i)=x(i)-z(i) goto 18 27 weiter=.false. tau=tau1/expand goto 2 20 weiter=.true. if (f2.lt.f) init=.false. if (f2.lt.f) goto 26 if (dabs(f2-f).gt.tolqu.and.dabs(amx*tau1).ge.xmin.and.dftn.ge. 1 dermin) goto 31 do 32 i=1,n x(i)=s(i) s(i)=x(i)-z(i) 32 g(i)=h(i) goto 18 31 f1=f2 dftau=dftau1 tau=tau1 goto 2 26 taudf=tau1*dftau1 if (nprint.gt.1) write (iout,1111) taudf 1111 format (/1x,'difference: ',1pe10.2) if (dabs(taudf).lt.tolqu.or.dftn.lt.dermin) goto 144 if (dftau1.lt.0.) goto 13 do 12 i=1,n d(i)=-d(i) h(i)=g(i) 12 s(i)=x(i) dftau=-dftau0 dftau0=-dftau1 f1=f f=f2 tau=tau1 goto 2 13 continue do 19 i=1,n h(i)=g(i) 19 s(i)=x(i) dftau0=dftau1 f=f2 tau=tau-tau1 2 continue if (nprint.gt.-1) write (iout,1050) 1050 format (//1x,' * * * * nonconvergence in dirmin, calculation ', 1 'terminated. * * * *') fail=.true. if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n) return 144 f=f2 goto 145 14 f=f1 145 do 15 i=1,n 15 s(i)=x(i)-z(i) 18 if (nprint.gt.1) write (iout,1060) ilin 1060 format (/1x,'dirmin converged in',i3,' iterations') return end ******************************************************************************* double precision function dirder(n,g,d) implicit real*8 (a-h,o-z) dimension d(n),g(n) c calculates the directional derivative. d is the direction vector, c g - gradient vector yy=0.0d0 do 1 i=1,n yy=yy+d(i)*g(i) 1 continue dirder=yy return end ******************************************************************************* double precision function fun(n,x,nfun) implicit real*8 (a-h,o-z) include 'DIMENSIONS.ZSCOPT' parameter (maxdim=max_paropt) dimension x(max_paropt) external fdum integer uiparm double precision f,obf print *,"FUN: n=",n nfun=nfun+1 call targetfun(n,x,nfun,f,uiparm,obf,fdum) fun=f return end ******************************************************************************* subroutine convrt(n,x,dir) implicit real*8 (a-h,o-z) dimension x(n) return end ******************************************************************************* subroutine dfpopt(inp,iout,ifirst,n,xopt) implicit real*8 (a-h,o-z) c c this subroutine reads in the iteration parameters, calls c the minimization procedures (davejr or dfp or both of them one by c one) and writes the final coordinates to file ares.xyz. c include 'DIMENSIONS.ZSCOPT' parameter (maxdim=max_paropt) parameter (maxvdim=maxdim*(maxdim+1)/2) common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim), 1 vy(maxdim),g(maxdim) logical conv,qua,rest,lint common /dirit/ nprint,maxexp,maxlin,maxd,expand,tolqu,dermin,amov, 1 xmin c integer option,odfp,onr dimension xopt(maxdim) data odfp /4h dfp/,onr /4h nr/, ncu /4h cub/,nrest /4hrest/ print *,'ce2 has been called' print *,"DFPOPT: n",n c c if (ifirst.eq.0) goto 1 c read (inp,1001) option,icub,irest c rest=.false. c if (irest.eq.nrest) rest=.true. c qua=.true. c if (icub.eq.ncu) qua=.false. c read (inp,1010) maxit,itolx,itolg,maxlin,maxexp,maxd,itolq,itolr, c 1 iderm,nprint,ixmin,nsave,nfile,iest c read (inp,1030) expand,amov c if (nsave.eq.0) nsave=10000 c nfile=4 c if (itolx.eq.0) itolx=3 c if (itolg.eq.0) itolg=1 c if (itolr.eq.0) itolr=6 c if (iderm.eq.0) iderm=7 c if (itolq.eq.0) itolq=4 c if (ixmin.eq.0) ixmin=5 c if (maxit.eq.0) maxit=5 c if (maxlin.eq.0) maxlin=10 c if (maxd.eq.0) maxd=10 c if (maxexp.eq.0) maxexp=10 c if (expand.eq.0.) expand=2.0d0 c if (amov.eq.0.) amov=.25d0 c tolx=10.0d0**(-itolx) c tolg=10.0d0**(-itolg) c xmin=10.0d0**(-ixmin) c tolqu=10.0d0**(-itolq) c dermin=10.0d0**(-iderm) c tolrec=10.0d0**(-itolr) c if (option.eq.odfp) goto 4 do i=1,n x(i)=xopt(i) enddo nprint=2 nsave=100000 nfile=89 maxit=100 maxlin=10 maxd=10 maxexp=10 expand=2.0d0 amov=.25d0 tolx=1.0d-3 tolg=1.0d-1 xmin=1.0d-5 tolqu=1.0d-3 dermin=1.0d-7 tolrec=1.0d-6 iest=0 qua=.false. rest=.false. 4 write (iout,2000) maxit,tolx,tolg,maxlin,maxexp,maxd,tolqu write (iout,2002) expand,amov,xmin,dermin,tolrec if (iest.gt.0) write (iout,2060) write (iout,2040) nsave,nfile if (rest) write (iout,2050) nfile if (qua) write (iout,2020) if (.not.qua) write (iout,2030) 1 call dfp(inp,iout,nfile,nsave,maxit,tolx,tolg,tolrec,qua,rest, 1 iest,lint,n,*2) do 300 i=1,n 300 xopt(i)=x(i) return 2 write (iout,2003) do 301 i=1,n 301 xopt(i)=x(i) 1001 format (3a4) 1030 format (8f10.5) 1010 format (14i5) 2000 format (//'welcome to dfp routine'//'the following ', 1 'options are used:'// 2 'maximum number of cycles:',i4/'tolerance on coordinates:', 3 1pe10.2/'tolerance on gradient:',1pe10.2/ 4 'maximum number of iterations in linesearch:',i4/ 5 'maximum number of unsuccesful function evaluations:',i4/ 6 'maximum number of steps in positive derivative search:',i4/ 7 'tolerance in linesearch:',1pe10.2) 2002 format ('expansion coefficient in linesearch:',1pe10.2/ 9 'maximum movement of atoms:',1pe10.2/ a 'minimum change of coordinates considered significant:',1pe10.2/ c 'minimum value of directional derivative considered different ', d 'from zero:',1pe10.2/ e 'lower bound of the scalar product (delta(x),delta(g)):', f 1pe10.2) 2003 format (/'not quite succesful... perhaps you will have a ' 1 'better time'/) c 2060 format (/'estimated partials used.') 2020 format (/'quadratic interpolation will be used in line', 1 ' minimization') 2030 format (/'cubic or quadratic interpolation will be used ', 1 'in line minimization') 2040 format (/'coordinates, gradients, and matrix v will be saved', 1 ' appr. every',i4,' energy and gradient evaluations on file',i3) 2050 format (//'restart to be done using information saved ', 1 'on file',i3/) return end ******************************************************************************* subroutine csave(lint) logical lint return end ******************************************************************************* subroutine psave(n,nvdim,nfun,nfile,kk,it,f,dftau0) implicit real*8 (a-h,o-z) include 'DIMENSIONS.ZSCOPT' parameter (maxdim=max_paropt) parameter (maxvdim=maxdim*(maxdim+1)/2) logical fail,init,qua,rest,lint common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim), 1 vy(maxdim),g(maxdim) common /dirit/ nprint common /restart/ restfile character*16 restfile common inp,iout,ipn,iprint rewind nfile nstep=nvdim/100 nstep1=nstep if (nstep*100.lt.nvdim) nstep1=nstep+1 write (nfile) kk,f,dftau0 do 222 i=1,n 222 write (nfile) g(i),y(i),x(i),s(i) do 111 i=1,nstep1 ns=(i-1)*100+1 ne=i*100 if (ne.gt.nvdim) ne=nvdim 111 write (nfile) (v(j),j=ns,ne) endfile nfile write (*,1010) it,nfun,f,restfile 1010 format (/1x,'cycle: ',i2,' energy evaluations: ',i5,' value:', 1 f20.10,'.'/' information for restart saved on file ',a16) return end ******************************************************************************* subroutine matout(ndim,n) implicit real*8 (a-h,o-z) include 'DIMENSIONS.ZSCOPT' parameter (maxdim=max_paropt) parameter (maxvdim=maxdim*(maxdim+1)/2) logical fail,init,qua,rest,lint common/bigmat/ a(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim), 1 vy(maxdim),g(maxdim) character*6 line6 / '------' / character*12 line12 / '------------' / dimension b(8) data iout /2/ do 1 i=1,n,6 nlim=min0(i+5,n) write (iout,1000) (k,k=i,nlim) write (iout,1020) line6,(line12,k=i,nlim) 1000 format (/7x,6(i7,5x)) 1020 format (a6,6a12) do 2 j=i,n jlim=min0(j,nlim) jlim1=jlim-i+1 do 3 k=i,jlim 3 b(k-i+1)=a(n*(k-1)-k*(k-1)/2+j) write (iout,1010) j,(b(k),k=1,jlim1) 2 continue 1 continue 1010 format (i3,3x,6(1pe11.4,1x)) return end