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