hdiff output

r31955/geopt.f 2017-02-22 08:30:08.510494983 +0000 r31954/geopt.f 2017-02-22 08:30:10.670524155 +0000
405:             IF (BULKT.AND.TWOD) NEXMODES=NATOMS+2405:             IF (BULKT.AND.TWOD) NEXMODES=NATOMS+2
406:             IF (PULLT.OR.EFIELDT) NEXMODES=4406:             IF (PULLT.OR.EFIELDT) NEXMODES=4
407: ! hk286407: ! hk286
408:             IF (GTHOMSONT .AND. (GTHOMMET .EQ. 5) ) NEXMODES = 3408:             IF (GTHOMSONT .AND. (GTHOMMET .EQ. 5) ) NEXMODES = 3
409:             IF (GTHOMSONT .AND. (GTHOMMET < 5) ) NEXMODES = 1409:             IF (GTHOMSONT .AND. (GTHOMMET < 5) ) NEXMODES = 1
410:             IF (TWOD) NEXMODES=NEXMODES+NATOMS410:             IF (TWOD) NEXMODES=NEXMODES+NATOMS
411:             IF (FREEZE) THEN411:             IF (FREEZE) THEN
412:                NEXMODES=3*NFREEZE412:                NEXMODES=3*NFREEZE
413:             ENDIF413:             ENDIF
414:             IF (RBAAT) THEN414:             IF (RBAAT) THEN
415:              IF(UNIAXT.AND..NOT.SANDBOXT) THEN415:              IF(UNIAXT) THEN
416:                IF (EFIELDT) THEN416:                IF (EFIELDT) THEN
417:                   NEXMODES = NATOMS/2 + 4417:                   NEXMODES = NATOMS/2 + 4
418:                ELSE418:                ELSE
419:                   NEXMODES = NATOMS/2 + 6419:                   NEXMODES = NATOMS/2 + 6
420:                ENDIF420:                ENDIF
421:              ELSE IF(SANDBOXT) THEN 
422:                NEXMODES=6 
423:                DO J1=1,NATOMS/2 
424:                  IF (UNIAXARRAY(J1)) THEN 
425:                   NEXMODES = NEXMODES + 1 
426:                  END IF 
427:                END DO 
428:              ELSE421:              ELSE
429:                IF (EFIELDT) THEN422:                IF (EFIELDT) THEN
430:                   NEXMODES = 4423:                   NEXMODES = 4
431:                ELSE424:                ELSE
432:                   NEXMODES = 6425:                   NEXMODES = 6
433:                ENDIF426:                ENDIF
434:              END IF427:              END IF
435:                IF (STOCKAAT) NEXMODES = NEXMODES + NATOMS/2428:                IF (STOCKAAT) NEXMODES = NEXMODES + NATOMS/2
436:             ENDIF429:             ENDIF
437:             IF (RINGPOLYMERT) THEN430:             IF (RINGPOLYMERT) THEN
468: !              ALLOCATE(ZWK(NOPT,NEXMODES+1))461: !              ALLOCATE(ZWK(NOPT,NEXMODES+1))
469:                ALLOCATE(ZWK(1,1))462:                ALLOCATE(ZWK(1,1))
470:                ALLOCATE(ZSAVE(NOPT,NOPT))463:                ALLOCATE(ZSAVE(NOPT,NOPT))
471:                ZSAVE(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT)464:                ZSAVE(1:NOPT,1:NOPT)=HESS(1:NOPT,1:NOPT)
472: 465: 
473:                IF (RIGIDINIT) THEN466:                IF (RIGIDINIT) THEN
474:                   RBAANORMALMODET = .TRUE.467:                   RBAANORMALMODET = .TRUE.
475:                   CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)468:                   CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)
476:                   RBAANORMALMODET = .FALSE.469:                   RBAANORMALMODET = .FALSE.
477:                   WRITE(*,*) "geopt> This combination of keywords (LOWESTFRQ) has not been checked - sn402"470:                   WRITE(*,*) "geopt> This combination of keywords (LOWESTFRQ) has not been checked - sn402"
478:                ELSEIF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN  ! hk286471:                ELSEIF (RBAAT.AND..NOT.PYGPERIODICT) THEN  ! hk286
479:                   RBAANORMALMODET = .TRUE.472:                   RBAANORMALMODET = .TRUE.
480:                   CALL NRMLMD (Q, EVALUES, .TRUE.)473:                   CALL NRMLMD (Q, EVALUES, .TRUE.)
481:                   RBAANORMALMODET = .FALSE.474:                   RBAANORMALMODET = .FALSE.
482:                   PRINT *, "geopt> This combination of keywords has not been checked carefully - hk286"475:                   PRINT *, "geopt> This combination of keywords has not been checked carefully - hk286"
483:                ELSE476:                ELSE
484:                   CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NEXMODES+1,ABSTOL,477:                   CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NEXMODES+1,ABSTOL,
485:      &                        NFOUND,EVALUES,478:      &                        NFOUND,EVALUES,
486:      &                        ZWK,NOPT,ISUPPZ,WORK,479:      &                        ZWK,NOPT,ISUPPZ,WORK,
487:      &                        LWORK, IWORK, ILWORK, INFO )480:      &                        LWORK, IWORK, ILWORK, INFO )
488:                ENDIF481:                ENDIF
901: ! Back into the main branch of the subroutine.894: ! Back into the main branch of the subroutine.
902: ! This is where we actually calculate the eigenvectors of the (mass-weighted) Hessian.895: ! This is where we actually calculate the eigenvectors of the (mass-weighted) Hessian.
903:             ELSE896:             ELSE
904:                IF (DUMPV) THEN ! use job type 'V' to get eigenvectors897:                IF (DUMPV) THEN ! use job type 'V' to get eigenvectors
905: ! hk286898: ! hk286
906:                   IF (RIGIDINIT) THEN899:                   IF (RIGIDINIT) THEN
907:                      RBAANORMALMODET = .TRUE.900:                      RBAANORMALMODET = .TRUE.
908:                      ! Eigenvectors will be calculated because DUMPV = .TRUE.901:                      ! Eigenvectors will be calculated because DUMPV = .TRUE.
909:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)902:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)
910:                      RBAANORMALMODET = .FALSE.903:                      RBAANORMALMODET = .FALSE.
911:                   ELSEIF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN904:                   ELSEIF (RBAAT.AND..NOT.PYGPERIODICT) THEN
912:                      RBAANORMALMODET = .TRUE.905:                      RBAANORMALMODET = .TRUE.
913:                      ! Last argument to NRMLMD is EIGENVECTORT. Set this to .TRUE. so eigenvectors are calculated906:                      ! Last argument to NRMLMD is EIGENVECTORT. Set this to .TRUE. so eigenvectors are calculated
914:                      CALL NRMLMD (Q, EVALUES, .TRUE.)907:                      CALL NRMLMD (Q, EVALUES, .TRUE.)
915:                      RBAANORMALMODET = .FALSE.908:                      RBAANORMALMODET = .FALSE.
916:                   ELSE909:                   ELSE
917:                      CALL DSYEV('V','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)910:                      CALL DSYEV('V','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)
918:                   ENDIF911:                   ENDIF
919:                ELSE   ! Then we don't need to calculate the eigenvectors. Eigenvalues only.912:                ELSE   ! Then we don't need to calculate the eigenvectors. Eigenvalues only.
920: ! hk286 - optional number of ENDHESS not yet implemented913: ! hk286 - optional number of ENDHESS not yet implemented
921:                   IF (RIGIDINIT) THEN914:                   IF (RIGIDINIT) THEN
922:                      RBAANORMALMODET = .TRUE.915:                      RBAANORMALMODET = .TRUE.
923:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)916:                      CALL GENRIGID_EIGENVALUES(Q, ATMASS, EVALUES, INFO)
924:                      RBAANORMALMODET = .FALSE.917:                      RBAANORMALMODET = .FALSE.
925:                   ELSE918:                   ELSE
926:                      IF ((NENDHESS.GE.NOPT).OR.(VARIABLES.AND.(NENDHESS.GE.NATOMS))) THEN919:                      IF ((NENDHESS.GE.NOPT).OR.(VARIABLES.AND.(NENDHESS.GE.NATOMS))) THEN
927: C                       CALL DSYEV('N','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)920: C                       CALL DSYEV('N','U',NOPT,HESS,NOPT,EVALUES,TEMPA,9*NATOMS,INFO)
928: C  csw34> changed the call used here to be the same as if NENDHESS < NOPT below921: C  csw34> changed the call used here to be the same as if NENDHESS < NOPT below
929: C         as this leads to a 20% speed increase!922: C         as this leads to a 20% speed increase!
930:                         ABSTOL=DLAMCH('Safe  minimum')923:                         ABSTOL=DLAMCH('Safe  minimum')
931:                         ALLOCATE(ZWK(1,1)) ! not referenced for job type 'N'924:                         ALLOCATE(ZWK(1,1)) ! not referenced for job type 'N'
932:                         IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN ! hk286925:                         IF (RBAAT.AND..NOT.PYGPERIODICT) THEN ! hk286
933:                            RBAANORMALMODET = .TRUE.926:                            RBAANORMALMODET = .TRUE.
934:                            CALL NRMLMD (Q, EVALUES, .TRUE.)927:                            CALL NRMLMD (Q, EVALUES, .TRUE.)
935:                            RBAANORMALMODET = .FALSE.928:                            RBAANORMALMODET = .FALSE.
936:                         ELSE929:                         ELSE
937:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NOPT,ABSTOL,930:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NOPT,ABSTOL,
938:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )931:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )
939:                         ENDIF932:                         ENDIF
940:                         DEALLOCATE(ZWK)933:                         DEALLOCATE(ZWK)
941:                      ELSE934:                      ELSE
942:                         ABSTOL=DLAMCH('Safe  minimum')935:                         ABSTOL=DLAMCH('Safe  minimum')
943:                         ALLOCATE(ZWK(1,1)) ! not referenced for job type 'N'936:                         ALLOCATE(ZWK(1,1)) ! not referenced for job type 'N'
944:                         IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN ! hk286937:                         IF (RBAAT.AND..NOT.PYGPERIODICT) THEN ! hk286
945:                            RBAANORMALMODET = .TRUE.938:                            RBAANORMALMODET = .TRUE.
946:                            CALL NRMLMD (Q, EVALUES, .TRUE.)939:                            CALL NRMLMD (Q, EVALUES, .TRUE.)
947:                            RBAANORMALMODET = .FALSE.940:                            RBAANORMALMODET = .FALSE.
948:                         ELSE941:                         ELSE
949:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NENDHESS,ABSTOL,942:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NENDHESS,ABSTOL,
950:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )943:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )
951:                         ENDIF944:                         ENDIF
952:                         DEALLOCATE(ZWK)945:                         DEALLOCATE(ZWK)
953:                      ENDIF946:                      ENDIF
954:                   ENDIF947:                   ENDIF
1375:                       IF (RINGPOLYMERT) THEN1368:                       IF (RINGPOLYMERT) THEN
1376:                          CALL MASSWTRP(NOPT,RPMASSES,RPDOF)1369:                          CALL MASSWTRP(NOPT,RPMASSES,RPDOF)
1377:                       ELSE1370:                       ELSE
1378:                          ! sn402: we're wasting a bit of time here when RBAAT is TRUE - we call this function1371:                          ! sn402: we're wasting a bit of time here when RBAAT is TRUE - we call this function
1379:                          ! but we then overwrite the mass-weighted hessian in NRMLMD1372:                          ! but we then overwrite the mass-weighted hessian in NRMLMD
1380:                          CALL MASSWT2(NATOMS,ATMASS,Q,VNEW,.TRUE.)1373:                          CALL MASSWT2(NATOMS,ATMASS,Q,VNEW,.TRUE.)
1381:                       ENDIF1374:                       ENDIF
1382:                    ENDIF1375:                    ENDIF
1383: ! Diagonalise1376: ! Diagonalise
1384: ! hk286 - normal modes calculation, for rigid body angle-axis the diagonalization is rather complex1377: ! hk286 - normal modes calculation, for rigid body angle-axis the diagonalization is rather complex
1385:                    IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN1378:                    IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
1386:                       RBAANORMALMODET = .TRUE.1379:                       RBAANORMALMODET = .TRUE.
1387:                       CALL NRMLMD (Q, DIAG, .TRUE.)1380:                       CALL NRMLMD (Q, DIAG, .TRUE.)
1388:                       RBAANORMALMODET = .FALSE.1381:                       RBAANORMALMODET = .FALSE.
1389:                    ELSE1382:                    ELSE
1390:                       CALL DSYEV('V','U',NOPT,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO)1383:                       CALL DSYEV('V','U',NOPT,HESS,SIZE(HESS,1),DIAG,TEMPA,9*NATOMS,INFO)
1391:                    ENDIF1384:                    ENDIF
1392: ! hk2861385: ! hk286
1393: ! If we're freezing atoms, some zero eiganvalue modes can creep in that1386: ! If we're freezing atoms, some zero eiganvalue modes can creep in that
1394: ! are NOT real i.e. one atoms moving and all others stationary. Here, we1387: ! are NOT real i.e. one atoms moving and all others stationary. Here, we
1395: ! check and remove these using ZT1388: ! check and remove these using ZT


r31955/key.f90 2017-02-22 08:30:08.734498011 +0000 r31954/key.f90 2017-02-22 08:30:10.898527234 +0000
 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T
 60:  60: 
 61: ! sy349 > for testing the flatpath after dneb 61: ! sy349 > for testing the flatpath after dneb
 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 63:       LOGICAL FLATPATHT 63:       LOGICAL FLATPATHT
 64:  64: 
 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 68:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 68:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 69:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 69:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
108:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &108:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
109:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &109:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
110:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2110:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2
111: 111: 
112: !     sf344112: !     sf344
113:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &113:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
114:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &114:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
115:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT115:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT
116:  116:  
117:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)117:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)
118:       LOGICAL, ALLOCATABLE :: uniaxarray(:) 
119:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)118:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)
120:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)119:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)
121:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)120:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)
122:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF121:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF
123:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION122:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION
124:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE123:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE
125:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR124:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR
126: !     CHRIS125: !     CHRIS
127: !     MCP126: !     MCP
128:       LOGICAL :: AMHT127:       LOGICAL :: AMHT


r31955/keywords.f 2017-02-22 08:30:08.978501305 +0000 r31954/keywords.f 2017-02-22 08:30:11.130530366 +0000
 60:          USE MULTIPOT 60:          USE MULTIPOT
 61:          ! includes toggle for switching frames - required for RBAA 61:          ! includes toggle for switching frames - required for RBAA
 62:          USE MODHESS, ONLY : RBAANORMALMODET 62:          USE MODHESS, ONLY : RBAANORMALMODET
 63:          ! jdf43> for MMEINITWRAPPER 63:          ! jdf43> for MMEINITWRAPPER
 64:          USE ISO_C_BINDING, ONLY: C_NULL_CHAR 64:          USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 65: ! AMBER12 65: ! AMBER12
 66:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ATOM, AMBER12_RESIDUE, AMBER12_ATOMS, 66:          USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ATOM, AMBER12_RESIDUE, AMBER12_ATOMS,
 67:      &                                    AMBER12_RESIDUES, AMBER12_GET_COORDS 67:      &                                    AMBER12_RESIDUES, AMBER12_GET_COORDS
 68:          USE CHIRALITY 68:          USE CHIRALITY
 69:          USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT 69:          USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT
 70:          use sandbox_module, only : num_atoms 
 71:  70: 
 72:          IMPLICIT NONE 71:          IMPLICIT NONE
 73:  72: 
 74:          DOUBLE PRECISION ::  Q(3*NATOMS) 73:          DOUBLE PRECISION ::  Q(3*NATOMS)
 75:  74: 
 76:          INTEGER NDUM, LUNIT, FUNIT, GETUNIT, MLPDATSTART, MLQDATSTART 75:          INTEGER NDUM, LUNIT, FUNIT, GETUNIT, MLPDATSTART, MLQDATSTART
 77:          INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST, NTYPEA, J1, J2, J3, J, I 76:          INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, IR, LAST, NTYPEA, J1, J2, J3, J, I
 78:          COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR, 77:          COMMON /BUFINF/ ITEM, NITEMS, LOC(132), LINE, SKIPBL, CLEAR, NCR,
 79:      &   NERROR, IR, ECHO, LAST, CAT 78:      &   NERROR, IR, ECHO, LAST, CAT
 80:          ! DOUBLE PRECISION ::AAA,AAB,ABB,PAA,PAB,PBB,QAA,QAB,QBB,ZAA,ZAB 79:          ! DOUBLE PRECISION ::AAA,AAB,ABB,PAA,PAB,PBB,QAA,QAB,QBB,ZAA,ZAB
