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