1 | !**************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the subroutines: hbond,chhb,ishybd,
|
---|
4 | ! ishybdo,nursat,interhbond
|
---|
5 | !
|
---|
6 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
7 | ! Shura Hayryan, Chin-Ku
|
---|
8 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
9 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
10 | !
|
---|
11 | ! **************************************************************
|
---|
12 |
|
---|
13 |
|
---|
14 | subroutine hbond(nml,mhb,ipr)
|
---|
15 | ! .................................................................
|
---|
16 | ! PURPOSE: find hydrogen bonds in molecule 'nml'
|
---|
17 | !
|
---|
18 | ! prints HBonds, if ipr > 0
|
---|
19 | !
|
---|
20 | ! OUTPUT: mhb - number of hyd.bds. of type i->i+4
|
---|
21 | !
|
---|
22 | ! to INCL.H:
|
---|
23 | !
|
---|
24 | ! ntyhb - number of different types of hyd. bds. found
|
---|
25 | ! nutyhb - number of hyd.bds. found for each type
|
---|
26 | ! ixtyhb - index for each type of hyd. bd. composed as
|
---|
27 | ! (atom idx. of H) * 1000 + atm.idx. of acceptor
|
---|
28 | !
|
---|
29 | ! CALLS: chhb,ishybd (ishybdo),nursat
|
---|
30 | !
|
---|
31 | !................................................................
|
---|
32 |
|
---|
33 | include 'INCL.H'
|
---|
34 |
|
---|
35 | integer nml
|
---|
36 |
|
---|
37 | integer mhb
|
---|
38 |
|
---|
39 | integer nursat
|
---|
40 | double precision atbase
|
---|
41 |
|
---|
42 | integer ipr, i2s, i, i14, i1, i1s, i2, ia, ixhb, id, ifivr, ii, ih
|
---|
43 | integer ims, io, iv, ivw, ix, k, jj, j, n, na, nd, ntlvr, iat
|
---|
44 |
|
---|
45 | !f2py intent(out) mhb
|
---|
46 | parameter (atbase=mxat)
|
---|
47 | logical ishb
|
---|
48 |
|
---|
49 |
|
---|
50 | do i=1,mxtyhb
|
---|
51 | nutyhb(i) = 0
|
---|
52 | ixtyhb(i) = 0
|
---|
53 | enddo
|
---|
54 | ntyhb=0
|
---|
55 | if (nml.eq.0) then
|
---|
56 | ntlvr = nvr
|
---|
57 | ifivr = ivrml1(1)
|
---|
58 | else
|
---|
59 | ntlvr=nvrml(nml)
|
---|
60 | if (ntlvr.eq.0) then
|
---|
61 | write (*,'(a,i4)')
|
---|
62 | & ' hbond> No variables defined in molecule #',nml
|
---|
63 | return
|
---|
64 | endif
|
---|
65 |
|
---|
66 | ifivr=ivrml1(nml)
|
---|
67 | ! Index of last moving set
|
---|
68 | i1s=imsml1(nml)+nmsml(nml)
|
---|
69 | endif
|
---|
70 | ! Loop over all variables
|
---|
71 | do io=ifivr+ntlvr-1,ifivr,-1
|
---|
72 | ! Get index of variable
|
---|
73 | iv=iorvr(io)
|
---|
74 | ! Index of next to last moving set
|
---|
75 | i2s=i1s-1
|
---|
76 | ! Index of moving set belonging to iv
|
---|
77 | i1s=imsvr1(iv)
|
---|
78 | ! Loop over all moving sets between the one belonging to iv and the
|
---|
79 | ! next to last one
|
---|
80 | do ims=i1s,i2s
|
---|
81 | ! First atom in moving set
|
---|
82 | i1=latms1(ims)
|
---|
83 | ! Last atom in moving set
|
---|
84 | i2=latms2(ims)
|
---|
85 | ! Loop over all atoms in moving set.
|
---|
86 | do i=i1,i2
|
---|
87 | ! Loop over van der Waals domains of atom i
|
---|
88 | do ivw=ivwat1(i),ivwat2(i)
|
---|
89 | ! Loop over atoms in van der Waals domain.
|
---|
90 | do j=lvwat1(ivw),lvwat2(ivw)
|
---|
91 |
|
---|
92 | call ishybd(i,j,ishb,ih,ia) ! Thornton criteria
|
---|
93 |
|
---|
94 | if (ishb) then
|
---|
95 |
|
---|
96 | ixhb=ih*atbase+ia
|
---|
97 |
|
---|
98 | do k=1,ntyhb
|
---|
99 | if (ixhb.eq.ixtyhb(k)) then
|
---|
100 | nutyhb(k)=nutyhb(k)+1
|
---|
101 | goto 1
|
---|
102 | endif
|
---|
103 | enddo
|
---|
104 |
|
---|
105 | if (ntyhb.lt.mxtyhb) then
|
---|
106 | ntyhb=ntyhb+1
|
---|
107 | nutyhb(ntyhb)=1
|
---|
108 | ixtyhb(ntyhb)=ixhb
|
---|
109 | else
|
---|
110 | write(*,*) ' hbond> increase parameter MXTYHB'
|
---|
111 | stop
|
---|
112 | endif
|
---|
113 |
|
---|
114 | endif ! have h.b.
|
---|
115 |
|
---|
116 | 1 enddo ! ... atoms j
|
---|
117 | enddo ! ... vdW-domains of i
|
---|
118 |
|
---|
119 | do i14=i14at1(i),i14at2(i) ! over 1-4 partners of 'i'
|
---|
120 | j=l14at(i14)
|
---|
121 |
|
---|
122 | call ishybd(i,j,ishb,ih,ia) ! Thornton criteria
|
---|
123 | ! call ishybdo(i,j,ishb,ih,ia)
|
---|
124 |
|
---|
125 | if (ishb) then
|
---|
126 |
|
---|
127 | ixhb=ih*atbase+ia
|
---|
128 |
|
---|
129 | do k=1,ntyhb
|
---|
130 | if (ixhb.eq.ixtyhb(k)) then
|
---|
131 | nutyhb(k)=nutyhb(k)+1
|
---|
132 | goto 2
|
---|
133 | endif
|
---|
134 | enddo
|
---|
135 |
|
---|
136 | if (ntyhb.lt.mxtyhb) then
|
---|
137 | ntyhb=ntyhb+1
|
---|
138 | nutyhb(ntyhb)=1
|
---|
139 | ixtyhb(ntyhb)=ixhb
|
---|
140 | else
|
---|
141 | write(*,*) ' hbond> increase parameter MXTYHB'
|
---|
142 | stop
|
---|
143 | endif
|
---|
144 |
|
---|
145 | endif ! have h.b.
|
---|
146 |
|
---|
147 | 2 enddo ! ... 1-4-partners of i
|
---|
148 |
|
---|
149 | enddo ! ... atoms i
|
---|
150 | enddo ! ... m.s.
|
---|
151 |
|
---|
152 | enddo ! ... variables
|
---|
153 |
|
---|
154 | mhb=0
|
---|
155 |
|
---|
156 | ! do inhb=1,ntyhb
|
---|
157 | ! mhb = mhb+nutyhb(inhb)
|
---|
158 | ! enddo
|
---|
159 |
|
---|
160 | if (ipr.gt.0) write(*,'(1x,a,/)') ' hbond> Hydrogen Bonds:'
|
---|
161 |
|
---|
162 | if (ntyhb.gt.0) then
|
---|
163 |
|
---|
164 | ii = 0
|
---|
165 | do i=1,ntyhb
|
---|
166 | jj = nutyhb(i)
|
---|
167 | do j = 1,jj
|
---|
168 | ii = ii + 1
|
---|
169 | ix = ixtyhb(ii)
|
---|
170 |
|
---|
171 | id =ix / atbase ! donor atom
|
---|
172 | nd = nursat(id) ! & residue
|
---|
173 |
|
---|
174 | ia = ix - id * atbase ! acceptor atom
|
---|
175 | na = nursat(ia)
|
---|
176 |
|
---|
177 | n = nd - na
|
---|
178 |
|
---|
179 | if (n.gt.4) mhb = mhb+1 ! only count these
|
---|
180 |
|
---|
181 | if (ipr.gt.0) then
|
---|
182 |
|
---|
183 | if (n.gt.0) then
|
---|
184 | write(*,'(1x,i3,a2,a4,a3,i3,1x,a4,a7,a4,a3,i3,1x,a4,a9,
|
---|
185 | & i2)')
|
---|
186 | & ii,') ',nmat(ia),' ( ',na,seq(na),' ) <-- ',nmat(id),
|
---|
187 | & ' ( ', nd,seq(nd),' ) = i +',n
|
---|
188 | else
|
---|
189 | write(*,'(1x,i3,a2,a4,a3,i3,1x,a4,a7,a4,a3,i3,1x,a4,a9,
|
---|
190 | & i2)')
|
---|
191 | & ii,') ',nmat(ia),' ( ',na,seq(na),' ) <-- ',nmat(id),
|
---|
192 | & ' ( ', nd,seq(nd),' ) = i -',abs(n)
|
---|
193 | endif
|
---|
194 |
|
---|
195 | call chhb(ia,id)
|
---|
196 | endif
|
---|
197 |
|
---|
198 | enddo
|
---|
199 | enddo
|
---|
200 | endif
|
---|
201 |
|
---|
202 | return
|
---|
203 | end
|
---|
204 | ! .....................................................................
|
---|
205 | ! Calculates hydrogen bonds between different chains.
|
---|
206 | !
|
---|
207 | ! @return number of intermolecular hydrogen bonds. Returns 0 if only
|
---|
208 | ! one molecule is present. The value is returned in the
|
---|
209 | ! variable mhb.
|
---|
210 | !
|
---|
211 | ! @author Jan H. Meinke <j.meinke@fz-juelich.de>
|
---|
212 | !
|
---|
213 | ! .....................................................................
|
---|
214 | subroutine interhbond(mhb)
|
---|
215 |
|
---|
216 | include 'INCL.H'
|
---|
217 |
|
---|
218 | !f2py intent(out) mhb
|
---|
219 |
|
---|
220 | logical ishb
|
---|
221 |
|
---|
222 | integer*4 mhb
|
---|
223 | integer iml, ires, jml, jat, jres
|
---|
224 |
|
---|
225 | integer ia, iat, ih
|
---|
226 |
|
---|
227 | mhb = 0
|
---|
228 | do iml = 1, ntlml
|
---|
229 | do jml = iml + 1, ntlml
|
---|
230 | mmhb(iml, jml) = 0
|
---|
231 | do ires= irsml1(iml), irsml2(iml)
|
---|
232 | do jres= irsml1(jml), irsml2(jml)
|
---|
233 | do iat = iatrs1(ires), iatrs2(ires)
|
---|
234 | do jat = iatrs1(jres), iatrs2(jres)
|
---|
235 | call ishybd(iat,jat,ishb,ih,ia)
|
---|
236 | if (ishb) then
|
---|
237 | mhb = mhb + 1
|
---|
238 | mmhb(iml, jml) = mmhb(iml, jml) + 1
|
---|
239 | endif
|
---|
240 | enddo ! jat
|
---|
241 | enddo ! iat
|
---|
242 | enddo ! jres
|
---|
243 | enddo ! ires
|
---|
244 | mmhb(jml, iml) = mmhb(iml, jml)
|
---|
245 | enddo ! jml
|
---|
246 | enddo ! iml
|
---|
247 |
|
---|
248 | end ! subroutine interhbond
|
---|
249 | ! ************************
|
---|
250 | subroutine chhb(i,j)
|
---|
251 |
|
---|
252 | include 'INCL.H'
|
---|
253 | integer i
|
---|
254 |
|
---|
255 | integer j
|
---|
256 |
|
---|
257 | integer ih
|
---|
258 |
|
---|
259 | integer ia, ib, id, ihb
|
---|
260 |
|
---|
261 | double precision valang
|
---|
262 |
|
---|
263 | double precision cdah, cdad, adab, adha, ahab, dah, dad
|
---|
264 |
|
---|
265 | ihb = ihbty(ityat(i),ityat(j))
|
---|
266 |
|
---|
267 | if (ihb.gt.0) then
|
---|
268 | ih=i
|
---|
269 | ia=j
|
---|
270 | else
|
---|
271 | ih=j
|
---|
272 | ia=i
|
---|
273 | endif
|
---|
274 |
|
---|
275 | dah=sqrt((xat(ih)-xat(ia))**2+(yat(ih)-yat(ia))**2+
|
---|
276 | & (zat(ih)-zat(ia))**2)
|
---|
277 |
|
---|
278 | id=iowat(ih)
|
---|
279 |
|
---|
280 | dad=sqrt((xat(id)-xat(ia))**2+(yat(id)-yat(ia))**2+
|
---|
281 | & (zat(id)-zat(ia))**2)
|
---|
282 | adha=valang(id,ih,ia)*crd
|
---|
283 |
|
---|
284 | ib=iowat(ia)
|
---|
285 |
|
---|
286 | ahab=valang(ih,ia,ib)*crd
|
---|
287 | adab=valang(id,ia,ib)*crd
|
---|
288 |
|
---|
289 | write(*,*) ' '
|
---|
290 | write(*,*) ' Dah: ',dah,' Dad: ',dad
|
---|
291 | write(*,*) ' Adha: ',adha,' Ahab: ',adha,' Adab: ',adab
|
---|
292 | write(*,*) ' '
|
---|
293 |
|
---|
294 | return
|
---|
295 | end
|
---|
296 | ! *************************************
|
---|
297 | subroutine ishybd(i,j,ishb,ih,ia)
|
---|
298 |
|
---|
299 |
|
---|
300 | ! ..........................................................
|
---|
301 | ! PURPOSE: checks for hydrogen bond between atoms 'i' & 'j'
|
---|
302 | ! according to geometric criteria
|
---|
303 | !
|
---|
304 | ! OUTPUT: logical 'ishb' - true, if have Hydrogen bond
|
---|
305 | ! ih - index of Hydrogen atom
|
---|
306 | ! ia - index of Acceptor atom
|
---|
307 | !
|
---|
308 | ! [I.K.McDonald,J.M.Thornton,Satisfying hydrogen bond
|
---|
309 | ! potential in proteins.J.Mol.Biol.238(5),777-793 (1994)]
|
---|
310 | !
|
---|
311 | ! D: hydrogen(=H) donor, A: acceptor, B: atom bound to A
|
---|
312 | !
|
---|
313 | ! Dis_HA <= 2.5 & Dis_DA <= 3.9 & Angle(D-H-A) > 90 &
|
---|
314 | ! Angle(H-A-B) > 90 & Angle(D-A-B) > 90
|
---|
315 | ! ..........................................................
|
---|
316 |
|
---|
317 | include 'INCL.H'
|
---|
318 | double precision cdah, cang, cahb, cdad
|
---|
319 | integer i
|
---|
320 |
|
---|
321 | integer j
|
---|
322 |
|
---|
323 | integer ih
|
---|
324 |
|
---|
325 | integer ia, ib, id, ihb
|
---|
326 |
|
---|
327 | double precision valang
|
---|
328 |
|
---|
329 | parameter (cdad=3.9d0,
|
---|
330 | & cdah=2.5d0,
|
---|
331 | & cang=110.d0)
|
---|
332 | ! # cang=90.d0)
|
---|
333 |
|
---|
334 | logical ishb
|
---|
335 |
|
---|
336 | cahb = cang * cdr
|
---|
337 | ishb = .false.
|
---|
338 |
|
---|
339 |
|
---|
340 | if (i.le.0.or.j.le.0) return
|
---|
341 |
|
---|
342 | ihb = ihbty(ityat(i),ityat(j))
|
---|
343 |
|
---|
344 | if (ihb.eq.0) then
|
---|
345 | return
|
---|
346 | elseif (ihb.gt.0) then
|
---|
347 | ih=i
|
---|
348 | ia=j
|
---|
349 | else
|
---|
350 | ih=j
|
---|
351 | ia=i
|
---|
352 | endif
|
---|
353 |
|
---|
354 | if (sqrt((xat(ih)-xat(ia))**2+(yat(ih)-yat(ia))**2+
|
---|
355 | & (zat(ih)-zat(ia))**2).gt.cdah) return
|
---|
356 |
|
---|
357 | id=iowat(ih)
|
---|
358 |
|
---|
359 | if (id.le.0.or.sqrt((xat(id)-xat(ia))**2+(yat(id)-yat(ia))**2+
|
---|
360 | & (zat(id)-zat(ia))**2).gt.cdad
|
---|
361 | & .or.valang(id,ih,ia).lt.cahb) return
|
---|
362 |
|
---|
363 | ib=iowat(ia)
|
---|
364 |
|
---|
365 | if (ib.gt.0.and.valang(ih,ia,ib).ge.cahb
|
---|
366 | & .and.valang(id,ia,ib).ge.cahb) ishb=.true.
|
---|
367 |
|
---|
368 | return
|
---|
369 | end
|
---|
370 | ! **************************************
|
---|
371 |
|
---|
372 | subroutine ishybdo(i,j,ishb,ih,ia)
|
---|
373 |
|
---|
374 | ! ..........................................................
|
---|
375 | ! PURPOSE: checks for hydrogen bond between atoms 'i' & 'j'
|
---|
376 | ! according to geometric criteria
|
---|
377 | !
|
---|
378 | ! OUTPUT: logical 'ishb' - true, if have Hydrogen bond
|
---|
379 | ! ih - index of Hydrogen atom
|
---|
380 | ! ia - index of Acceptor atom
|
---|
381 | !
|
---|
382 | ! D: hydrogen(=H) donor, A: acceptor
|
---|
383 | !
|
---|
384 | ! Dis_AH <= 2.5 & Angle(D-H-A) >= 160
|
---|
385 | ! ...........................................................
|
---|
386 |
|
---|
387 | include 'INCL.H'
|
---|
388 | integer i
|
---|
389 |
|
---|
390 | integer j
|
---|
391 |
|
---|
392 | integer ih
|
---|
393 |
|
---|
394 | double precision cahb, cdah, cang, valang
|
---|
395 |
|
---|
396 | integer ia, ihb, id
|
---|
397 | parameter (cdah=2.5d0,
|
---|
398 | & cang=140.d0)
|
---|
399 |
|
---|
400 | logical ishb
|
---|
401 |
|
---|
402 | cahb = cang * cdr
|
---|
403 |
|
---|
404 |
|
---|
405 | ishb = .false.
|
---|
406 |
|
---|
407 | if (i.le.0.or.j.le.0) return
|
---|
408 |
|
---|
409 | ihb = ihbty(ityat(i),ityat(j))
|
---|
410 |
|
---|
411 | if (ihb.eq.0) then
|
---|
412 | return
|
---|
413 | elseif (ihb.gt.0) then
|
---|
414 | ih=i
|
---|
415 | ia=j
|
---|
416 | else
|
---|
417 | ih=j
|
---|
418 | ia=i
|
---|
419 | endif
|
---|
420 |
|
---|
421 | if (sqrt((xat(ih)-xat(ia))**2+(yat(ih)-yat(ia))**2+
|
---|
422 | & (zat(ih)-zat(ia))**2).gt.cdah) return
|
---|
423 |
|
---|
424 | id=iowat(ih)
|
---|
425 |
|
---|
426 | if (id.gt.0.and.valang(id,ih,ia).ge.cahb) ishb=.true.
|
---|
427 |
|
---|
428 | return
|
---|
429 | end
|
---|
430 |
|
---|