hdiff output

r33499/congrad.f90 2017-11-18 14:30:14.672915855 +0000 r33498/congrad.f90 2017-11-18 14:30:15.116921720 +0000
 49: MYUNIT=6 49: MYUNIT=6
 50:  50: 
 51: EMAX=-1.0D200 51: EMAX=-1.0D200
 52: EEMAX(1:INTIMAGE+2)=-1.0D200 52: EEMAX(1:INTIMAGE+2)=-1.0D200
 53: JJMAX(1:INTIMAGE+2)=-1 53: JJMAX(1:INTIMAGE+2)=-1
 54: JMAX=-1 54: JMAX=-1
 55: IMAX=-1 55: IMAX=-1
 56: EMAXNOFF=-1.0D200 56: EMAXNOFF=-1.0D200
 57: JMAXNOFF=-1 57: JMAXNOFF=-1
 58: IMAXNOFF=-1 58: IMAXNOFF=-1
 59: IF (INTCONSTRAINTDEL.EQ.0.0D0) GOTO 531 
 60: ! 59: !
 61: !  Constraint energy and forces. 60: !  Constraint energy and forces.
 62: ! 61: !
 63: DO J2=1,NCONSTRAINT 62: DO J2=1,NCONSTRAINT
 64:    IF (.NOT.CONACTIVE(J2)) CYCLE 63:    IF (.NOT.CONACTIVE(J2)) CYCLE
 65:       CCLOCAL=CONCUTLOCAL(J2) 64:       CCLOCAL=CONCUTLOCAL(J2)
 66:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 65:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
 67:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 66:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
 68: ! 67: !
 69: ! For J1 we consider the line segment between image J1-1 and J1. 68: ! For J1 we consider the line segment between image J1-1 and J1.
 93:       ENDIF 92:       ENDIF
 94:       NI1=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1) 93:       NI1=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1)
 95:       NJ1=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1) 94:       NJ1=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1)
 96:       R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3) 95:       R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3)
 97:       R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3) 96:       R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3)
 98:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 97:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)
 99:       DUMMY=D2-CONDISTREFLOCAL(J2)   98:       DUMMY=D2-CONDISTREFLOCAL(J2)  
