hdiff output

r16977/eef1.fcm 2017-01-21 10:32:52.359188235 +0000 r16976/eef1.fcm 2017-01-21 10:32:54.419188235 +0000
  1: CHARMM Element source/fcm/eef1.fcm $Revision: 1.4 $  1: CHARMM Element source/fcm/eef1.fcm $Revision: 1.1.1.1 $
  2: ##IF ASPENER  2: ##IF ASPENER
  3: C  3: C
  4: C For the solvation calculations  4: C For the solvation calculations
  5: C  5: C
  6: C       SLVITC(atom types)= pointers to the correct entry in VOLM,GREF  6: C       SLVITC(atom types)= pointers to the correct entry in VOLM,GREF
  7: C       VOLM(atom types)=   volumes for each atom type  7: C       VOLM(atom types)=   volumes for each atom type
  8: C       GREF(atom types)=   reference solvation free energies  8: C       GREF(atom types)=   reference solvation free energies
  9: C       GFREE(atom types)=  solvation free energy of "free" group   9: C       GFREE(atom types)=  solvation free energy of "free" group 
 10: C 10: C
 11:       INTEGER SLVITC 11:         INTEGER SLVITC
 12:       INTEGER NSMTH 12:         REAL*8 VOLM,GREF,CTFSLV,GSOLV,GFREE,SIGW,FSIGW
 13:       COMMON /SLVI/ SLVITC(MAXATC),NSMTH 13:         REAL*8 HREF,CPREF,MYEXP
 14:  14:         REAL*8 VOLMI,GREFI,GFREEI,SIGWI,FSIGWI,VDWRI
 15:       REAL*8 VOLM,GREF,CTFSLV,GSOLV,GFREE,SIGW,FSIGW 15:         LOGICAL DOEEF1
 16:       REAL*8 HREF,CPREF,MYEXP 16:         COMMON /SLV/ SLVITC(MAXATC),VOLM(MAXATC*2),GREF(MAXATC*2),
 17:       REAL*8 VOLMI,GREFI,GFREEI,SIGWI,FSIGWI,VDWRI 
 18:       REAL*8 GREFI2,GFREEI2,HREFI,HREFI2,CPREFI,CPREFI2,SIGWI2 
 19:       REAL*8 FSIGWI2,WIDTH,LAM,DFDZ 
 20:       REAL*8 RDEBYE 
 21:       REAL*8 GALPHA,PSI0,GKAPPA,OFFST,TEMPR,VOLT,ACONS 
 22:       REAL*8 AEMPIR 
 23:       COMMON /SLVR/ VOLM(MAXATC*2),GREF(MAXATC*2), 
 24:      &     CTFSLV,GSOLV(MAXAIM),GFREE(MAXATC*2),SIGW(MAXATC*2), 17:      &     CTFSLV,GSOLV(MAXAIM),GFREE(MAXATC*2),SIGW(MAXATC*2),
 25:      &     FSIGW(MAXATC*2),HREF(MAXATC*2),CPREF(MAXATC*2),VOLMI(MAXAIM), 18:      &     FSIGW(MAXATC*2),HREF(MAXATC*2),CPREF(MAXATC*2),VOLMI(MAXAIM),
 26:      &     GREFI(MAXAIM),GFREEI(MAXAIM),SIGWI(MAXAIM),FSIGWI(MAXAIM), 19:      &     GREFI(MAXAIM),GFREEI(MAXAIM),SIGWI(MAXAIM),FSIGWI(MAXAIM),
 27:      &     VDWRI(MAXAIM),MYEXP(350), 20:      &     VDWRI(MAXAIM),MYEXP(350),
 28:      &     GREFI2(MAXA),GFREEI2(MAXA),HREFI(MAXA),HREFI2(MAXA), 21:      &     DOEEF1
 29:      &     CPREFI(MAXA),CPREFI2(MAXA),SIGWI2(MAXA),FSIGWI2(MAXA), 
 30:      &     WIDTH,LAM(MAXA),DFDZ(MAXA),AEMPIR, 
 31:      &     RDEBYE,GALPHA,PSI0,GKAPPA,OFFST,TEMPR,VOLT,ACONS 
 32:  
 33:       LOGICAL DOEEF1 
 34:       LOGICAL LMEMBR 
 35:       LOGICAL LDEBYE 
 36:       LOGICAL LGOUY,LVOLT 
 37:       COMMON /SLVL/ DOEEF1,LMEMBR,LDEBYE,LGOUY,LVOLT 
 38:  
 39:       CHARACTER*4 SLVNT,SLVNT2 
 40:       COMMON /SLVC/ SLVNT,SLVNT2 
 41:  
 42: ##IF SAVEFCM 
 43:       SAVE /SLVI/ 
 44:       SAVE /SLVR/ 
 45:       SAVE /SLVL/ 
 46:       SAVE /SLVC/ 
 47: ##ENDIF 
 48: ##ENDIF 22: ##ENDIF
 49: C 23: C


r16977/eef1.src 2017-01-21 10:32:52.647188235 +0000 r16976/eef1.src 2017-01-21 10:32:54.783188235 +0000
 35: ##INCLUDE '~/charmm_fcm/consta.fcm' 35: ##INCLUDE '~/charmm_fcm/consta.fcm'
 36: ##INCLUDE '~/charmm_fcm/psf.fcm' 36: ##INCLUDE '~/charmm_fcm/psf.fcm'
 37: ##INCLUDE '~/charmm_fcm/comand.fcm' 37: ##INCLUDE '~/charmm_fcm/comand.fcm'
 38: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 38: ##INCLUDE '~/charmm_fcm/exfunc.fcm'
 39: ##INCLUDE '~/charmm_fcm/param.fcm' 39: ##INCLUDE '~/charmm_fcm/param.fcm'
 40: ##INCLUDE '~/charmm_fcm/stream.fcm' 40: ##INCLUDE '~/charmm_fcm/stream.fcm'
 41: ##INCLUDE '~/charmm_fcm/number.fcm' 41: ##INCLUDE '~/charmm_fcm/number.fcm'
 42: ##INCLUDE '~/charmm_fcm/eef1.fcm' 42: ##INCLUDE '~/charmm_fcm/eef1.fcm'
 43:         CHARACTER*4 WRD 43:         CHARACTER*4 WRD
 44:         INTEGER UNPAR,I 44:         INTEGER UNPAR,I
 45: c Gouy-Chapman variables 45:         REAL*8 TEMPR
 46:         INTEGER VALENCE 
 47:         REAL*8 ANFR,AREA,CONC,RVAL 
 48:  46: 
 49:       WRD= NEXTA4(COMLYN,COMLEN) 47:       WRD= NEXTA4(COMLYN,COMLEN)
 50:       IF (WRD.eq.'PRIN') THEN 48:       IF (WRD.eq.'PRIN') THEN
 51:           CALL SLVPRINT(COMLYN,COMLEN) 49:           CALL SLVPRINT(COMLYN,COMLEN)
 52:       ELSE 50:       ELSE
 53:           DOEEF1=.true. 51:           DOEEF1=.true.
 54: C jmc          CTFSLV= 81.0d0                           ! 9^2 52: C jmc          CTFSLV= 81.0d0                           ! 9^2
 55: C jdb redefines this later. 53: C jdb redefines this later.
 56:           TEMPR= GTRMF(COMLYN,COMLEN,'TEMP',ROOMT) 54:           TEMPR= GTRMF(COMLYN,COMLEN,'TEMP',ROOMT)
 57:           LMEMBR = (INDXA(COMLYN,COMLEN,'MEMB') .GT. 0) 
 58:           RDEBYE= GTRMF(COMLYN,COMLEN,'IONI',ZERO) 
 59:           IF (RDEBYE.GT.1.E-10) THEN 
 60:               LDEBYE=.TRUE. 
 61: c The following gives Debye length in A, assuming water as solvent 
 62:               RDEBYE= SQRT(0.0316*TEMPR/RDEBYE) 
 63:               WRITE (OUTU,*) 'Debye Length is: ', RDEBYE 
 64:           ENDIF 
 65:           SLVNT= GTRMA(COMLYN,COMLEN,'SLVT') 
 66:           IF (SLVNT.eq.'    ') SLVNT= 'WATE'              !default water 
 67:           IF (LMEMBR) THEN 
 68:              SLVNT2= GTRMA(COMLYN,COMLEN,'SLV2') 
 69:              WIDTH= GTRMF(COMLYN,COMLEN,'WIDT',THIRTY) 
 70:              NSMTH= GTRMI(COMLYN,COMLEN,'NSMT',10) 
 71:              AEMPIR=0.85 
 72:              AEMPIR= GTRMF(COMLYN,COMLEN,'AEMP',AEMPIR) 
 73:              LGOUY = (INDXA(COMLYN,COMLEN,'GOUY') .GT. 0)      !oct03 
 74:              VOLT = GTRMF(COMLYN,COMLEN,'VOLT',ZERO) 
 75:              IF (VOLT.GT.0.0d0) LVOLT=.TRUE. 
 76:           IF (LGOUY.OR.LVOLT) THEN 
 77:              CONC= GTRMF(COMLYN,COMLEN,'CONC',ONE) 
 78:              VALENCE= GTRMI(COMLYN,COMLEN,'VALE',1) 
 79:              RVAL=FLOAT(VALENCE) 
 80:              GKAPPA= 5.622667*RVAL*SQRT(CONC/TEMPR)     !1/debye length 
 81:           ENDIF 
 82:           IF (LGOUY) THEN 
 83:              ANFR=0.3 
 84:              ANFR= GTRMF(COMLYN,COMLEN,'ANFR',ANFR) 
 85:              AREA=70.0 
 86:              AREA= GTRMF(COMLYN,COMLEN,'AREA',AREA) 
 87:              OFFST= GTRMF(COMLYN,COMLEN,'OFFS',THREE) 
 88:              PSI0=-2334.2*ANFR/AREA/SQRT(TEMPR*CONC) 
 89:              GALPHA= LOG(PSI0+SQRT(PSI0**2+1))/RVAL 
 90:              PSI0=1.7235D-4*TEMPR/RVAL*LOG(PSI0+SQRT(PSI0**2+1)) 
 91:              IF (PRNLEV .GT. 9) WRITE (6,*) 'Psi0 (V) ', PSI0 
 92:              GALPHA= (EXP(GALPHA)-1.0)/(EXP(GALPHA)+1.0) 
 93:           ENDIF 
 94:           IF (LVOLT) THEN 
 95:              ACONS= VOLT/(2.+40.*GKAPPA*WIDTH) 
 96:           ENDIF 
 97:           ENDIF   !LMEMBR 
 98:           COMLYN="READ CARD "//COMLYN 55:           COMLYN="READ CARD "//COMLYN
 99:           COMLEN=COMLEN+10 56:           COMLEN=COMLEN+10
