c$DEBUG c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c RGDMIN - THIS PROGRAM PREDICTS STRUCTURES OF DIMERS CONSISTING OF c A RARE-GAS ATOM AND A MOLECULE BY MEANS OF THE 'SIMPLE' MODEL c WHICH ATTEMPTS TO EMULATE THE GEOMETRICAL ANISOTROPY OF THE c DISPERSIVE INTERACTION RESPONSIBLE FOR MOST OF THE ATTRACTION. c Interaction is partitioned over atomic sites of the acceptor c molecule and assumeed to be described by London type r**(-6) c terms in which atomic polarisabilities are replaced by cubes c of covalent radii. Attraction is counterbalanced by hard c shell repulsion for which all atoms are assumed to be c surrounded by spheres of size determined by their van der c Waals radius. The sum of attractive and repulsive c terms so defined is minimised by reorienting the Rg-atom c relative to the molecule. c C ver. 7/93a ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c c c Program has following features: c c - automatic assignment of model parameters which can be overriden c - fixed geometry calculation c - full and limited minimisation (eg. for potential characterisation) c - Rg-atom translation vector can be either cartesian or polar c - input of translation vector from keyboard or file c - output to screen and two types of files c c RGDMIN was only run on PC's with Microsoft 3.31 Fortran, but should be c easily transferrable to other environments as programs it was derived c from were run on other machines and F77 compilers. c c............................................................................. c c RGDMIN can access up to four different files during operation. These are: c c 1. File containing molecular data with user provided name c 2. Optional file containing initial positions of B, user named c 3. Main output file containing results of minimisations, user named c 4. Ancillary output file containing details of individual iterations, c this has a fixed name MONITOR c c STRUCTURE OF THE DATA FILE: c c LINE 1 - comment line c c LINE 2 - NA,UNITIN,UNITOT,IMONIT (free format) c where: c NA is the number of centres in the acceptor molecule (atoms+ c dummy repulsive centres) c UNITIN determines whether the input atomic coordinates are in c Angstroms (=1), or bohrs (=0) (internally the program works c in atomic units) c UNITOT determines whether the output atomic coordinates are to c be in Angstr (=1), or bohr (=0). This parameter also defines c the units to be used for later input of translation vector c components for rg. position c IMONIT determines whether a brief (=1) or full (=2) MONITOR c listing is required. Brief listing summarises only the c minimised energies, full listing gives progress of c minimisation iteration by iteration and all occurrences of c repulsive surface clashes (in form: numbers of atoms in c an atom pair, sum of their v.d.W. radii,actual separation, and c current value of the repulsive term) c c LINE 3.. - for i=1,NA these lines are each to contain: c MM, x, y, z, [rcovcube, rvdwaal] (free format) c where: c MM is a two character, left justified atom name, this is case c sensitive and follows standard chemical nomenclature, c x,y,z are coordinates of the atom in units defined by UNITIN, c If the program recognises the atom name then the attractive and c repulsive parameter values are assigned internally. If MM is c not recognised then the two parameters r(cov)**3 (au) and c r(vdWaal) (angstr) are to be provided on the same line. This c procedure allows definition of dummy repulsive centres (in c which case r(cov)**3 should be set to zero), and also gives the c possibility for overriding internal parameter values. The c values of x,y,z for the Rg. are normally all set to zero. c c LAST LINE - line for the Rare-gas atom to be complexed, also in the form: c MM, x, y, z, [polarizability, rvdwaal] (free format) c c c INPUT OF INITIAL ORIENTATION OF THE RG-ATOM FROM SCREEN: c c I - selection of computation mode and type of orientation vector: c I=-1 or 1 cartesian translation vector c I=-2 or >=2 polar translation vector c I<0 calculation for fixed geometry c I>0 minimisation c I=1 or 2 full minimisation c I=3,4,5 limited minimisation: only one of R,theta,phi resp. is c fitted c P, Q, R /Angstr or bohr, degrees: c if I=+-1 then P,Q,R are cartesian components of c.m. A-B c separation vector, their polar counterparts used for c calculation c if I=+-2,3,4,5 then P,Q,R are polar coordinates R,Theta,phi c defining the c.m. separation vector, serving as minimisation C parameters par(1),par(2),par(3) c c c INPUT OF INITIAL ORIENTATION OF THE RG-ATOM FROM DISK: c c This should contain the same information as that input from screen with c the difference that a comment line should precede lines of orientation c data and that the control parameter I should be input in the same line as c P,Q,R. Two line blocks of the form: c comment line c I, P, Q, R c are therefore required and can be repeated as many times as necessary. c The comment line is not used by the program (can be empty) and is only c for legibility of the data file. This type of input can be used, for c example, to define a sphere of starting positions to ensure that no c minima are overlooked. c c_____________________________________________________________________________ c c implicit real*8 (a-h,o-z) parameter (natoms=11,natmax=50) common /blk0/vdwrad(natmax),x(natmax,3),q(natmax),radius, . pi,conv,na,nb,eattr,erep common /blkf/ifit common /blkder/iflag integer unitin,unitot real*8 par(6),r(3),xb(3),atwt(natmax) character*25 filnam character*76 coment,line character*10 names(5) character*2 namat(natoms),namrg(5),nameat(natmax),mm real*8 vdwaal(natoms),rqcov(natoms),rgpol(5),rgrad(5),rgmas(5), * atomas(natoms) c c...atomic data: c van der Waals radii/Angstrom: A.Bondi, J.Phys.Chem. 68,441(1964) c cubes of covalent radii/au derived from covalent radii from: c L.Pauling - The Nature of the Chemical Bond (except for value c of 0.505 for hydrogen which is derived by considering dispersion c in Ar-HCl) c data namat/'H ','C ','N ','O ','F ','Si','P ','S ','Cl','Br','I '/ data atomas/1.007825d0,12.d0,14.003241d0,15.9949150d0, * 18.9984046d0,27.976929d0,30.973765d0,31.9720737d0,34.968851d0, * 78.918329d0,126.904470d0/ data vdwaal/1.20,1.70,1.55,1.52,1.47,2.10,1.80,1.80,1.75, * 1.85,1.98/ data rqcov/0.505,3.1,2.7,2.7,2.5,10.8,9.0,7.6,6.6,10.0,15.9/ c c...rare gas data: c polarizabilities/au: C.G.Gray & K.E.Gubbins, Theory of Molecular Fluids, c Clarendon Press (1984) c van der Waals radii/Angstrom: also from Bondi but incremented by c 0.12 Angstrom. This is the difference between 1.88 and the c effective radius of 2 Angstroms derived for Ar in Rg-molecule c dimers on using to Bondi's radii for atoms in the molecule. c data namrg/'He','Ne','Ar','Kr','Xe'/ data rgrad/1.52,1.66,2.00,2.14,2.28/ data rgpol/1.39,2.67,11.08,16.76,27.1/ data rgmas/4.0026031d0,19.9924405d0,39.9623842d0, * 85.911503d0,131.904161d0/ c nrgas=5 pi=3.1415926535897932d0 conv=pi/180d0 bohr=0.5291771d0 c c c...Open input and output files, read molecular data c write (*,297) 297 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | RGDMIN - Minimisation of interaction energy for a dimer', * ' between a rare',T79,'|'/ * ' | gas atom and a molecule',T79,'|'/ * ' |',76(1H_),'|'/' version 7/1993',T64,'Zbigniew KISIEL'//) c 200 write(*,212)' Name of file containing the data: ' 212 format(1x/a,$) read(*,'(a)',err=200)filnam open(2,file=filnam,err=200) c 287 write(*,212)' Name of results file: ' read(*,'(a)',err=287)filnam open(3,file=filnam,status='unknown') open(4,file='monitor',status='unknown') write(3,211)' ' write(4,211)' ' write(3,199) 199 format(1x,78(1h_)) write(3,290) 290 format(1x/' MINIMISATION OF INTERACTION ENERGY OF RARE-GAS AT' * ,'OM AND A MOLECULE,'/' consisting of a sum of London type', * ' pairwise attraction and hard' * /' van der Waals shell repulsion. Polarizabilities of bound', * ' atoms are'/' simulated by cubes of covalent radii.') write(3,199) c c c...READ DATA FILE, comment and options c read(2,221)coment 221 format(a) write(*,375)coment write(3,375)coment 375 format(1x/1x,a//' atom',7x,'x',11x,'y',11x,'z', * 6x,'r.cov**3/au',3x,'r.vdW/Angstr'/) read(2,*,err=382)na,unitin,unitot,imonit if(unitin.ne.1.and.unitin.ne.0)unitin=1 if(unitot.ne.1.and.unitot.ne.0)unitot=1 if(imonit.ne.1.and.imonit.ne.2)imonit=1 c c...data on the molecule c ramax=0.d0 do 1 i=1,na read(2,221,err=380)line mm=line(1:2) line=line(3:) call cond(line) do 360 j=1,natoms if(namat(j).eq.mm)goto 352 360 continue c nameat(i)=mm read(line,353,err=380)(x(i,j),j=1,3),q(i),vdwrad(i) vdwrad(i)=vdwrad(i)/bohr goto 355 c 352 nameat(i)=namat(j) vdwrad(i)=vdwaal(j)/bohr atwt(i)=atomas(j) q(i)=rqcov(j) read(line,353,err=380)(x(i,j),j=1,3) 353 format(5f15.10) c 355 if(unitin.eq.1)then do 354 j=1,3 x(i,j)=x(i,j)/bohr 354 continue endif radius=dsqrt(x(i,1)**2+x(i,2)**2+x(i,3)**2)+vdwrad(i) if(radius.gt.ramax)ramax=radius c write(*,288)nameat(i),(x(i,j)*bohr,j=1,3),q(i),vdwrad(i)*bohr write(3,288)nameat(i),(x(i,j)*bohr,j=1,3),q(i),vdwrad(i)*bohr 288 format(1x,a2,3f12.6,f15.6,f12.6) 1 continue c c...data on the rare-gas atom c nb=1 i=na+nb read(2,221,err=383)line mm=line(1:2) line=line(3:) call cond(line) do 357 j=1,nrgas if(namrg(j).eq.mm)goto 356 357 continue c nameat(i)=mm read(line,353,err=383)(x(i,j),j=1,3),q(i),vdwrad(i) vdwrad(i)=vdwrad(i)/bohr goto 358 c 356 nameat(i)=namrg(j) vdwrad(i)=rgrad(j)/bohr q(i)=rgpol(j) atwt(i)=rgmas(j) read(line,353,err=383)(x(i,j),j=1,3) 358 if(unitin.eq.1)then do 359 j=1,3 x(i,j)=x(i,j)/bohr 359 continue endif write(*,288)nameat(i),(x(i,j)*bohr,j=1,3),q(i),vdwrad(i)*bohr write(3,288)nameat(i),(x(i,j)*bohr,j=1,3),q(i),vdwrad(i)*bohr c radius=ramax+0.2d0+vdwrad(i) write(3,295)radius*bohr write(*,295)radius*bohr 295 format(1x/' Standoff c.m. distance =',f10.5,' Angstr.') goto 302 c 380 write(*,'(1x//''**** ERROR on reading line for atom'',i3//)')i stop 382 write(*,'(1x//''**** ERROR on reading the second line''//)') stop 383 write(*,'(1x//''**** ERROR on reading the Rg-atom line''//)') stop c 302 close(2) c c c...determine whether initial orientation of Rg-atom will be input from disk c or from keyboard c 251 nrun=0 write(*,361)' Orientation of Rg-atom from keyboard (=0) or', * ' from disk (=1) ? ' 361 format(1x/2a,$) read(*,'(i5)',err=251)ikbdsk if(ikbdsk.ne.1.and.ikbdsk.ne.0)goto 251 if(ikbdsk.eq.0)goto 252 253 write(*,212)' Name of file containing reorientation data: ' read(*,'(a)',err=253)filnam open(2,file=filnam,err=253) c c c...beginning of loop for interaction energy minimisation for various c starting relative orientations of molecule and Rg-atom c c 201 if(ikbdsk.eq.0)goto 252 c c...initial orientation of Rg-atom from from disk c read(2,'(a)',err=393,end=256)filnam 393 read(2,*,err=393,end=256)iflag,(r(i),i=1,3) goto 254 256 close(2) goto 220 c c...initial orientation of Rg-atom from keyboard c 252 write(*,350) 350 format(/' TYPE OF CALCULATION, SPECIFY VALUE FROM -2 to 5:' * /' -ve: fixed calc only 1: cartesian Rg ', * 'translation vector ' * /' +ve: fit 2,3,4,5: polar Rg translation ', * 'vector' * /' 555: exit 3,4,5: only one of R,theta,', * 'phi fitted resp. '/' ..... ',$) read(*,*,err=252)iflag if(iflag.eq.555)goto 220 if(iflag.lt.-2.or.iflag.eq.0.or.iflag.gt.5)goto 252 351 write(*,211) * ' Vector from molecule to Rg (x,y,z or R,theta,phi):' 211 format(1x/1x,a) read(*,*,err=351)(r(i),i=1,3) c c...conversion to bohr and of Rg vector to polar coordinates (if necessaary) c 254 ifit=1 if(iflag.lt.0)then ifit=0 iter=0 iflag=-iflag endif if(iflag.ge.2.and.iflag.le.5)then par(1)=r(1) par(2)=r(2)*conv par(3)=r(3)*conv r(1)=par(1)*dsin(par(2))*dcos(par(3)) r(2)=par(1)*dsin(par(2))*dsin(par(3)) r(3)=par(1)*dcos(par(2)) endif c if(unitot.eq.0)goto 210 do 215 i=1,3 r(i)=r(i)/bohr 215 continue c 210 par(1)=dsqrt(r(1)**2+r(2)**2+r(3)**2) if(par(1).eq.0.d0)par(1)=1.d-10 par(2)=dacos(r(3)/par(1)) if(r(1).eq.0.0d0)r(1)=1.d-10 par(3)=datan2(r(2),r(1)) par(4)=0.d0 par(5)=0.d0 par(6)=0.d0 c write(*,199) write(*,281)' INITIAL',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv write(3,199) write(3,281)' INITIAL',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv if(imonit.ne.2)goto 291 write(4,199) if(imonit.ne.2)goto 291 write(4,281)' INITIAL',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv c c c...MINIMISATION c 291 nrun=nrun+1 if(ifit.eq.1)then 390 vtot=ESINTE(par) if(erep.gt.25.d0)goto 390 CALL DFPMIN(par,6,1.d-8,iter,fret) endif vtot=ESINTE(par) c c c...FINAL RESULTS, energy c write(3,280)' MINIMISED ENERGY =',vtot,vtot*219476.d0,iter,nrun 280 format(1x/a,f12.8,' <<',F11.2,'cm-1 NIT=',i3,' RUN=',i3,'>>') write(4,292)' MINIMISED ENERGY =',vtot, * vtot*219476.d0,iter,nrun 292 format( a,f12.8,' <<',f11.2,'cm-1 NIT=',i3,' RUN=',i3,'>>') r(1)=par(1)*dsin(par(2))*dcos(par(3)) r(2)=par(1)*dsin(par(2))*dsin(par(3)) r(3)=par(1)*dcos(par(2)) c c...coordinates and energy decomposition c write(*,373) write(3,373) 373 format(/1x,'CARTESIAN COORDINATES AND PAIRWISE ENERGIES FOR ', * 'THE RG-MOLECULE SYSTEM:'//' Atom',6x,'x',12x,'y',12x,'z', 8x, * 'mass',9x,'E.attr',6x,'E.rep'/) rm=1.d0 if(unitot.eq.1)rm=bohr do 370 i=1,3 xb(i)=x(na+1,i)+r(i) 370 continue c do 30 i=1,na sep=dsqrt((x(i,1)-xb(1))**2+(x(i,2)-xb(2))**2 * +(x(i,3)-xb(3))**2) peattr=-1.5d0*q(i)*q(na+1)/sep**6 sigij=vdwrad(i)+vdwrad(na+1) perep=25.d0*(1.d0+dtanh((sigij-sep-0.6d0)/0.1d0)) write(*,371)nameat(i),(x(i,j)*rm,j=1,3),atwt(i),peattr,perep write(3,371)nameat(i),(x(i,j)*rm,j=1,3),atwt(i),peattr,perep 30 continue write(*,206)nameat(na+1),(xb(j)*rm,j=1,3),atwt(na+1) write(3,206)nameat(na+1),(xb(j)*rm,j=1,3),atwt(na+1) 371 format(1x,a,3f13.6,f10.5,f14.7,f12.7) 206 format(1x,a,3f13.6,f10.5) c eint=eattr+erep write(*,245)eattr,erep,eattr*219476.D0,erep*219476.D0, * eint,eint*219476.D0 write(3,245)eattr,erep,eattr*219476.D0,erep*219476.D0, * eint,eint*219476.D0 245 format(38x,'totals:',f16.8,f14.8/38x,' ',f12.2,' cm-1', * f11.2,' cm-1'/38x,' E.int:',f16.8,f12.2,' cm-1') c c...translation vector c if(ifit.eq.1)then write(3,281)' MINIMISED',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv write(*,281)' MINIMISED',(r(i)*bohr,i=1,3),par(1)*bohr, * par(2)/conv,par(3)/conv endif 281 format(1x/a,' DISPLACEMENT VECTOR FOR THE RG-ATOM:'/ * ' cartesian: ',3f15.6/' polar: ', * 3f15.6) goto 201 c c 220 write(3,199) write(4,199) close(3) close(4) c stop end c c_____________________________________________________________________________ c subroutine cond(line) c c This routine conditions a line of numerical data to simulate free c format input from internal files (not allowed by the F77 standard) c character*76 line c nc=76 iput=1 inum=0 do 1 i=1,nc if(line(i:i).eq.'-'.or.line(i:i).eq.'.'.or. * (line(i:i).ge.'0'.and.line(i:i).le.'9'))then line(iput:iput)=line(i:i) iput=iput+1 inum=1 else if(inum.eq.1)then line(iput:iput)=',' iput=iput+1 inum=0 endif endif 1 continue if(iput.eq.nc+1)return do 2 i=iput,nc line(i:i)=' ' 2 continue c return end c c_____________________________________________________________________________ c FUNCTION ESINTE(par) c c Calculation of interaction energy between molecules A and B being c a sum of attractive London type pair terms and hard sphere repulsive c terms c implicit real*8 (a-h,o-z) parameter (natmax=50) common /blk0/vdwrad(natmax),x(natmax,3),q(natmax),radius, . pi,conv,na,nb,eattr,erep real*8 par(6) common /blkf/ifit c c...interaction energy c CALL ESUM(par,ifit) esinte=eattr+erep c return end c c_____________________________________________________________________________ c SUBROUTINE ESINTD(par,grad,f) c c Evaluation of derivatives GRAD of interaction energy w.r.t. the c minimisation parameters in PAR. F is the interaction energy at c PAR. c implicit real*8 (a-h,o-z) real*8 par(6),grad(6),delta(6) common /blkder/iflag data delta/0.0001d0,5*0.00017d0/ c do 1 n=1,3 par(n)=par(n)+delta(n) fd=esinte(par) grad(n)=(fd-f)/delta(n) par(n)=par(n)-delta(n) 1 continue grad(4)=0.d0 grad(5)=0.d0 grad(6)=0.d0 if(iflag.ge.3)then do 2 n=1,3 if(n.eq.iflag-2)goto 2 grad(n)=0.d0 2 continue endif c return end c c_____________________________________________________________________________ c SUBROUTINE ESUM(par,ifit) c c c Sum of interaction energy over pairs of sites c implicit real*8 (a-h,o-z) parameter (natmax=50) common /blk0/vdwrad(natmax),x(natmax,3),q(natmax),radius, . pi,conv,na,nb,eattr,erep real*8 par(6),r(3) C erep=0.d0 eattr=0.d0 c c...'dirty' reset of R that is much too long and normalisation of Theta c and Phi c par(1)=dabs(par(1)) if(ifit.eq.1.and.par(1).gt.radius)par(1)=radius CALL NORMAL(1,par(2)) CALL NORMAL(2,par(3)) r(1)=par(1)*dsin(par(2))*dcos(par(3)) r(2)=par(1)*dsin(par(2))*dsin(par(3)) r(3)=par(1)*dcos(par(2)) c c...loop over all atoms of molecule A c do 12 i=1,na c c...for every atom of A loop over all atoms of B c do 14 j=na+1, na+nb sepsq=0.0d0 do 15 k=1,3 term=r(k)+x(j,k)-x(i,k) sepsq=sepsq+term*term 15 continue sep=dsqrt(sepsq) c c...attractive term c eattr=eattr-1.5d0*q(i)*q(j)/sep**6 c c...repulsive term (from S.L.Price,A.J.Stone, J.Chem.Phys. 86,2859(1987)) c sigij=vdwrad(i)+vdwrad(j) if(sep.gt.2.d0*sigij)goto 14 erep=erep+25.d0*(1.d0+dtanh((sigij-sep-0.6d0)/0.1d0)) c 14 continue 12 continue c c...'dirty' reset of R that is too short c if(ifit.eq.1.and.erep.gt.25.d0)par(1)=par(1)+1.d0 c return end c c_____________________________________________________________________________ c SUBROUTINE DFPMIN(P,N,FTOL,ITER,FRET) C C Given a starting point P that is a vector of length N, the Broyden- C Fletcher-Goldfarb-Shanno varian of Davidon-Fletcher-Powell minimisation C is performed on a functin FUNC, using its gradient as calculated by a C routine DFUNC. The convergence requirement on the function value is C input as FTOL. Returned quantities are P (the location of the minimum), C ITER (the number of iterations that were performed), and FRET (the minimum C value of the function). The routine LINMIN is called to perform line C minimisations. C C Minimisation routines: DFPMIN, LINMIN, MNBRAK, BRENT and F1DIM C come from Numerical Recipes. Modified or added lines are indicated by a C 'ZK' in the comment field. In addition all standard function names C have been changed to explicit double precision form. C C IMPLICIT REAL*8 (A-H,O-Z) ZK PARAMETER (NMAX=6,ITMAX=200,EPS=1.d-10) common /mon/imonit ZK DIMENSION P(N),HESSIN(NMAX,NMAX),XI(NMAX),G(NMAX),DG(NMAX), * HDG(NMAX) pi=3.1415926535897932d0 ZK conv=pi/180d0 ZK bohr=0.5291771d0 ZK C C...calculate starting function value and gradient C FP=ESINTE(P) ZK CALL ESINTD(P,G,FP) ZK write(*,30)0,fp,p(1)*bohr,(p(i)/conv,i=2,3), ZK * g(1)*bohr,(g(i)/conv,i=2,3) ZK if(imonit.eq.2)write(4,30)0,fp,p(1)*bohr,(p(i)/conv,i=2,3), ZK * g(1)*bohr,(g(i)/conv,i=2,3) ZK C C...initialise the inverse Hessian to the unit matrix and set initial C line direction XI C DO 12 I=1,N DO 11 J=1,N HESSIN(I,J)=0.d0 11 CONTINUE HESSIN(I,I)=1.d0 XI(I)=-G(I) 12 CONTINUE C C...main loop over the iterations C DO 24 ITER=1,ITMAX CALL LINMIN(P,XI,N,FRET) if(imonit.eq.2)write(4,30)iter,fret,p(1)*bohr, ZK * (p(i)/conv,i=2,3),g(1)*bohr,(g(i)/conv,i=2,3) ZK write(*,30)iter,fret,p(1)*bohr,(p(i)/conv,i=2,3), ZK * g(1)*bohr,(g(i)/conv,i=2,3) ZK 30 format(' iteration no.',i2,36x,'E.int = ',f16.10/ ZK * 14x,'R,Theta,Phi: ',3f12.6/14x,'Derivatives: ',3f12.6) ZK C C...return if converged to a minimum IF(2.d0*DABS(FRET-FP).LE.FTOL*(DABS(FRET)+DABS(FP)+EPS)) * RETURN C C...save old function value and the old gradient FP=FRET DO 13 I=1,N DG(I)=G(I) 13 CONTINUE C C...get new function value and gradient FRET=ESINTE(P) ZK CALL ESINTD(P,G,FRET) ZK fcalc=fret ZK C C...compute difference of gradients DO 14 I=1,N DG(I)=G(I)-DG(I) 14 CONTINUE C C...and difference times the current matrix DO 16 I=1,N HDG(I)=0.d0 DO 15 J=1,N HDG(I)=HDG(I)+HESSIN(I,J)*DG(J) 15 CONTINUE 16 CONTINUE C C...calculate dot products for the denominators FAC=0.d0 FAE=0.d0 DO 17 I=1,N FAC=FAC+DG(I)*XI(I) FAE=FAE+DG(I)*HDG(I) 17 CONTINUE C C...and make the denominators multiplicative FAC=1.d0/FAC FAD=1.d0/FAE C C...the vector which makes BFGS different from DFP DO 18 I=1,N DG(I)=FAC*XI(I)-FAD*HDG(I) 18 CONTINUE C C...the BFGS updating formula DO 21 I=1,N DO 19 J=1,N HESSIN(I,J)=HESSIN(I,J)+FAC*XI(I)*XI(J) * -FAD*HDG(I)*HDG(J)+FAE*DG(I)*DG(J) 19 CONTINUE 21 CONTINUE C C...now calculate the next direction to go DO 23 I=1,N XI(I)=0.d0 DO 22 J=1,N XI(I)=XI(I)-HESSIN(I,J)*G(J) 22 CONTINUE 23 CONTINUE C C...and go back to another iteration 24 CONTINUE C PAUSE 'too many iterations in DFPMIN' C RETURN END C C_____________________________________________________________________________ C SUBROUTINE NORMAL(I,ANGLE) C C Normalisation of an angle (in radians) to range -PI to PI for I=1 and C 0 to 2*PI for I=2 C IMPLICIT REAL*8 (A-H,O-Z) DATA PI,TPI/3.141592654D0,6.283185307D0/ C IF(I.EQ.1)THEN A=DABS(ANGLE) IF(A.LE.PI)RETURN A=A-DINT(A/TPI)*TPI ANGLE=DSIGN(A,ANGLE) IF(ANGLE.LT.-PI)ANGLE=ANGLE+TPI IF(ANGLE.GT.PI)ANGLE=ANGLE-TPI RETURN ELSE IF(I.EQ.2)THEN A=DABS(ANGLE) A=A-DINT(A/TPI)*TPI ANGLE=DSIGN(A,ANGLE) IF(ANGLE.LT.0.D0)ANGLE=ANGLE+TPI RETURN ENDIF C RETURN END C C_____________________________________________________________________________ C SUBROUTINE LINMIN(P,XI,N,FRET) C C Given an N dimensional point P and an N dimensional direction XI, LINMIN C moves and resets P to where the function FUNC(P) takes on a minimum C along the direction XI from P, and replaces XI by the actual vector C displacement that P was moved. Also returns FRET the value of FUNC at C the returned location P. This is actually all accomplished by calling C the routines MNBRAK and BRENT. C IMPLICIT REAL*8 (A-H,O-Z) ZK PARAMETER (NMAX=6,TOL=3.d-4) EXTERNAL F1DIM DIMENSION P(N),XI(N) COMMON /F1COM/ PCOM(NMAX),XICOM(NMAX),NCOM C C...set up the common block C NCOM=N DO 11 J=1,N PCOM(J)=P(J) XICOM(J)=XI(J) 11 CONTINUE C C...initial guess for brackets C AX=0.d0 XX=1.d0 BX=2.d0 CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM) FRET=BRENT(AX,XX,BX,F1DIM,TOL,XMIN) C C...construct the vector results to return C DO 12 J=1,N XI(J)=XMIN*XI(J) P(J)=P(J)+XI(J) 12 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNC) C C Given a function FUNC, and given distinct initial points AX and BX, this C routine searches in the downhill direction (defined by the function as C evaluated at the initial points) and returns new points AX, BX, CX which C bracket a minimum of the function. Also returned are the function values C at the three points , FA, FB, and FC. C IMPLICIT REAL*8 (A-H,O-Z) ZK C C...The first parameter is the default ratio by which successive intervals C are magnified; the second is the maximum magnification allowed for a C parabolic-fit step C PARAMETER (GOLD=1.618034d0, GLIMIT=100.d0, TINY=1.d-20) C FA=FUNC(AX) FB=FUNC(BX) C C...Switch roles of A and B so we can go downhill in the direction from A to B C IF(FB.GT.FA)THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM ENDIF C C...First guess for C C CX=BX+GOLD*(BX-AX) FC=FUNC(CX) C C...'DO WHILE': keep returning here until we bracket C 1 IF(FB.GE.FC)THEN C C...compute U by parabolic extrapolation from A,B,C. TINY is used to prevent C any possible division by zero R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) U=BX-((BX-CX)*Q-(BX-AX)*R) * /(2.d0*DSIGN(DMAX1(DABS(Q-R),TINY),Q-R)) ULIM=BX+GLIMIT*(CX-BX) C C...got a minimum between B and C, which will exit IF((BX-U)*(U-CX).GT.0.d0)THEN FU=FUNC(U) IF(FU.LT.FC)THEN AX=BX FA=FB BX=U FB=FU GO TO 1 C C...got a minimum between A and U, which will exit ELSE IF(FU.GT.FB)THEN CX=U FC=FU GO TO 1 ENDIF C C...parabolic fit was no use. Use default magnification U=CX+GOLD*(CX-BX) FU=FUNC(U) C C...parabolic fit is between C and its allowed limit ELSE IF((CX-U)*(U-ULIM).GT.0.d0)THEN FU=FUNC(U) IF(FU.LT.FC)THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU FU=FUNC(U) ENDIF C C...limit parabolic U to maximum allowed value ELSE IF((U-ULIM)*(ULIM-CX).GE.0.d0)THEN U=ULIM FU=FUNC(U) C C...reject parabolic U, use default magnification ELSE U=CX+GOLD*(CX-BX) FU=FUNC(U) ENDIF C C...eliminate oldest point and continue AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GO TO 1 ENDIF C RETURN END C C_____________________________________________________________________________ C FUNCTION BRENT(AX,BX,CX,F,TOL,XMIN) C C Given a function F, and given a bracketing triplet of abscissas AX, BX, C CX (such that BX is between AX and CX, and F(BX) is less than both F(AX) C and F(CX)), this routine isolates the minimum to a fractional precision C of about TOL using Brent's method. The abscissa of the minimum is C returned as XMIN, and the minimum function value is returned as BRENT, the C returned function value. C IMPLICIT REAL*8 (A-H,O-Z) ZK C C...maximum allowed number of iterations; golden ratio; and a small number C which protects against trying to achieve fractional accuracy for a minimum C that happens to be exactly zero. C PARAMETER (ITMAX=100,CGOLD=.3819660d0,ZEPS=1.0d-10) C C...A and B must be in ascending order, though the input abscissas need not be C A=DMIN1(AX,CX) B=DMAX1(AX,CX) C C...initialisations C V=BX W=V X=V E=0.d0 FX=F(X) FV=FX FW=FX C C...Main program loop C DO 11 ITER=1,ITMAX XM=0.5d0*(A+B) TOL1=TOL*DABS(X)+ZEPS TOL2=2.d0*TOL1 C C...test for done here IF(DABS(X-XM).LE.(TOL2-.5d0*(B-A))) GOTO 3 C C...construct a trial parabolic fit IF(DABS(E).GT.TOL1) THEN R=(X-W)*(FX-FV) Q=(X-V)*(FX-FW) P=(X-V)*Q-(X-W)*R Q=2.d0*(Q-R) IF(Q.GT.0.d0) P=-P Q=DABS(Q) ETEMP=E E=D IF(DABS(P).GE.DABS(.5d0*Q*ETEMP).OR.P.LE.Q*(A-X).OR. * P.GE.Q*(B-X)) GOTO 1 C C...the above conditions determine the acceptability of the parabolic fit. C Here it is OK: take a parabolic step D=P/Q U=X+D IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=DSIGN(TOL1,XM-X) C C...skip over the golden section step GOTO 2 ENDIF C C...we arrive here for a golden section step, which we take into the larger C of the two segments 1 IF(X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF C C...take the golden section step D=CGOLD*E C C...arrive here with D computed either from parabolic fit, or else from golden C section 2 IF(DABS(D).GE.TOL1) THEN U=X+D ELSE U=X+DSIGN(TOL1,D) ENDIF C C...this is the one function evaluation per iteration, and we have to decide C what to do with our function evaluation FU=F(U) C C...housekeeping follows: IF(FU.LE.FX) THEN IF(U.GE.X) THEN A=X ELSE B=X ENDIF V=W FV=FW W=X FW=FX X=U FX=FU ELSE IF(U.LT.X) THEN A=U ELSE B=U ENDIF IF(FU.LE.FW .OR. W.EQ.X) THEN V=W FV=FW W=U FW=FU ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN V=U FV=FU ENDIF ENDIF C C...done with housekeeping, back for another iteration 11 CONTINUE PAUSE 'Brent exceed maximum iterations.' C C...Arrive here ready to exit with best values C 3 XMIN=X BRENT=FX C RETURN END C C_____________________________________________________________________________ C FUNCTION F1DIM(X) C C This function must accompany LINMIN C IMPLICIT REAL*8 (A-H,O-Z) ZK PARAMETER (NMAX=6) COMMON /F1COM/ PCOM(NMAX),XICOM(NMAX),NCOM DIMENSION XT(NMAX) C DO 11 J=1,NCOM XT(J)=PCOM(J)+X*XICOM(J) 11 CONTINUE F1DIM=ESINTE(XT) ZK C RETURN END C_____________________________________________________________________________ C_____________________________________________________________________________