hdiff output

r33542/congrad.f90 2017-12-03 16:30:13.683364607 +0000 r33541/congrad.f90 2017-12-03 16:30:14.223371667 +0000
 16: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, & 16: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
 17:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, & 17:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
 18:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, JMAXCON, & 18:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, JMAXCON, &
 19:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, & 19:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, &
 20:   &  FREEZENODEST, INTSPRINGACTIVET, INTMINFAC, NCONOFF, CONOFFLIST, CONOFFTRIED, INTCONMAX, ECON, EREP, ESPRING, & 20:   &  FREEZENODEST, INTSPRINGACTIVET, INTMINFAC, NCONOFF, CONOFFLIST, CONOFFTRIED, INTCONMAX, ECON, EREP, ESPRING, &
 21:   &  CONVERGECONTEST, CONVERGEREPTEST, FCONTEST, FREPTEST, QCIINTREPMINSEP, QCIAVDEV 21:   &  CONVERGECONTEST, CONVERGEREPTEST, FCONTEST, FREPTEST, QCIINTREPMINSEP, QCIAVDEV
 22: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG 22: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
 23: USE PORFUNCS 23: USE PORFUNCS
 24: IMPLICIT NONE 24: IMPLICIT NONE
 25:             25:            
 26: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,ISTAT,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF,NMAXINT,NMININT 26: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF
 27: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF, FMAX, FMIN, SEPARATION 27: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF, FMAX, FMIN, SEPARATION
 28: INTEGER OFFSET1, OFFSET2 28: INTEGER JJMAX(INTIMAGE+2)
  29: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
 29: 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
 30: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3) 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
 31: DOUBLE PRECISION DUMMY, REPGRAD(3), D12, DSQ2, DSQ1, DSQI 32: DOUBLE PRECISION DUMMY, REPGRAD(3), D12, DSQ2, DSQ1, DSQI
  33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
 32: LOGICAL NOINT 34: LOGICAL NOINT
 33: 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
 34: DOUBLE PRECISION GGGR((3*NATOMS)*(INTIMAGE+2)), EREP1,EREP2 
 35: LOGICAL IMGFREEZE(INTIMAGE) 36: LOGICAL IMGFREEZE(INTIMAGE)
 36: DOUBLE PRECISION DPLUS, SPGRAD(3), DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq,CUTMAX,DISTMAX 37: DOUBLE PRECISION DPLUS, SPGRAD(3), DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq,CUTMAX,DISTMAX
 37: DOUBLE PRECISION CONDMAX, CONREFMAX, CONCUTMAX, DUMMY2, CCLOCAL2, DVEC(INTIMAGE+1), DEVIATION(INTIMAGE+1) 38: DOUBLE PRECISION CONDMAX, CONREFMAX, CONCUTMAX, DUMMY2, CCLOCAL2, DVEC(INTIMAGE+1), DEVIATION(INTIMAGE+1)
 38: DOUBLE PRECISION GLOCAL1(3*NATOMS), GLOCAL2(3*NATOMS), XYZ1(3*NATOMS), XYZ2(3*NATOMS) 
 39:  39: 
 40: EEE(1:INTIMAGE+2)=0.0D0 40: EEE(1:INTIMAGE+2)=0.0D0
  41: CONE(1:INTIMAGE+2)=0.0D0
  42: REPE(1:INTIMAGE+2)=0.0D0
  43: REPEINT(1:INTIMAGE+2)=0.0D0
  44: NREPINT(1:INTIMAGE+2)=0
 41: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 45: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0
 42: ECON=0.0D0; EREP=0.0D0 46: ECON=0.0D0; EREP=0.0D0
 43: NMAXINT=0; NMININT=0 ! not actually used 47: NINTMIN=0
  48: NINTMIN2=0
 44: MYUNIT=6 49: MYUNIT=6
 45:  50: 
 46: EMAX=-1.0D200 51: EMAX=-1.0D200
 47: FMAX=-1.0D200 52: FMAX=-1.0D200
 48: FMIN=1.0D200 53: FMIN=1.0D200
 49: ! EEMAX(1:INTIMAGE+2)=-1.0D200 54: EEMAX(1:INTIMAGE+2)=-1.0D200
 50: ! JJMAX(1:INTIMAGE+2)=-1 55: JJMAX(1:INTIMAGE+2)=-1
 51: JMAX=-1 56: JMAX=-1
 52: IMAX=-1 57: IMAX=-1
 53: EMAXNOFF=-1.0D200 58: EMAXNOFF=-1.0D200
 54: JMAXNOFF=-1 59: JMAXNOFF=-1
 55: IMAXNOFF=-1 60: IMAXNOFF=-1
 56: IF (INTCONSTRAINTDEL.EQ.0.0D0) GOTO 531 61: IF (INTCONSTRAINTDEL.EQ.0.0D0) GOTO 531
 57: ! 62: !
 58: !  Constraint energy and forces. 63: !  Constraint energy and forces.
 59: ! 64: !
 60: DO J2=1,NCONSTRAINT 65: DO J2=1,NCONSTRAINT
