hdiff output

r33624/bfgsts.f 2018-01-19 09:34:32.870750079 +0000 r33623/bfgsts.f 2018-01-19 09:34:37.862790800 +0000
148:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH148:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
149:       SAVE149:       SAVE
150: 150: 
151:       ZWORK(:,1) = 0.0D0  ! sn402: Somehow we are arriving here with ZWORK already partially initialised.151:       ZWORK(:,1) = 0.0D0  ! sn402: Somehow we are arriving here with ZWORK already partially initialised.
152:       ! Surely this is a mistake? So I'm setting this to 0.152:       ! Surely this is a mistake? So I'm setting this to 0.
153: 153: 
154:       LWORK=33*3*NATOMS154:       LWORK=33*3*NATOMS
155:       ILWORK=33*3*NATOMS155:       ILWORK=33*3*NATOMS
156: 156: 
157:       IF ((ZSYM(1)(1:1).EQ.'W').AND.(.NOT.BFGSSTEP)) THEN157:       IF ((ZSYM(1)(1:1).EQ.'W').AND.(.NOT.BFGSSTEP)) THEN
 158:         IF (ZSYM(1)(1:2).NE.'WC') THEN
158:          PRINT*,'BFGSTS procedures have not been programmed for TIP potentials'159:          PRINT*,'BFGSTS procedures have not been programmed for TIP potentials'
159:          CALL FLUSH(6)160:          CALL FLUSH(6)
160:          STOP161:          STOP
 162:         ENDIF
161:       ENDIF163:       ENDIF
162:       FIXDSAVE=FIXD164:       FIXDSAVE=FIXD
163:       IF ((HINDEX.GT.1).AND.(.NOT.NOIT)) THEN165:       IF ((HINDEX.GT.1).AND.(.NOT.NOIT)) THEN
164:          WRITE(*,'(A)') 'For HINDEX > 1 you must use NOIT with BFGSTS'166:          WRITE(*,'(A)') 'For HINDEX > 1 you must use NOIT with BFGSTS'
165:          CALL FLUSH(6)167:          CALL FLUSH(6)
166:          STOP168:          STOP
167:       ENDIF169:       ENDIF
168: !170: !
169: ! Check for consistent convergence criteria on RMS gradient in MYLBFGS and EF part.171: ! Check for consistent convergence criteria on RMS gradient in MYLBFGS and EF part.
170: !     172: !     


