hdiff output

r30776/bfgsts.f 2016-07-15 07:30:13.599193387 +0100 r30775/bfgsts.f 2016-07-15 07:30:17.503245921 +0100
258:             CALL PVOPT(COORDS,ENERGY,GRAD)258:             CALL PVOPT(COORDS,ENERGY,GRAD)
259:          ENDIF259:          ENDIF
260:          IF ((.NOT.KNOWG).OR.(STEST.AND.(.NOT.KNOWH))) THEN260:          IF ((.NOT.KNOWG).OR.(STEST.AND.(.NOT.KNOWH))) THEN
261: ! hk286 - for RBAA, do this in moving frame261: ! hk286 - for RBAA, do this in moving frame
262:             RBAANORMALMODET = .TRUE.262:             RBAANORMALMODET = .TRUE.
263:             CALL POTENTIAL(COORDS,ENERGY,GRAD,.TRUE.,STEST,RMS,PTEST,.FALSE.)263:             CALL POTENTIAL(COORDS,ENERGY,GRAD,.TRUE.,STEST,RMS,PTEST,.FALSE.)
264:             RBAANORMALMODET = .FALSE.264:             RBAANORMALMODET = .FALSE.
265:          ENDIF265:          ENDIF
266:       ENDIF266:       ENDIF
267:       CALL DUMPP(COORDS,ENERGY)267:       CALL DUMPP(COORDS,ENERGY)
268:       IF ((.NOT.VARIABLES).AND.NOIT.AND.STEST.AND.(.NOT.MIEFT) .AND. (.NOT. NOTRANSROTT)) THEN268:       IF ((.NOT.VARIABLES).AND.NOIT.AND.STEST.AND.(.NOT.MIEFT)) THEN
269:          IF (RINGPOLYMERT.AND.(RPSYSTEM(1:4).EQ.'AECK')) THEN269:          IF (RINGPOLYMERT.AND.(RPSYSTEM(1:4).EQ.'AECK')) THEN
270:          ELSE IF (ZSYM(NATOMS).EQ.'SY') THEN270:          ELSE IF (ZSYM(NATOMS).EQ.'SY') THEN
271:             CALL SHIFTSTOCK(COORDS,NATOMS)271:             CALL SHIFTSTOCK(COORDS,NATOMS)
272:          ELSEIF (RBAAT) THEN272:          ELSEIF (RBAAT) THEN
273:             CALL SHIFTRIGID(COORDS,NATOMS)273:             CALL SHIFTRIGID(COORDS,NATOMS)
274:          ELSEIF (ZSYM(NATOMS).EQ.'TH') THEN274:          ELSEIF (ZSYM(NATOMS).EQ.'TH') THEN
275:             CALL SHIFTHTH(COORDS,NATOMS)275:             CALL SHIFTHTH(COORDS,NATOMS)
276: ! hk286276: ! hk286
277:          ELSEIF (GTHOMSONT) THEN277:          ELSEIF (GTHOMSONT) THEN
278:             CALL SHIFTHGTH(COORDS,NATOMS)278:             CALL SHIFTHGTH(COORDS,NATOMS)
1091:       USE COMMONS1091:       USE COMMONS
1092:       USE KEY1092:       USE KEY
1093:       USE VECCK1093:       USE VECCK
1094:       USE GENRIGID ! hk2861094:       USE GENRIGID ! hk286
1095:       IMPLICIT NONE1095:       IMPLICIT NONE
1096:       INTEGER J2, J3, NCHECK1096:       INTEGER J2, J3, NCHECK
1097:       DOUBLE PRECISION COORDS(*), VEC1(*), DUMMY1, DUMMY2, ROOTN, VDOT,1097:       DOUBLE PRECISION COORDS(*), VEC1(*), DUMMY1, DUMMY2, ROOTN, VDOT,
1098:      1                 CMX, CMY, CMZ, AMASS(NATOMS), TMASS, RMASS(NATOMS)1098:      1                 CMX, CMY, CMZ, AMASS(NATOMS), TMASS, RMASS(NATOMS)
1099:       LOGICAL OTEST1099:       LOGICAL OTEST
1100: 1100: 
1101:       IF (MIEFT .OR. NOTRANSROTT) THEN1101:       IF (MIEFT) THEN
1102:          IF (OTEST) CALL VECNORM(VEC1,NOPT)1102:          IF (OTEST) CALL VECNORM(VEC1,NOPT)
1103:          RETURN1103:          RETURN
1104:       ENDIF1104:       ENDIF
1105:       IF (VARIABLES.OR.MKTRAPT) THEN1105:       IF (VARIABLES.OR.MKTRAPT) THEN
1106:          IF (OTEST) CALL VECNORM(VEC1,NOPT)1106:          IF (OTEST) CALL VECNORM(VEC1,NOPT)
1107:          RETURN1107:          RETURN
1108:       ENDIF1108:       ENDIF
1109:       IF (RINGPOLYMERT.AND.(RPSYSTEM(1:4).EQ.'AECK')) THEN1109:       IF (RINGPOLYMERT.AND.(RPSYSTEM(1:4).EQ.'AECK')) THEN
1110:          IF (OTEST) CALL VECNORM(VEC1,NOPT)1110:          IF (OTEST) CALL VECNORM(VEC1,NOPT)
1111:          RETURN1111:          RETURN
1382: ! jwrm2> Make sure the z coordinate in 2D systems is fixed1382: ! jwrm2> Make sure the z coordinate in 2D systems is fixed
1383:       IF (TWOD .OR. GTHOMSONT) THEN1383:       IF (TWOD .OR. GTHOMSONT) THEN
1384:          DO J1 = 1,NATOMS1384:          DO J1 = 1,NATOMS
1385:             VECS(3*J1) = 0.0D01385:             VECS(3*J1) = 0.0D0
1386:          ENDDO1386:          ENDDO
1387:       ENDIF1387:       ENDIF
1388: ! jwrm2> end1388: ! jwrm2> end
1389: 1389: 
1390:       DUMMY1=-1.0D101390:       DUMMY1=-1.0D10
1391:       DO J1=1,NEVL1391:       DO J1=1,NEVL
1392:          IF (.NOT.FREEZE .AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN1392:          IF (.NOT.FREEZE .AND..NOT.MIEFT) THEN
1393:             CALL ORTHOGOPT(VEC1,COORDS,.TRUE.)1393:             CALL ORTHOGOPT(VEC1,COORDS,.TRUE.)
1394:          ELSE1394:          ELSE
1395:             CALL VECNORM(VEC1,NOPT)1395:             CALL VECNORM(VEC1,NOPT)
1396:          ENDIF1396:          ENDIF
1397:          CALL DSYMV('U',NOPT,1.0D0,HESS,SIZE(HESS,1),VEC1,1,0.0D0,VEC2,1)1397:          CALL DSYMV('U',NOPT,1.0D0,HESS,SIZE(HESS,1),VEC1,1,0.0D0,VEC2,1)
1398: C        CALL MATMUL(VEC2,HESS,VEC1,3*NATOMS,3*NATOMS,1,NOPT,NOPT,1)1398: C        CALL MATMUL(VEC2,HESS,VEC1,3*NATOMS,3*NATOMS,1,NOPT,NOPT,1)
1399:          DUMMY2=0.0D01399:          DUMMY2=0.0D0
1400:          DO J2=1,NOPT1400:          DO J2=1,NOPT
1401:             DUMMY2=DUMMY2+VEC2(J2)**21401:             DUMMY2=DUMMY2+VEC2(J2)**2
1402:          ENDDO1402:          ENDDO
1473: ! jwrm2> Make sure the z coordinate in 2D systems is fixed1473: ! jwrm2> Make sure the z coordinate in 2D systems is fixed
1474:       IF (TWOD .OR. GTHOMSONT) THEN1474:       IF (TWOD .OR. GTHOMSONT) THEN
1475:          DO J1 = 1,NATOMS1475:          DO J1 = 1,NATOMS
1476:             VECS(3*J1) = 0.0D01476:             VECS(3*J1) = 0.0D0
1477:          ENDDO1477:          ENDDO
1478:       ENDIF1478:       ENDIF
1479: ! jwrm2> end1479: ! jwrm2> end
1480: 1480: 
1481:       DUMMY1=-1.0D101481:       DUMMY1=-1.0D10
1482:       DO J1=1,NEVS1482:       DO J1=1,NEVS
1483:          IF (.NOT.FREEZE.AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN1483:          IF (.NOT.FREEZE.AND..NOT.MIEFT) THEN
1484:             CALL ORTHOGOPT(VEC1,COORDS,.TRUE.)1484:             CALL ORTHOGOPT(VEC1,COORDS,.TRUE.)
1485:          ELSE1485:          ELSE
1486:             CALL VECNORM(VEC1,NOPT)1486:             CALL VECNORM(VEC1,NOPT)
1487:          ENDIF1487:          ENDIF
1488:          CALL DSYMV('U',NOPT,1.0D0,HESS,SIZE(HESS,1),VEC1,1,0.0D0,VEC2,1)1488:          CALL DSYMV('U',NOPT,1.0D0,HESS,SIZE(HESS,1),VEC1,1,0.0D0,VEC2,1)
1489: 1489: 
1490:          DUMMY2=0.0D01490:          DUMMY2=0.0D0
1491:          DO J2=1,NOPT1491:          DO J2=1,NOPT
1492:             DUMMY2=DUMMY2+VEC2(J2)**21492:             DUMMY2=DUMMY2+VEC2(J2)**2
1493:          ENDDO1493:          ENDDO
1591: ! jwrm2> Make sure the z coordinate in 2D systems is fixed1591: ! jwrm2> Make sure the z coordinate in 2D systems is fixed
1592:       IF (TWOD .OR. GTHOMSONT) THEN1592:       IF (TWOD .OR. GTHOMSONT) THEN
1593:          DO J1 = 1,NATOMS1593:          DO J1 = 1,NATOMS
1594:             VECS(3*J1) = 0.0D01594:             VECS(3*J1) = 0.0D0
1595:          ENDDO1595:          ENDDO
1596:       ENDIF1596:       ENDIF
1597: ! jwrm2> end1597: ! jwrm2> end
1598: 1598: 
1599:       DUMMY4=-1.0D101599:       DUMMY4=-1.0D10
1600:       DO J1=1,NEVS1600:       DO J1=1,NEVS
1601:          IF (.NOT.FREEZE .AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN1601:          IF (.NOT.FREEZE .AND..NOT.MIEFT) THEN
1602:             CALL ORTHOGOPT(VEC1,COORDS,.TRUE.)1602:             CALL ORTHOGOPT(VEC1,COORDS,.TRUE.)
1603:          ELSE1603:          ELSE
1604:             CALL VECNORM(VEC1,NOPT)1604:             CALL VECNORM(VEC1,NOPT)
1605:          ENDIF1605:          ENDIF
1606: 1606: 
1607:          CALL DSYMV('U',NOPT,1.0D0,HESS,SIZE(HESS,1),VEC1,1,0.0D0,VEC2,1)1607:          CALL DSYMV('U',NOPT,1.0D0,HESS,SIZE(HESS,1),VEC1,1,0.0D0,VEC2,1)
1608: C        CALL MATMUL(VEC2,HESS,VEC1,3*NATOMS,3*NATOMS,1,NOPT,NOPT,1)1608: C        CALL MATMUL(VEC2,HESS,VEC1,3*NATOMS,3*NATOMS,1,NOPT,NOPT,1)
1609: 1609: 
1610:          DUMMY2=0.0D01610:          DUMMY2=0.0D0
1611:          DUMMY1=0.0D01611:          DUMMY1=0.0D0
1878:       FIXIMAGE=.TRUE.1878:       FIXIMAGE=.TRUE.
1879:       LWORK=33*3*NATOMS1879:       LWORK=33*3*NATOMS
1880:       ILWORK=33*3*NATOMS1880:       ILWORK=33*3*NATOMS
1881:       IF(ALLOCATED(VECCHK)) DEALLOCATE(VECCHK) 1881:       IF(ALLOCATED(VECCHK)) DEALLOCATE(VECCHK) 
1882:       IF (.NOT.ALLOCATED(VECCHK)) ALLOCATE(VECCHK(3*NATOMS,30))1882:       IF (.NOT.ALLOCATED(VECCHK)) ALLOCATE(VECCHK(3*NATOMS,30))
1883:       IF(ALLOCATED(ZWORK)) DEALLOCATE(ZWORK) 1883:       IF(ALLOCATED(ZWORK)) DEALLOCATE(ZWORK) 
1884:       IF (.NOT.ALLOCATED(ZWORK)) ALLOCATE(ZWORK(3*NATOMS,30))1884:       IF (.NOT.ALLOCATED(ZWORK)) ALLOCATE(ZWORK(3*NATOMS,30))
1885: 1885: 
1886:       IF (NOIT) THEN1886:       IF (NOIT) THEN
1887:          CALL POTENTIAL(COORDS,ENERGY,VEC1,.TRUE.,.TRUE.,DUMMY1,PTEST,.FALSE.) 1887:          CALL POTENTIAL(COORDS,ENERGY,VEC1,.TRUE.,.TRUE.,DUMMY1,PTEST,.FALSE.) 
1888:          IF ((.NOT.VARIABLES).AND.(.NOT.MIEFT) .AND. (.NOT. NOTRANSROTT)) THEN1888:          IF ((.NOT.VARIABLES).AND.(.NOT.MIEFT)) THEN
1889:             IF (ZSYM(NATOMS).EQ.'TH') THEN1889:             IF (ZSYM(NATOMS).EQ.'TH') THEN
1890:                CALL SHIFTHTH(COORDS,NATOMS)1890:                CALL SHIFTHTH(COORDS,NATOMS)
1891: ! hk2861891: ! hk286
1892:             ELSEIF (GTHOMSONT) THEN1892:             ELSEIF (GTHOMSONT) THEN
1893:                CALL SHIFTHGTH(COORDS,NATOMS)1893:                CALL SHIFTHGTH(COORDS,NATOMS)
1894:             ELSEIF (RBAAT) THEN1894:             ELSEIF (RBAAT) THEN
1895:                CALL SHIFTRIGID(COORDS,NATOMS)1895:                CALL SHIFTRIGID(COORDS,NATOMS)
1896:             ELSE1896:             ELSE
1897:                CALL SHIFTH(COORDS,.TRUE.,NOPT,NATOMS,ATMASS)1897:                CALL SHIFTH(COORDS,.TRUE.,NOPT,NATOMS,ATMASS)
1898:             ENDIF1898:             ENDIF


r30776/connect.f 2016-07-15 07:30:13.855196833 +0100 r30775/connect.f 2016-07-15 07:30:17.759249369 +0100
1537: C1537: C
1538: C  How many zero eigenvalues?1538: C  How many zero eigenvalues?
1539: C1539: C
1540:       IF (RTEST) THEN1540:       IF (RTEST) THEN
1541:          NZERO=21541:          NZERO=2
1542:          IF (JZ.NE.0.0D0) THEN1542:          IF (JZ.NE.0.0D0) THEN
1543:             NZERO=41543:             NZERO=4
1544:          ENDIF1544:          ENDIF
1545:       ELSE IF (MIEFT) THEN1545:       ELSE IF (MIEFT) THEN
1546:          NZERO=01546:          NZERO=0
1547:       ELSE IF (NOTRANSROTT) THEN 
1548:          NZERO=0 
1549:       ELSE IF (BULKT) THEN1547:       ELSE IF (BULKT) THEN
1550:          NZERO=31548:          NZERO=3
1551:          IF (TWOD) THEN1549:          IF (TWOD) THEN
1552:             NZERO=21550:             NZERO=2
1553:          ENDIF1551:          ENDIF
1554: C     ELSE IF ((FPGRP.EQ.'DXh'.OR.FPGRP.EQ.'CXv').AND.(ZSYMSAVE(1:1).NE.'W')) THEN1552: C     ELSE IF ((FPGRP.EQ.'DXh'.OR.FPGRP.EQ.'CXv').AND.(ZSYMSAVE(1:1).NE.'W')) THEN
1555: C        NZERO=51553: C        NZERO=5
1556:       ELSE IF (VARIABLES) THEN1554:       ELSE IF (VARIABLES) THEN
1557:       ELSE IF (FIELDT) THEN1555:       ELSE IF (FIELDT) THEN
1558:          NZERO=31556:          NZERO=3
3015:       INTEGER INDEX(NATOMS), INVERT, J5, J4, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE23013:       INTEGER INDEX(NATOMS), INVERT, J5, J4, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2
3016:       DOUBLE PRECISION TPOINTS1(3*NATOMS), TARGET(3*NATOMS), 3014:       DOUBLE PRECISION TPOINTS1(3*NATOMS), TARGET(3*NATOMS), 
3017:      1                 TPOINTS2(3*NATOMS), DIST, DMIN 3015:      1                 TPOINTS2(3*NATOMS), DIST, DMIN 
3018:       LOGICAL SUCCESS3016:       LOGICAL SUCCESS
3019: 3017: 
3020:       SUCCESS=.FALSE.3018:       SUCCESS=.FALSE.
3021:       INVERT=13019:       INVERT=1
3022:       IF (TWOD) THEN3020:       IF (TWOD) THEN
3023:          CALL STD_ORIENT2D(TPOINTS2,TPOINTS1,NORBIT1,1,NORBIT2,1)3021:          CALL STD_ORIENT2D(TPOINTS2,TPOINTS1,NORBIT1,1,NORBIT2,1)
3024:       ELSE IF (BULKT) THEN3022:       ELSE IF (BULKT) THEN
3025:          IF (.NOT.FREEZE .AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN3023:          IF (.NOT.FREEZE .AND..NOT.MIEFT) THEN
3026:             CALL CENTRE(TARGET)3024:             CALL CENTRE(TARGET)
3027:             CALL CENTRE(TPOINTS2)3025:             CALL CENTRE(TPOINTS2)
3028:          ENDIF3026:          ENDIF
3029:       ELSE3027:       ELSE
3030:          CALL STD_ORIENT(TPOINTS2,TPOINTS1,NORBIT1,1,NORBIT2,1)3028:          CALL STD_ORIENT(TPOINTS2,TPOINTS1,NORBIT1,1,NORBIT2,1)
3031:       ENDIF3029:       ENDIF
3032:       NCHOOSE1=13030:       NCHOOSE1=1
3033:       NCHOOSE2=13031:       NCHOOSE2=1
3034: 20    CONTINUE3032: 20    CONTINUE
3035:       IF (TWOD) THEN3033:       IF (TWOD) THEN


r30776/efol.f90 2016-07-15 07:30:14.103200177 +0100 r30775/efol.f90 2016-07-15 07:30:18.007252703 +0100
177: !        CALL  DIIS(NOPT,SVEC,ITER-1,QSAVE,VNEW,QTS,RMS,ENERGY,DONE)177: !        CALL  DIIS(NOPT,SVEC,ITER-1,QSAVE,VNEW,QTS,RMS,ENERGY,DONE)
178: !        IF (PTEST) WRITE(*,'(A,F20.10)') ' Energy before DIIS=',ENERGY178: !        IF (PTEST) WRITE(*,'(A,F20.10)') ' Energy before DIIS=',ENERGY
179: !        CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,PTEST,.FALSE.)179: !        CALL POTENTIAL(Q,ENERGY,VNEW,.TRUE.,.TRUE.,RMS,PTEST,.FALSE.)
180: !        IF (PTEST) WRITE(*,'(A,F20.10)') ' Energy after DIIS=',ENERGY180: !        IF (PTEST) WRITE(*,'(A,F20.10)') ' Energy after DIIS=',ENERGY
181: !     ENDIF181: !     ENDIF
182: !182: !
183: !  Transformation to mass weighted coordinates if required.183: !  Transformation to mass weighted coordinates if required.
184: !184: !
185:       IF (MASST) CALL MASSWT(NATOMS,ATMASS,QTS,VNEW,.TRUE.)185:       IF (MASST) CALL MASSWT(NATOMS,ATMASS,QTS,VNEW,.TRUE.)
186: 186: 
187:       IF ((.NOT.VARIABLES).AND.(.NOT.MIEFT).AND.(.NOT.NOTRANSROTT).AND.(.NOT.PHI4MODT).AND.(.NOT.(RINGPOLYMERT.AND.(RPSYSTEM(1:4).EQ.'AECK')))) THEN187:       IF ((.NOT.VARIABLES).AND.(.NOT.MIEFT).AND.(.NOT.PHI4MODT).AND.(.NOT.(RINGPOLYMERT.AND.(RPSYSTEM(1:4).EQ.'AECK')))) THEN
188:          IF (ZSYM(NATOMS).EQ.'SY') THEN188:          IF (ZSYM(NATOMS).EQ.'SY') THEN
189:             CALL SHIFTSTOCK(QTS,NATOMS)189:             CALL SHIFTSTOCK(QTS,NATOMS)
190:          ELSEIF (RBAAT) THEN190:          ELSEIF (RBAAT) THEN
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
302:          NZERO=(NATOMS/2)+6302:          NZERO=(NATOMS/2)+6
303:          DO J1=1,NZERO303:          DO J1=1,NZERO
304:             ZT(J1)=.FALSE.304:             ZT(J1)=.FALSE.
305:          ENDDO305:          ENDDO
306:       ELSE IF (PULLT.OR.EFIELDT) THEN306:       ELSE IF (PULLT.OR.EFIELDT) THEN
307:          NZERO=4307:          NZERO=4
308:          ZT(1)=.FALSE.308:          ZT(1)=.FALSE.
309:          ZT(2)=.FALSE.309:          ZT(2)=.FALSE.
310:          ZT(3)=.FALSE.310:          ZT(3)=.FALSE.
311:          ZT(4)=.FALSE.311:          ZT(4)=.FALSE.
312:       ELSE IF (MKTRAPT.OR.MIEFT.OR.NOTRANSROTT) THEN312:       ELSE IF (MKTRAPT.OR.MIEFT) THEN
313:          NZERO=0313:          NZERO=0
314:       ELSE IF ((EYTRAPT.OR.(ZSYM(1).EQ.'BE')).AND.(.NOT.TWOD)) THEN314:       ELSE IF ((EYTRAPT.OR.(ZSYM(1).EQ.'BE')).AND.(.NOT.TWOD)) THEN
315:          NZERO=3315:          NZERO=3
316:          ZT(1)=.FALSE.316:          ZT(1)=.FALSE.
317:          ZT(2)=.FALSE.317:          ZT(2)=.FALSE.
318:          ZT(3)=.FALSE.318:          ZT(3)=.FALSE.
319:       ELSE IF (ZSYM(1).EQ.'CK') THEN319:       ELSE IF (ZSYM(1).EQ.'CK') THEN
320:          DO J1=1,NATOMS+1320:          DO J1=1,NATOMS+1
321:             ZT(J1)=.FALSE.321:             ZT(J1)=.FALSE.
322:          ENDDO322:          ENDDO


r30776/geopt.f 2016-07-15 07:30:14.351203507 +0100 r30775/geopt.f 2016-07-15 07:30:18.267256202 +0100
397:                IF (BFGSTST.OR.(INR.EQ.2)) NEXMODES=NEXMODES+1397:                IF (BFGSTST.OR.(INR.EQ.2)) NEXMODES=NEXMODES+1
398:                IF (BFGSMINT) NEXMODES=NEXMODES+2398:                IF (BFGSMINT) NEXMODES=NEXMODES+2
399:             ENDIF399:             ENDIF
400:             IF (STEALTHYT) THEN400:             IF (STEALTHYT) THEN
401:                IF (ENERGY.LT.1.0D-8) THEN401:                IF (ENERGY.LT.1.0D-8) THEN
402:                   NEXMODES=NATOMS*3 - STM*2402:                   NEXMODES=NATOMS*3 - STM*2
403:                ELSE403:                ELSE
404:                   NEXMODES=3404:                   NEXMODES=3
405:                ENDIF405:                ENDIF
406:             ENDIF406:             ENDIF
407:             IF (VARIABLES.OR.MIEFT .OR. NOTRANSROTT) NEXMODES=0407:             IF (VARIABLES.OR.MIEFT) NEXMODES=0
408:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES408:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES
409: 409: 
410:             IF (LOWESTFRQT) THEN410:             IF (LOWESTFRQT) THEN
411: C411: C
412: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file412: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file
413: C  No mass-weighting here!413: C  No mass-weighting here!
414: C  'U' specifies that the upper triangle contains the Hessian.414: C  'U' specifies that the upper triangle contains the Hessian.
415: C  Need to save HESS before this call and restore afterwards.415: C  Need to save HESS before this call and restore afterwards.
416: C416: C
417:     417:     
1687:       RETURN1687:       RETURN
1688:       END1688:       END
1689: 1689: 
1690:       SUBROUTINE MAKENUMHESS(X,NATOMS)1690:       SUBROUTINE MAKENUMHESS(X,NATOMS)
1691: C1691: C
1692: C dae1692: C dae
1693: C1693: C
1694:       USE MODHESS1694:       USE MODHESS
1695:       USE MODCHARMM1695:       USE MODCHARMM
1696:       USE KEY,ONLY : FROZEN, CASTEP, AMHT, ONETEP, ENDNUMHESS2, QCHEM, VASP1696:       USE KEY,ONLY : FROZEN, CASTEP, AMHT, ONETEP, ENDNUMHESS2, QCHEM, VASP
1697:       USE SBMDATA,ONLY : SBMHESSATOM,SBMFIRSTDD 
1698:       USE COMMONS,ONLY : DEBUG, NOPT1697:       USE COMMONS,ONLY : DEBUG, NOPT
1699:       use porfuncs1698:       use porfuncs
1700:       IMPLICIT NONE1699:       IMPLICIT NONE
1701:       LOGICAL KNOWE, KNOWG, KNOWH1700:       LOGICAL KNOWE, KNOWG, KNOWH
1702:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH1701:       COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
1703: 1702: 
1704:       INTEGER I1,J1,NATOMS,ISTAT1703:       INTEGER I1,J1,NATOMS,ISTAT
1705:       DOUBLE PRECISION X(NOPT)1704:       DOUBLE PRECISION X(NOPT)
1706:       DOUBLE PRECISION DUM(NOPT),GRAD1(NOPT),GRAD2(NOPT),DELTA,RMS,ENERGY1705:       DOUBLE PRECISION DUM(NOPT),GRAD1(NOPT),GRAD2(NOPT),DELTA,RMS,ENERGY
1707: 1706: 
1725:       IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-11724:       IF(CHRMMT.AND.ACESOLV) NCHENCALLS=ACEUPSTEP-1
1726:  1725:  
1727: C1726: C
1728: C DAE having two potential calls for each hessian evaluation increases the accuracy1727: C DAE having two potential calls for each hessian evaluation increases the accuracy
1729: C compared to the analytical solution significantly (around 4sf), relative to the alternative1728: C compared to the analytical solution significantly (around 4sf), relative to the alternative
1730: C of evaluating the gradient once at the beginning then once for each element. But obviously it also1729: C of evaluating the gradient once at the beginning then once for each element. But obviously it also
1731: C doubles the number of potential evaluations ((3N)^2/2 for N atoms), so may be not worth it for1730: C doubles the number of potential evaluations ((3N)^2/2 for N atoms), so may be not worth it for
1732: C a long run on a big system.1731: C a long run on a big system.
1733: C1732: C
1734:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))1733:       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NOPT,NOPT))
1735: !      SBMFIRSTDD = .TRUE. 
1736:       DO I1=1,NOPT1734:       DO I1=1,NOPT
1737: !        write(*,*) 'working on atom', I1 
1738:          CALL FLUSH(6)1735:          CALL FLUSH(6)
1739:          IF (FROZEN((I1-1)/3+1)) THEN1736:          IF (FROZEN((I1-1)/3+1)) THEN
1740:             DO J1=I1,NOPT1737:             DO J1=I1,NOPT
1741:                HESS(I1,J1)=0.0D01738:                HESS(I1,J1)=0.0D0
1742:                HESS(J1,I1)=0.0D01739:                HESS(J1,I1)=0.0D0
1743:             ENDDO1740:             ENDDO
1744:          ELSE1741:          ELSE
1745:             DUM(I1)=X(I1)-DELTA1742:             DUM(I1)=X(I1)-DELTA
1746: !            SBMHESSATOM=int(I1/3) 
1747:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1743:             CALL POTENTIAL(DUM,ENERGY,GRAD1,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1748:             DUM(I1)=X(I1)+DELTA1744:             DUM(I1)=X(I1)+DELTA
1749:             CALL POTENTIAL(DUM,ENERGY,GRAD2,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)1745:             CALL POTENTIAL(DUM,ENERGY,GRAD2,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
1750:             DUM(I1)=X(I1)1746:             DUM(I1)=X(I1)
1751:             DO J1=I1,NOPT1747:             DO J1=I1,NOPT
1752:                HESS(I1,J1)=(GRAD2(J1)-GRAD1(J1))/(2.0D0*DELTA)1748:                HESS(I1,J1)=(GRAD2(J1)-GRAD1(J1))/(2.0D0*DELTA)
1753:                HESS(J1,I1)=HESS(I1,J1)1749:                HESS(J1,I1)=HESS(I1,J1)
1754:             ENDDO1750:             ENDDO
1755:          ENDIF1751:          ENDIF
1756:       ENDDO1752:       ENDDO
1757: !      SBMHESSATOM=-11753: 
1758:       IF (DEBUG) WRITE(*,'(A)') ' makenumhess> Hessian made'1754:       IF (DEBUG) WRITE(*,'(A)') ' makenumhess> Hessian made'
1759:       KNOWH=.TRUE.1755:       KNOWH=.TRUE.
1760:  1756:  
1761:       RETURN1757:       RETURN
1762:       END1758:       END
1763: 1759: 
1764:       SUBROUTINE MAKENUMHESSRP(X,NOPT)1760:       SUBROUTINE MAKENUMHESSRP(X,NOPT)
1765: C1761: C
1766: C dae1762: C dae
1767: C1763: C


r30776/inertia.f 2016-07-15 07:30:14.603206898 +0100 r30775/inertia.f 2016-07-15 07:30:18.519259593 +0100
117:                ENDIF117:                ENDIF
118:             ENDDO118:             ENDDO
119:          ENDDO119:          ENDDO
120:       ENDIF120:       ENDIF
121:       ! }}}121:       ! }}}
122:       RETURN122:       RETURN
123:       END123:       END
124: 124: 
125:       SUBROUTINE INERTIA2(DUMQ,ITX,ITY,ITZ)125:       SUBROUTINE INERTIA2(DUMQ,ITX,ITY,ITZ)
126:       USE COMMONS126:       USE COMMONS
127:       USE KEY,ONLY : FREEZE, RBAAT, MIEFT, NOTRANSROTT127:       USE KEY,ONLY : FREEZE, RBAAT, MIEFT
128:       IMPLICIT NONE128:       IMPLICIT NONE
129: ! subroutine parameters 129: ! subroutine parameters 
130:       DOUBLE PRECISION ITX,ITY,ITZ,DUMQ(3*NATOMS)130:       DOUBLE PRECISION ITX,ITY,ITZ,DUMQ(3*NATOMS)
131: ! local parameters 131: ! local parameters 
132:       INTEGER J1, J2, J3, J4132:       INTEGER J1, J2, J3, J4
133:       DOUBLE PRECISION IT(3,3), CMX, CMY, CMZ, VEC(3,3), MASST133:       DOUBLE PRECISION IT(3,3), CMX, CMY, CMZ, VEC(3,3), MASST
134: C Subroutine body {{{134: C Subroutine body {{{
135: 135: 
136:       IF (RBAAT) THEN136:       IF (RBAAT) THEN
137:          CALL RBINERTIA (DUMQ, ITX, ITY, ITZ)137:          CALL RBINERTIA (DUMQ, ITX, ITY, ITZ)
151: !        ELSE151: !        ELSE
152: !           CMX=CMX+DUMQ(3*(J1-1)+1)*ATMASS(J1)152: !           CMX=CMX+DUMQ(3*(J1-1)+1)*ATMASS(J1)
153: !           CMY=CMY+DUMQ(3*(J1-1)+2)*ATMASS(J1)153: !           CMY=CMY+DUMQ(3*(J1-1)+2)*ATMASS(J1)
154: !           CMZ=CMZ+DUMQ(3*(J1-1)+3)*ATMASS(J1)154: !           CMZ=CMZ+DUMQ(3*(J1-1)+3)*ATMASS(J1)
155: !           MASST=MASST+ATMASS(J1)155: !           MASST=MASST+ATMASS(J1)
156: !        ENDIF156: !        ENDIF
157:       ENDDO157:       ENDDO
158:       CMX=CMX/MASST158:       CMX=CMX/MASST
159:       CMY=CMY/MASST159:       CMY=CMY/MASST
160:       CMZ=CMZ/MASST160:       CMZ=CMZ/MASST
161:       IF (FREEZE .OR. MIEFT .OR. NOTRANSROTT ) THEN161:       IF (FREEZE .OR. MIEFT) THEN
162:          CMX=0.0D0162:          CMX=0.0D0
163:          CMY=0.0D0163:          CMY=0.0D0
164:          CMZ=0.0D0164:          CMZ=0.0D0
165:       ENDIF 165:       ENDIF 
166:       DO J1=1,NATOMS166:       DO J1=1,NATOMS
167:          DUMQ(3*(J1-1)+1)=DUMQ(3*(J1-1)+1)-CMX167:          DUMQ(3*(J1-1)+1)=DUMQ(3*(J1-1)+1)-CMX
168:          DUMQ(3*(J1-1)+2)=DUMQ(3*(J1-1)+2)-CMY168:          DUMQ(3*(J1-1)+2)=DUMQ(3*(J1-1)+2)-CMY
169:          DUMQ(3*(J1-1)+3)=DUMQ(3*(J1-1)+3)-CMZ169:          DUMQ(3*(J1-1)+3)=DUMQ(3*(J1-1)+3)-CMZ
170:       ENDDO170:       ENDDO
171: 171: 


r30776/key.f90 2016-07-15 07:30:14.851210235 +0100 r30775/key.f90 2016-07-15 07:30:18.767262928 +0100
 47:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 47:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 48:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 48:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 49:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 49:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 50:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 50:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 51:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 51:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 52:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 52:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 53:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 53:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 54:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 54:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 55:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 55:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 56:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, & 56:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, &
 57:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT 57:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT
 58:  58: 
 59: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 59: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 60:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 60:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 61:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 61:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 62:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 62:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 63:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 63:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 64:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads 64:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads
 65:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 65:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 66:  66: 
 67: ! hk286 > generalised THOMSON problem 67: ! hk286 > generalised THOMSON problem


r30776/keywords.f 2016-07-15 07:30:15.103213628 +0100 r30775/keywords.f 2016-07-15 07:30:19.019266321 +0100
451:          OSASAT=.FALSE.451:          OSASAT=.FALSE.
452:          RPRO=1.4D0452:          RPRO=1.4D0
453:          ODIHET=.FALSE.453:          ODIHET=.FALSE.
454: 454: 
455:          ! unres stuff455:          ! unres stuff
456:          UNRST=.FALSE.456:          UNRST=.FALSE.
457:          CONSECT=.FALSE.457:          CONSECT=.FALSE.
458:          STARTRES=0458:          STARTRES=0
459:          ENDRES=0459:          ENDRES=0
460:          NUMSEC=0460:          NUMSEC=0
461:  
462:          ! remove rotations and translations 
463:          NOTRANSROTT=.FALSE. 
464:  
465:          ! 461:          ! 
466:          ! AMH  stuff462:          ! AMH  stuff
467:          AMHT=.FALSE.463:          AMHT=.FALSE.
468: 464: 
469:          FREEZE=.FALSE.465:          FREEZE=.FALSE.
470:          FREEZERANGE=.FALSE.466:          FREEZERANGE=.FALSE.
471:          FREEZEGROUPT=.FALSE.467:          FREEZEGROUPT=.FALSE.
472:          FREEZEGROUPTYPE='GT'468:          FREEZEGROUPTYPE='GT'
473:          FREEZERES=.FALSE.469:          FREEZERES=.FALSE.
474:          NFREEZE=0470:          NFREEZE=0
1050:             IF (NITEMS.GT.3) THEN1046:             IF (NITEMS.GT.3) THEN
1051:                CALL READF(MALPHA2)1047:                CALL READF(MALPHA2)
1052:             ENDIF1048:             ENDIF
1053: ! 1049: ! 
1054: ! SAT: ALLPOINTS turns on printing of coordinates to file points for intermediate steps.1050: ! SAT: ALLPOINTS turns on printing of coordinates to file points for intermediate steps.
1055: ! This is the default.1051: ! This is the default.
1056: ! 1052: ! 
1057:          ELSE IF (WORD.EQ.'ALLPOINTS') THEN1053:          ELSE IF (WORD.EQ.'ALLPOINTS') THEN
1058:             PRINTPTS=.TRUE.1054:             PRINTPTS=.TRUE.
1059: 1055: 
1060: ! PCW: remove translations and rotations 
1061:  
1062:          ELSE IF (WORD.EQ.'NOTRANSROT') THEN 
1063:             NOTRANSROTT=.TRUE. 
1064:  
1065: ! davidg: introduced userpott here1056: ! davidg: introduced userpott here
1066:          ELSE IF (WORD.EQ.'USERPOT') THEN1057:          ELSE IF (WORD.EQ.'USERPOT') THEN
1067:             USERPOTT=.TRUE.1058:             USERPOTT=.TRUE.
1068:             RETURN1059:             RETURN
1069: ! MCP1060: ! MCP
1070:          ELSE IF (WORD.EQ.'AMH') THEN1061:          ELSE IF (WORD.EQ.'AMH') THEN
1071:             WRITE(6,*)'USING AMH ENERGIES FORCES'1062:             WRITE(6,*)'USING AMH ENERGIES FORCES'
1072:             WRITE(6,*)'CALCULATE ENERGY AND FORCE TABLES  '1063:             WRITE(6,*)'CALCULATE ENERGY AND FORCE TABLES  '
1073:             AMHT=.TRUE.1064:             AMHT=.TRUE.
1074:             WRITE(6,*)'AMH FLAG ', AMHT1065:             WRITE(6,*)'AMH FLAG ', AMHT


r30776/mindist.f 2016-07-15 07:30:15.355217019 +0100 r30775/mindist.f 2016-07-15 07:30:19.267269666 +0100
 19: C 19: C
 20: C  Finds the minimum distance between two geometries using LBFGS 20: C  Finds the minimum distance between two geometries using LBFGS
 21: C  Geometry R0 is reset to R1 after each optimisation step. Geometry 21: C  Geometry R0 is reset to R1 after each optimisation step. Geometry
 22: C  in R should not change, so neither should RA. RB is returned as the 22: C  in R should not change, so neither should RA. RB is returned as the
 23: C  closest geometry to RA. 23: C  closest geometry to RA.
 24: C 24: C
 25: C jmc As long as zsym isn't 'W' (in which case mind does something special) mind 25: C jmc As long as zsym isn't 'W' (in which case mind does something special) mind
 26: C doesn't care what atomic symbol we give it 26: C doesn't care what atomic symbol we give it
 27: C 27: C
 28:       SUBROUTINE MINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET) 28:       SUBROUTINE MINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET)
 29:       USE KEY,ONLY : MIEFT, NOTRANSROTT,FREEZE, PULLT, EFIELDT, GTHOMSONT ! hk286 29:       USE KEY,ONLY : MIEFT, FREEZE, PULLT, EFIELDT, GTHOMSONT ! hk286
 30:       IMPLICIT NONE 30:       IMPLICIT NONE
 31:       INTEGER J1, IG, I, ITER, J2, NCOUNT 31:       INTEGER J1, IG, I, ITER, J2, NCOUNT
 32:       DOUBLE PRECISION DPRAND 32:       DOUBLE PRECISION DPRAND
 33:       DOUBLE PRECISION P(3), DIST, DIST0, RG, RG0, DISTFUNC,  33:       DOUBLE PRECISION P(3), DIST, DIST0, RG, RG0, DISTFUNC, 
 34:      1                 MYROTMAT(3,3),  34:      1                 MYROTMAT(3,3), 
 35:      2                 OVEC(3), H1VEC(3), H2VEC(3), DSAVE, OMEGATOT(3,3), RA(*), RB(*) 35:      2                 OVEC(3), H1VEC(3), H2VEC(3), DSAVE, OMEGATOT(3,3), RA(*), RB(*)
 36:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: R, R0, R1, R1SAVE 36:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: R, R0, R1, R1SAVE
 37:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: RBSAVE 37:       DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: RBSAVE
 38:       INTEGER NSIZE, NATOMS 38:       INTEGER NSIZE, NATOMS
 39:       LOGICAL BULKT, MFLAG, TWOD, AGAIN, PRESERVET 39:       LOGICAL BULKT, MFLAG, TWOD, AGAIN, PRESERVET
102:             R(2,J1)=RA(3*(J1-1)+2)102:             R(2,J1)=RA(3*(J1-1)+2)
103:             R(3,J1)=RA(3*(J1-1)+3)103:             R(3,J1)=RA(3*(J1-1)+3)
104:             R0(1,J1)=RB(3*(J1-1)+1)104:             R0(1,J1)=RB(3*(J1-1)+1)
105:             R0(2,J1)=RB(3*(J1-1)+2)105:             R0(2,J1)=RB(3*(J1-1)+2)
106:             R0(3,J1)=RB(3*(J1-1)+3)106:             R0(3,J1)=RB(3*(J1-1)+3)
107:          ENDDO107:          ENDDO
108:       ENDIF108:       ENDIF
109: C109: C
110: C     put c.o.m. to origin110: C     put c.o.m. to origin
111: C111: C
112:       IF (.NOT.FREEZE .AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN112:       IF (.NOT.FREEZE .AND..NOT.MIEFT) THEN
113:          do ig=1,3113:          do ig=1,3
114:             rg=0.d0114:             rg=0.d0
115:             rg0=0.d0115:             rg0=0.d0
116:             do i=1,nsize116:             do i=1,nsize
117:                rg=rg+R(ig,i)117:                rg=rg+R(ig,i)
118:                rg0=rg0+R0(ig,i)118:                rg0=rg0+R0(ig,i)
119:             enddo119:             enddo
120:             rg=rg/nsize120:             rg=rg/nsize
121:             rg0=rg0/nsize121:             rg0=rg0/nsize
122:             do i=1,nsize122:             do i=1,nsize
157:          IF (NCOUNT.GT.0) THEN 157:          IF (NCOUNT.GT.0) THEN 
158:             PRINT*,'convergence failure in mind'158:             PRINT*,'convergence failure in mind'
159: c           STOP159: c           STOP
160:          ELSE160:          ELSE
161: C           WRITE(*,'(A,2F15.5,A,I6,A,I6)') 'Initial and final distances:',DIST0,DIST,' iterations=',ITER,' NCOUNT=',NCOUNT161: C           WRITE(*,'(A,2F15.5,A,I6,A,I6)') 'Initial and final distances:',DIST0,DIST,' iterations=',ITER,' NCOUNT=',NCOUNT
162:             DO J1=1,NSIZE162:             DO J1=1,NSIZE
163:                R0(1,J1)=R1(1,J1)163:                R0(1,J1)=R1(1,J1)
164:                R0(2,J1)=R1(2,J1)164:                R0(2,J1)=R1(2,J1)
165:                R0(3,J1)=R1(3,J1)165:                R0(3,J1)=R1(3,J1)
166:             ENDDO166:             ENDDO
167:             IF (.NOT.(TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT)) P(1)=2*(DPRAND()-0.5D0)/10.0D0167:             IF (.NOT.(TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT)) P(1)=2*(DPRAND()-0.5D0)/10.0D0
168:             IF (.NOT.(TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT)) P(2)=2*(DPRAND()-0.5D0)/10.0D0168:             IF (.NOT.(TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT)) P(2)=2*(DPRAND()-0.5D0)/10.0D0
169:             P(3)=2*(DPRAND()-0.5D0)/10.0D0169:             P(3)=2*(DPRAND()-0.5D0)/10.0D0
170:             GOTO 10170:             GOTO 10
171:          ENDIF171:          ENDIF
172:       ENDIF172:       ENDIF
173: C173: C
174: C  This block allows the second geometry to rotate out of the xy plane; not allowed!174: C  This block allows the second geometry to rotate out of the xy plane; not allowed!
175: C175: C
176: C     IF (TWOD) THEN176: C     IF (TWOD) THEN
177: C        IF (AGAIN) THEN177: C        IF (AGAIN) THEN
178: C           WRITE(*,'(A,2F20.10)') 'DIST,DSAVE=',DIST,DSAVE178: C           WRITE(*,'(A,2F20.10)') 'DIST,DSAVE=',DIST,DSAVE
258:             RA(3*(J1-1)+3)=R(3,J1)258:             RA(3*(J1-1)+3)=R(3,J1)
259:          ENDDO259:          ENDDO
260:       ENDIF260:       ENDIF
261: 261: 
262:       DO J1=1,3262:       DO J1=1,3
263:          DO J2=1,3263:          DO J2=1,3
264:             MYROTMAT(J2,J1)=OMEGATOT(J2,J1)264:             MYROTMAT(J2,J1)=OMEGATOT(J2,J1)
265:          ENDDO265:          ENDDO
266:       ENDDO266:       ENDDO
267: 267: 
268:       IF (PRESERVET.OR.MIEFT .OR. NOTRANSROTT) RB(1:3*NATOMS)=RBSAVE(1:3*NATOMS)268:       IF (PRESERVET.OR.MIEFT) RB(1:3*NATOMS)=RBSAVE(1:3*NATOMS)
269: 269: 
270:       DEALLOCATE(R, R0, R1, R1SAVE, RBSAVE)270:       DEALLOCATE(R, R0, R1, R1SAVE, RBSAVE)
271: 271: 
272:       RETURN272:       RETURN
273:       END273:       END
274: 274: 
275: c__________________________________________________________________________275: c__________________________________________________________________________
276: 276: 
277:         subroutine myinverse(A,B)277:         subroutine myinverse(A,B)
278:         IMPLICIT NONE278:         IMPLICIT NONE
728: 728: 
729:       PTEST=.TRUE.729:       PTEST=.TRUE.
730:       PTEST=.FALSE.730:       PTEST=.FALSE.
731:       ALPHA=1.0D0731:       ALPHA=1.0D0
732:       NFAIL=0732:       NFAIL=0
733:       ITER=0733:       ITER=0
734:       ITDONE=0734:       ITDONE=0
735:       ENERGY=DISTFUNC(X,R,R0,R1)735:       ENERGY=DISTFUNC(X,R,R0,R1)
736:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)736:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)
737: ! hk286 737: ! hk286 
738:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0738:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0
739:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0739:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0
740:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)740:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)
741:       DO J1=1,N741:       DO J1=1,N
742:          G(J1)=GSAVE(J1)742:          G(J1)=GSAVE(J1)
743:       ENDDO743:       ENDDO
744: 744: 
745:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,' LBFGS steps'745:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,' LBFGS steps'
746: 746: 
747: 10    CALL FLUSH(6)747: 10    CALL FLUSH(6)
748:       MFLAG=.FALSE.748:       MFLAG=.FALSE.
749:       IF (RMS.LE.EPS) THEN749:       IF (RMS.LE.EPS) THEN
859:       ENDIF859:       ENDIF
860: C860: C
861: C  Store the new search direction861: C  Store the new search direction
862: C862: C
863:       IF (ITER.GT.0) THEN863:       IF (ITER.GT.0) THEN
864:          DO I=1,N864:          DO I=1,N
865:             W(ISPT+POINT*N+I)= W(I)865:             W(ISPT+POINT*N+I)= W(I)
866:          ENDDO866:          ENDDO
867:       ENDIF867:       ENDIF
868: ! hk286868: ! hk286
869:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+1)=0.0D0869:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+1)=0.0D0
870:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+2)=0.0D0870:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+2)=0.0D0
871: 871: 
872: C     OVERLAP=DDOT(N,G,1,W,1)/SQRT(DDOT(N,G,1,G,1)*DDOT(N,W,1,W,1))872: C     OVERLAP=DDOT(N,G,1,W,1)/SQRT(DDOT(N,G,1,G,1)*DDOT(N,W,1,W,1))
873:       DOT1=SQRT(DDOT(N,G,1,G,1))873:       DOT1=SQRT(DDOT(N,G,1,G,1))
874:       DOT2=SQRT(DDOT(N,W,1,W,1))874:       DOT2=SQRT(DDOT(N,W,1,W,1))
875:       OVERLAP=0.0D0875:       OVERLAP=0.0D0
876:       IF (DOT1*DOT2.NE.0.0D0) OVERLAP=DDOT(N,G,1,W,1)/(DOT1*DOT2)876:       IF (DOT1*DOT2.NE.0.0D0) OVERLAP=DDOT(N,G,1,W,1)/(DOT1*DOT2)
877: 877: 
878: C     PRINT*,'OVERLAP,DIAG(1)=',OVERLAP,DIAG(1)878: C     PRINT*,'OVERLAP,DIAG(1)=',OVERLAP,DIAG(1)
879: C     PRINT*,'G . G=',DDOT(N,G,1,G,1)879: C     PRINT*,'G . G=',DDOT(N,G,1,G,1)
880: C     PRINT*,'W . W=',DDOT(N,W,1,W,1)880: C     PRINT*,'W . W=',DDOT(N,W,1,W,1)
960:       X(2)=0.0D0960:       X(2)=0.0D0
961:       X(3)=0.0D0961:       X(3)=0.0D0
962:       DO I=1,NSIZE962:       DO I=1,NSIZE
963:          DO J=1,3963:          DO J=1,3
964:             R0(J,I)=R1(J,I)964:             R0(J,I)=R1(J,I)
965:          ENDDO965:          ENDDO
966:       ENDDO966:       ENDDO
967: 967: 
968:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)968:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)
969: ! hk286969: ! hk286
970:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0970:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0
971:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.NOTRANSROTT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0971:       IF (TWOD.OR.PULLT.OR.EFIELDT.OR.MIEFT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0
972:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)972:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)
973:       DO J1=1,N973:       DO J1=1,N
974:          G(J1)=GSAVE(J1)974:          G(J1)=GSAVE(J1)
975:       ENDDO975:       ENDDO
976:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A,G15.5)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,976:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A,G15.5)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,
977:      1        ' LBFGS steps, step:',STP*SLENGTH977:      1        ' LBFGS steps, step:',STP*SLENGTH
978: C978: C
979: C     Compute the new step and gradient change979: C     Compute the new step and gradient change
980: C980: C
981: 30    NPT=POINT*N981: 30    NPT=POINT*N


