hdiff output

r32316/bfgsts.f 2017-04-22 17:30:12.326869573 +0100 r32315/bfgsts.f 2017-04-22 17:30:16.678925807 +0100
219:             GOTO 20219:             GOTO 20
220: 10          CLOSE(45)220: 10          CLOSE(45)
221:          ENDIF221:          ENDIF
222:          WRITE(*,'(A,F20.10)') ' Reaction vector read successfully. Eigenvalue=   ',EVALMIN222:          WRITE(*,'(A,F20.10)') ' Reaction vector read successfully. Eigenvalue=   ',EVALMIN
223:          READV=.FALSE. ! needed for runs that do bfgsts and then path223:          READV=.FALSE. ! needed for runs that do bfgsts and then path
224:          NS=100224:          NS=100
225:       ENDIF225:       ENDIF
226:       IF (FREEZE) THEN226:       IF (FREEZE) THEN
227:          DO J1=1,NATOMS227:          DO J1=1,NATOMS
228:             IF (FROZEN(J1)) THEN228:             IF (FROZEN(J1)) THEN
229:                IF (VARIABLES) THEN229:                VECS(3*(J1-1)+1)=0.0D0
230:                   VECS(J1)=0.0D0230:                VECS(3*(J1-1)+2)=0.0D0
231:                ELSE231:                VECS(3*(J1-1)+3)=0.0D0
232:                   VECS(3*(J1-1)+1)=0.0D0 
233:                   VECS(3*(J1-1)+2)=0.0D0 
234:                   VECS(3*(J1-1)+3)=0.0D0 
235:                ENDIF 
236:             ENDIF232:             ENDIF
237:          ENDDO233:          ENDDO
238:       ENDIF234:       ENDIF
239: ! jwrm2> Make sure the z coordinate in 2D systems is fixed235: ! jwrm2> Make sure the z coordinate in 2D systems is fixed
240:       IF (TWOD .OR. GTHOMSONT) THEN236:       IF (TWOD .OR. GTHOMSONT) THEN
241:          DO J1 = 1,NATOMS237:          DO J1 = 1,NATOMS
242:             VECS(3*J1) = 0.0D0238:             VECS(3*J1) = 0.0D0
243:          ENDDO239:          ENDDO
244:       ENDIF240:       ENDIF
245: ! jwrm2> end241: ! jwrm2> end
297:                      HESS(3*(J3-1)+2,J2)=0.0D0293:                      HESS(3*(J3-1)+2,J2)=0.0D0
298:                      HESS(3*(J3-1)+3,J2)=0.0D0294:                      HESS(3*(J3-1)+3,J2)=0.0D0
299:                   ENDDO295:                   ENDDO
300:                   HESS(J2,J2)=SHIFTV296:                   HESS(J2,J2)=SHIFTV
301:                ENDDO297:                ENDDO
302:             ENDIF298:             ENDIF
303:          ENDIF299:          ENDIF
304:       ENDIF300:       ENDIF
305: !301: !
306: !302: !
307:       IF (VARIABLES.AND.FREEZE) THEN 
308:          DO J1=1,NATOMS 
309:             IF (.NOT.FROZEN(J1))  CYCLE 
310:             HESS(J1,J1)=HESS(J1,J1)+SHIFTL(1) 
311:          ENDDO 
312:       ENDIF 
313:  
314:       IF (VARIABLES.AND.FOURDPBCT) CALL SHIFTFOURDPBCT303:       IF (VARIABLES.AND.FOURDPBCT) CALL SHIFTFOURDPBCT
315:       IF (VARIABLES.AND.THREEDPBCT) CALL SHIFTTHREEDPBCT304:       IF (VARIABLES.AND.THREEDPBCT) CALL SHIFTTHREEDPBCT
316:       IF (VARIABLES.AND.TWODPBCT) CALL SHIFTTWODPBCT305:       IF (VARIABLES.AND.TWODPBCT) CALL SHIFTTWODPBCT
317:       IF (VARIABLES.AND.ONEDPBCT) CALL SHIFTONEDPBCT306:       IF (VARIABLES.AND.ONEDPBCT) CALL SHIFTONEDPBCT
318:       IF (VARIABLES.AND.INVTONEDPBCT) CALL SHIFTONEDPBCT307:       IF (VARIABLES.AND.INVTONEDPBCT) CALL SHIFTONEDPBCT
319:       IF (VARIABLES.AND.INVTTWODPBCT) CALL SHIFTTWODPBCT308:       IF (VARIABLES.AND.INVTTWODPBCT) CALL SHIFTTWODPBCT
320: !      IF (VARIABLES.AND.EX1DT) CALL SHIFTEX1D309: !      IF (VARIABLES.AND.EX1DT) CALL SHIFTEX1D
321: 310: 
322: !311: !
323: !  This giant IF block computes the lowest eigenvectors312: !  This giant IF block computes the lowest eigenvectors
1375: C        VEC1(1)=1.0D01364: C        VEC1(1)=1.0D0
1376: 1365: 
1377:       ELSE1366:       ELSE
1378:          DO J1=1,NOPT1367:          DO J1=1,NOPT
1379:             VEC1(J1)=VECL(J1)1368:             VEC1(J1)=VECL(J1)
1380:          ENDDO1369:          ENDDO
1381:       ENDIF1370:       ENDIF
1382:       IF (FREEZE) THEN1371:       IF (FREEZE) THEN
1383:          DO J1=1,NATOMS1372:          DO J1=1,NATOMS
1384:             IF (FROZEN(J1)) THEN1373:             IF (FROZEN(J1)) THEN
1385:                IF (VARIABLES) THEN1374:                VEC1(3*(J1-1)+1)=0.0D0
1386:                   VEC1(J1)=0.0D01375:                VEC1(3*(J1-1)+2)=0.0D0
1387:                ELSE1376:                VEC1(3*(J1-1)+3)=0.0D0
1388:                   VEC1(3*(J1-1)+1)=0.0D0 
1389:                   VEC1(3*(J1-1)+2)=0.0D0 
1390:                   VEC1(3*(J1-1)+3)=0.0D0 
1391:                ENDIF 
1392:             ENDIF1377:             ENDIF
1393:          ENDDO1378:          ENDDO
1394:          CALL VECNORM(VEC1,NOPT)1379:          CALL VECNORM(VEC1,NOPT)
1395:       ENDIF1380:       ENDIF
1396: 1381: 
1397: ! jwrm2> Make sure the z coordinate in 2D systems is fixed1382: ! jwrm2> Make sure the z coordinate in 2D systems is fixed
1398:       IF (TWOD .OR. GTHOMSONT) THEN1383:       IF (TWOD .OR. GTHOMSONT) THEN
1399:          DO J1 = 1,NATOMS1384:          DO J1 = 1,NATOMS
1400:             VECS(3*J1) = 0.0D01385:             VECS(3*J1) = 0.0D0
1401:          ENDDO1386:          ENDDO


r32316/efol.f90 2017-04-22 17:30:12.546872421 +0100 r32315/efol.f90 2017-04-22 17:30:16.902928687 +0100
 98: ! 98: !
 99: !  Reset maximum step sizes in case this isn't the first call to EFOL. 99: !  Reset maximum step sizes in case this isn't the first call to EFOL.
100: !100: !
101:       IF (DEBUG) PRINT '(A,G20.10)',' efol> resetting maximum step sizes to ',MXSTP101:       IF (DEBUG) PRINT '(A,G20.10)',' efol> resetting maximum step sizes to ',MXSTP
102:       DO J1=1,NOPT102:       DO J1=1,NOPT
103:          STPMAX(J1)=MXSTP103:          STPMAX(J1)=MXSTP
104:       ENDDO104:       ENDDO
105: !105: !
106: !106: !
107: !  This call to gmetry is needed for REDO runs with TIP potentials, otherwise if107: !  This call to gmetry is needed for REDO runs with TIP potentials, otherwise if
108: !  an Euler angle goes bad the step is messed up.108: !  an Euler angle goes bad the step is f*cked.
109: !109: !
110: !     CALL GMETRY(0,VEC,QTS) 
111:       VEC(1:NOPT)=0.0D0110:       VEC(1:NOPT)=0.0D0
 111:       CALL GMETRY(0,VEC,QTS)
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') 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) 
191:             CALL SHIFTRIGID(QTS,NATOMS)191:             CALL SHIFTRIGID(QTS,NATOMS)
192:          ELSEIF (ZSYM(NATOMS).EQ.'TH') THEN192:          ELSEIF (ZSYM(NATOMS).EQ.'TH') THEN
193:             CALL SHIFTHTH(QTS,NATOMS)193:             CALL SHIFTHTH(QTS,NATOMS)
194: ! hk286194: ! hk286
195:          ELSEIF (GTHOMSONT) THEN195:          ELSEIF (GTHOMSONT) THEN
196:             CALL SHIFTHGTH(QTS,NATOMS)196:             CALL SHIFTHGTH(QTS,NATOMS)
197:          ELSE197:          ELSE
198:             CALL SHIFTH(QTS,.TRUE.,NOPT,NATOMS,ATMASS)198:             CALL SHIFTH(QTS,.TRUE.,NOPT,NATOMS,ATMASS)
199:          ENDIF199:          ENDIF
200:       ENDIF200:       ENDIF
201:       IF (FREEZE.AND.VARIABLES) THEN 
202:          DO J1=1,NATOMS 
203:             IF (FROZEN(J1)) HESS(J1,J1)=SHIFTL(1) 
204:          ENDDO 
205:       ENDIF 
206:  
207:       IF (VARIABLES.AND.FOURDPBCT) CALL SHIFTFOURDPBCT201:       IF (VARIABLES.AND.FOURDPBCT) CALL SHIFTFOURDPBCT
208:       IF (VARIABLES.AND.THREEDPBCT) CALL SHIFTTHREEDPBCT202:       IF (VARIABLES.AND.THREEDPBCT) CALL SHIFTTHREEDPBCT
209:       IF (VARIABLES.AND.TWODPBCT) CALL SHIFTTWODPBCT203:       IF (VARIABLES.AND.TWODPBCT) CALL SHIFTTWODPBCT
210:       IF (VARIABLES.AND.ONEDPBCT) CALL SHIFTONEDPBCT204:       IF (VARIABLES.AND.ONEDPBCT) CALL SHIFTONEDPBCT
211:       IF (VARIABLES.AND.INVTONEDPBCT) CALL SHIFTONEDPBCT205:       IF (VARIABLES.AND.INVTONEDPBCT) CALL SHIFTONEDPBCT
212:       IF (VARIABLES.AND.INVTTWODPBCT) CALL SHIFTTWODPBCT206:       IF (VARIABLES.AND.INVTTWODPBCT) CALL SHIFTTWODPBCT
213: !      IF (VARIABLES.AND.EX1DT) CALL SHIFTEX1D207: !      IF (VARIABLES.AND.EX1DT) CALL SHIFTEX1D
214:       IF (ISTCRT.NE.10) THEN208:       IF (ISTCRT.NE.10) THEN
215: !        IF (PTEST) WRITE(*,'(A,4F20.10)') 'ENERGY,EOLD,ENERGY-EOLD,DELE=',ENERGY,EOLD,ENERGY-EOLD,DELE209: !        IF (PTEST) WRITE(*,'(A,4F20.10)') 'ENERGY,EOLD,ENERGY-EOLD,DELE=',ENERGY,EOLD,ENERGY-EOLD,DELE
216:          IF ((ENERGY-EOLD.NE.0.0D0).AND.(EOLD.NE.0.0D0)) THEN210:          IF ((ENERGY-EOLD.NE.0.0D0).AND.(EOLD.NE.0.0D0)) THEN
1169:       IF (MASST) THEN1163:       IF (MASST) THEN
1170:          DO J1=1,NATOMS1164:          DO J1=1,NATOMS
1171:             AMASS=1/SQRT(ATMASS(J1))1165:             AMASS=1/SQRT(ATMASS(J1))
1172:             QTS(3*J1-2)=AMASS*QTS(3*J1-2)1166:             QTS(3*J1-2)=AMASS*QTS(3*J1-2)
1173:             QTS(3*J1-1)=AMASS*QTS(3*J1-1)1167:             QTS(3*J1-1)=AMASS*QTS(3*J1-1)
1174:             QTS(3*J1)=AMASS*QTS(3*J1)1168:             QTS(3*J1)=AMASS*QTS(3*J1)
1175:          ENDDO1169:          ENDDO
1176:       ENDIF1170:       ENDIF
1177: 1171: 
1178:       CALL FLUSH(6)1172:       CALL FLUSH(6)
1179: ! What was gmetry doing?? DJW1173:       CALL GMETRY(ITER,VEC,QTS)
1180: ! 
1181: !     CALL GMETRY(ITER,VEC,QTS) 
1182: !1174: !
1183: !     Print out the coordinates and distance matrix1175: !     Print out the coordinates and distance matrix
1184: !1176: !
1185:       IF (ADMT.AND.(MOD(ITER-1,NADM).EQ.0)) CALL ADM(QTS)1177:       IF (ADMT.AND.(MOD(ITER-1,NADM).EQ.0)) CALL ADM(QTS)
1186:       EOLD=ENERGY1178:       EOLD=ENERGY
1187:       DO J1=1,NOPT1179:       DO J1=1,NOPT
1188:          PSTEP(J1)=STEP(J1)1180:          PSTEP(J1)=STEP(J1)
1189:          PFOB(J1)=FOB(J1)1181:          PFOB(J1)=FOB(J1)
1190:          PZT(J1)=ZT(J1)1182:          PZT(J1)=ZT(J1)
1191:       ENDDO1183:       ENDDO


