hdiff output

r33514/congrad.f90 2017-11-23 14:30:19.580425705 +0000 r33513/congrad.f90 2017-11-23 14:30:20.996444533 +0000
 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, ECON, EREP, ESPRING, &
 21:   &  CONVERGECONTEST, CONVERGEREPTEST, FCONTEST, FREPTEST, QCIINTREPMINSEP 21:   &  CONVERGECONTEST, CONVERGEREPTEST
 22: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG 22: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
 23: USE PORFUNCS 23: USE PORFUNCS
 24: IMPLICIT NONE 24: IMPLICIT NONE
 25:             25:            
 26: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF 26: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF
 27: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF, FMAX, FMIN 27: DOUBLE PRECISION :: ETOTAL, RMS, EMAX, EMAXNOFF
 28: INTEGER JJMAX(INTIMAGE+2) 28: INTEGER JJMAX(INTIMAGE+2)
 29: DOUBLE PRECISION  EEMAX(INTIMAGE+2) 29: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
 30: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1 30: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3) 31: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
 32: DOUBLE PRECISION DUMMY, REPGRAD(3), D12, DSQ2, DSQ1, DSQI 32: DOUBLE PRECISION DUMMY, REPGRAD(3), D12, DSQ2, DSQ1, DSQI
 33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2) 33: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
 34: LOGICAL NOINT 34: LOGICAL NOINT
 35: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL 35: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL
 36: LOGICAL IMGFREEZE(INTIMAGE) 36: LOGICAL IMGFREEZE(INTIMAGE)
 37: DOUBLE PRECISION DPLUS, SPGRAD(3), DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq,CUTMAX,DISTMAX 37: DOUBLE PRECISION DPLUS, SPGRAD(3), DCUT, r1amr1bdr2amr2b,r1apr2bmr2amr1bsq,CUTMAX,DISTMAX
 42: REPE(1:INTIMAGE+2)=0.0D0 42: REPE(1:INTIMAGE+2)=0.0D0
 43: REPEINT(1:INTIMAGE+2)=0.0D0 43: REPEINT(1:INTIMAGE+2)=0.0D0
 44: NREPINT(1:INTIMAGE+2)=0 44: NREPINT(1:INTIMAGE+2)=0
 45: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 45: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0
 46: ECON=0.0D0; EREP=0.0D0 46: ECON=0.0D0; EREP=0.0D0
 47: NINTMIN=0 47: NINTMIN=0
 48: NINTMIN2=0 48: NINTMIN2=0
 49: MYUNIT=6 49: MYUNIT=6
 50:  50: 
 51: EMAX=-1.0D200 51: EMAX=-1.0D200
 52: FMAX=-1.0D200 
 53: FMIN=1.0D200 
 54: EEMAX(1:INTIMAGE+2)=-1.0D200 52: EEMAX(1:INTIMAGE+2)=-1.0D200
 55: JJMAX(1:INTIMAGE+2)=-1 53: JJMAX(1:INTIMAGE+2)=-1
 56: JMAX=-1 54: JMAX=-1
 57: IMAX=-1 55: IMAX=-1
 58: EMAXNOFF=-1.0D200 56: EMAXNOFF=-1.0D200
 59: JMAXNOFF=-1 57: JMAXNOFF=-1
 60: IMAXNOFF=-1 58: IMAXNOFF=-1
 61: IF (INTCONSTRAINTDEL.EQ.0.0D0) GOTO 531 59: IF (INTCONSTRAINTDEL.EQ.0.0D0) GOTO 531
 62: ! 60: !
 63: !  Constraint energy and forces. 61: !  Constraint energy and forces.
