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