hdiff output

r30873/connect_unconnected.f90 2016-07-28 07:30:11.213517269 +0100 r30872/connect_unconnected.f90 2016-07-28 07:30:12.817538887 +0100
 55: SUBROUTINE CONNECTUNCONNECTED(NAVAIL) 55: SUBROUTINE CONNECTUNCONNECTED(NAVAIL)
 56: USE PORFUNCS 56: USE PORFUNCS
 57: USE COMMONS 57: USE COMMONS
 58: !USE COMMONS, ONLY: UMIN, NMIN, NMINA, NMINB, NCPU, DMIN1, DMIN2, CONNECTLOWESTT, NATOMS, & 58: !USE COMMONS, ONLY: UMIN, NMIN, NMINA, NMINB, NCPU, DMIN1, DMIN2, CONNECTLOWESTT, NATOMS, &
 59: !&                  CONNECTETHRESHT, CONNECTDTHRESHT, REFMIN,  DEBUG, BOXLX, BOXLY, BOXLZ, & 59: !&                  CONNECTETHRESHT, CONNECTDTHRESHT, REFMIN,  DEBUG, BOXLX, BOXLY, BOXLZ, &
 60: !  &                BULKT, TWOD, RIGIDBODY, EMIN, NTS, LOCATIONA, LOCATIONB, NCONNMAX, NATT 60: !  &                BULKT, TWOD, RIGIDBODY, EMIN, NTS, LOCATIONA, LOCATIONB, NCONNMAX, NATT
 61: IMPLICIT NONE 61: IMPLICIT NONE
 62:  62: 
 63: INTEGER :: J1, J2, J3, J4, J5, J6, J7, J8, NPUSED, CLOSESTMIN 63: INTEGER :: J1, J2, J3, J4, J5, J6, J7, J8, NPUSED, CLOSESTMIN
 64: LOGICAL :: ISA(NMIN), ISB(NMIN), CHANGED, DEADTS(NTS), UNCONN(NMIN), DISTCALC 64: LOGICAL :: ISA(NMIN), ISB(NMIN), CHANGED, DEADTS(NTS), UNCONN(NMIN), DISTCALC
 65: INTEGER NOFF(NMIN) 65: INTEGER :: NCOL(NMIN), NVAL(NCONNMAX,NMIN), NDISTA(NMIN), NDISTB(NMIN), NCYCLE, DMIN
 66: INTEGER :: NCOL(NMIN), NVAL(2*NTS), NDISTA(NMIN), NDISTB(NMIN), NCYCLE, DMIN 
 67: INTEGER :: DMAX, NUNCONA, NUNCONB, NDEAD, NUNCONN, PAIRSTODO, NAVAIL 66: INTEGER :: DMAX, NUNCONA, NUNCONB, NDEAD, NUNCONN, PAIRSTODO, NAVAIL
 68: INTEGER, ALLOCATABLE :: UNCMIN(:),CONMIN(:) 67: INTEGER, ALLOCATABLE :: UNCMIN(:),CONMIN(:)
 69: DOUBLE PRECISION, ALLOCATABLE :: EUNCON(:), UNCONDIST(:), EUNDIFF(:), DISTCON(:) 68: DOUBLE PRECISION, ALLOCATABLE :: EUNCON(:), UNCONDIST(:), EUNDIFF(:), DISTCON(:)
 70: DOUBLE PRECISION :: DMATMC(2*NTS), KSUM(NMIN), CUT_UNDERFLOW, LOWESTX(3*NATOMS),SETX(3*NATOMS) 69: DOUBLE PRECISION :: DMATMC(NCONNMAX,NMIN), KSUM(NMIN), CUT_UNDERFLOW, LOWESTX(3*NATOMS),SETX(3*NATOMS)
 71: DOUBLE PRECISION :: DISTMIN, EREF, DISTANCE, DIST2, RMAT(3,3), REFX(3*NATOMS) 70: DOUBLE PRECISION :: DISTMIN, EREF, DISTANCE, DIST2, RMAT(3,3), REFX(3*NATOMS)
 72:  71: 
 73: ! 72: !
 74: ! Set up connectivity analysis 73: ! Set up connectivity analysis
 75: ! 74: !
 76: UNCONN(1:NMIN)=.FALSE. 75: UNCONN(1:NMIN)=.FALSE.
 77: ISA(1:NMIN)=.FALSE. 76: ISA(1:NMIN)=.FALSE.
 78: ISB(1:NMIN)=.FALSE. 77: ISB(1:NMIN)=.FALSE.
 79: DO J1=1,NMINA 78: DO J1=1,NMINA
 80:    ISA(LOCATIONA(J1))=.TRUE. 79:    ISA(LOCATIONA(J1))=.TRUE.
 81: ENDDO 80: ENDDO
 82: DO J1=1,NMINB 81: DO J1=1,NMINB
 83:    ISB(LOCATIONB(J1))=.TRUE. 82:    ISB(LOCATIONB(J1))=.TRUE.
 84: ENDDO 83: ENDDO
 85:  84: 
 86:  85: 
 87: CUT_UNDERFLOW=-300.0D0 86: CUT_UNDERFLOW=-300.0D0
 88: CALL RATECONST_SETUP(KSUM,DEADTS,NDEAD,.FALSE.,CUT_UNDERFLOW) 87: CALL RATECONST_SETUP(KSUM,DEADTS,NDEAD,.FALSE.,CUT_UNDERFLOW)
 89:  88: 
 90:  89: 
 91: CALL MAKEDLIN(DMATMC,NCOL,NOFF,NVAL,DEADTS,.TRUE.,ISA,ISB,KSUM) 90: CALL MAKED(DMATMC,NCOL,NVAL,DEADTS,.TRUE.,ISA,ISB,KSUM)
 92:  91: 
 93:  92: 
 94: ! 93: !
 95: ! Check connections to A set 94: ! Check connections to A set
 96: ! 95: !
 97: NDISTA(1:NMIN)=1000000 96: NDISTA(1:NMIN)=1000000
 98: DO J1=1,NMINA 97: DO J1=1,NMINA
 99:    NDISTA(LOCATIONA(J1))=0 98:    NDISTA(LOCATIONA(J1))=0
