hdiff output

r31259/amber_top_reader.f90 2016-10-10 11:30:09.085197270 +0100 r31258/amber_top_reader.f90 2016-10-10 11:30:12.273239631 +0100
 84:                ELSE 84:                ELSE
 85:                   GOTO 98 85:                   GOTO 98
 86:                ENDIF 86:                ENDIF
 87:                J4=J4+1 87:                J4=J4+1
 88:             ENDDO 88:             ENDDO
 89:          ENDDO 89:          ENDDO
 90:       ENDIF 90:       ENDIF
 91:  91: 
 92:    ENDDO reading 92:    ENDDO reading
 93: 99 CLOSE(MYUNIT2) 93: 99 CLOSE(MYUNIT2)
 94: !  DO J6=1,NBOND 94:    DO J6=1,NBOND
 95: !     WRITE(MYUNIT,'(A,I8,A,I8)') 'readtopology> Bond between',BONDS(J6,1),' and',BONDS(J6,2) 95:       WRITE(MYUNIT,'(A,I8,A,I8)') 'readtopology> Bond between',BONDS(J6,1),' and',BONDS(J6,2)
 96: !  ENDDO 96:    ENDDO
 97:  97: 
 98: END SUBROUTINE 98: END SUBROUTINE
 99:    99:   
100: SUBROUTINE READ_LINE(LINE,NWORDS,WORDSOUT)100: SUBROUTINE READ_LINE(LINE,NWORDS,WORDSOUT)
101:       CHARACTER(*), INTENT(IN) :: LINE101:       CHARACTER(*), INTENT(IN) :: LINE
102:       INTEGER, INTENT(IN) :: NWORDS102:       INTEGER, INTENT(IN) :: NWORDS
103:       CHARACTER(*), DIMENSION(NWORDS), INTENT(OUT) :: WORDSOUT103:       CHARACTER(*), DIMENSION(NWORDS), INTENT(OUT) :: WORDSOUT
104:       INTEGER:: J1,START_IND,END_IND,J2104:       INTEGER:: J1,START_IND,END_IND,J2
105:       CHARACTER(25) :: WORD105:       CHARACTER(25) :: WORD
106:       START_IND=0106:       START_IND=0


r31259/commons.f90 2016-10-10 11:30:09.349200787 +0100 r31258/commons.f90 2016-10-10 11:30:12.529243037 +0100
 32:      &        MYEUNIT, MYMUNIT, MYBUNIT, MYRUNIT, MYPUNIT, NFREEZETYPEA, & 32:      &        MYEUNIT, MYMUNIT, MYBUNIT, MYRUNIT, MYPUNIT, NFREEZETYPEA, &
 33:      &        TBPSTEPS, TBPCI, TBPBASIN, NTSITES, NRBGROUP, NZERO, PTMCDS_FRQ, PTMCDUMPENERFRQ, MONITORINT, NBLOCKS, & 33:      &        TBPSTEPS, TBPCI, TBPBASIN, NTSITES, NRBGROUP, NZERO, PTMCDS_FRQ, PTMCDUMPENERFRQ, MONITORINT, NBLOCKS, &
 34:      &        BINARY_EXAB_FRQ, NRESMIN, USERES, EXEQ, NONEDAPBC, STRUC, CHEMSHIFTITER, GRIDSIZE, MFETRUNS, BESTINVERT, GCNATOMS, & 34:      &        BINARY_EXAB_FRQ, NRESMIN, USERES, EXEQ, NONEDAPBC, STRUC, CHEMSHIFTITER, GRIDSIZE, MFETRUNS, BESTINVERT, GCNATOMS, &
 35:      &        GCINT, GCRELAX, MTARGETS, & 35:      &        GCINT, GCRELAX, MTARGETS, &
 36:      &        INTCONSEP, INTREPSEP, NCONSTRAINTON, CPREPSEP, CPCONSEP, NCONGEOM, & 36:      &        INTCONSEP, INTREPSEP, NCONSTRAINTON, CPREPSEP, CPCONSEP, NCONGEOM, &
 37:      &        NCPREPULSIVE, NCPCONSTRAINT, MAXCONUSE, INTCONSTEPS, INTRELSTEPS, INTSTEPS1, INTLJSTEPS, & 37:      &        NCPREPULSIVE, NCPCONSTRAINT, MAXCONUSE, INTCONSTEPS, INTRELSTEPS, INTSTEPS1, INTLJSTEPS, &
 38:      &        NTRAPPOW, MAXINTIMAGE, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, & 38:      &        NTRAPPOW, MAXINTIMAGE, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, &
 39:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, INTIMAGE, NREPULSIVE, & 39:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, INTIMAGE, NREPULSIVE, &
 40:      &        NNREPULSIVE, NCONSTRAINT, INTMUPDATE, DUMPINTEOSFREQ, DUMPINTXYZFREQ, & 40:      &        NNREPULSIVE, NCONSTRAINT, INTMUPDATE, DUMPINTEOSFREQ, DUMPINTXYZFREQ, &
 41:      &        LOCALPERMNEIGH, LOCALPERMMAXSEP, MAXNACTIVE, QCIPERMCHECKINT, & 41:      &        LOCALPERMNEIGH, LOCALPERMMAXSEP, MAXNACTIVE, QCIPERMCHECKINT, &
 42:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, DJWRBID, NHEXAMERS, QCIADDREP, QCIBONDS, QCISECOND, MQUNIT 42:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, DJWRBID, NHEXAMERS, QCIADDREP, QCIBONDS, QCISECOND
 43:   43:  
 44:       DOUBLE PRECISION RHO, GAMMA, SIG, SCEPS, SCC, TOLB, T12FAC, XMOVERENORM, RESIZE, QTSALLIS, & 44:       DOUBLE PRECISION RHO, GAMMA, SIG, SCEPS, SCC, TOLB, T12FAC, XMOVERENORM, RESIZE, QTSALLIS, &
 45:      &                 CQMAX, RADIUS, BQMAX,  MAXBFGS, DECAYPARAM, SYMTOL1, SYMTOL2, SYMTOL3, SYMTOL4, SYMTOL5, PGSYMTOLS(3),& 45:      &                 CQMAX, RADIUS, BQMAX,  MAXBFGS, DECAYPARAM, SYMTOL1, SYMTOL2, SYMTOL3, SYMTOL4, SYMTOL5, PGSYMTOLS(3),&
 46:      &                 ECONV, TOLD, TOLE, SYMREM(120,3,3), GMAX, CUTOFF, PCUT, EXPFAC, EXPD, CENTX, CENTY, CENTZ, & 46:      &                 ECONV, TOLD, TOLE, SYMREM(120,3,3), GMAX, CUTOFF, PCUT, EXPFAC, EXPD, CENTX, CENTY, CENTZ, &
 47:      &                 BOXLX, BOXLY, BOXLZ, BOX3D(3), PCUTOFF, SUPSTEP, SQUEEZER, SQUEEZED, COOPCUT, STOCKMU, STOCKLAMBDA, & 47:      &                 BOXLX, BOXLY, BOXLZ, BOX3D(3), PCUTOFF, SUPSTEP, SQUEEZER, SQUEEZED, COOPCUT, STOCKMU, STOCKLAMBDA, &
 48:      &                 TFAC(3), RMS, TEMPS, SACCRAT, CEIG, PNEWJUMP, EAMP, DISTFAC, ODDCHARGE, COULQ, COULSWAP, & 48:      &                 TFAC(3), RMS, TEMPS, SACCRAT, CEIG, PNEWJUMP, EAMP, DISTFAC, ODDCHARGE, COULQ, COULSWAP, &
 49:      &                 COULTEMP, APP, AMM, APM, XQP, XQM, ALPHAP, ALPHAM, ZSTAR, K_COMP, DGUESS, GUIDECUT, EFAC,&  49:      &                 COULTEMP, APP, AMM, APM, XQP, XQM, ALPHAP, ALPHAM, ZSTAR, K_COMP, DGUESS, GUIDECUT, EFAC,& 
 50:      &                 TRENORM, HISTMIN, HISTMAX, HISTFAC, EPSSPHERE, FINALCUTOFF, SHELLPROB, RINGROTSCALE, & 50:      &                 TRENORM, HISTMIN, HISTMAX, HISTFAC, EPSSPHERE, FINALCUTOFF, SHELLPROB, RINGROTSCALE, &
 51:      &                 HISTFACMUL, HPERCENT, AVOIDDIST, MAXERISE, MAXEFALL, TSTART, MATDIFF, STICKYSIG, SDTOL, & 51:      &                 HISTFACMUL, HPERCENT, AVOIDDIST, MAXERISE, MAXEFALL, TSTART, MATDIFF, STICKYSIG, SDTOL, &
 52:      &                 MinimalTemperature, MaximalTemperature, SwapProb, hdistconstraint, COREFRAC, TSTAR, & 52:      &                 MinimalTemperature, MaximalTemperature, SwapProb, hdistconstraint, COREFRAC, TSTAR, &
110:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&110:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
111:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &111:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &
112:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &112:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
113:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &113:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
114:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &114:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
115:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &115:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
116:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 116:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
117:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &117:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
118:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &118:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
119:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &119:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
120:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT, &120:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT, RIGIDMDT
121:      &        DUMPMQT 
122: !121: !
123:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 122:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
124:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)123:       DOUBLE PRECISION, ALLOCATABLE :: ATMASS(:)
125:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:) 124:       DOUBLE PRECISION, ALLOCATABLE :: SPECMASS(:) 
126: 125: 
127: ! csw34> FREEZEGROUP variables126: ! csw34> FREEZEGROUP variables
128: !127: !
129:       INTEGER :: GROUPCENTRE128:       INTEGER :: GROUPCENTRE
130:       DOUBLE PRECISION :: GROUPRADIUS129:       DOUBLE PRECISION :: GROUPRADIUS
131:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE130:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE
583:       LOGICAL :: DBYUKAWAT582:       LOGICAL :: DBYUKAWAT
584:       DOUBLE PRECISION :: LAMBDAYAA, LAMBDAYBB, LAMBDAYAB, YEPSFAC 583:       DOUBLE PRECISION :: LAMBDAYAA, LAMBDAYBB, LAMBDAYAB, YEPSFAC 
585: ! hk286 - RESTRAIN LENGTH584: ! hk286 - RESTRAIN LENGTH
586:       LOGICAL :: RESTRAINLT585:       LOGICAL :: RESTRAINLT
587:       DOUBLE PRECISION :: RESTRAINLK586:       DOUBLE PRECISION :: RESTRAINLK
588:       DOUBLE PRECISION, ALLOCATABLE :: RESTRAINLDIST(:)587:       DOUBLE PRECISION, ALLOCATABLE :: RESTRAINLDIST(:)
589:       INTEGER :: RESTRAINLNOPAIR588:       INTEGER :: RESTRAINLNOPAIR
590:       INTEGER, ALLOCATABLE :: RESTRAINLPAIRS(:,:)589:       INTEGER, ALLOCATABLE :: RESTRAINLPAIRS(:,:)
591: ! hk286590: ! hk286
592:       LOGICAL :: GTHOMSONT591:       LOGICAL :: GTHOMSONT
593:       DOUBLE PRECISION :: GThomsonZ, GThomsonC, GThomsonC2, VUNDULOID, GTHOMSONBIND592:       DOUBLE PRECISION :: GThomsonZ, GThomsonC, GThomsonC2, VUNDULOID
594:       DOUBLE PRECISION :: GTrefU, GTrefZ, GTmu, GTk, GTm, GTn, GTa, GTc593:       DOUBLE PRECISION :: GTrefU, GTrefZ, GTmu, GTk, GTm, GTn, GTa, GTc
595:       INTEGER :: GTHOMMET, GTHOMPOT, GTHOMSONBINN594:       INTEGER :: GTHOMMET, GTHOMPOT
596:       DOUBLE PRECISION :: GThomsonSigma, GThomsonRho595:       DOUBLE PRECISION :: GThomsonSigma, GThomsonRho
597: 596: 
598: ! dc550597: ! dc550
599:       LOGICAL :: SELECTMOVET598:       LOGICAL :: SELECTMOVET
600:       DOUBLE PRECISION, ALLOCATABLE :: SELTRANSSTEP(:), SELROTSCALE(:), SELMOVPROB(:)599:       DOUBLE PRECISION, ALLOCATABLE :: SELTRANSSTEP(:), SELROTSCALE(:), SELMOVPROB(:)
601:       INTEGER, ALLOCATABLE :: SELMOVFREQ(:), SELBEGIN(:), SELEND(:), SELSIZE(:)600:       INTEGER, ALLOCATABLE :: SELMOVFREQ(:), SELBEGIN(:), SELEND(:), SELSIZE(:)
602:       INTEGER :: SELMOVNO601:       INTEGER :: SELMOVNO
603:       602:       
604: !QUIP variables603: !QUIP variables
605:       CHARACTER(LEN=3) :: QUIPATOMTYPE604:       CHARACTER(LEN=3) :: QUIPATOMTYPE


r31259/congrad.f90 2016-10-10 11:30:09.609204235 +0100 r31258/congrad.f90 2016-10-10 11:30:12.785246433 +0100
 69:          G2(2)=(R2AY-R2BY)/D2 69:          G2(2)=(R2AY-R2BY)/D2
 70:          G2(3)=(R2AZ-R2BZ)/D2 70:          G2(3)=(R2AZ-R2BZ)/D2
 71:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3) 71:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
 72:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2) 72:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
 73:          EEE(J1)=EEE(J1)  +DUMMY 73:          EEE(J1)=EEE(J1)  +DUMMY
 74:          ECON=ECON        +DUMMY 74:          ECON=ECON        +DUMMY
 75:          CONE(J1)=CONE(J1)+DUMMY 75:          CONE(J1)=CONE(J1)+DUMMY
 76:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 76:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
 77:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 77:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
 78:       ENDIF 78:       ENDIF
 79: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1) 
 80:    ENDDO 79:    ENDDO
 81: ENDDO 80: ENDDO
 82:  81: 
 83: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes 82: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
 84: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes 83: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
 85:  84: 
 86: ! INTCONST=INTCONSTRAINREPCUT**13 85: ! INTCONST=INTCONSTRAINREPCUT**13
 87:  86: 
 88: DO J2=1,NNREPULSIVE 87: DO J2=1,NNREPULSIVE
 89: !  INTCONST=NREPCUT(J2)**13 88: !  INTCONST=NREPCUT(J2)**13
 92: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images 91: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images
 93:       IF (FREEZENODEST) THEN 92:       IF (FREEZENODEST) THEN
 94:          IF (J1.EQ.2) THEN 93:          IF (J1.EQ.2) THEN
 95:             IF (IMGFREEZE(1)) CYCLE 94:             IF (IMGFREEZE(1)) CYCLE
 96:          ELSE IF (J1.EQ.INTIMAGE+2) THEN 95:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
 97:             IF (IMGFREEZE(INTIMAGE)) CYCLE 96:             IF (IMGFREEZE(INTIMAGE)) CYCLE
 98:          ELSE 97:          ELSE
 99:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE 98:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
