hdiff output

r33306/finalq.f 2017-09-13 18:30:47.758440520 +0100 r33305/finalq.f 2017-09-13 18:30:50.454476350 +0100
167:                ENDDO167:                ENDDO
168:                CSMAV(1:3*NATOMS)=CSMAV(1:3*NATOMS)/CSMGPINDEX168:                CSMAV(1:3*NATOMS)=CSMAV(1:3*NATOMS)/CSMGPINDEX
169: !169: !
170: !  Check the CSM for the averaged structure. It should be zero if this structure has the170: !  Check the CSM for the averaged structure. It should be zero if this structure has the
171: !  right point group. Need to reset CSMIMAGES and CSMNORM temporarily. 171: !  right point group. Need to reset CSMIMAGES and CSMNORM temporarily. 
172: !  172: !  
173:                IF (PERMDIST) THEN173:                IF (PERMDIST) THEN
174:                   DO J2=1,CSMGPINDEX174:                   DO J2=1,CSMGPINDEX
175:                      XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)175:                      XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)
176:                      CALL CSMROT(XTEMP,DUMMY,1,J2)176:                      CALL CSMROT(XTEMP,DUMMY,1,J2)
177:                      CALL ALIGN_DECIDE(XTEMP,DUMMY,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY2,DIST2,RIGID,RMAT)177:                      CALL MINPERMDIST(XTEMP,DUMMY,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY2,DIST2,RIGID,RMAT)
178:                      CALL CSMROT(DUMMY,XTEMP,-1,J2) ! need to rotate the permuted rotated images back to the reference orientation178:                      CALL CSMROT(DUMMY,XTEMP,-1,J2) ! need to rotate the permuted rotated images back to the reference orientation
179:                      CSMIMAGES(1+3*NATOMS*(J2-1):3*NATOMS*J2)=XTEMP(1:3*NATOMS)179:                      CSMIMAGES(1+3*NATOMS*(J2-1):3*NATOMS*J2)=XTEMP(1:3*NATOMS)
180:                   ENDDO180:                   ENDDO
181:                ELSE181:                ELSE
182:                   DO J2=1,CSMGPINDEX182:                   DO J2=1,CSMGPINDEX
183:                      CSMIMAGES(1+3*NATOMS*(J2-1):3*NATOMS*J2)=CSMAV(1:3*NATOMS)183:                      CSMIMAGES(1+3*NATOMS*(J2-1):3*NATOMS*J2)=CSMAV(1:3*NATOMS)
184:                   ENDDO184:                   ENDDO
185:                ENDIF185:                ENDIF
186:                AA(1)=0.0D0; AA(2)=0.0D0; AA(3)=6.283185307D0 ! should give an identity matrix186:                AA(1)=0.0D0; AA(2)=0.0D0; AA(3)=6.283185307D0 ! should give an identity matrix
187:                SAVECSMNORM=CSMNORM187:                SAVECSMNORM=CSMNORM


r33306/genrigid.f90 2017-09-13 18:30:47.978443439 +0100 r33305/genrigid.f90 2017-09-13 18:30:50.678479328 +0100
563:      COM = COM / MASS563:      COM = COM / MASS
564:      XRIGIDCOORDS(3*J1-2:3*J1) = COM564:      XRIGIDCOORDS(3*J1-2:3*J1) = COM
565: 565: 
566:      DO J2 = 1, NSITEPERBODY(J1)566:      DO J2 = 1, NSITEPERBODY(J1)
567:         J9 = RIGIDGROUPS(J2, J1)567:         J9 = RIGIDGROUPS(J2, J1)
568:         PP1(3*J2-2:3*J2) = XCOORDS(3*J9-2:3*J9) - COM568:         PP1(3*J2-2:3*J2) = XCOORDS(3*J9-2:3*J9) - COM
569:         PP2(3*J2-2:3*J2) = SITESRIGIDBODY(J2,:,J1)569:         PP2(3*J2-2:3*J2) = SITESRIGIDBODY(J2,:,J1)
570:      ENDDO570:      ENDDO
571:      TEMPPERMDIST = PERMDIST571:      TEMPPERMDIST = PERMDIST
572:      PERMDIST = .FALSE.572:      PERMDIST = .FALSE.
573:      ! This could possibly be replaced with a call to ALIGN_DECIDE instead (might be faster and possibly slightly more 
574:      ! accurate). Need to test this! 
575:      CALL MINPERMDIST(PP1(1:3*NSITEPERBODY(J1)),PP2(1:3*NSITEPERBODY(J1)),NSITEPERBODY(J1),.FALSE., &573:      CALL MINPERMDIST(PP1(1:3*NSITEPERBODY(J1)),PP2(1:3*NSITEPERBODY(J1)),NSITEPERBODY(J1),.FALSE., &
576:           1.0D0,1.0D0,1.0D0,.FALSE.,.FALSE.,D,DIST2,.FALSE.,RMAT)574:           1.0D0,1.0D0,1.0D0,.FALSE.,.FALSE.,D,DIST2,.FALSE.,RMAT)
577:      PERMDIST = TEMPPERMDIST575:      PERMDIST = TEMPPERMDIST
578:      XRIGIDCOORDS(3*NRIGIDBODY+3*J1-2:3*NRIGIDBODY+3*J1) = rot_mx2aa(RMAT)576:      XRIGIDCOORDS(3*NRIGIDBODY+3*J1-2:3*NRIGIDBODY+3*J1) = rot_mx2aa(RMAT)
579: 577: 
580:      IF ( D/NSITEPERBODY(J1) > 0.1D0 ) THEN578:      IF ( D/NSITEPERBODY(J1) > 0.1D0 ) THEN
581:         !print *, 'not going so well...'579:         !print *, 'not going so well...'
582:         WRITE(MYUNIT, '(A, I3)')  'Warning: Genrigid > mapping looks bad for RB no ', J1 580:         WRITE(MYUNIT, '(A, I3)')  'Warning: Genrigid > mapping looks bad for RB no ', J1 
583:         WRITE(MYUNIT, '(A)')  'Warning: Genrigid >  Often it is the permutation of the RB members, e.g. Hs in NH3'581:         WRITE(MYUNIT, '(A)')  'Warning: Genrigid >  Often it is the permutation of the RB members, e.g. Hs in NH3'
584:      ENDIF582:      ENDIF


r33306/intlbfgs.f90 2017-09-13 18:30:48.210446523 +0100 r33305/intlbfgs.f90 2017-09-13 18:30:50.898482249 +0100
900: !        ENDIF900: !        ENDIF
901: !        WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT901: !        WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT
902:       ENDIF902:       ENDIF
903:    ENDIF903:    ENDIF
904: !904: !
905: ! End of add/subtract images block.905: ! End of add/subtract images block.
906: !906: !
907:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN907:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN
908:       LDEBUG=.FALSE.908:       LDEBUG=.FALSE.
909:       DO J2=2,INTIMAGE+2909:       DO J2=2,INTIMAGE+2
910:          CALL ALIGN_DECIDE(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &910:          CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
911:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DIST,DIST2,RIGID,RMAT)911:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DIST,DIST2,RIGID,RMAT)
912:       ENDDO912:       ENDDO
913:    ENDIF913:    ENDIF
914: 914: 
915:    IF (.NOT.SWITCHED) THEN915:    IF (.NOT.SWITCHED) THEN
916:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)916:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
917:       IF (QCIADDREP.GT.0) THEN917:       IF (QCIADDREP.GT.0) THEN
918:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)918:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
919:       ELSEIF (CHECKCONINT) THEN919:       ELSEIF (CHECKCONINT) THEN
920:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)920:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2229:       ENDDO2229:       ENDDO
2230:       NCONSTRAINT=J22230:       NCONSTRAINT=J2
2231:       WRITE(MYUNIT,'(A,I6,A)') ' checkperc> After allowing for frozen atoms there are ',NCONSTRAINT,' constraints'2231:       WRITE(MYUNIT,'(A,I6,A)') ' checkperc> After allowing for frozen atoms there are ',NCONSTRAINT,' constraints'
2232:       RETURN 2232:       RETURN 
2233:    ELSE2233:    ELSE
2234: !2234: !
2235: ! Put reference minima in optimal permutational alignment with reference minimum one.2235: ! Put reference minima in optimal permutational alignment with reference minimum one.
2236: !2236: !
2237:       DO J2=2,NCONGEOM2237:       DO J2=2,NCONGEOM
2238:          LDEBUG=.FALSE.2238:          LDEBUG=.FALSE.
2239:          CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),CONGEOM(J2,1:3*NATOMS),NATOMS,LDEBUG, &2239:          CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),CONGEOM(J2,1:3*NATOMS),NATOMS,LDEBUG, &
2240:   &                       BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,DIST2,RIGID,RMAT)2240:   &                       BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,DIST2,RIGID,RMAT)
2241:       ENDDO2241:       ENDDO
2242:    ENDIF2242:    ENDIF
2243:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))2243:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))
2244: ENDIF2244: ENDIF
2245: 2245: 
2246: INQUIRE(FILE='constraintfile',EXIST=CONFILET)2246: INQUIRE(FILE='constraintfile',EXIST=CONFILET)
2247: 2247: 
2248: 51   NCONSTRAINT=0 2248: 51   NCONSTRAINT=0 
2249: MAXCONDIST=-1.0D02249: MAXCONDIST=-1.0D0


