hdiff output

r33467/congrad.f90 2017-11-13 12:30:18.447371624 +0000 r33466/congrad.f90 2017-11-13 12:30:19.119380586 +0000
 10: ! 10: !
 11: !   You should have received a copy of the GNU General Public License 11: !   You should have received a copy of the GNU General Public License
 12: !   along with this program; if not, write to the Free Software 12: !   along with this program; if not, write to the Free Software
 13: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 13: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 14: ! 14: !
 15: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 15: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 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
 21: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG 21: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
 22: USE PORFUNCS 22: USE PORFUNCS
 23: IMPLICIT NONE 23: IMPLICIT NONE
 24:             24:            
 25: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF 25: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF
 26: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF 26: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX, EMAXNOFF
 27: INTEGER JJMAX(INTIMAGE+2) 27: INTEGER JJMAX(INTIMAGE+2)
 28: DOUBLE PRECISION  EEMAX(INTIMAGE+2) 28: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
 29: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1 29: 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) 30: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
 31: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI 31: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
 32: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2) 32: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
 33: LOGICAL NOINT 33: LOGICAL NOINT
 34: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL 34: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL
 35: LOGICAL IMGFREEZE(INTIMAGE) 35: LOGICAL IMGFREEZE(INTIMAGE)
 36: DOUBLE PRECISION DPLUS, SPGRAD(3), DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq 36: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3)
 37:  37: 
 38: EEE(1:INTIMAGE+2)=0.0D0 38: EEE(1:INTIMAGE+2)=0.0D0
 39: CONE(1:INTIMAGE+2)=0.0D0 39: CONE(1:INTIMAGE+2)=0.0D0
 40: REPE(1:INTIMAGE+2)=0.0D0 40: REPE(1:INTIMAGE+2)=0.0D0
 41: REPEINT(1:INTIMAGE+2)=0.0D0 41: REPEINT(1:INTIMAGE+2)=0.0D0
 42: NREPINT(1:INTIMAGE+2)=0 42: NREPINT(1:INTIMAGE+2)=0
 43: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 43: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0
 44: ECON=0.0D0; EREP=0.0D0 44: ECON=0.0D0; EREP=0.0D0
 45: NINTMIN=0 45: NINTMIN=0
 46: NINTMIN2=0 46: NINTMIN2=0
 47: MYUNIT=6 47: MYUNIT=6
 48:  48: 
 49: EMAX=-1.0D200 49: EMAX=-1.0D200
 50: EEMAX(1:INTIMAGE+2)=-1.0D200 50: EEMAX(1:INTIMAGE+2)=-1.0D200
 51: JJMAX(1:INTIMAGE+2)=-1 51: JJMAX(1:INTIMAGE+2)=-1
 52: JMAX=-1 52: JMAX=-1
 53: IMAX=-1 53: IMAX=-1
 54: EMAXNOFF=-1.0D200 54: EMAXNOFF=-1.0D200
 55: JMAXNOFF=-1 55: JMAXNOFF=-1
 56: IMAXNOFF=-1 56: IMAXNOFF=-1
  57: 
 57: ! 58: !
 58: !  Constraint energy and forces. 59: !  Constraint energy and forces.
 59: ! 60: !
 60: DO J2=1,NCONSTRAINT 61: DO J2=1,NCONSTRAINT
 61:    IF (.NOT.CONACTIVE(J2)) CYCLE 62:    IF (.NOT.CONACTIVE(J2)) CYCLE
 62:       CCLOCAL=CONCUTLOCAL(J2) 63:       CCLOCAL=CONCUTLOCAL(J2)
 63:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 64:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
 64:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 65:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
 65: ! 66: !
 66: ! For J1 we consider the line segment between image J1-1 and J1. 67: ! For J1 we consider the line segment between image J1-1 and J1.
 67: ! There are INTIMAGE+1 line segments in total, with an energy contribution 68: ! There are INTIMAGE+1 line segments in total, with an energy contribution
 68: ! and corresponding gradient terms for each.  69: ! and corresponding gradient terms for each. 
 69: ! A and B refer to atoms, 2 refers to image J1. 70: ! A and B refer to atoms, 2 refers to image J1.
 70: ! 71: !
 71:    DO J1=2,INTIMAGE+1 72: !  DO J1=2,INTIMAGE+1
 72: !  DO J1=1,INTIMAGE+2  ! checking for zero! 73:    DO J1=1,INTIMAGE+2  ! checking for zero!
 73:       IF (FREEZENODEST) THEN ! IMGFREEZE is not allocated otherwise! 
 74:          IF (J1.EQ.2) THEN 
 75:             IF (IMGFREEZE(1)) THEN 
 76: !              IF (J2.EQ.1) PRINT '(A)','J1=2 and IMGFREEZE(1)=T cycle' 
 77:                CYCLE 
 78:             ENDIF 
 79:          ELSE IF (J1.EQ.INTIMAGE+2) THEN 
 80:             IF (IMGFREEZE(INTIMAGE)) THEN 
 81: !              IF (J2.EQ.1) PRINT '(A)','J1=INTIMAGE+2 and IMGFREEZE(INTIMAGE)=T cycle' 
 82:                CYCLE 
 83:             ENDIF 
 84:          ELSE 
 85:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) THEN 
 86: !              IF (J2.EQ.1) PRINT '(A,I6,A)','J1=',J1,' IMGFREEZE(J1-2)=T and IMGFREEZE(J1-1)=T cycle' 
 87:                CYCLE 
 88:             ENDIF 
 89:          ENDIF 
 90:       ENDIF 
 91:       NI1=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1) 74:       NI1=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1)
 92:       NJ1=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1) 75:       NJ1=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1)
 93:       R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3) 76:       R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3)
 94:       R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3) 77:       R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3)
 95:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 78:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)
 96:       IF (ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL) THEN  79:       IF (ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL) THEN 
 97:          DUMMY=D2-CONDISTREFLOCAL(J2)   80:          DUMMY=D2-CONDISTREFLOCAL(J2)  
 98:          G2(1)=(R2AX-R2BX)/D2 81:          G2(1)=(R2AX-R2BX)/D2
 99:          G2(2)=(R2AY-R2BY)/D2 82:          G2(2)=(R2AY-R2BY)/D2
100:          G2(3)=(R2AZ-R2BZ)/D2 83:          G2(3)=(R2AZ-R2BZ)/D2
135:  & ' constraint ',JMAXNOFF, &118:  & ' constraint ',JMAXNOFF, &
136:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF119:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF
137: ELSEIF (JMAX.GT.0) THEN120: ELSEIF (JMAX.GT.0) THEN
138:    JMAXNOFF=JMAX121:    JMAXNOFF=JMAX
139:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'122:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'
140:    CONOFFTRIED(1:INTCONMAX)=.FALSE.123:    CONOFFTRIED(1:INTCONMAX)=.FALSE.
141:    NCONOFF=0124:    NCONOFF=0
142: ENDIF125: ENDIF
143: JMAXCON=JMAXNOFF126: JMAXCON=JMAXNOFF
144: 127: 
145: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes128: 
146: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes129: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
 130: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
147: 131: 
148: ! INTCONST=INTCONSTRAINREPCUT**13132: ! INTCONST=INTCONSTRAINREPCUT**13
149: 133: 
150: EMAX=-1.0D200134: EMAX=-1.0D200
151: EEMAX(1:INTIMAGE+2)=-1.0D200135: EEMAX(1:INTIMAGE+2)=-1.0D200
152: JJMAX(1:INTIMAGE+2)=-1136: JJMAX(1:INTIMAGE+2)=-1
153: JMAX=-1137: JMAX=-1
154: IMAX=-1138: IMAX=-1
155: 139: 
156: DO J2=1,NNREPULSIVE140: DO J2=1,NNREPULSIVE
160: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images144: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images
161:       IF (FREEZENODEST) THEN145:       IF (FREEZENODEST) THEN
162:          IF (J1.EQ.2) THEN146:          IF (J1.EQ.2) THEN
163:             IF (IMGFREEZE(1)) CYCLE147:             IF (IMGFREEZE(1)) CYCLE
164:          ELSE IF (J1.EQ.INTIMAGE+2) THEN148:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
165:             IF (IMGFREEZE(INTIMAGE)) CYCLE149:             IF (IMGFREEZE(INTIMAGE)) CYCLE
166:          ELSE150:          ELSE
167:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE151:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
168:          ENDIF152:          ENDIF
169:       ENDIF153:       ENDIF
170: !      IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN154:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
171: !!        WRITE(*, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)155: !        WRITE(*, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
172: !         STOP156:          STOP
173: !      ENDIF157:       ENDIF
174: 158: !     WRITE(*,'(A,2I8,6G20.10)') 'congrad> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)
175:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1) 
176:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)159:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
177:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1) 
178:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)160:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
179: 161:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
180:       G1(1:3)=XYZ(NI1+1:NI1+3)-XYZ(NJ1+1:NJ1+3)162:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
181:       G2(1:3)=XYZ(NI2+1:NI2+3)-XYZ(NJ2+1:NJ2+3)163:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2)
182: !164:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1
183: ! Squared distance between atoms A and B for theta=0 - distance in image 2165: !        D12=D2**12
184: !166:          D12=D2**2
185:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2 
186: ! 
187: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1 
188: ! 
189:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2 
190:       DCUT=NREPCUT(J2)**2 
191:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff 
192:       r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3) 
193: ! 
194: ! Is there an internal extremum? 
195: ! 
196:       r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b 
197:       IF (r1apr2bmr2amr1bsq.LT.1.0D-50) THEN 
198: !        D1=1.0D100; D2=1.0D100; 
199:          NOINT=.TRUE. 
200:          D1=SQRT(DSQ1) 
201:          D2=SQRT(DSQ2) 
202:          G2(1:3)=G2(1:3)/D2 
203:          G1(1:3)=G1(1:3)/D1 
204:       ELSE 
205:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq) 
206:       ENDIF 
207: ! 
208: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions. 
209: ! 
210:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1 
211: !        D12=DSQ2**6 
212:          D12=DSQ2 
213: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST) 
214: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)167: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
215:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)168:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
216:          EEE(J1)=EEE(J1)+DUMMY169:          EEE(J1)=EEE(J1)+DUMMY
217: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN170:          REPE(J1)=REPE(J1)+DUMMY
218: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &171:          EREP=EREP+DUMMY
219: ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY 
220: !        ENDIF 
221:          IF (DUMMY.GT.EMAX) THEN172:          IF (DUMMY.GT.EMAX) THEN
222:             IMAX=J1173:             IMAX=J1
223:             JMAX=J2174:             JMAX=J2
224:             EMAX=DUMMY175:             EMAX=DUMMY
225:          ENDIF176:          ENDIF
226:          IF (DUMMY.GT.EEMAX(J1)) THEN177:          IF (DUMMY.GT.EEMAX(J1)) THEN
227:             JJMAX(J1)=J2178:             JJMAX(J1)=J2
228:             EEMAX(J1)=DUMMY179:             EEMAX(J1)=DUMMY
229:          ENDIF180:          ENDIF
230:          REPE(J1)=REPE(J1)+DUMMY181: 
231:          EREP=EREP+DUMMY 
232: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)182: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
233:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)183:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
 184:          G2(1)=(R2AX-R2BX)/D2
 185:          G2(2)=(R2AY-R2BY)/D2
 186:          G2(3)=(R2AZ-R2BZ)/D2
234:          REPGRAD(1:3)=DUMMY*G2(1:3)187:          REPGRAD(1:3)=DUMMY*G2(1:3)
235: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), & 
236: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2) 
237:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)188:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
238:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)189:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
239:       ENDIF190:       ENDIF
 191: !     WRITE(MYUNIT,'(A,2I6,4G20.10)') 'J1,J2,D2,NREPCUT,EEE,REPE=',J1,J2,D2,NREPCUT(J2),EEE(J1),REPE(J1)
240: !192: !
241: ! For internal minima we are counting edges. 193: ! For internal minima we are counting edges. 
242: ! Edge J1 is between images J1-1 and J1, starting from J1=2.194: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
243: ! Energy contributions are shared evenly, except for195: ! Energy contributions are shared evenly, except for
244: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which196: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which
245: ! is assigned to image INTIMAGE+1. Gradients are set to zero for197: ! is assigned to image INTIMAGE+1. Gradients are set to zero for
246: ! the end images.198: ! the end images.
247: !199: !
248:       DUMMY=0.0D0200:       IF (J1.EQ.1) CYCLE
 201:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
 202:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
 203:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
 204:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
 205: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
 206:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-10) THEN
 207: !        WRITE(*, '(A,I6,A,2I6)') 'B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
 208: !        WRITE(*, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
 209: !        WRITE(*, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
 210: !        WRITE(*,'(A,7I10)') 'congrad> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)
 211: !        WRITE(*,'(A,7I10)') 'frames ',J1-1,J1
 212:       ELSE
 213:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
 214:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
 215:          IF (.NOT.NOINT) THEN
 216: !           WRITE(*,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), &
 217: ! &                        ' distance,cutoff=',DINT,NREPCUT(J2)
 218:             NINTMIN=NINTMIN+1
 219:          ENDIF
 220:       ENDIF
249:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN221:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
250:          NINTMIN2=NINTMIN2+1222:          NINTMIN2=NINTMIN2+1
251: !        D12=DSQI**6223: !        D12=DSQI**6
252:          D12=DSQI224:          D12=DSQI
253: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)225: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
254:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)226:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
255:          IF (J1.EQ.2) THEN227:          IF (J1.EQ.2) THEN
256:             EEE(J1)=EEE(J1)+DUMMY228:             EEE(J1)=EEE(J1)+DUMMY
257:             REPEINT(J1)=REPEINT(J1)+DUMMY229:             REPEINT(J1)=REPEINT(J1)+DUMMY
258:             NREPINT(J1)=NREPINT(J1)+1230:             NREPINT(J1)=NREPINT(J1)+1
340:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 312:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
341:                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))313:                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))
342:                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)314:                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)
343:                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)315:                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)
344:             ENDIF316:             ENDIF
345:          ENDDO317:          ENDDO
346: !     ENDIF318: !     ENDIF
347:    ENDDO319:    ENDDO
348: ENDIF320: ENDIF
349: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest spring  contribution for any image in image ',IMAX321: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest spring  contribution for any image in image ',IMAX
350: IF (DEBUG) WRITE(*, '(A,3G20.10)') 'congrad2> ECON,EREP,ESPRING=',ECON,EREP,ESPRING 
351: !322: !
352: ! Set gradients on frozen atoms to zero.323: ! Set gradients on frozen atoms to zero.
353: !324: !
354: IF (FREEZE) THEN325: IF (FREEZE) THEN
355:    DO J1=2,INTIMAGE+1  326:    DO J1=2,INTIMAGE+1  
356:       DO J2=1,NATOMS327:       DO J2=1,NATOMS
357:          IF (FROZEN(J2)) THEN328:          IF (FROZEN(J2)) THEN
358:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0329:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0
359:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0330:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0
360:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0331:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0
405: !  WRITE(*, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)376: !  WRITE(*, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1)
406:    IF (REPEINT(J1).LT.MININT) THEN377:    IF (REPEINT(J1).LT.MININT) THEN
407:       MININT=REPEINT(J1)378:       MININT=REPEINT(J1)
408:       NMININT=J1379:       NMININT=J1
409:    ENDIF380:    ENDIF
410:    IF (REPE(J1).GT.MAXINT) THEN381:    IF (REPE(J1).GT.MAXINT) THEN
411:       MAXINT=REPE(J1)382:       MAXINT=REPE(J1)
412:       NMAXINT=J1383:       NMAXINT=J1
413:    ENDIF384:    ENDIF
414: ENDDO385: ENDDO
415: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT386: IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') 'congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT
416: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> smallest internal energy=',MININT,' for image ',NMININT387: IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') 'congrad> smallest internal energy=',MININT,' for image ',NMININT
417: ! IF (DEBUG) WRITE(*, '(A,2I6)') ' congrad> number of internal minima=',NINTMIN,NINTMIN2388: IF (DEBUG) WRITE(*, '(A,2I6)') 'congrad> number of internal minima=',NINTMIN,NINTMIN2
418: 389: 
419: END SUBROUTINE CONGRAD390: END SUBROUTINE CONGRAD
420: 391: 
421: SUBROUTINE MINMAXD2(D2,D1,DINT,DSQ2,DSQ1,G1,G2,G1INT,G2INT,NOINT,DEBUG,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)392: SUBROUTINE MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
 393:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,DEBUG)
