hdiff output

r29843/minpermdist.f90 2016-03-16 18:33:11.690835940 +0000 r29842/minpermdist.f90 2016-03-16 18:33:12.474844003 +0000
 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the
 51: !  centre of coordinates of COORDSA will be the same as for COORDSB. 51: !  centre of coordinates of COORDSA will be the same as for COORDSB.
 52: ! 52: !
 53: !     ---------------------------------------------------------------------------------------------- 53: !     ----------------------------------------------------------------------------------------------
 54: ! jdf43>        Modified for generalised angle-axis 30/01/12 54: ! jdf43>        Modified for generalised angle-axis 30/01/12
 55: !     ---------------------------------------------------------------------------------------------- 55: !     ----------------------------------------------------------------------------------------------
 56:  56: 
 57: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST) 57: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 58: USE COMMONS,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, CHRMMT, MYUNIT, STOCKT, NFREEZE, & 58: USE COMMONS,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, CHRMMT, MYUNIT, STOCKT, NFREEZE, &
 59:   & AMBERT, CSMT, PERMDIST, PULLT, EFIELDT, OHCELLT, NTSITES, GEOMDIFFTOL, QCIPERMCHECK, & 59:   & AMBERT, CSMT, PERMDIST, PULLT, EFIELDT, OHCELLT, NTSITES, GEOMDIFFTOL, QCIPERMCHECK, &
 60:   & PERMOPT, PERMINVOPT, NOINVERSION, BESTPERM, BESTINVERT, GTHOMSONT, LOCALPERMDIST,  LPERMDIST, MKTRAPT 60:   & PERMOPT, PERMINVOPT, NOINVERSION, BESTPERM, BESTINVERT, GTHOMSONT, LOCALPERMDIST,  LPERMDIST
 61: USE PORFUNCS 61: USE PORFUNCS
 62: USE GENRIGID 62: USE GENRIGID
 63: IMPLICIT NONE 63: IMPLICIT NONE
 64:  64: 
 65: INTEGER, PARAMETER :: MAXIMUMTRIES=20 65: INTEGER, PARAMETER :: MAXIMUMTRIES=20
 66: INTEGER NATOMS, NPERM, PATOMS, NTRIES, ISTAT, OPNUM, NCHOOSEB1, NCHOOSEB2, NORBITB1, NORBITB2, I 66: INTEGER NATOMS, NPERM, PATOMS, NTRIES, ISTAT, OPNUM, NCHOOSEB1, NCHOOSEB2, NORBITB1, NORBITB2, I
 67: INTEGER J3, INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1 67: INTEGER J3, INVERT, NORBIT1, NORBIT2, NCHOOSE2, NDUMMY, LPERM(NATOMS), J1, J2, NCHOOSE1
 68: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) 68: DOUBLE PRECISION DIST2, COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS)
 69: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS) 69: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,WORSTRAD,RMAT(3,3),ENERGY, VNEW(3*NATOMS), DX, DY, DZ, RMS, DBEST, XBEST(3*NATOMS)
 70: DOUBLE PRECISION CMXA, CMXB, CMXC, PDISTANCE, CMX, CMY, CMZ 70: DOUBLE PRECISION CMXA, CMXB, CMXC, PDISTANCE, CMX, CMY, CMZ