r32316/fetchz.f 2017-04-22 17:30:12.778875415 +0100 r32315/fetchz.f 2017-04-22 17:30:17.122931533 +0100
 45:       INTEGER NTYPEA, NELEMENTS, NOXYGEN, NCHARGE 45:       INTEGER NTYPEA, NELEMENTS, NOXYGEN, NCHARGE
 46:       INTEGER, ALLOCATABLE :: ELEMENT_NUMS(:) 46:       INTEGER, ALLOCATABLE :: ELEMENT_NUMS(:)
 47:       COMMON /BIN/ EPSAB, EPSBB, SIGAB, SIGBB, NTYPEA 47:       COMMON /BIN/ EPSAB, EPSBB, SIGAB, SIGBB, NTYPEA
 48:       CHARACTER FNAME*80, TSTRING*80, CD1*1, CD2*2 48:       CHARACTER FNAME*80, TSTRING*80, CD1*1, CD2*2
 49:       COMMON /CAS/ AMAT, AINV, NELEMENTS, NTYPE 49:       COMMON /CAS/ AMAT, AINV, NELEMENTS, NTYPE
 50:       LOGICAL CONNECTT, DUMPPATH, READPATH, CALCRATES, STOPFIRST 50:       LOGICAL CONNECTT, DUMPPATH, READPATH, CALCRATES, STOPFIRST
 51:       INTEGER NCONNECT,I500,LUNIT,GETUNIT 51:       INTEGER NCONNECT,I500,LUNIT,GETUNIT
 52:       DOUBLE PRECISION TEMPERATURE, HRED, SCALE_FACT 52:       DOUBLE PRECISION TEMPERATURE, HRED, SCALE_FACT
 53:       COMMON /CONN/ STOPFIRST, CONNECTT, NCONNECT, DUMPPATH, READPATH, CALCRATES, TEMPERATURE, HRED 53:       COMMON /CONN/ STOPFIRST, CONNECTT, NCONNECT, DUMPPATH, READPATH, CALCRATES, TEMPERATURE, HRED
 54:       LOGICAL PATHT, DRAGT 54:       LOGICAL PATHT, DRAGT
 55:       INTEGER NPATHFRAME 55:       INTEGER NPATHFRAME, ISTAT
 56:       ! sn402: The name for this block appeared to be missing, so I've added it. 56:       ! sn402: The name for this block appeared to be missing, so I've added it.
 57:       COMMON /RUNTYPE/ DRAGT, PATHT, NPATHFRAME 57:       COMMON /RUNTYPE/ DRAGT, PATHT, NPATHFRAME
 58: C     DOUBLE PRECISION REPELTS(3*NATOMS,100), REPELPUSH 58: C     DOUBLE PRECISION REPELTS(3*NATOMS,100), REPELPUSH
 59: C     INTEGER NREPELTS, REPELFROM 59: C     INTEGER NREPELTS, REPELFROM
 60: C     LOGICAL REPELTST, REPEL 60: C     LOGICAL REPELTST, REPEL
 61: C     COMMON /OTS/ NREPELTS, REPELTST, REPELPUSH, REPEL, REPELFROM 61: C     COMMON /OTS/ NREPELTS, REPELTST, REPELPUSH, REPEL, REPELFROM
 62:       LOGICAL CUBIC 62:       LOGICAL CUBIC
 63:       COMMON /CUB/ CUBIC 63:       COMMON /CUB/ CUBIC
 64:       DOUBLE PRECISION CAPSRHO, EPS2, RAD, HEIGHT 64:       DOUBLE PRECISION CAPSRHO, EPS2, RAD, HEIGHT
 65:       COMMON /CAPS/ CAPSRHO, EPS2, RAD, HEIGHT 65:       COMMON /CAPS/ CAPSRHO, EPS2, RAD, HEIGHT
 66:       character(len=20) :: FTEMP 66:       character(len=20) :: str,str2,FTEMP,FINSTRING
 67:       character(len=20) :: FINSTRING 67:       CHARACTER(LEN=21) :: WW, RUBBISH
 68:       CHARACTER(LEN=21) :: WW 
 69:       CHARACTER(LEN=80) :: MYLINE 68:       CHARACTER(LEN=80) :: MYLINE
 70:       CHARACTER(5),ALLOCATABLE :: ELEMENT_NAMES(:)  69:       CHARACTER(5),ALLOCATABLE :: ELEMENT_NAMES(:) 
 71:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST 70:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST
 72:       LOGICAL END, SKIPBL, CLEAR, ECHO, CAT 71:       LOGICAL END, SKIPBL, CLEAR, ECHO, CAT
 73:       COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR, NERROR, IR, ECHO, LAST, CAT 72:       COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR, NERROR, IR, ECHO, LAST, CAT
 74:  73: 
 75:       CHARACTER(LEN=13) :: CUDAFILENAME1 74:       CHARACTER(LEN=13) :: CUDAFILENAME1
 76:       CHARACTER(LEN=34) :: CUDAFILENAME2 75:       CHARACTER(LEN=34) :: CUDAFILENAME2
 77:       CHARACTER(LEN=33) :: CUDAFILENAME3 76:       CHARACTER(LEN=33) :: CUDAFILENAME3
 78:       CHARACTER(LEN=20) :: CUDAFILENAME4 77:       CHARACTER(LEN=20) :: CUDAFILENAME4
1252: !           END DO1251: !           END DO
1253:          ENDIF1252:          ENDIF
1254: ! jwrm2> end1253: ! jwrm2> end
1255: 1254: 
1256: Cdavidg: changes to make lnatoms correct!1255: Cdavidg: changes to make lnatoms correct!
1257:          1256:          
1258:          IF (USERPOTT) THEN1257:          IF (USERPOTT) THEN
1259:             LNATOMS=NATOMS1258:             LNATOMS=NATOMS
1260:          ENDIF1259:          ENDIF
1261: 1260: 
1262:          LUNIT=GETUNIT() 
1263:          IF (TWOENDS.OR.CONNECTT.OR.NEBT.OR.NEWNEBT.OR.DRAGT.OR.GUESSPATHT1261:          IF (TWOENDS.OR.CONNECTT.OR.NEBT.OR.NEWNEBT.OR.DRAGT.OR.GUESSPATHT
1264:      $     .OR.MECCANOT.OR.MORPHT.OR.GREATCIRCLET.OR.GROWSTRINGT.OR.BHINTERPT.OR.BISECTT.OR.CLOSESTALIGNMENT) THEN1262:      $     .OR.MECCANOT.OR.MORPHT.OR.GREATCIRCLET.OR.GROWSTRINGT.OR.BHINTERPT.OR.BISECTT.OR.CLOSESTALIGNMENT) THEN
1265:             OPEN (UNIT=LUNIT,FILE=FINSTRING,STATUS='OLD') 1263:             OPEN (UNIT=7,FILE=FINSTRING,STATUS='OLD') 
1266: 1264: 
1267:             IF (GTHOMSONT) THEN1265:             IF (GTHOMSONT) THEN
1268: ! hk2861266: ! hk286
1269: !               DO J1=1,LNATOMS/2*31267: !               DO J1=1,LNATOMS/2*3
1270: !                  READ(LUNIT,*) FIN(3*(J1-1)+1),FIN(3*(J1-1)+2),FIN(3*(J1-1)+3)1268: !                  READ(7,*) FIN(3*(J1-1)+1),FIN(3*(J1-1)+2),FIN(3*(J1-1)+3)
1271: !               ENDDO1269: !               ENDDO
1272: !               CALL GTHOMSONCTOANG(FIN,THTEMP(1:3*NATOMS),NATOMS/2*3,1)1270: !               CALL GTHOMSONCTOANG(FIN,THTEMP(1:3*NATOMS),NATOMS/2*3,1)
1273: !               FIN(1:3*NATOMS)=THTEMP(1:3*NATOMS)1271: !               FIN(1:3*NATOMS)=THTEMP(1:3*NATOMS)
1274: 1272: 
1275: ! jwrm2> using TWOD framework for Gthomson coordinates1273: ! jwrm2> using TWOD framework for Gthomson coordinates
1276:                DO J1 = 1, NATOMS1274:                DO J1 = 1, NATOMS
1277:                   READ(LUNIT,*) FIN(3*(J1-1)+1), FIN(3*(J1-1)+2), FIN(3*(J1-1)+3)1275:                   READ(7,*) FIN(3*(J1-1)+1), FIN(3*(J1-1)+2), FIN(3*(J1-1)+3)
1278:                END DO1276:                END DO
1279:                CALL GTHOMSONCTOANG(FIN, THTEMP(1:3*NATOMS), NATOMS, 1)1277:                CALL GTHOMSONCTOANG(FIN, THTEMP(1:3*NATOMS), NATOMS, 1)
1280:                FIN(1:3*NATOMS) = THTEMP(1:3*NATOMS)1278:                FIN(1:3*NATOMS) = THTEMP(1:3*NATOMS)
1281: ! jwrm2> end1279: ! jwrm2> end
1282:             ELSE1280:             ELSE
1283:                DO J1=1,LNATOMS1281:                DO J1=1,LNATOMS
1284:                   READ(LUNIT,*) FIN(3*(J1-1)+1),FIN(3*(J1-1)+2),FIN(3*(J1-1)+3)1282:                   READ(7,*) FIN(3*(J1-1)+1),FIN(3*(J1-1)+2),FIN(3*(J1-1)+3)
1285:                ENDDO1283:                ENDDO
1286:             ENDIF1284:             ENDIF
1287:             CLOSE(LUNIT)1285:             CLOSE(7)
1288:             WRITE(*,'(A,I4,A,I4,A)') ' fetchz> Coordinates of second point read from file finish'1286:             WRITE(*,'(A,I4,A,I4,A)') ' fetchz> Coordinates of second point read from file finish'
1289: C           IF (ZSYM(1)(1:1).EQ.'W') THEN   !  WCOMMENT1287: C           IF (ZSYM(1)(1:1).EQ.'W') THEN   !  WCOMMENT
1290: C              OPEN (UNIT=7,FILE='Euler.finish',STATUS='OLD')1288: C              OPEN (UNIT=7,FILE='Euler.finish',STATUS='OLD')
1291: C              DO J1=1,LNATOMS1289: C              DO J1=1,LNATOMS
1292: C                 READ(7,*) FIN(3*LNATOMS+3*(J1-1)+1),FIN(3*LNATOMS+3*(J1-1)+2),FIN(3*LNATOMS+3*(J1-1)+3)1290: C                 READ(7,*) FIN(3*LNATOMS+3*(J1-1)+1),FIN(3*LNATOMS+3*(J1-1)+2),FIN(3*LNATOMS+3*(J1-1)+3)
1293: C              ENDDO1291: C              ENDDO
1294: C              CLOSE(7)1292: C              CLOSE(7)
1295: C              WRITE(*,'(A,I4,A,I4,A)') ' fetchz> Euler angles of second point read from file finish'1293: C              WRITE(*,'(A,I4,A,I4,A)') ' fetchz> Euler angles of second point read from file finish'
1296: C           ENDIF1294: C           ENDIF
1297:          ENDIF1295:          ENDIF
1679:          IF (T12FAC.LE.1.0D0) THEN1677:          IF (T12FAC.LE.1.0D0) THEN
1680:             WRITE(*,'(A,F12.4)') ' fetchz> Fixed initial uphill search direction, fraction of first collision time used=',T12FAC1678:             WRITE(*,'(A,F12.4)') ' fetchz> Fixed initial uphill search direction, fraction of first collision time used=',T12FAC
1681:          ELSE1679:          ELSE
1682:             WRITE(*,'(A,F12.4)') ' fetchz> Fixed initial uphill search direction, first colliding atoms will be placed half way'1680:             WRITE(*,'(A,F12.4)') ' fetchz> Fixed initial uphill search direction, first colliding atoms will be placed half way'
1683:          ENDIF1681:          ENDIF
1684:       ENDIF1682:       ENDIF
1685:       IF (READV) WRITE(*,'(A)') ' fetchz> Initial eigenvector will be read from vector.dump'1683:       IF (READV) WRITE(*,'(A)') ' fetchz> Initial eigenvector will be read from vector.dump'
1686:       IF (READSP) WRITE(*,'(A)') ' fetchz> stationary point info will be read'1684:       IF (READSP) WRITE(*,'(A)') ' fetchz> stationary point info will be read'
1687:       IF (DUMPSP) WRITE(*,'(A)') ' fetchz> stationary point info will be dumped'1685:       IF (DUMPSP) WRITE(*,'(A)') ' fetchz> stationary point info will be dumped'
1688:       IF (FREEZE) THEN1686:       IF (FREEZE) THEN
1689:          WRITE(*,'(A,I6,A)') ' fetchz> ', NFREEZE,' atoms or variables will be frozen:'1687:          WRITE(*,'(A,I6,A)') ' fetchz> ', NFREEZE,' atoms will be frozen:'
1690:          DO J1=1,NATOMS1688:          DO J1=1,NATOMS
1691:             IF (FROZEN(J1)) WRITE(*,'(I6)') J11689:             IF (FROZEN(J1)) WRITE(*,'(I6)') J1
1692:          ENDDO1690:          ENDDO
1693:       ENDIF1691:       ENDIF
1694: 1692: 
1695:       ! jbr36 classical rate calculations1693:       ! jbr36 classical rate calculations
1696:       IF (CLASSICALRATEST) THEN1694:       IF (CLASSICALRATEST) THEN
1697:          WRITE(*,'(A,F12.5,A,E20.10)') ' fetch> Classical rates and corrections will be calculated:'1695:          WRITE(*,'(A,F12.5,A,E20.10)') ' fetch> Classical rates and corrections will be calculated:'
1698:       ENDIF1696:       ENDIF
1699: 1697: 


