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
  23       integer i,j,iii,ires,ires_old,ishift,ibeg
  24       double precision dcj
  25       bfac=0.0d0
  26       do i=1,maxres
  27          iterter(i)=0
  28       enddo
  29       ibeg=1
  30       lsecondary=.false.
  31       nhfrag=0
  32       nbfrag=0
  33       do
  34         read (ipdbin,'(a80)',end=10) card
  35         if (card(:5).eq.'HELIX') then
  36          nhfrag=nhfrag+1
  37          lsecondary=.true.
  38          read(card(22:25),*) hfrag(1,nhfrag)
  39          read(card(34:37),*) hfrag(2,nhfrag)
  40         endif
  41         if (card(:5).eq.'SHEET') then
  42          nbfrag=nbfrag+1
  43          lsecondary=.true.
  44          read(card(24:26),*) bfrag(1,nbfrag)
  45          read(card(35:37),*) bfrag(2,nbfrag)
  46 crc----------------------------------------
  47 crc  to be corrected !!!
  48          bfrag(3,nbfrag)=bfrag(1,nbfrag)
  49          bfrag(4,nbfrag)=bfrag(2,nbfrag)
  50 crc----------------------------------------
  51         endif
  52         if (card(:3).eq.'END') then
  53           goto 10
  54         else if (card(:3).eq.'TER') then
  55 C End current chain
  56           ires_old=ires+2
  57           itype(ires_old-1)=ntyp1
  58           iterter(ires_old-1)=1
  59           itype(ires_old)=ntyp1
  60           iterter(ires_old)=1
  61           ibeg=2
  62           write (iout,*) "Chain ended",ires,ishift,ires_old
  63           if (unres_pdb) then
  64             do j=1,3
  65               dc(j,ires)=sccor(j,iii)
  66             enddo
  67           else
  68             call sccenter(ires,iii,sccor)
  69           endif
  70         endif
  71 C Fish out the ATOM cards.
  72         if (index(card(1:4),'ATOM').gt.0) then
  73           read (card(14:16),'(a3)') atom
  74           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  75 C Calculate the CM of the preceding residue.
  76             if (ibeg.eq.0) then
  77               if (unres_pdb) then
  78                 do j=1,3
  79                   dc(j,ires+nres)=sccor(j,iii)
  80                 enddo
  81               else
  82                 call sccenter(ires,iii,sccor)
  83               endif
  84             endif
  85 C Start new residue.
  86 c            write (iout,'(a80)') card
  87             read (card(23:26),*) ires
  88             read (card(18:20),'(a3)') res
  89             if (ibeg.eq.1) then
  90               ishift=ires-1
  91               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  92                 ishift=ishift-1
  93                 itype(1)=ntyp1
  94               endif
  95 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
  96               ibeg=0
  97             else if (ibeg.eq.2) then
  98 c Start a new chain
  99               ishift=-ires_old+ires-1
 100 c              write (iout,*) "New chain started",ires,ishift
 101               ibeg=0
 102             endif
 103             ires=ires-ishift
 104 c            write (2,*) "ires",ires," ishift",ishift
 105             if (res.eq.'ACE') then
 106               itype(ires)=10
 107             else
 108               itype(ires)=rescode(ires,res,0)
 109             endif
 110             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 111             read(card(61:66),*) bfac(ires)
 112 c            if(me.eq.king.or..not.out1file)
 113 c     &       write (iout,'(2i3,2x,a,3f8.3)')
 114 c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
 115             iii=1
 116             do j=1,3
 117               sccor(j,iii)=c(j,ires)
 118             enddo
 119           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
 120      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
 121             iii=iii+1
 122             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 123           endif
 124         endif
 125       enddo
 126    10 if(me.eq.king.or..not.out1file)
 127      & write (iout,'(a,i5)') ' Nres: ',ires
 128 C Calculate dummy residue coordinates inside the "chain" of a multichain
 129 C system
 130       nres=ires
 131       do i=2,nres-1
 132 c        write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
 133         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 134          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 135 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 136 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 137 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 138            if (unres_pdb) then
 139 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 140             print *,i,'tu dochodze'
 141             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 142             if (fail) then
 143               e2(1)=0.0d0
 144               e2(2)=1.0d0
 145               e2(3)=0.0d0
 146             endif !fail
 147             print *,i,'a tu?'
 148             do j=1,3
 149              c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 150             enddo
 151            else   !unres_pdb
 152            do j=1,3
 153              dcj=(c(j,i-2)-c(j,i-3))/2.0
 154             if (dcj.eq.0) dcj=1.23591524223
 155              c(j,i)=c(j,i-1)+dcj
 156              c(j,nres+i)=c(j,i)
 157            enddo
 158           endif   !unres_pdb
 159          else     !itype(i+1).eq.ntyp1
 160           if (unres_pdb) then
 161 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 162             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 163             if (fail) then
 164               e2(1)=0.0d0
 165               e2(2)=1.0d0
 166               e2(3)=0.0d0
 167             endif
 168             do j=1,3
 169               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 170             enddo
 171           else !unres_pdb
 172            do j=1,3
 173             dcj=(c(j,i+3)-c(j,i+2))/2.0
 174             if (dcj.eq.0) dcj=1.23591524223
 175             c(j,i)=c(j,i+1)-dcj
 176             c(j,nres+i)=c(j,i)
 177            enddo
 178           endif !unres_pdb
 179          endif !itype(i+1).eq.ntyp1
 180         endif  !itype.eq.ntyp1
 181       enddo
 182       write (iout,*) "After loop in readpbd"
 183 C Calculate the CM of the last side chain.
 184       if (unres_pdb) then
 185         do j=1,3
 186           dc(j,ires)=sccor(j,iii)
 187         enddo
 188       else
 189         call sccenter(ires,iii,sccor)
 190       endif
 191       nsup=nres
 192       nstart_sup=1
 193       if (itype(nres).ne.10) then
 194         nres=nres+1
 195         itype(nres)=ntyp1
 196         if (unres_pdb) then
 197 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 198           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 199           if (fail) then
 200             e2(1)=0.0d0
 201             e2(2)=1.0d0
 202             e2(3)=0.0d0
 203           endif
 204           do j=1,3
 205             c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 206           enddo
 207         else
 208         do j=1,3
 209           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 210             if (dcj.eq.0) dcj=1.23591524223
 211           c(j,nres)=c(j,nres-1)+dcj
 212           c(j,2*nres)=c(j,nres)
 213         enddo
 214         endif
 215       endif
 216       do i=2,nres-1
 217         do j=1,3
 218           c(j,i+nres)=dc(j,i)
 219         enddo
 220       enddo
 221       do j=1,3
 222         c(j,nres+1)=c(j,1)
 223         c(j,2*nres)=c(j,nres)
 224       enddo
 225       if (itype(1).eq.ntyp1) then
 226         nsup=nsup-1
 227         nstart_sup=2
 228         if (unres_pdb) then
 229 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 230           call refsys(2,3,4,e1,e2,e3,fail)
 231           if (fail) then
 232             e2(1)=0.0d0
 233             e2(2)=1.0d0
 234             e2(3)=0.0d0
 235           endif
 236           do j=1,3
 237             c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
 238           enddo
 239         else
 240         do j=1,3
 241           dcj=(c(j,4)-c(j,3))/2.0
 242           c(j,1)=c(j,2)-dcj
 243           c(j,nres+1)=c(j,1)
 244         enddo
 245         endif
 246       endif
 247 C Calculate internal coordinates.
 248       if(me.eq.king.or..not.out1file)then
 249        do ires=1,nres
 250         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 251      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 252      &    (c(j,nres+ires),j=1,3)
 253        enddo
 254       endif
 255       call flush(iout)
 256 c      write(iout,*)"before int_from_cart nres",nres
 257       call int_from_cart(.true.,.false.)
 258       do i=1,nres
 259         thetaref(i)=theta(i)
 260         phiref(i)=phi(i)
 261       enddo
 262       do i=1,nres-1
 263         do j=1,3
 264           dc(j,i)=c(j,i+1)-c(j,i)
 265           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 266         enddo
 267 c        write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
 268 c     &   vbld_inv(i+1)
 269       enddo
 270       do i=2,nres-1
 271         do j=1,3
 272           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 273           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 274         enddo
 275 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 276 c     &   vbld_inv(i+nres)
 277       enddo
 278       call sc_loc_geom(.false.)
 279       call int_from_cart1(.false.)
 280 c      call chainbuild
 281 C Copy the coordinates to reference coordinates
 282       do i=1,nres
 283         do j=1,3
 284           cref(j,i)=c(j,i)
 285           cref(j,i+nres)=c(j,i+nres)
 286         enddo
 287       enddo
 288   100 format (//'              alpha-carbon coordinates       ',
 289      &          '     centroid coordinates'/
 290      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 291      &                          10X,'X',11X,'Y',11X,'Z')
 292   110 format (a,'(',i3,')',6f12.5)
 293 cc enddiag
 294       do j=1,nbfrag
 295         do i=1,4
 296          bfrag(i,j)=bfrag(i,j)-ishift
 297         enddo
 298       enddo
 299
 300       do j=1,nhfrag
 301         do i=1,2
 302          hfrag(i,j)=hfrag(i,j)-ishift
 303         enddo
 304       enddo
 305       return
 306       end
 307 c---------------------------------------------------------------------------
 308       subroutine int_from_cart(lside,lprn)
 309       implicit none
 310       include 'DIMENSIONS'
 311 #ifdef MPI
 312       include "mpif.h"
 313 #endif
 314       include 'COMMON.LOCAL'
 315       include 'COMMON.VAR'
 316       include 'COMMON.CHAIN'
 317       include 'COMMON.INTERACT'
 318       include 'COMMON.IOUNITS'
 319       include 'COMMON.GEO'
 320       include 'COMMON.NAMES'
 321       include 'COMMON.CONTROL'
 322       include 'COMMON.SETUP'
 323       double precision dist,alpha,beta
 324       character*3 seq,atom,res
 325       character*80 card
 326       double precision sccor(3,50)
 327       integer rescode
 328       logical lside,lprn
 329       integer i,j,iti
 330       double precision di,cosfac2,sinfac2,cosfac,sinfac
 331 #ifdef MPI
 332       if(me.eq.king.or..not.out1file)then
 333 #endif
 334        if (lprn) then
 335         write (iout,'(/a)')
 336      &  'Internal coordinates calculated from crystal structure.'
 337         if (lside) then
 338           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 339      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 340      & '     Omega'
 341         else
 342           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 343      & '       Phi'
 344         endif
 345        endif
 346 #ifdef MPI
 347       endif
 348 #endif
 349       do i=1,nres-1
 350         iti=itype(i)
 351         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 352      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 353           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 354 ctest          stop
 355         endif
 356         vbld(i+1)=dist(i,i+1)
 357         vbld_inv(i+1)=1.0d0/vbld(i+1)
 358 c        write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
 359 c     &      vbld_inv(i+1)
 360         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 361         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 362       enddo
 363 c      if (unres_pdb) then
 364 c        if (itype(1).eq.21) then
 365 c          theta(3)=90.0d0*deg2rad
 366 c          phi(4)=180.0d0*deg2rad
 367 c          vbld(2)=3.8d0
 368 c          vbld_inv(2)=1.0d0/vbld(2)
 369 c        endif
 370 c        if (itype(nres).eq.21) then
 371 c          theta(nres)=90.0d0*deg2rad
 372 c          phi(nres)=180.0d0*deg2rad
 373 c          vbld(nres)=3.8d0
 374 c          vbld_inv(nres)=1.0d0/vbld(2)
 375 c        endif
 376 c      endif
 377 c      print *,"A TU2"
 378       if (lside) then
 379         do i=2,nres-1
 380           do j=1,3
 381             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 382      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 383           enddo
 384           iti=itype(i)
 385           di=dist(i,nres+i)
 386           vbld(i+nres)=di
 387           if (itype(i).ne.10) then
 388             vbld_inv(i+nres)=1.0d0/di
 389           else
 390             vbld_inv(i+nres)=0.0d0
 391           endif
 392           if (iti.ne.10) then
 393             alph(i)=alpha(nres+i,i,maxres2)
 394             omeg(i)=beta(nres+i,i,maxres2,i+1)
 395           endif
 396           if(me.eq.king.or..not.out1file)then
 397            if (lprn)
 398      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 399      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 400      &     rad2deg*alph(i),rad2deg*omeg(i)
 401           endif
 402         enddo
 403         if (lprn) then
 404            i=nres
 405            iti=itype(nres)
 406            write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 407      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 408      &     rad2deg*alph(i),rad2deg*omeg(i)
 409         endif
 410       else if (lprn) then
 411         do i=2,nres
 412           iti=itype(i)
 413           if(me.eq.king.or..not.out1file)
 414      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 415      &     rad2deg*theta(i),rad2deg*phi(i)
 416         enddo
 417       endif
 418       return
 419       end
 420 c-------------------------------------------------------------------------------
 421       subroutine sc_loc_geom(lprn)
 422       implicit none
 423       include 'DIMENSIONS'
 424 #ifdef MPI
 425       include "mpif.h"
 426 #endif
 427       include 'COMMON.LOCAL'
 428       include 'COMMON.VAR'
 429       include 'COMMON.CHAIN'
 430       include 'COMMON.INTERACT'
 431       include 'COMMON.IOUNITS'
 432       include 'COMMON.GEO'
 433       include 'COMMON.NAMES'
 434       include 'COMMON.CONTROL'
 435       include 'COMMON.SETUP'
 436       double precision x_prime(3),y_prime(3),z_prime(3)
 437       logical lprn
 438       integer i,j,it
 439       double precision xx,yy,zz,cosfac,cosfac2,sinfac,sinfac2
 440       do i=1,nres-1
 441         do j=1,3
 442           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 443         enddo
 444 c        write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
 445 c     &    " vbld",vbld_inv(i+1)
 446       enddo
 447       do i=2,nres-1
 448         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 449           do j=1,3
 450             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 451           enddo
 452 c          write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
 453 c     &    " vbld",vbld_inv(i+nres)
 454         else
 455           do j=1,3
 456             dc_norm(j,i+nres)=0.0d0
 457           enddo
 458         endif
 459       enddo
 460       do i=2,nres-1
 461         costtab(i+1) =dcos(theta(i+1))
 462         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 463         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 464         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 465         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 466         cosfac=dsqrt(cosfac2)
 467         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 468         sinfac=dsqrt(sinfac2)
 469         it=itype(i)
 470 c        write (iout,*) "i",i," costab",costtab(i+1),
 471 c     &    " sintab",sinttab(i+1)
 472 c        write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
 473 c        write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
 474         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 475 c
 476 C  Compute the axes of tghe local cartesian coordinates system; store in
 477 c   x_prime, y_prime and z_prime
 478 c
 479         do j=1,3
 480           x_prime(j) = 0.00
 481           y_prime(j) = 0.00
 482           z_prime(j) = 0.00
 483         enddo
 484         do j = 1,3
 485           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 486           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 487         enddo
 488 c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
 489 c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
 490         call vecpr(x_prime,y_prime,z_prime)
 491 c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
 492 c
 493 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 494 C to local coordinate system. Store in xx, yy, zz.
 495 c
 496         xx=0.0d0
 497         yy=0.0d0
 498         zz=0.0d0
 499         do j = 1,3
 500           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 501           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 502           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 503         enddo
 504
 505         xxref(i)=xx
 506         yyref(i)=yy
 507         zzref(i)=zz
 508         else
 509         xxref(i)=0.0d0
 510         yyref(i)=0.0d0
 511         zzref(i)=0.0d0
 512         endif
 513       enddo
 514       if (lprn) then
 515 #ifdef MPI
 516         if (me.eq.king.or..not.out1file) then
 517 #endif
 518         write (iout,*) "xxref,yyref,zzref"
 519         do i=2,nres
 520           write (iout,'(a3,i4,3f10.5)')
 521      &     restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
 522         enddo
 523 #ifdef MPI
 524         endif
 525 #endif
 526       endif
 527       return
 528       end
 529 c---------------------------------------------------------------------------
 530       subroutine sccenter(ires,nscat,sccor)
 531       implicit none
 532       include 'DIMENSIONS'
 533       include 'COMMON.CHAIN'
 534       integer i,j,ires,nscat
 535       double precision sccor(3,50)
 536       double precision sccmj
 537       do j=1,3
 538         sccmj=0.0D0
 539         do i=1,nscat
 540           sccmj=sccmj+sccor(j,i)
 541         enddo
 542         dc(j,ires)=sccmj/nscat
 543       enddo
 544       return
 545       end
 546 c---------------------------------------------------------------------------
 547       subroutine bond_regular
 548       implicit none
 549       include 'DIMENSIONS'
 550       include 'COMMON.VAR'
 551       include 'COMMON.LOCAL'
 552       include 'COMMON.INTERACT'
 553       include 'COMMON.CHAIN'
 554       integer i,i1,i2
 555       do i=1,nres-1
 556        vbld(i+1)=vbl
 557        vbld_inv(i+1)=vblinv
 558        vbld(i+1+nres)=dsc(iabs(itype(i+1)))
 559        vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
 560 c       print *,vbld(i+1),vbld(i+1+nres)
 561       enddo
 562 c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
 563       do i=1,nchain
 564         i1=chain_border(1,i)
 565         i2=chain_border(2,i)
 566         if (i1.gt.1) then
 567           vbld(i1)=vbld(i1)/2
 568           vbld_inv(i1)=vbld_inv(i1)*2
 569         endif
 570         if (i2.lt.nres) then
 571           vbld(i2+1)=vbld(i2+1)/2
 572           vbld_inv(i2+1)=vbld_inv(i2+1)*2
 573         endif
 574       enddo
 575       return
 576       end
 577 c---------------------------------------------------------------------------
 578       subroutine readpdb_template(k)
 579 C Read the PDB file for read_constr_homology with read2sigma
 580 C and convert the peptide geometry into virtual-chain geometry.
 581       implicit none
 582       include 'DIMENSIONS'
 583       include 'COMMON.LOCAL'
 584       include 'COMMON.VAR'
 585       include 'COMMON.CHAIN'
 586       include 'COMMON.INTERACT'
 587       include 'COMMON.IOUNITS'
 588       include 'COMMON.GEO'
 589       include 'COMMON.NAMES'
 590       include 'COMMON.CONTROL'
 591       include 'COMMON.FRAG'
 592       include 'COMMON.SETUP'
 593       integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
 594      &  ishift_pdb,ires_ca
 595       logical lprn /.false./,fail
 596       double precision e1(3),e2(3),e3(3)
 597       double precision dcj,efree_temp
 598       character*3 seq,res
 599       character*5 atom
 600       character*80 card
 601       double precision sccor(3,20)
 602       integer rescode,iterter(maxres)
 603       do i=1,maxres
 604          iterter(i)=0
 605       enddo
 606       ibeg=1
 607       ishift1=0
 608       ishift=0
 609 c      write (2,*) "UNRES_PDB",unres_pdb
 610       ires=0
 611       ires_old=0
 612       iii=0
 613       lsecondary=.false.
 614       nhfrag=0
 615       nbfrag=0
 616       do
 617         read (ipdbin,'(a80)',end=10) card
 618         if (card(:3).eq.'END') then
 619           goto 10
 620         else if (card(:3).eq.'TER') then
 621 C End current chain
 622           ires_old=ires+2
 623           itype(ires_old-1)=ntyp1
 624           iterter(ires_old-1)=1
 625           itype(ires_old)=ntyp1
 626           iterter(ires_old)=1
 627           ibeg=2
 628 c          write (iout,*) "Chain ended",ires,ishift,ires_old
 629           if (unres_pdb) then
 630             do j=1,3
 631               dc(j,ires)=sccor(j,iii)
 632             enddo
 633           else
 634             call sccenter(ires,iii,sccor)
 635           endif
 636         endif
 637 C Fish out the ATOM cards.
 638         if (index(card(1:4),'ATOM').gt.0) then
 639           read (card(12:16),*) atom
 640 c          write (iout,*) "! ",atom," !",ires
 641 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 642           read (card(23:26),*) ires
 643           read (card(18:20),'(a3)') res
 644 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 645 c     &      " ires_old",ires_old
 646 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 647 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 648           if (ires-ishift+ishift1.ne.ires_old) then
 649 C Calculate the CM of the preceding residue.
 650             if (ibeg.eq.0) then
 651               if (unres_pdb) then
 652                 do j=1,3
 653                   dc(j,ires)=sccor(j,iii)
 654                 enddo
 655               else
 656                 call sccenter(ires_old,iii,sccor)
 657               endif
 658               iii=0
 659             endif
 660 C Start new residue.
 661             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 662               ires=ires_old
 663               cycle
 664             else if (ibeg.eq.1) then
 665 c              write (iout,*) "BEG ires",ires
 666               ishift=ires-1
 667               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 668                 ishift=ishift-1
 669                 itype(1)=ntyp1
 670               endif
 671               ires=ires-ishift+ishift1
 672               ires_old=ires
 673 c              write (iout,*) "ishift",ishift," ires",ires,
 674 c     &         " ires_old",ires_old
 675 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
 676               ibeg=0
 677             else if (ibeg.eq.2) then
 678 c Start a new chain
 679               ishift=-ires_old+ires-1
 680               ires=ires_old+1
 681 c              write (iout,*) "New chain started",ires,ishift
 682               ibeg=0
 683             else
 684               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 685               ires=ires-ishift+ishift1
 686               ires_old=ires
 687             endif
 688             if (res.eq.'ACE' .or. res.eq.'NHE') then
 689               itype(ires)=10
 690             else
 691               itype(ires)=rescode(ires,res,0)
 692             endif
 693           else
 694             ires=ires-ishift+ishift1
 695           endif
 696 c          write (iout,*) "ires_old",ires_old," ires",ires
 697 c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 698 c            ishift1=ishift1+1
 699 c          endif
 700 c          write (2,*) "ires",ires," res ",res," ity",ity
 701           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 702      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 703             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 704 c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
 705 #ifdef DEBUG
 706             write (iout,'(2i3,2x,a,3f8.3)')
 707      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 708 #endif
 709             iii=iii+1
 710             do j=1,3
 711               sccor(j,iii)=c(j,ires)
 712             enddo
 713             if (ishift.ne.0) then
 714               ires_ca=ires+ishift-ishift1
 715             else
 716               ires_ca=ires
 717             endif
 718 c            write (*,*) card(23:27),ires,itype(ires)
 719           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 720      &             atom.ne.'N' .and. atom.ne.'C' .and.
 721      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 722      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 723 c            write (iout,*) "sidechain ",atom
 724             iii=iii+1
 725             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 726           endif
 727         endif
 728       enddo
 729    10 if(me.eq.king.or..not.out1file)
 730      & write (iout,'(a,i5)') ' Nres: ',ires
 731 C Calculate dummy residue coordinates inside the "chain" of a multichain
 732 C system
 733       nres=ires
 734       do i=2,nres-1
 735 c        write (iout,*) i,itype(i),itype(i+1)
 736         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 737          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 738 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 739 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 740 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 741            if (unres_pdb) then
 742 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 743             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 744             if (fail) then
 745               e2(1)=0.0d0
 746               e2(2)=1.0d0
 747               e2(3)=0.0d0
 748             endif !fail
 749             do j=1,3
 750              c(j,i)=c(j,i-1)-1.9d0*e2(j)
 751             enddo
 752            else   !unres_pdb
 753            do j=1,3
 754              dcj=(c(j,i-2)-c(j,i-3))/2.0
 755             if (dcj.eq.0) dcj=1.23591524223
 756              c(j,i)=c(j,i-1)+dcj
 757              c(j,nres+i)=c(j,i)
 758            enddo
 759           endif   !unres_pdb
 760          else     !itype(i+1).eq.ntyp1
 761           if (unres_pdb) then
 762 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 763             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 764             if (fail) then
 765               e2(1)=0.0d0
 766               e2(2)=1.0d0
 767               e2(3)=0.0d0
 768             endif
 769             do j=1,3
 770               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 771             enddo
 772           else !unres_pdb
 773            do j=1,3
 774             dcj=(c(j,i+3)-c(j,i+2))/2.0
 775             if (dcj.eq.0) dcj=1.23591524223
 776             c(j,i)=c(j,i+1)-dcj
 777             c(j,nres+i)=c(j,i)
 778            enddo
 779           endif !unres_pdb
 780          endif !itype(i+1).eq.ntyp1
 781         endif  !itype.eq.ntyp1
 782       enddo
 783 C Calculate the CM of the last side chain.
 784       if (unres_pdb) then
 785         do j=1,3
 786           dc(j,ires)=sccor(j,iii)
 787         enddo
 788       else
 789         call sccenter(ires,iii,sccor)
 790       endif
 791       nsup=nres
 792       nstart_sup=1
 793       if (itype(nres).ne.10) then
 794         nres=nres+1
 795         itype(nres)=ntyp1
 796         if (unres_pdb) then
 797 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 798           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 799           if (fail) then
 800             e2(1)=0.0d0
 801             e2(2)=1.0d0
 802             e2(3)=0.0d0
 803           endif
 804           do j=1,3
 805             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
 806           enddo
 807         else
 808         do j=1,3
 809           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 810             if (dcj.eq.0) dcj=1.23591524223
 811           c(j,nres)=c(j,nres-1)+dcj
 812           c(j,2*nres)=c(j,nres)
 813         enddo
 814       endif
 815       endif
 816       do i=2,nres-1
 817         do j=1,3
 818           c(j,i+nres)=dc(j,i)
 819         enddo
 820       enddo
 821       do j=1,3
 822         c(j,nres+1)=c(j,1)
 823         c(j,2*nres)=c(j,nres)
 824       enddo
 825       if (itype(1).eq.ntyp1) then
 826         nsup=nsup-1
 827         nstart_sup=2
 828         if (unres_pdb) then
 829 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 830           call refsys(2,3,4,e1,e2,e3,fail)
 831           if (fail) then
 832             e2(1)=0.0d0
 833             e2(2)=1.0d0
 834             e2(3)=0.0d0
 835           endif
 836           do j=1,3
 837             c(j,1)=c(j,2)-1.9d0*e2(j)
 838           enddo
 839         else
 840         do j=1,3
 841           dcj=(c(j,4)-c(j,3))/2.0
 842           c(j,1)=c(j,2)-dcj
 843           c(j,nres+1)=c(j,1)
 844         enddo
 845         endif
 846       endif
 847 C Copy the coordinates to reference coordinates
 848 c      do i=1,2*nres
 849 c        do j=1,3
 850 c          cref(j,i)=c(j,i)
 851 c        enddo
 852 c      enddo
 853 C Calculate internal coordinates.
 854       if (out_template_coord) then
 855       write (iout,'(/a)')
 856      &  "Cartesian coordinates of the reference structure"
 857       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 858      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 859       do ires=1,nres
 860         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 861      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 862      &    (c(j,ires+nres),j=1,3)
 863       enddo
 864       endif
 865 C Calculate internal coordinates.
 866       call int_from_cart(.true.,.true.)
 867       call sc_loc_geom(.false.)
 868       do i=1,nres
 869         thetaref(i)=theta(i)
 870         phiref(i)=phi(i)
 871       enddo
 872       do i=1,nres-1
 873         do j=1,3
 874           dc(j,i)=c(j,i+1)-c(j,i)
 875           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 876         enddo
 877       enddo
 878       do i=2,nres-1
 879         do j=1,3
 880           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 881           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 882         enddo
 883 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 884 c     &   vbld_inv(i+nres)
 885       enddo
 886       do i=1,nres
 887         do j=1,3
 888           cref(j,i)=c(j,i)
 889           cref(j,i+nres)=c(j,i+nres)
 890         enddo
 891       enddo
 892       do i=1,2*nres
 893         do j=1,3
 894           chomo(j,i,k)=c(j,i)
 895         enddo
 896       enddo
 897
 898       return
 899       end
 900