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