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