127:          DOUBLE PRECISION CAPSRHO, CAPSEPS2, CAPSRAD, HEIGHT126:          DOUBLE PRECISION CAPSRHO, CAPSEPS2, CAPSRAD, HEIGHT
128:          COMMON /CAPS/ CAPSRHO, CAPSEPS2, CAPSRAD, HEIGHT127:          COMMON /CAPS/ CAPSRHO, CAPSEPS2, CAPSRAD, HEIGHT
129:          CHARACTER(LEN=20) OSTRING, OTEMP128:          CHARACTER(LEN=20) OSTRING, OTEMP
130:          CHARACTER(LEN=20) :: PINFOSTRING129:          CHARACTER(LEN=20) :: PINFOSTRING
131:          CHARACTER(LEN=5) :: TEMPSTRING130:          CHARACTER(LEN=5) :: TEMPSTRING
132:          CHARACTER(LEN=9) UNSTRING131:          CHARACTER(LEN=9) UNSTRING
133:          CHARACTER(LEN=1) DUMMYCH132:          CHARACTER(LEN=1) DUMMYCH
134:          CHARACTER(LEN=100) TOPFILE,PARFILE133:          CHARACTER(LEN=100) TOPFILE,PARFILE
135:          DOUBLE PRECISION LJREPBB, LJATTBB, LJREPLL, LJATTLL, LJREPNN, LJATTNN,134:          DOUBLE PRECISION LJREPBB, LJATTBB, LJREPLL, LJATTLL, LJREPNN, LJATTNN,
136:      &   HABLN, HBBLN, HCBLN, HDBLN, EABLN, EBBLN, ECBLN, EDBLN, TABLN, TBBLN, TCBLN, TDBLN135:      &   HABLN, HBBLN, HCBLN, HDBLN, EABLN, EBBLN, ECBLN, EDBLN, TABLN, TBBLN, TCBLN, TDBLN
137:          INTEGER IGBNAB, NUNIAX ! sf344136:          INTEGER IGBNAB     ! sf344
138:          ! LOCAL AMH VARIABLES137:          ! LOCAL AMH VARIABLES
139:          INTEGER NRES_AMH, I_RES, GLY_COUNT138:          INTEGER NRES_AMH, I_RES, GLY_COUNT
140:          ! CHARACTER(LEN=5) TARFL139:          ! CHARACTER(LEN=5) TARFL
141:          ! DOUBLE PRECISION X, Y, Z140:          ! DOUBLE PRECISION X, Y, Z
142:          INTEGER :: GROUPCENTRE141:          INTEGER :: GROUPCENTRE
143:          DOUBLE PRECISION :: GROUPRADIUS,DISTGROUPX2,DISTGROUPY2,DISTGROUPZ2,DISTGROUPCENTRE142:          DOUBLE PRECISION :: GROUPRADIUS,DISTGROUPX2,DISTGROUPY2,DISTGROUPZ2,DISTGROUPCENTRE
144:          CHARACTER (LEN=2) :: FREEZEGROUPTYPE143:          CHARACTER (LEN=2) :: FREEZEGROUPTYPE
145:          LOGICAL :: FREEZEGROUPT, TURNOFFCHECKCHIRALITY, MLPDONE, MLPNORM,  MLQDONE, MLQNORM144:          LOGICAL :: FREEZEGROUPT, TURNOFFCHECKCHIRALITY, MLPDONE, MLPNORM,  MLQDONE, MLQNORM
146:          LOGICAL :: RES_IN_LIST145:          LOGICAL :: RES_IN_LIST
147:          DOUBLE PRECISION LPI146:          DOUBLE PRECISION LPI
878:          EYTRAPT=.FALSE.877:          EYTRAPT=.FALSE.
879:          MKTRAPT=.FALSE.878:          MKTRAPT=.FALSE.
880:          BFGSTSTOL=0.0001D0879:          BFGSTSTOL=0.0001D0
881:          EVPC=1.0D0880:          EVPC=1.0D0
882:          OHCELLT=.FALSE.881:          OHCELLT=.FALSE.
883:          CONDATT=.FALSE.882:          CONDATT=.FALSE.
884:          PAIRCOLOURT=.FALSE.883:          PAIRCOLOURT=.FALSE.
885:          NENDDUP=0884:          NENDDUP=0
886:          REVERSEUPHILLT=.FALSE.885:          REVERSEUPHILLT=.FALSE.
887:          ! sf344886:          ! sf344
888:          SANDBOXT=.FALSE. 
889:          EFIELDT=.FALSE.887:          EFIELDT=.FALSE.
890:          MACROIONT=.FALSE.888:          MACROIONT=.FALSE.
891:          NORMALMODET=.FALSE.889:          NORMALMODET=.FALSE.
892:          PYGPERIODICT=.FALSE.890:          PYGPERIODICT=.FALSE.
893:          PYT=.FALSE.891:          PYT=.FALSE.
894:          PYBINARYT=.FALSE.892:          PYBINARYT=.FALSE.
895:          MULTISITEPYT=.FALSE.893:          MULTISITEPYT=.FALSE.
896:          LJGSITET=.FALSE.894:          LJGSITET=.FALSE.
897:          LJSITE=.FALSE.895:          LJSITE=.FALSE.
898:          BLJSITE=.FALSE.896:          BLJSITE=.FALSE.
6095:             IF (NITEMS.GT.4) CALL READF(RPHQMAX)6093:             IF (NITEMS.GT.4) CALL READF(RPHQMAX)
6096:             IF (NITEMS.GT.5) CALL READI(RPHQBINS)6094:             IF (NITEMS.GT.5) CALL READI(RPHQBINS)
6097:             IF (NITEMS.LE.5) THEN6095:             IF (NITEMS.LE.5) THEN
6098:                PRINT '(A)','keywords> ERROR *** insufficient arguments for RPH keyword'6096:                PRINT '(A)','keywords> ERROR *** insufficient arguments for RPH keyword'
6099:                STOP6097:                STOP
6100:             ENDIF6098:             ENDIF
6101: 6099: 
6102:          ELSE IF (WORD .EQ. 'RPLINEAR') THEN6100:          ELSE IF (WORD .EQ. 'RPLINEAR') THEN
6103:             RPCYCLICT=.FALSE.6101:             RPCYCLICT=.FALSE.
6104:             print *, 'use linear polymer'6102:             print *, 'use linear polymer'
6105: 6103: ! 
6106: ! ---------------------------- SANDBOX potential --------------------- 
6107:  
6108:          ELSE IF (WORD .EQ. 'SANDBOX') THEN 
6109:            SANDBOXT = .TRUE. 
6110:            RBAAT = .TRUE. 
6111:            CALL SANDBOX_INPUT(6) 
6112:            IF (NITEMS.GT.1) THEN ! we have uniaxial particles as well in the mixture 
6113:              UNIAXT=.TRUE. 
6114:              CALL READI(nuniax) 
6115:              UNIAXARRAY(:)=.FALSE. 
6116:              DO i=1,nuniax 
6117:                 UNIAXARRAY(i)=.TRUE. 
6118: !                WRITE(*,*) 'uniaxial particle: ', i 
6119:              END DO 
6120:            END IF 
6121:            NTSITES=num_atoms   !num_atoms is the total number of sites 
6122:            ALLOCATE(RBSITE(NTSITES*2/NATOMS,3)) ! will have to change implementation  
6123:                                                 ! to allow for different number of sites in n-ary systems 
6124:  
6125: !            
6126: ! SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS6104: ! SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS
6127: ! 6105: ! 
6128: ! 6106: ! 
6129: ! Save candidate TS`s in SQVV run.6107: ! Save candidate TS`s in SQVV run.
6130: ! 6108: ! 
6131:          ELSE IF (WORD == 'SAVECANDIDATES') THEN6109:          ELSE IF (WORD == 'SAVECANDIDATES') THEN
6132:             SAVECANDIDATES=.TRUE.6110:             SAVECANDIDATES=.TRUE.
6133: ! 6111: ! 
6134: ! SCALE n sets the value of ISTCRT                             - default n=106112: ! SCALE n sets the value of ISTCRT                             - default n=10
6135: ! 6113: ! 


