hdiff output

r30496/keywords.f 2016-05-22 11:30:08.298418357 +0100 r30495/keywords.f 2016-05-22 11:30:08.682423433 +0100
1733:             IF (NITEMS.GT.2) CALL READI(MAXNSETS)1733:             IF (NITEMS.GT.2) CALL READI(MAXNSETS)
1734:             PRINT '(A,F15.5)','keywords> Distance tolerance for identifying atoms in the same orbit=',ORBITTOL1734:             PRINT '(A,F15.5)','keywords> Distance tolerance for identifying atoms in the same orbit=',ORBITTOL
1735:             PRINT '(A,I3)',' keywords> Maximum number of secondary sets in perm.allow file=', MAXNSETS1735:             PRINT '(A,I3)',' keywords> Maximum number of secondary sets in perm.allow file=', MAXNSETS
1736:          ENDIF1736:          ENDIF
1737:          IF (WORD.EQ.'PERMISOMER') PERMISOMER=.TRUE.1737:          IF (WORD.EQ.'PERMISOMER') PERMISOMER=.TRUE.
1738:          IF (WORD.EQ.'LPERMDIST') THEN1738:          IF (WORD.EQ.'LPERMDIST') THEN
1739:             LPERMDIST=.TRUE.1739:             LPERMDIST=.TRUE.
1740:             IF (NITEMS.GT.1) CALL READI(LOCALPERMNEIGH)1740:             IF (NITEMS.GT.1) CALL READI(LOCALPERMNEIGH)
1741:             IF (NITEMS.GT.2) CALL READF(LOCALPERMCUT)1741:             IF (NITEMS.GT.2) CALL READF(LOCALPERMCUT)
1742:             IF (NITEMS.GT.3) CALL READF(LOCALPERMCUT2)1742:             IF (NITEMS.GT.3) CALL READF(LOCALPERMCUT2)
1743:             IF (NITEMS.GT.4) CALL READF(ORBITTOL) 
1744:             PRINT '(A,F15.5)','keywords> Local permutational alignment: alignment threshold=',LOCALPERMCUT1743:             PRINT '(A,F15.5)','keywords> Local permutational alignment: alignment threshold=',LOCALPERMCUT
1745:             PRINT '(A,F15.5)','keywords> Local permutational alignment: alignment cutoff=   ',LOCALPERMCUT21744:             PRINT '(A,F15.5)','keywords> Local permutational alignment: alignment cutoff=   ',LOCALPERMCUT2
1746:             PRINT '(A,F15.5)',' keyword> Distance tolerance for distinguishing atoms in the same orbit=',ORBITTOL 
1747:          ENDIF1745:          ENDIF
1748:          INQUIRE(FILE='perm.allow',EXIST=PERMFILE)1746:          INQUIRE(FILE='perm.allow',EXIST=PERMFILE)
1749:          ALLOCATE(NPERMSIZE(3*NATOMS),PERMGROUP(3*NATOMS),NSETS(3*NATOMS),SETS(NATOMS,MAXNSETS))1747:          ALLOCATE(NPERMSIZE(3*NATOMS),PERMGROUP(3*NATOMS),NSETS(3*NATOMS),SETS(NATOMS,MAXNSETS))
1750: !1748: !
1751: !  The above dimensions were fixed at NATOMS because:1749: !  The above dimensions were fixed at NATOMS because:
1752: !  (a) Atoms were not allowed to appear in more than one group.1750: !  (a) Atoms were not allowed to appear in more than one group.
1753: !  (b) The maximum number of pair exchanges associated with a group is three.1751: !  (b) The maximum number of pair exchanges associated with a group is three.
1754: !1752: !
1755: ! However, for flexible water models we need to exchange all waters,1753: ! However, for flexible water models we need to exchange all waters,
1756: ! and we can exchange H's attached to the same O. The dimension required1754: ! and we can exchange H's attached to the same O. The dimension required