100: ENDDO  99: ENDDO 
101: NCYCLE=0100: NCYCLE=0
102: 5    CHANGED=.FALSE.101: 5    CHANGED=.FALSE.
103: NCYCLE=NCYCLE+1102: NCYCLE=NCYCLE+1
104: DMIN=100000103: DMIN=100000
105: DMAX=0104: DMAX=0
106: NUNCONA=0105: NUNCONA=0
107: DO J1=1,NMIN106: DO J1=1,NMIN
108:    IF (NDISTA(J1).EQ.0) CYCLE ! A minimum107:    IF (NDISTA(J1).EQ.0) CYCLE ! A minimum
109:    DO J2=1,NCOL(J1)108:    DO J2=1,NCOL(J1)
110:       IF (NDISTA(NVAL(NOFF(J1)+J2))+1.LT.NDISTA(J1)) THEN109:       IF (NDISTA(NVAL(J2,J1))+1.LT.NDISTA(J1)) THEN
111:          CHANGED=.TRUE.110:          CHANGED=.TRUE.
112:          NDISTA(J1)=NDISTA(NVAL(NOFF(J1)+J2))+1111:          NDISTA(J1)=NDISTA(NVAL(J2,J1))+1
113:       ENDIF112:       ENDIF
114:    ENDDO113:    ENDDO
115:    IF ((NDISTA(J1).GT.DMAX).AND.(NDISTA(J1).NE.1000000)) DMAX=NDISTA(J1)114:    IF ((NDISTA(J1).GT.DMAX).AND.(NDISTA(J1).NE.1000000)) DMAX=NDISTA(J1)
116:    IF (NDISTA(J1).LT.DMIN) DMIN=NDISTA(J1)115:    IF (NDISTA(J1).LT.DMIN) DMIN=NDISTA(J1)
117:    IF (NDISTA(J1).EQ.1000000) NUNCONA=NUNCONA+1116:    IF (NDISTA(J1).EQ.1000000) NUNCONA=NUNCONA+1
118: ENDDO117: ENDDO
119: IF (CHANGED) GOTO 5118: IF (CHANGED) GOTO 5
120: PRINT '(3(A,I8))','Dijkstra> steps to A region converged in ',NCYCLE-1, &119: PRINT '(3(A,I8))','Dijkstra> steps to A region converged in ',NCYCLE-1, &
121: &                    ' cycles; maximum=',DMAX,' disconnected=',NUNCONA120: &                    ' cycles; maximum=',DMAX,' disconnected=',NUNCONA
122: !121: !
128: ENDDO127: ENDDO
129: NCYCLE=0128: NCYCLE=0
130: 51    CHANGED=.FALSE.129: 51    CHANGED=.FALSE.
131: NCYCLE=NCYCLE+1130: NCYCLE=NCYCLE+1
132: DMIN=100000131: DMIN=100000
133: DMAX=0132: DMAX=0
134: NUNCONB=0133: NUNCONB=0
135: DO J1=1,NMIN134: DO J1=1,NMIN
136:    IF (NDISTB(J1).EQ.0) CYCLE ! B minimum135:    IF (NDISTB(J1).EQ.0) CYCLE ! B minimum
137:    DO J2=1,NCOL(J1)136:    DO J2=1,NCOL(J1)
138:       IF (NDISTB(NVAL(NOFF(J1)+J2))+1.LT.NDISTB(J1)) THEN137:       IF (NDISTB(NVAL(J2,J1))+1.LT.NDISTB(J1)) THEN
139:          CHANGED=.TRUE.138:          CHANGED=.TRUE.
140:          NDISTB(J1)=NDISTB(NVAL(NOFF(J1)+J2))+1139:          NDISTB(J1)=NDISTB(NVAL(J2,J1))+1
141:       ENDIF140:       ENDIF
142:    ENDDO141:    ENDDO
143:    IF ((NDISTB(J1).GT.DMAX).AND.(NDISTB(J1).NE.1000000)) DMAX=NDISTB(J1)142:    IF ((NDISTB(J1).GT.DMAX).AND.(NDISTB(J1).NE.1000000)) DMAX=NDISTB(J1)
144:    IF (NDISTB(J1).LT.DMIN) DMIN=NDISTB(J1)143:    IF (NDISTB(J1).LT.DMIN) DMIN=NDISTB(J1)
145:    IF (NDISTB(J1).EQ.1000000) NUNCONB=NUNCONB+1144:    IF (NDISTB(J1).EQ.1000000) NUNCONB=NUNCONB+1
146: ENDDO 145: ENDDO 
147: IF (CHANGED) GOTO 51146: IF (CHANGED) GOTO 51
148: PRINT '(3(A,I8))','Dijkstra> steps to B region converged in ',NCYCLE-1, &147: PRINT '(3(A,I8))','Dijkstra> steps to B region converged in ',NCYCLE-1, &
149: &                    ' cycles; maximum=',DMAX,' disconnected=',NUNCONB148: &                    ' cycles; maximum=',DMAX,' disconnected=',NUNCONB
150: !149: !


