hdiff output

r26189/build.csh 2017-01-21 10:33:53.698250000 +0000 r26188/build.csh 2017-01-21 10:33:55.786250000 +0000
822:      endif822:      endif
823:      set command = "$command T"823:      set command = "$command T"
824: # QMTYPE variable passed to Makefile to include SCC-DFTB libraries 824: # QMTYPE variable passed to Makefile to include SCC-DFTB libraries 
825:      set QMTYPE = DFTB825:      set QMTYPE = DFTB
826:    endif826:    endif
827:    if ( $compiler == pgi ) then827:    if ( $compiler == pgi ) then
828:      if($CTYPE == C36) then #compile.csh adjustments828:      if($CTYPE == C36) then #compile.csh adjustments
829:               export MPI_INCLUDE=/usr/local/shared/ubuntu-12.04/x86_64/openmpi-1.4.5-pgi12/include829:               export MPI_INCLUDE=/usr/local/shared/ubuntu-12.04/x86_64/openmpi-1.4.5-pgi12/include
830:                  export MPI_LIB=/usr/local/shared/ubuntu-12.04/x86_64/openmpi-1.4.5-pgi12/lib        830:                  export MPI_LIB=/usr/local/shared/ubuntu-12.04/x86_64/openmpi-1.4.5-pgi12/lib        
831:              #set command = "install.com gnu xxlarge PGF95 keepo keepf" 831:              #set command = "install.com gnu xxlarge PGF95 keepo keepf" 
832:              #set command = "install.com gnu xxlarge PGF95 keepo keepf E MPICH"  #mpif90832:              set command = "install.com gnu xxlarge PGF95 keepo keepf E MPICH"  #mpif90
833:          set command = "compile.csh" 
834:      else 833:      else 
835:              set command = "$command PGF90"834:              set command = "$command PGF90"
836:              if ( $level == opt ) then835:              if ( $level == opt ) then
837:                set command = "$command OPT"836:                set command = "$command OPT"
838:              endif837:              endif
839:              if ( $architecture == amd ) then838:              if ( $architecture == amd ) then
840:                set command = "$command AMD"839:                set command = "$command AMD"
841:              endif840:              endif
842:          endif841:          endif
843:    else if ( $compiler == ifort ) then842:    else if ( $compiler == ifort ) then


r26189/compile.csh 2017-01-21 10:33:52.506250000 +0000 r26188/compile.csh 2017-01-21 10:33:54.410250000 +0000
  1: #!/bin/bash -l  1: #!/bin/bash -l
  2: # jbr36  2: # jbr36
  3: # Charmm36 interface: this file is called from the OPTIM build.csh. The options below  3: # Charmm36 interface: this file ist not called from the OPTIM build.csh. The options below
  4: # are hard corded there for former interfaces   4: # are hard corded there. Do not forget to export the variables and load the module before running 
  5: # ./build.csh coptim c36 pgi  5: # ./build.csh coptim c36 pgi
  6: #  6: #
  7: # Default (ifort)  7: # Default (ifort)
  8: # tested on 'clust' with 'ifort/64/11.1/038' module  8: # tested on 'clust' with 'ifort/64/11.1/038' module
  9: #./install.com gnu large ifort keepo keepf > build.log  9: #./install.com gnu large ifort keepo keepf > build.log
 10:  10: 
 11: #adjust 11: #adjust
 12: #module load mpi/openmpi/64/intel12/1.4.5  12: #module load mpi/openmpi/64/intel12/1.4.5 
 13: module load mpi/openmpi/64/pgi12/1.4.5 13: module load mpi/openmpi/64/pgi12/1.4.5
 14: # from module display <modulename> 14: # from module display <modulename>


r26189/energy.src 2017-01-21 10:33:52.786250000 +0000 r26188/energy.src 2017-01-21 10:33:54.734250000 +0000
3645:      &               NATOM,NTRANS,IMTRNS,HEAP(BIMAG(IMATPT)),3645:      &               NATOM,NTRANS,IMTRNS,HEAP(BIMAG(IMATPT)),
3646:      &               HEAP(BIMAG(IMATTR)),NOROT,NATIM,IMINV)3646:      &               HEAP(BIMAG(IMATTR)),NOROT,NATIM,IMINV)
3647:           ENDIF3647:           ENDIF
3648:         ENDIF3648:         ENDIF
3649:       ENDIF                  !##PBOUND3649:       ENDIF                  !##PBOUND
3650: ##ENDIF  (not_noimag)3650: ##ENDIF  (not_noimag)
3651: ##ENDIF  (lrst)3651: ##ENDIF  (lrst)
3652: CCC3652: CCC
3653: ##ENDIF (ldm)3653: ##ENDIF (ldm)
3654: ##IF ENSEMBLE3654: ##IF ENSEMBLE
3655: CC jbr36 adjusted for glutamate mutase EVB from Bristol3655:       IF (QENSEXP) THEN
3656:       IF (QENSEXP.OR.QEVB.OR.QRPMD) THEN3656:          CALL ENSEXPAVG(DX,DY,DZ,EPROP(EPOT),NATOM)
3657:          CALL ENSEXPAVG(X,Y,Z,DX,DY,DZ,EPROP(EPOT),NATOM) 
3658:       ENDIF3657:       ENDIF
3659: ##ENDIF3658: ##ENDIF
3660: C---------------------------------------------------------------3659: C---------------------------------------------------------------
3661: ##IF ESTATS  (estats)3660: ##IF ESTATS  (estats)
3662:       IF (LESTAT) THEN3661:       IF (LESTAT) THEN
3663:        IENCAL = IENCAL + 13662:        IENCAL = IENCAL + 1
3664:        IF (IENCAL.LE.MXENCL) THEN !if not past end3663:        IF (IENCAL.LE.MXENCL) THEN !if not past end
3665:          IF (IENCAL.GT.ESTRCL) THEN !if after start3664:          IF (IENCAL.GT.ESTRCL) THEN !if after start
3666:            IF (MOD((IENCAL-ESTRCL),ESTMOD).EQ.0) THEN3665:            IF (MOD((IENCAL-ESTRCL),ESTMOD).EQ.0) THEN
3667:                CALL ESTACAL3666:                CALL ESTACAL


r26189/ensemble.fcm 2017-01-21 10:33:53.098250000 +0000 r26188/ensemble.fcm 2017-01-21 10:33:55.026250000 +0000
 18: C FLIPFLOP: tells ensemble which replicas we are trying to 18: C FLIPFLOP: tells ensemble which replicas we are trying to
 19: C           swap each time 19: C           swap each time
 20: C ENSATT: statistics of attempted swaps 20: C ENSATT: statistics of attempted swaps
 21: C ENSSUC: statistics of successful swaps 21: C ENSSUC: statistics of successful swaps
 22: C JREX: logical to tell whether replicas have been set up for 22: C JREX: logical to tell whether replicas have been set up for
 23: C       replica exchange. 23: C       replica exchange.
 24: C JRSWAP: is replica exchange actually turned on at the moment? 24: C JRSWAP: is replica exchange actually turned on at the moment?
 25: C ENSTEM: list of replica temperatures sorted in ascending order 25: C ENSTEM: list of replica temperatures sorted in ascending order
 26: C ENSDB: lookup table of beta_i-beta_j factors 26: C ENSDB: lookup table of beta_i-beta_j factors
 27: C ENSMYT: temperature of present replica  27: C ENSMYT: temperature of present replica 
 28: c  
 29: C 28: C
 30: C the following are for code which is still in testing... 29: C the following are for code which is still in testing...
 31: C ENSDX,ENSDY,ENSDZ,ENSH,ENSEXPU,QPARPRM,QENSEXP,SWAPC,JENSC, 30: C ENSDX,ENSDY,ENSDZ,ENSH,ENSEXPU,QPARPRM,QENSEXP,SWAPC,JENSC,
 32: C ENSEBETA, EXPOFF 31: C ENSEBETA, EXPOFF
 33: C 32: C
 34: c COUPRM: table where each row holds parameters for EVB gaussian 
 35: c   coupling matrix elements (constants & 1d/2d gaussians) 
 36: c 
 37: c  col 1: amplitude, A  (constant, 1d, or 2d gaussian) 
 38: c  col 2: center, R01   (1d or 2d gaussian) 
 39: c  col 3: width,  C1    (1d or 2d gaussian) 
 40: c  col 4: theta         (2d gaussian only) 
 41: c  col 5: center, R02   (2d gaussian only) 
 42: c  col 6: width,  C2    (2d gaussian only) 
 43: c  col 7: a             (2d gaussian only)   
 44: c  col 8: b             (2d gaussian only) 
 45: c  col 9: c             (2d gaussian only) 
 46: c 
 47: c 
 48: c 
 49: ##IF ENSEMBLE 33: ##IF ENSEMBLE
 50:       INTEGER WHOIAM,NENSEM,MAXENS,T2REPO,REP2TO,ENSBFL,ENSBUF,MAXSWP 34:       INTEGER WHOIAM,NENSEM,MAXENS,T2REPO,REP2TO,ENSBFL,ENSBUF,MAXSWP
 51: c-----PVG MPISETUP for RPMD AND EVB 
 52:       INTEGER NEVBS   
 53:       INTEGER NRPBDS  
 54:       INTEGER STATEID 
 55:       INTEGER BEADID  
 56:       INTEGER REPIDS(512,512) ! stores the beadID and stateID for each replica  
 57:       INTEGER interBeadGroup, interBeadGroupComm 
 58:       INTEGER interStateGroup, interStateGroupComm 
 59:       INTEGER NRPMDSEL 
 60: c-----PVG RPMD  
 61:       REAL*8  FCTEMP 
 62:       LOGICAL QRPMD, QEVBRP, QRPRXNCRD 
 63:       INTEGER FULLXYZFRQ, FULLXYZUNIT, CENTROIDXYZFRQ, CENTROIDXYZUNIT 
 64:       INTEGER RPRXNCRDUNIT 
 65:       INTEGER RPRXNSEL(4),RPRXNFRQ 
 66:       LOGICAL RPDISCRD, RPDIFFCRD 
 67: c-----END PVG RPMD  
 68:       INTEGER SWBUF 35:       INTEGER SWBUF
 69:       INTEGER ENSFRQ,ENSNSW 36:       INTEGER ENSFRQ,ENSNSW
 70:       INTEGER ENSDX,ENSDY,ENSDZ,ENSH,ENSEXPU 37:       INTEGER ENSDX,ENSDY,ENSDZ,ENSH,ENSEXPU
 71:       PARAMETER (MAXENS=512) 38:       PARAMETER (MAXENS=512)
 72:       PARAMETER (ENSBFL=120) 39:       PARAMETER (ENSBFL=120)
 73:       PARAMETER (MAXSWP=(MAXENS+1)*MAXENS/2) 40:       PARAMETER (MAXSWP=(MAXENS+1)*MAXENS/2)
 74:       INTEGER T2REP(MAXENS),REP2T(MAXENS),ENSATT(MAXSWP),ENSSUC(MAXSWP) 41:       INTEGER T2REP(MAXENS),REP2T(MAXENS),ENSATT(MAXSWP),ENSSUC(MAXSWP)
 75:       INTEGER ENSISW(MAXSWP),ENSJSW(MAXSWP) 42:       INTEGER ENSISW(MAXSWP),ENSJSW(MAXSWP)
 76:       INTEGER REPSWP(MAXENS) 43:       INTEGER REPSWP(MAXENS)
 77:       INTEGER ESEL(MAXENS*MAXENS,2),ESEL2(MAXENS*MAXENS,2) 44:       LOGICAL JREX,JRSWAP,JENSC,QPARPRM,QENSEXP,SWAPC,QEXPBW
 78:       INTEGER RPMDSEL(MAXENS) 
 79:       INTEGER EVBIIDX(MAXENS),EVBJIDX(MAXENS),NREAD 
 80:       LOGICAL JREX,JRSWAP,JENSC,QPARPRM,QENSEXP,SWAPC,QEXPBW,QEVB 
 81:       LOGICAL ONEDCRD,TWODCRD,DIFFCRD 
 82:       CHARACTER*4 COUTYPE(MAXENS) 
 83:  
 84:       REAL*8 ENSTEM(MAXENS),ENSDB(MAXENS),ENSMYT, 45:       REAL*8 ENSTEM(MAXENS),ENSDB(MAXENS),ENSMYT,
 85:      1       ENSEBETA,EXPOFF(MAXENS),COUPRM(MAXENS*MAXENS,9) 46:      1       ENSEBETA,EXPOFF(MAXENS)
 86:  
 87:       COMMON /ENSCOM/ENSTEM,ENSMYT,ENSDB,ENSEBETA,EXPOFF, 47:       COMMON /ENSCOM/ENSTEM,ENSMYT,ENSDB,ENSEBETA,EXPOFF,
 88:      1       WHOIAM,NENSEM,T2REPO,REP2TO,ENSBUF, 48:      1       WHOIAM,NENSEM,T2REPO,REP2TO,ENSBUF,
 89:      2       ENSFRQ,T2REP,REP2T,REPSWP,ENSISW,ENSJSW,ENSNSW, 49:      2       ENSFRQ,T2REP,REP2T,REPSWP,ENSISW,ENSJSW,ENSNSW,
 90:      3       ENSATT,ENSSUC,ENSDX,ENSDY,ENSDZ,ENSH,ENSEXPU, 50:      3       ENSATT,ENSSUC,
 91:      5       JREX,JRSWAP,JENSC,QPARPRM,QENSEXP,SWAPC,QEXPBW, 51:      4       ENSDX,ENSDY,ENSDZ,ENSH,ENSEXPU,
 92:      6       QEVB,COUPRM,ESEL,ESEL2,ONEDCRD,TWODCRD,EVBIIDX, 52:      5       JREX,JRSWAP,JENSC,QPARPRM,QENSEXP,SWAPC,QEXPBW
 93:      7       EVBJIDX,COUTYPE,NREAD,NRPBDS,NEVBS,STATEID,BEADID, 
 94:      8       REPIDS,interBeadGroup,interBeadGroupComm, 
 95:      9       interStateGroup,interStateGroupComm, 
 96:      &       FCTEMP,QRPMD,QEVBRP,RPMDSEL,NRPMDSEL,FULLXYZFRQ, 
 97:      &       FULLXYZUNIT,CENTROIDXYZFRQ,CENTROIDXYZUNIT, 
 98:      &       RPDISCRD, RPDIFFCRD, RPRXNCRDUNIT, QRPRXNCRD,  
 99:      &       RPRXNSEL,RPRXNFRQ 
100: ##ENDIF 53: ##ENDIF


r26189/ensemble.src 2017-01-21 10:33:53.438250000 +0000 r26188/ensemble.src 2017-01-21 10:33:55.398250000 +0000
  1: CHARMM Element source/machdep/hqbm.src $Revision: 1.3 $  1: CHARMM Element source/machdep/hqbm.src $Revision: 1.5 $
  2: C  2: C
  3: C ensemble code  3: C ensemble code
  4: C documentation in ensemble.doc  4: C documentation in ensemble.doc
  5: C  5: C
  6:  
  7:  
  8:  
  9: ##IF ENSEMBLE  6: ##IF ENSEMBLE
 10:       SUBROUTINE ENSINI  7:       SUBROUTINE ENSINI
 11: C----------------------------------------------------------------------  8: C----------------------------------------------------------------------
 12: C Initialize ensemble when charmm is started  9: C Initialize ensemble when charmm is started
 13: C---------------------------------------------------------------------- 10: C----------------------------------------------------------------------
 14: ##INCLUDE '~/charmm_fcm/impnon.fcm' 11: ##INCLUDE '~/charmm_fcm/impnon.fcm'
 15: ##INCLUDE '~/charmm_fcm/stream.fcm' 12: ##INCLUDE '~/charmm_fcm/stream.fcm'
 16: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 13: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
 17:       INCLUDE 'mpif.h' 14:       INCLUDE 'mpif.h'
 18:       INTEGER L 15:       INTEGER L
 22:       CALL MPI8INIT(IERROR) 19:       CALL MPI8INIT(IERROR)
 23:       CALL MPI8COMM_RANK(MPI_COMM_WORLD, WHOIAM, IERROR) 20:       CALL MPI8COMM_RANK(MPI_COMM_WORLD, WHOIAM, IERROR)
 24:       CALL MPI8COMM_SIZE(MPI_COMM_WORLD, NENSEM, IERROR) 21:       CALL MPI8COMM_SIZE(MPI_COMM_WORLD, NENSEM, IERROR)
 25:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR) 22:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR)
 26: ##ELSE 23: ##ELSE
 27:       CALL MPI_INIT(IERROR) 24:       CALL MPI_INIT(IERROR)
 28:       CALL MPI_COMM_RANK(MPI_COMM_WORLD, WHOIAM, IERROR) 25:       CALL MPI_COMM_RANK(MPI_COMM_WORLD, WHOIAM, IERROR)
 29:       CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NENSEM, IERROR) 26:       CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NENSEM, IERROR)
 30:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR) 27:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR)
 31: ##ENDIF 28: ##ENDIF
  29:       WRITE(TBUF, '(A,I5,A)')
  30:      1     ' ENSEMBLE>    REPLICA NODE ', WHOIAM, ' STARTED!'
  31:       CALL ENSPRN(OUTU,TBUF,LEN(TBUF))
 32: C     stop non-root processes from printing output 32: C     stop non-root processes from printing output
 33:       IF (WHOIAM.NE.0) THEN 33:       IF (WHOIAM.NE.0) THEN
 34:            PRNLEV = -1 34:            PRNLEV = -1
 35:            WRNLEV = -5 35:            WRNLEV = -5
 36:            IOLEV = -1 36:            IOLEV = -1
 37:       ENDIF 37:       ENDIF
 38: ce----PVG ERROR CHECKING 
 39: ce      ENDIF 
 40: ce      IF (WHOIAM.EQ.5) THEN 
 41: ce           PRNLEV = 5  
 42: ce      ENDIF 
 43: ce---PVG ERROR CHECKING 
 44:       WRITE(TBUF, '(A,I5,A)') 
 45:      1     ' ENSEMBLE>    REPLICA NODE ', WHOIAM, ' STARTED!' 
 46:       CALL ENSPRN(OUTU,TBUF,LEN(TBUF)) 
 47:       IF (NENSEM.GT.MAXENS) THEN 38:       IF (NENSEM.GT.MAXENS) THEN
 48:             IF (IOLEV.GT.0) WRITE (OUTU,'(A)') 39:             IF (IOLEV.GT.0) WRITE (OUTU,'(A)')
 49:      1            ' ENSINI> TOO MANY REPLICAS!', 40:      1            ' ENSINI> TOO MANY REPLICAS!',
 50:      2            ' ENSINI> INCREASE MAXENS IN ensemble.fcm' 41:      2            ' ENSINI> INCREASE MAXENS IN ensemble.fcm'
 51:             CALL WRNDIE(0, '<ENSINI>', 'UNRECOGNIZED SUBCOMMAND') 42:             CALL WRNDIE(0, '<ENSINI>', 'TOO MANY REPLICAS REQUESTED')
 52:       ENDIF     43:       ENDIF    
 53:       CALL SETMSI('WHOIAM',WHOIAM) 44:       CALL SETMSI('WHOIAM',WHOIAM)
 54:       CALL SETMSI('NENSEM',NENSEM) 45:       CALL SETMSI('NENSEM',NENSEM)
 55:       JREX = .FALSE. 46:       JREX = .FALSE.
 56:       JRSWAP = .FALSE. 47:       JRSWAP = .FALSE.
 57:       QENSEXP = .FALSE. 48:       QENSEXP = .FALSE.
 58:       QEVB = .FALSE. 
 59:       QRPMD = .FALSE. 
 60:       QEVBRP = .FALSE. 
 61:       ONEDCRD = .FALSE. 
 62:       TWODCRD = .FALSE. 
 63:       NREAD = 0 
 64: C     this is ignored: 49: C     this is ignored:
 65:       SWAPC = .TRUE. 50:       SWAPC = .TRUE.
 66:       T2REPO = -1 51:       T2REPO = -1
 67:       REP2TO = -1 52:       REP2TO = -1
 68:       ENSEXPU = -1 53:       ENSEXPU = -1
 69: C      ENSNT = 0 54: C      ENSNT = 0
 70:       IF (IOLEV.GT.0) THEN 55:       IF (IOLEV.GT.0) THEN
 71:           WRITE(OUTU,'(A)')  56:           WRITE(OUTU,'(A)') 
 72:      1    ' ENSINI> ALL REPLICAS PRESENTLY USE ONE TEMPERATURE', 57:      1    ' ENSINI> ALL REPLICAS PRESENTLY USE ONE TEMPERATURE',
 73:      2    ' ENSINI> USE "ENSEmble EXCHange"', 58:      2    ' ENSINI> USE "ENSEmble EXCHange"',
 74:      3    ' ENSINI> TO SET UP REPLICA EXCHANGE' 59:      3    ' ENSINI> TO SET UP REPLICA EXCHANGE.'
  60: ##IF STRINGM     
  61:      4    ,' ENSINI> USE ENSEmble STRIng [...]',
  62:      5    ' ENSINI> TO USE THE STRING METHOD.'
  63: ##ENDIF
 75:       ENDIF 64:       ENDIF
 76:       RETURN 65:       RETURN
 77:       END 66:       END
 78: C---------------------------------------------------------------------- 67: C----------------------------------------------------------------------
 79:       SUBROUTINE ENSFIN 68:       SUBROUTINE ENSFIN
 80: C---------------------------------------------------------------------- 69: C----------------------------------------------------------------------
 81: C clean up on exit 70: C clean up on exit
 82: C---------------------------------------------------------------------- 71: C----------------------------------------------------------------------
 83: ##INCLUDE '~/charmm_fcm/impnon.fcm' 72: ##INCLUDE '~/charmm_fcm/impnon.fcm'
 84: ##INCLUDE '~/charmm_fcm/stream.fcm' 73: ##INCLUDE '~/charmm_fcm/stream.fcm'
152: C     parse command line and call subcommands, namely141: C     parse command line and call subcommands, namely
153: C           OPEN142: C           OPEN
154: C           CLOSE143: C           CLOSE
155: C           SEED144: C           SEED
156: C     etc...145: C     etc...
157: C----------------------------------------------------------------------146: C----------------------------------------------------------------------
158: ##INCLUDE '~/charmm_fcm/impnon.fcm'147: ##INCLUDE '~/charmm_fcm/impnon.fcm'
159: ##INCLUDE '~/charmm_fcm/ensemble.fcm'148: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
160: ##INCLUDE '~/charmm_fcm/exfunc.fcm'149: ##INCLUDE '~/charmm_fcm/exfunc.fcm'
161: ##INCLUDE '~/charmm_fcm/stream.fcm'150: ##INCLUDE '~/charmm_fcm/stream.fcm'
 151: 
 152: ##INCLUDE '~/charmm_fcm/dimens.fcm'
 153: ##INCLUDE '~/charmm_fcm/coord.fcm'
 154: 
 155: 