r30496/lopermdist.f90 2016-05-22 11:30:08.486420839 +0100 r30495/lopermdist.f90 2016-05-22 11:30:08.894426228 +0100
 73: !  and is evaluated as a sum of local distances squared for permutable groups. 73: !  and is evaluated as a sum of local distances squared for permutable groups.
 74: !  We return to label 10 after every round of permutational/orientational alignment 74: !  We return to label 10 after every round of permutational/orientational alignment
 75: !  unless we have converged to the identity permutation. 75: !  unless we have converged to the identity permutation.
 76: ! 76: !
 77: !  The maximum number of pair exchanges associated with a group is two. 77: !  The maximum number of pair exchanges associated with a group is two.
 78: !  78: ! 
 79: DO J1=1,NATOMS 79: DO J1=1,NATOMS
 80:    NEWPERM(J1)=J1 80:    NEWPERM(J1)=J1
 81: ENDDO 81: ENDDO
 82: DSUM=0.0D0 82: DSUM=0.0D0
 83: LOCALPERMNEIGH=MIN(LOCALPERMNEIGH,NATOMS) 
 84:  83: 
 85: NDUMMY=1 84: NDUMMY=1
 86: DO J1=1,NPERMGROUP 85: DO J1=1,NPERMGROUP
 87:    PATOMS=NPERMSIZE(J1) 86:    PATOMS=NPERMSIZE(J1)
 88: !  PRINT '(A,I6,A,I6)','group ',J1,' size=',NPERMSIZE(J1) 
 89: !  PRINT '(A)','members:' 
 90: !  DO J2=1,PATOMS 
 91: !     PRINT '(20I6)',NEWPERM(PERMGROUP(NDUMMY+J2-1)) 
 92: !  ENDDO 
 93:    LDBEST(J1)=1.0D100 87:    LDBEST(J1)=1.0D100
 94:    TRIED(1:NATOMS)=0 88:    TRIED(1:NATOMS)=0
 95:    DO J2=1,PATOMS 89:    DO J2=1,PATOMS
 96:       LPERMBEST(J2)=J2 90:       LPERMBEST(J2)=J2
 97:    ENDDO 91:    ENDDO
 98:    XA=0.0D0; YA=0.0D0; ZA=0.0D0 92:    XA=0.0D0; YA=0.0D0; ZA=0.0D0
 99:    XB=0.0D0; YB=0.0D0; ZB=0.0D0 93:    XB=0.0D0; YB=0.0D0; ZB=0.0D0
