hdiff output

r29844/bfgsts.f 2016-03-16 18:32:33.886447181 +0000 r29843/bfgsts.f 2016-03-16 18:32:35.430463776 +0000
383: C              PRINT*,'Optimal and actual values of LWORK=',WORK(1),LWORK383: C              PRINT*,'Optimal and actual values of LWORK=',WORK(1),LWORK
384: C              PRINT*,'Optimal and actual values of ILWORK=',IWORK(1),ILWORK384: C              PRINT*,'Optimal and actual values of ILWORK=',IWORK(1),ILWORK
385:                EVALMIN=DIAG(1)385:                EVALMIN=DIAG(1)
386:                SOVER=0.0D0386:                SOVER=0.0D0
387:                DO J1=1,NOPT387:                DO J1=1,NOPT
388:                   VECS(J1)=ZWORK(J1,1)388:                   VECS(J1)=ZWORK(J1,1)
389:                   IF (ITER.GT.1) SOVER=SOVER+VECS(J1)*VECSP(J1)389:                   IF (ITER.GT.1) SOVER=SOVER+VECS(J1)*VECSP(J1)
390:                ENDDO390:                ENDDO
391:                IF (PTEST) THEN391:                IF (PTEST) THEN
392:                   IF (ITER.GT.1) THEN392:                   IF (ITER.GT.1) THEN
393:                      WRITE(*,'(A,G20.10,A,G20.10)') ' Smallest eigenvalue=',EVALMIN,' overlap with previous vector=',SOVER393:                      WRITE(*,'(A,F15.7,A,F15.7)') ' Smallest eigenvalue=',EVALMIN,' overlap with previous vector=',SOVER
394:                   ELSE394:                   ELSE
395:                      WRITE(*,'(A,G20.10,A,G20.10)') ' Smallest eigenvalue=',EVALMIN395:                      WRITE(*,'(A,F15.7,A,F15.7)') ' Smallest eigenvalue=',EVALMIN
396:                   ENDIF396:                   ENDIF
397:                ENDIF397:                ENDIF
398:                DO I=1,NOPT398:                DO I=1,NOPT
399:                   SCRATCH(I) = COORDS(I)399:                   SCRATCH(I) = COORDS(I)
400:                ENDDO400:                ENDDO
401:                DO I=1,HINDEX401:                DO I=1,HINDEX
402:                   XFOB(I)=0.0D0402:                   XFOB(I)=0.0D0
403:                   DO J = 1,NOPT403:                   DO J = 1,NOPT
404:                      XFOB(I)=XFOB(I)+GRAD(J)*ZWORK(J,I)404:                      XFOB(I)=XFOB(I)+GRAD(J)*ZWORK(J,I)
405:                   ENDDO405:                   ENDDO


r29844/minpermdist.f90 2016-03-16 18:32:34.086449240 +0000 r29843/minpermdist.f90 2016-03-16 18:32:35.642465241 +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, unless we 51: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we
 52: !  are doing an ion trap potential. 52: !  are doing an ion trap potential.
 53: ! 53: !
 54: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST) 54: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 55: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, STOCKT, GEOMDIFFTOL, AMBERT, & 55: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, STOCKT, GEOMDIFFTOL, AMBERT, &
 56:   &            NFREEZE, NABT, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, PULLT, EFIELDT, NTSITES, & 56:   &            NFREEZE, NABT, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, PULLT, EFIELDT, NTSITES, &
 57:   &            RIGIDBODY, PERMDIST, OHCELLT, LPERMDIST, EYTRAPT, MKTRAPT, LOCALPERMCUT, LOCALPERMCUT2, & 57:   &            RIGIDBODY, PERMDIST, OHCELLT, LPERMDIST, EYTRAPT, MKTRAPT, LOCALPERMCUT, LOCALPERMCUT2, &
 58:   &            LOCALPERMCUTINC, NOINVERSION, MIEFT, & 58:   &            LOCALPERMCUTINC, NOINVERSION, MIEFT, &
 59:   &            EDIFFTOL, GMAX, CONVR, ATOMMATCHDIST, NRANROT, GTHOMSONT, GTHOMMET, & ! hk286 59:   &            EDIFFTOL, GMAX, CONVR, ATOMMATCHDIST, NRANROT, GTHOMSONT, GTHOMMET, & ! hk286
 60:   &            PHI4MODT, MCPATHT, AMBER12T, VARIABLES, MKTRAPT 60:   &            PHI4MODT, MCPATHT, AMBER12T, VARIABLES
 61: USE COMMONS,ONLY : NOPT 61: USE COMMONS,ONLY : NOPT
 62: USE MODCHARMM,ONLY : CHRMMT 62: USE MODCHARMM,ONLY : CHRMMT
 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH
 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET
 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE
 66: USE GENRIGID 66: USE GENRIGID
 67: USE AMBER12_INTERFACE_MOD 67: USE AMBER12_INTERFACE_MOD
 68: USE CHIRALITY 68: USE CHIRALITY
 69: IMPLICIT NONE 69: IMPLICIT NONE
 70:  70: 