r30776/minpermdist.f90 2016-07-15 07:30:15.603220356 +0100 r30775/minpermdist.f90 2016-07-15 07:30:19.519273049 +0100
 48: !  +/- RMATCUMUL ROTA (COORDSA - CMA) = permutation(DUMMYA) 48: !  +/- RMATCUMUL ROTA (COORDSA - CMA) = permutation(DUMMYA)
 49: !  where +/- is given by the value of INVERT. 49: !  where +/- is given by the value of INVERT.
 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the
 51: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we 51: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we
 52: !  are doing an ion trap potential. 52: !  are doing an ion trap potential.
 53: ! 53: !
 54: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST) 54: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 55: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, STOCKT, GEOMDIFFTOL, AMBERT, & 55: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, STOCKT, GEOMDIFFTOL, AMBERT, &
 56:   &            NFREEZE, NABT, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, PULLT, EFIELDT, NTSITES, & 56:   &            NFREEZE, NABT, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, PULLT, EFIELDT, NTSITES, &
 57:   &            RIGIDBODY, PERMDIST, OHCELLT, LPERMDIST, EYTRAPT, MKTRAPT, LOCALPERMCUT, LOCALPERMCUT2, & 57:   &            RIGIDBODY, PERMDIST, OHCELLT, LPERMDIST, EYTRAPT, MKTRAPT, LOCALPERMCUT, LOCALPERMCUT2, &
 58:   &            LOCALPERMCUTINC, NOINVERSION, MIEFT, NOTRANSROTT, & 58:   &            LOCALPERMCUTINC, NOINVERSION, MIEFT, &
 59:   &            EDIFFTOL, GMAX, CONVR, ATOMMATCHDIST, NRANROT, GTHOMSONT, GTHOMMET, & ! hk286 59:   &            EDIFFTOL, GMAX, CONVR, ATOMMATCHDIST, NRANROT, GTHOMSONT, GTHOMMET, & ! hk286
 60:   &            PHI4MODT, MCPATHT, AMBER12T, VARIABLES, MKTRAPT, ALIGNRBST, QCIAMBERT 60:   &            PHI4MODT, MCPATHT, AMBER12T, VARIABLES, MKTRAPT, ALIGNRBST, QCIAMBERT
 61: USE COMMONS,ONLY : NOPT 61: USE COMMONS,ONLY : NOPT
 62: USE MODCHARMM,ONLY : CHRMMT 62: USE MODCHARMM,ONLY : CHRMMT
 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH
 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET
 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE
 66: USE GENRIGID 66: USE GENRIGID
 67: USE AMBER12_INTERFACE_MOD 67: USE AMBER12_INTERFACE_MOD
 68: USE CHIRALITY 68: USE CHIRALITY