162:       CHARACTER*(*) COMLYN156:       CHARACTER*(*) COMLYN
163:       INTEGER COMLEN157:       INTEGER COMLEN
164:       IF (INDXA(COMLYN,COMLEN,'OPEN') .GT. 0) THEN158:       IF (INDXA(COMLYN,COMLEN,'OPEN') .GT. 0) THEN
165:             CALL ENSOPN(COMLYN,COMLEN)159:             CALL ENSOPN(COMLYN,COMLEN)
166:       ELSE IF (INDXA(COMLYN,COMLEN,'CLOSE') .GT. 0) THEN160:       ELSE IF (INDXA(COMLYN,COMLEN,'CLOSE') .GT. 0) THEN
167:             CALL ENSCLO(COMLYN,COMLEN)161:             CALL ENSCLO(COMLYN,COMLEN)
168:       ELSE IF (INDXA(COMLYN,COMLEN,'SEED') .GT. 0) THEN162:       ELSE IF (INDXA(COMLYN,COMLEN,'SEED') .GT. 0) THEN
169:             CALL ENSEED(COMLYN,COMLEN)163:             CALL ENSEED(COMLYN,COMLEN)
170:       ELSE IF (INDXA(COMLYN,COMLEN,'PRSEED') .GT. 0) THEN164:       ELSE IF (INDXA(COMLYN,COMLEN,'PRSEED') .GT. 0) THEN
171:             CALL ENSDPR(COMLYN,COMLEN)165:             CALL ENSDPR(COMLYN,COMLEN)
172:       ELSE IF (INDXA(COMLYN,COMLEN,'EXPAVG') .GT. 0) THEN166:       ELSE IF (INDXA(COMLYN,COMLEN,'EXPAVG') .GT. 0) THEN
173:             CALL EAVGSETUP(COMLYN,COMLEN)167:             CALL EAVGSETUP(COMLYN,COMLEN)
174:       ELSE IF (INDXA(COMLYN,COMLEN,'EVB') .GT. 0) THEN 
175:             CALL EVBSETUP(COMLYN,COMLEN)       
176: c------PVG RPMD CODING 
177:       ELSE IF (INDXA(COMLYN,COMLEN,'RPMD') .GT. 0) THEN 
178:             CALL RPMDSETUP(COMLYN,COMLEN)       
179:       ELSE IF (INDXA(COMLYN,COMLEN,'RPRXNCRD') .GT. 0) THEN 
180:             CALL RPRXNCRD(COMLYN,COMLEN) 
181:       ELSE IF (INDXA(COMLYN,COMLEN,'SETEVBRPMDNODES') .GT. 0) THEN 
182:             CALL SETEVBRPMDNODES(COMLYN,COMLEN)       
183: c------END PVG RPMD CODING 
184:       ELSE IF (INDXA(COMLYN,COMLEN,'EXCH') .GT. 0) THEN168:       ELSE IF (INDXA(COMLYN,COMLEN,'EXCH') .GT. 0) THEN
185:             CALL ENSREX(COMLYN,COMLEN)169:             CALL ENSREX(COMLYN,COMLEN)
186:       ELSE IF (INDXA(COMLYN,COMLEN,'WRIT') .GT. 0) THEN170:       ELSE IF (INDXA(COMLYN,COMLEN,'WRIT') .GT. 0) THEN
187:             CALL ENSWRI(COMLYN,COMLEN)171:             CALL ENSWRI(COMLYN,COMLEN)
 172: c VO 5.08
 173: ##IF STRINGM
 174:       ELSE IF (INDXA(COMLYN,COMLEN,'STRI') .GT. 0) THEN
 175:             CALL parse_string_commands(comlyn,comlen) ! parse string commands
 176: ##ENDIF
 177: c VO 5.08
188:       ELSE IF (INDXA(COMLYN,COMLEN,'SWON') .GT. 0) THEN178:       ELSE IF (INDXA(COMLYN,COMLEN,'SWON') .GT. 0) THEN
189:             JRSWAP=.TRUE.179:             JRSWAP=.TRUE.
190:             IF (IOLEV.GT.0) 180:             IF (IOLEV.GT.0) 
191:      1          WRITE(OUTU,'(A)')181:      1          WRITE(OUTU,'(A)')
192:      2             'SWAPPING OF TEMPERATURE REPLICAS ENABLED'182:      2             'SWAPPING OF TEMPERATURE REPLICAS ENABLED'
193:       ELSE IF (INDXA(COMLYN,COMLEN,'SWOFF') .GT. 0) THEN183:       ELSE IF (INDXA(COMLYN,COMLEN,'SWOFF') .GT. 0) THEN
194:             JRSWAP=.FALSE.184:             JRSWAP=.FALSE.
195:             IF (IOLEV.GT.0) 185:             IF (IOLEV.GT.0) 
196:      1          WRITE(OUTU,'(A)') 186:      1          WRITE(OUTU,'(A)') 
197:      2             'SWAPPING OF TEMPERATURE REPLICAS DISABLED'187:      2             'SWAPPING OF TEMPERATURE REPLICAS DISABLED'
943: ##INCLUDE '~/charmm_fcm/image.fcm'933: ##INCLUDE '~/charmm_fcm/image.fcm'
944: ##INCLUDE '~/charmm_fcm/comand.fcm'934: ##INCLUDE '~/charmm_fcm/comand.fcm'
945: C##INCLUDE '~/charmm_fcm/control.fcm'935: C##INCLUDE '~/charmm_fcm/control.fcm'
946:       INCLUDE 'mpif.h'936:       INCLUDE 'mpif.h'
947:       REAL*8 XOLD(*),YOLD(*),ZOLD(*),XNEW(*),YNEW(*),ZNEW(*)937:       REAL*8 XOLD(*),YOLD(*),ZOLD(*),XNEW(*),YNEW(*),ZNEW(*)
948:       REAL*8 XCOMP(*),YCOMP(*),ZCOMP(*)938:       REAL*8 XCOMP(*),YCOMP(*),ZCOMP(*)
949:       INTEGER NATOM939:       INTEGER NATOM
950: C940: C
951:       REAL*8 TMP6(6),PIST(3),PISTMP(3)941:       REAL*8 TMP6(6),PIST(3),PISTMP(3)
952:       INTEGER IERROR,I,J,K,RI,RJ,PARTNER,TMPI,II,S942:       INTEGER IERROR,I,J,K,RI,RJ,PARTNER,TMPI,II,S
953:       INTEGER MYENEL,OLDISEED943:       INTEGER MYENEL
954:       PARAMETER (MYENEL=4)944:       PARAMETER (MYENEL=4)
955:       REAL*8 MYENE(MYENEL),ENSENE(MYENEL*MAXENS)945:       REAL*8 MYENE(MYENEL),ENSENE(MYENEL*MAXENS)
956:       REAL*8 EII,EIJ,EJJ,EJI,CDEL,RNUM,SCAFAC,OLDT,DB946:       REAL*8 EII,EIJ,EJJ,EJI,CDEL,RNUM,SCAFAC,OLDT,DB
957:       REAL*8 P_I,V_I,P_J,V_J947:       REAL*8 P_I,V_I,P_J,V_J
958:       REAL*8 PIVI,PIVJ,PJVI,PJVJ948:       REAL*8 PIVI,PIVJ,PJVI,PJVJ
959:       LOGICAL QSWP,SENDFIRST,DEBUG949:       LOGICAL QSWP,SENDFIRST,DEBUG
960: C950: C
961: C  If I is this node and J is potential swap partner, then:951: C  If I is this node and J is potential swap partner, then:
962: C  MYENE(1) = EI(XI)952: C  MYENE(1) = EI(XI)
963: C  MYENE(2) = EI(XJ)953: C  MYENE(2) = EI(XJ)
964: C954: C
965:       DEBUG=.false.955:       DEBUG=.false.
966: C     ROOT decides which nodes will attempt a swap956: C     ROOT decides which nodes will attempt a swap
967:       IF (WHOIAM.EQ.0) THEN957:       IF (WHOIAM.EQ.0) THEN
968:             OLDISEED = ISEED  ! so langevin will not be affected 
969:             WRITE(OUTU,'(A)') 958:             WRITE(OUTU,'(A)') 
970:      1        ' *****************************************************'959:      1        ' *****************************************************'
971:             WRITE(OUTU,'(A)') 960:             WRITE(OUTU,'(A)') 
972:      1        ' TESTING WHETHER REPLICAS NEED TO BE SWAPPED ...'961:      1        ' TESTING WHETHER REPLICAS NEED TO BE SWAPPED ...'
973: C962: C
974:             S=INT(RANDOM(ISEED)*(ENSNSW))+1963:             S=INT(RANDOM(ISEED)*(ENSNSW))+1
975:             I=ENSISW(S)964:             I=ENSISW(S)
976:             J=ENSJSW(S)965:             J=ENSJSW(S)
977:             WRITE(OUTU,'(A,I3,A,I3)') 966:             WRITE(OUTU,'(A,I3,A,I3)') 
978:      1        ' ATTEMPTING TO SWAP REPLICAS ', I, ' AND ', J 967:      1        ' ATTEMPTING TO SWAP REPLICAS ', I, ' AND ', J 
1147:          ELSE1136:          ELSE
1148:              CALL UPDATE(COMLYN,COMLEN,XCOMP,YCOMP,ZCOMP,WMAIN,1137:              CALL UPDATE(COMLYN,COMLEN,XCOMP,YCOMP,ZCOMP,WMAIN,
1149:      1              .FALSE.,.FALSE.,1138:      1              .FALSE.,.FALSE.,
1150:      2              .TRUE.,.FALSE.,.TRUE.,0,0,0,0,0,0,0)1139:      2              .TRUE.,.FALSE.,.TRUE.,0,0,0,0,0,0,0)
1151: CC           ... failure1140: CC           ... failure
1152: C            CALL ENSSCV(XCOMP,HEAP(ENSBUF),NATOM,PARTNER-1,SENDFIRST)1141: C            CALL ENSSCV(XCOMP,HEAP(ENSBUF),NATOM,PARTNER-1,SENDFIRST)
1153: C            CALL ENSSCV(YCOMP,HEAP(ENSBUF),NATOM,PARTNER-1,SENDFIRST)1142: C            CALL ENSSCV(YCOMP,HEAP(ENSBUF),NATOM,PARTNER-1,SENDFIRST)
1154: C            CALL ENSSCV(ZCOMP,HEAP(ENSBUF),NATOM,PARTNER-1,SENDFIRST)1143: C            CALL ENSSCV(ZCOMP,HEAP(ENSBUF),NATOM,PARTNER-1,SENDFIRST)
1155:          ENDIF1144:          ENDIF
1156:       ENDIF1145:       ENDIF
1157:       IF (WHOIAM.EQ.0) THEN 
1158:             ISEED = OLDISEED  ! so langevin will not be affected 
1159:       ENDIF 
1160: ##IF ALTIX1146: ##IF ALTIX
1161:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR)1147:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR)
1162: ##ELSE1148: ##ELSE
1163:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR)1149:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR)
1164: ##ENDIF1150: ##ENDIF
1165: C1151: C
1166:       RETURN1152:       RETURN
1167:       END1153:       END
1168: C1154: C
1169:       SUBROUTINE ENSSCV(ARRAY,TMPV,NATOM,PARTNER,SENDFIRST)1155:       SUBROUTINE ENSSCV(ARRAY,TMPV,NATOM,PARTNER,SENDFIRST)
1426:       CHARACTER*(*) COMLYN1412:       CHARACTER*(*) COMLYN
1427:       INTEGER COMLEN1413:       INTEGER COMLEN
1428: C1414: C
1429: ##INCLUDE '~/charmm_fcm/dimens.fcm'1415: ##INCLUDE '~/charmm_fcm/dimens.fcm'
1430: ##INCLUDE '~/charmm_fcm/exfunc.fcm'1416: ##INCLUDE '~/charmm_fcm/exfunc.fcm'
1431: ##INCLUDE '~/charmm_fcm/stream.fcm'1417: ##INCLUDE '~/charmm_fcm/stream.fcm'
1432: ##INCLUDE '~/charmm_fcm/psf.fcm'1418: ##INCLUDE '~/charmm_fcm/psf.fcm'
1433: ##INCLUDE '~/charmm_fcm/heap.fcm'1419: ##INCLUDE '~/charmm_fcm/heap.fcm'
1434: ##INCLUDE '~/charmm_fcm/ensemble.fcm'1420: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
1435: ##INCLUDE '~/charmm_fcm/number.fcm'1421: ##INCLUDE '~/charmm_fcm/number.fcm'
1436: ##INCLUDE '~/charmm_fcm/coord.fcm' 
1437: C local:1422: C local:
1438:       INTEGER BUFLEN,I,J,JJ,NSEL1423:       INTEGER BUFLEN,I
1439:       LOGICAL JOFFSET,JCOUPLE1424:       LOGICAL JOFFSET
1440: C1425: C
1441:       QENSEXP = .true.1426:       QENSEXP = .true.
1442:       BUFLEN = NATOM*NENSEM1427:       BUFLEN = NATOM*NENSEM
1443:       IF (IOLEV.GT.0) THEN1428:       IF (IOLEV.GT.0) THEN
1444:         WRITE(OUTU,'(A)') ' WILL DO EXPONENTIAL AVERAGING OF ENERGY'1429:         WRITE(OUTU,'(A)') ' WILL DO EXPONENTIAL AVERAGING OF ENERGY'
1445:       ENDIF1430:       ENDIF
1446: C     write output?1431: C     write output?
1447:       ENSEXPU=GTRMI(COMLYN,COMLEN,'UNIT',-1)1432:       ENSEXPU=GTRMI(COMLYN,COMLEN,'UNIT',-1)
1448:       IF ((IOLEV.GT.0).AND.(ENSEXPU.GT.0)) THEN1433:       IF ((IOLEV.GT.0).AND.(ENSEXPU.GT.0)) THEN
1449:         WRITE(OUTU,'(A,I5)') ' WRITING EXP ENERGY OUTPUT TO UNIT',1434:         WRITE(OUTU,'(A,I5)') ' WRITING EXP ENERGY OUTPUT TO UNIT',
1492:           ENDIF1477:           ENDIF
1493:       ENDIF1478:       ENDIF
1494: C     allocate space for MPI receive buffers1479: C     allocate space for MPI receive buffers
1495:       ENSDX=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDX')1480:       ENSDX=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDX')
1496:       ENSDY=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDY')1481:       ENSDY=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDY')
1497:       ENSDZ=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDZ')1482:       ENSDZ=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDZ')
1498:       ENSH=ALLHP(IREAL8(NENSEM),'ensemble.src','EAVGSETUP','ENSH')1483:       ENSH=ALLHP(IREAL8(NENSEM),'ensemble.src','EAVGSETUP','ENSH')
1499: C1484: C
1500:       RETURN1485:       RETURN
1501:       END1486:       END
1502:  
1503: C1487: C
1504: C----------------------------------------------------------------------1488: C----------------------------------------------------------------------
1505:       SUBROUTINE EVBSETUP(COMLYN,COMLEN)1489:       SUBROUTINE ENSEXPAVG(DXK,DYK,DZK,HK,NATOM)
1506: C-----------------------------------------------------------------------1490: C-----------------------------------------------------------------------
1507: C the nodes share their forces and energy and take the total energy to be1491: C the nodes share their forces and energy and and do 
1508: C the lowest eigenvalue of an NEVBS  by NEVBS  molecular mechanics matrix1492: C exp(-H)=exp(-HA)+exp(-HB) averaging
1509: C coded up by david glowacki, sept 20101493: C
1510: C ... setup & allocate storage at beginning1494: C ... calculate energy & forces at each time step
1511: C1495: C
1512: ##INCLUDE '~/charmm_fcm/impnon.fcm'1496: ##INCLUDE '~/charmm_fcm/impnon.fcm'
1513:       CHARACTER*(*) COMLYN 
1514:       INTEGER COMLEN 
1515: C1497: C
 1498:       REAL*8 DXK(*),DYK(*),DZK(*),HK
 1499:       INTEGER NATOM
