hdiff output

r29929/addperm.f 2016-02-09 15:30:11.638671532 +0000 r29928/addperm.f 2016-02-09 15:30:13.922701716 +0000
 24: C  One potential problem: the orders of the point groups are not 24: C  One potential problem: the orders of the point groups are not
 25: C  calculated for new stationary points. Of course, they are all one for 25: C  calculated for new stationary points. Of course, they are all one for
 26: C  this system. This affects the rates as well, so we can.t just copy 26: C  this system. This affects the rates as well, so we can.t just copy
 27: C  the corresponding quantities from the original transition state. 27: C  the corresponding quantities from the original transition state.
 28: C 28: C
 29:       SUBROUTINE ADDPERM(LPOINTSTS,LPLUS,LMINUS) 29:       SUBROUTINE ADDPERM(LPOINTSTS,LPLUS,LMINUS)
 30:       USE PORFUNCS 30:       USE PORFUNCS
 31:       USE COMMONS 31:       USE COMMONS
 32:       IMPLICIT NONE 32:       IMPLICIT NONE
 33:       INTEGER J1,J2,DOTS,ISTAT 33:       INTEGER J1,J2,DOTS,ISTAT
 34:       DOUBLE PRECISION LPOINTS(NOPT), LDUMMY,NIX,NIY,NIZ,LPOINTSTS(NOPT),LPLUS(NOPT), 34:       DOUBLE PRECISION LPOINTS(3*NATOMS), LDUMMY,NIX,NIY,NIZ,LPOINTSTS(3*NATOMS),LPLUS(3*NATOMS),
 35:      1                 LMINUS(NOPT), FRICTIONFAC 35:      1                 LMINUS(3*NATOMS), FRICTIONFAC
 36:  36: 
 37:       PRINT '(A)','WARNING- routine ADDPERM has not been tested' 37:       PRINT '(A)','WARNING- routine ADDPERM has not been tested'
 38:       IF (NTAG.LE.0) THEN 38:       IF (NTAG.LE.0) THEN
 39:          WRITE(*,'(A,I5)') 'WARNING - addperm called but NTAG=',NTAG 39:          WRITE(*,'(A,I5)') 'WARNING - addperm called but NTAG=',NTAG
 40:          RETURN 40:          RETURN
 41:       ENDIF 41:       ENDIF
 42:       DOTS=NTS 42:       DOTS=NTS
 43:  43: 
 44:       DO J1=2,NATOMS 44:       DO J1=2,NATOMS
 45:          DO J2=1,NOPT 45:          DO J2=1,3*NATOMS
 46:             LPOINTS(J2)=LPOINTSTS(J2) 46:             LPOINTS(J2)=LPOINTSTS(J2)
 47:          ENDDO 47:          ENDDO
 48: C 48: C
 49: C  Permute tagged atom NTAG and atom J1 49: C  Permute tagged atom NTAG and atom J1
 50: C 50: C
 51:          LDUMMY=LPOINTS(3*(NTAG-1)+1) 51:          LDUMMY=LPOINTS(3*(NTAG-1)+1)
 52:          LPOINTS(3*(NTAG-1)+1)=LPOINTS(3*(J1-1)+1) 52:          LPOINTS(3*(NTAG-1)+1)=LPOINTS(3*(J1-1)+1)
 53:          LPOINTS(3*(J1-1)+1)=LDUMMY 53:          LPOINTS(3*(J1-1)+1)=LDUMMY
 54:          LDUMMY=LPOINTS(3*(NTAG-1)+2) 54:          LDUMMY=LPOINTS(3*(NTAG-1)+2)
 55:          LPOINTS(3*(NTAG-1)+2)=LPOINTS(3*(J1-1)+2) 55:          LPOINTS(3*(NTAG-1)+2)=LPOINTS(3*(J1-1)+2)
 56:          LPOINTS(3*(J1-1)+2)=LDUMMY 56:          LPOINTS(3*(J1-1)+2)=LDUMMY
 57:          LDUMMY=LPOINTS(3*(NTAG-1)+3) 57:          LDUMMY=LPOINTS(3*(NTAG-1)+3)
 58:          LPOINTS(3*(NTAG-1)+3)=LPOINTS(3*(J1-1)+3) 58:          LPOINTS(3*(NTAG-1)+3)=LPOINTS(3*(J1-1)+3)
 59:          LPOINTS(3*(J1-1)+3)=LDUMMY 59:          LPOINTS(3*(J1-1)+3)=LDUMMY
 60:  60: 
 61:          CALL INERTIAWRAPPER(LPOINTS,NOPT,angleAxis,NIX,NIY,NIZ) 61:          CALL INERTIAWRAPPER(LPOINTS,NATOMS,angleAxis,NIX,NIY,NIZ)
 62:          DO J2=1,NTS 62:          DO J2=1,NTS
 63:             IF ((ABS(ETS(J2)-ETS(DOTS)).LT.EDIFFTOL).AND. 63:             IF ((ABS(ETS(J2)-ETS(DOTS)).LT.EDIFFTOL).AND.
 64:      1          (ABS(IXTS(J2)-NIX).LT.IDIFFTOL).AND. 64:      1          (ABS(IXTS(J2)-NIX).LT.IDIFFTOL).AND.
 65:      2          (ABS(IYTS(J2)-NIY).LT.IDIFFTOL).AND. 65:      2          (ABS(IYTS(J2)-NIY).LT.IDIFFTOL).AND.
 66:      3          (ABS(IZTS(J2)-NIZ).LT.IDIFFTOL)) THEN 66:      3          (ABS(IZTS(J2)-NIZ).LT.IDIFFTOL)) THEN
 67:                WRITE(*,'(A,I5,A,I5,A,I5)')  67:                WRITE(*,'(A,I5,A,I5,A,I5)') 
 68:      1                 'permuting atoms 1 and ',J1,' for ts ',DOTS,' gives old transition state ',J2 68:      1                 'permuting atoms 1 and ',J1,' for ts ',DOTS,' gives old transition state ',J2
 69:                GOTO 10 69:                GOTO 10
 70:             ENDIF 70:             ENDIF
 71:          ENDDO 71:          ENDDO
 76:          FVIBTS(NTS)=FVIBTS(DOTS) 76:          FVIBTS(NTS)=FVIBTS(DOTS)
 77:          IF (FRICTIONT) NEGEIG(NTS)=NEGEIG(DOTS) 77:          IF (FRICTIONT) NEGEIG(NTS)=NEGEIG(DOTS)
 78:          HORDERTS(NTS)=1          ! not valid in general ! see next line 78:          HORDERTS(NTS)=1          ! not valid in general ! see next line
 79:          IF ((NATOMS.NE.7).OR.(.NOT.TWOD)) THEN 79:          IF ((NATOMS.NE.7).OR.(.NOT.TWOD)) THEN
 80:             PRINT*,'addperm assumes all minima have point group order 1 - quit' 80:             PRINT*,'addperm assumes all minima have point group order 1 - quit'
 81:             STOP 81:             STOP
 82:          ENDIF 82:          ENDIF
 83:          IXTS(NTS)=NIX 83:          IXTS(NTS)=NIX
 84:          IYTS(NTS)=NIY 84:          IYTS(NTS)=NIY
 85:          IZTS(NTS)=NIZ 85:          IZTS(NTS)=NIZ
 86:          WRITE(UTS,REC=NTS) (LPOINTS(J2),J2=1,NOPT) 86:          WRITE(UTS,REC=NTS) (LPOINTS(J2),J2=1,3*NATOMS)
 87:          call flush(UTS) 87:          call flush(UTS)
 88: C 88: C
 89: C  Identify the connected minimum. 89: C  Identify the connected minimum.
 90: C 90: C
 91:          DO J2=1,NOPT 91:          DO J2=1,3*NATOMS
 92:             LPOINTS(J2)=LPLUS(J2) 92:             LPOINTS(J2)=LPLUS(J2)
 93:          ENDDO 93:          ENDDO
 94:          LDUMMY=LPOINTS(3*(NTAG-1)+1) 94:          LDUMMY=LPOINTS(3*(NTAG-1)+1)
 95:          LPOINTS(3*(NTAG-1)+1)=LPOINTS(3*(J1-1)+1) 95:          LPOINTS(3*(NTAG-1)+1)=LPOINTS(3*(J1-1)+1)
 96:          LPOINTS(3*(J1-1)+1)=LDUMMY 96:          LPOINTS(3*(J1-1)+1)=LDUMMY
 97:          LDUMMY=LPOINTS(3*(NTAG-1)+2) 97:          LDUMMY=LPOINTS(3*(NTAG-1)+2)
 98:          LPOINTS(3*(NTAG-1)+2)=LPOINTS(3*(J1-1)+2) 98:          LPOINTS(3*(NTAG-1)+2)=LPOINTS(3*(J1-1)+2)
 99:          LPOINTS(3*(J1-1)+2)=LDUMMY 99:          LPOINTS(3*(J1-1)+2)=LDUMMY