328:       IF (DEBUG) PRINT*, "msb50 minpermdist using intdistance", DISTANCE328:       IF (DEBUG) PRINT*, "msb50 minpermdist using intdistance", DISTANCE
329:     ENDIF329:     ENDIF
330:    DISTANCE=SQRT(DISTANCE)330:    DISTANCE=SQRT(DISTANCE)
331:    RETURN331:    RETURN
332: 332: 
333: ENDIF333: ENDIF
334: !334: !
335: !  Calculate original centres of mass. sn402 - this is actually the centre of geometry?335: !  Calculate original centres of mass. sn402 - this is actually the centre of geometry?
336: !336: !
337: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0337: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0
338: IF ((NFREEZE.LE.0).AND.(.NOT.MIEFT).AND.(.NOT.MKTRAPT) .AND. (.NOT. NOTRANSROTT)) THEN338: IF ((NFREEZE.LE.0).AND.(.NOT.MIEFT).AND.(.NOT.MKTRAPT)) THEN
339:    IF (STOCKT) THEN 339:    IF (STOCKT) THEN 
340:       NRB=(NATOMS/2)340:       NRB=(NATOMS/2)
341:       DO J1=1,NRB341:       DO J1=1,NRB
342:          CMAX=CMAX+COORDSA(3*(J1-1)+1)342:          CMAX=CMAX+COORDSA(3*(J1-1)+1)
343:          CMAY=CMAY+COORDSA(3*(J1-1)+2)343:          CMAY=CMAY+COORDSA(3*(J1-1)+2)
344:          CMAZ=CMAZ+COORDSA(3*(J1-1)+3)344:          CMAZ=CMAZ+COORDSA(3*(J1-1)+3)
345:       ENDDO345:       ENDDO
346:       CMAX=CMAX/NRB; CMAY=CMAY/NRB; CMAZ=CMAZ/NRB346:       CMAX=CMAX/NRB; CMAY=CMAY/NRB; CMAZ=CMAZ/NRB
347:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0347:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0
348:       DO J1=1,NRB348:       DO J1=1,NRB
376: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)376: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)
377: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)377: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)
378: DBEST=1.0D100378: DBEST=1.0D100
379: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)379: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
380: DBEST=DISTANCE**2380: DBEST=DISTANCE**2
381: 381: 
382: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE382: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE
383: 383: 
384: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then384: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then
385: ! translate the origin of structure A to the CoM of structure B.385: ! translate the origin of structure A to the CoM of structure B.
386: IF ((NFREEZE.GT.0).OR.BULKT.OR.MIEFT.OR.MKTRAPT.OR.NOTRANSROTT) THEN386: IF ((NFREEZE.GT.0).OR.BULKT.OR.MIEFT.OR.MKTRAPT) THEN
387:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)387:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
388: ELSE IF (STOCKT) THEN388: ELSE IF (STOCKT) THEN
389:    DO J1=1,(NATOMS/2)389:    DO J1=1,(NATOMS/2)
390:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX390:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX
391:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY391:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY
392:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ392:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ
393:    ENDDO393:    ENDDO
394: ELSE394: ELSE
395:    DO J1=1,NATOMS395:    DO J1=1,NATOMS
396:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX396:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX
453: ! because of the numerical cutoffs employed in MYORIENT. We could miss the453: ! because of the numerical cutoffs employed in MYORIENT. We could miss the
454: ! right orientation! 454: ! right orientation! 
455: !455: !
456: ! If we use MYORIENT to produce particular orientations then we end up aligning 456: ! If we use MYORIENT to produce particular orientations then we end up aligning 
457: ! COORDSA not with COORDSB but with the standard orientation of COORDSB in DUMMYB.457: ! COORDSA not with COORDSB but with the standard orientation of COORDSB in DUMMYB.
458: ! We now deal with this by tracking the complete transformation, including the458: ! We now deal with this by tracking the complete transformation, including the
459: ! contribution of MYORIENT using ROTB and ROTINVB.459: ! contribution of MYORIENT using ROTB and ROTINVB.
460: !460: !
461: 461: 
462: DISTANCE=0.0D0462: DISTANCE=0.0D0
463: IF (NFREEZE.LE.0 .AND. .NOT.MIEFT .AND. .NOT. NOTRANSROTT) THEN463: IF (NFREEZE.LE.0 .AND. .NOT.MIEFT) THEN
464:    IF (BULKT) THEN464:    IF (BULKT) THEN
465:       IF (BMTEST) BMCOORDSSV(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)465:       IF (BMTEST) BMCOORDSSV(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)
466:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;466:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;
467:       CALL BULKMINDIST(DUMMYB,DUMMYA,BMCOORDS,NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,.TRUE., TNMATCH, BMTEST)467:       CALL BULKMINDIST(DUMMYB,DUMMYA,BMCOORDS,NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,.TRUE., TNMATCH, BMTEST)
468:       IF (PITEST) THEN468:       IF (PITEST) THEN
469:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)469:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
470:          RMATBEST(1:3,1:3)=0.0D0470:          RMATBEST(1:3,1:3)=0.0D0
471:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0471:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
472:          DISTANCE=SQRT(DISTANCE)472:          DISTANCE=SQRT(DISTANCE)
473:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST, distance=',DISTANCE473:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST, distance=',DISTANCE
831: IF (NCHOOSEB2.LT.NORBITB2) GOTO 31831: IF (NCHOOSEB2.LT.NORBITB2) GOTO 31
832: IF (NCHOOSEB1.LT.NORBITB1) GOTO 66832: IF (NCHOOSEB1.LT.NORBITB1) GOTO 66
833: !833: !
834: !  Now try the enantiomer (or xz reflected structure for PULLT.OR.EFIELDT.OR.TWOD).834: !  Now try the enantiomer (or xz reflected structure for PULLT.OR.EFIELDT.OR.TWOD).
835: !  The tests for NCHOOSE1 and NCHOOSE2 appear to be redundant!835: !  The tests for NCHOOSE1 and NCHOOSE2 appear to be redundant!
836: !836: !
837: IF ((NCHOOSE2.EQ.NORBIT2).AND.(NCHOOSE1.EQ.NORBIT1).AND.(INVERT.EQ.1)) THEN837: IF ((NCHOOSE2.EQ.NORBIT2).AND.(NCHOOSE1.EQ.NORBIT1).AND.(INVERT.EQ.1)) THEN
838: !838: !
839: ! don't try inversion for bulk or charmm or amber or frozen atoms839: ! don't try inversion for bulk or charmm or amber or frozen atoms
840: !840: !
841:    IF (BULKT.OR.CHRMMT.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.QCIAMBERT.OR.(NFREEZE.GT.0).OR.NOINVERSION.OR.MIEFT.OR.NOTRANSROTT) GOTO 50 841:    IF (BULKT.OR.CHRMMT.OR.AMBERT.OR.NABT.OR.AMBER12T.OR.QCIAMBERT.OR.(NFREEZE.GT.0).OR.NOINVERSION.OR.MIEFT) GOTO 50 
842: !  IF (DEBUG) PRINT '(A)',' minpermdist> inverting geometry for comparison with target'842: !  IF (DEBUG) PRINT '(A)',' minpermdist> inverting geometry for comparison with target'
843:    INVERT=-1843:    INVERT=-1
844:    GOTO 60844:    GOTO 60
845: ENDIF845: ENDIF
846: IF (NROTDONE.LT.NRANROT) GOTO 11846: IF (NROTDONE.LT.NRANROT) GOTO 11
847: 847: 
848: 50 DISTANCE=DBEST848: 50 DISTANCE=DBEST
849: !849: !
850: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.850: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.
851: !  Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, 851: !  Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, 
852: !  undoing the reorientation to DUMMYB from MYORIENT. 852: !  undoing the reorientation to DUMMYB from MYORIENT. 
853: !  We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA) 853: !  We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA) 
854: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 854: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 
855: !  (aside from a possible permutation of the atom ordering)855: !  (aside from a possible permutation of the atom ordering)
856: !856: !
857:    IF ((NFREEZE.GT.0).OR.MIEFT.OR.MKTRAPT .OR. NOTRANSROTT) THEN857:    IF ((NFREEZE.GT.0).OR.MIEFT.OR.MKTRAPT) THEN
858:       XDUMMY=0.0D0858:       XDUMMY=0.0D0
859:       DO J1=1,NATOMS859:       DO J1=1,NATOMS
860:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &860:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
861:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &861:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
862:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2862:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
863:       ENDDO863:       ENDDO
864:    ELSE IF (STOCKT) THEN864:    ELSE IF (STOCKT) THEN
865:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)865:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)
866:       XDUMMY=0.0D0866:       XDUMMY=0.0D0
867:       DO J1=1,(NATOMS/2)867:       DO J1=1,(NATOMS/2)
891:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &891:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
892:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2892:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
893:       ENDDO893:       ENDDO
894:    ENDIF894:    ENDIF
895:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN895:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN
896:       PRINT '(2(A,F20.10))',' minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &896:       PRINT '(2(A,F20.10))',' minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &
897:   &                         ' should be ',SQRT(DISTANCE)897:   &                         ' should be ',SQRT(DISTANCE)
898:       STOP898:       STOP
899:    ENDIF899:    ENDIF
900: 900: 
901:    IF ((NFREEZE.GT.0).OR.MIEFT.OR.MKTRAPT.OR.NOTRANSROTT) THEN901:    IF ((NFREEZE.GT.0).OR.MIEFT.OR.MKTRAPT) THEN
902: !no rotation for NFREEZE .gt. 0902: !no rotation for NFREEZE .gt. 0
903:       RMATBEST(1:3,1:3)=0.0D0903:       RMATBEST(1:3,1:3)=0.0D0
904:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0904:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
905:    ELSE905:    ELSE
906:       RMATBEST=MATMUL(ROTINVBBEST,RMATBEST)906:       RMATBEST=MATMUL(ROTINVBBEST,RMATBEST)
907:    ENDIF907:    ENDIF
908: !!!!!!!!!!!!!!!!!!!!!!! DEBUG908: !!!!!!!!!!!!!!!!!!!!!!! DEBUG
909: !909: !
910: ! Test distance for COORDSA with permutation applied in BESTPERM910: ! Test distance for COORDSA with permutation applied in BESTPERM
911: !911: !


