hdiff output

r30713/commons.f90 2016-07-07 18:30:08.778560298 +0100 r30712/commons.f90 2016-07-07 18:30:09.654572137 +0100
346:       INTEGER, ALLOCATABLE :: INVATOMLISTS(:,:)346:       INTEGER, ALLOCATABLE :: INVATOMLISTS(:,:)
347:       LOGICAL :: KEEPLISTS347:       LOGICAL :: KEEPLISTS
348:       DOUBLE PRECISION :: NBRCUT1,NBRCUT2348:       DOUBLE PRECISION :: NBRCUT1,NBRCUT2
349:       INTEGER :: NVSITES_MAX, NNLIST_MAX349:       INTEGER :: NVSITES_MAX, NNLIST_MAX
350:       DOUBLE PRECISION, ALLOCATABLE :: VSITES(:,:) 350:       DOUBLE PRECISION, ALLOCATABLE :: VSITES(:,:) 
351: 351: 
352: !ds656> Stress tensor 352: !ds656> Stress tensor 
353:       LOGICAL :: STRESST353:       LOGICAL :: STRESST
354:       INTEGER :: STRESS_MODE354:       INTEGER :: STRESS_MODE
355:       DOUBLE PRECISION, ALLOCATABLE :: STRESS(:,:,:)355:       DOUBLE PRECISION, ALLOCATABLE :: STRESS(:,:,:)
356:  
357: !ds656> A saw-tooth temperature protocol 
358:       LOGICAL :: SAWTOOTH 
359:       INTEGER :: SAWTOOTH_NREJMAX 
360:       DOUBLE PRECISION :: SAWTOOTH_TMAX, SAWTOOTH_TFAC 
361:       356:       
362: !   allocatables357: !   allocatables
363: 358: 
364:       INTEGER,ALLOCATABLE,DIMENSION(:) :: MOVABLEATOMLIST         ! a list containing the movable atom indices359:       INTEGER,ALLOCATABLE,DIMENSION(:) :: MOVABLEATOMLIST         ! a list containing the movable atom indices
365:       LOGICAL,ALLOCATABLE,DIMENSION(:) :: MOVABLEATOMLISTLOGICAL  ! is atom i movable?360:       LOGICAL,ALLOCATABLE,DIMENSION(:) :: MOVABLEATOMLISTLOGICAL  ! is atom i movable?
366:       INTEGER,ALLOCATABLE,DIMENSION(:) :: ATOMSINBLOCK            ! for BLOCKMOVE, to split movableatoms into separate blocks361:       INTEGER,ALLOCATABLE,DIMENSION(:) :: ATOMSINBLOCK            ! for BLOCKMOVE, to split movableatoms into separate blocks
367:       INTEGER,ALLOCATABLE,DIMENSION(:) :: NSPECIES(:), NSPECIES_INI(:)             ! for multicomponent systems362:       INTEGER,ALLOCATABLE,DIMENSION(:) :: NSPECIES(:), NSPECIES_INI(:)             ! for multicomponent systems
368: 363: 
369:       INTEGER, ALLOCATABLE, DIMENSION(:) :: NT, NQ, BLOCK, JUMPINT, JUMPTO, NCORE, NSURFMOVES  ! dimension will be NPAR364:       INTEGER, ALLOCATABLE, DIMENSION(:) :: NT, NQ, BLOCK, JUMPINT, JUMPTO, NCORE, NSURFMOVES  ! dimension will be NPAR
370:       INTEGER, ALLOCATABLE, DIMENSION(:) :: HISTVALS, LHISTVALS  ! dimension declared HBINS in keyword365:       INTEGER, ALLOCATABLE, DIMENSION(:) :: HISTVALS, LHISTVALS  ! dimension declared HBINS in keyword


