hdiff output

r31037/bulkmindist.f90 2016-08-18 15:30:08.739992051 +0100 r31036/bulkmindist.f90 2016-08-18 15:30:11.516029611 +0100
 10: !   OPTIM is distributed in the hope that it will be useful, 10: !   OPTIM is distributed in the hope that it will be useful,
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE BULKMINDIST(DUMMYB,DUMMYA,XBEST, NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,RESETA, TNMATCH, BMTEST) 19: SUBROUTINE BULKMINDIST(DUMMYB,DUMMYA,XBEST, NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,RESETA, TNMATCH, BMTEST)
 20: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, GEOMDIFFTOL, ATOMMATCHFULL, BULKBOXT, GDSQ, GDSQT 20: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, GEOMDIFFTOL, ATOMMATCHFULL
 21:  21: 
 22: IMPLICIT NONE 22: IMPLICIT NONE
 23: INTEGER J1, NATOMS, NPMIN, NGMIN, J2, PERM(NATOMS), PBEST(NATOMS), NDUMMY, NMATCHED, PATOMS, J3, J4, NMBEST, ND1 23: INTEGER J1, NATOMS, NPMIN, NGMIN, J2, PERM(NATOMS), PBEST(NATOMS), NDUMMY, NMATCHED, PATOMS, J3, J4, NMBEST, ND1
 24: INTEGER NMATCHSV, J5, J6, LPERM(NATOMS), NREP, NREP2 24: INTEGER NMATCHSV, J5, J6, LPERM(NATOMS), NREP, NREP2
 25: DOUBLE PRECISION DUMMYB(3*NATOMS),DUMMYA(3*NATOMS),DISTANCE,BOXLX,BOXLY,BOXLZ,XSHIFT,YSHIFT,ZSHIFT,XTEMP(3*NATOMS) 25: DOUBLE PRECISION DUMMYB(3*NATOMS),DUMMYA(3*NATOMS),DISTANCE,BOXLX,BOXLY,BOXLZ,XSHIFT,YSHIFT,ZSHIFT,XTEMP(3*NATOMS)
 26: DOUBLE PRECISION XBEST(3*NATOMS), DMIN, DTOTAL, DIST 26: DOUBLE PRECISION XBEST(3*NATOMS), DMIN, DTOTAL, DIST, GDSQ
 27: DOUBLE PRECISION DIST2, LDISTANCE, WORSTRAD  27: DOUBLE PRECISION DIST2, LDISTANCE, WORSTRAD 
 28: LOGICAL TWOD,DEBUG,PITEST,SAMEMIN,RESETA, TNMATCH, BMTEST, BMTESTLOCAL  28: LOGICAL TWOD,DEBUG,PITEST,SAMEMIN,RESETA, TNMATCH, BMTEST, BMTESTLOCAL 
 29: COMMON /BULKSHIFT/ XSHIFT,YSHIFT,ZSHIFT 29: COMMON /BULKSHIFT/ XSHIFT,YSHIFT,ZSHIFT
 30: SAVE NMATCHSV 30: SAVE NMATCHSV
 31:  31: 
 32: IF (.NOT.TNMATCH) NMATCHSV=0  32: IF (.NOT.TNMATCH) NMATCHSV=0 
 33: ! 33: !
 34: ! Find smallest group of permutable atoms. 34: ! Find smallest group of permutable atoms.
 35: ! Translate first atom of group to all positions and then find nearest atom within 35: ! Translate first atom of group to all positions and then find nearest atom within
 36: ! the same group for every other atom. 36: ! the same group for every other atom.
 37: ! Keep the best translation/permutation, which corresponds to the smallest 37: ! Keep the best translation/permutation, which corresponds to the smallest
 38: ! minimum image distance. 38: ! minimum image distance.
 39: ! 39: !
 40: !PRINT*, NMATCHSV 40: !PRINT*, NMATCHSV
 41: TNMATCH=.TRUE. 41: TNMATCH=.TRUE.
 42: DISTANCE=1.0D100 42: DISTANCE=1.0D100
 43: PITEST=.FALSE. 43: PITEST=.FALSE.
 44: SAMEMIN=.TRUE. 44: SAMEMIN=.TRUE.
 45: BMTESTLOCAL=BMTEST 45: BMTESTLOCAL=BMTEST
 46: IF (.NOT.GDSQT) THEN 46: GDSQ=GEOMDIFFTOL**2/NATOMS ! because GEOMDIFFTOL is used for the total distance elsewhere in the program
 47:    GDSQ=GEOMDIFFTOL**2/NATOMS ! because GEOMDIFFTOL is used for the total distance elsewhere in the program 
 48: ENDIF 
 49: NPMIN=HUGE(1) 47: NPMIN=HUGE(1)
 50: ! PRINT *,'DUMMYA in bulkmindist:' 48: ! PRINT *,'DUMMYA in bulkmindist:'
 51: ! PRINT '(3G20.10)',DUMMYA(1:3*NATOMS) 49: ! PRINT '(3G20.10)',DUMMYA(1:3*NATOMS)
 52: ! PRINT *,'DUMMYB in bulkmindist:' 50: ! PRINT *,'DUMMYB in bulkmindist:'
 53: ! PRINT '(3G20.10)',DUMMYB(1:3*NATOMS) 51: ! PRINT '(3G20.10)',DUMMYB(1:3*NATOMS)
 54: DO J1=1,NPERMGROUP 52: DO J1=1,NPERMGROUP
 55:    IF (NPERMSIZE(J1).LT.NPMIN) THEN 53:    IF (NPERMSIZE(J1).LT.NPMIN) THEN
 56:       NPMIN=NPERMSIZE(J1) 54:       NPMIN=NPERMSIZE(J1)
 57:       NGMIN=J1 55:       NGMIN=J1
 58:    ENDIF 56:    ENDIF
 59: ENDDO 57: ENDDO
 60: ND1=0 58: ND1=0
 61: DO J1=1,NGMIN-1 59: DO J1=1,NGMIN-1
 62:    ND1=ND1+NPERMSIZE(J1) 60:    ND1=ND1+NPERMSIZE(J1)
 63: ENDDO 61: ENDDO
 64: ! IF (DEBUG) PRINT '(3(A,I6))',' bulkmindist> Smallest group of permutable atoms is number ',NGMIN,' with ',NPMIN,' members' 62: ! IF (DEBUG) PRINT '(3(A,I6))',' bulkmindist> Smallest group of permutable atoms is number ',NGMIN,' with ',NPMIN,' members'
 65: NREP=0 63: NREP=0
 66: NREP2=0 64: NREP2=0
 67: outer: DO J1=ND1+1,ND1+NPMIN 65: outer: DO J1=ND1+1,ND1+NPMIN
 68:    J2=PERMGROUP(J1) 66:    J2=PERMGROUP(J1)
 69:    XSHIFT=DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+1)-BOXLX*NINT((DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+1))/BOXLX) 67:    XSHIFT=DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1))+1)-BOXLX*NINT((DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1))+1))/BOXLX)
 70:    YSHIFT=DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+2)-BOXLY*NINT((DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+2))/BOXLY) 68:    YSHIFT=DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1))+2)-BOXLY*NINT((DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1))+2))/BOXLY)
 71:    IF (.NOT.TWOD) ZSHIFT=DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+3)-BOXLZ*NINT((DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+3))/BOXLZ) 69:    IF (.NOT.TWOD) ZSHIFT=DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1))+3)-BOXLZ*NINT((DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1))+3))/BOXLZ)
 72:    DO J2=1,NATOMS 70:    DO J2=1,NATOMS
 73:       XTEMP(3*(J2-1)+1)=DUMMYA(3*(J2-1)+1)-XSHIFT 71:       XTEMP(3*(J2-1)+1)=DUMMYA(3*(J2-1)+1)-XSHIFT
 74:       XTEMP(3*(J2-1)+2)=DUMMYA(3*(J2-1)+2)-YSHIFT 72:       XTEMP(3*(J2-1)+2)=DUMMYA(3*(J2-1)+2)-YSHIFT
 75:       IF (.NOT.TWOD) XTEMP(3*(J2-1)+3)=DUMMYA(3*(J2-1)+3)-ZSHIFT 73:       IF (.NOT.TWOD) XTEMP(3*(J2-1)+3)=DUMMYA(3*(J2-1)+3)-ZSHIFT
 76:    ENDDO 74:    ENDDO
 77:    NDUMMY=1 75:    NDUMMY=1
 78:    PERM(1:NATOMS)=-1 76:    PERM(1:NATOMS)=-1
 79:    NMATCHED=0 77:    NMATCHED=0
 80:    DTOTAL=0.0D0 78:    DTOTAL=0.0D0
 81:    DO J2=1,NPERMGROUP 79:    DO J2=1,NPERMGROUP
137:         ELSE IF (NMATCHED.EQ.NMATCHSV) THEN135:         ELSE IF (NMATCHED.EQ.NMATCHSV) THEN
138:           NREP=NREP+1136:           NREP=NREP+1
139:           IF (NREP.GT.10) THEN137:           IF (NREP.GT.10) THEN
140:             IF (.NOT.ATOMMATCHFULL) BMTESTLOCAL=.FALSE. 138:             IF (.NOT.ATOMMATCHFULL) BMTESTLOCAL=.FALSE. 
141:             CYCLE outer ! Match is always the same - give up, still do the PI test139:             CYCLE outer ! Match is always the same - give up, still do the PI test
142:           ENDIF140:           ENDIF
143:         ENDIF141:         ENDIF
144:        ELSE IF (J2.EQ.NPERMGROUP) THEN142:        ELSE IF (J2.EQ.NPERMGROUP) THEN
145:         NREP=0 143:         NREP=0 
146:        ENDIF  144:        ENDIF  
147:        IF (J2.EQ.NPERMGROUP.AND.NMATCHED.NE.NATOMS) CYCLE outer 145:       IF (J2.EQ.NPERMGROUP.AND.NMATCHED.NE.NATOMS) CYCLE outer 
148:       ENDIF  146:       ENDIF  
149:    ENDDO147:    ENDDO
150:    IF (SQRT(DTOTAL).GE.GEOMDIFFTOL) CYCLE outer 
151:    IF (SAMEMIN) THEN148:    IF (SAMEMIN) THEN
152:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> identical isomers identified for distance ',SQRT(DTOTAL)149:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> identical isomers identified for distance ',SQRT(DTOTAL)
153:    ELSE150:    ELSE
154:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> permutational isomers identified for distance ',SQRT(DTOTAL)151:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> permutational isomers identified for distance ',SQRT(DTOTAL)
155:    ENDIF152:    ENDIF
156:    PITEST=.TRUE.153:    PITEST=.TRUE.
157:    DISTANCE=DTOTAL154:    DISTANCE=DTOTAL
158:    IF (RESETA .AND. BULKBOXT) THEN155:    IF (RESETA) THEN
159:       DO J2=1,NATOMS 
160:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1) 
161:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2) 
162:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3) 
163:       ENDDO 
164:    ELSE IF (RESETA) THEN 
165:       DO J2=1,NATOMS156:       DO J2=1,NATOMS
166:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1)-BOXLX*NINT(XTEMP(3*(PERM(J2)-1)+1)/BOXLX)157:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1)-BOXLX*NINT(XTEMP(3*(PERM(J2)-1)+1)/BOXLX)
167:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2)-BOXLY*NINT(XTEMP(3*(PERM(J2)-1)+2)/BOXLY)158:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2)-BOXLY*NINT(XTEMP(3*(PERM(J2)-1)+2)/BOXLY)
168:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3)-BOXLZ*NINT(XTEMP(3*(PERM(J2)-1)+3)/BOXLZ)159:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3)-BOXLZ*NINT(XTEMP(3*(PERM(J2)-1)+3)/BOXLZ)
169:       ENDDO160:       ENDDO
170:    ENDIF161:    ENDIF
171: 162: 
172:    RETURN163:    RETURN
173: ENDDO outer164: ENDDO outer
174: 165: 


r31037/key.f90 2016-08-18 15:30:08.991995461 +0100 r31036/key.f90 2016-08-18 15:30:11.760032913 +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, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT 57:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT
 58:  
 59: ! sy349 > for testing the flatpath after dneb 
 60:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 
 61:       LOGICAL FLATPATHT 
 62:  58: 
 63: ! 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)
 64:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 60:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 65:       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)
 66:       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
 67:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 63:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 68:       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
 69:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 65:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 70:  66: 
 71: ! hk286 > generalised THOMSON problem 67: ! hk286 > generalised THOMSON problem
 98:      &        NEBRESEEDDEL2, INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 94:      &        NEBRESEEDDEL2, INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 99:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, & 95:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &
100:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, & 96:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &
101:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, & 97:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &
102:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, & 98:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &
103:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, & 99:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &
104:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &100:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &
105:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &101:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
106:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &102:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
107:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &103:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
108:      &         NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT104:      &        QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT
109: 105: 
110: !     sf344106: !     sf344
111:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &107:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
112:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &108:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
113:      &                     PYLOCALSTEP(2)109:      &                     PYLOCALSTEP(2)
114:  110:  
115:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)111:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)
116:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:)112:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:)
117:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF113:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF
118:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET114:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET


r31037/keywords.f 2016-08-18 15:30:09.247998924 +0100 r31036/keywords.f 2016-08-18 15:30:12.020036430 +0100
545:          DTHRESH=2.0D0545:          DTHRESH=2.0D0
546:          NSTEPNEB=1546:          NSTEPNEB=1
547:          NEBMAG=0547:          NEBMAG=0
548:          NMOVE=1548:          NMOVE=1
549:          NPATHFRAME=0549:          NPATHFRAME=0
550:          FRAMEEDIFF=0.0D0550:          FRAMEEDIFF=0.0D0
551:          FRAMESDIFF=0.0D0551:          FRAMESDIFF=0.0D0
552:          CUBIC=.FALSE.552:          CUBIC=.FALSE.
553:          BULKT=.FALSE.553:          BULKT=.FALSE.
554:          BULK_BOXVEC(:) = 1.d0554:          BULK_BOXVEC(:) = 1.d0
555:          BULKBOXT=.FALSE. 
556:          CUTT = .FALSE.555:          CUTT = .FALSE.
557:          POTENTIAL_CUTOFF = 1.D0556:          POTENTIAL_CUTOFF = 1.D0
558:          SDT=.FALSE.557:          SDT=.FALSE.
559:          SDOXYGEN=0558:          SDOXYGEN=0
560:          SDHYDROGEN=0559:          SDHYDROGEN=0
561:          SDCHARGE=0560:          SDCHARGE=0
562:          BOWMANT=.FALSE.561:          BOWMANT=.FALSE.
563:          TTM3T=.FALSE.562:          TTM3T=.FALSE.
564:          BOWMANPES=2563:          BOWMANPES=2
565:          BOWMANDIR='~/svn/OPTIM/source/Bowman/coef-3b/'564:          BOWMANDIR='~/svn/OPTIM/source/Bowman/coef-3b/'
586:          GSTANTYPE = 1585:          GSTANTYPE = 1
587:          REPARAMTOL = 0.75586:          REPARAMTOL = 0.75
588:          GSGROWTOL = 0.25587:          GSGROWTOL = 0.25
589:          GSCONV = 1.0D-3588:          GSCONV = 1.0D-3
590:          GSMXSTP = 0.1589:          GSMXSTP = 0.1
591:          STOCKT=.FALSE.590:          STOCKT=.FALSE.
592:          STOCKSPIN = .FALSE.591:          STOCKSPIN = .FALSE.
593:          STOCKZTOL = 1.0D-4592:          STOCKZTOL = 1.0D-4
594:          STOCKMAXSPIN = 20593:          STOCKMAXSPIN = 20
595:          GEOMDIFFTOL=1.0D-1594:          GEOMDIFFTOL=1.0D-1
596:          GDSQ=0.0D-1 
597:          GDSQT=.FALSE. 
598:          EDIFFTOL=1.0D-6595:          EDIFFTOL=1.0D-6
599:          NSECDIAG=2596:          NSECDIAG=2
600: 597: 
601:  
602:          ! MSEVB parameters598:          ! MSEVB parameters
603: 599: 
604:          shellsToCount = 3600:          shellsToCount = 3
605:          maxHbondLength = 2.5d0601:          maxHbondLength = 2.5d0
606:          minHbondAngle = 130.0d0602:          minHbondAngle = 130.0d0
607:          OOclash_sq = 4.41d0 ! 2.1^2603:          OOclash_sq = 4.41d0 ! 2.1^2
608:          printCoefficients = .FALSE.604:          printCoefficients = .FALSE.
609: 605: 
610:          NEBRESEEDT=.FALSE.606:          NEBRESEEDT=.FALSE.
611:          NEBRESEEDINT=100607:          NEBRESEEDINT=100
616:          NEBRESEEDPOW1=2612:          NEBRESEEDPOW1=2
617:          NEBRESEEDPOW2=10613:          NEBRESEEDPOW2=10
618:          ADDREPT=.FALSE.614:          ADDREPT=.FALSE.
619: 615: 
620:          INTLJT=.FALSE.616:          INTLJT=.FALSE.
621:          INTLJSTEPS=1000617:          INTLJSTEPS=1000
622:          INTLJTOL=1.0D-3618:          INTLJTOL=1.0D-3
623:          INTLJDEL=0.1D0619:          INTLJDEL=0.1D0
624:          INTLJEPS=1.0D0620:          INTLJEPS=1.0D0
625: 621: 
 622:          FREEZETOL=1.0D-3
