source/unres/src-HCD-5D/readpdb.F

   1       subroutine readpdb
   2 C Read the PDB file and convert the peptide geometry into virtual-chain
   3 C geometry.
   4       implicit none
   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.FRAG'
  15       include 'COMMON.SETUP'
  16       include 'COMMON.SBRIDGE'
  17       character*3 seq,atom,res
  18       character*80 card
  19       double precision sccor(3,50)
  20       double precision e1(3),e2(3),e3(3)
  21       integer rescode,iterter(maxres),cou
  22       logical fail,sccalc
  23       integer i,j,iii,ires,ires_old,ishift,ishift1,ibeg,ifree
  24       double precision dcj!,efree_temp
  25       logical zero
  26       bfac=0.0d0
  27       do i=1,maxres
  28          iterter(i)=0
  29       enddo
  30       ibeg=1
  31       ishift1=0
  32       lsecondary=.false.
  33       nhfrag=0
  34       nbfrag=0
  35       iii=0
  36       sccalc=.false.
  37       do
  38         read (ipdbin,'(a80)',end=10) card
  39 c        write (iout,'(a)') card
  40 c        call flush(iout)
  41         if (card(:5).eq.'HELIX') then
  42           nhfrag=nhfrag+1
  43           lsecondary=.true.
  44           read(card(22:25),*) hfrag(1,nhfrag)
  45           read(card(34:37),*) hfrag(2,nhfrag)
  46         endif
  47         if (card(:5).eq.'SHEET') then
  48           nbfrag=nbfrag+1
  49           lsecondary=.true.
  50           read(card(24:26),*) bfrag(1,nbfrag)
  51           read(card(35:37),*) bfrag(2,nbfrag)
  52 !rc----------------------------------------
  53 !rc  to be corrected !!!
  54           bfrag(3,nbfrag)=bfrag(1,nbfrag)
  55           bfrag(4,nbfrag)=bfrag(2,nbfrag)
  56 !rc----------------------------------------
  57         endif
  58         if (card(:3).eq.'END') then
  59           goto 10
  60         else if (card(:3).eq.'TER') then
  61 ! End current chain
  62           ires_old=ires+2
  63           itype(ires_old-1)=ntyp1
  64           iterter(ires_old-1)=1
  65           itype(ires_old)=ntyp1
  66           iterter(ires_old)=1
  67           ishift1=ishift1+1
  68           ibeg=2
  69           write (iout,*) "Chain ended",ires,ishift,ires_old
  70           if (unres_pdb) then
  71             do j=1,3
  72               dc(j,ires)=sccor(j,iii)
  73             enddo
  74           else
  75             call sccenter(ires,iii,sccor)
  76           endif
  77           iii=0
  78           sccalc=.true.
  79         endif
  80 ! Read free energy
  81 c        if (index(card,"FREE ENERGY").gt.0) then
  82 c          ifree=index(card,"FREE ENERGY")+12
  83 c          read(card(ifree:),*,err=1115,end=1115) efree_temp
  84 c 1115     continue
  85 c        endif
  86 ! Fish out the ATOM cards.
  87         if (index(card(1:4),'ATOM').gt.0) then
  88           sccalc=.false.
  89           read (card(12:16),*) atom
  90 c          write (2,'(a)') card
  91 c          write (iout,*) "ibeg",ibeg
  92 c          write (iout,*) "! ",atom," !",ires
  93 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
  94           read (card(23:26),*) ires
  95           read (card(18:20),'(a3)') res
  96 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
  97 c     &      " ires_old",ires_old
  98 c          write (iout,*) "ishift",ishift," ishift1",ishift1
  99 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 100           if (ires-ishift+ishift1.ne.ires_old) then
 101 ! Calculate the CM of the preceding residue.
 102 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
 103             if (ibeg.eq.0) then
 104 c              write (iout,*) "Calculating sidechain center iii",iii
 105 c              write (iout,*) "ires",ires
 106               if (unres_pdb) then
 107 c                write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3)
 108                 do j=1,3
 109                   dc(j,ires_old)=sccor(j,iii)
 110                 enddo
 111               else
 112                 call sccenter(ires_old,iii,sccor)
 113               endif
 114               iii=0
 115               sccalc=.true.
 116             endif
 117 ! Start new residue.
 118             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 119               ires=ires_old
 120               cycle
 121             else if (ibeg.eq.1) then
 122 c              write (iout,*) "BEG ires",ires
 123               ishift=ires-1
 124               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 125                 ishift=ishift-1
 126                 itype(1)=ntyp1
 127               endif
 128               ires=ires-ishift+ishift1
 129               ires_old=ires
 130 !              write (iout,*) "ishift",ishift," ires",ires,&
 131 !               " ires_old",ires_old
 132               ibeg=0
 133             else if (ibeg.eq.2) then
 134 ! Start a new chain
 135               ishift=-ires_old+ires-1 !!!!!
 136               ishift1=ishift1-1    !!!!!
 137 c              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
 138               ires=ires-ishift+ishift1
 139               ires_old=ires
 140               ibeg=0
 141             else
 142               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 143               ires=ires-ishift+ishift1
 144               ires_old=ires
 145             endif
 146             if (res.eq.'ACE' .or. res.eq.'NHE') then
 147               itype(ires)=10
 148             else
 149               itype(ires)=rescode(ires,res,0)
 150             endif
 151           else
 152             ires=ires-ishift+ishift1
 153           endif
 154 c          write (iout,*) "ires_old",ires_old," ires",ires
 155           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 156 !            ishift1=ishift1+1
 157           endif
 158 c          write (2,*) "ires",ires," res ",res!," ity"!,ity
 159           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 160      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 161             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 162             write (iout,*) "backbone ",atom
 163             write (iout,*) ires,res,(c(j,ires),j=1,3)
 164 #ifdef DEBUG
 165             write (iout,'(i6,i3,2x,a,3f8.3)')
 166      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 167 #endif
 168             iii=iii+1
 169             do j=1,3
 170               sccor(j,iii)=c(j,ires)
 171             enddo
 172 c            write (2,*) card(23:27),ires,itype(ires),iii
 173           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 174      &             atom.ne.'N' .and. atom.ne.'C' .and.
 175      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 176      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 177 !            write (iout,*) "sidechain ",atom
 178             iii=iii+1
 179             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 180 c            write (2,*) "iii",iii
 181           endif
 182         endif
 183       enddo
 184    10 if(me.eq.king.or..not.out1file)
 185      & write (iout,'(a,i5)') ' Nres: ',ires
 186 c      write (iout,*) "iii",iii
 187 C Calculate dummy residue coordinates inside the "chain" of a multichain
 188 C system
 189       nres=ires
 190 c      write (iout,*) "dc"
 191 c      do i=1,nres
 192 c        write (iout,'(i5,3f10.5)') i,(dc(j,i),j=1,3)
 193 c      enddo
 194       do i=2,nres-1
 195 c        write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
 196         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 197          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 198 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 199 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 200 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 201            if (unres_pdb) then
 202 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 203 c            print *,i,'tu dochodze'
 204             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 205             if (fail) then
 206               e2(1)=0.0d0
 207               e2(2)=1.0d0
 208               e2(3)=0.0d0
 209             endif !fail
 210 c            print *,i,'a tu?'
 211             do j=1,3
 212              c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 213             enddo
 214            else   !unres_pdb
 215            do j=1,3
 216              dcj=(c(j,i-2)-c(j,i-3))/2.0
 217             if (dcj.eq.0) dcj=1.23591524223
 218              c(j,i)=c(j,i-1)+dcj
 219              c(j,nres+i)=c(j,i)
 220              dC(j,i)=c(j,i)
 221            enddo
 222           endif   !unres_pdb
 223          else     !itype(i+1).eq.ntyp1
 224           if (unres_pdb) then
 225 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 226             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 227             if (fail) then
 228               e2(1)=0.0d0
 229               e2(2)=1.0d0
 230               e2(3)=0.0d0
 231             endif
 232             do j=1,3
 233               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 234             enddo
 235           else !unres_pdb
 236            do j=1,3
 237             dcj=(c(j,i+3)-c(j,i+2))/2.0
 238             if (dcj.eq.0) dcj=1.23591524223
 239             c(j,i)=c(j,i+1)-dcj
 240             c(j,nres+i)=c(j,i)
 241             dC(j,i)=c(j,i)
 242            enddo
 243           endif !unres_pdb
 244          endif !itype(i+1).eq.ntyp1
 245         endif  !itype.eq.ntyp1
 246       enddo
 247       write (iout,*) "After loop in readpbd"
 248 C Calculate the CM of the last side chain.
 249       if (.not. sccalc) then
 250       if (unres_pdb) then
 251         do j=1,3
 252           dc(j,ires)=sccor(j,iii)
 253         enddo
 254       else
 255 c        write (iout,*) "Calling sccenter iii",iii
 256         call sccenter(ires,iii,sccor)
 257       endif
 258       endif
 259       nsup=nres
 260       nstart_sup=1
 261       if (itype(nres).ne.10) then
 262         nres=nres+1
 263         itype(nres)=ntyp1
 264         if (unres_pdb) then
 265 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 266           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 267           if (fail) then
 268             e2(1)=0.0d0
 269             e2(2)=1.0d0
 270             e2(3)=0.0d0
 271           endif
 272           do j=1,3
 273             c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 274           enddo
 275         else
 276         do j=1,3
 277           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 278             if (dcj.eq.0) dcj=1.23591524223
 279           c(j,nres)=c(j,nres-1)+dcj
 280           c(j,2*nres)=c(j,nres)
 281         enddo
 282         endif
 283       endif
 284       do i=2,nres-1
 285         do j=1,3
 286           c(j,i+nres)=dc(j,i)
 287         enddo
 288       enddo
 289       do j=1,3
 290         c(j,nres+1)=c(j,1)
 291         c(j,2*nres)=c(j,nres)
 292       enddo
 293       if (itype(1).eq.ntyp1) then
 294         nsup=nsup-1
 295         nstart_sup=2
 296         if (unres_pdb) then
 297 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 298           call refsys(2,3,4,e1,e2,e3,fail)
 299           if (fail) then
 300             e2(1)=0.0d0
 301             e2(2)=1.0d0
 302             e2(3)=0.0d0
 303           endif
 304           do j=1,3
 305             c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
 306           enddo
 307         else
 308         do j=1,3
 309           dcj=(c(j,4)-c(j,3))/2.0
 310           c(j,1)=c(j,2)-dcj
 311           c(j,nres+1)=c(j,1)
 312         enddo
 313         endif
 314       endif
 315 C Calculate internal coordinates.
 316       if(me.eq.king.or..not.out1file)then
 317       write (iout,'(/a)')
 318      &  "Cartesian coordinates of the reference structure"
 319       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 320      & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 321       do ires=1,nres
 322         write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')
 323      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 324      &    (c(j,ires+nres),j=1,3)
 325       enddo
 326       call flush(iout)
 327       endif
 328       zero=.false.
 329       do ires=1,nres
 330         zero=zero.or.itype(ires).eq.0
 331       enddo
 332       if (zero) then
 333         write (iout,'(2a)') "Gaps in PDB coordinates detected;",
 334      &    " look for ZERO in the control output above."
 335         write (iout,'(2a)') "Repair the PDB file using MODELLER",
 336      &    " or other softwared and resubmit."
 337         call flush(iout)
 338         stop
 339       endif
 340 c      write(iout,*)"before int_from_cart nres",nres
 341       call int_from_cart(.true.,.false.)
 342       do i=1,nres
 343         thetaref(i)=theta(i)
 344         phiref(i)=phi(i)
 345       enddo
 346       dc(:,0)=c(:,1)
 347       do i=1,nres-1
 348         do j=1,3
 349           dc(j,i)=c(j,i+1)-c(j,i)
 350           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 351         enddo
 352 c        write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
 353 c     &   vbld_inv(i+1)
 354       enddo
 355       do i=2,nres-1
 356         do j=1,3
 357           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 358           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 359         enddo
 360 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 361 c     &   vbld_inv(i+nres)
 362       enddo
 363       call sc_loc_geom(.false.)
 364       call int_from_cart1(.false.)
 365 c      call chainbuild
 366 C Copy the coordinates to reference coordinates
 367       do i=1,nres
 368         do j=1,3
 369           cref(j,i)=c(j,i)
 370           cref(j,i+nres)=c(j,i+nres)
 371         enddo
 372       enddo
 373   100 format (//'              alpha-carbon coordinates       ',
 374      &          '     centroid coordinates'/
 375      1          '       ', 7X,'X',11X,'Y',11X,'Z',
 376      &                          10X,'X',11X,'Y',11X,'Z')
 377   110 format (a,'(',i4,')',6f12.5)
 378 cc enddiag
 379       do j=1,nbfrag
 380         do i=1,4
 381          bfrag(i,j)=bfrag(i,j)-ishift
 382         enddo
 383       enddo
 384
 385       do j=1,nhfrag
 386         do i=1,2
 387          hfrag(i,j)=hfrag(i,j)-ishift
 388         enddo
 389       enddo
 390       return
 391       end
 392