100:          ENDIF 99:          ENDIF
101:       ENDIF100:       ENDIF
102:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN101:    IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
103:          WRITE(MYUNIT, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)102:       WRITE(MYUNIT, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
104:          STOP103:       STOP
105:       ENDIF104:    ENDIF
106: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)105: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad2> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)
107:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)106:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
108:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)107:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
109:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)108:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
110:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)109:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
111:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)110:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)
112:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1111:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1
113: !        D12=D2**12112: !        D12=D2**12
114:          D12=D2**2113:          D12=D2**2
115: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)114: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
116:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)115:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
119:          EREP=EREP+DUMMY118:          EREP=EREP+DUMMY
120: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)119: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
121:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)120:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
122:          G2(1)=(R2AX-R2BX)/D2121:          G2(1)=(R2AX-R2BX)/D2
123:          G2(2)=(R2AY-R2BY)/D2122:          G2(2)=(R2AY-R2BY)/D2
124:          G2(3)=(R2AZ-R2BZ)/D2123:          G2(3)=(R2AZ-R2BZ)/D2
125:          REPGRAD(1:3)=DUMMY*G2(1:3)124:          REPGRAD(1:3)=DUMMY*G2(1:3)
126:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)125:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
127:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)126:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
128:       ENDIF127:       ENDIF
129: !     WRITE(MYUNIT,'(A,2I6,4G20.10)') 'J1,J2,D2,NREPCUT,EEE,REPE=',J1,J2,D2,NREPCUT(J2),EEE(J1),REPE(J1) 
130: !128: !
131: ! For internal minima we are counting edges. 129: ! For internal minima we are counting edges. 
132: ! Edge J1 is between images J1-1 and J1, starting from J1=2.130: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
133: ! Energy contributions are shared evenly, except for131: ! Energy contributions are shared evenly, except for
134: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which132: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which
135: ! is assigned to image INTIMAGE+1. Gradients are set to zero for133: ! is assigned to image INTIMAGE+1. Gradients are set to zero for
136: ! the end images.134: ! the end images.
137: !135: !
138:       IF (J1.EQ.1) CYCLE136:       IF (J1.EQ.1) CYCLE
139:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)137:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
175:             NREPINT(J1-1)=NREPINT(J1-1)+1173:             NREPINT(J1-1)=NREPINT(J1-1)+1
176:          ELSE IF (J1.EQ.INTIMAGE+2) THEN174:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
177:             EEE(J1-1)=EEE(J1-1)+DUMMY175:             EEE(J1-1)=EEE(J1-1)+DUMMY
178:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY176:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
179:             NREPINT(J1-1)=NREPINT(J1-1)+1177:             NREPINT(J1-1)=NREPINT(J1-1)+1
180:          ENDIF178:          ENDIF
181:          EREP=EREP+DUMMY179:          EREP=EREP+DUMMY
182: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)180: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
183:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)181:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
184:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)182:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)
185: !        WRITE(MYUNIT,'(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &183: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
186: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)184: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
187: !185: !
188: ! Gradient contributions for image J1-1186: ! Gradient contributions for image J1-1
189: !187: !
190:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)188:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
191:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)189:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
192:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)190:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)
193: !        WRITE(MYUNIT,'(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &191: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
194: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)192: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
195: !193: !
196: ! Gradient contributions for image J1194: ! Gradient contributions for image J1
197: !195: !
198:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)196:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
199:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)197:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
200:       ENDIF198:       ENDIF
201:    ENDDO199:    ENDDO
202: ENDDO200: ENDDO
203: !201: !
696:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)694:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
697:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)695:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
698:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)696:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
699:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)697:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
700:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)698:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
701: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN699: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
702:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-50) THEN700:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-50) THEN
703: !        WRITE(MYUNIT, '(A,I6,A,2I6)') 'A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)701: !        WRITE(MYUNIT, '(A,I6,A,2I6)') 'A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
704: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ702: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
705: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ703: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
706: !        WRITE(MYUNIT,'(A,7I10)') 'congrad2> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)704: !        WRITE(MYUNIT,'(A,7I10)') 'congrad> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)
707: !        WRITE(MYUNIT,'(A,7I10)') 'frames ',J1-1,J1705: !        WRITE(MYUNIT,'(A,7I10)') 'frames ',J1-1,J1
708:          D1=1.0D100; D2=1.0D100; NOINT=.TRUE.  ! to skip the next blocks706:          D1=1.0D100; D2=1.0D100; NOINT=.TRUE.  ! to skip the next blocks
709:       ELSE707:       ELSE
710:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &708:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
711:   &                    D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))709:   &                    D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
712:       ENDIF710:       ENDIF
713: !711: !
714: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.712: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
715: !713: !
716: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1714: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
1033:          G2(2)=(R2AY-R2BY)/D21031:          G2(2)=(R2AY-R2BY)/D2
1034:          G2(3)=(R2AZ-R2BZ)/D21032:          G2(3)=(R2AZ-R2BZ)/D2
1035:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)1033:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
1036:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)1034:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
1037:          EEE(J1)=EEE(J1)  +DUMMY1035:          EEE(J1)=EEE(J1)  +DUMMY
1038:          ECON=ECON        +DUMMY1036:          ECON=ECON        +DUMMY
1039:          CONE(J1)=CONE(J1)+DUMMY1037:          CONE(J1)=CONE(J1)+DUMMY
1040:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)1038:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
1041:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)1039:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
1042:       ENDIF1040:       ENDIF
1043: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1) 
1044:       IF (J2.GT.QCIBONDS) CYCLE1041:       IF (J2.GT.QCIBONDS) CYCLE
1045:       DO J3=J2+1,QCIBONDS1042:       DO J3=J2+1,QCIBONDS
1046:          IF (.NOT.CONACTIVE(J3)) CYCLE1043:          IF (.NOT.CONACTIVE(J3)) CYCLE
1047:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom1044:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom
1048:          IF (CONI(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom1045:          IF (CONI(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom
1049:          IF (CONJ(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom1046:          IF (CONJ(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom
1050:          IF (CONJ(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom1047:          IF (CONJ(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom
1051:          NI2=(3*NATOMS)*(J1-1)+3*(CONI(J3)-1)1048:          NI2=(3*NATOMS)*(J1-1)+3*(CONI(J3)-1)
1052:          NJ2=(3*NATOMS)*(J1-1)+3*(CONJ(J3)-1)1049:          NJ2=(3*NATOMS)*(J1-1)+3*(CONJ(J3)-1)
1053:          DO J4=1,QCIADDREP1050:          DO J4=1,QCIADDREP
1097:       IF (FREEZENODEST) THEN1094:       IF (FREEZENODEST) THEN
1098:          IF (J1.EQ.2) THEN1095:          IF (J1.EQ.2) THEN
1099:             IF (IMGFREEZE(1)) CYCLE1096:             IF (IMGFREEZE(1)) CYCLE
1100:          ELSE IF (J1.EQ.INTIMAGE+2) THEN1097:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
1101:             IF (IMGFREEZE(INTIMAGE)) CYCLE1098:             IF (IMGFREEZE(INTIMAGE)) CYCLE
1102:          ELSE1099:          ELSE
1103:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE1100:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
1104:          ENDIF1101:          ENDIF
1105:       ENDIF1102:       ENDIF
1106:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN1103:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
1107:          WRITE(MYUNIT, '(A,I6,A,2I6)') ' congrad3> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)1104:          WRITE(MYUNIT, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
1108:          STOP1105:          STOP
1109:       ENDIF1106:       ENDIF
1110: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad3> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)1107: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad2> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)
1111:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)1108:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
1112:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)1109:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
1113:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)1110:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
1114:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)1111:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
1115:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)1112:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)
1116:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J11113:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1
1117: !        D12=D2**121114: !        D12=D2**12
1118:          D12=D2**21115:          D12=D2**2
1119: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)1116: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
1120:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)1117:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
1123:          EREP=EREP+DUMMY1120:          EREP=EREP+DUMMY
1124: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)1121: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
1125:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)1122:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
1126:          G2(1)=(R2AX-R2BX)/D21123:          G2(1)=(R2AX-R2BX)/D2
1127:          G2(2)=(R2AY-R2BY)/D21124:          G2(2)=(R2AY-R2BY)/D2
1128:          G2(3)=(R2AZ-R2BZ)/D21125:          G2(3)=(R2AZ-R2BZ)/D2
1129:          REPGRAD(1:3)=DUMMY*G2(1:3)1126:          REPGRAD(1:3)=DUMMY*G2(1:3)
1130:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)1127:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
1131:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)1128:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
1132:       ENDIF1129:       ENDIF
1133: !     WRITE(MYUNIT,'(A,2I6,4G20.10)') 'J1,J2,D2,NREPCUT,EEE,REPE=',J1,J2,D2,NREPCUT(J2),EEE(J1),REPE(J1) 
1134: !1130: !
1135: ! For internal minima we are counting edges. 1131: ! For internal minima we are counting edges. 
1136: ! Edge J1 is between images J1-1 and J1, starting from J1=2.1132: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
1137: ! Energy contributions are shared evenly, except for1133: ! Energy contributions are shared evenly, except for
1138: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which1134: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which
1139: ! is assigned to image INTIMAGE+1. Gradients are set to zero for1135: ! is assigned to image INTIMAGE+1. Gradients are set to zero for
1140: ! the end images.1136: ! the end images.
1141: !1137: !
1142:       IF (J1.EQ.1) CYCLE1138:       IF (J1.EQ.1) CYCLE
1143:       IF (QCINOREPINT) CYCLE1139:       IF (QCINOREPINT) CYCLE
1144:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)1140:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
1145:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)1141:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
1146:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)1142:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
1147:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)1143:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
1148: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN1144: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
1149:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-10) THEN1145:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-10) THEN
1150: !        WRITE(MYUNIT, '(A,I6,A,2I6)') 'B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)1146: !        WRITE(MYUNIT, '(A,I6,A,2I6)') 'B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
1151: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ1147: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
1152: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ1148: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
1153: !        WRITE(MYUNIT,'(A,7I10)') 'congrad3> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)1149: !        WRITE(MYUNIT,'(A,7I10)') 'congrad> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)
1154: !        WRITE(MYUNIT,'(A,7I10)') 'frames ',J1-1,J11150: !        WRITE(MYUNIT,'(A,7I10)') 'frames ',J1-1,J1
1155:          NOINT=.TRUE.1151:          NOINT=.TRUE.
1156:       ELSE1152:       ELSE
1157:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &1153:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
1158:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))1154:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
1159:          IF (.NOT.NOINT) THEN1155:          IF (.NOT.NOINT) THEN
1160: !           WRITE(MYUNIT,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad3> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), &1156: !           WRITE(MYUNIT,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), &
1161: ! &                        ' distance,cutoff=',DINT,NREPCUT(J2)1157: ! &                        ' distance,cutoff=',DINT,NREPCUT(J2)
1162:             NINTMIN=NINTMIN+11158:             NINTMIN=NINTMIN+1
1163:          ENDIF1159:          ENDIF
1164:       ENDIF1160:       ENDIF
1165:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN1161:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
1166:          NINTMIN2=NINTMIN2+11162:          NINTMIN2=NINTMIN2+1
1167: !        D12=DSQI**61163: !        D12=DSQI**6
1168:          D12=DSQI1164:          D12=DSQI
1169: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)1165: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
1170:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)1166:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
1290: ENDIF1286: ENDIF
1291: !1287: !
1292: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,1288: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,
1293: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)1289: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)
1294: !1290: !
1295: ETOTAL=0.0D01291: ETOTAL=0.0D0
1296: MAXINT=-1.0D1001292: MAXINT=-1.0D100
1297: MININT=1.0D1001293: MININT=1.0D100
1298: DO J1=2,INTIMAGE+11294: DO J1=2,INTIMAGE+1
1299:    ETOTAL=ETOTAL+EEE(J1)1295:    ETOTAL=ETOTAL+EEE(J1)
1300: !  WRITE(MYUNIT, '(A,I6,A,3G20.10)') ' congrad3> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)1296: !  WRITE(MYUNIT, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)
1301:    IF (REPEINT(J1).LT.MININT) THEN1297:    IF (REPEINT(J1).LT.MININT) THEN
1302:       MININT=REPEINT(J1)1298:       MININT=REPEINT(J1)
1303:       NMININT=J11299:       NMININT=J1
1304:    ENDIF1300:    ENDIF
1305:    IF (REPE(J1).GT.MAXINT) THEN1301:    IF (REPE(J1).GT.MAXINT) THEN
1306:       MAXINT=REPE(J1)1302:       MAXINT=REPE(J1)
1307:       NMAXINT=J11303:       NMAXINT=J1
1308:    ENDIF1304:    ENDIF
1309: ENDDO1305: ENDDO
1310: 1306: 


r31259/finalio.f90 2016-10-10 11:30:09.865207635 +0100 r31258/finalio.f90 2016-10-10 11:30:13.053249945 +0100
515:         WRITE(MYUNIT2,65) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))515:         WRITE(MYUNIT2,65) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))
516: 65      FORMAT('AR ',3F20.10)516: 65      FORMAT('AR ',3F20.10)
517:      ELSE IF (TOSI.OR.WELCH) THEN517:      ELSE IF (TOSI.OR.WELCH) THEN
518:         DO J2=1,NATOMS518:         DO J2=1,NATOMS
519:            IF (ZSYM(J2).EQ.'PL') WRITE(MYUNIT2,'(A,3F20.10)') 'Na  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)519:            IF (ZSYM(J2).EQ.'PL') WRITE(MYUNIT2,'(A,3F20.10)') 'Na  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
520:            IF (ZSYM(J2).EQ.'MI') WRITE(MYUNIT2,'(A,3F20.10)') 'Cl  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)520:            IF (ZSYM(J2).EQ.'MI') WRITE(MYUNIT2,'(A,3F20.10)') 'Cl  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
521:         ENDDO521:         ENDDO
522:         ! hk286 - generalised Thomson problem, convert to Cartesians before printing522:         ! hk286 - generalised Thomson problem, convert to Cartesians before printing
523:      ELSE IF (GTHOMSONT) THEN523:      ELSE IF (GTHOMSONT) THEN
524:         CALL GTHOMSONANGTOC(DCOORDS(1:3*NATOMS), QMINP(J1, 1:3*NATOMS), NATOMS) 524:         CALL GTHOMSONANGTOC(DCOORDS(1:3*NATOMS), QMINP(J1, 1:3*NATOMS), NATOMS) 
525:         IF (GTHOMPOT.EQ.7) THEN525:         DO J2=1,NATOMS
526:            DO J2=1,NATOMS526:            WRITE(MYUNIT2,'(A,3F20.10)') 'C  ',(DCOORDS(3*(J2-1)+J3),J3=1,3)
527:               IF (J2.LE.GTHOMSONBINN) THEN527:         ENDDO
528:                  WRITE(MYUNIT2,'(A,3F20.10)') 'LA  ',(DCOORDS(3*(J2-1)+J3),J3=1,3) 
529:               ELSE 
530:                  WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(DCOORDS(3*(J2-1)+J3),J3=1,3) 
531:               ENDIF 
532:            ENDDO 
533:         ELSE 
534:            DO J2=1,NATOMS 
535:               WRITE(MYUNIT2,'(A,3F20.10)') 'C  ',(DCOORDS(3*(J2-1)+J3),J3=1,3) 
536:            ENDDO 
537:         ENDIF 
538:      ELSE IF ((BLJCLUSTER.OR.BLJCLUSTER_NOCUT).AND..NOT.QALCST) THEN528:      ELSE IF ((BLJCLUSTER.OR.BLJCLUSTER_NOCUT).AND..NOT.QALCST) THEN
539:         DO J2=1,NATOMS529:         DO J2=1,NATOMS
540:            IF (J2.LE.NTYPEA) THEN530:            IF (J2.LE.NTYPEA) THEN
541:               WRITE(MYUNIT2,'(A,3F20.10)') 'LA  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)531:               WRITE(MYUNIT2,'(A,3F20.10)') 'LA  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
542:            ELSE532:            ELSE
543:               WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)533:               WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
544:            ENDIF534:            ENDIF
545:         ENDDO535:         ENDDO
546:      ELSE IF (GLJT) THEN536:      ELSE IF (GLJT) THEN
547:         J3=1537:         J3=1


