hdiff output

r33340/ellipticintegral.f90 2017-09-23 12:30:15.648993832 +0100 r33339/ellipticintegral.f90 2017-09-23 12:30:16.989011646 +0100
 11:       TINY=1.5e-38 11:       TINY=1.5e-38
 12:       BIG=3.E37 12:       BIG=3.E37
 13:       THIRD=1.D0/3. 13:       THIRD=1.D0/3.
 14:       C1=1.D0/24. 14:       C1=1.D0/24.
 15:       C2=.1D0 15:       C2=.1D0
 16:       C3=3.D0/44. 16:       C3=3.D0/44.
 17:       C4=1.D0/14. 17:       C4=1.D0/14.
 18:  18: 
 19:       IF(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or. max(x,y,z).gt.BIG) THEN 19:       IF(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or. max(x,y,z).gt.BIG) THEN
 20:          PRINT *, x, y, z 20:          PRINT *, x, y, z
 21:          STOP ! pause 'invalid arguments in rf' 21:          pause 'invalid arguments in rf'
 22:       ENDIF 22:       ENDIF
 23:  23: 
 24:       xt=x 24:       xt=x
 25:       yt=y 25:       yt=y
 26:       zt=z 26:       zt=z
 27:  1    continue 27:  1    continue
 28:  28: 
 29:       sqrtx=sqrt(xt) 29:       sqrtx=sqrt(xt)
 30:       sqrty=sqrt(yt) 30:       sqrty=sqrt(yt)
 31:       sqrtz=sqrt(zt) 31:       sqrtz=sqrt(zt)
 60:       ERRTOL=.0015D0 60:       ERRTOL=.0015D0
 61:       TINY=1.e-25 61:       TINY=1.e-25
 62:       BIG=4.5E21 62:       BIG=4.5E21
 63:       C1=3.D0/14. 63:       C1=3.D0/14.
 64:       C2=1.D0/6. 64:       C2=1.D0/6.
 65:       C3=9.D0/22. 65:       C3=9.D0/22.
 66:       C4=3.D0/26. 66:       C4=3.D0/26.
 67:       C5=.25D0*C3 67:       C5=.25D0*C3
 68:       C6=1.5D0*C4 68:       C6=1.5D0*C4
 69:  69: 
 70:       if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y,z).gt.BIG) STOP ! pause 'invalid arguments in rd' 70:       if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y,z).gt.BIG)pause 'invalid arguments in rd'
 71:       xt=x 71:       xt=x
 72:       yt=y 72:       yt=y
 73:       zt=z 73:       zt=z
 74:       sum=0. 74:       sum=0.
 75:       fac=1. 75:       fac=1.
 76:  1    continue 76:  1    continue
 77:  77: 
 78:       sqrtx=sqrt(xt) 78:       sqrtx=sqrt(xt)
 79:       sqrty=sqrt(yt) 79:       sqrty=sqrt(yt)
 80:       sqrtz=sqrt(zt) 80:       sqrtz=sqrt(zt)


r33340/intlbfgs.f90 2017-09-23 12:30:15.876996861 +0100 r33339/intlbfgs.f90 2017-09-23 12:30:17.505018508 +0100
 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, & 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, &
 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, & 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, &
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, & 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, &
 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, & 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, &
 32:      & QCIADDREP, QCIXYZ, WHOLEDNEB, QCIIMAGE, FROZEN, QCIRESTART, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, & 32:      & QCIADDREP, QCIXYZ, WHOLEDNEB, QCIIMAGE, FROZEN
 33:      & PERMDIST, LOCALPERMCUT, QCILPERMDIST, QCIPDINT, QCIPERMCUT 
 34: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3 33: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3
 35: USE MODCHARMM, ONLY : CHRMMT 34: USE MODCHARMM, ONLY : CHRMMT
 36:  35: 
 37: IMPLICIT NONE  36: IMPLICIT NONE 
 38:  37: 
 39: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points 38: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points
 40: INTEGER D, U 39: INTEGER D, U
 41: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE 40: DOUBLE PRECISION DIST, DIST2, RMAT(3,3)
 42: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ 41: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ
 43: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF   42: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2
 44: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS) 43: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG
 45: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS) 
 46: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 44: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 47:  45: 
 48: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY 46: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY
 49: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF 
 50: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2 47: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2
 51: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE 48: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE
 52: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS) 49: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX
 53: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS) 50: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS)
 54: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, & 51: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, &
 55:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), & 52:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), &
 56:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS) 53:   &                 BESTWORST, WORST
 57: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA 54: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA
 58: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS) 55: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS)
 59: LOGICAL SWITCHED 56: LOGICAL SWITCHED
 60: DOUBLE PRECISION, POINTER :: X(:), G(:) 57: DOUBLE PRECISION, POINTER :: X(:), G(:)
 61: ! 58: !
 62: ! efk: for freezenodes 59: ! efk: for freezenodes
 63: ! 60: !
 64: DOUBLE PRECISION :: TESTG, TOTGNORM 61: DOUBLE PRECISION :: TESTG, TOTGNORM
 65: INTEGER :: IM 62: INTEGER :: IM
 66: ! 63: !
 67: ! Dimensions involving INTIMAGE 64: ! Dimensions involving INTIMAGE
 68: ! 65: !
 69: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), & 66: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), &
 70:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), & 67:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), &
 71:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:) 68:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:)
 72: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:) 69: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:)
 73: ! saved interpolation 70: ! saved interpolation
 74: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE 71: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE
 75: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:) 72: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:)
 76: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP) 73: LOGICAL READIMAGET
 77: INTEGER LUNIT, GETUNIT 74: INTEGER LUNIT, GETUNIT
 78: CHARACTER(LEN=2) SDUMMY 75: CHARACTER(LEN=2) SDUMMY
 79: INTEGER JMAXEEE,JMAXRMS 76: INTEGER JMAXEEE,JMAXRMS
 80: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT 77: DOUBLE PRECISION MAXEEE,MAXRMS
 81:  78: 
 82: WHOLEDNEB=.FALSE. 79: WHOLEDNEB=.FALSE.
  80: READIMAGET=.TRUE.
 83: READIMAGET=.FALSE. 81: READIMAGET=.FALSE.
 84: REMOVEIMAGE=.FALSE. 
 85: IF (QCIRESTART) READIMAGET=.TRUE. 
 86: IF (READIMAGET) THEN 82: IF (READIMAGET) THEN
 87:    LUNIT=GETUNIT() 83:    LUNIT=GETUNIT()
 88:    OPEN(UNIT=LUNIT,FILE='int.xyz',STATUS='OLD') 84:    OPEN(UNIT=LUNIT,FILE='restart.xyz',STATUS='OLD')
 89:    INTIMAGE=0 85:    INTIMAGE=0
 90: 653 CONTINUE 86: 653 CONTINUE
 91:    READ(LUNIT,*,END=654) NDUMMY 87:    READ(LUNIT,*,END=654) NDUMMY
 92:    READ(LUNIT,*)  88:    READ(LUNIT,*) 
 93:    DO J1=1,NATOMS 89:    DO J1=1,NATOMS
 94:       READ(LUNIT,*) SDUMMY, DUMMYX, DUMMYY, DUMMYZ 90:       READ(LUNIT,*) SDUMMY, DUMMYX, DUMMYY, DUMMYZ
 95: !     WRITE(*,'(A,I6,A2,3G20.10)') 'J1,sd,xd,yd,zd=',J1,SDUMMY, DUMMYX, DUMMYY, DUMMYZ 91: !     WRITE(*,'(A,I6,A2,3G20.10)') 'J1,sd,xd,yd,zd=',J1,SDUMMY, DUMMYX, DUMMYY, DUMMYZ
 96:    ENDDO 92:    ENDDO
 97:    INTIMAGE=INTIMAGE+1 93:    INTIMAGE=INTIMAGE+1
 98:    GOTO 653 94:    GOTO 653
 99: 654 CONTINUE 95: 654 CONTINUE
100:    INTIMAGE=INTIMAGE-2 96:    INTIMAGE=INTIMAGE-2
101:    WRITE(*,'(A,I10,A)') 'intlbfgs> Rereading ',INTIMAGE,' frames' 97:    WRITE(*,'(A,I10,A)') 'intlbfgs> Rereading ',INTIMAGE,' frames'
102:    CLOSE(LUNIT) 
103: ENDIF 98: ENDIF
104:  99: 
105: ALLOCATE(TRUEEE(INTIMAGE+2), &100: ALLOCATE(TRUEEE(INTIMAGE+2), &
106:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &101:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &
107:   &      GTMP(3*NATOMS*INTIMAGE), &102:   &      GTMP(3*NATOMS*INTIMAGE), &
108:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), &103:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), &
109:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &104:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &
110:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &105:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &
111:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))106:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))
112: 107: 
144: IF (.NOT.ALLOCATED(CONI)) THEN 139: IF (.NOT.ALLOCATED(CONI)) THEN 
145:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX))140:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX))
146:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))141:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
147: ENDIF142: ENDIF
148: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))143: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
149: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))144: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
150: !145: !
151: ! Initialise XYZ146: ! Initialise XYZ
152: !147: !
153: IF (READIMAGET) THEN148: IF (READIMAGET) THEN
154:    LUNIT=GETUNIT()149:    REWIND(LUNIT)
155:    OPEN(UNIT=LUNIT,FILE='int.xyz',STATUS='OLD') 
156:    DO J2=1,INTIMAGE+2150:    DO J2=1,INTIMAGE+2
157:       READ(LUNIT,*) NDUMMY151:       READ(LUNIT,*) NDUMMY
158:       READ(LUNIT,*) 152:       READ(LUNIT,*) 
159:       DO J1=1,NATOMS153:       DO J1=1,NATOMS
160:          READ(LUNIT,*) SDUMMY,XYZ(3*NATOMS*(J2-1)+3*(J1-1)+1),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+2),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+3)154:          READ(LUNIT,*) SDUMMY,XYZ(3*NATOMS*(J2-1)+3*(J1-1)+1),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+2),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+3)
161:       ENDDO155:       ENDDO
162:    ENDDO156:    ENDDO
163:    CLOSE(LUNIT)157:    CLOSE(LUNIT)
164: ELSE158: ELSE
165:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))159:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))
243:    WRITE(*,'(A)') ' intlbfgs> All atoms move less than threshold - skip to linear interpolation for end points'237:    WRITE(*,'(A)') ' intlbfgs> All atoms move less than threshold - skip to linear interpolation for end points'
244:    INTIMAGE=0238:    INTIMAGE=0
245:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))239:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))
246:    XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS))240:    XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS))
247:    DO J1=1,INTIMAGE+2241:    DO J1=1,INTIMAGE+2
248:       XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)242:       XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)
249:    ENDDO243:    ENDDO
250:    GOTO 678244:    GOTO 678
251: ENDIF245: ENDIF
252: 246: 
253: NACTIVE=0 
254: TURNONORDER(1:NATOMS)=0 
255: ATOMACTIVE(1:NATOMS)=.FALSE. 
256: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT**6 ! also needed for congrad.f90 potential 
257: IF (READIMAGET) THEN247: IF (READIMAGET) THEN
258:    LUNIT=GETUNIT()248:    NACTIVE=NATOMS
259:    OPEN(UNIT=LUNIT,FILE='QCIdump',STATUS='OLD')249:    DO J1=1,NATOMS
260: 250:       TURNONORDER(J1)=J1 ! fake initialisation
261:    READ(LUNIT,*) NACTIVE251:    ENDDO
262:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NACTIVE,' active atoms'252:    ATOMACTIVE(1:NATOMS)=.TRUE.
263:    WRITE(*,'(A)') ' intlbfgs> reading turnonorder'253:    CONACTIVE(1:NCONSTRAINT)=.TRUE.
264:    READ(LUNIT,*) TURNONORDER(1:NACTIVE) 
265:    WRITE(*,'(12I8)') TURNONORDER(1:NACTIVE) 
266:    WRITE(*,'(A)') ' intlbfgs> reading atomactive' 
267:    READ(LUNIT,*) ATOMACTIVE(1:NATOMS) 
268:    WRITE(*,'(12L5)') ATOMACTIVE(1:NATOMS)  
269:    READ(LUNIT,*) NCONSTRAINT 
270:    WRITE(*,'(A)') ' intlbfgs> reading conactive' 
271:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NCONSTRAINT,' constraints' 
272:    ALLOCATE(CONACTIVE(NCONSTRAINT)) 
273:    READ(LUNIT,*) CONACTIVE(1:NCONSTRAINT) 
274:    WRITE(*,'(12L5)') CONACTIVE(1:NCONSTRAINT)  
275:  
276:    ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT)) 
277:    ALLOCATE(CONCUTLOCAL(NCONSTRAINT)) 
278:    CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT) 
279:    CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT) 
280:    WRITE(*,'(A)') ' intlbfgs> CONCUT:' 
281:    WRITE(*,'(6G20.10)') CONCUT(1:NCONSTRAINT) 
282:    WRITE(*,'(A)') ' intlbfgs> CONDISTREF:' 
283:    WRITE(*,'(6G20.10)') CONDISTREF(1:NCONSTRAINT) 
284:    WRITE(*,'(A)') ' intlbfgs> CONI:' 
285:    WRITE(*,'(12I8)') CONI(1:NCONSTRAINT) 
286:    WRITE(*,'(A)') ' intlbfgs> CONJ:' 
287:    WRITE(*,'(12I8)') CONJ(1:NCONSTRAINT) 
288:  
289:    READ(LUNIT,*) NREPULSIVE,NNREPULSIVE,NREPMAX 
290:    USEFRAC=1.0D0 
291:    WRITE(*,'(A,3I10,G20.10)') 'intlbfgs> NREPULSIVE,NNREPULSIVE,NREPMAX=',NREPULSIVE,NNREPULSIVE,NREPMAX 
292:    IF (ALLOCATED(REPI)) DEALLOCATE(REPI) 
293:    IF (ALLOCATED(REPJ)) DEALLOCATE(REPJ) 
294:    IF (ALLOCATED(NREPI)) DEALLOCATE(NREPI) 
295:    IF (ALLOCATED(NREPJ)) DEALLOCATE(NREPJ) 
296:    IF (ALLOCATED(REPCUT)) DEALLOCATE(REPCUT) 
297:    IF (ALLOCATED(NREPCUT)) DEALLOCATE(NREPCUT) 
298:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX)) 
299:    READ(LUNIT,*) REPI(1:NREPULSIVE) 
300:    WRITE(*,'(A)') ' intlbfgs> read REPI:' 
301:    WRITE(*,'(12I8)') REPI(1:NREPULSIVE) 
302:    READ(LUNIT,*) REPJ(1:NREPULSIVE) 
303:    WRITE(*,'(A)') ' intlbfgs> read REPJ:' 
304:    WRITE(*,'(12I8)') REPJ(1:NREPULSIVE) 
305:    READ(LUNIT,*) NREPI(1:NNREPULSIVE) 
306:    WRITE(*,'(A)') ' intlbfgs> read NREPI:' 
307:    WRITE(*,'(12I8)') NREPI(1:NNREPULSIVE) 
308:    READ(LUNIT,*) NREPJ(1:NNREPULSIVE) 
309:    WRITE(*,'(A)') ' intlbfgs> read NREPJ:' 
310:    WRITE(*,'(12I8)') NREPJ(1:NNREPULSIVE) 
311:    READ(LUNIT,*) REPCUT(1:NREPULSIVE) 
312:    WRITE(*,'(A)') ' intlbfgs> read REPCUT:' 
313:    WRITE(*,'(6G20.10)') REPCUT(1:NREPULSIVE) 
314:    READ(LUNIT,*) NREPCUT(1:NNREPULSIVE) 
315:    WRITE(*,'(A)') ' intlbfgs> read NREPCUT:' 
316:    WRITE(*,'(6G20.10)') NREPCUT(1:NNREPULSIVE) 
317:    CLOSE(LUNIT) 
318:  
319:    GLAST(1:D)=G(1:D)254:    GLAST(1:D)=G(1:D)
320:    XSAVE(1:D)=X(1:D)255:    XSAVE(1:D)=X(1:D)
321:    GOTO 986256:    GOTO 986
322: ENDIF257: ENDIF
 258: NACTIVE=0
 259: TURNONORDER(1:NATOMS)=0
 260: ATOMACTIVE(1:NATOMS)=.FALSE.