133:                EMAXNOFF=DUMMY131:                EMAXNOFF=DUMMY
134:             ENDIF132:             ENDIF
135:          ENDIF133:          ENDIF
136:          IF (DUMMY.GT.EEMAX(J1)) THEN134:          IF (DUMMY.GT.EEMAX(J1)) THEN
137:             JJMAX(J1)=J2135:             JJMAX(J1)=J2
138:             EEMAX(J1)=DUMMY136:             EEMAX(J1)=DUMMY
139:          ENDIF137:          ENDIF
140:          CONE(J1)=CONE(J1)+DUMMY138:          CONE(J1)=CONE(J1)+DUMMY
141:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)139:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
142:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)140:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
143:          DUMMY2=MINVAL(REPGRAD) 
144:          IF (DUMMY2.LT.FMIN) FMIN=DUMMY 
145:          DUMMY2=MAXVAL(REPGRAD) 
146:          IF (DUMMY2.GT.FMAX) FMAX=DUMMY 
147:       ENDIF141:       ENDIF
148: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1)142: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1)
149:    ENDDO143:    ENDDO
150: ENDDO144: ENDDO
151: IF (-FMIN.GT.FMAX) FMAX=-FMIN 
152: FCONTEST=FMAX 
153: IF (JMAX.GT.0) THEN145: IF (JMAX.GT.0) THEN
154:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,3G15.5,A,G15.5)') ' congrad> Highest constraint for any image in ',IMAX, &146:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,3G15.5)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &
155:  & ' con ',JMAX, ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX, &147:  & ' constraint ',JMAX, &
156:  & ' max grad=',FMAX148:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX
157: 149: 
158: ENDIF150: ENDIF
159: !  WRITE(*,'(A,3L5)') 'congrad> here C atomactive 2113 2115 2123=',ATOMACTIVE(2113),ATOMACTIVE(2115),ATOMACTIVE(2123) 
160: CONVERGECONTEST=EMAX151: CONVERGECONTEST=EMAX
161: ! IF (JMAXNOFF.GT.0) THEN152: IF (JMAXNOFF.GT.0) THEN
162: !    WRITE(*,'(A,I6,A,I6,A,2I8,A,G20.10,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &153:    WRITE(*,'(A,I6,A,I6,A,2I8,A,G20.10,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &
163: !  & ' constraint ',JMAXNOFF, &154:  & ' constraint ',JMAXNOFF, &
164: !  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' off=',NCONOFF155:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' off=',NCONOFF
165: ! ELSEIF (JMAX.GT.0) THEN156: ELSEIF (JMAX.GT.0) THEN
166: !    JMAXNOFF=JMAX157:    JMAXNOFF=JMAX
167: !    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'158:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'
168: !    CONOFFTRIED(1:INTCONMAX)=.FALSE.159:    CONOFFTRIED(1:INTCONMAX)=.FALSE.
169: !    NCONOFF=0160:    NCONOFF=0
170: ! ENDIF161: ENDIF
171: ! JMAXCON=JMAXNOFF162: JMAXCON=JMAXNOFF
172: 531 CONTINUE163: 531 CONTINUE
173: 164: 
174: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes165: ! GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
175: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes166: ! GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
176: 167: 
177: ! INTCONST=INTCONSTRAINREPCUT**13168: ! INTCONST=INTCONSTRAINREPCUT**13
178: 169: 
179: EMAX=-1.0D200170: EMAX=-1.0D200
180: EEMAX(1:INTIMAGE+2)=-1.0D200171: EEMAX(1:INTIMAGE+2)=-1.0D200
181: JJMAX(1:INTIMAGE+2)=-1172: JJMAX(1:INTIMAGE+2)=-1
212: !203: !
213: ! Squared distance between atoms A and B for theta=0 - distance in image 2204: ! Squared distance between atoms A and B for theta=0 - distance in image 2
214: !205: !
215:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2206:       DSQ2=G2(1)**2 + G2(2)**2 + G2(3)**2
216: !207: !
217: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1208: ! Squared distance between atoms A and B for theta=Pi/2 - distance in image 1
218: !209: !
219:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2210:       DSQ1=G1(1)**2 + G1(2)**2 + G1(3)**2
220:       DCUT=NREPCUT(J2)**2211:       DCUT=NREPCUT(J2)**2
221:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff212:       IF ((DSQ1.GT.DCUT).AND.(DSQ2.GT.DCUT)) CYCLE ! don't look for an internal minimum if both repulsions outside cutoff
222:       IF (ABS(NREPI(J2)-NREPJ(J2)).GT.QCIINTREPMINSEP) THEN ! don't check for internal minimum in distance - atoms too close for chain crossing.213:       r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3)
223:          r1apr2bmr2amr1bsq=0.0D0 
224:       ELSE 
225: !        WRITE(*,'(A,I6,A,I6,A,2G20.10,A,G20.10)') 'distances in ',J1-1,' and ',J1,' are ',SQRT(DSQ1),SQRT(DSQ2),' cutoff=',NREPCUT(J2) 
226: !        WRITE(*,'(A,I6,A,2I6,A,6G15.7)') 'image ',J1-1,' atoms ',NREPI(J2),NREPJ(J2),' coords ',XYZ(NI1+1:NI1+3),XYZ(NJ1+1:NJ1+3) 
227: !        WRITE(*,'(A,I6,A,2I6,A,6G15.7)') 'image ',J1  ,' atoms ',NREPI(J2),NREPJ(J2),' coords ',XYZ(NI2+1:NI2+3),XYZ(NJ2+1:NJ2+3) 
228:          r1amr1bdr2amr2b=G1(1)*G2(1)+G1(2)*G2(2)+G1(3)*G2(3) 
229: !        WRITE(*,'(A,G20.10)') 'dot product (r1I-r1J).(r2I-r2J)=',r1amr1bdr2amr2b 
230: !214: !
231: ! Is the denominator of the d^2 internal minimum term close to zero?215: ! Is there an internal extremum?
232: !216: !
233:          r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b217:       r1apr2bmr2amr1bsq=DSQ1+DSQ2-2.0D0*r1amr1bdr2amr2b
234: !        WRITE(*,'(A,G20.10)') 'internal minimum test d(J1-1)^2+d(J1)^2-2(r1I-r1J).(r2I-r2J)=',r1apr2bmr2amr1bsq 
235:       ENDIF 
236:       IF (r1apr2bmr2amr1bsq.LT.1.0D-50) THEN218:       IF (r1apr2bmr2amr1bsq.LT.1.0D-50) THEN
237: !        D1=1.0D100; D2=1.0D100;219: !        D1=1.0D100; D2=1.0D100;
238:          NOINT=.TRUE.220:          NOINT=.TRUE.
239:          D1=SQRT(DSQ1)221:          D1=SQRT(DSQ1)
240:          D2=SQRT(DSQ2)222:          D2=SQRT(DSQ2)
241:          G2(1:3)=G2(1:3)/D2223:          G2(1:3)=G2(1:3)/D2
242:          G1(1:3)=G1(1:3)/D1224:          G1(1:3)=G1(1:3)/D1
243:       ELSE225:       ELSE
244: !        WRITE(*,'(A)') 'calling MINMAXD2R' 
245:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)226:          CALL MINMAXD2R(D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2),r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)
246:       ENDIF227:       ENDIF
247: !     WRITE(*,'(A,8G15.7)') 'D1,D2,G2,G1=',D1,D2,G2(1:3),G1(1:3) 
248: !228: !
249: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.229: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
250: !230: !
251: ! Multiply energy and gradient by NREPCUT**2/INTCONSTRAINTREP to put values on a common scale.231: ! Multiply energy and gradient by NREPCUT**2/INTCONSTRAINTREP to put values on a common scale.
252: !232: !
253: !     WRITE(*,'(A,2G20.10)') 'D2,NREPCUT=',D2,NREPCUT(J2) 
254:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1233:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
255: ! !        D12=DSQ2**6234: ! !        D12=DSQ2**6
256:          D12=DSQ2235:          D12=DSQ2
257: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)236: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
258: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)237: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
259: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)238: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
260: !        DUMMY=(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2239: !        DUMMY=(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2
261:          DUMMY=NREPCUT(J2)**2/D12+2.0D0*D2/NREPCUT(J2)-3.0D0  240:          DUMMY=NREPCUT(J2)**2/D12+2.0D0*D2/NREPCUT(J2)-3.0D0  
262:          EEE(J1)=EEE(J1)+DUMMY241:          EEE(J1)=EEE(J1)+DUMMY
263: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN242: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
280: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)259: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
281: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)260: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
282: !261: !
283: !        DUMMY=-2.0D0*(1.0D0/(D2*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2262: !        DUMMY=-2.0D0*(1.0D0/(D2*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2
284:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(D2*D12)-1.0D0/NREPCUT(J2))263:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(D2*D12)-1.0D0/NREPCUT(J2))
285:          REPGRAD(1:3)=DUMMY*G2(1:3)264:          REPGRAD(1:3)=DUMMY*G2(1:3)
286: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &265: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
287: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)266: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
288:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)267:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
289:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)268:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
290:          DUMMY2=MINVAL(REPGRAD) 
291:          IF (DUMMY2.LT.FMIN) THEN 
292:             FMIN=DUMMY2 
293: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
294:          ENDIF 
295:          DUMMY2=MAXVAL(REPGRAD) 
296:          IF (DUMMY2.GT.FMAX) THEN 
297:             FMAX=DUMMY2 
298: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
299:          ENDIF 
300: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
301:       ENDIF269:       ENDIF
302: !270: !
303: ! For internal minima we are counting edges. 271: ! For internal minima we are counting edges. 
304: ! Edge J1 is between images J1-1 and J1, starting from J1=2.272: ! Edge J1 is between images J1-1 and J1, starting from J1=2.
305: ! Energy contributions are shared evenly, except for273: ! Energy contributions are shared evenly, except for
306: ! edge 1, which was assigned to image 2, and edge INTIMAGE+1, which274: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which
307: ! was assigned to image INTIMAGE+1. Gradients are set to zero for the end images.275: ! is assigned to image INTIMAGE+1. Gradients are set to zero for
308: ! 20/11/17 - changed to turn off internal minimum checks involving the fixed endpoints. DJW276: ! the end images.
309: !277: !
310:       DUMMY=0.0D0278:       DUMMY=0.0D0
311: !     WRITE(*,'(A,2G20.10)') 'DINT,NREPCUT=',DINT,NREPCUT(J2)279:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN
312:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2)).AND.(J1.NE.2).AND.(J1.NE.INTIMAGE+2)) THEN 
313:          NINTMIN2=NINTMIN2+1280:          NINTMIN2=NINTMIN2+1
314: !        D12=DSQI**6281: !        D12=DSQI**6
315:          D12=DSQI282:          D12=DSQI
316: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)283: ! !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
317: !        DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)284: !        DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
318: !        DUMMY=INTMINFAC*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2285: !        DUMMY=INTMINFAC*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)*NREPCUT(J2)**2
319:          DUMMY=INTMINFAC*(NREPCUT(J2)**2/D12+2.0D0*DINT/NREPCUT(J2)-3.0D0)286:          DUMMY=INTMINFAC*(NREPCUT(J2)**2/D12+2.0D0*DINT/NREPCUT(J2)-3.0D0)
320:          IF (J1.EQ.2) THEN287:          IF (J1.EQ.2) THEN
321:             EEE(J1)=EEE(J1)+DUMMY288:             EEE(J1)=EEE(J1)+DUMMY
322:             REPEINT(J1)=REPEINT(J1)+DUMMY289:             REPEINT(J1)=REPEINT(J1)+DUMMY
348: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)315: ! !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
349: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)316: !        DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST)
350: !        DUMMY=-2.0D0*(1.0D0/(DINT*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2317: !        DUMMY=-2.0D0*(1.0D0/(DINT*D12)-1.0D0/INTCONST)*NREPCUT(J2)**2
351:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(DINT*D12)-1.0D0/NREPCUT(J2))318:          DUMMY=-2.0D0*(NREPCUT(J2)**2/(DINT*D12)-1.0D0/NREPCUT(J2))
352:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)319:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3)
353: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &320: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
354: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)321: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
355: !322: !
356: ! Gradient contributions for image J1-1323: ! Gradient contributions for image J1-1
357: !324: !
358:  
359:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)325:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
360:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)326:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
361:          DUMMY2=MINVAL(REPGRAD)327:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)
362:          IF (DUMMY2.LT.FMIN) THEN 
363:             FMIN=DUMMY2 
364: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
365:          ENDIF 
366:          DUMMY2=MAXVAL(REPGRAD) 
367:          IF (DUMMY2.GT.FMAX) THEN 
368:             FMAX=DUMMY2 
369: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
370:          ENDIF 
371: !        WRITE(*,'(A,2I6,5G20.10)') 'image J1-1 J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
372:  
373: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &328: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
374: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)329: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
375: !330: !
376: ! Gradient contributions for image J1331: ! Gradient contributions for image J1
377: !332: !
378:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3) 
379:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)333:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
380:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)334:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
381:          DUMMY2=MINVAL(REPGRAD) 
382:          IF (DUMMY2.LT.FMIN) THEN 
383:             FMIN=DUMMY2 
384: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
385:          ENDIF 
386:          DUMMY2=MAXVAL(REPGRAD) 
387:          IF (DUMMY2.GT.FMAX) THEN 
388:             FMAX=DUMMY2 
389: !           WRITE(*,'(A,2I6,5G20.10)') 'J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
390:          ENDIF 
391: !        WRITE(*,'(A,2I6,5G20.10)') 'image J1   J1,J2,REPGRAD,FMIN,FMAX=',J1,J2,REPGRAD(1:3),FMIN,FMAX 
392:       ENDIF335:       ENDIF
393:    ENDDO336:    ENDDO
394: ENDDO337: ENDDO
395: IF (-FMIN.GT.FMAX) FMAX=-FMIN 
396: FREPTEST=FMAX 
397: IF (JMAX.GT.0) THEN338: IF (JMAX.GT.0) THEN
398:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G15.5,A,2G15.5,A,G15.5)') ' congrad> Highest repulsion  for any image in ',IMAX, &  339:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G20.10,A,3G20.10)') ' congrad> Highest repulsive  contribution for any image in image ',IMAX, &  
399:  &  ' ind ',JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX),' value=',EMAX,' d,cutoff=',DISTMAX,CUTMAX, &340:  &  ' pair index ', &
400:  &  ' max grad=',FMAX341:  &                                JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX),' value=',EMAX,' d,cutoff=',DISTMAX,CUTMAX 
401: ENDIF342: ENDIF
402: CONVERGEREPTEST=EMAX343: CONVERGEREPTEST=EMAX
403: 654 CONTINUE344: 654 CONTINUE
404: 345: 
405: !346: !
406: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise347: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
407: ! energy terms between images except for the end points.348: ! energy terms between images except for the end points.
408: !349: !
409: ESPRING=0.0D0350: ESPRING=0.0D0
410: IF (KINT.NE.0.0D0) THEN351: IF (KINT.NE.0.0D0) THEN
424:       ENDDO365:       ENDDO
425:       DPLUS=SQRT(DPLUS)366:       DPLUS=SQRT(DPLUS)
426: !     IF (DPLUS.GT.IMSEPMAX) THEN367: !     IF (DPLUS.GT.IMSEPMAX) THEN
427: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2368: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2
428:          DUMMY=KINT*0.5D0*DPLUS**2369:          DUMMY=KINT*0.5D0*DPLUS**2
429:          IF (DUMMY.GT.EMAX) THEN370:          IF (DUMMY.GT.EMAX) THEN
430:             IMAX=J1371:             IMAX=J1
431:             EMAX=DUMMY372:             EMAX=DUMMY
432:          ENDIF373:          ENDIF
433:          ESPRING=ESPRING+DUMMY374:          ESPRING=ESPRING+DUMMY
434: !        IF (J1.EQ.1) THEN 
435: !           EEE(2)=EEE(2)+DUMMY 
436: !        ELSEIF (J1.EQ.INTIMAGE+1) THEN 
437: !           EEE(INTIMAGE+1)=EEE(INTIMAGE+1)+DUMMY 
438: !        ELSE 
439: !           EEE(J1)=EEE(J1)+DUMMY/2.0D0 
440: !           EEE(J1+1)=EEE(J1+1)+DUMMY/2.0D0 
441: !        ENDIF 
442: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS375: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS
443:          DUMMY=KINT376:          DUMMY=KINT
444:          DO J2=1,NATOMS377:          DO J2=1,NATOMS
445:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 378:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
446:                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))379:                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))
447:                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)380:                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)
448:                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)381:                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)
449:             ENDIF382:             ENDIF
450:          ENDDO383:          ENDDO
451: !     ENDIF384: !     ENDIF
514:    ENDIF447:    ENDIF
515:    IF (REPE(J1).GT.MAXINT) THEN448:    IF (REPE(J1).GT.MAXINT) THEN
516:       MAXINT=REPE(J1)449:       MAXINT=REPE(J1)
517:       NMAXINT=J1450:       NMAXINT=J1
518:    ENDIF451:    ENDIF
519: ENDDO452: ENDDO
520: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT453: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> largest  internal energy=',MAXINT,' for image ',NMAXINT
521: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> smallest internal energy=',MININT,' for image ',NMININT454: ! IF (DEBUG) WRITE(*, '(A,G20.10,A,2I6)') ' congrad> smallest internal energy=',MININT,' for image ',NMININT
522: ! IF (DEBUG) WRITE(*, '(A,2I6)') ' congrad> number of internal minima=',NINTMIN,NINTMIN2455: ! IF (DEBUG) WRITE(*, '(A,2I6)') ' congrad> number of internal minima=',NINTMIN,NINTMIN2
523: 456: 
524: ! STOP !!! DEBUG DJW 
525:  
526: END SUBROUTINE CONGRAD457: END SUBROUTINE CONGRAD
527: 458: 
528: SUBROUTINE MINMAXD2(D2,D1,DINT,DSQ2,DSQ1,G1,G2,G1INT,G2INT,NOINT,DEBUG,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)459: SUBROUTINE MINMAXD2(D2,D1,DINT,DSQ2,DSQ1,G1,G2,G1INT,G2INT,NOINT,DEBUG,r1amr1bdr2amr2b,r1apr2bmr2amr1bsq)
529: USE KEY, ONLY : CHECKCONINT460: USE KEY, ONLY : CHECKCONINT
530: IMPLICIT NONE461: IMPLICIT NONE
531: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT462: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1,DINT
532: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)463: DOUBLE PRECISION G1(3),G2(3),G1INT(3),G2INT(3)
533: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq464: DOUBLE PRECISION DSQ2, DSQ1, DSQI, r1apr2bmr2amr1bsq, r1amr1bsq, r2amr2bsq
534: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY, DUMMY2465: DOUBLE PRECISION r1amr1bdr2amr2b, r1amr1bdr2amr2bsq, DUMMY, DUMMY2
535: LOGICAL NOINT, DEBUG466: LOGICAL NOINT, DEBUG
883:       ENDIF814:       ENDIF
884:    ENDDO815:    ENDDO
885: ENDDO816: ENDDO
886: IF (JMAX.GT.0) THEN817: IF (JMAX.GT.0) THEN
887:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G20.10,A,3G20.10)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &818:    WRITE(*,'(A,I6,A,I6,A,2I6,A,G20.10,A,3G20.10)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &
888:  & ' constraint ',JMAX, &819:  & ' constraint ',JMAX, &
889:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX820:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' value=',EMAX,' d,ref,cutoff=',CONDMAX,CONREFMAX,CONCUTMAX
890: 821: 
891: ENDIF822: ENDIF
892: CONVERGECONTEST=EMAX/INTCONSTRAINTDEL823: CONVERGECONTEST=EMAX/INTCONSTRAINTDEL
893: ! IF (JMAXNOFF.GT.0) THEN824: IF (JMAXNOFF.GT.0) THEN
894: !    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &825:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &
895: !  & ' constraint ',JMAXNOFF, &826:  & ' constraint ',JMAXNOFF, &
896: !  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF827:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF
897: ! ELSEIF (JMAX.GT.0) THEN828: ELSEIF (JMAX.GT.0) THEN
898: !    JMAXNOFF=JMAX829:    JMAXNOFF=JMAX
899: !    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> *** Using highest constraint contribution for any image, setting NCONOFF to 0'830:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> *** Using highest constraint contribution for any image, setting NCONOFF to 0'
900: !    CONOFFTRIED(1:INTCONMAX)=.FALSE.831:    CONOFFTRIED(1:INTCONMAX)=.FALSE.
901: !    NCONOFF=0832:    NCONOFF=0
902: ! ENDIF833: ENDIF
903: JMAXCON=JMAXNOFF834: JMAXCON=JMAXNOFF
904: 835: 
905: ! INTCONST=INTCONSTRAINREPCUT**13836: ! INTCONST=INTCONSTRAINREPCUT**13
906: 837: 
907: EMAX=-1.0D200838: EMAX=-1.0D200
908: EEMAX(1:INTIMAGE+2)=-1.0D200839: EEMAX(1:INTIMAGE+2)=-1.0D200
909: JJMAX(1:INTIMAGE+2)=-1840: JJMAX(1:INTIMAGE+2)=-1
910: JMAX=-1841: JMAX=-1
911: IMAX=-1842: IMAX=-1
912: DO J2=1,NNREPULSIVE843: DO J2=1,NNREPULSIVE


r33514/intlbfgs.f90 2017-11-23 14:30:19.864429481 +0000 r33513/intlbfgs.f90 2017-11-23 14:30:21.288448415 +0000
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, & 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, &
 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, & 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, &
 32:      & QCIADDREP, QCIXYZ, WHOLEDNEB, QCIIMAGE, FROZEN, QCIRESTART, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, & 32:      & QCIADDREP, QCIXYZ, WHOLEDNEB, QCIIMAGE, FROZEN, QCIRESTART, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, &
 33:      & PERMDIST, LOCALPERMCUT, QCILPERMDIST, QCIPDINT, QCIPERMCUT, QCIAMBERT, BONDS, DOBACK, & 33:      & PERMDIST, LOCALPERMCUT, QCILPERMDIST, QCIPDINT, QCIPERMCUT, QCIAMBERT, BONDS, DOBACK, &
 34:      & QCIRESET, QCIRESETINT1, QCIRESETINT2, JMAXCON, NCONOFF, EREP, ECON, ESPRING, CONVERGECONTEST, CONVERGEREPTEST, & 34:      & QCIRESET, QCIRESETINT1, QCIRESETINT2, JMAXCON, NCONOFF, EREP, ECON, ESPRING, CONVERGECONTEST, CONVERGEREPTEST
 35:      & FCONTEST, FREPTEST 
 36: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3 35: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3
 37: USE MODCHARMM, ONLY : CHRMMT 36: USE MODCHARMM, ONLY : CHRMMT
 38: USE CHIRALITY 37: USE CHIRALITY
 39:  38: 
 40: IMPLICIT NONE  39: IMPLICIT NONE 
 41:  40: 
 42: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points 41: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points
 43: INTEGER D, U 42: INTEGER D, U
 44: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE, NEIGHBOUR_COORDS(12), CENTRE_COORDS(3) 43: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE, NEIGHBOUR_COORDS(12), CENTRE_COORDS(3)
 45: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ 44: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ
311:    READ(LUNIT,*) NACTIVE310:    READ(LUNIT,*) NACTIVE
312:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NACTIVE,' active atoms'311:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NACTIVE,' active atoms'
313:    READ(LUNIT,*) KINT, INTCONSTRAINTREP312:    READ(LUNIT,*) KINT, INTCONSTRAINTREP
314:    WRITE(*,'(A,2G20.10)') ' intlbfgs> spring constant and repulsive prefactor: ',KINT, INTCONSTRAINTREP313:    WRITE(*,'(A,2G20.10)') ' intlbfgs> spring constant and repulsive prefactor: ',KINT, INTCONSTRAINTREP
315:    WRITE(*,'(A)') ' intlbfgs> reading turnonorder'314:    WRITE(*,'(A)') ' intlbfgs> reading turnonorder'
316:    READ(LUNIT,*) TURNONORDER(1:NACTIVE)315:    READ(LUNIT,*) TURNONORDER(1:NACTIVE)
317: !  WRITE(*,'(12I8)') TURNONORDER(1:NACTIVE)316: !  WRITE(*,'(12I8)') TURNONORDER(1:NACTIVE)
318:    WRITE(*,'(A)') ' intlbfgs> reading atomactive'317:    WRITE(*,'(A)') ' intlbfgs> reading atomactive'
319:    READ(LUNIT,*) ATOMACTIVE(1:NATOMS)318:    READ(LUNIT,*) ATOMACTIVE(1:NATOMS)
320: !  WRITE(*,'(12L5)') ATOMACTIVE(1:NATOMS) 319: !  WRITE(*,'(12L5)') ATOMACTIVE(1:NATOMS) 
321: !  WRITE(*,'(A,3L5)') 'intlbfgs> here A atomactive 2113 2115 2123=',ATOMACTIVE(2113),ATOMACTIVE(2115),ATOMACTIVE(2123) 
322:    READ(LUNIT,*) NCONSTRAINT320:    READ(LUNIT,*) NCONSTRAINT
323:    WRITE(*,'(A)') ' intlbfgs> reading conactive'321:    WRITE(*,'(A)') ' intlbfgs> reading conactive'
324:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NCONSTRAINT,' constraints'322:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NCONSTRAINT,' constraints'
325:    ALLOCATE(CONACTIVE(NCONSTRAINT))323:    ALLOCATE(CONACTIVE(NCONSTRAINT))
326:    READ(LUNIT,*) CONACTIVE(1:NCONSTRAINT)324:    READ(LUNIT,*) CONACTIVE(1:NCONSTRAINT)
327: !  WRITE(*,'(12L5)') CONACTIVE(1:NCONSTRAINT) 325: !  WRITE(*,'(12L5)') CONACTIVE(1:NCONSTRAINT) 
328: 326: 
329:    ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))327:    ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))
330:    ALLOCATE(CONCUTLOCAL(NCONSTRAINT))328:    ALLOCATE(CONCUTLOCAL(NCONSTRAINT))
331:    CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)329:    CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)
661: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.659: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.
662: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!660: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!
663: !661: !
664: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.662: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.
665: !663: !
666: BESTWORST=1.0D100664: BESTWORST=1.0D100
667: 9876 CONTINUE665: 9876 CONTINUE
668: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)666: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
669: 667: 
670: 668: 
671: !                DIFF=1.0D-7669: !                DIFF=1.0D-6
672: !                PRINT*,'analytic and numerical gradients:'670: !                PRINT*,'analytic and numerical gradients:'
673: !                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)671: !                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
674: !!               DO J1=2,INTIMAGE+1672: !                DO J1=2,INTIMAGE+1
675: !                DO J1=INTIMAGE-1,INTIMAGE+1 
676: !                   DO J2=1,3*NATOMS673: !                   DO J2=1,3*NATOMS
677: !                      IF (.NOT.ATOMACTIVE((J2-1)/3+1)) CYCLE674: !                      IF (.NOT.ATOMACTIVE((J2-1)/3+1)) CYCLE
678: !                      J3=3*NATOMS*(J1-1)+J2675: !                      J3=3*NATOMS*(J1-1)+J2
679: !                      XYZ(J3)=XYZ(J3)+DIFF676: !                      XYZ(J3)=XYZ(J3)+DIFF
680: !                      CALL CONGRAD(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)677: !                      CALL CONGRAD2(NMAXINT,NMININT,EPLUS,XYZ,VPLUS,EEE,IMGFREEZE,RMS)
681: !                      XYZ(J3)=XYZ(J3)-2.0D0*DIFF678: !                      XYZ(J3)=XYZ(J3)-2.0D0*DIFF
682: !                      CALL CONGRAD(NMAXINT,NMININT,EMINUS,XYZ,VMINUS,EEE,IMGFREEZE,RMS)679: !                      CALL CONGRAD2(NMAXINT,NMININT,EMINUS,XYZ,VMINUS,EEE,IMGFREEZE,RMS)
683: !                      XYZ(J3)=XYZ(J3)+DIFF680: !                      XYZ(J3)=XYZ(J3)+DIFF
684: !                      IF ((ABS(GGG(J3)).NE.0.0D0).AND. &681: !                      IF ((ABS(GGG(J3)).NE.0.0D0).AND. &
685: !     &               (ABS(100.0D0*(GGG(J3)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GGG(J3)).GT.1.0D0)) THEN682: !     &               (ABS(100.0D0*(GGG(J3)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GGG(J3)).GT.1.0D0)) THEN
686: !                         WRITE(*,'(A,2I5,2G20.10,L5,A)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &683: !                         WRITE(*,'(A,2I5,2G20.10,L5,A)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &
687: !     & ATOMACTIVE((J2-1)/3+1),'   X'   684: !     & ATOMACTIVE((J2-1)/3+1),'   X'   
688: !                      ELSE685: !                      ELSE
689: !                         WRITE(*,'(A,2I5,2G20.10,L5)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &686: !                         WRITE(*,'(A,2I5,2G20.10,L5)') 'anal and num ',J1,J2,GGG(J3),(EPLUS-EMINUS)/(2.0D0*DIFF),  &
690: !     & ATOMACTIVE((J2-1)/3+1)  687: !     & ATOMACTIVE((J2-1)/3+1)  
691: !                      ENDIF688: !                      ENDIF
692: !                   ENDDO689: !                   ENDDO
740:       NDUMMY=NDUMMY+NPERMSIZE(J1)737:       NDUMMY=NDUMMY+NPERMSIZE(J1)
741:       ENDGROUP(J1)=NDUMMY-1738:       ENDGROUP(J1)=NDUMMY-1
742:    ENDDO739:    ENDDO
743: ENDIF740: ENDIF
744: 741: 
745: 567 CONTINUE742: 567 CONTINUE
746: 743: 
747: DO ! Main do loop with counter NITERDONE, initially set to one744: DO ! Main do loop with counter NITERDONE, initially set to one
748: 745: 
749: !746: !
750: ! Are we stuck? 747: ! Are we stuck? If so, try resetting problem atoms to previous image.
751: !748: !
752: IF (QCIRESET) THEN749: IF (QCIRESET) THEN
753: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN750: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN
754: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1751: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1
755:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN ! .AND.(NITERDONE-NLASTCHANGE.GT.QCIRESETINT1)) THEN752:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN ! .AND.(NITERDONE-NLASTCHANGE.GT.QCIRESETINT1)) THEN
756:       IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)753:       IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
757:       IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)754:       IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
758: !     TURNOFF=.TRUE.755:       TURNOFF=.TRUE.
759: !     IF (EREP.GT.2.0D0*ECON) THEN756: !     IF (EREP.GT.2.0D0*ECON) THEN
760: !        INTCONSTRAINTREP=INTCONSTRAINTREP/2.0D0757: !        INTCONSTRAINTREP=INTCONSTRAINTREP/2.0D0
761: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Interpolation seems to be stuck. Reducing repulsive prefactor to ',INTCONSTRAINTREP758: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Interpolation seems to be stuck. Reducing repulsive prefactor to ',INTCONSTRAINTREP
762: !        TURNOFF=.FALSE.759: !        TURNOFF=.FALSE.
763: !     ENDIF760: !     ENDIF
764: !     IF (ESPRING.GT.2.0D0*ECON) THEN761:       IF (ESPRING.GT.2.0D0*ECON) THEN
765: !        KINT=KINT/2.0D0762:          KINT=KINT/2.0D0
766: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Interpolation seems to be stuck. Reducing spring constant to ',KINT763:          WRITE(*,'(A,G20.10)') 'intlbfgs> Interpolation seems to be stuck. Reducing spring constant to ',KINT
767: !        TURNOFF=.FALSE.764:          TURNOFF=.FALSE.
768: !     ENDIF765:       ENDIF
769: !     IF (TURNOFF) THEN766:       IF (TURNOFF) THEN
770: !        WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Dump images. Turn off worst constraint ',JMAXCON, &767:          WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Dump images. Turn off worst constraint ',JMAXCON, &
771: ! &                        ' total off=',NCONOFF+1768:   &                        ' total off=',NCONOFF+1
772: !        IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN769:          IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN
773: !           WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'770:             WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'
774: !           WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF771:             WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF
775: !           WRITE(*,'(A)') 'CONOFFTRIED:'772:             WRITE(*,'(A)') 'CONOFFTRIED:'
776: !           WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT)773:             WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT)
777: !           STOP774:             STOP
778: !        ENDIF775:          ENDIF
779: !        NCONOFF=NCONOFF+1776:          NCONOFF=NCONOFF+1
780: !        CONOFFLIST(NCONOFF)=JMAXCON777:          CONOFFLIST(NCONOFF)=JMAXCON
781: !        CONACTIVE(JMAXCON)=.FALSE.778:          CONACTIVE(JMAXCON)=.FALSE.
782: !        CONOFFTRIED(JMAXCON)=.TRUE.779:          CONOFFTRIED(JMAXCON)=.TRUE.
783: !     ENDIF780:       ENDIF
784:  
785:       IF (MAX(CONVERGECONTEST,CONVERGEREPTEST).GT.MAXCONE) MAXCONE=MAXCONE*1.05D0 
786:       IF (MAX(FCONTEST,FREPTEST).GT.INTRMSTOL) INTRMSTOL=INTRMSTOL*1.05D0 
787:       WRITE(*,'(A,2G20.10)') 'intlbfgs> Interpolation seems to be stuck. Converge thresholds are now ',MAXCONE,INTRMSTOL 
788: 781: 
789:       NLASTGOODE=NITERDONE782:       NLASTGOODE=NITERDONE
790:       LASTGOODE=ETOTAL783:       LASTGOODE=ETOTAL
791: !     IF (NITERDONE.GT.650) STOP  !!!! DEBUG DJW784: !     IF (NITERDONE.GT.650) STOP  !!!! DEBUG DJW
792: !     STOP785: !     STOP
793:    ENDIF786:    ENDIF
794: ENDIF787: ENDIF
795: 788: 
796: !     IF (NITERDONE.GT.200) STOP  !!!! DEBUG DJW789: !     IF (NITERDONE.GT.200) STOP  !!!! DEBUG DJW
797: ! STOP790: ! STOP
1550:       WRITE(*,'(A,I8,A,G20.10,A,G20.10,A)') ' intlbfgs> Highest image ',JMAXEEE,' energy ',MAXEEE,' is ',ABS(MAXEEE-SUMEEE)/SIGMAEEE, &1543:       WRITE(*,'(A,I8,A,G20.10,A,G20.10,A)') ' intlbfgs> Highest image ',JMAXEEE,' energy ',MAXEEE,' is ',ABS(MAXEEE-SUMEEE)/SIGMAEEE, &
1551:   &                        ' sigma from the mean'1544:   &                        ' sigma from the mean'
1552: !       REMOVEIMAGE=.TRUE.1545: !       REMOVEIMAGE=.TRUE.
1553: !       IF (NITERDONE-NLASTGOODE.LT.20) REMOVEIMAGE=.FALSE.1546: !       IF (NITERDONE-NLASTGOODE.LT.20) REMOVEIMAGE=.FALSE.
1554: !    ENDIF1547: !    ENDIF
1555: 1548: 
1556:    IF (DEBUG) THEN1549:    IF (DEBUG) THEN
1557: !     WRITE(*,'(A,2G20.10)') ' intlbfgs> mean and sigma of image energies=',SUMEEE,SIGMAEEE1550: !     WRITE(*,'(A,2G20.10)') ' intlbfgs> mean and sigma of image energies=',SUMEEE,SIGMAEEE
1558: !     WRITE(*,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &1551: !     WRITE(*,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &
1559: ! &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS1552: ! &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS
1560:       WRITE(*,'(A,I6,5G20.10,I10,I6)') ' intlbfgs> steps: ',NITERDONE,CONVERGECONTEST,CONVERGEREPTEST,FCONTEST,FREPTEST,STEPTOT,NACTIVE,INTIMAGE+21553:       WRITE(*,'(A,I6,4G20.10,I10,I6)') ' intlbfgs> steps: ',NITERDONE,CONVERGECONTEST,CONVERGEREPTEST,RMS,STEPTOT,NACTIVE,INTIMAGE+2
1561:       CALL FLUSH(6)1554:       CALL FLUSH(6)
1562:    ENDIF1555:    ENDIF
1563: 1556: 
1564:    IF (.NOT.SWITCHED) THEN1557:    IF (.NOT.SWITCHED) THEN
1565: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN1558: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN
1566:       IF (.FALSE.) THEN ! no backtracking1559:       IF (.FALSE.) THEN ! no backtracking
1567:          WRITE(*,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE1560:          WRITE(*,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE
1568:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+11561:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+1
1569:          IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.1562:          IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.
1570: !1563: !
1664:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1657:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1665:          ELSE1658:          ELSE
1666:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1659:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1667:          ENDIF1660:          ENDIF
1668:       ENDIF1661:       ENDIF
1669:       LASTGOODE=ETOTAL1662:       LASTGOODE=ETOTAL
1670:    ENDIF1663:    ENDIF
1671: 1664: 
1672:    EXITSTATUS=01665:    EXITSTATUS=0
1673:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1666:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1674: !  IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.(NITERDONE>1).AND.(CONVERGECONTEST.LT.MAXCONE).AND.(CONVERGEREPTEST.LT.MAXCONE)) EXITSTATUS=1 1667:    PRINT *,'switched=',switched
1675:    IF ((.NOT.SWITCHED).AND.(FCONTEST.LT.INTRMSTOL).AND.(FREPTEST.LT.INTRMSTOL).AND.(NITERDONE>1) &1668:    PRINT *,'MAXRMS,INTRMSTOL=',MAXRMS,INTRMSTOL
1676:   &                   .AND.(CONVERGECONTEST.LT.MAXCONE).AND.(CONVERGEREPTEST.LT.MAXCONE)) EXITSTATUS=1 1669:    PRINT *,'CONVERGECONTEST,CONVERGEREPTEST,MAXCONE=',CONVERGECONTEST,CONVERGEREPTEST,MAXCONE
 1670:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.(NITERDONE>1).AND.(CONVERGECONTEST.LT.MAXCONE).AND.(CONVERGEREPTEST.LT.MAXCONE)) EXITSTATUS=1 
1677:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 1671:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 
1678:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21672:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1679:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1673:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
 1674:    PRINT *,'EXITSTATUS=',EXITSTATUS
1680: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX1675: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX
1681: 1676: 
1682:    IF (EXITSTATUS > 0) THEN  1677:    IF (EXITSTATUS > 0) THEN  
1683:       PRINT *,'here A'1678:       PRINT *,'here A'
1684:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1679:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1685: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771680: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1686: !        PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))1681: !        PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))
1687: !        IF (MAXEEE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771682: !        IF (MAXEEE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
 1683:       PRINT *,'here B nactive, natoms=',nactive,natoms
1688:          IF (NACTIVE.LT.NATOMS) THEN 1684:          IF (NACTIVE.LT.NATOMS) THEN 
1689:             ADDATOM=.TRUE.1685:             ADDATOM=.TRUE.
1690:             GOTO 7771686:             GOTO 777
1691:          ENDIF1687:          ENDIF
1692:          CALL MYCPU_TIME(FTIME,.FALSE.)1688:          CALL MYCPU_TIME(FTIME,.FALSE.)
1693:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1689:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1694:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1690:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1695:          CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1691:          CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1696:          CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1692:          CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1697:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1693:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
1883: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2))1879: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2))
1884: 1880: 
1885: FILENAME=XYZFILE1881: FILENAME=XYZFILE
1886: 1882: 
1887: ! IF (NITER.GT.0) THEN1883: ! IF (NITER.GT.0) THEN
1888: !    WRITE(DUMMYS,'(I8)') NITER1884: !    WRITE(DUMMYS,'(I8)') NITER
1889: !    FILENAME='int.' // TRIM(ADJUSTL(DUMMYS)) // '.xyz' ! so that vmd recognises the file type!1885: !    FILENAME='int.' // TRIM(ADJUSTL(DUMMYS)) // '.xyz' ! so that vmd recognises the file type!
1890: ! ENDIF1886: ! ENDIF
1891: LUNIT=GETUNIT()1887: LUNIT=GETUNIT()
1892: OPEN(UNIT=LUNIT,FILE=TRIM(ADJUSTL(FILENAME)),STATUS='replace')1888: OPEN(UNIT=LUNIT,FILE=TRIM(ADJUSTL(FILENAME)),STATUS='replace')
1893: !  WRITE(*,'(A,3L5)') 'intrwg> here B atomactive 2113 2115 2123=',ATOMACTIVE(2113),ATOMACTIVE(2115),ATOMACTIVE(2123) 
1894: DO J2=1,INTIMAGE+21889: DO J2=1,INTIMAGE+2
1895: !  WRITE(LUNIT,'(i4/)') NACTIVE1890: !  WRITE(LUNIT,'(i4/)') NACTIVE
1896:    WRITE(LUNIT,'(i4/)') NATOMS1891:    WRITE(LUNIT,'(i4/)') NATOMS
1897:    DO J3=1,NATOMS1892:    DO J3=1,NATOMS
1898:       IF (ATOMACTIVE(J3)) THEN1893:       IF (ATOMACTIVE(J3)) THEN
1899:          WRITE(LUNIT,'(A5,1X,3F20.10)') 'LA   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  1894:          WRITE(LUNIT,'(A5,1X,3F20.10)') 'LA   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  
1900:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  1895:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  
1901:       ELSE1896:       ELSE
1902:          WRITE(LUNIT,'(A5,1X,3F20.10)') 'DU   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  1897:          WRITE(LUNIT,'(A5,1X,3F20.10)') 'DU   ',XYZ((J2-1)*3*NATOMS+3*(J3-1)+1),XYZ((J2-1)*3*NATOMS+3*(J3-1)+2), &  
1903:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  1898:   &                                                                   XYZ((J2-1)*3*NATOMS+3*(J3-1)+3)  