r33306/io1.f 2017-09-13 18:30:48.438449553 +0100 r33305/io1.f 2017-09-13 18:30:51.178485972 +0100
 72: !         CLOSE(7) 72: !         CLOSE(7)
 73: !         OPEN(UNIT=7,FILE='coords',STATUS='OLD') 73: !         OPEN(UNIT=7,FILE='coords',STATUS='OLD')
 74: !js850>  this step (NATOMS=NATOMS/NPAR) is moved to earlier in the program (to 74: !js850>  this step (NATOMS=NATOMS/NPAR) is moved to earlier in the program (to
 75: !COUNTATOMS), so that the allocation is done with the correct NATOMS 75: !COUNTATOMS), so that the allocation is done with the correct NATOMS
 76: !        IF (MOD(NATOMS,NPAR).NE.0) THEN 76: !        IF (MOD(NATOMS,NPAR).NE.0) THEN
 77: !           WRITE(MYUNIT,'(A,I8,A,I8)') 'Number of atoms in coords file=',NATOMS, 77: !           WRITE(MYUNIT,'(A,I8,A,I8)') 'Number of atoms in coords file=',NATOMS,
 78: !    &               ' is not divisible by number of runs=',NPAR 78: !    &               ' is not divisible by number of runs=',NPAR
 79: !           STOP 79: !           STOP
 80: !        ENDIF 80: !        ENDIF
 81: !        NATOMS=NATOMS/NPAR 81: !        NATOMS=NATOMS/NPAR
 82:          IF (PERMOPT.OR.PERMINVOPT.OR.ALIGNT) THEN 82:          IF (PERMOPT.OR.PERMINVOPT) THEN
 83:             OPEN(UNIT=17,FILE='finish',STATUS='OLD') 83:             OPEN(UNIT=17,FILE='finish',STATUS='OLD')
 84:             READ(17,*) (FIN(J1),J1=1,3*NATOMS) 84:             READ(17,*) (FIN(J1),J1=1,3*NATOMS)
 85:             WRITE(MYUNIT,'(A)') 'Target coordinates read from file finish' 85:             WRITE(MYUNIT,'(A)') 'Target coordinates read from file finish'
 86:          ENDIF 86:          ENDIF
 87:          IF (TSALLIST.AND.(QTSALLIS.EQ.0)) QTSALLIS=1.0D0+1.0D0/(3*NATOMS) 87:          IF (TSALLIST.AND.(QTSALLIS.EQ.0)) QTSALLIS=1.0D0+1.0D0/(3*NATOMS)
 88:          IF (DFTBT.OR.TOSI.OR.WELCH) THEN 88:          IF (DFTBT.OR.TOSI.OR.WELCH) THEN
 89: !           REWIND(7) ! this seems to be needed now? 89: !           REWIND(7) ! this seems to be needed now?
 90:             DO JP=1,NPAR 90:             DO JP=1,NPAR
 91:                DO J1=1,NATOMS 91:                DO J1=1,NATOMS
 92:                   CALL INPUT(END, COORDS_UNIT) 92:                   CALL INPUT(END, COORDS_UNIT)


r33306/main.F 2017-09-13 18:30:48.654452425 +0100 r33305/main.F 2017-09-13 18:30:51.438489429 +0100
556:             COORDS(1:3*NATOMSALLOC,1)=ENDCOORDS(1,1:3*NATOMSALLOC)556:             COORDS(1:3*NATOMSALLOC,1)=ENDCOORDS(1,1:3*NATOMSALLOC)
557:             FINISH(1:3*NATOMSALLOC)=ENDCOORDS(2,1:3*NATOMSALLOC)557:             FINISH(1:3*NATOMSALLOC)=ENDCOORDS(2,1:3*NATOMSALLOC)
558:             DEALLOCATE(ENDCOORDS)558:             DEALLOCATE(ENDCOORDS)
559:          ELSE559:          ELSE
560: !560: !
561: ! align the two endpoints 561: ! align the two endpoints 
562: !562: !
563:             ALLOCATE(ENDCOORDS(2,3*NATOMSALLOC))563:             ALLOCATE(ENDCOORDS(2,3*NATOMSALLOC))
564:             ENDCOORDS(1,1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1)564:             ENDCOORDS(1,1:3*NATOMSALLOC)=COORDS(1:3*NATOMSALLOC,1)
565:             ENDCOORDS(2,1:3*NATOMSALLOC)=FINISH(1:3*NATOMSALLOC)565:             ENDCOORDS(2,1:3*NATOMSALLOC)=FINISH(1:3*NATOMSALLOC)
566:             CALL ALIGN_DECIDE(ENDCOORDS(1,1:3*NATOMSALLOC),ENDCOORDS(2,1:3*NATOMSALLOC),NATOMS,DEBUG, 566:             CALL MINPERMDIST(ENDCOORDS(1,1:3*NATOMSALLOC),ENDCOORDS(2,1:3*NATOMSALLOC),NATOMS,DEBUG, 
567:      &                       BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,DIST2,RIGID,RMAT)567:      &                       BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,DIST2,RIGID,RMAT)
568:             COORDS(1:3*NATOMSALLOC,1)=ENDCOORDS(1,1:3*NATOMSALLOC)568:             COORDS(1:3*NATOMSALLOC,1)=ENDCOORDS(1,1:3*NATOMSALLOC)
569:             FINISH(1:3*NATOMSALLOC)=ENDCOORDS(2,1:3*NATOMSALLOC)569:             FINISH(1:3*NATOMSALLOC)=ENDCOORDS(2,1:3*NATOMSALLOC)
570:             DEALLOCATE(ENDCOORDS)570:             DEALLOCATE(ENDCOORDS)
571:          ENDIF571:          ENDIF
572:       ENDIF572:       ENDIF
573: !573: !
574: ! If this is a CSM optimisation we now have to multiply the number of atoms by the number of574: ! If this is a CSM optimisation we now have to multiply the number of atoms by the number of
575: ! group operations and replicate some coordinates and allowed permutations.575: ! group operations and replicate some coordinates and allowed permutations.
576: !576: !


