source/wham/src-HCD/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 'DIMENSIONS.ZSCOPT'
   7       include 'COMMON.CONTROL'
   8       include 'COMMON.LOCAL'
   9       include 'COMMON.VAR'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.INTERACT'
  12       include 'COMMON.IOUNITS'
  13       include 'COMMON.GEO'
  14       include 'COMMON.NAMES'
  15       include 'COMMON.SBRIDGE'
  16       include 'COMMON.FRAG'
  17       character*3 seq,atom,res
  18       character*80 card
  19       double precision e1(3),e2(3),e3(3)
  20       double precision sccor(3,50)
  21       integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
  22       double precision dcj
  23       integer rescode,kkk,lll,icha,cou,kupa,iprzes
  24       logical lsecondary,sccalc,fail,zero
  25       integer iterter(maxres)
  26       double precision efree_temp
  27       iii=0
  28       ibeg=1
  29       ishift1=0
  30       sccalc=.false.
  31       bfac=0.0d0
  32       do i=1,maxres
  33          iterter(i)=0
  34       enddo
  35       ibeg=1
  36       ishift1=0
  37       lsecondary=.false.
  38       nhfrag=0
  39       nbfrag=0
  40       iii=0
  41       sccalc=.false.
  42       do
  43         read (ipdbin,'(a80)',end=10) card
  44 c        write (iout,'(a)') card
  45 c        call flush(iout)
  46         if (card(:5).eq.'HELIX') then
  47           nhfrag=nhfrag+1
  48           lsecondary=.true.
  49           read(card(22:25),*) hfrag(1,nhfrag)
  50           read(card(34:37),*) hfrag(2,nhfrag)
  51         endif
  52         if (card(:5).eq.'SHEET') then
  53           nbfrag=nbfrag+1
  54           lsecondary=.true.
  55           read(card(24:26),*) bfrag(1,nbfrag)
  56           read(card(35:37),*) bfrag(2,nbfrag)
  57 !rc----------------------------------------
  58 !rc  to be corrected !!!
  59           bfrag(3,nbfrag)=bfrag(1,nbfrag)
  60           bfrag(4,nbfrag)=bfrag(2,nbfrag)
  61 !rc----------------------------------------
  62         endif
  63         if (card(:3).eq.'END') then
  64           goto 10
  65         else if (card(:3).eq.'TER') then
  66 ! End current chain
  67           ires_old=ires+2
  68           itype(ires_old-1)=ntyp1
  69           iterter(ires_old-1)=1
  70           itype(ires_old)=ntyp1
  71           iterter(ires_old)=1
  72 c          ishift1=ishift1+1
  73           ibeg=2
  74           write (iout,*) "Chain ended",ires,ishift,ires_old,ibeg
  75           if (unres_pdb) then
  76             do j=1,3
  77               dc(j,ires)=sccor(j,iii)
  78             enddo
  79           else
  80             call sccenter(ires,iii,sccor)
  81           endif
  82           iii=0
  83           sccalc=.true.
  84         endif
  85 ! Read free energy
  86 c        if (index(card,"FREE ENERGY").gt.0) then
  87 c          ifree=index(card,"FREE ENERGY")+12
  88 c          read(card(ifree:),*,err=1115,end=1115) efree_temp
  89 c 1115     continue
  90 c        endif
  91 ! Fish out the ATOM cards.
  92         if (index(card(1:4),'ATOM').gt.0) then
  93           sccalc=.false.
  94           read (card(12:16),*) atom
  95 c          write (2,'(a)') card
  96 c          write (iout,*) "ibeg",ibeg
  97 c          write (iout,*) "! ",atom," !",ires
  98 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
  99           read (card(23:26),*) ires
 100           read (card(18:20),'(a3)') res
 101 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
 102 c     &      " ires_old",ires_old
 103 c         write (iout,*) "ishift",ishift," ishift1",ishift1
 104 c         write (iout,*) "IRES",ires-ishift+ishift1,ires_old
 105           if (ires-ishift+ishift1.ne.ires_old) then
 106 ! Calculate the CM of the preceding residue.
 107 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
 108             if (ibeg.eq.0) then
 109 c              write (iout,*) "Calculating sidechain center iii",iii
 110 c              write (iout,*) "ires",ires
 111               if (unres_pdb) then
 112 c                write (iout,'(i5,3f10.5)') ires,(sccor(j,iii),j=1,3)
 113                 do j=1,3
 114                   dc(j,ires_old)=sccor(j,iii)
 115                 enddo
 116               else
 117                 call sccenter(ires_old,iii,sccor)
 118               endif
 119               iii=0
 120               sccalc=.true.
 121             endif
 122 ! Start new residue.
 123 c            write (iout,*) "ibeg",ibeg
 124             if (res.eq.'Cl-' .or. res.eq.'Na+') then
 125               ires=ires_old
 126               cycle
 127             else if (ibeg.eq.1) then
 128 c              write (iout,*) "BEG ires",ires
 129               ishift=ires-1
 130               if (res.ne.'GLY' .and. res.ne. 'ACE') then
 131                 ishift=ishift-1
 132                 itype(1)=ntyp1
 133               endif
 134               ires=ires-ishift+ishift1
 135               ires_old=ires
 136 !              write (iout,*) "ishift",ishift," ires",ires,&
 137 !               " ires_old",ires_old
 138               ibeg=0
 139             else if (ibeg.eq.2) then
 140 ! Start a new chain
 141               ishift=-ires_old+ires-1 !!!!!
 142 c              ishift1=ishift1-1    !!!!!
 143 c              write (iout,*) "New chain started",ires,ires_old,ishift,
 144 c     &           ishift1
 145               ires=ires-ishift+ishift1
 146               write (iout,*) "New chain started ires",ires
 147               ires_old=ires
 148 c              ires=ires_old+1
 149               ibeg=0
 150             else
 151               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
 152               ires=ires-ishift+ishift1
 153               ires_old=ires
 154             endif
 155             if (res.eq.'ACE' .or. res.eq.'NHE') then
 156               itype(ires)=10
 157             else
 158               itype(ires)=rescode(ires,res,0)
 159             endif
 160           else
 161             ires=ires-ishift+ishift1
 162           endif
 163 c          write (iout,*) "ires_old",ires_old," ires",ires
 164           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 165 !            ishift1=ishift1+1
 166           endif
 167 c          write (2,*) "ires",ires," res ",res!," ity"!,ity
 168           if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
 169      &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
 170             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
 171 c            write (iout,*) "backbone ",atom
 172 c            write (iout,*) ires,res,(c(j,ires),j=1,3)
 173 #ifdef DEBUG
 174             write (iout,'(i6,i3,2x,a,3f8.3)')
 175      &      ires,itype(ires),res,(c(j,ires),j=1,3)
 176 #endif
 177             iii=iii+1
 178             do j=1,3
 179               sccor(j,iii)=c(j,ires)
 180             enddo
 181 c            write (2,*) card(23:27),ires,itype(ires),iii
 182           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
 183      &             atom.ne.'N' .and. atom.ne.'C' .and.
 184      &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
 185      &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
 186 !            write (iout,*) "sidechain ",atom
 187             iii=iii+1
 188             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 189 c            write (2,*) "iii",iii
 190           endif
 191         endif
 192       enddo
 193    10 write (iout,'(a,i5)') ' Nres: ',ires
 194 c      write (iout,*) "iii",iii
 195 C Calculate dummy residue coordinates inside the "chain" of a multichain
 196 C system
 197       nres=ires
 198 c      write (iout,*) "dc"
 199 c      do i=1,nres
 200 c        write (iout,'(i5,3f10.5)') i,(dc(j,i),j=1,3)
 201 c      enddo
 202       do i=2,nres-1
 203 c        write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
 204         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 205          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 206 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 207 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 208 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 209            if (unres_pdb) then
 210 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 211 c            print *,i,'tu dochodze'
 212             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 213             if (fail) then
 214               e2(1)=0.0d0
 215               e2(2)=1.0d0
 216               e2(3)=0.0d0
 217             endif !fail
 218 c            print *,i,'a tu?'
 219             do j=1,3
 220              c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 221             enddo
 222            else   !unres_pdb
 223            do j=1,3
 224              dcj=(c(j,i-2)-c(j,i-3))/2.0
 225             if (dcj.eq.0) dcj=1.23591524223
 226              c(j,i)=c(j,i-1)+dcj
 227              c(j,nres+i)=c(j,i)
 228              dC(j,i)=c(j,i)
 229            enddo
 230           endif   !unres_pdb
 231          else     !itype(i+1).eq.ntyp1
 232           if (unres_pdb) then
 233 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 234             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 235             if (fail) then
 236               e2(1)=0.0d0
 237               e2(2)=1.0d0
 238               e2(3)=0.0d0
 239             endif
 240             do j=1,3
 241               c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 242             enddo
 243           else !unres_pdb
 244            do j=1,3
 245             dcj=(c(j,i+3)-c(j,i+2))/2.0
 246             if (dcj.eq.0) dcj=1.23591524223
 247             c(j,i)=c(j,i+1)-dcj
 248             c(j,nres+i)=c(j,i)
 249             dC(j,i)=c(j,i)
 250            enddo
 251           endif !unres_pdb
 252          endif !itype(i+1).eq.ntyp1
 253         endif  !itype.eq.ntyp1
 254       enddo
 255       write (iout,*) "After loop in readpbd"
 256 C Calculate the CM of the last side chain.
 257       if (.not. sccalc) then
 258       if (unres_pdb) then
 259         do j=1,3
 260           dc(j,ires)=sccor(j,iii)
 261         enddo
 262       else
 263 c        write (iout,*) "Calling sccenter iii",iii
 264         call sccenter(ires,iii,sccor)
 265       endif
 266       endif
 267       nsup=nres
 268       nstart_sup=1
 269       if (itype(nres).ne.10) then
 270         nres=nres+1
 271         itype(nres)=ntyp1
 272         if (unres_pdb) then
 273 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 274           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 275           if (fail) then
 276             e2(1)=0.0d0
 277             e2(2)=1.0d0
 278             e2(3)=0.0d0
 279           endif
 280           do j=1,3
 281             c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 282           enddo
 283         else
 284         do j=1,3
 285           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 286             if (dcj.eq.0) dcj=1.23591524223
 287           c(j,nres)=c(j,nres-1)+dcj
 288           c(j,2*nres)=c(j,nres)
 289         enddo
 290         endif
 291       endif
 292       do i=2,nres-1
 293         do j=1,3
 294           c(j,i+nres)=dc(j,i)
 295         enddo
 296       enddo
 297       do j=1,3
 298         c(j,nres+1)=c(j,1)
 299         c(j,2*nres)=c(j,nres)
 300       enddo
 301       if (itype(1).eq.ntyp1) then
 302         nsup=nsup-1
 303         nstart_sup=2
 304         if (unres_pdb) then
 305 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 306           call refsys(2,3,4,e1,e2,e3,fail)
 307           if (fail) then
 308             e2(1)=0.0d0
 309             e2(2)=1.0d0
 310             e2(3)=0.0d0
 311           endif
 312           do j=1,3
 313             c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
 314           enddo
 315         else
 316         do j=1,3
 317           dcj=(c(j,4)-c(j,3))/2.0
 318           c(j,1)=c(j,2)-dcj
 319           c(j,nres+1)=c(j,1)
 320         enddo
 321         endif
 322       endif
 323 C Calculate internal coordinates.
 324       write (iout,'(/a)')
 325      &  "Cartesian coordinates of the reference structure"
 326       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 327      & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 328       do ires=1,nres
 329         write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')
 330      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 331      &    (c(j,ires+nres),j=1,3)
 332       zero=.false.
 333       enddo
 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 real*8 (a-h,o-z)
 402       include 'DIMENSIONS'
 403       include 'DIMENSIONS.ZSCOPT'
 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.SETUP'
 413       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
 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,20)
 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 write (iout,'(a,i5)') ' Nres: ',ires
 550 C Calculate dummy residue coordinates inside the "chain" of a multichain
 551 C system
 552       nres=ires
 553       do i=2,nres-1
 554 c        write (iout,*) i,itype(i),itype(i+1)
 555         if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
 556          if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
 557 C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 558 C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
 559 C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
 560            if (unres_pdb) then
 561 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 562             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
 563             if (fail) then
 564               e2(1)=0.0d0
 565               e2(2)=1.0d0
 566               e2(3)=0.0d0
 567             endif !fail
 568             do j=1,3
 569              c(j,i)=c(j,i-1)-1.9d0*e2(j)
 570             enddo
 571            else   !unres_pdb
 572            do j=1,3
 573              dcj=(c(j,i-2)-c(j,i-3))/2.0
 574             if (dcj.eq.0) dcj=1.23591524223
 575              c(j,i)=c(j,i-1)+dcj
 576              c(j,nres+i)=c(j,i)
 577            enddo
 578           endif   !unres_pdb
 579          else     !itype(i+1).eq.ntyp1
 580           if (unres_pdb) then
 581 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 582             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
 583             if (fail) then
 584               e2(1)=0.0d0
 585               e2(2)=1.0d0
 586               e2(3)=0.0d0
 587             endif
 588             do j=1,3
 589               c(j,i)=c(j,i+1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 590             enddo
 591           else !unres_pdb
 592            do j=1,3
 593             dcj=(c(j,i+3)-c(j,i+2))/2.0
 594             if (dcj.eq.0) dcj=1.23591524223
 595             c(j,i)=c(j,i+1)-dcj
 596             c(j,nres+i)=c(j,i)
 597            enddo
 598           endif !unres_pdb
 599          endif !itype(i+1).eq.ntyp1
 600         endif  !itype.eq.ntyp1
 601       enddo
 602 C Calculate the CM of the last side chain.
 603       if (unres_pdb) then
 604         do j=1,3
 605           dc(j,ires)=sccor(j,iii)
 606         enddo
 607       else
 608         call sccenter(ires,iii,sccor)
 609       endif
 610       nsup=nres
 611       nstart_sup=1
 612       if (itype(nres).ne.10) then
 613         nres=nres+1
 614         itype(nres)=ntyp1
 615         if (unres_pdb) then
 616 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 617           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 618           if (fail) then
 619             e2(1)=0.0d0
 620             e2(2)=1.0d0
 621             e2(3)=0.0d0
 622           endif
 623           do j=1,3
 624             c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
 625           enddo
 626         else
 627         do j=1,3
 628           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 629             if (dcj.eq.0) dcj=1.23591524223
 630           c(j,nres)=c(j,nres-1)+dcj
 631           c(j,2*nres)=c(j,nres)
 632         enddo
 633       endif
 634       endif
 635       do i=2,nres-1
 636         do j=1,3
 637           c(j,i+nres)=dc(j,i)
 638         enddo
 639       enddo
 640       do j=1,3
 641         c(j,nres+1)=c(j,1)
 642         c(j,2*nres)=c(j,nres)
 643       enddo
 644       if (itype(1).eq.ntyp1) then
 645         nsup=nsup-1
 646         nstart_sup=2
 647         if (unres_pdb) then
 648 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 649           call refsys(2,3,4,e1,e2,e3,fail)
 650           if (fail) then
 651             e2(1)=0.0d0
 652             e2(2)=1.0d0
 653             e2(3)=0.0d0
 654           endif
 655           do j=1,3
 656             c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
 657           enddo
 658         else
 659         do j=1,3
 660           dcj=(c(j,4)-c(j,3))/2.0
 661           c(j,1)=c(j,2)-dcj
 662           c(j,nres+1)=c(j,1)
 663         enddo
 664         endif
 665       endif
 666 C Copy the coordinates to reference coordinates
 667 c      do i=1,2*nres
 668 c        do j=1,3
 669 c          cref(j,i)=c(j,i)
 670 c        enddo
 671 c      enddo
 672 C Calculate internal coordinates.
 673       if (out_template_coord) then
 674       write (iout,'(/a)')
 675      &  "Cartesian coordinates of the reference structure"
 676       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 677      & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 678       do ires=1,nres
 679         write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')
 680      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 681      &    (c(j,ires+nres),j=1,3)
 682       enddo
 683       endif
 684       zero=.false.
 685       do ires=1,nres
 686         zero=zero.or.itype(ires).eq.0
 687       enddo
 688       if (zero) then
 689         write (iout,'(2a)') "Gaps in PDB coordinates detected;",
 690      &    " look for ZERO in the control output above."
 691         write (iout,'(2a)') "Repair the PDB file using MODELLER",
 692      &    " or other softwared and resubmit."
 693         call flush(iout)
 694         stop
 695       endif
 696 C Calculate internal coordinates.
 697       call int_from_cart(.true.,out_template_coord)
 698       call sc_loc_geom(.false.)
 699       do i=1,nres
 700         thetaref(i)=theta(i)
 701         phiref(i)=phi(i)
 702       enddo
 703       dc(:,0)=c(:,1)
 704       do i=1,nres-1
 705         do j=1,3
 706           dc(j,i)=c(j,i+1)-c(j,i)
 707           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 708         enddo
 709       enddo
 710       do i=2,nres-1
 711         do j=1,3
 712           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 713           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 714         enddo
 715 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 716 c     &   vbld_inv(i+nres)
 717       enddo
 718       do i=1,nres
 719         do j=1,3
 720           cref(j,i)=c(j,i)
 721           cref(j,i+nres)=c(j,i+nres)
 722         enddo
 723       enddo
 724       do i=1,2*nres
 725         do j=1,3
 726           chomo(j,i,k)=c(j,i)
 727         enddo
 728       enddo
 729
 730       return
 731       end
 732 c---------------------------------------------------------------------------
 733       subroutine int_from_cart(lside,lprn)
 734       implicit none
 735       include 'DIMENSIONS'
 736       include 'DIMENSIONS.ZSCOPT'
 737       include 'COMMON.LOCAL'
 738       include 'COMMON.VAR'
 739       include 'COMMON.CHAIN'
 740       include 'COMMON.INTERACT'
 741       include 'COMMON.IOUNITS'
 742       include 'COMMON.GEO'
 743       include 'COMMON.NAMES'
 744       character*3 seq,atom,res
 745       character*80 card
 746       double precision sccor(3,50)
 747       integer rescode
 748       double precision dist,alpha,beta,di
 749       integer i,j,iti
 750       logical lside,lprn
 751       if (lprn) then
 752         write (iout,'(/a)')
 753      &  'Internal coordinates calculated from crystal structure.'
 754         if (lside) then
 755           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 756      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 757      & '     Omega'
 758         else
 759           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 760      & '       Phi'
 761         endif
 762       endif
 763       do i=2,nres
 764         iti=itype(i)
 765 c        write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
 766         if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
 767      &    (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
 768           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 769 c          stop
 770         endif
 771         vbld(i)=dist(i-1,i)
 772         vbld_inv(i)=1.0d0/vbld(i)
 773         theta(i+1)=alpha(i-1,i,i+1)
 774         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 775       enddo
 776 c      if (itype(1).eq.ntyp1) then
 777 c        do j=1,3
 778 c          c(j,1)=c(j,2)+(c(j,3)-c(j,4))
 779 c        enddo
 780 c      endif
 781 c      if (itype(nres).eq.ntyp1) then
 782 c        do j=1,3
 783 c          c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
 784 c        enddo
 785 c      endif
 786       if (lside) then
 787         do i=2,nres-1
 788           do j=1,3
 789             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 790           enddo
 791           iti=itype(i)
 792           di=dist(i,nres+i)
 793            vbld(i+nres)=di
 794           if (itype(i).ne.10) then
 795             vbld_inv(i+nres)=1.0d0/di
 796           else
 797             vbld_inv(i+nres)=0.0d0
 798           endif
 799           if (iti.ne.10) then
 800             alph(i)=alpha(nres+i,i,maxres2)
 801             omeg(i)=beta(nres+i,i,maxres2,i+1)
 802           endif
 803           if (iti.ne.10) then
 804             alph(i)=alpha(nres+i,i,maxres2)
 805             omeg(i)=beta(nres+i,i,maxres2,i+1)
 806           endif
 807           if (lprn)
 808      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 809      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
 810      &    rad2deg*alph(i),rad2deg*omeg(i)
 811         enddo
 812       else if (lprn) then
 813         do i=2,nres
 814           iti=itype(i)
 815           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 816      &    rad2deg*theta(i),rad2deg*phi(i)
 817         enddo
 818       endif
 819       return
 820       end
 821 c---------------------------------------------------------------------------
 822       subroutine sccenter(ires,nscat,sccor)
 823       implicit none
 824       include 'DIMENSIONS'
 825       include 'COMMON.CHAIN'
 826       integer ires,nscat,i,j
 827       double precision sccor(3,50),sccmj
 828       do j=1,3
 829         sccmj=0.0D0
 830         do i=1,nscat
 831           sccmj=sccmj+sccor(j,i)
 832         enddo
 833         dc(j,ires)=sccmj/nscat
 834       enddo
 835       return
 836       end
 837 c---------------------------------------------------------------------------
 838       subroutine sc_loc_geom(lprn)
 839       implicit real*8 (a-h,o-z)
 840       include 'DIMENSIONS'
 841       include 'DIMENSIONS.ZSCOPT'
 842       include 'COMMON.LOCAL'
 843       include 'COMMON.VAR'
 844       include 'COMMON.CHAIN'
 845       include 'COMMON.INTERACT'
 846       include 'COMMON.IOUNITS'
 847       include 'COMMON.GEO'
 848       include 'COMMON.NAMES'
 849       include 'COMMON.CONTROL'
 850       include 'COMMON.SETUP'
 851       double precision x_prime(3),y_prime(3),z_prime(3)
 852       logical lprn
 853       do i=1,nres-1
 854         do j=1,3
 855           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 856         enddo
 857       enddo
 858       do i=2,nres-1
 859         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 860           do j=1,3
 861             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 862           enddo
 863         else
 864           do j=1,3
 865             dc_norm(j,i+nres)=0.0d0
 866           enddo
 867         endif
 868       enddo
 869       do i=2,nres-1
 870         costtab(i+1) =dcos(theta(i+1))
 871         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 872         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 873         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 874         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 875         cosfac=dsqrt(cosfac2)
 876         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 877         sinfac=dsqrt(sinfac2)
 878         it=itype(i)
 879         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 880 c
 881 C  Compute the axes of tghe local cartesian coordinates system; store in
 882 c   x_prime, y_prime and z_prime
 883 c
 884         do j=1,3
 885           x_prime(j) = 0.00
 886           y_prime(j) = 0.00
 887           z_prime(j) = 0.00
 888         enddo
 889         do j = 1,3
 890           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 891           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 892         enddo
 893 c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
 894 c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
 895         call vecpr(x_prime,y_prime,z_prime)
 896 c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
 897 c
 898 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 899 C to local coordinate system. Store in xx, yy, zz.
 900 c
 901         xx=0.0d0
 902         yy=0.0d0
 903         zz=0.0d0
 904         do j = 1,3
 905           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 906           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 907           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 908         enddo
 909
 910         xxref(i)=xx
 911         yyref(i)=yy
 912         zzref(i)=zz
 913         else
 914         xxref(i)=0.0d0
 915         yyref(i)=0.0d0
 916         zzref(i)=0.0d0
 917         endif
 918       enddo
 919       if (lprn) then
 920         write (iout,*) "xxref,yyref,zzref"
 921         do i=2,nres
 922           iti=itype(i)
 923           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
 924      &     zzref(i)
 925         enddo
 926       endif
 927       return
 928       end
 929 c---------------------------------------------------------------------------
 930       subroutine bond_regular
 931       implicit none
 932       include 'DIMENSIONS'
 933       include 'COMMON.VAR'
 934       include 'COMMON.LOCAL'
 935       include 'COMMON.INTERACT'
 936       include 'COMMON.CHAIN'
 937       integer i,i1,i2
 938       do i=1,nres-1
 939        vbld(i+1)=vbl
 940        vbld_inv(i+1)=vblinv
 941        vbld(i+1+nres)=dsc(iabs(itype(i+1)))
 942        vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
 943 c       print *,vbld(i+1),vbld(i+1+nres)
 944       enddo
 945 c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
 946       do i=1,nchain
 947         i1=chain_border(1,i)
 948         i2=chain_border(2,i)
 949         if (i1.gt.1) then
 950           vbld(i1)=vbld(i1)/2
 951           vbld_inv(i1)=vbld_inv(i1)*2
 952         endif
 953         if (i2.lt.nres) then
 954           vbld(i2+1)=vbld(i2+1)/2
 955           vbld_inv(i2+1)=vbld_inv(i2+1)*2
 956         endif
 957       enddo
 958       return
 959       end
 960