626: !623: !
627: ! QCI parameters624: ! QCI parameters
628: !625: !
629:          CONDATT=.FALSE.626:          CONDATT=.FALSE.
630:          QCIPOTT=.FALSE.627:          QCIPOTT=.FALSE.
631:          QCIPOT2T=.FALSE.628:          QCIPOT2T=.FALSE.
632:          QCIADDREP=0629:          QCIADDREP=0
633:          QCIADDREPCUT=1.0D0630:          QCIADDREPCUT=1.0D0
634:          QCIADDREPEPS=1.0D0631:          QCIADDREPEPS=1.0D0
635:          QCINOREPINT=.FALSE.632:          QCINOREPINT=.FALSE.
636:          MAXNACTIVE=0633:          MAXNACTIVE=0
637:  
638:          FREEZETOL=1.0D-3634:          FREEZETOL=1.0D-3
639:          FLATTESTT=.FALSE. 
640:          FLATEDIFF=1.0D-6 
641:          QCIPERMCHECK=.FALSE.635:          QCIPERMCHECK=.FALSE.
642:          QCIPERMCHECKINT=100636:          QCIPERMCHECKINT=100
643:          INTCONSTRAINTT=.FALSE.637:          INTCONSTRAINTT=.FALSE.
644:          INTCONSTRAINTTOL=0.1D0638:          INTCONSTRAINTTOL=0.1D0
645:          INTCONSTRAINTDEL=10.0D0639:          INTCONSTRAINTDEL=10.0D0
646:          INTCONSTRAINTREP=100.0D0640:          INTCONSTRAINTREP=100.0D0
647:          INTCONSTRAINREPCUT=1.7D0641:          INTCONSTRAINREPCUT=1.7D0
648:          INTFREEZET=.FALSE.642:          INTFREEZET=.FALSE.
649:          INTFREEZETOL=1.0D-3643:          INTFREEZETOL=1.0D-3
650:          INTFREEZEMIN=10644:          INTFREEZEMIN=10
770:          OHT=.FALSE.764:          OHT=.FALSE.
771:          IHT=.FALSE.765:          IHT=.FALSE.
772:          TDT=.FALSE.766:          TDT=.FALSE.
773:          D5HT=.FALSE.767:          D5HT=.FALSE.
774:          FOH=0.0D0768:          FOH=0.0D0
775:          FIH=0.0D0769:          FIH=0.0D0
776:          FTD=0.0D0770:          FTD=0.0D0
777:          FD5H=0.0D0771:          FD5H=0.0D0
778:          MAXERISE=1.0D-10772:          MAXERISE=1.0D-10
779:          XMAXERISE=1.0D-3773:          XMAXERISE=1.0D-3
780:          NEBMAXERISE=1.0D-3 
781:          INTEPSILON=1.0D-6774:          INTEPSILON=1.0D-6
782: 775: 
783:          EFIELD=0.0D0776:          EFIELD=0.0D0
784:          COLDFUSIONLIMIT=-1.0D6777:          COLDFUSIONLIMIT=-1.0D6
785:          BLNT=.FALSE.778:          BLNT=.FALSE.
786:          DUMPDATAT=.FALSE.779:          DUMPDATAT=.FALSE.
787:          LOWESTFRQT=.FALSE.780:          LOWESTFRQT=.FALSE.
788:          REDOPATH=.FALSE.781:          REDOPATH=.FALSE.
789:          REDOFRAC=0.5D0782:          REDOFRAC=0.5D0
790:          REDOK=0.0D0783:          REDOK=0.0D0
1643:             BULKT=.TRUE.1636:             BULKT=.TRUE.
1644:             IF (NITEMS.GT.1) THEN1637:             IF (NITEMS.GT.1) THEN
1645:                CALL READF(bulk_boxvec(1))1638:                CALL READF(bulk_boxvec(1))
1646:             ENDIF1639:             ENDIF
1647:             IF (NITEMS.GT.2) THEN1640:             IF (NITEMS.GT.2) THEN
1648:                CALL READF(bulk_boxvec(2))1641:                CALL READF(bulk_boxvec(2))
1649:             ENDIF1642:             ENDIF
1650:             IF (NITEMS.GT.3) THEN1643:             IF (NITEMS.GT.3) THEN
1651:                CALL READF(bulk_boxvec(3))1644:                CALL READF(bulk_boxvec(3))
1652:             ENDIF1645:             ENDIF
1653: ! 
1654: ! Allow particles move out of the box when dealing with bulk system. 
1655: ! 
1656:          ELSE IF (WORD.EQ.'BULKBOX') THEN 
1657:             BULKBOXT=.TRUE. 
1658: ! 1646: ! 
1659: ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC1647: ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1660: ! CADPAC tells the program to read derivative information in1648: ! CADPAC tells the program to read derivative information in
1661: ! CADPAC format.                                        - default FALSE1649: ! CADPAC format.                                        - default FALSE
1662: ! 1650: ! 
1663:          ELSE IF (WORD.EQ.'CADPAC') THEN1651:          ELSE IF (WORD.EQ.'CADPAC') THEN
1664:             CADPAC=.TRUE.1652:             CADPAC=.TRUE.
1665:             CALL READA(SYS)1653:             CALL READA(SYS)
1666:             DO J1=1,801654:             DO J1=1,80
1667:                IF (SYS(J1:J1).EQ.' ') THEN1655:                IF (SYS(J1:J1).EQ.' ') THEN
2838:             FIXATMS = .TRUE.2826:             FIXATMS = .TRUE.
2839: ! 2827: ! 
2840: ! Fix uphill direction until force changes sign.2828: ! Fix uphill direction until force changes sign.
2841: ! T12FAC is the fraction of the first collision time to be used in HSMOVE2829: ! T12FAC is the fraction of the first collision time to be used in HSMOVE
2842: ! 2830: ! 
2843:          ELSE IF (WORD.EQ.'FIXD') THEN2831:          ELSE IF (WORD.EQ.'FIXD') THEN
2844:             FIXD=.TRUE.2832:             FIXD=.TRUE.
2845:             IF (NITEMS.GT.1) CALL READF(T12FAC)2833:             IF (NITEMS.GT.1) CALL READF(T12FAC)
2846:             NMOVE=12834:             NMOVE=1
2847:             IF (NITEMS.GT.2) CALL READF(DTHRESH)2835:             IF (NITEMS.GT.2) CALL READF(DTHRESH)
2848: ! 
2849: ! FLATEDIFF: maximum energy difference when judging whether a dneb path is flat 
2850: ! 
2851:          ELSE IF (WORD.EQ.'FLATTEST') THEN 
2852:             FLATTESTT=.TRUE. 
2853:             IF (NITEMS.GT.1) CALL READF(FLATEDIFF) 
2854:  
2855: ! 2836: ! 
2856: ! FRACTIONAL: constant pressure calculation using fractional coordinates2837: ! FRACTIONAL: constant pressure calculation using fractional coordinates
2857: ! 2838: ! 
2858:          ELSE IF (WORD.EQ.'FRACTIONAL') THEN2839:          ELSE IF (WORD.EQ.'FRACTIONAL') THEN
2859:             FRACTIONAL=.TRUE.2840:             FRACTIONAL=.TRUE.
2860: 2841: 
2861: ! 2842: ! 
2862: ! Frozen atoms.2843: ! Frozen atoms.
2863: ! 2844: ! 
2864:          ELSE IF (WORD.EQ.'FREEZE') THEN2845:          ELSE IF (WORD.EQ.'FREEZE') THEN
3061: ! DTEST=.TRUE.3042: ! DTEST=.TRUE.
3062: ! CALL READF(PCUT)3043: ! CALL READF(PCUT)
3063: ! CALL READI(NDIIA)3044: ! CALL READI(NDIIA)
3064: ! CALL READI(NINTV)3045: ! CALL READI(NINTV)
3065: ! ENDIF3046: ! ENDIF
3066: ! IF (NDIIA.GT.NDIIS) THEN3047: ! IF (NDIIA.GT.NDIIS) THEN
3067: ! WRITE(*,'(A,I6)') ' NDIIA too large=',NDIIA3048: ! WRITE(*,'(A,I6)') ' NDIIA too large=',NDIIA
3068: ! STOP3049: ! STOP
3069: ! ENDIF3050: ! ENDIF
3070: ! 3051: ! 
3071:          ELSE IF (WORD.EQ.'GDSQ') THEN 
3072:             GDSQT=.TRUE. 
3073:             CALL READF(GDSQ) 
3074: ! GEOMDIFFTOL specifies the maximum displacement between identical permutational isomers in connect.3052: ! GEOMDIFFTOL specifies the maximum displacement between identical permutational isomers in connect.
3075: ! 3053: ! 
3076:          ELSE IF (WORD.EQ.'GEOMDIFFTOL') THEN3054:          ELSE IF (WORD.EQ.'GEOMDIFFTOL') THEN
3077:             CALL READF(GEOMDIFFTOL)3055:             CALL READF(GEOMDIFFTOL)
3078: ! 3056: ! 
3079: ! General LJ for mixed systems3057: ! General LJ for mixed systems
3080: ! 3058: ! 
3081:          ELSE IF (WORD.EQ.'GLJ') THEN3059:          ELSE IF (WORD.EQ.'GLJ') THEN
3082:             GLJT=.TRUE.3060:             GLJT=.TRUE.
3083:             NGLJ=NITEMS-1 ! the number of different species, since we specify how many of each on this line3061:             NGLJ=NITEMS-1 ! the number of different species, since we specify how many of each on this line
3692: ! The deafult is 3.3670: ! The deafult is 3.
3693: ! 3671: ! 
3694:          ELSE IF (WORD.EQ.'MAXCON') THEN3672:          ELSE IF (WORD.EQ.'MAXCON') THEN
3695:             CALL READI(MAXCONUSE)3673:             CALL READI(MAXCONUSE)
3696: ! 3674: ! 
3697: ! The maximum energy increase above which mylbfgs will reject a proposed step.3675: ! The maximum energy increase above which mylbfgs will reject a proposed step.
3698: ! 3676: ! 
3699:          ELSE IF (WORD.EQ.'MAXERISE') THEN3677:          ELSE IF (WORD.EQ.'MAXERISE') THEN
3700:             CALL READF(MAXERISE)3678:             CALL READF(MAXERISE)
3701:             IF (NITEMS.GT.1) CALL READF(XMAXERISE)3679:             IF (NITEMS.GT.1) CALL READF(XMAXERISE)
3702: 3680: ! 
3703: ! Maximum number of failures allowed in a minimisation before giving up.3681: ! Maximum number of failures allowed in a minimisation before giving up.
3704: ! 3682: ! 
3705:          ELSE IF (WORD.EQ.'MAXFAIL') THEN3683:          ELSE IF (WORD.EQ.'MAXFAIL') THEN
3706:             IF (NITEMS.GT.1) CALL READI(NFAILMAX)3684:             IF (NITEMS.GT.1) CALL READI(NFAILMAX)
3707: !3685: !
3708: ! Specify single connection attempt for largest gap in current Dijkstra selected chain.3686: ! Specify single connection attempt for largest gap in current Dijkstra selected chain.
3709: ! 3687: ! 
3710:          ELSE IF (WORD.EQ.'MAXGAP') THEN3688:          ELSE IF (WORD.EQ.'MAXGAP') THEN
3711:             MAXGAPT=.TRUE.3689:             MAXGAPT=.TRUE.
3712: ! 3690: ! 
4180:             IF (NITEMS.GT.2) CALL READI(NIMAGE)4158:             IF (NITEMS.GT.2) CALL READI(NIMAGE)
4181:             IF (NITEMS.GT.3) CALL READF(RMSNEB)4159:             IF (NITEMS.GT.3) CALL READF(RMSNEB)
4182:             IF (NITEMS.GT.3) CALL READI(NEBMAG)4160:             IF (NITEMS.GT.3) CALL READI(NEBMAG)
4183:          ELSE IF (WORD == 'NEBK') THEN4161:          ELSE IF (WORD == 'NEBK') THEN
4184:             CALL READF(NEBK)4162:             CALL READF(NEBK)
4185:             NEBKFINAL=NEBK4163:             NEBKFINAL=NEBK
4186:             NEBKINITIAL=NEBK4164:             NEBKINITIAL=NEBK
4187:             NEBFACTOR=1.01D04165:             NEBFACTOR=1.01D0
4188:             IF (NITEMS.GT.2) CALL READF(NEBKFINAL)4166:             IF (NITEMS.GT.2) CALL READF(NEBKFINAL)
4189:             IF (NITEMS.GT.3) CALL READF(NEBFACTOR)4167:             IF (NITEMS.GT.3) CALL READF(NEBFACTOR)
4190:  
4191:          ELSE IF (WORD.EQ.'NEBMAXERISE') THEN 
4192:             CALL READF(NEBMAXERISE) 
4193: !  
4194: ! 4168: ! 
4195: ! Read dneb guess images from file GUESSFILE, default name guess.xyz4169: ! Read dneb guess images from file GUESSFILE, default name guess.xyz
4196: ! 4170: ! 
4197:          ELSE IF ((WORD.EQ.'NEBREADGUESS').OR.(WORD.EQ.'NEBREADGUESSPERM')) THEN4171:          ELSE IF ((WORD.EQ.'NEBREADGUESS').OR.(WORD.EQ.'NEBREADGUESSPERM')) THEN
4198:             IF (WORD.EQ.'NEBREADGUESSPERM') PERMGUESS=.TRUE.4172:             IF (WORD.EQ.'NEBREADGUESSPERM') PERMGUESS=.TRUE.
4199:             READGUESS=.TRUE.4173:             READGUESS=.TRUE.
4200:             NEWNEBT=.TRUE.4174:             NEWNEBT=.TRUE.
4201:             FCD=.TRUE.4175:             FCD=.TRUE.
4202:             IF (NITEMS.GT.1) CALL READA(GUESSFILE)4176:             IF (NITEMS.GT.1) CALL READA(GUESSFILE)
4203:             FILENAME=GUESSFILE4177:             FILENAME=GUESSFILE
5780: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using5754: ! SQVV allows the first NIterSQVVGuessMax DNEB iterations to be done using
5781: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations5755: ! SQVV - switches to LBFGS minimisation after NIterSQVVGuessMax iterations
5782: ! or if the RMS force goes below SQVVGuessRMSTol.5756: ! or if the RMS force goes below SQVVGuessRMSTol.
5783: ! 5757: ! 
5784:          ELSE IF (WORD == 'SQVV') THEN5758:          ELSE IF (WORD == 'SQVV') THEN
5785:             SQVVGUESS=.TRUE.5759:             SQVVGUESS=.TRUE.
5786:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)5760:             IF (NITEMS.GT.1) CALL READI(NITERSQVVGUESSMAX)
5787:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)5761:             IF (NITEMS.GT.2) CALL READF(SQVVGUESSRMSTOL)
5788:          ELSE IF (WORD.EQ.'SSH') THEN5762:          ELSE IF (WORD.EQ.'SSH') THEN
5789:             SSHT=.TRUE.5763:             SSHT=.TRUE.
5790: !5764: 
5791: ! KLIM sets the radius in then wave vector space. 
5792: ! SCA sets the multiplier on the stealthy potential and gradients. 
5793: ! 
5794:          ELSE IF (WORD.EQ.'STEALTHY') THEN5765:          ELSE IF (WORD.EQ.'STEALTHY') THEN
5795:             STEALTHYT=.TRUE.5766:             STEALTHYT=.TRUE.
5796:             CALL READF(KLIM)5767:             CALL READF(KLIM)
5797:             IF (NITEMS.GT.2) THEN5768:             IF (NITEMS.GT.2) THEN
5798:                CALL READF(SCA)5769:                CALL READF(SCA)
5799:             ELSE5770:             ELSE
5800:                SCA=15771:                SCA=1
5801:             END IF5772:             END IF
5802: 5773: 
5803:          ELSE IF (WORD.EQ.'STEALTHYTEST') THEN5774:          ELSE IF (WORD.EQ.'STEALTHYTEST') THEN


