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