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