source/unres/src_CSA/readpdb.F

   1       subroutine readpdb
   2 C Read the PDB file and convert the peptide geometry into virtual-chain
   3 C geometry.
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6       include 'COMMON.LOCAL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.GEO'
  12       include 'COMMON.NAMES'
  13       include 'COMMON.CONTROL'
  14       include 'COMMON.DISTFIT'
  15       include 'COMMON.SETUP'
  16       character*3 seq,atom,res
  17       character*80 card
  18       dimension sccor(3,20)
  19       integer rescode
  20       ibeg=1
  21       lsecondary=.false.
  22       nhfrag=0
  23       nbfrag=0
  24       do i=1,10000
  25         read (ipdbin,'(a80)',end=10) card
  26         if (card(:5).eq.'HELIX') then
  27          nhfrag=nhfrag+1
  28          lsecondary=.true.
  29          read(card(22:25),*) hfrag(1,nhfrag)
  30          read(card(34:37),*) hfrag(2,nhfrag)
  31         endif
  32         if (card(:5).eq.'SHEET') then
  33          nbfrag=nbfrag+1
  34          lsecondary=.true.
  35          read(card(24:26),*) bfrag(1,nbfrag)
  36          read(card(35:37),*) bfrag(2,nbfrag)
  37 crc----------------------------------------
  38 crc  to be corrected !!!
  39          bfrag(3,nbfrag)=bfrag(1,nbfrag)
  40          bfrag(4,nbfrag)=bfrag(2,nbfrag)
  41 crc----------------------------------------
  42         endif
  43         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
  44 C Fish out the ATOM cards.
  45         if (index(card(1:4),'ATOM').gt.0) then
  46           read (card(14:16),'(a3)') atom
  47           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  48 C Calculate the CM of the preceding residue.
  49             if (ibeg.eq.0) then
  50               if (unres_pdb) then
  51                 do j=1,3
  52                   dc(j,ires+nres)=sccor(j,iii)
  53                 enddo
  54               else
  55                 call sccenter(ires,iii,sccor)
  56               endif
  57             endif
  58 C Start new residue.
  59             read (card(24:26),*) ires
  60             read (card(18:20),'(a3)') res
  61             if (ibeg.eq.1) then
  62               ishift=ires-1
  63               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  64                 ishift=ishift-1
  65                 itype(1)=21
  66               endif
  67               ibeg=0
  68             endif
  69             ires=ires-ishift
  70             if (res.eq.'ACE') then
  71               ity=10
  72             else
  73               itype(ires)=rescode(ires,res,0)
  74             endif
  75             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  76 c            if(me.eq.king.or..not.out1file)
  77 c     &       write (iout,'(2i3,2x,a,3f8.3)')
  78 c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
  79             iii=1
  80             do j=1,3
  81               sccor(j,iii)=c(j,ires)
  82             enddo
  83           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  84      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  85             iii=iii+1
  86             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  87           endif
  88         endif
  89       enddo
  90    10 if(me.eq.king.or..not.out1file)
  91      & write (iout,'(a,i5)') ' Nres: ',ires
  92 C Calculate the CM of the last side chain.
  93       if (unres_pdb) then
  94         do j=1,3
  95           dc(j,ires+nres)=sccor(j,iii)
  96         enddo
  97       else
  98         call sccenter(ires,iii,sccor)
  99       endif
 100       nres=ires
 101       nsup=nres
 102       nstart_sup=1
 103       if (itype(nres).ne.10) then
 104         nres=nres+1
 105         itype(nres)=21
 106         if (unres_pdb) then
 107           c(1,nres)=c(1,nres-1)+3.8d0
 108           c(2,nres)=c(2,nres-1)
 109           c(3,nres)=c(3,nres-1)
 110         else
 111         do j=1,3
 112           dcj=c(j,nres-2)-c(j,nres-3)
 113           c(j,nres)=c(j,nres-1)+dcj
 114           c(j,2*nres)=c(j,nres)
 115         enddo
 116         endif
 117       endif
 118       do i=2,nres-1
 119         do j=1,3
 120           c(j,i+nres)=dc(j,i)
 121         enddo
 122       enddo
 123       do j=1,3
 124         c(j,nres+1)=c(j,1)
 125         c(j,2*nres)=c(j,nres)
 126       enddo
 127       if (itype(1).eq.21) then
 128         nsup=nsup-1
 129         nstart_sup=2
 130         if (unres_pdb) then
 131           c(1,1)=c(1,2)-3.8d0
 132           c(2,1)=c(2,2)
 133           c(3,1)=c(3,2)
 134         else
 135         do j=1,3
 136           dcj=c(j,4)-c(j,3)
 137           c(j,1)=c(j,2)-dcj
 138           c(j,nres+1)=c(j,1)
 139         enddo
 140         endif
 141       endif
 142 C Calculate internal coordinates.
 143       if(me.eq.king.or..not.out1file)then
 144        write (iout,'(a)')
 145      &   "Backbone and SC coordinates as read from the PDB"
 146        do ires=1,nres
 147         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 148      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 149      &    (c(j,nres+ires),j=1,3)
 150        enddo
 151       endif
 152       call int_from_cart(.true.,.false.)
 153       call sc_loc_geom(.false.)
 154       do i=1,nres
 155         thetaref(i)=theta(i)
 156         phiref(i)=phi(i)
 157       enddo
 158       do i=1,nres-1
 159         do j=1,3
 160           dc(j,i)=c(j,i+1)-c(j,i)
 161           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 162         enddo
 163       enddo
 164       do i=2,nres-1
 165         do j=1,3
 166           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 167           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 168         enddo
 169 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 170 c     &   vbld_inv(i+nres)
 171       enddo
 172 c      call chainbuild
 173 C Copy the coordinates to reference coordinates
 174       do i=1,2*nres
 175         do j=1,3
 176           cref(j,i)=c(j,i)
 177         enddo
 178       enddo
 179
 180
 181       do j=1,nbfrag
 182         do i=1,4
 183          bfrag(i,j)=bfrag(i,j)-ishift
 184         enddo
 185       enddo
 186
 187       do j=1,nhfrag
 188         do i=1,2
 189          hfrag(i,j)=hfrag(i,j)-ishift
 190         enddo
 191       enddo
 192
 193       return
 194       end
 195 c---------------------------------------------------------------------------
 196       subroutine int_from_cart(lside,lprn)
 197       implicit real*8 (a-h,o-z)
 198       include 'DIMENSIONS'
 199 c      include "mpif.h"
 200       include 'COMMON.LOCAL'
 201       include 'COMMON.VAR'
 202       include 'COMMON.CHAIN'
 203       include 'COMMON.INTERACT'
 204       include 'COMMON.IOUNITS'
 205       include 'COMMON.GEO'
 206       include 'COMMON.NAMES'
 207       include 'COMMON.CONTROL'
 208       include 'COMMON.SETUP'
 209       character*3 seq,atom,res
 210       character*80 card
 211       dimension sccor(3,20)
 212       integer rescode
 213       logical lside,lprn
 214       if(me.eq.king.or..not.out1file)then
 215        if (lprn) then
 216         write (iout,'(/a)')
 217      &  'Internal coordinates calculated from crystal structure.'
 218         if (lside) then
 219           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 220      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 221      & '     Beta '
 222         else
 223           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 224      & '     Gamma'
 225         endif
 226        endif
 227       endif
 228       do i=1,nres-1
 229         iti=itype(i)
 230         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 231           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 232 ctest          stop
 233         endif
 234         vbld(i+1)=dist(i,i+1)
 235         vbld_inv(i+1)=1.0d0/vbld(i+1)
 236         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 237         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 238       enddo
 239 c      if (unres_pdb) then
 240 c        if (itype(1).eq.21) then
 241 c          theta(3)=90.0d0*deg2rad
 242 c          phi(4)=180.0d0*deg2rad
 243 c          vbld(2)=3.8d0
 244 c          vbld_inv(2)=1.0d0/vbld(2)
 245 c        endif
 246 c        if (itype(nres).eq.21) then
 247 c          theta(nres)=90.0d0*deg2rad
 248 c          phi(nres)=180.0d0*deg2rad
 249 c          vbld(nres)=3.8d0
 250 c          vbld_inv(nres)=1.0d0/vbld(2)
 251 c        endif
 252 c      endif
 253       if (lside) then
 254         do i=2,nres-1
 255           do j=1,3
 256             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 257      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 258           enddo
 259           iti=itype(i)
 260           di=dist(i,nres+i)
 261 C 10/03/12 Adam: Correction for zero SC-SC bond length
 262           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
 263      &     di=dsc(itype(i))
 264           vbld(i+nres)=di
 265           if (itype(i).ne.10) then
 266             vbld_inv(i+nres)=1.0d0/di
 267           else
 268             vbld_inv(i+nres)=0.0d0
 269           endif
 270           if (iti.ne.10) then
 271             alph(i)=alpha(nres+i,i,maxres2)
 272             omeg(i)=beta(nres+i,i,maxres2,i+1)
 273           endif
 274           if(me.eq.king.or..not.out1file)then
 275            if (lprn)
 276      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 277      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 278      &     rad2deg*alph(i),rad2deg*omeg(i)
 279           endif
 280         enddo
 281       else if (lprn) then
 282         do i=2,nres
 283           iti=itype(i)
 284           if(me.eq.king.or..not.out1file)
 285      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 286      &     rad2deg*theta(i),rad2deg*phi(i)
 287         enddo
 288       endif
 289       return
 290       end
 291 c-------------------------------------------------------------------------------
 292       subroutine sc_loc_geom(lprn)
 293       implicit real*8 (a-h,o-z)
 294       include 'DIMENSIONS'
 295 c      include "mpif.h"
 296       include 'COMMON.LOCAL'
 297       include 'COMMON.VAR'
 298       include 'COMMON.CHAIN'
 299       include 'COMMON.INTERACT'
 300       include 'COMMON.IOUNITS'
 301       include 'COMMON.GEO'
 302       include 'COMMON.NAMES'
 303       include 'COMMON.CONTROL'
 304       include 'COMMON.SETUP'
 305       double precision x_prime(3),y_prime(3),z_prime(3)
 306       logical lprn
 307       do i=1,nres-1
 308         do j=1,3
 309           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 310         enddo
 311       enddo
 312       do i=2,nres-1
 313         if (itype(i).ne.10) then
 314           do j=1,3
 315             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 316           enddo
 317         else
 318           do j=1,3
 319             dc_norm(j,i+nres)=0.0d0
 320           enddo
 321         endif
 322       enddo
 323       do i=2,nres-1
 324         costtab(i+1) =dcos(theta(i+1))
 325         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 326         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 327         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 328         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 329         cosfac=dsqrt(cosfac2)
 330         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 331         sinfac=dsqrt(sinfac2)
 332         it=itype(i)
 333         if (it.ne.10) then
 334 c
 335 C  Compute the axes of tghe local cartesian coordinates system; store in
 336 c   x_prime, y_prime and z_prime
 337 c
 338         do j=1,3
 339           x_prime(j) = 0.00
 340           y_prime(j) = 0.00
 341           z_prime(j) = 0.00
 342         enddo
 343         do j = 1,3
 344           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 345           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 346         enddo
 347         call vecpr(x_prime,y_prime,z_prime)
 348 c
 349 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 350 C to local coordinate system. Store in xx, yy, zz.
 351 c
 352         xx=0.0d0
 353         yy=0.0d0
 354         zz=0.0d0
 355         do j = 1,3
 356           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 357           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 358           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 359         enddo
 360
 361         xxref(i)=xx
 362         yyref(i)=yy
 363         zzref(i)=zz
 364         else
 365         xxref(i)=0.0d0
 366         yyref(i)=0.0d0
 367         zzref(i)=0.0d0
 368         endif
 369       enddo
 370       if (lprn) then
 371         do i=2,nres
 372           iti=itype(i)
 373           if(me.eq.king.or..not.out1file)
 374      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 375      &      yyref(i),zzref(i)
 376         enddo
 377       endif
 378       return
 379       end
 380 c---------------------------------------------------------------------------
 381       subroutine sccenter(ires,nscat,sccor)
 382       implicit real*8 (a-h,o-z)
 383       include 'DIMENSIONS'
 384       include 'COMMON.CHAIN'
 385       dimension sccor(3,20)
 386       do j=1,3
 387         sccmj=0.0D0
 388         do i=1,nscat
 389           sccmj=sccmj+sccor(j,i)
 390         enddo
 391         dc(j,ires)=sccmj/nscat
 392       enddo
 393       return
 394       end
 395 c---------------------------------------------------------------------------
 396       subroutine bond_regular
 397       implicit real*8 (a-h,o-z)
 398       include 'DIMENSIONS'
 399       include 'COMMON.VAR'
 400       include 'COMMON.LOCAL'
 401       include 'COMMON.CALC'
 402       include 'COMMON.INTERACT'
 403       include 'COMMON.CHAIN'
 404       do i=1,nres-1
 405        vbld(i+1)=vbl
 406        vbld_inv(i+1)=1.0d0/vbld(i+1)
 407        vbld(i+1+nres)=dsc(itype(i+1))
 408        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 409 c       print *,vbld(i+1),vbld(i+1+nres)
 410       enddo
 411       return
 412       end
 413