r33306/make_conpot.f90 2017-09-13 18:30:48.874455353 +0100 r33305/make_conpot.f90 2017-09-13 18:30:51.662492405 +0100
 14: !   along with this program; if not, write to the Free Software 14: !   along with this program; if not, write to the Free Software
 15: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 15: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 16: ! 16: !
 17: SUBROUTINE MAKE_CONPOT(NCPFIT,MINCOORDS) 17: SUBROUTINE MAKE_CONPOT(NCPFIT,MINCOORDS)
 18: USE COMMONS, ONLY : INTCONSEP, NREPMAX, NREPULSIVE, CONDISTREF, REPCON, INTCONSTRAINTREP, & 18: USE COMMONS, ONLY : INTCONSEP, NREPMAX, NREPULSIVE, CONDISTREF, REPCON, INTCONSTRAINTREP, &
 19:   & REPCUT, NCONSTRAINT, CONI, CONJ, CONDISTREFLOCAL, INTCONMAX, CONACTIVE, & 19:   & REPCUT, NCONSTRAINT, CONI, CONJ, CONDISTREFLOCAL, INTCONMAX, CONACTIVE, &
 20:   & INTCONSTRAINREPCUT, INTREPSEP, REPI, REPJ, INTCONSTRAINTTOL, REPCUT, NREPI, NREPJ, NREPCUT, & 20:   & INTCONSTRAINREPCUT, INTREPSEP, REPI, REPJ, INTCONSTRAINTTOL, REPCUT, NREPI, NREPJ, NREPCUT, &
 21:   & NCONGEOM, CONGEOM, NNREPULSIVE, PERIODIC, RIGID, TWOD, & 21:   & NCONGEOM, CONGEOM, NNREPULSIVE, PERIODIC, RIGID, TWOD, &
 22:   & INTFROZEN, FREEZE, INTFREEZET, INTFREEZETOL, INTFREEZEMIN, CONIFIX, CONJFIX, CONDISTREFFIX, REPIFIX, REPJFIX, & 22:   & INTFROZEN, FREEZE, INTFREEZET, INTFREEZETOL, INTFREEZEMIN, CONIFIX, CONJFIX, CONDISTREFFIX, REPIFIX, REPJFIX, &
 23:   & REPCUTFIX, NCONGEOM, NREPULSIVEFIX, CONDATT, NCONSTRAINTFIX, CONCUTLOCAL, CONCUTFIX, CONCUT, & 23:   & REPCUTFIX, NCONGEOM, NREPULSIVEFIX, CONDATT, NCONSTRAINTFIX, CONCUTLOCAL, CONCUTFIX, CONCUT, &
 24:   & NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,MYUNIT, FASTOVERLAPT, BNB_ALIGNT 24:   & NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,MYUNIT
 25:  25: 
 26: IMPLICIT NONE  26: IMPLICIT NONE 
 27: DOUBLE PRECISION DF, D, RMAT(3,3), DISTANCE, D2 27: DOUBLE PRECISION DF, D, RMAT(3,3), DISTANCE, D2
 28: INTEGER :: J2,ISTAT,J1,J3,J4,NCPFIT,J5,NQCIFREEZE,NDUMMY,LUNIT,GETUNIT 28: INTEGER :: J2,ISTAT,J1,J3,J4,NCPFIT,J5,NQCIFREEZE,NDUMMY,LUNIT,GETUNIT
 29: INTEGER NCONFORNEWATOM 29: INTEGER NCONFORNEWATOM
 30: DOUBLE PRECISION :: NDIST, MINCOORDS(NCPFIT,3*NATOMS), DMIN, LINTCONSTRAINTTOL, & 30: DOUBLE PRECISION :: NDIST, MINCOORDS(NCPFIT,3*NATOMS), DMIN, LINTCONSTRAINTTOL, &
 31:   &                 LXYZ(6*NATOMS) 31:   &                 LXYZ(6*NATOMS)
 32: DOUBLE PRECISION :: DMOVED(NATOMS) 32: DOUBLE PRECISION :: DMOVED(NATOMS)
 33: LOGICAL ADDREP(NATOMS) 33: LOGICAL ADDREP(NATOMS)
 34: INTEGER NDIST1(NATOMS), NCYCLE, NUNCON1, DLIST(NATOMS) 34: INTEGER NDIST1(NATOMS), NCYCLE, NUNCON1, DLIST(NATOMS)
 41: ! 41: !
 42: ! If this is not the first call, and we are being passed two minima, 42: ! If this is not the first call, and we are being passed two minima,
 43: ! then we are doing an interpolation metric for a new pair of minima. 43: ! then we are doing an interpolation metric for a new pair of minima.
 44: ! We should optimise the permutational isomers on reference minimum 1 44: ! We should optimise the permutational isomers on reference minimum 1
 45: ! and then do the overall alignment with newmindist, fixing the 45: ! and then do the overall alignment with newmindist, fixing the
 46: ! permutational isomers. This should put the permutational isomers 46: ! permutational isomers. This should put the permutational isomers
 47: ! in register with the constraints, which were calculated for all 47: ! in register with the constraints, which were calculated for all
 48: ! the reference minima after aligning with the first. 48: ! the reference minima after aligning with the first.
 49: ! 49: !
 50:    IF ((CALLED.OR.CONDATT).AND.(NCPFIT.EQ.2)) THEN 50:    IF ((CALLED.OR.CONDATT).AND.(NCPFIT.EQ.2)) THEN
 51:       IF (FASTOVERLAPT .OR. BNB_ALIGNT) THEN 51:       CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),MINCOORDS(1,1:3*NATOMS),NATOMS,DEBUG, &
 52:          WRITE(*,*) "make_conpot> Warning: this subroutine has not been tested with the new alignment routines." 
 53:          WRITE(*,*) "If they are not suitable for your use case, please change make_conpot to call MINPERMDIST instead of ALIGN_DECIDE" 
 54:       ENDIF 
 55:       CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),MINCOORDS(1,1:3*NATOMS),NATOMS,DEBUG, & 
 56:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,D2,RIGID,RMAT) 52:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,D2,RIGID,RMAT)
 57:       CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DEBUG, & 53:       CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DEBUG, &
 58:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,D2,RIGID,RMAT) 54:   &                    BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,D,D2,RIGID,RMAT)
 59:       ! sn402: could this be changed to ALIGN_DECIDE as well? 
 60:       CALL NEWMINDIST(MINCOORDS(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DISTANCE, & 55:       CALL NEWMINDIST(MINCOORDS(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DISTANCE, &
 61:   &                   PERIODIC,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT) 56:   &                   PERIODIC,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
 62:    ENDIF 57:    ENDIF
 63: ENDIF 58: ENDIF
 64:  59: 
 65: NQCIFREEZE=0 60: NQCIFREEZE=0
 66: IF (FREEZE) THEN 61: IF (FREEZE) THEN
 67:    WRITE(MYUNIT, '(A)') ' make_conpot> ERROR *** QCI has not been coded for frozen atoms yet' 62:    WRITE(MYUNIT, '(A)') ' make_conpot> ERROR *** QCI has not been coded for frozen atoms yet'
 68:    STOP 63:    STOP
 69: ENDIF 64: ENDIF


r33306/mc.F 2017-09-13 18:30:49.122458643 +0100 r33305/mc.F 2017-09-13 18:30:51.886495384 +0100
3232: !     WRITE(MYUNIT,'(A,I5,2G17.7,3I5,L8)') 'J1,POTEL,EBEST,JBEST,J1-JBEST,NRELAX,RES1=',3232: !     WRITE(MYUNIT,'(A,I5,2G17.7,3I5,L8)') 'J1,POTEL,EBEST,JBEST,J1-JBEST,NRELAX,RES1=',
3233: !    1                                J1,POTEL,EBEST(JP),JBEST(JP),J1-JBEST(JP),NRELAX,RES13233: !    1                                J1,POTEL,EBEST(JP),JBEST(JP),J1-JBEST(JP),NRELAX,RES1
3234:       RES2=.FALSE.3234:       RES2=.FALSE.
3235: !     IF ((.NOT.RES1).AND.AVOID) THEN3235: !     IF ((.NOT.RES1).AND.AVOID) THEN
3236: !     PRINT*,'RES1,AVOID,J1,JBEST(JP)=',RES1,AVOID,J1,JBEST(JP)3236: !     PRINT*,'RES1,AVOID,J1,JBEST(JP)=',RES1,AVOID,J1,JBEST(JP)
3237:       IF ((.NOT.RES1).AND.AVOID.AND.(J1.EQ.JBEST(JP)).AND.(NMSBSAVE.GT.0)) THEN ! best minimum has just changed.3237:       IF ((.NOT.RES1).AND.AVOID.AND.(J1.EQ.JBEST(JP)).AND.(NMSBSAVE.GT.0)) THEN ! best minimum has just changed.
3238:          FCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)3238:          FCOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)
3239:          savedloop: DO J2=1,MIN(NMSBSAVE,MAXSAVE)3239:          savedloop: DO J2=1,MIN(NMSBSAVE,MAXSAVE)
3240:             IF (DEBUG) THEN3240:             IF (DEBUG) THEN
3241:                XMSB(1:3*NATOMS)=MSBCOORDS(1:3*NATOMS,J2)3241:                XMSB(1:3*NATOMS)=MSBCOORDS(1:3*NATOMS,J2)
3242:                CALL ALIGN_DECIDE(FCOORDS,XMSB,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DISTANCE,DIST2,RIGID,RMAT)3242:                CALL MINPERMDIST(FCOORDS,XMSB,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DISTANCE,DIST2,RIGID,RMAT)
3243:                WRITE(MYUNIT,'(A,F10.3)')'newres> minimum distance is ',DISTANCE3243:                WRITE(MYUNIT,'(A,F10.3)')'newres> minimum distance is ',DISTANCE
3244:             ENDIF3244:             ENDIF
3245: 3245: 
3246: !3246: !
3247: !  Bipartite matching routine for permutations. Coordinates in FCOORDS do not change3247: !  Bipartite matching routine for permutations. Coordinates in FCOORDS do not change
3248: !  but the coordinates in XMSB do. DISTANCE is the distance in this case.3248: !  but the coordinates in XMSB do. DISTANCE is the distance in this case.
3249: !3249: !
3250: ! If the energy is lower no reseeding regardless of separation ? DJW3250: ! If the energy is lower no reseeding regardless of separation ? DJW
3251: !           PRINT '(A,2G20.10,L8)','POTEL,MSBE(J2)-ECONV,POTEL.LT.MSBE(J2)-ECONV=',3251: !           PRINT '(A,2G20.10,L8)','POTEL,MSBE(J2)-ECONV,POTEL.LT.MSBE(J2)-ECONV=',
3252: !    1                              POTEL,MSBE(J2)-ECONV,POTEL.LT.MSBE(J2)-ECONV3252: !    1                              POTEL,MSBE(J2)-ECONV,POTEL.LT.MSBE(J2)-ECONV
3253:             IF (POTEL.LT.MSBE(J2)-ECONV) CYCLE savedloop3253:             IF (POTEL.LT.MSBE(J2)-ECONV) CYCLE savedloop
3254:             XMSB(1:3*NATOMS)=MSBCOORDS(1:3*NATOMS,J2)3254:             XMSB(1:3*NATOMS)=MSBCOORDS(1:3*NATOMS,J2)
3255:             CALL ALIGN_DECIDE(FCOORDS,XMSB,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DISTANCE,DIST2,RIGID,RMAT)3255:             CALL MINPERMDIST(FCOORDS,XMSB,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DISTANCE,DIST2,RIGID,RMAT)
3256:             IF (DISTANCE.LT.AVOIDDIST) THEN3256:             IF (DISTANCE.LT.AVOIDDIST) THEN
3257:                RES2=.TRUE.3257:                RES2=.TRUE.
3258:                WRITE(MYUNIT,'(A,G20.10,A,I6,A,G20.10,A,F10.3)') 'newres> Minimum energy ',POTEL,3258:                WRITE(MYUNIT,'(A,G20.10,A,I6,A,G20.10,A,F10.3)') 'newres> Minimum energy ',POTEL,
3259:      &                                    ' is too close to saved structure ',3259:      &                                    ' is too close to saved structure ',
3260:      &                                 J2,' with energy ',MSBE(J2),' dist=',DISTANCE3260:      &                                 J2,' with energy ',MSBE(J2),' dist=',DISTANCE
3261:                GOTO 203261:                GOTO 20
3262:             ENDIF3262:             ENDIF
3263:          ENDDO savedloop3263:          ENDDO savedloop
3264:       ENDIF3264:       ENDIF
3265: 3265: 


r33306/potential.f90 2017-09-13 18:30:49.350461675 +0100 r33305/potential.f90 2017-09-13 18:30:52.110498361 +0100
972:          CSMGPSAVE=CSMGP972:          CSMGPSAVE=CSMGP
973:          CSMGP=CSMGUIDEGP973:          CSMGP=CSMGUIDEGP
974:          CSMGPINDEXSAVE=CSMGPINDEX974:          CSMGPINDEXSAVE=CSMGPINDEX
975:          CSMGPINDEX=CSMGUIDEGPINDEX975:          CSMGPINDEX=CSMGUIDEGPINDEX
976:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPGUIDE(1:3, 1:3, 1:2*CSMGPINDEX)976:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPGUIDE(1:3, 1:3, 1:2*CSMGPINDEX)
977:       END IF977:       END IF
978:       IF (PERMDIST) THEN978:       IF (PERMDIST) THEN
979:          DO J1=1, CSMGPINDEX979:          DO J1=1, CSMGPINDEX
980:             XTEMP(1:3*NATOMS)=X(1:3*NATOMS)980:             XTEMP(1:3*NATOMS)=X(1:3*NATOMS)
981:             CALL CSMROT(XTEMP, DUMMY, 1, J1)981:             CALL CSMROT(XTEMP, DUMMY, 1, J1)
982:             CALL ALIGN_DECIDE(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)982:             CALL MINPERMDIST(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)
983:             CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation983:             CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation
984:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)984:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)
985:          END DO985:          END DO
986:       ELSE986:       ELSE
987:          DO J1=1, CSMGPINDEX987:          DO J1=1, CSMGPINDEX
988:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=X(1:3*NATOMS)988:             CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=X(1:3*NATOMS)
989:          END DO989:          END DO
990:       END IF990:       END IF
991:       CALL CSMMIN(X, EREAL, CSMRMS, CSMIT)991:       CALL CSMMIN(X, EREAL, CSMRMS, CSMIT)
992: 992: 
1023:          END DO1023:          END DO
1024:          CSMAV(1:3*NATOMS)=CSMAV(1:3*NATOMS)/CSMGPINDEX1024:          CSMAV(1:3*NATOMS)=CSMAV(1:3*NATOMS)/CSMGPINDEX
1025: !1025: !
1026: !  Check the CSM for the averaged structure. It should be zero if this structure has the1026: !  Check the CSM for the averaged structure. It should be zero if this structure has the
1027: !  right point group. Need to reset CSMIMAGES and CSMNORM temporarily.1027: !  right point group. Need to reset CSMIMAGES and CSMNORM temporarily.
1028: !1028: !
1029:          IF (PERMDIST) THEN1029:          IF (PERMDIST) THEN
1030:             DO J1=1, CSMGPINDEX1030:             DO J1=1, CSMGPINDEX
1031:                XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)1031:                XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)
1032:                CALL CSMROT(XTEMP, DUMMY, 1, J1)1032:                CALL CSMROT(XTEMP, DUMMY, 1, J1)
1033:                CALL ALIGN_DECIDE(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, DUMMY2, DIST2, RIGID, RMAT)1033:                CALL MINPERMDIST(XTEMP, DUMMY, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, DUMMY2, DIST2, RIGID, RMAT)
1034:                CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation1034:                CALL CSMROT(DUMMY, XTEMP,-1, J1) ! need to rotate the permuted rotated images back to the reference orientation
1035:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)1035:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)
1036:             END DO1036:             END DO
1037:          ELSE1037:          ELSE
1038:             DO J1=1, CSMGPINDEX1038:             DO J1=1, CSMGPINDEX
1039:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=CSMAV(1:3*NATOMS)1039:                CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=CSMAV(1:3*NATOMS)
1040:             END DO1040:             END DO
1041:          END IF1041:          END IF
1042:          AA(1)=0.0D0; AA(2)=0.0D0; AA(3)=6.283185307D0 ! should give an identity matrix1042:          AA(1)=0.0D0; AA(2)=0.0D0; AA(3)=6.283185307D0 ! should give an identity matrix
1043:          SAVECSMPMAT(1:3, 1:3)=CSMPMAT(1:3, 1:3)1043:          SAVECSMPMAT(1:3, 1:3)=CSMPMAT(1:3, 1:3)
1053:          CSMIMAGES(1:3*NATOMS*CSMGPINDEX)=SAVECSMIMAGES(1:3*NATOMS*CSMGPINDEX)1053:          CSMIMAGES(1:3*NATOMS*CSMGPINDEX)=SAVECSMIMAGES(1:3*NATOMS*CSMGPINDEX)
1054:          WRITE(MYUNIT,'(A, 2G20.10)') 'potential> CSM values for reference structure and average=', EREAL, AVVAL1054:          WRITE(MYUNIT,'(A, 2G20.10)') 'potential> CSM values for reference structure and average=', EREAL, AVVAL
1055:       END IF ! DEBUG1055:       END IF ! DEBUG
1056:       IF (CSMDOGUIDET) THEN ! undo guiding changes1056:       IF (CSMDOGUIDET) THEN ! undo guiding changes
1057:          CSMGP=CSMGPSAVE1057:          CSMGP=CSMGPSAVE
1058:          CSMGPINDEX=CSMGPINDEXSAVE1058:          CSMGPINDEX=CSMGPINDEXSAVE
1059:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)1059:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)
1060:          CSMNORM=CSMNORMSAVE1060:          CSMNORM=CSMNORMSAVE
1061:       END IF1061:       END IF
1062: 1062: 
1063:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT .OR. ALIGNT) THEN1063:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN
1064: !  EREAL is the distance in this case1064: !  EREAL is the distance in this case
1065:       write(*,*) "Calling ALIGN_DECIDE from POTENTIAL, PERIODIC=",PERIODIC1065:       CALL MINPERMDIST(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)
1066:       CALL ALIGN_DECIDE(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT) 
1067:       IF(DEBUG) THEN 1066:       IF(DEBUG) THEN 
1068:         WRITE(MYUNIT,'(A)') 'Rotation matrix:'1067:         WRITE(MYUNIT,'(A)') 'Rotation matrix:'
1069:         WRITE(MYUNIT,*) RMAT(:,:)1068:         WRITE(MYUNIT,*) RMAT(:,:)
1070:       END IF1069:       END IF
1071: 1070:       
1072:    ELSE IF (MKTRAPT) THEN1071:    ELSE IF (MKTRAPT) THEN
1073:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)1072:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)
1074:    ELSE IF (BLJCLUSTER) THEN1073:    ELSE IF (BLJCLUSTER) THEN
1075:       CALL RAD(X, GRAD, EREAL, GRADT)1074:       CALL RAD(X, GRAD, EREAL, GRADT)
1076:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)1075:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)
1077: 1076: 
1078:    ELSE IF (BLJCLUSTER_NOCUT) THEN1077:    ELSE IF (BLJCLUSTER_NOCUT) THEN
1079: ! ds656> Binary LJ without cutoff1078: ! ds656> Binary LJ without cutoff
1080:       CALL RAD(X, GRAD, EREAL, GRADT)1079:       CALL RAD(X, GRAD, EREAL, GRADT)
1081:       CALL BLJ_CLUST(X, GRAD, EREAL, GRADT)1080:       CALL BLJ_CLUST(X, GRAD, EREAL, GRADT)


