hdiff output

r31072/commons.f90 2017-01-21 10:39:09.611814839 +0000 r31071/commons.f90 2017-01-21 10:39:11.031912530 +0000
 98:      &        BARRIERSHORT, FREEZE, RATESHORT, DUMMYRUNT, REWEIGHTT, REGROUPFREET, RFMULTIT, REGROUPFREEABT, READMINT, & 98:      &        BARRIERSHORT, FREEZE, RATESHORT, DUMMYRUNT, REWEIGHTT, REGROUPFREET, RFMULTIT, REGROUPFREEABT, READMINT, &
 99:      &        DUMPGROUPST, FREEPAIRT, KSHORTESTPATHST, KSHORT_FULL_PRINTT, DIJINITFLYT, BHINTERPT, ICINTERPT, & 99:      &        DUMPGROUPST, FREEPAIRT, KSHORTESTPATHST, KSHORT_FULL_PRINTT, DIJINITFLYT, BHINTERPT, ICINTERPT, &
100:      &        DUMMYTST, DOCKT, DSTAGE(6), USEPAIRST, LOWESTFRQT, BISECTT, NGTDISCONNECTALL, ANGLEAXIS2, TFOLDT, &100:      &        DUMMYTST, DOCKT, DSTAGE(6), USEPAIRST, LOWESTFRQT, BISECTT, NGTDISCONNECTALL, ANGLEAXIS2, TFOLDT, &
101:      &        SLURMT, INDEXCOSTFUNCTION, CVT, CVMINIMAT, DOST, IMFRQT, CLOSEFILEST, PULLT, FRICTIONT, ATOMMATCHFULL, &101:      &        SLURMT, INDEXCOSTFUNCTION, CVT, CVMINIMAT, DOST, IMFRQT, CLOSEFILEST, PULLT, FRICTIONT, ATOMMATCHFULL, &
102:      &        INTCONSTRAINTT, CHECKCONINT, INTLJT, INTERPCOSTFUNCTION, REMOVEUNCONNECTEDT, ATOMMATCHDIST, &102:      &        INTCONSTRAINTT, CHECKCONINT, INTLJT, INTERPCOSTFUNCTION, REMOVEUNCONNECTEDT, ATOMMATCHDIST, &
103:      &        DBPT, DBPTDT, DMBLPYT, EFIELDT, MSSTOCKT, NTIPT, PAHAT, PAPT, PATCHYDT, STOCKAAT, RBAAT, RBSYMT, TRAPT, SILANET, &103:      &        DBPT, DBPTDT, DMBLPYT, EFIELDT, MSSTOCKT, NTIPT, PAHAT, PAPT, PATCHYDT, STOCKAAT, RBAAT, RBSYMT, TRAPT, SILANET, &
104:      &        OHCELLT, INTFREEZET, LPERMDIST, PBST, RANDOMMETRICT, SSHT, ALLTST, USERPOTT, CHECKMINT, &104:      &        OHCELLT, INTFREEZET, LPERMDIST, PBST, RANDOMMETRICT, SSHT, ALLTST, USERPOTT, CHECKMINT, &
105:      &        CHECKTST, CHECKSPT, FROMLOWESTT, ADDMINXYZT, MACHINE, RATESCYCLET, NOINVERSION, NEWCONNECTIONST, NIMET, NIHEAM7T, &105:      &        CHECKTST, CHECKSPT, FROMLOWESTT, ADDMINXYZT, MACHINE, RATESCYCLET, NOINVERSION, NEWCONNECTIONST, NIMET, NIHEAM7T, &
106:      &        NIH2LEPST, DISTANCET, RATETARGETT, TARGETHIT, ALLOWABT, MICROTHERMT, RFKMCT, REGROUPKMCT, ONEREGROUPT, PHI4MODT, &106:      &        NIH2LEPST, DISTANCET, RATETARGETT, TARGETHIT, ALLOWABT, MICROTHERMT, RFKMCT, REGROUPKMCT, ONEREGROUPT, PHI4MODT, &
107:      &        PERSISTT, REGROUPPERSISTT, NOLABELST, SHANNONT, MAKEPAIRS, SKIPPAIRST, PERSISTAPPROXT, ALLCOMPONENTST, &107:      &        PERSISTT, REGROUPPERSISTT, NOLABELST, SHANNONT, MAKEPAIRS, SKIPPAIRST, PERSISTAPPROXT, ALLCOMPONENTST, &
108:      &        SHANNONRT, SHANNONZT, CUDAT, MLLJAT3, MLP3T, DIJPRUNET, PRINTSUMMARYT, MKTRAPT, MLPB3T, PRUNECYCLET, PAIRSIGNORET, &108:      &        SHANNONRT, SHANNONZT, CUDAT, MLLJAT3, MLP3T, DIJPRUNET, PRINTSUMMARYT, MKTRAPT, MLPB3T, PRUNECYCLET, PAIRSIGNORET
109:      &        NOTRANSROTT 
110: 109: 
111:       LOGICAL, ALLOCATABLE :: SHIFTABLE(:)110:       LOGICAL, ALLOCATABLE :: SHIFTABLE(:)
112:       CHARACTER(LEN=80) COORDSLIGANDSTR, COORDSCOMPLEXSTR, COORDSPROTEINSTR111:       CHARACTER(LEN=80) COORDSLIGANDSTR, COORDSCOMPLEXSTR, COORDSPROTEINSTR
113:       CHARACTER(LEN=80) EXEC,EXECGMIN112:       CHARACTER(LEN=80) EXEC,EXECGMIN
114:       CHARACTER(LEN=80) PATHNAME, MINNAME, ADDMINXYZNAME, ALLCOMPS113:       CHARACTER(LEN=80) PATHNAME, MINNAME, ADDMINXYZNAME, ALLCOMPS
115:       CHARACTER(LEN=150) COPYFILES114:       CHARACTER(LEN=150) COPYFILES
116:       CHARACTER(LEN=80) USEPAIRSFILE115:       CHARACTER(LEN=80) USEPAIRSFILE
117:       CHARACTER(LEN=80) MAKEPAIRSFILE116:       CHARACTER(LEN=80) MAKEPAIRSFILE
118:       CHARACTER(LEN=2) DIRECTION117:       CHARACTER(LEN=2) DIRECTION
119:       CHARACTER(LEN=5) UNCONNECTEDS118:       CHARACTER(LEN=5) UNCONNECTEDS


