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