source/unres/src_MD/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 #ifdef MPI
 216       include "mpif.h"
 217 #endif
 218       include 'COMMON.LOCAL'
 219       include 'COMMON.VAR'
 220       include 'COMMON.CHAIN'
 221       include 'COMMON.INTERACT'
 222       include 'COMMON.IOUNITS'
 223       include 'COMMON.GEO'
 224       include 'COMMON.NAMES'
 225       include 'COMMON.CONTROL'
 226       include 'COMMON.SETUP'
 227       character*3 seq,atom,res
 228       character*80 card
 229       dimension sccor(3,20)
 230       integer rescode
 231       logical lside,lprn
 232       if(me.eq.king.or..not.out1file)then
 233        if (lprn) then
 234         write (iout,'(/a)')
 235      &  'Internal coordinates calculated from crystal structure.'
 236         if (lside) then
 237           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 238      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 239      & '     Beta '
 240         else
 241           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 242      & '     Gamma'
 243         endif
 244        endif
 245       endif
 246       do i=1,nres-1
 247         iti=itype(i)
 248         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 249           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 250 ctest          stop
 251         endif
 252         vbld(i+1)=dist(i,i+1)
 253         vbld_inv(i+1)=1.0d0/vbld(i+1)
 254         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 255         if (i.gt.2) then
 256         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
 257         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
 258          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
 259         endif
 260         if (itype(i-1).ne.10) then
 261          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
 262          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
 263          omicron(2,i)=alpha(i-1+nres,i-1,i)
 264         endif
 265         if (itype(i).ne.10) then
 266          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
 267         endif
 268         endif
 269       enddo
 270 c      if (unres_pdb) then
 271 c        if (itype(1).eq.21) then
 272 c          theta(3)=90.0d0*deg2rad
 273 c          phi(4)=180.0d0*deg2rad
 274 c          vbld(2)=3.8d0
 275 c          vbld_inv(2)=1.0d0/vbld(2)
 276 c        endif
 277 c        if (itype(nres).eq.21) then
 278 c          theta(nres)=90.0d0*deg2rad
 279 c          phi(nres)=180.0d0*deg2rad
 280 c          vbld(nres)=3.8d0
 281 c          vbld_inv(nres)=1.0d0/vbld(2)
 282 c        endif
 283 c      endif
 284       if (lside) then
 285         do i=2,nres-1
 286           do j=1,3
 287             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 288      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 289           enddo
 290           iti=itype(i)
 291           di=dist(i,nres+i)
 292 C 10/03/12 Adam: Correction for zero SC-SC bond length
 293           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
 294      &     di=dsc(itype(i))
 295           vbld(i+nres)=di
 296           if (itype(i).ne.10) then
 297             vbld_inv(i+nres)=1.0d0/di
 298           else
 299             vbld_inv(i+nres)=0.0d0
 300           endif
 301           if (iti.ne.10) then
 302             alph(i)=alpha(nres+i,i,maxres2)
 303             omeg(i)=beta(nres+i,i,maxres2,i+1)
 304           endif
 305           if(me.eq.king.or..not.out1file)then
 306            if (lprn)
 307      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 308      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 309      &     rad2deg*alph(i),rad2deg*omeg(i)
 310           endif
 311         enddo
 312       else if (lprn) then
 313         do i=2,nres
 314           iti=itype(i)
 315           if(me.eq.king.or..not.out1file)
 316      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 317      &     rad2deg*theta(i),rad2deg*phi(i)
 318         enddo
 319       endif
 320       return
 321       end
 322 c-------------------------------------------------------------------------------
 323       subroutine sc_loc_geom(lprn)
 324       implicit real*8 (a-h,o-z)
 325       include 'DIMENSIONS'
 326 #ifdef MPI
 327       include "mpif.h"
 328 #endif
 329       include 'COMMON.LOCAL'
 330       include 'COMMON.VAR'
 331       include 'COMMON.CHAIN'
 332       include 'COMMON.INTERACT'
 333       include 'COMMON.IOUNITS'
 334       include 'COMMON.GEO'
 335       include 'COMMON.NAMES'
 336       include 'COMMON.CONTROL'
 337       include 'COMMON.SETUP'
 338       double precision x_prime(3),y_prime(3),z_prime(3)
 339       logical lprn
 340       do i=1,nres-1
 341         do j=1,3
 342           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 343         enddo
 344       enddo
 345       do i=2,nres-1
 346         if (itype(i).ne.10) then
 347           do j=1,3
 348             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 349           enddo
 350         else
 351           do j=1,3
 352             dc_norm(j,i+nres)=0.0d0
 353           enddo
 354         endif
 355       enddo
 356       do i=2,nres-1
 357         costtab(i+1) =dcos(theta(i+1))
 358         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 359         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 360         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 361         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 362         cosfac=dsqrt(cosfac2)
 363         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 364         sinfac=dsqrt(sinfac2)
 365         it=itype(i)
 366         if (it.ne.10) then
 367 c
 368 C  Compute the axes of tghe local cartesian coordinates system; store in
 369 c   x_prime, y_prime and z_prime
 370 c
 371         do j=1,3
 372           x_prime(j) = 0.00
 373           y_prime(j) = 0.00
 374           z_prime(j) = 0.00
 375         enddo
 376         do j = 1,3
 377           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 378           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 379         enddo
 380         call vecpr(x_prime,y_prime,z_prime)
 381 c
 382 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 383 C to local coordinate system. Store in xx, yy, zz.
 384 c
 385         xx=0.0d0
 386         yy=0.0d0
 387         zz=0.0d0
 388         do j = 1,3
 389           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 390           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 391           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 392         enddo
 393
 394         xxref(i)=xx
 395         yyref(i)=yy
 396         zzref(i)=zz
 397         else
 398         xxref(i)=0.0d0
 399         yyref(i)=0.0d0
 400         zzref(i)=0.0d0
 401         endif
 402       enddo
 403       if (lprn) then
 404         do i=2,nres
 405           iti=itype(i)
 406           if(me.eq.king.or..not.out1file)
 407      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 408      &      yyref(i),zzref(i)
 409         enddo
 410       endif
 411       return
 412       end
 413 c---------------------------------------------------------------------------
 414       subroutine sccenter(ires,nscat,sccor)
 415       implicit real*8 (a-h,o-z)
 416       include 'DIMENSIONS'
 417       include 'COMMON.CHAIN'
 418       dimension sccor(3,20)
 419       do j=1,3
 420         sccmj=0.0D0
 421         do i=1,nscat
 422           sccmj=sccmj+sccor(j,i)
 423         enddo
 424         dc(j,ires)=sccmj/nscat
 425       enddo
 426       return
 427       end
 428 c---------------------------------------------------------------------------
 429       subroutine bond_regular
 430       implicit real*8 (a-h,o-z)
 431       include 'DIMENSIONS'
 432       include 'COMMON.VAR'
 433       include 'COMMON.LOCAL'
 434       include 'COMMON.CALC'
 435       include 'COMMON.INTERACT'
 436       include 'COMMON.CHAIN'
 437       do i=1,nres-1
 438        vbld(i+1)=vbl
 439        vbld_inv(i+1)=1.0d0/vbld(i+1)
 440        vbld(i+1+nres)=dsc(itype(i+1))
 441        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 442 c       print *,vbld(i+1),vbld(i+1+nres)
 443       enddo
 444       return
 445       end