100:          LDUMMY=LPOINTS(3*(NTAG-1)+3)100:          LDUMMY=LPOINTS(3*(NTAG-1)+3)
101:          LPOINTS(3*(NTAG-1)+3)=LPOINTS(3*(J1-1)+3)101:          LPOINTS(3*(NTAG-1)+3)=LPOINTS(3*(J1-1)+3)
102:          LPOINTS(3*(J1-1)+3)=LDUMMY102:          LPOINTS(3*(J1-1)+3)=LDUMMY
103: 103: 
104:          CALL INERTIAWRAPPER(LPOINTS,NOPT,angleAxis,NIX,NIY,NIZ)104:          CALL INERTIAWRAPPER(LPOINTS,NATOMS,angleAxis,NIX,NIY,NIZ)
105:          DO J2=1,NMIN105:          DO J2=1,NMIN
106:             IF ((ABS(EMIN(PLUS(DOTS))-EMIN(J2)).LT.EDIFFTOL).AND.106:             IF ((ABS(EMIN(PLUS(DOTS))-EMIN(J2)).LT.EDIFFTOL).AND.
107:      1          (ABS(NIX-IXMIN(J2)).LT.IDIFFTOL).AND.107:      1          (ABS(NIX-IXMIN(J2)).LT.IDIFFTOL).AND.
108:      2          (ABS(NIY-IYMIN(J2)).LT.IDIFFTOL).AND.108:      2          (ABS(NIY-IYMIN(J2)).LT.IDIFFTOL).AND.
109:      3          (ABS(NIZ-IZMIN(J2)).LT.IDIFFTOL)) THEN109:      3          (ABS(NIZ-IZMIN(J2)).LT.IDIFFTOL)) THEN
110:                WRITE(*,'(A,I5,A,I5)') 'permuting atoms 1 and ',J1,' for plus minimum gives old minimum ',J2110:                WRITE(*,'(A,I5,A,I5)') 'permuting atoms 1 and ',J1,' for plus minimum gives old minimum ',J2
111:                PLUS(NTS)=J2111:                PLUS(NTS)=J2
112:                IF (ENSEMBLE.EQ.'T') THEN112:                IF (ENSEMBLE.EQ.'T') THEN
113:                   KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(J2)  / (2.0D0 * PI*HORDERTS(NTS))) +113:                   KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(J2)  / (2.0D0 * PI*HORDERTS(NTS))) +
114:      1                       (FVIBMIN(J2)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J2))/TEMPERATURE114:      1                       (FVIBMIN(J2)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J2))/TEMPERATURE
145:             PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))145:             PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
146:          ELSEIF (ENSEMBLE.EQ.'E') THEN146:          ELSEIF (ENSEMBLE.EQ.'E') THEN
147:             IF (TOTALE.GT.EMIN(NMIN)) THEN147:             IF (TOTALE.GT.EMIN(NMIN)) THEN
148:                PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))148:                PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
149:             ELSE149:             ELSE
150:                PFMIN(NMIN) = -1.0D250150:                PFMIN(NMIN) = -1.0D250
151:             ENDIF151:             ENDIF
152:          ENDIF152:          ENDIF
153:          PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN153:          PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN
154: 154: 
155:          WRITE(UMIN,REC=NMIN) (LPOINTS(J2),J2=1,NOPT)155:          WRITE(UMIN,REC=NMIN) (LPOINTS(J2),J2=1,3*NATOMS)
156:          CALL FLUSH(UMIN)156:          CALL FLUSH(UMIN)
157:          IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')157:          IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
158:          WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)158:          WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)
159:          CALL FLUSH(UMINDATA)159:          CALL FLUSH(UMINDATA)
160:          IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)160:          IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
161: 161: 
162:          PLUS(NTS)=NMIN162:          PLUS(NTS)=NMIN
163:          IF (ENSEMBLE.EQ.'T') THEN163:          IF (ENSEMBLE.EQ.'T') THEN
164:             KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN)  / (2.0D0 * PI*HORDERTS(NTS))) +164:             KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN)  / (2.0D0 * PI*HORDERTS(NTS))) +
165:      1                 (FVIBMIN(NMIN)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE165:      1                 (FVIBMIN(NMIN)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE
171:             ELSE171:             ELSE
172:                KPLUS(NTS)=-1.0D250172:                KPLUS(NTS)=-1.0D250
173:             ENDIF173:             ENDIF
174:          ENDIF174:          ENDIF
175:          IF (ZSYM(1:2).EQ.'CA') KPLUS(NTS)=KPLUS(NTS)+30.66356D0175:          IF (ZSYM(1:2).EQ.'CA') KPLUS(NTS)=KPLUS(NTS)+30.66356D0
176:          IF (PLUS(NTS).EQ.MINUS(NTS)) KPLUS(NTS)=KPLUS(NTS)+LOG(2.0D0)  ! degenenerate rearrangement176:          IF (PLUS(NTS).EQ.MINUS(NTS)) KPLUS(NTS)=KPLUS(NTS)+LOG(2.0D0)  ! degenenerate rearrangement
177: !        KSUM(NMIN)=0.0D0177: !        KSUM(NMIN)=0.0D0
178: !        IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(NMIN)=KPLUS(NTS)178: !        IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(NMIN)=KPLUS(NTS)
179: 11       CONTINUE179: 11       CONTINUE
180: 180: 
181:          DO J2=1,NOPT181:          DO J2=1,3*NATOMS
182:             LPOINTS(J2)=LMINUS(J2)182:             LPOINTS(J2)=LMINUS(J2)
183:          ENDDO183:          ENDDO
184:          LDUMMY=LPOINTS(3*(NTAG-1)+1)184:          LDUMMY=LPOINTS(3*(NTAG-1)+1)
185:          LPOINTS(3*(NTAG-1)+1)=LPOINTS(3*(J1-1)+1)185:          LPOINTS(3*(NTAG-1)+1)=LPOINTS(3*(J1-1)+1)
186:          LPOINTS(3*(J1-1)+1)=LDUMMY186:          LPOINTS(3*(J1-1)+1)=LDUMMY
187:          LDUMMY=LPOINTS(3*(NTAG-1)+2)187:          LDUMMY=LPOINTS(3*(NTAG-1)+2)
188:          LPOINTS(3*(NTAG-1)+2)=LPOINTS(3*(J1-1)+2)188:          LPOINTS(3*(NTAG-1)+2)=LPOINTS(3*(J1-1)+2)
189:          LPOINTS(3*(J1-1)+2)=LDUMMY189:          LPOINTS(3*(J1-1)+2)=LDUMMY
190:          LDUMMY=LPOINTS(3*(NTAG-1)+3)190:          LDUMMY=LPOINTS(3*(NTAG-1)+3)
191:          LPOINTS(3*(NTAG-1)+3)=LPOINTS(3*(J1-1)+3)191:          LPOINTS(3*(NTAG-1)+3)=LPOINTS(3*(J1-1)+3)
192:          LPOINTS(3*(J1-1)+3)=LDUMMY192:          LPOINTS(3*(J1-1)+3)=LDUMMY
193: 193: 
194:          CALL INERTIAWRAPPER(LPOINTS,NOPT,angleAxis,NIX,NIY,NIZ)194:          CALL INERTIAWRAPPER(LPOINTS,NATOMS,angleAxis,NIX,NIY,NIZ)
195:          DO J2=1,NMIN195:          DO J2=1,NMIN
196:             IF ((ABS(EMIN(MINUS(DOTS))-EMIN(J2)).LT.EDIFFTOL).AND.196:             IF ((ABS(EMIN(MINUS(DOTS))-EMIN(J2)).LT.EDIFFTOL).AND.
197:      1          (ABS(NIX-IXMIN(J2)).LT.IDIFFTOL).AND.197:      1          (ABS(NIX-IXMIN(J2)).LT.IDIFFTOL).AND.
198:      2          (ABS(NIY-IYMIN(J2)).LT.IDIFFTOL).AND.198:      2          (ABS(NIY-IYMIN(J2)).LT.IDIFFTOL).AND.
199:      3          (ABS(NIZ-IZMIN(J2)).LT.IDIFFTOL)) THEN199:      3          (ABS(NIZ-IZMIN(J2)).LT.IDIFFTOL)) THEN
200:                WRITE(*,'(A,I5,A,I5)') 'permuting atoms 1 and ',J1,' for minus minimum gives old minimum ',J2200:                WRITE(*,'(A,I5,A,I5)') 'permuting atoms 1 and ',J1,' for minus minimum gives old minimum ',J2
201:                MINUS(NTS)=J2201:                MINUS(NTS)=J2
202:                IF (ENSEMBLE.EQ.'T') THEN202:                IF (ENSEMBLE.EQ.'T') THEN
203:                   KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(J2) / (2.0D0 * PI*HORDERTS(NTS))) +203:                   KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(J2) / (2.0D0 * PI*HORDERTS(NTS))) +
204:      1                     (FVIBMIN(J2) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J2))/TEMPERATURE204:      1                     (FVIBMIN(J2) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J2))/TEMPERATURE
236:             PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))236:             PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
237:          ELSEIF (ENSEMBLE.EQ.'E') THEN237:          ELSEIF (ENSEMBLE.EQ.'E') THEN
238:             IF (TOTALE.GT.EMIN(NMIN)) THEN238:             IF (TOTALE.GT.EMIN(NMIN)) THEN
239:                PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))239:                PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
240:             ELSE240:             ELSE
241:                PFMIN(NMIN) = -1.0D250241:                PFMIN(NMIN) = -1.0D250
242:             ENDIF242:             ENDIF
243:          ENDIF243:          ENDIF
244:          PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN244:          PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN
245: 245: 
246:          WRITE(UMIN,REC=NMIN) (LPOINTS(J2),J2=1,NOPT)246:          WRITE(UMIN,REC=NMIN) (LPOINTS(J2),J2=1,3*NATOMS)
247:          call flush(UMIN)247:          call flush(UMIN)
248:          IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')248:          IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
249:          WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN),IZMIN(NMIN)249:          WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN),IZMIN(NMIN)
250:          CALL FLUSH(UMINDATA)250:          CALL FLUSH(UMINDATA)
251:          IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)251:          IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
252:          MINUS(NTS)=NMIN252:          MINUS(NTS)=NMIN
253:          IF (ENSEMBLE.EQ.'T') THEN253:          IF (ENSEMBLE.EQ.'T') THEN
254:             KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN) / (2.0D0 * PI*HORDERTS(NTS))) +254:             KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN) / (2.0D0 * PI*HORDERTS(NTS))) +
255:      1               (FVIBMIN(NMIN) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE255:      1               (FVIBMIN(NMIN) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE
256:             IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))256:             IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
326: C  calculated for new stationary points. 326: C  calculated for new stationary points. 
327: C327: C
328: C  No tagging assumed here.328: C  No tagging assumed here.
329: C329: C
330:       SUBROUTINE ADDPERM2(LDOTS,LPOINTSTS,LPLUS,LMINUS)330:       SUBROUTINE ADDPERM2(LDOTS,LPOINTSTS,LPLUS,LMINUS)
331:       USE PORFUNCS331:       USE PORFUNCS
332:       USE COMMONS332:       USE COMMONS
333:       IMPLICIT NONE333:       IMPLICIT NONE
334:       INTEGER J1,J2,DOTS,ISTAT,J3,J4,NDUMMY,NPERM,INVERSE, NPSIZE,LDOTS334:       INTEGER J1,J2,DOTS,ISTAT,J3,J4,NDUMMY,NPERM,INVERSE, NPSIZE,LDOTS
335:       INTEGER, ALLOCATABLE :: PERM(:), PERM2(:), REALPERM(:)335:       INTEGER, ALLOCATABLE :: PERM(:), PERM2(:), REALPERM(:)
336:       DOUBLE PRECISION LPOINTS(NOPT),NIX,NIY,NIZ,LPOINTSTS(NOPT),LPLUS(NOPT),336:       DOUBLE PRECISION LPOINTS(3*NATOMS),NIX,NIY,NIZ,LPOINTSTS(3*NATOMS),LPLUS(3*NATOMS),
337:      1                 LMINUS(NOPT), FRICTIONFAC, RMAT(3,3), DISTANCE, LOCALPOINTS2(NOPT), DIST2337:      1                 LMINUS(3*NATOMS), FRICTIONFAC, RMAT(3,3), DISTANCE, LOCALPOINTS2(3*NATOMS), DIST2
338: 338: 
339:       NPSIZE=PTFINISH-PTSTART+1339:       NPSIZE=PTFINISH-PTSTART+1
340:       PRINT '(A,3I6)','PTSTART,PTFINISH,NPSIZE=',PTSTART,PTFINISH,NPSIZE340:       PRINT '(A,3I6)','PTSTART,PTFINISH,NPSIZE=',PTSTART,PTFINISH,NPSIZE
341:       ALLOCATE(REALPERM(NATOMS),PERM2(NPSIZE),PERM(NPSIZE))341:       ALLOCATE(REALPERM(NATOMS),PERM2(NPSIZE),PERM(NPSIZE))
342:       DOTS=LDOTS342:       DOTS=LDOTS
343: 343: 
344:       DO J1=1,NATOMS344:       DO J1=1,NATOMS
345:          REALPERM(J1)=J1345:          REALPERM(J1)=J1
346:       ENDDO346:       ENDDO
347:       DO J1=1,NPSIZE347:       DO J1=1,NPSIZE
395: C395: C
396: 222   CONTINUE396: 222   CONTINUE
397:       DO J1=1,NATOMS397:       DO J1=1,NATOMS
398:          LPOINTS(3*(REALPERM(J1)-1)+1)=INVERSE*LPOINTSTS(3*(J1-1)+1)398:          LPOINTS(3*(REALPERM(J1)-1)+1)=INVERSE*LPOINTSTS(3*(J1-1)+1)
399:          LPOINTS(3*(REALPERM(J1)-1)+2)=INVERSE*LPOINTSTS(3*(J1-1)+2)399:          LPOINTS(3*(REALPERM(J1)-1)+2)=INVERSE*LPOINTSTS(3*(J1-1)+2)
400:          LPOINTS(3*(REALPERM(J1)-1)+3)=INVERSE*LPOINTSTS(3*(J1-1)+3)400:          LPOINTS(3*(REALPERM(J1)-1)+3)=INVERSE*LPOINTSTS(3*(J1-1)+3)
401:       ENDDO401:       ENDDO
402: 402: 
403:       tsloop: DO J3=1,NTS403:       tsloop: DO J3=1,NTS
404:          IF (ABS(ETS(J3)-ETS(DOTS)).LT.EDIFFTOL) THEN404:          IF (ABS(ETS(J3)-ETS(DOTS)).LT.EDIFFTOL) THEN
405:             READ(UTS,REC=J3) (LOCALPOINTS2(J4),J4=1,NOPT)405:             READ(UTS,REC=J3) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
406:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,406:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
407:      &                    RMAT,.FALSE.)407:      &                    RMAT,.FALSE.)
408:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN408:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN
409:                PRINT '(A,I6)','addperm2> permutation gives database TS ',J3409:                PRINT '(A,I6)','addperm2> permutation gives database TS ',J3
410:                GOTO 130410:                GOTO 130
411:             ENDIF411:             ENDIF
412:          ENDIF412:          ENDIF
413:       ENDDO tsloop413:       ENDDO tsloop
414:       PRINT '(A,2I6,A,I6)','addperm2> permutation gives new TS ',NTS+1414:       PRINT '(A,2I6,A,I6)','addperm2> permutation gives new TS ',NTS+1
415:       NTS=NTS+1415:       NTS=NTS+1
416:       IF (NTS.GT.MAXTS) CALL TSDOUBLE416:       IF (NTS.GT.MAXTS) CALL TSDOUBLE
417:       ETS(NTS)=ETS(DOTS)417:       ETS(NTS)=ETS(DOTS)
418:       FVIBTS(NTS)=FVIBTS(DOTS)418:       FVIBTS(NTS)=FVIBTS(DOTS)
419:       IF (FRICTIONT) NEGEIG(NTS)=NEGEIG(DOTS)419:       IF (FRICTIONT) NEGEIG(NTS)=NEGEIG(DOTS)
420:       HORDERTS(NTS)=1          ! not valid in general ! see next line420:       HORDERTS(NTS)=1          ! not valid in general ! see next line
421:       CALL INERTIAWRAPPER(LPOINTS,NOPT,ANGLEAXIS,NIX,NIY,NIZ)421:       CALL INERTIAWRAPPER(LPOINTS,NATOMS,ANGLEAXIS,NIX,NIY,NIZ)
422:       IXTS(NTS)=NIX422:       IXTS(NTS)=NIX
423:       IYTS(NTS)=NIY423:       IYTS(NTS)=NIY
424:       IZTS(NTS)=NIZ424:       IZTS(NTS)=NIZ
425:       WRITE(UTS,REC=NTS) (LPOINTS(J3),J3=1,NOPT)425:       WRITE(UTS,REC=NTS) (LPOINTS(J3),J3=1,3*NATOMS)
426:       CALL FLUSH(UTS)426:       CALL FLUSH(UTS)
427: C427: C
428: C  Identify the connected minimum.428: C  Identify the connected minimum.
429: C429: C
430:       DO J3=1,NATOMS430:       DO J3=1,NATOMS
431:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LPLUS(3*(J3-1)+1)431:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LPLUS(3*(J3-1)+1)
432:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LPLUS(3*(J3-1)+2)432:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LPLUS(3*(J3-1)+2)
433:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LPLUS(3*(J3-1)+3)433:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LPLUS(3*(J3-1)+3)
434:       ENDDO434:       ENDDO
435: 435: 
436:       DO J3=1,NMIN436:       DO J3=1,NMIN
437:          IF (ABS(EMIN(J3)-EMIN(PLUS(DOTS))).LT.EDIFFTOL) THEN437:          IF (ABS(EMIN(J3)-EMIN(PLUS(DOTS))).LT.EDIFFTOL) THEN
438:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,NOPT)438:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
439:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,439:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
440:      &                 RMAT,.FALSE.)440:      &                 RMAT,.FALSE.)
441:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN441:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN
442:                WRITE(*,'(A,I6)') 'permutation for plus minimum gives old minimum ',J3442:                WRITE(*,'(A,I6)') 'permutation for plus minimum gives old minimum ',J3
443:                PLUS(NTS)=J3443:                PLUS(NTS)=J3
444:                IF (ENSEMBLE.EQ.'T') THEN444:                IF (ENSEMBLE.EQ.'T') THEN
445:                   KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(J3)  / (2.0D0 * PI*HORDERTS(NTS))) +445:                   KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(J3)  / (2.0D0 * PI*HORDERTS(NTS))) +
446:      1                 (FVIBMIN(J3)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE446:      1                 (FVIBMIN(J3)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE
447:                   IF (FRICTIONT) KPLUS(NTS)=KPLUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))447:                   IF (FRICTIONT) KPLUS(NTS)=KPLUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
448:                   ELSE448:                   ELSE
457:                GOTO 11457:                GOTO 11
458:             ENDIF458:             ENDIF
459:          ENDIF459:          ENDIF
460:       ENDDO460:       ENDDO
461:       WRITE(*,'(A,I6)') 'permutation for plus minimum gives new minimum ',NMIN+1461:       WRITE(*,'(A,I6)') 'permutation for plus minimum gives new minimum ',NMIN+1
462:       NMIN=NMIN+1462:       NMIN=NMIN+1
463:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE463:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE
464:       EMIN(NMIN)=EMIN(PLUS(DOTS))464:       EMIN(NMIN)=EMIN(PLUS(DOTS))
465:       FVIBMIN(NMIN)=FVIBMIN(PLUS(DOTS))465:       FVIBMIN(NMIN)=FVIBMIN(PLUS(DOTS))
466:       HORDERMIN(NMIN)=1          ! not valid in general !466:       HORDERMIN(NMIN)=1          ! not valid in general !
467:       CALL INERTIAWRAPPER(LPOINTS,NOPT,ANGLEAXIS,NIX,NIY,NIZ)467:       CALL INERTIAWRAPPER(LPOINTS,NATOMS,ANGLEAXIS,NIX,NIY,NIZ)
468:       IXMIN(NMIN)=NIX468:       IXMIN(NMIN)=NIX
469:       IYMIN(NMIN)=NIY469:       IYMIN(NMIN)=NIY
470:       IZMIN(NMIN)=NIZ470:       IZMIN(NMIN)=NIZ
471:       GPFOLD(NMIN)=0.0D0471:       GPFOLD(NMIN)=0.0D0
472:       IF (ENSEMBLE.EQ.'T') THEN472:       IF (ENSEMBLE.EQ.'T') THEN
473:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))473:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
474:       ELSEIF (ENSEMBLE.EQ.'E') THEN474:       ELSEIF (ENSEMBLE.EQ.'E') THEN
475:          IF (TOTALE.GT.EMIN(NMIN)) THEN475:          IF (TOTALE.GT.EMIN(NMIN)) THEN
476:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))476:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
477:          ELSE477:          ELSE
478:             PFMIN(NMIN) = -1.0D250478:             PFMIN(NMIN) = -1.0D250
479:          ENDIF479:          ENDIF
480:       ENDIF480:       ENDIF
481:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN481:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN
482: 482: 
483:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,NOPT)483:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,3*NATOMS)
484:       CALL FLUSH(UMIN)484:       CALL FLUSH(UMIN)
485:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')485:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
486:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)486:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)
487:       CALL FLUSH(UMINDATA)487:       CALL FLUSH(UMINDATA)
488:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)488:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
489: 489: 
490:       PLUS(NTS)=NMIN490:       PLUS(NTS)=NMIN
491:       IF (ENSEMBLE.EQ.'T') THEN491:       IF (ENSEMBLE.EQ.'T') THEN
492:          KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN)  / (2.0D0 * PI*HORDERTS(NTS))) +492:          KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN)  / (2.0D0 * PI*HORDERTS(NTS))) +
493:      1           (FVIBMIN(NMIN)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE493:      1           (FVIBMIN(NMIN)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE
507: 11    CONTINUE507: 11    CONTINUE
508: 508: 
509:       DO J3=1,NATOMS509:       DO J3=1,NATOMS
510:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)510:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)
511:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)511:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)
512:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)512:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)
513:       ENDDO513:       ENDDO
514: 514: 
515:       DO J3=1,NMIN515:       DO J3=1,NMIN
516:          IF (ABS(EMIN(J3)-EMIN(MINUS(DOTS))).LT.EDIFFTOL) THEN516:          IF (ABS(EMIN(J3)-EMIN(MINUS(DOTS))).LT.EDIFFTOL) THEN
517:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,NOPT)517:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
518:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,518:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
519:      &                 RMAT,.FALSE.)519:      &                 RMAT,.FALSE.)
520:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN520:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN
521:                WRITE(*,'(A,I6)') 'permutation for minus minimum gives old minimum ',J3521:                WRITE(*,'(A,I6)') 'permutation for minus minimum gives old minimum ',J3
522:                MINUS(NTS)=J3522:                MINUS(NTS)=J3
523:                IF (ENSEMBLE.EQ.'T') THEN523:                IF (ENSEMBLE.EQ.'T') THEN
524:                   KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(J3) / (2.0D0 * PI*HORDERTS(NTS))) +524:                   KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(J3) / (2.0D0 * PI*HORDERTS(NTS))) +
525:      1               (FVIBMIN(J3) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE525:      1               (FVIBMIN(J3) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE
526:                   IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))526:                   IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
527:                ELSE527:                ELSE
537:                GOTO 12537:                GOTO 12
538:             ENDIF538:             ENDIF
539:          ENDIF539:          ENDIF
540:       ENDDO540:       ENDDO
541:       WRITE(*,'(A,I5,A,I5)') 'permutation for minus minimum gives new minimum ',NMIN+1541:       WRITE(*,'(A,I5,A,I5)') 'permutation for minus minimum gives new minimum ',NMIN+1
542:       NMIN=NMIN+1542:       NMIN=NMIN+1
543:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE543:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE
544:       EMIN(NMIN)=EMIN(MINUS(DOTS))544:       EMIN(NMIN)=EMIN(MINUS(DOTS))
545:       FVIBMIN(NMIN)=FVIBMIN(MINUS(DOTS))545:       FVIBMIN(NMIN)=FVIBMIN(MINUS(DOTS))
546:       HORDERMIN(NMIN)=1          ! not valid in general !546:       HORDERMIN(NMIN)=1          ! not valid in general !
547:       CALL INERTIAWRAPPER(LPOINTS,NOPT,ANGLEAXIS,NIX,NIY,NIZ)547:       CALL INERTIAWRAPPER(LPOINTS,NATOMS,ANGLEAXIS,NIX,NIY,NIZ)
548:       IXMIN(NMIN)=NIX548:       IXMIN(NMIN)=NIX
549:       IYMIN(NMIN)=NIY549:       IYMIN(NMIN)=NIY
550:       IZMIN(NMIN)=NIZ550:       IZMIN(NMIN)=NIZ
551:       GPFOLD(NMIN)=0.0D0551:       GPFOLD(NMIN)=0.0D0
552:       IF (ENSEMBLE.EQ.'T') THEN552:       IF (ENSEMBLE.EQ.'T') THEN
553:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))553:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
554:       ELSEIF (ENSEMBLE.EQ.'E') THEN554:       ELSEIF (ENSEMBLE.EQ.'E') THEN
555:          IF (TOTALE.GT.EMIN(NMIN)) THEN555:          IF (TOTALE.GT.EMIN(NMIN)) THEN
556:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))556:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
557:          ELSE557:          ELSE
558:             PFMIN(NMIN) = -1.0D250558:             PFMIN(NMIN) = -1.0D250
559:          ENDIF559:          ENDIF
560:       ENDIF560:       ENDIF
561:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN561:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN
562: 562: 
563:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,NOPT)563:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,3*NATOMS)
564:       CALL FLUSH(UMIN)564:       CALL FLUSH(UMIN)
565:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')565:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
566:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN),IZMIN(NMIN)566:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN),IZMIN(NMIN)
567:       CALL FLUSH(UMINDATA)567:       CALL FLUSH(UMINDATA)
568:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)568:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
569:       MINUS(NTS)=NMIN569:       MINUS(NTS)=NMIN
570:       IF (ENSEMBLE.EQ.'T') THEN570:       IF (ENSEMBLE.EQ.'T') THEN
571:          KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN) / (2.0D0 * PI*HORDERTS(NTS))) +571:          KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN) / (2.0D0 * PI*HORDERTS(NTS))) +
572:      1         (FVIBMIN(NMIN) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE572:      1         (FVIBMIN(NMIN) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE
573:          IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))573:          IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
650: C650: C
651:       SUBROUTINE ADDPERM3(LDOTS,LPOINTSTS,LPLUS,LMINUS)651:       SUBROUTINE ADDPERM3(LDOTS,LPOINTSTS,LPLUS,LMINUS)
652:       USE PORFUNCS652:       USE PORFUNCS
653:       USE COMMONS653:       USE COMMONS
654:       USE UTILS,ONLY : GETUNIT654:       USE UTILS,ONLY : GETUNIT
655:       IMPLICIT NONE655:       IMPLICIT NONE
656:       INTEGER J1,LDOTS,ISTAT,J3,J4,NDUMMY,NPERMTOTAL,INVERSE,K1,NDUMMY2,NPSIZEMAX,K2,DOTS656:       INTEGER J1,LDOTS,ISTAT,J3,J4,NDUMMY,NPERMTOTAL,INVERSE,K1,NDUMMY2,NPSIZEMAX,K2,DOTS
657:       LOGICAL CHANGED, LLPERMDIST657:       LOGICAL CHANGED, LLPERMDIST
658:       INTEGER, ALLOCATABLE :: NPSIZE(:), NPERM(:), J2(:)658:       INTEGER, ALLOCATABLE :: NPSIZE(:), NPERM(:), J2(:)
659:       INTEGER, ALLOCATABLE :: PERM(:,:), PERM2(:,:), REALPERM(:)659:       INTEGER, ALLOCATABLE :: PERM(:,:), PERM2(:,:), REALPERM(:)
660:       DOUBLE PRECISION LPOINTS(NOPT),NIX,NIY,NIZ,LPOINTSTS(NOPT),LPLUS(NOPT),660:       DOUBLE PRECISION LPOINTS(3*NATOMS),NIX,NIY,NIZ,LPOINTSTS(3*NATOMS),LPLUS(3*NATOMS),
661:      1                 LMINUS(NOPT), FRICTIONFAC, RMAT(3,3), DISTANCE, LOCALPOINTS2(NOPT), DIST2661:      1                 LMINUS(3*NATOMS), FRICTIONFAC, RMAT(3,3), DISTANCE, LOCALPOINTS2(3*NATOMS), DIST2
662: 662: 
663: !663: !
664: !  Need to combine all possible permutations of each set with one another:664: !  Need to combine all possible permutations of each set with one another:
665: !  a product of factorials of their size.665: !  a product of factorials of their size.
666: !  NPERMGROUP is the number of groups666: !  NPERMGROUP is the number of groups
667: !  NPERMSIZE(J1) is the number of permutable atoms in this group667: !  NPERMSIZE(J1) is the number of permutable atoms in this group
668: !668: !
669:       ALLOCATE(NPSIZE(NPERMGROUP),J2(NPERMGROUP),NPERM(NPERMGROUP))669:       ALLOCATE(NPSIZE(NPERMGROUP),J2(NPERMGROUP),NPERM(NPERMGROUP))
670:       ALLOCATE(REALPERM(NATOMS))670:       ALLOCATE(REALPERM(NATOMS))
671:       DO J1=1,NATOMS671:       DO J1=1,NATOMS
787:       PRINT '(20I6)',INVERSE*REALPERM(1:NATOMS)787:       PRINT '(20I6)',INVERSE*REALPERM(1:NATOMS)
788:       DO J1=1,NATOMS788:       DO J1=1,NATOMS
789:          LPOINTS(3*(REALPERM(J1)-1)+1)=INVERSE*LPOINTSTS(3*(J1-1)+1)789:          LPOINTS(3*(REALPERM(J1)-1)+1)=INVERSE*LPOINTSTS(3*(J1-1)+1)
790:          LPOINTS(3*(REALPERM(J1)-1)+2)=INVERSE*LPOINTSTS(3*(J1-1)+2)790:          LPOINTS(3*(REALPERM(J1)-1)+2)=INVERSE*LPOINTSTS(3*(J1-1)+2)
791:          LPOINTS(3*(REALPERM(J1)-1)+3)=INVERSE*LPOINTSTS(3*(J1-1)+3)791:          LPOINTS(3*(REALPERM(J1)-1)+3)=INVERSE*LPOINTSTS(3*(J1-1)+3)
792:       ENDDO792:       ENDDO
793: 793: 
794:       tsloop: DO J3=1,NTS794:       tsloop: DO J3=1,NTS
795: !        PRINT '(A,I6,2F20.10)','J3,ETS(J3),ETS(DOTS)=',J3,ETS(J3),ETS(DOTS)795: !        PRINT '(A,I6,2F20.10)','J3,ETS(J3),ETS(DOTS)=',J3,ETS(J3),ETS(DOTS)
796:          IF (ABS(ETS(J3)-ETS(DOTS)).LT.EDIFFTOL) THEN796:          IF (ABS(ETS(J3)-ETS(DOTS)).LT.EDIFFTOL) THEN
797:             READ(UTS,REC=J3) (LOCALPOINTS2(J4),J4=1,NOPT)797:             READ(UTS,REC=J3) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
798:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,798:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
799:      &                    RMAT,.FALSE.)799:      &                    RMAT,.FALSE.)
800: !           PRINT '(A,2F20.10)','DISTANCE,GEOMDIFFTOL=',DISTANCE,GEOMDIFFTOL800: !           PRINT '(A,2F20.10)','DISTANCE,GEOMDIFFTOL=',DISTANCE,GEOMDIFFTOL
801:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN801:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN
802:                PRINT '(A,I6)','addperm3> permutation gives database TS ',J3802:                PRINT '(A,I6)','addperm3> permutation gives database TS ',J3
803:                GOTO 130803:                GOTO 130
804:             ENDIF804:             ENDIF
805:          ENDIF805:          ENDIF
806:       ENDDO tsloop806:       ENDDO tsloop
807:       PRINT '(A,2I6,A,I6)','addperm3> permutation gives new TS ',NTS+1807:       PRINT '(A,2I6,A,I6)','addperm3> permutation gives new TS ',NTS+1
808:       NTS=NTS+1808:       NTS=NTS+1
809:       IF (NTS.GT.MAXTS) CALL TSDOUBLE809:       IF (NTS.GT.MAXTS) CALL TSDOUBLE
810:       ETS(NTS)=ETS(DOTS)810:       ETS(NTS)=ETS(DOTS)
811:       FVIBTS(NTS)=FVIBTS(DOTS)811:       FVIBTS(NTS)=FVIBTS(DOTS)
812:       IF (FRICTIONT) NEGEIG(NTS)=NEGEIG(DOTS)812:       IF (FRICTIONT) NEGEIG(NTS)=NEGEIG(DOTS)
813:       HORDERTS(NTS)=1          ! not valid in general ! see next line813:       HORDERTS(NTS)=1          ! not valid in general ! see next line
814:       CALL INERTIAWRAPPER(LPOINTS,NOPT,ANGLEAXIS,NIX,NIY,NIZ)814:       CALL INERTIAWRAPPER(LPOINTS,NATOMS,ANGLEAXIS,NIX,NIY,NIZ)
815:       IXTS(NTS)=NIX815:       IXTS(NTS)=NIX
816:       IYTS(NTS)=NIY816:       IYTS(NTS)=NIY
817:       IZTS(NTS)=NIZ817:       IZTS(NTS)=NIZ
818:       WRITE(UTS,REC=NTS) (LPOINTS(J3),J3=1,NOPT)818:       WRITE(UTS,REC=NTS) (LPOINTS(J3),J3=1,3*NATOMS)
819:       CALL FLUSH(UTS)819:       CALL FLUSH(UTS)
820: C820: C
821: C  Identify the connected minimum.821: C  Identify the connected minimum.
822: C822: C
823:       DO J3=1,NATOMS823:       DO J3=1,NATOMS
824:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LPLUS(3*(J3-1)+1)824:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LPLUS(3*(J3-1)+1)
825:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LPLUS(3*(J3-1)+2)825:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LPLUS(3*(J3-1)+2)
826:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LPLUS(3*(J3-1)+3)826:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LPLUS(3*(J3-1)+3)
827:       ENDDO827:       ENDDO
828: 828: 
829:       DO J3=1,NMIN829:       DO J3=1,NMIN
830: !        PRINT '(A,5I10)','addperm3> J3,DOTS,NTS,PLUS(DOTS),MINUS(DOTS)=',J3,DOTS,NTS,PLUS(DOTS),MINUS(DOTS)830: !        PRINT '(A,5I10)','addperm3> J3,DOTS,NTS,PLUS(DOTS),MINUS(DOTS)=',J3,DOTS,NTS,PLUS(DOTS),MINUS(DOTS)
831:          IF (ABS(EMIN(J3)-EMIN(PLUS(DOTS))).LT.EDIFFTOL) THEN831:          IF (ABS(EMIN(J3)-EMIN(PLUS(DOTS))).LT.EDIFFTOL) THEN
832:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,NOPT)832:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
833:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,833:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
834:      &                 RMAT,.FALSE.)834:      &                 RMAT,.FALSE.)
835:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN835:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN
836:                WRITE(*,'(A,I6)') 'permutation for plus minimum gives old minimum ',J3836:                WRITE(*,'(A,I6)') 'permutation for plus minimum gives old minimum ',J3
837:                PLUS(NTS)=J3837:                PLUS(NTS)=J3
838:                IF (ENSEMBLE.EQ.'T') THEN838:                IF (ENSEMBLE.EQ.'T') THEN
839:                   KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(J3)  / (2.0D0 * PI*HORDERTS(NTS))) +839:                   KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(J3)  / (2.0D0 * PI*HORDERTS(NTS))) +
840:      1                 (FVIBMIN(J3)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE840:      1                 (FVIBMIN(J3)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE
841:                   IF (FRICTIONT) KPLUS(NTS)=KPLUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))841:                   IF (FRICTIONT) KPLUS(NTS)=KPLUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
842:                   ELSE842:                   ELSE
851:                GOTO 11851:                GOTO 11
852:             ENDIF852:             ENDIF
853:          ENDIF853:          ENDIF
854:       ENDDO854:       ENDDO
855:       WRITE(*,'(A,I6)') 'permutation for plus minimum gives new minimum ',NMIN+1855:       WRITE(*,'(A,I6)') 'permutation for plus minimum gives new minimum ',NMIN+1
856:       NMIN=NMIN+1856:       NMIN=NMIN+1
857:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE857:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE
858:       EMIN(NMIN)=EMIN(PLUS(DOTS))858:       EMIN(NMIN)=EMIN(PLUS(DOTS))
859:       FVIBMIN(NMIN)=FVIBMIN(PLUS(DOTS))859:       FVIBMIN(NMIN)=FVIBMIN(PLUS(DOTS))
860:       HORDERMIN(NMIN)=1          ! not valid in general !860:       HORDERMIN(NMIN)=1          ! not valid in general !
861:       CALL INERTIAWRAPPER(LPOINTS,NOPT,ANGLEAXIS,NIX,NIY,NIZ)861:       CALL INERTIAWRAPPER(LPOINTS,NATOMS,ANGLEAXIS,NIX,NIY,NIZ)
862:       IXMIN(NMIN)=NIX862:       IXMIN(NMIN)=NIX
863:       IYMIN(NMIN)=NIY863:       IYMIN(NMIN)=NIY
864:       IZMIN(NMIN)=NIZ864:       IZMIN(NMIN)=NIZ
865:       GPFOLD(NMIN)=0.0D0865:       GPFOLD(NMIN)=0.0D0
866:       IF (ENSEMBLE.EQ.'T') THEN866:       IF (ENSEMBLE.EQ.'T') THEN
867:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))867:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
868:       ELSEIF (ENSEMBLE.EQ.'E') THEN868:       ELSEIF (ENSEMBLE.EQ.'E') THEN
869:          IF (TOTALE.GT.EMIN(NMIN)) THEN869:          IF (TOTALE.GT.EMIN(NMIN)) THEN
870:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))870:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
871:          ELSE871:          ELSE
872:             PFMIN(NMIN) = -1.0D250872:             PFMIN(NMIN) = -1.0D250
873:          ENDIF873:          ENDIF
874:       ENDIF874:       ENDIF
875:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN875:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN
876: 876: 
877:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,NOPT)877:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,3*NATOMS)
878:       CALL FLUSH(UMIN)878:       CALL FLUSH(UMIN)
879:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')879:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
880:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)880:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)
881:       CALL FLUSH(UMINDATA)881:       CALL FLUSH(UMINDATA)
882:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)882:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
883: 883: 
884:       PLUS(NTS)=NMIN884:       PLUS(NTS)=NMIN
885:       IF (ENSEMBLE.EQ.'T') THEN885:       IF (ENSEMBLE.EQ.'T') THEN
886:          KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN)  / (2.0D0 * PI*HORDERTS(NTS))) +886:          KPLUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN)  / (2.0D0 * PI*HORDERTS(NTS))) +
887:      1           (FVIBMIN(NMIN)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE887:      1           (FVIBMIN(NMIN)  - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE
901: 11    CONTINUE901: 11    CONTINUE
902: 902: 
903:       DO J3=1,NATOMS903:       DO J3=1,NATOMS
904:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)904:          LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)
905:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)905:          LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)
906:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)906:          LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)
907:       ENDDO907:       ENDDO
908: 908: 
909:       DO J3=1,NMIN909:       DO J3=1,NMIN
910:          IF (ABS(EMIN(J3)-EMIN(MINUS(DOTS))).LT.EDIFFTOL) THEN910:          IF (ABS(EMIN(J3)-EMIN(MINUS(DOTS))).LT.EDIFFTOL) THEN
911:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,NOPT)911:             READ(UMIN,REC=J3) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
912:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,912:             CALL MINPERMDIST(LPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
913:      &                 RMAT,.FALSE.)913:      &                 RMAT,.FALSE.)
914:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN914:             IF (DISTANCE.LT.GEOMDIFFTOL) THEN
915:                WRITE(*,'(A,I6)') 'permutation for minus minimum gives old minimum ',J3915:                WRITE(*,'(A,I6)') 'permutation for minus minimum gives old minimum ',J3
916:                MINUS(NTS)=J3916:                MINUS(NTS)=J3
917:                IF (ENSEMBLE.EQ.'T') THEN917:                IF (ENSEMBLE.EQ.'T') THEN
918:                   KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(J3) / (2.0D0 * PI*HORDERTS(NTS))) +918:                   KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(J3) / (2.0D0 * PI*HORDERTS(NTS))) +
919:      1               (FVIBMIN(J3) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE919:      1               (FVIBMIN(J3) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(J3))/TEMPERATURE
920:                   IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))920:                   IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
921:                ELSE921:                ELSE
931:                GOTO 12931:                GOTO 12
932:             ENDIF932:             ENDIF
933:          ENDIF933:          ENDIF
934:       ENDDO934:       ENDDO
935:       WRITE(*,'(A,I5,A,I5)') 'permutation for minus minimum gives new minimum ',NMIN+1935:       WRITE(*,'(A,I5,A,I5)') 'permutation for minus minimum gives new minimum ',NMIN+1
936:       NMIN=NMIN+1936:       NMIN=NMIN+1
937:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE937:       IF (NMIN.GT.MAXMIN) CALL MINDOUBLE
938:       EMIN(NMIN)=EMIN(MINUS(DOTS))938:       EMIN(NMIN)=EMIN(MINUS(DOTS))
939:       FVIBMIN(NMIN)=FVIBMIN(MINUS(DOTS))939:       FVIBMIN(NMIN)=FVIBMIN(MINUS(DOTS))
940:       HORDERMIN(NMIN)=1          ! not valid in general !940:       HORDERMIN(NMIN)=1          ! not valid in general !
941:       CALL INERTIAWRAPPER(LPOINTS,NOPT,ANGLEAXIS,NIX,NIY,NIZ)941:       CALL INERTIAWRAPPER(LPOINTS,NATOMS,ANGLEAXIS,NIX,NIY,NIZ)
942:       IXMIN(NMIN)=NIX942:       IXMIN(NMIN)=NIX
943:       IYMIN(NMIN)=NIY943:       IYMIN(NMIN)=NIY
944:       IZMIN(NMIN)=NIZ944:       IZMIN(NMIN)=NIZ
945:       GPFOLD(NMIN)=0.0D0945:       GPFOLD(NMIN)=0.0D0
946:       IF (ENSEMBLE.EQ.'T') THEN946:       IF (ENSEMBLE.EQ.'T') THEN
947:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))947:          PFMIN(NMIN) = -EMIN(NMIN)/TEMPERATURE - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
948:       ELSEIF (ENSEMBLE.EQ.'E') THEN948:       ELSEIF (ENSEMBLE.EQ.'E') THEN
949:          IF (TOTALE.GT.EMIN(NMIN)) THEN949:          IF (TOTALE.GT.EMIN(NMIN)) THEN
950:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))950:             PFMIN(NMIN) = (KAPPA-1)*LOG(TOTALE-EMIN(NMIN)) - FVIBMIN(NMIN)/2.0D0 - LOG(1.0D0*HORDERMIN(NMIN))
951:          ELSE951:          ELSE
952:             PFMIN(NMIN) = -1.0D250952:             PFMIN(NMIN) = -1.0D250
953:          ENDIF953:          ENDIF
954:       ENDIF954:       ENDIF
955:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN955:       PFMIN(NMIN)=PFMIN(NMIN)-PFMEAN
956: 956: 
957:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,NOPT)957:       WRITE(UMIN,REC=NMIN) (LPOINTS(J3),J3=1,3*NATOMS)
958:       CALL FLUSH(UMIN)958:       CALL FLUSH(UMIN)
959:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')959:       IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
960:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN),IZMIN(NMIN)960:       WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN),IZMIN(NMIN)
961:       CALL FLUSH(UMINDATA)961:       CALL FLUSH(UMINDATA)
962:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)962:       IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
963:       MINUS(NTS)=NMIN963:       MINUS(NTS)=NMIN
964:       IF (ENSEMBLE.EQ.'T') THEN964:       IF (ENSEMBLE.EQ.'T') THEN
965:          KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN) / (2.0D0 * PI*HORDERTS(NTS))) +965:          KMINUS(NTS)=LOG(1.0D0 * HORDERMIN(NMIN) / (2.0D0 * PI*HORDERTS(NTS))) +
966:      1         (FVIBMIN(NMIN) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE966:      1         (FVIBMIN(NMIN) - FVIBTS(NTS)) / 2.0D0 - (ETS(NTS) - EMIN(NMIN))/TEMPERATURE
967:          IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))967:          IF (FRICTIONT) KMINUS(NTS)=KMINUS(NTS)+LOG(FRICTIONFAC(NEGEIG(NTS)))
986: !     WRITE(TSTRING,'(I10)') NTS986: !     WRITE(TSTRING,'(I10)') NTS
987: !     TSTRING='path.ts' // TRIM(ADJUSTL(TSTRING)) // '.xyz'987: !     TSTRING='path.ts' // TRIM(ADJUSTL(TSTRING)) // '.xyz'
988: 988: 
989: !     DO J3=1,NATOMS989: !     DO J3=1,NATOMS
990: !        LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)990: !        LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)
991: !        LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)991: !        LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)
992: !        LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)992: !        LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)
993: !     ENDDO993: !     ENDDO
994: !     LUNIT=GETUNIT()994: !     LUNIT=GETUNIT()
995: !     OPEN(LUNIT,FILE=TRIM(ADJUSTL(TSTRING)),STATUS='UNKNOWN')995: !     OPEN(LUNIT,FILE=TRIM(ADJUSTL(TSTRING)),STATUS='UNKNOWN')
996: !     READ(UMIN,REC=PLUS(NTS)) (LOCALPOINTS2(J4),J4=1,NOPT)996: !     READ(UMIN,REC=PLUS(NTS)) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
997: !     WRITE(LUNIT,'(I10)') NATOMS997: !     WRITE(LUNIT,'(I10)') NATOMS
998: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for nts,plus(nts)',NTS,PLUS(NTS)998: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for nts,plus(nts)',NTS,PLUS(NTS)
999: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LOCALPOINTS2(3*(J4-1)+1),LOCALPOINTS2(3*(J4-1)+2),LOCALPOINTS2(3*(J4-1)+3),J4=1,NATOMS)999: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LOCALPOINTS2(3*(J4-1)+1),LOCALPOINTS2(3*(J4-1)+2),LOCALPOINTS2(3*(J4-1)+3),J4=1,NATOMS)
1000: !     WRITE(LUNIT,'(I10)') NATOMS1000: !     WRITE(LUNIT,'(I10)') NATOMS
1001: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for permuted LPLUS:'1001: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for permuted LPLUS:'
1002: !     DO J3=1,NATOMS1002: !     DO J3=1,NATOMS
1003: !        LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LPLUS(3*(J3-1)+1)1003: !        LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LPLUS(3*(J3-1)+1)
1004: !        LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LPLUS(3*(J3-1)+2)1004: !        LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LPLUS(3*(J3-1)+2)
1005: !        LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LPLUS(3*(J3-1)+3)1005: !        LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LPLUS(3*(J3-1)+3)
1006: !     ENDDO1006: !     ENDDO
1015: !     ENDDO1015: !     ENDDO
1016: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LPOINTS(3*(J4-1)+1),LPOINTS(3*(J4-1)+2),LPOINTS(3*(J4-1)+3),J4=1,NATOMS)1016: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LPOINTS(3*(J4-1)+1),LPOINTS(3*(J4-1)+2),LPOINTS(3*(J4-1)+3),J4=1,NATOMS)
1017: !     WRITE(LUNIT,'(I10)') NATOMS1017: !     WRITE(LUNIT,'(I10)') NATOMS
1018: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for permuted LMINUS:'1018: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for permuted LMINUS:'
1019: !     DO J3=1,NATOMS1019: !     DO J3=1,NATOMS
1020: !        LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)1020: !        LPOINTS(3*(REALPERM(J3)-1)+1)=INVERSE*LMINUS(3*(J3-1)+1)
1021: !        LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)1021: !        LPOINTS(3*(REALPERM(J3)-1)+2)=INVERSE*LMINUS(3*(J3-1)+2)
1022: !        LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)1022: !        LPOINTS(3*(REALPERM(J3)-1)+3)=INVERSE*LMINUS(3*(J3-1)+3)
1023: !     ENDDO1023: !     ENDDO
1024: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LPOINTS(3*(J4-1)+1),LPOINTS(3*(J4-1)+2),LPOINTS(3*(J4-1)+3),J4=1,NATOMS)1024: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LPOINTS(3*(J4-1)+1),LPOINTS(3*(J4-1)+2),LPOINTS(3*(J4-1)+3),J4=1,NATOMS)
1025: !     READ(UMIN,REC=MINUS(NTS)) (LOCALPOINTS2(J4),J4=1,NOPT)1025: !     READ(UMIN,REC=MINUS(NTS)) (LOCALPOINTS2(J4),J4=1,3*NATOMS)
1026: !     WRITE(LUNIT,'(I10)') NATOMS1026: !     WRITE(LUNIT,'(I10)') NATOMS
1027: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for nts,minus(nts)',NTS,MINUS(NTS)1027: !     WRITE(LUNIT,'(A,2I10)') 'coordinates for nts,minus(nts)',NTS,MINUS(NTS)
1028: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LOCALPOINTS2(3*(J4-1)+1),LOCALPOINTS2(3*(J4-1)+2),LOCALPOINTS2(3*(J4-1)+3),J4=1,NATOMS)1028: !     WRITE(LUNIT,'(A2,1X,3F20.10)') ('LA ',LOCALPOINTS2(3*(J4-1)+1),LOCALPOINTS2(3*(J4-1)+2),LOCALPOINTS2(3*(J4-1)+3),J4=1,NATOMS)
1029: !1029: !
1030: !     CLOSE(LUNIT)1030: !     CLOSE(LUNIT)
1031: C1031: C
1032: C end of debug printing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1032: C end of debug printing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1033: C1033: C
1034: 1034: 
1035: C1035: C