r31037/lbfgs.f90 2016-08-18 15:30:06.635963584 +0100 r31036/lbfgs.f90 2016-08-18 15:30:10.760019382 +0100
 28:      USE NEBUTILS 28:      USE NEBUTILS
 29:      USE GRADIENTS 29:      USE GRADIENTS
 30:      USE NEBOUTPUT 30:      USE NEBOUTPUT
 31:      USE PORFUNCS 31:      USE PORFUNCS
 32:      USE MODCHARMM, ONLY : CHRMMT,CHECKOMEGAT 32:      USE MODCHARMM, ONLY : CHRMMT,CHECKOMEGAT
 33:      USE KEY, ONLY : FREEZENODEST, FREEZETOL, DESMDEBUG, MAXNEBBFGS, NEBMUPDATE, NEBDGUESS, & 33:      USE KEY, ONLY : FREEZENODEST, FREEZETOL, DESMDEBUG, MAXNEBBFGS, NEBMUPDATE, NEBDGUESS, &
 34:           & DESMAXEJUMP,INTEPSILON, DESMAXAVGE, KADJUSTFRAC, KADJUSTFRQ, DNEBSWITCH, KADJUSTTOL, NEBRESEEDT, & 34:           & DESMAXEJUMP,INTEPSILON, DESMAXAVGE, KADJUSTFRAC, KADJUSTFRQ, DNEBSWITCH, KADJUSTTOL, NEBRESEEDT, &
 35:           & NEBRESEEDINT, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBKFINAL, NEBFACTOR, & 35:           & NEBRESEEDINT, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBKFINAL, NEBFACTOR, &
 36:           & NREPMAX, ORDERI, ORDERJ, EPSALPHA, NREPULSIVE, DISTREF, ADDREPT, NEBKINITIAL, & 36:           & NREPMAX, ORDERI, ORDERJ, EPSALPHA, NREPULSIVE, DISTREF, ADDREPT, NEBKINITIAL, &
 37:           & NEBRESEEDDEL1, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, REPPOW, & 37:           & NEBRESEEDDEL1, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, REPPOW, &
 38:           & BULKT, REDOTSIM, NREPMAX, COLDFUSIONLIMIT, DNEBEFRAC,VARIABLES, NEBMAXERISE 38:           & BULKT, REDOTSIM, NREPMAX, COLDFUSIONLIMIT, DNEBEFRAC,VARIABLES, MAXERISE
 39:      USE INTCOMMONS, ONLY : DESMINT, NNZ, KD, NINTC, PREVDIH, DIHINFO 39:      USE INTCOMMONS, ONLY : DESMINT, NNZ, KD, NINTC, PREVDIH, DIHINFO
 40:      USE COMMONS, ONLY: REDOPATHNEB 40:      USE COMMONS, ONLY: REDOPATHNEB
 41: ! hk286 41: ! hk286
 42:      USE GENRIGID 42:      USE GENRIGID
 43:  43: 
 44:      IMPLICIT NONE  44:      IMPLICIT NONE 
 45:  45: 
 46:      INTEGER,INTENT(IN) :: D,U  ! DIMENSIONALITY OF THE PROBLEM AND NUMBER OF UPDATES 46:      INTEGER,INTENT(IN) :: D,U  ! DIMENSIONALITY OF THE PROBLEM AND NUMBER OF UPDATES
 47:      INTEGER NPERSIST           ! number of persistent minima 47:      INTEGER NPERSIST           ! number of persistent minima
 48:      INTEGER PERSISTTHRESH      ! persistence threshold 48:      INTEGER PERSISTTHRESH      ! persistence threshold
263:      ENDIF263:      ENDIF
264:      NDECREASE=0264:      NDECREASE=0
265: 20   X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)265: 20   X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
266: 266: 
267:      IF (PREVGRAD.LT.DNEBSWITCH) THEN267:      IF (PREVGRAD.LT.DNEBSWITCH) THEN
268:         CALL OLDNEBGRADIENT268:         CALL OLDNEBGRADIENT
269:      ELSE269:      ELSE
270:         CALL NEBGRADIENT270:         CALL NEBGRADIENT
271:      ENDIF271:      ENDIF
272: 272: 
273:      DO WHILE ((ETOTAL-ENOLD).GT.NEBMAXERISE)273:      DO WHILE ((ETOTAL-ENOLD).GT.MAXERISE)
274:         X(1:D) = X(1:D) - STP(1:D)*SEARCHSTEP(POINT,1:D)274:         X(1:D) = X(1:D) - STP(1:D)*SEARCHSTEP(POINT,1:D)
275:         STP(1:D) = STP(1:D)/10275:         STP(1:D) = STP(1:D)/10
276:         X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)276:         X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
277:         IF (PREVGRAD.LT.DNEBSWITCH) THEN277:         IF (PREVGRAD.LT.DNEBSWITCH) THEN
278:            CALL OLDNEBGRADIENT278:            CALL OLDNEBGRADIENT
279:         ELSE279:         ELSE
280:            CALL NEBGRADIENT280:            CALL NEBGRADIENT
281:         END IF281:         END IF
282:      END DO282:      END DO
283: 283: 


r31037/minpermdist.f90 2016-08-18 15:30:09.496002280 +0100 r31036/minpermdist.f90 2016-08-18 15:30:12.632044710 +0100
325:       STOP325:       STOP
326:     ENDIF326:     ENDIF
327:     IF (INTDISTANCET) THEN327:     IF (INTDISTANCET) THEN
328:       CALL INTDISTANCE(COORDSB, COORDSA, DISTANCE, DEBUG)328:       CALL INTDISTANCE(COORDSB, COORDSA, DISTANCE, DEBUG)
329:       IF (DEBUG) PRINT*, "msb50 minpermdist using intdistance", DISTANCE329:       IF (DEBUG) PRINT*, "msb50 minpermdist using intdistance", DISTANCE
330:     ENDIF330:     ENDIF
331:    DISTANCE=SQRT(DISTANCE)331:    DISTANCE=SQRT(DISTANCE)
332:    RETURN332:    RETURN
333: 333: 
334: ENDIF334: ENDIF
335:  
336: !335: !
337: !  Calculate original centres of mass. sn402 - this is actually the centre of geometry?336: !  Calculate original centres of mass. sn402 - this is actually the centre of geometry?
338: !337: !
339: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0338: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0
340: IF ((NFREEZE.LE.0).AND.(.NOT.MIEFT).AND.(.NOT.MKTRAPT) .AND. (.NOT. NOTRANSROTT)) THEN339: IF ((NFREEZE.LE.0).AND.(.NOT.MIEFT).AND.(.NOT.MKTRAPT) .AND. (.NOT. NOTRANSROTT)) THEN
341:    IF (STOCKT) THEN 340:    IF (STOCKT) THEN 
342:       NRB=(NATOMS/2)341:       NRB=(NATOMS/2)
343:       DO J1=1,NRB342:       DO J1=1,NRB
344:          CMAX=CMAX+COORDSA(3*(J1-1)+1)343:          CMAX=CMAX+COORDSA(3*(J1-1)+1)
345:          CMAY=CMAY+COORDSA(3*(J1-1)+2)344:          CMAY=CMAY+COORDSA(3*(J1-1)+2)
862:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &861:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
863:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &862:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
864:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2863:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
865:       ENDDO864:       ENDDO
866:    ELSE IF (STOCKT) THEN865:    ELSE IF (STOCKT) THEN
867:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)866:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)
868:       XDUMMY=0.0D0867:       XDUMMY=0.0D0
869:       DO J1=1,(NATOMS/2)868:       DO J1=1,(NATOMS/2)
870:          XBEST(3*(J1-1)+1)=XBEST(3*(J1-1)+1)+CMBX869:          XBEST(3*(J1-1)+1)=XBEST(3*(J1-1)+1)+CMBX
871:          XBEST(3*(J1-1)+2)=XBEST(3*(J1-1)+2)+CMBY870:          XBEST(3*(J1-1)+2)=XBEST(3*(J1-1)+2)+CMBY
872:          XBEST(3*(J1-1)+3)=XBEST(3*(J1-1)+3)+CMBZ871:          xbest(3*(J1-1)+3)=XBEST(3*(J1-1)+3)+CMBZ
873:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &872:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
874:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &873:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
875:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2874:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
876:       ENDDO875:       ENDDO
877:    ELSEIF (BULKT) THEN876:    ELSEIF (BULKT) THEN
878:       XDUMMY=0.0D0877:       XDUMMY=0.0D0
879:       DO J1=1,NATOMS878:       DO J1=1,NATOMS
880:          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+ &879:          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+ &
881:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2) - BOXLY*NINT((COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))/BOXLY))**2880:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2) - BOXLY*NINT((COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))/BOXLY))**2
882:          IF (.NOT.TWOD) XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3) - &881:          IF (.NOT.TWOD) XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3) - &


