hdiff output

r33445/congrad.f90 2017-11-08 10:30:13.301331108 +0000 r33444/congrad.f90 2017-11-08 10:30:14.673349301 +0000
312:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 312:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
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))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))
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)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)
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)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)
316:             ENDIF316:             ENDIF
317:          ENDDO317:          ENDDO
318: !     ENDIF318: !     ENDIF
319:    ENDDO319:    ENDDO
320: ENDIF320: ENDIF
321: 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
 322: 
 323: 
322: !324: !
323: ! Set gradients on frozen atoms to zero.325: ! Set gradients on frozen atoms to zero.
324: !326: !
325: IF (FREEZE) THEN327: IF (FREEZE) THEN
326:    DO J1=2,INTIMAGE+1  328:    DO J1=2,INTIMAGE+1  
327:       DO J2=1,NATOMS329:       DO J2=1,NATOMS
328:          IF (FROZEN(J2)) THEN330:          IF (FROZEN(J2)) THEN
329:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0331:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0
330:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0332:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0
331:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0333:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0
586: !588: !
587: ! This version of congrad tests for an internal minimum in the589: ! This version of congrad tests for an internal minimum in the
588: ! constraint distances as well as the repulsions.590: ! constraint distances as well as the repulsions.
589: !591: !
590: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)592: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
591: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &593: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
592:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &594:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
593:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &595:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &
594:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &596:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &
595:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET, CHECKCONINT, JMAXCON, &597:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET, CHECKCONINT, JMAXCON, &
596:   &            NCONOFF, CONOFFTRIED, INTCONMAX, KINTENDS598:   &            NCONOFF, CONOFFTRIED, INTCONMAX
597: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG599: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
598: IMPLICIT NONE600: IMPLICIT NONE
599:            601:            
600: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2),JMAX,IMAX,JMAXNOFF,IMAXNOFF602: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2),JMAX,IMAX,JMAXNOFF,IMAXNOFF
601: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX, EMAXNOFF603: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX, EMAXNOFF
602: INTEGER JJMAX(INTIMAGE+2)604: INTEGER JJMAX(INTIMAGE+2)
603: DOUBLE PRECISION  EEMAX(INTIMAGE+2)605: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
604: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1606: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
605: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)607: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
606: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI608: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
636: !638: !
637: EMAX=-1.0D200639: EMAX=-1.0D200
638: EMAXNOFF=-1.0D200640: EMAXNOFF=-1.0D200
639: EEMAX(1:INTIMAGE+2)=-1.0D200641: EEMAX(1:INTIMAGE+2)=-1.0D200
640: JJMAX(1:INTIMAGE+2)=-1642: JJMAX(1:INTIMAGE+2)=-1
641: JMAX=-1643: JMAX=-1
642: IMAX=-1644: IMAX=-1
643: JMAXNOFF=-1645: JMAXNOFF=-1
644: IMAXNOFF=-1646: IMAXNOFF=-1
645: DO J2=1,NCONSTRAINT647: DO J2=1,NCONSTRAINT
646: !  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), & 
647: ! &                    ATOMACTIVE(CONI(J2)),ATOMACTIVE(CONJ(J2)) 
648: !  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), & 
649: ! &                    ATOMACTIVE(CONI(J2)),ATOMACTIVE(CONJ(J2)) 
650:    IF (.NOT.CONACTIVE(J2)) CYCLE648:    IF (.NOT.CONACTIVE(J2)) CYCLE
651:    CCLOCAL=CONCUTLOCAL(J2)649:    CCLOCAL=CONCUTLOCAL(J2)
652:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS650:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
653:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)651:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
654: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG652: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG
655:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN653:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN
656:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)654:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)
657:       STOP655:       STOP
658:    ENDIF656:    ENDIF
659: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG657: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG
1079:          DO J2=1,NATOMS1077:          DO J2=1,NATOMS
1080:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 1078:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
1081:                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))1079:                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))
1082:                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)1080:                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)
1083:                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)1081:                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)
1084:             ENDIF1082:             ENDIF
1085:          ENDDO1083:          ENDDO
1086: !     ENDIF1084: !     ENDIF
1087:    ENDDO1085:    ENDDO
1088: ENDIF1086: ENDIF
1089: IF (KINTENDS.GT.0.0D0) THEN 
1090: ! 
1091: ! Extra terms for the two fixed end points. 
1092: ! 
1093:    DO J1=2,INTIMAGE+1 ! sum over images 
1094:       NI1=(3*NATOMS)*(J1-1) 
1095:       NI2=(3*NATOMS)*(INTIMAGE+1) 
1096: ! 
1097: !  Spring between 1 and J1 
1098: ! 
1099:       DPLUS=0.0D0 
1100:       DO J2=1,NATOMS 
1101:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
1102:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(3*(J2-1)+1))**2 & 
1103:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(3*(J2-1)+2))**2 & 
1104:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(3*(J2-1)+3))**2 
1105:          ENDIF 
1106: !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS 
1107:       ENDDO 
1108:       DPLUS=SQRT(DPLUS) 
1109: !     IF (DPLUS.GT.IMSEPMAX) THEN 
1110: !        DUMMY=KINTENDS*0.5D0*(DPLUS-IMSEPMAX)**2 
1111:          DUMMY=KINTENDS*0.5D0*DPLUS**2 
1112:          DUMMY=DUMMY*(1.0D0/(J1-1))**2 ! (INTIMAGE/(J1-1))**2 
1113:          IF (DUMMY.GT.EMAX) THEN 
1114:             IMAX=J1 
1115:             EMAX=DUMMY 
1116:          ENDIF 
1117: !        DUMMY=KINTENDS*0.5D0*DPLUS**2 
1118:          ESPRING=ESPRING+DUMMY 
1119:          IF (DUMMY.GT.EEMAX(J1)) THEN 
1120:             EEMAX(J1)=DUMMY 
1121:          ENDIF 
1122:  
1123: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING 
1124: !        DUMMY=KINTENDS*(DPLUS-IMSEPMAX)/DPLUS 
1125:          DUMMY=KINTENDS 
1126:          DO J2=1,NATOMS 
1127:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
1128:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(3*(J2-1)+1:3*(J2-1)+3)) 
1129:                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) 
1130:             ENDIF 
1131:          ENDDO 
1132: !     ENDIF 
1133: ! 
1134: !  Spring between INTIMAGE+2 and J1 
1135: ! 
1136:       DPLUS=0.0D0 
1137:       DO J2=1,NATOMS 
1138:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
1139:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 & 
1140:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 & 
1141:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2 
1142:          ENDIF 
1143: !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS 
1144:       ENDDO 
1145:       DPLUS=SQRT(DPLUS) 
1146: !     IF (DPLUS.GT.IMSEPMAX) THEN 
1147: !        DUMMY=KINTENDS*0.5D0*(DPLUS-IMSEPMAX)**2 
1148:          DUMMY=KINTENDS*0.5D0*DPLUS**2 
1149:          DUMMY=DUMMY*(INTIMAGE/(INTIMAGE+2-J1))**2 
1150:          IF (DUMMY.GT.EMAX) THEN 
1151:             IMAX=J1 
1152:             EMAX=DUMMY 
1153:          ENDIF 
1154: !        DUMMY=KINTENDS*0.5D0*DPLUS**2 
1155:          ESPRING=ESPRING+DUMMY 
1156:          IF (DUMMY.GT.EEMAX(J1)) THEN 
1157:             EEMAX(J1)=DUMMY 
1158:          ENDIF 
1159:  
1160: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING 
1161: !        DUMMY=KINTENDS*(DPLUS-IMSEPMAX)/DPLUS 
1162:          DUMMY=KINTENDS 
1163:          DO J2=1,NATOMS 
1164:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
1165:                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)) 
1166:                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) 
1167:             ENDIF 
1168:          ENDDO 
1169: !     ENDIF 
1170:    ENDDO 
1171:  
1172: ENDIF 
1173: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest spring  contribution for any image in image ',IMAX1087: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest spring  contribution for any image in image ',IMAX
1174: IF (DEBUG) WRITE(*, '(A,3G20.10)') 'congrad2> ECON,EREP,ESPRING=',ECON,EREP,ESPRING1088:          IF (PRINTE) THEN
 1089:             WRITE(*, '(A,G20.10)') 'ESPRING=',ESPRING
 1090:          ENDIF
1175: !1091: !
1176: ! Set gradients on frozen atoms to zero.1092: ! Set gradients on frozen atoms to zero.
1177: !1093: !
1178: IF (FREEZE) THEN1094: IF (FREEZE) THEN
1179:    DO J1=2,INTIMAGE+1  1095:    DO J1=2,INTIMAGE+1  
1180:       DO J2=1,NATOMS1096:       DO J2=1,NATOMS
1181:          IF (FROZEN(J2)) THEN1097:          IF (FROZEN(J2)) THEN
1182:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D01098:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0
1183:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D01099:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0
1184:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D01100:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0