r33624/connect.f 2018-01-19 09:34:33.126752109 +0000 r33623/connect.f 2018-01-19 09:34:38.122792987 +0000
569:       IF (TSFRQDONE) THEN569:       IF (TSFRQDONE) THEN
570:          DO J1=1,NOPT570:          DO J1=1,NOPT
571:             FSAVETS(J1,NTS)=FRQSTS(J1)571:             FSAVETS(J1,NTS)=FRQSTS(J1)
572:          ENDDO572:          ENDDO
573:       ENDIF573:       ENDIF
574:    574:    
575:       IF (USEOLD) THEN575:       IF (USEOLD) THEN
576: C576: C
577: C  For water the paths files are in xyz format. We could convert instead of rerunning the path.577: C  For water the paths files are in xyz format. We could convert instead of rerunning the path.
578: C578: C
579:          IF (ZSYMSAVE(1:1).EQ.'W') THEN579:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
580:             CALL MKFNAMES(JUSE,FILTH,FILTHSTR,ITSTRING,EOFSSTRING)580:             CALL MKFNAMES(JUSE,FILTH,FILTHSTR,ITSTRING,EOFSSTRING)
581: C           CALL PATH(Q,ENERGY,VNEW,RMS,EVALMIN,VECS,.FALSE.,QPLUS,QMINUS,PTEST,ETS,581: C           CALL PATH(Q,ENERGY,VNEW,RMS,EVALMIN,VECS,.FALSE.,QPLUS,QMINUS,PTEST,ETS,
582: C    1             PLUSEN(JUSE),MINUSEN(JUSE),PATHLENGTH(JUSE),DISP(JUSE),GAMMA(JUSE),NTILDE(JUSE),582: C    1             PLUSEN(JUSE),MINUSEN(JUSE),PATHLENGTH(JUSE),DISP(JUSE),GAMMA(JUSE),NTILDE(JUSE),
583: C    2             FRQSTS,FRQSPLUS,FRQSMINUS,ITSTRING,EOFSSTRING,PATHFAILT)583: C    2             FRQSTS,FRQSPLUS,FRQSMINUS,ITSTRING,EOFSSTRING,PATHFAILT)
584:             TSUSED(JUSE)=.TRUE.584:             TSUSED(JUSE)=.TRUE.
585:             OPEN(UNIT=45,FILE=ITSTRING,STATUS='OLD')585:             OPEN(UNIT=45,FILE=ITSTRING,STATUS='OLD')
586:             READ(45,*)586:             READ(45,*)
587:             READ(45,*)587:             READ(45,*)
588: C           DO J1=1,NATOMS  !  WCOMMENT588: C           DO J1=1,NATOMS  !  WCOMMENT
589:             DO J1=1,NATOMS/2589:             DO J1=1,NATOMS/2
955:       DOUBLE PRECISION CMXA, CMYA, CMZA955:       DOUBLE PRECISION CMXA, CMYA, CMZA
956:       CHARACTER(LEN=5) ZSYMSAVE956:       CHARACTER(LEN=5) ZSYMSAVE
957:       CHARACTER(LEN=5), ALLOCATABLE, DIMENSION(:) :: ZSYML957:       CHARACTER(LEN=5), ALLOCATABLE, DIMENSION(:) :: ZSYML
958:       COMMON /SYS/ ZSYMSAVE958:       COMMON /SYS/ ZSYMSAVE
959:       INTEGER K1,K2,KD,NNZ,NINTB ! jmc for unres959:       INTEGER K1,K2,KD,NNZ,NINTB ! jmc for unres
960:       LOGICAL UNRST ! jmc for unres960:       LOGICAL UNRST ! jmc for unres
961: 961: 
962: C     ZSYMSAVE=ZSYM(1)962: C     ZSYMSAVE=ZSYM(1)
963: C     IMUL=1   !  WCOMMENT963: C     IMUL=1   !  WCOMMENT
964: C     IF (ZSYMSAVE(1:1).EQ.'W') IMUL=3   !  WCOMMENT964: C     IF (ZSYMSAVE(1:1).EQ.'W') IMUL=3   !  WCOMMENT
965:       IF (ZSYMSAVE(1:1).EQ.'W') THEN965:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
966:          NATOMSIMUL=(NATOMS/2)*3  !  number of atoms for water966:          NATOMSIMUL=(NATOMS/2)*3  !  number of atoms for water
967:          ALLOCATE(COORDS(9*(NATOMS/2)),PCOORDS(9*(NATOMS/2)),QW(9*(NATOMS/2)),ZSYML(3*(NATOMS/2)))967:          ALLOCATE(COORDS(9*(NATOMS/2)),PCOORDS(9*(NATOMS/2)),QW(9*(NATOMS/2)),ZSYML(3*(NATOMS/2)))
968:       ELSE968:       ELSE
969:          NATOMSIMUL=NATOMS969:          NATOMSIMUL=NATOMS
970:          ALLOCATE(COORDS(3*NATOMS),PCOORDS(3*NATOMS),ZSYML(NATOMS))970:          ALLOCATE(COORDS(3*NATOMS),PCOORDS(3*NATOMS),ZSYML(NATOMS))
971:          IF (UNRST) ALLOCATE(PEPCOORDS(3*NATOMS))971:          IF (UNRST) ALLOCATE(PEPCOORDS(3*NATOMS))
972:       ENDIF972:       ENDIF
973: 973: 
974:       WRITE(*,'(A)') 'Connected path found'974:       WRITE(*,'(A)') 'Connected path found'
975:       WRITE(*,'(A)')975:       WRITE(*,'(A)')
1031: C  two minima should be identical of course! To put the other frames in the path into the same orientation1031: C  two minima should be identical of course! To put the other frames in the path into the same orientation
1032: C  we need to save the overall rotation matrix from mind and apply the same rotation - this is done in1032: C  we need to save the overall rotation matrix from mind and apply the same rotation - this is done in
1033: C  ROTGEOM. 1033: C  ROTGEOM. 
1034: C1034: C
1035: C  For W4/W3/W2/W1 COORDS contains the Cartesian coordinates, so we don;t want1035: C  For W4/W3/W2/W1 COORDS contains the Cartesian coordinates, so we don;t want
1036: C  to convert them again in MIND. However, ZSYM(1) has been overwritten at this1036: C  to convert them again in MIND. However, ZSYM(1) has been overwritten at this
1037: C  point, so it won;t contain W any more. Need to turn RIGIDBODY off in the call for 1037: C  point, so it won;t contain W any more. Need to turn RIGIDBODY off in the call for 
1038: C  water.1038: C  water.
1039: C1039: C
1040:          DOREVERSE=.TRUE.1040:          DOREVERSE=.TRUE.
1041:          IF (ZSYMSAVE(1:1).EQ.'W') THEN1041:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1042:             CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,.FALSE.,DEBUG,RMAT)1042:             CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,.FALSE.,DEBUG,RMAT)
1043:          ELSE1043:          ELSE
1044:             CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,RIGIDBODY,DEBUG,RMAT)1044:             CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,RIGIDBODY,DEBUG,RMAT)
1045:          ENDIF1045:          ENDIF
1046: C        PRINT '(A,G20.10)','DISTPF=',DISTPF1046: C        PRINT '(A,G20.10)','DISTPF=',DISTPF
1047: C        PRINT '(A,9F15.10)','RMAT=',RMAT(1:3,1:3)1047: C        PRINT '(A,9F15.10)','RMAT=',RMAT(1:3,1:3)
1048: C        IF (DEBUG) PRINT '(A,G20.10)','Initial DISTPF=',DISTPF1048: C        IF (DEBUG) PRINT '(A,G20.10)','Initial DISTPF=',DISTPF
1049:          IF (DISTPF.LE.GEOMDIFFTOL) THEN ! there should be no need to reverse the path!1049:          IF (DISTPF.LE.GEOMDIFFTOL) THEN ! there should be no need to reverse the path!
1050:                                          ! unless both distances are less than DISTPF, sigh1050:                                          ! unless both distances are less than DISTPF, sigh
1051:             DOREVERSE=.FALSE.1051:             DOREVERSE=.FALSE.
1075:          IF (DOREVERSE) THEN1075:          IF (DOREVERSE) THEN
1076:             CALL REVERSE(84)1076:             CALL REVERSE(84)
1077:             REWIND(85)1077:             REWIND(85)
1078:             CALL REVERSEBLOCK(85,NATOMSIMUL+2)1078:             CALL REVERSEBLOCK(85,NATOMSIMUL+2)
1079:             READ(85,'(A100)') DUMMYST1079:             READ(85,'(A100)') DUMMYST
1080:             READ(85,'(A100)') DUMMYST21080:             READ(85,'(A100)') DUMMYST2
1081:             DO J4=1,NATOMSIMUL1081:             DO J4=1,NATOMSIMUL
1082:                READ(85,'(A2,4X,3F20.10)') ZSYML(J4),COORDS(3*(J4-1)+1),COORDS(3*(J4-1)+2),COORDS(3*(J4-1)+3)1082:                READ(85,'(A2,4X,3F20.10)') ZSYML(J4),COORDS(3*(J4-1)+1),COORDS(3*(J4-1)+2),COORDS(3*(J4-1)+3)
1083:             ENDDO1083:             ENDDO
1084: 1084: 
1085:             IF (ZSYMSAVE(1:1).EQ.'W') THEN1085:             IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1086:                CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,.FALSE.,DEBUG,RMAT)1086:                CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,.FALSE.,DEBUG,RMAT)
1087:             ELSE1087:             ELSE
1088:                CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,RIGIDBODY,DEBUG,RMAT)1088:                CALL NEWMINDIST(PCOORDS,COORDS,NATOMSIMUL,DISTPF,BULKT,TWOD,ZSYML(1),.FALSE.,RIGIDBODY,DEBUG,RMAT)
1089:             ENDIF1089:             ENDIF
1090: C           PRINT '(A,G20.10)','DISTPF=',DISTPF1090: C           PRINT '(A,G20.10)','DISTPF=',DISTPF
1091: C           PRINT '(A,9F15.10)','RMAT=',RMAT(1:3,1:3)1091: C           PRINT '(A,9F15.10)','RMAT=',RMAT(1:3,1:3)
1092: 1092: 
1093: C           IF (DEBUG) PRINT '(A,G20.10,A)','DISTPF for reversed path=',DISTPF ! DJW1093: C           IF (DEBUG) PRINT '(A,G20.10,A)','DISTPF for reversed path=',DISTPF ! DJW
1094:             IF ((DISTPF.GT.GEOMDIFFTOL).OR.(ABS(MINUSEN(J2)-EMIN(J1)).GT.EDIFFTOL)) THEN1094:             IF ((DISTPF.GT.GEOMDIFFTOL).OR.(ABS(MINUSEN(J2)-EMIN(J1)).GT.EDIFFTOL)) THEN
1095:                PRINT'(A)','ERROR - distance and energy tolerances are not consistent for reversed path'1095:                PRINT'(A)','ERROR - distance and energy tolerances are not consistent for reversed path'
1228:             IF (UNRST.AND.CALCDIHE) THEN1228:             IF (UNRST.AND.CALCDIHE) THEN
1229:                DO J2=1,3*NATOMS1229:                DO J2=1,3*NATOMS
1230:                   QMIN(J2)=QSAVE(J2,J1)1230:                   QMIN(J2)=QSAVE(J2,J1)
1231:                ENDDO1231:                ENDDO
1232:                CALL UNRESCALCDIHEREF(DIHE,ALLANG,QMIN)1232:                CALL UNRESCALCDIHEREF(DIHE,ALLANG,QMIN)
1233: C               WRITE(83,'(3F25.15)') EMIN(J1), DIHE, ALLANG1233: C               WRITE(83,'(3F25.15)') EMIN(J1), DIHE, ALLANG
1234:                WRITE(83,'(3F25.15)') EMIN(J1), DIHE1234:                WRITE(83,'(3F25.15)') EMIN(J1), DIHE
1235:             ELSE1235:             ELSE
1236:                WRITE(83,'(F25.15)') EMIN(J1)1236:                WRITE(83,'(F25.15)') EMIN(J1)
1237:             ENDIF1237:             ENDIF
1238:             IF (ZSYMSAVE(1:1).EQ.'W') THEN1238:             IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1239:                IF (ZSYMSAVE.EQ.'W4') IPOT=41239:                IF (ZSYMSAVE.EQ.'W4') IPOT=4
1240:                IF (ZSYMSAVE.EQ.'W3') IPOT=31240:                IF (ZSYMSAVE.EQ.'W3') IPOT=3
1241:                IF (ZSYMSAVE.EQ.'W2') IPOT=21241:                IF (ZSYMSAVE.EQ.'W2') IPOT=2
1242:                IF (ZSYMSAVE.EQ.'W1') IPOT=11242:                IF (ZSYMSAVE.EQ.'W1') IPOT=1
1243: C              DO J2=1,NATOMS  !  WCOMMENT1243: C              DO J2=1,NATOMS  !  WCOMMENT
1244:                DO J2=1,NATOMS/21244:                DO J2=1,NATOMS/2
1245:                   CALL CONVERT(QSAVE(3*(J2-1)+1,J1),QSAVE(3*(J2-1)+2,J1),QSAVE(3*(J2-1)+3,J1),1245:                   CALL CONVERT(QSAVE(3*(J2-1)+1,J1),QSAVE(3*(J2-1)+2,J1),QSAVE(3*(J2-1)+3,J1),
1246: C ! WCOMMENT1246: C ! WCOMMENT
1247: C    1                         QSAVE(3*(NATOMS+J2-1)+1,J1),QSAVE(3*(NATOMS+J2-1)+2,J1),QSAVE(3*(NATOMS+J2-1)+3,J1),1247: C    1                         QSAVE(3*(NATOMS+J2-1)+1,J1),QSAVE(3*(NATOMS+J2-1)+2,J1),QSAVE(3*(NATOMS+J2-1)+3,J1),
1248:      1                         QSAVE(3*(NATOMS/2+J2-1)+1,J1),QSAVE(3*(NATOMS/2+J2-1)+2,J1),QSAVE(3*(NATOMS/2+J2-1)+3,J1),1248:      1                         QSAVE(3*(NATOMS/2+J2-1)+1,J1),QSAVE(3*(NATOMS/2+J2-1)+2,J1),QSAVE(3*(NATOMS/2+J2-1)+3,J1),
1350:                ENDDO1350:                ENDDO
1351:                CALL SYMMETRY(HORDER,.FALSE.,Q,INERTIA)1351:                CALL SYMMETRY(HORDER,.FALSE.,Q,INERTIA)
1352:                WRITE(83,'(I6,1X,A4)') HORDER,FPGRP1352:                WRITE(83,'(I6,1X,A4)') HORDER,FPGRP
1353:                IF (.NOT.NOFRQS) WRITE(83,'(3G20.10)') (FSAVEMIN(J2,J1),J2=1,3*NATOMS)1353:                IF (.NOT.NOFRQS) WRITE(83,'(3G20.10)') (FSAVEMIN(J2,J1),J2=1,3*NATOMS)
1354:             ENDIF1354:             ENDIF
1355:             IF (J1.GT.1) CALL NEWMINDIST(QTS,QMIN,NATOMS,DISTPF,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)1355:             IF (J1.GT.1) CALL NEWMINDIST(QTS,QMIN,NATOMS,DISTPF,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)
1356:             WRITE(83,'(3F25.15)') (QMIN(J2),J2=1,NOPT)1356:             WRITE(83,'(3F25.15)') (QMIN(J2),J2=1,NOPT)
1357:             IF (J1.LT.NMIN) THEN1357:             IF (J1.LT.NMIN) THEN
1358:                J3=WHICHTS(J1)1358:                J3=WHICHTS(J1)
1359:                WRITE(83,'(F25.15)') TSEN(J3)1359:                WRITE(83,'(F25.15)') TSEN(J3)
1360:                IF (ZSYMSAVE(1:1).EQ.'W') THEN1360:                IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1361:                   IF (ZSYMSAVE.EQ.'W4') IPOT=41361:                   IF (ZSYMSAVE.EQ.'W4') IPOT=4
1362:                   IF (ZSYMSAVE.EQ.'W3') IPOT=31362:                   IF (ZSYMSAVE.EQ.'W3') IPOT=3
1363:                   IF (ZSYMSAVE.EQ.'W2') IPOT=21363:                   IF (ZSYMSAVE.EQ.'W2') IPOT=2
1364:                   IF (ZSYMSAVE.EQ.'W1') IPOT=11364:                   IF (ZSYMSAVE.EQ.'W1') IPOT=1
1365: C                 DO J2=1,NATOMS  !  WCOMMENT1365: C                 DO J2=1,NATOMS  !  WCOMMENT
1366:                   DO J2=1,NATOMS/21366:                   DO J2=1,NATOMS/2
1367:                      CALL CONVERT(QSAVETS(3*(J2-1)+1,J3),QSAVETS(3*(J2-1)+2,J3),QSAVETS(3*(J2-1)+3,J3),1367:                      CALL CONVERT(QSAVETS(3*(J2-1)+1,J3),QSAVETS(3*(J2-1)+2,J3),QSAVETS(3*(J2-1)+3,J3),
1368: C    1                            QSAVETS(3*(NATOMS+J2-1)+1,J3),QSAVETS(3*(NATOMS+J2-1)+2,J3),QSAVETS(3*(NATOMS+J2-1)+3,J3),1368: C    1                            QSAVETS(3*(NATOMS+J2-1)+1,J3),QSAVETS(3*(NATOMS+J2-1)+2,J3),QSAVETS(3*(NATOMS+J2-1)+3,J3),
1369:      1                            QSAVETS(3*(NATOMS/2+J2-1)+1,J3),QSAVETS(3*(NATOMS/2+J2-1)+2,J3),QSAVETS(3*(NATOMS/2+J2-1)+3,J3),1369:      1                            QSAVETS(3*(NATOMS/2+J2-1)+1,J3),QSAVETS(3*(NATOMS/2+J2-1)+2,J3),QSAVETS(3*(NATOMS/2+J2-1)+3,J3),
1370:      2                            OVEC,H1VEC,H2VEC)1370:      2                            OVEC,H1VEC,H2VEC)


r33624/dumpp.f 2018-01-19 09:34:33.402754309 +0000 r33623/dumpp.f 2018-01-19 09:34:38.362795018 +0000
 36:       WRITE(2,'(G30.20)') ENERGY 36:       WRITE(2,'(G30.20)') ENERGY
 37:       S=1.0D0 37:       S=1.0D0
 38:       IF ((ZSYM(NATOMS).EQ.'AR').AND.(ZSYM(1).NE.'CA')) S=3.4D0 38:       IF ((ZSYM(NATOMS).EQ.'AR').AND.(ZSYM(1).NE.'CA')) S=3.4D0
 39: C     IF (ZSYM(1)(1:1).EQ.'W') THEN 39: C     IF (ZSYM(1)(1:1).EQ.'W') THEN
 40: C        WRITE(1,'(I5)') 3*NATOMS 40: C        WRITE(1,'(I5)') 3*NATOMS
 41: C        WRITE(1,'(F20.10)') ENERGY 41: C        WRITE(1,'(F20.10)') ENERGY
 42: C     ENDIF 42: C     ENDIF
 43:       DO J1=1,NATOMS 43:       DO J1=1,NATOMS
 44:          IF (VARIABLES) THEN 44:          IF (VARIABLES) THEN
 45:             WRITE(1,'(G20.10)') QL(J1) 45:             WRITE(1,'(G20.10)') QL(J1)
 46:          ELSE IF (ZSYM(1)(1:1).EQ.'W') THEN 46:          ELSE IF (ZSYM(1)(1:1).EQ.'W'.and.ZSYM(1)(1:2).NE.'WC') THEN
 47:             IF (J1.LE.NATOMS/2) THEN  !  WCOMMENT new line 47:             IF (J1.LE.NATOMS/2) THEN  !  WCOMMENT new line
 48:             CALL CONVERT(QL(3*(J1-1)+1),QL(3*(J1-1)+2),QL(3*(J1-1)+3), 48:             CALL CONVERT(QL(3*(J1-1)+1),QL(3*(J1-1)+2),QL(3*(J1-1)+3),
 49:      1            QL(3*(NATOMS/2+J1-1)+1),QL(3*(NATOMS/2+J1-1)+2),QL(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC) 49:      1            QL(3*(NATOMS/2+J1-1)+1),QL(3*(NATOMS/2+J1-1)+2),QL(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC)
 50: C  ! WCOMMENT 50: C  ! WCOMMENT
 51: C    1            QL(3*(NATOMS+J1-1)+1),QL(3*(NATOMS+J1-1)+2),QL(3*(NATOMS+J1-1)+3),OVEC,H1VEC,H2VEC) 51: C    1            QL(3*(NATOMS+J1-1)+1),QL(3*(NATOMS+J1-1)+2),QL(3*(NATOMS+J1-1)+3),OVEC,H1VEC,H2VEC)
 52: C           WRITE(1,'(A3,3G20.10)') 'O  ',OVEC(1),OVEC(2),OVEC(3) 52: C           WRITE(1,'(A3,3G20.10)') 'O  ',OVEC(1),OVEC(2),OVEC(3)
 53: C           WRITE(1,'(A3,3G20.10)') 'H  ',H1VEC(1),H1VEC(2),H1VEC(3) 53: C           WRITE(1,'(A3,3G20.10)') 'H  ',H1VEC(1),H1VEC(2),H1VEC(3)
 54: C           WRITE(1,'(A3,3G20.10)') 'H  ',H2VEC(1),H2VEC(2),H2VEC(3) 54: C           WRITE(1,'(A3,3G20.10)') 'H  ',H2VEC(1),H2VEC(2),H2VEC(3)
 55:             WRITE(1,'(2X,3G20.10)') OVEC(1),OVEC(2),OVEC(3) 55:             WRITE(1,'(2X,3G20.10)') OVEC(1),OVEC(2),OVEC(3)
 56:             WRITE(1,'(2X,3G20.10)') H1VEC(1),H1VEC(2),H1VEC(3) 56:             WRITE(1,'(2X,3G20.10)') H1VEC(1),H1VEC(2),H1VEC(3)


r33624/efol.f90 2018-01-19 09:34:33.662756381 +0000 r33623/efol.f90 2018-01-19 09:34:38.630797286 +0000
109: !109: !
110: !     CALL GMETRY(0,VEC,QTS)110: !     CALL GMETRY(0,VEC,QTS)
111:       VEC(1:NOPT)=0.0D0111:       VEC(1:NOPT)=0.0D0
112:       MFLAG=.FALSE.112:       MFLAG=.FALSE.
113:       MINTEST=.FALSE.113:       MINTEST=.FALSE.
114:       TSTEST=.FALSE.114:       TSTEST=.FALSE.
115:       SDTEST=.FALSE.115:       SDTEST=.FALSE.
116:       NRTEST=.FALSE.116:       NRTEST=.FALSE.
117:       FRAME=1117:       FRAME=1
118: 118: 
119:       IF (ZSYMSAVE(1:1).EQ.'W') ALLOCATE(QW(9*(NATOMS/2)))119:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') ALLOCATE(QW(9*(NATOMS/2)))
120:       IF (INR.EQ.0) THEN120:       IF (INR.EQ.0) THEN
121:          IF (PTEST) WRITE(*,10) 121:          IF (PTEST) WRITE(*,10) 
122:          MINTEST=.TRUE.122:          MINTEST=.TRUE.
123:       ELSE IF (INR.EQ.1) THEN123:       ELSE IF (INR.EQ.1) THEN
124:          IF (PTEST) WRITE(*,20) 124:          IF (PTEST) WRITE(*,20) 
125:          NRTEST=.TRUE.125:          NRTEST=.TRUE.
126:       ELSE IF (INR.EQ.2) THEN126:       ELSE IF (INR.EQ.2) THEN
127:          IF (PTEST) WRITE(*,10) 127:          IF (PTEST) WRITE(*,10) 
128:          TSTEST=.TRUE.128:          TSTEST=.TRUE.
129:       ELSE IF (INR.EQ.3) THEN129:       ELSE IF (INR.EQ.3) THEN
1107:       ENDIF1107:       ENDIF
1108:       IF (TEST1.AND.TEST2) MFLAG=.TRUE.1108:       IF (TEST1.AND.TEST2) MFLAG=.TRUE.
1109:       IF (ITER.LT.NSTEPMIN) MFLAG=.FALSE.1109:       IF (ITER.LT.NSTEPMIN) MFLAG=.FALSE.
1110:       IF (PV.AND.(.NOT.PVFLAG)) MFLAG=.FALSE.1110:       IF (PV.AND.(.NOT.PVFLAG)) MFLAG=.FALSE.
1111: !1111: !
1112: !  Don't call symmetry if we're doing Fenske-Hall, or if the RMS force is too high.1112: !  Don't call symmetry if we're doing Fenske-Hall, or if the RMS force is too high.
1113: !1113: !
1114:       IF ((ZSYM(NATOMS).NE.'FH').AND.(.NOT.VARIABLES).AND.(.NOT.FIELDT).AND.(.NOT.RINGPOLYMERT) &1114:       IF ((ZSYM(NATOMS).NE.'FH').AND.(.NOT.VARIABLES).AND.(.NOT.FIELDT).AND.(.NOT.RINGPOLYMERT) &
1115:      &        .AND.(.NOT.PHI4MODT) &1115:      &        .AND.(.NOT.PHI4MODT) &
1116:      &        .AND.(PTEST).AND.(((RMS.LT.SYMCUT).OR.(ITER.EQ.NSTEPS)).OR.MFLAG)) THEN1116:      &        .AND.(PTEST).AND.(((RMS.LT.SYMCUT).OR.(ITER.EQ.NSTEPS)).OR.MFLAG)) THEN
1117:          IF (ZSYMSAVE(1:1).EQ.'W') THEN1117:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1118: !           DO J2=1,NOPT1118: !           DO J2=1,NOPT
1119: !              QSAVE(J2)=QTS(J2)1119: !              QSAVE(J2)=QTS(J2)
1120: !           ENDDO1120: !           ENDDO
1121: !           DO J2=1,NATOMS  !  WCOMMENT1121: !           DO J2=1,NATOMS  !  WCOMMENT
1122:             DO J2=1,NATOMS/21122:             DO J2=1,NATOMS/2
1123:                CALL CONVERT(QTS(3*(J2-1)+1),QTS(3*(J2-1)+2),QTS(3*(J2-1)+3),&1123:                CALL CONVERT(QTS(3*(J2-1)+1),QTS(3*(J2-1)+2),QTS(3*(J2-1)+3),&
1124: !    &                      QTS(3*(NATOMS+J2-1)+1),QTS(3*(NATOMS+J2-1)+2),QTS(3*(NATOMS+J2-1)+3),&1124: !    &                      QTS(3*(NATOMS+J2-1)+1),QTS(3*(NATOMS+J2-1)+2),QTS(3*(NATOMS+J2-1)+3),&
1125:      &                      QTS(3*(NATOMS/2+J2-1)+1),QTS(3*(NATOMS/2+J2-1)+2),QTS(3*(NATOMS/2+J2-1)+3),&1125:      &                      QTS(3*(NATOMS/2+J2-1)+1),QTS(3*(NATOMS/2+J2-1)+2),QTS(3*(NATOMS/2+J2-1)+3),&
1126:      &                      OVEC,H1VEC,H2VEC)1126:      &                      OVEC,H1VEC,H2VEC)
1127:                QW(9*(J2-1)+1)=OVEC(1)1127:                QW(9*(J2-1)+1)=OVEC(1)


r33624/fetchz.f 2018-01-19 09:34:33.942758614 +0000 r33623/fetchz.f 2018-01-19 09:34:38.890799486 +0000
1208:             NATOMS=(NATOMS/3)*2 ! The great big fudge1208:             NATOMS=(NATOMS/3)*2 ! The great big fudge
1209:             LNATOMS=NATOMS1209:             LNATOMS=NATOMS
1210:             Q(1:3*NATOMS)=THTEMP(1:3*NATOMS) ! coordinates saved as theta, phi1210:             Q(1:3*NATOMS)=THTEMP(1:3*NATOMS) ! coordinates saved as theta, phi
1211:             PRINT '(A,I6)',' fetchz> Number of atoms changed for the Thomson problem to ',NATOMS1211:             PRINT '(A,I6)',' fetchz> Number of atoms changed for the Thomson problem to ',NATOMS
1212:          ENDIF1212:          ENDIF
1213: 1213: 
1214:          IF (MACHINE) THEN1214:          IF (MACHINE) THEN
1215:              CALL READINPFILE(Q)1215:              CALL READINPFILE(Q)
1216:          ENDIF1216:          ENDIF
1217: 1217: 
1218:          IF (ZSYM(NATOMS)(1:1).EQ.'W') THEN1218:          IF (ZSYM(NATOMS)(1:1).EQ.'W'.and.ZSYM(NATOMS)(1:2).NE.'WC' ) THEN
1219:             RIGIDBODY=.TRUE.1219:             RIGIDBODY=.TRUE.
1220: C            OPEN(UNIT=18,FILE='Euler',STATUS='OLD')   !   WCOMMENT1220: C            OPEN(UNIT=18,FILE='Euler',STATUS='OLD')   !   WCOMMENT
1221: C            READ(18,*) (Q(3*LNATOMS+J1),J1=1,3*LNATOMS)1221: C            READ(18,*) (Q(3*LNATOMS+J1),J1=1,3*LNATOMS)
1222: C            CLOSE(18)1222: C            CLOSE(18)
1223:          ENDIF1223:          ENDIF
1224: 1224: 
1225: ! hk2861225: ! hk286
1226: !         IF (GTHOMSONT) THEN1226: !         IF (GTHOMSONT) THEN
1227: !            IF ((NATOMS/3)*1.0D0.NE.(NATOMS*1.0D0/3)) THEN1227: !            IF ((NATOMS/3)*1.0D0.NE.(NATOMS*1.0D0/3)) THEN
1228: !               PRINT '(A)',' fetchz> ERROR, the number of atoms needs to be divisible by 3 for GThomson fudge'1228: !               PRINT '(A)',' fetchz> ERROR, the number of atoms needs to be divisible by 3 for GThomson fudge'
1520:          ELSE IF ( ZSYM(NATOMS).EQ.'SB') THEN1520:          ELSE IF ( ZSYM(NATOMS).EQ.'SB') THEN
1521:             WRITE(*,'(A,I6,A,I5,A)')1521:             WRITE(*,'(A,I6,A,I5,A)')
1522:      1       ' SYSTEM ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' SBM sites'1522:      1       ' SYSTEM ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' SBM sites'
1523:          ELSE IF (ZSYM(NATOMS).EQ.'P6') THEN1523:          ELSE IF (ZSYM(NATOMS).EQ.'P6') THEN
1524:             WRITE(*,'(A,I4,A,I4,A)') 1524:             WRITE(*,'(A,I4,A,I4,A)') 
1525:      1       ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' Girifalco C60 molecules'1525:      1       ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' Girifalco C60 molecules'
1526:             WRITE(*,'(A,F15.8,A,F15.8,A,F15.8,A,F15.8)') 1526:             WRITE(*,'(A,F15.8,A,F15.8,A,F15.8,A,F15.8)') 
1527:      1       ' fetchz> Box lengths: x ',PARAM1,', y ',PARAM2,', z ',PARAM3,', cutoff ',PARAM41527:      1       ' fetchz> Box lengths: x ',PARAM1,', y ',PARAM2,', z ',PARAM3,', cutoff ',PARAM4
1528:          ELSE IF (ZSYM(NATOMS).EQ.'FH') THEN1528:          ELSE IF (ZSYM(NATOMS).EQ.'FH') THEN
1529:             WRITE(*,'(A,I4,A,I4,A)') ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' Fenske-Hall atoms'1529:             WRITE(*,'(A,I4,A,I4,A)') ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' Fenske-Hall atoms'
1530:          ELSE IF (ZSYM(NATOMS)(1:1).EQ.'W') THEN1530:          ELSE IF (ZSYM(NATOMS)(1:1).EQ.'W'.and.ZSYM(NATOMS)(1:2).NE.'WC') THEN
1531:             IF (ZSYM(NATOMS).EQ.'W4') WRITE(*,'(A,I4,A,I4,A)') 1531:             IF (ZSYM(NATOMS).EQ.'W4') WRITE(*,'(A,I4,A,I4,A)') 
1532:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIP4P water molecules'1532:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIP4P water molecules'
1533:             IF (ZSYM(NATOMS).EQ.'W3') WRITE(*,'(A,I4,A,I4,A)') 1533:             IF (ZSYM(NATOMS).EQ.'W3') WRITE(*,'(A,I4,A,I4,A)') 
1534:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIP3P water molecules'1534:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIP3P water molecules'
1535:             IF (ZSYM(NATOMS).EQ.'W2') WRITE(*,'(A,I4,A,I4,A)') 1535:             IF (ZSYM(NATOMS).EQ.'W2') WRITE(*,'(A,I4,A,I4,A)') 
1536:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIPS2 water molecules'1536:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIPS2 water molecules'
1537:             IF (ZSYM(NATOMS).EQ.'W1') WRITE(*,'(A,I4,A,I4,A)') 1537:             IF (ZSYM(NATOMS).EQ.'W1') WRITE(*,'(A,I4,A,I4,A)') 
1538:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIPS1 water molecules'1538:      2       ' fetchz> ',3*NATOMS,' coordinates will be optimised for ',NATOMS/2,' TIPS1 water molecules'
1539:          ELSE IF (ZSYM(NATOMS).EQ.'ME') THEN1539:          ELSE IF (ZSYM(NATOMS).EQ.'ME') THEN
1540:             WRITE(*,'(A,I4,A,I4,A)') ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' Mie atoms'1540:             WRITE(*,'(A,I4,A,I4,A)') ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' Mie atoms'
1595:      1          ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' CPMD atoms; bulk boundary conditions'1595:      1          ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,' CPMD atoms; bulk boundary conditions'
1596:          IF (CPMDC) WRITE(*,'(A,I6,A,I6,A)')1596:          IF (CPMDC) WRITE(*,'(A,I6,A,I6,A)')
1597:      1          ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,1597:      1          ' fetchz> ',NOPT,' Cartesian coordinates will be optimised for ',NATOMS,
1598:      1              ' CPMD atoms; cluster boundary conditions'1598:      1              ' CPMD atoms; cluster boundary conditions'
1599:          1599:          
1600:          IF (TWOD) WRITE(*,'(A)') ' fetchz> Two-dimensional flatland enforced'1600:          IF (TWOD) WRITE(*,'(A)') ' fetchz> Two-dimensional flatland enforced'
1601:          IF (DOUBLET) WRITE(*,'(A,F12.5,A,F12.5,A,F12.5)') 1601:          IF (DOUBLET) WRITE(*,'(A,F12.5,A,F12.5,A,F12.5)') 
1602:      1      ' fetchz> Double-well potential between first two atoms, barrier=',PARAM5,1602:      1      ' fetchz> Double-well potential between first two atoms, barrier=',PARAM5,
1603:      1                               ' first minimum at ',PARAM4,' second at ',PARAM4+2.0D0*PARAM61603:      1                               ' first minimum at ',PARAM4,' second at ',PARAM4+2.0D0*PARAM6
1604:       ENDIF1604:       ENDIF
1605:       IF (ZSYM(1)(1:1).EQ.'W') THEN ! must generalise for other rigid bodies somehow1605:       IF (ZSYM(1)(1:1).EQ.'W'.and.(ZSYM(1)(1:2).NE.'WC')) THEN ! must generalise for other rigid bodies somehow
1606:          DEALLOCATE (IATNUM,ATMASS)1606:          DEALLOCATE (IATNUM,ATMASS)
1607:          ALLOCATE (IATNUM(3*(NATOMS/2)),ATMASS(3*(NATOMS/2)))1607:          ALLOCATE (IATNUM(3*(NATOMS/2)),ATMASS(3*(NATOMS/2)))
1608:       ENDIF1608:       ENDIF
1609: C1609: C
1610: C DAE ATMASS is set in CHSETZSYMATMASS for CHARMM. PERTABLE would overwrite1610: C DAE ATMASS is set in CHSETZSYMATMASS for CHARMM. PERTABLE would overwrite
1611: C jmc similarly for unres in UNRSETZSYMATMASS...1611: C jmc similarly for unres in UNRSETZSYMATMASS...
1612: C1612: C
1613:       IF (.NOT.(CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.RINGPOLYMERT.OR.OPEPT)) CALL PERTABLE1613:       IF (.NOT.(CHRMMT.OR.UNRST.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.RINGPOLYMERT.OR.OPEPT)) CALL PERTABLE
1614:       IF (CASTEP) THEN1614:       IF (CASTEP) THEN
1615:          INQUIRE(FILE='castep.masses',EXIST=YESNO)1615:          INQUIRE(FILE='castep.masses',EXIST=YESNO)


r33624/gmetry.f 2018-01-19 09:34:34.210760772 +0000 r33623/gmetry.f 2018-01-19 09:34:39.150801701 +0000
 97:          ENDIF 97:          ENDIF
 98: 170   CONTINUE 98: 170   CONTINUE
 99:       If (IPRNT .GE. 20) then 99:       If (IPRNT .GE. 20) then
100:          WRITE (*,*) 'Cartesian coordinates from GMETRY.'100:          WRITE (*,*) 'Cartesian coordinates from GMETRY.'
101: C101: C
102:          WRITE (*,'(I3,3F12.6)') (i,Q(3*i-2),Q(3*i-1),Q(3*i),i=1,NATOMS)102:          WRITE (*,'(I3,3F12.6)') (i,Q(3*i-2),Q(3*i-1),Q(3*i),i=1,NATOMS)
103:       ENDIF103:       ENDIF
104: C104: C
105: C Check the theta angles if we are doing water.105: C Check the theta angles if we are doing water.
106: C106: C
107:       IF ((ZSYM(NATOMS)(1:1).EQ.'W').AND.(.NOT.ANGLEAXIS)) CALL TCHECK(VEC,Q)107:       IF (ZSYM(NATOMS)(1:1).EQ.'W'.AND.(.NOT.ANGLEAXIS)) then
 108:         IF (ZSYM(NATOMS)(1:2).NE.'WC') CALL TCHECK(VEC,Q)
 109:       ENDIF  
108:       RETURN110:       RETURN
109:       END111:       END


r33624/hybridmin.f 2018-01-19 09:34:34.474762903 +0000 r33623/hybridmin.f 2018-01-19 09:34:39.818808349 +0000
 50: C 50: C
 51:       INTEGER IWORK(33*3*NATOMS) 51:       INTEGER IWORK(33*3*NATOMS)
 52:       INTEGER ILWORK, LWORK, NFOUND, ISUPPZ(2*3*NATOMS) 52:       INTEGER ILWORK, LWORK, NFOUND, ISUPPZ(2*3*NATOMS)
 53:       INTEGER INFO, ISTAT, NEVSSAVE, NEVLSAVE 53:       INTEGER INFO, ISTAT, NEVSSAVE, NEVLSAVE
 54:       DOUBLE PRECISION WORK(33*3*NATOMS), ABSTOL, DIAG(3*NATOMS), DLAMCH, CEIGSAVE 54:       DOUBLE PRECISION WORK(33*3*NATOMS), ABSTOL, DIAG(3*NATOMS), DLAMCH, CEIGSAVE
 55:       LOGICAL KNOWE, KNOWG, KNOWH 55:       LOGICAL KNOWE, KNOWG, KNOWH
 56:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 56:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 57:  57: 
 58:       LWORK=33*3*NATOMS 58:       LWORK=33*3*NATOMS
 59:       ILWORK=33*3*NATOMS 59:       ILWORK=33*3*NATOMS
 60:       IF (ZSYM(1)(1:1).EQ.'W') THEN 60:       IF (ZSYM(1)(1:1).EQ.'W'.and.(ZSYM(1)(1:2).NE.'WC')) THEN
 61:          PRINT*,'HYBRIDMIN procedures have not been programmed for TIP potentials' 61:          PRINT*,'HYBRIDMIN procedures have not been programmed for TIP potentials'
 62:          CALL FLUSH(6) 62:          CALL FLUSH(6)
 63:          STOP 63:          STOP
 64:       ENDIF 64:       ENDIF
 65:       IF ((NUSEEV.GT.1).AND.(.NOT.NOIT.OR.NOHESS)) THEN 65:       IF ((NUSEEV.GT.1).AND.(.NOT.NOIT.OR.NOHESS)) THEN
 66:          WRITE(*,'(A)') 'For NUSEEV > 1 you must use a Hessian and NOIT with HYBRIDMIN' 66:          WRITE(*,'(A)') 'For NUSEEV > 1 you must use a Hessian and NOIT with HYBRIDMIN'
 67:          CALL FLUSH(6) 67:          CALL FLUSH(6)
 68:          STOP 68:          STOP
 69:       ENDIF 69:       ENDIF
 70: ! 70: !


r33624/keywords.f 2018-01-19 09:34:34.730764972 +0000 r33623/keywords.f 2018-01-19 09:34:40.066809548 +0000
  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/keywords.f' in revision 33624  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/keywords.f' in revision 33623


r33624/meccano.f90 2018-01-19 09:34:35.010767233 +0000 r33623/meccano.f90 2018-01-19 09:34:40.338811904 +0000
 46:    EINITIAL=E1 46:    EINITIAL=E1
 47:    EFINAL=E2 47:    EFINAL=E2
 48: ELSE 48: ELSE
 49:    CALL POTENTIAL(Q,EINITIAL,GRAD,.TRUE.,.FALSE.,RMSINITIAL,.FALSE.,.FALSE.) 49:    CALL POTENTIAL(Q,EINITIAL,GRAD,.TRUE.,.FALSE.,RMSINITIAL,.FALSE.,.FALSE.)
 50:    CALL POTENTIAL(FINAL,EFINAL,GRAD,.TRUE.,.FALSE.,RMSFINAL,.FALSE.,.FALSE.) 50:    CALL POTENTIAL(FINAL,EFINAL,GRAD,.TRUE.,.FALSE.,RMSFINAL,.FALSE.,.FALSE.)
 51: ENDIF 51: ENDIF
 52: IF (ITEST.AND.PTEST) WRITE(*,'(A,F20.10,A,F15.6,A,F20.10,A,F15.6)') & 52: IF (ITEST.AND.PTEST) WRITE(*,'(A,F20.10,A,F15.6,A,F20.10,A,F15.6)') &
 53:      &           ' Initial energy=',EINITIAL,' RMS force=',RMSINITIAL,' final energy=',EFINAL,' RMS=',RMSFINAL 53:      &           ' Initial energy=',EINITIAL,' RMS force=',RMSINITIAL,' final energy=',EFINAL,' RMS=',RMSFINAL
 54:  54: 
 55:  55: 
 56: IF (ZSYMSAVE(1:1).EQ.'W') THEN 56: IF (ZSYMSAVE(1:1).EQ.'W'.and.(ZSYMSAVE(1:2).NE.'WC')) THEN
 57:    PRINT '(A)','MECCANO not adapted for rigid bodies' 57:    PRINT '(A)','MECCANO not adapted for rigid bodies'
 58:    STOP 58:    STOP
 59: !  IADD=3*(NATOMS/2) 59: !  IADD=3*(NATOMS/2)
 60: !  LNOPT=NOPT/2      ! we only want to optimise the angular degrees of freedom for TIP 60: !  LNOPT=NOPT/2      ! we only want to optimise the angular degrees of freedom for TIP
 61: !  LNATOMS=NATOMS/2  ! for TIP potentials both centre-of-mass and Euler angles are now in odata 61: !  LNATOMS=NATOMS/2  ! for TIP potentials both centre-of-mass and Euler angles are now in odata
 62: ELSE 62: ELSE
 63:    IADD=0 63:    IADD=0
 64:    LNOPT=NOPT 64:    LNOPT=NOPT
 65:    LNATOMS=NATOMS 65:    LNATOMS=NATOMS
 66: ENDIF 66: ENDIF


r33624/mindist.f 2018-01-19 09:34:35.266769322 +0000 r33623/mindist.f 2018-01-19 09:34:40.606814224 +0000
 56: C 56: C
 57:       DO J1=1,3 57:       DO J1=1,3
 58:          DO J2=1,3 58:          DO J2=1,3
 59:             OMEGATOT(J2,J1)=0.0D0 59:             OMEGATOT(J2,J1)=0.0D0
 60:          ENDDO 60:          ENDDO
 61:          OMEGATOT(J1,J1)=1.0D0 61:          OMEGATOT(J1,J1)=1.0D0
 62:       ENDDO 62:       ENDDO
 63:  63: 
 64:       AGAIN=.FALSE. 64:       AGAIN=.FALSE.
 65:       DSAVE=0.0D0 65:       DSAVE=0.0D0
 66:       IF (ZUSE(1:1).EQ.'W') THEN 66:       IF (ZUSE(1:1).EQ.'W'.and.ZUSE(1:2).NE.'WC') THEN
 67:          ALLOCATE(R(3,3*(NATOMS/2)),R0(3,3*(NATOMS/2)),R1(3,3*(NATOMS/2)),R1SAVE(3,3*(NATOMS/2))) 67:          ALLOCATE(R(3,3*(NATOMS/2)),R0(3,3*(NATOMS/2)),R1(3,3*(NATOMS/2)),R1SAVE(3,3*(NATOMS/2)))
 68: C        NSIZE=3*NATOMS 68: C        NSIZE=3*NATOMS
 69:          NSIZE=3*(NATOMS/2) 69:          NSIZE=3*(NATOMS/2)
 70: C        DO J1=1,NATOMS 70: C        DO J1=1,NATOMS
 71:          DO J1=1,NATOMS/2 71:          DO J1=1,NATOMS/2
 72:             CALL CONVERT(RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3), 72:             CALL CONVERT(RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3),
 73: C    1              RA(3*(NATOMS+J1-1)+1),RA(3*(NATOMS+J1-1)+2),RA(3*(NATOMS+J1-1)+3),OVEC,H1VEC,H2VEC) 73: C    1              RA(3*(NATOMS+J1-1)+1),RA(3*(NATOMS+J1-1)+2),RA(3*(NATOMS+J1-1)+3),OVEC,H1VEC,H2VEC)
 74:      1              RA(3*(NATOMS/2+J1-1)+1),RA(3*(NATOMS/2+J1-1)+2),RA(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC) 74:      1              RA(3*(NATOMS/2+J1-1)+1),RA(3*(NATOMS/2+J1-1)+2),RA(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC)
 75:             R(1,(J1-1)*3+1)=OVEC(1) 75:             R(1,(J1-1)*3+1)=OVEC(1)
 76:             R(2,(J1-1)*3+1)=OVEC(2) 76:             R(2,(J1-1)*3+1)=OVEC(2)