r31037/ncutils.f90 2016-08-18 15:30:06.135956818 +0100 r31036/ncutils.f90 2016-08-18 15:30:10.248012455 +0100
 63:                          RETURN 63:                          RETURN
 64:                     ENDIF 64:                     ENDIF
 65:                ENDIF 65:                ENDIF
 66:           ENDDO 66:           ENDDO
 67:  67: 
 68:           NCISNEWTS=.TRUE. 68:           NCISNEWTS=.TRUE.
 69:      END FUNCTION NCISNEWTS 69:      END FUNCTION NCISNEWTS
 70:      70:     
 71:      SUBROUTINE ISNEWMIN(E,COORD,MINPOS,NEW,REDOPATH,PERMUTE,INVERT,INDEX,I) 71:      SUBROUTINE ISNEWMIN(E,COORD,MINPOS,NEW,REDOPATH,PERMUTE,INVERT,INDEX,I)
 72:           USE KEYUTILS 72:           USE KEYUTILS
 73:           USE KEY,ONLY : RIGIDBODY, BULKT, TWOD, PERMDIST, BULK_BOXVEC, BULKBOXT, FLATTESTT !msb50 73:           USE KEY,ONLY : RIGIDBODY, BULKT, TWOD, PERMDIST !msb50
 74:           USE NEWNEBMODULE 74:           USE COMMONS,ONLY : PARAM1,PARAM2,PARAM3,DEBUG
 75:           USE COMMONS,ONLY : PARAM1,PARAM2,PARAM3,DEBUG,NATOMS 
 76:           USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT 75:           USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT
 77:           IMPLICIT NONE 76:           IMPLICIT NONE
 78:            77:           
 79:           DOUBLE PRECISION,POINTER     :: E,COORD(:)  78:           DOUBLE PRECISION,POINTER     :: E,COORD(:) 
 80:           INTEGER,INTENT(OUT) :: MINPOS ! IF MINIMUM NEW IT IS EQUAL TO NMIN+1; OTHERWISE - POSITION IN MIN ARRAY 79:           INTEGER,INTENT(OUT) :: MINPOS ! IF MINIMUM NEW IT IS EQUAL TO NMIN+1; OTHERWISE - POSITION IN MIN ARRAY
 81:           LOGICAL,INTENT(OUT) :: NEW 80:           LOGICAL,INTENT(OUT) :: NEW
 82:           DOUBLE PRECISION QTEMP(NOPT) 81:           DOUBLE PRECISION QTEMP(NOPT)
 83:           INTEGER :: I, J1 82:           INTEGER :: I
 84:           INTEGER INVERT, INDEX(NATOMS), J2 83:           INTEGER INVERT, INDEX(NATOMS), J2
 85:           LOGICAL PERMUTE,SUCCESS,REDOPATH 84:           LOGICAL PERMUTE,SUCCESS,REDOPATH
 86:           DOUBLE PRECISION D, DIST2, RMAT(3,3), TSRDHE(NOPT) 85:           DOUBLE PRECISION D, DIST2, RMAT(3,3)
 87:  86: 
 88:           MINPOS=NMIN+1 87:           MINPOS=NMIN+1
 89:           NEW=.TRUE. 88:           NEW=.TRUE.
 90:  89:  
 91:           DO I=1,NMIN 90:           DO I=1,NMIN
 92:              PERMUTE=.FALSE. 91:              PERMUTE=.FALSE.
 93:              IF (DEBUG) PRINT '(A,I6,2G20.10)',' isnewmin> I,E,MI(I)%DATA%E=',I,E,MI(I)%DATA%E 92:              IF (DEBUG) PRINT '(A,I6,2G20.10)',' isnewmin> I,E,MI(I)%DATA%E=',I,E,MI(I)%DATA%E
 94:              IF (ABS(E-MI(I)%DATA%E) < EDIFFTOL) THEN 93:              IF (ABS(E-MI(I)%DATA%E) < EDIFFTOL) THEN
 95:                 CALL MINPERMDIST(MI(I)%DATA%X(1:NOPT), COORD, NATOMS, & 94:                 CALL MINPERMDIST(COORD,MI(I)%DATA%X(1:NOPT), NATOMS, &
 96:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT) 95:   &                                DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
 97:                 IF (DEBUG) PRINT '(A,G20.10)','isnewmin> minimum distance=',D 96:                 IF (DEBUG) PRINT '(A,G20.10)','isnewmin> minimum distance=',D
 98:                 IF (D<GEOMDIFFTOL) THEN 97:                 IF (D<GEOMDIFFTOL) THEN
 99:                    NEW=.FALSE. 98:                    NEW=.FALSE.
100:                    MINPOS=I 99:                    MINPOS=I
101: !                  IF (DEBUG) PRINT '(A,I6)','isnewmin> MINPOS=',I100: !                  IF (DEBUG) PRINT '(A,I6)','isnewmin> MINPOS=',I
102:                    RETURN101:                    RETURN
103:                 ENDIF102:                 ENDIF
104:              ENDIF103:              ENDIF
105:           ENDDO104:           ENDDO
194:           NMIN=NMIN+1193:           NMIN=NMIN+1
195:           MI(NMIN)%DATA%E => E194:           MI(NMIN)%DATA%E => E
196:           MI(NMIN)%DATA%X => COORD195:           MI(NMIN)%DATA%X => COORD
197: !         DO I=1,NMIN196: !         DO I=1,NMIN
198: !           PRINT '(A,I6,F20.10)',' addnewmin> I,energy=',I,MI(I)%DATA%E197: !           PRINT '(A,I6,F20.10)',' addnewmin> I,energy=',I,MI(I)%DATA%E
199: !         ENDDO198: !         ENDDO
200:           ALLOCATE( MI(NMIN)%DATA%D(NMIN-1),MI(NMIN)%DATA%NTRIES(NMIN-1),MI(NMIN)%DATA%C(1),MI(NMIN)%DATA%INTNTRIES(NMIN-1) )199:           ALLOCATE( MI(NMIN)%DATA%D(NMIN-1),MI(NMIN)%DATA%NTRIES(NMIN-1),MI(NMIN)%DATA%C(1),MI(NMIN)%DATA%INTNTRIES(NMIN-1) )
201:           IF (INTERPCOSTFUNCTION) ALLOCATE( MI(NMIN)%DATA%INTERP(NMIN-1))200:           IF (INTERPCOSTFUNCTION) ALLOCATE( MI(NMIN)%DATA%INTERP(NMIN-1))
202: 201: 
203:           DO I=1,NMIN-1202:           DO I=1,NMIN-1
204:                CALL MINPERMDIST(MI(I)%DATA%X(:),MI(NMIN)%DATA%X(:), NATOMS, &203:                CALL MINPERMDIST(MI(NMIN)%DATA%X(:), MI(I)%DATA%X(:), NATOMS, &
205:   &                             DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)204:   &                             DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
206:                IF (DEBUG) PRINT '(A,G20.10,2(A,I6))',' addnewmin> distance from MINPERMDIST/NEWMINDIST=',D,' for ',I,' and ',NMIN205:                IF (DEBUG) PRINT '(A,G20.10,2(A,I6))',' addnewmin> distance from MINPERMDIST/NEWMINDIST=',D,' for ',I,' and ',NMIN
207:                MI(NMIN)%DATA%D(I)=D  !  PASSING MIN(NMIN)%DATA%D(I) INSTEAD OF D DOES NOT WORK206:                MI(NMIN)%DATA%D(I)=D  !  PASSING MIN(NMIN)%DATA%D(I) INSTEAD OF D DOES NOT WORK
208:                                      !  perhaps because min is an intrisic function!207:                                      !  perhaps because min is an intrisic function!
209:                IF (D.LT.GEOMDIFFTOL) THEN208:                IF (D.LT.GEOMDIFFTOL) THEN
210:  
211:                   IF (ABS(MI(I)%DATA%E - MI(NMIN)%DATA%E) .LT. EDIFFTOL) THEN209:                   IF (ABS(MI(I)%DATA%E - MI(NMIN)%DATA%E) .LT. EDIFFTOL) THEN
212:                      PRINT '(A,I8,F12.6,A)',' addnewmin> Distance to minimum ',I,D,' is less than tolerance - this should never happen'210:                      PRINT '(A,I8,F12.6,A)',' addnewmin> Distance to minimum ',I,D,' is less than tolerance - this should never happen'
213:                      STOP ! DJW debug211:                      STOP ! DJW debug
214:                   ELSE212:                   ELSE
215:                      PRINT '(A,F12.6,A,2G20.10)',' addnewmin> Distance between minima ',D,' is less than tolerance' 213:                      PRINT '(A,F12.6,A,2G20.10)',' addnewmin> Distance between minima ',D,' is less than tolerance' 
216:                      PRINT '(A,2G20.10)',' addnewmin> Energies are ', MI(I)%DATA%E, MI(NMIN)%DATA%E 214:                      PRINT '(A,2G20.10)',' addnewmin> Energies are ', MI(I)%DATA%E, MI(NMIN)%DATA%E 
217:                   ENDIF215:                   ENDIF
218:                ENDIF216:                ENDIF
219:                IF (INTERPCOSTFUNCTION) THEN217:                IF (INTERPCOSTFUNCTION) THEN
220: !                   IF (INTCONSTRAINTT) THEN218: !                   IF (INTCONSTRAINTT) THEN


r31037/newmindist.f90 2016-08-18 15:30:09.744005636 +0100 r31036/newmindist.f90 2016-08-18 15:30:12.880048067 +0100
 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, NOTRANSROTT, RBAAT, PULLT, EFIELDT, EYTRAPT, MKTRAPT, USERPOTT, &
 36:   &            GTHOMSONT, VARIABLES, BULKBOXT! 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
 46: CHARACTER(LEN=5) ZUSE 46: CHARACTER(LEN=5) ZUSE
162:    &               + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2))**2 &162:    &               + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2))**2 &
163:    &               + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3))**2163:    &               + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3))**2
164:       ENDDO164:       ENDDO
165:    ENDIF165:    ENDIF
166:    DIST=SQRT(DIST)166:    DIST=SQRT(DIST)
167:    167:    
168:    RMAT(1:3,1:3)=0.0D0 ! rotation matrix is the identity168:    RMAT(1:3,1:3)=0.0D0 ! rotation matrix is the identity
169:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0169:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0
170:    RETURN170:    RETURN
171: ENDIF171: ENDIF
172: !172: IF (BULKT) THEN
173: ! If bulkboxt is on, the particles are allow to move out of the box, but 
174: ! structure still remains a box. 
175: ! 
176: IF (BULKT.AND.(.NOT.BULKBOXT)) THEN 
177:    BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3173:    BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3
178:    DO J1=1,NSIZE174:    DO J1=1,NSIZE
179:       XA(3*(J1-1)+1)=XA(3*(J1-1)+1)-BOXLX*NINT(XA(3*(J1-1)+1)/BOXLX)175:       XA(3*(J1-1)+1)=XA(3*(J1-1)+1)-BOXLX*NINT(XA(3*(J1-1)+1)/BOXLX)
180:       XA(3*(J1-1)+2)=XA(3*(J1-1)+2)-BOXLY*NINT(XA(3*(J1-1)+2)/BOXLY)176:       XA(3*(J1-1)+2)=XA(3*(J1-1)+2)-BOXLY*NINT(XA(3*(J1-1)+2)/BOXLY)
181:       IF (.NOT.TWOD) XA(3*(J1-1)+3)=XA(3*(J1-1)+3)-BOXLZ*NINT(XA(3*(J1-1)+3)/BOXLZ)177:       IF (.NOT.TWOD) XA(3*(J1-1)+3)=XA(3*(J1-1)+3)-BOXLZ*NINT(XA(3*(J1-1)+3)/BOXLZ)
182:    ENDDO178:    ENDDO
183:    DO J1=1,NSIZE179:    DO J1=1,NSIZE
184:       XB(3*(J1-1)+1)=XB(3*(J1-1)+1)-BOXLX*NINT(XB(3*(J1-1)+1)/BOXLX)180:       XB(3*(J1-1)+1)=XB(3*(J1-1)+1)-BOXLX*NINT(XB(3*(J1-1)+1)/BOXLX)
185:       XB(3*(J1-1)+2)=XB(3*(J1-1)+2)-BOXLY*NINT(XB(3*(J1-1)+2)/BOXLY)181:       XB(3*(J1-1)+2)=XB(3*(J1-1)+2)-BOXLY*NINT(XB(3*(J1-1)+2)/BOXLY)
186:       IF (.NOT.TWOD) XB(3*(J1-1)+3)=XB(3*(J1-1)+3)-BOXLZ*NINT(XB(3*(J1-1)+3)/BOXLZ)182:       IF (.NOT.TWOD) XB(3*(J1-1)+3)=XB(3*(J1-1)+3)-BOXLZ*NINT(XB(3*(J1-1)+3)/BOXLZ)
207:    CMXB=CMXB+XB(3*(J1-1)+1)203:    CMXB=CMXB+XB(3*(J1-1)+1)
208:    CMYB=CMYB+XB(3*(J1-1)+2)204:    CMYB=CMYB+XB(3*(J1-1)+2)
209:    CMZB=CMZB+XB(3*(J1-1)+3)205:    CMZB=CMZB+XB(3*(J1-1)+3)
210: ENDDO206: ENDDO
211: CMXB=CMXB/NSIZE; CMYB=CMYB/NSIZE; CMZB=CMZB/NSIZE207: CMXB=CMXB/NSIZE; CMYB=CMYB/NSIZE; CMZB=CMZB/NSIZE
212: DO J1=1,NSIZE208: DO J1=1,NSIZE
213:    XB(3*(J1-1)+1)=XB(3*(J1-1)+1)-CMXB209:    XB(3*(J1-1)+1)=XB(3*(J1-1)+1)-CMXB
214:    XB(3*(J1-1)+2)=XB(3*(J1-1)+2)-CMYB210:    XB(3*(J1-1)+2)=XB(3*(J1-1)+2)-CMYB
215:    XB(3*(J1-1)+3)=XB(3*(J1-1)+3)-CMZB211:    XB(3*(J1-1)+3)=XB(3*(J1-1)+3)-CMZB
216: ENDDO212: ENDDO
217: !ENDIF 
218: ! PRINT '(A,6F15.7)','CMA,CMB=',CMXA,CMYA,CMZA,CMXB,CMYB,CMZB213: ! PRINT '(A,6F15.7)','CMA,CMB=',CMXA,CMYA,CMZA,CMXB,CMYB,CMZB
219: 214: 
220: !  215: !  
221: ! Deal with TWOD non-bulk here.216: ! Deal with TWOD non-bulk here.
222: !217: !
223: IF (TWOD.AND.(.NOT.BULKT)) THEN218: IF (TWOD.AND.(.NOT.BULKT)) THEN
224:    CALL MINDIST(XA,XB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET)219:    CALL MINDIST(XA,XB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET)
225:    RMAT(1:3,1:3)=OMEGATOT(1:3,1:3)220:    RMAT(1:3,1:3)=OMEGATOT(1:3,1:3)
226: !  221: !  
227: !  To align the structures for minpermdist we use the orientation in XA222: !  To align the structures for minpermdist we use the orientation in XA
396: !391: !
397: !  Needs some thought for the angle/axis rigid body formulation.392: !  Needs some thought for the angle/axis rigid body formulation.
398: !393: !
399:       PRINT '(A)',' newmindist> WARNING *** back transformation not programmed yet for rigid bodies'394:       PRINT '(A)',' newmindist> WARNING *** back transformation not programmed yet for rigid bodies'
400: !395: !
401: !  Preserve centre of mass for EYTRAP, but not orientatation.396: !  Preserve centre of mass for EYTRAP, but not orientatation.
402: !397: !
403: !  ELSEIF (EYTRAPT.OR.MKTRAPT) THEN398: !  ELSEIF (EYTRAPT.OR.MKTRAPT) THEN
404:    ELSEIF (EYTRAPT) THEN399:    ELSEIF (EYTRAPT) THEN
405:       CALL NEWROTGEOM(NSIZE,RB,RMAT,0.0D0,0.0D0,0.0D0)400:       CALL NEWROTGEOM(NSIZE,RB,RMAT,0.0D0,0.0D0,0.0D0)
406:    !ELSEIF (BULKT) THEN 
407:    !   DO J1=1,NSIZE 
408:    !      RB(3*(J1-1)+1)=RB(3*(J1-1)+1)+VPBX+ZSHIFT 
409:    !      RB(3*(J1-1)+2)=RB(3*(J1-1)+2)+VPBY+YSHIFT 
410:    !      RB(3*(J1-1)+3)=RB(3*(J1-1)+3)+VPBZ+ZSHIFT 
411:    !      RB(3*(J1-1)+1)=RB(3*(J1-1)+1)+BOXLX*NINT(RB(3*(J1-1)+1)/BOXLX) 
412:    !      RB(3*(J1-1)+2)=RB(3*(J1-1)+2)+BOXLY*NINT(RB(3*(J1-1)+2)/BOXLY) 
413:    !      RB(3*(J1-1)+3)=RB(3*(J1-1)+3)+BOXLZ*NINT(RB(3*(J1-1)+3)/BOXLZ) 
414:    !   ENDDO 
415:    ELSE401:    ELSE
416: !402: !
417: !  Translate the RB coordinates to the centre of coordinates of RA.403: !  Translate the RB coordinates to the centre of coordinates of RA.
418: !404: !
419:       DO J1=1,NSIZE405:       DO J1=1,NSIZE
420:          RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-CMXB+CMXA+XSHIFT406:          RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-CMXB+CMXA+XSHIFT
421:          RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-CMYB+CMYA+YSHIFT407:          RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-CMYB+CMYA+YSHIFT
422:          RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT408:          RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT
423:       ENDDO409:       ENDDO
 410:       !write(*,*) "movetest"
 411:       !write(*,*) -CMXB+CMXA+XSHIFT, -CMYB+CMYA+YSHIFT, -CMZB+CMZA+ZSHIFT
424: 412: 
425:       !IF (BULKT) THEN413:       IF (BULKT) THEN
426:       !   BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3414:          DO J1=1,NSIZE
427:       !   DO J1=1,NSIZE415:             RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-BOXLX*NINT(RB(3*(J1-1)+1)/BOXLX)
428:       !      RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-BOXLX*NINT(RB(3*(J1-1)+1)/BOXLX)416:             RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-BOXLY*NINT(RB(3*(J1-1)+2)/BOXLY)
429:       !      RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-BOXLY*NINT(RB(3*(J1-1)+2)/BOXLY)417:             RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-BOXLZ*NINT(RB(3*(J1-1)+3)/BOXLZ)
430:       !      IF (.NOT.TWOD) RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-BOXLZ*NINT(RB(3*(J1-1)+3)/BOXLZ)418:          ENDDO
431:       !   ENDDO419:       ENDIF
432:       !ENDIF 
433:  
434: 420: 
435:       IF (.NOT.BULKT) THEN421:       IF (.NOT.BULKT) THEN
436:          IF (STOCKT) THEN422:          IF (STOCKT) THEN
437:             CALL NEWROTGEOMSTOCK(NATOMS,RB,RMAT,CMXA,CMYA,CMZA)423:             CALL NEWROTGEOMSTOCK(NATOMS,RB,RMAT,CMXA,CMYA,CMZA)
438:          ELSE424:          ELSE
439:             CALL NEWROTGEOM(NSIZE,RB,RMAT,CMXA,CMYA,CMZA)425:             CALL NEWROTGEOM(NSIZE,RB,RMAT,CMXA,CMYA,CMZA)
440:          ENDIF426:          ENDIF
441:       ENDIF427:       ENDIF
442:    ENDIF428:    ENDIF
443: ENDIF429: ENDIF


