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