r33445/intlbfgs.f90 2017-11-08 10:30:13.553334450 +0000 r33444/intlbfgs.f90 2017-11-08 10:30:14.909352431 +0000
 10: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 10: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 11: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 11: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 12: !   GNU General Public License for more details. 12: !   GNU General Public License for more details.
 13: ! 13: !
 14: !   You should have received a copy of the GNU General Public License 14: !   You should have received a copy of the GNU General Public License
 15: !   along with this program; if not, write to the Free Software 15: !   along with this program; if not, write to the Free Software
 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 17: ! 17: !
 18: SUBROUTINE INTLBFGS(QSTART,QFINISH) 18: SUBROUTINE INTLBFGS(QSTART,QFINISH)
 19: USE PORFUNCS 19: USE PORFUNCS
 20: USE KEY, ONLY : FREEZENODEST, FREEZETOL, MAXINTBFGS, CONVR, ATOMSTORES, & 20: USE KEY, ONLY : FREEZENODEST, FREEZETOL, MAXBFGS, CONVR, ATOMSTORES, &
 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, INTMUPDATE, INTDGUESS, & 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, INTMUPDATE, INTDGUESS, &
 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, CONOFFLIST, CONOFFTRIED,  & 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, CONOFFLIST, CONOFFTRIED,  &
 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, & 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, &
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, & 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, &
 37: USE CHIRALITY 37: USE CHIRALITY
 38:  38: 
 39: IMPLICIT NONE  39: IMPLICIT NONE 
 40:  40: 
 41: 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
 42: INTEGER D, U 42: INTEGER D, U
 43: 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)
 44: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ 44: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ
 45: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF 45: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF
 46: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS), STARTGROUP(NPERMGROUP), ENDGROUP(NPERMGROUP) 46: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS), STARTGROUP(NPERMGROUP), ENDGROUP(NPERMGROUP)
 47: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS), IDENTITY, IDONE 47: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS), IDENTITY
 48: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 48: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 49:  49: 
 50: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY 50: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY
 51: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF 51: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF
 52: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2 52: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2
 53: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND 53: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND
 54: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS), ACID, NLASTCHANGE 54: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS), ACID, NLASTCHANGE
 55: LOGICAL CHIRALSR, CHIRALSRP  55: LOGICAL CHIRALSR, CHIRALSRP 
 56: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS) 56: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS)
 57: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, & 57: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, &
 79: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP), CHIRALACTIVE(NATOMS) 79: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP), CHIRALACTIVE(NATOMS)
 80: INTEGER LUNIT, GETUNIT 80: INTEGER LUNIT, GETUNIT
 81: CHARACTER(LEN=2) SDUMMY 81: CHARACTER(LEN=2) SDUMMY
 82: INTEGER JMAXEEE,JMAXRMS 82: INTEGER JMAXEEE,JMAXRMS
 83: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT 83: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT
 84:  84: 
 85: WHOLEDNEB=.FALSE. 85: WHOLEDNEB=.FALSE.
 86: READIMAGET=.FALSE. 86: READIMAGET=.FALSE.
 87: REMOVEIMAGE=.FALSE. 87: REMOVEIMAGE=.FALSE.
 88:  88: 
 89: ! DO J1=1,NATOMS 89: DO J1=1,NATOMS
 90: !    WRITE(*,'(A,2I8)') 'intlbfgs> atom and residue: ',J1,ATOMSTORES(J1) 90:    WRITE(*,'(A,2I8)') 'intlbfgs> atom and residue: ',J1,ATOMSTORES(J1)
 91: ! ENDDO 91: ENDDO
 92: ! STOP 92: STOP
 93: NCONOFF=0 93: NCONOFF=0
 94: AABACK(1:NATOMS)=.FALSE. 94: AABACK(1:NATOMS)=.FALSE.
 95: BACKDONE=.FALSE. 95: BACKDONE=.FALSE.
 96: IF (DOBACK) THEN 96: IF (DOBACK) THEN
 97:    LUNIT=GETUNIT() 97:    LUNIT=GETUNIT()
 98:    OPEN(UNIT=LUNIT,FILE='aabk',STATUS='OLD') 98:    OPEN(UNIT=LUNIT,FILE='aabk',STATUS='OLD')
 99:    DO J1=1,NATOMS 99:    DO J1=1,NATOMS