r31259/Gthomson.f90 2016-10-10 11:30:08.569190412 +0100 r31258/Gthomson.f90 2016-10-10 11:30:11.757232776 +0100
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18:  18: 
 19:       SUBROUTINE GTHOMSON(X,V,ETHOMSON,GTEST) 19:       SUBROUTINE GTHOMSON(X,V,ETHOMSON,GTEST)
 20:  20: 
 21:       USE COMMONS, ONLY: GTHOMMET, GTHOMPOT, GThomsonRho, GThomsonSigma, NATOMS, GTHOMSONBINN, GTHOMSONBIND 21:       USE COMMONS, ONLY: GTHOMMET, GTHOMPOT, GThomsonRho, GThomsonSigma, NATOMS
 22:       IMPLICIT NONE 22:       IMPLICIT NONE
 23:  23: 
 24:       DOUBLE PRECISION, INTENT(INOUT)  :: X(*) 24:       DOUBLE PRECISION, INTENT(INOUT)  :: X(*)
 25:       DOUBLE PRECISION, INTENT(OUT) :: V(*), ETHOMSON 25:       DOUBLE PRECISION, INTENT(OUT) :: V(*), ETHOMSON
 26:       LOGICAL, INTENT(IN) :: GTEST 26:       LOGICAL, INTENT(IN) :: GTEST
 27:       INTEGER J1, J2 27:       INTEGER J1, J2
 28:       DOUBLE PRECISION DIST 28:       DOUBLE PRECISION DIST
 29:       DOUBLE PRECISION :: TMPCOORDS(3*NATOMS), DR(3), Gradient11(3), Gradient12(3), Gradient21(3), Gradient22(3) 29:       DOUBLE PRECISION :: TMPCOORDS(3*NATOMS), DR(3), Gradient11(3), Gradient12(3), Gradient21(3), Gradient22(3)
 30:  30: 
 31: ! jwrm2> Reset the range of the polar coordinates 31: ! jwrm2> Reset the range of the polar coordinates
 63:             ELSEIF (GTHOMPOT .EQ. 2) THEN  63:             ELSEIF (GTHOMPOT .EQ. 2) THEN 
 64:                ETHOMSON = ETHOMSON + 1/DIST**3 ! 1/R^3 64:                ETHOMSON = ETHOMSON + 1/DIST**3 ! 1/R^3
 65:             ELSEIF (GTHOMPOT .EQ. 3) THEN 65:             ELSEIF (GTHOMPOT .EQ. 3) THEN
 66:                ETHOMSON = ETHOMSON + 1/DIST*EXP(-1.0D0*DIST/GThomsonSigma)                !Yukawa 66:                ETHOMSON = ETHOMSON + 1/DIST*EXP(-1.0D0*DIST/GThomsonSigma)                !Yukawa
 67:             ELSEIF (GTHOMPOT .EQ. 4) THEN 67:             ELSEIF (GTHOMPOT .EQ. 4) THEN
 68:                ETHOMSON = ETHOMSON + (GThomsonSigma/DIST)**12 - (GThomsonSigma/DIST)**6  ! LJ 68:                ETHOMSON = ETHOMSON + (GThomsonSigma/DIST)**12 - (GThomsonSigma/DIST)**6  ! LJ
 69:             ELSEIF (GTHOMPOT .EQ. 5) THEN 69:             ELSEIF (GTHOMPOT .EQ. 5) THEN
 70:                ETHOMSON = ETHOMSON + (GThomsonSigma/DIST)**12                            !Repulsive LJ  70:                ETHOMSON = ETHOMSON + (GThomsonSigma/DIST)**12                            !Repulsive LJ 
 71:             ELSEIF (GTHOMPOT .EQ. 6) THEN 71:             ELSEIF (GTHOMPOT .EQ. 6) THEN
 72:                ETHOMSON = ETHOMSON + (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST)))**2 -1.0D0 ! Morse 72:                ETHOMSON = ETHOMSON + (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST)))**2 -1.0D0 ! Morse
 73:             ELSEIF (GTHOMPOT .EQ. 7) THEN 
 74:                IF ((J1.LE.GTHOMSONBINN).AND.(J2.GT.GTHOMSONBINN)) THEN 
 75:                   ETHOMSON = ETHOMSON + 1.0D0/SQRT(DIST**2+GTHOMSONBIND**2/NATOMS) ! softened 
 76:                ELSE 
 77:                   ETHOMSON = ETHOMSON + 1.0D0/DIST ! Coulomb 
 78:                ENDIF 
 79:             ENDIF 73:             ENDIF
 80:  74: 
  75: 
 81:             IF (GTEST) THEN 76:             IF (GTEST) THEN
 82:  77: 
 83:                IF (GTHOMPOT .EQ. 1) THEN 78:                IF (GTHOMPOT .EQ. 1) THEN
 84:                   V(3*J1-2) = V(3*J1-2) - 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient11) 79:                   V(3*J1-2) = V(3*J1-2) - 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient11)
 85:                   V(3*J1-1) = V(3*J1-1) - 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient12)                80:                   V(3*J1-1) = V(3*J1-1) - 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient12)               
 86:                ELSEIF (GTHOMPOT .EQ. 2) THEN  81:                ELSEIF (GTHOMPOT .EQ. 2) THEN 
 87:                   V(3*J1-2) = V(3*J1-2) - 3.0D0/DIST**5 * DOT_PRODUCT(DR,Gradient11) 82:                   V(3*J1-2) = V(3*J1-2) - 3.0D0/DIST**5 * DOT_PRODUCT(DR,Gradient11)
 88:                   V(3*J1-1) = V(3*J1-1) - 3.0D0/DIST**5 * DOT_PRODUCT(DR,Gradient12)                83:                   V(3*J1-1) = V(3*J1-1) - 3.0D0/DIST**5 * DOT_PRODUCT(DR,Gradient12)               
 89:                ELSEIF (GTHOMPOT .EQ. 3) THEN 84:                ELSEIF (GTHOMPOT .EQ. 3) THEN
 90:                   V(3*J1-2) = V(3*J1-2) - (1.0D0/DIST**3 + 1.0D0/GThomsonSigma/DIST**2) * EXP(-1.0D0*DIST/GThomsonSigma) & 85:                   V(3*J1-2) = V(3*J1-2) - (1.0D0/DIST**3 + 1.0D0/GThomsonSigma/DIST**2) * EXP(-1.0D0*DIST/GThomsonSigma) &
 97:                   V(3*J1-1) = V(3*J1-1) - ((12.0D0*GThomsonSigma**12/DIST**14) - (6.0D0*GThomsonSigma**6/DIST**8)) &  92:                   V(3*J1-1) = V(3*J1-1) - ((12.0D0*GThomsonSigma**12/DIST**14) - (6.0D0*GThomsonSigma**6/DIST**8)) & 
 98:                        * DOT_PRODUCT(DR,Gradient12) 93:                        * DOT_PRODUCT(DR,Gradient12)
 99:                ELSEIF (GTHOMPOT .EQ. 5) THEN 94:                ELSEIF (GTHOMPOT .EQ. 5) THEN
100:                   V(3*J1-2) = V(3*J1-2) - (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient11) ! Repulsive LJ 95:                   V(3*J1-2) = V(3*J1-2) - (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient11) ! Repulsive LJ
101:                   V(3*J1-1) = V(3*J1-1) - (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient12) ! Replusive LJ 96:                   V(3*J1-1) = V(3*J1-1) - (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient12) ! Replusive LJ
102:                ELSEIF (GTHOMPOT .EQ. 6) THEN 97:                ELSEIF (GTHOMPOT .EQ. 6) THEN
103:                   V(3*J1-2) = V(3*J1-2) + 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) * &  98:                   V(3*J1-2) = V(3*J1-2) + 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) * & 
104:                        EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient11) ! Morse 99:                        EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient11) ! Morse
105:                   V(3*J1-1) = V(3*J1-1) + 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) * &100:                   V(3*J1-1) = V(3*J1-1) + 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) * &
106:                        EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient12) ! Morse101:                        EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient12) ! Morse
107:                ELSEIF (GTHOMPOT .EQ. 7) THEN 
108:                   IF ((J1.LE.GTHOMSONBINN).AND.(J2.GT.GTHOMSONBINN)) THEN 
109:                      V(3*J1-2) = V(3*J1-2) - 1.0D0/(DIST**2+GTHOMSONBIND**2/NATOMS)**1.5D0 * DOT_PRODUCT(DR,Gradient11) 
110:                      V(3*J1-1) = V(3*J1-1) - 1.0D0/(DIST**2+GTHOMSONBIND**2/NATOMS)**1.5D0 * DOT_PRODUCT(DR,Gradient12)                
111:                   ELSE 
112:                      V(3*J1-2) = V(3*J1-2) - 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient11) 
113:                      V(3*J1-1) = V(3*J1-1) - 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient12)                
114:                   ENDIF 
115:                ENDIF102:                ENDIF
116:  103:  
117:                IF ( GTHOMMET .EQ. 1) THEN104:                IF ( GTHOMMET .EQ. 1) THEN
118:                   CALL GRADMETRICCYLINDER (X(3*J2-2:3*J2-1), Gradient21, Gradient22)105:                   CALL GRADMETRICCYLINDER (X(3*J2-2:3*J2-1), Gradient21, Gradient22)
119:                ELSEIF ( GTHOMMET .EQ. 2) THEN106:                ELSEIF ( GTHOMMET .EQ. 2) THEN
120:                   CALL GRADMETRICCATENOID (X(3*J2-2:3*J2-1), Gradient21, Gradient22)107:                   CALL GRADMETRICCATENOID (X(3*J2-2:3*J2-1), Gradient21, Gradient22)
121:                ELSEIF ( (GTHOMMET .EQ. 3) .OR. (GTHOMMET .EQ. 4) ) THEN108:                ELSEIF ( (GTHOMMET .EQ. 3) .OR. (GTHOMMET .EQ. 4) ) THEN
122:                   CALL GRADMETRICUNDULOID (X(3*J2-2:3*J2-1), Gradient21, Gradient22)109:                   CALL GRADMETRICUNDULOID (X(3*J2-2:3*J2-1), Gradient21, Gradient22)
123:                ELSEIF ( GTHOMMET .EQ. 5 ) THEN110:                ELSEIF ( GTHOMMET .EQ. 5 ) THEN
124:                   CALL GRADMETRICSPHERE (X(3*J2-2:3*J2-1), Gradient21, Gradient22)111:                   CALL GRADMETRICSPHERE (X(3*J2-2:3*J2-1), Gradient21, Gradient22)
143:                   V(3*J2-1) = V(3*J2-1) + ((12.0D0*GThomsonSigma**12/DIST**14) - (6.0D0*GThomsonSigma**6/DIST**8))&130:                   V(3*J2-1) = V(3*J2-1) + ((12.0D0*GThomsonSigma**12/DIST**14) - (6.0D0*GThomsonSigma**6/DIST**8))&
144:                        * DOT_PRODUCT(DR,Gradient22)131:                        * DOT_PRODUCT(DR,Gradient22)
145:                ELSEIF (GTHOMPOT .EQ. 5) THEN132:                ELSEIF (GTHOMPOT .EQ. 5) THEN
146:                   V(3*J2-2) = V(3*J2-2) + (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient21) ! Repulsive LJ133:                   V(3*J2-2) = V(3*J2-2) + (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient21) ! Repulsive LJ
147:                   V(3*J2-1) = V(3*J2-1) + (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient22) ! Repulsive LJ134:                   V(3*J2-1) = V(3*J2-1) + (12.0D0*GThomsonSigma**12/DIST**14)* DOT_PRODUCT(DR,Gradient22) ! Repulsive LJ
148:                ELSEIF (GTHOMPOT .EQ. 6) THEN135:                ELSEIF (GTHOMPOT .EQ. 6) THEN
149:                   V(3*J2-2) = V(3*J2-2) - 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) &136:                   V(3*J2-2) = V(3*J2-2) - 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) &
150:                        * EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient21)137:                        * EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient21)
151:                   V(3*J2-1) = V(3*J2-1) - 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) &138:                   V(3*J2-1) = V(3*J2-1) - 2.0D0 * (1.0D0 - EXP(GThomsonRho*(GThomsonSigma-DIST))) &
152:                        * EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient22) 139:                        * EXP(GThomsonRho*(GThomsonSigma-DIST)) * GThomsonRho / DIST * DOT_PRODUCT(DR,Gradient22) 
153:                ELSEIF (GTHOMPOT .EQ. 7) THEN 
154:                   IF ((J1.LE.GTHOMSONBINN).AND.(J2.GT.GTHOMSONBINN)) THEN 
155:                      V(3*J2-2) = V(3*J2-2) + 1.0D0/(DIST**2+GTHOMSONBIND**2/NATOMS)**1.5D0 * DOT_PRODUCT(DR,Gradient21) 
156:                      V(3*J2-1) = V(3*J2-1) + 1.0D0/(DIST**2+GTHOMSONBIND**2/NATOMS)**1.5D0 * DOT_PRODUCT(DR,Gradient22)                
157:                   ELSE 
158:                      V(3*J2-2) = V(3*J2-2) + 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient21)                
159:                      V(3*J2-1) = V(3*J2-1) + 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient22) 
160:                   ENDIF 
161:                ENDIF140:                ENDIF
162:                !PRINT *, 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient2), (1/DIST2 - 1/DIST)/0.00001141:                !PRINT *, 1.0D0/DIST**3 * DOT_PRODUCT(DR,Gradient2), (1/DIST2 - 1/DIST)/0.00001
163: 142: 
164:             ENDIF143:             ENDIF
165: 144: 
166:          ENDDO145:          ENDDO
167:       ENDDO146:       ENDDO
168: 147: 
169:       RETURN148:       RETURN
170:     END SUBROUTINE GTHOMSON149:     END SUBROUTINE GTHOMSON
963:          942:          
964:       PRINT *, "zmax ", GTHOMSONZ943:       PRINT *, "zmax ", GTHOMSONZ
965: 944: 
966:     END SUBROUTINE CONVERTUNDULOIDPARAMETERS945:     END SUBROUTINE CONVERTUNDULOIDPARAMETERS
967: 946: 
968: 947: 
969: !----------------------------------------------------------------948: !----------------------------------------------------------------
970: 949: 
971: SUBROUTINE GTHOMSONHESSIAN(X,HESS)950: SUBROUTINE GTHOMSONHESSIAN(X,HESS)
972: 951: 
973:   USE commons, ONLY: GTHOMMET, GTHOMPOT, GThomsonSigma, NATOMS, GTHOMSONBINN, GTHOMSONBIND952:   USE commons, ONLY: GTHOMMET, GTHOMPOT, GThomsonSigma, NATOMS
974:   IMPLICIT NONE953:   IMPLICIT NONE
975: 954: 
976:   INTEGER J1, J2, J3, J4955:   INTEGER J1, J2, J3, J4
977:   DOUBLE PRECISION, INTENT(IN)  :: X(*)956:   DOUBLE PRECISION, INTENT(IN)  :: X(*)
978:   DOUBLE PRECISION, INTENT(OUT) :: HESS(3*NATOMS, 3*NATOMS)957:   DOUBLE PRECISION, INTENT(OUT) :: HESS(3*NATOMS, 3*NATOMS)
979:   DOUBLE PRECISION DIST, D2VDR2, DVDR958:   DOUBLE PRECISION DIST, D2VDR2, DVDR
980:   DOUBLE PRECISION :: TMPCOORDS(3*NATOMS), DR(3), Gradient11(3), Gradient12(3), Gradient21(3), Gradient22(3) 959:   DOUBLE PRECISION :: TMPCOORDS(3*NATOMS), DR(3), Gradient11(3), Gradient12(3), Gradient21(3), Gradient22(3) 
981:   DOUBLE PRECISION :: HESS111(3), HESS112(3), HESS122(3)960:   DOUBLE PRECISION :: HESS111(3), HESS112(3), HESS122(3)
982:   DOUBLE PRECISION :: HESS211(3), HESS212(3), HESS222(3)961:   DOUBLE PRECISION :: HESS211(3), HESS212(3), HESS222(3)
983: 962: 
1019:            DVDR   = -( 1/DIST**3 + 1/GThomsonSigma/DIST**2 ) * EXP(-1.0D0*DIST/GThomsonSigma)998:            DVDR   = -( 1/DIST**3 + 1/GThomsonSigma/DIST**2 ) * EXP(-1.0D0*DIST/GThomsonSigma)
1020:         ELSEIF (GTHOMPOT .EQ. 4) THEN999:         ELSEIF (GTHOMPOT .EQ. 4) THEN
1021:            PRINT *, "POTENTIAL NOT IMPLEMENTED IN HESSIAN"1000:            PRINT *, "POTENTIAL NOT IMPLEMENTED IN HESSIAN"
1022:            STOP1001:            STOP
1023:         ELSEIF (GTHOMPOT .EQ. 5) THEN1002:         ELSEIF (GTHOMPOT .EQ. 5) THEN
1024:            PRINT *, "POTENTIAL NOT IMPLEMENTED IN HESSIAN"1003:            PRINT *, "POTENTIAL NOT IMPLEMENTED IN HESSIAN"
1025:            STOP1004:            STOP
1026:         ELSEIF (GTHOMPOT .EQ. 6) THEN1005:         ELSEIF (GTHOMPOT .EQ. 6) THEN
1027:            PRINT *, "POTENTIAL NOT IMPLEMENTED IN HESSIAN"1006:            PRINT *, "POTENTIAL NOT IMPLEMENTED IN HESSIAN"
1028:            STOP1007:            STOP
1029:         ELSEIF (GTHOMPOT .EQ. 7) THEN 
1030:            IF ((J1.LE.GTHOMSONBINN).AND.(J2.GT.GTHOMSONBINN)) THEN 
1031:               D2VDR2 =   3/(DIST**2+GTHOMSONBIND**2/NATOMS)**2.5D0 
1032:               DVDR   = - 1/(DIST**2+GTHOMSONBIND**2/NATOMS)**1.5D0 
1033:            ELSE 
1034:               D2VDR2 =   3/DIST**5 
1035:               DVDR   = - 1/DIST**3 
1036:            ENDIF 
1037:         ENDIF1008:         ENDIF
1038:         1009:         
1039:         HESS(3*J1-2,3*J1-2) = HESS(3*J1-2,3*J1-2) +  D2VDR2 * (DOT_PRODUCT(DR,Gradient11))**2 &1010:         HESS(3*J1-2,3*J1-2) = HESS(3*J1-2,3*J1-2) +  D2VDR2 * (DOT_PRODUCT(DR,Gradient11))**2 &
1040:              + DVDR * (DOT_PRODUCT(Gradient11,Gradient11)+ DOT_PRODUCT(DR,HESS111)) 1011:              + DVDR * (DOT_PRODUCT(Gradient11,Gradient11)+ DOT_PRODUCT(DR,HESS111)) 
1041:         HESS(3*J1-1,3*J1-1) = HESS(3*J1-1,3*J1-1) +  D2VDR2 * (DOT_PRODUCT(DR,Gradient12))**2 &1012:         HESS(3*J1-1,3*J1-1) = HESS(3*J1-1,3*J1-1) +  D2VDR2 * (DOT_PRODUCT(DR,Gradient12))**2 &
1042:              + DVDR * (DOT_PRODUCT(Gradient12,Gradient12)+ DOT_PRODUCT(DR,HESS122)) 1013:              + DVDR * (DOT_PRODUCT(Gradient12,Gradient12)+ DOT_PRODUCT(DR,HESS122)) 
1043:         HESS(3*J1-2,3*J1-1) = HESS(3*J1-2,3*J1-1) +  D2VDR2 * DOT_PRODUCT(DR,Gradient11) * DOT_PRODUCT(DR,Gradient12) &1014:         HESS(3*J1-2,3*J1-1) = HESS(3*J1-2,3*J1-1) +  D2VDR2 * DOT_PRODUCT(DR,Gradient11) * DOT_PRODUCT(DR,Gradient12) &
1044:              + DVDR * (DOT_PRODUCT(Gradient11,Gradient12) + DOT_PRODUCT(DR,HESS112)) 1015:              + DVDR * (DOT_PRODUCT(Gradient11,Gradient12) + DOT_PRODUCT(DR,HESS112)) 
1045: 1016: 
1046:         IF ( GTHOMMET .EQ. 1) THEN1017:         IF ( GTHOMMET .EQ. 1) THEN