160: !  WRITE(10,'(I6)') NATOMS/2160: !  WRITE(10,'(I6)') NATOMS/2
161: !  WRITE(10,'(A)') 'B initial'161: !  WRITE(10,'(A)') 'B initial'
162: !  DO J3=1,NATOMS/2162: !  DO J3=1,NATOMS/2
163: !      WRITE(10,'(A2,2X,3F20.10)') 'LA',COORDSB(3*(J3-1)+1),COORDSB(3*(J3-1)+2),COORDSB(3*(J3-1)+3)163: !      WRITE(10,'(A2,2X,3F20.10)') 'LA',COORDSB(3*(J3-1)+1),COORDSB(3*(J3-1)+2),COORDSB(3*(J3-1)+3)
164: !  ENDDO164: !  ENDDO
165: !  CLOSE(10)165: !  CLOSE(10)
166: !166: !
167: !  Calculate original centres of mass.167: !  Calculate original centres of mass.
168: !168: !
169: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0169: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0
170: IF ((NFREEZE.GT.0).OR.(MKTRAPT)) GOTO 11 ! don;t shift or reorient coordinates with frozen atoms170: IF (NFREEZE.GT.0) GOTO 11 ! don;t shift or reorient coordinates with frozen atoms
171: IF (STOCKT.OR.RIGID) THEN 171: IF (STOCKT.OR.RIGID) THEN 
172:    DO J1=1,NATOMS/2172:    DO J1=1,NATOMS/2
173:       CMAX=CMAX+COORDSA(3*(J1-1)+1)173:       CMAX=CMAX+COORDSA(3*(J1-1)+1)
174:       CMAY=CMAY+COORDSA(3*(J1-1)+2)174:       CMAY=CMAY+COORDSA(3*(J1-1)+2)
175:       CMAZ=CMAZ+COORDSA(3*(J1-1)+3)175:       CMAZ=CMAZ+COORDSA(3*(J1-1)+3)
176:    ENDDO176:    ENDDO
177:    CMAX=2*CMAX/NATOMS; CMAY=2*CMAY/NATOMS; CMAZ=2*CMAZ/NATOMS177:    CMAX=2*CMAX/NATOMS; CMAY=2*CMAY/NATOMS; CMAZ=2*CMAZ/NATOMS
178:    CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0178:    CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0
179:    DO J1=1,NATOMS/2179:    DO J1=1,NATOMS/2
180:       CMBX=CMBX+COORDSB(3*(J1-1)+1)180:       CMBX=CMBX+COORDSB(3*(J1-1)+1)
190:    ENDDO190:    ENDDO
191:    CMAX=CMAX/NATOMS; CMAY=CMAY/NATOMS; CMAZ=CMAZ/NATOMS191:    CMAX=CMAX/NATOMS; CMAY=CMAY/NATOMS; CMAZ=CMAZ/NATOMS
192:    CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0192:    CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0
193:    DO J1=1,NATOMS193:    DO J1=1,NATOMS
194:       CMBX=CMBX+COORDSB(3*(J1-1)+1)194:       CMBX=CMBX+COORDSB(3*(J1-1)+1)
195:       CMBY=CMBY+COORDSB(3*(J1-1)+2)195:       CMBY=CMBY+COORDSB(3*(J1-1)+2)
196:       CMBZ=CMBZ+COORDSB(3*(J1-1)+3)196:       CMBZ=CMBZ+COORDSB(3*(J1-1)+3)
197:    ENDDO197:    ENDDO
198:    CMBX=CMBX/NATOMS; CMBY=CMBY/NATOMS; CMBZ=CMBZ/NATOMS198:    CMBX=CMBX/NATOMS; CMBY=CMBY/NATOMS; CMBZ=CMBZ/NATOMS
199: ENDIF199: ENDIF
200: 11 CONTINUE 
201: BESTINVERT=1200: BESTINVERT=1
202: DO J1=1,NATOMS201: DO J1=1,NATOMS
203:    BESTPERM(J1)=J1202:    BESTPERM(J1)=J1
204: ENDDO203: ENDDO
 204: 11 CONTINUE