121:             CONREFMAX=CONDISTREFLOCAL(J2)126:             CONREFMAX=CONDISTREFLOCAL(J2)
122:             CONCUTMAX=CCLOCAL127:             CONCUTMAX=CCLOCAL
123:          ENDIF128:          ENDIF
124:          IF (DUMMY.GT.EMAXNOFF) THEN129:          IF (DUMMY.GT.EMAXNOFF) THEN
125:             IF (.NOT.CONOFFTRIED(J2)) THEN130:             IF (.NOT.CONOFFTRIED(J2)) THEN
126:                IMAXNOFF=J1131:                IMAXNOFF=J1
127:                JMAXNOFF=J2132:                JMAXNOFF=J2
128:                EMAXNOFF=DUMMY133:                EMAXNOFF=DUMMY
129:             ENDIF134:             ENDIF
130:          ENDIF135:          ENDIF
 136:          IF (DUMMY.GT.EEMAX(J1)) THEN
 137:             JJMAX(J1)=J2
 138:             EEMAX(J1)=DUMMY
 139:          ENDIF
 140:          CONE(J1)=CONE(J1)+DUMMY
131:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)141:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
132:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)142:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
133:          DUMMY2=MINVAL(REPGRAD)143:          DUMMY2=MINVAL(REPGRAD)
134:          IF (DUMMY2.LT.FMIN) FMIN=DUMMY144:          IF (DUMMY2.LT.FMIN) FMIN=DUMMY
135:          DUMMY2=MAXVAL(REPGRAD)145:          DUMMY2=MAXVAL(REPGRAD)
136:          IF (DUMMY2.GT.FMAX) FMAX=DUMMY146:          IF (DUMMY2.GT.FMAX) FMAX=DUMMY
137:       ENDIF147:       ENDIF
 148: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1)