r33306/quench.F 2017-09-13 18:30:49.578464706 +0100 r33305/quench.F 2017-09-13 18:30:52.478503265 +0100
177:          END IF177:          END IF
178:          POTEL=EREAL178:          POTEL=EREAL
179:          IF (.NOT.CFLAG) WRITE(MYUNIT,'(A,I7,A)') ' WARNING - compressed quench ',NQ(NP),'  did not converge'179:          IF (.NOT.CFLAG) WRITE(MYUNIT,'(A,I7,A)') ' WARNING - compressed quench ',NQ(NP),'  did not converge'
180:          WRITE(MYUNIT,'(A,I7,A,F20.10,A,I5,A,F15.7,A,I4,A,F12.2)') 'Comp Q ',NQ(NP),' energy=',180:          WRITE(MYUNIT,'(A,I7,A,F20.10,A,I5,A,F15.7,A,I4,A,F12.2)') 'Comp Q ',NQ(NP),' energy=',
181:      1              POTEL,' steps=',ITER,' RMS force=',RMS181:      1              POTEL,' steps=',ITER,' RMS force=',RMS
182:       ENDIF182:       ENDIF
183: 183: 
184:       IF (.NOT.PERCOLATET) COMPON=.FALSE.184:       IF (.NOT.PERCOLATET) COMPON=.FALSE.
185: 185: 
186: 186: 
187: 10    IF (PERMOPT.OR.PERMINVOPT.OR.DISTOPT.OR.ALIGNT) THEN ! lb415187: 10    IF (PERMOPT.OR.PERMINVOPT.OR.DISTOPT) THEN ! lb415
188:          !IF ( NQ(NP) .eq. 1) THEN188:          !IF ( NQ(NP) .eq. 1) THEN
189:          IF (DUMPT) THEN189:          IF (DUMPT) THEN
190:             IF (NP.EQ.1) WRITE(MYUNIT,'(A,4I6)') 'quench> initial NP,DUMPXYZUNIT=',NP,DUMPXYZUNIT(NP)190:             IF (NP.EQ.1) WRITE(MYUNIT,'(A,4I6)') 'quench> initial NP,DUMPXYZUNIT=',NP,DUMPXYZUNIT(NP)
191:             WRITE(DUMPXYZUNIT(NP),'(I6)') NDUMMY191:             WRITE(DUMPXYZUNIT(NP),'(I6)') NDUMMY
192:             WRITE(DUMPXYZUNIT(NP),'(A,I6)') 'quench> initial points before quench ',NQ(NP)192:             WRITE(DUMPXYZUNIT(NP),'(A,I6)') 'quench> initial points before quench ',NQ(NP)
193:             WRITE(DUMPXYZUNIT(NP),'(A,3G20.10)') ('LA ',P(3*(J2-1)+1),P(3*(J2-1)+2),P(3*(J2-1)+3),J2=1,NATOMS)193:             WRITE(DUMPXYZUNIT(NP),'(A,3G20.10)') ('LA ',P(3*(J2-1)+1),P(3*(J2-1)+2),P(3*(J2-1)+3),J2=1,NATOMS)
194:          ENDIF194:          ENDIF
195:          CALL POTENTIAL(P,GRAD,EREAL,.FALSE.,.FALSE.)195:          CALL POTENTIAL(P,GRAD,EREAL,.FALSE.,.FALSE.)
196:          CFLAG=.TRUE.196:          CFLAG=.TRUE.
197: !        ITER=1197: !        ITER=1
198: !        RMS=0.0D0198: !        RMS=0.0D0
199:          RMS=CSMRMS199:          RMS=CSMRMS
200:          ITER=CSMIT200:          ITER=CSMIT
201:          !ELSE201:          !ELSE
202:          !   CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP) ! minimize structure202:          !   CALL MYLBFGS(NOPT,MUPDATE,P,.FALSE.,GMAX,CFLAG,EREAL,MAXIT,ITER,.TRUE.,NP) ! minimize structure
203:          !   write(*,*) 'permdist mylbfgs', EREAL, ITER, RMS203:          !   write(*,*) 'permdist mylbfgs', EREAL, ITER, RMS
204:          !   POTEL=EREAL204:          !   POTEL=EREAL
205:          !   IF (.NOT.CFLAG) WRITE(MYUNIT,'(A,I7,A)') 'WARNING - Quench ',NQ(NP),'  did not converge'205:          !   IF (.NOT.CFLAG) WRITE(MYUNIT,'(A,I7,A)') 'WARNING - Quench ',NQ(NP),'  did not converge'
206:          !   DO II=1,NSAVE206:          !   DO II=1,NSAVE
207:          !      IF ( II .GE. NQ(NP) ) EXIT ! There's no need to check further, there's nothing207:          !      IF ( II .GE. NQ(NP) ) EXIT ! There's no need to check further, there's nothing
208:          !      CALL ALIGN_DECIDE(P,QMINP(II,:),NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY,DIST2,RIGID,RMAT)208:          !      CALL MINPERMDIST(P,QMINP(II,:),NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY,DIST2,RIGID,RMAT)
209:          !      write(*,*) DUMMY, 'dummy',ii209:          !      write(*,*) DUMMY, 'dummy',ii
210:          !      IF (DUMMY .LT. 0.5D0) THEN210:          !      IF (DUMMY .LT. 0.5D0) THEN
211:          !         !DO NOT ACCEPT THIS QUENCH211:          !         !DO NOT ACCEPT THIS QUENCH
212:          !         WRITE(MYUNIT,*) 'This quench ended in a known minimum. It won`t be counted.'212:          !         WRITE(MYUNIT,*) 'This quench ended in a known minimum. It won`t be counted.'
213:          !         RETURN213:          !         RETURN
214:          !      ENDIF214:          !      ENDIF
215:          !   ENDDO215:          !   ENDDO
216:          !ENDIF 216:          !ENDIF 
217:       ELSEIF (MODEL1T) THEN217:       ELSEIF (MODEL1T) THEN
218:          CALL MODEL1(P,GRAD,EREAL,QE,QX)218:          CALL MODEL1(P,GRAD,EREAL,QE,QX)
891:          ENDIF891:          ENDIF
892:       ENDIF892:       ENDIF
893: !893: !
894: !  If we are provided with target minimum coordinates in file coords.target then894: !  If we are provided with target minimum coordinates in file coords.target then
895: !  calculate the minimum distances. May be useful for algorithm development.895: !  calculate the minimum distances. May be useful for algorithm development.
896: !  If we get close, we don;t want to escape without a hit!896: !  If we get close, we don;t want to escape without a hit!
897: !897: !
898: !     IF (ALLOCATED(TCOORDS)) THEN898: !     IF (ALLOCATED(TCOORDS)) THEN
899: !        DO J1=1,NTARGETS899: !        DO J1=1,NTARGETS
900: !           TMPCOORDS(1:3*NATOMS)=TCOORDS(J1,1:3*NATOMS)900: !           TMPCOORDS(1:3*NATOMS)=TCOORDS(J1,1:3*NATOMS)
901: !           CALL ALIGN_DECIDE(P,TMPCOORDS,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY,DIST2,RIGID)901: !           CALL MINPERMDIST(P,TMPCOORDS,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY,DIST2,RIGID)
902: !           WRITE(MYUNIT, '(A,I5,A,F15.3,A,F15.3,A,F20.10)') 'for target structure ',J1,' dist=',DUMMY,' dist2=',DIST2,' V=',POTEL902: !           WRITE(MYUNIT, '(A,I5,A,F15.3,A,F15.3,A,F20.10)') 'for target structure ',J1,' dist=',DUMMY,' dist2=',DIST2,' V=',POTEL
903: !        ENDDO903: !        ENDDO
904: !        DO J1=1,MIN(NMSBSAVE,MAXSAVE)904: !        DO J1=1,MIN(NMSBSAVE,MAXSAVE)
905: !           TMPCOORDS(1:3*NATOMS)=MSBCOORDS(1:3*NATOMS,J1)905: !           TMPCOORDS(1:3*NATOMS)=MSBCOORDS(1:3*NATOMS,J1)
906: !           CALL ALIGN_DECIDE(P,TMPCOORDS,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY,DIST2,RIGID)906: !           CALL MINPERMDIST(P,TMPCOORDS,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DUMMY,DIST2,RIGID)
907: !           PRINT '(A,I5,A,F15.3,A,F15.3,A,F20.10)','for taboo  structure ',J1,' dist=',DUMMY,' dist2=',DIST2,' V=',POTEL907: !           PRINT '(A,I5,A,F15.3,A,F15.3,A,F20.10)','for taboo  structure ',J1,' dist=',DUMMY,' dist2=',DIST2,' V=',POTEL
908: !        ENDDO908: !        ENDDO
909: !     ENDIF909: !     ENDIF
910: !910: !
911: !  NORESET true does not set the configuration point to the quench geometry911: !  NORESET true does not set the configuration point to the quench geometry
912: !  A relaxed frozen core does not get saved, but the lowest minima are saved912: !  A relaxed frozen core does not get saved, but the lowest minima are saved
913: !  by GSAVEIT.913: !  by GSAVEIT.
914: !914: !
915:       IF (.NOT.NORESET) THEN915:       IF (.NOT.NORESET) THEN
916:          DO J1=1,3*(NATOMS-NSEED)916:          DO J1=1,3*(NATOMS-NSEED)