100:       DUMMY2=DUMMY**2 99:       DUMMY2=DUMMY**2
101:       CCLOCAL2=CCLOCAL**2100:       CCLOCAL2=CCLOCAL**2
102: !101: !
103: ! Reduced form for penalty function and gradient - multiply through by 1/(CCLOCAL**2*INTCONSTRAINTDEL)102: ! Reduced form for penalty function and gradient - multiply through by CCLOCAL**2/INTCONSTRAINTDEL
104: ! Should save CCLOCAL**2103: ! Should save CCLOCAL**2
105: ! Could save DUMMY**2-CCLOCAL**2 outside IF block and base the test on difference of squares104: ! Could save DUMMY**2-CCLOCAL**2 outside IF block and base the test on difference of squares
106: !105: !
107:       IF (DUMMY2.GT.CCLOCAL2) THEN 106:       IF (DUMMY2.GT.CCLOCAL2) THEN 
108:          G2(1)=(R2AX-R2BX)/D2107:          G2(1)=(R2AX-R2BX)/D2
109:          G2(2)=(R2AY-R2BY)/D2108:          G2(2)=(R2AY-R2BY)/D2
110:          G2(3)=(R2AZ-R2BZ)/D2109:          G2(3)=(R2AZ-R2BZ)/D2
111: !110: !
112: !        REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)111: !        REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
113: !        DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)112: !        DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
114: !113: !
115:          REPGRAD(1:3)=2*(DUMMY2-CCLOCAL2)*DUMMY*G2(1:3)/CCLOCAL2**2114:          REPGRAD(1:3)=2*(DUMMY2-CCLOCAL2)*DUMMY*G2(1:3)
116:          DUMMY=(DUMMY2-CCLOCAL2)**2/(2.0D0*CCLOCAL2**2)115:          DUMMY=(DUMMY2-CCLOCAL2)**2/2.0D0
117:          EEE(J1)=EEE(J1)  +DUMMY116:          EEE(J1)=EEE(J1)  +DUMMY
118:          ECON=ECON        +DUMMY117:          ECON=ECON        +DUMMY
119:          IF (DUMMY.GT.EMAX) THEN118:          IF (DUMMY.GT.EMAX) THEN
120:             IMAX=J1119:             IMAX=J1
121:             JMAX=J2120:             JMAX=J2
122:             EMAX=DUMMY121:             EMAX=DUMMY
123:             CONDMAX=D2122:             CONDMAX=D2
124:             CONREFMAX=CONDISTREFLOCAL(J2)123:             CONREFMAX=CONDISTREFLOCAL(J2)
125:             CONCUTMAX=CCLOCAL124:             CONCUTMAX=CCLOCAL
126:          ENDIF125:          ENDIF
136:             EEMAX(J1)=DUMMY135:             EEMAX(J1)=DUMMY
137:          ENDIF136:          ENDIF
138:          CONE(J1)=CONE(J1)+DUMMY137:          CONE(J1)=CONE(J1)+DUMMY
139:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)138:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
140:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)139:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
141:       ENDIF140:       ENDIF
142: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1)141: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1)
143:    ENDDO142:    ENDDO
144: ENDDO143: ENDDO
145: IF (JMAX.GT.0) THEN144: IF (JMAX.GT.0) THEN
146:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,3G15.5)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &145:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G20.10,A,3G20.10)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &
147:  & ' constraint ',JMAX, &146:  & ' constraint ',JMAX, &
148:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX147:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX
149: 148: 
150: ENDIF149: ENDIF
151: CONVERGECONTEST=EMAX150: CONVERGECONTEST=EMAX/INTCONSTRAINTDEL
152: IF (JMAXNOFF.GT.0) THEN151: IF (JMAXNOFF.GT.0) THEN
153:    WRITE(*,'(A,I6,A,I6,A,2I8,A,G20.10,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &152:    WRITE(*,'(A,I6,A,I6,A,2I8,A,G20.10,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &
154:  & ' constraint ',JMAXNOFF, &153:  & ' constraint ',JMAXNOFF, &
155:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' off=',NCONOFF154:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' off=',NCONOFF
156: ELSEIF (JMAX.GT.0) THEN155: ELSEIF (JMAX.GT.0) THEN
157:    JMAXNOFF=JMAX156:    JMAXNOFF=JMAX
158:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'157:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'
159:    CONOFFTRIED(1:INTCONMAX)=.FALSE.158:    CONOFFTRIED(1:INTCONMAX)=.FALSE.
160:    NCONOFF=0159:    NCONOFF=0
161: ENDIF160: ENDIF
162: JMAXCON=JMAXNOFF161: JMAXCON=JMAXNOFF
163: 531 CONTINUE 
164: 162: 
165: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes163: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
166: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes164: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
167: 165: 
168: ! INTCONST=INTCONSTRAINREPCUT**13166: ! INTCONST=INTCONSTRAINREPCUT**13
169: 167: 
170: EMAX=-1.0D200168: EMAX=-1.0D200
171: EEMAX(1:INTIMAGE+2)=-1.0D200169: EEMAX(1:INTIMAGE+2)=-1.0D200
172: JJMAX(1:INTIMAGE+2)=-1170: JJMAX(1:INTIMAGE+2)=-1
173: JMAX=-1171: JMAX=-1
174: IMAX=-1172: IMAX=-1
175: 173: 
176: IF (INTCONSTRAINTREP.EQ.0.0D0) GOTO 654 
177: DO J2=1,NNREPULSIVE174: DO J2=1,NNREPULSIVE
178: ! !  INTCONST=NREPCUT(J2)**13175: ! !  INTCONST=NREPCUT(J2)**13
179: !    INTCONST=NREPCUT(J2)**3176: !  INTCONST=NREPCUT(J2)**3
180:    DO J1=2,INTIMAGE+2177:    DO J1=2,INTIMAGE+2
181: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images178: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images
182:       IF (FREEZENODEST) THEN179:       IF (FREEZENODEST) THEN
183:          IF (J1.EQ.2) THEN180:          IF (J1.EQ.2) THEN
184:             IF (IMGFREEZE(1)) CYCLE181:             IF (IMGFREEZE(1)) CYCLE
185:          ELSE IF (J1.EQ.INTIMAGE+2) THEN182:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
186:             IF (IMGFREEZE(INTIMAGE)) CYCLE183:             IF (IMGFREEZE(INTIMAGE)) CYCLE
187:          ELSE184:          ELSE
188:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE185:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
189:          ENDIF186:          ENDIF
221:          D1=SQRT(DSQ1)218:          D1=SQRT(DSQ1)
222:          D2=SQRT(DSQ2)219:          D2=SQRT(DSQ2)
223:          G2(1:3)=G2(1:3)/D2220:          G2(1:3)=G2(1:3)/D2
224:          G1(1:3)=G1(1:3)/D1221:          G1(1:3)=G1(1:3)/D1
225:       ELSE222:       ELSE
226:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)223:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)
227:       ENDIF224:       ENDIF
228: !225: !
229: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.226: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
230: !227: !
231: ! Multiply energy and gradient by NREPCUT**2/INTCONSTRAINTREP to put values on a common scale. 
232: ! 
233:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1228:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
234: ! !        D12=DSQ2**6229: ! !        D12=DSQ2**6
235:          D12=DSQ2230:          D12=DSQ2
236: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)231: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
237: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)232: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
238: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)233: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
239: !        DUMMY=(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2234: !        DUMMY=(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2
240:          DUMMY=NREPCUT(J2)**2/D12+2.0D0*D2/NREPCUT(J2)-3.0D0  235:          DUMMY=DCUT/D12+2.0D0*D2/NREPCUT(J2)-3.0D0
241:          EEE(J1)=EEE(J1)+DUMMY236:          EEE(J1)=EEE(J1)+DUMMY
242: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN237: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
243: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &238: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
244: ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY239: ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
245: !        ENDIF240: !        ENDIF
246:          IF (DUMMY.GT.EMAX) THEN241:          IF (DUMMY.GT.EMAX) THEN
247:             IMAX=J1242:             IMAX=J1
248:             JMAX=J2243:             JMAX=J2
249:             CUTMAX=NREPCUT(J2)244:             CUTMAX=NREPCUT(J2)
250:             DISTMAX=D2245:             DISTMAX=D2
253:          IF (DUMMY.GT.EEMAX(J1)) THEN248:          IF (DUMMY.GT.EEMAX(J1)) THEN
254:             JJMAX(J1)=J2249:             JJMAX(J1)=J2
255:             EEMAX(J1)=DUMMY250:             EEMAX(J1)=DUMMY
256:          ENDIF251:          ENDIF
257:          REPE(J1)=REPE(J1)+DUMMY252:          REPE(J1)=REPE(J1)+DUMMY
258:          EREP=EREP+DUMMY253:          EREP=EREP+DUMMY
259: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)254: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
260: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)255: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
261: !256: !
262: !        DUMMY=-2.0D0*(1.0D0/(D2*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2257: !        DUMMY=-2.0D0*(1.0D0/(D2*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2
263:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(D2*D12)-1.0D0/NREPCUT(J2))258:          DUMMY=-2.0D0*(DCUT/(D2*D12)-1.0D0/NREPCUT(J2))
264:          REPGRAD(1:3)=DUMMY*G2(1:3)259:          REPGRAD(1:3)=DUMMY*G2(1:3)
265: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &260: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
266: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)261: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
267:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)262:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
268:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)263:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
269:       ENDIF264:       ENDIF
270: !265: !
271: ! For internal minima we are counting edges. 266: ! For internal minima we are counting edges. 
272: ! Edge J1 is between images J1-1 and J1, starting from J1=2.267: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
273: ! Energy contributions are shared evenly, except for268: ! Energy contributions are shared evenly, except for
276: ! the end images.271: ! the end images.
277: !272: !
278:       DUMMY=0.0D0273:       DUMMY=0.0D0
279:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN274:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
280:          NINTMIN2=NINTMIN2+1275:          NINTMIN2=NINTMIN2+1
281: !        D12=DSQI**6276: !        D12=DSQI**6
282:          D12=DSQI277:          D12=DSQI
283: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)278: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
284: !        DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)279: !        DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
285: !        DUMMY=INTMINFAC*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2280: !        DUMMY=INTMINFAC*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2
286:          DUMMY=INTMINFAC*(NREPCUT(J2)**2/D12+2.0D0*DINT/NREPCUT(J2)-3.0D0)281:          DUMMY=INTMINFAC*(DCUT/D12+2.0D0*DINT/NREPCUT(J2)-3.0D0)
287:          IF (J1.EQ.2) THEN282:          IF (J1.EQ.2) THEN
288:             EEE(J1)=EEE(J1)+DUMMY283:             EEE(J1)=EEE(J1)+DUMMY
289:             REPEINT(J1)=REPEINT(J1)+DUMMY284:             REPEINT(J1)=REPEINT(J1)+DUMMY
290:             NREPINT(J1)=NREPINT(J1)+1285:             NREPINT(J1)=NREPINT(J1)+1
291:          ELSE IF (J1.LT.INTIMAGE+2) THEN286:          ELSE IF (J1.LT.INTIMAGE+2) THEN
292:             EEE(J1)=EEE(J1)+DUMMY/2.0D0287:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
293:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0288:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
294:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0289:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0
295:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0290:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
296:             NREPINT(J1)=NREPINT(J1)+1291:             NREPINT(J1)=NREPINT(J1)+1
298:          ELSE IF (J1.EQ.INTIMAGE+2) THEN293:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
299:             EEE(J1-1)=EEE(J1-1)+DUMMY294:             EEE(J1-1)=EEE(J1-1)+DUMMY
300:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY295:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
301:             NREPINT(J1-1)=NREPINT(J1-1)+1296:             NREPINT(J1-1)=NREPINT(J1-1)+1
302:          ENDIF297:          ENDIF
303:          EREP=EREP+DUMMY298:          EREP=EREP+DUMMY
304:          IF (DUMMY.GT.EMAX) THEN299:          IF (DUMMY.GT.EMAX) THEN
305:             IMAX=J1300:             IMAX=J1
306:             JMAX=J2301:             JMAX=J2
307:             EMAX=DUMMY302:             EMAX=DUMMY
 303:             EMAX=DUMMY
