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