205: 205: 
206: INVERT=1206: INVERT=1
207: DBEST=1.0D100207: DBEST=1.0D100
208: 60 CONTINUE ! jump back here if INVERT changes sign.208: 60 CONTINUE ! jump back here if INVERT changes sign.
209:    NCHOOSEB1=0209:    NCHOOSEB1=0
210: 66 NCHOOSEB1=NCHOOSEB1+1210: 66 NCHOOSEB1=NCHOOSEB1+1
211:    NCHOOSEB2=0211:    NCHOOSEB2=0
212: 31 NCHOOSEB2=NCHOOSEB2+1212: 31 NCHOOSEB2=NCHOOSEB2+1
213:    NCHOOSE1=0213:    NCHOOSE1=0
214: 65 NCHOOSE1=NCHOOSE1+1214: 65 NCHOOSE1=NCHOOSE1+1
258:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)258:          COORDSA(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
259:          RMATBEST(1:3,1:3)=0.0D0259:          RMATBEST(1:3,1:3)=0.0D0
260:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0260:          RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
261:          DISTANCE=SQRT(DISTANCE)261:          DISTANCE=SQRT(DISTANCE)
262:          RETURN262:          RETURN
263:       ELSE263:       ELSE
264:          CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)264:          CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
265:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'minpermdist> after initial call to BULK/NEWMINDIST distance=',DISTANCE265:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'minpermdist> after initial call to BULK/NEWMINDIST distance=',DISTANCE
266:          DISTANCE=DISTANCE**2 ! minperdist returns the distance squared for historical reasons266:          DISTANCE=DISTANCE**2 ! minperdist returns the distance squared for historical reasons
267:       ENDIF267:       ENDIF
268:    ELSEIF (MKTRAPT) THEN 
269:       TMAT(1:3,1:3)=0.0D0 
270:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0 
271:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1; 
272:       ROTB(1:3,1:3)=0.0D0 
273:       ROTB(1,1)=1.0D0; ROTB(2,2)=1.0D0; ROTB(3,3)=1.0D0 
274:       ROTINVB(1:3,1:3)=0.0D0 
275:       ROTINVB(1,1)=1.0D0; ROTINVB(2,2)=1.0D0; ROTINVB(3,3)=1.0D0 
276:       ROTA(1:3,1:3)=0.0D0 
277:       ROTA(1,1)=1.0D0; ROTA(2,2)=1.0D0; ROTA(3,3)=1.0D0 
278:       ROTINVA(1:3,1:3)=0.0D0 
279:       ROTINVA(1,1)=1.0D0; ROTINVA(2,2)=1.0D0; ROTINVA(3,3)=1.0D0 
280:       RMAT(1:3,1:3)=0.0D0 
281:       RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0 
282:       RMATBEST(1:3,1:3)=0.0D0 
283:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0 
284:       CMX=0.0D0; CMY=0.0D0; CMZ=0.0D0 
285:       DUMMYA(1:3*NATOMS)=DUMMYC(1:3*NATOMS) 
286:       DISTANCE=0.0D0 
287:       DO J1=1,3*NATOMS 
288:          DISTANCE=DISTANCE+(DUMMYA(J1)-DUMMYB(J1))**2 
289:       ENDDO 
290:    ELSEIF (STOCKT) THEN268:    ELSEIF (STOCKT) THEN
291:       TMAT(1:3,1:3)=0.0D0269:       TMAT(1:3,1:3)=0.0D0
292:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0270:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0
293:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMYA,TMAT,0.0D0,0.0D0,0.0D0)271:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMYA,TMAT,0.0D0,0.0D0,0.0D0)
294:       DUMMY(1:3*NATOMS)=DUMMYA(1:3*NATOMS)272:       DUMMY(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
295:       CALL MYORIENT(DUMMYA,DUMMYC,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,NATOMS/2,DEBUG,ROTA,ROTINVA,STOCKT)273:       CALL MYORIENT(DUMMYA,DUMMYC,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,NATOMS/2,DEBUG,ROTA,ROTINVA,STOCKT)
296:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMY,ROTA,0.0D0,0.0D0,0.0D0)274:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMY,ROTA,0.0D0,0.0D0,0.0D0)
297:       DUMMYA(1:3*NATOMS)=DUMMY(1:3*NATOMS)275:       DUMMYA(1:3*NATOMS)=DUMMY(1:3*NATOMS)
298: 276: 
299:       DUMMY(1:3*NATOMS)=DUMMYB(1:3*NATOMS)277:       DUMMY(1:3*NATOMS)=DUMMYB(1:3*NATOMS)
347:    ROTB(1,1)=1.0D0; ROTB(2,2)=1.0D0; ROTB(3,3)=1.0D0325:    ROTB(1,1)=1.0D0; ROTB(2,2)=1.0D0; ROTB(3,3)=1.0D0
348:    ROTINVB(1:3,1:3)=0.0D0326:    ROTINVB(1:3,1:3)=0.0D0
349:    ROTINVB(1,1)=1.0D0; ROTINVB(2,2)=1.0D0; ROTINVB(3,3)=1.0D0327:    ROTINVB(1,1)=1.0D0; ROTINVB(2,2)=1.0D0; ROTINVB(3,3)=1.0D0
350:    ROTA(1:3,1:3)=0.0D0328:    ROTA(1:3,1:3)=0.0D0
351:    ROTA(1,1)=1.0D0; ROTA(2,2)=1.0D0; ROTA(3,3)=1.0D0329:    ROTA(1,1)=1.0D0; ROTA(2,2)=1.0D0; ROTA(3,3)=1.0D0
352:    ROTINVA(1:3,1:3)=0.0D0330:    ROTINVA(1:3,1:3)=0.0D0
353:    ROTINVA(1,1)=1.0D0; ROTINVA(2,2)=1.0D0; ROTINVA(3,3)=1.0D0331:    ROTINVA(1,1)=1.0D0; ROTINVA(2,2)=1.0D0; ROTINVA(3,3)=1.0D0
354:    RMAT(1:3,1:3)=0.0D0332:    RMAT(1:3,1:3)=0.0D0
355:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0333:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0
356:    CMX=0.0D0; CMY=0.0D0; CMZ=0.0D0334:    CMX=0.0D0; CMY=0.0D0; CMZ=0.0D0
357:    IF (.NOT.MKTRAPT) THEN335:    DO I=1,NATOMS
358:       DO I=1,NATOMS336:       CMX=CMX+DUMMYA(3*(I-1)+1)
359:          CMX=CMX+DUMMYA(3*(I-1)+1)337:       CMY=CMY+DUMMYA(3*(I-1)+2)
360:          CMY=CMY+DUMMYA(3*(I-1)+2)338:       CMZ=CMZ+DUMMYA(3*(I-1)+3)
361:          CMZ=CMZ+DUMMYA(3*(I-1)+3)339:    ENDDO
362:       ENDDO340:    CMX=CMX/NATOMS; CMY=CMY/NATOMS; CMZ=CMZ/NATOMS
363:       CMX=CMX/NATOMS; CMY=CMY/NATOMS; CMZ=CMZ/NATOMS341:    DO I=1,NATOMS
364:       DO I=1,NATOMS342:       DUMMYA(3*(I-1)+1)=DUMMYA(3*(I-1)+1)-CMX
365:          DUMMYA(3*(I-1)+1)=DUMMYA(3*(I-1)+1)-CMX343:       DUMMYA(3*(I-1)+2)=DUMMYA(3*(I-1)+2)-CMY
366:          DUMMYA(3*(I-1)+2)=DUMMYA(3*(I-1)+2)-CMY344:       DUMMYA(3*(I-1)+3)=DUMMYA(3*(I-1)+3)-CMZ
367:          DUMMYA(3*(I-1)+3)=DUMMYA(3*(I-1)+3)-CMZ345:    ENDDO
368:       ENDDO346:    CMX=0.0D0; CMY=0.0D0; CMZ=0.0D0
369:       CMX=0.0D0; CMY=0.0D0; CMZ=0.0D0347:    DO I=1,NATOMS
370:       DO I=1,NATOMS348:       CMX=CMX+DUMMYB(3*(I-1)+1)
371:          CMX=CMX+DUMMYB(3*(I-1)+1)349:       CMY=CMY+DUMMYB(3*(I-1)+2)
372:          CMY=CMY+DUMMYB(3*(I-1)+2)350:       CMZ=CMZ+DUMMYB(3*(I-1)+3)
373:          CMZ=CMZ+DUMMYB(3*(I-1)+3)351:    ENDDO
374:       ENDDO352:    CMX=CMX/NATOMS; CMY=CMY/NATOMS; CMZ=CMZ/NATOMS
375:       CMX=CMX/NATOMS; CMY=CMY/NATOMS; CMZ=CMZ/NATOMS353:    DO I=1,NATOMS
376:       DO I=1,NATOMS354:       DUMMYB(3*(I-1)+1)=DUMMYB(3*(I-1)+1)-CMX
377:          DUMMYB(3*(I-1)+1)=DUMMYB(3*(I-1)+1)-CMX355:       DUMMYB(3*(I-1)+2)=DUMMYB(3*(I-1)+2)-CMY
378:          DUMMYB(3*(I-1)+2)=DUMMYB(3*(I-1)+2)-CMY356:       DUMMYB(3*(I-1)+3)=DUMMYB(3*(I-1)+3)-CMZ
379:          DUMMYB(3*(I-1)+3)=DUMMYB(3*(I-1)+3)-CMZ357:    ENDDO
380:       ENDDO 
381:    ENDIF 
382: 358: 
383:    DISTANCE=0.0D0359:    DISTANCE=0.0D0
384:    DO J1=1,3*NATOMS360:    DO J1=1,3*NATOMS
385:       DISTANCE=DISTANCE+(DUMMYB(J1)-DUMMYA(J1))**2361:       DISTANCE=DISTANCE+(DUMMYB(J1)-DUMMYA(J1))**2
386:    ENDDO362:    ENDDO
387:    DISTANCE=SQRT(DISTANCE)363:    DISTANCE=SQRT(DISTANCE)
388:    IF (.NOT.(PERMOPT.OR.PERMINVOPT)) THEN364:    IF (.NOT.(PERMOPT.OR.PERMINVOPT)) THEN
389:       CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)365:       CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
390:       IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'minpermdist> after initial call to NEWMINDIST distance=',DISTANCE366:       IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'minpermdist> after initial call to NEWMINDIST distance=',DISTANCE
391:    ENDIF367:    ENDIF
633: ENDIF609: ENDIF
634: 610: 
635: 50 DISTANCE=DBEST611: 50 DISTANCE=DBEST
636: !612: !
637: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.613: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.
638: !  Rotate XBEST by ROTINVB to put in best correspondence with COORDSB, undoing the reorientation to DUMMYB from MYORIENT. 614: !  Rotate XBEST by ROTINVB to put in best correspondence with COORDSB, undoing the reorientation to DUMMYB from MYORIENT. 
639: !  We should get the same result for ROTINVB * RMATBEST * (COORDSA-CMA) 615: !  We should get the same result for ROTINVB * RMATBEST * (COORDSA-CMA) 
640: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 616: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 
641: !  (aside from a possible permutation of the atom ordering)617: !  (aside from a possible permutation of the atom ordering)
642: !618: !
643:    IF ((NFREEZE.GT.0).OR.MKTRAPT) THEN619:    IF (NFREEZE.GT.0) THEN
644:       XDUMMY=0.0D0620:       XDUMMY=0.0D0
645:       DO J1=1,NATOMS621:       DO J1=1,NATOMS
646:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &622:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
647:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &623:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
648:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2624:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
649:       ENDDO625:       ENDDO
650:    ELSE IF (STOCKT) THEN626:    ELSE IF (STOCKT) THEN
651:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)627:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)
652:       XDUMMY=0.0D0628:       XDUMMY=0.0D0
653:       DO J1=1,(NATOMS/2)629:       DO J1=1,(NATOMS/2)
676:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &652:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
677:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &653:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
678:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2654:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
679:       ENDDO655:       ENDDO
680:    ENDIF656:    ENDIF
681:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN657:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN
682:       WRITE(MYUNIT,'(2(A,F20.10))') 'minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &658:       WRITE(MYUNIT,'(2(A,F20.10))') 'minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &
683:   &                         ' should be ',SQRT(DISTANCE)659:   &                         ' should be ',SQRT(DISTANCE)
684:    ENDIF660:    ENDIF
685: 661: 
686:    IF ((NFREEZE.GT.0).OR.MKTRAPT) THEN662:    IF (NFREEZE.GT.0) THEN
687:       RMATBEST(1:3,1:3)=0.0D0663:       RMATBEST(1:3,1:3)=0.0D0
688:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0664:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
689:    ELSE665:    ELSE
690:       RMATBEST=MATMUL(ROTINVB,RMATBEST)666:       RMATBEST=MATMUL(ROTINVB,RMATBEST)
691:    ENDIF667:    ENDIF
692:    IF (DEBUG) THEN 668:    IF (DEBUG) THEN 
693:       WRITE(MYUNIT, '(A)') 'RMATBEST:'669:       WRITE(MYUNIT, '(A)') 'RMATBEST:'
694:       WRITE(MYUNIT, '(3F20.10)') RMATBEST(1:3,1:3)670:       WRITE(MYUNIT, '(3F20.10)') RMATBEST(1:3,1:3)
695:    ENDIF671:    ENDIF
696:    COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!672:    COORDSA(1:3*NATOMS)=XBEST(1:3*NATOMS) ! finally, best COORDSA should include permutations for DNEB input!
697:  
698: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!673: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
699: !  DO J1=1,(NATOMS/2)674: !  DO J1=1,(NATOMS/2)
700: !     XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-COORDSA(3*(J1-1)+1))**2+ &675: !     XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-COORDSA(3*(J1-1)+1))**2+ &
701: ! &                 (COORDSB(3*(J1-1)+2)-COORDSA(3*(J1-1)+2))**2+ &676: ! &                 (COORDSB(3*(J1-1)+2)-COORDSA(3*(J1-1)+2))**2+ &
702: ! &                 (COORDSB(3*(J1-1)+3)-COORDSA(3*(J1-1)+3))**2677: ! &                 (COORDSB(3*(J1-1)+3)-COORDSA(3*(J1-1)+3))**2
703: !  ENDDO678: !  ENDDO
704: !  PRINT '(A,F20.10)','XDUMMY=',XDUMMY679: !  PRINT '(A,F20.10)','XDUMMY=',XDUMMY
705: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!680: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
706: 681: 
707: DISTANCE=SQRT(DISTANCE)682: DISTANCE=SQRT(DISTANCE)


