hdiff output

r21974/cheq 2017-01-21 10:33:57.586250000 +0000 r21973/cheq 2017-01-21 10:34:01.318250000 +0000
  1: svn: warning: W195007: URL 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/source/cheq' refers to a directory  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/source/cheq' in revision 21973
  2: svn: E200009: Could not cat all targets because some targets are directories 
  3: svn: E200009: Illegal target for the requested operation 


r21974/cheq.mk 2017-01-21 10:33:56.778250000 +0000 r21973/cheq.mk 2017-01-21 10:34:00.342250000 +0000
  1: # cheq makefile  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/build/UNX/cheq.mk' in revision 21973
  2: # cheq library rules 
  3: OBJS_cheq= \ 
  4:         $(LIB)/cheq.a(fqcom.o) \ 
  5:         $(LIB)/cheq.a(fqener.o) \ 
  6:         $(LIB)/cheq.a(fqminut.o) 
  7: # 
  8: $(LIB)/cheq.a : $(OBJS_cheq) 
  9:         $(RANLIB) $(LIB)/cheq.a 
 10:         @echo cheq COMPLETED 
 11: # 
 12: # cheq source file rules 
 13: $(LIB)/cheq.a(fqcom.o)  : $(SRC)/cheq/fqcom.src 
 14:         $(FLX) $(SRC)/cheq/fqcom.src 
 15:         $(FC2) fqcom.f 
 16:         $(AR_COMMAND) $(LIB)/cheq.a fqcom.o 
 17:         $(REMOVE_F) fqcom.f 
 18:         $(REMOVE_O) fqcom.o 
 19: # 
 20: $(LIB)/cheq.a(fqener.o)  : $(SRC)/cheq/fqener.src 
 21:         $(FLX) $(SRC)/cheq/fqener.src 
 22:         $(FC2) fqener.f 
 23:         $(AR_COMMAND) $(LIB)/cheq.a fqener.o 
 24:         $(REMOVE_F) fqener.f 
 25:         $(REMOVE_O) fqener.o 
 26: # 
 27: $(LIB)/cheq.a(fqminut.o)  : $(SRC)/cheq/fqminut.src 
 28:         $(FLX) $(SRC)/cheq/fqminut.src 
 29:         $(FC2) fqminut.f 
 30:         $(AR_COMMAND) $(LIB)/cheq.a fqminut.o 
 31:         $(REMOVE_F) fqminut.f 
 32:         $(REMOVE_O) fqminut.o 
 33: # 
 34: # 
 35: # cheq dependency file 
 36: # 
 37: $(LIB)/cheq.a(fqcom.o) : cheqdyn.fcm 
 38: $(LIB)/cheq.a(fqcom.o) : comand.fcm 
 39: $(LIB)/cheq.a(fqcom.o) : consta.fcm 
 40: $(LIB)/cheq.a(fqcom.o) : coord.fcm 
 41: $(LIB)/cheq.a(fqcom.o) : dimens.fcm 
 42: $(LIB)/cheq.a(fqcom.o) : exfunc.fcm 
 43: $(LIB)/cheq.a(fqcom.o) : heap.fcm 
 44: $(LIB)/cheq.a(fqcom.o) : impnon.fcm 
 45: $(LIB)/cheq.a(fqcom.o) : number.fcm 
 46: $(LIB)/cheq.a(fqcom.o) : param.fcm 
 47: $(LIB)/cheq.a(fqcom.o) : psf.fcm 
 48: $(LIB)/cheq.a(fqcom.o) : rtf.fcm 
 49: $(LIB)/cheq.a(fqcom.o) : stack.fcm 
 50: $(LIB)/cheq.a(fqcom.o) : stream.fcm 
 51: # 
 52: # 
 53: $(LIB)/cheq.a(fqener.o) : cheqdyn.fcm 
 54: $(LIB)/cheq.a(fqener.o) : consta.fcm 
 55: $(LIB)/cheq.a(fqener.o) : derivq.fcm 
 56: $(LIB)/cheq.a(fqener.o) : dimens.fcm 
 57: $(LIB)/cheq.a(fqener.o) : energy.fcm 
 58: $(LIB)/cheq.a(fqener.o) : exfunc.fcm 
 59: $(LIB)/cheq.a(fqener.o) : fast.fcm 
 60: $(LIB)/cheq.a(fqener.o) : flucq.fcm 
 61: $(LIB)/cheq.a(fqener.o) : heap.fcm 
 62: $(LIB)/cheq.a(fqener.o) : image.fcm 
 63: $(LIB)/cheq.a(fqener.o) : impnon.fcm 
 64: $(LIB)/cheq.a(fqener.o) : inbnd.fcm 
 65: $(LIB)/cheq.a(fqener.o) : number.fcm 
 66: $(LIB)/cheq.a(fqener.o) : parallel.fcm 
 67: $(LIB)/cheq.a(fqener.o) : param.fcm 
 68: $(LIB)/cheq.a(fqener.o) : psf.fcm 
 69: $(LIB)/cheq.a(fqener.o) : stack.fcm 
 70: $(LIB)/cheq.a(fqener.o) : stream.fcm 
 71: # 
 72: # 
 73: $(LIB)/cheq.a(fqminut.o) : cheqdyn.fcm 
 74: $(LIB)/cheq.a(fqminut.o) : dimens.fcm 
 75: $(LIB)/cheq.a(fqminut.o) : impnon.fcm 
 76: # 


r21974/compile.csh 2017-01-21 10:33:57.014250000 +0000 r21973/compile.csh 2017-01-21 10:34:00.670250000 +0000
  1: #!/bin/tcsh  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/compile.csh' in revision 21973
  2:  
  3: # tested on 'clust' with 'ifort/64/11.1/038' module 
  4: #./install.com gnu large ifort keepo keepf > build.log 
  5:  
  6: # tested on 'clust' with 'pgi/64/7.2/5' module 
  7: ./install.com gnu large PGF90 keepo keepf AMD OPT > build.log 