r30873/getallpaths.f 2016-07-28 07:30:11.501521150 +0100 r30872/getallpaths.f 2016-07-28 07:30:13.117542931 +0100
297:             BADTRIPLE=.TRUE.297:             BADTRIPLE=.TRUE.
298:             GOTO 120298:             GOTO 120
299:          ENDIF299:          ENDIF
300:          DO J2=1,NTS300:          DO J2=1,NTS
301:             DISTANCE=1.0D100301:             DISTANCE=1.0D100
302:             IF (ABS(NEWETS-ETS(J2)).LT.EDIFFTOL) THEN302:             IF (ABS(NEWETS-ETS(J2)).LT.EDIFFTOL) THEN
303:                READ(UTS,REC=J2) (LOCALPOINTS2(J3),J3=1,NVARS)303:                READ(UTS,REC=J2) (LOCALPOINTS2(J3),J3=1,NVARS)
304:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 304:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 
305:      &                          RMAT,.FALSE.)305:      &                          RMAT,.FALSE.)
306:             ENDIF306:             ENDIF
307: !           PRINT '(A,I8,5G20.10)','J2,NEWETS,ETS(J2),diff,tol,DISTANCE=',J2,NEWETS,ETS(J2),ABS(NEWETS-ETS(J2)),EDIFFTOL,DISTANCE307: !           PRINT '(A,I8,3G20.10)','J2,NEWETS,ETS(J2),DISTANCE=',J2,NEWETS,ETS(J2),DISTANCE
308:             IF ((ABS(NEWETS-ETS(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN308:             IF ((ABS(NEWETS-ETS(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN
309:                IF (DEBUG) PRINT '(2(A,I6))','getallpaths> path ts ',J1,' is database ts ',J2309:                IF (DEBUG) PRINT '(2(A,I6))','getallpaths> path ts ',J1,' is database ts ',J2
310:                IF (ABS(NEWFVIBTS-FVIBTS(J2))/FVIBTS(J2).GT.1.0D-3) THEN310:                IF (ABS(NEWFVIBTS-FVIBTS(J2))/FVIBTS(J2).GT.1.0D-3) THEN
311:                   WRITE(*,'(A,F15.5,A,F15.5)') 'getallpaths> WARNING, NEWFVIBTS=',NEWFVIBTS,' should be ',FVIBTS(J2)311:                   WRITE(*,'(A,F15.5,A,F15.5)') 'getallpaths> WARNING, NEWFVIBTS=',NEWFVIBTS,' should be ',FVIBTS(J2)
312:                ENDIF312:                ENDIF
313:                IF (NEWHORDERTS.NE.HORDERTS(J2)) THEN313:                IF (NEWHORDERTS.NE.HORDERTS(J2)) THEN
314:                   WRITE(*,'(A,I6,A,I6)') 'getallpaths> ERROR, NEWHORDERTS=',NEWHORDERTS,' should be ',HORDERTS(J2)314:                   WRITE(*,'(A,I6,A,I6)') 'getallpaths> ERROR, NEWHORDERTS=',NEWHORDERTS,' should be ',HORDERTS(J2)
315:                   NEWHORDERTS=MAX(NEWHORDERTS,HORDERTS(J2))315:                   NEWHORDERTS=MAX(NEWHORDERTS,HORDERTS(J2))
316:                   WRITE(*,'(A,I6)') 'getallpaths> using maximum value: ',NEWHORDERTS316:                   WRITE(*,'(A,I6)') 'getallpaths> using maximum value: ',NEWHORDERTS
317:                ENDIF317:                ENDIF
477:                         IF (NTS.GT.MAXTS) CALL TSDOUBLE477:                         IF (NTS.GT.MAXTS) CALL TSDOUBLE
478:                         NEWTS=NTS478:                         NEWTS=NTS
479:                         ETS(NTS)=NEWETS479:                         ETS(NTS)=NEWETS
480:                         FVIBTS(NTS)=NEWFVIBTS480:                         FVIBTS(NTS)=NEWFVIBTS
481:                         HORDERTS(NTS)=NEWHORDERTS481:                         HORDERTS(NTS)=NEWHORDERTS
482:                         IXTS(NTS)=NEWIXTS482:                         IXTS(NTS)=NEWIXTS
483:                         IYTS(NTS)=NEWIYTS483:                         IYTS(NTS)=NEWIYTS
484:                         IZTS(NTS)=NEWIZTS484:                         IZTS(NTS)=NEWIZTS
485:                         NEGEIG(NTS)=NEWNEGEIG485:                         NEGEIG(NTS)=NEWNEGEIG
486:                         IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0486:                         IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0
487: !                       IF (DEBUG) THEN487:                         IF (DEBUG) THEN
488:                            WRITE(*,'(A,I6,A)') 'getallpaths> new intermediate ts recounted for different connections ',NTS488:                            WRITE(*,'(A,I6,A)') 'getallpaths> new intermediate ts recounted for different connections ',NTS
489: !                       ENDIF489:                         ENDIF
490:                         PLUS(NTS)=PLUSMIN490:                         PLUS(NTS)=PLUSMIN
491:                         WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NVARS)491:                         WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NVARS)
492:                         CALL FLUSH(UTS)492:                         CALL FLUSH(UTS)
493:                         GOTO 140493:                         GOTO 140
494:                      ELSE494:                      ELSE
495:                         FAILED=.TRUE.495:                         FAILED=.TRUE.
496: !                       IF (DEBUG) WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections'496:                         IF (DEBUG) WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections'
497:                         WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections' 
498:                      ENDIF497:                      ENDIF
499:                   ENDIF498:                   ENDIF
500:                ELSEIF (ABS(EPLUS-EMIN(MINUS(TSNUMBER))).LT.EDIFFTOL) THEN499:                ELSEIF (ABS(EPLUS-EMIN(MINUS(TSNUMBER))).LT.EDIFFTOL) THEN
501:                   IF (ABS(NEWEMIN-EMIN(PLUS(TSNUMBER))).LT.EDIFFTOL) THEN500:                   IF (ABS(NEWEMIN-EMIN(PLUS(TSNUMBER))).LT.EDIFFTOL) THEN
502:                      CALL INERTIAWRAPPER(NEWPOINTSMIN,NOPT,angleAxis,IXPLUS,IYPLUS,IZPLUS)501:                      CALL INERTIAWRAPPER(NEWPOINTSMIN,NOPT,angleAxis,IXPLUS,IYPLUS,IZPLUS)
503:                      CALL INERTIAWRAPPER(NEWPOINTSMINPLUS,NOPT,angleAxis,IXMINUS,IYMINUS,IZMINUS)502:                      CALL INERTIAWRAPPER(NEWPOINTSMINPLUS,NOPT,angleAxis,IXMINUS,IYMINUS,IZMINUS)
504:                      IF (DEBUG) THEN503:                      IF (DEBUG) THEN
505:                         PRINT '(A)','getallpaths> assigning + to - and - to +, energies are:'504:                         PRINT '(A)','getallpaths> assigning + to - and - to +, energies are:'
506:                         PRINT '(A,2G25.15)','getallpaths> new: ',EPLUS,NEWEMIN505:                         PRINT '(A,2G25.15)','getallpaths> new: ',EPLUS,NEWEMIN
507:                         PRINT '(A,2G25.15)','getallpaths> old: ',EMIN(MINUS(TSNUMBER)),EMIN(PLUS(TSNUMBER))506:                         PRINT '(A,2G25.15)','getallpaths> old: ',EMIN(MINUS(TSNUMBER)),EMIN(PLUS(TSNUMBER))
518:                         IF (NTS.GT.MAXTS) CALL TSDOUBLE517:                         IF (NTS.GT.MAXTS) CALL TSDOUBLE
519:                         NEWTS=NTS518:                         NEWTS=NTS
520:                         ETS(NTS)=NEWETS519:                         ETS(NTS)=NEWETS
521:                         FVIBTS(NTS)=NEWFVIBTS520:                         FVIBTS(NTS)=NEWFVIBTS
522:                         HORDERTS(NTS)=NEWHORDERTS521:                         HORDERTS(NTS)=NEWHORDERTS
523:                         IXTS(NTS)=NEWIXTS522:                         IXTS(NTS)=NEWIXTS
524:                         IYTS(NTS)=NEWIYTS523:                         IYTS(NTS)=NEWIYTS
525:                         IZTS(NTS)=NEWIZTS524:                         IZTS(NTS)=NEWIZTS
526:                         NEGEIG(NTS)=NEWNEGEIG525:                         NEGEIG(NTS)=NEWNEGEIG
527:                         IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0526:                         IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0
528: !                       IF (DEBUG) THEN527:                         IF (DEBUG) THEN
529:                            WRITE(*,'(A,I6,A)') 'getallpaths> new intermediate ts recounted for different connections ',NTS528:                            WRITE(*,'(A,I6,A)') 'getallpaths> new intermediate ts recounted for different connections ',NTS
530: !                       ENDIF529:                         ENDIF
531:                         PLUS(NTS)=PLUSMIN530:                         PLUS(NTS)=PLUSMIN
532:                         WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NVARS)531:                         WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NVARS)
533:                         CALL FLUSH(UTS)532:                         CALL FLUSH(UTS)
534:                         GOTO 140533:                         GOTO 140
535:                      ELSE534:                      ELSE
536:                         FAILED=.TRUE.535:                         FAILED=.TRUE.
537: !                       IF (DEBUG) WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections'536:                         IF (DEBUG) WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections'
538:                         WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections' 
539:                      ENDIF537:                      ENDIF
540:                   ENDIF538:                   ENDIF
541:                ELSE539:                ELSE
542:                   PRINT '(A,I6,A,G20.10)','getallpaths> WARNING failed consistency check for old ts ',540:                   PRINT '(A,I6,A,G20.10)','getallpaths> WARNING failed consistency check for old ts ',
543:      &                                     TSNUMBER,' energy=',ETS(TSNUMBER)541:      &                                     TSNUMBER,' energy=',ETS(TSNUMBER)
544:                   PRINT '(A,2G25.15)','getallpaths> +/- energies=',EPLUS,NEWEMIN542:                   PRINT '(A,2G25.15)','getallpaths> +/- energies=',EPLUS,NEWEMIN
545:                   PRINT '(A,I8,A,2G25.15)','getallpaths> ts ',TSNUMBER,543:                   PRINT '(A,I8,A,2G25.15)','getallpaths> ts ',TSNUMBER,
546:      &                                     ' minima are: ',EMIN(PLUS(TSNUMBER)),EMIN(MINUS(TSNUMBER))544:      &                                     ' minima are: ',EMIN(PLUS(TSNUMBER)),EMIN(MINUS(TSNUMBER))
547:                   IF (ALLTST) THEN545:                   IF (ALLTST) THEN
548:                      TSISOLD=.FALSE.546:                      TSISOLD=.FALSE.
550:                      IF (NTS.GT.MAXTS) CALL TSDOUBLE548:                      IF (NTS.GT.MAXTS) CALL TSDOUBLE
551:                      NEWTS=NTS549:                      NEWTS=NTS
552:                      ETS(NTS)=NEWETS550:                      ETS(NTS)=NEWETS
553:                      FVIBTS(NTS)=NEWFVIBTS551:                      FVIBTS(NTS)=NEWFVIBTS
554:                      HORDERTS(NTS)=NEWHORDERTS552:                      HORDERTS(NTS)=NEWHORDERTS
555:                      IXTS(NTS)=NEWIXTS553:                      IXTS(NTS)=NEWIXTS
556:                      IYTS(NTS)=NEWIYTS554:                      IYTS(NTS)=NEWIYTS
557:                      IZTS(NTS)=NEWIZTS555:                      IZTS(NTS)=NEWIZTS
558:                      NEGEIG(NTS)=NEWNEGEIG556:                      NEGEIG(NTS)=NEWNEGEIG
559:                      IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0557:                      IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0
560: !                    IF (DEBUG) THEN558:                      IF (DEBUG) THEN
561:                         WRITE(*,'(A,I6,A)') 'getallpaths> new intermediate ts recounted for different connections ',NTS559:                         WRITE(*,'(A,I6,A)') 'getallpaths> new intermediate ts recounted for different connections ',NTS
562: !                    ENDIF560:                      ENDIF
563:                      PLUS(NTS)=PLUSMIN561:                      PLUS(NTS)=PLUSMIN
564:                      WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NVARS)562:                      WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NVARS)
565:                      CALL FLUSH(UTS)563:                      CALL FLUSH(UTS)
566:                      GOTO 140564:                      GOTO 140
567:                   ELSE565:                   ELSE
568:                      FAILED=.TRUE.566:                      FAILED=.TRUE.
569: !                    IF (DEBUG) WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections'567:                      IF (DEBUG) WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections'
570:                      WRITE(*,'(A)') 'getallpaths> not recounting this ts for different connections' 
571:                   ENDIF568:                   ENDIF
572:                ENDIF569:                ENDIF
573:                IF ((.NOT.FAILED).AND.(.NOT.BULKT)) THEN570:                IF ((.NOT.FAILED).AND.(.NOT.BULKT)) THEN
574:                   IF (ABS(IXPLUS-IXMIN(PLUS(TSNUMBER))).GT.IDIFFTOL) THEN571:                   IF (ABS(IXPLUS-IXMIN(PLUS(TSNUMBER))).GT.IDIFFTOL) THEN
575:                      PRINT '(A,2G20.10)','getallpaths> WARNING failed consistency check for inertia +x: ',572:                      PRINT '(A,2G20.10)','getallpaths> WARNING failed consistency check for inertia +x: ',
576:      &                                    IXPLUS,IXMIN(PLUS(TSNUMBER))573:      &                                    IXPLUS,IXMIN(PLUS(TSNUMBER))
577:                   ENDIF574:                   ENDIF
578:                   IF (ABS(IYPLUS-IYMIN(PLUS(TSNUMBER))).GT.IDIFFTOL) THEN575:                   IF (ABS(IYPLUS-IYMIN(PLUS(TSNUMBER))).GT.IDIFFTOL) THEN
579:                      PRINT '(A,2G20.10)','getallpaths> WARNING failed consistency check for inertia +y: ',576:                      PRINT '(A,2G20.10)','getallpaths> WARNING failed consistency check for inertia +y: ',
580:      &                                    IYPLUS,IYMIN(PLUS(TSNUMBER))577:      &                                    IYPLUS,IYMIN(PLUS(TSNUMBER))


r30873/KMCcommit.f 2016-07-28 07:30:10.925513387 +0100 r30872/KMCcommit.f 2016-07-28 07:30:12.545535221 +0100
799:                DMATTMP(1:NCOL(J1))=DMAT0(1:NCOL(J1),J1)799:                DMATTMP(1:NCOL(J1))=DMAT0(1:NCOL(J1),J1)
800:                CALL SORT4(NCOL(J1),NCONNMAX,DMATTMP,NVALTMP)800:                CALL SORT4(NCOL(J1),NCONNMAX,DMATTMP,NVALTMP)
801:                NVAL(1:NCOL(J1),J1)=NVALTMP(1:NCOL(J1))801:                NVAL(1:NCOL(J1),J1)=NVALTMP(1:NCOL(J1))
802:                DMAT0(1:NCOL(J1),J1)=DMATTMP(1:NCOL(J1))802:                DMAT0(1:NCOL(J1),J1)=DMATTMP(1:NCOL(J1))
803:             ENDIF803:             ENDIF
804:          ENDDO804:          ENDDO
805:       ENDIF805:       ENDIF
806: 806: 
807:       RETURN807:       RETURN
808:       END808:       END
809: ! 
810: !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
811: ! 
812: !  Construct DMATMC using an economical storage format 
813: ! 
814: !  NCOL(M2)     = # connections for minimum M2 
815: !  NVAL(J1,M2)  = index of minimum involved in connection J1 from minimum M2 
816: !  DMAT0(J1,M2) = KMC-type probability of taking connection J1 from minimum  
817: !                 M2 to minimum NVAL(J1,M2) 
818: ! 
819: !  Degenerate rearrangements are excluded. 
820: ! 
821:       SUBROUTINE MAKEDLIN(DMAT0,NCOL,NOFF,NVAL,DEADTS,SORTT,ISA,ISB,KSUM) 
822:       USE COMMONS 
823:       IMPLICIT NONE 
824:       INTEGER NCOL(NMIN), J1, M1, M2, J2, NVAL(2*NTS), NVALTMP(NCONNMAX), OTHER, NOFF(NMIN) 
825:       DOUBLE PRECISION DMATTMP(NCONNMAX), DUMMY, KSUM(NMIN) 
826:       LOGICAL DEADTS(NTS), MATCHED, SORTT, ALLOWED, ISA(NMIN), ISB(NMIN) 
827: C     INTEGER, PARAMETER :: extra = selected_real_kind(20,200) 
828: C     REAL(KIND=extra) :: LDUMMY, DMAT0(NCONNMAX,NMIN) 
829:       DOUBLE PRECISION :: DMAT0(2*NTS) 
830:  
831: ! 
832: !  We are setting up DMAT0(NOFF(M2),NOFF(M2)+NCOL(M2)) 
833: !  First entry for M2 at position NOFF(M2)+1 
834: !  No entries if NCOL(M2)=0 
835: !  Cycle over connected transition states using pointers for minimum M2. 
836: ! 
837:       NCOL(1:NMIN)=0 
838:       NOFF(1:NMIN)=0 
839:  
840:       fromloop: DO M2=1,NMIN    
841:          IF (M2.GT.1) NOFF(M2)=NOFF(M2-1)+NCOL(M2-1) 
842:          IF (M2.GT.1) PRINT '(A,4I10)','M2,NOFF(M2-1),NCOL(M2-1),NOFF(M2)=',M2,NOFF(M2-1),NCOL(M2-1),NOFF(M2) 
843:          IF ((DIRECTION.EQ.'AB').AND.ISA(M2).AND.(.NOT.SORTT)) THEN 
844:             CYCLE fromloop ! no escape from A 
845:          ELSEIF ((DIRECTION.EQ.'BA').AND.ISB(M2).AND.(.NOT.SORTT)) THEN 
846:             CYCLE fromloop ! no escape from B 
847:          ENDIF 
848:          J1=TOPPOINTER(M2)  !  sets J1 to the ts connected to minimum M2 with the highest value 
849:          IF (J1.LE.0) CYCLE fromloop 
850: 11       CONTINUE 
851:          IF ((.NOT.DEADTS(J1)).AND.(PLUS(J1).NE.MINUS(J1))) THEN 
852:             IF (PLUS(J1).EQ.M2) THEN 
853:                OTHER=MINUS(J1) 
854:             ELSE  
855:                OTHER=PLUS(J1) 
856:             ENDIF 
857:             MATCHED=.FALSE. 
858:             ALLOWED=.TRUE. 
859:             IF (.NOT.SORTT) THEN ! assume this is the real KMC run if SORTT is .TRUE. 
860:                IF ((DIRECTION.EQ.'BA').AND.ISA(OTHER)) THEN 
861:                   ALLOWED=.FALSE. ! no contribution from transitions into A minima 
862:                ELSEIF ((DIRECTION.EQ.'AB').AND.ISB(OTHER)) THEN 
863:                   ALLOWED=.FALSE. ! no contribution from transitions into B minima 
864:                ENDIF 
865:             ENDIF 
866:             IF (ALLOWED) THEN 
867:                IF (PLUS(J1).EQ.M2) THEN  !  M2 M1   
868:                   matchcolp: DO M1=1,NCOL(M2) 
869:                      IF (NVAL(NOFF(M2)+M1).EQ.MINUS(J1)) THEN ! a previous ts also links this pair 
870:                         DMAT0(NOFF(M2)+M1)=MIN(DMAT0(NOFF(M2)+M1)+EXP(KPLUS(J1)-KSUM(PLUS(J1))),1.0D0) 
871:                         MATCHED=.TRUE. 
872:                         EXIT matchcolp 
873:                      ENDIF 
874:                   ENDDO matchcolp 
875:                   IF (.NOT.MATCHED) THEN ! this minimum has not been connected to from M1 before 
876:                      NCOL(M2)=NCOL(M2)+1 
877:                      NVAL(NOFF(M2)+NCOL(M2))=MINUS(J1) 
878:                      DMAT0(NOFF(M2)+NCOL(M2))=MIN(EXP(KPLUS(J1)-KSUM(PLUS(J1))),1.0D0) 
879:                   ENDIF 
880: !                 IF (DEBUG) PRINT '(A,3I6,G20.10)','KMCcommit> M2,NCOL,NVAL,DMAT0=', 
881: !    &                                               M2,NCOL(M2),NVAL(NOFF(M2)+NCOL(M2)),DMAT0(NOFF(M2)+NCOL(M2) 
882:                ELSE IF (MINUS(J1).EQ.M2) THEN  !  M1 M2  
883:                   matchcolm: DO M1=1,NCOL(M2) 
884:                      IF (NVAL(NOFF(M2)+M1).EQ.PLUS(J1)) THEN ! a previous ts also links this pair 
885:                         DMAT0(NOFF(M2)+M1)=MIN(DMAT0(NOFF(M2)+M1)+EXP(KMINUS(J1)-KSUM(MINUS(J1))),1.0D0) 
886:                         MATCHED=.TRUE. 
887:                         EXIT matchcolm 
888:                      ENDIF 
889:                   ENDDO matchcolm 
890:                   IF (.NOT.MATCHED) THEN ! this minimum has not been connected to from M1 before 
891:                      NCOL(M2)=NCOL(M2)+1 
892:                      NVAL(NOFF(M2)+NCOL(M2))=PLUS(J1) 
893:                      DMAT0(NOFF(M2)+NCOL(M2))=MIN(EXP(KMINUS(J1)-KSUM(MINUS(J1))),1.0D0) 
894:                   ENDIF 
895: !                 IF (DEBUG) PRINT '(A,3I6,G20.10)','KMCcommit> M2,NCOL,NVAL,DMAT0=', 
896: !    &                                               M2,NCOL(M2),NVAL(NOFF(M2)+NCOL(M2)),DMAT0(NOFF(M2)+NCOL(M2)) 
897:                ENDIF 
898:             ENDIF 
899:          ENDIF 
900:          IF (PLUS(J1).EQ.M2) THEN 
901:             J1=POINTERP(J1) 
902:          ELSE IF (MINUS(J1).EQ.M2) THEN 
903:             J1=POINTERM(J1) 
904:          ENDIF 
905:          IF (J1.GT.0) GOTO 11 
906:       ENDDO fromloop 
907: C 
908: C  Check row normalisation. 
909: C 
910:       DO J1=1,NMIN 
911:          DUMMY=0.0D0 
912:          DO J2=1,NCOL(J1) 
913:             DUMMY=DUMMY+DMAT0(NOFF(J1)+J2) 
914:             IF (DEBUG) WRITE(*,'(A,3I8,3G20.10)') 'J1,J2,NVAL,DMAT0,sum=',J1,J2,NVAL(NOFF(J1)+J2),DMAT0(NOFF(J1)+J2),DUMMY 
915:          ENDDO 
916:          IF (DEBUG.AND.(NCOL(J1).GT.0)) WRITE(*,'(A,2I8,3G20.10)') 'J1,ncol,sum=',J1,NCOL(J1),DUMMY 
917:          IF ((NCOL(J1).GT.0).AND.(ABS(DUMMY-1.0D0).GT.1.0D-10)) THEN 
918:          PRINT*,'KMCcommit> ERROR - J1,NCOL(J1),DUMMY=',J1,NCOL(J1),DUMMY 
919:                STOP 
920:          ENDIF 
921:       ENDDO 
922:  
923:       RETURN 
924:       END 


r30873/mergedb.f90 2016-07-28 07:30:11.793525085 +0100 r30872/mergedb.f90 2016-07-28 07:30:13.393546651 +0100
134:       WRITE(*,*)'mergedb> AVOIDING MOMENT OF INITERIA CALC FOR AMH'134:       WRITE(*,*)'mergedb> AVOIDING MOMENT OF INITERIA CALC FOR AMH'
135:    ELSE135:    ELSE
136:       CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,NEWIXMIN,NEWIYMIN,NEWIZMIN)136:       CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,NEWIXMIN,NEWIYMIN,NEWIZMIN)
137:       IF ((ABS(NEWIXMIN-IXMIN(INDEX)).GT.IDIFFTOL).OR. &137:       IF ((ABS(NEWIXMIN-IXMIN(INDEX)).GT.IDIFFTOL).OR. &
138:   &    (ABS(NEWIYMIN-IYMIN(INDEX)).GT.IDIFFTOL).OR. &138:   &    (ABS(NEWIYMIN-IYMIN(INDEX)).GT.IDIFFTOL).OR. &
139:   &    (ABS(NEWIZMIN-IZMIN(INDEX)).GT.IDIFFTOL)) THEN139:   &    (ABS(NEWIZMIN-IZMIN(INDEX)).GT.IDIFFTOL)) THEN
140:          WRITE(*,'(A)') 'mergedb> possible error - principal moments of inertia do not agree with input'140:          WRITE(*,'(A)') 'mergedb> possible error - principal moments of inertia do not agree with input'
141:          WRITE(*,'(A,3F20.10)') 'mergedb> values from coordinates: ',NEWIXMIN,NEWIYMIN,NEWIZMIN141:          WRITE(*,'(A,3F20.10)') 'mergedb> values from coordinates: ',NEWIXMIN,NEWIYMIN,NEWIZMIN
142:          WRITE(*,'(A,3F20.10)') 'mergedb> values from ' // TRIM(ADJUSTL(PATHNAME)) // '/min.data', &142:          WRITE(*,'(A,3F20.10)') 'mergedb> values from ' // TRIM(ADJUSTL(PATHNAME)) // '/min.data', &
143:    &                                     IXMIN(INDEX),IYMIN(INDEX),IZMIN(INDEX)143:    &                                     IXMIN(INDEX),IYMIN(INDEX),IZMIN(INDEX)
144:                      PRINT '(A,3G20.10)',' NEWEMIN=',NEWEMIN 
145:                      PRINT '(A,2I10)',' MINDB=',MINDB 
146:                       PRINT '(A)','LOCALPOINTS:' 
147:                       PRINT '(3F20.10)',LOCALPOINTS(1:NOPT) 
148:         STOP144:         STOP
149:       ENDIF145:       ENDIF
150:    ENDIF146:    ENDIF
151: !147: !
152: !  Is it a new minimum? Set up MINMAP.148: !  Is it a new minimum? Set up MINMAP.
153: !  Testing the new minima against themselves allows us to remove duplicates from 149: !  Testing the new minima against themselves allows us to remove duplicates from 
154: !  the database we are merging! 150: !  the database we are merging! 
155: !151: !
156: !  DO J2=1,NMIN152: !  DO J2=1,NMIN
157:    DO J2=1,NMIN+NEWMIN153:    DO J2=1,NMIN+NEWMIN
158:       IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN154:       IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN
159:          DISTANCE=1.0D100155:          DISTANCE=1.0D100
160:          READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NOPT)156:          READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NOPT)
161:          CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, &157:          CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, &
162:   &                       RMAT,.FALSE.)158:   &                       RMAT,.FALSE.)
163:          IF (.NOT.(ADDPT2.OR.ADDPT3)) THEN159:          IF (.NOT.(ADDPT2.OR.ADDPT3)) THEN
164:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.GT.GEOMDIFFTOL)) THEN160:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.GT.GEOMDIFFTOL)) THEN
165:                IF (PERMDIST.OR.LPERMDIST) THEN161:                IF (PERMDIST.OR.LPERMDIST) THEN
 162:                   PRINT '(A)',' likely error?'
 163:                   PRINT '(A,3G20.10)',' NEWEMIN,EMIN(J2),DISTANCE=',NEWEMIN,EMIN(J2),DISTANCE
 164:                   PRINT '(A,2I10)',' J2,MINDB=',J2,MINDB