r31955/ncutils.f90 2017-02-22 08:30:08.286491958 +0000 r31954/ncutils.f90 2017-02-22 08:30:10.442521076 +0000
597: ! as atoms in the following subroutine. 597: ! as atoms in the following subroutine. 
598: ! A good solution would probably be to define an unambiguous sense for the598: ! A good solution would probably be to define an unambiguous sense for the
599: ! dipoles so that this problem doesn't arise. It could also cause problems599: ! dipoles so that this problem doesn't arise. It could also cause problems
600: ! in recognising identical minima.600: ! in recognising identical minima.
601: !601: !
602: ! Merges path output files to produce full pathway for the rearrangement;602: ! Merges path output files to produce full pathway for the rearrangement;
603: ! frames in path are reversed as needed;603: ! frames in path are reversed as needed;
604: !604: !
605:      SUBROUTINE MERGEXYZEOFS  605:      SUBROUTINE MERGEXYZEOFS  
606:           USE KEY, ONLY: FILTH,UNRST,FILTHSTR,RIGIDBODY,BULKT,TWOD,STOCKT,STOCKAAT,RBAAT,PERMDIST, MIEFT, &606:           USE KEY, ONLY: FILTH,UNRST,FILTHSTR,RIGIDBODY,BULKT,TWOD,STOCKT,STOCKAAT,RBAAT,PERMDIST, MIEFT, &
607:   &                      AMHT,SEQ,NTSITES,NENDDUP, GTHOMSONT, PAPT, NFREEZE, VARIABLES, NOTRANSROTT, PYGPERIODICT,SANDBOXT ! hk286607:   &                      AMHT,SEQ,NTSITES,NENDDUP, GTHOMSONT, PAPT, NFREEZE, VARIABLES, NOTRANSROTT, PYGPERIODICT ! hk286
608:           USE KEYUTILS        ! frames in bits that are glued together are rotated accordingly;608:           USE KEYUTILS        ! frames in bits that are glued together are rotated accordingly;
609:           USE KEYCONNECT,ONLY : NCNSTEPS609:           USE KEYCONNECT,ONLY : NCNSTEPS
610:           USE COMMONS,ONLY : PARAM1,PARAM2,PARAM3,DEBUG610:           USE COMMONS,ONLY : PARAM1,PARAM2,PARAM3,DEBUG
611:           USE AMHGLOBALS, ONLY : NMRES611:           USE AMHGLOBALS, ONLY : NMRES
612: 612: 
613:           IMPLICIT NONE       ! prerequisites: chain of min/ts constructed; assumes path is dumping plus side of the path613:           IMPLICIT NONE       ! prerequisites: chain of min/ts constructed; assumes path is dumping plus side of the path
614:                               ! first, and there are no blank lines after last frame (!)614:                               ! first, and there are no blank lines after last frame (!)
615:                               ! does somewhat similar with EofS.ts files as well..615:                               ! does somewhat similar with EofS.ts files as well..
616:           DOUBLE PRECISION RMAT(3,3), Q2(4)616:           DOUBLE PRECISION RMAT(3,3), Q2(4)
617:           INTEGER :: I,J,K,EOF,J1,J2,INDEXTS !,FL617:           INTEGER :: I,J,K,EOF,J1,J2,INDEXTS !,FL
1293:              DUMMY=>DUMMY%NEXT1293:              DUMMY=>DUMMY%NEXT
1294:              I=I+11294:              I=I+1
1295:           ENDDO1295:           ENDDO
1296:           DEALLOCATE(LASTFRAME)1296:           DEALLOCATE(LASTFRAME)
1297:           CLOSE(50)1297:           CLOSE(50)
1298:           IF (RBAAT) CLOSE(55)1298:           IF (RBAAT) CLOSE(55)
1299:           CLOSE(41)1299:           CLOSE(41)
1300:      END SUBROUTINE MERGEXYZEOFS1300:      END SUBROUTINE MERGEXYZEOFS
1301: 1301: 
1302:      SUBROUTINE WRITEFRAME(C,S,Q)1302:      SUBROUTINE WRITEFRAME(C,S,Q)
1303:           USE KEY,ONLY : STOCKT, RBAAT, STOCKAAT, NTSITES, PAIRCOLOURT, GTHOMSONT, VARIABLES, PYA1BIN, PYGPERIODICT,SANDBOXT1303:           USE KEY,ONLY : STOCKT, RBAAT, STOCKAAT, NTSITES, PAIRCOLOURT, GTHOMSONT, VARIABLES, PYA1BIN, PYGPERIODICT
1304:           IMPLICIT NONE1304:           IMPLICIT NONE
1305: 1305: 
1306:           CHARACTER(LEN=132),INTENT(IN)         :: C1306:           CHARACTER(LEN=132),INTENT(IN)         :: C
1307:           CHARACTER(LEN=5),POINTER,DIMENSION(:) :: S1307:           CHARACTER(LEN=5),POINTER,DIMENSION(:) :: S
1308:           DOUBLE PRECISION,POINTER,DIMENSION(:) :: Q1308:           DOUBLE PRECISION,POINTER,DIMENSION(:) :: Q
1309:           DOUBLE PRECISION SITES(3*NTSITES), P(3), RM(3,3)1309:           DOUBLE PRECISION SITES(3*NTSITES), P(3), RM(3,3)
1310:           CHARACTER(LEN=2) ZSTRING(NATOMS)1310:           CHARACTER(LEN=2) ZSTRING(NATOMS)
1311: 1311: 
1312:           INTEGER :: J1312:           INTEGER :: J
1313:           DOUBLE PRECISION :: TMPCOORDS(9*NATOMS/2), EulerPhiDeg, EulerPsiDeg, EulerThetaDeg1313:           DOUBLE PRECISION :: TMPCOORDS(9*NATOMS/2), EulerPhiDeg, EulerPsiDeg, EulerThetaDeg
1534: ! format and updated for AMBER, NAB and AMH 30/5/11 DJW.1534: ! format and updated for AMBER, NAB and AMH 30/5/11 DJW.
1535: !1535: !
1536:      SUBROUTINE MAKEPATHINFO1536:      SUBROUTINE MAKEPATHINFO
1537:      USE SYMINF1537:      USE SYMINF
1538:      USE MODHESS1538:      USE MODHESS
1539:      USE MODCHARMM1539:      USE MODCHARMM
1540:      USE PORFUNCS1540:      USE PORFUNCS
1541:      USE MODUNRES1541:      USE MODUNRES
1542:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, RIGIDBODY, NOFRQS, PERMDIST, &1542:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, RIGIDBODY, NOFRQS, PERMDIST, &
1543:   &                 AMHT, SEQ, SDT, NRES_AMH_TEMP, AMBERT, NABT, MACROCYCLET, TTM3T, BOWMANT, &1543:   &                 AMHT, SEQ, SDT, NRES_AMH_TEMP, AMBERT, NABT, MACROCYCLET, TTM3T, BOWMANT, &
1544:   &                 HESSDUMPT,INSTANTONSTARTDUMPT, RBAAT, AMBER12T, VARIABLES, FRQCONV2, PYGPERIODICT, SANDBOXT1544:   &                 HESSDUMPT,INSTANTONSTARTDUMPT, RBAAT, AMBER12T, VARIABLES, FRQCONV2, PYGPERIODICT
1545: 1545: 
1546:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM, PARAM1, PARAM2, PARAM3, DEBUG1546:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM, PARAM1, PARAM2, PARAM3, DEBUG
1547: 1547: 
1548:      USE GENRIGID1548:      USE GENRIGID
1549: 1549: 
1550:      IMPLICIT NONE1550:      IMPLICIT NONE
1551:      DOUBLE PRECISION RMAT(3,3), DIST, DIST21551:      DOUBLE PRECISION RMAT(3,3), DIST, DIST2
1552: 1552: 
1553: !    LOCAL AMH VARIABLES1553: !    LOCAL AMH VARIABLES
1554:      INTEGER :: I_RES, GLY_COUNT1554:      INTEGER :: I_RES, GLY_COUNT
1775:               ENDIF1775:               ENDIF
1776:               WRITE(88,'(I6,1X,A4)') HORDER,FPGRP1776:               WRITE(88,'(I6,1X,A4)') HORDER,FPGRP
1777:               IF (.NOT.NOFRQS) THEN1777:               IF (.NOT.NOFRQS) THEN
1778:                  IF (RIGIDINIT) THEN1778:                  IF (RIGIDINIT) THEN
1779:                     CALL GENRIGID_EIGENVALUES(MI(DUMMY%I)%DATA%X, ATMASS, DIAG, INFO)1779:                     CALL GENRIGID_EIGENVALUES(MI(DUMMY%I)%DATA%X, ATMASS, DIAG, INFO)
1780: 1780: 
1781:                     IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN1781:                     IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
1782:                        CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)1782:                        CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
1783:                     ENDIF1783:                     ENDIF
1784: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame1784: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
1785:                  ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN1785:                  ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
1786:                     RBAANORMALMODET = .TRUE.1786:                     RBAANORMALMODET = .TRUE.
1787:                     CALL POTENTIAL(MI(DUMMY%I)%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)1787:                     CALL POTENTIAL(MI(DUMMY%I)%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
1788:                     CALL NRMLMD (MI(DUMMY%I)%DATA%X, DIAG, .FALSE.)1788:                     CALL NRMLMD (MI(DUMMY%I)%DATA%X, DIAG, .FALSE.)
1789:                     RBAANORMALMODET = .FALSE.1789:                     RBAANORMALMODET = .FALSE.
1790: !                   WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,NOPT)1790: !                   WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,NOPT)
1791:                  ELSE1791:                  ELSE
1792:                     IF (ENDNUMHESS) THEN1792:                     IF (ENDNUMHESS) THEN
1793:                         CALL MAKENUMHESS(MI(DUMMY%I)%DATA%X,NATOMS)1793:                         CALL MAKENUMHESS(MI(DUMMY%I)%DATA%X,NATOMS)
1794:                     ELSE1794:                     ELSE
1795:                        CALL POTENTIAL(MI(DUMMY%I)%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)1795:                        CALL POTENTIAL(MI(DUMMY%I)%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2103: ! hk286 - TS is recorded in rigid body coordinates2103: ! hk286 - TS is recorded in rigid body coordinates
2104: ! sn402 - but for some reason we have ATOMRIGIDCOORDT set to TRUE at this point (which suggests cartesian coordinates)2104: ! sn402 - but for some reason we have ATOMRIGIDCOORDT set to TRUE at this point (which suggests cartesian coordinates)
2105: ! In order for GENRIGID_EIGENVALUES to work, we switch it to FALSE for the duration of this call.2105: ! In order for GENRIGID_EIGENVALUES to work, we switch it to FALSE for the duration of this call.
2106:                     ATOMRIGIDCOORDT = .FALSE.2106:                     ATOMRIGIDCOORDT = .FALSE.
2107:                     CALL GENRIGID_EIGENVALUES(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X, ATMASS, DIAG, INFO)2107:                     CALL GENRIGID_EIGENVALUES(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X, ATMASS, DIAG, INFO)
2108:                     ATOMRIGIDCOORDT = .TRUE.2108:                     ATOMRIGIDCOORDT = .TRUE.
2109:                     IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN2109:                     IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
2110:                         CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)2110:                         CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
2111:                     ENDIF2111:                     ENDIF
2112: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame2112: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
2113:                   ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN2113:                   ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
2114:                     RBAANORMALMODET = .TRUE.2114:                     RBAANORMALMODET = .TRUE.
2115:                     CALL POTENTIAL(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)2115:                     CALL POTENTIAL(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2116:                     CALL NRMLMD (TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X, DIAG, .FALSE.)2116:                     CALL NRMLMD (TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X, DIAG, .FALSE.)
2117:                     RBAANORMALMODET = .FALSE.2117:                     RBAANORMALMODET = .FALSE.
2118: !                   WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,NOPT)2118: !                   WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,NOPT)
2119:                   ELSE2119:                   ELSE
2120:                     IF (ENDNUMHESS) THEN2120:                     IF (ENDNUMHESS) THEN
2121:                         CALL MAKENUMHESS(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X,NATOMS)2121:                         CALL MAKENUMHESS(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X,NATOMS)
2122:                     ELSE2122:                     ELSE
2123:                         CALL POTENTIAL(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)2123:                         CALL POTENTIAL(TS(MI(DUMMY%I)%DATA%CTS(DUMMY%J))%DATA%X,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2245: !2245: !
2246: !  Dump min to path.info2246: !  Dump min to path.info
2247: !2247: !
2248:      SUBROUTINE MAKEMINPATHINFO(MINIMUM)2248:      SUBROUTINE MAKEMINPATHINFO(MINIMUM)
2249:      USE SYMINF 2249:      USE SYMINF 
2250:      USE MODHESS2250:      USE MODHESS
2251:      USE MODCHARMM2251:      USE MODCHARMM
2252:      USE MODUNRES2252:      USE MODUNRES
2253:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, &2253:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, &
2254:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, FRQCONV2, &2254:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, FRQCONV2, &
2255:           PYGPERIODICT, SANDBOXT2255:           PYGPERIODICT
2256:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM2256:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM
2257:      USE PORFUNCS2257:      USE PORFUNCS
2258:      USE GENRIGID2258:      USE GENRIGID
2259:      IMPLICIT NONE2259:      IMPLICIT NONE
2260: 2260: 
2261:      CHARACTER(LEN=20) :: PINFOSTRING2261:      CHARACTER(LEN=20) :: PINFOSTRING
2262:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(3*NATOMS),RMS,DIAG(3*NATOMS),TEMPA(9*NATOMS),DUMQ(3*NATOMS)2262:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(3*NATOMS),RMS,DIAG(3*NATOMS),TEMPA(9*NATOMS),DUMQ(3*NATOMS)
2263:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT, MINIMUM2263:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT, MINIMUM
2264:      LOGICAL :: BTEST,KD,NNZ,NINTB,MINFRQDONE2264:      LOGICAL :: BTEST,KD,NNZ,NINTB,MINFRQDONE
2265:      DOUBLE PRECISION :: QMIN(3*NATOMS),FRQSMIN(3*NATOMS),EMIN,INERTIA(3,3)2265:      DOUBLE PRECISION :: QMIN(3*NATOMS),FRQSMIN(3*NATOMS),EMIN,INERTIA(3,3)
2501:            CALL SYMMETRY(HORDER,.FALSE.,QMIN,INERTIA)2501:            CALL SYMMETRY(HORDER,.FALSE.,QMIN,INERTIA)
2502:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP2502:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP
2503:            IF (.NOT.NOFRQS) THEN2503:            IF (.NOT.NOFRQS) THEN
2504:               IF (RIGIDINIT) THEN2504:               IF (RIGIDINIT) THEN
2505:                  CALL GENRIGID_EIGENVALUES(QMIN, ATMASS, DIAG, INFO)2505:                  CALL GENRIGID_EIGENVALUES(QMIN, ATMASS, DIAG, INFO)
2506: 2506: 
2507:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN2507:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
2508:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)2508:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
2509:                  ENDIF2509:                  ENDIF
2510: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame2510: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
2511:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN2511:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
2512:                  RBAANORMALMODET = .TRUE.2512:                  RBAANORMALMODET = .TRUE.
2513:                  CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)2513:                  CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2514:                  CALL NRMLMD (QMIN, DIAG, .FALSE.)2514:                  CALL NRMLMD (QMIN, DIAG, .FALSE.)
2515:                  RBAANORMALMODET = .FALSE.2515:                  RBAANORMALMODET = .FALSE.
2516:               ELSE2516:               ELSE
2517:                  IF (ENDNUMHESS) THEN2517:                  IF (ENDNUMHESS) THEN
2518:                     CALL MAKENUMHESS(QMIN,NATOMS)2518:                     CALL MAKENUMHESS(QMIN,NATOMS)
2519:                  ELSE2519:                  ELSE
2520:                     CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)2520:                     CALL POTENTIAL(QMIN,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2521:                  ENDIF2521:                  ENDIF
2611: !  Dump TS to path.info 2611: !  Dump TS to path.info 
2612: !2612: !
2613:      SUBROUTINE MAKETSPATHINFO(TSN)2613:      SUBROUTINE MAKETSPATHINFO(TSN)
2614: 2614: 
2615:      USE SYMINF 2615:      USE SYMINF 
2616:      USE MODHESS2616:      USE MODHESS
2617:      USE MODCHARMM2617:      USE MODCHARMM
2618:      USE MODUNRES2618:      USE MODUNRES
2619:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, &2619:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, &
2620:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, FRQCONV2, &2620:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, FRQCONV2, &
2621:           PYGPERIODICT, SANDBOXT2621:           PYGPERIODICT
2622:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM2622:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM
2623:      USE PORFUNCS2623:      USE PORFUNCS
2624:      USE GENRIGID2624:      USE GENRIGID
2625:      IMPLICIT NONE2625:      IMPLICIT NONE
2626: 2626: 
2627:      CHARACTER(LEN=20) :: PINFOSTRING2627:      CHARACTER(LEN=20) :: PINFOSTRING
2628:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(3*NATOMS),RMS,DIAG(3*NATOMS),TEMPA(9*NATOMS),DUMQ(3*NATOMS)2628:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(3*NATOMS),RMS,DIAG(3*NATOMS),TEMPA(9*NATOMS),DUMQ(3*NATOMS)
2629:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT,TSN2629:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT,TSN
2630:      LOGICAL :: BTEST,KD,NNZ,NINTB,TSFRQDONE2630:      LOGICAL :: BTEST,KD,NNZ,NINTB,TSFRQDONE
2631:      DOUBLE PRECISION :: QTS(3*NATOMS),FRQSTS(3*NATOMS),ETS,INERTIA(3,3)2631:      DOUBLE PRECISION :: QTS(3*NATOMS),FRQSTS(3*NATOMS),ETS,INERTIA(3,3)
2880:            IF (.NOT.NOFRQS) THEN2880:            IF (.NOT.NOFRQS) THEN
2881:               IF (RIGIDINIT) THEN2881:               IF (RIGIDINIT) THEN
2882: ! hk286 - TS is recorded in rigid body coordinates2882: ! hk286 - TS is recorded in rigid body coordinates
2883:                  ATOMRIGIDCOORDT = .FALSE.2883:                  ATOMRIGIDCOORDT = .FALSE.
2884:                  CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO)2884:                  CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO)
2885:                  ATOMRIGIDCOORDT = .TRUE.2885:                  ATOMRIGIDCOORDT = .TRUE.
2886:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN2886:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
2887:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)2887:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
2888:                  ENDIF2888:                  ENDIF
2889: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame2889: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
2890:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN2890:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
2891:                  RBAANORMALMODET = .TRUE.2891:                  RBAANORMALMODET = .TRUE.
2892:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)2892:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2893:                  CALL NRMLMD (QTS, DIAG, .FALSE.)2893:                  CALL NRMLMD (QTS, DIAG, .FALSE.)
2894:                  RBAANORMALMODET = .FALSE.2894:                  RBAANORMALMODET = .FALSE.
2895:               ELSE2895:               ELSE
2896:                  IF (ENDNUMHESS) THEN2896:                  IF (ENDNUMHESS) THEN
2897:                     CALL MAKENUMHESS(QTS,NATOMS)2897:                     CALL MAKENUMHESS(QTS,NATOMS)
2898:                  ELSE2898:                  ELSE
2899:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)2899:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
2900:                  ENDIF2900:                  ENDIF
2998: !2998: !
2999: !  Dump the latest min-sad-min triple to path.info in the usual format2999: !  Dump the latest min-sad-min triple to path.info in the usual format
3000: !  3000: !  
3001:      SUBROUTINE MAKEALLPATHINFO(QTS,QPLUS,QMINUS,ETS,EPLUS,EMINUS,FRQSTS,FRQSPLUS,FRQSMINUS)3001:      SUBROUTINE MAKEALLPATHINFO(QTS,QPLUS,QMINUS,ETS,EPLUS,EMINUS,FRQSTS,FRQSPLUS,FRQSMINUS)
3002:      USE SYMINF 3002:      USE SYMINF 
3003:      USE MODHESS3003:      USE MODHESS
3004:      USE MODCHARMM3004:      USE MODCHARMM
3005:      USE MODUNRES3005:      USE MODUNRES
3006:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, &3006:      USE KEY, ONLY: FILTH, FILTHSTR, UNRST, TWOD, BULKT, MACHINE, NOFRQS, AMBERT, NABT, AMHT, SEQ, TARFL, NRES_AMH_TEMP, SDT, &
3007:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, &3007:           AMBER12T, RBAAT, MACROCYCLET, GTHOMSONT, TTM3T, BOWMANT, HESSDUMPT, INSTANTONSTARTDUMPT,FREEZE,NONFREEZE, &
3008:   &       VARIABLES, FRQCONV2, PYGPERIODICT, SANDBOXT3008:   &       VARIABLES, FRQCONV2, PYGPERIODICT
3009:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM3009:      USE COMMONS, ONLY: ATMASS, NINTS, ZSYM
3010:      USE PORFUNCS3010:      USE PORFUNCS
3011:      USE GENRIGID3011:      USE GENRIGID
3012:      IMPLICIT NONE3012:      IMPLICIT NONE
3013: 3013: 
3014:      CHARACTER(LEN=20) :: PINFOSTRING3014:      CHARACTER(LEN=20) :: PINFOSTRING
3015:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(NOPT),RMS,DIAG(NOPT),TEMPA(9*NATOMS),DUMQ(NOPT)3015:      DOUBLE PRECISION :: DIHE,ALLANG,DISTPF,DUMMY1,GRAD(NOPT),RMS,DIAG(NOPT),TEMPA(9*NATOMS),DUMQ(NOPT)
3016:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT3016:      INTEGER :: HORDER,INFO,J2,K1,RECLEN,ISTAT,J1,LUNIT,GETUNIT
3017:      LOGICAL :: BTEST,KD,NNZ,NINTB,TSFRQDONE,MINFRQDONE3017:      LOGICAL :: BTEST,KD,NNZ,NINTB,TSFRQDONE,MINFRQDONE
3018:      DOUBLE PRECISION :: QTS(NOPT),QPLUS(NOPT),QMINUS(NOPT),FRQSTS(NOPT),FRQSPLUS(NOPT),FRQSMINUS(NOPT), &3018:      DOUBLE PRECISION :: QTS(NOPT),QPLUS(NOPT),QMINUS(NOPT),FRQSTS(NOPT),FRQSPLUS(NOPT),FRQSMINUS(NOPT), &
3260:            ENDIF3260:            ENDIF
3261:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP3261:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP
3262:            IF (.NOT.NOFRQS) THEN3262:            IF (.NOT.NOFRQS) THEN
3263:               IF (RIGIDINIT) THEN3263:               IF (RIGIDINIT) THEN
3264:                  CALL GENRIGID_EIGENVALUES(QPLUS, ATMASS, DIAG, INFO)3264:                  CALL GENRIGID_EIGENVALUES(QPLUS, ATMASS, DIAG, INFO)
3265: 3265: 
3266:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN3266:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
3267:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)3267:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
3268:                  ENDIF3268:                  ENDIF
3269: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame3269: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
3270:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN3270:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
3271:                  RBAANORMALMODET = .TRUE.3271:                  RBAANORMALMODET = .TRUE.
3272:                  CALL POTENTIAL(QPLUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)3272:                  CALL POTENTIAL(QPLUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
3273:                  CALL NRMLMD (QPLUS, DIAG, .FALSE.)3273:                  CALL NRMLMD (QPLUS, DIAG, .FALSE.)
3274:                  RBAANORMALMODET = .FALSE.3274:                  RBAANORMALMODET = .FALSE.
3275:               ELSE3275:               ELSE
3276:                  IF (ENDNUMHESS) THEN3276:                  IF (ENDNUMHESS) THEN
3277:                     CALL MAKENUMHESS(QPLUS,NATOMS)3277:                     CALL MAKENUMHESS(QPLUS,NATOMS)
3278:                  ELSE3278:                  ELSE
3279:                     CALL POTENTIAL(QPLUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)3279:                     CALL POTENTIAL(QPLUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
3280:                  ENDIF3280:                  ENDIF
3613:            IF (.NOT.NOFRQS) THEN3613:            IF (.NOT.NOFRQS) THEN
3614:               IF (RIGIDINIT) THEN3614:               IF (RIGIDINIT) THEN
3615: ! hk286 - TS is recorded in rigid body coordinates3615: ! hk286 - TS is recorded in rigid body coordinates
3616:                  ATOMRIGIDCOORDT = .FALSE.3616:                  ATOMRIGIDCOORDT = .FALSE.
3617:                  CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO)3617:                  CALL GENRIGID_EIGENVALUES(QTS, ATMASS, DIAG, INFO)
3618:                  ATOMRIGIDCOORDT = .TRUE.3618:                  ATOMRIGIDCOORDT = .TRUE.
3619:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN3619:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
3620:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)3620:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
3621:                  ENDIF3621:                  ENDIF
3622: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame3622: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
3623:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN3623:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
3624:                  RBAANORMALMODET = .TRUE.3624:                  RBAANORMALMODET = .TRUE.
3625:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)3625:                  CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
3626:                  CALL NRMLMD (QTS, DIAG, .FALSE.)3626:                  CALL NRMLMD (QTS, DIAG, .FALSE.)
3627:                  RBAANORMALMODET = .FALSE.3627:                  RBAANORMALMODET = .FALSE.
3628:               ELSE3628:               ELSE
3629:                  IF (ENDNUMHESS) THEN3629:                  IF (ENDNUMHESS) THEN
3630:                     CALL MAKENUMHESS(QTS,NATOMS)3630:                     CALL MAKENUMHESS(QTS,NATOMS)
3631:                  ELSE3631:                  ELSE
3632:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)3632:                     CALL POTENTIAL(QTS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
3633:                  ENDIF3633:                  ENDIF
3961:            ENDIF3961:            ENDIF
3962:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP3962:            WRITE(88,'(I6,1X,A4)') HORDER,FPGRP
3963:            IF (.NOT.NOFRQS) THEN3963:            IF (.NOT.NOFRQS) THEN
3964:               IF (RIGIDINIT) THEN3964:               IF (RIGIDINIT) THEN
3965:                  CALL GENRIGID_EIGENVALUES(QMINUS, ATMASS, DIAG, INFO)3965:                  CALL GENRIGID_EIGENVALUES(QMINUS, ATMASS, DIAG, INFO)
3966: 3966: 
3967:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN3967:                  IF (DIAG(1).LT.DIAG(DEGFREEDOMS)) THEN
3968:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)3968:                     CALL EIGENSORT_VAL_ASC(DIAG(1:DEGFREEDOMS),HESS(1:DEGFREEDOMS,1:DEGFREEDOMS),DEGFREEDOMS,DEGFREEDOMS)
3969:                  ENDIF3969:                  ENDIF
3970: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame3970: ! hk286 - compute normal modes for rigid body angle axis, go to moving frame
3971:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN3971:               ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
3972:                  RBAANORMALMODET = .TRUE.3972:                  RBAANORMALMODET = .TRUE.
3973:                  CALL POTENTIAL(QMINUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)3973:                  CALL POTENTIAL(QMINUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
3974:                  CALL NRMLMD (QMINUS, DIAG, .FALSE.)3974:                  CALL NRMLMD (QMINUS, DIAG, .FALSE.)
3975:                  RBAANORMALMODET = .FALSE.3975:                  RBAANORMALMODET = .FALSE.
3976:               ELSE3976:               ELSE
3977:                  IF (ENDNUMHESS) THEN3977:                  IF (ENDNUMHESS) THEN
3978:                     CALL MAKENUMHESS(QMINUS,NATOMS)3978:                     CALL MAKENUMHESS(QMINUS,NATOMS)
3979:                  ELSE3979:                  ELSE
3980:                     CALL POTENTIAL(QMINUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)3980:                     CALL POTENTIAL(QMINUS,DUMMY1,GRAD,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
3981:                  ENDIF3981:                  ENDIF


r31955/path.f 2017-02-22 08:30:09.202504330 +0000 r31954/path.f 2017-02-22 08:30:11.354533393 +0000
1892:          WRITE(88,'(I6,1X,A4)') HORDER,FPGRP1892:          WRITE(88,'(I6,1X,A4)') HORDER,FPGRP
1893: C        CALL H2OMODES(NATOMS,IPOT,Q,DIAG) ! WCOMMENT1893: C        CALL H2OMODES(NATOMS,IPOT,Q,DIAG) ! WCOMMENT
1894:          CALL H2OMODES(NATOMS/2,IPOT,Q,DIAG)1894:          CALL H2OMODES(NATOMS/2,IPOT,Q,DIAG)
1895: C        WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,6*NATOMS) ! WCOMMENT1895: C        WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,6*NATOMS) ! WCOMMENT
1896:          IF (.NOT.NOFRQS) WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS)1896:          IF (.NOT.NOFRQS) WRITE(88,'(3G20.10)') (DIAG(J2),J2=1,3*NATOMS)
1897: 1897: 
1898:       ! sn402: For the small number of potentials coded in the fully-rigid form (RBAAT), we can use the NRMLMD subroutine1898:       ! sn402: For the small number of potentials coded in the fully-rigid form (RBAAT), we can use the NRMLMD subroutine
1899:       ! directly to calculate the Hessian and calculate its eigenvalues. We write squared angular frequencies in whatever1899:       ! directly to calculate the Hessian and calculate its eigenvalues. We write squared angular frequencies in whatever
1900:       ! time unit is specified by FRQCONV (see documentation and comments in keywords.f)1900:       ! time unit is specified by FRQCONV (see documentation and comments in keywords.f)
1901: ! hk286 - compute potential for normal modes, notice the toggle1901: ! hk286 - compute potential for normal modes, notice the toggle
1902:       ELSE IF (RBAAT.AND..NOT.PYGPERIODICT.AND..NOT.SANDBOXT) THEN1902:       ELSE IF (RBAAT.AND..NOT.PYGPERIODICT) THEN
1903:          CALL POTENTIAL(Q,E,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)1903:          CALL POTENTIAL(Q,E,VNEW,.TRUE.,.TRUE.,RMS,.FALSE.,.FALSE.)
1904:          WRITE(88, '(G20.10)') E1904:          WRITE(88, '(G20.10)') E
1905:          ! The following line should definitely not be hard-coded in!1905:          ! The following line should definitely not be hard-coded in!
1906:          WRITE(88, '(A,A)') "1  ", "C1" ! TEMP Solution1906:          WRITE(88, '(A,A)') "1  ", "C1" ! TEMP Solution
1907: 1907: 
1908:          IF (.NOT.NOFRQS) THEN1908:          IF (.NOT.NOFRQS) THEN
1909:             RBAANORMALMODET = .TRUE.           1909:             RBAANORMALMODET = .TRUE.           
1910:             CALL NRMLMD (Q, DIAG, .FALSE.)1910:             CALL NRMLMD (Q, DIAG, .FALSE.)
1911:             RBAANORMALMODET = .FALSE.1911:             RBAANORMALMODET = .FALSE.
1912:             WRITE(88,'(3G20.10)') (DIAG(J2)*FRQCONV2,J2=1,3*NATOMS)  ! sn402: Added the frequency conversion1912:             WRITE(88,'(3G20.10)') (DIAG(J2)*FRQCONV2,J2=1,3*NATOMS)  ! sn402: Added the frequency conversion