r33514/key.f90 2017-11-23 14:30:20.140433151 +0000 r33513/key.f90 2017-11-23 14:30:21.588452407 +0000
 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, & 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, &
 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, & 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, &
 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, & 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, &
 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, & 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, &
 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, & 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, &
 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, & 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, &
 25:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS, & 25:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS, &
 26:      &        MLQIN, MLQSTART, MLQOUT, MLQDATA, NMLQ, & 26:      &        MLQIN, MLQSTART, MLQOUT, MLQDATA, NMLQ, &
 27:      &        QCIADDREP, QCIBONDS, QCISECOND, MAXNACTIVE, QCIIMAGE, NADDTARGET, NUMNN, MULTI_COUNT, MULTI_LAST, MULTI_STEP, & 27:      &        QCIADDREP, QCIBONDS, QCISECOND, MAXNACTIVE, QCIIMAGE, NADDTARGET, NUMNN, MULTI_COUNT, MULTI_LAST, MULTI_STEP, &
 28:      &        NDOF, RECCOUNT, MLPPROBPOS, PUSHOPTMAX, MLPNEIGH, QCICYCN, QCIPDINT, NPEAKS, NDISPLACEMENTS, NROTATIONS, & 28:      &        NDOF, RECCOUNT, MLPPROBPOS, PUSHOPTMAX, MLPNEIGH, QCICYCN, QCIPDINT, NPEAKS, NDISPLACEMENTS, NROTATIONS, &
 29:      &        MAX_ANGMOM, BNB_NSTEPS, QUIPZ, QCIRESETINT1, QCIRESETINT2, JMAXCON, NCONOFF, QCIINTREPMINSEP 29:      &        MAX_ANGMOM, BNB_NSTEPS, QUIPZ, QCIRESETINT1, QCIRESETINT2, JMAXCON, NCONOFF
 30:  30: 
 31:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, & 31:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, &
 32:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, & 32:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, &
 33:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, & 33:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, &
 34:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, & 34:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, &
 35:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, & 35:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, &
 36:      &        BFGSSTEP, EFOLSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, & 36:      &        BFGSSTEP, EFOLSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, &
 37:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, & 37:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, &
 38:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, & 38:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, &
 39:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,& 39:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,&
108:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &108:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &
109:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &109:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &
110:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &110:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &
111:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &111:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &
112:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &112:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &
113:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &113:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
114:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &114:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
115:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &115:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
116:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2, &116:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2, &
117:      &        TANHFAC, LJADDCUTOFF,LJADDREFNORM,MAXIMFACTOR, PUSHOPTCONV, QCICYCDIST, QCIPERMCUT, KERNELWIDTH, QUIPLATT(3,3), &117:      &        TANHFAC, LJADDCUTOFF,LJADDREFNORM,MAXIMFACTOR, PUSHOPTCONV, QCICYCDIST, QCIPERMCUT, KERNELWIDTH, QUIPLATT(3,3), &
118:      &        ECON, EREP, ESPRING, CONVERGECONTEST, CONVERGEREPTEST, FCONTEST, FREPTEST118:      &        ECON, EREP, ESPRING, CONVERGECONTEST, CONVERGEREPTEST
119: 119: 
120: !     sf344120: !     sf344
121:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &121:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
122:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &122:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
123:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT123:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT
124:  124:  
125:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)125:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)
126:       DOUBLE PRECISION, ALLOCATABLE :: SHIFTL(:)126:       DOUBLE PRECISION, ALLOCATABLE :: SHIFTL(:)
127:       LOGICAL, ALLOCATABLE :: uniaxarray(:)127:       LOGICAL, ALLOCATABLE :: uniaxarray(:)
128:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)128:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)