422: USE KEY, ONLY : CHECKCONINT394: USE KEY, ONLY : CHECKCONINT
423: IMPLICIT NONE395: IMPLICIT NONE
424: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT396: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT
425: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)397: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)
426: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq398: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq
427: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY, DUMMY2399: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY, DUMMY2
428: LOGICAL NOINT, DEBUG400: LOGICAL NOINT, DEBUG
429: 401: 
 402: G2(1)=r2ax - r2bx
 403: G2(2)=r2ay - r2by
 404: G2(3)=r2az - r2bz
 405: G1(1)=r1ax - r1bx
 406: G1(2)=r1ay - r1by
 407: G1(3)=r1az - r1bz
 408: r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3)
430: !409: !
431: ! Is there an internal extremum?410: ! Squared distance between atoms A and B for theta=0 - distance in image 2
432: ! THIS TEST IS NOT NEEDED IF CHECKCONINT IS FALSE. SHOULD SKIP.411: !
 412: DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2
 413: !
 414: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1
433: !415: !
434: ! IF (r1apr2bmr2amr1bsq.EQ.0.0D0) THEN ! now done in calling routine416: DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2
435: !417: !
436: DUMMY=(DSQ1-r1amr1bdr2amr2b)/r1apr2bmr2amr1bsq418: ! Is there an internal extremum?
 419: !
 420: r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b
 421: IF (r1apr2bmr2amr1bsq.EQ.0.0D0) THEN
 422:    DUMMY=2.0D0 ! just to skip the internal extremum part
 423: ELSE
 424:    DUMMY=(DSQ1-r1amr1bdr2amr2b)/r1apr2bmr2amr1bsq
 425: ENDIF
437: NOINT=.TRUE.426: NOINT=.TRUE.
438: IF ((DUMMY.GT.0.0D0).AND.(DUMMY.LT.1.0D0)) NOINT=.FALSE.427: IF ((DUMMY.GT.0.0D0).AND.(DUMMY.LT.1.0D0)) NOINT=.FALSE.
439: G1INT(1:3)=0.0D0428: G1INT(1:3)=0.0D0
440: G2INT(1:3)=0.0D0429: G2INT(1:3)=0.0D0
441: D2=SQRT(DSQ2)430: D2=SQRT(DSQ2)
442: D1=SQRT(DSQ1)431: D1=SQRT(DSQ1)
443: DSQI=1.0D10432: DSQI=1.0D10
444: DINT=1.0D10433: DINT=1.0D10
445: IF (.NOT.NOINT) THEN434: IF (CHECKCONINT.AND.(.NOT.NOINT)) THEN
446:    r1amr1bdr2amr2bsq=r1amr1bdr2amr2b**2435:    r1amr1bdr2amr2bsq=r1amr1bdr2amr2b**2
447:    DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2436:    DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2
448:    DSQI=MAX(-DUMMY2/r1apr2bmr2amr1bsq,0.0D0)437:    DSQI=MAX(-DUMMY2/r1apr2bmr2amr1bsq,0.0D0)
449:    DUMMY=r1apr2bmr2amr1bsq**2438:    DUMMY=r1apr2bmr2amr1bsq**2
450:    DINT=SQRT(DSQI)439:    DINT=SQRT(DSQI)
451:    IF (DINT.LE.0.0D0) THEN440:    IF (DINT.LE.0.0D0) THEN
452:       NOINT=.TRUE.441:       NOINT=.TRUE.
453:    ELSE442:    ELSE
454:      DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2443:      DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2
455:      DUMMY=DUMMY*DINT ! Convert derivatives of distance^2 to derivative of distance. 
456:      G1INT(1:3)= (DUMMY2*(G1(1:3) - G2(1:3)) + r1apr2bmr2amr1bsq*(G1(1:3)*DSQ2 -G2(1:3)*r1amr1bdr2amr2b))/DUMMY444:      G1INT(1:3)= (DUMMY2*(G1(1:3) - G2(1:3)) + r1apr2bmr2amr1bsq*(G1(1:3)*DSQ2 -G2(1:3)*r1amr1bdr2amr2b))/DUMMY
457:      G2INT(1:3)= (DUMMY2*(G2(1:3) - G1(1:3)) + r1apr2bmr2amr1bsq*(G2(1:3)*DSQ1 -G1(1:3)*r1amr1bdr2amr2b))/DUMMY445:      G2INT(1:3)= (DUMMY2*(G2(1:3) - G1(1:3)) + r1apr2bmr2amr1bsq*(G2(1:3)*DSQ1 -G1(1:3)*r1amr1bdr2amr2b))/DUMMY
458:    ENDIF446:    ENDIF
459: ENDIF447: ENDIF
460: !448: !
461: ! Convert derivatives of distance^2 to derivative of distance.449: ! Convert derivatives of distance^2 to derivative of distance.
462: ! We have cancelled a factor of two above and below!450: ! We have cancelled a factor of two above and below!
463: !451: !
464: G2(1:3)=G2(1:3)/D2452: G2(1:3)=G2(1:3)/D2
465: G1(1:3)=G1(1:3)/D1453: G1(1:3)=G1(1:3)/D1
466: ! IF (.NOT.NOINT) THEN454: IF (.NOT.NOINT) THEN
467: !    G1INT(1:3)=G1INT(1:3)/DINT455:    G1INT(1:3)=G1INT(1:3)/DINT
468: !    G2INT(1:3)=G2INT(1:3)/DINT456:    G2INT(1:3)=G2INT(1:3)/DINT
469: ! ENDIF457: ENDIF
470: 458: 
471: END SUBROUTINE MINMAXD2459: END SUBROUTINE MINMAXD2
472: 460: 
473: SUBROUTINE MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,DEBUG,INTCONSTRAINREPCUT,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)   461: SUBROUTINE MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
 462:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,DEBUG,INTCONSTRAINREPCUT)
474: IMPLICIT NONE463: IMPLICIT NONE
475: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT, DUMMY2464: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT, DUMMY2
476: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3),INTCONSTRAINREPCUT465: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3),INTCONSTRAINREPCUT
477: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq466: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq
478: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY467: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY
479: LOGICAL NOINT, DEBUG468: LOGICAL NOINT, DEBUG
480: 469: 
481: DUMMY=(DSQ1-r1amr1bdr2amr2b)/r1apr2bmr2amr1bsq470: G1(1)=r1ax - r1bx
 471: G1(2)=r1ay - r1by
 472: G1(3)=r1az - r1bz
 473: G2(1)=r2ax - r2bx
 474: G2(2)=r2ay - r2by
 475: G2(3)=r2az - r2bz
 476: r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3)
 477: !
 478: ! Squared distance between atoms A and B for theta=0 - distance in image 2
 479: !
 480: DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2
 481: !
 482: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1
 483: !
 484: DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2
 485: !
 486: ! Is there an internal extremum?
 487: !
 488: r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b
 489: IF (r1apr2bmr2amr1bsq.EQ.0.0D0) THEN
 490:    DUMMY=2.0D0 ! just to skip the internal solution
 491:    WRITE(*, '(A,G20.10)') 'r1apr2bmr2amr1bsq=',r1apr2bmr2amr1bsq
 492:    WRITE(*, '(A,3G20.10)') 'R1AX,R1AY,R1AZ=',R1AX,R1AY,R1AZ
 493:    WRITE(*, '(A,3G20.10)') 'R2AX,R2AY,R2AZ=',R2AX,R2AY,R2AZ
 494:    WRITE(*, '(A,3G20.10)') 'R1BX,R1BY,R1BZ=',R1BX,R1BY,R1BZ
 495:    WRITE(*, '(A,3G20.10)') 'R2BX,R2BY,R2BZ=',R2BX,R2BY,R2BZ
 496: ELSE
 497:    DUMMY=(DSQ1-r1amr1bdr2amr2b)/r1apr2bmr2amr1bsq
 498: ENDIF
482: NOINT=.TRUE.499: NOINT=.TRUE.
483: IF ((DUMMY.GT.0.0D0).AND.(DUMMY.LT.1.0D0)) NOINT=.FALSE.500: IF ((DUMMY.GT.0.0D0).AND.(DUMMY.LT.1.0D0)) NOINT=.FALSE.
484: G1INT(1:3)=0.0D0501: G1INT(1:3)=0.0D0
485: G2INT(1:3)=0.0D0502: G2INT(1:3)=0.0D0
486: D2=SQRT(DSQ2)503: D2=SQRT(DSQ2)
487: D1=SQRT(DSQ1)504: D1=SQRT(DSQ1)
488: DSQI=1.0D10505: DSQI=1.0D10
489: DINT=1.0D10506: DINT=1.0D10
490: IF (.NOT.NOINT) THEN507: IF (.NOT.NOINT) THEN
491:    r1amr1bdr2amr2bsq=r1amr1bdr2amr2b**2508:    r1amr1bdr2amr2bsq=r1amr1bdr2amr2b**2
492:    DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2509:    DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2
493:    DSQI=MAX(-DUMMY2/r1apr2bmr2amr1bsq,0.0D0)510:    DSQI=MAX(-DUMMY2/r1apr2bmr2amr1bsq,0.0D0)
494:    DUMMY=r1apr2bmr2amr1bsq**2511:    DUMMY=r1apr2bmr2amr1bsq**2
495:    DINT=SQRT(DSQI)512:    DINT=SQRT(DSQI)
496:    IF (DINT.LE.0.0D0) THEN513:    IF (DINT.LE.0.0D0) THEN
497:       NOINT=.TRUE.514:       NOINT=.TRUE.
498:    ELSEIF (DINT.LE.INTCONSTRAINREPCUT) THEN ! skip otherwise515:    ELSEIF (DINT.LE.INTCONSTRAINREPCUT) THEN ! skip otherwise
499:       DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2516:       DUMMY2=r1amr1bdr2amr2bsq - DSQ1*DSQ2
500:       DUMMY=DUMMY*DINT ! to convert derivatives of distance^2 to derivative of distance. 
501:       G1INT(1:3)= (DUMMY2*(G1(1:3) - G2(1:3)) + r1apr2bmr2amr1bsq*(G1(1:3)*DSQ2 -G2(1:3)*r1amr1bdr2amr2b))/DUMMY517:       G1INT(1:3)= (DUMMY2*(G1(1:3) - G2(1:3)) + r1apr2bmr2amr1bsq*(G1(1:3)*DSQ2 -G2(1:3)*r1amr1bdr2amr2b))/DUMMY
502:       G2INT(1:3)= (DUMMY2*(G2(1:3) - G1(1:3)) + r1apr2bmr2amr1bsq*(G2(1:3)*DSQ1 -G1(1:3)*r1amr1bdr2amr2b))/DUMMY518:       G2INT(1:3)= (DUMMY2*(G2(1:3) - G1(1:3)) + r1apr2bmr2amr1bsq*(G2(1:3)*DSQ1 -G1(1:3)*r1amr1bdr2amr2b))/DUMMY
503:    ENDIF519:    ENDIF
504: ENDIF520: ENDIF
505: !521: !
506: ! Convert derivatives of distance^2 to derivative of distance.522: ! Convert derivatives of distance^2 to derivative of distance.
507: ! We have cancelled a factor of two above and below!523: ! We have cancelled a factor of two above and below!
508: !524: !
509: G2(1:3)=G2(1:3)/D2525: G2(1:3)=G2(1:3)/D2
510: G1(1:3)=G1(1:3)/D1526: G1(1:3)=G1(1:3)/D1
 527: IF (.NOT.NOINT) THEN
 528: !  IF (DINT.EQ.0.0D0) THEN
 529: !     PRINT '(A,G20.10)','minmaxd2r> ERROR *** DINT=',DINT
 530: !     PRINT *,'original dummy=',((r1ax-r1bx)*(r1ax-r1bx-r2ax+r2bx)+ &
 531: ! &        (r1ay-r1by)*(r1ay-r1by-r2ay+r2by)+(r1az-r1bz)*(r1az-r1bz-r2az+r2bz))/r1apr2bmr2amr1bsq
 532: !     PRINT *,'r1amr1bdr2amr2b=',r1amr1bdr2amr2b
 533: !     PRINT *,'r1amr1bdr2amr2bsq=',r1amr1bdr2amr2bsq
 534: !     PRINT *,'DSQ1=',DSQ1
 535: !     PRINT *,'DSQ2=',DSQ2
 536: !     PRINT *,'DSQI=',DSQI
 537: !     PRINT *,'DUMMY=',DUMMY
 538: !     PRINT *,'G1INT=',G1INT(1:3)
 539: !     PRINT *,'G2INT=',G2INT(1:3)
 540: !     PRINT *,'R1AX,R1AY,R1AZ=',R1AX,R1AY,R1AZ
 541: !     PRINT *,'R2AX,R2AY,R2AZ=',R2AX,R2AY,R2AZ
 542: !     PRINT *,'R1BX,R1BY,R1BZ=',R1BX,R1BY,R1BZ
 543: !     PRINT *,'R2BX,R2BY,R2BZ=',R2BX,R2BY,R2BZ
 544: !     STOP
 545: !  ENDIF
 546:    G1INT(1:3)=G1INT(1:3)/DINT
 547:    G2INT(1:3)=G2INT(1:3)/DINT
 548: ENDIF
