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)=21
  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)=21
  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               ity=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.21) 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)=21
 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.21) 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.21)) 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       do j=1,3
 249       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 250       chain_rep(j,chain_length+nres,symetr)
 251      &=chain_rep(j,chain_length+nres,1)
 252       enddo
 253 c diagnostic
 254 c       write (iout,*) "spraw lancuchy",chain_length,symetr
 255 c       do i=1,4
 256 c         do kkk=1,chain_length
 257 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
 258 c         enddo
 259 c        enddo
 260 c enddiagnostic
 261 C makes copy of chains
 262         write (iout,*) "symetr", symetr
 263
 264       if (symetr.gt.1) then
 265        call permut(symetr)
 266        nperm=1
 267        do i=1,symetr
 268        nperm=nperm*i
 269        enddo
 270        do i=1,nperm
 271        write(iout,*) (tabperm(i,kkk),kkk=1,4)
 272        enddo
 273        do i=1,nperm
 274         cou=0
 275         do kkk=1,symetr
 276          icha=tabperm(i,kkk)
 277 c         write (iout,*) i,icha
 278          do lll=1,chain_length
 279           cou=cou+1
 280            if (cou.le.nres) then
 281            do j=1,3
 282             kupa=mod(lll,chain_length)
 283             iprzes=(kkk-1)*chain_length+lll
 284             if (kupa.eq.0) kupa=chain_length
 285 c            write (iout,*) "kupa", kupa
 286             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 287             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 288           enddo
 289           endif
 290          enddo
 291         enddo
 292        enddo
 293        endif
 294 C-koniec robienia kopii
 295 c diag
 296       do kkk=1,nperm
 297       write (iout,*) "nowa struktura", nperm
 298       do i=1,nres
 299         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 300      &cref(2,i,kkk),
 301      &cref(3,i,kkk),cref(1,nres+i,kkk),
 302      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 303       enddo
 304   100 format (//'              alpha-carbon coordinates       ',
 305      &          '     centroid coordinates'/
 306      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 307      &                          10X,'X',11X,'Y',11X,'Z')
 308   110 format (a,'(',i3,')',6f12.5)
 309
 310       enddo
 311 cc enddiag
 312       do j=1,nbfrag
 313         do i=1,4
 314          bfrag(i,j)=bfrag(i,j)-ishift
 315         enddo
 316       enddo
 317
 318       do j=1,nhfrag
 319         do i=1,2
 320          hfrag(i,j)=hfrag(i,j)-ishift
 321         enddo
 322       enddo
 323
 324       return
 325       end
 326 c---------------------------------------------------------------------------
 327       subroutine int_from_cart(lside,lprn)
 328       implicit real*8 (a-h,o-z)
 329       include 'DIMENSIONS'
 330 #ifdef MPI
 331       include "mpif.h"
 332 #endif
 333       include 'COMMON.LOCAL'
 334       include 'COMMON.VAR'
 335       include 'COMMON.CHAIN'
 336       include 'COMMON.INTERACT'
 337       include 'COMMON.IOUNITS'
 338       include 'COMMON.GEO'
 339       include 'COMMON.NAMES'
 340       include 'COMMON.CONTROL'
 341       include 'COMMON.SETUP'
 342       character*3 seq,atom,res
 343       character*80 card
 344       dimension sccor(3,20)
 345       integer rescode
 346       logical lside,lprn
 347 #ifdef MPI
 348       if(me.eq.king.or..not.out1file)then
 349 #endif
 350        if (lprn) then
 351         write (iout,'(/a)')
 352      &  'Internal coordinates calculated from crystal structure.'
 353         if (lside) then
 354           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 355      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 356      & '     Omega'
 357         else
 358           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 359      & '       Phi'
 360         endif
 361        endif
 362 #ifdef MPI
 363       endif
 364 #endif
 365       do i=1,nres-1
 366         iti=itype(i)
 367         if (iti.ne.21 .and. itype(i+1).ne.21 .and.
 368      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 369           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 370 ctest          stop
 371         endif
 372         vbld(i+1)=dist(i,i+1)
 373         vbld_inv(i+1)=1.0d0/vbld(i+1)
 374         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 375         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 376       enddo
 377 c      if (unres_pdb) then
 378 c        if (itype(1).eq.21) then
 379 c          theta(3)=90.0d0*deg2rad
 380 c          phi(4)=180.0d0*deg2rad
 381 c          vbld(2)=3.8d0
 382 c          vbld_inv(2)=1.0d0/vbld(2)
 383 c        endif
 384 c        if (itype(nres).eq.21) then
 385 c          theta(nres)=90.0d0*deg2rad
 386 c          phi(nres)=180.0d0*deg2rad
 387 c          vbld(nres)=3.8d0
 388 c          vbld_inv(nres)=1.0d0/vbld(2)
 389 c        endif
 390 c      endif
 391       if (lside) then
 392         do i=2,nres-1
 393           do j=1,3
 394             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 395      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 396           enddo
 397           iti=itype(i)
 398           di=dist(i,nres+i)
 399           vbld(i+nres)=di
 400           if (itype(i).ne.10) then
 401             vbld_inv(i+nres)=1.0d0/di
 402           else
 403             vbld_inv(i+nres)=0.0d0
 404           endif
 405           if (iti.ne.10) then
 406             alph(i)=alpha(nres+i,i,maxres2)
 407             omeg(i)=beta(nres+i,i,maxres2,i+1)
 408           endif
 409           if(me.eq.king.or..not.out1file)then
 410            if (lprn)
 411      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 412      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 413      &     rad2deg*alph(i),rad2deg*omeg(i)
 414           endif
 415         enddo
 416       else if (lprn) then
 417         do i=2,nres
 418           iti=itype(i)
 419           if(me.eq.king.or..not.out1file)
 420      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 421      &     rad2deg*theta(i),rad2deg*phi(i)
 422         enddo
 423       endif
 424       return
 425       end
 426 c-------------------------------------------------------------------------------
 427       subroutine sc_loc_geom(lprn)
 428       implicit real*8 (a-h,o-z)
 429       include 'DIMENSIONS'
 430 #ifdef MPI
 431       include "mpif.h"
 432 #endif
 433       include 'COMMON.LOCAL'
 434       include 'COMMON.VAR'
 435       include 'COMMON.CHAIN'
 436       include 'COMMON.INTERACT'
 437       include 'COMMON.IOUNITS'
 438       include 'COMMON.GEO'
 439       include 'COMMON.NAMES'
 440       include 'COMMON.CONTROL'
 441       include 'COMMON.SETUP'
 442       double precision x_prime(3),y_prime(3),z_prime(3)
 443       logical lprn
 444       do i=1,nres-1
 445         do j=1,3
 446           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 447         enddo
 448       enddo
 449       do i=2,nres-1
 450         if (itype(i).ne.10 .and. itype(i).ne.21) then
 451           do j=1,3
 452             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 453           enddo
 454         else
 455           do j=1,3
 456             dc_norm(j,i+nres)=0.0d0
 457           enddo
 458         endif
 459       enddo
 460       do i=2,nres-1
 461         costtab(i+1) =dcos(theta(i+1))
 462         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 463         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 464         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 465         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 466         cosfac=dsqrt(cosfac2)
 467         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 468         sinfac=dsqrt(sinfac2)
 469         it=itype(i)
 470         if (it.ne.10 .and. itype(i).ne.21) then
 471 c
 472 C  Compute the axes of tghe local cartesian coordinates system; store in
 473 c   x_prime, y_prime and z_prime
 474 c
 475         do j=1,3
 476           x_prime(j) = 0.00
 477           y_prime(j) = 0.00
 478           z_prime(j) = 0.00
 479         enddo
 480         do j = 1,3
 481           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 482           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 483         enddo
 484         call vecpr(x_prime,y_prime,z_prime)
 485 c
 486 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 487 C to local coordinate system. Store in xx, yy, zz.
 488 c
 489         xx=0.0d0
 490         yy=0.0d0
 491         zz=0.0d0
 492         do j = 1,3
 493           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 494           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 495           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 496         enddo
 497
 498         xxref(i)=xx
 499         yyref(i)=yy
 500         zzref(i)=zz
 501         else
 502         xxref(i)=0.0d0
 503         yyref(i)=0.0d0
 504         zzref(i)=0.0d0
 505         endif
 506       enddo
 507       if (lprn) then
 508         do i=2,nres
 509           iti=itype(i)
 510 #ifdef MPI
 511           if(me.eq.king.or..not.out1file)
 512      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 513      &      yyref(i),zzref(i)
 514 #else
 515           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
 516      &     zzref(i)
 517 #endif
 518         enddo
 519       endif
 520       return
 521       end
 522 c---------------------------------------------------------------------------
 523       subroutine sccenter(ires,nscat,sccor)
 524       implicit real*8 (a-h,o-z)
 525       include 'DIMENSIONS'
 526       include 'COMMON.CHAIN'
 527       dimension sccor(3,20)
 528       do j=1,3
 529         sccmj=0.0D0
 530         do i=1,nscat
 531           sccmj=sccmj+sccor(j,i)
 532         enddo
 533         dc(j,ires)=sccmj/nscat
 534       enddo
 535       return
 536       end
 537 c---------------------------------------------------------------------------
 538       subroutine bond_regular
 539       implicit real*8 (a-h,o-z)
 540       include 'DIMENSIONS'
 541       include 'COMMON.VAR'
 542       include 'COMMON.LOCAL'
 543       include 'COMMON.CALC'
 544       include 'COMMON.INTERACT'
 545       include 'COMMON.CHAIN'
 546       do i=1,nres-1
 547        vbld(i+1)=vbl
 548        vbld_inv(i+1)=1.0d0/vbld(i+1)
 549        vbld(i+1+nres)=dsc(itype(i+1))
 550        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 551 c       print *,vbld(i+1),vbld(i+1+nres)
 552       enddo
 553       return
 554       end
 555