323: IF (INTFREEZET) THEN261: IF (INTFREEZET) THEN
324:    DO J1=1,NATOMS262:    DO J1=1,NATOMS
325:       IF (INTFROZEN(J1)) THEN263:       IF (INTFROZEN(J1)) THEN
326: ! 264: ! 
327: ! linear interpolation 265: ! linear interpolation 
328: ! 266: ! 
329:          DO J2=2,INTIMAGE+1267:          DO J2=2,INTIMAGE+1
330:             XYZ((J2-1)*3*NATOMS+3*(J1-1)+1:(J2-1)*3*NATOMS+3*(J1-1)+3)= &268:             XYZ((J2-1)*3*NATOMS+3*(J1-1)+1:(J2-1)*3*NATOMS+3*(J1-1)+3)= &
331:   &            (INTIMAGE-J2+2)*XYZ(3*(J1-1)+1:3*(J1-1)+3)/(INTIMAGE+1) &269:   &            (INTIMAGE-J2+2)*XYZ(3*(J1-1)+1:3*(J1-1)+3)/(INTIMAGE+1) &
332:   &           +(J2-1)*XYZ(3*NATOMS*(INTIMAGE+1)+3*(J1-1)+1:3*NATOMS*(INTIMAGE+1)+3*(J1-1)+3)/(INTIMAGE+1)270:   &           +(J2-1)*XYZ(3*NATOMS*(INTIMAGE+1)+3*(J1-1)+1:3*NATOMS*(INTIMAGE+1)+3*(J1-1)+3)/(INTIMAGE+1)
364:    XSAVE(1:D)=X(1:D)302:    XSAVE(1:D)=X(1:D)
365:    GOTO 567303:    GOTO 567
366: ENDIF304: ENDIF
367: DO J1=1,NCONSTRAINT305: DO J1=1,NCONSTRAINT
368:    DF=SQRT((XYZ(3*(CONI(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+1))**2 &306:    DF=SQRT((XYZ(3*(CONI(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+1))**2 &
369:   &       +(XYZ(3*(CONI(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+2))**2 &307:   &       +(XYZ(3*(CONI(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+2))**2 &
370:   &       +(XYZ(3*(CONI(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+3))**2)&308:   &       +(XYZ(3*(CONI(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+3))**2)&
371:   &  +SQRT((XYZ(3*(CONJ(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+1))**2 &309:   &  +SQRT((XYZ(3*(CONJ(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+1))**2 &
372:   &       +(XYZ(3*(CONJ(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+2))**2 &310:   &       +(XYZ(3*(CONJ(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+2))**2 &
373:   &       +(XYZ(3*(CONJ(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+3))**2)311:   &       +(XYZ(3*(CONJ(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+3))**2)
 312: !  IF (J1.EQ.3505) THEN
 313: !     WRITE(*,'(A,3I10)') 'intlbfgs> J1,CONI(J1),CONJ(J1)=',J1,CONI(J1),CONJ(J1)
 314: !     WRITE(*,'(6G20.10)') XYZ(3*(CONI(J1)-1)+1:3*(CONI(J1)-1)+3), &
 315: ! &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+1:(INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+3)
 316: !     WRITE(*,'(6G20.10)') XYZ(3*(CONJ(J1)-1)+1:3*(CONJ(J1)-1)+3), &
 317: ! &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+1:(INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+3)
 318: !  ENDIF
374:    IF (DF.LT.DUMMY) THEN319:    IF (DF.LT.DUMMY) THEN
375:       NBEST=J1320:       NBEST=J1
376:       DUMMY=DF321:       DUMMY=DF
377:    ENDIF322:    ENDIF
378:    IF (DF.GT.DUMMY2) THEN323:    IF (DF.GT.DUMMY2) THEN
379:       NBEST2=J1324:       NBEST2=J1
380:       DUMMY2=DF325:       DUMMY2=DF
381:    ENDIF326:    ENDIF
382: ENDDO327: ENDDO
383: IF (DEBUG) WRITE(*,'(A,I6,A,2I6,A,F15.5)') ' intlbfgs> Smallest overall motion for constraint ',NBEST, ' atoms ', &328: IF (DEBUG) WRITE(*,'(A,I6,A,2I6,A,F15.5)') ' intlbfgs> Smallest overall motion for constraint ',NBEST, ' atoms ', &
491:       DMIN=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)436:       DMIN=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)
492:       NREPULSIVE=NREPULSIVE+1437:       NREPULSIVE=NREPULSIVE+1
493:       IF (NREPULSIVE.GT.NREPMAX) CALL REPDOUBLE438:       IF (NREPULSIVE.GT.NREPMAX) CALL REPDOUBLE
494:       REPI(NREPULSIVE)=J1439:       REPI(NREPULSIVE)=J1
495:       REPJ(NREPULSIVE)=CONJ(NBEST)440:       REPJ(NREPULSIVE)=CONJ(NBEST)
496:       REPCUT(NREPULSIVE)=DMIN441:       REPCUT(NREPULSIVE)=DMIN
497:       IF (DEBUG) WRITE(*,'(A,I6,A,I6,A,F15.5)') ' intlbfgs> Adding repulsion for new atom ',CONJ(NBEST),' with atom ',J1, &442:       IF (DEBUG) WRITE(*,'(A,I6,A,I6,A,F15.5)') ' intlbfgs> Adding repulsion for new atom ',CONJ(NBEST),' with atom ',J1, &
498:   &                                          ' cutoff=',DMIN443:   &                                          ' cutoff=',DMIN
499: 541   CONTINUE444: 541   CONTINUE
500:    ENDDO445:    ENDDO
501: ENDIF ! end of block to add constraints and repulsions for frozen atoms.446: ENDIF ! end of block to add constraints and repulstions for frozen atoms.
502: CALL MYCPU_TIME(FTIME,.FALSE.)447: CALL MYCPU_TIME(FTIME,.FALSE.)
503: WRITE(*,'(A,F10.1,A,I6)') ' intlbfgs> constrained potential finished, time=',FTIME-STIME,' number of repulsions=',NREPULSIVE448: WRITE(*,'(A,F10.1,A,I6)') ' intlbfgs> constrained potential finished, time=',FTIME-STIME,' number of repulsions=',NREPULSIVE
504: 986 CONTINUE449: 986 CONTINUE
505: STIME=FTIME450: STIME=FTIME
506: NSTEPSMAX=INTSTEPS1451: NSTEPSMAX=INTSTEPS1
507: !452: !
508: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.453: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.
509: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!454: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!
510: !455: !
511: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.456: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.
529:   &                       ETOTAL/INTIMAGE,COLDFUSIONLIMIT474:   &                       ETOTAL/INTIMAGE,COLDFUSIONLIMIT
530:    DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)475:    DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)
531:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &476:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
532:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)477:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
533:    INTIMAGE=INTIMAGESAVE478:    INTIMAGE=INTIMAGESAVE
534:    RETURN479:    RETURN
535: ENDIF480: ENDIF
536: 481: 
537: ! IF (DEBUG) WRITE(*,'(A6,A20,A20,A9,A9)') 'Iter','Energy per image','RMS Force','Step'482: ! IF (DEBUG) WRITE(*,'(A6,A20,A20,A9,A9)') 'Iter','Energy per image','RMS Force','Step'
538: 483: 
539: ! 
540: ! In this block PERMGROUP(NDUMMY+J2-1) counts through the atom indices specified in perm.allow, 
541: ! one group at a time. 
542: ! 
543: IF (PERMDIST) THEN 
544:    PERMUTABLE(1:NATOMS)=.FALSE. 
545:    INGROUP(1:NATOMS)=0 
546:    NDUMMY=1 
547:    DO J1=1,NPERMGROUP 
548:       GROUPACTIVE(J1)=.FALSE. 
549:       DO J2=1,NPERMSIZE(J1) 
550:          PERMUTABLE(PERMGROUP(NDUMMY+J2-1))=.TRUE. 
551:          INGROUP(PERMGROUP(NDUMMY+J2-1))=J1 
552:       ENDDO 
553:       NDUMMY=NDUMMY+NPERMSIZE(J1) 
554:    ENDDO 
555: ENDIF 
556:  
557: 567 CONTINUE484: 567 CONTINUE
558: 485: 
559: DO ! Main do loop with counter NITERDONE, initially set to one486: DO ! Main do loop with counter NITERDONE, initially set to one
560:  
561: ! 
562: !  Check permutational alignments. Maintain a list of the permutable groups where all 
563: !  members are active. See if we have any new complete groups. MUST update NDUMMY 
564: !  counter to step through permutable atom list. 
565: ! 
566: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN 
567:    NDUMMY=1 
568:    DO J1=1,NPERMGROUP 
569:       IF (GROUPACTIVE(J1)) GOTO 975 
570:       DO J2=1,NPERMSIZE(J1) 
571:          IF (.NOT.ATOMACTIVE(PERMGROUP(NDUMMY+J2-1))) GOTO 975 
572:       ENDDO 
573:       GROUPACTIVE(J1)=.TRUE. 
574:       WRITE(*,'(A,I6,A)') ' intlbfgs> All permutable atoms in group ',J1,' are active' 
575: 975   NDUMMY=NDUMMY+NPERMSIZE(J1) 
576:    ENDDO  
577:  
578: !  COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint 
579: !  COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint 
580: !  WRITE(*,'(A)') ' intlbfgs> checking alignment of endpoints - should be no permutations' 
581: !  CALL LOPERMDIST(COORDSB,COORDSC,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,.FALSE.,RMATBEST,0,NMOVE,PERM) 
582: !  WRITE(*,'(A,G20.10,A,I6)') ' intlbfgs> endpoint distance ',DISTANCE,' permutations=',NMOVE  
583:  
584:    COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint 
585:    COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint 
586:    SAVELOCALPERMCUT=LOCALPERMCUT 
587: ! 
588: ! needs to be a bit sloppier than usual to allow for distorted geometries. Otherwise we may get wrong assignments 
589: ! because we don't have enough neighbours. 
590: ! 
591:    LOCALPERMCUT=QCIPERMCUT ! 0.8D0  
592:    DO J3=1,INTIMAGE 
593:       np: DO J1=1,NPERMGROUP 
594:          IF (.NOT.GROUPACTIVE(J1)) CYCLE np 
595:          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint 
596:          CALL LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCES,DIST2,.FALSE.,RMATBEST,J1,NMOVES,PERMS) 
597:          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint 
598:          CALL LOPERMDIST(COORDSC,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCEF,DIST2,.FALSE.,RMATBEST,J1,NMOVEF,PERMF) 
599:          IF ((NMOVES.GT.0).AND.(NMOVEF.GT.0)) THEN 
600:             WRITE(*,'(A,I6,A,I6,A,2G20.10,A,2I6)') ' intlbfgs> image ',J3+1,' group ',J1,' start and finish distances =', & 
601:   &                                       DISTANCES,DISTANCEF,' permutations=',NMOVES, NMOVEF 
602:          ENDIF 
603:          IF ((NMOVES.GT.0).AND.(NMOVES.EQ.NMOVEF)) THEN 
604:             COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))  
605:             DO J2=1,NATOMS 
606:                IF (PERMS(J2).NE.J2) THEN 
607:                   IF (PERMF(J2).EQ.PERMS(J2)) THEN 
608:                      WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, & 
609:   &                                           ' consistent non-identity permutation for start and finish, move atom ', &   
610:   &                                           PERMS(J2),' to position ',J2   
611:                      COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+1) 
612:                      COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+2) 
613:                      COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+3) 
614:                   ELSE 
615:                      WRITE(*,'(A,I6,A,I6)') ' intlbfgs> inconsistent non-identity permutations for start and finish' 
616:                   ENDIF 
617:                ENDIF 
618:             ENDDO 
619:             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS) 
620:          ENDIF 
621:       ENDDO np 
622:    ENDDO 
623:    LOCALPERMCUT=SAVELOCALPERMCUT 
624:    IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER) 
625: !  STOP 
626: ENDIF 
627:  
628:  
629: !487: !
630: !  Add next atom to active set if ADDATOM is true. 488: !  Add next atom to active set if ADDATOM is true. 
631: !  Constraints to atoms already in the active set are turned on489: !  Constraints to atoms already in the active set are turned on
632: !  and short-range repulsions to active atoms that are not distance constrained are turned on.490: !  and short-range repulsions to active atoms that are not distance constrained are turned on.
633: !  *** OLD Find nearest atom to active set attached by a constraint491: !  *** OLD Find nearest atom to active set attached by a constraint
634: !  *** NEW Find atom with most constraints to active set492: !  *** NEW Find atom with most constraints to active set
635: !  Turn on constraint terms for this atom with all previous members of the active set493: !  Turn on constraint terms for this atom with all previous members of the active set
636: !  Add repulsions to non-constrained atoms in this set494: !  Add repulsions to non-constrained atoms in this set
637: !  NTOADD is the number of atoms to add to the active set in each pass. 1 seems best!495: !  NTOADD is the number of atoms to add to the active set in each pass. 1 seems best!
638: !496: !
773:       IF (DEBUG) PRINT '(2(A,I6))', ' intlbfgs> Number of frozen images=',NIMAGEFREEZE,' / ',INTIMAGE631:       IF (DEBUG) PRINT '(2(A,I6))', ' intlbfgs> Number of frozen images=',NIMAGEFREEZE,' / ',INTIMAGE
774:    ENDIF632:    ENDIF
775:    !  We now have the proposed step - update geometry and calculate new gradient633:    !  We now have the proposed step - update geometry and calculate new gradient
776:    NDECREASE=0634:    NDECREASE=0
777: 20 X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)635: 20 X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
778: 636: 
779: !  IF (.NOT.SWITCHED) THEN637: !  IF (.NOT.SWITCHED) THEN
780:    IF (.TRUE.) THEN638:    IF (.TRUE.) THEN
781: !     IF ((RMS.LT.INTRMSTOL*1.0D10).AND.(MOD(NITERDONE,10).EQ.0).AND.(NSTEPSMAX-NITERDONE.GT.100)) &639: !     IF ((RMS.LT.INTRMSTOL*1.0D10).AND.(MOD(NITERDONE,10).EQ.0).AND.(NSTEPSMAX-NITERDONE.GT.100)) &
782: ! &               CALL CHECKSEP(NMAXINT,NMININT,INTIMAGE,XYZ,(3*NATOMS),NATOMS)640: ! &               CALL CHECKSEP(NMAXINT,NMININT,INTIMAGE,XYZ,(3*NATOMS),NATOMS)
783:       IF (REMOVEIMAGE.OR.(MOD(NITERDONE,INTIMAGECHECK).EQ.0)) THEN641:       IF (MOD(NITERDONE,INTIMAGECHECK).EQ.0) THEN
784: 864      CONTINUE ! for adding more than one image at a time642: 864      CONTINUE ! for adding more than one image at a time
785:          DMAX=-1.0D0643:          DMAX=-1.0D0
786:          ADMAX=-1.0D0644:          ADMAX=-1.0D0
787:          DMIN=HUGE(1.0D0)645:          DMIN=HUGE(1.0D0)
788:          DO J1=1,INTIMAGE+1646:          DO J1=1,INTIMAGE+1
789:             DUMMY=0.0D0647:             DUMMY=0.0D0
790: !           DO J2=1,3*NATOMS648: !           DO J2=1,3*NATOMS
791: !              IF (ATOMACTIVE((J2-1)/3+1)) THEN649: !              IF (ATOMACTIVE((J2-1)/3+1)) THEN
792: !                 DUMMY=DUMMY+( XYZ((J1-1)*3*NATOMS+J2) - XYZ(J1*3*NATOMS+J2) )**2650: !                 DUMMY=DUMMY+( XYZ((J1-1)*3*NATOMS+J2) - XYZ(J1*3*NATOMS+J2) )**2
793: !              ENDIF651: !              ENDIF
820: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1678: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1
821:          ENDDO679:          ENDDO
822: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &680: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &
823: ! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE681: ! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE
824: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &682: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &
825: ! &                                                  DMAX,' for images ',JMAX,JMAX+1683: ! &                                                  DMAX,' for images ',JMAX,JMAX+1
826: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> smallest image separation is ', &684: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> smallest image separation is ', &
827: ! &                                                  DMIN,' for images ',JMIN,JMIN+1685: ! &                                                  DMIN,' for images ',JMIN,JMIN+1
828: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,G20.10)') ' intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)686: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,G20.10)') ' intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)
829: !        IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN687: !        IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN
830:          IF ((.NOT.REMOVEIMAGE).AND.((SQRT(ADMAX).GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE))) THEN688:          IF ((SQRT(ADMAX).GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN
831:             JMAX=JA1689:             JMAX=JA1
832:             WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> Add an image between ',JMAX,' and ',JMAX+1,' INTIMAGE=',INTIMAGE690:             WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> Add an image between ',JMAX,' and ',JMAX+1,' INTIMAGE=',INTIMAGE
833:             NITERUSE=0691:             NITERUSE=0
834:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))692:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
835:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))693:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
836:             DEALLOCATE(XYZ)694:             DEALLOCATE(XYZ)
837:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))695:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))
838:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)696:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)
839:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &697:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &
840:   &                                               + DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1)))/2.0D0698:   &                                               + DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1)))/2.0D0
904:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)762:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
905:             IF (QCIADDREP.GT.0) THEN763:             IF (QCIADDREP.GT.0) THEN
906:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)764:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
907:             ELSEIF (CHECKCONINT) THEN765:             ELSEIF (CHECKCONINT) THEN
908:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)766:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
909:             ELSE767:             ELSE
910:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)768:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
911:             ENDIF769:             ENDIF
912: !           GOTO 864770: !           GOTO 864
913:          ENDIF771:          ENDIF
914:          IF (REMOVEIMAGE.OR.((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1))) THEN772:          IF ((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1)) THEN
915:             IF (REMOVEIMAGE) JMIN=JMAXEEE 
916:             IF (JMIN.EQ.1) JMIN=2773:             IF (JMIN.EQ.1) JMIN=2
917:             WRITE(*,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN774:             WRITE(*,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN
918:             NITERUSE=0775:             NITERUSE=0
919:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))776:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
920:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))777:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
921:             DEALLOCATE(XYZ)778:             DEALLOCATE(XYZ)
922:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+1)))779:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+1)))
923:             XYZ(1:3*NATOMS*(JMIN-1))=DPTMP(1:3*NATOMS*(JMIN-1))780:             XYZ(1:3*NATOMS*(JMIN-1))=DPTMP(1:3*NATOMS*(JMIN-1))
924:             XYZ(3*NATOMS*(JMIN-1)+1:3*NATOMS*(INTIMAGE+1))=DPTMP(3*NATOMS*JMIN+1:3*NATOMS*(INTIMAGE+2))781:             XYZ(3*NATOMS*(JMIN-1)+1:3*NATOMS*(INTIMAGE+1))=DPTMP(3*NATOMS*JMIN+1:3*NATOMS*(INTIMAGE+2))
925: 782: 
982:             INTIMAGE=INTIMAGE-1839:             INTIMAGE=INTIMAGE-1
983:             D=(3*NATOMS)*INTIMAGE840:             D=(3*NATOMS)*INTIMAGE
984:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)841:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
985:             IF (QCIADDREP.GT.0) THEN842:             IF (QCIADDREP.GT.0) THEN
986:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)843:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
987:             ELSEIF (CHECKCONINT) THEN844:             ELSEIF (CHECKCONINT) THEN
988:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)845:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
989:             ELSE846:             ELSE
990:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)847:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
991:             ENDIF848:             ENDIF
992:             NLASTGOODE=NITERDONE 
993:             LASTGOODE=ETOTAL 
994: !           GOTO 864849: !           GOTO 864
995:          ENDIF850:          ENDIF
996:       ELSE851:       ELSE
997:          DMAX=-1.0D0852:          DMAX=-1.0D0
998:          ADMAX=-1.0D0853:          ADMAX=-1.0D0
999:          DMIN=HUGE(1.0D0)854:          DMIN=HUGE(1.0D0)
1000:          DUMMY2=0.0D0855:          DUMMY2=0.0D0
1001:          DO J1=1,INTIMAGE+1856:          DO J1=1,INTIMAGE+1
1002:             DUMMY=0.0D0857:             DUMMY=0.0D0
1003: !           DO J2=1,3*NATOMS858: !           DO J2=1,3*NATOMS
1042: !           KINT=MIN(1.0D6,KINT*1.1D0)897: !           KINT=MIN(1.0D6,KINT*1.1D0)
1043: !        ELSE898: !        ELSE
1044: !           KINT=MAX(1.0D-6,KINT/1.1D0)899: !           KINT=MAX(1.0D-6,KINT/1.1D0)
1045: !        ENDIF900: !        ENDIF
1046: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT901: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT
1047:       ENDIF902:       ENDIF
1048:    ENDIF903:    ENDIF
1049: !904: !
1050: ! End of add/subtract images block.905: ! End of add/subtract images block.
1051: !906: !
1052: !907:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN
1053: ! The new QCILPERMDIST check should be used instead of these lines908:       LDEBUG=.FALSE.
1054: !909:       DO J2=2,INTIMAGE+2
1055: !  IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN910:          CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
1056: !     LDEBUG=.FALSE.911:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
1057: !     DO J2=2,INTIMAGE+2912:       ENDDO
1058: !        CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &913:    ENDIF
1059: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT) 
1060: !     ENDDO 
1061: !  ENDIF 
1062: 914: 
1063:    IF (.NOT.SWITCHED) THEN915:    IF (.NOT.SWITCHED) THEN
1064:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)916:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1065:       IF (QCIADDREP.GT.0) THEN917:       IF (QCIADDREP.GT.0) THEN
1066:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)918:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1067:       ELSEIF (CHECKCONINT) THEN919:       ELSEIF (CHECKCONINT) THEN
1068:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)920:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1069:       ELSE921:       ELSE
1070:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)922:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1071:       ENDIF923:       ENDIF
1134:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)986:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
1135:       QCIIMAGE=INTIMAGE987:       QCIIMAGE=INTIMAGE
1136:       INTIMAGE=INTIMAGESAVE988:       INTIMAGE=INTIMAGESAVE
1137:       RETURN989:       RETURN
1138:    ENDIF990:    ENDIF
1139: 991: 
1140:    STEPTOT = SUM(STEPIMAGE)/INTIMAGE992:    STEPTOT = SUM(STEPIMAGE)/INTIMAGE
1141: 993: 
1142:    MAXRMS=-1.0D0994:    MAXRMS=-1.0D0
1143:    MAXEEE=-1.0D100995:    MAXEEE=-1.0D100
1144:    MINEEE=1.0D100 
1145:    SUMEEE=0.0D0 
1146:    SUMEEE2=0.0D0 
1147:    DO J1=2,INTIMAGE+1996:    DO J1=2,INTIMAGE+1
1148:       SUMEEE=SUMEEE+EEE(J1) 
1149:       SUMEEE2=SUMEEE2+EEE(J1)**2 
1150:       IF (EEE(J1).GT.MAXEEE) THEN997:       IF (EEE(J1).GT.MAXEEE) THEN
1151:          MAXEEE=EEE(J1)998:          MAXEEE=EEE(J1)
1152:          JMAXEEE=J1999:          JMAXEEE=J1
1153:       ENDIF1000:       ENDIF
1154:       IF (EEE(J1).LT.MINEEE) THEN 
1155:          MINEEE=EEE(J1) 
1156:       ENDIF 
1157:       DUMMY=0.0D01001:       DUMMY=0.0D0
1158:       DO J2=1,3*NATOMS1002:       DO J2=1,3*NATOMS
1159:          DUMMY=DUMMY+GGG(3*NATOMS*(J1-1)+J2)**21003:          DUMMY=DUMMY+GGG(3*NATOMS*(J1-1)+J2)**2
1160:       ENDDO1004:       ENDDO
1161:       IF (DUMMY.GT.MAXRMS) THEN1005:       IF (DUMMY.GT.MAXRMS) THEN
1162:          MAXRMS=DUMMY1006:          MAXRMS=DUMMY
1163:          JMAXRMS=J11007:          JMAXRMS=J1
1164:       ENDIF1008:       ENDIF
1165:    ENDDO1009:    ENDDO
1166:    MAXRMS=SQRT(MAXRMS/(3*NACTIVE))1010:    MAXRMS=SQRT(MAXRMS/(3*NACTIVE))
1167:    SUMEEE=SUMEEE/INTIMAGE 
1168:    SUMEEE2=SUMEEE2+EEE(J1)**2 
1169:    SUMEEE2=SUMEEE2/INTIMAGE 
1170:    SIGMAEEE=SQRT(SUMEEE2-SUMEEE**2) 
1171:    REMOVEIMAGE=.FALSE. 
1172: !  IF (ABS(MAXEEE-SUMEEE).GT.3.0D0*SIGMAEEE) THEN 
1173:    IF (MAXEEE.GT.1.0D2*MINEEE) THEN 
1174:       WRITE(*,'(A,I8,A,G20.10,A)') ' intlbfgs> warning *** highest image ',JMAXEEE,' energy ',MAXEEE,' is > 100 * higher than minimum '   
1175: !     REMOVEIMAGE=.TRUE. 
1176: !     IF (NITERDONE-NLASTGOODE.LT.40) REMOVEIMAGE=.FALSE. 
1177:    ENDIF 
1178: 1011: 
1179:    IF (DEBUG) THEN1012:    IF (DEBUG) THEN
1180:       WRITE(*,'(A,2G20.10)') ' intlbfgs> mean and sigma of image energies=',SUMEEE,SIGMAEEE 
1181:       WRITE(*,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &1013:       WRITE(*,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &
1182:   &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS1014:   &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS
1183:       CALL FLUSH(6)1015:       CALL FLUSH(6)
1184:    ENDIF1016:    ENDIF
1185: 1017: 
1186:    IF (.NOT.SWITCHED) THEN1018:    IF (.NOT.SWITCHED) THEN
1187: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN1019: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN
1188:       IF (.FALSE.) THEN ! no backtracking1020:       IF (.FALSE.) THEN ! no backtracking
1189:          WRITE(*,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE1021:          WRITE(*,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE
1190:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+11022:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+1
1303:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1135:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1304: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771136: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1305:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771137:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1306:          IF (NACTIVE.LT.NATOMS) THEN 1138:          IF (NACTIVE.LT.NATOMS) THEN 
1307:             ADDATOM=.TRUE.1139:             ADDATOM=.TRUE.
1308:             GOTO 7771140:             GOTO 777
1309:          ENDIF1141:          ENDIF
1310:          CALL MYCPU_TIME(FTIME,.FALSE.)1142:          CALL MYCPU_TIME(FTIME,.FALSE.)
1311:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1143:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1312:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1144:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1313:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)1145:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)
1314:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1146:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1315:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1147:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
1316:          DO J1=1,NATOMS1148:          DO J1=1,NATOMS
1317:             IF (.NOT.ATOMACTIVE(J1)) THEN1149:             IF (.NOT.ATOMACTIVE(J1)) THEN
1318:                WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> ERROR *** number of active atoms=',NACTIVE,' but atom ',J1,' is not active'1150:                WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> ERROR *** number of active atoms=',NACTIVE,' but atom ',J1,' is not active'
1319:             ENDIF1151:             ENDIF
1320:          ENDDO1152:          ENDDO
1321:          NSTEPSMAX=NITERDONE+INTCONSTEPS1153:          NSTEPSMAX=NITERDONE+INTCONSTEPS
1322:          SWITCHED=.TRUE.1154:          SWITCHED=.TRUE.
1323:          RMS=INTRMSTOL*10.0D0 ! to prevent premature convergence1155:          RMS=INTRMSTOL*10.0D0 ! to prevent premature convergence
1338:    777 CONTINUE1170:    777 CONTINUE
1339: !1171: !
1340: ! Compute the new step and gradient change1172: ! Compute the new step and gradient change
1341: !1173: !
1342:    NPT=POINT*D1174:    NPT=POINT*D
1343:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)1175:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)
1344:    GDIF(POINT,:)=G-GTMP1176:    GDIF(POINT,:)=G-GTMP
1345:    1177:    
1346:    POINT=POINT+1; IF (POINT==MUPDATE) POINT=01178:    POINT=POINT+1; IF (POINT==MUPDATE) POINT=0
1347: 1179: 
1348:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)1180:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)
1349:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1181:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1350: 1182: 
1351:    NITERDONE=NITERDONE+11183:    NITERDONE=NITERDONE+1
1352:    NITERUSE=NITERUSE+11184:    NITERUSE=NITERUSE+1
1353: 1185: 
1354:    IF (NITERDONE.GT.NSTEPSMAX) EXIT1186:    IF (NITERDONE.GT.NSTEPSMAX) EXIT
1355:    IF (NACTIVE.EQ.NATOMS) THEN1187:    IF (NACTIVE.EQ.NATOMS) THEN
1356:       IF (.NOT.SWITCHED) THEN1188:       IF (.NOT.SWITCHED) THEN
1357:          CALL MYCPU_TIME(FTIME,.FALSE.)1189:          CALL MYCPU_TIME(FTIME,.FALSE.)
1358:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1190:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1382: ENDIF1214: ENDIF
1383: IF (EXITSTATUS.EQ.1) THEN1215: IF (EXITSTATUS.EQ.1) THEN
1384:    WRITE(*,'(A,I6,A,G20.10,A,G15.8,A,I4)') ' intlbfgs> Converged after ',NITERDONE,' steps, energy/image=',ETOTAL/INTIMAGE, &1216:    WRITE(*,'(A,I6,A,G20.10,A,G15.8,A,I4)') ' intlbfgs> Converged after ',NITERDONE,' steps, energy/image=',ETOTAL/INTIMAGE, &
1385:   &                               ' RMS=',RMS,' images=',INTIMAGE1217:   &                               ' RMS=',RMS,' images=',INTIMAGE
1386: ELSEIF (EXITSTATUS.EQ.2) THEN1218: ELSEIF (EXITSTATUS.EQ.2) THEN
1387:    WRITE(*,'(A,I6,A,G20.10,A,G15.8,A,I4)') ' intlbfgs> After ',NITERDONE,' steps, energy/image=',ETOTAL/INTIMAGE, &1219:    WRITE(*,'(A,I6,A,G20.10,A,G15.8,A,I4)') ' intlbfgs> After ',NITERDONE,' steps, energy/image=',ETOTAL/INTIMAGE, &
1388:   &                               ' RMS=',RMS,' images=',INTIMAGE1220:   &                               ' RMS=',RMS,' images=',INTIMAGE
1389: ENDIF1221: ENDIF
1390: 678 CONTINUE1222: 678 CONTINUE
1391: 1223: 
1392: ! CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)1224: CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)
1393: ! CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1225: CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1394: 1226: 
1395: IF (DEBUG) WRITE(*,'(A,G20.10)') 'intlbfgs> WORST=',WORST1227: IF (DEBUG) WRITE(*,'(A,G20.10)') 'intlbfgs> WORST=',WORST
1396: 1228: 
1397: BESTWORST=WORST1229: BESTWORST=WORST
1398: BESTINTIMAGE=INTIMAGE1230: BESTINTIMAGE=INTIMAGE
1399: IF (ALLOCATED(QCIXYZ)) DEALLOCATE(QCIXYZ)1231: IF (ALLOCATED(QCIXYZ)) DEALLOCATE(QCIXYZ)
1400: ALLOCATE(QCIXYZ(3*NATOMS*(INTIMAGE+2)))1232: ALLOCATE(QCIXYZ(3*NATOMS*(INTIMAGE+2)))
1401: QCIXYZ(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))1233: QCIXYZ(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
1402: WRITE(*,'(A,I8,A,G20.10)') 'intlbfgs> retaining ',INTIMAGE,' QCI images, highest energy=',BESTWORST1234: WRITE(*,'(A,I8,A,G20.10)') 'intlbfgs> retaining ',INTIMAGE,' QCI images, highest energy=',BESTWORST
1403: 1235: 
1404: CALL INTRWG(NACTIVE,0,INTIMAGE,XYZ,TURNONORDER)1236: CALL INTRWG(NACTIVE,0,INTIMAGE,XYZ)
1405: CALL WRITEPROFILE(0,EEE,INTIMAGE)1237: CALL WRITEPROFILE(0,EEE,INTIMAGE)
1406: 1238: 
1407: DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)1239: DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)
1408: DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &1240: DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
1409:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)1241:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
1410: QCIIMAGE=INTIMAGE1242: QCIIMAGE=INTIMAGE
1411: INTIMAGE=INTIMAGESAVE1243: INTIMAGE=INTIMAGESAVE
1412: 1244: 
1413: END SUBROUTINE INTLBFGS1245: END SUBROUTINE INTLBFGS
1414: !1246: !
1468:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)1300:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)
1469:          GOTO 2461301:          GOTO 246
1470:       ENDIF1302:       ENDIF
1471:    ENDDO 1303:    ENDDO 
1472: 246 CONTINUE1304: 246 CONTINUE
1473: ENDDO1305: ENDDO
1474: IF (DEBUG) WRITE(*,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE1306: IF (DEBUG) WRITE(*,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE
1475: 1307: 
1476: END SUBROUTINE CHECKREP1308: END SUBROUTINE CHECKREP
1477: 1309: 
1478: SUBROUTINE INTRWG(NACTIVE,NITER,INTIMAGE,XYZ,TURNONORDER)1310: SUBROUTINE INTRWG(NACTIVE,NITER,INTIMAGE,XYZ)
1479: USE PORFUNCS1311: USE PORFUNCS
1480: USE KEY,ONLY: STOCKT,STOCKAAT, RBAAT, ATOMACTIVE, NCONSTRAINT, CONACTIVE, NREPULSIVE, NNREPULSIVE, REPI, REPJ, REPCUT, NREPCUT, &1312: USE KEY,ONLY: STOCKT,STOCKAAT, RBAAT, ATOMACTIVE
1481:   &           NREPMAX, NREPI, NREPJ 
1482: USE COMMONS, ONLY: NATOMS1313: USE COMMONS, ONLY: NATOMS
1483: IMPLICIT NONE1314: IMPLICIT NONE
1484: CHARACTER(LEN=10) :: XYZFILE   = 'int.xyz   '1315: CHARACTER(LEN=10) :: XYZFILE   = 'int.xyz   '
1485: CHARACTER(LEN=10) :: QCIFILE   = 'QCIdump   '1316: INTEGER,INTENT(IN) :: NITER
1486: INTEGER,INTENT(IN) :: NITER, TURNONORDER(NATOMS)1317: INTEGER :: J1,J2,INTIMAGE,J3,NACTIVE
1487: INTEGER :: J1,J2,INTIMAGE,J3,NACTIVE,LUNIT,GETUNIT 
1488: CHARACTER(LEN=80) :: FILENAME,DUMMYS1318: CHARACTER(LEN=80) :: FILENAME,DUMMYS
1489: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2))1319: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2))
1490: 1320: 
1491: FILENAME=XYZFILE1321: FILENAME=XYZFILE
1492: 1322: 
1493: ! IF (NITER.GT.0) THEN1323: IF (NITER.GT.0) THEN
1494: !    WRITE(DUMMYS,'(I8)') NITER1324:    WRITE(DUMMYS,'(I8)') NITER
1495: !    FILENAME='int.' // TRIM(ADJUSTL(DUMMYS)) // '.xyz' ! so that vmd recognises the file type!1325:    FILENAME='int.' // TRIM(ADJUSTL(DUMMYS)) // '.xyz' ! so that vmd recognises the file type!
1496: ! ENDIF1326: ENDIF
1497: LUNIT=GETUNIT()1327: OPEN(UNIT=993,FILE=FILENAME,STATUS='replace')
1498: OPEN(UNIT=LUNIT,FILE=TRIM(ADJUSTL(FILENAME)),STATUS='replace') 
1499: DO J2=1,INTIMAGE+21328: DO J2=1,INTIMAGE+2
1500: !  WRITE(LUNIT,'(i4/)') NACTIVE1329: !  WRITE(993,'(i4/)') NACTIVE
1501:    WRITE(LUNIT,'(i4/)') NATOMS1330:    WRITE(993,'(i4/)') NATOMS
1502:    DO J3=1,NATOMS1331:    DO J3=1,NATOMS
1503:       IF (ATOMACTIVE(J3)) THEN1332:       IF (ATOMACTIVE(J3)) THEN
1504:          WRITE(LUNIT,'(A5,1X,3F20.10)') 'LA   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  1333:          WRITE(993,'(A5,1X,3F20.10)') 'LA   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  
1505:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  1334:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  
1506:       ELSE1335:       ELSE
1507:          WRITE(LUNIT,'(A5,1X,3F20.10)') 'DU   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  1336:          WRITE(993,'(A5,1X,3F20.10)') 'DU   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  
1508:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  1337:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  
1509:       ENDIF1338:       ENDIF
1510:    ENDDO1339:    ENDDO
1511: ENDDO1340: ENDDO
1512: 1341: 
1513: WRITE(*,*) 'rwg> Interpolated image coordinates were saved to xyz file "'//TRIM(FILENAME)//'"'1342: WRITE(*,*) 'rwg> Interpolated image coordinates were saved to xyz file "'//TRIM(FILENAME)//'"'
1514: 1343: 
1515: CLOSE(LUNIT)1344: CLOSE(UNIT=993)
1516:  
1517: FILENAME=QCIFILE 
1518: LUNIT=GETUNIT() 
1519: OPEN(UNIT=LUNIT,FILE=TRIM(ADJUSTL(FILENAME)),STATUS='replace') 
1520:  
1521: WRITE(*,'(A,I10,A)') ' intlbfgs> dumping state for ',NACTIVE,' active atoms' 
1522: WRITE(LUNIT,'(I10)') NACTIVE 
1523: WRITE(*,'(A,I10,A)') ' intlbfgs> dumping turnonorder for ',NACTIVE,' active atoms' 
1524: WRITE(LUNIT,'(12I8)') TURNONORDER(1:NACTIVE) 
1525: WRITE(*,'(A)') ' intlbfgs> dumping atomactive' 
1526: WRITE(LUNIT,'(12L5)') ATOMACTIVE(1:NATOMS) 
1527: WRITE(LUNIT,'(I10)') NCONSTRAINT 
1528: WRITE(*,'(A,I10,A)') ' intlbfgs> dumping conactive for ',NCONSTRAINT,' constraints' 
1529: WRITE(LUNIT,'(12L5)') CONACTIVE(1:NCONSTRAINT) 
1530:  
1531:    WRITE(LUNIT,'(3I12,G20.10)') NREPULSIVE,NNREPULSIVE,NREPMAX 
1532:    WRITE(*,'(A,3I10,G20.10)') 'intlbfgs> dumping NREPULSIVE,NNREPULSIVE,NREPMAX=',NREPULSIVE,NNREPULSIVE,NREPMAX 
1533:  
1534:    WRITE(LUNIT,'(12I8)') REPI(1:NREPULSIVE) 
1535:    WRITE(*,'(A)') ' intlbfgs> dumped REPI:' 
1536:    WRITE(LUNIT,'(12I8)') REPJ(1:NREPULSIVE) 
1537:    WRITE(*,'(A)') ' intlbfgs> dumped REPJ:' 
1538:    WRITE(LUNIT,'(12I8)') NREPI(1:NNREPULSIVE) 
1539:    WRITE(*,'(A)') ' intlbfgs> dumped NREPI:' 
1540:    WRITE(LUNIT,'(12I8)') NREPJ(1:NNREPULSIVE) 
1541:    WRITE(*,'(A)') ' intlbfgs> dumped NREPJ:' 
1542:  
1543:    WRITE(LUNIT,'(6G20.10)') REPCUT(1:NREPULSIVE) 
1544:    WRITE(*,'(A)') ' intlbfgs> dumped REPCUT:' 
1545:    WRITE(LUNIT,'(6G20.10)') NREPCUT(1:NNREPULSIVE) 
1546:    WRITE(*,'(A)') ' intlbfgs> dumped NREPCUT:' 
1547:  
1548: CLOSE(LUNIT) 
1549:  
1550: END SUBROUTINE INTRWG1345: END SUBROUTINE INTRWG
1551: 1346: 
1552: SUBROUTINE WRITEPROFILE(NITER,EEE,INTIMAGE)1347: SUBROUTINE WRITEPROFILE(NITER,EEE,INTIMAGE)
1553: IMPLICIT NONE 1348: IMPLICIT NONE 
1554: INTEGER,INTENT(IN) :: NITER, INTIMAGE1349: INTEGER,INTENT(IN) :: NITER, INTIMAGE
1555: INTEGER :: I,LUNIT,GETUNIT1350: INTEGER :: I,UNIT
1556: DOUBLE PRECISION :: EEE(INTIMAGE+2)1351: DOUBLE PRECISION :: EEE(INTIMAGE+2)
1557: CHARACTER(LEN=20) :: FILENAME1352: CHARACTER(LEN=20) :: FILENAME
1558: 1353: 
1559: LUNIT=GETUNIT()1354: UNIT=992
1560: ! IF (NITER.GT.0) THEN1355: IF (NITER.GT.0) THEN
1561: !    WRITE(FILENAME,'(I8)') NITER1356:    WRITE(FILENAME,'(I8)') NITER
1562: !    FILENAME='int.EofS.' // TRIM(ADJUSTL(FILENAME))1357:    FILENAME='int.EofS.' // TRIM(ADJUSTL(FILENAME))
1563: ! ELSE   1358: ELSE   
1564:    FILENAME='int.EofS'1359:    FILENAME='int.EofS'
1565: ! ENDIF1360: ENDIF
1566: OPEN(UNIT=LUNIT,FILE=FILENAME,STATUS='replace')1361: OPEN(UNIT=UNIT,FILE=FILENAME,STATUS='replace')
1567: 1362: 
1568: WRITE(UNIT=LUNIT,FMT='(2g24.13)') EEE(1)1363: WRITE(UNIT=UNIT,FMT='(2g24.13)') EEE(1)
1569: DO I=2,INTIMAGE+11364: DO I=2,INTIMAGE+1
1570:    WRITE(UNIT=LUNIT,FMT='(2G24.13)') EEE(I)1365:    WRITE(UNIT=UNIT,FMT='(2G24.13)') EEE(I)
1571: ENDDO1366: ENDDO
1572: WRITE(UNIT=LUNIT,FMT='(2G24.13)') EEE(INTIMAGE+2)1367: WRITE(UNIT=UNIT,FMT='(2G24.13)') EEE(INTIMAGE+2)
1573: 1368: 
1574: CLOSE(LUNIT)1369: CLOSE(UNIT)
1575: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'1370: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'
1576: 1371: 
1577: END SUBROUTINE WRITEPROFILE1372: END SUBROUTINE WRITEPROFILE
1578: 1373: 
1579: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)1374: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)
1580: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &1375: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &
1581:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &1376:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &
1582:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, &1377:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, &
1583:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &1378:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &
1584:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, QCIADDREP1379:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, QCIADDREP
1690:             IF (.NOT.ATOMACTIVE(CONJ(J1))) THEN1485:             IF (.NOT.ATOMACTIVE(CONJ(J1))) THEN
1691:                IF (CONDISTREF(J1).LT.DUMMY2) THEN1486:                IF (CONDISTREF(J1).LT.DUMMY2) THEN
1692:                   DUMMY2=CONDISTREF(J1)1487:                   DUMMY2=CONDISTREF(J1)
1693:                   NEWATOM=CONJ(J1)1488:                   NEWATOM=CONJ(J1)
1694:                ENDIF1489:                ENDIF
1695:             ENDIF1490:             ENDIF
1696:          ENDIF1491:          ENDIF
1697:       ENDDO1492:       ENDDO
1698:       IF (DEBUG) WRITE(*,'(3(A,I6),A,F15.5)') ' intlbfgs> Choosing new active atom ',NEWATOM,' new constraints=', &1493:       IF (DEBUG) WRITE(*,'(3(A,I6),A,F15.5)') ' intlbfgs> Choosing new active atom ',NEWATOM,' new constraints=', &
1699:   &                                       NCONTOACTIVE(NEWATOM),' maximum=',NBEST,' shortest constraint=',DUMMY21494:   &                                       NCONTOACTIVE(NEWATOM),' maximum=',NBEST,' shortest constraint=',DUMMY2
1700: !     IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,ETOTAL)1495: !     IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)
1701: !     IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1496: !     IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1702: 1497: 
1703:           1498:           
1704:       IF (NEWATOM*NBEST.EQ.0) THEN ! sanity check1499:       IF (NEWATOM*NBEST.EQ.0) THEN ! sanity check
1705:          WRITE(*,'(A,I6,A,2I6)') ' intlbfgs> ERROR *** new active atom not set'1500:          WRITE(*,'(A,I6,A,2I6)') ' intlbfgs> ERROR *** new active atom not set'
1706:          STOP1501:          STOP
1707:       ELSE1502:       ELSE
1708: !1503: !
1709: !  We need a sorted list of up to 3 active atoms, sorted according to how well the1504: !  We need a sorted list of up to 3 active atoms, sorted according to how well the
1710: !  end point distance is preserved, even if they don't satisfy the constraint 1505: !  end point distance is preserved, even if they don't satisfy the constraint 


