hdiff output

r29771/bestpath.f90 2016-01-19 17:30:07.520881068 +0000 r29770/bestpath.f90 2016-01-19 17:30:09.192903446 +0000
  1: !   CONNECT module is an implementation of a connection algorithm for finding rearrangement pathways.  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/CONNECT/bestpath.f90' in revision 29770
  2: !   Copyright (C) 2003-2006 Semen A. Trygubenko and David J. Wales 
  3: !   This file is part of CONNECT module. CONNECT module is part of OPTIM. 
  4: ! 
  5: !   OPTIM is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   OPTIM is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: MODULE BESTPATHDIJKSTRA 
 20:      IMPLICIT NONE 
 21:  
 22:      CONTAINS 
 23:  
 24:      SUBROUTINE BESTPATH 
 25:           USE PORFUNCS 
 26:           USE CONNECTDATA 
 27:           USE CONNECTUTILS,ONLY: GETDISTANCE,SETDISTANCE,GETINTERP,NREDO,MAKEMINPATHINFO,MAKETSPATHINFO 
 28:           USE KEY, ONLY : INTCONSTRAINTT, INTLJT 
 29:           USE KEYDECIDE 
 30:           USE KEYCONNECT, ONLY : IMAGEMAX 
 31:           IMPLICIT NONE 
 32:  
 33:           ! use Dijkstra's algorithm to find the shortest path between the original two endpoints. 
 34:           ! The edge weights are zero for minima connected by a single ts, infinity for pairs 
 35:           ! where we have already tried the maximum allowed connections, and the minimum distance from 
 36:           ! mind (possibly squared) otherwise. 
 37:           ! DIJKSTRA should return the number of connections to try from this shortest path, and the 
 38:           ! corresponding pairs of minima. 
 39:           ! Implementation of Dijkstra's algorithm for finding shortest paths in an UNconnected database 
 40:           ! Weights edges of graph according to mi(i)%data%D(j) or mi(i)%data%INTERP(j) (i>j) 
 41:           ! Weight=0 if connected by a single ts 
 42:           ! Weight=infinity if the maximum tryconnect attempts have been done for this pair 
 43:           ! Weight=mind(i,j) otherwise 
 44: INTEGER :: I, NDUMMY, BEGIN, END, NSTEP, J, ISTAT 
 45: TYPE MINIMUM 
 46:    INTEGER :: PREVMIN                             ! PREVIOUS MIN/TS IN SHORTEST PATH 
 47:    DOUBLE PRECISION :: D                          ! DISTANCE FROM STARTING POINT 
 48:    LOGICAL :: PERMANENT                           ! ONCE IN PERMANENT SET, THE MINIMUM's 
 49:                                                   ! label cannot be changed 
 50: END TYPE MINIMUM 
 51: TYPE(MINIMUM), DIMENSION(:), ALLOCATABLE :: DIJMIN 
 52: TYPE PATHWAY 
 53:    INTEGER :: DIJMIN, DIJTS 
 54: END TYPE PATHWAY 
 55: DOUBLE PRECISION :: WEIGHT, MINWEIGHT, DUMMYDJW, DIST 
 56: INTEGER JS, JF, J1, J2, J3, J4, NSP, NPATHTS 
 57: INTEGER, DIMENSION(:), ALLOCATABLE :: PATHMINS, PATHTS 
 58: LOGICAL, DIMENSION(:), ALLOCATABLE :: PATHSP 
 59:  
 60: 975 CONTINUE 
 61:  
 62: ! Initialise connections 
 63:    ALLOCATE(DIJMIN(NMIN)) 
 64:    DIJMIN(1:NMIN)%PREVMIN=0 
 65:    DIJMIN(1:NMIN)%D=HUGE(1.0D0) 
 66:    DIJMIN(1:NMIN)%PERMANENT=.FALSE. 
 67:  
 68: ! Start algorithm at 'begin' minimum 
 69:  
 70:    I=1 
 71:    DIJMIN(I)%D=0.0D0 
 72:    DIJMIN(I)%PERMANENT=.TRUE. 
 73:        
 74:    MAIN: DO 
 75:       DO J=2, NMIN    ! LOOP OVER ALL MINIMA  
 76:          IF (DIJMIN(J)%PERMANENT) CYCLE          ! ONLY CONSIDER TEMPORARY MINIMA 
 77:          IF (GETDISTANCE(I,J)>HUGE(1.0D0)/2) THEN ! DON't raise a HUGE number to any power > 1!! 
 78:             WEIGHT=GETDISTANCE(I,J) 
 79:          ELSE 
 80:             IF (INDEXCOSTFUNCTION) THEN  
 81:                WEIGHT=GETDISTANCE(I,J) 
 82:                IF (WEIGHT.GT.0.0D0) THEN 
 83:                   WEIGHT=ABS(I-J) 
 84:                   IF (I.EQ.2) WEIGHT=ABS((NREDO+1)/2-J) 
 85:                   IF (J.EQ.2) WEIGHT=ABS((NREDO+1)/2-I) 
 86:                ENDIF 
 87:             ELSEIF (INTERPCOSTFUNCTION) THEN  
 88:                  DIST=GETDISTANCE(I,J) 
 89:                  WEIGHT=GETINTERP(I,J)/1.0D3+DIST 
 90:                  IF (DIST.GT.DIJKSTRADMAX) THEN 
 91:                     IF (EXPCOSTFUNCTION) THEN 
 92:                        WEIGHT=EXP(MIN(WEIGHT,300.0D0)) 
 93:                     ELSE 
 94:                        WEIGHT=COSTFUNCTIONPOWER*LOG(MAX(WEIGHT,1.0D-10))  
 95:                        WEIGHT=EXP(MIN(WEIGHT,300.0D0)) 
 96:                     ENDIF 
 97:                  ENDIF 
 98:             ELSEIF (EXPCOSTFUNCTION) THEN ! SAVES MEMORY AND CPU WHEN ENDPOINT SEPARATION IS VERY LARGE SAT 
 99:                  DIST=GETDISTANCE(I,J) 