r31072/inertia.f 2017-01-21 10:39:09.831830017 +0000 r31071/inertia.f 2017-01-21 10:39:11.383936636 +0000
 86:       DO J1=1,NOPT/3 86:       DO J1=1,NOPT/3
 87:          CMX=CMX+Q(3*(J1-1)+1)*MASS(J1) 87:          CMX=CMX+Q(3*(J1-1)+1)*MASS(J1)
 88:          CMY=CMY+Q(3*(J1-1)+2)*MASS(J1) 88:          CMY=CMY+Q(3*(J1-1)+2)*MASS(J1)
 89:          CMZ=CMZ+Q(3*(J1-1)+3)*MASS(J1) 89:          CMZ=CMZ+Q(3*(J1-1)+3)*MASS(J1)
 90:          MASST=MASST+MASS(J1) 90:          MASST=MASST+MASS(J1)
 91: C        WRITE(*,'(A,I5,4F20.10)') 'I,ATMASS(I),=',J1,MASS(J1),Q(3*(J1-1)+1),Q(3*(J1-1)+2),Q(3*(J1-1)+3) 91: C        WRITE(*,'(A,I5,4F20.10)') 'I,ATMASS(I),=',J1,MASS(J1),Q(3*(J1-1)+1),Q(3*(J1-1)+2),Q(3*(J1-1)+3)
 92:       ENDDO 92:       ENDDO
 93:       CMX=CMX/MASST 93:       CMX=CMX/MASST
 94:       CMY=CMY/MASST 94:       CMY=CMY/MASST
 95:       CMZ=CMZ/MASST 95:       CMZ=CMZ/MASST
 96:       IF (NOTRANSROTT) THEN 
 97:          CMX=0.0D0 
 98:          CMY=0.0D0 
 99:          CMZ=0.0D0 
100:       ENDIF 
101: C     PRINT*,'MASST,CMX,CMY,CMZ=',MASST,CMX,CMY,CMZ 96: C     PRINT*,'MASST,CMX,CMY,CMZ=',MASST,CMX,CMY,CMZ
102:       DO J1=1,NOPT/3 97:       DO J1=1,NOPT/3
103:          Q(3*(J1-1)+1)=Q(3*(J1-1)+1)-CMX 98:          Q(3*(J1-1)+1)=Q(3*(J1-1)+1)-CMX
104:          Q(3*(J1-1)+2)=Q(3*(J1-1)+2)-CMY 99:          Q(3*(J1-1)+2)=Q(3*(J1-1)+2)-CMY
105:          Q(3*(J1-1)+3)=Q(3*(J1-1)+3)-CMZ100:          Q(3*(J1-1)+3)=Q(3*(J1-1)+3)-CMZ
106:       ENDDO101:       ENDDO
107: 102: 
108:       DO J1=1,3103:       DO J1=1,3
109:          DO J2=1,3104:          DO J2=1,3
110:             IT(J1,J2)=0.0D0105:             IT(J1,J2)=0.0D0


r31072/keywords.f 2017-01-21 10:39:10.059845728 +0000 r31071/keywords.f 2017-01-21 10:39:11.607951976 +0000
391:       MAKEPAIRS = .FALSE.391:       MAKEPAIRS = .FALSE.
392:       SKIPPAIRST = .FALSE.392:       SKIPPAIRST = .FALSE.
393: 393: 
394: ! kr366: connect unconnected minima394: ! kr366: connect unconnected minima
395:       CONUNCONT=.FALSE.395:       CONUNCONT=.FALSE.
396:       CONNECTLOWESTT=.FALSE.396:       CONNECTLOWESTT=.FALSE.
397:       CONNECTETHRESHT=.FALSE.397:       CONNECTETHRESHT=.FALSE.
398:       CONNECTDTHRESHT=.FALSE.398:       CONNECTDTHRESHT=.FALSE.
399:       REFMIN=0399:       REFMIN=0
400:       NATT=2400:       NATT=2
401: ! 
402: ! Forbid overall translation and rotation in distance/alignment routines. 
403: ! 
404:       NOTRANSROTT=.FALSE. 
405: 401: 
406: C402: C
407: C SAT: Now, there it is! Unit 5 is a standard input unit, and we must not use it!403: C SAT: Now, there it is! Unit 5 is a standard input unit, and we must not use it!
408: C 5--->959404: C 5--->959
409: C405: C
410:       OPEN (959,FILE='pathdata',STATUS='OLD')406:       OPEN (959,FILE='pathdata',STATUS='OLD')
411: C407: C
412: C  Read from unit 959408: C  Read from unit 959
413: C409: C
414:       IR=959410:       IR=959
1637: C  in minpermdist. Needed for some twod examples.1633: C  in minpermdist. Needed for some twod examples.
1638: C1634: C
1639:       ELSE IF (WORD.EQ.'NOINVERSION') THEN1635:       ELSE IF (WORD.EQ.'NOINVERSION') THEN
1640:          NOINVERSION=.TRUE.1636:          NOINVERSION=.TRUE.
1641: C1637: C
1642: C  If NOPOINTS is specified then setup should not try to read the min.points or ts.points1638: C  If NOPOINTS is specified then setup should not try to read the min.points or ts.points
1643: C  files. Should be the default for post-DPS database kinetics analysis such as KMC.1639: C  files. Should be the default for post-DPS database kinetics analysis such as KMC.
1644: C1640: C
1645:       ELSE IF (WORD.EQ.'NOPOINTS') THEN1641:       ELSE IF (WORD.EQ.'NOPOINTS') THEN
1646:          NOPOINTS=.TRUE.1642:          NOPOINTS=.TRUE.
1647: ! 
1648: ! Forbid overall translation and rotation in distance/alignment calculations. 
1649: ! 
1650:       ELSE IF (WORD.EQ.'NOTRANSROT') THEN 
1651:             NOTRANSROTT=.TRUE. 
1652: 1643: 
1653:       ELSE IF (WORD.EQ.'NTIP') THEN1644:       ELSE IF (WORD.EQ.'NTIP') THEN
1654: 1645: 
1655:          CALL READI(TIPID)1646:          CALL READI(TIPID)
1656: 1647: 
1657:          IF (TIPID == 4) THEN1648:          IF (TIPID == 4) THEN
1658:             NRBSITES = 41649:             NRBSITES = 4
1659:          ELSE1650:          ELSE
1660:             PRINT *, 'TIPID NOT EQUAL TO 4 NOT YET DEFINED'1651:             PRINT *, 'TIPID NOT EQUAL TO 4 NOT YET DEFINED'
1661:             STOP1652:             STOP


