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