r31037/OPTIM.F 2016-08-18 15:30:08.491989890 +0100 r31036/OPTIM.F 2016-08-18 15:30:11.256026093 +0100
124:       CALL CPU_TIME(TSTART)124:       CALL CPU_TIME(TSTART)
125: C     IF (CONNECTT.AND.NEWNEBT) THEN125: C     IF (CONNECTT.AND.NEWNEBT) THEN
126: C        PRINT*,'WARNING - cannot use old connect with new neb, changing to old neb'126: C        PRINT*,'WARNING - cannot use old connect with new neb, changing to old neb'
127: C        NEWNEBT=.FALSE.127: C        NEWNEBT=.FALSE.
128: C        NEBT=.TRUE.128: C        NEBT=.TRUE.
129: C     ENDIF129: C     ENDIF
130:       IF ((FILTH2.EQ.0).AND.(FILTH.NE.0)) WRITE(FILTHSTR,'(I10)') FILTH ! otherwise FILTHSTR isn;t set correctly.130:       IF ((FILTH2.EQ.0).AND.(FILTH.NE.0)) WRITE(FILTHSTR,'(I10)') FILTH ! otherwise FILTHSTR isn;t set correctly.
131:       IF (REPELTST) ALLOCATE(REPELTS(NOPT,100)) ! PREVIOUS TS GEOMETRIES TO AVOID131:       IF (REPELTST) ALLOCATE(REPELTS(NOPT,100)) ! PREVIOUS TS GEOMETRIES TO AVOID
132:       IF (CHECKINDEX) ALLOCATE(VECCHK(NOPT,MAX(NUSEEV,HINDEX,1))) ! vectors to orthogonise to132:       IF (CHECKINDEX) ALLOCATE(VECCHK(NOPT,MAX(NUSEEV,HINDEX,1))) ! vectors to orthogonise to
133:       ALLOCATE(ZWORK(NOPT,MAX(NUSEEV,HINDEX,1)))                  ! partial eigenvectors storage133:       ALLOCATE(ZWORK(NOPT,MAX(NUSEEV,HINDEX,1)))                  ! partial eigenvectors storage
134:       IF (TWOENDS.OR.CONNECTT.OR.NEBT.OR.NEWNEBT.OR.DRAGT.OR.GUESSPATHT.OR.MECCANOT.OR.MORPHT.OR.GREATCIRCLET.OR.BHINTERPT134:       IF (TWOENDS.OR.CONNECTT.OR.NEWNEBT.OR.DRAGT.OR.GUESSPATHT.OR.MECCANOT.OR.MORPHT.OR.GREATCIRCLET.OR.BHINTERPT.OR.BISECTT 
135:      & .OR.BISECTT.OR.CLOSESTALIGNMENT.OR.GROWSTRINGT.OR.CLIMBERT) THEN135:      & .OR.CLOSESTALIGNMENT.OR.GROWSTRINGT.OR.CLIMBERT) THEN
136:          ALLOCATE(FIN(NOPT),START(NOPT))136:          ALLOCATE(FIN(NOPT),START(NOPT))
137:       ENDIF137:       ENDIF
138: 138: 
139:       NPCALL=0139:       NPCALL=0
140:       ECALL=0140:       ECALL=0
141:       FCALL=0141:       FCALL=0
142:       SCALL=0142:       SCALL=0
143:       ETIME=0143:       ETIME=0
144:       FTIME=0144:       FTIME=0
145:       STIME=0145:       STIME=0
220:       IF ((ZSYMSAVE.EQ.'QI').AND.(PARAM1.NE.0.0D0)) BULKT=.TRUE.220:       IF ((ZSYMSAVE.EQ.'QI').AND.(PARAM1.NE.0.0D0)) BULKT=.TRUE.
221:       IF ((ZSYMSAVE.EQ.'RM').AND.(PARAM1.NE.0.0D0)) BULKT=.TRUE.221:       IF ((ZSYMSAVE.EQ.'RM').AND.(PARAM1.NE.0.0D0)) BULKT=.TRUE.
222: C     LDUMMY=(NATOMS.EQ.64).AND.(ZSYMSAVE(1:1).EQ.'W') 222: C     LDUMMY=(NATOMS.EQ.64).AND.(ZSYMSAVE(1:1).EQ.'W') 
223:       BULKT=(BULKT.OR.223:       BULKT=(BULKT.OR.
224:      1      (ZSYMSAVE.EQ.'ME').OR.(ZSYMSAVE.EQ.'P6').OR.224:      1      (ZSYMSAVE.EQ.'ME').OR.(ZSYMSAVE.EQ.'P6').OR.
225:      2      (ZSYMSAVE.EQ.'SC').OR.(ZSYMSAVE.EQ.'MS').OR.225:      2      (ZSYMSAVE.EQ.'SC').OR.(ZSYMSAVE.EQ.'MS').OR.
226:      3      (ZSYMSAVE.EQ.'MP').OR.(ZSYMSAVE.EQ.'JM').OR.226:      3      (ZSYMSAVE.EQ.'MP').OR.(ZSYMSAVE.EQ.'JM').OR.
227:      3      (ZSYMSAVE.EQ.'DS').OR.227:      3      (ZSYMSAVE.EQ.'DS').OR.
228:      4      (ZSYMSAVE.EQ.'LP').OR.(ZSYMSAVE.EQ.'LS').OR.(ZSYMSAVE.EQ.'LC').OR.228:      4      (ZSYMSAVE.EQ.'LP').OR.(ZSYMSAVE.EQ.'LS').OR.(ZSYMSAVE.EQ.'LC').OR.
229:      5      (ZSYMSAVE.EQ.'LK'))229:      5      (ZSYMSAVE.EQ.'LK'))
230:  
231:       IF (BULKT.AND.NEWCONNECTT) THEN 
232:          CALL INITIALBOX(Q) 
233:          CALL INITIALBOX(FIN) 
234:       ENDIF 
235:  
236: C     NOSHIFT=(FIELDT.OR.BFGSMINT.OR.BSMIN.OR.RKMIN.OR.NOHESS.OR.BFGSSTEP)230: C     NOSHIFT=(FIELDT.OR.BFGSMINT.OR.BSMIN.OR.RKMIN.OR.NOHESS.OR.BFGSSTEP)
237: C     IF (INR.GE.0) NOSHIFT=.FALSE. ! changed INR default value in keywords to -1231: C     IF (INR.GE.0) NOSHIFT=.FALSE. ! changed INR default value in keywords to -1
238: C                                     This isn;t good enough - if we are doing a path and232: C                                     This isn;t good enough - if we are doing a path and
239: C                                     EFOL is called once to get the eigenvalues and eigenvectors233: C                                     EFOL is called once to get the eigenvalues and eigenvectors
240: C                                     then we need to SHIFT.234: C                                     then we need to SHIFT.
241: C  Shifting only needs to be turned off if we are doing VARIABLES or BFGSTS without NOIT;235: C  Shifting only needs to be turned off if we are doing VARIABLES or BFGSTS without NOIT;
242: C  account for this in BFGSTS.236: C  account for this in BFGSTS.
243: C     IF (NOIT) NOSHIFT=.FALSE.237: C     IF (NOIT) NOSHIFT=.FALSE.
244: C     IF (VARIABLES) NOSHIFT=.TRUE.238: C     IF (VARIABLES) NOSHIFT=.TRUE.
245: C     NOSHIFT=(FIELDT.OR.BFGSMINT.OR.BSMIN.OR.RKMIN.OR.NOHESS.OR.BFGSSTEP)239: C     NOSHIFT=(FIELDT.OR.BFGSMINT.OR.BSMIN.OR.RKMIN.OR.NOHESS.OR.BFGSSTEP)
554:             IF (RIGIDINIT) THEN548:             IF (RIGIDINIT) THEN
555: ! hk286               CALL CHARMMDUMP(Q,"test1.pdb",.FALSE.)549: ! hk286               CALL CHARMMDUMP(Q,"test1.pdb",.FALSE.)
556:                ATOMRIGIDCOORDT = .FALSE.550:                ATOMRIGIDCOORDT = .FALSE.
557:                CALL TRANSFORMCTORIGID (Q, XRIGIDCOORDS)551:                CALL TRANSFORMCTORIGID (Q, XRIGIDCOORDS)
558:                CALL POTENTIAL(XRIGIDCOORDS,EINITIAL,LGDUMMY,.TRUE.,.FALSE.,RMSINITIAL,.FALSE.,.FALSE.)552:                CALL POTENTIAL(XRIGIDCOORDS,EINITIAL,LGDUMMY,.TRUE.,.FALSE.,RMSINITIAL,.FALSE.,.FALSE.)
559:                ATOMRIGIDCOORDT = .TRUE.553:                ATOMRIGIDCOORDT = .TRUE.
560:             ELSE554:             ELSE
561:                CALL POTENTIAL(Q,EINITIAL,LGDUMMY,.TRUE.,.FALSE.,RMSINITIAL,.FALSE.,.FALSE.)555:                CALL POTENTIAL(Q,EINITIAL,LGDUMMY,.TRUE.,.FALSE.,RMSINITIAL,.FALSE.,.FALSE.)
562:             ENDIF556:             ENDIF
563:             WRITE(*,'(a,2(g20.10,a))') ' OPTIM> Initial energy=',EINITIAL,' RMS force=',RMSinitial557:             WRITE(*,'(a,2(g20.10,a))') ' OPTIM> Initial energy=',EINITIAL,' RMS force=',RMSinitial
564:             IF (UNRST) THEN558:            IF (UNRST) THEN
565:                 DO J1=1,NRES559:                 DO J1=1,NRES
566:                    C(1,J1)=FIN(6*(J1-1)+1)560:                    C(1,J1)=FIN(6*(J1-1)+1)
567:                    C(2,J1)=FIN(6*(J1-1)+2)561:                    C(2,J1)=FIN(6*(J1-1)+2)
568:                    C(3,J1)=FIN(6*(J1-1)+3)562:                    C(3,J1)=FIN(6*(J1-1)+3)
569:                    C(1,J1+NRES)=FIN(6*(J1-1)+4)563:                    C(1,J1+NRES)=FIN(6*(J1-1)+4)
570:                    C(2,J1+NRES)=FIN(6*(J1-1)+5)564:                    C(2,J1+NRES)=FIN(6*(J1-1)+5)
571:                    C(3,J1+NRES)=FIN(6*(J1-1)+6)565:                    C(3,J1+NRES)=FIN(6*(J1-1)+6)
572:                 ENDDO566:                 ENDDO
573:                 CALL UPDATEDC567:                 CALL UPDATEDC
574:                 CALL INT_FROM_CART(.TRUE.,.FALSE.)568:                 CALL INT_FROM_CART(.TRUE.,.FALSE.)
1033:      1         ' %=',ETIME*100.0D0/TFINISH1027:      1         ' %=',ETIME*100.0D0/TFINISH
1034:          WRITE(*,'(A,I10,A,F15.2,A,F5.1)') ' OPTIM> # of energy+gradient calls=        ',FCALL,' time=',FTIME,1028:          WRITE(*,'(A,I10,A,F15.2,A,F5.1)') ' OPTIM> # of energy+gradient calls=        ',FCALL,' time=',FTIME,
1035:      1         ' %=',FTIME*100.0D0/TFINISH1029:      1         ' %=',FTIME*100.0D0/TFINISH
1036:          WRITE(*,'(A,I10,A,F15.2,A,F5.1)') ' OPTIM> # of energy+gradient+Hessian calls=',SCALL,' time=',STIME,1030:          WRITE(*,'(A,I10,A,F15.2,A,F5.1)') ' OPTIM> # of energy+gradient+Hessian calls=',SCALL,' time=',STIME,
1037:      1         ' %=',STIME*100/TFINISH1031:      1         ' %=',STIME*100/TFINISH
1038:       ENDIF1032:       ENDIF
1039: 1033: 
1040: 1034: 
1041:       RETURN1035:       RETURN
1042:       END1036:       END
1043:  
1044:       SUBROUTINE INITIALBOX(QC) 
1045:       USE KEY, ONLY: BULK_BOXVEC 
1046:       USE COMMONS, ONLY: NATOMS 
1047:       IMPLICIT NONE 
1048:       DOUBLE PRECISION QC(NATOMS) 
1049:       INTEGER J1 
1050:        
1051:       DO J1=1, NATOMS 
1052:          QC(3*(J1-1)+1)=QC(3*(J1-1)+1)-BULK_BOXVEC(1)*NINT(QC(3*(J1-1)+1)/BULK_BOXVEC(1)) 
1053:          QC(3*(J1-1)+2)=QC(3*(J1-1)+2)-BULK_BOXVEC(2)*NINT(QC(3*(J1-1)+2)/BULK_BOXVEC(2)) 
1054:          QC(3*(J1-1)+3)=QC(3*(J1-1)+3)-BULK_BOXVEC(3)*NINT(QC(3*(J1-1)+3)/BULK_BOXVEC(3)) 
1055:       ENDDO 
1056:  
1057:       END 


r31037/output.f90 2016-08-18 15:30:06.883966942 +0100 r31036/output.f90 2016-08-18 15:30:11.008022738 +0100
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: MODULE NEBOUTPUT 19: MODULE NEBOUTPUT
 20:      IMPLICIT NONE 20:      IMPLICIT NONE
 21:      CONTAINS 21:      CONTAINS
 22:  22: 
 23:      SUBROUTINE TSLOCATOR(TSRESET) 23:      SUBROUTINE TSLOCATOR(TSRESET)
 24:           USE KEY,ONLY: BFGSTST,UNRST,NSTEPS,MACHINE, GROWSTRINGT, INTEPSILON, REDOTSIM, EDIFFTOL, & 24:           USE KEY,ONLY: BFGSTST,UNRST,NSTEPS,MACHINE, GROWSTRINGT, INTEPSILON, REDOTSIM, EDIFFTOL, &
 25:   &                     NONEBMAX, INTCONSTRAINTT, CUDAT, FLATPATHT, FLATTESTT, FLATEDIFF 25:   &                     NONEBMAX, INTCONSTRAINTT, CUDAT
 26:           USE GSDATA, ONLY: EVOLVESTRINGT 26:           USE GSDATA, ONLY: EVOLVESTRINGT
 27:           USE KEYOUTPUT 27:           USE KEYOUTPUT
 28:           USE MODCHARMM 28:           USE MODCHARMM
 29:           USE NEBDATA 29:           USE NEBDATA
 30:           USE KEYNEB,ONLY:NIMAGE,DEBUG 30:           USE KEYNEB,ONLY:NIMAGE,DEBUG
 31:           USE NEBTOCONNECT 31:           USE NEBTOCONNECT
 32:           USE CHARUTILS 32:           USE CHARUTILS
 33:           USE MODGUESS 33:           USE MODGUESS
 34:           USE MODMEC 34:           USE MODMEC
 35:           USE LINKEDLIST 35:           USE LINKEDLIST
 36:           USE MODEFOL 36:           USE MODEFOL
 37:           USE INTCOMMONS, ONLY : DESMINT, NINTC, NNZ, KD, INTNEWT 37:           USE INTCOMMONS, ONLY : DESMINT, NINTC, NNZ, KD, INTNEWT
 38:           USE COMMONS, ONLY : REDOPATH, REDOPATHNEB 38:           USE COMMONS, ONLY : REDOPATH, REDOPATHNEB
 39:           USE MODCUDABFGSTS, ONLY : CUDA_BFGSTS_WRAPPER 39:           USE MODCUDABFGSTS, ONLY : CUDA_BFGSTS_WRAPPER
 40:           USE GENRIGID, ONLY: DEGFREEDOMS, RIGIDINIT 40:           USE GENRIGID, ONLY: DEGFREEDOMS, RIGIDINIT
 41:  41: 
 42:           IMPLICIT NONE 42:           IMPLICIT NONE
 43:            43:           
 44:           INTEGER :: I,J,NT,NM,ITDONE=0,J1,RECLEN 44:           INTEGER :: I,J,NT,ITDONE=0,J1,RECLEN
 45:           INTEGER,PARAMETER :: MAXPRINTOUT = 50, ITMAX  = 30 45:           INTEGER,PARAMETER :: MAXPRINTOUT = 50, ITMAX  = 30
 46:           DOUBLE PRECISION :: EDUMMY,EVALMIN,EVALMAX,MAXE,VECSNORM 46:           DOUBLE PRECISION :: EDUMMY,EVALMIN,EVALMAX,MAXE,VECSNORM
 47:           LOGICAL :: TSCONVERGED,T,TSRESET 47:           LOGICAL :: TSCONVERGED,T,TSRESET
 48:           DOUBLE PRECISION,DIMENSION(NOPT) :: LGDUMMY, VECS, DIAG 48:           DOUBLE PRECISION,DIMENSION(NOPT) :: LGDUMMY, VECS, DIAG
 49:           INTEGER :: MLOC 49:           INTEGER :: MLOC
 50:           DOUBLE PRECISION :: TIME, TIME0 50:           DOUBLE PRECISION :: TIME, TIME0
 51:           DOUBLE PRECISION :: DPRAND 51:           DOUBLE PRECISION :: DPRAND
 52:           LOGICAL :: KNOWE, KNOWG, KNOWH ! JMC 52:           LOGICAL :: KNOWE, KNOWG, KNOWH ! JMC
 53:           COMMON /KNOWN/ KNOWE, KNOWG, KNOWH ! JMC 53:           COMMON /KNOWN/ KNOWE, KNOWG, KNOWH ! JMC
 54:           CHARACTER(LEN=256) :: FILENAME, METHSTR 54:           CHARACTER(LEN=256) :: FILENAME, METHSTR
 55:           INTEGER TSPOS(NIMAGE+2)!, MINPOS(NIMAGE+2) 55:           INTEGER TSPOS(NIMAGE+2)
 56:           INTEGER BIGN 56: 
 57:           !DOUBLE PRECISION :: SCA 
 58:           LOGICAL :: TMPINTNEWT, FAILED 57:           LOGICAL :: TMPINTNEWT, FAILED
 59:  58: 
 60:           NT = 0 59:           NT = 0
 61:           VECS(:) = 0 ! sn402: to avoid uninitialised value problems 60:           VECS(:) = 0 ! sn402: to avoid uninitialised value problems
 62:  61: 
 63:           IF (REDOPATHNEB) THEN 62:           IF (REDOPATHNEB) THEN
 64:              NT=1 63:              NT=1
 65:              MAXE=-1.0D100 64:              MAXE=-1.0D100
 66:              MLOC=REDOTSIM+1 65:              MLOC=REDOTSIM+1
 67:              PRINT '(A,F20.10)',' tslocator> transition state has energy ',EEE(REDOTSIM+1) 66:              PRINT '(A,F20.10)',' tslocator> transition state has energy ',EEE(REDOTSIM+1)
 75:                   MAXE=-1.0D100 74:                   MAXE=-1.0D100
 76:                   DO J1=2,NIMAGE+1 75:                   DO J1=2,NIMAGE+1
 77:                      IF (EEE(J1).GT.MAXE) THEN 76:                      IF (EEE(J1).GT.MAXE) THEN
 78:                         MLOC=J1 77:                         MLOC=J1
 79:                         MAXE=EEE(J1) 78:                         MAXE=EEE(J1)
 80:                      ENDIF 79:                      ENDIF
 81:                   ENDDO 80:                   ENDDO
 82:                   TSPOS(1)=MLOC 81:                   TSPOS(1)=MLOC
 83:              CASE('all','maxim') 82:              CASE('all','maxim')
 84:                   DO I=2,NIMAGE+1  83:                   DO I=2,NIMAGE+1 
 85:                        T=.FALSE. 84:                        T=.TRUE.
 86:                        IF (CANDIDATES=='maxim') then 85:                        IF (CANDIDATES=='maxim') then
 87:                             IF ( EEE(I-1)+EDIFFTOL*10.0D0 < EEE(I) .AND. EEE(I) > EEE(I+1)+EDIFFTOL*10.0D0 ) THEN 86:                             IF ( EEE(I-1)+EDIFFTOL*10.0D0 < EEE(I) .AND. EEE(I) > EEE(I+1)+EDIFFTOL*10.0D0 ) THEN
 88:                                  T=.TRUE. 87:                                  T=.TRUE.
 89:                             ELSE 88:                             ELSE
 90:                                  T=.FALSE. 89:                                  T=.FALSE.
 91:                             ENDIF 90:                             ENDIF
 92:                        ENDIF 91:                        ENDIF
 93: !                      PRINT '(A,I6,3F20.10,L5)','I,EEE(I-1),EEE(I),EEE(I+1),T=',I,EEE(I-1),EEE(I),EEE(I+1),T 92: !                      PRINT '(A,I6,3F20.10,L5)','I,EEE(I-1),EEE(I),EEE(I+1),T=',I,EEE(I-1),EEE(I),EEE(I+1),T
 94:                        IF (T) THEN 93:                        IF (T) THEN
 95:                             NT=NT+1 94:                             NT=NT+1