203: C              DO J2=1,3203: C              DO J2=1,3
204: C                 OSAVE(J2,J1)=0.0D0204: C                 OSAVE(J2,J1)=0.0D0
205: C                 OSAVESAVE(J2,J1)=OSAVE(J2,J1)205: C                 OSAVESAVE(J2,J1)=OSAVE(J2,J1)
206: C              ENDDO206: C              ENDDO
207: C              OSAVE(J1,J1)=1.0D0207: C              OSAVE(J1,J1)=1.0D0
208: C           ENDDO208: C           ENDDO
209: C           GOTO 10209: C           GOTO 10
210: C        ENDIF210: C        ENDIF
211: C     ENDIF211: C     ENDIF
212: 212: 
213:       IF (ZUSE(1:1).EQ.'W') THEN213:       IF (ZUSE(1:1).EQ.'W'.and.ZUSE(1:2).NE.'WC') THEN
214: C        DO J1=1,NATOMS  !  WCOMMENT214: C        DO J1=1,NATOMS  !  WCOMMENT
215:          DO J1=1,NATOMS/2215:          DO J1=1,NATOMS/2
216:             OVEC(1)=R(1,(J1-1)*3+1)216:             OVEC(1)=R(1,(J1-1)*3+1)
217:             OVEC(2)=R(2,(J1-1)*3+1)217:             OVEC(2)=R(2,(J1-1)*3+1)
218:             OVEC(3)=R(3,(J1-1)*3+1)218:             OVEC(3)=R(3,(J1-1)*3+1)
219:             H1VEC(1)=R(1,(J1-1)*3+2)219:             H1VEC(1)=R(1,(J1-1)*3+2)
220:             H1VEC(2)=R(2,(J1-1)*3+2)220:             H1VEC(2)=R(2,(J1-1)*3+2)
221:             H1VEC(3)=R(3,(J1-1)*3+2)221:             H1VEC(3)=R(3,(J1-1)*3+2)
222:             H2VEC(1)=R(1,(J1-1)*3+3)222:             H2VEC(1)=R(1,(J1-1)*3+3)
223:             H2VEC(2)=R(2,(J1-1)*3+3)223:             H2VEC(2)=R(2,(J1-1)*3+3)


