hdiff output

r33433/congrad.f90 2017-11-02 15:30:16.805495505 +0000 r33432/congrad.f90 2017-11-02 15:30:18.233514484 +0000
112:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &112:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &
113:  & ' constraint ',JMAX, &113:  & ' constraint ',JMAX, &
114:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF114:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF
115: ENDIF115: ENDIF
116: IF (JMAXNOFF.GT.0) THEN116: IF (JMAXNOFF.GT.0) THEN
117:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &117:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &
118:  & ' constraint ',JMAXNOFF, &118:  & ' constraint ',JMAXNOFF, &
119:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF119:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF
120: ELSEIF (JMAX.GT.0) THEN120: ELSEIF (JMAX.GT.0) THEN
121:    JMAXNOFF=JMAX121:    JMAXNOFF=JMAX
122:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image, setting NCONOFF to 0'122:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image'
123:    CONOFFTRIED(1:INTCONMAX)=.FALSE.123:    CONOFFTRIED(1:INTCONMAX)=.FALSE.
124:    NCONOFF=0 
125: ENDIF124: ENDIF
126: JMAXCON=JMAXNOFF125: JMAXCON=JMAXNOFF
127: 126: 
128: 127: 
129: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes128: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
130: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes129: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
131: 130: 
132: ! INTCONST=INTCONSTRAINREPCUT**13131: ! INTCONST=INTCONSTRAINREPCUT**13
133: 132: 
134: EMAX=-1.0D200133: EMAX=-1.0D200
826:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for any image in image ',IMAX, &825:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for any image in image ',IMAX, &
827:  & ' constraint ',JMAX, &826:  & ' constraint ',JMAX, &
828:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)827:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)
829: ENDIF828: ENDIF
830: IF (JMAXNOFF.GT.0) THEN829: IF (JMAXNOFF.GT.0) THEN
831:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &830:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, &
832:  & ' constraint ',JMAXNOFF, &831:  & ' constraint ',JMAXNOFF, &
833:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF832:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF
834: ELSEIF (JMAX.GT.0) THEN833: ELSEIF (JMAX.GT.0) THEN
835:    JMAXNOFF=JMAX834:    JMAXNOFF=JMAX
836:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> *** Using highest constraint contribution for any image, setting NCONOFF to 0' 
837:    CONOFFTRIED(1:INTCONMAX)=.FALSE.835:    CONOFFTRIED(1:INTCONMAX)=.FALSE.
838:    NCONOFF=0836:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> *** Using highest constraint contribution for any image'
839: ENDIF837: ENDIF
840: JMAXCON=JMAXNOFF838: JMAXCON=JMAXNOFF
841: 839: 
842: ! INTCONST=INTCONSTRAINREPCUT**13840: ! INTCONST=INTCONSTRAINREPCUT**13
843: 841: 
844: EMAX=-1.0D200842: EMAX=-1.0D200
845: EEMAX(1:INTIMAGE+2)=-1.0D200843: EEMAX(1:INTIMAGE+2)=-1.0D200
846: JJMAX(1:INTIMAGE+2)=-1844: JJMAX(1:INTIMAGE+2)=-1
847: JMAX=-1845: JMAX=-1
848: IMAX=-1846: IMAX=-1


r33433/intlbfgs.f90 2017-11-02 15:30:17.061498910 +0000 r33432/intlbfgs.f90 2017-11-02 15:30:18.497517990 +0000
 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 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
 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, &
 58:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), & 58:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), &
 59:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS) 59:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS)
 60: DOUBLE PRECISION, DIMENSION(INTMUPDATE)     :: RHO1,ALPHA 60: DOUBLE PRECISION, DIMENSION(INTMUPDATE)     :: RHO1,ALPHA
 61: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS) 61: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS)
 62: LOGICAL SWITCHED, AABACK(NATOMS), BACKDONE 62: LOGICAL SWITCHED, AABACK(NATOMS), BACKDONE
 63: DOUBLE PRECISION, POINTER :: X(:), G(:) 63: DOUBLE PRECISION, POINTER :: X(:), G(:)
 64: ! 64: !