100:    DMEAN(1:LOCALPERMNEIGH)=1.0D100 94:    DMEAN(1:LOCALPERMNEIGH)=1.0D100
101:    DO J2=1,PATOMS 95:    DO J2=1,PATOMS
102: !     TRIED(NEWPERM(PERMGROUP(NDUMMY+J2-1)))=-1 96:       TRIED(NEWPERM(PERMGROUP(NDUMMY+J2-1)))=-1
103:       TRIED(PERMGROUP(NDUMMY+J2-1))=-1 
104:       PDUMMYA(3*(J2-1)+1)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+1) 97:       PDUMMYA(3*(J2-1)+1)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+1)
105:       PDUMMYA(3*(J2-1)+2)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+2) 98:       PDUMMYA(3*(J2-1)+2)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+2)
106:       PDUMMYA(3*(J2-1)+3)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+3) 99:       PDUMMYA(3*(J2-1)+3)=DUMMYA(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+3)
107: !     PDUMMYB(3*(J2-1)+1)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+1)100:       PDUMMYB(3*(J2-1)+1)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+1)
108: !     PDUMMYB(3*(J2-1)+2)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+2)101:       PDUMMYB(3*(J2-1)+2)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+2)
109: !     PDUMMYB(3*(J2-1)+3)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+3)102:       PDUMMYB(3*(J2-1)+3)=DUMMYB(3*(NEWPERM(PERMGROUP(NDUMMY+J2-1))-1)+3)
110:       PDUMMYB(3*(J2-1)+1)=DUMMYB(3*(PERMGROUP(NDUMMY+J2-1)-1)+1) 
111:       PDUMMYB(3*(J2-1)+2)=DUMMYB(3*(PERMGROUP(NDUMMY+J2-1)-1)+2) 
112:       PDUMMYB(3*(J2-1)+3)=DUMMYB(3*(PERMGROUP(NDUMMY+J2-1)-1)+3) 
113:       XA=XA+PDUMMYA(3*(J2-1)+1)103:       XA=XA+PDUMMYA(3*(J2-1)+1)
114:       YA=YA+PDUMMYA(3*(J2-1)+2)104:       YA=YA+PDUMMYA(3*(J2-1)+2)
115:       ZA=ZA+PDUMMYA(3*(J2-1)+3)105:       ZA=ZA+PDUMMYA(3*(J2-1)+3)
116:       XB=XB+PDUMMYB(3*(J2-1)+1)106:       XB=XB+PDUMMYB(3*(J2-1)+1)
117:       YB=YB+PDUMMYB(3*(J2-1)+2)107:       YB=YB+PDUMMYB(3*(J2-1)+2)
118:       ZB=ZB+PDUMMYB(3*(J2-1)+3)108:       ZB=ZB+PDUMMYB(3*(J2-1)+3)
119:    ENDDO109:    ENDDO
120:    XA=XA/PATOMS; YA=YA/PATOMS; ZA=ZA/PATOMS110:    XA=XA/PATOMS; YA=YA/PATOMS; ZA=ZA/PATOMS
121:    XB=XB/PATOMS; YB=YB/PATOMS; ZB=ZB/PATOMS111:    XB=XB/PATOMS; YB=YB/PATOMS; ZB=ZB/PATOMS
122:    SPDUMMYA(1:3*PATOMS)=PDUMMYA(1:3*PATOMS)112:    SPDUMMYA(1:3*PATOMS)=PDUMMYA(1:3*PATOMS)
140: !130: !
141: ! Don't allow members of the same permutational group 131: ! Don't allow members of the same permutational group 
142: ! to appear as reference neighbours.132: ! to appear as reference neighbours.
143: !133: !
144:       IF (TRIED(J2).EQ.-1) THEN134:       IF (TRIED(J2).EQ.-1) THEN
145:          XDUMMY=1.0D9135:          XDUMMY=1.0D9
146:       ELSE136:       ELSE
147:          DA=(XA-DUMMYA(3*(NEWPERM(J2)-1)+1))**2 &137:          DA=(XA-DUMMYA(3*(NEWPERM(J2)-1)+1))**2 &
148:   &        +(YA-DUMMYA(3*(NEWPERM(J2)-1)+2))**2 &138:   &        +(YA-DUMMYA(3*(NEWPERM(J2)-1)+2))**2 &
149:   &        +(ZA-DUMMYA(3*(NEWPERM(J2)-1)+3))**2139:   &        +(ZA-DUMMYA(3*(NEWPERM(J2)-1)+3))**2
150: !        DB=(XB-DUMMYB(3*(NEWPERM(J2)-1)+1))**2 &140:          DB=(XB-DUMMYB(3*(NEWPERM(J2)-1)+1))**2 &
151: ! &        +(YB-DUMMYB(3*(NEWPERM(J2)-1)+2))**2 &141:   &        +(YB-DUMMYB(3*(NEWPERM(J2)-1)+2))**2 &
152: ! &        +(ZB-DUMMYB(3*(NEWPERM(J2)-1)+3))**2142:   &        +(ZB-DUMMYB(3*(NEWPERM(J2)-1)+3))**2
153:          DB=(XB-DUMMYB(3*(J2-1)+1))**2 & 
154:   &        +(YB-DUMMYB(3*(J2-1)+2))**2 & 
155:   &        +(ZB-DUMMYB(3*(J2-1)+3))**2 
156:          XDUMMY=(SQRT(DA)+SQRT(DB))/2.0D0143:          XDUMMY=(SQRT(DA)+SQRT(DB))/2.0D0
157:       ENDIF144:       ENDIF
158:       loop1: DO J3=1,J2145:       loop1: DO J3=1,J2
159:          IF (XDUMMY.LT.DMEAN(J3)) THEN146:          IF (XDUMMY.LT.DMEAN(J3)) THEN
160: !147: !
161: ! Move the rest down.148: ! Move the rest down.
162: !149: !
163:             DO J4=J2,J3+1,-1150:             DO J4=J2,J3+1,-1
164:                DMEAN(J4)=DMEAN(J4-1)151:                DMEAN(J4)=DMEAN(J4-1)
165:                SORTLIST(J4)=SORTLIST(J4-1)152:                SORTLIST(J4)=SORTLIST(J4-1)
166:             ENDDO153:             ENDDO
167:             DMEAN(J3)=XDUMMY154:             DMEAN(J3)=XDUMMY
168:             SORTLIST(J3)=J2155:             SORTLIST(J3)=J2
169:             EXIT loop1156:             EXIT loop1
170:          ENDIF157:          ENDIF
171:       ENDDO loop1158:       ENDDO loop1
172:    ENDDO outer1159:    ENDDO outer1
 160: !  IF (J1.EQ.16) THEN
 161: !     PRINT '(A)','SORTLIST:'
 162: !     PRINT '(20I5)',SORTLIST(1:NATOMS)
 163: !     PRINT '(A)','DMEAN:'
 164: !     PRINT '(10G10.4)',DMEAN(1:NATOMS)
 165: !  ENDIF
