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