r33306/ratio.f90 2017-09-13 18:30:49.794467576 +0100 r33305/ratio.f90 2017-09-13 18:30:52.694506122 +0100
  7: ! NULLMOVE also acts as the Metropolis test.  7: ! NULLMOVE also acts as the Metropolis test.
  8: SUBROUTINE NULLMOVE(NEWCOORDS,OLDCOORDS,ATEST,NULLMOVES,JP)  8: SUBROUTINE NULLMOVE(NEWCOORDS,OLDCOORDS,ATEST,NULLMOVES,JP)
  9: USE COMMONS, ONLY: NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, GEOMDIFFTOL, RIGID, TWOD, TEMP, MYUNIT, DEBUG, NPAR  9: USE COMMONS, ONLY: NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, GEOMDIFFTOL, RIGID, TWOD, TEMP, MYUNIT, DEBUG, NPAR
 10: IMPLICIT NONE 10: IMPLICIT NONE
 11: DOUBLE PRECISION  :: NEWCOORDS(3*NATOMS), OLDCOORDS(3*NATOMS), OLDCOORDSCP(3*NATOMS) 11: DOUBLE PRECISION  :: NEWCOORDS(3*NATOMS), OLDCOORDS(3*NATOMS), OLDCOORDSCP(3*NATOMS)
 12: DOUBLE PRECISION  :: RMAT(3,3), DISTANCE, DIST2, ENEW, EOLD, GRAD(3*NATOMS), DPRAND 12: DOUBLE PRECISION  :: RMAT(3,3), DISTANCE, DIST2, ENEW, EOLD, GRAD(3*NATOMS), DPRAND
 13: INTEGER           :: NULLMOVES(NPAR), JP 13: INTEGER           :: NULLMOVES(NPAR), JP
 14: LOGICAL           :: ATEST  14: LOGICAL           :: ATEST 
 15:  15: 
 16: OLDCOORDSCP(:)=OLDCOORDS(:) ! make local copy to avoid returning modified markov coordinates 16: OLDCOORDSCP(:)=OLDCOORDS(:) ! make local copy to avoid returning modified markov coordinates
 17: CALL ALIGN_DECIDE(NEWCOORDS,OLDCOORDSCP,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DISTANCE,DIST2,RIGID,RMAT) 17: CALL MINPERMDIST(NEWCOORDS,OLDCOORDSCP,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,DISTANCE,DIST2,RIGID,RMAT)
 18: CALL POTENTIAL(OLDCOORDSCP,GRAD,EOLD,.FALSE.,.FALSE.) ! possibly 18: CALL POTENTIAL(OLDCOORDSCP,GRAD,EOLD,.FALSE.,.FALSE.) ! possibly
 19: CALL POTENTIAL(NEWCOORDS,GRAD,ENEW,.FALSE.,.FALSE.) ! redundant 19: CALL POTENTIAL(NEWCOORDS,GRAD,ENEW,.FALSE.,.FALSE.) ! redundant
 20:  20: 
 21: IF (DISTANCE<GEOMDIFFTOL) THEN ! you didn't move to a new minimum 21: IF (DISTANCE<GEOMDIFFTOL) THEN ! you didn't move to a new minimum
 22:    ATEST=.FALSE. 22:    ATEST=.FALSE.
 23:    NULLMOVES(JP)=NULLMOVES(JP)+1 23:    NULLMOVES(JP)=NULLMOVES(JP)+1
 24:    IF (DEBUG) WRITE(MYUNIT,*)ATEST,"DISTFAIL, DISTANCE=",DISTANCE 24:    IF (DEBUG) WRITE(MYUNIT,*)ATEST,"DISTFAIL, DISTANCE=",DISTANCE
 25: ELSE IF ((EXP(-(ENEW-EOLD)/MAX(TEMP(JP),1.D-100))<DPRAND()).OR.(ENEW<-1.0D6)) THEN ! you failed the Metropolis test 25: ELSE IF ((EXP(-(ENEW-EOLD)/MAX(TEMP(JP),1.D-100))<DPRAND()).OR.(ENEW<-1.0D6)) THEN ! you failed the Metropolis test
 26:    ATEST=.FALSE. 26:    ATEST=.FALSE.
 27:    IF (DEBUG) WRITE(MYUNIT,*)ATEST,"EFAIL, EDIFF=",ENEW-EOLD 27:    IF (DEBUG) WRITE(MYUNIT,*)ATEST,"EFAIL, EDIFF=",ENEW-EOLD