260: 260: 
261: IF (NATOMS-NQCIFREEZE.LT.INTFREEZEMIN) THEN261: IF (NATOMS-NQCIFREEZE.LT.INTFREEZEMIN) THEN
262:    DO J1=NATOMS,NATOMS-INTFREEZEMIN+1,-1262:    DO J1=NATOMS,NATOMS-INTFREEZEMIN+1,-1
263:       INTFROZEN(DLIST(J1))=.FALSE.263:       INTFROZEN(DLIST(J1))=.FALSE.
264:    ENDDO264:    ENDDO
265:    NQCIFREEZE=MAX(0,NATOMS-INTFREEZEMIN)265:    NQCIFREEZE=MAX(0,NATOMS-INTFREEZEMIN)
266:    WRITE(*,'(A,I6,A)') ' intlbfgs> Freezing ',NQCIFREEZE,' atoms'266:    WRITE(*,'(A,I6,A)') ' intlbfgs> Freezing ',NQCIFREEZE,' atoms'
267: ENDIF267: ENDIF
268: 268: 
269: NLASTGOODE=0269: NLASTGOODE=0
270: NLASTCHANGE=0 
271: LASTGOODE=1.0D100270: LASTGOODE=1.0D100
272: 271: 
273: !272: !
274: ! Constraints are collected in a list and activated via the CONACTIVE(J1)273: ! Constraints are collected in a list and activated via the CONACTIVE(J1)
275: ! logical array. There will generally be of order NATOMS. However, the274: ! logical array. There will generally be of order NATOMS. However, the
276: ! repulsions will scale as NATOMS**2 and are treated differently. The275: ! repulsions will scale as NATOMS**2 and are treated differently. The
277: ! active repulsions are stored sequentially as atoms are added to the276: ! active repulsions are stored sequentially as atoms are added to the
278: ! growing list. This is done even if we have congeom or congeom.dat files277: ! growing list. This is done even if we have congeom or congeom.dat files
279: ! available. In this case we use the fixed list of possible constraints278: ! available. In this case we use the fixed list of possible constraints
280: ! via CHECKPERC, but the list of repulsions and cutoffs is recreated on279: ! via CHECKPERC, but the list of repulsions and cutoffs is recreated on
627: 567 CONTINUE626: 567 CONTINUE
628: 627: 
629: DO ! Main do loop with counter NITERDONE, initially set to one628: DO ! Main do loop with counter NITERDONE, initially set to one
630: 629: 
631: !630: !
632: ! Are we stuck? If so, try resetting problem atoms to previous image.631: ! Are we stuck? If so, try resetting problem atoms to previous image.
633: !632: !
634: IF (QCIRESET) THEN633: IF (QCIRESET) THEN
635: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN634: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN
636: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1635: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1
637:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1).AND.(NITERDONE-NLASTCHANGE.GT.QCIRESETINT1)) THEN636:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN
638:       CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)637:       CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
639:       CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)638:       CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
640:       WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Turn off worst constraint ',JMAXCON,' total off=',NCONOFF+1639:       WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Turn off worst constraint ',JMAXCON,' total off=',NCONOFF+1
641:       IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN640:       IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN
642:          WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'641:          WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'
643:          WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF642:          WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF
644:          WRITE(*,'(A)') 'CONOFFTRIED:'643:          WRITE(*,'(A)') 'CONOFFTRIED:'
645:          WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT)644:          WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT)
646:          STOP645:          STOP
647:       ENDIF646:       ENDIF
648:       NCONOFF=NCONOFF+1647:       NCONOFF=NCONOFF+1
649:       CONOFFLIST(NCONOFF)=JMAXCON648:       CONOFFLIST(NCONOFF)=JMAXCON
650:       CONACTIVE(JMAXCON)=.FALSE.649:       CONACTIVE(JMAXCON)=.FALSE.
651:       CONOFFTRIED(JMAXCON)=.TRUE.650:       CONOFFTRIED(JMAXCON)=.TRUE.
652:       NLASTGOODE=NITERDONE651:       NLASTGOODE=NITERDONE
653:       LASTGOODE=ETOTAL652:       LASTGOODE=ETOTAL
654: !     STOP  !!!! DEBUG DJW653: !     STOP
655:    ENDIF654:    ENDIF
656: ENDIF655: ENDIF
657: 656: 
658: !657: !
659: !  Check permutational alignments. Maintain a list of the permutable groups where all658: !  Check permutational alignments. Maintain a list of the permutable groups where all
660: !  members are active. See if we have any new complete groups. MUST update NDUMMY659: !  members are active. See if we have any new complete groups. MUST update NDUMMY
661: !  counter to step through permutable atom list.660: !  counter to step through permutable atom list.
662: !661: !
663: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN662: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN
664: 663: 
665:    PRINT *,'DOING CHIRALCHECK NOW' 
666:    CHIRALACTIVE(1:NATOMS)=.FALSE.664:    CHIRALACTIVE(1:NATOMS)=.FALSE.
667:    chicheck: DO J1=1,NATOMS665:    chicheck: DO J1=1,NATOMS
668:       IF (.NOT.ATOMACTIVE(J1)) CYCLE chicheck666:       IF (.NOT.ATOMACTIVE(J1)) CYCLE chicheck
669:       IF (NBONDED(J1).NE.4) CYCLE chicheck667:       IF (NBONDED(J1).NE.4) CYCLE chicheck
670:       DO J2=1,4668:       DO J2=1,4
671:          IF (.NOT.ATOMACTIVE(BONDEDLIST(J1,J2))) CYCLE chicheck669:          IF (.NOT.ATOMACTIVE(BONDEDLIST(J1,J2))) CYCLE chicheck
672:       ENDDO670:       ENDDO
673:       IF (J1.EQ.1070) WRITE(*,'(A,I6,A)') 'intlbfgs> Active four-coordinate atom ',J1,' has four active neighbours'671: !     WRITE(*,'(A,I6,A)') 'intlbfgs> Active four-coordinate atom ',J1,' has four active neighbours'
674:       CHIRALACTIVE(J1)=.TRUE.672:       CHIRALACTIVE(J1)=.TRUE.
675:       DO J3=1,INTIMAGE+2673:       DO J3=1,INTIMAGE+2
676:          CENTRE_COORDS(1)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+1)674:          CENTRE_COORDS(1)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+1)
677:          CENTRE_COORDS(2)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+2)675:          CENTRE_COORDS(2)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+2)
678:          CENTRE_COORDS(3)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+3)676:          CENTRE_COORDS(3)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+3)
679: 677: 
680:          DO J4=1,4678:          DO J4=1,4
681:             J2=BONDEDLIST(J1,J4)679:             J2=BONDEDLIST(J1,J4)
682:             NEIGHBOUR_COORDS(3*(J4-1)+1)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+1)680:             NEIGHBOUR_COORDS(3*(J4-1)+1)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+1)
683:             NEIGHBOUR_COORDS(3*(J4-1)+2)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+2)681:             NEIGHBOUR_COORDS(3*(J4-1)+2)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+2)
684:             NEIGHBOUR_COORDS(3*(J4-1)+3)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+3)682:             NEIGHBOUR_COORDS(3*(J4-1)+3)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+3)
685:          ENDDO683:          ENDDO
686:          CHIRALSR=CHIRALITY_SR(NEIGHBOUR_COORDS,CENTRE_COORDS)684:          CHIRALSR=CHIRALITY_SR(NEIGHBOUR_COORDS,CENTRE_COORDS)
687:          IF (J1.EQ.1070) WRITE(*,'(A,I6,A,I6,A,L5)') 'intlbfgs> Atom ',J1,' image ',J3,' chirality=',CHIRALSR685: !        WRITE(*,'(A,I6,A,I6,A,L5)') 'intlbfgs> Atom ',J1,' image ',J3,' chirality=',CHIRALSR
688:          IF ((J3.GT.1).AND.(CHIRALSR.NEQV.CHIRALSRP)) THEN686:          IF ((J3.GT.1).AND.(CHIRALSR.NEQV.CHIRALSRP)) THEN
689:             WRITE(*,'(A,I6,A,I6,A,I6)') 'intlbfgs> Atom ',J1,' image ',J3,' chirality CHANGED; use previous image coordinates'  687:             WRITE(*,'(A,I6,A,I6,A,I6)') 'intlbfgs> Atom ',J1,' image ',J3,' chirality CHANGED; use previous image coordinates'  
690:             NLASTCHANGE=NITERDONE 
691: !688: !
692: ! need to revert to whole aa coordinates, active atoms or not.689: ! need to revert to whole aa coordinates, active atoms or not.
693: !690: !
694:             ACID=ATOMSTORES(J1)691:             ACID=ATOMSTORES(J1)
695: 692: 
696:             DO J4=1,NATOMS693:             DO J4=1,NATOMS
697:                IF (ATOMSTORES(J4).NE.ACID) CYCLE694:                IF (ATOMSTORES(J4).NE.ACID) CYCLE
698:                IF (.NOT.ATOMACTIVE(J4)) CYCLE695:                WRITE(*,'(A,I6,A,I6,A,I6)') 'intlbfgs> Changing atom ',J4,' image ',J3
699:                WRITE(*,'(A,I6,A,I6,A,I6)') 'intlbfgs> Changing active atom ',J4,' image ',J3 
700:                XYZ(3*NATOMS*(J3-1)+3*(J4-1)+1)=XYZ(3*NATOMS*(J3-2)+3*(J4-1)+1)696:                XYZ(3*NATOMS*(J3-1)+3*(J4-1)+1)=XYZ(3*NATOMS*(J3-2)+3*(J4-1)+1)
701:                XYZ(3*NATOMS*(J3-1)+3*(J4-1)+2)=XYZ(3*NATOMS*(J3-2)+3*(J4-1)+2)697:                XYZ(3*NATOMS*(J3-1)+3*(J4-1)+2)=XYZ(3*NATOMS*(J3-2)+3*(J4-1)+2)
702:                XYZ(3*NATOMS*(J3-1)+3*(J4-1)+3)=XYZ(3*NATOMS*(J3-2)+3*(J4-1)+3)698:                XYZ(3*NATOMS*(J3-1)+3*(J4-1)+3)=XYZ(3*NATOMS*(J3-2)+3*(J4-1)+3)
703:             ENDDO699:             ENDDO
704:          ENDIF700:          ENDIF
705:          IF (J3.EQ.1) CHIRALSRP=CHIRALSR  ! just use result for fixed end point image 1701:          IF (J3.EQ.1) CHIRALSRP=CHIRALSR  ! just use result for fixed end point image 1
706:       ENDDO702:       ENDDO
707:    ENDDO chicheck703:    ENDDO chicheck
708: 704: 
709: !  STOP705: !  STOP
771:             DO J2=1,NATOMS767:             DO J2=1,NATOMS
772:                IF (PERMS(J2).NE.J2) THEN768:                IF (PERMS(J2).NE.J2) THEN
773:                   IF (PERMF(J2).EQ.PERMS(J2)) THEN769:                   IF (PERMF(J2).EQ.PERMS(J2)) THEN
774:                      WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &770:                      WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &
775:    &                 ' consistent non-identity lopermdist permutation for start and finish, move atom ',PERMS(J2),' to position ',J2   771:    &                 ' consistent non-identity lopermdist permutation for start and finish, move atom ',PERMS(J2),' to position ',J2   
776: !                WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &772: !                WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &
777: !  &           ' consistent non-identity lopermdist permutation for start and finish, would move atom ',PERMS(J2),' to position ',J2  773: !  &           ' consistent non-identity lopermdist permutation for start and finish, would move atom ',PERMS(J2),' to position ',J2  
778:                      COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+1)774:                      COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+1)
779:                      COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+2)775:                      COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+2)
780:                      COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+3)776:                      COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+3)
781:                      NLASTCHANGE=NITERDONE 
782:                   ELSE777:                   ELSE
783:                      WRITE(*,'(A,I6,A,I6)') ' intlbfgs> inconsistent non-identity lopermdist permutations for start and finish'778:                      WRITE(*,'(A,I6,A,I6)') ' intlbfgs> inconsistent non-identity lopermdist permutations for start and finish'
784:                   ENDIF779:                   ENDIF
785:                ENDIF780:                ENDIF
786:             ENDDO781:             ENDDO
787:             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS)782:             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS)
788:          ENDIF783:          ENDIF
789:        ENDDO784:        ENDDO
790:     ENDDO np785:     ENDDO np
791:     LOCALPERMCUT=SAVELOCALPERMCUT786:     LOCALPERMCUT=SAVELOCALPERMCUT
1516:    EXITSTATUS=01511:    EXITSTATUS=0
1517:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1512:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1518:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 1513:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 
1519:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 1514:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 
1520:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21515:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1521:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1516:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
1522: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX1517: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX
1523: 1518: 
1524:    IF (EXITSTATUS > 0) THEN  1519:    IF (EXITSTATUS > 0) THEN  
1525:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1520:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1526: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771521: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1527:          PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))1522:          PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))
1528:          IF (MAXEEE.GT.MAXCONE*MAX(0.3D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771523:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1529:          IF (NACTIVE.LT.NATOMS) THEN 1524:          IF (NACTIVE.LT.NATOMS) THEN 
1530:             ADDATOM=.TRUE.1525:             ADDATOM=.TRUE.
1531:             GOTO 7771526:             GOTO 777
1532:          ENDIF1527:          ENDIF
1533:          CALL MYCPU_TIME(FTIME,.FALSE.)1528:          CALL MYCPU_TIME(FTIME,.FALSE.)
1534:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1529:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1535:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1530:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1536:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1531:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1537:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1532:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1538:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1533:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
1814: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'1809: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'
1815: 1810: 
1816: END SUBROUTINE WRITEPROFILE1811: END SUBROUTINE WRITEPROFILE
1817: 1812: 
1818: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE,AABACK,BACKDONE)1813: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE,AABACK,BACKDONE)
1819: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &1814: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &
1820:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &1815:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &
1821:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, ATOMSTORES, QCIADDACIDT, &1816:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, ATOMSTORES, QCIADDACIDT, &
1822:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &1817:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &
1823:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, &1818:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, &
1824:   &             CONOFFTRIED, QCIADDREP, DOBACK, DOBACKALL, QCITRILAT1819:   &             CONOFFTRIED, QCIADDREP, DOBACK, DOBACKALL
1825: USE COMMONS, ONLY: NATOMS, DEBUG1820: USE COMMONS, ONLY: NATOMS, DEBUG
1826: IMPLICIT NONE1821: IMPLICIT NONE
1827: INTEGER INTIMAGE1822: INTEGER INTIMAGE
1828: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &1823: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &
1829:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE1824:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE
1830: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN1825: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN
1831: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID1826: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID
1832: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)1827: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)
1833: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID, AABACK(NATOMS), BACKDONE, FTEST1828: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID, AABACK(NATOMS), BACKDONE, FTEST
1834: DOUBLE PRECISION P1(3), P2(3), P3(3), SOL1(3), SOL2(3), R1, R2, R31829: DOUBLE PRECISION P1(3), P2(3), P3(3), SOL1(3), SOL2(3), R1, R2, R3
2293:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2288:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2294:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2289:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2295:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2290:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2296:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2291:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2297:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &2292:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
2298:   &            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)2293:   &            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)
2299: 2294: 
2300: !2295: !
2301: ! Alternative analytical solution from intersection of three spheres by trilateration2296: ! Alternative analytical solution from intersection of three spheres by trilateration
2302: !2297: !
2303:                IF (QCITRILAT) THEN2298:                P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)
2304:                   P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)2299:                P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3)
2305:                   P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3)2300:                P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3)
2306:                   P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3)2301:                R1=CONDIST(1)
2307:                   R1=CONDIST(1)2302:                R2=CONDIST(2)
2308:                   R2=CONDIST(2)2303:                R3=CONDIST(3)
2309:                   R3=CONDIST(3)2304:                CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)
2310:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST)2305:                IF (FTEST) THEN
2311:                   IF (FTEST) THEN2306: !                 WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1
2312: !                    WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J12307:                ELSE
 2308: !                 WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1
 2309: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3)
 2310: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3)
 2311: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
 2312:                   D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
 2313:   &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
 2314:   &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
 2315:                   D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 &
 2316:   &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 &
 2317:   &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2
 2318: !                 WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ
 2319:                   IF (D1SQ.LT.D2SQ) THEN
 2320:                      XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3)