100:           CALL OPNLGU(COMLYN,COMLEN,UNPAR) 57:           CALL OPNLGU(COMLYN,COMLEN,UNPAR)
101:           CALL XTRANE(COMLYN,COMLEN,'EEF1') 
102:           CALL RDSLVPAR(ATC,NATC,SLVITC,VOLM,GREF,MAXATC,COMLYN, 58:           CALL RDSLVPAR(ATC,NATC,SLVITC,VOLM,GREF,MAXATC,COMLYN,
103:      &    COMLEN,MXCMSZ,GFREE,UNPAR,TEMPR,SIGW,HREF,CPREF,SLVNT) 59:      &    COMLEN,MXCMSZ,GFREE,UNPAR,TEMPR,SIGW,HREF,CPREF)
104:           DO 10 I=1,NATC 60:           DO 10 I=1,NATC
105:           FSIGW(I)= TWOPI*SQRT(PI)*SIGW(I) 61:           FSIGW(I)= TWOPI*SQRT(PI)*SIGW(I)
106:           IF (SIGW(i).GT.RSMALL) THEN 62:           IF (SIGW(i).GT.RSMALL) THEN
107:           FSIGW(I)= ONE/FSIGW(I) 63:           FSIGW(I)= ONE/FSIGW(I)
108:           SIGW(I)= ONE/SIGW(I) 64:           SIGW(I)= ONE/SIGW(I)
109:           ENDIF 65:           ENDIF
110: 10        CONTINUE 66: 10        CONTINUE
111:           DO 20 I=1,NATOM 67:           DO 20 I=1,NATOM
112:           GREFI(I)= GREF(SLVITC(IAC(I))) 68:           GREFI(I)= GREF(SLVITC(IAC(I)))
113:           HREFI(I)= HREF(SLVITC(IAC(I))) 
114:           CPREFI(I)= CPREF(SLVITC(IAC(I))) 
115:           GFREEI(I)= GFREE(SLVITC(IAC(I))) 69:           GFREEI(I)= GFREE(SLVITC(IAC(I)))
116:           VOLMI(I)= VOLM(SLVITC(IAC(I))) 70:           VOLMI(I)= VOLM(SLVITC(IAC(I)))
117:           SIGWI(I)= SIGW(SLVITC(IAC(I))) 71:           SIGWI(I)= SIGW(SLVITC(IAC(I)))
118:           FSIGWI(I)= FSIGW(SLVITC(IAC(I))) 72:           FSIGWI(I)= FSIGW(SLVITC(IAC(I)))
119:           VDWRI(I)= VDWR(ITC(IAC(I))) 73:           VDWRI(I)= VDWR(ITC(IAC(I)))
120: 20        CONTINUE 74: 20        CONTINUE
121:           IF (LMEMBR) THEN 
122:           REWIND (UNIT=UNPAR) 
123:           CALL RDSLVPAR(ATC,NATC,SLVITC,VOLM,GREF,MAXATC,COMLYN, 
124:      &    COMLEN,MXCMSZ,GFREE,UNPAR,TEMPR,SIGW,HREF,CPREF,SLVNT2) 
125:           DO 110 I=1,NATC 
126:           FSIGW(I)= 2.d+0*3.14159d+0*SQRT(3.14159d+0)*SIGW(I) 
127:           IF (SIGW(I).GT.1.d-10) THEN 
128:           FSIGW(I)= 1.d0/FSIGW(I) 
129:           SIGW(I)= 1.d0/SIGW(I) 
130:           ENDIF 
131: 110       CONTINUE 
132:           DO 120 I=1,NATOM 
133:           GREFI2(I)= GREF(SLVITC(IAC(I))) 
134:           GFREEI2(I)= GFREE(SLVITC(IAC(I))) 
135:           HREFI2(I)= HREF(SLVITC(IAC(I))) 
136:           CPREFI2(I)= CPREF(SLVITC(IAC(I))) 
137: c          VOLMI(I)= VOLM(SLVITC(IAC(I)))        Volume,vdwr the same 
138: c          VDWRI(I)= VDWR(ITC(IAC(I))) 
139:           SIGWI2(I)= SIGW(SLVITC(IAC(I))) 
140:           FSIGWI2(I)= FSIGW(SLVITC(IAC(I))) 
141: 120       CONTINUE 
142:           ENDIF 
143:       ENDIF 75:       ENDIF
144:       RETURN 76:       RETURN
145:       END 77:       END
146: C 
147:       SUBROUTINE EEF1EN(EU,X,Y,Z,DX,DY,DZ,QECONT,ECONT,IFRSTA,NATOMX, 
148:      &                    IFRSTG,NGRPX,JNBG,INBLOG,INB14,IBLO14,LINTE, 
149:      &                    DD1,IUPT,QSECD 
150:      -                    ,QBSECD,HSBLC,MAXBHES  !## VIBBLOCK 
151:      -                    ) 
152: C 
153: C     EU   - Energy to be returned 
154: C     X,Y,Z - Coordinates for energy evaluation 
155: C     DX,DY,DZ - Forces. 
156: C     QECONT - Flag for analysis (0-no analysis,>0=analysis) 
157: C     ECONT(NATOMX) - Analysis array to be filled if QECONT>0. 
158: C     NATOMX - Number of atoms 
159: C 
160: C     Author: T.Lazaridis 
161: C     Second derivatives by I.A. 
162: C 
163: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
164: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
165: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
166: C 
167: ##INCLUDE '~/charmm_fcm/psf.fcm' 
168: ##INCLUDE '~/charmm_fcm/param.fcm' 
169: ##INCLUDE '~/charmm_fcm/eef1.fcm' 
170: ##INCLUDE '~/charmm_fcm/stack.fcm' 
171: ##INCLUDE '~/charmm_fcm/number.fcm' 
172: ##INCLUDE '~/charmm_fcm/vibblock.fcm'  !##VIBBLOCK 
173: ##INCLUDE '~/charmm_fcm/parallel.fcm' 
174: ##INCLUDE '~/charmm_fcm/stream.fcm' 
175: ##INCLUDE '~/charmm_fcm/actclus.fcm' 
176: C ARD 00-10-14 
177: ##INCLUDE '~/charmm_fcm/block.fcm' 
178:       REAL*8 EU 
179:       REAL*8 X(*),Y(*),Z(*) 
180:       REAL*8 DX(*),DY(*),DZ(*) 
181:       LOGICAL QECONT 
182:       REAL*8 ECONT(*) 
183:       INTEGER IFRSTA,NATOMX,IFRSTG,NGRPX 
184:       INTEGER JNBG(*),INBLOG(*),INB14(*),IBLO14(*) 
185:       LOGICAL LINTE 
186:       REAL*8 DD1(*) 
187:       INTEGER IUPT(*) 
188:       LOGICAL QSECD 
189: ##IF VIBBLOCK 
190:       INTEGER MAXBHES 
191:       LOGICAL QBSECD 
192:       REAL*8 HSBLC(MAXBHES,*) 
193: ##ENDIF 
194: C 
195:       INTEGER OLDUSD,FEFF 
196: C 
197:       OLDUSD=LSTUSD 
198:       FEFF=ALLSTK(IREAL8(NATOMX)) 
199:       CALL EEF1EN2(EU,X,Y,Z,DX,DY,DZ,QECONT,ECONT,IFRSTA,NATOMX, 
200:      &                    IFRSTG,NGRPX,JNBG,INBLOG,INB14,IBLO14,LINTE, 
201:      &                    DD1,IUPT,QSECD,STACK(FEFF) 
202:      -                    ,QBSECD,HSBLC,MAXBHES  !## VIBBLOCK 
203:      -                    ) 
204: C 
205:       LSTUSD=OLDUSD 
206:       RETURN 
207:       END 
208:  78: 
209: C***************************************************************** 79: C*****************************************************************
210: C jdb I have modified this subroutine to implement a switching 80: C jdb I have modified this subroutine to implement a switching
211: C function for the energy and forces.  This switching is based 81: C function for the energy and forces.  This switching is based
212: C on interatomic separation rather than group separation.  Note 82: C on interatomic separation rather than group separation.  Note
213: C that GROUP must still be specified in the data file since the 83: C that GROUP must still be specified in the data file since the
214: C nonbond updates are group based.  This discrepency may cause 84: C nonbond updates are group based.  This discrepency may cause
215: C problems -- I haven't been able to figure it out yet. 85: C problems -- I haven't been able to figure it out yet.
216: C This switching routine helps get rid of the discontinuities. 86: C This switching routine helps get rid of the discontinuities.
217: C Another cause of discontinuities is the use of MYEXP rather 87: C Another cause of discontinuities is the use of MYEXP rather
218: C than just calculating the exponentials.  I have also replaced 88: C than just calculating the exponentials.  I have also replaced
219: C this. 89: C this.
220:       SUBROUTINE EEF1EN2(EU,X,Y,Z,DX,DY,DZ,QECONT,ECONT,IFRSTA,NATOMX, 90:         SUBROUTINE EEF1EN(EU,X,Y,Z,DX,DY,DZ,QECONT,ECONT,IFRSTA,NATOMX,
221:      &               IFRSTG,NGRPX,JNBGP,INBLOGP,INB14P,IBLO14P,LINTE, 91:      &                  IFRSTG,NGRPX,JNBGP,INBLOGP,INB14P,IBLO14P,LINTE,
222:      &               DD1,IUPT,QSECD,F 92:      &                    DD1,IUPT,QSECD
223:      -               ,QBSECD,HSBLC,MAXBHES  !## VIBBLOCK 93:      -                    ,QBSECD,HSBLC,MAXBHES  !## VIBBLOCK
224:      -               ) 94:      -                    )
225: C jdb I added the suffix P to the passed variables JNBG, INBLOG, INB14, and 95: C jdb I added the suffix P to the passed variables JNBG, INBLOG, INB14, and
226: C IBLO14 to distinguish these passed variables from the ones declared in 96: C IBLO14 to distinguish these passed variables from the ones declared in
227: C inbnd.fcm 97: C inbnd.fcm
228: C 98: C
229: C     EU   - Energy to be returned 99: C     EU   - Energy to be returned
230: C     X,Y,Z - Coordinates for energy evaluation100: C     X,Y,Z - Coordinates for energy evaluation
231: C     DX,DY,DZ - Forces.101: C     DX,DY,DZ - Forces.
232: C     QECONT - Flag for analysis (0-no analysis,>0=analysis)102: C     QECONT - Flag for analysis (0-no analysis,>0=analysis)
233: C     ECONT(NATOMX) - Analysis array to be filled if QECONT>0.103: C     ECONT(NATOMX) - Analysis array to be filled if QECONT>0.
234: C     NATOMX - Number of atoms104: C     NATOMX - Number of atoms
235: C105: C
236: C     Author: T.Lazaridis106: C     Author: T.Lazaridis
237: C     Second derivatives by I.A.107: C     Second derivatives by I.A.
238: C108: C
239: ##INCLUDE '~/charmm_fcm/impnon.fcm'109: ##INCLUDE '~/charmm_fcm/impnon.fcm'
240: ##INCLUDE '~/charmm_fcm/dimens.fcm'110: ##INCLUDE '~/charmm_fcm/dimens.fcm'
241: ##INCLUDE '~/charmm_fcm/number.fcm' 
242: C 
243: ##INCLUDE '~/charmm_fcm/psf.fcm'111: ##INCLUDE '~/charmm_fcm/psf.fcm'
244: ##INCLUDE '~/charmm_fcm/param.fcm'112: ##INCLUDE '~/charmm_fcm/param.fcm'
245: ##INCLUDE '~/charmm_fcm/eef1.fcm'113: ##INCLUDE '~/charmm_fcm/eef1.fcm'
246: C##INCLUDE '~/charmm_fcm/stack.fcm'   ! bs360: was not commented originally114: ##INCLUDE '~/charmm_fcm/stack.fcm'
 115: ##INCLUDE '~/charmm_fcm/number.fcm'
247: ##INCLUDE '~/charmm_fcm/vibblock.fcm'  !##VIBBLOCK116: ##INCLUDE '~/charmm_fcm/vibblock.fcm'  !##VIBBLOCK
248: ##INCLUDE '~/charmm_fcm/parallel.fcm'117: ##INCLUDE '~/charmm_fcm/parallel.fcm'
249: ##INCLUDE '~/charmm_fcm/stream.fcm'   !bs360: added 
250: C jdb added next include so nonbond cutoff values could be accessed118: C jdb added next include so nonbond cutoff values could be accessed
251: ##INCLUDE '~/charmm_fcm/inbnd.fcm'119: ##INCLUDE '~/charmm_fcm/inbnd.fcm'
252: ##INCLUDE '~/charmm_fcm/actclus.fcm'120: ##INCLUDE '~/charmm_fcm/actclus.fcm'
253: C ARD 00-10-14121: C ARD 00-10-14
254: ##INCLUDE '~/charmm_fcm/block.fcm'122: ##INCLUDE '~/charmm_fcm/block.fcm'
255:       REAL*8 EU123:       REAL*8 EU
 124:       INTEGER NATOMX,NGRPX,IFRSTA,IFRSTG
256:       REAL*8 X(*),Y(*),Z(*)125:       REAL*8 X(*),Y(*),Z(*)
257:       REAL*8 DX(*),DY(*),DZ(*)126:       REAL*8 DX(*),DY(*),DZ(*)
258:       LOGICAL QECONT127:       LOGICAL QECONT
259:       REAL*8 ECONT(*)128:       REAL*8 ECONT(*)
260:       INTEGER IFRSTA,NATOMX,IFRSTG,NGRPX129: C      INTEGER ATOMON(*) !RJP
261:       INTEGER JNBGP(*),INBLOGP(*),INB14P(*),IBLO14P(*) 
262:       LOGICAL LINTE 
263:       REAL*8 DD1(*) 
264:       INTEGER IUPT(*) 
265:       LOGICAL QSECD 
266:       REAL*8 F(*) 
267: ##IF VIBBLOCK 
268:       INTEGER MAXBHES 
269:       LOGICAL QBSECD 
270:       REAL*8 HSBLC(MAXBHES,*) 
271: ##ENDIF 
272: C 
273:       INTEGER I,J,K130:       INTEGER I,J,K
274:       REAL*8 EXRI,EXRJ,XDISTI,XDISTJ,RSQ,ARSQ,RDIST,ARDIST,ARDIST3131:       REAL*8 EXRI,EXRJ,XDISTI,XDISTJ,RSQ,ARSQ,RDIST,ARDIST,ARDIST3
275:       REAL*8 GRI,VI,GRJ,VJ132:       REAL*8 GRI,VI,GRJ,VJ
276:       REAL*8 VGIJ,VGJI,FDERI,FDERJ,GSI,GSJ133:       REAL*8 VGIJ,VGJI,FDERI,FDERJ,GSI,GSJ
277:       REAL*8 SIGI,SIGJ,FSIGI,FSIGJ134:       REAL*8 SIGI,SIGJ,FSIGI,FSIGJ
278:       REAL*8 TMP,DPSI135:       REAL*8 TMP
279:       INTEGER NB,ITEMP,NPR,IS,IQ,JSLC,JQLC,JRS,JRSPR,NXI,NXIMAX,JSX136:       INTEGER NB,ITEMP,NPR,IS,IQ,JSLC,JQLC,JRS,JRSPR,NXI,NXIMAX,JSX
280:       INTEGER IRS,INBX,ICNT137:       INTEGER IRS,INBX
281:       LOGICAL LEXCL138:       LOGICAL LEXCL,LINTE
282:       REAL*8 DXI,DYI,DZI,DXIT,DYIT,DZIT,ACC139:       REAL*8 DXI,DYI,DZI,DXIT,DYIT,DZIT
283:       CHARACTER*1 STR1140:       CHARACTER*1 STR1
284:       REAL*8 FF 
285:       REAL*8 EGOUY,TMPE,PSI,CGGOUY,EVOLT 
286: 141: 
287: C jmc for switching function.142: C jmc for switching function.
288:       REAL*8 S,DS,DDS,DRDXI,DRDYI,DRDZI,DRDXJ,DRDYJ,DRDZJ,ETERM143:       REAL*8 S,DS,DDS,DRDXI,DRDYI,DRDZI,DRDXJ,DRDYJ,DRDZJ,ETERM
 144:       INTEGER JNBGP(*),INBLOGP(*),INB14P(*),IBLO14P(*)
289:       REAL*8 DXITI,DYITI,DZITI,DXITJ,DYITJ,DZITJ,CUTON,CUTOFF145:       REAL*8 DXITI,DYITI,DZITI,DXITJ,DYITJ,DZITJ,CUTON,CUTOFF
290: 146: 
 147: ##IF VIBBLOCK
 148:       INTEGER MAXBHES
 149:       LOGICAL QBSECD
 150:       REAL*8 HSBLC(MAXBHES,*)
 151: ##ENDIF
 152: 
291: C Local variables for second derivatives. I.A.153: C Local variables for second derivatives. I.A.
292:       REAL*8 EFDERI,EFDERJ,FCTI,FCTJ,ARDIST6,ETMP154:       REAL*8 EFDERI,EFDERJ,FCTI,FCTJ,ARDIST6,ETMP
293:       REAL*8 AXX,AYY,AZZ,AXY,AXZ,AYZ155:       REAL*8 AXX,AYY,AZZ,AXY,AXZ,AYZ
294:       INTEGER IADD,II,JJ156:       INTEGER IADD,II,JJ
295:       LOGICAL LSECD157:       LOGICAL LSECD
 158:       REAL*8 DD1(*)
 159:       INTEGER IUPT(*)
 160:       LOGICAL QSECD
296: C jmc161: C jmc
297:       REAL*8 D2RDX2,D2RDY2,D2RDZ2,D2RDXY,D2RDXZ,D2RDYZ162:       REAL*8 D2RDX2,D2RDY2,D2RDZ2,D2RDXY,D2RDXZ,D2RDYZ
298: 163: 
299: C -----------------------------164: C -----------------------------
300:       LSECD=QSECD   !set option flag for 2nd derivative, I.A.165:       LSECD=QSECD   !set option flag for 2nd derivative, I.A.
301: C166: C
302: C 
303: C     QC: bug fix 06/27/07 
304: ##IF VIBBLOCK167: ##IF VIBBLOCK
305:       IF(QBSECD) THEN168:       IF(QBSECD) THEN
306:       LSECD=.false.169:       LSECD=QBSECD
 170:       ELSE
 171: ##ENDIF
 172:       LSECD=QSECD
 173: ##IF VIBBLOCK