166:                   IF(DEBUG) THEN165:                   IF(DEBUG) THEN
167:                      PRINT '(A)',' likely error?' 
168:                      PRINT '(A,3G20.10)',' NEWEMIN,EMIN(J2),DISTANCE=',NEWEMIN,EMIN(J2),DISTANCE 
169:                      PRINT '(A,2I10)',' J2,MINDB=',J2,MINDB 
170:                       PRINT '(A)','LOCALPOINTS:'166:                       PRINT '(A)','LOCALPOINTS:'
171:                       PRINT '(3F20.10)',LOCALPOINTS(1:NOPT)167:                       PRINT '(3F20.10)',LOCALPOINTS(1:NOPT)
172:                       PRINT '(A)','LOCALPOINTS2:'168:                       PRINT '(A)','LOCALPOINTS2:'
173:                       PRINT '(3F20.10)',LOCALPOINTS2(1:NOPT)169:                       PRINT '(3F20.10)',LOCALPOINTS2(1:NOPT)
174:                   ENDIF170:                   ENDIF
175: !                 STOP !!! DJW171: !                 STOP !!! DJW
176:                ENDIF172:                ENDIF
177: !173: !
178: ! csw34> Added CHARMM structure dumping for debugging purposes174: ! csw34> Added CHARMM structure dumping for debugging purposes
179: !175: !


