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