307:       ENDIF174:       ENDIF
308: ##ENDIF175: ##ENDIF
309: 176: 
310:       EU = ZERO177:       EU = ZERO
311: 178: 
312: C jdb CUTOFF and CUTON are the shifting cutoffs used for EEF1.  They are 1.05 greater179: C jdb CUTOFF and CUTON are the shifting cutoffs used for EEF1.  They are 1.05 greater
313: C than CTOFNB and CTONNB since this gives good agreement with the energy ordering180: C than CTOFNB and CTONNB since this gives good agreement with the energy ordering
314: C from the original unshifted EEF1.  See the explanation in the docs files181: C from the original unshifted EEF1.  See the explanation in the docs files
315: C eef1gmin.doc for a more in depth explanation.182: C eef1gmin.doc for a more in depth explanation.
316: 183: 
317:       CUTON=CTONNB+1.05D0184:       CUTON=CTONNB+1.05D0
318:       CUTOFF=CTOFNB+1.05D0185:       CUTOFF=CTOFNB+1.05D0
319:       CTFSLV=CUTOFF*CUTOFF186:       CTFSLV=CUTOFF*CUTOFF
320: 187: 
321: C bs360: not clear 188:       IF(QECONT) THEN
322: C      IF(QECONT) THEN189:        DO 10 I=IFRSTA,NATOMX
323: C       DO 10 I=IFRSTA,NATOMX190:         ECONT(I)=ZERO
324: C        ECONT(I)=ZERO191: 10     CONTINUE
325: C10     CONTINUE192:       ENDIF
326: C      ENDIF 
327: C193: C
328: C Solvation calculation194: C Solvation calculation
329:       DO 20 I=IFRSTA,NATOMX195:       DO 20 I=IFRSTA,NATOMX
330:        GSOLV(I)= ZERO196:        GSOLV(I)= ZERO
331:        IF (LMEMBR) THEN 
332:             TMP = ABS(Z(I))*2.0/WIDTH 
333:             F(I)= TMP**NSMTH/(1.0+TMP**NSMTH) 
334:             LAM(I)= AEMPIR + (1.0-AEMPIR)*F(I) 
335: c            DFDZ(I)= (2.*Z(I)/ABS(Z(I))/WIDTH)* 
336:             DFDZ(I)= (2.*sign(ONE,Z(I))/WIDTH)* 
337:      &               FLOAT(NSMTH)*TMP**(NSMTH-1)/(1.0+TMP**NSMTH)**2 
338: c           The contribution to the z derivative due to self-energy 
339:             DZ(I)= DZ(I)+ DFDZ(I)*(GREFI(I)-GREFI2(I)) 
340:        ENDIF 
341: 20    CONTINUE197: 20    CONTINUE
342: C     LOOP OVER PAIRS OF GROUPS IN GROUPS LIST198: C     LOOP OVER PAIRS OF GROUPS IN GROUPS LIST
343: C199: C
344:       EGOUY=0.d0 
345:       EVOLT=0.d0 
346:       NB=0200:       NB=0
347:       ITEMP=0201:       ITEMP=0
348:       DO 160 IRS=IFRSTG,NGRPX202:       DO 160 IRS=IFRSTG,NGRPX
349:         NPR=INBLOGP(IRS)-ITEMP203:         NPR=INBLOGP(IRS)-ITEMP
350:         ITEMP=INBLOGP(IRS)204:         ITEMP=INBLOGP(IRS)
351:         IS=IGPBS(IRS)+1205:         IS=IGPBS(IRS)+1
352:         IQ=IGPBS(IRS+1)206:         IQ=IGPBS(IRS+1)
353: C Gouy-Chapman term 
354: C This calculation uses the full charge of ionic residues. They are 
355: C recognized only based on partial charge, so the partial charge 
356: C has to be unique in the topology file! 
357:          IF ((LGOUY.OR.LVOLT) .AND. .NOT. LINTE) THEN 
358:            DO 34 I=IS,IQ 
359:            IF (abs(CG(I)+0.9).LT.1d-6) THEN          !Lys or Nt 
360:               CGGOUY=+0.1d0 
361:            ELSE IF (abs(CG(I)+0.121).LT.1d-6) THEN   !ARG 
362:               CGGOUY=+0.379d0 
363:            ELSE IF (abs(CG(I)-0.45).LT.1d-6) THEN    !HSC/HSP,GLYP,PROP 
364:               CGGOUY=+0.95d0 
365:            ELSE IF (abs(CG(I)-1.0).LT.1d-6) THEN     !ASP,GLU,Ct 
366:               CGGOUY= 0.d0 
367:            ELSE 
368:               CGGOUY=CG(I) 
369:            ENDIF 
370:            IF (LGOUY) THEN 
371:              TMP=ABS(Z(I))-WIDTH/2.0-OFFST 
372:              IF (TMP .LE. 0.0) THEN 
373:                PSI=PSI0*23.06 
374:              ELSE 
375:                TMPE=EXP(-GKAPPA*TMP) 
376:                PSI=3.974d-3*TEMPR*LOG((1.+GALPHA*TMPE)/(1.-GALPHA*TMPE)) 
377:                DZ(I)=DZ(I)+CGGOUY*3.974d-3*TEMPR*(1.-GALPHA*TMPE)/ 
378:      & (1.+GALPHA*TMPE)*2.0*GALPHA/(1.-GALPHA*TMPE)**2*(-GKAPPA*TMPE* 
379:      &  SIGN(ONE,Z(I))) 
380:              ENDIF 
381:              EGOUY= EGOUY+PSI*CGGOUY 
382:            ENDIF 
383:            IF (LVOLT) THEN 
384:              IF (Z(I).LE.-WIDTH/2) THEN 
385:                 TMP=Z(I)+WIDTH/2 
386:                 TMPE=EXP(GKAPPA*TMP) 
387:                 PSI=ACONS*TMPE*23.06 
388:                 DPSI=PSI*GKAPPA 
389:              ELSE IF (Z(I).GE.WIDTH/2) THEN 
390:                 TMP=Z(I)-WIDTH/2 
391:                 TMPE=EXP(-GKAPPA*TMP) 
392:                 PSI=(VOLT-ACONS*TMPE)*23.06 
393:                 DPSI=GKAPPA*ACONS*TMPE*23.06 
394:              ELSE 
395:                 TMP=Z(I)+WIDTH/2 
396:                 PSI=ACONS*(40.*GKAPPA*TMP+1.0)*23.06 
397:                 DPSI=GKAPPA*ACONS*40.0*23.06 
398:              ENDIF 
399:              EVOLT=EVOLT+PSI*CGGOUY 
400:              DZ(I)=DZ(I)+DPSI*CGGOUY 
401:            ENDIF 
402: 34         CONTINUE 
403:          ENDIF 
404: C Gouy end 
405:         DO 150 JRSPR=1,NPR207:         DO 150 JRSPR=1,NPR
406:             NB=NB+1208:             NB=NB+1
407:             JRS=JNBGP(NB)209:             JRS=JNBGP(NB)
408:             LEXCL=(JRS.LT.0)210:             LEXCL=(JRS.LT.0)
409:             IF(LEXCL) JRS=-JRS211:             IF(LEXCL) JRS=-JRS
410:             JSLC=IGPBS(JRS)+1212:             JSLC=IGPBS(JRS)+1
411:             JQLC=IGPBS(JRS+1)213:             JQLC=IGPBS(JRS+1)
412: C214: C
413: C     PROCESS THIS GROUP PAIR215: C     PROCESS THIS GROUP PAIR
414: C216: C
424:                 ENDIF226:                 ENDIF
425:                 NXIMAX=IBLO14P(I)227:                 NXIMAX=IBLO14P(I)
426:                 IF (IS.EQ.JSLC) THEN228:                 IF (IS.EQ.JSLC) THEN
427:                   JSX= I+1229:                   JSX= I+1
428:                 ELSE230:                 ELSE
429:                   JSX= JSLC231:                   JSX= JSLC
430:                 ENDIF232:                 ENDIF
431:               ELSE233:               ELSE
432:                 JSX= JSLC234:                 JSX= JSLC
433:               ENDIF235:               ENDIF
 236:               GRI= GFREEI(I)
434:               VI= VOLMI(I)237:               VI= VOLMI(I)
435:               IF (LMEMBR) THEN238:               SIGI= SIGWI(I)
436:                 GRI= GFREEI2(I)+ F(I)*(GFREEI(I)-GFREEI2(I))239:               FSIGI= FSIGWI(I)
437:                 SIGI= 1./SIGWI2(I)+ F(I)*(1./SIGWI(I)-1./SIGWI2(I)) 
438:                 SIGI= 1./SIGI 
439:                 FSIGI= 1./FSIGWI2(I)+ F(I)*(1./FSIGWI(I)-1./FSIGWI2(I)) 
440:                 FSIGI= 1./FSIGI 
441:               ELSE 
442:                 GRI= GFREEI(I) 
443:                 SIGI= SIGWI(I) 
444:                 FSIGI= FSIGWI(I) 
445:               ENDIF 
446:               DO 80 J=JSX,JQLC240:               DO 80 J=JSX,JQLC
447:                 STR1= ATC(IAC(J))241:                 STR1= ATC(IAC(J))
448:                 IF (STR1.EQ.'H') GOTO 80242:                 IF (STR1.EQ.'H') GOTO 80
449:                 IF (LEXCL) THEN243:                 IF (LEXCL) THEN
450: C  CHECK ATOM EXCLUSION LIST FOR EXCLUSIONS244: C  CHECK ATOM EXCLUSION LIST FOR EXCLUSIONS
451:    40             IF(NXI.GT.NXIMAX) GOTO 60245:    40             IF(NXI.GT.NXIMAX) GOTO 60
452:                   IF (INB14P(NXI).LT.0) THEN246:                   IF (INB14P(NXI).LT.0) THEN
453:                     INBX=-INB14P(NXI)247:                     INBX=-INB14P(NXI)
454:                     IF(J.LE.INBX) GOTO 60248:                     IF(J.LE.INBX) GOTO 60
455:                   ELSE249:                   ELSE
473:                 ARSQ=ONE/RSQ267:                 ARSQ=ONE/RSQ
474: 268: 
475: C jdb Call the switching function.  S is returned as the value of the switching269: C jdb Call the switching function.  S is returned as the value of the switching
476: C function, and DS is the derivation with respect to R.  We use the nonbond270: C function, and DS is the derivation with respect to R.  We use the nonbond
477: C switching distances CUTOFF and CUTON271: C switching distances CUTOFF and CUTON
478:                 RDIST= SQRT(RSQ)272:                 RDIST= SQRT(RSQ)
479:                 CALL SWITCHFUNC(RDIST,CUTON,CUTOFF,S,DS,DDS)273:                 CALL SWITCHFUNC(RDIST,CUTON,CUTOFF,S,DS,DDS)
480: 274: 
481:                 ARDIST= ONE/RDIST275:                 ARDIST= ONE/RDIST
482:                 ARDIST3= TWO*ARDIST*ARDIST*ARDIST276:                 ARDIST3= TWO*ARDIST*ARDIST*ARDIST
483: C QC: bug fix 06/27/07277:                 IF(LSECD) ARDIST6= ARSQ*ARSQ*ARSQ !needed by 2nd deriv, I.A.
484: ##IF VIBBLOCK278:                 GRJ= GFREEI(J)
485: C     needed by 2nd deriv, I.A. 
486:                 IF(LSECD.OR.QBSECD) ARDIST6= ARSQ*ARSQ*ARSQ 
487: ##ELSE 
488:                 IF(LSECD) ARDIST6= ARSQ*ARSQ*ARSQ 
489: ##ENDIF 
490:                 VJ= VOLMI(J)279:                 VJ= VOLMI(J)
491:                 IF (LMEMBR) THEN280:                 SIGJ= SIGWI(J)
492:                    GRJ= GFREEI2(J)+ F(J)*(GFREEI(J)-GFREEI2(J))281:                 FSIGJ= FSIGWI(J)
493:                    SIGJ= 1./SIGWI2(J)+ F(J)*(1./SIGWI(J)-1./SIGWI2(J)) 
494:                    FSIGJ=1./FSIGWI2(J)+F(J)*(1./FSIGWI(J)-1./FSIGWI2(J)) 
495:                    SIGJ= 1./SIGJ 
496:                    FSIGJ= 1./FSIGJ 
497:                 ELSE 
498:                    GRJ= GFREEI(J) 
499:                    SIGJ= SIGWI(J) 
500:                    FSIGJ= FSIGWI(J) 
501:                 ENDIF 
502:                 XDISTI= ABS((RDIST-VDWRI(I))*SIGI) !reduced distance282:                 XDISTI= ABS((RDIST-VDWRI(I))*SIGI) !reduced distance
503:                 EXRI=DEXP(-XDISTI*XDISTI)283:                 EXRI=DEXP(-XDISTI*XDISTI)
504:                 XDISTJ= ABS((RDIST-VDWRI(J))*SIGJ) 284:                 XDISTJ= ABS((RDIST-VDWRI(J))*SIGJ) 
505:                 EXRJ=DEXP(-XDISTJ*XDISTJ)285:                 EXRJ=DEXP(-XDISTJ*XDISTJ)
506:                 VGJI= VJ*GRI286:                 VGJI= VJ*GRI
507:                 VGIJ= VI*GRJ287:                 VGIJ= VI*GRJ
508: 288: 
509: C ARD 00-10-14289: C ARD 00-10-14
510: C Moved EU sum into loop for Monte Carlo which uses only a partial290: C Moved EU sum into loop for Monte Carlo which uses only a partial
511: C INBLOG list.  Also put in NOFORC flag to skip force calculation.291: C INBLOG list.  Also put in NOFORC flag to skip force calculation.
512: C Note:  if IFRSTA>1 and NATOMX<NATOM not all GSOLV will be initialized.292: C Note:  if IFRSTA>1 and NATOMX<NATOM not all GSOLV will be initialized.
513: 293: 
514:                 GSI = VGJI*EXRI*ARSQ*FSIGI294:                 GSI = VGJI*EXRI*ARSQ*FSIGI
515:                 GSJ = VGIJ*EXRJ*ARSQ*FSIGJ295:                 GSJ = VGIJ*EXRJ*ARSQ*FSIGJ
516: C jdb multiply the solvation energy by the switching function.  Note that296: C jdb multiply the solvation energy by the switching function.  Note that
517: C this is applied termwise, and does not multiply the cumulative energy.297: C this is applied termwise, and does not multiply the cumulative energy.
518:                 GSOLV(I)= GSOLV(I) + S*(-GSI)298:                 GSOLV(I)= GSOLV(I) + S*(-GSI)
519:                 GSOLV(J)= GSOLV(J) + S*(-GSJ)299:                 GSOLV(J)= GSOLV(J) + S*(-GSJ)
520: C jmc                EU = EU - GSI - GSJ300: C jmc                EU = EU - GSI - GSJ
521: C We need to apply the switching function to EU here as well as to the individual 301: C Now (in c29) EU is cumulated here as opposed to summing over GSOLV at the end (c27),
 302: C so need to apply the switching function to EU here as well as to the individual 