r30776/ncutils.f90 2016-07-15 07:30:11.927170889 +0100 r30775/ncutils.f90 2016-07-15 07:30:16.619234027 +0100
595: ! as atoms in the following subroutine. 595: ! as atoms in the following subroutine. 
596: ! A good solution would probably be to define an unambiguous sense for the596: ! A good solution would probably be to define an unambiguous sense for the
597: ! dipoles so that this problem doesn't arise. It could also cause problems597: ! dipoles so that this problem doesn't arise. It could also cause problems
598: ! in recognising identical minima.598: ! in recognising identical minima.
599: !599: !
600: ! Merges path output files to produce full pathway for the rearrangement;600: ! Merges path output files to produce full pathway for the rearrangement;
601: ! frames in path are reversed as needed;601: ! frames in path are reversed as needed;
602: !602: !
603:      SUBROUTINE MERGEXYZEOFS  603:      SUBROUTINE MERGEXYZEOFS  
604:           USE KEY, ONLY: FILTH,UNRST,FILTHSTR,RIGIDBODY,BULKT,TWOD,STOCKT,STOCKAAT,RBAAT,PERMDIST, MIEFT, &604:           USE KEY, ONLY: FILTH,UNRST,FILTHSTR,RIGIDBODY,BULKT,TWOD,STOCKT,STOCKAAT,RBAAT,PERMDIST, MIEFT, &
605:   &                      AMHT,SEQ,NTSITES,NENDDUP, GTHOMSONT, PAPT, NFREEZE, VARIABLES, NOTRANSROTT ! hk286605:   &                      AMHT,SEQ,NTSITES,NENDDUP, GTHOMSONT, PAPT, NFREEZE, VARIABLES ! hk286
606:           USE KEYUTILS        ! frames in bits that are glued together are rotated accordingly;606:           USE KEYUTILS        ! frames in bits that are glued together are rotated accordingly;
607:           USE KEYCONNECT,ONLY : NCNSTEPS607:           USE KEYCONNECT,ONLY : NCNSTEPS
608:           USE COMMONS,ONLY : PARAM1,PARAM2,PARAM3,DEBUG608:           USE COMMONS,ONLY : PARAM1,PARAM2,PARAM3,DEBUG
609:           USE AMHGLOBALS, ONLY : NMRES609:           USE AMHGLOBALS, ONLY : NMRES
610: 610: 
611:           IMPLICIT NONE       ! prerequisites: chain of min/ts constructed; assumes path is dumping plus side of the path611:           IMPLICIT NONE       ! prerequisites: chain of min/ts constructed; assumes path is dumping plus side of the path
612:                               ! first, and there are no blank lines after last frame (!)612:                               ! first, and there are no blank lines after last frame (!)
613:                               ! does somewhat similar with EofS.ts files as well..613:                               ! does somewhat similar with EofS.ts files as well..
614:           DOUBLE PRECISION RMAT(3,3), Q2(4)614:           DOUBLE PRECISION RMAT(3,3), Q2(4)
615:           INTEGER :: I,J,K,EOF,J1,J2,INDEXTS !,FL615:           INTEGER :: I,J,K,EOF,J1,J2,INDEXTS !,FL
895:                            CALL MINPERMDIST(LASTFRAME,TMP%Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)895:                            CALL MINPERMDIST(LASTFRAME,TMP%Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
896:                         ENDIF896:                         ENDIF
897:                      ELSE897:                      ELSE
898:                         IF (RBAAT) THEN898:                         IF (RBAAT) THEN
899:                            CALL RBMINDIST(LASTFRAME,TMP%Q,NATOMS,DIST,Q2,DEBUG)899:                            CALL RBMINDIST(LASTFRAME,TMP%Q,NATOMS,DIST,Q2,DEBUG)
900:                            CALL QROTMAT (Q2, RMATBEST)900:                            CALL QROTMAT (Q2, RMATBEST)
901:                         ELSEIF (VARIABLES) THEN901:                         ELSEIF (VARIABLES) THEN
902:                         ELSE902:                         ELSE
903:                            CALL NEWMINDIST(LASTFRAME,TMP%Q,NATOMS,D,BULKT,TWOD,TMP%SYM(1),.TRUE.,RIGIDBODY,DEBUG,RMAT)903:                            CALL NEWMINDIST(LASTFRAME,TMP%Q,NATOMS,D,BULKT,TWOD,TMP%SYM(1),.TRUE.,RIGIDBODY,DEBUG,RMAT)
904:                            CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA)904:                            CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA)
905:                            IF(.NOT.MIEFT .AND. .NOT. NOTRANSROTT ) THEN905:                            IF(.NOT.MIEFT) THEN
906:                               DO J1=1,NATOMS906:                               DO J1=1,NATOMS
907:                                  J2=3*J1907:                                  J2=3*J1
908:                                  TMP%Q(J2-2)=TMP%Q(J2-2)+CMXFIX-CMXA908:                                  TMP%Q(J2-2)=TMP%Q(J2-2)+CMXFIX-CMXA
909:                                  TMP%Q(J2-1)=TMP%Q(J2-1)+CMYFIX-CMYA909:                                  TMP%Q(J2-1)=TMP%Q(J2-1)+CMYFIX-CMYA
910:                                  TMP%Q(J2)  =TMP%Q(J2)  +CMZFIX-CMZA910:                                  TMP%Q(J2)  =TMP%Q(J2)  +CMZFIX-CMZA
911:                               ENDDO911:                               ENDDO
912:                            ENDIF912:                            ENDIF
913:                         ENDIF913:                         ENDIF
914:                      ENDIF914:                      ENDIF
915:                   ENDIF915:                   ENDIF
986:                         IF ((I.EQ.1).AND.(K.EQ.1).AND.(NENDDUP.GT.0)) THEN986:                         IF ((I.EQ.1).AND.(K.EQ.1).AND.(NENDDUP.GT.0)) THEN
987:                            IF (DEBUG) PRINT '(A,3I6)','ncutils> writing NENDDUP extra start frames'987:                            IF (DEBUG) PRINT '(A,3I6)','ncutils> writing NENDDUP extra start frames'
988:                            DO J2=1,NENDDUP988:                            DO J2=1,NENDDUP
989:                               CALL WRITEFRAME(TMP%COMMENT,TMP%SYM,TMP%Q)989:                               CALL WRITEFRAME(TMP%COMMENT,TMP%SYM,TMP%Q)
990:                            ENDDO990:                            ENDDO
991:                         ENDIF991:                         ENDIF
992:                      ELSE992:                      ELSE
993:                         IF ((I>1.AND.K>1) .AND. (.NOT.(UNRST.OR.VARIABLES))) THEN993:                         IF ((I>1.AND.K>1) .AND. (.NOT.(UNRST.OR.VARIABLES))) THEN
994:                            IF (.NOT.BULKT) THEN994:                            IF (.NOT.BULKT) THEN
995:                               CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA)995:                               CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA)
996:                               IF(.NOT.MIEFT .AND. .NOT. NOTRANSROTT) THEN996:                               IF(.NOT.MIEFT) THEN
997:                                  DO J1=1,NATOMS997:                                  DO J1=1,NATOMS
998:                                     J2=3*J1998:                                     J2=3*J1
999:                                     TMP%Q(J2-2)=TMP%Q(J2-2)+CMXFIX-CMXA999:                                     TMP%Q(J2-2)=TMP%Q(J2-2)+CMXFIX-CMXA
1000:                                     TMP%Q(J2-1)=TMP%Q(J2-1)+CMYFIX-CMYA1000:                                     TMP%Q(J2-1)=TMP%Q(J2-1)+CMYFIX-CMYA
1001:                                     TMP%Q(J2)  =TMP%Q(J2)  +CMZFIX-CMZA1001:                                     TMP%Q(J2)  =TMP%Q(J2)  +CMZFIX-CMZA
1002:                                  ENDDO1002:                                  ENDDO
1003:                               ENDIF1003:                               ENDIF
1004:                            ELSE1004:                            ELSE
1005:                               DO J2=1,NATOMS1005:                               DO J2=1,NATOMS
1006:                                  TMP%Q(3*(J2-1)+1)=TMP%Q(3*(J2-1)+1)-BSX-PARAM1*NINT((TMP%Q(3*(J2-1)+1)-BSX)/PARAM1)1006:                                  TMP%Q(3*(J2-1)+1)=TMP%Q(3*(J2-1)+1)-BSX-PARAM1*NINT((TMP%Q(3*(J2-1)+1)-BSX)/PARAM1)
1082:                            CALL MINPERMDIST(LASTFRAME,TMP%Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)1082:                            CALL MINPERMDIST(LASTFRAME,TMP%Q,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
1083:                         ENDIF1083:                         ENDIF
1084:                      ELSE1084:                      ELSE
1085:                         IF (RBAAT) THEN1085:                         IF (RBAAT) THEN
1086:                            CALL RBMINDIST(LASTFRAME,TMP%Q,NATOMS,DIST,Q2,DEBUG)1086:                            CALL RBMINDIST(LASTFRAME,TMP%Q,NATOMS,DIST,Q2,DEBUG)
1087:                            CALL QROTMAT (Q2, RMATBEST)1087:                            CALL QROTMAT (Q2, RMATBEST)
1088:                         ELSEIF (VARIABLES) THEN1088:                         ELSEIF (VARIABLES) THEN
1089:                         ELSE1089:                         ELSE
1090:                            CALL NEWMINDIST(LASTFRAME,TMP%Q,NATOMS,D,BULKT,TWOD,TMP%SYM(1),.TRUE.,RIGIDBODY,DEBUG,RMAT)1090:                            CALL NEWMINDIST(LASTFRAME,TMP%Q,NATOMS,D,BULKT,TWOD,TMP%SYM(1),.TRUE.,RIGIDBODY,DEBUG,RMAT)
1091:                            CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA)1091:                            CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA)
1092:                            IF (NFREEZE.LE.0 .AND..NOT.MIEFT .AND. .NOT. NOTRANSROTT) THEN1092:                            IF (NFREEZE.LE.0 .AND..NOT.MIEFT) THEN
1093:                               DO J1=1,NATOMS1093:                               DO J1=1,NATOMS
1094:                                  J2=3*J11094:                                  J2=3*J1
1095:                                  TMP%Q(J2-2)=TMP%Q(J2-2)+CMXFIX-CMXA1095:                                  TMP%Q(J2-2)=TMP%Q(J2-2)+CMXFIX-CMXA
1096:                                  TMP%Q(J2-1)=TMP%Q(J2-1)+CMYFIX-CMYA1096:                                  TMP%Q(J2-1)=TMP%Q(J2-1)+CMYFIX-CMYA
1097:                                  TMP%Q(J2)  =TMP%Q(J2)  +CMZFIX-CMZA1097:                                  TMP%Q(J2)  =TMP%Q(J2)  +CMZFIX-CMZA
1098:                               ENDDO1098:                               ENDDO
1099:                            ENDIF1099:                            ENDIF
1100:                         ENDIF1100:                         ENDIF
1101:                      ENDIF1101:                      ENDIF
1102:                      IF (LOCALSTOCK) STOCKT=.TRUE.1102:                      IF (LOCALSTOCK) STOCKT=.TRUE.
1172:                         IF ((I.EQ.1).AND.(K.EQ.1).AND.(NENDDUP.GT.0)) THEN1172:                         IF ((I.EQ.1).AND.(K.EQ.1).AND.(NENDDUP.GT.0)) THEN
1173:                            IF (DEBUG) PRINT '(A,3I6)','ncutils> writing NENDDUP extra start frames'1173:                            IF (DEBUG) PRINT '(A,3I6)','ncutils> writing NENDDUP extra start frames'
1174:                            DO J2=1,NENDDUP1174:                            DO J2=1,NENDDUP
1175:                               CALL WRITEFRAME(TMP%COMMENT,TMP%SYM,TMP%Q)1175:                               CALL WRITEFRAME(TMP%COMMENT,TMP%SYM,TMP%Q)
1176:                            ENDDO1176:                            ENDDO
1177:                         ENDIF1177:                         ENDIF
1178:                      ELSE1178:                      ELSE
1179:                         IF ((I>1.AND.K>1).AND.(.NOT.(UNRST.OR.VARIABLES))) THEN1179:                         IF ((I>1.AND.K>1).AND.(.NOT.(UNRST.OR.VARIABLES))) THEN
1180:                            IF (.NOT.BULKT) THEN1180:                            IF (.NOT.BULKT) THEN
1181:                               CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA) 1181:                               CALL NEWROTGEOM(NATOMS,TMP%Q,RMAT,CMXA,CMYA,CMZA) 
1182:                               IF (NFREEZE.LE.0.AND..NOT.MIEFT .AND. .NOT.NOTRANSROTT) THEN1182:                               IF (NFREEZE.LE.0.AND..NOT.MIEFT) THEN
1183:                                  DO J1=1,NATOMS1183:                                  DO J1=1,NATOMS
1184:                                     TMP%Q(3*(J1-1)+1)=TMP%Q(3*(J1-1)+1)-CMXA+CMXFIX1184:                                     TMP%Q(3*(J1-1)+1)=TMP%Q(3*(J1-1)+1)-CMXA+CMXFIX
1185:                                     TMP%Q(3*(J1-1)+2)=TMP%Q(3*(J1-1)+2)-CMYA+CMYFIX1185:                                     TMP%Q(3*(J1-1)+2)=TMP%Q(3*(J1-1)+2)-CMYA+CMYFIX
1186:                                     TMP%Q(3*(J1-1)+3)=TMP%Q(3*(J1-1)+3)-CMZA+CMZFIX1186:                                     TMP%Q(3*(J1-1)+3)=TMP%Q(3*(J1-1)+3)-CMZA+CMZFIX
1187:                                  ENDDO1187:                                  ENDDO
1188:                               ENDIF1188:                               ENDIF
1189:                            ELSE1189:                            ELSE
1190:                               IF (NFREEZE.LE.0.AND..NOT.MIEFT.AND. .NOT.NOTRANSROTT ) THEN1190:                               IF (NFREEZE.LE.0.AND..NOT.MIEFT) THEN
1191:                                  DO J2=1,NATOMS1191:                                  DO J2=1,NATOMS
1192:                                     TMP%Q(3*(J2-1)+1)=TMP%Q(3*(J2-1)+1)-BSX-PARAM1*NINT((TMP%Q(3*(J2-1)+1)-BSX)/PARAM1)1192:                                     TMP%Q(3*(J2-1)+1)=TMP%Q(3*(J2-1)+1)-BSX-PARAM1*NINT((TMP%Q(3*(J2-1)+1)-BSX)/PARAM1)
1193:                                     TMP%Q(3*(J2-1)+2)=TMP%Q(3*(J2-1)+2)-BSY-PARAM2*NINT((TMP%Q(3*(J2-1)+2)-BSY)/PARAM2)1193:                                     TMP%Q(3*(J2-1)+2)=TMP%Q(3*(J2-1)+2)-BSY-PARAM2*NINT((TMP%Q(3*(J2-1)+2)-BSY)/PARAM2)
1194:                                     IF (.NOT.TWOD) TMP%Q(3*(J2-1)+3)= &1194:                                     IF (.NOT.TWOD) TMP%Q(3*(J2-1)+3)= &
1195:      &                                                  TMP%Q(3*(J2-1)+3)-BSZ-PARAM3*NINT((TMP%Q(3*(J2-1)+3)-BSZ)/PARAM3)1195:      &                                                  TMP%Q(3*(J2-1)+3)-BSZ-PARAM3*NINT((TMP%Q(3*(J2-1)+3)-BSZ)/PARAM3)
1196:                                  ENDDO1196:                                  ENDDO
1197:                               ENDIF1197:                               ENDIF
1198:                            ENDIF1198:                            ENDIF
1199:                         ENDIF1199:                         ENDIF
1200:                         CALL WRITEFRAME(TMP%COMMENT,TMP%SYM,TMP%Q)1200:                         CALL WRITEFRAME(TMP%COMMENT,TMP%SYM,TMP%Q)


