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