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