r32316/geopt.f 2017-04-22 17:30:13.006878360 +0100 r32315/geopt.f 2017-04-22 17:30:17.354934528 +0100
 28:       USE MODNEB,ONLY : NEWCONNECTT 28:       USE MODNEB,ONLY : NEWCONNECTT
 29:       USE AMHGLOBALS, ONLY : NMRES 29:       USE AMHGLOBALS, ONLY : NMRES
 30: ! hk286 30: ! hk286
 31:       USE GENRIGID 31:       USE GENRIGID
 32:       use internals_wrapper 32:       use internals_wrapper
 33:       USE MODCUDALBFGS, ONLY : CUDA_LBFGS_WRAPPER 33:       USE MODCUDALBFGS, ONLY : CUDA_LBFGS_WRAPPER
 34:       USE MODCUDABFGSTS, ONLY : CUDA_BFGSTS_WRAPPER 34:       USE MODCUDABFGSTS, ONLY : CUDA_BFGSTS_WRAPPER
 35:       USE CHIRALITY, ONLY: CIS_TRANS_CHECK, CHIRALITY_CHECK 35:       USE CHIRALITY, ONLY: CIS_TRANS_CHECK, CHIRALITY_CHECK
 36:  36: 
 37:       IMPLICIT NONE 37:       IMPLICIT NONE
 38:       INTEGER FCALL, INEG, ECALL, NPCALL, SCALL, ITDONE, INFO, NEXMODES, II1, I1, NUMHB, J1, HORDER, J3, NSTEP, J2, K1, K2, 38:       INTEGER FCALL, INEG, ECALL, NPCALL, SCALL, ITDONE, INFO, NEXMODES, I1, NUMHB, J1, HORDER, J3, NSTEP, J2, K1, K2,
 39:      &        IMAX, IDONE, ITS  39:      &        IMAX, IDONE, ITS 
 40:       DOUBLE PRECISION  DIAG(3*NATOMS), ENERGY, EREAL, EVALMAX, EVALMIN, ETIME, FTIME, STIME, Q(3*NATOMS), OMAX, 40:       DOUBLE PRECISION  DIAG(3*NATOMS), ENERGY, EREAL, EVALMAX, EVALMIN, ETIME, FTIME, STIME, Q(3*NATOMS), OMAX,
 41:      1                  RMS, RMS2, VECS(3*NATOMS), VNEW(3*NATOMS), IT(3,3), IV(3,3), AMASS, ATMASSSAVE(NATOMS), DUMMY, 41:      1                  RMS, RMS2, VECS(3*NATOMS), VNEW(3*NATOMS), IT(3,3), IV(3,3), AMASS, ATMASSSAVE(NATOMS), DUMMY,
 42:      2                  EVALUES(3*NATOMS),TEMPA(9*NATOMS), DUMQ(3*NATOMS), ITX, ITY, ITZ, RMSD, RGYR, PROD, DPRAND, TSDISP, 42:      2                  EVALUES(3*NATOMS),TEMPA(9*NATOMS), DUMQ(3*NATOMS), ITX, ITY, ITZ, RMSD, RGYR, PROD, DPRAND, TSDISP,
 43:      &                  EVSAVE(3*NATOMS), ESAVE, TSDISPSAVE, MINFRQ2, MINCURVE, GRAD(3*NATOMS), RPBN, QFAC, TMPC(3*NATOMS) 43:      &                  EVSAVE(3*NATOMS), ESAVE, TSDISPSAVE, MINFRQ2, MINCURVE, GRAD(3*NATOMS), RPBN, QFAC, TMPC(3*NATOMS)
 44:       DOUBLE PRECISION DIFF, ORDERPLUS, ORDERMINUS, DIST, DIST2, RMAT(3,3), QPATH(3*NATOMS) 44:       DOUBLE PRECISION DIFF, ORDERPLUS, ORDERMINUS, DIST, DIST2, RMAT(3,3), QPATH(3*NATOMS)
 45:       DOUBLE PRECISION, ALLOCATABLE :: ORDERDERIV(:),ORDER(:),ORDERSAVE(:) 45:       DOUBLE PRECISION, ALLOCATABLE :: ORDERDERIV(:),ORDER(:),ORDERSAVE(:)
 46:       DOUBLE PRECISION QSAVE(3*NATOMS) 46:       DOUBLE PRECISION QSAVE(3*NATOMS)
 47:       CHARACTER(LEN=80) FNAMEF 47:       CHARACTER(LEN=80) FNAMEF
 48:       CHARACTER(LEN=20) EFNAMEF 48:       CHARACTER(LEN=20) EFNAMEF