511: 549: 
512: END SUBROUTINE MINMAXD2R550: END SUBROUTINE MINMAXD2R
513: 551: 
514: !552: !
515: ! This version of congrad tests for an internal minimum in the553: ! This version of congrad tests for an internal minimum in the
516: ! constraint distances as well as the repulsions.554: ! constraint distances as well as the repulsions.
517: !555: !
518: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)556: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
519: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &557: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
520:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &558:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
529: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF567: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF
530: INTEGER JJMAX(INTIMAGE+2)568: INTEGER JJMAX(INTIMAGE+2)
531: DOUBLE PRECISION  EEMAX(INTIMAGE+2)569: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
532: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1570: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
533: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)571: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
534: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI572: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
535: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)573: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)
536: LOGICAL NOINT, LPRINT574: LOGICAL NOINT, LPRINT
537: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)575: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)
538: LOGICAL IMGFREEZE(INTIMAGE), PRINTE576: LOGICAL IMGFREEZE(INTIMAGE), PRINTE
539: DOUBLE PRECISION DPLUS, SPGRAD(3), CCLOCAL, DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq577: DOUBLE PRECISION DPLUS, SPGRAD(3), CCLOCAL
540: 578: 
541: PRINTE=.FALSE.579: PRINTE=.FALSE.
542: 111 CONTINUE580: 111 CONTINUE
543: 581: 
544: EEE(1:INTIMAGE+2)=0.0D0582: EEE(1:INTIMAGE+2)=0.0D0
545: CONE(1:INTIMAGE+2)=0.0D0583: CONE(1:INTIMAGE+2)=0.0D0
546: REPE(1:INTIMAGE+2)=0.0D0584: REPE(1:INTIMAGE+2)=0.0D0
547: NCONINT(1:INTIMAGE+2)=0585: NCONINT(1:INTIMAGE+2)=0
548: NREPINT(1:INTIMAGE+2)=0586: NREPINT(1:INTIMAGE+2)=0
549: REPEINT(1:INTIMAGE+2)=0.0D0587: REPEINT(1:INTIMAGE+2)=0.0D0
564: !602: !
565: EMAX=-1.0D200603: EMAX=-1.0D200
566: EMAXNOFF=-1.0D200604: EMAXNOFF=-1.0D200
567: EEMAX(1:INTIMAGE+2)=-1.0D200605: EEMAX(1:INTIMAGE+2)=-1.0D200
568: JJMAX(1:INTIMAGE+2)=-1606: JJMAX(1:INTIMAGE+2)=-1
569: JMAX=-1607: JMAX=-1
570: IMAX=-1608: IMAX=-1
571: JMAXNOFF=-1609: JMAXNOFF=-1
572: IMAXNOFF=-1610: IMAXNOFF=-1
573: DO J2=1,NCONSTRAINT611: DO J2=1,NCONSTRAINT
 612: !  IF (J2.EQ.5111) WRITE(*,'(A,I6,L5,2I6,2L5)') 'J2,CONACTIVE,CONI,CONJ,active i, active j=',J2,CONACTIVE(J2),CONI(J2),CONJ(J2), &
 613: ! &                    ATOMACTIVE(CONI(J2)),ATOMACTIVE(CONJ(J2))
 614: !  IF (J2.EQ.6378) WRITE(*,'(A,I6,L5,2I6,2L5)') 'J2,CONACTIVE,CONI,CONJ,active i, active j=',J2,CONACTIVE(J2),CONI(J2),CONJ(J2), &
 615: ! &                    ATOMACTIVE(CONI(J2)),ATOMACTIVE(CONJ(J2))
574:    IF (.NOT.CONACTIVE(J2)) CYCLE616:    IF (.NOT.CONACTIVE(J2)) CYCLE
575:    CCLOCAL=CONCUTLOCAL(J2)617:    CCLOCAL=CONCUTLOCAL(J2)
576:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS618:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
577:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)619:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
578: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG620: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG
579: !  IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN621:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN
580: !     WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)622:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)
581: !     STOP623:       STOP
582: !  ENDIF624:    ENDIF
583: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG625: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG
584:    DO J1=2,INTIMAGE+2626:    DO J1=2,INTIMAGE+2
585:       IF (FREEZENODEST) THEN ! IMGFREEZE is not allocated otherwise!627:       IF (FREEZENODEST) THEN ! IMGFREEZE is not allocated otherwise!
586:          IF (J1.EQ.2) THEN628:          IF (J1.EQ.2) THEN
587:             IF (IMGFREEZE(1)) THEN629:             IF (IMGFREEZE(1)) THEN
588: !              IF (J2.EQ.1) PRINT '(A)','J1=2 and IMGFREEZE(1)=T cycle'630: !              IF (J2.EQ.1) PRINT '(A)','J1=2 and IMGFREEZE(1)=T cycle'
589:                CYCLE631:                CYCLE
590:             ENDIF632:             ENDIF
591:          ELSE IF (J1.EQ.INTIMAGE+2) THEN633:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
592:             IF (IMGFREEZE(INTIMAGE)) THEN634:             IF (IMGFREEZE(INTIMAGE)) THEN
597:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) THEN639:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) THEN
598: !              IF (J2.EQ.1) PRINT '(A,I6,A)','J1=',J1,' IMGFREEZE(J1-2)=T and IMGFREEZE(J1-1)=T cycle'640: !              IF (J2.EQ.1) PRINT '(A,I6,A)','J1=',J1,' IMGFREEZE(J1-2)=T and IMGFREEZE(J1-1)=T cycle'
599:                CYCLE641:                CYCLE
600:             ENDIF642:             ENDIF
601:          ENDIF643:          ENDIF
602:       ENDIF644:       ENDIF
603:       NI1=(3*NATOMS)*(J1-2)+3*(CONI(J2)-1)645:       NI1=(3*NATOMS)*(J1-2)+3*(CONI(J2)-1)
604:       NI2=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1)646:       NI2=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1)
605:       NJ1=(3*NATOMS)*(J1-2)+3*(CONJ(J2)-1)647:       NJ1=(3*NATOMS)*(J1-2)+3*(CONJ(J2)-1)
606:       NJ2=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1)648:       NJ2=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1)
607: 649:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
608:       G1(1:3)=XYZ(NI1+1:NI1+3)-XYZ(NJ1+1:NJ1+3)650:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
609:       G2(1:3)=XYZ(NI2+1:NI2+3)-XYZ(NJ2+1:NJ2+3)651:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
610: !652:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
611: ! Squared distance between atoms A and B for theta=0 - distance in image 2653:       CALL MINMAXD2(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
612: !654:   &                 D2,D1,DINT,G1,G2,G1INT,G2INT,NOINT,.FALSE.)
613:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2 
614: ! 
615: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1 
616: ! 
617:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2 
618:       r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3) 
619: ! 
620: ! Is there an internal extremum? 
621: ! 
622:       r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b 
623:  
624:       IF ((.NOT.CHECKCONINT).OR.(r1apr2bmr2amr1bsq.LT.1.0D-50)) THEN 
625:          NOINT=.TRUE. 
626:          D1=SQRT(DSQ1) 
627:          D2=SQRT(DSQ2) 
628:          G2(1:3)=G2(1:3)/D2 
629:          G1(1:3)=G1(1:3)/D1 
630:       ELSE 
631:          CALL MINMAXD2(D2,D1,DINT,DSQ2,DSQ1,G1,G2,G1INT,G2INT,NOINT,DEBUG,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq) 
632:       ENDIF 
633: !655: !
634: ! Need to include both D2 and D1 contributions if they are both outside tolerance.656: ! Need to include both D2 and D1 contributions if they are both outside tolerance.
635: ! Otherwise we get discontinuities if they are very close and swap over.657: ! Otherwise we get discontinuities if they are very close and swap over.
636: !658: !
637: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)659: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)
638: !660: !
639: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.661: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.
640: !662: !
641: !     IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &663:       IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &
642: ! &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL664:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL
643:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 665:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 
644:          DUMMY=D2-CONDISTREFLOCAL(J2)  666:          DUMMY=D2-CONDISTREFLOCAL(J2)  
645:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)667:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
646:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)668:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
647:          IF (DUMMY.GT.EMAX) THEN669:          IF (DUMMY.GT.EMAX) THEN
648:             IMAX=J1670:             IMAX=J1
649:             JMAX=J2671:             JMAX=J2
650:             EMAX=DUMMY672:             EMAX=DUMMY
651:          ENDIF673:          ENDIF
652:          IF (DUMMY.GT.EMAXNOFF) THEN674:          IF (DUMMY.GT.EMAXNOFF) THEN
661:             EEMAX(J1)=DUMMY683:             EEMAX(J1)=DUMMY
662:          ENDIF684:          ENDIF
663:          EEE(J1)=EEE(J1)+DUMMY685:          EEE(J1)=EEE(J1)+DUMMY
664:          CONE(J1)=CONE(J1)+DUMMY686:          CONE(J1)=CONE(J1)+DUMMY
665:          ECON=ECON      +DUMMY687:          ECON=ECON      +DUMMY
666:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &688:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
667:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)689:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
668:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)690:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
669:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)691:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
670:       ENDIF692:       ENDIF
671: ! ! 
672: ! ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1. 
673: ! ! 
674: ! ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1. 
675: ! ! 
676: ! !     IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') & 
677: ! ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL 
678: !       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN   
679: !          DUMMY=D1-CONDISTREFLOCAL(J2)   
680: !          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3) 
681: !          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2) 
682: ! !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN 
683: ! !           WRITE(*, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', & 
684: ! ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY 
685: ! !        ENDIF 
686: !          IF (DUMMY.GT.EMAX) THEN 
687: !             IMAX=J1 
688: !             JMAX=J2 
689: !             EMAX=DUMMY 
690: !          ENDIF 
691: !          IF (DUMMY.GT.EMAXNOFF) THEN 
692: !             IF (.NOT.CONOFFTRIED(J2)) THEN 
693: !                IMAXNOFF=J1 
694: !                JMAXNOFF=J2 
695: !                EMAXNOFF=DUMMY 
696: !             ENDIF 
697: !          ENDIF 
698: !          IF (DUMMY.GT.EEMAX(J1-1)) THEN 
699: !             JJMAX(J1-1)=J2 
700: !             EEMAX(J1-1)=DUMMY 
701: !          ENDIF 
702: !          EEE(J1-1)=EEE(J1-1)+DUMMY 
703: !          CONE(J1-1)=CONE(J1-1)+DUMMY 
704: !          ECON=ECON      +DUMMY 
705: !          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), & 
706: !   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2) 
707: !          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 
708: !          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 
709: !       ENDIF 
710: !693: !
 694: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.
 695: !
 696: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.
 697: !
 698: !     IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &
 699: ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL
 700:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  
 701:          DUMMY=D1-CONDISTREFLOCAL(J2)  
 702:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)
 703:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
 704: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
 705: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
 706: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
 707: !        ENDIF
 708:          IF (DUMMY.GT.EMAX) THEN
 709:             IMAX=J1
 710:             JMAX=J2
 711:             EMAX=DUMMY
 712:          ENDIF
 713:          IF (DUMMY.GT.EMAXNOFF) THEN
 714:             IF (.NOT.CONOFFTRIED(J2)) THEN
 715:                IMAXNOFF=J1
 716:                JMAXNOFF=J2
 717:                EMAXNOFF=DUMMY
 718:             ENDIF
 719:          ENDIF
 720:          IF (DUMMY.GT.EEMAX(J1-1)) THEN
 721:             JJMAX(J1-1)=J2
 722:             EEMAX(J1-1)=DUMMY
 723:          ENDIF
 724:          EEE(J1-1)=EEE(J1-1)+DUMMY
 725:          CONE(J1-1)=CONE(J1-1)+DUMMY
 726:          ECON=ECON      +DUMMY
 727:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
 728:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
 729:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
 730:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
 731:       ENDIF
 732: !     WRITE(*,'(A,2L5,2G20.10)') 'before intmin block EEE(J1),EEE(J1-1)=',CHECKCONINT,NOINT,EEE(J1),EEE(J1-1)
711:       IF (CHECKCONINT.AND.(.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN733:       IF (CHECKCONINT.AND.(.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN
712:          DUMMY=DINT-CONDISTREFLOCAL(J2)  734:          DUMMY=DINT-CONDISTREFLOCAL(J2)  
713:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)735:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)
714:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)736:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
715:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)737:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
716:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)738:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)
717:          DUMMY=INTMINFAC*INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)739:          DUMMY=INTMINFAC*INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
718: !        WRITE(*,'(A,3I7,9F13.5)') 'J1,NI1,NJ1,INTMINFAC,INTCONSTRAINTDEL,DUMMY,GGG=',J1,NI1,NJ1,INTMINFAC,INTCONSTRAINTDEL, &740: !        WRITE(*,'(A,3I7,9F13.5)') 'J1,NI1,NJ1,INTMINFAC,INTCONSTRAINTDEL,DUMMY,GGG=',J1,NI1,NJ1,INTMINFAC,INTCONSTRAINTDEL, &
719: ! &            DUMMY, GGG(NI1+1:NI1+3),GGG(NJ1+1:NJ1+3)   741: ! &            DUMMY, GGG(NI1+1:NI1+3),GGG(NJ1+1:NJ1+3)   
720: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN742: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
761: !        WRITE(*, '(A,4I6,G15.5)') 'in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &783: !        WRITE(*, '(A,4I6,G15.5)') 'in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
762: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)784: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
763:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)785:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
764:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)786:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
765: !        WRITE(*, '(A,3I7,9G13.5)') 'J1,NI2,NJ2,INTMINFAC,INTCONSTRAINTDEL,DUMMY,GGG=',J1,NI2,NJ2,INTMINFAC,INTCONSTRAINTDEL, & 787: !        WRITE(*, '(A,3I7,9G13.5)') 'J1,NI2,NJ2,INTMINFAC,INTCONSTRAINTDEL,DUMMY,GGG=',J1,NI2,NJ2,INTMINFAC,INTCONSTRAINTDEL, & 
766: ! &            DUMMY, GGG(NI2+1:NI2+3),GGG(NJ2+1:NJ2+3)   788: ! &            DUMMY, GGG(NI2+1:NI2+3),GGG(NJ2+1:NJ2+3)   
767: !        WRITE(*,'(A,2G20.10)') 'in intmin block EEE(J1),EEE(J1-1)=',EEE(J1),EEE(J1-1)789: !        WRITE(*,'(A,2G20.10)') 'in intmin block EEE(J1),EEE(J1-1)=',EEE(J1),EEE(J1-1)
768:       ENDIF790:       ENDIF
769:    ENDDO791:    ENDDO
770: ENDDO792: ENDDO
 793: ! DO J2=1,INTIMAGE+2
 794: !    IF (JJMAX(J2).GT.0) THEN
 795: !       WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for image ',J2,' constraint ',JJMAX(J2),' atoms ', &
 796: !  &                                CONI(JJMAX(J2)),CONJ(JJMAX(J2))
 797: !    ENDIF
 798: ! ENDDO