r33514/keywords.f 2017-11-23 14:30:20.448437247 +0000 r33513/keywords.f 2017-11-23 14:30:21.876456232 +0000
654:          DOBACK=.FALSE.654:          DOBACK=.FALSE.
655:          DOBACKALL=.FALSE.655:          DOBACKALL=.FALSE.
656:          QCIRESET=.FALSE.656:          QCIRESET=.FALSE.
657:          QCIRESETINT1=300657:          QCIRESETINT1=300
658:          QCIRESETINT2=1000658:          QCIRESETINT2=1000
659:          QCIADDACIDT=.FALSE.659:          QCIADDACIDT=.FALSE.
660:          QCITRILAT=.FALSE.660:          QCITRILAT=.FALSE.
661:          QCIADDREPCUT=1.0D0661:          QCIADDREPCUT=1.0D0
662:          QCIADDREPEPS=1.0D0662:          QCIADDREPEPS=1.0D0
663:          QCINOREPINT=.FALSE.663:          QCINOREPINT=.FALSE.
664:          QCIINTREPMINSEP=20 
665:          MAXNACTIVE=0664:          MAXNACTIVE=0
666: 665: 
667:          FREEZETOL=1.0D-3666:          FREEZETOL=1.0D-3
668:          FLATTESTT=.FALSE.667:          FLATTESTT=.FALSE.
669:          FLATEDIFF=1.0D-6668:          FLATEDIFF=1.0D-6
670:          QCIPERMCHECK=.FALSE.669:          QCIPERMCHECK=.FALSE.
671:          QCIPERMCHECKINT=100670:          QCIPERMCHECKINT=100
672:          INTCONSTRAINTT=.FALSE.671:          INTCONSTRAINTT=.FALSE.
673:          INTCONSTRAINTTOL=0.1D0672:          INTCONSTRAINTTOL=0.1D0
674:          INTCONSTRAINTDEL=10.0D0673:          INTCONSTRAINTDEL=10.0D0
2186: ! input to OPTIM). The alignment is done via rigid-body movements and by considering2185: ! input to OPTIM). The alignment is done via rigid-body movements and by considering
2187: ! permutational isomerizations. The aligned `finish' coordinates are dumped and the2186: ! permutational isomerizations. The aligned `finish' coordinates are dumped and the
2188: ! minimized distance is printed.2187: ! minimized distance is printed.
2189: ! 2188: ! 
2190:          ELSE IF (WORD.EQ.'CLOSESTALIGNMENT') THEN2189:          ELSE IF (WORD.EQ.'CLOSESTALIGNMENT') THEN
2191:             CLOSESTALIGNMENT=.TRUE.2190:             CLOSESTALIGNMENT=.TRUE.
2192:             WRITE(*,*) 'Putting structures into closest alignment, then stopping'2191:             WRITE(*,*) 'Putting structures into closest alignment, then stopping'
2193: ! 2192: ! 
2194: ! Check for internal minimum in constraint terms for INTCONSTRAINT2193: ! Check for internal minimum in constraint terms for INTCONSTRAINT
2195: ! 2194: ! 
2196:          ELSE IF (WORD.EQ.'QCIINTREPMINSEP') THEN 
2197:             CALL READI(QCIINTREPMINSEP) 
2198:             WRITE(*,'(A,G20.10)') ' keyword> Minimum separation in atom index for internal minimum check in repulsion=', 
2199:      &                              QCIINTREPMINSEP   
2200: !  
2201: ! Check for internal minimum in constraint terms for INTCONSTRAINT 
2202: !  
2203:          ELSE IF ((WORD.EQ.'CONCONINT').OR.(WORD.EQ.'QCICONINT')) THEN2195:          ELSE IF ((WORD.EQ.'CONCONINT').OR.(WORD.EQ.'QCICONINT')) THEN
2204:             CHECKCONINT=.TRUE.2196:             CHECKCONINT=.TRUE.
2205:             WRITE(*,'(A,G20.10)') ' keyword> Turning on terms for internal minima in constraints'2197:             WRITE(*,'(A,G20.10)') ' keyword> Turning on terms for internal minima in constraints'
2206: ! 2198: ! 
2207: ! Scaling factor for internal minima in con or rep2199: ! Scaling factor for internal minima in con or rep
2208: ! 2200: ! 
2209:          ELSE IF ((WORD.EQ.'CONINT').OR.(WORD.EQ.'QCIINT')) THEN2201:          ELSE IF ((WORD.EQ.'CONINT').OR.(WORD.EQ.'QCIINT')) THEN
2210:             IF (NITEMS.GT.1) CALL READF(INTMINFAC)2202:             IF (NITEMS.GT.1) CALL READF(INTMINFAC)
2211:             WRITE(*,'(A,G20.10)') ' keyword> Internal minima terms will be scaled by a factor of ',INTMINFAC2203:             WRITE(*,'(A,G20.10)') ' keyword> Internal minima terms will be scaled by a factor of ',INTMINFAC
2212: !2204: !