r21974/fqcom.src 2017-01-21 10:33:57.794250000 +0000 r21973/fqcom.src 2017-01-21 10:34:01.634250000 +0000
  1: CHARMM Element source/flucq/fqcom.src 0.0  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/source/cheq/fqcom.src' in revision 21973
  2: ##IF CHEQ (cheq_main) 
  3: C----------------------------------------------------------------------- 
  4:       SUBROUTINE CHEQPREP(COMLYN,COMLEN) 
  5:  
  6: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
  7: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
  8: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
  9: ##INCLUDE '~/charmm_fcm/number.fcm' 
 10: C 
 11:       CHARACTER*(*) COMLYN 
 12:       INTEGER COMLEN 
 13: C 
 14: ##INCLUDE '~/charmm_fcm/coord.fcm' 
 15: ##INCLUDE '~/charmm_fcm/stream.fcm' 
 16: ##INCLUDE '~/charmm_fcm/psf.fcm' 
 17: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
 18: ##INCLUDE '~/charmm_fcm/heap.fcm' 
 19: C 
 20:       LOGICAL   WALL 
 21:       LOGICAL   MASQ 
 22:  
 23:       WALL=INDXA(COMLYN,COMLEN,'WALL').GT.0 
 24:  
 25:       IF ( (WALL) .AND. (CHEQALLHPCNTRWALL.EQ.1) ) THEN 
 26: C     allocate long-term heap space 
 27:       WALLPA1=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 28:      $              'WALLPA1') 
 29:       WALLPB1=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 30:      $              'WALLPB1') 
 31:       WALLPA2=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 32:      $              'WALLPA2') 
 33:       WALLPB2=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 34:      $             'WALLPB2') 
 35:       WALLPQ1=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 36:      $            'WALLPQ1') 
 37:       WALLPQ2=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 38:      $            'WALLPQ2') 
 39:       WALLPK=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 40:      $            'WALLPK') 
 41:       CHEQALLHPCNTRWALL = CHEQALLHPCNTRWALL + 1    ! ONE-TIME ASSIGNMENT ONLY 
 42:       QHPWALLP=.TRUE. 
 43:       ELSE 
 44:       ENDIF 
 45:  
 46:  
 47:       IF (CHEQALLHPCNTRMASQ.EQ.1) THEN 
 48:       CMASSQ=ALLHP(IREAL8(NATOM),'fqcom.src','CHEQPREP', 
 49:      $               'CMASSQ') 
 50:  
 51:       CHEQALLHPCNTRMASQ  = CHEQALLHPCNTRMASQ + 1 
 52:       QHPMASQ=.TRUE. 
 53:       ELSE 
 54:       ENDIF 
 55:  
 56:       CALL FQCOM(COMLYN,COMLEN) 
 57:  
 58:       RETURN 
 59:       END        ! CHEQPREP 
 60:  
 61: C----------------------------------------------------------------------- 
 62:       SUBROUTINE CHEQSTOP 
 63:  
 64: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
 65: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
 66: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
 67: C##INCLUDE '~/charmm_fcm/number.fcm' 
 68: C 
 69: C      INTEGER NATOM 
 70: C 
 71: C##INCLUDE '~/charmm_fcm/coord.fcm' 
 72: C##INCLUDE '~/charmm_fcm/stream.fcm' 
 73: ##INCLUDE '~/charmm_fcm/psf.fcm' 
 74: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
 75: C##INCLUDE '~/charmm_fcm/heap.fcm' 
 76:  
 77:  
 78:       IF (QHPWALLP) THEN 
 79: c   release heap space from CHEQ wall potential varibables 
 80:       CALL FREHP(WALLPA1,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 81:      $          'WALLPA1') 
 82:       CALL FREHP(WALLPB1,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 83:      $          'WALLPB1') 
 84:       CALL FREHP(WALLPA2,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 85:      $          'WALLPA2') 
 86:       CALL FREHP(WALLPB2,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 87:      $          'WALLPB2') 
 88:       CALL FREHP(WALLPQ1,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 89:      $          'WALLPQ1') 
 90:       CALL FREHP(WALLPQ2,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 91:      $          'WALLPQ2') 
 92:       CALL FREHP(WALLPK,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
 93:      $          'WALLPK') 
 94:  
 95:        ELSE 
 96:        ENDIF 
 97:  
 98:       IF (QHPMASQ) THEN 
 99:       CALL FREHP(CMASSQ,IREAL8(NATOM),'fqcom.src','CHEQSTOP', 
100:      $          'CMASSQ') 
101:       ENDIF 
102:  
103:       RETURN 
104:       END     ! CHEQSTOP 
105:  
106: C------------------------------------------------------------------------------ 
107:       SUBROUTINE FQCOM(COMLYN,COMLEN) 
108: C This subroutine is called from CHARMM (main) to process FLUQ commands 
109: C 
110: C written by Tom Cleveland 
111: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
112: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
113: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
114: ##INCLUDE '~/charmm_fcm/number.fcm' 
115: C 
116:       CHARACTER*(*) COMLYN 
117:       INTEGER COMLEN 
118: C 
119: ##INCLUDE '~/charmm_fcm/coord.fcm' 
120: ##INCLUDE '~/charmm_fcm/stream.fcm' 
121: ##INCLUDE '~/charmm_fcm/psf.fcm' 
122: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
123: ##INCLUDE '~/charmm_fcm/heap.fcm' 
124: C 
125:       INTEGER FLAGS(MAXAIM) 
126:       CHARACTER*4 QNRMTYP 
127:       INTEGER   I,WALN 
128:       REAL*8  CGMA,TSTA 
129:       REAL*8  QRA1,QRB1,QRA2,QRB2,QRQ1,QRQ2,QRK 
130:       LOGICAL   LUSED,QPRINT,QFLEX,QRESET 
131:       LOGICAL   QMASS,WALPP, WTYP 
132: C 
133: C 
134:       CALL SELCTA(COMLYN,COMLEN,FLAGS,X,Y,Z,WMAIN,.TRUE.) 
135:       QNRMTYP=GTRMA(COMLYN,COMLEN,'NORM') 
136:       IF (QNRMTYP.NE.'') THEN 
137:        CALL QPRTN(QNRMTYP,FLAGS,NATOM) 
138:       ENDIF 
139:       QCG=(INDXA(COMLYN,COMLEN,'ON').GT.0) 
140:       QCG=(.NOT.(INDXA(COMLYN,COMLEN,'OFF').GT.0)) 
141:       IF (.NOT.QCG) THEN 
142:        CGMODEL=0 
143:        QCG=.FALSE. 
144:        QCGINV=.FALSE. 
145:        QCGINVF=.FALSE. 
146:        QCG=.FALSE. 
147:        QNOCO=.FALSE. 
148:        QPOLAR1=.FALSE. 
149:        IPOLAR1=-1 
150:        QCGWATER=.FALSE. 
151:       ENDIF 
152:       QRESET=INDXA(COMLYN,COMLEN,'RESE').GT.0 
153:       IF (QRESET)  CALL FQRESET(NATOM) 
154:       QNOCO=(INDXA(COMLYN,COMLEN,'NOCO').GT. 0) 
155:       IF(QPRNETA) WRITE(*,*) "CGINV will write out Eta matrix" 
156:       CGMODEL=GTRMI(COMLYN,COMLEN,'CGMD',CGMODEL) 
157:       CGEQ=GTRMI(COMLYN,COMLEN,'CGEQ',0) 
158:       QCGINV=INDXA(COMLYN,COMLEN,'CGIN').GT.0 
159:       QCGPOL=INDXA(COMLYN,COMLEN,'POLT').GT.0 
160: c      QCGINVF=INDXA(COMLYN,COMLEN,'CGFC').GT.0 
161:       QCGWATER=INDXA(COMLYN,COMLEN,'WATE').GT.0 
162:       QCGSPC=INDXA(COMLYN,COMLEN,'SPC').GT.0 
163:       QFLEX=INDXA(COMLYN,COMLEN,'FLEX').GT.0 
164:       QT4P=INDXA(COMLYN,COMLEN,'TIP4').GT.0 
165:       IF (QCGWATER)  CALL TYPMOL('WATE',FLAGS) 
166:       IF (QCGSPC) CALL TYPMOL('SPC ',FLAGS) 
167:       IF (QT4P) THEN  
168:       GQT4P=.TRUE. 
169:       CALL TYPMOL('TIP4',FLAGS) 
170:       write(*,*) " GLOBAL QT4P = ", GQT4P 
171:        ENDIF  
172:       IF (QFLEX) CALL TYPMOL('FLEX',FLAGS) 
173:       QPOLAR1=INDXA(COMLYN, COMLEN, 'QPOL') .GT. 0 
174:       IPOLAR1=GTRMI(COMLYN,COMLEN,'IPOL',IPOLAR1) 
175:       QPRINT=(INDXA(COMLYN,COMLEN,'PRIN').GT.0) 
176:       QREF=(INDXA(COMLYN,COMLEN,'QREF').GT.0) 
177:  
178:       QMASS=INDXA(COMLYN,COMLEN,'QMAS').GT.0 
179:       IF (QMASS) THEN 
180:       CGMA=GTRMF(COMLYN,COMLEN,'CGMA',DEFQ) 
181:       TSTA=GTRMF(COMLYN,COMLEN,'TSTA',TCGDEF) 
182:       CALL ASSIGNQMAS(FLAGS,CGMA,TSTA,HEAP(CMASSQ)) 
183:       ENDIF 
184:       WTYP=INDXA(COMLYN,COMLEN,'WTYP').GT.0  ! Determine potential type (harmonic or other) 
185: c                                             1=Harmonic (preferred); 2=Wall Potential 
186:       IF (WTYP) THEN 
187:       PTYP=GTRMI(COMLYN,COMLEN,'PTYP',1) 
188:        IF (PTYP.EQ.2) THEN 
189:             QRA1=GTRMF(COMLYN,COMLEN,'QRA1',QRA1DEF) 
190:             QRB1=GTRMF(COMLYN,COMLEN,'QRB1',QRB1DEF) 
191:             QRA2=GTRMF(COMLYN,COMLEN,'QRA2',QRA2DEF) 
192:             QRB2=GTRMF(COMLYN,COMLEN,'QRB2',QRB2DEF) 
193:             QRQ1=GTRMF(COMLYN,COMLEN,'QRQ1',QRQ1DEF) 
194:             QRQ2=GTRMF(COMLYN,COMLEN,'QRQ2',QRQ2DEF) 
195:             QRK=GTRMF(COMLYN,COMLEN,'QRK',QRKDEF) 
196:             wallpN=GTRMI(COMLYN,COMLEN,'WALN',1) 
197:          ELSEIF (PTYP.EQ.1) THEN 
198:             QRQ1=GTRMF(COMLYN,COMLEN,'QRQ1',QRQ1DEF) 
199:             QRQ2=GTRMF(COMLYN,COMLEN,'QRQ2',QRQ2DEF) 
200:             QRK=GTRMF(COMLYN,COMLEN,'QRK',QRKDEF) 
201:          ELSE 
202:          ENDIF 
203:       CALL ASSWALLPP(FLAGS,QRA1,QRB1,QRA2,QRB2,QRQ1,QRQ2,QRK, 
204:      $               HEAP(wallpa1),HEAP(wallpb1),HEAP(wallpa2), 
205:      $               HEAP(wallpb2),HEAP(wallpq1),HEAP(wallpq2), 
206:      $               HEAP(wallpK),NATOM) 
207:       ENDIF 
208:  
209:  
210:       CALL XTRANE(COMLYN,COMLEN,'FQCOM') 
211:  
212:       IF (QREF) THEN 
213:        CALL FILLCGREF(FLAGS) 
214:       ENDIF 
215:  
216:       IF(QPRINT) THEN 
217:          IF (QCGINV) WRITE(OUTU,'(a)') "<FQCOM>: CG INVERSION REQUESTED" 
218:          IF (QCGINVF) WRITE(OUTU,'(a)') "<FQCOM>: CG FORCE REQUESTED" 
219:          WRITE(OUTU,'(a,i6)') "<FQCOM>: CGMODEL SET TO",CGMODEL 
220:          IF (QPOLAR1) WRITE(OUTU,'(2A,I5)') 
221:      $      "Single Molecule Polarization is requested for " 
222:      $      ,"Molecule #",IPOLAR1 
223:          IF (QNOCO) WRITE(OUTU,'(A)') 
224:      $      "Charge only, coordinates being fixed." 
225:          IF (QCGWATER) WRITE(OUTU,'(A)') 
226:      $      "Hardness ontribution to Force set to zero" 
227:       DO I=1,NATOM 
228:       WRITE (OUTU,20) I,QPARTBIN(I),QPART(I),QPRTPTR(I) 
229:       ENDDO 
230:        DO I = 1,NATOM 
231:        WRITE(*,*)"ATOM,MOLTYPE=",I,MOLTYPE(I) 
232:        ENDDO 
233:       ENDIF !QPRINT 
234:  
235:  
236: C 
237:  20   FORMAT(4i9) 
238:       RETURN 
239:       END !FQCOM 
240: C----------------------------------------------------------------------- 
241:       SUBROUTINE FQRESET(NATOM) 
242: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
243: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
244: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
245:  
246:       INTEGER I,NATOM 
247:        
248:       QNPART=0 
249:       QCG=.FALSE. 
250:       QPRNETA=.FALSE. 
251:       DO I=1,MAXAIM 
252:        QPARTBIN(I)=0 
253:        QPART(I)=0 
254:        MOLTYPE(I)=0 
255:        QPARTED(I)=.FALSE. 
256:       ENDDO 
257:       DO I=1,MAXGRP+1 
258:        QPRTPTR(I)=0 
259:       ENDDO 
260:  
261:       RETURN 
262:       END 
263: C----------------------------------------------------------------------- 
264:       SUBROUTINE QPRTN(QNRMTYP,FLAGS,NATOM) 
265: C  This subroutine is called from FQCOM to define the groups of atoms  
266: C  over which charge will be normalized. 
267: C 
268: C  QNRMTYP character       type of atom grouping 
269: C  FLAGS   integer NATOM    
270: C 
271: C  QNPART  integer scalar  total number of partitions      
272: C  QPRTPTR integer MAXPART pointers into QPART for end of current partition 
273: C  QPART   integer NATOM   atom numbers in partition order 
274: C  QPARTBIN integer NATOM   partition to which atom number belongs 
275: C       
276: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
277: ##INCLUDE '~/charmm_fcm/number.fcm' 
278: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
279: ##INCLUDE '~/charmm_fcm/stream.fcm' 
280: C 
281:       INTEGER FLAGS(*),NATOM 
282:       CHARACTER*4 QNRMTYP 
283:       INTEGER I,J,PARTCNTR,COUNTATOM 
284:       LOGICAL PARTFOUND 
285: C 
286: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
287: C 
288: C  assign new partition numbers to QPARTBIN depending on FLAGS array and 
289: C  type of partitioning QNRMTYP 
290: C 
291: c      write (*,*) '<QPRTN> before QPRTBY??? calls' 
292:       IF     (QNRMTYP.EQ.'BYRE') THEN 
293: C       CALL QPRTBYRES(FLAGS,QNPART) 
294:        CALL QPRTBYRES(FLAGS,QNPART,QPARTED) 
295:       ELSEIF (QNRMTYP.EQ.'BYAL') THEN 
296:        CALL QPRTBYALL(FLAGS,QNPART) 
297:       ELSEIF (QNRMTYP.EQ.'BYSE') THEN 
298:        CALL QPRTBYSEG(FLAGS,QNPART,QPARTED) 
299:       ELSEIF (QNRMTYP.EQ.'BYGR') THEN 
300:        CALL QPRTBYGRP(FLAGS,QNPART,QPARTED) 
301:       ELSEIF (QNRMTYP.EQ.'BYMO') THEN 
302:        CALL QPRTBYMOL(FLAGS) 
303:       ELSEIF (QNRMTYP.EQ.'NOFQ') THEN 
304:        CALL QPRTNOFQ(FLAGS,QPARTBIN,NATOM) 
305:       ELSEIF (QNRMTYP.EQ.'BYR1') THEN 
306:        CALL QPRTBYR1(FLAGS,QNPART) 
307:       ENDIF 
308:  
309: C  now merge the new elements from the FLAGS array with the old ones  
310: C  from the QPARTBIN array 
311: C 
312:       IF(QNRMTYP.NE.'NOFQ') THEN 
313:        CALL PARTMERGE(QPARTBIN,FLAGS,NATOM) 
314:       ENDIF 
315: C       
316: C  counts the number of partitions, compacts the partition bin array,  
317: C  and fills the partition occupancy and partition pointer array 
318: C 
319:       PARTCNTR=1 
320:       QPRTPTR(1)=0 
321:       COUNTATOM=0 
322:       DO I=1,QNPART 
323:        PARTFOUND=.FALSE. 
324:        DO J=1,NATOM 
325:          IF(QPARTBIN(J).EQ.I) THEN 
326:             COUNTATOM=COUNTATOM+1 
327:             QPARTBIN(J)=PARTCNTR 
328:             PARTFOUND=.TRUE. 
329:             QPART(COUNTATOM)=J 
330:           ENDIF 
331:        ENDDO 
332:        IF(PARTFOUND) THEN 
333:           QPRTPTR(PARTCNTR+1)=COUNTATOM 
334:           PARTCNTR=PARTCNTR+1 
335:        ENDIF 
336:       ENDDO 
337:       QNPART=PARTCNTR-1 
338:        
339:       IF (QNPART.GT.1) QMOLC=.TRUE. 
340:        
341:  
342:       RETURN 
343:       END !QNORMAL 
344: C----------------------------------------------------------------------- 
345:       SUBROUTINE QPRTBYRES(FLAGS,QNPART,QPARTED) 
346: C This subroutine is called from QPRTN and sets new partition  
347: C numbers for the FLAGS array by residue.  also sets new QNPART 
348: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
349: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
350: ##INCLUDE '~/charmm_fcm/psf.fcm' 
351: ##INCLUDE '~/charmm_fcm/stream.fcm'  
352:  
353:       INTEGER FLAGS(NATOM),QNPART 
354:       INTEGER I,J,PARTCTR,PARTATOMS,RESATOMS 
355:       LOGICAL FOUND,QPARTED(*) 
356: C 
357:       PARTCTR=1 
358:       DO I=1,NRES 
359:       FOUND=.FALSE. 
360:       PARTATOMS=0 
361:       DO J=(IBASE(I)+1),(IBASE(I+1)) 
362:            
363:           IF(FLAGS(J).NE.0) THEN 
364:            IF (.not.QPARTED(J)) THEN 
365:             QPARTED(J)=.TRUE. 
366:             FOUND=.TRUE. 
367:             FLAGS(J)=QNPART+PARTCTR 
368:             PARTATOMS=PARTATOMS+1 
369:            ELSE 
370:                if(prnlev.gt.1)write(outu,49) 
371:      $        '***CHEQ*** Attempting to assign an atom to more than',  
372:      2        'one normalization unit (RESIDUE); NOT A VALID SCHEME ' 
373:           CALL WRNDIE(-1,'<QPRTBYRES>', 
374:      $        'CHEQ: CHECK NORMALIZATION SPECIFICATIONS ') 
375:          ENDIF 
376:          ENDIF 
377: 49     format(/,3x,A,/,3x,A) 
378:        ENDDO 
379:         IF (FOUND) THEN 
380:           RESATOMS=IBASE(I+1)-IBASE(I) 
381:           IF(RESATOMS.GT.PARTATOMS) THEN  
382:             CALL WRNDIE(0,'<FQCOM>','All atoms of residue not selected') 
383:           ENDIF 
384:           PARTCTR=PARTCTR+1 
385:        ENDIF 
386:       ENDDO 
387:       QNPART=QNPART+PARTCTR-1 
388: C 
389:       RETURN 
390:       END 
391: C----------------------------------------------------------------------- 
392:       SUBROUTINE QPRTBYSEG(FLAGS,QNPART,QPARTED) 
393: C This subroutine is called from QPRTN and sets new partition  
394: C numbers for the FLAGS array by segment.  also sets new QNPART 
395: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
396: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
397: ##INCLUDE '~/charmm_fcm/psf.fcm' 
398: ##INCLUDE '~/charmm_fcm/stream.fcm'  
399:  
400:       INTEGER FLAGS(NATOM),QNPART 
401:       INTEGER I,J,PARTCTR,PARTATOMS,SEGATOMS 
402:       LOGICAL FOUND,QPARTED(*) 
403: C 
404:       PARTCTR=1 
405:       DO I=1,NSEG 
406:        FOUND=.FALSE. 
407:        PARTATOMS=0 
408:        DO J=(IBASE(NICTOT(I)+1)+1),(IBASE(NICTOT(I+1)+1)) 
409:          IF(FLAGS(J).NE.0) THEN 
410:           IF (.not.QPARTED(J)) THEN 
411:             QPARTED(J) = .TRUE. 
412:             FOUND=.TRUE. 
413:             FLAGS(J)=QNPART+PARTCTR 
414:             PARTATOMS=PARTATOMS+1 
415:           ELSE 
416:             if (PRNLEV.GT.1) write(outu,39) 
417:      $   '***CHEQ*** Attempting to assign an atom to more than', 
418:      2   'one normalization unit (SEGMENT): NOT A VALID SCHEME' 
419:         CALL WRNDIE(-1,'<QPRTBYSEG>', 
420:      $ 'CHEQ: CHECK NORMALIZATION SPECIFICATION!') 
421:           ENDIF 
422:           ENDIF 
423: 39     format(/,3x,A,/,3x,A) 
424:        ENDDO 
425:        IF (FOUND) THEN 
426:           SEGATOMS=IBASE(NICTOT(I+1)+1)-(IBASE(NICTOT(I)+1)) 
427:           IF (SEGATOMS.GT.PARTATOMS) THEN 
428:             CALL WRNDIE(0,'<FQCOM>','All atoms of segment not selected') 
429:           ENDIF 
430:           PARTCTR=PARTCTR+1 
431:        ENDIF 
432:       ENDDO 
433:       QNPART=QNPART+PARTCTR-1 
434: C 
435:       RETURN 
436:       END 
437: C----------------------------------------------------------------------- 
438:       SUBROUTINE QPRTBYGRP(FLAGS,QNPART,QPARTED) 
439: C This subroutine is called from QPRTN and sets new partition  
440: C numbers for the FLAGS array by segment.  also sets new QNPART 
441: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
442: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
443: ##INCLUDE '~/charmm_fcm/psf.fcm' 
444: ##INCLUDE '~/charmm_fcm/stream.fcm' 
445:       INTEGER FLAGS(NATOM),QNPART 
446:       INTEGER I,J,PARTCTR,PARTATOMS,GRPATOMS 
447:       LOGICAL FOUND,QPARTED(*) 
448: C 
449:       PARTCTR=1 
450:       DO I=1,NGRP 
451:        FOUND=.FALSE. 
452:        PARTATOMS=0 
453:        DO J=(IGPBS(I)+1),(IGPBS(I+1)) 
454:           IF(FLAGS(J).NE.0) THEN 
455:            IF (.not.QPARTED(J)) THEN 
456:             QPARTED(J) = .TRUE. 
457:             FOUND=.TRUE. 
458:             FLAGS(J)=QNPART+PARTCTR 
459:             PARTATOMS=PARTATOMS+1 
460:            ELSE 
461:              if(prnlev.gt.1)write(outu,59) 
462:      $        '***CHEQ*** Attempting to assign an atom to more than',  
463:      2        'one normalization unit (GROUP); NOT A VALID SCHEME ' 
464:               CALL WRNDIE(-1,'<QPRTBYGROUP>', 
465:      $        'CHEQ: CHECK NORMALIZATION SPECIFICATIONS ') 
466: 59     format(/,3x,A,/,3x,A) 
467:            ENDIF 
468:           ENDIF 
469:        ENDDO 
470:        IF (FOUND) THEN 
471:           GRPATOMS=IGPBS(I+1)-IGPBS(I) 
472:           IF (GRPATOMS.GT.PARTATOMS) THEN 
473:             CALL WRNDIE(0,'<FQCOM>','All atoms of group not selected') 
474:           ENDIF 
475:           PARTCTR=PARTCTR+1 
476:        ENDIF 
477:       ENDDO 
478:       QNPART=QNPART+PARTCTR-1 
479: C 
480:       RETURN 
481:       END 
482: C----------------------------------------------------------------------- 
483:       SUBROUTINE QPRTBYALL(FLAGS,QNPART) 
484: C  This subroutine is called from QPRTN and sets new partition  
485: C  number for trues from FLAGS array.  also sets new QNPART 
486: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
487: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
488: ##INCLUDE '~/charmm_fcm/psf.fcm' 
489: C 
490:       INTEGER FLAGS(*),QNPART 
491:       INTEGER I 
492: C 
493:       QNPART=QNPART+1 
494:       DO I=1,NATOM 
495:        IF(FLAGS(I).NE.0) FLAGS(I)=QNPART 
496:       ENDDO 
497:       RETURN 
498:       END 
499: C 
500: C----------------------------------------------------------------------- 
501:       SUBROUTINE QPRTBYMOL(FLAGS) 
502: C  This subroutine is called from QPRTN and sets new partition 
503: C  numbers for nonzero values from FLAGS array.  also set new QNPART 
504: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
505: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
506: ##INCLUDE '~/charmm_fcm/psf.fcm' 
507: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
508: ##INCLUDE '~/charmm_fcm/stream.fcm' 
509: C       
510:       INTEGER FLAGS(*) 
511:       INTEGER I,J,PARTCTR,PARTATOMS,MOLATOMS 
512:       LOGICAL FOUND 
513: C 
514:       PARTCTR=1 
515:       DO I=1,MOLNT 
516:        FOUND=.FALSE. 
517:        PARTATOMS=0 
518:        DO J=MOLPTR(I),(MOLPTR(I+1)-1) 
519:           IF(FLAGS(MOLBL(J)).NE.0) THEN 
520:            IF (.not.QPARTED(MOLBL(J))) THEN 
521:             QPARTED(MOLBL(J))=.TRUE. 
522:             FOUND=.TRUE. 
523:             FLAGS(MOLBL(J))=QNPART+PARTCTR 
524:             PARTATOMS=PARTATOMS+1 
525:            ELSE 
526:                      if(prnlev.gt.1)write(outu,29) 
527:      $      '***CHEQ*** Attempting to assign an atom to more than', 
528:      2      'one normalization unit (MOLECULE); NOT A VALID SCHEME ' 
529:           CALL WRNDIE(-1,'<QPRTBYMOL>', 
530:      $        'CHEQ: CHECK NORMALIZATION SPECIFICATIONS ') 
531:           ENDIF 
532:          ENDIF 
533: 29    format(/,3x,A,/,3x,A) 
534:        ENDDO 
535:        IF (FOUND) THEN 
536:           MOLATOMS=MOLPTR(I+1)-MOLPTR(I) 
537:           IF (MOLATOMS.GT.PARTATOMS) THEN 
538:            CALL WRNDIE(0,'<FQCOM>','All atoms of molecule not selected') 
539:           ENDIF 
540:           PARTCTR=PARTCTR+1 
541:        ENDIF 
542:       ENDDO 
543:       QNPART=QNPART+PARTCTR-1 
544: C 
545:       RETURN 
546:       END 
547: C----------------------------------------------------------------------- 
548:       SUBROUTINE QPRTNOFQ(FLAGS,QPARTBIN,NATOM) 
549: C  This subroutine is called from QPRTN and sets   
550: C   
551: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
552: C 
553:       INTEGER FLAGS(*),QPARTBIN(*),NATOM 
554:       INTEGER I 
555: C 
556:       DO I=1,NATOM 
557:        IF(FLAGS(I).NE.0) QPARTBIN(I)=0 
558:       ENDDO 
559:       RETURN 
560:       END 
561: C 
562: C----------------------------------------------------------------------- 
563:       SUBROUTINE QPRTBYR1(FLAGS,QNPART) 
564: C This subroutine is called from QPRTN and sets new partition 
565: C numbers for the FLAGS array by backbone and sidechain.  also sets new QNPART 
566: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
567: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
568: ##INCLUDE '~/charmm_fcm/psf.fcm' 
569: ##INCLUDE '~/charmm_fcm/rtf.fcm' 
570:  
571:       INTEGER FLAGS(NATOM),QNPART 
572:       INTEGER I,J,PARTCTR,PARTATOMS,RESATOMS 
573:  
574:       INTEGER ITHISSEG,NUMRES,FIRRES, LASRES,FIRTYP,LASTYP 
575:       LOGICAL FOUND 
576: C 
577:  
578:  
579:  
580:  
581: c   loop over NSEG and find which segment to define partitions over  
582: c    in this manner 
583:  
584:  
585:         DO I = 1,NSEG 
586:         J=(IBASE(NICTOT(I)+1)+1) 
587:          IF(FLAGS(J).NE.0) ITHISSEG = I 
588:        ENDDO 
589:  
590: C    Find number of residues in this segment 
591:  
592:         NUMRES = NICTOT(ITHISSEG+1) - ( NICTOT(ITHISSEG) )  
593:         FIRRES = NICTOT(ITHISSEG)+1 
594:         LASRES = NICTOT(ITHISSEG+1) 
595:  
596:  
597:        write(*,*) " SEGMENT TO APPLY NORM. METHOD = ", ITHISSEG 
598:        write(*,*) " Number of residues = ", NUMRES 
599:        write(*,*) " First , Last Residue Number =", FIRRES, LASRES 
600:        write(*,*) " NRES, QNPART  = ", NRES, QNPART 
601:  
602:       IF (NUMRES.GE.3) THEN 
603:  
604:       PARTCTR=1 
605:        FIRTYP=1 
606:  
607:           IF (FIRTYP.EQ.1) THEN   ! ACE 
608:  
609:           FOUND=.FALSE. 
610:           PARTATOMS=0 
611:           I = FIRRES 
612:            DO J = (IBASE(I)+1), (IBASE(I)+6) 
613:            IF(FLAGS(J).NE.0) THEN 
614:            FOUND=.TRUE. 
615:            FLAGS(J)=QNPART+PARTCTR 
616:            PARTATOMS=PARTATOMS+1 
617:            ENDIF 
618:            ENDDO 
619:            IF (FOUND) THEN 
620:            PARTCTR=PARTCTR+1 
621:            ENDIF 
622:  
623:            FOUND=.FALSE. 
624:            PARTATOMS=0 
625:            DO J = (IBASE(I)+7), (IBASE(I)+12) 
626:            IF(FLAGS(J).NE.0) THEN 
627:            FOUND=.TRUE. 
628:            FLAGS(J)=QNPART+PARTCTR 
629:            PARTATOMS=PARTATOMS+1 
630:            ENDIF 
631:            ENDDO 
632:  
633:  
634:            IF (FOUND) THEN 
635:            PARTCTR = PARTCTR + 1 
636:            ENDIF 
637:  
638:  
639:            FOUND=.FALSE. 
640:            PARTATOMS=0 
641:  
642:            DO J = (IBASE(I)+13), (IBASE(I+1)) !  (IBASE(I+1)-2) 
643:            IF(FLAGS(J).NE.0) THEN 
644:            FOUND=.TRUE. 
645:            FLAGS(J)=QNPART+PARTCTR 
646:            PARTATOMS=PARTATOMS+1 
647:            ENDIF 
648:            ENDDO 
649:            IF (FOUND) THEN 
650:            PARTCTR = PARTCTR + 1 
651:            ENDIF 
652:        
653:         ELSEIF (FIRTYP.EQ.2) THEN ! NTER 
654:  
655:      
656:          FOUND=.FALSE. 
657:           PARTATOMS=0 
658:           I = FIRRES 
659:            DO J = (IBASE(I)+1), (IBASE(I)+8) 
660:            IF(FLAGS(J).NE.0) THEN 
661:            FOUND=.TRUE. 
662:            FLAGS(J)=QNPART+PARTCTR 
663:            PARTATOMS=PARTATOMS+1 
664:            ENDIF 
665:            ENDDO 
666:            IF (FOUND) THEN 
667:            PARTCTR=PARTCTR+1 
668:            ENDIF 
669:  
670:            FOUND=.FALSE. 
671:            PARTATOMS=0 
672:            DO J = (IBASE(I)+9), (IBASE(I+1)) 
673:            IF(FLAGS(J).NE.0) THEN 
674:            FOUND=.TRUE. 
675:            FLAGS(J)=QNPART+PARTCTR 
676:            PARTATOMS=PARTATOMS+1 
677:            ENDIF 
678:            ENDDO 
679:            IF (FOUND) THEN 
680:            PARTCTR = PARTCTR + 1 
681:            ENDIF 
682:  
683:  
684:         ELSE  
685:         ENDIF  !  END CONDITIONAL FOR FIRTYP FOR FIRST RESIDUE 
686:  
687:       DO I=FIRRES+1,LASRES-1 
688:       FOUND=.FALSE. 
689:       PARTATOMS=0 
690:       DO J=(IBASE(I)+1),(IBASE(I)+6) 
691:           IF(FLAGS(J).NE.0) THEN 
692:             FOUND=.TRUE. 
693:             FLAGS(J)=QNPART+PARTCTR 
694:             PARTATOMS=PARTATOMS+1 
695:           ENDIF 
696:        ENDDO 
697:     
698:        IF (FOUND) THEN 
699:         PARTCTR = PARTCTR + 1 
700:        ENDIF 
701:        FOUND=.FALSE. 
702:        PARTATOMS=0 
703:            DO J = (IBASE(I)+7), (IBASE(I+1)) 
704:            IF(FLAGS(J).NE.0) THEN 
705:            FOUND=.TRUE. 
706:            FLAGS(J)=QNPART+PARTCTR 
707:            PARTATOMS=PARTATOMS+1 
708:            ENDIF 
709:            ENDDO 
710:            IF (FOUND) THEN 
711:            PARTCTR = PARTCTR + 1 
712:            ENDIF 
713:       ENDDO 
714:  
715:  
716:  
717: c    for first residue of this segment 
718: c     this is assumed to be patched with some patch residue 
719: c    the choices for the patch are hacked via hard-coding 
720: c    first patches accounted for are:   FIRTYP=1(ACE)   FIRTYP=2(NTER) 
721: c    last  patches accounted for are:   FIRTYP=1(CT3)   FIRTYP=2(CTER) 
722: c     these are taken to be common combinations of patch residues 
723: c    can be generalized later 
724:  
725:  
726:  
727:            IF (FIRTYP.EQ.1) THEN 
728:  
729:            I = LASRES      ! CT3 
730:  
731:            
732:            FOUND=.FALSE. 
733:            PARTATOMS=0 
734:            DO J= (IBASE(I)+1), (IBASE(I)+6) 
735:            IF(FLAGS(J).NE.0) THEN 
736:            FOUND=.TRUE. 
737:            FLAGS(J)=QNPART+PARTCTR 
738:            PARTATOMS=PARTATOMS+1 
739:            ENDIF 
740:            ENDDO 
741:  
742:            IF (FOUND) THEN 
743:            PARTCTR = PARTCTR + 1 
744:            ENDIF 
745:  
746:            FOUND=.FALSE. 
747:            PARTATOMS=0 
748:            DO J = (IBASE(I)+7), (IBASE(I+1)-6) 
749:            IF(FLAGS(J).NE.0) THEN 
750:            FOUND=.TRUE. 
751:            FLAGS(J)=QNPART+PARTCTR 
752:            PARTATOMS=PARTATOMS+1 
753:            ENDIF 
754:            ENDDO 
755:            IF (FOUND) THEN 
756:            PARTCTR = PARTCTR + 1 
757:            ENDIF 
758:  
759:            FOUND=.FALSE. 
760:            PARTATOMS=0 
761:            DO J = (IBASE(I+1)-5), (IBASE(I+1)) 
762:            IF(FLAGS(J).NE.0) THEN 
763:            FOUND=.TRUE. 
764:            FLAGS(J)=QNPART+PARTCTR 
765:            PARTATOMS=PARTATOMS+1 
766:            ENDIF 
767:            ENDDO 
768:            IF (FOUND) THEN 
769:            PARTCTR = PARTCTR + 1 
770:            ENDIF 
771:  
772:            ELSE  IF (FIRTYP.EQ.2) THEN 
773:            ELSE 
774:            ENDIF 
775:  
776:            
777:           QNPART=QNPART+PARTCTR-1 
778:  
779:  
780:           ELSEIF (NUMRES.EQ.1) THEN 
781:  
782: c       
783:           PARTCTR=1 
784:           FIRTYP = 1  
785:           IF (FIRTYP.EQ.1) THEN 
786:           FOUND=.FALSE. 
787:           PARTATOMS=0 
788:           I = FIRRES 
789:           DO J = (IBASE(I)+1), (IBASE(I)+6) 
790:           IF(FLAGS(J).NE.0) THEN 
791:            FOUND=.TRUE. 
792:            FLAGS(J)=QNPART+PARTCTR 
793:            PARTATOMS=PARTATOMS+1 
794:            ENDIF 
795:            ENDDO 
796:            IF (FOUND) THEN 
797:            PARTCTR = PARTCTR + 1 
798:            ENDIF 
799:  
800:  
801:           FOUND=.FALSE. 
802:           PARTATOMS=0 
803:           DO J = (IBASE(I)+7), (IBASE(I)+12) 
804:           IF(FLAGS(J).NE.0) THEN 
805:            FOUND=.TRUE. 
806:            FLAGS(J)=QNPART+PARTCTR 
807:            PARTATOMS=PARTATOMS+1 
808:            ENDIF 
809:            ENDDO 
810:            IF (FOUND) THEN 
811:            PARTCTR = PARTCTR + 1 
812:            ENDIF 
813:  
814:  
815:  
816:           FOUND=.FALSE. 
817:           PARTATOMS=0 
818:           DO J = (IBASE(I)+13), (IBASE(I+1)-6) 
819:           IF(FLAGS(J).NE.0) THEN 
820:            FOUND=.TRUE. 
821:            FLAGS(J)=QNPART+PARTCTR 
822:            PARTATOMS=PARTATOMS+1 
823:            ENDIF 
824:            ENDDO 
825:            IF (FOUND) THEN 
826:            PARTCTR = PARTCTR + 1 
827:            ENDIF 
828:  
829:  
830:          FOUND=.FALSE. 
831:           PARTATOMS=0 
832:           DO J = (IBASE(I+1)-5), (IBASE(I+1)) 
833:           IF(FLAGS(J).NE.0) THEN 
834:            FOUND=.TRUE. 
835:            FLAGS(J)=QNPART+PARTCTR 
836:            PARTATOMS=PARTATOMS+1 
837:            ENDIF 
838:            ENDDO 
839:            IF (FOUND) THEN 
840:            PARTCTR = PARTCTR + 1 
841:            ENDIF 
842:  
843:           ELSE 
844:           ENDIF  ! CONDITOINAL FOR FIRTYP 
845:  
846:           QNPART = QNPART + PARTCTR - 1 
847:  
848:  
849:  
850:           ELSEIF (NUMRES.EQ.2) THEN 
851:  
852:           PARTCTR=1 
853:           FIRTYP=1 
854:  
855:           IF (FIRTYP.EQ.1) THEN 
856:  
857:           I = FIRRES 
858:          FOUND=.FALSE. 
859:           PARTATOMS=0 
860:           I = FIRRES 
861:           DO J = (IBASE(I)+1), (IBASE(I)+6) 
862:           IF(FLAGS(J).NE.0) THEN 
863:            FOUND=.TRUE. 
864:            FLAGS(J)=QNPART+PARTCTR 
865:            PARTATOMS=PARTATOMS+1 
866:            ENDIF 
867:            ENDDO 
868:            IF (FOUND) THEN 
869:            PARTCTR = PARTCTR + 1 
870:            ENDIF 
871:  
872:  
873:           FOUND=.FALSE. 
874:           PARTATOMS=0 
875:           DO J = (IBASE(I)+7), (IBASE(I)+12) 
876:           IF(FLAGS(J).NE.0) THEN 
877:            FOUND=.TRUE. 
878:            FLAGS(J)=QNPART+PARTCTR 
879:            PARTATOMS=PARTATOMS+1 
880:            ENDIF 
881:            ENDDO 
882:            IF (FOUND) THEN 
883:            PARTCTR = PARTCTR + 1 
884:            ENDIF 
885:  
886:  
887:  
888:           FOUND=.FALSE. 
889:           PARTATOMS=0 
890:           DO J = (IBASE(I)+13), (IBASE(I+1)) 
891:           IF(FLAGS(J).NE.0) THEN 
892:            FOUND=.TRUE. 
893:            FLAGS(J)=QNPART+PARTCTR 
894:            PARTATOMS=PARTATOMS+1 
895:            ENDIF 
896:            ENDDO 
897:            IF (FOUND) THEN 
898:            PARTCTR = PARTCTR + 1 
899:            ENDIF 
900:   
901:            
902:  
903:             I = LASRES 
904:  
905:           FOUND=.FALSE. 
906:            PARTATOMS=0 
907:            DO J= (IBASE(I)+1), (IBASE(I)+6) 
908:            IF(FLAGS(J).NE.0) THEN 
909:            FOUND=.TRUE. 
910:            FLAGS(J)=QNPART+PARTCTR 
911:            PARTATOMS=PARTATOMS+1 
912:            ENDIF 
913:            ENDDO 
914:  
915:            IF (FOUND) THEN 
916:            PARTCTR = PARTCTR + 1 
917:            ENDIF 
918:  
919:            FOUND=.FALSE. 
920:            PARTATOMS=0 
921:            DO J = (IBASE(I)+7), (IBASE(I+1)-6) 
922:            IF(FLAGS(J).NE.0) THEN 
923:            FOUND=.TRUE. 
924:            FLAGS(J)=QNPART+PARTCTR 
925:            PARTATOMS=PARTATOMS+1 
926:            ENDIF 
927:            ENDDO 
928:            IF (FOUND) THEN 
929:            PARTCTR = PARTCTR + 1 
930:            ENDIF 
931:  
932:            FOUND=.FALSE. 
933:            PARTATOMS=0 
934:            DO J = (IBASE(I+1)-5), (IBASE(I+1)) 
935:            IF(FLAGS(J).NE.0) THEN 
936:            FOUND=.TRUE. 
937:            FLAGS(J)=QNPART+PARTCTR 
938:            PARTATOMS=PARTATOMS+1 
939:            ENDIF 
940:            ENDDO 
941:            IF (FOUND) THEN 
942:            PARTCTR = PARTCTR + 1 
943:            ENDIF 
944:  
945:            ELSE 
946:            ENDIF ! CONDITIONAL on FIRTYP 
947:  
948:           QNPART = QNPART + PARTCTR - 1 
949:  
950:           ELSE 
951:           ENDIF  ! CONDITIONAL on NUMRES 
952:  
953:           RETURN 
954:           END           
955:  
956:  
957: C        DO I = 1,NRES,NRES-1 
958: C           write(*,*) " I , RTRTYP = ", I,RTRTYP(I) 
959: C           FOUND=.FALSE. 
960: C           PARTATOMS=0 
961: C           DO J=(IBASE(I)+1),(IBASE(I+1)) 
962: C             IF(FLAGS(J).NE.0) THEN 
963: C             FOUND=.TRUE. 
964: C             FLAGS(J)=QNPART+PARTCTR 
965: C             PARTATOMS=PARTATOMS+1 
966: C             ENDIF 
967: C           ENDDO 
968: C           IF (FOUND) THEN 
969: C           PARTCTR=PARTCTR+1 
970: C           ENDIF 
971: C        ENDDO 
972: C 
973: C         QNPART=QNPART+PARTCTR-1 
974: C 
975: C           write(*,*)  " QNPART = ", QNPART 
976: C 
977: C       ELSEIF(NRES.EQ.1) THEN 
978: C 
979: C       write(*,*) " NRES, RTRTYP = ", NRES,RTRTYP(NRES) 
980: CC       DO I = 1,NRES 
981: C       DO J = IBASE(I)+1, IBASE(I+1) 
982: C        write(*,*) " IAC, TYPE = ", IAC(J),TYPE(J) 
983: C       ENDDO 
984: C       ENDDO 
985: C       write(*,*) " QNPART = ", QNPART 
986: C 
987: C       ELSEIF(NRES.EQ.2) THEN 
988: C 
989: C       write(*,*) " NRES = ", NRES 
990: C 
991: C 
992: C       write(*,*) " QNPART = ", QNPART 
993: C 
994: C       ENDIF 
995:  
996:          
997: C        IF (FOUND) THEN 
998: C          RESATOMS=IBASE(I+1)-IBASE(I) 
999: C          IF(RESATOMS.GT.PARTATOMS) THEN 
1000: C            CALL WRNDIE(0,'<FQCOM>','All atoms of residue not selected') 
1001: C          ENDIF 
1002: C          PARTCTR=PARTCTR+1 
1003: C       ENDIF 
1004: C      ENDDO 
1005: CCC      QNPART=QNPART+PARTCTR-1 
1006: CCC 
1007: CCC      RETURN 
1008: CCC      END 
1009: C 
1010: C----------------------------------------------------------------------- 
1011:  
1012:       SUBROUTINE PARTMERGE(QPRTBIN,NEWPART,NATOM) 
1013: C This subroutine is called from QPRTN and merges the new parts of 
1014: C FLAGS (i.e. nonzero values) into the QPRTBIN array 
1015: C 
1016: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1017: C 
1018:       INTEGER NEWPART(*),QPRTBIN(*),NATOM 
1019:       INTEGER I 
1020: C 
1021:       DO I=1,NATOM 
1022:        IF(NEWPART(I).NE.0) QPRTBIN(I)=NEWPART(I) 
1023:       ENDDO 
1024:       RETURN 
1025:       END 
1026: C 
1027: C----------------------------------------------------------------------- 
1028: c zz12-02-95 
1029:       SUBROUTINE MOLSRT(NATOM,MOLT,K,MOLPTR,MOLBL,IB,JB,NBOND, 
1030:      &    QCGRP,NGRP,IGPBS,BLIST) 
1031: C 
1032: C This routine sorts out molecules defined by bonds. Atoms with same 
1033: C MOLT are linked to each other through chemical bonds 
1034: C 
1035: C  Zhou 12-02-95 
1036: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1037: ##INCLUDE '~/charmm_fcm/stream.fcm' 
1038:       INTEGER NATOM,IB(*),JB(*),NBOND,MOLT(*),K,BLIST(13,*),MOLPTR(*), 
1039:      $      MOLBL(*),IGPBS(*),NGRP 
1040:       LOGICAL QCGRP 
1041: C  Local variable 
1042:       INTEGER I,J,M,IS,IQ 
1043:       IF (QCGRP) THEN 
1044: C use groups as molecules 
1045:          M=1 
1046:          DO I=1,NGRP 
1047:            MOLPTR(I)=M 
1048:            IS=IGPBS(I)+1 
1049:            IQ=IGPBS(I+1) 
1050:            DO J=IS,IQ 
1051:              MOLT(J)=I 
1052:              MOLBL(M)=J 
1053:              M=M+1 
1054:            ENDDO 
1055:          ENDDO 
1056:          K=NGRP 
1057: C This is to make sure the last molecule has the right # of atoms 
1058:          MOLPTR(K+1)=M 
1059:          IF(M.NE.(NATOM+1)) CALL WRNDIE(-3,'MOLSRT>', 
1060:      $      'Some atoms have not assigned a mol. #.') 
1061: C test printout 
1062: C         WRITE(OUTU,'(A)')"I,MOLT,MOLPTR,MOLBL:" 
1063: C         DO I=1,NATOM 
1064: C           WRITE(OUTU,'(4(2x,I9))')I,MOLT(I),MOLPTR(I),MOLBL(I) 
1065: C         ENDDO 
1066:          RETURN 
1067:       ENDIF 
1068:       DO I=1,NATOM 
1069:          BLIST(1,I)=1 
1070:          MOLPTR(I)=0 
1071:          MOLBL(I)=0 
1072:       ENDDO 
1073:       DO I=1,NBOND 
1074:          IF (IB(I).GT.JB(I)) THEN 
1075:             J=JB(I) 
1076:             K=IB(I) 
1077:          ELSE 
1078:             J=IB(I) 
1079:             K=JB(I) 
1080:          ENDIF 
1081:          IF ((J.LT.1) .OR. (J.GT.NATOM) .OR. (K.LT.1) .OR. (K.GT.NATOM)) 
1082:      $        GOTO 21 
1083:          BLIST(1,J)=BLIST(1,J)+1 
1084:          IF (BLIST(1,J).GT.13) THEN 
1085:             WRITE(OUTU,900)J,(BLIST(M,J),M=2,13),K 
1086:  900        FORMAT(/,1X,'<MOLSRT>: More than 12 atoms bonded to atom ', 
1087:      $       I8,':',/,10X,7(2X,I8),/) 
1088:             CALL WRNDIE(-3,'MOLSRT>', 
1089:      $      'More than 12 bonds found for an atom') 
1090:          ENDIF 
1091:          BLIST(BLIST(1,J),J)=K 
1092: 21       CONTINUE 
1093:       ENDDO 
1094: C test printout 
1095: C      DO I=1,NATOM 
1096: C        WRITE(OUTU,'(7(2X,I5))')I,(BLIST(J,I),J=2,BLIST(1,I)) 
1097: C      ENDDO 
1098:       DO I=NATOM,1,-1 
1099:          MOLT(I)=I 
1100:          IF (BLIST(1,I).GT.1) THEN 
1101:             DO J=2,BLIST(1,I) 
1102:                M=MOLT(BLIST(J,I)) 
1103:                DO K=I+1,NATOM 
1104:                   IF (MOLT(K).EQ.M) MOLT(K)=I 
1105:                ENDDO 
1106:             ENDDO 
1107:          ENDIF 
1108:       ENDDO 
1109: C eliminate gaps between molecular labels 
1110:       BLIST(1,1)=1 
1111:       K=2 
1112:       DO I=2,NATOM 
1113:          IF (MOLT(I).GT.MOLT(I-1)) THEN 
1114:             BLIST(1,MOLT(I))=K 
1115:             K=K+1 
1116:           ENDIF 
1117:       ENDDO 
1118:       K=K-1 
1119:       DO I=1,NATOM 
1120:          MOLT(I)=BLIST(1,MOLT(I)) 
1121:       ENDDO 
1122: C Rearrange atom so all atoms with the same MOLT will be in consecutive 
1123: C block in MOLBL. MOLPTR points to the first atom of a molecule in MOLBL 
1124:       M=1 
1125:       DO I=1,K 
1126:          MOLPTR(I)=M 
1127:          DO J=1,NATOM 
1128:             IF (MOLT(J).EQ.I) THEN 
1129:                MOLBL(M)=J 
1130:                M=M+1 
1131:             ENDIF 
1132:          ENDDO 
1133:       ENDDO 
1134: C This is to make sure the last molecule has the right # of atoms 
1135:       MOLPTR(K+1)=M 
1136:       IF(M.NE.(NATOM+1)) CALL WRNDIE(-3,'MOLSRT>', 
1137:      $   'Some atoms have not assigned a mol. #.') 
1138: C test printout 
1139: C      WRITE(OUTU,'(/,A,I8,/,A)')'MOLNT = ',K,'I,MOLT,MOLPTR,MOLBL' 
1140: C      WRITE(OUTU,'(4(2X,I9))')(I,MOLT(I),MOLPTR(I),MOLBL(I),I=1,NATOM) 
1141: C 
1142:       RETURN 
1143:       END 
1144: C----------------------------------------------------------------------- 
1145:        SUBROUTINE TYPMOL(MODE,FLAGS) 
1146: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1147: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1148: ##INCLUDE '~/charmm_fcm/param.fcm' 
1149: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
1150: ##INCLUDE '~/charmm_fcm/stack.fcm' 
1151: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1152: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1153: ##INCLUDE '~/charmm_fcm/comand.fcm' 
1154: ##INCLUDE '~/charmm_fcm/coord.fcm' 
1155:        CHARACTER*4 MODE 
1156:        INTEGER I,IW,OLDUSD,MARK,J,J1,JMD,JMT,K 
1157:        INTEGER IWN(5),FLAGS(*) 
1158:  
1159:        write(*,*) 'TYPMOL: Setting selection to type ',MODE 
1160:        DO I=1,NATOM 
1161:          IF (FLAGS(I).GT.0) THEN  
1162:            MOLTYPE(I)=0 
1163:          ENDIF 
1164:        ENDDO 
1165:        IF (MODE.EQ.'WATE') THEN 
1166:           MARK=-1 
1167:        ELSEIF (MODE.EQ.'SPC ') THEN 
1168:           MARK=1 
1169:        ELSEIF (MODE.EQ.'TIP4') THEN 
1170:           MARK=3 
1171:        ELSE   
1172:          RETURN 
1173:        END IF 
1174: C Just in case 
1175:        OLDUSD=LSTUSD 
1176:        IW=ALLSTK(INTEG4(NATOM)) 
1177:        CALL SELCTA(COMLYN,COMLEN,STACK(IW),X,Y,Z,WMAIN,.FALSE.) 
1178:  
1179:        IF (MODE.EQ.'SPC') THEN 
1180:  
1181:        DO I=1,NATOM 
1182:          IF(FLAGS(I).GT.0) THEN 
1183:            IF ((STACK(IW+I-1).EQ.1).AND.(MOLTYPE(I).EQ.0)) THEN 
1184:              JMT=MOLT(I) 
1185:              IF ((MOLPTR(JMT+1)-MOLPTR(JMT)).EQ.3) THEN  !must be water! 
1186:                K=1 
1187:                DO JMD=MOLPTR(JMT),MOLPTR(JMT+1)-1 
1188:                  J=MOLBL(JMD) 
1189:                  IWN(K)=J 
1190:                  K=K+1 
1191:                ENDDO 
1192:                IF (ITC(IAC(IWN(1))).EQ.ITC(IAC(IWN(2)))) THEN 
1193:                  MOLTYPE(IWN(1))=2*MARK 
1194:                  MOLTYPE(IWN(2))=2*MARK 
1195:                  MOLTYPE(IWN(3))=MARK 
1196:                ELSE IF (ITC(IAC(IWN(1))).EQ.ITC(IAC(IWN(3)))) THEN 
1197:                  MOLTYPE(IWN(1))=2*MARK 
1198:                  MOLTYPE(IWN(2))=MARK 
1199:                  MOLTYPE(IWN(3))=2*MARK 
1200:                ELSE IF (ITC(IAC(IWN(2))).EQ.ITC(IAC(IWN(3)))) THEN 
1201:                  MOLTYPE(IWN(1))=MARK 
1202:                  MOLTYPE(IWN(2))=2*MARK 
1203:                  MOLTYPE(IWN(3))=2*MARK 
1204:                ELSE 
1205:                  CALLWRNDIE(-3,'<MOLTYP>','It is not Water.') 
1206:                END IF 
1207:              END IF 
1208:            END IF 
1209:          END IF 
1210:        ENDDO 
1211:        
1212:       ELSEIF (MODE.EQ.'TIP4') THEN 
1213:  
1214:        DO I=1,NATOM 
1215:          IF(FLAGS(I).GT.0) THEN 
1216:            IF ((STACK(IW+I-1).EQ.1).AND.(MOLTYPE(I).EQ.0)) THEN 
1217:              JMT=MOLT(I) 
1218:              IF ((MOLPTR(JMT+1)-MOLPTR(JMT)).EQ.4) THEN  !must be 4-point water! 
1219:                K=1 
1220:                DO JMD=MOLPTR(JMT),MOLPTR(JMT+1)-1 
1221:                  J=MOLBL(JMD) 
1222:                  IWN(K)=J 
1223:                  K=K+1 
1224:                ENDDO 
1225: c       assume tip4p with following order of atoms:  OH2, OM, H1, H2 
1226: c       OM is the lone pair which carries the charge and undergoes 
1227: c       fluctuating charge dynamics 
1228: c       hardcoded for now 
1229: c       S. Patel: (08/28/02) 
1230:  
1231:                   MOLTYPE(IWN(1))=9 
1232:                   MOLTYPE(IWN(2))=3 
1233:                   MOLTYPE(IWN(3))=4 
1234:                   MOLTYPE(IWN(4))=4 
1235:  
1236: C               IF (ITC(IAC(IWN(1))).EQ.ITC(IAC(IWN(2)))) THEN 
1237: C                 MOLTYPE(IWN(1))=2*MARK 
1238: C                 MOLTYPE(IWN(2))=2*MARK 
1239: C                 MOLTYPE(IWN(3))=MARK 
1240: C               ELSE IF (ITC(IAC(IWN(1))).EQ.ITC(IAC(IWN(3)))) THEN 
1241: C                 MOLTYPE(IWN(1))=2*MARK 
1242: C                 MOLTYPE(IWN(2))=MARK 
1243: C                 MOLTYPE(IWN(3))=2*MARK 
1244: C               ELSE IF (ITC(IAC(IWN(2))).EQ.ITC(IAC(IWN(3)))) THEN 
1245: C                 MOLTYPE(IWN(1))=MARK 
1246: C                 MOLTYPE(IWN(2))=2*MARK 
1247: C                 MOLTYPE(IWN(3))=2*MARK 
1248: C               ELSE 
1249: C                 CALLWRNDIE(-3,'<MOLTYP>','It is not Water.') 
1250: C               END IF 
1251:              END IF 
1252:            END IF 
1253:          END IF 
1254:        ENDDO 
1255:  
1256:        ELSE 
1257:        ENDIF 
1258:  
1259:        LSTUSD=OLDUSD 
1260:        RETURN 
1261:        END 
1262: C----------------------------------------------------------------------- 
1263:       SUBROUTINE FILLCGREF(FLAGS) 
1264: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1265: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1266: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1267: ##INCLUDE '~/charmm_fcm/coord.fcm' 
1268: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1269:  
1270:       INTEGER FLAGS(*) 
1271:       INTEGER I 
1272:  
1273:       DO I=1,NATOM 
1274:        IF(FLAGS(I).NE.0) CGREF(I)=WMAIN(I) 
1275:       ENDDO 
1276:  
1277:       RETURN 
1278:       END 
1279: C--------------------------------------------------------- 
1280:       SUBROUTINE ASSIGNQMAS(FLAGS,QMA,TSTA,PMASSQ) 
1281: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1282: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1283: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1284: ##INCLUDE '~/charmm_fcm/coord.fcm' 
1285: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1286: ##INCLUDE '~/charmm_fcm/consta.fcm' 
1287:  
1288:        INTEGER FLAGS(*) 
1289:        INTEGER I 
1290:        REAL*8  QMA,TSTA,PMASSQ(NATOM) 
1291:  
1292:        DO I=1,NATOM 
1293:           IF(FLAGS(I).NE.0) THEN  
1294:               PMASSQ(I)=QMA/(TIMFAC*TIMFAC) 
1295:           ENDIF 
1296:        ENDDO 
1297:        RETURN 
1298:        END 
1299: C----------------------------------------------------------------------- 
1300:       SUBROUTINE ASSWALLPP(FLAGS,QRA1,QRB1,QRA2,QRB2,QRQ1,QRQ2,QRK, 
1301:      $                     wallpa1,wallpb1,wallpa2,wallpb2,wallpq1, 
1302:      $                     wallpq2,wallpK,NATOM) 
1303: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1304: C 
1305:        INTEGER FLAGS(*) 
1306:        INTEGER I,NATOM 
1307:        REAL*8 QRA1,QRB1,QRA2,QRB2,QRQ1,QRQ2,QRK 
1308:        REAL*8 wallpa1(NATOM), wallpb1(NATOM), wallpa2(NATOM), 
1309:      $        wallpb2(NATOM), wallpq1(NATOM), wallpq2(NATOM), 
1310:      $        wallpK(NATOM) 
1311: C 
1312: C 
1313:        DO I=1,NATOM 
1314:           IF(FLAGS(I).NE.0) THEN 
1315:           wallpa1(I)=QRA1 
1316:           wallpb1(I)=QRB1 
1317:           wallpa2(I)=QRA2 
1318:           wallpb2(I)=QRB2 
1319:           wallpQ1(I)=QRQ1 
1320:           wallpQ2(I)=QRQ2 
1321:           wallpK(I)=QRK 
1322:           ENDIF 
1323:        ENDDO 
1324:       RETURN 
1325:       END 
1326:  
1327:  
1328:  
1329: C -- 
1330:  
1331:       SUBROUTINE CHEQRESETNORM 
1332:  
1333: C     reset the QPARTED array used in checking for assignmetn of atoms 
1334: c     to more than one normalization unit 
1335:  
1336: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1337: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1338: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
1339: ##INCLUDE '~/charmm_fcm/number.fcm' 
1340: ##INCLUDE '~/charmm_fcm/coord.fcm' 
1341: ##INCLUDE '~/charmm_fcm/stream.fcm' 
1342: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1343: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1344:  
1345:  
1346:        INTEGER I 
1347:  
1348:         DO I = 1, NATOM 
1349:           QPARTED(I)=.FALSE. 
1350:         ENDDO 
1351:  
1352: ##ELSE (cheq_main) 
1353:       SUBROUTINE NULL_FQCOM 
1354: ##ENDIF (cheq_main) 
1355:       RETURN 
1356:       END 