r33306/symmetrycsm.f90 2017-09-13 18:30:50.014470500 +0100 r33305/symmetrycsm.f90 2017-09-13 18:30:52.914509048 +0100
 81: ! 81: !
 82:  82: 
 83: CSMBEST=1.0D100 83: CSMBEST=1.0D100
 84: OLCOORDS(1:3*NATOMS)=LCOORDS(3*NATOMS) 84: OLCOORDS(1:3*NATOMS)=LCOORDS(3*NATOMS)
 85: CSMSTEP=0.5D0 85: CSMSTEP=0.5D0
 86: DO J2=1,CSMSTEPS 86: DO J2=1,CSMSTEPS
 87:    IF (PERMDIST) THEN 87:    IF (PERMDIST) THEN
 88:       DO J1=1,CSMGPINDEX 88:       DO J1=1,CSMGPINDEX
 89:          XTEMP(1:3*NATOMS)=LCOORDS(1:3*NATOMS) 89:          XTEMP(1:3*NATOMS)=LCOORDS(1:3*NATOMS)
 90:          CALL CSMROT(XTEMP,DUMMY,1,J1) 90:          CALL CSMROT(XTEMP,DUMMY,1,J1)
 91:          CALL ALIGN_DECIDE(XTEMP,DUMMY,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,EREAL,DIST2,RIGID,RMAT) 91:          CALL MINPERMDIST(XTEMP,DUMMY,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,EREAL,DIST2,RIGID,RMAT)
 92:          CALL CSMROT(DUMMY,XTEMP,-1,J1) ! need to rotate the permuted rotated images back to the reference orientation 92:          CALL CSMROT(DUMMY,XTEMP,-1,J1) ! need to rotate the permuted rotated images back to the reference orientation
 93:          CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS) 93:          CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)
 94:       ENDDO 94:       ENDDO
 95:    ELSE 95:    ELSE
 96:       DO J1=1,CSMGPINDEX 96:       DO J1=1,CSMGPINDEX
 97:          CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=LCOORDS(1:3*NATOMS) 97:          CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=LCOORDS(1:3*NATOMS)
 98:       ENDDO 98:       ENDDO
 99:    ENDIF 99:    ENDIF