r29929/checkspodata.f 2016-02-09 15:30:11.862674500 +0000 r29928/checkspodata.f 2016-02-09 15:30:14.146704677 +0000
 59:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN') 59:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN')
 60:          WRITE(LUNIT,'(3F20.10)') (LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2), 60:          WRITE(LUNIT,'(3F20.10)') (LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2),
 61:      &                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS) 61:      &                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS)
 62:          CLOSE(LUNIT) 62:          CLOSE(LUNIT)
 63:       ELSE IF (AMHT) THEN 63:       ELSE IF (AMHT) THEN
 64:          FPOO='start.'//TRIM(ADJUSTL(CONNSTR)) !  64:          FPOO='start.'//TRIM(ADJUSTL(CONNSTR)) ! 
 65:          LUNIT=GETUNIT() 65:          LUNIT=GETUNIT()
 66:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN') 66:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN')
 67:          WRITE(LUNIT,'(3G25.15)') LOCALPOINTS1(1:3*NATOMS) 67:          WRITE(LUNIT,'(3G25.15)') LOCALPOINTS1(1:3*NATOMS)
 68:          CLOSE(LUNIT) 68:          CLOSE(LUNIT)
 69:       ELSE IF (MLP3T.OR.MLPB3T) THEN 
 70:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) 
 71:          OPEN(2,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 
 72:          WRITE(2,'(A2,2X,F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(J2),J2=1,NOPT) 
 73:          CLOSE(2) 
 74:       ELSE 69:       ELSE
 75:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug 70:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug
 76:          LUNIT=GETUNIT() 71:          LUNIT=GETUNIT()
 77:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 72:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND')
 78:          WRITE(LUNIT,'(A2,2X,3F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2), 73:          WRITE(LUNIT,'(A2,2X,3F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2),
 79:      &                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS) 74:      &                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS)
 80:          CLOSE(LUNIT) 75:          CLOSE(LUNIT)
 81:       ENDIF 76:       ENDIF
 82:  77: 
 83:       RETURN 78:       RETURN


r29929/cycle2.f90 2016-02-09 15:30:12.106677735 +0000 r29928/cycle2.f90 2016-02-09 15:30:14.370707634 +0000
195:             IF (CONNDOING.GT.NMIN) THEN195:             IF (CONNDOING.GT.NMIN) THEN
196:                PRINT '(A)','cycle2> Completed sweep through all available minima'196:                PRINT '(A)','cycle2> Completed sweep through all available minima'
197:                NOMORE=.TRUE.197:                NOMORE=.TRUE.
198:             ENDIF198:             ENDIF
199:          ENDIF199:          ENDIF
200:       ENDIF200:       ENDIF
201:    ELSEIF (CHECKSPT) THEN201:    ELSEIF (CHECKSPT) THEN
202:       NSPCHECK=NSPCHECK+1202:       NSPCHECK=NSPCHECK+1
203:       IF (NSPCHECK.LE.CHECKSPF) THEN203:       IF (NSPCHECK.LE.CHECKSPF) THEN
204:          IF (CHECKMINT) THEN204:          IF (CHECKMINT) THEN
205:             READ(UMIN,REC=NSPCHECK) SPOINTS(1:NOPT)205:             READ(UMIN,REC=NSPCHECK) SPOINTS(1:3*NATOMS)
206:          ELSE206:          ELSE
207:             READ(UTS,REC=NSPCHECK) SPOINTS(1:NOPT)207:             READ(UTS,REC=NSPCHECK) SPOINTS(1:3*NATOMS)
208:          ENDIF208:          ENDIF
209:       ENDIF209:       ENDIF
210:    ELSE210:    ELSE
211:       CALL GETPAIR(NAVAIL,NUSED,MINS(J3),MINF(J3),SPOINTS,FPOINTS)211:       CALL GETPAIR(NAVAIL,NUSED,MINS(J3),MINF(J3),SPOINTS,FPOINTS)
212:    ENDIF212:    ENDIF
213:    IF (ADDMINXYZT) THEN 213:    IF (ADDMINXYZT) THEN 
214:       IF (ADDMINXYZCHECK.LE.NMINADDXYZ) CALL ADDMINXYZODATA(J3,SPOINTS) 214:       IF (ADDMINXYZCHECK.LE.NMINADDXYZ) CALL ADDMINXYZODATA(J3,SPOINTS) 
215:    ELSEIF (CHECKSPT) THEN 215:    ELSEIF (CHECKSPT) THEN 
216:       IF (NSPCHECK.LE.CHECKSPF) CALL CHECKSPODATA(J3,SPOINTS) 216:       IF (NSPCHECK.LE.CHECKSPF) CALL CHECKSPODATA(J3,SPOINTS) 
217:    ELSEIF (NEWCONNECTIONST) THEN 217:    ELSEIF (NEWCONNECTIONST) THEN 
468:          ADDMINXYZCHECK=ADDMINXYZCHECK+1468:          ADDMINXYZCHECK=ADDMINXYZCHECK+1
469:          IF (ADDMINXYZCHECK.LE.NMINADDXYZ) THEN469:          IF (ADDMINXYZCHECK.LE.NMINADDXYZ) THEN
470:             READ(XYZUNIT,'(I6)') NDUMMY470:             READ(XYZUNIT,'(I6)') NDUMMY
471:             READ(XYZUNIT,'(A)') SDUMMY471:             READ(XYZUNIT,'(A)') SDUMMY
472:             READ(XYZUNIT,*) (SDUMMY,SPOINTS(3*(J1-1)+1:3*(J1-1)+3),J1=1,NATOMS)472:             READ(XYZUNIT,*) (SDUMMY,SPOINTS(3*(J1-1)+1:3*(J1-1)+3),J1=1,NATOMS)
473:          ENDIF473:          ENDIF
474:       ELSEIF (CHECKSPT) THEN474:       ELSEIF (CHECKSPT) THEN
475:          NSPCHECK=NSPCHECK+1475:          NSPCHECK=NSPCHECK+1
476:          IF (NSPCHECK.LE.CHECKSPF) THEN476:          IF (NSPCHECK.LE.CHECKSPF) THEN
477:             IF (CHECKMINT) THEN477:             IF (CHECKMINT) THEN
478:                READ(UMIN,REC=NSPCHECK) SPOINTS(1:NOPT)478:                READ(UMIN,REC=NSPCHECK) SPOINTS(1:3*NATOMS)
479:             ELSE479:             ELSE
480:                READ(UTS,REC=NSPCHECK) SPOINTS(1:NOPT)480:                READ(UTS,REC=NSPCHECK) SPOINTS(1:3*NATOMS)
481:             ENDIF481:             ENDIF
482:          ENDIF482:          ENDIF
483:       ELSE483:       ELSE
484:          CALL GETPAIR(NAVAIL,NUSED,MINS(NEWJOB),MINF(NEWJOB),SPOINTS,FPOINTS)484:          CALL GETPAIR(NAVAIL,NUSED,MINS(NEWJOB),MINF(NEWJOB),SPOINTS,FPOINTS)
485:       ENDIF485:       ENDIF
486: !     WRITE(*,'(2(A,I6),A,F12.1,A,F12.3,A,I8)') 'cycle2> connecting minima ',MINS(NEWJOB),' and ', &486: !     WRITE(*,'(2(A,I6),A,F12.1,A,F12.3,A,I8)') 'cycle2> connecting minima ',MINS(NEWJOB),' and ', &
487: ! &      MINF(NEWJOB), ' distance=',DISTANCE,' |Pfold diff|=',ABS(GPFOLD(MINS(NEWJOB))-GPFOLD(MINF(NEWJOB))),' rejects=',NREJ487: ! &      MINF(NEWJOB), ' distance=',DISTANCE,' |Pfold diff|=',ABS(GPFOLD(MINS(NEWJOB))-GPFOLD(MINF(NEWJOB))),' rejects=',NREJ
488:       IF (CHECKSPT) THEN 488:       IF (CHECKSPT) THEN 
489:          IF (NSPCHECK.LE.CHECKSPF) CALL CHECKSPODATA(NEWJOB,SPOINTS) 489:          IF (NSPCHECK.LE.CHECKSPF) CALL CHECKSPODATA(NEWJOB,SPOINTS) 
490:       ELSEIF (NEWCONNECTIONST) THEN 490:       ELSEIF (NEWCONNECTIONST) THEN 


r29929/getallpaths.f 2016-02-09 15:30:12.334680741 +0000 r29928/getallpaths.f 2016-02-09 15:30:14.590710542 +0000
 23: C  This subroutine analyses a path.info in the new min-sad-min min-sad-min format 23: C  This subroutine analyses a path.info in the new min-sad-min min-sad-min format
 24: C  as generated with OPTIM keyword DUMPALLPATHS. 24: C  as generated with OPTIM keyword DUMPALLPATHS.
 25: C 25: C
 26:       SUBROUTINE GETALLPATHS 26:       SUBROUTINE GETALLPATHS
 27:       USE PORFUNCS 27:       USE PORFUNCS
 28:       USE COMMONS 28:       USE COMMONS
 29:       USE UTILS,ONLY : GETUNIT 29:       USE UTILS,ONLY : GETUNIT
 30:       IMPLICIT NONE 30:       IMPLICIT NONE
 31:  31: 
 32:       INTEGER J1, J2, ISTAT, NMINOLD, TSNUMBER, J3, NCOUNT, NMINSAVE, NTSSAVE, J4, J5, LUNIT, PLUSMIN 32:       INTEGER J1, J2, ISTAT, NMINOLD, TSNUMBER, J3, NCOUNT, NMINSAVE, NTSSAVE, J4, J5, LUNIT, PLUSMIN
 33:       DOUBLE PRECISION LOCALPOINTS(NOPT), ENERGY, NEWEMIN, NEWETS, DISTANCE, RMAT(3,3), 33:       DOUBLE PRECISION LOCALPOINTS(3*NATOMS), ENERGY, NEWEMIN, NEWETS, DISTANCE, RMAT(3,3),
 34:      1                 LPOINTSTS(NOPT), LPLUS(NOPT), LMINUS(NOPT), LOCALPOINTS2(NOPT) 34:      1                 LPOINTSTS(3*NATOMS), LPLUS(3*NATOMS), LMINUS(3*NATOMS), LOCALPOINTS2(3*NATOMS)
 35:       DOUBLE PRECISION DUMMY, DIST2, ELAPSED, TNEW 35:       DOUBLE PRECISION DUMMY, DIST2, ELAPSED, TNEW
 36:       DOUBLE PRECISION NEWFVIBMIN, NEWFVIBTS, NEWNEGEIG, NEWPOINTSMIN(NOPT), NEWPOINTSMINPLUS(NOPT), EPLUS, 36:       DOUBLE PRECISION NEWFVIBMIN, NEWFVIBTS, NEWNEGEIG, NEWPOINTSMIN(3*NATOMS), NEWPOINTSMINPLUS(3*NATOMS), EPLUS,
 37:      1                 NEWPOINTSTS(NOPT), NEWIXMIN,  NEWIYMIN, NEWIZMIN, IXPLUS, IYPLUS, IZPLUS, 37:      1                 NEWPOINTSTS(3*NATOMS), NEWIXMIN,  NEWIYMIN, NEWIZMIN, IXPLUS, IYPLUS, IZPLUS,
 38:      2                 NEWIXTS,  NEWIYTS, NEWIZTS, IXMINUS, IYMINUS, IZMINUS, FRICTIONFAC, TEMPD(PAIRDISTMAX) 38:      2                 NEWIXTS,  NEWIYTS, NEWIZTS, IXMINUS, IYMINUS, IZMINUS, FRICTIONFAC, TEMPD(PAIRDISTMAX)
 39:       INTEGER NEWHORDERMIN, NEWHORDERTS, NEWMIN, NEWTS, NTRIPLES, TEMPL(PAIRDISTMAX), NFRQS, NVARS 39:       INTEGER NEWHORDERMIN, NEWHORDERTS, NEWMIN, NEWTS, NTRIPLES, TEMPL(PAIRDISTMAX), NFRQS, NVARS
 40:       LOGICAL TSISOLD, FAILED, MINPOLD, MINMOLD, BADTRIPLE 40:       LOGICAL TSISOLD, FAILED, MINPOLD, MINMOLD, BADTRIPLE
 41:       CHARACTER(LEN=1) DUMMYSTRING 41:       CHARACTER(LEN=1) DUMMYSTRING
 42:  42: 
 43:       NMINOLD=NMIN 43:       NMINOLD=NMIN
 44:       IF (MACHINE) THEN 44:       IF (MACHINE) THEN
 45:            OPEN(1,FILE='path.info',STATUS='OLD',form='unformatted') 45:            OPEN(1,FILE='path.info',STATUS='OLD',form='unformatted')
 46:       ELSE 46:       ELSE
 47:            OPEN(1,FILE='path.info',STATUS='OLD') 47:            OPEN(1,FILE='path.info',STATUS='OLD')
 48:       ENDIF 48:       ENDIF
 49:       NFRQS=3*(NATOMS-NGLY) 49:       NFRQS=3*(NATOMS-NGLY)
 50:       NVARS=NOPT 50:       NVARS=3*NATOMS
 51:       IF (PHI4MODT) NFRQS=NATOMS 51:       IF (PHI4MODT) NFRQS=NATOMS
 52:       IF (PHI4MODT) NVARS=NATOMS 52:       IF (PHI4MODT) NVARS=NATOMS
 53:       IF (MLP3T.OR.MLPB3T) NFRQS=NATOMS 53:       IF (MLP3T.OR.MLPB3T) NFRQS=NATOMS
 54:       IF (MLP3T.OR.MLPB3T) NVARS=NATOMS 54:       IF (MLP3T.OR.MLPB3T) NVARS=NATOMS
 55: C 55: C
 56: C  Nasty things can happen if an OPTIM job crashes and leaves a path.info file incomplete. 56: C  Nasty things can happen if an OPTIM job crashes and leaves a path.info file incomplete.
 57: C  A transition state could be written without the necessary minima.  57: C  A transition state could be written without the necessary minima. 
 58: C  To avoid this problem, we now parse the path.info file first and only read in complete triples. 58: C  To avoid this problem, we now parse the path.info file first and only read in complete triples.
 59: C 59: C
 60: C  To avoid non-Morse points and higher index saddles getting into the database we now skip 60: C  To avoid non-Morse points and higher index saddles getting into the database we now skip
 76:          IF (MOD(NATOMS,3).EQ.0) THEN 76:          IF (MOD(NATOMS,3).EQ.0) THEN
 77:             NTRIPLES=NCOUNT/(3*(2*(NATOMS/3)+2)) 77:             NTRIPLES=NCOUNT/(3*(2*(NATOMS/3)+2))
 78:             J1=NTRIPLES*(3*(2*(NATOMS/3)+2)) 78:             J1=NTRIPLES*(3*(2*(NATOMS/3)+2))
 79:          ELSE 79:          ELSE
 80:             NTRIPLES=NCOUNT/(3*(2*(NATOMS/3+1)+2)) 80:             NTRIPLES=NCOUNT/(3*(2*(NATOMS/3+1)+2))
 81:             J1=NTRIPLES*(3*(2*(NATOMS/3+1)+2)) 81:             J1=NTRIPLES*(3*(2*(NATOMS/3+1)+2))
 82:          ENDIF 82:          ENDIF
 83:          IF (DEBUG) PRINT '(2(A,I8))','getallpaths> number of triples=',NTRIPLES,' number of trailing lines=',NCOUNT-J1 83:          IF (DEBUG) PRINT '(2(A,I8))','getallpaths> number of triples=',NTRIPLES,' number of trailing lines=',NCOUNT-J1
 84: ! hk286 84: ! hk286
 85:       ELSEIF ( (.NOT. NOFRQS) .AND. RIGIDINIT) THEN 85:       ELSEIF ( (.NOT. NOFRQS) .AND. RIGIDINIT) THEN
 86:          NTRIPLES=NCOUNT/(NOPT+DEGFREEDOMS+6) 86:          NTRIPLES=NCOUNT/(3*NATOMS+DEGFREEDOMS+6)
 87:          J1=NTRIPLES*(NOPT+DEGFREEDOMS+6) 87:          J1=NTRIPLES*(3*NATOMS+DEGFREEDOMS+6)
 88:          IF (DEBUG) PRINT '(2(A,I8))','getallpaths> number of triples=',NTRIPLES,' number of trailing lines=',NCOUNT-J1 88:          IF (DEBUG) PRINT '(2(A,I8))','getallpaths> number of triples=',NTRIPLES,' number of trailing lines=',NCOUNT-J1
 89:       ELSE 89:       ELSE
 90:          NTRIPLES=NCOUNT/(3*(2*NATOMS-NGLY+2)) 90:          NTRIPLES=NCOUNT/(3*(2*NATOMS-NGLY+2))
 91:          J1=NTRIPLES*3*(2*NATOMS-NGLY+2) 91:          J1=NTRIPLES*3*(2*NATOMS-NGLY+2)
 92:          IF (DEBUG) PRINT '(2(A,I8))','getallpaths> number of triples=',NTRIPLES,' number of trailing lines=',NCOUNT-J1 92:          IF (DEBUG) PRINT '(2(A,I8))','getallpaths> number of triples=',NTRIPLES,' number of trailing lines=',NCOUNT-J1
 93:       ENDIF 93:       ENDIF
 94:  94: 
 95:       TSISOLD=.TRUE. 95:       TSISOLD=.TRUE.
 96:       DO J1=1,NTRIPLES  96:       DO J1=1,NTRIPLES 
 97:          IF (DEBUG) PRINT '(A,I6,A,2I10)','getallpaths> doing triple number ',J1,' number of minima and ts=',NMIN,NTS 97:          IF (DEBUG) PRINT '(A,I6,A,2I10)','getallpaths> doing triple number ',J1,' number of minima and ts=',NMIN,NTS
116:          NEWEMIN=ENERGY116:          NEWEMIN=ENERGY
117:          IF (NOFRQS) THEN117:          IF (NOFRQS) THEN
118: !            DUMMY=4.675754133D0 ! 2 ln(2pi) +1118: !            DUMMY=4.675754133D0 ! 2 ln(2pi) +1
119:             DUMMY = 1.D0119:             DUMMY = 1.D0
120:          ELSE120:          ELSE
121:             IF (MACHINE) THEN121:             IF (MACHINE) THEN
122:                READ(1) (FRQS(J2),J2=1,NFRQS)122:                READ(1) (FRQS(J2),J2=1,NFRQS)
123:             ELSE123:             ELSE
124:                IF (RIGIDINIT) THEN ! hk286124:                IF (RIGIDINIT) THEN ! hk286
125:                   READ(1,*) (FRQS(J2),J2=1,DEGFREEDOMS)125:                   READ(1,*) (FRQS(J2),J2=1,DEGFREEDOMS)
126:                   FRQS(DEGFREEDOMS+1:NOPT) = 1.0D0126:                   FRQS(DEGFREEDOMS+1:3*NATOMS) = 1.0D0
127:                ELSE127:                ELSE
128:                   READ(1,*) (FRQS(J2),J2=1,NFRQS)128:                   READ(1,*) (FRQS(J2),J2=1,NFRQS)
129:                ENDIF129:                ENDIF
130:             ENDIF130:             ENDIF
131:             DUMMY=0.0D0131:             DUMMY=0.0D0
132:             DO J2=NFSTART,NFFINISH132:             DO J2=NFSTART,NFFINISH
133:                IF (FRQS(J2).LE.EVCUT) THEN133:                IF (FRQS(J2).LE.EVCUT) THEN
134:                   PRINT '(A,I8,A,G20.10)','getallpaths> SKIPPING - vibrational frequency ',J2,' of + minimum is ',FRQS(J2)134:                   PRINT '(A,I8,A,G20.10)','getallpaths> SKIPPING - vibrational frequency ',J2,' of + minimum is ',FRQS(J2)
135:                   BADTRIPLE=.TRUE.135:                   BADTRIPLE=.TRUE.
136:                ELSE136:                ELSE
140:          ENDIF140:          ENDIF
141:          NEWHORDERMIN=HORDER141:          NEWHORDERMIN=HORDER
142:          NEWFVIBMIN=DUMMY142:          NEWFVIBMIN=DUMMY
143: 143: 
144:          IF (MACHINE) THEN144:          IF (MACHINE) THEN
145:             READ(1) (NEWPOINTSMIN(J2),J2=1,NVARS)  145:             READ(1) (NEWPOINTSMIN(J2),J2=1,NVARS)  
146:          ELSE146:          ELSE
147:             READ(1,*) (NEWPOINTSMIN(J2),J2=1,NVARS)  147:             READ(1,*) (NEWPOINTSMIN(J2),J2=1,NVARS)  
148:          ENDIF148:          ENDIF
149:          LOCALPOINTS(1:NVARS)=NEWPOINTSMIN(1:NVARS)149:          LOCALPOINTS(1:NVARS)=NEWPOINTSMIN(1:NVARS)
150:          CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,NEWIXMIN,NEWIYMIN,NEWIZMIN)150:          CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,angleAxis,NEWIXMIN,NEWIYMIN,NEWIZMIN)
151:          MINPOLD=.TRUE.151:          MINPOLD=.TRUE.
152:          DO J2=1,NMIN152:          DO J2=1,NMIN
153:             DISTANCE=1.0D100153:             DISTANCE=1.0D100
154:             IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN154:             IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN
155:                READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NVARS)  155:                READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NVARS)  
156:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 156:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 
157:      &                          RMAT,.FALSE.)157:      &                          RMAT,.FALSE.)
158:             ENDIF158:             ENDIF
159: 159: 
160:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN160:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN
268:          ENDIF268:          ENDIF
269: C269: C
270: C  Now we store the transition state coordinates.270: C  Now we store the transition state coordinates.
271: C271: C
272:          IF (MACHINE) THEN272:          IF (MACHINE) THEN
273:             READ(1) (NEWPOINTSTS(J2),J2=1,NVARS)  273:             READ(1) (NEWPOINTSTS(J2),J2=1,NVARS)  
274:          ELSE274:          ELSE
275:             READ(1,*) (NEWPOINTSTS(J2),J2=1,NVARS)275:             READ(1,*) (NEWPOINTSTS(J2),J2=1,NVARS)
276:          ENDIF276:          ENDIF
277:          LOCALPOINTS(1:NVARS)=NEWPOINTSTS(1:NVARS)277:          LOCALPOINTS(1:NVARS)=NEWPOINTSTS(1:NVARS)
278:          CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,ANGLEAXIS,NEWIXTS,NEWIYTS,NEWIZTS)278:          CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,ANGLEAXIS,NEWIXTS,NEWIYTS,NEWIZTS)
279:          NEWFVIBTS=DUMMY279:          NEWFVIBTS=DUMMY
280:          NEWHORDERTS=HORDER280:          NEWHORDERTS=HORDER
281: 281: 
282:          IF (NEWETS.LT.NEWEMIN) PRINT '(2(A,G20.10),A)', 'getallpaths> WARNING *** New TS (energy = ',NEWETS,') is lower in energy 282:          IF (NEWETS.LT.NEWEMIN) PRINT '(2(A,G20.10),A)', 'getallpaths> WARNING *** New TS (energy = ',NEWETS,') is lower in energy 
283:      &than connected minimum (energy = ',NEWEMIN,')'283:      &than connected minimum (energy = ',NEWEMIN,')'
284: C284: C
285: C  Now check for new transition states.285: C  Now check for new transition states.
286: C286: C
287:          TSISOLD=.FALSE.287:          TSISOLD=.FALSE.
288: !        IF ((DIJINITT.OR.DIJINITFLYT).AND.(NEWETS.GT.TSTHRESH)) THEN ! reject this ts288: !        IF ((DIJINITT.OR.DIJINITFLYT).AND.(NEWETS.GT.TSTHRESH)) THEN ! reject this ts
350:          NEWEMIN=ENERGY350:          NEWEMIN=ENERGY
351:          IF (NOFRQS) THEN351:          IF (NOFRQS) THEN
352: !            DUMMY=4.675754133D0 ! 2 ln(2pi) +1352: !            DUMMY=4.675754133D0 ! 2 ln(2pi) +1
353:             DUMMY= 1.D0353:             DUMMY= 1.D0
354:          ELSE354:          ELSE
355:             IF (MACHINE) THEN355:             IF (MACHINE) THEN
356:                READ(1) (FRQS(J2),J2=1,NFRQS)356:                READ(1) (FRQS(J2),J2=1,NFRQS)
357:             ELSE357:             ELSE
358:                IF (RIGIDINIT) THEN ! hk286358:                IF (RIGIDINIT) THEN ! hk286
359:                   READ(1,*) (FRQS(J2),J2=1,DEGFREEDOMS)359:                   READ(1,*) (FRQS(J2),J2=1,DEGFREEDOMS)
360:                   FRQS(DEGFREEDOMS+1:NOPT) = 1.0D0360:                   FRQS(DEGFREEDOMS+1:3*NATOMS) = 1.0D0
361:                ELSE361:                ELSE
362:                   READ(1,*) (FRQS(J2),J2=1,NFRQS)362:                   READ(1,*) (FRQS(J2),J2=1,NFRQS)
363:                ENDIF363:                ENDIF
364:             ENDIF364:             ENDIF
365:             DUMMY=0.0D0365:             DUMMY=0.0D0
366:             DO J2=NFSTART,NFFINISH366:             DO J2=NFSTART,NFFINISH
367:                IF (FRQS(J2).LE.EVCUT) THEN367:                IF (FRQS(J2).LE.EVCUT) THEN
368:                   PRINT '(A,I8,A,G20.10)','getallpaths> SKIPPING - vibrational frequency ',J2,' of - minimum is ',FRQS(J2)368:                   PRINT '(A,I8,A,G20.10)','getallpaths> SKIPPING - vibrational frequency ',J2,' of - minimum is ',FRQS(J2)
369:                   BADTRIPLE=.TRUE.369:                   BADTRIPLE=.TRUE.
370:                ELSE370:                ELSE
376:          NEWFVIBMIN=DUMMY376:          NEWFVIBMIN=DUMMY
377:          IF (NEWETS.LT.NEWEMIN) PRINT '(2(A,G20.10),A)', 'getallpaths> WARNING *** New TS (energy = ',NEWETS,') is lower in energy 377:          IF (NEWETS.LT.NEWEMIN) PRINT '(2(A,G20.10),A)', 'getallpaths> WARNING *** New TS (energy = ',NEWETS,') is lower in energy 
378:      &than connected minimum (energy = ',NEWEMIN,')'378:      &than connected minimum (energy = ',NEWEMIN,')'
379: 379: 
380:          IF (MACHINE) THEN380:          IF (MACHINE) THEN
381:               READ(1) (NEWPOINTSMIN(J2),J2=1,NVARS)  381:               READ(1) (NEWPOINTSMIN(J2),J2=1,NVARS)  
382:          ELSE382:          ELSE
383:               READ(1,*) (NEWPOINTSMIN(J2),J2=1,NVARS)  383:               READ(1,*) (NEWPOINTSMIN(J2),J2=1,NVARS)  
384:          ENDIF384:          ENDIF
385:          LOCALPOINTS(1:NVARS)=NEWPOINTSMIN(1:NVARS)385:          LOCALPOINTS(1:NVARS)=NEWPOINTSMIN(1:NVARS)
386:          CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,NEWIXMIN,NEWIYMIN,NEWIZMIN)386:          CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,angleAxis,NEWIXMIN,NEWIYMIN,NEWIZMIN)
387:          MINMOLD=.TRUE.387:          MINMOLD=.TRUE.
388:          DO J2=1,NMIN388:          DO J2=1,NMIN
389:             DISTANCE=1.0D100389:             DISTANCE=1.0D100
390:             IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN390:             IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN
391:                READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NVARS)391:                READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NVARS)
392:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 392:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 
393:      &                          RMAT,.FALSE.)393:      &                          RMAT,.FALSE.)
394:             ENDIF394:             ENDIF
395: 395: 
396:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN396:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN
448: C  Now we know the id of the MINUS minimum we can finish the bookkeeping for this ts.448: C  Now we know the id of the MINUS minimum we can finish the bookkeeping for this ts.
449: C449: C
450:          IF (TSISOLD) THEN450:          IF (TSISOLD) THEN
451:             IF (TSNUMBER.GT.0) THEN451:             IF (TSNUMBER.GT.0) THEN
452:                FAILED=.FALSE.452:                FAILED=.FALSE.
453: C453: C
454: C  Consistency check for minima.454: C  Consistency check for minima.
455: C455: C
456:                IF (ABS(EPLUS-EMIN(PLUS(TSNUMBER))).LT.EDIFFTOL) THEN456:                IF (ABS(EPLUS-EMIN(PLUS(TSNUMBER))).LT.EDIFFTOL) THEN
457:                   IF (ABS(NEWEMIN-EMIN(MINUS(TSNUMBER))).LT.EDIFFTOL) THEN457:                   IF (ABS(NEWEMIN-EMIN(MINUS(TSNUMBER))).LT.EDIFFTOL) THEN
458:                      CALL INERTIAWRAPPER(NEWPOINTSMINPLUS,NOPT,angleAxis,IXPLUS,IYPLUS,IZPLUS)458:                      CALL INERTIAWRAPPER(NEWPOINTSMINPLUS,NATOMS,angleAxis,IXPLUS,IYPLUS,IZPLUS)
459:                      CALL INERTIAWRAPPER(NEWPOINTSMIN,NOPT,angleAxis,IXMINUS,IYMINUS,IZMINUS)459:                      CALL INERTIAWRAPPER(NEWPOINTSMIN,NATOMS,angleAxis,IXMINUS,IYMINUS,IZMINUS)
460:                      IF (DEBUG) THEN460:                      IF (DEBUG) THEN
461:                         PRINT '(A)','getallpaths> assigning + to + and - to -, energies are:'461:                         PRINT '(A)','getallpaths> assigning + to + and - to -, energies are:'
462:                         PRINT '(A,2G25.15)','getallpaths> new: ',EPLUS,NEWEMIN462:                         PRINT '(A,2G25.15)','getallpaths> new: ',EPLUS,NEWEMIN
463:                         PRINT '(A,2G25.15)','getallpaths> old: ',EMIN(PLUS(TSNUMBER)),EMIN(MINUS(TSNUMBER))463:                         PRINT '(A,2G25.15)','getallpaths> old: ',EMIN(PLUS(TSNUMBER)),EMIN(MINUS(TSNUMBER))
464:                      ENDIF464:                      ENDIF
465:                   ELSE465:                   ELSE
466:                      PRINT '(A,I6,A,G20.10)','getallpaths> WARNING failed consistency check for old ts ',466:                      PRINT '(A,I6,A,G20.10)','getallpaths> WARNING failed consistency check for old ts ',
467:      &                                        TSNUMBER,' energy=',ETS(TSNUMBER)467:      &                                        TSNUMBER,' energy=',ETS(TSNUMBER)
468:                      PRINT '(A,2G25.15)','getallpaths> +/- energies=',EPLUS,NEWEMIN468:                      PRINT '(A,2G25.15)','getallpaths> +/- energies=',EPLUS,NEWEMIN
469:                      PRINT '(A,I8,A,2G25.15)','getallpaths> ts ',TSNUMBER,469:                      PRINT '(A,I8,A,2G25.15)','getallpaths> ts ',TSNUMBER,
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:                      ENDIF497:                      ENDIF
498:                   ENDIF498:                   ENDIF
499:                ELSEIF (ABS(EPLUS-EMIN(MINUS(TSNUMBER))).LT.EDIFFTOL) THEN499:                ELSEIF (ABS(EPLUS-EMIN(MINUS(TSNUMBER))).LT.EDIFFTOL) THEN
500:                   IF (ABS(NEWEMIN-EMIN(PLUS(TSNUMBER))).LT.EDIFFTOL) THEN500:                   IF (ABS(NEWEMIN-EMIN(PLUS(TSNUMBER))).LT.EDIFFTOL) THEN
501:                      CALL INERTIAWRAPPER(NEWPOINTSMIN,NOPT,angleAxis,IXPLUS,IYPLUS,IZPLUS)501:                      CALL INERTIAWRAPPER(NEWPOINTSMIN,NATOMS,angleAxis,IXPLUS,IYPLUS,IZPLUS)
502:                      CALL INERTIAWRAPPER(NEWPOINTSMINPLUS,NOPT,angleAxis,IXMINUS,IYMINUS,IZMINUS)502:                      CALL INERTIAWRAPPER(NEWPOINTSMINPLUS,NATOMS,angleAxis,IXMINUS,IYMINUS,IZMINUS)
503:                      IF (DEBUG) THEN503:                      IF (DEBUG) THEN
504:                         PRINT '(A)','getallpaths> assigning + to - and - to +, energies are:'504:                         PRINT '(A)','getallpaths> assigning + to - and - to +, energies are:'
505:                         PRINT '(A,2G25.15)','getallpaths> new: ',EPLUS,NEWEMIN505:                         PRINT '(A,2G25.15)','getallpaths> new: ',EPLUS,NEWEMIN
506:                         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))
507:                      ENDIF507:                      ENDIF
508:                   ELSE508:                   ELSE
509:                      PRINT '(A,I6,A,G20.10)','getallpaths> WARNING failed consistency check for old ts ',509:                      PRINT '(A,I6,A,G20.10)','getallpaths> WARNING failed consistency check for old ts ',
510:      &                                        TSNUMBER,' energy=',ETS(TSNUMBER)510:      &                                        TSNUMBER,' energy=',ETS(TSNUMBER)
511:                      PRINT '(A,2G20.10,A,I8,A,2G20.10)','getallpaths> +/- energies=',EPLUS,NEWEMIN511:                      PRINT '(A,2G20.10,A,I8,A,2G20.10)','getallpaths> +/- energies=',EPLUS,NEWEMIN
512:                      PRINT '(A,I8,A,2G25.15)','getallpaths>  ts ',TSNUMBER,512:                      PRINT '(A,I8,A,2G25.15)','getallpaths>  ts ',TSNUMBER,


