source/unres/src_MD-M/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       integer rescode
  21       logical fail
  22       ibeg=1
  23       lsecondary=.false.
  24       nhfrag=0
  25       nbfrag=0
  26       do
  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') then
  46           goto 10
  47         else if (card(:3).eq.'TER') then
  48 C End current chain
  49           ires_old=ires+1
  50           itype(ires_old)=ntyp1
  51           ibeg=2
  52 c          write (iout,*) "Chain ended",ires,ishift,ires_old
  53           if (unres_pdb) then
  54             do j=1,3
  55               dc(j,ires)=sccor(j,iii)
  56             enddo
  57           else
  58             call sccenter(ires,iii,sccor)
  59           endif
  60         endif
  61 C Fish out the ATOM cards.
  62         if (index(card(1:4),'ATOM').gt.0) then
  63           read (card(14:16),'(a3)') atom
  64           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  65 C Calculate the CM of the preceding residue.
  66             if (ibeg.eq.0) then
  67               if (unres_pdb) then
  68                 do j=1,3
  69                   dc(j,ires+nres)=sccor(j,iii)
  70                 enddo
  71               else
  72                 call sccenter(ires,iii,sccor)
  73               endif
  74             endif
  75 C Start new residue.
  76 c            write (iout,'(a80)') card
  77             read (card(24:26),*) ires
  78             read (card(18:20),'(a3)') res
  79             if (ibeg.eq.1) then
  80               ishift=ires-1
  81               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  82                 ishift=ishift-1
  83                 itype(1)=ntyp1
  84               endif
  85 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
  86               ibeg=0
  87             else if (ibeg.eq.2) then
  88 c Start a new chain
  89               ishift=-ires_old+ires-1
  90 c              write (iout,*) "New chain started",ires,ishift
  91               ibeg=0
  92             endif
  93             ires=ires-ishift
  94 c            write (2,*) "ires",ires," ishift",ishift
  95             if (res.eq.'ACE') then
  96               itype(ires)=10
  97             else
  98               itype(ires)=rescode(ires,res,0)
  99             endif
 100             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 101             if(me.eq.king.or..not.out1file)
 102      &       write (iout,'(2i3,2x,a,3f8.3)')
 103      &       ires,itype(ires),res,(c(j,ires),j=1,3)
 104             iii=1
 105             do j=1,3
 106               sccor(j,iii)=c(j,ires)
 107             enddo
 108           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
 109      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
 110             iii=iii+1
 111             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 112           endif
 113         endif
 114       enddo
 115    10 if(me.eq.king.or..not.out1file)
 116      & write (iout,'(a,i5)') ' Nres: ',ires
 117 C Calculate dummy residue coordinates inside the "chain" of a multichain
 118 C system
 119       nres=ires
 120       do i=2,nres-1
 121 c        write (iout,*) i,itype(i)
 122         if (itype(i).eq.ntyp1) then
 123 c          write (iout,*) "dummy",i,itype(i)
 124           do j=1,3
 125             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
 126 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
 127             dc(j,i)=c(j,i)
 128           enddo
 129         endif
 130       enddo
 131 C Calculate the CM of the last side chain.
 132       if (unres_pdb) then
 133         do j=1,3
 134           dc(j,ires)=sccor(j,iii)
 135         enddo
 136       else
 137         call sccenter(ires,iii,sccor)
 138       endif
 139       nsup=nres
 140       nstart_sup=1
 141       if (itype(nres).ne.10) then
 142         nres=nres+1
 143         itype(nres)=ntyp1
 144         if (unres_pdb) then
 145 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 146           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 147           if (fail) then
 148             e2(1)=0.0d0
 149             e2(2)=1.0d0
 150             e2(3)=0.0d0
 151           endif
 152           do j=1,3
 153             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
 154           enddo
 155         else
 156         do j=1,3
 157           dcj=c(j,nres-2)-c(j,nres-3)
 158           c(j,nres)=c(j,nres-1)+dcj
 159           c(j,2*nres)=c(j,nres)
 160         enddo
 161         endif
 162       endif
 163       do i=2,nres-1
 164         do j=1,3
 165           c(j,i+nres)=dc(j,i)
 166         enddo
 167       enddo
 168       do j=1,3
 169         c(j,nres+1)=c(j,1)
 170         c(j,2*nres)=c(j,nres)
 171       enddo
 172       if (itype(1).eq.ntyp1) then
 173         nsup=nsup-1
 174         nstart_sup=2
 175         if (unres_pdb) then
 176 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 177           call refsys(2,3,4,e1,e2,e3,fail)
 178           if (fail) then
 179             e2(1)=0.0d0
 180             e2(2)=1.0d0
 181             e2(3)=0.0d0
 182           endif
 183           do j=1,3
 184             c(j,1)=c(j,2)-3.8d0*e2(j)
 185           enddo
 186         else
 187         do j=1,3
 188           dcj=c(j,4)-c(j,3)
 189           c(j,1)=c(j,2)-dcj
 190           c(j,nres+1)=c(j,1)
 191         enddo
 192         endif
 193       endif
 194 C Calculate internal coordinates.
 195       if(me.eq.king.or..not.out1file)then
 196        do ires=1,nres
 197         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 198      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 199      &    (c(j,nres+ires),j=1,3)
 200        enddo
 201       endif
 202       call int_from_cart(.true.,.false.)
 203       call sc_loc_geom(.true.)
 204       do i=1,nres
 205         thetaref(i)=theta(i)
 206         phiref(i)=phi(i)
 207       enddo
 208       do i=1,nres-1
 209         do j=1,3
 210           dc(j,i)=c(j,i+1)-c(j,i)
 211           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 212         enddo
 213       enddo
 214       do i=2,nres-1
 215         do j=1,3
 216           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 217           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 218         enddo
 219 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 220 c     &   vbld_inv(i+nres)
 221       enddo
 222 c      call chainbuild
 223 C Copy the coordinates to reference coordinates
 224 C Splits to single chain if occurs
 225       kkk=1
 226       lll=0
 227       cou=1
 228       do i=1,nres
 229       lll=lll+1
 230 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 231       if (i.gt.1) then
 232       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
 233       chain_length=lll-1
 234       kkk=kkk+1
 235 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 236       lll=1
 237       endif
 238       endif
 239         do j=1,3
 240           cref(j,i,cou)=c(j,i)
 241           cref(j,i+nres,cou)=c(j,i+nres)
 242           if (i.le.nres) then
 243           chain_rep(j,lll,kkk)=c(j,i)
 244           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
 245           endif
 246          enddo
 247       enddo
 248       write (iout,*) chain_length
 249       if (chain_length.eq.0) chain_length=nres
 250       do j=1,3
 251       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 252       chain_rep(j,chain_length+nres,symetr)
 253      &=chain_rep(j,chain_length+nres,1)
 254       enddo
 255 c diagnostic
 256 c       write (iout,*) "spraw lancuchy",chain_length,symetr
 257 c       do i=1,4
 258 c         do kkk=1,chain_length
 259 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
 260 c         enddo
 261 c        enddo
 262 c enddiagnostic
 263 C makes copy of chains
 264         write (iout,*) "symetr", symetr
 265
 266       if (symetr.gt.1) then
 267        call permut(symetr)
 268        nperm=1
 269        do i=1,symetr
 270        nperm=nperm*i
 271        enddo
 272        do i=1,nperm
 273        write(iout,*) (tabperm(i,kkk),kkk=1,4)
 274        enddo
 275        do i=1,nperm
 276         cou=0
 277         do kkk=1,symetr
 278          icha=tabperm(i,kkk)
 279 c         write (iout,*) i,icha
 280          do lll=1,chain_length
 281           cou=cou+1
 282            if (cou.le.nres) then
 283            do j=1,3
 284             kupa=mod(lll,chain_length)
 285             iprzes=(kkk-1)*chain_length+lll
 286             if (kupa.eq.0) kupa=chain_length
 287 c            write (iout,*) "kupa", kupa
 288             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 289             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 290           enddo
 291           endif
 292          enddo
 293         enddo
 294        enddo
 295        endif
 296 C-koniec robienia kopii
 297 c diag
 298       do kkk=1,nperm
 299       write (iout,*) "nowa struktura", nperm
 300       do i=1,nres
 301         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 302      &cref(2,i,kkk),
 303      &cref(3,i,kkk),cref(1,nres+i,kkk),
 304      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 305       enddo
 306   100 format (//'              alpha-carbon coordinates       ',
 307      &          '     centroid coordinates'/
 308      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 309      &                          10X,'X',11X,'Y',11X,'Z')
 310   110 format (a,'(',i3,')',6f12.5)
 311
 312       enddo
 313 cc enddiag
 314       do j=1,nbfrag
 315         do i=1,4
 316          bfrag(i,j)=bfrag(i,j)-ishift
 317         enddo
 318       enddo
 319
 320       do j=1,nhfrag
 321         do i=1,2
 322          hfrag(i,j)=hfrag(i,j)-ishift
 323         enddo
 324       enddo
 325
 326       return
 327       end
 328 c---------------------------------------------------------------------------
 329       subroutine int_from_cart(lside,lprn)
 330       implicit real*8 (a-h,o-z)
 331       include 'DIMENSIONS'
 332 #ifdef MPI
 333       include "mpif.h"
 334 #endif
 335       include 'COMMON.LOCAL'
 336       include 'COMMON.VAR'
 337       include 'COMMON.CHAIN'
 338       include 'COMMON.INTERACT'
 339       include 'COMMON.IOUNITS'
 340       include 'COMMON.GEO'
 341       include 'COMMON.NAMES'
 342       include 'COMMON.CONTROL'
 343       include 'COMMON.SETUP'
 344       character*3 seq,atom,res
 345       character*80 card
 346       dimension sccor(3,20)
 347       integer rescode
 348       logical lside,lprn
 349 #ifdef MPI
 350       if(me.eq.king.or..not.out1file)then
 351 #endif
 352        if (lprn) then
 353         write (iout,'(/a)')
 354      &  'Internal coordinates calculated from crystal structure.'
 355         if (lside) then
 356           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 357      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 358      & '     Omega'
 359         else
 360           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 361      & '       Phi'
 362         endif
 363        endif
 364 #ifdef MPI
 365       endif
 366 #endif
 367       do i=1,nres-1
 368         iti=itype(i)
 369         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 370      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 371           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 372 ctest          stop
 373         endif
 374         vbld(i+1)=dist(i,i+1)
 375         vbld_inv(i+1)=1.0d0/vbld(i+1)
 376         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 377         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 378       enddo
 379 c      if (unres_pdb) then
 380 c        if (itype(1).eq.21) then
 381 c          theta(3)=90.0d0*deg2rad
 382 c          phi(4)=180.0d0*deg2rad
 383 c          vbld(2)=3.8d0
 384 c          vbld_inv(2)=1.0d0/vbld(2)
 385 c        endif
 386 c        if (itype(nres).eq.21) then
 387 c          theta(nres)=90.0d0*deg2rad
 388 c          phi(nres)=180.0d0*deg2rad
 389 c          vbld(nres)=3.8d0
 390 c          vbld_inv(nres)=1.0d0/vbld(2)
 391 c        endif
 392 c      endif
 393       if (lside) then
 394         do i=2,nres-1
 395           do j=1,3
 396             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 397      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 398           enddo
 399           iti=itype(i)
 400           di=dist(i,nres+i)
 401           vbld(i+nres)=di
 402           if (itype(i).ne.10) then
 403             vbld_inv(i+nres)=1.0d0/di
 404           else
 405             vbld_inv(i+nres)=0.0d0
 406           endif
 407           if (iti.ne.10) then
 408             alph(i)=alpha(nres+i,i,maxres2)
 409             omeg(i)=beta(nres+i,i,maxres2,i+1)
 410           endif
 411           if(me.eq.king.or..not.out1file)then
 412            if (lprn)
 413      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 414      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 415      &     rad2deg*alph(i),rad2deg*omeg(i)
 416           endif
 417         enddo
 418       else if (lprn) then
 419         do i=2,nres
 420           iti=itype(i)
 421           if(me.eq.king.or..not.out1file)
 422      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 423      &     rad2deg*theta(i),rad2deg*phi(i)
 424         enddo
 425       endif
 426       return
 427       end
 428 c-------------------------------------------------------------------------------
 429       subroutine sc_loc_geom(lprn)
 430       implicit real*8 (a-h,o-z)
 431       include 'DIMENSIONS'
 432 #ifdef MPI
 433       include "mpif.h"
 434 #endif
 435       include 'COMMON.LOCAL'
 436       include 'COMMON.VAR'
 437       include 'COMMON.CHAIN'
 438       include 'COMMON.INTERACT'
 439       include 'COMMON.IOUNITS'
 440       include 'COMMON.GEO'
 441       include 'COMMON.NAMES'
 442       include 'COMMON.CONTROL'
 443       include 'COMMON.SETUP'
 444       double precision x_prime(3),y_prime(3),z_prime(3)
 445       logical lprn
 446       do i=1,nres-1
 447         do j=1,3
 448           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 449         enddo
 450       enddo
 451       do i=2,nres-1
 452         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 453           do j=1,3
 454             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 455           enddo
 456         else
 457           do j=1,3
 458             dc_norm(j,i+nres)=0.0d0
 459           enddo
 460         endif
 461       enddo
 462       do i=2,nres-1
 463         costtab(i+1) =dcos(theta(i+1))
 464         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 465         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 466         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 467         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 468         cosfac=dsqrt(cosfac2)
 469         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 470         sinfac=dsqrt(sinfac2)
 471         it=itype(i)
 472         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 473 c
 474 C  Compute the axes of tghe local cartesian coordinates system; store in
 475 c   x_prime, y_prime and z_prime
 476 c
 477         do j=1,3
 478           x_prime(j) = 0.00
 479           y_prime(j) = 0.00
 480           z_prime(j) = 0.00
 481         enddo
 482         do j = 1,3
 483           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 484           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 485         enddo
 486         call vecpr(x_prime,y_prime,z_prime)
 487 c
 488 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 489 C to local coordinate system. Store in xx, yy, zz.
 490 c
 491         xx=0.0d0
 492         yy=0.0d0
 493         zz=0.0d0
 494         do j = 1,3
 495           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 496           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 497           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 498         enddo
 499
 500         xxref(i)=xx
 501         yyref(i)=yy
 502         zzref(i)=zz
 503         else
 504         xxref(i)=0.0d0
 505         yyref(i)=0.0d0
 506         zzref(i)=0.0d0
 507         endif
 508       enddo
 509       if (lprn) then
 510         do i=2,nres
 511           iti=itype(i)
 512 #ifdef MPI
 513           if(me.eq.king.or..not.out1file)
 514      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 515      &      yyref(i),zzref(i)
 516 #else
 517           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
 518      &     zzref(i)
 519 #endif
 520         enddo
 521       endif
 522       return
 523       end
 524 c---------------------------------------------------------------------------
 525       subroutine sccenter(ires,nscat,sccor)
 526       implicit real*8 (a-h,o-z)
 527       include 'DIMENSIONS'
 528       include 'COMMON.CHAIN'
 529       dimension sccor(3,20)
 530       do j=1,3
 531         sccmj=0.0D0
 532         do i=1,nscat
 533           sccmj=sccmj+sccor(j,i)
 534         enddo
 535         dc(j,ires)=sccmj/nscat
 536       enddo
 537       return
 538       end
 539 c---------------------------------------------------------------------------
 540       subroutine bond_regular
 541       implicit real*8 (a-h,o-z)
 542       include 'DIMENSIONS'
 543       include 'COMMON.VAR'
 544       include 'COMMON.LOCAL'
 545       include 'COMMON.CALC'
 546       include 'COMMON.INTERACT'
 547       include 'COMMON.CHAIN'
 548       do i=1,nres-1
 549        vbld(i+1)=vbl
 550        vbld_inv(i+1)=1.0d0/vbld(i+1)
 551        vbld(i+1+nres)=dsc(iabs(itype(i+1)))
 552        vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
 553 c       print *,vbld(i+1),vbld(i+1+nres)
 554       enddo
 555       return
 556       end
 557