138:    ENDDO149:    ENDDO
139: ENDDO150: ENDDO
140: IF (-FMIN.GT.FMAX) FMAX=-FMIN151: IF (-FMIN.GT.FMAX) FMAX=-FMIN
141: FCONTEST=FMAX152: FCONTEST=FMAX
142: IF (JMAX.GT.0) THEN153: IF (JMAX.GT.0) THEN
143:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,3G15.5,A,G15.5)') ' congrad> Highest constraint for any image in ',IMAX, &154:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,3G15.5,A,G15.5)') ' congrad> Highest constraint for any image in ',IMAX, &
144:  & ' con ',JMAX, ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX, &155:  & ' con ',JMAX, ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX, &
145:  & ' max grad=',FMAX156:  & ' max grad=',FMAX
146: 157: 
147: ENDIF158: ENDIF
158: ! ENDIF169: ! ENDIF
159: ! JMAXCON=JMAXNOFF170: ! JMAXCON=JMAXNOFF
160: 531 CONTINUE171: 531 CONTINUE
161: 172: 
162: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes173: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
163: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes174: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
164: 175: 
165: ! INTCONST=INTCONSTRAINREPCUT**13176: ! INTCONST=INTCONSTRAINREPCUT**13
166: 177: 
167: EMAX=-1.0D200178: EMAX=-1.0D200
168: FMAX=-1.0D200179: EEMAX(1:INTIMAGE+2)=-1.0D200
169: FMIN=1.0D200180: JJMAX(1:INTIMAGE+2)=-1
170: JMAX=-1181: JMAX=-1
171: IMAX=-1182: IMAX=-1
172: 183: 
173: GGGR(1:3*NATOMS*(INTIMAGE+2))=0.0D0 
174: IF (INTCONSTRAINTREP.EQ.0.0D0) GOTO 654184: IF (INTCONSTRAINTREP.EQ.0.0D0) GOTO 654
175: DO J1=2,INTIMAGE+1185: DO J1=2,INTIMAGE+2
176: ! DO J1=2,INTIMAGE+2 ! we don't do anything for INTIMAGE+2 any more 
177: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images186: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images
178:    IF (FREEZENODEST) THEN187:    IF (FREEZENODEST) THEN
179:       IF (J1.EQ.2) THEN188:       IF (J1.EQ.2) THEN
180:          IF (IMGFREEZE(1)) CYCLE189:          IF (IMGFREEZE(1)) CYCLE
181: !     ELSE IF (J1.EQ.INTIMAGE+2) THEN190:       ELSE IF (J1.EQ.INTIMAGE+2) THEN
182: !        IF (IMGFREEZE(INTIMAGE)) CYCLE191:          IF (IMGFREEZE(INTIMAGE)) CYCLE
183:       ELSE192:       ELSE
184:          IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE193:          IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
185:       ENDIF194:       ENDIF
186:    ENDIF195:    ENDIF
187:    OFFSET2=(3*NATOMS)*(J1-1) 
188:    OFFSET1=(3*NATOMS)*(J1-2) 
189:    XYZ1(1:3*NATOMS)=XYZ(OFFSET1+1:OFFSET1+3*NATOMS) 
190:    XYZ2(1:3*NATOMS)=XYZ(OFFSET2+1:OFFSET2+3*NATOMS) 
191:    GLOCAL1(1:3*NATOMS)=0.0D0 
192:    GLOCAL2(1:3*NATOMS)=0.0D0 
193:    EREP1=0.0D0 
194:    EREP2=0.0D0 
195:    DO J2=1,NNREPULSIVE196:    DO J2=1,NNREPULSIVE
196: ! !  INTCONST=NREPCUT(J2)**13197: ! !  INTCONST=NREPCUT(J2)**13
197: !    INTCONST=NREPCUT(J2)**3198: !    INTCONST=NREPCUT(J2)**3
198: !      IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN199: !      IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
199: !!        WRITE(*, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)200: !!        WRITE(*, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
200: !         STOP201: !         STOP
201: !      ENDIF202: !      ENDIF
202: 203: 
203: !     NI1=OFFSET1+3*(NREPI(J2)-1)204:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
204: !     NI2=OFFSET2+3*(NREPI(J2)-1)205:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
205: !     NJ1=OFFSET1+3*(NREPJ(J2)-1)206:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
206: !     NJ2=OFFSET2+3*(NREPJ(J2)-1)207:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
207:  
208:       NI1=3*(NREPI(J2)-1) 
209:       NI2=3*(NREPI(J2)-1) 
210:       NJ1=3*(NREPJ(J2)-1) 
211:       NJ2=3*(NREPJ(J2)-1) 
212:  
213: !     G1(1:3)=XYZ(NI1+1:NI1+3)-XYZ(NJ1+1:NJ1+3) 
214: !     G2(1:3)=XYZ(NI2+1:NI2+3)-XYZ(NJ2+1:NJ2+3) 
215: 208: 
216:       G1(1:3)=XYZ1(NI1+1:NI1+3)-XYZ1(NJ1+1:NJ1+3)209:       G1(1:3)=XYZ(NI1+1:NI1+3)-XYZ(NJ1+1:NJ1+3)
217:       G2(1:3)=XYZ2(NI2+1:NI2+3)-XYZ2(NJ2+1:NJ2+3)210:       G2(1:3)=XYZ(NI2+1:NI2+3)-XYZ(NJ2+1:NJ2+3)
218: !211: !
219: ! Squared distance between atoms A and B for theta=0 - distance in image 2212: ! Squared distance between atoms A and B for theta=0 - distance in image 2
220: !213: !
221:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2214:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2
222: !215: !
223: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1216: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1
224: !217: !
225:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2218:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2
226:       DCUT=NREPCUT(J2)**2219:       DCUT=NREPCUT(J2)**2
227:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff220:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff
228:       IF (ABS(NREPI(J2)-NREPJ(J2)).GT.QCIINTREPMINSEP) THEN ! don't check for internal minimum in distance - atoms too close for chain crossing.221:       IF (ABS(NREPI(J2)-NREPJ(J2)).GT.QCIINTREPMINSEP) THEN ! don't check for internal minimum in distance - atoms too close for chain crossing.
229:          r1apr2bmr2amr1bsq=0.0D0222:          r1apr2bmr2amr1bsq=0.0D0
230:       ELSE223:       ELSE
231: !        WRITE(*,'(A,I6,A,I6,A,2G20.10,A,G20.10)') 'distances in ',J1-1,' and ',J1,' are ',SQRT(DSQ1),SQRT(DSQ2),' cutoff=',NREPCUT(J2)224: !        WRITE(*,'(A,I6,A,I6,A,2G20.10,A,G20.10)') 'distances in ',J1-1,' and ',J1,' are ',SQRT(DSQ1),SQRT(DSQ2),' cutoff=',NREPCUT(J2)
232: !        WRITE(*,'(A,I6,A,2I6,A,6G15.7)') 'image ',J1-1,' atoms ',NREPI(J2),NREPJ(J2),' coords ',XYZ(NI1+1:NI1+3),XYZ(NJ1+1:NJ1+3)225: !        WRITE(*,'(A,I6,A,2I6,A,6G15.7)') 'image ',J1-1,' atoms ',NREPI(J2),NREPJ(J2),' coords ',XYZ(NI1+1:NI1+3),XYZ(NJ1+1:NJ1+3)
233: !        WRITE(*,'(A,I6,A,2I6,A,6G15.7)') 'image ',J1  ,' atoms ',NREPI(J2),NREPJ(J2),' coords ',XYZ(NI2+1:NI2+3),XYZ(NJ2+1:NJ2+3)226: !        WRITE(*,'(A,I6,A,2I6,A,6G15.7)') 'image ',J1  ,' atoms ',NREPI(J2),NREPJ(J2),' coords ',XYZ(NI2+1:NI2+3),XYZ(NJ2+1:NJ2+3)
234:          r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3)227:          r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3)
235:          r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b228: !        WRITE(*,'(A,G20.10)') 'dot product (r1I-r1J).(r2I-r2J)=',r1amr1bdr2amr2b
236:       ENDIF 
237: !229: !
238: ! Is the denominator of the d^2 internal minimum term close to zero?230: ! Is the denominator of the d^2 internal minimum term close to zero?
239: !231: !
 232:          r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b
 233: !        WRITE(*,'(A,G20.10)') 'internal minimum test d(J1-1)^2+d(J1)^2-2(r1I-r1J).(r2I-r2J)=',r1apr2bmr2amr1bsq
 234:       ENDIF
 235: 