771: IF (JMAX.GT.0) THEN799: IF (JMAX.GT.0) THEN
772:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for any image in image ',IMAX, &800:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for any image in image ',IMAX, &
773:  & ' constraint ',JMAX, &801:  & ' constraint ',JMAX, &
774:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)802:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)
775: ENDIF803: ENDIF
776: IF (JMAXNOFF.GT.0) THEN804: IF (JMAXNOFF.GT.0) THEN
777:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &805:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &
778:  & ' constraint ',JMAXNOFF, &806:  & ' constraint ',JMAXNOFF, &
779:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF807:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF
780: ELSEIF (JMAX.GT.0) THEN808: ELSEIF (JMAX.GT.0) THEN
798:    DO J1=2,INTIMAGE+2826:    DO J1=2,INTIMAGE+2
799:       IF (FREEZENODEST) THEN827:       IF (FREEZENODEST) THEN
800:          IF (J1.EQ.2) THEN828:          IF (J1.EQ.2) THEN
801:             IF (IMGFREEZE(1)) CYCLE829:             IF (IMGFREEZE(1)) CYCLE
802:          ELSE IF (J1.EQ.INTIMAGE+2) THEN830:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
803:             IF (IMGFREEZE(INTIMAGE)) CYCLE831:             IF (IMGFREEZE(INTIMAGE)) CYCLE
804:          ELSE832:          ELSE
805:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE833:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
806:          ENDIF834:          ENDIF
807:       ENDIF835:       ENDIF
808: !     IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN836:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
809: !        WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)837:          WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
810: !        STOP838:          STOP
811: !     ENDIF839:       ENDIF
812:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)840:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
813:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)841:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
814:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)842:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
815:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)843:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
816: 844:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
817:       G1(1:3)=XYZ(NI1+1:NI1+3)-XYZ(NJ1+1:NJ1+3) 845:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
818:       G2(1:3)=XYZ(NI2+1:NI2+3)-XYZ(NJ2+1:NJ2+3) 846:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
819: !847:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
820: ! Squared distance between atoms A and B for theta=0 - distance in image 2848: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
821: !849:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-50) THEN
822:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2850: !        WRITE(*, '(A,I6,A,2I6)') 'A repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
823: !851: !        WRITE(*, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
824: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1852: !        WRITE(*, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
825: !853: !        WRITE(*,'(A,7I10)') 'congrad2> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)
826:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2854: !        WRITE(*,'(A,7I10)') 'frames ',J1-1,J1
827:       DCUT=NREPCUT(J2)**2855:          D1=1.0D100; D2=1.0D100; NOINT=.TRUE.  ! to skip the next blocks
828:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff 
829:       r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3) 
830: ! 
831: ! Is there an internal extremum? 
832: ! 
833:       r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b 
834:  
835:       IF (r1apr2bmr2amr1bsq.LT.1.0D-50) THEN 
836: !        D1=1.0D100; D2=1.0D100;  
837:          NOINT=.TRUE.   
838:          D1=SQRT(DSQ1) 
839:          D2=SQRT(DSQ2) 
840:          G2(1:3)=G2(1:3)/D2 
841:          G1(1:3)=G1(1:3)/D1 
842:       ELSE856:       ELSE
843:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)857:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
 858:   &                    D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
844:       ENDIF859:       ENDIF
845: !?????????????????????????????????????????????????????????????????????????????? 
846: ! 
847: !  WHY ARE WE DOING both D2 and D1? Isn't this double counting? See CONGRAD routine above 
848: ! 
849: !860: !
850: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.861: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
851: !862: !
 863: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
 864: 
 865: !     IF ((NREPJ(J2).EQ.2024).OR.(NREPI(J2).EQ.2024)) WRITE(*, '(A,2I6,3G20.10,I6)') &
 866: ! &    'j,i,D2,NREPCUT,INTCONST,DUMMY,J1=',NREPJ(J2),NREPI(J2),D2,NREPCUT(J2),INTCONST,J1
 867: 
852:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1868:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
853: !        D12=DSQ2**6869: !        D12=DSQ2**6
854:          D12=DSQ2870:          D12=DSQ2
855: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)871: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
856: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)872: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
857:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)873:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
 874: !        IF ((NREPJ(J2).EQ.2024).OR.(NREPI(J2).EQ.2024)) WRITE(*, '(A,2I6,6G20.10)') &
 875: ! &        'j,i,INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',NREPJ(J2),NREPI(J2),INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 
 876: !        IF ((NREPJ(J2).EQ.83).AND.(NREPI(J2).EQ.357)) WRITE(*, '(A,6G20.10)') &
 877: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 
 878: !        IF ((NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.357)) WRITE(*, '(A,6G20.10)') &
 879: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 
858:          EEE(J1)=EEE(J1)+DUMMY880:          EEE(J1)=EEE(J1)+DUMMY
859: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN881:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
860: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &882:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
861: ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY883:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
862: !        ENDIF884:          ENDIF
863:          IF (DUMMY.GT.EMAX) THEN885:          IF (DUMMY.GT.EMAX) THEN
864:             IMAX=J1886:             IMAX=J1
865:             JMAX=J2887:             JMAX=J2
866:             EMAX=DUMMY888:             EMAX=DUMMY
867:          ENDIF889:          ENDIF
868:          IF (DUMMY.GT.EEMAX(J1)) THEN890:         IF (DUMMY.GT.EEMAX(J1)) THEN
869:             JJMAX(J1)=J2891:             JJMAX(J1)=J2
870:             EEMAX(J1)=DUMMY892:             EEMAX(J1)=DUMMY
871:          ENDIF893:          ENDIF
872:          REPE(J1)=REPE(J1)+DUMMY894:          REPE(J1)=REPE(J1)+DUMMY
873:          EREP=EREP+DUMMY895:          EREP=EREP+DUMMY
874: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)896: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
875:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)897:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
876:          REPGRAD(1:3)=DUMMY*G2(1:3)898:          REPGRAD(1:3)=DUMMY*G2(1:3)
877: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &899: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
878: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)900: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
879:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)901:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
880:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)902:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
881:       ENDIF903:       ENDIF
882:       DUMMY=0.0D0904:       DUMMY=0.0D0
883: !905: !
884: ! SURELY THIS BLOCK JUST DUPLICATES THE TERMS FOR THE D2 BLOCK?906: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.
885: !907: !
886: ! !908: !     IF ((D1.LT.INTCONSTRAINREPCUT).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1
887: ! ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.909:       IF ((D1.LT.NREPCUT(J2)).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1
888: ! !910: !        D12=DSQ1**6
889: ! !     IF ((D1.LT.INTCONSTRAINREPCUT).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1911:          D12=DSQ1
890: !       IF ((D1.LT.NREPCUT(J2)).AND.(J1.GT.2)) THEN ! terms for image J1-1 - non-zero derivatives only for J1-1912: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
891: ! !        D12=DSQ1**6913: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*NREPCUT(J2))/INTCONST)
892: !          D12=DSQ1914:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D1-3.0D0*NREPCUT(J2))/INTCONST)
893: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*INTCONSTRAINREPCUT)/INTCONST)915:          EEE(J1-1)=EEE(J1-1)+DUMMY
894: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*NREPCUT(J2))/INTCONST)916:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
895: !          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D1-3.0D0*NREPCUT(J2))/INTCONST)917:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R2 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
896: !          EEE(J1-1)=EEE(J1-1)+DUMMY918:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
897: ! !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN919:          ENDIF
898: ! !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R2 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &920:          IF (DUMMY.GT.EMAX) THEN
899: ! ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY921:             IMAX=J1
900: ! !        ENDIF922:             JMAX=J2
901: !          IF (DUMMY.GT.EMAX) THEN923:             EMAX=DUMMY
902: !             IMAX=J1924:          ENDIF
903: !             JMAX=J2925:          IF (DUMMY.GT.EEMAX(J1-1)) THEN
904: !             EMAX=DUMMY926:             JJMAX(J1-1)=J2
905: !          ENDIF927:             EEMAX(J1-1)=DUMMY
906: !          IF (DUMMY.GT.EEMAX(J1-1)) THEN928:          ENDIF
907: !             JJMAX(J1-1)=J2929:          REPE(J1-1)=REPE(J1-1)+DUMMY
908: !             EEMAX(J1-1)=DUMMY930:          EREP=EREP+DUMMY
909: !          ENDIF931: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)
910: !          REPE(J1-1)=REPE(J1-1)+DUMMY932:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)
911: !          EREP=EREP+DUMMY933:          REPGRAD(1:3)=DUMMY*G1(1:3)
912: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)934: !        WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
913: !          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)935: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
914: !          REPGRAD(1:3)=DUMMY*G1(1:3)936:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
915: ! !        WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &937:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
916: ! ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)938:       ENDIF
917: !          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 
918: !          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 
919: !       ENDIF 
920:       DUMMY=0.0D0939:       DUMMY=0.0D0
921: !     IF ((.NOT.NOINT).AND.(DINT.LT.INTCONSTRAINREPCUT)) THEN940: !     IF ((.NOT.NOINT).AND.(DINT.LT.INTCONSTRAINREPCUT)) THEN
922:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN941:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
923: !        D12=DSQI**6942: !        D12=DSQI**6
924:          D12=DSQI943:          D12=DSQI
925: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*INTCONSTRAINREPCUT)/INTCONST)944: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
926: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)945: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
927:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)946:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
928:          EREP=EREP+DUMMY947:          EREP=EREP+DUMMY
929: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN948: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &
930: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &949: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY
931: ! &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY 
932: !        ENDIF950: !        ENDIF
 951:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
 952:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
 953:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
 954:          ENDIF
933:          IF (DUMMY.GT.EMAX) THEN955:          IF (DUMMY.GT.EMAX) THEN
934:             IMAX=J1956:             IMAX=J1
935:             JMAX=J2957:             JMAX=J2
936:             EMAX=DUMMY958:             EMAX=DUMMY
937:          ENDIF959:          ENDIF
938:          IF (DUMMY.GT.EEMAX(J1)) THEN960:          IF (DUMMY.GT.EEMAX(J1)) THEN
939:             JJMAX(J1)=J2961:             JJMAX(J1)=J2
940:             EEMAX(J1)=DUMMY962:             EEMAX(J1)=DUMMY
941:          ENDIF963:          ENDIF
942:          IF (DUMMY.GT.EEMAX(J1-1)) THEN964:          IF (DUMMY.GT.EEMAX(J1-1)) THEN
955:             NREPINT(J1)=NREPINT(J1)+1977:             NREPINT(J1)=NREPINT(J1)+1
956:             NREPINT(J1-1)=NREPINT(J1-1)+1978:             NREPINT(J1-1)=NREPINT(J1-1)+1
957:          ELSE IF (J1.EQ.INTIMAGE+2) THEN979:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
958:             EEE(J1-1)=EEE(J1-1)+DUMMY980:             EEE(J1-1)=EEE(J1-1)+DUMMY
959:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY981:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY
960:             NREPINT(J1-1)=NREPINT(J1-1)+1982:             NREPINT(J1-1)=NREPINT(J1-1)+1
961:          ENDIF983:          ENDIF
962: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)984: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
963:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)985:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
964:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)986:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)
 987: !        WRITE(*, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
 988: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
965:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)989:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
966:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)990:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
967:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)991:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)
 992: !        WRITE(*, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
 993: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
968:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)994:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
969:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)995:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
970:       ENDIF996:       ENDIF
971:    ENDDO997:    ENDDO
972: ENDDO998: ENDDO
 999: ! DO J2=1,INTIMAGE+2
 1000: !    IF (JJMAX(J2).GT.0) THEN
 1001: !       WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest repulsive  contribution for image ',J2,' pair index ', &
 1002: !  &                                JJMAX(J2),' atoms ', &
 1003: !  &                                NREPI(JJMAX(J2)),NREPJ(JJMAX(J2))
 1004: !    ENDIF
 1005: ! ENDDO