r31955/potential.f 2017-02-22 08:30:09.434507462 +0000 r31954/potential.f 2017-02-22 08:30:11.586536525 +0000
1866:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1866:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1867:             END IF1867:             END IF
1868: 1868: 
1869: 1869: 
1870:          ELSE IF (STOCKAAT) THEN1870:          ELSE IF (STOCKAAT) THEN
1871:             CALL STOCKGHAA(COORDS,VNEW,ENERGY,GTEST,SSTEST)1871:             CALL STOCKGHAA(COORDS,VNEW,ENERGY,GTEST,SSTEST)
1872:             IF (PTEST) THEN1872:             IF (PTEST) THEN
1873:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1873:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1874:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1874:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1875:             ENDIF1875:             ENDIF
1876:          ELSE IF (SANDBOXT) THEN 
1877:             CALL SANDBOX (COORDS,VNEW,ENERGY,GTEST,SSTEST) 
1878:             IF (SSTEST) THEN 
1879:                CALL SANDBOXSECDER(COORDS,SSTEST) 
1880:             END IF 
1881:  
1882:             IF (PTEST) THEN 
1883:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         ' 
1884:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         ' 
1885:             ENDIF 
1886:         
1887:          ELSE IF (SILANET) THEN1876:          ELSE IF (SILANET) THEN
1888:             CALL SILANE(COORDS,VNEW,ENERGY,GTEST,SSTEST)1877:             CALL SILANE(COORDS,VNEW,ENERGY,GTEST,SSTEST)
1889:             IF (PTEST) THEN1878:             IF (PTEST) THEN
1890:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1879:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1891:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1880:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1892:             END IF1881:             END IF
1893: 1882: 
1894:          ELSE IF (VASP) THEN1883:          ELSE IF (VASP) THEN
1895:             FNAME='POSCAR'1884:             FNAME='POSCAR'
1896:             OPEN(UNIT=7,FILE=FNAME,STATUS='OLD')1885:             OPEN(UNIT=7,FILE=FNAME,STATUS='OLD')


r31955/rbperm.f90 2017-02-22 08:30:09.698511031 +0000 r31954/rbperm.f90 2017-02-22 08:30:12.186544665 +0000
 32: !     centre of coordinates of COORDSA will be the same as for COORDSB. 32: !     centre of coordinates of COORDSA will be the same as for COORDSB.
 33: ! 33: !
 34: !     ---------------------------------------------------------------------------------------------- 34: !     ----------------------------------------------------------------------------------------------
 35:   35:  
 36:       SUBROUTINE RBMINPERMDIST(COORDSB,COORDSA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,SITESB,SITESA) 36:       SUBROUTINE RBMINPERMDIST(COORDSB,COORDSA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,SITESB,SITESA)
 37:  37: 
 38: !     returns DISTANCE as the actual distance, rather than the squared distance 38: !     returns DISTANCE as the actual distance, rather than the squared distance
 39:  39: 
 40:       USE COMMONS, ONLY : NATOMS, NRBSITES  40:       USE COMMONS, ONLY : NATOMS, NRBSITES 
 41:       USE KEY,ONLY : NTSITES, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, GEOMDIFFTOL, BESTPERM, EFIELDT, PAPT, PTSTSTT, & 41:       USE KEY,ONLY : NTSITES, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, GEOMDIFFTOL, BESTPERM, EFIELDT, PAPT, PTSTSTT, &
 42:      &               RBOPS, NRBGROUP, RBSYMT, DBPTDT, MULTISITEPYT, PYGPERIODICT, SANDBOXT 42:      &               RBOPS, NRBGROUP, RBSYMT, DBPTDT, MULTISITEPYT, PYGPERIODICT 
 43:  43: 
 44:       IMPLICIT NONE 44:       IMPLICIT NONE
 45:  45: 
 46:       INTEGER, PARAMETER :: MAXIMUMTRIES = 100 46:       INTEGER, PARAMETER :: MAXIMUMTRIES = 100
 47:       INTEGER            :: NPERM, PATOMS, NTRIES, NSIZE, J1, J2, J3, NRB !jwrm2> unused, INFO, JMAX, LOCMAX(1) 47:       INTEGER            :: NPERM, PATOMS, NTRIES, NSIZE, J1, J2, J3, NRB !jwrm2> unused, INFO, JMAX, LOCMAX(1)
 48:       INTEGER            :: INVERT, NORBIT1, NORBIT2, PERM(NATOMS), NCHOOSE2, NDUMMY, LPERM(NATOMS), NCHOOSE1 48:       INTEGER            :: INVERT, NORBIT1, NORBIT2, PERM(NATOMS), NCHOOSE2, NDUMMY, LPERM(NATOMS), NCHOOSE1
 49:       INTEGER            :: NEWPERM(NATOMS), ALLPERM(NATOMS) 49:       INTEGER            :: NEWPERM(NATOMS), ALLPERM(NATOMS)
 50:       DOUBLE PRECISION   :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DIST2 !jwrm2> unused, DISTWP, TEMPA(9*NATOMS)   50:       DOUBLE PRECISION   :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), DISTANCE, DIST2 !jwrm2> unused, DISTWP, TEMPA(9*NATOMS)  
 51:       DOUBLE PRECISION   :: DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) !jwrm2> unused, DUMMYWP(3*NATOMS) 51:       DOUBLE PRECISION   :: DUMMYA(3*NATOMS), DUMMYB(3*NATOMS), DUMMY(3*NATOMS) !jwrm2> unused, DUMMYWP(3*NATOMS)
 52:       DOUBLE PRECISION   :: XA(3*NTSITES),  XB(3*NTSITES), XTMP(3*NATOMS) !jwrm2> unused, XBS(3*NTSITES) 52:       DOUBLE PRECISION   :: XA(3*NTSITES),  XB(3*NTSITES), XTMP(3*NATOMS) !jwrm2> unused, XBS(3*NTSITES)
 68:       COORDSCOMA(:) = COORDSA(1:3*NATOMS/2) 68:       COORDSCOMA(:) = COORDSA(1:3*NATOMS/2)
 69:       COORDSCOMB(:) = COORDSB(1:3*NATOMS/2) 69:       COORDSCOMB(:) = COORDSB(1:3*NATOMS/2)
 70:  70: 
 71:       IF (DEBUG) THEN 71:       IF (DEBUG) THEN
 72:          WRITE(*,'(3F20.10)')COORDSA(1:3*NATOMS) 72:          WRITE(*,'(3F20.10)')COORDSA(1:3*NATOMS)
 73:          WRITE(*,*) 73:          WRITE(*,*)
 74:          WRITE(*,'(3F20.10)')COORDSB(1:3*NATOMS) 74:          WRITE(*,'(3F20.10)')COORDSB(1:3*NATOMS)
 75:       ENDIF 75:       ENDIF
 76:  76: 
 77: !     do not call minpermdistrbCOM for multisite PY  77: !     do not call minpermdistrbCOM for multisite PY 
 78:       IF(MULTISITEPYT.OR.PYGPERIODICT.OR.SANDBOXT) GOTO 100 78:       IF(MULTISITEPYT.OR.PYGPERIODICT) GOTO 100
 79:         NATOMS = NATOMS/2 79:         NATOMS = NATOMS/2
 80:  80: 
 81:         CALL MINPERMDISTRBCOM(COORDSCOMB,COORDSCOMA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT) 81:         CALL MINPERMDISTRBCOM(COORDSCOMB,COORDSCOMA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT)
 82:  82: 
 83:         NATOMS = 2*NATOMS 83:         NATOMS = 2*NATOMS
 84:  84: 
 85:         IF (SQRT(DISTANCE) <= GEOMDIFFTOL) THEN 85:         IF (SQRT(DISTANCE) <= GEOMDIFFTOL) THEN
 86:          DISTANCE = SQRT(DISTANCE) 86:          DISTANCE = SQRT(DISTANCE)
 87:          IF (DEBUG) PRINT '(A)',' rbpermdist> minpermdistrbcom suggests identical' 87:          IF (DEBUG) PRINT '(A)',' rbpermdist> minpermdistrbcom suggests identical'
 88:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0 88:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0