r30713/keywords.f 2016-07-07 18:30:09.074564298 +0100 r30712/keywords.f 2016-07-07 18:30:09.950576137 +0100
529:       SPANSWAPST = .FALSE.529:       SPANSWAPST = .FALSE.
530:       ENPERMST = .FALSE.530:       ENPERMST = .FALSE.
531:       KEEPLISTS = .FALSE.531:       KEEPLISTS = .FALSE.
532:       NBRCUT1 = 0.0D0532:       NBRCUT1 = 0.0D0
533:       NBRCUT2 = 0.0D0533:       NBRCUT2 = 0.0D0
534: 534: 
535: !     ds656> Stress tensor calculation535: !     ds656> Stress tensor calculation
536:       STRESST=.FALSE.536:       STRESST=.FALSE.
537:       STRESS_MODE=1537:       STRESS_MODE=1
538: 538: 
539: !     ds656> Saw-tooth protocol 
540:       SAWTOOTH=.FALSE. 
541:       SAWTOOTH_NREJMAX=0 
542:       SAWTOOTH_TMAX=0.3D0 
543:       SAWTOOTH_TFAC=1.0D0 
544:  
545: !     ds656> Generalised LJ with Yukawa539: !     ds656> Generalised LJ with Yukawa
546:       GLJY = .FALSE.540:       GLJY = .FALSE.
547:       GLJ_EXP = 6541:       GLJ_EXP = 6
548:       YUK_A = 0.0D0542:       YUK_A = 0.0D0
549:       YUK_XI = 1.0D0543:       YUK_XI = 1.0D0
550:       544:       
551:       CHRMMT=.FALSE.545:       CHRMMT=.FALSE.
552:       CHARMMTYPE=1546:       CHARMMTYPE=1
553:       CHARMMDFTBT=.FALSE.547:       CHARMMDFTBT=.FALSE.
554:       ACESOLV=.FALSE.548:       ACESOLV=.FALSE.
1620:             SPANSWAPST=.TRUE.1614:             SPANSWAPST=.TRUE.
1621:          ENDIF1615:          ENDIF
1622:       ELSE IF (WORD.EQ.'ENPERMS') THEN1616:       ELSE IF (WORD.EQ.'ENPERMS') THEN
1623:          ENPERMST=.TRUE.1617:          ENPERMST=.TRUE.
1624: ! ds656> Generalised LJ with Yukawa1618: ! ds656> Generalised LJ with Yukawa
1625:       ELSE IF (WORD .EQ. 'GLJY') THEN1619:       ELSE IF (WORD .EQ. 'GLJY') THEN
1626:          GLJY = .TRUE.1620:          GLJY = .TRUE.
1627:          CALL READI(GLJ_EXP)1621:          CALL READI(GLJ_EXP)
1628:          CALL READF(YUK_A)1622:          CALL READF(YUK_A)
1629:          CALL READF(YUK_XI)1623:          CALL READF(YUK_XI)
1630: ! ds656> Saw-tooth protocol 
1631:       ELSE IF (WORD .EQ. 'SAWTOOTH') THEN 
1632:          SAWTOOTH=.TRUE. 
1633:          IF(NITEMS.GT.1) CALL READF(SAWTOOTH_TMAX) 
1634:          IF(NITEMS.GT.2) CALL READF(SAWTOOTH_TFAC) 
1635:          IF(NITEMS.GT.3) CALL READI(SAWTOOTH_NREJMAX) 
1636: ! ds656> Stress tensor calculation1624: ! ds656> Stress tensor calculation
1637:       ELSE IF (WORD .EQ. 'STRESS') THEN1625:       ELSE IF (WORD .EQ. 'STRESS') THEN
1638:          STRESST=.TRUE.1626:          STRESST=.TRUE.
1639:          IF(NITEMS.GT.1) CALL READI(STRESS_MODE)1627:          IF(NITEMS.GT.1) CALL READI(STRESS_MODE)
1640: ! BLN 1628: ! BLN 
1641:       ELSE IF (WORD.EQ.'BLN') THEN1629:       ELSE IF (WORD.EQ.'BLN') THEN
1642:          BLNT=.TRUE.1630:          BLNT=.TRUE.
1643:          CALL READF(RK_R)1631:          CALL READF(RK_R)
1644:          CALL READF(RK_THETA)1632:          CALL READF(RK_THETA)
1645:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),1633:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),