r31072/mindist.f 2017-01-21 10:39:10.291861676 +0000 r31071/mindist.f 2017-01-21 10:39:11.855968959 +0000
 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 COMMONS,ONLY : FREEZE, PULLT, GTHOMSONT, GTHOMMET, NOTRANSROTT ! hk286 29:       USE COMMONS,ONLY : FREEZE, PULLT, GTHOMSONT, GTHOMMET ! 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                 ROTMAT(3,3),  34:      1                 ROTMAT(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
 95:             R(3,J1)=RA(3*(J1-1)+3) 95:             R(3,J1)=RA(3*(J1-1)+3)
 96:             R0(1,J1)=RB(3*(J1-1)+1) 96:             R0(1,J1)=RB(3*(J1-1)+1)
 97:             R0(2,J1)=RB(3*(J1-1)+2) 97:             R0(2,J1)=RB(3*(J1-1)+2)
 98:             R0(3,J1)=RB(3*(J1-1)+3) 98:             R0(3,J1)=RB(3*(J1-1)+3)
 99:          ENDDO 99:          ENDDO
100:       ENDIF100:       ENDIF
101: C101: C
102: C     put c.o.m. to origin102: C     put c.o.m. to origin
103: C103: C
104: ! hk286104: ! hk286
105:       IF ((.NOT.FREEZE).AND.(.NOT.GTHOMSONT) .AND. (.NOT. NOTRANSROTT)) THEN105:       IF ((.NOT.FREEZE).AND.(.NOT.GTHOMSONT)) THEN
106:          do ig=1,3106:          do ig=1,3
107:             rg=0.d0107:             rg=0.d0
108:             rg0=0.d0108:             rg0=0.d0
109:             do i=1,nsize109:             do i=1,nsize
110:                rg=rg+R(ig,i)110:                rg=rg+R(ig,i)
111:                rg0=rg0+R0(ig,i)111:                rg0=rg0+R0(ig,i)
112:             enddo112:             enddo
113:             rg=rg/nsize113:             rg=rg/nsize
114:             rg0=rg0/nsize114:             rg0=rg0/nsize
115:             do i=1,nsize115:             do i=1,nsize
151:             PRINT*,'convergence failure in mind'151:             PRINT*,'convergence failure in mind'
152: c           STOP152: c           STOP
153:          ELSE153:          ELSE
154: C           WRITE(*,'(A,2F15.5,A,I6,A,I6)') 'Initial and final distances:',DIST0,DIST,' iterations=',ITER,' NCOUNT=',NCOUNT154: C           WRITE(*,'(A,2F15.5,A,I6,A,I6)') 'Initial and final distances:',DIST0,DIST,' iterations=',ITER,' NCOUNT=',NCOUNT
155:             DO J1=1,NSIZE155:             DO J1=1,NSIZE
156:                R0(1,J1)=R1(1,J1)156:                R0(1,J1)=R1(1,J1)
157:                R0(2,J1)=R1(2,J1)157:                R0(2,J1)=R1(2,J1)
158:                R0(3,J1)=R1(3,J1)158:                R0(3,J1)=R1(3,J1)
159:             ENDDO159:             ENDDO
160:             ! hk286160:             ! hk286
161:             IF (.NOT.(NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5)))) P(1)=2*(DPRAND()-0.5D0)/10.0D0161:             IF (.NOT.(TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5)))) P(1)=2*(DPRAND()-0.5D0)/10.0D0
162:             IF (.NOT.(NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5)))) P(2)=2*(DPRAND()-0.5D0)/10.0D0162:             IF (.NOT.(TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5)))) P(2)=2*(DPRAND()-0.5D0)/10.0D0
163:             P(3)=2*(DPRAND()-0.5D0)/10.0D0163:             P(3)=2*(DPRAND()-0.5D0)/10.0D0
164:             GOTO 10164:             GOTO 10
165:          ENDIF165:          ENDIF
166:       ENDIF166:       ENDIF
167: C167: C
168: C  This block allows the second geometry to rotate out of the xy plane; not allowed!168: C  This block allows the second geometry to rotate out of the xy plane; not allowed!
169: C169: C
170: C     IF (TWOD) THEN170: C     IF (TWOD) THEN
171: C        IF (AGAIN) THEN171: C        IF (AGAIN) THEN
172: C           WRITE(*,'(A,2F20.10)') 'DIST,DSAVE=',DIST,DSAVE172: C           WRITE(*,'(A,2F20.10)') 'DIST,DSAVE=',DIST,DSAVE
252:             RA(3*(J1-1)+3)=R(3,J1)252:             RA(3*(J1-1)+3)=R(3,J1)
253:          ENDDO253:          ENDDO
254:       ENDIF254:       ENDIF
255: 255: 
256:       DO J1=1,3256:       DO J1=1,3
257:          DO J2=1,3257:          DO J2=1,3
258:             ROTMAT(J2,J1)=OMEGATOT(J2,J1)258:             ROTMAT(J2,J1)=OMEGATOT(J2,J1)
259:          ENDDO259:          ENDDO
260:       ENDDO260:       ENDDO
261: 261: 
262:       IF (PRESERVET .OR. NOTRANSROTT) RB(1:3*NATOMS)=RBSAVE(1:3*NATOMS)262:       IF (PRESERVET) RB(1:3*NATOMS)=RBSAVE(1:3*NATOMS)
263: 263: 
264:       DEALLOCATE(R, R0, R1, R1SAVE, RBSAVE)264:       DEALLOCATE(R, R0, R1, R1SAVE, RBSAVE)
265: 265: 
266:       RETURN266:       RETURN
267:       END267:       END
268: 268: 
269: c__________________________________________________________________________269: c__________________________________________________________________________
270: 270: 
271:         subroutine inverse(A,B)271:         subroutine inverse(A,B)
272:         IMPLICIT NONE272:         IMPLICIT NONE
694:       RETURN694:       RETURN
695:       END695:       END
696: 696: 
697: C        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION697: C        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
698: C                          JORGE NOCEDAL698: C                          JORGE NOCEDAL
699: C                        *** July 1990 ***699: C                        *** July 1990 ***
700: C700: C
701: C        Line search removed plus small modifications, DJW 2001701: C        Line search removed plus small modifications, DJW 2001
702: C702: C
703:       SUBROUTINE MMYLBFGS(X,EPS,MFLAG,ENERGY,ITMAX,ITDONE,R,R0,R1)703:       SUBROUTINE MMYLBFGS(X,EPS,MFLAG,ENERGY,ITMAX,ITDONE,R,R0,R1)
704:       USE COMMONS, ONLY : PULLT, TWOD, GTHOMSONT, GTHOMMET, NOTRANSROTT ! hk286704:       USE COMMONS, ONLY : PULLT, TWOD, GTHOMSONT, GTHOMMET ! hk286
705:       use porfuncs705:       use porfuncs
706:       IMPLICIT NONE706:       IMPLICIT NONE
707:       INTEGER N,M,J1,J2,ITMAX,ITDONE,NFAIL,J707:       INTEGER N,M,J1,J2,ITMAX,ITDONE,NFAIL,J
708:       PARAMETER (N=3,M=5)  !  MMUPDATE is actually ignored708:       PARAMETER (N=3,M=5)  !  MMUPDATE is actually ignored
709:       DOUBLE PRECISION X(N),G(3),DIAG(N),W(N*(2*M+1)+2*M),SLENGTH,DDOT,OVERLAP,DISTFUNC709:       DOUBLE PRECISION X(N),G(3),DIAG(N),W(N*(2*M+1)+2*M),SLENGTH,DDOT,OVERLAP,DISTFUNC
710:       DOUBLE PRECISION EPS,DUMMY1,ENERGY,ENEW,RMS,ALPHA,GSAVE(3)710:       DOUBLE PRECISION EPS,DUMMY1,ENERGY,ENEW,RMS,ALPHA,GSAVE(3)
711:       DOUBLE PRECISION GNORM,STP,YS,YY,SQ,YR,BETA,ROTMAT(3,3),OMEGATOT(3,3),OTEMP(3,3)711:       DOUBLE PRECISION GNORM,STP,YS,YY,SQ,YR,BETA,ROTMAT(3,3),OMEGATOT(3,3),OTEMP(3,3)
712:       INTEGER ITER,POINT,ISPT,IYPT,BOUND,NPT,CP,I,INMC,IYCN,ISCN,NDECREASE712:       INTEGER ITER,POINT,ISPT,IYPT,BOUND,NPT,CP,I,INMC,IYCN,ISCN,NDECREASE
713:       LOGICAL MFLAG, PTEST713:       LOGICAL MFLAG, PTEST
714:       INTEGER NSIZE, ISTAT714:       INTEGER NSIZE, ISTAT
722: 722: 
723:       PTEST=.TRUE.723:       PTEST=.TRUE.
724:       PTEST=.FALSE.724:       PTEST=.FALSE.
725:       ALPHA=1.0D0725:       ALPHA=1.0D0
726:       NFAIL=0726:       NFAIL=0
727:       ITER=0727:       ITER=0
728:       ITDONE=0728:       ITDONE=0
729:       ENERGY=DISTFUNC(X,R,R0,R1)729:       ENERGY=DISTFUNC(X,R,R0,R1)
730:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)730:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)
731:       ! hk286731:       ! hk286
732:       IF (NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0732:       IF (TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0
733:       IF (NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0733:       IF (TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0
734:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)734:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)
735:       DO J1=1,N735:       DO J1=1,N
736:          G(J1)=GSAVE(J1)736:          G(J1)=GSAVE(J1)
737:       ENDDO737:       ENDDO
738: 738: 
739:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,' LBFGS steps'739:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,' LBFGS steps'
740: 740: 
741: 10    CALL FLUSH(6)741: 10    CALL FLUSH(6)
742:       MFLAG=.FALSE.742:       MFLAG=.FALSE.
743:       IF (RMS.LE.EPS) THEN743:       IF (RMS.LE.EPS) THEN
853:       ENDIF853:       ENDIF
854: C854: C
855: C  Store the new search direction855: C  Store the new search direction
856: C856: C
857:       IF (ITER.GT.0) THEN857:       IF (ITER.GT.0) THEN
858:          DO I=1,N858:          DO I=1,N
859:             W(ISPT+POINT*N+I)= W(I)859:             W(ISPT+POINT*N+I)= W(I)
860:          ENDDO860:          ENDDO
861:       ENDIF861:       ENDIF
862:       ! hk286862:       ! hk286
863:       IF (NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+1)=0.0D0863:       IF (TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+1)=0.0D0
864:       IF (NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+2)=0.0D0864:       IF (TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) W(ISPT+POINT*N+2)=0.0D0
865: 865: 
866: C     OVERLAP=DDOT(N,G,1,W,1)/SQRT(DDOT(N,G,1,G,1)*DDOT(N,W,1,W,1))866: C     OVERLAP=DDOT(N,G,1,W,1)/SQRT(DDOT(N,G,1,G,1)*DDOT(N,W,1,W,1))
867:       DOT1=SQRT(DDOT(N,G,1,G,1))867:       DOT1=SQRT(DDOT(N,G,1,G,1))
868:       DOT2=SQRT(DDOT(N,W,1,W,1))868:       DOT2=SQRT(DDOT(N,W,1,W,1))
869:       OVERLAP=0.0D0869:       OVERLAP=0.0D0
870:       IF (DOT1*DOT2.NE.0.0D0) OVERLAP=DDOT(N,G,1,W,1)/(DOT1*DOT2)870:       IF (DOT1*DOT2.NE.0.0D0) OVERLAP=DDOT(N,G,1,W,1)/(DOT1*DOT2)
871: 871: 
872: C     PRINT*,'OVERLAP,DIAG(1)=',OVERLAP,DIAG(1)872: C     PRINT*,'OVERLAP,DIAG(1)=',OVERLAP,DIAG(1)
873: C     PRINT*,'G . G=',DDOT(N,G,1,G,1)873: C     PRINT*,'G . G=',DDOT(N,G,1,G,1)
874: C     PRINT*,'W . W=',DDOT(N,W,1,W,1)874: C     PRINT*,'W . W=',DDOT(N,W,1,W,1)
954:       X(2)=0.0D0954:       X(2)=0.0D0
955:       X(3)=0.0D0955:       X(3)=0.0D0
956:       DO I=1,NSIZE956:       DO I=1,NSIZE
957:          DO J=1,3957:          DO J=1,3
958:             R0(J,I)=R1(J,I)958:             R0(J,I)=R1(J,I)
959:          ENDDO959:          ENDDO
960:       ENDDO960:       ENDDO
961: 961: 
962:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)962:       CALL DDISTFUNC(X,GSAVE,R,R0,R1)
963:       ! hk286963:       ! hk286
964:       IF (NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0964:       IF (TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(1)=0.0D0
965:       IF (NOTRANSROTT.OR.TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0965:       IF (TWOD.OR.PULLT.OR.(GTHOMSONT.AND.(GTHOMMET<5))) GSAVE(2)=0.0D0
966:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)966:       RMS=SQRT((GSAVE(1)**2+GSAVE(2)**2+GSAVE(3)**2)/3)
967:       DO J1=1,N967:       DO J1=1,N
968:          G(J1)=GSAVE(J1)968:          G(J1)=GSAVE(J1)
969:       ENDDO969:       ENDDO
970:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A,G15.5)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,970:       IF (PTEST) WRITE(*,'(A,2F20.10,A,I6,A,G15.5)') ' Distance and RMS force=',ENERGY,RMS,' after ',ITDONE,
971:      1        ' LBFGS steps, step:',STP*SLENGTH971:      1        ' LBFGS steps, step:',STP*SLENGTH
972: C972: C
973: C     Compute the new step and gradient change973: C     Compute the new step and gradient change
974: C974: C
975: 30    NPT=POINT*N975: 30    NPT=POINT*N


r31072/minpermdist.f90 2017-01-21 10:39:10.519877347 +0000 r31071/minpermdist.f90 2017-01-21 10:39:12.576018459 +0000
 49: !  +/- RMATCUMUL ROTA (COORDSA - CMA) = permutation(DUMMYA) 49: !  +/- RMATCUMUL ROTA (COORDSA - CMA) = permutation(DUMMYA)
 50: !  where +/- is given by the value of INVERT. 50: !  where +/- is given by the value of INVERT.
 51: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 51: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the
 52: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we 52: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we
 53: !  are doing an ion trap potential. 53: !  are doing an ion trap potential.
 54: ! 54: !
 55: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST,USEINT) 55: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST,USEINT)
 56: USE COMMONS,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, AMBERT, NFREEZE, CHARMMT, RBAAT, PULLT, & 56: USE COMMONS,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, AMBERT, NFREEZE, CHARMMT, RBAAT, PULLT, &
 57:   &               ANGLEAXIS, PERMISOMER, PERMDIST, ZSYM, INTCONSTRAINTT, INTLJT, OHCELLT, ATOMMATCHDIST, LPERMDIST, & 57:   &               ANGLEAXIS, PERMISOMER, PERMDIST, ZSYM, INTCONSTRAINTT, INTLJT, OHCELLT, ATOMMATCHDIST, LPERMDIST, &
 58:   &               TRAPT, NRANROT, MACROCYCLET, MCYCLEPERIOD, MCYCLEREPEATS, NOINVERSION, GTHOMSONT, & 58:   &               TRAPT, NRANROT, MACROCYCLET, MCYCLEPERIOD, MCYCLEREPEATS, NOINVERSION, GTHOMSONT, &
 59:   &               RATETARGETT, MLP3T, NOPT, MKTRAPT, MLPB3T, NOTRANSROTT 59:   &               RATETARGETT, MLP3T, NOPT, MKTRAPT, MLPB3T
 60: USE PORFUNCS  60: USE PORFUNCS 
 61: USE UTILS,ONLY : GETUNIT 61: USE UTILS,ONLY : GETUNIT
 62:  62: 
 63: IMPLICIT NONE 63: IMPLICIT NONE
 64:  64: 
 65: INTEGER :: MAXIMUMTRIES=10 65: INTEGER :: MAXIMUMTRIES=10
 66: INTEGER NATOMS, NPERM, PATOMS, NTRIES, BESTINVERT, I, LUNIT 66: INTEGER NATOMS, NPERM, PATOMS, NTRIES, BESTINVERT, I, LUNIT
 67: INTEGER INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1, OPNUM, J3, NROTDONE 67: INTEGER INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1, OPNUM, J3, NROTDONE
 68: INTEGER NORBITB1, NORBITB2, NCHOOSEB1, NCHOOSEB2 68: INTEGER NORBITB1, NORBITB2, NCHOOSEB1, NCHOOSEB2
 69:  69: 