r29843/newmindist.f90 2016-03-16 18:33:11.890837997 +0000 r29842/newmindist.f90 2016-03-16 18:33:12.686846182 +0000
 24: !   Kearsley, Acta Cryst. A, 45, 208-210, 1989. 24: !   Kearsley, Acta Cryst. A, 45, 208-210, 1989.
 25: ! 25: !
 26: ! jmc As long as zsym isn't 'W' (in which case mind does something special) mind 26: ! jmc As long as zsym isn't 'W' (in which case mind does something special) mind
 27: ! doesn't care what atomic symbol we give it. 27: ! doesn't care what atomic symbol we give it.
 28: ! 28: !
 29: !     ---------------------------------------------------------------------------------------------- 29: !     ----------------------------------------------------------------------------------------------
 30: ! jdf43>        Modified for general angle-axis 30/01/12 30: ! jdf43>        Modified for general angle-axis 30/01/12
 31: !     ---------------------------------------------------------------------------------------------- 31: !     ----------------------------------------------------------------------------------------------
 32:  32: 
 33: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT) 33: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT)
 34: USE COMMONS,ONLY : MYUNIT, MULLERBROWNT, BOXLX, BOXLY, BOXLZ, STOCKT, CSMT, MKTRAPT 34: USE COMMONS,ONLY : MYUNIT, MULLERBROWNT, BOXLX, BOXLY, BOXLZ, STOCKT, CSMT
 35: IMPLICIT NONE 35: IMPLICIT NONE
 36: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN 36: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN
 37: INTEGER,PARAMETER :: LWORK=12 37: INTEGER,PARAMETER :: LWORK=12
 38: DOUBLE PRECISION RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), & 38: DOUBLE PRECISION RA(3*NATOMS), RB(3*NATOMS), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), &
 39:   &              DIAG(4), TEMPA(LWORK), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, & 39:   &              DIAG(4), TEMPA(LWORK), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, &
 40:   &              NCMXB, NCMYB, NCMZB 40:   &              NCMXB, NCMYB, NCMZB
 41: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:) 41: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
 42: LOGICAL BULKT, TWOD, RIGIDBODY, PRESERVET, DEBUG 42: LOGICAL BULKT, TWOD, RIGIDBODY, PRESERVET, DEBUG
 43: CHARACTER(LEN=5) ZUSE 43: CHARACTER(LEN=5) ZUSE
 44: INTEGER NCIT 44: INTEGER NCIT
 45: DOUBLE PRECISION XSHIFT, YSHIFT, ZSHIFT, XSHIFTNEW, YSHIFTNEW, ZSHIFTNEW 45: DOUBLE PRECISION XSHIFT, YSHIFT, ZSHIFT, XSHIFTNEW, YSHIFTNEW, ZSHIFTNEW
 46: DOUBLE PRECISION MYROTMAT(3,3), OMEGATOT(3,3) 46: DOUBLE PRECISION MYROTMAT(3,3), OMEGATOT(3,3)
 47: COMMON /MINDOM/ MYROTMAT, OMEGATOT 47: COMMON /MINDOM/ MYROTMAT, OMEGATOT
 48:  48: 
 49: IF (RIGIDBODY) THEN 49: IF (RIGIDBODY) THEN
 50:    CALL RBMINDIST(RA,RB,NATOMS,DIST,RMAT,DEBUG) 50:    CALL RBMINDIST(RA,RB,NATOMS,DIST,RMAT,DEBUG)
 51:    RETURN 51:    RETURN
 52: ELSEIF (MKTRAPT) THEN 
 53:    DIST=0.0D0 
 54:    DO J1=1,3*NATOMS 
 55:       DIST=DIST+(RA(J1)-RB(J1))**2 
 56:    ENDDO 
 57:    DIST=SQRT(DIST) 
 58:    RETURN 
 59: !  52: ! 
 60: ! Convert rigid body coordinates to Cartesians for rigid bodies.  53: ! Convert rigid body coordinates to Cartesians for rigid bodies. 
 61: ! 54: !
 62: ELSE IF (ZUSE(1:1).EQ.'W') THEN 55: ELSE IF (ZUSE(1:1).EQ.'W') THEN
 63:    ALLOCATE(XA(3*3*(NATOMS/2)),XB(3*3*(NATOMS/2))) 56:    ALLOCATE(XA(3*3*(NATOMS/2)),XB(3*3*(NATOMS/2)))
 64:    NSIZE=3*(NATOMS/2) 57:    NSIZE=3*(NATOMS/2)
 65:    DO J1=1,NATOMS/2 58:    DO J1=1,NATOMS/2
 66:       CALL CONVERT(RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3), & 59:       CALL CONVERT(RA(3*(J1-1)+1),RA(3*(J1-1)+2),RA(3*(J1-1)+3), &
 67:      &        RA(3*(NATOMS/2+J1-1)+1),RA(3*(NATOMS/2+J1-1)+2),RA(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC) 60:      &        RA(3*(NATOMS/2+J1-1)+1),RA(3*(NATOMS/2+J1-1)+2),RA(3*(NATOMS/2+J1-1)+3),OVEC,H1VEC,H2VEC)
 68:       XA(3*(J1-1)+1+1)=OVEC(1) 61:       XA(3*(J1-1)+1+1)=OVEC(1)