r33624/morph.f 2018-01-19 09:34:35.494771184 +0000 r33623/morph.f 2018-01-19 09:34:40.854816369 +0000
 40:       LOGICAL MFLAG, PVFLAG, RESET, POTCALL, PTEST, PUSH, PULL, FINISHED 40:       LOGICAL MFLAG, PVFLAG, RESET, POTCALL, PTEST, PUSH, PULL, FINISHED
 41:       COMMON /PVF/ PVFLAG 41:       COMMON /PVF/ PVFLAG
 42:       CHARACTER(LEN=5) ZSYMSAVE 42:       CHARACTER(LEN=5) ZSYMSAVE
 43:       COMMON /SYS/ ZSYMSAVE 43:       COMMON /SYS/ ZSYMSAVE
 44:       LOGICAL KNOWE, KNOWG, KNOWH 44:       LOGICAL KNOWE, KNOWG, KNOWH
 45:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 45:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 46:       COMMON /MORPHDATA/ PUSH, PULL 46:       COMMON /MORPHDATA/ PUSH, PULL
 47:  47: 
 48:       IF (DEBUG) OPEN(UNIT=995,FILE='morph.xyz',STATUS='UNKNOWN') 48:       IF (DEBUG) OPEN(UNIT=995,FILE='morph.xyz',STATUS='UNKNOWN')
 49:       CALL MYCPU_TIME(STARTTIME,.TRUE.) 49:       CALL MYCPU_TIME(STARTTIME,.TRUE.)
 50:       IF (ZSYM(1)(1:1).EQ.'W') THEN 50:       IF (ZSYM(1)(1:1).EQ.'W'.and.ZSYM(1)(1:2).NE.'WC') THEN
 51:          PRINT*,'MORPH procedures have not been programmed for TIP potentials' 51:          PRINT*,'MORPH procedures have not been programmed for TIP potentials'
 52:          STOP 52:          STOP
 53:       ENDIF 53:       ENDIF
 54:       PRINT '(A,F8.2,A,F12.2,A,F10.2)',' morph> Interpolation using maximum step ',MORPHMXSTP,' energy maximum ', 54:       PRINT '(A,F8.2,A,F12.2,A,F10.2)',' morph> Interpolation using maximum step ',MORPHMXSTP,' energy maximum ',
 55:      &                                   MORPHEMAX,' maximum rise ',MORPHERISE 55:      &                                   MORPHEMAX,' maximum rise ',MORPHERISE
 56:       PULL=.TRUE. 56:       PULL=.TRUE.
 57:       PUSH=.FALSE. 57:       PUSH=.FALSE.
 58:       NPU=0 58:       NPU=0
 59:       FINISHED=.FALSE. 59:       FINISHED=.FALSE.
 60:       MNBFGSMAX1SAVE=MNBFGSMAX1 60:       MNBFGSMAX1SAVE=MNBFGSMAX1


r33624/newmindist.f90 2018-01-19 09:34:35.758773341 +0000 r33623/newmindist.f90 2018-01-19 09:34:41.110818600 +0000
 84:  84: 
 85: ! WRITE(*,*) NATOMS 85: ! WRITE(*,*) NATOMS
 86: ! WRITE(*,*) 'RA starting geometry' 86: ! WRITE(*,*) 'RA starting geometry'
 87: ! WRITE(*,'(A3,3G20.10)') ('LA ',RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3),J1=1,NATOMS) 87: ! WRITE(*,'(A3,3G20.10)') ('LA ',RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3),J1=1,NATOMS)
 88: ! WRITE(*,*) NATOMS 88: ! WRITE(*,*) NATOMS
 89: ! WRITE(*,*) 'RB starting geometry' 89: ! WRITE(*,*) 'RB starting geometry'
 90: ! WRITE(*,'(A3,3G20.10)') ('LA ',RB(3*(J1-1)+1),RB(3*(J1-1)+2),RB(3*(J1-1)+3),J1=1,NATOMS) 90: ! WRITE(*,'(A3,3G20.10)') ('LA ',RB(3*(J1-1)+1),RB(3*(J1-1)+2),RB(3*(J1-1)+3),J1=1,NATOMS)
 91: !  91: ! 
 92: ! Convert rigid body coordinates to Cartesians for rigid bodies.  92: ! Convert rigid body coordinates to Cartesians for rigid bodies. 
 93: ! 93: !
 94: IF (ZUSE(1:1).EQ.'W') THEN 94: IF (ZUSE(1:1).EQ.'W'.and.ZUSE(1:2).NE.'WC') THEN
 95:    ALLOCATE(XA(3*3*(NATOMS/2)),XB(3*3*(NATOMS/2))) 95:    ALLOCATE(XA(3*3*(NATOMS/2)),XB(3*3*(NATOMS/2)))
 96:    NSIZE=3*(NATOMS/2) 96:    NSIZE=3*(NATOMS/2)
 97:    DO J1=1,NATOMS/2 97:    DO J1=1,NATOMS/2
 98:       CALL CONVERT(RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3), & 98:       CALL CONVERT(RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3), &
 99:      &        RA(3*(NATOMS/2+J1-1)+1),RA(3*(NATOMS/2+J1-1)+2),RA(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC) 99:      &        RA(3*(NATOMS/2+J1-1)+1),RA(3*(NATOMS/2+J1-1)+2),RA(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC)
100:       XA(9*(J1-1)+0+1)=OVEC(1)100:       XA(9*(J1-1)+0+1)=OVEC(1)
101:       XA(9*(J1-1)+0+2)=OVEC(2)101:       XA(9*(J1-1)+0+2)=OVEC(2)
102:       XA(9*(J1-1)+0+3)=OVEC(3)102:       XA(9*(J1-1)+0+3)=OVEC(3)
103:       XA(9*(J1-1)+3+1)=H1VEC(1)103:       XA(9*(J1-1)+3+1)=H1VEC(1)
104:       XA(9*(J1-1)+3+2)=H1VEC(2)104:       XA(9*(J1-1)+3+2)=H1VEC(2)
359:    RMAT(1,2)=2*(Q2*Q3+Q1*Q4)359:    RMAT(1,2)=2*(Q2*Q3+Q1*Q4)
360:    RMAT(1,3)=2*(Q2*Q4-Q1*Q3)360:    RMAT(1,3)=2*(Q2*Q4-Q1*Q3)
361:    RMAT(2,1)=2*(Q2*Q3-Q1*Q4)361:    RMAT(2,1)=2*(Q2*Q3-Q1*Q4)
362:    RMAT(2,2)=Q1**2+Q3**2-Q2**2-Q4**2362:    RMAT(2,2)=Q1**2+Q3**2-Q2**2-Q4**2
363:    RMAT(2,3)=2*(Q3*Q4+Q1*Q2)363:    RMAT(2,3)=2*(Q3*Q4+Q1*Q2)
364:    RMAT(3,1)=2*(Q2*Q4+Q1*Q3)364:    RMAT(3,1)=2*(Q2*Q4+Q1*Q3)
365:    RMAT(3,2)=2*(Q3*Q4-Q1*Q2)365:    RMAT(3,2)=2*(Q3*Q4-Q1*Q2)
366:    RMAT(3,3)=Q1**2+Q4**2-Q2**2-Q3**2366:    RMAT(3,3)=Q1**2+Q4**2-Q2**2-Q3**2
367: ENDIF367: ENDIF
368: IF (.NOT.(PRESERVET.OR.MIEFT.OR.NOTRANSROTT)) THEN ! never translate with substrate!368: IF (.NOT.(PRESERVET.OR.MIEFT.OR.NOTRANSROTT)) THEN ! never translate with substrate!
369:    IF (ZUSE(1:1).EQ.'W') THEN369:    IF (ZUSE(1:1).EQ.'W'.and.ZUSE(1:2).NE.'WC') THEN
370: !370: !
371: !  Translate the XB coordinates to the centre of coordinates of XA.371: !  Translate the XB coordinates to the centre of coordinates of XA.
372: !372: !
373:       DO J1=1,NSIZE373:       DO J1=1,NSIZE
374:          XB(3*(J1-1)+1)=XB(3*(J1-1)+1)+CMXA+XSHIFT374:          XB(3*(J1-1)+1)=XB(3*(J1-1)+1)+CMXA+XSHIFT
375:          XB(3*(J1-1)+2)=XB(3*(J1-1)+2)+CMYA+YSHIFT375:          XB(3*(J1-1)+2)=XB(3*(J1-1)+2)+CMYA+YSHIFT
376:          XB(3*(J1-1)+3)=XB(3*(J1-1)+3)+CMZA+ZSHIFT376:          XB(3*(J1-1)+3)=XB(3*(J1-1)+3)+CMZA+ZSHIFT
377:       ENDDO377:       ENDDO
378: !378: !
379: !  Rotate XB coordinates about new centre of mass379: !  Rotate XB coordinates about new centre of mass


r33624/oldneb.f 2018-01-19 09:34:36.058775800 +0000 r33623/oldneb.f 2018-01-19 09:34:41.370820874 +0000
 37:       CHARACTER FNAME*80 37:       CHARACTER FNAME*80
 38:       INTEGER IADD, LNOPT, LNATOMS 38:       INTEGER IADD, LNOPT, LNATOMS
 39:       COMMON /INTS/ IADD, LNOPT, LNATOMS 39:       COMMON /INTS/ IADD, LNOPT, LNATOMS
 40:       CHARACTER(LEN=5) ZSYMSAVE 40:       CHARACTER(LEN=5) ZSYMSAVE
 41:       COMMON /SYS/ ZSYMSAVE 41:       COMMON /SYS/ ZSYMSAVE
 42:  42: 
 43: ! 43: !
 44: ! The declaration of POINTS provides more space than is actually needed for TIP potential. 44: ! The declaration of POINTS provides more space than is actually needed for TIP potential.
 45: ! 9*NATOMS*NIMAGEMAX/2 is needed in reality. 45: ! 9*NATOMS*NIMAGEMAX/2 is needed in reality.
 46: ! 46: !
 47:       IF (ZSYMSAVE(1:1).EQ.'W') THEN 47:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
 48:          IADD=3*(NATOMS/2) 48:          IADD=3*(NATOMS/2)
 49:          LNOPT=NOPT/2      ! we only want to optimise the angular degrees of freedom for TIP 49:          LNOPT=NOPT/2      ! we only want to optimise the angular degrees of freedom for TIP
 50:          LNATOMS=NATOMS/2  ! for TIP potentials both centre-of-mass and Euler angles are now in odata 50:          LNATOMS=NATOMS/2  ! for TIP potentials both centre-of-mass and Euler angles are now in odata
 51:       ELSE 51:       ELSE
 52:          IADD=0 52:          IADD=0
 53:          LNOPT=NOPT 53:          LNOPT=NOPT
 54:          LNATOMS=NATOMS 54:          LNATOMS=NATOMS
 55:       ENDIF 55:       ENDIF
 56:       NMAGNIFY=0 56:       NMAGNIFY=0
 57:       IF (NIMAGE.GT.NIMAGEMAX) THEN 57:       IF (NIMAGE.GT.NIMAGEMAX) THEN