1516: ##INCLUDE '~/charmm_fcm/dimens.fcm'1500: ##INCLUDE '~/charmm_fcm/dimens.fcm'
1517: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
1518: ##INCLUDE '~/charmm_fcm/stream.fcm'1501: ##INCLUDE '~/charmm_fcm/stream.fcm'
1519: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1520: ##INCLUDE '~/charmm_fcm/heap.fcm'1502: ##INCLUDE '~/charmm_fcm/heap.fcm'
1521: ##INCLUDE '~/charmm_fcm/ensemble.fcm'1503: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
1522: ##INCLUDE '~/charmm_fcm/number.fcm' 
1523: ##INCLUDE '~/charmm_fcm/coord.fcm' 
1524: C local:1504: C local:
1525:       INTEGER BUFLEN,I,J,II,JJ,KK,NG,GI,NSEL,CEIDX,GIDX,MAXII,MAXJJ1505:       INTEGER BUFLEN
1526:       INTEGER EMSLCT(MAXAIM),NCOUPLE,ITMP,STATEIDX 
1527:       LOGICAL JOFFSET,JCOUPLE,FOUND,EVBCONS,EVBGAUS 
1528:       REAL*8 SMALLA, SMALLB, SMALLC,SHFTVAL 
1529:       CHARACTER*4 WRD 
1530: C 
1531:       QEVB = .true. 
1532:       BUFLEN = NATOM*NEVBS  
1533:  
1534: C     write output? 
1535:       ENSEXPU=GTRMI(COMLYN,COMLEN,'UNIT',-1) 
1536:       IF ((IOLEV.GT.0).AND.(ENSEXPU.GT.0)) THEN 
1537:         WRITE(OUTU,'(A,I5)') ' WRITING EVB ENERGY OUTPUT TO UNIT', 
1538:      1  ENSEXPU 
1539:       ENDIF 
1540:  
1541: c     read offsets, specify default offsets of zero 
1542:       IF (IOLEV.GT.0) THEN 
1543:         WRITE(OUTU,'(A)') ' ' 
1544:         WRITE(OUTU,'(A)')' READING ENERGY SHIFTS (KEYWORD SHFT)' 
1545:         WRITE(OUTU,'(A)')' IF UNSPECIFIED, THE DEFAULT SHIFT IS ZERO' 
1546:       ENDIF 
1547:       DO I=1,NEVBS  
1548:         EXPOFF(I) = 0.0 
1549:       ENDDO 
1550:  
1551: c---- get, but don't remove the next word 
1552:       WRD=CURRA4(COMLYN,COMLEN) 
1553:        
1554:       DO WHILE (WRD.EQ.'SHFT') 
1555:  
1556: c------ if the word was 'SHFT', remove it and read the following input 
1557:         WRD=NEXTA4(COMLYN,COMLEN) 
1558:         STATEIDX=NEXTI(COMLYN,COMLEN) 
1559:         SHFTVAL=NEXTF(COMLYN,COMLEN)         
1560:        
1561:         IF(STATEIDX.LT.1)THEN 
1562:           WRITE(OUTU,'(A,I3,A)')' STATE INDEX ', STATEIDX ,' FOR COUPLI 
1563:      &NG IS LESS THAN 1' 
1564:           CALL WRNDIE(-1,'ENSEMBLE>','CHECK YOUR INPUT!') 
1565:         ELSEIF(STATEIDX.GT.NEVBS) THEN 
1566:           WRITE(OUTU,'(A,I3,A,I3)')' SHIFT STATE INDEX ', STATEIDX ,' IS 
1567:      &GREATER THAN THE MAXIMUM ALLOWABLE STATE INDEX ', NEVBS   
1568:           CALL WRNDIE(-1,'ENSEMBLE>','CHECK YOUR INPUT!')           
1569:         ELSE 
1570:           EXPOFF(STATEIDX)=SHFTVAL 
1571:           WRITE(OUTU,'(A,I3,A,F12.5)') ' STATE ', STATEIDX,' WILL BE  
1572:      &SHIFTED BY ', EXPOFF(STATEIDX) 
1573:         ENDIF 
1574:          
1575:         WRD=CURRA4(COMLYN,COMLEN)        
1576:  
1577:       ENDDO 
1578:        
1579:       IF (IOLEV.GT.0) THEN 
1580:         WRITE(OUTU,'(A)') ' ' 
1581:         WRITE(OUTU,'(A)') ' -----------------------------------------' 
1582:         WRITE(OUTU,'(A)') ' ENERGY SHIFT SUMMARY: ' 
1583:         DO I=1,NEVBS 
1584:           WRITE(OUTU,'(A,I3,A,F12.5)') ' STATE = ', I,' SHIFT = ',  
1585:      &EXPOFF 
1586:      &(I) 
1587:         ENDDO 
1588:         WRITE(OUTU,'(A)') ' -----------------------------------------' 
1589:         WRITE(OUTU,'(A)') ' ' 
1590:       ENDIF 
1591:        
1592: c---- davidglo, read in the matrix coupling elements 
1593:       NCOUPLE=NEVBS*(NEVBS-1)/2 
1594:       MAXII=NEVBS   
1595:       MAXJJ=NEVBS   
1596:        
1597:       IF (IOLEV.GT.0) THEN 
1598:         WRITE(OUTU,'(A)') ' USING MULTI-REF MM (EVB)' 
1599:         WRITE(OUTU,'(A,I3,A)') ' ENSEMBLE EVB SETTING ALL ',NCOUPLE,'  
1600:      &ELEM 
1601:      &ENTS OF THE COUPLING MATRIX TO A DEFAULT VALUE OF ZERO' 
1602:         IF (IOLEV.GT.0) THEN 
1603:           WRITE(OUTU,'(A,I3,A,I3,A)')' ENERGY WILL BE THE LOWEST  
1604:      &EIGENVALUE OF A ',NEVBS,' BY ',NEVBS,' COUPLING MATRIX' 
1605:         ENDIF 
1606:       ENDIF 
1607: c---- nread keeps track of how many different times 'couple' has been read 
1608:       NREAD=0     
1609:       DO WHILE ( INDXA(COMLYN, COMLEN, 'COUPLE') .GT. 0 )  
1610:  
1611:         WRITE(OUTU,'(A)') ' ' 
1612:         NREAD=NREAD+1 
1613:         II = NEXTI(COMLYN,COMLEN) 
1614:         JJ = NEXTI(COMLYN,COMLEN) 
1615:          
1616: c        write(*,*)'II,jj', II, JJ 
1617:  
1618: c------ check upper bounds of INDICES 
1619:         IF(II.GT.MAXII) THEN 
1620:           WRITE(OUTU,'(A,I3,A,I3)')'THE MATRIX ELEMENT INDEX ', II,' IS  
1621:      &LARGER THAN THE MAX STATE INDEX ',MAXII 
1622:           CALL WRNDIE(-1,'ENSEMBLE>','MATRIX ELEMENT INDEX TOO LARGE!') 
1623:         ELSEIF(JJ.GT.MAXJJ) THEN 
1624:           WRITE(OUTU,'(A,I3,A,I3)')'THE MATRIX ELEMENT INDEX ', JJ,' IS  
1625:      &LARGER THAN THE MAX STATE INDEX ',MAXJJ 
1626:           CALL WRNDIE(-1,'ENSEMBLE>','MATRIX ELEMENT INDEX TOO LARGE!')         
1627:  
1628: c------ check lower bounds of indices 
1629:         ELSEIF(II.LT.1) THEN 
1630:           WRITE(OUTU,'(A,I3,A)')'THE MATRIX ELEMENT INDEX ', II,' IS LES 
1631:      &S THAN ONE' 
1632:           CALL WRNDIE(-1,'ENSEMBLE>','MATRIX ELEMENT LESS THAN ONE!') 
1633:         ELSEIF(JJ.LT.1) THEN 
1634:           WRITE(OUTU,'(A,I3,A)')'THE MATRIX ELEMENT INDEX ', JJ,' IS LES 
1635:      &S THAN ONE' 
1636:           CALL WRNDIE(-1,'ENSEMBLE>','MATRIX ELEMENT LESS THAN ONE!') 
1637:  
1638: c------ read the coupling elements in upper tridiagonal form, which means that ii < jj 
1639:         ELSEIF(II.EQ.JJ) THEN 
1640:           WRITE(OUTU,'(A,I3,A,I3,A)')'THE SPECIFIED MATRIX INDICES ', II 
1641:      &,' AND ', JJ, 'ARE EQUAL' 
1642:           CALL WRNDIE(-1,'ENSEMBLE>','MATRIX ELEMENT INDICES ARE EQUAL') 
1643:         ENDIF 
1644:         IF(II.GT.JJ) THEN 
1645:           ITMP=II 
1646:           II=JJ 
1647:           JJ=ITMP 
1648:         ENDIF 
1649:  
1650: c------ EVBIIDX & EVBJIDX are two arrays that keep track of the coupling element  
1651: c------ to which the relevant parameters belong 
1652:  
1653:  
1654:         EVBIIDX(NREAD)=II 
1655:         EVBJIDX(NREAD)=JJ 
1656:         EVBCONS=.FALSE. 
1657:         EVBGAUS=.FALSE. 
1658:         ONEDCRD=.FALSE. 
1659:         TWODCRD=.FALSE. 
1660:         IF (IOLEV.GT.0) THEN 
1661:           WRITE(OUTU,'(A,I3,A,I3,A)')' READING TERMS FOR INCLUSION IN  
1662:      &EVB COUPLING ELEMENT ',II,'-',JJ,' ...' 
1663:         ENDIF 
1664:          
1665: c------ determine whether it's gonna be a constant or gaussian we're reading 
1666:         WRD=NEXTA4(COMLYN,COMLEN) 
1667:         IF(WRD.EQ.'CONS') THEN 
1668:           EVBCONS=.TRUE. 
1669:         ELSEIF (WRD.EQ.'GAUS') THEN 
1670:           EVBGAUS=.TRUE. 
1671:         ENDIF 
1672:  
1673:         IF (EVBCONS) THEN 
1674:          COUPRM(NREAD,1) = NEXTF(COMLYN,COMLEN) 
1675:          IF (IOLEV.GT.0) THEN 
1676:            WRITE(OUTU,'(A,F12.5)')' SPECIFIED TERM IS A CONSTANT WITH  
1677:      &VALUE ',COUPRM(NREAD,1) 
1678:          ENDIF 
1679:          COUTYPE(NREAD)='CONS' 
1680:           
1681: c------ if it's a gaussian we're reading - need to check if ONED or TWOD or DIFF 
1682:         ELSEIF(EVBGAUS) THEN 
1683:          WRD=NEXTA4(COMLYN,COMLEN) 
1684:          IF(WRD.EQ.'ONED') THEN 
1685:            ONEDCRD=.TRUE. 
1686:          ELSEIF(WRD.EQ.'TWOD') THEN 
1687:            TWODCRD=.TRUE. 
1688:          ELSEIF(WRD.EQ.'DIFF') THEN 
1689:            DIFFCRD=.TRUE. 
1690:          ENDIF 
1691:  
1692:          IF(ONEDCRD) THEN 
1693: C          ONEDCRD=.TRUE. 
1694:            COUTYPE(NREAD)='ONED' 
1695:            IF (IOLEV.GT.0) THEN 
1696:              WRITE(OUTU,'(A)')' SPECIFIED TERM IS A ONED GAUSSIAN' 
1697:              WRITE(OUTU,'(A)')' IT DEPENDS ON A DISTANCE R1:' 
1698:              WRITE(OUTU,'(A)')' A*EXP(-(R1-R01)**2/(2C**2))' 
1699:            ENDIF 
1700:          ELSEIF(TWODCRD) THEN 
1701:            IF (IOLEV.GT.0) THEN 
1702: C          TWODCRD=.TRUE. 
1703:            COUTYPE(NREAD)='TWOD' 
1704:              WRITE(OUTU,'(A)')' SPECIFIED TERM IS A TWOD GAUSSIAN' 
1705:              WRITE(OUTU,'(A)')' IT DEPENDS ON TWO DISTANCES R1 & R2:' 
1706:              WRITE(OUTU,'(A)')' A*EXP(-(a[R1-R01]**2 +  
1707:      &2b[R1-R01][R2-R02] + c[R2-R02]**2))' 
1708:              WRITE(OUTU,'(A)')' where:' 
1709:              WRITE(OUTU,'(A)')' a = [COS(THETA)]**2/(2*C1**2) +  
1710:      &[SIN(THETA)]**2/(2*C2**2) ' 
1711:              WRITE(OUTU,'(A)')' b = [SIN(THETA)]**2/(2*C1**2) +  
1712:      &[COS(THETA)]**2/(2*C2**2) ' 
1713:              WRITE(OUTU,'(A)')' c = SIN(2*THETA)/(4*C1**2) +  
1714:      &[SIN(2*THETA)]**2/(4*C2**2) ' 
1715:            ENDIF 
1716:          ELSEIF(DIFFCRD) THEN 
1717: C          DIFFCRD=.TRUE. 
1718:            COUTYPE(NREAD)='DIFF' 
1719:            IF (IOLEV.GT.0) THEN 
1720:              WRITE(OUTU,'(A)')' SPECIFIED TERM IS A DIFFERENCE  
1721:      &GAUSSIAN' 
1722:              WRITE(OUTU,'(A)')' IT DEPENDS ON THE DIFF OF TWO  
1723:      &DISTANCES:' 
1724:              WRITE(OUTU,'(A)')' A*EXP(-((DIFF)-DIFF0)**2/(2C**2)) WHERE  
1725:      &DIFF=R1-R2' 
1726:            ENDIF 
1727:          ELSE 
1728:            IF (IOLEV.GT.0) THEN 
1729:              WRITE(OUTU,'(A)')' UNRECOGNIZED INPUT OPTION FOR GAUSSIAN  
1730:      &COUPLING ELEMENT... EXITING' 
1731:              CALL WRNDIE(-1,'ENSEMBLE>','FOR EVB WITH GAUSSIAN COUPLING 
1732:      &YOU MUST SPECIFY KEYWORD ONED OR TWOD') 
1733:            ENDIF 
1734:          ENDIF 
1735:           
1736: c------- now read the parameters for the first gaussian, used by ONED, DIFF, TWOD 
1737:          IF (IOLEV.GT.0) THEN 
1738:            WRITE(OUTU,'(A)') ' READING A, R01, C1 PARAMETERS' 
1739:          ENDIF 
1740:          COUPRM(NREAD,1) = NEXTF(COMLYN,COMLEN) 
1741:          IF (IOLEV.GT.0) THEN 
1742:            WRITE(OUTU,'(A,F12.5)')' A   : ',COUPRM(NREAD,1) 
1743:          ENDIF 
1744:          COUPRM(NREAD,2) = NEXTF(COMLYN,COMLEN) 
1745:          IF (IOLEV.GT.0) THEN 
1746:            WRITE(OUTU,'(A,F12.5)')' R01 : ',COUPRM(NREAD,2) 
1747:          ENDIF 
1748:          COUPRM(NREAD,3) = NEXTF(COMLYN,COMLEN) 
1749:          IF (IOLEV.GT.0) THEN 
1750:            WRITE(OUTU,'(A,F12.5)')' C1  : ',COUPRM(NREAD,3) 
1751:          ENDIF 
1752: c------- read in the first atom selection; this is always required 
1753:          EMSLCT=0 
1754:          IF(ONEDCRD.OR.TWODCRD.OR.DIFFCRD) THEN 
1755:            IF (IOLEV.GT.0) THEN 
1756:              WRITE(OUTU,'(A)') ' READING ATOM SELECTION TO DEFINE R1' 
1757:            ENDIF 
1758:            CALL SELCTA(COMLYN,COMLEN,EMSLCT,X,Y,Z,WMAIN,.TRUE.) 
1759:            NSEL = 0 
1760:            DO II=1,MAXAIM 
1761:              IF (EMSLCT(II) .EQ. 1) THEN 
1762:                NSEL = NSEL + 1 
1763:                ESEL(NREAD,NSEL) = II 
1764:              ENDIF 
1765:            ENDDO 
1766:            IF(NSEL.LT.2.OR.NSEL.GT.2) THEN 
1767:               WRITE(OUTU,'(A)') ' EVB WITH GAUSSIAN COUPLING ELEMENTS RE 
1768:      &QUIRES SPECIFICATION OF 2 ATOMS' 
1769:              CALL WRNDIE(-1,'ENSEMBLE>','FOR EVB WITH GAUSSIAN COUPLING 
1770:      & YOU MUST SELECT 2 ATOMS') 
1771:            ELSE 
1772:              IF (IOLEV.GT.0) THEN 
1773:                WRITE(OUTU,'(A)') ' ATOM SELECTIONS:' 
1774:              ENDIF 
1775:              DO KK=1,NSEL 
1776:                IF (IOLEV.GT.0) THEN 
1777:                  WRITE(OUTU,'(A,I3)')' ATOM ',ESEL(NREAD,KK) 
1778:                ENDIF 
1779:              ENDDO 
1780:            ENDIF 
1781:          ENDIF 
1782: c------- reset the elements of EMSLCT to zero 
1783:          EMSLCT=0 
1784: c------- if there's no second selection, set ESEL2 elements to zero 
1785:          IF(ONEDCRD) THEN 
1786:           ESEL2(NREAD,1)=0 
1787:           ESEL2(NREAD,2)=0 
1788:          ELSEIF(DIFFCRD) THEN 
1789: c-------- read in the second atom selection 
1790:           IF (IOLEV.GT.0) THEN 
1791:             WRITE(OUTU,'(A)') ' ATOM SELECTION TO DEFINE R2:' 
1792:           ENDIF 
1793:           CALL SELCTA(COMLYN,COMLEN,EMSLCT,X,Y,Z,WMAIN,.TRUE.) 
1794:           NSEL = 0 
1795:           DO II=1,MAXAIM 
1796:             IF (EMSLCT(II) .EQ. 1) THEN 
1797:               NSEL = NSEL + 1 
1798:               ESEL2(NREAD,NSEL) = II 
1799:             ENDIF 
1800:           ENDDO 
1801:           IF(NSEL.LT.2.OR.NSEL.GT.2) THEN 
1802:             WRITE(OUTU,'(A)') ' EVB WITH GAUSSIAN COUPLING ELEMENTS RE 
1803:      &QUIRES SPECIFICATION OF 2 ATOMS' 
1804:             CALL WRNDIE(-1,'ENSEMBLE>','FOR EVB WITH GAUSSIAN COUPLING 
1805:      & YOU MUST SELECT 2 ATOMS') 
1806:           ELSE 
1807:             WRITE(OUTU,'(A)')' ATOM SELECTIONS:' 
1808:             DO KK=1,NSEL 
1809:               IF (IOLEV.GT.0) THEN 
1810:                 WRITE(OUTU,'(A,I3)')' ATOM ',ESEL2(NREAD,KK) 
1811:               ENDIF 
1812:             ENDDO 
1813:          ENDIF 
1814:            
1815: c------- now read the additional TWOD gaussian parameters: theta, r02, c2 
1816:          ELSEIF(TWODCRD) THEN 
1817:           IF (IOLEV.GT.0) THEN 
1818:             WRITE(OUTU,'(A)') ' READING THETA, R02, C2 PARAMETERS' 
1819:           ENDIF 
1820:           COUPRM(NREAD,4) = NEXTF(COMLYN,COMLEN) 
1821:           IF (IOLEV.GT.0) THEN 
1822:             WRITE(OUTU,'(A,F12.5)') ' THETA : ',COUPRM(NREAD,4) 
1823:           ENDIF 
1824:           COUPRM(NREAD,5) = NEXTF(COMLYN,COMLEN) 
1825:           IF (IOLEV.GT.0) THEN 
1826:             WRITE(OUTU,'(A,F12.5)') ' R02   : ',COUPRM(NREAD,5) 
1827:           ENDIF 
1828:           COUPRM(NREAD,6) = NEXTF(COMLYN,COMLEN) 
1829:           IF (IOLEV.GT.0) THEN 
1830:             WRITE(OUTU,'(A,F12.5)') ' C2    : ',COUPRM(NREAD,6) 
1831:           ENDIF 
1832: c-------- how to calculate the values of a, b, and c 
1833: c         SMALLA = COS(THETA)**2/(2*SIGR1**2) 
1834: c     &          + SIN(THETA)**2/(2*SIGR2**2) 
1835: c         SMALLB = -1*SIN(2.0D0*THETA)/(4.0D0*SIGR1**2)  
1836: c     &          +  SIN(2*THETA)/(4*SIGR2**2) 
1837: c         SMALLC = SIN(THETA)**2/(2*SIGR1**2) 
1838: c     &          + COS(THETA)**2/(2*SIGR2**2) 
1839:  
1840: c-------- calculate a, b, c parameters, held in COUPRM(NREAD,7-9) 
1841:           SMALLA=((COS(COUPRM(NREAD,4)))**2)/(2.0D0*COUPRM(NREAD,3)**2) 
1842:      &          +((SIN(COUPRM(NREAD,4)))**2)/(2.0D0*COUPRM(NREAD,6)**2) 
1843:  
1844:           SMALLB=-1.0D0*SIN(2.0D0*COUPRM(NREAD,4))/ 
1845:      &                (4.0D0*COUPRM(NREAD,3)**2) 
1846:      &          +SIN(2.0D0*COUPRM(NREAD,4))/(4.0D0*COUPRM(NREAD,6)**2) 
1847:  
1848:           SMALLC=((SIN(COUPRM(NREAD,4)))**2)/(2.0D0*COUPRM(NREAD,3)**2) 
1849:      &          +((COS(COUPRM(NREAD,4)))**2)/(2.0D0*COUPRM(NREAD,6)**2) 
1850:       
1851:           COUPRM(NREAD,7)=SMALLA 
1852:           COUPRM(NREAD,8)=SMALLB 
1853:           COUPRM(NREAD,9)=SMALLC 
1854:       
1855:           IF (IOLEV.GT.0) THEN 
1856:             WRITE(OUTU,'(A,F12.5)')' DERIVED PARAMETER a =',SMALLA 
1857:             WRITE(OUTU,'(A,F12.5)')' DERIVED PARAMETER b =',SMALLB 
1858:             WRITE(OUTU,'(A,F12.5)')' DERIVED PARAMETER c =',SMALLC 
1859:           ENDIF 
1860: c-------- read in the second atom selection 
1861:           IF (IOLEV.GT.0) THEN 
1862:             WRITE(OUTU,'(A)') ' ATOM SELECTION TO DEFINE R2:' 
1863:           ENDIF 
1864:           CALL SELCTA(COMLYN,COMLEN,EMSLCT,X,Y,Z,WMAIN,.TRUE.) 
1865:           NSEL = 0 
1866:           DO II=1,MAXAIM 
1867:             IF (EMSLCT(II) .EQ. 1) THEN 
1868:               NSEL = NSEL + 1 
1869:               ESEL2(NREAD,NSEL) = II 
1870:             ENDIF 
1871:           ENDDO 
1872:           IF(NSEL.LT.2.OR.NSEL.GT.2) THEN 
1873:             WRITE(OUTU,'(A)') ' EVB WITH GAUSSIAN COUPLING ELEMENTS RE 
1874:      &QUIRES SPECIFICATION OF 2 ATOMS' 
1875:             CALL WRNDIE(-1,'ENSEMBLE>','FOR EVB WITH GAUSSIAN COUPLING 
1876:      & YOU MUST SELECT 2 ATOMS') 
1877:           ELSE 
1878:             WRITE(OUTU,'(A)')' ATOM SELECTIONS:' 
1879:             DO KK=1,NSEL 
1880:               WRITE(OUTU,'(A,I3)')' ATOM ',ESEL2(NREAD,KK) 
1881:             ENDDO 
1882:           ENDIF 
1883:          ENDIF 
1884:         ENDIF 
1885:       ENDDO  
1886:  
1887: c      DO I=1,NREAD 
1888: c        WRITE(*,*)'I,EVBIIDX,EVBJIDX,COUTYPE',I,EVBIIDX(I),EVBJIDX(I),CO 
1889: c     &UTYPE(I) 
1890: c      ENDDO 
1891:  
1892: c---- write out the matrix coupling elements to check input is correct 
1893:       IF (IOLEV.GT.0) THEN 
1894:         WRITE(OUTU,'(A)') ' ' 
1895:         WRITE(OUTU,'(A)') ' -----------------------------------------' 
1896:         WRITE(OUTU,'(A)') ' COUPLING MATRIX ELEMENT SUMMARY: ' 
1897:       ENDIF 
1898:       DO II=1,NEVBS  
1899:         DO JJ=II+1,NEVBS 
1900:           FOUND=.FALSE. 
1901:           IF (IOLEV.GT.0) THEN 
1902:             WRITE(OUTU,'(A)') ' ' 
1903:             WRITE(OUTU,'(A,I3,A,I3,A)')' ELEMENT WHICH COUPLES  
1904:      &STATES ', II,' & ', JJ,' IS A SUM OF THE FOLLOWING TERMS: ' 
1905:           ENDIF 
1906:           DO J=1,NREAD 
1907:             IF((II).EQ.EVBIIDX(J).AND.(JJ).EQ.EVBJIDX(J)) THEN 
1908:               FOUND=.TRUE. 
1909:               IF(COUTYPE(J).EQ.'CONS') THEN 
1910:                 IF (IOLEV.GT.0) THEN 
1911:                   WRITE(OUTU,'(A,F12.5)')' >> CONSTANT: ',COUPRM(J,1) 
1912:                 ENDIF 
1913:               ELSEIF(COUTYPE(J).EQ.'ONED') THEN 
1914:                 IF (IOLEV.GT.0) THEN 
1915:                   WRITE(OUTU,'(A,F12.5,F12.5,F12.5)')' >> ONED  
1916:      &GAUSSIAN: ',COUPRM(J,1),COUPRM(J,2),COUPRM(J,3) 
1917:                   WRITE(OUTU,'(A,I3,A,I3)')'    DEPENDS ON DISTANCE: ', 
1918:      &ESEL(J,1),' - ',ESEL(J,2) 
1919:                 ENDIF 
1920:               ELSEIF(COUTYPE(J).EQ.'DIFF') THEN 
1921:                 IF (IOLEV.GT.0) THEN 
1922:                   WRITE(OUTU,'(A,F12.5,F12.5,F12.5)')' >> DIFF  
1923:      &GAUSSIAN: ',COUPRM(J,1),COUPRM(J,2),COUPRM(J,3) 
1924:                   WRITE(OUTU,'(A,I3,A,I3,A,I3,A,I3)')'    DEPENDS ON  
1925:      &TWO DISTANCES: ',ESEL(J,1),' - ',ESEL(J,2),' AND ',ESEL2(J,1), 
1926:      &' - ',ESEL2(J,2)          
1927:                 ENDIF        
1928:               ELSEIF(COUTYPE(J).EQ.'TWOD') THEN 
1929:                 IF (IOLEV.GT.0) THEN 
1930:                   WRITE(OUTU,'(A,F12.5,F12.5,F12.5,F12.5,F12.5,F12.5)')' 
1931:      &>> TWOD GAUSSIAN: ',COUPRM(J,1),COUPRM(J,2),COUPRM(J,3), 
1932:      &COUPRM(J,4),COUPRM(J,5),COUPRM(J,6) 
1933:                   WRITE(OUTU,'(A,I3,A,I3,A,I3,A,I3)')'    DEPENDS ON  
1934:      &TWO DISTANCES: ',ESEL(J,1),' - ',ESEL(J,2),' AND ',ESEL2(J,1), 
1935:      &' - ',ESEL2(J,2)                
1936:                 ENDIF 
1937:               ENDIF      
1938:             ENDIF 
1939:           ENDDO 
1940:  
1941:           IF(.NOT.FOUND) THEN 
1942:             WRITE(OUTU,'(A)')' DEFAULT VALUE: 0.0000 '             
1943:           ENDIF 
1944:  
1945:         ENDDO 
1946:       ENDDO 
1947:       WRITE(OUTU,'(A)') ' ---------------------------------------------' 
1948:       WRITE(OUTU,'(A)') ' ' 
1949:          
1950: c---- davidglo, end block for reading matrix coupling elements 
1951:  
1952: C     allocate space for MPI receive buffers 
1953:       ENSDX=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDX') 
1954:       ENSDY=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDY') 
1955:       ENSDZ=ALLHP(IREAL8(BUFLEN),'ensemble.src','EAVGSETUP','ENSDZ') 
1956:       ENSH=ALLHP(IREAL8(NENSEM),'ensemble.src','EAVGSETUP','ENSH') 
1957: C1506: C
 1507:       BUFLEN = NATOM*NENSEM
 1508:       CALL ENSEXPAVG2(DXK,DYK,DZK,HK,NATOM,HEAP(ENSDX),HEAP(ENSDY),
 1509:      1        HEAP(ENSDZ),HEAP(ENSH),BUFLEN)
1958:       RETURN1510:       RETURN
1959:       END1511:       END
1960: 1512: C
1961: 1513:       SUBROUTINE ENSEXPAVG2(DXK,DYK,DZK,HK,NATOM,
1962: C----------------------------------------------------------------------1514:      1                      RCVDX,RCVDY,RCVDZ,RCVH,BUFLEN)
1963:       SUBROUTINE SETEVBRPMDNODES(COMLYN,COMLEN) 
1964: C-----------------------------------------------------------------------1515: C-----------------------------------------------------------------------
1965: C This subroutine manages the mpi process assignment of system replicas1516: C the nodes share their forces and energy and and do exp(-H)=exp(-HA)+exp(-HB)
1966: C when using ring polymer molecular dynamics (RPMD) with 1517: C averaging
1967: C empirical valence bond (EVB). For a simulation with NRPBDS1518: C
1968: C RPMD beads and NEVBS evb states the mpi processes are divided up1519: C ... calculate energy & forces at each time step
1969: C as follows (for a 4 RPMD bead and 2 evb state simulation: 
1970: C 
1971: C MPI Process:   EVB state:   BEAD: 
1972: C      0           1           1           
1973: C      1           1           2 
1974: C      2           1           3 
1975: C      3           1           4 
1976: C      4           2           1 
1977: C      5           2           2 
1978: C      6           2           3 
1979: C      7           2           4 
1980: C1520: C
1981:  
1982: C Writen by Patrick von Glehn Sept 2013 
1983: ##INCLUDE '~/charmm_fcm/impnon.fcm'1521: ##INCLUDE '~/charmm_fcm/impnon.fcm'
1984: ##INCLUDE '~/charmm_fcm/dimens.fcm'1522: C
1985: ##INCLUDE '~/charmm_fcm/exfunc.fcm'1523:       REAL*8 DXK(*),DYK(*),DZK(*),RCVDX(*),RCVDY(*),RCVDZ(*)
 1524:       REAL*8 RCVH(*),HK
 1525:       INTEGER NATOM,BUFLEN
1986: ##INCLUDE '~/charmm_fcm/stream.fcm'1526: ##INCLUDE '~/charmm_fcm/stream.fcm'
1987: ##INCLUDE '~/charmm_fcm/heap.fcm' 
1988: ##INCLUDE '~/charmm_fcm/ensemble.fcm'1527: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
 1528: ##INCLUDE '~/charmm_fcm/contrl.fcm'
1989: ##INCLUDE '~/charmm_fcm/number.fcm'1529: ##INCLUDE '~/charmm_fcm/number.fcm'
1990:       INCLUDE 'mpif.h'1530:       INCLUDE 'mpif.h'
1991:       CHARACTER*(*) COMLYN1531: C     local
1992:       INTEGER COMLEN1532:       INTEGER IERROR,I,M
1993:       integer, allocatable, dimension(:) :: interStateGroupRanks(:)1533:       REAL*8 EXPSUM,EXP2SUM,DXI,DYI,DZI,MINE
1994:       integer, allocatable, dimension(:) :: interBeadGroupRanks(:)1534:       REAL*8 DXI2,DYI2,DZI2
1995:       integer, allocatable, dimension(:) :: rcvInterBeadRanks(:)1535:       REAL*8 TMPH(1)
1996:       integer, allocatable, dimension(:) :: rcvInterStateRanks(:)1536:       REAL*8 EXPWT(MAXENS)
1997:       integer ierr,ierror1537: C
1998:       integer orig_group1538:       TMPH(1) = HK
1999:       integer XX 
2000:  
2001:       INTEGER REPLICAID, NREP, II, JJ 
2002:  
2003:       QEVBRP = .True. 
2004:  
2005: c-----initialize to default values 
2006:       NEVBS=1 
2007:       NRPBDS=1 
2008:  
2009:       IF ( INDXA(COMLYN, COMLEN, 'NEVBSTATES') .GT. 0 ) THEN 
2010:         NEVBS=NEXTF(COMLYN,COMLEN) 
2011:         IF (WHOIAM .EQ. 0) THEN 
2012:           WRITE(OUTU,'(A,I5)')'NUMBER OF EVB STATES =',NEVBS 
2013:         ENDIF 
2014:       ENDIF 
2015:       IF ( INDXA(COMLYN, COMLEN, 'NRPMDBEADS') .GT. 0 ) THEN 
2016:         NRPBDS=NEXTF(COMLYN,COMLEN) 
2017:         IF (WHOIAM .EQ. 0) THEN 
2018:           WRITE(OUTU,'(A,I5)')'NUMBER OF RING POLYEMER MOLECULAR DYNAMICS 
2019:      &BEADS =',NRPBDS 
2020:         ENDIF 
2021:       ENDIF 
2022:  
2023: c-----check that the number of replicas = number of mpi processes requested 
2024:       NREP = NRPBDS*NEVBS 
2025:       IF (WHOIAM .EQ. 0) THEN 
2026:         WRITE(OUTU,'(A,I5)')'NUMBER OF MPI PROCESSES CALLED:',NENSEM 
2027:         WRITE(OUTU,'(A,I5)')'NUMBER OF REPLICAS:',NREP 
2028:         IF (NREP.NE.NENSEM) THEN 
2029:           CALL WRNDIE(-1,'ENSEMBLE>','NUMBER OF MPI PROCESSES REQUESTED 
2030:      &DOES NOT EQUAL NUMBER OF REPLICAS SPECIFIED IN INPUT FILE')  
2031:         ENDIF 
2032:       ENDIF 
2033:  
2034: c-----build 2D array REPIDS(STATEID,BEADID) 
2035:       REPLICAID=1 
2036:       DO II=1, NEVBS 
2037:         DO JJ=1, NRPBDS 
2038:           REPIDS(REPLICAID,1)=II 
2039:           REPIDS(REPLICAID,2)=JJ     
2040:           IF (WHOIAM .EQ. 0) THEN 
2041:             WRITE(OUTU,'(A,I5,A,I5,A,I5)')' REPLICA ',REPLICAID,'  
2042:      &EVBSTATE=',REPIDS(REPLICAID,1),' BEADID= ', 
2043:      &REPIDS(REPLICAID,2) 
2044:           ENDIF 
2045:           REPLICAID=REPLICAID + 1 
2046:         ENDDO 
2047:       ENDDO 
2048:  
2049:       STATEID=REPIDS(WHOIAM+1,1) 
2050:       CALL SETMSI('STATEID',STATEID) 
2051:       BEADID=REPIDS(WHOIAM+1,2) 
2052:       CALL SETMSI('BEADID',BEADID) 
2053: ce---PVG error checking 
2054: ce      WRITE(OUTU,'(A,I,A,I,A,I)')'I AM NODE',WHOIAM,' MY EVB STATE IS', 
2055: ce     &STATEID,' MY BEAD NUMBER IS',BEADID 
2056: ce---PVG end error checking           
2057:  
2058:  
2059: ##IF ALTIX1539: ##IF ALTIX
2060:       CALL MPI8COMM_RANK(MPI_COMM_WORLD, WHOIAM, IERROR)1540:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR)
2061:       CALL MPI8COMM_GROUP(MPI_COMM_WORLD, orig_group, IERROR)1541:       CALL MPI8ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1),
 1542:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
 1543:       CALL MPI8ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1),
 1544:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
 1545:       CALL MPI8ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1),
 1546:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
 1547:       CALL MPI8ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1),
 1548:      1       1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
2062: ##ELSE1549: ##ELSE
2063:       CALL MPI_COMM_RANK(MPI_COMM_WORLD, WHOIAM, IERROR)1550:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR)
2064:       CALL MPI_COMM_GROUP(MPI_COMM_WORLD, orig_group, IERROR)1551:       CALL MPI_ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1),
 1552:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
 1553:       CALL MPI_ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1),
 1554:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
 1555:       CALL MPI_ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1),
 1556:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
 1557:       CALL MPI_ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1),
 1558:      1       1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR)