166: !  WRITE(LUNIT,'(I6)') NATOMS/2166: !  WRITE(LUNIT,'(I6)') NATOMS/2
167: !  WRITE(LUNIT,'(A)') 'B initial'167: !  WRITE(LUNIT,'(A)') 'B initial'
168: !  DO J3=1,NATOMS/2168: !  DO J3=1,NATOMS/2
169: !      WRITE(LUNIT,'(A2,2X,3F20.10)') 'LA',COORDSB(3*(J3-1)+1),COORDSB(3*(J3-1)+2),COORDSB(3*(J3-1)+3)169: !      WRITE(LUNIT,'(A2,2X,3F20.10)') 'LA',COORDSB(3*(J3-1)+1),COORDSB(3*(J3-1)+2),COORDSB(3*(J3-1)+3)
170: !  ENDDO170: !  ENDDO
171: !  CLOSE(LUNIT)171: !  CLOSE(LUNIT)
172: !172: !
173: !  Calculate original centres of mass.173: !  Calculate original centres of mass.
174: !174: !
175: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0175: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0
176: IF ((NFREEZE.LE.0).AND.(.NOT.MKTRAPT).AND.(.NOT. NOTRANSROTT)) THEN176: IF ((NFREEZE.LE.0).AND.(.NOT.MKTRAPT)) THEN
177:    IF (RBAAT) THEN177:    IF (RBAAT) THEN
178:       DO J1=1,NATOMS/2178:       DO J1=1,NATOMS/2
179:          CMAX=CMAX+COORDSA(3*(J1-1)+1)179:          CMAX=CMAX+COORDSA(3*(J1-1)+1)
180:          CMAY=CMAY+COORDSA(3*(J1-1)+2)180:          CMAY=CMAY+COORDSA(3*(J1-1)+2)
181:          CMAZ=CMAZ+COORDSA(3*(J1-1)+3)181:          CMAZ=CMAZ+COORDSA(3*(J1-1)+3)
182:       ENDDO182:       ENDDO
183:       CMAX=2*CMAX/NATOMS; CMAY=2*CMAY/NATOMS; CMAZ=2*CMAZ/NATOMS183:       CMAX=2*CMAX/NATOMS; CMAY=2*CMAY/NATOMS; CMAZ=2*CMAZ/NATOMS
184:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0184:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0
185:       DO J1=1,NATOMS/2185:       DO J1=1,NATOMS/2
186:          CMBX=CMBX+COORDSB(3*(J1-1)+1)186:          CMBX=CMBX+COORDSB(3*(J1-1)+1)
207: !207: !
208: ! It is possible for the standard orientation to result in a distance that is worse than208: ! It is possible for the standard orientation to result in a distance that is worse than
209: ! the starting distance. Hence we need to set XBEST here.209: ! the starting distance. Hence we need to set XBEST here.
210: !210: !
211: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)211: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)
212: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)212: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)
213: DBEST=1.0D100213: DBEST=1.0D100
214: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)214: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
215: DBEST=DISTANCE**2215: DBEST=DISTANCE**2
216: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE216: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE
217: 217: IF ((NFREEZE.GT.0).OR.BULKT.OR.MKTRAPT) THEN
218: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then 
219: ! translate the origin of structure A to the CoM of structure B. 
220: IF ((NFREEZE.GT.0).OR.BULKT.OR.MKTRAPT.OR.NOTRANSROTT) THEN 
221:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)218:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
222: ELSE219: ELSE
223:    DO J1=1,NATOMS220:    DO J1=1,NATOMS
224:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX221:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX
225:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY222:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY
226:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ223:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ
227:    ENDDO224:    ENDDO
228: ENDIF225: ENDIF
229: BESTINVERT=1226: BESTINVERT=1
230: DO J1=1,NATOMS227: DO J1=1,NATOMS
281: ! Problems can occur if we don't use all the atoms specified by NORBIT1 and NORBIT2278: ! Problems can occur if we don't use all the atoms specified by NORBIT1 and NORBIT2
282: ! because of the numerical cutoffs employed in MYORIENT. We could miss the279: ! because of the numerical cutoffs employed in MYORIENT. We could miss the
283: ! right orientation! 280: ! right orientation! 
284: !281: !
285: ! If we use MYORIENT to produce particular orientations then we end up aligning 282: ! If we use MYORIENT to produce particular orientations then we end up aligning 
286: ! COORDSA not with COORDSB but with the standard orientation of COORDSB in DUMMYB.283: ! COORDSA not with COORDSB but with the standard orientation of COORDSB in DUMMYB.
287: ! We now deal with this by tracking the complete transformation, including the284: ! We now deal with this by tracking the complete transformation, including the
288: ! contribution of MYORIENT using ROTB and ROTINVB.285: ! contribution of MYORIENT using ROTB and ROTINVB.
289: !286: !
290: DISTANCE=0.0D0 287: DISTANCE=0.0D0 
291: IF (NFREEZE.LE.0..AND. .NOT. NOTRANSROTT) THEN288: IF (NFREEZE.LE.0) THEN
292:    IF (BULKT) THEN289:    IF (BULKT) THEN
293:       IF (BMTEST) BMCOORDSSV(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)290:       IF (BMTEST) BMCOORDSSV(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)
294:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;291:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;
295:       CALL BULKMINDIST(DUMMYB,DUMMYA,BMCOORDS,NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,.TRUE.,TNMATCH, BMTEST)292:       CALL BULKMINDIST(DUMMYB,DUMMYA,BMCOORDS,NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,.TRUE.,TNMATCH, BMTEST)
296:       IF (PITEST) THEN293:       IF (PITEST) THEN
297:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)294:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
298:          RMATBEST(1:3,1:3)=0.0D0295:          RMATBEST(1:3,1:3)=0.0D0
299:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0296:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
300:          DISTANCE=SQRT(DISTANCE)297:          DISTANCE=SQRT(DISTANCE)
301:          RETURN298:          RETURN
616: IF ((NCHOOSE2.EQ.NORBIT2).AND.(NCHOOSE1.EQ.NORBIT1).AND.(INVERT.EQ.1)) THEN613: IF ((NCHOOSE2.EQ.NORBIT2).AND.(NCHOOSE1.EQ.NORBIT1).AND.(INVERT.EQ.1)) THEN
617: !614: !
618: ! don't try inversion for bulk or charmm or amber or frozen atoms615: ! don't try inversion for bulk or charmm or amber or frozen atoms
619: ! Inversion is allowed for all macrocycles in AMBER/CHARMM616: ! Inversion is allowed for all macrocycles in AMBER/CHARMM
620: !617: !
621: !  IF ((BULKT.OR.CHARMMT.OR.AMBERT.OR.(NFREEZE.GT.0)).AND.MACROCYCLET.EQV..FALSE.) GOTO 50 618: !  IF ((BULKT.OR.CHARMMT.OR.AMBERT.OR.(NFREEZE.GT.0)).AND.MACROCYCLET.EQV..FALSE.) GOTO 50 
622: !619: !
623: !  Bug fix DJW 12/7/12. When MACROCYCLET was false all inversion tests were turned off!!620: !  Bug fix DJW 12/7/12. When MACROCYCLET was false all inversion tests were turned off!!
624: !621: !
625:    IF (NOINVERSION.OR.BULKT.OR.(CHARMMT.AND.(.NOT.MACROCYCLET)).OR.(AMBERT.AND.(.NOT.MACROCYCLET)) &622:    IF (NOINVERSION.OR.BULKT.OR.(CHARMMT.AND.(.NOT.MACROCYCLET)).OR.(AMBERT.AND.(.NOT.MACROCYCLET)) &
626:   &     .OR.(NFREEZE.GT.0).OR.NOTRANSROTT) GOTO 50 623:   &     .OR.(NFREEZE.GT.0)) GOTO 50 
627: !  IF (DEBUG) PRINT '(A)','minpermdist> inverting geometry for comparison with target'624: !  IF (DEBUG) PRINT '(A)','minpermdist> inverting geometry for comparison with target'
628:    INVERT=-1625:    INVERT=-1
629:    GOTO 60626:    GOTO 60
630: ENDIF627: ENDIF
631: MCYCLESTEP=MCYCLESTEP+1 628: MCYCLESTEP=MCYCLESTEP+1 
632: IF (MACROCYCLET.AND.(MCYCLESTEP.LE.MCYCLEREPEATS)) GOTO 70 !End of macrocycle loop629: IF (MACROCYCLET.AND.(MCYCLESTEP.LE.MCYCLEREPEATS)) GOTO 70 !End of macrocycle loop
633: IF (NROTDONE.LT.NRANROT) GOTO 11630: IF (NROTDONE.LT.NRANROT) GOTO 11
634: 631: 
635: 50 DISTANCE=DBEST632: 50 DISTANCE=DBEST
636: !633: !
637: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.634: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.
638: !  Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, 635: !  Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, 
639: !  undoing the reorientation to DUMMYB from MYORIENT. 636: !  undoing the reorientation to DUMMYB from MYORIENT. 
640: !  We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA)637: !  We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA)
641: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 638: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 
642: !  (aside from a possible permutation of the atom ordering)639: !  (aside from a possible permutation of the atom ordering)
643: !  640: !  
644:    IF ((NFREEZE.GT.0).OR.RBAAT.OR.MKTRAPT .OR. NOTRANSROTT) THEN641:    IF ((NFREEZE.GT.0).OR.RBAAT.OR.MKTRAPT) THEN
645:       XDUMMY=0.0D0642:       XDUMMY=0.0D0
646:       DO J1=1,NATOMS643:       DO J1=1,NATOMS
647:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &644:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
648:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &645:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
649:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2646:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
650:       ENDDO647:       ENDDO
651:    ELSEIF (BULKT) THEN648:    ELSEIF (BULKT) THEN
652:       XDUMMY=0.0D0649:       XDUMMY=0.0D0
653:       DO J1=1,NATOMS650:       DO J1=1,NATOMS
654:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1) - BOXLX*NINT((COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))/BOXLX))**2+ &651:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1) - BOXLX*NINT((COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))/BOXLX))**2+ &
671:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL .AND. (.NOT. RBAAT)) THEN668:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL .AND. (.NOT. RBAAT)) THEN
672:       PRINT '(2(A,G20.10))','minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &669:       PRINT '(2(A,G20.10))','minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &
673:   &                         ' should be ',SQRT(DISTANCE)670:   &                         ' should be ',SQRT(DISTANCE)
674:       PRINT '(A)','transformed XBEST:'671:       PRINT '(A)','transformed XBEST:'
675:       PRINT '(3F20.10)',XBEST(1:3*NATOMS)672:       PRINT '(3F20.10)',XBEST(1:3*NATOMS)
676:       PRINT '(A)','COORDSB:'673:       PRINT '(A)','COORDSB:'
677:       PRINT '(3F20.10)',COORDSB(1:3*NATOMS)674:       PRINT '(3F20.10)',COORDSB(1:3*NATOMS)
678:       STOP675:       STOP
679:    ENDIF676:    ENDIF
680: 677: 
681:    IF ((NFREEZE.GT.0).OR.MKTRAPT.OR.NOTRANSROTT) THEN678:    IF ((NFREEZE.GT.0).OR.MKTRAPT) THEN
682: !no rotation for NFREEZE .gt. 0679: !no rotation for NFREEZE .gt. 0
683:       RMATBEST(1:3,1:3)=0.0D0680:       RMATBEST(1:3,1:3)=0.0D0
684:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0681:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
685:    ELSE682:    ELSE
686:       RMATBEST=MATMUL(ROTINVBBEST,RMATBEST)683:       RMATBEST=MATMUL(ROTINVBBEST,RMATBEST)
687:    ENDIF684:    ENDIF
688: !685: !
689: ! For TRAP potentials we need to preserve the centre of mass because686: ! For TRAP potentials we need to preserve the centre of mass because
690: ! the origin is a fixed point. Try using RMAT best and rotating COORDSA687: ! the origin is a fixed point. Try using RMAT best and rotating COORDSA
691: ! around the origin.688: ! around the origin.