r29929/getnewpath.f 2016-02-09 15:30:12.550683591 +0000 r29928/getnewpath.f 2016-02-09 15:30:14.810713448 +0000
 21: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 21: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 22: C 22: C
 23: C  This subroutine analyses a new path segment from an OPTIM connect run.  23: C  This subroutine analyses a new path segment from an OPTIM connect run. 
 24: C 24: C
 25:       SUBROUTINE GETNEWPATH(MINS,MINF) 25:       SUBROUTINE GETNEWPATH(MINS,MINF)
 26:       USE PORFUNCS 26:       USE PORFUNCS
 27:       USE COMMONS 27:       USE COMMONS
 28:       IMPLICIT NONE 28:       IMPLICIT NONE
 29:  29: 
 30:       INTEGER J1, J2, ISTAT, MINS, MINF, NMINOLD, J3 30:       INTEGER J1, J2, ISTAT, MINS, MINF, NMINOLD, J3
 31:       DOUBLE PRECISION LOCALPOINTS(NOPT), ENERGY, NEWEMIN, NEWETS, DISTANCE, RMAT(3,3), 31:       DOUBLE PRECISION LOCALPOINTS(3*NATOMS), ENERGY, NEWEMIN, NEWETS, DISTANCE, RMAT(3,3),
 32:      1                 LPOINTSTS(NOPT), LPLUS(NOPT), LMINUS(NOPT), LOCALPOINTS2(NOPT), DIST2 32:      1                 LPOINTSTS(3*NATOMS), LPLUS(3*NATOMS), LMINUS(3*NATOMS), LOCALPOINTS2(3*NATOMS), DIST2
 33:       DOUBLE PRECISION DUMMY 33:       DOUBLE PRECISION DUMMY
 34:       DOUBLE PRECISION NEWFVIBMIN, NEWFVIBTS, NEWPOINTSMIN(NOPT), NEWPOINTSMINPREV(NOPT), 34:       DOUBLE PRECISION NEWFVIBMIN, NEWFVIBTS, NEWPOINTSMIN(3*NATOMS), NEWPOINTSMINPREV(3*NATOMS),
 35:      1                 NEWPOINTSTS(NOPT), NEWIXMIN,  NEWIYMIN, NEWIZMIN, 35:      1                 NEWPOINTSTS(3*NATOMS), NEWIXMIN,  NEWIYMIN, NEWIZMIN,
 36:      2                 NEWIXTS,  NEWIYTS, NEWIZTS, NEWNEGEIG, FRICTIONFAC 36:      2                 NEWIXTS,  NEWIYTS, NEWIZTS, NEWNEGEIG, FRICTIONFAC
 37:       INTEGER NEWHORDERMIN, NEWHORDERTS, NEWMIN, NEWTS 37:       INTEGER NEWHORDERMIN, NEWHORDERTS, NEWMIN, NEWTS
 38:       LOGICAL TSISOLD 38:       LOGICAL TSISOLD
 39:  39: 
 40:       NMINOLD=NMIN 40:       NMINOLD=NMIN
 41:  41: 
 42:       if (machine) then 42:       if (machine) then
 43:            OPEN(1,FILE='path.info',STATUS='OLD',form='unformatted') 43:            OPEN(1,FILE='path.info',STATUS='OLD',form='unformatted')
 44:       else 44:       else
 45:            OPEN(1,FILE='path.info',STATUS='OLD') 45:            OPEN(1,FILE='path.info',STATUS='OLD')
 83:                WRITE(*,'(A,2I6)') 'getnewpath> ERROR, order of initial minimum does not agree with database: ', 83:                WRITE(*,'(A,2I6)') 'getnewpath> ERROR, order of initial minimum does not agree with database: ',
 84:      1                                         HORDER,HORDERMIN(MINS) 84:      1                                         HORDER,HORDERMIN(MINS)
 85:                HORDER=HORDERMIN(MINS) 85:                HORDER=HORDERMIN(MINS)
 86:                WRITE(*,'(A,I6)') 'getnewpath> Using the original value of ',HORDER 86:                WRITE(*,'(A,I6)') 'getnewpath> Using the original value of ',HORDER
 87:             ENDIF 87:             ENDIF
 88:          ENDIF 88:          ENDIF
 89:          NEWHORDERMIN=HORDER 89:          NEWHORDERMIN=HORDER
 90:          NEWFVIBMIN=DUMMY 90:          NEWFVIBMIN=DUMMY
 91:  91: 
 92:          if (machine) then 92:          if (machine) then
 93:               READ(1) (NEWPOINTSMIN(J2),J2=1,NOPT)   93:               READ(1) (NEWPOINTSMIN(J2),J2=1,3*NATOMS)  
 94:          else 94:          else
 95:               READ(1,*) (NEWPOINTSMIN(J2),J2=1,NOPT)   95:               READ(1,*) (NEWPOINTSMIN(J2),J2=1,3*NATOMS)  
 96:          endif 96:          endif
 97:          LOCALPOINTS(1:NOPT)=NEWPOINTSMIN(1:NOPT) 97:          LOCALPOINTS(1:3*NATOMS)=NEWPOINTSMIN(1:3*NATOMS)
 98:          CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,ANGLEAXIS,NEWIXMIN,NEWIYMIN,NEWIZMIN) 98:          CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,ANGLEAXIS,NEWIXMIN,NEWIYMIN,NEWIZMIN)
 99:          DO J2=1,NMIN 99:          DO J2=1,NMIN