973: IF (JMAX.GT.0) THEN1006: IF (JMAX.GT.0) THEN
974:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest repulsive  contribution for any image in image ',IMAX, &1007:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest repulsive  contribution for any image in image ',IMAX, &
975:  &  ' pair index ', &1008:  &  ' pair index ', &
976:  &                                JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX)1009:  &                                JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX)
977: ENDIF1010: ENDIF
978: !1011: !
979: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise1012: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
980: ! energy terms between images except for the end points.1013: ! energy terms between images except for the end points.
981: !1014: !
982: ESPRING=0.0D01015: ESPRING=0.0D0
1018:          DO J2=1,NATOMS1051:          DO J2=1,NATOMS
1019:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 1052:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
1020:                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))1053:                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))
1021:                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)1054:                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)
1022:                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)1055:                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)
1023:             ENDIF1056:             ENDIF
1024:          ENDDO1057:          ENDDO
1025: !     ENDIF1058: !     ENDIF
1026:    ENDDO1059:    ENDDO
1027: ENDIF1060: ENDIF
1028: ! IF (KINTENDS.GT.0.0D0) THEN1061: IF (KINTENDS.GT.0.0D0) THEN
1029: ! !1062: !
1030: ! ! Extra terms for the two fixed end points.1063: ! Extra terms for the two fixed end points.
1031: ! !1064: !
1032: !    DO J1=2,INTIMAGE+1 ! sum over images1065:    DO J1=2,INTIMAGE+1 ! sum over images
1033: !       NI1=(3*NATOMS)*(J1-1)1066:       NI1=(3*NATOMS)*(J1-1)
1034: !       NI2=(3*NATOMS)*(INTIMAGE+1)1067:       NI2=(3*NATOMS)*(INTIMAGE+1)
1035: ! !1068: !
1036: ! !  Spring between 1 and J11069: !  Spring between 1 and J1
1037: ! !1070: !
1038: !       DPLUS=0.0D01071:       DPLUS=0.0D0
1039: !       DO J2=1,NATOMS1072:       DO J2=1,NATOMS
1040: !          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN1073:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN
1041: !             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(3*(J2-1)+1))**2 &1074:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(3*(J2-1)+1))**2 &
1042: !   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(3*(J2-1)+2))**2 &1075:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(3*(J2-1)+2))**2 &
1043: !   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(3*(J2-1)+3))**21076:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(3*(J2-1)+3))**2
1044: !          ENDIF1077:          ENDIF
1045: ! !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS1078: !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS
1046: !       ENDDO1079:       ENDDO
1047: !       DPLUS=SQRT(DPLUS)1080:       DPLUS=SQRT(DPLUS)
1048: ! !     IF (DPLUS.GT.IMSEPMAX) THEN1081: !     IF (DPLUS.GT.IMSEPMAX) THEN
1049: ! !        DUMMY=KINTENDS*0.5D0*(DPLUS-IMSEPMAX)**21082: !        DUMMY=KINTENDS*0.5D0*(DPLUS-IMSEPMAX)**2
1050: !          DUMMY=KINTENDS*0.5D0*DPLUS**21083:          DUMMY=KINTENDS*0.5D0*DPLUS**2
1051: !          DUMMY=DUMMY*(1.0D0/(J1-1))**2 ! (INTIMAGE/(J1-1))**21084:          DUMMY=DUMMY*(1.0D0/(J1-1))**2 ! (INTIMAGE/(J1-1))**2
1052: !          IF (DUMMY.GT.EMAX) THEN1085:          IF (DUMMY.GT.EMAX) THEN
1053: !             IMAX=J11086:             IMAX=J1
1054: !             EMAX=DUMMY1087:             EMAX=DUMMY
1055: !          ENDIF1088:          ENDIF
1056: ! !        DUMMY=KINTENDS*0.5D0*DPLUS**21089: !        DUMMY=KINTENDS*0.5D0*DPLUS**2
1057: !          ESPRING=ESPRING+DUMMY1090:          ESPRING=ESPRING+DUMMY
1058: !          IF (DUMMY.GT.EEMAX(J1)) THEN1091:          IF (DUMMY.GT.EEMAX(J1)) THEN
1059: !             EEMAX(J1)=DUMMY1092:             EEMAX(J1)=DUMMY
1060: !          ENDIF1093:          ENDIF
1061: ! 1094: 
1062: ! !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING1095: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING
1063: ! !        DUMMY=KINTENDS*(DPLUS-IMSEPMAX)/DPLUS1096: !        DUMMY=KINTENDS*(DPLUS-IMSEPMAX)/DPLUS
1064: !          DUMMY=KINTENDS1097:          DUMMY=KINTENDS
1065: !          DO J2=1,NATOMS1098:          DO J2=1,NATOMS
1066: !             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN1099:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN
1067: !                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(3*(J2-1)+1:3*(J2-1)+3))1100:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(3*(J2-1)+1:3*(J2-1)+3))
1068: !                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)1101:                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)
1069: !             ENDIF1102:             ENDIF
1070: !          ENDDO1103:          ENDDO
1071: ! !     ENDIF1104: !     ENDIF
1072: ! !1105: !
1073: ! !  Spring between INTIMAGE+2 and J11106: !  Spring between INTIMAGE+2 and J1
1074: ! !1107: !
1075: !       DPLUS=0.0D01108:       DPLUS=0.0D0
1076: !       DO J2=1,NATOMS1109:       DO J2=1,NATOMS
1077: !          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN1110:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN
1078: !             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &1111:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &
1079: !   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &1112:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &
1080: !   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**21113:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2
1081: !          ENDIF1114:          ENDIF
1082: ! !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS1115: !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS
1083: !       ENDDO1116:       ENDDO
1084: !       DPLUS=SQRT(DPLUS)1117:       DPLUS=SQRT(DPLUS)
1085: ! !     IF (DPLUS.GT.IMSEPMAX) THEN1118: !     IF (DPLUS.GT.IMSEPMAX) THEN
1086: ! !        DUMMY=KINTENDS*0.5D0*(DPLUS-IMSEPMAX)**21119: !        DUMMY=KINTENDS*0.5D0*(DPLUS-IMSEPMAX)**2
1087: !          DUMMY=KINTENDS*0.5D0*DPLUS**21120:          DUMMY=KINTENDS*0.5D0*DPLUS**2
1088: !          DUMMY=DUMMY*(INTIMAGE/(INTIMAGE+2-J1))**21121:          DUMMY=DUMMY*(INTIMAGE/(INTIMAGE+2-J1))**2
1089: !          IF (DUMMY.GT.EMAX) THEN1122:          IF (DUMMY.GT.EMAX) THEN
1090: !             IMAX=J11123:             IMAX=J1
1091: !             EMAX=DUMMY1124:             EMAX=DUMMY
1092: !          ENDIF1125:          ENDIF
1093: ! !        DUMMY=KINTENDS*0.5D0*DPLUS**21126: !        DUMMY=KINTENDS*0.5D0*DPLUS**2
1094: !          ESPRING=ESPRING+DUMMY1127:          ESPRING=ESPRING+DUMMY
1095: !          IF (DUMMY.GT.EEMAX(J1)) THEN1128:          IF (DUMMY.GT.EEMAX(J1)) THEN
1096: !             EEMAX(J1)=DUMMY1129:             EEMAX(J1)=DUMMY
1097: !          ENDIF1130:          ENDIF
1098: ! 1131: 
1099: ! !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING1132: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING
1100: ! !        DUMMY=KINTENDS*(DPLUS-IMSEPMAX)/DPLUS1133: !        DUMMY=KINTENDS*(DPLUS-IMSEPMAX)/DPLUS
1101: !          DUMMY=KINTENDS1134:          DUMMY=KINTENDS
1102: !          DO J2=1,NATOMS1135:          DO J2=1,NATOMS
1103: !             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN1136:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN
1104: !                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))1137:                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))
1105: !                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)1138:                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)
1106: !             ENDIF1139:             ENDIF
1107: !          ENDDO1140:          ENDDO
1108: ! !     ENDIF1141: !     ENDIF
1109: !    ENDDO1142:    ENDDO
1110: ! 1143: 
1111: ! ENDIF1144: ENDIF
1112: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest spring  contribution for any image in image ',IMAX1145: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest spring  contribution for any image in image ',IMAX
1113: IF (DEBUG) WRITE(*, '(A,3G20.10)') 'congrad2> ECON,EREP,ESPRING=',ECON,EREP,ESPRING1146: IF (DEBUG) WRITE(*, '(A,3G20.10)') 'congrad2> ECON,EREP,ESPRING=',ECON,EREP,ESPRING
1114: !1147: !
1115: ! Set gradients on frozen atoms to zero.1148: ! Set gradients on frozen atoms to zero.
1116: !1149: !
1117: IF (FREEZE) THEN1150: IF (FREEZE) THEN
1118:    DO J1=2,INTIMAGE+1  1151:    DO J1=2,INTIMAGE+1  
1119:       DO J2=1,NATOMS1152:       DO J2=1,NATOMS
1120:          IF (FROZEN(J2)) THEN1153:          IF (FROZEN(J2)) THEN
1121:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D01154:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0
1176:    ENDIF1209:    ENDIF
1177:    IF (CONEINT(J1)+REPEINT(J1).GT.MAXINT) THEN1210:    IF (CONEINT(J1)+REPEINT(J1).GT.MAXINT) THEN
1178:       MAXINT=CONEINT(J1)+REPEINT(J1)1211:       MAXINT=CONEINT(J1)+REPEINT(J1)
1179:       NMAXINT=J11212:       NMAXINT=J1
1180:    ENDIF1213:    ENDIF
1181: ENDDO1214: ENDDO
1182: ! IF (DEBUG) PRINT '(A,G20.10,A,2I6)',' congrad2> largest  internal energy=',MAXINT,' for image ',NMAXINT1215: ! IF (DEBUG) PRINT '(A,G20.10,A,2I6)',' congrad2> largest  internal energy=',MAXINT,' for image ',NMAXINT
1183: ! IF (DEBUG) PRINT '(A,G20.10,A,2I6)',' congrad2> smallest internal energy=',MININT,' for image ',NMININT1216: ! IF (DEBUG) PRINT '(A,G20.10,A,2I6)',' congrad2> smallest internal energy=',MININT,' for image ',NMININT
1184: IF (INTIMAGE.EQ.0) ETOTAL=EEE(1)+EEE(2)1217: IF (INTIMAGE.EQ.0) ETOTAL=EEE(1)+EEE(2)
1185: 1218: 
1186: ! IF ((RMS.LT.1.0D-50).AND.(ETOTAL.GT.0.1D0).AND.(INTIMAGE.GT.0.0D0)) THEN1219: IF ((RMS.LT.1.0D-50).AND.(ETOTAL.GT.0.1D0).AND.(INTIMAGE.GT.0.0D0)) THEN
1187: !    WRITE(*, '(A,2G20.10)') 'ETOTAL,RMS=',ETOTAL,RMS1220:    WRITE(*, '(A,2G20.10)') 'ETOTAL,RMS=',ETOTAL,RMS
1188: !    WRITE(*, '(A,G20.10)') 'ECON=',ECON1221:    WRITE(*, '(A,G20.10)') 'ECON=',ECON
1189: !    WRITE(*, '(A,G20.10)') 'EREP=',EREP1222:    WRITE(*, '(A,G20.10)') 'EREP=',EREP
1190: !    IF (PRINTE) STOP1223:    IF (PRINTE) STOP
1191: !    PRINTE=.TRUE.1224:    PRINTE=.TRUE.
1192: !    GOTO 1111225:    GOTO 111
1193: ! ENDIF1226: ENDIF
1194: 1227: 
1195: END SUBROUTINE CONGRAD21228: END SUBROUTINE CONGRAD2
1196: 1229: 
1197: SUBROUTINE INTMINONLY(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DINT,NOINT)1230: SUBROUTINE INTMINONLY(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DINT,NOINT)
1198: IMPLICIT NONE1231: IMPLICIT NONE
1199: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DINT,DUMMY1232: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DINT,DUMMY
1200: DOUBLE PRECISION DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq, r1amr1bdr2amr2b, r1amr1bdr2amr2bsq1233: DOUBLE PRECISION DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq, r1amr1bdr2amr2b, r1amr1bdr2amr2bsq
1201: LOGICAL NOINT1234: LOGICAL NOINT
1202: !1235: !
1203: ! Is there an internal extremum?1236: ! Is there an internal extremum?
1246:            1279:            
1247: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,J3,J4,J5,NINTMIN,NINTMIN2,MYUNIT1280: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,J3,J4,J5,NINTMIN,NINTMIN2,MYUNIT
1248: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS1281: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
1249: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D11282: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
1250: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)1283: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
1251: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI1284: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
1252: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)1285: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
1253: LOGICAL NOINT1286: LOGICAL NOINT
1254: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL1287: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL
1255: LOGICAL IMGFREEZE(INTIMAGE)1288: LOGICAL IMGFREEZE(INTIMAGE)
1256: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), X1, X2, Y1, Y2, Z1, Z2, DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq1289: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), X1, X2, Y1, Y2, Z1, Z2
1257: 1290: 
1258: EEE(1:INTIMAGE+2)=0.0D01291: EEE(1:INTIMAGE+2)=0.0D0
1259: CONE(1:INTIMAGE+2)=0.0D01292: CONE(1:INTIMAGE+2)=0.0D0
1260: REPE(1:INTIMAGE+2)=0.0D01293: REPE(1:INTIMAGE+2)=0.0D0
1261: REPEINT(1:INTIMAGE+2)=0.0D01294: REPEINT(1:INTIMAGE+2)=0.0D0
1262: NREPINT(1:INTIMAGE+2)=01295: NREPINT(1:INTIMAGE+2)=0
1263: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D01296: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0
1264: ECON=0.0D0; EREP=0.0D01297: ECON=0.0D0; EREP=0.0D0
1265: MYUNIT=61298: MYUNIT=6
1266: 1299: 
1393: ! Edge J1 is between images J1-1 and J1, starting from J1=2.1426: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
1394: ! Energy contributions are shared evenly, except for1427: ! Energy contributions are shared evenly, except for
1395: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which1428: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which
1396: ! is assigned to image INTIMAGE+1. Gradients are set to zero for1429: ! is assigned to image INTIMAGE+1. Gradients are set to zero for
1397: ! the end images.1430: ! the end images.
1398: !1431: !
1399:       IF (J1.EQ.1) CYCLE1432:       IF (J1.EQ.1) CYCLE
1400:       IF (QCINOREPINT) CYCLE1433:       IF (QCINOREPINT) CYCLE
1401:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)1434:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
1402:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)1435:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
1403: 1436:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
1404:       G1(1:3)=XYZ(NI1+1:NI1+3)-XYZ(NJ1+1:NJ1+3)1437:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
1405:       G2(1:3)=XYZ(NI2+1:NI2+3)-XYZ(NJ2+1:NJ2+3) 
1406: ! 
1407: ! Squared distance between atoms A and B for theta=0 - distance in image 2 
1408: ! 
1409:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2 
1410: ! 
1411: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1 
1412: ! 
1413:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2 
1414:       DCUT=NREPCUT(J2)**2 
1415:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff 
1416:       r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3) 
1417: ! 
1418: ! Is there an internal extremum? 
1419: ! 
1420:       r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b 
1421:  
1422: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN1438: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
1423:       IF (r1apr2bmr2amr1bsq.LT.1.0D-10) THEN1439:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-10) THEN
1424: !        D1=1.0D100; D2=1.0D100;1440: !        WRITE(*, '(A,I6,A,2I6)') 'B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2)
 1441: !        WRITE(*, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ
 1442: !        WRITE(*, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ
 1443: !        WRITE(*,'(A,7I10)') 'congrad3> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2)
 1444: !        WRITE(*,'(A,7I10)') 'frames ',J1-1,J1
1425:          NOINT=.TRUE.1445:          NOINT=.TRUE.
1426:          D1=SQRT(DSQ1) 
1427:          D2=SQRT(DSQ2) 
1428:          G2(1:3)=G2(1:3)/D2 
1429:          G1(1:3)=G1(1:3)/D1 
1430:       ELSE1446:       ELSE
1431:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq) 1447:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
 1448:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
1432:          IF (.NOT.NOINT) THEN1449:          IF (.NOT.NOINT) THEN
1433: !           WRITE(*,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad3> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), &1450: !           WRITE(*,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad3> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), &
1434: ! &                        ' distance,cutoff=',DINT,NREPCUT(J2)1451: ! &                        ' distance,cutoff=',DINT,NREPCUT(J2)
1435:             NINTMIN=NINTMIN+11452:             NINTMIN=NINTMIN+1
1436:          ENDIF1453:          ENDIF
1437:       ENDIF1454:       ENDIF
1438:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN1455:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
1439:          NINTMIN2=NINTMIN2+11456:          NINTMIN2=NINTMIN2+1
1440: !        D12=DSQI**61457: !        D12=DSQI**6
1441:          D12=DSQI1458:          D12=DSQI


r33467/intlbfgs.f90 2017-11-13 12:30:18.675374667 +0000 r33466/intlbfgs.f90 2017-11-13 12:30:19.347383627 +0000
 65: ! efk: for freezenodes 65: ! efk: for freezenodes
 66: ! 66: !
 67: DOUBLE PRECISION :: TESTG, TOTGNORM 67: DOUBLE PRECISION :: TESTG, TOTGNORM
 68: INTEGER :: IM 68: INTEGER :: IM
 69: ! 69: !
 70: ! Dimensions involving INTIMAGE 70: ! Dimensions involving INTIMAGE
 71: ! 71: !
 72: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), & 72: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), &
 73:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), & 73:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), &
 74:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:) 74:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:)
 75: DOUBLE PRECISION, ALLOCATABLE :: VPLUS(:), VMINUS(:)    75: ! DOUBLE PRECISION, ALLOCATABLE :: VPLUS(:), VMINUS(:)   
 76: DOUBLE PRECISION  EPLUS, EMINUS, DIFF    76: ! DOUBLE PRECISION  EPLUS, EMINUS, DIFF   
 77: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:) 77: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:)
 78: ! saved interpolation 78: ! saved interpolation
 79: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE 79: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE
 80: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:) 80: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:)
 81: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP), CHIRALACTIVE(NATOMS) 81: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP), CHIRALACTIVE(NATOMS)
 82: INTEGER LUNIT, GETUNIT 82: INTEGER LUNIT, GETUNIT
 83: CHARACTER(LEN=2) SDUMMY 83: CHARACTER(LEN=2) SDUMMY
 84: INTEGER JMAXEEE,JMAXRMS 84: INTEGER JMAXEEE,JMAXRMS
 85: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT 85: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT
 86:  86: 
