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