100:             DISTANCE=1.0D100100:             DISTANCE=1.0D100
101:             IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN101:             IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN
102:                READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NOPT)102:                READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,3*NATOMS)
103:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 103:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, 
104:      &                          RMAT,.FALSE.)104:      &                          RMAT,.FALSE.)
105:             ENDIF105:             ENDIF
106: 106: 
107:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN107:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN
108:                NEWMIN=J2108:                NEWMIN=J2
109:                IF (DEBUG) PRINT '(2(A,I6))','getnewpath> path minimum ',J1,' is database minimum ',J2109:                IF (DEBUG) PRINT '(2(A,I6))','getnewpath> path minimum ',J1,' is database minimum ',J2
110:                IF (ABS(NEWFVIBMIN-FVIBMIN(J2))/FVIBMIN(J2).GT.1.0D-4) THEN110:                IF (ABS(NEWFVIBMIN-FVIBMIN(J2))/FVIBMIN(J2).GT.1.0D-4) THEN
111:                   WRITE(*,'(A,F15.5,A,F15.5)') 'getnewpath> WARNING, NEWFVIBMIN=',NEWFVIBMIN,' should be ',FVIBMIN(J2)111:                   WRITE(*,'(A,F15.5,A,F15.5)') 'getnewpath> WARNING, NEWFVIBMIN=',NEWFVIBMIN,' should be ',FVIBMIN(J2)
112:                ENDIF112:                ENDIF
137:          IXMIN(NMIN)=NEWIXMIN137:          IXMIN(NMIN)=NEWIXMIN
138:          IYMIN(NMIN)=NEWIYMIN138:          IYMIN(NMIN)=NEWIYMIN
139:          IZMIN(NMIN)=NEWIZMIN139:          IZMIN(NMIN)=NEWIZMIN
140:          GPFOLD(NMIN)=0.0D0140:          GPFOLD(NMIN)=0.0D0
141:          WRITE(*,'(A,I6,A)') 'getnewpath> new minimum ',NMIN,' writing parameters to file min.data and points to points.min'141:          WRITE(*,'(A,I6,A)') 'getnewpath> new minimum ',NMIN,' writing parameters to file min.data and points to points.min'
142:          TOPPOINTER(NMIN)=0142:          TOPPOINTER(NMIN)=0
143:          IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')143:          IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
144:          WRITE(UMINDATA,'(2F20.10,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)144:          WRITE(UMINDATA,'(2F20.10,I6,3F20.10)') EMIN(NMIN), FVIBMIN(NMIN), HORDERMIN(NMIN), IXMIN(NMIN), IYMIN(NMIN), IZMIN(NMIN)
145:          CALL FLUSH(UMINDATA)145:          CALL FLUSH(UMINDATA)
146:          IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)146:          IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
147:          WRITE(UMIN,REC=NMIN) (NEWPOINTSMIN(J2),J2=1,NOPT)147:          WRITE(UMIN,REC=NMIN) (NEWPOINTSMIN(J2),J2=1,3*NATOMS)
148: 148: 
149:          CALL FLUSH(UMIN)149:          CALL FLUSH(UMIN)
150: 130      CONTINUE150: 130      CONTINUE
151: C151: C
152: C  If TSISOLD is .FALSE. the previous ts was new. Now we know the id of the MINUS152: C  If TSISOLD is .FALSE. the previous ts was new. Now we know the id of the MINUS
153: C  minimum we can finish the bookkeeping for this ts.153: C  minimum we can finish the bookkeeping for this ts.
154: C154: C
155:          IF (.NOT.TSISOLD) THEN155:          IF (.NOT.TSISOLD) THEN
156:             MINUS(NTS)=NEWMIN156:             MINUS(NTS)=NEWMIN
157:             IF (CLOSEFILEST) OPEN(UNIT=UTSDATA,FILE='ts.data',STATUS='OLD',POSITION='APPEND')157:             IF (CLOSEFILEST) OPEN(UNIT=UTSDATA,FILE='ts.data',STATUS='OLD',POSITION='APPEND')
228: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(PLUS(NTS))=KPLUS(NTS)228: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(PLUS(NTS))=KPLUS(NTS)
229: !           ELSE229: !           ELSE
230: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(PLUS(NTS)) =LOG(EXP(KSUM(PLUS(NTS))-KMEAN) + EXP(KPLUS(NTS) -KMEAN)) + KMEAN230: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(PLUS(NTS)) =LOG(EXP(KSUM(PLUS(NTS))-KMEAN) + EXP(KPLUS(NTS) -KMEAN)) + KMEAN
231: !           ENDIF231: !           ENDIF
232: !           IF (KSUM(MINUS(NTS)).EQ.0.0D0) THEN232: !           IF (KSUM(MINUS(NTS)).EQ.0.0D0) THEN
233: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(MINUS(NTS))=KMINUS(NTS)233: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(MINUS(NTS))=KMINUS(NTS)
234: !           ELSE234: !           ELSE
235: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(MINUS(NTS))=LOG(EXP(KSUM(MINUS(NTS))-KMEAN) + EXP(KMINUS(NTS)-KMEAN)) + KMEAN235: !              IF (PLUS(NTS).NE.MINUS(NTS)) KSUM(MINUS(NTS))=LOG(EXP(KSUM(MINUS(NTS))-KMEAN) + EXP(KMINUS(NTS)-KMEAN)) + KMEAN
236: !           ENDIF236: !           ENDIF
237: 237: 
238:             DO J2=1,NOPT238:             DO J2=1,3*NATOMS
239:                LPOINTSTS(J2)=NEWPOINTSTS(J2)239:                LPOINTSTS(J2)=NEWPOINTSTS(J2)
240:                LPLUS(J2)=NEWPOINTSMINPREV(J2)240:                LPLUS(J2)=NEWPOINTSMINPREV(J2)
241:                LMINUS(J2)=NEWPOINTSMIN(J2)241:                LMINUS(J2)=NEWPOINTSMIN(J2)
242:             ENDDO242:             ENDDO
243:             IF (ADDPT) CALL ADDPERM(LPOINTSTS,LPLUS,LMINUS) 243:             IF (ADDPT) CALL ADDPERM(LPOINTSTS,LPLUS,LMINUS) 
244:             IF (ADDPT2) CALL ADDPERM2(NTS,LPOINTSTS,LPLUS,LMINUS) 244:             IF (ADDPT2) CALL ADDPERM2(NTS,LPOINTSTS,LPLUS,LMINUS) 
245:             IF (ADDPT3) CALL ADDPERM3(NTS,LPOINTSTS,LPLUS,LMINUS)245:             IF (ADDPT3) CALL ADDPERM3(NTS,LPOINTSTS,LPLUS,LMINUS)
246:          ENDIF246:          ENDIF
247:          NEWPOINTSMINPREV(1:NOPT)=NEWPOINTSMIN(1:NOPT)247:          NEWPOINTSMINPREV(1:3*NATOMS)=NEWPOINTSMIN(1:3*NATOMS)
248: C248: C
249: C  Read TS data - if the above minimum is the last one then there won;t be a TS here249: C  Read TS data - if the above minimum is the last one then there won;t be a TS here
250: C  and we must exit.250: C  and we must exit.
251: C251: C
252:          if (machine) then252:          if (machine) then
253:               READ(1,END=110) ENERGY253:               READ(1,END=110) ENERGY
254:               READ(1) HORDER254:               READ(1) HORDER
255:          else255:          else
256:               READ(1,*,END=110) ENERGY, HORDER256:               READ(1,*,END=110) ENERGY, HORDER
257:          endif257:          endif
276:                ELSE276:                ELSE
277:                   DUMMY=DUMMY+LOG(FRQS(J2))277:                   DUMMY=DUMMY+LOG(FRQS(J2))
278:                ENDIF278:                ENDIF
279:             ENDDO279:             ENDDO
280:             NEWNEGEIG = FRQS(3*(NATOMS-NGLY))280:             NEWNEGEIG = FRQS(3*(NATOMS-NGLY))
281:          ENDIF281:          ENDIF
282: C282: C
283: C  Now we store the transition state coordinates.283: C  Now we store the transition state coordinates.
284: C284: C
285:          IF (MACHINE) THEN285:          IF (MACHINE) THEN
286:             READ(1) (NEWPOINTSTS(J2),J2=1,NOPT)  286:             READ(1) (NEWPOINTSTS(J2),J2=1,3*NATOMS)  
287:          ELSE287:          ELSE
288:             READ(1,*) (NEWPOINTSTS(J2),J2=1,NOPT)  288:             READ(1,*) (NEWPOINTSTS(J2),J2=1,3*NATOMS)  
289:          ENDIF289:          ENDIF
290:          LOCALPOINTS(1:NOPT)=NEWPOINTSTS(1:NOPT)290:          LOCALPOINTS(1:3*NATOMS)=NEWPOINTSTS(1:3*NATOMS)
291:          CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,ANGLEAXIS,NEWIXTS,NEWIYTS,NEWIZTS)291:          CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,ANGLEAXIS,NEWIXTS,NEWIYTS,NEWIZTS)
292:          NEWFVIBTS=DUMMY292:          NEWFVIBTS=DUMMY
293:          NEWHORDERTS=HORDER293:          NEWHORDERTS=HORDER
294: C294: C
295: C  Now check for new transition states.295: C  Now check for new transition states.
296: C296: C
297:          TSISOLD=.FALSE. 297:          TSISOLD=.FALSE. 
298: !298: !
299: ! We should also reject if both barriers are > MAXBARRIER, but we don't know the minus energy299: ! We should also reject if both barriers are > MAXBARRIER, but we don't know the minus energy
300: ! yet. We can allow such transition states into the database, but exclude them from300: ! yet. We can allow such transition states into the database, but exclude them from
301: ! pathway analysis.301: ! pathway analysis.
302: !302: !
303:          IF (DIJINITT.AND.(NEWETS.GT.TSTHRESH)) THEN ! reject this ts303:          IF (DIJINITT.AND.(NEWETS.GT.TSTHRESH)) THEN ! reject this ts
304:             PRINT '(A,I8,A,G20.10,A)','getnewpath> Transition state ',J1,' energy ',NEWETS,' lies above the threshold'304:             PRINT '(A,I8,A,G20.10,A)','getnewpath> Transition state ',J1,' energy ',NEWETS,' lies above the threshold'
305:             TSISOLD=.TRUE.305:             TSISOLD=.TRUE.
306:             GOTO 120306:             GOTO 120
307:          ENDIF307:          ENDIF
308:          DO J2=1,NTS308:          DO J2=1,NTS
309:             DISTANCE=1.0D100309:             DISTANCE=1.0D100
310:             IF (ABS(NEWETS-ETS(J2)).LT.EDIFFTOL) THEN310:             IF (ABS(NEWETS-ETS(J2)).LT.EDIFFTOL) THEN
311:                READ(UTS,REC=J2) (LOCALPOINTS2(J3),J3=1,NOPT)311:                READ(UTS,REC=J2) (LOCALPOINTS2(J3),J3=1,3*NATOMS)
312:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,312:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
313:      &                          RMAT,.FALSE.)313:      &                          RMAT,.FALSE.)
314:             ENDIF314:             ENDIF
315:             IF ((ABS(NEWETS-ETS(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN315:             IF ((ABS(NEWETS-ETS(J2)).LT.EDIFFTOL).AND.(DISTANCE.LT.GEOMDIFFTOL)) THEN
316:                IF (DEBUG) PRINT '(2(A,I6))','getnewpath> path ts ',J1,' is database ts ',J2316:                IF (DEBUG) PRINT '(2(A,I6))','getnewpath> path ts ',J1,' is database ts ',J2
317:                IF (ABS(NEWFVIBTS-FVIBTS(J2))/FVIBTS(J2).GT.1.0D-4) THEN317:                IF (ABS(NEWFVIBTS-FVIBTS(J2))/FVIBTS(J2).GT.1.0D-4) THEN
318:                   WRITE(*,'(A,F15.5,A,F15.5)') 'getnewpath> WARNING, NEWFVIBTS=',NEWFVIBTS,' should be ',FVIBTS(J2)318:                   WRITE(*,'(A,F15.5,A,F15.5)') 'getnewpath> WARNING, NEWFVIBTS=',NEWFVIBTS,' should be ',FVIBTS(J2)
319:                ENDIF319:                ENDIF
320:                IF (NEWHORDERTS.NE.HORDERTS(J2)) THEN320:                IF (NEWHORDERTS.NE.HORDERTS(J2)) THEN
321:                   WRITE(*,'(A,I6,A,I6)') 'getnewpath> ERROR, NEWHORDERTS=',NEWHORDERTS,' should be ',HORDERTS(J2)321:                   WRITE(*,'(A,I6,A,I6)') 'getnewpath> ERROR, NEWHORDERTS=',NEWHORDERTS,' should be ',HORDERTS(J2)
332:          ETS(NTS)=NEWETS332:          ETS(NTS)=NEWETS
333:          FVIBTS(NTS)=NEWFVIBTS333:          FVIBTS(NTS)=NEWFVIBTS
334:          HORDERTS(NTS)=NEWHORDERTS334:          HORDERTS(NTS)=NEWHORDERTS
335:          IXTS(NTS)=NEWIXTS335:          IXTS(NTS)=NEWIXTS
336:          IYTS(NTS)=NEWIYTS336:          IYTS(NTS)=NEWIYTS
337:          IZTS(NTS)=NEWIZTS337:          IZTS(NTS)=NEWIZTS
338:          NEGEIG(NTS)=NEWNEGEIG338:          NEGEIG(NTS)=NEWNEGEIG
339:          IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0339:          IF (DIJKSTRAT .OR. KSHORTESTPATHST) TSATTEMPT(NTS)=0
340:          WRITE(*,'(A,I6,A)') 'getnewpath> new intermediate ts ',NTS,' writing parameters to file ts.data and updating rates'340:          WRITE(*,'(A,I6,A)') 'getnewpath> new intermediate ts ',NTS,' writing parameters to file ts.data and updating rates'
341:          PLUS(NTS)=NEWMIN341:          PLUS(NTS)=NEWMIN
342:          WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,NOPT)342:          WRITE(UTS,REC=NTS) (NEWPOINTSTS(J2),J2=1,3*NATOMS)
343: 343: 
344:          CALL FLUSH(UTS)344:          CALL FLUSH(UTS)
345: 120      CONTINUE345: 120      CONTINUE
346:       ENDDO346:       ENDDO
347: 110   CLOSE(1)347: 110   CLOSE(1)
348: C348: C
349: C  We already know the parameters for the final minimum if MINF>0, so check that they agree!349: C  We already know the parameters for the final minimum if MINF>0, so check that they agree!
350: C  350: C  
351:       IF (MINF.GT.0) THEN351:       IF (MINF.GT.0) THEN
352:             DISTANCE=1.0D100352:             DISTANCE=1.0D100
353:             IF (ABS(NEWEMIN-EMIN(MINF)).LT.EDIFFTOL) THEN353:             IF (ABS(NEWEMIN-EMIN(MINF)).LT.EDIFFTOL) THEN
354:                READ(UMIN,REC=MINF) (LOCALPOINTS2(J3),J3=1,NOPT)354:                READ(UMIN,REC=MINF) (LOCALPOINTS2(J3),J3=1,3*NATOMS)
355:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,355:                CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY,
356:      &                          RMAT,.FALSE.)356:      &                          RMAT,.FALSE.)
357:             ENDIF357:             ENDIF
358:          IF ((ABS(NEWEMIN-EMIN(MINF)).GT.EDIFFTOL).OR.(DISTANCE.GT.GEOMDIFFTOL)) THEN358:          IF ((ABS(NEWEMIN-EMIN(MINF)).GT.EDIFFTOL).OR.(DISTANCE.GT.GEOMDIFFTOL)) THEN
359:             WRITE(*,'(A,2F20.10,A,2F20.10)') 'getnewpath> ERROR, energy and DISTANCE=',ENERGY,359:             WRITE(*,'(A,2F20.10,A,2F20.10)') 'getnewpath> ERROR, energy and DISTANCE=',ENERGY,
360:      1               DISTANCE,' should be ',EMIN(MINF),0.0D0360:      1               DISTANCE,' should be ',EMIN(MINF),0.0D0
361:             STOP361:             STOP
362:          ENDIF362:          ENDIF
363:          IF (ABS((DUMMY-FVIBMIN(MINF))/DUMMY).GT.1.0D-4) THEN363:          IF (ABS((DUMMY-FVIBMIN(MINF))/DUMMY).GT.1.0D-4) THEN
364:             WRITE(*,'(A,F15.5,A,F15.5)') 'getnewpath> WARNING - frequencies of final minimum ',DUMMY,364:             WRITE(*,'(A,F15.5,A,F15.5)') 'getnewpath> WARNING - frequencies of final minimum ',DUMMY,


r29929/inertia.f 2016-02-09 15:30:12.786686704 +0000 r29928/inertia.f 2016-02-09 15:30:15.022716250 +0000
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19:  19: 
 20: ! ################################################################################ 20: ! ################################################################################
 21:  21: 
 22: ! Wrapper for inertia calculation to avoid explicit coordinate transformations 22: ! Wrapper for inertia calculation to avoid explicit coordinate transformations
 23: ! where possible 23: ! where possible
 24:  24: 
 25:       subroutine inertiaWrapper(Q, nCoordSites, angleAxis, ITX, ITY, ITZ) 25:       subroutine inertiaWrapper(Q, nCoordSites, angleAxis, ITX, ITY, ITZ)
 26:  26: 
 27:       USE COMMONS, ONLY: RBAAT, NOPT, MLP3T, MLPB3T 27:       USE COMMONS, ONLY: RBAAT
 28:  28: 
 29:       use rigidBodymod 29:       use rigidBodymod
 30:  30: 
 31:       implicit none 31:       implicit none
 32:  32: 
 33: ! Subroutine arguments 33: ! Subroutine arguments
 34:  34: 
 35:       integer, intent(IN) :: nCoordSites 35:       integer, intent(IN) :: nCoordSites
 36:       real (kind=kind(0.0d0)), intent(INOUT) :: Q(NOPT) 36:       real (kind=kind(0.0d0)), intent(INOUT) :: Q(3*nCoordSites)
 37:       logical, intent(IN) :: angleAxis 37:       logical, intent(IN) :: angleAxis
 38:       real (kind=kind(0.0d0)), intent(OUT) :: ITX, ITY, ITZ 38:       real (kind=kind(0.0d0)), intent(OUT) :: ITX, ITY, ITZ
 39:  39: 
 40: ! Local variables 40: ! Local variables
 41:  41: 
 42:       integer :: numCartPoints 42:       integer :: numCartPoints
 43:       real (kind=kind(0.0d0)), allocatable :: cartCoords(:) 43:       real (kind=kind(0.0d0)), allocatable :: cartCoords(:)
 44:  44: 
 45:       IF (RBAAT) THEN 45:       IF (RBAAT) THEN
 46:          CALL RBINERTIA(Q, ITX, ITY, ITZ) 46:          CALL RBINERTIA(Q, ITX, ITY, ITZ)
 47:          RETURN 47:          RETURN
 48:       ENDIF 48:       ENDIF
 49:  49: 
 50:       if (angleAxis) then 50:       if (angleAxis) then
 51:          numCartPoints = (nCoordSites/2)*rbPotential%nPhysicalSites 51:          numCartPoints = (nCoordSites/2)*rbPotential%nPhysicalSites
 52:          allocate (cartCoords(numCartPoints*3)) 52:          allocate (cartCoords(numCartPoints*3))
 53:          call systemToCartesians(nCoordSites/2, Q, cartCoords) 53:          call systemToCartesians(nCoordSites/2, Q, cartCoords)
 54:          call inertia(cartCoords,ITX,ITY,ITZ) 54:          call inertia(cartCoords,numCartPoints,ITX,ITY,ITZ)
 55:          deallocate (cartCoords) 55:          deallocate (cartCoords)
 56:       ELSEIF (MLP3T.OR.MLPB3T) THEN 
 57:          ITX=1.0D0; ITY=1.0D0; ITZ=1.0D0 
 58:       else 56:       else
 59:          call inertia(Q,ITX,ITY,ITZ) 57:          call inertia(Q,nCoordSites,ITX,ITY,ITZ)
 60:       endif 58:       endif
 61: !     PRINT *,'ccords in inertiawrapper:' 59: !     PRINT *,'ccords in inertiawrapper:'
 62: !     PRINT '(3F20.10)',Q(1:NOPT) 60: !     PRINT '(3F20.10)',Q(1:3*nCoordSites)
 63: !     PRINT '(A,3F20.10)','moments of inertia: ',ITX,ITY,ITZ 61: !     PRINT '(A,3F20.10)','moments of inertia: ',ITX,ITY,ITZ
 64:  62: 
 65:       end 63:       end
 66:  64: 
 67: C************************************************************************* 65: C*************************************************************************
 68: C 66: C
 69: C  Build inertia tensor and get principal values. 67: C  Build inertia tensor and get principal values.
 70: C 68: C
 71:       SUBROUTINE INERTIA(Q,ITX,ITY,ITZ) 69:       SUBROUTINE INERTIA(Q,nCartPoints,ITX,ITY,ITZ)
 72:       USE COMMONS 70:       USE COMMONS
 73:       IMPLICIT NONE 71:       IMPLICIT NONE
  72:       integer, intent(IN) :: nCartPoints
 74:       INTEGER J1, J2, J3 73:       INTEGER J1, J2, J3
 75:       DOUBLE PRECISION IT(3,3), Q(NOPT), CMX, CMY, CMZ, VEC(3,3), ITX, ITY, ITZ, MASST 74:       DOUBLE PRECISION IT(3,3), Q(3*nCartPoints), CMX, CMY, CMZ, VEC(3,3), ITX, ITY, ITZ, MASST
 76:  75: 
 77:       if (size(MASS).ne.NOPT/3) then 76:       if (size(MASS).ne.nCartPoints) then
 78:          print *, 'inertia> Size of MASS not equal to number of points' 77:          print *, 'inertia> Size of MASS not equal to number of points'
 79:          stop 78:          stop
 80:       endif 79:       endif
 81:  80: 
 82:       CMX=0.0D0 81:       CMX=0.0D0
 83:       CMY=0.0D0 82:       CMY=0.0D0
 84:       CMZ=0.0D0 83:       CMZ=0.0D0
 85:       MASST=0.0D0 84:       MASST=0.0D0
 86:       DO J1=1,NOPT/3 85:       DO J1=1,nCartPoints
 87:          CMX=CMX+Q(3*(J1-1)+1)*MASS(J1) 86:          CMX=CMX+Q(3*(J1-1)+1)*MASS(J1)
 88:          CMY=CMY+Q(3*(J1-1)+2)*MASS(J1) 87:          CMY=CMY+Q(3*(J1-1)+2)*MASS(J1)
 89:          CMZ=CMZ+Q(3*(J1-1)+3)*MASS(J1) 88:          CMZ=CMZ+Q(3*(J1-1)+3)*MASS(J1)
 90:          MASST=MASST+MASS(J1) 89:          MASST=MASST+MASS(J1)
 91: C        WRITE(*,'(A,I5,4F20.10)') 'I,ATMASS(I),=',J1,MASS(J1),Q(3*(J1-1)+1),Q(3*(J1-1)+2),Q(3*(J1-1)+3) 90: C        WRITE(*,'(A,I5,4F20.10)') 'I,ATMASS(I),=',J1,MASS(J1),Q(3*(J1-1)+1),Q(3*(J1-1)+2),Q(3*(J1-1)+3)
 92:       ENDDO 91:       ENDDO
 93:       CMX=CMX/MASST 92:       CMX=CMX/MASST
 94:       CMY=CMY/MASST 93:       CMY=CMY/MASST
 95:       CMZ=CMZ/MASST 94:       CMZ=CMZ/MASST
 96: C     PRINT*,'MASST,CMX,CMY,CMZ=',MASST,CMX,CMY,CMZ 95: C     PRINT*,'MASST,CMX,CMY,CMZ=',MASST,CMX,CMY,CMZ
 97:       DO J1=1,NOPT/3 96:       DO J1=1,nCartPoints
 98:          Q(3*(J1-1)+1)=Q(3*(J1-1)+1)-CMX 97:          Q(3*(J1-1)+1)=Q(3*(J1-1)+1)-CMX
 99:          Q(3*(J1-1)+2)=Q(3*(J1-1)+2)-CMY 98:          Q(3*(J1-1)+2)=Q(3*(J1-1)+2)-CMY
100:          Q(3*(J1-1)+3)=Q(3*(J1-1)+3)-CMZ 99:          Q(3*(J1-1)+3)=Q(3*(J1-1)+3)-CMZ
101:       ENDDO100:       ENDDO
102: 101: 
103:       DO J1=1,3102:       DO J1=1,3
104:          DO J2=1,3103:          DO J2=1,3
105:             IT(J1,J2)=0.0D0104:             IT(J1,J2)=0.0D0
106:             DO J3=1,NOPT/3105:             DO J3=1,nCartPoints
107:                IT(J1,J2)=IT(J1,J2)-Q(3*(J3-1)+J1)*Q(3*(J3-1)+J2)*MASS(J3)106:                IT(J1,J2)=IT(J1,J2)-Q(3*(J3-1)+J1)*Q(3*(J3-1)+J2)*MASS(J3)
108:             ENDDO107:             ENDDO
109:             IF (J1.EQ.J2) THEN108:             IF (J1.EQ.J2) THEN
110:                DO J3=1,NOPT/3109:                DO J3=1,nCartPoints
111:                   IT(J1,J2)=IT(J1,J2)+(Q(3*(J3-1)+1)**2+Q(3*(J3-1)+2)**2+Q(3*(J3-1)+3)**2)*MASS(J3)110:                   IT(J1,J2)=IT(J1,J2)+(Q(3*(J3-1)+1)**2+Q(3*(J3-1)+2)**2+Q(3*(J3-1)+3)**2)*MASS(J3)
112:                ENDDO111:                ENDDO
113:             ENDIF112:             ENDIF
114:          ENDDO113:          ENDDO
115:       ENDDO114:       ENDDO
116: 115: 
117:       CALL EIG(IT,VEC,3,3,0)116:       CALL EIG(IT,VEC,3,3,0)
118:       ITX=IT(1,1)117:       ITX=IT(1,1)
119:       ITY=IT(2,2)118:       ITY=IT(2,2)
120:       ITZ=IT(3,3)119:       ITZ=IT(3,3)
121: !120: !
122: !  Must not move frozen atoms!121: !  Must not move frozen atoms!
123: !122: !
124:       DO J1=1,NOPT/3123:       DO J1=1,nCartPoints
125:          Q(3*(J1-1)+1)=Q(3*(J1-1)+1)+CMX124:          Q(3*(J1-1)+1)=Q(3*(J1-1)+1)+CMX
126:          Q(3*(J1-1)+2)=Q(3*(J1-1)+2)+CMY125:          Q(3*(J1-1)+2)=Q(3*(J1-1)+2)+CMY
127:          Q(3*(J1-1)+3)=Q(3*(J1-1)+3)+CMZ126:          Q(3*(J1-1)+3)=Q(3*(J1-1)+3)+CMZ
128:       ENDDO127:       ENDDO
129: 128: 
130:       RETURN129:       RETURN
131:       END130:       END


r29929/mergedb.f90 2016-02-09 15:30:13.010689668 +0000 r29928/mergedb.f90 2016-02-09 15:30:15.254719314 +0000
 21: !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 21: !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 22: ! 22: !
 23: !  This subroutine merges the database information in directory PATHNAME. 23: !  This subroutine merges the database information in directory PATHNAME.
 24: ! 24: !
 25: SUBROUTINE MERGEDB 25: SUBROUTINE MERGEDB
 26: USE COMMONS,ONLY : NATOMS, IYTS, IZTS, UTSDATA, UTS, HORDERMIN, TOPPOINTER, HORDERTS, PLUS, MINUS, GPFOLD, & 26: USE COMMONS,ONLY : NATOMS, IYTS, IZTS, UTSDATA, UTS, HORDERMIN, TOPPOINTER, HORDERTS, PLUS, MINUS, GPFOLD, &
 27:    &              MAXMIN, MAXTS, FVIBTS, EMIN, FVIBMIN, IXMIN, IYMIN, IZMIN, NEGEIG, PATHNAME, ETS, DEBUG, & 27:    &              MAXMIN, MAXTS, FVIBTS, EMIN, FVIBMIN, IXMIN, IYMIN, IZMIN, NEGEIG, PATHNAME, ETS, DEBUG, &
 28:    &              NMIN, UNRST, CHARMMT, IDIFFTOL, EDIFFTOL, UMINDATA, UMIN, NTS, IXTS, NMINA, NMINB, & 28:    &              NMIN, UNRST, CHARMMT, IDIFFTOL, EDIFFTOL, UMINDATA, UMIN, NTS, IXTS, NMINA, NMINB, &
 29:    &              LOCATIONA, LOCATIONB, ANGLEAXIS, PERMDIST, BOXLX, BOXLY, BOXLZ, GEOMDIFFTOL, TWOD, & 29:    &              LOCATIONA, LOCATIONB, ANGLEAXIS, PERMDIST, BOXLX, BOXLY, BOXLZ, GEOMDIFFTOL, TWOD, &
 30:    &              RIGIDBODY, BULKT, ZSYM, PERMISOMER, IMFRQT, CLOSEFILEST, AMHT, PAIRDISTMAX, DIJINITT, & 30:    &              RIGIDBODY, BULKT, ZSYM, PERMISOMER, IMFRQT, CLOSEFILEST, AMHT, PAIRDISTMAX, DIJINITT, &
 31:    &              DIJINITCONTT, ALLTST, ADDPT2, ADDPT3, NOPT 31:    &              DIJINITCONTT, ALLTST, ADDPT2, ADDPT3
 32: USE UTILS,ONLY : GETUNIT 32: USE UTILS,ONLY : GETUNIT
 33:  33: 
 34: USE PORFUNCS 34: USE PORFUNCS
 35: IMPLICIT NONE 35: IMPLICIT NONE
 36:  36: 
 37: INTEGER J1, J2, ISTAT, NMINOLD, NTSOLD, NMINDB, NDUMMY, J3, LUNIT, J4, LPLUNIT, LPDUNIT, PLUNIT, PDUNIT 37: INTEGER J1, J2, ISTAT, NMINOLD, NTSOLD, NMINDB, NDUMMY, J3, LUNIT, J4, LPLUNIT, LPDUNIT, PLUNIT, PDUNIT
 38: DOUBLE PRECISION LOCALPOINTS(NOPT), NEWEMIN, NEWETS 38: DOUBLE PRECISION LOCALPOINTS(3*NATOMS), NEWEMIN, NEWETS
 39: DOUBLE PRECISION NEWFVIBMIN, NEWFVIBTS, NEWPOINTSMIN(NOPT), NEWNEGEIG, LPLUS(NOPT), LMINUS(NOPT), & 39: DOUBLE PRECISION NEWFVIBMIN, NEWFVIBTS, NEWPOINTSMIN(3*NATOMS), NEWNEGEIG, LPLUS(3*NATOMS), LMINUS(3*NATOMS), &
 40:   &  NEWPOINTSTS(NOPT), NEWIXMIN,  NEWIYMIN, NEWIZMIN, LPAIRDIST(PAIRDISTMAX), & 40:   &  NEWPOINTSTS(3*NATOMS), NEWIXMIN,  NEWIYMIN, NEWIZMIN, LPAIRDIST(PAIRDISTMAX), &
 41:   &  NEWIXTS,  NEWIYTS, NEWIZTS, DISTANCE, DIST2, LOCALPOINTS2(NOPT), RMAT(3,3) 41:   &  NEWIXTS,  NEWIYTS, NEWIZTS, DISTANCE, DIST2, LOCALPOINTS2(3*NATOMS), RMAT(3,3)
 42: INTEGER NEWHORDERMIN, NEWHORDERTS, NEWMIN, NEWTS, INDEX, LPAIRLIST(PAIRDISTMAX) 42: INTEGER NEWHORDERMIN, NEWHORDERTS, NEWMIN, NEWTS, INDEX, LPAIRLIST(PAIRDISTMAX)
 43: INTEGER, ALLOCATABLE :: MINMAP(:), IDUM(:), TSMAP(:), LOCATIONDB(:) 43: INTEGER, ALLOCATABLE :: MINMAP(:), IDUM(:), TSMAP(:), LOCATIONDB(:)
 44: CHARACTER(LEN=130) FNAME, PDFILE, PLFILE, LPDFILE, LPLFILE 44: CHARACTER(LEN=130) FNAME, PDFILE, PLFILE, LPDFILE, LPLFILE
 45: LOGICAL YESNO 45: LOGICAL YESNO
 46: INTEGER MINDB, TSDB 46: INTEGER MINDB, TSDB
 47: INTEGER :: NMINDBMAX=10 47: INTEGER :: NMINDBMAX=10
 48: INTEGER :: NTSDBMAX=10 48: INTEGER :: NTSDBMAX=10
 49: LOGICAL I_OPENED 49: LOGICAL I_OPENED
 50: CHARACTER(LEN=80) I_NAME, I_ACTION 50: CHARACTER(LEN=80) I_NAME, I_ACTION
 51:  51: 
 52: ! 52: !
 53: ! First check for new minima in <pathname>/min.data and <pathname>/points.min 53: ! First check for new minima in <pathname>/min.data and <pathname>/points.min
 54: ! 54: !
 55: NMINOLD=NMIN 55: NMINOLD=NMIN
 56: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/min.data' 56: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/min.data'
 57: OPEN(UNIT=1,FILE=TRIM(ADJUSTL(FNAME)),STATUS='OLD') 57: OPEN(UNIT=1,FILE=TRIM(ADJUSTL(FNAME)),STATUS='OLD')
 58: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/points.min' 58: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/points.min'
 59: OPEN(UNIT=2,FILE=TRIM(ADJUSTL(FNAME)),ACCESS='DIRECT',FORM='UNFORMATTED',STATUS='OLD',RECL=8*NOPT) 59: OPEN(UNIT=2,FILE=TRIM(ADJUSTL(FNAME)),ACCESS='DIRECT',FORM='UNFORMATTED',STATUS='OLD',RECL=8*3*NATOMS)
 60: NEWMIN=0 60: NEWMIN=0
 61: MINDB=0 61: MINDB=0
 62: ALLOCATE(MINMAP(NMINDBMAX)) 62: ALLOCATE(MINMAP(NMINDBMAX))
 63: ! 63: !
 64: ! check for pairlist and pairdist files 64: ! check for pairlist and pairdist files
 65: ! 65: !
 66: PLUNIT=-1 66: PLUNIT=-1
 67: PDUNIT=-1 67: PDUNIT=-1
 68: IF (DIJINITT.OR.DIJINITCONTT) THEN 68: IF (DIJINITT.OR.DIJINITCONTT) THEN
 69:    INQUIRE(FILE='pairlist',EXIST=YESNO) 69:    INQUIRE(FILE='pairlist',EXIST=YESNO)
121:    NEWHORDERMIN=HORDERMIN(INDEX)121:    NEWHORDERMIN=HORDERMIN(INDEX)
122:    IF (DIJINITT.OR.DIJINITCONTT) THEN122:    IF (DIJINITT.OR.DIJINITCONTT) THEN
123:       READ(PLUNIT,'(10I10)') (LPAIRLIST(J4),J4=1,PAIRDISTMAX)123:       READ(PLUNIT,'(10I10)') (LPAIRLIST(J4),J4=1,PAIRDISTMAX)
124: !     PRINT '(A)','mergedb> read pairlist values:'124: !     PRINT '(A)','mergedb> read pairlist values:'
125: !     WRITE(*,'(10I10)') (LPAIRLIST(J4),J4=1,PAIRDISTMAX)125: !     WRITE(*,'(10I10)') (LPAIRLIST(J4),J4=1,PAIRDISTMAX)
126:       READ(PDUNIT,'(10G20.10)') (LPAIRDIST(J4),J4=1,PAIRDISTMAX)126:       READ(PDUNIT,'(10G20.10)') (LPAIRDIST(J4),J4=1,PAIRDISTMAX)
127:    ENDIF127:    ENDIF
128: !128: !
129: !  Read in points and check for agreement with moments of inertia as in setup129: !  Read in points and check for agreement with moments of inertia as in setup
130: !130: !
131:    READ(2,REC=MINDB) (NEWPOINTSMIN(J2),J2=1,NOPT)  131:    READ(2,REC=MINDB) (NEWPOINTSMIN(J2),J2=1,3*NATOMS)  
132:    LOCALPOINTS(1:NOPT)=NEWPOINTSMIN(1:NOPT)132:    LOCALPOINTS(1:3*NATOMS)=NEWPOINTSMIN(1:3*NATOMS)
133:    IF (AMHT) THEN133:    IF (AMHT) THEN
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,NATOMS,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:         STOP144:         STOP
145:       ENDIF145:       ENDIF
146:    ENDIF146:    ENDIF
147: !147: !
148: !  Is it a new minimum? Set up MINMAP.148: !  Is it a new minimum? Set up MINMAP.
149: !  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 
150: !  the database we are merging! 150: !  the database we are merging! 
151: !151: !
152: !  DO J2=1,NMIN152: !  DO J2=1,NMIN
153:    DO J2=1,NMIN+NEWMIN153:    DO J2=1,NMIN+NEWMIN
154:       IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN154:       IF (ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL) THEN
155:          DISTANCE=1.0D100155:          DISTANCE=1.0D100
156:          READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,NOPT)156:          READ(UMIN,REC=J2) (LOCALPOINTS2(J3),J3=1,3*NATOMS)
157:          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, &
158:   &                       RMAT,.FALSE.)158:   &                       RMAT,.FALSE.)
159:          IF (.NOT.(ADDPT2.OR.ADDPT3)) THEN159:          IF (.NOT.(ADDPT2.OR.ADDPT3)) THEN
160:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.GT.GEOMDIFFTOL)) THEN160:             IF ((ABS(NEWEMIN-EMIN(J2)).LT.EDIFFTOL).AND.(DISTANCE.GT.GEOMDIFFTOL)) THEN
161:                PRINT '(A)',' likely error?'161:                PRINT '(A)',' likely error?'
162:                PRINT '(A,3G20.10)',' NEWEMIN,EMIN(J2),DISTANCE=',NEWEMIN,EMIN(J2),DISTANCE162:                PRINT '(A,3G20.10)',' NEWEMIN,EMIN(J2),DISTANCE=',NEWEMIN,EMIN(J2),DISTANCE
163:                PRINT '(A,2I10)',' J2,MINDB=',J2,MINDB163:                PRINT '(A,2I10)',' J2,MINDB=',J2,MINDB
164:                IF(DEBUG) THEN164:                IF(DEBUG) THEN
165:                    PRINT '(A)','LOCALPOINTS:'165:                    PRINT '(A)','LOCALPOINTS:'
166:                    PRINT '(3F20.10)',LOCALPOINTS(1:NOPT)166:                    PRINT '(3F20.10)',LOCALPOINTS(1:3*NATOMS)
167:                    PRINT '(A)','LOCALPOINTS2:'167:                    PRINT '(A)','LOCALPOINTS2:'
168:                    PRINT '(3F20.10)',LOCALPOINTS2(1:NOPT)168:                    PRINT '(3F20.10)',LOCALPOINTS2(1:3*NATOMS)
169:                ENDIF169:                ENDIF
170: !170: !
171: ! csw34> Added CHARMM structure dumping for debugging purposes171: ! csw34> Added CHARMM structure dumping for debugging purposes
172: !172: !
173:                IF (CHARMMT) THEN173:                IF (CHARMMT) THEN
174:                   PRINT '(A)','Dumping both sets of coordinates to .crd files'174:                   PRINT '(A)','Dumping both sets of coordinates to .crd files'
175:                   CALL CHARMMDUMP(LOCALPOINTS,'localpoints.crd')175:                   CALL CHARMMDUMP(LOCALPOINTS,'localpoints.crd')
176:                   CALL CHARMMDUMP(LOCALPOINTS2,'localpoints2.crd')176:                   CALL CHARMMDUMP(LOCALPOINTS2,'localpoints2.crd')
177:                ENDIF 177:                ENDIF 
178: !           STOP !!! DJW178: !           STOP !!! DJW
201:    MINMAP(MINDB)=NMIN+NEWMIN201:    MINMAP(MINDB)=NMIN+NEWMIN
202:    GPFOLD(NMIN+NEWMIN)=0.0D0202:    GPFOLD(NMIN+NEWMIN)=0.0D0
203:    TOPPOINTER(NMIN+NEWMIN)=0203:    TOPPOINTER(NMIN+NEWMIN)=0
204:    IF (DEBUG) PRINT '(2(A,I10))','mergedb> new minimum number ',NMIN+NEWMIN,' number of new minima=',NEWMIN204:    IF (DEBUG) PRINT '(2(A,I10))','mergedb> new minimum number ',NMIN+NEWMIN,' number of new minima=',NEWMIN
205:    IF (DEBUG) WRITE(*,'(A,I10,A)') 'mergedb> new minimum ',NMIN+NEWMIN,&205:    IF (DEBUG) WRITE(*,'(A,I10,A)') 'mergedb> new minimum ',NMIN+NEWMIN,&
206:   &                ' writing parameters to file min.data and points to points.min'206:   &                ' writing parameters to file min.data and points to points.min'
207:    IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')207:    IF (CLOSEFILEST) OPEN(UNIT=UMINDATA,FILE='min.data',STATUS='UNKNOWN',POSITION='APPEND')
208:    WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(INDEX),FVIBMIN(INDEX),HORDERMIN(INDEX),IXMIN(INDEX),IYMIN(INDEX),IZMIN(INDEX)208:    WRITE(UMINDATA,'(2F25.15,I6,3F20.10)') EMIN(INDEX),FVIBMIN(INDEX),HORDERMIN(INDEX),IXMIN(INDEX),IYMIN(INDEX),IZMIN(INDEX)
209:    CALL FLUSH(UMINDATA)209:    CALL FLUSH(UMINDATA)
210:    IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)210:    IF (CLOSEFILEST) CLOSE(UNIT=UMINDATA)
211:    WRITE(UMIN,REC=INDEX) (NEWPOINTSMIN(J2),J2=1,NOPT)211:    WRITE(UMIN,REC=INDEX) (NEWPOINTSMIN(J2),J2=1,3*NATOMS)
212: !212: !
213: ! pairlist and pairdist if applicable213: ! pairlist and pairdist if applicable
214: !214: !
215:    IF (DIJINITT.OR.DIJINITCONTT) THEN215:    IF (DIJINITT.OR.DIJINITCONTT) THEN
216:       WRITE(LPLUNIT,'(10I10)') (LPAIRLIST(J4),J4=1,PAIRDISTMAX)216:       WRITE(LPLUNIT,'(10I10)') (LPAIRLIST(J4),J4=1,PAIRDISTMAX)
217:       WRITE(LPDUNIT,'(10G20.10)') (LPAIRDIST(J4),J4=1,PAIRDISTMAX)217:       WRITE(LPDUNIT,'(10G20.10)') (LPAIRDIST(J4),J4=1,PAIRDISTMAX)
218:    ENDIF218:    ENDIF
219:    219:    
220: 130 CONTINUE220: 130 CONTINUE
221: ENDDO221: ENDDO
231: IF (NMIN.GT.MAXMIN) CALL MINDOUBLE231: IF (NMIN.GT.MAXMIN) CALL MINDOUBLE
232: MINDB=MINDB-1232: MINDB=MINDB-1
233: PRINT '(A,I10,A,I10,A)','mergedb> merged ',MINDB,' minima - ',NEWMIN,' are new'233: PRINT '(A,I10,A,I10,A)','mergedb> merged ',MINDB,' minima - ',NEWMIN,' are new'
234: !234: !
235: !  Now for the transition states.235: !  Now for the transition states.
236: !236: !
237: NTSOLD=NTS237: NTSOLD=NTS
238: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/ts.data'238: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/ts.data'
239: OPEN(UNIT=1,FILE=TRIM(ADJUSTL(FNAME)),STATUS='OLD')239: OPEN(UNIT=1,FILE=TRIM(ADJUSTL(FNAME)),STATUS='OLD')
240: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/points.ts'240: FNAME=TRIM(ADJUSTL(PATHNAME)) // '/points.ts'
241: OPEN(UNIT=2,FILE=TRIM(ADJUSTL(FNAME)),ACCESS='DIRECT',FORM='UNFORMATTED',STATUS='OLD',RECL=8*NOPT)241: OPEN(UNIT=2,FILE=TRIM(ADJUSTL(FNAME)),ACCESS='DIRECT',FORM='UNFORMATTED',STATUS='OLD',RECL=8*3*NATOMS)
242: NEWTS=0242: NEWTS=0
243: TSDB=0243: TSDB=0
244: NDUMMY=0244: NDUMMY=0
245: ALLOCATE(TSMAP(NTSDBMAX))245: ALLOCATE(TSMAP(NTSDBMAX))
246: DO246: DO
247:    TSDB=TSDB+1247:    TSDB=TSDB+1
248:    IF (TSDB.GT.NTSDBMAX) THEN248:    IF (TSDB.GT.NTSDBMAX) THEN
249:       ALLOCATE(IDUM(NTSDBMAX))249:       ALLOCATE(IDUM(NTSDBMAX))
250:       IDUM(1:NTSDBMAX)=TSMAP(1:NTSDBMAX)250:       IDUM(1:NTSDBMAX)=TSMAP(1:NTSDBMAX)
251:       DEALLOCATE(TSMAP)251:       DEALLOCATE(TSMAP)
262:    ELSE262:    ELSE
263:       READ(1,*,END=40) ETS(INDEX),FVIBTS(INDEX),HORDERTS(INDEX),PLUS(INDEX),MINUS(INDEX),IXTS(INDEX),IYTS(INDEX),IZTS(INDEX)263:       READ(1,*,END=40) ETS(INDEX),FVIBTS(INDEX),HORDERTS(INDEX),PLUS(INDEX),MINUS(INDEX),IXTS(INDEX),IYTS(INDEX),IZTS(INDEX)
264:    ENDIF264:    ENDIF
265:    NEWETS=ETS(INDEX)265:    NEWETS=ETS(INDEX)
266:    NEWFVIBTS=FVIBTS(INDEX)266:    NEWFVIBTS=FVIBTS(INDEX)
267:    NEWHORDERTS=HORDERTS(INDEX)267:    NEWHORDERTS=HORDERTS(INDEX)
268:    NEWNEGEIG=NEGEIG(INDEX)268:    NEWNEGEIG=NEGEIG(INDEX)
269: !269: !
270: !  Read in points and check for agreement with moments of inertia as in setup270: !  Read in points and check for agreement with moments of inertia as in setup
271: !271: !
272:    READ(2,REC=TSDB) (NEWPOINTSTS(J2),J2=1,NOPT)  272:    READ(2,REC=TSDB) (NEWPOINTSTS(J2),J2=1,3*NATOMS)  
273:    LOCALPOINTS(1:NOPT)=NEWPOINTSTS(1:NOPT)273:    LOCALPOINTS(1:3*NATOMS)=NEWPOINTSTS(1:3*NATOMS)
274:    IF (AMHT) THEN274:    IF (AMHT) THEN
275:       WRITE(*,*)'mergedb> AVOIDING MOMENT OF INITERIA CALC FOR AMH'275:       WRITE(*,*)'mergedb> AVOIDING MOMENT OF INITERIA CALC FOR AMH'
276:    ELSE276:    ELSE
277:       CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,NEWIXTS,NEWIYTS,NEWIZTS)277:       CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,angleAxis,NEWIXTS,NEWIYTS,NEWIZTS)
278:       IF ((ABS(NEWIXTS-IXTS(INDEX)).GT.IDIFFTOL).OR. &278:       IF ((ABS(NEWIXTS-IXTS(INDEX)).GT.IDIFFTOL).OR. &
279:   &       (ABS(NEWIYTS-IYTS(INDEX)).GT.IDIFFTOL).OR. &279:   &       (ABS(NEWIYTS-IYTS(INDEX)).GT.IDIFFTOL).OR. &
280:   &       (ABS(NEWIZTS-IZTS(INDEX)).GT.IDIFFTOL)) THEN280:   &       (ABS(NEWIZTS-IZTS(INDEX)).GT.IDIFFTOL)) THEN
281:          WRITE(*,'(A)') 'mergedb> possible error - principal moments of inertia do not agree with input'281:          WRITE(*,'(A)') 'mergedb> possible error - principal moments of inertia do not agree with input'
282:          WRITE(*,'(A,3F20.10)') 'mergedb> values from coordinates: ',NEWIXTS,NEWIYTS,NEWIZTS282:          WRITE(*,'(A,3F20.10)') 'mergedb> values from coordinates: ',NEWIXTS,NEWIYTS,NEWIZTS
283:          WRITE(*,'(A,3F20.10)') 'mergedb> values from ' // TRIM(ADJUSTL(PATHNAME)) // '/ts.data', &283:          WRITE(*,'(A,3F20.10)') 'mergedb> values from ' // TRIM(ADJUSTL(PATHNAME)) // '/ts.data', &
284:    &                                        IXMIN(INDEX),IYMIN(INDEX),IZMIN(INDEX)284:    &                                        IXMIN(INDEX),IYMIN(INDEX),IZMIN(INDEX)
285:          IF (.NOT.BULKT) STOP285:          IF (.NOT.BULKT) STOP
286:       ENDIF286:       ENDIF
287:    ENDIF287:    ENDIF
288: !288: !
289: !  Is it a new ts? Set up TSMAP.289: !  Is it a new ts? Set up TSMAP.
290: !290: !
291: !  DO J2=1,NTS291: !  DO J2=1,NTS
292:    DO J2=1,NTS+NEWTS292:    DO J2=1,NTS+NEWTS
293:       IF (ABS(NEWETS-ETS(J2)).LT.EDIFFTOL) THEN293:       IF (ABS(NEWETS-ETS(J2)).LT.EDIFFTOL) THEN
294:          DISTANCE=1.0D100294:          DISTANCE=1.0D100
295:          READ(UTS,REC=J2) (LOCALPOINTS2(J3),J3=1,NOPT)295:          READ(UTS,REC=J2) (LOCALPOINTS2(J3),J3=1,3*NATOMS)
296:          CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, &296:          CALL MINPERMDIST(LOCALPOINTS,LOCALPOINTS2,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGIDBODY, &
297:   &                       RMAT,.FALSE.)297:   &                       RMAT,.FALSE.)
298:          IF (DISTANCE.LT.GEOMDIFFTOL) THEN298:          IF (DISTANCE.LT.GEOMDIFFTOL) THEN
299: 299: 
300:             IF (DEBUG) PRINT '(2(A,I10))','mergedb> ts ',TSDB,' is database ts ',J2300:             IF (DEBUG) PRINT '(2(A,I10))','mergedb> ts ',TSDB,' is database ts ',J2
301:             IF (ABS(NEWFVIBTS-FVIBTS(J2))/FVIBTS(J2).GT.1.0D-3) THEN301:             IF (ABS(NEWFVIBTS-FVIBTS(J2))/FVIBTS(J2).GT.1.0D-3) THEN
302:                WRITE(*,'(A,F15.5,A,F15.5)') 'mergedb> WARNING, NEWFVIBTS=',NEWFVIBTS,' should be ',FVIBTS(J2)302:                WRITE(*,'(A,F15.5,A,F15.5)') 'mergedb> WARNING, NEWFVIBTS=',NEWFVIBTS,' should be ',FVIBTS(J2)
303:             ENDIF303:             ENDIF
304:             IF (NEWHORDERTS.NE.HORDERTS(J2)) THEN304:             IF (NEWHORDERTS.NE.HORDERTS(J2)) THEN
305:             WRITE(*,'(A,I6,A,I6)') 'mergedb> ERROR, NEWHORDERTS=',NEWHORDERTS,' should be ',HORDERTS(J2)305:             WRITE(*,'(A,I6,A,I6)') 'mergedb> ERROR, NEWHORDERTS=',NEWHORDERTS,' should be ',HORDERTS(J2)
342: 342: 
343:    IF (IMFRQT) THEN343:    IF (IMFRQT) THEN
344:       WRITE(UTSDATA,'(2F25.15,3I10,4F20.10)') ETS(INDEX),FVIBTS(INDEX),HORDERTS(INDEX),MINMAP(PLUS(INDEX)),MINMAP(MINUS(INDEX)), &344:       WRITE(UTSDATA,'(2F25.15,3I10,4F20.10)') ETS(INDEX),FVIBTS(INDEX),HORDERTS(INDEX),MINMAP(PLUS(INDEX)),MINMAP(MINUS(INDEX)), &
345:         &                                          IXTS(INDEX),IYTS(INDEX),IZTS(INDEX),NEGEIG(INDEX)345:         &                                          IXTS(INDEX),IYTS(INDEX),IZTS(INDEX),NEGEIG(INDEX)
346:    ELSE346:    ELSE
347:       WRITE(UTSDATA,'(2F25.15,3I10,3F20.10)') ETS(INDEX),FVIBTS(INDEX),HORDERTS(INDEX),MINMAP(PLUS(INDEX)),MINMAP(MINUS(INDEX)), &347:       WRITE(UTSDATA,'(2F25.15,3I10,3F20.10)') ETS(INDEX),FVIBTS(INDEX),HORDERTS(INDEX),MINMAP(PLUS(INDEX)),MINMAP(MINUS(INDEX)), &
348:         &                                          IXTS(INDEX),IYTS(INDEX),IZTS(INDEX)348:         &                                          IXTS(INDEX),IYTS(INDEX),IZTS(INDEX)
349:    ENDIF349:    ENDIF
350:    CALL FLUSH(UTSDATA)350:    CALL FLUSH(UTSDATA)
351:    IF (CLOSEFILEST) CLOSE(UNIT=UTSDATA)351:    IF (CLOSEFILEST) CLOSE(UNIT=UTSDATA)
352:    WRITE(UTS,REC=INDEX) (NEWPOINTSTS(J2),J2=1,NOPT)352:    WRITE(UTS,REC=INDEX) (NEWPOINTSTS(J2),J2=1,3*NATOMS)
353:    353:    
354: 140 CONTINUE354: 140 CONTINUE
355: ENDDO355: ENDDO
356: 40  CONTINUE356: 40  CONTINUE
357: NTS=NTSOLD+NEWTS357: NTS=NTSOLD+NEWTS
358: 358: 
359: IF (ADDPT2.OR.ADDPT3) THEN359: IF (ADDPT2.OR.ADDPT3) THEN
360:    DO J1=1,NEWTS360:    DO J1=1,NEWTS
361:       READ(UTS,REC=NTSOLD+NEWTS-NEWTS+J1) (NEWPOINTSTS(J2),J2=1,NOPT)361:       READ(UTS,REC=NTSOLD+NEWTS-NEWTS+J1) (NEWPOINTSTS(J2),J2=1,3*NATOMS)
362:       READ(UMIN,REC=PLUS(NTSOLD+NEWTS-NEWTS+J1)) (LPLUS(J2),J2=1,NOPT)362:       READ(UMIN,REC=PLUS(NTSOLD+NEWTS-NEWTS+J1)) (LPLUS(J2),J2=1,3*NATOMS)
363:       READ(UMIN,REC=MINUS(NTSOLD+NEWTS-NEWTS+J1)) (LMINUS(J2),J2=1,NOPT)363:       READ(UMIN,REC=MINUS(NTSOLD+NEWTS-NEWTS+J1)) (LMINUS(J2),J2=1,3*NATOMS)
364:       IF (ADDPT2) CALL ADDPERM2(NTSOLD+NEWTS-NEWTS+J1,NEWPOINTSTS,LPLUS,LMINUS)364:       IF (ADDPT2) CALL ADDPERM2(NTSOLD+NEWTS-NEWTS+J1,NEWPOINTSTS,LPLUS,LMINUS)
365:       IF (ADDPT3) CALL ADDPERM3(NTSOLD+NEWTS-NEWTS+J1,NEWPOINTSTS,LPLUS,LMINUS)365:       IF (ADDPT3) CALL ADDPERM3(NTSOLD+NEWTS-NEWTS+J1,NEWPOINTSTS,LPLUS,LMINUS)
366:    ENDDO366:    ENDDO
367: ENDIF367: ENDIF
368: 368: 
369: CLOSE(1)369: CLOSE(1)
370: CLOSE(2)370: CLOSE(2)
371: IF (NTS.GT.MAXTS) CALL TSDOUBLE371: IF (NTS.GT.MAXTS) CALL TSDOUBLE
372: TSDB=TSDB-1372: TSDB=TSDB-1
373: PRINT '(A,I10,A,I10,A)','mergedb> merged ',TSDB,' ts - ',NEWTS,' are new'373: PRINT '(A,I10,A,I10,A)','mergedb> merged ',TSDB,' ts - ',NEWTS,' are new'


