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