hdiff output

r29153/commons.f90 2015-11-17 23:33:35.336978654 +0000 r29152/commons.f90 2015-11-17 23:33:37.813011859 +0000
 71:      &                 RHOCH10, RHOC20H, RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ, SYMFCTR, & 71:      &                 RHOCH10, RHOC20H, RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ, SYMFCTR, &
 72:      &                 GBANISOTROPYR,GBWELLDEPTHR,PARAMa1,PARAMb1,PARAMc1,PARAMa2,PARAMB2,PARAMc2,PSIGMA0(2),PEPSILON0,& 72:      &                 GBANISOTROPYR,GBWELLDEPTHR,PARAMa1,PARAMb1,PARAMc1,PARAMa2,PARAMB2,PARAMc2,PSIGMA0(2),PEPSILON0,&
 73:      &                 PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2),PYA11(3),PYA21(3),PYA12(3),PYA22(3), & 73:      &                 PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2),PYA11(3),PYA21(3),PYA12(3),PYA22(3), &
 74:      &                 PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, PYCFTHRESH, LJSITECOORDS(3), & 74:      &                 PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, PYCFTHRESH, LJSITECOORDS(3), &
 75:      &                 MSTART,MFINISH,MBSTART1,MBFINISH1,MBSTART2,MBFINISH2,MBHEIGHT1,MBHEIGHT2,ME1,ME2,ME3, & 75:      &                 MSTART,MFINISH,MBSTART1,MBFINISH1,MBSTART2,MBFINISH2,MBHEIGHT1,MBHEIGHT2,ME1,ME2,ME3, &
 76:      &                 BSPTQMAX, BSPTQMIN, PFORCE, CSMNORM, CSMGUIDENORM, CSMEPS, PERCCUT, & 76:      &                 BSPTQMAX, BSPTQMIN, PFORCE, CSMNORM, CSMGUIDENORM, CSMEPS, PERCCUT, &
 77:      &                 LOWESTE, PERTSTEP, GCPLUS, & 77:      &                 LOWESTE, PERTSTEP, GCPLUS, &
 78:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, & 78:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, &
 79:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, & 79:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, &
 80:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 80:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 81:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL 81:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF
 82:  82: 
 83:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, & 83:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, &
 84:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, & 84:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, &
 85:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, JMT, LJCOULT, SETCENT, & 85:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, JMT, LJCOULT, SETCENT, &
 86:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, & 86:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, &
 87:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, & 87:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, &
 88:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, & 88:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, &
 89:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, & 89:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, &
 90:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, & 90:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, &
 91:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, &  91:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, & 
104:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF,&104:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF,&
105:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&105:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
106:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &106:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &
107:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &107:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
108:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &108:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
109:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &109:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
110:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &110:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
111:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 111:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
112:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, &112:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, &
113:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &113:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
114:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET114:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT
115: !115: !
116:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 116:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
117:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)117:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)
118:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 118:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 
119: 119: 
120: ! csw34> FREEZEGROUP variables120: ! csw34> FREEZEGROUP variables
121: !121: !
122:       INTEGER :: GROUPCENTRE122:       INTEGER :: GROUPCENTRE
123:       DOUBLE PRECISION :: GROUPRADIUS123:       DOUBLE PRECISION :: GROUPRADIUS
124:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE124:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE


r29153/congrad.f90 2015-11-17 23:33:35.528981232 +0000 r29152/congrad.f90 2015-11-17 23:33:38.001014387 +0000
 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 CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 19: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 20: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, & 20: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
 21:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, & 21:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
 22:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, & 22:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, &
 23:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, & 23:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, &
 24:   &  NATOMS, DEBUG, MYUNIT, INTMINFAC, FREEZENODEST, INTSPRINGACTIVET 24:   &  NATOMS, DEBUG, MYUNIT
 25: USE PORFUNCS 25: USE PORFUNCS
 26: IMPLICIT NONE 26: IMPLICIT NONE
 27:             27:            
 28: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2 28: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT
 29: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS 29: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
 30: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1 30: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3) 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
 32: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI 32: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
 33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2) 33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
 34: LOGICAL NOINT 34: LOGICAL NOINT
 35: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL 35: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL
 36: LOGICAL IMGFREEZE(INTIMAGE) 36: LOGICAL IMGFREEZE(INTIMAGE)
 37: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3) 37: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3)
 38:  38: 
 39: EEE(1:INTIMAGE+2)=0.0D0 39: EEE(1:INTIMAGE+2)=0.0D0
 40: CONE(1:INTIMAGE+2)=0.0D0 40: CONE(1:INTIMAGE+2)=0.0D0
 41: REPE(1:INTIMAGE+2)=0.0D0 41: REPE(1:INTIMAGE+2)=0.0D0
 42: REPEINT(1:INTIMAGE+2)=0.0D0 42: REPEINT(1:INTIMAGE+2)=0.0D0
 43: NREPINT(1:INTIMAGE+2)=0 43: NREPINT(1:INTIMAGE+2)=0
 44: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 44: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0
 45: ECON=0.0D0; EREP=0.0D0 45: ECON=0.0D0; EREP=0.0D0
 46: NINTMIN=0 
 47: NINTMIN2=0 
 48:  46: 
 49: ! 47: !
 50: !  Constraint energy and forces. 48: !  Constraint energy and forces.
 51: ! 49: !
 52: DO J2=1,NCONSTRAINT 50: DO J2=1,NCONSTRAINT
 53:    IF (.NOT.CONACTIVE(J2)) CYCLE 51:    IF (.NOT.CONACTIVE(J2)) CYCLE
 54:       CCLOCAL=CONCUTLOCAL(J2) 52:       CCLOCAL=CONCUTLOCAL(J2)
 55:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 53:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
 56:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 54:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
 57: ! 55: !
 84: ENDDO 82: ENDDO
 85:  83: 
 86: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes 84: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
 87: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes 85: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
 88:  86: 
 89: ! INTCONST=INTCONSTRAINREPCUT**13 87: ! INTCONST=INTCONSTRAINREPCUT**13
 90:  88: 
 91: DO J2=1,NNREPULSIVE 89: DO J2=1,NNREPULSIVE
 92: !  INTCONST=NREPCUT(J2)**13 90: !  INTCONST=NREPCUT(J2)**13
 93:    INTCONST=NREPCUT(J2)**3 91:    INTCONST=NREPCUT(J2)**3
 94:    DO J1=2,INTIMAGE+2 92: !  DO J1=2,INTIMAGE+2
 95: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images 93:    DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images
 96:       IF (FREEZENODEST) THEN 
 97:          IF (J1.EQ.2) THEN 
 98:             IF (IMGFREEZE(1)) CYCLE 
 99:          ELSE IF (J1.EQ.INTIMAGE+2) THEN 
100:             IF (IMGFREEZE(INTIMAGE)) CYCLE 
101:          ELSE 
102:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE 
103:          ENDIF 
104:       ENDIF 
105: !  IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN 
106: !     PRINT '(A,I6,A,2I6)',' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2) 
107: !     STOP 
108: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad2> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6) 
109:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1) 94:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
110:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1) 95:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
111:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3) 96:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
112:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3) 97:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
113:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 98:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)
114:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1 99:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1
115: !        D12=D2**12100: !        D12=D2**12
116:          D12=D2**2101:          D12=D2**2
117: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)102: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
118:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)103:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
142:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)127:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
143:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)128:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
144: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN129: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
145: !        PRINT '(A,I6,A,2I6)','B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)130: !        PRINT '(A,I6,A,2I6)','B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
146: !        PRINT '(A,I6)','image number ',J1131: !        PRINT '(A,I6)','image number ',J1
147: !        PRINT '(A,6F15.10)','R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ132: !        PRINT '(A,6F15.10)','R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
148: !        PRINT '(A,6F15.10)','R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ133: !        PRINT '(A,6F15.10)','R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
149: !     ENDIF134: !     ENDIF
150:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &135:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
151:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))136:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
152:       IF (.NOT.NOINT) THEN 
153: !        WRITE(MYUNIT,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), & 
154: ! &                     ' distance,cutoff=',DINT,NREPCUT(J2) 
155:          NINTMIN=NINTMIN+1 
156:       ENDIF 
157:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN137:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
158:          NINTMIN2=NINTMIN2+1 
159: !        D12=DSQI**6138: !        D12=DSQI**6
160:          D12=DSQI139:          D12=DSQI
161: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)140: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
162:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)141:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
163:          IF (J1.EQ.2) THEN142:          IF (J1.EQ.2) THEN
164:             EEE(J1)=EEE(J1)+DUMMY143:             EEE(J1)=EEE(J1)+DUMMY
165:             REPEINT(J1)=REPEINT(J1)+DUMMY144:             REPEINT(J1)=REPEINT(J1)+DUMMY
166:             NREPINT(J1)=NREPINT(J1)+1145:             NREPINT(J1)=NREPINT(J1)+1
167:          ELSE IF (J1.LT.INTIMAGE+2) THEN146:          ELSE IF (J1.LT.INTIMAGE+2) THEN
168:             EEE(J1)=EEE(J1)+DUMMY/2.0D0147:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
169:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0148:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
170:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0149:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0
171:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0150:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
172:             NREPINT(J1)=NREPINT(J1)+1151:             NREPINT(J1)=NREPINT(J1)+1
173:             NREPINT(J1-1)=NREPINT(J1-1)+1152:             NREPINT(J1-1)=NREPINT(J1-1)+1
174:          ELSE IF (J1.EQ.INTIMAGE+2) THEN153:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
175:             EEE(J1-1)=EEE(J1-1)+DUMMY154:             EEE(J1-1)=EEE(J1-1)+DUMMY
176:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY155:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
177:             NREPINT(J1-1)=NREPINT(J1-1)+1156:             NREPINT(J1-1)=NREPINT(J1-1)+1
178:          ENDIF157:          ENDIF
179:          EREP=EREP+DUMMY158:          EREP=EREP+DUMMY
180: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)159: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
181:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)160:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
182:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)161:          REPGRAD(1:3)=DUMMY*G1INT(1:3)
183: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &162: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
184: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)163: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
185: !164: !
186: ! Gradient contributions for image J1-1165: ! Gradient contributions for image J1-1
187: !166: !
188:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)167:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
189:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)168:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
190:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)169:          REPGRAD(1:3)=DUMMY*G2INT(1:3)
191: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &170: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
192: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)171: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
193: !172: !
194: ! Gradient contributions for image J1173: ! Gradient contributions for image J1
195: !174: !
196:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)175:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
197:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)176:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
198:       ENDIF177:       ENDIF
199:    ENDDO178:    ENDDO
200: ENDDO179: ENDDO
205: ESPRING=0.0D0184: ESPRING=0.0D0
206: IF (KINT.NE.0.0D0) THEN185: IF (KINT.NE.0.0D0) THEN
207:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1186:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1
208:       NI1=(3*NATOMS)*(J1-1)187:       NI1=(3*NATOMS)*(J1-1)
209:       NI2=(3*NATOMS)*J1188:       NI2=(3*NATOMS)*J1
210: !189: !
211: !  Edge between J1 and J1+1190: !  Edge between J1 and J1+1
212: !191: !
213:       DPLUS=0.0D0192:       DPLUS=0.0D0
214:       DO J2=1,NATOMS193:       DO J2=1,NATOMS
215:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 194:          DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &
216:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &195:   &                 +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &
217:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &196:   &                 +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2
218:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2 
219:          ENDIF 
220:       ENDDO197:       ENDDO
221:       DPLUS=SQRT(DPLUS)198:       DPLUS=SQRT(DPLUS)
222:       IF (DPLUS.GT.IMSEPMAX) THEN199:       IF (DPLUS.GT.IMSEPMAX) THEN
223: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2200: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2
224:          DUMMY=KINT*0.5D0*DPLUS**2201:          DUMMY=KINT*0.5D0*DPLUS**2
225:          ESPRING=ESPRING+DUMMY202:          ESPRING=ESPRING+DUMMY
226: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS203: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS
227:          DUMMY=KINT204:          DUMMY=KINT
228:          DO J2=1,NATOMS205:          DO J2=1,NATOMS
229:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 206:             SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))
230:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))207:             GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)
231:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)208:             GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3)
232:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3) 
233:             ENDIF 
234:          ENDDO209:          ENDDO
235:       ENDIF210:       ENDIF
236:    ENDDO211:    ENDDO
237: ENDIF212: ENDIF
238: !213: !
239: ! Set gradients on frozen atoms to zero.214: ! Set gradients on frozen atoms to zero.
240: !215: !
241: IF (FREEZE) THEN216: IF (FREEZE) THEN
242:    DO J1=2,INTIMAGE+1  217:    DO J1=2,INTIMAGE+1  
243:       DO J2=1,NATOMS218:       DO J2=1,NATOMS
282: ENDIF257: ENDIF
283: !258: !
284: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,259: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,
285: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)260: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)
286: !261: !
287: ETOTAL=0.0D0262: ETOTAL=0.0D0
288: MAXINT=-1.0D100263: MAXINT=-1.0D100
289: MININT=1.0D100264: MININT=1.0D100
290: DO J1=2,INTIMAGE+1265: DO J1=2,INTIMAGE+1
291:    ETOTAL=ETOTAL+EEE(J1)266:    ETOTAL=ETOTAL+EEE(J1)
292: !  WRITE(MYUNIT, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)267:    WRITE(MYUNIT, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)
293:    IF (REPEINT(J1).LT.MININT) THEN268:    IF (REPEINT(J1).LT.MININT) THEN
294:       MININT=REPEINT(J1)269:       MININT=REPEINT(J1)
295:       NMININT=J1270:       NMININT=J1
296:    ENDIF271:    ENDIF
297:    IF (REPE(J1).GT.MAXINT) THEN272:    IF (REPE(J1).GT.MAXINT) THEN
298:       MAXINT=REPE(J1)273:       MAXINT=REPE(J1)
299:       NMAXINT=J1274:       NMAXINT=J1
300:    ENDIF275:    ENDIF
301: ENDDO276: ENDDO
302: IF (DEBUG) WRITE(MYUNIT, '(A,G20.10,A,2I6)') 'congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT277: IF (DEBUG) WRITE(MYUNIT, '(A,G20.10,A,2I6)') ' congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT
303: IF (DEBUG) WRITE(MYUNIT, '(A,G20.10,A,2I6)') 'congrad> smallest internal energy=',MININT,' for image ',NMININT278: IF (DEBUG) WRITE(MYUNIT, '(A,G20.10,A,2I6)') ' congrad> smallest internal energy=',MININT,' for image ',NMININT
304: IF (DEBUG) WRITE(MYUNIT, '(A,2I6)') 'congrad> number of internal minima=',NINTMIN,NINTMIN2 
305: 279: 
306: END SUBROUTINE CONGRAD280: END SUBROUTINE CONGRAD
307: 281: 
308: SUBROUTINE MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &282: SUBROUTINE MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
309:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,DEBUG)283:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,DEBUG)
310: IMPLICIT NONE284: IMPLICIT NONE
311: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT285: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT
312: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)286: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)
313: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq287: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq
314: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY288: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY
500: 474: 
501: !475: !
502: ! This version of congrad tests for an internal minimum in the476: ! This version of congrad tests for an internal minimum in the
503: ! constraint distances as well as the repulsions.477: ! constraint distances as well as the repulsions.
504: !478: !
505: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)479: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
506: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &480: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
507:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &481:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
508:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &482:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &
509:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &483:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &
510:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, NATOMS, DEBUG, MYUNIT, INTMINFAC, INTSPRINGACTIVET 484:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, NATOMS, DEBUG, MYUNIT
511: IMPLICIT NONE485: IMPLICIT NONE
512:            486:            
513: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)487: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)
514: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS488: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
515: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1489: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
516: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)490: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
517: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI491: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
518: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)492: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)
519: LOGICAL NOINT, LPRINT493: LOGICAL NOINT, LPRINT
520: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)494: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)
586:       CALL MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &560:       CALL MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
587:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,.FALSE.)561:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,.FALSE.)
588: !562: !
589: ! Need to include both D2 and D1 contributions if they are both outside tolerance.563: ! Need to include both D2 and D1 contributions if they are both outside tolerance.
590: ! Otherwise we get discontinuities if they are very close and swap over.564: ! Otherwise we get discontinuities if they are very close and swap over.
591: !565: !
592: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)566: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)
593: !567: !
594: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.568: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.
595: !569: !
596:       IF (LPRINT) WRITE(MYUNIT, '(A,I6,5G15.5)') &570: !     IF (LPRINT) PRINT '(A,I6,5G15.5)', &
597:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL571: ! &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL
598:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 572:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 
599:          DUMMY=D2-CONDISTREFLOCAL(J2)  573:          DUMMY=D2-CONDISTREFLOCAL(J2)  
600:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)574:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
601:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)575:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
602:          EEE(J1)=EEE(J1)+DUMMY576:          EEE(J1)=EEE(J1)+DUMMY
603:          CONE(J1)=CONE(J1)+DUMMY577:          CONE(J1)=CONE(J1)+DUMMY
604:          ECON=ECON      +DUMMY578:          ECON=ECON      +DUMMY
605:          IF (LPRINT) WRITE(MYUNIT, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &579: !        IF (LPRINT) PRINT '(A,4I6,G15.5)','min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
606:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)580: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
607:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)581:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
608:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)582:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
609:       ENDIF583:       ENDIF
610: !584: !
611: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.585: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.
612: !586: !
613: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.587: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.
614: !588: !
615: !     IF (LPRINT) WRITE(MYUNIT, '(A,I6,5G15.5)') &589: !     IF (LPRINT) PRINT '(A,I6,5G15.5)', &
616: ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL590: ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL
617:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  591:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  
618:          DUMMY=D1-CONDISTREFLOCAL(J2)  592:          DUMMY=D1-CONDISTREFLOCAL(J2)  
619:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)593:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)
620:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)594:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
621:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN595:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
622:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &596:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
623:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY597:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
624:          ENDIF598:          ENDIF
625:          EEE(J1-1)=EEE(J1-1)+DUMMY599:          EEE(J1-1)=EEE(J1-1)+DUMMY
626:          CONE(J1-1)=CONE(J1-1)+DUMMY600:          CONE(J1-1)=CONE(J1-1)+DUMMY
627:          ECON=ECON      +DUMMY601:          ECON=ECON      +DUMMY
628:          IF (LPRINT) WRITE(MYUNIT, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &602: !        IF (LPRINT) PRINT '(A,4I6,G15.5)','max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
629:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)603: ! &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
630:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)604:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
631:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)605:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
632:       ENDIF606:       ENDIF
633:       IF ((.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN607:       IF ((.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN
634:          DUMMY=DINT-CONDISTREFLOCAL(J2)  608:          DUMMY=DINT-CONDISTREFLOCAL(J2)  
635:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)609:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)
636:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)610:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
637:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)611:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
638:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)612:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)
639:          DUMMY=INTMINFAC*INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)613:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
640:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN614:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
641:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &615:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
642:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY616:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
643:          ENDIF617:          ENDIF
644:          ECON=ECON+DUMMY618:          ECON=ECON+DUMMY
 619: !        IF ((J2.EQ.304).OR.(J2.EQ.309).OR.(J2.EQ.311)) THEN
 620: !           PRINT '(A,G20.10)','adding INT term DUMMY=',DUMMY
 621: !        ENDIF