240:       IF (r1apr2bmr2amr1bsq.LT.1.0D-50) THEN236:       IF (r1apr2bmr2amr1bsq.LT.1.0D-50) THEN
241:          NOINT=.TRUE.237:          NOINT=.TRUE.
242:          D1=SQRT(DSQ1)238:          D1=SQRT(DSQ1)
243:          D2=SQRT(DSQ2)239:          D2=SQRT(DSQ2)
244:          G2(1:3)=G2(1:3)/D2240:          G2(1:3)=G2(1:3)/D2
245:          G1(1:3)=G1(1:3)/D1241:          G1(1:3)=G1(1:3)/D1
246:       ELSE242:       ELSE
247:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)243:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)
248:       ENDIF244:       ENDIF
 245: !     WRITE(*,'(A,8G15.7)') 'D1,D2,G2,G1=',D1,D2,G2(1:3),G1(1:3)
249: !246: !
250: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.247: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
251: !248: !
252: ! Multiply energy and gradient by NREPCUT**2/INTCONSTRAINTREP to put values on a common scale.249: ! Multiply energy and gradient by NREPCUT**2/INTCONSTRAINTREP to put values on a common scale.
253: !250: !
254: !     IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1251:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
255:       IF (D2.LT.NREPCUT(J2)) THEN ! terms for image J1 - non-zero derivatives only for J1252: ! !        D12=DSQ2**6
256: !        D12=DSQ2**6 
257:          D12=DSQ2253:          D12=DSQ2
258: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)254: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
259: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)255: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
260: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)256: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
261: !        DUMMY=(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2257: !        DUMMY=(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2
262:          DUMMY=NREPCUT(J2)**2/D12+2.0D0*D2/NREPCUT(J2)-3.0D0  258:          DUMMY=NREPCUT(J2)**2/D12+2.0D0*D2/NREPCUT(J2)-3.0D0  
263: !        EEE(J1)=EEE(J1)+DUMMY259:          EEE(J1)=EEE(J1)+DUMMY
264:          EREP1=EREP1+DUMMY260: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
 261: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
 262: ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
 263: !        ENDIF
265:          IF (DUMMY.GT.EMAX) THEN264:          IF (DUMMY.GT.EMAX) THEN
266:             IMAX=J1265:             IMAX=J1
267:             JMAX=J2266:             JMAX=J2
268:             CUTMAX=NREPCUT(J2)267:             CUTMAX=NREPCUT(J2)
269:             DISTMAX=D2268:             DISTMAX=D2
270:             EMAX=DUMMY269:             EMAX=DUMMY
271:          ENDIF270:          ENDIF
 271:          IF (DUMMY.GT.EEMAX(J1)) THEN
 272:             JJMAX(J1)=J2
 273:             EEMAX(J1)=DUMMY
 274:          ENDIF
 275:          REPE(J1)=REPE(J1)+DUMMY
272:          EREP=EREP+DUMMY276:          EREP=EREP+DUMMY
273: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)277: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
274: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)278: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
 279: !