100:       READ(LUNIT,*,END=861) NDUMMY100:       READ(LUNIT,*,END=861) NDUMMY
101:       AABACK(NDUMMY)=.TRUE.101:       AABACK(NDUMMY)=.TRUE.
102:    ENDDO102:    ENDDO
131:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:INTMUPDATE,(3*NATOMS)*INTIMAGE), &131:   &      DIAG(3*NATOMS*INTIMAGE), STP(3*NATOMS*INTIMAGE), SEARCHSTEP(0:INTMUPDATE,(3*NATOMS)*INTIMAGE), &
132:   &      GDIF(0:INTMUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &132:   &      GDIF(0:INTMUPDATE,(3*NATOMS)*INTIMAGE),GLAST((3*NATOMS)*INTIMAGE), XSAVE((3*NATOMS)*INTIMAGE), &
133:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &133:   &      XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), CHECKG((3*NATOMS)*INTIMAGE), IMGFREEZE(INTIMAGE), &
134:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))134:   &      EEE(INTIMAGE+2), STEPIMAGE(INTIMAGE))
135: 135: 
136: SWITCHED=.FALSE.136: SWITCHED=.FALSE.
137: INTIMAGESAVE=INTIMAGE137: INTIMAGESAVE=INTIMAGE
138: NBACKTRACK=1138: NBACKTRACK=1
139: CALL MYCPU_TIME(STIME,.FALSE.)139: CALL MYCPU_TIME(STIME,.FALSE.)
140: WRITE(*,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1140: WRITE(*,'(A,I6)') ' intlbfgs> Maximum number of steps for constraint potential phase is ',INTSTEPS1
141: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',INTMUPDATE,' maximum step size=',MAXINTBFGS141: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',INTMUPDATE,' maximum step size=',MAXBFGS
142: ADDATOM=.FALSE.142: ADDATOM=.FALSE.
143: NFAIL=0143: NFAIL=0
144: IMGFREEZE(1:INTIMAGE)=.FALSE.144: IMGFREEZE(1:INTIMAGE)=.FALSE.
145: D=(3*NATOMS)*INTIMAGE145: D=(3*NATOMS)*INTIMAGE
146: U=INTMUPDATE146: U=INTMUPDATE
147: NITERDONE=1147: NITERDONE=1
148: NITERUSE=1148: NITERUSE=1
149: NQDONE=0149: NQDONE=0
150: 150: 
151: !151: !
634: 634: 
635: !635: !
636: ! Are we stuck? If so, try resetting problem atoms to previous image.636: ! Are we stuck? If so, try resetting problem atoms to previous image.
637: !637: !
638: IF (QCIRESET) THEN638: IF (QCIRESET) THEN
639: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN639: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN
640: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1640: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1
641:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN ! .AND.(NITERDONE-NLASTCHANGE.GT.QCIRESETINT1)) THEN641:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN ! .AND.(NITERDONE-NLASTCHANGE.GT.QCIRESETINT1)) THEN
642:       CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)642:       CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
643:       CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)643:       CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
644:       WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Dump images. Turn off worst constraint ',JMAXCON, &644:       WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Turn off worst constraint ',JMAXCON,' total off=',NCONOFF+1
645:   &                        ' total off=',NCONOFF+1 
646:       IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN645:       IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN
647:          WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'646:          WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'
648:          WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF647:          WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF
649:          WRITE(*,'(A)') 'CONOFFTRIED:'648:          WRITE(*,'(A)') 'CONOFFTRIED:'
650:          WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT)649:          WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT)
651:          STOP650:          STOP
652:       ENDIF651:       ENDIF
653:       NCONOFF=NCONOFF+1652:       NCONOFF=NCONOFF+1
654:       CONOFFLIST(NCONOFF)=JMAXCON653:       CONOFFLIST(NCONOFF)=JMAXCON
655:       CONACTIVE(JMAXCON)=.FALSE.654:       CONACTIVE(JMAXCON)=.FALSE.
656:       CONOFFTRIED(JMAXCON)=.TRUE.655:       CONOFFTRIED(JMAXCON)=.TRUE.
657:       NLASTGOODE=NITERDONE656:       NLASTGOODE=NITERDONE
658:       LASTGOODE=ETOTAL657:       LASTGOODE=ETOTAL
659: !     IF (NITERDONE.GT.3123) STOP  !!!! DEBUG DJW658:       IF (NITERDONE.GT.1780) STOP  !!!! DEBUG DJW
660: !     STOP 
661:    ENDIF659:    ENDIF
662: ENDIF660: ENDIF
663: 661: 
664: !662: !
665: !  Check permutational alignments. Maintain a list of the permutable groups where all663: !  Check permutational alignments. Maintain a list of the permutable groups where all
666: !  members are active. See if we have any new complete groups. MUST update NDUMMY664: !  members are active. See if we have any new complete groups. MUST update NDUMMY
667: !  counter to step through permutable atom list.665: !  counter to step through permutable atom list.
668: !666: !
669: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN667: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN
670: 668: 
966: 964: 
967: !  We should apply the maximum LBFGS step to each image separately.965: !  We should apply the maximum LBFGS step to each image separately.
968: !  However, using different scale factors for different images leads to huge966: !  However, using different scale factors for different images leads to huge
969: !  discontinuities! Now take the minimum scale factor for all images. DJW 26/11/07967: !  discontinuities! Now take the minimum scale factor for all images. DJW 26/11/07
970: 968: 
971:    STPMIN=1.0D0969:    STPMIN=1.0D0
972:    DO J2=1,INTIMAGE970:    DO J2=1,INTIMAGE
973:       STEPIMAGE(J2) = SQRT(DOT_PRODUCT(SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2), &971:       STEPIMAGE(J2) = SQRT(DOT_PRODUCT(SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2), &
974:   &                                    SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2)))972:   &                                    SEARCHSTEP(POINT,(3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2)))
975:       DUMMY=STEPIMAGE(J2)973:       DUMMY=STEPIMAGE(J2)
976:       IF (STEPIMAGE(J2) > MAXINTBFGS) THEN974:       IF (STEPIMAGE(J2) > MAXBFGS) THEN
977:            STP((3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2) = MAXINTBFGS/STEPIMAGE(J2)975:            STP((3*NATOMS)*(J2-1)+1:(3*NATOMS)*J2) = MAXBFGS/STEPIMAGE(J2)
978:            STPMIN=MIN(STPMIN,STP((3*NATOMS)*(J2-1)+1))976:            STPMIN=MIN(STPMIN,STP((3*NATOMS)*(J2-1)+1))
979:       ENDIF977:       ENDIF
980: !     WRITE(*,'(A,I8,3G20.10)') ' image,initial step size,STP,prod=',J2,DUMMY,STP(3*NATOMS*(J2-1)+1), &978: !     WRITE(*,'(A,I8,3G20.10)') ' image,initial step size,STP,prod=',J2,DUMMY,STP(3*NATOMS*(J2-1)+1), &
981: ! &                                   STEPIMAGE(J2)*STP(3*NATOMS*(J2-1)+1)   979: ! &                                   STEPIMAGE(J2)*STP(3*NATOMS*(J2-1)+1)   
982:    ENDDO980:    ENDDO
983:    STP(1:D)=STPMIN981:    STP(1:D)=STPMIN
984: ! EFK: decide whether to freeze some nodes982: ! EFK: decide whether to freeze some nodes
985:    IF (FREEZENODEST) THEN983:    IF (FREEZENODEST) THEN
986:       TOTGNORM=SQRT(DOT_PRODUCT(G(1:(3*NATOMS)*INTIMAGE),G(1:(3*NATOMS)*INTIMAGE))/INTIMAGE)984:       TOTGNORM=SQRT(DOT_PRODUCT(G(1:(3*NATOMS)*INTIMAGE),G(1:(3*NATOMS)*INTIMAGE))/INTIMAGE)
987:       NIMAGEFREEZE=0985:       NIMAGEFREEZE=0
1573: ! Compute the new step and gradient change1571: ! Compute the new step and gradient change
1574: !1572: !
1575:    NPT=POINT*D1573:    NPT=POINT*D
1576:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)1574:    SEARCHSTEP(POINT,:) = STP*SEARCHSTEP(POINT,:)
1577:    GDIF(POINT,:)=G-GTMP1575:    GDIF(POINT,:)=G-GTMP
1578:    1576:    
1579:    POINT=POINT+1; IF (POINT==INTMUPDATE) POINT=01577:    POINT=POINT+1; IF (POINT==INTMUPDATE) POINT=0
1580: 1578: 
1581:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1579:    IF (DUMPINTXYZ.AND.MOD(NITERDONE,DUMPINTXYZFREQ)==0) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1582:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1580:    IF (DUMPINTEOS.AND.MOD(NITERDONE,DUMPINTEOSFREQ)==0) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1583: !  IF (NITERDONE.GT.2921) CALL INTRWG2(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF) !!! DEBUG DJW1581:    IF (NITERDONE.GT.1780) CALL INTRWG2(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1584: 1582: 
1585:    NITERDONE=NITERDONE+11583:    NITERDONE=NITERDONE+1
1586:    NITERUSE=NITERUSE+11584:    NITERUSE=NITERUSE+1
1587: 1585: 
1588:    IF (NITERDONE.GT.NSTEPSMAX) EXIT1586:    IF (NITERDONE.GT.NSTEPSMAX) EXIT
1589:    IF (NACTIVE.EQ.NATOMS) THEN1587:    IF (NACTIVE.EQ.NATOMS) THEN
1590:       IF (.NOT.SWITCHED) THEN1588:       IF (.NOT.SWITCHED) THEN
1591:          CALL MYCPU_TIME(FTIME,.FALSE.)1589:          CALL MYCPU_TIME(FTIME,.FALSE.)
1592:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1590:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1593:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1591:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1834:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, &1832:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, &
1835:   &             CONOFFTRIED, QCIADDREP, DOBACK, DOBACKALL, QCITRILAT1833:   &             CONOFFTRIED, QCIADDREP, DOBACK, DOBACKALL, QCITRILAT
1836: USE COMMONS, ONLY: NATOMS, DEBUG1834: USE COMMONS, ONLY: NATOMS, DEBUG
1837: IMPLICIT NONE1835: IMPLICIT NONE
1838: INTEGER INTIMAGE1836: INTEGER INTIMAGE
1839: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &1837: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &
1840:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE1838:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE
1841: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN1839: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN
1842: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID1840: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID
1843: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)1841: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)
1844: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID, AABACK(NATOMS), BACKDONE, FTEST, IDONE1842: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID, AABACK(NATOMS), BACKDONE, FTEST
1845: DOUBLE PRECISION P1(3), P2(3), P3(3), SOL1(3), SOL2(3), R1, R2, R31843: DOUBLE PRECISION P1(3), P2(3), P3(3), SOL1(3), SOL2(3), R1, R2, R3
1846: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE01844: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE0
1847: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE1845: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE
1848: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), XSAVED(3,INTIMAGE+2), XSAVEC(3,INTIMAGE+2), XSAVE0(3,INTIMAGE+2),FRAC,RAN1, &1846: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), XSAVED(3,INTIMAGE+2), XSAVEC(3,INTIMAGE+2), XSAVE0(3,INTIMAGE+2),FRAC,RAN1, &
1849:   &              RMS,EEE(INTIMAGE+2),GGG((3*NATOMS)*(INTIMAGE+2)),ETOTAL,DS,DF,DNORM,D1SQ,D2SQ1847:   &              RMS,EEE(INTIMAGE+2),GGG((3*NATOMS)*(INTIMAGE+2)),ETOTAL,DS,DF,DNORM,D1SQ,D2SQ
1850: 1848: 
1851: 1849: 
1852: NTOADD=11850: NTOADD=1
1853: NADDED=01851: NADDED=0
1854: CHOSENACID=.FALSE.1852: CHOSENACID=.FALSE.
2002: !2000: !
2003: !  We need a sorted list of up to 3 active atoms, sorted according to how well the2001: !  We need a sorted list of up to 3 active atoms, sorted according to how well the
2004: !  end point distance is preserved, even if they don't satisfy the constraint 2002: !  end point distance is preserved, even if they don't satisfy the constraint 
2005: !  condition. We want three atoms to use for a local axis system in the interpolation.2003: !  condition. We want three atoms to use for a local axis system in the interpolation.
2006: !2004: !
2007: !  Try sorting on the shortest average distances in the endpoint structures instead, to avoid2005: !  Try sorting on the shortest average distances in the endpoint structures instead, to avoid
2008: !  problems with distant atoms acidentally having a well-preserved distance.2006: !  problems with distant atoms acidentally having a well-preserved distance.
2009: !2007: !
2010:          NDFORNEWATOM=02008:          NDFORNEWATOM=0
2011:          BESTPRESERVEDD(1:NATOMS)=1.0D1002009:          BESTPRESERVEDD(1:NATOMS)=1.0D100
2012: !          DO J1=1,NATOMS2010:          DO J1=1,NATOMS
2013: !             IF (ABS(J1-NEWATOM).GT.INTCONSEP) CYCLE2011:             IF (ABS(J1-NEWATOM).GT.INTCONSEP) CYCLE
2014: !             IF (.NOT.ATOMACTIVE(J1)) CYCLE2012:             IF (.NOT.ATOMACTIVE(J1)) CYCLE
2015: !             DS=SQRT((XYZ(3*(NEWATOM-1)+1)-XYZ(3*(J1-1)+1))**2 &2013:             DS=SQRT((XYZ(3*(NEWATOM-1)+1)-XYZ(3*(J1-1)+1))**2 &
2016: !   &                +(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(J1-1)+2))**2 &2014:   &                +(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(J1-1)+2))**2 &
2017: !   &                +(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(J1-1)+3))**2) 2015:   &                +(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(J1-1)+3))**2) 
2018: !             DF=SQRT((XYZ((INTIMAGE+1)*3*NATOMS+3*(NEWATOM-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+1))**2 &2016:             DF=SQRT((XYZ((INTIMAGE+1)*3*NATOMS+3*(NEWATOM-1)+1)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+1))**2 &
2019: !   &                +(XYZ((INTIMAGE+1)*3*NATOMS+3*(NEWATOM-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+2))**2 &2017:   &                +(XYZ((INTIMAGE+1)*3*NATOMS+3*(NEWATOM-1)+2)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+2))**2 &
2020: !   &                +(XYZ((INTIMAGE+1)*3*NATOMS+3*(NEWATOM-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+3))**2) 2018:   &                +(XYZ((INTIMAGE+1)*3*NATOMS+3*(NEWATOM-1)+3)-XYZ((INTIMAGE+1)*3*NATOMS+3*(J1-1)+3))**2) 
2021: !             IF (DS.GT.INTCONCUT) CYCLE2019:             IF (DS.GT.INTCONCUT) CYCLE
2022: !             IF (DF.GT.INTCONCUT) CYCLE2020:             IF (DF.GT.INTCONCUT) CYCLE
2023: !             DUMMY=ABS(DS-DF)2021:             DUMMY=ABS(DS-DF)
2024: !             NDFORNEWATOM=NDFORNEWATOM+12022:             NDFORNEWATOM=NDFORNEWATOM+1
2025: !             DO J2=1,NDFORNEWATOM 2023:             DO J2=1,NDFORNEWATOM 
2026: !                IF (DUMMY.LT.BESTPRESERVEDD(J2)) THEN2024:                IF (DUMMY.LT.BESTPRESERVEDD(J2)) THEN
2027: ! !                 WRITE(*,'(A,I6,G12.4,I6,G12.4)') 'J1,DUMMY < J2,BESTPRESERVEDD: ',J1,DUMMY,J2,BESTPRESERVEDD(J2)2025: !                 WRITE(*,'(A,I6,G12.4,I6,G12.4)') 'J1,DUMMY < J2,BESTPRESERVEDD: ',J1,DUMMY,J2,BESTPRESERVEDD(J2)
2028: !                   DO J3=NDFORNEWATOM,J2+1,-1 2026:                   DO J3=NDFORNEWATOM,J2+1,-1 
2029: ! !                    WRITE(*,'(A,I6,A,I6,A,G12.4)') ' moving diff and list from ',J3-1,' to ',J3, &2027: !                    WRITE(*,'(A,I6,A,I6,A,G12.4)') ' moving diff and list from ',J3-1,' to ',J3, &
2030: ! !&                                               ' DIFF=',BESTPRESERVEDD(J3-1)2028: !&                                               ' DIFF=',BESTPRESERVEDD(J3-1)
2031: !                      BESTPRESERVEDD(J3)=BESTPRESERVEDD(J3-1)2029:                      BESTPRESERVEDD(J3)=BESTPRESERVEDD(J3-1)
2032: !                      BESTPRESERVEDN(J3)=BESTPRESERVEDN(J3-1)2030:                      BESTPRESERVEDN(J3)=BESTPRESERVEDN(J3-1)
2033: !                   ENDDO2031:                   ENDDO
2034: !                   BESTPRESERVEDD(J2)=DUMMY2032:                   BESTPRESERVEDD(J2)=DUMMY
2035: ! !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTPRESERVEDD element ',J2,' to ',DUMMY2033: !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTPRESERVEDD element ',J2,' to ',DUMMY
2036: !                   BESTPRESERVEDN(J2)=J12034:                   BESTPRESERVEDN(J2)=J1
2037: ! !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTPRESERVEDN element ',J2,' to ',J12035: !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTPRESERVEDN element ',J2,' to ',J1
2038: !                   GOTO 6532036:                   GOTO 653
2039: !                ENDIF2037:                ENDIF
2040: !             ENDDO2038:             ENDDO
2041: ! 653         CONTINUE2039: 653         CONTINUE
2042: !          ENDDO2040:          ENDDO
2043: !          IF (DEBUG) THEN2041:          IF (DEBUG) THEN
2044: !             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> New active atom ',NEWATOM,' best preserved distances:'2042:             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> New active atom ',NEWATOM,' best preserved distances:'
2045: !             WRITE(*,'(20I6)') BESTPRESERVEDN(1:MIN(10,NDFORNEWATOM))2043:             WRITE(*,'(20I6)') BESTPRESERVEDN(1:MIN(10,NDFORNEWATOM))
2046: !             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> sorted differences:'2044:             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> sorted differences:'
2047: !             WRITE(*,'(10G12.4)') BESTPRESERVEDD(1:MIN(10,NDFORNEWATOM))2045:             WRITE(*,'(10G12.4)') BESTPRESERVEDD(1:MIN(10,NDFORNEWATOM))
2048: !          ENDIF2046:          ENDIF
2049:          IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.2047:          IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.
2050: 2048: 
2051:          NCFORNEWATOM=02049:          NCFORNEWATOM=0
2052:          BESTCLOSESTD(1:NATOMS)=1.0D1002050:          BESTCLOSESTD(1:NATOMS)=1.0D100
2053:          DO J1=1,NATOMS2051:          DO J1=1,NATOMS
2054:             IF (ABS(J1-NEWATOM).GT.INTCONSEP) CYCLE2052:             IF (ABS(J1-NEWATOM).GT.INTCONSEP) CYCLE
2055:             IF (.NOT.ATOMACTIVE(J1)) CYCLE2053:             IF (.NOT.ATOMACTIVE(J1)) CYCLE
2056:             DS=SQRT((XYZ(3*(NEWATOM-1)+1)-XYZ(3*(J1-1)+1))**2 &2054:             DS=SQRT((XYZ(3*(NEWATOM-1)+1)-XYZ(3*(J1-1)+1))**2 &
2057:   &                +(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(J1-1)+2))**2 &2055:   &                +(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(J1-1)+2))**2 &
2058:   &                +(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(J1-1)+3))**2) 2056:   &                +(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(J1-1)+3))**2) 
2075:                   BESTCLOSESTD(J2)=DUMMY2073:                   BESTCLOSESTD(J2)=DUMMY
2076: !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTCLOSESTD element ',J2,' to ',DUMMY2074: !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTCLOSESTD element ',J2,' to ',DUMMY
2077:                   BESTCLOSESTN(J2)=J12075:                   BESTCLOSESTN(J2)=J1
2078: !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTCLOSESTN element ',J2,' to ',J12076: !                 WRITE(*,'(A,I6,A,G12.4)') ' setting BESTCLOSESTN element ',J2,' to ',J1
2079:                   GOTO 6592077:                   GOTO 659
2080:                ENDIF2078:                ENDIF
2081:             ENDDO2079:             ENDDO
2082: 659         CONTINUE2080: 659         CONTINUE
2083:          ENDDO2081:          ENDDO
2084:          IF (DEBUG) THEN2082:          IF (DEBUG) THEN
2085:             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> New active atom ',NEWATOM,' closest average distances in endpoints:'2083:             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> New active atom ',NEWATOM,' shortest average distances in endpoints:'
2086:             WRITE(*,'(20I6)') BESTCLOSESTN(1:MIN(10,NCFORNEWATOM))2084:             WRITE(*,'(20I6)') BESTCLOSESTN(1:MIN(10,NCFORNEWATOM))
2087:             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> sorted average distances:'2085:             WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> sorted differences:'
2088:             WRITE(*,'(10G12.4)') BESTCLOSESTD(1:MIN(10,NCFORNEWATOM))2086:             WRITE(*,'(10G12.4)') BESTCLOSESTD(1:MIN(10,NCFORNEWATOM))
2089:          ENDIF2087:          ENDIF
2090: !2088: !
2091: !  Maintain a sorted list of active atoms that are constrained to the new atom, sorted2089: !  Maintain a sorted list of active atoms that are constrained to the new atom, sorted
2092: !  according to their distance.2090: !  according to their distance.
2093: !2091: !
2094:          NCONFORNEWATOM=02092:          NCONFORNEWATOM=0
2095:          CONDIST(1:NATOMS)=1.0D1002093:          CONDIST(1:NATOMS)=1.0D100
2096:          IF (DEBUG) WRITE(*,'(3(A,I6))') ' intlbfgs> New active atom is number ',NEWATOM,' total=',NACTIVE+1, &2094:          IF (DEBUG) WRITE(*,'(3(A,I6))') ' intlbfgs> New active atom is number ',NEWATOM,' total=',NACTIVE+1, &
2097:  &                        ' steps=',NITERDONE2095:  &                        ' steps=',NITERDONE
2257:          ENDIF2255:          ENDIF
2258: 2256: 
2259:          TURNONORDER(NACTIVE)=NEWATOM2257:          TURNONORDER(NACTIVE)=NEWATOM
2260: !2258: !
2261: ! Initial guess for new active atom position. This is crucial for success in INTCONSTRAINT schemes!2259: ! Initial guess for new active atom position. This is crucial for success in INTCONSTRAINT schemes!
2262: !2260: !
2263:          ESAVED=1.0D1002261:          ESAVED=1.0D100
2264:          ESAVE0=1.0D1002262:          ESAVE0=1.0D100
2265:          ESAVEC=1.0D1002263:          ESAVEC=1.0D100
2266:          FTEST=.TRUE.2264:          FTEST=.TRUE.
2267:          IDONE=.FALSE. 
2268:          IF (NCONFORNEWATOM.GE.3) THEN2265:          IF (NCONFORNEWATOM.GE.3) THEN
2269:             IDONE=.TRUE. 
2270: !2266: !
2271: ! Move the new atom consistently in the local environment of its three nearest actively constrained atoms.2267: ! Move the new atom consistently in the local environment of its three nearest actively constrained atoms.
2272: ! Make a local orthogonal coordinate system and use constant components in this basis.2268: ! Make a local orthogonal coordinate system and use constant components in this basis.
2273: !2269: !
2274:             N1=NCONFORNEWATOM-2; N2=NCONFORNEWATOM-1; N3=NCONFORNEWATOM2270:             IF (DEBUG) WRITE(*,'(A,3I6)') ' intlbfgs> initial guess from closest three constrained active atoms, ',CONLIST(1:3)
2275: !           N1=1; N2=2; N3=32271:             VEC1(1:3)=XYZ(3*(CONLIST(2)-1)+1:3*(CONLIST(2)-1)+3)-XYZ(3*(CONLIST(1)-1)+1:3*(CONLIST(1)-1)+3)
2276:             IF (DEBUG) WRITE(*,'(A,3I6)') ' intlbfgs> initial guess from furthest three constrained active atoms, ',CONLIST(N1:N3) 
2277:             VEC1(1:3)=XYZ(3*(CONLIST(N2)-1)+1:3*(CONLIST(N2)-1)+3)-XYZ(3*(CONLIST(N1)-1)+1:3*(CONLIST(N1)-1)+3) 
2278:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)2272:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
2279:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY2273:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
2280:             VEC2(1:3)=XYZ(3*(CONLIST(N3)-1)+1:3*(CONLIST(N3)-1)+3)-XYZ(3*(CONLIST(N1)-1)+1:3*(CONLIST(N1)-1)+3)2274:             VEC2(1:3)=XYZ(3*(CONLIST(3)-1)+1:3*(CONLIST(3)-1)+3)-XYZ(3*(CONLIST(1)-1)+1:3*(CONLIST(1)-1)+3)
2281:             DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)2275:             DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
2282:             VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)2276:             VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
2283:             DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)2277:             DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
2284:             IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2278:             IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2285:             VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2279:             VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2286:             VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2280:             VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2287:             VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2281:             VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2288:             C1=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(N1)-1)+1))*VEC1(1)+ &2282:             C1=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC1(1)+ &
2289:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(N1)-1)+2))*VEC1(2)+ &2283:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC1(2)+ &
2290:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(N1)-1)+3))*VEC1(3)2284:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC1(3)
2291:             C2=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(N1)-1)+1))*VEC2(1)+ &2285:             C2=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC2(1)+ &
2292:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(N1)-1)+2))*VEC2(2)+ &2286:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC2(2)+ &
2293:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(N1)-1)+3))*VEC2(3)2287:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC2(3)
2294:             C3=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(N1)-1)+1))*VEC3(1)+ &2288:             C3=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC3(1)+ &
2295:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(N1)-1)+2))*VEC3(2)+ &2289:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC3(2)+ &
2296:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(N1)-1)+3))*VEC3(3)2290:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC3(3)
2297: 2291: 
2298:             DO J1=2,INTIMAGE+12292:             DO J1=2,INTIMAGE+1
2299:                VEC1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(N2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N2)-1)+3) &2293:                VEC1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3) &
2300:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+3)2294:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)
2301:                DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)2295:                DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
2302:                IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY2296:                IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
2303:                VEC2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(N3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N3)-1)+3) &2297:                VEC2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3) &
2304:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+3)2298:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)
2305:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)2299:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
2306:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)2300:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
2307:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)2301:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
2308:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2302:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2309:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2303:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2310:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2304:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2311:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2305:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2312:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &2306:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
2313:   &               XYZ((J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+3)+C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)2307:   &            XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)+C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)
2314: 2308: 
2315: !2309: !
2316: ! Alternative analytical solution from intersection of three spheres by trilateration2310: ! Alternative analytical solution from intersection of three spheres by trilateration
2317: !2311: !
2318:                IF (QCITRILAT) THEN2312:                IF (QCITRILAT) THEN
2319:                   P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N1)-1)+3)2313:                   P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)
2320:                   P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(N2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N2)-1)+3)2314:                   P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3)
2321:                   P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(N3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(N3)-1)+3)2315:                   P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3)
2322:                   R1=CONDIST(N1)2316:                   R1=CONDIST(1)
2323:                   R2=CONDIST(N2)2317:                   R2=CONDIST(2)
2324:                   R3=CONDIST(N3)2318:                   R3=CONDIST(3)
2325:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)2319:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)
2326:                   IF (FTEST) THEN2320:                   IF (FTEST) THEN
2327: !                    WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J12321: !                    WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1
2328:                   ELSE2322:                   ELSE
2329: !                    WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J12323: !                    WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1
2330: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3)2324: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3)
2331: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3)2325: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3)
2332: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2326: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2333: !                    D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &2327:                      D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
2334: ! &                      +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &2328:   &                      +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
2335: ! &                      +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**22329:   &                      +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
2336: !                    D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &2330:                      D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
2337: ! &                      +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &2331:   &                      +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
2338: ! &                      +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**22332:   &                      +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
2339: ! 
2340: ! Try minmum distance from previous solution 
2341: ! 
2342:                      D1SQ=(SOL1(1)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2343:   &                      +(SOL1(2)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2344:   &                      +(SOL1(3)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2345:                      D2SQ=(SOL2(1)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2346:   &                      +(SOL2(2)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2347:   &                      +(SOL2(3)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2348: !                    WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ2333: !                    WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ
2349:                      IF (D1SQ.LT.D2SQ) THEN2334:                      IF (D1SQ.LT.D2SQ) THEN
2350:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3)2335:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3)
2351:                      ELSE2336:                      ELSE
2352:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)2337:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)
2353:                      ENDIF2338:                      ENDIF
2354:                   ENDIF2339:                   ENDIF
2355:                ENDIF2340:                ENDIF
2356: 2341: 
2357:             ENDDO2342:             ENDDO
2362:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2347:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2363:             ELSE2348:             ELSE
2364:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2349:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2365:             ENDIF2350:             ENDIF
2366:             ESAVE0=ETOTAL2351:             ESAVE0=ETOTAL
2367:             DO J1=2,INTIMAGE+12352:             DO J1=2,INTIMAGE+1
2368:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2353:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2369:             ENDDO2354:             ENDDO
2370:          ENDIF2355:          ENDIF
2371: 2356: 
2372: !          IF ((.NOT.IDONE).AND.(NDFORNEWATOM.GE.3)) THEN2357:          IF (NDFORNEWATOM.GE.3) THEN
2373: !             IDONE=.TRUE.2358: !
2374: ! !2359: ! Choose three atoms from the BESTPRESERVEDN list at random with bias towards the 
2375: ! ! Choose three atoms from the BESTPRESERVEDN list at random with bias towards the 2360: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate
2376: ! ! start of the list. Let the relative weight for position i be 1/i**2 and calculate2361: ! the sum to normalise.
2377: ! ! the sum to normalise.2362: !
2378: ! !2363:             DUMMY=0.0D0
2379: !             DUMMY=0.0D02364:             DO J1=1,NDFORNEWATOM
2380: !             DO J1=1,NDFORNEWATOM2365: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)
2381: ! !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)2366: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTPRESERVEDD(J1))
2382: ! !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTPRESERVEDD(J1))2367:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)
2383: !                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)2368:             ENDDO
2384: !             ENDDO2369:             N1=0; N2=0; N3=0
2385: !             N1=0; N2=0; N3=02370:             DO WHILE (N3.EQ.0)
2386: !             DO WHILE (N3.EQ.0)2371:                DUMMY2=0.0D0
2387: !                DUMMY2=0.0D02372:                RAN1=DPRAND()*DUMMY
2388: !                RAN1=DPRAND()*DUMMY2373:                DO J1=1,NDFORNEWATOM
2389: !                DO J1=1,NDFORNEWATOM2374: !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*J1)
2390: ! !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*J1)2375: !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*BESTPRESERVEDD(J1))
2391: ! !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*BESTPRESERVEDD(J1))2376:                   DUMMY2=DUMMY2+1.0D0/(1.0D0*J1**2)
2392: !                   DUMMY2=DUMMY2+1.0D0/(1.0D0*J1**2)2377:                   IF (DUMMY2.GE.RAN1) THEN
2393: !                   IF (DUMMY2.GE.RAN1) THEN2378:                      IF ((J1.EQ.N1).OR.(J1.EQ.N2)) EXIT ! already chosen
2394: !                      IF ((J1.EQ.N1).OR.(J1.EQ.N2)) EXIT ! already chosen2379:                      IF (N1.EQ.0) THEN
2395: !                      IF (N1.EQ.0) THEN2380:                         N1=J1
2396: !                         N1=J12381:                         EXIT
2397: !                         EXIT2382:                      ENDIF
2398: !                      ENDIF2383:                      IF (N2.EQ.0) THEN
2399: !                      IF (N2.EQ.0) THEN2384:                         N2=J1
2400: !                         N2=J12385:                         EXIT
2401: !                         EXIT2386:                      ENDIF
2402: !                      ENDIF2387:                      N3=J1
2403: !                      N3=J12388:                      EXIT
2404: !                      EXIT2389:                   ENDIF
2405: !                   ENDIF2390:                ENDDO
2406: !                ENDDO2391:             ENDDO
2407: !             ENDDO2392:             IF (DEBUG) WRITE(*,'(A,3I6,A)') ' intlbfgs> choosing positions ',N1,N2,N3,' in best preserved list'
2408: !             IF (DEBUG) WRITE(*,'(A,3I6,A)') ' intlbfgs> choosing positions ',N1,N2,N3,' in best preserved list'2393:             IF (DEBUG) WRITE(*,'(A,3I6)') ' intlbfgs> atoms are ',BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)
2409: !             IF (DEBUG) WRITE(*,'(A,3I6)') ' intlbfgs> atoms are ',BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)2394: !           IF (DEBUG) WRITE(*,'(A,3I6,A)') ' intlbfgs> full list has length ',NDFORNEWATOM
2410: ! !           IF (DEBUG) WRITE(*,'(A,3I6,A)') ' intlbfgs> full list has length ',NDFORNEWATOM2395: !           IF (DEBUG) WRITE(*,'(20I6)') BESTPRESERVEDN(1:NDFORNEWATOM)
2411: ! !           IF (DEBUG) WRITE(*,'(20I6)') BESTPRESERVEDN(1:NDFORNEWATOM) 
2412: !  
2413: ! ! 
2414: ! ! Move the new atom consistently in the local environment of the three active atoms with the 
2415: ! ! best preserved absolute distances or the shortest average distances in the end points. 
2416: ! ! Check the energies and compare linear interpolation as well, then choose the interpolation 
2417: ! ! with the lowest energy. 
2418: ! ! Make a local orthogonal coordinate system and use constant components in this basis. 
2419: ! ! 
2420: !             VEC1(1:3)=XYZ(3*(BESTPRESERVEDN(N2)-1)+1:3*(BESTPRESERVEDN(N2)-1)+3) & 
2421: !   &                  -XYZ(3*(BESTPRESERVEDN(N1)-1)+1:3*(BESTPRESERVEDN(N1)-1)+3) 
2422: !             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2) 
2423: !             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY 
2424: !             VEC2(1:3)=XYZ(3*(BESTPRESERVEDN(N3)-1)+1:3*(BESTPRESERVEDN(N3)-1)+3) & 
2425: !   &                  -XYZ(3*(BESTPRESERVEDN(N1)-1)+1:3*(BESTPRESERVEDN(N1)-1)+3) 
2426: !             DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3) 
2427: !             VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3) 
2428: !             DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2) 
2429: !             IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY 
2430: !             VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2) 
2431: !             VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1) 
2432: !             VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1) 
2433: !             C1=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(BESTPRESERVEDN(N1)-1)+1))*VEC1(1)+ & 
2434: !   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(BESTPRESERVEDN(N1)-1)+2))*VEC1(2)+ & 
2435: !   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(BESTPRESERVEDN(N1)-1)+3))*VEC1(3) 
2436: !             C2=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(BESTPRESERVEDN(N1)-1)+1))*VEC2(1)+ & 
2437: !   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(BESTPRESERVEDN(N1)-1)+2))*VEC2(2)+ & 
2438: !   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(BESTPRESERVEDN(N1)-1)+3))*VEC2(3) 
2439: !             C3=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(BESTPRESERVEDN(N1)-1)+1))*VEC3(1)+ & 
2440: !   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(BESTPRESERVEDN(N1)-1)+2))*VEC3(2)+ & 
2441: !   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(BESTPRESERVEDN(N1)-1)+3))*VEC3(3) 
2442: !             DO J1=2,INTIMAGE+1 
2443: !                VEC1(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+3) & 
2444: !   &                     -XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3) 
2445: !                DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2) 
2446: !                IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY 
2447: !                VEC2(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+3) & 
2448: !   &                     -XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3) 
2449: !                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3) 
2450: !                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3) 
2451: !                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2) 
2452: !                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY 
2453: !                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2) 
2454: !                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1) 
2455: !                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1) 
2456: !                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= & 
2457: !   &            XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)+ & 
2458: !   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)+0.01D0*(DPRAND()-0.5D0)*2.0D0 
2459: ! !              WRITE(*,'(A,I6,3G20.10)') 'intlbfgs> J1,C1,C2,C3=',J1,C1,C2,C3 
2460: ! !              WRITE(*,'(A,9G20.10)') 'intlbfgs> VEC1,2,3=',VEC1(1:3),VEC2(1:3),VEC3(1:3) 
2461: ! !              WRITE(*,'(A,6I6)') 'intlbfgs> N1,N2,N3,Bestpreserved N1,N2,N3=',N1,N2,N3, & 
2462: ! ! &                 BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3) 
2463: !  
2464: ! ! 
2465: ! ! Alternative analytical solution from intersection of three spheres by trilateration not available - we haven't saved the distances, 
2466: ! ! just the differences in the endpoints 
2467: ! ! 
2468: !             ENDDO 
2469: !  
2470: !             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list 
2471: !             IF (QCIADDREP.GT.0) THEN 
2472: !                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2473: !             ELSEIF (CHECKCONINT) THEN 
2474: !                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2475: !             ELSE 
2476: !                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2477: !             ENDIF 
2478: !             ESAVED=ETOTAL 
2479: !             DO J1=2,INTIMAGE+1 
2480: !                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3) 
2481: !             ENDDO 
2482: !          ENDIF 
2483: 2396: 
2484:          IF ((.NOT.IDONE).AND.(NCFORNEWATOM.GE.3)) THEN2397: !
2485:             IDONE=.TRUE.2398: ! Move the new atom consistently in the local environment of the three active atoms with the
 2399: ! best preserved absolute distances or the shortest average distances in the end points.
 2400: ! Check the energies and compare linear interpolation as well, then choose the interpolation
 2401: ! with the lowest energy.
 2402: ! Make a local orthogonal coordinate system and use constant components in this basis.
 2403: !
 2404:             VEC1(1:3)=XYZ(3*(BESTPRESERVEDN(N2)-1)+1:3*(BESTPRESERVEDN(N2)-1)+3) &
 2405:   &                  -XYZ(3*(BESTPRESERVEDN(N1)-1)+1:3*(BESTPRESERVEDN(N1)-1)+3)
 2406:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
 2407:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
 2408:             VEC2(1:3)=XYZ(3*(BESTPRESERVEDN(N3)-1)+1:3*(BESTPRESERVEDN(N3)-1)+3) &
 2409:   &                  -XYZ(3*(BESTPRESERVEDN(N1)-1)+1:3*(BESTPRESERVEDN(N1)-1)+3)
 2410:             DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
 2411:             VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
 2412:             DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
 2413:             IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
 2414:             VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
 2415:             VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
 2416:             VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
 2417:             C1=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(BESTPRESERVEDN(N1)-1)+1))*VEC1(1)+ &
 2418:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(BESTPRESERVEDN(N1)-1)+2))*VEC1(2)+ &
 2419:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(BESTPRESERVEDN(N1)-1)+3))*VEC1(3)
 2420:             C2=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(BESTPRESERVEDN(N1)-1)+1))*VEC2(1)+ &
 2421:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(BESTPRESERVEDN(N1)-1)+2))*VEC2(2)+ &
 2422:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(BESTPRESERVEDN(N1)-1)+3))*VEC2(3)
 2423:             C3=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(BESTPRESERVEDN(N1)-1)+1))*VEC3(1)+ &
 2424:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(BESTPRESERVEDN(N1)-1)+2))*VEC3(2)+ &
 2425:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(BESTPRESERVEDN(N1)-1)+3))*VEC3(3)
 2426:             DO J1=2,INTIMAGE+1
 2427:                VEC1(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+3) &
 2428:   &                     -XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)
 2429:                DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
 2430:                IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
 2431:                VEC2(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+3) &
 2432:   &                     -XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)
 2433:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
 2434:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
 2435:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
 2436:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
 2437:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
 2438:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
 2439:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
 2440:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
 2441:   &            XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)+ &
 2442:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)+0.01D0*(DPRAND()-0.5D0)*2.0D0
 2443: !              WRITE(*,'(A,I6,3G20.10)') 'intlbfgs> J1,C1,C2,C3=',J1,C1,C2,C3
 2444: !              WRITE(*,'(A,9G20.10)') 'intlbfgs> VEC1,2,3=',VEC1(1:3),VEC2(1:3),VEC3(1:3)
 2445: !              WRITE(*,'(A,6I6)') 'intlbfgs> N1,N2,N3,Bestpreserved N1,N2,N3=',N1,N2,N3, &
 2446: ! &                 BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)
 2447: 
 2448: !
 2449: ! Alternative analytical solution from intersection of three spheres by trilateration
 2450: !
 2451:                IF (QCITRILAT) THEN
 2452:                   P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)
 2453:                   P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+3)
 2454:                   P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+3)
 2455:                   R1=BESTPRESERVEDD(N1)
 2456:                   R2=BESTPRESERVEDD(N2)
 2457:                   R3=BESTPRESERVEDD(N3)
 2458:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)
 2459:                   IF (FTEST) THEN
 2460: !                    WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1
 2461:                   ELSE
 2462: !                    WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1
 2463: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3)
 2464: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3)
 2465: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
 2466:                      D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
 2467:      &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
 2468:      &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
 2469:                      D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
 2470:      &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
 2471:      &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
 2472: !                    WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ
 2473:                      IF (D1SQ.LT.D2SQ) THEN
 2474:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3)
 2475:                      ELSE
 2476:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)
 2477:                      ENDIF
 2478:                   ENDIF
 2479:                ENDIF
 2480: 
 2481:             ENDDO
 2482: 
 2483:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
 2484:             IF (QCIADDREP.GT.0) THEN
 2485:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2486:             ELSEIF (CHECKCONINT) THEN
 2487:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2488:             ELSE
 2489:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2490:             ENDIF
 2491:             ESAVED=ETOTAL
 2492:             DO J1=2,INTIMAGE+1
 2493:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
 2494:             ENDDO
 2495:          ENDIF
 2496: 
 2497:          IF (NCFORNEWATOM.GE.3) THEN