522: C components GSOLV!!! 303: C components GSOLV!!! 
523: C EU is now switched energy.304: C EU is now switched energy.
524:                 EU=EU-S*(GSI+GSJ)305:                 EU=EU-S*(GSI+GSJ)
525: 306: 
526:                 IF (.NOT. NOFORC) THEN                         !##BLOCK307:                 IF (.NOT. NOFORC) THEN                         !##BLOCK
527: 308: 
528:                 FDERI= EXRI*(XDISTI*SIGI+ARDIST)*ARDIST3*FSIGI309:                 FDERI= EXRI*(XDISTI*SIGI+ARDIST)*ARDIST3*FSIGI
529:                 FDERJ= EXRJ*(XDISTJ*SIGJ+ARDIST)*ARDIST3*FSIGJ310:                 FDERJ= EXRJ*(XDISTJ*SIGJ+ARDIST)*ARDIST3*FSIGJ
530:                 TMP= FDERI*VGJI+FDERJ*VGIJ311:                 TMP= FDERI*VGJI+FDERJ*VGIJ
531: 312: 
569: C jmc original code350: C jmc original code
570: c               DXIT= DXI*TMP351: c               DXIT= DXI*TMP
571: c               DYIT= DYI*TMP352: c               DYIT= DYI*TMP
572: C               DZIT= DZI*TMP353: C               DZIT= DZI*TMP
573: C               DX(I)= DX(I)+ DXIT354: C               DX(I)= DX(I)+ DXIT
574: C               DY(I)= DY(I)+ DYIT355: C               DY(I)= DY(I)+ DYIT
575: C               DZ(I)= DZ(I)+ DZIT356: C               DZ(I)= DZ(I)+ DZIT
576: c               DX(J)= DX(J)- DXIT357: c               DX(J)= DX(J)- DXIT
577: C               DY(J)= DY(J)- DYIT358: C               DY(J)= DY(J)- DYIT
578: C               DZ(J)= DZ(J)- DZIT359: C               DZ(J)= DZ(J)- DZIT
579: c The contributions to z derivatives360: 
580:                 IF (LMEMBR) THEN 
581:                    DZ(I)= DZ(I) - S * ( VJ*ARSQ*FSIGI*EXRI*DFDZ(I)* 
582:      &             (GFREEI(I)-GFREEI2(I)+ GRI*SIGI*(2.*XDISTI**2 -1.0)* 
583:      &             (1./SIGWI(I)-1./SIGWI2(I))) ) 
584:                    DZ(J)= DZ(J) - S * ( VI*ARSQ*FSIGJ*EXRJ*DFDZ(J)* 
585:      &             (GFREEI(J)-GFREEI2(J)+ GRJ*SIGJ*(2.*XDISTJ**2 -1.0)* 
586:      &             (1./SIGWI(J)-1./SIGWI2(J))) ) 
587:                 ENDIF 
588: C 
589: C jmc for 2nd derivs361: C jmc for 2nd derivs
590: C signs as for d2R/dXidXi.  Reverse sign for d2R/dXidXj etc.362: C signs as for d2R/dXidXi.  Reverse sign for d2R/dXidXj etc.
591:                 D2RDX2=-DXI*DXI/(RSQ*RDIST)+ARDIST363:                 D2RDX2=-DXI*DXI/(RSQ*RDIST)+ARDIST
592:                 D2RDY2=-DYI*DYI/(RSQ*RDIST)+ARDIST364:                 D2RDY2=-DYI*DYI/(RSQ*RDIST)+ARDIST
593:                 D2RDZ2=-DZI*DZI/(RSQ*RDIST)+ARDIST365:                 D2RDZ2=-DZI*DZI/(RSQ*RDIST)+ARDIST
594: C jmc D2RDXY = d2R/dxidyi = d2R/dxjdyj, similarly for xz and yz.366: C jmc D2RDXY = d2R/dxidyi = d2R/dxjdyj, similarly for xz and yz.
595: C Reverse sign for d2r/dxidyj  etc.367: C Reverse sign for d2r/dxidyj  etc.
596:                 D2RDXY=-DXI*DYI/(RSQ*RDIST)368:                 D2RDXY=-DXI*DYI/(RSQ*RDIST)
597:                 D2RDXZ=-DXI*DZI/(RSQ*RDIST)369:                 D2RDXZ=-DXI*DZI/(RSQ*RDIST)
598:                 D2RDYZ=-DYI*DZI/(RSQ*RDIST)370:                 D2RDYZ=-DYI*DZI/(RSQ*RDIST)
807: C##ENDIF  (lsecd_main1)579: C##ENDIF  (lsecd_main1)
808:    80         CONTINUE580:    80         CONTINUE
809:    90       CONTINUE581:    90       CONTINUE
810:  150    CONTINUE582:  150    CONTINUE
811: 160   CONTINUE583: 160   CONTINUE
812: 584: 
813: C AvdV 12/5/01585: C AvdV 12/5/01
814: ##IF PARALLEL586: ##IF PARALLEL
815: ##IFN VIBRAN587: ##IFN VIBRAN
816:       call gcomb(gsolv,natomx)588:       call gcomb(gsolv,natomx)
 589: C Communication of the atomon (integer) array can be done with
 590: C the bitwise or routine, instead of adapting gcomb for integer routines.
 591: C      call gbor(atomon,natomx)
817: ##ENDIF592: ##ENDIF
818: ##ENDIF593: ##ENDIF
819: 594: 
820: C ARD 00-10-14595: C ARD 00-10-14
821: C Moved addition of GREFI contribution down here due to EU move to loop.596: C Moved addition of GREFI contribution down here due to EU move to loop.
822: C AvdV 12/05/01597: C AvdV 12/05/01
823: C Add the GREFI only to one node for parallel runs to avoid double counting.598: C Add the GREFI only to one node for parallel runs to avoid double counting.
824: ##IF PARALLEL599: ##IF PARALLEL
825: ##IFN VIBRAN600: ##IFN VIBRAN
826:       if (mynod.eq.0) then601:       if (mynod.eq.0) then
827: ##ENDIF602: ##ENDIF
828: ##ENDIF603: ##ENDIF
829:       IF (.NOT. LINTE) THEN604:       IF (.NOT. LINTE) THEN
830:         DO 208 I=IFRSTA,NATOMX605:         DO 208 I=IFRSTA,NATOMX
831:           IF (LMEMBR) THEN606:           GSOLV(I) = GSOLV(I) + GREFI(I)
832:             GSOLV(I) = GSOLV(I) + GREFI2(I)+ F(I)*(GREFI(I)-GREFI2(I))607:           IF (ACAFLG(I).EQ.1) EU = EU + GREFI(I)  !if atom is active
833:           ELSE 
834:             GSOLV(I) = GSOLV(I) + GREFI(I) 
835:           ENDIF 
836: Cbs360          IF(QLBYCC) THEN 
837:           IF(LBYCC) THEN 
838:             IF (ACAFLG(I).EQ.1) THEN  !if atom is active 
839:               IF (LMEMBR) THEN 
840:                   EU = EU + GREFI2(I)+ F(I)*(GREFI(I)-GREFI2(I)) 
841:               ELSE 
842:                   EU = EU + GREFI(I) 
843:               ENDIF 
844:             ENDIF 
845:           ELSE 
846:             IF (LMEMBR) THEN 
847:                EU = EU + GREFI2(I)+ F(I)*(GREFI(I)-GREFI2(I)) 
848:             ELSE 
849:                EU = EU + GREFI(I) 
850:             ENDIF 
851:           ENDIF 
852: 208     CONTINUE608: 208     CONTINUE
853:       ENDIF609:       ENDIF
854: ##IF PARALLEL610: ##IF PARALLEL
855: ##IFN VIBRAN611: ##IFN VIBRAN
856:       endif612:       endif
857: ##ENDIF613: ##ENDIF
858: ##ENDIF614: ##ENDIF
859: 615: 
860: C AvdV 12/05/01616: C AvdV 12/05/01
861: C The following code is fine for both parallel and serial617: C The following code is fine for both parallel and serial
862:       IF (QECONT) THEN618:       IF (QECONT) THEN
863:          DO 209 I=IFRSTA,NATOMX619:           DO 209 I=IFRSTA,NATOMX
864: C added if test--RJP620: C added if test--RJP
865: C bs360           IF(QLBYCC) THEN621:            IF (ACAFLG(I).EQ.1) THEN !if atom is active
866:            IF (LBYCC) THEN622:             ECONT(I)=GSOLV(I)
867:               IF (ACAFLG(I).EQ.1) ECONT(I)=GSOLV(I) !if atom is active 
868:            ELSE 
869:               ECONT(I)=ECONT(I)+GSOLV(I) 
870:            ENDIF623:            ENDIF
871: 209      CONTINUE624: 209       CONTINUE
872:       ENDIF 
873:       IF (LGOUY) THEN 
874:         EU=EU+EGOUY 
875:         IF (PRNLEV .GT. 9) WRITE (6,*) 'Egouy is ', EGOUY 
876:       ENDIF625:       ENDIF
877:       IF (LVOLT) THEN626: 
878:         EU=EU+EVOLT 
879:         IF (PRNLEV .GT. 9) WRITE (6,*) 'Evolt is ', EVOLT 
880:       ENDIF 
881: C 
882:       RETURN627:       RETURN
883:       END628:       END
884: 629: 
885:         SUBROUTINE SLVPRINT(COMLYN,COMLEN)630:         SUBROUTINE SLVPRINT(COMLYN,COMLEN)
886: C Calls EEF1EN and then prints out the solvation energies of each atom631: C Calls EEF1EN and then prints out the solvation energies of each atom
887: C632: C
888: ##INCLUDE '~/charmm_fcm/impnon.fcm'633: ##INCLUDE '~/charmm_fcm/impnon.fcm'
889: ##INCLUDE '~/charmm_fcm/dimens.fcm'634: ##INCLUDE '~/charmm_fcm/dimens.fcm'
890: ##INCLUDE '~/charmm_fcm/psf.fcm'635: ##INCLUDE '~/charmm_fcm/psf.fcm'
891: ##INCLUDE '~/charmm_fcm/coord.fcm'636: ##INCLUDE '~/charmm_fcm/coord.fcm'
894: ##INCLUDE '~/charmm_fcm/param.fcm'639: ##INCLUDE '~/charmm_fcm/param.fcm'
895: ##INCLUDE '~/charmm_fcm/econt.fcm'640: ##INCLUDE '~/charmm_fcm/econt.fcm'
896: ##INCLUDE '~/charmm_fcm/exfunc.fcm'641: ##INCLUDE '~/charmm_fcm/exfunc.fcm'
897: ##INCLUDE '~/charmm_fcm/heap.fcm'642: ##INCLUDE '~/charmm_fcm/heap.fcm'
898: ##INCLUDE '~/charmm_fcm/bases.fcm'643: ##INCLUDE '~/charmm_fcm/bases.fcm'
899: ##INCLUDE '~/charmm_fcm/inbnd.fcm'644: ##INCLUDE '~/charmm_fcm/inbnd.fcm'
900: ##INCLUDE '~/charmm_fcm/eef1.fcm'645: ##INCLUDE '~/charmm_fcm/eef1.fcm'
901: ##INCLUDE '~/charmm_fcm/parallel.fcm'646: ##INCLUDE '~/charmm_fcm/parallel.fcm'
902: ##INCLUDE '~/charmm_fcm/stack.fcm' !RJP647: ##INCLUDE '~/charmm_fcm/stack.fcm' !RJP
903: C648: C
904:         REAL*8 EU,TOTALI,TOTARO,TOTPOL,TMP,F,HRF,CPRF649:         REAL*8 EU,TOTALI,TOTARO,TOTPOL
905:         REAL*8 HSOLV,CPSOLV,HALI,HARO,HPOL,CPALI,CPARO,CPPOL,GRF650:         REAL*8 HSOLV,CPSOLV,HALI,HARO,HPOL,CPALI,CPARO,CPPOL,GRF
906:         CHARACTER*(*) COMLYN651:         CHARACTER*(*) COMLYN
907:         INTEGER ISEG,IRES,I,COMLEN652:         INTEGER ISEG,IRES,I,COMLEN
 653: C        INTEGER ATOMON,OLDLST  !RJP
908:         CHARACTER*2 STR2654:         CHARACTER*2 STR2
909: 655: C allocate temporary memory !RJP
 656: C        OLDLST=LSTUSD
 657: C        ATOMON=ALLSTK(INTEG4(NATOM))
910:         CALL EEF1EN(EU,X,Y,Z,DX,DY,DZ,.FALSE.,ECONT,1,NATOM,1,NGRP,658:         CALL EEF1EN(EU,X,Y,Z,DX,DY,DZ,.FALSE.,ECONT,1,NATOM,1,NGRP,
911:      &        HEAP(BNBND(JNBG)),HEAP(BNBND(INBLOG)),659:      &        HEAP(BNBND(JNBG)),HEAP(BNBND(INBLOG)),
912:      &        HEAP(BNBND(INB14)),HEAP(BNBND(IBLO14)),.FALSE.,660:      &        HEAP(BNBND(INB14)),HEAP(BNBND(IBLO14)),.FALSE.,
913:      &        0,0,.FALSE.661:      &        0,0,.FALSE.
914:      -               ,.false.,0,0  !## VIBBLOCK662:      -               ,.false.,0,0  !## VIBBLOCK
915:      -                )663:      -                )
 664: 
 665: C  free workspace
 666: C        CALL FRESTK(LSTUSD-OLDLST)