r33340/key.f90 2017-09-23 12:30:16.096999786 +0100 r33339/key.f90 2017-09-23 12:30:17.725021435 +0100
 18:      &        NTRAPPOW, MAXINTIMAGE, CHECKDID, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, & 18:      &        NTRAPPOW, MAXINTIMAGE, CHECKDID, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, &
 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, & 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, &
 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, & 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, &
 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, & 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, &
 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, & 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, &
 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, & 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, &
 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, & 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, &
 25:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS, & 25:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS, &
 26:      &        MLQIN, MLQSTART, MLQOUT, MLQDATA, NMLQ, & 26:      &        MLQIN, MLQSTART, MLQOUT, MLQDATA, NMLQ, &
 27:      &        QCIADDREP, QCIBONDS, QCISECOND, MAXNACTIVE, QCIIMAGE, NADDTARGET, NUMNN, MULTI_COUNT, MULTI_LAST, MULTI_STEP, & 27:      &        QCIADDREP, QCIBONDS, QCISECOND, MAXNACTIVE, QCIIMAGE, NADDTARGET, NUMNN, MULTI_COUNT, MULTI_LAST, MULTI_STEP, &
 28:      &        NDOF, RECCOUNT, MLPPROBPOS, PUSHOPTMAX, MLPNEIGH, QCICYCN, QCIPDINT 28:      &        NDOF, RECCOUNT, MLPPROBPOS, PUSHOPTMAX, MLPNEIGH, QCICYCN
 29:  29: 
 30:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, & 30:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, &
 31:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, & 31:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, &
 32:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, & 32:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, &
 33:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, & 33:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, &
 34:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, & 34:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, &
 35:      &        BFGSSTEP, EFOLSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, & 35:      &        BFGSSTEP, EFOLSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, &
 36:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, & 36:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, &
 37:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, & 37:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, &
 38:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,& 38:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,&
 52:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 52:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 53:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 53:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 54:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 54:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 55:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 55:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 56:      &        PBST, SSHT, GAUSSIAN03, GAUSSIAN09, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 56:      &        PBST, SSHT, GAUSSIAN03, GAUSSIAN09, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 57:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 57:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 58:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 58:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 59:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 59:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 60:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T, & 60:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T, &
 61:      &        MBPOLT, MULTIJOB_MACHINET, DUMPDATA_MACHINET, PLUSSIDET, MINUSSIDET, PUSHOPTT, MLPVB3NNT, GAUSSIAN16, QCICYCLEST, & 61:      &        MBPOLT, MULTIJOB_MACHINET, DUMPDATA_MACHINET, PLUSSIDET, MINUSSIDET, PUSHOPTT, MLPVB3NNT, GAUSSIAN16, QCICYCLEST, &
 62:      &        QCIDNEBT, QCIRESTART, QCILPERMDIST 62:      &        QCIDNEBT
 63:  
 64:  63: 
 65: ! sy349 > for testing the flatpath after dneb 64: ! sy349 > for testing the flatpath after dneb
 66:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 65:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 67:       LOGICAL FLATPATHT 66:       LOGICAL FLATPATHT
 68:  67: 
 69: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 68: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 70:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 69:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 71:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 70:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 72:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 71:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 73:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 72:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
105:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &104:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &
106:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &105:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &
107:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &106:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &
108:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &107:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &
109:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &108:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &
110:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &109:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &
111:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &110:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
112:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &111:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
113:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &112:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
114:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2, &113:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2, &
115:      &        TANHFAC, LJADDCUTOFF,LJADDREFNORM,MAXIMFACTOR, PUSHOPTCONV, QCICYCDIST, QCIPERMCUT114:      &        TANHFAC, LJADDCUTOFF,LJADDREFNORM,MAXIMFACTOR, PUSHOPTCONV, QCICYCDIST
116: 115: 
117: !     sf344116: !     sf344
118:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &117:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
119:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &118:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
120:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT119:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT
121:  120:  
122:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)121:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)
123:       DOUBLE PRECISION, ALLOCATABLE :: SHIFTL(:)122:       DOUBLE PRECISION, ALLOCATABLE :: SHIFTL(:)
124:       LOGICAL, ALLOCATABLE :: uniaxarray(:)123:       LOGICAL, ALLOCATABLE :: uniaxarray(:)
125:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)124:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)


r33340/keywords.f 2017-09-23 12:30:16.325002818 +0100 r33339/keywords.f 2017-09-23 12:30:17.965024625 +0100
  1: ! OPTIM: A prograa for optimizing geometries and calculating reaction pathways  1: ! OPTIM: A program for optimizing geometries and calculating reaction pathways
  2: ! Copyright (C) 1999-2006 David J. Wales  2: ! Copyright (C) 1999-2006 David J. Wales
  3: ! This file is part of OPTIM.  3: ! This file is part of OPTIM.
  4: !   4: ! 
  5: ! OPTIM is free software; you can redistribute it and/or modify  5: ! OPTIM is free software; you can redistribute it and/or modify
  6: ! it under the terms of the GNU General Public License as published by  6: ! it under the terms of the GNU General Public License as published by
  7: ! the Free Software Foundation; either version 2 of the License, or  7: ! the Free Software Foundation; either version 2 of the License, or
  8: ! at your option) any later version.  8: ! at your option) any later version.
  9: !   9: ! 
 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