100:                  WEIGHT=DIST 
101:                  IF (DIST.GT.DIJKSTRADMAX) WEIGHT=EXP(DIST) 
102:             ELSE ! Compare squares to favour more small jumps over big ones DJW 
103:                  DIST=GETDISTANCE(I,J) 
104:                  WEIGHT=DIST 
105:                  IF (DIST.GT.DIJKSTRADMAX) WEIGHT=DIST**COSTFUNCTIONPOWER 
106:             ENDIF 
107:          ENDIF 
108:          JS=MAX(I,J) 
109:          JF=MIN(I,J) 
110:  
111:          IF ((DIJMIN(I)%D + WEIGHT) < DIJMIN(J)%D) THEN    ! UPDATE D FOR EACH CONNECTED MINIMUM 
112: !           PRINT '(A,2I8,3G20.10)','decide> I,J,DIJMIN(J)%D,DIJMIN(I)%D,WEIGHT=',I,J,DIJMIN(J)%D,DIJMIN(I)%D,WEIGHT 
113:             DIJMIN(J)%D=DIJMIN(I)%D + WEIGHT 
114:             DIJMIN(J)%PREVMIN=I 
115:          ENDIF 
116: !        PRINT '(A,2I8,3G20.10)','decide> I,J,DIJMIN(J)%D,DIJMIN(I)%D,WEIGHT=',I,J,DIJMIN(J)%D,DIJMIN(I)%D,WEIGHT 
117:       ENDDO 
118:              
119:       MINWEIGHT=HUGE(MINWEIGHT) 
120:       NDUMMY=I ! SAT: IF THIS IS THE FIRST CYCLE OF THE MAIN DO LOOP AND ALL MINIMA HAVE 
121:       ! HUGE weight ndummy will never be set which is a bit unfortunate given 
122:       ! that we use it immediately afterwards in the IF test ... 
123:       DO J=2,NMIN                           ! FIND TEMPORARY MINIMUM WITH LOWEST D 
124: !        PRINT '(A,I8,L5,2G20.10,I6)','decide> J,DIJMIN(J)%PERMANENT,DIJMIN(J)%D,MINWEIGHT,NDUMMY=', & 
125: ! &            J,DIJMIN(J)%PERMANENT,DIJMIN(J)%D,MINWEIGHT,NDUMMY 
126:          IF (DIJMIN(J)%PERMANENT) CYCLE     ! ONLY CONSIDER TEMPORARY MINIMA 
127:          IF (DIJMIN(J)%D < MINWEIGHT) THEN  ! FIND MINIMUM WEIGHT  
128:             MINWEIGHT=DIJMIN(J)%D 
129:             NDUMMY=J 
130:          ENDIF 
131:       ENDDO 
132:  
133: ! If no temporary minima with total weight less than HUGE are present, we will never exit the outer do loop! 
134:       IF (NDUMMY==I) THEN 
135:          WRITE(*,'(A)') ' bestpath> There appear to be no connections left to try in DIJKSTRA' 
136:          IMAGEMAX=MAX(1.1D0*IMAGEMAX,DBLE(IMAGEMAX+1)) 
137:          IF (IMAGEMAX.GT.10000) STOP 
138:          WRITE(*,'(A,I6,A)') ' bestpath> Increasing maximum number of images to ',IMAGEMAX,' and clearing weights' 
139:          DO J1=1,NMIN-1 
140:             DO J2=J1+1,NMIN 
141:                IF (MI(J2)%DATA%D(J1).GT.HUGE(1.0D0)/1.0D2) THEN 
142: ! 
143: ! Perhaps we should calculate the metric again? On the other hand, 
144: ! we just need a weight that enables us to continue! 
145: ! 
146:                   IF (INTERPCOSTFUNCTION) MI(J2)%DATA%INTERP(J1)=10.0D0 
147:                   MI(J2)%DATA%D(J1)=10.0D0 
148:                   MI(J2)%DATA%INTNTRIES(J1)=1 
149:                ENDIF 
150:             ENDDO 
151:          ENDDO 
152: !        CALL TSUMMARY 
153:          DEALLOCATE(DIJMIN) 
154:          CALL FLUSH(6) 
155:          GOTO 975 
156: !        STOP 
157:       ENDIF 
158:  
159:       I=NDUMMY 
160:       DIJMIN(I)%PERMANENT=.TRUE. 
161:       IF (I==2) EXIT    
162:    ENDDO MAIN 
163:  
164:    IF (MINWEIGHT.GT.HUGE(MINWEIGHT)/2) THEN 
165:       WRITE(*,'(A)') ' bestpath> There appear to be no connections left to try in DIJKSTRA' 
166:       IMAGEMAX=MAX(1.1D0*IMAGEMAX,DBLE(IMAGEMAX+1)) 
167:       IF (IMAGEMAX.GT.10000) STOP 
168:       WRITE(*,'(A,I6,A)') ' bestpath> Increasing maximum number of images to ',IMAGEMAX,' and clearing weights' 
169:       DO J1=1,NMIN-1 
170:          DO J2=J1+1,NMIN 
171:             IF (MI(J2)%DATA%D(J1).GT.HUGE(1.0D0)/1.0D2) THEN 
172: ! 
173: ! Perhaps we should calculate the metric again? On the other hand, 
174: ! we just need a weight that enables us to continue! 
175: ! 
176:                IF (INTERPCOSTFUNCTION) MI(J2)%DATA%INTERP(J1)=10.0D0 
177:                MI(J2)%DATA%D(J1)=10.0D0 
178:                MI(J2)%DATA%INTNTRIES(J1)=1 
179:             ENDIF 
180:          ENDDO 
181:       ENDDO 
182: !     CALL TSUMMARY 
183:       DEALLOCATE(DIJMIN) 
184:       CALL FLUSH(6) 
185:       GOTO 975 
186: !     STOP 
187:    ENDIF 
188:  
189:    IF (ABS(MINWEIGHT).LT.1.0D-100) THEN 
190:       WRITE(*,'(A)') ' WARNING - zero weight path already exists in DIJKSTRA' 
191:    ENDIF 
192:  
193: ! Should have found shortest path now. Get it, backwards from end 
194:  
195:   BEGIN=1 
196:   END=2 
197:   I=END 
198:   NSTEP=0 
199:   NDIJPAIRS=0 
200:   DO  
201:      IF (I==BEGIN) EXIT 
202:      NSTEP=NSTEP + 1 
203:      IF (GETDISTANCE(I,DIJMIN(I)%PREVMIN)>0.0D0) THEN 
204:         NDIJPAIRS=NDIJPAIRS+1 
205:      ENDIF 
206:      I=DIJMIN(I)%PREVMIN 
207:   ENDDO 
208:  
209:      ALLOCATE(DIJPAIR(NDIJPAIRS,2)) 
210:      ALLOCATE(DIJPAIRDIST(NDIJPAIRS)) 
211: !    IF (INDEXCOSTFUNCTION) THEN  
212: !       MINWEIGHT=MINWEIGHT 
213: !    ELSEIF (EXPCOSTFUNCTION) THEN ! SAT 
214: !         MINWEIGHT=LOG(MINWEIGHT) 
215: !    ELSEIF (INTERPCOSTFUNCTION) THEN  
216: !         MINWEIGHT=MINWEIGHT 
217: !    ELSE 
218: !         MINWEIGHT=MINWEIGHT**(1.0D0/COSTFUNCTIONPOWER) 
219: !    ENDIF 
220:      WRITE(*,'(A,I6,A,I6,A,G15.5)') ' bestpath> Shortest path in Dijkstra has ',nstep,' steps with ',NDIJPAIRS, & 
221:                                  ' missing connections, weight=',MINWEIGHT 
222:      WRITE(*,'(A)') ' bestpath> The unconnected minima in the chain and their distances are:' 
223:   I=END 
224:   NSTEP=0 
225:   NDIJPAIRS=0 
226:   ALLOCATE(PATHMINS(NMIN)) 
227:   ALLOCATE(PATHTS(NTS)) 
228:   ALLOCATE(PATHSP(NMIN+NTS)) 
229:   NSP=0 
230:   NPATHTS=0 
231:   DO  
232:      IF (I==BEGIN) EXIT 
233:      NSTEP=NSTEP + 1 
234:      NSP=NSP + 1 
235: ! add min to pathlist 
236:      PATHMINS(NSTEP)=I 
237:      PATHSP(NSP)=.TRUE. 
238:      DUMMYDJW=GETDISTANCE(I,DIJMIN(I)%PREVMIN) 
239: !    WRITE(*,'(I6,F12.2)',ADVANCE="NO") I,DUMMYDJW 
240:      IF (DUMMYDJW>0.0D0) THEN 
241: !        PRINT *, I,DUMMYDJW,DIJMIN(I)%PREVMIN 
242:         WRITE(*,'(I6,F12.2,I6)',ADVANCE="NO") I,DUMMYDJW,DIJMIN(I)%PREVMIN 
243:         NDIJPAIRS=NDIJPAIRS+1 
244:         DIJPAIR(NDIJPAIRS,1)=I 
245:         DIJPAIR(NDIJPAIRS,2)=DIJMIN(I)%PREVMIN 
246:         DIJPAIRDIST(NDIJPAIRS)=DUMMYDJW 
247:      ELSE 
248:        NPATHTS=NPATHTS + 1 
249:        NSP=NSP + 1 
250:        PATHSP(NSP)=.FALSE. 
251:        DO J4=1,NTS 
252:           IF (((TS(J4)%data%P).EQ.I).AND.((TS(J4)%data%M).EQ.DIJMIN(I)%PREVMIN)) THEN 
253:              PATHTS(NPATHTS)=J4 
254:           ELSEIF (((TS(J4)%data%M).EQ.I).AND.((TS(J4)%data%P).EQ.DIJMIN(I)%PREVMIN)) THEN 
255:              PATHTS(NPATHTS)=J4 
256:           ENDIF 
257:        ENDDO 
258:      ENDIF 
259:      I=DIJMIN(I)%PREVMIN 
260:   ENDDO 
261:   PATHMINS(NSTEP+1)=I 
262:   PATHSP(NSP+1)=.TRUE. 
263:   NSTEP=NSTEP +1 
264:   NSP=NSP + 1 
265:   NSTEP=0 
266:   NPATHTS=0 
267:   DO J3=1,NSP 
268:      IF (PATHSP(J3)) THEN 
269:         NSTEP=NSTEP + 1 
270:         CALL MAKEMINPATHINFO(PATHMINS(NSTEP)) 
271:      ELSE 
272:         NPATHTS=NPATHTS + 1 
273:         CALL MAKETSPATHINFO(PATHTS(NPATHTS)) 
274:      ENDIF 
275:   ENDDO 
276: ! WRITE(*,'(I6)') BEGIN 
277:   WRITE(*,'(A)') ' ' 
278:  
279:   DEALLOCATE(PATHMINS) 
280:   DEALLOCATE(PATHTS) 
281:   DEALLOCATE(PATHSP) 
282:   DEALLOCATE(DIJMIN) 
283:  
284:      END SUBROUTINE BESTPATH 
285:  
286: END MODULE BESTPATHDIJKSTRA  


