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