r31259/intlbfgs.f90 2016-10-10 11:30:10.129211142 +0100 r31258/intlbfgs.f90 2016-10-10 11:30:13.313253451 +0100
102:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &102:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &
103:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &103:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &
104:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))104:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))
105: ALLOCATE(BESTXYZ((3*NATOMS)*(INTIMAGE+2)),BESTEEE(INTIMAGE+2))105: ALLOCATE(BESTXYZ((3*NATOMS)*(INTIMAGE+2)),BESTEEE(INTIMAGE+2))
106: 106: 
107: SWITCHED=.FALSE.107: SWITCHED=.FALSE.
108: INTIMAGESAVE=INTIMAGE108: INTIMAGESAVE=INTIMAGE
109: NBACKTRACK=1109: NBACKTRACK=1
110: CALL MYCPU_TIME(STIME,.FALSE.)110: CALL MYCPU_TIME(STIME,.FALSE.)
111: WRITE(MYUNIT,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1111: WRITE(MYUNIT,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1
112: WRITE(MYUNIT,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',MUPDATE,' maximum step size=',MAXBFGS112: WRITE(MYUNIT,'(A,I6)') ' intlbfgs> Updates: ',MUPDATE
113: ADDATOM=.FALSE.113: ADDATOM=.FALSE.
114: NFAIL=0114: NFAIL=0
115: IMGFREEZE(1:INTIMAGE)=.FALSE.115: IMGFREEZE(1:INTIMAGE)=.FALSE.
116: D=(3*NATOMS)*INTIMAGE116: D=(3*NATOMS)*INTIMAGE
117: U=MUPDATE117: U=MUPDATE
118: NITERDONE=1118: NITERDONE=1
119: NITERUSE=1119: NITERUSE=1
120: NQDONE=0120: NQDONE=0
121: 121: 
122: IF ( D<=0 ) THEN122: IF ( D<=0 ) THEN
667:             IF (DUMMY.GT.DMAX) THEN667:             IF (DUMMY.GT.DMAX) THEN
668:                DMAX=DUMMY668:                DMAX=DUMMY
669:                JMAX=J1669:                JMAX=J1
670:             ENDIF670:             ENDIF
671:             IF (DUMMY.LT.DMIN) THEN671:             IF (DUMMY.LT.DMIN) THEN
672:                DMIN=DUMMY672:                DMIN=DUMMY
673:                JMIN=J1673:                JMIN=J1
674:             ENDIF674:             ENDIF
675: !            IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &675: !            IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &
676: !  &                                                  J1,' and ',J1+1,' is ',DUMMY676: !  &                                                  J1,' and ',J1+1,' is ',DUMMY
677: !            IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &677: !!           IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &
678: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1678: !! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1
679:          ENDDO679:          ENDDO
680:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &680:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &
681:   &                                                  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
682:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &682:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &
683:   &                                                  DMAX,' for images ',JMAX,JMAX+1683:   &                                                  DMAX,' for images ',JMAX,JMAX+1
684:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> smallest image separation is ', &684:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> smallest image separation is ', &
685:   &                                                  DMIN,' for images ',JMIN,JMIN+1685:   &                                                  DMIN,' for images ',JMIN,JMIN+1
686:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)686:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)
687: !        IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN687: !        IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN
688:          IF ((SQRT(ADMAX).GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN688:          IF ((SQRT(ADMAX).GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN
878:             IF (DUMMY.GT.DMAX) THEN878:             IF (DUMMY.GT.DMAX) THEN
879:                DMAX=DUMMY879:                DMAX=DUMMY
880:                JMAX=J1880:                JMAX=J1
881:             ENDIF881:             ENDIF
882:             IF (DUMMY.LT.DMIN) THEN882:             IF (DUMMY.LT.DMIN) THEN
883:                DMIN=DUMMY883:                DMIN=DUMMY
884:                JMIN=J1884:                JMIN=J1
885:             ENDIF885:             ENDIF
886: !            IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &886: !            IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &
887: !  &                                                  J1,' and ',J1+1,' is ',DUMMY887: !  &                                                  J1,' and ',J1+1,' is ',DUMMY
888: !            IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &888: !!           IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &
889: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1889: !! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1
890:          ENDDO890:          ENDDO
891:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &891:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &
892:   &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE892:   &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE
893:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &893:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &
894:   &                                                  DMAX,' for images ',JMAX,JMAX+1894:   &                                                  DMAX,' for images ',JMAX,JMAX+1
895:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)895:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)
896: !        IF (SQRT(ADMAX).GT.IMSEPMAX) THEN896: !        IF (SQRT(ADMAX).GT.IMSEPMAX) THEN
897: !           KINT=MIN(1.0D6,KINT*1.1D0)897: !           KINT=MIN(1.0D6,KINT*1.1D0)
898: !        ELSE898: !        ELSE
899: !           KINT=MAX(1.0D-6,KINT/1.1D0)899: !           KINT=MAX(1.0D-6,KINT/1.1D0)
900: !        ENDIF900: !        ENDIF
901: !        WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT901:          WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT
902:       ENDIF902:       ENDIF
903:    ENDIF903:    ENDIF
904: !904: !
905: ! End of add/subtract images block.905: ! End of add/subtract images block.
906: !906: !
907:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN907:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN
908:       LDEBUG=.FALSE.908:       LDEBUG=.FALSE.
909:       DO J2=2,INTIMAGE+2909:       DO J2=2,INTIMAGE+2
910:          CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &910:          CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
911:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DIST,DIST2,RIGID,RMAT)911:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DIST,DIST2,RIGID,RMAT)
912:       ENDDO912:       ENDDO
913:    ENDIF913:    ENDIF
914: 914: 
915:    IF (.NOT.SWITCHED) THEN915:    IF (.NOT.SWITCHED) THEN
916:       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)
 917: !     IF (INTIMAGE.GT.300) THEN
 918: !        WRITE(MYUNIT,'(A,2L5)') 'atom 375 intfrozen and atomactive: ',INTFROZEN(375),ATOMACTIVE(375)
 919: !        WRITE(MYUNIT,'(A,2L5)') 'atom 384 intfrozen and atomactive: ',INTFROZEN(384),ATOMACTIVE(384)
 920: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 284:',XYZ(3*400*283+3*374+1:3*400*283+3*374+3)
 921: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 285:',XYZ(3*400*284+3*374+1:3*400*284+3*374+3)
 922: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 286:',XYZ(3*400*285+3*374+1:3*400*285+3*374+3)
 923: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 284:',XYZ(3*400*283+3*383+1:3*400*283+3*383+3)
 924: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 285:',XYZ(3*400*284+3*383+1:3*400*284+3*383+3)
 925: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 286:',XYZ(3*400*285+3*383+1:3*400*285+3*383+3)
 926: !     ENDIF
917:       IF (QCIADDREP.GT.0) THEN927:       IF (QCIADDREP.GT.0) THEN
918:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)928:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
919:       ELSEIF (CHECKCONINT) THEN929:       ELSEIF (CHECKCONINT) THEN
920:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)930:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
921:       ELSE931:       ELSE
922:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)932:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
923:       ENDIF933:       ENDIF
 934: !     IF (INTIMAGE.GT.300) THEN
 935: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 284 GGG:',GGG(3*400*283+3*374+1:3*400*283+3*374+3)
 936: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 285 GGG:',GGG(3*400*284+3*374+1:3*400*284+3*374+3)
 937: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 286 GGG:',GGG(3*400*285+3*374+1:3*400*285+3*374+3)
 938: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 284 GGG:',GGG(3*400*283+3*383+1:3*400*283+3*383+3)
 939: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 285 GGG:',GGG(3*400*284+3*383+1:3*400*284+3*383+3)
 940: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 286 GGG:',GGG(3*400*285+3*383+1:3*400*285+3*383+3)
 941: !     ENDIF
924: 942: 
925:       IF ((ETOTAL-EOLD.LT.1.0D100).OR.ADDATOM) THEN ! MAXERISE effectively set to 1.0D100 here943:       IF ((ETOTAL-EOLD.LT.1.0D100).OR.ADDATOM) THEN ! MAXERISE effectively set to 1.0D100 here
926:          EOLD=ETOTAL944:          EOLD=ETOTAL
927:          GLAST(1:D)=G(1:D)945:          GLAST(1:D)=G(1:D)
928:          XSAVE(1:D)=X(1:D)946:          XSAVE(1:D)=X(1:D)
929:       ELSE947:       ELSE
930:          NDECREASE=NDECREASE+1948:          NDECREASE=NDECREASE+1
931:          IF (NDECREASE.GT.5) THEN949:          IF (NDECREASE.GT.5) THEN
932:             NFAIL=NFAIL+1950:             NFAIL=NFAIL+1
933:             WRITE(*,'(A,I6)') ' intlbfgs> WARNING *** in lbfgs cannot find a lower energy, NFAIL=',NFAIL951:             WRITE(*,'(A,I6)') ' intlbfgs> WARNING *** in lbfgs cannot find a lower energy, NFAIL=',NFAIL
970:       ENDIF988:       ENDIF
971:       ETOTAL=USEFRAC*ETOTALTMP+(1.0D0-USEFRAC)*ETOTAL989:       ETOTAL=USEFRAC*ETOTALTMP+(1.0D0-USEFRAC)*ETOTAL
972:       G(1:D)=USEFRAC*MYGTMP(1:D)+(1.0D0-USEFRAC)*G(1:D)990:       G(1:D)=USEFRAC*MYGTMP(1:D)+(1.0D0-USEFRAC)*G(1:D)
973:       RMS=SUM(G(1:D)**2)991:       RMS=SUM(G(1:D)**2)
974:       RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))992:       RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))
975:       EEE(1:INTIMAGE+2)=USEFRAC*EEETMP(1:INTIMAGE+2)+(1.0D0-USEFRAC)*EEE(1:INTIMAGE+2)993:       EEE(1:INTIMAGE+2)=USEFRAC*EEETMP(1:INTIMAGE+2)+(1.0D0-USEFRAC)*EEE(1:INTIMAGE+2)
976:       WORST=-1.0D100994:       WORST=-1.0D100
977:       DO J4=2,INTIMAGE+1995:       DO J4=2,INTIMAGE+1
978:          IF (EEE(J4).GT.WORST) WORST=EEE(J4)996:          IF (EEE(J4).GT.WORST) WORST=EEE(J4)
979:       ENDDO997:       ENDDO
980:       IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I8)') 'intlbfgs> Highest QCI image energy=',WORST,' images=',INTIMAGE998:       WRITE(MYUNIT,'(A,G20.10,A,I8)') 'intlbfgs> Highest QCI image energy=',WORST,' images=',INTIMAGE
981:    ENDIF999:    ENDIF
982:    IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN1000:    IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN
983:       WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/INTIMAGE,COLDFUSIONLIMIT1001:       WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/INTIMAGE,COLDFUSIONLIMIT
984:       DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)1002:       DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)
985:       DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &1003:       DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
986:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)1004:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
987:       INTIMAGE=INTIMAGESAVE1005:       INTIMAGE=INTIMAGESAVE
988:       RETURN1006:       RETURN
989:    ENDIF1007:    ENDIF
990: 1008: 
1117:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1135:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1118:          ELSE1136:          ELSE
1119:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1137:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1120:          ENDIF1138:          ENDIF
1121:       ENDIF1139:       ENDIF
1122:       LASTGOODE=ETOTAL1140:       LASTGOODE=ETOTAL
1123:    ENDIF1141:    ENDIF
1124: 1142: 
1125:    EXITSTATUS=01143:    EXITSTATUS=0
1126:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1144:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1127:    1145: !  IF ((.NOT.SWITCHED).AND.(RMS<=INTRMSTOL*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))).AND.NITERDONE>1) EXITSTATUS=1 
1128:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 1146:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 
 1147: !  IF (SWITCHED.AND.(RMS<=CQMAX).AND.NITERDONE>1) EXITSTATUS=1 
1129:    IF (SWITCHED.AND.(MAXRMS<=CQMAX).AND.NITERDONE>1) EXITSTATUS=1 1148:    IF (SWITCHED.AND.(MAXRMS<=CQMAX).AND.NITERDONE>1) EXITSTATUS=1 
1130:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21149:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
 1150: !  IF (SQRT(ADMAX).GT.IMSEPMAX) EXITSTATUS=0 ! prevent converge if largest atomic displacement is too big