229:                 ENDIF229:                 ENDIF
230:              ELSE IF (BSMIN.OR.RKMIN) THEN230:              ELSE IF (BSMIN.OR.RKMIN) THEN
231:                 CALL ODESD(NSTEPS,Q,MFLAG,ITDONE,.TRUE.)231:                 CALL ODESD(NSTEPS,Q,MFLAG,ITDONE,.TRUE.)
232:              ENDIF232:              ENDIF
233:           ENDIF233:           ENDIF
234: 234: 
235: C235: C
236: C Using xmylbfgs to locate the lowest eigenvalue and the corresponding eigenvector236: C Using xmylbfgs to locate the lowest eigenvalue and the corresponding eigenvector
237: C237: C
238:       ELSE IF (EIGENONLY) THEN238:       ELSE IF (EIGENONLY) THEN
239:          TSTART=TSTART-1.0d0239:       TSTART=TSTART-1.0d0
240:           PRINT '(A)', ' geopt> locate the lowest eigenvector with lbfgs'240:           PRINT '(A)', ' geopt> locate the lowest eigenvector with lbfgs'
241:           DO J1=1,NOPT241:           DO J1=1,NOPT
242:              VECS(J1)=DPRAND()*2-1.0D0242:              VECS(J1)=DPRAND()*2-1.0D0
243:           ENDDO243:           ENDDO
244:           IF (READV) THEN244:           IF (READV) THEN
245:              PRINT '(A)', ' geopt> read vector.dump as initial vector'245:              PRINT '(A)', ' geopt> read vector.dump as initial vector'
246:              IF(.NOT.DUMPV) OPEN(44,FILE="vector.dump",STATUS='UNKNOWN')246:              IF(.NOT.DUMPV) OPEN(44,FILE="vector.dump",STATUS='UNKNOWN')
247:              REWIND(44)247:              REWIND(44)
248:              READ(44,*,END=111) EVALMIN248:              READ(44,*,END=111) EVALMIN
249:              READ(44,*,END=111) (VECS(J1),J1=1,NOPT)249:              READ(44,*,END=111) (VECS(J1),J1=1,NOPT)
250:              IF(.NOT.DUMPV) CLOSE(44)250:              IF(.NOT.DUMPV) CLOSE(44)
251: 111          CONTINUE251: 111          CONTINUE
252:           ENDIF252:           ENDIF
253:           IF (FREEZE) THEN253:           IF (FREEZE) THEN
254:              DO J1=1,NATOMS254:              DO J1=1,NATOMS
255:                 IF (FROZEN(J1)) THEN255:                 IF (FROZEN(J1)) THEN
256:                    IF (VARIABLES) THEN256:                    VECS(3*(J1-1)+1)=0.0D0
257:                       VECS(J1)=0.0D0257:                    VECS(3*(J1-1)+2)=0.0D0
258:                    ELSE258:                    VECS(3*(J1-1)+3)=0.0D0
259:                       VECS(3*(J1-1)+1)=0.0D0 
260:                       VECS(3*(J1-1)+2)=0.0D0 
261:                       VECS(3*(J1-1)+3)=0.0D0 
262:                    ENDIF 
263:                 ENDIF259:                 ENDIF
264:              ENDDO260:              ENDDO
265:           ENDIF261:           ENDIF
266:           CALL VECNORM(VECS,NOPT)262:           CALL VECNORM(VECS,NOPT)
267:           SOVER=0.0d0263:           SOVER=0.0d0
268:           CALL BEIG(ITDONE,Q,ENERGY,VECS,EVALMIN,J1,SOVER,.TRUE.,MFLAG)264:           CALL BEIG(ITDONE,Q,ENERGY,VECS,EVALMIN,J1,SOVER,.TRUE.,MFLAG)
269:           IF (DUMPV) THEN265:           IF (DUMPV) THEN
270:             REWIND(44)266:             REWIND(44)
271:             WRITE(44,'(E20.10)') EVALMIN267:             WRITE(44,'(E20.10)') EVALMIN
272:             WRITE(44,'(3F20.10)') (VECS(J1),J1=1,NOPT)268:             WRITE(44,'(3F20.10)') (VECS(J1),J1=1,NOPT)
1276: ! are written correctly to points.final etc.1272: ! are written correctly to points.final etc.
1277: ! hk2861273: ! hk286
1278:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN1274:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN
1279:          ! DUMQ will hold the cartesian coordinates. We don't change ATOMRIGIDCOORDT because we're only changing1275:          ! DUMQ will hold the cartesian coordinates. We don't change ATOMRIGIDCOORDT because we're only changing
1280:          ! coordinates of a dummy array, not Q itself.1276:          ! coordinates of a dummy array, not Q itself.
1281:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, DUMQ, Q(1:DEGFREEDOMS))1277:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, DUMQ, Q(1:DEGFREEDOMS))
1282:       ELSE IF (GTHOMSONT) THEN1278:       ELSE IF (GTHOMSONT) THEN
1283: !jwrm2> With GTHOMSON, want Cartesian coordinates for moment of inertia calculation1279: !jwrm2> With GTHOMSON, want Cartesian coordinates for moment of inertia calculation
1284:          CALL GTHOMSONANGTOC(DUMQ(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)1280:          CALL GTHOMSONANGTOC(DUMQ(1:3*NATOMS), Q(1:3*NATOMS), NATOMS)
1285:       ELSE1281:       ELSE
1286:          DO II1=1,NOPT1282:          DO I1=1,3*NATOMS
1287:             DUMQ(II1)=Q(II1)1283:             DUMQ(I1)=Q(I1)
1288:          ENDDO1284:          ENDDO
1289:       ENDIF1285:       ENDIF
1290: 1286: 
1291:       INERTIAT=.TRUE.1287:       INERTIAT=.TRUE.
1292: C     IF ((UNRST.OR.CHRMMT).AND.INERTIAT) THEN1288: C     IF ((UNRST.OR.CHRMMT).AND.INERTIAT) THEN
1293: 1289: 
1294:       ! The next block is used for dumping to a min.data.info file. The frequencies calculated1290:       ! The next block is used for dumping to a min.data.info file. The frequencies calculated
1295:       ! previously will be used, unless NOFRQS is set.1291:       ! previously will be used, unless NOFRQS is set.
1296:       IF ((DUMPDATAT).AND.INERTIAT) THEN1292:       IF ((DUMPDATAT).AND.INERTIAT) THEN
1297: C1293: C
1298: C  Pathsample uses unit masses, so we need to change to unit mass temporarily here,1294: C  Pathsample uses unit masses, so we need to change to unit mass temporarily here,
1299: C  and then put the right values back, otherwise the frequencies will be wrong1295: C  and then put the right values back, otherwise the frequencies will be wrong
1300: C  when we call geopt again from newconnect in a bhinterp run!!!!1296: C  when we call geopt again from newconnect in a bhinterp run!!!!
1301: C1297: C
1302:          DO J1=1,NATOMS1298:          DO J1=1,NATOMS
1303:             ATMASSSAVE(J1)=ATMASS(J1)1299:             ATMASSSAVE(J1)=ATMASS(J1)
1304:             ATMASS(J1)=1.0D0 1300:             ATMASS(J1)=1.0D0 
1305:          ENDDO1301:          ENDDO
1306: 1302: 
1307:          ! Get the inertia tensor1303:          ! Get the inertia tensor
1308:          IF (VARIABLES) THEN1304:          CALL INERTIA2(DUMQ,ITX,ITY,ITZ)
1309:             ITX=1.0D0; ITY=1.0D0; ITZ=1.0D0 
1310:          ELSE 
1311:             CALL INERTIA2(DUMQ,ITX,ITY,ITZ) 
1312:          ENDIF 
1313:          WRITE(*,'(A,4X,F20.10)') ' geopt> X component of inertia tensor=',ITX1305:          WRITE(*,'(A,4X,F20.10)') ' geopt> X component of inertia tensor=',ITX
1314:          WRITE(*,'(A,4X,F20.10)') ' geopt> Y component of inertia tensor=',ITY1306:          WRITE(*,'(A,4X,F20.10)') ' geopt> Y component of inertia tensor=',ITY
1315:          WRITE(*,'(A,4X,F20.10)') ' geopt> Z component of inertia tensor=',ITZ1307:          WRITE(*,'(A,4X,F20.10)') ' geopt> Z component of inertia tensor=',ITZ
1316:          HORDER=11308:          HORDER=1
1317:          DO J1=1,NATOMS1309:          DO J1=1,NATOMS
1318:             ATMASS(J1)=ATMASSSAVE(J1)1310:             ATMASS(J1)=ATMASSSAVE(J1)
1319:          ENDDO1311:          ENDDO
1320:       ENDIF1312:       ENDIF
1321: 1313: 
1322: ! jwrm2> With GTHOMSON, coordinates are expected to be polar for symmetry1314: ! jwrm2> With GTHOMSON, coordinates are expected to be polar for symmetry
1717:       ENDDO1709:       ENDDO
1718: 1710: 
1719:       RETURN1711:       RETURN
1720:       END1712:       END
1721: !1713: !
1722: ! Somewhat better eigenvalues for twice the cpu time. Uses Richardson extrapolation.1714: ! Somewhat better eigenvalues for twice the cpu time. Uses Richardson extrapolation.
1723: !1715: !
1724:       SUBROUTINE MAKENUMHESS2(X,NATOMS)1716:       SUBROUTINE MAKENUMHESS2(X,NATOMS)
1725:       USE MODHESS1717:       USE MODHESS
1726:       USE MODCHARMM1718:       USE MODCHARMM
1727:       USE KEY,ONLY : FROZEN, CASTEP, AMHT, ONETEP, ENDNUMHESSDELTA, VARIABLES1719:       USE KEY,ONLY : FROZEN, CASTEP, AMHT, ONETEP, ENDNUMHESSDELTA
1728:       USE COMMONS,ONLY : DEBUG, NOPT1720:       USE COMMONS,ONLY : DEBUG, NOPT
1729:       use porfuncs1721:       use porfuncs
1730:       IMPLICIT NONE1722:       IMPLICIT NONE
1731:       LOGICAL KNOWE, KNOWG, KNOWH1723:       LOGICAL KNOWE, KNOWG, KNOWH
1732:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH1724:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
1733: 1725: 
1734:       INTEGER I1,J1,NATOMS,ISTAT1726:       INTEGER I1,J1,NATOMS,ISTAT
1735:       DOUBLE PRECISION X(NOPT)1727:       DOUBLE PRECISION X(NOPT)
1736:       DOUBLE PRECISION DUM(NOPT),GRAD1(NOPT),GRAD2(NOPT),DELTA,RMS,ENERGY1728:       DOUBLE PRECISION DUM(NOPT),GRAD1(NOPT),GRAD2(NOPT),DELTA,RMS,ENERGY
1737:       DOUBLE PRECISION GRAD1B(NOPT),GRAD2B(NOPT)1729:       DOUBLE PRECISION GRAD1B(NOPT),GRAD2B(NOPT)
1742:       ENDDO1734:       ENDDO
1743: 1735: 
1744:       DELTA=ENDNUMHESSDELTA1736:       DELTA=ENDNUMHESSDELTA
1745:       PRINT '(A,G20.10)',' makenumhess2> displacement=',DELTA1737:       PRINT '(A,G20.10)',' makenumhess2> displacement=',DELTA
1746: C1738: C
1747:       IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-11739:       IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-1
1748: 1740: 
1749:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))1741:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))
1750:       DO I1=1,NOPT1742:       DO I1=1,NOPT
1751:          CALL FLUSH(6)1743:          CALL FLUSH(6)
1752:          IF (VARIABLES.AND.FROZEN(I1)) THEN1744:          IF (FROZEN((I1-1)/3+1)) THEN
1753:             DO J1=I1,NOPT 
1754:                HESS(I1,J1)=0.0D0 
1755:                HESS(J1,I1)=0.0D0 
1756:             ENDDO 
1757:          ELSE IF ((.NOT.VARIABLES).AND.FROZEN((I1-1)/3+1)) THEN 
1758:             DO J1=I1,NOPT1745:             DO J1=I1,NOPT
1759:                HESS(I1,J1)=0.0D01746:                HESS(I1,J1)=0.0D0
1760:                HESS(J1,I1)=0.0D01747:                HESS(J1,I1)=0.0D0
1761:             ENDDO1748:             ENDDO
1762:          ELSE1749:          ELSE
1763:             DUM(I1)=X(I1)-DELTA1750:             DUM(I1)=X(I1)-DELTA
1764:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1751:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1765:             DUM(I1)=X(I1)-2.0D0*DELTA1752:             DUM(I1)=X(I1)-2.0D0*DELTA
1766:             CALL POTENTIAL(DUM,ENERGY,GRAD1B,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1753:             CALL POTENTIAL(DUM,ENERGY,GRAD1B,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1767:             DUM(I1)=X(I1)+DELTA1754:             DUM(I1)=X(I1)+DELTA
1781: 1768: 
1782:       RETURN1769:       RETURN
1783:       END1770:       END
1784: 1771: 
1785:       SUBROUTINE MAKENUMHESS(X,NATOMS)1772:       SUBROUTINE MAKENUMHESS(X,NATOMS)
1786: C1773: C
1787: C dae1774: C dae
1788: C1775: C
1789:       USE MODHESS1776:       USE MODHESS
1790:       USE MODCHARMM1777:       USE MODCHARMM
1791:       USE KEY,ONLY : FROZEN, CASTEP, AMHT, ONETEP, ENDNUMHESS2, QCHEM, VASP, VARIABLES1778:       USE KEY,ONLY : FROZEN, CASTEP, AMHT, ONETEP, ENDNUMHESS2, QCHEM, VASP
1792:       USE SBMDATA,ONLY : SBMHESSATOM,SBMFIRSTDD1779:       USE SBMDATA,ONLY : SBMHESSATOM,SBMFIRSTDD
1793:       USE COMMONS,ONLY : DEBUG, NOPT1780:       USE COMMONS,ONLY : DEBUG, NOPT
1794:       use porfuncs1781:       use porfuncs
1795:       IMPLICIT NONE1782:       IMPLICIT NONE
1796:       LOGICAL KNOWE, KNOWG, KNOWH1783:       LOGICAL KNOWE, KNOWG, KNOWH
1797:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH1784:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
1798: 1785: 
1799:       INTEGER I1,J1,NATOMS,ISTAT1786:       INTEGER I1,J1,NATOMS,ISTAT
1800:       DOUBLE PRECISION X(NOPT)1787:       DOUBLE PRECISION X(NOPT)
1801:       DOUBLE PRECISION DUM(NOPT),GRAD1(NOPT),GRAD2(NOPT),DELTA,RMS,ENERGY1788:       DOUBLE PRECISION DUM(NOPT),GRAD1(NOPT),GRAD2(NOPT),DELTA,RMS,ENERGY
1824: C compared to the analytical solution significantly (around 4sf), relative to the alternative1811: C compared to the analytical solution significantly (around 4sf), relative to the alternative
1825: C of evaluating the gradient once at the beginning then once for each element. But obviously it also1812: C of evaluating the gradient once at the beginning then once for each element. But obviously it also
1826: C doubles the number of potential evaluations ((3N)^2/2 for N atoms), so may be not worth it for1813: C doubles the number of potential evaluations ((3N)^2/2 for N atoms), so may be not worth it for
1827: C a long run on a big system.1814: C a long run on a big system.
1828: C1815: C
1829:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))1816:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))
1830: !      SBMFIRSTDD = .TRUE.1817: !      SBMFIRSTDD = .TRUE.
1831:       DO I1=1,NOPT1818:       DO I1=1,NOPT
1832: !        write(*,*) 'working on atom', I11819: !        write(*,*) 'working on atom', I1
1833:          CALL FLUSH(6)1820:          CALL FLUSH(6)
1834:          IF (VARIABLES.AND.FROZEN(I1)) THEN1821:          IF (FROZEN((I1-1)/3+1)) THEN
1835:             DO J1=I1,NOPT 
1836:                HESS(I1,J1)=0.0D0 
1837:                HESS(J1,I1)=0.0D0 
1838:             ENDDO 
1839:          ELSE IF ((.NOT.VARIABLES).AND.FROZEN((I1-1)/3+1)) THEN 
1840:             DO J1=I1,NOPT1822:             DO J1=I1,NOPT
1841:                HESS(I1,J1)=0.0D01823:                HESS(I1,J1)=0.0D0
1842:                HESS(J1,I1)=0.0D01824:                HESS(J1,I1)=0.0D0
1843:             ENDDO1825:             ENDDO
1844:          ELSE1826:          ELSE
1845:             DUM(I1)=X(I1)-DELTA1827:             DUM(I1)=X(I1)-DELTA
1846: !            SBMHESSATOM=int(I1/3)1828: !            SBMHESSATOM=int(I1/3)
1847:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1829:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1848:             DUM(I1)=X(I1)+DELTA1830:             DUM(I1)=X(I1)+DELTA
1849:             CALL POTENTIAL(DUM,ENERGY,GRAD2,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1831:             CALL POTENTIAL(DUM,ENERGY,GRAD2,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)


r32316/grad.f90 2017-04-22 17:30:10.654847975 +0100 r32315/grad.f90 2017-04-22 17:30:14.822901822 +0100
 29:           USE NEBDATA 29:           USE NEBDATA
 30:           USE KEYNEB,ONLY: NIMAGE 30:           USE KEYNEB,ONLY: NIMAGE
 31:           USE TANGENT 31:           USE TANGENT
 32:           USE SPFUNCTS 32:           USE SPFUNCTS
 33:           USE NEBUTILS 33:           USE NEBUTILS
 34:           USE MODUNRES, ONLY: C,NRES 34:           USE MODUNRES, ONLY: C,NRES
 35:           USE GENRIGID, ONLY: RIGIDINIT, DEGFREEDOMS, RB_DISTANCE 35:           USE GENRIGID, ONLY: RIGIDINIT, DEGFREEDOMS, RB_DISTANCE
 36:           USE KEY, ONLY: UNRST, FROZEN, FREEZE, NEBK, BULKT, NEBRESEEDT, NEBRESEEDEMAX, NEBRESEEDBMAX, & 36:           USE KEY, ONLY: UNRST, FROZEN, FREEZE, NEBK, BULKT, NEBRESEEDT, NEBRESEEDEMAX, NEBRESEEDBMAX, &
 37:   &                      RBAAT, NREPULSIVE, NEBRESEEDINT, INTCONFRAC, INTCONSTRAINTT, & 37:   &                      RBAAT, NREPULSIVE, NEBRESEEDINT, INTCONFRAC, INTCONSTRAINTT, &
 38:   &                      ORDERI,ORDERJ,EPSALPHA, REPPOW,PERMDIST, TWOD, RIGIDBODY, NEBRESEEDDEL1, & 38:   &                      ORDERI,ORDERJ,EPSALPHA, REPPOW,PERMDIST, TWOD, RIGIDBODY, NEBRESEEDDEL1, &
 39:   &                      DISTREF, ADDREPT, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, CHECKCONINT, VARIABLES 39:   &                      DISTREF, ADDREPT, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, CHECKCONINT
 40:           USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3, DEBUG 40:           USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3, DEBUG
 41:           IMPLICIT NONE 41:           IMPLICIT NONE
 42:             42:            
 43:           INTEGER :: J1,J2,K,J3,J4,J5,JOUT,JIN,NSHRINK,NSHRINKATOM(NATOMS),J6,NATTACH,NATTACHS,NMAXINT,NMININT,LUNIT 43:           INTEGER :: J1,J2,K,J3,J4,J5,JOUT,JIN,NSHRINK,NSHRINKATOM(NATOMS),J6,NATTACH,NATTACHS,NMAXINT,NMININT,LUNIT
 44:           DOUBLE PRECISION :: PERPCOMP=0.1D0 44:           DOUBLE PRECISION :: PERPCOMP=0.1D0
 45:           DOUBLE PRECISION :: DIHEDIST,TMPRMS,QINT(NINTS*(NIMAGE+2)),DIFFM(NINTS),DIFFP(NINTS) ! JMC 45:           DOUBLE PRECISION :: DIHEDIST,TMPRMS,QINT(NINTS*(NIMAGE+2)),DIFFM(NINTS),DIFFP(NINTS) ! JMC
 46:           DOUBLE PRECISION DIST 46:           DOUBLE PRECISION DIST
 47:           DOUBLE PRECISION :: PI=3.141592653589793D0 47:           DOUBLE PRECISION :: PI=3.141592653589793D0
 48:           DOUBLE PRECISION TEMP1(NOPT), TEMP2(NOPT), GGGSAVE(NOPT*(NIMAGE+2)), DUMMY, DUMMY2, SPRING(NOPT*(NIMAGE+2)), & 48:           DOUBLE PRECISION TEMP1(NOPT), TEMP2(NOPT), GGGSAVE(NOPT*(NIMAGE+2)), DUMMY, DUMMY2, SPRING(NOPT*(NIMAGE+2)), &
 49:   &                        TRUEPERP(NOPT*(NIMAGE+2)), IMCOORDS(NOPT), RSITE(NOPT), RMAT(3,3), D, DIST2, LX(3), LV(3), & 49:   &                        TRUEPERP(NOPT*(NIMAGE+2)), IMCOORDS(NOPT), RSITE(NOPT), RMAT(3,3), D, DIST2, LX(3), LV(3), &
576:                 ENDDO576:                 ENDDO
577:              ENDDO577:              ENDDO
578:           ENDIF578:           ENDIF
579: !579: !
580: ! Set gradients on frozen atoms to zero.580: ! Set gradients on frozen atoms to zero.
581: !581: !
582: 555       IF (FREEZE) THEN582: 555       IF (FREEZE) THEN
583:              DO J1=2,NIMAGE+1  583:              DO J1=2,NIMAGE+1  
584:                 DO J2=1,NATOMS584:                 DO J2=1,NATOMS
585:                    IF (FROZEN(J2)) THEN585:                    IF (FROZEN(J2)) THEN
586:                       IF (VARIABLES) THEN586:                       GGG(NOPT*(J1-1)+3*(J2-1)+1)=0.0D0
587:                          GGG(NOPT*(J1-1)+J2)=0.0D0587:                       GGG(NOPT*(J1-1)+3*(J2-1)+2)=0.0D0
588:                       ELSE588:                       GGG(NOPT*(J1-1)+3*(J2-1)+3)=0.0D0
589:                          GGG(NOPT*(J1-1)+3*(J2-1)+1)=0.0D0 
590:                          GGG(NOPT*(J1-1)+3*(J2-1)+2)=0.0D0 
591:                          GGG(NOPT*(J1-1)+3*(J2-1)+3)=0.0D0 
592:                       ENDIF 
593:                    ENDIF589:                    ENDIF
594:                 ENDDO590:                 ENDDO
595:              ENDDO591:              ENDDO
596:           ENDIF592:           ENDIF
597:           RMS=0.0D0593:           RMS=0.0D0
598:           DO J1=2,NIMAGE+1594:           DO J1=2,NIMAGE+1
599:              DO J2=1,NOPT595:              DO J2=1,NOPT
600:                 RMS=RMS+GGG(NOPT*(J1-1)+J2)**2596:                 RMS=RMS+GGG(NOPT*(J1-1)+J2)**2
601:              ENDDO597:              ENDDO
602:           ENDDO598:           ENDDO


r32316/hybridmin.f 2017-04-22 17:30:13.226881203 +0100 r32315/hybridmin.f 2017-04-22 17:30:17.578937423 +0100
111: 10          CLOSE(45)111: 10          CLOSE(45)
112:             DIAG(1)=EVALMIN112:             DIAG(1)=EVALMIN
113:          ENDIF113:          ENDIF
114:          WRITE(*,'(A,F20.10)') ' hybridmin> Reaction vector read successfully. Eigenvalue=   ',EVALMIN114:          WRITE(*,'(A,F20.10)') ' hybridmin> Reaction vector read successfully. Eigenvalue=   ',EVALMIN
115:          NS=100115:          NS=100
116:          READV=.FALSE. ! needed for runs that do bfgsts then path116:          READV=.FALSE. ! needed for runs that do bfgsts then path
117:       ENDIF117:       ENDIF
118:       IF (FREEZE) THEN118:       IF (FREEZE) THEN
119:          DO J1=1,NATOMS119:          DO J1=1,NATOMS
120:             IF (FROZEN(J1)) THEN120:             IF (FROZEN(J1)) THEN
121:                IF (VARIABLES) THEN121:                VECS(3*(J1-1)+1)=0.0D0
122:                   VECS(J1)=0.0D0122:                VECS(3*(J1-1)+2)=0.0D0
123:                ELSE123:                VECS(3*(J1-1)+3)=0.0D0
124:                   VECS(3*(J1-1)+1)=0.0D0 
125:                   VECS(3*(J1-1)+2)=0.0D0 
126:                   VECS(3*(J1-1)+3)=0.0D0 
127:                ENDIF 
128:             ENDIF124:             ENDIF
129:          ENDDO125:          ENDDO
130:       ENDIF126:       ENDIF
131: 127: 
132:       STEST=.TRUE.128:       STEST=.TRUE.
133:       KNOWG=.FALSE.129:       KNOWG=.FALSE.
134:       IF (NOHESS) STEST=.FALSE.130:       IF (NOHESS) STEST=.FALSE.
135:       IF (READV.AND.(ITER.EQ.1).AND.(NUSEEV.LE.1)) STEST=.FALSE.131:       IF (READV.AND.(ITER.EQ.1).AND.(NUSEEV.LE.1)) STEST=.FALSE.
136: 90    IF (PTEST) PRINT*132: 90    IF (PTEST) PRINT*
137:       NUP=NUSEEV133:       NUP=NUSEEV


r32316/keywords.f 2017-04-22 17:30:13.470884355 +0100 r32315/keywords.f 2017-04-22 17:30:17.798940263 +0100
6668: ! 6668: ! 
6669:          ELSE IF (WORD .EQ. 'VALUES') THEN6669:          ELSE IF (WORD .EQ. 'VALUES') THEN
6670:             CALL READI(NVALUES)6670:             CALL READI(NVALUES)
6671: ! 6671: ! 
6672: ! VARIABLES - keyword at the end of the list of options after which6672: ! VARIABLES - keyword at the end of the list of options after which
6673: ! the general variables follow. NZERO is the number of zero6673: ! the general variables follow. NZERO is the number of zero
6674: ! eigenvalues, default 0.6674: ! eigenvalues, default 0.
6675: ! 6675: ! 
6676:          ELSE IF (WORD.EQ.'VARIABLES') THEN6676:          ELSE IF (WORD.EQ.'VARIABLES') THEN
6677:             VARIABLES=.TRUE.6677:             VARIABLES=.TRUE.
6678:             ZSYM(1:NATOMS)=' ' ! otherwise valgrind reports uninitialised 
6679: 6678: 
6680:             IF (DUMMY_FRQCONV.EQ.0) THEN 6679:             IF (DUMMY_FRQCONV.EQ.0) THEN 
6681:                FRQCONV = 1.0D06680:                FRQCONV = 1.0D0
6682:             ELSE6681:             ELSE
6683:                FRQCONV = DUMMY_FRQCONV6682:                FRQCONV = DUMMY_FRQCONV
6684:             ENDIF6683:             ENDIF
6685:             FRQCONV2 = FRQCONV*FRQCONV6684:             FRQCONV2 = FRQCONV*FRQCONV
6686:             IF (NITEMS.GT.1) CALL READI(NZERO) 
6687: 6685: 
6688:             RETURN6686:             RETURN
6689: ! 6687: ! 
6690: ! 6688: ! 
6691: ! VASP tells the program to read derivative information in6689: ! VASP tells the program to read derivative information in
6692: ! VASP format.                                        - default FALSE6690: ! VASP format.                                        - default FALSE
6693: ! 6691: ! 
6694:          ELSE IF (WORD.EQ.'VASP') THEN6692:          ELSE IF (WORD.EQ.'VASP') THEN
6695:             VASP=.TRUE.6693:             VASP=.TRUE.
6696:             BULKT=.TRUE.6694:             BULKT=.TRUE.


r32316/lbfgs.f90 2017-04-22 17:30:10.874850820 +0100 r32315/lbfgs.f90 2017-04-22 17:30:15.058904866 +0100
 83:      INTPTEST = .FALSE. 83:      INTPTEST = .FALSE.
 84:      IF (DESMDEBUG) MOREPRINTING = .TRUE. 84:      IF (DESMDEBUG) MOREPRINTING = .TRUE.
 85:      ADDATOM=.FALSE. 85:      ADDATOM=.FALSE.
 86:      NFAIL=0 86:      NFAIL=0
 87:      IF (FREEZENODEST) IMGFREEZE(1:NIMAGE)=.FALSE. 87:      IF (FREEZENODEST) IMGFREEZE(1:NIMAGE)=.FALSE.
 88:      BESTE=HUGE(1.0D0) 88:      BESTE=HUGE(1.0D0)
 89:  89: 
 90:      CALL CHECKINPUT 90:      CALL CHECKINPUT
 91:      IF (DEBUG) CALL DUMPFILES("b") 91:      IF (DEBUG) CALL DUMPFILES("b")
 92:      IF (DUMPNEBXYZ) CALL RWG("w",.False.,0) 92:      IF (DUMPNEBXYZ) CALL RWG("w",.False.,0)
 93:      IF (DUMPNEBEOS) DVEC(1:NIMAGE+1)=0.0D0 ! otherwise uninitialised 
 94:      IF (DUMPNEBEOS) CALL WRITEPROFILE(0) 93:      IF (DUMPNEBEOS) CALL WRITEPROFILE(0)
 95:      IF (MOREPRINTING) THEN 94:      IF (MOREPRINTING) THEN
 96:           WRITE(*,'(A)') '                  Iter   Energy per image      RMS Force       Av.Dev    Path     Step per image' 95:           WRITE(*,'(A)') '                  Iter   Energy per image      RMS Force       Av.Dev    Path     Step per image'
 97:      ENDIF 96:      ENDIF
 98:      IF (PREVGRAD.LT.DNEBSWITCH) THEN 97:      IF (PREVGRAD.LT.DNEBSWITCH) THEN
 99:         CALL OLDNEBGRADIENT 98:         CALL OLDNEBGRADIENT
100:      ELSE 99:      ELSE
101:         CALL NEBGRADIENT100:         CALL NEBGRADIENT
102:      ENDIF101:      ENDIF
103: 102: 


r32316/mylbfgs.f 2017-04-22 17:30:13.694887240 +0100 r32315/mylbfgs.f 2017-04-22 17:30:18.030943262 +0100
 52:       USE GENRIGID 52:       USE GENRIGID
 53:       IMPLICIT NONE 53:       IMPLICIT NONE
 54:       INTEGER N,M,J1,J2,ITMAX,ITDONE,NFAIL,NCOUNT 54:       INTEGER N,M,J1,J2,ITMAX,ITDONE,NFAIL,NCOUNT
 55: C     DOUBLE PRECISION X(*),G(3*NATOMS),DIAG(N),W(N*(2*M+1)+2*M),SLENGTH,DDOT,VEC2(3*NATOMS),OVERLAP 55: C     DOUBLE PRECISION X(*),G(3*NATOMS),DIAG(N),W(N*(2*M+1)+2*M),SLENGTH,DDOT,VEC2(3*NATOMS),OVERLAP
 56:       DOUBLE PRECISION X(3*NATOMS),SLENGTH,DDOT,VEC2(3*NATOMS),OVERLAP,DISTF,DISTS,GAMMA 56:       DOUBLE PRECISION X(3*NATOMS),SLENGTH,DDOT,VEC2(3*NATOMS),OVERLAP,DISTF,DISTS,GAMMA
 57:       DOUBLE PRECISION, TARGET :: G(3*NATOMS) 57:       DOUBLE PRECISION, TARGET :: G(3*NATOMS)
 58: ! 58: !
 59: ! It is necessary to dynamically allocate DIAG and W for SAVE to work with pgi and ifort compilers. 59: ! It is necessary to dynamically allocate DIAG and W for SAVE to work with pgi and ifort compilers.
 60: ! 60: !
 61:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DIAG, W 61:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DIAG, W
 62:       DOUBLE PRECISION DUMMY1,ENERGY,ENEW,RMS,EREAL,REALRMS,RVEC(3*NATOMS),ALPHA,GSAVE(N),DUMMY2 62:       DOUBLE PRECISION DUMMY1,ENERGY,ENEW,RMS,EREAL,REALRMS,RVEC(3*NATOMS),ALPHA,GSAVE(3*NATOMS),DUMMY2
 63:       LOGICAL DIAGCO, RESET, PTEST, NODUMP, PROJECT 63:       LOGICAL DIAGCO, RESET, PTEST, NODUMP, PROJECT
 64:       DOUBLE PRECISION GNORM,STP,YS,YY,SQ,YR,BETA 64:       DOUBLE PRECISION GNORM,STP,YS,YY,SQ,YR,BETA
 65:       LOGICAL NOCOOR, NODERV ! internals stuff 65:       LOGICAL NOCOOR, NODERV ! internals stuff
 66:       DOUBLE PRECISION DELTACART(3*NATOMS) 66:       DOUBLE PRECISION DELTACART(3*NATOMS)
 67:       DOUBLE PRECISION CART(3*NATOMS),OLDQ(N),NEWQ(N),OLDGINT(N) 67:       DOUBLE PRECISION CART(3*NATOMS),OLDQ(N),NEWQ(N),OLDGINT(N)
 68:       DOUBLE PRECISION, TARGET :: GINT(N),XINT(N), DELTAQ(N) 68:       DOUBLE PRECISION, TARGET :: GINT(N),XINT(N), DELTAQ(N)
 69:       DOUBLE PRECISION OLDCART(3*NATOMS),TMPINT(NINTS) ! JMC 69:       DOUBLE PRECISION OLDCART(3*NATOMS),TMPINT(NINTS) ! JMC
 70:       DOUBLE PRECISION GLAST(N),XSAVE(N) 70:       DOUBLE PRECISION GLAST(3*NATOMS),XSAVE(N)
 71:       INTEGER ITER,POINT,ISPT,IYPT,BOUND,NPT,CP,I,INMC,IYCN,ISCN,I1,NREV 71:       INTEGER ITER,POINT,ISPT,IYPT,BOUND,NPT,CP,I,INMC,IYCN,ISCN,I1,NREV
 72:       LOGICAL MFLAG,FAILED 72:       LOGICAL MFLAG,FAILED
 73:       DOUBLE PRECISION GSQSCALE, GSTHRESH, DOT1, DOT2 73:       DOUBLE PRECISION GSQSCALE, GSTHRESH, DOT1, DOT2
 74:       INTEGER NSPECIAL, NALLOW, NINFO, FRAME 74:       INTEGER NSPECIAL, NALLOW, NINFO, FRAME
 75:       LOGICAL PVFLAG 75:       LOGICAL PVFLAG
 76:       COMMON /PVF/ PVFLAG 76:       COMMON /PVF/ PVFLAG
 77:       COMMON /G2/ GSTHRESH, GSQSCALE, NSPECIAL, NALLOW, NINFO 77:       COMMON /G2/ GSTHRESH, GSQSCALE, NSPECIAL, NALLOW, NINFO
 78:       LOGICAL PATHT, DRAGT 78:       LOGICAL PATHT, DRAGT
 79:       INTEGER NPATHFRAME, NDECREASE, ISTAT 79:       INTEGER NPATHFRAME, NDECREASE, ISTAT
 80:       COMMON /RUNTYPE/ DRAGT, PATHT, NPATHFRAME 80:       COMMON /RUNTYPE/ DRAGT, PATHT, NPATHFRAME
271:                   GSAVE(J1)=GSAVE(J1)-2.0D0*ZWORK(J1,J2)*DUMMY1271:                   GSAVE(J1)=GSAVE(J1)-2.0D0*ZWORK(J1,J2)*DUMMY1
272:                   RMS=RMS+GSAVE(J1)**2272:                   RMS=RMS+GSAVE(J1)**2
273:                ENDDO273:                ENDDO
274:                RMS=SQRT(RMS/N)274:                RMS=SQRT(RMS/N)
275:             ENDDO275:             ENDDO
276:          ENDIF276:          ENDIF
277:       ELSE277:       ELSE
278:          RMS=REALRMS278:          RMS=REALRMS
279:       ENDIF279:       ENDIF
280: 280: 
281:       G(1:N)=GSAVE(1:N)281:       G(1:3*NATOMS)=GSAVE(1:3*NATOMS)
282:       GLAST(1:N)=GSAVE(1:N)282:       GLAST(1:3*NATOMS)=GSAVE(1:3*NATOMS)
283:       ! Write initial energy and coordinates to file283:       ! Write initial energy and coordinates to file
284:       IF (.NOT.(NODUMP)) CALL DUMPP(X,ENERGY)284:       IF (.NOT.(NODUMP)) CALL DUMPP(X,ENERGY)
285: 285: 
286: 286: 
287:       IF (TWOENDS.AND.(.NOT.TTDONE).AND.(FORCE.NE.0.0D0)) THEN287:       IF (TWOENDS.AND.(.NOT.TTDONE).AND.(FORCE.NE.0.0D0)) THEN
288:          IF (DOINTERNALSTRANSFORM) THEN288:          IF (DOINTERNALSTRANSFORM) THEN
289:             PRINT*,'ERROR - TWOENDS not available for internal coordinates'289:             PRINT*,'ERROR - TWOENDS not available for internal coordinates'
290:             STOP290:             STOP
291:          ENDIF291:          ENDIF
292:          DUMMY1=0.0D0292:          DUMMY1=0.0D0
340: C  For CHARMM internal coordinate minimisation we project the uphill direction340: C  For CHARMM internal coordinate minimisation we project the uphill direction
341: C  out of the Cartesian gradient every time. If MYLBFGS is reset on each call then341: C  out of the Cartesian gradient every time. If MYLBFGS is reset on each call then
342: C  the total step should be a linear combination of gradients that all have this342: C  the total step should be a linear combination of gradients that all have this
343: C  component removed. However - if we don't reset then on previous calls the343: C  component removed. However - if we don't reset then on previous calls the
344: C  uphill direction will be different! Suggests that we should reset in bfgsts.344: C  uphill direction will be different! Suggests that we should reset in bfgsts.
345: C345: C
346:          DO J2=1,NUP346:          DO J2=1,NUP
347:             IF (FREEZE) THEN347:             IF (FREEZE) THEN
348:                DO J1=1,NATOMS348:                DO J1=1,NATOMS
349:                   IF (.NOT.FROZEN(J1)) CYCLE349:                   IF (.NOT.FROZEN(J1)) CYCLE
350:                   IF (VARIABLES) THEN350:                   ZWORK(3*(J1-1)+1,J2)=0.0D0
351:                      ZWORK(J1,J2)=0.0D0351:                   ZWORK(3*(J1-1)+2,J2)=0.0D0
352:                   ELSE352:                   ZWORK(3*(J1-1)+3,J2)=0.0D0
353:                      ZWORK(3*(J1-1)+1,J2)=0.0D0 
354:                      ZWORK(3*(J1-1)+2,J2)=0.0D0 
355:                      ZWORK(3*(J1-1)+3,J2)=0.0D0 
356:                   ENDIF 
357:                ENDDO353:                ENDDO
358:             ENDIF354:             ENDIF
359: ! jwrm2> Make sure the z coordinates in 2D systems remain zero355: ! jwrm2> Make sure the z coordinates in 2D systems remain zero
360:             IF (TWOD .OR. GTHOMSONT) THEN356:             IF (TWOD .OR. GTHOMSONT) THEN
361:               DO J1 = 1, NATOMS357:               DO J1 = 1, NATOMS
362:                 ZWORK(3*J1, J2) = 0.0D0358:                 ZWORK(3*J1, J2) = 0.0D0
363:               END DO359:               END DO
364:             END IF360:             END IF
365: ! jwrm2> end361: ! jwrm2> end
366:             DUMMY2=0.0D0362:             DUMMY2=0.0D0
764: !     PRINT '(A,6G20.10)','W=',W(1:3)760: !     PRINT '(A,6G20.10)','W=',W(1:3)
765: !     PRINT '(A,6G20.10)','ZWORK=',ZWORK(1:3,1)761: !     PRINT '(A,6G20.10)','ZWORK=',ZWORK(1:3,1)
766: 762: 
767: !     COMMENT vr274: I did not change this if block to the pointer stuff, not used with internals763: !     COMMENT vr274: I did not change this if block to the pointer stuff, not used with internals
768:       IF (PROJECT.AND.(.NOT.TWOENDS).AND.(.NOT.DOINTERNALSTRANSFORM).AND.(.NOT.REVERSEUPHILLT)) THEN764:       IF (PROJECT.AND.(.NOT.TWOENDS).AND.(.NOT.DOINTERNALSTRANSFORM).AND.(.NOT.REVERSEUPHILLT)) THEN
769:          DO J2=1,NUP765:          DO J2=1,NUP
770: 766: 
771:             IF (FREEZE) THEN ! this may duplicate the block in the projection of the gradient767:             IF (FREEZE) THEN ! this may duplicate the block in the projection of the gradient
772:                DO J1=1,NATOMS768:                DO J1=1,NATOMS
773:                   IF (.NOT.FROZEN(J1)) CYCLE769:                   IF (.NOT.FROZEN(J1)) CYCLE
774:                   IF (VARIABLES) THEN770:                   ZWORK(3*(J1-1)+1,J2)=0.0D0
775:                      ZWORK(J1,J2)=0.0D0771:                   ZWORK(3*(J1-1)+2,J2)=0.0D0
776:                   ELSE772:                   ZWORK(3*(J1-1)+3,J2)=0.0D0
777:                      ZWORK(3*(J1-1)+1,J2)=0.0D0 
778:                      ZWORK(3*(J1-1)+2,J2)=0.0D0 
779:                      ZWORK(3*(J1-1)+3,J2)=0.0D0 
780:                   ENDIF 
781:                ENDDO773:                ENDDO
782:             ENDIF774:             ENDIF
783: ! jwrm2> Make sure the z coordinates in 2D systems remain zero775: ! jwrm2> Make sure the z coordinates in 2D systems remain zero
784:             IF (TWOD .OR. GTHOMSONT) THEN776:             IF (TWOD .OR. GTHOMSONT) THEN
785:                DO J1 = 1, NATOMS777:                DO J1 = 1, NATOMS
786:                   ZWORK(3*J1, J2) = 0.0D0778:                   ZWORK(3*J1, J2) = 0.0D0
787:                END DO779:                END DO
788:              END IF780:              END IF
789: ! jwrm2> end781: ! jwrm2> end
790: 782: 
1045: C        CALL POTENTIAL(X,EMINUS,GDUM,.TRUE.,.FALSE.,RMSDUM,.FALSE.,.FALSE.)1037: C        CALL POTENTIAL(X,EMINUS,GDUM,.TRUE.,.FALSE.,RMSDUM,.FALSE.,.FALSE.)
1046: C        EMINUS=DDOT(3*NATOMS,GDUM,1,GDUM,1)1038: C        EMINUS=DDOT(3*NATOMS,GDUM,1,GDUM,1)
1047: C        X(J1)=DUMMY11039: C        X(J1)=DUMMY1
1048: C        IF (100.0D0*ABS((VEC2(J1)-(EPLUS-EMINUS)/(2.0D0*SHIT))/VEC2(J1)).GT.1.0D0) 1040: C        IF (100.0D0*ABS((VEC2(J1)-(EPLUS-EMINUS)/(2.0D0*SHIT))/VEC2(J1)).GT.1.0D0) 
1049: C    1     WRITE(*,'(A,I5,5G20.10)') 'J1,anal,num,%,+,-=',J1,VEC2(J1),(EPLUS-EMINUS)/(2.0D0*SHIT),1041: C    1     WRITE(*,'(A,I5,5G20.10)') 'J1,anal,num,%,+,-=',J1,VEC2(J1),(EPLUS-EMINUS)/(2.0D0*SHIT),
1050: C    2                  100.0D0*ABS((VEC2(J1)-(EPLUS-EMINUS)/(2.0D0*SHIT))/VEC2(J1)),EPLUS,EMINUS1042: C    2                  100.0D0*ABS((VEC2(J1)-(EPLUS-EMINUS)/(2.0D0*SHIT))/VEC2(J1)),EPLUS,EMINUS
1051: C     ENDDO1043: C     ENDDO
1052: C     FIXIMAGE=FIXSAVE1044: C     FIXIMAGE=FIXSAVE
1053: C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1045: C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1054: 1046: 
1055:       G(1:N)=GSAVE(1:N)1047:       G(1:3*NATOMS)=GSAVE(1:3*NATOMS)
1056: 1048: 
1057: C  vr274: I didn't fully understand this if/elseif block, not changed to pointer stuff1049: C  vr274: I didn't fully understand this if/elseif block, not changed to pointer stuff
1058: 1050: 
1059:       IF (TWOENDS.AND.(.NOT.TTDONE)) THEN1051:       IF (TWOENDS.AND.(.NOT.TTDONE)) THEN
1060:          IF (DOINTERNALSTRANSFORM) THEN1052:          IF (DOINTERNALSTRANSFORM) THEN
1061:             PRINT*,'ERROR - cant use twoends with internals'1053:             PRINT*,'ERROR - cant use twoends with internals'
1062:             STOP1054:             STOP
1063:          ENDIF1055:          ENDIF
1064:          IF (FORCE.NE.0.0D0) THEN1056:          IF (FORCE.NE.0.0D0) THEN
1065:             DUMMY1=0.0D01057:             DUMMY1=0.0D0
1238: C           OLDQ(1:N)=OLDQ(1:N)+DELTAQ(1:N)1230: C           OLDQ(1:N)=OLDQ(1:N)+DELTAQ(1:N)
1239:             OLDQ(1:N)=XINT(1:N)1231:             OLDQ(1:N)=XINT(1:N)
1240:          ELSEIF (UNRST) THEN1232:          ELSEIF (UNRST) THEN
1241: !           TEST1(1:N)=X(1:N)1233: !           TEST1(1:N)=X(1:N)
1242:             CALL geom_to_var(N,X(1:N)) ! testing!!! - to put X back into register with the common block internals (and g)1234:             CALL geom_to_var(N,X(1:N)) ! testing!!! - to put X back into register with the common block internals (and g)
1243: !           CALL geom_to_var(N,TEST1(1:N))1235: !           CALL geom_to_var(N,TEST1(1:N))
1244: !           do j1=1,N1236: !           do j1=1,N
1245: !           if (abs((TEST1(j1)-x(j1))/x(j1))*100.0d0.gt.1.0D-6) print *,'hello coords ',J11237: !           if (abs((TEST1(j1)-x(j1))/x(j1))*100.0d0.gt.1.0D-6) print *,'hello coords ',J1
1246: !           enddo1238: !           enddo
1247:          ENDIF1239:          ENDIF
1248:          GLAST(1:N)=GSAVE(1:N) 1240:          GLAST(1:3*NATOMS)=GSAVE(1:3*NATOMS) 
1249:          ! Write energy and coordinates for this step to file1241:          ! Write energy and coordinates for this step to file
1250:          IF (.NOT.(NODUMP)) THEN ! jmc dumpp dumps x but x is in internals...1242:          IF (.NOT.(NODUMP)) THEN ! jmc dumpp dumps x but x is in internals...
1251:             IF (INTMINT) THEN1243:             IF (INTMINT) THEN
1252:                IF (UNRST) THEN1244:                IF (UNRST) THEN
1253: !                 TMPINT(1:N)=X(1:N) ! the Cartesians should not have changed from when they were set before the 1245: !                 TMPINT(1:N)=X(1:N) ! the Cartesians should not have changed from when they were set before the 
1254:                                      ! call to potential above, so commented out these three lines1246:                                      ! call to potential above, so commented out these three lines
1255: !                 CALL var_to_geom(N,TMPINT)1247: !                 CALL var_to_geom(N,TMPINT)
1256: !                 CALL chainbuild1248: !                 CALL chainbuild
1257:                   DO I1=1,nres1249:                   DO I1=1,nres
1258:                      CART(6*(I1-1)+1)=c(1,I1)1250:                      CART(6*(I1-1)+1)=c(1,I1)


r32316/Natb.f 2017-04-22 17:30:11.866863634 +0100 r32315/Natb.f 2017-04-22 17:30:15.898915725 +0100
 52: C          print *,'ANION' 52: C          print *,'ANION'
 53: C       endif 53: C       endif
 54:  54: 
 55:         CALL entots(P,grad,nsize,e0) 55:         CALL entots(P,grad,nsize,e0)
 56:  56: 
 57:         end 57:         end
 58: c______________________________________________________________________ 58: c______________________________________________________________________
 59:  59: 
 60:         SUBROUTINE entots(P,DP,nbat,e0) 60:         SUBROUTINE entots(P,DP,nbat,e0)
 61:         IMPLICIT DOUBLE PRECISION(A-H,O-Z) 61:         IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 62:         IMPLICIT INTEGER(I-N) 
 63:         PARAMETER(ANM=41901.052,RKB=3.165D-6,erepmax=100.d0) 62:         PARAMETER(ANM=41901.052,RKB=3.165D-6,erepmax=100.d0)
 64:  63: 
 65:         DIMENSION P(3*NBAT) 64:         DIMENSION P(3*NBAT)
 66:         DIMENSION GRAD(3*NBAT),DP(3*NBAT) 65:         DIMENSION GRAD(3*NBAT),DP(3*NBAT)
 67:  66: 
 68:         DIMENSION HDER(3*NBAT,NBAT*(NBAT+1)/2),R(NBAT*(NBAT+1)/2) 67:         DIMENSION HDER(3*NBAT,NBAT*(NBAT+1)/2),R(NBAT*(NBAT+1)/2)
 69:         DIMENSION HA(NBAT*(NBAT+1)/2) 68:         DIMENSION HA(NBAT*(NBAT+1)/2)
 70:         DIMENSION VECT(NBAT*(NBAT+1)/2),RINV(NBAT*(NBAT+1)/2) 69:         DIMENSION VECT(NBAT*(NBAT+1)/2),RINV(NBAT*(NBAT+1)/2)
 71:         DIMENSION H(NBAT*(NBAT+1)/2),DIF(3*NBAT*(NBAT+1)/2) 70:         DIMENSION H(NBAT*(NBAT+1)/2),DIF(3*NBAT*(NBAT+1)/2)
 72:         DIMENSION DBSZ(3*NBAT*(NBAT+1)/2),UNIT(3*NBAT*(NBAT+1)/2) 71:         DIMENSION DBSZ(3*NBAT*(NBAT+1)/2),UNIT(3*NBAT*(NBAT+1)/2)
537:            enddo536:            enddo
538:            537:            
539:            e0=e0+erep538:            e0=e0+erep
540:            539:            
541:         return540:         return
542:         end541:         end
543: 542: 
544: c-----------------------------------------------------543: c-----------------------------------------------------
545:         subroutine ftss(r,f,df)544:         subroutine ftss(r,f,df)
546:         implicit DOUBLE PRECISION(a-h,o-z)545:         implicit DOUBLE PRECISION(a-h,o-z)
547:         IMPLICIT INTEGER(I-N) 
548:         dimension a(10),b(10),c(10)546:         dimension a(10),b(10),c(10)
549:         character*2 metal547:         character*2 metal
550:         common/METAL/metal548:         common/METAL/metal
551:         if (metal.eq.'NA') then549:         if (metal.eq.'NA') then
552:            a(1)=-0.0000600345550:            a(1)=-0.0000600345
553:            a(2)=0.7366633997551:            a(2)=0.7366633997
554:            a(3)=-21.2536277450552:            a(3)=-21.2536277450
555:            b(1)=7.8227026687553:            b(1)=7.8227026687
556:            b(2)=5.0784652072554:            b(2)=5.0784652072
557:            b(3)=-5.7014389921555:            b(3)=-5.7014389921
589:            xn=8.d0587:            xn=8.d0
590:            fexp=aa*dexp(-bb*r)588:            fexp=aa*dexp(-bb*r)
591:            f=(r**xn)*fexp589:            f=(r**xn)*fexp
592:            df=(xn*r**(xn-1.D0)-bb*(r**xn))*fexp590:            df=(xn*r**(xn-1.D0)-bb*(r**xn))*fexp
593:         endif591:         endif
594:         return592:         return
595:         end593:         end
596: c-----------------------------------------------------594: c-----------------------------------------------------
597:         subroutine ftsz(r,f,df)595:         subroutine ftsz(r,f,df)
598:         implicit DOUBLE PRECISION(a-h,o-z)596:         implicit DOUBLE PRECISION(a-h,o-z)
599:         IMPLICIT INTEGER(I-N) 
600:         dimension a(10),b(10),c(10)597:         dimension a(10),b(10),c(10)
601:         character*2 metal598:         character*2 metal
602:         common/METAL/metal599:         common/METAL/metal
603:         if (metal.eq.'NA') then600:         if (metal.eq.'NA') then
604:            a(1)=0.0002818237601:            a(1)=0.0002818237
605:            a(2)=-0.0025316896602:            a(2)=-0.0025316896
606:            a(3)=67.9894048289603:            a(3)=67.9894048289
607:            b(1)=9.2162096989604:            b(1)=9.2162096989
608:            b(2)= 4.0958902363605:            b(2)= 4.0958902363
609:            b(3)= -5.2076863938606:            b(3)= -5.2076863938


r32316/newconnect.f90 2017-04-22 17:30:10.430845092 +0100 r32315/newconnect.f90 2017-04-22 17:30:14.606899036 +0100
 60:           INTEGER INVERT, INDEX(NA), IMATCH, NMINSAVE, NMINSAVE2, ISTAT, MYJS, MYJF 60:           INTEGER INVERT, INDEX(NA), IMATCH, NMINSAVE, NMINSAVE2, ISTAT, MYJS, MYJF
 61:           CHARACTER(LEN=80) FNAMEF 61:           CHARACTER(LEN=80) FNAMEF
 62:           CHARACTER(LEN=20) EFNAME 62:           CHARACTER(LEN=20) EFNAME
 63:           DOUBLE PRECISION BHENERGY 63:           DOUBLE PRECISION BHENERGY
 64:           COMMON /BHINTE/ BHENERGY 64:           COMMON /BHINTE/ BHENERGY
 65:           DOUBLE PRECISION DIST2, RMAT(3,3) 65:           DOUBLE PRECISION DIST2, RMAT(3,3)
 66:  66: 
 67:           FINISHED=.FALSE. 67:           FINISHED=.FALSE.
 68:           ALLOCATE(EI,EF,Q(LNOPT),FIN(LNOPT)) 68:           ALLOCATE(EI,EF,Q(LNOPT),FIN(LNOPT))
 69:           EI=EII;EF=EFF;Q=QQ;FIN=FINFIN; 69:           EI=EII;EF=EFF;Q=QQ;FIN=FINFIN;
 70:           TSREDO(1:LNOPT)=-1.0D0 ! can be uninitialised otherwise 
 71:  70: 
 72:           MOREPRINTING=PTEST 71:           MOREPRINTING=PTEST
 73:           IF (MOREPRINTING) THEN 72:           IF (MOREPRINTING) THEN
 74:                CALL ALLKEYCONNECTPRINT 73:                CALL ALLKEYCONNECTPRINT
 75:                PRINT* 74:                PRINT*
 76:           ENDIF 75:           ENDIF
 77:           INQUIRE(FILE='redopoints',EXIST=YESNO) 76:           INQUIRE(FILE='redopoints',EXIST=YESNO)
 78:           IF (YESNO) THEN 77:           IF (YESNO) THEN
 79:              IF (REDOPATH.AND.(.NOT.REDOPATHXYZ)) THEN 78:              IF (REDOPATH.AND.(.NOT.REDOPATHXYZ)) THEN
 80:                 PRINT '(A)',' newconnect> Transition state coordinates will be read from file redopoints' 79:                 PRINT '(A)',' newconnect> Transition state coordinates will be read from file redopoints'


r32316/nnutils.f90 2017-04-22 17:30:11.422858016 +0100 r32315/nnutils.f90 2017-04-22 17:30:15.286907815 +0100
1186:                IF (.NOT.FILTH==0) THEN1186:                IF (.NOT.FILTH==0) THEN
1187:                     FILENAME=TRIM(FILENAME)//'.'//TRIM(ADJUSTL(FILTHSTR))1187:                     FILENAME=TRIM(FILENAME)//'.'//TRIM(ADJUSTL(FILTHSTR))
1188:                ENDIF1188:                ENDIF
1189:                OPEN(UNIT=UNIT,FILE=FILENAME,STATUS='replace')1189:                OPEN(UNIT=UNIT,FILE=FILENAME,STATUS='replace')
1190: !         ENDIF1190: !         ENDIF
1191: 1191: 
1192:           DUMMY=0.0D01192:           DUMMY=0.0D0
1193:           WRITE(UNIT=UNIT,FMT='(2g24.13)') dummy,eee(1)1193:           WRITE(UNIT=UNIT,FMT='(2g24.13)') dummy,eee(1)
1194:           DO I=2,NIMAGE+11194:           DO I=2,NIMAGE+1
1195:                DUMMY = DUMMY + DVEC(I-1)1195:                DUMMY = DUMMY + DVEC(I-1)
1196:                PRINT *,'I,DUMMY,EEE,DVEC=',I,DUMMY,EEE(I),DVEC(I-1) 
1197:                WRITE(UNIT=UNIT,FMT='(2g24.13)') dummy,eee(i)1196:                WRITE(UNIT=UNIT,FMT='(2g24.13)') dummy,eee(i)
1198:           ENDDO1197:           ENDDO
1199:           DUMMY = DUMMY + DVEC(NIMAGE+1)1198:           DUMMY = DUMMY + DVEC(NIMAGE+1)
1200:           WRITE(UNIT=UNIT,FMT='(2g24.13)') dummy,eee(Nimage+2)1199:           WRITE(UNIT=UNIT,FMT='(2g24.13)') dummy,eee(Nimage+2)
1201: 1200: 
1202: !         IF (.NOT.PRESENT(UNITIN)) THEN1201: !         IF (.NOT.PRESENT(UNITIN)) THEN
1203:                CLOSE(UNIT)1202:                CLOSE(UNIT)
1204: !         ENDIF1203: !         ENDIF
1205:           PRINT *, 'writeprofile> NEB profile was saved to file "'//trim(filename)//'"'1204:           PRINT *, 'writeprofile> NEB profile was saved to file "'//trim(filename)//'"'
1206:      END SUBROUTINE WRITEPROFILE1205:      END SUBROUTINE WRITEPROFILE


r32316/oldneb.f 2017-04-22 17:30:13.914890090 +0100 r32315/oldneb.f 2017-04-22 17:30:18.274946413 +0100
812:       RETURN812:       RETURN
813:       END813:       END
814: C814: C
815: C****************************************************************************************815: C****************************************************************************************
816: C816: C
817:       SUBROUTINE MAKEGRADR(X,G,ETOTAL,EREAL,RMSREAL,QSTART,DISCONT,LLGMAX)817:       SUBROUTINE MAKEGRADR(X,G,ETOTAL,EREAL,RMSREAL,QSTART,DISCONT,LLGMAX)
818:       USE COMMONS818:       USE COMMONS
819:       USE MODTWOEND819:       USE MODTWOEND
820:       USE MODNEB820:       USE MODNEB
821:       USE MODHESS821:       USE MODHESS
822:       USE KEY,ONLY : NEBK, FROZEN, FREEZE, VARIABLES822:       USE KEY,ONLY : NEBK, FROZEN, FREEZE
823:       use porfuncs823:       use porfuncs
824:       IMPLICIT NONE824:       IMPLICIT NONE
825:       INTEGER J1,J2825:       INTEGER J1,J2
826:       DOUBLE PRECISION DUMMY1,DUMMY2,EREAL(NIMAGEMAX),ETOTAL,TANVEC(3*NATOMS,NIMAGEMAX),WPLUS,WMINUS,EMAX,826:       DOUBLE PRECISION DUMMY1,DUMMY2,EREAL(NIMAGEMAX),ETOTAL,TANVEC(3*NATOMS,NIMAGEMAX),WPLUS,WMINUS,EMAX,
827:      1                 XX(3*NATOMS),GG(3*NATOMS),RMS,G(3*NATOMS*NIMAGEMAX),X(3*NATOMS*NIMAGEMAX),LLGMAX(3*NATOMS),827:      1                 XX(3*NATOMS),GG(3*NATOMS),RMS,G(3*NATOMS*NIMAGEMAX),X(3*NATOMS*NIMAGEMAX),LLGMAX(3*NATOMS),
828:      2                 QSTART(3*NATOMS),RMSREAL(NIMAGEMAX),DVEC(NIMAGEMAX+1),EINITIAL,EFINAL,SEPARATION,PI2828:      2                 QSTART(3*NATOMS),RMSREAL(NIMAGEMAX),DVEC(NIMAGEMAX+1),EINITIAL,EFINAL,SEPARATION,PI2
829:       PARAMETER (PI2=6.283185307D0)829:       PARAMETER (PI2=6.283185307D0)
830:       COMMON /NEBRMS/ RMS,EINITIAL,EFINAL,SEPARATION830:       COMMON /NEBRMS/ RMS,EINITIAL,EFINAL,SEPARATION
831:       COMMON /DVNEB/ DVEC831:       COMMON /DVNEB/ DVEC
832:       INTEGER IADD, LNOPT, LNATOMS, ISTAT832:       INTEGER IADD, LNOPT, LNATOMS, ISTAT
1052:             G(J2+(LNOPT+IADD)*(J1-1))=G(J2+(LNOPT+IADD)*(J1-1))-(DUMMY1-DUMMY2)*TANVEC(J2,J1)1052:             G(J2+(LNOPT+IADD)*(J1-1))=G(J2+(LNOPT+IADD)*(J1-1))-(DUMMY1-DUMMY2)*TANVEC(J2,J1)
1053:          ENDDO1053:          ENDDO
1054:       ENDDO1054:       ENDDO
1055: !   1055: !   
1056: ! Set gradients on frozen atoms to zero.1056: ! Set gradients on frozen atoms to zero.
1057: !   1057: !   
1058:       IF (FREEZE) THEN1058:       IF (FREEZE) THEN
1059:          DO J1=1,NIMAGE1059:          DO J1=1,NIMAGE
1060:             DO J2=1,NATOMS1060:             DO J2=1,NATOMS
1061:                IF (FROZEN(J2)) THEN1061:                IF (FROZEN(J2)) THEN
1062:                   IF (VARIABLES) THEN1062:                   G(LNOPT*(J1-1)+3*(J2-1)+1)=0.0D0
1063:                      G(LNOPT*(J1-1)+J2)=0.0D01063:                   G(LNOPT*(J1-1)+3*(J2-1)+2)=0.0D0
1064:                   ELSE1064:                   G(LNOPT*(J1-1)+3*(J2-1)+3)=0.0D0
1065:                      G(LNOPT*(J1-1)+3*(J2-1)+1)=0.0D0 
1066:                      G(LNOPT*(J1-1)+3*(J2-1)+2)=0.0D0 
1067:                      G(LNOPT*(J1-1)+3*(J2-1)+3)=0.0D0 
1068:                   ENDIF 
1069:                ENDIF1065:                ENDIF
1070:             ENDDO1066:             ENDDO
1071:          ENDDO1067:          ENDDO
1072:       ENDIF1068:       ENDIF
1073: 1069: 
1074:       RETURN1070:       RETURN
1075:       END1071:       END
1076: C1072: C
1077: C  Put a specified number of images between specified starting and1073: C  Put a specified number of images between specified starting and
1078: C  finishing points using suitable interpolation.1074: C  finishing points using suitable interpolation.


r32316/oldnebgradient.f90 2017-04-22 17:30:11.650860851 +0100 r32315/oldnebgradient.f90 2017-04-22 17:30:15.654912529 +0100
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE OLDNEBGRADIENT 19: SUBROUTINE OLDNEBGRADIENT
 20: USE KEYGRAD 20: USE KEYGRAD
 21: USE NEBDATA 21: USE NEBDATA
 22: USE KEYNEB,ONLY: NIMAGE 22: USE KEYNEB,ONLY: NIMAGE
 23: USE TANGENT 23: USE TANGENT
 24: USE NEBUTILS 24: USE NEBUTILS
 25: USE KEY, ONLY: UNRST, FROZEN, FREEZE, NEBK, BULKT, TWOD, VARIABLES 25: USE KEY, ONLY: UNRST, FROZEN, FREEZE, NEBK, BULKT, TWOD
 26: USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3 26: USE COMMONS, ONLY : PARAM1, PARAM2, PARAM3
 27: USE GRADIENTS 27: USE GRADIENTS
 28: IMPLICIT NONE 28: IMPLICIT NONE
 29: INTEGER J1, J2 29: INTEGER J1, J2
 30: DOUBLE PRECISION DUMMY1, DUMMY2 30: DOUBLE PRECISION DUMMY1, DUMMY2
 31:  31: 
 32: ! 32: !
 33: ! Image 1 is START 33: ! Image 1 is START
 34: ! Movable images are numbers 2 to NIMAGE+1 34: ! Movable images are numbers 2 to NIMAGE+1
 35: ! Image NIMAGE+2 is FINISH 35: ! Image NIMAGE+2 is FINISH
118: !  PRINT '(A,I6)','oldnebgrad> GGG g perp + with parallel spring derivatives for image ',J1118: !  PRINT '(A,I6)','oldnebgrad> GGG g perp + with parallel spring derivatives for image ',J1
119: !  PRINT '(6G20.10)',GGG(1+NOPT*(J1-1):J1*NOPT)119: !  PRINT '(6G20.10)',GGG(1+NOPT*(J1-1):J1*NOPT)
120: ENDDO120: ENDDO
121: !   121: !   
122: ! Set gradients on frozen atoms to zero.122: ! Set gradients on frozen atoms to zero.
123: !   123: !   
124: IF (FREEZE) THEN124: IF (FREEZE) THEN
125:    DO J1=1,NIMAGE+2125:    DO J1=1,NIMAGE+2
126:       DO J2=1,NATOMS126:       DO J2=1,NATOMS
127:          IF (FROZEN(J2)) THEN127:          IF (FROZEN(J2)) THEN
128:             IF (VARIABLES) THEN128:             GGG(NOPT*(J1-1)+3*(J2-1)+1)=0.0D0
129:                GGG(NOPT*(J1-1)+J2)=0.0D0129:             GGG(NOPT*(J1-1)+3*(J2-1)+2)=0.0D0
130:             ELSE130:             GGG(NOPT*(J1-1)+3*(J2-1)+3)=0.0D0
131:                GGG(NOPT*(J1-1)+3*(J2-1)+1)=0.0D0 
132:                GGG(NOPT*(J1-1)+3*(J2-1)+2)=0.0D0 
133:                GGG(NOPT*(J1-1)+3*(J2-1)+3)=0.0D0 
134:             ENDIF 
135:          ENDIF131:          ENDIF
136:       ENDDO132:       ENDDO
137:    ENDDO133:    ENDDO
138: ENDIF134: ENDIF
139: 135: 
140: IF (ORT) THEN ! REMOVE OVERALL ROTATION/TRANSLATION BY FREEZING 6 DEGREES OF FREEDOM136: IF (ORT) THEN ! REMOVE OVERALL ROTATION/TRANSLATION BY FREEZING 6 DEGREES OF FREEDOM
141:    G(1::NOPT)=0.0D0137:    G(1::NOPT)=0.0D0
142:    G(2::NOPT)=0.0D0138:    G(2::NOPT)=0.0D0
143:    G(3::NOPT)=0.0D0139:    G(3::NOPT)=0.0D0
144:    G(4::NOPT)=0.0D0140:    G(4::NOPT)=0.0D0


r32316/OPTIM.F 2017-04-22 17:30:12.086866478 +0100 r32315/OPTIM.F 2017-04-22 17:30:16.450923866 +0100
452:             WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV452:             WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV
453:          ELSE453:          ELSE
454:             WRITE(*,'(A)') ' OPTIM> z rotational ev shifting'454:             WRITE(*,'(A)') ' OPTIM> z rotational ev shifting'
455:          ENDIF455:          ENDIF
456:          SHIFTL(1)=SHIFTV456:          SHIFTL(1)=SHIFTV
457:          SHIFTL(2)=SHIFTV457:          SHIFTL(2)=SHIFTV
458:          SHIFTL(3)=SHIFTV458:          SHIFTL(3)=SHIFTV
459:          SHIFTL(4)=SHIFTV459:          SHIFTL(4)=SHIFTV
460:          SHIFTL(5)=SHIFTV460:          SHIFTL(5)=SHIFTV
461:          SHIFTL(6)=SHIFTV461:          SHIFTL(6)=SHIFTV
462:       ELSE IF (MIEFT.OR.NOTRANSROTT) THEN462:       ELSE IF (MIEFT.OR.MLP3T.OR.MLPB3T.OR.MLPVB3T.OR.NOTRANSROTT) THEN
463:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'463:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'
464:       ELSE IF (VARIABLES) THEN 
465:          WRITE(*,'(A,I6,A)') ' OPTIM> ',NZERO,' zero eigenvalues shifted' 
466:          DO J1=1,NZERO 
467:             SHIFTL(J1)=SHIFTV 
468:          ENDDO 
469:       ELSE464:       ELSE
470:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV465:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV
471:          DO J1=1,6466:          DO J1=1,6
472:             SHIFTL(J1)=SHIFTV467:             SHIFTL(J1)=SHIFTV
473:          ENDDO468:          ENDDO
474:       ENDIF469:       ENDIF
475: 470: 
476:       !js850> 471:       !js850> 
477:       IF (ZSYM(NATOMS).EQ.'SS') THEN472:       IF (ZSYM(NATOMS).EQ.'SS') THEN
478:          ALLOCATE(SOFT_SPHERE_RADII(NATOMS))473:          ALLOCATE(SOFT_SPHERE_RADII(NATOMS))
1003:       IF (ALLOCATED(ZWORK)) DEALLOCATE(ZWORK)998:       IF (ALLOCATED(ZWORK)) DEALLOCATE(ZWORK)
1004:       IF (ALLOCATED(HESS)) DEALLOCATE(HESS)999:       IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
1005:       IF (ALLOCATED(FIN)) DEALLOCATE(FIN)1000:       IF (ALLOCATED(FIN)) DEALLOCATE(FIN)
1006:       IF (ALLOCATED(START)) DEALLOCATE(START)1001:       IF (ALLOCATED(START)) DEALLOCATE(START)
1007:       IF (ALLOCATED(QW)) DEALLOCATE(QW)1002:       IF (ALLOCATED(QW)) DEALLOCATE(QW)
1008:       IF (ALLOCATED(MIEF_SIG)) DEALLOCATE(MIEF_SIG)1003:       IF (ALLOCATED(MIEF_SIG)) DEALLOCATE(MIEF_SIG)
1009:       IF (ALLOCATED(MIEF_EPS)) DEALLOCATE(MIEF_EPS)1004:       IF (ALLOCATED(MIEF_EPS)) DEALLOCATE(MIEF_EPS)
1010:       IF (ALLOCATED(MIEF_SITES)) DEALLOCATE(MIEF_SITES)1005:       IF (ALLOCATED(MIEF_SITES)) DEALLOCATE(MIEF_SITES)
1011:       IF (ALLOCATED(MIEF_U_RCUT)) DEALLOCATE(MIEF_U_RCUT)1006:       IF (ALLOCATED(MIEF_U_RCUT)) DEALLOCATE(MIEF_U_RCUT)
1012:       IF (ALLOCATED(MIEF_DUDR_RCUT)) DEALLOCATE(MIEF_DUDR_RCUT)1007:       IF (ALLOCATED(MIEF_DUDR_RCUT)) DEALLOCATE(MIEF_DUDR_RCUT)
1013:       IF (ALLOCATED(MLPDAT)) DEALLOCATE(MLPDAT) 1008:       
1014:       IF (ALLOCATED(MLPOUTCOME)) DEALLOCATE(MLPOUTCOME)  
1015:       1009:       
1016:       IF (DUMPMAG) CLOSE(321)1010:       IF (DUMPMAG) CLOSE(321)
1017:       CALL FLUSH(6) 1011:       CALL FLUSH(6) 
1018:       1012:       
1019:       RETURN1013:       RETURN
1020: 1014: 
1021:       END1015:       END
1022: 1016: 
1023:       SUBROUTINE TSUMMARY1017:       SUBROUTINE TSUMMARY
1024:       USE KEY,ONLY : TSTART1018:       USE KEY,ONLY : TSTART


r32316/potential.f 2017-04-22 17:30:14.146893088 +0100 r32315/potential.f 2017-04-22 17:30:18.494949254 +0100
2080: 2080: 
2081:             ! read gradients and freeze them2081:             ! read gradients and freeze them
2082:             LUNIT=GETUNIT()2082:             LUNIT=GETUNIT()
2083:             OPEN(LUNIT,FILE='xmolout',STATUS='OLD') !assuming that the output is called like this!2083:             OPEN(LUNIT,FILE='xmolout',STATUS='OLD') !assuming that the output is called like this!
2084: 2084: 
2085:             READ(LUNIT,'(A80)',END=812) FNAME2085:             READ(LUNIT,'(A80)',END=812) FNAME
2086:             READ(LUNIT,'(A80)',END=812) FNAME2086:             READ(LUNIT,'(A80)',END=812) FNAME
2087:             DO J1=1,NATOMS2087:             DO J1=1,NATOMS
2088:                READ(LUNIT,*) STRING,XTEMP,YTEMP,ZTEMP,2088:                READ(LUNIT,*) STRING,XTEMP,YTEMP,ZTEMP,
2089:      1         VNEW(3*(J1-1)+1),VNEW(3*(J1-1)+2),VNEW(3*(J1-1)+3),NDUMMY2089:      1         VNEW(3*(J1-1)+1),VNEW(3*(J1-1)+2),VNEW(3*(J1-1)+3),NDUMMY
2090:                IF (FROZEN(J1)) THEN2090:                IF(FROZEN(J1)) THEN
2091:                   IF (VARIABLES) THEN2091:                   VNEW(3*(J1-1)+1)=0.D0
2092:                      VNEW(J1)=0.D02092:                   VNEW(3*(J1-1)+2)=0.D0
2093:                   ELSE2093:                   VNEW(3*(J1-1)+3)=0.D0
2094:                      VNEW(3*(J1-1)+1)=0.D0 
2095:                      VNEW(3*(J1-1)+2)=0.D0 
2096:                      VNEW(3*(J1-1)+3)=0.D0 
2097:                   ENDIF 
2098:                ENDIF2094:                ENDIF
2099:             ENDDO2095:             ENDDO
2100: 2096: 
2101: 812         CONTINUE2097: 812         CONTINUE
2102:             CLOSE(LUNIT)2098:             CLOSE(LUNIT)
2103:             VNEW(1:3*NATOMS)=VNEW(1:3*NATOMS)/23.06054D0   !kcal/Angstrom to eV/Angstrom.2099:             VNEW(1:3*NATOMS)=VNEW(1:3*NATOMS)/23.06054D0   !kcal/Angstrom to eV/Angstrom.
2104:             ! WRITE(*,*) 'VNEW: ',VNEW2100:             ! WRITE(*,*) 'VNEW: ',VNEW
2105: 2101: 
2106: 2102: 
2107: 2103: 
4039:                STOP4035:                STOP
4040:             ENDIF4036:             ENDIF
4041:          ENDIF4037:          ENDIF
4042: 4038: 
4043:          IF (GTEST) THEN4039:          IF (GTEST) THEN
4044:             ! PRINT '(A,L5)',' potential> FREEZE,VNEW=',FREEZE4040:             ! PRINT '(A,L5)',' potential> FREEZE,VNEW=',FREEZE
4045:             ! PRINT '(3G20.10)',VNEW(1:3*NATOMS)4041:             ! PRINT '(3G20.10)',VNEW(1:3*NATOMS)
4046:             IF (FREEZE) THEN4042:             IF (FREEZE) THEN
4047:                DO J1=1,NATOMS4043:                DO J1=1,NATOMS
4048:                   IF (FROZEN(J1)) THEN4044:                   IF (FROZEN(J1)) THEN
4049:                      IF (VARIABLES) THEN4045:                      VNEW(3*(J1-1)+1)=0.0D0
4050:                         VNEW(J1)=0.0D04046:                      VNEW(3*(J1-1)+2)=0.0D0
4051:                      ELSE4047:                      VNEW(3*(J1-1)+3)=0.0D0
4052:                         VNEW(3*(J1-1)+1)=0.0D0 
4053:                         VNEW(3*(J1-1)+2)=0.0D0 
4054:                         VNEW(3*(J1-1)+3)=0.0D0 
4055:                      ENDIF 
4056:                   ENDIF4048:                   ENDIF
4057:                ENDDO4049:                ENDDO
4058:             ENDIF4050:             ENDIF
4059: 4051: 
4060:             ! sn402: other half of hack to allow genrigid to work with arbitrary potential4052:             ! sn402: other half of hack to allow genrigid to work with arbitrary potential
4061:             IF (RIGIDMOLECULEST .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN4053:             IF (RIGIDMOLECULEST .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.)) THEN
4062:                 IF (STEST) THEN4054:                 IF (STEST) THEN
4063:                      CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)4055:                      CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)
4064:                          HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D04056:                          HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0
4065:                          HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D04057:                          HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0
4158:                   ! We need to have correct zero eigenvalues for path.info, for example.4150:                   ! We need to have correct zero eigenvalues for path.info, for example.
4159:                   ! 4151:                   ! 
4160:                   ! HESS(J2,J2)=SHIFTV4152:                   ! HESS(J2,J2)=SHIFTV
4161:                ENDIF4153:                ENDIF
4162:             ENDDO4154:             ENDDO
4163:          ENDIF4155:          ENDIF
4164:          IF (STEST) THEN4156:          IF (STEST) THEN
4165:             IF (FREEZE) THEN4157:             IF (FREEZE) THEN
4166:                DO J1=1,NATOMS4158:                DO J1=1,NATOMS
4167:                   IF (FROZEN(J1)) THEN4159:                   IF (FROZEN(J1)) THEN
4168:                      IF (VARIABLES) THEN4160:                      DO J2=1,NOPT
4169:                         DO J2=1,NOPT4161:                         HESS(3*(J1-1)+1,J2)=0.0D0
4170:                            HESS(J1,J2)=0.0D04162:                         HESS(3*(J1-1)+2,J2)=0.0D0
4171:                            HESS(J2,J1)=0.0D04163:                         HESS(3*(J1-1)+3,J2)=0.0D0
4172:                         ENDDO4164:                         HESS(J2,3*(J1-1)+1)=0.0D0
4173:                      ELSE4165:                         HESS(J2,3*(J1-1)+2)=0.0D0
4174:                         DO J2=1,NOPT4166:                         HESS(J2,3*(J1-1)+3)=0.0D0
4175:                            HESS(3*(J1-1)+1,J2)=0.0D04167:                      ENDDO
4176:                            HESS(3*(J1-1)+2,J2)=0.0D0 
4177:                            HESS(3*(J1-1)+3,J2)=0.0D0 
4178:                            HESS(J2,3*(J1-1)+1)=0.0D0 
4179:                            HESS(J2,3*(J1-1)+2)=0.0D0 
4180:                            HESS(J2,3*(J1-1)+3)=0.0D0 
4181:                         ENDDO 
4182:                      ENDIF 
4183:                   ENDIF4168:                   ENDIF
4184:                ENDDO4169:                ENDDO
4185:             ENDIF4170:             ENDIF
4186:          ENDIF4171:          ENDIF
4187:          IF (INVERTPT) THEN4172:          IF (INVERTPT) THEN
4188:             ENERGY=-ENERGY4173:             ENERGY=-ENERGY
4189:             IF (GTEST) VNEW(1:NOPT)=-VNEW(1:NOPT)4174:             IF (GTEST) VNEW(1:NOPT)=-VNEW(1:NOPT)
4190:             IF (SSTEST) HESS(1:NOPT,1:NOPT)=-HESS(1:NOPT,1:NOPT)4175:             IF (SSTEST) HESS(1:NOPT,1:NOPT)=-HESS(1:NOPT,1:NOPT)
4191:          ENDIF4176:          ENDIF
4192: 4177: 