r33514/lopermdistcon.f90.save 2017-11-23 14:30:20.716440810 +0000 r33513/lopermdistcon.f90.save 2017-11-23 14:30:22.148459854 +0000
  1: !     Copyright (C) 1999-2008 David J. Wales  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/lopermdistcon.f90.save' in revision 33513
  2: !  This file is part of OPTIM. 
  3: ! 
  4: !  OPTIM is free software; you can redistribute it and/or modify 
  5: !  it under the terms of the GNU General Public License as published by 
  6: !  the Free Software Foundation; either version 2 of the License, or 
  7: !  (at your option) any later version. 
  8: ! 
  9: !  OPTIM is distributed in the hope that it will be useful, 
 10: !  but WITHOUT ANY WARRANTY; without even the implied warranty of 
 11: !  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 12: !  GNU General Public License for more details. 
 13: ! 
 14: !  You should have received a copy of the GNU General Public License 
 15: !  along with this program; if not, write to the Free Software 
 16: !  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 17: ! 
 18: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 19: ! 
 20: !  This routine uses local optimal alignment for each group of permutable atoms. 
 21: !  It is intended for use with CHARMM and AMBER. 
 22: !  Overall alignment is based on the locally constrained atoms from QCI. 
 23: !  CONLIST(J1,1:NCONATOM(J1)) contains the list fo constrained atoms for J1, total NCONATOM(J1) 
 24: ! 
 25: !  COORDSA becomes the optimal alignment of the optimal permutation 
 26: !  isomer, but without the permutations. DISTANCE is the residual square distance 
 27: !  for the best alignment with respect to permutation as well as 
 28: !  orientation and centre of mass. 
 29: ! 
 30: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 
 31: !  centre of coordinates of COORDSA will be the same as for COORDSB. 
 32: ! 
 33: SUBROUTINE LOPERMDISTCON(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST,DOGROUP,NMOVE,NEWPERM) 
 34:  
 35: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, GEOMDIFFTOL, & 
 36:   &            NFREEZE, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, NTSITES, & 
 37:   &            LOCALPERMCUT, STOCKT, GMAX, CONVR, EDIFFTOL, INTCONSTRAINTT, & 
 38:   &            LOCALPERMNEIGH, LOCALPERMCUT2, ATOMACTIVE 
 39: USE MODCHARMM,ONLY : CHRMMT 
 40: IMPLICIT NONE 
 41:  
 42: INTEGER, PARAMETER :: MAXIMUMTRIES=10 
 43: INTEGER NATOMS, NPERM, PATOMS, NRB, OPNUM,  NORBIT1, NORBIT2, NCHOOSE2, NCHOOSE1, NTRIES, NORBITB1, NORBITB2, NMOVE, NDMEAN 
 44: INTEGER J3, J4, NDUMMY, LPERM(NATOMS), J1, J2, NOTHER, LPERMBEST(NATOMS), NCHOOSEB1, NCHOOSEB2, & 
 45:         LPERMBESTATOM(NATOMS) 
 46: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), & 
 47:   &              BESTA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS), DIST, DSUM 
 48: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS) 
 49: DOUBLE PRECISION CMXA, CMXB, CMXC, QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES) 
 50: DOUBLE PRECISION ROTA(3,3), ROTINVA(3,3), ROTB(3,3), ROTINVB(3,3), RMATBEST(3,3), TMAT(3,3) 
 51: DOUBLE PRECISION PVEC(3), RTEMP1(3,3), RTEMP2(3,3) 
 52: LOGICAL DEBUG, TWOD, RIGID, BULKT, PITEST, AOK, BOK, ADDED, PERMUTABLE(NATOMS), USEATOM 
 53: DOUBLE PRECISION PDUMMYA(3*NATOMS), PDUMMYB(3*NATOMS), LDISTANCE, DUMMYC(3*NATOMS), XDUMMY, DUMMYD(3*NATOMS), & 
 54:    &             LDBEST(NPERMGROUP), LDBESTATOM 
 55: DOUBLE PRECISION SPDUMMYA(3*NATOMS), SPDUMMYB(3*NATOMS), AINIT, BINIT 
 56: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS) 
 57: DOUBLE PRECISION TIME0, TIME1 
 58: DOUBLE PRECISION, ALLOCATABLE :: TEMPA(:), TEMPB(:) 
 59: CHARACTER(LEN=5) ZSYMSAVE 
 60: COMMON /SYS/ ZSYMSAVE 
 61: DOUBLE PRECISION XA, XB, YA, YB, ZA, ZB, DMEAN(NATOMS), DA, DB 
 62: INTEGER TRIED(NATOMS), DLIST(NATOMS), SORTLIST(NATOMS), NDUMMY2, INGROUP(NATOMS), NADDED, DOGROUP 
 63:  
 64: !IF (DOGROUP.GT.0) THEN 
 65: !   WRITE(*,'(A,I6,A)') ' lopermdist> Checking only group ',DOGROUP,' and using only QCI active atoms' 
 66: !ENDIF 
 67:  
 68: DBEST=1.0D100 
 69: PERMUTABLE(1:NATOMS)=.FALSE. 
 70: NDUMMY=1 
 71: DO J1=1,NPERMGROUP 
 72:    DO J2=1,NPERMSIZE(J1) 
 73:       PERMUTABLE(PERMGROUP(NDUMMY+J2-1))=.TRUE. 
 74:       INGROUP(PERMGROUP(NDUMMY+J2-1))=J1 
 75:    ENDDO 
 76:    NDUMMY=NDUMMY+NPERMSIZE(J1) 
 77: ENDDO 
 78:  
 79: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS) 
 80: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS) 
 81: ! 
 82: !  Bipartite matching routine for permutations. Coordinates in DUMMYB do not change 
 83: !  but the coordinates in DUMMYA do. DISTANCE is the distance^2 in this case, 
 84: !  and is evaluated as a sum of local distances squared for permutable groups. 
 85: !  We return to label 10 after every round of permutational/orientational alignment 
 86: !  unless we have converged to the identity permutation. 
 87: ! 
 88: !  The maximum number of pair exchanges associated with a group is two. 
 89: !  
 90: DO J1=1,NATOMS 
 91:    NEWPERM(J1)=J1 
 92: ENDDO 
 93: DSUM=0.0D0 
 94: LOCALPERMNEIGH=MIN(LOCALPERMNEIGH,NATOMS) 
 95:  
 96: NDUMMY=1 
 97: DO J1=1,NPERMGROUP 
 98:    IF (DOGROUP.GT.0) THEN 
 99:       IF (J1.GT.DOGROUP) EXIT 
