source/unres/src-5hdiag-tmp/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 c        write (iout,*) "i",i," costab",costtab(i+1),
 458 c     &    " sintab",sinttab(i+1)
 459 c        write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
 460 c        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 c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
 476 c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
 477         call vecpr(x_prime,y_prime,z_prime)
 478 c        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 c---------------------------------------------------------------------------
 550       subroutine readpdb_template(k)
 551 C Read the PDB file for read_constr_homology with read2sigma
 552 C and convert the peptide geometry into virtual-chain geometry.
 553       implicit real*8 (a-h,o-z)
 554       include 'DIMENSIONS'
 555       include 'COMMON.LOCAL'
 556       include 'COMMON.VAR'
 557       include 'COMMON.CHAIN'
 558       include 'COMMON.INTERACT'
 559       include 'COMMON.IOUNITS'
 560       include 'COMMON.GEO'
 561       include 'COMMON.NAMES'
 562       include 'COMMON.CONTROL'
 563       include 'COMMON.DISTFIT'
 564       include 'COMMON.SETUP'
 565       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
 566      &  ishift_pdb
 567       logical lprn /.false./,fail
 568       double precision e1(3),e2(3),e3(3)
 569       double precision dcj,efree_temp
 570       character*3 seq,res
 571       character*5 atom
 572       character*80 card
 573       double precision sccor(3,20)
 574       integer rescode,iterter(maxres)
 575       do i=1,maxres
 576          iterter(i)=0
 577       enddo
 578       ibeg=1
 579       ishift1=0
 580       ishift=0
 581 c      write (2,*) "UNRES_PDB",unres_pdb
 582       ires=0
 583       ires_old=0
 584       iii=0
 585       lsecondary=.false.
 586       nhfrag=0
 587       nbfrag=0
 588       do
 589         read (ipdbin,'(a80)',end=10) card
 590         if (card(:3).eq.'END') then
 591           goto 10
 592         else if (card(:3).eq.'TER') then
 593 C End current chain
 594           ires_old=ires+2
 595           itype(ires_old-1)=ntyp1
 596           iterter(ires_old-1)=1
 597           itype(ires_old)=ntyp1
 598           iterter(ires_old)=1
 599           ibeg=2
 600 c          write (iout,*) "Chain ended",ires,ishift,ires_old
 601           if (unres_pdb) then
 602             do j=1,3
 603               dc(j,ires)=sccor(j,iii)
 604             enddo
 605           else
 606             call sccenter(ires,iii,sccor)
 607           endif
 608         endif
 609 C Fish out the ATOM cards.
 610         if (index(card(1:4),'ATOM').gt.0) then
 611           read (card(12:16),*) atom
 612 c          write (iout,*) "! ",atom," !",ires
 613 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 614           read (card(23:26),*) ires
 615           read (card(18:20),'(a3)') res
 616 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 617 c     &      " ires_old",ires_old
 618 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 619 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 620           if (ires-ishift+ishift1.ne.ires_old) then
 621 C Calculate the CM of the preceding residue.
 622             if (ibeg.eq.0) then
 623               if (unres_pdb) then
 624                 do j=1,3
 625                   dc(j,ires)=sccor(j,iii)
 626                 enddo
 627               else
 628                 call sccenter(ires_old,iii,sccor)
 629               endif
 630               iii=0
 631             endif
 632 C Start new residue.
 633             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 634               ires=ires_old
 635               cycle
 636             else if (ibeg.eq.1) then
 637 c              write (iout,*) "BEG ires",ires
 638               ishift=ires-1
 639               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 640                 ishift=ishift-1
 641                 itype(1)=ntyp1
 642               endif
 643               ires=ires-ishift+ishift1
 644               ires_old=ires
 645 c              write (iout,*) "ishift",ishift," ires",ires,
 646 c     &         " ires_old",ires_old
 647 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
 648               ibeg=0
 649             else if (ibeg.eq.2) then
 650 c Start a new chain
 651               ishift=-ires_old+ires-1
 652               ires=ires_old+1
 653 c              write (iout,*) "New chain started",ires,ishift
 654               ibeg=0
 655             else
 656               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 657               ires=ires-ishift+ishift1
 658               ires_old=ires
 659             endif
 660             if (res.eq.'ACE' .or. res.eq.'NHE') then
 661               itype(ires)=10
 662             else
 663               itype(ires)=rescode(ires,res,0)
 664             endif
 665           else
 666             ires=ires-ishift+ishift1
 667           endif
 668 c          write (iout,*) "ires_old",ires_old," ires",ires
 669 c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 670 c            ishift1=ishift1+1
 671 c          endif
 672 c          write (2,*) "ires",ires," res ",res," ity",ity
 673           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 674      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 675             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 676 c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
 677 #ifdef DEBUG
 678             write (iout,'(2i3,2x,a,3f8.3)')
 679      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 680 #endif
 681             iii=iii+1
 682             do j=1,3
 683               sccor(j,iii)=c(j,ires)
 684             enddo
 685             if (ishift.ne.0) then
 686               ires_ca=ires+ishift-ishift1
 687             else
 688               ires_ca=ires
 689             endif
 690 c            write (*,*) card(23:27),ires,itype(ires)
 691           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 692      &             atom.ne.'N' .and. atom.ne.'C' .and.
 693      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 694      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 695 c            write (iout,*) "sidechain ",atom
 696             iii=iii+1
 697             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 698           endif
 699         endif
 700       enddo
 701    10 if(me.eq.king.or..not.out1file)
 702      & write (iout,'(a,i5)') ' Nres: ',ires
 703 C Calculate dummy residue coordinates inside the "chain" of a multichain
 704 C system
 705       nres=ires
 706       do i=2,nres-1
 707 c        write (iout,*) i,itype(i),itype(i+1)
 708         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 709          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 710 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 711 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 712 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 713            if (unres_pdb) then
 714 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 715             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 716             if (fail) then
 717               e2(1)=0.0d0
 718               e2(2)=1.0d0
 719               e2(3)=0.0d0
 720             endif !fail
 721             do j=1,3
 722              c(j,i)=c(j,i-1)-1.9d0*e2(j)
 723             enddo
 724            else   !unres_pdb
 725            do j=1,3
 726              dcj=(c(j,i-2)-c(j,i-3))/2.0
 727             if (dcj.eq.0) dcj=1.23591524223
 728              c(j,i)=c(j,i-1)+dcj
 729              c(j,nres+i)=c(j,i)
 730            enddo
 731           endif   !unres_pdb
 732          else     !itype(i+1).eq.ntyp1
 733           if (unres_pdb) then
 734 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 735             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 736             if (fail) then
 737               e2(1)=0.0d0
 738               e2(2)=1.0d0
 739               e2(3)=0.0d0
 740             endif
 741             do j=1,3
 742               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 743             enddo
 744           else !unres_pdb
 745            do j=1,3
 746             dcj=(c(j,i+3)-c(j,i+2))/2.0
 747             if (dcj.eq.0) dcj=1.23591524223
 748             c(j,i)=c(j,i+1)-dcj
 749             c(j,nres+i)=c(j,i)
 750            enddo
 751           endif !unres_pdb
 752          endif !itype(i+1).eq.ntyp1
 753         endif  !itype.eq.ntyp1
 754       enddo
 755 C Calculate the CM of the last side chain.
 756       if (unres_pdb) then
 757         do j=1,3
 758           dc(j,ires)=sccor(j,iii)
 759         enddo
 760       else
 761         call sccenter(ires,iii,sccor)
 762       endif
 763       nsup=nres
 764       nstart_sup=1
 765       if (itype(nres).ne.10) then
 766         nres=nres+1
 767         itype(nres)=ntyp1
 768         if (unres_pdb) then
 769 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 770           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 771           if (fail) then
 772             e2(1)=0.0d0
 773             e2(2)=1.0d0
 774             e2(3)=0.0d0
 775           endif
 776           do j=1,3
 777             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
 778           enddo
 779         else
 780         do j=1,3
 781           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 782             if (dcj.eq.0) dcj=1.23591524223
 783           c(j,nres)=c(j,nres-1)+dcj
 784           c(j,2*nres)=c(j,nres)
 785         enddo
 786       endif
 787       endif
 788       do i=2,nres-1
 789         do j=1,3
 790           c(j,i+nres)=dc(j,i)
 791         enddo
 792       enddo
 793       do j=1,3
 794         c(j,nres+1)=c(j,1)
 795         c(j,2*nres)=c(j,nres)
 796       enddo
 797       if (itype(1).eq.ntyp1) then
 798         nsup=nsup-1
 799         nstart_sup=2
 800         if (unres_pdb) then
 801 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 802           call refsys(2,3,4,e1,e2,e3,fail)
 803           if (fail) then
 804             e2(1)=0.0d0
 805             e2(2)=1.0d0
 806             e2(3)=0.0d0
 807           endif
 808           do j=1,3
 809             c(j,1)=c(j,2)-1.9d0*e2(j)
 810           enddo
 811         else
 812         do j=1,3
 813           dcj=(c(j,4)-c(j,3))/2.0
 814           c(j,1)=c(j,2)-dcj
 815           c(j,nres+1)=c(j,1)
 816         enddo
 817         endif
 818       endif
 819 C Copy the coordinates to reference coordinates
 820 c      do i=1,2*nres
 821 c        do j=1,3
 822 c          cref(j,i)=c(j,i)
 823 c        enddo
 824 c      enddo
 825 C Calculate internal coordinates.
 826       if (out_template_coord) then
 827       write (iout,'(/a)')
 828      &  "Cartesian coordinates of the reference structure"
 829       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 830      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 831       do ires=1,nres
 832         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 833      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 834      &    (c(j,ires+nres),j=1,3)
 835       enddo
 836       endif
 837 C Calculate internal coordinates.
 838       call int_from_cart(.true.,.false.)
 839       call sc_loc_geom(.false.)
 840       do i=1,nres
 841         thetaref(i)=theta(i)
 842         phiref(i)=phi(i)
 843       enddo
 844       do i=1,nres-1
 845         do j=1,3
 846           dc(j,i)=c(j,i+1)-c(j,i)
 847           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 848         enddo
 849       enddo
 850       do i=2,nres-1
 851         do j=1,3
 852           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 853           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 854         enddo
 855 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 856 c     &   vbld_inv(i+nres)
 857       enddo
 858       do i=1,nres
 859         do j=1,3
 860           cref(j,i)=c(j,i)
 861           cref(j,i+nres)=c(j,i+nres)
 862         enddo
 863       enddo
 864       do i=1,2*nres
 865         do j=1,3
 866           chomo(j,i,k)=c(j,i)
 867         enddo
 868       enddo
 869
 870       return
 871       end
 872