916: C667: C
917: ##IF PARALLEL668: ##IF PARALLEL
918: ##IFN VIBRAN669: ##IFN VIBRAN
919: 670: 
920:         CALL GCOMB(EU,1)671:         CALL GCOMB(EU,1)
921: C AvdV 12/5/01672: C AvdV 12/5/01
922: C Removed the gcomb statement: this is already done in eef1en673: C Removed the gcomb statement: this is already done in eef1en
923: ccc        CALL GCOMB(GSOLV,NATOM)674: ccc        CALL GCOMB(GSOLV,NATOM)
924:         IF (MYNOD.NE.0) RETURN675:         IF (MYNOD.NE.0) RETURN
925: ##ENDIF676: ##ENDIF
926: ##ENDIF677: ##ENDIF
927: C678: C
928:         IF (IOLEV.LE.0) RETURN 
929: C 
930:         WRITE (OUTU,100) EU679:         WRITE (OUTU,100) EU
931: 100        FORMAT ('SLVPRINT> TOTAL SOLVATION FREE ENERGY: ',F15.3,680: 100        FORMAT ('SLVPRINT> TOTAL SOLVATION FREE ENERGY: ',F15.3,
932:      &   ' Kcal/mol')681:      &   ' Kcal/mol')
933:         WRITE (OUTU,*) 'SLVPRINT>                            G         H682:         WRITE (OUTU,*) 'SLVPRINT>                            G         H
934:      &       CP'683:      &       CP'
935:         TOTALI= 0.0d+0684:         TOTALI= 0.0d+0
936:         TOTARO= 0.0d+0685:         TOTARO= 0.0d+0
937:         TOTPOL= 0.0d+0686:         TOTPOL= 0.0d+0
938:         HALI= 0.0d+0687:         HALI= 0.0d+0
939:         HARO= 0.0d+0688:         HARO= 0.0d+0
940:         HPOL= 0.0d+0689:         HPOL= 0.0d+0
941:         CPALI= 0.0d+0690:         CPALI= 0.0d+0
942:         CPARO= 0.0d+0691:         CPARO= 0.0d+0
943:         CPPOL= 0.0d+0692:         CPPOL= 0.0d+0
944:         DO 300 ISEG=1,NSEG693:         DO 300 ISEG=1,NSEG
945:         DO 200 IRES=NICTOT(ISEG)+1,NICTOT(ISEG+1)694:         DO 200 IRES=NICTOT(ISEG)+1,NICTOT(ISEG+1)
946:         DO 5 I=IBASE(IRES)+1,IBASE(IRES+1)695:         DO 5 I=IBASE(IRES)+1,IBASE(IRES+1)
947:         IF (AMASS(I).LT.1.1) GOTO 5696:         IF (AMASS(I).LT.1.1) GOTO 5
948:         IF (LMEMBR) THEN697:         GRF= GREF(SLVITC(IAC(I)))
949:             TMP = ABS(Z(I))*2.0/WIDTH 
950:             F= TMP**NSMTH/(1.0+TMP**NSMTH) 
951:             GRF = GREFI2(I)+ F*(GREFI(I)-GREFI2(I)) 
952:             HRF = HREFI2(I)+ F*(HREFI(I)-HREFI2(I)) 
953:             CPRF = CPREFI2(I)+ F*(CPREFI(I)-CPREFI2(I)) 
954:         ELSE 
955:             GRF= GREFI(I) 
956:             HRF= HREFI(I) 
957:             CPRF= CPREFI(I) 
958:         ENDIF 
959:         IF (DABS(GRF).GT.1.d-10) THEN698:         IF (DABS(GRF).GT.1.d-10) THEN
960:           HSOLV= HRF*GSOLV(I)/GRF699:           HSOLV= HREF(SLVITC(IAC(I)))*GSOLV(I)/GRF
961:           CPSOLV= CPRF*GSOLV(I)/GRF700:           CPSOLV= CPREF(SLVITC(IAC(I)))*GSOLV(I)/GRF
962:         ELSE701:         ELSE
963:           HSOLV= 0.d0702:           HSOLV= 0.d0
964:           CPSOLV= 0.d0703:           CPSOLV= 0.d0
965:         ENDIF704:         ENDIF
966:         WRITE(OUTU,'(I5,1X,A,1X,A,1X,A,1X,A,1X,A4,1X,3f11.3)')705:         WRITE(OUTU,'(I5,1X,A,1X,A,1X,A,1X,A,1X,A4,1X,3f11.3)')
967:      $       I,SEGID(ISEG)(1:idleng),RESID(IRES)(1:idleng),706:      $       I,SEGID(ISEG)(1:idleng),RESID(IRES)(1:idleng),
968:      $       RES(IRES)(1:idleng),TYPE(I)(1:idleng),ATC(IAC(I)),707:      $       RES(IRES)(1:idleng),TYPE(I)(1:idleng),ATC(IAC(I)),
969:      $       GSOLV(I),HSOLV,CPSOLV708:      $       GSOLV(I),HSOLV,CPSOLV
970:         STR2= ATC(IAC(I))                        !Want first character709:         STR2= ATC(IAC(I))                        !Want first character
971:         IF (STR2.EQ.'CH') THEN710:         IF (STR2.EQ.'CH') THEN
992: 103     FORMAT ('SLVPRINT> TOTAL POLAR     : ',F15.3,' Kcal/mol')731: 103     FORMAT ('SLVPRINT> TOTAL POLAR     : ',F15.3,' Kcal/mol')
993:         WRITE (OUTU,105) 'SLVPRINT> ENTHALPY (ALIPHATIC,AROMATIC,POLAR)'732:         WRITE (OUTU,105) 'SLVPRINT> ENTHALPY (ALIPHATIC,AROMATIC,POLAR)'
994:      &  ,HALI+HARO+HPOL,HALI,HARO,HPOL733:      &  ,HALI+HARO+HPOL,HALI,HARO,HPOL
995:         WRITE (OUTU,105) 'SLVPRINT>       CP (ALIPHATIC,AROMATIC,POLAR)'734:         WRITE (OUTU,105) 'SLVPRINT>       CP (ALIPHATIC,AROMATIC,POLAR)'
996:      &  ,CPALI+CPARO+CPPOL,CPALI,CPARO,CPPOL735:      &  ,CPALI+CPARO+CPPOL,CPALI,CPARO,CPPOL
997: 105     FORMAT (A,F15.3,' (',3F15.3,' )')736: 105     FORMAT (A,F15.3,' (',3F15.3,' )')
998:         RETURN737:         RETURN
999:         END738:         END
1000: 739: 
1001:         SUBROUTINE RDSLVPAR(ATC,NATC,SLVITC,VOLM,GREF,MAXATC,COMLYN,740:         SUBROUTINE RDSLVPAR(ATC,NATC,SLVITC,VOLM,GREF,MAXATC,COMLYN,
1002:      &     COMLEN,MXCMSZ,GFREE,UNPAR,TEMPR,SIGW,HREF,CPREF,SLVNT)741:      &     COMLEN,MXCMSZ,GFREE,UNPAR,TEMPR,SIGW,HREF,CPREF)
1003: C This reads the file solvpar.inp742: C This reads the file solvpar.inp
1004: C Modeled it after PARRDR743: C Modeled it after PARRDR
1005: C        GREF(i=1,atom types)= reference solvation free energy of an atom744: C        GREF(i=1,atom types)= reference solvation free energy of an atom
1006: C        GFREE(i=1,atom types)= free group solvation free energy745: C        GFREE(i=1,atom types)= free group solvation free energy
1007: C        VOLM(i=1,atom types)= volume of each atom746: C        VOLM(i=1,atom types)= volume of each atom
1008: C        SLVITC(atom type)= like ITC for NB terms, points to the correct index747: C        SLVITC(atom type)= like ITC for NB terms, points to the correct index
1009: C                            in GREF,VOLM, i.e., the Gref of atom J is748: C                            in GREF,VOLM, i.e., the Gref of atom J is
1010: C                                GREF(SLVITC(IAC(J)))749: C                                GREF(SLVITC(IAC(J)))
1011: C750: C
1012: ##INCLUDE '~/charmm_fcm/impnon.fcm'751: ##INCLUDE '~/charmm_fcm/impnon.fcm'
1018: C757: C
1019:         INTEGER NAT,MAXATC,NATC,MXCMSZ,COMLEN,I,J,UNPAR758:         INTEGER NAT,MAXATC,NATC,MXCMSZ,COMLEN,I,J,UNPAR
1020:         INTEGER SLVITC(*)759:         INTEGER SLVITC(*)
1021:         CHARACTER*(*) COMLYN760:         CHARACTER*(*) COMLYN
1022:         CHARACTER*(*) ATC(*)761:         CHARACTER*(*) ATC(*)
1023:         CHARACTER*4 WRD,AI762:         CHARACTER*4 WRD,AI
1024:         LOGICAL EOF763:         LOGICAL EOF
1025:         REAL*8 VOLM(MAXATC*2),GREF(MAXATC*2),GFREE(MAXATC*2)764:         REAL*8 VOLM(MAXATC*2),GREF(MAXATC*2),GFREE(MAXATC*2)
1026:         REAL*8 SIGW(MAXATC*2),HREF(MAXATC*2),CPREF(MAXATC*2)765:         REAL*8 SIGW(MAXATC*2),HREF(MAXATC*2),CPREF(MAXATC*2)
1027:         REAL*8 V,G,GF,H,CP,RATIO,TEMPR,SG766:         REAL*8 V,G,GF,H,CP,RATIO,TEMPR,SG
1028:         CHARACTER*4 SLVNT 
1029: C767: C
1030: ccc        IF(MYNOD.GT.0) GOTO 100      !##PARALLEL768:         IF(MYNOD.GT.0) GOTO 100      !##PARALLEL
1031: C769: C
1032:         NAT=0770:         NAT=0
1033: C Get input line771: C Get input line
1034:         EOF= .FALSE.772:         EOF= .FALSE.
1035: 5       CALL RDCMND(COMLYN,MXCMSZ,COMLEN,UNPAR,EOF,.TRUE.,773: 10      CALL RDCMND(COMLYN,MXCMSZ,COMLEN,UNPAR,EOF,.FALSE.,
1036:      &            .FALSE.,'SETUPSLV> ') 
1037:         CALL TRIME(COMLYN,COMLEN) 
1038:         IF (EOF) THEN      !only one solvent is there 
1039:            REWIND (UNIT=UNPAR) 
1040:            GOTO 10 
1041:         ENDIF 
1042:         WRD= NEXTA4(COMLYN,COMLEN) 
1043:         IF (WRD.NE.SLVNT) GOTO 5 
1044: C We have found the right solvent. Proceed with parameters 
1045: 10      EOF= .FALSE. 
1046:         CALL RDCMND(COMLYN,MXCMSZ,COMLEN,UNPAR,EOF,.TRUE., 
1047:      &            .FALSE.,'SETUPSLV> ')774:      &            .FALSE.,'SETUPSLV> ')
1048:         CALL TRIME(COMLYN,COMLEN)775:         CALL TRIME(COMLYN,COMLEN)
1049:         IF (EOF) GOTO 100776:         IF (EOF) GOTO 100
1050:         IF (.NOT.(COMLEN.GT.0)) GOTO 10777:         IF (.NOT.(COMLEN.GT.0)) GOTO 10
1051: C Data processing778: C Data processing
1052:         WRD= NEXTA4(COMLYN,COMLEN)779:         WRD= NEXTA4(COMLYN,COMLEN)
1053:         IF (WRD.EQ.'END') GOTO 100 
1054:         AI=WRD780:         AI=WRD
1055:         V=NEXTF(COMLYN,COMLEN)781:         V=NEXTF(COMLYN,COMLEN)
1056:         G=NEXTF(COMLYN,COMLEN)782:         G=NEXTF(COMLYN,COMLEN)
1057:         GF=NEXTF(COMLYN,COMLEN)783:         GF=NEXTF(COMLYN,COMLEN)
1058:         H=NEXTF(COMLYN,COMLEN)784:         H=NEXTF(COMLYN,COMLEN)
1059:         CP=NEXTF(COMLYN,COMLEN)785:         CP=NEXTF(COMLYN,COMLEN)
1060:         SG=NEXTF(COMLYN,COMLEN)786:         SG=NEXTF(COMLYN,COMLEN)
1061:         IF (DABS(G).LT.RSMALL) THEN787:         IF (DABS(G).LT.RSMALL) THEN
1062:           RATIO= ZERO788:           RATIO= ZERO
1063:         ELSE789:         ELSE
1119: C845: C
1120:       REAL*8 EIMSLV,ESLV846:       REAL*8 EIMSLV,ESLV
1121:       INTEGER BNBND(*),BIMAG(*),BIMAGX(*)847:       INTEGER BNBND(*),BIMAG(*),BIMAGX(*)
1122:       INTEGER IGPBS(*),IGPTYP(*),IAC(*)848:       INTEGER IGPBS(*),IGPTYP(*),IAC(*)
1123:       REAL*8 X(*),Y(*),Z(*),DX(*),DY(*),DZ(*)849:       REAL*8 X(*),Y(*),Z(*),DX(*),DY(*),DZ(*)
1124:       LOGICAL QECONT850:       LOGICAL QECONT
1125:       REAL*8 ECONT(*)851:       REAL*8 ECONT(*)
1126:       INTEGER NATOMX !RJP852:       INTEGER NATOMX !RJP
1127: C853: C
1128:       INTEGER I,IMX854:       INTEGER I,IMX
 855: C      INTEGER ATOMON,OLDLST !RJP
1129: C856: C
1130: ##INCLUDE '~/charmm_fcm/inbnd.fcm'857: ##INCLUDE '~/charmm_fcm/inbnd.fcm'
1131: ##INCLUDE '~/charmm_fcm/image.fcm'858: ##INCLUDE '~/charmm_fcm/image.fcm'
1132: ##INCLUDE '~/charmm_fcm/timer.fcm'859: ##INCLUDE '~/charmm_fcm/timer.fcm'
1133: ##INCLUDE '~/charmm_fcm/heap.fcm'860: ##INCLUDE '~/charmm_fcm/heap.fcm'
1134: ##INCLUDE '~/charmm_fcm/stream.fcm'861: ##INCLUDE '~/charmm_fcm/stream.fcm'
1135: ##INCLUDE '~/charmm_fcm/stack.fcm' !RJP862: ##INCLUDE '~/charmm_fcm/stack.fcm' !RJP
1136: C863: C
1137:       INTEGER BDUMMY(SNBNDT+1)864:       INTEGER BDUMMY(SNBNDT+1)
1138: C     Set up a dummy data structure for enbond.865: C     Set up a dummy data structure for enbond.
1166:         BDUMMY(INBLO)=BIMAGX(IMBLOS)893:         BDUMMY(INBLO)=BIMAGX(IMBLOS)
1167:         BDUMMY(JNB)=BIMAGX(IMJNBS)894:         BDUMMY(JNB)=BIMAGX(IMJNBS)
1168:         NNNB =HEAP(BIMAGX(NIMNBS))895:         NNNB =HEAP(BIMAGX(NIMNBS))
1169:         BDUMMY(INBLOG)=BIMAGX(IMBLOX)896:         BDUMMY(INBLOG)=BIMAGX(IMBLOX)
1170:         BDUMMY(JNBG)=BIMAGX(IMJNBX)897:         BDUMMY(JNBG)=BIMAGX(IMJNBX)
1171:         NNNBG=HEAP(BIMAGX(NIMNBX))898:         NNNBG=HEAP(BIMAGX(NIMNBX))
1172:         NNB14= 0899:         NNB14= 0
1173: C900: C
1174:         CALL SETBND(BDUMMY)901:         CALL SETBND(BDUMMY)
1175: C902: C
 903: C        OLDLST=LSTUSD
 904: C        ATOMON=ALLSTK(INTEG4(NATOMX))