r29771/fetchz.f 2016-01-19 17:30:08.104888885 +0000 r29770/fetchz.f 2016-01-19 17:30:09.776911261 +0000
1783:          WRITE(*,'(A,F12.5,A,E20.10)') 1783:          WRITE(*,'(A,F12.5,A,E20.10)') 
1784:      1          ' fetchz> Rate constants will be calculated for temperature ',TEMPERATURE,' and Plank;s constant=',HRED1784:      1          ' fetchz> Rate constants will be calculated for temperature ',TEMPERATURE,' and Plank;s constant=',HRED
1785:          IF (READPATH)  WRITE(*,'(A)') ' fetchz> Rate constants will be calculated for pathway saved in path.info'1785:          IF (READPATH)  WRITE(*,'(A)') ' fetchz> Rate constants will be calculated for pathway saved in path.info'
1786:       ENDIF1786:       ENDIF
1787: 1787: 
1788:       IF (PATHT) WRITE(*,'(A,I6,A)')' fetchz> Pathways will be calculated saving ',NPATHFRAME,' frames on each side'1788:       IF (PATHT) WRITE(*,'(A,I6,A)')' fetchz> Pathways will be calculated saving ',NPATHFRAME,' frames on each side'
1789:       if (machine) write(*,'(a)') ' fetchz> Will use binary files for communication'1789:       if (machine) write(*,'(a)') ' fetchz> Will use binary files for communication'
1790:       if (machine) write(*,'(a)')' fetchz> WARNING Reading binary files is not fully supported. Thouroughly tested for CHARMM only'1790:       if (machine) write(*,'(a)')' fetchz> WARNING Reading binary files is not fully supported. Thouroughly tested for CHARMM only'
1791:       IF (DUMPPATH) WRITE(*,'(A)')' fetchz> Pathway information will be printed to path.info'1791:       IF (DUMPPATH) WRITE(*,'(A)')' fetchz> Pathway information will be printed to path.info'
1792:       IF (DUMPDATAT) WRITE(*,'(A)')' fetchz> Stationary point information will be printed to min.data'1792:       IF (DUMPDATAT) WRITE(*,'(A)')' fetchz> Stationary point information will be printed to min.data'
1793:       IF (DUMPBESTPATH) WRITE(*,'(A)')' fetchz> Best pathway will be printed to path.info' 
1794:       IF (ORDERPARAMT) THEN1793:       IF (ORDERPARAMT) THEN
1795:          IF (.NOT.DUMPDATAT) THEN1794:          IF (.NOT.DUMPDATAT) THEN
1796:             PRINT '(A)','fetchz> WARNING - order parameter calculation requested, but DUMPDATA not set in odata'1795:             PRINT '(A)','fetchz> WARNING - order parameter calculation requested, but DUMPDATA not set in odata'
1797:          ELSE1796:          ELSE
1798:             WRITE(*,'(A,I8)')' fetchz> number of order parameters and derivatives to be printed=',NORDER1797:             WRITE(*,'(A,I8)')' fetchz> number of order parameters and derivatives to be printed=',NORDER
1799:             DO J1=1,NORDER1798:             DO J1=1,NORDER
1800:                WRITE(*,'(A,A,I8)')' fetchz> order parameter and additional info=',WHICHORDER(J1),ORDERNUM(J1)1799:                WRITE(*,'(A,A,I8)')' fetchz> order parameter and additional info=',WHICHORDER(J1),ORDERNUM(J1)
1801:             ENDDO1800:             ENDDO
1802:          ENDIF1801:          ENDIF
1803:       ENDIF1802:       ENDIF


r29771/ido.f90 2016-01-19 17:30:07.712883638 +0000 r29770/ido.f90 2016-01-19 17:30:09.384906015 +0000
293:      !     p => ts(i)%data%X293:      !     p => ts(i)%data%X
294:      !     nullify(ts(i)%data%X)294:      !     nullify(ts(i)%data%X)
295:      !     deallocate(p)295:      !     deallocate(p)
296:      !enddo296:      !enddo
297: END SUBROUTINE DEINITIALISE297: END SUBROUTINE DEINITIALISE
298: 298: 
299: SUBROUTINE OUTPUT299: SUBROUTINE OUTPUT
300:      USE PORFUNCS300:      USE PORFUNCS
301:      USE KEYCONNECT301:      USE KEYCONNECT
302:      USE CONNECTUTILS302:      USE CONNECTUTILS
303:      USE BESTPATHDIJKSTRA303:      USE KEY, ONLY: DUMPSP, MACHINE, DUMPALLPATHS
304:      USE KEY, ONLY: DUMPSP, MACHINE, DUMPALLPATHS, DUMPBESTPATH 
305:      IMPLICIT NONE304:      IMPLICIT NONE
306: 305: 
307:      INTEGER :: I306:      INTEGER :: I
308:      INTEGER,POINTER :: FINALPATHTS(:)307:      INTEGER,POINTER :: FINALPATHTS(:)
309:      DOUBLE PRECISION :: ESAVE,NTILDEAV308:      DOUBLE PRECISION :: ESAVE,NTILDEAV
310:      DOUBLE PRECISION :: ETSMAX,SMAX,DMAX,NTILDEMAX ! WHERE E IS ENERGY, S IS THE INTEGRATED PATHLENGTH AND D IS THE DISTANCE309:      DOUBLE PRECISION :: ETSMAX,SMAX,DMAX,NTILDEMAX ! WHERE E IS ENERGY, S IS THE INTEGRATED PATHLENGTH AND D IS THE DISTANCE
311:      DOUBLE PRECISION :: ETSMIN,SMIN,DMIN,NTILDEMIN ! BETWEEN MINIMA, AND NTILDE IS THE COOPERATIVITY INDEX310:      DOUBLE PRECISION :: ETSMIN,SMIN,DMIN,NTILDEMIN ! BETWEEN MINIMA, AND NTILDE IS THE COOPERATIVITY INDEX
312:      DOUBLE PRECISION :: BARPLUSMAX,BARPLUSMIN,BARMINUSMAX,BARMINUSMIN311:      DOUBLE PRECISION :: BARPLUSMAX,BARPLUSMIN,BARMINUSMAX,BARMINUSMIN
313:      DOUBLE PRECISION :: BAIMAX,BAIMIN,PAIMAX,PAIMIN ! BARRIER AND PATHLENGTH ASYMMETRY INDECES: {0..1}; 0 - ABSOLUTE SYMMETRY312:      DOUBLE PRECISION :: BAIMAX,BAIMIN,PAIMAX,PAIMIN ! BARRIER AND PATHLENGTH ASYMMETRY INDECES: {0..1}; 0 - ABSOLUTE SYMMETRY
314: 313: 
380:                I=I+1379:                I=I+1
381:           ENDDO380:           ENDDO
382: 381: 
383:           WRITE(*,'(1x,a,i7)') "Number of cycles               =",NConDone382:           WRITE(*,'(1x,a,i7)') "Number of cycles               =",NConDone
384:           IF (DUMPSP) CALL DODUMPSP383:           IF (DUMPSP) CALL DODUMPSP
385:           CALL MERGEXYZEOFS384:           CALL MERGEXYZEOFS
386:           IF (DUMPPATH) CALL MAKEPATHINFO385:           IF (DUMPPATH) CALL MAKEPATHINFO
387:      ELSE386:      ELSE
388:           PRINT *,'Number of allowed connection steps reached - exit newconnect'387:           PRINT *,'Number of allowed connection steps reached - exit newconnect'
389:           IF (DUMPSP) CALL DODUMPSP388:           IF (DUMPSP) CALL DODUMPSP
390:           IF (DUMPBESTPATH) CALL BESTPATH 
391:           CALL TSUMMARY389:           CALL TSUMMARY
392:           STOP390:           STOP
393:      ENDIF391:      ENDIF
394: 392: 
395:      IF (MACHINE) CALL DUMPDB(FINISHED,FINALPATHTS) ! DON't pour the water of wisdom on the dry sand393:      IF (MACHINE) CALL DUMPDB(FINISHED,FINALPATHTS) ! DON't pour the water of wisdom on the dry sand
396: END SUBROUTINE OUTPUT394: END SUBROUTINE OUTPUT
397: 395: 
398: END MODULE IDOMODULE396: END MODULE IDOMODULE


r29771/key.f90 2016-01-19 17:30:08.804897046 +0000 r29770/key.f90 2016-01-19 17:30:09.972913885 +0000
 44:      &        CHECKCONINT, ATOMMATCHFULL, NIH2LEPST,NIHEAM7T,NIHLEPST, NIHPAIRONLYT, & 44:      &        CHECKCONINT, ATOMMATCHFULL, NIH2LEPST,NIHEAM7T,NIHLEPST, NIHPAIRONLYT, &
 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, & 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, &
 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, DUMPBESTPATH 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T
 55:  55: 
 56: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 56: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 57:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 57:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 58:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 58:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 59:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 59:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 60:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 60:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 61:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads 61:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads
 62:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 62:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 63:  63: 
 64: ! hk286 > generalised THOMSON problem 64: ! hk286 > generalised THOMSON problem