r30873/shannon.f90 2016-07-28 07:30:12.249531232 +0100 r30872/shannon.f90 2016-07-28 07:30:13.641549993 +0100
 20: SUBROUTINE SHANNON 20: SUBROUTINE SHANNON
 21: USE PORFUNCS 21: USE PORFUNCS
 22: USE UTILS 22: USE UTILS
 23: USE SAVESTATE 23: USE SAVESTATE
 24: USE COMMONS, ONLY : DEBUG, KAPPA, TOTALE, TEMPERATURE, ENSEMBLE, HORDERMIN, FVIBMIN, TOTALE, NMIN, EMIN, & 24: USE COMMONS, ONLY : DEBUG, KAPPA, TOTALE, TEMPERATURE, ENSEMBLE, HORDERMIN, FVIBMIN, TOTALE, NMIN, EMIN, &
 25:   &                 SHANNONTMIN, SHANNONTMAX, SHANNONTINC, EINC, NPEQ, NMINA, NMINB, LOCATIONA, LOCATIONB, & 25:   &                 SHANNONTMIN, SHANNONTMAX, SHANNONTINC, EINC, NPEQ, NMINA, NMINB, LOCATIONA, LOCATIONB, &
 26:   &                 TSATTEMPT, MAXTS, NTS, MINGROUP, FVIBTS, HORDERTS, NEGEIG, REGROUPFREETHRESH, GPFOLD, & 26:   &                 TSATTEMPT, MAXTS, NTS, MINGROUP, FVIBTS, HORDERTS, NEGEIG, REGROUPFREETHRESH, GPFOLD, &
 27:   &                 IXMIN, IYMIN, IZMIN, POINTERP, POINTERM,  MINUS, PLUS, TOPPOINTER, KMINUS,  KPLUS, & 27:   &                 IXMIN, IYMIN, IZMIN, POINTERP, POINTERM,  MINUS, PLUS, TOPPOINTER, KMINUS,  KPLUS, &
 28:   &                 ETS, PFMIN, ZSYM, PI, PFMEAN, PFTOTALA, PFTOTALB, FRICTIONT, RATEAB, RATEBA, & 28:   &                 ETS, PFMIN, ZSYM, PI, PFMEAN, PFTOTALA, PFTOTALB, FRICTIONT, RATEAB, RATEBA, &
 29:   &                 NRATESCYCLETEMPS, RATESCYCLETEMPS, NCONNMAX, SHANNONRT, & 29:   &                 NRATESCYCLETEMPS, RATESCYCLETEMPS, NCONNMAX, SHANNONRT, &
 30:   &                 SHANNONZT, EDIFFTOL 30:   &                 SHANNONZT
 31: IMPLICIT NONE 31: IMPLICIT NONE
 32: INTEGER J1, LUNIT, NGMIN, NHIGHESTPEQ(NPEQ), J2, J3, J4, NSTEPS, NDUMMY, NCONNMAXSAVE, NGMINCOPIES 32: INTEGER J1, LUNIT, NGMIN, NHIGHESTPEQ(NPEQ), J2, J3, J4, NSTEPS, NDUMMY, NCONNMAXSAVE
 33: DOUBLE PRECISION PFNORM1, PFNORM2, BARRIER(NMIN), DUMMY, FRUSTRATION, HIGHESTPEQ(NPEQ) 33: DOUBLE PRECISION PFNORM1, PFNORM2, BARRIER(NMIN), DUMMY, FRUSTRATION, HIGHESTPEQ(NPEQ)
 34: DOUBLE PRECISION DUMMY2, RATESUM1, RATESUM2, PRENORM 34: DOUBLE PRECISION DUMMY2, RATESUM1, RATESUM2, PRENORM
 35: DOUBLE PRECISION GLOBALMIN, SENTROPY, TEMPSAVE 35: DOUBLE PRECISION GLOBALMIN, SENTROPY, TEMPSAVE
 36: DOUBLE PRECISION FRICTIONFAC, PSUM 36: DOUBLE PRECISION FRICTIONFAC, PSUM
 37: INTEGER, ALLOCATABLE :: MINMAP(:) 37: INTEGER, ALLOCATABLE :: MINMAP(:)
 38:  38: 
 39: ALLOCATE(TSATTEMPT(MAXTS)) 39: ALLOCATE(TSATTEMPT(MAXTS))
 40: TEMPSAVE=TEMPERATURE 40: TEMPSAVE=TEMPERATURE
 41: ! 41: !
 42: !  Save state. 42: !  Save state.
 50: NMINASAVE=NMINA; NMINBSAVE=NMINB; NMINSAVE=NMIN; NTSSAVE=NTS; LOCATIONASAVE(1:NMINA)=LOCATIONA(1:NMINA) 50: NMINASAVE=NMINA; NMINBSAVE=NMINB; NMINSAVE=NMIN; NTSSAVE=NTS; LOCATIONASAVE(1:NMINA)=LOCATIONA(1:NMINA)
 51: LOCATIONBSAVE(1:NMINB)=LOCATIONB(1:NMINB); EMINSAVE(1:NMIN)=EMIN(1:NMIN); PFMINSAVE(1:NMIN)=PFMIN(1:NMIN) 51: LOCATIONBSAVE(1:NMINB)=LOCATIONB(1:NMINB); EMINSAVE(1:NMIN)=EMIN(1:NMIN); PFMINSAVE(1:NMIN)=PFMIN(1:NMIN)
 52: ETSSAVE(1:NTS)=ETS(1:NTS); KPLUSSAVE(1:NTS)=KPLUS(1:NTS); KMINUSSAVE(1:NTS)=KMINUS(1:NTS) 52: ETSSAVE(1:NTS)=ETS(1:NTS); KPLUSSAVE(1:NTS)=KPLUS(1:NTS); KMINUSSAVE(1:NTS)=KMINUS(1:NTS)
 53: TOPPOINTERSAVE(1:NMIN)=TOPPOINTER(1:NMIN); PLUSSAVE(1:NTS)=PLUS(1:NTS); MINUSSAVE(1:NTS)=MINUS(1:NTS) 53: TOPPOINTERSAVE(1:NMIN)=TOPPOINTER(1:NMIN); PLUSSAVE(1:NTS)=PLUS(1:NTS); MINUSSAVE(1:NTS)=MINUS(1:NTS)
 54: POINTERMSAVE(1:NTS)=POINTERM(1:NTS); POINTERPSAVE(1:NTS)=POINTERP(1:NTS) 54: POINTERMSAVE(1:NTS)=POINTERM(1:NTS); POINTERPSAVE(1:NTS)=POINTERP(1:NTS)
 55: PFMEANSAVE=PFMEAN; PFTOTALASAVE=PFTOTALA; PFTOTALBSAVE=PFTOTALB 55: PFMEANSAVE=PFMEAN; PFTOTALASAVE=PFTOTALA; PFTOTALBSAVE=PFTOTALB
 56: FVIBMINSAVE(1:NMIN)=FVIBMIN(1:NMIN); HORDERMINSAVE(1:NMIN)=HORDERMIN(1:NMIN) 56: FVIBMINSAVE(1:NMIN)=FVIBMIN(1:NMIN); HORDERMINSAVE(1:NMIN)=HORDERMIN(1:NMIN)
 57: IXMINSAVE(1:NMIN)=IXMIN(1:NMIN); IYMINSAVE(1:NMIN)=IYMIN(1:NMIN); IZMINSAVE(1:NMIN)=IZMIN(1:NMIN) 57: IXMINSAVE(1:NMIN)=IXMIN(1:NMIN); IYMINSAVE(1:NMIN)=IYMIN(1:NMIN); IZMINSAVE(1:NMIN)=IZMIN(1:NMIN)
 58:  58: 
 59: CALL GETBARRIER2(BARRIER,NGMIN,GLOBALMIN) 59: CALL GETBARRIER2(BARRIER,NGMIN,GLOBALMIN)
 60: NGMINCOPIES=0 
 61: DO J1=1,NMIN 60: DO J1=1,NMIN
 62:    IF (ABS(EMIN(J1)-GLOBALMIN).LT.EDIFFTOL) NGMINCOPIES=NGMINCOPIES+1 
 63:    IF (BARRIER(J1).LT.0.0D0) THEN 61:    IF (BARRIER(J1).LT.0.0D0) THEN
 64:      BARRIER(J1)=0.0D0 62:      BARRIER(J1)=0.0D0
 65:      CYCLE 63:      CYCLE
 66:    ENDIF 64:    ENDIF
 67:    IF (J1.EQ.NGMIN) CYCLE 65:    IF (J1.EQ.NGMIN) CYCLE
 68:    IF (ABS(EMIN(J1)-GLOBALMIN).LT.EDIFFTOL) CYCLE ! allow for degeneracy 
 69:    BARRIER(J1)=BARRIER(J1)/ABS(EMIN(J1)-GLOBALMIN) 66:    BARRIER(J1)=BARRIER(J1)/ABS(EMIN(J1)-GLOBALMIN)
 70: ENDDO 67: ENDDO
 71: PRINT '(A,I8,A)','shannon> ',NGMINCOPIES,' entries for global minimum detected in min.data' 
 72: ! 68: !
 73: ! BARRIER contains barrier divided by PE difference from lowest minimum. 69: ! BARRIER contains barrier divided by PE difference from lowest minimum.
 74: ! 70: !
 75: !  Calculate partition functions for minima. Note that the total partition function 71: !  Calculate partition functions for minima. Note that the total partition function
 76: !  is not needed, just the totals for A and B. Since A and B sets are fixed here 72: !  is not needed, just the totals for A and B. Since A and B sets are fixed here
 77: !  we don;t need to change the totals. 73: !  we don;t need to change the totals.
 78: ! 74: !
 79:        75:       
 80: GLOBALMIN=HUGE(1.0D0) 76: GLOBALMIN=HUGE(1.0D0)
 81: DO J1=1,NMIN 77: DO J1=1,NMIN
