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