108:              PRINT '(1x,a)', 'No maximum in profile - using highest image'107:              PRINT '(1x,a)', 'No maximum in profile - using highest image'
109:              NT=1108:              NT=1
110:              IF (EEE(2).GT.EEE(NIMAGE+1)) THEN109:              IF (EEE(2).GT.EEE(NIMAGE+1)) THEN
111:                 TSPOS(1)=2110:                 TSPOS(1)=2
112:              ELSE111:              ELSE
113:                 TSPOS(1)=NIMAGE+1112:                 TSPOS(1)=NIMAGE+1
114:              ENDIF113:              ENDIF
115:           ENDIF114:           ENDIF
116: 115: 
117:           WRITE(*,'(1X,A,I4,A)',advance='No') 'Following ',NT,' images are candidates for TS:'116:           WRITE(*,'(1X,A,I4,A)',advance='No') 'Following ',NT,' images are candidates for TS:'
118:           !open(UNIT=323, FILE="energydiff") 
119:           DO J=1,NT117:           DO J=1,NT
120:              WRITE(*,'(i5)',advance='No') TSPOS(J)-1118:              WRITE(*,'(i5)',advance='No') TSPOS(J)-1
121:           !   write(323,*) EEE(TSPOS(J))-EEE(1), EEE(TSPOS(J))-EEE(NIMAGE+2) 
122:           ENDDO119:           ENDDO
123:           PRINT *,' '120:           PRINT *,' '
124: 121:           
125: !122:           IF (OPTIMIZETS) THEN
126: !sy349: This block is to test whether it is a flat path after the dneb calculation. If 
127: !the energy difference between highest-energy image and two end points is 
128: !smaller than FLATEDIFF. If it is a flat path, it will regard the highest-energy image directly 
129: !as transition state, and these two end points will be regarded connected. 
130: ! 
131:           IF (FLATTESTT) THEN 
132:  
133:              FLATPATHT=.FALSE. 
134:              BIGN=TSPOS(1) 
135:              DO J=2, NT                 
136:                 IF (EEE(BIGN)<EEE(TSPOS(J))) BIGN=TSPOS(J) 
137:              ENDDO 
138:  
139:              IF (ABS(EEE(1)-EEE(NIMAGE+2)).LT.EDIFFTOL) THEN 
140:                 IF (ABS(EEE(BIGN)-EEE(1)).LT.FLATEDIFF & 
141:                    &.OR.ABS(EEE(BIGN)-EEE(NIMAGE+2)).LT.FLATEDIFF) THEN 
142:                    FLATPATHT=.TRUE. 
143:                    IF (TSRESET) NTSFOUND=0 
144:                    NTSFOUND=NTSFOUND+1 
145:                    ALLOCATE(TSFOUND(NTSFOUND)%E,TSFOUND(NTSFOUND)%COORD(NOPT),& 
146:                            &TSFOUND(NTSFOUND)%EVALMIN,TSFOUND(NTSFOUND)%VECS(NOPT)) 
147:                    TSFOUND(NTSFOUND)%VECS=VECS(1:NOPT) 
148:                    TSFOUND(NTSFOUND)%COORD=XYZ(NOPT*(BIGN-1)+1:NOPT*BIGN) 
149:                    TSFOUND(NTSFOUND)%E=EEE(BIGN) 
150:                    TSFOUND(NTSFOUND)%EVALMIN=EVALMIN 
151:                 ENDIF 
152:              ENDIF 
153:  
154:              IF (DEBUG .AND. FLATPATHT) THEN 
155:                  write(*,*) "Index of first candidate is ", TSPOS(1)  ! sn402 
156:                  write(*,*) "Image energies are", EEE(:) !sn402 
157:              ENDIF 
158:           ENDIF 
159:  
160:           IF (OPTIMIZETS.AND. .NOT.FLATPATHT) THEN 
161:              IF (DEBUG) THEN123:              IF (DEBUG) THEN
162:                  write(*,*) "Index of first candidate is ", TSPOS(1)  ! sn402124:                  write(*,*) "Index of first candidate is ", TSPOS(1)  ! sn402
163:                  write(*,*) "Image energies are", EEE(:) !sn402125:                  write(*,*) "Image energies are", EEE(:) !sn402
164:              ENDIF126:              ENDIF
165:              IF (TSRESET) NTSFOUND=0127:              IF (TSRESET) NTSFOUND=0
166:              CALL MYCPU_TIME(STARTTIME,.FALSE.)128:              CALL MYCPU_TIME(STARTTIME,.FALSE.)
167:              DO J=1,NT129:              DO J=1,NT
168:                 !IF (FLATPATHT(J)) CYCLE 
169:                 CALL MYCPU_TIME(TIME0,.FALSE.)130:                 CALL MYCPU_TIME(TIME0,.FALSE.)
170:                 EDUMMY=EEE(TSPOS(J))131:                 EDUMMY=EEE(TSPOS(J))
171:                 LGDUMMY(1:NOPT)=TRUEGRAD((TSPOS(J)-1)*NOPT+1:TSPOS(J)*NOPT)132:                 LGDUMMY(1:NOPT)=TRUEGRAD((TSPOS(J)-1)*NOPT+1:TSPOS(J)*NOPT)
172:                 KNOWE=.TRUE.133:                 KNOWE=.TRUE.
173:                 KNOWG=.TRUE.134:                 KNOWG=.TRUE.
174:                 IF (REDOPATH) THEN135:                 IF (REDOPATH) THEN
175:                    KNOWG = .FALSE.136:                    KNOWG = .FALSE.
176:                    KNOWE = .FALSE.137:                    KNOWE = .FALSE.
177:                 ENDIF138:                 ENDIF
178:                 IF (BFGSTST) THEN139:                 IF (BFGSTST) THEN
287:                    WRITE(METHSTR,'(A)') 'ES'248:                    WRITE(METHSTR,'(A)') 'ES'
288:                 ELSE249:                 ELSE
289:                    WRITE(METHSTR,'(A)') 'GS'250:                    WRITE(METHSTR,'(A)') 'GS'
290:                 ENDIF251:                 ENDIF
291:              ELSE252:              ELSE
292:                 WRITE(METHSTR,'(A)') 'DNEB'253:                 WRITE(METHSTR,'(A)') 'DNEB'
293:              ENDIF              254:              ENDIF              
294: 255: 
295:              WRITE(*, '(1x,a,f10.2)',advance='yes') trim(METHSTR)//' run yielded '//trim(adjustl(IntStr))// &256:              WRITE(*, '(1x,a,f10.2)',advance='yes') trim(METHSTR)//' run yielded '//trim(adjustl(IntStr))// &
296:                           &' true transition state(s) time=',EndTime-StartTime257:                           &' true transition state(s) time=',EndTime-StartTime
297:           ENDIF258:         ENDIF
298:           259:         IF (SAVECANDIDATES) THEN
299:           IF (SAVECANDIDATES) THEN260:            DO J=1,NTSFOUND
300:              DO J=1,NTSFOUND261:               IF (DESMINT) THEN
301:                 IF (DESMINT) THEN262:                  INQUIRE(IOLENGTH=RECLEN) (XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)))
302:                    INQUIRE(IOLENGTH=RECLEN) (XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)))263:               ELSE                       
303:                 ELSE                       264:                  INQUIRE(IOLENGTH=RECLEN) (XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)))
304:                    INQUIRE(IOLENGTH=RECLEN) (XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)))265:               ENDIF
305:                 ENDIF266:               WRITE(FILENAME,'(i10)') J
306:                 WRITE(FILENAME,'(i10)') J267:               FILENAME='points'//trim(adjustl(filename))//'.out'
307:                 FILENAME='points'//trim(adjustl(filename))//'.out'268:               OPEN(UNIT=40,FILE=FILENAME,STATUS='unknown',form='unformatted',access='direct',recl=reclen)
308:                 OPEN(UNIT=40,FILE=FILENAME,STATUS='unknown',form='unformatted',access='direct',recl=reclen)269: 
309: 270:               IF (DESMINT) THEN
310:                 IF (DESMINT) THEN271:                  WRITE(40,REC=1) ( XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)) )
311:                    WRITE(40,REC=1) ( XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)) )272:               ELSE
312:                 ELSE273:                  WRITE(40,REC=1) ( XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)) )
313:                    WRITE(40,REC=1) ( XYZ(NOPT*(TSPOS(J)-1)+1:NOPT*TSPOS(J)) )274:               ENDIF
314:                 ENDIF275: 
315: 276:               CLOSE(40)
316:                 CLOSE(40)277:            ENDDO
317:              ENDDO278:         ENDIF
318:           ENDIF 
319: 279: 
320:           RETURN280:         RETURN
321: 281: 
322:       END SUBROUTINE TSLOCATOR282:       END SUBROUTINE TSLOCATOR
323: 283: 
324: SUBROUTINE CONTSLOCATOR284: SUBROUTINE CONTSLOCATOR
325: USE KEY,ONLY: BFGSTST,UNRST,NSTEPS,MACHINE, GROWSTRINGT, INTEPSILON, REDOTSIM285: USE KEY,ONLY: BFGSTST,UNRST,NSTEPS,MACHINE, GROWSTRINGT, INTEPSILON, REDOTSIM
326: USE GSDATA, ONLY: EVOLVESTRINGT286: USE GSDATA, ONLY: EVOLVESTRINGT
327: USE KEYOUTPUT287: USE KEYOUTPUT
328: USE MODCHARMM288: USE MODCHARMM
329: USE NEBDATA289: USE NEBDATA
330: USE KEYNEB,ONLY:NIMAGE,DEBUG290: USE KEYNEB,ONLY:NIMAGE,DEBUG


r31037/stealthy.f90 2016-08-18 15:30:09.996009045 +0100 r31036/stealthy.f90 2016-08-18 15:30:13.140051583 +0100
  1: SUBROUTINE STEALTHY(XS,TIDU,ENERGY,GTEST,STEST)  1: SUBROUTINE STEALTHY(XS,TIDU,ENERGY,GTEST,STEST)
  2:     USE COMMONS  2:     USE COMMONS
  3:     USE KEY  3:     USE KEY
  4:     USE MODHESS  4:     USE MODHESS
  5:     IMPLICIT NONE  5:     IMPLICIT NONE
  6:     LOGICAL GTEST,STEST  6:     LOGICAL GTEST,STEST
  7:     INTEGER M, M1, M2, M3, J1, J2, I1, I2, L1, L2, KI, I  7:     INTEGER M, M1, M2, M3, J1, J2, I1, I2, L1, L2, KI
  8:     DOUBLE PRECISION XS(3*NATOMS), KX(3), KY(3), KZ(3), KN(3), TIDU(3*NATOMS), &  8:     DOUBLE PRECISION XS(3*NATOMS), KX(3), KY(3), KZ(3), KN(3), TIDU(3*NATOMS), &
  9:                      CKR(NATOMS), SKR(NATOMS)  9:                      CKR(NATOMS), SKR(NATOMS)
 10:     DOUBLE PRECISION PI, KR ,C, S, ENERGY, VF, IM, KM, FORPI, STRMS 10:     DOUBLE PRECISION PI, KR ,C, S, ENERGY, VF, IM, KM, FORPI
 11:  11: 
 12:     IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS, 3*NATOMS)) 12:     IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
 13:  13: 
 14: !   initialize 14: !   initialize
 15:     FORPI=1.0 15:     FORPI=1.0
 16:     PI=4*DATAN(FORPI) 16:     PI=4*DATAN(FORPI)
  17:     !PI=3.141592653589793
 17:     KM=KLIM*KLIM 18:     KM=KLIM*KLIM
 18:     KI=INT(KLIM) 19:     KI=INT(KLIM)
 19:     ENERGY=0 20:     ENERGY=0
 20:     STM=0 21:     STM=0
 21:  22: 
 22:     IF (GTEST) THEN 23:     IF (GTEST) THEN
 23:        TIDU=0 24:        TIDU=0
 24:     END IF 25:     END IF
 25:  26: 
 26:     IF (STEST) THEN 27:     IF (STEST) THEN
 27:        HESS=0 28:        HESS=0
 28:     END IF 29:     END IF
 29:  30: 
  31: !   get the initial box
  32: 
 30: !   calculate the wave vector 33: !   calculate the wave vector
 31:     KX(1)=2*PI/BULK_BOXVEC(1) 34:     KX(1)=2*PI/BULK_BOXVEC(1)
 32:     KX(2)=0 35:     KX(2)=0
 33:     KX(3)=0 36:     KX(3)=0
 34:     KY(1)=0 37:     KY(1)=0
 35:     KY(2)=2*PI/BULK_BOXVEC(2) 38:     KY(2)=2*PI/BULK_BOXVEC(2)
 36:     KY(3)=0 39:     KY(3)=0
 37:     KZ(1)=0 40:     KZ(1)=0
 38:     KZ(2)=0 41:     KZ(2)=0
 39:     KZ(3)=2*PI/BULK_BOXVEC(3) 42:     KZ(3)=2*PI/BULK_BOXVEC(3)
 40:     VF=BULK_BOXVEC(1)*BULK_BOXVEC(2)*BULK_BOXVEC(3) 43:     VF=BULK_BOXVEC(1)*BULK_BOXVEC(2)*BULK_BOXVEC(3)
 41: !   WRITE(*,*) KI(1), KI(2), KI(3) 44: !   WRITE(*,*) KI(1), KI(2), KI(3)
 42:  45: 
 43: !   periodic boundary condition 46: !   periodic boundary condition
 44:     IF (BULKT.AND. .NOT. BULKBOXT) THEN 47:     IF (BULKT) THEN
 45:        DO J1=1, NATOMS 48:        DO J1=1, NATOMS
 46:           J2=3*(J1-1) 49:           J2=3*(J1-1)
 47:           XS(J2+1)=XS(J2+1) - BULK_BOXVEC(1)*DNINT(XS(J2+1)/BULK_BOXVEC(1)) 50:           XS(J2+1)=XS(J2+1) - BULK_BOXVEC(1)*DNINT(XS(J2+1)/BULK_BOXVEC(1))
 48:           XS(J2+2)=XS(J2+2) - BULK_BOXVEC(2)*DNINT(XS(J2+2)/BULK_BOXVEC(2)) 51:           XS(J2+2)=XS(J2+2) - BULK_BOXVEC(2)*DNINT(XS(J2+2)/BULK_BOXVEC(2))
 49:           XS(J2+3)=XS(J2+3) - BULK_BOXVEC(3)*DNINT(XS(J2+3)/BULK_BOXVEC(3)) 52:           XS(J2+3)=XS(J2+3) - BULK_BOXVEC(3)*DNINT(XS(J2+3)/BULK_BOXVEC(3))
 50:        END DO 53:        END DO
 51:     END IF 54:     END IF
 52: !    write(*,*) XS 55: !    write(*,*) XS
 53:   56:  
 54: !   sum of k 57: !   sum of k
 67:              KN(1)=KX(1)*M1 + KY(1)*M2 + KZ(1)*M3 70:              KN(1)=KX(1)*M1 + KY(1)*M2 + KZ(1)*M3
 68:              KN(2)=KX(2)*M1 + KY(2)*M2 + KZ(2)*M3 71:              KN(2)=KX(2)*M1 + KY(2)*M2 + KZ(2)*M3
 69:              KN(3)=KX(3)*M1 + KY(3)*M2 + KZ(3)*M3 72:              KN(3)=KX(3)*M1 + KY(3)*M2 + KZ(3)*M3
 70:  73: 
 71: !            calculation of Re[n(k)] and Im[n(k)] 74: !            calculation of Re[n(k)] and Im[n(k)]
 72:              DO I1=1, NATOMS 75:              DO I1=1, NATOMS
 73:                 I2=3*(I1-1) 76:                 I2=3*(I1-1)
 74:                 KR=KN(1)*XS(I2+1) + KN(2)*XS(I2+2) + KN(3)*XS(I2+3) 77:                 KR=KN(1)*XS(I2+1) + KN(2)*XS(I2+2) + KN(3)*XS(I2+3)
 75:                 CKR(I1)=COS(KR) 78:                 CKR(I1)=COS(KR)
 76:                 SKR(I1)=SIN(KR) 79:                 SKR(I1)=SIN(KR)
 77:                 !if (M2==5) write(*,'(T1,4(I,1X),F,1X,F)') M1, M2, M3, I1, CKR(I1), SKR(I1) 80:                 !if (M1==0.AND.M2==-3) write(*,'(T1,4(I,1X),F,1X,F)') M1, M2, M3, I1, CKR(I1), SKR(I1)
 78:  81: 
 79:                 C=C+CKR(I1) 82:                 C=C+CKR(I1)
 80:                 S=S+SKR(I1) 83:                 S=S+SKR(I1)
 81: !                if (M1==0.AND.M2==-3) write(*,'(T1,F,1X,F)') C, S 84: !                if (M1==0.AND.M2==-3) write(*,'(T1,F,1X,F)') C, S
 82:              END DO 85:              END DO
 83:              !write(*,*) M1, M2, M3 86:              !write(*,*) M1, M2, M3
 84:              !write(*,*) CKR 87:              !write(*,*) CKR(1:6)
 85:              !write(*,*) SKR 88:              !write(*,*) SKR(1:6)
 86:              !if (M1==0.AND.M2==-8) then 89:              !if (M1==0.AND.M2==-8) then
 87:              !   write(*,*) CKR(1:6) 90:              !   write(*,*) CKR(1:6)
 88:              !   write(*,*) SKR(1:6) 91:              !   write(*,*) SKR(1:6)
 89:              !end if 92:              !end if
 90:              !write(*,'(T1,3(I,1X),F,1X,F)') M1, M2, M3, C, S 93: !             write(*,'(T1,3(I,1X),F,1X,F)') M1, M2, M3, C, S
 91:  94: 
 92: !            E=n(k)^2/Vf 95: !            E=n(k)^2/Vf
 93:              ENERGY=ENERGY+ SCA*(C*C+S*S)/VF 96:              ENERGY=ENERGY+ SCA*(C*C+S*S)/VF
 94:  97: 
 95: !            gradient 98: !            gradient
 96:              IF (GTEST) THEN 99:              IF (GTEST) THEN
 97:                 DO L1=1, NATOMS100:                 DO L1=1, NATOMS
 98:                    L2=3*(L1-1)101:                    L2=3*(L1-1)
 99:                    IM=( C*SKR(L1) - S*CKR(L1) )*SCA102:                    IM=( C*SKR(L1) - S*CKR(L1) )*SCA