2313:                   ELSE2321:                   ELSE
2314: !                    WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J12322:                      XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3)
2315: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3) 
2316: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3) 
2317: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3) 
2318:                      D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2319:   &                      +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2320:   &                      +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2321:                      D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2322:   &                      +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2323:   &                      +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2324: !                    WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ 
2325:                      IF (D1SQ.LT.D2SQ) THEN 
2326:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3) 
2327:                      ELSE 
2328:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3) 
2329:                      ENDIF 
2330:                   ENDIF2323:                   ENDIF
2331:                ENDIF2324:                ENDIF
2332: 2325: 
2333:             ENDDO2326:             ENDDO
2334:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2327:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2335:             IF (QCIADDREP.GT.0) THEN2328:             IF (QCIADDREP.GT.0) THEN
2336:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2329:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2337:             ELSEIF (CHECKCONINT) THEN2330:             ELSEIF (CHECKCONINT) THEN
2338:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2331:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2339:             ELSE2332:             ELSE
2340:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2333:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2341:             ENDIF2334:             ENDIF
2342:             ESAVE0=ETOTAL2335:             ESAVE0=ETOTAL
2343:             DO J1=2,INTIMAGE+12336:             DO J1=2,INTIMAGE+1
2344:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2337:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2345:             ENDDO2338:             ENDDO
2346:          ENDIF2339:          ENDIF
2347: 2340:          FTEST=.TRUE.  !!!!!!  DEBUG DJW
2348:          IF (NDFORNEWATOM.GE.3) THEN2341:          IF (FTEST.AND.(NDFORNEWATOM.GE.3)) THEN
2349: !2342: !
2350: ! Choose three atoms from the BESTPRESERVEDN list at random with bias towards the 2343: ! Choose three atoms from the BESTPRESERVEDN list at random with bias towards the 
2351: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate2344: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate
2352: ! the sum to normalise.2345: ! the sum to normalise.
2353: !2346: !
2354:             DUMMY=0.0D02347:             DUMMY=0.0D0
2355:             DO J1=1,NDFORNEWATOM2348:             DO J1=1,NDFORNEWATOM
2356: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)2349: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)
2357: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTPRESERVEDD(J1))2350: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTPRESERVEDD(J1))
2358:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)2351:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)
2428:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2421:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2429:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2422:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2430:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2423:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2431:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &2424:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
2432:   &            XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)+ &2425:   &            XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)+ &
2433:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)+0.01D0*(DPRAND()-0.5D0)*2.0D02426:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)+0.01D0*(DPRAND()-0.5D0)*2.0D0
2434: !              WRITE(*,'(A,I6,3G20.10)') 'intlbfgs> J1,C1,C2,C3=',J1,C1,C2,C32427: !              WRITE(*,'(A,I6,3G20.10)') 'intlbfgs> J1,C1,C2,C3=',J1,C1,C2,C3
2435: !              WRITE(*,'(A,9G20.10)') 'intlbfgs> VEC1,2,3=',VEC1(1:3),VEC2(1:3),VEC3(1:3)2428: !              WRITE(*,'(A,9G20.10)') 'intlbfgs> VEC1,2,3=',VEC1(1:3),VEC2(1:3),VEC3(1:3)
2436: !              WRITE(*,'(A,6I6)') 'intlbfgs> N1,N2,N3,Bestpreserved N1,N2,N3=',N1,N2,N3, &2429: !              WRITE(*,'(A,6I6)') 'intlbfgs> N1,N2,N3,Bestpreserved N1,N2,N3=',N1,N2,N3, &
2437: ! &                 BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)2430: ! &                 BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)
2438:  
2439: ! 
2440: ! Alternative analytical solution from intersection of three spheres by trilateration 
2441: ! 
2442:                IF (QCITRILAT) THEN 
2443:                   P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3) 
2444:                   P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N2)-1)+3) 
2445:                   P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N3)-1)+3) 
2446:                   R1=BESTPRESERVEDD(N1) 
2447:                   R2=BESTPRESERVEDD(N2) 
2448:                   R3=BESTPRESERVEDD(N3) 
2449:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST) 
2450:                   IF (FTEST) THEN 
2451: !                    WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1 
2452:                   ELSE 
2453: !                    WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1 
2454: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3) 
2455: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3) 
2456: !                    WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3) 
2457:                      D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2458:      &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2459:      &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2460:                      D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2461:      &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2462:      &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2463: !                    WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ 
2464:                      IF (D1SQ.LT.D2SQ) THEN 
2465:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3) 
2466:                      ELSE 
2467:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3) 
2468:                      ENDIF 
2469:                   ENDIF 
2470:                ENDIF 
2471:  
2472:             ENDDO2431:             ENDDO
2473: 2432: 
2474:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2433:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2475:             IF (QCIADDREP.GT.0) THEN2434:             IF (QCIADDREP.GT.0) THEN
2476:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2435:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2477:             ELSEIF (CHECKCONINT) THEN2436:             ELSEIF (CHECKCONINT) THEN
2478:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2437:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2479:             ELSE2438:             ELSE
2480:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2439:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2481:             ENDIF2440:             ENDIF
2482:             ESAVED=ETOTAL2441:             ESAVED=ETOTAL
2483:             DO J1=2,INTIMAGE+12442:             DO J1=2,INTIMAGE+1
2484:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2443:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2485:             ENDDO2444:             ENDDO
2486:          ENDIF2445:          ENDIF
2487: 2446: 
2488:          IF (NCFORNEWATOM.GE.3) THEN2447:          IF (FTEST.AND.(NCFORNEWATOM.GE.3)) THEN
2489: !2448: !
2490: ! Choose three atoms from the BESTCLOSEST list at random with bias towards the2449: ! Choose three atoms from the BESTCLOSEST list at random with bias towards the
2491: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate2450: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate
2492: ! the sum to normalise.2451: ! the sum to normalise.
2493: !2452: !
2494:             DUMMY=0.0D02453:             DUMMY=0.0D0
2495:             DO J1=1,NCFORNEWATOM2454:             DO J1=1,NCFORNEWATOM
2496: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)2455: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)
2497: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTCLOSESTD(J1))2456: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTCLOSESTD(J1))
2498:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)2457:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)
2552:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)2511:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
2553:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)2512:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
2554:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)2513:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
2555:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2514:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2556:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2515:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2557:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2516:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2558:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2517:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2559:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &2518:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
2560:   &            XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+3)+ &2519:   &            XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+3)+ &
2561:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)2520:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)
2562:  
2563: ! 
2564: ! Alternative analytical solution from intersection of three spheres by trilateration 
2565: ! 
2566:                IF (QCITRILAT) THEN 
2567:                   P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+3) 
2568:                   P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N2)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N2)-1)+3) 
2569:                   P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N3)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N3)-1)+3) 
2570:                   R1=BESTCLOSESTD(N1) 
2571:                   R2=BESTCLOSESTD(N2) 
2572:                   R3=BESTCLOSESTD(N3) 
2573:                   CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST) 
2574:                   IF (FTEST) THEN 
2575: !                 WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1 
2576:                   ELSE 
2577: !                 WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1 
2578: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3) 
2579: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3) 
2580: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3) 
2581:                      D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2582:      &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2583:      &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2584:                      D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2585:      &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2586:      &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2587: !                 WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ 
2588:                      IF (D1SQ.LT.D2SQ) THEN 
2589:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3) 
2590:                      ELSE 
2591:                         XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3) 
2592:                      ENDIF 
2593:                   ENDIF 
2594:                ENDIF 
2595:             ENDDO2521:             ENDDO
2596: 2522: 
2597:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2523:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2598:             IF (QCIADDREP.GT.0) THEN2524:             IF (QCIADDREP.GT.0) THEN
2599:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2525:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2600:             ELSEIF (CHECKCONINT) THEN2526:             ELSEIF (CHECKCONINT) THEN
2601:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2527:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2602:             ELSE2528:             ELSE
2603:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2529:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2604:             ENDIF2530:             ENDIF
2607:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2533:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2608:             ENDDO2534:             ENDDO
2609:          ENDIF2535:          ENDIF
2610: !2536: !
2611: ! Standard linear interpolation, with constraint distance scaled by FRAC.2537: ! Standard linear interpolation, with constraint distance scaled by FRAC.
2612: ! Works for FRAC as small as 0.1 with repulsion turned off.2538: ! Works for FRAC as small as 0.1 with repulsion turned off.
2613: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements2539: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements
2614: ! in the two end points.2540: ! in the two end points.
2615: !2541: !
2616:          ETOTAL=1.0D62542:          ETOTAL=1.0D6
2617:          FRAC=1.0D02543:          IF (FTEST) THEN
2618:          DO J1=2,INTIMAGE+12544:             FRAC=1.0D0
2619:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1)  &2545:             DO J1=2,INTIMAGE+1
 2546:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1)  &
