1 subroutine dfp(inp,iout,nfile,nsave,maxit,tolx,tolg,tolrec,qua,
4 c this subroutine performs minimization of a certain objective function,
5 c f, by the davidon-fletcher-powell method. diagonal and upper-diagonal
6 c elements of the matrix v are stored in the vector v. linesearch is
7 c included in the subroutine dirmin. the variables are converted into
8 c their working form (vector x), and vice-versa by the subroutine convrt
9 c calculation terminates normally if both tests on delta(x) and delta(g)
10 c are satisfied. maximum number of cycles is maxit.
12 c ************** written by a.liwo, institute of chemistry
13 c ************** university of gdansk
14 c ************** ul.sobieskiego 18
15 c ************** 80952 gdansk, poland
17 c ************** first r-32 version (fortran iv): august 1986
18 c ************** ibm pc adaptation (fortran 77): may 1988
19 c ************** unix version (fortran 77): november 1990
21 implicit real*8 (a-h,o-z)
22 include 'DIMENSIONS.ZSCOPT'
23 parameter (maxdim=max_paropt)
24 parameter (maxvdim=maxdim*(maxdim+1)/2)
25 logical fail,init,qua,rest,lint
26 common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
27 1 vy(maxdim),g(maxdim)
31 call convrt(n,x,.true.)
35 if (.not.rest) goto 20
37 print *,'reading restart information from the restart file'
41 if (nstep*100.lt.nvdim) nstep1=nstep+1
42 read(nfile) kkk,f,dftau0
44 222 read (nfile) g(i),y(i),x(i),s(i)
48 if (ne.gt.nvdim) ne=nvdim
49 111 read (nfile) (v(j),j=ns,ne)
51 print *,'restart information has been read in'
76 call dirmin(iout,iest,nfun,n,f,dftau0,d,y,x,s,fail,init,qua)
77 if (nprint.gt.2) write (iout,1040) (i,x(i),i=1,n)
80 if (nprint.gt.-1) write (iout,1100) it,nfun
81 if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
82 1100 format (//1x,' * * * * * * * * no point remarkably lower in ene',
83 1 'rgy can be found.'/1x,'calculation terminated in cycle',i3,
84 2 ' after',i6,' function & gradient evaluations')
85 call convrt(n,x,.false.)
97 if (dxi.gt.xmax) xmax=dxi
98 if (gi.gt.gmax) gmax=gi
100 163 gnorm=gnorm+gi*gi
104 if(nprint.gt.0)write(iout,1000)it,kk,xmax,xnorm,gmax,gnorm,f,nfun
105 if (nprint.gt.2) write (iout,1070) (i,y(i),i=1,n)
106 1070 format (/5x,' * * * * gradients'/100(1x,i5,5x,1pe12.4/))
107 if (nfun.lt.isave) goto 23
108 call convrt (n,x,.false.)
109 call psave(n,nvdim,nfun,nfile,kk,it,f,dftau0)
112 23 if (xmax.gt.tolx.or.gmax.gt.tolg) goto 17
113 if (nprint.gt.-1) write (iout,1010) it,nfun
114 if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n)
124 if (dabs(sy).gt.tolrec) goto 162
125 if (nprint.gt.-1) write (iout,1030)
126 1030 format (//' * * * delta(x) is perpendicular to delta(g). ',
127 1 'the program will try gradient direction * * *')
132 ind=(2*n-j)*(j-1)/2+i
133 10 vyi=vyi+v(ind)*y(j)
137 ind=(2*n-i)*(i-1)/2+j
138 11 vyi=vyi+v(ind)*y(j)
142 12 yvy=yvy+y(i)*vy(i)
149 v(ind)=v(ind)+si*s(j)-vyi*vy(j)
151 if (nprint.gt.4) call matout(nvdim,n)
153 if (v((2*n-i)*(i-1)/2+i).gt.0.) goto 18
154 if (nprint.gt.-1) write (iout,1080) i
157 1080 format (/1x,i3,'-th diagonal element of the matrix v is negative.'
158 1 ,' restart with gradient direction follows.')
162 ind=(2*n-j)*(j-1)/2+i
167 ind=(2*n-i)*(i-1)/2+j
173 call dirmin(iout,iest,nfun,n,f,dftau0,d,y,x,s,fail,init,qua)
175 if (.not.init) goto 19
176 if (nprint.gt.-1) write (iout,1110)
177 1110 format (/1x,'energy from dirmin differs but slightly from the'
178 1 ,' initial value. restart with gradient direction follows.')
187 if (xmax.lt.dxi) xmax=dxi
188 if (gmax.lt.gi) gmax=gi
190 161 gnorm=gnorm+gi*gi
193 if (nprint.gt.0) write (iout,1000) it,kk,xmax,xnorm,gmax,gnorm,f,
195 if (nprint.gt.2) write (iout,1040) (i,x(i),i=1,n)
196 1040 format (/1x,' * * * * parameters:'/100(1x,i5,5x,1pe14.6/))
197 if (nprint.gt.2) write (iout,1070) (i,y(i),i=1,n)
198 1000 format (/1x,'iteration:',i4,' direction:',i4,' max. element of ',
199 1 'delta(x): ',1pe14.6,' rms of delta(x): ',1pe14.6/' maximum ',
200 2 'element of gradient: ',1pe14.6,' rms of gradient: ',1pe14.6,
201 3 ' energy: ',1pe14.6,' fun. eval:',i6)
203 if (nfun.lt.isave) goto 22
204 call convrt(n,x,.false.)
205 call psave(n,nvdim,nfun,nfile,kk,it,f,dftau0)
208 22 if (xmax.gt.tolx.or.gmax.gt.tolg) goto 6
209 if (nprint.gt.-1) write (iout,1010) it,nfun
210 if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n)
211 if (nprint.gt.-1) write (iout,1111)
212 1111 format (//11x,'matrix v')
213 if (nprint.gt.-1) call matout(nvdim,n)
214 1010 format (//1x,'* * * * * * convergence has been achieved in ',
215 1 'cycle',i4,' after',i6,' function & gradient evaluations')
216 call convrt(n,x,.false.)
222 c if this point is reached, convergence hasn't been achieved.
223 if (nprint.gt.-1) write (iout,1020) maxit,nfun
224 if (nprint.gt.-1) write (iout,1070) (i,y(i),i=1,n)
225 if (nprint.gt.-1) write (iout,1111)
226 if (nprint.gt.-1) call matout(nvdim,n)
227 call convrt(n,x,.false.)
230 1020 format (//1x,' * * * * * * * convergence has not been achieved',
231 1 ' * * * * * * * '/1x,'number of cycles done:',i3,' number of ',
232 2 'function & gradient evaluations:',i6)
235 *******************************************************************************
236 subroutine grad(n,iest,g,x)
237 implicit real*8 (a-h,o-z)
238 include 'DIMENSIONS.ZSCOPT'
239 dimension g(max_paropt),x(max_paropt)
240 integer uiparm,nf,iest
242 double precision obg(max_paropt)
243 call targetgrad(n,x,nf,g,uiparm,obg,fdum)
246 *******************************************************************************
247 subroutine dirmin(iout,iest,nfun,n,f,dftau0,d,g,x,s,fail,init,qua)
248 implicit real*8 (a-h,o-z)
249 include 'DIMENSIONS.ZSCOPT'
250 parameter (maxdim=max_paropt)
251 logical fail,init,qua,weiter
252 common /dirit/ nprint,maxexp,maxlin,maxd,expand,tolqu,dermin,amov,
254 dimension x(n),s(n),d(n),g(n),z(maxdim),h(maxdim)
255 * write (2,*) 'dirmin has been called'
262 if (pom.gt.amx) amx=pom
264 tau=dmin1(1.0d0,amov/amx)
269 * write (2,'(2hh:/(i3,1pe15.5))' ) (i,h(i),i=1,n)
270 if (nprint.gt.1) write (iout,1020) dftau0
272 * write (2,'(2hh:/(i3,1pe15.5))' ) (i,h(i),i=1,n)
279 call grad(n,iest,g,x)
281 if (nprint.gt.1) write (iout,1000) ilin,tau,iexp,f1,f,dftau
282 1000 format (1x,'dirmin iter: ',i4,' tau=',1pe12.4,
283 1 ' contr. step:',i4/' actual f value: ',1pe14.6,' old value: ',
284 2 1pe14.6,' derivative: ',1pe14.6)
285 if (nprint.gt.4) write (iout,1080) (i,x(i),i=1,n)
286 if (nprint.gt.5) write (iout,1070) (i,g(i),i=1,n)
287 if (dabs(f1-f).gt.tolqu.and.amx*tau.ge.xmin) goto 21
294 21 if (f1.lt.f.or.dftau.gt.0.) goto 6
296 if (nprint.gt.-1) write(iout,1010) f
300 if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
303 1010 format (//' * * * * * dirmin failed. least function value is: ',
304 1 1pe14.6,' * * * * *')
305 1080 format (/' * * * * parameters:'/100(1x,i5,5x,1pe14.6/))
306 1070 format (/' * * * * gradients:'/100(1x,i5,5x,1pe14.6/))
307 1020 format (/1x,'initial value of directional derivative: ',1pe14.6)
308 6 if (dftau.gt.0.) goto 10
310 if (dabs(dftau/amx).le.dermin) goto 14
318 if (nprint.gt.-1) write (iout,1030)
319 1030 format (//1x,' * * * * * error in dirmin. positive directional ',
320 1 'derivative cannot be found * * * * *')
321 if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
324 10 term=f1-f-dftau0*tau
325 if (qua.or.dabs(dftau0/amx).lt.0.01.and.dabs(dftau/amx).lt.0.01)
327 difd=tau*(dftau-dftau0)
330 aa=(tau*difd-2.*term)/ttt
331 if (aa.eq.0.) goto 30
333 delta=4.*bb*bb-12.*aa*dftau0
334 if (delta.lt.0.) goto 30
336 tau1=(-2.*bb-delta)/(6.*aa)
337 tau2=(-2.*bb+delta)/(6.*aa)
338 if (nprint.gt.1) write (iout,1113) ilin,tau1,tau2
339 1113 format (/1x,'dirmin iter.:',i3,' cubic interpolation',
340 1 ' tau1=',1pe12.4,' tau2= ',1pe12.4)
341 if (tau1.gt.0.0d0.and.tau1.le.tau) goto 16
342 if (tau2.le.0.0d0.or.tau2.gt.tau) goto 30
345 30 tau1=-dftau0*tau*tau/(2.*term)
346 if (nprint.gt.1) write (iout,1040) ilin,tau1
347 if (tau1.lt.tau.and.tau1.gt.0.) goto 16
348 if (dabs(amx*tau1).lt.xmin) goto 145
359 1040 format (/1x,'dirmin iteration:',i4,' tau after quadr. interpo',
360 1 'lation: ',1pe12.4)
362 11 x(i)=s(i)+tau1*d(i)
363 if (nprint.gt.4) write (iout,1080) (i,x(i),i=1,n)
364 if (nprint.gt.5) write (iout,1070) (i,g(i),i=1,n)
366 call grad(n,iest,g,x)
368 dftn=dabs(dftau1/amx)
369 if (nprint.gt.1) write (iout,1112) f2,dftau1
370 1112 format (1x,'function value: ',1pe14.6,' derivative: ',1pe14.6)
371 if (f2.lt.f.or.dftau1.gt.0.) goto 20
372 if (dabs(f2-f).gt.tolqu.and.dabs(amx*tau1).gt.xmin.and.dftn.ge.
383 if (f2.lt.f) init=.false.
385 if (dabs(f2-f).gt.tolqu.and.dabs(amx*tau1).ge.xmin.and.dftn.ge.
397 if (nprint.gt.1) write (iout,1111) taudf
398 1111 format (/1x,'difference: ',1pe10.2)
399 if (dabs(taudf).lt.tolqu.or.dftn.lt.dermin) goto 144
400 if (dftau1.lt.0.) goto 13
419 if (nprint.gt.-1) write (iout,1050)
420 1050 format (//1x,' * * * * nonconvergence in dirmin, calculation ',
421 1 'terminated. * * * *')
423 if (nprint.gt.-1) write (iout,1070) (i,g(i),i=1,n)
430 18 if (nprint.gt.1) write (iout,1060) ilin
431 1060 format (/1x,'dirmin converged in',i3,' iterations')
434 *******************************************************************************
435 double precision function dirder(n,g,d)
436 implicit real*8 (a-h,o-z)
438 c calculates the directional derivative. d is the direction vector,
439 c g - gradient vector
447 *******************************************************************************
448 double precision function fun(n,x,nfun)
449 implicit real*8 (a-h,o-z)
450 include 'DIMENSIONS.ZSCOPT'
451 parameter (maxdim=max_paropt)
452 dimension x(max_paropt)
455 double precision f,obf
458 call targetfun(n,x,nfun,f,uiparm,obf,fdum)
462 *******************************************************************************
463 subroutine convrt(n,x,dir)
464 implicit real*8 (a-h,o-z)
468 *******************************************************************************
469 subroutine dfpopt(inp,iout,ifirst,n,xopt)
470 implicit real*8 (a-h,o-z)
472 c this subroutine reads in the iteration parameters, calls
473 c the minimization procedures (davejr or dfp or both of them one by
474 c one) and writes the final coordinates to file ares.xyz.
476 include 'DIMENSIONS.ZSCOPT'
477 parameter (maxdim=max_paropt)
478 parameter (maxvdim=maxdim*(maxdim+1)/2)
479 common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
480 1 vy(maxdim),g(maxdim)
481 logical conv,qua,rest,lint
482 common /dirit/ nprint,maxexp,maxlin,maxd,expand,tolqu,dermin,amov,
485 integer option,odfp,onr
486 dimension xopt(maxdim)
487 data odfp /4h dfp/,onr /4h nr/, ncu /4h cub/,nrest /4hrest/
488 print *,'ce2 has been called'
489 print *,"DFPOPT: n",n
491 c if (ifirst.eq.0) goto 1
492 c read (inp,1001) option,icub,irest
494 c if (irest.eq.nrest) rest=.true.
496 c if (icub.eq.ncu) qua=.false.
497 c read (inp,1010) maxit,itolx,itolg,maxlin,maxexp,maxd,itolq,itolr,
498 c 1 iderm,nprint,ixmin,nsave,nfile,iest
499 c read (inp,1030) expand,amov
500 c if (nsave.eq.0) nsave=10000
502 c if (itolx.eq.0) itolx=3
503 c if (itolg.eq.0) itolg=1
504 c if (itolr.eq.0) itolr=6
505 c if (iderm.eq.0) iderm=7
506 c if (itolq.eq.0) itolq=4
507 c if (ixmin.eq.0) ixmin=5
508 c if (maxit.eq.0) maxit=5
509 c if (maxlin.eq.0) maxlin=10
510 c if (maxd.eq.0) maxd=10
511 c if (maxexp.eq.0) maxexp=10
512 c if (expand.eq.0.) expand=2.0d0
513 c if (amov.eq.0.) amov=.25d0
514 c tolx=10.0d0**(-itolx)
515 c tolg=10.0d0**(-itolg)
516 c xmin=10.0d0**(-ixmin)
517 c tolqu=10.0d0**(-itolq)
518 c dermin=10.0d0**(-iderm)
519 c tolrec=10.0d0**(-itolr)
520 c if (option.eq.odfp) goto 4
543 4 write (iout,2000) maxit,tolx,tolg,maxlin,maxexp,maxd,tolqu
544 write (iout,2002) expand,amov,xmin,dermin,tolrec
545 if (iest.gt.0) write (iout,2060)
546 write (iout,2040) nsave,nfile
547 if (rest) write (iout,2050) nfile
548 if (qua) write (iout,2020)
549 if (.not.qua) write (iout,2030)
550 1 call dfp(inp,iout,nfile,nsave,maxit,tolx,tolg,tolrec,qua,rest,
561 2000 format (//'welcome to dfp routine'//'the following ',
562 1 'options are used:'//
563 2 'maximum number of cycles:',i4/'tolerance on coordinates:',
564 3 1pe10.2/'tolerance on gradient:',1pe10.2/
565 4 'maximum number of iterations in linesearch:',i4/
566 5 'maximum number of unsuccesful function evaluations:',i4/
567 6 'maximum number of steps in positive derivative search:',i4/
568 7 'tolerance in linesearch:',1pe10.2)
569 2002 format ('expansion coefficient in linesearch:',1pe10.2/
570 9 'maximum movement of atoms:',1pe10.2/
571 a 'minimum change of coordinates considered significant:',1pe10.2/
572 c 'minimum value of directional derivative considered different ',
573 d 'from zero:',1pe10.2/
574 e 'lower bound of the scalar product (delta(x),delta(g)):',
576 2003 format (/'not quite succesful... perhaps you will have a '
579 2060 format (/'estimated partials used.')
580 2020 format (/'quadratic interpolation will be used in line',
582 2030 format (/'cubic or quadratic interpolation will be used ',
583 1 'in line minimization')
584 2040 format (/'coordinates, gradients, and matrix v will be saved',
585 1 ' appr. every',i4,' energy and gradient evaluations on file',i3)
586 2050 format (//'restart to be done using information saved ',
590 *******************************************************************************
591 subroutine csave(lint)
595 *******************************************************************************
596 subroutine psave(n,nvdim,nfun,nfile,kk,it,f,dftau0)
597 implicit real*8 (a-h,o-z)
598 include 'DIMENSIONS.ZSCOPT'
599 parameter (maxdim=max_paropt)
600 parameter (maxvdim=maxdim*(maxdim+1)/2)
601 logical fail,init,qua,rest,lint
602 common/bigmat/ v(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
603 1 vy(maxdim),g(maxdim)
604 common /dirit/ nprint
605 common /restart/ restfile
606 character*16 restfile
607 common inp,iout,ipn,iprint
611 if (nstep*100.lt.nvdim) nstep1=nstep+1
612 write (nfile) kk,f,dftau0
614 222 write (nfile) g(i),y(i),x(i),s(i)
618 if (ne.gt.nvdim) ne=nvdim
619 111 write (nfile) (v(j),j=ns,ne)
621 write (*,1010) it,nfun,f,restfile
622 1010 format (/1x,'cycle: ',i2,' energy evaluations: ',i5,' value:',
623 1 f20.10,'.'/' information for restart saved on file ',a16)
626 *******************************************************************************
627 subroutine matout(ndim,n)
628 implicit real*8 (a-h,o-z)
629 include 'DIMENSIONS.ZSCOPT'
630 parameter (maxdim=max_paropt)
631 parameter (maxvdim=maxdim*(maxdim+1)/2)
632 logical fail,init,qua,rest,lint
633 common/bigmat/ a(maxvdim),s(maxdim),x(maxdim),y(maxdim),d(maxdim),
634 1 vy(maxdim),g(maxdim)
635 character*6 line6 / '------' /
636 character*12 line12 / '------------' /
641 write (iout,1000) (k,k=i,nlim)
642 write (iout,1020) line6,(line12,k=i,nlim)
643 1000 format (/7x,6(i7,5x))
644 1020 format (a6,6a12)
649 3 b(k-i+1)=a(n*(k-1)-k*(k-1)/2+j)
650 write (iout,1010) j,(b(k),k=1,jlim1)
653 1010 format (i3,3x,6(1pe11.4,1x))