r29929/newconnodata.f 2016-02-09 15:30:13.238692677 +0000 r29928/newconnodata.f 2016-02-09 15:30:15.458722009 +0000
 74:          LUNIT=GETUNIT() 74:          LUNIT=GETUNIT()
 75:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 75:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND')
 76:          WRITE(LUNIT,'(A2,2X,F20.10)') ('  ',LOCALPOINTS1(J2),J2=1,NOPT) 76:          WRITE(LUNIT,'(A2,2X,F20.10)') ('  ',LOCALPOINTS1(J2),J2=1,NOPT)
 77:          CLOSE(LUNIT) 77:          CLOSE(LUNIT)
 78:       ELSE IF (PHI4MODT) THEN 78:       ELSE IF (PHI4MODT) THEN
 79:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug 79:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug
 80:          LUNIT=GETUNIT() 80:          LUNIT=GETUNIT()
 81:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 81:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND')
 82:          WRITE(LUNIT,'(A2,2X,F20.10)') ('  ',LOCALPOINTS1(J2),J2=1,NATOMS) 82:          WRITE(LUNIT,'(A2,2X,F20.10)') ('  ',LOCALPOINTS1(J2),J2=1,NATOMS)
 83:          CLOSE(LUNIT) 83:          CLOSE(LUNIT)
 84:       ELSE IF (MLP3T.OR.MLPB3T) THEN 
 85:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) 
 86:          OPEN(2,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 
 87:          WRITE(2,'(A2,2X,F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(J2),J2=1,NOPT) 
 88:          CLOSE(2) 
 89:       ELSE 84:       ELSE
 90:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug 85:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug
 91:          LUNIT=GETUNIT() 86:          LUNIT=GETUNIT()
 92:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 87:          OPEN(LUNIT,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND')
 93:          WRITE(LUNIT,'(A2,2X,3F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2), 88:          WRITE(LUNIT,'(A2,2X,3F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2),
 94:      &                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS) 89:      &                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS)
 95:          CLOSE(LUNIT) 90:          CLOSE(LUNIT)
 96:       ENDIF 91:       ENDIF
 97:  92: 
 98:       RETURN 93:       RETURN


r29929/orderodata.f 2016-02-09 15:30:13.446695427 +0000 r29928/orderodata.f 2016-02-09 15:30:15.686725019 +0000
 46:                  SAVEPOINTS(J2)=LOCALPOINTS1(J2) 46:                  SAVEPOINTS(J2)=LOCALPOINTS1(J2)
 47:               ENDDO 47:               ENDDO
 48:               CALL CHARMMDUMP(SAVEPOINTS,'input.crd.'//TRIM(ADJUSTL(CONNSTR))) 48:               CALL CHARMMDUMP(SAVEPOINTS,'input.crd.'//TRIM(ADJUSTL(CONNSTR)))
 49:          endif 49:          endif
 50:       ELSE IF (UNRST) THEN 50:       ELSE IF (UNRST) THEN
 51:          DO J2=1,3*NATOMS 51:          DO J2=1,3*NATOMS
 52:             SAVEPOINTS(J2)=LOCALPOINTS1(J2) 52:             SAVEPOINTS(J2)=LOCALPOINTS1(J2)
 53:          ENDDO 53:          ENDDO
 54:          WRITE(UNSTRING,'(A)') 'coords.'//TRIM(ADJUSTL(CONNSTR)) 54:          WRITE(UNSTRING,'(A)') 'coords.'//TRIM(ADJUSTL(CONNSTR))
 55:          CALL MYUNRESDUMP(SAVEPOINTS,UNSTRING) 55:          CALL MYUNRESDUMP(SAVEPOINTS,UNSTRING)
 56:       ELSE IF (MLP3T.OR.MLPB3T) THEN 
 57:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) 
 58:          OPEN(2,FILE=TRIM(ADJUSTL(FPOO)),STATUS='UNKNOWN',POSITION='APPEND') 
 59:          WRITE(2,'(A2,2X,F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(J2),J2=1,NOPT) 
 60:          CLOSE(2) 
 61:       ELSE 56:       ELSE
 62:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug 57:          FPOO='odata.'//TRIM(ADJUSTL(CONNSTR)) ! workaround for Sun compiler bug
 63:          OPEN(2,FILE=TRIM(ADJUSTL(FPOO)),STATUS='OLD',POSITION='APPEND') 58:          OPEN(2,FILE=TRIM(ADJUSTL(FPOO)),STATUS='OLD',POSITION='APPEND')
 64: C        PRINT '(A,A)','workaround for Sun compiler bug ZSYMBOL(1)=',ZSYMBOL(1) 59: C        PRINT '(A,A)','workaround for Sun compiler bug ZSYMBOL(1)=',ZSYMBOL(1)
 65:          WRITE(2,'(A2,2X,3F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2), 60:          WRITE(2,'(A2,2X,3F20.10)') (ZSYMBOL(J2),LOCALPOINTS1(3*(J2-1)+1),LOCALPOINTS1(3*(J2-1)+2),
 66:      1                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS) 61:      1                               LOCALPOINTS1(3*(J2-1)+3),J2=1,NATOMS)
 67:          CLOSE(2) 62:          CLOSE(2)
 68:       ENDIF 63:       ENDIF
 69:  64: 
 70:       RETURN 65:       RETURN


r29929/setup.f 2016-02-09 15:30:13.682698540 +0000 r29928/setup.f 2016-02-09 15:30:15.918728075 +0000
 74: !        the 'NOLABELS' option can be used to read just the coordinates 74: !        the 'NOLABELS' option can be used to read just the coordinates
 75:             IF (NOLABELST) THEN 75:             IF (NOLABELST) THEN
 76:                READ(1,'(A25,F20.10),A71') CDUMMY,EMIN(1) 76:                READ(1,'(A25,F20.10),A71') CDUMMY,EMIN(1)
 77:                READ(1,*) (LOCALPOINTS(J2),J2=1,NOPT) 77:                READ(1,*) (LOCALPOINTS(J2),J2=1,NOPT)
 78:             ELSE 78:             ELSE
 79:                READ(1,'(A25,F20.10,A71)') CDUMMY,EMIN(1),CDUMMY 79:                READ(1,'(A25,F20.10,A71)') CDUMMY,EMIN(1),CDUMMY
 80:                DO J1=1,NATOMS 80:                DO J1=1,NATOMS
 81:                   READ(1,*) CDUMMY,LOCALPOINTS(3*J1-2),LOCALPOINTS(3*J1-1),LOCALPOINTS(3*J1) 81:                   READ(1,*) CDUMMY,LOCALPOINTS(3*J1-2),LOCALPOINTS(3*J1-1),LOCALPOINTS(3*J1)
 82:                ENDDO 82:                ENDDO
 83:             ENDIF 83:             ENDIF
 84:             CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,IXM,IYM,IZM) 84:             CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,angleAxis,IXM,IYM,IZM)
 85:             WRITE(LUNIT,REC=NDUMMY) (LOCALPOINTS(J2),J2=1,NOPT) 85:             WRITE(LUNIT,REC=NDUMMY) (LOCALPOINTS(J2),J2=1,NOPT)
 86:             WRITE(LUNIT2,'(2F25.15,I6,3F20.10)') EMIN(1), FVIBMIN(1), HORDERMIN(1),IXM,IYM,IZM 86:             WRITE(LUNIT2,'(2F25.15,I6,3F20.10)') EMIN(1), FVIBMIN(1), HORDERMIN(1),IXM,IYM,IZM
 87:          ENDDO 87:          ENDDO
 88: 776      NDUMMY=NDUMMY-1 88: 776      NDUMMY=NDUMMY-1
 89:          WRITE(*,'(A,I10)')'setup> number of minima written=',NDUMMY 89:          WRITE(*,'(A,I10)')'setup> number of minima written=',NDUMMY
 90:          CLOSE(LUNIT) 90:          CLOSE(LUNIT)
 91:          CLOSE(LUNIT2) 91:          CLOSE(LUNIT2)
 92:          CLOSE(1) 92:          CLOSE(1)
 93:          STOP 93:          STOP
 94:       ENDIF 94:       ENDIF