321:       IF (DEBUG) PRINT*, "msb50 minpermdist using intdistance", DISTANCE321:       IF (DEBUG) PRINT*, "msb50 minpermdist using intdistance", DISTANCE
322:     ENDIF322:     ENDIF
323:    DISTANCE=SQRT(DISTANCE)323:    DISTANCE=SQRT(DISTANCE)
324:    RETURN324:    RETURN
325: 325: 
326: ENDIF326: ENDIF
327: !327: !
328: !  Calculate original centres of mass. sn402 - this is actually the centre of geometry?328: !  Calculate original centres of mass. sn402 - this is actually the centre of geometry?
329: !329: !
330: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0330: CMAX=0.0D0; CMAY=0.0D0; CMAZ=0.0D0
331: IF ((NFREEZE.LE.0).AND.(.NOT.MIEFT).AND.(.NOT.MKTRAPT)) THEN331: IF (NFREEZE .LE. 0 .AND..NOT.MIEFT) THEN
332:    IF (STOCKT) THEN 332:    IF (STOCKT) THEN 
333:       NRB=(NATOMS/2)333:       NRB=(NATOMS/2)
334:       DO J1=1,NRB334:       DO J1=1,NRB
335:          CMAX=CMAX+COORDSA(3*(J1-1)+1)335:          CMAX=CMAX+COORDSA(3*(J1-1)+1)
336:          CMAY=CMAY+COORDSA(3*(J1-1)+2)336:          CMAY=CMAY+COORDSA(3*(J1-1)+2)
337:          CMAZ=CMAZ+COORDSA(3*(J1-1)+3)337:          CMAZ=CMAZ+COORDSA(3*(J1-1)+3)
338:       ENDDO338:       ENDDO
339:       CMAX=CMAX/NRB; CMAY=CMAY/NRB; CMAZ=CMAZ/NRB339:       CMAX=CMAX/NRB; CMAY=CMAY/NRB; CMAZ=CMAZ/NRB
340:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0340:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0
341:       DO J1=1,NRB341:       DO J1=1,NRB
368: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)368: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)
369: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)369: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)
370: DBEST=1.0D100370: DBEST=1.0D100
371: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)371: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
372: DBEST=DISTANCE**2372: DBEST=DISTANCE**2
373: 373: 
374: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE374: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE
375: 375: 
376: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then376: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then
377: ! translate the origin of structure A to the CoM of structure B.377: ! translate the origin of structure A to the CoM of structure B.
378: IF ((NFREEZE.GT.0).OR.BULKT.OR.MIEFT.OR.MKTRAPT) THEN378: IF ((NFREEZE.GT.0).OR.BULKT.OR.MIEFT) THEN
379:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)379:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
380: ELSE IF (STOCKT) THEN380: ELSE IF (STOCKT) THEN
381:    DO J1=1,(NATOMS/2)381:    DO J1=1,(NATOMS/2)
382:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX382:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX
383:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY383:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY
384:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ384:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ
385:    ENDDO385:    ENDDO
386: ELSE386: ELSE
387:    DO J1=1,NATOMS387:    DO J1=1,NATOMS
388:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX388:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX
504:              RMATBEST(1:3,1:3)=0.0D0504:              RMATBEST(1:3,1:3)=0.0D0
505:              RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0505:              RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
506:              RETURN506:              RETURN
507:             ENDIF507:             ENDIF
508:          ENDIF508:          ENDIF
509:       ELSE509:       ELSE
510:          CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)510:          CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
511:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULK/NEWMINDIST distance=',DISTANCE511:          IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> after initial call to BULK/NEWMINDIST distance=',DISTANCE
512:          DISTANCE=DISTANCE**2 ! minperdist returns the distance squared for historical reasons512:          DISTANCE=DISTANCE**2 ! minperdist returns the distance squared for historical reasons
513:       ENDIF513:       ENDIF
514:    ELSEIF (MKTRAPT) THEN 
515:       TMAT(1:3,1:3)=0.0D0 
516:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0 
517:       NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1; 
518:       ROTB(1:3,1:3)=0.0D0 
519:       ROTB(1,1)=1.0D0; ROTB(2,2)=1.0D0; ROTB(3,3)=1.0D0 
520:       ROTINVB(1:3,1:3)=0.0D0 
521:       ROTINVB(1,1)=1.0D0; ROTINVB(2,2)=1.0D0; ROTINVB(3,3)=1.0D0 
522:       ROTA(1:3,1:3)=0.0D0 
523:       ROTA(1,1)=1.0D0; ROTA(2,2)=1.0D0; ROTA(3,3)=1.0D0 
524:       ROTINVA(1:3,1:3)=0.0D0 
525:       ROTINVA(1,1)=1.0D0; ROTINVA(2,2)=1.0D0; ROTINVA(3,3)=1.0D0 
526:       RMAT(1:3,1:3)=0.0D0 
527:       RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0 
528:       CMX=0.0D0; CMY=0.0D0; CMZ=0.0D0 
529:       DUMMYA(1:3*NATOMS)=INVERT*COORDSA(1:3*NATOMS) 
530:       DISTANCE=0.0D0 
531:       DO J1=1,3*NATOMS 
532:          DISTANCE=DISTANCE+(DUMMYA(J1)-DUMMYB(J1))**2 
533:       ENDDO 
534:    ELSEIF (STOCKT) THEN514:    ELSEIF (STOCKT) THEN
535:       TMAT(1:3,1:3)=0.0D0515:       TMAT(1:3,1:3)=0.0D0
536:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0516:       TMAT(1,1)=INVERT*1.0D0; TMAT(2,2)=INVERT*1.0D0; TMAT(3,3)=INVERT*1.0D0
537:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMYA,TMAT,0.0D0,0.0D0,0.0D0)517:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMYA,TMAT,0.0D0,0.0D0,0.0D0)
538:       DUMMY(1:3*NATOMS)=DUMMYA(1:3*NATOMS)518:       DUMMY(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
539:       CALL MYORIENT(DUMMYA,DUMMYC,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,NATOMS/2,DEBUG,ROTA,ROTINVA,STOCKT)519:       CALL MYORIENT(DUMMYA,DUMMYC,NORBIT1,NCHOOSE1,NORBIT2,NCHOOSE2,NATOMS/2,DEBUG,ROTA,ROTINVA,STOCKT)
540:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMY,ROTA,0.0D0,0.0D0,0.0D0)520:       CALL NEWROTGEOMSTOCK(NATOMS,DUMMY,ROTA,0.0D0,0.0D0,0.0D0)
541:       DUMMYA(1:3*NATOMS)=DUMMY(1:3*NATOMS)521:       DUMMYA(1:3*NATOMS)=DUMMY(1:3*NATOMS)
542: 522: 
543:       DUMMY(1:3*NATOMS)=DUMMYB(1:3*NATOMS)523:       DUMMY(1:3*NATOMS)=DUMMYB(1:3*NATOMS)
550:    ELSE530:    ELSE
551:       IF ((PULLT.OR.EFIELDT.OR.TWOD).AND.(INVERT.EQ.-1)) THEN ! reflect in xz plane531:       IF ((PULLT.OR.EFIELDT.OR.TWOD).AND.(INVERT.EQ.-1)) THEN ! reflect in xz plane
552:          DO J1=1,NATOMS532:          DO J1=1,NATOMS
553:             DUMMYC(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)533:             DUMMYC(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)
554:             DUMMYC(3*(J1-1)+2)=-DUMMYA(3*(J1-1)+2)534:             DUMMYC(3*(J1-1)+2)=-DUMMYA(3*(J1-1)+2)
555:             DUMMYC(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)535:             DUMMYC(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)
556:          ENDDO536:          ENDDO
557:       ELSE537:       ELSE
558:          DUMMYC(1:3*NATOMS)=INVERT*DUMMYA(1:3*NATOMS)538:          DUMMYC(1:3*NATOMS)=INVERT*DUMMYA(1:3*NATOMS)
559:       ENDIF 539:       ENDIF 
560:       IF ((NRANROT.GT.0).AND.(NROTDONE.LE.NRANROT).AND.(NROTDONE.GT.0).AND.(.NOT.MKTRAPT)) THEN540:       IF ((NRANROT.GT.0).AND.(NROTDONE.LE.NRANROT).AND.(NROTDONE.GT.0)) THEN
561:          IF (DEBUG) PRINT '(A,I6,A,G20.10)',' minpermdist> Trying random starting orientation number ',NROTDONE, &541:          IF (DEBUG) PRINT '(A,I6,A,G20.10)',' minpermdist> Trying random starting orientation number ',NROTDONE, &
562:   &                                         ' minimum distance=',SQRT(DBEST)542:   &                                         ' minimum distance=',SQRT(DBEST)
563:          NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;543:          NORBIT1=1; NORBIT2=1; NORBITB1=1; NORBITB2=1;
564:          ROTB(1:3,1:3)=0.0D0544:          ROTB(1:3,1:3)=0.0D0
565:          ROTB(1,1)=1.0D0; ROTB(2,2)=1.0D0; ROTB(3,3)=1.0D0545:          ROTB(1,1)=1.0D0; ROTB(2,2)=1.0D0; ROTB(3,3)=1.0D0
566:          ROTINVB(1:3,1:3)=0.0D0546:          ROTINVB(1:3,1:3)=0.0D0
567:          ROTINVB(1,1)=1.0D0; ROTINVB(2,2)=1.0D0; ROTINVB(3,3)=1.0D0547:          ROTINVB(1,1)=1.0D0; ROTINVB(2,2)=1.0D0; ROTINVB(3,3)=1.0D0
568:          ROTA(1:3,1:3)=0.0D0548:          ROTA(1:3,1:3)=0.0D0
569:          ROTA(1,1)=1.0D0; ROTA(2,2)=1.0D0; ROTA(3,3)=1.0D0549:          ROTA(1,1)=1.0D0; ROTA(2,2)=1.0D0; ROTA(3,3)=1.0D0
570:          ROTINVA(1:3,1:3)=0.0D0550:          ROTINVA(1:3,1:3)=0.0D0
838: 818: 
839: 50 DISTANCE=DBEST819: 50 DISTANCE=DBEST
840: !820: !
841: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.821: !  XBEST contains the best alignment of A coordinates for the orientation of B coordinates in DUMMYB.
842: !  Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, 822: !  Rotate XBEST by ROTINVBBEST to put in best correspondence with COORDSB, 
843: !  undoing the reorientation to DUMMYB from MYORIENT. 823: !  undoing the reorientation to DUMMYB from MYORIENT. 
844: !  We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA) 824: !  We should get the same result for ROTINVBBEST * RMATBEST * (COORDSA-CMA) 
845: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 825: !  where RMATBEST = +/- RMATCUMUL * ROTA for the best alignment 
846: !  (aside from a possible permutation of the atom ordering)826: !  (aside from a possible permutation of the atom ordering)
847: !827: !
848:    IF ((NFREEZE.GT.0).OR.MIEFT.OR.MKTRAPT) THEN828:    IF (NFREEZE.GT.0 .OR. MIEFT) THEN
849:       XDUMMY=0.0D0829:       XDUMMY=0.0D0
850:       DO J1=1,NATOMS830:       DO J1=1,NATOMS
851:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &831:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
852:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &832:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
853:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2833:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
854:       ENDDO834:       ENDDO
855:    ELSE IF (STOCKT) THEN835:    ELSE IF (STOCKT) THEN
856:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)836:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)
857:       XDUMMY=0.0D0837:       XDUMMY=0.0D0
858:       DO J1=1,(NATOMS/2)838:       DO J1=1,(NATOMS/2)
882:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &862:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
883:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2863:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
884:       ENDDO864:       ENDDO
885:    ENDIF865:    ENDIF
886:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN866:    IF (ABS(SQRT(XDUMMY)-SQRT(DISTANCE)).GT.GEOMDIFFTOL) THEN
887:       PRINT '(2(A,F20.10))',' minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &867:       PRINT '(2(A,F20.10))',' minpermdist> ERROR *** distance between transformed XBEST and COORDSB=',SQRT(XDUMMY), &
888:   &                         ' should be ',SQRT(DISTANCE)868:   &                         ' should be ',SQRT(DISTANCE)
889:       STOP869:       STOP
890:    ENDIF870:    ENDIF
891: 871: 
892:    IF ((NFREEZE.GT.0).OR.MIEFT.OR.MKTRAPT) THEN872:    IF (NFREEZE.GT.0 .OR. MIEFT) THEN
893: !no rotation for NFREEZE .gt. 0873: !no rotation for NFREEZE .gt. 0
894:       RMATBEST(1:3,1:3)=0.0D0874:       RMATBEST(1:3,1:3)=0.0D0
895:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0875:       RMATBEST(1,1)=1.0D0; RMATBEST(2,2)=1.0D0; RMATBEST(3,3)=1.0D0
896:    ELSE876:    ELSE
897:       RMATBEST=MATMUL(ROTINVBBEST,RMATBEST)877:       RMATBEST=MATMUL(ROTINVBBEST,RMATBEST)
898:    ENDIF878:    ENDIF
899: !!!!!!!!!!!!!!!!!!!!!!! DEBUG879: !!!!!!!!!!!!!!!!!!!!!!! DEBUG
900: !880: !
901: ! Test distance for COORDSA with permutation applied in BESTPERM881: ! Test distance for COORDSA with permutation applied in BESTPERM
902: !882: !
909: !  CALL NEWMINDIST(COORDSB,DUMMYA,NATOMS,XDUMMY,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)889: !  CALL NEWMINDIST(COORDSB,DUMMYA,NATOMS,XDUMMY,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
910: !  WRITE(*,'(2(A,G20.10))') ' minpermdist> distance check for permuted COORDSA and original COORDSB=',XDUMMY, &890: !  WRITE(*,'(2(A,G20.10))') ' minpermdist> distance check for permuted COORDSA and original COORDSB=',XDUMMY, &
911: ! &      ' should be ',SQRT(DISTANCE)891: ! &      ' should be ',SQRT(DISTANCE)
912: !!!!!!!!!!!!!!!!!!!!!!! DEBUG892: !!!!!!!!!!!!!!!!!!!!!!! DEBUG
913: 893: 
914: !894: !
915: ! For TRAP potentials we need to preserve the centre of mass because895: ! For TRAP potentials we need to preserve the centre of mass because
916: ! the origin is a fixed point. Try using RMAT best and rotating COORDSA896: ! the origin is a fixed point. Try using RMAT best and rotating COORDSA
917: ! around the origin.897: ! around the origin.
918: !898: !
919:    IF (EYTRAPT) THEN899:    IF (EYTRAPT.OR.MKTRAPT) THEN
920:       XDUMMY=0.0D0900:       XDUMMY=0.0D0
921:       DO J3=1,NATOMS901:       DO J3=1,NATOMS
922:          DUMMYC(3*(J3-1)+1)=COORDSA(3*(BESTPERM(J3)-1)+1)902:          DUMMYC(3*(J3-1)+1)=COORDSA(3*(BESTPERM(J3)-1)+1)
923:          DUMMYC(3*(J3-1)+2)=COORDSA(3*(BESTPERM(J3)-1)+2)903:          DUMMYC(3*(J3-1)+2)=COORDSA(3*(BESTPERM(J3)-1)+2)
924:          DUMMYC(3*(J3-1)+3)=COORDSA(3*(BESTPERM(J3)-1)+3)904:          DUMMYC(3*(J3-1)+3)=COORDSA(3*(BESTPERM(J3)-1)+3)
925:          XDUMMY=XDUMMY+(COORDSB(3*(J3-1)+1)-DUMMYC(3*(J3-1)+1))**2+ &905:          XDUMMY=XDUMMY+(COORDSB(3*(J3-1)+1)-DUMMYC(3*(J3-1)+1))**2+ &
926:   &                    (COORDSB(3*(J3-1)+2)-DUMMYC(3*(J3-1)+2))**2+ &906:   &                    (COORDSB(3*(J3-1)+2)-DUMMYC(3*(J3-1)+2))**2+ &
927:   &                    (COORDSB(3*(J3-1)+3)-DUMMYC(3*(J3-1)+3))**2907:   &                    (COORDSB(3*(J3-1)+3)-DUMMYC(3*(J3-1)+3))**2
928:       ENDDO908:       ENDDO
929: !     PRINT '(A,G20.10)',' minpermdist> with permutations distance is ',SQRT(XDUMMY)909: !     PRINT '(A,G20.10)',' minpermdist> with permutations distance is ',SQRT(XDUMMY)
997: ! IF (DEBUG) PRINT '(A)',' minpermdist> Overall permutation for COORDSA (second argument):'977: ! IF (DEBUG) PRINT '(A)',' minpermdist> Overall permutation for COORDSA (second argument):'
998: ! IF (DEBUG) PRINT '(20I6)',BESTPERM(1:NATOMS)978: ! IF (DEBUG) PRINT '(20I6)',BESTPERM(1:NATOMS)
999: ! PRINT '(I6)',NATOMS979: ! PRINT '(I6)',NATOMS
1000: ! PRINT '(A)','coordsa in minpermdist:'980: ! PRINT '(A)','coordsa in minpermdist:'
1001: ! PRINT '(A,3F20.10)',('LA ',COORDSA(3*(J1-1)+1),COORDSA(3*(J1-1)+2),COORDSA(3*(J1-1)+3),J1=1,NATOMS)981: ! PRINT '(A,3F20.10)',('LA ',COORDSA(3*(J1-1)+1),COORDSA(3*(J1-1)+2),COORDSA(3*(J1-1)+3),J1=1,NATOMS)
1002: ! PRINT '(I6)',NATOMS982: ! PRINT '(I6)',NATOMS
1003: ! PRINT '(A)','coordsb in minpermdist:'983: ! PRINT '(A)','coordsb in minpermdist:'
1004: ! PRINT '(A,3F20.10)',('LB ',COORDSB(3*(J1-1)+1),COORDSB(3*(J1-1)+2),COORDSB(3*(J1-1)+3),J1=1,NATOMS)984: ! PRINT '(A,3F20.10)',('LB ',COORDSB(3*(J1-1)+1),COORDSB(3*(J1-1)+2),COORDSB(3*(J1-1)+3),J1=1,NATOMS)
1005: 985: 
1006: DISTANCE=SQRT(DISTANCE) ! now changed to return distance, not distance^2 22/11/10 DJW986: DISTANCE=SQRT(DISTANCE) ! now changed to return distance, not distance^2 22/11/10 DJW
1007:  
1008: RETURN987: RETURN
1009: END SUBROUTINE MINPERMDIST988: END SUBROUTINE MINPERMDIST
1010: 989: 
1011: SUBROUTINE RANROT(COORDS,ROT,ROTINV,NATOMS)990: SUBROUTINE RANROT(COORDS,ROT,ROTINV,NATOMS)
1012: ! rewritten by js850 to use the new unbiased rotations991: ! rewritten by js850 to use the new unbiased rotations
1013: !992: !
1014: ! subtract the center of mass from coords993: ! subtract the center of mass from coords
1015: ! rotate coords by a uniform random rotation.994: ! rotate coords by a uniform random rotation.
1016: ! return the rotation ROT and the inverse ROTINV995: ! return the rotation ROT and the inverse ROTINV
1017: USE KEY, ONLY : GTHOMSONT, GTHOMMET996: USE KEY, ONLY : GTHOMSONT, GTHOMMET


r29844/MKtrap.f 2016-03-16 18:32:33.678445043 +0000 r29843/MKtrap.f 2016-03-16 18:32:34.930457919 +0000
  1:       SUBROUTINE MKTRAP(N,X,VNEW,ENERGY,SSTEST)  1:       SUBROUTINE MKTRAP(N,X,VNEW,ENERGY,SSTEST)
  2:       USE MODHESS  2:       USE MODHESS
  3:       IMPLICIT NONE  3:       IMPLICIT NONE
  4:       INTEGER N,I,J,J1,J2,J3,J4,J5,J6  4:       INTEGER N,I,J,J1,J2,J3,J4,J5,J6
  5:       DOUBLE PRECISION V0,DELTA,OMEGA,M,B,A,VW,C1,SUM,DIST,  5:       DOUBLE PRECISION V0,DELTA,OMEGA,M,B,A,VW,C1,SUM,DIST,
  6:      & SUM2_11,SUM2_12,SUM2_13,SUM2_22,SUM2_23,SUM2_33,VNEW(3*N),ENERGY,SUM1_1,SUM1_2,SUM1_3  6:      & SUM2_11,SUM2_12,SUM2_13,SUM2_22,SUM2_23,SUM2_33,VNEW(3*N),ENERGY,SUM1_1,SUM1_2,SUM1_3
  7:   7: 
  8:       PARAMETER(OMEGA=0.37D-6)  8:       PARAMETER(OMEGA=0.37D-6)
  9:       PARAMETER(M=1.6D4)  9:       PARAMETER(M=1.6D4)
 10:       PARAMETER(V0=0.33D-6) 10:       PARAMETER(V0=0.33D-6)
 11: !     PARAMETER(DELTA=0.001575D0) ! weak wall 11:       PARAMETER(DELTA=0.001575D0) ! weak wall
 12:       PARAMETER(DELTA=0.004725D0) ! strong wall 12: !     PARAMETER(DELTA=0.004725D0) ! strong wall
 13:       PARAMETER(B=0.0139D-4) 13:       PARAMETER(B=0.0139D-4)
 14:       PARAMETER(A=5.29D-6) ! 5.29 micrometres 14:       PARAMETER(A=5.29D-6) ! 5.29 micrometres
 15:  15: 
 16:       DOUBLE PRECISION X(3*N) 16:       DOUBLE PRECISION X(3*N)
 17:       DOUBLE PRECISION R(N,N),R3(N,N),R5(N,N),V(N) 17:       DOUBLE PRECISION R(N,N),R3(N,N),R5(N,N),V(N)
 18:       DOUBLE PRECISION V1(N,3) 18:       DOUBLE PRECISION V1(N,3)
 19:       DOUBLE PRECISION V2(N,9) 19:       DOUBLE PRECISION V2(N,9)
 20:       LOGICAL SSTEST 20:       LOGICAL SSTEST
 21:  21: 
 22:       VW=V0*DELTA 22:       VW=V0*DELTA
 23:       VNEW(1:3*N)=0.0D0 
 24:       IF (SSTEST) HESS(1:3*N,1:3*N)=0.0D0 23:       IF (SSTEST) HESS(1:3*N,1:3*N)=0.0D0
 25:  24: 
 26: C     N IS THE TOTAL NUMBER OF ATOMS. 25: C     N IS THE TOTAL NUMBER OF ATOMS.
 27: C     V0,VW,OMEGA,M,B ARE PARAMETERS CHARACTERIZING THE POTENTIAL. 26: C     V0,VW,OMEGA,M,B ARE PARAMETERS CHARACTERIZING THE POTENTIAL.
 28: C     A IS THE LENGTH UNIT. 27: C     A IS THE LENGTH UNIT.
 29: C     X(3*N) CONTAIN THE COORDINATES OF THE ATOMS IN THE 3 SPATIAL DIRECTIONS X,Y,Z. 28: C     X(3*N) CONTAIN THE COORDINATES OF THE ATOMS IN THE 3 SPATIAL DIRECTIONS X,Y,Z.
 30: C     R(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)  29: C     R(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS) 
 31: C     R3(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**3  30: C     R3(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**3 
 32: C     R5(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**5  31: C     R5(N,N) IS THE MATRIX WHOSE IJ-TH ELEMENT IS 1/(DISTANCE BETWEEN ITH AND JTH ATOMS)**5 
 33: C     V(I) IS THE POTENTIAL ENERGY OF THE ITH ATOM. 32: C     V(I) IS THE POTENTIAL ENERGY OF THE ITH ATOM.
 37: C     V2(I,1) IS \PARTIAL_XX V(I). 36: C     V2(I,1) IS \PARTIAL_XX V(I).
 38: C     V2(I,2) IS \PARTIAL_XY V(I). 37: C     V2(I,2) IS \PARTIAL_XY V(I).
 39: C     V2(I,3) IS \PARTIAL_XZ V(I). 38: C     V2(I,3) IS \PARTIAL_XZ V(I).
 40: C     V2(I,4) IS \PARTIAL_YX V(I). 39: C     V2(I,4) IS \PARTIAL_YX V(I).
 41: C     V2(I,5) IS \PARTIAL_YY V(I). 40: C     V2(I,5) IS \PARTIAL_YY V(I).
 42: C     V2(I,6) IS \PARTIAL_YZ V(I). 41: C     V2(I,6) IS \PARTIAL_YZ V(I).
 43: C     V2(I,7) IS \PARTIAL_ZX V(I). 42: C     V2(I,7) IS \PARTIAL_ZX V(I).
 44: C     V2(I,8) IS \PARTIAL_ZY V(I). 43: C     V2(I,8) IS \PARTIAL_ZY V(I).
 45: C     V2(I,9) IS \PARTIAL_ZZ V(I). 44: C     V2(I,9) IS \PARTIAL_ZZ V(I).
 46:  45: 
 47:       C1=0.5D0*(OMEGA-M*OMEGA**2-V0) 46:       C1=0.5*(OMEGA-M*OMEGA**2-V0)
 48:  47: 
 49: C     COMPUTING THE DISTANCE MATRIX. 48: C     COMPUTING THE DISTANCE MATRIX.
 50:       CALL RMATMK(N,X,R,R3,R5) 49:       CALL RMATMK(N,X,R,R3,R5)
 51:  50: 
 52:       ENERGY=0.0D0 51:       ENERGY=0.0D0
 53:       DO I=1,N 52:       DO I=1,N
 54:  53: 
 55: c     THE POTENTIAL.  54: c     THE POTENTIAL. 
 56:        SUM=0.0D0 55:        SUM=0.0D0
 57:        DO J=I+1,I+N-1 56:        DO J=I+1,I+N-1
115:         SUM2_22=SUM2_22+((X(3*(I-1)+2)-X(3*(J1-1)+2))**2)*R5(I,J1)114:         SUM2_22=SUM2_22+((X(3*(I-1)+2)-X(3*(J1-1)+2))**2)*R5(I,J1)
116:         SUM2_23=SUM2_23+115:         SUM2_23=SUM2_23+
117:      &  ((X(3*(I-1)+2)-X(3*(J1-1)+2))*(X(3*(I-1)+3)-X(3*(J1-1)+3)))116:      &  ((X(3*(I-1)+2)-X(3*(J1-1)+2))*(X(3*(I-1)+3)-X(3*(J1-1)+3)))
118:      &  *R5(I,J1)117:      &  *R5(I,J1)
119:         SUM2_33=SUM2_33+((X(3*(I-1)+3)-X(3*(J1-1)+3))**2)*R5(I,J1)118:         SUM2_33=SUM2_33+((X(3*(I-1)+3)-X(3*(J1-1)+3))**2)*R5(I,J1)
120:        ENDDO 119:        ENDDO 
121:        V2(I,1)=2*(C1+VW)-B*SUM+3*B*SUM2_11120:        V2(I,1)=2*(C1+VW)-B*SUM+3*B*SUM2_11
122:        V2(I,2)=3*B*SUM2_12121:        V2(I,2)=3*B*SUM2_12
123:        V2(I,3)=3*B*SUM2_13122:        V2(I,3)=3*B*SUM2_13
124:        V2(I,5)=2*(C1-VW)-B*SUM+3*B*SUM2_22123:        V2(I,5)=2*(C1-VW)-B*SUM+3*B*SUM2_22
125:        V2(I,6)=3*B*SUM2_23124:        V2(I,8)=V2(I,6)
126:        V2(I,9)=2*V0-B*SUM+3*B*SUM2_33125:        V2(I,9)=2*V0-B*SUM+3*B*SUM2_33
127:        HESS(3*(I-1)+1,3*(I-1)+1)=V2(I,1)126:        HESS(3*(I-1)+1,3*(I-1)+1)=V2(I,1)
128:        HESS(3*(I-1)+1,3*(I-1)+2)=V2(I,2)127:        HESS(3*(I-1)+1,3*(I-1)+2)=V2(I,2)
129:        HESS(3*(I-1)+1,3*(I-1)+3)=V2(I,3)128:        HESS(3*(I-1)+1,3*(I-1)+3)=V2(I,3)
130:        HESS(3*(I-1)+2,3*(I-1)+1)=V2(I,2)129:        HESS(3*(I-1)+2,3*(I-1)+1)=V2(I,2)
131:        HESS(3*(I-1)+2,3*(I-1)+2)=V2(I,5)130:        HESS(3*(I-1)+2,3*(I-1)+2)=V2(I,5)
132:        HESS(3*(I-1)+2,3*(I-1)+3)=V2(I,6)131:        HESS(3*(I-1)+2,3*(I-1)+3)=V2(I,6)
133:        HESS(3*(I-1)+3,3*(I-1)+1)=V2(I,3)132:        HESS(3*(I-1)+3,3*(I-1)+1)=V2(I,3)
134:        HESS(3*(I-1)+3,3*(I-1)+2)=V2(I,6)133:        HESS(3*(I-1)+3,3*(I-1)+2)=V2(I,6)
135:        HESS(3*(I-1)+3,3*(I-1)+3)=V2(I,9)134:        HESS(3*(I-1)+3,3*(I-1)+3)=V2(I,9)


r29844/newmindist.f90 2016-03-16 18:32:34.310451542 +0000 r29843/newmindist.f90 2016-03-16 18:32:35.846467338 +0000
 48: DOUBLE PRECISION XSHIFT, YSHIFT, ZSHIFT, XSHIFTNEW, YSHIFTNEW, ZSHIFTNEW, BOXLX, BOXLY, BOXLZ 48: DOUBLE PRECISION XSHIFT, YSHIFT, ZSHIFT, XSHIFTNEW, YSHIFTNEW, ZSHIFTNEW, BOXLX, BOXLY, BOXLZ
 49: DOUBLE PRECISION ENERGY, VNEW(3*NATOMS), RMS 49: DOUBLE PRECISION ENERGY, VNEW(3*NATOMS), RMS
 50:  50: 
 51: ! hk286 51: ! hk286
 52: IF (GTHOMSONT) THEN 52: IF (GTHOMSONT) THEN
 53:    CALL GTHOMSONNEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT) 53:    CALL GTHOMSONNEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT)
 54:    RETURN 54:    RETURN
 55: ENDIF 55: ENDIF
 56:  56: 
 57: ! DJW 57: ! DJW
 58: IF (VARIABLES.OR.MKTRAPT) THEN 58: IF (VARIABLES) THEN
 59:    DIST=0.0D0 59:    DIST=0.0D0
 60:    DO J1=1,NOPT 60:    DO J1=1,NOPT
 61:       DIST=DIST+(RA(J1)-RB(J1))**2 61:       DIST=DIST+(RA(J1)-RB(J1))**2
 62:    ENDDO 62:    ENDDO
 63:    DIST=SQRT(DIST) 63:    DIST=SQRT(DIST)
 64:    RETURN 64:    RETURN
 65: ENDIF 65: ENDIF
 66:  66: 
 67: IF (RBAAT) THEN 67: IF (RBAAT) THEN
 68:    CALL RBMINDIST(RA,RB,NATOMS,DIST,RMAT,DEBUG) 68:    CALL RBMINDIST(RA,RB,NATOMS,DIST,RMAT,DEBUG)
387:   &                    RB(3*(NATOMS/2+J1-1)+1),RB(3*(NATOMS/2+J1-1)+2),RB(3*(NATOMS/2+J1-1)+3))387:   &                    RB(3*(NATOMS/2+J1-1)+1),RB(3*(NATOMS/2+J1-1)+2),RB(3*(NATOMS/2+J1-1)+3))
388:       ENDDO388:       ENDDO
389:    ELSEIF (RIGIDBODY) THEN389:    ELSEIF (RIGIDBODY) THEN
390: !390: !
391: !  Needs some thought for the angle/axis rigid body formulation.391: !  Needs some thought for the angle/axis rigid body formulation.
392: !392: !
393:       PRINT '(A)',' newmindist> WARNING *** back transformation not programmed yet for rigid bodies'393:       PRINT '(A)',' newmindist> WARNING *** back transformation not programmed yet for rigid bodies'
394: !394: !
395: !  Preserve centre of mass for EYTRAP, but not orientatation.395: !  Preserve centre of mass for EYTRAP, but not orientatation.
396: !396: !
397: !  ELSEIF (EYTRAPT.OR.MKTRAPT) THEN397:    ELSEIF (EYTRAPT.OR.MKTRAPT) THEN
398:    ELSEIF (EYTRAPT) THEN 
399:       CALL NEWROTGEOM(NSIZE,RB,RMAT,0.0D0,0.0D0,0.0D0)398:       CALL NEWROTGEOM(NSIZE,RB,RMAT,0.0D0,0.0D0,0.0D0)
400:    ELSE399:    ELSE
401: !400: !
402: !  Translate the RB coordinates to the centre of coordinates of RA.401: !  Translate the RB coordinates to the centre of coordinates of RA.
403: !402: !
404:       DO J1=1,NSIZE403:       DO J1=1,NSIZE
405:          RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-CMXB+CMXA+XSHIFT404:          RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-CMXB+CMXA+XSHIFT
406:          RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-CMYB+CMYA+YSHIFT405:          RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-CMYB+CMYA+YSHIFT
407:          RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT406:          RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT
408:       ENDDO407:       ENDDO