117: 100   CONTINUE117: 100   CONTINUE
118: 118: 
119:       COORDSBS(1:3*NATOMS) = COORDSB(1:3*NATOMS) ! to trace error, see at the end119:       COORDSBS(1:3*NATOMS) = COORDSB(1:3*NATOMS) ! to trace error, see at the end
120:       COORDSAS(1:3*NATOMS) = COORDSA(1:3*NATOMS) ! to trace error, see at the end120:       COORDSAS(1:3*NATOMS) = COORDSA(1:3*NATOMS) ! to trace error, see at the end
121: !      print *, 'COORDSB'121: !      print *, 'COORDSB'
122: !      print *, COORDSB122: !      print *, COORDSB
123: !      print *, 'COORDSA'123: !      print *, 'COORDSA'
124: !      print *, COORDSA124: !      print *, COORDSA
125: 125: 
126:       NRB   = (NATOMS/2)126:       NRB   = (NATOMS/2)
127:       IF (DBPTDT.OR.SANDBOXT) THEN127:       IF (DBPTDT) THEN
128:          NSIZE = NTSITES128:          NSIZE = NTSITES
129:       ELSE 129:       ELSE
130:          NSIZE = NRB*NRBSITES130:          NSIZE = NRB*NRBSITES
131:       ENDIF131:       ENDIF
 132: 
132: !      CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0133: !      CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0
133: !      DO J1 = 1, NATOMS/2134: !      DO J1 = 1, NATOMS/2
134: !         CMAX = CMAX + COORDSA(3*(J1-1)+1)135: !         CMAX = CMAX + COORDSA(3*(J1-1)+1)
135: !         CMAY = CMAY + COORDSA(3*(J1-1)+2)136: !         CMAY = CMAY + COORDSA(3*(J1-1)+2)
136: !         CMAZ = CMAZ + COORDSA(3*(J1-1)+3)137: !         CMAZ = CMAZ + COORDSA(3*(J1-1)+3)
137: !       ENDDO138: !       ENDDO
138:       CMAX = 2*CMAX/NATOMS; CMAY = 2*CMAY/NATOMS; CMAZ = 2*CMAZ/NATOMS139:       CMAX = 2*CMAX/NATOMS; CMAY = 2*CMAY/NATOMS; CMAZ = 2*CMAZ/NATOMS
139:       CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0140:       CMBX = 0.0D0; CMBY = 0.0D0; CMBZ = 0.0D0
140:       DO J1 = 1, NATOMS/2141:       DO J1 = 1, NATOMS/2
141:          CMBX = CMBX + COORDSB(3*(J1-1)+1)142:          CMBX = CMBX + COORDSB(3*(J1-1)+1)
649: 650: 
650:       END SUBROUTINE RBMINPERMDIST651:       END SUBROUTINE RBMINPERMDIST
651: 652: 
652: !     ----------------------------------------------------------------------------------------------653: !     ----------------------------------------------------------------------------------------------
653: 654: 
654:       SUBROUTINE SITEPOS(Q, X)655:       SUBROUTINE SITEPOS(Q, X)
655: 656: 
656: !     GENERATE SITE POSITIONS FROM THE RIGID-BODY COORDINATES657: !     GENERATE SITE POSITIONS FROM THE RIGID-BODY COORDINATES
657: 658: 
658:       USE COMMONS, ONLY: NATOMS, NRBSITES, RBSITE659:       USE COMMONS, ONLY: NATOMS, NRBSITES, RBSITE
659:       USE KEY, ONLY: NTSITES, DBPTDT, SANDBOXT660:       USE KEY, ONLY: NTSITES, DBPTDT
660:       use sandbox_module 
661: 661: 
662:       IMPLICIT NONE662:       IMPLICIT NONE
663: 663: 
664:       INTEGER :: J1, J2, J3, J4664:       INTEGER :: J1, J2, J3, J4
665:       DOUBLE PRECISION :: Q(3*NATOMS), X(3*NTSITES), R(3), P(3), RM(3,3)665:       DOUBLE PRECISION :: Q(3*NATOMS), X(3*NTSITES), R(3), P(3), RM(3,3)
666:       DOUBLE PRECISION :: RBSITETD(4,3), SIGTD(4), SIGDB(2)666:       DOUBLE PRECISION :: RBSITETD(4,3), SIGTD(4), SIGDB(2)
667:       type(sandbox_molecule), pointer :: mol 
668:       type(sandbox_site), pointer :: site 
669:       integer :: atom_num, atom_index, site_num, mol_num 
670:       double precision :: dummy(3, 3) 
671: 667: 
672:       IF (DBPTDT) CALL DEFTD(RBSITETD,SIGTD,SIGDB)668:       IF (DBPTDT) CALL DEFTD(RBSITETD,SIGTD,SIGDB)
673: 669: 
674:       IF (SANDBOXT) THEN 
675:         J1=0     
676:         do mol_num = 1, size(molecules) 
677:             mol => molecules(mol_num) 
678:  
679:             call rmdrvt(q(mol%p_index: mol%p_index + 2), RM, dummy, dummy, dummy, .false.) 
680: !            WRITE(*,*) 'size(molecules)=', size(molecules(mol_num)%sites) 
681:  
682:             call update_sandbox_molecule(.false., mol, q(mol%r_index: mol%r_index + 2),q(mol%p_index: mol%p_index + 2)) 
683:  
684:             do atom_num=1,size(molecules(mol_num)%sites) 
685:                J1=J1+1 
686:                X(3*J1-2:3*J1)=molecules(mol_num)%sites(atom_num)%r 
687:             enddo 
688:         end do 
689: !        WRITE(*,*) X(:) 
690:         RETURN  
691:  
692:       END IF 
693:  
694:       DO J1 = 1, NATOMS/2670:       DO J1 = 1, NATOMS/2
695: 671: 
696:          J2   = 3*J1672:          J2   = 3*J1
697:          R(:) = Q(J2-2:J2)673:          R(:) = Q(J2-2:J2)
698:          J2   = 3*NATOMS/2 + J2674:          J2   = 3*NATOMS/2 + J2
699:          P(:) = Q(J2-2:J2)675:          P(:) = Q(J2-2:J2)
700: 676: 
701:          CALL ROTMAT(P, RM)677:          CALL ROTMAT(P, RM)
702: 678: 
703:          IF (DBPTDT) THEN679:          IF (DBPTDT) THEN


r31955/rigidb.f90 2017-02-22 08:30:09.962514611 +0000 r31954/rigidb.f90 2017-02-22 08:30:12.410547654 +0000
622:                EV(J1,NZERO+I)  = 1.D0 622:                EV(J1,NZERO+I)  = 1.D0 
623:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2623:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2
624:             624:             
625:             ELSE 625:             ELSE 
626:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH))626:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH))
627:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))627:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH))
628:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 &628:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 &
629:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH)629:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH)
630:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2  630:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2  
631:             ENDIF631:             ENDIF
632:          ELSE IF (SANDBOXT.AND.UNIAXARRAY(I)) THEN   !sandbox potential with arbitrary indices of rotationally symmetric rigid bodies632:        
633:             IF (THETA == 0.D0) THEN 
634:  
635:                EV(J1,NZERO+I)  = 1.D0 
636:                NRMFCT(NZERO+I) = NRMFCT(NZERO+I) + EV(J1,NZERO+I)**2 
637:  
638:             ELSE 
639:                EV(J1-2,NZERO+I) = 0.5D0*Q(J1-1) + Q(J1-2)*Q(J1)/THETA2 - 0.5D0*Q(J1-2)*Q(J1)/(THETA*TAN(THETAH)) 
640:                EV(J1-1,NZERO+I) = - 0.5D0*Q(J1-2) + Q(J1-1)*Q(J1)/THETA2 - 0.5D0*Q(J1-1)*Q(J1)/(THETA*TAN(THETAH)) 
641:                EV(J1,NZERO+I)   = THETAH*SIN(THETA) + Q(J1)*Q(J1)/THETA2 & 
642:                                 + 0.5D0*(THETA*COS(THETA)-Q(J1)*Q(J1)/THETA)/TAN(THETAH) 
643:                NRMFCT(NZERO+I)  = NRMFCT(NZERO+I) + EV(J1-2,NZERO+I)**2 + EV(J1-1,NZERO+I)**2 + EV(J1,NZERO+I)**2 
644:             ENDIF 
645:          ENDIF633:          ENDIF
646: 634: 
647:       ENDDO635:       ENDDO
648: ! sf344> set the number of zeros properly for SANDBOX, if only some of the building blocks have rotational symmetry636: 
649:       IF(UNIAXT) THEN  
650:          DO I=1,NATOMS/2 
651:            IF (SANDBOXT.AND.UNIAXARRAY(I)) THEN 
652:              NZERO = NZERO + 1 
653:            END IF 
654:          END DO 
655:       END IF 
656:       IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN637:       IF (GBT.OR.GBDT.OR.(PYGT.AND.UNIAXT).OR.STOCKAAT.OR.((PYT.OR.PYGPERIODICT.OR.PYBINARYT.OR.MULTISITEPYT).AND.UNIAXT)) THEN
657:          NZERO = NZERO + NATOMS/2638:          NZERO = NZERO + NATOMS/2
658:       ELSE639:       ELSE
659:          NZERO = NZERO640:          NZERO = NZERO
660:       ENDIF641:       ENDIF
661: 642: 
662:       DO J = 1, NZERO643:       DO J = 1, NZERO
663: 644: 
664:          NRMFCT(J) = DSQRT(NRMFCT(J))645:          NRMFCT(J) = DSQRT(NRMFCT(J))
665:          EV(:,J)   = EV(:,J)/NRMFCT(J)646:          EV(:,J)   = EV(:,J)/NRMFCT(J)


r31955/sandbox.f90 2017-02-22 08:30:10.214517994 +0000 r31954/sandbox.f90 2017-02-22 08:30:12.638550735 +0000
  1: module sandbox_module  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/sandbox.f90' in revision 31954
  2:     implicit none 
  3:  
  4:     type sandbox_site 
  5:         type(sandbox_molecule), pointer :: molecule 
  6:         integer :: class    ! interaction class 
  7:         double precision :: r_molframe(3), p_molframe(3), rotmat_molframe(3, 3) ! molecule frame position and orientation 
  8:         double precision :: r(3), p(3), rotmat(3, 3)  ! lab frame position and orientation 
  9:         double precision :: rotmat_deriv(3, 3, 3)   ! dR^i/dp^I = deriv with resp to molecule angle-axis 
 10:  
 11:         logical :: is_aniso 
 12:  
 13:         logical :: has_pole 
 14:         double precision :: mu_hat0(3), mu_hat(3), mu_hat_deriv(3, 3)   ! mu_hat derivatives with resp to p^i 
 15:         character(len=20):: name 
 16:     end type sandbox_site 
 17:  
 18:     type sandbox_molecule 
 19:         type(sandbox_site), allocatable :: sites(:) 
 20:         double precision :: r(3), p(3) 
 21:         double precision :: rotmat(3, 3), rotmat_deriv(3, 3, 3) 
 22:         integer :: r_index, p_index 
 23:         character :: atom_symbol 
 24:     end type sandbox_molecule 
 25:      
 26:     type interaction 
 27:         logical :: exists 
 28:         integer :: func_id  ! see table below for potential <-> func_id correspondence 
 29:         double precision, allocatable :: params(:) 
 30:     end type interaction 
 31:  
 32:     type interaction_class 
 33:         logical :: has_pole, is_aniso 
 34:     end type interaction_class 
 35:  
 36:     double precision, parameter :: z_hat(3) = (/ 0.D0, 0.D0, 1.D0 /) 
 37:  
 38: ! inserts for compatibility with OPTIM 
 39:     double precision, allocatable :: vt(:), coords(:) 
 40:     double precision :: step(1), ostep(1), astep(1), perccut, radius 
 41:     logical :: mpit, modulart, percolatet, omove(1), tmove(1) 
 42:     integer :: nsave, mynode 
 43:  
 44:     type(sandbox_molecule), allocatable, target :: molecules(:) 
 45:     type(interaction), allocatable, target :: interactions(:, :) 
 46:     type(interaction_class), allocatable :: classes(:) 
 47:  
 48:     integer :: num_mols, x_size, num_atoms   ! length of x array 
 49:  
 50:     ! to add a new potential, add a sandbox_[potential name] subroutine, update pairwise_sandbox, and update 
 51:     ! sandbox_input 
 52:  
 53:     ! table of potentials and their integer IDs func_id: 
 54:     ! lj 1 
 55:     ! coulomb 2 
 56:     ! dipole 3 
 57:     ! chiro 4 
 58:     ! morse 5 
 59:     ! coulomb with distance dependent dielectric 6 
 60:  
 61:  
 62: contains 
 63:     subroutine initialize_sandbox_molecule(mol, mol_index) 
 64:         implicit none 
 65:         type(sandbox_molecule), target, intent(inout) :: mol 
 66:         integer, intent(in) :: mol_index 
 67:  
 68:         integer :: site_index 
 69:         type(sandbox_site), pointer :: site 
 70:         double precision :: dummy(3, 3) 
 71:  
 72:         ! initialize molecule values 
 73:         mol%r_index = 3 * mol_index - 2 
 74:         mol%p_index = mol%r_index + 3 * num_mols 
 75:  
 76:         do site_index = 1, size(mol%sites) 
 77:             site => mol%sites(site_index) 
 78:             site%molecule => mol    ! give all sites a pointer to their parent 
 79:  
 80:             ! if the site is anisotropic, need to do some extra bookkeeping 
 81:             if(classes(site%class)%is_aniso) then 
 82:                 site%is_aniso = .true. 
 83:  
 84:                 ! compute rotation matrix for site 
 85:                 call rmdrvt(site%p_molframe, site%rotmat_molframe, dummy, dummy, dummy, .false.) 
 86:  
 87:                 ! if site has a pole, turn on a marker to keep track of mu_hat 
 88:                 if(classes(site%class)%has_pole) then 
 89:                     site%has_pole = .true. 
 90:                     site%mu_hat0 = matmul(site%rotmat_molframe, z_hat) 
 91:                 else 
 92:                     site%has_pole = .false. 
 93:                 end if 
 94:             end if 
 95:         end do 
 96:     end subroutine initialize_sandbox_molecule 
 97:  
 98:     subroutine update_sandbox_molecule(gtest, mol, r, p) 
 99:         implicit none 