r30776/newmindist.f90 2016-07-15 07:30:15.859223800 +0100 r30775/newmindist.f90 2016-07-15 07:30:19.771276440 +0100
 25: ! 25: !
 26: ! jmc As long as zsym isn't 'W' (in which case mind does something special) mind 26: ! jmc As long as zsym isn't 'W' (in which case mind does something special) mind
 27: ! doesn't care what atomic symbol we give it. 27: ! doesn't care what atomic symbol we give it.
 28: ! 28: !
 29: !  If PRESERVET is false we put RB into best correspondence with RA. This involves 29: !  If PRESERVET is false we put RB into best correspondence with RA. This involves
 30: !  a translation to the same centre of coordinates, followed by a rotation about that 30: !  a translation to the same centre of coordinates, followed by a rotation about that
 31: !  centre. 31: !  centre.
 32: ! 32: !
 33: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT) 33: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT)
 34: USE COMMONS,ONLY : PARAM1, PARAM2, PARAM3, NOPT 34: USE COMMONS,ONLY : PARAM1, PARAM2, PARAM3, NOPT
 35: USE KEY,ONLY : STOCKT, NFREEZE, MIEFT, NOTRANSROTT, RBAAT, PULLT, EFIELDT, EYTRAPT, MKTRAPT, USERPOTT, & 35: USE KEY,ONLY : STOCKT, NFREEZE, MIEFT, RBAAT, PULLT, EFIELDT, EYTRAPT, MKTRAPT, USERPOTT, &
 36:   &            GTHOMSONT, VARIABLES! hk286 36:   &            GTHOMSONT, VARIABLES! hk286
 37:  37: 
 38: IMPLICIT NONE 38: IMPLICIT NONE
 39:  39: 
 40: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN 40: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN
 41: DOUBLE PRECISION RA(NOPT), RB(NOPT), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), & 41: DOUBLE PRECISION RA(NOPT), RB(NOPT), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), &
 42:   &              DIAG(4), TEMPA(3*NOPT), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, & 42:   &              DIAG(4), TEMPA(3*NOPT), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, &
 43:   &              MYROTMAT(3,3), OMEGATOT(3,3) 43:   &              MYROTMAT(3,3), OMEGATOT(3,3)
 44: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:) 44: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
 45: LOGICAL BULKT, TWOD, RIGIDBODY, PRESERVET, DEBUG, UPRETURN 45: LOGICAL BULKT, TWOD, RIGIDBODY, PRESERVET, DEBUG, UPRETURN