r29843/potential.f90 2016-03-16 18:33:12.086840012 +0000 r29842/potential.f90 2016-03-16 18:33:12.886848239 +0000
264:          CALL QCIPOT2(EREAL,X,GRAD)264:          CALL QCIPOT2(EREAL,X,GRAD)
265:       ELSE265:       ELSE
266:          CALL QCIPOT(EREAL,X,GRAD)266:          CALL QCIPOT(EREAL,X,GRAD)
267:       ENDIF267:       ENDIF
268: 268: 
269:    ELSE IF (TOSI) THEN269:    ELSE IF (TOSI) THEN
270:       CALL RAD(X, GRAD, EREAL, GRADT)270:       CALL RAD(X, GRAD, EREAL, GRADT)
271:       IF (EVAPREJECT) RETURN271:       IF (EVAPREJECT) RETURN
272:       CALL TOSIFUMI(X, GRAD, EREAL, GRADT, SECT)272:       CALL TOSIFUMI(X, GRAD, EREAL, GRADT, SECT)
273: 273: 
 274:    ELSE IF (MKTRAPT) THEN
 275:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)
 276: 
274:    ELSE IF (WELCH) THEN277:    ELSE IF (WELCH) THEN
275:       CALL RAD(X, GRAD, EREAL, GRADT)278:       CALL RAD(X, GRAD, EREAL, GRADT)
276:       CALL WEL(X, GRAD, EREAL, GRADT, SECT)279:       CALL WEL(X, GRAD, EREAL, GRADT, SECT)
277: 280: 
278:    ELSE IF (SCT) THEN281:    ELSE IF (SCT) THEN
279:       CALL RAD(X, GRAD, EREAL, GRADT)282:       CALL RAD(X, GRAD, EREAL, GRADT)
280:       CALL SC(X, GRAD, EREAL, GRADT)283:       CALL SC(X, GRAD, EREAL, GRADT)
281:       IF (CPMD) EREAL=EREAL+1.0D6284:       IF (CPMD) EREAL=EREAL+1.0D6
282: 285: 
283:    ELSE IF (MSCT) THEN286:    ELSE IF (MSCT) THEN
811:          CSMGP=CSMGPSAVE814:          CSMGP=CSMGPSAVE
812:          CSMGPINDEX=CSMGPINDEXSAVE815:          CSMGPINDEX=CSMGPINDEXSAVE
813:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)816:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)
814:          CSMNORM=CSMNORMSAVE817:          CSMNORM=CSMNORMSAVE
815:       END IF818:       END IF
816: 819: 
817:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN820:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN
818: !  EREAL is the distance in this case821: !  EREAL is the distance in this case
819:       CALL MINPERMDIST(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)822:       CALL MINPERMDIST(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)
820: 823: 
821:    ELSE IF (MKTRAPT) THEN 
822:       CALL MKTRAP(NATOMS, X, GRAD, EREAL) 
823:    ELSE IF (BLJCLUSTER) THEN824:    ELSE IF (BLJCLUSTER) THEN
824:       CALL RAD(X, GRAD, EREAL, GRADT)825:       CALL RAD(X, GRAD, EREAL, GRADT)
825:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)826:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)
826: 827: 
827:    ELSE IF (BLJCLUSTER_NOCUT) THEN828:    ELSE IF (BLJCLUSTER_NOCUT) THEN
828: ! ds656> Binary LJ without cutoff829: ! ds656> Binary LJ without cutoff
829:       CALL RAD(X, GRAD, EREAL, GRADT)830:       CALL RAD(X, GRAD, EREAL, GRADT)
830:       CALL BLJ_CLUST(X, GRAD, EREAL, GRADT)831:       CALL BLJ_CLUST(X, GRAD, EREAL, GRADT)
831: 832: 
832:    ELSE IF (MLJT) THEN833:    ELSE IF (MLJT) THEN