182: HIGHESTPEQ(1:NPEQ)=-1.0D0178: HIGHESTPEQ(1:NPEQ)=-1.0D0
183: PSUM=0.0D0179: PSUM=0.0D0
184: DO J1=1,NMIN180: DO J1=1,NMIN
185:    DUMMY=EXP(PFMIN(J1)-PFMIN(1)-LOG(PFNORM2))181:    DUMMY=EXP(PFMIN(J1)-PFMIN(1)-LOG(PFNORM2))
186:    DUMMY2=EXP(PFMIN(J1)-PFMIN(1)-LOG(PFNORM2))182:    DUMMY2=EXP(PFMIN(J1)-PFMIN(1)-LOG(PFNORM2))
187:    IF (DEBUG) WRITE(*,'(F20.10,I6,2G20.10)') EMIN(J1),HORDERMIN(J1),  &183:    IF (DEBUG) WRITE(*,'(F20.10,I6,2G20.10)') EMIN(J1),HORDERMIN(J1),  &
188:   &              EXP(-FVIBMIN(J1)/2.0D0-LOG(1.0D0*HORDERMIN(J1))-LOG(PFNORM1)+ FVIBMIN(1)/2.0D0 + LOG(1.0D0*HORDERMIN(1))),  &184:   &              EXP(-FVIBMIN(J1)/2.0D0-LOG(1.0D0*HORDERMIN(J1))-LOG(PFNORM1)+ FVIBMIN(1)/2.0D0 + LOG(1.0D0*HORDERMIN(1))),  &
189:   &              DUMMY185:   &              DUMMY
190:    IF(BARRIER(J1).GT.0.0D0.OR.J1.EQ.NGMIN) SENTROPY=SENTROPY+DUMMY*LOG(DUMMY)186:    IF(BARRIER(J1).GT.0.0D0.OR.J1.EQ.NGMIN) SENTROPY=SENTROPY+DUMMY*LOG(DUMMY)
191:    IF(BARRIER(J1).GT.0.0D0.OR.J1.EQ.NGMIN) PSUM=PSUM+DUMMY187:    IF(BARRIER(J1).GT.0.0D0.OR.J1.EQ.NGMIN) PSUM=PSUM+DUMMY
192:    IF (ABS(EMIN(J1)-GLOBALMIN).GT.EDIFFTOL) THEN188:    IF (J1.NE.NGMIN) FRUSTRATION=FRUSTRATION+DUMMY*BARRIER(J1)
193:       FRUSTRATION=FRUSTRATION+DUMMY*BARRIER(J1)189:    IF (J1.NE.NGMIN.AND.BARRIER(J1).GT.0.0D0) PRENORM=PRENORM+DUMMY
194:       IF (BARRIER(J1).GT.0.0D0) PRENORM=PRENORM+DUMMY 
195:    ENDIF 
196: !190: !
197: ! Maintain ordered list of the NPEQ highest probability minima at this temperature191: ! Maintain ordered list of the NPEQ highest probability minima at this temperature
198: !192: !
199:    sortloop: DO J3=1,NPEQ ! sort to find the largest barriers to product193:    sortloop: DO J3=1,NPEQ ! sort to find the largest barriers to product
200:       IF (DUMMY.GT.HIGHESTPEQ(J3)) THEN194:       IF (DUMMY.GT.HIGHESTPEQ(J3)) THEN
201:          DO J4=NPEQ,J3+1,-1195:          DO J4=NPEQ,J3+1,-1
202:             NHIGHESTPEQ(J4)=NHIGHESTPEQ(J4-1)196:             NHIGHESTPEQ(J4)=NHIGHESTPEQ(J4-1)
203:             HIGHESTPEQ(J4)=HIGHESTPEQ(J4-1)197:             HIGHESTPEQ(J4)=HIGHESTPEQ(J4-1)
204:          ENDDO198:          ENDDO
205:          NHIGHESTPEQ(J3)=J1199:          NHIGHESTPEQ(J3)=J1


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0