138:    XB(1:3*NSIZE)=RB(1:3*NSIZE)138:    XB(1:3*NSIZE)=RB(1:3*NSIZE)
139: ELSE139: ELSE
140:    ALLOCATE(XA(3*NATOMS),XB(3*NATOMS))140:    ALLOCATE(XA(3*NATOMS),XB(3*NATOMS))
141:    NSIZE=NATOMS141:    NSIZE=NATOMS
142:    XA(1:3*NATOMS)=RA(1:3*NATOMS)142:    XA(1:3*NATOMS)=RA(1:3*NATOMS)
143:    XB(1:3*NATOMS)=RB(1:3*NATOMS)143:    XB(1:3*NATOMS)=RB(1:3*NATOMS)
144: ENDIF144: ENDIF
145: !145: !
146: ! If there are frozen atoms then just calculate the distance and return.146: ! If there are frozen atoms then just calculate the distance and return.
147: !147: !
148: IF (NFREEZE.GT.0 .OR. MIEFT .OR. NOTRANSROTT) THEN148: IF (NFREEZE.GT.0 .OR. MIEFT) THEN
149:    DIST=0.0D0149:    DIST=0.0D0
150:    IF (BULKT) THEN150:    IF (BULKT) THEN
151:       BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3151:       BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3
152:       DO J1=1,NSIZE152:       DO J1=1,NSIZE
153:          DIST=DIST + (XA(3*(J1-1)+1)-XB(3*(J1-1)+1) - BOXLX*NINT((XA(3*(J1-1)+1)-XB(3*(J1-1)+1))/BOXLX))**2 &153:          DIST=DIST + (XA(3*(J1-1)+1)-XB(3*(J1-1)+1) - BOXLX*NINT((XA(3*(J1-1)+1)-XB(3*(J1-1)+1))/BOXLX))**2 &
154:    &               + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2) - BOXLY*NINT((XA(3*(J1-1)+2)-XB(3*(J1-1)+2))/BOXLY))**2 &154:    &               + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2) - BOXLY*NINT((XA(3*(J1-1)+2)-XB(3*(J1-1)+2))/BOXLY))**2 &
155:    &               + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3) - BOXLZ*NINT((XA(3*(J1-1)+3)-XB(3*(J1-1)+3))/BOXLZ))**2155:    &               + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3) - BOXLZ*NINT((XA(3*(J1-1)+3)-XB(3*(J1-1)+3))/BOXLZ))**2
156:       ENDDO156:       ENDDO
157:    ELSE157:    ELSE
158:       ! Note: this is the site-by-site cartesian distance, which is equivalent to the angle-axis158:       ! Note: this is the site-by-site cartesian distance, which is equivalent to the angle-axis
353:    RMAT(1,1)=Q1**2+Q2**2-Q3**2-Q4**2353:    RMAT(1,1)=Q1**2+Q2**2-Q3**2-Q4**2
354:    RMAT(1,2)=2*(Q2*Q3+Q1*Q4)354:    RMAT(1,2)=2*(Q2*Q3+Q1*Q4)
355:    RMAT(1,3)=2*(Q2*Q4-Q1*Q3)355:    RMAT(1,3)=2*(Q2*Q4-Q1*Q3)
356:    RMAT(2,1)=2*(Q2*Q3-Q1*Q4)356:    RMAT(2,1)=2*(Q2*Q3-Q1*Q4)
357:    RMAT(2,2)=Q1**2+Q3**2-Q2**2-Q4**2357:    RMAT(2,2)=Q1**2+Q3**2-Q2**2-Q4**2
358:    RMAT(2,3)=2*(Q3*Q4+Q1*Q2)358:    RMAT(2,3)=2*(Q3*Q4+Q1*Q2)
359:    RMAT(3,1)=2*(Q2*Q4+Q1*Q3)359:    RMAT(3,1)=2*(Q2*Q4+Q1*Q3)
360:    RMAT(3,2)=2*(Q3*Q4-Q1*Q2)360:    RMAT(3,2)=2*(Q3*Q4-Q1*Q2)
361:    RMAT(3,3)=Q1**2+Q4**2-Q2**2-Q3**2361:    RMAT(3,3)=Q1**2+Q4**2-Q2**2-Q3**2
362: ENDIF362: ENDIF
363: IF (.NOT.(PRESERVET.OR.MIEFT.OR.NOTRANSROTT)) THEN ! never translate with substrate!363: 
 364: IF (.NOT.(PRESERVET.OR.MIEFT)) THEN ! never translate with substrate!
