source/cluster/wham/src-HCD-5D/readpdb.F

   1       subroutine readpdb(lprint)
   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 'COMMON.CONTROL'
   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.SBRIDGE'
  15       character*3 seq,atom,res
  16       character*80 card
  17       double precision sccor(3,50)
  18       integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
  19       double precision dcj
  20       integer rescode,kkk,lll,icha,cou,kupa,iprzes
  21       logical lprint
  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       if (lprint) then
 176       write (iout,100)
 177       do ires=1,nres
 178         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 179      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 180      &    (c(j,nres+ires),j=1,3)
 181       enddo
 182       endif
 183       call int_from_cart(.true.,.false.)
 184       call flush(iout)
 185       do i=1,nres-1
 186         do j=1,3
 187           dc(j,i)=c(j,i+1)-c(j,i)
 188           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 189         enddo
 190       enddo
 191       do i=2,nres-1
 192         do j=1,3
 193           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 194           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 195         enddo
 196 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 197 c     &   vbld_inv(i+nres)
 198       enddo
 199 c      call chainbuild
 200 C Copy the coordinates to reference coordinates
 201       do i=1,nres
 202         do j=1,3
 203           cref(j,i)=c(j,i)
 204           cref(j,i+nres)=c(j,i+nres)
 205         enddo
 206       enddo
 207   100 format ('Residue    alpha-carbon coordinates    ',
 208      &          '     centroid coordinates'/
 209      1          '         ', 6X,'X',7X,'Y',7X,'Z',
 210      &                          12X,'X',7X,'Y',7X,'Z')
 211   110 format (a,'(',i3,')',6f12.5)
 212
 213       ishift_pdb=ishift
 214       return
 215       end
 216 c---------------------------------------------------------------------------
 217       subroutine int_from_cart(lside,lprn)
 218       implicit none
 219       include 'DIMENSIONS'
 220       include 'COMMON.LOCAL'
 221       include 'COMMON.VAR'
 222       include 'COMMON.CHAIN'
 223       include 'COMMON.INTERACT'
 224       include 'COMMON.IOUNITS'
 225       include 'COMMON.GEO'
 226       include 'COMMON.NAMES'
 227       character*3 seq,atom,res
 228       character*80 card
 229       double precision sccor(3,50)
 230       integer rescode
 231       double precision dist,alpha,beta,di
 232       integer i,j,iti
 233       logical lside,lprn
 234       if (lprn) then
 235         write (iout,'(/a)')
 236      &  'Internal coordinates calculated from crystal structure.'
 237         if (lside) then
 238           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 239      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 240      & '     Omega'
 241         else
 242           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 243      & '       Phi'
 244         endif
 245       endif
 246       do i=2,nres
 247         iti=itype(i)
 248 c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
 249         if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
 250      &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
 251           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 252           stop
 253         endif
 254         vbld(i)=dist(i-1,i)
 255         vbld_inv(i)=1.0d0/vbld(i)
 256         theta(i+1)=alpha(i-1,i,i+1)
 257         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 258       enddo
 259 c      if (itype(1).eq.ntyp1) then
 260 c        do j=1,3
 261 c          c(j,1)=c(j,2)+(c(j,3)-c(j,4))
 262 c        enddo
 263 c      endif
 264 c      if (itype(nres).eq.ntyp1) then
 265 c        do j=1,3
 266 c          c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
 267 c        enddo
 268 c      endif
 269       if (lside) then
 270         do i=2,nres-1
 271           do j=1,3
 272             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 273           enddo
 274           iti=itype(i)
 275           di=dist(i,nres+i)
 276            vbld(i+nres)=di
 277           if (itype(i).ne.10) then
 278             vbld_inv(i+nres)=1.0d0/di
 279           else
 280             vbld_inv(i+nres)=0.0d0
 281           endif
 282           if (iti.ne.10) then
 283             alph(i)=alpha(nres+i,i,maxres2)
 284             omeg(i)=beta(nres+i,i,maxres2,i+1)
 285           endif
 286           if (iti.ne.10) then
 287             alph(i)=alpha(nres+i,i,maxres2)
 288             omeg(i)=beta(nres+i,i,maxres2,i+1)
 289           endif
 290           if (lprn)
 291      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 292      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 293      &    rad2deg*alph(i),rad2deg*omeg(i)
 294         enddo
 295       else if (lprn) then
 296         do i=2,nres
 297           iti=itype(i)
 298           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 299      &    rad2deg*theta(i),rad2deg*phi(i)
 300         enddo
 301       endif
 302       return
 303       end
 304 c---------------------------------------------------------------------------
 305       subroutine sccenter(ires,nscat,sccor)
 306       implicit none
 307       include 'DIMENSIONS'
 308       include 'COMMON.CHAIN'
 309       integer ires,nscat,i,j
 310       double precision sccor(3,50),sccmj
 311       do j=1,3
 312         sccmj=0.0D0
 313         do i=1,nscat
 314           sccmj=sccmj+sccor(j,i)
 315         enddo
 316         dc(j,ires)=sccmj/nscat
 317       enddo
 318       return
 319       end
 320 c---------------------------------------------------------------------------
 321       subroutine sc_loc_geom(lprn)
 322       implicit real*8 (a-h,o-z)
 323       include 'DIMENSIONS'
 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 'COMMON.LOCAL'
 436       include 'COMMON.VAR'
 437       include 'COMMON.CHAIN'
 438       include 'COMMON.INTERACT'
 439       include 'COMMON.IOUNITS'
 440       include 'COMMON.GEO'
 441       include 'COMMON.NAMES'
 442       include 'COMMON.CONTROL'
 443       include 'COMMON.SETUP'
 444       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
 445       logical lprn /.false./,fail
 446       double precision e1(3),e2(3),e3(3)
 447       double precision dcj,efree_temp
 448       character*3 seq,res
 449       character*5 atom
 450       character*80 card
 451       double precision sccor(3,20)
 452       integer rescode,iterter(maxres)
 453       do i=1,maxres
 454          iterter(i)=0
 455       enddo
 456       ibeg=1
 457       ishift1=0
 458       ishift=0
 459 c      write (2,*) "UNRES_PDB",unres_pdb
 460       ires=0
 461       ires_old=0
 462       iii=0
 463       lsecondary=.false.
 464       nhfrag=0
 465       nbfrag=0
 466       do
 467         read (ipdbin,'(a80)',end=10) card
 468         if (card(:3).eq.'END') then
 469           goto 10
 470         else if (card(:3).eq.'TER') then
 471 C End current chain
 472           ires_old=ires+2
 473           itype(ires_old-1)=ntyp1
 474           iterter(ires_old-1)=1
 475           itype(ires_old)=ntyp1
 476           iterter(ires_old)=1
 477           ibeg=2
 478 c          write (iout,*) "Chain ended",ires,ishift,ires_old
 479           if (unres_pdb) then
 480             do j=1,3
 481               dc(j,ires)=sccor(j,iii)
 482             enddo
 483           else
 484             call sccenter(ires,iii,sccor)
 485           endif
 486         endif
 487 C Fish out the ATOM cards.
 488         if (index(card(1:4),'ATOM').gt.0) then
 489           read (card(12:16),*) atom
 490 c          write (iout,*) "! ",atom," !",ires
 491 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
 492           read (card(23:26),*) ires
 493           read (card(18:20),'(a3)') res
 494 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 495 c     &      " ires_old",ires_old
 496 c          write (iout,*) "ishift",ishift," ishift1",ishift1
 497 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 498           if (ires-ishift+ishift1.ne.ires_old) then
 499 C Calculate the CM of the preceding residue.
 500             if (ibeg.eq.0) then
 501               if (unres_pdb) then
 502                 do j=1,3
 503                   dc(j,ires)=sccor(j,iii)
 504                 enddo
 505               else
 506                 call sccenter(ires_old,iii,sccor)
 507               endif
 508               iii=0
 509             endif
 510 C Start new residue.
 511             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 512               ires=ires_old
 513               cycle
 514             else if (ibeg.eq.1) then
 515 c              write (iout,*) "BEG ires",ires
 516               ishift=ires-1
 517               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 518                 ishift=ishift-1
 519                 itype(1)=ntyp1
 520               endif
 521               ires=ires-ishift+ishift1
 522               ires_old=ires
 523 c              write (iout,*) "ishift",ishift," ires",ires,
 524 c     &         " ires_old",ires_old
 525 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
 526               ibeg=0
 527             else if (ibeg.eq.2) then
 528 c Start a new chain
 529               ishift=-ires_old+ires-1
 530               ires=ires_old+1
 531 c              write (iout,*) "New chain started",ires,ishift
 532               ibeg=0
 533             else
 534               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 535               ires=ires-ishift+ishift1
 536               ires_old=ires
 537             endif
 538             if (res.eq.'ACE' .or. res.eq.'NHE') then
 539               itype(ires)=10
 540             else
 541               itype(ires)=rescode(ires,res,0)
 542             endif
 543           else
 544             ires=ires-ishift+ishift1
 545           endif
 546 c          write (iout,*) "ires_old",ires_old," ires",ires
 547 c          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 548 c            ishift1=ishift1+1
 549 c          endif
 550 c          write (2,*) "ires",ires," res ",res," ity",ity
 551           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 552      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 553             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 554 c            write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
 555 #ifdef DEBUG
 556             write (iout,'(2i3,2x,a,3f8.3)')
 557      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 558 #endif
 559             iii=iii+1
 560             do j=1,3
 561               sccor(j,iii)=c(j,ires)
 562             enddo
 563             if (ishift.ne.0) then
 564               ires_ca=ires+ishift-ishift1
 565             else
 566               ires_ca=ires
 567             endif
 568 c            write (*,*) card(23:27),ires,itype(ires)
 569           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 570      &             atom.ne.'N' .and. atom.ne.'C' .and.
 571      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 572      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 573 c            write (iout,*) "sidechain ",atom
 574             iii=iii+1
 575             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 576           endif
 577         endif
 578       enddo
 579    10 write (iout,'(a,i5)') ' Nres: ',ires
 580 C Calculate dummy residue coordinates inside the "chain" of a multichain
 581 C system
 582       nres=ires
 583       do i=2,nres-1
 584 c        write (iout,*) i,itype(i),itype(i+1)
 585         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 586          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 587 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 588 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 589 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 590            if (unres_pdb) then
 591 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 592             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 593             if (fail) then
 594               e2(1)=0.0d0
 595               e2(2)=1.0d0
 596               e2(3)=0.0d0
 597             endif !fail
 598             do j=1,3
 599              c(j,i)=c(j,i-1)-1.9d0*e2(j)
 600             enddo
 601            else   !unres_pdb
 602            do j=1,3
 603              dcj=(c(j,i-2)-c(j,i-3))/2.0
 604             if (dcj.eq.0) dcj=1.23591524223
 605              c(j,i)=c(j,i-1)+dcj
 606              c(j,nres+i)=c(j,i)
 607            enddo
 608           endif   !unres_pdb
 609          else     !itype(i+1).eq.ntyp1
 610           if (unres_pdb) then
 611 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 612             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 613             if (fail) then
 614               e2(1)=0.0d0
 615               e2(2)=1.0d0
 616               e2(3)=0.0d0
 617             endif
 618             do j=1,3
 619               c(j,i)=c(j,i+1)-1.9d0*e2(j)
 620             enddo
 621           else !unres_pdb
 622            do j=1,3
 623             dcj=(c(j,i+3)-c(j,i+2))/2.0
 624             if (dcj.eq.0) dcj=1.23591524223
 625             c(j,i)=c(j,i+1)-dcj
 626             c(j,nres+i)=c(j,i)
 627            enddo
 628           endif !unres_pdb
 629          endif !itype(i+1).eq.ntyp1
 630         endif  !itype.eq.ntyp1
 631       enddo
 632 C Calculate the CM of the last side chain.
 633       if (unres_pdb) then
 634         do j=1,3
 635           dc(j,ires)=sccor(j,iii)
 636         enddo
 637       else
 638         call sccenter(ires,iii,sccor)
 639       endif
 640       nsup=nres
 641       nstart_sup=1
 642       if (itype(nres).ne.10) then
 643         nres=nres+1
 644         itype(nres)=ntyp1
 645         if (unres_pdb) then
 646 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 647           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 648           if (fail) then
 649             e2(1)=0.0d0
 650             e2(2)=1.0d0
 651             e2(3)=0.0d0
 652           endif
 653           do j=1,3
 654             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
 655           enddo
 656         else
 657         do j=1,3
 658           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 659             if (dcj.eq.0) dcj=1.23591524223
 660           c(j,nres)=c(j,nres-1)+dcj
 661           c(j,2*nres)=c(j,nres)
 662         enddo
 663       endif
 664       endif
 665       do i=2,nres-1
 666         do j=1,3
 667           c(j,i+nres)=dc(j,i)
 668         enddo
 669       enddo
 670       do j=1,3
 671         c(j,nres+1)=c(j,1)
 672         c(j,2*nres)=c(j,nres)
 673       enddo
 674       if (itype(1).eq.ntyp1) then
 675         nsup=nsup-1
 676         nstart_sup=2
 677         if (unres_pdb) then
 678 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 679           call refsys(2,3,4,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
 685           do j=1,3
 686             c(j,1)=c(j,2)-1.9d0*e2(j)
 687           enddo
 688         else
 689         do j=1,3
 690           dcj=(c(j,4)-c(j,3))/2.0
 691           c(j,1)=c(j,2)-dcj
 692           c(j,nres+1)=c(j,1)
 693         enddo
 694         endif
 695       endif
 696 C Copy the coordinates to reference coordinates
 697 c      do i=1,2*nres
 698 c        do j=1,3
 699 c          cref(j,i)=c(j,i)
 700 c        enddo
 701 c      enddo
 702 C Calculate internal coordinates.
 703       if (out_template_coord) then
 704       write (iout,'(/a)')
 705      &  "Cartesian coordinates of the reference structure"
 706       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 707      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 708       do ires=1,nres
 709         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 710      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 711      &    (c(j,ires+nres),j=1,3)
 712       enddo
 713       endif
 714 C Calculate internal coordinates.
 715 c      call int_from_cart1(.false.)
 716       call int_from_cart(.true.,.true.)
 717       call sc_loc_geom(.true.)
 718       do i=1,nres
 719         thetaref(i)=theta(i)
 720         phiref(i)=phi(i)
 721       enddo
 722       do i=1,nres-1
 723         do j=1,3
 724           dc(j,i)=c(j,i+1)-c(j,i)
 725           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 726         enddo
 727       enddo
 728       do i=2,nres-1
 729         do j=1,3
 730           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 731           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 732         enddo
 733 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 734 c     &   vbld_inv(i+nres)
 735       enddo
 736       do i=1,nres
 737         do j=1,3
 738           cref(j,i)=c(j,i)
 739           cref(j,i+nres)=c(j,i+nres)
 740         enddo
 741       enddo
 742       do i=1,2*nres
 743         do j=1,3
 744           chomo(j,i,k)=c(j,i)
 745         enddo
 746       enddo
 747
 748       return
 749       end
 750
 751