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