364:    IF (ZUSE(1:1).EQ.'W') THEN365:    IF (ZUSE(1:1).EQ.'W') THEN
365: !366: !
366: !  Translate the XB coordinates to the centre of coordinates of XA.367: !  Translate the XB coordinates to the centre of coordinates of XA.
367: !368: !
368:       DO J1=1,NSIZE369:       DO J1=1,NSIZE
369:          XB(3*(J1-1)+1)=XB(3*(J1-1)+1)+CMXA+XSHIFT370:          XB(3*(J1-1)+1)=XB(3*(J1-1)+1)+CMXA+XSHIFT
370:          XB(3*(J1-1)+2)=XB(3*(J1-1)+2)+CMYA+YSHIFT371:          XB(3*(J1-1)+2)=XB(3*(J1-1)+2)+CMYA+YSHIFT
371:          XB(3*(J1-1)+3)=XB(3*(J1-1)+3)+CMZA+ZSHIFT372:          XB(3*(J1-1)+3)=XB(3*(J1-1)+3)+CMZA+ZSHIFT
372:       ENDDO373:       ENDDO
373: !374: !


r30776/OPTIM.F 2016-07-15 07:30:12.195174504 +0100 r30775/OPTIM.F 2016-07-15 07:30:16.863237310 +0100
445:             WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV445:             WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV
446:          ELSE446:          ELSE
447:             WRITE(*,'(A)') ' OPTIM> z rotational ev shifting'447:             WRITE(*,'(A)') ' OPTIM> z rotational ev shifting'
448:          ENDIF448:          ENDIF
449:          SHIFTL(1)=SHIFTV449:          SHIFTL(1)=SHIFTV
450:          SHIFTL(2)=SHIFTV450:          SHIFTL(2)=SHIFTV
451:          SHIFTL(3)=SHIFTV451:          SHIFTL(3)=SHIFTV
452:          SHIFTL(4)=SHIFTV452:          SHIFTL(4)=SHIFTV
453:          SHIFTL(5)=SHIFTV453:          SHIFTL(5)=SHIFTV
454:          SHIFTL(6)=SHIFTV454:          SHIFTL(6)=SHIFTV
455:       ELSE IF (MIEFT.OR.MLP3T.OR.MLPB3T.OR.NOTRANSROTT) THEN455:       ELSE IF (MIEFT.OR.MLP3T.OR.MLPB3T) THEN
456:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'456:          WRITE(*,'(A,G20.10)') ' OPTIM> No eigenvalues shifted'
457:       ELSE457:       ELSE
458:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV458:          WRITE(*,'(A,G20.10)') ' OPTIM> Using translational/rotational ev shift=',SHIFTV
459:          DO J1=1,6459:          DO J1=1,6
460:             SHIFTL(J1)=SHIFTV460:             SHIFTL(J1)=SHIFTV
461:          ENDDO461:          ENDDO
462:       ENDIF462:       ENDIF
463: 463: 
464:       !js850> 464:       !js850> 
465:       IF (ZSYM(NATOMS).EQ.'SS') THEN465:       IF (ZSYM(NATOMS).EQ.'SS') THEN


