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