hdiff output

r31741/bulkmindist.f90 2017-01-21 10:35:23.126250000 +0000 r31740/bulkmindist.f90 2017-01-21 10:35:23.586250000 +0000
  9: !  9: !
 10: !   OPTIM is distributed in the hope that it will be useful, 10: !   OPTIM is distributed in the hope that it will be useful,
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19:  
 20: ! Atom-matching method for aligning periodic structures in cartesian coordinates, and identifying permutational isomers. 
 21: ! Input: DUMMYB and DUMMYA are the two structures to be aligned 
 22: ! Output:  
 23: !         DUMMYB will be left unchanged. If a permutational isomer is identified, it is returned in DUMMYA and PITEST=TRUE 
 24: !         Otherwise, DUMMYA is left unchanged and PITEST=FALSE 
 25: !         If BMTEST is true, we are using the atom-matching method to find the best alignment of the two structures. In this 
 26: !         case XBEST will contain structure A with the best-identified translation and permutation applied. 
 27: SUBROUTINE BULKMINDIST(DUMMYB,DUMMYA,XBEST, NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,RESETA, TNMATCH, BMTEST) 19: SUBROUTINE BULKMINDIST(DUMMYB,DUMMYA,XBEST, NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,RESETA, TNMATCH, BMTEST)
 28: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, GEOMDIFFTOL, ATOMMATCHFULL, BULKBOXT, GDSQ, GDSQT 20: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, GEOMDIFFTOL, ATOMMATCHFULL, BULKBOXT, GDSQ, GDSQT
 29:  21: 
 30: IMPLICIT NONE 22: IMPLICIT NONE
 31: INTEGER J1, NATOMS, NPMIN, NGMIN, J2, PERM(NATOMS), PBEST(NATOMS), NDUMMY, NMATCHED, PATOMS, J3, J4, NMBEST, ND1, ND2 23: INTEGER J1, NATOMS, NPMIN, NGMIN, J2, PERM(NATOMS), PBEST(NATOMS), NDUMMY, NMATCHED, PATOMS, J3, J4, NMBEST, ND1
 32: INTEGER NMATCHSV, J5, J6, LPERM(NATOMS), NREP, NREP2, J7, J8, DEGS(3*NATOMS), TEMP_PERM(NATOMS) 24: INTEGER NMATCHSV, J5, J6, LPERM(NATOMS), NREP, NREP2
 33: DOUBLE PRECISION DUMMYB(3*NATOMS),DUMMYA(3*NATOMS),DISTANCE,BOXLX,BOXLY,BOXLZ,XSHIFT,YSHIFT,ZSHIFT,XTEMP(3*NATOMS) 25: DOUBLE PRECISION DUMMYB(3*NATOMS),DUMMYA(3*NATOMS),DISTANCE,BOXLX,BOXLY,BOXLZ,XSHIFT,YSHIFT,ZSHIFT,XTEMP(3*NATOMS)
 34: DOUBLE PRECISION XBEST(3*NATOMS), DMIN, DTOTAL, DIST, TEMP_DIST 26: DOUBLE PRECISION XBEST(3*NATOMS), DMIN, DTOTAL, DIST
 35: DOUBLE PRECISION DIST2, LDISTANCE, WORSTRAD , DUMMYE, VNEW(3*NATOMS), RMS 27: DOUBLE PRECISION DIST2, LDISTANCE, WORSTRAD 
 36: LOGICAL TWOD,DEBUG,PITEST,SAMEMIN,RESETA, TNMATCH, BMTEST, BMTESTLOCAL  28: LOGICAL TWOD,DEBUG,PITEST,SAMEMIN,RESETA, TNMATCH, BMTEST, BMTESTLOCAL 
 37: COMMON /BULKSHIFT/ XSHIFT,YSHIFT,ZSHIFT 29: COMMON /BULKSHIFT/ XSHIFT,YSHIFT,ZSHIFT
 38: SAVE NMATCHSV 30: SAVE NMATCHSV
 39:  31: 
 40: IF(BMTEST .AND. ANY(NSETS(:).GT.0)) THEN 
 41:    WRITE(*,*) "bulkmindist> Error. Atom matching method does not currently support secondary sets of permutable atoms." 
 42:    STOP 
 43: ENDIF 
 44:  
 45: IF (.NOT.TNMATCH) NMATCHSV=0  32: IF (.NOT.TNMATCH) NMATCHSV=0 
 46: ! 33: !
 47: ! Find smallest group of permutable atoms. 34: ! Find smallest group of permutable atoms.
 48: ! Translate first atom of group to all positions and then find nearest atom within 35: ! Translate first atom of group to all positions and then find nearest atom within
 49: ! the same group for every other atom. 36: ! the same group for every other atom.
 50: ! Keep the best translation/permutation, which corresponds to the smallest 37: ! Keep the best translation/permutation, which corresponds to the smallest
 51: ! minimum image distance. 38: ! minimum image distance.
 52: ! 39: !
 53: !PRINT*, NMATCHSV 40: !PRINT*, NMATCHSV
 54: TNMATCH=.TRUE. 41: TNMATCH=.TRUE.
 55: DISTANCE=1.0D100 42: DISTANCE=1.0D100
 56: PITEST=.FALSE. 43: PITEST=.FALSE.
 57: SAMEMIN=.TRUE. 44: SAMEMIN=.TRUE.
 58: BMTESTLOCAL=BMTEST 45: BMTESTLOCAL=BMTEST
 59: IF (.NOT.GDSQT) THEN 46: IF (.NOT.GDSQT) THEN
 60:    GDSQ=GEOMDIFFTOL**2/NATOMS ! because GEOMDIFFTOL is used for the total distance elsewhere in the program 47:    GDSQ=GEOMDIFFTOL**2/NATOMS ! because GEOMDIFFTOL is used for the total distance elsewhere in the program
 61: ENDIF 48: ENDIF
 62: NPMIN=HUGE(1) 49: NPMIN=HUGE(1)
 63:  50: ! PRINT *,'DUMMYA in bulkmindist:'
  51: ! PRINT '(3G20.10)',DUMMYA(1:3*NATOMS)
  52: ! PRINT *,'DUMMYB in bulkmindist:'
  53: ! PRINT '(3G20.10)',DUMMYB(1:3*NATOMS)
 64: DO J1=1,NPERMGROUP 54: DO J1=1,NPERMGROUP
 65:    IF (NPERMSIZE(J1).LT.NPMIN) THEN 55:    IF (NPERMSIZE(J1).LT.NPMIN) THEN
 66:       NPMIN=NPERMSIZE(J1) 56:       NPMIN=NPERMSIZE(J1)
 67:       NGMIN=J1 57:       NGMIN=J1
 68:    ENDIF 58:    ENDIF
 69: ENDDO 59: ENDDO
 70: ND1=0 60: ND1=0
 71: ND2=1 
 72: DO J1=1,NGMIN-1 61: DO J1=1,NGMIN-1
 73:    ND1=ND1+NPERMSIZE(J1)  ! ND1 is an offset index for the PERMGROUP array, indexing the start of the smallest permutational group. 62:    ND1=ND1+NPERMSIZE(J1)
 74: ENDDO 63: ENDDO
 75:  IF (DEBUG) PRINT '(3(A,I6))',' bulkmindist> Smallest group of permutable atoms is number ',NGMIN,' with ',NPMIN,' members' 64: ! IF (DEBUG) PRINT '(3(A,I6))',' bulkmindist> Smallest group of permutable atoms is number ',NGMIN,' with ',NPMIN,' members'
 76:  
 77: NREP=0 65: NREP=0
 78: NREP2=0 66: NREP2=0
 79: outer: DO J1=ND1+1,ND1+NPMIN  ! Loop over the smallest permutable group 67: outer: DO J1=ND1+1,ND1+NPMIN
 80:    ! Consider renaming J2 in the next line - it's confusing to use the same name for a loop index and a normal variable. 68:    J2=PERMGROUP(J1)
 81:    J2=PERMGROUP(J1)  ! Get the atom index for this cycle (PERMGROUP contains the indices of the atoms in all permutable groups, in group order) 
 82:    IF(DEBUG) WRITE(*,*) "Mapping atom ", J2, "in DUMMYA onto atom ", PERMGROUP(ND1+1), "in DUMMYB " 
 83:    ! Map the current atom (J2) onto the first atom in this group (PERMGROUP(ND1+1). First, calculate the (negative of the) required shift... 
 84:    XSHIFT=DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+1)-BOXLX*NINT((DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+1))/BOXLX) 69:    XSHIFT=DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+1)-BOXLX*NINT((DUMMYA(3*(J2-1)+1)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+1))/BOXLX)
 85:    YSHIFT=DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+2)-BOXLY*NINT((DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+2))/BOXLY) 70:    YSHIFT=DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+2)-BOXLY*NINT((DUMMYA(3*(J2-1)+2)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+2))/BOXLY)
 86:    IF (.NOT.TWOD) ZSHIFT=DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+3)-BOXLZ*NINT((DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+3))/BOXLZ) 71:    IF (.NOT.TWOD) ZSHIFT=DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+3)-BOXLZ*NINT((DUMMYA(3*(J2-1)+3)-DUMMYB(3*(PERMGROUP(ND1+1)-1)+3))/BOXLZ)
 87:    ! ...and then apply this shift to all atoms in DUMMYA. 
 88:    DO J2=1,NATOMS 72:    DO J2=1,NATOMS
 89:       XTEMP(3*(J2-1)+1)=DUMMYA(3*(J2-1)+1)-XSHIFT 73:       XTEMP(3*(J2-1)+1)=DUMMYA(3*(J2-1)+1)-XSHIFT
 90:       XTEMP(3*(J2-1)+2)=DUMMYA(3*(J2-1)+2)-YSHIFT 74:       XTEMP(3*(J2-1)+2)=DUMMYA(3*(J2-1)+2)-YSHIFT
 91:       IF (.NOT.TWOD) XTEMP(3*(J2-1)+3)=DUMMYA(3*(J2-1)+3)-ZSHIFT 75:       IF (.NOT.TWOD) XTEMP(3*(J2-1)+3)=DUMMYA(3*(J2-1)+3)-ZSHIFT
 92:    ENDDO 76:    ENDDO
 93:  
 94:    ! Now, we will find the best matches between the remaining atoms in the shifted DUMMYA and the unshifted DUMMYB. If all atoms match closely, we 
 95:    ! have a translation-permutation isomer. Otherwise, we continue superimposing atoms to find the translation which gives the greatest number 
 96:    ! of exact matches of atoms. This is returned as XBEST. 
 97:    NDUMMY=1 77:    NDUMMY=1
 98:    PERM(1:NATOMS)=-1 78:    PERM(1:NATOMS)=-1
 99:    NMATCHED=0 79:    NMATCHED=0
