source/unres/src_Eshel/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,maxres
  27         itype(i)=21
  28         do j=1,3
  29           c(j,i)=0.0d0
  30           c(j,i+nres)=0.0d0
  31         enddo
  32       enddo
  33       do i=1,10000
  34         read (ipdbin,'(a80)',end=10) card
  35         if (card(:5).eq.'HELIX') then
  36          nhfrag=nhfrag+1
  37          lsecondary=.true.
  38          read(card(22:25),*) hfrag(1,nhfrag)
  39          read(card(34:37),*) hfrag(2,nhfrag)
  40         endif
  41         if (card(:5).eq.'SHEET') then
  42          nbfrag=nbfrag+1
  43          lsecondary=.true.
  44          read(card(24:26),*) bfrag(1,nbfrag)
  45          read(card(35:37),*) bfrag(2,nbfrag)
  46 crc----------------------------------------
  47 crc  to be corrected !!!
  48          bfrag(3,nbfrag)=bfrag(1,nbfrag)
  49          bfrag(4,nbfrag)=bfrag(2,nbfrag)
  50 crc----------------------------------------
  51         endif
  52         if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
  53 C Fish out the ATOM cards.
  54         if (index(card(1:4),'ATOM').gt.0) then
  55           read (card(14:16),'(a3)') atom
  56           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  57 C Calculate the CM of the preceding residue.
  58             if (ibeg.eq.0) then
  59               if (unres_pdb) then
  60                 do j=1,3
  61                   dc(j,ires+nres)=sccor(j,iii)
  62                 enddo
  63               else
  64                 call sccenter(ires,iii,sccor)
  65               endif
  66             endif
  67 C Start new residue.
  68             read (card(24:26),*) ires
  69             read (card(18:20),'(a3)') res
  70             if (ibeg.eq.1) then
  71               ishift=ires-1
  72               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  73                 ishift=ishift-1
  74                 itype(1)=21
  75               endif
  76               ibeg=0
  77             endif
  78             ires=ires-ishift
  79             if (res.eq.'ACE') then
  80               ity=10
  81             else
  82               itype(ires)=rescode(ires,res,0)
  83             endif
  84             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  85 c            if(me.eq.king.or..not.out1file)
  86 c     &       write (iout,'(2i3,2x,a,3f8.3)')
  87 c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
  88             iii=1
  89             do j=1,3
  90               sccor(j,iii)=c(j,ires)
  91             enddo
  92           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  93      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  94             iii=iii+1
  95             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  96           endif
  97         endif
  98       enddo
  99    10 if(me.eq.king.or..not.out1file)
 100      & write (iout,'(a,i5)') ' Nres: ',ires
 101 C Calculate the CM of the last side chain.
 102       if (unres_pdb) then
 103         do j=1,3
 104           dc(j,ires+nres)=sccor(j,iii)
 105         enddo
 106       else if (.not.catrace) then
 107         call sccenter(ires,iii,sccor)
 108       endif
 109       nres=ires
 110       nsup=nres
 111       nstart_sup=1
 112       if (itype(nres).ne.10) then
 113         nres=nres+1
 114         itype(nres)=21
 115         if (unres_pdb) then
 116 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 117           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 118           if (fail) then
 119             e2(1)=0.0d0
 120             e2(2)=1.0d0
 121             e2(3)=0.0d0
 122           endif
 123           do j=1,3
 124             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
 125           enddo
 126         else if (.not.catrace) then
 127         do j=1,3
 128           dcj=c(j,nres-2)-c(j,nres-3)
 129           c(j,nres)=c(j,nres-1)+dcj
 130           c(j,2*nres)=c(j,nres)
 131         enddo
 132         endif
 133       endif
 134       do i=2,nres-1
 135         do j=1,3
 136           c(j,i+nres)=dc(j,i)
 137         enddo
 138       enddo
 139       do j=1,3
 140         c(j,nres+1)=c(j,1)
 141         c(j,2*nres)=c(j,nres)
 142       enddo
 143       if (itype(1).eq.21) then
 144         nsup=nsup-1
 145         nstart_sup=2
 146         if (unres_pdb) then
 147 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 148           call refsys(2,3,4,e1,e2,e3,fail)
 149           if (fail) then
 150             e2(1)=0.0d0
 151             e2(2)=1.0d0
 152             e2(3)=0.0d0
 153           endif
 154           do j=1,3
 155             c(j,1)=c(j,2)-3.8d0*e2(j)
 156           enddo
 157         else if (.not.catrace) then
 158         do j=1,3
 159           dcj=c(j,4)-c(j,3)
 160           c(j,1)=c(j,2)-dcj
 161           c(j,nres+1)=c(j,1)
 162         enddo
 163         endif
 164       endif
 165 C Calculate internal coordinates.
 166       if(me.eq.king.or..not.out1file)then
 167        write (iout,'(a)')
 168      &   "Backbone and SC coordinates as read from the PDB"
 169        do ires=1,nres
 170         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 171      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 172      &    (c(j,nres+ires),j=1,3)
 173        enddo
 174       endif
 175       call int_from_cart(.not.catrace,.false.)
 176       if (.not.catrace) call sc_loc_geom(.false.)
 177       do i=1,nres
 178         thetaref(i)=theta(i)
 179         phiref(i)=phi(i)
 180       enddo
 181       do i=1,nres-1
 182         do j=1,3
 183           dc(j,i)=c(j,i+1)-c(j,i)
 184           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 185         enddo
 186       enddo
 187       do i=2,nres-1
 188         do j=1,3
 189           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 190           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 191         enddo
 192 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 193 c     &   vbld_inv(i+nres)
 194       enddo
 195 c      call chainbuild
 196 C Copy the coordinates to reference coordinates
 197       do i=1,2*nres
 198         do j=1,3
 199           cref(j,i)=c(j,i)
 200         enddo
 201       enddo
 202
 203
 204       do j=1,nbfrag
 205         do i=1,4
 206          bfrag(i,j)=bfrag(i,j)-ishift
 207         enddo
 208       enddo
 209
 210       do j=1,nhfrag
 211         do i=1,2
 212          hfrag(i,j)=hfrag(i,j)-ishift
 213         enddo
 214       enddo
 215
 216       return
 217       end
 218 c---------------------------------------------------------------------------
 219       subroutine int_from_cart(lside,lprn)
 220       implicit real*8 (a-h,o-z)
 221       include 'DIMENSIONS'
 222 #ifdef MPI
 223       include "mpif.h"
 224 #endif
 225       include 'COMMON.LOCAL'
 226       include 'COMMON.VAR'
 227       include 'COMMON.CHAIN'
 228       include 'COMMON.INTERACT'
 229       include 'COMMON.IOUNITS'
 230       include 'COMMON.GEO'
 231       include 'COMMON.NAMES'
 232       include 'COMMON.CONTROL'
 233       include 'COMMON.SETUP'
 234       character*3 seq,atom,res
 235       character*80 card
 236       dimension sccor(3,20)
 237       integer rescode
 238       logical lside,lprn
 239       if(me.eq.king.or..not.out1file)then
 240        if (lprn) then
 241         write (iout,'(/a)')
 242      &  'Internal coordinates calculated from crystal structure.'
 243         if (lside) then
 244           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 245      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 246      & '     Beta '
 247         else
 248           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 249      & '     Gamma'
 250         endif
 251        endif
 252       endif
 253       do i=1,nres-1
 254         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
 255         iti=itype(i)
 256         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 257           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 258 ctest          stop
 259         endif
 260         vbld(i+1)=dist(i,i+1)
 261         vbld_inv(i+1)=1.0d0/vbld(i+1)
 262         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 263         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 264       enddo
 265 c      if (unres_pdb) then
 266 c        if (itype(1).eq.21) then
 267 c          theta(3)=90.0d0*deg2rad
 268 c          phi(4)=180.0d0*deg2rad
 269 c          vbld(2)=3.8d0
 270 c          vbld_inv(2)=1.0d0/vbld(2)
 271 c        endif
 272 c        if (itype(nres).eq.21) then
 273 c          theta(nres)=90.0d0*deg2rad
 274 c          phi(nres)=180.0d0*deg2rad
 275 c          vbld(nres)=3.8d0
 276 c          vbld_inv(nres)=1.0d0/vbld(2)
 277 c        endif
 278 c      endif
 279       if (lside) then
 280         do i=2,nres-1
 281           if (itype(i).eq.ntyp1) cycle
 282           do j=1,3
 283             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 284      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 285           enddo
 286           iti=itype(i)
 287           di=dist(i,nres+i)
 288 C 10/03/12 Adam: Correction for zero SC-SC bond length
 289           if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
 290      &     di=dsc(itype(i))
 291           vbld(i+nres)=di
 292           if (itype(i).ne.10) then
 293             vbld_inv(i+nres)=1.0d0/di
 294           else
 295             vbld_inv(i+nres)=0.0d0
 296           endif
 297           if (iti.ne.10) then
 298             alph(i)=alpha(nres+i,i,maxres2)
 299             omeg(i)=beta(nres+i,i,maxres2,i+1)
 300           endif
 301           if(me.eq.king.or..not.out1file)then
 302            if (lprn)
 303      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 304      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 305      &     rad2deg*alph(i),rad2deg*omeg(i)
 306           endif
 307         enddo
 308       else if (lprn) then
 309         do i=2,nres
 310           iti=itype(i)
 311           if(me.eq.king.or..not.out1file)
 312      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 313      &     rad2deg*theta(i),rad2deg*phi(i)
 314         enddo
 315       endif
 316       return
 317       end
 318 c-------------------------------------------------------------------------------
 319       subroutine sc_loc_geom(lprn)
 320       implicit real*8 (a-h,o-z)
 321       include 'DIMENSIONS'
 322 #ifdef MPI
 323       include "mpif.h"
 324 #endif
 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