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