701:          QCIAMBERT=.FALSE.701:          QCIAMBERT=.FALSE.
702:          INTMINT=.FALSE.702:          INTMINT=.FALSE.
703:          INTSPRINGACTIVET=.TRUE.703:          INTSPRINGACTIVET=.TRUE.
704:          INTMINFAC=1.0D0704:          INTMINFAC=1.0D0
705:          QCIRADSHIFTT=.FALSE.705:          QCIRADSHIFTT=.FALSE.
706:          QCIRADSHIFT=1.0D0706:          QCIRADSHIFT=1.0D0
707:          QCICYCLEST=.FALSE.707:          QCICYCLEST=.FALSE.
708:          QCICYCDIST=0.0708:          QCICYCDIST=0.0
709:          QCICYCN=100709:          QCICYCN=100
710:          QCIDNEBT=.FALSE.710:          QCIDNEBT=.FALSE.
711:          QCIRESTART=.FALSE. 
712:          QCILPERMDIST=.FALSE. 
713:          QCIPDINT=1000 
714:          QCIPERMCUT=0.8D0 
715: 711: 
716:          CONPOTT=.FALSE.712:          CONPOTT=.FALSE.
717:          CPCONSTRAINTTOL=0.1D0713:          CPCONSTRAINTTOL=0.1D0
718:          CPCONSTRAINTDEL=1.0D5714:          CPCONSTRAINTDEL=1.0D5
719:          CPCONSTRAINTREP=1.0D0715:          CPCONSTRAINTREP=1.0D0
720:          CPCONSTRAINREPCUT=20.0D0716:          CPCONSTRAINREPCUT=20.0D0
721:          CPCONFRAC=1.0D-4717:          CPCONFRAC=1.0D-4
722:          CPREPSEP=0718:          CPREPSEP=0
723:          CPCONSEP=10000719:          CPCONSEP=10000
724:          CHECKOVERLAPT=.FALSE.720:          CHECKOVERLAPT=.FALSE.
1197: ! Store coordinates in Q1193: ! Store coordinates in Q
1198:             IF (NITEMS .GT. 1) THEN1194:             IF (NITEMS .GT. 1) THEN
1199:                CALL READA(AMBERSTR)1195:                CALL READA(AMBERSTR)
1200:                IF (FILTH2 .NE. 0) THEN1196:                IF (FILTH2 .NE. 0) THEN
1201:                   WRITE(OTEMP, *) FILTH21197:                   WRITE(OTEMP, *) FILTH2
1202:                   WRITE(OSTRING,'(A)') TRIM(ADJUSTL(AMBERSTR))//'.'//TRIM(ADJUSTL(OTEMP))1198:                   WRITE(OSTRING,'(A)') TRIM(ADJUSTL(AMBERSTR))//'.'//TRIM(ADJUSTL(OTEMP))
1203:                   WRITE(*,*) 'ostring=', OSTRING1199:                   WRITE(*,*) 'ostring=', OSTRING
1204:                ELSE1200:                ELSE
1205:                   WRITE(OSTRING,'(A)') TRIM(ADJUSTL(AMBERSTR))1201:                   WRITE(OSTRING,'(A)') TRIM(ADJUSTL(AMBERSTR))
1206:                END IF1202:                END IF
1207:                WRITE(*,'(A,A)') ' keywords> input coordinates for AMBER12 system will be read from ', OSTRING1203:                WRITE(*,'(A)') ' keywords> input coordinates for AMBER12 system will be read from ', OSTRING
1208:                OPEN(UNIT=3827, FILE=TRIM(ADJUSTL(OSTRING)), STATUS='UNKNOWN')1204:                OPEN(UNIT=3827, FILE=TRIM(ADJUSTL(OSTRING)), STATUS='UNKNOWN')
1209:                IF (NITEMS == 3) THEN1205:                IF (NITEMS == 3) THEN
1210:                   CALL READA(AMBERSTR)1206:                   CALL READA(AMBERSTR)
1211:                   IF (TRIM(ADJUSTL(AMBERSTR)) == 'inpcrd') THEN1207:                   IF (TRIM(ADJUSTL(AMBERSTR)) == 'inpcrd') THEN
1212:                      WRITE(*,'(A)') ' keywords> reading in AMBER restart format'1208:                      WRITE(*,'(A)') ' keywords> reading in AMBER restart format'
1213:                      CALL AMBER12_GET_COORDS(NATOMS, Q(1:3*NATOMS))1209:                      CALL AMBER12_GET_COORDS(NATOMS, Q(1:3*NATOMS))
1214:                   ELSE1210:                   ELSE
1215:                      WRITE(*,'(A)') ' keywords> reading in xyz format'1211:                      WRITE(*,'(A)') ' keywords> reading in xyz format'
1216:                      DO I = 1, NATOMS1212:                      DO I = 1, NATOMS
1217:                         READ(3827, *) Q(3*I-2:3*I)1213:                         READ(3827, *) Q(3*I-2:3*I)
2179: ! 2175: ! 
2180: ! Check for internal minimum in constraint terms for INTCONSTRAINT2176: ! Check for internal minimum in constraint terms for INTCONSTRAINT
2181: ! 2177: ! 
2182:          ELSE IF ((WORD.EQ.'CONINT').OR.(WORD.EQ.'QCIINT')) THEN2178:          ELSE IF ((WORD.EQ.'CONINT').OR.(WORD.EQ.'QCIINT')) THEN
2183:             CHECKCONINT=.TRUE.2179:             CHECKCONINT=.TRUE.
2184:             IF (NITEMS.GT.1) CALL READF(INTMINFAC)2180:             IF (NITEMS.GT.1) CALL READF(INTMINFAC)
2185:             WRITE(*,'(A,G20.10)') ' keyword> Internal minima terms will be scaled by a factor of ',INTMINFAC2181:             WRITE(*,'(A,G20.10)') ' keyword> Internal minima terms will be scaled by a factor of ',INTMINFAC
2186: !2182: !
2187: ! Maximum active atoms in QCI procedure.2183: ! Maximum active atoms in QCI procedure.
2188: !2184: !
2189:       ELSE IF (WORD.EQ.'QCIRESTART') THEN 
2190:          QCIRESTART=.TRUE. 
2191: ! 
2192: ! Maximum active atoms in QCI procedure. 
2193: ! 
2194:       ELSE IF (WORD.EQ.'QCIMAXACTIVE') THEN2185:       ELSE IF (WORD.EQ.'QCIMAXACTIVE') THEN
2195:          CALL READI(MAXNACTIVE)2186:          CALL READI(MAXNACTIVE)
2196: 2187: 
2197: 2188: 
2198:       ELSE IF (WORD.EQ.'QCICYCLES') THEN2189:       ELSE IF (WORD.EQ.'QCICYCLES') THEN
2199:          QCICYCLEST=.TRUE.2190:          QCICYCLEST=.TRUE.
2200:          QCIDNEBT=.FALSE.2191:          QCIDNEBT=.FALSE.
2201:          CALL READF(QCICYCDIST)2192:          CALL READF(QCICYCDIST)
2202:          CALL READI(QCICYCN)2193:          CALL READI(QCICYCN)
2203:          WRITE(*,'(A,F8.2,A,I8,A)') ' keyword> DNEB for distance <' , QCICYCDIST , ' and after ', QCICYCN,' cycles'   2194:          WRITE(*,'(A,F8.2,A,I8,A)') ' keyword> DNEB for distance <' , QCICYCDIST , ' and after ', QCICYCN,' cycles'   
3527: ! 3518: ! 
3528:          ELSE IF (WORD.EQ.'INDEX') THEN3519:          ELSE IF (WORD.EQ.'INDEX') THEN
3529:             CALL READI(HINDEX)3520:             CALL READI(HINDEX)
3530: !3521: !
3531: ! Radial shift to make space for new atoms.3522: ! Radial shift to make space for new atoms.
3532: !3523: !
3533:          ELSE IF (WORD.EQ.'QCIRADSHIFT') THEN3524:          ELSE IF (WORD.EQ.'QCIRADSHIFT') THEN
3534:             QCIRADSHIFTT=.TRUE.3525:             QCIRADSHIFTT=.TRUE.
3535:             IF (NITEMS.GT.1) CALL READF(QCIRADSHIFT)3526:             IF (NITEMS.GT.1) CALL READF(QCIRADSHIFT)
3536:             WRITE(*,'(A,G20.10)') ' keyword> Shifting unconstrained atoms away from added atoms by ',QCIRADSHIFT3527:             WRITE(*,'(A,G20.10)') ' keyword> Shifting unconstrained atoms away from added atoms by ',QCIRADSHIFT
3537: !3528:          ELSE IF (WORD.EQ.'QCIPERMCHECK') THEN
3538: ! Use new QCILPERMDIST check instead of this flag3529:             QCIPERMCHECK=.TRUE.
3539: !3530:             CALL READI(QCIPERMCHECKINT)
3540: !        ELSE IF (WORD.EQ.'QCIPERMCHECK') THEN 
3541: !           QCIPERMCHECK=.TRUE. 
3542: !           CALL READI(QCIPERMCHECKINT) 
3543: ! 3531: ! 
3544: ! Images for INTCONSTRAINT3532: ! Images for INTCONSTRAINT
3545: ! 3533: ! 
3546:          ELSE IF ((WORD.EQ.'INTIMAGE').OR.(WORD.EQ.'QCIIMAGE')) THEN3534:          ELSE IF ((WORD.EQ.'INTIMAGE').OR.(WORD.EQ.'QCIIMAGE')) THEN
3547:             IF (NITEMS.GT.1) CALL READF(IMSEPMIN)3535:             IF (NITEMS.GT.1) CALL READF(IMSEPMIN)
3548:             IF (NITEMS.GT.2) CALL READF(IMSEPMAX)3536:             IF (NITEMS.GT.2) CALL READF(IMSEPMAX)
3549:             IF (NITEMS.GT.3) CALL READI(INTIMAGE)3537:             IF (NITEMS.GT.3) CALL READI(INTIMAGE)
3550:             IF (NITEMS.GT.4) CALL READI(MAXINTIMAGE)3538:             IF (NITEMS.GT.4) CALL READI(MAXINTIMAGE)
3551:             IF (NITEMS.GT.5) CALL READI(INTNTRIESMAX)3539:             IF (NITEMS.GT.5) CALL READI(INTNTRIESMAX)
3552:             IF (NITEMS.GT.6) CALL READI(INTIMAGEINCR)3540:             IF (NITEMS.GT.6) CALL READI(INTIMAGEINCR)
3553:             IF (NITEMS.GT.7) CALL READI(INTIMAGECHECK)3541:             IF (NITEMS.GT.7) CALL READI(INTIMAGECHECK)
3554: !3542: !
3555: ! Maximum distance for constrained atoms3543: ! Maximum distance for constrained atoms
3556: !3544: !
3557:          ELSE IF ((WORD.EQ.'INTCONCUT').OR.(WORD.EQ.'QCICONCUT')) THEN3545:          ELSE IF ((WORD.EQ.'INTCONCUT').OR.(WORD.EQ.'QCICONCUT')) THEN
3558:             CALL READF(INTCONCUT)3546:             CALL READF(INTCONCUT)
3559: ! 
3560: ! QCI LOPERMDIST checks for images 
3561: ! 
3562:          ELSE IF ((WORD.EQ.'INTLPERMDIST').OR.(WORD.EQ.'QCILPERMDIST')) THEN 
3563:             QCILPERMDIST=.TRUE. 
3564:             IF (NITEMS.GT.1) CALL READI(QCIPDINT)   ! interval for check 
3565:             IF (NITEMS.GT.2) CALL READI(QCIPERMCUT) ! maximum allowed distance for valid alignment 
3566:             WRITE(*,'(A,I8,A,G20.10)') ' keywords> QCI active permutational alignment check every ',QCIPDINT,' steps, tolerance=',  
3567:      &                                  QCIPERMCUT      
3568: ! 3547: ! 
3569: ! Use constraint potential for initial interpolation in each cycle.3548: ! Use constraint potential for initial interpolation in each cycle.
3570: ! 3549: ! 
3571:          ELSE IF ((WORD.EQ.'INTCONSTRAINT').OR.(WORD.EQ.'QCI')) THEN3550:          ELSE IF ((WORD.EQ.'INTCONSTRAINT').OR.(WORD.EQ.'QCI')) THEN
3572:             INTCONSTRAINTT=.TRUE.3551:             INTCONSTRAINTT=.TRUE.
3573:             IF (NITEMS.GT.1) CALL READF(INTCONSTRAINTTOL)3552:             IF (NITEMS.GT.1) CALL READF(INTCONSTRAINTTOL)
3574:             IF (NITEMS.GT.2) CALL READF(INTCONSTRAINTDEL)3553:             IF (NITEMS.GT.2) CALL READF(INTCONSTRAINTDEL)
3575:             IF (NITEMS.GT.3) CALL READF(INTCONSTRAINTREP)3554:             IF (NITEMS.GT.3) CALL READF(INTCONSTRAINTREP)
3576:             IF (NITEMS.GT.4) CALL READF(INTCONSTRAINREPCUT)3555:             IF (NITEMS.GT.4) CALL READF(INTCONSTRAINREPCUT)
3577:             IF (NITEMS.GT.5) CALL READF(INTCONFRAC)3556:             IF (NITEMS.GT.5) CALL READF(INTCONFRAC)