100:    DTOTAL=0.0D0 80:    DTOTAL=0.0D0
101:    permgroups: DO J2=1,NPERMGROUP  ! For each permgroup (not just the smallest)  WARNING: J2 is redefined here! 81:    DO J2=1,NPERMGROUP
102:       IF(J2.GT.1) NDUMMY=NDUMMY+NPERMSIZE(J2-1) 82:       IF(J2.GT.1) NDUMMY=NDUMMY+NPERMSIZE(J2-1)
103:       PATOMS=NPERMSIZE(J2) 83:       PATOMS=NPERMSIZE(J2)
104:       loop1: DO J3=1,PATOMS    ! for each atom in fixed structure B in group J2 84:       loop1: DO J3=1,PATOMS    ! for each atom in fixed structure B in group J2
105:          DMIN=1.0D100 85:          DMIN=1.0D100
106:          loop2: DO J4=1,PATOMS ! which is the closest atom in the same group for the structure in XTEMP (shifted A)? 86:          loop2: DO J4=1,PATOMS ! which is the closest atom in the same group for the structure in XTEMP (shifted A)?
107:             DIST=(XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+1)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+1) & 87:             DIST=(XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+1)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+1) &
108:  &  - BOXLX*NINT((XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+1)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+1))/BOXLX))**2 & 88:  &  - BOXLX*NINT((XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+1)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+1))/BOXLX))**2 &
109:             &  + (XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+2)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+2) & 89:             &  + (XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+2)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+2) &
110:  &  - BOXLY*NINT((XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+2)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+2))/BOXLY))**2  90:  &  - BOXLY*NINT((XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+2)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+2))/BOXLY))**2 
111:             IF (.NOT.TWOD) DIST=DIST+(XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+3)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+3) & 91:             IF (.NOT.TWOD) DIST=DIST+(XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+3)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+3) &
112:  &  - BOXLZ*NINT((XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+3)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+3))/BOXLZ))**2 92:  &  - BOXLZ*NINT((XTEMP(3*(PERMGROUP(NDUMMY+J4-1)-1)+3)-DUMMYB(3*(PERMGROUP(NDUMMY+J3-1)-1)+3))/BOXLZ))**2
113:             IF (DIST.LT.DMIN) THEN  ! Found a new atom which is the best match in the current cycle 93:             IF (DIST.LT.DMIN) THEN
114:                DMIN=DIST 94:                DMIN=DIST
115:                PERM(PERMGROUP(NDUMMY+J3-1))=PERMGROUP(NDUMMY+J4-1) 95:                PERM(PERMGROUP(NDUMMY+J3-1))=PERMGROUP(NDUMMY+J4-1)
116:                IF (DIST.LT.GDSQ) THEN  ! Is it a genuine match as defined by GEOMDIFFTOL? (See above) 96:                IF (DIST.LT.GDSQ) THEN
117:                  PRINT '(A,I6,A,I6,A,G20.10)',' match found between atom ',PERMGROUP(NDUMMY+J3-1), & 97: !                 PRINT '(A,I6,A,I6,A,G20.10)',' match found between atom ',PERMGROUP(NDUMMY+J3-1), &
118:  &                                            ' and ',PERMGROUP(NDUMMY+J4-1),' DIST=',DIST 98: ! &                                            ' and ',PERMGROUP(NDUMMY+J4-1),' DIST=',DIST
119:                   NMATCHED=NMATCHED+1 99:                   NMATCHED=NMATCHED+1
120:                   DTOTAL=DTOTAL+DMIN100:                   DTOTAL=DTOTAL+DMIN
121:                   IF (PERM(PERMGROUP(NDUMMY+J3-1)).NE.PERMGROUP(NDUMMY+J3-1)) SAMEMIN=.FALSE.101:                   IF (PERM(PERMGROUP(NDUMMY+J3-1)).NE.PERMGROUP(NDUMMY+J3-1)) SAMEMIN=.FALSE.
122:                   CYCLE loop1  ! We've matched this atom in the unshifted structure, so move on to the next.102:                   CYCLE loop1
123:                ENDIF103:                ENDIF
124:             ENDIF104:             ENDIF
125:          ENDDO loop2  ! We've tried matching atom J3 with every permutationally related atom in the shifted structure.105:          ENDDO loop2
126: !       DTOTAL=DTOTAL+DMIN106: !       DTOTAL=DTOTAL+DMIN
127:        ! If we've got this far, it means we haven't hit the CYCLE loop1 statement, and hence we did not find a match. 
128: !       PRINT '(A,I6,A,G20.10,A,I6)',' match failed for atom ',PERMGROUP(NDUMMY+J3-1),' DMIN=',DMIN,' J1=',J1107: !       PRINT '(A,I6,A,G20.10,A,I6)',' match failed for atom ',PERMGROUP(NDUMMY+J3-1),' DMIN=',DMIN,' J1=',J1
129: 108:       IF (.NOT.BMTESTLOCAL) CYCLE outer ! If we reached here then we don't have a permutational isomer because
130:       ! If we reached here then the atom specified in the J3 loop does not have a partner, so this choice of atoms to superimpose does not109:                       ! the atom specified in the J3 loop does not have a partner.
131:       ! give us exactly-aligned structures. 
132:       ! If BMTESTLOCAL is FALSE, we are not trying to find the best alignment of the two structures by atom matching (either because  
133:       ! ATOMMATCHDIST is not set, or has already failed). In this case, we are only interested in whether the two structures are  
134:       ! permutational isomers so we can skip the next section. Instead, we try the next pair of atoms to superimpose. 
135:       IF (.NOT.BMTESTLOCAL) CYCLE outer  
136:  
137:       IF ((J3-NMATCHED).GT.(NATOMS-NMATCHSV)) CYCLE outer ! Match cannot be better than the previous best110:       IF ((J3-NMATCHED).GT.(NATOMS-NMATCHSV)) CYCLE outer ! Match cannot be better than the previous best
138:       IF ((.NOT.ATOMMATCHFULL).AND.NMATCHSV.GT.0.AND.(J3-NMATCHED).GT.INT(NATOMS/2.0D0)) THEN  ! Conditions to give up on ATOMMATCHDIST111:       IF ((.NOT.ATOMMATCHFULL).AND.NMATCHSV.GT.0.AND.(J3-NMATCHED).GT.INT(NATOMS/2.0D0)) THEN
139:           NREP2=NREP2+1112:           NREP2=NREP2+1
140:           IF (NREP2.GE.5) THEN  113:           IF (NREP2.GE.5) THEN  
141:             BMTESTLOCAL=.FALSE.114:             BMTESTLOCAL=.FALSE.
142:             CYCLE outer ! Match is not good enough - give up115:             CYCLE outer ! Match is not good enough - give up
143:           ENDIF 116:           ENDIF 
144:       ELSE117:       ELSE
145:           NREP2=0118:           NREP2=0
146:       ENDIF119:       ENDIF
147:       ENDDO loop1 ! We have finished determining the number of matches associated with this choice of atoms to superimpose120:       ENDDO loop1 
148:       IF (DEBUG) PRINT '(A, I6, A, I6)',' bulkmindist> number of matching atoms=', NMATCHED, ' in permgroups ', J2   121:       IF (DEBUG) PRINT '(A, I6, A, I6)',' bulkmindist> number of matching atoms=', NMATCHED, ' in permgroups ', J2   
149:       IF (BMTESTLOCAL.AND.((NATOMS-NMATCHED).GT.0)) THEN  ! We are still looking for better matches (i.e. translation vectors which match more122:       IF (BMTESTLOCAL.AND.((NATOMS-NMATCHED).GT.0)) THEN
150:                                                           ! atoms). There are still atoms left to match.123:        IF (J2.EQ.NPERMGROUP.AND.NMATCHED.GE.NMATCHSV) THEN ! We have cycled over all atoms. Record best match. 
151:          IF (J2.EQ.NPERMGROUP.AND.NMATCHED.GE.NMATCHSV) THEN ! We have checked all permutational groups for matches, and found an improvement124:         CALL MINPERM(NATOMS, DUMMYB, XTEMP, BOXLX, BOXLY, BOXLZ, .TRUE., LPERM, LDISTANCE, DIST2, WORSTRAD)
152:                                                              ! compared to the previous best.125:         IF (LDISTANCE.LT.DISTANCE) THEN
153: 126:           NREP=0
154:             ! See whether a permutational alignment improves the match and record the new best match.127:           PRINT*, 'From minperm', LDISTANCE, DIST2
155: 128:           DISTANCE=LDISTANCE
156:             ! This obsolete call to MINPERM will find the best permutational alignment of ALL atoms, rather than each of the129:           NMATCHSV=NMATCHED
157:             ! permutational groups.130:           IF (RESETA) THEN
158: !            CALL MINPERM(NATOMS, DUMMYB, XTEMP, BOXLX, BOXLY, BOXLZ, .TRUE., LPERM, LDISTANCE, DIST2, WORSTRAD)131:            DO J6=1,NATOMS
159: 132:             XBEST(3*(J6-1)+1)=XTEMP(3*(LPERM(J6)-1)+1)-BOXLX*NINT(XTEMP(3*(LPERM(J6)-1)+1)/BOXLX)
160:             ! Initialise some variables for the permutational alignment133:             XBEST(3*(J6-1)+2)=XTEMP(3*(LPERM(J6)-1)+2)-BOXLY*NINT(XTEMP(3*(LPERM(J6)-1)+2)/BOXLY)
161:             LDISTANCE = 0.0D0134:             IF (.NOT.TWOD) XBEST(3*(J6-1)+3)=XTEMP(3*(LPERM(J6)-1)+3)-BOXLZ*NINT(XTEMP(3*(LPERM(J6)-1)+3)/BOXLZ)
162:             TEMP_DIST = 0.0D0135:            ENDDO
163:             DO J7=1, NATOMS136:           ENDIF
164:                LPERM(J7) = J7137:         ELSE IF (NMATCHED.EQ.NMATCHSV) THEN
165:             ENDDO138:           NREP=NREP+1
166:             DO J7=1,NPERMGROUP 139:           IF (NREP.GT.10) THEN
167:                ! ND2 is the index of the start of this permutation group.140:             IF (.NOT.ATOMMATCHFULL) BMTESTLOCAL=.FALSE. 
168:                ! PERMGROUP(ND2:ND2+NPERMSIZE(J7)-1) is an array containing the indices of all atoms belonging to this group.141:             CYCLE outer ! Match is always the same - give up, still do the PI test
169:                ! We build up a new array DEGS(:NPERMSIZE(J7)) which contains the indices of the corresponding degrees of freedom142:           ENDIF
170:                ! (i.e. indices for the coordinate arrays that correspond to this permutational group)143:         ENDIF
171:                DO J8=1,NPERMSIZE(J7)144:        ELSE IF (J2.EQ.NPERMGROUP) THEN
172:                   DEGS(3*J8-2) = 3*(PERMGROUP(ND2+J8-1))-2145:         NREP=0 
173:                   DEGS(3*J8-1) = 3*(PERMGROUP(ND2+J8-1))-1146:        ENDIF  
174:                   DEGS(3*J8) = 3*(PERMGROUP(ND2+J8-1))147:        IF (J2.EQ.NPERMGROUP.AND.NMATCHED.NE.NATOMS) CYCLE outer 
175:                ENDDO 
176:   
177:                ! Now do the actual permutational alignment 
178:                ! DUMMYB(DEGS(1:NPERMSIZE(J7))) extracts from DUMMYB the coordinates of this group only. 
179:                CALL MINPERM(NPERMSIZE(J7), DUMMYB(DEGS(1:3*NPERMSIZE(J7))), XTEMP(DEGS(1:3*NPERMSIZE(J7))), & 
180: &                           BOXLX, BOXLY, BOXLZ, .TRUE., TEMP_PERM(:NPERMSIZE(J7)), TEMP_DIST, DIST2, WORSTRAD) 
181:                DO J8=1,NPERMSIZE(J7) 
182:                   LPERM(PERMGROUP(ND2+J8-1)) = PERMGROUP(ND2+TEMP_PERM(J8)-1) 
183:                ENDDO 
184:                LDISTANCE = LDISTANCE + TEMP_DIST 
185:                ND2 = ND2+NPERMSIZE(J7) 
186:             ENDDO 
187:             ND2 = 1 ! Reset for next time we go through this block 
188:  
189:             IF(DEBUG) THEN 
190:                TEMP_DIST = 0.0D0 
191:                DO J7=1,NATOMS      
192:                   TEMP_DIST=TEMP_DIST+ (XTEMP(3*LPERM(J7)-2)-DUMMYB(3*J7-2) - BOXLX*NINT((XTEMP(3*LPERM(J7)-2)-DUMMYB(3*J7-2))/BOXLX))**2 & 
193: &                                    + (XTEMP(3*LPERM(J7)-1)-DUMMYB(3*J7-1) - BOXLY*NINT((XTEMP(3*LPERM(J7)-1)-DUMMYB(3*J7-1))/BOXLY))**2  
194:                   IF (.NOT.TWOD) TEMP_DIST=TEMP_DIST+(XTEMP(3*LPERM(J7))-DUMMYB(3*J7) - BOXLZ*NINT((XTEMP(3*LPERM(J7))-DUMMYB(3*J7))/BOXLZ))**2 
195:                ENDDO 
196:                IF(ABS(TEMP_DIST-LDISTANCE)/LDISTANCE.GT.1.0D-4) THEN 
197:                   WRITE(*,*) "bulkmindist> Warning: distance from MINPERM and recalculated distance do not agree:" 
198:                   WRITE(*,*) LDISTANCE, TEMP_DIST 
199:                ENDIF 
200:             ENDIF 
201:  
202:             IF (LDISTANCE.LT.DISTANCE) THEN  ! Performing the alignment improved the distance. Update the best match. 
203:                NREP=0 
204:                ! We used to print DIST2 as well, but it is now meaningless (it only refers to one of the permutational groups) 
205:                IF(DEBUG) PRINT*, 'bulkmindist> Distance from minperm', LDISTANCE!, DIST2 
206:                DISTANCE=LDISTANCE 
207:                NMATCHSV=NMATCHED 
208:                IF (RESETA) THEN  ! (I think this is always TRUE for this subroutine) 
209:                   ! XBEST contains the best-so-far aligned structure (translation+permutation) 
210:                   DO J6=1,NATOMS 
211:                      ! XTEMP is the most recent shifted structure, LPERM is the best permutation as obtained by MINPERM 
212:                      XBEST(3*(J6-1)+1)=XTEMP(3*(LPERM(J6)-1)+1)-BOXLX*NINT(XTEMP(3*(LPERM(J6)-1)+1)/BOXLX) 
213:                      XBEST(3*(J6-1)+2)=XTEMP(3*(LPERM(J6)-1)+2)-BOXLY*NINT(XTEMP(3*(LPERM(J6)-1)+2)/BOXLY) 
214:                      IF (.NOT.TWOD) XBEST(3*(J6-1)+3)=XTEMP(3*(LPERM(J6)-1)+3)-BOXLZ*NINT(XTEMP(3*(LPERM(J6)-1)+3)/BOXLZ) 
215:                   ENDDO 
216:                ENDIF 
217:  
218:             ELSE IF (NMATCHED.EQ.NMATCHSV) THEN  ! We have been through a whole cycle of outer (i.e. tried a different superposition) 
219:                                                  ! without any change in the number of matches 
220:                NREP=NREP+1  ! Increment the fail counter. Give up when it passes 10, unless ATOMMATCHFULL is set. 
221:                IF (NREP.GT.10) THEN 
222:                   IF (.NOT.ATOMMATCHFULL) BMTESTLOCAL=.FALSE.  
223:                   write(*,*) "Giving up now; ATOMMATCHDIST just set BMTESTLOCAL to FALSE" 
224:                   CYCLE outer ! Match is always the same - give up, still do the PI test 
225:                ENDIF 
226:             ENDIF 
227:          ELSE IF (J2.EQ.NPERMGROUP) THEN  ! we have reached the end of this permutational group, no improvement on previous best 
228:             NREP=0  
229:          ENDIF   
230:          IF (J2.EQ.NPERMGROUP.AND.NMATCHED.NE.NATOMS) CYCLE outer ! We have finished this group and not yet found a full match. 
231:                                                                   ! Move on to the next  
232:       ENDIF  148:       ENDIF  
233:    ENDDO  permgroups149:    ENDDO
234:    IF (SQRT(DTOTAL).GE.GEOMDIFFTOL) CYCLE outer  ! We don't have an overall permutational isomer. Try next superposition.150:    IF (SQRT(DTOTAL).GE.GEOMDIFFTOL) CYCLE outer
235:    IF (SAMEMIN) THEN  ! no permutation was necessary to identify the isomers: the two structures are translationally related.151:    IF (SAMEMIN) THEN
236:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> identical isomers identified for distance ',SQRT(DTOTAL)152:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> identical isomers identified for distance ',SQRT(DTOTAL)
237:    ELSE153:    ELSE
238:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> permutational isomers identified for distance ',SQRT(DTOTAL)154:       IF (DEBUG) PRINT '(A,G20.10)',' bulkmindist> permutational isomers identified for distance ',SQRT(DTOTAL)
239:    ENDIF155:    ENDIF
240:    PITEST=.TRUE.   ! We have a permutational isomer. This variable gets used by MINPERMDIST on return.156:    PITEST=.TRUE.
241:    DISTANCE=DTOTAL ! Update the distance (should be 0)157:    DISTANCE=DTOTAL
242:    ! Now if RESETA is TRUE (which it always is), update DUMMYA to contain the coordinates of the isomer 
243:    IF (RESETA .AND. BULKBOXT) THEN158:    IF (RESETA .AND. BULKBOXT) THEN
244:       ! Don't put the coordinates back in the box first 
245:       DO J2=1,NATOMS159:       DO J2=1,NATOMS
246:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1)160:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1)
247:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2)161:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2)
248:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3)162:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3)
249:       ENDDO163:       ENDDO
250:    ELSE IF (RESETA) THEN164:    ELSE IF (RESETA) THEN
251:       ! Do put the coordinates back in the box. 
252:       DO J2=1,NATOMS165:       DO J2=1,NATOMS
253:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1)-BOXLX*NINT(XTEMP(3*(PERM(J2)-1)+1)/BOXLX)166:          DUMMYA(3*(J2-1)+1)=XTEMP(3*(PERM(J2)-1)+1)-BOXLX*NINT(XTEMP(3*(PERM(J2)-1)+1)/BOXLX)
254:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2)-BOXLY*NINT(XTEMP(3*(PERM(J2)-1)+2)/BOXLY)167:          DUMMYA(3*(J2-1)+2)=XTEMP(3*(PERM(J2)-1)+2)-BOXLY*NINT(XTEMP(3*(PERM(J2)-1)+2)/BOXLY)
255:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3)-BOXLZ*NINT(XTEMP(3*(PERM(J2)-1)+3)/BOXLZ)168:          IF (.NOT.TWOD) DUMMYA(3*(J2-1)+3)=XTEMP(3*(PERM(J2)-1)+3)-BOXLZ*NINT(XTEMP(3*(PERM(J2)-1)+3)/BOXLZ)
256:       ENDDO169:       ENDDO
257:    ENDIF170:    ENDIF
258: 171: 
259:    RETURN ! Don't need to do any more if we found an isomer!172:    RETURN
260: ENDDO outer  173: ENDDO outer
261: 174: 
262: ! If we got this far, we tried all the possible superpositions and didn't find a PI isomer. XBEST should contain the best translation-permutational 
263: ! alignment of DUMMYA that we found. DUMMYB is unchanged, as usual. DUMMYA is also unchanged, in this case. 
264: IF (DEBUG) PRINT '(A, G20.10)',' bulkmindist> distance=', DISTANCE175: IF (DEBUG) PRINT '(A, G20.10)',' bulkmindist> distance=', DISTANCE
265: IF (DEBUG) PRINT '(A)',' bulkmindist> structures are not permutational isomers'176: IF (DEBUG) PRINT '(A)',' bulkmindist> structures are not permutational isomers'
266: 177: 
267: RETURN178: RETURN
268: 179: 
269: END SUBROUTINE BULKMINDIST180: END SUBROUTINE BULKMINDIST
270: !181: !
271: ! Apply Oh point group operation number OPNUM to coordinates in182: ! Apply Oh point group operation number OPNUM to coordinates in
272: ! vector X of dimension 3*NLOCAL, returning the result in 183: ! vector X of dimension 3*NLOCAL, returning the result in 
273: ! vector Y.184: ! vector Y.