2486: !2498: !
2487: ! Choose three atoms from the BESTCLOSEST list at random with bias towards the2499: ! Choose three atoms from the BESTCLOSEST list at random with bias towards the
2488: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate2500: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate
2489: ! the sum to normalise.2501: ! the sum to normalise.
2490: !2502: !
2491:             DUMMY=0.0D02503:             DUMMY=0.0D0
2492: !             DO J1=1,NCFORNEWATOM2504:             DO J1=1,NCFORNEWATOM
2493: ! !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)2505: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)
2494: ! !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTCLOSESTD(J1))2506: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTCLOSESTD(J1))
2495: !                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)2507:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)
2496: !             ENDDO2508:             ENDDO
2497: !             N1=0; N2=0; N3=02509:             N1=0; N2=0; N3=0
2498: !             DO WHILE (N3.EQ.0)2510:             DO WHILE (N3.EQ.0)
2499: !                DUMMY2=0.0D02511:                DUMMY2=0.0D0
2500: !                RAN1=DPRAND()*DUMMY2512:                RAN1=DPRAND()*DUMMY
2501: !                DO J1=1,NCFORNEWATOM2513:                DO J1=1,NCFORNEWATOM
2502: ! !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*J1)2514: !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*J1)
2503: ! !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*BESTCLOSESTD(J1))2515: !                 DUMMY2=DUMMY2+1.0D0/(1.0D0*BESTCLOSESTD(J1))
2504: !                   DUMMY2=DUMMY2+1.0D0/(1.0D0*J1**2)2516:                   DUMMY2=DUMMY2+1.0D0/(1.0D0*J1**2)
2505: !                   IF (DUMMY2.GE.RAN1) THEN2517:                   IF (DUMMY2.GE.RAN1) THEN
2506: !                      IF ((J1.EQ.N1).OR.(J1.EQ.N2)) EXIT ! already chosen2518:                      IF ((J1.EQ.N1).OR.(J1.EQ.N2)) EXIT ! already chosen
2507: !                      IF (N1.EQ.0) THEN2519:                      IF (N1.EQ.0) THEN
2508: !                         N1=J12520:                         N1=J1
2509: !                         EXIT2521:                         EXIT
2510: !                      ENDIF2522:                      ENDIF
2511: !                      IF (N2.EQ.0) THEN2523:                      IF (N2.EQ.0) THEN
2512: !                         N2=J12524:                         N2=J1
2513: !                         EXIT2525:                         EXIT
2514: !                      ENDIF2526:                      ENDIF
2515: !                      N3=J12527:                      N3=J1
2516: !                      EXIT2528:                      EXIT
2517: !                   ENDIF2529:                   ENDIF
2518: !                ENDDO2530:                ENDDO
2519: !             ENDDO2531:             ENDDO
2520:             N1=1; N2=2; N3=3 
2521:             IF (DEBUG) WRITE(*,'(A,3I6,A)') ' intlbfgs> choosing positions ',N1,N2,N3,' in closest list'2532:             IF (DEBUG) WRITE(*,'(A,3I6,A)') ' intlbfgs> choosing positions ',N1,N2,N3,' in closest list'
2522: 2533: 
2523:             VEC1(1:3)=XYZ(3*(BESTCLOSESTN(N2)-1)+1:3*(BESTCLOSESTN(N2)-1)+3)-XYZ(3*(BESTCLOSESTN(N1)-1)+1:3*(BESTCLOSESTN(N1)-1)+3)2534:             VEC1(1:3)=XYZ(3*(BESTCLOSESTN(N2)-1)+1:3*(BESTCLOSESTN(N2)-1)+3)-XYZ(3*(BESTCLOSESTN(N1)-1)+1:3*(BESTCLOSESTN(N1)-1)+3)
2524:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)2535:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
2525:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY2536:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
2526:             VEC2(1:3)=XYZ(3*(BESTCLOSESTN(N3)-1)+1:3*(BESTCLOSESTN(N3)-1)+3)-XYZ(3*(BESTCLOSESTN(N1)-1)+1:3*(BESTCLOSESTN(N1)-1)+3)2537:             VEC2(1:3)=XYZ(3*(BESTCLOSESTN(N3)-1)+1:3*(BESTCLOSESTN(N3)-1)+3)-XYZ(3*(BESTCLOSESTN(N1)-1)+1:3*(BESTCLOSESTN(N1)-1)+3)
2527:             DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)2538:             DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
2528:             VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)2539:             VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
2529:             DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)2540:             DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
2530:             IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2541:             IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2569:                   R2=BESTCLOSESTD(N2)2580:                   R2=BESTCLOSESTD(N2)
2570:                   R3=BESTCLOSESTD(N3)2581:                   R3=BESTCLOSESTD(N3)
2571:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)2582:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)
2572:                   IF (FTEST) THEN2583:                   IF (FTEST) THEN
2573: !                 WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J12584: !                 WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1
2574:                   ELSE2585:                   ELSE
2575: !                 WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J12586: !                 WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1
2576: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3)2587: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3)
2577: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3)2588: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3)
2578: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2589: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2579: !                    D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &2590:                      D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
2580: !    &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &2591:      &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
2581: !    &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**22592:      &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
2582: !                    D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &2593:                      D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
2583: !    &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &2594:      &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
2584: !    &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**22595:      &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
2585: ! 
2586: ! Try minmum distance from previous solution 
2587: ! 
2588:                      D1SQ=(SOL1(1)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2589:   &                      +(SOL1(2)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2590:   &                      +(SOL1(3)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2591:                      D2SQ=(SOL2(1)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2592:   &                      +(SOL2(2)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2593:   &                      +(SOL2(3)-XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2594: !                 WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ2596: !                 WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ
2595:                      IF (D1SQ.LT.D2SQ) THEN2597:                      IF (D1SQ.LT.D2SQ) THEN
2596:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3)2598:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3)
2597:                      ELSE2599:                      ELSE
2598:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)2600:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)
2599:                      ENDIF2601:                      ENDIF
2600:                   ENDIF2602:                   ENDIF
2601:                ENDIF2603:                ENDIF
2602:             ENDDO2604:             ENDDO
2603: 2605: 
2613:             DO J1=2,INTIMAGE+12615:             DO J1=2,INTIMAGE+1
2614:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2616:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2615:             ENDDO2617:             ENDDO
2616:          ENDIF2618:          ENDIF
2617: !2619: !
2618: ! Standard linear interpolation, with constraint distance scaled by FRAC.2620: ! Standard linear interpolation, with constraint distance scaled by FRAC.
2619: ! Works for FRAC as small as 0.1 with repulsion turned off.2621: ! Works for FRAC as small as 0.1 with repulsion turned off.
2620: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements2622: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements
2621: ! in the two end points.2623: ! in the two end points.
2622: !2624: !
2623:          ETOTAL=1.0D1002625:          ETOTAL=1.0D6
2624:          IF (.NOT.IDONE) THEN2626:          FRAC=1.0D0
2625:             FRAC=1.0D02627:          DO J1=2,INTIMAGE+1
2626:             DO J1=2,INTIMAGE+12628:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1)  &
2627:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1)  & 
2628:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &2629:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &
2629:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)2630:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)
2630:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &2631:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &
2631:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &2632:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &
2632:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)2633:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)
2633:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &2634:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &
2634:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &2635:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &
2635:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)2636:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)
2636:             ENDDO2637:          ENDDO
2637:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2638:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2638:             IF (QCIADDREP.GT.0) THEN2639:          IF (QCIADDREP.GT.0) THEN
2639:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2640:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2640:             ELSEIF (CHECKCONINT) THEN2641:          ELSEIF (CHECKCONINT) THEN
2641:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2642:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2642:             ELSE2643:          ELSE
2643:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2644:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2644:             ENDIF 
2645:          ENDIF2645:          ENDIF
2646: 2646: 
2647:          IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &2647:          IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &
2648:   &              ESAVE0,ESAVED,ESAVEC,ETOTAL2648:   &              ESAVE0,ESAVED,ESAVEC,ETOTAL
2649:     2649:     
2650:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN2650:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN
2651:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'2651:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'
2652:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN2652:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN
2653:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'2653:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'
2654:             DO J1=2,INTIMAGE+12654:             DO J1=2,INTIMAGE+1