r30713/mc.F 2016-07-07 18:30:09.370568298 +0100 r30712/mc.F 2016-07-07 18:30:10.238580030 +0100
 37:       USE AMH_INTERFACES, ONLY:E_WRITE 37:       USE AMH_INTERFACES, ONLY:E_WRITE
 38:       USE ROTAMER 38:       USE ROTAMER
 39:  39: 
 40:       IMPLICIT NONE 40:       IMPLICIT NONE
 41: #ifdef MPI 41: #ifdef MPI
 42:       INCLUDE 'mpif.h' 42:       INCLUDE 'mpif.h'
 43:       INTEGER MPIERR 43:       INTEGER MPIERR
 44:       LOGICAL HITANY 44:       LOGICAL HITANY
 45: #endif 45: #endif
 46:  46: 
 47:       INTEGER J1, NSUCCESS(NPAR), NFAIL(NPAR), NFAILT(NPAR), NSUCCESST(NPAR), J2, NSTEPS, JP, REJSTREAK, 47:       INTEGER J1, NSUCCESS(NPAR), NFAIL(NPAR), NFAILT(NPAR), NSUCCESST(NPAR), J2, NSTEPS, JP, 
 48:      1        UNT, ITERATIONS, NSUPERCOUNT, NQTOT, JACCPREV, NREN, NLAST, NSTEPREN, BRUN,QDONE,JBEST(NPAR),GCJBEST(NPAR), 48:      1        UNT, ITERATIONS, NSUPERCOUNT, NQTOT, JACCPREV, NREN, NLAST, NSTEPREN, BRUN,QDONE,JBEST(NPAR),GCJBEST(NPAR),
 49:      2        NRMS, NDONE, I, RNDSEED, J, NTOT, ITRAJ, NEACCEPT, J3, J4, ISTAT, LOCALCOUNT, QNEWRES, QGCBH, GCNATOMSBEST(NPAR), 49:      2        NRMS, NDONE, I, RNDSEED, J, NTOT, ITRAJ, NEACCEPT, J3, J4, ISTAT, LOCALCOUNT, QNEWRES, QGCBH, GCNATOMSBEST(NPAR),
 50:      &        GCJESTP(NPAR), GCNATOMSBESTP(NPAR), GCJBESTP(NPAR) 50:      &        GCJESTP(NPAR), GCNATOMSBESTP(NPAR), GCJBESTP(NPAR)
 51:       INTEGER :: NSYMCALL=0, NULLMOVES(NPAR), NULLMOVEST(NPAR), NSWAPS=0 51:       INTEGER :: NSYMCALL=0, NULLMOVES(NPAR), NULLMOVEST(NPAR), NSWAPS=0
 52:       DOUBLE PRECISION POTEL, SCALEFAC, RANDOM, DPRAND, DPRAND_UNIVERSAL, SAVECOORDS(3*NATOMSALLOC), TEMPCOORDS(3*NATOMSALLOC), 52:       DOUBLE PRECISION POTEL, SCALEFAC, RANDOM, DPRAND, DPRAND_UNIVERSAL, SAVECOORDS(3*NATOMSALLOC), TEMPCOORDS(3*NATOMSALLOC),
 53:      1                 TIME, SPOTEL(NSUPER), SCOORDS(3*NATOMSALLOC,NSUPER), SCREENC(3*NATOMSALLOC), 53:      1                 TIME, SPOTEL(NSUPER), SCOORDS(3*NATOMSALLOC,NSUPER), SCREENC(3*NATOMSALLOC),
 54:      2                 EPPREV(NPAR), QSTART, QFINISH, RANNJ, RMIN, RMINO, RCOORDS(3*NATOMSALLOC),ELASTSYM(NPAR), 54:      2                 EPPREV(NPAR), QSTART, QFINISH, RANNJ, RMIN, RMINO, RCOORDS(3*NATOMSALLOC),ELASTSYM(NPAR),
 55:      3                 RCOORDSO(3*NATOMSALLOC), RVAT(NATOMSALLOC), RVATO(NATOMSALLOC), EPSSAVE, EBEST(NPAR), GCEBEST(NPAR),  55:      3                 RCOORDSO(3*NATOMSALLOC), RVAT(NATOMSALLOC), RVATO(NATOMSALLOC), EPSSAVE, EBEST(NPAR), GCEBEST(NPAR), 
 56:      4                 GCEBESTP(NPAR), 56:      4                 GCEBESTP(NPAR),
 57:      4                 BESTCOORDS(3*NATOMSALLOC,NPAR), RMSD, VINIT, CTE, TEMPTRAJ(0:NPAR-1), GCBESTCOORDS(3*NATOMSALLOC,NPAR), 57:      4                 BESTCOORDS(3*NATOMSALLOC,NPAR), RMSD, VINIT, CTE, TEMPTRAJ(0:NPAR-1), GCBESTCOORDS(3*NATOMSALLOC,NPAR),
 58:      5                 T, BETA(0:NPAR-1), GRAD(3*NATOMSALLOC), E, W, GCBESTCOORDSP(3*NATOMSALLOC,NPAR), 58:      5                 T, BETA(0:NPAR-1), GRAD(3*NATOMSALLOC), E, W, GCBESTCOORDSP(3*NATOMSALLOC,NPAR),
 59:      &                 DUMMY1, DUMMY2, DUMMY3, INTE, OPOTEL, GCBESTVAT(NATOMSALLOC,NPAR), GCBESTVATP(NATOMSALLOC,NPAR),  59:      &                 DUMMY1, DUMMY2, DUMMY3, INTE, OPOTEL, GCBESTVAT(NATOMSALLOC,NPAR), GCBESTVATP(NATOMSALLOC,NPAR), 
 60:      &                 XMASS, YMASS, ZMASS, CMMAX 60:      &                 XMASS, YMASS, ZMASS, CMMAX
 61:       LOGICAL CHANGEDE, RES1, RES2, COMPON, QUENCHCOMP, MARKOVSTEP 61:       LOGICAL CHANGEDE, RES1, RES2, COMPON, QUENCHCOMP
 62:       LOGICAL CHIRALFAIL,AMIDEFAIL, LOGDUMMY, DISTOK 62:       LOGICAL CHIRALFAIL,AMIDEFAIL, LOGDUMMY, DISTOK
 63:       CHARACTER FNAME*9 63:       CHARACTER FNAME*9
 64:       CHARACTER (LEN= 3)  ISTR 64:       CHARACTER (LEN= 3)  ISTR
 65:       CHARACTER (LEN=20) QUENCHNUM, QUNAME 65:       CHARACTER (LEN=20) QUENCHNUM, QUNAME
 66:       CHARACTER (LEN=20) BESTNAME 66:       CHARACTER (LEN=20) BESTNAME
 67:       LOGICAL  CHANGED 67:       LOGICAL  CHANGED
 68: c  AMH  68: c  AMH 
 69:       INTEGER :: GLY_COUNT,III,I2,I500,SNAPCOUNT 69:       INTEGER :: GLY_COUNT,III,I2,I500,SNAPCOUNT
 70:       DOUBLE PRECISION prcord(NATOMSALLOC,3,3,3) 70:       DOUBLE PRECISION prcord(NATOMSALLOC,3,3,3)
 71:       DOUBLE PRECISION :: mctemp 71:       DOUBLE PRECISION :: mctemp