645:          IF (J1.EQ.2) THEN622:          IF (J1.EQ.2) THEN
646:             EEE(J1)=EEE(J1)+DUMMY623:             EEE(J1)=EEE(J1)+DUMMY
647:             CONEINT(J1)=CONEINT(J1)+DUMMY624:             CONEINT(J1)=CONEINT(J1)+DUMMY
648:             NCONINT(J1)=NCONINT(J1)+1625:             NCONINT(J1)=NCONINT(J1)+1
649:          ELSE IF (J1.LT.INTIMAGE+2) THEN626:          ELSE IF (J1.LT.INTIMAGE+2) THEN
650:             EEE(J1)=EEE(J1)+DUMMY/2.0D0627:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
651:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0628:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
652:             CONEINT(J1)=CONEINT(J1)+DUMMY/2.0D0629:             CONEINT(J1)=CONEINT(J1)+DUMMY/2.0D0
653:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY/2.0D0630:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY/2.0D0
654:             NCONINT(J1)=NCONINT(J1)+1631:             NCONINT(J1)=NCONINT(J1)+1
655:             NCONINT(J1-1)=NCONINT(J1-1)+1632:             NCONINT(J1-1)=NCONINT(J1-1)+1
656:          ELSE IF (J1.EQ.INTIMAGE+2) THEN633:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
657:             EEE(J1-1)=EEE(J1-1)+DUMMY634:             EEE(J1-1)=EEE(J1-1)+DUMMY
658:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY635:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY
659:             NCONINT(J1-1)=NCONINT(J1-1)+1636:             NCONINT(J1-1)=NCONINT(J1-1)+1
660:          ENDIF637:          ENDIF
661: !        WRITE(MYUNIT, '(A,4I6,G15.5)') 'in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &638: !        PRINT '(A,4I6,G15.5)','in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
662: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)639: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
663:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)640:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
664:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)641:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
665:       ENDIF642:       ENDIF
666:    ENDDO643:    ENDDO
667: ENDDO644: ENDDO
668: 645: 
669: ! INTCONST=INTCONSTRAINREPCUT**13646: ! INTCONST=INTCONSTRAINREPCUT**13
670: 647: 
671: DO J2=1,NNREPULSIVE648: DO J2=1,NNREPULSIVE
672: !  INTCONST=NREPCUT(J2)**13649: !  INTCONST=NREPCUT(J2)**13
673:    INTCONST=NREPCUT(J2)**3650:    INTCONST=NREPCUT(J2)**3
674:    DO J1=2,INTIMAGE+2651:    DO J1=2,INTIMAGE+2
675:       IF (FREEZENODEST) THEN652:       If (FREEZENODEST) THEN
676:          IF (J1.EQ.2) THEN653:          IF (J1.EQ.2) THEN
677:             IF (IMGFREEZE(1)) CYCLE654:             IF (IMGFREEZE(1)) CYCLE
678:          ELSE IF (J1.EQ.INTIMAGE+2) THEN655:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
679:             IF (IMGFREEZE(INTIMAGE)) CYCLE656:             IF (IMGFREEZE(INTIMAGE)) CYCLE
680:          ELSE657:          ELSE
681:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE658:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
682:          ENDIF659:          ENDIF
683:       ENDIF660:       ENDIF
684: !     IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN661: !  IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
685: !        PRINT '(A,I6,A,2I6)',' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)662: !     PRINT '(A,I6,A,2I6)',' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
686: !        STOP663: !     STOP
687: !     ENDIF664: !  ENDIF
688: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad2> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6) 
689:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)665:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
690:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)666:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
691:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)667:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
692:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)668:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
693:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)669:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
694:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)670:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
695:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)671:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
696:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)672:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
697:       IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN673:       IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
698: !        WRITE(MYUNIT, '(A,I6,A,2I6)') 'A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)674:          WRITE(MYUNIT, '(A,I6,A,2I6)') 'A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
699: !        WRITE(MYUNIT, '(A,I6)') 'image number ',J1675:          WRITE(MYUNIT, '(A,I6)') 'image number ',J1
700: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ676:          WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
701: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ677:          WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
702:       ENDIF678:       ENDIF
703:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &679:       CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
704:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))680:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
705: !     IF ((NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.357)) THEN681: !     IF ((NREPI(J2).EQ.135).AND.(NREPJ(J2).EQ.192)) THEN
706: !        WRITE(MYUNIT, '(A,3G20.10)') ' congrad2> R1AX,R1AY,R1AZ=',R1AX,R1AY,R1AZ682: !        PRINT '(A,3G20.10)',' congrad2> R1AX,R1AY,R1AZ=',R1AX,R1AY,R1AZ
707: !        WRITE(MYUNIT, '(A,3G20.10)') ' congrad2> R1BX,R1BY,R1BZ=',R1BX,R1BY,R1BZ683: !        PRINT '(A,3G20.10)',' congrad2> R1BX,R1BY,R1BZ=',R1BX,R1BY,R1BZ
708: !        WRITE(MYUNIT, '(A,3G20.10)') ' congrad2> R2AX,R2AY,R2AZ=',R2AX,R2AY,R2AZ684: !        PRINT '(A,3G20.10)',' congrad2> R2AX,R2AY,R2AZ=',R2AX,R2AY,R2AZ
709: !        WRITE(MYUNIT, '(A,3G20.10)') ' congrad2> R2BX,R2BY,R2BZ=',R2BX,R2BY,R2BZ685: !        PRINT '(A,3G20.10)',' congrad2> R2BX,R2BY,R2BZ=',R2BX,R2BY,R2BZ
710: !        WRITE(MYUNIT, '(A,I6,A,2I6)') ' congrad2> J1=',J1,' edge between images: ',J1-1,J1686: !        PRINT '(A,I6,A,2I6)',' congrad2> J1=',J1,' edge between images: ',J1-1,J1
711: !        WRITE(MYUNIT, '(A,L5,3G20.10)') ' congrad2> NOINT,D2,D1,DINT=',NOINT,D2,D1,DINT687: !        PRINT '(A,L5,3G20.10)',' congrad2> NOINT,D2,D1,DINT=',NOINT,D2,D1,DINT
712: !     ENDIF688: !     ENDIF
713: !     IF (.NOT.NOINT.AND.(NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.357)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
714: ! &                                                 ' congrad2> im 83 357 D2,D1,DINT=',J1,D2,D1,DINT 
715: !     IF (.NOT.NOINT.AND.(NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.359)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
716: ! &                                                 ' congrad2> im 83 359 D2,D1,DINT=',J1,D2,D1,DINT 
717: !     IF (.NOT.NOINT.AND.(NREPI(J2).EQ.86).AND.(NREPJ(J2).EQ.357)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
718: ! &                                                 ' congrad2> im 86 357 D2,D1,DINT=',J1,D2,D1,DINT 
719: !     IF (.NOT.NOINT.AND.(NREPI(J2).EQ.86).AND.(NREPJ(J2).EQ.359)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
720: ! &                                                 ' congrad2> im 86 359 D2,D1,DINT=',J1,D2,D1,DINT 
721: !     IF (.NOT.NOINT.AND.(NREPJ(J2).EQ.83).AND.(NREPI(J2).EQ.357)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
722: ! &                                                 ' congrad2> im 83 357 D2,D1,DINT=',J1,D2,D1,DINT 
723: !     IF (.NOT.NOINT.AND.(NREPJ(J2).EQ.83).AND.(NREPI(J2).EQ.359)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
724: ! &                                                 ' congrad2> im 83 359 D2,D1,DINT=',J1,D2,D1,DINT 
725: !     IF (.NOT.NOINT.AND.(NREPJ(J2).EQ.86).AND.(NREPI(J2).EQ.357)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
726: ! &                                                 ' congrad2> im 86 357 D2,D1,DINT=',J1,D2,D1,DINT 
727: !     IF (.NOT.NOINT.AND.(NREPJ(J2).EQ.86).AND.(NREPI(J2).EQ.359)) WRITE(MYUNIT, '(A,I5,3G20.10)') & 
728: ! &                                                 ' congrad2> im 86 359 D2,D1,DINT=',J1,D2,D1,DINT 
729:       DUMMY=0.0D0 689:       DUMMY=0.0D0 
730: !690: !
731: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.691: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
732: !692: !
733: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1693: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
734:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1694:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
735: !        D12=DSQ2**6695: !        D12=DSQ2**6
736:          D12=DSQ2696:          D12=DSQ2
737: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)697: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
738: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)698: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
739:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)699:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
740: !        IF ((NREPJ(J2).EQ.83).AND.(NREPI(J2).EQ.357)) WRITE(MYUNIT, '(A,6G20.10)') & 
741: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY  
742: !        IF ((NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.357)) WRITE(MYUNIT, '(A,6G20.10)') & 
743: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY  
744:          EEE(J1)=EEE(J1)+DUMMY700:          EEE(J1)=EEE(J1)+DUMMY
745:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN701:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
746:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &702:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
747:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY703:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
748:          ENDIF704:          ENDIF
749:          REPE(J1)=REPE(J1)+DUMMY705:          REPE(J1)=REPE(J1)+DUMMY
750:          EREP=EREP+DUMMY706:          EREP=EREP+DUMMY
751: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)707: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
752:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)708:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
753:          REPGRAD(1:3)=DUMMY*G2(1:3)709:          REPGRAD(1:3)=DUMMY*G2(1:3)
754: !        WRITE(MYUNIT, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &710: !        PRINT '(A,4I6,G15.5)','min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
755: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)711: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
756:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)712:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
757:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)713:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
758:       ENDIF714:       ENDIF
759:       DUMMY=0.0D0715:       DUMMY=0.0D0
760: !716: !
761: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.717: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.
762: !718: !
763: !     IF ((D1.LT.INTCONSTRAINREPCUT).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1719: !     IF ((D1.LT.INTCONSTRAINREPCUT).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1
764:       IF ((D1.LT.NREPCUT(J2)).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1720:       IF ((D1.LT.NREPCUT(J2)).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1
770:          EEE(J1-1)=EEE(J1-1)+DUMMY726:          EEE(J1-1)=EEE(J1-1)+DUMMY
771:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN727:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
772:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'R2 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &728:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'R2 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
773:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY729:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
774:          ENDIF730:          ENDIF
775:          REPE(J1-1)=REPE(J1-1)+DUMMY731:          REPE(J1-1)=REPE(J1-1)+DUMMY
776:          EREP=EREP+DUMMY732:          EREP=EREP+DUMMY
777: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)733: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)
778:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)734:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)
779:          REPGRAD(1:3)=DUMMY*G1(1:3)735:          REPGRAD(1:3)=DUMMY*G1(1:3)
780: !        WRITE(MYUNIT, '(A,4I6,G15.5)') 'max J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &736: !        PRINT '(A,4I6,G15.5)','max J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
781: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)737: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
782:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)738:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
783:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)739:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
784:       ENDIF740:       ENDIF
785:       DUMMY=0.0D0741:       DUMMY=0.0D0
786: !     IF ((.NOT.NOINT).AND.(DINT.LT.INTCONSTRAINREPCUT)) THEN742: !     IF ((.NOT.NOINT).AND.(DINT.LT.INTCONSTRAINREPCUT)) THEN
787:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN743:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
788: !        D12=DSQI**6744: !        D12=DSQI**6
789:          D12=DSQI745:          D12=DSQI
790: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*INTCONSTRAINREPCUT)/INTCONST)746: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
791: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)747: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
792:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)748:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
793:          EREP=EREP+DUMMY749:          EREP=EREP+DUMMY
794: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &750: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &
795: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY751: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY
 752: !        IF (((NREPI(J2).EQ.143).AND.(NREPJ(J2).EQ.191)).OR.  &
 753: ! &          ((NREPJ(J2).EQ.143).AND.(NREPI(J2).EQ.191))) THEN
 754: !            PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &
 755: ! &                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY
796: !        ENDIF756: !        ENDIF
797:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN757:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
798:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &758:             WRITE(MYUNIT, '(A,2I6,2L5,G20.10)') 'R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
799:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY759:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
800:          ENDIF760:          ENDIF
801:          IF (J1.EQ.2) THEN761:          IF (J1.EQ.2) THEN
802:             EEE(J1)=EEE(J1)+DUMMY762:             EEE(J1)=EEE(J1)+DUMMY
803:             REPEINT(J1)=REPEINT(J1)+DUMMY763:             REPEINT(J1)=REPEINT(J1)+DUMMY
804:             NREPINT(J1)=NREPINT(J1)+1764:             NREPINT(J1)=NREPINT(J1)+1
805:          ELSE IF (J1.LT.INTIMAGE+2) THEN765:          ELSE IF (J1.LT.INTIMAGE+2) THEN
809:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0769:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
810:             NREPINT(J1)=NREPINT(J1)+1770:             NREPINT(J1)=NREPINT(J1)+1
811:             NREPINT(J1-1)=NREPINT(J1-1)+1771:             NREPINT(J1-1)=NREPINT(J1-1)+1
812:          ELSE IF (J1.EQ.INTIMAGE+2) THEN772:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
813:             EEE(J1-1)=EEE(J1-1)+DUMMY773:             EEE(J1-1)=EEE(J1-1)+DUMMY
814:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY774:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
815:             NREPINT(J1-1)=NREPINT(J1-1)+1775:             NREPINT(J1-1)=NREPINT(J1-1)+1
816:          ENDIF776:          ENDIF
817: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)777: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
818:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)778:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
819:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)779:          REPGRAD(1:3)=DUMMY*G1INT(1:3)
820: !        WRITE(MYUNIT, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &780: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
821: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)781: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
822:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)782:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
823:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)783:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
824:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)784:          REPGRAD(1:3)=DUMMY*G2INT(1:3)
825: !        WRITE(MYUNIT, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &785: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
826: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)786: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
827:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)787:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
828:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)788:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
829:       ENDIF789:       ENDIF
830: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad2> C J1,J2,GGG(1:6)=',J1,J2,GGG(1:6) 
831:    ENDDO790:    ENDDO
832: ENDDO791: ENDDO
833: !792: !
834: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise793: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
835: ! energy terms between images except for the end points.794: ! energy terms between images except for the end points.
836: !795: !
837: ESPRING=0.0D0796: ESPRING=0.0D0
838: IF (KINT.NE.0.0D0) THEN797: IF (KINT.NE.0.0D0) THEN
839:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1798:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1
840:       NI1=(3*NATOMS)*(J1-1)799:       NI1=(3*NATOMS)*(J1-1)
841:       NI2=(3*NATOMS)*J1800:       NI2=(3*NATOMS)*J1
842: !801: !
843: !  Edge between J1 and J1+1802: !  Edge between J1 and J1+1
844: !803: !
845:       DPLUS=0.0D0804:       DPLUS=0.0D0
846:       DO J2=1,NATOMS805:       DO J2=1,NATOMS
847:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 806:          DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &
848:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &807:   &                 +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &
849:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &808:   &                 +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2
850:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2 
851:          ENDIF 
852: !        WRITE(MYUNIT,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS 
853:       ENDDO809:       ENDDO
854:       DPLUS=SQRT(DPLUS)810:       DPLUS=SQRT(DPLUS)
855:       IF (DPLUS.GT.IMSEPMAX) THEN811:       IF (DPLUS.GT.IMSEPMAX) THEN
856:          DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2812:          DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2
857: !        DUMMY=KINT*0.5D0*DPLUS**2813: !        DUMMY=KINT*0.5D0*DPLUS**2
858:          ESPRING=ESPRING+DUMMY814:          ESPRING=ESPRING+DUMMY
859: !        WRITE(MYUNIT,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING 
860:          DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS815:          DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS
861: !        DUMMY=KINT816: !        DUMMY=KINT
862:          DO J2=1,NATOMS817:          DO J2=1,NATOMS
863:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 818:             SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))
864:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))819:             GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)
865:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)820:             GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3)
866:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3) 
867:             ENDIF 
868:          ENDDO821:          ENDDO
869:       ENDIF822:       ENDIF
870:    ENDDO823:    ENDDO
871: ENDIF824: ENDIF
872:          IF (PRINTE) THEN825:          IF (PRINTE) THEN
873:             WRITE(MYUNIT, '(A,G20.10)') 'ESPRING=',ESPRING826:             WRITE(MYUNIT, '(A,G20.10)') 'ESPRING=',ESPRING
874:          ENDIF827:          ENDIF
875: !828: !
876: ! Set gradients on frozen atoms to zero.829: ! Set gradients on frozen atoms to zero.
877: !830: !
907:    GGG(1:(3*NATOMS))=0.0D0860:    GGG(1:(3*NATOMS))=0.0D0
908:    GGG((INTIMAGE+1)*(3*NATOMS)+1:(INTIMAGE+2)*(3*NATOMS))=0.0D0861:    GGG((INTIMAGE+1)*(3*NATOMS)+1:(INTIMAGE+2)*(3*NATOMS))=0.0D0
909: ENDIF862: ENDIF
910: RMS=0.0D0863: RMS=0.0D0
911: RMSIMAGE(1:INTIMAGE+2)=0.0D0864: RMSIMAGE(1:INTIMAGE+2)=0.0D0
912: DO J1=2,INTIMAGE+1865: DO J1=2,INTIMAGE+1
913:    DO J2=1,(3*NATOMS)866:    DO J2=1,(3*NATOMS)
914:       RMSIMAGE(J1)=RMSIMAGE(J1)+GGG((3*NATOMS)*(J1-1)+J2)**2867:       RMSIMAGE(J1)=RMSIMAGE(J1)+GGG((3*NATOMS)*(J1-1)+J2)**2
915:    ENDDO868:    ENDDO
916:    RMS=RMS+RMSIMAGE(J1)869:    RMS=RMS+RMSIMAGE(J1)
917:    IF (LPRINT) WRITE(MYUNIT, '(A,I6,2G20.10)') ' congrad2> J1,EEE,RMSIMAGE=',J1,EEE(J1),RMSIMAGE(J1)870:    IF (LPRINT) WRITE(MYUNIT, '(A,I6,2G20.10,L5)') ' congrad2> J1,EEE,RMSIMAGE,freeze=', &
 871:   &                                                  J1,EEE(J1),RMSIMAGE(J1),IMGFREEZE(J1)
918: ENDDO872: ENDDO
919: IF (INTIMAGE.NE.0) THEN873: IF (INTIMAGE.NE.0) THEN
920:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))874:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))
921: ENDIF875: ENDIF
922: !876: !
923: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,877: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,
924: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)878: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)
925: !879: !
926: ETOTAL=0.0D0880: ETOTAL=0.0D0
927: MAXINT=-1.0D100881: MAXINT=-1.0D100


