C diff_acc C C Evaluates the 'diffusion accessibility' of the surface C atoms in a protein structure. The diffusion accessibility C is obtained by solving the steady state equations of C diffusion (Poisson's eqn) with the protein acting as an C adsorber and an outer sphere acting as a source. C The extremely large system of linear equations is solved C using a Jacobian iterative method with SOR (successive C overrelaxation) and Chebyshev "checkerboard" acceleration. C The number of cycles required for convergence is roughly C twice the number of grid points in the longest direction. C C The diffusion accessibility of each atom is evaluated C according to the maximum flux vector magnitude in its C vicinity. The values are expressed as a percentage of C the highest value for any atom (actually the max is set C at 99 instead of 100 to give an output that C reads better in the f6.2 field format of the pdb file). C The values are output in the 'temperature factor' column of C a PDB file to allow coloring by temperature factor in a C variety of molecular graphics programs. C C For algorithmic details: C Yeates, T. O. (1995). Evaluating the Long- C Range Accessibility of Protein Surfaces. J. Mol. Biol. C 249, 804-815. C C T. O. Yeates C copyright 1994 C UCLA Dept of Chemistry and Biochemistry C This program may be distributed freely as long as proper C attribution is retained. C C version 2. - Jan., 2004 C version 3. - June 2013 C C C implicit none C integer*4 maxsize, maxatoms parameter (maxsize= 16 000 000) parameter (maxatoms= 360 000) real*4 pot(0:maxsize), output2(maxatoms) integer*4 icode(0:maxsize) real*4 x(3), range(3), rm, rgt, d, cm(3) integer*4 ngr(3), ind(3), jatt(maxatoms) character*60 filein, fileout character resnam*3, atnam*4, junk2*4, chain*2 character*2 refnam(20) real*4 standard_rad(20) character*3 rt(maxatoms) character*30 junk real*4 xch, output, grid, r2, pi, sp, wopt, rmstest, rnormtest real*4 vdwsq, dsq, sq1, sq2, sq3, grad2, acc_max, scale, occ real*4 gradmax, dx, dy, dz, pad, probe_radius real*4 sample_radius, atom_radius, buffer, fudge real*4 diels, dielp, delta, w, bigatom, vol, x0, xx, yy, xtmp real*4 rad_gyr, unitflux, vdw, zz, ytmp, ztmp integer*4 ntypes, nchar, i, ntest, j, k, nsp, ngrtot, ndel, nprot integer*4 iat, iattype, np0, np1, np5, nc0, nc1, nc integer*4 indgrad1, indgrad2, indgrad3, indgrad4, indgrad5 integer*4 indgrad6, indgrad7, indgrad8, ii, jj integer*4 index, indtmp, kk, nat1, nat2, nat3, nat4, natoms integer*4 nres, ncyc, ncyc2 common diels,dielp,ntest,ncyc2,w(1000),delta,rmstest,rnormtest, % x0(3),grid,xch(4,200 000),nat1,nat2,nat3,nat4, % output(200 000), natoms parameter (pad=2.) parameter (probe_radius= 1.4) parameter (sample_radius= 1.5) parameter (atom_radius=1.7) parameter (buffer= 1.5) parameter (fudge=0.5) C C C ***** Files and grid parameters C C write (6,*) ' Input coordinate filename ' read (5,901) filein C write (6,901) filein open (1,file=filein,status='old') C open (2,file='wwwtmp/tmpoutput.pdb',status='unknown',form='formatted') open (2,file='tmpoutput.pdb',status='unknown',form='formatted') C write (6,*) ' Filename for output of diffusion flux ' read (5,901) fileout C write(6,901) fileout open (4,file=fileout,status='new',form='formatted') 901 format (a60) C write (6,*) ' Input grid spacing in Angstroms ' C write (6,*) ' (0.75 suggested [or coarser to increase speed])' C read (5,*) grid C ***** Experience shows that 0.75A is a good grid spacing - the program C automatically adjusts to coarser spacing if required grid=0.75 C write (6,*) grid C *** use fixed atomic radii *** refnam(1)="C " standard_rad(1)=1.9 refnam(2)="N " standard_rad(2)=1.7 refnam(3)="O " standard_rad(3)=1.4 refnam(4)="S " standard_rad(4)=1.9 refnam(5)="P " standard_rad(5)=1.8 refnam(6)="A " standard_rad(6)=1.5 refnam(7)="FE" standard_rad(7)=2.0 refnam(8)="MG" standard_rad(8)=2.0 ntypes = 8 C *** READ ATOMIC RADII FROM A FILE *** write (6,*) '0' nchar=1 C ntypes=1 C open (3,file='radius.dat',status='old') C bigatom=0. C do while (.true.) C read (3,803,end=85) refnam(ntypes), standard_rad(ntypes) C write (6,803) refnam(ntypes), standard_rad(ntypes) C if (standard_rad(ntypes) .gt. bigatom) bigatom=standard_rad(ntypes) C ntypes=ntypes+1 C end do C 85 continue C ntypes=ntypes-1 C 803 format (a2,1x,f8.4) C *** read atoms write (6,*) '4' natoms=0 do while (.true.) 83 read(1,904,err=83,end=10) junk2,atnam,resnam,chain,nres,x C *** currently excludes heteroatoms, such as ligands if (junk2 .ne. 'ATOM') goto 83 natoms=natoms+1 if (natoms .gt. MAXATOMS) stop ' error: too many atoms ' write (2,904) junk2, atnam, resnam, chain, nres, x end do 10 continue 904 format (a4, 9x, a4, a3, a2, i4, 4x, 3f8.3) rewind (2) C C *** read coordinates for center of mass *** C write (6,*) '9' do i=1,3 cm(i)=0. end do do iat=1,natoms read(2,902) junk, x do i=1,3 cm(i)=cm(i)+x(i) end do end do 902 format (a30,3f8.3) do i=1,3 cm(i)=cm(i)/float(natoms) end do rewind(2) C write (6,903) (cm(i),i=1,3) 903 format (/,' CM : ',/,3f10.3) C C *** find largest distance from center, and radius of gyration *** C rm=0. rgt=0. do i=1,natoms d=0. read(2,902) junk, x do j=1,3 d=d+(x(j)-cm(j))**2 end do if (d .gt. rm) rm=d rgt=rgt+d end do rm=sqrt(rm)+probe_radius+bigatom rad_gyr=sqrt(rgt/natoms)+probe_radius+bigatom C write (6,*) ' Molecular radius (max, gyr, extended boundary) =', C % rm, rad_gyr, rm*buffer rm=rm*buffer unitflux=rm/rad_gyr/(rm-rad_gyr) rewind(2) C write (6,*) '12' C C *** get size and origin of grid *** C note that the potential and dielectric grids extend beyond the charge C (and residual) grid. C vol=1. ntest=1 do i=1,3 x0(i)=cm(i)-rm-pad range(i)=2.*(rm+pad) vol=vol*range(i) ngr(i)=2+int(range(i)/grid) ntest=ntest*ngr(i) end do C write(6,*) ' Origin of (potential) map at ', x0 if (ntest .gt. maxsize) then grid=(vol/maxsize)**0.3333 do while (.true.) ntest=1 do i=1,3 ngr(i)=2+int(range(i)/grid) ntest=ntest*ngr(i) end do if (ntest .le. maxsize) goto 20 grid=grid*1.02 end do 20 continue C write (6,*) ' Grid size reset to ', grid end if C write (6,*) ' Grid size = ', (ngr(i),i=1,3) C ngrtot=ngr(1)*ngr(2)*ngr(3) C write (6,*) '18' C C *** spectral radius *** C pi=4.*atan(1.) sp = cos(pi/(ngr(1)*fudge)) wopt=2./(1.+sqrt(1.-sp*sp)) C write(6,*) ' Spectral radius is approximately ',sp C write(6,*) ' Optimal (asymtotic) over-relaxation factor ',wopt C write(6,*) ' The relaxation factor will be varied according to ' C write(6,*) ' the Chebyshev acceleration for iteration on an ' C write(6,*) ' alternating checkerboard.' C write(6,*) ' ' C C *** # of cycles *** C C write(6,*) ' ' C write(6,*) ' Number of cycles (half the number of half-iterations) ' C write(6,*) ' Number of cycles for iteration ' ncyc=ngr(1) C write(6,*) ncyc C ncyc2=2*ncyc ncyc2=ncyc C C *** omega *** C write (6,*) '23' w(1)=1. w(2)=1./(1.-sp*sp/2.) do j=3,ncyc2 w(j)=1./(1.-w(j-1)*sp*sp/4.) end do C write(6,*) ' 1st 20 SOR factors ' C write (6,101) (w(i),i=1,20) C write(6,*) ' Last 20 SOR factors ' C write (6,101) (w(i),i=ncyc2-19,ncyc2) 101 Format (10f7.4) C C *** termination conditions *** C write (6,*) '29' C type *,' ' C type *,' RMS and NORM (of residual) for termination ' C accept *, rmstest,rnormtest C write(6,*) rmstest,rnormtest C type *,' ' rmstest=0. rnormtest=0. C C *** set up boundary conditions for potential *** C write (6,*) '34' do i=0,ngrtot-1 pot(i)=1. icode(i)=-1 end do r2=rm**2 nsp=0 do k=0,ngr(3)-1 do j=0,ngr(2)-1 do i=0,ngr(1)-1 zz=x0(3)+k*grid-cm(3) yy=x0(2)+j*grid-cm(2) xx=x0(1)+i*grid-cm(1) index=i+j*ngr(1)+k*ngr(1)*ngr(2) if ((xx*xx+yy*yy+zz*zz) .lt. r2) then pot(index)=0.5 icode(index)=0 nsp=nsp+1 end if end do end do end do C write (6,*) ' # pts inside spherical boundary = ', nsp 992 format(I30) C C *** reread coordinates for xyz and radii *** C write (6,*) '44' ndel=int(vdw/grid+0.5) vdwsq=vdw**2 nprot=0 do iat=1,natoms read(2,904) junk2, atnam, resnam, chain, nres, % (xch(i,iat),i=1,3) do i=1,3 x(i)=xch(i,iat) end do iattype=0 do i=1,ntypes if (refnam(i)(1:nchar) .eq. atnam(1:nchar)) then iattype=i jatt(iat)=i goto 87 end if end do 87 continue if (iattype .eq. 0) then C write (6,804) atnam iattype=1 C stop end if xch(4,iat)=standard_rad(iattype) rt(iat)=resnam 804 format (' WARNING : unknown atom type ', a4, ' -- using default') vdw=xch(4,iat)+probe_radius ndel=int(vdw/grid+0.5) vdwsq=vdw**2 C write (6,*) vdw, ndel, vdwsq do i=1,3 ind(i)=nint((x(i)-x0(i))/grid) end do do kk=-ndel,ndel sq3=(x(3)-x0(3)-(ind(3)+kk)*grid)**2 do jj=-ndel,ndel sq2=(x(2)-x0(2)-(ind(2)+jj)*grid)**2 do ii=-ndel,ndel sq1=(x(1)-x0(1)-(ind(1)+ii)*grid)**2 dsq=sq1+sq2+sq3 if (dsq .lt. vdwsq) then indtmp = ind(1)+ii + ngr(1)*(ind(2)+jj) + % ngr(1)*ngr(2)*(ind(3)+kk) nprot=nprot+1 icode(indtmp)=iat pot(indtmp)=0. end if end do end do end do output(iat)=0. end do C write (6,*) ' # of points inside atoms = ',nprot C write (6,*) '51' do i=0,ngrtot-1 if (pot(i) .eq. 0.) np0=np0+1 if (pot(i) .eq. 1.) np1=np1+1 if (pot(i) .eq. 0.5) np5=np5+1 if (icode(i) .eq. 0) nc0=nc0+1 if (icode(i) .eq. -1) nc1=nc1+1 if (icode(i) .ge. 1) nc=nc+1 end do C write (6,*) ' npprot, npext, npfree, ncfree, ncext, ncatoms ' C write (6,*) np0, np1, np5, nc0, nc1, nc call solve (pot,icode,ngr(1),ngr(2),ngr(3)) do iat=1,natoms gradmax=-999. do i=1,3 x(i)=xch(i,iat) end do vdw=xch(4,iat)+probe_radius+sample_radius ndel=int(vdw/grid+0.5) vdwsq=vdw**2 do i=1,3 ind(i)=nint((x(i)-x0(i))/grid) end do do kk=-ndel-1,ndel sq3=(x(3)-x0(3)-(ind(3)+kk+0.5)*grid)**2 do jj=-ndel-1,ndel sq2=(x(2)-x0(2)-(ind(2)+jj+0.5)*grid)**2 do ii=-ndel-1,ndel sq1=(x(1)-x0(1)-(ind(1)+ii+0.5)*grid)**2 dsq=sq1+sq2+sq3 if (dsq .lt. vdwsq) then indgrad1 = ind(1)+ii + ngr(1)*(ind(2)+jj) + % ngr(1)*ngr(2)*(ind(3)+kk) indgrad2 = ind(1)+ii+1 + ngr(1)*(ind(2)+jj) + % ngr(1)*ngr(2)*(ind(3)+kk) indgrad3 = ind(1)+ii + ngr(1)*(ind(2)+jj+1) + % ngr(1)*ngr(2)*(ind(3)+kk) indgrad4 = ind(1)+ii+1 + ngr(1)*(ind(2)+jj+1) + % ngr(1)*ngr(2)*(ind(3)+kk) indgrad5 = ind(1)+ii + ngr(1)*(ind(2)+jj) + % ngr(1)*ngr(2)*(ind(3)+kk+1) indgrad6 = ind(1)+ii+1 + ngr(1)*(ind(2)+jj) + % ngr(1)*ngr(2)*(ind(3)+kk+1) indgrad7 = ind(1)+ii + ngr(1)*(ind(2)+jj+1) + % ngr(1)*ngr(2)*(ind(3)+kk+1) indgrad8 = ind(1)+ii+1 + ngr(1)*(ind(2)+jj+1) + % ngr(1)*ngr(2)*(ind(3)+kk+1) dx=0.25*(pot(indgrad2)+pot(indgrad4)+ % pot(indgrad6)+pot(indgrad8)- % (pot(indgrad1)+pot(indgrad3)+ % pot(indgrad5)+pot(indgrad7)))/grid dy=0.25*(pot(indgrad3)+pot(indgrad4)+ % pot(indgrad7)+pot(indgrad8)- % (pot(indgrad1)+pot(indgrad2)+ % pot(indgrad5)+pot(indgrad6)))/grid dz=0.25*(pot(indgrad5)+pot(indgrad6)+ % pot(indgrad7)+pot(indgrad8)- % (pot(indgrad1)+pot(indgrad2)+ % pot(indgrad3)+pot(indgrad4)))/grid grad2=(dx*dx+dy*dy+dz*dz) if (grad2 .gt. gradmax) gradmax=grad2 end if end do end do end do output2(iat)=sqrt(gradmax) C end do C set scale for outputing diffusion accessibility acc_max=0. do iat=1,natoms if (output2(iat) .gt. acc_max) acc_max=output2(iat) end do scale=99./acc_max do iat=1,natoms output2(iat)=output2(iat)*scale*1.40 if (output2(iat) .gt. 99.) output2(iat)=99. C tmp=100.*output2(iat)/unitflux C if (tmp .lt. 25.) tmp=25. C if (tmp .gt. 80.) tmp=80. C output2(iat)=tmp end do rewind (2) do iat=1,natoms read (2,902) junk, xtmp, ytmp, ztmp occ=1. write (4,802) junk, xtmp, ytmp, ztmp, occ, output2(iat) end do 802 format (a30, 3f8.3, 2f6.2) 701 format (i6, 2f12.4,2x,a3, i3) C write (6,*) '100' stop end C C ********************************************************************* C subroutine solve (pot,icode,ng1,ng2,ng3) C common diels,dielp,ntest,ncyc2,w(1000),delta,rmstest,rnormtest, % x0(3),grid,xch(4,200 000),nat1,nat2,nat3,nat4, % output(200 000), natoms real*4 pot(0:ng1-1,0:ng2-1,0:ng3-1) real*4 xch, output integer*4 icode(0:ng1-1,0:ng2-1,0:ng3-1) C integer*4 ind(3) C C ********* CYCLE ******************** C tot=ntest C type *,' Cycle #, SOR factor, RMS, Norm :' C write(6,*) ' Half-cycle #, SOR factor:' C write(6,*) ' Iteration cycle:' C write(6,*) ' ' do icyc=1,ncyc2 C sum1=0. rnorm=0. sor2=w(icyc) do k=1,ng3-2 do j=1,ng2-2 ipar=mod(k+j+icyc,2)+1 do i = ipar,ng1-2,2 if (icode(i,j,k) .eq. 0) then C bc=-6.*pot(i,j,k)+pot(i+1,j,k)+pot(i-1,j,k) * +pot(i,j+1,k)+pot(i,j-1,k) * +pot(i,j,k+1)+pot(i,j,k-1) C C ***** Calculate Residual (b-Ax) ***** C resid=-bc rnorm=rnorm+abs(resid) sum1=sum1+resid*resid C C ***** Iterate Potential ***** C pot(i,j,k)=pot(i,j,k)-sor2*resid/6. end if end do end do end do root=sqrt(2.*sum1/tot) C type *, icyc-1,'/2 ',w(icyc),root,rnorm if (10*(icyc/10) .eq. icyc) then C write(6,*) ' ',icyc-1,' ',w(icyc) C write(6,*) ' ',icyc-1 end if if (root .le. rmstest .or. (rnorm) .le. rnormtest) * goto 30 C end do C 30 Continue C write(6,*) ' ' C write(6,*) ' END OF ITERATION ' C write(6,*) ' ' do i=1,natoms output(i)=0. end do C do k=1,ng3-2 do j=1,ng2-2 do i = 1,ng1-2 iat=icode(i,j,k) if (iat .gt. 0) then C bc=-6.*pot(i,j,k)+pot(i+1,j,k)+pot(i-1,j,k) * +pot(i,j+1,k)+pot(i,j-1,k) * +pot(i,j,k+1)+pot(i,j,k-1) C C ***** Accumulate Flux ***** C output(iat)=output(iat)+bc C end if end do end do end do C relic code related to P-B electrostatics C C ***** OUTPUT ************** C Cty write(2) ng1,ng2,ng3, x0, grid Cty write(3) ng1-2,ng2-2,ng3-2, (x0(i)+grid,i=1,3), grid Cty do k=0,ng3-1 Cty write(2) ((pot(i,j,k),i=0,ng1-1),j=0,ng2-1) Cty end do Cty do k=1,ng3-2 Cty write(3) ((diel(i,j,k),i=1,ng1-2),j=1,ng2-2) Cty end do CtyC CtyC ******* CALCULATE FINAL RESIDUAL and 'self energy' ***** CtyC Cty sum1=0. Cty rnorm=0. Cty self=0. Cty do k=1,ng3-2 Cty do j=1,ng2-2 Cty do i = 1,ng1-2 Cty e=diel(i,j,k) Cty e1=diel(i+1,j,k)+e Cty e2=diel(i-1,j,k)+e Cty e3=diel(i,j+1,k)+e Cty e4=diel(i,j-1,k)+e Cty e5=diel(i,j,k+1)+e Cty e6=diel(i,j,k-1)+e Cty et=e1+e2+e3+e4+e5+e6 Cty if (e .ne. dielp) et=et+delta CtyC CtyC ***** Compute Calculated Charge ***** CtyC Cty bc=(et*pot(i,j,k)-e1*pot(i+1,j,k)-e2*pot(i-1,j,k) Cty * -e3*pot(i,j+1,k)-e4*pot(i,j-1,k) Cty * -e5*pot(i,j,k+1)-e6*pot(i,j,k-1))*0.5 CtyC CtyC ***** Calculate Residual (b-Ax) ***** CtyC Cty resid=b(i,j,k)-bc Cty rnorm=rnorm+abs(resid) Cty sum1=sum1+resid*resid CtyC self=self+b(i,j,k)*pot(i,j,k) CtyC Cty end do Cty end do Cty end do Cty do iat=nat3,nat4 Cty charge=xch(4,iat) Cty if (charge .ne. 0.) then Cty do i=1,3 Cty ind(i)=int((xch(i,iat)-x0(i))/grid) Cty end do Cty do kk=0,1 Cty fac3=xch(3,iat)-x0(3)-grid*(ind(3)+kk) Cty do jj=0,1 Cty fac2=xch(2,iat)-x0(2)-grid*(ind(2)+jj) Cty do ii=0,1 Cty fac1=xch(1,iat)-x0(1)-grid*(ind(1)+ii) Cty factor=(grid-abs(fac1))*(grid-abs(fac2))*(grid-abs(fac3)) Cty % /(grid**3) Cty self=self+factor*charge*pot(ind(1)+ii,ind(2)+jj,ind(3)+kk) Cty end do Cty end do Cty end do Cty end if Cty end do Cty Cty root=sqrt(sum1/tot) Cty type *,' RMS following cycle ',ncyc2/2,' : ',root Cty type *,' NORM following cycle ',ncyc2/2,' :',rnorm Cty type *,' Self energy = ', self*4.*3.14159*332./grid CtyC return end