173: 166: 
174: 71 CONTINUE167: 71 CONTINUE
175:    PDUMMYA(1:3*PATOMS)=SPDUMMYA(1:3*PATOMS)168:    PDUMMYA(1:3*PATOMS)=SPDUMMYA(1:3*PATOMS)
176:    PDUMMYB(1:3*PATOMS)=SPDUMMYB(1:3*PATOMS)169:    PDUMMYB(1:3*PATOMS)=SPDUMMYB(1:3*PATOMS)
177: 170: 
178:    LDBESTATOM=1.0D100171:    LDBESTATOM=1.0D100
179:    NOTHER=0172:    NOTHER=0
180:    DO J2=1,NATOMS173:    DO J2=1,NATOMS
181:       IF (TRIED(J2).EQ.1) THEN174:       IF (TRIED(J2).EQ.1) THEN
182:          NOTHER=NOTHER+1175:          NOTHER=NOTHER+1
226:             PRINT '(A,I6,A)',' lopermdist> ERROR *** Partner atom ',DLIST(NOTHER),' has already been tried'219:             PRINT '(A,I6,A)',' lopermdist> ERROR *** Partner atom ',DLIST(NOTHER),' has already been tried'
227:             STOP220:             STOP
228:          ENDIF221:          ENDIF
229:       ENDDO222:       ENDDO
230:    ENDIF223:    ENDIF
231:    224:    
232:    DO J2=1,NOTHER225:    DO J2=1,NOTHER
233:       PDUMMYA(3*(PATOMS+J2-1)+1)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+1)226:       PDUMMYA(3*(PATOMS+J2-1)+1)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+1)
234:       PDUMMYA(3*(PATOMS+J2-1)+2)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+2)227:       PDUMMYA(3*(PATOMS+J2-1)+2)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+2)
235:       PDUMMYA(3*(PATOMS+J2-1)+3)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+3)228:       PDUMMYA(3*(PATOMS+J2-1)+3)=DUMMYA(3*(NEWPERM(DLIST(J2))-1)+3)
236: !     PDUMMYB(3*(PATOMS+J2-1)+1)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+1)229:       PDUMMYB(3*(PATOMS+J2-1)+1)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+1)
237: !     PDUMMYB(3*(PATOMS+J2-1)+2)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+2)230:       PDUMMYB(3*(PATOMS+J2-1)+2)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+2)
238: !     PDUMMYB(3*(PATOMS+J2-1)+3)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+3)231:       PDUMMYB(3*(PATOMS+J2-1)+3)=DUMMYB(3*(NEWPERM(DLIST(J2))-1)+3)
239:       PDUMMYB(3*(PATOMS+J2-1)+1)=DUMMYB(3*(DLIST(J2)-1)+1) 
240:       PDUMMYB(3*(PATOMS+J2-1)+2)=DUMMYB(3*(DLIST(J2)-1)+2) 
241:       PDUMMYB(3*(PATOMS+J2-1)+3)=DUMMYB(3*(DLIST(J2)-1)+3) 
242:    ENDDO232:    ENDDO
 233: !  IF ((J1.EQ.16).OR.(J1.EQ.12)) THEN
 234: !  PRINT '(4(A,I6))',' lopermdist> For group ',J1,' size ',PATOMS,' aligning with ',NOTHER,' other atoms'
 235: !  PRINT '(A)',' DLIST:'
 236: !  PRINT '(20I6)',DLIST(1:NOTHER)
 237: !  ENDIF