100:         logical, intent(in) :: gtest 
101:         type(sandbox_molecule), target, intent(inout) :: mol 
102:         double precision, intent(inout) :: r(3), p(3) 
103:  
104:         double precision, parameter :: pi = 3.14159265359 
105:         double precision :: pmod 
106:         type(sandbox_site), pointer :: site 
107:         integer :: site_index, k 
108:         double precision :: dummy(3, 3), dR_dp(3, 3, 3) 
109:  
110:         double precision :: swo(4, 3, 3) 
111:  
112:         mol%r(:) = r(:) 
113:  
114:         ! make sure that 0 < |p| < 2*pi 
115:         pmod = sqrt(dot_product(p, p)) 
116:         if(pmod > 2 * pi) p(:) = p(:) / pmod * mod(pmod, 2 * pi) 
117:         mol%p(:) = p(:) 
118:  
119:         ! recompute the rotation matrices and derivatives 
120:         call rmdrvt(mol%p, mol%rotmat, & 
121:             & mol%rotmat_deriv(1, :, :), mol%rotmat_deriv(2, :, :), mol%rotmat_deriv(3, :, :), gtest) 
122:  
123:         ! loop over sites 
124:         do site_index = 1, size(mol%sites) 
125:             site => mol%sites(site_index) 
126:  
127:             site%r(:) = mol%r(:) + matmul(mol%rotmat, site%r_molframe(:))   ! reset site position 
128:  
129:             if(site%is_aniso) then 
130:                 ! if site has a pole, update the pole 
131:                 if(site%has_pole) then 
132:                     ! mu_hat = R^I . mu_hat0 = R^I . R^i0 . z_hat 
133:                     site%mu_hat = matmul(mol%rotmat, site%mu_hat0) 
134:  
135:                     ! dmu/dp^I_k = dR^I/dp^I_k . mu_hat0 
136:                     do k = 1, 3 
137:                         site%mu_hat_deriv(k, :) = matmul(mol%rotmat_deriv(k, :, :), site%mu_hat0) 
138:                     end do 
139:                 end if 
140:             end if 
141:         end do 
142:  
143:     end subroutine update_sandbox_molecule 
144:  
145:     subroutine pairwise_sandbox(gtest, sitei, sitej, energy_contrib, grad_contrib) 
146:         use rotations 
147:         use vec3 
148:         implicit none 
149:         logical, intent(in) :: gtest 
150:         type(sandbox_site), intent(in) :: sitei, sitej 
151:         double precision, intent(out) :: energy_contrib, grad_contrib(12) 
152:  
153:         integer :: classi, classj 
154:         type(interaction), pointer :: this_interaction 
155:         integer, pointer :: func_id 
156:         double precision, pointer :: params(:) 
157:         double precision :: dr_dp(3, 3) 
158:         integer :: k 
159:  
160:         energy_contrib = 0.D0 
161:         grad_contrib(:) = 0.D0 
162:  
163:         classi = min(sitei%class, sitej%class) 
164:         classj = max(sitei%class, sitej%class) 
165:  
166:         if(interactions(classi, classj)%exists) then 
167:  
168:         func_id => interactions(classi, classj)%func_id 
169:         params => interactions(classi, classj)%params 
170:  
171:             if(func_id == 1) then 
172:                 ! lj 
173:                 ! sigma_0 eps_0 rep, att 
174:                 call sandbox_lj(gtest, params(1), params(2), params(3), params(4), sitei%r - sitej%r, energy_contrib, grad_contrib) 
175:             elseif(func_id == 2) then 
176:                 ! coulomb 
177:                 ! sigma_0 k_c*q*q 
178:                 call sandbox_coulomb(gtest, params(1), params(2), sitei%r - sitej%r, energy_contrib, grad_contrib) 
179:             elseif(func_id == 3) then 
180:                 ! dipole 
181:                 ! sigma_0 mu^2 
182:                 call sandbox_dipole(gtest, params(1), params(2), sitei%r - sitej%r, sitei%mu_hat, & 
183:                     & sitej%mu_hat, sitei%mu_hat_deriv, sitej%mu_hat_deriv, energy_contrib, grad_contrib) 
184:             elseif(func_id == 4) then 
185:                 ! chiro 
186:                 ! sigma_0 mu^2 cos_gamma sin_gamma 
187:                 call sandbox_chiro(gtest, params(1), params(2), params(3), params(4), sitei%r - sitej%r, sitei%mu_hat, & 
188:                     & sitej%mu_hat, sitei%mu_hat_deriv, sitej%mu_hat_deriv, energy_contrib, grad_contrib) 
189:             elseif(func_id == 5) then 
190:                 ! morse 
191:                 ! eps_0 beta r_equilibrium 
192:                 call sandbox_morse(gtest, params(1), params(2), params(3), & 
193:                     sitei%r - sitej%r, energy_contrib, grad_contrib) 
194:             elseif(func_id == 6) then 
195:                 ! coulomb 
196:                 ! sigma_0 k_c*q*q 
197:                 call sandbox_coulomb_ddiel(gtest, params(1), params(2), sitei%r - sitej%r, energy_contrib, grad_contrib) 
198:             else 
199:                 print *, "pairwise_sandbox> don't have potential with func_id='", func_id, & 
200:                     & "' for interaction of classes", classi, classj 
201:                 return 
202:             end if 
203:  
204:             if(gtest) then 
205:                 ! if site is not at molecule origin, need to compute dr^i/dp^I 
206:                 if(any(sitei%r_molframe /= 0.D0)) then 
207:                     do k = 1, 3 
208:                         ! dr^i/dp^I_k = R^I_k . r^i0 
209:                         dr_dp(k, :) = matmul(sitei%molecule%rotmat_deriv(k, :, :), sitei%r_molframe) 
210:                     end do 
211:  
212:                     ! assign dU/dp^I = dr^i/dp^I . dU/dr^i 
213:                     grad_contrib(7: 9) = grad_contrib(7: 9) + matmul(dr_dp, grad_contrib(1: 3)) 
214:                 end if 
215:  
216:                 ! repeat for molecule J 
217:                 if(any(sitej%r_molframe /= 0.D0)) then 
218:                     do k = 1, 3 
219:                         dr_dp(k, :) = matmul(sitej%molecule%rotmat_deriv(k, :, :), sitej%r_molframe) 
220:                     end do 
221:  
222:                     grad_contrib(10: 12) = grad_contrib(10: 12) + matmul(dr_dp, grad_contrib(4: 6)) 
223:                 end if 
224:             end if 
225:         end if 
226:  
227:     end subroutine pairwise_sandbox 
228:  
229:  
230: ! --- begin sandbox potentials 
231:     subroutine sandbox_lj(gtest, sigma_0, eps_0, rep, att, rij, energy, grad) 
232:         ! gtest: whether gradients should be calculated 
233:         ! sigma_0: length scale 
234:         ! eps_0: energy_scale 
235:         ! rep: repulsive coefficient, i.e., A in A / r^12 (for normal LJ, this is 1.0) 
236:         ! att: attractive coefficient, i.e., B in B / r^6 (for normal LJ, this is 1.0) 
237:         ! rij: r_j - r_i 
238:  
239:         logical, intent(in) :: gtest 
240:         double precision, intent(in) :: sigma_0, eps_0, rep, att, rij(3) 
241:         double precision, intent(out) :: energy, grad(12) 
242:  
243:         double precision :: rijmod, rij_hat(3) ! |rij| and normal vector pointing rij 
244:  
245:         ! r = (|rij| / sigma_0), r_pm1 = r to power minus 1 = r^-1, etc 
246:         double precision :: r_pm1, r_pm6, r_pm7, r_pm12, r_pm13 
247:  
248:         energy = 0.D0 
249:         if(gtest) grad(:) = 0.D0 
250:  
251:         rijmod = sqrt(dot_product(rij, rij)) 
252:  
253:         r_pm1 = sigma_0 / rijmod 
254:         r_pm6 = r_pm1 * r_pm1 * r_pm1 * r_pm1 * r_pm1 * r_pm1  
255:         r_pm12 = r_pm6 * r_pm6 
256:  
257:         energy = 4.0 * eps_0 * (rep * r_pm12 - att * r_pm6) 
258:  
259:         if(gtest) then 
260:             rij_hat(:) = rij / rijmod 
261:             r_pm7 = r_pm6 * r_pm1 
262:             r_pm13 = r_pm12 * r_pm1 
263:  
264:             grad(1: 3) = 24.0 * eps_0 / sigma_0 * (-2.0 * rep * r_pm13 + att * r_pm7) * rij_hat 
265:             grad(4: 6) = -grad(1: 3) 
266:             grad(7: 9) = 0.D0 
267:             grad(10: 12) = 0.D0 
268:         end if 
269:  
270:     end subroutine sandbox_lj 
271:  
272:     subroutine sandbox_coulomb(gtest, sigma_0, kqq, rij, energy, grad) 
273:         ! gtest: if gradients are computed 
274:         ! sigma_0: length scale for rij 
275:         ! kqq: Coulomb's constant * charge_1 * charge_2 
276:         ! rij: separation between sites 
277:  
278:         ! energy: pairwise energy 
279:         ! gradient: (1:3) = radial for site i, (4:6) = radial j, (7:9) = angular for i, (10:12), angular j 
280:  
281:         logical, intent(in) :: gtest 
282:         double precision, intent(in) :: sigma_0, kqq, rij(3) 
283:         double precision, intent(out) :: energy, grad(12) 
284:  
285:         double precision :: rijmod, rij_hat(3) ! |rij| and normal vector pointing rij 
286:  
287:         ! r = (|rij| / sigma_0), r_pm1 = r to power minus 1 = r^-1, etc 
288:         double precision :: r_pm1, r_pm2  
289:  
290:         energy = 0.D0 
291:  
292:         if(gtest) grad(:) = 0.D0 
293:  
294:         rijmod = sqrt(dot_product(rij, rij)) 
295:  
296:         r_pm1 = sigma_0 / rijmod 
297:  
298:         energy = kqq * r_pm1 
299:  
300:         if(gtest) then 
301:             rij_hat(:) = rij / rijmod 
302:             r_pm2 = r_pm1 * r_pm1 
303:  
304:             grad(1: 3) = -kqq / sigma_0 * r_pm2 * rij_hat(:)  
305:             grad(4: 6) = -grad(1: 3) 
306:             grad(7: 9) = 0.D0 
307:             grad(10: 12) = 0.D0 
308:         end if 
309:  
310:     end subroutine sandbox_coulomb 
311:  
312:     subroutine sandbox_coulomb_ddiel(gtest, sigma_0, kqq, rij, energy, grad) 
313:         ! gtest: if gradients are computed 
314:         ! sigma_0: length scale for rij 
315:         ! kqq: Coulomb's constant * charge_1 * charge_2 
316:         ! rij: separation between sites 
317:  
318:         ! energy: pairwise energy 
319:         ! gradient: (1:3) = radial for site i, (4:6) = radial j, (7:9) = angular for i, (10:12), angular j 
320:  
321:         logical, intent(in) :: gtest 
322:         double precision, intent(in) :: sigma_0, kqq, rij(3) 
323:         double precision, intent(out) :: energy, grad(12) 
324:  
325:         double precision :: rijmod, rij_hat(3) ! |rij| and normal vector pointing rij 
326:         double precision :: eps, epsilon_r ! epsilon_r is a distance dependent dielectric (eps*rij) 
327:  
328:  
329:         ! r = (|rij| / sigma_0), r_pm1 = r to power minus 1 = r^-1, etc 
330:         double precision :: r_pm1, r_pm2, eps_r_pm1 
331:  
332:         energy = 0.D0 
333:  
334:         eps = 4.D0 
335:  
336:         if(gtest) grad(:) = 0.D0 
337:  
338:         rijmod = sqrt(dot_product(rij, rij)) 
339:  
340:         r_pm1 = 1.0 / rijmod 
341:  
342: ! distance dependent dielectric 
343:         epsilon_r = eps * rijmod 
344:  
345:         eps_r_pm1 = 1 / epsilon_r 
346:  
347:         energy = kqq * r_pm1 * r_pm1 * 1/eps 
348:  
349:         if(gtest) then 
350:             rij_hat(:) = rij / rijmod 
351:             r_pm2 = r_pm1 * r_pm1 
352:  
353:             grad(1: 3) = -2.0D0 * kqq * r_pm2 * rij_hat(:) * r_pm1/eps 
354:             grad(4: 6) = -grad(1: 3) 
355:             grad(7: 9) = 0.D0 
356:             grad(10: 12) = 0.D0 
357:         end if 
358:  
359:     end subroutine sandbox_coulomb_ddiel 
360:  
361:     subroutine sandbox_dipole(gtest, sigma_0, mu_p2, rij, mu_hati, mu_hatj, mu_hat_derivi, mu_hat_derivj, energy, grad) 
362:         logical, intent(in) :: gtest 
363:         double precision, intent(in) :: sigma_0, mu_p2, rij(3), mu_hati(3), mu_hatj(3), & 
364:             & mu_hat_derivi(3, 3), mu_hat_derivj(3, 3) 
365:         double precision, intent(out) :: energy, grad(12) 
366:  
367:         double precision :: rijmod, rij_hat(3) 
368:         double precision :: r_pm1, r_pm2, r_pm3, r_pm6 
369:         double precision :: mu_dot_ri, mu_dot_rj, mu_dot_mu 
370:         integer :: k 
371:  
372:         energy = 0.D0 
373:         if(gtest) grad(:) = 0.D0 
374:  
375:         rijmod = sqrt(dot_product(rij, rij)) 
376:         if(rijmod == 0.D0) then 
377:             rij_hat(:) = 0.D0 
378:         else 
379:             rij_hat(:) = rij / rijmod 
380:         end if 
381:  
382:         r_pm1 = sigma_0 / rijmod 
383:         r_pm2 = r_pm1 * r_pm1 
384:         r_pm3 = r_pm2 * r_pm1 
385:         r_pm6 = r_pm2 * r_pm2 * r_pm2 
386:  
387:         mu_dot_ri = dot_product(mu_hati, rij_hat) 
388:         mu_dot_rj = dot_product(mu_hatj, rij_hat) 
389:         mu_dot_mu = dot_product(mu_hati, mu_hatj) 
390:  
391:         energy = mu_p2 * r_pm3 * (mu_dot_mu - 3.0 * mu_dot_ri * mu_dot_rj) 
392:  
393:         if(gtest) then 
394:             grad(1: 3) = -3.0 * mu_p2 * r_pm3 / rijmod * & 
395:                 & ((mu_dot_mu - 5.0 * mu_dot_ri * mu_dot_rj) * rij_hat(:) + mu_dot_ri * mu_hatj + mu_dot_rj * mu_hati) 
396:             grad(4: 6) = -grad(1: 3) 
397:  
398:             do k = 1, 3 
399:                 grad(6 + k) = mu_p2 * r_pm3 * (dot_product(mu_hat_derivi(k, :), mu_hatj) & 
400:                     & - 3.0 * dot_product(mu_hat_derivi(k, :), rij_hat) * mu_dot_rj) 
401:                 grad(9 + k) = mu_p2 * r_pm3 * (dot_product(mu_hat_derivj(k, :), mu_hati) & 
402:                     & - 3.0 * dot_product(mu_hat_derivj(k, :), rij_hat) * mu_dot_ri) 
403:             end do 
404:         end if 
405:  
406:     end subroutine sandbox_dipole 
407:  
408:     subroutine sandbox_chiro(gtest, sigma_0, mu_p2, cos_gamma, sin_gamma, rij, mu_hati, mu_hatj, & 
409:       & mu_hat_derivi, mu_hat_derivj, energy, grad) 
410:         use vec3, only: vec_cross 
411:         implicit none 
412:         logical, intent(in) :: gtest 
413:         double precision, intent(in) :: sigma_0, mu_p2, sin_gamma, cos_gamma, rij(3), mu_hati(3), mu_hatj(3), & 
414:             & mu_hat_derivi(3, 3), mu_hat_derivj(3, 3) 
415:         double precision, intent(out) :: energy, grad(12) 
416:  
417:         double precision :: rijmod, rij_hat(3) 
418:         double precision :: r_pm1, r_pm3, r_pm4 ! r is in reduced units 
419:         double precision :: mu_dot, mu_cross(3)  ! mu_hati . mu_hatj and mu_hati x mu_hatj 
420:  
421:         integer :: k 
422:  
423:         energy = 0.D0 
424:         if(gtest) grad(:) = 0.D0 
425:  
426:         rijmod = sqrt(dot_product(rij, rij)) 
427:         if(rijmod == 0.D0) then 
428:             rij_hat(:) = 0.D0 
429:         else 
430:             rij_hat(:) = rij / rijmod 
431:         end if 
432:          
433:         r_pm1 = sigma_0 / rijmod 
434:         r_pm3 = r_pm1 * r_pm1 * r_pm1 
435:  
436:         mu_dot = dot_product(mu_hati, mu_hatj) 
437:         mu_cross(:) = vec_cross(mu_hati, mu_hatj) 
438:  
439:         energy = -mu_p2 * r_pm3 * (cos_gamma * mu_dot + sin_gamma * dot_product(mu_cross, rij_hat)) 
440:  
441:         if(gtest) then 
442:             r_pm4 = r_pm3 * r_pm1 
443:  
444:             grad(1: 3) = -mu_p2 * r_pm4 / sigma_0 * (cos_gamma * (-3.0 * mu_dot * rij_hat(:)) + & 
445:                 & sin_gamma * (mu_cross(:) - 4.0 * dot_product(mu_cross, rij_hat) * rij_hat(:))) 
446:             grad(4: 6) = -grad(1: 3) 
447:  
448:             do k = 1, 3 
449:                 grad(6 + k) = -mu_p2 * r_pm3 * (cos_gamma * dot_product(mu_hat_derivi(k, :), mu_hatj) &  
450:                     & + sin_gamma * dot_product(vec_cross(mu_hat_derivi(k, :), mu_hatj), rij_hat)) 
451:  
452:                 grad(9 + k) = -mu_p2 * r_pm3 * (cos_gamma * dot_product(mu_hati, mu_hat_derivj(k, :)) &  
453:                     & + sin_gamma * dot_product(vec_cross(mu_hati, mu_hat_derivj(k, :)), rij_hat)) 
454:             end do 
455:         end if 
456:  
457:     end subroutine sandbox_chiro 
458:  
459:     subroutine sandbox_morse(gtest, beta, requil, eps_0, rij, energy, grad) 
460:         ! beta -- range parameter 
461:         ! requil -- equilibrium bond distance 
462:         ! eps_0 -- energy scale 
463:  
464:         implicit none 
465:         logical, intent(in) :: gtest 
466:         double precision, intent(in) :: beta, requil, eps_0, rij(3) 
467:         double precision, intent(out) :: energy, grad(12) 
468:  
469:         double precision :: rijmod, rij_hat(3) 
470:         double precision :: exp_term    ! = exp{-beta * (r - re)} 
471:         double precision :: one_m_exp   ! = 1 - exp{-beta * (r - re)} 
472:  
473:         energy = 0.D0 
474:         if(gtest) grad(:) = 0.D0 
475:  
476:         rijmod = sqrt(dot_product(rij, rij)) 
477:         if(rijmod == 0.D0) then 
478:             rij_hat(:) = 0.D0 
479:         else 
480:             rij_hat(:) = rij / rijmod 
481:         end if 
482:  
483:         exp_term = exp(-beta * (rijmod - requil)) 
484:         one_m_exp = 1.d0 - exp_term 
485:  
486:         energy = eps_0 * (one_m_exp * one_m_exp - 1.d0)  ! = eps * {[1 - exp{-beta * (r - re)}]^2 - 1} 
487:  
488:         if(gtest) then 
489:             grad(1: 3) = 2.d0 * eps_0 * beta * one_m_exp * exp_term * rij_hat(:) 
490:             grad(4: 6) = -grad(1: 3) 
491:             grad(7: 12) = 0.d0 
492:         end if 
493:  
494:     end subroutine sandbox_morse 
495: ! --- end sandbox potentials 
496:  
497:  
498: ! --- begin sandbox utilities 
499:     function lower_case(word) 
500:         character(len=*), intent(in) :: word 
501:         character(len=len(word)) :: lower_case 
502:  
503:         character(len=26), parameter :: uc = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' 
504:         character(len=26), parameter :: lc = 'abcdefghijklmnopqrstuvwxyz' 
505:  
506:         integer :: i, ind 
507:  
508:         do i = 1, len(word) 
509:             ind = index(uc, word(i: i)) 
510:             if (ind /= 0) then 
511:                 lower_case(i: i) = lc(ind: ind) 
512:             else 
513:                 lower_case(i: i) = word(i: i) 
514:             end if 
515:         end do 
516:  
517:     end function lower_case 
518:  
519:     function first_char(u) result(c) 
520:         implicit none 
521:         integer, intent(in) :: u 
522:         character :: c 
523:         integer :: stat 
524:  
525:         do 
526:             read(u, '(a)', advance='no', iostat=stat) c 
527:             if(stat /= 0) then 
528:                 c = '' 
529:                 backspace(u) 
530:                 return 
531:             elseif(c == ' ') then 
532:                 cycle 
533:             else 
534:                 backspace(u) 
535:                 return 
536:             end if 
537:         end do 
538:     end function first_char 
539:  
540:     subroutine print_matrix(m) 
541:         implicit none 
542:         double precision, intent(in) :: m(3, 3) 
543:  
544:         integer :: i 
545:          
546:         do i = 1, 3 
547:             print *, m(i, 1), m(i, 2), m(i, 3) 
548:         end do 
549:  
550:     end subroutine print_matrix 
551:  
552:     function max_distance(mols) result(d) 
553:         ! compute the maximum distance of a molecule from the origin 
554:  
555:         use vec3, only: vec_len 
556:  
557:         implicit none 
558:         type(sandbox_molecule), intent(in) :: mols(:) 
559:         double precision :: d 
560:  
561:         double precision :: this_d 
562:         integer :: i 
563:  
564:         d = 0.D0 
565:  
566:         do i = 1, size(mols) 
567:             this_d = vec_len(mols(i)%r) 
568:             if(this_d > d) d = this_d 
569:         end do 
570:  
571:     end function max_distance 
572: ! --- end sandbox utilities 
573:  
574: end module sandbox_module 
575:  
576: subroutine sandbox_input(out_unit) 
577:     use commons, only: natoms 
578:     use key, only: uniaxarray 
579:     use sandbox_module 
580:     implicit none     
581:  
582:     integer, intent(in) :: out_unit ! probably GMIN_out 
583:  
584:     integer :: sand_unit, stat 
585:     type(sandbox_molecule), pointer :: mol 
586:     type(sandbox_site), pointer :: site 
587:     integer :: arity, num_these_mols, num_sites, num_classes 
588:     character(len=20) :: nary_check, label(3) 
589:     double precision :: number_fraction 
590:     integer :: getunit, mol_index, mol_type_index, site_index 
591:     character(len=20) :: buffer 
592:  
593:     integer :: classi, classj, tmp 
594:     character(len=20) :: potential_name 
595:     type(interaction), pointer :: this_int 
596:     character(len=20), pointer :: func 
597:     double precision, pointer :: params(:) 
598:  
599:     character :: c, atom_symbol 
600:     integer :: i, j 
601:  
602:     if(.not.allocated(uniaxarray)) allocate(uniaxarray(natoms)) 
603:     if(.not.allocated(vt)) allocate(vt(3*natoms)) 
604:  
605:     x_size = 3 * natoms 
606:     num_mols = natoms / 2 
607:     allocate(molecules(num_mols)) 
608:     WRITE(*,*) 'SANDBOX - ', size(molecules) 
609:     num_atoms = 0 
610:     num_classes = 0 
611:  
612:     sand_unit = getunit() 
613:     open(unit=sand_unit, file='rbdata', status='old') 
614:  
615:     ! determine if first line is an integer num_sites or a string 'n-ary' 
616:     read(sand_unit, *) buffer 
617:     read(buffer, *, iostat=stat) num_sites 
618:     if(stat == 0) then 
619:         ! it was an integer, so do single site 
620:         arity = 1 
621:         num_these_mols = num_mols 
622:     else 
623:         ! it was hopefully a string 
624:         read(buffer, *, iostat=stat) nary_check 
625:  
626:         ! sanity check: should have copied a string 'n-ary' 
627:         if(.not.(stat == 0 .and. nary_check == 'n-ary')) then 
628:             print *, 'woah, son: got to put in some proper input for rb!' 
629:             return 
630:         end if 
631:  
632:         ! it was a string 'n-ary', so do n-ary mixture 
633:         read(sand_unit, *) arity 
634:         write(out_unit, *) 'sandbox_input> arity of mixture:', arity 
635:     end if 
636:  
637:     mol_index = 1 
638:     do mol_type_index = 1, arity 
639:         read(sand_unit, *)  ! read blank line 
640:         if(arity == 1) then 
641:             write(out_unit, *) 'sandbox_input>', num_these_mols, 'rigid bodies' 
642:         elseif(arity > 1) then 
643:             atom_symbol = '' 
644:             read(sand_unit, *) number_fraction, atom_symbol  ! number fraction of this molecule type 
645:             read(sand_unit, *) num_sites    ! number of sites in this molecule 
646:             if(atom_symbol == '') atom_symbol = 'O' 
647:  
648:             num_these_mols = nint(num_mols * number_fraction) 
649:             write(out_unit, *) 'adding', num_these_mols, 'molecules of type', mol_type_index 
650:         end if 
651:  
652:         if(num_these_mols == 0) then    ! the number fraction is too small to give even one molecule 
653:             read(sand_unit, *)  ! just skip the line 
654:         else 
655:             num_atoms=num_atoms+num_these_mols*num_sites 
656:             allocate(molecules(mol_index)%sites(num_sites)) ! prepare space for the sites in the molecule 
657:  
658:             site_index = 1 
659:             do 
660:                 if(site_index > num_sites) exit 
661:  
662:                 ! skip blank lines and lines beginning with # 
663:                 c = first_char(sand_unit) 
664:                 if(c == '#' .or. c == '') then 
665:                     read(sand_unit, *, iostat=stat) 
666:                     if(stat /= 0) exit 
667:                     cycle 
668:                 end if 
669:  
670:                 ! if the line wasn't blank or comment, start assigning values 
671:                 molecules(mol_index)%atom_symbol = atom_symbol 
672:                 site => molecules(mol_index)%sites(site_index) 
673:                 read(sand_unit, *) label(1), site%class, site%name, site%r_molframe(:), label(3), site%p_molframe(:) 
674:  
675:                 ! update number of classes 
676:                 if(site%class > num_classes) num_classes = site%class 
677:  
678:                 ! increment the site index 
679:                 site_index = site_index + 1 
680:             end do 
681:             mol_index = mol_index + 1 
682:  
683:             ! copy the first molecule into the rest of the mols of this type 
684:             do mol_index = mol_index, mol_index + (num_these_mols - 2) 
685:                 allocate(molecules(mol_index)%sites(num_sites)) 
686:                 molecules(mol_index) = molecules(mol_index - 1) 
687:             end do 
688:         end if 
689:     end do 
690:  
691:     ! read in blank lines until reach 'matrix' 
692:     do 
693:         read(sand_unit, *, iostat=stat) buffer 
694:         if(stat /= 0) then 
695:             print *, "sandbox> where's the interation matrix?" 
696:             return 
697:         end if 
698:         if(buffer == 'matrix') exit 
699:     end do 
700:  
701:     ! assign the interaction matrix 
702:     allocate(classes(num_classes)) 
703:     classes%is_aniso = .false. 
704:     classes%has_pole = .false. 
705:     allocate(interactions(num_classes, num_classes)) 
706:     interactions%exists = .false. 
707:     do 
708:         ! skip blank lines and lines beginning with # 
709:         c = first_char(sand_unit) 
710:         if(c == '#' .or. c == '') then 
711:             read(sand_unit, *, iostat=stat) 
712:             if(stat /= 0) then 
713:                 exit 
714:             end if 
715:             cycle 
716:         end if 
717:  
718:         read(sand_unit, *, iostat=stat) classi, classj, potential_name 
719:         if(stat /= 0) exit 
720:  
721:         if(classi > num_classes .or. classj > num_classes) then 
722:             print *, 'sandbox_input> youve specified interactions for classes with no sites' 
723:             exit 
724:         end if 
725:  
726:         backspace(sand_unit)    ! get ready to read same line again 
727:  
728:         ! keep interaction information in upper-right matrix 
729:         if(classi > classj) then 
730:             tmp = classi 
731:             classi = classj 
732:             classj = tmp 
733:         end if 
734:  
735:         ! assign existence 
736:         interactions(classi, classj)%exists = .true. 
737:         call locase(potential_name) 
738:  
739:         ! read in input params and do some internal processing if necessary 
740:         if(potential_name == 'lj') then 
741:             ! lj sigma_0 eps_0 rep att 
742:             interactions(classi, classj)%func_id = 1 
743:             allocate(interactions(classi, classj)%params(4)) 
744:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 4) 
745:         else if(potential_name == 'coulomb') then 
746:             ! coulomb sigma_0 k_c q q -> sigma_0 kcqq 
747:             interactions(classi, classj)%func_id = 2 
748:             allocate(interactions(classi, classj)%params(4)) 
749:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 4) 
750:             interactions(classi, classj)%params(2) = interactions(classi, classj)%params(2) * & 
751:                 & interactions(classi, classj)%params(3) * interactions(classi, classj)%params(4) 
752:         else if(potential_name == 'dipole') then 
753:             ! dipole sigma_0 mu -> sigma_0 mu^2 
754:             interactions(classi, classj)%func_id = 3 
755:             allocate(interactions(classi, classj)%params(2)) 
756:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 2) 
757:             interactions(classi, classj)%params(2) = interactions(classi, classj)%params(2) ** 2 
758:             classes(classi)%is_aniso = .true. 
759:             classes(classj)%is_aniso = .true. 
760:             classes(classi)%has_pole = .true. 
761:             classes(classj)%has_pole = .true. 
762:         else if(potential_name == 'chiro') then 
763:             ! chiro sigma_0 mu gamma -> sigma_0 mu^2 cos_gamma sin_gamma 
764:             interactions(classi, classj)%func_id = 4 
765:             allocate(interactions(classi, classj)%params(4)) 
766:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 3) 
767:             interactions(classi, classj)%params(2) = interactions(classi, classj)%params(2) ** 2 
768:             interactions(classi, classj)%params(4) = sin(interactions(classi, classj)%params(3)) 
769:             interactions(classi, classj)%params(3) = cos(interactions(classi, classj)%params(3)) 
770:             classes(classi)%is_aniso = .true. 
771:             classes(classj)%is_aniso = .true. 
772:             classes(classi)%has_pole = .true. 
773:             classes(classj)%has_pole = .true. 
774:         else if(potential_name == 'morse') then 
775:             ! morse eps_0 beta r_equilirbium -> eps beta re 
776:             interactions(classi, classj)%func_id = 5 
777:             allocate(interactions(classi, classj)%params(5)) 
778:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 3) 
779:         else if(potential_name == 'coulomb_ddiel') then 
780:             ! coulomb_ddiel sigma_0 k_c q q -> sigma_0 kcqq 
781:             interactions(classi, classj)%func_id = 6 
782:             allocate(interactions(classi, classj)%params(4)) 
783:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 4) 
784:             interactions(classi, classj)%params(2) = interactions(classi, classj)%params(2) * & 
785:                 & interactions(classi, classj)%params(3) * interactions(classi, classj)%params(4) 
786:         else 
787:             print *, 'sandbox_input> entered unknown interaction name!' 
788:             exit 
789:         end if 
790:     end do 
791:  
792:     ! initialize molecules 
793:     do mol_index = 1, size(molecules) 
794:         call initialize_sandbox_molecule(molecules(mol_index), mol_index) 
795:     end do 
796:  
797: end subroutine sandbox_input 
798:  
799: subroutine sandbox(x, grad, energy, gtest, stest) 
800:     use key, only: twod, frozen 
801:     use sandbox_module 
802:     implicit none 
803:     double precision, intent(inout) :: x(x_size) 
804:     logical, intent(in) :: gtest, stest 
805:     double precision, intent(out) :: energy, grad(x_size) 
806:  
807:     type(sandbox_molecule), pointer :: moli, molj 
808:     type(sandbox_site), pointer :: sitei, sitej 
809:     integer moli_index, molj_index, sitei_index, sitej_index 
810:     double precision :: energy_contrib, grad_contrib(12) 
811:  
812:     energy = 0.D0 
813:     vt(:) = 0.D0 
814:     if(gtest) grad(:) = 0.D0 
815:  
816:     ! update the values for all the molecules 
817:     do moli_index = 1, num_mols 
818:         moli => molecules(moli_index) 
819:         call update_sandbox_molecule(gtest, moli, x(moli%r_index: moli%r_index + 2), & 
820:             & x(moli%p_index: moli%p_index + 2)) 
821:     end do 
822:  
823:     ! outer loop over molecules 
824:     do moli_index = 1, num_mols - 1 
825:         moli => molecules(moli_index) 
826:  
827:         ! inner loop over molecules 
828:         do molj_index = moli_index + 1, num_mols 
829:             molj => molecules(molj_index) 
830:  
831:             ! loop over sites in outer molecule 
832:             do sitei_index = 1, size(moli%sites) 
833:                 sitei => moli%sites(sitei_index) 
834:  
835:                 ! loop over sites in inner molecule 
836:                 do sitej_index = 1, size(molj%sites) 
837:                     sitej => molj%sites(sitej_index) 
838:  
839:                     energy_contrib = 0.D0 
840:                     grad_contrib(:) = 0.D0 
841:  
842:                     ! compute the energy and gradient for this pair of sites 
843:                     call pairwise_sandbox(gtest, sitei, sitej, energy_contrib, grad_contrib) 
844:  
845:                     ! add the energy to total and pairwise 
846:                     energy = energy + energy_contrib 
847:                     vt(moli_index) = vt(moli_index) + energy_contrib 
848:                     vt(molj_index) = vt(molj_index) + energy_contrib 
849:  
850:                     if(gtest) then  ! add gradient contributions 
851:                         grad(moli%r_index: moli%r_index + 2) = grad(moli%r_index: moli%r_index + 2) + grad_contrib(1: 3) 
852:                         grad(molj%r_index: molj%r_index + 2) = grad(molj%r_index: molj%r_index + 2) + grad_contrib(4: 6) 
853:                         grad(moli%p_index: moli%p_index + 2) = grad(moli%p_index: moli%p_index + 2) + grad_contrib(7: 9) 
854:                         grad(molj%p_index: molj%p_index + 2) = grad(molj%p_index: molj%p_index + 2) + grad_contrib(10: 12) 
855:                     end if 
856:  
857:                 end do 
858:             end do 
859:         end do 
860:     end do 
861:  
862:     ! 2D: set all z gradients to zero 
863:     if(twod) then 
864:         do moli_index = 1, num_mols 
865:             grad(molecules(moli_index)%r_index + 2) = 0.D0 
866:         end do 
867:     end if 
868:     ! freeze keyword: set gradients to zero 
869:     if(gtest .and. any(frozen)) then 
870: !        STOP 
871:         do moli_index = 1, num_mols 
872:             moli => molecules(moli_index) 
873:             if(frozen(moli_index)) grad(moli%r_index: moli%r_index + 2) = 0.d0 
874:             if(frozen(moli_index + num_mols)) grad(moli%p_index: moli%p_index + 2) = 0.d0 
875:         end do 
876:     end if 
877:  
878: end subroutine sandbox 
879:  
880: subroutine sandbox_output 
881:     use commons, only: natoms  
882:     use sandbox_module 
883:     implicit none 
884:  
885:     integer :: sandout_unit, coords_unit 
886:     integer :: min_num, mol_num, site_num, atom_num, atom_index 
887:     type(sandbox_molecule), pointer :: mol 
888:     type(sandbox_site), pointer :: site 
889:     character(len=25) :: sandout_name, coords_name, min_name, node 
890:     double precision :: rotmat(3, 3), dummy(3, 3) 
891:  
892:     integer :: getunit 
893:  
894:       integer,allocatable :: FF(:),INTEFF(:) ! NSAVE 
895:       INTEGER, ALLOCATABLE :: NPCALL_QMIN(:) ! NSAVE 
896:       DOUBLE PRECISION, ALLOCATABLE :: QMIN(:), INTEQMIN(:) ! NSAVE 
897:       DOUBLE PRECISION, ALLOCATABLE :: QMINP(:,:), INTEQMINP(:,:) 
898:       INTEGER, ALLOCATABLE :: QMINT(:,:), QMINNATOMS(:) 
899:  
900:  
901:     sandout_unit = getunit() 
902:  
903:     ! open sandout file for writing 
904:     if(mpit) then 
905:         write(node, *) mynode + 1 
906:         sandout_name = 'sandout.' // trim(adjustl(node)) // '.xyz' 
907:     else 
908:         sandout_name = 'sandout.xyz' 
909:     end if 
910:     open(unit=sandout_unit, file=sandout_name, status='unknown') 
911:  
912:     ! loop over saved minima 
913:     do min_num = 1, nsave 
914:         ! put number of atoms and comment line. for now this we'll do just one atom per molecule 
915:         !write(sandout_unit, '(I0)') num_mols 
916:         write(sandout_unit, '(I0)') num_atoms 
917:         write(sandout_unit, '(A,1x,I0,A,1x,F20.10,1x,A,1x,I0)') 'energy of minimum', min_num, '=', qmin(min_num), 'first found at step', ff(min_num) 
918:  
919:         ! loop over molecules 
920:         do mol_num = 1, size(molecules) 
921:             mol => molecules(mol_num) 
922:  
923:             call rmdrvt(qminp(min_num, mol%p_index: mol%p_index + 2), rotmat, dummy, dummy, dummy, .false.) 
924:  
925:             call update_sandbox_molecule(.false., mol, qminp(min_num, mol%r_index: mol%r_index + 2),qminp(min_num,mol%p_index: mol%p_index + 2)) 
926:  
927:             !write(sandout_unit, '(A6,3F10.5,A12,1X,3F10.5)') mol%atom_symbol, qminp(min_num, mol%r_index: mol%r_index + 2), & 
928:              !   & 'atom_vector', matmul(rotmat, z_hat) 
929:             do atom_num=1,size(molecules(mol_num)%sites) 
930:                write(sandout_unit, '(A6,3F10.5,A12,1X,3F10.5)')molecules(mol_num)%sites(atom_num)%name,molecules(mol_num)%sites(atom_num)%r 
931:             enddo 
932:         end do 
933:  
934:         ! output coords. files 
935:         coords_unit = getunit() 
936:         write(min_name, '(i3)') min_num 
937:  
938:         if(mpit) then 
939:             coords_name = 'coords.' // trim(adjustl(min_name)) // '.' // trim(adjustl(node)) 
940:         else 
941:             coords_name = 'coords.' // trim(adjustl(min_name)) 
942:         end if 
943:  
944:         open(coords_unit, file=coords_name, status='unknown') 
945:         do atom_num = 1, natoms 
946:             atom_index = 3 * (atom_num - 1) + 1 
947:             write(coords_unit, *) qminp(min_num, atom_index), qminp(min_num, atom_index + 1), qminp(min_num, atom_index + 2) 
948:         end do 
949:         close(coords_unit) 
950:  
951:     end do 
952:  
953:     close(sandout_unit) 
954:              
955: end subroutine sandbox_output 
956:  
957: subroutine sandbox_takestep(np) 
958:     use key, only: myunit,twod, frozen 
959:     use sandbox_module 
960:     use vec3, only: vec_len, vec_random 
961:     use rotations, only: rot_takestep_aa 
962:  
963:     implicit none 
964:     integer, intent(in) :: np 
965:  
966:     double precision, target :: x(size(coords(:))) 
967:     double precision, pointer :: r(:), p(:) 
968:     double precision :: step_size, min_vt, dice_roll, twod_step(3) 
969:     type(sandbox_molecule), pointer :: moli, molj 
970:     integer :: moli_index, molj_index 
971:     logical :: stop_pulling ! for angular moves 
972:  
973:     min_vt = minval(vt) 
974:     x(:)=coords(:) 
975:     ! update ellipsoids before taking any steps 
976:     do moli_index = 1, num_mols 
977:         moli => molecules(moli_index) 
978:         r => x(moli%r_index: moli%r_index + 2) 
979:         p => x(moli%p_index: moli%p_index + 2) 
980:  
981:         if(.not. percolatet) then 
982:             if(vec_len(moli%r) > radius) then 
983:                 write(myunit, *) 'sandbox_takestep> initial coord outside container, brining in' 
984:                 r(:) = r(:) - sqrt(radius) * nint(r(:) / sqrt(radius)) 
985:             end if 
986:         end if 
987:  
988:         call update_sandbox_molecule(.false., moli, r, p) 
989:     end do 
990:  
991:     do moli_index = 1, num_mols 
992:         moli => molecules(moli_index) 
993:         r => x(moli%r_index: moli%r_index + 2) 
994:         p => x(moli%p_index: moli%p_index + 2) 
995:  
996:         if(twod) then 
997:             ! translational step 
998:             do 
999:                 call random_number(twod_step(1)) 
1000:                 call random_number(twod_step(2)) 
1001:                 twod_step(3) = 0.D0 
1002:                 if(vec_len(twod_step) <= 1.0) exit 
1003:             end do 
1004:             r(:) = r(:) + (step(np) * twod_step(:)) 
1005:         else 
1006:          if(.not.(frozen(moli_index) .or. frozen(moli_index + num_mols))) then 
1007:             ! angular move: accept step with probability 1 - exp(-beta * delta_energy), where beta = 1/kT = astep 
1008:             if(astep(np) > 0.D0) then 
1009:                 call random_number(dice_roll) 
1010:                 if(dice_roll > exp(-astep(np) * (vt(moli_index) - min_vt))) then 
1011:                     r = max_distance(molecules) * vec_random() 
1012:  
1013:                     ! if using percolate, then pull the molecule in until it is within percolation distance of another 
1014:                     if(percolatet) then 
1015:                         stop_pulling = .false. 
1016:                         do 
1017:                             do molj_index = 1, num_mols 
1018:                                 if(moli_index == molj_index) cycle 
1019:                                 molj => molecules(molj_index) 
1020:  
1021:                                 if(vec_len(moli%r - molj%r) < perccut) then 
1022:                                     stop_pulling = .true. 
1023:                                     exit 
1024:                                 end if 
1025:                             end do 
1026:  
1027:                             if(stop_pulling) exit 
1028:                             r(:) = 0.95 * r(:) 
1029:                         end do 
1030:                     else ! if using radius, pull in until the molecule is inside the container 
1031:                         do 
1032:                             if(vec_len(r) < radius) exit 
1033:                             r(:) = 0.95 * r(:) 
1034:                         end do 
1035:                     end if 
1036:  
1037:                 end if 
1038:             end if 
1039:          end if  
1040:             ! translational move: uniformly distributed in a sphere of radius step(np) 
1041:          if(.not. frozen(moli_index)) then 
1042:             if(tmove(np)) then 
1043:                 call random_number(step_size) 
1044:                 step_size = sqrt(step(np) * step_size)  ! need to take square root to get uniform volume distribution 
1045:                 r(:) = r(:) + step_size * vec_random() 
1046:             end if 
1047:          end if 
1048:         end if ! 2d 
1049:  
1050:         ! orientational move 
1051:          if(.not. frozen(moli_index + num_mols)) then 
1052:             if(omove(np)) then 
1053:               call random_number(step_size) 
1054:               call rot_takestep_aa(p(:), step_size * ostep(np)) 
1055:             end if 
1056:          end if 
1057:     end do 
1058:  
1059:     ! copy local results back 
1060:     coords(:) = x(:) 
1061:  
1062: end subroutine sandbox_takestep 
1063:  
1064: subroutine sandboxsecder(OLDX,STEST) 
1065: use commons 
1066: use modhess 
1067: implicit none 
1068:  
1069: DOUBLE PRECISION  :: V(3*NATOMS),EGB,X(3*NATOMS),OLDX(3*NATOMS),VTEMP(2,3*NATOMS),ksi 
1070: INTEGER    :: i,j 
1071: LOGICAL    :: GTEST,STEST 
1072:  
1073: ksi=0.00001D0 
1074: X(:)=OLDX(:) 
1075:  
1076: IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS,3*NATOMS)) 
1077: HESS(:,:)=0.0D0 
1078: V(:)=0.0D0 
1079: VTEMP(:,:)=0.0D0 
1080:  
1081: DO i=1,3*NATOMS 
1082:  
1083:          X(i)=X(i)-ksi 
1084:  
1085:          call sandbox(X,V,EGB,.TRUE.,.FALSE.) 
1086:          VTEMP(1,:)=V(:) 
1087:  
1088:          X(i)=X(i)+2.0D0*ksi 
1089:  
1090:          call sandbox(X,V,EGB,.TRUE.,.FALSE.) 
1091:  
1092:          VTEMP(2,:)=V(:) 
1093:  
1094:                 DO j=i,3*NATOMS 
1095:                         HESS(i,j)=(VTEMP(2,j)-VTEMP(1,j))/(2.0D0*ksi) 
1096:                         HESS(j,i)=HESS(i,j) 
1097:                 END DO 
1098: END DO 
1099:  
1100: end subroutine sandboxsecder 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0