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(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       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       return
 373       end
 374 c---------------------------------------------------------------------------
 375       subroutine int_from_cart(lside,lprn)
 376       implicit real*8 (a-h,o-z)
 377       include 'DIMENSIONS'
 378 #ifdef MPI
 379       include "mpif.h"
 380 #endif
 381       include 'COMMON.LOCAL'
 382       include 'COMMON.VAR'
 383       include 'COMMON.CHAIN'
 384       include 'COMMON.INTERACT'
 385       include 'COMMON.IOUNITS'
 386       include 'COMMON.GEO'
 387       include 'COMMON.NAMES'
 388       include 'COMMON.CONTROL'
 389       include 'COMMON.SETUP'
 390       character*3 seq,atom,res
 391       character*80 card
 392       dimension sccor(3,20)
 393       integer rescode
 394       logical lside,lprn
 395 #ifdef MPI
 396       if(me.eq.king.or..not.out1file)then
 397 #endif
 398        if (lprn) then
 399         write (iout,'(/a)')
 400      &  'Internal coordinates calculated from crystal structure.'
 401         if (lside) then
 402           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 403      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 404      & '     Omega'
 405         else
 406           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 407      & '       Phi'
 408         endif
 409        endif
 410 #ifdef MPI
 411       endif
 412 #endif
 413       do i=1,nres-1
 414         iti=itype(i)
 415         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 416      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 417           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 418 ctest          stop
 419         endif
 420         vbld(i+1)=dist(i,i+1)
 421         vbld_inv(i+1)=1.0d0/vbld(i+1)
 422         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 423         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 424       enddo
 425 c      if (unres_pdb) then
 426 c        if (itype(1).eq.21) then
 427 c          theta(3)=90.0d0*deg2rad
 428 c          phi(4)=180.0d0*deg2rad
 429 c          vbld(2)=3.8d0
 430 c          vbld_inv(2)=1.0d0/vbld(2)
 431 c        endif
 432 c        if (itype(nres).eq.21) then
 433 c          theta(nres)=90.0d0*deg2rad
 434 c          phi(nres)=180.0d0*deg2rad
 435 c          vbld(nres)=3.8d0
 436 c          vbld_inv(nres)=1.0d0/vbld(2)
 437 c        endif
 438 c      endif
 439       if (lside) then
 440         do i=2,nres-1
 441           do j=1,3
 442             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 443      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 444           enddo
 445           iti=itype(i)
 446           di=dist(i,nres+i)
 447           vbld(i+nres)=di
 448           if (itype(i).ne.10) then
 449             vbld_inv(i+nres)=1.0d0/di
 450           else
 451             vbld_inv(i+nres)=0.0d0
 452           endif
 453           if (iti.ne.10) then
 454             alph(i)=alpha(nres+i,i,maxres2)
 455             omeg(i)=beta(nres+i,i,maxres2,i+1)
 456           endif
 457           if(me.eq.king.or..not.out1file)then
 458            if (lprn)
 459      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 460      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 461      &     rad2deg*alph(i),rad2deg*omeg(i)
 462           endif
 463         enddo
 464       else if (lprn) then
 465         do i=2,nres
 466           iti=itype(i)
 467           if(me.eq.king.or..not.out1file)
 468      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 469      &     rad2deg*theta(i),rad2deg*phi(i)
 470         enddo
 471       endif
 472       return
 473       end
 474 c-------------------------------------------------------------------------------
 475       subroutine sc_loc_geom(lprn)
 476       implicit real*8 (a-h,o-z)
 477       include 'DIMENSIONS'
 478 #ifdef MPI
 479       include "mpif.h"
 480 #endif
 481       include 'COMMON.LOCAL'
 482       include 'COMMON.VAR'
 483       include 'COMMON.CHAIN'
 484       include 'COMMON.INTERACT'
 485       include 'COMMON.IOUNITS'
 486       include 'COMMON.GEO'
 487       include 'COMMON.NAMES'
 488       include 'COMMON.CONTROL'
 489       include 'COMMON.SETUP'
 490       double precision x_prime(3),y_prime(3),z_prime(3)
 491       logical lprn
 492       do i=1,nres-1
 493         do j=1,3
 494           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 495         enddo
 496       enddo
 497       do i=2,nres-1
 498         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 499           do j=1,3
 500             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 501           enddo
 502         else
 503           do j=1,3
 504             dc_norm(j,i+nres)=0.0d0
 505           enddo
 506         endif
 507       enddo
 508       do i=2,nres-1
 509         costtab(i+1) =dcos(theta(i+1))
 510         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 511         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 512         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 513         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 514         cosfac=dsqrt(cosfac2)
 515         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 516         sinfac=dsqrt(sinfac2)
 517         it=itype(i)
 518         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 519 c
 520 C  Compute the axes of tghe local cartesian coordinates system; store in
 521 c   x_prime, y_prime and z_prime
 522 c
 523         do j=1,3
 524           x_prime(j) = 0.00
 525           y_prime(j) = 0.00
 526           z_prime(j) = 0.00
 527         enddo
 528         do j = 1,3
 529           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 530           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 531         enddo
 532         call vecpr(x_prime,y_prime,z_prime)
 533 c
 534 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 535 C to local coordinate system. Store in xx, yy, zz.
 536 c
 537         xx=0.0d0
 538         yy=0.0d0
 539         zz=0.0d0
 540         do j = 1,3
 541           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 542           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 543           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 544         enddo
 545
 546         xxref(i)=xx
 547         yyref(i)=yy
 548         zzref(i)=zz
 549         else
 550         xxref(i)=0.0d0
 551         yyref(i)=0.0d0
 552         zzref(i)=0.0d0
 553         endif
 554       enddo
 555       if (lprn) then
 556         do i=2,nres
 557           iti=itype(i)
 558 #ifdef MPI
 559           if(me.eq.king.or..not.out1file)
 560      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 561      &      yyref(i),zzref(i)
 562 #else
 563           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
 564      &     zzref(i)
 565 #endif
 566         enddo
 567       endif
 568       return
 569       end
 570 c---------------------------------------------------------------------------
 571       subroutine sccenter(ires,nscat,sccor)
 572       implicit real*8 (a-h,o-z)
 573       include 'DIMENSIONS'
 574       include 'COMMON.CHAIN'
 575       dimension sccor(3,20)
 576       do j=1,3
 577         sccmj=0.0D0
 578         do i=1,nscat
 579           sccmj=sccmj+sccor(j,i)
 580         enddo
 581         dc(j,ires)=sccmj/nscat
 582       enddo
 583       return
 584       end
 585 c---------------------------------------------------------------------------
 586       subroutine bond_regular
 587       implicit real*8 (a-h,o-z)
 588       include 'DIMENSIONS'
 589       include 'COMMON.VAR'
 590       include 'COMMON.LOCAL'
 591       include 'COMMON.CALC'
 592       include 'COMMON.INTERACT'
 593       include 'COMMON.CHAIN'
 594       do i=1,nres-1
 595        vbld(i+1)=vbl
 596        vbld_inv(i+1)=1.0d0/vbld(i+1)
 597        vbld(i+1+nres)=dsc(iabs(itype(i+1)))
 598        vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
 599 c       print *,vbld(i+1),vbld(i+1+nres)
 600       enddo
 601       return
 602       end
 603