hdiff output

r32044/sqnm.f90 2017-03-07 20:36:28.472524484 +0000 r32043/sqnm.f90 2017-03-07 20:36:28.752528151 +0000
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/sqnm.f90' in revision 32043
  2: !   Copyright (C) 1999-2017 David J. Wales 
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: !                   SUBROUTINE SLBFGS 
 20: !Written by jk669 Feb 2017, based on code from mymylbfgs.f 
 21: ! 
 22: !   Stabilized L-BFGS algorithm (Schaefer et al) added 2017: 
 23: !   Stabilized quasi-Newton optimization of noisy potential energy surfaces 
 24: !        J. Chem. Phys. 142. 034112 (2015) doi:10.1063/1.4905665 
 25: !       
 26: !This subroutine is called from mylbfgs.f90 only. 
 27: !Requires LAPACK functions: dsyev, dsysv 
 28: !Dependencies: Requires calls to potential and topology_reader, and uses some global variables from Commons 
 29: ! 
 30: module sqnm_func  
 31:     implicit none 
 32: contains 
 33:     !dot product 
 34:     pure double precision function dot(numcoords, v1, v2) 
 35:         implicit none 
 36:         integer, intent(in) :: numcoords 
 37:         double precision, intent(in) :: v1(numcoords), v2(numcoords) 
 38:         double precision :: temp 
 39:         integer :: i 
 40:         temp=0.0d0 
 41:         do i=1, numcoords 
 42:             temp=temp+v1(i)*v2(i) 
 43:         end do  
 44:         dot=temp 
 45:     end function dot 
 46:  
 47:     !vector norm 
 48:     pure double precision function norm(numcoords, v1) 
 49:         implicit none 
 50:         integer, intent(in) :: numcoords 
 51:         double precision, intent(in) :: v1(numcoords) 
 52:         norm = sqrt(dot(numcoords,v1,v1)) 
 53:     end function norm 
 54: end module sqnm_func 
 55:  
 56: subroutine sqnm(numcoords,xcoords,maxrmsgrad,nhistmax,converget,energy,itermax,iter) 
 57:     use COMMONS, only: MAXERISE,SQNM_DEBUGT,AMBERT,AMBER12T,BONDS,SQNM_DEBUGRUN,SQNM_DEBUGLEVEL,SQNM_WRITEMAX,SQNM_BIOT,RMS 
 58:     use sqnm_func 
 59:     implicit none 
 60:     !arguments 
 61:     integer, intent(in)             :: numcoords  
 62:     double precision, intent(inout) :: xcoords(numcoords)  
 63:     double precision, intent(in)    :: maxrmsgrad 
 64:     integer, intent(in)             :: nhistmax  
 65:     logical, intent(inout)          :: converget  
 66:     double precision, intent(inout) :: energy  
 67:     integer, intent(in)             :: itermax  
 68:     integer, intent(inout)          :: iter 
 69:  
 70: !The subroutine is invoked as  
 71: !   slbfgs(numcoords,xcoords,maxrmsgrad,converget,energy,itermax,iter) 
 72: !where 
 73: !   numcoords is a positive integer equal to the number of coordinates. This should be equal to NATOMS*3. It is the responisibility 
 74: !           of the calling function to verify this.  
 75: !   xcoords is a vector of length numcoords containing the coordinates of the initial position on input, and containing the  
 76: !           optimized position upon exit if converget==.true., or the best positon obtained thus far if converget==.false. 
 77: !   maxrmsgrad is a cut-off for the gradient. The rms of the gradient is compared to maxrmsgrad to determine if the algorithm 
 78: !           has found a local minimum. The algorithm is converged if rms(grad)<maxrmsgrad. 
 79: !   nhistmax is a positive integer containing the maximum number of historic position and gradient vectors that will be stored for   
 80: !           determination of teh significant subspace. 
 81: !   converget is a logical which should be false on input (since the system is not converged) and is false on output UNLESS the  
 82: !           gradient at the new position is less than maxrmsgrad. 
 83: !   energy is a single double precision value that contains the energy of the function at the initial position on input, and the  
 84: !           energy of the optimal position on output.  
 85: !   itermax is a positive integer containing the maximum number of minimization steps that may be performed during the algorithm. 
 86: !   iter is the actual number of steps performed during the minimization algorithm. Note that any preliminary steepest descent 
 87: !           steps used to move the initial position closer to a minimum are NOT included in this count.  
 88: ! 
 89:     !**********parameters********** 
 90:     double precision, parameter :: alphascale=1.1d0 !preconditioned gradient scaling 
 91:     double precision, parameter :: minalpha=1.0d-15 !used to prevent alpha from getting too small, which may cause underflow 
 92:     double precision, parameter :: initsteepstepsize=1.0d-3 !initial steepest descent stepsize 
 93:     double precision, parameter :: cutoffratio=1.0d-4 !cutoff ratio for small eigenvalues for finding significant subspace 
 94:     double precision, parameter :: anglecutoff=0.2d0 !cos(angle) cutoff for changing preconditioning with factor alpha  
 95:     integer, parameter :: maxdescents=15 !maximum number of steepest descent steps allowed 
 96:     integer, parameter :: maxeeval=100 !maximum number potential calls allowed in a line search 
 97:     !**********local variables********** 
 98:     integer :: i, j !counters 
 99:     integer :: nhist=0; !current history length 