275: !        DUMMY=-2.0D0*(1.0D0/(D2*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2280: !        DUMMY=-2.0D0*(1.0D0/(D2*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2
276:  
277:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(D2*D12)-1.0D0/NREPCUT(J2))281:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(D2*D12)-1.0D0/NREPCUT(J2))
278:          REPGRAD(1:3)=DUMMY*G2(1:3)282:          REPGRAD(1:3)=DUMMY*G2(1:3)
279: !        GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)283: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
280: !        GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)284: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
281:          GLOCAL2(NI2+1:NI2+3)=GLOCAL2(NI2+1:NI2+3)+REPGRAD(1:3)285:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
282:          GLOCAL2(NJ2+1:NJ2+3)=GLOCAL2(NJ2+1:NJ2+3)-REPGRAD(1:3)286:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
 287:          DUMMY2=MINVAL(REPGRAD)
 288:          IF (DUMMY2.LT.FMIN) THEN
 289:             FMIN=DUMMY2
 290:          ENDIF
 291:          DUMMY2=MAXVAL(REPGRAD)
 292:          IF (DUMMY2.GT.FMAX) THEN
 293:             FMAX=DUMMY2
 294:          ENDIF
283:       ENDIF295:       ENDIF
284: !296: !
285: ! For internal minima we are counting edges. 297: ! For internal minima we are counting edges. 
286: ! Edge J1 is between images J1-1 and J1, starting from J1=2.298: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
287: ! Energy contributions are shared evenly, except for299: ! Energy contributions are shared evenly, except for
288: ! edge 1, which was assigned to image 2, and edge INTIMAGE+1, which300: ! edge 1, which was assigned to image 2, and edge INTIMAGE+1, which
289: ! was assigned to image INTIMAGE+1. Gradients are set to zero for the end images.301: ! was assigned to image INTIMAGE+1. Gradients are set to zero for the end images.
290: ! 20/11/17 - changed to turn off internal minimum checks involving the fixed endpoints. DJW302: ! 20/11/17 - changed to turn off internal minimum checks involving the fixed endpoints. DJW
291: !303: !
292:       DUMMY=0.0D0304:       DUMMY=0.0D0
293:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2)).AND.(J1.NE.2)) THEN305: !     WRITE(*,'(A,2G20.10)') 'DINT,NREPCUT=',DINT,NREPCUT(J2)
 306:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2)).AND.(J1.NE.2).AND.(J1.NE.INTIMAGE+2)) THEN
 307:          NINTMIN2=NINTMIN2+1