184:       CLOSE(4)184:       CLOSE(4)
185: C185: C
186: C  Dump pathway as an xyz file186: C  Dump pathway as an xyz file
187: C187: C
188:       IF (FILTH.EQ.0) THEN188:       IF (FILTH.EQ.0) THEN
189:          WRITE(FNAME,'(A12)') 'neb.path.xyz'189:          WRITE(FNAME,'(A12)') 'neb.path.xyz'
190:       ELSE 190:       ELSE 
191:          WRITE(FNAME,'(A)') 'neb.path.xyz.'//TRIM(ADJUSTL(FILTHSTR))191:          WRITE(FNAME,'(A)') 'neb.path.xyz.'//TRIM(ADJUSTL(FILTHSTR))
192:       ENDIF192:       ENDIF
193:       OPEN(UNIT=3,FILE=FNAME,STATUS='UNKNOWN')193:       OPEN(UNIT=3,FILE=FNAME,STATUS='UNKNOWN')
194:       IF (ZSYM(1)(1:1).EQ.'W') THEN194:       IF (ZSYM(1)(1:1).EQ.'W'.and.ZSYM(1)(1:2).NE.'WC') THEN
195:          WRITE(3,'(I6)') 3*LNATOMS195:          WRITE(3,'(I6)') 3*LNATOMS
196:          WRITE(3,'(F20.10)') EINITIAL196:          WRITE(3,'(F20.10)') EINITIAL
197:          DO J1=1,LNATOMS197:          DO J1=1,LNATOMS
198:             CALL CONVERT(Q(3*(J1-1)+1),Q(3*(J1-1)+2),Q(3*(J1-1)+3),198:             CALL CONVERT(Q(3*(J1-1)+1),Q(3*(J1-1)+2),Q(3*(J1-1)+3),
199:      1       Q(3*(LNATOMS+J1-1)+1),Q(3*(LNATOMS+J1-1)+2),Q(3*(LNATOMS+J1-1)+3),OVEC,H1VEC,H2VEC)199:      1       Q(3*(LNATOMS+J1-1)+1),Q(3*(LNATOMS+J1-1)+2),Q(3*(LNATOMS+J1-1)+3),OVEC,H1VEC,H2VEC)
200:             WRITE(3,'(A2,4X,3F20.10)') 'O  ',OVEC(1),OVEC(2),OVEC(3)200:             WRITE(3,'(A2,4X,3F20.10)') 'O  ',OVEC(1),OVEC(2),OVEC(3)
201:             WRITE(3,'(A2,4X,3F20.10)') 'H  ',H1VEC(1),H1VEC(2),H1VEC(3)201:             WRITE(3,'(A2,4X,3F20.10)') 'H  ',H1VEC(1),H1VEC(2),H1VEC(3)
202:             WRITE(3,'(A2,4X,3F20.10)') 'H  ',H2VEC(1),H2VEC(2),H2VEC(3)202:             WRITE(3,'(A2,4X,3F20.10)') 'H  ',H2VEC(1),H2VEC(2),H2VEC(3)
203:          ENDDO203:          ENDDO
204:          DO J1=1,NIMAGE204:          DO J1=1,NIMAGE
1091:       COMMON /INTS/ IADD, LNOPT, LNATOMS1091:       COMMON /INTS/ IADD, LNOPT, LNATOMS
1092:       INTEGER J1,J2,J31092:       INTEGER J1,J2,J3
1093:       DOUBLE PRECISION VECS(3*NATOMS), STEP, GUESS(3*NATOMS),DPRAND,TEMP(3),Q(3*NATOMS),1093:       DOUBLE PRECISION VECS(3*NATOMS), STEP, GUESS(3*NATOMS),DPRAND,TEMP(3),Q(3*NATOMS),
1094:      1                 POINTS(5*NATOMS*NIMAGEMAX), OVEC1(3),H1VEC1(3),H2VEC1(3), C2VEC(3),1094:      1                 POINTS(5*NATOMS*NIMAGEMAX), OVEC1(3),H1VEC1(3),H2VEC1(3), C2VEC(3),
1095:      2                 OVEC2(3),H1VEC2(3),H2VEC2(3),1095:      2                 OVEC2(3),H1VEC2(3),H2VEC2(3),
1096:      3                 HSVEL(3*NATOMS),1096:      3                 HSVEL(3*NATOMS),
1097:      4                 DUMMY, EFINAL, EINITIAL,SEPARATION,RMS,DIST1097:      4                 DUMMY, EFINAL, EINITIAL,SEPARATION,RMS,DIST
1098:       COMMON /NEBRMS/ RMS,EINITIAL,EFINAL,SEPARATION1098:       COMMON /NEBRMS/ RMS,EINITIAL,EFINAL,SEPARATION
1099:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: FINCART,STARTCART,POINTSCART1099:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: FINCART,STARTCART,POINTSCART
1100: 1100: 
1101:       IF (ZSYM(1)(1:1).EQ.'W') THEN1101:       IF (ZSYM(1)(1:1).EQ.'W'.and.ZSYM(1)(1:2).NE.'WC') THEN
1102:          ALLOCATE(FINCART(9*LNATOMS),STARTCART(9*LNATOMS),POINTSCART(9*LNATOMS))1102:          ALLOCATE(FINCART(9*LNATOMS),STARTCART(9*LNATOMS),POINTSCART(9*LNATOMS))
1103:       ELSE1103:       ELSE
1104:          ALLOCATE(FINCART(3*NATOMS),STARTCART(3*NATOMS),POINTSCART(3*NATOMS))1104:          ALLOCATE(FINCART(3*NATOMS),STARTCART(3*NATOMS),POINTSCART(3*NATOMS))
1105:       ENDIF1105:       ENDIF
1106:       IF (.NOT.ALLOCATED(LANV)) ALLOCATE(LANV(3*LNATOMS+IADD))1106:       IF (.NOT.ALLOCATED(LANV)) ALLOCATE(LANV(3*LNATOMS+IADD))
1107: C1107: C
1108: C  Bulk principal image code assumes that PARAM1, PARAM2 and PARAM3 are the box lengths.1108: C  Bulk principal image code assumes that PARAM1, PARAM2 and PARAM3 are the box lengths.
1109: C1109: C
1110:       DO J1=1,LNATOMS+IADD/31110:       DO J1=1,LNATOMS+IADD/3
1111:          J3=3*(J1-1)1111:          J3=3*(J1-1)
1164:             DIST=DIST+(FIN(J1)-START(J1)-PARAM3*LANV(J1))**21164:             DIST=DIST+(FIN(J1)-START(J1)-PARAM3*LANV(J1))**2
1165:          ENDDO1165:          ENDDO
1166:       ENDIF1166:       ENDIF
1167:       DIST=SQRT(DIST)1167:       DIST=SQRT(DIST)
1168: C     PRINT *,'oldneb> DIST=',DIST1168: C     PRINT *,'oldneb> DIST=',DIST
1169: C1169: C
1170: C  Assign initial image positions.1170: C  Assign initial image positions.
1171: C1171: C
1172:       SEPARATION=DIST/(NIMAGE+1)1172:       SEPARATION=DIST/(NIMAGE+1)
1173: C     NEBK=(ABS(EFINAL)+ABS(EINITIAL))/(100*(NIMAGE+1)*SEPARATION**2)1173: C     NEBK=(ABS(EFINAL)+ABS(EINITIAL))/(100*(NIMAGE+1)*SEPARATION**2)
1174:       IF (ZSYM(1)(1:1).EQ.'W') NEBK=NEBK*0.0D01174:       IF (ZSYM(1)(1:1).EQ.'W'.and.ZSYM(1)(1:2).NE.'WC') NEBK=NEBK*0.0D0
1175:       IF (ZSYM(1)(1:1).EQ.'W') THEN1175:       IF (ZSYM(1)(1:1).EQ.'W'.and.ZSYM(1)(1:2).NE.'WC') THEN
1176: C1176: C
1177: C  Note that interpolating the Cartesian coordinates will generally break the geometry of the individual water molecules,1177: C  Note that interpolating the Cartesian coordinates will generally break the geometry of the individual water molecules,
1178: C  so there won't be a perfect conversion back to centre-of-mass/Euler angle coordinates.1178: C  so there won't be a perfect conversion back to centre-of-mass/Euler angle coordinates.
1179: C1179: C
1180:          DO J1=1,NIMAGE1180:          DO J1=1,NIMAGE
1181:             DO J2=1,LNATOMS1181:             DO J2=1,LNATOMS
1182:                POINTSCART(9*(J2-1)+1)=STARTCART(9*(J2-1)+1)+(FINCART(9*(J2-1)+1)-STARTCART(9*(J2-1)+1))*J1*SEPARATION/DIST1182:                POINTSCART(9*(J2-1)+1)=STARTCART(9*(J2-1)+1)+(FINCART(9*(J2-1)+1)-STARTCART(9*(J2-1)+1))*J1*SEPARATION/DIST
1183:                POINTSCART(9*(J2-1)+2)=STARTCART(9*(J2-1)+2)+(FINCART(9*(J2-1)+2)-STARTCART(9*(J2-1)+2))*J1*SEPARATION/DIST1183:                POINTSCART(9*(J2-1)+2)=STARTCART(9*(J2-1)+2)+(FINCART(9*(J2-1)+2)-STARTCART(9*(J2-1)+2))*J1*SEPARATION/DIST
1184:                POINTSCART(9*(J2-1)+3)=STARTCART(9*(J2-1)+3)+(FINCART(9*(J2-1)+3)-STARTCART(9*(J2-1)+3))*J1*SEPARATION/DIST1184:                POINTSCART(9*(J2-1)+3)=STARTCART(9*(J2-1)+3)+(FINCART(9*(J2-1)+3)-STARTCART(9*(J2-1)+3))*J1*SEPARATION/DIST
1185:                POINTSCART(9*(J2-1)+4)=STARTCART(9*(J2-1)+4)+(FINCART(9*(J2-1)+4)-STARTCART(9*(J2-1)+4))*J1*SEPARATION/DIST1185:                POINTSCART(9*(J2-1)+4)=STARTCART(9*(J2-1)+4)+(FINCART(9*(J2-1)+4)-STARTCART(9*(J2-1)+4))*J1*SEPARATION/DIST


r33624/OPTIM.F 2018-01-19 09:34:32.618748094 +0000 r33623/OPTIM.F 2018-01-19 09:34:37.586788491 +0000
167: 167: 
168:       IF (AMBERT .OR. NABT) CALL SET_CHECK_CHIRAL(Q,NATOMS,GOODSTRUCTURE1,CHIARRAY1)168:       IF (AMBERT .OR. NABT) CALL SET_CHECK_CHIRAL(Q,NATOMS,GOODSTRUCTURE1,CHIARRAY1)
169: 169: 
170:       IF (CHECKDT) CALL CHECKDRVTS(Q)170:       IF (CHECKDT) CALL CHECKDRVTS(Q)
171: 171: 
172:       IF ((INTCONSTRAINTT.OR.INTLJT).AND.(NCONGEOM.GE.2)) THEN172:       IF ((INTCONSTRAINTT.OR.INTLJT).AND.(NCONGEOM.GE.2)) THEN
173: !173: !
174: ! Set up all the constraints and repulsions for zero frozen atoms.174: ! Set up all the constraints and repulsions for zero frozen atoms.
175: !175: !
176:          IF (.NOT.ALLOCATED(CONI)) THEN176:          IF (.NOT.ALLOCATED(CONI)) THEN
177:             ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),177:             ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX))
178:      &               CONOFFTRIED(INTCONMAX)) 
179:             CONOFFTRIED(1:INTCONMAX)=.FALSE. 
180:             ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))178:             ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
181:          ENDIF179:          ENDIF
182:          INTFREEZETOLSAVE=INTFREEZETOL180:          INTFREEZETOLSAVE=INTFREEZETOL
183:          INTFREEZETOL=-1.0D0181:          INTFREEZETOL=-1.0D0
184:          PRINT *,'OPTIM> NCONGEOM=',NCONGEOM182:          PRINT *,'OPTIM> NCONGEOM=',NCONGEOM
185:          CALL MAKE_CONPOT(NCONGEOM,CONGEOM) 183:          CALL MAKE_CONPOT(NCONGEOM,CONGEOM) 
186:          INTFREEZETOL=INTFREEZETOLSAVE184:          INTFREEZETOL=INTFREEZETOLSAVE
187:       ENDIF185:       ENDIF
188: 186: 
189: 187: 
222:       IF ((ZSYMSAVE.EQ.'RM').AND.(PARAM1.NE.0.0D0)) BULKT=.TRUE.220:       IF ((ZSYMSAVE.EQ.'RM').AND.(PARAM1.NE.0.0D0)) BULKT=.TRUE.
223: C     LDUMMY=(NATOMS.EQ.64).AND.(ZSYMSAVE(1:1).EQ.'W') 221: C     LDUMMY=(NATOMS.EQ.64).AND.(ZSYMSAVE(1:1).EQ.'W') 
224:       BULKT=(BULKT.OR.222:       BULKT=(BULKT.OR.
225:      1      (ZSYMSAVE.EQ.'ME').OR.(ZSYMSAVE.EQ.'P6').OR.223:      1      (ZSYMSAVE.EQ.'ME').OR.(ZSYMSAVE.EQ.'P6').OR.
226:      2      (ZSYMSAVE.EQ.'SC').OR.(ZSYMSAVE.EQ.'MS').OR.224:      2      (ZSYMSAVE.EQ.'SC').OR.(ZSYMSAVE.EQ.'MS').OR.
227:      3      (ZSYMSAVE.EQ.'MP').OR.(ZSYMSAVE.EQ.'JM').OR.225:      3      (ZSYMSAVE.EQ.'MP').OR.(ZSYMSAVE.EQ.'JM').OR.
228:      3      (ZSYMSAVE.EQ.'DS').OR.226:      3      (ZSYMSAVE.EQ.'DS').OR.
229:      4      (ZSYMSAVE.EQ.'LP').OR.(ZSYMSAVE.EQ.'LS').OR.(ZSYMSAVE.EQ.'LC').OR.227:      4      (ZSYMSAVE.EQ.'LP').OR.(ZSYMSAVE.EQ.'LS').OR.(ZSYMSAVE.EQ.'LC').OR.
230:      5      (ZSYMSAVE.EQ.'LK'))228:      5      (ZSYMSAVE.EQ.'LK'))
231: 229: 
232:       IF (BULKT.AND.NEWCONNECTT.AND.(.NOT.QUIPT)) THEN230:       IF (BULKT.AND.NEWCONNECTT) THEN
233:          CALL INITIALBOX(Q)231:          CALL INITIALBOX(Q)
234:          CALL INITIALBOX(FIN)232:          CALL INITIALBOX(FIN)
235:       ENDIF233:       ENDIF
236: 234: 
237: C     NOSHIFT=(FIELDT.OR.BFGSMINT.OR.BSMIN.OR.RKMIN.OR.NOHESS.OR.BFGSSTEP)235: C     NOSHIFT=(FIELDT.OR.BFGSMINT.OR.BSMIN.OR.RKMIN.OR.NOHESS.OR.BFGSSTEP)
238: C     IF (INR.GE.0) NOSHIFT=.FALSE. ! changed INR default value in keywords to -1236: C     IF (INR.GE.0) NOSHIFT=.FALSE. ! changed INR default value in keywords to -1
239: C                                     This isn;t good enough - if we are doing a path and237: C                                     This isn;t good enough - if we are doing a path and
240: C                                     EFOL is called once to get the eigenvalues and eigenvectors238: C                                     EFOL is called once to get the eigenvalues and eigenvectors
241: C                                     then we need to SHIFT.239: C                                     then we need to SHIFT.
242: C  Shifting only needs to be turned off if we are doing VARIABLES or BFGSTS without NOIT;240: C  Shifting only needs to be turned off if we are doing VARIABLES or BFGSTS without NOIT;
250: C248: C
251: C  Jump back here for MULTIJOB runs.249: C  Jump back here for MULTIJOB runs.
252: C250: C
253:       MULTIINR=INR251:       MULTIINR=INR
254: 963   CONTINUE252: 963   CONTINUE
255: C253: C
256: C  Resize the system if required.254: C  Resize the system if required.
257: C255: C
258:       IF (RESIZE.NE.1.0D0) THEN256:       IF (RESIZE.NE.1.0D0) THEN
259:          PRINT*,'Scaling coordinates by ',RESIZE257:          PRINT*,'Scaling coordinates by ',RESIZE
260:          IF (ZSYMSAVE(1:1).EQ.'W') THEN258:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
261:             DO J1=1,3*(NATOMS/2)259:             DO J1=1,3*(NATOMS/2)
262:                Q(J1)=Q(J1)*RESIZE260:                Q(J1)=Q(J1)*RESIZE
263:             ENDDO261:             ENDDO
264:          ELSE262:          ELSE
265:             DO J1=1,NOPT263:             DO J1=1,NOPT
266:                Q(J1)=Q(J1)*RESIZE264:                Q(J1)=Q(J1)*RESIZE
267:             ENDDO265:             ENDDO
268:          ENDIF266:          ENDIF
269:       ENDIF267:       ENDIF
270:       IF (CASTEP.AND.PRESSURE) CALL CLATMIN(Q,VNEW)268:       IF (CASTEP.AND.PRESSURE) CALL CLATMIN(Q,VNEW)
290: C     DO J1=1,NATOMS288: C     DO J1=1,NATOMS
291: C        IF (IATNUM(J1).NE.0) J2=J2+1289: C        IF (IATNUM(J1).NE.0) J2=J2+1
292: C     ENDDO290: C     ENDDO
293: C     NREAL=J2291: C     NREAL=J2
294: C292: C
295: C     CALL GMETRY(0,VEC,Q)293: C     CALL GMETRY(0,VEC,Q)
296: C294: C
297: C  Don't call symmetry if we're doing Fenske-Hall295: C  Don't call symmetry if we're doing Fenske-Hall
298: C296: C
299:       IF ((ZSYMSAVE.NE.'FH').AND.(.NOT.VARIABLES).AND.(.NOT.AMBERT).AND.(.NOT.CHRMMT).AND.297:       IF ((ZSYMSAVE.NE.'FH').AND.(.NOT.VARIABLES).AND.(.NOT.AMBERT).AND.(.NOT.CHRMMT).AND.
300:      1 (.NOT.NABT).AND.(.NOT.UNRST).AND.(.NOT.RINGPOLYMERT).AND.(.NOT.AMBER12T).AND.(.NOT.BULKT)) THEN298:      1 (.NOT.NABT).AND.(.NOT.UNRST).AND.(.NOT.RINGPOLYMERT).AND.(.NOT.AMBER12T)) THEN
301: C299: C
302: C  For W1/W2/W3/W4 potentials Q contains the centre of mass coordinates300: C  For W1/W2/W3/W4 potentials Q contains the centre of mass coordinates
303: C  followed by the Euler angles and the number of molecules is NATOMS/2301: C  followed by the Euler angles and the number of molecules is NATOMS/2
304: C302: C
305:          IF (ZSYMSAVE(1:1).EQ.'W') THEN303:          IF (ZSYMSAVE(1:1).EQ.'W' .AND. ZSYMSAVE(1:2) .NE. 'WC') THEN
306:             ALLOCATE(QW(9*(NATOMS/2)))304:             ALLOCATE(QW(9*(NATOMS/2)))
307: C           DO J2=1,NOPT+IADD305: C           DO J2=1,NOPT+IADD
308: C              QSAVE(J2)=Q(J2)306: C              QSAVE(J2)=Q(J2)
309: C           ENDDO307: C           ENDDO
310: C           DO J2=1,NATOMS  !  WCOMMENT308: C           DO J2=1,NATOMS  !  WCOMMENT
311:             DO J2=1,NATOMS/2309:             DO J2=1,NATOMS/2
312:                CALL CONVERT(Q(3*(J2-1)+1),Q(3*(J2-1)+2),Q(3*(J2-1)+3),310:                CALL CONVERT(Q(3*(J2-1)+1),Q(3*(J2-1)+2),Q(3*(J2-1)+3),
313: C    1                      Q(3*(NATOMS+J2-1)+1),Q(3*(NATOMS+J2-1)+2),Q(3*(NATOMS+J2-1)+3),311: C    1                      Q(3*(NATOMS+J2-1)+1),Q(3*(NATOMS+J2-1)+2),Q(3*(NATOMS+J2-1)+3),
314:      1                      Q(3*(NATOMS/2+J2-1)+1),Q(3*(NATOMS/2+J2-1)+2),Q(3*(NATOMS/2+J2-1)+3),312:      1                      Q(3*(NATOMS/2+J2-1)+1),Q(3*(NATOMS/2+J2-1)+2),Q(3*(NATOMS/2+J2-1)+3),
315:      2                      OVEC,H1VEC,H2VEC)313:      2                      OVEC,H1VEC,H2VEC)
499: 497: 
500: 555   CONTINUE ! JUMP BACK TO HERE AFTER REOPTIMISATION OF BAD ENDPOINTS498: 555   CONTINUE ! JUMP BACK TO HERE AFTER REOPTIMISATION OF BAD ENDPOINTS
501:       IF (RPHT) THEN499:       IF (RPHT) THEN
502:          CALL RPH500:          CALL RPH
503:       ELSE IF (MCPATHT) THEN501:       ELSE IF (MCPATHT) THEN
504:          CALL MCPATH502:          CALL MCPATH
505:       ELSE IF (MCPATH2T) THEN503:       ELSE IF (MCPATH2T) THEN
506:          CALL MCPATH2504:          CALL MCPATH2
507:       ELSE IF (CALCRATES.AND.READPATH) THEN505:       ELSE IF (CALCRATES.AND.READPATH) THEN
508:          CALL RATES(NATOMS,NINTS) ! JMC IF UNRES, THEN HAVE NON-ZERO NINTS, OTHERWISE JUST PASS 0506:          CALL RATES(NATOMS,NINTS) ! JMC IF UNRES, THEN HAVE NON-ZERO NINTS, OTHERWISE JUST PASS 0
 507:       
 508: 
 509: 