262: 262: 
263:       IF (NACCEPT.EQ.0) NACCEPT=NSTEPS+1263:       IF (NACCEPT.EQ.0) NACCEPT=NSTEPS+1
264:       NRMS=0264:       NRMS=0
265:       NLAST=0265:       NLAST=0
266:       STAY=.FALSE.266:       STAY=.FALSE.
267:       JACCPREV=0267:       JACCPREV=0
268:       NQTOT=0268:       NQTOT=0
269:       RMINO=1.0D100269:       RMINO=1.0D100
270:       RMIN=1.0D100270:       RMIN=1.0D100
271:       NREN=NRENORM271:       NREN=NRENORM
272:       REJSTREAK=0 
273: #ifdef MPI272: #ifdef MPI
274: #else273: #else
275:       DO JP=1,NPAR 274:       DO JP=1,NPAR 
276: #endif275: #endif
277:          TMOVE(JP)=.TRUE.276:          TMOVE(JP)=.TRUE.
278:          OMOVE(JP)=.TRUE.277:          OMOVE(JP)=.TRUE.
279:          NSUCCESS(JP)=0278:          NSUCCESS(JP)=0
280:          NFAIL(JP)=0279:          NFAIL(JP)=0
281:          NSUCCESST(JP)=0280:          NSUCCESST(JP)=0
282:          NFAILT(JP)=0281:          NFAILT(JP)=0
1454:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)1453:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)
1455:                   IF (DEBUG) THEN1454:                   IF (DEBUG) THEN
1456:                      WRITE(MYUNIT,'(A)') 'Taboo list:'1455:                      WRITE(MYUNIT,'(A)') 'Taboo list:'
1457:                      WRITE(MYUNIT,'(6F20.10)') (ESAVE(J2,1),J2=1,NT(1))1456:                      WRITE(MYUNIT,'(6F20.10)') (ESAVE(J2,1),J2=1,NT(1))
1458:                      WRITE(MYUNIT,73) JP,J1,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)1457:                      WRITE(MYUNIT,73) JP,J1,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
1459: 73                   FORMAT('JP,J1,POTEL,EPREV,NSUC,NFAIL=',I2,I6,2F15.7,2I6,' TABOO')1458: 73                   FORMAT('JP,J1,POTEL,EPREV,NSUC,NFAIL=',I2,I6,2F15.7,2I6,' TABOO')
1460:                   ENDIF1459:                   ENDIF
1461:                   GOTO 231460:                   GOTO 23
1462:                ENDIF1461:                ENDIF
1463: 1462: 
1464: !ds656> For the saw-tooth temperature protocol, step size is 
1465: !       dynamically adjusted to help escape from the current  
1466: !       Markov minimum, but not too far... 
1467:                IF(SAWTOOTH) THEN 
1468:                   ! Energy criterion may not be enough... 
1469:                   IF(ABS(POTEL-EPREV(JP)) < ECONV) THEN 
1470:                      STEP(JP)=1.10*STEP(JP) 
1471:                   ELSE 
1472:                      STEP(JP)=0.99*STEP(JP) 
1473:                   ENDIF 
1474:                ENDIF 
1475: !ds656> print some energies to test BASWAP1463: !ds656> print some energies to test BASWAP
1476:                IF(BASWAPTEST) THEN1464:                IF(BASWAPTEST) THEN
1477:                   !CALL BGUPTA(COORDS,GRAD,EASWAPQ,.FALSE.)1465:                   !CALL BGUPTA(COORDS,GRAD,EASWAPQ,.FALSE.)
1478:                   IF(NSWAPS == 1) THEN1466:                   IF(NSWAPS == 1) THEN
1479:                      WRITE(MYUNIT, '(2(A,1X,E12.6,1X))') 1467:                      WRITE(MYUNIT, '(2(A,1X,E12.6,1X))') 
1480:      1                    "mc> 1 swap: de_A=",DE1,"and de_B=",DE21468:      1                    "mc> 1 swap: de_A=",DE1,"and de_B=",DE2
1481:                   ENDIF1469:                   ENDIF
1482:                   WRITE(MYUNIT, '(2(A,1X,E12.6,1X),A)') 1470:                   WRITE(MYUNIT, '(2(A,1X,E12.6,1X),A)') 
1483:      2                 "mc> BA swaps: dE=", EASWAP-EPREV(JP), "before and dE=",1471:      2                 "mc> BA swaps: dE=", EASWAP-EPREV(JP), "before and dE=",
1484:      3                 POTEL-EPREV(JP), "after quench."1472:      3                 POTEL-EPREV(JP), "after quench."
1699:               CALL NEWRES(J1,JP,JBEST,EBEST,BESTCOORDS,EPPREV,POTEL,ITERATIONS,TIME,RCOORDS,RMIN,RVAT,BRUN,SCREENC,QDONE,1687:               CALL NEWRES(J1,JP,JBEST,EBEST,BESTCOORDS,EPPREV,POTEL,ITERATIONS,TIME,RCOORDS,RMIN,RVAT,BRUN,SCREENC,QDONE,
1700:      &                    JACCPREV,NSUCCESS,NFAIL,NFAILT,NSUCCESST,RES1,RES2)1688:      &                    JACCPREV,NSUCCESS,NFAIL,NFAILT,NSUCCESST,RES1,RES2)
1701:               IF (RES1.OR.RES2) QNEWRES=01689:               IF (RES1.OR.RES2) QNEWRES=0
1702: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DJW1690: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DJW
1703: !             IF (CENT.AND.(.NOT.SEEDT)) CALL CENTRE2(COORDS(1:3*NATOMS,JP))1691: !             IF (CENT.AND.(.NOT.SEEDT)) CALL CENTRE2(COORDS(1:3*NATOMS,JP))
1704: !             COORDSO(1:3*(NATOMS-NSEED),JP)=COORDS(1:3*(NATOMS-NSEED),JP)1692: !             COORDSO(1:3*(NATOMS-NSEED),JP)=COORDS(1:3*(NATOMS-NSEED),JP)
1705: !             WRITE(MYUNIT,'(A,2G20.10)'),'mc> coordso changed: ',COORDSO(1,JP),COORDS(1,JP)     1693: !             WRITE(MYUNIT,'(A,2G20.10)'),'mc> coordso changed: ',COORDSO(1,JP),COORDS(1,JP)     
1706: !             VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)1694: !             VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
1707: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DJW1695: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DJW
1708:             ENDIF1696:             ENDIF
1709:              
1710:             MARKOVSTEP=.FALSE. !<ds656 
1711:              
1712:             IF (STAY) THEN1697:             IF (STAY) THEN
1713:             ELSE IF (EVAP .and. .not.evapreject) THEN1698:             ELSE IF (EVAP .and. .not.evapreject) THEN
1714:                NFAIL(JP)=NFAIL(JP)+11699:                NFAIL(JP)=NFAIL(JP)+1
1715:                CALL MYRESET(JP,NATOMS,NPAR,NSEED)1700:                CALL MYRESET(JP,NATOMS,NPAR,NSEED)
1716:                IF (DEBUG) THEN1701:                IF (DEBUG) THEN
1717:                   WRITE(MYUNIT,33) JP,J1,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)1702:                   WRITE(MYUNIT,33) JP,J1,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
1718: 33                FORMAT('JP,J1,POTEL,EPREV,NSUC,NFAIL=',I2,I6,2F15.7,2I6,' EVAP,REJ')1703: 33                FORMAT('JP,J1,POTEL,EPREV,NSUC,NFAIL=',I2,I6,2F15.7,2I6,' EVAP,REJ')
1719:                ENDIF1704:                ENDIF
1720:             ELSE1705:             ELSE
1721: !     csw34> A series of tests start here to check if a structure should1706: !     csw34> A series of tests start here to check if a structure should
1991:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,1976:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,
1992:      &                                                 ' but Markov energy=',EPREV(JP)1977:      &                                                 ' but Markov energy=',EPREV(JP)
1993:                         IF (.NOT.(CHEMSHIFT2)) STOP1978:                         IF (.NOT.(CHEMSHIFT2)) STOP
1994:                      ENDIF1979:                      ENDIF
1995:                   ENDIF1980:                   ENDIF
1996:                ENDIF1981:                ENDIF
1997:                1982:                
1998:                IF (RATIOT) THEN1983:                IF (RATIOT) THEN
1999:                   CALL NULLMOVE(SCREENC,COORDSO(:,JP),ATEST,NULLMOVES,JP)1984:                   CALL NULLMOVE(SCREENC,COORDSO(:,JP),ATEST,NULLMOVES,JP)
2000:                ENDIF1985:                ENDIF
 1986: 