100:                    TIDU(L2+1)=TIDU(L2+1) - 2*KN(1)*IM/VF103:                    TIDU(L2+1)=TIDU(L2+1) - 2*KN(1)*IM/VF
106: !            hessian matrix109: !            hessian matrix
107:              IF (STEST) THEN110:              IF (STEST) THEN
108:                 CALL STHESS(XS, KN, VF, CKR, C, SKR, S)111:                 CALL STHESS(XS, KN, VF, CKR, C, SKR, S)
109:              END IF112:              END IF
110: 113: 
111:              STM=STM+1114:              STM=STM+1
112: 115: 
113:           END DO116:           END DO
114:        END DO117:        END DO
115:     END DO118:     END DO
116:     !write(*,*) "STM", STM 
117:  
118:     !DO I=1, NATOMS 
119:     !   STRMS=STRMS+TIDU(I)**2 
120:     !ENDDO 
121:  
122:     !STRMS=DSQRT(STRMS/NATOMS) 
123:  
124:     !write(*,*) "STRMS", STRMS 
125: 119: 
126:     IF (STEST.AND.STEALTV) THEN120:     IF (STEST.AND.STEALTV) THEN
127:        CALL STEALTHY_EIGEN_TEST(ENERGY, TIDU)121:        CALL STEALTHY_EIGEN_TEST(ENERGY, TIDU)
128:     ENDIF122:     ENDIF
129: 123: 
 124: !   change the box back to cubic
 125: 
130: END SUBROUTINE126: END SUBROUTINE
131: 127: 
132: SUBROUTINE STHESS(X, K, V, CK, CSUM, SK, SSUM)128: SUBROUTINE STHESS(X, K, V, CK, CSUM, SK, SSUM)
133:    USE COMMONS129:    USE COMMONS
134:    USE KEY130:    USE KEY
135:    USE MODHESS131:    USE MODHESS
136:    IMPLICIT NONE132:    IMPLICIT NONE
137:    INTEGER I1, I2, J1, J2133:    INTEGER I1, I2, J1, J2
138:    DOUBLE PRECISION X(3*NATOMS), K(3), RIJ(3), CK(NATOMS), SK(NATOMS)134:    DOUBLE PRECISION X(3*NATOMS), K(3), RIJ(3), CK(NATOMS), SK(NATOMS)
139:    DOUBLE PRECISION KRIJ, V, CSUM, SSUM, CI, SI, CF, CSKR135:    DOUBLE PRECISION KRIJ, V, CSUM, SSUM, CI, SI, CF, CSKR


r31037/tryconnect.f90 2016-08-18 15:30:06.387960230 +0100 r31036/tryconnect.f90 2016-08-18 15:30:10.504015919 +0100
 25:           USE PORFUNCS 25:           USE PORFUNCS
 26:           USE CONNECTDATA 26:           USE CONNECTDATA
 27:           USE KEYCONNECT 27:           USE KEYCONNECT
 28:           USE CONNECTUTILS 28:           USE CONNECTUTILS
 29:           USE NEBTOCONNECT 29:           USE NEBTOCONNECT
 30:           USE KEYNEB, ONLY : NIMAGE, NITERMAX, READGUESS 30:           USE KEYNEB, ONLY : NIMAGE, NITERMAX, READGUESS
 31:           USE KEY, ONLY : UNRST, FILTH, FILTHSTR, DUMPALLPATHS, TWOD, MAXTSENERGY, RIGIDBODY, & 31:           USE KEY, ONLY : UNRST, FILTH, FILTHSTR, DUMPALLPATHS, TWOD, MAXTSENERGY, RIGIDBODY, &
 32:   &                   MAXMAXBARRIER, DIJKSTRALOCAL, INTNTRIESMAX, & 32:   &                   MAXMAXBARRIER, DIJKSTRALOCAL, INTNTRIESMAX, &
 33:   &                   PERMDIST, MAXBARRIER, GROWSTRINGT, BULKT, AMBERT, NABT, AMBER12T, FREEZE, FROZEN, NFREEZE, & 33:   &                   PERMDIST, MAXBARRIER, GROWSTRINGT, BULKT, AMBERT, NABT, AMBER12T, FREEZE, FROZEN, NFREEZE, &
 34:   &                   PUSHOFF, MAXBFGS, MIN1REDO, MIN2REDO, REDOKADD, D1INIT, D2INIT, REDOE1, REDOE2, & 34:   &                   PUSHOFF, MAXBFGS, MIN1REDO, MIN2REDO, REDOKADD, D1INIT, D2INIT, REDOE1, REDOE2, &
 35:   &                   AMHT, SEQ, INTIMAGE, NEBKFINAL, INTIMAGEINCR, MAXINTIMAGE, SETCHIRAL, REDOTS, VARIABLES, & 35:   &                   AMHT, SEQ, INTIMAGE, NEBKFINAL, INTIMAGEINCR, MAXINTIMAGE, SETCHIRAL, REDOTS, VARIABLES, QCIIMAGE
 36:   &                   FLATPATHT, QCIIMAGE 
 37:           USE MODGUESS 36:           USE MODGUESS
 38:           USE MODUNRES 37:           USE MODUNRES
 39:           USE MODCHARMM, ONLY : CHRMMT,NCHENCALLS,CHECKOMEGAT,CHECKCHIRALT,NOCISTRANS,MINOMEGA 38:           USE MODCHARMM, ONLY : CHRMMT,NCHENCALLS,CHECKOMEGAT,CHECKCHIRALT,NOCISTRANS,MINOMEGA
 40:           USE MODAMBER9, ONLY : NOCISTRANSRNA, NOCISTRANSDNA, GOODSTRUCTURE1, GOODSTRUCTURE2, CISARRAY1, CISARRAY2,CHIARRAY1,CHIARRAY2 39:           USE MODAMBER9, ONLY : NOCISTRANSRNA, NOCISTRANSDNA, GOODSTRUCTURE1, GOODSTRUCTURE2, CISARRAY1, CISARRAY2,CHIARRAY1,CHIARRAY2
 41:           USE MODMEC 40:           USE MODMEC
 42:           USE KEYUTILS 41:           USE KEYUTILS
 43:           USE COMMONS, ONLY : NINTS, PARAM1, PARAM2, PARAM3, ZSYM, EVDISTTHRESH, REDOPATHNEB, DEBUG 42:           USE COMMONS, ONLY : NINTS, PARAM1, PARAM2, PARAM3, ZSYM, EVDISTTHRESH, REDOPATHNEB, DEBUG
 44:           USE PORFUNCS 43:           USE PORFUNCS
 45:           USE AMHGLOBALS, ONLY : NMRES 44:           USE AMHGLOBALS, ONLY : NMRES
 46:           USE GSDATA, ONLY : GSITERDENSITY, GSCURITERD=>ITERD 45:           USE GSDATA, ONLY : GSITERDENSITY, GSCURITERD=>ITERD
 47: ! hk286 46: ! hk286
 48:           USE GENRIGID 47:           USE GENRIGID
 49:           USE CHIRALITY, ONLY: CIS_TRANS_CHECK, CHIRALITY_CHECK 48:           USE CHIRALITY, ONLY: CIS_TRANS_CHECK, CHIRALITY_CHECK
 50:           IMPLICIT NONE 49:           IMPLICIT NONE
 51:           DOUBLE PRECISION RMAT(3,3), DISTFAC 50:           DOUBLE PRECISION RMAT(3,3), DISTFAC
 52:  51: 
 53:           INTEGER,INTENT(IN) :: JS,JF 52:           INTEGER,INTENT(IN) :: JS,JF
 54:           LOGICAL,INTENT(IN) :: USEINT, USEINTLJ 53:           LOGICAL,INTENT(IN) :: USEINT, USEINTLJ
 55:           LOGICAL, ALLOCATABLE :: FOUNDBEFORE(:), DOAGAIN(:) 54:           LOGICAL, ALLOCATABLE :: FOUNDBEFORE(:), DOAGAIN(:)
 56:            55:           
 57:           INTEGER         :: I,I1,UNIQUE=0,MINPLUSPOS,MINMINUSPOS,J1,NC1,NC2,GLY_COUNT,LINTIMAGE,SAMEAS,ISTAT 56:           INTEGER         :: I,UNIQUE=0,MINPLUSPOS,MINMINUSPOS,J1,NC1,NC2,GLY_COUNT,LINTIMAGE,SAMEAS,ISTAT
 58:           DOUBLE PRECISION         :: EDUMMY,EDUMMY2,TMPTS(NOPT),SAVEPUSHOFF,SAVEMAXBFGS 57:           DOUBLE PRECISION         :: EDUMMY,EDUMMY2,TMPTS(NOPT),SAVEPUSHOFF,SAVEMAXBFGS
 59:           DOUBLE PRECISION,POINTER :: QPLUS(:),QMINUS(:),EPLUS,EMINUS 58:           DOUBLE PRECISION,POINTER :: QPLUS(:),QMINUS(:),EPLUS,EMINUS
 60:           LOGICAL         :: PLUSNEW,MINUSNEW,PATHFAILT,AMIDEFAIL,CHIRALFAIL,RERUN 59:           LOGICAL         :: PLUSNEW,MINUSNEW,PATHFAILT,AMIDEFAIL,CHIRALFAIL,RERUN
 61:           CHARACTER       :: ITSTRING*80, EOFSSTRING*80 60:           CHARACTER       :: ITSTRING*80, EOFSSTRING*80
 62:           DOUBLE PRECISION :: STARTINT(NINTS), FINISHINT(NINTS),DUM, LDUMMY 61:           DOUBLE PRECISION :: STARTINT(NINTS), FINISHINT(NINTS),DUM, LDUMMY
 63:           LOGICAL REDOPATH, REDOPATHXYZ, PERMUTE, EXISTS, NORERUN, PERMTEST 62:           LOGICAL REDOPATH, REDOPATHXYZ, PERMUTE, EXISTS, NORERUN, PERMTEST
 64:           DOUBLE PRECISION TSREDO(NOPT), ETSPREV, ETSDUM, DIST2, QP(NOPT), QM(NOPT), LGDUMMY(NOPT) 63:           DOUBLE PRECISION TSREDO(NOPT), ETSPREV, ETSDUM, DIST2, QP(NOPT), QM(NOPT), LGDUMMY(NOPT)
 65:           DOUBLE PRECISION NEBKFINALSAVE 64:           DOUBLE PRECISION NEBKFINALSAVE
 66:           INTEGER INVERT, INDEX(NATOMS), J2, IMATCH 65:           INTEGER INVERT, INDEX(NATOMS), J2, IMATCH
 67:           CHARACTER(LEN=2) ZDUM 66:           CHARACTER(LEN=2) ZDUM
 72:           INTEGER POSITION, NMINSAVE, NMINORIG 71:           INTEGER POSITION, NMINSAVE, NMINORIG
 73:           DOUBLE PRECISION DIST1P, DIST1M, DIST2P, DIST2M 72:           DOUBLE PRECISION DIST1P, DIST1M, DIST2P, DIST2M
 74:           LOGICAL PATHFAILED 73:           LOGICAL PATHFAILED
 75:           INTEGER :: NREDOPATHTRIES1=1 74:           INTEGER :: NREDOPATHTRIES1=1
 76:           INTEGER :: NREDOPATHTRIES2=1 75:           INTEGER :: NREDOPATHTRIES2=1
 77:           DOUBLE PRECISION :: REDOSTRETCH=5.0D0 76:           DOUBLE PRECISION :: REDOSTRETCH=5.0D0
 78:           TYPE (MINFOUNDTYPE) :: MYMINFOUND(NMINMAX) 77:           TYPE (MINFOUNDTYPE) :: MYMINFOUND(NMINMAX)
 79:           TYPE (TSFOUNDTYPE) :: MYTSFOUND(NTSMAX) 78:           TYPE (TSFOUNDTYPE) :: MYTSFOUND(NTSMAX)
 80: ! hk286 79: ! hk286
 81:           DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(NOPT) 80:           DOUBLE PRECISION :: XRIGIDCOORDS(DEGFREEDOMS), XCOORDS(NOPT)
 82:           DOUBLE PRECISION :: DISTFLAT, PATHLENGTH(3), EOFS(3), FLATSUM2, FLATSUM4, MINIM 
 83:  81: 
 84: ! 82: !
 85: !  We want to return here to rerun the path for transition states read in with REDOPATH 83: !  We want to return here to rerun the path for transition states read in with REDOPATH
 86: !  in case the connection fails. 84: !  in case the connection fails.
 87: ! 85: !
 88:           NC1=0 ! counter for changes of PUSHOFF value 86:           NC1=0 ! counter for changes of PUSHOFF value
 89:           NC2=0 ! counter for changes of BFGSSTEP value 87:           NC2=0 ! counter for changes of BFGSSTEP value
 90:           SAVEPUSHOFF=PUSHOFF 88:           SAVEPUSHOFF=PUSHOFF
 91:           SAVEMAXBFGS=MAXBFGS 89:           SAVEMAXBFGS=MAXBFGS
 92:           NMINORIG=NMIN 90:           NMINORIG=NMIN
