1 | ! *******************************************************************
|
---|
2 | ! SMMP version of Anders Irback's force field, to be called the Lund
|
---|
3 | ! force field. This file contains the function enylun, which in turn
|
---|
4 | ! calls all the terms in the energy function. The terms Bias (ebias),
|
---|
5 | ! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the
|
---|
6 | ! Excluded volume (eexvol and eloexv) are also implemented in this
|
---|
7 | ! file.
|
---|
8 | !
|
---|
9 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
10 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
11 | !
|
---|
12 | subroutine set_local_constr(ires,frst,lst,root)
|
---|
13 | include 'INCL.H'
|
---|
14 | include 'incl_lund.h'
|
---|
15 | integer iat1, iat2
|
---|
16 |
|
---|
17 | integer ires,frst,lst,root
|
---|
18 | do iat1=iCa(ires)+frst,iCa(ires)+lst
|
---|
19 | do iat2=iat1+1,iCa(ires)+lst
|
---|
20 | matcon(iat2-iat1,iat1)=0
|
---|
21 | matcon(iat1-iat2,iat2)=0
|
---|
22 | enddo
|
---|
23 | enddo
|
---|
24 | iat1=iCa(ires)+root
|
---|
25 | do iat2=iCa(ires)+root,iCa(ires)+lst
|
---|
26 | matcon(iat2-iat1,iat1)=0
|
---|
27 | matcon(iat1-iat2,iat2)=0
|
---|
28 | enddo
|
---|
29 | end subroutine set_local_constr
|
---|
30 |
|
---|
31 | subroutine init_lundff
|
---|
32 | include 'INCL.H'
|
---|
33 | include 'incl_lund.h'
|
---|
34 | double precision eunit, csacc
|
---|
35 |
|
---|
36 | integer i, j, iml, iat1, iat2, k, iat3, l, iat4, m, iat3p, irs
|
---|
37 | integer iatoff, iatmrg, ilp, lci
|
---|
38 |
|
---|
39 | character mynm*4
|
---|
40 | logical prlvr
|
---|
41 |
|
---|
42 | print *,'initializing Lund forcefield'
|
---|
43 | ! Some parameters in the Lund force field.
|
---|
44 | ! The correspondence between internal energy scale and kcal/mol
|
---|
45 | eunit=1.3315d0
|
---|
46 | ! eunit=1.0
|
---|
47 | ! Bias
|
---|
48 | kbias=100.0d0*eunit
|
---|
49 | ! print *,'Bias'
|
---|
50 | ! Hydrogen bonds
|
---|
51 | epshb1=3.1d0*eunit
|
---|
52 | epshb2=2.0d0*eunit
|
---|
53 | sighb=2.0d0
|
---|
54 | cthb=4.5d0
|
---|
55 | cthb2=cthb*cthb
|
---|
56 | powa=0.5d0
|
---|
57 | powb=0.5d0
|
---|
58 | blhb=-30.0d0*(((sighb/cthb)**10-(sighb/cthb)**12))/cthb2
|
---|
59 | alhb=-(5*((sighb/cthb)**12)-6*((sighb/cthb)**10))-blhb*cthb2
|
---|
60 | sighb2=sighb*sighb
|
---|
61 | cdon=1.0d0
|
---|
62 | cacc=(1.0d0/1.23d0)**powb
|
---|
63 | csacc=(1.0d0/1.25d0)**powb
|
---|
64 | ! print *,'Hydrogen bonds'
|
---|
65 | ! Hydrophobicity
|
---|
66 | ! print *,'Hydrophobicity with nhptyp = ',nhptyp
|
---|
67 |
|
---|
68 | hpstrg(1)=0.0d0*eunit
|
---|
69 | hpstrg(2)=0.1d0*eunit
|
---|
70 | hpstrg(3)=0.1d0*eunit
|
---|
71 | hpstrg(4)=0.1d0*eunit
|
---|
72 | hpstrg(5)=0.9d0*eunit
|
---|
73 | hpstrg(6)=2.8d0*eunit
|
---|
74 | hpstrg(7)=0.1d0*eunit
|
---|
75 | hpstrg(8)=2.8d0*eunit
|
---|
76 | hpstrg(9)=3.2d0*eunit
|
---|
77 |
|
---|
78 | do i=1,mxrs
|
---|
79 | do j=1,6
|
---|
80 | ihpat(i,j)=0
|
---|
81 | enddo
|
---|
82 | nhpat(i)=0
|
---|
83 | enddo
|
---|
84 | do i=irsml1(1),irsml2(ntlml)
|
---|
85 | mynm=seq(i)
|
---|
86 | call tolost(mynm)
|
---|
87 | if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
|
---|
88 | & .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
|
---|
89 | & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
|
---|
90 | prlvr=.true. ! residue i is a proline variant
|
---|
91 | print *, 'proline variant ',mynm,i
|
---|
92 | else
|
---|
93 | prlvr=.false.
|
---|
94 | endif
|
---|
95 |
|
---|
96 | if (mynm.eq.'ala') then
|
---|
97 | nhpat(i)=1
|
---|
98 | ihpat(i,1)=iCa(i)+2
|
---|
99 | else if (mynm.eq.'val') then
|
---|
100 | nhpat(i)=3
|
---|
101 | ihpat(i,1)=iCa(i)+2
|
---|
102 | ihpat(i,2)=iCa(i)+4
|
---|
103 | ihpat(i,3)=iCa(i)+8
|
---|
104 | else if (prlvr) then
|
---|
105 | nhpat(i)=3
|
---|
106 | ihpat(i,1)=iCa(i)+2
|
---|
107 | ihpat(i,2)=iCa(i)+5
|
---|
108 | ihpat(i,3)=iCa(i)+8
|
---|
109 | else if (mynm.eq.'leu') then
|
---|
110 | nhpat(i)=4
|
---|
111 | ihpat(i,1)=iCa(i)+2
|
---|
112 | ihpat(i,2)=iCa(i)+5
|
---|
113 | ihpat(i,3)=iCa(i)+7
|
---|
114 | ihpat(i,4)=iCa(i)+11
|
---|
115 | else if (mynm.eq.'ile') then
|
---|
116 | nhpat(i)=4
|
---|
117 | ihpat(i,1)=iCa(i)+2
|
---|
118 | ihpat(i,2)=iCa(i)+4
|
---|
119 | ihpat(i,3)=iCa(i)+8
|
---|
120 | ihpat(i,4)=iCa(i)+11
|
---|
121 | else if (mynm.eq.'met') then
|
---|
122 | nhpat(i)=4
|
---|
123 | ihpat(i,1)=iCa(i)+2
|
---|
124 | ihpat(i,2)=iCa(i)+5
|
---|
125 | ihpat(i,3)=iCa(i)+8
|
---|
126 | ihpat(i,4)=iCa(i)+9
|
---|
127 | else if (mynm.eq.'phe') then
|
---|
128 | nhpat(i)=6
|
---|
129 | ihpat(i,1)=iCa(i)+5
|
---|
130 | ihpat(i,2)=iCa(i)+6
|
---|
131 | ihpat(i,3)=iCa(i)+8
|
---|
132 | ihpat(i,4)=iCa(i)+10
|
---|
133 | ihpat(i,5)=iCa(i)+12
|
---|
134 | ihpat(i,6)=iCa(i)+14
|
---|
135 | else if (mynm.eq.'tyr') then
|
---|
136 | nhpat(i)=6
|
---|
137 | ihpat(i,1)=iCa(i)+5
|
---|
138 | ihpat(i,2)=iCa(i)+6
|
---|
139 | ihpat(i,3)=iCa(i)+8
|
---|
140 | ihpat(i,4)=iCa(i)+10
|
---|
141 | ihpat(i,5)=iCa(i)+13
|
---|
142 | ihpat(i,6)=iCa(i)+15
|
---|
143 | else if (mynm.eq.'trp') then
|
---|
144 | nhpat(i)=6
|
---|
145 | ihpat(i,1)=iCa(i)+10
|
---|
146 | ihpat(i,2)=iCa(i)+11
|
---|
147 | ihpat(i,3)=iCa(i)+13
|
---|
148 | ihpat(i,4)=iCa(i)+15
|
---|
149 | ihpat(i,5)=iCa(i)+17
|
---|
150 | ihpat(i,6)=iCa(i)+19
|
---|
151 | endif
|
---|
152 | enddo
|
---|
153 | ! print *,'Hydrophobicity'
|
---|
154 |
|
---|
155 | ! Excluded volume and local pair excluded volume terms
|
---|
156 | exvk=0.1*eunit
|
---|
157 | exvcut=4.3
|
---|
158 | exvcut2=exvcut*exvcut
|
---|
159 |
|
---|
160 | sigsa(1)=1.0d0 ! hydrogen
|
---|
161 | sigsa(2)=1.0d0 ! hydrogen
|
---|
162 | sigsa(3)=1.0d0 ! hydrogen
|
---|
163 | sigsa(4)=1.0d0 ! hydrogen
|
---|
164 | sigsa(5)=1.0d0 ! hydrogen
|
---|
165 | sigsa(6)=1.0d0 ! hydrogen
|
---|
166 | sigsa(7)=1.75d0 ! carbon
|
---|
167 | sigsa(8)=1.75d0 ! carbon
|
---|
168 | sigsa(9)=1.75d0 ! carbon
|
---|
169 | sigsa(10)=1.42d0 ! oxygen
|
---|
170 | sigsa(11)=1.42d0 ! oxygen
|
---|
171 | sigsa(12)=1.42d0 ! oxygen
|
---|
172 | sigsa(13)=1.55d0 ! nitrogen
|
---|
173 | sigsa(14)=1.55d0 ! nitrogen
|
---|
174 | sigsa(15)=1.55d0 ! nitrogen
|
---|
175 | sigsa(16)=1.77d0 ! sulfur
|
---|
176 | sigsa(17)=1.0d0 ! hydrogen
|
---|
177 | sigsa(18)=1.75d0 ! carbon
|
---|
178 |
|
---|
179 | do i=1,mxtyat
|
---|
180 | do j=1,mxtyat
|
---|
181 | sig2lcp(i,j)=(sigsa(i)+sigsa(j))**2
|
---|
182 | asalcp(i,j)=-7*((sig2lcp(i,j)/exvcut2)**6.0d0)
|
---|
183 | bsalcp(i,j)=6*((sig2lcp(i,j)/exvcut2)**6.0d0)/exvcut2
|
---|
184 | enddo
|
---|
185 | enddo
|
---|
186 | ! print *,'Local pair excluded volume constants'
|
---|
187 |
|
---|
188 | exvlam=0.75d0
|
---|
189 | exvcutg=exvcut*exvlam
|
---|
190 | exvcutg2=exvcutg*exvcutg
|
---|
191 |
|
---|
192 | do i=1,mxtyat
|
---|
193 | do j=1,mxtyat
|
---|
194 | sig2exv(i,j)=(exvlam*(sigsa(i)+sigsa(j)))**2
|
---|
195 | asaexv(i,j)=-7*((sig2exv(i,j)/exvcutg2)**6.0)
|
---|
196 | bsaexv(i,j)=6*((sig2exv(i,j)/exvcutg2)**6.0)/exvcutg2
|
---|
197 | enddo
|
---|
198 | enddo
|
---|
199 | ! print *,'General excluded volume constants'
|
---|
200 |
|
---|
201 | ! Initialization of the connections matrix matcon(i,j). The index
|
---|
202 | ! i runs from -mxconr to +mxconr, and j from 1 to mxat.
|
---|
203 | ! matcon(i2-i1,i1) = 0, if the distance between atoms i1 and i2 is fixed
|
---|
204 | ! = 2, if atoms i1 and i2 are separated by 3 covalent
|
---|
205 | ! bonds and their distance can change
|
---|
206 | ! = 1, for all other pairs
|
---|
207 | ! if abs(i2-i1) > mxconr, the atoms are assumed to be separated by
|
---|
208 | ! many bonds, and with no restriction on their distances. On a protein
|
---|
209 | ! molecule made of natural amino acids, atoms with indices separated
|
---|
210 | ! by more than 35 can not be connected by three covalent bonds.
|
---|
211 |
|
---|
212 | do i=1,mxat
|
---|
213 | do j=-mxconr,mxconr
|
---|
214 | matcon(j,i)=1
|
---|
215 | enddo
|
---|
216 | matcon(0,i)=0
|
---|
217 | enddo
|
---|
218 | ! continued...
|
---|
219 | do iml=1,ntlml
|
---|
220 | do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
|
---|
221 | do j=1,nbdat(iat1)
|
---|
222 | iat2=ibdat(j,iat1)
|
---|
223 | matcon(iat2-iat1,iat1)=0 ! directly bonded atoms
|
---|
224 | matcon(iat1-iat2,iat2)=0 !
|
---|
225 | do k=1,nbdat(iat2)
|
---|
226 | iat3=ibdat(k,iat2)
|
---|
227 | matcon(iat3-iat1,iat1)=0 !
|
---|
228 | matcon(iat1-iat3,iat3)=0 ! 2 covalent bonds
|
---|
229 | do l=1,nbdat(iat3)
|
---|
230 | iat4=ibdat(l,iat3)
|
---|
231 | if (iat4.ne.iat2) then
|
---|
232 | matcon(iat4-iat1,iat1)=2 !
|
---|
233 | matcon(iat1-iat4,iat4)=2 ! 3 covalent bonds
|
---|
234 | endif
|
---|
235 | do m=1,nbdat(iat2)
|
---|
236 | iat3p=ibdat(m,iat2)
|
---|
237 | if (iat3p.ne.iat3) then
|
---|
238 | matcon(iat4-iat3p,iat3p)=2 ! 3 covalent bonds
|
---|
239 | matcon(iat3p-iat4,iat4)=2 !
|
---|
240 | endif
|
---|
241 | enddo
|
---|
242 | enddo
|
---|
243 | enddo
|
---|
244 | do k=1,nbdat(iat1)
|
---|
245 | iat3=ibdat(k,iat1)
|
---|
246 | matcon(iat3-iat2,iat2)=0 ! also 2 bonds iat2-iat1-iat3
|
---|
247 | matcon(iat2-iat3,iat3)=0 !
|
---|
248 | enddo
|
---|
249 | enddo
|
---|
250 | enddo
|
---|
251 |
|
---|
252 | ! print *,'going to initialize connections for first residue'
|
---|
253 | ! print *,'iN,iCa,iC =',iN(irsml1(iml)),
|
---|
254 | ! # iCa(irsml1(iml)),iC(irsml1(iml))
|
---|
255 | do iat1=iN(irsml1(iml))+1,iCa(irsml1(iml))-1
|
---|
256 | ! print *,'connections for iat1 = ',iat1
|
---|
257 | matcon(iat1-iN(irsml1(iml)),iN(irsml1(iml)))=0
|
---|
258 | matcon(iN(irsml1(iml))-iat1,iat1)=0
|
---|
259 | matcon(iat1-iCa(irsml1(iml)),iCa(irsml1(iml)))=0
|
---|
260 | matcon(iCa(irsml1(iml))-iat1,iat1)=0
|
---|
261 | do iat2=iat1,iCa(irsml1(iml))-1
|
---|
262 | matcon(iat2-iat1,iat1)=0
|
---|
263 | matcon(iat1-iat2,iat2)=0
|
---|
264 | enddo
|
---|
265 | matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2
|
---|
266 | matcon(iCa(irsml1(iml))+1-iat1,iat1)=2
|
---|
267 | matcon(iat1-iCa(irsml1(iml))-2,iCa(irsml1(iml))+2)=2
|
---|
268 | matcon(iCa(irsml1(iml))+2-iat1,iat1)=2
|
---|
269 | matcon(iat1-iC(irsml1(iml)),iC(irsml1(iml)))=2
|
---|
270 | matcon(iC(irsml1(iml))-iat1,iat1)=2
|
---|
271 | enddo
|
---|
272 |
|
---|
273 | ! Below: for certain residues, some atoms separated by 3 or more bonds
|
---|
274 | ! do not change distance. So, the connection matrix term for such pairs
|
---|
275 | ! should be zero.
|
---|
276 |
|
---|
277 | do irs=irsml1(iml),irsml2(iml)
|
---|
278 | if (irs.eq.irsml1(iml)) then
|
---|
279 | iatoff=1
|
---|
280 | else
|
---|
281 | iatoff=0
|
---|
282 | endif
|
---|
283 | if (irs.eq.irsml2(iml)) then
|
---|
284 | iatmrg=2
|
---|
285 | else
|
---|
286 | iatmrg=0
|
---|
287 | endif
|
---|
288 | mynm=seq(irs)
|
---|
289 | call tolost(mynm)
|
---|
290 | if ((mynm.eq.'pro')) then
|
---|
291 | prlvr=.true. ! residue i is a proline variant
|
---|
292 | else
|
---|
293 | prlvr=.false.
|
---|
294 | endif
|
---|
295 | if ((mynm.eq.'asn')) then
|
---|
296 | call set_local_constr(irs,5,9,2)
|
---|
297 | else if ((mynm.eq.'gln')) then
|
---|
298 | call set_local_constr(irs,8,12,5)
|
---|
299 | else if ((mynm.eq.'arg')) then
|
---|
300 | call set_local_constr(irs,11,19,8)
|
---|
301 | else if ((mynm.eq.'his')) then
|
---|
302 | call set_local_constr(irs,5,12,2)
|
---|
303 | else if (mynm.eq.'phe') then
|
---|
304 | call set_local_constr(irs,5,15,2)
|
---|
305 | else if (mynm.eq.'tyr') then
|
---|
306 | call set_local_constr(irs,5,16,2)
|
---|
307 | iat1=iCa(irs)+12 ! H_h
|
---|
308 | iat2=iCa(irs)+13 ! C_e2
|
---|
309 | iat3=iCa(irs)+8 ! C_e1
|
---|
310 | matcon(iat2-iat1,iat1)=2
|
---|
311 | matcon(iat1-iat2,iat2)=2
|
---|
312 | matcon(iat3-iat1,iat1)=2
|
---|
313 | matcon(iat1-iat3,iat3)=2
|
---|
314 | else if (mynm.eq.'trp') then
|
---|
315 | call set_local_constr(irs,5,19,2)
|
---|
316 | else if (prlvr) then
|
---|
317 | ! Proline. Many more distances are fixed because of the fixed
|
---|
318 | ! phi angle
|
---|
319 | call set_local_constr(irs,-3,11,-3)
|
---|
320 | matcon(iCa(irs-1)-iCa(irs)-8,iCa(irs)+8)=0
|
---|
321 | matcon(iCa(irs)+8-iCa(irs-1),iCa(irs-1))=0
|
---|
322 | endif
|
---|
323 | enddo
|
---|
324 |
|
---|
325 | ! Distances fixed because of the constant omega angle
|
---|
326 | do irs=irsml1(iml),irsml2(iml)-1
|
---|
327 | matcon(iCa(irs+1)-iC(irs)-1,iC(irs)+1)=0 ! O_i -- Ca_{i+1}
|
---|
328 | matcon(-iCa(irs+1)+iC(irs)+1,iCa(irs+1))=0 ! O_i -- Ca_{i+1}
|
---|
329 |
|
---|
330 | matcon(iN(irs+1)-iC(irs),iC(irs)+1)=0 ! O_i -- H_{i+1}
|
---|
331 | matcon(-iN(irs+1)+iC(irs),iN(irs+1)+1)=0 ! O_i -- H_{i+1}
|
---|
332 |
|
---|
333 | matcon(iCa(irs+1)-iCa(irs),iCa(irs))=0 ! Ca_i -- Ca_{i+1}
|
---|
334 | matcon(-iCa(irs+1)+iCa(irs),iCa(irs+1))=0 ! Ca_i -- Ca_{i+1}
|
---|
335 |
|
---|
336 | matcon(iN(irs+1)+1-iCa(irs),iCa(irs))=0 ! Ca_i -- H_{i+1}
|
---|
337 | matcon(-iN(irs+1)-1+iCa(irs),iN(irs+1)+1)=0 ! Ca_i -- H_{i+1}
|
---|
338 | enddo
|
---|
339 |
|
---|
340 | enddo
|
---|
341 | ! finished initializing matrix conmat
|
---|
342 | ! print *,'Connections matrix'
|
---|
343 | ! Local pair excluded volume
|
---|
344 | do i=1,mxml
|
---|
345 | ilpst(i)=1
|
---|
346 | ilpnd(i)=0
|
---|
347 | enddo
|
---|
348 | do i=1,50*mxrs
|
---|
349 | lcp1(i)=0
|
---|
350 | lcp2(i)=0
|
---|
351 | enddo
|
---|
352 | ilp=0
|
---|
353 | do iml=1,ntlml
|
---|
354 | do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
|
---|
355 | do iat2=iat1+1,iatrs2(irsml2(iml))
|
---|
356 | if ((iat2-iat1.le.mxconr).and.
|
---|
357 | & matcon(iat2-iat1,iat1).eq.2) then
|
---|
358 | ilp=ilp+1
|
---|
359 | lcp1(ilp)=iat1
|
---|
360 | lcp2(ilp)=iat2
|
---|
361 | endif
|
---|
362 | enddo
|
---|
363 | enddo
|
---|
364 |
|
---|
365 | ilpnd(iml)=ilp
|
---|
366 | if (iml.lt.ntlml) then
|
---|
367 | ilpst(iml+1)=ilp+1
|
---|
368 | endif
|
---|
369 | print *,'molecule ',iml,' lc pair range ',ilpst(iml),ilpnd(iml)
|
---|
370 | ! print *,'local pair list'
|
---|
371 | do lci=ilpst(iml),ilpnd(iml)
|
---|
372 | iat1=lcp1(lci)
|
---|
373 | iat2=lcp2(lci)
|
---|
374 | ! print *,lci,iat1,iat2,ityat(iat1),ityat(iat2)
|
---|
375 | enddo
|
---|
376 | enddo
|
---|
377 |
|
---|
378 | print *,'finished initializing Lund force field'
|
---|
379 | end
|
---|
380 |
|
---|
381 | integer function ihptype(iaa)
|
---|
382 | include 'INCL.H'
|
---|
383 | integer iaa, ityp
|
---|
384 |
|
---|
385 | character mynm*4
|
---|
386 | mynm=seq(iaa)
|
---|
387 | ityp=-1
|
---|
388 | call tolost(mynm)
|
---|
389 | if (mynm.eq.'ala') then
|
---|
390 | ityp=1
|
---|
391 | else if ((mynm.eq.'val').or.(mynm.eq.'leu').or.(mynm.eq.'ile')
|
---|
392 | & .or.(mynm.eq.'met').or.(mynm.eq.'pro')) then
|
---|
393 | ityp=2
|
---|
394 | else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp'))
|
---|
395 | & then
|
---|
396 | ityp=3
|
---|
397 | endif
|
---|
398 | ihptype=ityp
|
---|
399 | return
|
---|
400 | end
|
---|
401 |
|
---|
402 | real*8 function ebiasrs(irsd)
|
---|
403 | include 'INCL.H'
|
---|
404 | include 'incl_lund.h'
|
---|
405 | double precision q1, q2, et, xij, yij, zij
|
---|
406 |
|
---|
407 | integer i, iat1, irsd, j, iat2
|
---|
408 |
|
---|
409 | dimension q1(2),q2(2)
|
---|
410 | data q1/-0.2,0.2/
|
---|
411 | data q2/0.42,-0.42/
|
---|
412 |
|
---|
413 | et=0.0
|
---|
414 | do i=0,1
|
---|
415 | iat1=iN(irsd)+i
|
---|
416 | do j=0,1
|
---|
417 | iat2=iC(irsd)+j
|
---|
418 | xij=xat(iat1)-xat(iat2)
|
---|
419 | yij=yat(iat1)-yat(iat2)
|
---|
420 | zij=zat(iat1)-zat(iat2)
|
---|
421 | et=et+q1(i+1)*q2(j+1)/sqrt(xij*xij+yij*yij+zij*zij)
|
---|
422 | enddo
|
---|
423 | enddo
|
---|
424 | ebiasrs=kbias*et
|
---|
425 | return
|
---|
426 | end
|
---|
427 | ! Evaluates backbone backbone hydrogen bond strength for residues
|
---|
428 | ! i and j, taking the donor from residue i and acceptor from residue j
|
---|
429 | real*8 function ehbmmrs(i,j)
|
---|
430 | include 'INCL.H'
|
---|
431 | include 'incl_lund.h'
|
---|
432 | double precision rdon2, racc2, evlu
|
---|
433 |
|
---|
434 | integer i, j
|
---|
435 |
|
---|
436 | double precision r2,r4,r6,dx,dy,dz,ca,cb
|
---|
437 | integer d1,d2,a1,a2
|
---|
438 | d1=iN(i)
|
---|
439 | d2=d1+1
|
---|
440 | a1=iC(j)+1 ! dipoles are numbered from -ve to +ve charge
|
---|
441 | a2=iC(j) ! so, acceptor 1 is the Oxygen, and a2, the carbon
|
---|
442 | dx=xat(a1)-xat(d2)
|
---|
443 | dy=yat(a1)-yat(d2)
|
---|
444 | dz=zat(a1)-zat(d2)
|
---|
445 | r2=dx*dx+dy*dy+dz*dz
|
---|
446 | if (r2.gt.cthb2) then
|
---|
447 | ! print *,'hbmm = 0 ',cthb2,r2,a1,a2,d1,d2
|
---|
448 | ! print *,'a1,a2,d1,d2,r2 = ',a1,a2,d1,d2,r2,sighb2,cthb
|
---|
449 | ehbmmrs=0
|
---|
450 | return
|
---|
451 | endif
|
---|
452 | ca=(xat(d2)-xat(d1))*dx+(yat(d2)-yat(d1))*dy+(zat(d2)-zat(d1))*dz
|
---|
453 | cb=(xat(a2)-xat(a1))*dx+(yat(a2)-yat(a1))*dy+(zat(a2)-zat(a1))*dz
|
---|
454 | if (powa.gt.0.and.ca.le.0) then
|
---|
455 | ! print *,'hbmm, returning 0 because of angle a'
|
---|
456 | ehbmmrs=0
|
---|
457 | return
|
---|
458 | endif
|
---|
459 | if (powb.gt.0.and.cb.le.0) then
|
---|
460 | ! print *,'hbmm, returning 0 because of angle b'
|
---|
461 | ehbmmrs=0
|
---|
462 | return
|
---|
463 | endif
|
---|
464 | r6=sighb2/r2
|
---|
465 | r4=r6*r6
|
---|
466 | r6=r6*r4
|
---|
467 | rdon2=(xat(d2)-xat(d1))**2+(yat(d2)-yat(d1))**2+
|
---|
468 | & (zat(d2)-zat(d1))**2
|
---|
469 | racc2=(xat(a2)-xat(a1))**2+(yat(a2)-yat(a1))**2+
|
---|
470 | & (zat(a2)-zat(a1))**2
|
---|
471 | evlu=((ca*ca/(r2*rdon2))**(0.5*powa))*
|
---|
472 | & ((cb*cb/(r2*racc2))**(0.5*powb))
|
---|
473 | evlu=evlu*(r6*(5*r6-6*r4)+alhb+blhb*r2)
|
---|
474 | ! print *,'found hbmm contribution ',i,j,epshb1,evlu
|
---|
475 | ehbmmrs=epshb1*evlu
|
---|
476 | return
|
---|
477 | end
|
---|
478 | real*8 function enylun(nml)
|
---|
479 | ! nml = 1 .. ntlml. No provision exists to handle out of range values
|
---|
480 | ! for nml inside this function.
|
---|
481 | include 'INCL.H'
|
---|
482 | include 'incl_lund.h'
|
---|
483 | double precision ebiasrs, shbm1, shbm2, etmp, ehbmmrs, ehp, etmp1
|
---|
484 | double precision xij, yij, zij, r2, r6
|
---|
485 |
|
---|
486 | integer istres, nml, indres, i, j, ihpi, ihptype, ihpj, i1, i2
|
---|
487 | integer iat1, iat2, iatt1, iatt2
|
---|
488 |
|
---|
489 | character mynm*4
|
---|
490 | logical prlvr
|
---|
491 | eyhb=0.0 ! backbone-backbone and sidechain-backbone HB
|
---|
492 | eyel=0.0 ! Bias, or local electrostatic term
|
---|
493 | eysl=0.0 ! Hydrophobicity, replaces the solvent term of SMMP
|
---|
494 | eyvr=0.0 ! Local pair excluded volume, in a sense a variable potential
|
---|
495 | eyvw=0.0 ! atom-atom repulsion, excluded volume
|
---|
496 | ! atom-atom repulsion is calculated on a system wide basis, instead of
|
---|
497 | ! molecule by molecule for efficiency. Look into function exvlun.
|
---|
498 |
|
---|
499 | istres=irsml1(nml)
|
---|
500 | indres=irsml2(nml)
|
---|
501 |
|
---|
502 | ! First, all terms that can be calculated on a residue by residue basis
|
---|
503 | do i=istres,indres
|
---|
504 | mynm=seq(i)
|
---|
505 | call tolost(mynm)
|
---|
506 | if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
|
---|
507 | & .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
|
---|
508 | & .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
|
---|
509 | prlvr=.true. ! residue i is a proline variant
|
---|
510 | else
|
---|
511 | prlvr=.false.
|
---|
512 | endif
|
---|
513 |
|
---|
514 | ! Bias, or local electrostatic term. Excluded from the list are
|
---|
515 | ! residues at the ends of the chain, glycine and all proline variants
|
---|
516 | if ((i.ne.istres).and.(i.ne.indres).and.
|
---|
517 | & .not.prlvr.and.mynm.ne.'gly') then
|
---|
518 | eyel=eyel+ebiasrs(i)
|
---|
519 | endif
|
---|
520 | ! Backbone--backbone hydrogen bonds
|
---|
521 | shbm1=1.0
|
---|
522 | shbm2=1.0
|
---|
523 | if ((i.eq.istres).or.(i.eq.indres)) shbm1=0.5
|
---|
524 | ! Residue i contributes the donor, and j, the acceptor, so both i and
|
---|
525 | ! j run over the whole set of amino acids.
|
---|
526 | ! No terms for residue i, if it is a proline variant.
|
---|
527 | if (.not.prlvr) then
|
---|
528 | do j=istres,indres
|
---|
529 | if ((j.lt.(i-2)).or.(j.gt.(i+1))) then
|
---|
530 | shbm2=1.0
|
---|
531 | if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
|
---|
532 | etmp=ehbmmrs(i,j)
|
---|
533 | eyhb=eyhb+shbm1*shbm2*etmp
|
---|
534 | endif
|
---|
535 | enddo
|
---|
536 | endif
|
---|
537 | ! Hydrophobicity, only if residue i is hydrophobic to start with
|
---|
538 | ihpi=ihptype(i)
|
---|
539 | if (ihpi.ge.0) then
|
---|
540 | ! Unlike hydrogen bonds, the hydrophobicity potential is symmetric
|
---|
541 | ! in i and j. So, the loop for j runs from i+1 to the end.
|
---|
542 |
|
---|
543 | do j=i+1,indres
|
---|
544 | ihpj=ihptype(j)
|
---|
545 | if (ihpj.ge.0) then
|
---|
546 | etmp=ehp(i,j,ihpi,ihpj)
|
---|
547 | if (j.eq.(i+1)) etmp=0
|
---|
548 | if (j.eq.(i+2)) etmp=0.5*etmp
|
---|
549 | ! print *, 'hydrophobicity contribution ',etmp,i,j
|
---|
550 | eysl=eysl+etmp
|
---|
551 | endif
|
---|
552 | enddo
|
---|
553 | endif
|
---|
554 | enddo
|
---|
555 |
|
---|
556 | ! Terms that are not calculated residue by residue ...
|
---|
557 |
|
---|
558 | ! Local pair or third-neighbour excluded volume
|
---|
559 | ! Numerically this is normally the term with largest positive
|
---|
560 | ! contribution to the energy in an equilibrated stystem.
|
---|
561 |
|
---|
562 | i1=ilpst(nml)
|
---|
563 | i2=ilpnd(nml)
|
---|
564 | etmp=0.0
|
---|
565 | do i=i1,i2
|
---|
566 | etmp1=0.0
|
---|
567 | iat1=lcp1(i)
|
---|
568 | iat2=lcp2(i)
|
---|
569 | iatt1=ityat(iat1)
|
---|
570 | iatt2=ityat(iat2)
|
---|
571 | xij=xat(iat1)-xat(iat2)
|
---|
572 | yij=yat(iat1)-yat(iat2)
|
---|
573 | zij=zat(iat1)-zat(iat2)
|
---|
574 | r2=(xij*xij + yij*yij + zij*zij)
|
---|
575 | if (r2.le.exvcut2) then
|
---|
576 | r6=sig2lcp(iatt1,iatt2)/r2
|
---|
577 | r6=r6*r6*r6
|
---|
578 | etmp1=(r6*r6+asalcp(iatt1,iatt2)+bsalcp(iatt1,iatt2)*r2)
|
---|
579 | if (etmp1.ge.0) then
|
---|
580 | ! print *,'local pair contribution ',iat1,iat2,iatt1,
|
---|
581 | ! & iatt2,sqrt(sig2lcp(iatt1,iatt2)),etmp1
|
---|
582 | endif
|
---|
583 | etmp=etmp+etmp1
|
---|
584 | endif
|
---|
585 | ! print *,'pair : ',iat1,iat2,' contribution ',etmp1
|
---|
586 | ! print *,exvcut2,r2
|
---|
587 | enddo
|
---|
588 | eyvr=exvk*etmp
|
---|
589 | eysm=eyel+eyhb+eyvr+eysl
|
---|
590 | enylun=eysm
|
---|
591 | end
|
---|
592 | real*8 function eninlun()
|
---|
593 | include 'INCL.H'
|
---|
594 | eysmi=0.0
|
---|
595 | eyeli=0.0
|
---|
596 | eyvwi=0.0
|
---|
597 | eyhbi=0.0
|
---|
598 | eninlun=0.0
|
---|
599 | return
|
---|
600 | end
|
---|
601 | real*8 function ehp(i1,i2,ihp1,ihp2)
|
---|
602 | include 'INCL.H'
|
---|
603 | include 'incl_lund.h'
|
---|
604 | double precision b2, a2, r2min, xij, yij, zij, dtmp, ssum
|
---|
605 |
|
---|
606 | integer ihp1, ihp2, ni, i1, nj, i2, i, k, j, l
|
---|
607 |
|
---|
608 | dimension r2min(12)
|
---|
609 |
|
---|
610 | b2=20.25
|
---|
611 | a2=12.25
|
---|
612 | ! ihp1=ihptype(i1)
|
---|
613 | ! ihp2=ihptype(i2)
|
---|
614 | if ((ihp1.le.0).or.(ihp2.le.0)) then
|
---|
615 | ehp=0.0
|
---|
616 | return
|
---|
617 | endif
|
---|
618 | ni=nhpat(i1)
|
---|
619 | nj=nhpat(i2)
|
---|
620 | do i=1,ni+nj
|
---|
621 | r2min(i)=100000000 ! quite far
|
---|
622 | enddo
|
---|
623 | do i=1,ni
|
---|
624 | k=ihpat(i1,i)
|
---|
625 | do j=1,nj
|
---|
626 | l=ihpat(i2,j)
|
---|
627 | xij=xat(k)-xat(l)
|
---|
628 | yij=yat(k)-yat(l)
|
---|
629 | zij=zat(k)-zat(l)
|
---|
630 | dtmp=(xij*xij + yij*yij + zij*zij)
|
---|
631 | if (dtmp.le.r2min(i)) r2min(i)=dtmp
|
---|
632 | if (dtmp.le.r2min(ni+j)) r2min(ni+j)=dtmp
|
---|
633 | enddo
|
---|
634 | enddo
|
---|
635 | ssum=0
|
---|
636 | do i=1,ni+nj
|
---|
637 | if (r2min(i).le.b2) then
|
---|
638 | if (r2min(i).lt.a2) then
|
---|
639 | ssum=ssum+1
|
---|
640 | else
|
---|
641 | ssum=ssum+(b2-r2min(i))/(b2-a2)
|
---|
642 | endif
|
---|
643 | endif
|
---|
644 | ! hpstrgth=hpstrg((ihp1-1)*nhptyp+ihp2)
|
---|
645 | ! print *, 'hp diagnosis ',ni,nj,ssum/(ni+nj),hpstrgth
|
---|
646 | enddo
|
---|
647 | ehp=-hpstrg((ihp1-1)*nhptyp+ihp2)*ssum/(ni+nj)
|
---|
648 | return
|
---|
649 | end
|
---|
650 |
|
---|
651 | real*8 function exvlun(nml)
|
---|
652 | include 'INCL.H'
|
---|
653 | include 'incl_lund.h'
|
---|
654 | ! For multi-chain systems it makes little sense to split the calculation
|
---|
655 | ! of this term into an 'interaction part' and a contribution from
|
---|
656 | ! individual molecules. So, normally this should always be called with
|
---|
657 | ! argument nml=0. Only for diagnostic reasons, you might want to find
|
---|
658 | ! the contribution from one molecule in a multi-chain system assuming
|
---|
659 | ! there was no other molecule.
|
---|
660 | double precision xmin, ymin, zmin, xmax, ymax, zmax, sizex, sizey
|
---|
661 | double precision sizez, shiftx, shifty, shiftz, etmp, xij, yij
|
---|
662 | double precision zij, r2, r6, etmp1
|
---|
663 |
|
---|
664 | integer nml, istat, indat, i, incell, ndx, ndy, ndz, nxy, ncell
|
---|
665 | integer nocccl, locccl, j, mx, my, mz, icellj, icell, jj, isort
|
---|
666 | integer icl, lcell, ix, iz, iy, nex, ney, nez, nsame, nngbr, jx
|
---|
667 | integer jy, jz, jcl, ii, ngbr, i1, iat1, i2, iat2, iatt1, iatt2
|
---|
668 |
|
---|
669 | dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell)
|
---|
670 | dimension icell(mxat)
|
---|
671 | integer :: jx1, jy1, jz1
|
---|
672 |
|
---|
673 | if (nml.eq.0) then
|
---|
674 | istat=iatrs1(irsml1(1))
|
---|
675 | indat=iatrs2(irsml2(ntlml))
|
---|
676 | else
|
---|
677 | istat=iatrs1(irsml1(nml))
|
---|
678 | indat=iatrs2(irsml2(nml))
|
---|
679 | endif
|
---|
680 |
|
---|
681 | eyvw=0.0
|
---|
682 | ! The beginning part of this implementation is very similar to the
|
---|
683 | ! assignment of cells to the atoms during calculation of solvent
|
---|
684 | ! accessible surface area. So, much of that part is similar. But
|
---|
685 | ! unlike the accessible surface calculations, this term is symmetric
|
---|
686 | ! in any two participating atoms. So, the part after the assignment
|
---|
687 | ! of cells differs even in the structure of the code.
|
---|
688 |
|
---|
689 | do i=1,mxcell
|
---|
690 | incell(i)=0
|
---|
691 | enddo
|
---|
692 | ! print *,'evaluating general excluded volume :',istat,',',indat
|
---|
693 | ! Find minimal containing box
|
---|
694 | xmin=xat(istat)
|
---|
695 | ymin=yat(istat)
|
---|
696 | zmin=zat(istat)
|
---|
697 | xmax=xmin
|
---|
698 | ymax=ymin
|
---|
699 | zmax=zmin
|
---|
700 |
|
---|
701 | do i=istat,indat
|
---|
702 | if (xat(i).le.xmin) then
|
---|
703 | xmin=xat(i)
|
---|
704 | else if (xat(i).ge.xmax) then
|
---|
705 | xmax=xat(i)
|
---|
706 | endif
|
---|
707 | if (yat(i).le.ymin) then
|
---|
708 | ymin=yat(i)
|
---|
709 | else if (yat(i).ge.ymax) then
|
---|
710 | ymax=yat(i)
|
---|
711 | endif
|
---|
712 | if (zat(i).le.zmin) then
|
---|
713 | zmin=zat(i)
|
---|
714 | else if (zat(i).ge.zmax) then
|
---|
715 | zmax=zat(i)
|
---|
716 | endif
|
---|
717 | enddo
|
---|
718 |
|
---|
719 | sizex=xmax-xmin
|
---|
720 | sizey=ymax-ymin
|
---|
721 | sizez=zmax-zmin
|
---|
722 | ! Number of cells along each directions that fit into the box.
|
---|
723 | ndx=int(sizex/exvcutg)+1
|
---|
724 | ndy=int(sizey/exvcutg)+1
|
---|
725 | ndz=int(sizez/exvcutg)+1
|
---|
726 |
|
---|
727 | nxy=ndx*ndy
|
---|
728 | ncell=nxy*ndz
|
---|
729 | ! print *,'Number of cells along x,y,z = ',ndx,',',ndy,',',ndz
|
---|
730 | if (ncell.ge.mxcell) then
|
---|
731 | print *,'exvlun> required number of cells',ncell,
|
---|
732 | & ' exceeded the limit ',mxcell
|
---|
733 | print *,'recompile with a higher mxcell.'
|
---|
734 | stop
|
---|
735 | endif
|
---|
736 | ! Expand box to contain an integral number of cells along each direction
|
---|
737 | shiftx=(dble(ndx)*exvcutg-sizex)/2.0
|
---|
738 | shifty=(dble(ndy)*exvcutg-sizey)/2.0
|
---|
739 | shiftz=(dble(ndz)*exvcutg-sizez)/2.0
|
---|
740 | xmin=xmin-shiftx
|
---|
741 | ymin=ymin-shifty
|
---|
742 | zmin=zmin-shiftz
|
---|
743 | xmax=xmax+shiftx
|
---|
744 | ymax=ymax+shifty
|
---|
745 | zmax=zmax+shiftz
|
---|
746 |
|
---|
747 | ! Set occupied cells to zero. Note that the maximum number of occupied
|
---|
748 | ! cells is the number of atoms in the system.
|
---|
749 | nocccl=0
|
---|
750 | do i=1,mxat
|
---|
751 | locccl(i)=0
|
---|
752 | enddo
|
---|
753 |
|
---|
754 | ! Put atoms in cells
|
---|
755 | do j=istat,indat
|
---|
756 | mx=min(int(max((xat(j)-xmin)/exvcutg,0.0d0)),ndx-1)
|
---|
757 | my=min(int(max((yat(j)-ymin)/exvcutg,0.0d0)),ndy-1)
|
---|
758 | mz=min(int(max((zat(j)-zmin)/exvcutg,0.0d0)),ndz-1)
|
---|
759 | icellj=mx+my*ndx+mz*nxy+1
|
---|
760 | icell(j)=icellj
|
---|
761 | if (icellj.gt.mxcell) then
|
---|
762 | print *,'exvlun> bad cell index ',icellj,' for atom ',j
|
---|
763 | stop
|
---|
764 | else
|
---|
765 | if (incell(icellj).eq.0) then
|
---|
766 | ! previously unoccupied cell
|
---|
767 | nocccl=nocccl+1
|
---|
768 | locccl(nocccl)=icellj
|
---|
769 | endif
|
---|
770 | incell(icellj)=incell(icellj)+1
|
---|
771 | endif
|
---|
772 | enddo
|
---|
773 | ! print *,'finished assigning cells. nocccl = ',nocccl
|
---|
774 | ! Cummulative occupancy of i'th cell
|
---|
775 | do i=1,ncell
|
---|
776 | incell(i+1)=incell(i+1)+incell(i)
|
---|
777 | enddo
|
---|
778 | ! print *,'finished making cumulative cell sums'
|
---|
779 | ! Sorting atoms by their cell index
|
---|
780 | do i=istat,indat
|
---|
781 | j=icell(i)
|
---|
782 | jj=incell(j)
|
---|
783 | isort(jj)=i
|
---|
784 | incell(j)=jj-1
|
---|
785 | enddo
|
---|
786 | ! print *,'sorted atoms by cell index'
|
---|
787 | etmp=0.0
|
---|
788 | do icl=1,nocccl
|
---|
789 | ! loop through occupied cells
|
---|
790 | lcell=locccl(icl)
|
---|
791 | ix=mod(lcell-1,ndx)
|
---|
792 | iz = (lcell - 1) / nxy
|
---|
793 | iy = ((lcell - 1) - iz * nxy) / ndx
|
---|
794 |
|
---|
795 | c print *,'icl=',icl,'absolute index of cell = ',lcell
|
---|
796 | c print *,'iz,iy,ix = ',iz,iy,ix
|
---|
797 | c find all atoms in current cell and all its forward-going neighbours
|
---|
798 | nex=ix+1
|
---|
799 | ney=iy+1
|
---|
800 | nez=iz+1
|
---|
801 | nsame=0
|
---|
802 | nngbr=0
|
---|
803 | do jx=ix,nex
|
---|
804 | if (jx.ge.ndx) then
|
---|
805 | jx1=0
|
---|
806 | else
|
---|
807 | jx1=jx
|
---|
808 | endif
|
---|
809 | do jy=iy,ney
|
---|
810 | if (jy.ge.ndy) then
|
---|
811 | jy1=0
|
---|
812 | else
|
---|
813 | jy1=jy
|
---|
814 | endif
|
---|
815 | do jz=iz,nez
|
---|
816 | if (jz.ge.ndz) then
|
---|
817 | jz1=0
|
---|
818 | else
|
---|
819 | jz1=jz
|
---|
820 | endif
|
---|
821 | jcl=jx1 + ndx*jy1 + nxy*jz1 + 1
|
---|
822 | ! write (logString, *)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
|
---|
823 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
824 | ! count the total number of neighbours
|
---|
825 | nngbr=nngbr+1
|
---|
826 | if (jx.eq.ix.and.jy.eq.iy.and.jz.eq.iz) then
|
---|
827 | ! count how many neighbours are from the same cell
|
---|
828 | nsame=nsame+1
|
---|
829 | endif
|
---|
830 | ngbr(nngbr)=isort(ii)
|
---|
831 | enddo
|
---|
832 | enddo
|
---|
833 | enddo
|
---|
834 | enddo
|
---|
835 | ! A few more cells need to be searched, so that we cover 13 of the 26
|
---|
836 | ! neighbouring cells.
|
---|
837 | ! 1
|
---|
838 | jx=ix+1
|
---|
839 | if (jx.ge.ndx) jx=0
|
---|
840 | jy=iy
|
---|
841 | jz=iz-1
|
---|
842 | if (jz.lt.0) jz=ndz-1
|
---|
843 | jcl=jx+ndx*jy+nxy*jz+1
|
---|
844 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
845 | nngbr=nngbr+1
|
---|
846 | ngbr(nngbr)=isort(ii)
|
---|
847 | enddo
|
---|
848 | ! 2
|
---|
849 | jx=ix
|
---|
850 | jy=iy-1
|
---|
851 | if (jy.lt.0) jy =ndy-1
|
---|
852 | jz=iz+1
|
---|
853 | if (jz.ge.ndz) jz=0
|
---|
854 | jcl=jx+ndx*jy+nxy*jz+1
|
---|
855 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
856 | nngbr=nngbr+1
|
---|
857 | ngbr(nngbr)=isort(ii)
|
---|
858 | enddo
|
---|
859 | ! 3
|
---|
860 | jx=ix-1
|
---|
861 | if (jx.lt.0)jx =ndx-1
|
---|
862 | jy=iy+1
|
---|
863 | if (jy.ge.ndy) jy=0
|
---|
864 | jz=iz
|
---|
865 | jcl=jx+ndx*jy+nxy*jz+1
|
---|
866 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
867 | nngbr=nngbr+1
|
---|
868 | ngbr(nngbr)=isort(ii)
|
---|
869 | enddo
|
---|
870 | ! 4
|
---|
871 | jx=ix+1
|
---|
872 | if (jx.ge.ndx) jx=0
|
---|
873 | jy=iy+1
|
---|
874 | if (jy.ge.ndy) jy=0
|
---|
875 | jz=iz-1
|
---|
876 | if (jz.lt.0) jz=ndz-1
|
---|
877 | jcl=jx+ndx*jy+nxy*jz+1
|
---|
878 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
879 | nngbr=nngbr+1
|
---|
880 | ngbr(nngbr)=isort(ii)
|
---|
881 | enddo
|
---|
882 | ! 5
|
---|
883 | jx=ix+1
|
---|
884 | if (jx.ge.ndx) jx = 0
|
---|
885 | jy=iy-1
|
---|
886 | if (jy.lt.0) jy=ndy-1
|
---|
887 | jz=iz+1
|
---|
888 | if (jz.ge.ndz) jz=0
|
---|
889 | jcl=jx+ndx*jy+nxy*jz+1
|
---|
890 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
891 | nngbr=nngbr+1
|
---|
892 | ngbr(nngbr)=isort(ii)
|
---|
893 | enddo
|
---|
894 | ! 6
|
---|
895 | jx=ix+1
|
---|
896 | if (jx.ge.ndx) jx=0
|
---|
897 | jy=iy-1
|
---|
898 | if (jy.lt.0) jy=ndy-1
|
---|
899 | jz=iz-1
|
---|
900 | if (jz.lt.0) jz=ndy-1
|
---|
901 | jcl=jx+ndx*jy+nxy*jz+1
|
---|
902 | do ii=incell(jcl)+1,incell(jcl+1)
|
---|
903 | nngbr=nngbr+1
|
---|
904 | ngbr(nngbr)=isort(ii)
|
---|
905 | enddo
|
---|
906 |
|
---|
907 | ! print *,'atoms in same cell ',nsame
|
---|
908 | ! print *,'atoms in neighbouring cells ',nngbr
|
---|
909 | do i1=1,nsame
|
---|
910 | ! Over all atoms from the original cell
|
---|
911 | iat1=ngbr(i1)
|
---|
912 | do i2=i1,nngbr
|
---|
913 | ! Over all atoms in the original+neighbouring cells
|
---|
914 | iat2=ngbr(i2)
|
---|
915 | xij=xat(iat1)-xat(iat2)
|
---|
916 | yij=yat(iat1)-yat(iat2)
|
---|
917 | zij=zat(iat1)-zat(iat2)
|
---|
918 | r2=(xij*xij+yij*yij+zij*zij)
|
---|
919 |
|
---|
920 | if (r2.le.exvcutg2) then
|
---|
921 | if (abs(iat2-iat1).gt.mxconr.or.
|
---|
922 | & matcon(iat2-iat1,iat1).eq.1) then
|
---|
923 | iatt1=ityat(iat1)
|
---|
924 | iatt2=ityat(iat2)
|
---|
925 | r6=sig2exv(iatt1,iatt2)/r2
|
---|
926 | r6=r6*r6*r6
|
---|
927 | etmp1=r6*r6+asaexv(iatt1,iatt2)
|
---|
928 | & +bsaexv(iatt1,iatt2)*r2
|
---|
929 | etmp=etmp+etmp1
|
---|
930 | ! if (etmp1.ge.2000) then
|
---|
931 | ! print *,'contribution ',iat1,iat2,etmp1
|
---|
932 | ! call outpdb(1,'EXAMPLES/clash.pdb')
|
---|
933 | ! stop
|
---|
934 | ! endif
|
---|
935 | endif
|
---|
936 | endif
|
---|
937 | enddo
|
---|
938 | enddo
|
---|
939 | enddo
|
---|
940 | ! irs=1
|
---|
941 | ! do iat=iatrs1(irs),iatrs2(irs)
|
---|
942 | ! do j=-mxconr,mxconr
|
---|
943 | ! print *,iat,j,':',matcon(j,iat)
|
---|
944 | ! enddo
|
---|
945 | ! enddo
|
---|
946 | ! irs=irsml2(1)
|
---|
947 | ! do iat=iatrs1(irs),iatrs2(irs)
|
---|
948 | ! do j=-mxconr,mxconr
|
---|
949 | ! print *,iat,j,':',matcon(j,iat)
|
---|
950 | ! enddo
|
---|
951 | ! enddo
|
---|
952 |
|
---|
953 | eyvw=exvk*etmp
|
---|
954 | exvlun=eyvw
|
---|
955 | return
|
---|
956 | end
|
---|
957 |
|
---|
958 | real*8 function exvbrfc()
|
---|
959 | ! Brute force excluded volume evaluation
|
---|
960 | include 'INCL.H'
|
---|
961 | include 'incl_lund.h'
|
---|
962 |
|
---|
963 | double precision etmp, etmp1, xij, yij, zij, r2, r6
|
---|
964 |
|
---|
965 | integer iat1, iat2, iatt1, iatt2
|
---|
966 |
|
---|
967 | etmp=0.0
|
---|
968 | etmp1=0.0
|
---|
969 | print *,'max connection radius is ',mxconr
|
---|
970 | do iat1=iatrs1(irsml1(1)),iatrs2(irsml2(ntlml))
|
---|
971 | do iat2=iat1+1,iatrs2(irsml2(ntlml))
|
---|
972 | xij=xat(iat1)-xat(iat2)
|
---|
973 | yij=yat(iat1)-yat(iat2)
|
---|
974 | zij=zat(iat1)-zat(iat2)
|
---|
975 | r2=(xij*xij+yij*yij+zij*zij)
|
---|
976 |
|
---|
977 | if (r2.le.exvcutg2) then
|
---|
978 | if (abs(iat2-iat1).gt.mxconr.or.
|
---|
979 | & matcon(iat2-iat1,iat1).eq.1) then
|
---|
980 | iatt1=ityat(iat1)
|
---|
981 | iatt2=ityat(iat2)
|
---|
982 | r6=sig2exv(iatt1,iatt2)/r2
|
---|
983 | r6=r6*r6*r6
|
---|
984 | etmp1=r6*r6+asaexv(iatt1,iatt2)
|
---|
985 | & +bsaexv(iatt1,iatt2)*r2
|
---|
986 | etmp=etmp+etmp1
|
---|
987 | else
|
---|
988 | ! print *,'atoms ', iat1,' and ',iat2,' were close',
|
---|
989 | ! # 'but matcon is ',matcon(iat2-iat1,iat1)
|
---|
990 | endif
|
---|
991 | endif
|
---|
992 | enddo
|
---|
993 | enddo
|
---|
994 | exvbrfc=etmp*exvk
|
---|
995 | return
|
---|
996 | end
|
---|