r33340/lopermdist.f90 2017-09-23 12:30:16.545005743 +0100 r33339/lopermdist.f90 2017-09-23 12:30:18.193027657 +0100
 22: !  Overall alignment is based on the transformation for the best preserved local group. 22: !  Overall alignment is based on the transformation for the best preserved local group.
 23: ! 23: !
 24: !  COORDSA becomes the optimal alignment of the optimal permutation 24: !  COORDSA becomes the optimal alignment of the optimal permutation
 25: !  isomer, but without the permutations. DISTANCE is the residual square distance 25: !  isomer, but without the permutations. DISTANCE is the residual square distance
 26: !  for the best alignment with respect to permutation as well as 26: !  for the best alignment with respect to permutation as well as
 27: !  orientation and centre of mass. 27: !  orientation and centre of mass.
 28: ! 28: !
 29: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 29: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the
 30: !  centre of coordinates of COORDSA will be the same as for COORDSB. 30: !  centre of coordinates of COORDSA will be the same as for COORDSB.
 31: ! 31: !
 32: SUBROUTINE LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST,DOGROUP,NMOVE,NEWPERM) 32: SUBROUTINE LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 33:  33: 
 34: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, & 34: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, &
 35:   &            NFREEZE, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, NTSITES, & 35:   &            NFREEZE, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, NTSITES, &
 36:   &            LOCALPERMCUT, STOCKT, GMAX, CONVR, EDIFFTOL, INTCONSTRAINTT, & 36:   &            LOCALPERMCUT, STOCKT, GMAX, CONVR, EDIFFTOL, INTCONSTRAINTT, &
 37:   &            LOCALPERMNEIGH, LOCALPERMCUT2, ATOMACTIVE 37:   &            LOCALPERMNEIGH, LOCALPERMCUT2
 38: USE MODCHARMM,ONLY : CHRMMT 38: USE MODCHARMM,ONLY : CHRMMT
 39: IMPLICIT NONE 39: IMPLICIT NONE
 40:  40: 
 41: INTEGER, PARAMETER :: MAXIMUMTRIES=10 41: INTEGER, PARAMETER :: MAXIMUMTRIES=10
 42: INTEGER NATOMS, NPERM, PATOMS, NRB, OPNUM,  NORBIT1, NORBIT2, NCHOOSE2, NCHOOSE1, NTRIES, NORBITB1, NORBITB2, NMOVE 42: INTEGER NATOMS, NPERM, PATOMS, NRB, OPNUM,  NORBIT1, NORBIT2, NCHOOSE2, NCHOOSE1, NTRIES, NORBITB1, NORBITB2
 43: INTEGER J3, J4, NDUMMY, LPERM(NATOMS), J1, J2, NOTHER, LPERMBEST(NATOMS), NCHOOSEB1, NCHOOSEB2, & 43: INTEGER J3, J4, NDUMMY, LPERM(NATOMS), J1, J2, NOTHER, LPERMBEST(NATOMS), NCHOOSEB1, NCHOOSEB2, &
 44:         LPERMBESTATOM(NATOMS) 44:         LPERMBESTATOM(NATOMS)
 45: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), & 45: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), &
 46:   &              BESTA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DIST, DSUM 46:   &              BESTA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DIST, DSUM
 47: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS) 47: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS)
 48: DOUBLE PRECISION CMXA, CMXB, CMXC, QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES) 48: DOUBLE PRECISION CMXA, CMXB, CMXC, QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES)
 49: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), RMATBEST(3,3), TMAT(3,3) 49: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), RMATBEST(3,3), TMAT(3,3)
 50: DOUBLE PRECISION PVEC(3), RTEMP1(3,3), RTEMP2(3,3) 50: DOUBLE PRECISION PVEC(3), RTEMP1(3,3), RTEMP2(3,3)
 51: LOGICAL DEBUG, TWOD, RIGID, BULKT, PITEST, AOK, BOK, ADDED, PERMUTABLE(NATOMS) 51: LOGICAL DEBUG, TWOD, RIGID, BULKT, PITEST, AOK, BOK, ADDED, PERMUTABLE(NATOMS)
 52: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY, DUMMYD(3*NATOMS), & 52: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY, DUMMYD(3*NATOMS), &
 53:    &             LDBEST(NPERMGROUP), LDBESTATOM 53:    &             LDBEST(NPERMGROUP), LDBESTATOM
 54: DOUBLE PRECISION SPDUMMYA(3*NATOMS), SPDUMMYB(3*NATOMS), AINIT, BINIT 54: DOUBLE PRECISION SPDUMMYA(3*NATOMS), SPDUMMYB(3*NATOMS), AINIT, BINIT
 55: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS) 55: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS)
 56: DOUBLE PRECISION TIME0, TIME1 56: DOUBLE PRECISION TIME0, TIME1
 57: DOUBLE PRECISION, ALLOCATABLE :: TEMPA(:), TEMPB(:) 57: DOUBLE PRECISION, ALLOCATABLE :: TEMPA(:), TEMPB(:)
 58: CHARACTER(LEN=5) ZSYMSAVE 58: CHARACTER(LEN=5) ZSYMSAVE
 59: COMMON /SYS/ ZSYMSAVE 59: COMMON /SYS/ ZSYMSAVE
 60: DOUBLE PRECISION XA, XB, YA, YB, ZA, ZB, DMEAN(NATOMS), DA, DB 60: DOUBLE PRECISION XA, XB, YA, YB, ZA, ZB, DMEAN(NATOMS), DA, DB
 61: INTEGER TRIED(NATOMS), DLIST(NATOMS), SORTLIST(NATOMS), NDUMMY2, INGROUP(NATOMS), NADDED, DOGROUP 61: INTEGER TRIED(NATOMS), DLIST(NATOMS), SORTLIST(NATOMS), NDUMMY2, INGROUP(NATOMS), NADDED
 62:  62: 
 63: IF (DEBUG.AND.(DOGROUP.EQ.0)) THEN 63: IF (DEBUG) THEN
 64:    IF (CHRMMT) CALL UPDATENBONDS(COORDSA) 64:    IF (CHRMMT) CALL UPDATENBONDS(COORDSA)
 65:    CALL POTENTIAL(COORDSA,AINIT,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 65:    CALL POTENTIAL(COORDSA,AINIT,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
 66:    WRITE(*,'(2(A,G25.15))') ' initial energy for structure A=             ',AINIT,' RMS=',RMS 66:    WRITE(*,'(2(A,G25.15))') ' initial energy for structure A=             ',AINIT,' RMS=',RMS
 67:    IF (RMS-MAX(GMAX,CONVR).GT.1.0D-6) THEN 67:    IF (RMS-MAX(GMAX,CONVR).GT.1.0D-6) THEN
 68:       WRITE(*,'(A)') ' lopermdist> WARNING *** RMS for structure A is outside tolerance' 68:       WRITE(*,'(A)') ' lopermdist> WARNING *** RMS for structure A is outside tolerance'
 69:    ENDIF 69:    ENDIF
 70:    IF (CHRMMT) CALL UPDATENBONDS(COORDSB) 70:    IF (CHRMMT) CALL UPDATENBONDS(COORDSB)
 71:    CALL POTENTIAL(COORDSB,BINIT,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 71:    CALL POTENTIAL(COORDSB,BINIT,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
 72:    PRINT '(2(A,G25.15))',' initial energy for structure B=             ',BINIT,' RMS=',RMS 72:    PRINT '(2(A,G25.15))',' initial energy for structure B=             ',BINIT,' RMS=',RMS
 73:    IF (RMS-MAX(GMAX,CONVR).GT.1.0D-6) THEN 73:    IF (RMS-MAX(GMAX,CONVR).GT.1.0D-6) THEN
 74:       PRINT '(A)',' lopermdist> WARNING *** RMS for structure B is outside tolerance' 74:       PRINT '(A)',' lopermdist> WARNING *** RMS for structure B is outside tolerance'
 75:    ENDIF 75:    ENDIF
 76: ENDIF 76: ENDIF
 77:  77: 
 78: !IF (DOGROUP.GT.0) THEN 
 79: !   WRITE(*,'(A,I6,A)') ' lopermdist> Checking only group ',DOGROUP,' and using only QCI active atoms' 
 80: !ENDIF 
 81:  
 82: DBEST=1.0D100 78: DBEST=1.0D100
 83: PERMUTABLE(1:NATOMS)=.FALSE. 79: PERMUTABLE(1:NATOMS)=.FALSE.
 84: NDUMMY=1 80: NDUMMY=1
 85: DO J1=1,NPERMGROUP 81: DO J1=1,NPERMGROUP
 86:    DO J2=1,NPERMSIZE(J1) 82:    DO J2=1,NPERMSIZE(J1)
 87:       PERMUTABLE(PERMGROUP(NDUMMY+J2-1))=.TRUE. 83:       PERMUTABLE(PERMGROUP(NDUMMY+J2-1))=.TRUE.
 88:       INGROUP(PERMGROUP(NDUMMY+J2-1))=J1 84:       INGROUP(PERMGROUP(NDUMMY+J2-1))=J1
 89:    ENDDO 85:    ENDDO
 90:    NDUMMY=NDUMMY+NPERMSIZE(J1) 86:    NDUMMY=NDUMMY+NPERMSIZE(J1)
 91: ENDDO 87: ENDDO
102: !  The maximum number of pair exchanges associated with a group is two. 98: !  The maximum number of pair exchanges associated with a group is two.
103: !  99: ! 
104: DO J1=1,NATOMS100: DO J1=1,NATOMS
105:    NEWPERM(J1)=J1101:    NEWPERM(J1)=J1
106: ENDDO102: ENDDO
107: DSUM=0.0D0103: DSUM=0.0D0
108: LOCALPERMNEIGH=MIN(LOCALPERMNEIGH,NATOMS)104: LOCALPERMNEIGH=MIN(LOCALPERMNEIGH,NATOMS)
109: 105: 
110: NDUMMY=1106: NDUMMY=1
111: DO J1=1,NPERMGROUP107: DO J1=1,NPERMGROUP
112:    IF (DOGROUP.GT.0) THEN 
113:       IF (J1.GT.DOGROUP) EXIT 
114:       IF (J1.LT.DOGROUP) GOTO 864 ! for QCI test single groups - MUST increment NDUMMY on line 864!! 
115:    ENDIF 
116:    PATOMS=NPERMSIZE(J1)108:    PATOMS=NPERMSIZE(J1)
117: !  PRINT '(A,I6,A,I6)','group ',J1,' size=',NPERMSIZE(J1)109: !  PRINT '(A,I6,A,I6)','group ',J1,' size=',NPERMSIZE(J1)
118: !  PRINT '(A)','members:'110: !  PRINT '(A)','members:'
119: !  DO J2=1,PATOMS111: !  DO J2=1,PATOMS
120: !     PRINT '(20I6)',NEWPERM(PERMGROUP(NDUMMY+J2-1))112: !     PRINT '(20I6)',NEWPERM(PERMGROUP(NDUMMY+J2-1))
121: !  ENDDO113: !  ENDDO
122:    LDBEST(J1)=1.0D100114:    LDBEST(J1)=1.0D100
123:    TRIED(1:NATOMS)=0115:    TRIED(1:NATOMS)=0
124:    DO J2=1,PATOMS116:    DO J2=1,PATOMS
125:       LPERMBEST(J2)=J2117:       LPERMBEST(J2)=J2
165: ! DMEAN, SORTLIST, TRIED, PERMUTABLE, and DLIST entries refer to original157: ! DMEAN, SORTLIST, TRIED, PERMUTABLE, and DLIST entries refer to original
166: ! atom labels. Use NEWPERM to find where they are in coordinate lists.158: ! atom labels. Use NEWPERM to find where they are in coordinate lists.
167: !159: !
168:    outer1: DO J2=1,NATOMS160:    outer1: DO J2=1,NATOMS
169: !161: !
170: ! Don't allow members of the same permutational group 162: ! Don't allow members of the same permutational group 
171: ! to appear as reference neighbours.163: ! to appear as reference neighbours.
172: !164: !
173:       IF (TRIED(J2).EQ.-1) THEN165:       IF (TRIED(J2).EQ.-1) THEN
174:          XDUMMY=1.0D9166:          XDUMMY=1.0D9
175:       ELSE IF ((DOGROUP.GT.0).AND.(.NOT.ATOMACTIVE(J2))) THEN ! only use active atoms for QCI single group 
176:          XDUMMY=1.0D9 
177:       ELSE167:       ELSE
178:          DA=(XA-DUMMYA(3*(NEWPERM(J2)-1)+1))**2 &168:          DA=(XA-DUMMYA(3*(NEWPERM(J2)-1)+1))**2 &
179:   &        +(YA-DUMMYA(3*(NEWPERM(J2)-1)+2))**2 &169:   &        +(YA-DUMMYA(3*(NEWPERM(J2)-1)+2))**2 &
180:   &        +(ZA-DUMMYA(3*(NEWPERM(J2)-1)+3))**2170:   &        +(ZA-DUMMYA(3*(NEWPERM(J2)-1)+3))**2
181: !        DB=(XB-DUMMYB(3*(NEWPERM(J2)-1)+1))**2 &171: !        DB=(XB-DUMMYB(3*(NEWPERM(J2)-1)+1))**2 &
182: ! &        +(YB-DUMMYB(3*(NEWPERM(J2)-1)+2))**2 &172: ! &        +(YB-DUMMYB(3*(NEWPERM(J2)-1)+2))**2 &
183: ! &        +(ZB-DUMMYB(3*(NEWPERM(J2)-1)+3))**2173: ! &        +(ZB-DUMMYB(3*(NEWPERM(J2)-1)+3))**2
184:          DB=(XB-DUMMYB(3*(J2-1)+1))**2 &174:          DB=(XB-DUMMYB(3*(J2-1)+1))**2 &
185:   &        +(YB-DUMMYB(3*(J2-1)+2))**2 &175:   &        +(YB-DUMMYB(3*(J2-1)+2))**2 &
186:   &        +(ZB-DUMMYB(3*(J2-1)+3))**2176:   &        +(ZB-DUMMYB(3*(J2-1)+3))**2
312: !  DO J2=1,PATOMS302: !  DO J2=1,PATOMS
313: !     PRINT '(2I6)',PERMGROUP(NDUMMY+J2-1),NEWPERM(PERMGROUP(NDUMMY+J2-1))303: !     PRINT '(2I6)',PERMGROUP(NDUMMY+J2-1),NEWPERM(PERMGROUP(NDUMMY+J2-1))
314: !  ENDDO304: !  ENDDO
315:    305:    
316: !  PRINT '(I6,A)',NOTHER,' other atoms:'306: !  PRINT '(I6,A)',NOTHER,' other atoms:'
317: !  PRINT '(A)','original and new other atom labels:'307: !  PRINT '(A)','original and new other atom labels:'
318: !  DO J2=1,NOTHER308: !  DO J2=1,NOTHER
319: !     PRINT '(2I6)',DLIST(J2),NEWPERM(DLIST(J2))309: !     PRINT '(2I6)',DLIST(J2),NEWPERM(DLIST(J2))
320: !  ENDDO310: !  ENDDO
321: 311: 
322: !  PRINT '(A,3I6,3G20.10)','J1,PATOMS,NOTHER,root LDBEST(J1),root LDISTANCE=',J1,PATOMS,NOTHER,SQRT(LDBEST(J1)),SQRT(LDISTANCE)312: !  PRINT '(A,3I6,3G20.10)','J1,PATOMS,NOTHER,LDBEST(J1),LDISTANCE=',J1,PATOMS,NOTHER,LDBEST(J1),LDISTANCE
323: !  PRINT '(A,20I6)','LPERM after MINPERM: ',LPERM(1:PATOMS+NOTHER)313: !  PRINT '(A,20I6)','LPERM after MINPERM: ',LPERM(1:PATOMS+NOTHER)
324: 314: 
325:    LDISTANCE=LDISTANCE315:    LDISTANCE=LDISTANCE
326:    DO J2=1,PATOMS316:    DO J2=1,PATOMS
327:       IF (LPERM(J2).GT.PATOMS) THEN317:       IF (LPERM(J2).GT.PATOMS) THEN
328:          LDISTANCE=1.0D300318:          LDISTANCE=1.0D300
329: !        IF (DEBUG) PRINT '(A,I6,A,I6,A)',' lopermdist> For group ',J1,' with ',NOTHER,' neighbours - neighbours mix in' 319: !        IF (DEBUG) PRINT '(A,I6,A,I6,A)',' lopermdist> For group ',J1,' with ',NOTHER,' neighbours - neighbours mix in' 
330: !        IF (DEBUG) PRINT '(A,I6,A,I6)',' lopermdist> atom ',J2,' lperm value is ',LPERM(J2)320: !        IF (DEBUG) PRINT '(A,I6,A,I6)',' lopermdist> atom ',J2,' lperm value is ',LPERM(J2)
331:          EXIT321:          EXIT
332:       ENDIF322:       ENDIF
347:    ENDDO337:    ENDDO
348: !338: !
349: ! Save the best permutation and local distance for this subset of atoms.339: ! Save the best permutation and local distance for this subset of atoms.
350: ! NEWPERM and coordinates are only reset after all the cycles over orbits and NEWMINDIST.340: ! NEWPERM and coordinates are only reset after all the cycles over orbits and NEWMINDIST.
351: ! Hence we need to track a cumulative permutation and save the best current values.341: ! Hence we need to track a cumulative permutation and save the best current values.
352: !342: !
353:    IF (LDISTANCE.LT.LDBESTATOM) THEN343:    IF (LDISTANCE.LT.LDBESTATOM) THEN
354:       LDBESTATOM=LDISTANCE344:       LDBESTATOM=LDISTANCE
355:       LPERMBESTATOM(1:PATOMS)=LPERM(1:PATOMS)345:       LPERMBESTATOM(1:PATOMS)=LPERM(1:PATOMS)
356:    ENDIF346:    ENDIF
357: !  PRINT '(A,2G20.10)','root LDISTANCE,root LDBESTATOM=',SQRT(LDISTANCE),SQRT(LDBESTATOM)347: !  PRINT '(A,2G20.10)','LDISTANCE,LDBESTATOM=',LDISTANCE,LDBESTATOM
358: 348: 
359: !  PRINT '(A,4I6,2G20.10)','NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,root LDISTANCE,root LDBEST=', &349: !  PRINT '(A,4I6,2G20.10)','NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,LDISTANCE,LDBEST=', &
360: ! &                         NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,SQRT(LDISTANCE),SQRT(LDBEST(J1))350: ! &                         NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,LDISTANCE,LDBEST(J1)
361: 351: 
362:    IF (NCHOOSE2.LT.NORBIT2) GOTO 30352:    IF (NCHOOSE2.LT.NORBIT2) GOTO 30
363:    IF (NCHOOSE1.LT.NORBIT1) GOTO 65353:    IF (NCHOOSE1.LT.NORBIT1) GOTO 65
364:    IF (NCHOOSEB2.LT.NORBITB2) GOTO 31354:    IF (NCHOOSEB2.LT.NORBITB2) GOTO 31
365:    IF (NCHOOSEB1.LT.NORBITB1) GOTO 66355:    IF (NCHOOSEB1.LT.NORBITB1) GOTO 66
366: 356: 
367: !  PRINT '(A,2G20.10)','root LDBESTATOM,LOCALPERMCUT=',SQRT(LDBESTATOM),LOCALPERMCUT 357: !  PRINT '(A,2G20.10)','LDBESTATOM,LOCALPERMCUT=',LDBESTATOM,LOCALPERMCUT
368:    IF (SQRT(LDBESTATOM).GT.LOCALPERMCUT) THEN358:    IF (SQRT(LDBESTATOM).GT.LOCALPERMCUT) THEN
369:       IF (DEBUG) THEN359: !     IF (DEBUG) THEN
370: !        PRINT '(A,G15.5,A,I6)',' lopermdist> Best distance ',SQRT(LDBESTATOM), &360: !        PRINT '(A,G15.5,A,I6)',' lopermdist> Best distance ',SQRT(LDBESTATOM), &
371: ! &                                     ' is too large for atom ',DLIST(NOTHER)361: ! &                                     ' is too large for atom ',DLIST(NOTHER)
372:       ENDIF362: !     ENDIF
373:       TRIED(DLIST(NOTHER))=-1363:       TRIED(DLIST(NOTHER))=-1
374:       IF (NADDED.GT.1) THEN364:       IF (NADDED.GT.1) THEN
375:          IF (DEBUG) THEN365: !        IF (DEBUG) THEN
376: !           PRINT '(A)',' lopermdist> and partner atoms:'366: !           PRINT '(A)',' lopermdist> and partner atoms:'
377: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)367: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)
378: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:'368: !        ENDIF
379: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1) 
380:          ENDIF 
381:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1369:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1
382:       ENDIF370:       ENDIF
383:       GOTO 71371:       GOTO 71
384:    ELSE372:    ELSE
385:       IF (DEBUG) PRINT '(A,G20.10,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), &373: !     IF (DEBUG) PRINT '(A,G20.10,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), &
386:   &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours' 374: ! &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours' 
387:       TRIED(DLIST(NOTHER))=1375:       TRIED(DLIST(NOTHER))=1
388:       IF (NADDED.GT.1) THEN376:       IF (NADDED.GT.1) THEN
389: !        IF (DEBUG) THEN377: !        IF (DEBUG) THEN
390: !           PRINT '(A)',' lopermdist> and partner atoms:'378: !           PRINT '(A)',' lopermdist> and partner atoms:'
391: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)379: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)
392: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:' 
393: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1) 
394: !        ENDIF380: !        ENDIF
395:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1381:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1
396:       ENDIF382:       ENDIF
397: !     IF (DEBUG) THEN383: !     IF (DEBUG) THEN
398: !        PRINT '(A)',' all other atoms:'384: !        PRINT '(A)',' all other atoms:'
399: !        PRINT '(20I5)',DLIST(1:NOTHER)385: !        PRINT '(20I5)',DLIST(1:NOTHER)
400: !        IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(DLIST(1:NOTHER)) 
401: !     ENDIF386: !     ENDIF
402:       LDBEST(J1)=LDBESTATOM387:       LDBEST(J1)=LDBESTATOM
403:       LPERMBEST(1:PATOMS)=LPERMBESTATOM(1:PATOMS)388:       LPERMBEST(1:PATOMS)=LPERMBESTATOM(1:PATOMS)
404: !     PRINT '(A,2G20.10)','Updating permutation: sqrt(LDBEST)=',SQRT(LDBEST(J1))389: !     PRINT '(A,2G20.10)','Updating permutation: sqrt(LDBEST)=',SQRT(LDBEST(J1))
405: !     PRINT '(A,10I6)','LPERMBEST: ',LPERMBEST(1:PATOMS)390: !     PRINT '(A,10I6)','LPERMBEST: ',LPERMBEST(1:PATOMS)
406:    ENDIF391:    ENDIF
407: !392: !
408: ! Add the next eligible atom and try alignment again.393: ! Add the next eligible atom and try alignment again.
409: ! Stop if we already have LOCALPERMNEIGH neighbours.394: ! Stop if we already have LOCALPERMNEIGH neighbours.
410: !395: !
444:          DO J3=1,NSETS(J1)429:          DO J3=1,NSETS(J1)
445:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3)430:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3)
446:          ENDDO431:          ENDDO
447:       ENDDO432:       ENDDO
448:    ENDIF433:    ENDIF
449: !434: !
450: ! Save current optimal permutation in NEWPERM435: ! Save current optimal permutation in NEWPERM
451: !436: !
452:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS)437:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS)
453:    DSUM=DSUM+SQRT(LDBEST(J1))438:    DSUM=DSUM+SQRT(LDBEST(J1))
454:    PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' best distance=',SQRT(LDBEST(J1)),' total=',DSUM439: !  PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' best distance=',SQRT(LDBEST(J1)),' total=',DSUM
455: !  PRINT '(A)','best permutation is now'440: !  PRINT '(A)','best permutation is now'
456: !  PRINT '(20I6)',NEWPERM(1:NATOMS)441: !  PRINT '(20I6)',NEWPERM(1:NATOMS)
457: 442: 
458: !443: !
459: ! Update NDUMMY, the cumulative offset for PERMGROUP444: ! Update NDUMMY, the cumulative offset for PERMGROUP
460: !445: !
461: 864   NDUMMY=NDUMMY+NPERMSIZE(J1)446:    NDUMMY=NDUMMY+NPERMSIZE(J1)
462: ENDDO  !  end of loop over groups of permutable atoms447: ENDDO  !  end of loop over groups of permutable atoms
463:  
464: IF (DOGROUP.GT.0) THEN 
465:    DISTANCE=SQRT(LDBEST(DOGROUP)) 
466:    NMOVE=0 
467:    DO J2=1,NATOMS 
468:       IF (NEWPERM(J2).NE.J2) THEN 
469: !        WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2 
470:          NMOVE=NMOVE+1 
471:       ENDIF 
472:    ENDDO 
473: !  PRINT '(A,I6,A,I6)',' lopermdist> Total permutations for optimal alignment of group ',DOGROUP,' is ',NMOVE 
474:    RETURN 
475: ENDIF 
476: NMOVE=0 
477: DO J2=1,NATOMS 
478:    IF (NEWPERM(J2).NE.J2) THEN 
479: !     WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2 
480:       NMOVE=NMOVE+1 
481:    ENDIF 
482: ENDDO 
483: ! PRINT '(A,I6)',' lopermdist> Total permutations for optimal alignment=',NMOVE 
484: !448: !
485: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA449: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA
486: ! to the current best alignment. 450: ! to the current best alignment. 
487: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current451: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current
488: ! NEWPERM.452: ! NEWPERM.
489: !453: !
490: !454: !
491: ! Putting the ALLPERM(J1)=J1 into the second loop causes pgf90 to miscompile!!455: ! Putting the ALLPERM(J1)=J1 into the second loop causes pgf90 to miscompile!!
492: !456: !
493: DO J1=1,NATOMS457: DO J1=1,NATOMS
556: CALL NEWMINDIST(DUMMYB,XBEST,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)520: CALL NEWMINDIST(DUMMYB,XBEST,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
557: IF (DEBUG) PRINT '(A,G20.10)',' lopermdist> after overall alignment distance=',DISTANCE521: IF (DEBUG) PRINT '(A,G20.10)',' lopermdist> after overall alignment distance=',DISTANCE
558: RMATBEST(1:3,1:3)=RMAT(1:3,1:3)522: RMATBEST(1:3,1:3)=RMAT(1:3,1:3)
559: 523: 
560: COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!524: COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!
561: 525: 
562: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!526: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
563: ! IF (DEBUG) PRINT '(A)',' lopermdist> Overall permutation for COORDSA (second argument):'527: ! IF (DEBUG) PRINT '(A)',' lopermdist> Overall permutation for COORDSA (second argument):'
564: ! IF (DEBUG) PRINT '(20I6)',BESTPERM(1:NATOMS)528: ! IF (DEBUG) PRINT '(20I6)',BESTPERM(1:NATOMS)
565: 529: 
566: IF (DEBUG.AND.(DOGROUP.EQ.0)) THEN530: IF (DEBUG) THEN
567:    IF (CHRMMT) CALL UPDATENBONDS(COORDSA)531:    IF (CHRMMT) CALL UPDATENBONDS(COORDSA)
568:    CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)532:    CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
569:    PRINT '(2(A,G25.15))',' lopermdist> final   energy for structure A=             ',ENERGY,' RMS=',RMS533:    PRINT '(2(A,G25.15))',' lopermdist> final   energy for structure A=             ',ENERGY,' RMS=',RMS
570:    IF (ABS(ENERGY-AINIT).GT.EDIFFTOL) THEN534:    IF (ABS(ENERGY-AINIT).GT.EDIFFTOL) THEN
571:       PRINT '(A)',' minpermdist> WARNING *** energy change for structure A is outside tolerance'535:       PRINT '(A)',' minpermdist> WARNING *** energy change for structure A is outside tolerance'
572:    ENDIF536:    ENDIF
573:    IF (CHRMMT) CALL UPDATENBONDS(COORDSB)537:    IF (CHRMMT) CALL UPDATENBONDS(COORDSB)
574:    CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)538:    CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
575:    PRINT '(2(A,G25.15))',' lopermdist> final   energy for structure B=             ',ENERGY,' RMS=',RMS539:    PRINT '(2(A,G25.15))',' lopermdist> final   energy for structure B=             ',ENERGY,' RMS=',RMS
576:    IF (ABS(ENERGY-BINIT).GT.EDIFFTOL) THEN540:    IF (ABS(ENERGY-BINIT).GT.EDIFFTOL) THEN


