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