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