294: !        D12=DSQI**6308: !        D12=DSQI**6
295:          D12=DSQI309:          D12=DSQI
296: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)310: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
297: !        DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)311: !        DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
298: !        DUMMY=INTMINFAC*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2312: !        DUMMY=INTMINFAC*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2
299:          DUMMY=INTMINFAC*(NREPCUT(J2)**2/D12+2.0D0*DINT/NREPCUT(J2)-3.0D0)313:          DUMMY=INTMINFAC*(NREPCUT(J2)**2/D12+2.0D0*DINT/NREPCUT(J2)-3.0D0)
300: !        IF (J1.EQ.2) THEN314:          IF (J1.EQ.2) THEN
301: !           EEE(J1)=EEE(J1)+DUMMY315:             EEE(J1)=EEE(J1)+DUMMY
302: !        ELSE ! IF (J1.LT.INTIMAGE+2) THEN316:             REPEINT(J1)=REPEINT(J1)+DUMMY
303: !           EEE(J1)=EEE(J1)+DUMMY/2.0D0317:             NREPINT(J1)=NREPINT(J1)+1
304: !           EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0318:          ELSE IF (J1.LT.INTIMAGE+2) THEN
305: !        ENDIF319:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
306:          EREP2=EREP2+DUMMY320:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
 321:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0
 322:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
 323:             NREPINT(J1)=NREPINT(J1)+1
 324:             NREPINT(J1-1)=NREPINT(J1-1)+1
 325:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
 326:             EEE(J1-1)=EEE(J1-1)+DUMMY
 327:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
 328:             NREPINT(J1-1)=NREPINT(J1-1)+1
 329:          ENDIF
307:          EREP=EREP+DUMMY330:          EREP=EREP+DUMMY
308:          IF (DUMMY.GT.EMAX) THEN331:          IF (DUMMY.GT.EMAX) THEN
309:             IMAX=J1332:             IMAX=J1
310:             JMAX=J2333:             JMAX=J2
311:             EMAX=DUMMY334:             EMAX=DUMMY
312:             DISTMAX=DINT335:             DISTMAX=DINT
313:             CUTMAX=NREPCUT(J2)336:             CUTMAX=NREPCUT(J2)
314:          ENDIF337:          ENDIF
 338:          IF (DUMMY.GT.EEMAX(J1)) THEN
 339:             JJMAX(J1)=J2
 340:             EEMAX(J1)=DUMMY
 341:          ENDIF
315: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)342: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
316: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)343: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
317: !        DUMMY=-2.0D0*(1.0D0/(DINT*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2344: !        DUMMY=-2.0D0*(1.0D0/(DINT*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2
318:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(DINT*D12)-1.0D0/NREPCUT(J2))345:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(DINT*D12)-1.0D0/NREPCUT(J2))
319:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)346:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)
320: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &347: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
321: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)348: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
322: !349: !
323: ! Gradient contributions for image J1-1350: ! Gradient contributions for image J1-1
324: !351: !
325: 352: 
326: !        GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)353:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
327: !        GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)354:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
328:          GLOCAL1(NI1+1:NI1+3)=GLOCAL1(NI1+1:NI1+3)+REPGRAD(1:3) 
329:          GLOCAL1(NJ1+1:NJ1+3)=GLOCAL1(NJ1+1:NJ1+3)-REPGRAD(1:3) 
330:          DUMMY2=MINVAL(REPGRAD)355:          DUMMY2=MINVAL(REPGRAD)
 356:          IF (DUMMY2.LT.FMIN) THEN
 357:             FMIN=DUMMY2
 358: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX
 359:          ENDIF
 360:          DUMMY2=MAXVAL(REPGRAD)
 361:          IF (DUMMY2.GT.FMAX) THEN
 362:             FMAX=DUMMY2
 363: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX
 364:          ENDIF
 365: !        WRITE(*,'(A,2I6,5G20.10)') 'image J1-1 J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX
 366: 
 367: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
 368: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
