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