r33340/minpermdist.f90 2017-09-23 12:30:16.765008670 +0100 r33339/minpermdist.f90 2017-09-23 12:30:18.421030688 +0100
 62: USE MODCHARMM,ONLY : CHRMMT 62: USE MODCHARMM,ONLY : CHRMMT
 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH
 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET
 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE
 66: USE GENRIGID 66: USE GENRIGID
 67: USE AMBER12_INTERFACE_MOD 67: USE AMBER12_INTERFACE_MOD
 68: USE CHIRALITY 68: USE CHIRALITY
 69: IMPLICIT NONE 69: IMPLICIT NONE
 70:  70: 
 71: INTEGER :: MAXIMUMTRIES=10 71: INTEGER :: MAXIMUMTRIES=10
 72: INTEGER NATOMS, NPERM, PATOMS, NTRIES, NRB, OPNUM, BESTINVERT, I, LOPERM(NATOMS) 72: INTEGER NATOMS, NPERM, PATOMS, NTRIES, NRB, OPNUM, BESTINVERT, I
 73: INTEGER J3, INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1, NROTDONE, NORBITB1, NORBITB2, & 73: INTEGER J3, INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1, NROTDONE, NORBITB1, NORBITB2, &
 74:   &     NCHOOSEB1, NCHOOSEB2 74:   &     NCHOOSEB1, NCHOOSEB2
 75: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), & 75: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), &
 76:   &              DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DX, DY, DZ 76:   &              DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DX, DY, DZ
 77: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), RMS, DBEST, XBEST(3*NATOMS) 77: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), RMS, DBEST, XBEST(3*NATOMS)
 78: DOUBLE PRECISION MAXE1, MAXE2, DISTANCE1, SAVECUT, DIST, AINIT, BINIT 78: DOUBLE PRECISION MAXE1, MAXE2, DISTANCE1, SAVECUT, DIST, AINIT, BINIT
 79: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES), CMX, CMY, CMZ 79: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES), CMX, CMY, CMZ
 80: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), ROTINVBBEST(3,3), ROTABEST(3,3), RMATBEST(3,3), TMAT(3,3) 80: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), ROTINVBBEST(3,3), ROTABEST(3,3), RMATBEST(3,3), TMAT(3,3)
 81: DOUBLE PRECISION CMAX, CMAY, CMAZ, CMBX, CMBY, CMBZ, RMATCUMUL(3,3) 81: DOUBLE PRECISION CMAX, CMAY, CMAZ, CMBX, CMBY, CMBZ, RMATCUMUL(3,3)
 82: DOUBLE PRECISION REFXZ(3,3) 82: DOUBLE PRECISION REFXZ(3,3)
 83: LOGICAL DEBUG, TWOD, RIGID, BULKT, PITEST, TNMATCH, BMTEST, LDB 83: LOGICAL DEBUG, TWOD, RIGID, BULKT, PITEST, TNMATCH, BMTEST, LDB
 84: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY 84: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY
 85: DOUBLE PRECISION BMDIST, BMCOORDS(3*NATOMS), BMCOORDSSV(3*NATOMS) 85: DOUBLE PRECISION BMDIST, BMCOORDS(3*NATOMS), BMCOORDSSV(3*NATOMS)
 86: DOUBLE PRECISION TEMPCOORDSA(DEGFREEDOMS), TEMPCOORDSB(DEGFREEDOMS) ! sn402 86: DOUBLE PRECISION TEMPCOORDSA(DEGFREEDOMS), TEMPCOORDSB(DEGFREEDOMS) ! sn402
 87: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS), NMOVE 87: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS)
 88: CHARACTER(LEN=5) ZSYMSAVE 88: CHARACTER(LEN=5) ZSYMSAVE
 89: COMMON /SYS/ ZSYMSAVE 89: COMMON /SYS/ ZSYMSAVE
 90:  90: 
 91: LDB=.FALSE. 91: LDB=.FALSE.
 92: ! hk286 92: ! hk286
 93: IF (GTHOMSONT) THEN 93: IF (GTHOMSONT) THEN
 94:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST) 94:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 95:    RETURN 95:    RETURN
 96: ELSEIF (VARIABLES) THEN 96: ELSEIF (VARIABLES) THEN
 97:    DISTANCE=0.0D0 97:    DISTANCE=0.0D0