308:             DISTMAX=DINT304:             DISTMAX=DINT
309:             CUTMAX=NREPCUT(J2)305:             CUTMAX=NREPCUT(J2)
310:          ENDIF306:          ENDIF
311:          IF (DUMMY.GT.EEMAX(J1)) THEN307:          IF (DUMMY.GT.EEMAX(J1)) THEN
312:             JJMAX(J1)=J2308:             JJMAX(J1)=J2
313:             EEMAX(J1)=DUMMY309:             EEMAX(J1)=DUMMY
314:          ENDIF310:          ENDIF
315: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)311: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
316: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)312: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
317: !        DUMMY=-2.0D0*(1.0D0/(DINT*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2313: !        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))314:          DUMMY=-2.0D0*(DCUT/(DINT*D12)-1.0D0/NREPCUT(J2))
319:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)315:          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), &316: !        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)317: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
322: !318: !
323: ! Gradient contributions for image J1-1319: ! Gradient contributions for image J1-1
324: !320: !
325:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)321:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
326:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)322:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
327:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)323:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)
328: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &324: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
333:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)329:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
334:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)330:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
335:       ENDIF331:       ENDIF
336:    ENDDO332:    ENDDO
337: ENDDO333: ENDDO
338: IF (JMAX.GT.0) THEN334: IF (JMAX.GT.0) THEN
339:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G20.10,A,3G20.10)') ' congrad> Highest repulsive  contribution for any image in image ',IMAX, &  335:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G20.10,A,3G20.10)') ' congrad> Highest repulsive  contribution for any image in image ',IMAX, &  
340:  &  ' pair index ', &336:  &  ' pair index ', &
341:  &                                JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX),' value=',EMAX,' d,cutoff=',DISTMAX,CUTMAX 337:  &                                JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX),' value=',EMAX,' d,cutoff=',DISTMAX,CUTMAX 
342: ENDIF338: ENDIF
343: CONVERGEREPTEST=EMAX339: CONVERGEREPTEST=EMAX/INTCONSTRAINTDEL
344: 654 CONTINUE 
345: 340: 
346: !341: !
347: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise342: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
348: ! energy terms between images except for the end points.343: ! energy terms between images except for the end points.
349: !344: !
350: ESPRING=0.0D0345: ESPRING=0.0D0
351: IF (KINT.NE.0.0D0) THEN346: IF (KINT.NE.0.0D0) THEN
352:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1347:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1
353:       NI1=(3*NATOMS)*(J1-1)348:       NI1=(3*NATOMS)*(J1-1)
354:       NI2=(3*NATOMS)*J1349:       NI2=(3*NATOMS)*J1


