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       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) phi(i+1)=beta(i-2,i-1,i,i+1)
 254       enddo
 255 c      if (unres_pdb) then
 256 c        if (itype(1).eq.21) then
 257 c          theta(3)=90.0d0*deg2rad
 258 c          phi(4)=180.0d0*deg2rad
 259 c          vbld(2)=3.8d0
 260 c          vbld_inv(2)=1.0d0/vbld(2)
 261 c        endif
 262 c        if (itype(nres).eq.21) then
 263 c          theta(nres)=90.0d0*deg2rad
 264 c          phi(nres)=180.0d0*deg2rad
 265 c          vbld(nres)=3.8d0
 266 c          vbld_inv(nres)=1.0d0/vbld(2)
 267 c        endif
 268 c      endif
 269       if (lside) then
 270         do i=2,nres-1
 271           do j=1,3
 272             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 273      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 274           enddo
 275           iti=itype(i)
 276           di=dist(i,nres+i)
 277 C 9/29/12 Adam: Correction for zero SC-SC bond length
 278           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
 279      &     di=dsc(itype(i))
 280           vbld(i+nres)=di
 281           if (itype(i).ne.10) then
 282             vbld_inv(i+nres)=1.0d0/di
 283           else
 284             vbld_inv(i+nres)=0.0d0
 285           endif
 286           if (iti.ne.10) then
 287             alph(i)=alpha(nres+i,i,maxres2)
 288             omeg(i)=beta(nres+i,i,maxres2,i+1)
 289           endif
 290           if(me.eq.king.or..not.out1file)then
 291            if (lprn)
 292      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 293      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 294      &     rad2deg*alph(i),rad2deg*omeg(i)
 295           endif
 296         enddo
 297       else if (lprn) then
 298         do i=2,nres
 299           iti=itype(i)
 300           if(me.eq.king.or..not.out1file)
 301      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 302      &     rad2deg*theta(i),rad2deg*phi(i)
 303         enddo
 304       endif
 305       return
 306       end
 307 c-------------------------------------------------------------------------------
 308       subroutine sc_loc_geom(lprn)
 309       implicit real*8 (a-h,o-z)
 310       include 'DIMENSIONS'
 311 c      include "mpif.h"
 312       include 'COMMON.LOCAL'
 313       include 'COMMON.VAR'
 314       include 'COMMON.CHAIN'
 315       include 'COMMON.INTERACT'
 316       include 'COMMON.IOUNITS'
 317       include 'COMMON.GEO'
 318       include 'COMMON.NAMES'
 319       include 'COMMON.CONTROL'
 320       include 'COMMON.SETUP'
 321       double precision x_prime(3),y_prime(3),z_prime(3)
 322       logical lprn
 323       do i=1,nres-1
 324         do j=1,3
 325           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 326         enddo
 327       enddo
 328       do i=2,nres-1
 329         if (itype(i).ne.10) then
 330           do j=1,3
 331             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 332           enddo
 333         else
 334           do j=1,3
 335             dc_norm(j,i+nres)=0.0d0
 336           enddo
 337         endif
 338       enddo
 339       do i=2,nres-1
 340         costtab(i+1) =dcos(theta(i+1))
 341         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 342         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 343         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 344         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 345         cosfac=dsqrt(cosfac2)
 346         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 347         sinfac=dsqrt(sinfac2)
 348         it=itype(i)
 349         if (it.ne.10) then
 350 c
 351 C  Compute the axes of tghe local cartesian coordinates system; store in
 352 c   x_prime, y_prime and z_prime
 353 c
 354         do j=1,3
 355           x_prime(j) = 0.00
 356           y_prime(j) = 0.00
 357           z_prime(j) = 0.00
 358         enddo
 359         do j = 1,3
 360           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 361           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 362         enddo
 363         call vecpr(x_prime,y_prime,z_prime)
 364 c
 365 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 366 C to local coordinate system. Store in xx, yy, zz.
 367 c
 368         xx=0.0d0
 369         yy=0.0d0
 370         zz=0.0d0
 371         do j = 1,3
 372           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 373           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 374           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 375         enddo
 376
 377         xxref(i)=xx
 378         yyref(i)=yy
 379         zzref(i)=zz
 380         else
 381         xxref(i)=0.0d0
 382         yyref(i)=0.0d0
 383         zzref(i)=0.0d0
 384         endif
 385       enddo
 386       if (lprn) then
 387         do i=2,nres
 388           iti=itype(i)
 389           if(me.eq.king.or..not.out1file)
 390      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 391      &      yyref(i),zzref(i)
 392         enddo
 393       endif
 394       return
 395       end
 396 c---------------------------------------------------------------------------
 397       subroutine sccenter(ires,nscat,sccor)
 398       implicit real*8 (a-h,o-z)
 399       include 'DIMENSIONS'
 400       include 'COMMON.CHAIN'
 401       dimension sccor(3,20)
 402       do j=1,3
 403         sccmj=0.0D0
 404         do i=1,nscat
 405           sccmj=sccmj+sccor(j,i)
 406         enddo
 407         dc(j,ires)=sccmj/nscat
 408       enddo
 409       return
 410       end
 411 c---------------------------------------------------------------------------
 412       subroutine bond_regular
 413       implicit real*8 (a-h,o-z)
 414       include 'DIMENSIONS'
 415       include 'COMMON.VAR'
 416       include 'COMMON.LOCAL'
 417       include 'COMMON.CALC'
 418       include 'COMMON.INTERACT'
 419       include 'COMMON.CHAIN'
 420       do i=1,nres-1
 421        vbld(i+1)=vbl
 422        vbld_inv(i+1)=1.0d0/vbld(i+1)
 423        vbld(i+1+nres)=dsc(itype(i+1))
 424        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 425 c       print *,vbld(i+1),vbld(i+1+nres)
 426       enddo
 427       return
 428       end