r31072/newmindist.f90 2017-01-21 10:39:10.767894392 +0000 r31071/newmindist.f90 2017-01-21 10:39:12.812034223 +0000
 21: !   closest geometry to RA if PRESERVET is .FALSE. 21: !   closest geometry to RA if PRESERVET is .FALSE.
 22: ! 22: !
 23: !   New analytic method based on quaterions from 23: !   New analytic method based on quaterions from
 24: !   Kearsley, Acta Cryst. A, 45, 208-210, 1989. 24: !   Kearsley, Acta Cryst. A, 45, 208-210, 1989.
 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: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,ANGLEAXIS,DEBUG,RMAT) 29: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,ANGLEAXIS,DEBUG,RMAT)
 30: use rigidbodymod 30: use rigidbodymod
 31: USE COMMONS,ONLY : FREEZE, BOXLX, BOXLY, BOXLZ, RBAAT, PULLT, TRAPT, USERPOTT, GTHOMSONT, GTHOMMET, MKTRAPT, NOTRANSROTT ! hk286 31: USE COMMONS,ONLY : FREEZE, BOXLX, BOXLY, BOXLZ, RBAAT, PULLT, TRAPT, USERPOTT, GTHOMSONT, GTHOMMET, MKTRAPT ! hk286
 32: IMPLICIT NONE 32: IMPLICIT NONE
 33: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN 33: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN
 34: DOUBLE PRECISION RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), & 34: DOUBLE PRECISION RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), &
 35:   &              DIAG(4), TEMPA(9*NATOMS), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, & 35:   &              DIAG(4), TEMPA(9*NATOMS), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, &
 36:   &              NCMXB, NCMYB, NCMZB, ROTMAT(3,3), OMEGATOT(3,3) 36:   &              NCMXB, NCMYB, NCMZB, ROTMAT(3,3), OMEGATOT(3,3)
 37: DOUBLE PRECISION :: COCA(3), COCB(3) 37: DOUBLE PRECISION :: COCA(3), COCB(3)
 38: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:), TEMPCOORDS(:,:) 38: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:), TEMPCOORDS(:,:)
 39: COMMON /MINDOM/ ROTMAT, OMEGATOT 39: COMMON /MINDOM/ ROTMAT, OMEGATOT
 40: LOGICAL, INTENT(IN) :: BULKT, TWOD, PRESERVET, DEBUG, ANGLEAXIS 40: LOGICAL, INTENT(IN) :: BULKT, TWOD, PRESERVET, DEBUG, ANGLEAXIS
 41: LOGICAL UPRETURN 41: LOGICAL UPRETURN