129: ENDIF129: ENDIF
130: 130: 
131: ALLOCATE(TRUEEE(INTIMAGE+2), &131: ALLOCATE(TRUEEE(INTIMAGE+2), &
132:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &132:   &      EEETMP(INTIMAGE+2), MYGTMP(3*NATOMS*INTIMAGE), &
133:   &      GTMP(3*NATOMS*INTIMAGE), &133:   &      GTMP(3*NATOMS*INTIMAGE), &
134:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:INTMUPDATE,(3*NATOMS)*INTIMAGE), &134:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:INTMUPDATE,(3*NATOMS)*INTIMAGE), &
135:   &      GDIF(0:INTMUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &135:   &      GDIF(0:INTMUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &
136:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &136:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &
137:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))137:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))
138: 138: 
139: ALLOCATE(VPLUS((3*NATOMS)*(INTIMAGE+2)),VMINUS((3*NATOMS)*(INTIMAGE+2)))  139: ! ALLOCATE(VPLUS((3*NATOMS)*(INTIMAGE+2)),VMINUS((3*NATOMS)*(INTIMAGE+2)))  
140: 140: 
141: SWITCHED=.FALSE.141: SWITCHED=.FALSE.
142: INTIMAGESAVE=INTIMAGE142: INTIMAGESAVE=INTIMAGE
143: NBACKTRACK=1143: NBACKTRACK=1
144: CALL MYCPU_TIME(STIME,.FALSE.)144: CALL MYCPU_TIME(STIME,.FALSE.)
145: WRITE(*,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1145: WRITE(*,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1
146: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',INTMUPDATE,' maximum step size=',MAXINTBFGS146: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',INTMUPDATE,' maximum step size=',MAXINTBFGS
147: ADDATOM=.FALSE.147: ADDATOM=.FALSE.
148: NFAIL=0148: NFAIL=0
149: IMGFREEZE(1:INTIMAGE)=.FALSE.149: IMGFREEZE(1:INTIMAGE)=.FALSE.
573: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.573: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.
574: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!574: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!
575: !575: !
576: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.576: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.
577: !577: !
578: BESTWORST=1.0D100578: BESTWORST=1.0D100
579: 9876 CONTINUE579: 9876 CONTINUE
580: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)580: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
581: 581: 
582: 582: 
583: !               DIFF=1.0D-4583: !                DIFF=1.0D-7
584: !               PRINT*,'analytic and numerical gradients:'584: !                PRINT*,'analytic and numerical gradients:'
585: !               CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)585: !                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
586: !               DO J1=2,INTIMAGE+1586: !                DO J1=2,INTIMAGE+1
587: !                  DO J2=1,3*NATOMS587: !                   DO J2=1,3*NATOMS
588: !                     J3=3*NATOMS*(J1-1)+J2588: !                      J3=3*NATOMS*(J1-1)+J2
589: !                     XYZ(J3)=XYZ(J3)+DIFF589: !                      XYZ(J3)=XYZ(J3)+DIFF
590: !                     CALL CONGRAD2(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)590: !                      CALL CONGRAD2(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)
591: !                     XYZ(J3)=XYZ(J3)-2.0D0*DIFF591: !                      XYZ(J3)=XYZ(J3)-2.0D0*DIFF
592: !                     CALL CONGRAD2(NMAXINT,NMININT,EMINUS,XYZ,VMINUS,EEE,IMGFREEZE,RMS)592: !                      CALL CONGRAD2(NMAXINT,NMININT,EMINUS,XYZ,VMINUS,EEE,IMGFREEZE,RMS)
593: !                     XYZ(J3)=XYZ(J3)+DIFF593: !                      XYZ(J3)=XYZ(J3)+DIFF
594: !                     IF ((ABS(GGG(J3)).NE.0.0D0).AND. &594: !                      IF ((ABS(GGG(J3)).NE.0.0D0).AND. &
595: !    &               (ABS(100.0D0*(GGG(J3)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GGG(J3)).GT.1.0D0)) THEN595: !     &               (ABS(100.0D0*(GGG(J3)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GGG(J3)).GT.1.0D0)) THEN
596: !                        WRITE(*,'(A,2I5,2F20.10,L5,A)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &596: !                         WRITE(*,'(A,2I5,2F20.10,L5,A)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &
597: !    & ATOMACTIVE((J2-1)/3+1),'   X'   597: !     & ATOMACTIVE((J2-1)/3+1),'   X'   
598: !                     ELSE598: !                      ELSE
599: !                        WRITE(*,'(A,2I5,2F20.10,L5)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &599: !                         WRITE(*,'(A,2I5,2F20.10,L5)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &
600: !    & ATOMACTIVE((J2-1)/3+1)  600: !     & ATOMACTIVE((J2-1)/3+1)  
601: !                     ENDIF601: !                      ENDIF
602: !                  ENDDO602: !                   ENDDO
603: !               ENDDO603: !                ENDDO
604: ! !604: ! 
605: !               STOP605: !                STOP
606: 606: 
607: 607: 
608: IF (QCIADDREP.GT.0) THEN608: IF (QCIADDREP.GT.0) THEN
609:    CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)609:    CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
610: !610: !
611: ! Don't do the CONINT part of CONGRAD2 if CONINT isn't set. CONGRAD seems to be611: ! Don't do the CONINT part of CONGRAD2 if CONINT isn't set. CONGRAD seems to be
612: ! dong something different at the moment. Focus on CONGRAD2612: ! dong something different at the moment. Focus on CONGRAD2
613: !613: !
614: ELSEIF (CHECKCONINT) THEN614: ! ELSEIF (CHECKCONINT) THEN
615:    CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
616: ELSE615: ELSE
617:    CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)616:    CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 617: ! ELSE
 618:    ! CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
618: ENDIF619: ENDIF
619: EOLD=ETOTAL620: EOLD=ETOTAL
620: GLAST(1:D)=G(1:D)621: GLAST(1:D)=G(1:D)
621: XSAVE(1:D)=X(1:D)622: XSAVE(1:D)=X(1:D)
622: 623: 
623: IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN624: IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN
624:    WRITE(*,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=', &625:    WRITE(*,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=', &
625:   &                       ETOTAL/INTIMAGE,COLDFUSIONLIMIT626:   &                       ETOTAL/INTIMAGE,COLDFUSIONLIMIT
626:    DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST,CONOFFTRIED)627:    DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST,CONOFFTRIED)
627:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &628:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
754:    DO J1=1,NPERMGROUP755:    DO J1=1,NPERMGROUP
755:       IF (GROUPACTIVE(J1)) GOTO 975756:       IF (GROUPACTIVE(J1)) GOTO 975
756:       DO J2=1,NPERMSIZE(J1)757:       DO J2=1,NPERMSIZE(J1)
757:          IF (.NOT.ATOMACTIVE(PERMGROUP(NDUMMY+J2-1))) GOTO 975758:          IF (.NOT.ATOMACTIVE(PERMGROUP(NDUMMY+J2-1))) GOTO 975
758:       ENDDO759:       ENDDO
759:       GROUPACTIVE(J1)=.TRUE.760:       GROUPACTIVE(J1)=.TRUE.
760:       IF (DEBUG) WRITE(*,'(A,I6,A)') ' intlbfgs> All permutable atoms in group ',J1,' are active'761:       IF (DEBUG) WRITE(*,'(A,I6,A)') ' intlbfgs> All permutable atoms in group ',J1,' are active'
761: 975   NDUMMY=NDUMMY+NPERMSIZE(J1)762: 975   NDUMMY=NDUMMY+NPERMSIZE(J1)
762:    ENDDO 763:    ENDDO 
763: 764: 
764: !  IF (NITERDONE.EQ.1) THEN765:    IF (NITERDONE.EQ.1) THEN
765: !     COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint766:       COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint
766: !     COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint767:       COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint
767: !     WRITE(*,'(A)') ' intlbfgs> checking alignment of endpoints - should be no permutations'768:       WRITE(*,'(A)') ' intlbfgs> checking alignment of endpoints - should be no permutations'
768: !     CALL LOPERMDIST(COORDSB,COORDSC,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,.FALSE.,RMATBEST,0,NMOVE,PERM)769:       CALL LOPERMDIST(COORDSB,COORDSC,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,.FALSE.,RMATBEST,0,NMOVE,PERM)
769: !     WRITE(*,'(A,G20.10,A,I6)') ' intlbfgs> endpoint distance ',DISTANCE,' permutations=',NMOVE 770:       WRITE(*,'(A,G20.10,A,I6)') ' intlbfgs> endpoint distance ',DISTANCE,' permutations=',NMOVE 
770: !  ENDIF771:    ENDIF
771: 772: 
772:    SAVELOCALPERMCUT=LOCALPERMCUT773:    SAVELOCALPERMCUT=LOCALPERMCUT
773:    LOCALPERMCUT=QCIPERMCUT ! 0.8D0 774:    LOCALPERMCUT=QCIPERMCUT ! 0.8D0 
774:    np: DO J1=1,NPERMGROUP775:    np: DO J1=1,NPERMGROUP
775:       IF (.NOT.GROUPACTIVE(J1)) CYCLE np776:       IF (.NOT.GROUPACTIVE(J1)) CYCLE np
776:       DO J3=1,INTIMAGE777:       DO J3=1,INTIMAGE
777: !        WRITE(*,'(A,I6,A,I6)') 'intlbfgs> Doing group ',J1,' image ',J3+1778: !        WRITE(*,'(A,I6,A,I6)') 'intlbfgs> Doing group ',J1,' image ',J3+1
 779: 
778: !          COORDSB(1:3*NATOMS)=XYZ(3*NATOMS*(J3-1)+1:3*NATOMS*J3) ! coordinates for intervening image J3-1, which is the starting endpoint for J3=1780: !          COORDSB(1:3*NATOMS)=XYZ(3*NATOMS*(J3-1)+1:3*NATOMS*J3) ! coordinates for intervening image J3-1, which is the starting endpoint for J3=1
779: !          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))  ! coordinates for intervening image J3, which is image J3+1 including starting endpoint781: !          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))  ! coordinates for intervening image J3, which is image J3+1 including starting endpoint
780: !          CALL QCIOPTPERM(COORDSB,COORDSA,NATOMS,DEBUG,J1,PERM,STARTGROUP,J3)782: !          CALL QCIOPTPERM(COORDSB,COORDSA,NATOMS,DEBUG,J1,PERM,STARTGROUP,J3)
781: !          IF (.NOT.IDENTITY) THEN783: !          IF (.NOT.IDENTITY) THEN
782: !             DO J2=1,NATOMS784: !             DO J2=1,NATOMS
783: ! !           IF (PERM(J2).NE.J2) WRITE(*,'(4(A,I6))') ' intlbfgs> image ',J3+1,' qcipermopt would move atom ',J2,' to position ',PERM(J2), &  785: ! !           IF (PERM(J2).NE.J2) WRITE(*,'(4(A,I6))') ' intlbfgs> image ',J3+1,' qcipermopt would move atom ',J2,' to position ',PERM(J2), &  
784: ! !  &   ' group ',J1  786: ! !  &   ' group ',J1  
785: !                IF (PERM(J2).NE.J2) WRITE(*,'(4(A,I6))') ' intlbfgs> image ',J3+1,' qcipermopt move atom ',J2,' to position ',PERM(J2), &787: !                IF (PERM(J2).NE.J2) WRITE(*,'(4(A,I6))') ' intlbfgs> image ',J3+1,' qcipermopt move atom ',J2,' to position ',PERM(J2), &
786: !    &   ' group ',J1  788: !    &   ' group ',J1  
787: !                COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERM(J2)-1)+1)789: !                COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERM(J2)-1)+1)
1161:             EEE(1:INTIMAGE+3)=0.0D01163:             EEE(1:INTIMAGE+3)=0.0D0
1162:             STEPIMAGE(1:INTIMAGE+1)=0.0D01164:             STEPIMAGE(1:INTIMAGE+1)=0.0D0
1163: 1165: 
1164:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))1166:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))
1165:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))1167:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))
1166:             INTIMAGE=INTIMAGE+11168:             INTIMAGE=INTIMAGE+1
1167:             D=(3*NATOMS)*INTIMAGE1169:             D=(3*NATOMS)*INTIMAGE
1168:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1170:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1169:             IF (QCIADDREP.GT.0) THEN1171:             IF (QCIADDREP.GT.0) THEN
1170:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1172:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1171:             ELSEIF (CHECKCONINT) THEN1173: !           ELSEIF (CHECKCONINT) THEN
1172:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1173:             ELSE1174:             ELSE
1174:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1175:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 1176: !           ELSE
 1177: !              CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1175:             ENDIF1178:             ENDIF
1176: !           GOTO 8641179: !           GOTO 864
1177:          ENDIF1180:          ENDIF
1178:          IF (REMOVEIMAGE.OR.((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1))) THEN1181:          IF (REMOVEIMAGE.OR.((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1))) THEN
1179:             IF (REMOVEIMAGE) JMIN=JMAXEEE1182:             IF (REMOVEIMAGE) JMIN=JMAXEEE
1180:             IF (JMIN.EQ.1) JMIN=21183:             IF (JMIN.EQ.1) JMIN=2
1181:             WRITE(*,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN1184:             WRITE(*,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN
1182:             NITERUSE=01185:             NITERUSE=0
1183:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))1186:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
1184:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))1187:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
1241:             EEE(1:INTIMAGE+1)=0.0D01244:             EEE(1:INTIMAGE+1)=0.0D0
1242:             STEPIMAGE(1:INTIMAGE-1)=0.0D01245:             STEPIMAGE(1:INTIMAGE-1)=0.0D0
1243: 1246: 
1244:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))1247:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))
1245:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))1248:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))
1246:             INTIMAGE=INTIMAGE-11249:             INTIMAGE=INTIMAGE-1
1247:             D=(3*NATOMS)*INTIMAGE1250:             D=(3*NATOMS)*INTIMAGE
1248:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1251:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1249:             IF (QCIADDREP.GT.0) THEN1252:             IF (QCIADDREP.GT.0) THEN
1250:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1253:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1251:             ELSEIF (CHECKCONINT) THEN1254: !           ELSEIF (CHECKCONINT) THEN
1252:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1253:             ELSE1255:             ELSE
1254:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1256:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 1257: !           ELSE
 1258: !              CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1255:             ENDIF1259:             ENDIF