r29771/keywords.f 2016-01-19 17:30:09.000900875 +0000 r29770/keywords.f 2016-01-19 17:30:10.168916507 +0000
494:          FINC=1.0D0494:          FINC=1.0D0
495:          RMSTWO=0.001D0495:          RMSTWO=0.001D0
496:          NTWO=100496:          NTWO=100
497:          NTWOITER=25497:          NTWOITER=25
498:          TWOEVAL=0.0D0498:          TWOEVAL=0.0D0
499:          PATHT=.FALSE.499:          PATHT=.FALSE.
500:          STOPFIRST=.FALSE.500:          STOPFIRST=.FALSE.
501:          CONNECTT=.FALSE.501:          CONNECTT=.FALSE.
502:          CPPNEBT=.FALSE.502:          CPPNEBT=.FALSE.
503:          DUMPPATH=.FALSE.503:          DUMPPATH=.FALSE.
504:          DUMPBESTPATH=.FALSE. 
505:          DUMPALLPATHS=.FALSE.504:          DUMPALLPATHS=.FALSE.
506:          HESSDUMPT=.FALSE.505:          HESSDUMPT=.FALSE.
507:          HESSREADT=.FALSE.506:          HESSREADT=.FALSE.
508:          INSTANTONOPTT=.FALSE.507:          INSTANTONOPTT=.FALSE.
509:          INSTANTONRATET=.FALSE.508:          INSTANTONRATET=.FALSE.
510:          INSTANTONSTARTDUMPT=.FALSE.509:          INSTANTONSTARTDUMPT=.FALSE.
511:          NIMAGEINST=1510:          NIMAGEINST=1
512:          DISTORTINST=0.4D0511:          DISTORTINST=0.4D0
513:          DELTAINST=1.D-2512:          DELTAINST=1.D-2
514:          READPATH=.FALSE.513:          READPATH=.FALSE.
2372:             IF (FILTH.EQ.0) THEN2371:             IF (FILTH.EQ.0) THEN
2373:                WRITE(PINFOSTRING,'(A9)') 'path.info'2372:                WRITE(PINFOSTRING,'(A9)') 'path.info'
2374:             ELSE2373:             ELSE
2375:                WRITE(PINFOSTRING,'(A)') 'path.info.'//TRIM(ADJUSTL(FILTHSTR))2374:                WRITE(PINFOSTRING,'(A)') 'path.info.'//TRIM(ADJUSTL(FILTHSTR))
2376:             ENDIF2375:             ENDIF
2377:             IF (MACHINE) THEN2376:             IF (MACHINE) THEN
2378:                OPEN(UNIT=88,FILE=PINFOSTRING,STATUS='UNKNOWN',FORM='UNFORMATTED')2377:                OPEN(UNIT=88,FILE=PINFOSTRING,STATUS='UNKNOWN',FORM='UNFORMATTED')
2379:             ELSE2378:             ELSE
2380:                OPEN(UNIT=88,FILE=PINFOSTRING,STATUS='UNKNOWN')2379:                OPEN(UNIT=88,FILE=PINFOSTRING,STATUS='UNKNOWN')
2381:             ENDIF2380:             ENDIF
2382: !  
2383: ! DUMPBESTPATH only writes the minima and TS on the best path as found in NEWCONNECT 
2384: ! to path.info. The format is identical to DUMPALLPATHS. 
2385: !  
2386:          ELSE IF (WORD.EQ.'DUMPBESTPATH') THEN 
2387:             DUMPBESTPATH=.TRUE. 
2388:             IF (FILTH.EQ.0) THEN 
2389:                WRITE(PINFOSTRING,'(A9)') 'path.info' 
2390:             ELSE 
2391:                WRITE(PINFOSTRING,'(A)') 'path.info.'//TRIM(ADJUSTL(FILTHSTR)) 
2392:             ENDIF 
2393:             IF (MACHINE) THEN 
2394:                OPEN(UNIT=88,FILE=PINFOSTRING,STATUS='UNKNOWN',FORM='UNFORMATTED') 
2395:             ELSE 
2396:                OPEN(UNIT=88,FILE=PINFOSTRING,STATUS='UNKNOWN') 
2397:             ENDIF 
2398: ! 2381: ! 
2399: ! Creates a file in pathsample min.data format for the minimum found2382: ! Creates a file in pathsample min.data format for the minimum found
2400: ! following a minimisation. Useful for a DPS initial path run in2383: ! following a minimisation. Useful for a DPS initial path run in
2401: ! creating entries for the two endpoints.2384: ! creating entries for the two endpoints.
2402: ! Can also be used with BHINTERP alone to generate a list of entries2385: ! Can also be used with BHINTERP alone to generate a list of entries
2403: ! for interpolated minima.2386: ! for interpolated minima.
2404: ! 2387: ! 
2405:          ELSE IF (WORD.EQ.'DUMPDATA') THEN2388:          ELSE IF (WORD.EQ.'DUMPDATA') THEN
2406:             DUMPDATAT=.TRUE.2389:             DUMPDATAT=.TRUE.
2407:             IF (FILTH.EQ.0) THEN2390:             IF (FILTH.EQ.0) THEN