243: !238: !
244: ! Save PDUMMYA and PDUMMYB for cycling over possible orbits in MYORIENT alignment.239: ! Save PDUMMYA and PDUMMYB for cycling over possible orbits in MYORIENT alignment.
245: !240: !
246:    SPDUMMYA(3*PATOMS+1:3*(PATOMS+NOTHER))=PDUMMYA(3*PATOMS+1:3*(PATOMS+NOTHER))241:    SPDUMMYA(3*PATOMS+1:3*(PATOMS+NOTHER))=PDUMMYA(3*PATOMS+1:3*(PATOMS+NOTHER))
247:    SPDUMMYB(3*PATOMS+1:3*(PATOMS+NOTHER))=PDUMMYB(3*PATOMS+1:3*(PATOMS+NOTHER))242:    SPDUMMYB(3*PATOMS+1:3*(PATOMS+NOTHER))=PDUMMYB(3*PATOMS+1:3*(PATOMS+NOTHER))
248:    NCHOOSEB1=0243:    NCHOOSEB1=0
249: 66 NCHOOSEB1=NCHOOSEB1+1244: 66 NCHOOSEB1=NCHOOSEB1+1
250:    NCHOOSEB2=0245:    NCHOOSEB2=0
251: 31 NCHOOSEB2=NCHOOSEB2+1246: 31 NCHOOSEB2=NCHOOSEB2+1
252:    NCHOOSE1=0247:    NCHOOSE1=0
268: ! Optimimise permutational isomer for the standard orientation for the263: ! Optimimise permutational isomer for the standard orientation for the
269: ! current choice of atoms from the possible orbits.264: ! current choice of atoms from the possible orbits.
270: !265: !
271: ! MINPERM does not change PDUMMYB and PDUMMYA.266: ! MINPERM does not change PDUMMYB and PDUMMYA.
272: !267: !
273: ! Note that LDISTANCE is actually the distance squared. LDBEST also has dimensions of268: ! Note that LDISTANCE is actually the distance squared. LDBEST also has dimensions of
274: ! length squared.269: ! length squared.
275: !270: !
276:    LDISTANCE=0.0D0271:    LDISTANCE=0.0D0
277:    CALL MINPERM(PATOMS+NOTHER, PDUMMYB, PDUMMYA, BOXLX, BOXLY, BOXLZ, BULKT, LPERM, LDISTANCE, DIST2, WORSTRAD) 272:    CALL MINPERM(PATOMS+NOTHER, PDUMMYB, PDUMMYA, BOXLX, BOXLY, BOXLZ, BULKT, LPERM, LDISTANCE, DIST2, WORSTRAD) 
278: !  PRINT '(A,I6,A,I6)','for group ',J1,' size=',NPERMSIZE(J1)273: !  PRINT '(A,I6,3G20.10)','J1,LDBEST(J1),LDISTANCE=',J1,LDBEST(J1),LDISTANCE
279: !  PRINT '(A)','original and new atom labels:' 
280: !  DO J2=1,PATOMS 
281: !     PRINT '(2I6)',PERMGROUP(NDUMMY+J2-1),NEWPERM(PERMGROUP(NDUMMY+J2-1)) 
282: !  ENDDO 
283:  
284: !  PRINT '(I6,A)',NOTHER,' other atoms:' 
285: !  PRINT '(A)','original and new other atom labels:' 
286: !  DO J2=1,NOTHER 
287: !     PRINT '(2I6)',DLIST(J2),NEWPERM(DLIST(J2)) 
288: !  ENDDO 
289:  
290: !  PRINT '(A,3I6,3G20.10)','J1,PATOMS,NOTHER,LDBEST(J1),LDISTANCE=',J1,PATOMS,NOTHER,LDBEST(J1),LDISTANCE 
291: !  PRINT '(A,20I6)','LPERM after MINPERM: ',LPERM(1:PATOMS+NOTHER) 
292: 274: 
293:    LDISTANCE=LDISTANCE275:    LDISTANCE=LDISTANCE
294:    DO J2=1,PATOMS276:    DO J2=1,PATOMS
295:       IF (LPERM(J2).GT.PATOMS) THEN277:       IF (LPERM(J2).GT.PATOMS) THEN
296:          LDISTANCE=1.0D300278:          LDISTANCE=1.0D300
297: !        IF (DEBUG) PRINT '(A,I6,A,I6,A)',' lopermdist> For group ',J1,' with ',NOTHER,' neighbours - neighbours mix in' 279: !        IF (DEBUG) PRINT '(A,I6,A,I6,A)',' lopermdist> For group ',J1,' with ',NOTHER,' neighbours - neighbours mix in' 
298: !        IF (DEBUG) PRINT '(A,I6,A,I6)',' lopermdist> atom ',J2,' lperm value is ',LPERM(J2) 
299:          EXIT280:          EXIT
300:       ENDIF281:       ENDIF
301:    ENDDO282:    ENDDO
302: 283: 
 284: !  IF (J1.EQ.16) THEN
 285: !  PRINT '(I6)',PATOMS+NOTHER
 286: !  PRINT '(A,8I6,G20.10)',' PDUMMYB for NO1,NO2,NC1,NC2,NOB1,NOB2,NCB1,NCB2,distance ', &
 287: ! &                                     NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,NORBITB1,NORBITB2,NCHOOSEB1,NCHOOSEB2,SQRT(LDISTANCE)
 288: !  PRINT '(A,3G20.10)',('LA ',PDUMMYB(3*(J2-1)+1:3*(J2-1)+3),J2=1,PATOMS+NOTHER)
 289: !  PRINT '(I6)',PATOMS+NOTHER
 290: !  PRINT '(A,8I6,G20.10)',' PDUMMYA for NO1,NO2,NC1,NC2,NOB1,NOB2,NCB1,NCB2,distance ', &
 291: ! &                                     NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,NORBITB1,NORBITB2,NCHOOSEB1,NCHOOSEB2,SQRT(LDISTANCE)
 292: !  PRINT '(A,3G20.10)',('LA ',PDUMMYA(3*(J2-1)+1:3*(J2-1)+3),J2=1,PATOMS+NOTHER)
 293: !  ENDIF
 294: 