100: !100: !
101: ! At this point we have the best permutations for each point group operation,101: ! At this point we have the best permutations for each point group operation,
200: ! !200: ! !
201: ! NDONE=0201: ! NDONE=0
202: ! 5 CONTINUE202: ! 5 CONTINUE
203: ! CALL CENTRE2(CSMAV)203: ! CALL CENTRE2(CSMAV)
204: ! NDONE=NDONE+1204: ! NDONE=NDONE+1
205: ! IF (PERMDIST) THEN205: ! IF (PERMDIST) THEN
206: !    CSMVALUE=0.0D0206: !    CSMVALUE=0.0D0
207: !    DO J1=1,CSMGPINDEX207: !    DO J1=1,CSMGPINDEX
208: !       XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)208: !       XTEMP(1:3*NATOMS)=CSMAV(1:3*NATOMS)
209: !       CALL CSMROT(XTEMP,DUMMY,1,J1)209: !       CALL CSMROT(XTEMP,DUMMY,1,J1)
210: !       CALL ALIGN_DECIDE(XTEMP,DUMMY,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,EREAL,DIST2,RIGID,RMAT)210: !       CALL MINPERMDIST(XTEMP,DUMMY,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,PERIODIC,TWOD,EREAL,DIST2,RIGID,RMAT)
211: !       CSMVALUE=CSMVALUE+EREAL**2211: !       CSMVALUE=CSMVALUE+EREAL**2
212: !       CALL CSMROT(DUMMY,XTEMP,-1,J1) ! need to rotate the permuted rotated images back to the reference orientation212: !       CALL CSMROT(DUMMY,XTEMP,-1,J1) ! need to rotate the permuted rotated images back to the reference orientation
213: !       CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)213: !       CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=XTEMP(1:3*NATOMS)
214: !    ENDDO214: !    ENDDO
215: ! ELSE215: ! ELSE
216: !    DO J1=1,CSMGPINDEX216: !    DO J1=1,CSMGPINDEX
217: !       CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=CSMAV(1:3*NATOMS)217: !       CSMIMAGES(1+3*NATOMS*(J1-1):3*NATOMS*J1)=CSMAV(1:3*NATOMS)
218: !    ENDDO218: !    ENDDO
219: ! ENDIF219: ! ENDIF
220: ! CSMNORM=0.0D0220: ! CSMNORM=0.0D0


r33306/takestep.f 2017-09-13 18:30:50.234473424 +0100 r33305/takestep.f 2017-09-13 18:30:53.138512022 +0100
 45: !        WRITE(177,'(I6)') NATOMS 45: !        WRITE(177,'(I6)') NATOMS
 46: !        WRITE(MYUNIT,'(A)') ' start of takestep coords ' 46: !        WRITE(MYUNIT,'(A)') ' start of takestep coords '
 47: !        WRITE(MYUNIT,'(A2,2X,3G20.10)') ('LA',COORDS(3*(J2-1)+1,NP),COORDS(3*(J2-1)+2,NP), 47: !        WRITE(MYUNIT,'(A2,2X,3G20.10)') ('LA',COORDS(3*(J2-1)+1,NP),COORDS(3*(J2-1)+2,NP),
 48: !    &                                 COORDS(3*(J2-1)+3,NP),J2=1,NATOMS) 48: !    &                                 COORDS(3*(J2-1)+3,NP),J2=1,NATOMS)
 49: !        WRITE(177,'(I6)') NATOMS 49: !        WRITE(177,'(I6)') NATOMS
 50: !        WRITE(MYUNIT,'(A)') ' start of takestep coordso' 50: !        WRITE(MYUNIT,'(A)') ' start of takestep coordso'
 51: !        WRITE(MYUNIT,'(A2,2X,3G20.10)') ('LA',COORDSO(3*(J2-1)+1,NP),COORDSO(3*(J2-1)+2,NP), 51: !        WRITE(MYUNIT,'(A2,2X,3G20.10)') ('LA',COORDSO(3*(J2-1)+1,NP),COORDSO(3*(J2-1)+2,NP),
 52: !    &                                 COORDSO(3*(J2-1)+3,NP),J2=1,NATOMS) 52: !    &                                 COORDSO(3*(J2-1)+3,NP),J2=1,NATOMS)
 53: !        CLOSE(177) 53: !        CLOSE(177)
 54:  54: 
 55:       IF ((PERMOPT.OR.PERMINVOPT.OR.DISTOPT.OR.ALIGNT).AND.PERIODIC) RETURN  55:       IF ((PERMOPT.OR.PERMINVOPT.OR.DISTOPT).AND.PERIODIC) RETURN 
 56:       IF ((PERMOPT.OR.PERMINVOPT.OR.ALIGNT).AND.MKTRAPT) RETURN  56:       IF ((PERMOPT.OR.PERMINVOPT).AND.MKTRAPT) RETURN 
 57:        57:       
 58:       IF (MODEL1T) THEN 58:       IF (MODEL1T) THEN
 59:          RANDOM=(DPRAND()-0.5D0)*2.0D0 59:          RANDOM=(DPRAND()-0.5D0)*2.0D0
 60:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM 60:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM
 61:          RETURN 61:          RETURN
 62:       ENDIF 62:       ENDIF
 63:       IF (MLQT.OR.MLPVB3T.OR.MLP3T.OR.ORBITALS) THEN 63:       IF (MLQT.OR.MLPVB3T.OR.MLP3T.OR.ORBITALS) THEN
 64:          DO J1=1,NATOMS 64:          DO J1=1,NATOMS
 65:             IF (FROZEN(J1)) CYCLE 65:             IF (FROZEN(J1)) CYCLE
 66:             RANDOM=(DPRAND()-0.5D0)*2.0D0 66:             RANDOM=(DPRAND()-0.5D0)*2.0D0