r29153/GMINdump.f90 2015-11-17 23:33:35.144976083 +0000 r29152/GMINdump.f90 2015-11-17 23:33:37.629009392 +0000
 27:    CALL FLUSH(MYUNIT) 27:    CALL FLUSH(MYUNIT)
 28: !   CALL SYSTEM('cp "GMIN.dump."//trim(adjustl(istr)) "GMIN.dump."//trim(adjustl(istr))//".save"') 28: !   CALL SYSTEM('cp "GMIN.dump."//trim(adjustl(istr)) "GMIN.dump."//trim(adjustl(istr))//".save"')
 29:    OPEN(MYUNIT2,FILE="GMIN.dump."//trim(adjustl(istr)),STATUS='UNKNOWN') 29:    OPEN(MYUNIT2,FILE="GMIN.dump."//trim(adjustl(istr)),STATUS='UNKNOWN')
 30: ELSE 30: ELSE
 31:    MYUNIT2=GETUNIT() 31:    MYUNIT2=GETUNIT()
 32:    WRITE (ISTR, '(A)') '' 32:    WRITE (ISTR, '(A)') ''
 33:    CALL SYSTEM('cp GMIN.dump GMIN.dump.save') 33:    CALL SYSTEM('cp GMIN.dump GMIN.dump.save')
 34:    OPEN(MYUNIT2,FILE='GMIN.dump',STATUS='UNKNOWN') 34:    OPEN(MYUNIT2,FILE='GMIN.dump',STATUS='UNKNOWN')
 35: ENDIF 35: ENDIF
 36:  36: 
 37:   WRITE(MYUNIT2, '(A)') 'steps completed NQ(JP) in mc' 37:   WRITE(MYUNIT2, '(A)') 'steps completed J1 in mc'
 38:   WRITE(MYUNIT2, '(I8)') NDONE 38:   WRITE(MYUNIT2, '(I8)') NDONE
 39:   WRITE(MYUNIT2, '(I8)') NSAVE 39:   WRITE(MYUNIT2, '(I8)') NSAVE
 40: ! Number of potential energy calls done so far 40: ! Number of potential energy calls done so far
 41:   WRITE(MYUNIT2, '(I20)') NPCALL  41:   WRITE(MYUNIT2, '(I20)') NPCALL 
 42:   WRITE(MYUNIT2, '(A)') 'COORDS' 42:   WRITE(MYUNIT2, '(A)') 'COORDS'
 43:   WRITE(MYUNIT2, '(A,I8)') 'run number ',JP 43:   WRITE(MYUNIT2, '(A,I8)') 'run number ',JP
 44:   WRITE(MYUNIT2, '(I8)') NATOMS 44:   WRITE(MYUNIT2, '(I8)') NATOMS
 45:   CALL WRITE_MARKOV_COORDS(MYUNIT2, '(3F25.15)', JP) 45:   CALL WRITE_MARKOV_COORDS(MYUNIT2, '(3F25.15)', JP)
 46: !  WRITE(MYUNIT2, '(3F25.15)') COORDSO(1:3*NATOMS,JP) 46: !  WRITE(MYUNIT2, '(3F25.15)') COORDSO(1:3*NATOMS,JP)
 47:   WRITE(MYUNIT2, '(A)') 'STEP, ASTEP, OSTEP, TEMP:' 47:   WRITE(MYUNIT2, '(A)') 'STEP, ASTEP, OSTEP, TEMP:'


r29153/intlbfgs.f90 2015-11-17 23:33:35.712983703 +0000 r29152/intlbfgs.f90 2015-11-17 23:33:38.193016958 +0000
 10: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 10: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 11: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 11: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 12: !   GNU General Public License for more details. 12: !   GNU General Public License for more details.
 13: ! 13: !
 14: !   You should have received a copy of the GNU General Public License 14: !   You should have received a copy of the GNU General Public License
 15: !   along with this program; if not, write to the Free Software 15: !   along with this program; if not, write to the Free Software
 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 17: ! 17: !
 18: SUBROUTINE INTLBFGS(QSTART,QFINISH) 18: SUBROUTINE INTLBFGS(QSTART,QFINISH)
 19: USE PORFUNCS 19: USE PORFUNCS
 20: USE COMMONS, ONLY : FREEZENODEST, FREEZETOL, MAXBFGS, CHRMMT, MYUNIT, CQMAX, & 20: USE COMMONS, ONLY : FREEZENODEST, MAXBFGS, CHRMMT, MYUNIT, CQMAX, &
 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, MUPDATE, INTDGUESS, & 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, MUPDATE, INTDGUESS, &
 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, & 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, &
 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, & 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, &
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, NATOMS, DEBUG, STEP, MCSTEPS 30:      & CONCUT, CONCUTLOCAL, NATOMS, DEBUG, STEP, MCSTEPS
 32: IMPLICIT NONE  32: IMPLICIT NONE 
 33:  33: 
 34: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points 34: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points
 35: INTEGER D, U 35: INTEGER D, U
 36: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP 36: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP
 37: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE 37: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE
 38: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS) 38: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS)
 39: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 39: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 40:  40: 
 41: DOUBLE PRECISION DUMMY, DPRAND 41: DOUBLE PRECISION DUMMY, DPRAND
 42: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM 42: INTEGER POINT,NPT,J3,J4,NACTIVE,NBEST,NEWATOM
 43: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE 43: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE
 44: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX 44: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX
 45: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS) 45: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS)
 46: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, & 46: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, &
 47:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), & 47:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), &
 48:   &                 BESTWORST, WORST 48:   &                 BESTWORST, WORST
 49: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA 49: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA
 50: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS) 50: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS)
 51: LOGICAL SWITCHED 51: LOGICAL SWITCHED
 52: DOUBLE PRECISION, POINTER :: X(:), G(:) 52: DOUBLE PRECISION, POINTER :: X(:), G(:)
 53: ! 53: !
 54: ! efk: for freezenodes 
 55: ! 
 56: DOUBLE PRECISION :: TESTG, TOTGNORM 
 57: INTEGER :: IM 
 58: ! 
 59: ! Dimensions involving INTIMAGE 54: ! Dimensions involving INTIMAGE
 60: ! 55: !
 61: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), & 56: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), &
 62:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), & 57:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), &
 63:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:) 58:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:)
 64: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:) 59: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:)
 65: ! saved interpolation 60: ! saved interpolation
 66: DOUBLE PRECISION, ALLOCATABLE :: BESTXYZ(:), BESTEEE(:) 61: DOUBLE PRECISION, ALLOCATABLE :: BESTXYZ(:), BESTEEE(:)
 67: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE 62: INTEGER BESTINTIMAGE, NSTEPS
 68: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:) 63: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:)
 69:  64: 
 70: ALLOCATE(TRUEEE(INTIMAGE+2), & 65: ALLOCATE(TRUEEE(INTIMAGE+2), &
 71:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), & 66:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &
 72:   &      GTMP(3*NATOMS*INTIMAGE), & 67:   &      GTMP(3*NATOMS*INTIMAGE), &
 73:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), & 68:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), &
 74:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), & 69:   &      GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &
 75:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), & 70:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &
 76:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE)) 71:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))
 77: ALLOCATE(BESTXYZ((3*NATOMS)*(INTIMAGE+2)),BESTEEE(INTIMAGE+2)) 72: ALLOCATE(BESTXYZ((3*NATOMS)*(INTIMAGE+2)),BESTEEE(INTIMAGE+2))
 78:  73: 
 79: SWITCHED=.FALSE. 74: SWITCHED=.FALSE.
 80: INTIMAGESAVE=INTIMAGE 75: INTIMAGESAVE=INTIMAGE
 81: NBACKTRACK=1 76: NBACKTRACK=1
 82: CALL MYCPU_TIME(STIME,.FALSE.) 77: CALL MYCPU_TIME(STIME,.FALSE.)
 83: WRITE(MYUNIT,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1 78: WRITE(MYUNIT,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1
 84: WRITE(MYUNIT,'(A,I6)') ' intlbfgs> Updates: ',MUPDATE 
 85: ADDATOM=.FALSE. 79: ADDATOM=.FALSE.
 86: NFAIL=0 80: NFAIL=0
 87: IMGFREEZE(1:INTIMAGE)=.FALSE. 81: IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.
 88: D=(3*NATOMS)*INTIMAGE 82: D=(3*NATOMS)*INTIMAGE
 89: U=MUPDATE 83: U=MUPDATE
 90: NITERDONE=1 84: NITERDONE=1
 91: NITERUSE=1 
 92: NQDONE=0 85: NQDONE=0
 93:  86: 
 94: IF ( D<=0 ) THEN 87: IF ( D<=0 ) THEN
 95:    WRITE(MYUNIT,*) 'd is not positive, d=',d 88:    WRITE(MYUNIT,*) 'd is not positive, d=',d
 96:    STOP 89:    STOP
 97: ENDIF 90: ENDIF
 98: IF ( U<=0 ) THEN 91: IF ( U<=0 ) THEN
 99:    WRITE(MYUNIT,*) 'u is not positive, u=',u 92:    WRITE(MYUNIT,*) 'u is not positive, u=',u
100:    STOP 93:    STOP
101: ENDIF 94: ENDIF
114: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))107: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
115: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))108: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
116: !109: !
117: ! Initialise XYZ110: ! Initialise XYZ
118: !111: !
119: XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))112: XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))
120: XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS))113: XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS))
121: DO J1=1,INTIMAGE+2114: DO J1=1,INTIMAGE+2
122:    XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)115:    XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)
123: ENDDO116: ENDDO
124:       WRITE(MYUNIT,'(A)') 'intlbfgs> here Z' 
125:       WRITE(MYUNIT,'(6G20.10)') XYZ(3*(398-1)+1:3*(398-1)+3), & 
126:   &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(398-1)+1:(INTIMAGE+1)*3*NATOMS+3*(398-1)+3) 
127:       WRITE(MYUNIT,'(6G20.10)') XYZ(3*(400-1)+1:3*(400-1)+3), & 
128:   &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(400-1)+1:(INTIMAGE+1)*3*NATOMS+3*(400-1)+3) 
129:       WRITE(MYUNIT,'(6G20.10)') QSTART(1:6) 
130:       WRITE(MYUNIT,'(6G20.10)') QFINISH(1:6) 
131: 117: 
132: NQCIFREEZE=0118: NQCIFREEZE=0
133: IF (FREEZE) THEN119: IF (FREEZE) THEN
134:    WRITE(MYUNIT,'(A)') ' intlbfgs> ERROR *** QCI has not been coded for frozen atoms yet'120:    WRITE(MYUNIT,'(A)') ' intlbfgs> ERROR *** QCI has not been coded for frozen atoms yet'
135:    STOP     121:    STOP     
136: ENDIF122: ENDIF
137: IF (ALLOCATED(INTFROZEN)) DEALLOCATE(INTFROZEN)123: IF (ALLOCATED(INTFROZEN)) DEALLOCATE(INTFROZEN)
138: ALLOCATE(INTFROZEN(NATOMS))124: ALLOCATE(INTFROZEN(NATOMS))
139: INTFROZEN(1:NATOMS)=.FALSE.125: INTFROZEN(1:NATOMS)=.FALSE.
140: DLIST(1:NATOMS)=-1126: DLIST(1:NATOMS)=-1
141: DMOVED(1:NATOMS)=1.0D100127: DMOVED(1:NATOMS)=1.0D100
142: IF (INTFREEZET) THEN128: IF (INTFREEZET) THEN
143:    DUMMY=INTFREEZETOL**2129:    DUMMY=INTFREEZETOL**2
144:    WRITE(MYUNIT,'(A,6G20.10)') ' intlbfgs> INTFREEZETOL,DUMMY=',INTFREEZETOL,DUMMY 
145:    DO J1=1,NATOMS130:    DO J1=1,NATOMS
146:       DF=(XYZ(3*(J1-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+1))**2 &131:       DF=(XYZ(3*(J1-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+1))**2 &
147:   &     +(XYZ(3*(J1-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+2))**2 &132:   &     +(XYZ(3*(J1-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+2))**2 &
148:   &     +(XYZ(3*(J1-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+3))**2133:   &     +(XYZ(3*(J1-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+3))**2
149:       IF (DF.LT.DUMMY) THEN134:       IF (DF.LT.DUMMY) THEN
150:          NQCIFREEZE=NQCIFREEZE+1135:          NQCIFREEZE=NQCIFREEZE+1
151:          INTFROZEN(J1)=.TRUE.136:          INTFROZEN(J1)=.TRUE.
152:          IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6)') &137: !        IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6)') &
153:   &            ' intlbfgs> atom ',J1,' moves less than threshold: dist^2=',DF,' total=',NQCIFREEZE138: ! &            ' intlbfgs> atom ',J1,' moves less than threshold: dist^2=',DF,' total=',NQCIFREEZE
154:       ENDIF139:       ENDIF
155:       sortd: DO J2=1,J1140:       sortd: DO J2=1,J1
156:          IF (DF.LT.DMOVED(J2)) THEN141:          IF (DF.LT.DMOVED(J2)) THEN
157:             DO J3=J1,J2+1,-1142:             DO J3=J1,J2+1,-1
158:                DMOVED(J3)=DMOVED(J3-1)143:                DMOVED(J3)=DMOVED(J3-1)
159:                DLIST(J3)=DLIST(J3-1)144:                DLIST(J3)=DLIST(J3-1)
160:             ENDDO145:             ENDDO
161:             DMOVED(J2)=DF146:             DMOVED(J2)=DF
162:             DLIST(J2)=J1147:             DLIST(J2)=J1
163:             EXIT sortd148:             EXIT sortd
164:          ENDIF149:          ENDIF
165:       ENDDO sortd150:       ENDDO sortd
166:    ENDDO151:    ENDDO
167:    WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6)') ' intlbfgs> Total number of atoms moving less than threshold=',NQCIFREEZE152:    WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6)') ' intlbfgs> Total number of atoms moving less than threshold=',NQCIFREEZE
168: ENDIF153: ENDIF
169: 154: 
170:       WRITE(MYUNIT,'(6G20.10)') XYZ(3*(398-1)+1:3*(398-1)+3), & 
171:   &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(398-1)+1:(INTIMAGE+1)*3*NATOMS+3*(398-1)+3) 
172:       WRITE(MYUNIT,'(6G20.10)') XYZ(3*(400-1)+1:3*(400-1)+3), & 
173:   &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(400-1)+1:(INTIMAGE+1)*3*NATOMS+3*(400-1)+3) 
174:  
175: IF (NATOMS-NQCIFREEZE.LT.INTFREEZEMIN) THEN155: IF (NATOMS-NQCIFREEZE.LT.INTFREEZEMIN) THEN
176:    DO J1=NATOMS,NATOMS-INTFREEZEMIN+1,-1156:    DO J1=NATOMS,NATOMS-INTFREEZEMIN+1,-1
177:       INTFROZEN(DLIST(J1))=.FALSE.157:       INTFROZEN(DLIST(J1))=.FALSE.
178:    ENDDO158:    ENDDO
179:    NQCIFREEZE=NATOMS-INTFREEZEMIN159:    NQCIFREEZE=NATOMS-INTFREEZEMIN
180:    WRITE(MYUNIT,'(A,I6,A)') ' intlbfgs> Freezing ',NQCIFREEZE,' atoms'160:    WRITE(MYUNIT,'(A,I6,A)') ' intlbfgs> Freezing ',NQCIFREEZE,' atoms'
181: ENDIF161: ENDIF
182: 162: 
183: NLASTGOODE=0163: NLASTGOODE=0
184: LASTGOODE=1.0D100164: LASTGOODE=1.0D100
261:    XSAVE(1:D)=X(1:D)241:    XSAVE(1:D)=X(1:D)
262:    GOTO 567242:    GOTO 567
263: ENDIF243: ENDIF
264: DO J1=1,NCONSTRAINT244: DO J1=1,NCONSTRAINT
265:    DF=SQRT((XYZ(3*(CONI(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+1))**2 &245:    DF=SQRT((XYZ(3*(CONI(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+1))**2 &
266:   &       +(XYZ(3*(CONI(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+2))**2 &246:   &       +(XYZ(3*(CONI(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+2))**2 &
267:   &       +(XYZ(3*(CONI(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+3))**2)&247:   &       +(XYZ(3*(CONI(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+3))**2)&
268:   &  +SQRT((XYZ(3*(CONJ(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+1))**2 &248:   &  +SQRT((XYZ(3*(CONJ(J1)-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+1))**2 &
269:   &       +(XYZ(3*(CONJ(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+2))**2 &249:   &       +(XYZ(3*(CONJ(J1)-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+2))**2 &
270:   &       +(XYZ(3*(CONJ(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+3))**2)250:   &       +(XYZ(3*(CONJ(J1)-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+3))**2)
271:    IF (J1.EQ.3505) THEN 
272:       WRITE(MYUNIT,'(A,3I10)') 'intlbfgs> J1,CONI(J1),CONJ(J1)=',J1,CONI(J1),CONJ(J1) 
273:       WRITE(MYUNIT,'(6G20.10)') XYZ(3*(CONI(J1)-1)+1:3*(CONI(J1)-1)+3), & 
274:   &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+1:(INTIMAGE+1)*3*NATOMS+3*(CONI(J1)-1)+3) 
275:       WRITE(MYUNIT,'(6G20.10)') XYZ(3*(CONJ(J1)-1)+1:3*(CONJ(J1)-1)+3), & 
276:   &                             XYZ((INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+1:(INTIMAGE+1)*3*NATOMS+3*(CONJ(J1)-1)+3) 
277:    ENDIF 
278:    IF (DF.LT.DUMMY) THEN251:    IF (DF.LT.DUMMY) THEN
279:       NBEST=J1252:       NBEST=J1
280:       DUMMY=DF253:       DUMMY=DF
281:    ENDIF254:    ENDIF
282: ENDDO255: ENDDO
283: IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,2I6,A,F15.5)') ' intlbfgs> Smallest overall motion for constraint ',NBEST,' atoms ', &256: IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,2I6,A,F15.5)') ' intlbfgs> Smallest overall motion for constraint ',NBEST,' atoms ', &
284:   &                           CONI(NBEST),CONJ(NBEST),' distance=',DUMMY257:   &                           CONI(NBEST),CONJ(NBEST),' distance=',DUMMY
285: 258: 
286: TURNONORDER(1:NATOMS)=0259: TURNONORDER(1:NATOMS)=0
287: NTRIES(1:NATOMS)=1260: NTRIES(1:NATOMS)=1
440: !  *** NEW Find atom with most constraints to active set413: !  *** NEW Find atom with most constraints to active set
441: !  Turn on constraint terms for this atom with all previous members of the active set414: !  Turn on constraint terms for this atom with all previous members of the active set
442: !  Add repulsions to non-constrained atoms in this set415: !  Add repulsions to non-constrained atoms in this set
443: !  NTOADD is the number of atoms to add to the active set in each pass. 1 seems best!416: !  NTOADD is the number of atoms to add to the active set in each pass. 1 seems best!
444: !417: !
445:    IF (ADDATOM.AND.(NACTIVE.LT.NATOMS)) THEN418:    IF (ADDATOM.AND.(NACTIVE.LT.NATOMS)) THEN
446:       CALL DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)419:       CALL DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)
447:       NLASTGOODE=NITERDONE420:       NLASTGOODE=NITERDONE
448:       LASTGOODE=ETOTAL421:       LASTGOODE=ETOTAL
449:    ENDIF422:    ENDIF
450:    GTMP(1:D)=0.0D0423:    CALL MAKESTEP(NITERDONE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)
451:    CALL MAKESTEP(NITERUSE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA) 
452: !424: !
453: ! If the number of images has changed since G was declared then G is not the same425: ! If the number of images has changed since G was declared then G is not the same
454: ! size as Gtmp and Dot_Product cannot be used.426: ! size as Gtmp and Dot_Product cannot be used.
455: !427: !
456: !  IF (Dot_Product(G,Gtmp)/SQRT( Dot_Product(G,G)*Dot_Product(Gtmp,Gtmp) ) > 0.0D0) THEN428: !  IF (Dot_Product(G,Gtmp)/SQRT( Dot_Product(G,G)*Dot_Product(Gtmp,Gtmp) ) > 0.0D0) THEN
457: !429: !
458: !  Separate sqrt;s to avoid overflow.430: !  Separate sqrt;s to avoid overflow.
459: !431: !
460:    IF (DDOT(D,G,1,GTMP,1)/MAX(1.0D-100,SQRT( DDOT(D,G,1,G,1))*SQRT(DDOT(D,GTMP,1,GTMP,1)) ) > 0.0D0) THEN432:    IF (DDOT(D,G,1,GTMP,1)/MAX(1.0D-100,SQRT( DDOT(D,G,1,G,1))*SQRT(DDOT(D,GTMP,1,GTMP,1)) ) > 0.0D0) THEN
461:         IF (DEBUG) WRITE(MYUNIT,*) 'Search direction has positive projection onto gradient - reversing step'433:         IF (DEBUG) WRITE(MYUNIT,*) 'Search direction has positive projection onto gradient - reversing step'
470: 442: 
471:    STPMIN=1.0D0443:    STPMIN=1.0D0
472:    DO J2=1,INTIMAGE444:    DO J2=1,INTIMAGE
473:       STEPIMAGE(J2) = SQRT(DOT_PRODUCT(SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2), &445:       STEPIMAGE(J2) = SQRT(DOT_PRODUCT(SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2), &
474:   &                                    SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2)))446:   &                                    SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2)))
475:       DUMMY=STEPIMAGE(J2)447:       DUMMY=STEPIMAGE(J2)
476:       IF (STEPIMAGE(J2) > MAXBFGS) THEN448:       IF (STEPIMAGE(J2) > MAXBFGS) THEN
477:            STP((3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2) = MAXBFGS/STEPIMAGE(J2)449:            STP((3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2) = MAXBFGS/STEPIMAGE(J2)
478:            STPMIN=MIN(STPMIN,STP((3*NATOMS)*(J2-1)+1))450:            STPMIN=MIN(STPMIN,STP((3*NATOMS)*(J2-1)+1))
479:       ENDIF451:       ENDIF
480: !     WRITE(MYUNIT,'(A,I8,3G20.10)') ' image,initial step size,STP,prod=',J2,DUMMY,STP(3*NATOMS*(J2-1)+1), & 
481: ! &                                   STEPIMAGE(J2)*STP(3*NATOMS*(J2-1)+1)    
482:    ENDDO452:    ENDDO
483:    STP(1:D)=STPMIN453:    STP(1:D)=STPMIN
484: ! EFK: decide whether to freeze some nodes454: 
485:    IF (FREEZENODEST) THEN 
486:       TOTGNORM=SQRT(DOT_PRODUCT(G(1:(3*NATOMS)*INTIMAGE),G(1:(3*NATOMS)*INTIMAGE))/INTIMAGE) 
487:       NIMAGEFREEZE=0 
488:       DO IM=1,INTIMAGE 
489:          TESTG=SQRT(DOT_PRODUCT(G((3*NATOMS)*(IM-1)+1:(3*NATOMS)*IM),G((3*NATOMS)*(IM-1)+1:(3*NATOMS)*IM))) 
490:          IMGFREEZE(IM)=.FALSE. 
491:          IF (TOTGNORM.NE.0.0D0) THEN 
492: !           IF (TESTG/TOTGNORM.LT.FREEZETOL) THEN 
493:             IF (TESTG/SQRT(3.0D0*NATOMS).LT.FREEZETOL) THEN 
494: !              IF (DEBUG) PRINT '(A,I6,3G20.10)', ' intlbfgs> Freezing image: ',IM,TESTG,FREEZETOL,TOTGNORM 
495:                IMGFREEZE(IM)=.TRUE. 
496:                STEPIMAGE(IM)=0.0D0 
497:                NIMAGEFREEZE=NIMAGEFREEZE+1 
498:                STP((3*NATOMS)*(IM-1)+1:(3*NATOMS)*IM)=0.0D0 
499:             ENDIF 
500:          ENDIF 
501:       ENDDO 
502:       IF (DEBUG) PRINT '(2(A,I6))', ' intlbfgs> Number of frozen images=',NIMAGEFREEZE,' / ',INTIMAGE 
503:    ENDIF 
504:    !  We now have the proposed step - update geometry and calculate new gradient455:    !  We now have the proposed step - update geometry and calculate new gradient
505:    NDECREASE=0456:    NDECREASE=0
506: 20 X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)457: 20 X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
507: 458: 
508: !  IF (.NOT.SWITCHED) THEN459: !  IF (.NOT.SWITCHED) THEN
509:    IF (.TRUE.) THEN460:    IF (.TRUE.) THEN
510: !     IF ((RMS.LT.INTRMSTOL*1.0D10).AND.(MOD(NITERDONE,10).EQ.0).AND.(NSTEPSMAX-NITERDONE.GT.100)) &461: !     IF ((RMS.LT.INTRMSTOL*1.0D10).AND.(MOD(NITERDONE,10).EQ.0).AND.(NSTEPSMAX-NITERDONE.GT.100)) &
511: ! &               CALL CHECKSEP(NMAXINT,NMININT,INTIMAGE,XYZ,(3*NATOMS),NATOMS)462: ! &               CALL CHECKSEP(NMAXINT,NMININT,INTIMAGE,XYZ,(3*NATOMS),NATOMS)
512:       IF (MOD(NITERDONE,INTIMAGECHECK).EQ.0) THEN463:       IF (MOD(NITERDONE,INTIMAGECHECK).EQ.0) THEN
513: 864      CONTINUE ! for adding more than one image at a time464: 864      CONTINUE ! for adding more than one image at a time
522:             ENDDO473:             ENDDO
523:             DUMMY=SQRT(DUMMY)474:             DUMMY=SQRT(DUMMY)
524:             IF (DUMMY.GT.DMAX) THEN475:             IF (DUMMY.GT.DMAX) THEN
525:                DMAX=DUMMY476:                DMAX=DUMMY
526:                JMAX=J1477:                JMAX=J1
527:             ENDIF478:             ENDIF
528:             IF (DUMMY.LT.DMIN) THEN479:             IF (DUMMY.LT.DMIN) THEN
529:                DMIN=DUMMY480:                DMIN=DUMMY
530:                JMIN=J1481:                JMIN=J1
531:             ENDIF482:             ENDIF
532:             IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &483: !           IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &
533:   &                                                  J1,' and ',J1+1,' is ',DUMMY484: ! &                                                  J1,' and ',J1+1,' is ',DUMMY
534:          ENDDO485:          ENDDO
535:          IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN486:          IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN
536:             WRITE(MYUNIT,'(A,I6,A,I6,A,I6)') ' intlbfgs> Add an image between ',JMAX,' and ',JMAX+1,' INTIMAGE=',INTIMAGE487:             WRITE(MYUNIT,'(A,I6,A,I6,A,I6)') ' intlbfgs> Add an image between ',JMAX,' and ',JMAX+1,' INTIMAGE=',INTIMAGE
537:             NITERUSE=0 
538:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))488:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
539:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))489:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
540:             DEALLOCATE(XYZ)490:             DEALLOCATE(XYZ)
541:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))491:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))
542:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)492:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)
543:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &493:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &
544:   &                                               + DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1)))/2.0D0494:   &                                               + DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1)))/2.0D0
545:             XYZ(3*NATOMS*(JMAX+1)+1:3*NATOMS*(INTIMAGE+3))=DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+2))495:             XYZ(3*NATOMS*(JMAX+1)+1:3*NATOMS*(INTIMAGE+3))=DPTMP(3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+2))
546: !496: !
547: ! Save step-taking memories in SEARCHSTEP and GDIF.497: ! Save step-taking memories in SEARCHSTEP and GDIF.
553:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)503:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
554:             DEALLOCATE(SEARCHSTEP)504:             DEALLOCATE(SEARCHSTEP)
555:             ALLOCATE(SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))505:             ALLOCATE(SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))
556:             DO J1=0,MUPDATE506:             DO J1=0,MUPDATE
557:                IF (JMAX.GT.1) SEARCHSTEP(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))507:                IF (JMAX.GT.1) SEARCHSTEP(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))
558:                IF (JMAX.LT.INTIMAGE+1) SEARCHSTEP(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &508:                IF (JMAX.LT.INTIMAGE+1) SEARCHSTEP(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &
559:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)509:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)
560:                SEARCHSTEP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &510:                SEARCHSTEP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &
561:   &                             D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))511:   &                             D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))
562:             ENDDO512:             ENDDO
563: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
564:             SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1))=0.0D0 
565: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
566:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)513:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
567:             DEALLOCATE(GDIF)514:             DEALLOCATE(GDIF)
568:             ALLOCATE(GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))515:             ALLOCATE(GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1)))
569:             DO J1=0,MUPDATE516:             DO J1=0,MUPDATE
570:                IF (JMAX.GT.1) GDIF(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))517:                IF (JMAX.GT.1) GDIF(J1,1:3*NATOMS*(JMAX-1))=D2TMP(J1,1:3*NATOMS*(JMAX-1))
571:                IF (JMAX.LT.INTIMAGE+1) GDIF(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &518:                IF (JMAX.LT.INTIMAGE+1) GDIF(J1,3*NATOMS*JMAX+1:3*NATOMS*(INTIMAGE+1))= &
572:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)519:   &                 D2TMP(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*INTIMAGE)
573:                GDIF(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &520:                GDIF(J1,3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX)= &
574:   &                       D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))521:   &                       D2TMP(J1,3*NATOMS*(MIN(JMAX,INTIMAGE)-1)+1:3*NATOMS*MIN(JMAX,INTIMAGE))
575:             ENDDO522:             ENDDO
576: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
577:             GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE+1))=0.0D0 
578: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
579:             DEALLOCATE(D2TMP)523:             DEALLOCATE(D2TMP)
580: 524: 
581:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &525:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &
582:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)526:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)
583:             ALLOCATE(TRUEEE(INTIMAGE+3), &527:             ALLOCATE(TRUEEE(INTIMAGE+3), &
584:   &                  EEETMP(INTIMAGE+3), MYGTMP(3*NATOMS*(INTIMAGE+1)), &528:   &                  EEETMP(INTIMAGE+3), MYGTMP(3*NATOMS*(INTIMAGE+1)), &
585:   &                  GTMP(3*NATOMS*(INTIMAGE+1)), &529:   &                  GTMP(3*NATOMS*(INTIMAGE+1)), &
586:   &                  DIAG(3*NATOMS*(INTIMAGE+1)), STP(3*NATOMS*(INTIMAGE+1)), &530:   &                  DIAG(3*NATOMS*(INTIMAGE+1)), STP(3*NATOMS*(INTIMAGE+1)), &
587:   &                  GLAST((3*NATOMS)*(INTIMAGE+1)), &531:   &                  GLAST((3*NATOMS)*(INTIMAGE+1)), &
588:   &                  XSAVE((3*NATOMS)*(INTIMAGE+1)), CHECKG((3*NATOMS)*(INTIMAGE+1)), IMGFREEZE(INTIMAGE+1), &532:   &                  XSAVE((3*NATOMS)*(INTIMAGE+1)), CHECKG((3*NATOMS)*(INTIMAGE+1)), IMGFREEZE(INTIMAGE+1), &
607:             D=(3*NATOMS)*INTIMAGE551:             D=(3*NATOMS)*INTIMAGE
608:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)552:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
609:             IF (CHECKCONINT) THEN553:             IF (CHECKCONINT) THEN
610:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)554:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
611:             ELSE555:             ELSE
612:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)556:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
613:             ENDIF557:             ENDIF
614: !           GOTO 864558: !           GOTO 864
615:          ELSEIF ((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1)) THEN559:          ELSEIF ((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1)) THEN
616:             IF (JMIN.EQ.1) JMIN=2560:             IF (JMIN.EQ.1) JMIN=2
617:             WRITE(MYUNIT,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN561: !           WRITE(MYUNIT,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN
618:             NITERUSE=0 
619:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))562:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
620:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))563:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
621:             DEALLOCATE(XYZ)564:             DEALLOCATE(XYZ)
622:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+1)))565:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+1)))
623:             XYZ(1:3*NATOMS*(JMIN-1))=DPTMP(1:3*NATOMS*(JMIN-1))566:             XYZ(1:3*NATOMS*(JMIN-1))=DPTMP(1:3*NATOMS*(JMIN-1))
624:             XYZ(3*NATOMS*(JMIN-1)+1:3*NATOMS*(INTIMAGE+1))=DPTMP(3*NATOMS*JMIN+1:3*NATOMS*(INTIMAGE+2))567:             XYZ(3*NATOMS*(JMIN-1)+1:3*NATOMS*(INTIMAGE+1))=DPTMP(3*NATOMS*JMIN+1:3*NATOMS*(INTIMAGE+2))
625: 568: 
626:             DEALLOCATE(DPTMP)569:             DEALLOCATE(DPTMP)
627: !570: !
628: ! Save step-taking memories in SEARCHSTEP and GDIF.571: ! Save step-taking memories in SEARCHSTEP and GDIF.
631: !574: !
632:             ALLOCATE(D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE))575:             ALLOCATE(D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE))
633:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)576:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
634:             DEALLOCATE(SEARCHSTEP)577:             DEALLOCATE(SEARCHSTEP)
635:             ALLOCATE(SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))578:             ALLOCATE(SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))
636:             DO J1=0,MUPDATE579:             DO J1=0,MUPDATE
637:                SEARCHSTEP(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))580:                SEARCHSTEP(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))
638:                SEARCHSTEP(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &581:                SEARCHSTEP(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &
639:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)582:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)
640:             ENDDO583:             ENDDO
641: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
642:             SEARCHSTEP(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1))=0.0D0 
643: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
644:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)584:             D2TMP(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)=GDIF(0:MUPDATE,1:(3*NATOMS)*INTIMAGE)
645:             DEALLOCATE(GDIF)585:             DEALLOCATE(GDIF)
646:             ALLOCATE(GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))586:             ALLOCATE(GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1)))
647:             DO J1=0,MUPDATE587:             DO J1=0,MUPDATE
648:                GDIF(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))588:                GDIF(J1,1:3*NATOMS*(JMIN-2))=D2TMP(J1,1:3*NATOMS*(JMIN-2))
649:                GDIF(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &589:                GDIF(J1,3*NATOMS*(JMIN-2)+1:3*NATOMS*(INTIMAGE-1))= &
650:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)590:   &                     D2TMP(J1,3*NATOMS*(JMIN-1)+1:3*NATOMS*INTIMAGE)
651:             ENDDO591:             ENDDO
652: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
653:             GDIF(0:MUPDATE,1:(3*NATOMS)*(INTIMAGE-1))=0.0D0 
654: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
655:             DEALLOCATE(D2TMP)592:             DEALLOCATE(D2TMP)
656: 593: 
657:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &594:             DEALLOCATE(TRUEEE,EEETMP,MYGTMP,GTMP,GGG, &
658:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)595:   &                    DIAG,STP,GLAST,XSAVE,EEE,STEPIMAGE,CHECKG,IMGFREEZE)
659:             ALLOCATE(TRUEEE(INTIMAGE+1),&596:             ALLOCATE(TRUEEE(INTIMAGE+1),&
660:   &                  EEETMP(INTIMAGE+1), MYGTMP(3*NATOMS*(INTIMAGE-1)), &597:   &                  EEETMP(INTIMAGE+1), MYGTMP(3*NATOMS*(INTIMAGE-1)), &
661:   &                  GTMP(3*NATOMS*(INTIMAGE-1)), &598:   &                  GTMP(3*NATOMS*(INTIMAGE-1)), &
662:   &                  DIAG(3*NATOMS*(INTIMAGE-1)), STP(3*NATOMS*(INTIMAGE-1)), &599:   &                  DIAG(3*NATOMS*(INTIMAGE-1)), STP(3*NATOMS*(INTIMAGE-1)), &
663:   &                  GLAST((3*NATOMS)*(INTIMAGE-1)), &600:   &                  GLAST((3*NATOMS)*(INTIMAGE-1)), &
664:   &                  XSAVE((3*NATOMS)*(INTIMAGE-1)), CHECKG((3*NATOMS)*(INTIMAGE-1)), IMGFREEZE(INTIMAGE-1), &601:   &                  XSAVE((3*NATOMS)*(INTIMAGE-1)), CHECKG((3*NATOMS)*(INTIMAGE-1)), IMGFREEZE(INTIMAGE-1), &
913:       ENDIF850:       ENDIF
914:       EXIT851:       EXIT
915:    ENDIF852:    ENDIF
916:    777 CONTINUE853:    777 CONTINUE
917: !854: !
918: ! Compute the new step and gradient change855: ! Compute the new step and gradient change
919: !856: !
920:    NPT=POINT*D857:    NPT=POINT*D
921:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)858:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)
922:    GDIF(POINT,:)=G-GTMP859:    GDIF(POINT,:)=G-GTMP
923:     
924:    POINT=POINT+1; IF (POINT==MUPDATE) POINT=0860:    POINT=POINT+1; IF (POINT==MUPDATE) POINT=0
925: 861: 
926:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL RWG(NITERDONE,INTIMAGE,XYZ)862:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL RWG(NITERDONE,INTIMAGE,XYZ)
927:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)863:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)
928: 864: 
929:    NITERDONE=NITERDONE+1865:    NITERDONE=NITERDONE+1
930:    NITERUSE=NITERUSE+1 
931: 866: 
932:    IF (NITERDONE.GT.NSTEPSMAX) EXIT867:    IF (NITERDONE.GT.NSTEPSMAX) EXIT
933:    IF (NACTIVE.EQ.NATOMS) THEN868:    IF (NACTIVE.EQ.NATOMS) THEN
934:       IF (.NOT.SWITCHED) THEN869:       IF (.NOT.SWITCHED) THEN
935:          CALL MYCPU_TIME(FTIME,.FALSE.)870:          CALL MYCPU_TIME(FTIME,.FALSE.)
936:          WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &871:          WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
937:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME872:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
938:          WRITE(MYUNIT,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'873:          WRITE(MYUNIT,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
939:          DO J1=1,NATOMS874:          DO J1=1,NATOMS
940:             IF (.NOT.ATOMACTIVE(J1)) THEN875:             IF (.NOT.ATOMACTIVE(J1)) THEN
970: ! Linear interpolation for constraint potential and real potential separately.905: ! Linear interpolation for constraint potential and real potential separately.
971: ! Constraint potential need not be flat if we have done some steps with both906: ! Constraint potential need not be flat if we have done some steps with both
972: ! potentials turned on.907: ! potentials turned on.
973: !908: !
974: 678 CONTINUE909: 678 CONTINUE
975: 910: 
976: CALL RWG(NITERDONE,INTIMAGE,XYZ)911: CALL RWG(NITERDONE,INTIMAGE,XYZ)
977: CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)912: CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)
978: NQDONE=NQDONE+1913: NQDONE=NQDONE+1
979: 914: 
980: WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> WORST=',WORST 
981: WRITE(MYUNIT,'(A,2I8)') 'intlbfgs> NQDONE,MCSTEPS=',NQDONE,MCSTEPS(1) 
982: IF (WORST.EQ.0.0D0) GOTO 8765915: IF (WORST.EQ.0.0D0) GOTO 8765
983: IF (NQDONE.EQ.MCSTEPS(1)) GOTO 8765916: IF (NQDONE.EQ.MCSTEPS(1)) GOTO 8765
984: 917: 
985: !918: !
986: ! Accept/reject this QCI set until BH steps exceeded, or worst energy is zero.919: ! Accept/reject this QCI set until BH steps exceeded, or worst energy is zero.
987: ! Reset everything necessary and go back to 9876 CONTINUE - start of QCI optimisation.920: ! Reset everything necessary and go back to 9876 CONTINUE - start of QCI optimisation.
988: !921: !
989: IF (WORST.GT.BESTWORST) THEN922: IF (WORST.GT.BESTWORST) THEN
990:    WRITE(MYUNIT,'(A)') 'intlbfgs> rejecting step - resetting' 
991:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &923:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
992:   &         DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)924:   &         DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
993:    NULLIFY(X,G)925:    NULLIFY(X,G)
994:    INTIMAGE=BESTINTIMAGE926:    INTIMAGE=BESTINTIMAGE
995:    D=(3*NATOMS)*INTIMAGE927:    D=(3*NATOMS)*INTIMAGE
996: 928: 
997:    ALLOCATE(TRUEEE(INTIMAGE+2), &929:    ALLOCATE(TRUEEE(INTIMAGE+2), &
998:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &930:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &
999:   &      GTMP(3*NATOMS*INTIMAGE), &931:   &      GTMP(3*NATOMS*INTIMAGE), &
1000:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), &932:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE), &
1969: ! IF (DEBUG) WRITE(MYUNIT,'(A,F15.5)') ' checkperc> Final constraint tolerance parameter ',LINTCONSTRAINTTOL1901: ! IF (DEBUG) WRITE(MYUNIT,'(A,F15.5)') ' checkperc> Final constraint tolerance parameter ',LINTCONSTRAINTTOL
1970: 1902: 
1971: ! WRITE(MYUNIT,'(A,I6,3(A,F15.5))') ' checkperc> Total distance constraints=',NCONSTRAINT, &1903: ! WRITE(MYUNIT,'(A,I6,3(A,F15.5))') ' checkperc> Total distance constraints=',NCONSTRAINT, &
1972: !   &                    ' shortest=',MINCONDIST,' longest=',MAXCONDIST,' tolerance=',LINTCONSTRAINTTOL1904: !   &                    ' shortest=',MINCONDIST,' longest=',MAXCONDIST,' tolerance=',LINTCONSTRAINTTOL
1973: 1905: 
1974: CALLED=.TRUE.1906: CALLED=.TRUE.
1975: 1907: 
1976: END SUBROUTINE CHECKPERC1908: END SUBROUTINE CHECKPERC
1977: 1909: 
1978: SUBROUTINE MAKESTEP(NITERDONE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)1910: SUBROUTINE MAKESTEP(NITERDONE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)
1979: USE COMMONS, ONLY : MUPDATE, INTDGUESS, NATOMS, MYUNIT1911: USE COMMONS, ONLY : MUPDATE, INTDGUESS, NATOMS
1980: IMPLICIT NONE1912: IMPLICIT NONE
1981: INTEGER NITERDONE, POINT, BOUND, NPT, D, CP, INTIMAGE, I1913: INTEGER NITERDONE, POINT, BOUND, NPT, D, CP, INTIMAGE, I
1982: DOUBLE PRECISION DIAG(3*NATOMS*INTIMAGE),SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE),G((3*NATOMS)*INTIMAGE), &1914: DOUBLE PRECISION DIAG(3*NATOMS*INTIMAGE),SEARCHSTEP(0:MUPDATE,(3*NATOMS)*INTIMAGE),G((3*NATOMS)*INTIMAGE), &
1983:   &  GTMP(3*NATOMS*INTIMAGE), GNORM, STP(3*NATOMS*INTIMAGE), YS, GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE), YY, &1915:   &  GTMP(3*NATOMS*INTIMAGE), GNORM, STP(3*NATOMS*INTIMAGE), YS, GDIF(0:MUPDATE,(3*NATOMS)*INTIMAGE), YY, &
1984:   &  SQ, YR, BETA1916:   &  SQ, YR, BETA
1985: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA1917: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA
1986: LOGICAL CHANGEIMAGE 
1987: SAVE1918: SAVE
1988: 1919: 
1989: ! WRITE(MYUNIT,'(A,I8)') 'makestep> NITERDONE=',NITERDONE 
1990: ! WRITE(MYUNIT,'(A,3G20.10)') 'makestep> DIAG: ',DIAG(118:120) 
1991: ! WRITE(MYUNIT,'(A,3G20.10)') 'makestep> INTDGUESS: ',INTDGUESS 
1992: ! WRITE(MYUNIT,'(A,3G20.10)') 'makestep> G: ',G(118:120) 
1993: ! WRITE(MYUNIT,'(A,4I8)') 'NPT,D,NPT/D,POINT=',NPT,D,NPT/D,POINT 
1994: ! WRITE(MYUNIT,'(A,3G20.10)') 'makestep> GDIF: ',GDIF(NPT/D,118:120) 
1995: ! WRITE(MYUNIT,'(A,3G20.10)') 'makestep> SEARCHSTEP: ',SEARCHSTEP(NPT/D,118:120) 
1996:  
1997: MAIN: IF (NITERDONE==1) THEN1920: MAIN: IF (NITERDONE==1) THEN
1998:      POINT = 01921:      POINT = 0
1999:      DIAG(1:D)=INTDGUESS1922:      DIAG(1:D)=INTDGUESS
2000:      SEARCHSTEP(0,1:D)= -G(1:D)*INTDGUESS            ! NR STEP FOR DIAGONAL INVERSE HESSIAN1923:      SEARCHSTEP(0,1:D)= -G(1:D)*INTDGUESS            ! NR STEP FOR DIAGONAL INVERSE HESSIAN
2001:      GTMP(1:D)        = SEARCHSTEP(0,1:D)1924:      GTMP(1:D)        = SEARCHSTEP(0,1:D)
2002:      GNORM            = MAX(SQRT(DOT_PRODUCT(G(1:D),G(1:D))),1.0D-100)1925:      GNORM            = MAX(SQRT(DOT_PRODUCT(G(1:D),G(1:D))),1.0D-100)
2003:      STP(1:D)         = MIN(1.0D0/GNORM, GNORM) ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS1926:      STP(1:D)         = MIN(1.0D0/GNORM, GNORM) ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS
2004: ELSE MAIN1927: ELSE MAIN
2005:      BOUND=NITERDONE-11928:      BOUND=NITERDONE-1
2006:      IF (NITERDONE.GT.MUPDATE) BOUND=MUPDATE1929:      IF (NITERDONE.GT.MUPDATE) BOUND=MUPDATE
2030:           ALPHA(CP+1) = RHO1(CP+1) * SQ1953:           ALPHA(CP+1) = RHO1(CP+1) * SQ
2031:           GTMP(1:D)        = -ALPHA(CP+1)*GDIF(CP,1:D) + GTMP(1:D)1954:           GTMP(1:D)        = -ALPHA(CP+1)*GDIF(CP,1:D) + GTMP(1:D)
2032:      ENDDO1955:      ENDDO
2033:               1956:               
2034:      GTMP(1:D)=DIAG(1)*GTMP(1:D)1957:      GTMP(1:D)=DIAG(1)*GTMP(1:D)
2035: 1958: 
2036:      DO I=1,BOUND1959:      DO I=1,BOUND
2037:           YR= DOT_PRODUCT( GDIF(CP,1:D) , GTMP )1960:           YR= DOT_PRODUCT( GDIF(CP,1:D) , GTMP )
2038:           BETA= RHO1(CP+1)*YR1961:           BETA= RHO1(CP+1)*YR
2039:           BETA= ALPHA(CP+1)-BETA1962:           BETA= ALPHA(CP+1)-BETA
2040: !         WRITE(MYUNIT,'(A,I8,4G20.10)') 'makestep> I,YR,BETA,RHO1,ALPHA=',I,YR,BETA,RHO1(CP+1),ALPHA(CP+1) 
2041:           GTMP(1:D) = BETA*SEARCHSTEP(CP,1:D) + GTMP(1:D)1963:           GTMP(1:D) = BETA*SEARCHSTEP(CP,1:D) + GTMP(1:D)
2042:           CP=CP+11964:           CP=CP+1
2043: !         IF (CP==M) CP=01965: !         IF (CP==M) CP=0
2044:           IF (CP==MUPDATE) CP=01966:           IF (CP==MUPDATE) CP=0
2045:      ENDDO1967:      ENDDO
2046:               1968:               
2047:      STP(1:D) = 1.0D01969:      STP(1:D) = 1.0D0
2048: ENDIF MAIN1970: ENDIF MAIN
2049: 1971: 
2050: !  Store the new search direction1972: !  Store the new search direction