r33499/intlbfgs.f90 2017-11-18 14:30:14.896918814 +0000 r33498/intlbfgs.f90 2017-11-18 14:30:15.340924678 +0000
659: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.659: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.
660: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!660: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!
661: !661: !
662: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.662: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.
663: !663: !
664: BESTWORST=1.0D100664: BESTWORST=1.0D100
665: 9876 CONTINUE665: 9876 CONTINUE
666: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)666: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
667: 667: 
668: 668: 
669: !                DIFF=1.0D-6669: !               DIFF=1.0D-7
670: !                PRINT*,'analytic and numerical gradients:'670: !               PRINT*,'analytic and numerical gradients:'
671: !                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)671: !               CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
672: !                DO J1=2,INTIMAGE+1672: !               DO J1=2,INTIMAGE+1
673: !                   DO J2=1,3*NATOMS673: !                  DO J2=1,3*NATOMS
674: !                      IF (.NOT.ATOMACTIVE((J2-1)/3+1)) CYCLE674: !                     J3=3*NATOMS*(J1-1)+J2
675: !                      J3=3*NATOMS*(J1-1)+J2675: !                     XYZ(J3)=XYZ(J3)+DIFF
676: !                      XYZ(J3)=XYZ(J3)+DIFF676: !                     CALL CONGRAD2(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)
677: !                      CALL CONGRAD2(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)677: !                     XYZ(J3)=XYZ(J3)-2.0D0*DIFF
678: !                      XYZ(J3)=XYZ(J3)-2.0D0*DIFF678: !                     CALL CONGRAD2(NMAXINT,NMININT,EMINUS,XYZ,VMINUS,EEE,IMGFREEZE,RMS)
679: !                      CALL CONGRAD2(NMAXINT,NMININT,EMINUS,XYZ,VMINUS,EEE,IMGFREEZE,RMS)679: !                     XYZ(J3)=XYZ(J3)+DIFF
680: !                      XYZ(J3)=XYZ(J3)+DIFF680: !                     IF ((ABS(GGG(J3)).NE.0.0D0).AND. &
681: !                      IF ((ABS(GGG(J3)).NE.0.0D0).AND. &681: !    &               (ABS(100.0D0*(GGG(J3)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GGG(J3)).GT.1.0D0)) THEN
682: !     &               (ABS(100.0D0*(GGG(J3)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GGG(J3)).GT.1.0D0)) THEN682: !                        WRITE(*,'(A,2I5,2F20.10,L5,A)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &
683: !                         WRITE(*,'(A,2I5,2G20.10,L5,A)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &683: !    & ATOMACTIVE((J2-1)/3+1),'   X'   
684: !     & ATOMACTIVE((J2-1)/3+1),'   X'   684: !                     ELSE
685: !                      ELSE685: !                        WRITE(*,'(A,2I5,2F20.10,L5)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &
686: !                         WRITE(*,'(A,2I5,2G20.10,L5)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &686: !    & ATOMACTIVE((J2-1)/3+1)  
687: !     & ATOMACTIVE((J2-1)/3+1)  687: !                     ENDIF
688: !                      ENDIF688: !                  ENDDO
689: !                   ENDDO689: !               ENDDO
690: !                ENDDO690: ! !
691: !  !691: !               STOP
692: !                STOP 
693: 692: 
694: 693: 
695: IF (QCIADDREP.GT.0) THEN694: IF (QCIADDREP.GT.0) THEN
696:    CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)695:    CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
697: !696: !
698: ! Don't do the CONINT part of CONGRAD2 if CONINT isn't set. CONGRAD seems to be697: ! Don't do the CONINT part of CONGRAD2 if CONINT isn't set. CONGRAD seems to be
699: ! dong something different at the moment. Focus on CONGRAD2698: ! dong something different at the moment. Focus on CONGRAD2
700: !699: !
701: ELSEIF (CHECKCONINT) THEN700: ELSEIF (CHECKCONINT) THEN
702:    CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)701:    CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1657:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1656:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1658:          ELSE1657:          ELSE
1659:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1658:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1660:          ENDIF1659:          ENDIF
1661:       ENDIF1660:       ENDIF
1662:       LASTGOODE=ETOTAL1661:       LASTGOODE=ETOTAL
1663:    ENDIF1662:    ENDIF
1664: 1663: 
1665:    EXITSTATUS=01664:    EXITSTATUS=0
1666:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1665:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1667:    PRINT *,'switched=',switched 
1668:    PRINT *,'MAXRMS,INTRMSTOL=',MAXRMS,INTRMSTOL 
1669:    PRINT *,'CONVERGECONTEST,CONVERGEREPTEST,MAXCONE=',CONVERGECONTEST,CONVERGEREPTEST,MAXCONE 
1670:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.(NITERDONE>1).AND.(CONVERGECONTEST.LT.MAXCONE).AND.(CONVERGEREPTEST.LT.MAXCONE)) EXITSTATUS=1 1666:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.(NITERDONE>1).AND.(CONVERGECONTEST.LT.MAXCONE).AND.(CONVERGEREPTEST.LT.MAXCONE)) EXITSTATUS=1 
1671:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 1667:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 
1672:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21668:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1673:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1669:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
1674:    PRINT *,'EXITSTATUS=',EXITSTATUS 
1675: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX1670: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX
1676: 1671: 
1677:    IF (EXITSTATUS > 0) THEN  1672:    IF (EXITSTATUS > 0) THEN  
1678:       PRINT *,'here A' 
1679:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1673:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1680: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771674: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1681: !        PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))1675: !        PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))
1682: !        IF (MAXEEE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771676: !        IF (MAXEEE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1683:       PRINT *,'here B nactive, natoms=',nactive,natoms 
1684:          IF (NACTIVE.LT.NATOMS) THEN 1677:          IF (NACTIVE.LT.NATOMS) THEN 
1685:             ADDATOM=.TRUE.1678:             ADDATOM=.TRUE.
1686:             GOTO 7771679:             GOTO 777
1687:          ENDIF1680:          ENDIF
1688:          CALL MYCPU_TIME(FTIME,.FALSE.)1681:          CALL MYCPU_TIME(FTIME,.FALSE.)
1689:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1682:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1690:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1683:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1691:          CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1684:          CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1692:          CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1685:          CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1693:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1686:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0