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