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 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 c            if(me.eq.king.or..not.out1file)
 111 c     &       write (iout,'(2i3,2x,a,3f8.3)')
 112 c     &       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 c      write(iout,*)"before int_from_cart nres",nres
 255       call int_from_cart(.true.,.false.)
 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 c        write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
 266 c     &   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 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 274 c     &   vbld_inv(i+nres)
 275       enddo
 276       call sc_loc_geom(.false.)
 277       call int_from_cart1(.false.)
 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 c        write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
 354 c     &      vbld_inv(i+1)
 355         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 356         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 357       enddo
 358 c      if (unres_pdb) then
 359 c        if (itype(1).eq.21) then
 360 c          theta(3)=90.0d0*deg2rad
 361 c          phi(4)=180.0d0*deg2rad
 362 c          vbld(2)=3.8d0
 363 c          vbld_inv(2)=1.0d0/vbld(2)
 364 c        endif
 365 c        if (itype(nres).eq.21) then
 366 c          theta(nres)=90.0d0*deg2rad
 367 c          phi(nres)=180.0d0*deg2rad
 368 c          vbld(nres)=3.8d0
 369 c          vbld_inv(nres)=1.0d0/vbld(2)
 370 c        endif
 371 c      endif
 372 c      print *,"A TU2"
 373       if (lside) then
 374         do i=2,nres-1
 375           do j=1,3
 376             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 377      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 378           enddo
 379           iti=itype(i)
 380           di=dist(i,nres+i)
 381           vbld(i+nres)=di
 382           if (itype(i).ne.10) then
 383             vbld_inv(i+nres)=1.0d0/di
 384           else
 385             vbld_inv(i+nres)=0.0d0
 386           endif
 387           if (iti.ne.10) then
 388             alph(i)=alpha(nres+i,i,maxres2)
 389             omeg(i)=beta(nres+i,i,maxres2,i+1)
 390           endif
 391           if(me.eq.king.or..not.out1file)then
 392            if (lprn)
 393      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 394      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 395      &     rad2deg*alph(i),rad2deg*omeg(i)
 396           endif
 397         enddo
 398         if (lprn) then
 399            i=nres
 400            iti=itype(nres)
 401            write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 402      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 403      &     rad2deg*alph(i),rad2deg*omeg(i)
 404         endif
 405       else if (lprn) then
 406         do i=2,nres
 407           iti=itype(i)
 408           if(me.eq.king.or..not.out1file)
 409      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 410      &     rad2deg*theta(i),rad2deg*phi(i)
 411         enddo
 412       endif
 413       return
 414       end
 415 c-------------------------------------------------------------------------------
 416       subroutine sc_loc_geom(lprn)
 417       implicit real*8 (a-h,o-z)
 418       include 'DIMENSIONS'
 419 #ifdef MPI
 420       include "mpif.h"
 421 #endif
 422       include 'COMMON.LOCAL'
 423       include 'COMMON.VAR'
 424       include 'COMMON.CHAIN'
 425       include 'COMMON.INTERACT'
 426       include 'COMMON.IOUNITS'
 427       include 'COMMON.GEO'
 428       include 'COMMON.NAMES'
 429       include 'COMMON.CONTROL'
 430       include 'COMMON.SETUP'
 431       double precision x_prime(3),y_prime(3),z_prime(3)
 432       logical lprn
 433       do i=1,nres-1
 434         do j=1,3
 435           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 436         enddo
 437 c        write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
 438 c     &    " vbld",vbld_inv(i+1)
 439       enddo
 440       do i=2,nres-1
 441         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 442           do j=1,3
 443             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 444           enddo
 445 c          write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
 446 c     &    " vbld",vbld_inv(i+nres)
 447         else
 448           do j=1,3
 449             dc_norm(j,i+nres)=0.0d0
 450           enddo
 451         endif
 452       enddo
 453       do i=2,nres-1
 454         costtab(i+1) =dcos(theta(i+1))
 455         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 456         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 457         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 458         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 459         cosfac=dsqrt(cosfac2)
 460         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 461         sinfac=dsqrt(sinfac2)
 462         it=itype(i)
 463 c        write (iout,*) "i",i," costab",costtab(i+1),
 464 c     &    " sintab",sinttab(i+1)
 465 c        write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
 466 c        write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
 467         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 468 c
 469 C  Compute the axes of tghe local cartesian coordinates system; store in
 470 c   x_prime, y_prime and z_prime
 471 c
 472         do j=1,3
 473           x_prime(j) = 0.00
 474           y_prime(j) = 0.00
 475           z_prime(j) = 0.00
 476         enddo
 477         do j = 1,3
 478           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 479           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 480         enddo
 481 c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
 482 c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
 483         call vecpr(x_prime,y_prime,z_prime)
 484 c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
 485 c
 486 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 487 C to local coordinate system. Store in xx, yy, zz.
 488 c
 489         xx=0.0d0
 490         yy=0.0d0
 491         zz=0.0d0
 492         do j = 1,3
 493           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 494           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 495           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 496         enddo
 497
 498         xxref(i)=xx
 499         yyref(i)=yy
 500         zzref(i)=zz
 501         else
 502         xxref(i)=0.0d0
 503         yyref(i)=0.0d0
 504         zzref(i)=0.0d0
 505         endif
 506       enddo
 507       if (lprn) then
 508 #ifdef MPI
 509         if (me.eq.king.or..not.out1file) then
 510 #endif
 511         write (iout,*) "xxref,yyref,zzref"
 512         do i=2,nres
 513           write (iout,'(a3,i4,3f10.5)')
 514      &     restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
 515         enddo
 516 #ifdef MPI
 517         endif
 518 #endif
 519       endif
 520       return
 521       end
 522 c---------------------------------------------------------------------------
 523       subroutine sccenter(ires,nscat,sccor)
 524       implicit real*8 (a-h,o-z)
 525       include 'DIMENSIONS'
 526       include 'COMMON.CHAIN'
 527       dimension sccor(3,50)
 528       do j=1,3
 529         sccmj=0.0D0
 530         do i=1,nscat
 531           sccmj=sccmj+sccor(j,i)
 532         enddo
 533         dc(j,ires)=sccmj/nscat
 534       enddo
 535       return
 536       end
 537 c---------------------------------------------------------------------------
 538       subroutine bond_regular
 539       implicit real*8 (a-h,o-z)
 540       include 'DIMENSIONS'
 541       include 'COMMON.VAR'
 542       include 'COMMON.LOCAL'
 543       include 'COMMON.CALC'
 544       include 'COMMON.INTERACT'
 545       include 'COMMON.CHAIN'
 546       do i=1,nres-1
 547        vbld(i+1)=vbl
 548        vbld_inv(i+1)=vblinv
 549        vbld(i+1+nres)=dsc(iabs(itype(i+1)))
 550        vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
 551 c       print *,vbld(i+1),vbld(i+1+nres)
 552       enddo
 553 c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
 554       do i=1,nchain
 555         i1=chain_border(1,i)
 556         i2=chain_border(2,i)
 557         if (i1.gt.1) then
 558           vbld(i1)=vbld(i1)/2
 559           vbld_inv(i1)=vbld_inv(i1)*2
 560         endif
 561         if (i2.lt.nres) then
 562           vbld(i2+1)=vbld(i2+1)/2
 563           vbld_inv(i2+1)=vbld_inv(i2+1)*2
 564         endif
 565       enddo
 566       return
 567       end
 568 c---------------------------------------------------------------------------
 569       subroutine readpdb_template(k)
 570 C Read the PDB file for read_constr_homology with read2sigma
 571 C and convert the peptide geometry into virtual-chain geometry.
 572       implicit real*8 (a-h,o-z)
 573       include 'DIMENSIONS'
 574       include 'COMMON.LOCAL'
 575       include 'COMMON.VAR'
 576       include 'COMMON.CHAIN'
 577       include 'COMMON.INTERACT'
 578       include 'COMMON.IOUNITS'
 579       include 'COMMON.GEO'
 580       include 'COMMON.NAMES'
 581       include 'COMMON.CONTROL'
 582       include 'COMMON.DISTFIT'
 583       include 'COMMON.SETUP'
 584       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
 585      &  ishift_pdb
 586       logical lprn /.false./,fail
 587       double precision e1(3),e2(3),e3(3)
 588       double precision dcj,efree_temp
 589       character*3 seq,res
 590       character*5 atom
 591       character*80 card
 592       double precision sccor(3,20)
 593       integer rescode,iterter(maxres)
 594       do i=1,maxres
 595          iterter(i)=0
 596       enddo
 597       ibeg=1
 598       ishift1=0
 599       ishift=0
 600 c      write (2,*) "UNRES_PDB",unres_pdb
 601       ires=0
 602       ires_old=0
 603       iii=0
 604       lsecondary=.false.
 605       nhfrag=0
 606       nbfrag=0
 607       do
 608         read (ipdbin,'(a80)',end=10) card
 609         if (card(:3).eq.'END') then
 610           goto 10
 611         else if (card(:3).eq.'TER') then
 612 C End current chain
 613           ires_old=ires+2
 614           itype(ires_old-1)=ntyp1
 615           iterter(ires_old-1)=1
 616           itype(ires_old)=ntyp1
 617           iterter(ires_old)=1
 618           ibeg=2
 619 c          write (iout,*) "Chain ended",ires,ishift,ires_old
 620           if (unres_pdb) then
 621             do j=1,3
 622               dc(j,ires)=sccor(j,iii)
 623             enddo
 624           else
 625             call sccenter(ires,iii,sccor)
 626           endif
 627         endif
 628 C Fish out the ATOM cards.
 629         if (index(card(1:4),'ATOM').gt.0) then
 630           read (card(12:16),*) atom
 631 c          write (iout,*) "! ",atom," !",ires
 632 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 633           read (card(23:26),*) ires
 634           read (card(18:20),'(a3)') res
 635 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 636 c     &      " ires_old",ires_old
 637 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 638 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 639           if (ires-ishift+ishift1.ne.ires_old) then
 640 C Calculate the CM of the preceding residue.
 641             if (ibeg.eq.0) then
 642               if (unres_pdb) then
 643                 do j=1,3
 644                   dc(j,ires)=sccor(j,iii)
 645                 enddo
 646               else
 647                 call sccenter(ires_old,iii,sccor)
 648               endif
 649               iii=0
 650             endif
 651 C Start new residue.
 652             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 653               ires=ires_old
 654               cycle
 655             else if (ibeg.eq.1) then
 656 c              write (iout,*) "BEG ires",ires
 657               ishift=ires-1
 658               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 659                 ishift=ishift-1
 660                 itype(1)=ntyp1
 661               endif
 662               ires=ires-ishift+ishift1
 663               ires_old=ires
 664 c              write (iout,*) "ishift",ishift," ires",ires,
 665 c     &         " ires_old",ires_old
 666 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
 667               ibeg=0
 668             else if (ibeg.eq.2) then
 669 c Start a new chain
 670               ishift=-ires_old+ires-1
 671               ires=ires_old+1
 672 c              write (iout,*) "New chain started",ires,ishift
 673               ibeg=0
 674             else
 675               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 676               ires=ires-ishift+ishift1
 677               ires_old=ires
 678             endif
 679             if (res.eq.'ACE' .or. res.eq.'NHE') then
 680               itype(ires)=10
 681             else
 682               itype(ires)=rescode(ires,res,0)
 683             endif
 684           else
 685             ires=ires-ishift+ishift1
 686           endif
 687 c          write (iout,*) "ires_old",ires_old," ires",ires
 688 c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 689 c            ishift1=ishift1+1
 690 c          endif
 691 c          write (2,*) "ires",ires," res ",res," ity",ity
 692           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 693      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 694             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 695 c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
 696 #ifdef DEBUG
 697             write (iout,'(2i3,2x,a,3f8.3)')
 698      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 699 #endif
 700             iii=iii+1
 701             do j=1,3
 702               sccor(j,iii)=c(j,ires)
 703             enddo
 704             if (ishift.ne.0) then
 705               ires_ca=ires+ishift-ishift1
 706             else
 707               ires_ca=ires
 708             endif
 709 c            write (*,*) card(23:27),ires,itype(ires)
 710           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 711      &             atom.ne.'N' .and. atom.ne.'C' .and.
 712      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 713      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 714 c            write (iout,*) "sidechain ",atom
 715             iii=iii+1
 716             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 717           endif
 718         endif
 719       enddo
 720    10 if(me.eq.king.or..not.out1file)
 721      & write (iout,'(a,i5)') ' Nres: ',ires
 722 C Calculate dummy residue coordinates inside the "chain" of a multichain
 723 C system
 724       nres=ires
 725       do i=2,nres-1
 726 c        write (iout,*) i,itype(i),itype(i+1)
 727         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 728          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 729 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 730 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 731 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 732            if (unres_pdb) then
 733 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 734             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 735             if (fail) then
 736               e2(1)=0.0d0
 737               e2(2)=1.0d0
 738               e2(3)=0.0d0
 739             endif !fail
 740             do j=1,3
 741              c(j,i)=c(j,i-1)-1.9d0*e2(j)
 742             enddo
 743            else   !unres_pdb
 744            do j=1,3
 745              dcj=(c(j,i-2)-c(j,i-3))/2.0
 746             if (dcj.eq.0) dcj=1.23591524223
 747              c(j,i)=c(j,i-1)+dcj
 748              c(j,nres+i)=c(j,i)
 749            enddo
 750           endif   !unres_pdb
 751          else     !itype(i+1).eq.ntyp1
 752           if (unres_pdb) then
 753 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 754             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 755             if (fail) then
 756               e2(1)=0.0d0
 757               e2(2)=1.0d0
 758               e2(3)=0.0d0
 759             endif
 760             do j=1,3
 761               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 762             enddo
 763           else !unres_pdb
 764            do j=1,3
 765             dcj=(c(j,i+3)-c(j,i+2))/2.0
 766             if (dcj.eq.0) dcj=1.23591524223
 767             c(j,i)=c(j,i+1)-dcj
 768             c(j,nres+i)=c(j,i)
 769            enddo
 770           endif !unres_pdb
 771          endif !itype(i+1).eq.ntyp1
 772         endif  !itype.eq.ntyp1
 773       enddo
 774 C Calculate the CM of the last side chain.
 775       if (unres_pdb) then
 776         do j=1,3
 777           dc(j,ires)=sccor(j,iii)
 778         enddo
 779       else
 780         call sccenter(ires,iii,sccor)
 781       endif
 782       nsup=nres
 783       nstart_sup=1
 784       if (itype(nres).ne.10) then
 785         nres=nres+1
 786         itype(nres)=ntyp1
 787         if (unres_pdb) then
 788 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 789           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 790           if (fail) then
 791             e2(1)=0.0d0
 792             e2(2)=1.0d0
 793             e2(3)=0.0d0
 794           endif
 795           do j=1,3
 796             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
 797           enddo
 798         else
 799         do j=1,3
 800           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 801             if (dcj.eq.0) dcj=1.23591524223
 802           c(j,nres)=c(j,nres-1)+dcj
 803           c(j,2*nres)=c(j,nres)
 804         enddo
 805       endif
 806       endif
 807       do i=2,nres-1
 808         do j=1,3
 809           c(j,i+nres)=dc(j,i)
 810         enddo
 811       enddo
 812       do j=1,3
 813         c(j,nres+1)=c(j,1)
 814         c(j,2*nres)=c(j,nres)
 815       enddo
 816       if (itype(1).eq.ntyp1) then
 817         nsup=nsup-1
 818         nstart_sup=2
 819         if (unres_pdb) then
 820 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 821           call refsys(2,3,4,e1,e2,e3,fail)
 822           if (fail) then
 823             e2(1)=0.0d0
 824             e2(2)=1.0d0
 825             e2(3)=0.0d0
 826           endif
 827           do j=1,3
 828             c(j,1)=c(j,2)-1.9d0*e2(j)
 829           enddo
 830         else
 831         do j=1,3
 832           dcj=(c(j,4)-c(j,3))/2.0
 833           c(j,1)=c(j,2)-dcj
 834           c(j,nres+1)=c(j,1)
 835         enddo
 836         endif
 837       endif
 838 C Copy the coordinates to reference coordinates
 839 c      do i=1,2*nres
 840 c        do j=1,3
 841 c          cref(j,i)=c(j,i)
 842 c        enddo
 843 c      enddo
 844 C Calculate internal coordinates.
 845       if (out_template_coord) then
 846       write (iout,'(/a)')
 847      &  "Cartesian coordinates of the reference structure"
 848       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 849      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 850       do ires=1,nres
 851         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 852      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 853      &    (c(j,ires+nres),j=1,3)
 854       enddo
 855       endif
 856 C Calculate internal coordinates.
 857       call int_from_cart(.true.,.true.)
 858       call sc_loc_geom(.false.)
 859       do i=1,nres
 860         thetaref(i)=theta(i)
 861         phiref(i)=phi(i)
 862       enddo
 863       do i=1,nres-1
 864         do j=1,3
 865           dc(j,i)=c(j,i+1)-c(j,i)
 866           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 867         enddo
 868       enddo
 869       do i=2,nres-1
 870         do j=1,3
 871           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 872           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 873         enddo
 874 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 875 c     &   vbld_inv(i+nres)
 876       enddo
 877       do i=1,nres
 878         do j=1,3
 879           cref(j,i)=c(j,i)
 880           cref(j,i+nres)=c(j,i+nres)
 881         enddo
 882       enddo
 883       do i=1,2*nres
 884         do j=1,3
 885           chomo(j,i,k)=c(j,i)
 886         enddo
 887       enddo
 888
 889       return
 890       end
 891