r33445/key.f90 2017-11-08 10:30:13.781337474 +0000 r33444/key.f90 2017-11-08 10:30:15.137355454 +0000
 93:      &        APP, AMM, APM, XQP,XQM, ALPHAP, ALPHAM, PHIG, SHIFTV, DGUESS, XDGUESS, TIMELIMIT, TSTART, & 93:      &        APP, AMM, APM, XQP,XQM, ALPHAP, ALPHAM, PHIG, SHIFTV, DGUESS, XDGUESS, TIMELIMIT, TSTART, &
 94:      &        INTEPSILON, FRAMEEDIFF, MORPHERISE, MORPHEMAX, MORPHMXSTP, MAXTSENERGY, MAXBARRIER, GCMXSTP, GCSCALE, & 94:      &        INTEPSILON, FRAMEEDIFF, MORPHERISE, MORPHEMAX, MORPHMXSTP, MAXTSENERGY, MAXBARRIER, GCMXSTP, GCSCALE, &
 95:      &        FREEZETOL, DESMAXEJUMP, DESMAXAVGE, STOCKMU, STOCKLAMBDA, EDIFFTOL, GEOMDIFFTOL, DNEBSWITCH, & 95:      &        FREEZETOL, DESMAXEJUMP, DESMAXAVGE, STOCKMU, STOCKLAMBDA, EDIFFTOL, GEOMDIFFTOL, DNEBSWITCH, &
 96:      &        BHACCREJ, BHSTEPSIZE, BHCONV, BHTEMP, BHDISTTHRESH, BHMAXENERGY, BHK, DIJKSTRALOCAL, KADJUSTFRAC, KADJUSTTOL, & 96:      &        BHACCREJ, BHSTEPSIZE, BHCONV, BHTEMP, BHDISTTHRESH, BHMAXENERGY, BHK, DIJKSTRALOCAL, KADJUSTFRAC, KADJUSTTOL, &
 97:      &        BHSFRAC, HMCEIG, HMMXSTP, HMEVMAX, BBRGAM, BBREPS, BBRSIGMA1, BBRSIGMA2, BBRALPHA, BBRCONV, STOCKZTOL, & 97:      &        BHSFRAC, HMCEIG, HMMXSTP, HMEVMAX, BBRGAM, BBREPS, BBRSIGMA1, BBRSIGMA2, BBRALPHA, BBRCONV, STOCKZTOL, &
 98:      &        GBKAPPA, GBKAPPRM, GBMU, GBNU, GBSIGNOT, GBEPSNOT, GBCHI, GBCHIPRM, & 98:      &        GBKAPPA, GBKAPPRM, GBMU, GBNU, GBSIGNOT, GBEPSNOT, GBCHI, GBCHIPRM, &
 99:      &        PYA1(3), PYA2(3), PYSIGNOT, PYEPSNOT, NEBKINITIAL, NEBKFINAL, NEBFACTOR, & 99:      &        PYA1(3), PYA2(3), PYSIGNOT, PYEPSNOT, NEBKINITIAL, NEBKFINAL, NEBFACTOR, &