2001: ! Accept or reject step. If the quench did not converge then allow a1987: ! Accept or reject step. If the quench did not converge then allow a
2002: ! potenial move, but count it as a rejection in terms of NSUCCESS and1988: ! potenial move, but count it as a rejection in terms of NSUCCESS and
2003: ! NFAIL. This way we will accept a lower minimum if found, but the steps won;t become so big.1989: ! NFAIL. This way we will accept a lower minimum if found, but the steps won;t become so big.
2004: ! However, for TIP5P some cold fusion events that had not actually reached the threshold for1990: ! However, for TIP5P some cold fusion events that had not actually reached the threshold for
2005: ! rejection were actually accepted. Must prevent this!1991: ! rejection were actually accepted. Must prevent this!
2006:                IF (ATEST) THEN1992:                IF (ATEST) THEN
2007: #ifdef MPI1993: #ifdef MPI
2008:                   IF (DEBUG) THEN1994:                   IF (DEBUG) THEN
2009:                      WRITE(MYUNIT,334) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)1995:                      WRITE(MYUNIT,334) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
2010: 334                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' ACC')1996: 334                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' ACC')
2027: !                    IF (RESTART) WRITE(MYUNIT,'(A,I6,A)') ' relaxation time set to ',NRELAX,' steps'2013: !                    IF (RESTART) WRITE(MYUNIT,'(A,I6,A)') ' relaxation time set to ',NRELAX,' steps'
2028:                      JACCPREV=J12014:                      JACCPREV=J1
2029:                   ENDIF2015:                   ENDIF
2030:                   IF (QDONE.EQ.1) THEN2016:                   IF (QDONE.EQ.1) THEN
2031:                      NSUCCESS(JP)=NSUCCESS(JP)+12017:                      NSUCCESS(JP)=NSUCCESS(JP)+1
2032:                   ELSE2018:                   ELSE
2033:                      NFAIL(JP)=NFAIL(JP)+12019:                      NFAIL(JP)=NFAIL(JP)+1
2034:                   ENDIF2020:                   ENDIF
2035:                   EPPREV(JP)=EPREV(JP)2021:                   EPPREV(JP)=EPREV(JP)
2036:                   EPREV(JP)=POTEL2022:                   EPREV(JP)=POTEL
2037:                   REJSTREAK=0 
2038:                   MARKOVSTEP=.TRUE. 
2039: ! jwrm2> if we're using percolate with compression, we need to know whether EPREV was calculated with compression on2023: ! jwrm2> if we're using percolate with compression, we need to know whether EPREV was calculated with compression on
2040:                   IF (DEBUG .AND. PERCOLATET) PERCCOMPMARKOV = QUENCHCOMP2024:                   IF (DEBUG .AND. PERCOLATET) PERCCOMPMARKOV = QUENCHCOMP
2041: ! jwrm2> end2025: ! jwrm2> end
2042: ! HBONDMATRIX 2026: ! HBONDMATRIX 
2043: ! csw34> If quench is part of existing group, update the Markov energy 2027: ! csw34> If quench is part of existing group, update the Markov energy 
2044:                   IF ((AMBERT.OR.AMBER12T.OR.(CUDAT.AND.(CUDAPOT.EQ.'A'))).AND.HBONDMATRIX.AND.HBONDMATCH) THEN2028:                   IF ((AMBERT.OR.AMBER12T.OR.(CUDAT.AND.(CUDAPOT.EQ.'A'))).AND.HBONDMATRIX.AND.HBONDMATCH) THEN
2045:                      HBONDMARKOV(HBONDMATCHN)=POTEL  2029:                      HBONDMARKOV(HBONDMATCHN)=POTEL  
2046: ! csw34> Update the group max energy if necessary2030: ! csw34> Update the group max energy if necessary
2047:                      IF (POTEL.GT.HBONDMAXE(HBONDMATCHN)) THEN2031:                      IF (POTEL.GT.HBONDMAXE(HBONDMATCHN)) THEN
2048:                         WRITE(MYUNIT,'(A,I3,A,G20.10)') ' HBONDMATRIX> New highest energy structure for group ',2032:                         WRITE(MYUNIT,'(A,I3,A,G20.10)') ' HBONDMATRIX> New highest energy structure for group ',
2066:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))2050:                            HBONDGROUPNAME='group'//TRIM(ADJUSTL(NHBONDGROUPSCHAR))
2067:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)2051:                            CALL WRITEHBONDMATRIX(HBONDMAT(:,:),HBONDNRES,HBONDGROUPNAME)
2068:                         ENDIF2052:                         ENDIF
2069:                      ENDIF2053:                      ENDIF
2070:                   ENDIF2054:                   ENDIF
2071: ! END OF HBONDMATRIX2055: ! END OF HBONDMATRIX
2072:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)2056:                   COORDSO(1:3*NATOMS,JP)=COORDS(1:3*NATOMS,JP)
2073:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)2057:                   VATO(1:NATOMS,JP)=VAT(1:NATOMS,JP)
2074:                ELSE2058:                ELSE
2075:                   NFAIL(JP)=NFAIL(JP)+12059:                   NFAIL(JP)=NFAIL(JP)+1
2076:                   REJSTREAK=REJSTREAK+1 
2077:                   MARKOVSTEP=.FALSE. 
2078:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)2060:                   CALL MYRESET(JP,NATOMS,NPAR,NSEED)
2079:                   IF (DEBUG) THEN2061:                   IF (DEBUG) THEN
2080:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)2062:                      WRITE(MYUNIT,36) JP,RANDOM,POTEL,EPREV(JP),NSUCCESS(JP),NFAIL(JP)
2081: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')2063: 36                   FORMAT('JP,RAN,POTEL,EPREV,NSUC,NFAIL=',I2,3F15.7,2I6,' REJ')
2082:                   ENDIF2064:                   ENDIF
2083:                ENDIF2065:                ENDIF
2084:             ENDIF2066:             ENDIF
2085: !           WRITE(MYUNIT,'(A,4F20.10)') 'Q4 values ',QS,QF,QSTART,QFINISH2067: !           WRITE(MYUNIT,'(A,4F20.10)') 'Q4 values ',QS,QF,QSTART,QFINISH
2086: !2068: !
2087: !  Dump coordinates and energy if another run is attempting to jump to this one.2069: !  Dump coordinates and energy if another run is attempting to jump to this one.
2088: !2070: !
2089:             IF ((NPAR.GT.1).AND.(.NOT.NEWJUMP)) CALL DUMPJ(JP,JUMPTO,NPAR,COORDS(1:3*NATOMS,1:NPAR),NATOMS,EPREV)2071:             IF ((NPAR.GT.1).AND.(.NOT.NEWJUMP)) CALL DUMPJ(JP,JUMPTO,NPAR,COORDS(1:3*NATOMS,1:NPAR),NATOMS,EPREV)
2090: !2072: !
2091: !  If RESTART then reseed if we haven t accepted a step in twice the relaxation time.2073: !  If RESTART then reseed if we haven t accepted a step in twice the relaxation time.
2092: !2074: !
2093:             IF (RESTART.AND.(J1-JACCPREV.GT.1.1D0*NRELAX)) CALL REST(ITERATIONS,TIME,J1,RCOORDS,RMIN,RVAT,JACCPREV)2075:             IF (RESTART.AND.(J1-JACCPREV.GT.1.1D0*NRELAX)) CALL REST(ITERATIONS,TIME,J1,RCOORDS,RMIN,RVAT,JACCPREV)
2094: !2076: !
2095: !  Check the acceptance ratio.2077: !  Check the acceptance ratio.
2096: ! 2078: ! 
2097:             IF ((MOD(J1,NACCEPT).EQ.0).AND.(NSEED.EQ.0).AND.(.NOT.STAY).AND.(.NOT.SAWTOOTH)) THEN2079:             IF ((MOD(J1,NACCEPT).EQ.0).AND.(NSEED.EQ.0).AND.(.NOT.STAY)) THEN
2098:                IF (RATIOT) THEN2080:                IF (RATIOT) THEN
2099:                   CALL FIXRATIO(NQ(JP),NULLMOVES,NSUCCESS,JP,NULLMOVEST,NSUCCESST)2081:                   CALL FIXRATIO(NQ(JP),NULLMOVES,NSUCCESS,JP,NULLMOVEST,NSUCCESST)
2100:                ELSE2082:                ELSE
2101:                   CALL ACCREJ(NSUCCESS,NFAIL,JP,NSUCCESST,NFAILT)2083:                   CALL ACCREJ(NSUCCESS,NFAIL,JP,NSUCCESST,NFAILT)
2102:                ENDIF2084:                ENDIF
2103:             ENDIF2085:             ENDIF
2104: !2086: 
2105: ! ds656> Saw-tooth temperature protocol2087:             TEMP(JP)=TEMP(JP)*SCALEFAC
2106: ! 
2107:             IF(SAWTOOTH) THEN  
2108:                IF(MARKOVSTEP) THEN 
2109:                   ! If the current trial minimum has been accepted, 
2110:                   ! then multiply the temperature by a factor... 
2111:                   IF(ABS(TEMP(JP)-SAWTOOTH_TMAX) > 1.0D-10) THEN 
2112:                      TEMP(JP)=TEMP(JP)*SAWTOOTH_TFAC 
2113:                      WRITE(MYUNIT,'(A,F13.10,A,F20.10,A,F9.6)')  
2114:      &                    'mc> New T= ',TEMP(JP),' E= ',EPREV(JP), 
2115:      &                    ' STEP= ',STEP(JP) 
2116:                   ENDIF 
2117:                ELSEIF(REJSTREAK>0.AND.MOD(REJSTREAK,SAWTOOTH_NREJMAX).EQ.0) THEN 
2118:                   ! The number of successive rejections has exceeded 
2119:                   ! the allowed maximum => trapped! So we reset the  
2120:                   ! temperature to the maximum value and reset 
2121:                   ! the rejection count. 
2122:                   IF(ABS(TEMP(JP)-SAWTOOTH_TMAX) < 1.0D-10) THEN 
2123:                      SAWTOOTH_TMAX=1.25D0*SAWTOOTH_TMAX 
2124:                   ENDIF 
2125:                   TEMP(JP)=SAWTOOTH_TMAX 
2126:                   WRITE(MYUNIT,'(A,F13.10,A,F20.10,A,F9.6)')  
2127:      &                 'mc> Reset T= ',TEMP(JP),' E= ',EPREV(JP), 
2128:      &                 ' STEP= ',STEP(JP) 
2129:                   REJSTREAK=0 
2130:                ENDIF 
2131:             ELSE 
2132:                ! Do what has always been done... 
2133:                TEMP(JP)=TEMP(JP)*SCALEFAC 
2134:             ENDIF 
2135:                 
2136: #ifdef MPI2088: #ifdef MPI
2137: #else2089: #else
2138:             IF (HIT) GOTO 372090:             IF (HIT) GOTO 37
2139:             IF (EPREV(JP)<EBEST(JP)) EBEST(JP)=EPREV(JP)2091:             IF (EPREV(JP)<EBEST(JP)) EBEST(JP)=EPREV(JP)
2140:             IF ((REPMATCHT).AND.(ABS(MAXVAL(EBEST(:))-MINVAL(EBEST(:)))<ECONV)) GOTO 372092:             IF ((REPMATCHT).AND.(ABS(MAXVAL(EBEST(:))-MINVAL(EBEST(:)))<ECONV)) GOTO 37
2141: #endif2093: #endif
2142:             IF (DUMPINT.GT.0) THEN2094:             IF (DUMPINT.GT.0) THEN
2143:                IF (MOD(J1,DUMPINT).EQ.0) THEN2095:                IF (MOD(J1,DUMPINT).EQ.0) THEN
2144:                   CALL DUMPSTATE(NQ(JP),EBEST,BESTCOORDS,JBEST,JP)2096:                   CALL DUMPSTATE(NQ(JP),EBEST,BESTCOORDS,JBEST,JP)
2145:                ENDIF2097:                ENDIF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0