r21974/fqener.src 2017-01-21 10:33:58.058250000 +0000 r21973/fqener.src 2017-01-21 10:34:01.978250000 +0000
  1: CHARMM Element source/flucq/fqene.src 0.0  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/source/cheq/fqener.src' in revision 21973
  2: ##IF CHEQ (cheq_main) 
  3: C 
  4:       SUBROUTINE QCGSET(NATOM,X,Y,Z,DX,DY,DZ,CG,INBLO,IBLO14,INB14,JNB, 
  5:      &                  IAC,LELEC,OUTU) 
  6: C----------------------------------------------------------------------- 
  7: C  This routine sets the charges, forces and second derivatives at the 
  8: C  beginning of the energy call. 
  9: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
 10: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
 11: ##INCLUDE '~/charmm_fcm/number.fcm' 
 12: C 
 13:       INTEGER NATOM,IAC(*),OUTU,JJ,thisisjunk 
 14:       INTEGER INBLO(*),IBLO14(*),INB14(*),JNB(*) 
 15:       REAL*8 X(*),Y(*),Z(*),DX(*),DY(*),DZ(*),CG(*) 
 16:       LOGICAL LELEC 
 17: C 
 18: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
 19: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
 20: ##INCLUDE '~/charmm_fcm/heap.fcm' 
 21: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
 22: ##INCLUDE '~/charmm_fcm/fast.fcm' 
 23: ##IF PARALLEL 
 24: ##INCLUDE '~/charmm_fcm/parallel.fcm' 
 25: ##ENDIF 
 26: C 
 27:       INTEGER I,ATMP,DATMP,NDA 
 28:       INTEGER ETASIZE 
 29:       INTEGER inter,Rpoint,delR,temp 
 30:  
 31:       IF (LELEC) THEN 
 32:          IF (QCGPOL) THEN 
 33:           write(OUTU,*)" CALING POLT" 
 34:           ETASIZE=NATOM*NATOM 
 35:             ATMP=ALLHP(IREAL8(ETASIZE),'fqener.src','QCGSET','ATMP') 
 36:              NDA=1 
 37:           DATMP=ALLHP(IREAL8(NDA*NDA),'fqener.src','QCGSET','DATMP') 
 38:  
 39:           inter =ALLHP(IREAL8(3*NATOM),'fqener.src','QCGSET','inter') 
 40:           Rpoint=ALLHP(IREAL8(3*NATOM),'fqener.src','QCGSET','Rpoint') 
 41:           delR  =ALLHP(IREAL8(3*NATOM),'fqener.src','QCGSET','delR') 
 42:           temp  =ALLHP(IREAL8(3*NATOM),'fqener.src','QCGSET','temp') 
 43:  
 44:           CALL FQPOLTEN(NATOM,X,Y,Z,DX,DY,DZ,CG,ECH,EHA,QMOLC,MOLT, 
 45:      $  HEAP(ATMP),HEAP(DATMP),DCH,QCGINVF,JNB, 
 46:      $  INBLO,INB14,IBLO14,IAC,QPARTBIN,MOLTYPE,QPRNETA, 
 47:      $  HEAP(ICCNBA),HEAP(ICCNBB),HEAP(ICCNBC),HEAP(ICCNBD),IACNB, 
 48:      $           NITCC2,LOWTP,OUTU, 
 49:      &  HEAP(inter),HEAP(Rpoint),HEAP(delR),HEAP(temp)) 
 50:  
 51:       CALL FREHP(inter, IREAL8(3*NATOM),'fqener.src','QCGSET','inter') 
 52:       CALL FREHP(Rpoint,IREAL8(3*NATOM),'fqener.src','QCGSET','Rpoint') 
 53:       CALL FREHP(delR,  IREAL8(3*NATOM),'fqener.src','QCGSET','delR') 
 54:       CALL FREHP(temp,  IREAL8(3*NATOM),'fqener.src','QCGSET','temp') 
 55:  
 56:        ENDIF 
 57:          IF (QCGINV) THEN 
 58:             ETASIZE=NATOM*NATOM 
 59:             ATMP=ALLHP(IREAL8(ETASIZE),'fqener.src','QCGSET','ATMP') 
 60:             IF (QCGINVF) THEN 
 61:                NDA=NATOM 
 62:             ELSE 
 63:                NDA=1 
 64:             ENDIF 
 65:             DATMP=ALLHP(IREAL8(NDA*NDA),'fqener.src','QCGSET','DATMP') 
 66:  
 67:  
 68:             CALL CGINV(NATOM,X,Y,Z,DX,DY,DZ,CG,ECH,EHA,QMOLC,MOLT, 
 69:      $  HEAP(ATMP),HEAP(DATMP),DCH,QCGINVF,JNB, 
 70:      $  INBLO,INB14,IBLO14,IAC,QPARTBIN,MOLTYPE,QPRNETA, 
 71:      $  HEAP(ICCNBA),HEAP(ICCNBB),HEAP(ICCNBC),HEAP(ICCNBD),IACNB, 
 72:      $           NITCC2,LOWTP,OUTU) 
 73:  
 74:             DO I=1,NATOM 
 75:                CG(I)=DCH(I) 
 76:                DCH(I)=ECH(I) 
 77:                DDCH(I)=ZERO 
 78:             ENDDO 
 79:             CALL FREHP(ATMP,IREAL8(ETASIZE),'fqener.src','QCGSET', 
 80:      $                'ATMP') 
 81:             CALL FREHP(DATMP,IREAL8(NDA*NDA),'fqener.src','QCGSET', 
 82:      $                'DATMP') 
 83:             QCGINV=.FALSE. 
 84:          ELSE 
 85:             DO I=1,NATOM 
 86:                DCH(I)=ECH(I) 
 87:                DDCH(I)=ZERO 
 88:             ENDDO 
 89:          ENDIF 
 90:       ELSE 
 91:          DO I=1,NATOM 
 92:             DCH(I)=ZERO 
 93:             DDCH(I)=ZERO 
 94:          ENDDO 
 95:       ENDIF 
 96:  
 97:  
 98:       RETURN 
 99:       END !SUB QCGSET 