r31741/minpermdist.f90 2017-01-21 10:35:23.350250000 +0000 r31740/minpermdist.f90 2017-01-21 10:35:23.810250000 +0000
456: ! right orientation! 456: ! right orientation! 
457: !457: !
458: ! If we use MYORIENT to produce particular orientations then we end up aligning 458: ! If we use MYORIENT to produce particular orientations then we end up aligning 
459: ! COORDSA not with COORDSB but with the standard orientation of COORDSB in DUMMYB.459: ! COORDSA not with COORDSB but with the standard orientation of COORDSB in DUMMYB.
460: ! We now deal with this by tracking the complete transformation, including the460: ! We now deal with this by tracking the complete transformation, including the
461: ! contribution of MYORIENT using ROTB and ROTINVB.461: ! contribution of MYORIENT using ROTB and ROTINVB.
462: !462: !
463: 463: 
464: DISTANCE=0.0D0464: DISTANCE=0.0D0
465: IF (NFREEZE.LE.0 .AND. .NOT.MIEFT .AND. .NOT. NOTRANSROTT) THEN465: IF (NFREEZE.LE.0 .AND. .NOT.MIEFT .AND. .NOT. NOTRANSROTT) THEN
466:    IF (BULKT) THEN  ! we will not be doing any rotations.466:    IF (BULKT) THEN
467:       IF (BMTEST) BMCOORDSSV(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)467:       IF (BMTEST) BMCOORDSSV(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)
468:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;468:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;
469:       ! BULKMINDIST is used for a quick deterministic check for permutational isomers for periodic systems.  
470:       ! If BMTEST is TRUE (i.e. ATOMMATCHDIST or ATOMMATCHFULL) then this routine also attempts to find the translational-permutational  
471:       ! alignment that maximises the number of exact atom matches between the two structures. This is not necessarily the alignment with  
472:       ! the minimum cartesian distance. 
473:       ! BULKMINDIST returns DUMMYB unchanged and BMCOORDS as the best trans+perm alignment of structure A. If a permutational isomer is 
474:       ! detected, DUMMYA contains its coordinates and PITEST is TRUE. Otherwise, PITEST is FALSE and DUMMYA should be unchanged. 
475:       write(*,*) "Before call to BULKMINDIST" 
476:       CALL BULKMINDIST(DUMMYB,DUMMYA,BMCOORDS,NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,.TRUE., TNMATCH, BMTEST)469:       CALL BULKMINDIST(DUMMYB,DUMMYA,BMCOORDS,NATOMS,DISTANCE,TWOD,DEBUG,BOXLX,BOXLY,BOXLZ,PITEST,.TRUE., TNMATCH, BMTEST)
477:       IF (PITEST) THEN470:       IF (PITEST) THEN
478:          ! The two structures are permutational isomers. The aligned structure A is contained in DUMMYA 
479:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)471:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
480:          ! No rotations (periodic system) 
481:          RMATBEST(1:3,1:3)=0.0D0472:          RMATBEST(1:3,1:3)=0.0D0
482:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0473:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
483:          DISTANCE=SQRT(DISTANCE)474:          DISTANCE=SQRT(DISTANCE)
484:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST, distance=',DISTANCE475:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST, distance=',DISTANCE
485:          IF (DEBUG) PRINT '(A)',' minpermdist> permutation-inversion isomers identified by BULKMINDIST' 
486:          ! As required, DUMMYA contains the best alignment of structure A and DISTANCE contains the optimal distance (should be approx. 0) 
487:          ! We can now return from MINPERMDIST. 
488:          RETURN476:          RETURN
489:       ELSE IF (BMTEST) THEN477:       ELSE IF (BMTEST) THEN
490:          ! We are using an atom-matching method to get the best trans-perm alignment. 
491:          ! If OHCELLT is TRUE, we will go through this block 48 times, once for each symmetry operation of the Oh group - because of the 
492:          ! GOTO statement below. 
493:          ! If OHCELLT is FALSE, we should only go through it once. 
494:  
495:          IF (SQRT(DISTANCE).LT.BMDIST) THEN478:          IF (SQRT(DISTANCE).LT.BMDIST) THEN
496:             ! If OHCELLT is FALSE, we should always end up here (since BMDIST is initialised to a HUGE value and we only go through once). 
497:             BMDIST=SQRT(DISTANCE)479:             BMDIST=SQRT(DISTANCE)
498:             IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST distance=', BMDIST480:             IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST distance=', BMDIST
499:          ELSE481:          ELSE
500:             ! This should only ever happen if applying an Oh symmetry operation makes the alignment less good than the previous pass. 
501:             BMCOORDS(1:3*NATOMS)=BMCOORDSSV(1:3*NATOMS)482:             BMCOORDS(1:3*NATOMS)=BMCOORDSSV(1:3*NATOMS)
502:             IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST distance has not improved'483:             IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULKMINDIST distance has not improved'
503:          ENDIF484:          ENDIF
504: !         IF ((OHCELLT.AND.OPNUM.EQ.0).OR.(.NOT.OHCELLT)) THEN485: !         IF ((OHCELLT.AND.OPNUM.EQ.0).OR.(.NOT.OHCELLT)) THEN
505: !           CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX',.FALSE.,RIGID,DEBUG,RMAT)486: !           CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX',.FALSE.,RIGID,DEBUG,RMAT)
506: !         IF (BMDIST.GT.DISTANCE) THEN487: !         IF (BMDIST.GT.DISTANCE) THEN
507: !! Only reach here if atom-matching is really awful488: !! Only reach here if atom-matching is really awful
508: !! No bipartite matching etc.. for distance489: !! No bipartite matching etc.. for distance
509: !            IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> shortest distance so far from NEWMINDIST distance=', DISTANCE490: !            IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> shortest distance so far from NEWMINDIST distance=', DISTANCE
510: !            PRINT '(A,G20.10)',' minpermdist> WARNING, atom-matching aborted, NEWMINDIST distance=', DISTANCE491: !            PRINT '(A,G20.10)',' minpermdist> WARNING, atom-matching aborted, NEWMINDIST distance=', DISTANCE
511: !          BMTEST=.FALSE.492: !          BMTEST=.FALSE.
512: !!         OPNUM=48 493: !!         OPNUM=48 
513: !         ENDIF494: !         ENDIF
514: !         ENDIF495: !         ENDIF
515:  
516:          ! The following statment should always be TRUE. Looking at the preceding IF block, then either BMDIST.GE.SQRT(DISTANCE),  
517:          ! in which case we set BMDIST=SQRT(DISTANCE) so BMDIST.LT.DISTANCE, or BMDIST.LT.SQRT(DISTANCE) in which case  
518:          ! BMDIST.LT.DISTANCE is implied - unless DISTANCE<1. I'm not sure what happens in that case! This may be a bug. 
519:          IF (BMDIST.LT.DISTANCE) THEN496:          IF (BMDIST.LT.DISTANCE) THEN
520:             DISTANCE=BMDIST497:             DISTANCE=BMDIST
521:             IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> shortest distance so far from BULKMINDIST distance=', DISTANCE498:             IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> shortest distance so far from BULKMINDIST distance=', DISTANCE
522:             IF (OHCELLT.AND.(OPNUM.LT.48)) THEN499:             IF (OHCELLT.AND.(OPNUM.LT.48)) THEN
523:              ! Apply the next symmetry operation from the Oh point group (which is appropriate for cubic simulation boxes) 
524:              ! and try BULKMINDIST again. 
525:              GOTO 25500:              GOTO 25
526:             ELSE501:             ELSE
527:              ! If we're ignoring the Oh symmetry operations, or if we've checked all of them, then we have obtained the best 
528:              ! answer possible from this atom matching method. So we set DUMMYA equal to BMCOORDS (the best tr-perm alignment  
529:              ! of structure A) and recalculate the distance from B to compare with the distance obtained by BULKMINDIST. 
530:              ! Then we return from MINPERMDIST. 
531:              DUMMYA(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)502:              DUMMYA(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)
532:              COORDSA(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)503:              COORDSA(1:3*NATOMS)=BMCOORDS(1:3*NATOMS)
533:              DISTANCE=0.0D0504:              DISTANCE=0.0D0
534:              DO J3=1,NATOMS505:              DO J3=1,NATOMS
535:                   DISTANCE=DISTANCE + (COORDSA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1) &506:                   DISTANCE=DISTANCE + (COORDSA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1) &
536:   &                       -BOXLX*NINT((COORDSA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1))/BOXLX))**2 &507:   &                       -BOXLX*NINT((COORDSA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1))/BOXLX))**2 &
537:   &                                 + (COORDSA(3*(J3-1)+2)-DUMMYB(3*(J3-1)+2) &508:   &                                 + (COORDSA(3*(J3-1)+2)-DUMMYB(3*(J3-1)+2) &
538:   &                       -BOXLY*NINT((COORDSA(3*(J3-1)+2)-DUMMYB(3*(J3-1)+2))/BOXLY))**2509:   &                       -BOXLY*NINT((COORDSA(3*(J3-1)+2)-DUMMYB(3*(J3-1)+2))/BOXLY))**2
539:                   IF (.NOT.TWOD) DISTANCE=DISTANCE+(COORDSA(3*(J3-1)+3)-DUMMYB(3*(J3-1)+3)&510:                   IF (.NOT.TWOD) DISTANCE=DISTANCE+(COORDSA(3*(J3-1)+3)-DUMMYB(3*(J3-1)+3)&
540:   &                       -BOXLZ*NINT((COORDSA(3*(J3-1)+3)-DUMMYB(3*(J3-1)+3))/BOXLZ))**2511:   &                       -BOXLZ*NINT((COORDSA(3*(J3-1)+3)-DUMMYB(3*(J3-1)+3))/BOXLZ))**2
541:              ENDDO512:              ENDDO
542:              DISTANCE=SQRT(DISTANCE)513:              DISTANCE=SQRT(DISTANCE)
543:              IF (DEBUG) PRINT*, ' minpermdist> Recalculated distance=', DISTANCE514:              IF (DEBUG) PRINT*, ' minpermdist> Recalculated distance=', DISTANCE
544:  
545:              RMATBEST(1:3,1:3)=0.0D0515:              RMATBEST(1:3,1:3)=0.0D0
546:              RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0516:              RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
547:  
548:  
549:              IF (DEBUG) THEN 
550:                 IF (CHRMMT) CALL UPDATENBONDS(COORDSA)  ! But we probably shouldn't have BULKT and CHRMMT anyway 
551:                 IF (RIGIDINIT) THEN 
552:                    CALL GENRIGID_POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
553:                 ELSE 
554:                    CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
555:                 ENDIF 
556:                 PRINT '(2(A,F25.15))',' final   energy for structure A=             ',ENERGY,' RMS=',RMS 
557:                 IF (ABS(ENERGY-AINIT).GT.2*EDIFFTOL) THEN 
558:                    PRINT '(A)',' minpermdist> ERROR *** energy change for structure A is outside tolerance' 
559:                    STOP 
560:                 ENDIF 
561:                 IF (CHRMMT) CALL UPDATENBONDS(COORDSB) 
562:                 IF (RIGIDINIT) THEN 
563:                    CALL GENRIGID_POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
564:                 ELSE 
565:                    CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
566:                 ENDIF 
567:                 PRINT '(2(A,F25.15))',' final   energy for structure B=             ',ENERGY,' RMS=',RMS 
568:                 IF (ABS(ENERGY-BINIT).GT.2*EDIFFTOL) THEN 
569:                    PRINT '(A)',' minpermdist> ERROR *** energy change for structure B is outside tolerance' 
570:                    STOP 
571:                 ENDIF 
572:              ENDIF 
573:  
574:              RETURN517:              RETURN
575:             ENDIF518:             ENDIF
576:          ELSE 
577:             WRITE(*,*) "minpermdist> Error in BMTEST (atom-matching) block. Stopping now" 
578:             STOP 
579:          ENDIF519:          ENDIF
580:       ELSE520:       ELSE
581:          CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)521:          CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
582:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULK/NEWMINDIST distance=',DISTANCE522:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULK/NEWMINDIST distance=',DISTANCE
583:          DISTANCE=DISTANCE**2 ! minperdist returns the distance squared for historical reasons523:          DISTANCE=DISTANCE**2 ! minperdist returns the distance squared for historical reasons
584:       ENDIF524:       ENDIF
585:    ELSEIF (MKTRAPT) THEN525:    ELSEIF (MKTRAPT) THEN
586:       TMAT(1:3,1:3)=0.0D0526:       TMAT(1:3,1:3)=0.0D0
587:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0527:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0
588:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;528:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;
700: !  The maximum number of pair exchanges associated with a group is two.640: !  The maximum number of pair exchanges associated with a group is two.
701: !641: !
702: NTRIES=0642: NTRIES=0
703: !643: !
704: !  RMATCUMUL contains the accumulated rotation matrix that relates the original 644: !  RMATCUMUL contains the accumulated rotation matrix that relates the original 
705: !  DUMMYA obtained from COORDSA to the final one.645: !  DUMMYA obtained from COORDSA to the final one.
706: !646: !
707: RMATCUMUL(1:3,1:3)=0.0D0647: RMATCUMUL(1:3,1:3)=0.0D0
708: RMATCUMUL(1,1)=1.0D0; RMATCUMUL(2,2)=1.0D0; RMATCUMUL(3,3)=1.0D0648: RMATCUMUL(1,1)=1.0D0; RMATCUMUL(2,2)=1.0D0; RMATCUMUL(3,3)=1.0D0
709: 10 CONTINUE649: 10 CONTINUE
710:  
711: NTRIES=NTRIES+1650: NTRIES=NTRIES+1
712: 651: 
713: NDUMMY=1652: NDUMMY=1
714: DO J1=1,NATOMS653: DO J1=1,NATOMS
715:    NEWPERM(J1)=J1654:    NEWPERM(J1)=J1
716: ENDDO655: ENDDO
717: !656: !
718: ! ALLPERM saves the permutation from the previous cycle.657: ! ALLPERM saves the permutation from the previous cycle.
719: ! NEWPERM contains the permutation for this cycle, relative to the identity.658: ! NEWPERM contains the permutation for this cycle, relative to the identity.
720: ! SAVEPERM is temporary storage for NEWPERM.659: ! SAVEPERM is temporary storage for NEWPERM.


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0