107:             ENDDO107:             ENDDO
108: 777         NDUMMY=NDUMMY-1108: 777         NDUMMY=NDUMMY-1
109:             PRINT '(A,I10)','setup> number of minima extracted=',NDUMMY109:             PRINT '(A,I10)','setup> number of minima extracted=',NDUMMY
110:             CLOSE(LUNIT)110:             CLOSE(LUNIT)
111:          ELSEIF (WHICHMIN.LE.0) THEN111:          ELSEIF (WHICHMIN.LE.0) THEN
112:             NDUMMY=0112:             NDUMMY=0
113:             PRINT '(A)', 'setup> extracting all minima '113:             PRINT '(A)', 'setup> extracting all minima '
114:             DO 114:             DO 
115:                NDUMMY=NDUMMY+1115:                NDUMMY=NDUMMY+1
116:                READ(UMIN,REC=NDUMMY,ERR=877) (LOCALPOINTS(J2),J2=1,NOPT)116:                READ(UMIN,REC=NDUMMY,ERR=877) (LOCALPOINTS(J2),J2=1,NOPT)
117:                IF (MLP3T.OR.MLPB3T) THEN117:                WRITE(1,'(3F25.15)') (LOCALPOINTS(J2),J2=1,NOPT)
118:                   WRITE(1,'(F25.15)') (LOCALPOINTS(J2),J2=1,NOPT) 
119:                ELSE 
120:                   WRITE(1,'(3F25.15)') (LOCALPOINTS(J2),J2=1,NOPT) 
121:                ENDIF 
122:             ENDDO118:             ENDDO
123: 877         NDUMMY=NDUMMY-1119: 877         NDUMMY=NDUMMY-1
124:             PRINT '(A,I10)','setup> number of minima extracted=',NDUMMY120:             PRINT '(A,I10)','setup> number of minima extracted=',NDUMMY
125:          ELSE121:          ELSE
126:             PRINT '(A,I10)', 'setup> extracting minimum ',WHICHMIN122:             PRINT '(A,I10)', 'setup> extracting minimum ',WHICHMIN
127:             READ(UMIN,REC=WHICHMIN) (LOCALPOINTS(J2),J2=1,NOPT)123:             READ(UMIN,REC=WHICHMIN) (LOCALPOINTS(J2),J2=1,NOPT)
128: 124: 
129:            IF (AMHT) THEN125:            IF (AMHT) THEN
130:               CALL AMHDUMP(LOCALPOINTS,'amhmin.pdb')126:               CALL AMHDUMP(LOCALPOINTS,'amhmin.pdb')
131:   127:   
202:             ENDDO198:             ENDDO
203: 778         NDUMMY=NDUMMY-1199: 778         NDUMMY=NDUMMY-1
204:             CLOSE(LUNIT)200:             CLOSE(LUNIT)
205:             PRINT '(A,I10)','setup> number of ts extracted=',NDUMMY201:             PRINT '(A,I10)','setup> number of ts extracted=',NDUMMY
206:          ELSEIF (WHICHTS.LE.0) THEN202:          ELSEIF (WHICHTS.LE.0) THEN
207:             NDUMMY=0203:             NDUMMY=0
208:             PRINT '(A)', 'setup> extracting all ts '204:             PRINT '(A)', 'setup> extracting all ts '
209:             DO205:             DO
210:                NDUMMY=NDUMMY+1206:                NDUMMY=NDUMMY+1
211:                READ(UTS,REC=NDUMMY,ERR=878) (LOCALPOINTS(J2),J2=1,NOPT)207:                READ(UTS,REC=NDUMMY,ERR=878) (LOCALPOINTS(J2),J2=1,NOPT)
212:                IF (MLP3T.OR.MLPB3T) THEN208:                WRITE(1,'(3F25.15)') (LOCALPOINTS(J2),J2=1,NOPT)
213:                   WRITE(1,'(F25.15)') (LOCALPOINTS(J2),J2=1,NOPT) 
214:                ELSE 
215:                   WRITE(1,'(3F25.15)') (LOCALPOINTS(J2),J2=1,NOPT) 
216:                ENDIF 
217:             ENDDO209:             ENDDO
218: 878         NDUMMY=NDUMMY-1210: 878         NDUMMY=NDUMMY-1
219:             PRINT '(A,I10)','setup> number of ts extracted=',NDUMMY211:             PRINT '(A,I10)','setup> number of ts extracted=',NDUMMY
220:          ELSE212:          ELSE
221:             PRINT '(A,I10)', 'setup> extracting ts ',WHICHTS213:             PRINT '(A,I10)', 'setup> extracting ts ',WHICHTS
222:             READ(UTS,REC=WHICHTS) (LOCALPOINTS(J2),J2=1,NOPT)214:             READ(UTS,REC=WHICHTS) (LOCALPOINTS(J2),J2=1,NOPT)
223: 215: 
224: !           IF (AMHT) THEN216: !           IF (AMHT) THEN
225: !              CALL AMHDUMP(LOCALPOINTS,'amhts.pdb')217: !              CALL AMHDUMP(LOCALPOINTS,'amhts.pdb')
226: !           ELSE218: !           ELSE
546:       ENDIF538:       ENDIF
547: C539: C
548: C  Check that the necessary coordinates are in fact present. 540: C  Check that the necessary coordinates are in fact present. 
549: C541: C
550:       IF ((.NOT.NOPOINTS).AND.(NATTEMPT.GT.0)) THEN542:       IF ((.NOT.NOPOINTS).AND.(NATTEMPT.GT.0)) THEN
551:           IF (AMHT) THEN543:           IF (AMHT) THEN
552:             WRITE(*,*) 'setup> AVOIDING MOMENT OF INITERIA CALC FOR AMH'544:             WRITE(*,*) 'setup> AVOIDING MOMENT OF INITERIA CALC FOR AMH'
553:           ELSE545:           ELSE
554:            DO J1=1,NMIN546:            DO J1=1,NMIN
555:             READ(UMIN,REC=J1) (LOCALPOINTS(J2),J2=1,NOPT)547:             READ(UMIN,REC=J1) (LOCALPOINTS(J2),J2=1,NOPT)
556:             CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,IXM,IYM,IZM)548:             CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,angleAxis,IXM,IYM,IZM)
557:             IF (PRINTT) WRITE(*,'(2F20.10,I6,3F20.10)') EMIN(J1),FVIBMIN(J1),HORDERMIN(J1),IXM,IYM,IZM549: C           IF (PRINTT) WRITE(*,'(2F20.10,I6,3F20.10)') EMIN(J1),FVIBMIN(J1),HORDERMIN(J1),IXM,IYM,IZM
558:             IF ((ABS(IXM-IXMIN(J1)).GT.IDIFFTOL).OR.550:             IF ((ABS(IXM-IXMIN(J1)).GT.IDIFFTOL).OR.
559:      1          (ABS(IYM-IYMIN(J1)).GT.IDIFFTOL).OR.551:      1          (ABS(IYM-IYMIN(J1)).GT.IDIFFTOL).OR.
560:      1          (ABS(IZM-IZMIN(J1)).GT.IDIFFTOL)) THEN552:      1          (ABS(IZM-IZMIN(J1)).GT.IDIFFTOL)) THEN
561:                WRITE(*,'(A,I10)') 'setup> WARNING - principal moments of inertia do not agree with input for minimum ',J1553:                WRITE(*,'(A,I10)') 'setup> WARNING - principal moments of inertia do not agree with input for minimum ',J1
562:                WRITE(*,'(A,3F20.10)') 'setup> values from coordinates: ',IXM,IYM,IZM554:                WRITE(*,'(A,3F20.10)') 'setup> values from coordinates: ',IXM,IYM,IZM
563:                WRITE(*,'(A,2F20.10,I6,3F20.10)') 'setup> values from min.data: ', 555:                WRITE(*,'(A,2F20.10,I6,3F20.10)') 'setup> values from min.data: ', 
564:      &                     EMIN(J1),FVIBMIN(J1),HORDERMIN(J1),IXMIN(J1),IYMIN(J1),IZMIN(J1)556:      &                     EMIN(J1),FVIBMIN(J1),HORDERMIN(J1),IXMIN(J1),IYMIN(J1),IZMIN(J1)
 557:                PRINT '(A)','LOCALPOINTS:'
 558:                PRINT '(3G20.10)',LOCALPOINTS(1:NOPT)