100:      &        DBEPSBB, DBEPSAB, DBSIGBB, DBSIGAB, DBPMU, YEPS, YKAPPA, MAXMAXBARRIER, PAPALP, PAPEPS, PAPS, PAPCD, &100:      &        DBEPSBB, DBEPSAB, DBSIGBB, DBSIGAB, DBPMU, YEPS, YKAPPA, MAXMAXBARRIER, PAPALP, PAPEPS, PAPS, PAPCD, &
101:      &        RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &101:      &        RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, RHOCH10, RHOC20H, RHOCH20, &
102:      &        EPS11, EPS22, EPS12, MRHO, MRHO11, MRHO22, MRHO12, REQ11, REQ22, REQ12, &102:      &        EPS11, EPS22, EPS12, MRHO, MRHO11, MRHO22, MRHO12, REQ11, REQ22, REQ12, &
103:      &        ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ, KINT, KINTENDS, &103:      &        ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ, KINT, &
104:      &        LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, EVPC, &104:      &        LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, EVPC, &
105:      &        BISECTMAXENERGY, BISECTMINDIST, BLFACTOR, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBRESEEDDEL1, &105:      &        BISECTMAXENERGY, BISECTMINDIST, BLFACTOR, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBRESEEDDEL1, &
106:      &        NEBRESEEDDEL2, INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &106:      &        NEBRESEEDDEL2, INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
107:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &107:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &
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, &