r30776/RPH.f90 2016-07-15 07:30:13.343191178 +0100 r30775/RPH.f90 2016-07-15 07:30:17.111240646 +0100
310:    NEXMODES=3*NFREEZE310:    NEXMODES=3*NFREEZE
311: ENDIF311: ENDIF
312: IF (RBAAT) THEN312: IF (RBAAT) THEN
313:    IF (EFIELDT) THEN313:    IF (EFIELDT) THEN
314:       NEXMODES = 4314:       NEXMODES = 4
315:    ELSE315:    ELSE
316:       NEXMODES = 6316:       NEXMODES = 6
317:    ENDIF317:    ENDIF
318:    IF (STOCKAAT) NEXMODES = NEXMODES + NATOMS/2318:    IF (STOCKAAT) NEXMODES = NEXMODES + NATOMS/2
319: ENDIF319: ENDIF
320: IF(MIEFT.OR.NOTRANSROTT) NEXMODES = 0320: IF(MIEFT) NEXMODES = 0
321: CALL GETPROD(PROD,NEXMODES,EVALUES)321: CALL GETPROD(PROD,NEXMODES,EVALUES)
322: 322: 
323: KAPPA=3*NATOMS-6323: KAPPA=3*NATOMS-6
324: IF (TWOD) THEN324: IF (TWOD) THEN
325:    KAPPA=2*NATOMS-3325:    KAPPA=2*NATOMS-3
326:    IF (BULKT) KAPPA=2*NATOMS-2326:    IF (BULKT) KAPPA=2*NATOMS-2
327:       ELSE IF (PULLT) THEN327:       ELSE IF (PULLT) THEN
328:    KAPPA=3*NATOMS-4328:    KAPPA=3*NATOMS-4
329: ELSE IF (BULKT) THEN329: ELSE IF (BULKT) THEN
330:    KAPPA=3*NATOMS-3330:    KAPPA=3*NATOMS-3
331: ELSE IF (MKTRAPT.OR.MIEFT.OR.NOTRANSROTT) THEN331: ELSE IF (MKTRAPT.OR.MIEFT) THEN
332:    KAPPA=3*NATOMS332:    KAPPA=3*NATOMS
333: ELSE IF (EYTRAPT) THEN333: ELSE IF (EYTRAPT) THEN
334:    KAPPA=3*NATOMS-3334:    KAPPA=3*NATOMS-3
335: ELSE IF (FREEZE) THEN335: ELSE IF (FREEZE) THEN
336:    KAPPA=3*NATOMS-3*NFREEZE336:    KAPPA=3*NATOMS-3*NFREEZE
337: ELSE IF (RIGIDINIT) THEN 337: ELSE IF (RIGIDINIT) THEN 
338:    KAPPA=DEGFREEDOMS-6338:    KAPPA=DEGFREEDOMS-6
339: ELSE IF (UNRST) THEN339: ELSE IF (UNRST) THEN
340:    KAPPA=NINTS340:    KAPPA=NINTS
341: ENDIF341: ENDIF


r30776/secdiag.f90 2016-07-15 07:30:16.111227203 +0100 r30775/secdiag.f90 2016-07-15 07:30:20.031279939 +0100
 80: ! 80: !
 81: !  Must read VEC into LOCALV because we are going to play with the vector in 81: !  Must read VEC into LOCALV because we are going to play with the vector in
 82: !  question. This would mess up cases where we need to retry with a smaller 82: !  question. This would mess up cases where we need to retry with a smaller
 83: !  step size, because we cannot undo the previous step cleanly if the magnitude 83: !  step size, because we cannot undo the previous step cleanly if the magnitude
 84: !  of VEC is changed! 1/7/04 DJW 84: !  of VEC is changed! 1/7/04 DJW
 85: ! 85: !
 86:  86: 
 87:       LOCALV(1:NOPT)=VEC(1:NOPT) 87:       LOCALV(1:NOPT)=VEC(1:NOPT)
 88: !     PRINT '(A)','secdiag> vec before orthogopt' 88: !     PRINT '(A)','secdiag> vec before orthogopt'
 89: !     PRINT '(6G20.10)',LOCALV(1:NOPT) 89: !     PRINT '(6G20.10)',LOCALV(1:NOPT)
 90:       IF (NFREEZE.LT.3.AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN 90:       IF (NFREEZE.LT.3.AND..NOT.MIEFT) THEN
 91:         CALL ORTHOGOPT(LOCALV,COORDS,.TRUE.) 91:         CALL ORTHOGOPT(LOCALV,COORDS,.TRUE.)
 92:       ELSE 92:       ELSE
 93:          CALL VECNORM(LOCALV,NOPT) 93:          CALL VECNORM(LOCALV,NOPT)
 94:       ENDIF 94:       ENDIF
 95: !     PRINT '(A)','secdiag> vec after orthogopt' 95: !     PRINT '(A)','secdiag> vec after orthogopt'
 96: !     PRINT '(6G20.10)',LOCALV(1:NOPT) 96: !     PRINT '(6G20.10)',LOCALV(1:NOPT)
 97:  97: 
 98:       IF (FREEZE) THEN 98:       IF (FREEZE) THEN
 99:          DO J1=1,NATOMS 99:          DO J1=1,NATOMS
100:             IF (FROZEN(J1)) THEN100:             IF (FROZEN(J1)) THEN
188:                GL(J1)=2.d0 * (GRAD1(J1)-GRAD(J1))/(ZETA*VECL**2)-2.0D0*DIAG4*LOCALV(J1)/VECL**2188:                GL(J1)=2.d0 * (GRAD1(J1)-GRAD(J1))/(ZETA*VECL**2)-2.0D0*DIAG4*LOCALV(J1)/VECL**2
189:             ELSEIF (NSECDIAG.EQ.2) THEN189:             ELSEIF (NSECDIAG.EQ.2) THEN
190:                ! second order central finite differences190:                ! second order central finite differences
191:                GL(J1)=(GRAD1(J1)-GRAD2(J1))/(ZETA*VECL**2)-2.0D0*DIAG2*LOCALV(J1)/VECL**2191:                GL(J1)=(GRAD1(J1)-GRAD2(J1))/(ZETA*VECL**2)-2.0D0*DIAG2*LOCALV(J1)/VECL**2
192:             ELSE192:             ELSE
193:                ! second order central finite differences193:                ! second order central finite differences
194:                GL(J1)=(GRAD1(J1)-GRAD2(J1))/(ZETA*VECL**2)-2.0D0*DIAG*LOCALV(J1)/VECL**2194:                GL(J1)=(GRAD1(J1)-GRAD2(J1))/(ZETA*VECL**2)-2.0D0*DIAG*LOCALV(J1)/VECL**2
195:             ENDIF195:             ENDIF
196: !           WRITE(*,'(A,I4,4G16.7)') 'secdiag> J1,GRAD1,GRAD2,LOCALV,GL=',J1,GRAD1(J1),GRAD2(J1),LOCALV(J1),GL(J1)196: !           WRITE(*,'(A,I4,4G16.7)') 'secdiag> J1,GRAD1,GRAD2,LOCALV,GL=',J1,GRAD1(J1),GRAD2(J1),LOCALV(J1),GL(J1)
197:          ENDDO197:          ENDDO
198:          IF (NFREEZE.LT.3.AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) CALL ORTHOGOPT(GL,COORDS,.FALSE.)198:          IF (NFREEZE.LT.3.AND..NOT.MIEFT) CALL ORTHOGOPT(GL,COORDS,.FALSE.)
199: !        PRINT *,'secdiag> before proj stuff GL:'199: !        PRINT *,'secdiag> before proj stuff GL:'
200: !        PRINT '(3F20.10)',GL(1:NOPT)200: !        PRINT '(3F20.10)',GL(1:NOPT)
201: !        CALL ORTHOGOPT(GL,COORDS,.FALSE.) ! seems to do some good for MSEVB201: !        CALL ORTHOGOPT(GL,COORDS,.FALSE.) ! seems to do some good for MSEVB
202: !202: !
203: !  Project out any component of the gradient along LOCALV (which is a unit vector)203: !  Project out any component of the gradient along LOCALV (which is a unit vector)
204: !  This is a big improvement for DFTB.  js850> The parallel components have already204: !  This is a big improvement for DFTB.  js850> The parallel components have already
205: !  been projected out, this simply improves the numerical precision205: !  been projected out, this simply improves the numerical precision
206: !206: !
207:          PROJ=0.0D0207:          PROJ=0.0D0
208:          DO J1=1,NOPT208:          DO J1=1,NOPT


r30776/xmylbfgs.f 2016-07-15 07:30:16.359230528 +0100 r30775/xmylbfgs.f 2016-07-15 07:30:20.287283382 +0100
 90:       ITER=0 90:       ITER=0
 91:       ITDONE=0 91:       ITDONE=0
 92:       NFAIL=0 92:       NFAIL=0
 93:       EVPERCENT=0.0D0 93:       EVPERCENT=0.0D0
 94:       FAILED=.FALSE. 94:       FAILED=.FALSE.
 95:       IF (FILTH.EQ.0) THEN 95:       IF (FILTH.EQ.0) THEN
 96:          WRITE(VECSTRING,'(A)') 'vector.dump.xlbfgs' 96:          WRITE(VECSTRING,'(A)') 'vector.dump.xlbfgs'
 97:       ELSE 97:       ELSE
 98:          WRITE(VECSTRING,'(A)') 'vector.dump.xlbfgs.'//TRIM(ADJUSTL(FILTHSTR)) 98:          WRITE(VECSTRING,'(A)') 'vector.dump.xlbfgs.'//TRIM(ADJUSTL(FILTHSTR))
 99:       ENDIF 99:       ENDIF
100:       IF (NFREEZE.LT.3.AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN100:       IF (NFREEZE.LT.3.AND..NOT.MIEFT) THEN
101:          CALL ORTHOGOPT(X,COORDS,.TRUE.)101:          CALL ORTHOGOPT(X,COORDS,.TRUE.)
102:       ELSE102:       ELSE
103:          CALL VECNORM(X,N)103:          CALL VECNORM(X,N)
104:       ENDIF104:       ENDIF
105:       VECSAVE(1:NOPT)=X(1:NOPT)105:       VECSAVE(1:NOPT)=X(1:NOPT)
106:       CALL SECDIAG(X,COORDS,EENERGY,ggrad,G,ENERGY,.TRUE.,RMS)106:       CALL SECDIAG(X,COORDS,EENERGY,ggrad,G,ENERGY,.TRUE.,RMS)
107:       CURVSAVE=ENERGY107:       CURVSAVE=ENERGY
108:       IF (PTEST) WRITE(*,'(A,2G20.10,A,I6,A)') 'xmylbfgs> Eigenvalue and RMS=',ENERGY,RMS,' after ',ITDONE,' steps'108:       IF (PTEST) WRITE(*,'(A,2G20.10,A,I6,A)') 'xmylbfgs> Eigenvalue and RMS=',ENERGY,RMS,' after ',ITDONE,' steps'
109: !     PRINT '(A)','X for this eigenvalue:'109: !     PRINT '(A)','X for this eigenvalue:'
110: !     PRINT '(3G20.10)',X(1:15)110: !     PRINT '(3G20.10)',X(1:15)
413:       DO I=1,N413:       DO I=1,N
414:          W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)414:          W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
415:          W(IYPT+NPT+I)= G(I)-W(I)415:          W(IYPT+NPT+I)= G(I)-W(I)
416:       ENDDO416:       ENDDO
417: C417: C
418: C  Does it help to keep X normalised? 418: C  Does it help to keep X normalised? 
419: C  Calling ORTHOGOPT will change X slightly and can therefore change419: C  Calling ORTHOGOPT will change X slightly and can therefore change
420: C  the eigenvalue. It could therefore cause NFAIL messages if the eigenvalue420: C  the eigenvalue. It could therefore cause NFAIL messages if the eigenvalue
421: C  isn;t allowed to rise slightly.421: C  isn;t allowed to rise slightly.
422: C422: C
423:       IF (NFREEZE.LT.3.AND..NOT.MIEFT .AND. (.NOT. NOTRANSROTT)) THEN423:       IF (NFREEZE.LT.3.AND..NOT.MIEFT) THEN
424:          CALL ORTHOGOPT(X,COORDS,.TRUE.)424:          CALL ORTHOGOPT(X,COORDS,.TRUE.)
425:       ELSE425:       ELSE
426:          CALL VECNORM(X,N)426:          CALL VECNORM(X,N)
427:       ENDIF427:       ENDIF
428: 428: 
429:       POINT=POINT+1429:       POINT=POINT+1
430:       IF (POINT.EQ.M) POINT=0430:       IF (POINT.EQ.M) POINT=0
431:       FIXIMAGE=.FALSE.431:       FIXIMAGE=.FALSE.
432:       IF ((FIXAFTER.GT.0).AND.(ITDONE.GE.FIXAFTER)) FIXIMAGE=.TRUE.432:       IF ((FIXAFTER.GT.0).AND.(ITDONE.GE.FIXAFTER)) FIXIMAGE=.TRUE.
433:       GOTO 10433:       GOTO 10


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0