r29153/keywords.f 2015-11-17 23:33:35.908986321 +0000 r29152/keywords.f 2015-11-17 23:33:38.389019580 +0000
158:       CHFREQSC=1158:       CHFREQSC=1
159:       CHFREQBB=1159:       CHFREQBB=1
160:       FTRANS=1160:       FTRANS=1
161:       FROT=1161:       FROT=1
162: 162: 
163: !163: !
164: ! QCI parameters164: ! QCI parameters
165: !165: !
166:          CONDATT=.FALSE.166:          CONDATT=.FALSE.
167:          QCIPOTT=.FALSE.167:          QCIPOTT=.FALSE.
168:          QCIPOT2T=.FALSE. 
169:          FREEZETOL=1.0D-3 
170:          INTCONSTRAINTT=.FALSE.168:          INTCONSTRAINTT=.FALSE.
171:          INTCONSTRAINTTOL=0.1D0169:          INTCONSTRAINTTOL=0.1D0
172:          INTCONSTRAINTDEL=10.0D0170:          INTCONSTRAINTDEL=10.0D0
173:          INTCONSTRAINTREP=100.0D0171:          INTCONSTRAINTREP=100.0D0
174:          INTCONSTRAINREPCUT=1.7D0172:          INTCONSTRAINREPCUT=1.7D0
175:          INTFREEZET=.FALSE.173:          INTFREEZET=.FALSE.
176:          INTFREEZETOL=1.0D-3174:          INTFREEZETOL=1.0D-3
177:          INTFREEZEMIN=10175:          INTFREEZEMIN=10
178:          INTCONFRAC=0.9D0176:          INTCONFRAC=0.9D0
179:          INTCONSEP=15177:          INTCONSEP=15
532:       CHNMIN=0.D0530:       CHNMIN=0.D0
533:       CHNMAX=HUGE(1.0D0)531:       CHNMAX=HUGE(1.0D0)
534:       CHMDT=.FALSE.532:       CHMDT=.FALSE.
535:       CHMDFREQ=HUGE(1)533:       CHMDFREQ=HUGE(1)
536:       CURRENTIMP=0534:       CURRENTIMP=0
537:       BOXT=.FALSE.535:       BOXT=.FALSE.
538:       SPHERET=.FALSE.536:       SPHERET=.FALSE.
539:       RMST=.FALSE.537:       RMST=.FALSE.
540:       NEWCONFT=.FALSE.538:       NEWCONFT=.FALSE.
541:       INTMINT=.FALSE.539:       INTMINT=.FALSE.
542:       INTSPRINGACTIVET=.FALSE. 
543:       INTMINFAC=1.0D0 
544:       DAESTAT=.FALSE.540:       DAESTAT=.FALSE.
545:       MAKEOLIGOT=.FALSE.541:       MAKEOLIGOT=.FALSE.
546:       MAKEOLIGOSTART=.FALSE.542:       MAKEOLIGOSTART=.FALSE.
547:       TRANSXYT=.FALSE.543:       TRANSXYT=.FALSE.
548:       ROTZT=.FALSE.544:       ROTZT=.FALSE.
549:       NREPEAT=0545:       NREPEAT=0
550:       NFIXSEG=0546:       NFIXSEG=0
551:       OHCELLT=.FALSE.547:       OHCELLT=.FALSE.
552: 548: 
553: ! khs26> AMBER12 stuff549: ! khs26> AMBER12 stuff
699:       PERMINVOPT=.FALSE.695:       PERMINVOPT=.FALSE.
700: 696: 
701:       GAMMA=1.0D0697:       GAMMA=1.0D0
702:       TUNNELT=.FALSE.698:       TUNNELT=.FALSE.
703:       699:       
704:       TWOD=.FALSE.700:       TWOD=.FALSE.
705:       COMPRESST=.FALSE.701:       COMPRESST=.FALSE.
706: 702: 
707:       MUPDATE=4703:       MUPDATE=4
708:       DGUESS=0.1D0704:       DGUESS=0.1D0
709:       INTDGUESS=0.001D0 
710:       BFGS=.FALSE.705:       BFGS=.FALSE.
711:       LBFGST=.TRUE.706:       LBFGST=.TRUE.
712:       CONJG=.FALSE.707:       CONJG=.FALSE.
713:       TNT=.FALSE.708:       TNT=.FALSE.
714:       TOLB=0.1D0709:       TOLB=0.1D0
715:       DBRENTT=.FALSE.710:       DBRENTT=.FALSE.
716:       GUIDECUT=0.0001D0711:       GUIDECUT=0.0001D0
717:       CPMD=.FALSE.712:       CPMD=.FALSE.
718:       DL_POLY=.FALSE.713:       DL_POLY=.FALSE.
719:       EFAC=0.0D0714:       EFAC=0.0D0
1892:          ELSE IF (NITEMS.EQ.3) THEN1887:          ELSE IF (NITEMS.EQ.3) THEN
1893:             CALL READF(CENTX)1888:             CALL READF(CENTX)
1894:             CALL READF(CENTY)1889:             CALL READF(CENTY)
1895:          ELSE IF (NITEMS.EQ.4) THEN1890:          ELSE IF (NITEMS.EQ.4) THEN
1896:             CALL READF(CENTX)1891:             CALL READF(CENTX)
1897:             CALL READF(CENTY)1892:             CALL READF(CENTY)
1898:             CALL READF(CENTZ)1893:             CALL READF(CENTZ)
1899:          ENDIF1894:          ENDIF
1900:       ELSE IF (WORD.EQ.'QCIPOT') THEN1895:       ELSE IF (WORD.EQ.'QCIPOT') THEN
1901:          QCIPOTT=.TRUE.1896:          QCIPOTT=.TRUE.
1902:       ELSE IF (WORD.EQ.'QCIPOT2') THEN 
1903:          QCIPOTT=.TRUE. 
1904:          QCIPOT2T=.TRUE. 
1905: 1897: 
1906: !       csw34> QUCENTRE moves the centre of coordinates to (0,0,0)1898: !       csw34> QUCENTRE moves the centre of coordinates to (0,0,0)
1907: !       before each step is taken. 1899: !       before each step is taken. 
1908:       ELSE IF (WORD.EQ.'QUCENTRE') THEN1900:       ELSE IF (WORD.EQ.'QUCENTRE') THEN
1909:          QUCENTRE=.TRUE.1901:          QUCENTRE=.TRUE.
1910: !1902: !
1911: !  Conjugate gradient optimisation instead of default LBFGS1903: !  Conjugate gradient optimisation instead of default LBFGS
1912: !1904: !
1913:       ELSE IF (WORD.EQ.'CG') THEN1905:       ELSE IF (WORD.EQ.'CG') THEN
1914:          LBFGST=.FALSE.1906:          LBFGST=.FALSE.
2501:          CHRIGIDROTT=.TRUE.2493:          CHRIGIDROTT=.TRUE.
2502:          CALL READF(PROT)2494:          CALL READF(PROT)
2503:          CALL READF(ROTMAX)2495:          CALL READF(ROTMAX)
2504:          CALL READI(FROT)2496:          CALL READI(FROT)
2505: 2497: 
2506: !2498: !
2507: ! Check for internal minimum in constraint terms for INTCONSTRAINT2499: ! Check for internal minimum in constraint terms for INTCONSTRAINT
2508: !2500: !
2509:          ELSE IF ((WORD.EQ.'CONINT').OR.(WORD.EQ.'QCIINT')) THEN2501:          ELSE IF ((WORD.EQ.'CONINT').OR.(WORD.EQ.'QCIINT')) THEN
2510:             CHECKCONINT=.TRUE.2502:             CHECKCONINT=.TRUE.
2511:             IF (NITEMS.GT.1) CALL READF(INTMINFAC) 
2512:             WRITE(MYUNIT,'(A,G20.10)') ' keyword> Internal minima terms will be scaled by a factor of ',INTMINFAC  
2513: !2503: !
2514: ! Absolute distance to allow before turning on constraint potential.2504: ! Absolute distance to allow before turning on constraint potential.
2515: !2505: !
2516:          ELSE IF (WORD.EQ.'CONCUTABS') THEN2506:          ELSE IF (WORD.EQ.'CONCUTABS') THEN
2517:             CONCUTABST=.TRUE.2507:             CONCUTABST=.TRUE.
2518:             CONCUTFRACT=.FALSE.2508:             CONCUTFRACT=.FALSE.
2519:             IF (NITEMS.GT.1) CALL READF(CONCUTABS)2509:             IF (NITEMS.GT.1) CALL READF(CONCUTABS)
2520: !2510: !
2521: ! Fraction of constraint distance to allow before turning on constraint potential.2511: ! Fraction of constraint distance to allow before turning on constraint potential.
2522: !2512: !
2929: !2919: !
2930: ! Tiffany's DFTB potential for carbon, transplanted from OPTIM.2920: ! Tiffany's DFTB potential for carbon, transplanted from OPTIM.
2931: !2921: !
2932:       ELSE IF (WORD.EQ.'DFTBC') THEN2922:       ELSE IF (WORD.EQ.'DFTBC') THEN
2933:          DFTBCT=.TRUE.2923:          DFTBCT=.TRUE.
2934: !2924: !
2935: !  Initial guess for diagonal matrix elements in LBFGS.2925: !  Initial guess for diagonal matrix elements in LBFGS.
2936: !2926: !
2937:       ELSE IF (WORD.EQ.'DGUESS') THEN2927:       ELSE IF (WORD.EQ.'DGUESS') THEN
2938:          CALL READF(DGUESS)2928:          CALL READF(DGUESS)
2939:          IF (NITEMS.GT.2) CALL READF(INTDGUESS) 
2940: 2929: 
2941: !      ELSE IF (WORD.EQ.'DIELEC') THEN2930: !      ELSE IF (WORD.EQ.'DIELEC') THEN
2942: !         CALL READF(XX)2931: !         CALL READF(XX)
2943: !         DPARAM=XX2932: !         DPARAM=XX
2944: !         WRITE(MYUNIT,'(A,F9.5)') ' Dielectric constant = ',DPARAM2933: !         WRITE(MYUNIT,'(A,F9.5)') ' Dielectric constant = ',DPARAM
2945: !2934: !
2946: !  NOT DOCUMENTED - INTENTIONAL2935: !  NOT DOCUMENTED - INTENTIONAL
2947: !2936: !
2948:       ELSE IF (WORD.EQ.'DIFFRACT') THEN2937:       ELSE IF (WORD.EQ.'DIFFRACT') THEN
2949:          DIFFRACTT=.TRUE.2938:          DIFFRACTT=.TRUE.
3213: !3202: !
3214: !  Frozen atoms.3203: !  Frozen atoms.
3215: !3204: !
3216:       ELSE IF (WORD.EQ.'FREEZE') THEN3205:       ELSE IF (WORD.EQ.'FREEZE') THEN
3217:          FREEZE=.TRUE.3206:          FREEZE=.TRUE.
3218:          DO J1=1,NITEMS-13207:          DO J1=1,NITEMS-1
3219:             NFREEZE=NFREEZE+13208:             NFREEZE=NFREEZE+1
3220:             CALL READI(NDUMMY)3209:             CALL READI(NDUMMY)
3221:             FROZEN(NDUMMY)=.TRUE.3210:             FROZEN(NDUMMY)=.TRUE.
3222:          ENDDO3211:          ENDDO
3223:       ELSE IF (WORD.EQ.'FREEZENODES') THEN 
3224:          FREEZENODEST=.TRUE. 
3225:          CALL READF(FREEZETOL) 
3226: 3212: 
3227: !3213: !
3228: ! sf344> unfreeze everything at the final quenches3214: ! sf344> unfreeze everything at the final quenches
3229: !3215: !
3230:       ELSE IF (WORD.EQ.'UNFREEZEFINALQ') THEN3216:       ELSE IF (WORD.EQ.'UNFREEZEFINALQ') THEN
3231:         UNFREEZEFINALQ=.TRUE.3217:         UNFREEZEFINALQ=.TRUE.
3232: !3218: !
3233: ! csw343219: ! csw34
3234: ! Frozen residues (to be converted to frozen atoms)3220: ! Frozen residues (to be converted to frozen atoms)
3235: !3221: !
3911:             CALL READF(XX)3897:             CALL READF(XX)
3912:             EXPFAC=XX3898:             EXPFAC=XX
3913:          ENDIF3899:          ENDIF
3914:          IF (NITEMS.GT.3) THEN3900:          IF (NITEMS.GT.3) THEN
3915:             CALL READF(XX)3901:             CALL READF(XX)
3916:             EXPD=XX3902:             EXPD=XX
3917:          ENDIF3903:          ENDIF
3918: 3904: 
3919:       ELSE IF (WORD.EQ.'INTMIN') THEN3905:       ELSE IF (WORD.EQ.'INTMIN') THEN
3920:          INTMINT=.TRUE.3906:          INTMINT=.TRUE.
 3907: !        IF (NITEMS.GT.1) THEN
 3908: !           CALL READF(IMINCUT)
 3909: !        ENDIF
 3910: 
