hdiff output

r33372/bhinterp.f90 2017-10-04 18:30:16.540295608 +0100 r33371/bhinterp.f90 2017-10-04 18:30:20.816352013 +0100
157:       IF ((SFRAC.LE.0).OR.(SFRAC.GE.1.D0)) THEN157:       IF ((SFRAC.LE.0).OR.(SFRAC.GE.1.D0)) THEN
158:          PRINT *,'bhinterp> could not find acceptable interpolated minimum'158:          PRINT *,'bhinterp> could not find acceptable interpolated minimum'
159:          RETURN159:          RETURN
160:       ENDIF160:       ENDIF
161:       IF(MOD(NTRY,2).EQ.0) NCOUNT=NCOUNT+1161:       IF(MOD(NTRY,2).EQ.0) NCOUNT=NCOUNT+1
162:       PRINT '(A,F8.3)','bhinterp> retry interpolation, SFRAC= ',SFRAC162:       PRINT '(A,F8.3)','bhinterp> retry interpolation, SFRAC= ',SFRAC
163:       GOTO 10163:       GOTO 10
164:    ENDIF164:    ENDIF
165: ENDIF165: ENDIF
166: 166: 
167: WRITE(*,*) "sn402: Changed BHINTERP so that it calls ALIGN_DECIDE instead of MINPERMDIST" 
168: WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
169: WRITE(*,*) "Please check carefully with FASTOVERLAP set that this part of the code is working as you expect, then remove these messages!" 
170: !msb50 changed order CSTART, COORDS because this is only for distance so shouldn't matter167: !msb50 changed order CSTART, COORDS because this is only for distance so shouldn't matter
171: CALL ALIGN_DECIDE(CSTART,COORDS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)168: CALL MINPERMDIST(CSTART,COORDS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
172: DSTART=DIST169: DSTART=DIST
173: CALL ALIGN_DECIDE(CFINISH,COORDS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)170: CALL MINPERMDIST(CFINISH,COORDS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
174: DFINISH=DIST171: DFINISH=DIST
175: 172: 
176: ESPREV=0.5D0*BHK*(DSTART**2+DFINISH**2)173: ESPREV=0.5D0*BHK*(DSTART**2+DFINISH**2)
177: ! PRINT '(A,5G18.8,I8)','EREAL,ESPREV,DS,DF,SFRAC,ITDONE=',EREAL,ESPREV,DSTART,DFINISH,SFRAC,ITDONE174: ! PRINT '(A,5G18.8,I8)','EREAL,ESPREV,DS,DF,SFRAC,ITDONE=',EREAL,ESPREV,DSTART,DFINISH,SFRAC,ITDONE
178: 175: 
179: IF (.NOT.MFLAG) THEN176: IF (.NOT.MFLAG) THEN
180:    PRINT '(2(A,G20.10),A,I8)','bhinterp> WARNING - initial quench failed to converge, energy=',EREAL, &177:    PRINT '(2(A,G20.10),A,I8)','bhinterp> WARNING - initial quench failed to converge, energy=',EREAL, &
181:   &  ' spring energy=',ESPREV,' lbfgs steps=',ITDONE178:   &  ' spring energy=',ESPREV,' lbfgs steps=',ITDONE
182:    PRINT '(A,2G15.5)',' bhinterp> Initial distances: ',DSTART,DFINISH179:    PRINT '(A,2G15.5)',' bhinterp> Initial distances: ',DSTART,DFINISH
183: ENDIF180: ENDIF
188: !185: !
189: IF (MIN(DSTART,DFINISH).LT.1.1D0) THEN186: IF (MIN(DSTART,DFINISH).LT.1.1D0) THEN
190:    IF (BHDEBUG) PRINT '(A,G20.10)',' bhinterp> one of the initial distances is suspiciously small: reoptimise'187:    IF (BHDEBUG) PRINT '(A,G20.10)',' bhinterp> one of the initial distances is suspiciously small: reoptimise'
191:    DSTARTSAVE=DSTART; DFINISHSAVE=DFINISH188:    DSTARTSAVE=DSTART; DFINISHSAVE=DFINISH
192:    GMAXSAVE=GMAX; GMAX=BHCONV189:    GMAXSAVE=GMAX; GMAX=BHCONV
193:    NCOUNT=0190:    NCOUNT=0
194:    DO WHILE (GMAX.GT.GMAXSAVE) 191:    DO WHILE (GMAX.GT.GMAXSAVE) 
195:       GMAX=GMAX/10.0D0192:       GMAX=GMAX/10.0D0
196:       CALL MYLBFGS(NCOORDS,MUPDATE,COORDS,.FALSE.,MFLAG,ENERGY,RMS2,EREAL,RMS,BFGSSTEPS, &193:       CALL MYLBFGS(NCOORDS,MUPDATE,COORDS,.FALSE.,MFLAG,ENERGY,RMS2,EREAL,RMS,BFGSSTEPS, &
197:      &                      .TRUE.,ITDONE,PTEST,VNEW,.TRUE.,.FALSE.)194:      &                      .TRUE.,ITDONE,PTEST,VNEW,.TRUE.,.FALSE.)
198:       WRITE(*,*) "sn402: Changed BHINTERP so that it calls ALIGN_DECIDE instead of MINPERMDIST"195:       CALL MINPERMDIST(COORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
199:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
200:       WRITE(*,*) "Please check carefully with FASTOVERLAP set that this part of the code is working as you expect, then remove these messages!" 
201:       CALL ALIGN_DECIDE(COORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT) 
202:       DSTART=DIST196:       DSTART=DIST
203:       CALL ALIGN_DECIDE(COORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)197:       CALL MINPERMDIST(COORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
204:       DFINISH=DIST198:       DFINISH=DIST
205:       IF (BHDEBUG) PRINT '(A,G20.10,A,2G15.5)',' bhinterp> for RMS condition ',GMAX,' distances: ',DSTART,DFINISH199:       IF (BHDEBUG) PRINT '(A,G20.10,A,2G15.5)',' bhinterp> for RMS condition ',GMAX,' distances: ',DSTART,DFINISH
206:       IF ((ABS(DSTARTSAVE-DSTART)*100/DSTARTSAVE.LT.10.0D0).AND. &200:       IF ((ABS(DSTARTSAVE-DSTART)*100/DSTARTSAVE.LT.10.0D0).AND. &
207:           (ABS(DFINISHSAVE-DFINISH)*100/DFINISHSAVE.LT.10.0D0)) EXIT ! guess that we really have converged201:           (ABS(DFINISHSAVE-DFINISH)*100/DFINISHSAVE.LT.10.0D0)) EXIT ! guess that we really have converged
208:       IF ((DSTART.LT.GEOMDIFFTOL).OR.(DFINISH.LT.GEOMDIFFTOL)) EXIT ! we have hit an end point202:       IF ((DSTART.LT.GEOMDIFFTOL).OR.(DFINISH.LT.GEOMDIFFTOL)) EXIT ! we have hit an end point
209:     203:     
210:       DSTARTSAVE=DSTART; DFINISHSAVE=DFINISH204:       DSTARTSAVE=DSTART; DFINISHSAVE=DFINISH
211:    END DO205:    END DO
212:    GMAX=GMAXSAVE206:    GMAX=GMAXSAVE
213: ENDIF207: ENDIF
295:    ELSE289:    ELSE
296:       DO J2=1,NCOORDS290:       DO J2=1,NCOORDS
297:          COORDS(J2)=COORDS(J2)+BHSTEPSIZE*(DPRAND()-0.5D0)*2.0D0291:          COORDS(J2)=COORDS(J2)+BHSTEPSIZE*(DPRAND()-0.5D0)*2.0D0
298:       ENDDO292:       ENDDO
299:    ENDIF293:    ENDIF
300:    KNOWE=.FALSE.294:    KNOWE=.FALSE.
301:    GMAXSAVE=GMAX; GMAX=BHCONV ! mylbfgs now uses GMAX so that we can change this parameter via changep295:    GMAXSAVE=GMAX; GMAX=BHCONV ! mylbfgs now uses GMAX so that we can change this parameter via changep
302:    CALL MYLBFGS(NCOORDS,MUPDATE,COORDS,.FALSE.,MFLAG,ENERGY,RMS2,EREAL,RMS,BFGSSTEPS, &296:    CALL MYLBFGS(NCOORDS,MUPDATE,COORDS,.FALSE.,MFLAG,ENERGY,RMS2,EREAL,RMS,BFGSSTEPS, &
303:      &                     .TRUE.,ITDONE,PTEST,VNEW,.TRUE.,.FALSE.)297:      &                     .TRUE.,ITDONE,PTEST,VNEW,.TRUE.,.FALSE.)
304:    GMAX=GMAXSAVE298:    GMAX=GMAXSAVE
305:       WRITE(*,*) "sn402: Changed BHINTERP so that it calls ALIGN_DECIDE instead of MINPERMDIST"299:    CALL MINPERMDIST(COORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
306:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
307:       WRITE(*,*) "Please check carefully with FASTOVERLAP set that this part of the code is working as you expect, then remove these messages!" 
308:    CALL ALIGN_DECIDE(COORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT) 
309:    DSTART=DIST300:    DSTART=DIST
310:    CALL ALIGN_DECIDE(COORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)301:    CALL MINPERMDIST(COORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
311:    DFINISH=DIST302:    DFINISH=DIST
312:    ESTRING=0.5D0*BHK*(DSTART**2+DFINISH**2)303:    ESTRING=0.5D0*BHK*(DSTART**2+DFINISH**2)
313:    ETOTAL=EREAL+ESTRING304:    ETOTAL=EREAL+ESTRING
314: !  PRINT '(A,4G20.10)',' bhinterp> E, E+ES, distances=',EREAL,ETOTAL,DSTART,DFINISH305: !  PRINT '(A,4G20.10)',' bhinterp> E, E+ES, distances=',EREAL,ETOTAL,DSTART,DFINISH
315: ! WRITE(99,*) 864306: ! WRITE(99,*) 864
316: ! WRITE(99,'(A,4G20.10)') 'initial interp, E, E+ES, dists=',EREAL,ESPREV,DSTART,DFINISH307: ! WRITE(99,'(A,4G20.10)') 'initial interp, E, E+ES, dists=',EREAL,ESPREV,DSTART,DFINISH
317: ! DO K=1,NATOMS308: ! DO K=1,NATOMS
318: !    WRITE(99,'(A,3G20.10)') 'LA ',COORDS(3*(K-1)+1),COORDS(3*(K-1)+2),COORDS(3*(K-1)+3)309: !    WRITE(99,'(A,3G20.10)') 'LA ',COORDS(3*(K-1)+1),COORDS(3*(K-1)+2),COORDS(3*(K-1)+3)
319: ! ENDDO310: ! ENDDO
320: !311: !
384: 375: 
385: !  IF ((DSTART.LT.DENDPT).AND.(DFINISH.LT.DENDPT)) EXIT ! stop if both distances are < end point distance376: !  IF ((DSTART.LT.DENDPT).AND.(DFINISH.LT.DENDPT)) EXIT ! stop if both distances are < end point distance
386:    377:    
387: ENDDO ! main loop over BH steps378: ENDDO ! main loop over BH steps
388: !379: !
389: !  Reoptimise best minimum to global RMS tolerance. Here we actually use GMAX not BHCONV!380: !  Reoptimise best minimum to global RMS tolerance. Here we actually use GMAX not BHCONV!
390: !381: !
391: KNOWE=.FALSE.382: KNOWE=.FALSE.
392: CALL MYLBFGS(NCOORDS,MUPDATE,BESTCOORDS,.FALSE.,MFLAG,ENERGY,RMS2,EREAL,RMS,BFGSSTEPS, &383: CALL MYLBFGS(NCOORDS,MUPDATE,BESTCOORDS,.FALSE.,MFLAG,ENERGY,RMS2,EREAL,RMS,BFGSSTEPS, &
393:   &          .TRUE.,ITDONE,PTEST,VNEW,.TRUE.,.FALSE.)384:   &          .TRUE.,ITDONE,PTEST,VNEW,.TRUE.,.FALSE.)
394:       WRITE(*,*) "sn402: Changed BHINTERP so that it calls ALIGN_DECIDE instead of MINPERMDIST"385: CALL MINPERMDIST(BESTCOORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
395:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
396:       WRITE(*,*) "Please check carefully with FASTOVERLAP set that this part of the code is working as you expect, then remove these messages!" 
397: CALL ALIGN_DECIDE(BESTCOORDS,CSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT) 
398: DSTART=DIST386: DSTART=DIST
399: CALL ALIGN_DECIDE(BESTCOORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)387: CALL MINPERMDIST(BESTCOORDS,CFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
400: DFINISH=DIST388: DFINISH=DIST
401: !389: !
402: ! Don't allow the interpolated minimum to be either of the endpoints!390: ! Don't allow the interpolated minimum to be either of the endpoints!
403: !391: !
404: IF (BHDEBUG) PRINT '(A,2G20.10)',' bhinterp> DSTART,DFINISH=',DSTART,DFINISH392: IF (BHDEBUG) PRINT '(A,2G20.10)',' bhinterp> DSTART,DFINISH=',DSTART,DFINISH
405: IF ((DSTART.LT.GEOMDIFFTOL).OR.(DFINISH.LT.GEOMDIFFTOL)) THEN393: IF ((DSTART.LT.GEOMDIFFTOL).OR.(DFINISH.LT.GEOMDIFFTOL)) THEN
406:    SUCCESS=.FALSE.394:    SUCCESS=.FALSE.
407:    IF (BHDEBUG) PRINT '(A)',' bhinterp> Tight quench optimised to an end point'395:    IF (BHDEBUG) PRINT '(A)',' bhinterp> Tight quench optimised to an end point'
408: ELSE396: ELSE
409: ! sf344> AMBER stuff397: ! sf344> AMBER stuff


r33372/bisect.f90 2017-10-04 18:30:16.756298458 +0100 r33371/bisect.f90 2017-10-04 18:30:21.060355232 +0100
164:    ENDIF164:    ENDIF
165:    IF(MOD(NBISECT,2).EQ.0) NCOUNT=NCOUNT+1165:    IF(MOD(NBISECT,2).EQ.0) NCOUNT=NCOUNT+1
166:    PRINT '(A,F8.3)',' bisect> retry interpolation, SFRAC= ',SFRAC166:    PRINT '(A,F8.3)',' bisect> retry interpolation, SFRAC= ',SFRAC
167:    GOTO 12167:    GOTO 12
168: ENDIF168: ENDIF
169: IF ((.NOT.MFLAG).OR.(EREAL.GT.BISECTMAXENERGY)) THEN169: IF ((.NOT.MFLAG).OR.(EREAL.GT.BISECTMAXENERGY)) THEN
170:    TRIED(NDO)=.TRUE.170:    TRIED(NDO)=.TRUE.
171:    GOTO 11171:    GOTO 11
172: ENDIF172: ENDIF
173: 173: 
174: WRITE(*,*) "sn402: Changed BISECT so that it calls ALIGN_DECIDE instead of MINPERMDIST"174: CALL MINPERMDIST(COORDS,TEMPSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
175: WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
176: WRITE(*,*) "Please set FASTOVERLAP and check carefully that this part of the code is working as you expect, then remove these messages!" 
177: CALL ALIGN_DECIDE(COORDS,TEMPSTART,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT) 
178: DSTART=DIST175: DSTART=DIST
179: CALL ALIGN_DECIDE(COORDS,TEMPFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)176: CALL MINPERMDIST(COORDS,TEMPFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
180: DFINISH=DIST177: DFINISH=DIST
181: ! PRINT '(A,4G18.8,I8)','EREAL,DS,DF,SFRAC,ITDONE=',EREAL,DSTART,DFINISH,SFRAC,ITDONE178: ! PRINT '(A,4G18.8,I8)','EREAL,DS,DF,SFRAC,ITDONE=',EREAL,DSTART,DFINISH,SFRAC,ITDONE
182: IF (.NOT.MFLAG) THEN179: IF (.NOT.MFLAG) THEN
183:    PRINT '(A,G20.10,A,I8)',' bisect> WARNING - initial quench failed to converge, energy=',EREAL, &180:    PRINT '(A,G20.10,A,I8)',' bisect> WARNING - initial quench failed to converge, energy=',EREAL, &
184:   &  ' lbfgs steps=',ITDONE181:   &  ' lbfgs steps=',ITDONE
185:    PRINT '(A,2G15.5)',' bisect> distances from bracketing minima: ',DSTART,DFINISH182:    PRINT '(A,2G15.5)',' bisect> distances from bracketing minima: ',DSTART,DFINISH
186: ENDIF183: ENDIF
187: IF (BISECTDEBUG) PRINT '(A,G20.10,A,I8)',' bisect> quench energy=',EREAL,' lbfgs steps=',ITDONE184: IF (BISECTDEBUG) PRINT '(A,G20.10,A,I8)',' bisect> quench energy=',EREAL,' lbfgs steps=',ITDONE
188: IF (BISECTDEBUG) PRINT '(A,2G15.5)',' bisect> distances from bracketing minima: ',DSTART,DFINISH185: IF (BISECTDEBUG) PRINT '(A,2G15.5)',' bisect> distances from bracketing minima: ',DSTART,DFINISH
189: 186: 
224:    IF (BISECTDEBUG) PRINT '(A,G20.10)',' bisect> Initial guess minimised to second end point - change fraction to ',SFRAC221:    IF (BISECTDEBUG) PRINT '(A,G20.10)',' bisect> Initial guess minimised to second end point - change fraction to ',SFRAC
225:    HITF=.TRUE.222:    HITF=.TRUE.
226:    GOTO 12223:    GOTO 12
227: ENDIF224: ENDIF
228: !225: !
229: ! Now check that the minimum isn;t identical to any of the others. Don't add it to the list if so.226: ! Now check that the minimum isn;t identical to any of the others. Don't add it to the list if so.
230: !227: !
231: DO J1=1,NMIN228: DO J1=1,NMIN
232:    IF (ABS(EREAL-EMIN(J1)).GT.EDIFFTOL) CYCLE229:    IF (ABS(EREAL-EMIN(J1)).GT.EDIFFTOL) CYCLE
233:    CTEMP(1:3*NATOMS)=MINCOORDS(1:3*NATOMS,J1)230:    CTEMP(1:3*NATOMS)=MINCOORDS(1:3*NATOMS,J1)
234:    CALL ALIGN_DECIDE(COORDS,CTEMP,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)231:    CALL MINPERMDIST(COORDS,CTEMP,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
235:    IF (DIST.LT.GEOMDIFFTOL) THEN232:    IF (DIST.LT.GEOMDIFFTOL) THEN
236:       PRINT '(A,I8,A)',' bisect> quench minimum is the same as number ',J1,' not adding it to the list'233:       PRINT '(A,I8,A)',' bisect> quench minimum is the same as number ',J1,' not adding it to the list'
237:       TRIED(NDO)=.TRUE.234:       TRIED(NDO)=.TRUE.
238:       GOTO 11235:       GOTO 11
239:    ENDIF236:    ENDIF
240: ENDDO237: ENDDO
241: 238: 
242: !239: !
243: ! If we have reached this point then we need to add the new minimum to the list and 240: ! If we have reached this point then we need to add the new minimum to the list and 
244: ! then go back to the beginning. Must set the TRIED values appropriately. Connections to 241: ! then go back to the beginning. Must set the TRIED values appropriately. Connections to 


r33372/connect.f 2017-10-04 18:30:16.980301412 +0100 r33371/connect.f 2017-10-04 18:30:21.628363046 +0100
1795:       ENDDO1795:       ENDDO
1796: C 1796: C 
1797: C The declared dimension of VECS is 3*NATOMS, so the above DO loop is OK but we need to get the dimension right below.1797: C The declared dimension of VECS is 3*NATOMS, so the above DO loop is OK but we need to get the dimension right below.
1798: C1798: C
1799:       IF (UNRST) THEN1799:       IF (UNRST) THEN
1800:          CALL VECNORM(VECS,NINTS)1800:          CALL VECNORM(VECS,NINTS)
1801:       ELSE1801:       ELSE
1802:          CALL VECNORM(VECS,NOPT)1802:          CALL VECNORM(VECS,NOPT)
1803:       ENDIF1803:       ENDIF
1804:       FIXDSAVE=FIXD1804:       FIXDSAVE=FIXD
1805:       CALL ALIGN_DECIDE(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTPF,DIST2,RIGIDBODY,RMAT)1805:       CALL MINPERMDIST(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTPF,DIST2,RIGIDBODY,RMAT)
1806:       EINITIAL=EMIN(JSTART)1806:       EINITIAL=EMIN(JSTART)
1807:       EFINAL=EMIN(JFINISH)1807:       EFINAL=EMIN(JFINISH)
1808:       DO J1=1,NOPT1808:       DO J1=1,NOPT
1809:          QSTART(J1)=Q(J1)1809:          QSTART(J1)=Q(J1)
1810:       ENDDO1810:       ENDDO
1811:       WRITE(*,'(A,I5,A,I5,A,F15.5)')' Minima ',JSTART,' and ',JFINISH,' unconnected: minimum distance=',DISTPF1811:       WRITE(*,'(A,I5,A,I5,A,F15.5)')' Minima ',JSTART,' and ',JFINISH,' unconnected: minimum distance=',DISTPF
1812:       IF (DISTPF.LT.GEOMDIFFTOL) THEN1812:       IF (DISTPF.LT.GEOMDIFFTOL) THEN
1813:          WRITE(*,'(A,I6,A,I6,A)') ' Minima ',JSTART,' and ',JFINISH,' appear to be identical - quit'1813:          WRITE(*,'(A,I6,A,I6,A)') ' Minima ',JSTART,' and ',JFINISH,' appear to be identical - quit'
1814:          STOP1814:          STOP
1815:       ENDIF1815:       ENDIF


r33372/gay-berne.f90 2017-10-04 18:30:17.204304366 +0100 r33371/gay-berne.f90 2017-10-04 18:30:21.864365837 +0100
2087:         WRITE(*,'(A)') 'ERROR - number of atoms determined from OPTIM does not match '// &2087:         WRITE(*,'(A)') 'ERROR - number of atoms determined from OPTIM does not match '// &
2088: &                        ' with the number of atoms determined from the path.xyz file!!!'2088: &                        ' with the number of atoms determined from the path.xyz file!!!'
2089:         STOP2089:         STOP
2090:    END IF2090:    END IF
2091:    DO J2=1,REALNATOMS2091:    DO J2=1,REALNATOMS
2092:         READ(133,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &2092:         READ(133,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &
2093:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))2093:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))
2094:    END DO2094:    END DO
2095:    COORDSDUMMY(:)=COORDSNEW(:)2095:    COORDSDUMMY(:)=COORDSNEW(:)
2096: !   COORDSDUMMY(50)=COORDSNEW(50)*10.0D02096: !   COORDSDUMMY(50)=COORDSNEW(50)*10.0D0
2097: WRITE(*,*) "sn402: Changed PYREALIGNXYZ so that it calls ALIGN_DECIDE instead of MINPERMDIST"2097:    CALL MINPERMDIST(COORDSREF,COORDSDUMMY,NATOMS,.false.,1.0,1.0,1.0,.false.,.false.,D,DIST2,.false.,RMAT)
2098: WRITE(*,*) "I haven't tested this change and suspect that it may not work (especially if the coordinates are non-cartesian)." 2098:    D=D**2 ! since MINPERMDIST now returns the distance
2099: WRITE(*,*) "Please set FASTOVERLAP and check carefully that this part of the code is working as you expect, then remove these messages!" 
2100:    CALL ALIGN_DECIDE(COORDSREF,COORDSDUMMY,NATOMS,.false.,1.0,1.0,1.0,.false.,.false.,D,DIST2,.false.,RMAT) 
2101:    D=D**2 ! since ALIGN_DECIDE now returns the distance 
2102:    COORDSNEW(:)=COORDSDUMMY(:)2099:    COORDSNEW(:)=COORDSDUMMY(:)
2103:    WRITE(*,*) 'writing frame ', J1+12100:    WRITE(*,*) 'writing frame ', J1+1
2104:    WRITE(134,*) REALNATOMS2101:    WRITE(134,*) REALNATOMS
2105:    WRITE(134,*) dummyline2102:    WRITE(134,*) dummyline
2106:    DO J2=1,REALNATOMS2103:    DO J2=1,REALNATOMS
2107:      WRITE(134,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &2104:      WRITE(134,'(A1,6X,6G20.10)') ATOMLABELS(J2),COORDSNEW(3*J2-2),COORDSNEW(3*J2-1),COORDSNEW(3*J2), &
2108:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))2105:         & COORDSNEW(3*(J2+REALNATOMS)-2),COORDSNEW(3*(J2+REALNATOMS)-1),COORDSNEW(3*(J2+REALNATOMS))
2109:    END DO2106:    END DO
2110:    COORDSREF(:)=COORDSNEW(:)2107:    COORDSREF(:)=COORDSNEW(:)
2111:  END DO2108:  END DO
2186:         !WRITE(MYUNIT,*) COORDSSTORE(:,3)2183:         !WRITE(MYUNIT,*) COORDSSTORE(:,3)
2187:         !WRITE(MYUNIT,*) COORDSSTORE(:,2)2184:         !WRITE(MYUNIT,*) COORDSSTORE(:,2)
2188:         !WRITE(MYUNIT,*) COORDSSTORE(:,4)2185:         !WRITE(MYUNIT,*) COORDSSTORE(:,4)
2189:         DO J1=1,32186:         DO J1=1,3
2190:          COORDSDUMMY(3*(RANDOM1-1)+J1)=COORDSSTORE(J1,2)2187:          COORDSDUMMY(3*(RANDOM1-1)+J1)=COORDSSTORE(J1,2)
2191:          COORDSDUMMY(3*(RANDOM2-1)+J1)=COORDSSTORE(J1,1)2188:          COORDSDUMMY(3*(RANDOM2-1)+J1)=COORDSSTORE(J1,1)
2192:          COORDSDUMMY(3*NATOMS/2+3*(RANDOM1-1)+J1)=COORDSSTORE(J1,4)2189:          COORDSDUMMY(3*NATOMS/2+3*(RANDOM1-1)+J1)=COORDSSTORE(J1,4)
2193:          COORDSDUMMY(3*NATOMS/2+3*(RANDOM2-1)+J1)=COORDSSTORE(J1,3)2190:          COORDSDUMMY(3*NATOMS/2+3*(RANDOM2-1)+J1)=COORDSSTORE(J1,3)
2194:         END DO2191:         END DO
2195: 2192: 
2196:    IF(MOD(J2,1000)==0)   then2193:    IF(MOD(J2,1000)==0)   then 
2197: WRITE(*,*) "sn402: Changed PYRANDOMSWAP so that it calls ALIGN_DECIDE instead of MINPERMDIST"2194:    CALL MINPERMDIST(COORDSREF,COORDSDUMMY,NATOMS,.false.,1.0,1.0,1.0,.false.,.false.,D,D,.false.,RMAT)
2198: WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
2199: WRITE(*,*) "Please set FASTOVERLAP and check carefully that this part of the code is working as you expect, then remove these messages!"  
2200:    CALL ALIGN_DECIDE(COORDSREF,COORDSDUMMY,NATOMS,.false.,1.0,1.0,1.0,.false.,.false.,D,D,.false.,RMAT) 
2201:    D=D**22195:    D=D**2
2202:         IF(D<DBEST) THEN2196:         IF(D<DBEST) THEN
2203:           DBEST=D2197:           DBEST=D
2204:           BESTCOORDS(:)=COORDSDUMMY(:)2198:           BESTCOORDS(:)=COORDSDUMMY(:)
2205: !        ELSE IF (D-DBEST<0.1D0.AND.DPRAND()<0.5D0) THEN2199: !        ELSE IF (D-DBEST<0.1D0.AND.DPRAND()<0.5D0) THEN
2206: !          DBEST=D2200: !          DBEST=D
2207: !          BESTCOORDS(:)=COORDSDUMMY(:)2201: !          BESTCOORDS(:)=COORDSDUMMY(:)
2208:         END IF2202:         END IF
2209:    END IF2203:    END IF
2210: END DO2204: END DO


r33372/geopt.f 2017-10-04 18:30:17.428307320 +0100 r33371/geopt.f 2017-10-04 18:30:22.084368740 +0100
627:                      Q(3*(J2-1)+3)=QSAVE(3*(J2-1)+3)+(TSDISP*VECS(3*(J2-1)+3))/SQRT(ATMASS(J2))627:                      Q(3*(J2-1)+3)=QSAVE(3*(J2-1)+3)+(TSDISP*VECS(3*(J2-1)+3))/SQRT(ATMASS(J2))
628:                   ENDDO628:                   ENDDO
629:                   DO I1=1,NORDER629:                   DO I1=1,NORDER
630:                      IF ((CHRMMT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN630:                      IF ((CHRMMT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN
631:                         CALL GETDIHE(Q,ORDER(I1),ORDERNUM(I1))631:                         CALL GETDIHE(Q,ORDER(I1),ORDERNUM(I1))
632:                      ELSEIF ((AMBERT.OR.NABT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN632:                      ELSEIF ((AMBERT.OR.NABT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN
633:                         IF (I1.EQ.1) CALL AMBERDIHEDR(Q,NATOMS,5,7,9,15,ORDER(I1))633:                         IF (I1.EQ.1) CALL AMBERDIHEDR(Q,NATOMS,5,7,9,15,ORDER(I1))
634:                         IF (I1.EQ.2) CALL AMBERDIHEDR(Q,NATOMS,7,9,15,17,ORDER(I1))634:                         IF (I1.EQ.2) CALL AMBERDIHEDR(Q,NATOMS,7,9,15,17,ORDER(I1))
635:                      ELSE 635:                      ELSE 
636: !636: !
637: ! align_decide selects a subroutine that will align the second coordinate array with the first. We don't want to change Q637: ! minpermdist aligns second coordinate array with first. We don't want to change Q
638: ! so choose the order carefully.638: ! so choose the order carefully.
639: !639: !
640:                         IF (I1.EQ.1) THEN640:                         IF (I1.EQ.1) THEN
641:                            CALL ALIGN_DECIDE(Q,POINTSDECA,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDER(1),DIST2,641:                            CALL MINPERMDIST(Q,POINTSDECA,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDER(1),DIST2,
642:      &                                                RIGIDBODY,RMAT)642:      &                                                RIGIDBODY,RMAT)
643:                         ELSE IF (I1.EQ.2) THEN643:                         ELSE IF (I1.EQ.2) THEN
644:                            CALL ALIGN_DECIDE(Q,POINTSICOS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDER(2),DIST2,644:                            CALL MINPERMDIST(Q,POINTSICOS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDER(2),DIST2,
645:      &                                                RIGIDBODY,RMAT)645:      &                                                RIGIDBODY,RMAT)
646:                         ELSE IF (I1.GT.2) THEN646:                         ELSE IF (I1.GT.2) THEN
647:                            DO K1=1,3647:                            DO K1=1,3
648:                               DO K2=1,NATOMS648:                               DO K2=1,NATOMS
649:                                  QTEMP(K1,K2)=Q(3*(K2-1)+K1)649:                                  QTEMP(K1,K2)=Q(3*(K2-1)+K1)
650:                               ENDDO650:                               ENDDO
651:                            ENDDO651:                            ENDDO
652:                            CALL QORDER_LJ(QTEMP,Q4,Q6)652:                            CALL QORDER_LJ(QTEMP,Q4,Q6)
653:                            ORDER(3)=Q6653:                            ORDER(3)=Q6
654: !                          ORDER(4)=Q6654: !                          ORDER(4)=Q6
675:                      ELSEIF ((AMBERT.OR.NABT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN675:                      ELSEIF ((AMBERT.OR.NABT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN
676:                         IF (I1.EQ.1) CALL AMBERDIHEDR(Q,NATOMS,5,7,9,15,ORDERPLUS)676:                         IF (I1.EQ.1) CALL AMBERDIHEDR(Q,NATOMS,5,7,9,15,ORDERPLUS)
677:                         IF (I1.EQ.2) CALL AMBERDIHEDR(Q,NATOMS,7,9,15,17,ORDERPLUS)677:                         IF (I1.EQ.2) CALL AMBERDIHEDR(Q,NATOMS,7,9,15,17,ORDERPLUS)
678:                      ELSE 678:                      ELSE 
679: !679: !
680: ! Fix the permutation to save time and avoid discontinuities.680: ! Fix the permutation to save time and avoid discontinuities.
681: !681: !
682:                         PERMDISTLOCAL=PERMDIST682:                         PERMDISTLOCAL=PERMDIST
683:                         PERMDIST=.FALSE.683:                         PERMDIST=.FALSE.
684:                         IF (I1.EQ.1) THEN684:                         IF (I1.EQ.1) THEN
685:                            CALL ALIGN_DECIDE(POINTSDECA,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERPLUS,DIST2,685:                            CALL MINPERMDIST(POINTSDECA,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERPLUS,DIST2,
686:      &                                                RIGIDBODY,RMAT)686:      &                                                RIGIDBODY,RMAT)
687:                         ELSE IF (I1.EQ.2) THEN687:                         ELSE IF (I1.EQ.2) THEN
688:                            CALL ALIGN_DECIDE(POINTSICOS,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERPLUS,DIST2,688:                            CALL MINPERMDIST(POINTSICOS,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERPLUS,DIST2,
689:      &                                                RIGIDBODY,RMAT)689:      &                                                RIGIDBODY,RMAT)
690:                         ELSE IF (I1.GT.2) THEN690:                         ELSE IF (I1.GT.2) THEN
691:                            DO K1=1,3691:                            DO K1=1,3
692:                               DO K2=1,NATOMS692:                               DO K2=1,NATOMS
693:                                  QTEMP(K1,K2)=Q(3*(K2-1)+K1)693:                                  QTEMP(K1,K2)=Q(3*(K2-1)+K1)
694:                               ENDDO694:                               ENDDO
695:                            ENDDO695:                            ENDDO
696:                            CALL QORDER_LJ(QTEMP,Q4,Q6)696:                            CALL QORDER_LJ(QTEMP,Q4,Q6)
697:                            IF (I1.EQ.3) ORDERPLUS=Q6697:                            IF (I1.EQ.3) ORDERPLUS=Q6
698: !                          IF (I1.EQ.4) ORDERPLUS=Q6698: !                          IF (I1.EQ.4) ORDERPLUS=Q6
709:                      ELSEIF ((AMBERT.OR.NABT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN709:                      ELSEIF ((AMBERT.OR.NABT).AND.(WHICHORDER(I1).EQ.'DIHE')) THEN
710:                         IF (I1.EQ.1) CALL AMBERDIHEDR(Q,NATOMS,5,7,9,15,ORDERMINUS)710:                         IF (I1.EQ.1) CALL AMBERDIHEDR(Q,NATOMS,5,7,9,15,ORDERMINUS)
711:                         IF (I1.EQ.2) CALL AMBERDIHEDR(Q,NATOMS,7,9,15,17,ORDERMINUS)711:                         IF (I1.EQ.2) CALL AMBERDIHEDR(Q,NATOMS,7,9,15,17,ORDERMINUS)
712:                      ELSE712:                      ELSE
713: !713: !
714: ! Fix the permutation to save time and avoid discontinuities.714: ! Fix the permutation to save time and avoid discontinuities.
715: !715: !
716:                         PERMDISTLOCAL=PERMDIST716:                         PERMDISTLOCAL=PERMDIST
717:                         PERMDIST=.FALSE.717:                         PERMDIST=.FALSE.
718:                         IF (I1.EQ.1) THEN718:                         IF (I1.EQ.1) THEN
719:                            CALL ALIGN_DECIDE(POINTSDECA,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERMINUS,DIST2,719:                            CALL MINPERMDIST(POINTSDECA,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERMINUS,DIST2,
720:      &                                                RIGIDBODY,RMAT)720:      &                                                RIGIDBODY,RMAT)
721:                         ELSE IF (I1.EQ.2) THEN721:                         ELSE IF (I1.EQ.2) THEN
722:                            CALL ALIGN_DECIDE(POINTSICOS,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERMINUS,DIST2,722:                            CALL MINPERMDIST(POINTSICOS,Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ORDERMINUS,DIST2,
723:      &                                                RIGIDBODY,RMAT)723:      &                                                RIGIDBODY,RMAT)
724:                         ELSE IF (I1.GT.2) THEN724:                         ELSE IF (I1.GT.2) THEN
725:                            DO K1=1,3725:                            DO K1=1,3
726:                               DO K2=1,NATOMS726:                               DO K2=1,NATOMS
727:                                  QTEMP(K1,K2)=Q(3*(K2-1)+K1)727:                                  QTEMP(K1,K2)=Q(3*(K2-1)+K1)
728:                               ENDDO728:                               ENDDO
729:                            ENDDO729:                            ENDDO
730:                            CALL QORDER_LJ(QTEMP,Q4,Q6)730:                            CALL QORDER_LJ(QTEMP,Q4,Q6)
731:                            IF (I1.EQ.3) ORDERMINUS=Q6731:                            IF (I1.EQ.3) ORDERMINUS=Q6
732: !                          IF (I1.EQ.4) ORDERMINUS=Q6732: !                          IF (I1.EQ.4) ORDERMINUS=Q6
812:                ENDDO812:                ENDDO
813:                DEALLOCATE(ORDER,ORDERDERIV)813:                DEALLOCATE(ORDER,ORDERDERIV)
814:                IDONE=IDONE+1814:                IDONE=IDONE+1
815:                IF ((NEXMODES.EQ.7).AND.(IDONE.LT.1000)) THEN815:                IF ((NEXMODES.EQ.7).AND.(IDONE.LT.1000)) THEN
816:                   IF (PATHT) THEN816:                   IF (PATHT) THEN
817:                      READ(9124,*,END=444) 817:                      READ(9124,*,END=444) 
818:                      READ(9124,*,END=444) 818:                      READ(9124,*,END=444) 
819:                      DO J2=1,NATOMS819:                      DO J2=1,NATOMS
820:                         READ(9124,*,END=444) ADUMMY,Q(3*(J2-1)+1),Q(3*(J2-1)+2),Q(3*(J2-1)+3)820:                         READ(9124,*,END=444) ADUMMY,Q(3*(J2-1)+1),Q(3*(J2-1)+2),Q(3*(J2-1)+3)
821:                      ENDDO821:                      ENDDO
822:                      CALL ALIGN_DECIDE(Q,QSAVE,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)822:                      CALL MINPERMDIST(Q,QSAVE,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
823:                      DUMMY=0.D0823:                      DUMMY=0.D0
824:                      DO J2=1,NATOMS824:                      DO J2=1,NATOMS
825:                         DUMMY=DUMMY+ATMASS(J2)*( (Q(3*(J2-1)+1)-QSAVE(3*(J2-1)+1))**2 825:                         DUMMY=DUMMY+ATMASS(J2)*( (Q(3*(J2-1)+1)-QSAVE(3*(J2-1)+1))**2 
826:      &                                          +(Q(3*(J2-1)+2)-QSAVE(3*(J2-1)+2))**2826:      &                                          +(Q(3*(J2-1)+2)-QSAVE(3*(J2-1)+2))**2
827:      &                                          +(Q(3*(J2-1)+3)-QSAVE(3*(J2-1)+3))**2 )827:      &                                          +(Q(3*(J2-1)+3)-QSAVE(3*(J2-1)+3))**2 )
828: C                        DUMMY=DUMMY+( (Q(3*(J2-1)+1)-QSAVE(3*(J2-1)+1))**2828: C                        DUMMY=DUMMY+( (Q(3*(J2-1)+1)-QSAVE(3*(J2-1)+1))**2
829: C     &                                          +(Q(3*(J2-1)+2)-QSAVE(3*(J2-1)+2))**2829: C     &                                          +(Q(3*(J2-1)+2)-QSAVE(3*(J2-1)+2))**2
830: C     &                                          +(Q(3*(J2-1)+3)-QSAVE(3*(J2-1)+3))**2 )830: C     &                                          +(Q(3*(J2-1)+3)-QSAVE(3*(J2-1)+3))**2 )
831:                      ENDDO831:                      ENDDO
832:                      TSDISP=SQRT(DUMMY)832:                      TSDISP=SQRT(DUMMY)
1856: C compared to the analytical solution significantly (around 4sf), relative to the alternative1856: C compared to the analytical solution significantly (around 4sf), relative to the alternative
1857: C of evaluating the gradient once at the beginning then once for each element. But obviously it also1857: C of evaluating the gradient once at the beginning then once for each element. But obviously it also
1858: C doubles the number of potential evaluations ((3N)^2/2 for N atoms), so may be not worth it for1858: C doubles the number of potential evaluations ((3N)^2/2 for N atoms), so may be not worth it for
1859: C a long run on a big system.1859: C a long run on a big system.
1860: C1860: C
1861:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))1861:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))
1862: !      SBMFIRSTDD = .TRUE.1862: !      SBMFIRSTDD = .TRUE.
1863:       DO I1=1,NOPT1863:       DO I1=1,NOPT
1864: !        write(*,*) 'working on atom', I11864: !        write(*,*) 'working on atom', I1
1865:          CALL FLUSH(6)1865:          CALL FLUSH(6)
1866:          IF (VARIABLES) THEN1866:          IF (VARIABLES.AND.FROZEN(I1)) THEN
1867:             IF (FROZEN(I1)) THEN1867:             DO J1=I1,NOPT
1868:                DO J1=I1,NOPT1868:                HESS(I1,J1)=0.0D0
1869:                   HESS(I1,J1)=0.0D01869:                HESS(J1,I1)=0.0D0
1870:                   HESS(J1,I1)=0.0D01870:             ENDDO
1871:                ENDDO 
1872:             ENDIF 
1873:          ELSE IF ((.NOT.VARIABLES).AND.FROZEN((I1-1)/3+1)) THEN1871:          ELSE IF ((.NOT.VARIABLES).AND.FROZEN((I1-1)/3+1)) THEN
1874:             DO J1=I1,NOPT1872:             DO J1=I1,NOPT
1875:                HESS(I1,J1)=0.0D01873:                HESS(I1,J1)=0.0D0
1876:                HESS(J1,I1)=0.0D01874:                HESS(J1,I1)=0.0D0
1877:             ENDDO1875:             ENDDO
1878:          ELSE1876:          ELSE
1879:             DUM(I1)=X(I1)-DELTA1877:             DUM(I1)=X(I1)-DELTA
1880: !            SBMHESSATOM=int(I1/3)1878: !            SBMHESSATOM=int(I1/3)
1881:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1879:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1882:             DUM(I1)=X(I1)+DELTA1880:             DUM(I1)=X(I1)+DELTA


r33372/greatcirc.f 2017-10-04 18:30:17.652310275 +0100 r33371/greatcirc.f 2017-10-04 18:30:22.304371642 +0100
 67:          PRINT '(A,I10,A,I10,A)', 'ERROR, dimension of W=',SIZE(W,1), 67:          PRINT '(A,I10,A,I10,A)', 'ERROR, dimension of W=',SIZE(W,1),
 68:      $        ' but N*(2*M+1)+2*M=',N*(2*M+1)+2*M,' in mylbfgs' 68:      $        ' but N*(2*M+1)+2*M=',N*(2*M+1)+2*M,' in mylbfgs'
 69:          STOP 69:          STOP
 70:       ENDIF 70:       ENDIF
 71:  71: 
 72:       IF (N.NE.3*NATOMS+1) THEN 72:       IF (N.NE.3*NATOMS+1) THEN
 73:          PRINT*,'ERROR - N and 3*NATOMS are different in GCLBFGS: ',N,3*NATOMS 73:          PRINT*,'ERROR - N and 3*NATOMS are different in GCLBFGS: ',N,3*NATOMS
 74:          STOP 74:          STOP
 75:       ENDIF 75:       ENDIF
 76:  76: 
 77:       WRITE(*,*) "sn402: Changed GCLBFGS so that it calls ALIGN_DECIDE instead of MINPERMDIST" 77:       CALL MINPERMDIST(RS,RF,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTF,DIST2,RIGIDBODY,RMAT)
 78:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
 79:       WRITE(*,*) "Please set FASTOVERLAP and check carefully that this part of the code is working as you expect." 
 80:       WRITE(*,*) "Then remove these messages!" 
 81:       CALL ALIGN_DECIDE(RS,RF,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTF,DIST2,RIGIDBODY,RMAT) 
 82:  78: 
 83:       ALPHA=1.0D0 79:       ALPHA=1.0D0
 84:       NFAIL=0 80:       NFAIL=0
 85:       FRAME=1 81:       FRAME=1
 86:       IF (RESET) ITER=0 82:       IF (RESET) ITER=0
 87:       ITDONE=0 83:       ITDONE=0
 88:       IF (RESET.AND.PTEST) WRITE(*,'(A)') ' Resetting LBFGS minimiser' 84:       IF (RESET.AND.PTEST) WRITE(*,'(A)') ' Resetting LBFGS minimiser'
 89:       IF ((.NOT.RESET).AND.PTEST) WRITE(*,'(A)') 85:       IF ((.NOT.RESET).AND.PTEST) WRITE(*,'(A)')
 90:      $     ' Not resetting LBFGS minimiser' 86:      $     ' Not resetting LBFGS minimiser'
 91: 1     FIXIMAGE=.FALSE. 87: 1     FIXIMAGE=.FALSE.


r33372/intlbfgs.f90 2017-10-04 18:30:17.872313178 +0100 r33371/intlbfgs.f90 2017-10-04 18:30:22.528374597 +0100
1051:    ENDIF1051:    ENDIF
1052: !1052: !
1053: ! End of add/subtract images block.1053: ! End of add/subtract images block.
1054: !1054: !
1055: !1055: !
1056: ! The new QCILPERMDIST check should be used instead of these lines1056: ! The new QCILPERMDIST check should be used instead of these lines
1057: !1057: !
1058: !  IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN1058: !  IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN
1059: !     LDEBUG=.FALSE.1059: !     LDEBUG=.FALSE.
1060: !     DO J2=2,INTIMAGE+21060: !     DO J2=2,INTIMAGE+2
1061: !        CALL ALIGN_DECIDE(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &1061: !        CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
1062: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)1062: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
1063: !     ENDDO1063: !     ENDDO
1064: !  ENDIF1064: !  ENDIF
1065: 1065: 
1066:    IF (.NOT.SWITCHED) THEN1066:    IF (.NOT.SWITCHED) THEN
1067:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1067:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1068:       IF (QCIADDREP.GT.0) THEN1068:       IF (QCIADDREP.GT.0) THEN
1069:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1069:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1070:       ELSEIF (CHECKCONINT) THEN1070:       ELSEIF (CHECKCONINT) THEN
1071:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1071:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2385:       ENDDO2385:       ENDDO
2386:       NCONSTRAINT=J22386:       NCONSTRAINT=J2
2387:       WRITE(*,'(A,I6,A)') ' checkperc> After allowing for frozen atoms there are ',NCONSTRAINT,' constraints'2387:       WRITE(*,'(A,I6,A)') ' checkperc> After allowing for frozen atoms there are ',NCONSTRAINT,' constraints'
2388:       RETURN 2388:       RETURN 
2389:    ELSE2389:    ELSE
2390: !2390: !
2391: ! Put reference minima in optimal permutational alignment with reference minimum one.2391: ! Put reference minima in optimal permutational alignment with reference minimum one.
2392: !2392: !
2393:       DO J2=2,NCONGEOM2393:       DO J2=2,NCONGEOM
2394:          LDEBUG=.FALSE.2394:          LDEBUG=.FALSE.
2395:          CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),CONGEOM(J2,1:3*NATOMS),NATOMS,LDEBUG, &2395:          CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),CONGEOM(J2,1:3*NATOMS),NATOMS,LDEBUG, &
2396:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)2396:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
2397:       ENDDO2397:       ENDDO
2398:    ENDIF2398:    ENDIF
2399:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))2399:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))
2400: ENDIF2400: ENDIF
2401: 2401: 
2402: INQUIRE(FILE='constraintfile',EXIST=CONFILET)2402: INQUIRE(FILE='constraintfile',EXIST=CONFILET)
2403: 2403: 
2404: 51   NCONSTRAINT=0 2404: 51   NCONSTRAINT=0 
2405: MAXCONDIST=-1.0D02405: MAXCONDIST=-1.0D0


r33372/make_conpot.f90 2017-10-04 18:30:18.096316132 +0100 r33371/make_conpot.f90 2017-10-04 18:30:22.760377655 +0100
 41: ! 41: !
 42: ! If this is not the first call, and we are being passed two minima, 42: ! If this is not the first call, and we are being passed two minima,
 43: ! then we are doing an interpolation metric for a new pair of minima. 43: ! then we are doing an interpolation metric for a new pair of minima.
 44: ! We should optimise the permutational isomers on reference minimum 1 44: ! We should optimise the permutational isomers on reference minimum 1
 45: ! and then do the overall alignment with newmindist, fixing the 45: ! and then do the overall alignment with newmindist, fixing the
 46: ! permutational isomers. This should put the permutational isomers 46: ! permutational isomers. This should put the permutational isomers
 47: ! in register with the constraints, which were calculated for all 47: ! in register with the constraints, which were calculated for all
 48: ! the reference minima after aligning with the first. 48: ! the reference minima after aligning with the first.
 49: ! 49: !
 50:    IF ((CALLED.OR.CONDATT).AND.(NCPFIT.EQ.2)) THEN 50:    IF ((CALLED.OR.CONDATT).AND.(NCPFIT.EQ.2)) THEN
 51:       CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),MINCOORDS(1,1:3*NATOMS),NATOMS,DEBUG, & 51:       CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),MINCOORDS(1,1:3*NATOMS),NATOMS,DEBUG, &
 52:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT) 52:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT)
 53:       CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DEBUG, & 53:       CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DEBUG, &
 54:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT) 54:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT)
 55:       CALL NEWMINDIST(MINCOORDS(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DISTANCE, & 55:       CALL NEWMINDIST(MINCOORDS(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DISTANCE, &
 56:   &                   BULKT,TWOD,'AX    ',.FALSE.,RIGIDBODY,DEBUG,RMAT) 56:   &                   BULKT,TWOD,'AX    ',.FALSE.,RIGIDBODY,DEBUG,RMAT)
 57:    ENDIF 57:    ENDIF
 58: ENDIF 58: ENDIF
 59:  59: 
 60: NQCIFREEZE=0 60: NQCIFREEZE=0
 61: IF (FREEZE) THEN 61: IF (FREEZE) THEN
 62:    WRITE(*, '(A)') ' make_conpot> ERROR *** QCI has not been coded for frozen atoms yet' 62:    WRITE(*, '(A)') ' make_conpot> ERROR *** QCI has not been coded for frozen atoms yet'
 63:    STOP 63:    STOP


r33372/mcpath2.f90 2017-10-04 18:30:18.544322044 +0100 r33371/mcpath2.f90 2017-10-04 18:30:23.236383936 +0100
165: CLOSE(LUNIT)165: CLOSE(LUNIT)
166: !166: !
167: ! Just in case things aren't permutationally aligned and should be, let's check this here.167: ! Just in case things aren't permutationally aligned and should be, let's check this here.
168: !168: !
169: DUMMY2=-1.0D0169: DUMMY2=-1.0D0
170: SLENGTH(1)=0.0D0170: SLENGTH(1)=0.0D0
171: DO J1=2,NPATH171: DO J1=2,NPATH
172:    NTRIES=0172:    NTRIES=0
173: 863 CONTINUE173: 863 CONTINUE
174:    NTRIES=NTRIES+1174:    NTRIES=NTRIES+1
175:    CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &175:    CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
176:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)176:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
177:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1,' after alignment is ',DUMMY177:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1,' after alignment is ',DUMMY
178:    IF (DEBUG) PRINT '(3(A,I6),A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1, &178:    IF (DEBUG) PRINT '(3(A,I6),A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1, &
179:   &                 ' after alignment attempt ',NTRIES,' is ',DUMMY179:   &                 ' after alignment attempt ',NTRIES,' is ',DUMMY
180:    IF ((DUMMY.GT.2.0D0*MAXBFGS).AND.(NTRIES.LT.10)) GOTO 863180:    IF ((DUMMY.GT.2.0D0*MAXBFGS).AND.(NTRIES.LT.10)) GOTO 863
181:    PRINT '(3(A,I6),A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1,' after alignment attempt ',NTRIES,' is ',DUMMY181:    PRINT '(3(A,I6),A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1,' after alignment attempt ',NTRIES,' is ',DUMMY
182:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY182:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY
183:    IF (DUMMY.GT.1.0D0) THEN183:    IF (DUMMY.GT.1.0D0) THEN
184:       PRINT '(A,I6)','frame ',J1-1184:       PRINT '(A,I6)','frame ',J1-1
185:       PRINT '(3G20.10)',PATHFRAMES(J1-1,1:3*NATOMS)185:       PRINT '(3G20.10)',PATHFRAMES(J1-1,1:3*NATOMS)
192: PRINT '(A,G20.10)','mcpath2> largest distance between frames in path.xyz is now',DUMMY2192: PRINT '(A,G20.10)','mcpath2> largest distance between frames in path.xyz is now',DUMMY2
193: !193: !
194: ! Should now be safe to turn permutations off!194: ! Should now be safe to turn permutations off!
195: !195: !
196: PERMDISTSAVE=.FALSE. 196: PERMDISTSAVE=.FALSE. 
197: LPERMDISTSAVE=.FALSE.197: LPERMDISTSAVE=.FALSE.
198: PERMDIST=.FALSE. 198: PERMDIST=.FALSE. 
199: LPERMDIST=.FALSE.199: LPERMDIST=.FALSE.
200: DUMMY2=-1.0D0200: DUMMY2=-1.0D0
201: DO J1=2,NPATH201: DO J1=2,NPATH
202:    CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &202:    CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
203:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)203:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
204:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1,' is now ',DUMMY204:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath2> Distance between frames ',J1,' and ',J1-1,' is now ',DUMMY
205:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY205:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY
206: ENDDO206: ENDDO
207: PRINT '(A,G20.10)','mcpath2> check: largest distance between frames in path.xyz is ',DUMMY2207: PRINT '(A,G20.10)','mcpath2> check: largest distance between frames in path.xyz is ',DUMMY2
208: 208: 
209: IF (DUMMY2.GT.1.0D0) STOP !!! debug DJW209: IF (DUMMY2.GT.1.0D0) STOP !!! debug DJW
210: 210: 
211: NSTRUCTREF=NPATH211: NSTRUCTREF=NPATH
212: 212: 
263:             DO K2=1,NATOMS263:             DO K2=1,NATOMS
264:                QTEMP(K1,K2)=COORDSREF(J1,3*(K2-1)+K1)264:                QTEMP(K1,K2)=COORDSREF(J1,3*(K2-1)+K1)
265:             ENDDO265:             ENDDO
266:          ENDDO266:          ENDDO
267:          CALL GETORDER(QTEMP,QORDER,NATOMS)267:          CALL GETORDER(QTEMP,QORDER,NATOMS)
268:          PRINT '(A,G20.10)','mcpath2> Q=',QORDER268:          PRINT '(A,G20.10)','mcpath2> Q=',QORDER
269:       ENDIF269:       ENDIF
270:    ENDIF270:    ENDIF
271:    IF (EOFS(J1-1)>EOFS(J1).AND.EOFS(J1)<=EOFS(J1+1)) THEN271:    IF (EOFS(J1-1)>EOFS(J1).AND.EOFS(J1)<=EOFS(J1+1)) THEN
272:       IF (VREF(J1).GT.VREF(HIGHESTREF)) HIGHESTREF=J1272:       IF (VREF(J1).GT.VREF(HIGHESTREF)) HIGHESTREF=J1
273:       CALL ALIGN_DECIDE(COORDSREF(J1-1,1:3*NATOMS),COORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG, &273:       CALL MINPERMDIST(COORDSREF(J1-1,1:3*NATOMS),COORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG, &
274:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)274:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
275:       PRINT '(A,I10,2(A,G20.10))','mcpath2> minimum identified for frame ',J1,' energy=', &275:       PRINT '(A,I10,2(A,G20.10))','mcpath2> minimum identified for frame ',J1,' energy=', &
276:   &         VREF(J1)276:   &         VREF(J1)
277:       IF (.NOT.LASTTS) THEN277:       IF (.NOT.LASTTS) THEN
278:          PRINT '(A)','mcpath2> ERROR *** previous stationary point identified was not a transition state'278:          PRINT '(A)','mcpath2> ERROR *** previous stationary point identified was not a transition state'
279:          STOP279:          STOP
280:       ENDIF280:       ENDIF
281:       LASTTS=.FALSE.281:       LASTTS=.FALSE.
282:       LASTMIN=.TRUE.282:       LASTMIN=.TRUE.
283:       IF (.TRUE.) THEN283:       IF (.TRUE.) THEN
284:          DO K1=1,3284:          DO K1=1,3
285:             DO K2=1,NATOMS285:             DO K2=1,NATOMS
286:                QTEMP(K1,K2)=COORDSREF(J1,3*(K2-1)+K1)286:                QTEMP(K1,K2)=COORDSREF(J1,3*(K2-1)+K1)
287:             ENDDO287:             ENDDO
288:          ENDDO288:          ENDDO
289:          CALL GETORDER(QTEMP,QORDER,NATOMS)289:          CALL GETORDER(QTEMP,QORDER,NATOMS)
290:          PRINT '(A,G20.10)','mcpath2> Q=',QORDER290:          PRINT '(A,G20.10)','mcpath2> Q=',QORDER
291:       ENDIF291:       ENDIF
292:    ENDIF292:    ENDIF
293: ENDDO293: ENDDO
294: CALL ALIGN_DECIDE(COORDSREF(NPATH-1,1:3*NATOMS),COORDSREF(NPATH,1:3*NATOMS),NATOMS,DEBUG, &294: CALL MINPERMDIST(COORDSREF(NPATH-1,1:3*NATOMS),COORDSREF(NPATH,1:3*NATOMS),NATOMS,DEBUG, &
295:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)295:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
296:       PRINT '(A,I10,2(A,G20.10))','mcpath2> minimum assumed for frame ',NPATH,' energy=', &296:       PRINT '(A,I10,2(A,G20.10))','mcpath2> minimum assumed for frame ',NPATH,' energy=', &
297:   &         VREF(NPATH)297:   &         VREF(NPATH)
298: !298: !
299: ! This interval will depend upon the order parameter. 299: ! This interval will depend upon the order parameter. 
300: !300: !
301: QORDERHISTINT=(MAX(MCPATHQMAX,MCPATHQMIN)-MIN(MCPATHQMAX,MCPATHQMIN))/MCPATHBINS 301: QORDERHISTINT=(MAX(MCPATHQMAX,MCPATHQMIN)-MIN(MCPATHQMAX,MCPATHQMIN))/MCPATHBINS 
302: DO J1=1,MCPATHBINS302: DO J1=1,MCPATHBINS
303: ! these Q values correspond to the bottom of the Q bins303: ! these Q values correspond to the bottom of the Q bins
304:    BINLABELQORDER(J1)=MIN(MCPATHQMAX,MCPATHQMIN)+QORDERHISTINT*(J1-1.0D0) 304:    BINLABELQORDER(J1)=MIN(MCPATHQMAX,MCPATHQMIN)+QORDERHISTINT*(J1-1.0D0) 
317:    GOTO 963317:    GOTO 963
318: ENDIF318: ENDIF
319: 319: 
320: PRINT '(4(A,I10))','mcpath2> Total number of reference structures=',NSTRUCTREF320: PRINT '(4(A,I10))','mcpath2> Total number of reference structures=',NSTRUCTREF
321: CALL FLUSH(6)321: CALL FLUSH(6)
322: !322: !
323: ! Align path frames if necessary. Should not be executed, since PERMDISTSAVE and LPERMDISTSAVE323: ! Align path frames if necessary. Should not be executed, since PERMDISTSAVE and LPERMDISTSAVE
324: ! were turned off above!324: ! were turned off above!
325: !325: !
326: IF (PERMDISTSAVE.OR.LPERMDISTSAVE) THEN326: IF (PERMDISTSAVE.OR.LPERMDISTSAVE) THEN
327:    CALL ALIGN_DECIDE(COORDSREF(1,1:3*NATOMS),PATHFRAMES(1,1:3*NATOMS),NATOMS,DEBUG, &327:    CALL MINPERMDIST(COORDSREF(1,1:3*NATOMS),PATHFRAMES(1,1:3*NATOMS),NATOMS,DEBUG, &
328:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)328:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
329:    PRINT '(A,F20.10)','mcpath2> first distance=',DUMMY329:    PRINT '(A,F20.10)','mcpath2> first distance=',DUMMY
330:    DO J1=2,NPATH330:    DO J1=2,NPATH
331:       CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &331:       CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
332:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)332:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
333:       PRINT '(A,I10,F20.10)','mcpath2> frame and distance=',J1,DUMMY333:       PRINT '(A,I10,F20.10)','mcpath2> frame and distance=',J1,DUMMY
334:    ENDDO334:    ENDDO
335: !335: !
336: ! Need to resave COORDSREF to get consistent permutational alignment.336: ! Need to resave COORDSREF to get consistent permutational alignment.
337: !337: !
338:    DO J1=1,NSTRUCTREF 338:    DO J1=1,NSTRUCTREF 
339:       COORDSREF(J1,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)339:       COORDSREF(J1,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)
340:    ENDDO340:    ENDDO
341:    PERMDIST=.FALSE.341:    PERMDIST=.FALSE.
1100: ! Check that we are within MCPATHDMAX of at least one reference geometry.1100: ! Check that we are within MCPATHDMAX of at least one reference geometry.
1101: ! If not, reject and recount. We also reject steps outside the current block.1101: ! If not, reject and recount. We also reject steps outside the current block.
1102: !1102: !
1103:   RECOUNT=.TRUE.1103:   RECOUNT=.TRUE.
1104:   DISTMIN=1.0D1001104:   DISTMIN=1.0D100
1105:   IMAGEMIN=11105:   IMAGEMIN=1
1106:   DISTMAX=-1.0D1001106:   DISTMAX=-1.0D100
1107:   IMAGEMAX=11107:   IMAGEMAX=1
1108: 1108: 
1109: ! DO J1=1,LSTRUCTREF1109: ! DO J1=1,LSTRUCTREF
1110: !    CALL ALIGN_DECIDE(MCCOORDS,LCOORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1110: !    CALL MINPERMDIST(MCCOORDS,LCOORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1111: ! &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)1111: ! &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)
1112: !    IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY31112: !    IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY3
1113: !    IF (DUMMY3.LT.DISTMIN) THEN1113: !    IF (DUMMY3.LT.DISTMIN) THEN
1114: !       DISTMIN=DUMMY31114: !       DISTMIN=DUMMY3
1115: !       IMAGEMIN=J11115: !       IMAGEMIN=J1
1116: !    ENDIF1116: !    ENDIF
1117: !    IF (DUMMY3.GT.DISTMAX) THEN1117: !    IF (DUMMY3.GT.DISTMAX) THEN
1118: !       DISTMAX=DUMMY31118: !       DISTMAX=DUMMY3
1119: !       IMAGEMAX=J11119: !       IMAGEMAX=J1
1120: !    ENDIF1120: !    ENDIF
1121: ! ENDDO1121: ! ENDDO
1122: 1122: 
1123: !1123: !
1124: ! First try previous closest frame and neighbours.1124: ! First try previous closest frame and neighbours.
1125: !1125: !
1126: ! IF (IMCSTEP.GT.1) THEN1126: ! IF (IMCSTEP.GT.1) THEN
1127: !    DO J1=MAX(MAX(1,STARTSAMPLE-MCPATHBLOCK+MCPATHOVER),IMAGEMINO+(STARTSAMPLE-1)-1), &1127: !    DO J1=MAX(MAX(1,STARTSAMPLE-MCPATHBLOCK+MCPATHOVER),IMAGEMINO+(STARTSAMPLE-1)-1), &
1128: ! &        MIN(MIN(NPATH,ENDSAMPLE+MCPATHBLOCK-MCPATHOVER),IMAGEMINO+(STARTSAMPLE-1)+1)1128: ! &        MIN(MIN(NPATH,ENDSAMPLE+MCPATHBLOCK-MCPATHOVER),IMAGEMINO+(STARTSAMPLE-1)+1)
1129: !       CALL ALIGN_DECIDE(MCCOORDS,COORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1129: !       CALL MINPERMDIST(MCCOORDS,COORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1130: ! &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)1130: ! &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)
1131: !       IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY31131: !       IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY3
1132: !       PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY31132: !       PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY3
1133: !       IF (DUMMY3.LT.DISTMIN) THEN1133: !       IF (DUMMY3.LT.DISTMIN) THEN
1134: !          DISTMIN=DUMMY31134: !          DISTMIN=DUMMY3
1135: !          IMAGEMIN=J11135: !          IMAGEMIN=J1
1136: !       ENDIF1136: !       ENDIF
1137: !       IF (DUMMY3.GT.DISTMAX) THEN1137: !       IF (DUMMY3.GT.DISTMAX) THEN
1138: !          DISTMAX=DUMMY31138: !          DISTMAX=DUMMY3
1139: !          IMAGEMAX=J11139: !          IMAGEMAX=J1
1141: !    ENDDO1141: !    ENDDO
1142: ! ENDIF1142: ! ENDIF
1143: !1143: !
1144: ! Check all frames every NCHECKINT steps.1144: ! Check all frames every NCHECKINT steps.
1145: !1145: !
1146:   NCHECKINT=11146:   NCHECKINT=1
1147:   IF ((MOD(IMCSTEP,NCHECKINT).EQ.0).OR.(IMCSTEP.EQ.1)) THEN1147:   IF ((MOD(IMCSTEP,NCHECKINT).EQ.0).OR.(IMCSTEP.EQ.1)) THEN
1148:      DISTMIN=1.0D1001148:      DISTMIN=1.0D100
1149:      DISTMAX=-1.0D1001149:      DISTMAX=-1.0D100
1150:      DO J1=MAX(1,STARTSAMPLE-MCPATHBLOCK+MCPATHOVER),MIN(NPATH,ENDSAMPLE+MCPATHBLOCK-MCPATHOVER)1150:      DO J1=MAX(1,STARTSAMPLE-MCPATHBLOCK+MCPATHOVER),MIN(NPATH,ENDSAMPLE+MCPATHBLOCK-MCPATHOVER)
1151:         CALL ALIGN_DECIDE(MCCOORDS,COORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1151:         CALL MINPERMDIST(MCCOORDS,COORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1152:   &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)1152:   &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)
1153:         IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY31153:         IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY3
1154: !       PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY31154: !       PRINT '(A,I10,A,G20.10)','mcpath2> distance for local reference structure ',J1,' is ',DUMMY3
1155:         IF (DUMMY3.LT.DISTMIN) THEN1155:         IF (DUMMY3.LT.DISTMIN) THEN
1156:            DISTMIN=DUMMY31156:            DISTMIN=DUMMY3
1157:            IMAGEMIN2=J11157:            IMAGEMIN2=J1
1158:         ENDIF1158:         ENDIF
1159:         IF (DUMMY3.GT.DISTMAX) THEN1159:         IF (DUMMY3.GT.DISTMAX) THEN
1160:            DISTMAX=DUMMY31160:            DISTMAX=DUMMY3
1161:            IMAGEMAX=J11161:            IMAGEMAX=J1


r33372/mcpath.f90 2017-10-04 18:30:18.320319087 +0100 r33371/mcpath.f90 2017-10-04 18:30:23.000380824 +0100
130: CLOSE(LUNIT)130: CLOSE(LUNIT)
131: !131: !
132: ! Just in case things aren't permutationally aligned and should be, let's check this here.132: ! Just in case things aren't permutationally aligned and should be, let's check this here.
133: !133: !
134: DUMMY2=-1.0D0134: DUMMY2=-1.0D0
135: SLENGTH(1)=0.0D0135: SLENGTH(1)=0.0D0
136: DO J1=2,NPATH136: DO J1=2,NPATH
137:    NTRIES=0137:    NTRIES=0
138: 163 CONTINUE138: 163 CONTINUE
139:    NTRIES=NTRIES+1139:    NTRIES=NTRIES+1
140:    CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &140:    CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
141:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)141:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
142:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1,' after alignment is ',DUMMY142:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1,' after alignment is ',DUMMY
143:    IF (DEBUG) PRINT '(3(A,I6),A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1, &143:    IF (DEBUG) PRINT '(3(A,I6),A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1, &
144:   &                 ' after alignment attempt ',NTRIES,' is ',DUMMY144:   &                 ' after alignment attempt ',NTRIES,' is ',DUMMY
145:    IF ((DUMMY.GT.2.0D0*MAXBFGS).AND.(NTRIES.LT.10000)) GOTO 163145:    IF ((DUMMY.GT.2.0D0*MAXBFGS).AND.(NTRIES.LT.10000)) GOTO 163
146:    PRINT '(3(A,I6),A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1,' after alignment attempt ',NTRIES,' is ',DUMMY146:    PRINT '(3(A,I6),A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1,' after alignment attempt ',NTRIES,' is ',DUMMY
147:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY147:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY
148:    IF (DUMMY.GT.1.0D0) THEN148:    IF (DUMMY.GT.1.0D0) THEN
149:       PRINT '(A,I6)','frame ',J1-1149:       PRINT '(A,I6)','frame ',J1-1
150:       PRINT '(3G20.10)',PATHFRAMES(J1-1,1:3*NATOMS)150:       PRINT '(3G20.10)',PATHFRAMES(J1-1,1:3*NATOMS)
157: PRINT '(A,G20.10)','mcpath> largest distance between frames in path.xyz is now',DUMMY2157: PRINT '(A,G20.10)','mcpath> largest distance between frames in path.xyz is now',DUMMY2
158: !158: !
159: ! Should now be safe to turn permutations off!159: ! Should now be safe to turn permutations off!
160: !160: !
161: PERMDISTSAVE=.FALSE. 161: PERMDISTSAVE=.FALSE. 
162: LPERMDISTSAVE=.FALSE.162: LPERMDISTSAVE=.FALSE.
163: PERMDIST=.FALSE. 163: PERMDIST=.FALSE. 
164: LPERMDIST=.FALSE.164: LPERMDIST=.FALSE.
165: DUMMY2=-1.0D0165: DUMMY2=-1.0D0
166: DO J1=2,NPATH166: DO J1=2,NPATH
167:    CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &167:    CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
168:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)168:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
169:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1,' is now ',DUMMY169:    IF (DEBUG) PRINT '(A,I6,A,I6,A,G20.10)','mcpath> Distance between frames ',J1,' and ',J1-1,' is now ',DUMMY
170:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY170:    IF (DUMMY.GT.DUMMY2) DUMMY2=DUMMY
171: ENDDO171: ENDDO
172: PRINT '(A,G20.10)','mcpath> largest distance between frames in path.xyz is now',DUMMY2172: PRINT '(A,G20.10)','mcpath> largest distance between frames in path.xyz is now',DUMMY2
173: 173: 
174: IF (DUMMY2.GT.1.0D0) STOP !!! debug DJW174: IF (DUMMY2.GT.1.0D0) STOP !!! debug DJW
175: 175: 
176: NSTRUCTREF=2176: NSTRUCTREF=2
177: DO J1=2,NPATH-1177: DO J1=2,NPATH-1
217:    PRINT '(A,G20.10)','mcpath> Q=',QORDER217:    PRINT '(A,G20.10)','mcpath> Q=',QORDER
218: ENDIF218: ENDIF
219: 219: 
220: DO J1=2,NPATH-1220: DO J1=2,NPATH-1
221:    IF (EOFS(J1-1)+EDIFFTOL*10.0D0<EOFS(J1).AND.EOFS(J1)>EOFS(J1+1)+EDIFFTOL*10.0D0) THEN221:    IF (EOFS(J1-1)+EDIFFTOL*10.0D0<EOFS(J1).AND.EOFS(J1)>EOFS(J1+1)+EDIFFTOL*10.0D0) THEN
222:       NSTRUCTREF=NSTRUCTREF+1222:       NSTRUCTREF=NSTRUCTREF+1
223:       REFSTRUCTFRAME(NSTRUCTREF)=J1223:       REFSTRUCTFRAME(NSTRUCTREF)=J1
224:       VREF(NSTRUCTREF)=EOFS(J1) 224:       VREF(NSTRUCTREF)=EOFS(J1) 
225:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF225:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF
226:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)226:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)
227:       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &227:       CALL MINPERMDIST(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &
228:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)228:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
229:       PRINT '(A,I10,2(A,G20.10))','mcpath> transition state identified for frame ',J1,' energy=', &229:       PRINT '(A,I10,2(A,G20.10))','mcpath> transition state identified for frame ',J1,' energy=', &
230:   &         VREF(NSTRUCTREF),' distance to previous minimum=',DUMMY230:   &         VREF(NSTRUCTREF),' distance to previous minimum=',DUMMY
231:       NTS=NTS+1231:       NTS=NTS+1
232:       IF (.NOT.LASTMIN) THEN232:       IF (.NOT.LASTMIN) THEN
233:          PRINT '(A)','mcpath> ERROR *** previous stationary point identified was not a minimum'233:          PRINT '(A)','mcpath> ERROR *** previous stationary point identified was not a minimum'
234:          STOP234:          STOP
235:       ENDIF235:       ENDIF
236:       LASTMIN=.FALSE.236:       LASTMIN=.FALSE.
237:       LASTTS=.TRUE.237:       LASTTS=.TRUE.
244:          CALL GETORDER(QTEMP,QORDER,NATOMS)244:          CALL GETORDER(QTEMP,QORDER,NATOMS)
245:          PRINT '(A,G20.10)','mcpath> Q=',QORDER245:          PRINT '(A,G20.10)','mcpath> Q=',QORDER
246:       ENDIF246:       ENDIF
247:    ENDIF247:    ENDIF
248:    IF ((EOFS(J1-1)>EOFS(J1).AND.EOFS(J1)<=EOFS(J1+1)).AND.(ABS(EOFS(J1)-VREF(NSTRUCTREF)).GT.EDIFFTOL)) THEN248:    IF ((EOFS(J1-1)>EOFS(J1).AND.EOFS(J1)<=EOFS(J1+1)).AND.(ABS(EOFS(J1)-VREF(NSTRUCTREF)).GT.EDIFFTOL)) THEN
249:       NSTRUCTREF=NSTRUCTREF+1249:       NSTRUCTREF=NSTRUCTREF+1
250:       REFSTRUCTFRAME(NSTRUCTREF)=J1250:       REFSTRUCTFRAME(NSTRUCTREF)=J1
251:       VREF(NSTRUCTREF)=EOFS(J1) 251:       VREF(NSTRUCTREF)=EOFS(J1) 
252:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF252:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF
253:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)253:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)
254:       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &254:       CALL MINPERMDIST(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &
255:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)255:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
256:       PRINT '(A,I10,2(A,G20.10))','mcpath> minimum identified for frame ',J1,' energy=', &256:       PRINT '(A,I10,2(A,G20.10))','mcpath> minimum identified for frame ',J1,' energy=', &
257:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY257:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY
258:       NMIN=NMIN+1258:       NMIN=NMIN+1
259:       IF (.NOT.LASTTS) THEN259:       IF (.NOT.LASTTS) THEN
260:          PRINT '(A)','mcpath> ERROR *** previous stationary point identified was not a transition state'260:          PRINT '(A)','mcpath> ERROR *** previous stationary point identified was not a transition state'
261:          STOP261:          STOP
262:       ENDIF262:       ENDIF
263:       LASTTS=.FALSE.263:       LASTTS=.FALSE.
264:       LASTMIN=.TRUE.264:       LASTMIN=.TRUE.
297: !    PERMDISTSAVE=.FALSE. 297: !    PERMDISTSAVE=.FALSE. 
298: !    LPERMDISTSAVE=.FALSE.298: !    LPERMDISTSAVE=.FALSE.
299: !    PERMDIST=.FALSE. 299: !    PERMDIST=.FALSE. 
300: !    LPERMDIST=.FALSE.300: !    LPERMDIST=.FALSE.
301: ! ENDIF301: ! ENDIF
302: IF (ABS(EOFS(NPATH)-VREF(NSTRUCTREF)).GT.EDIFFTOL) THEN302: IF (ABS(EOFS(NPATH)-VREF(NSTRUCTREF)).GT.EDIFFTOL) THEN
303:    NSTRUCTREF=NSTRUCTREF+1303:    NSTRUCTREF=NSTRUCTREF+1
304:    REFSTRUCTFRAME(NSTRUCTREF)=NPATH304:    REFSTRUCTFRAME(NSTRUCTREF)=NPATH
305:    VREF(NSTRUCTREF)=EOFS(NPATH) ! last minimum305:    VREF(NSTRUCTREF)=EOFS(NPATH) ! last minimum
306:    COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(NPATH,1:3*NATOMS)306:    COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(NPATH,1:3*NATOMS)
307:       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &307:       CALL MINPERMDIST(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &
308:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)308:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
309:    PRINT '(A,I10,2(A,G20.10))','mcpath> minimum assumed    for frame ',NPATH,' energy=', &309:    PRINT '(A,I10,2(A,G20.10))','mcpath> minimum assumed    for frame ',NPATH,' energy=', &
310:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY310:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY
311:    NMIN=NMIN+1311:    NMIN=NMIN+1
312:    IF (.NOT.LASTTS) THEN312:    IF (.NOT.LASTTS) THEN
313:       PRINT '(A)','mcpath> ERROR *** previous stationary point identified was not a transition state'313:       PRINT '(A)','mcpath> ERROR *** previous stationary point identified was not a transition state'
314:       STOP314:       STOP
315:    ENDIF315:    ENDIF
316:    IF (.TRUE.) THEN316:    IF (.TRUE.) THEN
317:       DO K1=1,3317:       DO K1=1,3
332: !332: !
333: ! Check whether any of the reference structures are separated by distances greater than MCPATHDMAX333: ! Check whether any of the reference structures are separated by distances greater than MCPATHDMAX
334: ! If so, add more references! We might want to tie the mcpath simulation closer to the path than334: ! If so, add more references! We might want to tie the mcpath simulation closer to the path than
335: ! the separation between two stationary points.335: ! the separation between two stationary points.
336: !336: !
337: 654 CONTINUE337: 654 CONTINUE
338: DO J1=1,NSTRUCTREF338: DO J1=1,NSTRUCTREF
339: !339: !
340: ! Check distance to original frame is zero!340: ! Check distance to original frame is zero!
341: !341: !
342:    CALL ALIGN_DECIDE(COORDSREF(J1,1:3*NATOMS),PATHFRAMES(REFSTRUCTFRAME(J1),1:3*NATOMS),NATOMS,DEBUG, &342:    CALL MINPERMDIST(COORDSREF(J1,1:3*NATOMS),PATHFRAMES(REFSTRUCTFRAME(J1),1:3*NATOMS),NATOMS,DEBUG, &
343:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)343:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
344:    PRINT '(A,I6,A,I6,A,F20.10)','mcpath> distance between reference structure ',J1,' and corresponding path frame ', &344:    PRINT '(A,I6,A,I6,A,F20.10)','mcpath> distance between reference structure ',J1,' and corresponding path frame ', &
345:   &                 REFSTRUCTFRAME(J1),' is ',DUMMY345:   &                 REFSTRUCTFRAME(J1),' is ',DUMMY
346: !346: !
347: ! Find closest references forwards and backwards.347: ! Find closest references forwards and backwards.
348: !348: !
349:    DPLUS=1.0D100349:    DPLUS=1.0D100
350:    DMINUS=1.0D100350:    DMINUS=1.0D100
351:    DO J2=1,NSTRUCTREF351:    DO J2=1,NSTRUCTREF
352:       IF (J2.EQ.J1) CYCLE352:       IF (J2.EQ.J1) CYCLE
353: 353: 
354:       CALL ALIGN_DECIDE(COORDSREF(J1,1:3*NATOMS),COORDSREF(J2,1:3*NATOMS),NATOMS,DEBUG, &354:       CALL MINPERMDIST(COORDSREF(J1,1:3*NATOMS),COORDSREF(J2,1:3*NATOMS),NATOMS,DEBUG, &
355:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)355:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
356:       IF (DEBUG) PRINT '(2(A,I10),A,G20.10,A,I6,A,I6)','mcpath> distance between reference structures ',J1,' and ',J2,' is ',DUMMY, &356:       IF (DEBUG) PRINT '(2(A,I10),A,G20.10,A,I6,A,I6)','mcpath> distance between reference structures ',J1,' and ',J2,' is ',DUMMY, &
357:   &   ' frames: ',REFSTRUCTFRAME(J1),' and ',REFSTRUCTFRAME(J2)357:   &   ' frames: ',REFSTRUCTFRAME(J1),' and ',REFSTRUCTFRAME(J2)
358:       IF (REFSTRUCTFRAME(J1).LT.REFSTRUCTFRAME(J2)) THEN358:       IF (REFSTRUCTFRAME(J1).LT.REFSTRUCTFRAME(J2)) THEN
359:          IF (DUMMY.LT.DPLUS) THEN359:          IF (DUMMY.LT.DPLUS) THEN
360:             NCLOSEP=J2360:             NCLOSEP=J2
361:             DPLUS=DUMMY361:             DPLUS=DUMMY
362:          ENDIF362:          ENDIF
363:       ELSE363:       ELSE
364:          IF (DUMMY.LT.DMINUS) THEN364:          IF (DUMMY.LT.DMINUS) THEN
378:    IF (DPLUS.GT.MCPATHADDREF) THEN ! add another reference structure between J1 and NCLOSEP378:    IF (DPLUS.GT.MCPATHADDREF) THEN ! add another reference structure between J1 and NCLOSEP
379:       DUMMY2=1.0D100379:       DUMMY2=1.0D100
380:       NDUMMY=-1380:       NDUMMY=-1
381:       rloopp: DO J2=MIN(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEP))+1,MAX(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEP))-1 381:       rloopp: DO J2=MIN(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEP))+1,MAX(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEP))-1 
382: !382: !
383: ! Check whether this frame is already a reference.383: ! Check whether this frame is already a reference.
384: !384: !
385:          DO J3=1,NSTRUCTREF385:          DO J3=1,NSTRUCTREF
386:             IF (REFSTRUCTFRAME(J3).EQ.J2) CYCLE rloopp386:             IF (REFSTRUCTFRAME(J3).EQ.J2) CYCLE rloopp
387:          ENDDO387:          ENDDO
388:          CALL ALIGN_DECIDE(COORDSREF(J1,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &388:          CALL MINPERMDIST(COORDSREF(J1,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &
389:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)389:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
390:          DUMMY3=DUMMY390:          DUMMY3=DUMMY
391:          CALL ALIGN_DECIDE(COORDSREF(NCLOSEP,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &391:          CALL MINPERMDIST(COORDSREF(NCLOSEP,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &
392:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)392:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
393:          DUMMY3=ABS(DUMMY3-DUMMY)393:          DUMMY3=ABS(DUMMY3-DUMMY)
394: !        PRINT '(A,G20.10,A,I6)','mcpath> Difference of distances is ',DUMMY3,' for frame ',J2394: !        PRINT '(A,G20.10,A,I6)','mcpath> Difference of distances is ',DUMMY3,' for frame ',J2
395:          IF (DUMMY3.LT.DUMMY2) THEN395:          IF (DUMMY3.LT.DUMMY2) THEN
396:             DUMMY2=DUMMY3396:             DUMMY2=DUMMY3
397:             NDUMMY=J2397:             NDUMMY=J2
398:             IF (DEBUG) PRINT '(A,G20.10,A,I6)','mcpath> Smallest difference of distances so far is ',DUMMY2,' for frame ',NDUMMY398:             IF (DEBUG) PRINT '(A,G20.10,A,I6)','mcpath> Smallest difference of distances so far is ',DUMMY2,' for frame ',NDUMMY
399:          ENDIF399:          ENDIF
400:       ENDDO rloopp400:       ENDDO rloopp
401: 401: 
446:    IF (DMINUS.GT.MCPATHADDREF) THEN ! add another reference structure between J1 and NCLOSEM446:    IF (DMINUS.GT.MCPATHADDREF) THEN ! add another reference structure between J1 and NCLOSEM
447:       DUMMY2=1.0D100447:       DUMMY2=1.0D100
448:       NDUMMY=-1448:       NDUMMY=-1
449:       rloopm: DO J2=MIN(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEM))+1,MAX(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEM))-1 449:       rloopm: DO J2=MIN(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEM))+1,MAX(REFSTRUCTFRAME(J1),REFSTRUCTFRAME(NCLOSEM))-1 
450: !450: !
451: ! Check whether this frame is already a reference.451: ! Check whether this frame is already a reference.
452: !452: !
453:          DO J3=1,NSTRUCTREF453:          DO J3=1,NSTRUCTREF
454:             IF (REFSTRUCTFRAME(J3).EQ.J2) CYCLE rloopm454:             IF (REFSTRUCTFRAME(J3).EQ.J2) CYCLE rloopm
455:          ENDDO455:          ENDDO
456:          CALL ALIGN_DECIDE(COORDSREF(J1,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &456:          CALL MINPERMDIST(COORDSREF(J1,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &
457:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)457:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
458:          DUMMY3=DUMMY458:          DUMMY3=DUMMY
459:          CALL ALIGN_DECIDE(COORDSREF(NCLOSEM,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &459:          CALL MINPERMDIST(COORDSREF(NCLOSEM,1:3*NATOMS),PATHFRAMES(J2,1:3*NATOMS),NATOMS,DEBUG, &
460:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)460:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
461:          DUMMY3=ABS(DUMMY3-DUMMY)461:          DUMMY3=ABS(DUMMY3-DUMMY)
462: !        PRINT '(A,G20.10,A,I6)','mcpath> Difference of distances is ',DUMMY3,' for frame ',J2462: !        PRINT '(A,G20.10,A,I6)','mcpath> Difference of distances is ',DUMMY3,' for frame ',J2
463:          IF (DUMMY3.LT.DUMMY2) THEN463:          IF (DUMMY3.LT.DUMMY2) THEN
464:             DUMMY2=DUMMY3464:             DUMMY2=DUMMY3
465:             NDUMMY=J2465:             NDUMMY=J2
466:             IF (DEBUG) PRINT '(A,G20.10,A,I6)','mcpath> Smallest difference of distances so far is ',DUMMY2,' for frame ',NDUMMY466:             IF (DEBUG) PRINT '(A,G20.10,A,I6)','mcpath> Smallest difference of distances so far is ',DUMMY2,' for frame ',NDUMMY
467:          ENDIF467:          ENDIF
468:       ENDDO rloopm468:       ENDDO rloopm
469: 469: 
516:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(NDUMMY,1:3*NATOMS)516:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(NDUMMY,1:3*NATOMS)
517:       GOTO 654517:       GOTO 654
518:    ENDIF518:    ENDIF
519: ENDDO519: ENDDO
520: PRINT '(A,I6)','mcpath> Total number of reference structures is now: ',NSTRUCTREF520: PRINT '(A,I6)','mcpath> Total number of reference structures is now: ',NSTRUCTREF
521: !521: !
522: ! Align path frames if necessary. Should not be executed, since PERMDISTSAVE and LPERMDISTSAVE522: ! Align path frames if necessary. Should not be executed, since PERMDISTSAVE and LPERMDISTSAVE
523: ! were turned off above!523: ! were turned off above!
524: !524: !
525: IF (PERMDISTSAVE.OR.LPERMDISTSAVE) THEN525: IF (PERMDISTSAVE.OR.LPERMDISTSAVE) THEN
526:    CALL ALIGN_DECIDE(COORDSREF(1,1:3*NATOMS),PATHFRAMES(1,1:3*NATOMS),NATOMS,DEBUG, &526:    CALL MINPERMDIST(COORDSREF(1,1:3*NATOMS),PATHFRAMES(1,1:3*NATOMS),NATOMS,DEBUG, &
527:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)527:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
528:    PRINT '(A,F20.10)','mcpath> first distance=',DUMMY528:    PRINT '(A,F20.10)','mcpath> first distance=',DUMMY
529:    DO J1=2,NPATH529:    DO J1=2,NPATH
530:       CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &530:       CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
531:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)531:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
532:       PRINT '(A,I10,F20.10)','mcpath> frame and distance=',J1,DUMMY532:       PRINT '(A,I10,F20.10)','mcpath> frame and distance=',J1,DUMMY
533:    ENDDO533:    ENDDO
534: !534: !
535: ! Need to resave COORDSREF to get consistent permutational alignment.535: ! Need to resave COORDSREF to get consistent permutational alignment.
536: !536: !
537:    DO J1=1,NSTRUCTREF 537:    DO J1=1,NSTRUCTREF 
538:       COORDSREF(J1,1:3*NATOMS)=PATHFRAMES(REFSTRUCTFRAME(J1),1:3*NATOMS)538:       COORDSREF(J1,1:3*NATOMS)=PATHFRAMES(REFSTRUCTFRAME(J1),1:3*NATOMS)
539:    ENDDO539:    ENDDO
540:    PERMDIST=.FALSE.540:    PERMDIST=.FALSE.
547: ! This was not working!! DISTFRAME was never set.547: ! This was not working!! DISTFRAME was never set.
548: !548: !
549: IF (.NOT.ALLOCATED(DISTFRAME))   ALLOCATE(DISTFRAME(NPATH,NSTRUCTREF))549: IF (.NOT.ALLOCATED(DISTFRAME))   ALLOCATE(DISTFRAME(NPATH,NSTRUCTREF))
550: !550: !
551: ! Shepard interpolation using just the pe at551: ! Shepard interpolation using just the pe at
552: ! each reference structure and the distance.552: ! each reference structure and the distance.
553: !553: !
554: ! IF (NSTRUCTREF.EQ.NMIN+NTS) THEN554: ! IF (NSTRUCTREF.EQ.NMIN+NTS) THEN
555:    DO J2=1,NSTRUCTREF555:    DO J2=1,NSTRUCTREF
556:       DO J1=1,NPATH556:       DO J1=1,NPATH
557:          CALL ALIGN_DECIDE(COORDSREF(J2,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &557:          CALL MINPERMDIST(COORDSREF(J2,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
558:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTFRAME(J1,J2),DIST2,RIGIDBODY,RMAT)558:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTFRAME(J1,J2),DIST2,RIGIDBODY,RMAT)
559:       ENDDO559:       ENDDO
560:    ENDDO560:    ENDDO
561: ! ELSE561: ! ELSE
562: !    DO J1=1,NPATH562: !    DO J1=1,NPATH
563: !       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &563: !       CALL MINPERMDIST(COORDSREF(NSTRUCTREF,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
564: !   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTFRAME(J1,NSTRUCTREF),DIST2,RIGIDBODY,RMAT)564: !   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTFRAME(J1,NSTRUCTREF),DIST2,RIGIDBODY,RMAT)
565: !    ENDDO565: !    ENDDO
566: ! ENDIF566: ! ENDIF
567: 567: 
568: BIASSTEP=1.0D0568: BIASSTEP=1.0D0
569: CALL SUMSQ(BIASFAC,DISTFRAME,NPATH,NSTRUCTREF,SUM1,VREF,EOFS,MAXDEV,IMAX)569: CALL SUMSQ(BIASFAC,DISTFRAME,NPATH,NSTRUCTREF,SUM1,VREF,EOFS,MAXDEV,IMAX)
570: BIASBEST=BIASFAC570: BIASBEST=BIASFAC
571: SUMBEST=SUM1571: SUMBEST=SUM1
572: NDUMMY=0572: NDUMMY=0
573:  DO WHILE (NDUMMY<51)573:  DO WHILE (NDUMMY<51)
655: ! However, we could reoptimise the bias exponent for different655: ! However, we could reoptimise the bias exponent for different
656: ! ts regions.656: ! ts regions.
657: ! Visualise new bias potential along the reaction pathway.657: ! Visualise new bias potential along the reaction pathway.
658: !658: !
659: OPEN(LUNIT,FILE='Bias.new',STATUS='UNKNOWN')659: OPEN(LUNIT,FILE='Bias.new',STATUS='UNKNOWN')
660: DO J1=1,NPATH 660: DO J1=1,NPATH 
661:    DUMMY=0.0D0661:    DUMMY=0.0D0
662:    DUMMY2=0.0D0662:    DUMMY2=0.0D0
663:    DMIN=1.0D100663:    DMIN=1.0D100
664:    DO J2=1,NSTRUCTREF664:    DO J2=1,NSTRUCTREF
665:       CALL ALIGN_DECIDE(COORDSREF(J2,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &665:       CALL MINPERMDIST(COORDSREF(J2,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
666:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)666:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)
667:       IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)667:       IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)
668:      668:      
669:    ENDDO669:    ENDDO
670:    DO J2=1,NSTRUCTREF670:    DO J2=1,NSTRUCTREF
671:       DUMMY=DUMMY+VREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))671:       DUMMY=DUMMY+VREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
672:       DUMMY2=DUMMY2+       EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))672:       DUMMY2=DUMMY2+       EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
673: !     PRINT '(A,2I6,5G15.5)','J1,J2, VREF,DSTRUCT,DUMMY,DUMMY2,DMIN=',J1,J2,VREF(J2),DSTRUCT(J2),DUMMY,DUMMY2,DMIN673: !     PRINT '(A,2I6,5G15.5)','J1,J2, VREF,DSTRUCT,DUMMY,DUMMY2,DMIN=',J1,J2,VREF(J2),DSTRUCT(J2),DUMMY,DUMMY2,DMIN
674:    ENDDO674:    ENDDO
675: 675: 
777:    LVREF(3)=VREF(ENDSAMPLE)777:    LVREF(3)=VREF(ENDSAMPLE)
778:    LREFSTRUCTFRAME(3)=REFSTRUCTFRAME(ENDSAMPLE)778:    LREFSTRUCTFRAME(3)=REFSTRUCTFRAME(ENDSAMPLE)
779:    DO K1=1,3779:    DO K1=1,3
780:       DO K2=1,NATOMS780:       DO K2=1,NATOMS
781:          QTEMP(K1,K2)=LCOORDSREF(3,3*(K2-1)+K1)781:          QTEMP(K1,K2)=LCOORDSREF(3,3*(K2-1)+K1)
782:       ENDDO782:       ENDDO
783:    ENDDO783:    ENDDO
784:    CALL GETORDER(QTEMP,QORDER,NATOMS)784:    CALL GETORDER(QTEMP,QORDER,NATOMS)
785:    PRINT '(3(A,G20.10),A,I6)','mcpath> Local reference 3 energy=',DUMMY,' should be ',VREF(ENDSAMPLE),' Q=',QORDER, &785:    PRINT '(3(A,G20.10),A,I6)','mcpath> Local reference 3 energy=',DUMMY,' should be ',VREF(ENDSAMPLE),' Q=',QORDER, &
786:   &                        ' frame=',REFSTRUCTFRAME(ENDSAMPLE)786:   &                        ' frame=',REFSTRUCTFRAME(ENDSAMPLE)
787:          CALL ALIGN_DECIDE(LCOORDSREF(1,1:3*NATOMS),LCOORDSREF(2,1:3*NATOMS),NATOMS,DEBUG, &787:          CALL MINPERMDIST(LCOORDSREF(1,1:3*NATOMS),LCOORDSREF(2,1:3*NATOMS),NATOMS,DEBUG, &
788:   &                  PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)788:   &                  PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
789:    PRINT '(A,G20.10)','minimum distance between reference structures 1 and 2=',DUMMY789:    PRINT '(A,G20.10)','minimum distance between reference structures 1 and 2=',DUMMY
790:       CALL ALIGN_DECIDE(LCOORDSREF(3,1:3*NATOMS),LCOORDSREF(2,1:3*NATOMS),NATOMS,DEBUG, &790:       CALL MINPERMDIST(LCOORDSREF(3,1:3*NATOMS),LCOORDSREF(2,1:3*NATOMS),NATOMS,DEBUG, &
791:   &                  PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)791:   &                  PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
792:    PRINT '(A,G20.10)','minimum distance between reference structures 3 and 2=',DUMMY792:    PRINT '(A,G20.10)','minimum distance between reference structures 3 and 2=',DUMMY
793:    LSTRUCTREF=3793:    LSTRUCTREF=3
794:    LNPATH=REFSTRUCTFRAME(ENDSAMPLE)-REFSTRUCTFRAME(STARTSAMPLE)+1794:    LNPATH=REFSTRUCTFRAME(ENDSAMPLE)-REFSTRUCTFRAME(STARTSAMPLE)+1
795:    PRINT '(A,I6)','mcpath> Setting local number of path frames to ',LNPATH795:    PRINT '(A,I6)','mcpath> Setting local number of path frames to ',LNPATH
796:    IF (NSTRUCTREF.GT.NMIN+NTS) THEN796:    IF (NSTRUCTREF.GT.NMIN+NTS) THEN
797:       DO J1=1,NSTRUCTREF797:       DO J1=1,NSTRUCTREF
798:          IF (J1.EQ.2*DOTS) CYCLE ! don't add the transition state again!798:          IF (J1.EQ.2*DOTS) CYCLE ! don't add the transition state again!
799:          IF ((REFSTRUCTFRAME(J1).LT.REFSTRUCTFRAME(ENDSAMPLE)).AND.(REFSTRUCTFRAME(J1).GT.REFSTRUCTFRAME(STARTSAMPLE))) THEN799:          IF ((REFSTRUCTFRAME(J1).LT.REFSTRUCTFRAME(ENDSAMPLE)).AND.(REFSTRUCTFRAME(J1).GT.REFSTRUCTFRAME(STARTSAMPLE))) THEN
800:             LSTRUCTREF=LSTRUCTREF+1800:             LSTRUCTREF=LSTRUCTREF+1
874: ! IF (MCPATHTS.EQ.0) THEN874: ! IF (MCPATHTS.EQ.0) THEN
875: IF ((NTS.GT.1).AND.(MCPATHTS.GE.0)) THEN875: IF ((NTS.GT.1).AND.(MCPATHTS.GE.0)) THEN
876: !876: !
877: ! We need to create LDISTFRAME for use in BIASFAC reoptimisation877: ! We need to create LDISTFRAME for use in BIASFAC reoptimisation
878: !878: !
879:    ALLOCATE(LDISTFRAME(LNPATH,LSTRUCTREF))879:    ALLOCATE(LDISTFRAME(LNPATH,LSTRUCTREF))
880:    ALLOCATE(LEOFS(LNPATH))880:    ALLOCATE(LEOFS(LNPATH))
881:    LEOFS(1:LNPATH)=EOFS(REFSTRUCTFRAME(STARTSAMPLE):REFSTRUCTFRAME(ENDSAMPLE))881:    LEOFS(1:LNPATH)=EOFS(REFSTRUCTFRAME(STARTSAMPLE):REFSTRUCTFRAME(ENDSAMPLE))
882:    DO J1=1,LNPATH882:    DO J1=1,LNPATH
883:       DO J2=1,LSTRUCTREF883:       DO J2=1,LSTRUCTREF
884:          CALL ALIGN_DECIDE(LCOORDSREF(J2,1:3*NATOMS), &884:          CALL MINPERMDIST(LCOORDSREF(J2,1:3*NATOMS), &
885:   &                    PATHFRAMES(J1+REFSTRUCTFRAME(STARTSAMPLE)-1,1:3*NATOMS),NATOMS,DEBUG, &885:   &                    PATHFRAMES(J1+REFSTRUCTFRAME(STARTSAMPLE)-1,1:3*NATOMS),NATOMS,DEBUG, &
886:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,LDISTFRAME(J1,J2),DIST2,RIGIDBODY,RMAT)886:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,LDISTFRAME(J1,J2),DIST2,RIGIDBODY,RMAT)
887:       ENDDO887:       ENDDO
888:    ENDDO888:    ENDDO
889:    BIASSTEP=1.0D0889:    BIASSTEP=1.0D0
890:    BIASSTEP=1.0D0890:    BIASSTEP=1.0D0
891:    CALL SUMSQ(BIASFAC,LDISTFRAME,LNPATH,LSTRUCTREF,SUM1,LVREF,LEOFS,MAXDEV,IMAX)891:    CALL SUMSQ(BIASFAC,LDISTFRAME,LNPATH,LSTRUCTREF,SUM1,LVREF,LEOFS,MAXDEV,IMAX)
892:    BIASBEST=BIASFAC892:    BIASBEST=BIASFAC
893:    SUMBEST=SUM1893:    SUMBEST=SUM1
894:    NDUMMY=0894:    NDUMMY=0
918:    WRITE(FNAME,'(I6)') DOTS918:    WRITE(FNAME,'(I6)') DOTS
919:    FNAME='Bias.new.ts' // TRIM(ADJUSTL(FNAME))919:    FNAME='Bias.new.ts' // TRIM(ADJUSTL(FNAME))
920:    OPEN(LUNIT,FILE=FNAME,STATUS='UNKNOWN')920:    OPEN(LUNIT,FILE=FNAME,STATUS='UNKNOWN')
921:    PRINT '(A,I6)','DOTS=',DOTS921:    PRINT '(A,I6)','DOTS=',DOTS
922:    PRINT '(A,I6)','REFSTRUCTFRAME(STARTSAMPLE)=',REFSTRUCTFRAME(STARTSAMPLE)922:    PRINT '(A,I6)','REFSTRUCTFRAME(STARTSAMPLE)=',REFSTRUCTFRAME(STARTSAMPLE)
923:    DO J1=1,LNPATH923:    DO J1=1,LNPATH
924:       DUMMY=0.0D0924:       DUMMY=0.0D0
925:       DUMMY2=0.0D0925:       DUMMY2=0.0D0
926:       DMIN=1.0D100926:       DMIN=1.0D100
927:       DO J2=1,LSTRUCTREF927:       DO J2=1,LSTRUCTREF
928:          CALL ALIGN_DECIDE(LCOORDSREF(J2,1:3*NATOMS), &928:          CALL MINPERMDIST(LCOORDSREF(J2,1:3*NATOMS), &
929:   &                    PATHFRAMES(J1+REFSTRUCTFRAME(STARTSAMPLE)-1,1:3*NATOMS),NATOMS,DEBUG, &929:   &                    PATHFRAMES(J1+REFSTRUCTFRAME(STARTSAMPLE)-1,1:3*NATOMS),NATOMS,DEBUG, &
930:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)930:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)
931:          IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)931:          IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)
932:       ENDDO932:       ENDDO
933:       DO J2=1,LSTRUCTREF933:       DO J2=1,LSTRUCTREF
934:          DUMMY=DUMMY+LVREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))934:          DUMMY=DUMMY+LVREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
935:          DUMMY2=DUMMY2+        EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))935:          DUMMY2=DUMMY2+        EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
936: !     PRINT '(A,2I6,5G15.5)','J1,J2,LVREF,DSTRUCT,DUMMY,DUMMY2,DMIN=',J1,J2,LVREF(J2),DSTRUCT(J2),DUMMY,DUMMY2,DMIN936: !     PRINT '(A,2I6,5G15.5)','J1,J2,LVREF,DSTRUCT,DUMMY,DUMMY2,DMIN=',J1,J2,LVREF(J2),DSTRUCT(J2),DUMMY,DUMMY2,DMIN
937:       ENDDO937:       ENDDO
938: 938: 
977: ENDIF977: ENDIF
978: 978: 
979: CALL POTENTIAL(MCCOORDS,VOLD,GRAD,.FALSE.,.FALSE.,RMS,.FALSE.,.FALSE.)979: CALL POTENTIAL(MCCOORDS,VOLD,GRAD,.FALSE.,.FALSE.,RMS,.FALSE.,.FALSE.)
980: PRINT '(A,G20.10)','mcpath> Initial configuration energy             = ',VOLD980: PRINT '(A,G20.10)','mcpath> Initial configuration energy             = ',VOLD
981: 981: 
982: IF (MCBIAST) THEN982: IF (MCBIAST) THEN
983:    DUMMY=0.0D0983:    DUMMY=0.0D0
984:    DUMMY2=0.0D0984:    DUMMY2=0.0D0
985:    DMIN=1.0D100985:    DMIN=1.0D100
986:    DO J2=1,LSTRUCTREF986:    DO J2=1,LSTRUCTREF
987:       CALL ALIGN_DECIDE(LCOORDSREF(J2,1:3*NATOMS),MCCOORDS,NATOMS,DEBUG, &987:       CALL MINPERMDIST(LCOORDSREF(J2,1:3*NATOMS),MCCOORDS,NATOMS,DEBUG, &
988:   &                  PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)988:   &                  PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)
989:       IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)989:       IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)
990:    ENDDO990:    ENDDO
991:    DO J2=1,LSTRUCTREF991:    DO J2=1,LSTRUCTREF
992:       DUMMY=DUMMY+LVREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))992:       DUMMY=DUMMY+LVREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
993:       DUMMY2=DUMMY2+        EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))993:       DUMMY2=DUMMY2+        EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
994: !     PRINT '(A,I6,3G20.10)','J2,DSTRUCT,DUMMY,DUMMY2=',J2,DSTRUCT(J2),DUMMY,DUMMY2994: !     PRINT '(A,I6,3G20.10)','J2,DSTRUCT,DUMMY,DUMMY2=',J2,DSTRUCT(J2),DUMMY,DUMMY2
995: !     PRINT '(A,2G20.10)','MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)=',MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)995: !     PRINT '(A,2G20.10)','MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)=',MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)
996:    ENDDO996:    ENDDO
997:    WOLD=-DUMMY/MAX(1.0D-10,DUMMY2)997:    WOLD=-DUMMY/MAX(1.0D-10,DUMMY2)
1092: ! Check that we are within MCPATHDMAX of at least one reference geometry.1092: ! Check that we are within MCPATHDMAX of at least one reference geometry.
1093: ! If not, reject and recount.1093: ! If not, reject and recount.
1094: !1094: !
1095:   RECOUNT=.TRUE.1095:   RECOUNT=.TRUE.
1096:   DISTMIN=1.0D1001096:   DISTMIN=1.0D100
1097:   IMAGEMIN=11097:   IMAGEMIN=1
1098:   IMAGEMINF=11098:   IMAGEMINF=1
1099:   DISTMAX=-1.0D1001099:   DISTMAX=-1.0D100
1100:   IMAGEMAX=11100:   IMAGEMAX=1
1101:   DO J1=1,LSTRUCTREF1101:   DO J1=1,LSTRUCTREF
1102:      CALL ALIGN_DECIDE(MCCOORDS,LCOORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1102:      CALL MINPERMDIST(MCCOORDS,LCOORDSREF(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1103:   &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)1103:   &       TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)
1104:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for local reference structure ',J1,' is ',DUMMY31104:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for local reference structure ',J1,' is ',DUMMY3
1105:      IF (DUMMY3.LT.DISTMIN) THEN1105:      IF (DUMMY3.LT.DISTMIN) THEN
1106:         DISTMIN=DUMMY31106:         DISTMIN=DUMMY3
1107:         IMAGEMIN=J11107:         IMAGEMIN=J1
1108:      ENDIF1108:      ENDIF
1109:      IF (DUMMY3.GT.DISTMAX) THEN1109:      IF (DUMMY3.GT.DISTMAX) THEN
1110:         DISTMAX=DUMMY31110:         DISTMAX=DUMMY3
1111:         IMAGEMAX=J11111:         IMAGEMAX=J1
1112:      ENDIF1112:      ENDIF
1113:   ENDDO1113:   ENDDO
1114: 1114: 
1115:   IF (IMCSTEP.EQ.1) THEN1115:   IF (IMCSTEP.EQ.1) THEN
1116:      IMAGEMINO=IMAGEMIN ! initialise this information at first step1116:      IMAGEMINO=IMAGEMIN ! initialise this information at first step
1117:      DISTMINF=1.0D1001117:      DISTMINF=1.0D100
1118:      IMAGEMINF=11118:      IMAGEMINF=1
1119:      DO J1=1,NPATH1119:      DO J1=1,NPATH
1120:         CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1120:         CALL MINPERMDIST(MCCOORDS,PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1121:      &       TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)1121:      &       TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
1122:         IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',J1,' is ',DUMMY1122:         IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',J1,' is ',DUMMY
1123:         IF (DUMMY.LT.DISTMINF) THEN1123:         IF (DUMMY.LT.DISTMINF) THEN
1124:            DISTMINF=DUMMY1124:            DISTMINF=DUMMY
1125:            IMAGEMINF=J11125:            IMAGEMINF=J1
1126:         ENDIF1126:         ENDIF
1127:      ENDDO1127:      ENDDO
1128:      PRINT '(A,I6,A,G20.10)','mcpath> Initial configuration is closest to path frame ',IMAGEMINF,' distance=',DISTMINF1128:      PRINT '(A,I6,A,G20.10)','mcpath> Initial configuration is closest to path frame ',IMAGEMINF,' distance=',DISTMINF
1129:      IF ((MCPATHTS.GE.0).AND.((IMAGEMINF.LT.REFSTRUCTFRAME(STARTSAMPLE)).OR.(IMAGEMINF.GT.REFSTRUCTFRAME(ENDSAMPLE)))) THEN1129:      IF ((MCPATHTS.GE.0).AND.((IMAGEMINF.LT.REFSTRUCTFRAME(STARTSAMPLE)).OR.(IMAGEMINF.GT.REFSTRUCTFRAME(ENDSAMPLE)))) THEN
1130:         PRINT '(A)','mcath> ERROR *** Initial configuration lies outside the frames defined by the two bracketing stationary points'1130:         PRINT '(A)','mcath> ERROR *** Initial configuration lies outside the frames defined by the two bracketing stationary points'
1157: 1157: 
1158:         GOTO 9611158:         GOTO 961
1159:      ENDIF1159:      ENDIF
1160:      IMAGEMINFO=IMAGEMINF1160:      IMAGEMINFO=IMAGEMINF
1161:   ELSE1161:   ELSE
1162:      IMAGEMINSTART=IMAGEMINFO1162:      IMAGEMINSTART=IMAGEMINFO
1163:      NDIRECTION=11163:      NDIRECTION=1
1164: ! 1164: ! 
1165: ! Get distance for previous best frame, initialised at IMAGEMINFO1165: ! Get distance for previous best frame, initialised at IMAGEMINFO
1166: !1166: !
1167:      CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(IMAGEMINSTART,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1167:      CALL MINPERMDIST(MCCOORDS,PATHFRAMES(IMAGEMINSTART,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1168:   &    TWOD,DISTMINF,DIST2,RIGIDBODY,RMAT)1168:   &    TWOD,DISTMINF,DIST2,RIGIDBODY,RMAT)
1169: !    IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART,' is ',DISTMINF1169: !    IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART,' is ',DISTMINF
1170: 263 CONTINUE1170: 263 CONTINUE
1171:      IF ((IMAGEMINSTART+NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN1171:      IF ((IMAGEMINSTART+NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN
1172:         CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(IMAGEMINSTART+NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1172:         CALL MINPERMDIST(MCCOORDS,PATHFRAMES(IMAGEMINSTART+NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1173:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)1173:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
1174:      ELSE1174:      ELSE
1175:         DUMMY=1.0D1001175:         DUMMY=1.0D100
1176:      ENDIF1176:      ENDIF
1177:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART+NDIRECTION,' is ',DUMMY1177:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART+NDIRECTION,' is ',DUMMY
1178:      IF (DUMMY.LT.DISTMINF) THEN1178:      IF (DUMMY.LT.DISTMINF) THEN
1179:         DISTMINF=DUMMY1179:         DISTMINF=DUMMY
1180:         IMAGEMINSTART=IMAGEMINSTART+NDIRECTION1180:         IMAGEMINSTART=IMAGEMINSTART+NDIRECTION
1181: !1181: !
1182: ! Closest frame has changed. Need to go back and check again so that we identify the1182: ! Closest frame has changed. Need to go back and check again so that we identify the
1183: ! closest frame consistently and avoid getting stuck! It needs to be a local minimum at least!1183: ! closest frame consistently and avoid getting stuck! It needs to be a local minimum at least!
1184: !1184: !
1185:         GOTO 2631185:         GOTO 263
1186:      ENDIF1186:      ENDIF
1187:      IF ((IMAGEMINSTART-NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN1187:      IF ((IMAGEMINSTART-NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN
1188:         CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(IMAGEMINSTART-NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1188:         CALL MINPERMDIST(MCCOORDS,PATHFRAMES(IMAGEMINSTART-NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1189:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)1189:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
1190:           ELSE1190:           ELSE
1191:         DUMMY=1.0D1001191:         DUMMY=1.0D100
1192:      ENDIF1192:      ENDIF
1193:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART-NDIRECTION,' is ',DUMMY1193:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART-NDIRECTION,' is ',DUMMY
1194:      IF (DUMMY.LT.DISTMINF) THEN1194:      IF (DUMMY.LT.DISTMINF) THEN
1195:         DISTMINF=DUMMY1195:         DISTMINF=DUMMY
1196:         IMAGEMINSTART=IMAGEMINSTART-NDIRECTION1196:         IMAGEMINSTART=IMAGEMINSTART-NDIRECTION
1197:         NDIRECTION=-NDIRECTION1197:         NDIRECTION=-NDIRECTION
1198: !1198: !
1208:      PRINT '(A,I6,A,I6,A,G15.5,A,G15.5,A,I6)','mcpath> WARNING *** Distance to local reference structure ',IMAGEMIN,' frame ', &1208:      PRINT '(A,I6,A,I6,A,G15.5,A,G15.5,A,I6)','mcpath> WARNING *** Distance to local reference structure ',IMAGEMIN,' frame ', &
1209:   &              LREFSTRUCTFRAME(IMAGEMIN),' of ',DISTMIN,' is less than ',DISTMINF,' for current frame ',IMAGEMINF1209:   &              LREFSTRUCTFRAME(IMAGEMIN),' of ',DISTMIN,' is less than ',DISTMINF,' for current frame ',IMAGEMINF
1210: !1210: !
1211: ! Try for a local minimum from this reference structure.1211: ! Try for a local minimum from this reference structure.
1212: !1212: !
1213:      IMAGEMINSTART=LREFSTRUCTFRAME(IMAGEMIN)1213:      IMAGEMINSTART=LREFSTRUCTFRAME(IMAGEMIN)
1214:      NDIRECTION=11214:      NDIRECTION=1
1215:      DISTMINF=DISTMIN1215:      DISTMINF=DISTMIN
1216: 764 CONTINUE1216: 764 CONTINUE
1217:      IF ((IMAGEMINSTART+NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN1217:      IF ((IMAGEMINSTART+NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN
1218:         CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(IMAGEMINSTART+NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1218:         CALL MINPERMDIST(MCCOORDS,PATHFRAMES(IMAGEMINSTART+NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1219:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)1219:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
1220:      ELSE1220:      ELSE
1221:         DUMMY=1.0D1001221:         DUMMY=1.0D100
1222:      ENDIF1222:      ENDIF
1223:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART+NDIRECTION,' is ',DUMMY1223:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART+NDIRECTION,' is ',DUMMY
1224:      IF (DUMMY.LT.DISTMINF) THEN1224:      IF (DUMMY.LT.DISTMINF) THEN
1225:         DISTMINF=DUMMY1225:         DISTMINF=DUMMY
1226:         IMAGEMINSTART=IMAGEMINSTART+NDIRECTION1226:         IMAGEMINSTART=IMAGEMINSTART+NDIRECTION
1227: !1227: !
1228: ! Closest frame has changed. Need to go back and check again so that we identify the1228: ! Closest frame has changed. Need to go back and check again so that we identify the
1229: ! closest frame consistently and avoid getting stuck! It needs to be a local minimum at least!1229: ! closest frame consistently and avoid getting stuck! It needs to be a local minimum at least!
1230: !1230: !
1231:         GOTO 7641231:         GOTO 764
1232:      ENDIF1232:      ENDIF
1233:      IF ((IMAGEMINSTART-NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN1233:      IF ((IMAGEMINSTART-NDIRECTION.GT.0).AND.(IMAGEMINSTART+NDIRECTION.LE.NPATH)) THEN
1234:         CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(IMAGEMINSTART-NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1234:         CALL MINPERMDIST(MCCOORDS,PATHFRAMES(IMAGEMINSTART-NDIRECTION,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1235:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)1235:   &                      TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
1236:           ELSE1236:           ELSE
1237:         DUMMY=1.0D1001237:         DUMMY=1.0D100
1238:      ENDIF1238:      ENDIF
1239:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART-NDIRECTION,' is ',DUMMY1239:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for path frame ',IMAGEMINSTART-NDIRECTION,' is ',DUMMY
1240:      IF (DUMMY.LT.DISTMINF) THEN1240:      IF (DUMMY.LT.DISTMINF) THEN
1241:         DISTMINF=DUMMY1241:         DISTMINF=DUMMY
1242:         IMAGEMINSTART=IMAGEMINSTART-NDIRECTION1242:         IMAGEMINSTART=IMAGEMINSTART-NDIRECTION
1243:         NDIRECTION=-NDIRECTION1243:         NDIRECTION=-NDIRECTION
1244: !1244: !
1245: ! Closest frame has changed. Need to go back and check again so that we identify the1245: ! Closest frame has changed. Need to go back and check again so that we identify the
1246: ! closest frame consistently and avoid getting stuck!1246: ! closest frame consistently and avoid getting stuck!
1247: !1247: !
1248:         GOTO 7641248:         GOTO 764
1249:      ENDIF1249:      ENDIF
1250:      IMAGEMINF=IMAGEMINSTART1250:      IMAGEMINF=IMAGEMINSTART
1251:      PRINT '(A,I6,A,G15.5)','mcpath> WARNING *** Resetting closest frame to ',IMAGEMINF,' distance ',DISTMINF1251:      PRINT '(A,I6,A,G15.5)','mcpath> WARNING *** Resetting closest frame to ',IMAGEMINF,' distance ',DISTMINF
1252: 1252: 
1253: !    CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(LREFSTRUCTFRAME(IMAGEMIN),1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1253: !    CALL MINPERMDIST(MCCOORDS,PATHFRAMES(LREFSTRUCTFRAME(IMAGEMIN),1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1254: ! &    TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)1254: ! &    TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
1255: !    IF (LREFSTRUCTFRAME(IMAGEMIN)-1.GT.0) THEN1255: !    IF (LREFSTRUCTFRAME(IMAGEMIN)-1.GT.0) THEN
1256: !       CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(LREFSTRUCTFRAME(IMAGEMIN)-1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1256: !       CALL MINPERMDIST(MCCOORDS,PATHFRAMES(LREFSTRUCTFRAME(IMAGEMIN)-1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1257: ! &                      TWOD,DUMMY2,DIST2,RIGIDBODY,RMAT)1257: ! &                      TWOD,DUMMY2,DIST2,RIGIDBODY,RMAT)
1258: !    ELSE1258: !    ELSE
1259: !       DUMMY2=1.0D1001259: !       DUMMY2=1.0D100
1260: !    ENDIF1260: !    ENDIF
1261: !    IF (LREFSTRUCTFRAME(IMAGEMIN)+1.LE.NPATH) THEN1261: !    IF (LREFSTRUCTFRAME(IMAGEMIN)+1.LE.NPATH) THEN
1262: !       CALL ALIGN_DECIDE(MCCOORDS,PATHFRAMES(LREFSTRUCTFRAME(IMAGEMIN)+1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1262: !       CALL MINPERMDIST(MCCOORDS,PATHFRAMES(LREFSTRUCTFRAME(IMAGEMIN)+1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1263: ! &                      TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)1263: ! &                      TWOD,DUMMY3,DIST2,RIGIDBODY,RMAT)
1264: !    ELSE1264: !    ELSE
1265: !       DUMMY3=1.0D1001265: !       DUMMY3=1.0D100
1266: !    ENDIF1266: !    ENDIF
1267: !    IF ((DUMMY.LT.DUMMY2).AND.(DUMMY.LT.DUMMY3)) THEN1267: !    IF ((DUMMY.LT.DUMMY2).AND.(DUMMY.LT.DUMMY3)) THEN
1268: !       PRINT '(A,I6,A)','mcpath> WARNING *** Resetting to frame ',LREFSTRUCTFRAME(IMAGEMIN),' for proposed step'1268: !       PRINT '(A,I6,A)','mcpath> WARNING *** Resetting to frame ',LREFSTRUCTFRAME(IMAGEMIN),' for proposed step'
1269: !       DISTMINF=DISTMIN1269: !       DISTMINF=DISTMIN
1270: !       IMAGEMINF=LREFSTRUCTFRAME(IMAGEMIN)1270: !       IMAGEMINF=LREFSTRUCTFRAME(IMAGEMIN)
1271: !    ELSE1271: !    ELSE
1272: !       PRINT '(A)','mcpath> WARNING *** Cannot reset - proposed frame is not a local minimum'1272: !       PRINT '(A)','mcpath> WARNING *** Cannot reset - proposed frame is not a local minimum'
1281:   ENDIF1281:   ENDIF
1282:   IF ((MCPATHTS.GE.0).AND.((IMAGEMINF.LT.REFSTRUCTFRAME(STARTSAMPLE)).OR.(IMAGEMINF.GT.REFSTRUCTFRAME(ENDSAMPLE)))) THEN1282:   IF ((MCPATHTS.GE.0).AND.((IMAGEMINF.LT.REFSTRUCTFRAME(STARTSAMPLE)).OR.(IMAGEMINF.GT.REFSTRUCTFRAME(ENDSAMPLE)))) THEN
1283:      IF (DEBUG) PRINT '(A)','mcpath> reject: closest frame is outside the range defined by the two neighbouring minima'1283:      IF (DEBUG) PRINT '(A)','mcpath> reject: closest frame is outside the range defined by the two neighbouring minima'
1284:      RECOUNT=.TRUE.1284:      RECOUNT=.TRUE.
1285:      IF (IMCSTEP.GT.MCPATHEQUIL) IREJFRAME=IREJFRAME+11285:      IF (IMCSTEP.GT.MCPATHEQUIL) IREJFRAME=IREJFRAME+1
1286:   ENDIF1286:   ENDIF
1287: !1287: !
1288: ! Check rejection condition on neighbouring reference strucures.1288: ! Check rejection condition on neighbouring reference strucures.
1289: !1289: !
1290:   DO J1=1,NCHECKSTRUCT1290:   DO J1=1,NCHECKSTRUCT
1291:      CALL ALIGN_DECIDE(MCCOORDS,LCOORDSCHECK(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &1291:      CALL MINPERMDIST(MCCOORDS,LCOORDSCHECK(J1,1:3*NATOMS),NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT, &
1292:   &       TWOD,DISTCHECK,DIST2,RIGIDBODY,RMAT)1292:   &       TWOD,DISTCHECK,DIST2,RIGIDBODY,RMAT)
1293:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for check reference structure ',J1,' is ',DISTCHECK1293:      IF (DEBUG) PRINT '(A,I10,A,G20.10)','mcpath> distance for check reference structure ',J1,' is ',DISTCHECK
1294:      IF (DISTCHECK.LT.DISTMIN) THEN ! reject step1294:      IF (DISTCHECK.LT.DISTMIN) THEN ! reject step
1295:         IREJCHECK=IREJCHECK+11295:         IREJCHECK=IREJCHECK+1
1296:         RECOUNT=.TRUE.1296:         RECOUNT=.TRUE.
1297:         PRINT '(A)','mcpath> Rejecting step because configuration is nearer to check structure ',J11297:         PRINT '(A)','mcpath> Rejecting step because configuration is nearer to check structure ',J1
1298:         PRINT '(A,F15.5,A,F15.5)','mcpath> Distance is ',DISTCHECK,' compared to ',DISTMIN1298:         PRINT '(A,F15.5,A,F15.5)','mcpath> Distance is ',DISTCHECK,' compared to ',DISTMIN
1299:         EXIT1299:         EXIT
1300:      ENDIF1300:      ENDIF
1301:   ENDDO1301:   ENDDO
1326:      VNEW=VOLD1326:      VNEW=VOLD
1327:      WNEW=WOLD1327:      WNEW=WOLD
1328:   ELSE1328:   ELSE
1329:      CALL POTENTIAL(MCCOORDS,VNEW,GRAD,.FALSE.,.FALSE.,RMS,.FALSE.,.FALSE.)1329:      CALL POTENTIAL(MCCOORDS,VNEW,GRAD,.FALSE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1330: !    PRINT '(A,G20.10)','mcpath> after potential call VNEW=',VNEW1330: !    PRINT '(A,G20.10)','mcpath> after potential call VNEW=',VNEW
1331:      IF (MCBIAST) THEN1331:      IF (MCBIAST) THEN
1332:         DUMMY=0.0D01332:         DUMMY=0.0D0
1333:         DUMMY2=0.0D01333:         DUMMY2=0.0D0
1334:         DMIN=1.0D1001334:         DMIN=1.0D100
1335:         DO J2=1,LSTRUCTREF1335:         DO J2=1,LSTRUCTREF
1336:            CALL ALIGN_DECIDE(LCOORDSREF(J2,1:3*NATOMS),MCCOORDS,NATOMS,DEBUG, &1336:            CALL MINPERMDIST(LCOORDSREF(J2,1:3*NATOMS),MCCOORDS,NATOMS,DEBUG, &
1337:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)1337:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSTRUCT(J2),DIST2,RIGIDBODY,RMAT)
1338:            IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)1338:            IF (DSTRUCT(J2).LT.DMIN) DMIN=DSTRUCT(J2)
1339:         ENDDO1339:         ENDDO
1340:         DO J2=1,LSTRUCTREF1340:         DO J2=1,LSTRUCTREF
1341:            DUMMY=DUMMY+LVREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))1341:            DUMMY=DUMMY+LVREF(J2)*EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
1342:            DUMMY2=DUMMY2+       EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))1342:            DUMMY2=DUMMY2+       EXP(BIASFAC/MAX(1.0D-2,DSTRUCT(J2))-BIASFAC/MAX(1.0D-2,DMIN))
1343: !          PRINT '(A,I6,3G20.10)','J2,DSTRUCT,DUMMY,DUMMY2=',J2,DSTRUCT(J2),DUMMY,DUMMY21343: !          PRINT '(A,I6,3G20.10)','J2,DSTRUCT,DUMMY,DUMMY2=',J2,DSTRUCT(J2),DUMMY,DUMMY2
1344: !          PRINT '(A,2G20.10)','MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)=',MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)1344: !          PRINT '(A,2G20.10)','MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)=',MAX(1.0D-2,DSTRUCT(J2)),MAX(1.0D-2,DMIN)
1345:         ENDDO1345:         ENDDO
1346:         WNEW=-DUMMY/MAX(1.0D-10,DUMMY2)1346:         WNEW=-DUMMY/MAX(1.0D-10,DUMMY2)


r33372/morph.f 2017-10-04 18:30:18.764324944 +0100 r33371/morph.f 2017-10-04 18:30:23.456386838 +0100
 57:       PUSH=.FALSE. 57:       PUSH=.FALSE.
 58:       NPU=0 58:       NPU=0
 59:       FINISHED=.FALSE. 59:       FINISHED=.FALSE.
 60:       MNBFGSMAX1SAVE=MNBFGSMAX1 60:       MNBFGSMAX1SAVE=MNBFGSMAX1
 61:       MNBFGSMAX2SAVE=MNBFGSMAX2 61:       MNBFGSMAX2SAVE=MNBFGSMAX2
 62:  62: 
 63:       OPEN(UNIT=44,FILE='morph.points',STATUS='UNKNOWN') 63:       OPEN(UNIT=44,FILE='morph.points',STATUS='UNKNOWN')
 64:  64: 
 65:       NUP=1 ! mylbfgs uses this for projection of the search direction 65:       NUP=1 ! mylbfgs uses this for projection of the search direction
 66:       QSTART(1:3*NATOMS)=COORDS(1:3*NATOMS) 66:       QSTART(1:3*NATOMS)=COORDS(1:3*NATOMS)
 67:       WRITE(*,*) "sn402: Changed MORPH so that it calls ALIGN_DECIDE instead of MINPERMDIST" 67:       CALL MINPERMDIST(QSTART,QFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTF,DIST2,RIGIDBODY,RMAT)
 68:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
 69:       WRITE(*,*) "Please set FASTOVERLAP and check carefully that this part of the code is working as you expect." 
 70:       WRITE(*,*) "Then remove these messages!" 
 71:       CALL ALIGN_DECIDE(QSTART,QFINISH,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DISTF,DIST2,RIGIDBODY,RMAT) 
 72:       IF (DEBUG) PRINT '(A,F20.10)',' morph> initial distance between start and finish=',DISTF 68:       IF (DEBUG) PRINT '(A,F20.10)',' morph> initial distance between start and finish=',DISTF
 73: !     QSTART(1:3*NATOMS)=COORDS(1:3*NATOMS) 69: !     QSTART(1:3*NATOMS)=COORDS(1:3*NATOMS)
 74:       NREPEL=0 70:       NREPEL=0
 75:       PRINTPTS=.TRUE. 71:       PRINTPTS=.TRUE.
 76:       ITER=1 72:       ITER=1
 77:       PTEST=.FALSE. 73:       PTEST=.FALSE.
 78:       IF (DEBUG) PTEST=.TRUE. 74:       IF (DEBUG) PTEST=.TRUE.
 79: 90    CONTINUE 75: 90    CONTINUE
 80:       IF (POTCALL) THEN 76:       IF (POTCALL) THEN
 81:          IF (PV) THEN 77:          IF (PV) THEN


r33372/ncutils.f90 2017-10-04 18:30:15.172277560 +0100 r33371/ncutils.f90 2017-10-04 18:30:19.472334283 +0100
 35:  35: 
 36:           IF (NTS==0) THEN 36:           IF (NTS==0) THEN
 37:                ISNEWTS=.TRUE. 37:                ISNEWTS=.TRUE.
 38:                RETURN 38:                RETURN
 39:           ENDIF 39:           ENDIF
 40:  40: 
 41:           DO I=1,NTS 41:           DO I=1,NTS
 42:                IF(ABS(TSTOCHECK%E-TS(I)%DATA%E) < EDIFFTOL) THEN 42:                IF(ABS(TSTOCHECK%E-TS(I)%DATA%E) < EDIFFTOL) THEN
 43:                     !print *, "Energy of the TS found is the same as of the TS #",i,"; checking distance..." 43:                     !print *, "Energy of the TS found is the same as of the TS #",i,"; checking distance..."
 44:                     IF (PERMDIST) THEN 44:                     IF (PERMDIST) THEN
 45:                        CALL ALIGN_DECIDE(TSTOCHECK%COORD,TS(I)%DATA%X, NATOMS, & 45:                        CALL MINPERMDIST(TSTOCHECK%COORD,TS(I)%DATA%X, NATOMS, &
 46:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT) 46:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
 47:                     ELSE 47:                     ELSE
 48:                        CALL NEWMINDIST(TSTOCHECK%COORD,TS(I)%DATA%X,NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT) 48:                        CALL NEWMINDIST(TSTOCHECK%COORD,TS(I)%DATA%X,NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT)
 49:                     ENDIF 49:                     ENDIF
 50:                     !print *, "The distance is",d 50:                     !print *, "The distance is",d
 51:                     IF (D<GEOMDIFFTOL) THEN 51:                     IF (D<GEOMDIFFTOL) THEN
 52:                          !print *, "This TS is already known." 52:                          !print *, "This TS is already known."
 53:                          ISNEWTS=.FALSE. 53:                          ISNEWTS=.FALSE.
 54:                          IF (I > NTSOLD) THEN ! TWO OR MORE TS RETURNED FROM _ONE_ NEB RUN ARE THE SAME 54:                          IF (I > NTSOLD) THEN ! TWO OR MORE TS RETURNED FROM _ONE_ NEB RUN ARE THE SAME
 55:                               PRINT *, 'Shortcut was found! - TS is same as already saved TS #',i  55:                               PRINT *, 'Shortcut was found! - TS is same as already saved TS #',i 
 76:           INTEGER INVERT, INDEX(NATOMS), J2 76:           INTEGER INVERT, INDEX(NATOMS), J2
 77:           LOGICAL PERMUTE,SUCCESS,REDOPATH 77:           LOGICAL PERMUTE,SUCCESS,REDOPATH
 78:           DOUBLE PRECISION D, DIST2, RMAT(3,3) 78:           DOUBLE PRECISION D, DIST2, RMAT(3,3)
 79:  79: 
 80:           MINPOS=NMIN+1 80:           MINPOS=NMIN+1
 81:           NEW=.TRUE. 81:           NEW=.TRUE.
 82:           DO I=1,NMIN 82:           DO I=1,NMIN
 83:              PERMUTE=.FALSE. 83:              PERMUTE=.FALSE.
 84:              IF ( ABS(E-MI(I)%DATA%E) < EDIFFTOL) THEN 84:              IF ( ABS(E-MI(I)%DATA%E) < EDIFFTOL) THEN
 85:                     IF (PERMDIST) THEN 85:                     IF (PERMDIST) THEN
 86:                        CALL ALIGN_DECIDE(COORD,MI(I)%DATA%X(1:3*NATOMS), NATOMS, & 86:                        CALL MINPERMDIST(COORD,MI(I)%DATA%X(1:3*NATOMS), NATOMS, &
 87:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT) 87:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
 88:                     ELSE 88:                     ELSE
 89:                        CALL NEWMINDIST(COORD,MI(I)%DATA%X(1:3*NATOMS),NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT) 89:                        CALL NEWMINDIST(COORD,MI(I)%DATA%X(1:3*NATOMS),NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT)
 90:                     ENDIF 90:                     ENDIF
 91:                     IF (D<GEOMDIFFTOL) THEN 91:                     IF (D<GEOMDIFFTOL) THEN
 92:                        NEW=.FALSE. 92:                        NEW=.FALSE.
 93:                        MINPOS=I 93:                        MINPOS=I
 94:                        RETURN 94:                        RETURN
 95:                     ENDIF 95:                     ENDIF
 96:              ENDIF 96:              ENDIF
177:           IF (NMIN==MINRACKSIZE) CALL REALLOCATEMINRACK177:           IF (NMIN==MINRACKSIZE) CALL REALLOCATEMINRACK
178: 178: 
179:           NMIN=NMIN+1179:           NMIN=NMIN+1
180:           MI(NMIN)%DATA%E => E180:           MI(NMIN)%DATA%E => E
181:           MI(NMIN)%DATA%X => COORD181:           MI(NMIN)%DATA%X => COORD
182:           ALLOCATE( MI(NMIN)%DATA%D(NMIN-1),MI(NMIN)%DATA%NTRIES(NMIN-1),MI(NMIN)%DATA%C(1) )182:           ALLOCATE( MI(NMIN)%DATA%D(NMIN-1),MI(NMIN)%DATA%NTRIES(NMIN-1),MI(NMIN)%DATA%C(1) )
183:           IF (INTERPCOSTFUNCTION) ALLOCATE( MI(NMIN)%DATA%INTERP(NMIN-1))183:           IF (INTERPCOSTFUNCTION) ALLOCATE( MI(NMIN)%DATA%INTERP(NMIN-1))
184: 184: 
185:           DO I=1,NMIN-1185:           DO I=1,NMIN-1
186:                IF (PERMDIST) THEN186:                IF (PERMDIST) THEN
187:                   CALL ALIGN_DECIDE(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:), NATOMS, &187:                   CALL MINPERMDIST(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:), NATOMS, &
188:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)188:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
189:                   D=SQRT(D)189:                   D=SQRT(D)
190:                ELSE190:                ELSE
191:                   CALL NEWMINDIST(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:), NATOMS, D,BULKT,TWOD,'AX   ',.True., &191:                   CALL NEWMINDIST(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:), NATOMS, D,BULKT,TWOD,'AX   ',.True., &
192:   &                               RIGIDBODY,DEBUG,RMAT)192:   &                               RIGIDBODY,DEBUG,RMAT)
193:                ENDIF193:                ENDIF
194:                IF (DEBUG) PRINT '(A,G20.10,2(A,I6))',' addnewmin> distance from MINPERMDIST/NEWMINDIST=',D,' for ',I,' and ',NMIN194:                IF (DEBUG) PRINT '(A,G20.10,2(A,I6))',' addnewmin> distance from MINPERMDIST/NEWMINDIST=',D,' for ',I,' and ',NMIN
195:                MI(NMIN)%DATA%D(I)=D  !  PASSING MIN(NMIN)%DATA%D(I) INSTEAD OF D DOES NOT WORK195:                MI(NMIN)%DATA%D(I)=D  !  PASSING MIN(NMIN)%DATA%D(I) INSTEAD OF D DOES NOT WORK
196:                                      !  perhaps because min is an intrisic function!196:                                      !  perhaps because min is an intrisic function!
197:                IF (INTERPCOSTFUNCTION) MI(NMIN)%DATA%INTERP(I)=INTERPVALUE(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:),D)  197:                IF (INTERPCOSTFUNCTION) MI(NMIN)%DATA%INTERP(I)=INTERPVALUE(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:),D)  
1568:           INTEGER J1, J21568:           INTEGER J1, J2
1569: 1569: 
1570:           IF (ABS(MI(I)%DATA%E-MI(J)%DATA%E) < EDIFFTOL) THEN1570:           IF (ABS(MI(I)%DATA%E-MI(J)%DATA%E) < EDIFFTOL) THEN
1571:              WRITE(*,'(/1x,a,2i5,a)') "Energies of the minima in the pair ",i,j," are the same - checking distance ..."1571:              WRITE(*,'(/1x,a,2i5,a)') "Energies of the minima in the pair ",i,j," are the same - checking distance ..."
1572: 1572: 
1573: !            IF (BULKT) THEN1573: !            IF (BULKT) THEN
1574: !               CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT)1574: !               CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT)
1575: !            ELSEIF (PERMDIST) THEN1575: !            ELSEIF (PERMDIST) THEN
1576: 1576: 
1577:              IF (PERMDIST) THEN1577:              IF (PERMDIST) THEN
1578:                 CALL ALIGN_DECIDE(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS, &1578:                 CALL MINPERMDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS, &
1579:   &                              DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1579:   &                              DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1580:                 D=SQRT(D)1580:                 D=SQRT(D)
1581:              ELSE1581:              ELSE
1582:                 CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT)1582:                 CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,'AX   ',.True.,RIGIDBODY,DEBUG,RMAT)
1583:              ENDIF1583:              ENDIF
1584:              IF (D < GEOMDIFFTOL) THEN1584:              IF (D < GEOMDIFFTOL) THEN
1585:                   WRITE(*,'(1x,a)') "Distance is zero: this should not happen"1585:                   WRITE(*,'(1x,a)') "Distance is zero: this should not happen"
1586:                   CALL TSUMMARY1586:                   CALL TSUMMARY
1587:                   CALL EXIT(10)1587:                   CALL EXIT(10)
1588:              ENDIF1588:              ENDIF
1589:           ELSE1589:           ELSE
1590: !            IF (BULKT) THEN1590: !            IF (BULKT) THEN
1591: !               CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,"AX   ",.False.,RIGIDBODY,DEBUG,RMAT)1591: !               CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,"AX   ",.False.,RIGIDBODY,DEBUG,RMAT)
1592: !            ELSEIF (PERMDIST) THEN1592: !            ELSEIF (PERMDIST) THEN
1593:              IF (PERMDIST) THEN1593:              IF (PERMDIST) THEN
1594:                 CALL ALIGN_DECIDE(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS, &1594:                 CALL MINPERMDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS, &
1595:   &                              DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1595:   &                              DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1596:                 D=SQRT(D)1596:                 D=SQRT(D)
1597:              ELSE1597:              ELSE
1598:                 CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,"AX   ",.False.,RIGIDBODY,DEBUG,RMAT)1598:                 CALL NEWMINDIST(MI(I)%DATA%X,MI(J)%DATA%X,NATOMS,D,BULKT,TWOD,"AX   ",.False.,RIGIDBODY,DEBUG,RMAT)
1599:              ENDIF1599:              ENDIF
1600:           ENDIF1600:           ENDIF
1601: 1601: 
1602:      END SUBROUTINE CHECKPAIR1602:      END SUBROUTINE CHECKPAIR
1603: 1603: 
1604:      FUNCTION GETDISTANCE(I,J)1604:      FUNCTION GETDISTANCE(I,J)


r33372/newconnect.f90 2017-10-04 18:30:14.952274659 +0100 r33371/newconnect.f90 2017-10-04 18:30:19.216330904 +0100
156:              EI=EII;EF=EFF;Q=QQ;FIN=FINFIN;156:              EI=EII;EF=EFF;Q=QQ;FIN=FINFIN;
157:              PRINT '(3(A,I6))',' newconnect> redopoints file contains ',(NREDO-1)/2,' transition states and ',(NREDO+1)/2,' minima'157:              PRINT '(3(A,I6))',' newconnect> redopoints file contains ',(NREDO-1)/2,' transition states and ',(NREDO+1)/2,' minima'
158:              IF (MOD(NREDO,2).EQ.0) THEN158:              IF (MOD(NREDO,2).EQ.0) THEN
159:                 PRINT '(A,I6,A)',' newconnect> WARNING *** # of coordinate sets read from redopoints is ', &159:                 PRINT '(A,I6,A)',' newconnect> WARNING *** # of coordinate sets read from redopoints is ', &
160:   &                          NREDO,' which is not an odd number!'160:   &                          NREDO,' which is not an odd number!'
161:              ENDIF161:              ENDIF
162:  162:  
163:              PRINT '(A,F20.10)',' newconnect> start and finish minima overwritten from redopoints file'163:              PRINT '(A,F20.10)',' newconnect> start and finish minima overwritten from redopoints file'
164:              WRITE(*,'(a,2(g20.10,a))') ' newconnect> Initial energy=',EI,' RMS force=',RMSINITIAL164:              WRITE(*,'(a,2(g20.10,a))') ' newconnect> Initial energy=',EI,' RMS force=',RMSINITIAL
165:              WRITE(*,'(a,2(g20.10,a))') ' newconnect> Final   energy=',EF,' RMS force=',RMSFINAL165:              WRITE(*,'(a,2(g20.10,a))') ' newconnect> Final   energy=',EF,' RMS force=',RMSFINAL
166:              CALL ALIGN_DECIDE(QQ,FINFIN,NA,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ENDPOINTSEP,DIST2,RIGIDBODY,RMAT)166:              CALL MINPERMDIST(QQ,FINFIN,NA,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,ENDPOINTSEP,DIST2,RIGIDBODY,RMAT)
167:              WRITE(*,'(A,G20.10)') ' newconnect> separation=',ENDPOINTSEP167:              WRITE(*,'(A,G20.10)') ' newconnect> separation=',ENDPOINTSEP
168:           ENDIF168:           ENDIF
169: 169: 
170:           CALL INITIALISE(NA,EI,Q,EF,FIN,ENDPOINTSEP)170:           CALL INITIALISE(NA,EI,Q,EF,FIN,ENDPOINTSEP)
171: 171: 
172:           IF (REDOPATHXYZ) PRINT '(A)',' newconnect> Redo run will use available path.<n>.xyz files'172:           IF (REDOPATHXYZ) PRINT '(A)',' newconnect> Redo run will use available path.<n>.xyz files'
173:           IF (REDOPATH.AND.(.NOT.REDOPATHXYZ)) THEN173:           IF (REDOPATH.AND.(.NOT.REDOPATHXYZ)) THEN
174:              REWIND(99)174:              REWIND(99)
175:              IF (AMHT) THEN175:              IF (AMHT) THEN
176:                 NCOUNT=0176:                 NCOUNT=0


r33372/newneb.f90 2017-10-04 18:30:15.616283420 +0100 r33371/newneb.f90 2017-10-04 18:30:19.912340088 +0100
374: !                     IF (.NOT.WHOLEDNEB) THEN374: !                     IF (.NOT.WHOLEDNEB) THEN
375:                       IF (.FALSE.) THEN375:                       IF (.FALSE.) THEN
376:                          IF (ALLOCATED(INTNEBIMAGES)) DEALLOCATE(INTNEBIMAGES)376:                          IF (ALLOCATED(INTNEBIMAGES)) DEALLOCATE(INTNEBIMAGES)
377:                          ALLOCATE(INTNEBIMAGES(NIMAGE*NOPT))377:                          ALLOCATE(INTNEBIMAGES(NIMAGE*NOPT))
378:                          PERMDISTSAVE=PERMDIST378:                          PERMDISTSAVE=PERMDIST
379:                          PERMDIST=.FALSE.379:                          PERMDIST=.FALSE.
380:                          LDEBUG = DEBUG380:                          LDEBUG = DEBUG
381:                          DEBUG = .FALSE. 381:                          DEBUG = .FALSE. 
382:                          CALL GENRIGID_IMAGE_RIGIDTOC(NIMAGE, XYZ)382:                          CALL GENRIGID_IMAGE_RIGIDTOC(NIMAGE, XYZ)
383:                          DEBUG = LDEBUG383:                          DEBUG = LDEBUG
384:                    WRITE(*,*) "sn402: Changed NEWNEB so that it calls ALIGN_DECIDE instead of MINPERMDIST"384:                          CALL MINPERMDIST(XYZ(NOPT+1:2*NOPT),XYZ(1:3*NATOMS),NATOMS,DEBUG, &
385:                    WRITE(*,*) "At the time of writing, this WRITE statement is inside an IF(FALSE) block, so it hasn't been tested" 
386:                    WRITE(*,*) "If you are reading this, please check carefully that this part of the code is working as you expect," 
387:                    WRITE(*,*) "then remove this message!" 
388:                          CALL ALIGN_DECIDE(XYZ(NOPT+1:2*NOPT),XYZ(1:3*NATOMS),NATOMS,DEBUG, & 
389:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)385:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
390:                          CALL ALIGN_DECIDE(XYZ(NOPT*NIMAGE+1:NOPT*(NIMAGE+1)), &386:                          CALL MINPERMDIST(XYZ(NOPT*NIMAGE+1:NOPT*(NIMAGE+1)), &
391:   &                       XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2)),NATOMS,DEBUG, &387:   &                       XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2)),NATOMS,DEBUG, &
392:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)388:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
393:                          PERMDIST=PERMDISTSAVE389:                          PERMDIST=PERMDISTSAVE
394: 390: 
395:                          TOTALDIST=0.0D0391:                          TOTALDIST=0.0D0
396:                          DO J2=1,NIMAGE+1392:                          DO J2=1,NIMAGE+1
397:                             XDUMMY=0.0D0393:                             XDUMMY=0.0D0
398:                             DO J5=1,3*NATOMS394:                             DO J5=1,3*NATOMS
399:                                XDUMMY=XDUMMY+( XYZ((J2-1)*3*NATOMS+J5) - XYZ(J2*3*NATOMS+J5) )**2395:                                XDUMMY=XDUMMY+( XYZ((J2-1)*3*NATOMS+J5) - XYZ(J2*3*NATOMS+J5) )**2
400:                             ENDDO396:                             ENDDO


r33372/nnutils.f90 2017-10-04 18:30:15.840286372 +0100 r33371/nnutils.f90 2017-10-04 18:30:20.144343149 +0100
509:            QS(1:3*NATOMS)=XYZ(1:NOPT)509:            QS(1:3*NATOMS)=XYZ(1:NOPT)
510:            QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)510:            QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)
511: !511: !
512: !  Now change the whole XYZ array to be based on the permutational isomer that512: !  Now change the whole XYZ array to be based on the permutational isomer that
513: !  minimises the overall distance. The difference is in the initial straight line513: !  minimises the overall distance. The difference is in the initial straight line
514: !  guess, which may be for different permutational isomers.514: !  guess, which may be for different permutational isomers.
515: !  Use the permutation that gives the lowest value for the highest image.515: !  Use the permutation that gives the lowest value for the highest image.
516: !516: !
517:            QS2(1:3*NATOMS)=XYZ(1:NOPT)517:            QS2(1:3*NATOMS)=XYZ(1:NOPT)
518:            QF2(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)518:            QF2(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)
519:            CALL ALIGN_DECIDE(QS2,QF2,NATOMS, &519:            CALL MINPERMDIST(QS2,QF2,NATOMS, &
520:   &                         DEBUG,BULK_BOXVEC(1),BULK_BOXVEC(2),BULK_BOXVEC(3),BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)520:   &                         DEBUG,BULK_BOXVEC(1),BULK_BOXVEC(2),BULK_BOXVEC(3),BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
521:            EWORST2=-HUGE(1.0D0)521:            EWORST2=-HUGE(1.0D0)
522:            DO J1=1,NIMAGE522:            DO J1=1,NIMAGE
523:               COORDS(1:3*NATOMS)=(J1*QF2(1:3*NATOMS)+(NIMAGE+1-J1)*QS2(1:3*NATOMS))/(NIMAGE+1)523:               COORDS(1:3*NATOMS)=(J1*QF2(1:3*NATOMS)+(NIMAGE+1-J1)*QS2(1:3*NATOMS))/(NIMAGE+1)
524:               CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)524:               CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)
525:               LEIMAGE2(J1)=E1525:               LEIMAGE2(J1)=E1
526:               IF (E1.GT.EWORST2) THEN526:               IF (E1.GT.EWORST2) THEN
527:                  JWORST2=J1527:                  JWORST2=J1
528:                  EWORST2=E1528:                  EWORST2=E1
529:               ENDIF529:               ENDIF
1369: !                   CALL TSUMMARY1369: !                   CALL TSUMMARY
1370: !                   STOP1370: !                   STOP
1371: !              ENDIF1371: !              ENDIF
1372:                DO J2=1,NIMAGE+21372:                DO J2=1,NIMAGE+2
1373: !1373: !
1374: ! Here we skip two lines, allowing for the second line to be blank.1374: ! Here we skip two lines, allowing for the second line to be blank.
1375: !1375: !
1376:                   READ(993,*)1376:                   READ(993,*)
1377:                   READ(993,*)1377:                   READ(993,*)
1378:                   READ(993,*) (DUMMYS,XYZ( (J2-1)*NOPT+J1),XYZ((J2-1)*NOPT+J1+1),XYZ((J2-1)*NOPT+J1+2),J1=1,NOPT,3)1378:                   READ(993,*) (DUMMYS,XYZ( (J2-1)*NOPT+J1),XYZ((J2-1)*NOPT+J1+1),XYZ((J2-1)*NOPT+J1+2),J1=1,NOPT,3)
1379:                   IF (PERMGUESS.AND.(J2.GT.1)) THEN1379:                   IF (PERMGUESS.AND.(J2.GT.1)) THEN 
1380:                    WRITE(*,*) "sn402: Changed NNUTILS::RWG so that it calls ALIGN_DECIDE instead of MINPERMDIST for PERMGUESS"1380:                      CALL MINPERMDIST(XYZ((J2-2)*NOPT+1:(J2-1)*NOPT),XYZ((J2-1)*NOPT+1:J2*NOPT),NATOMS,DEBUG, &
1381:                    WRITE(*,*) "I haven't tested this change and am not certain whether it''s sensible."  
1382:                    WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove this message!"  
1383:                      CALL ALIGN_DECIDE(XYZ((J2-2)*NOPT+1:(J2-1)*NOPT),XYZ((J2-1)*NOPT+1:J2*NOPT),NATOMS,DEBUG, & 
1384:   &                    BULK_BOXVEC(1),BULK_BOXVEC(2),BULK_BOXVEC(3),BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1381:   &                    BULK_BOXVEC(1),BULK_BOXVEC(2),BULK_BOXVEC(3),BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1385:                   ENDIF1382:                   ENDIF
1386:                ENDDO1383:                ENDDO
1387:                PRINT '(A,I6,A)','nnutils> ',NIMAGE+2,' images read from file ' // TRIM(ADJUSTL(GUESSFILE))1384:                PRINT '(A,I6,A)','nnutils> ',NIMAGE+2,' images read from file ' // TRIM(ADJUSTL(GUESSFILE))
1388:           END SELECT1385:           END SELECT
1389:           CLOSE(UNIT=993)1386:           CLOSE(UNIT=993)
1390:      END SUBROUTINE RWG1387:      END SUBROUTINE RWG
1391: 1388: 
1392:      SUBROUTINE SAVEBANDCOORD1389:      SUBROUTINE SAVEBANDCOORD
1393:           USE NEBDATA,ONLY: XYZ,NOPT1390:           USE NEBDATA,ONLY: XYZ,NOPT


r33372/OPTIM.F 2017-10-04 18:30:16.060289276 +0100 r33371/OPTIM.F 2017-10-04 18:30:20.368346103 +0100
 85:       DOUBLE PRECISION :: XRIGIDCOORDS(3*NATOMS), XRIGIDGRAD(3*NATOMS) 85:       DOUBLE PRECISION :: XRIGIDCOORDS(3*NATOMS), XRIGIDGRAD(3*NATOMS)
 86:       INTEGER GETUNIT 86:       INTEGER GETUNIT
 87:       COMMON /OEPATH/ ETS,EPLUS,EMINUS 87:       COMMON /OEPATH/ ETS,EPLUS,EMINUS
 88: ! Print OPTIM version in the output 88: ! Print OPTIM version in the output
 89:       !VERSIONTEMP=25661 89:       !VERSIONTEMP=25661
 90:       !WRITE(*, '(A,I5)') ' OPTIM> version r',VERSIONTEMP 90:       !WRITE(*, '(A,I5)') ' OPTIM> version r',VERSIONTEMP
 91: C 91: C
 92: C  Dynamic memory allocation 92: C  Dynamic memory allocation
 93: C 93: C
 94:       ALLOCATE (FROZEN(NATOMS),ZSYM(NATOMS),NR(NATOMS)) 94:       ALLOCATE (FROZEN(NATOMS),ZSYM(NATOMS),NR(NATOMS))
 95:       IF (DEBUG) WRITE(*,*) ' OPTIM> allocated FROZEN and ZSYM with dimension NATOMS=',NATOMS 95:       IF (DEBUG) WRITE(*,*) ' OPTIM> allocated ZSYM with dimension NATOMS=',NATOMS
 96:       ALLOCATE (FROZENRES(NATOMS)) 96:       ALLOCATE (FROZENRES(NATOMS))
 97: C      STPMAX(:)=0.0D0 97: C      STPMAX(:)=0.0D0
 98:  98: 
 99:       FILTH=F1 ; FILTH2=F2 ; FILTHSTR=TRIM(ADJUSTL(FLSTRING)) 99:       FILTH=F1 ; FILTH2=F2 ; FILTHSTR=TRIM(ADJUSTL(FLSTRING))
100:       KNOWE=.FALSE. ; KNOWG=.FALSE. ; KNOWH=.FALSE.100:       KNOWE=.FALSE. ; KNOWG=.FALSE. ; KNOWH=.FALSE.
101:       RBATOMSMAX=10101:       RBATOMSMAX=10
102: 102: 
103:       CALL KEYWORDS(Q)103:       CALL KEYWORDS(Q)
104:       ALLOCATE (STPMAX(NOPT))104:       ALLOCATE (STPMAX(NOPT))
105:       IF (UNRST.AND.(RKMIN.OR.BSMIN.OR.(INR.GT.-1))) THEN105:       IF (UNRST.AND.(RKMIN.OR.BSMIN.OR.(INR.GT.-1))) THEN
677:             IF (AMBER12T.AND.NOCISTRANS) THEN ! check the finish geometry for cis/trans wrt start677:             IF (AMBER12T.AND.NOCISTRANS) THEN ! check the finish geometry for cis/trans wrt start
678:                CALL CIS_TRANS_CHECK(FIN,GOODENDPOINT)678:                CALL CIS_TRANS_CHECK(FIN,GOODENDPOINT)
679:                IF (.NOT.GOODENDPOINT) THEN679:                IF (.NOT.GOODENDPOINT) THEN
680:                   WRITE(*,'(A)') 680:                   WRITE(*,'(A)') 
681:      &     ' tryconnect> Cis-trans isomerisation of a peptide bond detected in finish end point (wrt. the original structure), quit' 681:      &     ' tryconnect> Cis-trans isomerisation of a peptide bond detected in finish end point (wrt. the original structure), quit' 
682:                   STOP682:                   STOP
683:                ELSE683:                ELSE
684:                   WRITE(*,'(A)') ' tryconnect> finish endpoint cis/trans consistent with starting geometry'684:                   WRITE(*,'(A)') ' tryconnect> finish endpoint cis/trans consistent with starting geometry'
685:                ENDIF685:                ENDIF
686:             ENDIF686:             ENDIF
687:             CALL ALIGN_DECIDE(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)687:             CALL MINPERMDIST(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
688:             IF (PERMDISTINIT) PERMDIST=.FALSE.688:             IF (PERMDISTINIT) PERMDIST=.FALSE.
689:             IF (ATOMMATCHINIT) ATOMMATCHDIST=.FALSE.689:             IF (ATOMMATCHINIT) ATOMMATCHDIST=.FALSE.
690:             IF (BISECTT) THEN690:             IF (BISECTT) THEN
691:                CALL BISECT_OPT(NATOMS,EINITIAL,Q,EFINAL,FIN,DIST)691:                CALL BISECT_OPT(NATOMS,EINITIAL,Q,EFINAL,FIN,DIST)
692:             ELSE692:             ELSE
693:                IF (ALLOCATED(SAVES)) DEALLOCATE(SAVES)693:                IF (ALLOCATED(SAVES)) DEALLOCATE(SAVES)
694:                IF (ALLOCATED(SAVEF)) DEALLOCATE(SAVEF)694:                IF (ALLOCATED(SAVEF)) DEALLOCATE(SAVEF)
695:                ALLOCATE(SAVES(NOPT),SAVEF(NOPT))695:                ALLOCATE(SAVES(NOPT),SAVEF(NOPT))
696:                SAVES(1:NOPT)=Q(1:NOPT)696:                SAVES(1:NOPT)=Q(1:NOPT)
697:                SAVEF(1:NOPT)=FIN(1:NOPT)697:                SAVEF(1:NOPT)=FIN(1:NOPT)
744:                DIST=DIST+(Q(J1)-FIN(J1))**2744:                DIST=DIST+(Q(J1)-FIN(J1))**2
745:             ENDDO745:             ENDDO
746:             DIST=SQRT(DIST)746:             DIST=SQRT(DIST)
747:             NNNIMAGE=NINT(MIN(MECIMDENS*DIST,MECMAXIMAGES*1.0D0)) ! IMAGE DENSITY TIMES DISTANCE747:             NNNIMAGE=NINT(MIN(MECIMDENS*DIST,MECMAXIMAGES*1.0D0)) ! IMAGE DENSITY TIMES DISTANCE
748:             IF (NNNIMAGE < 1       ) NNNIMAGE=1748:             IF (NNNIMAGE < 1       ) NNNIMAGE=1
749:             NITERMAX=NINT(MIN(NNNIMAGE*MECITDENS,MECMAXIT*1.0D0)) ! NUMBER OF IMAGES TIMES ITERATION DENSITY749:             NITERMAX=NINT(MIN(NNNIMAGE*MECITDENS,MECMAXIT*1.0D0)) ! NUMBER OF IMAGES TIMES ITERATION DENSITY
750:             CALL MECCANO(.FALSE.,.TRUE.,ENERGY,.FALSE.,Q,FIN,E1,E2,RMSINITIAL,RMSFINAL)750:             CALL MECCANO(.FALSE.,.TRUE.,ENERGY,.FALSE.,Q,FIN,E1,E2,RMSINITIAL,RMSFINAL)
751:          ENDIF751:          ENDIF
752:       ELSE IF (NEWNEBT.OR.NEBT.OR.GROWSTRINGT) THEN752:       ELSE IF (NEWNEBT.OR.NEBT.OR.GROWSTRINGT) THEN
753:          IF (NEBMIND) THEN753:          IF (NEBMIND) THEN
754:             CALL ALIGN_DECIDE(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)754:             CALL MINPERMDIST(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
755:             REALSTR=WR(DIST,3)755:             REALSTR=WR(DIST,3)
756:             WRITE(*,'(a,f12.2)')' OPTIM> Structures were put in the closest coincidence, distance = '//trim(RealStr)756:             WRITE(*,'(a,f12.2)')' OPTIM> Structures were put in the closest coincidence, distance = '//trim(RealStr)
757:          ENDIF757:          ENDIF
758:          IF (UNRST) THEN758:          IF (UNRST) THEN
759:              DO J1=1,NRES759:              DO J1=1,NRES
760:                 C(1,J1)=Q(6*(J1-1)+1)760:                 C(1,J1)=Q(6*(J1-1)+1)
761:                 C(2,J1)=Q(6*(J1-1)+2)761:                 C(2,J1)=Q(6*(J1-1)+2)
762:                 C(3,J1)=Q(6*(J1-1)+3)762:                 C(3,J1)=Q(6*(J1-1)+3)
763:                 C(1,J1+NRES)=Q(6*(J1-1)+4)763:                 C(1,J1+NRES)=Q(6*(J1-1)+4)
764:                 C(2,J1+NRES)=Q(6*(J1-1)+5)764:                 C(2,J1+NRES)=Q(6*(J1-1)+5)
887:             ENDIF887:             ENDIF
888: ! to get the instanton rates.888: ! to get the instanton rates.
889:             CALL GEOPT(FNAMEF,EFNAME,Q,VECS,MFLAG,ENERGY,EVALMIN,VNEW) 889:             CALL GEOPT(FNAMEF,EFNAME,Q,VECS,MFLAG,ENERGY,EVALMIN,VNEW) 
890:          ENDIF890:          ENDIF
891: 891: 
892: 892: 
893: 893: 
894: C     ELSE IF (DRAGT) THEN894: C     ELSE IF (DRAGT) THEN
895: C        CALL DRAG895: C        CALL DRAG
896:       ELSE IF (CLOSESTALIGNMENT) THEN896:       ELSE IF (CLOSESTALIGNMENT) THEN
897:          CALL ALIGN_DECIDE(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)897:          CALL MINPERMDIST(Q,FIN,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
898:          WRITE(*,'(a,F10.5)') ' OPTIM> Minimized distance between aligned structures = ',DIST898:          WRITE(*,'(a,F10.5)') ' OPTIM> Minimized distance between aligned structures = ',DIST
899:          IF (CHRMMT) THEN 899:          IF (CHRMMT) THEN 
900:             CALL CHARMMDUMP(FIN,'aligned.crd',MACHINE)900:             CALL CHARMMDUMP(FIN,'aligned.crd',MACHINE)
901:          ELSE IF (AMBERT .OR. NABT .OR. AMBER12T) THEN901:          ELSE IF (AMBERT .OR. NABT .OR. AMBER12T) THEN
902:             ! Formats used come from the AMBER routine minrit902:             ! Formats used come from the AMBER routine minrit
903:             OPEN(UNIT=3,FILE='aligned.rst', STATUS='UNKNOWN')903:             OPEN(UNIT=3,FILE='aligned.rst', STATUS='UNKNOWN')
904:             WRITE(3,'(20A4)') 'MOL'904:             WRITE(3,'(20A4)') 'MOL'
905:             IF (NATOMS.GT.100000) THEN905:             IF (NATOMS.GT.100000) THEN
906:                WRITE(3,'(I6)') NATOMS906:                WRITE(3,'(I6)') NATOMS
907:             ELSE907:             ELSE


r33372/RPH.f90 2017-10-04 18:30:16.324292757 +0100 r33371/RPH.f90 2017-10-04 18:30:20.592349058 +0100
136:    PRINT '(A,G20.10)','RPH> for first minimum Q=',QORDER136:    PRINT '(A,G20.10)','RPH> for first minimum Q=',QORDER
137: ENDIF137: ENDIF
138: 138: 
139: DO J1=2,NPATH-1139: DO J1=2,NPATH-1
140:    IF (EOFS(J1-1)+EDIFFTOL*10.0D0<EOFS(J1).AND.EOFS(J1)>EOFS(J1+1)+EDIFFTOL*10.0D0) THEN140:    IF (EOFS(J1-1)+EDIFFTOL*10.0D0<EOFS(J1).AND.EOFS(J1)>EOFS(J1+1)+EDIFFTOL*10.0D0) THEN
141:       NSTRUCTREF=NSTRUCTREF+1141:       NSTRUCTREF=NSTRUCTREF+1
142:       REFSTRUCTFRAME(NSTRUCTREF)=J1142:       REFSTRUCTFRAME(NSTRUCTREF)=J1
143:       VREF(NSTRUCTREF)=EOFS(J1) 143:       VREF(NSTRUCTREF)=EOFS(J1) 
144:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF144:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF
145:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)145:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)
146:       WRITE(*,*) "sn402: Changed RPH so that it calls ALIGN_DECIDE instead of MINPERMDIST"146:       CALL MINPERMDIST(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &
147:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
148:       WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove these messages!" 
149:       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, & 
150:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)147:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
151:       PRINT '(A,I10,2(A,G20.10))','RPH> transition state identified for frame ',J1,' energy=', &148:       PRINT '(A,I10,2(A,G20.10))','RPH> transition state identified for frame ',J1,' energy=', &
152:   &         VREF(NSTRUCTREF),' distance to previous minimum=',DUMMY149:   &         VREF(NSTRUCTREF),' distance to previous minimum=',DUMMY
153:       NTS=NTS+1150:       NTS=NTS+1
154:       IF (.NOT.LASTMIN) THEN151:       IF (.NOT.LASTMIN) THEN
155:          PRINT '(A)','RPH> WARNING *** previous stationary point identified was not a minimum'152:          PRINT '(A)','RPH> WARNING *** previous stationary point identified was not a minimum'
156:          STOP153:          STOP
157:       ENDIF154:       ENDIF
158:       LASTMIN=.FALSE.155:       LASTMIN=.FALSE.
159:       LASTTS=.TRUE.156:       LASTTS=.TRUE.
166:          CALL GETORDER(QTEMP,QORDER,NATOMS)163:          CALL GETORDER(QTEMP,QORDER,NATOMS)
167:          PRINT '(A,G20.10)','RPH> for this transition state Q=',QORDER164:          PRINT '(A,G20.10)','RPH> for this transition state Q=',QORDER
168:       ENDIF165:       ENDIF
169:    ENDIF166:    ENDIF
170:    IF ((EOFS(J1-1)>EOFS(J1).AND.EOFS(J1)<=EOFS(J1+1)).AND.(ABS(EOFS(J1)-VREF(NSTRUCTREF)).GT.EDIFFTOL)) THEN167:    IF ((EOFS(J1-1)>EOFS(J1).AND.EOFS(J1)<=EOFS(J1+1)).AND.(ABS(EOFS(J1)-VREF(NSTRUCTREF)).GT.EDIFFTOL)) THEN
171:       NSTRUCTREF=NSTRUCTREF+1168:       NSTRUCTREF=NSTRUCTREF+1
172:       REFSTRUCTFRAME(NSTRUCTREF)=J1169:       REFSTRUCTFRAME(NSTRUCTREF)=J1
173:       VREF(NSTRUCTREF)=EOFS(J1) 170:       VREF(NSTRUCTREF)=EOFS(J1) 
174:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF171:       IF (VREF(NSTRUCTREF).GT.VREF(HIGHESTREF)) HIGHESTREF=NSTRUCTREF
175:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)172:       COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(J1,1:3*NATOMS)
176:       WRITE(*,*) "sn402: Changed RPH so that it calls ALIGN_DECIDE instead of MINPERMDIST"173:       CALL MINPERMDIST(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &
177:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
178:       WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove these messages!" 
179:       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, & 
180:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)174:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
181:       PRINT '(A,I10,2(A,G20.10))','RPH> minimum identified for frame ',J1,' energy=', &175:       PRINT '(A,I10,2(A,G20.10))','RPH> minimum identified for frame ',J1,' energy=', &
182:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY176:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY
183:       NMIN=NMIN+1177:       NMIN=NMIN+1
184:       IF (.NOT.LASTTS) THEN178:       IF (.NOT.LASTTS) THEN
185:          PRINT '(A)','RPH> WARNING *** previous stationary point identified was not a transition state'179:          PRINT '(A)','RPH> WARNING *** previous stationary point identified was not a transition state'
186:       ENDIF180:       ENDIF
187:       LASTTS=.FALSE.181:       LASTTS=.FALSE.
188:       LASTMIN=.TRUE.182:       LASTMIN=.TRUE.
189:       IF (.TRUE.) THEN183:       IF (.TRUE.) THEN
201:    PERMDISTSAVE=.FALSE. 195:    PERMDISTSAVE=.FALSE. 
202:    LPERMDISTSAVE=.FALSE.196:    LPERMDISTSAVE=.FALSE.
203:    PERMDIST=.FALSE. 197:    PERMDIST=.FALSE. 
204:    LPERMDIST=.FALSE.198:    LPERMDIST=.FALSE.
205: ENDIF199: ENDIF
206: IF (ABS(EOFS(NPATH)-VREF(NSTRUCTREF)).GT.EDIFFTOL) THEN200: IF (ABS(EOFS(NPATH)-VREF(NSTRUCTREF)).GT.EDIFFTOL) THEN
207:    NSTRUCTREF=NSTRUCTREF+1201:    NSTRUCTREF=NSTRUCTREF+1
208:    REFSTRUCTFRAME(NSTRUCTREF)=NPATH202:    REFSTRUCTFRAME(NSTRUCTREF)=NPATH
209:    VREF(NSTRUCTREF)=EOFS(NPATH) ! last minimum203:    VREF(NSTRUCTREF)=EOFS(NPATH) ! last minimum
210:    COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(NPATH,1:3*NATOMS)204:    COORDSREF(NSTRUCTREF,1:3*NATOMS)=PATHFRAMES(NPATH,1:3*NATOMS)
211:       WRITE(*,*) "sn402: Changed RPH so that it calls ALIGN_DECIDE instead of MINPERMDIST"205:       CALL MINPERMDIST(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, &
212:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
213:       WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove these messages!" 
214:       CALL ALIGN_DECIDE(COORDSREF(NSTRUCTREF-1,1:3*NATOMS),COORDSREF(NSTRUCTREF,1:3*NATOMS),NATOMS,DEBUG, & 
215:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)206:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
216:       PRINT '(A,I10,2(A,G20.10))','RPH> minimum assumed    for frame ',J1,' energy=', &207:       PRINT '(A,I10,2(A,G20.10))','RPH> minimum assumed    for frame ',J1,' energy=', &
217:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY208:   &         VREF(NSTRUCTREF),' distance to previous transition state=',DUMMY
218:    NMIN=NMIN+1209:    NMIN=NMIN+1
219:    IF (.NOT.LASTTS) THEN210:    IF (.NOT.LASTTS) THEN
220:       PRINT '(A)','RPH> WARNING *** previous stationary point identified was not a transition state'211:       PRINT '(A)','RPH> WARNING *** previous stationary point identified was not a transition state'
221:    ENDIF212:    ENDIF
222:    IF (.TRUE.) THEN213:    IF (.TRUE.) THEN
223:       DO K1=1,3214:       DO K1=1,3
224:          DO K2=1,NATOMS215:          DO K2=1,NATOMS
228:       CALL GETORDER(QTEMP,QORDER,NATOMS)219:       CALL GETORDER(QTEMP,QORDER,NATOMS)
229:       PRINT '(A,G20.10)','RPH> for this minimum          Q=',QORDER220:       PRINT '(A,G20.10)','RPH> for this minimum          Q=',QORDER
230:    ENDIF221:    ENDIF
231: ENDIF222: ENDIF
232: PRINT '(4(A,I10))','RPH> Total number of reference structures=',NSTRUCTREF,' with ',NTS,' ts and ',NMIN,' minima'223: PRINT '(4(A,I10))','RPH> Total number of reference structures=',NSTRUCTREF,' with ',NTS,' ts and ',NMIN,' minima'
233: CALL FLUSH(6)224: CALL FLUSH(6)
234: !225: !
235: ! work out total distance and slice width - use this for SLICELABEL226: ! work out total distance and slice width - use this for SLICELABEL
236: !227: !
237: DISTTOTAL=0.0D0228: DISTTOTAL=0.0D0
238:       WRITE(*,*) "sn402: Changed RPH so that it calls ALIGN_DECIDE instead of MINPERMDIST"229: CALL MINPERMDIST(COORDSREF(1,1:3*NATOMS),PATHFRAMES(1,1:3*NATOMS),NATOMS,DEBUG, &
239:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
240:       WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove these messages!" 
241: CALL ALIGN_DECIDE(COORDSREF(1,1:3*NATOMS),PATHFRAMES(1,1:3*NATOMS),NATOMS,DEBUG, & 
242:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)230:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
243: PRINT '(A,F20.10)','RPH> first distance=',DUMMY231: PRINT '(A,F20.10)','RPH> first distance=',DUMMY
244: DISTTOTAL=DUMMY232: DISTTOTAL=DUMMY
245: DO J1=2,NPATH233: DO J1=2,NPATH
246:       WRITE(*,*) "sn402: Changed RPH so that it calls ALIGN_DECIDE instead of MINPERMDIST"234:    CALL MINPERMDIST(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, &
247:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
248:       WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove these messages!" 
249:    CALL ALIGN_DECIDE(PATHFRAMES(J1-1,1:3*NATOMS),PATHFRAMES(J1,1:3*NATOMS),NATOMS,DEBUG, & 
250:   &                 PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)235:   &                 PARAM1,PARAM2,PARAM3,BULKT,TWOD,DUMMY,DIST2,RIGIDBODY,RMAT)
251:    PRINT '(A,I10,F20.10)','RPH> frame and distance=',J1,DUMMY236:    PRINT '(A,I10,F20.10)','RPH> frame and distance=',J1,DUMMY
252:    DISTTOTAL=DISTTOTAL+DUMMY237:    DISTTOTAL=DISTTOTAL+DUMMY
253: ENDDO238: ENDDO
254: !239: !
255: ! Need to resave COORDSREF to get consistent permutational alignment.240: ! Need to resave COORDSREF to get consistent permutational alignment.
256: !241: !
257: DO J1=1,NSTRUCTREF 242: DO J1=1,NSTRUCTREF 
258:    COORDSREF(J1,1:3*NATOMS)=PATHFRAMES(REFSTRUCTFRAME(J1),1:3*NATOMS)243:    COORDSREF(J1,1:3*NATOMS)=PATHFRAMES(REFSTRUCTFRAME(J1),1:3*NATOMS)
259: ENDDO244: ENDDO
362: NEXMODES=NEXMODES+1347: NEXMODES=NEXMODES+1
363: ALLOCATE(VSAVE(NOPT,NEXMODES))348: ALLOCATE(VSAVE(NOPT,NEXMODES))
364: !349: !
365: ! Now step through intervening slice geometries with additional projection along the gradient350: ! Now step through intervening slice geometries with additional projection along the gradient
366: ! (or the frame displacement vector?).351: ! (or the frame displacement vector?).
367: !352: !
368: DISTTOTAL2=0.0D0353: DISTTOTAL2=0.0D0
369: ATFRAME=2354: ATFRAME=2
370: DO J1=1,RPHSLICES355: DO J1=1,RPHSLICES
371:    DO WHILE (DISTTOTAL2.LT.SLICELABEL(J1)) 356:    DO WHILE (DISTTOTAL2.LT.SLICELABEL(J1)) 
372:       WRITE(*,*) "sn402: Changed RPH so that it calls ALIGN_DECIDE instead of MINPERMDIST"357:       CALL MINPERMDIST(PATHFRAMES(ATFRAME-1,1:3*NATOMS),PATHFRAMES(ATFRAME,1:3*NATOMS),NATOMS,DEBUG, &
373:       WRITE(*,*) "I haven't tested this change and am not certain whether it's sensible."  
374:       WRITE(*,*) "Please check carefully that this part of the code is working as you expect, then remove these messages!" 
375:       CALL ALIGN_DECIDE(PATHFRAMES(ATFRAME-1,1:3*NATOMS),PATHFRAMES(ATFRAME,1:3*NATOMS),NATOMS,DEBUG, & 
376:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSAVE,DIST2,RIGIDBODY,RMAT)358:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DSAVE,DIST2,RIGIDBODY,RMAT)
377: 359: 
378:       DISTTOTAL2=DISTTOTAL2+DSAVE360:       DISTTOTAL2=DISTTOTAL2+DSAVE
379:       PRINT '(A,2I10,2F20.10)','RPH> slice, frame, distance, and total=',J1,ATFRAME,DSAVE,DISTTOTAL2361:       PRINT '(A,2I10,2F20.10)','RPH> slice, frame, distance, and total=',J1,ATFRAME,DSAVE,DISTTOTAL2
380:       ATFRAME=ATFRAME+1362:       ATFRAME=ATFRAME+1
381:       IF (ATFRAME.GT.NPATH) PRINT '(A)','WARNING *** Reached last frame but still in slices loop'363:       IF (ATFRAME.GT.NPATH) PRINT '(A)','WARNING *** Reached last frame but still in slices loop'
382:    ENDDO364:    ENDDO
383:    FRACTION=(SLICELABEL(J1)-DISTTOTAL2+DSAVE)/DSAVE365:    FRACTION=(SLICELABEL(J1)-DISTTOTAL2+DSAVE)/DSAVE
384:    IF ((FRACTION.LT.0.0D0).OR.(FRACTION.GT.1.0)) THEN366:    IF ((FRACTION.LT.0.0D0).OR.(FRACTION.GT.1.0)) THEN
385:       PRINT '(A,I6,A,I6,A,G20.10,A,2F20.10)','RPH> Fractional distance for frames ',ATFRAME,' and ',ATFRAME-1,' is ',FRACTION367:       PRINT '(A,I6,A,I6,A,G20.10,A,2F20.10)','RPH> Fractional distance for frames ',ATFRAME,' and ',ATFRAME-1,' is ',FRACTION


r33372/tryconnect.f90 2017-10-04 18:30:15.392280464 +0100 r33371/tryconnect.f90 2017-10-04 18:30:19.692337187 +0100
103: 103: 
104:           IF (CHRMMT) NCHENCALLS = 999 ! update non-bonded list on next call to potential.104:           IF (CHRMMT) NCHENCALLS = 999 ! update non-bonded list on next call to potential.
105: !105: !
106: !  Subroutine CHECKPAIR puts the endpoints into optimal alignment.106: !  Subroutine CHECKPAIR puts the endpoints into optimal alignment.
107: !107: !
108:           IF (.NOT.REDOPATH) THEN108:           IF (.NOT.REDOPATH) THEN
109:              CALL CHECKPAIR(JS,JF,PERMTEST)109:              CALL CHECKPAIR(JS,JF,PERMTEST)
110:           ELSE110:           ELSE
111:             PERMTEST=.FALSE.111:             PERMTEST=.FALSE.
112:             IF (ABS(MI(JS)%DATA%E-MI(JF)%DATA%E) < EDIFFTOL) PERMTEST=.TRUE. ! must initialise this logical112:             IF (ABS(MI(JS)%DATA%E-MI(JF)%DATA%E) < EDIFFTOL) PERMTEST=.TRUE. ! must initialise this logical
113:             CALL ALIGN_DECIDE(TSREDO,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)113:             CALL MINPERMDIST(TSREDO,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
114:             D1INIT=D114:             D1INIT=D
115:             CALL ALIGN_DECIDE(TSREDO,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)115:             CALL MINPERMDIST(TSREDO,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
116:             D2INIT=D116:             D2INIT=D
117: 117: 
118:              IF (DEBUG) PRINT '(A,2F20.10)',' tryconnect> Initial distances of transition state to minima are :',D1INIT,D2INIT118:              IF (DEBUG) PRINT '(A,2F20.10)',' tryconnect> Initial distances of transition state to minima are :',D1INIT,D2INIT
119:           ENDIF119:           ENDIF
120: 120: 
121:           IF (GUESSPATHT) THEN121:           IF (GUESSPATHT) THEN
122:             IF (UNRST) THEN122:             IF (UNRST) THEN
123:                DO J1=1,NRES123:                DO J1=1,NRES
124:                   C(1,J1)=MI(JS)%DATA%X(6*(J1-1)+1)124:                   C(1,J1)=MI(JS)%DATA%X(6*(J1-1)+1)
125:                   C(2,J1)=MI(JS)%DATA%X(6*(J1-1)+2)125:                   C(2,J1)=MI(JS)%DATA%X(6*(J1-1)+2)
151:             IF (NIMAGE < 2       ) PRINT*,'WARNING - Nimage is < 2'151:             IF (NIMAGE < 2       ) PRINT*,'WARNING - Nimage is < 2'
152:             NIMAGE=NINTERP152:             NIMAGE=NINTERP
153: !           NIterMax = Nimage*IterDensity ! try zero neb iterations if we have a GUESSPATH path153: !           NIterMax = Nimage*IterDensity ! try zero neb iterations if we have a GUESSPATH path
154:             NITERMAX = 0154:             NITERMAX = 0
155:             IF (NINTERP.LT.2) THEN ! NO IMAGES FROM GUESSPATH - REVERT TO USUAL SCHEME155:             IF (NINTERP.LT.2) THEN ! NO IMAGES FROM GUESSPATH - REVERT TO USUAL SCHEME
156:                IF (.NOT.(NCONDONE==1 .AND. FCD)) THEN156:                IF (.NOT.(NCONDONE==1 .AND. FCD)) THEN
157: !                   NIMAGE=IMAGEDENSITY*MI(JF)%DATA%D(JS) &157: !                   NIMAGE=IMAGEDENSITY*MI(JF)%DATA%D(JS) &
158: !                       +IMAGEINCR*IMAGEDENSITY*MI(JF)%DATA%D(JS)*(MI(JF)%DATA%NTRIES(JS)-1)158: !                       +IMAGEINCR*IMAGEDENSITY*MI(JF)%DATA%D(JS)*(MI(JF)%DATA%NTRIES(JS)-1)
159:                   DISTFAC=MI(JF)%DATA%D(JS)159:                   DISTFAC=MI(JF)%DATA%D(JS)
160:                   IF (DIJKSTRALOCAL.NE.1.0D0) THEN160:                   IF (DIJKSTRALOCAL.NE.1.0D0) THEN
161:                      CALL ALIGN_DECIDE(MI(JS)%DATA%X,MI(JF)%DATA%X,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD, &161:                      CALL MINPERMDIST(MI(JS)%DATA%X,MI(JF)%DATA%X,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD, &
162:      &                                   DISTFAC,DIST2,RIGIDBODY,RMAT)162:      &                                   DISTFAC,DIST2,RIGIDBODY,RMAT)
163:                   ENDIF163:                   ENDIF
164:                   IF (.NOT.(USEINT.OR.USEINTLJ)) THEN164:                   IF (.NOT.(USEINT.OR.USEINTLJ)) THEN
165:                      NIMAGE=DISTFAC*(IMAGEDENSITY+IMAGEINCR*(MI(JF)%DATA%NTRIES(JS)-1))165:                      NIMAGE=DISTFAC*(IMAGEDENSITY+IMAGEINCR*(MI(JF)%DATA%NTRIES(JS)-1))
166:                      IF (NIMAGE >= IMAGEMAX) THEN166:                      IF (NIMAGE >= IMAGEMAX) THEN
167:                         NIMAGE = IMAGEMAX167:                         NIMAGE = IMAGEMAX
168:                         MI(JF)%DATA%NTRIES(JS)=NTRIESMAX ! no point trying again with the same number of images168:                         MI(JF)%DATA%NTRIES(JS)=NTRIESMAX ! no point trying again with the same number of images
169:                      ENDIF169:                      ENDIF
170:                      IF (NIMAGE < 2       ) NIMAGE = 2170:                      IF (NIMAGE < 2       ) NIMAGE = 2
171:                      NITERMAX = NIMAGE*ITERDENSITY171:                      NITERMAX = NIMAGE*ITERDENSITY
1039: 1039: 
1040: 333            CONTINUE1040: 333            CONTINUE
1041:                CALL ISNEWMIN(EPLUS,QPLUS,MINPLUSPOS,PLUSNEW,REDOPATH,PERMUTE,INVERT,INDEX,IMATCH)1041:                CALL ISNEWMIN(EPLUS,QPLUS,MINPLUSPOS,PLUSNEW,REDOPATH,PERMUTE,INVERT,INDEX,IMATCH)
1042:                CALL ISNEWMIN(EMINUS,QMINUS,MINMINUSPOS,MINUSNEW,REDOPATH,PERMUTE,INVERT,INDEX,IMATCH)1042:                CALL ISNEWMIN(EMINUS,QMINUS,MINMINUSPOS,MINUSNEW,REDOPATH,PERMUTE,INVERT,INDEX,IMATCH)
1043: !1043: !
1044: ! The above check will not discover the case when the plus minimum is new, and is the same1044: ! The above check will not discover the case when the plus minimum is new, and is the same
1045: ! as the minus minimum.1045: ! as the minus minimum.
1046: !1046: !
1047:                IF (PLUSNEW.AND.MINUSNEW) THEN1047:                IF (PLUSNEW.AND.MINUSNEW) THEN
1048:                   IF (ABS(EMINUS-EPLUS) < EDIFFTOL) THEN1048:                   IF (ABS(EMINUS-EPLUS) < EDIFFTOL) THEN
1049:                      CALL ALIGN_DECIDE(QMINUS,QPLUS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1049:                      CALL MINPERMDIST(QMINUS,QPLUS,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1050:                      IF (D<GEOMDIFFTOL) THEN1050:                      IF (D<GEOMDIFFTOL) THEN
1051:                         MINUSNEW=.FALSE.1051:                         MINUSNEW=.FALSE.
1052:                         MINMINUSPOS=MINPLUSPOS1052:                         MINMINUSPOS=MINPLUSPOS
1053:                      ENDIF1053:                      ENDIF
1054:                   ENDIF1054:                   ENDIF
1055:                ENDIF1055:                ENDIF
1056: 1056: 
1057:                EDUMMY=TS(I)%DATA%E1057:                EDUMMY=TS(I)%DATA%E
1058:                TMPTS(1:NOPT)=TS(I)%DATA%X(1:NOPT)1058:                TMPTS(1:NOPT)=TS(I)%DATA%X(1:NOPT)
1059:                IF (DUMPALLPATHS) CALL MAKEALLPATHINFO(TMPTS,QPLUS,QMINUS,EDUMMY,EPLUS,EMINUS,FRQSTS,FRQSPLUS,FRQSMINUS)1059:                IF (DUMPALLPATHS) CALL MAKEALLPATHINFO(TMPTS,QPLUS,QMINUS,EDUMMY,EPLUS,EMINUS,FRQSTS,FRQSPLUS,FRQSMINUS)
1096:                   PUSHOFF=SAVEPUSHOFF1096:                   PUSHOFF=SAVEPUSHOFF
1097:                ENDIF1097:                ENDIF
1098:           ENDDO1098:           ENDDO
1099:         ENDIF  1099:         ENDIF  
1100: 1100: 
1101:           !IF(ALLOCATED(FLATPATHT)) DEALLOCATE(FLATPATHT)1101:           !IF(ALLOCATED(FLATPATHT)) DEALLOCATE(FLATPATHT)
1102: !1102: !
1103: ! Allow for new pathway calculation with different PUSHOFF and MAXBFGS1103: ! Allow for new pathway calculation with different PUSHOFF and MAXBFGS
1104: !1104: !
1105:           IF (REDOPATH.AND.(.NOT.PATHFAILT)) THEN1105:           IF (REDOPATH.AND.(.NOT.PATHFAILT)) THEN
1106:              CALL ALIGN_DECIDE(QP,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1106:              CALL MINPERMDIST(QP,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1107:              DIST1P=D1107:              DIST1P=D
1108: 1108: 
1109:              CALL ALIGN_DECIDE(QM,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1109:              CALL MINPERMDIST(QM,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1110:              DIST1M=D1110:              DIST1M=D
1111: 1111: 
1112:              CALL ALIGN_DECIDE(QM,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1112:              CALL MINPERMDIST(QM,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1113:              DIST2M=D1113:              DIST2M=D
1114: 1114: 
1115:              CALL ALIGN_DECIDE(QP,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1115:              CALL MINPERMDIST(QP,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1116:              DIST2P=D1116:              DIST2P=D
1117: 1117: 
1118:              PATHFAILED=.FALSE.1118:              PATHFAILED=.FALSE.
1119:              IF ((DIST1P.GT.GEOMDIFFTOL).AND.(DIST1M.GT.GEOMDIFFTOL)) THEN1119:              IF ((DIST1P.GT.GEOMDIFFTOL).AND.(DIST1M.GT.GEOMDIFFTOL)) THEN
1120:                 PRINT '(A)',' tryconnect> path failed to match first minimum'1120:                 PRINT '(A)',' tryconnect> path failed to match first minimum'
1121:                 PATHFAILED=.TRUE.1121:                 PATHFAILED=.TRUE.
1122:              ENDIF1122:              ENDIF
1123:              IF ((DIST2P.GT.GEOMDIFFTOL).AND.(DIST2M.GT.GEOMDIFFTOL)) THEN1123:              IF ((DIST2P.GT.GEOMDIFFTOL).AND.(DIST2M.GT.GEOMDIFFTOL)) THEN
1124:                 PRINT '(A)',' tryconnect> path failed to match second minimum'1124:                 PRINT '(A)',' tryconnect> path failed to match second minimum'
1125:                 PATHFAILED=.TRUE.1125:                 PATHFAILED=.TRUE.


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0