135: !     ENDIF135: !     ENDIF
136: 136: 
137:       IF (WENZEL) THEN137:       IF (WENZEL) THEN
138: 11       RANDOM=(DPRAND()-0.5D0)*2.0D0138: 11       RANDOM=(DPRAND()-0.5D0)*2.0D0
139:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM139:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM
140:          IF ((COORDS(1,NP).GT.1.0D0).OR.(COORDS(1,NP).LT.0.0D0)) GOTO 11140:          IF ((COORDS(1,NP).GT.1.0D0).OR.(COORDS(1,NP).LT.0.0D0)) GOTO 11
141: 12       RANDOM=(DPRAND()-0.5D0)*2.0D0141: 12       RANDOM=(DPRAND()-0.5D0)*2.0D0
142:          COORDS(2,NP)=COORDSO(2,NP)+STEP(NP)*RANDOM142:          COORDS(2,NP)=COORDSO(2,NP)+STEP(NP)*RANDOM
143:          IF ((COORDS(2,NP).GT.1.0D0).OR.(COORDS(2,NP).LT.0.0D0)) GOTO 12143:          IF ((COORDS(2,NP).GT.1.0D0).OR.(COORDS(2,NP).LT.0.0D0)) GOTO 12
144:          RETURN144:          RETURN
145:       ELSE IF (PERMOPT.OR.PERMINVOPT.OR.ALIGNT.OR.(CSMT.AND.(.NOT.SYMMETRIZECSM))) THEN145:       ELSE IF (PERMOPT.OR.PERMINVOPT.OR.(CSMT.AND.(.NOT.SYMMETRIZECSM))) THEN
146: 146: 
147:           P(:)=ROT_SMALL_RANDOM_AA(STEP(NP))147:           P(:)=ROT_SMALL_RANDOM_AA(STEP(NP))
148:           M(:,:)=ROT_AA2MX(P(:))148:           M(:,:)=ROT_AA2MX(P(:))
149: 149: 
150:           IF (RIGID) NATOMS=NATOMS/2150:           IF (RIGID) NATOMS=NATOMS/2
151:           DO J1=1,NATOMS151:           DO J1=1,NATOMS
152:              J2=3*J1152:              J2=3*J1
153:              J3=3*(J1+NATOMS)153:              J3=3*(J1+NATOMS)
154:              COORDS(J2-2:J2,NP)=MATMUL(M(:,:),COORDS(J2-2:J2,NP))154:              COORDS(J2-2:J2,NP)=MATMUL(M(:,:),COORDS(J2-2:J2,NP))
155:              IF (RIGID) COORDS(J3-2:J3,NP)=ROT_ROTATE_AA(COORDS(J3-2:J3,NP),P(:))155:              IF (RIGID) COORDS(J3-2:J3,NP)=ROT_ROTATE_AA(COORDS(J3-2:J3,NP),P(:))
407:          ENDIF407:          ENDIF
408: !408: !
409: !  Angular move block.409: !  Angular move block.
410: !  If NORESET is .TRUE. then VAT won;t be set, so we should skip this block.410: !  If NORESET is .TRUE. then VAT won;t be set, so we should skip this block.
411: !411: !
412: !        IF (J1.EQ.JMAX) WRITE(MYUNIT,'(A,I6,4F15.5)') 'JMAX,VAT,ASTEP(NP),VMIN,prod=',JMAX,VAT(J1,NP), 412: !        IF (J1.EQ.JMAX) WRITE(MYUNIT,'(A,I6,4F15.5)') 'JMAX,VAT,ASTEP(NP),VMIN,prod=',JMAX,VAT(J1,NP), 
413: !    &                                    ASTEP(NP),VMIN,ASTEP(NP)*VMIN413: !    &                                    ASTEP(NP),VMIN,ASTEP(NP)*VMIN
414:          IF (((VAT(J1,NP).GT.ASTEP(NP)*VMIN).AND.(J1.EQ.JMAX)).AND.(.NOT.BLNT).AND.!(.NOT.RIGID).AND.414:          IF (((VAT(J1,NP).GT.ASTEP(NP)*VMIN).AND.(J1.EQ.JMAX)).AND.(.NOT.BLNT).AND.!(.NOT.RIGID).AND.
415:      &         (.NOT.DIFFRACTT).AND.(.NOT.GAUSST).AND.(.NOT.PERCOLATET).AND.(.NOT.MLPVB3T).AND.(.NOT.ORBITALS)415:      &         (.NOT.DIFFRACTT).AND.(.NOT.GAUSST).AND.(.NOT.PERCOLATET).AND.(.NOT.MLPVB3T).AND.(.NOT.ORBITALS)
416:      &        .AND.(.NOT.NORESET).AND.(.NOT.PERIODIC).AND.(.NOT.THOMSONT).AND.(.NOT.ONEDAPBCT).AND.(.NOT.ONEDPBCT)416:      &        .AND.(.NOT.NORESET).AND.(.NOT.PERIODIC).AND.(.NOT.THOMSONT).AND.(.NOT.ONEDAPBCT).AND.(.NOT.ONEDPBCT)
417:      &      .AND.(.NOT.TWODPBCT).AND.(.NOT.THREEDAPBCT).AND.(.NOT.THREEDPBCT).AND.(.NOT.QCIPOTT).AND.(.NOT.MLP3T)417:      &      .AND.(.NOT.TWODPBCT).AND.(.NOT.THREEDAPBCT).AND.(.NOT.THREEDPBCT).AND.(.NOT.QCIPOTT).AND.(.NOT.MLP3T).AND.(.NOT.MLPB3T)
418:      &      .AND.(.NOT.MLPB3T)418:      &       .AND.(.NOT.MLQT).AND.(.NOT.TWODAPBCT).AND.(.NOT.((NCORE(NP).GT.0).AND.(J1.GT.NATOMS-NCORE(NP))))) THEN
419:      &       .AND.(.NOT.MLQT).AND.(.NOT.TWODAPBCT).AND.(.NOT.((NCORE(NP).GT.0).AND.(J1.GT.NATOMS-NCORE(NP)))) 
420:      &       .AND.(.NOT.RIGIDINIT)) THEN 
421: 419: 
422:             IF (DEBUG) WRITE(MYUNIT,'(A,I4,A,F12.4,A,F12.4,A,I4,A,F12.4)') 'angular move for atom ',J1, 420:             IF (DEBUG) WRITE(MYUNIT,'(A,I4,A,F12.4,A,F12.4,A,I4,A,F12.4)') 'angular move for atom ',J1, 
423:      &           ' V=',VMAX,' Vmin=',VMIN,' next most weakly bound atom is ',JMAX2,' V=',VMAX2421:      &           ' V=',VMAX,' Vmin=',VMIN,' next most weakly bound atom is ',JMAX2,' V=',VMAX2
424: 422: 
425:            THETA=DPRAND()*PI423:            THETA=DPRAND()*PI
426:            PHI=DPRAND()*PI*2.0D0424:            PHI=DPRAND()*PI*2.0D0
427: !425: !
428: !  Evaporation is judged from the origin, not the centre of mass. We don't want the426: !  Evaporation is judged from the origin, not the centre of mass. We don't want the
429: !  angular move to cause evaporation. Obviously this will cause problems if we have a cluster that drifts427: !  angular move to cause evaporation. Obviously this will cause problems if we have a cluster that drifts
430: !  away from the origin up to the container radius.  428: !  away from the origin up to the container radius.  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0