100:     double precision :: grad(numcoords) !current gradient 
101:     double precision :: pregrad(numcoords) !preconditioned gradient in subspace G 
102:     double precision :: gradhist(numcoords,nhistmax) !history of gradients 
103:     double precision :: xhist(numcoords,nhistmax) !history of gradients 
104:     integer :: lwork !used for LAPACK function dsyev 
105:     integer :: ndim !number of dimensions of the significant subspace.  
106:     integer, save :: nbonds !number of bonds in the molecule (Used for biomolecules) 
107:     integer :: eeval ! number of potential calls in the linesearch  
108:     double precision :: maxlambda !the value of 1/(maximum eigenvalue of the overlap matrix) 
109:     double precision :: tempval !exactly what you think. 
110:     double precision :: gradperp(numcoords) !portion of gradient perpendicular to preconditioned gradient 
111:     double precision :: newxcoords(numcoords) !newly calculated x-coordinates 
112:     double precision :: oldxcoords(numcoords) !previous x-coordinates 
113:     double precision :: oldgrad(numcoords) !previous gradient. 
114:     double precision :: newgrad(numcoords) ! gradient of energy at newly calculated x-coordinates 
115:     double precision :: residuetemp(numcoords) !temporary vector for calculating residues 
116:     double precision :: invproj(numcoords,numcoords) !inverse projection from subspace to find orthogonal complement 
117:     double precision :: enew !energy at newly calculated x-coordinates 
118:     double precision :: ethresh !Energy threshhold to determine if a step is accepted.  
119:     double precision :: alpha !preconditioned gradient scaling 
120:     double precision :: angle !angle between gradient and preconditioned gradient 
121:     double precision :: gradrms !Root mean squared of the gradient 
122:     double precision :: steepstepsize !steepest descent step size 
123:     logical :: acceptstep !Set based on whether the last step was accepted or not. 
124:     logical :: bio = .false. !Whether the molecule being minimized is a biomolecule (true) or not (false). 
125:     double precision :: gradstr(numcoords) !the bond-stretching component of the gradient for biomolecules 
126:     integer :: atom1, atom2 !the atom numbers used for finding atom displacements for biomolecules  
127:     double precision :: atom_disp(3) !Displacement between atom1 and atom2 for biomolecules 
128:     integer :: signflip !number of projections which have flipped since the last iteration for biomolecules 
129:     double precision :: alpha_bio !steepest descent stepsize for biomolecule gradients  
130:     integer :: infoc !information integer from subroutines and LAPACK functions.  
131:     integer :: totalfunceval ! number of calls to the potential energy function. 
132:     integer, save :: runs=0 !keeps track of how many times this subroutine is called during the program. 
133:     logical :: debug_run !If set to true, this prints debug information for the run.  
134:     integer :: writemax !maximum number of lines of something to pring when debugging 
135:  
136:     !allocatable local variables. All variables here have at least one dimension of length nhist, ndim, or nbond, which  
137:     !change or are unknown. 
138:     double precision, allocatable, dimension(:,:) :: bondvect !sparce vectors containing atom distances of bonded atoms. 
139:     double precision, allocatable, dimension(:)   :: oldcoeff !sparce vectors containing atom distances of bonded atoms.   
140:     double precision, allocatable, dimension(:)   :: bondcoeff !vector with coefficients of bond vectors to find gradstr 
141:     double precision, allocatable, dimension(:,:) :: bondmat !Matrix containging bondvect^t bondvect 
142:     double precision, allocatable, dimension(:)   :: ipiv !vector of length nbonds for dsysv 
143:     double precision, allocatable, dimension(:,:) :: graddisp !gradient displacement vectors 
144:     double precision, allocatable, dimension(:,:) :: xdisp !xcoord displacement vectors 
145:     double precision, allocatable, dimension(:)   :: norms !the norms of the xdisp vectors 
146:     double precision, allocatable, dimension(:,:) :: overlap !Overlap matrix of normalized displacement vectors 
147:     double precision, allocatable, dimension(:)   :: work !used for LAPACK function dsyev 
148:     double precision, allocatable, dimension(:)   :: eigenvalues !contains eigenvalues of overlap matrix or projected hessian 
149:     double precision, allocatable, dimension(:,:) :: subspan !matrix with subspace-spanning vectors as columns 
150:     double precision, allocatable, dimension(:,:) :: gradspan !matrix with curvature-describing vectors as columns 
151:     double precision, allocatable, dimension(:,:) :: sigsubspan !matrix with significant subspace-spanning vectors as columns 
152:     double precision, allocatable, dimension(:,:) :: siggradspan !matrix with significant curvature-describing vectors as columns 
153:     double precision, allocatable, dimension(:,:) :: SO !approximation of shape operator 
154:     double precision, allocatable, dimension(:,:) :: principals !holds principal directions of the projected hessian 
155:     double precision, allocatable, dimension(:)   :: residues !contains residues calculated by Weinstein's criterion 
156:  
157:     !increment number of runs to keep track of calls to the function during GMIN. 
158:     runs=runs+1 
159:  
160:     !Determine if run is being debugged 
161:     debug_run=.false. 
162:     if(SQNM_DEBUGT .and. (SQNM_DEBUGRUN==0 .or. SQNM_DEBUGRUN+1==runs)) then  
163:         debug_run=.true. 
164:         writemax=SQNM_WRITEMAX 
165:     end if  
166:     if((runs > SQNM_DEBUGRUN+1) .and. SQNM_DEBUGT ) STOP !stop after debug info is printed 
167:  
168:     !check input for any potential issues (responisbility for correctness here lies with calling subroutine).  
169:     if (numcoords<=0 .or.itermax<=0 .or. maxrmsgrad<=0.0d0 ) then  
170:         if(debug_run .and. SQNM_DEBUGLEVEL>=0) write(*, '(a,i5,a)') & 
171:             'WARNING: subroutine slbfgs in slbfgs.f90 received bad input on run ', runs, & 
172:             '. No minimization occurred.' 
173:         return 
174:     end if 
175:  
176:     if(debug_run .and. SQNM_DEBUGLEVEL>=0) write(*, '(a,i5,a)') 'STARTING RUN :', runs, '.' 
177:  
178:     !initialize variables. 
179:     converget=.false. 
180:     iter=1 
181:     nhist=0 
182:     ethresh=MAXERISE !from COMMONS 
183:     lwork=4*numcoords+64 !size of work vector for LAPACK function dsyev 
184:     totalfunceval=0 
185:     call potential(xcoords,grad,energy,numcoords,.FALSE.) 
186:     allocate(work(lwork)) 
187:     steepstepsize=min(initsteepstepsize,1.0d0/maxval(abs(grad))) 
188:     totalfunceval=totalfunceval+1 !update number of energy function calls 
189:     if (debug_run .and. SQNM_DEBUGLEVEL>=0) then  
190:         write(*, '(A)') 'xcoords and grad before steepest descent steps' 
191:         do i=1,min(writemax,numcoords) 
192:             write(*, '(2(f20.15))') xcoords(i), grad(i) 
193:         end do  
194:     end if 
195:     !In GMIN, atoms are moved randomly, and initial gradients may be very large. Do steepest descent steps to reduce energy. 
196:     !perform at least one steepest descent step, and store more recent previous position and gradient in the histories. 
197:     !Steepest descent is performed with a constant stepsize to prevent numerous function calls for linesearch.  
198:     steepest_descent: do while (iter<maxdescents) 
199:         if (debug_run .and. SQNM_DEBUGLEVEL>=1) then 
200:             print *, 'Iter: ', iter, ' Init E: ', energy, '. Init gradrms: ', norm(numcoords, grad), 'stepsize', steepstepsize  
201:         end if 
202:         eeval=0 
203:         ! store position and gradient in history before linsearch, since this subroutine changes values.  
204:         xhist(:,1)=xcoords(:) 
205:         gradhist(:,1)=grad(:) 
206:         !update the stepsize. Try a larger step if possible, or way larger if step is very small. 
207:         if (steepstepsize<1.0d-8 .and. iter>=3) then  
208:             steepstepsize=5.0D-2 
209:         else  
210:             steepstepsize=alphascale*steepstepsize  
211:         end if  
212:         tempval=energy !store energy value to ensure it's lower 
213:         call linesearch(numcoords,xcoords,energy,grad,-grad,steepstepsize,ethresh,maxeeval,infoc,eeval) 
214:         totalfunceval=totalfunceval+eeval !update number of energy function calls 
215:         gradrms=norm(numcoords,grad) 
216:         if (RMS<maxrmsgrad) then !we already converged! Great work, everyone. 
217:             converget=.true. 
218:             if(debug_run .and. SQNM_DEBUGLEVEL>=0 ) then 
219:                 write(*, '(a, i4,a,i4)') 'Run ', runs, ' has converged during steepest descent, iteration ', iter 
220:                 write(*,'(a,f12.8,a,f12.6)') 'Final rms of gradient = ', gradrms, ', Final energy: ', energy 
221:                 print *, 'Total potential energy function evaluations: ', totalfunceval 
222:             end if 
223:             return 
224:         else if (tempval+ethresh<energy.or.(gradrms>1.0d10.and.iter>=itermax/2).or.(gradrms>1.0d15.and.iter>=2)) then  
225:             !bad things have happened; the energy has increased after the stepsearch. 
226:             !the gradient was probably way too massive coming in, and a small enough step could not be taken. 
227:             !Escape to prevent overflow.  
228:             if(debug_run .and. SQNM_DEBUGLEVEL>=0 ) then 
229:                 print *, 'ERROR: Run did not converge due to stepsearch being unable to locate a minimum. ' 
230:                 write(*, '(a, i4,a,i4)') ' Stepsearch terminated on run ', runs, ', stepsearch iteration ', iter 
231:                 write(*,'(a,f12.8,a,f12.6)') 'RMS of gradient = ', gradrms, ', Energy: ', energy 
232:                 print *, 'Total potential energy function evaluations: ', totalfunceval 
233:             end if 
234:             return 
235:         end if 
236:         iter=iter+1 
237:     end do steepest_descent 
238:  
239:     acceptstep=.true. !since we did a steepest descent step, we had a successful step, and it is stored in history 
240:     nhist=1 !since previous position stored in steepest descent portion 
241:     alpha=min(initsteepstepsize,1.0d0/maxval(abs(grad))) !works well emperically for LJ clusters and oligopeptides 
242:     if (debug_run .and. SQNM_DEBUGLEVEL>=1) then 
243:         write(*, '(a, i4, a, f15.8, a, f15.8)') 'Number of steepest descent steps taken: ', iter, &  
244:             '. Energy after Steepest Descents: ', energy, '. Gradient RMS after  Steepest Descents: ', gradrms 
245:     end if 
246:     if (converget) then 
247:         if(debug_run .and. SQNM_DEBUGLEVEL>=0) write(*, '(a,i4,a,i3,a)') 'Run ', runs, ' converged after ', iter,& 
248:          ' steps. No slbfgs needed.' 
249:         return 
250:     end if  
251:     if (debug_run .and. SQNM_DEBUGLEVEL>=1) then  
252:         write(*, '(A)') 'xcoords and grad after steepest descent steps' 
253:         do i=1,min(writemax,numcoords) 
254:             write(*, '(2(f20.15))') xcoords(i), grad(i) 
255:         end do  
256:     end if 
257:     !if we are using an AMBER force field, we will assume that we are dealing with a biomolecule, and use the  
258:     !procedure for biomolecules. For this, we need the bond information. Also make sure that the bio portion was not  
259:     ! turned off.  
260:     biocheck: if((AMBERT .or. AMBER12T) .and. (SQNM_BIOT)) then  
261:         !this function call sets nbonds to the number of bonds, and then the global variable BONDS  
262:         !is a nbonds*2 integer array containing pairs of bonded atoms.  
263:         bio = .true. 
264:         if (.not. allocated(BONDS)) then  
265:             call topology_reader(nbonds) 
266:         end if  
267:         if (.not. allocated(bondvect)) then 
268:             allocate(bondvect(numcoords,nbonds)) 
269:             allocate(bondcoeff(nbonds)) 
270:             allocate(oldcoeff(nbonds)) 
271:             allocate(bondmat(nbonds,nbonds)) 
272:         end if  
273:         bondcoeff(:)=0.0d0 
274:         !set bond steepest descent stepsize.  
275:         alpha_bio = steepstepsize  
276:         if (debug_run .and. SQNM_DEBUGLEVEL>=1) then  
277:             write(*,'(a,i5,a)') 'The molecule has been identified as a biomolecule. ', nbonds, ' bonds were found.' 
278:             write(*,'(a,i5,a)') 'The molecule has has ', numcoords/3, ' atoms.' 
279:         end if  
280:     else if ((AMBERT .or. AMBER12T) .and. .not. (SQNM_BIOT)) then biocheck 
281:         if (debug_run .and. SQNM_DEBUGLEVEL>=1) then  
282:             write(*,'(2a)') 'The minimization of this molecule is using the AMBER potential but is not being treated as a ', & 
283:             'biomolcule due to the presence of the BIOOFF keyword in the data file.' 
284:         end if  
285:     end if biocheck 
286:  
287: !******************************************************************************************************** 
288: !Main slbfgs do loop. This is broken down into numerous parts:                                          * 
289: !    If we have a biomolecule,                                                                          * 
290: !       Compute bond stretching components.                                                             * 
291: !       Find bond vector coefficients                                                                   * 
292: !       Adjust alpha_bio                                                                                * 
293: !       Adjust gradient and position based on the stretching components.                                * 
294: !    Find preconditioned gradient.                                                                      * 
295: !       find x-displacements and grad displacements from histories.                                     * 
296: !           find x-displacement norms and normalize x-displacement vectors.                             * 
297: !       find overlap matrix.                                                                            * 
298: !           get eigenvectors and eigenvalues of overlap matrix                                          * 
299: !           find number of dimensions of significant subspace based on eigenvalues                      * 
300: !       find significant subspace of tangent space on energy manifold at position                       * 
301: !           calculate spanning vectors of subspace                                                      * 
302: !           determine significance of vectors from overlap eigenvalues                                  * 
303: !       obtain curvature information from significant subspace                                          * 
304: !           find 'normalized' gradients in tangent subspace                                             * 
305: !           produce symmetrized shape operator                                                          * 
306: !           project principal directions of shape operator onto energy manifold                         * 
307: !       correct curvatures due to deviations of calculated shape operator from second fundamental form  * 
308: !           calculate Weinstein residues                                                                * 
309: !           adjust curvatures                                                                           * 
310: !       precondition the gradient                                                                       * 
311: !           find rejection of gradient from projection onto preconditioned gradient vector              * 
312: !           adjust rejection by factor of alpha                                                         * 
313: !    x(iter+1)=x-preconditioned gradient.                                                               * 
314: !    Find energy at x(iter+1), and determine if minimization occurred.                                  * 
315: !       If energy decreased or alpha < alphastart/10, accept step                                       * 
316: !           update x-coords and gradient histories                                                      * 
317: !           if nhist<nhistmax, nhist++                                                                  * 
318: !           adjust alpha                                                                                * 
319: !       Else, reject step (x(iter+1)=x(iter)                                                            * 
320: !           alpha=alpha/2                                                                               * 
321: !           purge history                                                                               * 
322: !   iter++                                                                                              * 
323: !Loop until converged or maxiter is reached.                                                            * 
324: !******************************************************************************************************** 
325:     mainloop: do while (.not. converget .and. iter .le. itermax) 
326:         if(debug_run .and. SQNM_DEBUGLEVEL>=0) then  
327:             write(*, '(a, i4, a, f15.8, a, f15.8)') 'Number of iterations done: ', iter, &  
328:             '. Energy: ', energy, '. Gradient RMS: ', gradrms 
329:             print *, 'Bio On:',SQNM_BIOT,'alpha_bio: ', alpha_bio, ' bio: ',bio,' acceptstep: ',acceptstep,' nhist : ',nhist 
330:             print *, 'run: ',runs,' iter: ',iter,' nhist: ',nhist,' Energy: ',energy,' RMSgrad: ',gradrms,' alpha: ',alpha  
331:             write(*, '(A)') 'xcoords and grad at loop start' 
332:             do i=1,min(writemax,numcoords) 
333:                 write(*, '(2(f24.20))') xcoords(i), grad(i) 
334:             end do  
335:         end if  
336:  
337:         oldgrad(:)=grad(:) 
338:         oldxcoords(:)=xcoords(:) 
339:         !FOR BIOMOLECULES, WE MUST FIRST FIND BOND SRETCHING COMPONENTS 
340:         bio_optim: if (bio .and. acceptstep) then 
341:             !save old positions in case new ones are higher in energy.  
342:             !FIND ATOM DISTANCES FOR BOND VECTORS  
343:             do i=1,nbonds 
344:                 !get bonded atom distances. 
345:                 !bondvect is sparce with six elements in each column, corresponding to the bonded atom displacements.  
346:                 !Because displacement is antisymmetric, each column sum should be 0. 
347:                 atom1=BONDS(i,1) 
348:                 atom2=BONDS(i,2) 
349:                 !displacement is xyz(atom2)-xyz(atom1) 
350:                 atom_disp(:)=xcoords(3*atom2-2:3*atom2)-xcoords(3*atom1-2:3*atom1) 
351:                 bondvect(3*atom1-2:3*atom1,i)=atom_disp(:) 
352:                 bondvect(3*atom2-2:3*atom2,i)=-atom_disp(:) !by antisymmetry 
353:             end do 
354:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then 
355:                 print *, 'bond vectors' 
356:                 do i=1,min(writemax,nbonds)  
357:                     write(*, '(*(f13.9))') (bondvect(i,j),j=1,min(writemax,nbonds)) 
358:                 end do 
359:             end if  
360:             !FIND BOND VECTOR COEFFICIENTS  
361:             !This is done by solving (bondvect^t * bondvect)*coeffs = bondvect^t * grad 
362:             !find bondgrad=bondvect^t * grad and bondmat=(bondvect^t * bondvect) 
363:             do i=1,nbonds 
364:                 bondcoeff(i)=dot(numcoords,bondvect(:,i),grad(:)) 
365:                 bondmat(i,i)=dot(numcoords,bondvect(:,i),bondvect(:,i)) 
366:                 do j=i+1,nbonds 
367:                     !compute only lower half since matrix is symmetric. 
368:                     bondmat(j,i)=dot(numcoords,bondvect(:,i),bondvect(:,j)) 
369:                 end do 
370:             end do 
371:             !count number of projections bondcoeff which have not changed sign from last time by comparing to the old coefficient 
372:             signflip=0 
373:             do i=1,nbonds 
374:                 if ((oldcoeff(i)>=0.0d0 .and. bondcoeff(i)>=0.0d0) .or. (oldcoeff(i)<0.0d0 .and. bondcoeff(i)<0.0d0)) then 
375:                     signflip=signflip+1 
376:                 end if  
377:             end do 
378:             oldcoeff(:)=bondcoeff(:) 
379:  
380:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then 
381:                 print *, 'bond matrix' 
382:                 do i=1,min(writemax,nbonds)  
383:                     write(*, '(*(f13.9))') (bondmat(i,j),j=1,min(writemax,nbonds)) 
384:                 end do 
385:                 print *, 'bond coeffs' 
386:                 do i=1,min(writemax,nbonds) 
387:                     write(*,'(f16.10)') bondcoeff(i) 
388:                 end do 
389:                 print *, 'Number of sign flips: ', signflip, '; alpha_bio = ', alpha_bio 
390:             end if  
391:  
392:  
393:             !Solve linear system Bondmat*coeff=bondvect.  
394:             !Note that a LAPACK function for symmetric matrices is called since we calculated the lower 
395:             !triangular portion of bondmat. Solution should now be in bondcoeff.  
396:             allocate(ipiv(nbonds)) 
397:             call dsysv('L',nbonds,1,bondmat,nbonds,ipiv,bondcoeff,nbonds,work,lwork,infoc) 
398:             deallocate(ipiv) 
399:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then 
400:                 print *, 'bond coeffs AFTER solving' 
401:                 do i=1,min(writemax,nbonds)  
402:                     write(*,'(f16.10)') bondcoeff(i) 
403:                 end do 
404:             end if  
405:  
406:  
407:             !FIND BOND STRETCHING COMPONENTS OF THE GRDAIENT, gradstr 
408:             gradstr(:)=0 
409:             do i=1,nbonds 
410:                 gradstr(:)=gradstr(:)+bondcoeff(i)*bondvect(:,i) 
411:             end do 
412:  
413:             !MINIMIZE ALONG BOND-STRETCHING COMPONENTS AND SEND REST OF GREADIENT TO PRECONDITIONING 
414:             !adjust gradient to non-bonded portions  
415:             grad(:)=grad(:)-gradstr(:) 
416:             !scale bonded gradient component steepest descent stepsize based on number of sign flips  
417:             if(signflip > 2*nbonds/3) then  
418:                 alpha_bio=min(alpha_bio*alphascale,1.0d0) 
419:             else 
420:                 alpha_bio=alpha_bio/alphascale 
421:             end if 
422:             !perform simple steepest descent step along bonded portions of gradient  
423:             xcoords(:)=xcoords(:)-alpha_bio*gradstr(:) 
424:             if (debug_run .and. SQNM_DEBUGLEVEL>=1) then 
425:                 print *, 'New alpha_bio: ', alpha_bio 
426:                 write(*, '(A)') '   old xcoords         xcoords             old grad            new grad            grad_str' 
427:                 do i=1,min(writemax,numcoords) 
428:                     write(*, '(5(f20.15))') oldxcoords(i), xcoords(i), oldgrad(i), grad(i), gradstr(i) 
429:                 end do  
430:             end if 
431:         end if bio_optim 
432:  
433:  
434:         !**********FINDING PRECONDITIONED GRADIENT********** 
435:         !FIND COORDINATE AND GRADIANT DISPLACEMENTS   
436:         preconditioning: if (nhist==0) then !last step failed; perform a steepest descent-like step 
437:             pregrad(:)=alpha*grad 
438:         else preconditioning 
439:             !FIND POSITION AND GRADIENT FINITE DIFFERENCES 
440:             allocate(xdisp(numcoords,nhist)) 
441:             allocate(graddisp(numcoords,nhist))  
442:             !for first column of displacements, use current position and then first historical position 
443:             xdisp(:,1)=xcoords(:)-xhist(:,1) 
444:             graddisp(:,1)=grad(:)-gradhist(:,1) 
445:             !for remaining columns, just use historical positions. 
446:             do i=2,nhist 
447:                 xdisp(:,i)=xhist(:,i-1)-xhist(:,i) 
448:                 graddisp(:,i)=gradhist(:,i-1)-gradhist(:,i) 
449:             end do   
450:  
451:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then  
452:                 write(*, '(A)') 'xhist' 
453:                 do i=1,min(writemax,numcoords) 
454:                     write(*, '(*(f20.15))') (xhist(i,j),j=1,nhist) 
455:                 end do  
456:                 write(*, '(A)') 'gradhist' 
457:                 do i=1,min(writemax,numcoords) 
458:                     write(*, '(*(f20.15))') (gradhist(i,j),j=1,nhist) 
459:                 end do  
460:                 write(*, '(A)') 'xdisp' 
461:                 do i=1,min(writemax,numcoords) 
462:                     write(*, '(*(f16.10))') (xdisp(i,j),j=1,nhist) 
463:                 end do 
464:                 write(*, '(A)') 'graddisp' 
465:                 do i=1,min(writemax,numcoords) 
466:                     write(*, '(*(f16.10))') (graddisp(i,j),j=1,nhist) 
467:                 end do    
468:             end if  
469:  
470:             !FIND NORMS AND NORMALIZE X-DISPLACEMENT VECTORS 
471:             allocate(norms(nhist)) 
472:             do i=1, nhist  
473:                 norms(i)=norm(numcoords,xdisp(:,i)) 
474:                 if (norms(i)<1.0d-100) then !floating point errors caused underflow. Return to avoid division by 0. 
475:                     if (debug_run  .and. SQNM_DEBUGLEVEL>=0) print *,'Norm too small for division on run ',runs,', iter,',iter,& 
476:                         'hist ', nhist, ' Bombing out...' 
477:                     converget=.false. 
478:                     return 
479:                 end if 
480:                 !non-normalized x-displacements will not be needed again, so we can normalize in original matrix. 
481:                 xdisp(:,i)=xdisp(:,i)/norms(i) 
482:             end do 
483:  
484:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then  
485:                 write(*, '(A)') 'displacement norms' 
486:                 do i=1,min(writemax,nhist) 
487:                     write(*, '(*(f16.10))') (norms(i)) 
488:                 end do  
489:                 write(*, '(A)') 'xdisp - normalized' 
490:                 do i=1,min(writemax,numcoords) 
491:                     write(*, '(*(f16.10))') (xdisp(i,j),j=1,min(writemax,nhist)) 
492:                 end do   
493:             end if 
494:  
495:             !FIND OVERLAP MATRIX 
496:             allocate(overlap(nhist, nhist)) 
497:             do i=1, nhist 
498:                 !diagonal elements 
499:                 overlap(i,i)=dot(numcoords,xdisp(:,i),xdisp(:,i)) 
500:                 do j=i+1, nhist 
501:                     !compute only lower half since matrix is symmetric. 
502:                     overlap(j,i)=dot(numcoords,xdisp(:,i),xdisp(:,j)) 
503:                 end do 
504:             end do 
505:  
506:             if (debug_run .and. SQNM_DEBUGLEVEL>=3) then  
507:                 write(*, '(A)') 'overlap' 
508:                 do i=1,nhist 
509:                     write(*, '(*(f16.10))') (overlap(i,j),j=1,nhist) 
510:                 end do 
511:             end if 
512:  
513:             !GET EIGENVALUES AND EIGENVECTORS OF OVERLAP MATRIX 
514:             !allocate vectors for eigenvectors and for work for dsyev 
515:             allocate(eigenvalues(nhist)) 
516:             !dsyev will place eigenvalues of overlap into the eigenvalues vector in order of ascending value,  
517:             !and the columns of overlap will be replaces with the corresponding eigenvectors. 
518:             !Note that the LAPACK function for symmetric matrices is called since we calculated the lower 
519:             !triangular portion of overlap. 
520:             call dsyev('V',"L",nhist,overlap,nhist,eigenvalues,work,lwork,infoc)  
521:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then  
522:                 write(*, '(A)') 'eigenvalues' 
523:                 do i=1,nhist 
524:                     write(*, '(*(f16.10))') (eigenvalues(i)) 
525:                 end do 
526:                 write(*, '(A)') 'overlap matrix eigenvectors' 
527:                 do i=1,nhist 
528:                     write(*, '(*(f16.10))') (overlap(i,j),j=1,nhist) 
529:                 end do    
530:             end if  
531:  
532:             !FIND THE DIMENSIONS OF THE SIGNIFICANT SUBSPACE. This is done by comparing the ratio of all eigenvalues to  
533:             !the largest eigenvalue. This is used to find the dimension of the significant subspace: 
534:             ndim=1 !minimum dimension if nhist>0 
535:             maxlambda=1.0D0/eigenvalues(nhist) !eigenvalues are in ascending order, so eigenvalues(nhist) is the largest. 
536:             do i=1, nhist-1 
537:                 if (maxlambda*eigenvalues(i)>=cutoffratio) then !equivalent to checking eigenvalues(i)/eigenvalues(nhist) 
538:                     ndim=ndim+1 
539:                 end if  
540:             end do 
541:  
542:             !FIND SPANNING VECTORS OF SUBSPACE. ALSO FIND NORMALIZED GRADIENTS 
543:             allocate(subspan(numcoords,nhist)) 
544:             allocate(gradspan(numcoords,nhist)) 
545:             do i = 1, nhist 
546:                 subspan(:,i)=0.0D0 
547:                 gradspan(:,i)=0.0D0 
548:                 do j=1, nhist 
549:                     subspan(:,i)=subspan(:,i)+overlap(j,i)*xdisp(:,j) 
550:                     gradspan(:,i)=gradspan(:,i)+overlap(j,i)*graddisp(:,j)/norms(j) 
551:                 end do 
552:                 !note that any eigenvalues should be at least zero since the overlap matrix is positive semi-definite;  
553:                 !however, due to floating point error, we may have very slightly negative or zero eigenvalues. Take max().  
554:                 tempval=1.0D0/sqrt(max(eigenvalues(i),1.0d-10)) 
555:                 subspan(:,i)=subspan(:,i)*tempval 
556:                 gradspan(:,i)=gradspan(:,i)*tempval 
557:             end do 
558:  
559:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then  
560:                 write(*, '(A)') 'subspance span' 
561:                 do i=1,min(writemax,numcoords) 
562:                     write(*, '(*(f16.10))') (subspan(i,j),j=1,min(writemax,nhist)) 
563:                 end do  
564:                 write(*, '(A)') 'subspance gradients' 
565:                 do i=1,min(writemax,numcoords) 
566:                     write(*, '(*(f16.10))') (gradspan(i,j),j=1,min(writemax,nhist)) 
567:                 end do  
568:             end if  
569:  
570:             !SELECT VECTORS AND GRADIENTS IN THE SIGNIFICANT SUBSPACE 
571:             !there will be ndim of these vectors, out of nhist, with ndim<=nhist 
572:             allocate(sigsubspan(numcoords,ndim)) 
573:             allocate(siggradspan(numcoords,ndim)) 
574:             !Note that eigenvalues were listed in ascending order before, so we reverse the direction of vector selection. 
575:             !The first vectors correspond with the largest eigenvalue and so on. 
576:             do i=1, ndim 
577:                 sigsubspan(:,i)=subspan(:,nhist-i+1) 
578:                 siggradspan(:,i)=gradspan(:,nhist-i+1) 
579:             end do 
580:  
581:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then  
582:                 write(*, '(A)') 'significant subspance span' 
583:                 do i=1,min(writemax,numcoords) 
584:                     write(*, '(*(f16.10))') (sigsubspan(i,j),j=1,min(writemax,ndim)) 
585:                 end do  
586:                 write(*, '(A)') 'significant subspance gradients' 
587:                 do i=1,min(writemax,numcoords) 
588:                     write(*, '(*(f16.10))') (siggradspan(i,j),j=1,min(writemax,ndim)) 
589:                 end do  
590:             end if  
591:  
592:             !PRODUCE SYMMETRIC SHAPE OPERATOR 
593:             allocate(SO(ndim,ndim)) 
594:             SO(:,:)=0.0D0 
595:             do i=1, ndim 
596:                 SO(i,i)=dot(numcoords,siggradspan(:,i),sigsubspan(:,i)) !diagonals of shape operator 
597:                 do j=i+1, ndim 
598:                     !matrix is symmetric; find only lower triangular  
599:                     SO(j,i)=0.5D0*(dot(numcoords,siggradspan(:,i),sigsubspan(:,j))+ & 
600:                         dot(numcoords,siggradspan(:,j),sigsubspan(:,i)))  
601:                 end do 
602:             end do 
603:  
604:             if (debug_run .and. SQNM_DEBUGLEVEL>=5) then  
605:                 write(*, '(A)') 'shape operator matrix' 
606:                 do i=1,min(writemax,ndim) 
607:                     write(*, '(*(f16.10))') (SO(i,j),j=1,min(writemax,ndim)) 
608:                 end do    
609:             end if  
610:  
611:             !PERFORM EIGENVALUE DECOMFC=ifort cmake -DCMAKE_BUILD_TYPE=Debug -DWITH_AMBER=1 ../../sourcePOSITION OF SHAPE OPERATOR 
612:             !we no longer need the eigenvectors from the overlap matrix, so we reuse the vector. 
613:             deallocate(eigenvalues) 
614:             allocate(eigenvalues(ndim)) 
615:             !Note that the LAPACK function for symmetric matrices is called since we calculated the lower 
616:             !triangular portion of the shape operator. 
617:             call dsyev('V',"L",ndim,SO,ndim,eigenvalues,work,lwork,infoc)   
618:  
619:             if (debug_run .and. SQNM_DEBUGLEVEL>=5) then  
620:                 write(*, '(A)') 'shape operator eigenvalues' 
621:                 do i=1,min(writemax,ndim) 
622:                     write(*, '(f16.10)') (eigenvalues(i)) 
623:                 end do 
624:                 write(*, '(A)') 'shape operator matrix eigenvectors' 
625:                 do i=1,min(writemax,ndim) 
626:                     write(*, '(*(f16.10))') (SO(i,j),j=1,min(writemax,ndim)) 
627:                 end do    
628:             end if  
629:  
630:             !FIND PROJECTED PRINCIPAL DIRECTIONS 
631:             allocate(principals(numcoords,ndim)) 
632:             do i=1, ndim 
633:                 principals(:,i)=0.0D0 
634:                 do j=1, ndim 
635:                     principals(:,i)=principals(:,i)+SO(j,i)*sigsubspan(:,j) 
636:                 end do 
637:             end do  
638:  
639:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then  
640:                 write(*, '(A)') 'principals' 
641:                 do i=1,min(writemax,numcoords) 
642:                     write(*, '(*(f16.10))') (principals(i,j),j=1,min(writemax,ndim)) 
643:                 end do   
644:             end if 
645:  
646:             !ADJUST CURVATURES BY FINDING RESIDUES 
647:             allocate(residues(ndim)) 
648:             !calculate residues 
649:             do i=1, ndim 
650:                 residues(i)=0.0D0 
651:                 residuetemp(:)=0.0D0 
652:                 do j =1, ndim 
653:                     residuetemp(:)=residuetemp(:)+SO(j,i)*siggradspan(:,j) 
654:                 end do  
655:                 residuetemp(:)=residuetemp(:)-eigenvalues(i)*principals(:,i) 
656:                 residues(i)=norm(numcoords,residuetemp) 
657:             end do  
658:             !adjust egienvalues 
659:             do i=1,ndim 
660:                 eigenvalues(i)=sqrt(eigenvalues(i)**2+residues(i)**2) 
661:             end do 
662:  
663:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then  
664:                 write(*, '(2A)') 'residues     ', ' Adjusted eigenvalues' 
665:                 do i=1,min(writemax,ndim) 
666:                     write(*, '(i3,a,f12.5,A,g15.7)') i,': ', residues(i), ' ', eigenvalues(i) 
667:                 end do  
668:             end if 
669:  
670:  
671:             !FORM SUBSPACE PORTION OF PRECONDITIONED GRADIENT 
672:             pregrad(:)=0.0D0 
673:             do i=1, ndim 
674:                 tempval=dot(numcoords,grad,principals(:,i))/eigenvalues(i) 
675:                 pregrad(:)=pregrad(:)+tempval*principals(:,i) 
676:             end do 
677:             if (debug_run .and. SQNM_DEBUGLEVEL>=3) then  
678:                 write(*, '(A)') 'Preconditioned Gradient in Subspace' 
679:                 do i=1,min(writemax,numcoords) 
680:                     write(*, '(*(f16.10))') pregrad(i) 
681:                 end do 
682:             end if 
683:             !FIND COS(ANGLE) BETWEEN GRADIENT AND PRECONDITIONED GRADIENT TO ADJUST ALPHA LATER. 
684:             !ALSO FIND ORTHOGONAL COMPLEMENT OF GRADIENT: gradperp=(I-P')grad; invproj=(I-P') 
685:             angle=dot(numcoords,pregrad,grad)/(norm(numcoords,pregrad)*norm(numcoords,grad)) 
686:             !Note P'=principals*principals^T; can use dot products for each entry of invproj; add 1 to diagonals to get (I-P') 
687:             do i=1, numcoords 
688:                 invproj(i,i)=1.0d0-dot(ndim,principals(i,:),principals(i,:)) 
689:                 do j=i+1, numcoords 
690:                     invproj(j,i)=-dot(ndim,principals(i,:),principals(j,:)) 
691:                     invproj(i,j)=invproj(j,i) 
692:                 end do 
693:             end do 
694:  
695:             !Find rejection of the gradient projected onto the significant subspace  
696:             do i=1, numcoords 
697:                 gradperp(i)=0.0d0 
698:                 do j=1,numcoords 
699:                 gradperp(i)=gradperp(i)+invproj(j,i)*grad(j) 
700:                 end do 
701:             end do 
702:             pregrad(:)=pregrad(:)+alpha*gradperp(:) !full preconditioned gradient 
703:  
704:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then  
705:                 print *, 'alpha: ', alpha 
706:                 print *, 'Preconditioned Gradient, ', 'grad, ', 'gradperp' 
707:                 do i=1,min(writemax,numcoords) 
708:                     print *, pregrad(i), grad(i), gradperp(i) 
709:                 end do 
710:             end if 
711:  
712:             !Iteration Cleanup 
713:             deallocate(xdisp) 
714:             deallocate(graddisp) 
715:             deallocate(overlap) 
716:             deallocate(subspan) 
717:             deallocate(gradspan) 
718:             deallocate(norms) 
719:             deallocate(eigenvalues) 
720:             deallocate(residues) 
721:             deallocate(siggradspan) 
722:             deallocate(sigsubspan) 
723:             deallocate(SO) 
724:             deallocate(principals) 
725:         end if preconditioning 
726:         !**********DONE FINDING PRECONDITIONED GRADIENT********** 
727:  
728:         if (nhist==0) then !step was a pure steepest step. Do a line search. 
729:             newxcoords(:)=xcoords(:) 
730:             newgrad(:)=grad(:) 
731:             enew=energy 
732:             !try making alpha bigger 
733:             alpha=max(1.0d-3,alpha) 
734:             eeval=0 
735:             call linesearch(numcoords,newxcoords,enew,newgrad,-grad,alpha,ethresh,maxeeval,infoc,eeval) 
736:             totalfunceval=totalfunceval+eeval !update number of energy function calls 
737:         else 
738:             newxcoords(:)=xcoords(:)-pregrad(:) 
739:             call potential(newxcoords,newgrad,enew,numcoords,.FALSE.) 
740:             totalfunceval=totalfunceval+1 !update number of energy function calls 
741:         end if 
742:  
743:         !**********CHECK FOR STEP ACCEPTANCE********** 
744:         accept_check: if(enew>energy+ethresh) then 
745:             acceptstep=.false. !step rejected 
746:             if (bio) then  
747:                 grad(:)=oldgrad(:) 
748:                 xcoords(:)=oldxcoords(:) 
749:             end if  
750:             nhist=0 
751:             if (alpha>steepstepsize/1.0d1) alpha=max(alpha/1.5d0,minalpha) 
752:  
753:             if (debug_run .and. SQNM_DEBUGLEVEL>=1) then  
754:                 write(*,'(a,g12.5,a,f15.7,a,f12.6,a,f15.12,a,f15.7)') 'Step REJECTED. New angle:', angle, & 
755:                     ', rms of gradient: ', RMS, ', Energy: ', energy, ', alpha: ', alpha, ', newE: ', enew 
756:             end if 
757:         else accept_check 
758:             if (debug_run .and. SQNM_DEBUGLEVEL>=1) then  
759:                 write(*,'(a,g12.5,a,f15.7,a,f12.6,a,f15.12,a,f14.10)') 'Step ACCEPTED. New angle:', angle, & 
760:                     ', rms of gradient: ', RMS, ',New Energy: ', enew, ', alpha: ', alpha, ', E decrease: ', enew-energy 
761:             end if 
762:             acceptstep=.true. !step accepted! 
763:             !adjust histories, and current coordinates, gradient, and energy. 
764:             energy=enew 
765:             !first cycle old coordinates by shifting columns. Note that if there is an xcoord that is  
766:             !nhistmax old, it is deleted. Same for gradients 
767:             do i=nhistmax-1,1,-1 
768:                 xhist(:,i+1)=xhist(:,i) 
769:                 gradhist(:,i+1)=gradhist(:,i) 
770:             end do 
771:             !now insert last previous x-coord and grads: 
772:             xhist(:,1)=xcoords(:) 
773:             gradhist(:,1)=grad(:) 
774:              
775:             grad(:)=newgrad(:) 
776:             xcoords(:)=newxcoords(:) 
777:             !increment history length if allowed  
778:             if (nhist<nhistmax) nhist=nhist+1 
779:  
780:             !algorithm seems to occasionally get stuck with small steps when nhist is high; this resets it 
781:             if (nhist==nhistmax) nhist=0 
782:  
783:             !Adjust alpha 
784:             !alpha = alpha * alphascale 
785:             if (angle>anglecutoff) then 
786:                 alpha=alpha*alphascale 
787:             else 
788:                 alpha=max(alpha/alphascale, minalpha) 
789:             end if 
790:         end if accept_check 
791:         !**********DONE CHECKING FOR STEP ACCEPTANCE********** 
792:  
793:         !check for convergence.  
794:         gradrms=norm(numcoords,grad) 
795:         if (RMS < maxrmsgrad) then  
796:             !converged! Yay!  
797:         if(debug_run .and. SQNM_DEBUGLEVEL>=0) then 
798:             write(*, '(a, i4,a,i4)') 'Run ', runs, ' has converged! num iterations: ', iter 
799:             write(*,'(a,f12.8,a,f12.6)') 'Final rms of gradient = ', gradrms, ', Final energy: ', energy 
800:             print *, 'Total potential energy function evaluations: ', totalfunceval 
801:         end if 
802:             converget=.true. 
803:             return 
804:         end if 
805:         !increment iterations. If not converged, loop on.  
806:         iter=iter+1 
807:     end do mainloop 
808:  
809:  
810:     !algorithm has not converged. 
811:    if(debug_run .and. SQNM_DEBUGLEVEL>=0) then 
812:         write(*, '(a, i4,a,i4)') 'Run ', runs, ' failed to converge.' 
813:         write(*,'(a,f12.8,a,f12.6)') 'Final rms of gradient = ', gradrms, ', Final energy: ', energy 
814:         print *, 'Total potential energy function evaluations: ', totalfunceval 
815:    end if 
816: end subroutine sqnm 
817:  
818: !================================================================================================================================ 
819: !================================================================================================================================== 
820: !     ****************************** 
821: !     LINE SEARCH ROUTINE linesearch 
822: !     ****************************** 
823:  
824: SUBROUTINE linesearch(N,X,F,G,S,STP,FTOL,MAXFEV,INFO,NFEV) 
825:     use sqnm_func, only: dot 
826:     implicit none 
827:     integer, intent(in)             :: N,MAXFEV 
828:     integer, intent(out)            :: INFO 
829:     integer, intent(inout)          :: NFEV 
830:     double precision,intent(inout)  :: F,STP 
831:     double precision, intent(in)    :: FTOL 
832:     double precision, intent(inout) :: X(N),G(N) 
833:     double precision, intent(in)    :: S(N) 
834:  
835: !          SUBROUTINE linesearch  
836: ! Taken from Nocedal's L-BFGS algorith (MCSRCH) at http://users.iems.northwestern.edu/~nocedal/lbfgs.html. 
837: ! J. Nocedal. Updating Quasi-Newton Matrices with Limited Storage (1980), Mathematics of Computation 35, pp. 773-782. 
838: ! Note in Nocedal's implemntation, it is called MCSRCH, not linesearch. This does not include any of the actual LBFGS algorithm 
839: !  
840: ! Modifications were made to bring up to Fortran 90 standards.  GOTO/CONTINUE commands were removed, and no returns to the  
841: ! parent function are present since the energies and gradients are obtained from potential.f90, and several do loops were  
842: ! removed in favour of operations on arrays. No interaction with the calling function occurs. Typically, things in lowercase 
843: ! were added by James Klimavicz March 2017 
844: ! 
845: ! THE PURPOSE OF linesearch IS TO FIND A STEP WHICH SATISFIES A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION. 
846: ! AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF UNCERTAINTY  
847: ! IS INITIALLY CHOSEN SO THAT IT CONTAINS A MINIMIZER OF THE MODIFIED FUNCTION 
848: !     F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S). 
849: ! IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, THEN THE  
850: ! INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT CONTAINS A MINIMIZER OF F(X+STP*S). 
851: ! THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES THE SUFFICIENT DECREASE CONDITION 
852: !     F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S), 
853: ! AND THE CURVATURE CONDITION 
854: !     ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S). 
855: ! IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES BOTH 
856: ! CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING ERRORS PREVENT 
857: ! FURTHER PROGRESS. IN THIS CASE STP ONLY SATISFIES THE SUFFICIENT DECREASE CONDITION. 
858: !  
859: ! THE SUBROUTINE STATEMENT IS SUBROUTINE  
860: !       MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA) WHERE 
861: ! N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER OF VARIABLES. 
862: ! X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS X + STP*S. 
863: ! F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S. 
864: ! G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT OF F AT X + STP*S. 
865: ! S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE SEARCH DIRECTION. 
866: ! STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT STP CONTAINS  
867: !     THE FINAL ESTIMATE. 
868: ! FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL  
869: !     DERIVATIVE CONDITION ARE SATISFIED. 
870: ! STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. 
871: ! MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV BY THE END OF  
872: ! AN ITERATION. 
873: ! INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: 
874: !   INFO = 0  IMPROPER INPUT PARAMETERS. 
875: !   INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION HOLD. 
876: !   INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY IS AT MOST XTOL. 
877: !   INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV. 
878: !   INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN. 
879: !   INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX. 
880: !   INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS. THERE MAY NOT BE A STEP WHICH SATISFIES THE 
881: !          SUFFICIENT DECREASE AND CURVATURE CONDITIONS. TOLERANCES MAY BE TOO SMALL. 
882: !   info = 7  search direction was not a descent direction 
883: ! NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF CALLS TO FCN. 
884: ! WA IS A WORK ARRAY OF LENGTH N. 
885: ! SUBPROGRAMS CALLED: MCSTEP 
886: ! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 JORGE J. MORE', DAVID J. THUENTE 
887:  
888:     double precision, parameter :: GTOL = 1.0d-1 !controls the accuracy of the line search routine MCSRCH. 
889:     double precision, parameter :: STPMIN = 1.0d-20 
890:     double precision, parameter :: STPMAX = 1.0d5 
891:     double precision, parameter :: XTOL = epsilon(1.0d0)  
892:  
893:     INTEGER :: INFOC=6 
894:     LOGICAL :: BRACKT,STAGE1 
895:     DOUBLE PRECISION :: DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,FINIT,FTEST1,FM,FX,FXM,FY,FYM,STX,STY 
896:     DOUBLE PRECISION :: STMIN,STMAX,WIDTH,WIDTH1, WA(N) 
897:     double precision, parameter :: P5 = 0.5d0 
898:     double precision, parameter :: P66 = 2.0d0/3.0d0 
899:     double precision, parameter :: XTRAPF=4.0d0 
900:     double precision, parameter :: zero=0.0d0 
901:  
902:     INFO=0 
903:  
904:     ! CHECK THE INPUT PARAMETERS FOR ERRORS. 
905:  
906:     IF (N<= 0.OR.STP<=ZERO.OR.FTOL<ZERO.OR.GTOL<ZERO.OR.XTOL<ZERO.OR.STPMIN<ZERO.OR.STPMAX<STPMIN.OR.MAXFEV<0) RETURN 
907:  
908:     ! COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION AND CHECK THAT S IS A DESCENT DIRECTION. 
909:     DGINIT=dot(N,G,S) 
910:         IF (DGINIT .GE. ZERO) then 
911:             info = 8         
912:             RETURN 
913:         ENDIF 
914:  
915:     !INITIALIZE LOCAL VARIABLES. 
916:     BRACKT = .FALSE. 
917:     STAGE1 = .TRUE. 
918:     NFEV = 0 
919:     FINIT = F 
920:     DGTEST = FTOL*DGINIT 
921:     WIDTH = STPMAX - STPMIN 
922:     WIDTH1 = WIDTH/P5 
923:     WA(:) = X(:) 
924:  
925:     ! THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. 
926:     ! THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF THE INTERVAL  
927:     ! OF UNCERTAINTY. THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. 
928:     STX = ZERO 
929:     FX = FINIT 
930:     DGX = DGINIT 
931:     STY = ZERO 
932:     FY = FINIT 
933:     DGY = DGINIT 
934:  
935:     ! START OF ITERATION. 
936:     do while (info==0) 
937:         ! SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND TO THE PRESENT INTERVAL OF UNCERTAINTY. 
938:         IF (BRACKT) THEN 
939:             STMIN = MIN(STX,STY) 
940:             STMAX = MAX(STX,STY) 
941:         ELSE 
942:             STMIN = STX 
943:             STMAX = STP + XTRAPF*(STP - STX) 
944:         END IF 
945:  
946:         ! FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. 
947:          STP = MAX(STP,STPMIN) 
948:          STP = MIN(STP,STPMAX) 
949:  
950:         ! IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET STP BE THE LOWEST POINT OBTAINED SO FAR. 
951:         IF ((BRACKT .AND. (STP <= STMIN .OR. STP >= STMAX)) .OR. NFEV >= MAXFEV-1 .OR. INFOC == 0 .OR. & 
952:             (BRACKT .AND. STMAX-STMIN <= XTOL*STMAX)) STP = STX 
953:  
954:         ! EVALUATE THE FUNCTION AND GRADIENT AT STP AND COMPUTE THE DIRECTIONAL DERIVATIVE. 
955:         X(:) = WA(:)+STP*S(:) 
956:         call potential(X,G,F,N,.FALSE.) 
957:         NFEV = NFEV + 1 
958:         DG = dot(N,G,S) 
959:  
960:         FTEST1 = FINIT + STP*DGTEST 
961:  
962:         ! TEST FOR CONVERGENCE. 
963:         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX)) .OR. INFOC .EQ. 0) INFO = 6 
964:         IF (STP .EQ. STPMAX .AND. F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5 
965:         IF (STP .EQ. STPMIN .AND. (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4 
966:         IF (NFEV .GE. MAXFEV) INFO = 3 
967:         IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2 
968:         IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1 
969:  
970:         ! CHECK FOR TERMINATION. 
971:         IF (INFO .NE. 0) RETURN 
972:  
973:         ! IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. 
974:         IF (STAGE1 .AND. F .LE. FTEST1 .AND. DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE. 
975:  
976:         ! A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED 
977:         ! FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN 
978:         ! OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. 
979:         IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN 
980:             ! DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. 
981:             FM = F - STP*DGTEST 
982:             FXM = FX - STX*DGTEST 
983:             FYM = FY - STY*DGTEST 
984:             DGM = DG - DGTEST 
985:             DGXM = DGX - DGTEST 
986:             DGYM = DGY - DGTEST 
987:             ! CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY AND TO COMPUTE THE NEW STEP. 
988:             CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,BRACKT,STMIN,STMAX,INFOC) 
989:  
990:             ! RESET THE FUNCTION AND GRADIENT VALUES FOR F. 
991:             FX = FXM + STX*DGTEST 
992:             FY = FYM + STY*DGTEST 
993:             DGX = DGXM + DGTEST 
994:             DGY = DGYM + DGTEST 
995:         ELSE 
996:             ! CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY AND TO COMPUTE THE NEW STEP. 
997:             CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,BRACKT,STMIN,STMAX,INFOC) 
998:         END IF 
999:  
1000:         ! FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE INTERVAL OF UNCERTAINTY. 
1001:         IF (BRACKT) THEN 
1002:             IF (ABS(STY-STX) .GE. P66*WIDTH1) STP = STX + P5*(STY - STX) 
1003:             WIDTH1 = WIDTH 
1004:             WIDTH = ABS(STY-STX) 
1005:         END IF 
1006:          
1007:     end do 
1008: end subroutine linesearch 
1009:  
1010: !=========================================================================================================== 
1011: SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,STPMIN,STPMAX,INFO) 
1012: implicit none 
1013:     INTEGER INFO 
1014:     DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX 
1015:     LOGICAL BRACKT,BOUND 
1016:  
1017: ! THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR 
1018: ! A MINIMIZER OF THE FUNCTION. 
1019: ! THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS 
1020: ! ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A 
1021: ! MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY WITH ENDPOINTS STX AND STY. 
1022: ! THE SUBROUTINE STATEMENT IS SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,STPMIN,STPMAX,INFO) WHERE 
1023: ! STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP, THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED 
1024: !     SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE 
1025: !     SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY. 
1026: ! STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP, THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF 
1027: !     THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY. 
1028: ! STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP, THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP. 
1029: !     IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP. 
1030: !     BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED 
1031: !     THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE. 
1032: ! STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. 
1033: ! INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED 
1034: ! ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS. 
1035: ! ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 JORGE J. MORE', DAVID J. THUENTE 
1036:  
1037:     DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA 
1038:     INFO = 0 
1039:  
1040:     ! CHECK THE INPUT PARAMETERS FOR ERRORS. 
1041:     IF ((BRACKT .AND. (STP <= MIN(STX,STY) .OR. STP >= MAX(STX,STY))) .OR. DX*(STP-STX) >= 0.0 .OR. STPMAX < STPMIN) RETURN 
1042:  
1043:     ! DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. 
1044:     SGND = DP*(DX/ABS(DX)) 
1045:  
1046:     ! FIRST CASE. A HIGHER FUNCTION VALUE. THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER TO STX THAN THE QUADRATIC STEP,  
1047:     ! THE CUBIC STEP IS TAKEN, ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. 
1048:     cases: IF (FP .GT. FX) THEN 
1049:         INFO = 1 
1050:         BOUND = .TRUE. 
1051:         THETA = 3*(FX - FP)/(STP - STX) + DX + DP 
1052:         S = MAX(ABS(THETA),ABS(DX),ABS(DP)) 
1053:         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) 
1054:         IF (STP .LT. STX) GAMMA = -GAMMA 
1055:         P = (GAMMA - DX) + THETA 
1056:         Q = ((GAMMA - DX) + GAMMA) + DP 
1057:         R = P/Q 
1058:         STPC = STX + R*(STP - STX) 
1059:         STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX) 
1060:         IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN 
1061:             STPF = STPC 
1062:         ELSE 
1063:             STPF = STPC + (STPQ - STPC)/2 
1064:         END IF 
1065:         BRACKT = .TRUE. 
1066:  
1067:     ! SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC 
1068:     ! STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. 
1069:     ELSE IF (SGND .LT. 0.0) THEN cases 
1070:         INFO = 2 
1071:         BOUND = .FALSE. 
1072:         THETA = 3*(FX - FP)/(STP - STX) + DX + DP 
1073:         S = MAX(ABS(THETA),ABS(DX),ABS(DP)) 
1074:         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) 
1075:         IF (STP .GT. STX) GAMMA = -GAMMA 
1076:         P = (GAMMA - DP) + THETA 
1077:         Q = ((GAMMA - DP) + GAMMA) + DX 
1078:         R = P/Q 
1079:         STPC = STP + R*(STX - STP) 
1080:         STPQ = STP + (DP/(DP-DX))*(STX - STP) 
1081:         IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN 
1082:             STPF = STPC 
1083:         ELSE 
1084:             STPF = STPQ 
1085:         END IF 
1086:         BRACKT = .TRUE. 
1087:  
1088:     ! THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. THE CUBIC 
1089:     ! STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC IS BEYOND STP.  
1090:     ! OTHERWISE THE CUBIC STEP IS DEFINED TO BE EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO COMPUTED AND IF THE  
1091:     ! MINIMUM IS BRACKETED THEN THE THE STEP CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. 
1092:     ELSE IF (ABS(DP) .LT. ABS(DX)) THEN cases 
1093:         INFO = 3 
1094:         BOUND = .TRUE. 
1095:         THETA = 3*(FX - FP)/(STP - STX) + DX + DP 
1096:         S = MAX(ABS(THETA),ABS(DX),ABS(DP)) 
1097:         !THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND TO INFINITY IN THE DIRECTION OF THE STEP. 
1098:         GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S))) 
1099:         IF (STP .GT. STX) GAMMA = -GAMMA 
1100:         P = (GAMMA - DP) + THETA 
1101:         Q = (GAMMA + (DX - DP)) + GAMMA 
1102:         R = P/Q 
1103:         IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN 
1104:             STPC = STP + R*(STX - STP) 
1105:         ELSE IF (STP .GT. STX) THEN 
1106:             STPC = STPMAX 
1107:         ELSE 
1108:             STPC = STPMIN 
1109:         END IF 
1110:         STPQ = STP + (DP/(DP-DX))*(STX - STP) 
1111:         IF (BRACKT) THEN 
1112:             IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN 
1113:                 STPF = STPC 
1114:             ELSE 
1115:                 STPF = STPQ 
1116:             END IF 
1117:         ELSE 
1118:             IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN 
1119:                 STPF = STPC 
1120:             ELSE 
1121:                 STPF = STPQ 
1122:             END IF 
1123:         END IF 
1124:  
1125:     ! FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES 
1126:     ! NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. 
1127:     ELSE cases 
1128:         INFO = 4 
1129:         BOUND = .FALSE. 
1130:         IF (BRACKT) THEN 
1131:             THETA = 3*(FP - FY)/(STY - STP) + DY + DP 
1132:             S = MAX(ABS(THETA),ABS(DY),ABS(DP)) 
1133:             GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S)) 
1134:             IF (STP .GT. STY) GAMMA = -GAMMA 
1135:             P = (GAMMA - DP) + THETA 
1136:             Q = ((GAMMA - DP) + GAMMA) + DY 
1137:             R = P/Q 
1138:             STPC = STP + R*(STY - STP) 
1139:             STPF = STPC 
1140:         ELSE IF (STP .GT. STX) THEN 
1141:             STPF = STPMAX 
1142:         ELSE 
1143:             STPF = STPMIN 
1144:         END IF 
1145:     END IF cases 
1146:  
1147:     ! UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. 
1148:     IF (FP .GT. FX) THEN 
1149:         STY = STP 
1150:         FY = FP 
1151:         DY = DP 
1152:     ELSE 
1153:         IF (SGND .LT. 0.0) THEN 
1154:             STY = STX 
1155:             FY = FX 
1156:             DY = DX 
1157:             END IF 
1158:         STX = STP 
1159:         FX = FP 
1160:         DX = DP 
1161:     END IF 
1162:  
1163:     ! COMPUTE THE NEW STEP AND SAFEGUARD IT. 
1164:     STPF = MIN(STPMAX,STPF) 
1165:     STPF = MAX(STPMIN,STPF) 
1166:     STP = STPF 
1167:     IF (BRACKT .AND. BOUND) THEN 
1168:         IF (STY .GT. STX) THEN 
1169:             STP = MIN(STX+0.66*(STY-STX),STP) 
1170:         ELSE 
1171:             STP = MAX(STX+0.66*(STY-STX),STP) 
1172:         END IF 
1173:     END IF 
1174:     RETURN 
1175: END subroutine mcstep  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0