2620:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &2547:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &
2621:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)2548:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)
2622:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &2549:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &
2623:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &2550:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &
2624:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)2551:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)
2625:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &2552:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &
2626:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &2553:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &
2627:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)2554:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)
2628:          ENDDO2555:             ENDDO
2629:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2556:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2630:          IF (QCIADDREP.GT.0) THEN2557:             IF (QCIADDREP.GT.0) THEN
2631:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2558:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2632:          ELSEIF (CHECKCONINT) THEN2559:             ELSEIF (CHECKCONINT) THEN
2633:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2560:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2561:             ELSE
 2562:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 2563:             ENDIF
 2564:             IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &
 2565:   &                 ESAVE0,ESAVED,ESAVEC,ETOTAL
2634:          ELSE2566:          ELSE
2635:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2567:             IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> Using trilateration scheme for interpolation'
2636:          ENDIF2568:          ENDIF
2637:  
2638:          IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', & 
2639:   &              ESAVE0,ESAVED,ESAVEC,ETOTAL 
2640:     2569:     
2641:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN2570:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN
2642:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'2571:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'
2643:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN2572:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN
2644:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'2573:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'
2645:             DO J1=2,INTIMAGE+12574:             DO J1=2,INTIMAGE+1
2646:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVEC(1:3,J1)2575:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVEC(1:3,J1)
2647:             ENDDO2576:             ENDDO
2648:             ETOTAL=ESAVEC2577:             ETOTAL=ESAVEC
2649:          ELSE IF (ESAVED.LT.ESAVE0) THEN2578:          ELSE IF (ESAVED.LT.ESAVE0) THEN