1131:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1151:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
1132: !  WRITE(MYUNIT,'(A,2G20.10,3I8)') 'MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX 
1133: 1152: 
1134:    IF (EXITSTATUS > 0) THEN  1153:    IF (EXITSTATUS > 0) THEN  
1135:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1154:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1136: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771155: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1137:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771156:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1138:          IF (NACTIVE.LT.NATOMS) THEN 1157:          IF (NACTIVE.LT.NATOMS) THEN 
1139:             ADDATOM=.TRUE.1158:             ADDATOM=.TRUE.
1140:             GOTO 7771159:             GOTO 777
1141:          ENDIF1160:          ENDIF
1142:          CALL MYCPU_TIME(FTIME,.FALSE.)1161:          CALL MYCPU_TIME(FTIME,.FALSE.)
1223: ! Linear interpolation for constraint potential and real potential separately.1242: ! Linear interpolation for constraint potential and real potential separately.
1224: ! Constraint potential need not be flat if we have done some steps with both1243: ! Constraint potential need not be flat if we have done some steps with both
1225: ! potentials turned on.1244: ! potentials turned on.
1226: !1245: !
1227: 678 CONTINUE1246: 678 CONTINUE
1228: 1247: 
1229: CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)1248: CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)
1230: CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)1249: CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)
1231: NQDONE=NQDONE+11250: NQDONE=NQDONE+1
1232: 1251: 
1233: IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> WORST=',WORST1252: WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> WORST=',WORST
1234: IF (DEBUG) WRITE(MYUNIT,'(A,2I8)') 'intlbfgs> NQDONE,MCSTEPS=',NQDONE,MCSTEPS(1)1253: WRITE(MYUNIT,'(A,2I8)') 'intlbfgs> NQDONE,MCSTEPS=',NQDONE,MCSTEPS(1)
1235: IF (WORST.EQ.0.0D0) GOTO 87651254: IF (WORST.EQ.0.0D0) GOTO 8765
1236: IF (NQDONE.EQ.MCSTEPS(1)) GOTO 87651255: IF (NQDONE.EQ.MCSTEPS(1)) GOTO 8765
1237: 1256: 
1238: !1257: !
1239: ! Accept/reject this QCI set until BH steps exceeded, or worst energy is zero.1258: ! Accept/reject this QCI set until BH steps exceeded, or worst energy is zero.
1240: ! Reset everything necessary and go back to 9876 CONTINUE - start of QCI optimisation.1259: ! Reset everything necessary and go back to 9876 CONTINUE - start of QCI optimisation.
1241: !1260: !
1242: IF (WORST.GT.BESTWORST) THEN1261: IF (WORST.GT.BESTWORST) THEN
1243:    WRITE(MYUNIT,'(A)') 'intlbfgs> rejecting step - resetting'1262:    WRITE(MYUNIT,'(A)') 'intlbfgs> rejecting step - resetting'
1244:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &1263:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
1750: 543         CONTINUE1769: 543         CONTINUE
1751:          ENDDO1770:          ENDDO
1752:          ATOMACTIVE(NEWATOM)=.TRUE.1771:          ATOMACTIVE(NEWATOM)=.TRUE.
1753:          NACTIVE=NACTIVE+11772:          NACTIVE=NACTIVE+1
1754: 1773: 
1755: !1774: !
1756: ! Freeze atoms that became active more than NACTIVE-MAXNACTIVE events ago.1775: ! Freeze atoms that became active more than NACTIVE-MAXNACTIVE events ago.
1757: ! For example, with MAXNACTIVE=5 and 40 active atoms, we would freeze those 1776: ! For example, with MAXNACTIVE=5 and 40 active atoms, we would freeze those 
1758: ! turned on first, second, up to the 35th in the TURNONORDER list.1777: ! turned on first, second, up to the 35th in the TURNONORDER list.
1759: !1778: !
1760:          IF (DEBUG) WRITE(MYUNIT,'(A,I6)') 'doaddatom> Number of active atoms is now ',NACTIVE1779:          WRITE(MYUNIT,'(A,I6)') 'doaddatom> Number of active atoms is now ',NACTIVE
1761:          IF (NACTIVE.GT.MAXNACTIVE) THEN1780:          IF (NACTIVE.GT.MAXNACTIVE) THEN
1762:             WRITE(MYUNIT,'(A)') 'doaddatom> TURNONORDER:'1781:             WRITE(MYUNIT,'(A)') 'doaddatom> TURNONORDER:'
1763:             WRITE(MYUNIT,'(5I6)') TURNONORDER(1:NACTIVE-1)1782:             WRITE(MYUNIT,'(5I6)') TURNONORDER(1:NACTIVE-1)
1764:             NDUMMY=TURNONORDER(NACTIVE-MAXNACTIVE)1783:             NDUMMY=TURNONORDER(NACTIVE-MAXNACTIVE)
1765:             IF (INTFROZEN(NDUMMY)) THEN1784:             IF (INTFROZEN(NDUMMY)) THEN
1766:                IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,2I6)') ' doaddatom> Not turning off frozen active atom ',NDUMMY,' already frozen'1785:                IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,2I6)') ' doaddatom> Not turning off frozen active atom ',NDUMMY,' already frozen'
1767:             ELSE1786:             ELSE
1768:                IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,2I6)') ' doaddatom> Freezing active atom ',NDUMMY1787:                IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,2I6)') ' doaddatom> Freezing active atom ',NDUMMY
1769:                INTFROZEN(NDUMMY)=.TRUE.1788:                INTFROZEN(NDUMMY)=.TRUE.
1770: !1789: !
2251: IF (QCIAMBERT) THEN             2270: IF (QCIAMBERT) THEN             
2252:    CALL TOPOLOGY_READER(NBOND)  2271:    CALL TOPOLOGY_READER(NBOND)  
2253: !2272: !
2254: !  kr366> assume we use two endpoints and topology for amber constraints2273: !  kr366> assume we use two endpoints and topology for amber constraints
2255: !  get number of bonds and bonds from topology2274: !  get number of bonds and bonds from topology
2256: !  loop through all bonds and add them to constraint list2275: !  loop through all bonds and add them to constraint list
2257: !2276: !
2258:    DO J2=1,NBOND                !loop through all bonds and add them to constraint list2277:    DO J2=1,NBOND                !loop through all bonds and add them to constraint list
2259:       IF (INTFROZEN(BONDS(J2,1)).AND.INTFROZEN(BONDS(J2,2))) CYCLE ! no constraints between intfrozen atoms2278:       IF (INTFROZEN(BONDS(J2,1)).AND.INTFROZEN(BONDS(J2,2))) CYCLE ! no constraints between intfrozen atoms
2260:       NCONSTRAINT=NCONSTRAINT+12279:       NCONSTRAINT=NCONSTRAINT+1
2261:       IF (DEBUG) WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',BONDS(J2,1),BONDS(J2,2), &2280:       WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',BONDS(J2,1),BONDS(J2,2), &
2262:   &                     '  total=',NCONSTRAINT2281:   &                     '  total=',NCONSTRAINT
2263:       DS=SQRT((LXYZ(3*(BONDS(J2,1)-1)+1)-LXYZ(3*(BONDS(J2,2)-1)+1))**2 &2282:       DS=SQRT((LXYZ(3*(BONDS(J2,1)-1)+1)-LXYZ(3*(BONDS(J2,2)-1)+1))**2 &
2264:   &          +(LXYZ(3*(BONDS(J2,1)-1)+2)-LXYZ(3*(BONDS(J2,2)-1)+2))**2 &2283:   &          +(LXYZ(3*(BONDS(J2,1)-1)+2)-LXYZ(3*(BONDS(J2,2)-1)+2))**2 &
2265:   &          +(LXYZ(3*(BONDS(J2,1)-1)+3)-LXYZ(3*(BONDS(J2,2)-1)+3))**2)2284:   &          +(LXYZ(3*(BONDS(J2,1)-1)+3)-LXYZ(3*(BONDS(J2,2)-1)+3))**2)
2266:       DF=SQRT((LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+1)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+1))**2 &2285:       DF=SQRT((LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+1)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+1))**2 &
2267:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+2)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+2))**2 &2286:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+2)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+2))**2 &
2268:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+3)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+3))**2)2287:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+3)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+3))**2)
2269:       IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2288:       IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE
2270:       CONI(NCONSTRAINT)=MIN(BONDS(J2,1),BONDS(J2,2))2289:       CONI(NCONSTRAINT)=MIN(BONDS(J2,1),BONDS(J2,2))
2271:       CONJ(NCONSTRAINT)=MAX(BONDS(J2,2),BONDS(J2,2))2290:       CONJ(NCONSTRAINT)=MAX(BONDS(J2,2),BONDS(J2,2))
2272:       CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02291:       CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2273:       CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02292:       CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2274:       IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2293:       IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2275:       IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2294:       IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2276:    ENDDO2295:    ENDDO
2277:    WRITE(MYUNIT,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for ',CONI(NCONSTRAINT), &2296:    WRITE(MYUNIT,'(A,2I6,A,3F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &
2278:   &             CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &2297:   &             CONJ(NCONSTRAINT),' values are ',DS,DF,DUMMY,' fraction=',2*ABS(DS-DF)/(DS+DF), &
2279:   &            ' # bond constraints=',NCONSTRAINT2298:   &            ' # bond constraints=',NCONSTRAINT
2280:    QCIBONDS=NCONSTRAINT2299:    QCIBONDS=NCONSTRAINT
2281: !2300: !
2282: ! Add constraints for second-nearest neighbours - should correspond to bond angles2301: ! Add constraints for second-nearest neighbours - should correspond to bond angles
2283: !2302: !
2284:    DO J2=1,NBOND2303:    DO J2=1,NBOND
2285:       inloop: DO J3=J2+1,NBOND2304:       inloop: DO J3=J2+1,NBOND
2286:         IF (BONDS(J2,1).EQ.BONDS(J3,1)) THEN2305:         IF (BONDS(J2,1).EQ.BONDS(J3,1)) THEN
2287:            ATOM1=BONDS(J2,2)2306:            ATOM1=BONDS(J2,2)
2288:            ATOM2=BONDS(J3,2)2307:            ATOM2=BONDS(J3,2)
2293:            ATOM1=BONDS(J2,1)2312:            ATOM1=BONDS(J2,1)
2294:            ATOM2=BONDS(J3,2)2313:            ATOM2=BONDS(J3,2)
2295:         ELSEIF (BONDS(J2,2).EQ.BONDS(J3,2)) THEN2314:         ELSEIF (BONDS(J2,2).EQ.BONDS(J3,2)) THEN
2296:            ATOM1=BONDS(J2,1)2315:            ATOM1=BONDS(J2,1)
2297:            ATOM2=BONDS(J3,1)2316:            ATOM2=BONDS(J3,1)
2298:         ELSE2317:         ELSE
2299:            CYCLE inloop2318:            CYCLE inloop
2300:         ENDIF2319:         ENDIF
2301:         IF (INTFROZEN(ATOM1).AND.INTFROZEN(ATOM2)) CYCLE ! no constraints between intfrozen atoms2320:         IF (INTFROZEN(ATOM1).AND.INTFROZEN(ATOM2)) CYCLE ! no constraints between intfrozen atoms
2302:         NCONSTRAINT=NCONSTRAINT+12321:         NCONSTRAINT=NCONSTRAINT+1
2303: !       WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for second neighbours ',ATOM1,ATOM2, &2322:         WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for second neighbours ',ATOM1,ATOM2, &
2304: ! &                     '  total=',NCONSTRAINT2323:   &                     '  total=',NCONSTRAINT
2305:          DS=SQRT((LXYZ(3*(ATOM1-1)+1)-LXYZ(3*(ATOM2-1)+1))**2 &2324:          DS=SQRT((LXYZ(3*(ATOM1-1)+1)-LXYZ(3*(ATOM2-1)+1))**2 &
2306:   &             +(LXYZ(3*(ATOM1-1)+2)-LXYZ(3*(ATOM2-1)+2))**2 &2325:   &             +(LXYZ(3*(ATOM1-1)+2)-LXYZ(3*(ATOM2-1)+2))**2 &
2307:   &             +(LXYZ(3*(ATOM1-1)+3)-LXYZ(3*(ATOM2-1)+3))**2)2326:   &             +(LXYZ(3*(ATOM1-1)+3)-LXYZ(3*(ATOM2-1)+3))**2)
2308:          DF=SQRT((LXYZ(3*NATOMS+3*(ATOM1-1)+1)-LXYZ(3*NATOMS+3*(ATOM2-1)+1))**2 &2327:          DF=SQRT((LXYZ(3*NATOMS+3*(ATOM1-1)+1)-LXYZ(3*NATOMS+3*(ATOM2-1)+1))**2 &
2309:   &             +(LXYZ(3*NATOMS+3*(ATOM1-1)+2)-LXYZ(3*NATOMS+3*(ATOM2-1)+2))**2 &2328:   &             +(LXYZ(3*NATOMS+3*(ATOM1-1)+2)-LXYZ(3*NATOMS+3*(ATOM2-1)+2))**2 &
2310:   &             +(LXYZ(3*NATOMS+3*(ATOM1-1)+3)-LXYZ(3*NATOMS+3*(ATOM2-1)+3))**2)2329:   &             +(LXYZ(3*NATOMS+3*(ATOM1-1)+3)-LXYZ(3*NATOMS+3*(ATOM2-1)+3))**2)
2311:          IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2330:          IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE
2312:          CONI(NCONSTRAINT)=MIN(ATOM1,ATOM2)2331:          CONI(NCONSTRAINT)=MIN(ATOM1,ATOM2)
2313:          CONJ(NCONSTRAINT)=MAX(ATOM1,ATOM2)2332:          CONJ(NCONSTRAINT)=MAX(ATOM1,ATOM2)
2314:          CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02333:          CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2315:          CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02334:          CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2316:          IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2335:          IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2317:          IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2336:          IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2318: !        WRITE(MYUNIT,'(A,2I6,A,2F12.2,A,F12.4,A,2I8)') ' intlbfgs> constrain distance for ',CONI(NCONSTRAINT), &2337:          WRITE(MYUNIT,'(A,2I6,A,3F12.2,A,F12.4,A,2I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &
2319: ! &             CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &2338:   &             CONJ(NCONSTRAINT),' values are ',DS,DF,DUMMY,' fraction=',2*ABS(DS-DF)/(DS+DF), &
2320: ! &            ' # second neighbour constraints, total=',QCISECOND,NCONSTRAINT2339:   &            ' # second neighbour constraints and total=',QCISECOND,NCONSTRAINT
2321:       ENDDO inloop2340:       ENDDO inloop
2322:    ENDDO2341:    ENDDO
2323:    QCISECOND=NCONSTRAINT-QCIBONDS2342:    QCISECOND=NCONSTRAINT-QCIBONDS
2324:    IF (CONFILET) THEN 2343: ENDIF
2325:       LUNIT=GETUNIT()2344: IF (CONFILET) THEN 
2326:       OPEN(LUNIT,FILE='constraintfile',STATUS='OLD') 
2327: ! 
2328: !  Additional amber constraints, e.g. cis/trans 
2329: ! 
2330:       DO  
2331:          READ(LUNIT,*,END=534)  J2, J3 
2332: ! 
2333: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to  
2334: ! turn this off 
2335: ! 
2336:          IF (J3-J2.GT.INTCONSEP) CYCLE  
2337:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms 
2338:          NCONSTRAINT=NCONSTRAINT+1 
2339:          WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding extra constraint for atoms ',J2,J3,'  total=',NCONSTRAINT 
2340:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 & 
2341:   &             +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 & 
2342:   &             +(LXYZ(3*(J2-1)+3)-LXYZ(3*(J3-1)+3))**2)  
2343:          DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 & 
2344:   &             +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 & 
2345:   &             +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2)  
2346:          IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE 
2347:          CONI(NCONSTRAINT)=J2 
2348:          CONJ(NCONSTRAINT)=J3 
2349:          CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0 
2350:          CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0 
2351:          IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT) 
2352:          IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT) 
2353:          WRITE(MYUNIT,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for ',CONI(NCONSTRAINT), & 
2354:   &                     CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), & 
2355:   &                  ' # constraints=',NCONSTRAINT 
2356:       ENDDO 
2357: 534   CONTINUE 
2358:       WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, & 
2359:   &                                     ' shortest=',MINCONDIST,' longest=',MAXCONDIST 
2360:       CLOSE(LUNIT) 
2361:    ENDIF 
2362: ELSE IF (CONFILET) THEN  
2363:     LUNIT=GETUNIT()2345:     LUNIT=GETUNIT()
2364:     OPEN(LUNIT,FILE='constraintfile',STATUS='OLD')2346:     OPEN(LUNIT,FILE='constraintfile',STATUS='OLD')
2365: !2347: !
2366: !  Add constraint for this distance to the list.2348: !  Add constraint for this distance to the list.
2367: !2349: !
2368:     DO 2350:     DO 
2369:        READ(LUNIT,*,END=531)  J2, J32351:        READ(LUNIT,*,END=531)  J2, J3
2370: !2352: !
2371: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to 2353: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to 
2372: ! turn this off2354: ! turn this off
2381:        DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 &2363:        DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 &
2382:   &           +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 &2364:   &           +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 &
2383:   &           +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2) 2365:   &           +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2) 
2384:        IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2366:        IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE
2385:        CONI(NCONSTRAINT)=J22367:        CONI(NCONSTRAINT)=J2
2386:        CONJ(NCONSTRAINT)=J32368:        CONJ(NCONSTRAINT)=J3
2387:        CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02369:        CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2388:        CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02370:        CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2389:        IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2371:        IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2390:        IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2372:        IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2391:        WRITE(MYUNIT,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for ',CONI(NCONSTRAINT), &2373:        WRITE(MYUNIT,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &
2392:   &                 CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &2374:   &                 CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &
2393:   &                ' # constraints=',NCONSTRAINT2375:   &                ' # constraints=',NCONSTRAINT
2394:     ENDDO2376:     ENDDO
2395: 531 CONTINUE2377: 531 CONTINUE
2396:     CLOSE(LUNIT) 
2397:     WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &2378:     WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &
2398:   &                               ' shortest=',MINCONDIST,' longest=',MAXCONDIST2379:   &                               ' shortest=',MINCONDIST,' longest=',MAXCONDIST
2399: 2380: 
2400: ELSE IF (NCONGEOM.LT.2) THEN 2381: ELSE IF (NCONGEOM.LT.2) THEN 
2401:    DO J2=1,NATOMS2382:    DO J2=1,NATOMS
2402:       DO J3=J2+1,NATOMS2383:       DO J3=J2+1,NATOMS
2403: 2384: 
2404:          IF (J3-J2.GT.INTCONSEP) CYCLE ! forbid constraints corresponding to atoms distant in sequence2385:          IF (J3-J2.GT.INTCONSEP) CYCLE ! forbid constraints corresponding to atoms distant in sequence
2405:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms2386:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms
2406:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &2387:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &
2419: !2400: !
2420:             NCONSTRAINT=NCONSTRAINT+12401:             NCONSTRAINT=NCONSTRAINT+1
2421:             WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',J2,J3,'  total=',NCONSTRAINT2402:             WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',J2,J3,'  total=',NCONSTRAINT
2422:             IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2403:             IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE
2423:             CONI(NCONSTRAINT)=J22404:             CONI(NCONSTRAINT)=J2
2424:             CONJ(NCONSTRAINT)=J32405:             CONJ(NCONSTRAINT)=J3
2425:             CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02406:             CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2426:             CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02407:             CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2427:             IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2408:             IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2428:             IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2409:             IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2429: !           IF (DEBUG) PRINT '(A,2I6,A,2F12.2,A,F12.4,A,I8)',' intlbfgs> constrain distance for ',CONI(NCONSTRAINT), &2410: !           IF (DEBUG) PRINT '(A,2I6,A,2F12.2,A,F12.4,A,I8)',' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &
2430: ! &                 CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &2411: ! &                 CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &
2431: ! &                ' # constraints=',NCONSTRAINT2412: ! &                ' # constraints=',NCONSTRAINT
2432:          ENDIF2413:          ENDIF
2433:       ENDDO2414:       ENDDO
2434:    ENDDO2415:    ENDDO
2435:    IF (DEBUG) WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &2416:    IF (DEBUG) WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &
2436:   &                                     ' shortest=',MINCONDIST,' longest=',MAXCONDIST2417:   &                                     ' shortest=',MINCONDIST,' longest=',MAXCONDIST
2437: ELSE2418: ELSE
2438:    DO J2=1,NATOMS2419:    DO J2=1,NATOMS
2439:       DO J3=J2+1,NATOMS2420:       DO J3=J2+1,NATOMS
2470: 753      CONTINUE2451: 753      CONTINUE
2471:       ENDDO2452:       ENDDO
2472:    ENDDO2453:    ENDDO
2473:    CONIFIX(1:NCONSTRAINT)=CONI(1:NCONSTRAINT)2454:    CONIFIX(1:NCONSTRAINT)=CONI(1:NCONSTRAINT)
2474:    CONJFIX(1:NCONSTRAINT)=CONJ(1:NCONSTRAINT)2455:    CONJFIX(1:NCONSTRAINT)=CONJ(1:NCONSTRAINT)
2475:    CONDISTREFFIX(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)2456:    CONDISTREFFIX(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)
2476:    CONCUTFIX(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)2457:    CONCUTFIX(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)
2477:    NCONSTRAINTFIX=NCONSTRAINT2458:    NCONSTRAINTFIX=NCONSTRAINT
2478: ENDIF2459: ENDIF
2479: 2460: 
 2461: write(myunit,'(A,2I6)') 'QCIADDREP,QCIBONDS=',QCIADDREP,QCIBONDS
 2462: 
2480: IF (QCIADDREP.GT.0) THEN2463: IF (QCIADDREP.GT.0) THEN
2481:    DMIN=1.0D1002464:    DMIN=1.0D100
2482:    DO J2=1,QCIBONDS2465:    DO J2=1,QCIBONDS
2483: !2466: !
2484: ! end point 12467: ! end point 1
2485: !2468: !
2486:       NI1=3*(CONI(J2)-1)2469:       NI1=3*(CONI(J2)-1)
2487:       NJ1=3*(CONJ(J2)-1)2470:       NJ1=3*(CONJ(J2)-1)
2488:       DO J3=J2+1,QCIBONDS2471:       DO J3=J2+1,QCIBONDS
2489:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom2472:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom


r31259/keywords.f 2016-10-10 11:30:10.413214915 +0100 r31258/keywords.f 2016-10-10 11:30:13.573256904 +0100
172: ! QCI parameters172: ! QCI parameters
173: !173: !
174:          CONDATT=.FALSE.174:          CONDATT=.FALSE.
175:          QCIPOTT=.FALSE.175:          QCIPOTT=.FALSE.
176:          QCIPOT2T=.FALSE.176:          QCIPOT2T=.FALSE.
177:          QCIADDREP=0177:          QCIADDREP=0
178:          QCIADDREPCUT=1.0D0178:          QCIADDREPCUT=1.0D0
179:          QCIADDREPEPS=1.0D0179:          QCIADDREPEPS=1.0D0
180:          QCINOREPINT=.FALSE.180:          QCINOREPINT=.FALSE.
181:          MAXNACTIVE=0181:          MAXNACTIVE=0
182:          FREEZENODEST=.FALSE. 
183:          FREEZETOL=1.0D-3182:          FREEZETOL=1.0D-3
184:          QCIPERMCHECK=.FALSE.183:          QCIPERMCHECK=.FALSE.
185:          QCIPERMCHECKINT=100184:          QCIPERMCHECKINT=100
186:          INTCONSTRAINTT=.FALSE.185:          INTCONSTRAINTT=.FALSE.
187:          INTCONSTRAINTTOL=0.1D0186:          INTCONSTRAINTTOL=0.1D0
188:          INTCONSTRAINTDEL=10.0D0187:          INTCONSTRAINTDEL=10.0D0
189:          INTCONSTRAINTREP=100.0D0188:          INTCONSTRAINTREP=100.0D0
190:          INTCONSTRAINREPCUT=1.7D0189:          INTCONSTRAINREPCUT=1.7D0
191:          INTFREEZET=.FALSE.190:          INTFREEZET=.FALSE.
192:          INTFREEZETOL=1.0D-3191:          INTFREEZETOL=1.0D-3
1165: !1164: !
1166:       MLP3T=.FALSE.1165:       MLP3T=.FALSE.
1167:       MLPB3T=.FALSE.1166:       MLPB3T=.FALSE.
1168:       MLPB3NEWT=.FALSE.1167:       MLPB3NEWT=.FALSE.
1169:       MLPNEWREG=.FALSE.1168:       MLPNEWREG=.FALSE.
1170:       MLPDONE=.FALSE.1169:       MLPDONE=.FALSE.
1171:       MLPNORM=.FALSE.1170:       MLPNORM=.FALSE.
1172:       MLPLAMBDA=0.0D01171:       MLPLAMBDA=0.0D0
1173: 1172: 
1174:       LJADDT=.FALSE.1173:       LJADDT=.FALSE.
1175:  
1176:       DUMPMQT=.FALSE. 
1177:       1174:       
1178:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)1175:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)
1179:       1176:       
1180: !      OPEN (5,FILE='data',STATUS='OLD')1177: !      OPEN (5,FILE='data',STATUS='OLD')
1181: 1178: 
1182: !190   CALL INPUT(END,5)1179: !190   CALL INPUT(END,5)
1183: 190   CALL INPUT(END, DATA_UNIT)1180: 190   CALL INPUT(END, DATA_UNIT)
1184:       IF (.NOT. END) THEN1181:       IF (.NOT. END) THEN
1185:         CALL READU(WORD)1182:         CALL READU(WORD)
1186:       ENDIF1183:       ENDIF
2220:             STOP2217:             STOP
2221:          ENDIF2218:          ENDIF
2222:    2219:    
2223:          LUNIT=GETUNIT()2220:          LUNIT=GETUNIT()
2224:          OPEN(LUNIT,FILE='MLPdata',STATUS='OLD')2221:          OPEN(LUNIT,FILE='MLPdata',STATUS='OLD')
2225:          ALLOCATE(MLPDAT(MLPDATA,MLPIN),MLPOUTCOME(MLPDATA),MLPMEAN(MLPIN))2222:          ALLOCATE(MLPDAT(MLPDATA,MLPIN),MLPOUTCOME(MLPDATA),MLPMEAN(MLPIN))
2226:          MLPMEAN(1:MLPIN)=0.0D02223:          MLPMEAN(1:MLPIN)=0.0D0
2227:          DO J1=1,MLPDATA2224:          DO J1=1,MLPDATA
2228:             READ(LUNIT,*) MLPOUTCOME(J1),(DUMMY,J2=1,MLPSTART-1),MLPDAT(J1,1:MLPIN)2225:             READ(LUNIT,*) MLPOUTCOME(J1),(DUMMY,J2=1,MLPSTART-1),MLPDAT(J1,1:MLPIN)
2229:             MLPOUTCOME(J1)=MLPOUTCOME(J1)+1 ! to shift the range to 1 to 42226:             MLPOUTCOME(J1)=MLPOUTCOME(J1)+1 ! to shift the range to 1 to 4
 2227:             WRITE(MYUNIT,'(9G20.10,I8)') MLPOUTCOME(J1),MLPDAT(J1,1:MLPIN)
