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