1256:             NLASTGOODE=NITERDONE1260:             NLASTGOODE=NITERDONE
1257:             LASTGOODE=ETOTAL1261:             LASTGOODE=ETOTAL
1258: !           GOTO 8641262: !           GOTO 864
1259:          ENDIF1263:          ENDIF
1260:       ELSE1264:       ELSE
1261:          DMAX=-1.0D01265:          DMAX=-1.0D0
1262:          ADMAX=-1.0D01266:          ADMAX=-1.0D0
1263:          DMIN=HUGE(1.0D0)1267:          DMIN=HUGE(1.0D0)
1264:          DUMMY2=0.0D01268:          DUMMY2=0.0D0
1321: !     DO J2=2,INTIMAGE+21325: !     DO J2=2,INTIMAGE+2
1322: !        CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &1326: !        CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
1323: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)1327: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
1324: !     ENDDO1328: !     ENDDO
1325: !  ENDIF1329: !  ENDIF
1326: 1330: 
1327:    IF (.NOT.SWITCHED) THEN1331:    IF (.NOT.SWITCHED) THEN
1328:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1332:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1329:       IF (QCIADDREP.GT.0) THEN1333:       IF (QCIADDREP.GT.0) THEN
1330:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1334:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1331:       ELSEIF (CHECKCONINT) THEN1335: !     ELSEIF (CHECKCONINT) THEN
1332:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1333:       ELSE1336:       ELSE
1334:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1337:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 1338: !     ELSE
 1339: !        CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1335:       ENDIF1340:       ENDIF
1336: 1341: 
1337:       IF ((ETOTAL-EOLD.LT.1.0D100).OR.ADDATOM) THEN ! MAXERISE effectively set to 1.0D100 here1342:       IF ((ETOTAL-EOLD.LT.1.0D100).OR.ADDATOM) THEN ! MAXERISE effectively set to 1.0D100 here
1338:          EOLD=ETOTAL1343:          EOLD=ETOTAL
1339:          GLAST(1:D)=G(1:D)1344:          GLAST(1:D)=G(1:D)
1340:          XSAVE(1:D)=X(1:D)1345:          XSAVE(1:D)=X(1:D)
1341:       ELSE1346:       ELSE
1342:          NDECREASE=NDECREASE+11347:          NDECREASE=NDECREASE+1
1343:          IF (NDECREASE.GT.5) THEN1348:          IF (NDECREASE.GT.5) THEN
1344:             NFAIL=NFAIL+11349:             NFAIL=NFAIL+1
1364:             CALL POTENTIAL(XYZ((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),GGG((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),EEE(J4), &1369:             CALL POTENTIAL(XYZ((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),GGG((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),EEE(J4), &
1365:   &                                    .TRUE.,.FALSE.)1370:   &                                    .TRUE.,.FALSE.)
1366:             ETOTALTMP=ETOTALTMP+EEE(J4)1371:             ETOTALTMP=ETOTALTMP+EEE(J4)
1367:          ENDDO1372:          ENDDO
1368:       ENDIF1373:       ENDIF
1369:       EEETMP(1:INTIMAGE+2)=EEE(1:INTIMAGE+2)1374:       EEETMP(1:INTIMAGE+2)=EEE(1:INTIMAGE+2)
1370:       MYGTMP(1:D)=G(1:D)1375:       MYGTMP(1:D)=G(1:D)
1371:       IF (USEFRAC.LT.1.0D0) THEN1376:       IF (USEFRAC.LT.1.0D0) THEN
1372:          IF (QCIADDREP.GT.0) THEN1377:          IF (QCIADDREP.GT.0) THEN
1373:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1378:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1374:          ELSEIF (CHECKCONINT) THEN1379: !        ELSEIF (CHECKCONINT) THEN
1375:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1376:          ELSE1380:          ELSE
1377:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1381:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 1382: !        ELSE
 1383: !           CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1378:          ENDIF1384:          ENDIF
1379:       ELSE1385:       ELSE
1380:          ETOTAL=0.0D01386:          ETOTAL=0.0D0
1381:          G(1:D)=0.0D01387:          G(1:D)=0.0D0
1382:       ENDIF1388:       ENDIF
1383:       ETOTAL=USEFRAC*ETOTALTMP+(1.0D0-USEFRAC)*ETOTAL1389:       ETOTAL=USEFRAC*ETOTALTMP+(1.0D0-USEFRAC)*ETOTAL
1384:       G(1:D)=USEFRAC*MYGTMP(1:D)+(1.0D0-USEFRAC)*G(1:D)1390:       G(1:D)=USEFRAC*MYGTMP(1:D)+(1.0D0-USEFRAC)*G(1:D)
1385:       RMS=SUM(G(1:D)**2)1391:       RMS=SUM(G(1:D)**2)
1386:       RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))1392:       RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))
1387:       EEE(1:INTIMAGE+2)=USEFRAC*EEETMP(1:INTIMAGE+2)+(1.0D0-USEFRAC)*EEE(1:INTIMAGE+2)1393:       EEE(1:INTIMAGE+2)=USEFRAC*EEETMP(1:INTIMAGE+2)+(1.0D0-USEFRAC)*EEE(1:INTIMAGE+2)
1540:             DO J1=1,NATOMS1546:             DO J1=1,NATOMS
1541:                IF (ATOMACTIVE(J1)) WRITE(*,'(A,I6)') ' active atom ',J11547:                IF (ATOMACTIVE(J1)) WRITE(*,'(A,I6)') ' active atom ',J1
1542:             ENDDO1548:             ENDDO
1543:             STOP1549:             STOP
1544:          ENDIF1550:          ENDIF
1545:          ADDATOM=.TRUE.1551:          ADDATOM=.TRUE.
1546: 1552: 
1547:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1553:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1548:          IF (QCIADDREP.GT.0) THEN1554:          IF (QCIADDREP.GT.0) THEN
1549:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1555:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1550:          ELSEIF (CHECKCONINT) THEN1556: !        ELSEIF (CHECKCONINT) THEN
1551:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1552:          ELSE1557:          ELSE
1553:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1558:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 1559: !        ELSE
 1560: !           CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1554:          ENDIF1561:          ENDIF
1555:       ENDIF1562:       ENDIF
1556:       LASTGOODE=ETOTAL1563:       LASTGOODE=ETOTAL
1557:    ENDIF1564:    ENDIF
1558: 1565: 
1559:    EXITSTATUS=01566:    EXITSTATUS=0
1560:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1567:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1561:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 1568:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 
1562:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 1569:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 
1563:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21570:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1699: 1706: 
1700: NNREPULSIVE=NNSTART1707: NNREPULSIVE=NNSTART
1701: DO JJ=NSTART,NREPULSIVE1708: DO JJ=NSTART,NREPULSIVE
1702:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**21709:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**2
1703:    NI=REPI(JJ)1710:    NI=REPI(JJ)
1704:    NJ=REPJ(JJ)1711:    NJ=REPJ(JJ)
1705:    DO KK=1,INTIMAGE+2 ! first check for standard distances within threshold1712:    DO KK=1,INTIMAGE+2 ! first check for standard distances within threshold
1706:       LDIST=(XYZ((KK-1)*NOPT+3*(NI-1)+1)-XYZ((KK-1)*NOPT+3*(NJ-1)+1))**2 &1713:       LDIST=(XYZ((KK-1)*NOPT+3*(NI-1)+1)-XYZ((KK-1)*NOPT+3*(NJ-1)+1))**2 &
1707:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+2)-XYZ((KK-1)*NOPT+3*(NJ-1)+2))**2 &1714:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+2)-XYZ((KK-1)*NOPT+3*(NJ-1)+2))**2 &
1708:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+3)-XYZ((KK-1)*NOPT+3*(NJ-1)+3))**21715:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+3)-XYZ((KK-1)*NOPT+3*(NJ-1)+3))**2
 1716: !     IF ((REPI(JJ).EQ.2024).OR.(REPJ(JJ).EQ.2024)) THEN
 1717: !        WRITE(*,'(A,2I12,2I6,3G20.10)') 'checkrep> KK,JJ,repi,repj,ldist,compare,repcut=', &
 1718: ! &                                    KK,JJ,repi(JJ),repj(JJ),ldist,compare,repcut(jj)
 1719: !     ENDIF
1709:       IF (LDIST.LT.COMPARE) THEN1720:       IF (LDIST.LT.COMPARE) THEN
1710:          NNREPULSIVE=NNREPULSIVE+11721:          NNREPULSIVE=NNREPULSIVE+1
1711:          NREPI(NNREPULSIVE)=NI1722:          NREPI(NNREPULSIVE)=NI
1712:          NREPJ(NNREPULSIVE)=NJ1723:          NREPJ(NNREPULSIVE)=NJ
1713:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)1724:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)
1714: !        IF ((REPI(JJ).EQ.2024).OR.(REPJ(JJ).EQ.2024)) THEN1725: !        IF ((REPI(JJ).EQ.2024).OR.(REPJ(JJ).EQ.2024)) THEN
1715: !           WRITE(*,'(A)') 'checkrep> KK,JJ,ldist < compare : active'1726: !           WRITE(*,'(A)') 'checkrep> KK,JJ,ldist < compare : active'
1716: !        ENDIF1727: !        ENDIF
1717:          GOTO 2461728:          GOTO 246
1718:       ENDIF1729:       ENDIF
1719:    ENDDO 1730:    ENDDO 
1720: !1731:    COMPARE=CHECKREPCUTOFF*REPCUT(JJ)
1721: ! We don't check for internal minima in repulsions in congrad now unless both distances are1732:    DO KK=2,INTIMAGE+2 ! now check internal minima within threshold
1722: ! within threshold.1733:       DMIN=1.0D10
1723: !1734:       NI2=NOPT*(KK-2)+3*(NI-1)
1724: !  COMPARE=CHECKREPCUTOFF*REPCUT(JJ)1735:       NI1=NOPT*(KK-1)+3*(NI-1)
1725: !  DO KK=2,INTIMAGE+2 ! now check internal minima within threshold1736:       NJ2=NOPT*(KK-2)+3*(NJ-1)
1726: !     DMIN=1.0D101737:       NJ1=NOPT*(KK-1)+3*(NJ-1)
1727: !     NI2=NOPT*(KK-2)+3*(NI-1)1738:       R1AX=XYZ(NI2+1); R1AY=XYZ(NI2+2); R1AZ=XYZ(NI2+3)
1728: !     NI1=NOPT*(KK-1)+3*(NI-1)1739:       R1BX=XYZ(NJ2+1); R1BY=XYZ(NJ2+2); R1BZ=XYZ(NJ2+3)
1729: !     NJ2=NOPT*(KK-2)+3*(NJ-1)1740:       R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3)
1730: !     NJ1=NOPT*(KK-1)+3*(NJ-1)1741:       R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3)
1731: !     R1AX=XYZ(NI2+1); R1AY=XYZ(NI2+2); R1AZ=XYZ(NI2+3)1742:       CALL INTMINONLY(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DMIN,NOINT)
1732: !     R1BX=XYZ(NJ2+1); R1BY=XYZ(NJ2+2); R1BZ=XYZ(NJ2+3)1743: 
1733: !     R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3)1744:       IF (NOINT) CYCLE
1734: !     R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3)1745:       IF (DMIN.LT.COMPARE) THEN
1735: !     CALL INTMINONLY(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DMIN,NOINT)1746:          NNREPULSIVE=NNREPULSIVE+1
1736: 1747:          NREPI(NNREPULSIVE)=NI
1737: !     IF (NOINT) CYCLE1748:          NREPJ(NNREPULSIVE)=NJ
1738: !     IF (DMIN.LT.COMPARE) THEN1749:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)
1739: !        NNREPULSIVE=NNREPULSIVE+11750:          GOTO 246
1740: !        NREPI(NNREPULSIVE)=NI1751:       ENDIF
1741: !        NREPJ(NNREPULSIVE)=NJ1752:    ENDDO 
1742: !        NREPCUT(NNREPULSIVE)=REPCUT(JJ) 
1743: !        GOTO 246 
1744: !     ENDIF 
1745: !  ENDDO  
1746: 246 CONTINUE1753: 246 CONTINUE
1747: ENDDO1754: ENDDO
1748: IF (DEBUG) WRITE(*,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE1755: IF (DEBUG) WRITE(*,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE
1749: 1756: 
1750: END SUBROUTINE CHECKREP1757: END SUBROUTINE CHECKREP
1751: 1758: 
1752: SUBROUTINE INTRWG(NACTIVE,NITER,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1759: SUBROUTINE INTRWG(NACTIVE,NITER,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1753: USE PORFUNCS1760: USE PORFUNCS
1754: USE KEY,ONLY: STOCKT,STOCKAAT, RBAAT, ATOMACTIVE, NCONSTRAINT, CONACTIVE, NREPULSIVE, NNREPULSIVE, REPI, REPJ, REPCUT, NREPCUT, &1761: USE KEY,ONLY: STOCKT,STOCKAAT, RBAAT, ATOMACTIVE, NCONSTRAINT, CONACTIVE, NREPULSIVE, NNREPULSIVE, REPI, REPJ, REPCUT, NREPCUT, &
1755:   &           NREPMAX, NREPI, NREPJ, INTFROZEN, CONOFFLIST,CONOFFTRIED1762:   &           NREPMAX, NREPI, NREPJ, INTFROZEN, CONOFFLIST,CONOFFTRIED
2384:                      ELSE2391:                      ELSE
2385:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)2392:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)
2386:                      ENDIF2393:                      ENDIF
2387:                   ENDIF2394:                   ENDIF
2388:                ENDIF2395:                ENDIF
2389: 2396: 
2390:             ENDDO2397:             ENDDO
2391:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2398:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2392:             IF (QCIADDREP.GT.0) THEN2399:             IF (QCIADDREP.GT.0) THEN
2393:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2400:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2394:             ELSEIF (CHECKCONINT) THEN2401: !           ELSEIF (CHECKCONINT) THEN
2395:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2396:             ELSE2402:             ELSE
2397:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2403:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2404: !           ELSE
 2405: !              CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2398:             ENDIF2406:             ENDIF
2399:             ESAVE0=ETOTAL2407:             ESAVE0=ETOTAL
2400:             DO J1=2,INTIMAGE+12408:             DO J1=2,INTIMAGE+1
2401:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2409:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2402:             ENDDO2410:             ENDDO
2403:          ENDIF2411:          ENDIF
2404: 2412: 
2405: !          IF ((.NOT.IDONE).AND.(NDFORNEWATOM.GE.3)) THEN2413: !          IF ((.NOT.IDONE).AND.(NDFORNEWATOM.GE.3)) THEN
2406: !             IDONE=.TRUE.2414: !             IDONE=.TRUE.
2407: ! !2415: ! !
2630:                      ELSE2638:                      ELSE
2631:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)2639:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)
2632:                      ENDIF2640:                      ENDIF
2633:                   ENDIF2641:                   ENDIF
2634:                ENDIF2642:                ENDIF
2635:             ENDDO2643:             ENDDO
2636: 2644: 
2637:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2645:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2638:             IF (QCIADDREP.GT.0) THEN2646:             IF (QCIADDREP.GT.0) THEN
2639:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2647:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2640:             ELSEIF (CHECKCONINT) THEN2648: !           ELSEIF (CHECKCONINT) THEN
2641:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2642:             ELSE2649:             ELSE
2643:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2650:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2651: !           ELSE
 2652: !              CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2644:             ENDIF2653:             ENDIF