2230:             DO J2=1,MLPIN2228:             DO J2=1,MLPIN
2231:                MLPMEAN(J2)=MLPMEAN(J2)+ABS(MLPDAT(J1,J2))2229:                MLPMEAN(J2)=MLPMEAN(J2)+ABS(MLPDAT(J1,J2))
2232:             ENDDO2230:             ENDDO
2233:          ENDDO2231:          ENDDO
2234:          CLOSE(LUNIT)2232:          CLOSE(LUNIT)
2235:          IF (MLPNORM) THEN2233:          IF (MLPNORM) THEN
2236:             MLPMEAN(1:MLPIN)=MLPMEAN(1:MLPIN)/MLPDATA 2234:             MLPMEAN(1:MLPIN)=MLPMEAN(1:MLPIN)/MLPDATA 
2237:             WRITE(MYUNIT,'(A)') 'keyword> Rescaling inputs by mean absolute values:'2235:             WRITE(MYUNIT,'(A)') 'keyword> Rescaling inputs by mean absolute values:'
2238:             WRITE(MYUNIT,'(6G20.10)') MLPMEAN(1:MLPIN)2236:             WRITE(MYUNIT,'(6G20.10)') MLPMEAN(1:MLPIN)
2239:             DO J1=1,MLPIN2237:             DO J1=1,MLPIN
3334:          DUMPT=.TRUE.3332:          DUMPT=.TRUE.
3335: 3333: 
3336:       ELSE IF (WORD.EQ.'DUMPINT') THEN3334:       ELSE IF (WORD.EQ.'DUMPINT') THEN
3337:          CALL READI(DUMPINT)3335:          CALL READI(DUMPINT)
3338: 3336: 
3339:       ELSE IF (WORD.EQ.'DUMP_MARKOV') THEN3337:       ELSE IF (WORD.EQ.'DUMP_MARKOV') THEN
3340:          DUMP_MARKOV=.TRUE.3338:          DUMP_MARKOV=.TRUE.
3341:          CALL READI(DUMP_MARKOV_NWAIT)3339:          CALL READI(DUMP_MARKOV_NWAIT)
3342:          CALL READI(DUMP_MARKOV_NFREQ)3340:          CALL READI(DUMP_MARKOV_NFREQ)
3343: 3341: 
3344:       ELSE IF (WORD.EQ.'DUMPMQ') THEN 
3345:          DUMPMQT=.TRUE. 
3346:          MQUNIT=GETUNIT() 
3347:          OPEN(MQUNIT,FILE='dumpMQ',STATUS='UNKNOWN') 
3348:  
3349:       ELSE IF (WORD.EQ.'DZUGUTOV') THEN3342:       ELSE IF (WORD.EQ.'DZUGUTOV') THEN
3350:          DZTEST=.TRUE.3343:          DZTEST=.TRUE.
3351:          CALL READF(DZP1)3344:          CALL READF(DZP1)
3352:          CALL READF(DZP2)3345:          CALL READF(DZP2)
3353:          CALL READF(DZP3)3346:          CALL READF(DZP3)
3354:          CALL READF(DZP4)3347:          CALL READF(DZP4)
3355:          CALL READF(DZP5)3348:          CALL READF(DZP5)
3356:          CALL READF(DZP6)3349:          CALL READF(DZP6)
3357:          CALL READF(DZP7)3350:          CALL READF(DZP7)
3358: 3351: 
3792:             ENDDO3785:             ENDDO
3793:             WRITE(MYUNIT,'(A)') ' '3786:             WRITE(MYUNIT,'(A)') ' '
3794:          ENDDO3787:          ENDDO
3795:       ELSE IF (WORD.EQ.'GROUND') THEN3788:       ELSE IF (WORD.EQ.'GROUND') THEN
3796:          GROUND=.TRUE.3789:          GROUND=.TRUE.
3797: 3790: 
3798: !--------------------------------!3791: !--------------------------------!
3799: ! hk286 > Generalised Thomson    !3792: ! hk286 > Generalised Thomson    !
3800: !--------------------------------!3793: !--------------------------------!
3801: 3794: 
3802:       ELSE IF (WORD .EQ. 'GTHOMSONBIN') THEN 
3803:          GTHOMSONT = .TRUE. 
3804:          THOMSONT=.TRUE. 
3805:          RIGID = .TRUE. 
3806:          GTHOMPOT=7 
3807:          CALL READI(GTHOMSONBINN) 
3808:          CALL READF(GTHOMSONBIND) 
3809:       ELSE IF (WORD .EQ. 'GTHOMSON') THEN3795:       ELSE IF (WORD .EQ. 'GTHOMSON') THEN
3810:          GTHOMSONT = .TRUE.3796:          GTHOMSONT = .TRUE.
3811:          THOMSONT=.TRUE.3797:          THOMSONT=.TRUE.
3812:          RIGID = .TRUE.3798:          RIGID = .TRUE.
3813:          CALL READI(GTHOMMET)3799:          CALL READI(GTHOMMET)
3814:          CALL READF(GTHOMSONZ)3800:          CALL READF(GTHOMSONZ)
3815:          IF (NITEMS.GT.3) THEN3801:          IF (NITEMS.GT.3) THEN
3816:             CALL READF(GTHOMSONC)3802:             CALL READF(GTHOMSONC)
3817:          ENDIF3803:          ENDIF
3818:          IF (NITEMS.GT.4) THEN3804:          IF (NITEMS.GT.4) THEN
5419: ! NOT DOCUMENTED.5405: ! NOT DOCUMENTED.
5420: !5406: !
5421:       ELSE IF (WORD.EQ.'CALCQ') THEN5407:       ELSE IF (WORD.EQ.'CALCQ') THEN
5422:          CALCQT=.TRUE.5408:          CALCQT=.TRUE.
5423: 5409: 
5424: !5410: !
5425: !Use topology information for QCI constraints for AMBER5411: !Use topology information for QCI constraints for AMBER
5426: !5412: !
5427:       ELSE IF (WORD.EQ.'QCIAMBER') THEN5413:       ELSE IF (WORD.EQ.'QCIAMBER') THEN
5428:          QCIAMBERT=.TRUE.5414:          QCIAMBERT=.TRUE.
5429:          WRITE(MYUNIT,'(A)') ' keyword> Use topology file for constraints in QCI'5415:          WRITE(MYUNIT,'(A)') ' keyword> Use topology file for constrtaints in QCI'
5430: 5416: 
5431: !5417: !
5432: ! Distance cut-off for Coulomb interactions in AMBER (the PNM hand-coded version).5418: ! Distance cut-off for Coulomb interactions in AMBER (the PNM hand-coded version).
5433: !5419: !
5434: !      ELSE IF (WORD.EQ.'QCUTOFF') THEN5420: !      ELSE IF (WORD.EQ.'QCUTOFF') THEN
5435: !         AMCUT=.TRUE.5421: !         AMCUT=.TRUE.
5436: !         CALL READF(REALQCUTOFF)5422: !         CALL READF(REALQCUTOFF)
5437: !         QCUTOFF=1.1D0*REALQCUTOFF5423: !         QCUTOFF=1.1D0*REALQCUTOFF
5438: !5424: !
5439: ! NOT DOCUMENTED5425: ! NOT DOCUMENTED