r33433/key.f90 2017-11-02 15:30:17.313502259 +0000 r33432/key.f90 2017-11-02 15:30:18.729521073 +0000
 53:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 53:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 54:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 54:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 55:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 55:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 56:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 56:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 57:      &        PBST, SSHT, GAUSSIAN03, GAUSSIAN09, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 57:      &        PBST, SSHT, GAUSSIAN03, GAUSSIAN09, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 58:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 58:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 59:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 59:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 60:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 60:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 61:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T, & 61:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T, &
 62:      &        MBPOLT, MULTIJOB_MACHINET, DUMPDATA_MACHINET, PLUSSIDET, MINUSSIDET, PUSHOPTT, MLPVB3NNT, GAUSSIAN16, QCICYCLEST, & 62:      &        MBPOLT, MULTIJOB_MACHINET, DUMPDATA_MACHINET, PLUSSIDET, MINUSSIDET, PUSHOPTT, MLPVB3NNT, GAUSSIAN16, QCICYCLEST, &
 63:      &        QCIDNEBT, QCIRESTART, QCILPERMDIST, FASTOVERLAPT, BNB_ALIGNT, QUIPT, QCIADDACIDT, DOBACK, QCIRESET, DOBACKALL, & 63:      &        QCIDNEBT, QCIRESTART, QCILPERMDIST, FASTOVERLAPT, BNB_ALIGNT, QUIPT, QCIADDACIDT, DOBACK, QCIRESET, DOBACKALL
 64:      &        QCITRILAT 
 65:  64: 
 66:  65: 
 67: ! sy349 > for testing the flatpath after dneb 66: ! sy349 > for testing the flatpath after dneb
 68:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 67:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 69:       LOGICAL FLATPATHT 68:       LOGICAL FLATPATHT
 70:  69: 
 71: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 70: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 72:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 71:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 73:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 72:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 74:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 73:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1