3921: !3911: !
3922: ! Images for INTCONSTRAINT3912: ! Images for INTCONSTRAINT
3923: !3913: !
3924:          ELSE IF ((WORD.EQ.'INTIMAGE').OR.(WORD.EQ.'QCIIMAGE')) THEN3914:          ELSE IF ((WORD.EQ.'INTIMAGE').OR.(WORD.EQ.'QCIIMAGE')) THEN
3925:             IF (NITEMS.GT.1) CALL READF(IMSEPMIN)3915:             IF (NITEMS.GT.1) CALL READF(IMSEPMIN)
3926:             IF (NITEMS.GT.2) CALL READF(IMSEPMAX)3916:             IF (NITEMS.GT.2) CALL READF(IMSEPMAX)
3927:             IF (NITEMS.GT.3) CALL READI(INTIMAGE)3917:             IF (NITEMS.GT.3) CALL READI(INTIMAGE)
3928:             IF (NITEMS.GT.4) CALL READI(MAXINTIMAGE)3918:             IF (NITEMS.GT.4) CALL READI(MAXINTIMAGE)
3929:             IF (NITEMS.GT.5) CALL READI(INTNTRIESMAX)3919:             IF (NITEMS.GT.5) CALL READI(INTNTRIESMAX)
3930:             IF (NITEMS.GT.6) CALL READI(INTIMAGEINCR)3920:             IF (NITEMS.GT.6) CALL READI(INTIMAGEINCR)
4061:                ENDDO4051:                ENDDO
4062:                CLOSE(LUNIT)4052:                CLOSE(LUNIT)
4063:                PRINT '(A,I6,A)',' keyword> Read ',NCONGEOM,' reference geometries from congeom file'4053:                PRINT '(A,I6,A)',' keyword> Read ',NCONGEOM,' reference geometries from congeom file'
4064:                IF (NCONGEOM.LT.2) PRINT '(A)',' WARNING *** insufficient reference geometries - using end point minima'4054:                IF (NCONGEOM.LT.2) PRINT '(A)',' WARNING *** insufficient reference geometries - using end point minima'
4065:             ENDIF4055:             ENDIF
4066: !4056: !
4067: ! Use the quasi-continuous metric for connection attempts, instead of distance.4057: ! Use the quasi-continuous metric for connection attempts, instead of distance.
4068: !4058: !
4069:             INTERPCOSTFUNCTION=.TRUE.4059:             INTERPCOSTFUNCTION=.TRUE.
4070: !4060: !
4071: ! Only include spring gradients for active atoms. 
4072: ! 
4073:       ELSE IF (WORD.EQ.'INTSPRINGACTIVE') THEN 
4074:          INTSPRINGACTIVET=.TRUE. 
4075: ! 
4076: ! KINT: force constant for springs in INTCONSTRAINT calculations.4061: ! KINT: force constant for springs in INTCONSTRAINT calculations.
4077: ! Default zero.4062: ! Default zero.
4078: !4063: !
4079:          ELSE IF (WORD.EQ.'KINT') THEN4064:          ELSE IF (WORD.EQ.'KINT') THEN
4080:             CALL READF(KINT)4065:             CALL READF(KINT)
4081:  4066:  
4082:       ELSE IF (WORD.EQ.'JM') THEN4067:       ELSE IF (WORD.EQ.'JM') THEN
4083:          JMT=.TRUE.4068:          JMT=.TRUE.
4084: 4069: 
4085:       ELSE IF (WORD.EQ.'JUMPMOVE') THEN4070:       ELSE IF (WORD.EQ.'JUMPMOVE') THEN