r31259/make_conpot.f90 2016-10-10 11:30:10.661218225 +0100 r31258/make_conpot.f90 2016-10-10 11:30:13.829260302 +0100
 28: INTEGER :: J2,ISTAT,J1,J3,J4,NCPFIT,J5,NQCIFREEZE,NDUMMY,LUNIT,GETUNIT 28: INTEGER :: J2,ISTAT,J1,J3,J4,NCPFIT,J5,NQCIFREEZE,NDUMMY,LUNIT,GETUNIT
 29: INTEGER NCONFORNEWATOM 29: INTEGER NCONFORNEWATOM
 30: DOUBLE PRECISION :: NDIST, MINCOORDS(NCPFIT,3*NATOMS), DMIN, LINTCONSTRAINTTOL, & 30: DOUBLE PRECISION :: NDIST, MINCOORDS(NCPFIT,3*NATOMS), DMIN, LINTCONSTRAINTTOL, &
 31:   &                 LXYZ(6*NATOMS) 31:   &                 LXYZ(6*NATOMS)
 32: DOUBLE PRECISION :: DMOVED(NATOMS) 32: DOUBLE PRECISION :: DMOVED(NATOMS)
 33: LOGICAL ADDREP(NATOMS) 33: LOGICAL ADDREP(NATOMS)
 34: INTEGER NDIST1(NATOMS), NCYCLE, NUNCON1, DLIST(NATOMS) 34: INTEGER NDIST1(NATOMS), NCYCLE, NUNCON1, DLIST(NATOMS)
 35: LOGICAL :: YESNO, CALLED=.FALSE. 35: LOGICAL :: YESNO, CALLED=.FALSE.
 36: SAVE CALLED 36: SAVE CALLED
 37:  37: 
 38: WRITE(MYUNIT,'(A,I6,2L5,I6)') 'make_conpot> NCONGEOM,CALLED,CONDATT,NCPFIT=',NCONGEOM,CALLED,CONDATT,NCPFIT 
 39:  
 40: IF (NCONGEOM.GE.2) THEN 38: IF (NCONGEOM.GE.2) THEN
 41: ! 39: !
 42: ! If this is not the first call, and we are being passed two minima, 40: ! If this is not the first call, and we are being passed two minima,
 43: ! then we are doing an interpolation metric for a new pair of minima. 41: ! then we are doing an interpolation metric for a new pair of minima.
 44: ! We should optimise the permutational isomers on reference minimum 1 42: ! We should optimise the permutational isomers on reference minimum 1
 45: ! and then do the overall alignment with newmindist, fixing the 43: ! and then do the overall alignment with newmindist, fixing the
 46: ! permutational isomers. This should put the permutational isomers 44: ! permutational isomers. This should put the permutational isomers
 47: ! in register with the constraints, which were calculated for all 45: ! in register with the constraints, which were calculated for all
 48: ! the reference minima after aligning with the first. 46: ! the reference minima after aligning with the first.
 49: ! 47: !
198: ! The rest of this code is for initial setup. It isn't needed if CONDATT is true.196: ! The rest of this code is for initial setup. It isn't needed if CONDATT is true.
199: !197: !
200: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT198: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT
201: IF (ALLOCATED(CONDISTREFLOCAL)) DEALLOCATE(CONDISTREFLOCAL)199: IF (ALLOCATED(CONDISTREFLOCAL)) DEALLOCATE(CONDISTREFLOCAL)
202: ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))200: ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))
203: CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)201: CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)
204: IF (ALLOCATED(CONCUTLOCAL)) DEALLOCATE(CONCUTLOCAL)202: IF (ALLOCATED(CONCUTLOCAL)) DEALLOCATE(CONCUTLOCAL)
205: ALLOCATE(CONCUTLOCAL(NCONSTRAINT))203: ALLOCATE(CONCUTLOCAL(NCONSTRAINT))
206: CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)204: CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)
207: NREPULSIVEFIX=0205: NREPULSIVEFIX=0
208: WRITE(MYUNIT,*) 'make_conpot> INTCONSTRAINTREP,NREPULSIVEFIX=',INTCONSTRAINTREP,NREPULSIVEFIX 
209: IF (INTCONSTRAINTREP.EQ.0.0D0) GOTO 963206: IF (INTCONSTRAINTREP.EQ.0.0D0) GOTO 963
210: !207: !
211: ! Add repulsions to non-constrained atoms.208: ! Add repulsions to non-constrained atoms.
212: ! Note that we do not limit the number of constraints per site in this209: ! Note that we do not limit the number of constraints per site in this
213: ! routine, unlike NEB/lbfgs.f90, where the result will depend on the210: ! routine, unlike NEB/lbfgs.f90, where the result will depend on the
214: ! order in which the constraints are turned on. 211: ! order in which the constraints are turned on. 
215: !212: !
216: NDUMMY=1213: NDUMMY=1
217: DO J1=1,NATOMS214: DO J1=1,NATOMS
218: !215: !


r31259/mc.F 2016-10-10 11:30:10.945221986 +0100 r31258/mc.F 2016-10-10 11:30:14.089263761 +0100
 87:       integer :: J6, L, M, SUMSQUAREDIFF, SUMSQUAREDIFF1D, HBONDDIFF 87:       integer :: J6, L, M, SUMSQUAREDIFF, SUMSQUAREDIFF1D, HBONDDIFF
 88:       INTEGER GETUNIT 88:       INTEGER GETUNIT
 89:       INTEGER DUMPUNIQUEUNIT 89:       INTEGER DUMPUNIQUEUNIT
 90:  90: 
 91: ! khs26> Energy decomposition for AMBER 12 91: ! khs26> Energy decomposition for AMBER 12
 92:       TYPE(POT_ENE_REC_C) :: AMBER12_ENERGY_DECOMP 92:       TYPE(POT_ENE_REC_C) :: AMBER12_ENERGY_DECOMP
 93:  93: 
 94: ! ds656 > energy before and after BGUPTA swap (without relaxation) 94: ! ds656 > energy before and after BGUPTA swap (without relaxation)
 95:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0 95:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0
 96:  96: 
 97:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6 
 98:  
 99:       QNEWRES=0 97:       QNEWRES=0
100:       QGCBH=0 98:       QGCBH=0
101:  99: 
102: ! sf344> selectively unfreeze the last added building blocks (MODULARNINCR)100: ! sf344> selectively unfreeze the last added building blocks (MODULARNINCR)
103: !        for sequential MC loops in modular global optimization101: !        for sequential MC loops in modular global optimization
104:       IF(MODULART) THEN102:       IF(MODULART) THEN
105:         FROZEN(:)=.TRUE.103:         FROZEN(:)=.TRUE.
106:            IF(MODULARCURRENTN==MODULARNMIN) THEN   ! unfreeze every atom in the first step of MODULARNMIN molecules104:            IF(MODULARCURRENTN==MODULARNMIN) THEN   ! unfreeze every atom in the first step of MODULARNMIN molecules
107:                 J6=0105:                 J6=0
108:            ELSE106:            ELSE
2144: #else2142: #else
2145:             IF (HIT) GOTO 372143:             IF (HIT) GOTO 37
2146:             IF (EPREV(JP)<EBEST(JP)) EBEST(JP)=EPREV(JP)2144:             IF (EPREV(JP)<EBEST(JP)) EBEST(JP)=EPREV(JP)
2147:             IF ((REPMATCHT).AND.(ABS(MAXVAL(EBEST(:))-MINVAL(EBEST(:)))<ECONV)) GOTO 372145:             IF ((REPMATCHT).AND.(ABS(MAXVAL(EBEST(:))-MINVAL(EBEST(:)))<ECONV)) GOTO 37
2148: #endif2146: #endif
2149:             IF (DUMPINT.GT.0) THEN2147:             IF (DUMPINT.GT.0) THEN
2150:                IF (MOD(J1,DUMPINT).EQ.0) THEN2148:                IF (MOD(J1,DUMPINT).EQ.0) THEN
2151:                   CALL DUMPSTATE(NQ(JP),EBEST,BESTCOORDS,JBEST,JP)2149:                   CALL DUMPSTATE(NQ(JP),EBEST,BESTCOORDS,JBEST,JP)
2152:                ENDIF2150:                ENDIF
2153:             ENDIF2151:             ENDIF
2154: ! 
2155: ! Dump current energy, Q4 and Q6 parameters for minimum in Markov chain. 
2156: ! 
2157:             IF (DUMPMQT) THEN  
2158:                DO J2=1,NATOMS 
2159:                   QTEMP(1:3,J2)=COORDSO(3*(J2-1)+1:3*(J2-1)+3,JP) 
2160:                ENDDO 
2161:                CALL POTENTIAL(QTEMP,GRAD,DUMMY,.TRUE.,.FALSE.) 
2162:                WRITE(MQUNIT,'(A,G20.10)') 'energy from coordso=',DUMMY 
2163:                CALL QORDER_LJ(QTEMP,Q4,Q6) 
2164:                WRITE(MQUNIT,'(A,I10,A,G20.10,A,G20.10,A,G20.10)') 'MC step ',J1,' Markov energy ', EPREV(JP), ' Q4 ',Q4,' Q6 ',Q6 
2165:                DO J2=1,NATOMS 
2166:                   WRITE(MQUNIT,'(3G20.10)') QTEMP(1:3,J2) 
2167:                ENDDO 
2168:             ENDIF 
2169:             2152:             
2170:             IF (DUMP_MARKOV) THEN 2153:             IF (DUMP_MARKOV) THEN 
2171:                ! ds656> Dump current Markov state to XYZ for post-processing 2154:                ! ds656> Dump current Markov state to XYZ for post-processing 
2172:                ! (e.g. characterisation of equilibrium structure)2155:                ! (e.g. characterisation of equilibrium structure)
2173:                J2=J1-1-DUMP_MARKOV_NWAIT2156:                J2=J1-1-DUMP_MARKOV_NWAIT
2174:                IF(J2.GE.0 .AND. 2157:                IF(J2.GE.0 .AND. 
2175:      &              MOD(J2,DUMP_MARKOV_NFREQ).EQ.0) THEN2158:      &              MOD(J2,DUMP_MARKOV_NFREQ).EQ.0) THEN
2176:                   WRITE(DUMPXYZUNIT(JP),'(I4)') NATOMS2159:                   WRITE(DUMPXYZUNIT(JP),'(I4)') NATOMS
2177:                   WRITE(DUMPXYZUNIT(JP),'(A,I9,A,F20.10)') 2160:                   WRITE(DUMPXYZUNIT(JP),'(A,I9,A,F20.10)') 
2178:      &                 'MC step ',J1,' Markov energy ', EPREV(JP)2161:      &                 'MC step ',J1,' Markov energy ', EPREV(JP)
2277:       ENDDO2260:       ENDDO
2278: #ifdef MPI2261: #ifdef MPI
2279:       IF (MPIT) THEN2262:       IF (MPIT) THEN
2280: ! sf344> Apparently MPI_FINALIZE can only be called once during a run, so we call it at the end of main.F2263: ! sf344> Apparently MPI_FINALIZE can only be called once during a run, so we call it at the end of main.F
2281: !         CALL MPI_FINALIZE(MPIERR)2264: !         CALL MPI_FINALIZE(MPIERR)
2282: !         IF (DEBUG) WRITE(MYUNIT, '(A,I6)') 'In mc after MPI_FINALIZE MPIERR=',MPIERR2265: !         IF (DEBUG) WRITE(MYUNIT, '(A,I6)') 'In mc after MPI_FINALIZE MPIERR=',MPIERR
2283:       ENDIF2266:       ENDIF
2284: #endif2267: #endif
2285: 2268: 
2286: 37    CONTINUE2269: 37    CONTINUE
2287: ! 
2288: ! Dump current energy, Q4 and Q6 paramaters for minimum in Markov chain. 
2289: ! 
2290:             IF (DUMPMQT) THEN  
2291:                DO J2=1,NATOMS 
2292:                   QTEMP(1:3,J2)=COORDS(3*(J2-1)+1:3*(J2-1)+3,JP) 
2293:                ENDDO 
2294:                CALL POTENTIAL(QTEMP,GRAD,DUMMY,.TRUE.,.FALSE.) 
2295:                WRITE(MQUNIT,'(A,G20.10)') 'energy from coords=',DUMMY 
2296:                CALL QORDER_LJ(QTEMP,Q4,Q6) 
2297:                WRITE(MQUNIT,'(A,I10,A,G20.10,A,G20.10,A,G20.10)') 'MC step ',J1,' Markov energy ', EPREV(JP), ' Q4 ',Q4,' Q6 ',Q6 
2298:                DO J2=1,NATOMS 
2299:                   WRITE(MQUNIT,'(3G20.10)') QTEMP(1:3,J2) 
2300:                ENDDO 
2301:             ENDIF 
2302: 2270: 
2303: #ifdef MPI2271: #ifdef MPI
2304: !jdf43>Wait for all runs to get to 37 then collate statistics2272: !jdf43>Wait for all runs to get to 37 then collate statistics
2305:       CALL BHPT_IO(TIME-TSTART,NQTOT-1,NPCALL-1,NTOT,NEACCEPT)2273:       CALL BHPT_IO(TIME-TSTART,NQTOT-1,NPCALL-1,NTOT,NEACCEPT)
2306: 2274: 
2307:       WRITE(MYUNIT,10) NSUCCESST(JP)*1.0D0/MAX(1.0D0,1.0D0*(NSUCCESST(JP)+NFAILT(JP))),STEP(JP),ASTEP(JP),TEMP(JP)2275:       WRITE(MYUNIT,10) NSUCCESST(JP)*1.0D0/MAX(1.0D0,1.0D0*(NSUCCESST(JP)+NFAILT(JP))),STEP(JP),ASTEP(JP),TEMP(JP)
2308: 10    FORMAT('Acceptance ratio for run=',F12.5,' Step=',F12.5,' Angular step factor=',F12.5,' T=',F12.5)2276: 10    FORMAT('Acceptance ratio for run=',F12.5,' Step=',F12.5,' Angular step factor=',F12.5,' T=',F12.5)
2309:       IF (RIGID) WRITE(MYUNIT,'(A,F12.5)') 'Rigid body orientational step=',OSTEP(JP)2277:       IF (RIGID) WRITE(MYUNIT,'(A,F12.5)') 'Rigid body orientational step=',OSTEP(JP)
2310: 2278: 
2311: #else2279: #else


