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