331: !369: !
332: ! Gradient contributions for image J1370: ! Gradient contributions for image J1
333: !371: !
334:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)372:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)
335: !        GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)373:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
336: !        GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)374:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
337:          GLOCAL2(NI2+1:NI2+3)=GLOCAL2(NI2+1:NI2+3)+REPGRAD(1:3)375:          DUMMY2=MINVAL(REPGRAD)
338:          GLOCAL2(NJ2+1:NJ2+3)=GLOCAL2(NJ2+1:NJ2+3)-REPGRAD(1:3)376:          IF (DUMMY2.LT.FMIN) THEN
 377:             FMIN=DUMMY2
 378: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX
 379:          ENDIF
 380:          DUMMY2=MAXVAL(REPGRAD)
 381:          IF (DUMMY2.GT.FMAX) THEN
 382:             FMAX=DUMMY2
 383: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX
 384:          ENDIF
 385: !        WRITE(*,'(A,2I6,5G20.10)') 'image J1   J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX
339:       ENDIF386:       ENDIF
340:    ENDDO387:    ENDDO
341:    GGGR(OFFSET1+1:OFFSET1+3*NATOMS)=GGGR(OFFSET1+1:OFFSET1+3*NATOMS)+GLOCAL1(1:3*NATOMS) 
342:    GGGR(OFFSET2+1:OFFSET2+3*NATOMS)=GGGR(OFFSET2+1:OFFSET2+3*NATOMS)+GLOCAL2(1:3*NATOMS) 
343:    EEE(J1)=EEE(J1)+EREP1 
344:    IF (J1.EQ.2) THEN 
345:       EEE(J1)=EEE(J1)+EREP2 
346:    ELSE ! IF (J1.LT.INTIMAGE+2) THEN 
347:       EEE(J1)=EEE(J1)+EREP2/2.0D0 
348:       EEE(J1-1)=EEE(J1-1)+EREP2/2.0D0 
349:    ENDIF 
350: ENDDO388: ENDDO
351: FMIN=MINVAL(GGGR(3*NATOMS+1:3*NATOMS*(INTIMAGE+1))) 
352: FMAX=MAXVAL(GGGR(3*NATOMS+1:3*NATOMS*(INTIMAGE+1))) 
353: IF (-FMIN.GT.FMAX) FMAX=-FMIN389: IF (-FMIN.GT.FMAX) FMAX=-FMIN
354: FREPTEST=FMAX390: FREPTEST=FMAX
355:  
356: GGG(1:3*NATOMS*(INTIMAGE+2))=GGG(1:3*NATOMS*(INTIMAGE+2))+GGGR(1:3*NATOMS*(INTIMAGE+2)) 
357:  
358: IF (JMAX.GT.0) THEN391: IF (JMAX.GT.0) THEN
359:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,2G15.5,A,G15.5)') ' congrad> Highest repulsion  for any image in ',IMAX, &  392:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,2G15.5,A,G15.5)') ' congrad> Highest repulsion  for any image in ',IMAX, &  
360:  &  ' ind ',JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX),' value=',EMAX,' d,cutoff=',DISTMAX,CUTMAX, &393:  &  ' ind ',JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX),' value=',EMAX,' d,cutoff=',DISTMAX,CUTMAX, &
361:  &  ' max grad=',FMAX394:  &  ' max grad=',FMAX
362: ENDIF395: ENDIF
363: CONVERGEREPTEST=EMAX396: CONVERGEREPTEST=EMAX
364: 654 CONTINUE397: 654 CONTINUE
365: 398: 
366: !399: !
367: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise400: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
451:       ENDDO484:       ENDDO
452:    ENDDO485:    ENDDO
453: ENDIF486: ENDIF
454: !487: !
455: ! Set gradients to zero for start and finish images.488: ! Set gradients to zero for start and finish images.
456: !489: !
457: GGG(1:(3*NATOMS))=0.0D0490: GGG(1:(3*NATOMS))=0.0D0
458: GGG((INTIMAGE+1)*(3*NATOMS)+1:(INTIMAGE+2)*(3*NATOMS))=0.0D0491: GGG((INTIMAGE+1)*(3*NATOMS)+1:(INTIMAGE+2)*(3*NATOMS))=0.0D0
459: RMS=0.0D0492: RMS=0.0D0
460: DO J1=2,INTIMAGE+1493: DO J1=2,INTIMAGE+1
 494:    RMSIM(J1)=0.0D0
461:    DO J2=1,(3*NATOMS)495:    DO J2=1,(3*NATOMS)
462:       RMS=RMS+GGG((3*NATOMS)*(J1-1)+J2)**2496:       RMS=RMS+GGG((3*NATOMS)*(J1-1)+J2)**2
 497:       RMSIM(J1)=RMSIM(J1)+GGG((3*NATOMS)*(J1-1)+J2)**2