r29771/ncutils.f90 2016-01-19 17:30:07.912886314 +0000 r29770/ncutils.f90 2016-01-19 17:30:09.580908637 +0000
2185:         ENDIF2185:         ENDIF
2186:         PREVIOUSTS=TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X2186:         PREVIOUSTS=TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X
2187:         DUMMY=>DUMMY%NEXT2187:         DUMMY=>DUMMY%NEXT
2188:         I=I+12188:         I=I+1
2189:      ENDDO2189:      ENDDO
2190:      CALL FLUSH(88)2190:      CALL FLUSH(88)
2191:      CLOSE(88)2191:      CLOSE(88)
2192: 2192: 
2193:      END SUBROUTINE MAKEPATHINFO2193:      END SUBROUTINE MAKEPATHINFO
2194: 2194: 
2195: ! 
2196: !  Dump min to path.info 
2197: ! 
2198:      SUBROUTINE MAKEMINPATHINFO(MINIMUM) 
2199:      USE SYMINF  
2200:      USE MODHESS 
2201:      USE MODCHARMM 
2202:      USE MODUNRES 
2203:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, & 
2204:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, METRICTENSOR ! hk286 
2205:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM 
2206:      USE PORFUNCS 
2207:      USE GENRIGID 
2208:      IMPLICIT NONE 
2209:  
2210:      CHARACTER(LEN=20) :: PINFOSTRING 
2211:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(3*NATOMS),RMS,DIAG(3*NATOMS),TEMPA(9*NATOMS),DUMQ(3*NATOMS) 
2212:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT, MINIMUM 
2213:      LOGICAL :: BTEST,KD,NNZ,NINTB,MINFRQDONE 
2214:      DOUBLE PRECISION :: QMIN(3*NATOMS),FRQSMIN(3*NATOMS),EMIN,INERTIA(3,3) 
2215:  
2216: !    LOCAL AMH VARIABLES 
2217:      INTEGER :: I_RES, GLY_COUNT 
2218:  
2219:      DOUBLE PRECISION :: TMPCOORDS(9*NATOMS/2), TMPHESS(2*NATOMS,2*NATOMS) 
2220:      DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(3*NATOMS) 
2221:  
2222:      LOGICAL KNOWE, KNOWG, KNOWH 
2223:      COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 
2224:  
2225:      MINFRQDONE=.FALSE. ! ASSUME THAT WE NEVER KNOW THE FREQUENCIES 
2226:  
2227:      EMIN=(mi(MINIMUM)%data%E) 
2228:      QMIN=(mi(MINIMUM)%data%X) 
2229:  
2230:      IF (UNRST.AND.CALCDIHE) THEN 
2231:         CALL UNRESCALCDIHEREF(DIHE,ALLANG,QMIN) 
2232:         WRITE(88,'(3F25.15)') EMIN, DIHE 
2233:      ELSE 
2234:         WRITE(88,'(F25.15)') EMIN 
2235:      ENDIF 
2236:      IF ((.NOT.MINFRQDONE).OR.TWOD.OR.CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T) THEN 
2237:         IF (UNRST) THEN ! JMC UPDATE COORDS 
2238:            DO K1=1,NRES  
2239:               C(1:3,K1)=QMIN(6*(K1-1)+1:6*(K1-1)+3) 
2240:               C(1:3,K1+NRES)=QMIN(6*(K1-1)+4:6*(K1-1)+6) 
2241:            ENDDO 
2242:            CALL UPDATEDC() 
2243:            CALL INT_FROM_CART(.TRUE.,.FALSE.) 
2244:            CALL CHAINBUILD() 
2245:         ENDIF 
2246:         IF (CHRMMT) THEN 
2247:            NCHENCALLS = 999 
2248:            HORDER=1  
2249:            FPGRP='C1' 
2250:            IF (MACHINE) THEN 
2251:               WRITE(88) HORDER,FPGRP 
2252:            ELSE 
2253:               WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2254:            ENDIF 
2255:            IF (.NOT.NOFRQS) THEN 
2256:               IF (RIGIDINIT) THEN 
2257:                  CALL GENRIGID_EIGENVALUES(QMIN, ATMASS, DIAG, INFO) 
2258:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN 
2259:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS) 
2260:                  ENDIF 
2261:                  IF (MACHINE) THEN 
2262:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2263:                  ELSE 
2264:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2265:                  ENDIF 
2266:               ELSE 
2267:                  IF (ENDNUMHESS) THEN 
2268:                     CALL MAKENUMHESS(QMIN,NATOMS) 
2269:                  ELSE 
2270:                     CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2271:                  ENDIF 
2272:                  CALL MASSWT2(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2273:                  IF (HESSDUMPT) THEN 
2274:                     LUNIT=GETUNIT() 
2275:                     OPEN(LUNIT,FILE='minhess.plus',STATUS='UNKNOWN',POSITION ='APPEND') 
2276:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2277:                     CLOSE(LUNIT) 
2278:                  ENDIF 
2279:                  CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2280:                  IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2281: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2282:                  IF (INSTANTONSTARTDUMPT) THEN 
2283:                     LUNIT=555 
2284:                     OPEN(LUNIT,file='qmrate_reactant.plus.txt', action='WRITE') 
2285:                     WRITE(LUNIT,'(a)') "Energy and Hessian eigenvalues of reactant.plus" 
2286:                     WRITE(LUNIT,*) NATOMS,NATOMS*3 
2287:                     WRITE(LUNIT,*) DUMMY1 
2288:                     WRITE(LUNIT,*) "Coordinates" 
2289:                     WRITE(LUNIT,*) QMIN 
2290:                     WRITE(LUNIT,*) "Hessian Eigenvalues" 
2291:                     WRITE(LUNIT,*) DIAG 
2292:                     WRITE(LUNIT,*) "Masses in amu (M(12C)=12)" 
2293:                     WRITE(LUNIT,*) ATMASS 
2294:                     CLOSE(LUNIT) 
2295:                  ENDIF 
2296:                  IF (MACHINE) THEN 
2297:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2298:                  ELSE 
2299:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2300:                  ENDIF 
2301:               ENDIF 
2302:            ENDIF 
2303:         ELSEIF (AMBER12T.OR.AMBERT.OR.NABT) THEN 
2304:            IF (.NOT.MACROCYCLET) THEN 
2305:               HORDER=1  
2306:               FPGRP='C1' 
2307:            ELSE 
2308:               CALL SYMMETRY(HORDER,.FALSE.,QMIN,INERTIA) 
2309:            ENDIF 
2310:            IF (MACHINE) THEN 
2311:               WRITE(88) HORDER,FPGRP 
2312:            ELSE 
2313:               WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2314:            ENDIF 
2315:            IF (.NOT.NOFRQS) THEN 
2316:               IF (RIGIDINIT) THEN 
2317:                  CALL GENRIGID_EIGENVALUES(QMIN, ATMASS, DIAG, INFO) 
2318:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN 
2319:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS) 
2320:                  ENDIF 
2321:                  IF (MACHINE) THEN 
2322:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2323:                  ELSE 
2324:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2325:                  ENDIF 
2326:               ELSE 
2327:                  IF (ENDNUMHESS.OR.AMBERT.OR.AMBER12T) THEN 
2328:                     CALL MAKENUMHESS(QMIN,NATOMS) 
2329:                  ELSE 
2330:                     CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2331:                  ENDIF 
2332:                  CALL MASSWT2(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2333:                  IF (HESSDUMPT) THEN 
2334:                     LUNIT=GETUNIT() 
2335:                     OPEN(LUNIT,FILE='minhess.plus',STATUS='UNKNOWN',POSITION ='APPEND') 
2336:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2337:                     CLOSE(LUNIT) 
2338:                  ENDIF 
2339:                  IF (FREEZE) THEN 
2340:                     CALL SWEEP_ZERO() 
2341:                     DIAG=0.0D0 
2342:                     CALL DSYEV('N','U',3*NONFREEZE,NONFROZENHESS,SIZE(NONFROZENHESS,1),DIAG(1:3*NONFREEZE),TEMPA,9*NONFREEZE,INFO) 
2343:                  ELSE 
2344:                     CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2345:                  ENDIF 
2346:                  IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2347: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2348:                  IF (INSTANTONSTARTDUMPT) THEN 
2349:                     LUNIT=555 
2350:                     OPEN(LUNIT,file='qmrate_reactant.plus.txt', action='WRITE') 
2351:                     WRITE(LUNIT,'(a)') "Energy and Hessian eigenvalues of reactant.plus" 
2352:                     WRITE(LUNIT,*) NATOMS,NATOMS*3 
2353:                     WRITE(LUNIT,*) DUMMY1 
2354:                     WRITE(LUNIT,*) "Coordinates" 
2355:                     WRITE(LUNIT,*) QMIN 
2356:                     WRITE(LUNIT,*) "Hessian Eigenvalues" 
2357:                     WRITE(LUNIT,*) DIAG 
2358:                     WRITE(LUNIT,*) "Masses in amu (M(12C)=12)" 
2359:                     WRITE(LUNIT,*) ATMASS 
2360:                     CLOSE(LUNIT) 
2361:                  ENDIF 
2362:                  IF (MACHINE) THEN 
2363:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2364:                  ELSE 
2365:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2366:                  ENDIF 
2367:               ENDIF 
2368:            ENDIF 
2369:         ELSEIF (UNRST) THEN 
2370:            HORDER=1 
2371:            FPGRP='C1' 
2372:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2373:            IF (.NOT.NOFRQS) THEN 
2374:               IF (ENDNUMHESS) THEN 
2375:                  CALL MAKENUMINTHESS(NINTS,NATOMS) 
2376:                  CALL GETSTUFF(KD,NNZ,NINTB) 
2377:                  CALL INTSECDET(QMIN,3*NATOMS,KD,NNZ,NINTB,DIAG) 
2378:               ELSE 
2379:                  CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2380:               ENDIF 
2381:               WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS) 
2382:            ENDIF 
2383:         ELSEIF (AMHT) THEN 
2384:            WRITE(88,'(I6,1X,A4)') 1,' C1' 
2385:            IF (.NOT.NOFRQS) THEN 
2386:               IF (ENDNUMHESS) THEN 
2387:                  CALL MAKENUMHESS(QMIN,NATOMS) 
2388:               ELSE 
2389:                  CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2390:               ENDIF 
2391:               CALL MASSWT(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2392:                  IF (HESSDUMPT) THEN 
2393:                     LUNIT=GETUNIT() 
2394:                     OPEN(LUNIT,FILE='minhess.plus',STATUS='UNKNOWN',POSITION ='APPEND') 
2395:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2396:                     CLOSE(LUNIT) 
2397:                  ENDIF 
2398:               CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2399:               IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2400:               WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS) 
2401: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2402:                     IF (INSTANTONSTARTDUMPT) THEN 
2403: !                      CALL POTENTIAL(QPLUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2404:                       LUNIT=555 
2405:                       open(LUNIT,file='qmrate_reactant.plus.txt', action='write') 
2406:                       write(LUNIT,'(a)') "Energy and Hessian eigenvalues of reactant.plus" 
2407:                       write(LUNIT,*) NATOMS,NATOMS*3 
2408:                       write(LUNIT,*) DUMMY1 
2409:                       write(LUNIT,*) "Coordinates" 
2410:                       write(LUNIT,*) QMIN 
2411:                       write(LUNIT,*) "Hessian Eigenvalues" 
2412:                       write(LUNIT,*) DIAG 
2413:                       write(LUNIT,*) "Masses in amu (M(12C)=12)" 
2414:                       write(LUNIT,*) ATMASS 
2415:                       close(LUNIT) 
2416:                     ENDIF 
2417:            ENDIF              
2418: ! hk286 
2419:         ELSEIF (GTHOMSONT) THEN 
2420:            CALL GTHOMSONANGTOC(TMPCOORDS,QMIN,NATOMS) 
2421:            CALL SYMMETRY(HORDER,.FALSE.,TMPCOORDS,INERTIA) 
2422:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2423:            IF (.NOT.NOFRQS) THEN 
2424:               IF (ENDNUMHESS) THEN 
2425:                  CALL MAKENUMHESS(QMIN,NATOMS) 
2426:               ELSE 
2427:                  CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2428:               ENDIF 
2429:               CALL MASSWT(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2430:               IF (HESSDUMPT) THEN 
2431:                  LUNIT=GETUNIT() 
2432:                  OPEN(LUNIT,FILE='minhess.plus',STATUS='UNKNOWN',POSITION ='APPEND') 
2433:                  WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2434:                  CLOSE(LUNIT) 
2435:               ENDIF 
2436:               CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2437:               IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2438:            ENDIF 
2439:            WRITE(88,'(3G20.10)') (1.0D10, J2=1, NATOMS) 
2440:            IF (SDT.OR.TTM3T) THEN 
2441:               WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,2*NATOMS) 
2442:            ELSEIF (BOWMANT) THEN 
2443:               WRITE(88,'(3G20.10)') (DIAG(J2)*2625.47D26,J2=1,2*NATOMS) 
2444:            ELSE 
2445:               WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,2*NATOMS) 
2446:            ENDIF 
2447:             
2448:         ELSE 
2449:            CALL SYMMETRY(HORDER,.FALSE.,QMIN,INERTIA) 
2450:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2451:            IF (.NOT.NOFRQS) THEN 
2452:               ! sn402: Currently there are two different methods implemented for finding the normal modes of 
2453:               ! local rigid bodies. GENRIGID_NORMALMODES makes use of the metric tensor formulation and so should 
2454:               ! in principle be more accurate. Eventually this should be made the default (or indeed only) option 
2455:               ! and the keyword METRICTENSOR should be removed. 
2456:               IF (RIGIDINIT) THEN 
2457:                  IF(METRICTENSOR) THEN 
2458:                      write(*,*) "ATMASS going in:" 
2459:                      write(*,*) ATMASS 
2460:                      CALL GENRIGID_NORMALMODES(QMIN, ATMASS, DIAG, INFO) 
2461:                  ELSE 
2462:                      CALL GENRIGID_EIGENVALUES(QMIN, ATMASS, DIAG, INFO) 
2463:                  ENDIF 
2464:  
2465:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN 
2466:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS) 
2467:                  ENDIF 
2468: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame 
2469:               ELSE IF (RBAAT) THEN 
2470:                  RBAANORMALMODET = .TRUE. 
2471:                  CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2472:                  CALL NRMLMD (QMIN, DIAG, .FALSE.) 
2473:                  RBAANORMALMODET = .FALSE. 
2474:               ELSE 
2475:                  IF (ENDNUMHESS) THEN 
2476:                     CALL MAKENUMHESS(QMIN,NATOMS) 
2477:                  ELSE 
2478:                     CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2479:                  ENDIF 
2480:                  CALL MASSWT(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2481:                  IF (HESSDUMPT) THEN 
2482:                     LUNIT=GETUNIT() 
2483:                     OPEN(LUNIT,FILE='minhess.plus',STATUS='UNKNOWN',POSITION ='APPEND') 
2484:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2485:                     CLOSE(LUNIT) 
2486:                  ENDIF 
2487:                  CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2488:                  IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2489: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2490:                  IF (INSTANTONSTARTDUMPT) THEN 
2491:                     LUNIT=555 
2492:                     OPEN(LUNIT,file='qmrate_reactant.plus.txt', action='WRITE') 
2493:                     WRITE(LUNIT,'(a)') "Energy and Hessian eigenvalues of reactant.plus" 
2494:                     WRITE(LUNIT,*) NATOMS,NATOMS*3 
2495:                     WRITE(LUNIT,*) DUMMY1 
2496:                     WRITE(LUNIT,*) "Coordinates" 
2497:                     WRITE(LUNIT,*) QMIN 
2498:                     WRITE(LUNIT,*) "Hessian Eigenvalues" 
2499:                     WRITE(LUNIT,*) DIAG 
2500:                     WRITE(LUNIT,*) "Masses in amu (M(12C)=12)" 
2501:                     WRITE(LUNIT,*) ATMASS 
2502:                     CLOSE(LUNIT) 
2503:                  ENDIF 
2504:               ENDIF 
2505:               IF (SDT.OR.TTM3T) THEN 
2506:                  WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2507:               ELSEIF (BOWMANT) THEN 
2508:                  WRITE(88,'(3G20.10)') (DIAG(J2)*2625.47D26,J2=1,3*NATOMS) 
2509:               ELSEIF (RIGIDINIT) THEN 
2510:                  IF (MACHINE) THEN 
2511: !                    WRITE(88) (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2512:                     WRITE(88) (DIAG(J2),J2=1,DEGFREEDOMS) 
2513:                  ELSE 
2514:                     WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,DEGFREEDOMS) 
2515: !                    WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2516:                  ENDIF 
2517:               ELSE 
2518:                  WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS) 
2519:               ENDIF 
2520:            ENDIF 
2521:         ENDIF 
2522:      ELSE 
2523:         CALL SYMMETRY(HORDER,.FALSE.,QMIN,INERTIA) 
2524:         WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2525:      ENDIF 
2526:      IF (MACHINE) THEN 
2527:         IF (GTHOMSONT) THEN 
2528:            CALL GTHOMSONANGTOC(TMPCOORDS, QMIN, NATOMS) 
2529:            WRITE(88) (TMPCOORDS, J2=1, 3*NATOMS) 
2530:         ELSE 
2531:            WRITE(88) (QMIN,J2=1,3*NATOMS) 
2532:         ENDIF 
2533:      ELSEIF (AMHT) THEN 
2534:  
2535: !  THIS IS FOR PLACE HOLDING C-BETAS FOR GLYCINE IN AMH 
2536:         GLY_COUNT = 0 
2537:  
2538:         DO J2=1, NRES_AMH_TEMP 
2539:            IF (SEQ(J2).EQ.8) THEN 
2540: !             WRITE(2,*)SEQ(J2) , J2 
2541:                WRITE(88,*)QMIN(9*(J2-1)+1-GLY_COUNT*3), & 
2542:                QMIN(9*(J2-1)+2-GLY_COUNT*3),QMIN(9*(J2-1)+3-GLY_COUNT*3) 
2543:               WRITE(88,*)QMIN(9*(J2-1)+1-GLY_COUNT*3), & 
2544:                QMIN(9*(J2-1)+2-GLY_COUNT*3),QMIN(9*(J2-1)+3-GLY_COUNT*3) 
2545:               WRITE(88,*)QMIN(9*(J2-1)+4-GLY_COUNT*3), & 
2546:                 QMIN(9*(J2-1)+5-GLY_COUNT*3),QMIN(9*(J2-1)+6-GLY_COUNT*3) 
2547:               GLY_COUNT = GLY_COUNT + 1 
2548:            ELSE 
2549: !            WRITE(2,*)SEQ(J2) , J2 
2550:              WRITE(88,*)QMIN(9*(J2-1)+1-GLY_COUNT*3), & 
2551:                QMIN(9*(J2-1)+2-GLY_COUNT*3),QMIN(9*(J2-1)+3-GLY_COUNT*3) 
2552:              WRITE(88,*)QMIN(9*(J2-1)+4-GLY_COUNT*3), & 
2553:                QMIN(9*(J2-1)+5-GLY_COUNT*3),QMIN(9*(J2-1)+6-GLY_COUNT*3) 
2554:             WRITE(88,*)QMIN(9*(J2-1)+7-GLY_COUNT*3), & 
2555:                QMIN(9*(J2-1)+8-GLY_COUNT*3),QMIN(9*(J2-1)+9-GLY_COUNT*3) 
2556:            ENDIF 
2557:         ENDDO 
2558:      ELSE 
2559:          IF (GTHOMSONT) THEN 
2560:             CALL GTHOMSONANGTOC(TMPCOORDS, QMIN, NATOMS) 
2561:             WRITE(88,'(3F25.15)') (TMPCOORDS(J2), J2=1, 3*NATOMS) 
2562:          ELSE                        
2563:             WRITE(88,'(3F25.15)') (QMIN(J2),J2=1,3*NATOMS) 
2564:          ENDIF 
2565:      ENDIF 
2566:      END SUBROUTINE MAKEMINPATHINFO 
2567:  
2568: ! 
2569: !  Dump TS to path.info  
2570: ! 
2571:      SUBROUTINE MAKETSPATHINFO(TSN) 
2572:  
2573:      USE SYMINF  
2574:      USE MODHESS 
2575:      USE MODCHARMM 
2576:      USE MODUNRES 
2577:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, & 
2578:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, METRICTENSOR ! hk286 
2579:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM 
2580:      USE PORFUNCS 
2581:      USE GENRIGID 
2582:      IMPLICIT NONE 
2583:  
2584:      CHARACTER(LEN=20) :: PINFOSTRING 
2585:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(3*NATOMS),RMS,DIAG(3*NATOMS),TEMPA(9*NATOMS),DUMQ(3*NATOMS) 
2586:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT,TSN 
2587:      LOGICAL :: BTEST,KD,NNZ,NINTB,TSFRQDONE 
2588:      DOUBLE PRECISION :: QTS(3*NATOMS),FRQSTS(3*NATOMS),ETS,INERTIA(3,3) 
2589:  
2590: !    LOCAL AMH VARIABLES 
2591:      INTEGER :: I_RES, GLY_COUNT 
2592:  
2593:      DOUBLE PRECISION :: TMPCOORDS(9*NATOMS/2), TMPHESS(2*NATOMS,2*NATOMS) 
2594:      DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(3*NATOMS) 
2595:  
2596:      LOGICAL KNOWE, KNOWG, KNOWH 
2597:      COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 
2598:  
2599:      TSFRQDONE=.FALSE.  ! ASSUME THAT WE NEVER KNOW THE FREQUENCIES 
2600:       
2601:      ETS=(TS(TSN)%data%E) 
2602:      QTS=(TS(TSN)%data%X) 
2603:  
2604:      IF (MACHINE) THEN 
2605:           WRITE(88) ETS 
2606:      ELSE 
2607:           WRITE(88,'(F25.15)') ETS 
2608:      ENDIF 
2609:      IF ((.NOT.TSFRQDONE).OR.TWOD.OR.CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T) THEN 
2610:         IF (UNRST) THEN ! JMC UPDATE COORDS 
2611:            DO K1=1,NRES  
2612:               C(1:3,K1)=QTS(6*(K1-1)+1:6*(K1-1)+3) 
2613:               C(1:3,K1+NRES)=QTS(6*(K1-1)+4:6*(K1-1)+6) 
2614:            ENDDO 
2615:            CALL UPDATEDC() 
2616:            CALL INT_FROM_CART(.TRUE.,.FALSE.) 
2617:            CALL CHAINBUILD() 
2618:         ENDIF 
2619:         IF (CHRMMT) THEN 
2620:            NCHENCALLS = 999 
2621:            HORDER=1 
2622:            FPGRP='C1' 
2623:            IF (MACHINE) THEN 
2624:               WRITE(88) HORDER,FPGRP 
2625:            ELSE 
2626:               WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2627:            ENDIF 
2628:            IF (.NOT.NOFRQS) THEN 
2629:               IF (RIGIDINIT) THEN 
2630: ! hk286 - TS is recorded in rigid body coordinates 
2631:                  ATOMRIGIDCOORDT = .FALSE. 
2632:                  CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO) 
2633:                  ATOMRIGIDCOORDT = .TRUE. 
2634:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN 
2635:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS) 
2636:                  ENDIF 
2637:                  IF (MACHINE) THEN 
2638:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2639:                  ELSE 
2640:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2641:                  ENDIF 
2642:               ELSE 
2643:                  IF (ENDNUMHESS) THEN 
2644:                     CALL MAKENUMHESS(QTS,NATOMS) 
2645:                  ELSE 
2646:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2647:                  ENDIF 
2648:                  CALL MASSWT2(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2649:                  IF (HESSDUMPT) THEN 
2650:                     LUNIT=GETUNIT() 
2651:                     OPEN(LUNIT,FILE='tshess',STATUS='UNKNOWN',POSITION ='APPEND') 
2652:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2653:                     CLOSE(LUNIT) 
2654:                  ENDIF 
2655:                  CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2656:                  IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2657: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2658:                  IF (INSTANTONSTARTDUMPT) THEN 
2659:                     LUNIT=555 
2660:                     OPEN(LUNIT,file='qmrate_ts.txt', action='WRITE') 
2661:                     WRITE(LUNIT,'(a)') "Energy and Hessian eigenvalues of transition state" 
2662:                     WRITE(LUNIT,*) NATOMS,NATOMS*3 
2663:                     WRITE(LUNIT,*) DUMMY1 
2664:                     WRITE(LUNIT,*) "Coordinates" 
2665:                     WRITE(LUNIT,*) QTS 
2666:                     WRITE(LUNIT,*) "Hessian Eigenvalues" 
2667:                     WRITE(LUNIT,*) DIAG 
2668:                     WRITE(LUNIT,*) "Masses in amu (M(12C)=12)" 
2669:                     WRITE(LUNIT,*) ATMASS 
2670:                     CLOSE(LUNIT) 
2671:                  ENDIF 
2672:                  IF (MACHINE) THEN 
2673:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2674:                  ELSE 
2675:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2676:                  ENDIF 
2677:               ENDIF 
2678:            ENDIF 
2679:         ELSE IF (AMBER12T.OR.AMBERT.OR.NABT) THEN 
2680:            IF (.NOT.MACROCYCLET) THEN 
2681:               HORDER=1 
2682:               FPGRP='C1' 
2683:            ELSE 
2684:               CALL SYMMETRY(HORDER,.FALSE.,QTS,INERTIA) 
2685:            ENDIF 
2686:            IF (MACHINE) THEN 
2687:               WRITE(88) HORDER,FPGRP 
2688:            ELSE 
2689:               WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2690:            ENDIF 
2691:            IF (.NOT.NOFRQS) THEN 
2692:               IF (RIGIDINIT) THEN 
2693: ! h286 - TS is saved in rigid body coordinates 
2694:                  ATOMRIGIDCOORDT = .FALSE. 
2695:                  CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO) 
2696:                  ATOMRIGIDCOORDT = .TRUE. 
2697:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN 
2698:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS) 
2699:                  ENDIF 
2700:                  IF (MACHINE) THEN 
2701:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2702:                  ELSE 
2703:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2704:                  ENDIF 
2705:               ELSE 
2706:                  IF (ENDNUMHESS.OR.AMBERT.OR.AMBER12T) THEN 
2707:                     CALL MAKENUMHESS(QTS,NATOMS) 
2708:                  ELSE 
2709:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2710:                  ENDIF 
2711:                  CALL MASSWT2(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2712:                  IF (HESSDUMPT) THEN 
2713:                     LUNIT=GETUNIT() 
2714:                     OPEN(LUNIT,FILE='tshess',STATUS='UNKNOWN',POSITION ='APPEND') 
2715:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2716:                     CLOSE(LUNIT) 
2717:                  ENDIF 
2718:                  IF (FREEZE) THEN 
2719:                     CALL SWEEP_ZERO() 
2720:                     DIAG=0.0D0 
2721:                     CALL DSYEV('N','U',3*NONFREEZE,NONFROZENHESS,SIZE(NONFROZENHESS,1),DIAG(1:3*NONFREEZE),TEMPA,9*NONFREEZE,INFO) 
2722:                  ELSE 
2723:                     CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2724:                  ENDIF 
2725:                  IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2726: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2727:                  IF (INSTANTONSTARTDUMPT) THEN 
2728:                     LUNIT=555 
2729:                     OPEN(LUNIT,file='qmrate_ts.txt', action='WRITE') 
2730:                     WRITE(LUNIT,'(a)') "Energy and Hessian eigenvalues of transition state" 
2731:                     WRITE(LUNIT,*) NATOMS,NATOMS*3 
2732:                     WRITE(LUNIT,*) DUMMY1 
2733:                     WRITE(LUNIT,*) "Coordinates" 
2734:                     WRITE(LUNIT,*) QTS 
2735:                     WRITE(LUNIT,*) "Hessian Eigenvalues" 
2736:                     WRITE(LUNIT,*) DIAG 
2737:                     WRITE(LUNIT,*) "Masses in amu (M(12C)=12)" 
2738:                     WRITE(LUNIT,*) ATMASS 
2739:                     CLOSE(LUNIT) 
2740:                  ENDIF 
2741:                  IF (MACHINE) THEN 
2742:                     WRITE(88) (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2743:                  ELSE 
2744:                     WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2745:                  ENDIF 
2746:               ENDIF 
2747:            ENDIF 
2748:         ELSEIF (UNRST) THEN 
2749:            HORDER=1 
2750:            FPGRP='C1' 
2751:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2752:            IF (.NOT.NOFRQS) THEN 
2753:               IF (ENDNUMHESS) THEN 
2754:                  CALL MAKENUMINTHESS(NINTS,NATOMS) 
2755:                  CALL GETSTUFF(KD,NNZ,NINTB) 
2756:                  CALL INTSECDET(QTS,3*NATOMS,KD,NNZ,NINTB,DIAG) 
2757:               ELSE 
2758:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2759:               ENDIF 
2760:               DO J2=1,NINTS-1 
2761:                  IF (DIAG(J2).LT.0.0D0) PRINT *,'Higher order saddle found in pathway - ts eigenvalue ',DIAG(J2) 
2762:               END DO 
2763:               WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS) 
2764:            ENDIF 
2765:         ELSEIF (AMHT) THEN 
2766:            WRITE(88,'(I6,1X,A4)') 1,' C1' 
2767:            IF (.NOT.NOFRQS) THEN 
2768:               IF (ENDNUMHESS) THEN 
2769:                  CALL MAKENUMHESS(QTS,NATOMS) 
2770:               ELSE 
2771:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2772:               ENDIF 
2773:               CALL MASSWT(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2774:                  IF (HESSDUMPT) THEN 
2775:                     LUNIT=GETUNIT() 
2776:                     OPEN(LUNIT,FILE='tshess',STATUS='UNKNOWN',POSITION ='APPEND') 
2777:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2778:                     CLOSE(LUNIT) 
2779:                  ENDIF 
2780:               CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2781:               IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2782:               WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS) 
2783: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2784:                     IF (INSTANTONSTARTDUMPT) THEN 
2785: !                      CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2786:                       LUNIT=555 
2787:                       open(LUNIT,file='qmrate_ts.txt', action='write') 
2788:                       write(LUNIT,'(a)') "Energy and Hessian eigenvalues of transition state" 
2789:                       write(LUNIT,*) NATOMS,NATOMS*3 
2790:                       write(LUNIT,*) DUMMY1 
2791:                       write(LUNIT,*) "Coordinates" 
2792:                       write(LUNIT,*) QTS 
2793:                       write(LUNIT,*) "Hessian Eigenvalues" 
2794:                       write(LUNIT,*) DIAG 
2795:                       write(LUNIT,*) "Masses in amu (M(12C)=12)" 
2796:                       write(LUNIT,*) ATMASS 
2797:                       close(LUNIT) 
2798:                     ENDIF 
2799:            ENDIF 
2800:  
2801:         ELSEIF (GTHOMSONT) THEN 
2802:            CALL GTHOMSONANGTOC(TMPCOORDS, QTS, NATOMS) 
2803:            CALL SYMMETRY(HORDER,.FALSE.,TMPCOORDS,INERTIA) 
2804:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2805:            IF (.NOT.NOFRQS) THEN 
2806:               IF (ENDNUMHESS) THEN 
2807:                  CALL MAKENUMHESS(QTS,NATOMS) 
2808:               ELSE 
2809:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2810:               ENDIF 
2811:               CALL MASSWT(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2812:                  IF (HESSDUMPT) THEN 
2813:                     LUNIT=GETUNIT() 
2814:                     OPEN(LUNIT,FILE='tshess',STATUS='UNKNOWN',POSITION ='APPEND') 
2815:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2816:                     CLOSE(LUNIT) 
2817:                  ENDIF 
2818:               CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2819:               IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2820:               IF (DIAG(3*NATOMS) < 0.0) THEN 
2821:                  DIAG(2*NATOMS) = DIAG(3*NATOMS) 
2822:               ENDIF 
2823:  
2824:            ENDIF 
2825:            WRITE(88,'(3G20.10)') (1.0D10, J2=1, NATOMS) 
2826:            IF (SDT.OR.TTM3T) THEN 
2827:               WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,2*NATOMS) 
2828:            ELSEIF (BOWMANT) THEN 
2829:               WRITE(88,'(3G20.10)') (DIAG(J2)*2625.47D26,J2=1,2*NATOMS) 
2830:            ELSE 
2831:               WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,2*NATOMS) 
2832:            ENDIF 
2833:  
2834:         ELSE 
2835:            CALL SYMMETRY(HORDER,.FALSE.,QTS,INERTIA) 
2836:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2837:            IF (.NOT.NOFRQS) THEN 
2838:               IF (RIGIDINIT) THEN 
2839: ! hk286 - TS is recorded in rigid body coordinates 
2840:                  ATOMRIGIDCOORDT = .FALSE. 
2841:                  IF(METRICTENSOR) THEN 
2842:                      CALL GENRIGID_NORMALMODES(QTS, ATMASS, DIAG, INFO) 
2843:                  ELSE 
2844:                      CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO) 
2845:                  ENDIF 
2846:                  ATOMRIGIDCOORDT = .TRUE. 
2847:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN 
2848:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS) 
2849:                  ENDIF 
2850: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame 
2851:               ELSE IF (RBAAT) THEN 
2852:                  RBAANORMALMODET = .TRUE. 
2853:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2854:                  CALL NRMLMD (QTS, DIAG, .FALSE.) 
2855:                  RBAANORMALMODET = .FALSE. 
2856:               ELSE 
2857:                  IF (ENDNUMHESS) THEN 
2858:                     CALL MAKENUMHESS(QTS,NATOMS) 
2859:                  ELSE 
2860:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.) 
2861:                  ENDIF 
2862:                  CALL MASSWT(NATOMS,ATMASS,DUMQ,GRAD,.TRUE.) 
2863:                  IF (HESSDUMPT) THEN 
2864:                     LUNIT=GETUNIT() 
2865:                     OPEN(LUNIT,FILE='tshess',STATUS='UNKNOWN',POSITION ='APPEND') 
2866:                     WRITE(LUNIT,'(6G20.10)') HESS(1:3*NATOMS,1:3*NATOMS) 
2867:                     CLOSE(LUNIT) 
2868:                  ENDIF 
2869:                  CALL DSYEV('N','U',3*NATOMS,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO) 
2870:                  IF (DIAG(1).LT.DIAG(3*NATOMS)) CALL EIGENSORT_VAL_ASC(DIAG,HESS,3*NATOMS,3*NATOMS) 
2871: ! jbr36 - writes the first input for qm rate calculations from classical rates 
2872:                  IF (INSTANTONSTARTDUMPT) THEN 
2873:                     LUNIT=555 
2874:                     OPEN(LUNIT,file='qmrate_ts.txt', action='WRITE') 
2875:                     WRITE(LUNIT,'(a)') "Energy and Hessian eigenvalues of transition state" 
2876:                     WRITE(LUNIT,*) NATOMS,NATOMS*3 
2877:                     WRITE(LUNIT,*) DUMMY1 
2878:                     WRITE(LUNIT,*) "Coordinates" 
2879:                     WRITE(LUNIT,*) QTS 
2880:                     WRITE(LUNIT,*) "Hessian Eigenvalues" 
2881:                     WRITE(LUNIT,*) DIAG 
2882:                     WRITE(LUNIT,*) "Masses in amu (M(12C)=12)" 
2883:                     WRITE(LUNIT,*) ATMASS 
2884:                     CLOSE(LUNIT) 
2885:                  ENDIF 
2886:               ENDIF 
2887:               IF (SDT.OR.TTM3T) THEN 
2888:                  WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,3*NATOMS) 
2889:               ELSEIF (BOWMANT) THEN 
2890:                  WRITE(88,'(3G20.10)') (DIAG(J2)*2625.47D26,J2=1,3*NATOMS) 
2891:               ELSEIF (RIGIDINIT) THEN 
2892:                  IF (MACHINE) THEN 
2893: !                    WRITE(88) (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2894:                     WRITE(88) (DIAG(J2),J2=1,DEGFREEDOMS) 
2895:                  ELSE 
2896:                     WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,DEGFREEDOMS) 
2897: !                    WRITE(88,'(3G20.10)') (DIAG(J2)*4.184D26,J2=1,DEGFREEDOMS) 
2898:                  ENDIF 
2899:               ELSE 
2900:                  WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS) 
2901:               ENDIF 
2902:            ENDIF 
2903:         ENDIF 
2904:      ELSE 
2905:         CALL SYMMETRY(HORDER,.FALSE.,QTS,INERTIA) 
2906:         WRITE(88,'(I6,1X,A4)') HORDER,FPGRP 
2907:      ENDIF 
2908:      IF (MACHINE) THEN 
2909:         IF (GTHOMSONT) THEN 
2910:            CALL GTHOMSONANGTOC(TMPCOORDS, QTS, NATOMS) 
2911:            WRITE(88) (TMPCOORDS(J2), J2=1, 3*NATOMS) 
2912:         ELSEIF (RIGIDINIT) THEN 
2913:            CALL TRANSFORMRIGIDTOC (1, NRIGIDBODY, XCOORDS, QTS(1:DEGFREEDOMS)) 
2914:            WRITE(88) (XCOORDS(J2),J2=1,NOPT) 
2915:         ELSE            
2916:            WRITE(88) (QTS(J2),J2=1,NOPT) 
2917:         ENDIF 
2918:      ELSEIF (AMHT) THEN 
2919: !       READ SEQUENCE 
2920:  
2921: !  THIS IS FOR PLACE HOLDING C-BETAS FOR GLYCINE IN AMH 
2922:         GLY_COUNT = 0 
2923:  
2924:         DO J2=1, NRES_AMH_TEMP 
2925:            IF (SEQ(J2).EQ.8) THEN 
2926: !             WRITE(2,*)SEQ(J2) , J2 
2927:                WRITE(88,*)QTS(9*(J2-1)+1-GLY_COUNT*3), & 
2928:                 QTS(9*(J2-1)+2-GLY_COUNT*3),QTS(9*(J2-1)+3-GLY_COUNT*3) 
2929:               WRITE(88,*)QTS(9*(J2-1)+1-GLY_COUNT*3), & 
2930:                 QTS(9*(J2-1)+2-GLY_COUNT*3),QTS(9*(J2-1)+3-GLY_COUNT*3) 
2931:               WRITE(88,*)QTS(9*(J2-1)+4-GLY_COUNT*3), & 
2932:                 QTS(9*(J2-1)+5-GLY_COUNT*3),QTS(9*(J2-1)+6-GLY_COUNT*3) 
2933:               GLY_COUNT = GLY_COUNT + 1 
2934:            ELSE 
2935: !            WRITE(2,*)SEQ(J2) , J2 
2936:              WRITE(88,*)QTS(9*(J2-1)+1-GLY_COUNT*3), & 
2937:                QTS(9*(J2-1)+2-GLY_COUNT*3),QTS(9*(J2-1)+3-GLY_COUNT*3) 
2938:              WRITE(88,*)QTS(9*(J2-1)+4-GLY_COUNT*3), & 
2939:                QTS(9*(J2-1)+5-GLY_COUNT*3),QTS(9*(J2-1)+6-GLY_COUNT*3) 
2940:              WRITE(88,*)QTS(9*(J2-1)+7-GLY_COUNT*3), & 
2941:                QTS(9*(J2-1)+8-GLY_COUNT*3),QTS(9*(J2-1)+9-GLY_COUNT*3) 
2942:            ENDIF 
2943:          ENDDO 
2944:      ELSE 
2945:         IF (GTHOMSONT) THEN 
2946:            CALL GTHOMSONANGTOC(TMPCOORDS, QTS, NATOMS) 
2947:            WRITE(88,'(3F25.15)') (TMPCOORDS(J2), J2=1, 3*NATOMS) 
2948:         ELSEIF (RIGIDINIT) THEN 
2949:            CALL TRANSFORMRIGIDTOC (1, NRIGIDBODY, XCOORDS, QTS(1:DEGFREEDOMS)) 
2950:            WRITE(88,'(3F25.15)') (XCOORDS(J2),J2=1,3*NATOMS) 
2951:         ELSE                       
2952:            WRITE(88,'(3F25.15)') (QTS(J2),J2=1,NOPT) 
2953:         ENDIF 
2954:      ENDIF 
2955:      END SUBROUTINE MAKETSPATHINFO 
2956:  
2957: ! --------------------------------------------------------------------------------------------------------------------2195: ! --------------------------------------------------------------------------------------------------------------------
2958: 2196: 
2959: !2197: !
2960: !  Dump the latest min-sad-min triple to path.info in the usual format2198: !  Dump the latest min-sad-min triple to path.info in the usual format
2961: !  2199: !  
2962:      SUBROUTINE MAKEALLPATHINFO(QTS,QPLUS,QMINUS,ETS,EPLUS,EMINUS,FRQSTS,FRQSPLUS,FRQSMINUS)2200:      SUBROUTINE MAKEALLPATHINFO(QTS,QPLUS,QMINUS,ETS,EPLUS,EMINUS,FRQSTS,FRQSPLUS,FRQSMINUS)
2963:      USE SYMINF 2201:      USE SYMINF 
2964:      USE MODHESS2202:      USE MODHESS
2965:      USE MODCHARMM2203:      USE MODCHARMM
2966:      USE MODUNRES2204:      USE MODUNRES


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0