2065: ##ENDIF1559: ##ENDIF
2066: 1560: C     first calculate mean to use as an offset to avoid numerical
2067:       allocate(interStateGroupRanks(NEVBS))1561: C     errors in exp(-beta U) for very small (large negative) U.
2068:       allocate(rcvInterStateRanks(NEVBS))1562:       MINE=RCVH(1)+EXPOFF(1)
2069:       allocate(rcvInterBeadRanks(NRPBDS))1563:       DO I=2,NENSEM
2070:       allocate(interBeadGroupRanks(NRPBDS))1564:          MINE = MIN(MINE,RCVH(I)+EXPOFF(I))
2071:  
2072: c-----Set up the MPI group containing all of the states 
2073: c-----associated with a the same RPMD bead as a particular replica 
2074:       DO II=1,NEVBS 
2075:          interStateGroupRanks(II) = (II-1)*NRPBDS + MOD(WHOIAM,NRPBDS) 
2076:       ENDDO1565:       ENDDO
2077: ##IF ALTIX1566: C     convert entries HK(I) to EXP(-beta*HK(I))
2078:       call MPI8COMM_GROUP(MPI_COMM_WORLD, orig_group, ierr)1567:       EXPSUM=0.0
2079:       call MPI8GROUP_INCL(orig_group, NEVBS, interStateGroupRanks,1568:       EXP2SUM=0.0
2080:      &                  interStateGroup, ierr)1569:       DO I=1,NENSEM
2081:       call MPI8COMM_CREATE(MPI_COMM_WORLD, interStateGroup,1570:          EXPWT(I)=DEXP(-ENSEBETA*(RCVH(I)+EXPOFF(I)-MINE))
2082:      &                  interStateGroupComm, ierr)1571:          EXPSUM=EXPSUM+EXPWT(I)
2083:       call MPI8ALLGATHER(WHOIAM,1,MPI_INTEGER,rcvInterStateRanks(1),1572:          IF (QEXPBW) EXP2SUM=EXP2SUM+(EXPWT(I))**2
2084:      &       1,MPI_INTEGER,interStateGroupComm,IERROR) 
2085: ##ELSE 
2086:       call MPI_COMM_GROUP(MPI_COMM_WORLD, orig_group, ierr) 
2087:       call MPI_GROUP_INCL(orig_group, NEVBS, interStateGroupRanks, 
2088:      &                  interStateGroup, ierr) 
2089:       call MPI_COMM_CREATE(MPI_COMM_WORLD, interStateGroup, 
2090:      &                  interStateGroupComm, ierr) 
2091:       call MPI_ALLGATHER(WHOIAM,1,MPI_INTEGER,rcvInterStateRanks(1), 
2092:      &       1,MPI_INTEGER,interStateGroupComm,IERROR) 
2093: ##ENDIF 
2094:  
2095: c----Set up the MPI group containing each replica with the same state 
2096: c----as a particular replica 
2097:  
2098:       DO II=1,NEVBS 
2099:         DO JJ=1,NRPBDS 
2100:           if (WHOIAM .GE. (II-1)*NRPBDS) then 
2101:             interBeadGroupRanks(JJ) = (II-1)*NRPBDS + JJ - 1 
2102:           endif 
2103:         ENDDO 
2104:       ENDDO 
2105:  
2106: ##IF ALTIX 
2107:  
2108:       call MPI8GROUP_INCL(orig_group, NRPBDS, interBeadGroupRanks, 
2109:      &                 interBeadGroup, ierr) 
2110:       call MPI8COMM_CREATE(MPI_COMM_WORLD, interBeadGroup, 
2111:      &                 interBeadGroupComm, ierr) 
2112:       call MPI8ALLGATHER(WHOIAM,1,MPI_INTEGER,rcvInterBeadRanks(1), 
2113:      &        1,MPI_INTEGER,interBeadGroupComm,IERROR) 
2114: ##ELSE 
2115:       call MPI_GROUP_INCL(orig_group, NRPBDS, interBeadGroupRanks, 
2116:      &                 interBeadGroup, ierr) 
2117:       call MPI_COMM_CREATE(MPI_COMM_WORLD, interBeadGroup, 
2118:      &                 interBeadGroupComm, ierr) 
2119:       call MPI_ALLGATHER(WHOIAM,1,MPI_INTEGER,rcvInterBeadRanks(1), 
2120:      &        1,MPI_INTEGER,interBeadGroupComm,IERROR) 
2121: ##ENDIF 
2122:  
2123: ce----PVG error checking 
2124:       write (*,'(A,I3,X,A,X)',advance='no') 'NODE:', WHOIAM, 
2125:      &'rcvinterBeadRanks' 
2126:       do XX = 1, NRPBDS 
2127:         write (*,'(I5)',advance='no') rcvInterBeadRanks(XX) 
2128:       end do 
2129:       write (*,'(X,A)',advance='no') 'rcvinterStateRanks' 
2130:         do XX = 1, nevbs-1 
2131:         write (*,'(I5)',advance='no') rcvInterStateRanks(XX) 
2132:       end do 
2133:       write (*,'(I5)') rcvInterStateRanks(nevbs) 
2134: ce----PVG error checking 
2135:  
2136: ce----PVG error checking 
2137: ce      write (*,'(A,I)') 'interBeadGroupComm =', interBeadGroupComm 
2138: ce     write (*,'(A,I)') 'interStateGroupComm =', interStateGroupComm 
2139: ce----PVG error checking 
2140: ce      write (*,*) 'rank',WHOIAM,'sub evb interBeadGroupComm =', 
2141: ce     &interBeadGroupComm 
2142: ce      write (*,*) 'rank',WHOIAM,'sub evb interStateGroupComm =', 
2143: ce     &interStateGroupComm 
2144:       
2145: ce----PVG END ERROR CHECKING 
2146:  
2147: ce----PVG END ERROR CHECKING 
2148:  
2149:  
2150:       deallocate(interStateGroupRanks) 
2151:       deallocate(rcvInterStateRanks) 
2152:       deallocate(rcvInterBeadRanks) 
2153:       deallocate(interBeadGroupRanks) 
2154:  
2155:  
2156:       RETURN 
2157:       END 
2158:  
2159: C---------------------------------------------------------------------- 
2160:       SUBROUTINE RPMDSETUP(COMLYN,COMLEN) 
2161: C----------------------------------------------------------------------- 
2162: C 
2163: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
2164:       CHARACTER*(*) COMLYN 
2165:       INTEGER COMLEN 
2166: C include global variables from these files: 
2167: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
2168: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
2169: ##INCLUDE '~/charmm_fcm/stream.fcm' 
2170: ##INCLUDE '~/charmm_fcm/psf.fcm' 
2171: ##INCLUDE '~/charmm_fcm/heap.fcm' 
2172: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 
2173: ##INCLUDE '~/charmm_fcm/number.fcm' 
2174: ##INCLUDE '~/charmm_fcm/coord.fcm' 
2175:  
2176:       INTEGER BUFLEN,II,JJ,KK,CTR 
2177:       INTEGER EMSLCT(MAXAIM) 
2178:       CHARACTER*6 WRD 
2179:       REAL*8 R 
2180:       REAL*8 TMPX,TMPY 
2181:  
2182: C     local variable declarations: 
2183:        
2184:       WRITE(*,*)'READING RING POLYMER MOLECULAR DYNAMICS VARIABLES' 
2185:  
2186:       QRPMD=.TRUE. 
2187:  
2188: c------------------------------------ 
2189: c-----PARSING RPMD COMMAND LINE------ 
2190: c------------------------------------ 
2191:  
2192: c-----Get the frequency at which to write to a catonated xyz of the 
2193: c-----full system (all beads). Default is not to print out (0) 
2194:       FULLXYZFRQ=GTRMI(COMLYN,COMLEN,'FULLXYZFRQ',0) 
2195:       IF (PRNLEV .GE. 3) WRITE(OUTU,'(A,I5)')'Frequency of printing out 
2196:      &full system xyz (with all ring polymer beads)=', 
2197:      &FULLXYZFRQ 
2198:       FULLXYZUNIT=GTRMI(COMLYN,COMLEN,'FULLXYZUNIT',0) 
2199:       IF (FULLXYZFRQ .GT. 0 .and. FULLXYZUNIT .LE. 0) THEN 
2200:         CALL WRNDIE(-1,'ENSEMBLE>','IF FULLXYZFRQ IS GREATER THAN 0 
2201:      &FULLXYZUNIT MUST BE SPECIFIED') 
2202:       ENDIF         
2203:       IF (PRNLEV .GE. 3 .and. FULLXYZFRQ .GT. 0) THEN 
2204:         WRITE(OUTU,'(A,I5)')'full system xyz will be printed to unit', 
2205:      &FULLXYZUNIT 
2206:       ENDIF 
2207:        
2208: c-----Get the temperature that defines the RPMD spring force constant 
2209:       FCTEMP=GTRMF(COMLYN,COMLEN,'FCTEMP',-1) 
2210:       IF (PRNLEV .GE. 3) WRITE(OUTU,'(A,F10.3)')'Temperature used to  
2211:      &define ring polymer spring force constants=',FCTEMP 
2212:       IF (FCTEMP.LE.0) THEN 
2213:         CALL WRNDIE(-1,'ENSEMBLE>','FCTEMP MUST BE GREATER THAN 0') 
2214:       ENDIF 
2215:        
2216:       IF ( INDXA(COMLYN, COMLEN, 'FCTEMP') .GT. 0 ) THEN 
2217:         FCTEMP=NEXTF(COMLYN,COMLEN) 
2218:         WRITE(OUTU,'(A,F12.3,A)')'TEMPERATURE USED TO DEFINE RING  
2219:      &POLYMER FORCE CONSTANT = ',FCTEMP,' K' 
2220:         IF (FCTEMP.LE.0) THEN 
2221:         CALL WRNDIE(-1,'ENSEMBLE>','FCTEMP MUST BE GREATER THAN 0') 
2222:         ENDIF 
2223:       ENDIF 
2224:        
2225: c------- read in the atom selection for the atoms to be treated as ring polymers; this is always required 
2226:       EMSLCT=0 
2227:       WRITE(OUTU,'(A)') ' READING ATOM SELECTION TO DEFINE ATOMS TO BE  
2228:      &TREATED AS RING POLYMERS' 
2229:       CALL SELCTA(COMLYN,COMLEN,EMSLCT,X,Y,Z,WMAIN,.TRUE.) 
2230:       NRPMDSEL = 0 
2231:       DO II=1,MAXAIM 
2232:         IF (EMSLCT(II) .EQ. 1) THEN 
2233:           NRPMDSEL = NRPMDSEL + 1 
2234:           RPMDSEL(NRPMDSEL) = II 
2235:         ENDIF 
2236:       ENDDO 
2237:       IF(NRPMDSEL.LT.1) THEN 
2238:         WRITE(OUTU,'(A)') 'RPMD REQUIRES THE SELECTION OF AT LEAST ONE 
2239:      &ATOM' 
2240:         CALL WRNDIE(-1,'ENSEMBLE>','RPMD REQUIRES THE SELECTION OF  
2241:      &AT LEAST ONE ATOM TO BE TREATED AS A RING POLYMER') 
2242:       ELSEIF(NRPMDSEL.GT.10000) THEN 
2243:         WRITE(OUTU,'(A,I5,A)') 'MAX NUMBER OF ALLOWED RPMD ATOMS  
2244:      &(',10000,') EXCEEDED' 
2245:         CALL WRNDIE(-1,'ENSEMBLE>','MAX NUMBER OF ALLOWED RPMD ATOMS 
2246:      &EXCEEDED. This can be increased in ensemble.src') 
2247:       ELSE 
2248:         WRITE(OUTU,'(A)') 'RING POLYMER ATOM SELECTIONS:' 
2249:           DO KK=1,NRPMDSEL 
2250:             WRITE(OUTU,'(A,I5)')' ATOM ',RPMDSEL(KK) 
2251:           ENDDO 
2252:       ENDIF 
2253:        
2254:  
2255: c      IF ( INDXA(COMLYN, COMLEN, 'NRPBDS') .GT. 0 ) THEN 
2256: c        NRPBDS=NEXTI(COMLYN,COMLEN) 
2257: c        IF (NRPBDS.LT.1) THEN 
2258: c          CALL WRNDIE(-1,'ENSEMBLE>','NRPBDS MUST BE GREATER THAN 0') 
2259: c        ENDIF 
2260: c        WRITE(OUTU,'(A,I5)')'NUMBER OF RING POLYMER BEADS:',NRPBDS 
2261: cc        WRITE(OUTU,*)'NUMBER OF RING POLYMER BEADS:',NRPBDS 
2262: c      ENDIF 
2263: c      IF ( INDXA(COMLYN, COMLEN, 'FULLXYZFRQ') .GT. 0 ) THEN 
2264: c        FULLXYZFRQ=NEXTI(COMLYN,COMLEN) 
2265: c        IF (FULLXYZFRQ.LT.0) THEN 
2266: c          CALL WRNDIE(-1,'ENSEMBLE>','FULLXYZFRQ MUST NOT BE  
2267: c     &NEGATIVE') 
2268: c        ENDIF 
2269: c        WRITE(OUTU,'(A,I5)')'FULLXYZFRQ:',FULLXYZFRQ 
2270: c      ENDIF 
2271:  
2272:  
2273:  
2274: c------- reset the elements of EMSLCT to zero 
2275:       EMSLCT=0 
2276:        
2277: c----PROCEDURE: INITIALIZE_BEAD_COORDS 
2278: c----set bead coordinates in a circle so that they are not all on  
2279: c----top of each other 
2280:  
2281: c      R = 0.001 
2282: c      DO JJ=1,NRPMDSEL 
2283: c        CTR=1 
2284: c        IF (MOD(NRPBDS,2).EQ.1) THEN 
2285: c          DO II=1,NRPBDS/2+2 
2286: c            IF(BEADID.EQ.CTR)THEN 
2287: c              TMPX = X(RPMDSEL(JJ)) - R + (II-1)*(4*R/REAL(NRPBDS+1)) 
2288: c              TMPY = SQRT(R*R-(TMPX - X(RPMDSEL(JJ)))* 
2289: c     &        (TMPX - X(RPMDSEL(JJ)))) + Y(RPMDSEL(JJ)) 
2290: c              X(RPMDSEL(JJ)) = TMPX 
2291: c              Y(RPMDSEL(JJ)) = TMPY 
2292: c            ENDIF 
2293: c            CTR = CTR + 1 
2294: c          ENDDO 
2295: c          DO II=NRPBDS/2,2,-1 
2296: c            IF(BEADID.EQ.CTR) THEN 
2297: c              TMPX = X(RPMDSEL(JJ)) - R + (II-1)*(4*R/REAL(NRPBDS+1)) 
2298: c              TMPY = -SQRT(R*R-(TMPX-X(RPMDSEL(JJ)))* 
2299: c     &        (TMPX-X(RPMDSEL(JJ)))) + Y(RPMDSEL(JJ)) 
2300: c              X(RPMDSEL(JJ)) = TMPX 
2301: c              Y(RPMDSEL(JJ)) = TMPY 
2302: c            ENDIF 
2303: c            CTR = CTR + 1 
2304: c          ENDDO 
2305: c        ELSE 
2306: c          DO II=1,NRPBDS/2+1 
2307: c            IF(BEADID.EQ.CTR)THEN 
2308: c              TMPX = X(RPMDSEL(JJ)) - R + (II-1)*(4*R/REAL(NRPBDS)) 
2309: c              TMPY = SQRT(R*R-(TMPX - X(RPMDSEL(JJ)))* 
2310: c     &        (TMPX - X(RPMDSEL(JJ)))) + Y(RPMDSEL(JJ)) 
2311: c              X(RPMDSEL(JJ)) = TMPX 
2312: c              Y(RPMDSEL(JJ)) = TMPY 
2313: c            ENDIF 
2314: c            CTR = CTR + 1 
2315: c          ENDDO 
2316: c          DO II=NRPBDS/2,2,-1 
2317: c            IF(BEADID.EQ.CTR) THEN 
2318: c              TMPX = X(RPMDSEL(JJ)) - R + (II-1)*(4*R/REAL(NRPBDS)) 
2319: c              TMPY = -SQRT(R*R-(TMPX-X(RPMDSEL(JJ)))* 
2320: c     &        (TMPX-X(RPMDSEL(JJ)))) + Y(RPMDSEL(JJ)) 
2321: c              X(RPMDSEL(JJ)) = TMPX 
2322: c              Y(RPMDSEL(JJ)) = TMPY 
2323: c            ENDIF 
2324: c            CTR = CTR + 1 
2325: c          ENDDO 
2326: c        ENDIF 
2327: c      ENDDO 
2328: c-----END PROCEDURE INITIALIZE_BEAD_COORDS 
2329:        
2330:       RETURN 
2331:       END 
2332:  
2333: C---------------------------------------------------------------------- 
2334:       SUBROUTINE RPRXNCRD(COMLYN,COMLEN) 
2335: C----------------------------------------------------------------------- 
2336: C 
2337: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
2338:       CHARACTER*(*) COMLYN 
2339:       INTEGER COMLEN 
2340: C include global variables from these files: 
2341: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
2342: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
2343: ##INCLUDE '~/charmm_fcm/stream.fcm' 
2344: ##INCLUDE '~/charmm_fcm/psf.fcm' 
2345: ##INCLUDE '~/charmm_fcm/heap.fcm' 
2346: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 
2347: ##INCLUDE '~/charmm_fcm/number.fcm' 
2348: ##INCLUDE '~/charmm_fcm/coord.fcm' 
2349:  
2350:       CHARACTER*4 WRD 
2351:       INTEGER EMSLCT(NATOM) 
2352:       INTEGER NSEL, II, KK 
2353:        
2354:  
2355:       QRPRXNCRD = .TRUE. 
2356:        
2357:       RPDISCRD = .FALSE. 
2358:       RPDIFFCRD = .FALSE. 
2359:       
2360:  
2361: c-----Assign unit to write reaction coordinate trace to 
2362:       RPRXNCRDUNIT=GTRMI(COMLYN,COMLEN,'UNIT',0) 
2363:       IF (RPRXNCRDUNIT.LE.0) THEN 
2364:         CALL WRNDIE(-1,'ENSEMBLE>','ASSIGN A UNIT TO WRITE THE REACTION  
2365:      &COORDINATE TRACE TO, FOR EXAMPLE: ENSEMBLE RPRXNCRD UNIT 57 DIS  
2366:      &SELE... END') 
2367:       ENDIF 
2368:  
2369: c-----Set frequency with which to write out the reaction coordinate 
2370: c-----default is every timestep      
2371:       RPRXNFRQ=GTRMI(COMLYN,COMLEN,'FREQ',1) 
2372:  
2373: c-----Determine if the reaction coordinate is distance or 
2374: c-----a difference of two distances 
2375:       WRD=NEXTA4(COMLYN,COMLEN) 
2376:       IF(WRD.EQ.'DIS') THEN 
2377:         RPDISCRD=.TRUE. 
2378:       ELSEIF(WRD.EQ.'DIFF') THEN 
2379:         RPDIFFCRD=.TRUE. 
2380:       ELSE 
2381:         CALL WRNDIE(-1,'ENSEMBLE>','SPECIFY WHETHER THE REACTION  
2382:      &COORDINATE IS A DISTANCE (DIS) OR A DIFFERENCE OF  
2383:      &DISTANCES (DIFF)') 
2384:       ENDIF 
2385:        
2386:       IF(RPDISCRD) THEN 
2387:         IF (IOLEV.GT.0) THEN 
2388:           WRITE(OUTU,'(A)')' A SINGLE DISTANCE WILL BE FOLLOWED AS THE  
2389:      &REACTION COORDINATE. THE CENTROIDS OF ANY RING POLYMER 
2390:      &ATOMS WILL BE USED IN CALCULATING THIS REACTION COORDINATE' 
2391:         ENDIF 
2392:       ELSEIF(RPDIFFCRD) THEN 
2393:         IF (IOLEV.GT.0) THEN 
2394:           WRITE(OUTU,'(A)')' A DIFFERENCE OF DISTANCES R1 - R2 WILL BE 
2395:      &FOLLOWED AS THE REACTION COORDINATE.' 
2396:           WRITE(OUTU,'(A)')'THE CENTROID OF ANY RING POLYMER 
2397:      &ATOMS WILL BE USED IN CALCULATING THIS REACTION COORDINATE' 
2398:         ENDIF 
2399:       ENDIF 
2400:   
2401: c-----read in the first atom selection; this is always required 
2402:       EMSLCT=0 
2403:  
2404:       IF (IOLEV.GT.0) THEN 
2405:         WRITE(OUTU,'(A)') ' READING ATOM SELECTION TO DEFINE R1' 
2406:       ENDIF 
2407:         CALL SELCTA(COMLYN,COMLEN,EMSLCT,X,Y,Z,WMAIN,.TRUE.) 
2408:         NSEL = 0 
2409:          DO II=1,NATOM 
2410:            IF (EMSLCT(II) .EQ. 1) THEN 
2411:              NSEL = NSEL + 1 
2412:              RPRXNSEL(NSEL) = II 
2413:            ENDIF 
2414:          ENDDO 
2415:       IF(NSEL.LT.2.OR.NSEL.GT.2) THEN 
2416:         CALL WRNDIE(-1,'ENSEMBLE>','EACH COMPONENT OF THE REACTION  
2417:      &COORDINATE REQUIRES THE SELECTION OF OF EXACTLY TWO ATOMS') 
2418:       ELSE 
2419:         IF (IOLEV.GT.0) THEN 
2420:           WRITE(OUTU,'(A,I5,A,I5)') 'R1 IS THE DISTANCE BETWEEN  
2421:      &ATOMS',RPRXNSEL(1),' AND ',RPRXNSEL(2) 
2422:         ENDIF 
2423:       ENDIF 
2424:        
2425:        
2426:       IF (RPDIFFCRD) THEN 
2427:         IF (IOLEV.GT.0) THEN 
2428:           WRITE(OUTU,'(A)') ' READING ATOM SELECTION TO DEFINE R2' 
2429:         ENDIF 
2430:           CALL SELCTA(COMLYN,COMLEN,EMSLCT,X,Y,Z,WMAIN,.TRUE.) 
2431:            DO II=1,NATOM 
2432:              IF (EMSLCT(II) .EQ. 1) THEN 
2433:                NSEL = NSEL + 1 
2434:                RPRXNSEL(NSEL) = II 
2435:              ENDIF 
2436:            ENDDO         
2437:         IF(NSEL.LT.4.OR.NSEL.GT.4) THEN 
2438:           CALL WRNDIE(-1,'ENSEMBLE>','EACH COMPONENT OF THE REACTION  
2439:      &COORDINATE REQUIRES THE SELECTION OF OF EXACTLY TWO ATOMS 
2440:      &, IF DIFFCRD IS USED EXACTLY 4 ATOMS MUST BE SPECIFIED IN TOTAL') 
2441:         ELSE 
2442:           IF (IOLEV.GT.0) THEN 
2443:             WRITE(OUTU,'(A,I5,A,I5)') 'R2 IS THE DISTANCE BETWEEN  
2444:      &ATOMS',RPRXNSEL(3),' AND ',RPRXNSEL(4) 
2445:           ENDIF 
2446:         ENDIF 
2447:       ENDIF 
2448:  
2449:  
2450:       RETURN 
2451:       END 
2452:  
2453:  
2454:  
2455: C 
2456: C---------------------------------------------------------------------- 
2457:       SUBROUTINE ENSEXPAVG(QX,QY,QZ,DXK,DYK,DZK,HK,NATOM,AMASS) 
2458: C----------------------------------------------------------------------- 
2459: C the nodes share their forces and energy and and do  
2460: C exp(-H)=exp(-HA)+exp(-HB) averaging 
2461: C 
2462: C ... calculate energy & forces at each time step 
2463: C 
2464: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
2465: C 
2466:       REAL*8 QX(*),QY(*),QZ(*),AMASS(*) 
2467:       REAL*8 DXK(*),DYK(*),DZK(*),HK 
2468:       INTEGER NATOM  
2469:    
2470: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
2471: ##INCLUDE '~/charmm_fcm/stream.fcm' 
2472: ##INCLUDE '~/charmm_fcm/heap.fcm' 
2473: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 
2474: C local: 
2475:       INTEGER BUFLEN 
2476: C 
2477:       BUFLEN = NATOM*NENSEM 
2478:       IF(QENSEXP) THEN 
2479:         CALL ENSEXPAVG2(DXK,DYK,DZK,HK,NATOM,HEAP(ENSDX), 
2480:      1                  HEAP(ENSDY),HEAP(ENSDZ),HEAP(ENSH),BUFLEN) 
2481:       ELSE 
2482:         IF(QEVB) THEN 
2483:           CALL EVB(QX,QY,QZ,DXK,DYK,DZK,HK,NATOM,BUFLEN) 
2484:         ENDIF 
2485:         IF(QRPMD) THEN 
2486:           CALL RPMD(QX,QY,QZ,DXK,DYK,DZK,HK,BUFLEN) 
2487:         ENDIF 
2488:       ENDIF 
2489:       RETURN 
2490:       END 
2491:   
2492:  
2493:  
2494:       
2495:       SUBROUTINE RPMD(QX,QY,QZ,DXK,DYK,DZK,HK,BUFLEN) 
2496: C----------------------------------------------------------------------- 
2497: C The nodes share their forces and energy; the total energy is calculated 
2498: c as the smallest eigenvalue of an NENSEM by NENSEM matrix of mm energies. 
2499: c Gradients are calculated on the EVB surface using Feynmann-Hellman;  
2500: c n-state evb functionality coded by david glowacki, sept 2010 
2501: C 
2502: C 
2503: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
2504: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
2505: ##INCLUDE '~/charmm_fcm/psf.fcm' 
2506: ##INCLUDE '~/charmm_fcm/stream.fcm' 
2507: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 
2508: ##INCLUDE '~/charmm_fcm/contrl.fcm' 
2509: ##INCLUDE '~/charmm_fcm/number.fcm' 
2510: ##INCLUDE '~/charmm_fcm/consta.fcm' 
2511: ##INCLUDE '~/charmm_fcm/energy.fcm' 
2512:       INCLUDE 'mpif.h' 
2513:  
2514: C     arguments: 
2515: C     qx,qy,qz are the coordinates 
2516:       REAL*8 QX(*),QY(*),QZ(*) 
2517:       REAL*8 DXK(*),DYK(*),DZK(*),RCVDX(NATOM*NRPBDS) 
2518:       REAL*8 RCVDY(NATOM*NRPBDS),RCVDZ(NATOM*NRPBDS) 
2519:       REAL*8 RCVH(NRPBDS),HK 
2520: C     local: 
2521:       REAL*8 TMPH(1) 
2522:       integer tempprnlev 
2523:       INTEGER BUFLEN,IERROR 
2524:       LOGICAL QAVFORCE 
2525:       LOGICAL SCRATCHLOGICAL 
2526:       REAL*8 FORCESUM, FORCEAVE 
2527:       REAL*8 COORSUM, COORAVE 
2528:       INTEGER II, JJ, KK, CTR 
2529:       REAL*8 BETA, BETAN, OMEGA, OMEGAN 
2530:       REAL*8 DX1, DX2, DX3, DY1, DY2, DY3, DZ1, DZ2, DZ3 
2531:       REAL*8 RPMDX(NRPMDSEL),RPMDY(NRPMDSEL),RPMDZ(NRPMDSEL) 
2532:       REAL*8 RCVRPX(NRPMDSEL*NRPBDS) 
2533:       REAL*8 RCVRPY(NRPMDSEL*NRPBDS) 
2534:       REAL*8 RCVRPZ(NRPMDSEL*NRPBDS) 
2535:       REAL*8 CENTROIDX(NATOM) 
2536:       REAL*8 CENTROIDY(NATOM) 
2537:       REAL*8 CENTROIDZ(NATOM) 
2538:       REAL*8,parameter :: aukboltz = 3.166829d-6 
2539:       REAL*8, parameter :: auhbar   = 1.d0 
2540:       REAL*8, parameter :: au_to_kcalmol = 627.509646392d0 
2541:       REAL*8, parameter :: kcalmol_to_au = 0.001593601d0 
2542:       REAL*8, parameter :: ang_to_bohr   = 1.889726128d0 
2543:       REAL*8, parameter :: bohr_to_ang   = 0.5291772108d0 
2544:       REAL*8, parameter :: atomicmass_to_me = 1822.88839d0 
2545:       character*4 :: basename 
2546:       character*5 :: basename2 
2547:       character*4 :: extension 
2548:       character(100):: filename 
2549:       character*1 :: kstring 
2550:       character*1 :: kstring2 
2551:       REAL*8 TMPQX(NATOM), TMPQY(NATOM), TMPQZ(NATOM) 
2552:       REAL*8 TMPQX2(NATOM), TMPQY2(NATOM), TMPQZ2(NATOM) 
2553:       REAL*8 DEBUGX,DEBUGY,DEBUGZ 
2554:       REAL*8 RCDIS1, RCDIS2 
2555:       REAL*8 SPRINGEN 
2556:                    
2557: ce     CALL MPI_BARRIER(interBeadGroupComm, IERROR) 
2558: ce     call MPI_ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION, 
2559: ce    &RCVDX(1),NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2560: ce     &rcvArray1(1),NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2561: ce     call MPI_BARRIER(MPI_COMM_WORLD,IERROR) 
2562: ce      write(*,*)rcvArray1 
2563:  
2564:        
2565:       DO II=1,NRPMDSEL 
2566:         RPMDX(II)=QX(RPMDSEL(II)) 
2567:         RPMDY(II)=QY(RPMDSEL(II)) 
2568:         RPMDZ(II)=QZ(RPMDSEL(II)) 
2569:       ENDDO 
2570:       DO II=1,NATOM 
2571:         TMPQX(II) = QX(II) 
2572:         TMPQY(II) = QY(II) 
2573:         TMPQZ(II) = QZ(II) 
2574:         TMPQX2(II) = QX(II) 
2575:         TMPQY2(II) = QY(II) 
2576:         TMPQZ2(II) = QZ(II) 
2577:       ENDDO 
2578:       TMPH(1) = HK 
2579: ##IF ALTIX 
2580:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR) 
2581:       CALL MPI8ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1), 
2582:      1       NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2583:       CALL MPI8ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1), 
2584:      1       NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2585:       CALL MPI8ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1), 
2586:      1       NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2587:       CALL MPI8ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1), 
2588:      1       1,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2589:       CALL MPI_ALLGATHER(RPMDX(1),NRPMDSEL,MPI_DOUBLE_PRECISION, 
2590:      1RCVRPX(1),NRPMDSEL,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2591:       CALL MPI_ALLGATHER(RPMDY(1),NRPMDSEL,MPI_DOUBLE_PRECISION, 
2592:      1RCVRPY(1),NRPMDSEL,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2593:       CALL MPI_ALLGATHER(RPMDZ(1),NRPMDSEL,MPI_DOUBLE_PRECISION, 
2594:      1RCVRPZ(1),NRPMDSEL,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2595:       CALL MPI8BCAST(TMPQX(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2596:      1MPI_COMM_WORLD,IERROR) 
2597:       CALL MPI8BCAST(TMPQY(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2598:      1MPI_COMM_WORLD,IERROR) 
2599:       CALL MPI8BCAST(TMPQZ(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2600:      1MPI_COMM_WORLD,IERROR) 
2601:       CALL MPI8BCAST(TMPQX2(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2602:      1interStateGroupComm,IERROR) 
2603:       CALL MPI8BCAST(TMPQY2(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2604:      1interStateGroupComm,IERROR) 
2605:       CALL MPI8BCAST(TMPQZ2(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2606:      1interStateGroupComm,IERROR) 
2607: ##ELSE 
2608:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR) 
2609:       CALL MPI_ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1), 
2610:      1       NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2611:       CALL MPI_ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1), 
2612:      1       NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2613:       CALL MPI_ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1), 
2614:      1       NATOM,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2615:       CALL MPI_ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1), 
2616:      1       1,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2617:       CALL MPI_ALLGATHER(RPMDX(1),NRPMDSEL,MPI_DOUBLE_PRECISION, 
2618:      1RCVRPX(1),NRPMDSEL,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2619:       CALL MPI_ALLGATHER(RPMDY(1),NRPMDSEL,MPI_DOUBLE_PRECISION, 
2620:      1RCVRPY(1),NRPMDSEL,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2621:       CALL MPI_ALLGATHER(RPMDZ(1),NRPMDSEL,MPI_DOUBLE_PRECISION, 
2622:      1RCVRPZ(1),NRPMDSEL,MPI_DOUBLE_PRECISION,interBeadGroupComm,IERROR) 
2623:       CALL MPI_BCAST(TMPQX(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2624:      1MPI_COMM_WORLD,IERROR) 
2625:       CALL MPI_BCAST(TMPQY(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2626:      1MPI_COMM_WORLD,IERROR) 
2627:       CALL MPI_BCAST(TMPQZ(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2628:      1MPI_COMM_WORLD,IERROR) 
2629:       CALL MPI_BCAST(TMPQX2(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2630:      1interStateGroupComm,IERROR) 
2631:       CALL MPI_BCAST(TMPQY2(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2632:      1interStateGroupComm,IERROR) 
2633:       CALL MPI_BCAST(TMPQZ2(1),NATOM,MPI_DOUBLE_PRECISION,0, 
2634:      1interStateGroupComm,IERROR) 
2635: ##ENDIF 
2636: ce      IF (IOLEV.GE.0) THEN 
2637: ce      ENDIF 
2638: ce      write(*,'(A,I,A,36F10.3)')'node',WHOIAM,'recvddxk',RCVDX 
2639: ce      write(*,'(A,I,A,36F10.3)')'node',WHOIAM,'recvddyk',RCVDY 
2640: ce      write(*,'(A,I,A,36F10.3)')'node',WHOIAM,'recvddzk',RCVDZ 
2641:  
2642: c-----for non quantum atoms, average the forces from the different beads 
2643:  
2644:  
2645:    
2646:       DO II=1,NATOM 
2647:         QAVFORCE = .True. 
2648:         DO JJ=1,NRPMDSEL 
2649:           IF(II.EQ.RPMDSEL(JJ)) THEN 
2650:             QAVFORCE = .False. 
2651:           ENDIF 
2652:         ENDDO 
2653:         IF(QAVFORCE) THEN 
2654:           FORCESUM = 0 
2655:           DO KK=1,NRPBDS  
2656:             FORCESUM = FORCESUM + RCVDX(II + (KK-1)*NATOM) 
2657:           ENDDO 
2658:           FORCEAVE = FORCESUM/dble(NRPBDS) 
2659:           DXK(II) = FORCEAVE 
2660: c---------set non rpmd atom position to that of the equivalent atom on bead 1 
2661:           QX(II) = TMPQX(II) 
2662:           CENTROIDX(II)=QX(II) 
2663:         ELSE 
2664: c---------for rpmd atoms average position to get centroid position 
2665:           COORSUM = 0 
2666:           COORAVE = 0 
2667:           DO JJ=1,NRPMDSEL 
2668:             IF(RPMDSEL(JJ).EQ.II) THEN 
2669:               DO KK=1,NRPBDS  
2670:                 COORSUM = COORSUM + RCVRPX(JJ + NRPMDSEL*(KK-1)) 
2671:               ENDDO 
2672:             ENDIF 
2673:           ENDDO 
2674:           COORAVE = COORSUM/dble(NRPBDS) 
2675:           CENTROIDX(II) = COORAVE 
2676:         ENDIF 
2677:       ENDDO 
2678:        
2679:  
2680:  
2681:       DO II=1,NATOM 
2682:         QAVFORCE = .True. 
2683:         DO JJ=1,NRPMDSEL 
2684:           IF(II.EQ.RPMDSEL(JJ)) THEN 
2685:             QAVFORCE = .False. 
2686:           ENDIF 
2687:         ENDDO 
2688:         IF(QAVFORCE) THEN 
2689:           FORCESUM = 0 
2690:           DO KK=1,NRPBDS  
2691:             FORCESUM = FORCESUM + RCVDY(II + (KK-1)*NATOM) 
2692:           ENDDO 
2693:           FORCEAVE = FORCESUM/dble(NRPBDS) 
2694:           DYK(II) = FORCEAVE 
2695: c---------set non rpmd atom position to that of the equivalent atom on bead 1 
2696:           QY(II) = TMPQY(II) 
2697:           CENTROIDY(II) = QY(II) 
2698:         ELSE 
2699: c---------for rpmd atoms average position to get centroid position 
2700:           COORSUM = 0 
2701:           COORAVE = 0 
2702:           DO JJ=1,NRPMDSEL 
2703:             IF(RPMDSEL(JJ).EQ.II) THEN 
2704:               DO KK=1,NRPBDS  
2705:                 COORSUM = COORSUM + RCVRPY(JJ + NRPMDSEL*(KK-1)) 
2706:               ENDDO 
2707:             ENDIF 
2708:           ENDDO 
2709:           COORAVE = COORSUM/dble(NRPBDS) 
2710:           CENTROIDY(II) = COORAVE 
2711:         ENDIF 
2712:       ENDDO 
2713:  
2714:       DO II=1,NATOM 
2715:         QAVFORCE = .True. 
2716:         DO JJ=1,NRPMDSEL 
2717:           IF(II.EQ.RPMDSEL(JJ)) THEN 
2718:             QAVFORCE = .False. 
2719:           ENDIF 
2720:         ENDDO 
2721:         IF(QAVFORCE) THEN 
2722:           FORCESUM = 0 
2723:           DO KK=1,NRPBDS  
2724:             FORCESUM = FORCESUM + RCVDZ(II + (KK-1)*NATOM) 
2725:           ENDDO 
2726:           FORCEAVE = FORCESUM/dble(NRPBDS) 
2727:           DZK(II) = FORCEAVE 
2728: c---------set non rpmd atom position to that of the equivalent atom on bead 1 
2729:           QZ(II) = TMPQZ(II) 
2730:           CENTROIDZ(II) = QZ(II) 
2731:         ELSE 
2732: c---------for rpmd atoms average position to get centroid position 
2733:           COORSUM = 0 
2734:           COORAVE = 0 
2735:           DO JJ=1,NRPMDSEL 
2736:             IF(RPMDSEL(JJ).EQ.II) THEN 
2737:               DO KK=1,NRPBDS  
2738:                 COORSUM = COORSUM + RCVRPZ(JJ + NRPMDSEL*(KK-1)) 
2739:               ENDDO 
2740:             ENDIF 
2741:           ENDDO 
2742:           COORAVE = COORSUM/dble(NRPBDS) 
2743:           CENTROIDZ(II) = COORAVE 
2744:         ENDIF 
2745:       ENDDO 
2746:              
2747: c---------make sure each replica in the same beadgroup have the same coordinates 
2748:           DO JJ=1,NRPMDSEL 
2749:             QX(RPMDSEL(JJ)) = TMPQX2(RPMDSEL(JJ))  
2750:             QY(RPMDSEL(JJ)) = TMPQY2(RPMDSEL(JJ))  
2751:             QZ(RPMDSEL(JJ)) = TMPQZ2(RPMDSEL(JJ))  
2752:           ENDDO 
2753:  
2754:  
2755: cc--------write out centroid coordinates 
2756: cc      IF(MOD(MDSTEP,1000).EQ.0) THEN 
2757:       IF(MOD(MDSTEP,1000).EQ.0) THEN 
2758:         IF(WHOIAM.EQ.0) THEN 
2759:           open(60,file = 'centroids.xyz', 
2760:      &status='unknown', access='append') 
2761:           write(60,*),NATOM 
2762:           write(60,*),'rpmd centroids','timestep=',MDSTEP 
2763:           DO II=1,NATOM 
2764:             write(60,'(A,F10.6,F10.6,F10.6)')TYPE(II),CENTROIDX(II), 
2765:      &CENTROIDY(II),CENTROIDZ(II) 
2766:           ENDDO 
2767:           close(60) 
2768:         ENDIF 
2769:       ENDIF 
2770:  
2771: c----------------------------------------------------------------- 
2772: c--------write out all coordinates including all RPMD beads------- 
2773: c----------------------------------------------------------------- 
2774:       IF(FULLXYZFRQ .GT. 0 .and. MOD(MDSTEP,FULLXYZFRQ).EQ.0) THEN 
2775:         IF(WHOIAM.EQ.0) THEN 
2776:           write(FULLXYZUNIT,*),NATOM+NRPMDSEL*(NRPBDS-1) 
2777:           write(FULLXYZUNIT,*),'rpmd all atoms','timestep=',MDSTEP 
2778:           DO JJ=1,NRPMDSEL 
2779:             DO II=1,NRPBDS 
2780:               write(FULLXYZUNIT,'(A,F10.6,F10.6,F10.6)') 
2781:      &TYPE(RPMDSEL(JJ)),RCVRPX(JJ+(II-1)*NRPMDSEL), 
2782:      &RCVRPY(JJ+(II-1)*NRPMDSEL),RCVRPZ(JJ+(II-1)*NRPMDSEL) 
2783:             ENDDO 
2784:           ENDDO 
2785:           DO II=1,NATOM 
2786:             SCRATCHLOGICAL = .TRUE. 
2787:             DO JJ=1,NRPMDSEL 
2788:               IF (RPMDSEL(JJ).EQ.II) SCRATCHLOGICAL = .FALSE. 
2789:             ENDDO 
2790:             IF (SCRATCHLOGICAL .EQV. .TRUE.) THEN 
2791:               write(FULLXYZUNIT,'(A,F10.6,F10.6,F10.6)')TYPE(II), 
2792:      &TMPQX(II),TMPQY(II),TMPQZ(II) 
2793:             ENDIF 
2794:           ENDDO 
2795:         ENDIF 
2796:       ENDIF 
2797:  
2798: cc--------write out reaction coordinate (specific to malonaldehyde 
2799: c      IF(WHOIAM.EQ.0) THEN 
2800: c        open(61,file = 'trace.dat',status='unknown',access='append') 
2801: c        RCDIS1=SQRT((CENTROIDX(1)-CENTROIDX(2))**2 
2802: c     &               +(CENTROIDY(1)-CENTROIDY(2))**2 
2803: c     &               +(CENTROIDZ(1)-CENTROIDZ(2))**2) 
2804: c        RCDIS2=SQRT((CENTROIDX(1)-CENTROIDX(5))**2 
2805: c     &               +(CENTROIDY(1)-CENTROIDY(5))**2 
2806: c     &               +(CENTROIDZ(1)-TMPQZ(5))**2) 
2807: c        write(61,*)MDSTEP,RCDIS1-RCDIS2 
2808: c        close(61) 
2809: c      ENDIF     
2810: c-------------------------------------------- 
2811: c--------write out reaction coordinate------- 
2812: c-------------------------------------------- 
2813:       IF(QRPRXNCRD .and. MOD(MDSTEP,RPRXNFRQ).EQ.0) THEN 
2814:         IF(WHOIAM.EQ.0) THEN 
2815:           IF(RPDISCRD) THEN 
2816:             RCDIS1=SQRT((CENTROIDX(RPRXNSEL(1))- 
2817:      &             CENTROIDX(RPRXNSEL(2)))**2 
2818:      &             +(CENTROIDY(RPRXNSEL(1))-CENTROIDY(RPRXNSEL(2)))**2 
2819:      &             +(CENTROIDZ(RPRXNSEL(1))-CENTROIDZ(RPRXNSEL(2)))**2) 
2820:             WRITE(RPRXNCRDUNIT,'(I10,F10.5)')MDSTEP,RCDIS1 
2821:           ELSEIF(RPDIFFCRD) THEN 
2822:             RCDIS1=SQRT((CENTROIDX(RPRXNSEL(1))- 
2823:      &             CENTROIDX(RPRXNSEL(2)))**2 
2824:      &             +(CENTROIDY(RPRXNSEL(1))-CENTROIDY(RPRXNSEL(2)))**2 
2825:      &             +(CENTROIDZ(RPRXNSEL(1))-CENTROIDZ(RPRXNSEL(2)))**2) 
2826:             RCDIS2=SQRT((CENTROIDX(RPRXNSEL(3)) 
2827:      &             -CENTROIDX(RPRXNSEL(4)))**2 
2828:      &             +(CENTROIDY(RPRXNSEL(3))-CENTROIDY(RPRXNSEL(4)))**2 
2829:      &             +(CENTROIDZ(RPRXNSEL(3))-CENTROIDZ(RPRXNSEL(4)))**2) 
2830:             WRITE(RPRXNCRDUNIT,'(I10,F10.5)')MDSTEP,RCDIS1-RCDIS2 
2831:           ENDIF 
2832:         ENDIF 
2833:       ENDIF     
2834:          
2835:  
2836:  
2837:  
2838: c-----calculate the forces due to the spring terms for the  
2839: c-----RPMD selected atoms 
2840:  
2841:       beta  = 1.d0/(aukboltz*FCTEMP) 
2842:       betan = beta/dble(NRPBDS) 
2843:       omega = 1.d0/ (betan * auhbar) 
2844:       omegan = omega*omega 
2845:       SPRINGEN = 0 
2846:       DO II=1,NRPMDSEL 
2847:         IF(BEADID.EQ.1) THEN 
2848:           DX1 = RCVRPX(II) 
2849:           DX2 = RCVRPX(II + NRPMDSEL*(NRPBDS-1)) 
2850:           DX3 = RCVRPX(II + NRPMDSEL) 
2851:           DY1 = RCVRPY(II) 
2852:           DY2 = RCVRPY(II + NRPMDSEL*(NRPBDS-1)) 
2853:           DY3 = RCVRPY(II + NRPMDSEL) 
2854:           DZ1 = RCVRPZ(II) 
2855:           DZ2 = RCVRPZ(II + NRPMDSEL*(NRPBDS-1)) 
2856:           DZ3 = RCVRPZ(II + NRPMDSEL) 
2857:         ELSEIF(BEADID.EQ.NRPBDS) THEN 
2858:           DX1 = RCVRPX(II + NRPMDSEL*(NRPBDS-1)) 
2859:           DX2 = RCVRPX(II + NRPMDSEL*(BEADID-2)) 
2860:           DX3 = RCVRPX(II) 
2861:           DY1 = RCVRPY(II + NRPMDSEL*(NRPBDS-1)) 
2862:           DY2 = RCVRPY(II + NRPMDSEL*(BEADID-2)) 
2863:           DY3 = RCVRPY(II) 
2864:           DZ1 = RCVRPZ(II + NRPMDSEL*(NRPBDS-1)) 
2865:           DZ2 = RCVRPZ(II + NRPMDSEL*(BEADID-2)) 
2866:           DZ3 = RCVRPZ(II) 
2867:         ELSE 
2868:           DX1 = RCVRPX(II + NRPMDSEL*(BEADID-1)) 
2869:           DX2 = RCVRPX(II + NRPMDSEL*(BEADID-2)) 
2870:           DX3 = RCVRPX(II + NRPMDSEL*(BEADID)) 
2871:           DY1 = RCVRPY(II + NRPMDSEL*(BEADID-1)) 
2872:           DY2 = RCVRPY(II + NRPMDSEL*(BEADID-2)) 
2873:           DY3 = RCVRPY(II + NRPMDSEL*(BEADID)) 
2874:           DZ1 = RCVRPZ(II + NRPMDSEL*(BEADID-1)) 
2875:           DZ2 = RCVRPZ(II + NRPMDSEL*(BEADID-2)) 
2876:           DZ3 = RCVRPZ(II + NRPMDSEL*(BEADID))           
2877:         ENDIF 
2878: ce------debugging 
2879:         DEBUGX = DXK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang 
2880:         DEBUGY = DYK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang 
2881:         DEBUGZ = DZK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang 
2882: ce------debugging         
2883:           DXK(RPMDSEL(II)) = DXK(RPMDSEL(II)) + 
2884:      &(AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DX1-DX2-DX3)* 
2885:      &ang_to_bohr)*ang_to_bohr*au_to_kcalmol 
2886:           DYK(RPMDSEL(II)) = DYK(RPMDSEL(II)) + 
2887:      &(AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DY1-DY2-DY3)* 
2888:      &ang_to_bohr)*ang_to_bohr*au_to_kcalmol 
2889:           DZK(RPMDSEL(II)) = DZK(RPMDSEL(II)) + 
2890:      &(AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DZ1-DZ2-DZ3)* 
2891:      &ang_to_bohr)*ang_to_bohr*au_to_kcalmol 
2892:  
2893: ce-------PVG debugging 
2894:       IF(WHOIAM.EQ.0) THEN 
2895:          open(59, file = "debugging.txt",status='unknown',  
2896:      &access='append') 
2897:          write(59,*)AMASS(RPMDSEL(II)) 
2898:          close(59) 
2899:       ENDIF 
2900: ce------PVG end debugging  
2901:      
2902:  
2903:           SPRINGEN = SPRINGEN + 0.5d0 *  
2904:      &(AMASS(RPMDSEL(II))*atomicmass_to_me * omegan * 
2905:      &((DX1-DX2)*(DX1-DX2)+(DY1-DY2)*(DY1-DY2)+(DZ1-DZ2)*(DZ1-DZ2))* 
2906:      &ang_to_bohr*ang_to_bohr)*au_to_kcalmol 
2907:           basename = 'springen' 
2908:           extension = '.dat' 
2909:           write( kstring, '(i1)' ) BEADID 
2910:           filename = basename//kstring//extension 
2911:  
2912:            
2913:        ENDDO 
2914: c       IF(STATEID.EQ.1)THEN 
2915: c         IF(MOD(MDSTEP,1000).EQ.0) THEN 
2916: c         open(59, file = filename,status='unknown', access='append') 
2917: c         write(59,'(F20.5$)') SPRINGEN 
2918: c         write(59,'(A)') ' ' 
2919: c         close(59)    
2920: c         ENDIF 
2921: c       ENDIF 
2922:         
2923: c       IF(STATEID.EQ.1) THEN 
2924: c         IF(MOD(MDSTEP,1000).EQ.0) THEN 
2925: c           tempprnlev = PRNLEV 
2926: c           PRNLEV = 5 
2927: c           basename = 'beadenergy' 
2928: c           extension = '.log' 
2929: c           write( kstring, '(i1)' ) BEADID 
2930: c           filename = basename//kstring//extension 
2931: c           OPEN(59, file = filename,status='unknown', 
2932: c     &access='append') 
2933: c           CALL PRINTE(59, EPROP, ETERM,'DYNA','DYN', 
2934: c     &                .TRUE., ZERO, ZERO, ZERO, .TRUE.) 
2935: c           CLOSE(59) 
2936: c           PRNLEV = tempprnlev 
2937: c         ENDIF 
2938: c       ENDIF 
2939:            
2940:  
2941: ce-------PVG DEBUGGING 
2942: c          basename = 'debuggin' 
2943: c          extension = '.xyz' 
2944: c          write( kstring, '(i1)' ) BEADID  
2945: c          write( kstring2, '(i1)' ) STATEID 
2946: c          filename = basename//kstring//kstring2//extension 
2947: c      open(59, file = filename,status='unknown') 
2948: c          DO CTR=1,SIZE(RCVRPX) 
2949: c            write(59,*),'RCVRPX',RCVRPX(CTR)*ang_to_bohr 
2950: c          ENDDO     
2951: c          write(59,*)'mass',AMASS(RPMDSEL(II))*atomicmass_to_me 
2952: c          write(59,*)'omegan',omegan 
2953: c          write(59,*)'2.d0*DZ1-DZ2-DZ3',ang_to_bohr*(2.d0*DZ1-DZ2-DZ3) 
2954: c          write(59,*)'beadid',BEADID 
2955: c          write(59,*)'dz1',DZ1*ang_to_bohr 
2956: c          write(59,*)'dz2',DZ2*ang_to_bohr 
2957: c          write(59,*)'dz3',DZ3*ang_to_bohr 
2958: c          write(59,*)'2.d0*DZ1-DZ2-DZ3',(2.d0*DZ1-DZ2-DZ3)*ang_to_bohr 
2959: c          write(59,*)'dvdr before springs',DEBUGZ 
2960: c          write(59,*)'dvdr due to springs',AMASS(RPMDSEL(II))* 
2961: c     &atomicmass_to_me*omegan*(2.d0*DZ1-DZ2-DZ3)*ang_to_bohr 
2962: c          write(59,*)'dvdr after springs',AMASS(RPMDSEL(II))* 
2963: c     &atomicmass_to_me*omegan*(2.d0*DZ1-DZ2-DZ3)*ang_to_bohr + 
2964: c     &DEBUGZ 
2965: c          write(59,*)'force in charmm units',DZK(RPMDSEL(II)) 
2966: c          write(59,*)DZK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang 
2967: c          write(59,*)DZK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang + 
2968: c     &AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DZ1-DZ2-DZ3)* 
2969: c     &ang_to_bohr 
2970: c          write(59,*)(DZK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang + 
2971: c     &AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DZ1-DZ2-DZ3)* 
2972: c     &ang_to_bohr)*ang_to_bohr*au_to_kcalmol 
2973: c         write(59,*)'mass',AMASS(RPMDSEL(II))*atomicmass_to_me 
2974: c         write(59,*)'omegan',omegan 
2975: c         write(59,*)'2.d0*DY1-DY2-DY3',ang_to_bohr*(2.d0*DY1-DY2-DY3) 
2976: c         write(59,*)'beadid',BEADID 
2977: c         write(59,*)'dy1',DY1*ang_to_bohr 
2978: c         write(59,*)'dy2',DY2*ang_to_bohr 
2979: c         write(59,*)'dy3',DY3*ang_to_bohr 
2980: c         write(59,*)'2.d0*DY1-DY2-DY3',(2.d0*DY1-DY2-DY3)*ang_to_bohr 
2981: c         write(59,*)'dvdr before springs',DEBUGY 
2982: c         write(59,*)'dvdr due to springs',AMASS(RPMDSEL(II))* 
2983: c     &atomicmass_to_me*omegan*(2.d0*DY1-DY2-DY3)*ang_to_bohr 
2984: c          write(59,*)'dvdr after springs',AMASS(RPMDSEL(II))* 
2985: c     &atomicmass_to_me*omegan*(2.d0*DY1-DY2-DY3)*ang_to_bohr + 
2986: c     &DEBUGY 
2987: c         write(59,*)'force in charmm units',DYK(RPMDSEL(II)) 
2988: c         write(59,*)DYK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang 
2989: c         write(59,*)DYK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang + 
2990: c     &AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DY1-DY2-DY3)* 
2991: c     &ang_to_bohr 
2992: c          write(59,*)(DYK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang + 
2993: c     &AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DY1-DY2-DY3)* 
2994: c     &ang_to_bohr)*ang_to_bohr*au_to_kcalmol 
2995: c          write(59,*)'mass',AMASS(RPMDSEL(II))*atomicmass_to_me 
2996: c          write(59,*)'omegan',omegan 
2997: c          write(59,*)'hbar',auhbar    
2998: c          write(59,*)'betan',betan   
2999: c          write(59,*)'2.d0*DX1-DX2-DX3',ang_to_bohr*(2.d0*DX1-DX2-DX3) 
3000: c          write(59,*)'beadid',BEADID 
3001: c          write(59,*)'dx1',DX1*ang_to_bohr 
3002: c          write(59,*)'dx2',DX2*ang_to_bohr 
3003: c          write(59,*)'dx3',DX3*ang_to_bohr 
3004: c          write(59,*)'2.d0*DX1-DX2-DX3',(2.d0*DX1-DX2-DX3)*ang_to_bohr 
3005: c          write(59,*)'dvdr before springs',DEBUGX 
3006: c          write(59,*)'dvdr due to springs',AMASS(RPMDSEL(II))* 
3007: c     &atomicmass_to_me*omegan*(2.d0*DX1-DX2-DX3)*ang_to_bohr 
3008: c          write(59,*)'dvdr after springs',AMASS(RPMDSEL(II))* 
3009: c     &atomicmass_to_me*omegan*(2.d0*DX1-DX2-DX3)*ang_to_bohr + 
3010: c     &DEBUGX 
3011: c          write(59,*)'force in charmm units',DXK(RPMDSEL(II)) 
3012: c          write(59,*)DXK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang 
3013: c          write(59,*)DXK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang + 
3014: c     &AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DX1-DX2-DX3)* 
3015: c     &ang_to_bohr 
3016: c          write(59,*)(DXK(RPMDSEL(II))*kcalmol_to_au*bohr_to_ang + 
3017: c     &AMASS(RPMDSEL(II))*atomicmass_to_me*omegan*(2.d0*DX1-DX2-DX3)* 
3018: c     &ang_to_bohr)*ang_to_bohr*au_to_kcalmol 
3019: c          close(59) 
3020: ce        IF(beadid.EQ.4) THEN 
3021: ce          open(59,file = 'debuggin.txt',status='unknown') 
3022: ce          write(59,*)'beadid',beadid,'DYK',DYK(RPMDSEL(II)) 
3023: ce          write(59,*)'dx1',DX1*ang_to_bohr 
3024: ce          close(59)     
3025: ce        ENDIF  
3026:             
3027:  
3028: ce----PVG ERROR CHECKING PRINT OUT GEOMETRIES OF EACH REPLICA 
3029: c      basename = 'bead' 
3030: c      basename2 = 'state' 
3031: c      extension = '.xyz' 
3032: c      write( kstring, '(i1)' )  BEADID 
3033: c      write( kstring2, '(i1)' ) STATEID 
3034: cc      write(*,'(A)') kstring 
3035: c      filename = basename//kstring//basename2//kstring2//extension 
3036: c      open(61, file = filename,status='unknown') 
3037: c      write(61,'(i5)')NATOM 
3038: c      write(61,'("* title")') 
3039: c      do II = 1, NATOM 
3040: c        write(61,'(a10,F10.5,F10.5,F10.5)')"H",QX(II),QY(II),QZ(II) 
3041: c      enddo 
3042: c      close(61) 
3043: c 
3044: cce----PVG ERROR CHECKING PRINT OUT FORCE FOR EACH REPLICA 
3045: c      basename = 'bead' 
3046: c      basename2 = 'state' 
3047: c      extension = '.force' 
3048: c      write( kstring, '(i1)' )  BEADID 
3049: c      write( kstring2, '(i1)' ) STATEID 
3050: cc      write(*,'(A)') kstring 
3051: c      filename = basename//kstring//basename2//kstring2//extension 
3052: c      open(61, file = filename,status='unknown') 
3053: c      write(61,'(i5)')NATOM 
3054: c      write(61,'("* title")') 
3055: c      do II = 1, NATOM 
3056: c        write(61,'(a10,F10.5,F10.5,F10.5)')"H",DXK(II),DYK(II),DZK(II) 
3057: c      enddo 
3058: c      close(61) 
3059:  
3060: ce-------PVG DEBUGGING END       
3061:  
3062:  
3063:            
3064:       RETURN 
3065:       END 
3066:  
3067:        
3068: C 
3069:       SUBROUTINE ENSEXPAVG2(DXK,DYK,DZK,HK,NATOM, 
3070:      1                      RCVDX,RCVDY,RCVDZ,RCVH,BUFLEN) 
3071: C----------------------------------------------------------------------- 
3072: C the nodes share their forces and energy and and do exp(-H)=exp(-HA)+exp(-HB) 
3073: C averaging 
3074: C 
3075: C ... calculate energy & forces at each time step 
3076: C 
3077: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
3078: ##INCLUDE '~/charmm_fcm/stream.fcm' 
3079: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 
3080: ##INCLUDE '~/charmm_fcm/contrl.fcm' 
3081: ##INCLUDE '~/charmm_fcm/number.fcm' 
3082:       INCLUDE 'mpif.h' 
3083:  
3084: C     arguments: 
3085:       REAL*8 DXK(*),DYK(*),DZK(*),RCVDX(*),RCVDY(*),RCVDZ(*) 
3086:       REAL*8 RCVH(*),HK 
3087:       INTEGER NATOM,BUFLEN 
3088:  
3089: C     local: 
3090:       INTEGER IERROR,I,M 
3091:       REAL*8 EXPSUM,EXP2SUM,DXI,DYI,DZI,MINE 
3092:       REAL*8 DXI2,DYI2,DZI2 
3093:       REAL*8 TMPH(1) 
3094:       REAL*8 EXPWT(MAXENS) 
3095: C 
3096:       TMPH(1) = HK 
3097: ##IF ALTIX 
3098:       CALL MPI8BARRIER(MPI_COMM_WORLD, IERROR) 
3099:       CALL MPI8ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1), 
3100:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3101:       CALL MPI8ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1), 
3102:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3103:       CALL MPI8ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1), 
3104:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3105:       CALL MPI8ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1), 
3106:      1       1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3107: ##ELSE 
3108:       CALL MPI_BARRIER(MPI_COMM_WORLD, IERROR) 
3109:       CALL MPI_ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1), 
3110:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3111:       CALL MPI_ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1), 
3112:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3113:       CALL MPI_ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1), 
3114:      1       NATOM,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3115:       CALL MPI_ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1), 
3116:      1       1,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERROR) 
3117: ##ENDIF 
3118: C     first calculate mean to use as an offset to avoid numerical 
3119: C     errors in exp(-beta U) for very small (large negative) U. 
3120:       MINE=RCVH(1)+EXPOFF(1) 
3121:       DO I=2,NENSEM 
3122:          MINE = MIN(MINE,RCVH(I)+EXPOFF(I)) 
3123:       ENDDO 
3124: C     convert entries HK(I) to EXP(-beta*HK(I)) 
3125:       EXPSUM=0.0 
3126:       EXP2SUM=0.0 
3127:       DO I=1,NENSEM 
3128:          EXPWT(I)=DEXP(-ENSEBETA*(RCVH(I)+EXPOFF(I)-MINE)) 
3129:          EXPSUM=EXPSUM+EXPWT(I) 
3130:          IF (QEXPBW) EXP2SUM=EXP2SUM+(EXPWT(I))**2 
3131:       ENDDO1573:       ENDDO
3132: C     calculate forces1574: C     calculate forces
3133:       IF (QEXPBW) THEN1575:       IF (QEXPBW) THEN
3134: C        write (outu,*) 'DOING EXPBW'1576: C        write (outu,*) 'DOING EXPBW'
3135: C        call flush(outu)1577: C        call flush(outu)
3136:         DO I=1,NATOM1578:         DO I=1,NATOM
3137:            DXI = 0.01579:            DXI = 0.0
3138:            DYI = 0.01580:            DYI = 0.0
3139:            DZI = 0.01581:            DZI = 0.0
3140:            DXI2 = 0.01582:            DXI2 = 0.0
3176:            DO I=1,NENSEM1618:            DO I=1,NENSEM
3177:                WRITE(ENSEXPU,'(1X,F12.5$)') RCVH(I)1619:                WRITE(ENSEXPU,'(1X,F12.5$)') RCVH(I)
3178:            ENDDO1620:            ENDDO
3179:            WRITE(ENSEXPU,'(1X,F12.5)') HK1621:            WRITE(ENSEXPU,'(1X,F12.5)') HK
3180:         ENDIF1622:         ENDIF
3181:       ENDIF1623:       ENDIF
3182: C     1624: C     
3183:       RETURN1625:       RETURN
3184:       END1626:       END
3185: C1627: C
3186:  
3187:       SUBROUTINE EVB(QX,QY,QZ,DXK,DYK,DZK,HK,NATOM,BUFLEN) 
3188: c     1                RCVDX,RCVDY,RCVDZ,RCVH,BUFLEN) 
3189: C----------------------------------------------------------------------- 
3190: C The nodes share their forces and energy; the total energy is calculated 
3191: c as the smallest eigenvalue of an NENSEM by NENSEM matrix of mm energies. 
3192: c Gradients are calculated on the EVB surface using Feynmann-Hellman;  
3193: c n-state evb functionality coded by david glowacki, sept 2010 
3194: C 
3195: C 
3196: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
3197: ##INCLUDE '~/charmm_fcm/stream.fcm' 
3198: ##INCLUDE '~/charmm_fcm/ensemble.fcm' 
3199: ##INCLUDE '~/charmm_fcm/contrl.fcm' 
3200: ##INCLUDE '~/charmm_fcm/number.fcm' 
3201:       INCLUDE 'mpif.h' 
3202:  
3203: C     arguments: 
3204: C     qx,qy,qz are the coordinates 
3205:       REAL*8 QX(*),QY(*),QZ(*) 
3206:       REAL*8 DXK(*),DYK(*),DZK(*) 
3207:       REAL*8 RCVDX(NEVBS*NATOM),RCVDY(NEVBS*NATOM),RCVDZ(NEVBS*NATOM) 
3208:       REAL*8 RCVH(NEVBS),HK 
3209:       INTEGER NATOM,BUFLEN 
3210:  
3211: C     local: 
3212:       INTEGER IERROR,I,M,J,CTR,II,JJ,NG,GI,CEIDX 
3213:       REAL*8 TMPH(1) 
3214: cccc  EVB coupling matrix 
3215:       REAL*8 EVBMAT(MAXENS,MAXENS) 
3216: cccc  array to hold upper triagonal evb matrix 
3217:       REAL*8 UPTRI(MAXENS*(MAXENS+1)/2) 
3218: cccc  eigenvector matrix, eigenvalue array, and array to hold 
3219: cccc  the eigenvector corresponding to the smallest eigenvalue 
3220:       REAL*8 REVECS(MAXENS,MAXENS),EVALS(MAXENS),SMEVEC(MAXENS) 
3221: Cccc  davidglo scratch arrays for EVB diagonalization 
3222:       REAL*8 S1(MAXENS),S2(MAXENS),S3(MAXENS),S4(MAXENS) 
3223:       REAL*8 S5(MAXENS),S6(MAXENS),S7(MAXENS) 
3224: cccc  davidglo, evb x,y,z forces 
3225:       REAL*8 DXEVB,DYEVB,DZEVB 
3226: cccc  davidglo, distance & parameter labels for gaussian coupling element 
3227:       REAL*8 DIS1(MAXENS*(MAXENS+1)/2),DIS2(MAXENS*(MAXENS+1)/2) 
3228:       REAL*8 dDIS1dX,dDIS1dY,dDIS1dZ,dDIS2dX,dDIS2dY,dDIS2dZ 
3229:       REAL*8 DIFFDIS,dDIFFDISdX,dDIFFDISdY,dDIFFDISdZ 
3230:       REAL*8 dGdDIS1,dGdDIS2 
3231:       REAL*8 GA,GR01,GC1,GR02,GAU,SMALLA,SMALLB,SMALLC 
3232: cccc  davidglo, gradient coupling matrices 
3233:       REAL*8 dMdX(MAXENS,MAXENS),dMdY(MAXENS,MAXENS) 
3234:       REAL*8 dMdZ(MAXENS,MAXENS),diffMX,diffMY,diffMZ 
3235: C 
3236:  
3237: ce----PVG error checking 
3238: ce      write (*,'(A,I)') 'interBeadGroupComm =', interBeadGroupComm 
3239: ce      write (*,'(A,I)') 'interStateGroupComm =', interStateGroupComm 
3240: ce----PVG error checking 
3241: ce      write (*,*) 'rank',WHOIAM,'sub evb interBeadGroupComm =', 
3242: ce     &interBeadGroupComm 
3243: ce      write (*,*) 'rank',WHOIAM,'sub evb interStateGroupComm =', 
3244: ce     &interStateGroupComm 
3245:       
3246: ce----PVG END ERROR CHECKING 
3247: ce----PVG END ERROR CHECKING 
3248:  
3249:  
3250:  
3251:  
3252:  
3253:       TMPH(1) = HK 
3254: ##IF ALTIX 
3255:       CALL MPI8BARRIER(interStateGroupComm, IERROR) 
3256:       CALL MPI8ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1), 
3257:      1       NATOM,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3258:       CALL MPI8ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1), 
3259:      1       NATOM,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3260:       CALL MPI8ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1), 
3261:      1       NATOM,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3262:       CALL MPI8ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1), 
3263:      1       1,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3264: ##ELSE 
3265:       CALL MPI_BARRIER(interStateGroupComm, IERROR) 
3266:       CALL MPI_ALLGATHER(DXK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDX(1), 
3267:      1       NATOM,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3268:       CALL MPI_ALLGATHER(DYK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDY(1), 
3269:      1       NATOM,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3270:       CALL MPI_ALLGATHER(DZK(1),NATOM,MPI_DOUBLE_PRECISION,RCVDZ(1), 
3271:      1       NATOM,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3272:       CALL MPI_ALLGATHER(TMPH(1),1,MPI_DOUBLE_PRECISION,RCVH(1), 
3273:      1       1,MPI_DOUBLE_PRECISION,interStateGroupComm,IERROR) 
3274: ##ENDIF 
3275:  
3276: c---- set up the EVB matrix for an NEVBS-state system 
3277: c---- the diagonal elements are the energies + the user specified offset (rchv+expoff) 
3278: c---- the off-diagonal elements are the coupling elements 
3279:  
3280:       DO II=1,NEVBS  
3281:         EVBMAT(II,II)=0.0 
3282:         DO JJ=II+1,NEVBS  
3283:           EVBMAT(II,JJ)=0.0 
3284:           EVBMAT(JJ,II)=0.0 
3285:         ENDDO 
3286:       ENDDO 
3287:  
3288: ce      WRITE(*,*)'NREAD',NREAD 
3289: ce      DO I=1,NREAD 
3290: ce        WRITE(*,*)'I,EVBIIDX(I),EVBJIDX(I)',I,EVBIIDX(I),EVBJIDX(I) 
3291: ce      ENDDO 
3292:  
3293: ce      WRITE(*,*)'NENSEM',NENSEM      
3294:  
3295:       DO II=1,NEVBS  
3296:         EVBMAT(II,II)=RCVH(II)+EXPOFF(II) 
3297: ce------error checking code 
3298: ce        WRITE(*,*)' II,II:',II,II,RCVH(II-1),EXPOFF(II-1) 
3299: ce------end error checking 
3300:         DO JJ=II+1,NEVBS  
3301: c  ------ initialize the off diagonals to a consistent state 
3302:           EVBMAT(II,JJ)=0.0D0 
3303: c-------- see whether there's an entry for coupling element II,JJ 
3304:           DO J=1,NREAD 
3305: c---------- "II-1" and "JJ-1" arise because the node indices stored in  
3306: c---------- EVBIIDX and EVBJIDX begin with zero; but Fortran arrays begin with 1 
3307:             IF((II).EQ.EVBIIDX(J).AND.(JJ).EQ.EVBJIDX(J)) THEN 
3308:  
3309: c------------ what to do if the entry is a constant 
3310:               IF(COUTYPE(J).EQ.'CONS') THEN 
3311:                 EVBMAT(II,JJ)=EVBMAT(II,JJ) + COUPRM(J,1) 
3312:                  
3313: c------------ what to do if the entry is a gaussian 
3314:               ELSEIF(COUTYPE(J).EQ.'ONED'.OR.COUTYPE(J).EQ.'TWOD'.OR. 
3315:      &        COUTYPE(J).EQ.'DIFF') THEN               
3316:                 GA=COUPRM(J,1) 
3317:                 GR01=COUPRM(J,2) 
3318:                 GC1=COUPRM(J,3) 
3319:                 DIS1(J)=SQRT((QX(ESEL(J,1))-QX(ESEL(J,2)))**2 
3320:      &                      +(QY(ESEL(J,1))-QY(ESEL(J,2)))**2 
3321:      &                      +(QZ(ESEL(J,1))-QZ(ESEL(J,2)))**2) 
3322: ce--------------error checking 
3323: ce                WRITE(*,'(A,I3,A,E)'),"1986: DIS1(",J,")=",DIS1(J) 
3324: ce--------------end error checking 
3325:  
3326: c-------------- what to do if the entry is a ONED gaussian 
3327:                 IF(COUTYPE(J).EQ.'ONED') THEN 
3328:                   GAU=GA*EXP(((DIS1(J)-GR01)**2)/(-2.0D0*GC1**2)) 
3329:                   EVBMAT(II,JJ)=EVBMAT(II,JJ) + GAU 
3330:                    
3331: c-------------- what to do if the entry is a TWOD gaussian 
3332:                 ELSEIF(COUTYPE(J).EQ.'TWOD') THEN 
3333:                   DIS2(J)=SQRT((QX(ESEL2(J,1))-QX(ESEL2(J,2)))**2 
3334:      &                        +(QY(ESEL2(J,1))-QY(ESEL2(J,2)))**2 
3335:      &                        +(QZ(ESEL2(J,1))-QZ(ESEL2(J,2)))**2) 
3336: ce----------------error checking  
3337: ce                  WRITE(*,'(A,I3,A,E)'),"1993: DIS2(",J,")=",DIS2(J) 
3338: ce----------------end error checking  
3339:                   GR02=COUPRM(J,5) 
3340:                   SMALLA=COUPRM(J,7) 
3341:                   SMALLB=COUPRM(J,8) 
3342:                   SMALLC=COUPRM(J,9) 
3343:                   GAU=GA*EXP(-1.0D0*(SMALLA*(DIS1(J)-GR01)**2 
3344:      &              +2.0D0*SMALLB*(DIS1(J)-GR01)*(DIS2(J)-GR02) 
3345:      &              +SMALLC*(DIS2(J)-GR02)**2)) 
3346:                   EVBMAT(II,JJ)=EVBMAT(II,JJ) + GAU 
3347: ce--------------- error checking 
3348: ce                  WRITE(*,*)' GA', GA      
3349: ce                  WRITE(*,*)' DIS1', DIS1(J) 
3350: ce                  WRITE(*,*)' DIS2', DIS2(J) 
3351: ce                  WRITE(*,*)' GR01', GR01 
3352: ce                  WRITE(*,*)' GR02', GR02 
3353: ce                  WRITE(*,*)' SMALLA', SMALLA 
3354: ce                  WRITE(*,*)' SMALLB', SMALLB 
3355: ce                  WRITE(*,*)' SMALLC', SMALLC 
3356: ce                  WRITE(*,*)' EVBMAT(II,JJ)',EVBMAT(II,JJ) 
3357: ce--------------- end error checking 
3358: c----------------- what to fo if the entry is a DIFF gaussian 
3359:                 ELSEIF(COUTYPE(J).EQ.'DIFF') THEN 
3360:                   DIS2(J)=SQRT((QX(ESEL2(J,1))-QX(ESEL2(J,2)))**2 
3361:      &                        +(QY(ESEL2(J,1))-QY(ESEL2(J,2)))**2 
3362:      &                        +(QZ(ESEL2(J,1))-QZ(ESEL2(J,2)))**2) 
3363:                   GAU=GA*EXP(((DIS1(J)-DIS2(J)-GR01)**2) 
3364:      &            /(-2.0D0*GC1**2)) 
3365:                   EVBMAT(II,JJ)=EVBMAT(II,JJ) + GAU 
3366:                 ENDIF 
3367:               ENDIF      
3368:             ENDIF 
3369:           ENDDO 
3370:           EVBMAT(JJ,II)=EVBMAT(II,JJ) 
3371: c-------- davidglo write output to evbfile 
3372: c          WRITE(ENSEXPU,'(A,I1,A,I1,F10.5$)')' V_',I-1,'-',J-1,EVBMAT(II 
3373: c     &,JJ) 
3374:         ENDDO 
3375:       ENDDO 
3376:  
3377: c---- put EVBMAT in an upper triangular vector & diagonalize 
3378:  
3379:       CTR=1 
3380:       DO I=1,NEVBS 
3381:         DO J=I,NEVBS 
3382:           UPTRI(CTR)=EVBMAT(I,J) 
3383: ce------PVG ERROR CHECKING 
3384: ce          WRITE(*,*) 'EVBMAT(',I,J,') =', EVBMAT(I,J) 
3385: ce------PVG ERROR CHECKING 
3386:           CTR=CTR+1 
3387:         ENDDO 
3388:       ENDDO 
3389:  
3390:       CALL DIAGQ(NEVBS,NEVBS,UPTRI,REVECS,S1,S2,S3,S4,EVALS, 
3391:      1S5,S6,S7,0) 
3392:  
3393: ce----error checking 
3394: ce      WRITE(*,*)'EVALS AFTER DIAGQ',EVALS(1),EVALS(2) 
3395: ce----end error checking 
3396:  
3397: c---- search for the smallest eigenvalue, and assign it to the variable  
3398: c---- HK, which is the total energy. For a two state system: 
3399: c---- HK=0.5*(V1+V2-SQRT((V1-V2)**2+4*ENSEBETA**2)) 
3400:  
3401:       HK=EVALS(1) 
3402:       CTR=1 
3403:       DO I=1,NEVBS 
3404:         IF(EVALS(I).LT.HK) THEN 
3405:           HK=EVALS(I) 
3406:           CTR=I 
3407:         ENDIF 
3408:       ENDDO 
3409:  
3410: c---- get the eigenvector corresponding to the smallest eigenvalue 
3411:       DO I=1,NEVBS 
3412:         SMEVEC(I)=REVECS(I,CTR) 
3413: ce------error checking 
3414: ce        write(*,*)'I,SMEVEC(I)',I,SMEVEC(I) 
3415: ce------end error checking 
3416:       ENDDO      
3417:  
3418: c---- get the gradients using the formula from  
3419: c---- P. Lancaster, Numerische Mathematik 6, 377-387 (1964). This is identical 
3420: c---- to using the Hellman-Feynmann approach. 
3421: c---- If M is EVBMAT, its eigenvector matrix is U, and its eigenvalue 
3422: c---- matrix is D (i.e., D=(U^T)*M*U), then dD/dp=(U^T)*(dM/dp)*(U),  
3423: c---- where M depends on some parameter p. In this case, p respresents the  
3424: c---- cartesian directions - x,y,z, so that we have dM/dx, dM/dy, and dM/dz. 
3425: c---- For a two state system, this method gives results identical 
3426: c---- to those obtained with an analytic formula, but is extendable to  
3427: c---- multi-state systems where analytic forms are not available 
3428: c**** note: 
3429: c---- Analytic 2-state gradients (e.g., DXK) would be calculated as follows: 
3430: c---- DXK(I)=0.5*(RCVDX(I)+RCVDX(NATOM+I)-TMPFAC*(RCVDX(I)-RCVDX(NATOM+I))) 
3431: c---- where TMPFAC=(V1-V2)/SQRT((V1-V2)**2+4*ENSEBETA**2) 
3432:  
3433: c---- set up the dM/dp matrices - one for each cartesian direction 
3434:       DO I=1,NATOM 
3435:         DO II=1,NEVBS 
3436: c-------- diagonals are the force at each state (dV/dx,dV/dy,dV/dz) 
3437:           dMdX(II,II)=RCVDX((II-1)*NATOM+I) 
3438:           dMdY(II,II)=RCVDY((II-1)*NATOM+I) 
3439:           dMdZ(II,II)=RCVDZ((II-1)*NATOM+I) 
3440: ce--------error checking 
3441: ce          WRITE(*,'(A,I3,I3,A,E)'),"2071: dMdX(",II,II,")= ", 
3442: ce     &dMdX(II,II) 
3443: ce--------end error checking 
3444: c-------  the off-diagonals are the gradients of the coupling elements 
3445:           DO JJ=II+1,NEVBS 
3446:             dMdX(II,JJ)=0.0D0 
3447:             dMdY(II,JJ)=0.0D0 
3448:             dMdZ(II,JJ)=0.0D0 
3449: c---------- see whether there's an entry for coupling element II,JJ 
3450:             DO J=1,NREAD 
3451:               IF((II).EQ.EVBIIDX(J).AND.(JJ).EQ.EVBJIDX(J)) THEN 
3452: c-------------- if we have gaussian coupling elements, then some of the  
3453: c-------------- off-diagonals of the dM/dp matrices are not zero 
3454:                 IF(COUTYPE(J).EQ.'ONED'.OR.COUTYPE(J).EQ.'TWOD') THEN 
3455:                   IF(ESEL(J,1).EQ.I.OR.ESEL(J,2).EQ.I) THEN 
3456: c------------------ dDIS1/dx, dDIS1/dy, & dDIS1/dz elements w.r.t DIS1:  
3457:                     dDIS1dX=(QX(ESEL(J,1))-QX(ESEL(J,2)))/DIS1(J) 
3458:                     dDIS1dY=(QY(ESEL(J,1))-QY(ESEL(J,2)))/DIS1(J) 
3459:                     dDIS1dZ=(QZ(ESEL(J,1))-QZ(ESEL(J,2)))/DIS1(J) 
3460:                     IF(ESEL(J,2).EQ.I)THEN 
3461:                       dDIS1dX=-1.0D0*dDIS1dX 
3462:                       dDIS1dY=-1.0D0*dDIS1dY 
3463:                       dDIS1dZ=-1.0D0*dDIS1dZ 
3464:                     ENDIF 
3465: c---------------- what to do if the entry is a ONED gaussian 
3466:                     IF(COUTYPE(J).EQ.'ONED') THEN 
3467: c-------------------- get the parameters for 1d gaussian 
3468:                       GA=COUPRM(J,1) 
3469:                       GR01=COUPRM(J,2) 
3470:                       GC1=COUPRM(J,3) 
3471:                       GAU=GA*EXP(((DIS1(J)-GR01)**2)/(-2.0D0*GC1**2)) 
3472: ce------------------- error checking 
3473: ce                    WRITE(*,'(A,I3,A,I3,A,I3,A,I3,A,E)'),"2111: I= ",I, 
3474: ce     &" GI= ",GI," II= ",II," JJ= ",JJ," GAU= ",GAU 
3475: ce------------------- end error checking 
3476: c-------------------- calculate the chain rule gradients  
3477:                       diffMX=dDIS1dX*(GR01-DIS1(J))*GAU/(GC1**2) 
3478:                       diffMY=dDIS1dY*(GR01-DIS1(J))*GAU/(GC1**2) 
3479:                       diffMZ=dDIS1dZ*(GR01-DIS1(J))*GAU/(GC1**2) 
3480: c-------------------- update the matrix elements 
3481:                       dMdX(II,JJ)=dMdX(II,JJ)+diffMX 
3482:                       dMdY(II,JJ)=dMdY(II,JJ)+diffMY 
3483:                       dMdZ(II,JJ)=dMdZ(II,JJ)+diffMZ 
3484: ce--------------------error checking 
3485: ce                    WRITE(*,'(A,I3,A,I3,A,I3,A,I3,A,E)'),"2121: I= ",I, 
3486: ce     &" GI= ",GI," II= ",II," JJ= ",JJ," diffMZ ",diffMZ 
3487: ce------------------- end error checking 
3488:  
3489: c------------------ what to do if the entry is a TWOD gaussian 
3490:                     ELSEIF(COUTYPE(J).EQ.'TWOD') THEN 
3491: c-------------------- get the parameters for 2d gaussian 
3492:                       GA=COUPRM(J,1) 
3493:                       GR01=COUPRM(J,2) 
3494:                       GR02=COUPRM(J,5) 
3495:                       SMALLA=COUPRM(J,7) 
3496:                       SMALLB=COUPRM(J,8) 
3497:                       SMALLC=COUPRM(J,9) 
3498:                       GAU=GA*EXP(-1.0D0*(SMALLA*(DIS1(J)-GR01)**2 
3499:      &               +2.0D0*SMALLB*(DIS1(J)-GR01)*(DIS2(J)-GR02) 
3500:      &               +SMALLC*(DIS2(J)-GR02)**2)) 
3501: ce--------------------error checking 
3502: ce                    WRITE(*,'(A,I3)'),"2232: I= ",I 
3503: ce                    WRITE(*,*)'X1',QX(ESEL(J,1)),'X2',QX(ESEL(J,2)) 
3504: ce                    WRITE(*,*)'Y1',QY(ESEL(J,1)),'Y2',QY(ESEL(J,2)) 
3505: ce                    WRITE(*,*)'Z1',QZ(ESEL(J,1)),'Z2',QZ(ESEL(J,2)) 
3506: ce                    WRITE(*,*)'DIS1 = ',DIS1(J) 
3507: ce                    WRITE(*,*)'GR01 = ',GR01 
3508: ce                    WRITE(*,*)'DIS2 = ',DIS2(J) 
3509: ce                    WRITE(*,*)'GR02 = ',GR02 
3510: ce                    WRITE(*,*)'GAU = ',GAU 
3511: ce--------------------end error checking 
3512: c-------------------- calculate the chain rule gradients: 
3513:                       dGdDIS1=-2.0D0*SMALLA*(DIS1(J)-GR01) 
3514:      &                      -2.0D0*SMALLB*(DIS2(J)-GR02) 
3515:                       diffMX=GAU*dGdDIS1*dDIS1dX 
3516:                       diffMY=GAU*dGdDIS1*dDIS1dY 
3517:                       diffMZ=GAU*dGdDIS1*dDIS1dZ 
3518: ce--------------------error checking 
3519: ce                    WRITE(*,*)'diffMX',diffMX 
3520: ce                    WRITE(*,*)'diffMY',diffMY 
3521: ce                    WRITE(*,*)'diffMZ',diffMZ 
3522: ce--------------------end error checking 
3523: c-------------------- update the matrix elements 
3524:                       dMdX(II,JJ)=dMdX(II,JJ)+diffMX 
3525:                       dMdY(II,JJ)=dMdY(II,JJ)+diffMY 
3526:                       dMdZ(II,JJ)=dMdZ(II,JJ)+diffMZ 
3527:                     ENDIF 
3528:                   ENDIF 
3529:                 ENDIF 
3530:                 IF(COUTYPE(J).EQ.'TWOD') THEN 
3531:                   IF(ESEL2(J,1).EQ.I.OR.ESEL2(J,2).EQ.I) THEN 
3532: c------------------ dDIS2/dx, dDIS2/dy, & dDIS2/dz elements w.r.t DIS2: 
3533:                     dDIS2dX=(QX(ESEL2(J,1))-QX(ESEL2(J,2)))/DIS2(J) 
3534:                     dDIS2dY=(QY(ESEL2(J,1))-QY(ESEL2(J,2)))/DIS2(J) 
3535:                     dDIS2dZ=(QZ(ESEL2(J,1))-QZ(ESEL2(J,2)))/DIS2(J) 
3536:                     IF(ESEL2(J,2).EQ.I)THEN 
3537:                       dDIS2dX=-1.0D0*dDIS2dX 
3538:                       dDIS2dY=-1.0D0*dDIS2dY 
3539:                       dDIS2dZ=-1.0D0*dDIS2dZ 
3540:                     ENDIF 
3541: c------------------ get the parameters for 2d gaussian 
3542:                     GA=COUPRM(J,1) 
3543:                     GR01=COUPRM(J,2) 
3544:                     GR02=COUPRM(J,5) 
3545:                     SMALLA=COUPRM(J,7) 
3546:                     SMALLB=COUPRM(J,8) 
3547:                     SMALLC=COUPRM(J,9) 
3548:                     GAU=GA*EXP(-1.0D0*(SMALLA*(DIS1(J)-GR01)**2 
3549:      &               +2.0D0*SMALLB*(DIS1(J)-GR01)*(DIS2(J)-GR02) 
3550:      &               +SMALLC*(DIS2(J)-GR02)**2)) 
3551: ce--------------------error checking 
3552: ce                    WRITE(*,'(A,I3)'),"2280: I= ",I 
3553: ce                    WRITE(*,*)'X1',QX(ESEL2(J,1)),'X2',QX(ESEL2(J,2)) 
3554: ce                    WRITE(*,*)'Y1',QY(ESEL2(J,1)),'Y2',QY(ESEL2(J,2)) 
3555: ce                    WRITE(*,*)'Z1',QZ(ESEL2(J,1)),'Z2',QZ(ESEL2(J,2)) 
3556: ce                    WRITE(*,*)'DIS1 = ',DIS1(J) 
3557: ce                    WRITE(*,*)'GR01 = ',GR01 
3558: ce                    WRITE(*,*)'DIS2 = ',DIS2(J) 
3559: ce                    WRITE(*,*)'GR02 = ',GR02 
3560: ce                    WRITE(*,*)'GAU = ',GAU 
3561: ce--------------------end error checking 
3562: c------------------ calculate the chain rule gradients: 
3563:                     dGdDIS2=-2.0D0*SMALLB*(DIS1(J)-GR01) 
3564:      &                    -2.0D0*SMALLC*(DIS2(J)-GR02) 
3565:                     diffMX=GAU*dGdDIS2*dDIS2dX 
3566:                     diffMY=GAU*dGdDIS2*dDIS2dY 
3567:                     diffMZ=GAU*dGdDIS2*dDIS2dZ 
3568: ce------------------error checking 
3569: ce                  WRITE(*,*)'diffMX',diffMX 
3570: ce                  WRITE(*,*)'diffMY',diffMY 
3571: ce                  WRITE(*,*)'diffMZ',diffMZ 
3572: ce------------------end error checking 
3573: c------------------ update the matrix elements 
3574:                     dMdX(II,JJ)=dMdX(II,JJ)+diffMX 
3575:                     dMdY(II,JJ)=dMdY(II,JJ)+diffMY 
3576:                     dMdZ(II,JJ)=dMdZ(II,JJ)+diffMZ 
3577:                   ENDIF 
3578:                 ELSEIF(COUTYPE(J).EQ.'DIFF') THEN 
3579:                   IF(ESEL(J,1).EQ.I.OR.ESEL(J,2).EQ.I) THEN 
3580: c------------------ dDIFFDIS/dx, dDIFFDIS1/dy, & dDIFFDIS1/dz elements w.r.t DIFFDIS:  
3581:                     dDIFFDISdX=(QX(ESEL(J,1))-QX(ESEL(J,2)))/DIS1(J) 
3582:                     dDIFFDISdY=(QY(ESEL(J,1))-QY(ESEL(J,2)))/DIS1(J) 
3583:                     dDIFFDISdZ=(QZ(ESEL(J,1))-QZ(ESEL(J,2)))/DIS1(J) 
3584:                     IF(ESEL(J,2).EQ.I) THEN 
3585:                       dDIFFDISdX=-1.0D0*dDIFFDISdX 
3586:                       dDIFFDISdY=-1.0D0*dDIFFDISdY 
3587:                       dDIFFDISdZ=-1.0D0*dDIFFDISdZ 
3588:                     ENDIF 
3589:                     GA=COUPRM(J,1) 
3590:                     GR01=COUPRM(J,2) 
3591:                     GC1=COUPRM(J,3) 
3592:                     DIFFDIS=DIS1(J)-DIS2(J) 
3593:                     GAU=GA*EXP(((DIFFDIS-GR01)**2)/(-2.0D0*GC1**2)) 
3594: c-------------------calculate the chain rule gradients  
3595:                     diffMX=dDIFFDISdX*(GR01-DIFFDIS)*GAU/(GC1**2) 
3596: ce----------------error checking 
3597: ce                  WRITE(*,'(A,E)'),"HELLO PVG diffMX=",diffMX 
3598: ce----------------end error checking 
3599:                     diffMY=dDIFFDISdY*(GR01-DIFFDIS)*GAU/(GC1**2) 
3600:                     diffMZ=dDIFFDISdZ*(GR01-DIFFDIS)*GAU/(GC1**2) 
3601: c-------------------update the matrix elements 
3602:                     dMdX(II,JJ)=dMdX(II,JJ)+diffMX 
3603:                     dMdY(II,JJ)=dMdY(II,JJ)+diffMY 
3604:                     dMdZ(II,JJ)=dMdZ(II,JJ)+diffMZ 
3605:                   ENDIF 
3606:                   IF(ESEL2(J,1).EQ.I.OR.ESEL2(J,2).EQ.I) THEN 
3607: c------------------ dDIFFDIS/dx, dDIFFDIS1/dy, & dDIFFDIS1/dz elements w.r.t DIFFDIS:  
3608:                     dDIFFDISdX=(QX(ESEL2(J,1))-QX(ESEL2(J,2)))/DIS2(J) 
3609:                     dDIFFDISdY=(QY(ESEL2(J,1))-QY(ESEL2(J,2)))/DIS2(J) 
3610:                     dDIFFDISdZ=(QZ(ESEL2(J,1))-QZ(ESEL2(J,2)))/DIS2(J) 
3611:                     IF(ESEL2(J,1).EQ.I) THEN 
3612:                       dDIFFDISdX=-1.0D0*dDIFFDISdX 
3613:                       dDIFFDISdY=-1.0D0*dDIFFDISdY 
3614:                       dDIFFDISdZ=-1.0D0*dDIFFDISdZ 
3615:                     ENDIF 
3616:                     GA=COUPRM(J,1) 
3617:                     GR01=COUPRM(J,2) 
3618:                     GC1=COUPRM(J,3) 
3619:                     DIFFDIS=DIS1(J)-DIS2(J) 
3620:                     GAU=GA*EXP(((DIFFDIS-GR01)**2)/(-2.0D0*GC1**2)) 
3621: c-------------------calculate the chain rule gradients  
3622:                     diffMX=dDIFFDISdX*(GR01-DIFFDIS)*GAU/(GC1**2) 
3623:                     diffMY=dDIFFDISdY*(GR01-DIFFDIS)*GAU/(GC1**2) 
3624:                     diffMZ=dDIFFDISdZ*(GR01-DIFFDIS)*GAU/(GC1**2) 
3625: c-------------------update the matrix elements 
3626:                     dMdX(II,JJ)=dMdX(II,JJ)+diffMX 
3627:                     dMdY(II,JJ)=dMdY(II,JJ)+diffMY 
3628:                     dMdZ(II,JJ)=dMdZ(II,JJ)+diffMZ 
3629:                   ENDIF 
3630:                 ENDIF 
3631:               ENDIF 
3632:             ENDDO 
3633:           ENDDO 
3634: c-------- assign the other elements of the symmetric matrix 
3635:           DO JJ=II+1,NEVBS 
3636:             dMdX(JJ,II)=dMdX(II,JJ)   
3637:             dMdY(JJ,II)=dMdY(II,JJ)   
3638:             dMdZ(JJ,II)=dMdZ(II,JJ) 
3639:           ENDDO  
3640:         ENDDO 
3641: c------ now calculate the forces acting on the atom on the EVB surface. The 
3642: c------ procedure below effectively amounts to calculating (v^T)*(dM/dp)*v, where 
3643: c------ v is the eigenvector corresponding to the smallest eigenvalue (SMEVEC) 
3644:         DXEVB=0.0D0 
3645:         DYEVB=0.0D0 
3646:         DZEVB=0.0D0 
3647:         DO II=1,NEVBS 
3648:           DO JJ=1,NEVBS 
3649:             DXEVB=DXEVB+SMEVEC(II)*SMEVEC(JJ)*dMdX(II,JJ) 
3650:             DYEVB=DYEVB+SMEVEC(II)*SMEVEC(JJ)*dMdY(II,JJ) 
3651:             DZEVB=DZEVB+SMEVEC(II)*SMEVEC(JJ)*dMdZ(II,JJ) 
3652:           ENDDO 
3653:         ENDDO 
3654:         DXK(I)=DXEVB 
3655:         DYK(I)=DYEVB 
3656:         DZK(I)=DZEVB 
3657:       ENDDO 
3658:  
3659: c---- davidglo - end forces loop 
3660: c---- write output ... 
3661:       IF ((IOLEV.GE.0).AND.(ENSEXPU.GT.0)) THEN 
3662:         IF(MOD(MDSTEP,NPRINT).EQ.0) THEN 
3663:          WRITE(ENSEXPU,'(I10$)') MDSTEP 
3664:          DO I=1,NEVBS 
3665:            WRITE(ENSEXPU,'(1X,F12.5$)') RCVH(I)+EXPOFF(I) 
3666:          ENDDO 
3667:          WRITE(ENSEXPU,'(1X,F12.5)') HK 
3668:         ENDIF 
3669:       ENDIF 
3670: C      
3671:       RETURN 
3672:       END 
3673: C 
3674:  
3675:       SUBROUTINE ERDCMD(COMLYN,MXCMS2,COMLEN,UNIT,EOF,1628:       SUBROUTINE ERDCMD(COMLYN,MXCMS2,COMLEN,UNIT,EOF,
3676:      $                  QCMPRS,QPRINT,ECHOST)1629:      $                  QCMPRS,QPRINT,ECHOST)
3677: C-----------------------------------------------------------------------1630: C-----------------------------------------------------------------------
3678: C     parallel input version of rdcmnd1631: C     parallel input version of rdcmnd
3679: ##INCLUDE '~/charmm_fcm/impnon.fcm'1632: ##INCLUDE '~/charmm_fcm/impnon.fcm'
3680: ##INCLUDE '~/charmm_fcm/dimens.fcm'1633: ##INCLUDE '~/charmm_fcm/dimens.fcm'
3681: ##INCLUDE '~/charmm_fcm/exfunc.fcm'1634: ##INCLUDE '~/charmm_fcm/exfunc.fcm'
3682: C1635: C
3683:       INTEGER MXCMS2,COMLEN,UNIT1636:       INTEGER MXCMS2,COMLEN,UNIT
3684:       CHARACTER*(*) COMLYN1637:       CHARACTER*(*) COMLYN
3836: ##INCLUDE '~/charmm_fcm/dimens.fcm'1789: ##INCLUDE '~/charmm_fcm/dimens.fcm'
3837: ##INCLUDE '~/charmm_fcm/ensemble.fcm'1790: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
3838: C1791: C
3839:       INCLUDE 'mpif.h'          !##.not.TERRA1792:       INCLUDE 'mpif.h'          !##.not.TERRA
3840:       INTEGER STATUS1793:       INTEGER STATUS
3841: C1794: C
3842:       IF(NENSEM.EQ.1) RETURN1795:       IF(NENSEM.EQ.1) RETURN
3843: ##IF ALTIX1796: ##IF ALTIX
3844:       CALL MPI8BCAST(ARRAY,8*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)1797:       CALL MPI8BCAST(ARRAY,8*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
3845: ##ELSE1798: ##ELSE
 1799: ##IF INTEGER8
 1800:       CALL MPI_BCAST(ARRAY,8*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
 1801: ##ELSE
3846:       CALL MPI_BCAST(ARRAY,4*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)1802:       CALL MPI_BCAST(ARRAY,4*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
3847: ##ENDIF1803: ##ENDIF
 1804: ##ENDIF
3848:       RETURN1805:       RETURN
3849:       END1806:       END
3850: C1807: C
3851:       SUBROUTINE PSND8(ARRAY, LENGTH)1808:       SUBROUTINE PSND8(ARRAY, LENGTH)
3852: C-----------------------------------------------------------------------1809: C-----------------------------------------------------------------------
3853: C This routine performs on node 0 broadcast to all other nodes1810: C This routine performs on node 0 broadcast to all other nodes
3854: C and receive from node 0 on all other nodes.1811: C and receive from node 0 on all other nodes.
3855: C Usually called after read on node 0. For Real*8 precision arrays.1812: C Usually called after read on node 0. For Real*8 precision arrays.
3856: C1813: C
3857: ##INCLUDE '~/charmm_fcm/impnon.fcm'1814: ##INCLUDE '~/charmm_fcm/impnon.fcm'
3961: C1918: C
3962:       comm4 = comm1919:       comm4 = comm
3963:       datatype4 = datatype1920:       datatype4 = datatype
3964:       count4 = count1921:       count4 = count
3965:       dest4 = dest1922:       dest4 = dest
3966:       tag4 = tag1923:       tag4 = tag
3967: C1924: C
3968:       Call MPI_SEND(buf, count4, datatype4, dest4, tag4, comm4, ierror4)1925:       Call MPI_SEND(buf, count4, datatype4, dest4, tag4, comm4, ierror4)
3969: C1926: C
3970:       ierror = ierror41927:       ierror = ierror4
 1928: C
 1929:       Return
 1930:       End
 1931: C
 1932:       Subroutine MPI8RECV(buf, count, datatype, source , tag, comm, 
 1933:      $    status, ierror)
 1934: C-----------------------------------------------------------------------
 1935: ##INCLUDE '~/charmm_fcm/impnon.fcm'
 1936:       INCLUDE 'mpif.h'
 1937:       INTEGER count, datatype, source, tag, comm, ierror, 
 1938:      $     status(*), ii
 1939:       INTEGER*4 count4, datatype4, source4, tag4, comm4, ierror4,
3971:      $     status4(MPI_STATUS_SIZE)1940:      $     status4(MPI_STATUS_SIZE)
3972:       REAL*8 buf(*)1941:       REAL*8 buf(*)
3973: C1942: C
3974:       comm4 = comm1943:       comm4 = comm
3975:       datatype4 = datatype1944:       datatype4 = datatype
3976:       count4 = count1945:       count4 = count
3977:       source4 = source1946:       source4 = source
3978:       tag4 = tag1947:       tag4 = tag
3979: C1948: C
3980:       Call MPI_RECV(buf, count4, datatype4, source4, tag4, comm4, 1949:       Call MPI_RECV(buf, count4, datatype4, source4, tag4, comm4, 
4057:      $     sendcount4,sendtype4,2026:      $     sendcount4,sendtype4,
4058:      $     recvbuf,recvcount4,recvtype4,2027:      $     recvbuf,recvcount4,recvtype4,
4059:      $     root4,comm4,ierror4)2028:      $     root4,comm4,ierror4)
4060: C2029: C
4061:       ierror = ierror42030:       ierror = ierror4
4062: C2031: C
4063:       Return2032:       Return
4064:       END2033:       END
4065: C2034: C
4066: ##ENDIF2035: ##ENDIF
 2036: c
 2037: c
 2038:       SUBROUTINE GCOMB(X,N)
 2039: C-----------------------------------------------------------------------
 2040: C
 2041: C     Global COMBine routine.
 2042: C
 2043: ##INCLUDE '~/charmm_fcm/impnon.fcm'
 2044: C
 2045:       REAL*8 X(*)
 2046: C
 2047: ##INCLUDE '~/charmm_fcm/dimens.fcm'
 2048: C##INCLUDE '~/charmm_fcm/parallel.fcm'
 2049: ##INCLUDE '~/charmm_fcm/ensemble.fcm'
 2050: C
 2051: ##INCLUDE '~/charmm_fcm/exfunc.fcm'
 2052: ##INCLUDE '~/charmm_fcm/stack.fcm'
 2053: ##IF TIMESTAMP
 2054: ##INCLUDE '~/charmm_fcm/timer.fcm'
 2055: ##ENDIF
 2056: C
 2057: ##IF CMPI
 2058:       INTEGER STATUS,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD
 2059:       PARAMETER(MPI_DOUBLE_PRECISION=0,MPI_SUM=0,MPI_COMM_WORLD=0)
 2060: ##ELSE
 2061:       INCLUDE  'mpif.h'
 2062:       INTEGER STATUS
 2063: ##ENDIF
 2064: C
 2065:       INTEGER N, W, I
 2066:       INTEGER TAG      !##TIMESTAMP
 2067: C
 2068: ##IF TIMESTAMP
 2069:       IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
 2070:      $     'GCOMB called from CHARMM'//CHAR(0),N)
 2071:       IF (ICMRW.EQ.1) THEN
 2072:           CALL GRDFMM(X,N,'GCOMB    ')
 2073:           RETURN
 2074:       ELSE IF (ICMRW.EQ.2) THEN
 2075:           CALL WRTFMR(X,'GCOMB    ',TAG,N)
 2076:           RETURN
 2077:       ENDIF
 2078: ##ENDIF
 2079: C
 2080:       IF(NENSEM.EQ.1) RETURN
 2081: ##IF DDI
 2082:       CALL DDI_GSUMF(3334,X,N)
 2083: ##ELIF CMPI
 2084:       W=ALLSTK(IREAL8(N))
 2085:       CALL CIMPI_ALLRED(X,STACK(W),N,MPI_DOUBLE_PRECISION,MPI_SUM,
 2086:      @                   MPI_COMM_WORLD,STATUS)
 2087:       CALL FRESTK(IREAL8(N))
 2088: ##ELIF CMPIx
 2089: C     Just checking...
 2090:       W=ALLSTK(IREAL8(N))
 2091:       JPARPT(0)=0
 2092:       DO I = 1, NENSEM
 2093:          JPARPT(I)=N*I/NENSEM
 2094:       ENDDO
 2095: ##IF GENCOMM
 2096:       CALL AGVDGS(X,JPARPT,STACK(W),NEMSEM,WHOIAM,IPPMAP)
 2097:       CALL AGVDGB(X,JPARPT,NENSEM,WHOIAM,IPPMAP)
 2098: ##ELSE
 2099:       CALL AVDGSUM(X,JPARPT,STACK(W),NODDIM,MYNOD,IPPMAP)
 2100:       CALL AVDGBR(X,JPARPT,NODDIM,MYNOD,IPPMAP)
 2101: ##ENDIF
 2102:       CALL FRESTK(IREAL8(N))
 2103: ##ELSE
 2104: C     Uses compliant calls.
 2105:       W=ALLSTK(IREAL8(N))
 2106: ##IF ALTIX_MPI IBMAIX_MPI64 IA64_MPI
 2107:       CALL MPI8ALLREDUCE(X,STACK(W),N,MPI_DOUBLE_PRECISION,MPI_SUM,
 2108:      @                   MPI_COMM_WORLD,STATUS)
 2109: ##ELSE
 2110:       CALL MPI_ALLREDUCE(X,STACK(W),N,MPI_DOUBLE_PRECISION,MPI_SUM,
 2111:      @                   MPI_COMM_WORLD,STATUS)
 2112: ##ENDIF
 2113:       CALL COPYR8(STACK(W),X,N)
 2114:       CALL FRESTK(IREAL8(N))
 2115: C
 2116: ##ENDIF
 2117: C
 2118:       RETURN
 2119:       END
 2120: C
4067: ##ELSE2121: ##ELSE
4068: C====================== NO ENSEMBLES ... ==============================2122: C====================== NO ENSEMBLES ... ==============================
4069:       SUBROUTINE ENSCMD(COMLYN,COMLEN)2123:       SUBROUTINE ENSCMD(COMLYN,COMLEN)
4070: C----------------------------------------------------------------------2124: C----------------------------------------------------------------------
4071:       CHARACTER*(*) COMLYN2125:       CHARACTER*(*) COMLYN
4072:       INTEGER COMLEN2126:       INTEGER COMLEN
4073: ##INCLUDE '~/charmm_fcm/stream.fcm'2127: ##INCLUDE '~/charmm_fcm/stream.fcm'
4074:       WRITE (OUTU, '(A)') 'NO ENSEMBLE CODE COMPILED'2128:       WRITE (OUTU, '(A)') 'NO ENSEMBLE CODE COMPILED'
4075:       RETURN2129:       RETURN
4076:       END2130:       END
4077: C2131: C
4078:       SUBROUTINE ENSFIN2132:       SUBROUTINE ENSFIN
4079:       RETURN2133:       RETURN
4080:       END2134:       END
4081: C2135: C
4082:       SUBROUTINE ENSCHK2136:       SUBROUTINE ENSCHK
4083:       RETURN2137:       RETURN
4084:       END2138:       END
4085: ##ENDIF2139: ##ENDIF
4086:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0