100:       IF (J1.LT.DOGROUP) GOTO 864 ! for QCI test single groups - MUST increment NDUMMY on line 864!! 
101:    ENDIF 
102:    PATOMS=NPERMSIZE(J1) 
103: !  PRINT '(A,I6,A,I6)','group ',J1,' size=',NPERMSIZE(J1) 
104: !  PRINT '(A)','members:' 
105: !  DO J2=1,PATOMS 
106: !     PRINT '(20I6)',NEWPERM(PERMGROUP(NDUMMY+J2-1)) 
107: !  ENDDO 
108:    LDBEST(J1)=1.0D100 
109:    TRIED(1:NATOMS)=0 
110:    DO J2=1,PATOMS 
111:       LPERMBEST(J2)=J2 
112:    ENDDO 
113:    XA=0.0D0; YA=0.0D0; ZA=0.0D0 
114:    XB=0.0D0; YB=0.0D0; ZB=0.0D0 
115:    DMEAN(1:LOCALPERMNEIGH)=1.0D100 
116:    DO J2=1,PATOMS 
117: !     TRIED(NEWPERM(PERMGROUP(NDUMMY+J2-1)))=-1 
118:       TRIED(PERMGROUP(NDUMMY+J2-1))=-1 
119:       PDUMMYA(3*(J2-1)+1)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+1) 
120:       PDUMMYA(3*(J2-1)+2)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+2) 
121:       PDUMMYA(3*(J2-1)+3)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+3) 
122: !     PDUMMYB(3*(J2-1)+1)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+1) 
123: !     PDUMMYB(3*(J2-1)+2)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+2) 
124: !     PDUMMYB(3*(J2-1)+3)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+3) 
125:       PDUMMYB(3*(J2-1)+1)=DUMMYB(3*(PERMGROUP(NDUMMY+J2-1)-1)+1) 
126:       PDUMMYB(3*(J2-1)+2)=DUMMYB(3*(PERMGROUP(NDUMMY+J2-1)-1)+2) 
127:       PDUMMYB(3*(J2-1)+3)=DUMMYB(3*(PERMGROUP(NDUMMY+J2-1)-1)+3) 
128:       XA=XA+PDUMMYA(3*(J2-1)+1) 
129:       YA=YA+PDUMMYA(3*(J2-1)+2) 
130:       ZA=ZA+PDUMMYA(3*(J2-1)+3) 
131:       XB=XB+PDUMMYB(3*(J2-1)+1) 
132:       YB=YB+PDUMMYB(3*(J2-1)+2) 
133:       ZB=ZB+PDUMMYB(3*(J2-1)+3) 
134:    ENDDO 
135:    XA=XA/PATOMS; YA=YA/PATOMS; ZA=ZA/PATOMS 
136:    XB=XB/PATOMS; YB=YB/PATOMS; ZB=ZB/PATOMS 
137:    SPDUMMYA(1:3*PATOMS)=PDUMMYA(1:3*PATOMS) 
138:    SPDUMMYB(1:3*PATOMS)=PDUMMYB(1:3*PATOMS) 
139: ! 
140: ! TRIED(J2) is 0 if atom J2 is eligible to be a neighbour, but has not 
141: ! yet been tried. It is -1 if it is ineligible, or has been tried and 
142: ! broke the alignment. It is +1 if it has been tried and did not break 
143: ! the alignment. It is -1 for atoms already in the set of permutable 
144: ! atoms in question. We add neighbours one at a time in order of  
145: ! increasing distance from primary permutable set 
146: ! and test whether they break the alignment. 
147: ! 
148:    DMEAN(1:NATOMS)=1.0D10 
149:    NDMEAN=0 
150: ! 
151: ! Use active atoms that are constrained to all the atoms in the permutaing group DOGROUP. 
152: ! COMMONCON(J1,1:NCONCOMMON(J1)) are the entries. 
153: ! DMEAN, SORTLIST, TRIED, PERMUTABLE, and DLIST entries refer to original 
154: ! atom labels. Use NEWPERM to find where they are in coordinate lists. 
155: ! 
156:    DO J2=1,NCONCOMMON(J1) 
157:       J3=COMMONCON(J1,J2) 
158:       IF (DOGROUP.GT.0) THEN 
159:          IF (ATOMACTIVE(J3)) THEN 
160:             NDMEAN=NDMEAN+1 
161:             SORTLIST(NDMEAN)=J3 
162:             DMEAN(NDMEAN)=0.0D0  ! don't care about distances in this set up. Use all eligible neighbours. 
163:          ENDIF 
164:       ENDIF 
165:    ENDDO 
166:  
167: !  UTH DJW 
168:  
169: 71 CONTINUE 
170:    PDUMMYA(1:3*PATOMS)=SPDUMMYA(1:3*PATOMS) 
171:    PDUMMYB(1:3*PATOMS)=SPDUMMYB(1:3*PATOMS) 
172:  
173:    LDBESTATOM=1.0D100 
174:    NOTHER=0 
175:    DO J2=1,NATOMS 
176:       IF (TRIED(J2).EQ.1) THEN 
177:          NOTHER=NOTHER+1 
178:          DLIST(NOTHER)=J2 
179:       ENDIF 
180:    ENDDO 
181:    ADDED=.FALSE. 
182:    outer2: DO J2=1,NATOMS 
183: !     WRITE(*,'(A,I6,2G20.10)') 'lopermdist> J2,dmean,localpermcut2=',J2,DMEAN(J2),LOCALPERMCUT2 
184:       IF (DMEAN(J2).GT.LOCALPERMCUT2) THEN 
185: !        PRINT '(A)',' lopermdist> No more atoms within cutoff' 
186:          GOTO 91 
187:       ENDIF 
188:       IF (TRIED(SORTLIST(J2)).EQ.0) THEN 
189:          ADDED=.TRUE. 
190:          NOTHER=NOTHER+1 
191:          IF (NOTHER+PATOMS.GT.NATOMS) THEN 
192:             PRINT '(A,I6)', & 
193:   & ' lopermdist> ERROR *** number of neighbours plus number of permutable atoms exceeds total for group ',J1 
194:             STOP 
195:          ENDIF 
196:          DLIST(NOTHER)=SORTLIST(J2) 
197:          EXIT outer2 
198:       ENDIF 
199:    ENDDO outer2 
200:  
201:    NADDED=1 
202:    IF (PERMUTABLE(DLIST(NOTHER))) THEN 
203: !     IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(NOTHER),' belongs to permutable set ', & 
204: !  &                                INGROUP(DLIST(NOTHER)) 
205:       NDUMMY2=1 
206:       DO J2=1,INGROUP(DLIST(NOTHER))-1 
207:          NDUMMY2=NDUMMY2+NPERMSIZE(J2) 
208:       ENDDO 
209:       DO J2=1,NPERMSIZE(INGROUP(DLIST(NOTHER))) 
210:          IF (PERMGROUP(NDUMMY2+J2-1).EQ.DLIST(NOTHER-NADDED+1)) CYCLE 
211:          IF (TRIED(PERMGROUP(NDUMMY2+J2-1)).EQ.0) THEN 
212:             NOTHER=NOTHER+1 
213:             NADDED=NADDED+1 
214:             IF (NOTHER+PATOMS.GT.NATOMS) THEN 
215:                PRINT '(A,I6)', & 
216:      ' lopermdist> ERROR *** number of neighbours plus number of permutable atoms exceeds total for group ',J1 
217:                STOP 
218:             ENDIF 
219:             DLIST(NOTHER)=PERMGROUP(NDUMMY2+J2-1) 
220: !           IF (DEBUG) PRINT '(A,I6)',' lopermdist> Adding partner atom ',DLIST(NOTHER) 
221:          ELSE 
222:             PRINT '(A,I6,A)',' lopermdist> ERROR *** Partner atom ',DLIST(NOTHER),' has already been tried' 
223:             STOP 
224:          ENDIF 
225:       ENDDO 
226:    ENDIF 
227:     
228:    DO J2=1,NOTHER 
229:       PDUMMYA(3*(PATOMS+J2-1)+1)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+1) 
230:       PDUMMYA(3*(PATOMS+J2-1)+2)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+2) 
231:       PDUMMYA(3*(PATOMS+J2-1)+3)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+3) 
232: !     PDUMMYB(3*(PATOMS+J2-1)+1)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+1) 
233: !     PDUMMYB(3*(PATOMS+J2-1)+2)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+2) 
234: !     PDUMMYB(3*(PATOMS+J2-1)+3)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+3) 
235:       PDUMMYB(3*(PATOMS+J2-1)+1)=DUMMYB(3*(DLIST(J2)-1)+1) 
236:       PDUMMYB(3*(PATOMS+J2-1)+2)=DUMMYB(3*(DLIST(J2)-1)+2) 
237:       PDUMMYB(3*(PATOMS+J2-1)+3)=DUMMYB(3*(DLIST(J2)-1)+3) 
238:    ENDDO 
239: ! 
240: ! Save PDUMMYA and PDUMMYB for cycling over possible orbits in MYORIENT alignment. 
241: ! 
242:    SPDUMMYA(3*PATOMS+1:3*(PATOMS+NOTHER))=PDUMMYA(3*PATOMS+1:3*(PATOMS+NOTHER)) 
243:    SPDUMMYB(3*PATOMS+1:3*(PATOMS+NOTHER))=PDUMMYB(3*PATOMS+1:3*(PATOMS+NOTHER)) 
244:    NCHOOSEB1=0 
245: 66 NCHOOSEB1=NCHOOSEB1+1 
246:    NCHOOSEB2=0 
247: 31 NCHOOSEB2=NCHOOSEB2+1 
248:    NCHOOSE1=0 
249: 65 NCHOOSE1=NCHOOSE1+1 
250:    NCHOOSE2=0 
251: 30 NCHOOSE2=NCHOOSE2+1 
252: ! 
253: ! Reset the coordinates of the PATOMS+NOTHER atoms in PDUMMYA and PDUMMYB 
254: ! to the subset of atoms from COORDSA and COORDSB. 
255: ! 
256:    PDUMMYA(1:3*(PATOMS+NOTHER))=SPDUMMYA(1:3*(PATOMS+NOTHER)) 
257:    PDUMMYB(1:3*(PATOMS+NOTHER))=SPDUMMYB(1:3*(PATOMS+NOTHER)) 
258:  
259:    CALL MYORIENT(PDUMMYA,DUMMY,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,PATOMS+NOTHER,DEBUG,ROTA,ROTINVA,STOCKT) 
260:    PDUMMYA(1:3*(PATOMS+NOTHER))=DUMMY(1:3*(PATOMS+NOTHER)) 
261:    CALL MYORIENT(PDUMMYB,DUMMY,NORBITB1,NCHOOSEB1,NORBITB2,NCHOOSEB2,PATOMS+NOTHER,DEBUG,ROTB,ROTINVB,STOCKT) 
262:    PDUMMYB(1:3*(PATOMS+NOTHER))=DUMMY(1:3*(PATOMS+NOTHER)) 
263:  
264: ! 
265: ! Optimimise permutational isomer for the standard orientation for the 
266: ! current choice of atoms from the possible orbits. 
267: ! 
268: ! MINPERM does not change PDUMMYB and PDUMMYA. 
269: ! 
270: ! Note that LDISTANCE is actually the distance squared. LDBEST also has dimensions of 
271: ! length squared. 
272: ! 
273:    LDISTANCE=0.0D0 
274:    CALL MINPERM(PATOMS+NOTHER, PDUMMYB, PDUMMYA, BOXLX, BOXLY, BOXLZ, BULKT, LPERM, LDISTANCE, DIST2, WORSTRAD)  
275: !  PRINT '(A,I6,A,I6)','for group ',J1,' size=',NPERMSIZE(J1) 
276: !  PRINT '(A)','original and new atom labels:' 
277: !  DO J2=1,PATOMS 
278: !     PRINT '(2I6)',PERMGROUP(NDUMMY+J2-1),NEWPERM(PERMGROUP(NDUMMY+J2-1)) 
279: !  ENDDO 
280:     
281: !  PRINT '(I6,A)',NOTHER,' other atoms:' 
282: !  PRINT '(A)','original and new other atom labels:' 
283: !  DO J2=1,NOTHER 
284: !     PRINT '(2I6)',DLIST(J2),NEWPERM(DLIST(J2)) 
285: !  ENDDO 
286:  
287: !  PRINT '(A,3I6,3G20.10)','J1,PATOMS,NOTHER,root LDBEST(J1),root LDISTANCE=',J1,PATOMS,NOTHER,SQRT(LDBEST(J1)),SQRT(LDISTANCE) 
288: !  PRINT '(A,20I6)','LPERM after MINPERM: ',LPERM(1:PATOMS+NOTHER) 
289:  
290:    LDISTANCE=LDISTANCE 
291:    DO J2=1,PATOMS 
292:       IF (LPERM(J2).GT.PATOMS) THEN 
293:          LDISTANCE=1.0D300 
294: !        IF (DEBUG) PRINT '(A,I6,A,I6,A)',' lopermdist> For group ',J1,' with ',NOTHER,' neighbours - neighbours mix in'  
295: !        IF (DEBUG) PRINT '(A,I6,A,I6)',' lopermdist> atom ',J2,' lperm value is ',LPERM(J2) 
296:          EXIT 
297:       ENDIF 
298:    ENDDO 
299:  
300:  
301:    DO J2=1,NOTHER 
302:       IF (LPERM(PATOMS+J2).NE.PATOMS+J2) THEN 
303: !        IF (DEBUG) PRINT '(A,I6,A,I6)',' lopermdist> Atom ',DLIST(J2),' also needs to permute to ',LPERM(PATOMS+J2) 
304:          IF (PERMUTABLE(DLIST(J2))) THEN 
305: !           IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(J2),' belongs to permutable set ', & 
306: !  &                                INGROUP(DLIST(J2)) 
307:          ELSE 
308: !           IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(J2),' is NOT permutable!' 
309:             LDISTANCE=1.0D300 
310:          ENDIF 
311:       ENDIF 
312:    ENDDO 
313: ! 
314: ! Save the best permutation and local distance for this subset of atoms. 
315: ! NEWPERM and coordinates are only reset after all the cycles over orbits and NEWMINDIST. 
316: ! Hence we need to track a cumulative permutation and save the best current values. 
317: ! 
318:    IF (LDISTANCE.LT.LDBESTATOM) THEN 
319:       LDBESTATOM=LDISTANCE 
320:       LPERMBESTATOM(1:PATOMS)=LPERM(1:PATOMS) 
321:    ENDIF 
322: !  PRINT '(A,2G20.10)','root LDISTANCE,root LDBESTATOM=',SQRT(LDISTANCE),SQRT(LDBESTATOM) 
323:  
324: !  PRINT '(A,4I6,2G20.10)','NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,root LDISTANCE,root LDBEST=', & 
325: ! &                         NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,SQRT(LDISTANCE),SQRT(LDBEST(J1)) 
326:  
327:    IF (NCHOOSE2.LT.NORBIT2) GOTO 30 
328:    IF (NCHOOSE1.LT.NORBIT1) GOTO 65 
329:    IF (NCHOOSEB2.LT.NORBITB2) GOTO 31 
330:    IF (NCHOOSEB1.LT.NORBITB1) GOTO 66 
331:  
332: !  PRINT '(A,2G20.10)','root LDBESTATOM,LOCALPERMCUT=',SQRT(LDBESTATOM),LOCALPERMCUT  
333:    IF (SQRT(LDBESTATOM).GT.LOCALPERMCUT) THEN 
334: !     IF (DEBUG) THEN 
335: !        PRINT '(A,G15.5,A,I6,A,G15.5)',' lopermdist> Best distance ',SQRT(LDBESTATOM), & 
336: ! &                                     ' is too large for atom ',DLIST(NOTHER),' cutoff=',LOCALPERMCUT 
337: !     ENDIF 
338:       TRIED(DLIST(NOTHER))=-1 
339:       IF (NADDED.GT.1) THEN 
340: !        IF (DEBUG) THEN 
341: !           PRINT '(A)',' lopermdist> and partner atoms:' 
342: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1) 
343: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:' 
344: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1) 
345: !        ENDIF 
346:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1 
347:       ENDIF 
348:       GOTO 71 
349:    ELSE 
350: !     IF (DEBUG) PRINT '(A,G20.10,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), & 
351: ! &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours'  
352:       TRIED(DLIST(NOTHER))=1 
353:       IF (NADDED.GT.1) THEN 
354: !        IF (DEBUG) THEN 
355: !           PRINT '(A)',' lopermdist> and partner atoms:' 
356: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1) 
357: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:' 
358: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1) 
359: !        ENDIF 
360:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1 
361:       ENDIF 
362: !     IF (DEBUG) THEN 
363: !        PRINT '(A)',' all other atoms:' 
364: !        PRINT '(20I5)',DLIST(1:NOTHER) 
365: !        IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(DLIST(1:NOTHER)) 
366: !     ENDIF 
367:       LDBEST(J1)=LDBESTATOM 
368:       LPERMBEST(1:PATOMS)=LPERMBESTATOM(1:PATOMS) 
369: !     PRINT '(A,2G20.10)','Updating permutation: sqrt(LDBEST)=',SQRT(LDBEST(J1)) 
370: !     PRINT '(A,10I6)','LPERMBEST: ',LPERMBEST(1:PATOMS) 
371:    ENDIF 
372: ! 
373: ! Add the next eligible atom and try alignment again. 
374: ! Stop if we already have LOCALPERMNEIGH neighbours. 
375: ! 
376:    IF (NOTHER.LT.LOCALPERMNEIGH) GOTO 71 
377:  
378: 91 CONTINUE ! jump here when there are no atoms left to try. 
379:  
380: !  IF (DEBUG) PRINT '(2(A,I6),A,G15.5)',' lopermdist> For group ',J1,' maximum neighbours=', & 
381: ! &                                      NOTHER,' distance=',SQRT(LDBEST(J1)) 
382: ! 
383: ! We now have the best permutation for group J1 and standard orientations 
384: ! based upon all atoms belonging to the two possible orbits that appear 
385: ! for the standard alignment. 
386: ! 
387:    LPERM(1:PATOMS)=LPERMBEST(1:PATOMS) 
388: ! 
389: ! Fill SAVEPERM with NEWPERM, which contains the current best permutation 
390: ! after the previous pass through J1 
391: ! 
392:    SAVEPERM(1:NATOMS)=NEWPERM(1:NATOMS) 
393: ! 
394: ! Update best permutation for atoms in subset J1, specified by PERMGROUP 
395: ! with offset NDUMMY (updated below after each pass through J1) 
396: ! 
397:    DO J2=1,PATOMS 
398:       SAVEPERM(PERMGROUP(NDUMMY+J2-1))=NEWPERM(PERMGROUP(NDUMMY+LPERMBEST(J2)-1)) 
399: !     PRINT '(2(A,I6))',' lopermdist> Atom ',NEWPERM(PERMGROUP(NDUMMY+J2-1)), & 
400: ! &                     ' moves to position ',PERMGROUP(NDUMMY+LPERMBEST(J2)-1) 
401:    ENDDO 
402: !   
403: ! Update permutation of associated atoms, if any. 
404: ! We must do this as we go along, because these atoms could move in more than 
405: ! one permutational group now. 
406: !  
407:    IF (NSETS(J1).GT.0) THEN 
408:       DO J2=1,PATOMS 
409:          DO J3=1,NSETS(J1) 
410:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3) 
411:          ENDDO 
412:       ENDDO 
413:    ENDIF 
414: ! 
415: ! Save current optimal permutation in NEWPERM 
416: ! 
417:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS) 
418:    DSUM=DSUM+SQRT(LDBEST(J1)) 
419: !  PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' best distance=',SQRT(LDBEST(J1)),' total=',DSUM 
420: !  PRINT '(A)','best permutation is now' 
421: !  PRINT '(20I6)',NEWPERM(1:NATOMS) 
422:  
423: ! 
424: ! Update NDUMMY, the cumulative offset for PERMGROUP 
425: ! 
426: 864   NDUMMY=NDUMMY+NPERMSIZE(J1) 
427: ENDDO  !  end of loop over groups of permutable atoms 
428:  
429: IF (DOGROUP.GT.0) THEN 
430:    DISTANCE=SQRT(LDBEST(DOGROUP)) 
431:    NMOVE=0 
432:    DO J2=1,NATOMS 
433:       IF (NEWPERM(J2).NE.J2) THEN 
434: !        WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2 
435:          NMOVE=NMOVE+1 
436:       ENDIF 
437:    ENDDO 
438: !  PRINT '(A,I6,A,I6)',' lopermdist> Total permutations (not applied!) for optimal alignment of group ',DOGROUP,' is ',NMOVE 
439:    RETURN 
440: ENDIF 
441: NMOVE=0 
442: ! DO J2=1,NATOMS 
443: !    IF (NEWPERM(J2).NE.J2) THEN 
444: !       WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2 
445: !       NMOVE=NMOVE+1 
446: !    ENDIF 
447: ! ENDDO 
448:   IF (DEBUG) PRINT '(A,I6)',' lopermdist> Total permutations for optimal alignment (will be applied)=',NMOVE 
449: ! 
450: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA 
451: ! to the current best alignment.  
452: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current 
453: ! NEWPERM. 
454: ! 
455: ! 
456: ! Putting the ALLPERM(J1)=J1 into the second loop causes pgf90 to miscompile!! 
457: ! 
458: DO J1=1,NATOMS 
459:    ALLPERM(J1)=J1 
460: ENDDO 
461: DO J1=1,NATOMS 
462:    SAVEPERM(J1)=ALLPERM(NEWPERM(J1)) 
463: ENDDO 
464: ALLPERM(1:NATOMS)=SAVEPERM(1:NATOMS) 
465: ! 
466: ! At this point DUMMYA should not have changed from COORDSA, so we are 
467: ! putting COORDSA in DUMMY 
468: ! 
469: DUMMY(1:3*NATOMS)=DUMMYA(1:3*NATOMS) 
470: NPERM=0 
471: ! 
472: ! Update coordinates in DUMMYA to current best overall permutation using NEWPERM. 
473: ! We are doing this to operate with NEWPERMDIST in the next block. 
474: ! 
475: DO J3=1,NATOMS 
476:    DUMMYA(3*(J3-1)+1)=DUMMY(3*(NEWPERM(J3)-1)+1) 
477:    DUMMYA(3*(J3-1)+2)=DUMMY(3*(NEWPERM(J3)-1)+2) 
478:    DUMMYA(3*(J3-1)+3)=DUMMY(3*(NEWPERM(J3)-1)+3) 
479:  
480: !  IF (DEBUG) WRITE(*,'(A,I5,A,I5)') ' lopermdist> Overall permutations after MYORIENT alignment:' 
481:    IF (J3.NE.NEWPERM(J3)) THEN 
482: !     IF (DEBUG) WRITE(*,'(A,I5,A,I5)') ' lopermdist> Moving position ',NEWPERM(J3),' to ',J3 
483:       NPERM=NPERM+1 
484:    ENDIF 
485: ENDDO 
486:  
487: DISTANCE=DSUM 
488: ! IF (DEBUG) WRITE(*,'(A,G20.10)') ' lopermdist> After myorient block sum of distances=',DISTANCE 
489: ! 
490: ! Save current best overall distance, permuted version of COORDSA, and permutation. 
491: ! 
492: DBEST=DISTANCE 
493: XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS) 
494: BESTPERM(1:NATOMS)=ALLPERM(1:NATOMS) 
495: ! 
496: ! At this point NEWPERM, ALLPERM, SAVEPERM, BESTPERM 
497: ! are all the same! 
498: ! 
499: ! PRINT '(A)',' lopermdist> NEWPERM, ALLPERM, SAVEPERM, BESTPERM:' 
500: ! PRINT '(4I6)',(NEWPERM(J1),ALLPERM(J1),SAVEPERM(J1),BESTPERM(J1),J1=1,NATOMS) 
501:  
502: !!!!!!!!!!!!!!!!!!!!!!! DEBUG 
503: ! 
504: ! Test distance for COORDSA with permutation applied in BESTPERM 
505: ! 
506: !  DO J1=1,NATOMS 
507: !     DUMMYA(3*(J1-1)+1)=COORDSA(3*(BESTPERM(J1)-1)+1) 
508: !     DUMMYA(3*(J1-1)+2)=COORDSA(3*(BESTPERM(J1)-1)+2) 
509: !     DUMMYA(3*(J1-1)+3)=COORDSA(3*(BESTPERM(J1)-1)+3) 
510: !  ENDDO 
511:  
512: !  CALL NEWMINDIST(COORDSB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT) 
513: !  CALL NEWMINDIST(COORDSB,XBEST,NATOMS,XDUMMY,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT) 
514: !  IF (DEBUG) WRITE(*,'(A,2G20.10)') &  
515: ! &   ' lopermdist> distance check for permuted COORDSA and original COORDSB=',XDUMMY,DISTANCE 
516: !!!!!!!!!!!!!!!!!!!!!!! DEBUG 
517: ! 
518: ! Now align and reorient the permuted coordinates in COORDSA  
519: ! Try using the best locally aligned group of atoms 
520: ! 
521: CALL NEWMINDIST(DUMMYB,XBEST,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT) 
522: IF (DEBUG) PRINT '(A,G20.10)',' lopermdist> after overall alignment distance=',DISTANCE 
523: RMATBEST(1:3,1:3)=RMAT(1:3,1:3) 
524:  
525: COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input! 
526:  
527: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
528: ! IF (DEBUG) PRINT '(A)',' lopermdist> Overall permutation for COORDSA (second argument):' 
529: ! IF (DEBUG) PRINT '(20I6)',BESTPERM(1:NATOMS) 
530:  
531: IF (DEBUG.AND.(DOGROUP.EQ.0)) THEN 
532:    IF (CHRMMT) CALL UPDATENBONDS(COORDSA) 
533:    CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
534:    PRINT '(2(A,G25.15))',' lopermdist> final   energy for structure A=             ',ENERGY,' RMS=',RMS 
535:    IF (ABS(ENERGY-AINIT).GT.EDIFFTOL) THEN 
536:       PRINT '(A)',' minpermdist> WARNING *** energy change for structure A is outside tolerance' 
537:    ENDIF 
538:    IF (CHRMMT) CALL UPDATENBONDS(COORDSB) 
539:    CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
540:    PRINT '(2(A,G25.15))',' lopermdist> final   energy for structure B=             ',ENERGY,' RMS=',RMS 
541:    IF (ABS(ENERGY-BINIT).GT.EDIFFTOL) THEN 
542:       PRINT '(A)',' minpermdist> WARNING *** energy change for structure B is outside tolerance' 
543:    ENDIF 
544: ENDIF 
545:  
546: RETURN 
547: END SUBROUTINE LOPERMDISTCON 
548:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0