r29153/main.F 2015-11-17 23:33:36.100988900 +0000 r29152/main.F 2015-11-17 23:33:38.581022157 +0000
 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS
 35:       USE GAUSS_MOD, ONLY: KEGEN 35:       USE GAUSS_MOD, ONLY: KEGEN
 36:       IMPLICIT NONE 36:       IMPLICIT NONE
 37:       !EXTERNAL READ_CMD_ARGS 37:       !EXTERNAL READ_CMD_ARGS
 38: #ifdef MPI 38: #ifdef MPI
 39:       INCLUDE 'mpif.h' 39:       INCLUDE 'mpif.h'
 40: #endif 40: #endif
 41:       INTEGER J1,J2, JP, MPIERR, NDUMMY3,NPTOTAL,VERSIONTEMP,GETUNIT,LUNIT 41:       INTEGER J1,J2, JP, MPIERR, NDUMMY3,NPTOTAL,VERSIONTEMP,GETUNIT,LUNIT
 42:       DOUBLE PRECISION, ALLOCATABLE :: SCREENC(:) 42:       DOUBLE PRECISION, ALLOCATABLE :: SCREENC(:)
 43:       DOUBLE PRECISION POTEL, DUMMY, INTFREEZETOLSAVE, RMAT(3,3), DUMMY2, DIST2, LINTCONSTRAINTTOL 43:       DOUBLE PRECISION POTEL, DUMMY, INTFREEZETOLSAVE, RMAT(3,3), DUMMY2, DIST2, LINTCONSTRAINTTOL
 44:       DOUBLE PRECISION, ALLOCATABLE :: TMPCOORDS(:), ENDCOORDS(:,:) 44:       DOUBLE PRECISION, ALLOCATABLE :: TMPCOORDS(:)
 45:       INTEGER, ALLOCATABLE :: NDUMMY(:), NDUMMY2(:,:) 45:       INTEGER, ALLOCATABLE :: NDUMMY(:), NDUMMY2(:,:)
 46:       LOGICAL LOPEN, YESNO 46:       LOGICAL LOPEN, YESNO
 47:       INTEGER J6 47:       INTEGER J6
 48:  48: 
 49:       CHARACTER(LEN=130) ISTR,JSTR 49:       CHARACTER(LEN=130) ISTR,JSTR
 50:       CHARACTER(LEN=40) :: atom1,atom2,atompair 50:       CHARACTER(LEN=40) :: atom1,atom2,atompair
 51:  51: 
 52:       CHARACTER(LEN=13) :: CUDAFILENAME1 52:       CHARACTER(LEN=13) :: CUDAFILENAME1
 53:       CHARACTER(LEN=19) :: CUDAFILENAME2 53:       CHARACTER(LEN=19) :: CUDAFILENAME2
 54:       CHARACTER(LEN=21) :: CUDAFILENAME3 54:       CHARACTER(LEN=21) :: CUDAFILENAME3
135: 135: 
136:       ALLOCATE(FINISH(3*NATOMSALLOC))136:       ALLOCATE(FINISH(3*NATOMSALLOC))
137: 137: 
138:       INQUIRE(UNIT=1,OPENED=LOPEN)138:       INQUIRE(UNIT=1,OPENED=LOPEN)
139:       IF (LOPEN) THEN139:       IF (LOPEN) THEN
140:          WRITE(*,'(A,I2,A)') 'main> A ERROR *** Unit ', 1, ' is not free '140:          WRITE(*,'(A,I2,A)') 'main> A ERROR *** Unit ', 1, ' is not free '
141:          STOP141:          STOP
142:       ENDIF142:       ENDIF
143: 143: 
144:       CALL KEYWORD144:       CALL KEYWORD
 145:       IF (QCIPOTT) THEN
 146:          ALLOCATE(TMPCOORDS(3*NATOMSALLOC))
 147:          TMPCOORDS(1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1)
 148:          IF (ALLOCATED(INTFROZEN)) DEALLOCATE(INTFROZEN)
 149:          ALLOCATE(INTFROZEN(NATOMS))
 150:          INTFROZEN(1:NATOMS)=.FALSE.
 151:          CALL CHECKPERC(TMPCOORDS,LINTCONSTRAINTTOL,0,1)
 152:          DEALLOCATE(TMPCOORDS)
 153:       ENDIF
145: !     write out the atom indices to check:154: !     write out the atom indices to check:
146: !        IF(TWISTT) THEN155: !        IF(TWISTT) THEN
147: !              WRITE(MYUNIT,*)'NWISTGROUPS dihedrals to be constrained..'156: !              WRITE(MYUNIT,*)'NWISTGROUPS dihedrals to be constrained..'
148: !              DO J1=1, NTWISTGROUPS157: !              DO J1=1, NTWISTGROUPS
149: !                WRITE(MYUNIT,*) 'Group J1  indices:'158: !                WRITE(MYUNIT,*) 'Group J1  indices:'
150: !                WRITE(MYUNIT,*) TWIST_ATOMS(1:4, J1)159: !                WRITE(MYUNIT,*) TWIST_ATOMS(1:4, J1)
151: !                WRITE(MYUNIT,*) 'Group J1  k value, ref.dihedral:'160: !                WRITE(MYUNIT,*) 'Group J1  k value, ref.dihedral:'
152: !                WRITE(MYUNIT,*) TWIST_K(J1), TWIST_THETA0(J1)161: !                WRITE(MYUNIT,*) TWIST_K(J1), TWIST_THETA0(J1)
153: !              END DO162: !              END DO
154: !        END IF163: !        END IF
298: 307: 
299:       ALLOCATE(FF(NSAVE),QMIN(MAX(NSAVE,1)))308:       ALLOCATE(FF(NSAVE),QMIN(MAX(NSAVE,1)))
300:       ALLOCATE(QMINP(NSAVE,3*NATOMSALLOC))309:       ALLOCATE(QMINP(NSAVE,3*NATOMSALLOC))
301:       ALLOCATE(QMINNATOMS(NSAVE))310:       ALLOCATE(QMINNATOMS(NSAVE))
302:       ALLOCATE(QMINT(NSAVE,NATOMSALLOC))311:       ALLOCATE(QMINT(NSAVE,NATOMSALLOC))
303:       ALLOCATE(NPCALL_QMIN(NSAVE))312:       ALLOCATE(NPCALL_QMIN(NSAVE))
304:       IF (MONITORT) THEN313:       IF (MONITORT) THEN
305:          ALLOCATE(LOWESTC(3*NATOMSALLOC))314:          ALLOCATE(LOWESTC(3*NATOMSALLOC))
306:       ENDIF315:       ENDIF
307: 316: 
 317: !
 318: ! Allocations for QCI
 319: !
 320:       IF (INTCONSTRAINTT.AND.(NCONGEOM.GE.2)) THEN
 321: !
 322: ! Set up all the constraints and repulsions for zero frozen atoms.
 323: !
 324:          IF (.NOT.ALLOCATED(CONI)) THEN
 325:             ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX))
 326:             ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
 327:          ENDIF
 328:          CALL MINPERMDIST(COORDS,FINISH,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY2,DIST2,RIGID,RMAT)
 329:          INTFREEZETOLSAVE=INTFREEZETOL
 330:          INTFREEZETOL=-1.0D0
 331:          CALL MAKE_CONPOT(NCONGEOM,CONGEOM)
 332:          INTFREEZETOL=INTFREEZETOLSAVE
 333:       ENDIF
308: !        csw34> ALLOCATE the interaction energy tracking arrays if A9INTE in data334: !        csw34> ALLOCATE the interaction energy tracking arrays if A9INTE in data
309:       IF (A9INTET.AND.AMBERT) THEN 335:       IF (A9INTET.AND.AMBERT) THEN 
310:          ALLOCATE(INTEFF(NSAVEINTE),INTEQMIN(NSAVEINTE))336:          ALLOCATE(INTEFF(NSAVEINTE),INTEQMIN(NSAVEINTE))
311:          ALLOCATE(INTEQMINP(NSAVEINTE,3*NATOMS))337:          ALLOCATE(INTEQMINP(NSAVEINTE,3*NATOMS))
312:          INTEQMIN(:)=1.0D10 338:          INTEQMIN(:)=1.0D10 
313:          INTEQMINP(1:NSAVEINTE,1:3*NATOMS)=0.0D0 ! to prevent reading from uninitialised memory339:          INTEQMINP(1:NSAVEINTE,1:3*NATOMS)=0.0D0 ! to prevent reading from uninitialised memory
314:          INTEFF(1:NSAVEINTE)=0 ! to prevent reading from uninitialised memory340:          INTEFF(1:NSAVEINTE)=0 ! to prevent reading from uninitialised memory
315:       ENDIF  341:       ENDIF  
316: 342: 
317:       IF (GAUSST) THEN343:       IF (GAUSST) THEN
412:                IF (RMST) OPEN(MYRUNIT,FILE='rmsd',STATUS='UNKNOWN',FORM='FORMATTED')438:                IF (RMST) OPEN(MYRUNIT,FILE='rmsd',STATUS='UNKNOWN',FORM='FORMATTED')
413:             ENDIF439:             ENDIF
414:             IF (A9INTET) THEN440:             IF (A9INTET) THEN
415:                OPEN(UNIT=3998,FILE='intE.dat',STATUS='UNKNOWN',FORM='FORMATTED')441:                OPEN(UNIT=3998,FILE='intE.dat',STATUS='UNKNOWN',FORM='FORMATTED')
416:                OPEN(UNIT=3999,FILE='bestintE.dat',STATUS='UNKNOWN',FORM='FORMATTED')442:                OPEN(UNIT=3999,FILE='bestintE.dat',STATUS='UNKNOWN',FORM='FORMATTED')
417:             ENDIF443:             ENDIF
418:          ENDIF444:          ENDIF
419:       ENDIF445:       ENDIF
420:       CALL FLUSH(6)446:       CALL FLUSH(6)
421:       CALL IO1447:       CALL IO1
422:       IF (QCIPOTT) THEN 
423:          ALLOCATE(TMPCOORDS(3*NATOMSALLOC)) 
424:          TMPCOORDS(1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1) 
425:          IF (ALLOCATED(INTFROZEN)) DEALLOCATE(INTFROZEN) 
426:          ALLOCATE(INTFROZEN(NATOMS)) 
427:          INTFROZEN(1:NATOMS)=.FALSE. 
428: !        WRITE(MYUNIT, *) NATOMSALLOC 
429: !        WRITE(MYUNIT, *) 'A tmpcoords before' 
430: !        DO J1=1,NATOMSALLOC 
431: !           WRITE(MYUNIT,'(A,3G20.10)') 'O ',TMPCOORDS(3*(J1-1)+1:3*(J1-1)+3) 
432: !        ENDDO 
433:          CALL CHECKPERC(TMPCOORDS,LINTCONSTRAINTTOL,0,1) 
434: !        WRITE(MYUNIT, *) NATOMSALLOC 
435: !        WRITE(MYUNIT, *) 'A tmpcoords after' 
436: !        DO J1=1,NATOMSALLOC 
437: !           WRITE(MYUNIT,'(A,3G20.10)') 'O ',TMPCOORDS(3*(J1-1)+1:3*(J1-1)+3) 
438: !        ENDDO 
439:          DEALLOCATE(TMPCOORDS) 
440:       ENDIF 
441: ! 
442: ! Allocations for QCI 
443: ! 
444:       IF (INTCONSTRAINTT.AND.(NCONGEOM.GE.2)) THEN 
445: ! 
446: ! Set up all the constraints and repulsions for zero frozen atoms. 
447: ! 
448:          IF (.NOT.ALLOCATED(CONI)) THEN 
449:             ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX)) 
450:             ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX)) 
451:          ENDIF 
452:  
453:          ALLOCATE(TMPCOORDS(3*NATOMSALLOC)) 
454:          TMPCOORDS(1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1) 
455:          DEALLOCATE(TMPCOORDS) 
456:          INTFREEZETOLSAVE=INTFREEZETOL 
457:          INTFREEZETOL=-1.0D0 
458:          CALL MAKE_CONPOT(NCONGEOM,CONGEOM) 
459:          INTFREEZETOL=INTFREEZETOLSAVE 
460: ! 
461: ! Now align the two endpoints to the first reference minimum. 
462: ! 
463:           
464:          ALLOCATE(ENDCOORDS(2,3*NATOMSALLOC)) 
465:          ENDCOORDS(1,1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1) 
466:          ENDCOORDS(2,1:3*NATOMSALLOC)=FINISH(1:3*NATOMSALLOC) 
467:          CALL MAKE_CONPOT(2,ENDCOORDS) 
468:          COORDS(1:3*NATOMSALLOC,1)=ENDCOORDS(1,1:3*NATOMSALLOC) 
469:          FINISH(1:3*NATOMSALLOC)=ENDCOORDS(2,1:3*NATOMSALLOC) 
470:          DEALLOCATE(ENDCOORDS) 
471:       ENDIF 
472: !448: !
473: ! If this is a CSM optimisation we now have to multiply the number of atoms by the number of449: ! If this is a CSM optimisation we now have to multiply the number of atoms by the number of
474: ! group operations and replicate some coordinates and allowed permutations.450: ! group operations and replicate some coordinates and allowed permutations.
475: !451: !
476:       IF (CSMT) THEN452:       IF (CSMT) THEN
477:          CALL CSMINIT453:          CALL CSMINIT
478:          IF (SYMMETRIZECSM) THEN454:          IF (SYMMETRIZECSM) THEN
479:             IF (CSMMAXIT.EQ.0) CSMMAXIT=MAXIT455:             IF (CSMMAXIT.EQ.0) CSMMAXIT=MAXIT
480:          ELSE456:          ELSE
481:             CSMMAXIT=MAXIT457:             CSMMAXIT=MAXIT