303:    DO J2=1,NOTHER295:    DO J2=1,NOTHER
304:       IF (LPERM(PATOMS+J2).NE.PATOMS+J2) THEN296:       IF (LPERM(PATOMS+J2).NE.PATOMS+J2) THEN
305: !        IF (DEBUG) PRINT '(A,I6,A,I6)',' lopermdist> Atom ',DLIST(J2),' also needs to permute to ',LPERM(PATOMS+J2)297: !        IF (DEBUG) PRINT '(A,I6,A)',' lopermdist> Atom ',DLIST(J2),' also needs to permute'
306:          IF (PERMUTABLE(DLIST(J2))) THEN298:          IF (PERMUTABLE(DLIST(J2))) THEN
307: !           IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(J2),' belongs to permutable set ', &299: !           IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(J2),' belongs to permutable set ', &
308: !  &                                INGROUP(DLIST(J2))300: !  &                                INGROUP(DLIST(J2))
309:          ELSE301:          ELSE
310: !           IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(J2),' is NOT permutable!'302: !           IF (DEBUG) PRINT '(2(A,I6))',' lopermdist> Atom ',DLIST(J2),' is NOT permutable!'
311:             LDISTANCE=1.0D300303:             LDISTANCE=1.0D300
312:          ENDIF304:          ENDIF
313:       ENDIF305:       ENDIF
314:    ENDDO306:    ENDDO
315: !307: !
340:       TRIED(DLIST(NOTHER))=-1332:       TRIED(DLIST(NOTHER))=-1
341:       IF (NADDED.GT.1) THEN333:       IF (NADDED.GT.1) THEN
342: !        IF (DEBUG) THEN334: !        IF (DEBUG) THEN
343: !           PRINT '(A)',' lopermdist> and partner atoms:'335: !           PRINT '(A)',' lopermdist> and partner atoms:'
344: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)336: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)
345: !        ENDIF337: !        ENDIF
346:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1338:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1
347:       ENDIF339:       ENDIF
348:       GOTO 71340:       GOTO 71
349:    ELSE341:    ELSE
350: !     IF (DEBUG) PRINT '(A,G20.10,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), &342: !     IF (DEBUG) PRINT '(A,F12.2,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), &
351: ! &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours' 343: ! &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours' 
352:       TRIED(DLIST(NOTHER))=1344:       TRIED(DLIST(NOTHER))=1
353:       IF (NADDED.GT.1) THEN345:       IF (NADDED.GT.1) THEN
354: !        IF (DEBUG) THEN346: !        IF (DEBUG) THEN
355: !           PRINT '(A)',' lopermdist> and partner atoms:'347: !           PRINT '(A)',' lopermdist> and partner atoms:'
356: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)348: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)
357: !        ENDIF349: !        ENDIF
358:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1350:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1
359:       ENDIF351:       ENDIF
360:       LDBEST(J1)=LDBESTATOM352:       LDBEST(J1)=LDBESTATOM
363: !     PRINT '(A,10I6)','LPERMBEST: ',LPERMBEST(1:PATOMS)355: !     PRINT '(A,10I6)','LPERMBEST: ',LPERMBEST(1:PATOMS)
364:    ENDIF356:    ENDIF
365: !357: !
366: ! Add the next eligible atom and try alignment again.358: ! Add the next eligible atom and try alignment again.
367: ! Stop if we already have LOCALPERMNEIGH neighbours.359: ! Stop if we already have LOCALPERMNEIGH neighbours.
368: !360: !
369:    IF (NOTHER.LT.LOCALPERMNEIGH) GOTO 71361:    IF (NOTHER.LT.LOCALPERMNEIGH) GOTO 71
370: 362: 
371: 91 CONTINUE ! jump here when there are no atoms left to try.363: 91 CONTINUE ! jump here when there are no atoms left to try.
372: 364: 
373: !  IF (DEBUG) PRINT '(2(A,I6),A,G15.5)',' lopermdist> For group ',J1,' maximum neighbours=', &365:    IF (DEBUG) PRINT '(2(A,I6),A,G15.5)',' lopermdist> For group ',J1,' maximum neighbours=', &
374: ! &                                      NOTHER,' distance=',SQRT(LDBEST(J1))366:   &                                      NOTHER,' distance=',SQRT(LDBEST(J1))
375: !367: !
376: ! We now have the best permutation for group J1 and standard orientations368: ! We now have the best permutation for group J1 and standard orientations
377: ! based upon all atoms belonging to the two possible orbits that appear369: ! based upon all atoms belonging to the two possible orbits that appear
378: ! for the standard alignment.370: ! for the standard alignment.
379: !371: !
380:    LPERM(1:PATOMS)=LPERMBEST(1:PATOMS)372:    LPERM(1:PATOMS)=LPERMBEST(1:PATOMS)
381: !373: !
382: ! Fill SAVEPERM with NEWPERM, which contains the current best permutation374: ! Fill SAVEPERM with NEWPERM, which contains the current best permutation
383: ! after the previous pass through J1375: ! after the previous pass through J1
384: !376: !
402:          DO J3=1,NSETS(J1)394:          DO J3=1,NSETS(J1)
403:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3)395:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3)
404:          ENDDO396:          ENDDO
405:       ENDDO397:       ENDDO
406:    ENDIF398:    ENDIF
407: !399: !
408: ! Save current optimal permutation in NEWPERM400: ! Save current optimal permutation in NEWPERM
409: !401: !
410:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS)402:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS)
411:    DSUM=DSUM+SQRT(LDBEST(J1))403:    DSUM=DSUM+SQRT(LDBEST(J1))
412: !  PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' best distance=',SQRT(LDBEST(J1)),' total=',DSUM404: !  PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' after myorient distance=',SQRT(LDBEST(J1)),' total=',DSUM
413: !  PRINT '(A)','best permutation is now' 
414: !  PRINT '(20I6)',NEWPERM(1:NATOMS) 
415: 405: 
416: !406: !
417: ! Update NDUMMY, the cumulative offset for PERMGROUP407: ! Update NDUMMY, the cumulative offset for PERMGROUP
418: !408: !
419:    NDUMMY=NDUMMY+NPERMSIZE(J1)409:    NDUMMY=NDUMMY+NPERMSIZE(J1)
420: ENDDO  !  end of loop over groups of permutable atoms410: ENDDO  !  end of loop over groups of permutable atoms
421: !411: !
422: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA412: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA
423: ! to the current best alignment. 413: ! to the current best alignment. 
424: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current414: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current
468: ! At this point NEWPERM, ALLPERM, SAVEPERM, BESTPERM458: ! At this point NEWPERM, ALLPERM, SAVEPERM, BESTPERM
469: ! are all the same!459: ! are all the same!
470: !460: !
471: ! PRINT '(A)',' lopermdist> NEWPERM, ALLPERM, SAVEPERM, BESTPERM:'461: ! PRINT '(A)',' lopermdist> NEWPERM, ALLPERM, SAVEPERM, BESTPERM:'
472: ! PRINT '(4I6)',(NEWPERM(J1),ALLPERM(J1),SAVEPERM(J1),BESTPERM(J1),J1=1,NATOMS)462: ! PRINT '(4I6)',(NEWPERM(J1),ALLPERM(J1),SAVEPERM(J1),BESTPERM(J1),J1=1,NATOMS)
473: 463: 
474: !!!!!!!!!!!!!!!!!!!!!!! DEBUG464: !!!!!!!!!!!!!!!!!!!!!!! DEBUG
475: !465: !
476: ! Test distance for COORDSA with permutation applied in BESTPERM466: ! Test distance for COORDSA with permutation applied in BESTPERM
477: !467: !
478: !  DO J1=1,NATOMS468:    DO J1=1,NATOMS
479: !     DUMMYA(3*(J1-1)+1)=COORDSA(3*(BESTPERM(J1)-1)+1)469:       DUMMYA(3*(J1-1)+1)=COORDSA(3*(BESTPERM(J1)-1)+1)
480: !     DUMMYA(3*(J1-1)+2)=COORDSA(3*(BESTPERM(J1)-1)+2)470:       DUMMYA(3*(J1-1)+2)=COORDSA(3*(BESTPERM(J1)-1)+2)
481: !     DUMMYA(3*(J1-1)+3)=COORDSA(3*(BESTPERM(J1)-1)+3)471:       DUMMYA(3*(J1-1)+3)=COORDSA(3*(BESTPERM(J1)-1)+3)
482: !  ENDDO472:    ENDDO
483: 473: 
484: !  CALL NEWMINDIST(COORDSB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)474:    CALL NEWMINDIST(COORDSB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
485: !  CALL NEWMINDIST(COORDSB,XBEST,NATOMS,XDUMMY,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)475:    CALL NEWMINDIST(COORDSB,XBEST,NATOMS,XDUMMY,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
486: !  IF (DEBUG) WRITE(*,'(A,2G20.10)') & 476:    IF (DEBUG) WRITE(*,'(A,2G20.10)') & 
487: ! &   ' lopermdist> distance check for permuted COORDSA and original COORDSB=',XDUMMY,DISTANCE477:   &   ' lopermdist> distance check for permuted COORDSA and original COORDSB=',XDUMMY,DISTANCE
488: !!!!!!!!!!!!!!!!!!!!!!! DEBUG478: !!!!!!!!!!!!!!!!!!!!!!! DEBUG
489: !479: !
490: ! Now align and reorient the permuted coordinates in COORDSA 480: ! Now align and reorient the permuted coordinates in COORDSA 
491: ! Try using the best locally aligned group of atoms481: ! Try using the best locally aligned group of atoms
492: !482: !
493: CALL NEWMINDIST(DUMMYB,XBEST,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)483: CALL NEWMINDIST(DUMMYB,XBEST,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
494: IF (DEBUG) PRINT '(A,G20.10)',' lopermdist> after overall alignment distance=',DISTANCE484: IF (DEBUG) PRINT '(A,G20.10)',' lopermdist> after overall alignment distance=',DISTANCE
495: RMATBEST(1:3,1:3)=RMAT(1:3,1:3)485: RMATBEST(1:3,1:3)=RMAT(1:3,1:3)
496: 486: 
497: COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!487: COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0