r29844/potential.f 2016-03-16 18:32:34.518453682 +0000 r29843/potential.f 2016-03-16 18:32:36.046469395 +0000
1362: !              CALL MKTRAP(NATOMS, COORDS, VNEW, ENERGY,.TRUE.)1362: !              CALL MKTRAP(NATOMS, COORDS, VNEW, ENERGY,.TRUE.)
1363: !              PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS)1363: !              PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS)
1364: !              HESSDUM(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT)1364: !              HESSDUM(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT)
1365: !              DO J1=1,NATOMS1365: !              DO J1=1,NATOMS
1366: !                 COORDS(J1)=COORDS(J1)+DIFF1366: !                 COORDS(J1)=COORDS(J1)+DIFF
1367: !                 CALL MKTRAP(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.)1367: !                 CALL MKTRAP(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.)
1368: !                 COORDS(J1)=COORDS(J1)-2.0D0*DIFF1368: !                 COORDS(J1)=COORDS(J1)-2.0D0*DIFF
1369: !                 CALL MKTRAP(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.)1369: !                 CALL MKTRAP(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.)
1370: !                 COORDS(J1)=COORDS(J1)+DIFF1370: !                 COORDS(J1)=COORDS(J1)+DIFF
1371: !                 IF ((ABS(VNEW(J1)).NE.0.0D0).AND.ABS(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.-1.0D0) THEN1371: !                 IF ((ABS(VNEW(J1)).NE.0.0D0).AND.ABS(100.0D0*(VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.-1.0D0) THEN
1372: !                    WRITE(*,'(I5,2G25.15)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)1372: !                    WRITE(*,'(I5,2F20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
1373: !                 ENDIF1373: !                 ENDIF
1374: !              ENDDO1374: !              ENDDO
1375: !              PRINT*,'analytic and numerical second derivatives:'1375: !              PRINT*,'analytic and numerical second derivatives:'
1376: !              DO J1=1,NATOMS1376: !              DO J1=1,NATOMS
1377: !                 COORDS(J1)=COORDS(J1)+DIFF1377: !                 COORDS(J1)=COORDS(J1)+DIFF
1378: !                 CALL MKTRAP(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.)1378: !                 CALL MKTRAP(NATOMS,COORDS,VPLUS,EPLUS,.FALSE.)
1379: !                 COORDS(J1)=COORDS(J1)-2.0D0*DIFF1379: !                 COORDS(J1)=COORDS(J1)-2.0D0*DIFF
1380: !                 CALL MKTRAP(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.)1380: !                 CALL MKTRAP(NATOMS,COORDS,VMINUS,EMINUS,.FALSE.)
1381: !                 COORDS(J1)=COORDS(J1)+DIFF1381: !                 COORDS(J1)=COORDS(J1)+DIFF
1382: !                 DO J2=1,NATOMS1382: !                 DO J2=1,NATOMS
1383: !                    IF ((ABS(HESS(J1,J2)).GT.1.0D-100)) THEN1383: !                    IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND.
1384: !                       IF ((ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D-3)) THEN1384: !    &                   (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D-3)) THEN
1385: !                          WRITE(*,'(2I5,2G25.15,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'1385: !                    WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'
1386: !                       ELSE 
1387: !                          WRITE(*,'(2I5,2G25.15,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF) 
1388: !                       ENDIF 
1389: !                    ELSE1386: !                    ELSE
1390: !                       WRITE(*,'(2I5,2G25.15,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)1387: !                       WRITE(*,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)
1391: !                    ENDIF1388: !                    ENDIF
1392: !                 ENDDO1389: !                 ENDDO
1393: !              ENDDO1390: !              ENDDO
1394: !          PRINT*,'coords in potential:' 
1395: !          WRITE(*,'(3G25.15)') (COORDS(J1),J1=1,NOPT) 
1396: !              STOP1391: !              STOP
1397: 1392: 
1398:             IF (PTEST) THEN1393:             IF (PTEST) THEN
1399:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY ! ,' hartree'1394:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY ! ,' hartree'
1400:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY ! ,' hartree'1395:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY ! ,' hartree'
1401:             ENDIF1396:             ENDIF
1402:          ELSE IF (ZSYM(NATOMS).EQ.'C1') THEN1397:          ELSE IF (ZSYM(NATOMS).EQ.'C1') THEN
1403:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'1398:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'
1404:             CALL C10(COORDS,NATOMS,VNEW,ENERGY,GTEST,SSTEST)1399:             CALL C10(COORDS,NATOMS,VNEW,ENERGY,GTEST,SSTEST)
1405:             IF (PTEST) THEN1400:             IF (PTEST) THEN
3771:      &         TEMP(5)*SQRT(1.0D0*(NINTS))3766:      &         TEMP(5)*SQRT(1.0D0*(NINTS))
3772:             ELSE IF ( RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN !hk2863767:             ELSE IF ( RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN !hk286
3773:                CALL VSTAT(VNEW(1:DEGFREEDOMS),TEMP,DEGFREEDOMS,DEGFREEDOMS)3768:                CALL VSTAT(VNEW(1:DEGFREEDOMS),TEMP,DEGFREEDOMS,DEGFREEDOMS)
3774:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',3769:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',
3775:      &         TEMP(5)*SQRT(1.0D0*(DEGFREEDOMS))3770:      &         TEMP(5)*SQRT(1.0D0*(DEGFREEDOMS))
3776:                IF (AACONVERGENCET .AND. (TEMP(5) < 5.0D0 * GMAX)) THEN3771:                IF (AACONVERGENCET .AND. (TEMP(5) < 5.0D0 * GMAX)) THEN
3777:                   CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, TEMP(5))3772:                   CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, TEMP(5))
3778:                END IF3773:                END IF
3779:             ELSE3774:             ELSE
3780:                CALL VSTAT(VNEW,TEMP,NOPT,NOPT)3775:                CALL VSTAT(VNEW,TEMP,NOPT,NOPT)
3781:                IF (PTEST) WRITE(*,'(A,43X,G15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',3776:                IF (PTEST) WRITE(*,'(A,43X,F15.10,2X,A,G15.10)') ' potential> RMS force: ',TEMP(5),' |gradient|=',
3782:      &         TEMP(5)*SQRT(1.0D0*(NOPT))3777:      &         TEMP(5)*SQRT(1.0D0*(NOPT))
3783:             ENDIF3778:             ENDIF
3784:             RMS=TEMP(5)3779:             RMS=TEMP(5)
3785:             IF(DEBUG.AND.(RMS.NE.RMS)) THEN3780:             IF(DEBUG.AND.(RMS.NE.RMS)) THEN
3786:                WRITE(*,'(A)' ) 'potential> WARNING - RMS force is NaN - if using AMBER igb=1, can be due to negative Born radii'3781:                WRITE(*,'(A)' ) 'potential> WARNING - RMS force is NaN - if using AMBER igb=1, can be due to negative Born radii'
3787:             END IF3782:             END IF
3788:             ! PRINT '(A,G20.10)',' potential> RMS=',RMS3783:             ! PRINT '(A,G20.10)',' potential> RMS=',RMS
3789:             IF (CPMD.AND.(RMS.EQ.0.0D0)) RMS=1.0D0  !  to prevent convergence when CPMD SCF fails3784:             IF (CPMD.AND.(RMS.EQ.0.0D0)) RMS=1.0D0  !  to prevent convergence when CPMD SCF fails
3790:          ENDIF3785:          ENDIF
3791:          ! 3786:          ! 
3877:             IF (GTEST) VNEW(1:NOPT)=-VNEW(1:NOPT)3872:             IF (GTEST) VNEW(1:NOPT)=-VNEW(1:NOPT)
3878:             IF (SSTEST) HESS(1:NOPT,1:NOPT)=-HESS(1:NOPT,1:NOPT)3873:             IF (SSTEST) HESS(1:NOPT,1:NOPT)=-HESS(1:NOPT,1:NOPT)
3879:          ENDIF3874:          ENDIF
3880: 3875: 
3881:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN3876:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN
3882:         ! sn402: transform back to rigid coords if required3877:         ! sn402: transform back to rigid coords if required
3883:           COORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)3878:           COORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)
3884:           COORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D03879:           COORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D0
3885:       ENDIF3880:       ENDIF
3886: 3881: 
3887:          ! WRITE(*,'(A,G30.20)') 'Energy in POTENTIAL:',ENERGY3882:          ! WRITE(*,'(A,F30.20)') 'Energy in POTENTIAL:',ENERGY
3888:          ! PRINT*,'GTEST,SSTEST=',GTEST,SSTEST3883:          ! PRINT*,'GTEST,SSTEST=',GTEST,SSTEST
3889:          ! WRITE(*,'(A,G20.10)') 'RMS in potential=',RMS3884:          ! WRITE(*,'(A,F20.10)') 'RMS in potential=',RMS
3890:          ! CALL FLUSH(6)3885:          ! CALL FLUSH(6)
3891:          ! PRINT*,'coords in potential:'3886:          ! PRINT*,'coords in potential:'
3892:          ! WRITE(*,'(6G25.15)') (COORDS(J1),J1=1,NOPT)3887:          ! WRITE(*,'(6F20.10)') (COORDS(J1),J1=1,NOPT)
3893:          ! PRINT*,'PARAMS'3888:          ! PRINT*,'PARAMS'
3894:          ! WRITE(*,'(3F20.10)') PARAM1,PARAM2,PARAM33889:          ! WRITE(*,'(3F20.10)') PARAM1,PARAM2,PARAM3
3895:          ! IF (GTEST) PRINT*,'grad:'3890:          ! IF (GTEST) PRINT*,'grad:'
3896:          ! IF (GTEST) WRITE(*,'(6F15.5)') (VNEW(J1),J1=1,NOPT)3891:          ! IF (GTEST) WRITE(*,'(6F15.5)') (VNEW(J1),J1=1,NOPT)
3897:          ! IF (SSTEST) PRINT*,'hess:'3892:          ! IF (SSTEST) PRINT*,'hess:'
3898:          ! IF (SSTEST) WRITE(*,'(6F15.5)') ((HESS(J1,J2),J1=1,NOPT),J2=1,NOPT)3893:          ! IF (SSTEST) WRITE(*,'(6F15.5)') ((HESS(J1,J2),J1=1,NOPT),J2=1,NOPT)
3899: 3894: 
3900:          CALL MYCPU_TIME(TIME,.FALSE.)3895:          CALL MYCPU_TIME(TIME,.FALSE.)
3901:          IF (SSTEST) THEN3896:          IF (SSTEST) THEN
3902:             SCALL=SCALL+13897:             SCALL=SCALL+1


r29844/secdiag.f90 2016-03-16 18:32:34.734455902 +0000 r29843/secdiag.f90 2016-03-16 18:32:36.250471492 +0000
 31: ! 31: !
 32: !           d lambda(x,v) / dv = 2* v * H - 2 * lambda * v 32: !           d lambda(x,v) / dv = 2* v * H - 2 * lambda * v
 33: ! 33: !
 34: !   Again expanding the hessian in powers alows us to compute this without 34: !   Again expanding the hessian in powers alows us to compute this without
 35: !   knowing the actual hessian 35: !   knowing the actual hessian
 36: ! 36: !
 37: ! 37: !
 38:       SUBROUTINE SECDIAG(VEC,COORDS,ENERGY,grad,GL,DIAG,GTEST,XRMS) 38:       SUBROUTINE SECDIAG(VEC,COORDS,ENERGY,grad,GL,DIAG,GTEST,XRMS)
 39:       USE COMMONS 39:       USE COMMONS
 40:       USE KEY 40:       USE KEY
 41:       USE MODCHARMM 41:       USE modcharmm
 42:       USE PORFUNCS 42:       use porfuncs
 43:       IMPLICIT NONE 43:       IMPLICIT NONE
 44:       INTEGER J1 44:       INTEGER J1
 45:       LOGICAL GTEST, FPLUS, FMINUS 45:       LOGICAL GTEST, FPLUS, FMINUS
 46:       DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS) ! the point at which to compute the curvature (this should be intent(IN), but I couldn't easily guarantee that orthogopt didn't modify it) 46:       double precision, intent(INOUT) :: coords(3*natoms) ! the point at which to compute the curvature (this should be intent(IN), but I couldn't easily guarantee that orthogopt didn't modify it)
 47:       DOUBLE PRECISION, INTENT(IN) :: VEC(3*NATOMS) ! the direction along which to compute the curvature (the current best estimate for the lowest eigenvector) 47:       double precision, intent(IN) :: vec(3*natoms) ! the direction along which to compute the curvature (the current best estimate for the lowest eigenvector)
 48:       DOUBLE PRECISION, INTENT(IN) :: ENERGY ! the energy at point coords 48:       double precision, intent(IN) :: energy ! the energy at point coords
 49:       DOUBLE PRECISION, INTENT(IN) :: GRAD(3*NATOMS) ! the gradient at coords (will be nonsense unless nsecdiag > 2) 49:       double precision, intent(IN) :: grad(3*natoms) ! the gradient at coords (will be nonsense unless nsecdiag > 2)
 50:       DOUBLE PRECISION, INTENT(OUT) :: DIAG ! the curvature 50:       double precision, intent(OUT) :: diag ! the curvature
 51:       DOUBLE PRECISION, INTENT(OUT) :: GL(3*NATOMS) ! the gradient of the curvature 51:       double precision, intent(OUT) :: GL(3*natoms) ! the gradient of the curvature
 52:       DOUBLE PRECISION, INTENT(OUT) :: XRMS ! the rms of GL 52:       double precision, intent(out) :: xrms ! the rms of GL
 53:       DOUBLE PRECISION DUMMY3(3*NATOMS), DIFF, DIAG2, DIAG3, & 53:       DOUBLE PRECISION DUMMY3(3*NATOMS), DIFF, DIAG2, DIAG3, &
 54:                        EPLUS,EMINUS,GRAD1(3*NATOMS),LOCALV(3*NATOMS), & 54:                        EPLUS,EMINUS,GRAD1(3*NATOMS),LOCALV(3*NATOMS), &
 55:                        RMS,GRAD2(3*NATOMS), VECL, ZETA, PROJ 55:                        RMS,GRAD2(3*NATOMS), VECL, ZETA, PROJ
 56:       double precision diag4 ! the curvature computed with the first order forward finite differences method 56:       double precision diag4 ! the curvature computed with the first order forward finite differences method
 57:  57: 
 58:       DIAG2 = 0.0D0 58:       DIAG2 = 0.0D0
 59:       DIAG3 = 0.0D0 59:       DIAG3 = 0.0D0
 60: !      EPLUS = 0.0D0 60: !      EPLUS = 0.0D0
 61:       EMINUS = 0.0D0 61:       EMINUS = 0.0D0
 62: !      XRMS = 0.0D0 62: !      XRMS = 0.0D0


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0