509:       ELSE IF (GUESSPATHT.AND.(.NOT.CONNECTT)) THEN510:       ELSE IF (GUESSPATHT.AND.(.NOT.CONNECTT)) THEN
510:          CALL NEWMINDIST(FIN,Q,NATOMS,DISTSF,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)511:          CALL NEWMINDIST(FIN,Q,NATOMS,DISTSF,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)
511:          WRITE(*,'(A,F12.2)') ' OPTIM> distance between start and finish=',DISTSF512:          WRITE(*,'(A,F12.2)') ' OPTIM> distance between start and finish=',DISTSF
512:          IF (UNRST) THEN513:          IF (UNRST) THEN
513:             DO J1=1,NRES514:             DO J1=1,NRES
514:                C(1,J1)=Q(6*(J1-1)+1)515:                C(1,J1)=Q(6*(J1-1)+1)
515:                C(2,J1)=Q(6*(J1-1)+2)516:                C(2,J1)=Q(6*(J1-1)+2)
516:                C(3,J1)=Q(6*(J1-1)+3)517:                C(3,J1)=Q(6*(J1-1)+3)
517:                C(1,J1+NRES)=Q(6*(J1-1)+4)518:                C(1,J1+NRES)=Q(6*(J1-1)+4)
518:                C(2,J1+NRES)=Q(6*(J1-1)+5)519:                C(2,J1+NRES)=Q(6*(J1-1)+5)