r33433/keywords.f 2017-11-02 15:30:17.973511028 +0000 r33432/keywords.f 2017-11-02 15:30:18.993524582 +0000
650:          CONDATT=.FALSE.650:          CONDATT=.FALSE.
651:          QCIPOTT=.FALSE.651:          QCIPOTT=.FALSE.
652:          QCIPOT2T=.FALSE.652:          QCIPOT2T=.FALSE.
653:          QCIADDREP=0653:          QCIADDREP=0
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. 
661:          QCIADDREPCUT=1.0D0660:          QCIADDREPCUT=1.0D0
662:          QCIADDREPEPS=1.0D0661:          QCIADDREPEPS=1.0D0
663:          QCINOREPINT=.FALSE.662:          QCINOREPINT=.FALSE.
664:          MAXNACTIVE=0663:          MAXNACTIVE=0
665: 664: 
666:          FREEZETOL=1.0D-3665:          FREEZETOL=1.0D-3
667:          FLATTESTT=.FALSE.666:          FLATTESTT=.FALSE.
668:          FLATEDIFF=1.0D-6667:          FLATEDIFF=1.0D-6
669:          QCIPERMCHECK=.FALSE.668:          QCIPERMCHECK=.FALSE.
670:          QCIPERMCHECKINT=100669:          QCIPERMCHECKINT=100
3620:                READ(LUNIT,*) CONCUTFIX(1:NCONSTRAINTFIX)3619:                READ(LUNIT,*) CONCUTFIX(1:NCONSTRAINTFIX)
3621:                READ(LUNIT,*) NREPULSIVEFIX3620:                READ(LUNIT,*) NREPULSIVEFIX
3622:                ALLOCATE(REPIFIX(NREPULSIVEFIX),REPJFIX(NREPULSIVEFIX),REPCUTFIX(NREPULSIVEFIX))3621:                ALLOCATE(REPIFIX(NREPULSIVEFIX),REPJFIX(NREPULSIVEFIX),REPCUTFIX(NREPULSIVEFIX))
3623:                READ(LUNIT,*) REPIFIX(1:NREPULSIVEFIX)3622:                READ(LUNIT,*) REPIFIX(1:NREPULSIVEFIX)
3624:                READ(LUNIT,*) REPJFIX(1:NREPULSIVEFIX)3623:                READ(LUNIT,*) REPJFIX(1:NREPULSIVEFIX)
3625:                READ(LUNIT,*) REPCUTFIX(1:NREPULSIVEFIX)3624:                READ(LUNIT,*) REPCUTFIX(1:NREPULSIVEFIX)
3626:                CLOSE(LUNIT)3625:                CLOSE(LUNIT)
3627:                PRINT '(A)',' keyword> Constraint potential parameters read from file congeom.dat'3626:                PRINT '(A)',' keyword> Constraint potential parameters read from file congeom.dat'
3628:                INTCONMAX=NCONSTRAINTFIX3627:                INTCONMAX=NCONSTRAINTFIX
3629:                NREPMAX=NREPULSIVEFIX3628:                NREPMAX=NREPULSIVEFIX
3630:                ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),3629:                ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),CONOFFTRIED(INTCONMAX))
3631:      &                  CONOFFTRIED(INTCONMAX)) 
3632:                CONOFFTRIED(1:INTCONMAX)=.FALSE.3630:                CONOFFTRIED(1:INTCONMAX)=.FALSE.
3633:                ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))3631:                ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
3634:                ALLOCATE(CONACTIVE(NCONSTRAINTFIX))3632:                ALLOCATE(CONACTIVE(NCONSTRAINTFIX))
3635:             ELSE3633:             ELSE
3636:                INQUIRE(FILE='congeom',EXIST=CONFILE)3634:                INQUIRE(FILE='congeom',EXIST=CONFILE)
3637:                NCONGEOM=03635:                NCONGEOM=0
3638:                IF (.NOT.CONFILE) THEN3636:                IF (.NOT.CONFILE) THEN
3639:                   PRINT '(A)',' keyword> WARNING *** no congeom file found. Will use end point minima only.'3637:                   PRINT '(A)',' keyword> WARNING *** no congeom file found. Will use end point minima only.'
3640:                ELSE3638:                ELSE
3641:                   LUNIT=GETUNIT()3639:                   LUNIT=GETUNIT()
3654:                   CLOSE(LUNIT)3652:                   CLOSE(LUNIT)
3655:                   PRINT '(A,I6,A)',' keyword> Read ',NCONGEOM,' reference geometries from congeom file'3653:                   PRINT '(A,I6,A)',' keyword> Read ',NCONGEOM,' reference geometries from congeom file'
3656:                   IF (NCONGEOM.LT.2) PRINT '(A)',' WARNING *** insufficient reference geometries - using end point minima'3654:                   IF (NCONGEOM.LT.2) PRINT '(A)',' WARNING *** insufficient reference geometries - using end point minima'
3657:                ENDIF3655:                ENDIF
3658:             ENDIF3656:             ENDIF
3659: ! 3657: ! 
3660: ! DO NOT Use the quasi-continuous metric for connection attempts, instead of distance.3658: ! DO NOT Use the quasi-continuous metric for connection attempts, instead of distance.
3661: ! 3659: ! 
3662:             INTERPCOSTFUNCTION=.FALSE.3660:             INTERPCOSTFUNCTION=.FALSE.
3663: !3661: !
3664: ! Use trilateration in placing next atom for QCI 
3665: ! 
3666:       ELSE IF (WORD.EQ.'QCITRILAT') THEN 
3667:          QCITRILAT=.TRUE. 
3668: ! 
3669: ! Do the backbone first3662: ! Do the backbone first
3670: !3663: !
3671:       ELSE IF (WORD.EQ.'QCIDOBACK') THEN3664:       ELSE IF (WORD.EQ.'QCIDOBACK') THEN
3672:          DOBACK=.TRUE.3665:          DOBACK=.TRUE.
3673: !3666: !
3674: ! Do the backbone first and add all backbone atoms in each aa in one go3667: ! Do the backbone first and add all backbone atoms in each aa in one go
3675: !3668: !
3676:       ELSE IF (WORD.EQ.'QCIDOBACKALL') THEN3669:       ELSE IF (WORD.EQ.'QCIDOBACKALL') THEN
3677:          DOBACKALL=.TRUE.3670:          DOBACKALL=.TRUE.
3678:          DOBACK=.TRUE.3671:          DOBACK=.TRUE.


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0