1176:         CALL EEF1EN(ESLV,X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATIM,1,NIMGRP,905:         CALL EEF1EN(ESLV,X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATIM,1,NIMGRP,
1177:      &        HEAP(BDUMMY(JNBG)),HEAP(BDUMMY(INBLOG)),906:      &        HEAP(BDUMMY(JNBG)),HEAP(BDUMMY(INBLOG)),
1178:      &        HEAP(BDUMMY(INB14)),HEAP(BDUMMY(IBLO14)),.TRUE.,907:      &        HEAP(BDUMMY(INB14)),HEAP(BDUMMY(IBLO14)),.TRUE.,
1179:      &        0,0,.FALSE.908:      &        0,0,.FALSE.
1180:      -               ,.false.,0,0  !## VIBBLOCK909:      -               ,.false.,0,0  !## VIBBLOCK
1181:      -                )910:      -                )
1182: 911: 
 912: C        CALL FRESTK(LSTUSD-OLDLST)
1183:         ESLV=ESLV*HALF913:         ESLV=ESLV*HALF
1184:         DO I=1,IMX914:         DO I=1,IMX
1185:           DX(I)=DX(I)*HALF915:           DX(I)=DX(I)*HALF
1186:           DY(I)=DY(I)*HALF916:           DY(I)=DY(I)*HALF
1187:           DZ(I)=DZ(I)*HALF917:           DZ(I)=DZ(I)*HALF
1188:         ENDDO918:         ENDDO
1189:       ENDIF919:       ENDIF
1190: C920: C
1191: C     compute image nonbonded energies921: C     compute image nonbonded energies
1192: C922: C
1196: C926: C
1197:         BDUMMY(JNB)=BIMAGX(IMJNB)927:         BDUMMY(JNB)=BIMAGX(IMJNB)
1198:         BDUMMY(INBLO)=BIMAGX(IMBLO)928:         BDUMMY(INBLO)=BIMAGX(IMBLO)
1199:         NNNB =HEAP(BIMAGX(NIMNB))929:         NNNB =HEAP(BIMAGX(NIMNB))
1200:         BDUMMY(INBLOG)=BIMAGX(IMBLOG)930:         BDUMMY(INBLOG)=BIMAGX(IMBLOG)
1201:         BDUMMY(JNBG)=BIMAGX(IMJNBG)931:         BDUMMY(JNBG)=BIMAGX(IMJNBG)
1202:         NNNBG=HEAP(BIMAGX(NIMNBG))932:         NNNBG=HEAP(BIMAGX(NIMNBG))
1203: C933: C
1204:         CALL SETBND(BDUMMY)934:         CALL SETBND(BDUMMY)
1205: C935: C
 936: C        OLDLST=LSTUSD
 937: C        ATOMON=ALLSTK(INTEG4(NATOMX))
1206:         CALL EEF1EN(EIMSLV,X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATIM,1,NIMGRP,938:         CALL EEF1EN(EIMSLV,X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATIM,1,NIMGRP,
1207:      &        HEAP(BDUMMY(JNBG)),HEAP(BDUMMY(INBLOG)),939:      &        HEAP(BDUMMY(JNBG)),HEAP(BDUMMY(INBLOG)),
1208:      &        HEAP(BDUMMY(INB14)),HEAP(BDUMMY(IBLO14)),.TRUE.,940:      &        HEAP(BDUMMY(INB14)),HEAP(BDUMMY(IBLO14)),.TRUE.,
1209:      &        0,0,.FALSE.941:      &        0,0,.FALSE.
1210:      -               ,.false.,0,0  !## VIBBLOCK942:      -               ,.false.,0,0  !## VIBBLOCK
1211:      -                )943:      -                )
1212: 944: 
 945: C        CALL FRESTK(LSTUSD-OLDLST)
1213:         IF (TIMER.GT.1) THEN946:         IF (TIMER.GT.1) THEN
1214:           IF(PRNLEV.GE.2) WRITE(OUTU,'(1X,A)')947:           IF(PRNLEV.GE.2) WRITE(OUTU,'(1X,A)')
1215:      $      'IMAGE EEF1 TIMES:'948:      $      'IMAGE EEF1 TIMES:'
1216:           CALL TIMRE949:           CALL TIMRE
1217:           CALL TIMRB950:           CALL TIMRB
1218:         ENDIF951:         ENDIF
1219:       ENDIF952:       ENDIF
1220: C953: C
1221:       EIMSLV=EIMSLV+ESLV954:       EIMSLV=EIMSLV+ESLV
1222: C955: C


r16977/enbfast.src 2017-01-21 10:32:52.943188235 +0000 r16976/enbfast.src 2017-01-21 10:32:55.083188235 +0000
1146: ##IF LDM1146: ##IF LDM
1147: ##INCLUDE '~/charmm_fcm/lambda.fcm'1147: ##INCLUDE '~/charmm_fcm/lambda.fcm'
1148:       REAL*8 FALPHA1148:       REAL*8 FALPHA
1149:       INTEGER ILDM, JLDM1149:       INTEGER ILDM, JLDM
1150: ##ENDIF ! LDM1150: ##ENDIF ! LDM
1151: ##IF FLUCQ1151: ##IF FLUCQ
1152:       LOGICAL QFLUC1152:       LOGICAL QFLUC
1153:       REAL*8 FQCFOR(*)1153:       REAL*8 FQCFOR(*)
1154: ##ENDIF1154: ##ENDIF
1155: C1155: C
1156: ##IF ASPENER 
1157: ##INCLUDE '~/charmm_fcm/eef1.fcm' 
1158:       REAL*8 LAMD,IR2 
1159: ##ENDIF ! ASPENER 
1160: C 
1161:       REAL*8 ETEMP1,ETEMP2,EELPR,ENBPR1156:       REAL*8 ETEMP1,ETEMP2,EELPR,ENBPR
1162:       REAL*8 C2ONNB,C2OFNB,RUL3,RUL121157:       REAL*8 C2ONNB,C2OFNB,RUL3,RUL12
1163:       REAL*8 CA,CB,TR61158:       REAL*8 CA,CB,TR6
1164:       INTEGER IRS,IS,IQ,I,JRS,JS,JQ,J1159:       INTEGER IRS,IS,IQ,I,JRS,JS,JQ,J
1165:       INTEGER NAT,NB,ITEMP,NPR,JRSPR,NI,NJ1160:       INTEGER NAT,NB,ITEMP,NPR,JRSPR,NI,NJ
1166:       INTEGER LEX14,NXI,NXIMAX,JSX,INBX,IVECT,IACI1161:       INTEGER LEX14,NXI,NXIMAX,JSX,INBX,IVECT,IACI
1167:       REAL*8 DXI,DYI,DZI,DXIT,DYIT,DZIT,DXJT,DYJT,DZJT,DXIC,DYIC,DZIC1162:       REAL*8 DXI,DYI,DZI,DXIT,DYIT,DZIT,DXJT,DYJT,DZJT,DXIC,DYIC,DZIC
1168:       REAL*8 CGF,FUNCT,CGT,S,R2,DF,DFN,DFI,DFJ,DFF1163:       REAL*8 CGF,FUNCT,CGT,S,R2,DF,DFN,DFI,DFJ,DFF
1169:       REAL*8 DFVDW, DFELEC, DFZ1164:       REAL*8 DFVDW, DFELEC
1170:       REAL*8 SCENT,RIJL,RIJU1165:       REAL*8 SCENT,RIJL,RIJU
1171:       LOGICAL LST2,LEXCL,LSWITR,LVGRP,LSWIT,LCSW,LRSW1166:       LOGICAL LST2,LEXCL,LSWITR,LVGRP,LSWIT,LCSW,LRSW
1172: ##IF MC1167: ##IF MC
1173:       INTEGER LCENTR1168:       INTEGER LCENTR
1174: ##ENDIF1169: ##ENDIF
1175: ##IF MTS1170: ##IF MTS
1176: ##INCLUDE '~/charmm_fcm/tbmts.fcm'1171: ##INCLUDE '~/charmm_fcm/tbmts.fcm'
1177:       REAL*8 RR1,RR2,RR31172:       REAL*8 RR1,RR2,RR3
1178:       REAL*8 SWF,SWFE1173:       REAL*8 SWF,SWFE
1179: ##ENDIF1174: ##ENDIF
1371:                DYI=Y(I)-Y(J)1366:                DYI=Y(I)-Y(J)
1372:                DZI=Z(I)-Z(J)1367:                DZI=Z(I)-Z(J)
1373: CCC Changed by Georgios to prevent overflow of non-interacting blocks1368: CCC Changed by Georgios to prevent overflow of non-interacting blocks
1374: CCC that are on top of each other; Nov. 7 951369: CCC that are on top of each other; Nov. 7 95
1375:                S=MAX(RSMALL,DXI*DXI+DYI*DYI+DZI*DZI)1370:                S=MAX(RSMALL,DXI*DXI+DYI*DYI+DZI*DZI)
1376: CCC1371: CCC
1377: CCC  original code             S=DXI*DXI+DYI*DYI+DZI*DZI1372: CCC  original code             S=DXI*DXI+DYI*DYI+DZI*DZI
1378: CCC1373: CCC
1379:                R2=ONE/S1374:                R2=ONE/S
1380: C1375: C
1381: C               IF (LRSW) THEN 
1382: C                  EELPR=CGT*CGX(J)*R2 
1383: C                  DFELEC=R2*MINTWO*EELPR 
1384: C               ELSE 
1385: C                  EELPR=CGT*CGX(J)*SQRT(R2) 
1386: C                  DFELEC=-R2*EELPR 
1387: C               ENDIF 
1388: C               IF(LEX14.GT.0) THEN 
1389: C                  DFELEC=DFELEC*E14FAC 
1390: C                  EELPR=EELPR*E14FAC 
1391: C               ENDIF 
1392: C 
1393: ##IF ASPENER 
1394:                IF (LRSW .AND. LMEMBR) THEN 
1395:                   LAMD=SQRT(LAM(I)*LAM(J)) 
1396:                   IR2=SQRT(R2) 
1397:                   EELPR=CGT*CGX(J)*IR2*IR2**LAMD 
1398:                   DFELEC=-(ONE+LAMD)*R2*EELPR 
1399:                   DFZ= EELPR*LOG(IR2)*(ONE-AEMPIR)/TWO/LAMD 
1400:                ELSE IF (LRSW .AND. LDEBYE) THEN 
1401:                   IR2=SQRT(R2) 
1402:                   EELPR=CGT*CGX(J)*R2*EXP(-1./IR2/RDEBYE) 
1403:                   DFELEC=(R2*MINTWO-IR2/RDEBYE)*EELPR 
1404:                ELSE IF (LRSW) THEN 
1405: ##ELSE 
1406:                IF (LRSW) THEN1376:                IF (LRSW) THEN
1407: ##ENDIF !ASPENER 
1408:                   EELPR=CGT*CGX(J)*R21377:                   EELPR=CGT*CGX(J)*R2
1409:                   DFELEC=R2*MINTWO*EELPR1378:                   DFELEC=R2*MINTWO*EELPR
1410:                ELSE1379:                ELSE
1411:                   EELPR=CGT*CGX(J)*SQRT(R2)1380:                   EELPR=CGT*CGX(J)*SQRT(R2)
1412:                   DFELEC=-R2*EELPR1381:                   DFELEC=-R2*EELPR
1413:                ENDIF1382:                ENDIF
1414:                IF(LEX14.GT.0) THEN1383:                IF(LEX14.GT.0) THEN
1415:                   DFELEC=DFELEC*E14FAC1384:                   DFELEC=DFELEC*E14FAC
1416:                   EELPR=EELPR*E14FAC1385:                   EELPR=EELPR*E14FAC
1417:                   IF (LMEMBR) DFZ=DFZ*E14FAC   !##ASPENER 
1418:                ENDIF1386:                ENDIF
1419: 1387: C
1420:  
1421:                IF (LVGRP) THEN1388:                IF (LVGRP) THEN
1422:                   IVECT=LOWTP(MAX(IACNB(J),IACI))+IACNB(J)+IACI+LEX141389:                   IVECT=LOWTP(MAX(IACNB(J),IACI))+IACNB(J)+IACI+LEX14
1423:                   TR6=R2*R2*R21390:                   TR6=R2*R2*R2
1424:                   IF (.NOT.QETEN) THEN                      !##GOMODEL1391:                   IF (.NOT.QETEN) THEN                      !##GOMODEL
1425:                      CA=CCNBA(IVECT)*TR6*TR61392:                      CA=CCNBA(IVECT)*TR6*TR6
1426:                      CB=CCNBB(IVECT)*TR61393:                      CB=CCNBB(IVECT)*TR6
1427:                      ENBPR=CA-CB1394:                      ENBPR=CA-CB
1428:                      DFVDW=-SIX*R2*(ENBPR+CA)1395:                      DFVDW=-SIX*R2*(ENBPR+CA)
1429: ##IF GOMODEL1396: ##IF GOMODEL
1430:                   ELSE1397:                   ELSE
1610:                     DY(I)=DY(I)+DYIT*DOCFI1577:                     DY(I)=DY(I)+DYIT*DOCFI
1611:                     DZ(I)=DZ(I)+DZIT*DOCFI1578:                     DZ(I)=DZ(I)+DZIT*DOCFI
1612:                     DX(J)=DX(J)-DXIT*DOCFJ1579:                     DX(J)=DX(J)-DXIT*DOCFJ
1613:                     DY(J)=DY(J)-DYIT*DOCFJ1580:                     DY(J)=DY(J)-DYIT*DOCFJ
1614:                     DZ(J)=DZ(J)-DZIT*DOCFJ1581:                     DZ(J)=DZ(J)-DZIT*DOCFJ
1615:                   ELSE1582:                   ELSE
1616: ##ENDIF !dock1583: ##ENDIF !dock
1617: ##ENDIF (middle_block)1584: ##ENDIF (middle_block)
1618:                   DX(I)=DX(I)+DXIT1585:                   DX(I)=DX(I)+DXIT
1619:                   DY(I)=DY(I)+DYIT1586:                   DY(I)=DY(I)+DYIT
 1587:                   DZ(I)=DZ(I)+DZIT
