1 | ! **************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the subroutines: enysol,tessel
|
---|
4 | !
|
---|
5 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
6 | ! Shura Hayryan, Chin-Ku
|
---|
7 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
8 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
9 | !
|
---|
10 | ! **************************************************************
|
---|
11 |
|
---|
12 |
|
---|
13 | real*8 function enysol(nmol)
|
---|
14 |
|
---|
15 | include 'INCL.H'
|
---|
16 | include 'mpif.h'
|
---|
17 |
|
---|
18 | ! --------------------------------------------------------------
|
---|
19 | !
|
---|
20 | ! Double Cubic Lattice algorithm for calculating the
|
---|
21 | ! solvation energy of proteins using
|
---|
22 | ! solvent accessible area method.
|
---|
23 | !
|
---|
24 | ! if nmol == 0 do solvation energy over all residues.
|
---|
25 | ! CALLS: nursat
|
---|
26 | !
|
---|
27 | ! -------------------------------------------------------------
|
---|
28 | ! TODO: Check the solvent energy for multiple molecules
|
---|
29 | double precision startwtime, avr_x, avr_y, avr_z, xmin, ymin, zmin
|
---|
30 | double precision xmax, ymax, zmax, rmax, diamax, sizex, sizey
|
---|
31 | double precision sizez, shiftx, shifty, shiftz, boxpp, trad, dx
|
---|
32 | double precision dy, dz, dd, akrad, dr, xyz, radb, radb2, sdd, sdr
|
---|
33 | double precision area, volume, eyslsum, endwtime
|
---|
34 |
|
---|
35 | integer nmol, nrslow, nrshi, nlow, nup, i, numat, inbox, j, ndx
|
---|
36 | integer ndy, ndz, nqxy, ncbox, mx, my, mz, nboxj, numbox, jj
|
---|
37 | integer indsort, iboxmin, iboxmax, ibox, iz, iy, ix, lbn, nsx, nsy
|
---|
38 | integer nsz, nex, ney, nez, jcnt, jz, jy, jx, jbox, ii, look, ia
|
---|
39 | integer jbi, nnei, ib, jtk, il, lst, ilk, ik, icount, jres, nursat
|
---|
40 | integer ierror, nhx, mhx, nbt, mbt
|
---|
41 |
|
---|
42 | dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat)
|
---|
43 | dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox)
|
---|
44 | logical surfc(mxpoint)
|
---|
45 | integer root
|
---|
46 |
|
---|
47 | ! common/ressurf/surfres(mxrs)
|
---|
48 | real*8 tsurfres(mxrs)
|
---|
49 | startwtime = MPI_Wtime()
|
---|
50 | root = 0
|
---|
51 | enysolct = enysolct + 1
|
---|
52 | eyslh = 0.0
|
---|
53 | eyslp = 0.0
|
---|
54 | if (nmol.eq.0) then
|
---|
55 | nrslow=irsml1(1)
|
---|
56 | nrshi=irsml2(ntlml)
|
---|
57 | else
|
---|
58 | nrslow = irsml1(nmol)
|
---|
59 | nrshi = irsml2(nmol)
|
---|
60 | endif
|
---|
61 | nlow = iatrs1(nrslow)
|
---|
62 | nup = iatrs2(nrshi)
|
---|
63 | do i=nrslow,nrshi
|
---|
64 | surfres(i) = 0.0d0
|
---|
65 | end do
|
---|
66 |
|
---|
67 | numat= nup - nlow + 1
|
---|
68 |
|
---|
69 | do i=1,mxbox+1
|
---|
70 | inbox(i)=0
|
---|
71 | end do
|
---|
72 |
|
---|
73 | asa=0.0d0
|
---|
74 | vdvol=0.0d0
|
---|
75 | eysl=0.0d0
|
---|
76 |
|
---|
77 | avr_x=xat(nlow)
|
---|
78 | avr_y=yat(nlow)
|
---|
79 | avr_z=zat(nlow)
|
---|
80 | xmin=xat(nlow)
|
---|
81 | ymin=yat(nlow)
|
---|
82 | zmin=zat(nlow)
|
---|
83 | xmax=xmin
|
---|
84 | ymax=ymin
|
---|
85 | zmax=zmin
|
---|
86 |
|
---|
87 | rmax=rvdw(nlow)
|
---|
88 |
|
---|
89 | do j=nlow+1,nup
|
---|
90 | if(xat(j).le.xmin) then
|
---|
91 | xmin=xat(j)
|
---|
92 | else if(xat(j).ge.xmax) then
|
---|
93 | xmax=xat(j)
|
---|
94 | end if
|
---|
95 | avr_x=avr_x+xat(j)
|
---|
96 | if(yat(j).le.ymin) then
|
---|
97 | ymin=yat(j)
|
---|
98 | else if(yat(j).ge.ymax) then
|
---|
99 | ymax=yat(j)
|
---|
100 | end if
|
---|
101 | avr_y=avr_y+yat(j)
|
---|
102 | if(zat(j).le.zmin) then
|
---|
103 | zmin=zat(j)
|
---|
104 | else if(zat(j).ge.zmax) then
|
---|
105 | zmax=zat(j)
|
---|
106 | end if
|
---|
107 | avr_z=avr_z+zat(j)
|
---|
108 | if(rvdw(j).ge.rmax) rmax=rvdw(j)
|
---|
109 | end do
|
---|
110 |
|
---|
111 | avr_x=avr_x/dble(numat)
|
---|
112 | avr_y=avr_y/dble(numat)
|
---|
113 | avr_z=avr_z/dble(numat)
|
---|
114 | diamax=2.d0*rmax
|
---|
115 |
|
---|
116 | ! The sizes of the big box
|
---|
117 |
|
---|
118 | sizex=xmax-xmin
|
---|
119 | sizey=ymax-ymin
|
---|
120 | sizez=zmax-zmin
|
---|
121 |
|
---|
122 | ! How many maximal diameters in each size ?
|
---|
123 |
|
---|
124 | ndx=sizex/diamax + 1
|
---|
125 | ndy=sizey/diamax + 1
|
---|
126 | ndz=sizez/diamax + 1
|
---|
127 |
|
---|
128 | ! We may need the number of quadratic boxes in (X,Y) plane
|
---|
129 |
|
---|
130 | nqxy=ndx*ndy
|
---|
131 |
|
---|
132 | ! The number of cubic boxes of the size "diamax"
|
---|
133 |
|
---|
134 | ncbox=nqxy*ndz
|
---|
135 | if(ncbox.ge.mxbox) then
|
---|
136 | print *,'enysol> bad ncbox',ncbox
|
---|
137 | stop
|
---|
138 | end if
|
---|
139 |
|
---|
140 | ! Let us shift the borders to home all boxes
|
---|
141 |
|
---|
142 | shiftx=(dble(ndx)*diamax-sizex)/2.d0
|
---|
143 | shifty=(dble(ndy)*diamax-sizey)/2.d0
|
---|
144 | shiftz=(dble(ndz)*diamax-sizez)/2.d0
|
---|
145 | xmin=xmin-shiftx
|
---|
146 | ymin=ymin-shifty
|
---|
147 | zmin=zmin-shiftz
|
---|
148 | xmax=xmax+shiftx
|
---|
149 | ymax=ymax+shifty
|
---|
150 | zmax=zmax+shiftz
|
---|
151 |
|
---|
152 | ! Finding the box of each atom
|
---|
153 |
|
---|
154 | do j=nlow,nup
|
---|
155 | mx=min(int(max((xat(j)-xmin)/diamax,0.0d0)),ndx-1)
|
---|
156 | my=min(int(max((yat(j)-ymin)/diamax,0.0d0)),ndy-1)
|
---|
157 | mz=min(int(max((zat(j)-zmin)/diamax,0.0d0)),ndz-1)
|
---|
158 | nboxj=mx+my*ndx+mz*nqxy+1
|
---|
159 | numbox(j)=nboxj
|
---|
160 | if (nboxj.gt.mxbox) then
|
---|
161 | write (logString, '(a)') 'enysol> bad mxboxe-2'
|
---|
162 | write (logString, *) 'diagnostics ...'
|
---|
163 | write (logString, *) 'atom ',j
|
---|
164 | write (logString, *) 'position ',xat(j),yat(j),zat(j)
|
---|
165 | write (logString, *) 'box indices ',mx,my,mz
|
---|
166 | write (logString, *) 'resulting boxindex and limit ',nboxj,
|
---|
167 | & mxbox
|
---|
168 |
|
---|
169 | stop
|
---|
170 | else
|
---|
171 | inbox(nboxj)=inbox(nboxj)+1
|
---|
172 | end if
|
---|
173 | end do
|
---|
174 |
|
---|
175 | ! Summation over the boxes
|
---|
176 |
|
---|
177 | do i=1,ncbox
|
---|
178 | inbox(i+1)=inbox(i+1)+inbox(i)
|
---|
179 | end do
|
---|
180 |
|
---|
181 |
|
---|
182 | ! Sorting the atoms by the their box numbers
|
---|
183 |
|
---|
184 | do i=nlow,nup
|
---|
185 | j=numbox(i)
|
---|
186 | jj=inbox(j)
|
---|
187 | indsort(jj)=i
|
---|
188 | inbox(j)=jj-1
|
---|
189 | end do
|
---|
190 |
|
---|
191 | ! Getting started
|
---|
192 | ! We have to loop over ncbox boxes and have no processors available
|
---|
193 | boxpp = 1.0 * ncbox / no
|
---|
194 | iboxmin = boxpp * myrank
|
---|
195 | iboxmax = boxpp * (myrank + 1) - 1
|
---|
196 | if (myrank.eq.(no - 1)) iboxmax = ncbox
|
---|
197 | if (myrank.eq.-1) then
|
---|
198 | write (logString, *) 'enysol> Summary:', enysolct, ncbox,
|
---|
199 | & boxpp, ndx, ndy, ndz
|
---|
200 | endif
|
---|
201 | do ibox = iboxmin + 1, iboxmax + 1
|
---|
202 | iz = (ibox - 1) / nqxy
|
---|
203 | iy = (ibox - 1 - iz * nqxy) / ndx
|
---|
204 | ix = ibox - 1 - iy*ndx - iz*nqxy
|
---|
205 | ! ibox=ix+iy*ndx+iz*nqxy + 1
|
---|
206 |
|
---|
207 | ! Does this box contain atoms ?
|
---|
208 |
|
---|
209 | lbn=inbox(ibox+1)-inbox(ibox)
|
---|
210 |
|
---|
211 | if(lbn.gt.0) then ! There are some atoms
|
---|
212 | nsx=max(ix-1,0)
|
---|
213 | nsy=max(iy-1,0)
|
---|
214 | nsz=max(iz-1,0)
|
---|
215 | nex=min(ix+1,ndx-1)
|
---|
216 | ney=min(iy+1,ndy-1)
|
---|
217 | nez=min(iz+1,ndz-1)
|
---|
218 |
|
---|
219 | ! Atoms in the boxes around
|
---|
220 |
|
---|
221 | jcnt=1
|
---|
222 | do jz=nsz,nez
|
---|
223 | do jy=nsy,ney
|
---|
224 | do jx=nsx,nex
|
---|
225 | jbox=jx+jy*ndx+jz*nqxy+1
|
---|
226 | do ii=inbox(jbox)+1,inbox(jbox+1)
|
---|
227 | if(rvdw(indsort(ii)).gt.0.0d0) then
|
---|
228 | look(jcnt)=indsort(ii)
|
---|
229 | jcnt=jcnt+1
|
---|
230 | end if
|
---|
231 | end do
|
---|
232 | end do
|
---|
233 | end do
|
---|
234 | end do
|
---|
235 |
|
---|
236 | do ia=inbox(ibox)+1,inbox(ibox+1)
|
---|
237 | jbi=indsort(ia)
|
---|
238 | trad=rvdw(jbi)
|
---|
239 | if(trad.gt.0.0) then
|
---|
240 | nnei=0
|
---|
241 | do ib=1,jcnt-1
|
---|
242 | jtk=look(ib)
|
---|
243 | if(jtk.ne.jbi)then
|
---|
244 | dx=(xat(jbi)-xat(jtk))/trad
|
---|
245 | dy=(yat(jbi)-yat(jtk))/trad
|
---|
246 | dz=(zat(jbi)-zat(jtk))/trad
|
---|
247 | dd=dx*dx+dy*dy+dz*dz
|
---|
248 | akrad=rvdw(jtk)/trad
|
---|
249 | dr=1.0d0+akrad
|
---|
250 | dr=dr*dr
|
---|
251 | !c if contact
|
---|
252 | if(dd.le.dr) then
|
---|
253 | nnei=nnei+1
|
---|
254 | xyz(nnei,1)=dx
|
---|
255 | xyz(nnei,2)=dy
|
---|
256 | xyz(nnei,3)=dz
|
---|
257 | radb(nnei)=akrad
|
---|
258 | radb2(nnei)=akrad*akrad
|
---|
259 | end if
|
---|
260 | end if
|
---|
261 | end do
|
---|
262 | !c
|
---|
263 | do il=1,npnt
|
---|
264 | surfc(il)=.false.
|
---|
265 | end do
|
---|
266 |
|
---|
267 | ! Check overlap
|
---|
268 |
|
---|
269 | lst=1
|
---|
270 | do il=1,npnt
|
---|
271 | sdd = 0.0d0
|
---|
272 | do ilk=1,3
|
---|
273 | sdd = sdd +(xyz(lst,ilk)+spoint(il,ilk))**2
|
---|
274 | end do
|
---|
275 | if(sdd.gt.radb2(lst)) then
|
---|
276 | do ik=1,nnei
|
---|
277 | sdd =0.0d0
|
---|
278 | do ilk=1,3
|
---|
279 | sdd = sdd +(xyz(ik,ilk)+spoint(il,ilk))**2
|
---|
280 | end do
|
---|
281 | if(sdd.le.radb2(ik)) then
|
---|
282 | lst=ik
|
---|
283 | go to 99
|
---|
284 | end if
|
---|
285 | end do
|
---|
286 | 99 continue
|
---|
287 |
|
---|
288 | if(ik.gt.nnei)then
|
---|
289 | surfc(il)=.true.
|
---|
290 | end if
|
---|
291 | end if
|
---|
292 | end do
|
---|
293 |
|
---|
294 | icount=0
|
---|
295 | dx=0.0d0
|
---|
296 | dy=0.0d0
|
---|
297 | dz=0.0d0
|
---|
298 | do il=1,npnt
|
---|
299 | if(surfc(il)) then
|
---|
300 | icount=icount+1
|
---|
301 | dx=dx+spoint(il,1)
|
---|
302 | dy=dy+spoint(il,2)
|
---|
303 | dz=dz+spoint(il,3)
|
---|
304 | end if
|
---|
305 | end do
|
---|
306 | sdr=4.d0*pi*trad*trad/dble(npnt)
|
---|
307 | area = sdr*dble(icount)
|
---|
308 | volume = sdr/3.0d0*(trad*dble(icount)
|
---|
309 | & +(xat(jbi)-avr_x)*dx
|
---|
310 | & +(yat(jbi)-avr_y)*dy
|
---|
311 | & +(zat(jbi)-avr_z)*dz)
|
---|
312 |
|
---|
313 | asa=asa+area
|
---|
314 | vdvol=vdvol+volume
|
---|
315 | eysl=eysl+area*sigma(jbi)
|
---|
316 | ! Separate hydrophilic (h) and hyrdophobic (p) contributions to eysl
|
---|
317 | if (sigma(jbi).lt.0) then
|
---|
318 | eyslp = eyslp + area * sigma(jbi)
|
---|
319 | asap = asap + area
|
---|
320 | endif
|
---|
321 | if (sigma(jbi).gt.0) then
|
---|
322 | eyslh = eyslh + area * sigma(jbi)
|
---|
323 | asah = asah + area
|
---|
324 | endif
|
---|
325 | ! Measure how much a residue is solvent accessible:
|
---|
326 | jres = nursat(jbi)
|
---|
327 | surfres(jres) = surfres(jres) + area
|
---|
328 | end if
|
---|
329 | end do
|
---|
330 | end if
|
---|
331 | ! end do
|
---|
332 | ! end do
|
---|
333 | end do
|
---|
334 | call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION,
|
---|
335 | & MPI_SUM,my_mpi_comm, ierror)
|
---|
336 | write (logString, *) 'enysol>', myrank, eysl, eyslsum
|
---|
337 | tsurfres = surfres
|
---|
338 | call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION,
|
---|
339 | & MPI_SUM,my_mpi_comm, ierror)
|
---|
340 | eysl = eyslsum
|
---|
341 |
|
---|
342 | endwtime = MPI_Wtime()
|
---|
343 | if (myrank.le.-1) then
|
---|
344 | write (logString, *) 'enysol>',myrank,enysolct,
|
---|
345 | & iboxmin + 1, iboxmax + 1,
|
---|
346 | & endwtime - startwtime, "s"
|
---|
347 | endif
|
---|
348 | !
|
---|
349 | if (isolscl) then
|
---|
350 | nhx=0
|
---|
351 | mhx=0
|
---|
352 | nbt=0
|
---|
353 | mbt=0
|
---|
354 | call helix(nhx,mhx,nbt,mbt)
|
---|
355 | eysl=((nhx+nbt)*eysl)/(irsml2(ntlml)-irsml1(1))
|
---|
356 | endif
|
---|
357 |
|
---|
358 | enysol = eysl
|
---|
359 |
|
---|
360 | return
|
---|
361 | end
|
---|
362 | ! *********************
|
---|
363 | subroutine tessel
|
---|
364 | include 'INCL.H'
|
---|
365 | integer i
|
---|
366 |
|
---|
367 | character lin*80
|
---|
368 |
|
---|
369 | ! Skipping comment lines, which begin with '!'
|
---|
370 |
|
---|
371 | read(20,'(a)') lin
|
---|
372 | do while(lin(1:1).eq.'!')
|
---|
373 | read (20,'(a)') lin
|
---|
374 | end do
|
---|
375 |
|
---|
376 | ! The first non-comment line is the number of the surface points
|
---|
377 |
|
---|
378 | read(lin(1:5),'(i5)') npnt
|
---|
379 | ! write (logString, '(a,i5)') 'the number of points---->',npnt
|
---|
380 |
|
---|
381 | ! Read the surface points
|
---|
382 |
|
---|
383 | do i=1,npnt
|
---|
384 | read(20,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
|
---|
385 |
|
---|
386 | ! write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
|
---|
387 | end do
|
---|
388 |
|
---|
389 | return
|
---|
390 |
|
---|
391 | end
|
---|
392 |
|
---|