100:  
101:       SUBROUTINE QNORM(DCH,QNPART,QPRTPTR,QPART,PMASSQ, 
102:      $                 MINNORM,QT4P) 
103: C----------------------------------------------------------------------- 
104: C  This subroutine keeps the charge for the given partitions constant 
105: C  over a partition by setting the net derivative with respect to 
106: C  charge to zero. 
107: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
108: ##INCLUDE '~/charmm_fcm/number.fcm' 
109:  
110:       INTEGER QNPART,QPRTPTR(*),QPART(*) 
111:       REAL*8  DCH(*),PMASSQ(*) 
112:  
113:       INTEGER I,J,NAT(10),frstatom,IFLAG 
114:       REAL*8 PARTDCH,TMASS,MASINV 
115:       LOGICAL MINNORM,QT4P 
116:  
117:  
118:            IFLAG=2 
119:          IF (IFLAG.EQ.1) THEN 
120:           NAT(1)=4 
121:           NAT(2)=2 
122:           NAT(3)=4 
123:           NAT(4)=2 
124:           NAT(5)=4 
125:           NAT(6)=4 
126:           NAT(7)=2 
127:           NAT(8)=4 
128:           NAT(9)=4 
129:           NAT(10)=3 
130:  
131:        frstatom=0 
132:        DO I=1,10 
133:            PARTDCH=ZERO 
134:            DO J=frstatom+1,frstatom+NAT(I) 
135:              PARTDCH=PARTDCH + DCH(J) 
136:            ENDDO 
137:            PARTDCH=PARTDCH/(NAT(I)) 
138:            DO J=frstatom+1,frstatom+NAT(I) 
139:             DCH(J)= DCH(J)-PARTDCH 
140:            ENDDO 
141:             frstatom= J - 1 
142:           ENDDO 
143:  
144:           ELSEIF(IFLAG.EQ.2) THEN 
145:  
146:  
147:        IF (MINNORM) THEN 
148:  
149:  
150:  
151:             IF (QT4P ) THEN 
152:  
153:         DO I=1,QNPART 
154:         PARTDCH=ZERO 
155:  
156:         IF ( ((QPRTPTR(I+1)-QPRTPTR(I)).EQ.4) ) THEN ! .AND. 
157: C     2      (MOLTYPE(QPRTPTR(I)+1).EQ.9) ) THEN 
158:  
159:         DO J=QPRTPTR(I)+2,QPRTPTR(I+1) 
160:           PARTDCH=PARTDCH+DCH(QPART(J)) 
161:        ENDDO !J 
162:         PARTDCH=PARTDCH/(QPRTPTR(I+1)-QPRTPTR(I)-1) 
163: c      Write(6,*)' Correction for partition ', i, ' = ', partdch 
164:         DO J=QPRTPTR(I)+2,QPRTPTR(I+1) 
165:           DCH(QPART(J))=DCH(QPART(J))-PARTDCH 
166:         ENDDO !J 
167:  
168:       ELSE 
169:  
170:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
171:           PARTDCH=PARTDCH+DCH(QPART(J)) 
172:        ENDDO !J 
173:         PARTDCH=PARTDCH/(QPRTPTR(I+1)-QPRTPTR(I)) 
174: c      Write(6,*)' Correction for partition ', i, ' = ', partdch 
175:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
176:           DCH(QPART(J))=DCH(QPART(J))-PARTDCH 
177:         ENDDO !J 
178:  
179:        ENDIF 
180:  
181:       ENDDO !I 
182:  
183:  
184:  
185:         ELSE 
186:  
187:         DO I=1,QNPART 
188:         PARTDCH=ZERO 
189:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
190:           PARTDCH=PARTDCH+DCH(QPART(J)) 
191:        ENDDO !J 
192:         PARTDCH=PARTDCH/(QPRTPTR(I+1)-QPRTPTR(I)) 
193: c      Write(6,*)' Correction for partition ', i, ' = ', partdch 
194:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
195:           DCH(QPART(J))=DCH(QPART(J))-PARTDCH 
196:         ENDDO !J 
197:       ENDDO !I 
198:  
199:  
200:          ENDIF  ! CONDITIONAL ON QT4P 
201:  
202:        ELSE 
203:  
204:  
205:  
206:         IF (QT4P) THEN 
207:  
208:       DO I=1,QNPART 
209:         TMASS=ZERO 
210:         PARTDCH=ZERO 
211:  
212:         IF ( ((QPRTPTR(I+1)-QPRTPTR(I)).EQ.4)) THEN  !   .AND. 
213: C     2      (MOLTYPE(QPRTPTR(I)+1).EQ.9) ) THEN 
214:  
215:         DO J=QPRTPTR(I)+2,QPRTPTR(I+1) ! Starte with OM site on TIP4P hardwired topology 
216:           MASINV=ONE/PMASSQ(QPART(J)) 
217:           PARTDCH=PARTDCH+DCH(QPART(J))*MASINV 
218: C           PARTDCH=PARTDCH+DCH(QPART(J)) 
219:           TMASS = TMASS + MASINV 
220:        ENDDO !J 
221:         PARTDCH=PARTDCH/TMASS 
222:         DO J=QPRTPTR(I)+2,QPRTPTR(I+1) ! Start with OM site on TIP4P hardwired topology 
223:           DCH(QPART(J))=DCH(QPART(J))-PARTDCH 
224:         ENDDO !J 
225:  
226:         ELSE 
227:  
228:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
229:           MASINV=ONE/PMASSQ(QPART(J)) 
230:           PARTDCH=PARTDCH+DCH(QPART(J))*MASINV 
231:           TMASS = TMASS + MASINV 
232:        ENDDO !J 
233:         PARTDCH=PARTDCH/TMASS 
234:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
235:           DCH(QPART(J))=DCH(QPART(J))-PARTDCH 
236:         ENDDO !J 
237:  
238:           ENDIF 
239:  
240:       ENDDO !I 
241:  
242:        ELSE 
243:  
244:       DO I=1,QNPART 
245:         TMASS=ZERO 
246:         PARTDCH=ZERO 
247:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
248:           MASINV=ONE/PMASSQ(QPART(J)) 
249:           PARTDCH=PARTDCH+DCH(QPART(J))*MASINV 
250: C           PARTDCH=PARTDCH+DCH(QPART(J)) 
251:           TMASS = TMASS + MASINV 
252:        ENDDO !J 
253: C        PARTDCH=PARTDCH/(QPRTPTR(I+1)-QPRTPTR(I)) 
254:         PARTDCH=PARTDCH/TMASS 
255:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
256: C          DCH(QPART(J))=DCH(QPART(J))-PMASSQ(QPART(J))*PARTDCH 
257:           DCH(QPART(J))=DCH(QPART(J))-PARTDCH 
258:         ENDDO !J 
259:       ENDDO !I 
260:  
261: C             ENDIF 
262:  
263:        ENDIF  !  CONDITONIAL on Q4TP 
264:  
265:  
266:  
267:        ENDIF  ! MINNORM CONDITIONAL 
268:  
269:         ELSEIF(IFLAG.EQ.3) THEN 
270:  
271: C 
272:  
273:                ENDIF 
274:  
275:       RETURN 
276:       END !SUB QNORM 
277:  
278:       SUBROUTINE QBENER(NATOM,CG,X,Y,Z,DX,DY,DZ,IBLO14,INB14,IAC, 
279:      &                  QSECD,LEWALD,LELEC,pwallpa1,pwallpb1,pwallpa2, 
280:      &                  pwallpb2,pwallpq1,pwallpq2,pwallpK) 
281: C----------------------------------------------------------------------- 
282: C  This subroutine calculates the electrostatic energy contribution of 
283: C  the pairs on the non-bonded exclusion list. 
284: C 
285: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
286: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
287: ##INCLUDE '~/charmm_fcm/exfunc.fcm' 
288: ##INCLUDE '~/charmm_fcm/number.fcm' 
289: ##INCLUDE '~/charmm_fcm/parallel.fcm' 
290: C 
291:       INTEGER NATOM,IAC(*),IBLO14(*),INB14(*) 
292:       REAL*8 CG(*) 
293:       REAL*8 X(*),Y(*),Z(*),DX(*),DY(*),DZ(*) 
294:       LOGICAL QSECD,LEWALD,LELEC 
295: C 
296: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
297: ##INCLUDE '~/charmm_fcm/energy.fcm' 
298: ##INCLUDE '~/charmm_fcm/param.fcm' 
299: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
300: ##INCLUDE '~/charmm_fcm/stream.fcm' 
301: ##INCLUDE '~/charmm_fcm/consta.fcm' 
302:  
303:  
304:       INTEGER I,J,JJ,NBINDX,JCTR,JCOUNT,thisisjunk,KK,KJ 
305:       INTEGER IHIJ,IMASS,JMASS,ISTEP,IFLAG 
306:       REAL*8 HIJ,HIJ1,CGT,FI,CGF,ECSUM,ECSUM1,ECSUM2,ENPB,S 
307:       REAL*8 CRXI,CRYI,CRZI,DXIJ,DYIJ,DZIJ,FXI,FYI,FZI 
308:       REAL*8 FJ,FXJ,FYJ,FZJ,CGTJ,sum1  
309:  
310:       REAL*8 term1,term2,term3,swchfunc,myeterm 
311:       REAL*8 term4,myforce,term2f,BRACK,QGP,ETILD 
312:       REAL*8 B(12) 
313: c 
314:       REAL*8 pwallpa1(NATOM),pwallpb1(NATOM), pwallpa2(NATOM), 
315:      &       pwallpb2(NATOM),pwallpq1(NATOM), pwallpq2(NATOM), 
316:      &       pwallpk(NATOM) 
317:  
318: 433     FORMAT(' QBENER: ',A) 
319: 435     FORMAT(20I5) 
320:  
321:  
322: ##IF PARALLEL 
323:  
324:        IF (MYNOD.EQ.0) THEN 
325: ##ENDIF 
326:  
327:           sum1 = 0.0 
328:  
329:       ECSUM = ZERO 
330:       ECSUM1=ZERO 
331:       ECSUM2 = ZERO 
332: C  main looping structure 
333:       IF(CGMODEL.EQ.0) THEN 
334: C No CG-XYZ mixing 
335:         DO I=1,NATOM 
336:           ECSUM=ECSUM+ECH(I)*CG(I) 
337:           HIJ=EHA(I)+EHA(I) 
338:           ECSUM1=ECSUM1+CG(I)*HIJ*CG(I) 
339:           DCH(I)=DCH(I)+HIJ*CG(I) 
340:         ENDDO 
341:         DO I=1,NATOM-1 
342:           IF (I.EQ.1) THEN 
343:             NBINDX=1 
344:           ELSE 
345:             NBINDX=IBLO14(I-1)+1 
346:           ENDIF 
347:           JCTR=IBLO14(I) 
348:           DO JJ=NBINDX,JCTR 
349:              J=IABS(INB14(JJ)) 
350:             IF (J.GT.0) THEN 
351:               S=(X(I)-X(J))*(X(I)-X(J))+(Y(I)-Y(J))*(Y(I)-Y(J)) 
352:      $         +(Z(I)-Z(J))*(Z(I)-Z(J)) 
353:               HIJ=CCELEC*HIJ/SQRT((1.0+HIJ*HIJ*S)) 
354:               DCH(I)=DCH(I)+HIJ*CG(J) 
355:               DCH(J)=DCH(J)+HIJ*CG(I) 
356:               IF(QREF) THEN 
357:                 DCH(I)=DCH(I)+HIJ*CGREF(J) 
358:                 DCH(J)=DCH(J)+HIJ*CGREF(I) 
359:               ENDIF 
360:             ENDIF 
361:           ENDDO !JJ 
362:         ENDDO ! I 
363:       ELSE IF (CGMODEL.EQ.1)  THEN 
364:         NBINDX=1 
365: C electrochemical potential term and diagonal hardness term, 
366: C hardness term 
367:         DO I=1,NATOM 
368:          IF (QPARTBIN(I).NE.0) THEN 
369:           IF (MOLTYPE(I).LE.0) THEN 
370:             ECSUM=ECSUM+ECH(I)*CG(I) 
371:             HIJ=EHA(I)+EHA(I) 
372:             ECSUM1=ECSUM1+CG(I)*HIJ*CG(I)*CCELEC 
373:                  DCH(I)=DCH(I)+HIJ*CG(I)*CCELEC 
374:  
375:             IF(QREF) THEN 
376:               DCH(I)=DCH(I)+HIJ*CGREF(I) 
377:             ENDIF 
378:           ELSEIF (MOLTYPE(I).GT.0) THEN   ! for water molecules 
379:             IF (MOLTYPE(I).EQ.1) THEN 
380:               HIJ=367.0 
381:             ELSEIF (MOLTYPE(I).EQ.2) THEN 
382:               HIJ=392.2 
383:             ELSEIF (MOLTYPE(I).EQ.3) THEN 
384:               HIJ = 371.6   !  TIP4P OM SITE 
385:             ELSEIF (MOLTYPE(I).EQ.4) THEN 
386:               HIJ = 353.0   ! TIP4P HYDROGEN SITE 
387:              ELSEIF (MOLTYPE(I).EQ.9) THEN 
388:               HIJ = 0.0   !   TIP4P oxygen site (no charge) 
389:             ELSE 
390:             ENDIF 
391:             ECSUM=ECSUM+ECH(I)*CG(I) 
392:             ECSUM1=ECSUM1+CG(I)*HIJ*CG(I) 
393:             DCH(I)=DCH(I)+HIJ*CG(I) 
394:             IF(QREF) THEN 
395:               DCH(I)=DCH(I)+HIJ*CGREF(I) 
396:             ENDIF 
397:           ENDIF 
398:  
399: C___________________________________________________________________________ 
400: c      add in wall potential energy 
401:  
402:        IF (PTYP.EQ.1) THEN 
403: C    Harmonic potential 
404:          IF (CG(I).LT.pwallpQ2(I)) THEN 
405:  
406:          myeterm = pwallpK(I)*(CG(I)-pwallpQ2(I))*(CG(I)-pwallpQ2(I)) 
407:          myforce = 1.0*2.0*pwallpK(I)*(CG(I)-pwallpQ2(I)) 
408:  
409:          ECSUM = ECSUM + myeterm 
410:          DCH(I) = DCH(I) + myforce 
411:  
412:  
413:       ELSE IF ((CG(I).GE.pwallpQ2(I)).and.(CG(I).LE.pwallpQ1(I))) THEN 
414:  
415:          myeterm = 0.0 
416:          myforce = 0.0 
417:  
418:          ECSUM = ECSUM + myeterm 
419:          DCH(I) = DCH(I) + myforce 
420:  
421:          ELSE IF ( CG(I).GT.pwallpQ1(I) ) THEN 
422:  
423:          myeterm = pwallpK(I)*(CG(I)-pwallpQ1(I))*(CG(I)-pwallpQ1(I)) 
424:          myforce = 1.0*2.0*pwallpK(I)*(CG(I)-pwallpQ1(I)) 
425:  
426:          ECSUM = ECSUM + myeterm 
427:          DCH(I) = DCH(I) + myforce 
428:  
429:          ELSE 
430:          ENDIF 
431:  
432:  
433:        ELSE IF (PTYP.EQ.2) THEN 
434:  
435: c---------------------------------------------------------------------- 
436: C  Nth order wall potential 
437:  
438:         IF ((CG(I).ge.pwallpa1(I)).and.(CG(I).lt.pwallpb1(I))) THEN 
439:  
440:         term1=(1.0)/((pwallpa1(I)-pwallpb1(I))**3) 
441:  
442:         term2=(CG(I)-pwallpa1(I))**2 
443:  
444:         term2f=CG(I)-pwallpa1(I) 
445:         term3= 2.0*CG(I) + pwallpa1(I) - 3.0*pwallpb1(I) 
446:  
447:         term4=1.0/((pwallpQ1(I)-CG(I))) 
448:  
449:         swchfunc=term1*term2*term3 
450:         myeterm = pwallpK(I)/(abs(CG(I)-pwallpQ1(I)))**wallpN 
451:         myeterm = swchfunc*myeterm 
452:         ECSUM = ECSUM + myeterm 
453:  
454:         myforce = term2f + term3 + (0.5*wallpN*term2f*term3*term4) 
455:         myforce = myforce*2.0*term1*term2f*pwallpK(I)*(term4**wallpN) 
456:  
457:         DCH(I) = DCH(I) + myforce 
458:  
459:         ELSE IF ((CG(I).ge.pwallpb1(I)).and.(CG(I).le.pwallpQ1(I))) THEN 
460:         myeterm = pwallpK(I)/(abs(CG(I)-pwallpQ1(I)))**wallpN 
461:         ECSUM = ECSUM + myeterm 
462:  
463:         myforce = wallpN*pwallpK(I)/((pwallpQ1(I)-CG(I))**(wallpN+1)) 
464:         DCH(I) = DCH(I) + myforce 
465:  
466:  
467:         ELSE IF (CG(I).gt.pwallpQ1(I)) THEN 
468:  
469:       myeterm=pwallpK(I)/(abs(CG(I)-pwallpQ1(I)))**wallpN 
470:       myforce=-1.0*wallpN*pwallpK(I)/((CG(I)-pwallpQ1(I))**(wallpN+1)) 
471:  
472:         ECSUM = ECSUm + myeterm 
473:         DCH(I) = DCH(I) + myforce 
474:  
475:       ELSE IF ((CG(I).le.pwallpa2(I)).and.(CG(I).gt.pwallpb2(I))) THEN 
476:  
477:         term1=(1.0)/((pwallpa2(I)-pwallpb2(I))**3) 
478:         term2=(CG(I)-pwallpa2(I))**2 
479:         term2f=CG(I)-pwallpa2(I) 
480:         term3= 2.0*CG(I) + pwallpa2(I) - 3.0*pwallpb2(I) 
481:         term4= 1.0/(CG(I)-pwallpQ2(I)) 
482:         swchfunc=term1*term2*term3 
483:         myeterm = pwallpK(I)/(abs(CG(I)-pwallpQ2(I)))**wallpN 
484:         myeterm = swchfunc*myeterm 
485:         ECSUM = ECSUM + myeterm 
486:  
487:         myforce = term2f + term3 - (0.5*wallpN*term2f*term3*term4) 
488:         myforce = myforce*2.0*term1*term2f*pwallpK(I)*(term4**wallpN) 
489:        DCH(I) = DCH(I) + myforce 
490:  
491:        ELSE IF ((CG(I).le.pwallpb2(I)).and.(CG(I).ge.pwallpQ2(I))) THEN 
492:  
493:        myeterm = pwallpK(I)/(abs(CG(I)-pwallpQ2(I)))**wallpN 
494:        ECSUM = ECSUM + myeterm 
495:  
496:        myforce=-1.0*wallpN*pwallpK(I)/((CG(I)-pwallpQ2(I))**(wallpN+1)) 
497:        DCH(I)=DCH(I) + myforce 
498:  
499:        ELSE IF (CG(I).lt.pwallpQ2(I)) THEN 
500:  
501:        myeterm =  pwallpK(I)/(abs(CG(I)-pwallpQ2(I)))**wallpN 
502:        myforce = wallpN*pwallpK(I)/((pwallpQ2(I)-CG(I))**(wallpN+1)) 
503:        ECSUM = ECSUM + myeterm 
504:        DCH(I) = DCH(I) + myforce 
505:  
506:        ELSE 
507:        ENDIF 
508:  
509:        ELSE 
510:        ENDIF 
511:  
512: C--------------------------------------------------------------------- 
513: C_____________________________________________________________________ 
514:          ENDIF 
515:         ENDDO 
516:  
517: C off-diagonal terms for hardness term 
518:         DO I=1,NATOM-1 
519:         IF (QPARTBIN(I).NE.0) THEN 
520: C  General case 
521:         IF (MOLTYPE(I).EQ.0) THEN 
522:           IF (I.EQ.1) THEN 
523:             NBINDX=1 
524:           ELSE 
525:             NBINDX=IBLO14(I-1)+1 
526:           ENDIF 
527:           JCOUNT=0 
528:           JCTR=IBLO14(I) 
529:           CGT=CCELEC*CG(I) 
530:           IF (CGT.NE.ZERO) THEN 
531:             CRXI=X(I) 
532:             CRYI=Y(I) 
533:             CRZI=Z(I) 
534:             DO JJ=NBINDX,JCTR 
535:                J=IABS(INB14(JJ)) 
536:               IF (J.GT.0) THEN 
537:                 JCOUNT=JCOUNT+1 
538:                 CGTJ=CCELEC*CG(J) 
539:                 DXIJ=CRXI-X(J) 
540:                 DYIJ=CRYI-Y(J) 
541:                 DZIJ=CRZI-Z(J) 
542:                 S=DXIJ*DXIJ+DYIJ*DYIJ+DZIJ*DZIJ 
543:                 HIJ1=EHA(I)+EHA(J) 
544:                 HIJ=1.00*HIJ1/SQRT((1.0+HIJ1*HIJ1*S)) ! SCALED OFF DIAG. ETA 
545:                 ECSUM1=ECSUM1+2*CG(I)*HIJ*CG(J)*CCELEC 
546:                 DCH(I)=DCH(I)+HIJ*CG(J)*CCELEC 
547:                 DCH(J)=DCH(J)+HIJ*CG(I)*CCELEC 
548:                 IF(QREF) THEN 
549:                   DCH(I)=DCH(I)+HIJ*CGREF(J) 
550:                   DCH(J)=DCH(J)+HIJ*CGREF(I) 
551:                 ENDIF 
552: C FORCE FROM CHARGE FLUCTUATING ON NUCLEI 
553:                IF (QCGMINE.OR.(.NOT.QCGWATER)) THEN 
554:                   FI=CG(J)*HIJ*HIJ*HIJ 
555:                   DX(I)=DX(I)-FI*DXIJ*CGT 
556:                   DY(I)=DY(I)-FI*DYIJ*CGT 
557:                   DZ(I)=DZ(I)-FI*DZIJ*CGT 
558:                   FJ=CG(I)*HIJ*HIJ*HIJ 
559:                   DX(J)=DX(J)+FJ*DXIJ*CGTJ 
560:                   DY(J)=DY(J)+FJ*DYIJ*CGTJ 
561:                   DZ(J)=DZ(J)+FJ*DZIJ*CGTJ 
562:                 ENDIF 
563:                ELSE 
564:               ENDIF 
565:             ENDDO !JJ 
566:           ELSE 
567:             CRXI=X(I) 
568:             CRYI=Y(I) 
569:             CRZI=Z(I) 
570:             DO JJ=NBINDX,JCTR 
571:                J=IABS(INB14(JJ)) 
572:               IF (J.GT.0) THEN 
573:                 CGTJ=CCELEC*CG(J) 
574:                 DXIJ=CRXI-X(J) 
575:                 DYIJ=CRYI-Y(J) 
576:                 DZIJ=CRZI-Z(J) 
577:                 S=DXIJ*DXIJ+DYIJ*DYIJ+DZIJ*DZIJ 
578:                 HIJ=EHA(I)+EHA(J) 
579:                 HIJ=HIJ/SQRT((1.0+HIJ*HIJ*S)) 
580:                 DCH(I)=DCH(I)+HIJ*CG(J)*CCELEC 
581:                 IF(QREF) THEN 
582:                   DCH(I)=DCH(I)+HIJ*CGREF(J) 
583:                   DCH(J)=DCH(J)+HIJ*CGREF(I) 
584:                 ENDIF 
585:                 IF (QCGMINE.OR.(.NOT.QCGWATER).AND.(CGTJ.NE.ZERO)) THEN 
586:                   FJ=CG(I)*HIJ*HIJ*HIJ 
587:                   DX(J)=DX(J)+FJ*DXIJ*CGTJ 
588:                   DY(J)=DY(J)+FJ*DYIJ*CGTJ 
589:                   DZ(J)=DZ(J)+FJ*DZIJ*CGTJ 
590:                 ENDIF 
591:               ENDIF 
592:             ENDDO !JJ 
593:           ENDIF !(CGT.NE.ZERO) 
594:         ELSE IF ( (MOLTYPE(I).GT.0).AND.(MOLTYPE(I).LE.2)) THEN  ! for SPC water molecules 
595:           CGT=CCELEC*CG(I) 
596:           IF (I.EQ.1) THEN 
597:             NBINDX=1 
598:           ELSE 
599:             NBINDX=IBLO14(I-1)+1 
600:           ENDIF 
601:           JCTR=IBLO14(I) 
602:           IF (CGT.NE.ZERO) THEN 
603: C            CRXI=X(I) 
604: C            CRYI=Y(I) 
605: C            CRZI=Z(I) 
606:             DO JJ=NBINDX,JCTR 
607:                J=IABS(INB14(JJ)) 
608:               IF (J.GT.0) THEN 
609:                 IF (MOLTYPE(I).EQ.1.OR.MOLTYPE(J).EQ.1) THEN 
610:                    HIJ = 276.0 
611:                 ELSE 
612:                     HIJ = 196.0 
613:                 ENDIF 
614:                 ECSUM1=ECSUM1+2*CG(I)*HIJ*CG(J) 
615:                 DCH(I)=DCH(I)+HIJ*CG(J) 
616:                 DCH(J)=DCH(J)+HIJ*CG(I) 
617:                 IF(QREF) THEN 
618:                   DCH(I)=DCH(I)+HIJ*CGREF(J) 
619:                   DCH(J)=DCH(J)+HIJ*CGREF(I) 
620:                 ENDIF 
621:               ENDIF 
622:             ENDDO !JJ 
623:             NBINDX=NBINDX+IBLO14(I) 
624:           ELSE 
625:             DO JJ=NBINDX,JCTR 
626:                J=IABS(INB14(JJ)) 
627:               IF (J.GT.0) THEN 
628:                 S=(X(I)-X(J))*(X(I)-X(J))+(Y(I)-Y(J))*(Y(I)-Y(J)) 
629:      $         +(Z(I)-Z(J))*(Z(I)-Z(J)) 
630: C The type stuff is for SPC test only 
631:                 IF (MOLTYPE(I).EQ.1) THEN 
632:                   HIJ=276.0 
633:                 ELSE IF (ITC(IAC(J)).EQ.ITC(IAC(1))) THEN 
634:                   HIJ=276.0 
635:                 ELSE 
636:                   HIJ=196.0 
637:                 ENDIF 
638:                 ECSUM1=ECSUM1+2*CG(I)*HIJ*CG(J) 
639:                 DCH(I)=DCH(I)+HIJ*CG(J) 
640:                 IF(QREF) THEN 
641:                   DCH(I)=DCH(I)+HIJ*CGREF(J) 
642:                   DCH(J)=DCH(J)+HIJ*CGREF(I) 
643:                 ENDIF 
644:               ENDIF 
645:             ENDDO !JJ 
646:           ENDIF 
647: c 
648: c  TIP4P STUFF 
649:         ELSE IF ( (MOLTYPE(I).GT.2).AND.(MOLTYPE(I).LE.9)) THEN  ! for TIP4P water molecules 
650:           CGT=CCELEC*CG(I) 
651:           IF (I.EQ.1) THEN 
652:             NBINDX=1 
653:           ELSE 
654:             NBINDX=IBLO14(I-1)+1 
655:           ENDIF 
656:           JCTR=IBLO14(I) 
657:           IF (CGT.NE.ZERO) THEN 
658: C            CRXI=X(I) 
659: C            CRYI=Y(I) 
660: C            CRZI=Z(I) 
661:             DO JJ=NBINDX,JCTR 
662:                J=IABS(INB14(JJ)) 
663:               IF (J.GT.0) THEN 
664:  
665:                 IF (MOLTYPE(J).NE.9) THEN 
666:                 IF (MOLTYPE(I).EQ.3.OR.MOLTYPE(J).EQ.3) THEN 
667:                    HIJ = 286.4 
668:                 ELSE 
669:                     HIJ = 203.6 
670:                 ENDIF 
671:                 ECSUM1=ECSUM1+2*CG(I)*HIJ*CG(J) 
672:                 DCH(I)=DCH(I)+HIJ*CG(J) 
673:                 DCH(J)=DCH(J)+HIJ*CG(I) 
674:                 IF(QREF) THEN 
675:                   DCH(I)=DCH(I)+HIJ*CGREF(J) 
676:                   DCH(J)=DCH(J)+HIJ*CGREF(I) 
677:                 ENDIF 
678:               ENDIF 
679:             ENDIF 
680:             ENDDO !JJ 
681:             ENDIF 
682:  
683:         ELSE ! Other Water 
684:           IF (I.EQ.1) THEN 
685:             NBINDX=1 
686:           ELSE 
687:             NBINDX=IBLO14(I-1)+1 
688:           ENDIF 
689:           JCTR=IBLO14(I) 
690:           CGT=CCELEC*CG(I) 
691:           IF (CGT.NE.ZERO) THEN 
692:             CRXI=X(I) 
693:             CRYI=Y(I) 
694:             CRZI=Z(I) 
695:             DO JJ=NBINDX,JCTR 
696:                J=IABS(INB14(JJ)) 
697:               IF (J.GT.0) THEN 
698:                 DXIJ=CRXI-X(J) 
699:                 DYIJ=CRYI-Y(J) 
700:                 DZIJ=CRZI-Z(J) 
701:                 S=DXIJ*DXIJ+DYIJ*DYIJ+DZIJ*DZIJ 
702:                 HIJ1=EHA(I)+EHA(J) 
703:                 HIJ=HIJ1/SQRT((1.0+HIJ1*HIJ1*S)) 
704:                 ECSUM1=ECSUM1+2*CG(I)*HIJ*CG(J)*CCELEC 
705:                 DCH(I)=DCH(I)+HIJ*CG(J)*CCELEC 
706:                 DCH(J)=DCH(J)+HIJ*CG(I)*CCELEC 
707:                 IF(QREF) THEN 
708:                   DCH(I)=DCH(I)+HIJ*CGREF(J)*CCELEC 
709:                   DCH(J)=DCH(J)+HIJ*CGREF(I)*CCELEC 
710:                 ENDIF 
711:               ENDIF 
712:             ENDDO !JJ 
713:           ELSE 
714:             DO JJ=NBINDX,JCTR 
715:                J=IABS(INB14(JJ)) 
716:               IF (J.GT.0) THEN 
717:                 S=(X(I)-X(J))*(X(I)-X(J))+(Y(I)-Y(J))*(Y(I)-Y(J)) 
718:      $         +(Z(I)-Z(J))*(Z(I)-Z(J)) 
719:                 HIJ=EHA(I)+EHA(J) 
720:                 HIJ=HIJ/SQRT((1.0+HIJ*HIJ*S)) 
721:                 ECSUM1=ECSUM1+2*CG(I)*HIJ*CG(J)*CCELEC 
722:                 DCH(I)=DCH(I)+HIJ*CG(J)*CCELEC 
723:                 IF(QREF) THEN 
724:                   DCH(I)=DCH(I)+HIJ*CGREF(J)*CCELEC 
725:                 ENDIF 
726:               ENDIF 
727:             ENDDO !JJ 
728:           ENDIF !(CGT.NE.ZERO) 
729:         ENDIF  ! MOLTYPE 
730:         ENDIF  ! QPARTBIN 
731:         ENDDO 
732:       ELSE !(CGMODEL.EQ. (anything but 0 or 1) ) 
733:         WRITE(OUTU,'(A,I5)')"CGMODEL out of range: ",CGMODEL 
734:         CALL DIE 
735:       ENDIF !(CGMODEL.EQ. 
736:       ECSUM=(ECSUM+0.5*ECSUM1) 
737:       IF (CGMODEL.GE.1) THEN 
738: C don't double count the simple coulombic interactions 
739:         ETERM(ELEC)=ETERM(ELEC)+ECSUM 
740:         EPROP(CGPOT)=ETERM(ELEC) 
741:         IF(LEWALD.AND.LELEC.AND.QETERM(ELEC)) 
742:      &           EPROP(CGPOT)=EPROP(CGPOT)+ETERM(IMELEC) 
743:         CALL SETMSR('ESEL',ECSUM) 
744:       ENDIF 
745:  
746:  
747: ##IF PARALLEL 
748:       ELSE    ! (IF MYNOD.EQ.0) 
749:       EPROP(CGPOT)=0.0 
750:       ENDIF   ! (IF MYNOD.EQ.0) 
751:  
752: ##ENDIF  ! (IF PARALLEL) 
753:  
754: 69      format(22(f10.5,3x)) 
755:         IFLAG = 0 
756:  
757:  
758:       RETURN 
759:       END ! SUBROUTINE QBENER 
760:  
761:       SUBROUTINE QNOCOD(DX,DY,DZ,NATOM) 
762: C----------------------------------------------------------------------- 
763: C  This subroutine will set all atom first derivatives (or any three 
764: C  arrays of length NATOM) to zero. 
765: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
766: ##INCLUDE '~/charmm_fcm/number.fcm' 
767:       INTEGER I,NATOM 
768:       REAL*8 DX(*),DY(*),DZ(*) 
769: C 
770:       DO I=1,NATOM 
771:         DX(I)=ZERO 
772:         DY(I)=ZERO 
773:         DZ(I)=ZERO 
774:       ENDDO 
775:       RETURN 
776:       END !SUB QNOCOD 
777:  
778:       SUBROUTINE DDCG(BNBND,NATOM,CG,X,Y,Z) 
779: C----------------------------------------------------------------------- 
780: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
781: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
782:       INTEGER BNBND(*) 
783:       INTEGER NATOM 
784:       REAL*8  CG(*) 
785:       REAL*8  X(*),Y(*),Z(*) 
786: C 
787: ##INCLUDE '~/charmm_fcm/heap.fcm' 
788: ##INCLUDE '~/charmm_fcm/inbnd.fcm' 
789: ##INCLUDE '~/charmm_fcm/stack.fcm' 
790: ##INCLUDE '~/charmm_fcm/stream.fcm' 
791: C 
792: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
793: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
794:       REAL*8 HIJ 
795: C 
796:       IF(NATOM.LE.0) RETURN 
797: C     Make sure we have the current non-bond flags and cut-offs. This 
798: C     will handle any changes made due to pert or tsm. 
799:       CALL GETBND(BNBND,.TRUE.) 
800: C 
801:       CALL NBDDCG(NATOM,HEAP(BNBND(JNB)),HEAP(BNBND(INBLO)), 
802:      &     HEAP(BNBND(IBLO14)),HEAP(BNBND(INB14)),CG,X,Y,Z,CTOFNB) 
803:       RETURN 
804:       END 
805:  
806:       SUBROUTINE NBDDCG(NATOM,JNB,INBLO,IBLO14,INB14,CG,X,Y,Z,CTOFNB) 
807: C----------------------------------------------------------------------- 
808: C  need to check these derivatives.  I'm not sure they are correct. 
809: C 
810: C     NATOM  - number of atoms 
811: C     JNB    - nonbond pair list  (INBLO(NATOM)) 
812: C     INBLO  - pointers into JNB  (NATOM) 
813: C     CG     - charges  (NATOM) 
814: C 
815: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
816: ##INCLUDE '~/charmm_fcm/stream.fcm' 
817: C 
818: ##INCLUDE '~/charmm_fcm/heap.fcm' 
819: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
820: C 
821: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
822: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
823:       REAL*8 HIJ 
824: C 
825:       INTEGER NATOM 
826:       INTEGER JNB(*) 
827:       INTEGER INBLO(*) 
828:       INTEGER IBLO14(*) 
829:       INTEGER INB14(*) 
830:       REAL*8 CG(*) 
831:       REAL*8 X(*),Y(*),Z(*) 
832:       REAL*8 CTOFNB 
833: C 
834:       REAL*8 C2OFNB 
835:       REAL*8 CRXI,CRYI,CRZI,DXI,DYI,DZI 
836:       REAL*8 S 
837:       INTEGER I,J,I1,J1,J2,NB,NPR,JPR,ITEMP 
838:       INTEGER NBINDX,JJ,JCTR 
839: C 
840:       DO I=1,NATOM 
841:         DDCH(I)=0.0 
842:       ENDDO 
843:       NB=0 
844:       C2OFNB=CTOFNB*CTOFNB 
845: C 
846:       ITEMP=0 
847:       DO 40 I=1,NATOM 
848:         NPR=INBLO(I)-ITEMP 
849:         ITEMP=INBLO(I) 
850:         IF (NPR.EQ.0) GOTO 30 
851: C 
852:         CRXI=X(I) 
853:         CRYI=Y(I) 
854:         CRZI=Z(I) 
855: C 
856:         DO 20 JPR=1,NPR 
857:           NB=NB+1 
858:           IF (JNB(NB).LT.0) THEN 
859:             J=-JNB(NB) 
860:           ELSE 
861:             J=JNB(NB) 
862:           ENDIF 
863: C 
864:           DXI=CRXI-X(J) 
865:           DYI=CRYI-Y(J) 
866:           DZI=CRZI-Z(J) 
867:           S=DXI*DXI+DYI*DYI+DZI*DZI 
868: C 
869:           IF (S.LT.C2OFNB) THEN 
870:             HIJ=1.0/SQRT(S) 
871:             DDCH(I)=DDCH(I)+HIJ*DCH(J) 
872:             DDCH(J)=DDCH(J)+HIJ*DCH(I) 
873:           ENDIF 
874: C 
875:    20   CONTINUE 
876:    30   CONTINUE 
877:    40 CONTINUE 
878: C 
879: C Calculate the contributions to DDCH from nonbond exclusion list 
880:       NBINDX=1 
881:       DO I=1,NATOM-1 
882:         IF (I.EQ.1) THEN 
883:           NBINDX=1 
884:         ELSE 
885:           NBINDX=IBLO14(I-1)+1 
886:         ENDIF 
887:         JCTR=IBLO14(I) 
888:         DO JJ=NBINDX,JCTR 
889:            J=IABS(INB14(JJ))  !!! SAP 
890: C          J=INB14(JJ) 
891:           S=(X(I)-X(J))*(X(I)-X(J))+(Y(I)-Y(J))*(Y(I)-Y(J)) 
892:      $     +(Z(I)-Z(J))*(Z(I)-Z(J)) 
893:           HIJ=EHA(I)+EHA(J) 
894:           HIJ=HIJ/SQRT((1.0+HIJ*HIJ*S)) 
895:           DDCH(I)=DDCH(I)+HIJ*DCH(J) 
896: c          DDCH(J)=DDCH(J)+HIJ*DCH(I) 
897:         ENDDO 
898:       ENDDO 
899:       DO I=1,NATOM 
900: C        DDCH(I)=-2.0*DDCH(I) 
901:         DDCH(I)=2.0*DDCH(I) 
902:       ENDDO 
903:  
904: C      WRITE(OUTU,*)(DCH(I),DDCH(I),I=1,NATOM) 
905:  
906:       RETURN 
907:       END 
908:  
909:        SUBROUTINE CGINV(NATOM,X,Y,Z,DX,DY,DZ,CG,ECH,EHA,QMOL, 
910:      $      MOLT,A,DA,B,QCGINVF,JNB,INBLO,INB14,IBLO14,IAC,QPARTBIN, 
911:      $      MOLTYPE,QPRNETA,CCNBA,CCNBB,CCNBC,CCNBD,IACNB,NITCC2,LOWTP, 
912:      $      OUTU) 
913: C---------------------------------------------------------------------- 
914: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
915: ##INCLUDE '~/charmm_fcm/consta.fcm' 
916: ##INCLUDE '~/charmm_fcm/flucq.fcm' !##FLUCQ 
917:  
918:       INTEGER NB,ITEMP,NATOM,JNB(*),INB14(*),INBLO(*),IBLO14(*) 
919:       REAL*8 X(*),Y(*),Z(*),CG(*),ECH(*),EHA(*),A(NATOM,*),B(*) 
920:       REAL*8 U(6,6),V(6) 
921:       REAL*8 DX(*),DY(*),DZ(*),DA(NATOM,*) 
922:       INTEGER N,MOLT(*),MOLTYPE(*) 
923:       INTEGER OUTU 
924: ##IF BLOCK 
925: C will need to change these declarations to get block working with 
926: C charge inversion 
927:       INTEGER IBLOCK(1) 
928:       REAL*8  BLCOE(1),BLCOV(1) 
929: ##ENDIF ! BLOCK 
930:       REAL*8  CCNBA(*),CCNBB(*),CCNBC(*),CCNBD(*) 
931:       INTEGER IACNB(*),NITCC2,LOWTP(*) 
932:       LOGICAL QMOL,QCGINVF,QPRNETA 
933:       REAL*8 RIJ2,RIJ,HIJ,RIJ1,HIJ1 
934:       INTEGER QPARTBIN(*),IAC(*),NPR,JPR,JCTR,JJ,NBINDX 
935:       REAL*8 EVDW,ECG 
936:       LOGICAL LUSED 
937: C 
938: C       DIMENSION INDX(6) 
939:        INTEGER I,J,K 
940:        INTEGER NMOL 
941:  
942:        DO I=1,NATOM 
943:          DO J=1,NATOM 
944:            A(I,J)=0.0 
945:          ENDDO 
946:          IF (MOLTYPE(I).EQ.0) THEN 
947:            A(I,I)=2.0*EHA(I) 
948:          ELSE IF (MOLTYPE(I).GT.0) THEN  ! for SPC water molecules 
949:            IF (MOLTYPE(I).EQ.1) THEN 
950:              A(I,I)=367.0/CCELEC 
951:            ELSE 
952:              A(I,I)=392.2/CCELEC 
953:            ENDIF 
954:          ENDIF 
955:        ENDDO 
956: c 
957:  
958:             CALL ENBFS8(EVDW,ECG,.TRUE.,.FALSE.,1,NATOM,CG,JNB,INBLO, 
959:      &           CCNBA, CCNBB, CCNBC, CCNBD, IACNB, NITCC2,LOWTP, 
960: ##IF BLOCK 
961:      $           IBLOCK,BLCOE,BLCOV, 
962: ##ENDIF 
963:      &           A,.TRUE., 
964:      $           QFLUC,FQCFOR,  !##FLUCQ 
965:      $           LUSED) 
966:  
967:  
968:  
969:       NBINDX=1 
970:       DO I=1,NATOM-1 
971:         IF (MOLTYPE(I).EQ.0) THEN 
972:           IF (I.EQ.1) THEN 
973:             NBINDX=1 
974:           ELSE 
975:             NBINDX=IBLO14(I-1)+1 
976:           ENDIF 
977:           JCTR=IBLO14(I) 
978:           DO JJ=NBINDX,JCTR 
979:              J=IABS(INB14(JJ)) 
980:             IF (J.GT.0) THEN 
981: C              write(*,*) "CGINV: i,j",i,j 
982:               RIJ2=(X(I)-X(J))*(X(I)-X(J)) 
983:      $               +(Y(I)-Y(J))*(Y(I)-Y(J)) 
984:      $               +(Z(I)-Z(J))*(Z(I)-Z(J)) 
985:               HIJ=EHA(I)+EHA(J) 
986:               HIJ=HIJ/SQRT(1.0+HIJ*HIJ*RIJ2) 
987:               A(I,J)=HIJ 
988:               A(J,I)=HIJ 
989:  
990:             ELSE 
991:               HIJ = 1.0/SQRT(RIJ2) 
992:               A(I,J) = HIJ 
993:               A(J,I) = HIJ 
994:              ENDIF 
995:  
996:           ENDDO 
997:         ELSE IF (MOLTYPE(I).GT.0) THEN  ! for SPC water molecules 
998:           IF (I.EQ.1) THEN 
999:             NBINDX=1 
1000:           ELSE 
1001:             NBINDX=IBLO14(I-1)+1 
1002:           ENDIF 
1003:           JCTR=IBLO14(I) 
1004:           DO JJ=NBINDX,JCTR 
1005:              J=IABS(INB14(JJ)) 
1006:             IF (J.GT.0) THEN 
1007:               IF (MOLTYPE(I).EQ.1.OR.MOLTYPE(J).EQ.1) THEN 
1008:                 A(I,J)=276.0/CCELEC 
1009:                 A(J,I)=276.0/CCELEC 
1010:               ELSE 
1011:                 A(I,J)=196.0/CCELEC 
1012:                 A(J,I)=196.0/CCELEC 
1013:               ENDIF 
1014:             ENDIF 
1015:           ENDDO !JJ 
1016:           NBINDX=NBINDX+IBLO14(I) 
1017:         ELSE 
1018:           write(6,*) 'CGINV: other waters not supported yet' 
1019:         ENDIF 
1020:       ENDDO 
1021:  
1022:       IF(QPRNETA) THEN 
1023:         DO I=1,NATOM 
1024:           write(6,*) 'Row ',I 
1025:           WRITE(6,60) (A(I,J), J=1,NATOM) 
1026:         ENDDO 
1027:       ENDIF 
1028:  60   FORMAT(12(F10.5,1x)) 
1029:  
1030:        IF (QMOL) THEN 
1031: C assumes partitions are consecutive!!! 
1032: c         write(*,*) "CGINV: charge normalization over partitions" 
1033:          NMOL=0 
1034:          DO I=1,NATOM 
1035:            IF (QPARTBIN(I).NE.QPARTBIN(I-1)) NMOL=NMOL+1 
1036:          ENDDO 
1037:          K=1 
1038:          B(1)=CG(1) 
1039:          DO I=2,NATOM 
1040:            IF (QPARTBIN(I).EQ.QPARTBIN(I-1)) THEN 
1041:              DO J=1,NATOM 
1042:                A(I,J)=A(K,J)-A(I,J) 
1043:              ENDDO 
1044:              B(I)=(ECH(I)-ECH(K))/CCELEC 
1045:              B(K)=B(K)+CG(I) 
1046:            ELSE 
1047:              K=I 
1048:              B(I)=CG(I) 
1049:            ENDIF 
1050:          ENDDO 
1051:          DO J=1,NATOM 
1052:            IF(QPARTBIN(J).EQ.QPARTBIN(1)) THEN 
1053:              A(1,J)=1.0 
1054:            ELSE 
1055:              A(1,J)=0.0 
1056:            ENDIF 
1057:          ENDDO 
1058:          DO I=2,NATOM 
1059:            IF (QPARTBIN(I).NE.QPARTBIN(I-1)) THEN 
1060:              DO J=1,NATOM 
1061:                IF(QPARTBIN(J).EQ.QPARTBIN(I)) THEN 
1062:                  A(I,J)=1.0 
1063:                ELSE 
1064:                  A(I,J)=0.0 
1065:                ENDIF 
1066:              ENDDO 
1067:            ENDIF 
1068:          ENDDO 
1069:        ELSE 
1070: c         write(*,*) "CGINV: charge normalization over all atoms" 
1071:          DO I=2,NATOM 
1072:            DO J=1,NATOM 
1073:               A(I,J)=A(1,J)-A(I,J) 
1074:            ENDDO 
1075:          ENDDO 
1076:          DO I=1,NATOM 
1077:            B(I)=(ECH(I)-ECH(1))/CCELEC 
1078:            A(1,I)=1.0 
1079:          ENDDO 
1080:        ENDIF 
1081:       IF(QPRNETA) THEN 
1082:         write(6,*) 'A matrix going into GAUSSJ' 
1083:         DO I=1,NATOM 
1084:           write(6,*) 'Row ',I 
1085:           WRITE(6,60) (A(I,J), J=1,NATOM) 
1086:         ENDDO 
1087:       ENDIF 
1088:       QPRNETA=.FALSE. 
1089:        CALL GAUSSJ(A,NATOM,NATOM,B,1,1) 
1090: C So A is the inverse matrix now 
1091:  
1092:        RETURN 
1093:        END 
1094:  
1095:       SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) 
1096: C---------------------------------------------------------------------- 
1097: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1098: ##INCLUDE '~/charmm_fcm/number.fcm' 
1099: C 
1100:       INTEGER N,NP,M,MP 
1101:       REAL*8 A(NP,NP),B(NP,MP) 
1102: C 
1103:       INTEGER NMAX 
1104:       PARAMETER (NMAX=1000) 
1105:       INTEGER IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) 
1106: C 
1107:       INTEGER I,J,K,L,LL,IROW,ICOL 
1108:       REAL*8  BIG,DUM,PIVINV 
1109: C 
1110:       IF(N.GT.NMAX) CALL WRNDIE(-4,'<GAUSSJ>','N is too large.') 
1111: C 
1112:       DO J=1,N 
1113:         IPIV(J)=0 
1114:       ENDDO 
1115:       DO I=1,N 
1116:         BIG=ZERO 
1117:         DO J=1,N 
1118:           IF(IPIV(J).NE.1)THEN 
1119:             DO K=1,N 
1120:               IF (IPIV(K).EQ.0) THEN 
1121:                 IF (ABS(A(J,K)).GE.BIG)THEN 
1122:                   BIG=ABS(A(J,K)) 
1123:                   IROW=J 
1124:                   ICOL=K 
1125:                 ENDIF 
1126:               ELSE IF (IPIV(K).GT.1) THEN 
1127:                  CALL WRNDIE(-4,'<GAUSSJ>','Singular matrix.1') 
1128:               ENDIF 
1129:             ENDDO 
1130:           ENDIF 
1131:         ENDDO 
1132:         IPIV(ICOL)=IPIV(ICOL)+1 
1133:         IF (IROW.NE.ICOL) THEN 
1134:           DO L=1,N 
1135:             DUM=A(IROW,L) 
1136:             A(IROW,L)=A(ICOL,L) 
1137:             A(ICOL,L)=DUM 
1138:           ENDDO 
1139:           DO L=1,M 
1140:             DUM=B(IROW,L) 
1141:             B(IROW,L)=B(ICOL,L) 
1142:             B(ICOL,L)=DUM 
1143:           ENDDO 
1144:         ENDIF 
1145:         INDXR(I)=IROW 
1146:         INDXC(I)=ICOL 
1147:         IF (A(ICOL,ICOL).EQ.0.)  
1148:      $       CALL WRNDIE(-4,'<GAUSSJ>','Singular matrix.2') 
1149:         PIVINV=1./A(ICOL,ICOL) 
1150:         A(ICOL,ICOL)=1. 
1151:         DO L=1,N 
1152:           A(ICOL,L)=A(ICOL,L)*PIVINV 
1153:         ENDDO 
1154:         DO L=1,M 
1155:           B(ICOL,L)=B(ICOL,L)*PIVINV 
1156:         ENDDO 
1157:         DO LL=1,N 
1158:           IF(LL.NE.ICOL)THEN 
1159:             DUM=A(LL,ICOL) 
1160:             A(LL,ICOL)=0. 
1161:             DO L=1,N 
1162:               A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 
1163:             ENDDO 
1164:             DO L=1,M 
1165:               B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 
1166:             ENDDO 
1167:           ENDIF 
1168:         ENDDO 
1169:       ENDDO 
1170:       DO L=N,1,-1 
1171:         IF(INDXR(L).NE.INDXC(L))THEN 
1172:           DO K=1,N 
1173:             DUM=A(K,INDXR(L)) 
1174:             A(K,INDXR(L))=A(K,INDXC(L)) 
1175:             A(K,INDXC(L))=DUM 
1176:           ENDDO 
1177:         ENDIF 
1178:       ENDDO 
1179:       RETURN 
1180:       END 
1181:  
1182:       SUBROUTINE CGCOPY(CG,CGX,NATOM,NATIM) 
1183: C---------------------------------------------------------------------- 
1184: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1185:       REAL*8 CG(*),CGX(*) 
1186:       INTEGER NATOM,NATIM,I,NATX 
1187:       IF (NATIM.GT.NATOM) THEN 
1188:          NATX=NATIM 
1189:       ELSE 
1190:          NATX=NATOM 
1191:       ENDIF 
1192:       DO I=1,NATX 
1193:          CGX(I)=CG(I) 
1194:       ENDDO 
1195:       RETURN 
1196:       END 
1197:  
1198:       SUBROUTINE CGCPY(CGX) 
1199: C---------------------------------------------------------------------- 
1200: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1201: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1202: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1203:       REAL*8 CGX(*) 
1204:       INTEGER I 
1205:       DO I=1,NATOM 
1206:          CG(I)=CGX(I) 
1207:       ENDDO 
1208:       RETURN 
1209:       END 
1210:  
1211:       SUBROUTINE QCGTRANSI(DCH,NATOM,NATIM,IMATTR) 
1212: C---------------------------------------------------------------------- 
1213: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1214: ##INCLUDE '~/charmm_fcm/number.fcm' 
1215:       REAL*8 DCH(*) 
1216:       INTEGER NATOM,NATIM,IMATTR(*) 
1217:       INTEGER I 
1218:  
1219:        DO I=NATOM+1,NATIM 
1220:          DCH(IMATTR(I)) = DCH(IMATTR(I)) + DCH(I) 
1221:          DCH(I)=ZERO 
1222:        ENDDO 
1223:  
1224:        RETURN 
1225:        END 
1226:  
1227:       SUBROUTINE CPCGIMG(IMATTRX) 
1228: C---------------------------------------------------------------------- 
1229: C THIS ROUTINE COPIES CG FROM MAIN SET TO IMAGE ATOMS 
1230: C FOR USE IN CHARGE DYNAMICS 
1231: C 
1232: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1233: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1234: ##INCLUDE '~/charmm_fcm/image.fcm' 
1235: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1236: C 
1237:       INTEGER IMATTRX(*) 
1238: C 
1239:       INTEGER I,J 
1240:       IF (NATIM.GT.NATOM) THEN 
1241:         DO I=NATOM+1,NATIM 
1242:            J=IMATTRX(I) 
1243:            CG(I)=CG(J) 
1244:         ENDDO 
1245:       ENDIF 
1246:       RETURN 
1247:       END 
1248:  
1249:       SUBROUTINE PRCGIMG(IMATTRX) 
1250: C---------------------------------------------------------------------- 
1251: C THIS ROUTINE COPIES CG FROM MAIN SET TO IMAGE ATOMS 
1252: C FOR USE IS CHARGE DYNAMICS 
1253: C 
1254: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1255: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1256: ##INCLUDE '~/charmm_fcm/image.fcm' 
1257: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1258: C 
1259:       INTEGER IMATTRX(*) 
1260: C 
1261:       INTEGER I,J 
1262:       IF (NATIM.GT.NATOM) THEN 
1263:         WRITE(6,*)"IMAGE CHARGES>>>> comparison" 
1264:         DO I=NATOM+1,NATIM 
1265:            J=IMATTRX(I) 
1266:            write(6,'(2I10,2F20.6)')I,J,CG(I),CG(J) 
1267:         ENDDO 
1268:       ENDIF 
1269:       RETURN 
1270:       END 
1271:  
1272:       SUBROUTINE FQPOLTEN(NATOM,X,Y,Z,DX,DY,DZ,CG,ECH,EHA,QMOL, 
1273:      $      MOLT,A,DA,B,QCGINVF,JNB,INBLO,INB14,IBLO14,IAC,QPARTBIN, 
1274:      $      MOLTYPE,QPRNETA,CCNBA,CCNBB,CCNBC,CCNBD,IACNB,NITCC2,LOWTP, 
1275:      $      OUTU,INTER,R,delR,temp) 
1276: c------------------------------------------------------------------------ 
1277: C    subroutine to calcuate molecualr polarizability tensor 
1278: c    from atomic hardness parameters 
1279: c   USE ONLY FOR SINGLE MOLECULES WITH NET ZERO CHARGE!!!!! 
1280: C    (modify normalization scheme for use with charged species 
1281: c    or more exotic normalization schemes 
1282: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1283: ##INCLUDE '~/charmm_fcm/consta.fcm' 
1284:       INTEGER NB,ITEMP,NATOM,JNB(*),INB14(*),INBLO(*),IBLO14(*) 
1285:       REAL*8 X(*),Y(*),Z(*),CG(*),ECH(*),EHA(*),A(NATOM,*),B(*) 
1286:       REAL*8 ISICK(3,3),BSICK(3,1) 
1287:       REAL*8 U(6,6),V(6) 
1288:       REAL*8 DX(*),DY(*),DZ(*),DA(NATOM,*) 
1289:       INTEGER N,MOLT(*),MOLTYPE(*) 
1290:       INTEGER OUTU,KK,LL,JJ,NATMM1,L,M,TIPFLAG 
1291: ##IF BLOCK 
1292: C will need to change these declarations to get block working with 
1293: C charge inversion 
1294:       INTEGER IBLOCK(1) 
1295:       REAL*8  BLCOE(1),BLCOV(1) 
1296: ##ENDIF ! BLOCK 
1297:       REAL*8  CCNBA(*),CCNBB(*),CCNBC(*),CCNBD(*) 
1298:       INTEGER IACNB(*),NITCC2,LOWTP(*) 
1299:       LOGICAL QMOL,QCGINVF,QPRNETA 
1300:       REAL*8 RIJ2,RIJ,HIJ,RIJ1,HIJ1 
1301:       INTEGER QPARTBIN(*),IAC(*),NPR,JPR,JCTR,NBINDX 
1302:       REAL*8 EVDW,ECG, alpha(3,3),inter(3,NATOM),R(3,NATOM) 
1303:       REAL*8 delR(NATOM,3), temp(3,NATOM) 
1304:       LOGICAL LUSED 
1305:       Character*2 comp 
1306: C 
1307:       INTEGER I,J,K 
1308:       INTEGER NMOL 
1309:  
1310:        TIPFLAG=0 
1311:  
1312: c       set up the eta matrix 
1313:  
1314:         DO I=1,NATOM 
1315:          DO J=1,NATOM 
1316:            A(I,J)=0.0 
1317:          ENDDO 
1318:        ENDDO 
1319:  
1320:         DO I=1,3 
1321:           DO J=1,3 
1322:          ISICK(I,J) = 0.0 
1323:          alpha(I,J) = 0.0 
1324:           ENDDO 
1325:         ENDDO 
1326:  
1327:         DO I = 1,3 
1328:          DO J = 1,NATOM 
1329:            inter(I,J) = 0.0 
1330:          ENDDO 
1331:         ENDDO 
1332:  
1333:        DO I=1,NATOM 
1334:          DO J=1,NATOM 
1335:            A(I,J)=0.0 
1336:          ENDDO 
1337:          IF (MOLTYPE(I).EQ.0) THEN 
1338:            A(I,I)=2.0*EHA(I) 
1339:          ELSE IF (MOLTYPE(I).GT.0) THEN  ! for SPC / TIP4P water molecules 
1340:            IF (MOLTYPE(I).EQ.1) THEN 
1341:              A(I,I)=367.0/CCELEC 
1342:            ELSEIF (MOLTYPE(I).EQ.2) THEN 
1343:              A(I,I)=392.2/CCELEC 
1344:             ELSEIF (MOLTYPE(I).EQ.3) THEN 
1345:              A(I,I) = 371.6/CCELEC   !  TIP4P OM SITE 
1346:             ELSEIF (MOLTYPE(I).EQ.4) THEN 
1347:              A(I,I) = 353.0/CCELEC   ! TIP4P HYDROGEN SITE 
1348:              ELSEIF (MOLTYPE(I).EQ.9) THEN 
1349:              A(I,I) = 0.0   !   TIP4P oxygen site (no charge) 
1350:            ELSE 
1351:            ENDIF 
1352:          ENDIF 
1353:        ENDDO 
1354:  
1355:       NBINDX=1 
1356:       DO I=1,NATOM-1 
1357:         IF (MOLTYPE(I).EQ.0) THEN 
1358:           IF (I.EQ.1) THEN 
1359:             NBINDX=1 
1360:           ELSE 
1361:             NBINDX=IBLO14(I-1)+1 
1362:           ENDIF 
1363:           JCTR=IBLO14(I) 
1364:           DO JJ=NBINDX,JCTR 
1365:  
1366:              J=IABS(INB14(JJ)) !!! SAP 
1367:             IF (J.GT.0) THEN 
1368:               RIJ2=(X(I)-X(J))*(X(I)-X(J)) 
1369:      $               +(Y(I)-Y(J))*(Y(I)-Y(J)) 
1370:      $               +(Z(I)-Z(J))*(Z(I)-Z(J)) 
1371:               HIJ=EHA(I)+EHA(J) 
1372:               HIJ=HIJ/SQRT(1.0+HIJ*HIJ*RIJ2) 
1373:               A(I,J)=HIJ 
1374:               A(J,I)=HIJ 
1375:             ELSE 
1376:               HIJ = 1.0/SQRT(RIJ2) 
1377:               A(I,J) = HIJ 
1378:               A(J,I) = HIJ 
1379:              ENDIF 
1380:  
1381:           ENDDO 
1382:  
1383:         ELSE IF ( (MOLTYPE(I).GT.0).AND.(MOLTYPE(I).LE.2) ) THEN 
1384: C     for SPC/TIP4P water molecules 
1385:           IF (I.EQ.1) THEN 
1386:             NBINDX=1 
1387:           ELSE 
1388:             NBINDX=IBLO14(I-1)+1 
1389:           ENDIF 
1390:           JCTR=IBLO14(I) 
1391:           DO JJ=NBINDX,JCTR 
1392:              J=IABS(INB14(JJ)) 
1393:             IF (J.GT.0) THEN 
1394:               IF (MOLTYPE(I).EQ.1.OR.MOLTYPE(J).EQ.1) THEN 
1395:                 A(I,J)=276.0/CCELEC 
1396:                 A(J,I)=276.0/CCELEC 
1397:               ELSE 
1398:                 A(I,J)=196.0/CCELEC 
1399:                 A(J,I)=196.0/CCELEC 
1400:               ENDIF 
1401:             ENDIF 
1402:           ENDDO !JJ 
1403:           NBINDX=NBINDX+IBLO14(I) 
1404:  
1405:         ELSE IF ( (MOLTYPE(I).GT.2).AND.(MOLTYPE(I).Lt.9)) THEN 
1406: C     for TIP4P water molecu 
1407:         TIPFLAG=1 
1408:           IF (I.EQ.1) THEN 
1409:             NBINDX=1 
1410:           ELSE 
1411:             NBINDX=IBLO14(I-1)+1 
1412:           ENDIF 
1413:           JCTR=IBLO14(I) 
1414: C            CRXI=X(I) 
1415: C            CRYI=Y(I) 
1416: C            CRZI=Z(I) 
1417:             DO JJ=NBINDX,JCTR 
1418:                J=IABS(INB14(JJ))   !!! SAP 
1419: C              J=INB14(JJ) 
1420:               IF (J.GT.0) THEN 
1421:  
1422:                 IF (MOLTYPE(J).NE.9)  THEN 
1423:                 IF (MOLTYPE(I).EQ.3.OR.MOLTYPE(J).EQ.3) THEN 
1424:                    A(I,J) = 286.4/CCELEC 
1425:                    A(J,I) = 286.4/CCELEC 
1426:                 ELSE 
1427:                     A(I,J) = 203.6/CCELEC 
1428:                     A(J,I) = 203.6/CCELEC 
1429:                 ENDIF 
1430:               ENDIF 
1431:             ENDIF 
1432:             ENDDO !JJ 
1433: C            ENDIF 
1434:  
1435:         ELSE 
1436:           write(6,*) 'CGINV: other waters not supported yet' 
1437:         ENDIF 
1438:       ENDDO 
1439:  
1440:           IF (TIPFLAG.EQ.1) THEN 
1441:           DO L=1,3 
1442:             DO M=1,3 
1443:               ISICK(L,M) = A(L+1,M+1) 
1444:             ENDDO 
1445:          ENDDO 
1446:          ELSE 
1447:          ENDIF 
1448:  
1449: C     compute inverse of eta matrix 
1450:  
1451:       IF (TIPFLAG.EQ.0) THEN 
1452:        CALL GAUSSJ(A,NATOM,NATOM,B,1,1) 
1453:        ELSEIF (TIPFLAG.EQ.1) THEN 
1454:          CALL GAUSSJ(ISICK,3,3,BSICK,1,1) 
1455:       ELSE 
1456:       ENDIF 
1457:  
1458: c     A is the inverse matrix now 
1459:  
1460: c     now set up the coordinates matrix   `R` 
1461:  
1462:             IF (TIPFLAG.EQ.0) THEN 
1463:             DO JJ=1,NATOM 
1464:              R(1,JJ)= X(JJ) 
1465:              R(2,JJ)= Y(JJ) 
1466:              R(3,JJ)= Z(JJ) 
1467:             ENDDO 
1468:  
1469:             ELSE 
1470:  
1471:             DO JJ=1,3 
1472:               R(1,JJ) = X(JJ+1) 
1473:               R(2,JJ) = Y(JJ+1) 
1474:               R(3,JJ) = Z(JJ+1) 
1475:             ENDDO 
1476:  
1477:             ENDIF 
1478:  
1479: c     now set up the delR matrix 
1480:  
1481:           DO JJ=1,3 
1482:               DO KK=1,NATOM 
1483:                temp(JJ,KK)=R(JJ,KK)-R(JJ,1) 
1484:                delR(KK,JJ)=temp(JJ,KK) 
1485:               ENDDO 
1486:           ENDDO 
1487:  
1488:         IF (TIPFLAG.EQ.0) THEN 
1489:  
1490: c     multiply  R*inv(eta) 
1491:  
1492:           DO JJ=1,3 
1493:             DO KK=1,NATOM 
1494:                DO LL=1,NATOM 
1495:               inter(JJ,KK)  = inter(JJ,KK)+R(JJ,LL)*A(LL,KK) 
1496:               ENDDO 
1497:              ENDDO 
1498:            ENDDO 
1499:  
1500:           DO JJ=1,3 
1501:             DO KK=1,3 
1502:                DO LL=1,NATOM 
1503:                alpha(JJ,KK) = alpha(JJ,KK)+inter(JJ,LL)*R(KK,LL) 
1504:               ENDDO 
1505:              ENDDO 
1506:            ENDDO 
1507:  
1508:        ELSEIF (TIPFLAG.EQ.1) THEN 
1509:  
1510: c    multiply  R*inv(eta) 
1511:  
1512:           DO JJ=1,3 
1513:             DO KK=1,3 ! NATOM 
1514:                DO LL=1,3 ! NATOM 
1515:               inter(JJ,KK) = inter(JJ,KK) + R(JJ,LL)*ISICK(LL,KK) 
1516:               ENDDO 
1517:              ENDDO 
1518:            ENDDO 
1519:  
1520:           DO JJ=1,3 
1521:             DO KK=1,3 
1522:                DO LL=1,3 ! NATOM 
1523:               alpha(JJ,KK)  = alpha(JJ,KK)+inter(JJ,LL)*R(KK,LL) 
1524:               ENDDO 
1525:              ENDDO 
1526:            ENDDO 
1527:  
1528:  
1529:         ELSE 
1530:            ENDIF 
1531:  
1532:          write(OUTU,*) " " 
1533:          write(OUTU,*)" MOLECULAR POLARIZABILITY TENSOR ", TIPFLAG 
1534:          write(OUTU,*) "          X            Y            Z " 
1535:          write(OUTU,*) "  _____________________________________" 
1536:          DO I = 1,3 
1537:           if (I.eq.1) comp='X|' 
1538:           if (I.eq.2) comp='Y|' 
1539:           if (I.eq.3) comp='Z|' 
1540:           write(6,22) comp,(alpha(I,J),J=1,3) 
1541:           ENDDO 
1542:         write(OUTU,*) " " 
1543:  
1544: 22    format    (1x,a2,2x,3(f10.5,2x)) 
1545: 23    format    (4(f10.5,2x)) 
1546:           return 
1547:            end 
1548:  
1549: C  routine to get forces on charges arising from the GB model for solvation 
1550:  
1551:       SUBROUTINE  GBFONCHEQ(alpha,Eps_GB,X,Y,Z) 
1552: C----------------------------------------------------------------------- 
1553: C     This subroutine calculates the forces on charges arising from the 
1554: C     generalized born expression for the solvation energy.  Coming into 
1555: C     the routine, the Born radii, alpha(i), should be calculated, of if 
1556: C     needed, updated in a multiple time step approach, and available. 
1557: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1558: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1559: ##INCLUDE '~/charmm_fcm/consta.fcm' 
1560: ##INCLUDE '~/charmm_fcm/number.fcm' 
1561: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
1562: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1563: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1564: ##INCLUDE '~/charmm_fcm/inbnd.fcm' 
1565: ##INCLUDE '~/charmm_fcm/parallel.fcm' 
1566:  
1567:        Integer I, J 
1568:        Real*8  alpha(*),XDIFF,YDIFF,ZDIFF,RIJ2,term1,term2,K 
1569:        Real*8  X(*), Y(*), Z(*), Eps_GB 
1570:  
1571: c     make sure units are consistent 
1572:  
1573: c     prefactor with solute and solvent dielectrics 
1574:       K = -CCELEC* ( (1.0/EPS) - (1.0/Eps_GB) ) 
1575:  
1576: C     looping structure not the most efficient for now 
1577:       DO I = 1, NATOM 
1578:          DO J = 1, NATOM 
1579:             IF (I.NE.J) THEN 
1580:                XDIFF = X(I) - X(J) 
1581:                YDIFF = Y(I) - Y(J) 
1582:                ZDIFF = Z(I) - Z(J) 
1583:                RIJ2 = XDIFF**2 + YDIFF**2 + ZDIFF**2 
1584:                term1 = alpha(I)*alpha(J) 
1585: c     THE 8.0  is for P6 = 8.0 in GBMV methods 
1586:                term2 = RIJ2/(8.0*term1)   
1587:                DCH(I) = DCH(I) + ( -1.0*0.5 * K * CG(J)/ 
1588:      $              (sqrt( RIJ2 + term1*DEXP(-term2) ) ) ) 
1589:             ELSE 
1590:                DCH(I) = DCH(I) + ( -1.0*0.5*K*( CG(I)/alpha(I) )  ) 
1591:             ENDIF 
1592:          ENDDO !  J LOOP 
1593:       ENDDO    !  I LOOP 
1594:       RETURN 
1595:       END 
1596:  
1597:       SUBROUTINE CHECKETA(QCHEQRDPRM) 
1598: C----------------------------------------------------------------------- 
1599: C     Routine to check to see the state of hardness parameters 
1600: C     for CHEQ code. 
1601: C     IF indeed one wishes to use CHEQ, then they have to have 
1602: C     read in the correct CHEQ model paramater.s 
1603: C     IF the parameters haven't been read, but the user tries to 
1604: C     run any CHEQ method, there should be an alert and die. 
1605: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1606: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1607: ##INCLUDE '~/charmm_fcm/number.fcm' 
1608: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
1609: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1610: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1611: ##INCLUDE '~/charmm_fcm/parallel.fcm' 
1612:  
1613:       INTEGER I 
1614:       LOGICAL FOUND,QCHEQRDPRM 
1615:  
1616:       QCHEQRDPRM = .TRUE. 
1617:  
1618: ##IF PARALLEL 
1619:       IF (MYNOD.EQ.0) THEN 
1620: ##ENDIF 
1621:  
1622:       FOUND = .FALSE. 
1623:       DO I = 1, NATOM 
1624:          IF (MOLTYPE(I).NE.9) THEN 
1625:             IF (EHA(I).LE.0.0) FOUND = .TRUE. 
1626:          ENDIF 
1627:       ENDDO 
1628:  
1629:       IF (FOUND) THEN 
1630:          QCHEQRDPRM = .FALSE. 
1631:       ELSE 
1632:          QCHEQRDPRM = .TRUE. 
1633:       ENDIF 
1634:  
1635: ##IF PARALLEL 
1636:       ENDIF 
1637: ##ENDIF 
1638:  
1639:       RETURN 
1640:       END 
1641:  
1642:       SUBROUTINE CHECKQNORM(QCHEQNORM) 
1643: C----------------------------------------------------------------------- 
1644: C     Routine to check that there are no CHEQ partitions defined 
1645: c     with zero atoms; this indicates some error and must be corrected 
1646: c     by the user before any valid CHEQ method is used. 
1647: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
1648: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
1649: ##INCLUDE '~/charmm_fcm/number.fcm' 
1650: ##INCLUDE '~/charmm_fcm/derivq.fcm' 
1651: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
1652: ##INCLUDE '~/charmm_fcm/psf.fcm' 
1653: ##INCLUDE '~/charmm_fcm/stream.fcm' 
1654: ##INCLUDE '~/charmm_fcm/parallel.fcm' 
1655:  
1656:       LOGICAL FOUND,QCHEQNORM 
1657:  
1658:       QCHEQNORM = .TRUE. 
1659:  
1660: ##IF PARALLEL 
1661:       IF (MYNOD.EQ.0) THEN 
1662: ##ENDIF 
1663:  
1664:       FOUND = .FALSE. 
1665:       IF (QNPART.LE.0) FOUND = .TRUE. 
1666:  
1667:       IF (FOUND) THEN 
1668:           QCHEQNORM = .FALSE. 
1669:       ELSE 
1670:           QCHEQNORM = .TRUE. 
1671:       ENDIF 
1672:  
1673:       FOUND = .FALSE. 
1674:  
1675: ##IF PARALLEL 
1676:       ENDIF 
1677: ##ENDIF 
1678: C 
1679: ##ELSE (cheq_main) 
1680:       SUBROUTINE NULL_QCGSET 
1681: ##ENDIF (cheq_main) 
1682:       RETURN 
1683:       END 


r21974/fqminut.src 2017-01-21 10:33:58.390250000 +0000 r21973/fqminut.src 2017-01-21 10:34:02.330250000 +0000
  1: C  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/CHARMM35/source/cheq/fqminut.src' in revision 21973
  2: ##IF CHEQ (fqminut_main) 
  3: CCC The unused dummy variables are here in the following two routines 
  4: CCC for later development and consistency with their original routines. 
  5:       SUBROUTINE GETVR1Q(VARB,CG) 
  6: C     Put the coordinates and crystal variables into the variable array. 
  7: C 
  8: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
  9: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
 10: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
 11: C 
 12:       REAL*8      VARB(*), CG(*) 
 13: C 
 14:       INTEGER     I, IP, J 
 15: C 
 16:       IP = 0 
 17:  
 18:       DO I=1,QNPART 
 19:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
 20:           IP=IP+1 
 21:           VARB(IP)=CG(QPART(J)) 
 22:         ENDDO 
 23:       ENDDO 
 24: C 
 25:       RETURN 
 26:       END 
 27:        
 28: C----------------------------------------------------------------------- 
 29:       SUBROUTINE PUTVR1Q(VARB,CG,CG_FIX,CGT) 
 30: C     Fill the coordinate and crystal arrays with the variables. 
 31: C 
 32: ##INCLUDE '~/charmm_fcm/impnon.fcm' 
 33: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
 34: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
 35: C 
 36:       REAL*8      VARB(*),CG(*) 
 37:       REAL*8      CG_FIX(*),CGT(*) 
 38: C 
 39: C 
 40:       INTEGER     I, IP, J 
 41:       REAL*8  QDIF 
 42: C 
 43:       IP = 0 
 44: C 
 45:       DO I=1,QNPART 
 46:         CGT(I)=0.0 
 47:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
 48:           IP=IP+1 
 49:           CGT(I)=CGT(I)+VARB(IP) 
 50:           CG(QPART(J))=VARB(IP) 
 51:         ENDDO 
 52: !        QDIF=(CG_FIX(I)-CGT(I))/((QPRTPTR(I)+1)-QPRTPTR(I+1)) 
 53: !        DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
 54: !           CG(QPART(J))=CG(QPART(J))-QDIF 
 55: !        ENDDO 
 56:       ENDDO 
 57: C 
 58:       RETURN 
 59:       END 
 60:  
 61: C----------------------------------------------------------------------- 
 62:       SUBROUTINE MOLCGT(NATOM,CG,CG_FIX) 
 63: c  This routine computes the total charge on each  
 64: c  partition over which charge is nomalized. 
 65:       REAL*8 CG(*),CG_FIX(*) 
 66:       INTEGER NATOM 
 67: ##INCLUDE '~/charmm_fcm/dimens.fcm' 
 68: ##INCLUDE '~/charmm_fcm/cheqdyn.fcm' 
 69:       INTEGER I, I1,J,J1 
 70:       DO I=1,QNPART 
 71:         CG_FIX(I)=0.0 
 72:         DO J=QPRTPTR(I)+1,QPRTPTR(I+1) 
 73:           CG_FIX(I)=CG_FIX(I)+CG(QPART(J)) 
 74:         ENDDO 
 75: C      write(6,'(A,I5,F15.4)')"CGFIX",I,CG_FIX(I) 
 76:       ENDDO 
 77:  
 78:       RETURN 
 79: C 
 80: ##ELSE (fqminut_main) 
 81:       subroutine fqminut_dummy() 
 82: ##ENDIF (fqminut_main) 
 83:       END 


r21974/install.com 2017-01-21 10:33:57.258250000 +0000 r21973/install.com 2017-01-21 10:34:01.002250000 +0000
141: set xlf95 = 0141: set xlf95 = 0
142: set osx = 0142: set osx = 0
143: set ifc = 0143: set ifc = 0
144: set efc = 0144: set efc = 0
145: set amd64 = 0145: set amd64 = 0
146: set x86_64 = 0146: set x86_64 = 0
147: set gfortran = 0147: set gfortran = 0
148: set pathscale = 0148: set pathscale = 0
149: set ifort = 0149: set ifort = 0
150: set g95 = 0150: set g95 = 0
151: set pgf90 = 0151: set pgf95 = 0
152: set cfort = 0152: set cfort = 0
153: set f90 = 0153: set f90 = 0
154: set nih = 0154: set nih = 0
155: set tsri = 0155: set tsri = 0
156: set polyrate = 0156: set polyrate = 0
157: set full = 1157: set full = 1
158: set longint = 0158: set longint = 0
159: set itanium = 0159: set itanium = 0
160: set DEBUG = 0160: set DEBUG = 0
161: # PJ 06/2005161: # PJ 06/2005
294:     if ( $opt == 'F2C') set f2c = 1294:     if ( $opt == 'F2C') set f2c = 1
295:     if ( $opt == 'F77') set f77 = 1295:     if ( $opt == 'F77') set f77 = 1
296:     if ( $opt == 'IFC' || $opt == 'ifc') set ifc = 1296:     if ( $opt == 'IFC' || $opt == 'ifc') set ifc = 1
297:     if ( $opt == 'IFORT' || $opt == 'ifort') set ifort = 1297:     if ( $opt == 'IFORT' || $opt == 'ifort') set ifort = 1
298:     if ( $opt == 'G95' || $opt == 'g95') set g95 = 1298:     if ( $opt == 'G95' || $opt == 'g95') set g95 = 1
299:     if ( $opt == 'G77' || $opt == 'g77') set g77 = 1299:     if ( $opt == 'G77' || $opt == 'g77') set g77 = 1
300:     if ( $opt == 'EFC'|| $opt == 'efc' ) then300:     if ( $opt == 'EFC'|| $opt == 'efc' ) then
301:       set efc = 1301:       set efc = 1
302:       set longint = 1302:       set longint = 1
303:     endif303:     endif
304:     if ( $opt == 'PGF90' || $opt == 'pgf90' ) set pgf90 = 1304:     if ( $opt == 'PGF95' || $opt == 'pgf95' ) set pgf95 = 1
305:     if ( $opt == 'FORT') set cfort = 1305:     if ( $opt == 'FORT') set cfort = 1
306:     if ( $opt == 'F90'|| $opt == 'f90' ) set f90 = 1306:     if ( $opt == 'F90'|| $opt == 'f90' ) set f90 = 1
307:     if ( $opt == 'GFORTRAN'|| $opt == 'gfortran' ) set gfortran = 1307:     if ( $opt == 'GFORTRAN'|| $opt == 'gfortran' ) set gfortran = 1
308:     if ( $opt == 'X86_64'|| $opt == 'x86_64' ) set x86_64 = 1308:     if ( $opt == 'X86_64'|| $opt == 'x86_64' ) set x86_64 = 1
309:     if ( $opt == 'AMD64'|| $opt == 'amd64' ) set amd64 = 1309:     if ( $opt == 'AMD64'|| $opt == 'amd64' ) set amd64 = 1
310:     if ( $opt == 'PS'|| $opt == 'ps' ) set pathscale = 1310:     if ( $opt == 'PS'|| $opt == 'ps' ) set pathscale = 1
311:     if ( $opt == 'NIH'|| $opt == 'nih' ) set nih = 1311:     if ( $opt == 'NIH'|| $opt == 'nih' ) set nih = 1
312:     if ( $opt == 'NODISP'  ) set nodsp = 1312:     if ( $opt == 'NODISP'  ) set nodsp = 1
313:     if ( ( $opt == 'TSRI'|| $opt == 'tsri' ) && $full == 0 ) set tsri = 1313:     if ( ( $opt == 'TSRI'|| $opt == 'tsri' ) && $full == 0 ) set tsri = 1
314:     if ( $opt == 'big_e'|| $opt == 'BIG_E' ) set bigendian = 1314:     if ( $opt == 'big_e'|| $opt == 'BIG_E' ) set bigendian = 1
482:   else if ( $g77 == 1 ) then482:   else if ( $g77 == 1 ) then
483:     setenv GNU_G77 YES483:     setenv GNU_G77 YES
484:   else if ( $amd64 == 1 ) then484:   else if ( $amd64 == 1 ) then
485:     setenv AMD64 YES485:     setenv AMD64 YES
486:   else if ( $pathscale == 1 ) then486:   else if ( $pathscale == 1 ) then
487:     setenv PATHSCALE YES487:     setenv PATHSCALE YES
488:     if ( $x86_64 == 1 ) then488:     if ( $x86_64 == 1 ) then
489:       set longint = 1489:       set longint = 1
490:       setenv X86_64 YES490:       setenv X86_64 YES
491:     endif491:     endif
492:   else if ( $pgf90 == 1 ) then492:   else if ( $pgf95 == 1 ) then
493:     setenv PGI_F90 YES493:     setenv PGI_F95 YES
494:     if ( $x86_64 == 1 ) then494:     if ( $x86_64 == 1 ) then
495:       set longint = 1495:       set longint = 1
496:       setenv X86_64 YES496:       setenv X86_64 YES
497:     endif497:     endif
498:   else if ( $cfort == 1 ) then498:   else if ( $cfort == 1 ) then
499:     set longint=1499:     set longint=1
500:     setenv COMPAQ_FORT YES500:     setenv COMPAQ_FORT YES
501:   else501:   else
502:     setenv GFORTRAN YES502:     setenv GFORTRAN YES
503:     if ( $x86_64 == 1 ) then503:     if ( $x86_64 == 1 ) then
917: #    ----- IFORT ---------------------------------         917: #    ----- IFORT ---------------------------------         
918:         else if ( $ifort == 1 ) then918:         else if ( $ifort == 1 ) then
919:           if ( $mpiset == 1 && $scali == 0 ) then919:           if ( $mpiset == 1 && $scali == 0 ) then
920:             sed -e 's@GNU-Linux-compiler@mpif90 -O3 -tpp7 -axW@g' \920:             sed -e 's@GNU-Linux-compiler@mpif90 -O3 -tpp7 -axW@g' \
921:                     $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$921:                     $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$
922:           else922:           else
923:             sed -e 's@GNU-Linux-compiler@ifort -O3 -tpp7 -axW@g' \923:             sed -e 's@GNU-Linux-compiler@ifort -O3 -tpp7 -axW@g' \
924:                    $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$924:                    $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$
925:           endif925:           endif
926: #    ----- PGF77 ---------------------------------         926: #    ----- PGF77 ---------------------------------         
927:         else if ( $pgf90 == 1 ) then927:         else if ( $pgf95 == 1 ) then
928:           if ( $mpiset == 1 && $scali == 0 ) then928:           if ( $mpiset == 1 && $scali == 0 ) then
929:             sed -e \929:             sed -e \
930: #    's@GNU-Linux-compiler@mpif90 -O2 -Munroll -tp p7 -Mnoframe@g' \930:     's@GNU-Linux-compiler@mpif90 -O2 -Munroll -tp p7 -Mnoframe@g' \
931:     's@GNU-Linux-compiler@mpif90 -O2 -Munroll -Mnoframe@g' \ 
932:                    $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$931:                    $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$
933:           else932:           else
934:             sed -e \933:             sed -e \
935: #                's@GNU-Linux-compiler@pgf90 -O2 -Munroll -tp p7 -Mnoframe@g' \934:                 's@GNU-Linux-compiler@pgf95 -O2 -Munroll -tp p7 -Mnoframe@g' \
936:                 's@GNU-Linux-compiler@pgf90 -O2 -Munroll -Mnoframe@g' \ 
937:                     $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$935:                     $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$
938:           endif936:           endif
939:         else if ( $cfort == 1 ) then937:         else if ( $cfort == 1 ) then
940:             sed -e \938:             sed -e \
941:           's@GNU-Linux-compiler@fort -i8 -O4 -math_library fast -tune host@g' \939:           's@GNU-Linux-compiler@fort -i8 -O4 -math_library fast -tune host@g' \
942:                     $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$940:                     $chmtool/gmscomp_$chmhost > $chmtool/gmscomp_$$
943: #   ----- G95 ---------------------------------         941: #   ----- G95 ---------------------------------         
944:         else if ( $g95 == 1) then942:         else if ( $g95 == 1) then
945:           if ( $mpiset == 1 ) then943:           if ( $mpiset == 1 ) then
946:               sed -e 's@GNU-Linux-compiler@mpif90 -O3 @g' \944:               sed -e 's@GNU-Linux-compiler@mpif90 -O3 @g' \
1199:     breaksw1197:     breaksw
1200:   case gnu:1198:   case gnu:
1201:     set forcom = gfortran1199:     set forcom = gfortran
1202:     if ( $g95 == 1 ) set forcom = g951200:     if ( $g95 == 1 ) set forcom = g95
1203:     if ( $f77 == 1 ) set forcom = f771201:     if ( $f77 == 1 ) set forcom = f77
1204:     if ( $f2c == 1 ) set forcom = fort771202:     if ( $f2c == 1 ) set forcom = fort77
1205:     if ( $g77 == 1 ) set forcom = g771203:     if ( $g77 == 1 ) set forcom = g77
1206:     if ( $ifc == 1 ) set forcom = ifc1204:     if ( $ifc == 1 ) set forcom = ifc
1207:     if ( $efc == 1 ) set forcom = efc1205:     if ( $efc == 1 ) set forcom = efc
1208:     if ( $ifort == 1 ) set forcom = ifort1206:     if ( $ifort == 1 ) set forcom = ifort
1209:     if ( $pgf90 == 1 ) set forcom = pgf901207:     if ( $pgf95 == 1 ) set forcom = pgf95
1210:     if ( $cfort == 1 ) set forcom = fort1208:     if ( $cfort == 1 ) set forcom = fort
1211:     if ( $gfortran == 1 ) set forcom = "gfortran -static"1209:     if ( $gfortran == 1 ) set forcom = "gfortran -static"
1212:     if ( $pathscale == 1 ) set forcom = pathf901210:     if ( $pathscale == 1 ) set forcom = pathf90
1213:     if (! -e prefx_gnu  )   $forcom -o prefx_$$    prefx.f1211:     if (! -e prefx_gnu  )   $forcom -o prefx_$$    prefx.f
1214:     echo "GNU"      >! $chmbuild/pref$$.dat1212:     echo "GNU"      >! $chmbuild/pref$$.dat
1215:     if($g95 == 1) echo "G95" >> $chmbuild/pref$$.dat1213:     if($g95 == 1) echo "G95" >> $chmbuild/pref$$.dat
1216:     echo "UNIX"     >> $chmbuild/pref$$.dat1214:     echo "UNIX"     >> $chmbuild/pref$$.dat
1217:     echo "SCALAR"   >> $chmbuild/pref$$.dat1215:     echo "SCALAR"   >> $chmbuild/pref$$.dat
1218:     if ( $pvmset == 1 ) echo "PVMC"       >> $chmbuild/pref$$.dat1216:     if ( $pvmset == 1 ) echo "PVMC"       >> $chmbuild/pref$$.dat
1219:     if ( $socket == 1 ) echo "SOCKET"     >> $chmbuild/pref$$.dat1217:     if ( $socket == 1 ) echo "SOCKET"     >> $chmbuild/pref$$.dat
1532:     echo " "1530:     echo " "
1533:     echo " install.com> Unrecognized CHARMM size specified."1531:     echo " install.com> Unrecognized CHARMM size specified."
1534:     echo "              As default, MEDIUM will be used."1532:     echo "              As default, MEDIUM will be used."
1535:     echo "MEDIUM"   >> $chmbuild/pref$$.dat1533:     echo "MEDIUM"   >> $chmbuild/pref$$.dat
1536:     breaksw1534:     breaksw
1537: endsw1535: endsw
1538: #1536: #
1539: # generic preprocessing keys to build CHARMM1537: # generic preprocessing keys to build CHARMM
1540:   echo "EXPAND"     >> $chmbuild/pref$$.dat1538:   echo "EXPAND"     >> $chmbuild/pref$$.dat
1541:   echo "PUTFCM"     >> $chmbuild/pref$$.dat1539:   echo "PUTFCM"     >> $chmbuild/pref$$.dat
1542:   echo "FCMDIR=$chmsrc/fcm" >> $chmbuild/pref$$.dat1540:   echo "FCMDIR=fcm" >> $chmbuild/pref$$.dat
1543: #1541: #
1544: # specify parallel platform preprocessor keys1542: # specify parallel platform preprocessor keys
1545: # and supported features to non-parallel platforms1543: # and supported features to non-parallel platforms
1546: if ( $chm_host == alphamp ) then1544: if ( $chm_host == alphamp ) then
1547:   echo "NOCORREL"   >> $chmbuild/pref$$.dat1545:   echo "NOCORREL"   >> $chmbuild/pref$$.dat
1548:   echo "NOVIBRAN"   >> $chmbuild/pref$$.dat1546:   echo "NOVIBRAN"   >> $chmbuild/pref$$.dat
1549:   echo "NOST2"      >> $chmbuild/pref$$.dat1547:   echo "NOST2"      >> $chmbuild/pref$$.dat
1550: else if ( $chm_host == t3e ) then1548: else if ( $chm_host == t3e ) then
1551:   echo "NOCORREL"   >> $chmbuild/pref$$.dat1549:   echo "NOCORREL"   >> $chmbuild/pref$$.dat
1552:   echo "NOVIBRAN"   >> $chmbuild/pref$$.dat1550:   echo "NOVIBRAN"   >> $chmbuild/pref$$.dat
1870: #1868: #
1871: # PS check??1869: # PS check??
1872: if( $chm_host == gnu && ( $mpiset == 1 || $ensemble ==  1 ) && $scali == 0) then1870: if( $chm_host == gnu && ( $mpiset == 1 || $ensemble ==  1 ) && $scali == 0) then
1873:     mv Makefile_$$ Makefile.tmp_$$1871:     mv Makefile_$$ Makefile.tmp_$$
1874:     sed -e 's/g77/mpif90 -Impi/g' -e 's/ gcc/ mpicc/g' \1872:     sed -e 's/g77/mpif90 -Impi/g' -e 's/ gcc/ mpicc/g' \
1875:         -e 's/ fort/ mpif90 -fc=fort -Impi/g' -e 's/ccc/mpicc -cc=ccc/g' \1873:         -e 's/ fort/ mpif90 -fc=fort -Impi/g' -e 's/ccc/mpicc -cc=ccc/g' \
1876:         -e 's/ifc/mpif90 -Impi/g' \1874:         -e 's/ifc/mpif90 -Impi/g' \
1877:         -e 's/ gfortran/ mpif90 -Impi/g' \1875:         -e 's/ gfortran/ mpif90 -Impi/g' \
1878:         -e 's/ifort/mpif90 -Impi/g' \1876:         -e 's/ifort/mpif90 -Impi/g' \
1879:         -e 's/pathf90/mpif90 -Impi/g' \1877:         -e 's/pathf90/mpif90 -Impi/g' \
1880:         -e 's/pgf90/mpif90 -Impi/g' < Makefile.tmp_$$ > Makefile_$$1878:         -e 's/pgf95/mpif90 -Impi/g' < Makefile.tmp_$$ > Makefile_$$
1881:     /bin/rm Makefile.tmp_$$1879:     /bin/rm Makefile.tmp_$$
1882: else if( $chm_host == hpitanium && $mpiset == 1 ) then1880: else if( $chm_host == hpitanium && $mpiset == 1 ) then
1883:     mv Makefile_$$ Makefile.tmp_$$1881:     mv Makefile_$$ Makefile.tmp_$$
1884:     sed -e 's/f90/mpif90 -Impi/g' -e 's/cc/mpicc/g' \1882:     sed -e 's/f90/mpif90 -Impi/g' -e 's/cc/mpicc/g' \
1885:         < Makefile.tmp_$$ > Makefile_$$1883:         < Makefile.tmp_$$ > Makefile_$$
1886:     /bin/rm Makefile.tmp_$$1884:     /bin/rm Makefile.tmp_$$
1887: else if( $chm_host == itanium && $mpiset == 1 ) then1885: else if( $chm_host == itanium && $mpiset == 1 ) then
1888:     mv Makefile_$$ Makefile.tmp_$$1886:     mv Makefile_$$ Makefile.tmp_$$
1889:     sed -e 's/ifort/mpif90 -Impi/g' -e 's/icc/mpicc/g' \1887:     sed -e 's/ifort/mpif90 -Impi/g' -e 's/icc/mpicc/g' \
1890:         < Makefile.tmp_$$ > Makefile_$$1888:         < Makefile.tmp_$$ > Makefile_$$
2122:       setenv F77 absoft2120:       setenv F77 absoft
2123:     else if ( $f2c == 1 ) then2121:     else if ( $f2c == 1 ) then
2124:       # this will need work2122:       # this will need work
2125:       setenv F77 f2c2123:       setenv F77 f2c
2126:     else if ( $g77 == 1 ) then2124:     else if ( $g77 == 1 ) then
2127:       setenv F77 g772125:       setenv F77 g77
2128:     else if ( $ifc == 1 ) then2126:     else if ( $ifc == 1 ) then
2129:       setenv F77 ifc2127:       setenv F77 ifc
2130:     else if ( $efc == 1 ) then2128:     else if ( $efc == 1 ) then
2131:       setenv F77 efc2129:       setenv F77 efc
2132:     else if ( $pgf90 == 1 ) then2130:     else if ( $pgf95 == 1 ) then
2133:       setenv F77 pgf902131:       setenv F77 pgf95
2134:     else if ( $g95 == 1 ) then2132:     else if ( $g95 == 1 ) then
2135:       setenv F77 g952133:       setenv F77 g95
2136:     else if ( $gfortran == 1 ) then2134:     else if ( $gfortran == 1 ) then
2137:       setenv F77 gfortran2135:       setenv F77 gfortran
2138:     else if ( $cfort == 1 ) then2136:     else if ( $cfort == 1 ) then
2139:       # this will need work2137:       # this will need work
2140:       setenv F77 fort2138:       setenv F77 fort
2141:     else2139:     else
2142:       setenv F77 gfortran2140:       setenv F77 gfortran
2143:     endif2141:     endif
2259: if ( 1 == `grep EMAP pref.dat | wc -l` ) then2257: if ( 1 == `grep EMAP pref.dat | wc -l` ) then
2260:  set emapmk = '-f emap.mk'2258:  set emapmk = '-f emap.mk'
2261:  awk '/dynamc/ { print "        $(LIB)/emap.a \\" } \2259:  awk '/dynamc/ { print "        $(LIB)/emap.a \\" } \
2262:       /./      {print $0}' \2260:       /./      {print $0}' \
2263:                 objlibs.mk > objlibs.mk.tmp$$2261:                 objlibs.mk > objlibs.mk.tmp$$
2264:  if ( 1 != `grep emap objlibs.mk | wc -l` ) mv objlibs.mk.tmp$$ objlibs.mk2262:  if ( 1 != `grep emap objlibs.mk | wc -l` ) mv objlibs.mk.tmp$$ objlibs.mk
2265: else2263: else
2266:  set emapmk = ''2264:  set emapmk = ''
2267: endif2265: endif
2268: 2266: 
2269: # Setup CHEQ processing 
2270: if ( 1 == `grep CHEQ pref.dat | wc -l` ) then 
2271:  set cheq = '-f cheq.mk' 
2272:  awk '/cff/ { print "   $(LIB)/cheq.a \\" } \ 
2273:       /./      {print $0}' \ 
2274:                 Makefile > Makefile.tmp$$ 
2275:  mv Makefile.tmp$$  Makefile 
2276: else 
2277:  set cheq = '' 
2278: endif 
2279:  
2280: # Check for socket compile - some machines don't support this2267: # Check for socket compile - some machines don't support this
2281: # include only when really needed2268: # include only when really needed
2282: if ( 1 == `grep SOCKET pref.dat | wc -l` ) then2269: if ( 1 == `grep SOCKET pref.dat | wc -l` ) then
2283:      sed -e 's@cc @cc -Dchmsocket @' \2270:      sed -e 's@cc @cc -Dchmsocket @' \
2284:       Makefile > Makefile.tmp$$2271:       Makefile > Makefile.tmp$$
2285:  mv Makefile.tmp$$  Makefile2272:  mv Makefile.tmp$$  Makefile
2286: endif2273: endif
2287: 2274: 
2288: #2275: #
2289: # POLYRATE installation control2276: # POLYRATE installation control
2317: else2304: else
2318:    set prate = ''2305:    set prate = ''
2319: endif2306: endif
2320: #2307: #
2321: 2308: 
2322: if ( ! $?MAKE_COMMAND ) then2309: if ( ! $?MAKE_COMMAND ) then
2323:   setenv MAKE_COMMAND make2310:   setenv MAKE_COMMAND make
2324: endif2311: endif
2325: 2312: 
2326: set go_make = "$MAKE_COMMAND -f Makefile -f charmm.mk -f adumb.mk -f cadint.mk"2313: set go_make = "$MAKE_COMMAND -f Makefile -f charmm.mk -f adumb.mk -f cadint.mk"
2327: set go_make = "$go_make -f cff.mk -f cheq.mk -f correl.mk -f dimb.mk -f dynamc.mk"2314: set go_make = "$go_make -f cff.mk -f correl.mk -f dimb.mk -f dynamc.mk"
2328: set go_make = "$go_make -f energy.mk -f gamint.mk -f gukint.mk"2315: set go_make = "$go_make -f energy.mk -f gamint.mk -f gukint.mk"
2329: set go_make = "$go_make -f gener.mk -f image.mk -f io.mk -f machdep.mk "2316: set go_make = "$go_make -f gener.mk -f image.mk -f io.mk -f machdep.mk "
2330: set go_make = "$go_make -f mc.mk -f mmff.mk -f manip.mk -f minmiz.mk"2317: set go_make = "$go_make -f mc.mk -f mmff.mk -f manip.mk -f minmiz.mk"
2331: set go_make = "$go_make -f misc.mk -f mndint.mk -f molvib.mk -f mscale.mk"2318: set go_make = "$go_make -f misc.mk -f mndint.mk -f molvib.mk -f mscale.mk"
2332: set go_make = "$go_make -f nbonds.mk -f pert.mk -f quantum.mk -f rxncor.mk"2319: set go_make = "$go_make -f nbonds.mk -f pert.mk -f quantum.mk -f rxncor.mk"
2333: set go_make = "$go_make -f shapes.mk -f solvation.mk -f util.mk "2320: set go_make = "$go_make -f shapes.mk -f solvation.mk -f util.mk "
2334: set go_make = "$go_make -f vibran.mk -f zerom.mk"2321: set go_make = "$go_make -f vibran.mk -f zerom.mk"
2335: 2322: 
2336: set go_make = "$go_make $cadpac $emapmk $flucq $gamess $graphics $mbond"2323: set go_make = "$go_make $cadpac $emapmk $flucq $gamess $graphics $mbond"
2337: set go_make = "$go_make $mndo97 $squantm $pipf $prate $sccdftb"2324: set go_make = "$go_make $mndo97 $squantm $pipf $prate $sccdftb"


r21974/Makefile 2017-01-21 10:33:59.374250000 +0000 r21973/Makefile 2017-01-21 10:34:03.294250000 +0000
  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/Makefile' in revision 21974  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/Makefile' in revision 21973


r21974/Makefile_gnu 2017-01-21 10:33:56.466250000 +0000 r21973/Makefile_gnu 2017-01-21 10:34:00.014250000 +0000
 18: # GNU options for MBO(N)D code: 18: # GNU options for MBO(N)D code:
 19: # note that -fno-automatic -finit-local-zero are dangerous 19: # note that -fno-automatic -finit-local-zero are dangerous
 20: # options because they mask incorrect behaviour. They are 20: # options because they mask incorrect behaviour. They are
 21: # here only to simulate SGI compiler behaviour to facilitate 21: # here only to simulate SGI compiler behaviour to facilitate
 22: # initial porting from SGI (with -static) 22: # initial porting from SGI (with -static)
 23: # and eventually should be removed and offending code fixed. 23: # and eventually should be removed and offending code fixed.
 24: # NOTE: BLAS library is needed by MBO(N) code 24: # NOTE: BLAS library is needed by MBO(N) code
 25: # BLAS for Linux/PPro from: http://www.cs.utk.edu/~ghenry/distrib 25: # BLAS for Linux/PPro from: http://www.cs.utk.edu/~ghenry/distrib
 26: # BLAS dependency has been removed as of c27b2 and c28a2 26: # BLAS dependency has been removed as of c27b2 and c28a2
 27: ifdef MBOND 27: ifdef MBOND
 28: #MBONDFLG = -fno-automatic -finit-local-zero -O -Wall -Wsurprising -W -DCHARMM 28: MBONDFLG = -fno-automatic -finit-local-zero -O -Wall -Wsurprising -W -DCHARMM
 29: MBONDFLG = -O -Wall -Wsurprising -W -DCHARMM 
 30: INCLUDE = -I../../source/moldyn 29: INCLUDE = -I../../source/moldyn
 31: endif 30: endif
 32:  31: 
 33: # options for compilation with APBS 32: # options for compilation with APBS
 34: ifdef APBS 33: ifdef APBS
 35: # these must be defined: APBS_LIB, IAPBS_LIB and APBS_BLAS 34: # these must be defined: APBS_LIB, IAPBS_LIB and APBS_BLAS
 36: ADDLIB := $(ADDLIB) $(IAPBS_LIB)/iapbs.a \ 35: ADDLIB := $(ADDLIB) $(IAPBS_LIB)/iapbs.a \
 37:         -L$(APBS_LIB) -lapbsmainroutines -lapbs -lmaloc $(APBS_BLAS) 36:         -L$(APBS_LIB) -lapbsmainroutines -lapbs -lmaloc $(APBS_BLAS)
 38:  37: 
 39: endif 38: endif
 67: LD := g77 $(FFLAGS) 66: LD := g77 $(FFLAGS)
 68: endif 67: endif
 69:  68: 
 70: # Absoft f77 v3.4 69: # Absoft f77 v3.4
 71: ifdef ABSOFT_F77_V34 70: ifdef ABSOFT_F77_V34
 72: FC := f77 -N9 -O -U -f 71: FC := f77 -N9 -O -U -f
 73: LD := f77 -N9 -O -U -f 72: LD := f77 -N9 -O -U -f
 74: endif 73: endif
 75:  74: 
 76: # pgi pgf95 7.0 75: # pgi pgf95 7.0
 77: ifdef PGI_F90 76: ifdef PGI_F95
 78: ifdef BIG_ENDIAN 77: ifdef BIG_ENDIAN
 79: ENDIAN:= -Mbyteswapio 78: ENDIAN:= -Mbyteswapio
 80: endif 79: endif
 81: ifdef X86_64 80: ifdef X86_64
 82: #FC = pgf95 -mcmodel=medium -tp x64 $(ENDIAN) $(I8DUM1) 81: FC = pgf95 -mcmodel=medium -tp x64 $(ENDIAN) $(I8DUM1)
 83: #LD = pgf95 -mcmodel=medium -tp x64 -Munroll -Mnoframe $(ENDIAN) $(I8DUM1) 82: LD = pgf95 -mcmodel=medium -tp x64 -Munroll -Mnoframe $(ENDIAN) $(I8DUM1)
 84: FC = pgf90 -Mextend -O3 -Munroll -Mnoframe 
 85: LD = pgf90 -Mextend -O3 -Munroll -Mnoframe 
 86: else 83: else
 87: #FC = pgf95 -tp p7 $(ENDIAN) 84: FC = pgf95 -tp p7 $(ENDIAN)
 88: #LD = pgf95 -tp p7 -Munroll -Mnoframe $(ENDIAN) 85: LD = pgf95 -tp p7 -Munroll -Mnoframe $(ENDIAN)
 89: FC = pgf90 -Mextend -O3 -Munroll -Mnoframe 
 90: LD = pgf90 -Mextend -O3 -Munroll -Mnoframe 
 91: endif 86: endif
 92: endif 87: endif
 93:  88: 
 94: # Absoft  f77 4.4 (f90 v6.0): 89: # Absoft  f77 4.4 (f90 v6.0):
 95: # options: f lower case symbolic names, N3 include record length in unformatted 90: # options: f lower case symbolic names, N3 include record length in unformatted
 96: # files, N26 big endian unformatted files, B108 append underscore to external 91: # files, N26 big endian unformatted files, B108 append underscore to external
 97: # procedure names, B100 optimize for PentiumII, N86 optimize address expr., 92: # procedure names, B100 optimize for PentiumII, N86 optimize address expr.,
 98: # U enable loop unrolling 93: # U enable loop unrolling
 99: # 94: #
100: ifdef ABSOFT_F77_V44 95: ifdef ABSOFT_F77_V44
278:                 ifdef GFORTRAN273:                 ifdef GFORTRAN
279:                   FC0 = $(FC) -c -O0274:                   FC0 = $(FC) -c -O0
280:                   FC1 = $(FC) -c -O1275:                   FC1 = $(FC) -c -O1
281:                   FC2 = $(FC) -c -O2276:                   FC2 = $(FC) -c -O2
282:                   FC3 = $(FC) -c -O3277:                   FC3 = $(FC) -c -O3
283:                   FCR = $(FC) -c -g 278:                   FCR = $(FC) -c -g 
284:                   FCD = $(FC) -c -g -O0 279:                   FCD = $(FC) -c -g -O0 
285:                   FCRD = $(FC) -c -g -V -O0 280:                   FCRD = $(FC) -c -g -V -O0 
286:                   LDD = $(FC)281:                   LDD = $(FC)
287:                 else282:                 else
288: #                  FC0 = $(FC) -c -finit-local-zero -fno-automatic283:                   FC0 = $(FC) -c -finit-local-zero -fno-automatic
289:                   FC0 = $(FC) -c  
290:                 endif284:                 endif
291:               endif285:               endif
292:             endif286:             endif
293:           endif287:           endif
294:         endif288:         endif
295:       endif289:       endif
296:     endif290:     endif
297:   endif291:   endif
298: endif292: endif
299: 293: 


r21974/SVNREV 2017-01-21 10:33:59.714250000 +0000 r21973/SVNREV 2017-01-21 10:34:03.582250000 +0000
  1: 21969  1: 21954


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0