r32316/secdiag.f90 2017-04-22 17:30:14.366895928 +0100 r32315/secdiag.f90 2017-04-22 17:30:18.718952149 +0100
 92:         CALL ORTHOGOPT(LOCALV,COORDS,.TRUE.) 92:         CALL ORTHOGOPT(LOCALV,COORDS,.TRUE.)
 93:       ELSE 93:       ELSE
 94:          CALL VECNORM(LOCALV,NOPT) 94:          CALL VECNORM(LOCALV,NOPT)
 95:       ENDIF 95:       ENDIF
 96: !     PRINT '(A)','secdiag> vec after orthogopt' 96: !     PRINT '(A)','secdiag> vec after orthogopt'
 97: !     PRINT '(6G20.10)',LOCALV(1:NOPT) 97: !     PRINT '(6G20.10)',LOCALV(1:NOPT)
 98:  98: 
 99:       IF (FREEZE) THEN 99:       IF (FREEZE) THEN
100:          DO J1=1,NATOMS100:          DO J1=1,NATOMS
101:             IF (FROZEN(J1)) THEN101:             IF (FROZEN(J1)) THEN
102:                IF (VARIABLES) THEN102:                LOCALV(3*(J1-1)+1)=0.0D0
103:                   LOCALV(J1)=0.0D0103:                LOCALV(3*(J1-1)+2)=0.0D0
104:                ELSE104:                LOCALV(3*(J1-1)+3)=0.0D0
105:                   LOCALV(3*(J1-1)+1)=0.0D0 
106:                   LOCALV(3*(J1-1)+2)=0.0D0 
107:                   LOCALV(3*(J1-1)+3)=0.0D0 
108:                ENDIF 
109:             ENDIF105:             ENDIF
110:          ENDDO106:          ENDDO
111:       ENDIF107:       ENDIF
112: 108: 
113:       VECL=1.0D0109:       VECL=1.0D0
114:       ZETA=DIFF110:       ZETA=DIFF
115: 111: 
116:       ! compute the energy and gradient at coords + zeta*vec112:       ! compute the energy and gradient at coords + zeta*vec
117:       DO J1=1,NOPT113:       DO J1=1,NOPT
118:          DUMMY3(J1)=COORDS(J1)+ZETA*LOCALV(J1)114:          DUMMY3(J1)=COORDS(J1)+ZETA*LOCALV(J1)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0