608: !                    WRITE(*,'(A)') ' tryconnect> All of TS found appear to be new'606: !                    WRITE(*,'(A)') ' tryconnect> All of TS found appear to be new'
609: !               ENDIF607: !               ENDIF
610: !          ELSEIF (UNIQUE < NTSFOUND) THEN608: !          ELSEIF (UNIQUE < NTSFOUND) THEN
611: !               WRITE(CHR,'(i7)') unique 609: !               WRITE(CHR,'(i7)') unique 
612: !               WRITE(*,'(1X,A)') trim(adjustl(chr))//' of TS found appear to be new.'610: !               WRITE(*,'(1X,A)') trim(adjustl(chr))//' of TS found appear to be new.'
613: !          ELSEIF (UNIQUE ==0 .AND..NOT.NTSFOUND==0) THEN611: !          ELSEIF (UNIQUE ==0 .AND..NOT.NTSFOUND==0) THEN
614: !               WRITE(*,'(1X,A)') ' tryconnect> All of TS found are already known'612: !               WRITE(*,'(1X,A)') ' tryconnect> All of TS found are already known'
615: !          ENDIF613: !          ENDIF
616: 614: 
617:           ! path run for all unique ts615:           ! path run for all unique ts
618: ! 
619: !sy349: if the path got from dneb is regarded as flat, the highest-energy image 
620: !will be directly regarded as the transition state, and the two points will be 
621: !regarded connected. The judgement of whether the path is flat is in NEB/output.f90 
622: ! 
623:           IF (FLATPATHT) THEN 
624:              WRITE(*,'(A45,I3,A45,I3)') 'tryconnect > flat path found between minimum ', JS, ' and ', JF 
625:              CALL NEWCONNECTION(JS,JF,NTS) 
626:              CALL SETDISTANCE(JS,JF,0.0D0) 
627:              IF (INTERPCOSTFUNCTION) CALL SETINTERP(JS,JF,0.0D0) 
628:              CALL MKFNAMES(NTS,FILTH,FILTHSTR,ITSTRING,EOFSSTRING) 
629:  
630:              EOFS(1)=MI(JS)%DATA%E 
631:              EOFS(2)=TS(NTS)%DATA%E 
632:              EOFS(3)=MI(JF)%DATA%E 
633:  
634:              PATHLENGTH(1)=0.0D0 
635:              DISTFLAT=0.0D0 
636:              IF (BULKT) THEN 
637:                 DO I1=1,NATOMS 
638:                    DISTFLAT=DISTFLAT & 
639:    &               +(MI(JS)%DATA%X(3*(I1-1)+1)-TS(NTS)%DATA%X(3*(I1-1)+1) & 
640:    &            -PARAM1*NINT((MI(JS)%DATA%X(3*(I1-1)+1)-TS(NTS)%DATA%X(3*(I1-1)+1))/PARAM1))**2 & 
641:    &               +(MI(JS)%DATA%X(3*(I1-1)+2)-TS(NTS)%DATA%X(3*(I1-1)+2) & 
642:    &            -PARAM2*NINT((MI(JS)%DATA%X(3*(I1-1)+2)-TS(NTS)%DATA%X(3*(I1-1)+2))/PARAM2))**2 & 
643:    &               +(MI(JS)%DATA%X(3*(I1-1)+3)-TS(NTS)%DATA%X(3*(I1-1)+3) & 
644:    &            -PARAM3*NINT((MI(JS)%DATA%X(3*(I1-1)+3)-TS(NTS)%DATA%X(3*(I1-1)+3))/PARAM3))**2 
645:                 ENDDO 
646:                 PATHLENGTH(2)=SQRT(DISTFLAT) 
647:  
648:                 DISTFLAT=0.0D0 
649:                 DO I1=1,NATOMS 
650:                    DISTFLAT=DISTFLAT & 
651:    &               +(MI(JF)%DATA%X(3*(I1-1)+1)-TS(NTS)%DATA%X(3*(I1-1)+1) & 
652:    &            -PARAM1*NINT((MI(JF)%DATA%X(3*(I1-1)+1)-TS(NTS)%DATA%X(3*(I1-1)+1))/PARAM1))**2 & 
653:    &               +(MI(JF)%DATA%X(3*(I1-1)+2)-TS(NTS)%DATA%X(3*(I1-1)+2) & 
654:    &            -PARAM2*NINT((MI(JF)%DATA%X(3*(I1-1)+2)-TS(NTS)%DATA%X(3*(I1-1)+2))/PARAM2))**2 & 
655:    &               +(MI(JF)%DATA%X(3*(I1-1)+3)-TS(NTS)%DATA%X(3*(I1-1)+3) & 
656:    &            -PARAM3*NINT((MI(JF)%DATA%X(3*(I1-1)+3)-TS(NTS)%DATA%X(3*(I1-1)+3))/PARAM3))**2 
657:                 ENDDO 
658:                 PATHLENGTH(3)=SQRT(DISTFLAT)+PATHLENGTH(2) 
659:  
660:                 DISTFLAT=0.0D0 
661:                 DO I1=1,NATOMS 
662:                    DISTFLAT=DISTFLAT & 
663:    &               +(MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1) & 
664:    &            -PARAM1*NINT((MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1))/PARAM1))**2 & 
665:    &               +(MI(JS)%DATA%X(3*(I1-1)+2)-MI(JF)%DATA%X(3*(I1-1)+2) & 
666:    &            -PARAM2*NINT((MI(JS)%DATA%X(3*(I1-1)+2)-MI(JF)%DATA%X(3*(I1-1)+2))/PARAM2))**2 & 
667:    &               +(MI(JS)%DATA%X(3*(I1-1)+3)-MI(JF)%DATA%X(3*(I1-1)+3) & 
668:    &            -PARAM3*NINT((MI(JS)%DATA%X(3*(I1-1)+3)-MI(JF)%DATA%X(3*(I1-1)+3))/PARAM3))**2 
669:                 ENDDO 
670:                 FLATSUM2=DISTFLAT 
671:  
672:                 DISTFLAT=0.0D0 
673:                 DO I1=1,NATOMS 
674:                    DISTFLAT=DISTFLAT & 
675:    &               +((MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1) & 
676:    &            -PARAM1*NINT((MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1))/PARAM1))**2 & 
677:    &               +(MI(JS)%DATA%X(3*(I1-1)+2)-MI(JF)%DATA%X(3*(I1-1)+2) & 
678:    &            -PARAM2*NINT((MI(JS)%DATA%X(3*(I1-1)+2)-MI(JF)%DATA%X(3*(I1-1)+2))/PARAM2))**2 & 
679:    &               +(MI(JS)%DATA%X(3*(I1-1)+3)-MI(JF)%DATA%X(3*(I1-1)+3) & 
680:    &            -PARAM3*NINT((MI(JS)%DATA%X(3*(I1-1)+3)-MI(JF)%DATA%X(3*(I1-1)+3))/PARAM3))**2)**2 
681:                 ENDDO 
682:                 FLATSUM4=DISTFLAT 
683:              ELSE 
684:                 DO I1=1,NOPT 
685:                    DISTFLAT=DISTFLAT+(MI(JS)%DATA%X(I1)-TS(NTS)%DATA%X(I1))**2 
686:                 ENDDO 
687:                 PATHLENGTH(2)=SQRT(DISTFLAT) 
688:  
689:                 DISTFLAT=0.0D0 
690:                 DO I1=1,NOPT 
691:                    DISTFLAT=DISTFLAT+(MI(JF)%DATA%X(I1)-TS(NTS)%DATA%X(I1))**2 
692:                 ENDDO 
693:                 PATHLENGTH(3)=SQRT(DISTFLAT)+PATHLENGTH(2) 
694:                  
695:                 DISTFLAT=0.0D0 
696:                 DO I1=1,NOPT 
697:                    DISTFLAT=DISTFLAT+(MI(JS)%DATA%X(I1)-MI(JF)%DATA%X(I1))**2 
698:                 ENDDO 
699:                 FLATSUM2=DISTFLAT 
700:  
701:                 DISTFLAT=0.0D0 
702:                 DO I1=1,NATOMS 
703:                    DISTFLAT=DISTFLAT & 
704:    &                  +((MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1))**2 & 
705:    &                  +(MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1))**2 & 
706:    &                  +(MI(JS)%DATA%X(3*(I1-1)+1)-MI(JF)%DATA%X(3*(I1-1)+1))**2)**2 
707:                 ENDDO 
708:                 FLATSUM4=DISTFLAT 
709:              ENDIF 
710:               
711:              OPEN(UNIT=250,FILE=EOFSSTRING,STATUS='UNKNOWN') 
712:              DO I1=1,3 
713:                 WRITE(250,'(2G20.10,I6)') PATHLENGTH(I1), EOFS(I1), I1 
714:              ENDDO 
715:              CLOSE(250) 
716:  
717:              OPEN (UNIT=251,FILE=ITSTRING,STATUS='UNKNOWN') 
718:  
719:              WRITE(251,'(I6)') NATOMS 
720:              WRITE(251,'(A,G25.15)') 'Energy=', EOFS(1) 
721:              DO I1=1,NATOMS 
722:                 WRITE(251,'(A2,4X,3G20.10)') & 
723:    &               ZSYM(I1),MI(JS)%DATA%X(3*(I1-1)+1),MI(JS)%DATA%X(3*(I1-1)+2), & 
724:    &               MI(JS)%DATA%X(3*(I1-1)+3) 
725:              ENDDO 
726:  
727:              WRITE(251,'(I6)') NATOMS 
728:              WRITE(251,'(A,G25.15)') 'Energy=', EOFS(2) 
729:              DO I1=1,NATOMS 
730:                 WRITE(251,'(A2,4X,3G20.10)') & 
731:    &                 ZSYM(I1),TS(NTS)%DATA%X(3*(I1-1)+1),TS(NTS)%DATA%X(3*(I1-1)+2), & 
732:    &                 TS(NTS)%DATA%X(3*(I1-1)+3) 
733:              ENDDO 
734:  
735:              WRITE(251,'(I6)') NATOMS 
736:              WRITE(251,'(A,G25.15)') 'Energy=', EOFS(3) 
737:              DO I1=1,NATOMS 
738:                 WRITE(251,'(A2,4X,3G20.10)') & 
739:    &                 ZSYM(I1),MI(JF)%DATA%X(3*(I1-1)+1),MI(JF)%DATA%X(3*(I1-1)+2), & 
740:    &                 MI(JF)%DATA%X(3*(I1-1)+3) 
741:              ENDDO 
742:              CLOSE(251) 
743:               
744:              TS(NTS)%DATA%SLENGTH=PATHLENGTH(3) 
745:              TS(NTS)%DATA%DISP=SQRT(FLATSUM2) 
746:              TS(NTS)%DATA%GAMMA=FLATSUM4*NATOMS/FLATSUM2**2 
747:              TS(NTS)%DATA%NTILDE=FLATSUM2**2/FLATSUM4 
748:  
749:              IF (DUMPALLPATHS) CALL MAKEALLPATHINFO(TS(NTS)%DATA%X,MI(JS)%DATA%X,MI(JF)%DATA%X,EOFS(2),EOFS(1),EOFS(3),FRQSTS,FRQSPLUS,FRQSMINUS) 
750:  
751:          ELSE  
752:           DO I=NTS-UNIQUE+1,NTS616:           DO I=NTS-UNIQUE+1,NTS
753:                WRITE(CHR,'(i5)') I617:                WRITE(CHR,'(i5)') I
754:                PRINT '(/1x,a)', '>>>>>  Path run for ts '//trim(adjustl(chr))//' ...'618:                PRINT '(/1x,a)', '>>>>>  Path run for ts '//trim(adjustl(chr))//' ...'
755:                IF (DOAGAIN(I-NTS+UNIQUE).AND.(.NOT.REDOPATH)) THEN619:                IF (DOAGAIN(I-NTS+UNIQUE).AND.(.NOT.REDOPATH)) THEN
756:                   PUSHOFF=PUSHOFF/10.0D0620:                   PUSHOFF=PUSHOFF/10.0D0
757:                   PRINT '(A,G20.10)',' tryconnect> Trying this transition state path again with pushoff=',PUSHOFF621:                   PRINT '(A,G20.10)',' tryconnect> Trying this transition state path again with pushoff=',PUSHOFF
758:                ENDIF622:                ENDIF
759:                ALLOCATE( QPLUS(NOPT),QMINUS(NOPT),EPLUS,EMINUS )623:                ALLOCATE( QPLUS(NOPT),QMINUS(NOPT),EPLUS,EMINUS )
760:                CALL MKFNAMES(I,FILTH,FILTHSTR,ITSTRING,EOFSSTRING)624:                CALL MKFNAMES(I,FILTH,FILTHSTR,ITSTRING,EOFSSTRING)
761:                EDUMMY=TS(I)%DATA%E625:                EDUMMY=TS(I)%DATA%E
1099:                          CALL ADDNEWMIN(EMINUS,QMINUS)963:                          CALL ADDNEWMIN(EMINUS,QMINUS)
1100:                     ENDIF964:                     ENDIF
1101:                     CALL NEWCONNECTION(MINPLUSPOS,MINMINUSPOS,I)965:                     CALL NEWCONNECTION(MINPLUSPOS,MINMINUSPOS,I)
1102:                     CALL SETDISTANCE(MINPLUSPOS,MINMINUSPOS,0.0D0)966:                     CALL SETDISTANCE(MINPLUSPOS,MINMINUSPOS,0.0D0)
1103:                     IF (INTERPCOSTFUNCTION) CALL SETINTERP(MINPLUSPOS,MINMINUSPOS,0.0D0)967:                     IF (INTERPCOSTFUNCTION) CALL SETINTERP(MINPLUSPOS,MINMINUSPOS,0.0D0)
1104:                ENDIF968:                ENDIF
1105:                IF (DOAGAIN(I-NTS+UNIQUE).AND.(.NOT.REDOPATH)) THEN969:                IF (DOAGAIN(I-NTS+UNIQUE).AND.(.NOT.REDOPATH)) THEN
1106:                   PUSHOFF=SAVEPUSHOFF970:                   PUSHOFF=SAVEPUSHOFF
1107:                ENDIF971:                ENDIF
1108:           ENDDO972:           ENDDO
1109:         ENDIF   
1110:            
1111:           !IF(ALLOCATED(FLATPATHT)) DEALLOCATE(FLATPATHT) 
1112: !973: !
1113: ! Allow for new pathway calculation with different PUSHOFF and MAXBFGS974: ! Allow for new pathway calculation with different PUSHOFF and MAXBFGS
1114: !975: !
1115:           IF (REDOPATH) THEN976:          IF (REDOPATH) THEN
1116:              CALL MINPERMDIST(QP,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)977:             CALL MINPERMDIST(QP,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1117:              DIST1P=D978:             DIST1P=D
1118: 979: 
1119:              CALL MINPERMDIST(QM,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)980:             CALL MINPERMDIST(QM,MIN1REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1120:              DIST1M=D981:             DIST1M=D
1121: 982: 
1122:              CALL MINPERMDIST(QM,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)983:             CALL MINPERMDIST(QM,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1123:              DIST2M=D984:             DIST2M=D
1124: 985: 
1125:              CALL MINPERMDIST(QP,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)986:             CALL MINPERMDIST(QP,MIN2REDO,NATOMS,DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1126:              DIST2P=D987:             DIST2P=D
1127: 988: 
1128:              PATHFAILED=.FALSE.989:             PATHFAILED=.FALSE.
1129:              IF ((DIST1P.GT.GEOMDIFFTOL).AND.(DIST1M.GT.GEOMDIFFTOL)) THEN990:             IF ((DIST1P.GT.GEOMDIFFTOL).AND.(DIST1M.GT.GEOMDIFFTOL)) THEN
1130:                 PRINT '(A)',' tryconnect> path failed to match first minimum'991:                PRINT '(A)',' tryconnect> path failed to match first minimum'
1131:                 PATHFAILED=.TRUE.992:                PATHFAILED=.TRUE.
1132:              ENDIF993:             ENDIF
1133:              IF ((DIST2P.GT.GEOMDIFFTOL).AND.(DIST2M.GT.GEOMDIFFTOL)) THEN994:             IF ((DIST2P.GT.GEOMDIFFTOL).AND.(DIST2M.GT.GEOMDIFFTOL)) THEN
1134:                 PRINT '(A)',' tryconnect> path failed to match second minimum'995:                PRINT '(A)',' tryconnect> path failed to match second minimum'
1135:                 PATHFAILED=.TRUE.996:                PATHFAILED=.TRUE.
1136:              ENDIF997:             ENDIF
1137:              IF (PATHFAILED) THEN998:             IF (PATHFAILED) THEN
1138:                 NC1=NC1+1999:                NC1=NC1+1
1139:                 IF (NC1.GT.2*NREDOPATHTRIES1) THEN1000:                IF (NC1.GT.2*NREDOPATHTRIES1) THEN
1140:                    NC2=NC2+11001:                   NC2=NC2+1
1141:                    NC1=01002:                   NC1=0
1142:                 ENDIF1003:                ENDIF
1143:                 IF (NREDOPATHTRIES1.EQ.0) THEN1004:                IF (NREDOPATHTRIES1.EQ.0) THEN
1144:                    PUSHOFF=SAVEPUSHOFF1005:                   PUSHOFF=SAVEPUSHOFF
1145:                 ELSEIF (NC1.GT.NREDOPATHTRIES1) THEN1006:                ELSEIF (NC1.GT.NREDOPATHTRIES1) THEN
1146:                    PUSHOFF=SAVEPUSHOFF+(NC1-NREDOPATHTRIES1)*(REDOSTRETCH*SAVEPUSHOFF-SAVEPUSHOFF)/NREDOPATHTRIES11007:                   PUSHOFF=SAVEPUSHOFF+(NC1-NREDOPATHTRIES1)*(REDOSTRETCH*SAVEPUSHOFF-SAVEPUSHOFF)/NREDOPATHTRIES1
1147:                 ELSE1008:                ELSE
1148:                    PUSHOFF=SAVEPUSHOFF-NC1*(SAVEPUSHOFF-SAVEPUSHOFF/REDOSTRETCH)/NREDOPATHTRIES11009:                   PUSHOFF=SAVEPUSHOFF-NC1*(SAVEPUSHOFF-SAVEPUSHOFF/REDOSTRETCH)/NREDOPATHTRIES1
1149:                 ENDIF1010:                ENDIF
1150:                 IF (NREDOPATHTRIES2.EQ.0) THEN1011:                IF (NREDOPATHTRIES2.EQ.0) THEN
1151:                    MAXBFGS=SAVEMAXBFGS1012:                   MAXBFGS=SAVEMAXBFGS
1152:                 ELSEIF (NC2.GT.NREDOPATHTRIES2) THEN1013:                ELSEIF (NC2.GT.NREDOPATHTRIES2) THEN
1153:                    MAXBFGS=SAVEMAXBFGS+(NC2-NREDOPATHTRIES2)*(REDOSTRETCH*SAVEMAXBFGS-SAVEMAXBFGS)/NREDOPATHTRIES21014:                   MAXBFGS=SAVEMAXBFGS+(NC2-NREDOPATHTRIES2)*(REDOSTRETCH*SAVEMAXBFGS-SAVEMAXBFGS)/NREDOPATHTRIES2
1154:                 ELSE1015:                ELSE
1155:                    MAXBFGS=SAVEMAXBFGS-NC2*(SAVEMAXBFGS-SAVEMAXBFGS/REDOSTRETCH)/NREDOPATHTRIES21016:                   MAXBFGS=SAVEMAXBFGS-NC2*(SAVEMAXBFGS-SAVEMAXBFGS/REDOSTRETCH)/NREDOPATHTRIES2
1156:                 ENDIF1017:                ENDIF
1157:                 IF (NC2.LE.2*NREDOPATHTRIES2) THEN1018:                IF (NC2.LE.2*NREDOPATHTRIES2) THEN
1158:                    PRINT '(2(A,F15.5))',' tryconnect> Redo path with pushoff=',PUSHOFF,' maxbfgs=',MAXBFGS1019:                   PRINT '(2(A,F15.5))',' tryconnect> Redo path with pushoff=',PUSHOFF,' maxbfgs=',MAXBFGS
1159:                    RERUN=.TRUE.1020:                   RERUN=.TRUE.
1160:                    GOTO 101021:                   GOTO 10
1161:                 ENDIF1022:                ENDIF
1162:              ENDIF1023:             ENDIF
1163:            ENDIF1024:           ENDIF
 1025: 
 1026:           RERUN=.FALSE.
 1027:           PUSHOFF=SAVEPUSHOFF
 1028:           MAXBFGS=SAVEMAXBFGS
 1029:           DEALLOCATE(FOUNDBEFORE,DOAGAIN)
1164: 1030: 
1165:            RERUN=.FALSE. 
1166:            PUSHOFF=SAVEPUSHOFF 
1167:            MAXBFGS=SAVEMAXBFGS 
1168:            DEALLOCATE(FOUNDBEFORE,DOAGAIN) 
1169:   
1170:      END SUBROUTINE TRYCONNECT1031:      END SUBROUTINE TRYCONNECT
1171: 1032: 
1172: END MODULE TRYCONNECTMODULE1033: END MODULE TRYCONNECTMODULE


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0