565:                IXMIN(J1)=IXM559:                IXMIN(J1)=IXM
566:                IYMIN(J1)=IYM560:                IYMIN(J1)=IYM
567:                IZMIN(J1)=IZM561:                IZMIN(J1)=IZM
568: !              STOP  562: !              STOP  
569:             ENDIF563:             ENDIF
570:            ENDDO564:            ENDDO
571:           ENDIF565:           ENDIF
572:          IF (PRINTT) WRITE(*,'(A,I10,A)') 'setup> points for ',NMIN,' minima read from file points.min'566:          IF (PRINTT) WRITE(*,'(A,I10,A)') 'setup> points for ',NMIN,' minima read from file points.min'
573:       ENDIF567:       ENDIF
574: !568: !
940: !              WRITE(*,'(A,I10,A,4I10)') 'setup> for ts ',J1,' +,-,p+,p-:',PLUS(J1),MINUS(J1),POINTERP(J1),POINTERM(J1)934: !              WRITE(*,'(A,I10,A,4I10)') 'setup> for ts ',J1,' +,-,p+,p-:',PLUS(J1),MINUS(J1),POINTERP(J1),POINTERM(J1)
941: !           ENDDO935: !           ENDDO
942: !        ENDIF936: !        ENDIF
943:     937:     
944:          IF ((.NOT.NOPOINTS).AND.(NATTEMPT.GT.0).AND.(.NOT.DUMMYTST)) THEN938:          IF ((.NOT.NOPOINTS).AND.(NATTEMPT.GT.0).AND.(.NOT.DUMMYTST)) THEN
945:             IF (AMHT) THEN939:             IF (AMHT) THEN
946:               WRITE(*,*) 'setup> AVOIDING MOMENT OF INITERIA CALC FOR AMH'940:               WRITE(*,*) 'setup> AVOIDING MOMENT OF INITERIA CALC FOR AMH'
947:             ELSE941:             ELSE
948:                DO J1=1,NTS942:                DO J1=1,NTS
949:                   READ(UTS,REC=J1) (LOCALPOINTS(J2),J2=1,NOPT)943:                   READ(UTS,REC=J1) (LOCALPOINTS(J2),J2=1,NOPT)
950:                   CALL INERTIAWRAPPER(LOCALPOINTS,NOPT,angleAxis,IXM,IYM,IZM)944:                   CALL INERTIAWRAPPER(LOCALPOINTS,NATOMS,angleAxis,IXM,IYM,IZM)
951: !                 IF (DEBUG) WRITE(*,'(A,I10,2F17.7,3I6,3F15.5)') 'setup> ',J1,ETS(J1),FVIBTS(J1),HORDERTS(J1),945: !                 IF (DEBUG) WRITE(*,'(A,I10,2F17.7,3I6,3F15.5)') 'setup> ',J1,ETS(J1),FVIBTS(J1),HORDERTS(J1),
952: !    1                                                            PLUS(J1),MINUS(J1),IXM,IYM,IZM946: !    1                                                            PLUS(J1),MINUS(J1),IXM,IYM,IZM
953:                   IF ((ABS(IXM-IXTS(J1)).GT.IDIFFTOL).OR.947:                   IF ((ABS(IXM-IXTS(J1)).GT.IDIFFTOL).OR.
954:      1                (ABS(IYM-IYTS(J1)).GT.IDIFFTOL).OR.948:      1                (ABS(IYM-IYTS(J1)).GT.IDIFFTOL).OR.
955:      1                (ABS(IZM-IZTS(J1)).GT.IDIFFTOL)) THEN949:      1                (ABS(IZM-IZTS(J1)).GT.IDIFFTOL)) THEN
956:                      WRITE(*,'(A,I10)') 'setup> WARNING - principal moments of inertia do not agree with input for ts ',J1950:                      WRITE(*,'(A,I10)') 'setup> WARNING - principal moments of inertia do not agree with input for ts ',J1
957:                      WRITE(*,'(A)') 'values in ts.data:'951:                      WRITE(*,'(A)') 'values in ts.data:'
958:                      WRITE(*,'(3F20.10)') IXTS(J1),IYTS(J1),IZTS(J1)952:                      WRITE(*,'(3F20.10)') IXTS(J1),IYTS(J1),IZTS(J1)
959:                      WRITE(*,'(A)') 'values recalculated from points.ts:'953:                      WRITE(*,'(A)') 'values recalculated from points.ts:'
960:                      WRITE(*,'(3F20.10)') IXM,IYM,IZM954:                      WRITE(*,'(3F20.10)') IXM,IYM,IZM


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0