r33445/keywords.f 2017-11-08 10:30:14.041340921 +0000 r33444/keywords.f 2017-11-08 10:30:15.373358583 +0000
698:          CONCUTABST=.TRUE.698:          CONCUTABST=.TRUE.
699:          CONCUTFRAC=0.1D0699:          CONCUTFRAC=0.1D0
700:          CONCUTFRACT=.FALSE.700:          CONCUTFRACT=.FALSE.
701:          CHECKREPINTERVAL=10701:          CHECKREPINTERVAL=10
702:          CHECKREPCUTOFF=2.0702:          CHECKREPCUTOFF=2.0
703:          DUMPINTXYZ=.FALSE.703:          DUMPINTXYZ=.FALSE.
704:          DUMPINTEOS=.FALSE.704:          DUMPINTEOS=.FALSE.
705:          DUMPINTXYZFREQ=100705:          DUMPINTXYZFREQ=100
706:          DUMPINTEOSFREQ=100706:          DUMPINTEOSFREQ=100
707:          KINT=0.0D0707:          KINT=0.0D0
708:          KINTENDS=0.0D0 
709:          QCIAMBERT=.FALSE.708:          QCIAMBERT=.FALSE.
710:          INTMINT=.FALSE.709:          INTMINT=.FALSE.
711:          INTSPRINGACTIVET=.TRUE.710:          INTSPRINGACTIVET=.TRUE.
712:          INTMINFAC=1.0D0711:          INTMINFAC=1.0D0
713:          QCIRADSHIFTT=.FALSE.712:          QCIRADSHIFTT=.FALSE.
714:          QCIRADSHIFT=1.0D0713:          QCIRADSHIFT=1.0D0
715:          QCICYCLEST=.FALSE.714:          QCICYCLEST=.FALSE.
716:          QCICYCDIST=0.0715:          QCICYCDIST=0.0
717:          QCICYCN=100716:          QCICYCN=100
718:          QCIDNEBT=.FALSE.717:          QCIDNEBT=.FALSE.
3823:             KTWNT=.TRUE.3822:             KTWNT=.TRUE.
3824:             IF (NITEMS.GT.1) THEN3823:             IF (NITEMS.GT.1) THEN
3825:                CALL READF(KTWN)3824:                CALL READF(KTWN)
3826:             ENDIF3825:             ENDIF
3827: ! 3826: ! 
3828: ! KINT: force constant for springs in INTCONSTRAINT calculations.3827: ! KINT: force constant for springs in INTCONSTRAINT calculations.
3829: ! Default zero.3828: ! Default zero.
3830: ! 3829: ! 
3831:          ELSE IF (WORD.EQ.'KINT') THEN3830:          ELSE IF (WORD.EQ.'KINT') THEN
3832:             CALL READF(KINT)3831:             CALL READF(KINT)
3833:             IF (NITEMS.GT.2) CALL READF(KINTENDS) 
3834: ! 3832: ! 
3835: ! LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL3833: ! LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL
3836: ! 3834: ! 
3837: ! Use Lanczos to diagonalize the Hamiltonian. Defaults for the three3835: ! Use Lanczos to diagonalize the Hamiltonian. Defaults for the three
3838: ! associated parameters are ACCLAN=1.0D-8 SHIFTLAN=1.0D-2 CUTLAN=-1.0D0.3836: ! associated parameters are ACCLAN=1.0D-8 SHIFTLAN=1.0D-2 CUTLAN=-1.0D0.
3839: ! 3837: ! 
3840:          ELSE IF (WORD.EQ.'LANCZOS') THEN3838:          ELSE IF (WORD.EQ.'LANCZOS') THEN
3841:             LANCZOST=.TRUE.3839:             LANCZOST=.TRUE.
3842:             IF (NITEMS.GT.1) THEN3840:             IF (NITEMS.GT.1) THEN
3843:                CALL READF(ACCLAN)3841:                CALL READF(ACCLAN)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0