r33624/path.f 2018-01-19 09:34:36.338778115 +0000 r33623/path.f 2018-01-19 09:34:41.650823322 +0000
713: 713: 
714:       DO J1=1,NOPT714:       DO J1=1,NOPT
715:          QMINUS(J1)=Q(J1)715:          QMINUS(J1)=Q(J1)
716:       ENDDO716:       ENDDO
717:       EMINUS=ENERGY717:       EMINUS=ENERGY
718: C718: C
719: C  The total number of energies and coordinates is NSTEPPLUS + NSTEPMINUS + 1 for the transition state.719: C  The total number of energies and coordinates is NSTEPPLUS + NSTEPMINUS + 1 for the transition state.
720: C  The rest of this subroutine is post-processing 720: C  The rest of this subroutine is post-processing 
721: C  Minimise IO if PRINTPTS is .FALSE.721: C  Minimise IO if PRINTPTS is .FALSE.
722: C722: C
723:       IF (ZSYMSAVE(1:1).EQ.'W') THEN  !  WCOMMENT723:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN  !  WCOMMENT
724:          ALLOCATE(QFRAMEP(9*(NATOMS/2),NFMAX),QFRAMEM(9*(NATOMS/2),NFMAX))724:          ALLOCATE(QFRAMEP(9*(NATOMS/2),NFMAX),QFRAMEM(9*(NATOMS/2),NFMAX))
725:       ELSE725:       ELSE
726:          ALLOCATE(QFRAMEP(NOPT,NFMAX),QFRAMEM(NOPT,NFMAX))726:          ALLOCATE(QFRAMEP(NOPT,NFMAX),QFRAMEM(NOPT,NFMAX))
727: !        ALLOCATE(QFRAMEP(3*NATOMS,NFMAX),QFRAMEM(3*NATOMS,NFMAX))727: !        ALLOCATE(QFRAMEP(3*NATOMS,NFMAX),QFRAMEM(3*NATOMS,NFMAX))
728:       ENDIF728:       ENDIF
729:       IF (.NOT.PRINTPTS) THEN 729:       IF (.NOT.PRINTPTS) THEN 
730:          ALLOCATE(EOFS(3), PATHLENGTH(3), EOFSFRAMEP(3), EOFSFRAMEM(3))730:          ALLOCATE(EOFS(3), PATHLENGTH(3), EOFSFRAMEP(3), EOFSFRAMEM(3))
731:          NSTEPPLUS=1731:          NSTEPPLUS=1
732:          NSTEPMINUS=1732:          NSTEPMINUS=1
733:          NFPLUS=1733:          NFPLUS=1
762:       REWIND(2)762:       REWIND(2)
763:       J1=0763:       J1=0
764:       DO 764:       DO 
765:          READ(2,*,END=975) DUMMY765:          READ(2,*,END=975) DUMMY
766:          J1=J1+1766:          J1=J1+1
767:       ENDDO767:       ENDDO
768: 975   IF (DEBUG) PRINT '(A,I6)',' path> number of entries in EofS file=',J1768: 975   IF (DEBUG) PRINT '(A,I6)',' path> number of entries in EofS file=',J1
769:       REWIND(2)769:       REWIND(2)
770:       ALLOCATE(EOFS(J1), PATHLENGTH(J1), EOFSFRAMEP(J1), EOFSFRAMEM(J1))770:       ALLOCATE(EOFS(J1), PATHLENGTH(J1), EOFSFRAMEP(J1), EOFSFRAMEM(J1))
771:       READ(2,*) EOFS(NSTEPPLUS+1)771:       READ(2,*) EOFS(NSTEPPLUS+1)
772:       IF (ZSYMSAVE(1:1).EQ.'W') THEN  !  WCOMMENT772:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN  !  WCOMMENT
773:          ALLOCATE(Q1(9*(NATOMS/2)),Q2(9*(NATOMS/2)),QW(9*(NATOMS/2)))773:          ALLOCATE(Q1(9*(NATOMS/2)),Q2(9*(NATOMS/2)),QW(9*(NATOMS/2)))
774: C        READ(1,*) (Q1(J1),J1=1,9*NATOMS)774: C        READ(1,*) (Q1(J1),J1=1,9*NATOMS)
775:          READ(1,*) (Q1(J1),J1=1,9*(NATOMS/2))775:          READ(1,*) (Q1(J1),J1=1,9*(NATOMS/2))
776:       ELSE776:       ELSE
777:          ALLOCATE(Q1(NOPT),Q2(NOPT))777:          ALLOCATE(Q1(NOPT),Q2(NOPT))
778:          READ(1,*) (Q1(J1),J1=1,NOPT)778:          READ(1,*) (Q1(J1),J1=1,NOPT)
779:       ENDIF779:       ENDIF
780:       PATHLENGTH(NSTEPPLUS+1)=0.0D0780:       PATHLENGTH(NSTEPPLUS+1)=0.0D0
781: C     781: C     
782: C  The number of frames for each side of the path, NPATHFRAME, is now treated 782: C  The number of frames for each side of the path, NPATHFRAME, is now treated 
801:          PMINUS = 1.0d0801:          PMINUS = 1.0d0
802:       else802:       else
803:          PMINUS=MIN(MIN(NPATHFRAME,NFMAX)*1.0D0/(1.0D0*(NSTEPMINUS-1)),1.0D0)803:          PMINUS=MIN(MIN(NPATHFRAME,NFMAX)*1.0D0/(1.0D0*(NSTEPMINUS-1)),1.0D0)
804:       endif804:       endif
805: 805: 
806:       IF (PTEST) WRITE(*,'(A,F10.4,A,F10.4,A)') ' Frames will be dumped to points.path.xyz with probability ',806:       IF (PTEST) WRITE(*,'(A,F10.4,A,F10.4,A)') ' Frames will be dumped to points.path.xyz with probability ',
807:      1                          PPLUS,'/',PMINUS,' steps on the plus/minus sides'807:      1                          PPLUS,'/',PMINUS,' steps on the plus/minus sides'
808:       NFPLUS=0808:       NFPLUS=0
809:       DO J1=1,NSTEPPLUS809:       DO J1=1,NSTEPPLUS
810:          READ(2,*) EOFS(NSTEPPLUS+1-J1)810:          READ(2,*) EOFS(NSTEPPLUS+1-J1)
811:          IF (ZSYMSAVE(1:1).EQ.'W') THEN  811:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN  
812: C           READ(1,*) (Q2(J2),J2=1,9*NATOMS)  !  WCOMMENT812: C           READ(1,*) (Q2(J2),J2=1,9*NATOMS)  !  WCOMMENT
813:             READ(1,*) (Q2(J2),J2=1,9*(NATOMS/2))813:             READ(1,*) (Q2(J2),J2=1,9*(NATOMS/2))
814:          ELSE814:          ELSE
815:             READ(1,*) (Q2(J2),J2=1,NOPT)815:             READ(1,*) (Q2(J2),J2=1,NOPT)
816:          ENDIF816:          ENDIF
817: 817: 
818:          ! Only print the frame if ETEST is TRUE.818:          ! Only print the frame if ETEST is TRUE.
819:          ETEST=.FALSE.819:          ETEST=.FALSE.
820:          IF ((J1.EQ.1).OR.(J1.EQ.NSTEPPLUS)) THEN ! always take the end points820:          IF ((J1.EQ.1).OR.(J1.EQ.NSTEPPLUS)) THEN ! always take the end points
821:             ETEST=.TRUE.821:             ETEST=.TRUE.
826:          ! This is where we test whether to print the frame:826:          ! This is where we test whether to print the frame:
827:          IF (((DPRAND().LE.PPLUS).OR.(J1.EQ.NSTEPPLUS)).AND.(NFPLUS.LT.NFMAX).AND.ETEST) THEN827:          IF (((DPRAND().LE.PPLUS).OR.(J1.EQ.NSTEPPLUS)).AND.(NFPLUS.LT.NFMAX).AND.ETEST) THEN
828:             IF ((J1.EQ.NSTEPPLUS).OR.(NFPLUS.LT.NFMAX-1)) THEN ! save room for the last endpoint828:             IF ((J1.EQ.NSTEPPLUS).OR.(NFPLUS.LT.NFMAX-1)) THEN ! save room for the last endpoint
829:                NFPLUS=NFPLUS+1829:                NFPLUS=NFPLUS+1
830: 830: 
831:                ! Save the energy of the frame831:                ! Save the energy of the frame
832:                EOFSFRAMEP(NFPLUS)=EOFS(NSTEPPLUS+1-J1)832:                EOFSFRAMEP(NFPLUS)=EOFS(NSTEPPLUS+1-J1)
833:                LASTE=EOFS(NSTEPPLUS+1-J1)833:                LASTE=EOFS(NSTEPPLUS+1-J1)
834: 834: 
835:                ! Save the coordinates of the frame835:                ! Save the coordinates of the frame
836:                IF (ZSYMSAVE(1:1).EQ.'W') THEN836:                IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
837: C                 DO J2=1,9*NATOMS  !  WCOMMENT837: C                 DO J2=1,9*NATOMS  !  WCOMMENT
838:                   DO J2=1,9*(NATOMS/2)838:                   DO J2=1,9*(NATOMS/2)
839:                      QFRAMEP(J2,NFPLUS)=Q2(J2)839:                      QFRAMEP(J2,NFPLUS)=Q2(J2)
840:                   ENDDO840:                   ENDDO
841:                ELSE IF (ZSYMSAVE(1:2).EQ.'CD') THEN841:                ELSE IF (ZSYMSAVE(1:2).EQ.'CD') THEN
842:                    DO J2=1,NATOMS/2842:                    DO J2=1,NATOMS/2
843:                       CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),843:                       CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),
844:      1                              Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),844:      1                              Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),
845:      2                              CAPSCOORDS2,RAD,HEIGHT)845:      2                              CAPSCOORDS2,RAD,HEIGHT)
846:                       DO J3=1,18846:                       DO J3=1,18
860:          IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN860:          IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN
861:             CALL RB_DISTANCE(TEMP,Q1(1:DEGFREEDOMS),Q2(1:DEGFREEDOMS),DUMMYA,DUMMYB,.FALSE.)861:             CALL RB_DISTANCE(TEMP,Q1(1:DEGFREEDOMS),Q2(1:DEGFREEDOMS),DUMMYA,DUMMYB,.FALSE.)
862:             TEMP = TEMP**2862:             TEMP = TEMP**2
863:          ELSE IF (BULKT.AND..NOT.VASP) THEN863:          ELSE IF (BULKT.AND..NOT.VASP) THEN
864:             DO J2=1,NATOMS864:             DO J2=1,NATOMS
865:                TEMP=TEMP+MINIM(Q2(3*(J2-1)+1),Q1(3*(J2-1)+1),PARAM1)**2865:                TEMP=TEMP+MINIM(Q2(3*(J2-1)+1),Q1(3*(J2-1)+1),PARAM1)**2
866:      1                  +MINIM(Q2(3*(J2-1)+2),Q1(3*(J2-1)+2),PARAM2)**2866:      1                  +MINIM(Q2(3*(J2-1)+2),Q1(3*(J2-1)+2),PARAM2)**2
867:                IF (.NOT.TWOD) TEMP=TEMP+MINIM(Q2(3*(J2-1)+3),Q1(3*(J2-1)+3),PARAM3)**2867:                IF (.NOT.TWOD) TEMP=TEMP+MINIM(Q2(3*(J2-1)+3),Q1(3*(J2-1)+3),PARAM3)**2
868:             ENDDO868:             ENDDO
869:          ELSE869:          ELSE
870:             IF (ZSYMSAVE(1:1).EQ.'W') THEN870:             IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
871: C              DO J2=1,9*NATOMS  !  WCOMMENT871: C              DO J2=1,9*NATOMS  !  WCOMMENT
872:                DO J2=1,9*(NATOMS/2)872:                DO J2=1,9*(NATOMS/2)
873:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2873:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2
874:                ENDDO874:                ENDDO
875:             ELSE IF (ZSYMSAVE.EQ.'CD') THEN875:             ELSE IF (ZSYMSAVE.EQ.'CD') THEN
876:                DO J2=1,NATOMS/2876:                DO J2=1,NATOMS/2
877:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),877:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),
878:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),878:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),
879:      2                          CAPSCOORDS2,RAD,HEIGHT)879:      2                          CAPSCOORDS2,RAD,HEIGHT)
880:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),880:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),
888: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only888: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only
889:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)889:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)
890:             ELSE890:             ELSE
891:                DO J2=1,NOPT891:                DO J2=1,NOPT
892:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2892:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2
893:                ENDDO893:                ENDDO
894:             ENDIF894:             ENDIF
895:          ENDIF895:          ENDIF
896:          ! Compute the total pathlength from the TS to this frame896:          ! Compute the total pathlength from the TS to this frame
897:          PATHLENGTH(NSTEPPLUS+1-J1)=PATHLENGTH(NSTEPPLUS+2-J1)-SQRT(TEMP)897:          PATHLENGTH(NSTEPPLUS+1-J1)=PATHLENGTH(NSTEPPLUS+2-J1)-SQRT(TEMP)
898:          IF (ZSYMSAVE(1:1).EQ.'W') THEN898:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
899: C           DO J2=1,9*NATOMS  !  WCOMMENT899: C           DO J2=1,9*NATOMS  !  WCOMMENT
900:             DO J2=1,9*(NATOMS/2)900:             DO J2=1,9*(NATOMS/2)
901:                Q1(J2)=Q2(J2)901:                Q1(J2)=Q2(J2)
902:             ENDDO902:             ENDDO
903:          ELSE903:          ELSE
904:             ! Copy the current frame to the "previous frame" array904:             ! Copy the current frame to the "previous frame" array
905:             DO J2=1,NOPT905:             DO J2=1,NOPT
906:                Q1(J2)=Q2(J2)906:                Q1(J2)=Q2(J2)
907:             ENDDO907:             ENDDO
908:          ENDIF908:          ENDIF
909:       ENDDO  ! Loop over + direction frames909:       ENDDO  ! Loop over + direction frames
910:       IF (PTEST) WRITE(*,'(A,I6,A,I6,A)') ' Transition state will be frame number ',NFPLUS+1910:       IF (PTEST) WRITE(*,'(A,I6,A,I6,A)') ' Transition state will be frame number ',NFPLUS+1
911: 911: 
912:       ! Now set up for the - side of the path912:       ! Now set up for the - side of the path
913:       IF (ZSYMSAVE(1:1).EQ.'W') THEN913:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
914: C        DO J1=1,NATOMS  !  WCOMMENT914: C        DO J1=1,NATOMS  !  WCOMMENT
915:          DO J1=1,NATOMS/2915:          DO J1=1,NATOMS/2
916:             CALL CONVERT(QINIT(3*(J1-1)+1),QINIT(3*(J1-1)+2),QINIT(3*(J1-1)+3),916:             CALL CONVERT(QINIT(3*(J1-1)+1),QINIT(3*(J1-1)+2),QINIT(3*(J1-1)+3),
917: C    1                   QINIT(3*(NATOMS+J1-1)+1),QINIT(3*(NATOMS+J1-1)+2),QINIT(3*(NATOMS+J1-1)+3),917: C    1                   QINIT(3*(NATOMS+J1-1)+1),QINIT(3*(NATOMS+J1-1)+2),QINIT(3*(NATOMS+J1-1)+3),
918:      1                   QINIT(3*(NATOMS/2+J1-1)+1),QINIT(3*(NATOMS/2+J1-1)+2),QINIT(3*(NATOMS/2+J1-1)+3),918:      1                   QINIT(3*(NATOMS/2+J1-1)+1),QINIT(3*(NATOMS/2+J1-1)+2),QINIT(3*(NATOMS/2+J1-1)+3),
919:      2                   OVEC,H1VEC,H2VEC)919:      2                   OVEC,H1VEC,H2VEC)
920:             Q1(9*(J1-1)+1)=OVEC(1)920:             Q1(9*(J1-1)+1)=OVEC(1)
921:             Q1(9*(J1-1)+2)=OVEC(2)921:             Q1(9*(J1-1)+2)=OVEC(2)
922:             Q1(9*(J1-1)+3)=OVEC(3)922:             Q1(9*(J1-1)+3)=OVEC(3)
923:             Q1(9*(J1-1)+4)=H1VEC(1)923:             Q1(9*(J1-1)+4)=H1VEC(1)
929:          ENDDO929:          ENDDO
930:       ELSE930:       ELSE
931:          DO J1=1,NOPT931:          DO J1=1,NOPT
932:             Q1(J1)=QINIT(J1)932:             Q1(J1)=QINIT(J1)
933:          ENDDO933:          ENDDO
934:       ENDIF934:       ENDIF
935: 935: 
936:       NFMINUS=0936:       NFMINUS=0
937:       DO J1=1,NSTEPMINUS937:       DO J1=1,NSTEPMINUS
938:          READ(2,*) EOFS(NSTEPPLUS+1+J1)938:          READ(2,*) EOFS(NSTEPPLUS+1+J1)
939:          IF (ZSYMSAVE(1:1).EQ.'W') THEN939:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
940: C           READ(1,*) (Q2(J2),J2=1,9*NATOMS)  !  WCOMMENT940: C           READ(1,*) (Q2(J2),J2=1,9*NATOMS)  !  WCOMMENT
941:             READ(1,*) (Q2(J2),J2=1,9*(NATOMS/2))941:             READ(1,*) (Q2(J2),J2=1,9*(NATOMS/2))
942:          ELSE942:          ELSE
943:             READ(1,*) (Q2(J2),J2=1,NOPT)943:             READ(1,*) (Q2(J2),J2=1,NOPT)
944:          ENDIF944:          ENDIF
945: 945: 
946:          ETEST=.FALSE.946:          ETEST=.FALSE.
947:          IF ((J1.EQ.1).OR.(J1.EQ.NSTEPMINUS)) THEN ! always take the end points947:          IF ((J1.EQ.1).OR.(J1.EQ.NSTEPMINUS)) THEN ! always take the end points
948:             ETEST=.TRUE.948:             ETEST=.TRUE.
949:          ELSE949:          ELSE
950:             IF (ABS(EOFS(NSTEPPLUS+1+J1)-LASTE).GE.FRAMEEDIFF) ETEST=.TRUE. ! energy must have changed enough950:             IF (ABS(EOFS(NSTEPPLUS+1+J1)-LASTE).GE.FRAMEEDIFF) ETEST=.TRUE. ! energy must have changed enough
951:          ENDIF951:          ENDIF
952:          IF (ETEST) LASTE=EOFS(NSTEPPLUS+1+J1)952:          IF (ETEST) LASTE=EOFS(NSTEPPLUS+1+J1)
953:          RANDOM=DPRAND()953:          RANDOM=DPRAND()
954:          IF (((RANDOM.LE.PMINUS).OR.(J1.EQ.NSTEPMINUS)).AND.(NFMINUS.LE.NFMAX).AND.ETEST) THEN954:          IF (((RANDOM.LE.PMINUS).OR.(J1.EQ.NSTEPMINUS)).AND.(NFMINUS.LE.NFMAX).AND.ETEST) THEN
955:             IF ((J1.EQ.NSTEPMINUS).OR.(NFMINUS.LT.NFMAX-1)) THEN ! save a space for the stationary point at the end955:             IF ((J1.EQ.NSTEPMINUS).OR.(NFMINUS.LT.NFMAX-1)) THEN ! save a space for the stationary point at the end
956:                NFMINUS=NFMINUS+1956:                NFMINUS=NFMINUS+1
957:                EOFSFRAMEM(NFMINUS)=EOFS(NSTEPPLUS+1+J1)957:                EOFSFRAMEM(NFMINUS)=EOFS(NSTEPPLUS+1+J1)
958:                LASTE=EOFS(NSTEPPLUS+1+J1)958:                LASTE=EOFS(NSTEPPLUS+1+J1)
959:                IF (ZSYMSAVE(1:1).EQ.'W') THEN959:                IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
960: C                 DO J2=1,9*NATOMS  !  WCOMMENT960: C                 DO J2=1,9*NATOMS  !  WCOMMENT
961:                      DO J2=1,9*(NATOMS/2)961:                      DO J2=1,9*(NATOMS/2)
962:                      QFRAMEM(J2,NFMINUS)=Q2(J2)962:                      QFRAMEM(J2,NFMINUS)=Q2(J2)
963:                   ENDDO963:                   ENDDO
964:                ELSE IF (ZSYMSAVE.EQ.'CD') THEN964:                ELSE IF (ZSYMSAVE.EQ.'CD') THEN
965:                    DO J2=1,NATOMS/2965:                    DO J2=1,NATOMS/2
966:                       CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),966:                       CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),
967:      1                              Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),967:      1                              Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),
968:      2                              CAPSCOORDS2,RAD,HEIGHT)968:      2                              CAPSCOORDS2,RAD,HEIGHT)
969:                       DO J3=1,18969:                       DO J3=1,18
983:          IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN983:          IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN
984:             CALL RB_DISTANCE(TEMP,Q1(1:DEGFREEDOMS),Q2(1:DEGFREEDOMS),DUMMYA,DUMMYB,.FALSE.)984:             CALL RB_DISTANCE(TEMP,Q1(1:DEGFREEDOMS),Q2(1:DEGFREEDOMS),DUMMYA,DUMMYB,.FALSE.)
985:             TEMP = TEMP**2985:             TEMP = TEMP**2
986:          ELSE IF (BULKT.AND..NOT.VASP) THEN986:          ELSE IF (BULKT.AND..NOT.VASP) THEN
987:             DO J2=1,NATOMS987:             DO J2=1,NATOMS
988:                TEMP=TEMP+MINIM(Q2(3*(J2-1)+1),Q1(3*(J2-1)+1),PARAM1)**2988:                TEMP=TEMP+MINIM(Q2(3*(J2-1)+1),Q1(3*(J2-1)+1),PARAM1)**2
989:      1                  +MINIM(Q2(3*(J2-1)+2),Q1(3*(J2-1)+2),PARAM2)**2989:      1                  +MINIM(Q2(3*(J2-1)+2),Q1(3*(J2-1)+2),PARAM2)**2
990:                IF (.NOT.TWOD) TEMP=TEMP+MINIM(Q2(3*(J2-1)+3),Q1(3*(J2-1)+3),PARAM3)**2990:                IF (.NOT.TWOD) TEMP=TEMP+MINIM(Q2(3*(J2-1)+3),Q1(3*(J2-1)+3),PARAM3)**2
991:             ENDDO991:             ENDDO
992:          ELSE992:          ELSE
993:             IF (ZSYMSAVE(1:1).EQ.'W') THEN993:             IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
994: C              DO J2=1,9*NATOMS  !  WCOMMENT994: C              DO J2=1,9*NATOMS  !  WCOMMENT
995:                DO J2=1,9*(NATOMS/2)995:                DO J2=1,9*(NATOMS/2)
996:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2996:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2
997:                ENDDO997:                ENDDO
998:             ELSE IF (ZSYMSAVE.EQ.'CD') THEN998:             ELSE IF (ZSYMSAVE.EQ.'CD') THEN
999:                DO J2=1,NATOMS/2999:                DO J2=1,NATOMS/2
1000:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),1000:                   CALL CAPSIDIO(Q2(3*(J2-1)+1),Q2(3*(J2-1)+2),Q2(3*(J2-1)+3),
1001:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),1001:      1                          Q2(3*(NATOMS/2+J2-1)+1),Q2(3*(NATOMS/2+J2-1)+2),Q2(3*(NATOMS/2+J2-1)+3),
1002:      2                          CAPSCOORDS2,RAD,HEIGHT)1002:      2                          CAPSCOORDS2,RAD,HEIGHT)
1003:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),1003:                   CALL CAPSIDIO(Q1(3*(J2-1)+1),Q1(3*(J2-1)+2),Q1(3*(J2-1)+3),
1010:             ELSE IF((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT) THEN1010:             ELSE IF((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT) THEN
1011: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only1011: C           calculate path lengths wrt. Cartesian coordinate of the centre and two coordinates along the symmetry axis only
1012:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)1012:                CALL UNIAXGETPATHLENGTH(Q1,Q2,TEMP)
1013:             ELSE1013:             ELSE
1014:                DO J2=1,NOPT1014:                DO J2=1,NOPT
1015:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**21015:                   TEMP=TEMP+(Q2(J2)-Q1(J2))**2
1016:                ENDDO1016:                ENDDO
1017:             ENDIF1017:             ENDIF
1018:          ENDIF1018:          ENDIF
1019:          PATHLENGTH(NSTEPPLUS+1+J1)=PATHLENGTH(NSTEPPLUS+J1)+SQRT(TEMP)1019:          PATHLENGTH(NSTEPPLUS+1+J1)=PATHLENGTH(NSTEPPLUS+J1)+SQRT(TEMP)
1020:          IF (ZSYMSAVE(1:1).EQ.'W') THEN1020:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1021: C           DO J2=1,9*NATOMS  !  WCOMMENT1021: C           DO J2=1,9*NATOMS  !  WCOMMENT
1022:             DO J2=1,9*(NATOMS/2)1022:             DO J2=1,9*(NATOMS/2)
1023:                Q1(J2)=Q2(J2)1023:                Q1(J2)=Q2(J2)
1024:             ENDDO1024:             ENDDO
1025:          ELSE1025:          ELSE
1026:             DO J2=1,NOPT1026:             DO J2=1,NOPT
1027:                Q1(J2)=Q2(J2)1027:                Q1(J2)=Q2(J2)
1028:             ENDDO1028:             ENDDO
1029:          ENDIF1029:          ENDIF
1030:       ENDDO1030:       ENDDO
1033:       ! Write the data we have just obtained to the EofS file for this path1033:       ! Write the data we have just obtained to the EofS file for this path
1034:       OPEN(UNIT=3,FILE=EOFSSTRING,STATUS='UNKNOWN')1034:       OPEN(UNIT=3,FILE=EOFSSTRING,STATUS='UNKNOWN')
1035:       WRITE(3,'(2G25.15,I6)') (PATHLENGTH(J1),EOFS(J1),J1,J1=1,NSTEPPLUS+NSTEPMINUS+1)1035:       WRITE(3,'(2G25.15,I6)') (PATHLENGTH(J1),EOFS(J1),J1,J1=1,NSTEPPLUS+NSTEPMINUS+1)
1036:       CLOSE(3)1036:       CLOSE(3)
1037:       1037:       
1038:       SUM2=0.0D01038:       SUM2=0.0D0
1039:       SUM4=0.0D01039:       SUM4=0.0D0
1040: C     NDUMMY=1 !  WCOMMENT1040: C     NDUMMY=1 !  WCOMMENT
1041: C     IF (ZSYMSAVE(1:1).EQ.'W') NDUMMY=9  !  WCOMMENT 9 should have been 3?1041: C     IF (ZSYMSAVE(1:1).EQ.'W') NDUMMY=9  !  WCOMMENT 9 should have been 3?
1042:       NATOMSIMUL=NATOMS1042:       NATOMSIMUL=NATOMS
1043:       IF (ZSYMSAVE(1:1).EQ.'W') NATOMSIMUL=3*(NATOMS/2)1043:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') NATOMSIMUL=3*(NATOMS/2)
1044:       IF (ZSYMSAVE(1:1).EQ.'W') THEN1044:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1045:          DO J1=1,NATOMSIMUL  !  WCOMMENT1045:          DO J1=1,NATOMSIMUL  !  WCOMMENT
1046:             SUM2=SUM2+1046:             SUM2=SUM2+
1047:      1           (QFRAMEP(3*(J1-1)+1,NFPLUS)-QFRAMEM(3*(J1-1)+1,NFMINUS))**21047:      1           (QFRAMEP(3*(J1-1)+1,NFPLUS)-QFRAMEM(3*(J1-1)+1,NFMINUS))**2
1048:      2          +(QFRAMEP(3*(J1-1)+2,NFPLUS)-QFRAMEM(3*(J1-1)+2,NFMINUS))**21048:      2          +(QFRAMEP(3*(J1-1)+2,NFPLUS)-QFRAMEM(3*(J1-1)+2,NFMINUS))**2
1049:      3          +(QFRAMEP(3*(J1-1)+3,NFPLUS)-QFRAMEM(3*(J1-1)+3,NFMINUS))**21049:      3          +(QFRAMEP(3*(J1-1)+3,NFPLUS)-QFRAMEM(3*(J1-1)+3,NFMINUS))**2
1050:             SUM4=SUM4+1050:             SUM4=SUM4+
1051:      1          ((QFRAMEP(3*(J1-1)+1,NFPLUS)-QFRAMEM(3*(J1-1)+1,NFMINUS))**21051:      1          ((QFRAMEP(3*(J1-1)+1,NFPLUS)-QFRAMEM(3*(J1-1)+1,NFMINUS))**2
1052:      2          +(QFRAMEP(3*(J1-1)+2,NFPLUS)-QFRAMEM(3*(J1-1)+2,NFMINUS))**21052:      2          +(QFRAMEP(3*(J1-1)+2,NFPLUS)-QFRAMEM(3*(J1-1)+2,NFMINUS))**2
1053:      3          +(QFRAMEP(3*(J1-1)+3,NFPLUS)-QFRAMEM(3*(J1-1)+3,NFMINUS))**2)**21053:      3          +(QFRAMEP(3*(J1-1)+3,NFPLUS)-QFRAMEM(3*(J1-1)+3,NFMINUS))**2)**2
1054:          ENDDO1054:          ENDDO
1355: C           WRITE(3,'(A2,4X,3g20.10)') 'C4  ',CAPSCOORDS2(16),CAPSCOORDS2(17),CAPSCOORDS2(18)1355: C           WRITE(3,'(A2,4X,3g20.10)') 'C4  ',CAPSCOORDS2(16),CAPSCOORDS2(17),CAPSCOORDS2(18)
1356: C        ENDDO1356: C        ENDDO
1357: C     ELSE1357: C     ELSE
1358: C        WRITE(3,'(I6)') NATOMS1358: C        WRITE(3,'(I6)') NATOMS
1359: C        WRITE(3,'(g20.10)') EOFS(1)1359: C        WRITE(3,'(g20.10)') EOFS(1)
1360: C        WRITE(3,'(A2,4X,3g20.10)') (ZSYM(J1),QPLUS(3*(J1-1)+1),QPLUS(3*(J1-1)+2),QPLUS(3*(J1-1)+3),J1=1,NATOMS)1360: C        WRITE(3,'(A2,4X,3g20.10)') (ZSYM(J1),QPLUS(3*(J1-1)+1),QPLUS(3*(J1-1)+2),QPLUS(3*(J1-1)+3),J1=1,NATOMS)
1361: C     ENDIF1361: C     ENDIF
1362: 1362: 
1363:       ! Write + frame coordinates for this path1363:       ! Write + frame coordinates for this path
1364:       DO J1=NFPLUS,1,-11364:       DO J1=NFPLUS,1,-1
1365:          IF (ZSYMSAVE(1:1).EQ.'W') THEN1365:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1366: C           WRITE(3,'(I6)') 3*NATOMS  !  WCOMMENT1366: C           WRITE(3,'(I6)') 3*NATOMS  !  WCOMMENT
1367:             WRITE(3,'(I6)') 3*(NATOMS/2)1367:             WRITE(3,'(I6)') 3*(NATOMS/2)
1368:             WRITE(3,'(A,g25.15)') 'Energy=',EOFSFRAMEP(J1)1368:             WRITE(3,'(A,g25.15)') 'Energy=',EOFSFRAMEP(J1)
1369: C           IF (J1.EQ.NFPLUS) THEN 1369: C           IF (J1.EQ.NFPLUS) THEN 
1370: C              WRITE(3,'(g20.10)') EOFS(1)1370: C              WRITE(3,'(g20.10)') EOFS(1)
1371: C           ELSE1371: C           ELSE
1372: C              WRITE(3,'(A)') ' '1372: C              WRITE(3,'(A)') ' '
1373: C           ENDIF1373: C           ENDIF
1374:             WRITE(3,'(A2,4X,3g20.10)') 1374:             WRITE(3,'(A2,4X,3g20.10)') 
1375:      1   ('O  ',QFRAMEP(9*(J2-1)+1,J1),QFRAMEP(9*(J2-1)+2,J1),QFRAMEP(9*(J2-1)+3,J1),1375:      1   ('O  ',QFRAMEP(9*(J2-1)+1,J1),QFRAMEP(9*(J2-1)+2,J1),QFRAMEP(9*(J2-1)+3,J1),
1503:      &         (ZSYM(J2),QFRAMEP(3*(J2-1)+1,J1),QFRAMEP(3*(J2-1)+2,J1),QFRAMEP(3*(J2-1)+3,J1),J2=1,NATOMS)1503:      &         (ZSYM(J2),QFRAMEP(3*(J2-1)+1,J1),QFRAMEP(3*(J2-1)+2,J1),QFRAMEP(3*(J2-1)+3,J1),J2=1,NATOMS)
1504:             ENDIF1504:             ENDIF
1505:          ENDIF1505:          ENDIF
1506:       ENDDO1506:       ENDDO
1507: 1507: 
1508: 1508: 
1509:       ! Now write coords of the transition state to file.1509:       ! Now write coords of the transition state to file.
1510: !1510: !
1511: !  ZSYMSAVE=W,CD  1511: !  ZSYMSAVE=W,CD  
1512: !  STOCKT RBAAT AMHT RINGPOLYMERT ELSE1512: !  STOCKT RBAAT AMHT RINGPOLYMERT ELSE
1513:       IF (ZSYMSAVE(1:1).EQ.'W') THEN1513:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1514: C        WRITE(3,'(I6)') 3*NATOMS ! WCOMMENT1514: C        WRITE(3,'(I6)') 3*NATOMS ! WCOMMENT
1515:          WRITE(3,'(I6)') 3*(NATOMS/2)1515:          WRITE(3,'(I6)') 3*(NATOMS/2)
1516:          WRITE(3,'(A,g25.15)') 'Energy=',EOFS(NSTEPPLUS+1)1516:          WRITE(3,'(A,g25.15)') 'Energy=',EOFS(NSTEPPLUS+1)
1517: C        DO J1=1,NATOMS ! WCOMMENT1517: C        DO J1=1,NATOMS ! WCOMMENT
1518:          DO J1=1,NATOMS/21518:          DO J1=1,NATOMS/2
1519:             CALL CONVERT(QINIT(3*(J1-1)+1),QINIT(3*(J1-1)+2),QINIT(3*(J1-1)+3),1519:             CALL CONVERT(QINIT(3*(J1-1)+1),QINIT(3*(J1-1)+2),QINIT(3*(J1-1)+3),
1520: C    1                   QINIT(3*(NATOMS+J1-1)+1),QINIT(3*(NATOMS+J1-1)+2),QINIT(3*(NATOMS+J1-1)+3),1520: C    1                   QINIT(3*(NATOMS+J1-1)+1),QINIT(3*(NATOMS+J1-1)+2),QINIT(3*(NATOMS+J1-1)+3),
1521:      1                   QINIT(3*(NATOMS/2+J1-1)+1),QINIT(3*(NATOMS/2+J1-1)+2),QINIT(3*(NATOMS/2+J1-1)+3),1521:      1                   QINIT(3*(NATOMS/2+J1-1)+1),QINIT(3*(NATOMS/2+J1-1)+2),QINIT(3*(NATOMS/2+J1-1)+3),
1522:      2                   OVEC,H1VEC,H2VEC)1522:      2                   OVEC,H1VEC,H2VEC)
1523:             WRITE(3,'(A2,4X,3g20.10)') 'O  ',OVEC(1),OVEC(2),OVEC(3)1523:             WRITE(3,'(A2,4X,3g20.10)') 'O  ',OVEC(1),OVEC(2),OVEC(3)
1621:             enddo1621:             enddo
1622:          ELSEIF (VARIABLES) THEN1622:          ELSEIF (VARIABLES) THEN
1623:             WRITE(3,'(G20.10)') QINIT(1:NOPT)1623:             WRITE(3,'(G20.10)') QINIT(1:NOPT)
1624:          ELSE1624:          ELSE
1625:             WRITE(3,'(A2,4X,3G20.10)') (ZSYM(J1),QINIT(3*(J1-1)+1),QINIT(3*(J1-1)+2),QINIT(3*(J1-1)+3),J1=1,NATOMS)1625:             WRITE(3,'(A2,4X,3G20.10)') (ZSYM(J1),QINIT(3*(J1-1)+1),QINIT(3*(J1-1)+2),QINIT(3*(J1-1)+3),J1=1,NATOMS)
1626:          endif1626:          endif
1627:       ENDIF1627:       ENDIF
1628: 1628: 
1629:       ! Now - path frames1629:       ! Now - path frames
1630:       DO J1=1,NFMINUS1630:       DO J1=1,NFMINUS
1631:          IF (ZSYMSAVE(1:1).EQ.'W') THEN1631:          IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1632: C           WRITE(3,'(I6)') 3*NATOMS  !  WCOMMENT1632: C           WRITE(3,'(I6)') 3*NATOMS  !  WCOMMENT
1633:             WRITE(3,'(I6)') 3*(NATOMS/2)1633:             WRITE(3,'(I6)') 3*(NATOMS/2)
1634:             WRITE(3,'(A,g25.15)') 'Energy=',EOFSFRAMEM(J1)1634:             WRITE(3,'(A,g25.15)') 'Energy=',EOFSFRAMEM(J1)
1635: C           IF (J1.EQ.NFMINUS) THEN1635: C           IF (J1.EQ.NFMINUS) THEN
1636: C              WRITE(3,'(g20.10)') EOFS(NSTEPPLUS+NSTEPMINUS+1)1636: C              WRITE(3,'(g20.10)') EOFS(NSTEPPLUS+NSTEPMINUS+1)
1637: C           ELSE1637: C           ELSE
1638: C              WRITE(3,'(A)') ' '1638: C              WRITE(3,'(A)') ' '
1639: C           ENDIF1639: C           ENDIF
1640:             WRITE(3,'(A2,4X,3g20.10)') 1640:             WRITE(3,'(A2,4X,3g20.10)') 
1641:      1            ('O  ',QFRAMEM(9*(J2-1)+1,J1),QFRAMEM(9*(J2-1)+2,J1),QFRAMEM(9*(J2-1)+3,J1),1641:      1            ('O  ',QFRAMEM(9*(J2-1)+1,J1),QFRAMEM(9*(J2-1)+2,J1),QFRAMEM(9*(J2-1)+3,J1),
1932: 1932: 
1933:       ! We have identified a stationary point (minimum or ts) and now wish to write its information to a path.info file.1933:       ! We have identified a stationary point (minimum or ts) and now wish to write its information to a path.info file.
1934:       ! We know the energy and coordinates of the point, we need to determine the point group and the normal mode frequencies.1934:       ! We know the energy and coordinates of the point, we need to determine the point group and the normal mode frequencies.
1935:       ! This is done here.1935:       ! This is done here.
1936:       ! (Note, sometimes we may have already calculated frequencies but if so they have probably been shifted, so we don't1936:       ! (Note, sometimes we may have already calculated frequencies but if so they have probably been shifted, so we don't
1937:       ! use them. We calculate new ones instead)1937:       ! use them. We calculate new ones instead)
1938: 1938: 
1939: 1939: 
1940:       ! First, we deal with a couple of the more special cases. Both of these should be largely obsolete now: the GENRIGID1940:       ! First, we deal with a couple of the more special cases. Both of these should be largely obsolete now: the GENRIGID
1941:       ! versions of rigid body potentials (including the TIP potentials) are more likely to be maintained.1941:       ! versions of rigid body potentials (including the TIP potentials) are more likely to be maintained.
1942:       IF (ZSYMSAVE(1:1).EQ.'W') THEN1942:       IF (ZSYMSAVE(1:1).EQ.'W'.and.ZSYMSAVE(1:2).NE.'WC') THEN
1943:          IF (ZSYMSAVE.EQ.'W4') IPOT=41943:          IF (ZSYMSAVE.EQ.'W4') IPOT=4
1944:          IF (ZSYMSAVE.EQ.'W3') IPOT=31944:          IF (ZSYMSAVE.EQ.'W3') IPOT=3
1945:          IF (ZSYMSAVE.EQ.'W2') IPOT=21945:          IF (ZSYMSAVE.EQ.'W2') IPOT=2
1946:          IF (ZSYMSAVE.EQ.'W1') IPOT=11946:          IF (ZSYMSAVE.EQ.'W1') IPOT=1
1947: C        DO J2=1,NATOMS1947: C        DO J2=1,NATOMS
1948:          DO J2=1,NATOMS/2 ! WCOMMENT1948:          DO J2=1,NATOMS/2 ! WCOMMENT
1949:             CALL CONVERT(Q(3*(J2-1)+1),Q(3*(J2-1)+2),Q(3*(J2-1)+3),1949:             CALL CONVERT(Q(3*(J2-1)+1),Q(3*(J2-1)+2),Q(3*(J2-1)+3),
1950: C    1                Q(3*(NATOMS+J2-1)+1),Q(3*(NATOMS+J2-1)+2),Q(3*(NATOMS+J2-1)+3),1950: C    1                Q(3*(NATOMS+J2-1)+1),Q(3*(NATOMS+J2-1)+2),Q(3*(NATOMS+J2-1)+3),
1951:      1                Q(3*(NATOMS/2+J2-1)+1),Q(3*(NATOMS/2+J2-1)+2),Q(3*(NATOMS/2+J2-1)+3),1951:      1                Q(3*(NATOMS/2+J2-1)+1),Q(3*(NATOMS/2+J2-1)+2),Q(3*(NATOMS/2+J2-1)+3),
1952:      2                OVEC,H1VEC,H2VEC)1952:      2                OVEC,H1VEC,H2VEC)


r33624/pertable.f 2018-01-19 09:34:36.586780166 +0000 r33623/pertable.f 2018-01-19 09:34:41.882825352 +0000
 90:       IF (IPRNT .GE. 10) WRITE (*,9000) 90:       IF (IPRNT .GE. 10) WRITE (*,9000)
 91:  9000 FORMAT ('PERTABLE: Mass and at. nr. lookup'/ 91:  9000 FORMAT ('PERTABLE: Mass and at. nr. lookup'/
 92:      1     'Line Symbol AtNr   At. Mass') 92:      1     'Line Symbol AtNr   At. Mass')
 93:       DO 10 I = 1, NATOMS 93:       DO 10 I = 1, NATOMS
 94:          IF (ZSYM(I) .EQ. 'X ') THEN 94:          IF (ZSYM(I) .EQ. 'X ') THEN
 95:             IATNUM(I) = 0 95:             IATNUM(I) = 0
 96:             ATMASS(I) = 0.0 96:             ATMASS(I) = 0.0
 97:          ELSE 97:          ELSE
 98:             DO 20 J = 1, MXTABL 98:             DO 20 J = 1, MXTABL
 99:                IF (ZSYM(I) .EQ. ATSYM(J)) THEN 99:                IF (ZSYM(I) .EQ. ATSYM(J)) THEN
100:                   IF (ZSYM(I)(1:1).EQ.'W') THEN ! to generalise this to other rigid 100:                   IF (ZSYM(I)(1:1).EQ.'W'.and.ZSYM(I)(1:2).NE.'WC') THEN ! to generalise this to other rigid 
101:                      IF (I.LE.NATOMS/2) THEN101:                      IF (I.LE.NATOMS/2) THEN
102:                         IATNUM(3*(I-1)+1)=8        ! bodies will require site/mass specification !102:                         IATNUM(3*(I-1)+1)=8        ! bodies will require site/mass specification !
103:                         IATNUM(3*(I-1)+2)=1103:                         IATNUM(3*(I-1)+2)=1
104:                         IATNUM(3*(I-1)+3)=1104:                         IATNUM(3*(I-1)+3)=1
105:                         ATMASS(3*(I-1)+1)=16.0D0105:                         ATMASS(3*(I-1)+1)=16.0D0
106:                         ATMASS(3*(I-1)+2)=1.0D0106:                         ATMASS(3*(I-1)+2)=1.0D0
107:                         ATMASS(3*(I-1)+3)=1.0D0107:                         ATMASS(3*(I-1)+3)=1.0D0
108:                      ENDIF108:                      ENDIF
109:                   ELSE109:                   ELSE
110:                      IATNUM(I) = J110:                      IATNUM(I) = J


r33624/potential.f 2018-01-19 09:34:36.866782479 +0000 r33623/potential.f 2018-01-19 09:34:42.238828487 +0000
  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/potential.f' in revision 33624  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/potential.f' in revision 33623


r33624/wca.f 2018-01-19 09:34:37.106784475 +0000 r33623/wca.f 2018-01-19 09:34:42.494830757 +0000
  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/wca.f' in revision 33624  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/wca.f' in revision 33623


r33624/wca_tail_modified.f 2018-01-19 09:34:37.330786349 +0000 r33623/wca_tail_modified.f 2018-01-19 09:34:42.782833304 +0000
  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/wca_tail_modified.f' in revision 33624  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/wca_tail_modified.f' in revision 33623


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0