r29153/make_conpot.f90 2015-11-17 23:33:36.292991475 +0000 r29152/make_conpot.f90 2015-11-17 23:33:38.765024629 +0000
 64: ALLOCATE(INTFROZEN(NATOMS)) 64: ALLOCATE(INTFROZEN(NATOMS))
 65: INTFROZEN(1:NATOMS)=.FALSE. 65: INTFROZEN(1:NATOMS)=.FALSE.
 66: DLIST(1:NATOMS)=-1 66: DLIST(1:NATOMS)=-1
 67: DMOVED(1:NATOMS)=1.0D100 67: DMOVED(1:NATOMS)=1.0D100
 68: IF (INTFREEZET) THEN 68: IF (INTFREEZET) THEN
 69:    IF (NCPFIT.GT.1) THEN 69:    IF (NCPFIT.GT.1) THEN
 70:       DO J1=1,NATOMS 70:       DO J1=1,NATOMS
 71:          DF=SQRT((MINCOORDS(1,3*(J1-1)+1)-MINCOORDS(2,3*(J1-1)+1))**2 & 71:          DF=SQRT((MINCOORDS(1,3*(J1-1)+1)-MINCOORDS(2,3*(J1-1)+1))**2 &
 72:   &             +(MINCOORDS(1,3*(J1-1)+2)-MINCOORDS(2,3*(J1-1)+2))**2 & 72:   &             +(MINCOORDS(1,3*(J1-1)+2)-MINCOORDS(2,3*(J1-1)+2))**2 &
 73:   &             +(MINCOORDS(1,3*(J1-1)+3)-MINCOORDS(2,3*(J1-1)+3))**2) 73:   &             +(MINCOORDS(1,3*(J1-1)+3)-MINCOORDS(2,3*(J1-1)+3))**2)
 74:          IF (J1.EQ.NATOMS) THEN 
 75:             WRITE(MYUNIT,'(A,6G20.10)') 'mincoords atom 400: ',MINCOORDS(1,1198:1200),MINCOORDS(2,1198:1200) 
 76:             WRITE(MYUNIT,'(A,6G20.10)') 'DF,INTFREEZETOL=',DF,INTFREEZETOL 
 77:          ENDIF 
 78:  
 79:          IF (DF.LT.INTFREEZETOL) THEN 74:          IF (DF.LT.INTFREEZETOL) THEN
 80:             NQCIFREEZE=NQCIFREEZE+1 75:             NQCIFREEZE=NQCIFREEZE+1
 81:             INTFROZEN(J1)=.TRUE. 76:             INTFROZEN(J1)=.TRUE.
 82:             IF (DEBUG) WRITE(MYUNIT, '(A,I6,A,F12.6,A,I6)') ' make_conpot> atom ',J1, & 77: !           IF (DEBUG) WRITE(MYUNIT, '(A,I6,A,F12.6,A,I6)') ' make_conpot> atom ',J1, &
 83:   &                          ' moves less than threshold: distance=',DF,' total=',NQCIFREEZE 78: ! &                          ' moves less than threshold: distance=',DF,' total=',NQCIFREEZE
 84:          ENDIF 79:          ENDIF
 85:          sortd: DO J2=1,J1 80:          sortd: DO J2=1,J1
 86:             IF (DF.LT.DMOVED(J2)) THEN 81:             IF (DF.LT.DMOVED(J2)) THEN
 87:                DO J3=J1,J2+1,-1 82:                DO J3=J1,J2+1,-1
 88:                   DMOVED(J3)=DMOVED(J3-1) 83:                   DMOVED(J3)=DMOVED(J3-1)
 89:                   DLIST(J3)=DLIST(J3-1) 84:                   DLIST(J3)=DLIST(J3-1)
 90:                ENDDO 85:                ENDDO
 91:                DMOVED(J2)=DF 86:                DMOVED(J2)=DF
 92:                DLIST(J2)=J1 87:                DLIST(J2)=J1
 93:                EXIT sortd 88:                EXIT sortd
137: 132: 
138: !133: !
139: ! Fixed repulsions based on congeom file entries134: ! Fixed repulsions based on congeom file entries
140: ! Just need to adjust the list based on any frozen atoms and check135: ! Just need to adjust the list based on any frozen atoms and check
141: ! to make sure a new pair of minima don't have repulsive atoms within136: ! to make sure a new pair of minima don't have repulsive atoms within
142: ! the current cutoff.137: ! the current cutoff.
143: ! 138: ! 
144: IF (NCONGEOM.GE.2) THEN139: IF (NCONGEOM.GE.2) THEN
145:    IF (CALLED.OR.CONDATT) THEN140:    IF (CALLED.OR.CONDATT) THEN
146:       J2=0141:       J2=0
147:       WRITE(MYUNIT,'(A,I8)') 'NCPFIT=',NCPFIT 
148:       DO J1=1,NREPULSIVEFIX142:       DO J1=1,NREPULSIVEFIX
149: !143: !
150: ! If called with two minima check that REPCUTFIX doesn't exceed the separation in 144: ! If called with two minima check that REPCUTFIX doesn't exceed the separation in 
151: ! either minimum.145: ! either minimum.
152: !146: !
153:           
154:          IF (NCPFIT.EQ.2) THEN147:          IF (NCPFIT.EQ.2) THEN
155:             DF=MIN(SQRT((MINCOORDS(1,3*(REPIFIX(J1)-1)+1)-MINCOORDS(1,3*(REPJFIX(J1)-1)+1))**2+ &148:             DF=MIN(SQRT((MINCOORDS(1,3*(REPIFIX(J1)-1)+1)-MINCOORDS(1,3*(REPJFIX(J1)-1)+1))**2+ &
156:   &                     (MINCOORDS(1,3*(REPIFIX(J1)-1)+2)-MINCOORDS(1,3*(REPJFIX(J1)-1)+2))**2+ &149:   &                     (MINCOORDS(1,3*(REPIFIX(J1)-1)+2)-MINCOORDS(1,3*(REPJFIX(J1)-1)+2))**2+ &
157:   &                     (MINCOORDS(1,3*(REPIFIX(J1)-1)+3)-MINCOORDS(1,3*(REPJFIX(J1)-1)+3))**2),&150:   &                     (MINCOORDS(1,3*(REPIFIX(J1)-1)+3)-MINCOORDS(1,3*(REPJFIX(J1)-1)+3))**2),&
158:                    SQRT((MINCOORDS(2,3*(REPIFIX(J1)-1)+1)-MINCOORDS(2,3*(REPJFIX(J1)-1)+1))**2+ &151:                    SQRT((MINCOORDS(2,3*(REPIFIX(J1)-1)+1)-MINCOORDS(2,3*(REPJFIX(J1)-1)+1))**2+ &
159:   &                     (MINCOORDS(2,3*(REPIFIX(J1)-1)+2)-MINCOORDS(2,3*(REPJFIX(J1)-1)+2))**2+ &152:   &                     (MINCOORDS(2,3*(REPIFIX(J1)-1)+2)-MINCOORDS(2,3*(REPJFIX(J1)-1)+2))**2+ &
160:   &                     (MINCOORDS(2,3*(REPIFIX(J1)-1)+3)-MINCOORDS(2,3*(REPJFIX(J1)-1)+3))**2))153:   &                     (MINCOORDS(2,3*(REPIFIX(J1)-1)+3)-MINCOORDS(2,3*(REPJFIX(J1)-1)+3))**2))
161: !           WRITE(MYUNIT,'(A,I8,2G20.10)') 'J1,DF,REPCUTFIX(J1)=',J1,DF,REPCUTFIX(J1) 
162:             IF (DF.LT.REPCUTFIX(J1)) THEN154:             IF (DF.LT.REPCUTFIX(J1)) THEN
163:                WRITE(MYUNIT, '(A,2I6,2(A,G15.5))') ' make_conpot> Reducing repulsive cutoff for atoms ', &155:                WRITE(MYUNIT, '(A,2I6,2(A,G15.5))') ' make_conpot> Reducing repulsive cutoff for atoms ', &
164:   &                       REPIFIX(J1),REPJFIX(J1),' from ',REPCUTFIX(J1),' to ',DF-1.0D-3156:   &                       REPIFIX(J1),REPJFIX(J1),' from ',REPCUTFIX(J1),' to ',DF-1.0D-3
165:                REPCUTFIX(J1)=DF-1.0D-3157:                REPCUTFIX(J1)=DF-1.0D-3
166:             ENDIF158:             ENDIF
167:          ENDIF159:          ENDIF
168:          IF (INTFROZEN(REPIFIX(J1)).AND.INTFROZEN(REPJFIX(J1))) CYCLE160:          IF (INTFROZEN(REPIFIX(J1)).AND.INTFROZEN(REPJFIX(J1))) CYCLE
169:          J2=J2+1161:          J2=J2+1
170:          REPI(J2)=REPIFIX(J1)162:          REPI(J2)=REPIFIX(J1)
171:          REPJ(J2)=REPJFIX(J1)163:          REPJ(J2)=REPJFIX(J1)
172:          REPCUT(J2)=REPCUTFIX(J1)164:          REPCUT(J2)=REPCUTFIX(J1)
173:       ENDDO165:       ENDDO
174:       NREPULSIVE=J2166:       NREPULSIVE=J2
175:       WRITE(MYUNIT, '(A,I6,A)') ' make_conpot> After allowing for frozen atoms there are ',NREPULSIVE,' possible repulsions'167: !     WRITE(MYUNIT, '(A,I6,A)') ' make_conpot> After allowing for frozen atoms there are ',NREPULSIVE,' possible repulsions'
176: !     STOP !!! DJW 
177:       NREPI(1:NREPULSIVE)=REPI(1:NREPULSIVE)168:       NREPI(1:NREPULSIVE)=REPI(1:NREPULSIVE)
178:       NREPJ(1:NREPULSIVE)=REPJ(1:NREPULSIVE)169:       NREPJ(1:NREPULSIVE)=REPJ(1:NREPULSIVE)
179:       NNREPULSIVE=NREPULSIVE170:       NNREPULSIVE=NREPULSIVE
180:       NREPCUT(1:NREPULSIVE)=REPCUT(1:NREPULSIVE)171:       NREPCUT(1:NREPULSIVE)=REPCUT(1:NREPULSIVE)
181:       IF (ALLOCATED(CONACTIVE)) DEALLOCATE(CONACTIVE)172:       IF (ALLOCATED(CONACTIVE)) DEALLOCATE(CONACTIVE)
182:       ALLOCATE(CONACTIVE(NCONSTRAINT))173:       ALLOCATE(CONACTIVE(NCONSTRAINT))
183:       CONACTIVE(1:NCONSTRAINT)=.TRUE. 174:       CONACTIVE(1:NCONSTRAINT)=.TRUE. 
184:       IF (ALLOCATED(CONDISTREFLOCAL)) DEALLOCATE(CONDISTREFLOCAL)175:       IF (ALLOCATED(CONDISTREFLOCAL)) DEALLOCATE(CONDISTREFLOCAL)
185:       ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))176:       ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))
186:       CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)177:       CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)


r29153/mc.F 2015-11-17 23:33:36.484994050 +0000 r29152/mc.F 2015-11-17 23:33:38.957027202 +0000
2031: 2031: 
2032:             TEMP(JP)=TEMP(JP)*SCALEFAC2032:             TEMP(JP)=TEMP(JP)*SCALEFAC
2033: #ifdef MPI2033: #ifdef MPI
2034: #else2034: #else
2035:             IF (HIT) GOTO 372035:             IF (HIT) GOTO 37
2036:             IF (EPREV(JP)<EBEST(JP)) EBEST(JP)=EPREV(JP)2036:             IF (EPREV(JP)<EBEST(JP)) EBEST(JP)=EPREV(JP)
2037:             IF ((REPMATCHT).AND.(ABS(MAXVAL(EBEST(:))-MINVAL(EBEST(:)))<ECONV)) GOTO 372037:             IF ((REPMATCHT).AND.(ABS(MAXVAL(EBEST(:))-MINVAL(EBEST(:)))<ECONV)) GOTO 37
2038: #endif2038: #endif
2039:             IF (DUMPINT.GT.0) THEN2039:             IF (DUMPINT.GT.0) THEN
2040:                IF (MOD(J1,DUMPINT).EQ.0) THEN2040:                IF (MOD(J1,DUMPINT).EQ.0) THEN
2041:                   CALL DUMPSTATE(NQ(JP),EBEST,BESTCOORDS,JBEST,JP)2041:                   CALL DUMPSTATE(J1,EBEST,BESTCOORDS,JBEST,JP)
2042:                ENDIF2042:                ENDIF
2043:             ENDIF2043:             ENDIF
2044: !2044: !
2045: ! Accept/reject step over grand canonical relxation blocks.2045: ! Accept/reject step over grand canonical relxation blocks.
2046: !2046: !
2047:             IF (GCBHT.AND.(QGCBH.EQ.GCRELAX)) THEN2047:             IF (GCBHT.AND.(QGCBH.EQ.GCRELAX)) THEN
2048: !2048: !
2049: ! did the GC potential decrease overall in the GC relax block?2049: ! did the GC potential decrease overall in the GC relax block?
2050: ! this is accept/reject for potential compared with previous best value in previous relaxation block2050: ! this is accept/reject for potential compared with previous best value in previous relaxation block
2051: ! (not the value in the Markov chain, the lowest value in the previous block)2051: ! (not the value in the Markov chain, the lowest value in the previous block)


r29153/minpermdist.f90 2015-11-17 23:33:36.676996622 +0000 r29152/minpermdist.f90 2015-11-17 23:33:39.213030635 +0000
100:    ENDIF100:    ENDIF
101:    RETURN101:    RETURN
102: ENDIF102: ENDIF
103: 103: 
104: IF (.NOT.PERMDIST) THEN104: IF (.NOT.PERMDIST) THEN
105:    IF (RIGID) THEN105:    IF (RIGID) THEN
106:       CALL RBMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,QBEST,DEBUG)106:       CALL RBMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,QBEST,DEBUG)
107:       CALL QROTMAT(QBEST,RMATBEST)107:       CALL QROTMAT(QBEST,RMATBEST)
108:    ELSE108:    ELSE
109:       CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGID,DEBUG,RMAT)109:       CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGID,DEBUG,RMAT)
110:       RMATBEST = RMAT110:        RMATBEST = RMAT
111:    ENDIF111:    ENDIF
112:    RETURN112:    RETURN
113: ENDIF113: ENDIF
114: 114: 
115: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!115: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116: ! CALL OCHARMM(DUMMYA,VNEW,ENERGY,.FALSE.,.FALSE.)116: ! CALL OCHARMM(DUMMYA,VNEW,ENERGY,.FALSE.,.FALSE.)
117: ! CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)117: ! CALL POTENTIAL(DUMMYA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.)
118: ! PRINT '(2(A,F25.15))',' Initial energy=',ENERGY,' RMS=',RMS118: ! PRINT '(2(A,F25.15))',' Initial energy=',ENERGY,' RMS=',RMS
119: ! PRINT '(2(A,F25.15))',' for coordinates:'119: ! PRINT '(2(A,F25.15))',' for coordinates:'
120: ! PRINT '(3F25.15)',DUMMYA(1:3*NATOMS)120: ! PRINT '(3F25.15)',DUMMYA(1:3*NATOMS)