2645:             ESAVEC=ETOTAL2654:             ESAVEC=ETOTAL
2646:             DO J1=2,INTIMAGE+12655:             DO J1=2,INTIMAGE+1
2647:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2656:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2648:             ENDDO2657:             ENDDO
2649:          ENDIF2658:          ENDIF
2650: !2659: !
2651: ! Standard linear interpolation, with constraint distance scaled by FRAC.2660: ! Standard linear interpolation, with constraint distance scaled by FRAC.
2652: ! Works for FRAC as small as 0.1 with repulsion turned off.2661: ! Works for FRAC as small as 0.1 with repulsion turned off.
2653: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements2662: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements
2663:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &2672:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &
2664:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &2673:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &
2665:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)2674:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)
2666:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &2675:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &
2667:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &2676:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &
2668:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)2677:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)
2669:             ENDDO2678:             ENDDO
2670:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2679:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2671:             IF (QCIADDREP.GT.0) THEN2680:             IF (QCIADDREP.GT.0) THEN
2672:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2681:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2673:             ELSEIF (CHECKCONINT) THEN2682: !           ELSEIF (CHECKCONINT) THEN
2674:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2675:             ELSE2683:             ELSE
2676:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2684:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2685: !           ELSE
 2686: !              CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2677:             ENDIF2687:             ENDIF
2678:          ENDIF2688:          ENDIF
2679: 2689: 
2680:          IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &2690:          IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &
2681:   &              ESAVE0,ESAVED,ESAVEC,ETOTAL2691:   &              ESAVE0,ESAVED,ESAVEC,ETOTAL
2682:     2692:     
2683:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN2693:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN
2684:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'2694:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'
2685:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN2695:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN
2686:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'2696:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'
2750: ! Turn frozen images off for new added atom.2760: ! Turn frozen images off for new added atom.
2751: !2761: !
2752: !     IF (DEBUG) WRITE(*,'(A)') ' intlbfgs> turning off frozen images'2762: !     IF (DEBUG) WRITE(*,'(A)') ' intlbfgs> turning off frozen images'
2753: !     IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.2763: !     IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.
2754:       CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2764:       CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2755: !2765: !
2756: ! need a new gradient since the active atom has changed !2766: ! need a new gradient since the active atom has changed !
2757: !2767: !
2758:       IF (QCIADDREP.GT.0) THEN2768:       IF (QCIADDREP.GT.0) THEN
2759:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2769:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2760:       ELSEIF (CHECKCONINT) THEN2770: !     ELSEIF (CHECKCONINT) THEN
2761:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2762:       ELSE2771:       ELSE
2763:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2772:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2773: !     ELSE
 2774: !        CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2764:       ENDIF2775:       ENDIF
2765: 2776: 
2766: END SUBROUTINE DOADDATOM2777: END SUBROUTINE DOADDATOM
2767: 2778: 
2768: SUBROUTINE CHECKPERC(LXYZ,LINTCONSTRAINTTOL,NQCIFREEZE,NCPFIT)2779: SUBROUTINE CHECKPERC(LXYZ,LINTCONSTRAINTTOL,NQCIFREEZE,NCPFIT)
2769: USE KEY, ONLY : ATOMACTIVE, NCONSTRAINT, INTFROZEN, CONI, CONJ, CONDISTREF, INTCONMAX, INTCONSTRAINTTOL, &2780: USE KEY, ONLY : ATOMACTIVE, NCONSTRAINT, INTFROZEN, CONI, CONJ, CONDISTREF, INTCONMAX, INTCONSTRAINTTOL, &
2770:   &             INTCONSEP, NCONGEOM, CONGEOM, CONIFIX, CONJFIX, CONDISTREFFIX, INTCONCUT, &2781:   &             INTCONSEP, NCONGEOM, CONGEOM, CONIFIX, CONJFIX, CONDISTREFFIX, INTCONCUT, &
2771:   &             NCONSTRAINTFIX, BULKT, TWOD, RIGIDBODY, CONDATT, CONCUT, CONCUTFIX, &2782:   &             NCONSTRAINTFIX, BULKT, TWOD, RIGIDBODY, CONDATT, CONCUT, CONCUTFIX, &
2772:   &             BONDS, QCIAMBERT, QCIADDREP, QCIADDREPCUT, QCIBONDS, QCISECOND2783:   &             BONDS, QCIAMBERT, QCIADDREP, QCIADDREPCUT, QCIBONDS, QCISECOND
2773: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM32784: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3


r33467/myorient.f 2017-11-13 12:30:18.895377599 +0000 r33466/myorient.f 2017-11-13 12:30:19.575386667 +0000
302: C        ENDDO302: C        ENDDO
303: C     ENDDO303: C     ENDDO
304: C     PRINT '(A)','myorient> ROTINV * T1:'304: C     PRINT '(A)','myorient> ROTINV * T1:'
305: C     PRINT '(6F20.10)',TVEC(1:12)305: C     PRINT '(6F20.10)',TVEC(1:12)
306: 306: 
307:       RETURN307:       RETURN
308:       END308:       END
309: 309: 
310:       SUBROUTINE ROTXZ(NATOMS,JDO,T1,PROJ,DIST,XVEC,YVEC,ZVEC,CUTOFF1)310:       SUBROUTINE ROTXZ(NATOMS,JDO,T1,PROJ,DIST,XVEC,YVEC,ZVEC,CUTOFF1)
311:       IMPLICIT NONE311:       IMPLICIT NONE
312:       INTEGER NATOMS, JDO, J1, J3, J2312:       INTEGER NATOMS, JDO, J1, J3
313:       DOUBLE PRECISION T1(3*NATOMS), PROJ, DIST(NATOMS), RVEC(3), COST, SINT, RDOTN, TX, TY, TZ, CUTOFF1313:       DOUBLE PRECISION T1(3*NATOMS), PROJ, DIST(NATOMS), RVEC(3), COST, SINT, RDOTN, TX, TY, TZ, CUTOFF1
314:       DOUBLE PRECISION XVEC(3), YVEC(3), ZVEC(3)314:       DOUBLE PRECISION XVEC(3), YVEC(3), ZVEC(3)
315: 315: 
316: C     PRINT '(I6)',NATOMS316: C     PRINT '(I6)',NATOMS
317: C     PRINT '(A,I6,G20.10)','before rotation'317: C     PRINT '(A,I6,G20.10)','before rotation'
318: C     WRITE(*,'(A2,2X,3G20.10)') ('LA',T1(3*(J3-1)+1),T1(3*(J3-1)+2),T1(3*(J3-1)+3),J3=1,NATOMS)318: C     WRITE(*,'(A2,2X,3G20.10)') ('LA',T1(3*(J3-1)+1),T1(3*(J3-1)+2),T1(3*(J3-1)+3),J3=1,NATOMS)
319: 319: 
320:       IF (ABS(T1(3*(JDO-1)+2)).LT.1.0D-8) THEN320:       IF (ABS(T1(3*(JDO-1)+2)).LT.1.0D-8) THEN
321:          IF (T1(3*(JDO-1)+1).LT.0.0D0) THEN ! rotate about the z axis DO NOT INVERT!!321:          IF (T1(3*(JDO-1)+1).LT.0.0D0) THEN ! rotate about the z axis DO NOT INVERT!!
322:             DO J1=1,NATOMS322:             DO J1=1,NATOMS
330:          GOTO 20330:          GOTO 20
331:       ENDIF331:       ENDIF
332: 332: 
333:       COST=T1(3*(JDO-1)+1)/DIST(JDO)333:       COST=T1(3*(JDO-1)+1)/DIST(JDO)
334:       SINT=T1(3*(JDO-1)+2)/DIST(JDO)334:       SINT=T1(3*(JDO-1)+2)/DIST(JDO)
335: C     PRINT '(A,4G20.10)','T1(3*(JDO-1)+2),COST,SINT,1=',T1(3*(JDO-1)+2),COST,SINT,COST**2+SINT**2335: C     PRINT '(A,4G20.10)','T1(3*(JDO-1)+2),COST,SINT,1=',T1(3*(JDO-1)+2),COST,SINT,COST**2+SINT**2
336:       IF (ABS(COST**2+SINT**2-1.0D0).GT.2.0D-3) THEN336:       IF (ABS(COST**2+SINT**2-1.0D0).GT.2.0D-3) THEN
337:          WRITE(*, '(A,G20.10)') 'ERROR - in ROTXZ cos**2+sin**2=',COST**2+SINT**2337:          WRITE(*, '(A,G20.10)') 'ERROR - in ROTXZ cos**2+sin**2=',COST**2+SINT**2
338:          STOP338:          STOP
339:       ENDIF339:       ENDIF
340: !     RVEC(1)=0.0D0340:       RVEC(1)=0.0D0
341: !     RVEC(2)=0.0D0341:       RVEC(2)=0.0D0
342: !     RVEC(3)=1.0D0342:       RVEC(3)=1.0D0
343:       DO J1=1,NATOMS343:       DO J1=1,NATOMS
344: !        IF (DIST(J1).NE.0.0D0) THEN344:          IF (DIST(J1).NE.0.0D0) THEN
345:             J2=3*(J1-1)345:             RDOTN=T1(3*(J1-1)+1)*RVEC(1)+T1(3*(J1-1)+2)*RVEC(2)+T1(3*(J1-1)+3)*RVEC(3)
346: !           RDOTN=T1(J2+3)*(1-COST)346:             TX=T1(3*(J1-1)+1)*COST + RVEC(1)*RDOTN*(1-COST)+(T1(3*(J1-1)+2)*RVEC(3)-T1(3*(J1-1)+3)*RVEC(2))*SINT
347:             TX=T1(J2+1)*COST + T1(J2+2)*SINT347:             TY=T1(3*(J1-1)+2)*COST + RVEC(2)*RDOTN*(1-COST)+(T1(3*(J1-1)+3)*RVEC(1)-T1(3*(J1-1)+1)*RVEC(3))*SINT
348:             TY=T1(J2+2)*COST - T1(J2+1)*SINT348:             TZ=T1(3*(J1-1)+3)*COST + RVEC(3)*RDOTN*(1-COST)+(T1(3*(J1-1)+1)*RVEC(2)-T1(3*(J1-1)+2)*RVEC(1))*SINT
349:             T1(J2+1)=TX349:             T1(3*(J1-1)+1)=TX
350:             T1(J2+2)=TY350:             T1(3*(J1-1)+2)=TY
351: !        ENDIF351:             T1(3*(J1-1)+3)=TZ
 352:          ENDIF
352:       ENDDO353:       ENDDO
353:       RDOTN=XVEC(3)354:       RDOTN=XVEC(1)*RVEC(1)+XVEC(2)*RVEC(2)+XVEC(3)*RVEC(3)
354:       TX=XVEC(1)*COST+XVEC(2)*SINT355:       TX=XVEC(1)*COST + RVEC(1)*RDOTN*(1-COST)+(XVEC(2)*RVEC(3)-XVEC(3)*RVEC(2))*SINT
355:       TY=XVEC(2)*COST-XVEC(1)*SINT356:       TY=XVEC(2)*COST + RVEC(2)*RDOTN*(1-COST)+(XVEC(3)*RVEC(1)-XVEC(1)*RVEC(3))*SINT
356:       XVEC(1)=TX; XVEC(2)=TY357:       TZ=XVEC(3)*COST + RVEC(3)*RDOTN*(1-COST)+(XVEC(1)*RVEC(2)-XVEC(2)*RVEC(1))*SINT
357:       RDOTN=YVEC(3)358:       XVEC(1)=TX; XVEC(2)=TY; XVEC(3)=TZ
358:       TX=YVEC(1)*COST + YVEC(2)*SINT359:       RDOTN=YVEC(1)*RVEC(1)+YVEC(2)*RVEC(2)+YVEC(3)*RVEC(3)
359:       TY=YVEC(2)*COST - YVEC(1)*SINT360:       TX=YVEC(1)*COST + RVEC(1)*RDOTN*(1-COST)+(YVEC(2)*RVEC(3)-YVEC(3)*RVEC(2))*SINT
360:       YVEC(1)=TX; YVEC(2)=TY361:       TY=YVEC(2)*COST + RVEC(2)*RDOTN*(1-COST)+(YVEC(3)*RVEC(1)-YVEC(1)*RVEC(3))*SINT
361:       RDOTN=ZVEC(3)362:       TZ=YVEC(3)*COST + RVEC(3)*RDOTN*(1-COST)+(YVEC(1)*RVEC(2)-YVEC(2)*RVEC(1))*SINT
362:       TX=ZVEC(1)*COST + ZVEC(2)*SINT363:       YVEC(1)=TX; YVEC(2)=TY; YVEC(3)=TZ
363:       TY=ZVEC(2)*COST - ZVEC(1)*SINT364:       RDOTN=ZVEC(1)*RVEC(1)+ZVEC(2)*RVEC(2)+ZVEC(3)*RVEC(3)
364:       ZVEC(1)=TX; ZVEC(2)=TY365:       TX=ZVEC(1)*COST + RVEC(1)*RDOTN*(1-COST)+(ZVEC(2)*RVEC(3)-ZVEC(3)*RVEC(2))*SINT
 366:       TY=ZVEC(2)*COST + RVEC(2)*RDOTN*(1-COST)+(ZVEC(3)*RVEC(1)-ZVEC(1)*RVEC(3))*SINT
 367:       TZ=ZVEC(3)*COST + RVEC(3)*RDOTN*(1-COST)+(ZVEC(1)*RVEC(2)-ZVEC(2)*RVEC(1))*SINT
 368:       ZVEC(1)=TX; ZVEC(2)=TY; ZVEC(3)=TZ
365: 369: 
366: 20    CONTINUE370: 20    CONTINUE
367: 371: 
368:       PROJ=0.0D0372:       PROJ=0.0D0
369: !373: !
370: ! For possibly sloppy alignments with LOCALPERMDIST it may be important to374: ! For possibly sloppy alignments with LOCALPERMDIST it may be important to
371: ! increase the cutoff and use CUTOFF1 ! DJW 27/8/09375: ! increase the cutoff and use CUTOFF1 ! DJW 27/8/09
372: !376: !
373:       DO J1=1,NATOMS377:       DO J1=1,NATOMS
 378: C        IF (T1(3*(J1-1)+3).GT.1.0D-2) PROJ=PROJ+T1(3*(J1-1)+1)
374:          IF (T1(3*(J1-1)+3).GT.CUTOFF1) PROJ=PROJ+T1(3*(J1-1)+1)379:          IF (T1(3*(J1-1)+3).GT.CUTOFF1) PROJ=PROJ+T1(3*(J1-1)+1)
375:       ENDDO380:       ENDDO
376: C     PRINT '(I6)',NATOMS381: C     PRINT '(I6)',NATOMS
377: C     PRINT '(A,I6,G20.10)','after rotation atom,PROJ=',JDO,PROJ382: C     PRINT '(A,I6,G20.10)','after rotation atom,PROJ=',JDO,PROJ
378: C     WRITE(*,'(A2,2X,3G20.10)') ('LA',T1(3*(J3-1)+1),T1(3*(J3-1)+2),T1(3*(J3-1)+3),J3=1,NATOMS)383: C     WRITE(*,'(A2,2X,3G20.10)') ('LA',T1(3*(J3-1)+1),T1(3*(J3-1)+2),T1(3*(J3-1)+3),J3=1,NATOMS)
379: 384: 
380:       RETURN385:       RETURN
381:       END386:       END


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0