463:    ENDDO498:    ENDDO
 499:    RMSIM(J1)=SQRT(RMSIM(J1)/(3*NATOMS))
464: ENDDO500: ENDDO
465: IF (INTIMAGE.NE.0) THEN501: IF (INTIMAGE.NE.0) THEN
466:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))502:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))
467: ENDIF503: ENDIF
468: !504: !
469: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,505: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,
470: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)506: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)
471: !507: !
472: ETOTAL=0.0D0508: ETOTAL=0.0D0
 509: MAXINT=-1.0D100
 510: MININT=1.0D100
473: DO J1=2,INTIMAGE+1511: DO J1=2,INTIMAGE+1
474:    ETOTAL=ETOTAL+EEE(J1)512:    ETOTAL=ETOTAL+EEE(J1)
 513: !  WRITE(*, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)
 514:    IF (REPEINT(J1).LT.MININT) THEN
 515:       MININT=REPEINT(J1)
 516:       NMININT=J1
 517:    ENDIF
 518:    IF (REPE(J1).GT.MAXINT) THEN
 519:       MAXINT=REPE(J1)
 520:       NMAXINT=J1
 521:    ENDIF
475: ENDDO522: ENDDO
 523: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT
 524: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> smallest internal energy=',MININT,' for image ',NMININT
 525: ! IF (DEBUG) WRITE(*, '(A,2I6)') ' congrad> number of internal minima=',NINTMIN,NINTMIN2
 526: 
 527: ! STOP !!! DEBUG DJW
476: 528: 
477: END SUBROUTINE CONGRAD529: END SUBROUTINE CONGRAD
478: 530: 
479: SUBROUTINE MINMAXD2(D2,D1,DINT,DSQ2,DSQ1,G1,G2,G1INT,G2INT,NOINT,DEBUG,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)531: SUBROUTINE MINMAXD2(D2,D1,DINT,DSQ2,DSQ1,G1,G2,G1INT,G2INT,NOINT,DEBUG,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)
480: USE KEY, ONLY : CHECKCONINT532: USE KEY, ONLY : CHECKCONINT
481: IMPLICIT NONE533: IMPLICIT NONE
482: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT534: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT
483: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)535: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)
484: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq536: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq
485: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY, DUMMY2537: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY, DUMMY2


r33542/intlbfgs.f90 2017-12-03 16:30:13.963368213 +0000 r33541/intlbfgs.f90 2017-12-03 16:30:14.487375154 +0000
661: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.661: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.
662: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!662: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!
663: !663: !
664: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.664: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.
665: !665: !
666: BESTWORST=1.0D100666: BESTWORST=1.0D100
667: 9876 CONTINUE667: 9876 CONTINUE
668: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)668: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
669: 669: 
670: 670: 
671: !                DIFF=1.0D-6671: !                DIFF=1.0D-7
672: !                PRINT*,'analytic and numerical gradients:'672: !                PRINT*,'analytic and numerical gradients:'
673: !                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)673: !                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
674: !!               DO J1=2,INTIMAGE+1674: !!               DO J1=2,INTIMAGE+1
675: !                DO J1=INTIMAGE-1,INTIMAGE+1675: !                DO J1=INTIMAGE-1,INTIMAGE+1
676: !                   DO J2=1,3*NATOMS676: !                   DO J2=1,3*NATOMS
677: !                      IF (.NOT.ATOMACTIVE((J2-1)/3+1)) CYCLE677: !                      IF (.NOT.ATOMACTIVE((J2-1)/3+1)) CYCLE
678: !                      J3=3*NATOMS*(J1-1)+J2678: !                      J3=3*NATOMS*(J1-1)+J2
679: !                      XYZ(J3)=XYZ(J3)+DIFF679: !                      XYZ(J3)=XYZ(J3)+DIFF
680: !                      CALL CONGRAD(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)680: !                      CALL CONGRAD(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)
681: !                      XYZ(J3)=XYZ(J3)-2.0D0*DIFF681: !                      XYZ(J3)=XYZ(J3)-2.0D0*DIFF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0