263:             IF (.NOT.TWOD) COORDSA(3*(J1-1)+3)=COORDSA(3*(J1-1)+3)+DZ263:             IF (.NOT.TWOD) COORDSA(3*(J1-1)+3)=COORDSA(3*(J1-1)+3)+DZ
264:          ENDDO264:          ENDDO
265:       ENDIF265:       ENDIF
266:       RMATBEST = RMAT266:       RMATBEST = RMAT
267:    ELSE267:    ELSE
268:       CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)268:       CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)
269:       RMATBEST = RMAT      ! hk286 - this line keeps getting removed by someone else. If this is not suitable for you, please let me know.269:       RMATBEST = RMAT      ! hk286 - this line keeps getting removed by someone else. If this is not suitable for you, please let me know.
270:    ENDIF270:    ENDIF
271:    RETURN271:    RETURN
272: ELSEIF (LPERMDIST) THEN272: ELSEIF (LPERMDIST) THEN
273:    CALL LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST,0,NMOVE,LOPERM)273:    CALL LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
274:    IF (DEBUG) THEN274:    IF (DEBUG) THEN
275:       IF (CHRMMT) CALL UPDATENBONDS(COORDSA)275:       IF (CHRMMT) CALL UPDATENBONDS(COORDSA)
276:       IF (RIGIDINIT ) THEN276:       IF (RIGIDINIT ) THEN
277:          CALL GENRIGID_POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)277:          CALL GENRIGID_POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
278:       ELSE278:       ELSE
279:          CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)279:          CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
280:       ENDIF280:       ENDIF
281:       PRINT '(2(A,G25.15))',' minpermdist> final   energy for structure A=             ',ENERGY,' RMS=',RMS281:       PRINT '(2(A,G25.15))',' minpermdist> final   energy for structure A=             ',ENERGY,' RMS=',RMS
282:       IF (ABS(ENERGY-AINIT).GT.2*EDIFFTOL) THEN282:       IF (ABS(ENERGY-AINIT).GT.2*EDIFFTOL) THEN
283:          PRINT '(A)',' minpermdist> ERROR *** energy change for structure A is outside tolerance - QCI/DNEB endpoint alignment?'283:          PRINT '(A)',' minpermdist> ERROR *** energy change for structure A is outside tolerance - QCI/DNEB endpoint alignment?'


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0