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