1620:                   DX(J)=DX(J)-DXIT1588:                   DX(J)=DX(J)-DXIT
1621:                   DY(J)=DY(J)-DYIT1589:                   DY(J)=DY(J)-DYIT
1622: ##IF ASPENER1590:                   DZ(J)=DZ(J)-DZIT
1623:                   IF (LRSW .AND. LMEMBR) THEN 
1624:                     DZ(I)=DZ(I)+DZIT+DFZ*LAM(J)*DFDZ(I) 
1625:                     DZ(J)=DZ(J)-DZIT+DFZ*LAM(I)*DFDZ(J) 
1626:                   ELSE 
1627: ##ENDIF 
1628:                     DZ(I)=DZ(I)+DZIT 
1629:                     DZ(J)=DZ(J)-DZIT 
1630:                   ENDIF              !##ASPENER 
1631: C 
1632: ##IF BLOCK1591: ##IF BLOCK
1633: ##IF DOCK1592: ##IF DOCK
1634:                   ENDIF1593:                   ENDIF
1635: ##ENDIF1594: ##ENDIF
1636:                ENDIF1595:                ENDIF
1637: ##ENDIF1596: ##ENDIF
1638: C1597: C
1639: ##IF FLUCQ1598: ##IF FLUCQ
1640: C Add in electrostatic energy for this pair to FlucQ arrays1599: C Add in electrostatic energy for this pair to FlucQ arrays
1641:                IF (QFLUC) THEN1600:                IF (QFLUC) THEN
1731: C T.S.1690: C T.S.
1732:                   DFELEC=0.1691:                   DFELEC=0.
1733:                   DFVDW=0.1692:                   DFVDW=0.
1734: C1693: C
1735:                   DXI=X(I)-X(J)1694:                   DXI=X(I)-X(J)
1736:                   DYI=Y(I)-Y(J)1695:                   DYI=Y(I)-Y(J)
1737:                   DZI=Z(I)-Z(J)1696:                   DZI=Z(I)-Z(J)
1738:                   S=MAX(RSMALL,DXI*DXI+DYI*DYI+DZI*DZI)1697:                   S=MAX(RSMALL,DXI*DXI+DYI*DYI+DZI*DZI)
1739: CCCCC                  S=DXI*DXI+DYI*DYI+DZI*DZI1698: CCCCC                  S=DXI*DXI+DYI*DYI+DZI*DZI
1740:                   R2=ONE/S1699:                   R2=ONE/S
1741: ##IF ASPENER 
1742:                   IF (LRSW .AND. LMEMBR) THEN 
1743:                     LAMD=SQRT(LAM(I)*LAM(J)) 
1744:                     IR2=SQRT(R2) 
1745:                     EELPR=CGT*CGX(J)*IR2*IR2**LAMD 
1746:                     DFELEC=-(ONE+LAMD)*R2*EELPR 
1747:                     DFZ= EELPR*LOG(IR2)*(ONE-AEMPIR)/TWO/LAMD 
1748:                   ELSE IF (LRSW .AND. LDEBYE) THEN 
1749:                     IR2=SQRT(R2) 
1750:                     EELPR=CGT*CGX(J)*R2*EXP(-1./IR2/RDEBYE) 
1751:                     DFELEC=(R2*MINTWO-IR2/RDEBYE)*EELPR 
1752:                   ELSE IF (LRSW) THEN 
1753: ##ELSE 
1754:                   IF (LRSW) THEN1700:                   IF (LRSW) THEN
1755: ##ENDIF ! ASPENER 
1756:                      EELPR=CGT*CGX(J)*R21701:                      EELPR=CGT*CGX(J)*R2
1757:                      DFELEC=R2*MINTWO*EELPR1702:                      DFELEC=R2*MINTWO*EELPR
1758:                   ELSE1703:                   ELSE
1759:                      EELPR=CGT*CGX(J)*SQRT(R2)1704:                      EELPR=CGT*CGX(J)*SQRT(R2)
1760:                      DFELEC=-R2*EELPR1705:                      DFELEC=-R2*EELPR
1761:                   ENDIF1706:                   ENDIF
1762: C1707: C
1763: C                  IF (LRSW) THEN 
1764: C                     EELPR=CGT*CGX(J)*R2 
1765: C                     DFELEC=R2*MINTWO*EELPR 
1766: C                  ELSE 
1767: C                     EELPR=CGT*CGX(J)*SQRT(R2) 
1768: C                     DFELEC=-R2*EELPR 
1769: C                  ENDIF 
1770: C 
1771:                   IF(LVGRP) THEN1708:                   IF(LVGRP) THEN
1772:                      IVECT=LOWTP(MAX(IACNB(J),IACI))+IACNB(J)+IACI1709:                      IVECT=LOWTP(MAX(IACNB(J),IACI))+IACNB(J)+IACI
1773:                      TR6=R2*R2*R21710:                      TR6=R2*R2*R2
1774:                      IF (.NOT.QETEN) THEN                      !##GOMODEL1711:                      IF (.NOT.QETEN) THEN                      !##GOMODEL
1775:                         CA=CCNBA(IVECT)*TR6*TR61712:                         CA=CCNBA(IVECT)*TR6*TR6
1776:                         CB=CCNBB(IVECT)*TR61713:                         CB=CCNBB(IVECT)*TR6
1777:                         ENBPR=CA-CB1714:                         ENBPR=CA-CB
1778:                         DFVDW=-SIX*R2*(ENBPR+CA)1715:                         DFVDW=-SIX*R2*(ENBPR+CA)
1779: ##IF GOMODEL1716: ##IF GOMODEL
1780:                      ELSE1717:                      ELSE
1959:                     DY(I)=DY(I)+DYIT*DOCFI1896:                     DY(I)=DY(I)+DYIT*DOCFI
1960:                     DZ(I)=DZ(I)+DZIT*DOCFI1897:                     DZ(I)=DZ(I)+DZIT*DOCFI
1961:                     DX(J)=DX(J)-DXIT*DOCFJ1898:                     DX(J)=DX(J)-DXIT*DOCFJ
1962:                     DY(J)=DY(J)-DYIT*DOCFJ1899:                     DY(J)=DY(J)-DYIT*DOCFJ
1963:                     DZ(J)=DZ(J)-DZIT*DOCFJ1900:                     DZ(J)=DZ(J)-DZIT*DOCFJ
1964:                   ELSE1901:                   ELSE
1965: ##ENDIF !dock1902: ##ENDIF !dock
1966: ##ENDIF (bl_b)1903: ##ENDIF (bl_b)
1967:                   DX(I)=DX(I)+DXIT1904:                   DX(I)=DX(I)+DXIT
1968:                   DY(I)=DY(I)+DYIT1905:                   DY(I)=DY(I)+DYIT
 1906:                   DZ(I)=DZ(I)+DZIT
1969:                   DX(J)=DX(J)-DXIT1907:                   DX(J)=DX(J)-DXIT
1970:                   DY(J)=DY(J)-DYIT1908:                   DY(J)=DY(J)-DYIT
1971: ##IF ASPENER1909:                   DZ(J)=DZ(J)-DZIT
1972:                   IF (LRSW .AND. LMEMBR) THEN 
1973:                      IF (LSWITR) THEN 
1974:                   DZ(I)=DZ(I)+DZIT+DFZ*LAM(J)*DFDZ(I)*FUNCT 
1975:                   DZ(J)=DZ(J)-DZIT+DFZ*LAM(I)*DFDZ(J)*FUNCT 
1976:                      ELSE 
1977:                   DZ(I)=DZ(I)+DZIT+DFZ*LAM(J)*DFDZ(I) 
1978:                   DZ(J)=DZ(J)-DZIT+DFZ*LAM(I)*DFDZ(J) 
1979:                      ENDIF 
1980:                   ELSE 
1981: ##ENDIF 
1982:                     DZ(I)=DZ(I)+DZIT 
1983:                     DZ(J)=DZ(J)-DZIT 
1984:                   ENDIF              !##ASPENER 
1985: ##IF BLOCK1910: ##IF BLOCK
1986: ##IF DOCK1911: ##IF DOCK
1987:                    ENDIF1912:                    ENDIF
1988: ##ENDIF1913: ##ENDIF
1989:                   ENDIF1914:                   ENDIF
1990: ##ENDIF1915: ##ENDIF
1991: C1916: C
1992: ##IF FLUCQ1917: ##IF FLUCQ
1993: C If FlucQ is interactive, sum the interaction energies (scaled by the1918: C If FlucQ is interactive, sum the interaction energies (scaled by the
1994: C switching function if necessary) into the FlucQ arrays1919: C switching function if necessary) into the FlucQ arrays


r16977/energy.src 2017-01-21 10:32:51.159188235 +0000 r16976/energy.src 2017-01-21 10:32:53.551188235 +0000
158:       INTEGER I,IUPT,ATFRST,ATLAST158:       INTEGER I,IUPT,ATFRST,ATLAST
159:       REAL*8  ECSUM,ECHARMBF159:       REAL*8  ECSUM,ECHARMBF
160:       LOGICAL TBM2,TBM4,TBHY3,TBHY4160:       LOGICAL TBM2,TBM4,TBHY3,TBHY4
161:       LOGICAL QFARRAY,QDIM4161:       LOGICAL QFARRAY,QDIM4
162:       LOGICAL QOK162:       LOGICAL QOK
163:       LOGICAL QOKANGLE                       !## CFF163:       LOGICAL QOKANGLE                       !## CFF
164:       REAL*8 ct3,lrcn                        !## LRVDW164:       REAL*8 ct3,lrcn                        !## LRVDW
165: 165: 
166:       call timer_start(T_energy)             !## NEWTIMER166:       call timer_start(T_energy)             !## NEWTIMER
167:       call timer_start( T_inte)              !## NEWTIMER167:       call timer_start( T_inte)              !## NEWTIMER
168:   168: 
169: ##IF MTS (mts_setflags)169: ##IF MTS (mts_setflags)
170:       IF(.NOT.TBMTS) THEN170:       IF(.NOT.TBMTS) THEN
171: C        If MTS is not in use, then do all energy terms - BRB 01/06/98171: C        If MTS is not in use, then do all energy terms - BRB 01/06/98
172:          ENE1=.TRUE.  ! flag to do fast terms172:          ENE1=.TRUE.  ! flag to do fast terms
173:          ENE2=.TRUE.  ! flag to do medium terms173:          ENE2=.TRUE.  ! flag to do medium terms
174:          ENE3=.TRUE.  ! flag to do slow terms174:          ENE3=.TRUE.  ! flag to do slow terms
175:       ENDIF175:       ENDIF
176: ##ENDIF (mts_setflags)176: ##ENDIF (mts_setflags)
177: C177: C
178: C     Check if we have enough time and no interruptions178: C     Check if we have enough time and no interruptions
1220: ##ENDIF1220: ##ENDIF
1221: ##IF PARALLEL1221: ##IF PARALLEL
1222: C get energies1222: C get energies
1223:         TMERI(TIMINTE) = TMERI(TIMINTE) + ECLOCK()-TIMMER1223:         TMERI(TIMINTE) = TMERI(TIMINTE) + ECLOCK()-TIMMER
1224:         TIMMER = ECLOCK()1224:         TIMMER = ECLOCK()
1225: ##ENDIF1225: ##ENDIF
1226: 1226: 
1227:         call timer_stop(  T_ephi)                !##NEWTIMER1227:         call timer_stop(  T_ephi)                !##NEWTIMER
1228:         call timer_stpstrt( T_inte,T_nonbon)     !##NEWTIMER1228:         call timer_stpstrt( T_inte,T_nonbon)     !##NEWTIMER
1229: C1229: C
1230: C--------------------------------------------------------------------- 
1231: C ASPENER is moved here, before ENBOND because it calculates LAMD 
1232: C which is needed in the electrostatic calculation in IMM1 
1233: ##IF ASPENER 
1234:       IF(ENE3) THEN                                          !##MTS 
1235:         IF(SOLVENT_ENTERED.AND.QETERM(ASP)) THEN 
1236:           CALL ASPENR(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT) 
1237:           IF (TIMER.GT.1) CALL WRTTIM('ASP surface energy times:') 
1238:         ELSEIF (DOEEF1.AND.QETERM(ASP)) THEN 
1239:           CALL EEF1EN(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATOM, 
1240:      &        1,NGRP,HEAP(BNBND(JNBG)),HEAP(BNBND(INBLOG)), 
1241:      &        HEAP(BNBND(INB14)),HEAP(BNBND(IBLO14)),.FALSE., 
1242:      &        DD1,STACK(IUPT),QSECD 
1243:      -       ,.false.,0,0    !## VIBBLOCK 
1244:      -        ) 
1245:  
1246:           IF (TIMER.GT.1) CALL WRTTIM('EEF1 energy times:') 
1247:  
1248: ##IFN NOIMAGES 
1249:       if(.not.qboun) then   !##PBOUND 
1250:       if(.not.lbypbcb) then   !##PBCUBES 
1251:        IF(NTRANS.GT.0) THEN 
1252: C  added "NATOM" to argument list - RJP 
1253:           CALL EEF1IM(EIMSLV,BNBND,BIMAG,BIMAG, 
1254:      $      IGPBS,IGPTYP,IAC,DX,DY,DZ,X,Y,Z,QECONT,ECONT,NATOM) 
1255: C EIMSLV is the amount of solvation free energy of primary atoms 
1256: C excluded by image atoms 
1257:           ETERM(ASP)= ETERM(ASP)+ EIMSLV 
1258:        ENDIF 
1259:       endif                 !##PBCUBES 
1260:       endif                 !##PBOUND 
1261: ##ENDIF 
1262:         ENDIF 
1263:       ENDIF                                                  !##MTS 
1264: ##ENDIF ! ASPENER 
1265:  
1266: C 
1267: C================= NONBOND ENERGY TERMS ================================1230: C================= NONBOND ENERGY TERMS ================================
1268: C . Non-bond terms.1231: C . Non-bond terms.
1269:       IF (NATOM .GT. 0) THEN1232:       IF (NATOM .GT. 0) THEN
1270: C-----------------------------------------------------------------------1233: C-----------------------------------------------------------------------
1271: ##IF MTS1234: ##IF MTS
1272:        IF(TBMTS) THEN1235:        IF(TBMTS) THEN
1273:          IF(ENE1.AND.(TBHY1.AND.(.NOT.TBHY2))) TBHY3=.TRUE.1236:          IF(ENE1.AND.(TBHY1.AND.(.NOT.TBHY2))) TBHY3=.TRUE.
1274:          IF(ENE2) THEN1237:          IF(ENE2) THEN
1275:            IF((TBHY1.AND.TBHY2).OR.SLFG) TBHY3=.TRUE.1238:            IF((TBHY1.AND.TBHY2).OR.SLFG) TBHY3=.TRUE.
1276:          ENDIF1239:          ENDIF
1619:       IF (QFLUC) THEN1582:       IF (QFLUC) THEN
1620:          IF (QETERM(FQPOL)) THEN1583:          IF (QETERM(FQPOL)) THEN
1621:             CALL FQENER(ETERM(FQPOL),X,Y,Z,QETERM(FQPOL))1584:             CALL FQENER(ETERM(FQPOL),X,Y,Z,QETERM(FQPOL))
1622:          ENDIF1585:          ENDIF
1623:       ENDIF1586:       ENDIF
1624: ##ENDIF ! FLUCQ1587: ##ENDIF ! FLUCQ
1625: C-----------------------------------------------------------------------1588: C-----------------------------------------------------------------------
1626: ##IF OVERLAP1589: ##IF OVERLAP
1627:       IF(QETERM(QOVLAP).AND.QOLAP)CALL OLAPENER(ETERM(QOVLAP),DX,DY,DZ)1590:       IF(QETERM(QOVLAP).AND.QOLAP)CALL OLAPENER(ETERM(QOVLAP),DX,DY,DZ)
1628: ##ENDIF1591: ##ENDIF
1629: C1592: ##IF ASPENER
 1593:       IF(ENE3) THEN                                          !##MTS
 1594:         IF(SOLVENT_ENTERED.AND.QETERM(ASP)) THEN
 1595:           CALL ASPENR(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT)
 1596:           IF (TIMER.GT.1) CALL WRTTIM('ASP surface energy times:')
 1597:         ELSEIF (DOEEF1.AND.QETERM(ASP)) THEN
 1598:           CALL EEF1EN(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATOM,
 1599:      &        1,NGRP,HEAP(BNBND(JNBG)),HEAP(BNBND(INBLOG)),
 1600:      &        HEAP(BNBND(INB14)),HEAP(BNBND(IBLO14)),.FALSE.,
 1601:      &        DD1,STACK(IUPT),QSECD
 1602:      -       ,.false.,0,0    !## VIBBLOCK
 1603:      -        )
 1604: 
 1605:           IF (TIMER.GT.1) CALL WRTTIM('EEF1 energy times:')
 1606: 
 1607: ##IFN NOIMAGES
 1608:       if(.not.qboun) then   !##PBOUND
 1609:       if(.not.lbypbcb) then   !##PBCUBES
 1610:        IF(NTRANS.GT.0) THEN
 1611: C  added "NATOM" to argument list - RJP
 1612:           CALL EEF1IM(EIMSLV,BNBND,BIMAG,BIMAG,
 1613:      $      IGPBS,IGPTYP,IAC,DX,DY,DZ,X,Y,Z,QECONT,ECONT,NATOM)
 1614: C EIMSLV is the amount of solvation free energy of primary atoms
 1615: C excluded by image atoms
 1616:           ETERM(ASP)= ETERM(ASP)+ EIMSLV
 1617:        ENDIF
 1618:       endif                 !##PBCUBES
 1619:       endif                 !##PBOUND
 1620: ##ENDIF
 1621:         ENDIF
 1622:       ENDIF                                                  !##MTS
 1623: ##ENDIF ! ASPENER