r31259/minpermdist.f90 2016-10-10 11:30:11.205225441 +0100 r31258/minpermdist.f90 2016-10-10 11:30:14.349267211 +0100
 49: !  where +/- is given by the value of INVERT. 49: !  where +/- is given by the value of INVERT.
 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the
 51: !  centre of coordinates of COORDSA will be the same as for COORDSB. 51: !  centre of coordinates of COORDSA will be the same as for COORDSB.
 52: ! 52: !
 53: !     ---------------------------------------------------------------------------------------------- 53: !     ----------------------------------------------------------------------------------------------
 54: ! jdf43>        Modified for generalised angle-axis 30/01/12 54: ! jdf43>        Modified for generalised angle-axis 30/01/12
 55: !     ---------------------------------------------------------------------------------------------- 55: !     ----------------------------------------------------------------------------------------------
 56:  56: 
 57: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST) 57: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 58: USE COMMONS,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, CHRMMT, MYUNIT, STOCKT, NFREEZE, & 58: USE COMMONS,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, CHRMMT, MYUNIT, STOCKT, NFREEZE, &
 59:   & AMBERT, CSMT, PERMDIST, PULLT, EFIELDT, OHCELLT, NTSITES, GEOMDIFFTOL, QCIPERMCHECK, QCIAMBERT, & 59:   & AMBERT, CSMT, PERMDIST, PULLT, EFIELDT, OHCELLT, NTSITES, GEOMDIFFTOL, QCIPERMCHECK, &
 60:   & PERMOPT, PERMINVOPT, NOINVERSION, BESTPERM, BESTINVERT, GTHOMSONT, LOCALPERMDIST,  LPERMDIST, MKTRAPT, AMBER12T 60:   & PERMOPT, PERMINVOPT, NOINVERSION, BESTPERM, BESTINVERT, GTHOMSONT, LOCALPERMDIST,  LPERMDIST, MKTRAPT
 61: USE PORFUNCS 61: USE PORFUNCS
 62: USE GENRIGID 62: USE GENRIGID
 63: IMPLICIT NONE 63: IMPLICIT NONE
 64:  64: 
 65: INTEGER, PARAMETER :: MAXIMUMTRIES=20 65: INTEGER, PARAMETER :: MAXIMUMTRIES=20
 66: INTEGER NATOMS, NPERM, PATOMS, NTRIES, ISTAT, OPNUM, NCHOOSEB1, NCHOOSEB2, NORBITB1, NORBITB2, I 66: INTEGER NATOMS, NPERM, PATOMS, NTRIES, ISTAT, OPNUM, NCHOOSEB1, NCHOOSEB2, NORBITB1, NORBITB2, I
 67: INTEGER J3, INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1 67: INTEGER J3, INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1
 68: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) 68: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS)
 69: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS) 69: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS)
 70: DOUBLE PRECISION CMXA, CMXB, CMXC, PDISTANCE, CMX, CMY, CMZ 70: DOUBLE PRECISION CMXA, CMXB, CMXC, PDISTANCE, CMX, CMY, CMZ
397: !     DUMMYA(1:3*NATOMS)=DUMMY(1:3*NATOMS)397: !     DUMMYA(1:3*NATOMS)=DUMMY(1:3*NATOMS)
398: !     CALL MYORIENT(DUMMYB,DUMMY,NORBITB1,NCHOOSEB1,NORBITB2,NCHOOSEB2,NATOMS,DEBUG,ROTB,ROTINVB,STOCKT)398: !     CALL MYORIENT(DUMMYB,DUMMY,NORBITB1,NCHOOSEB1,NORBITB2,NCHOOSEB2,NATOMS,DEBUG,ROTB,ROTINVB,STOCKT)
399: !     DUMMYB(1:3*NATOMS)=DUMMY(1:3*NATOMS)399: !     DUMMYB(1:3*NATOMS)=DUMMY(1:3*NATOMS)
400: !     DISTANCE=0.0D0400: !     DISTANCE=0.0D0
401: !     DO J1=1,3*NATOMS401: !     DO J1=1,3*NATOMS
402: !        DISTANCE=DISTANCE+(DUMMYA(J1)-DUMMYB(J1))**2402: !        DISTANCE=DISTANCE+(DUMMYA(J1)-DUMMYB(J1))**2
403: !     ENDDO403: !     ENDDO
404: !  IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'minpermdist> B distance=',SQRT(DISTANCE)404: !  IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'minpermdist> B distance=',SQRT(DISTANCE)
405: 405: 
406: ENDIF406: ENDIF
407:    IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A)') & 
408:   &       ' minpermdist> after initial call to MYORIENT distance=',SQRT(DISTANCE), ' for ',NATOMS,' atoms' 
409:    IF (DEBUG) WRITE(MYUNIT,'(A,6I8)') ' minpermdist> size of orbits, selected atoms, invert: ', & 
410:   &       NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,INVERT 
411:  
412: IF (DEBUG) WRITE(MYUNIT,'(A,4I8)') 'minpermdist> size of orbits and selected atoms: ',NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2407: IF (DEBUG) WRITE(MYUNIT,'(A,4I8)') 'minpermdist> size of orbits and selected atoms: ',NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2
413: 408: 
414: ! CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)409: ! CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
415: ! PRINT '(2(A,F25.15))',' New initial energy=',ENERGY,' RMS=',RMS410: ! PRINT '(2(A,F25.15))',' New initial energy=',ENERGY,' RMS=',RMS
416: ! PRINT '(2(A,F25.15))',' for coordinates:'411: ! PRINT '(2(A,F25.15))',' for coordinates:'
417: ! PRINT '(3F25.15)',DUMMYA(1:3*NATOMS)412: ! PRINT '(3F25.15)',DUMMYA(1:3*NATOMS)
418: !413: !
419: !  Bipartite matching routine for permutations. Coordinates in DUMMYB do not change414: !  Bipartite matching routine for permutations. Coordinates in DUMMYB do not change
420: !  but the coordinates in DUMMYA do. DISTANCE is the distance in this case.415: !  but the coordinates in DUMMYA do. DISTANCE is the distance in this case.
421: !  We return to label 10 after every round of permutational/orientational alignment416: !  We return to label 10 after every round of permutational/orientational alignment
624: IF (NCHOOSEB1.LT.NORBITB1) GOTO 66619: IF (NCHOOSEB1.LT.NORBITB1) GOTO 66
625: 620: 
626: !621: !
627: !  Now try the enantiomer (or xz reflected structure for TWOD.OR.PULLT.OR.EFIELDT).622: !  Now try the enantiomer (or xz reflected structure for TWOD.OR.PULLT.OR.EFIELDT).
628: !623: !
629: IF ((NCHOOSE2.EQ.NORBIT2).AND.(NCHOOSE1.EQ.NORBIT1).AND.(INVERT.EQ.1)) THEN624: IF ((NCHOOSE2.EQ.NORBIT2).AND.(NCHOOSE1.EQ.NORBIT1).AND.(INVERT.EQ.1)) THEN
630: !625: !
631: ! don't try inversion for bulk or charmm or amber or frozen atoms or CSM or PERMOPT626: ! don't try inversion for bulk or charmm or amber or frozen atoms or CSM or PERMOPT
632: ! Need keyword PERMINVOPT to do permutation-inversions.627: ! Need keyword PERMINVOPT to do permutation-inversions.
633: !628: !
634:    IF (NOINVERSION.OR.BULKT.OR.CHRMMT.OR.AMBERT.OR.AMBER12T.OR.QCIAMBERT.OR.(NFREEZE.GT.0).OR.CSMT.OR.PERMOPT) GOTO 50 629:    IF (NOINVERSION.OR.BULKT.OR.CHRMMT.OR.AMBERT.OR.(NFREEZE.GT.0).OR.CSMT.OR.PERMOPT) GOTO 50 
635:    IF (DEBUG) WRITE(MYUNIT,'(A)') 'minpermdist> inverting geometry for comparison with target'630:    IF (DEBUG) WRITE(MYUNIT,'(A)') 'minpermdist> inverting geometry for comparison with target'
636:    INVERT=-1631:    INVERT=-1
637:    GOTO 60632:    GOTO 60
638: ENDIF633: ENDIF
639: 634: 
640: 50 DISTANCE=DBEST635: 50 DISTANCE=DBEST
641: !636: !
642: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.637: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.
643: !  Rotate XBEST by ROTINVB to put in best correspondence with COORDSB, undoing the reorientation to DUMMYB from MYORIENT. 638: !  Rotate XBEST by ROTINVB to put in best correspondence with COORDSB, undoing the reorientation to DUMMYB from MYORIENT. 
644: !  We should get the same result for ROTINVB * RMATBEST * (COORDSA-CMA) 639: !  We should get the same result for ROTINVB * RMATBEST * (COORDSA-CMA) 


r31259/MLP3.f90 2016-10-10 11:30:08.829193874 +0100 r31258/MLP3.f90 2016-10-10 11:30:12.017236227 +0100
 17: ! w^1_{ij} at MLPHIDDEN*MLPIN + (i-1)*MLPHIDDEN+j 17: ! w^1_{ij} at MLPHIDDEN*MLPIN + (i-1)*MLPHIDDEN+j
 18: !   up to MLPHIDDEN*MLPIN + MLPOUT*MLPHIDDEN 18: !   up to MLPHIDDEN*MLPIN + MLPOUT*MLPHIDDEN
 19: ! 19: !
 20:  20: 
 21: MLPOFFSET=MLPHIDDEN*MLPIN 21: MLPOFFSET=MLPHIDDEN*MLPIN
 22: ENERGY=0.0D0 22: ENERGY=0.0D0
 23: V(1:NMLP)=0.0D0 23: V(1:NMLP)=0.0D0
 24: ! HESS(1:NMLP,1:NMLP)=0.0D0 24: ! HESS(1:NMLP,1:NMLP)=0.0D0
 25: DO J1=1,MLPDATA 25: DO J1=1,MLPDATA
 26:    MLPOUTJ1=MLPOUTCOME(J1) 26:    MLPOUTJ1=MLPOUTCOME(J1)
 27: !  WRITE(MYUNIT,'(A,3I10)') 'J1,MLPDATA,MLPOUTCOME(J1)=',J1,MLPDATA,MLPOUTCOME(J1) 
 28:    DO J2=1,MLPHIDDEN 27:    DO J2=1,MLPHIDDEN
 29:       DUMMY1=0.0D0 28:       DUMMY1=0.0D0
 30:       DO J3=1,MLPIN 29:       DO J3=1,MLPIN
 31:          DUMMY1=DUMMY1+X((J2-1)*MLPIN+J3)*MLPDAT(J1,J3) 30:          DUMMY1=DUMMY1+X((J2-1)*MLPIN+J3)*MLPDAT(J1,J3)
 32:       ENDDO 31:       ENDDO
 33:       TANHSUM(J2)=TANH(DUMMY1)  32:       TANHSUM(J2)=TANH(DUMMY1) 
 34:       DYW1G(J2)=TANHSUM(J2) 33:       DYW1G(J2)=TANHSUM(J2)
 35:       SECH2(J2)=1.0D0/COSH(DUMMY1)**2  34:       SECH2(J2)=1.0D0/COSH(DUMMY1)**2 
 36:    ENDDO 35:    ENDDO
 37:    DUMMY3=0.0D0 36:    DUMMY3=0.0D0
 41:          DO J3=1,MLPIN 40:          DO J3=1,MLPIN
 42:             DYW2G(J4,J2,J3)=X( MLPOFFSET + (J4-1)*MLPHIDDEN + J2 ) * MLPDAT(J1,J3)*SECH2(J2) 41:             DYW2G(J4,J2,J3)=X( MLPOFFSET + (J4-1)*MLPHIDDEN + J2 ) * MLPDAT(J1,J3)*SECH2(J2)
 43:          ENDDO 42:          ENDDO
 44:          DUMMY2=DUMMY2+X(MLPOFFSET+(J4-1)*MLPHIDDEN+J2)*TANHSUM(J2) 43:          DUMMY2=DUMMY2+X(MLPOFFSET+(J4-1)*MLPHIDDEN+J2)*TANHSUM(J2)
 45:       ENDDO 44:       ENDDO
 46:       Y(J4)=DUMMY2 45:       Y(J4)=DUMMY2
 47:       DUMMY3=DUMMY3+EXP(DUMMY2) 46:       DUMMY3=DUMMY3+EXP(DUMMY2)
 48:    ENDDO   47:    ENDDO  
 49:    DO J4=1,MLPOUT 48:    DO J4=1,MLPOUT
 50:       PROB(J4)=EXP(Y(J4))/DUMMY3 49:       PROB(J4)=EXP(Y(J4))/DUMMY3
 51: !     WRITE(MYUNIT,'(A,8G15.5)') 'J4,Y(J4),DUMMY3,PROB(J4)=',J4,Y(J4),DUMMY3,PROB(J4) 
 52:    ENDDO 50:    ENDDO
 53: !  WRITE(MYUNIT,'(A,I10)') 'MLPOUTJ1=',MLPOUTJ1 
 54:    PMLPOUTJ1=PROB(MLPOUTJ1) 51:    PMLPOUTJ1=PROB(MLPOUTJ1)
 55: !  WRITE(MYUNIT,'(A,G20.10)') 'PMLPOUTJ1',PMLPOUTJ1 52: !  IF (DEBUG) THEN
 56:    IF (DEBUG) THEN 53: !     WRITE(MYUNIT,'(A,I8,A)') 'MLP3> data point ',J1,' outputs and probabilities:'
 57:       WRITE(MYUNIT,'(A,I8,A)') 'MLP3> data point ',J1,' outputs and probabilities:' 54: !     WRITE(MYUNIT,'(8G15.5)') Y(1:MLPOUT),PROB(1:MLPOUT)
 58:       WRITE(MYUNIT,'(8G15.5)') Y(1:MLPOUT),PROB(1:MLPOUT) 55: !  ENDIF
 59:    ENDIF 
 60:    ENERGY=ENERGY-LOG(PMLPOUTJ1) 56:    ENERGY=ENERGY-LOG(PMLPOUTJ1)
 61:    IF (GTEST) THEN 57:    IF (GTEST) THEN
 62: ! 58: !
 63: ! We only need the probability derivative for the probability corresponding to the correct outcome for this data point 59: ! We only need the probability derivative for the probability corresponding to the correct outcome for this data point
 64: ! 60: !
 65:       DPCW1BG(1:MLPOUT,1:MLPHIDDEN)=0.0D0 61:       DPCW1BG(1:MLPOUT,1:MLPHIDDEN)=0.0D0
 66:       DO J2=1,MLPHIDDEN 62:       DO J2=1,MLPHIDDEN
 67:          DO J4=1,MLPOUT 63:          DO J4=1,MLPOUT
 68:             DPCW1BG(J4,J2)=DPCW1BG(J4,J2)-PMLPOUTJ1*PROB(J4)*DYW1G(J2) 64:             DPCW1BG(J4,J2)=DPCW1BG(J4,J2)-PMLPOUTJ1*PROB(J4)*DYW1G(J2)
 69:          ENDDO 65:          ENDDO


r31259/quench.F 2016-10-10 11:30:11.457228784 +0100 r31258/quench.F 2016-10-10 11:30:14.621270807 +0100
303: !           PEPSILON1(:)=PEPSILON1(:)/10000.0D0303: !           PEPSILON1(:)=PEPSILON1(:)/10000.0D0
304: !            CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)304: !            CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
305: !           WRITE(MYUNIT,*) 'second iteration: increasing epsilon_rep values by a factor of 100' 305: !           WRITE(MYUNIT,*) 'second iteration: increasing epsilon_rep values by a factor of 100' 
306: !           PEPSILON1(:)=PEPSILON1(:)*100.0D0306: !           PEPSILON1(:)=PEPSILON1(:)*100.0D0
307: !            CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)307: !            CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
308: !            WRITE(MYUNIT,*) 'third iteration: increasing epsilon_rep values by a factor of 100' 308: !            WRITE(MYUNIT,*) 'third iteration: increasing epsilon_rep values by a factor of 100' 
309: !           PEPSILON1(:)=PEPSILON1(:)*100.0D0309: !           PEPSILON1(:)=PEPSILON1(:)*100.0D0
310: !            CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)310: !            CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
311: !          END IF311: !          END IF
312:          ELSEIF ((.NOT.QCIPOTT).AND.INTCONSTRAINTT) THEN312:          ELSEIF ((.NOT.QCIPOTT).AND.INTCONSTRAINTT) THEN
 313:             WRITE(MYUNIT,'(A,2L5)') 'QCIPOTT,INTCONSTRAINTT=',QCIPOTT,INTCONSTRAINTT
313:             CALL CPU_TIME(T_LBFGS_START)314:             CALL CPU_TIME(T_LBFGS_START)
314:             CALL INTLBFGS(P,FINISH)315:             CALL INTLBFGS(P,FINISH)
315:             CALL CPU_TIME(T_LBFGS_END)316:             CALL CPU_TIME(T_LBFGS_END)
316:             STOP317:             STOP
317:          ELSE318:          ELSE
318: ! Added timing for LBFGS call319: ! Added timing for LBFGS call
319:             CALL CPU_TIME(T_LBFGS_START)320:             CALL CPU_TIME(T_LBFGS_START)
320:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)321:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
321:             CALL CPU_TIME(T_LBFGS_END)322:             CALL CPU_TIME(T_LBFGS_END)
322:             IF (FEBHT.AND.DEBUG) WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START323:             IF (FEBHT.AND.DEBUG) WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0