r29843/takestep.f 2016-03-16 18:33:12.278841988 +0000 r29842/takestep.f 2016-03-16 18:33:13.078850214 +0000
 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).AND.PERIODIC) RETURN  55:       IF ((PERMOPT.OR.PERMINVOPT.OR.DISTOPT).AND.PERIODIC) RETURN 
 56:       IF ((PERMOPT.OR.PERMINVOPT).AND.MKTRAPT) RETURN  
 57:        56:       
 58:       IF (MODEL1T) THEN 57:       IF (MODEL1T) THEN
 59:          RANDOM=(DPRAND()-0.5D0)*2.0D0 58:          RANDOM=(DPRAND()-0.5D0)*2.0D0
 60:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM 59:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM
 61:          RETURN 60:          RETURN
 62:       ENDIF 61:       ENDIF
 63: ! 62: !
 64: !  Calling CENTRE if NORESET is .TRUE. can lead to problems with COORDSO containing an atom 63: !  Calling CENTRE if NORESET is .TRUE. can lead to problems with COORDSO containing an atom
 65: !  outside the permitted radius. Then it may be impossible to take a step that keeps all the 64: !  outside the permitted radius. Then it may be impossible to take a step that keeps all the
 66: !  atoms inside. 65: !  atoms inside.


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0