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