180: ! jwrm2> Try this instead180: ! jwrm2> Try this instead
181:       DO J1=1,NATOMS181:       DO J1=1,NATOMS
182:          RB(3*(J1-1)+1)=XB(3*(J1-1)+1)+CMXA182:          RB(3*(J1-1)+1)=XB(3*(J1-1)+1)+CMXA
183:          RB(3*(J1-1)+2)=XB(3*(J1-1)+2)+CMYA183:          RB(3*(J1-1)+2)=XB(3*(J1-1)+2)+CMYA
184:          RB(3*(J1-1)+3)=XB(3*(J1-1)+3)+CMZA184:          RB(3*(J1-1)+3)=XB(3*(J1-1)+3)+CMZA
185:       ENDDO185:       ENDDO
186:    ENDIF186:    ENDIF
187:    RETURN187:    RETURN
188: ENDIF188: ENDIF
189: 189: 
190: IF (NOTRANSROTT) THEN 
191:    DIST=0.0D0 
192:    DO J1=1,NSIZE 
193:       DIST=DIST + (XA(3*(J1-1)+1)-XB(3*(J1-1)+1))**2 & 
194:    &            + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2))**2 & 
195:    &            + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3))**2 
196:    ENDDO 
197:    DIST=SQRT(DIST) 
198:  
199:    RMAT(1:3,1:3)=0.0D0 ! rotation matrix is the identity 
200:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0 
201:    RETURN 
202: ENDIF 
203:  
204: XSHIFT=0.0D0; YSHIFT=0.0D0; ZSHIFT=0.0D0190: XSHIFT=0.0D0; YSHIFT=0.0D0; ZSHIFT=0.0D0
205: NCIT=0191: NCIT=0
206: IF (BULKT) THEN 192: IF (BULKT) THEN 
207: ! 1  NCIT=NCIT+1193: ! 1  NCIT=NCIT+1
208: !    IF (NCIT.GT.1000) THEN194: !    IF (NCIT.GT.1000) THEN
209: !       PRINT '(A)','inertia> WARNING - iterative calculation of centre of mass shift did not converge'195: !       PRINT '(A)','inertia> WARNING - iterative calculation of centre of mass shift did not converge'
210: !    ENDIF196: !    ENDIF
211: !    XSHIFTNEW=0.0D0197: !    XSHIFTNEW=0.0D0
212: !    YSHIFTNEW=0.0D0198: !    YSHIFTNEW=0.0D0
213: !    ZSHIFTNEW=0.0D0199: !    ZSHIFTNEW=0.0D0
303:    RMAT(3,1)=2*(Q2*Q4+Q1*Q3)289:    RMAT(3,1)=2*(Q2*Q4+Q1*Q3)
304:    RMAT(3,2)=2*(Q3*Q4-Q1*Q2)290:    RMAT(3,2)=2*(Q3*Q4-Q1*Q2)
305:    RMAT(3,3)=Q1**2+Q4**2-Q2**2-Q3**2291:    RMAT(3,3)=Q1**2+Q4**2-Q2**2-Q3**2
306: ENDIF292: ENDIF
307: !293: !
308: ! Superpose RB if required294: ! Superpose RB if required
309: ! For the angle-axis formulation we are not translating RB but simply295: ! For the angle-axis formulation we are not translating RB but simply
310: ! setting RB equal to XB, so the origin will be the centre of coordinates296: ! setting RB equal to XB, so the origin will be the centre of coordinates
311: ! including all the rigid body sites. This might need to be changed.297: ! including all the rigid body sites. This might need to be changed.
312: !298: !
313: IF (.NOT.PRESERVET.OR.NOTRANSROTT) THEN299: IF (.NOT.PRESERVET) THEN
314:    CALL NEWROTGEOM(NSIZE,XB,RMAT,CMXA,CMYA,CMZA)300:    CALL NEWROTGEOM(NSIZE,XB,RMAT,CMXA,CMYA,CMZA)
315: 301: 
316:    IF (ANGLEAXIS) THEN302:    IF (ANGLEAXIS) THEN
317:       ALLOCATE(TEMPCOORDS(NSIZE,3))303:       ALLOCATE(TEMPCOORDS(NSIZE,3))
318:       TEMPCOORDS = RESHAPE(XB,SHAPE=(/NSIZE,3/),ORDER=(/2,1/))304:       TEMPCOORDS = RESHAPE(XB,SHAPE=(/NSIZE,3/),ORDER=(/2,1/))
319:       CALL SYSTEMTOAA(NATOMS/2,TEMPCOORDS,RB)305:       CALL SYSTEMTOAA(NATOMS/2,TEMPCOORDS,RB)
320:       DEALLOCATE(TEMPCOORDS)306:       DEALLOCATE(TEMPCOORDS)
321:    ELSE IF (ZUSE(1:1).EQ.'W') THEN307:    ELSE IF (ZUSE(1:1).EQ.'W') THEN
322:       DO J1=1,NATOMS/2308:       DO J1=1,NATOMS/2
323:          CALL CONVERT2(XB(9*J1-8:9*J1-6),XB(9*J1-5:9*J1-3),XB(9*J1-2:9*J1), &309:          CALL CONVERT2(XB(9*J1-8:9*J1-6),XB(9*J1-5:9*J1-3),XB(9*J1-2:9*J1), &


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0