1630: ##IF ASPMEMB1624: ##IF ASPMEMB
1631:       IF(ENE3) THEN                                          !##MTS1625:       IF(ENE3) THEN                                          !##MTS
1632:         IF(SOLVMEMB_ENTRD.AND.QETERM(ASP)) THEN1626:         IF(SOLVMEMB_ENTRD.AND.QETERM(ASP)) THEN
1633:           CALL ASPENERMB(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT)1627:           CALL ASPENERMB(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT)
1634:         ENDIF1628:         ENDIF
1635:       ENDIF                                                  !##MTS 1629:       ENDIF                                                  !##MTS 
1636: ##ENDIF ! ASPMEMB1630: ##ENDIF ! ASPMEMB
1637: C1631: C
1638: ##IF PARALLEL1632: ##IF PARALLEL
1639:       TMERI(TIMEXTE)=TMERI(TIMEXTE)+ECLOCK()-TIMMER1633:       TMERI(TIMEXTE)=TMERI(TIMEXTE)+ECLOCK()-TIMMER


r16977/eutil.src 2017-01-21 10:32:51.791188235 +0000 r16976/eutil.src 2017-01-21 10:32:53.831188235 +0000
 19: ##INCLUDE '~/charmm_fcm/energy.fcm' 19: ##INCLUDE '~/charmm_fcm/energy.fcm'
 20: C 20: C
 21:       REAL*8  X(*),Y(*),Z(*) 21:       REAL*8  X(*),Y(*),Z(*)
 22:       REAL*8  XREF(*),YREF(*),ZREF(*) 22:       REAL*8  XREF(*),YREF(*),ZREF(*)
 23:       INTEGER MASSW 23:       INTEGER MASSW
 24: C 24: C
 25:       INTEGER IX,IY,IZ,IXR,IYR,IZR 25:       INTEGER IX,IY,IZ,IXR,IYR,IZR
 26:       LOGICAL LMASS,QOK 26:       LOGICAL LMASS,QOK
 27:       INTEGER NSHKIT,OLDUSD 27:       INTEGER NSHKIT,OLDUSD
 28: C 28: C
 29:       INTEGER ia 
 30: C 
 31:       IF(QHOLO) THEN 29:       IF(QHOLO) THEN
 32:         LMASS=(MASSW.EQ.1) 30:         LMASS=(MASSW.EQ.1)
 33:         OLDUSD=LSTUSD 31:         OLDUSD=LSTUSD
 34:         IX=ALLSTK(IREAL8(NATOM)) 32:         IX=ALLSTK(IREAL8(NATOM))
 35:         IY=ALLSTK(IREAL8(NATOM)) 33:         IY=ALLSTK(IREAL8(NATOM))
 36:         IZ=ALLSTK(IREAL8(NATOM)) 34:         IZ=ALLSTK(IREAL8(NATOM))
 37:         CALL COPYR8(X,STACK(IX),NATOM) 35:         CALL COPYR8(X,STACK(IX),NATOM)
 38:         CALL COPYR8(Y,STACK(IY),NATOM) 36:         CALL COPYR8(Y,STACK(IY),NATOM)
 39:         CALL COPYR8(Z,STACK(IZ),NATOM) 37:         CALL COPYR8(Z,STACK(IZ),NATOM)
 40:         IXR=ALLSTK(IREAL8(NATOM)) 38:         IXR=ALLSTK(IREAL8(NATOM))


r16977/intere.src 2017-01-21 10:32:52.091188235 +0000 r16976/intere.src 2017-01-21 10:32:54.127188235 +0000
422:        CALL INTER3('Normal Atom Nonbond List',ISLCT,JSLCT,NATOM,NGRP,422:        CALL INTER3('Normal Atom Nonbond List',ISLCT,JSLCT,NATOM,NGRP,
423:      &            HEAP(BNBNDC(INBLO)),HEAP(BNBNDC(JNB)),423:      &            HEAP(BNBNDC(INBLO)),HEAP(BNBNDC(JNB)),
424:      &            HEAP(BNBND(INBLO)),HEAP(BNBND(JNB)),424:      &            HEAP(BNBND(INBLO)),HEAP(BNBND(JNB)),
425:      &            HEAP(BNBNDC(INBLOG))425:      &            HEAP(BNBNDC(INBLOG))
426:      $        ,lbycbim                        !##IMCUBES426:      $        ,lbycbim                        !##IMCUBES
427:      $        ,lbypbcb                        !##PBCUBES427:      $        ,lbypbcb                        !##PBCUBES
428:      $        )428:      $        )
429:       ENDIF429:       ENDIF
430:       EndIf  !##GRID430:       EndIf  !##GRID
431: C431: C
432: C EEF1EN added by Lazaridis, Apr 99. May2004: Moved before ENBOND432:       CALL ENBOND(ETERM(VDW),ETERM(ELEC),BNBNDC,
 433:      $       1,NATOM,CG,RSCLF,NGRP,IGPBS,IGPTYP,IAC,
 434:      $       DX,DY,DZ,X,Y,Z,QECONT,ECONT,.FALSE.,0,0,0,.FALSE.,
 435:      $       QETERM(VDW),QETERM(ELEC),
 436:      $       QFLUC,FQCFOR,  !##FLUCQ
 437:      $       QETERM(ST2),NST2,ETERM(ST2),.FALSE.
 438:      -               ,.false.,0,0  !## VIBBLOCK
 439:      -                )
 440: 
 441: C EEF1EN added by Lazaridis, Apr 99
433: ##IF ASPENER442: ##IF ASPENER
434:         IF (DOEEF1.AND.QETERM(ASP)) THEN443:         IF (DOEEF1.AND.QETERM(ASP)) THEN
435: c added last 3 arguments for 2nd derivatives. I.A.444: c added last 3 arguments for 2nd derivatives. I.A.
436:           CALL EEF1EN(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATOM,445:           CALL EEF1EN(ETERM(ASP),X,Y,Z,DX,DY,DZ,QECONT,ECONT,1,NATOM,
437:      &        1,NGRP,HEAP(BNBNDC(JNBG)),HEAP(BNBNDC(INBLOG)),446:      &        1,NGRP,HEAP(BNBNDC(JNBG)),HEAP(BNBNDC(INBLOG)),
438:      &        HEAP(BNBND(INB14)),HEAP(BNBND(IBLO14)),.TRUE.,447:      &        HEAP(BNBND(INB14)),HEAP(BNBND(IBLO14)),.TRUE.,
439:      &        0,0,.FALSE.448:      &        0,0,.FALSE.
440:      -               ,.false.,0,0  !## VIBBLOCK449:      -               ,.false.,0,0  !## VIBBLOCK
441:      -                )450:      -                )
442: 451: 
443:         ENDIF452:         ENDIF
444: ##ENDIF ! ASPENER453: ##ENDIF ! ASPENER
445: C 
446:       CALL ENBOND(ETERM(VDW),ETERM(ELEC),BNBNDC, 
447:      $       1,NATOM,CG,RSCLF,NGRP,IGPBS,IGPTYP,IAC, 
448:      $       DX,DY,DZ,X,Y,Z,QECONT,ECONT,.FALSE.,0,0,0,.FALSE., 
449:      $       QETERM(VDW),QETERM(ELEC), 
450:      $       QFLUC,FQCFOR,  !##FLUCQ 
451:      $       QETERM(ST2),NST2,ETERM(ST2),.FALSE. 
452:      -               ,.false.,0,0  !## VIBBLOCK 
453:      -                ) 
454: C 
455: ##IF GENBORN (gb1)454: ##IF GENBORN (gb1)
456: C  Generalized Born Solvation energy term455: C  Generalized Born Solvation energy term
457: C456: C
458: C  First call set-up and get alphas and sigmas for this config457: C  First call set-up and get alphas and sigmas for this config
459: C458: C
460:         IF (QGenBorn) THEN459:         IF (QGenBorn) THEN
461:            If ( NTrans .gt. 0 ) Call WrnDie(-3,'<ENERGY>',460:            If ( NTrans .gt. 0 ) Call WrnDie(-3,'<ENERGY>',
462:      *          'Images not implemented w/ Generalized Born.')461:      *          'Images not implemented w/ Generalized Born.')
463: ##IF QBLOCK (gb2)462: ##IF QBLOCK (gb2)
464: ##IF GBBLCK (gb3)463: ##IF GBBLCK (gb3)


r16977/toph19_eef1.1.inp 2017-01-21 10:32:53.255188235 +0000 r16976/toph19_eef1.1.inp 2017-01-21 10:32:55.371188235 +0000
586: RESI HSC     0.00000                ! Doubly protonated histidine.  Ring charges586: RESI HSC     0.00000                ! Doubly protonated histidine.  Ring charges
587: GROU                                587: GROU                                
588: ATOM N    NH1    -0.35                ! Hayes and Kollman jacs 98:3335 (1976).588: ATOM N    NH1    -0.35                ! Hayes and Kollman jacs 98:3335 (1976).
589: ATOM H    H       0.25                !  TL modif: neutralize589: ATOM H    H       0.25                !  TL modif: neutralize
590: ATOM CA   CH1E    0.10590: ATOM CA   CH1E    0.10
591: GROU591: GROU
592: ATOM CB   CH2E    0.00592: ATOM CB   CH2E    0.00
593: ATOM CG   CR      0.05593: ATOM CG   CR      0.05
594: ATOM CD2  CR1E    0.05594: ATOM CD2  CR1E    0.05
595: GROU595: GROU
596: !bs360: changed back to NH1596: ATOM ND1  NC2    -0.55
597: !ATOM ND1  NC2    -0.55 
598: ATOM ND1  NH1    -0.55 
599: ATOM HD1  H       0.45597: ATOM HD1  H       0.45
600: GROU598: GROU
601: ATOM CE1  CR1E    0.10599: ATOM CE1  CR1E    0.10
602: GROU600: GROU
603: !ATOM NE2  NC2    -0.55601: ATOM NE2  NC2    -0.55
604: ATOM NE2  NH1    -0.55 
605: ATOM HE2  H       0.45602: ATOM HE2  H       0.45
606: GROU603: GROU
607: ATOM C    C       0.55604: ATOM C    C       0.55
608: ATOM O    O      -0.55605: ATOM O    O      -0.55
609: BOND N    CA        CA   C         C    +N        C    O         N    H606: BOND N    CA        CA   C         C    +N        C    O         N    H
610: BOND CA   CB        CB   CG        CG   ND1       CG   CD2       ND1  HD1607: BOND CA   CB        CB   CG        CG   ND1       CG   CD2       ND1  HD1
611: BOND ND1  CE1       CD2  NE2       CE1  NE2          NE2        HE2608: BOND ND1  CE1       CD2  NE2       CE1  NE2          NE2        HE2
612: DIHE -C   N    CA   C         N    CA   C    +N        CA   C    +N   +CA609: DIHE -C   N    CA   C         N    CA   C    +N        CA   C    +N   +CA
613: DIHE N    CA   CB   CG        CA   CB   CG   ND1610: DIHE N    CA   CB   CG        CA   CB   CG   ND1
614: IMPH N    -C   CA   H         C    CA   +N   O         CA   N    C    CB611: IMPH N    -C   CA   H         C    CA   +N   O         CA   N    C    CB


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0