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