r29153/newmindist.f90 2015-11-17 23:33:36.860999094 +0000 r29152/newmindist.f90 2015-11-17 23:33:39.433033584 +0000
294:       RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT294:       RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT
295:    ENDDO295:    ENDDO
296:    IF (.NOT.BULKT) THEN296:    IF (.NOT.BULKT) THEN
297:       IF (STOCKT) THEN297:       IF (STOCKT) THEN
298:          CALL NEWROTGEOMSTOCK(NATOMS,RB,RMAT,CMXA,CMYA,CMZA)298:          CALL NEWROTGEOMSTOCK(NATOMS,RB,RMAT,CMXA,CMYA,CMZA)
299:       ELSE299:       ELSE
300:          CALL NEWROTGEOM(NSIZE,RB,RMAT,CMXA,CMYA,CMZA)300:          CALL NEWROTGEOM(NSIZE,RB,RMAT,CMXA,CMYA,CMZA)
301:       ENDIF301:       ENDIF
302:    ENDIF302:    ENDIF
303: ENDIF303: ENDIF
304: ! WRITE(MYUNIT,'(A,6G20.10)') 'CMA,CMB=',CMXA,CMYA,CMZA,CMXB,CMYB,CMZB 
305: 304: 
306: DEALLOCATE(XA,XB)305: DEALLOCATE(XA,XB)
307: 306: 
308: END SUBROUTINE NEWMINDIST307: END SUBROUTINE NEWMINDIST
309: 308: 
310: SUBROUTINE NEWROTGEOM(NATOMS,COORDS,ROTMAT,CX,CY,CZ)309: SUBROUTINE NEWROTGEOM(NATOMS,COORDS,ROTMAT,CX,CY,CZ)
311: IMPLICIT NONE310: IMPLICIT NONE
312: INTEGER I, J, K, NATOMS311: INTEGER I, J, K, NATOMS
313: DOUBLE PRECISION COORDS(*), R1, R0(3), ROTMAT(3,3), CX, CY, CZ312: DOUBLE PRECISION COORDS(*), R1, R0(3), ROTMAT(3,3), CX, CY, CZ
314: 313: 


r29153/potential.f90 2015-11-17 23:33:37.049001615 +0000 r29152/potential.f90 2015-11-17 23:33:39.621036105 +0000
253: !     END DO253: !     END DO
254: !254: !
255: 255: 
256:    ELSE IF (MORSET) THEN256:    ELSE IF (MORSET) THEN
257: !      CALL RAD(X, GRAD, EREAL, GRADT)257: !      CALL RAD(X, GRAD, EREAL, GRADT)
258: !      CALL MORSE(X, GRAD, EREAL, GRADT)258: !      CALL MORSE(X, GRAD, EREAL, GRADT)
259:       CALL MORSEGH(X, GRAD, EREAL, GRADT, SECT)259:       CALL MORSEGH(X, GRAD, EREAL, GRADT, SECT)
260: 260: 
261:    ELSE IF (QCIPOTT) THEN261:    ELSE IF (QCIPOTT) THEN
262: !     IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREPQCIPOT(X,3*NATOMS,0,1)262: !     IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREPQCIPOT(X,3*NATOMS,0,1)
263:       IF (QCIPOT2T) THEN263:       CALL QCIPOT(EREAL,X,GRAD)
264:          CALL QCIPOT2(EREAL,X,GRAD) 
265:       ELSE 
266:          CALL QCIPOT(EREAL,X,GRAD) 
267:       ENDIF 
268: 264: 
269:    ELSE IF (TOSI) THEN265:    ELSE IF (TOSI) THEN
270:       CALL RAD(X, GRAD, EREAL, GRADT)266:       CALL RAD(X, GRAD, EREAL, GRADT)
271:       IF (EVAPREJECT) RETURN267:       IF (EVAPREJECT) RETURN
272:       CALL TOSIFUMI(X, GRAD, EREAL, GRADT, SECT)268:       CALL TOSIFUMI(X, GRAD, EREAL, GRADT, SECT)
273: 269: 
274:    ELSE IF (WELCH) THEN270:    ELSE IF (WELCH) THEN
275:       CALL RAD(X, GRAD, EREAL, GRADT)271:       CALL RAD(X, GRAD, EREAL, GRADT)
276:       CALL WEL(X, GRAD, EREAL, GRADT, SECT)272:       CALL WEL(X, GRAD, EREAL, GRADT, SECT)
277: 273: 


r29153/qcipot.f90 2015-11-17 23:33:37.241004189 +0000 r29152/qcipot.f90 2015-11-17 23:33:39.817038740 +0000
 30: DOUBLE PRECISION D2,R2AX,R2AY,R2AZ,R2BX,R2BY,R2BZ 30: DOUBLE PRECISION D2,R2AX,R2AY,R2AZ,R2BX,R2BY,R2BZ
 31: DOUBLE PRECISION G2(3) 31: DOUBLE PRECISION G2(3)
 32: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12 32: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12
 33: DOUBLE PRECISION XYZ(3*NATOMS), GGG(3*NATOMS), CCLOCAL 33: DOUBLE PRECISION XYZ(3*NATOMS), GGG(3*NATOMS), CCLOCAL
 34:  34: 
 35: GGG(1:3*NATOMS)=0.0D0 35: GGG(1:3*NATOMS)=0.0D0
 36: ECON=0.0D0; EREP=0.0D0 36: ECON=0.0D0; EREP=0.0D0
 37: ! 37: !
 38: !  Constraint energy and forces. 38: !  Constraint energy and forces.
 39: ! 39: !
 40: ! WRITE(MYUNIT,'(A,I8)') 'qcipot> NCONSTRAINT=',NCONSTRAINT 40: WRITE(MYUNIT,'(A,I8)') 'qcipot> NCONSTRAINT=',NCONSTRAINT
 41: DO J2=1,NCONSTRAINT 41: DO J2=1,NCONSTRAINT
 42: !  WRITE(MYUNIT,'(A,I8,L5)') 'qcipot> J2,CONACTIVE=',J2,CONACTIVE(J2) 42: !  WRITE(MYUNIT,'(A,I8,L5)') 'qcipot> J2,CONACTIVE=',J2,CONACTIVE(J2)
 43:    IF (.NOT.CONACTIVE(J2)) CYCLE 43:    IF (.NOT.CONACTIVE(J2)) CYCLE
 44:    CCLOCAL=CONCUTLOCAL(J2) 44:    CCLOCAL=CONCUTLOCAL(J2)
 45:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 45:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
 46:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 46:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
 47:    NI1=3*(CONI(J2)-1) 47:    NI1=3*(CONI(J2)-1)
 48:    NJ1=3*(CONJ(J2)-1) 48:    NJ1=3*(CONJ(J2)-1)
 49:    R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3) 49:    R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3)
 50:    R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3) 50:    R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3)
127:       NREPCUT(NNREPULSIVE)=REPCUT(JJ)127:       NREPCUT(NNREPULSIVE)=REPCUT(JJ)
128:       GOTO 246128:       GOTO 246
129:    ENDIF129:    ENDIF
130:    COMPARE=CHECKREPCUTOFF*REPCUT(JJ)130:    COMPARE=CHECKREPCUTOFF*REPCUT(JJ)
131: 246 CONTINUE131: 246 CONTINUE
132: ENDDO132: ENDDO
133: IF (DEBUG) WRITE(MYUNIT,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE133: IF (DEBUG) WRITE(MYUNIT,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE
134: 134: 
135: END SUBROUTINE CHECKREPQCIPOT135: END SUBROUTINE CHECKREPQCIPOT
136: 136: 
137: ! 
138: ! QCI potential with no flat region for constraints. 
139: ! 
140: SUBROUTINE QCIPOT2(ECON,XYZ,GGG) 
141: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, & 
142:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, INTCONSTRAINTREP, CONDISTREFLOCAL, & 
143:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT, IMSEPMAX, & 
144:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, & 
145:   &  NATOMS, DEBUG, MYUNIT 
146: USE PORFUNCS 
147: IMPLICIT NONE 
148:             
149: INTEGER :: J2,NI2,NI1,NJ2,NJ1 
150: DOUBLE PRECISION :: ECON, EREP, RMS 
151: DOUBLE PRECISION D2,R2AX,R2AY,R2AZ,R2BX,R2BY,R2BZ 
152: DOUBLE PRECISION G2(3) 
153: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12 
154: DOUBLE PRECISION XYZ(3*NATOMS), GGG(3*NATOMS), CCLOCAL 
155:  
156: GGG(1:3*NATOMS)=0.0D0 
157: ECON=0.0D0; EREP=0.0D0 
158: ! 
159: !  Constraint energy and forces. 
160: ! 
161: ! WRITE(MYUNIT,'(A,I8)') 'qcipot> NCONSTRAINT=',NCONSTRAINT 
162: DO J2=1,NCONSTRAINT 
163: !  WRITE(MYUNIT,'(A,I8,L5)') 'qcipot> J2,CONACTIVE=',J2,CONACTIVE(J2) 
164:    IF (.NOT.CONACTIVE(J2)) CYCLE 
165: !  CCLOCAL=CONCUTLOCAL(J2) 
166: !  IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 
167: !  IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 
168:    NI1=3*(CONI(J2)-1) 
169:    NJ1=3*(CONJ(J2)-1) 
170:    R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3) 
171:    R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3) 
172:    D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 
173: !  WRITE(MYUNIT,'(A,I8,3G20.10)') 'J2,D2,CONDISTREFLOCAL,CCLOCAL=',J2,D2,CONDISTREFLOCAL(J2),CCLOCAL 
174:    DUMMY=D2-CONDISTREFLOCAL(J2)   
175:    G2(1)=(R2AX-R2BX)/D2 
176:    G2(2)=(R2AY-R2BY)/D2 
177:    G2(3)=(R2AZ-R2BZ)/D2 
178:    REPGRAD(1:3)=2*INTCONSTRAINTDEL*DUMMY**2*DUMMY*G2(1:3) 
179:    DUMMY=INTCONSTRAINTDEL*DUMMY**4/2.0D0 
180:    ECON=ECON        +DUMMY 
181: !  WRITE(MYUNIT,'(A,I8,2G20.10)') 'J2,DUMMY,ECON=',J2,DUMMY,ECON 
182:    GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 
183:    GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 
184: !  WRITE(MYUNIT,'(A,6G20.10)') 'con grads=',GGG(NI1+1:NI1+3),GGG(NJ1+1:NJ1+3) 
185: ENDDO 
186:  
187: DO J2=1,NNREPULSIVE 
188:    INTCONST=NREPCUT(J2)**3 
189:    NI2=3*(NREPI(J2)-1) 
190:    NJ2=3*(NREPJ(J2)-1) 
191:    R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3) 
192:    R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3) 
193:    D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 
194:    IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1 
195: !     D12=D2**12 
196:       D12=D2**2 
197: !     DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST) 
198:       DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST) 
199: !     WRITE(MYUNIT,'(A,I8,2G20.10)') 'J2,DUMMY,EREP=',J2,DUMMY,EREP 
200:       EREP=EREP+DUMMY 
201: !     DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST) 
202:       DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST) 
203:       G2(1)=(R2AX-R2BX)/D2 
204:       G2(2)=(R2AY-R2BY)/D2 
205:       G2(3)=(R2AZ-R2BZ)/D2 
206:       REPGRAD(1:3)=DUMMY*G2(1:3) 
207:       GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3) 
208:       GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3) 
209: !     WRITE(MYUNIT,'(A,6G20.10)') 'rep grads=',GGG(NI2+1:NI2+3),GGG(NJ2+1:NJ2+3) 
210:    ENDIF 
211: ENDDO 
212: ! 
213: ! Set gradients on frozen atoms to zero. 
214: ! 
215: IF (FREEZE) THEN 
216:    DO J2=1,NATOMS 
217:       IF (FROZEN(J2)) THEN 
218:          GGG(3*(J2-1)+1)=0.0D0 
219:          GGG(3*(J2-1)+2)=0.0D0 
220:          GGG(3*(J2-1)+3)=0.0D0 
221:       ENDIF 
222:    ENDDO 
223: ENDIF 
224:  
225: END SUBROUTINE QCIPOT2 


r29153/quench.F 2015-11-17 23:33:37.429006709 +0000 r29152/quench.F 2015-11-17 23:33:40.009041309 +0000
312:             WRITE(MYUNIT,'(A,2L5)') 'QCIPOTT,INTCONSTRAINTT=',QCIPOTT,INTCONSTRAINTT312:             WRITE(MYUNIT,'(A,2L5)') 'QCIPOTT,INTCONSTRAINTT=',QCIPOTT,INTCONSTRAINTT
313:             CALL CPU_TIME(T_LBFGS_START)313:             CALL CPU_TIME(T_LBFGS_START)
314:             CALL INTLBFGS(P,FINISH)314:             CALL INTLBFGS(P,FINISH)
315:             CALL CPU_TIME(T_LBFGS_END)315:             CALL CPU_TIME(T_LBFGS_END)
316:             STOP316:             STOP
317:          ELSE317:          ELSE
318: ! Added timing for LBFGS call318: ! Added timing for LBFGS call
319:             CALL CPU_TIME(T_LBFGS_START)319:             CALL CPU_TIME(T_LBFGS_START)
320:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)320:             CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP)
321:             CALL CPU_TIME(T_LBFGS_END)321:             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_START322:             IF (FEBHT) THEN
 323:                WRITE(MYUNIT, '(A, F10.3)') 'Time to minimise (s):', T_LBFGS_END - T_LBFGS_START
 324:             ENDIF
323:          ENDIF325:          ENDIF
324:          IF (EVAPREJECT) RETURN326:          IF (EVAPREJECT) RETURN
325:       ENDIF327:       ENDIF
326:       !328:       !
327:       IF (FEBHT .AND. CFLAG) THEN329:       IF (FEBHT .AND. CFLAG) THEN
328:       ! Calculate the free energy330:       ! Calculate the free energy
329:          NUM_ZERO_EVS=6331:          NUM_ZERO_EVS=6
330:          IF (NATOMS.EQ.2) NUM_ZERO_EVS=5332:          IF (NATOMS.EQ.2) NUM_ZERO_EVS=5
331:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)333:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
332:          ALLOCATE(HESS(3*NATOMS,3*NATOMS))334:          ALLOCATE(HESS(3*NATOMS,3*NATOMS))
336:          IF (SPARSET) THEN338:          IF (SPARSET) THEN
337:              IF (DEBUG) THEN339:              IF (DEBUG) THEN
338:                 WRITE(MYUNIT, *) 'Using sparse cholesky decomposition to calculate determinant.'340:                 WRITE(MYUNIT, *) 'Using sparse cholesky decomposition to calculate determinant.'
339:              END IF341:              END IF
340:              SHIFT_ARRAY(:) = 1.0D0342:              SHIFT_ARRAY(:) = 1.0D0
341:              CALL CPU_TIME(T_DIAG_START)343:              CALL CPU_TIME(T_DIAG_START)
342:              CALL SHIFT_HESS_ZEROS(P, SHIFT_ARRAY)344:              CALL SHIFT_HESS_ZEROS(P, SHIFT_ARRAY)
343:              CALL FILTER_ZEROS(HESS, ZERO_THRESH)345:              CALL FILTER_ZEROS(HESS, ZERO_THRESH)
344:              CALL GET_DETERMINANT(3*NATOMS, LOG_PROD_EV, NQ(NP))346:              CALL GET_DETERMINANT(3*NATOMS, LOG_PROD_EV, NQ(NP))
345:              CALL CPU_TIME(T_DIAG_END)347:              CALL CPU_TIME(T_DIAG_END)
346:              IF (DEBUG) WRITE(MYUNIT, '(A, F10.3)') 'Time to diagonalise with sparse routines (s):', T_DIAG_END - T_DIAG_START348:              WRITE(MYUNIT, '(A, F10.3)') 'Time to diagonalise with sparse routines (s):', T_DIAG_END - T_DIAG_START
347:              IF (LOG_PROD_EV < SMALL_DOUBLE) THEN349:              IF (LOG_PROD_EV < SMALL_DOUBLE) THEN
348:                 WRITE(MYUNIT, '(A)') 'Log product from SuiteSparse is -infinity. The hessian ' //350:                 WRITE(MYUNIT, '(A)') 'Log product from SuiteSparse is -infinity. The hessian ' //
349:      &                               'is probably not positive definite. If using a threshold ' //351:      &                               'is probably not positive definite. If using a threshold ' //
350:      &                               'on zero values, try reducing it.'352:      &                               'on zero values, try reducing it.'
351: !353: !
352: ! Changed to convergence failure so that the run can continue, as for ts location. DJW354: ! Changed to convergence failure so that the run can continue, as for ts location. DJW
353: !355: !
354:                 CFLAG=.FALSE.356:                 CFLAG=.FALSE.
355:                 RETURN357:                 RETURN
356:              END IF 358:              END IF 
361: !             CALL MASSWT()363: !             CALL MASSWT()
362: !             CALL SHIFT_HESS_ZEROS(P, (/ 1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0 /))364: !             CALL SHIFT_HESS_ZEROS(P, (/ 1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0, 1.0D0 /))
363: !             CALL DSYEV('N','L',3*NATOMS,HESS,3*NATOMS,EVALUES,WORK,LWORK,INFO)             365: !             CALL DSYEV('N','L',3*NATOMS,HESS,3*NATOMS,EVALUES,WORK,LWORK,INFO)             
364: !             PRINT *, "DSYEV:", SUM(LOG(EVALUES))366: !             PRINT *, "DSYEV:", SUM(LOG(EVALUES))
365: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!367: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
366:          ELSE368:          ELSE
367: #endif /* __SPARSE */369: #endif /* __SPARSE */
368:             CALL CPU_TIME(T_DIAG_START)370:             CALL CPU_TIME(T_DIAG_START)
369:             CALL DSYEV('N','L',3*NATOMS,HESS,3*NATOMS,EVALUES,WORK,LWORK,INFO)371:             CALL DSYEV('N','L',3*NATOMS,HESS,3*NATOMS,EVALUES,WORK,LWORK,INFO)
370:             CALL CPU_TIME(T_DIAG_END)372:             CALL CPU_TIME(T_DIAG_END)
371:             IF (DEBUG) WRITE(MYUNIT, '(A, F10.3)') 'Time to diagonalise with DSYEV (s):', T_DIAG_END - T_DIAG_START373:             WRITE(MYUNIT, '(A, F10.3)') 'Time to diagonalise with DSYEV (s):', T_DIAG_END - T_DIAG_START
372:             IF (DEBUG) THEN374:             IF (DEBUG) THEN
373:                WRITE(MYUNIT, *) 'Using DSYEV to calculate determinant.'375:                WRITE(MYUNIT, *) 'Using DSYEV to calculate determinant.'
374:             END IF376:             END IF
375:             IF (DEBUG) THEN377:             IF (DEBUG) THEN
376:                 WRITE(MYUNIT, '(A)') "Eigenvalues:"378:                 WRITE(MYUNIT, '(A)') "Eigenvalues:"
377:                 WRITE(MYUNIT, '(3F12.3)') EVALUES379:                 WRITE(MYUNIT, '(3F12.3)') EVALUES
378:             END IF380:             END IF
379: ! Check that we don't have any proper negative eigenvalues (the magnitude should drop).381: ! Check that we don't have any proper negative eigenvalues (the magnitude should drop).
380: ! We do it this way, since some of the zeros come out negative as well, so you can't 382: ! We do it this way, since some of the zeros come out negative as well, so you can't 
381: ! just check the sign.383: ! just check the sign.


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0