source/unres/src_MD-M/readpdb.F

   1       subroutine readpdb
   2 C Read the PDB file and convert the peptide geometry into virtual-chain
   3 C geometry.
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6       include 'COMMON.LOCAL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.GEO'
  12       include 'COMMON.NAMES'
  13       include 'COMMON.CONTROL'
  14       include 'COMMON.DISTFIT'
  15       include 'COMMON.SETUP'
  16       integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
  17      &  ishift_pdb
  18       logical lprn /.true./,fail
  19       double precision e1(3),e2(3),e3(3)
  20       double precision dcj,efree_temp
  21       character*3 seq,res
  22       character*5 atom
  23       character*80 card
  24       double precision sccor(3,20)
  25       integer rescode
  26       efree_temp=0.0d0
  27       ibeg=1
  28       ishift1=0
  29       ishift=0
  30 c      write (2,*) "UNRES_PDB",unres_pdb
  31       ires=0
  32       ires_old=0
  33       nres=0
  34       iii=0
  35       lsecondary=.false.
  36       nhfrag=0
  37       nbfrag=0
  38       do i=1,100000
  39         read (ipdbin,'(a80)',end=10) card
  40 c        write (iout,'(a)') card
  41         if (card(:5).eq.'HELIX') then
  42          nhfrag=nhfrag+1
  43          lsecondary=.true.
  44          read(card(22:25),*) hfrag(1,nhfrag)
  45          read(card(34:37),*) hfrag(2,nhfrag)
  46         endif
  47         if (card(:5).eq.'SHEET') then
  48          nbfrag=nbfrag+1
  49          lsecondary=.true.
  50          read(card(24:26),*) bfrag(1,nbfrag)
  51          read(card(35:37),*) bfrag(2,nbfrag)
  52 crc----------------------------------------
  53 crc  to be corrected !!!
  54          bfrag(3,nbfrag)=bfrag(1,nbfrag)
  55          bfrag(4,nbfrag)=bfrag(2,nbfrag)
  56 crc----------------------------------------
  57         endif
  58         if (card(:3).eq.'END') then
  59           goto 10
  60         else if (card(:3).eq.'TER') then
  61 C End current chain
  62           ires_old=ires+1
  63           ishift1=ishift1+1
  64           itype(ires_old)=ntyp1
  65           ibeg=2
  66 c          write (iout,*) "Chain ended",ires,ishift,ires_old
  67           if (unres_pdb) then
  68             do j=1,3
  69               dc(j,ires)=sccor(j,iii)
  70             enddo
  71           else
  72             call sccenter(ires,iii,sccor)
  73           endif
  74           iii=0
  75         endif
  76 c Read free energy
  77         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
  78 C Fish out the ATOM cards.
  79         if (index(card(1:4),'ATOM').gt.0) then
  80           read (card(12:16),*) atom
  81 c          write (iout,*) "! ",atom," !",ires
  82 c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
  83           read (card(23:26),*) ires
  84           read (card(18:20),'(a3)') res
  85 c          write (iout,*) "ires",ires,ires-ishift+ishift1,
  86 c     &      " ires_old",ires_old
  87 c          write (iout,*) "ishift",ishift," ishift1",ishift1
  88 c          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
  89           if (ires-ishift+ishift1.ne.ires_old) then
  90 C Calculate the CM of the preceding residue.
  91 c            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
  92             if (ibeg.eq.0) then
  93 c              write (iout,*) "Calculating sidechain center iii",iii
  94               if (unres_pdb) then
  95                 do j=1,3
  96                   dc(j,ires+nres)=sccor(j,iii)
  97                 enddo
  98               else
  99                 call sccenter(ires_old,iii,sccor)
 100               endif
 101               iii=0
 102             endif
 103 C 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 c              write (iout,*) "ishift",ishift," ires",ires,
 117 c     &         " ires_old",ires_old
 118               ibeg=0
 119             else if (ibeg.eq.2) then
 120 c Start a new chain
 121 c              ishift=-ires_old+ires-1
 122 c              ishift1=ishift1+1
 123 c              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 c          write (iout,*) "ires_old",ires_old," ires",ires
 141           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 142 c            ishift1=ishift1+1
 143           endif
 144 c          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 c            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 (*,*) card(23:27),ires,itype(ires)
 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 c            write (iout,*) "sidechain ",atom
 163             iii=iii+1
 164             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
 165           endif
 166         endif
 167       enddo
 168    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
 169       if (ires.eq.0) return
 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         if (itype(i).eq.ntyp1) then
 176 c          write (iout,*) "dummy",i,itype(i)
 177           do j=1,3
 178             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
 179 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
 180             dc(j,i)=c(j,i)
 181           enddo
 182         endif
 183       enddo
 184 C Calculate the CM of the last side chain.
 185       if (iii.gt.0)  then
 186       if (unres_pdb) then
 187         do j=1,3
 188           dc(j,ires)=sccor(j,iii)
 189         enddo
 190       else
 191         call sccenter(ires,iii,sccor)
 192       endif
 193       endif
 194 c      nres=ires
 195       nsup=nres
 196       nstart_sup=1
 197       if (itype(nres).ne.10) then
 198         nres=nres+1
 199         itype(nres)=ntyp1
 200         if (unres_pdb) then
 201 C 2/15/2013 by Adam: corrected insertion of the last dummy residue
 202           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
 203           if (fail) then
 204             e2(1)=0.0d0
 205             e2(2)=1.0d0
 206             e2(3)=0.0d0
 207           endif
 208           do j=1,3
 209             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
 210           enddo
 211         else
 212         do j=1,3
 213           dcj=c(j,nres-2)-c(j,nres-3)
 214           c(j,nres)=c(j,nres-1)+dcj
 215           c(j,2*nres)=c(j,nres)
 216         enddo
 217         endif
 218       endif
 219       do i=2,nres-1
 220         do j=1,3
 221           c(j,i+nres)=dc(j,i)
 222         enddo
 223       enddo
 224       do j=1,3
 225         c(j,nres+1)=c(j,1)
 226         c(j,2*nres)=c(j,nres)
 227       enddo
 228       if (itype(1).eq.ntyp1) then
 229         nsup=nsup-1
 230         nstart_sup=2
 231         if (unres_pdb) then
 232 C 2/15/2013 by Adam: corrected insertion of the first dummy residue
 233           call refsys(2,3,4,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,1)=c(j,2)-3.8d0*e2(j)
 241           enddo
 242         else
 243         do j=1,3
 244           dcj=c(j,4)-c(j,3)
 245           c(j,1)=c(j,2)-dcj
 246           c(j,nres+1)=c(j,1)
 247         enddo
 248         endif
 249       endif
 250 C Copy the coordinates to reference coordinates
 251 c      do i=1,2*nres
 252 c        do j=1,3
 253 c          cref(j,i)=c(j,i)
 254 c        enddo
 255 c      enddo
 256 C Calculate internal coordinates.
 257       if (lprn) then
 258       write (iout,'(/a)')
 259      &  "Cartesian coordinates of the reference structure"
 260       write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
 261      & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
 262       do ires=1,nres
 263         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
 264      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
 265      &    (c(j,ires+nres),j=1,3)
 266       enddo
 267       endif
 268 C Calculate internal coordinates.
 269       if(me.eq.king.or..not.out1file)then
 270        write (iout,'(a)')
 271      &   "Backbone and SC coordinates as read from the PDB"
 272        do ires=1,nres
 273         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 274      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 275      &    (c(j,nres+ires),j=1,3)
 276        enddo
 277       endif
 278       call int_from_cart(.true.,.false.)
 279       call sc_loc_geom(.true.)
 280 c wczesbiej bylo false
 281       do i=1,nres
 282         thetaref(i)=theta(i)
 283         phiref(i)=phi(i)
 284       enddo
 285       do i=1,nres-1
 286         do j=1,3
 287           dc(j,i)=c(j,i+1)-c(j,i)
 288           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 289         enddo
 290       enddo
 291       do i=2,nres-1
 292         do j=1,3
 293           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 294           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 295         enddo
 296 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 297 c     &   vbld_inv(i+nres)
 298       enddo
 299 c      call chainbuild
 300 C Copy the coordinates to reference coordinates
 301 C Splits to single chain if occurs
 302
 303 c      do i=1,2*nres
 304 c        do j=1,3
 305 c          cref(j,i,cou)=c(j,i)
 306 c        enddo
 307 c      enddo
 308 c
 309       kkk=1
 310       lll=0
 311       cou=1
 312       do i=1,nres
 313       lll=lll+1
 314 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 315       if (i.gt.1) then
 316       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
 317       chain_length=lll-1
 318       kkk=kkk+1
 319 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 320       lll=1
 321       endif
 322       endif
 323         do j=1,3
 324           cref(j,i,cou)=c(j,i)
 325           cref(j,i+nres,cou)=c(j,i+nres)
 326           if (i.le.nres) then
 327           chain_rep(j,lll,kkk)=c(j,i)
 328           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
 329           endif
 330          enddo
 331       enddo
 332       write (iout,*) chain_length
 333       if (chain_length.eq.0) chain_length=nres
 334       do j=1,3
 335       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 336       chain_rep(j,chain_length+nres,symetr)
 337      &=chain_rep(j,chain_length+nres,1)
 338       enddo
 339 c diagnostic
 340 c       write (iout,*) "spraw lancuchy",chain_length,symetr
 341 c       do i=1,4
 342 c         do kkk=1,chain_length
 343 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
 344 c         enddo
 345 c        enddo
 346 c enddiagnostic
 347 C makes copy of chains
 348         write (iout,*) "symetr", symetr
 349
 350       if (symetr.gt.1) then
 351        call permut(symetr)
 352        nperm=1
 353        do i=1,symetr
 354        nperm=nperm*i
 355        enddo
 356        do i=1,nperm
 357        write(iout,*) (tabperm(i,kkk),kkk=1,4)
 358        enddo
 359        do i=1,nperm
 360         cou=0
 361         do kkk=1,symetr
 362          icha=tabperm(i,kkk)
 363 c         write (iout,*) i,icha
 364          do lll=1,chain_length
 365           cou=cou+1
 366            if (cou.le.nres) then
 367            do j=1,3
 368             kupa=mod(lll,chain_length)
 369             iprzes=(kkk-1)*chain_length+lll
 370             if (kupa.eq.0) kupa=chain_length
 371 c            write (iout,*) "kupa", kupa
 372             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 373             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 374           enddo
 375           endif
 376          enddo
 377         enddo
 378        enddo
 379        endif
 380 C-koniec robienia kopii
 381 c diag
 382       do kkk=1,nperm
 383       write (iout,*) "nowa struktura", nperm
 384       do i=1,nres
 385         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 386      &cref(2,i,kkk),
 387      &cref(3,i,kkk),cref(1,nres+i,kkk),
 388      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 389       enddo
 390   100 format (//'              alpha-carbon coordinates       ',
 391      &          '     centroid coordinates'/
 392      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 393      &                          10X,'X',11X,'Y',11X,'Z')
 394   110 format (a,'(',i3,')',6f12.5)
 395
 396       enddo
 397 cc enddiag
 398       do j=1,nbfrag
 399         do i=1,4
 400          bfrag(i,j)=bfrag(i,j)-ishift
 401         enddo
 402       enddo
 403
 404       do j=1,nhfrag
 405         do i=1,2
 406          hfrag(i,j)=hfrag(i,j)-ishift
 407         enddo
 408       enddo
 409       ishift_pdb=ishift
 410       return
 411       end
 412 c---------------------------------------------------------------------------
 413       subroutine int_from_cart(lside,lprn)
 414       implicit real*8 (a-h,o-z)
 415       include 'DIMENSIONS'
 416 #ifdef MPI
 417       include "mpif.h"
 418 #endif
 419       include 'COMMON.LOCAL'
 420       include 'COMMON.VAR'
 421       include 'COMMON.CHAIN'
 422       include 'COMMON.INTERACT'
 423       include 'COMMON.IOUNITS'
 424       include 'COMMON.GEO'
 425       include 'COMMON.NAMES'
 426       include 'COMMON.CONTROL'
 427       include 'COMMON.SETUP'
 428       character*3 seq,res
 429 c      character*5 atom
 430       character*80 card
 431       dimension sccor(3,20)
 432       integer rescode
 433       logical lside,lprn
 434       if(me.eq.king.or..not.out1file)then
 435        if (lprn) then
 436         write (iout,'(/a)')
 437      &  'Internal coordinates calculated from crystal structure.'
 438         if (lside) then
 439           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 440      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 441      & '     Beta '
 442         else
 443           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 444      & '     Gamma'
 445         endif
 446        endif
 447       endif
 448       do i=1,nres-1
 449         iti=itype(i)
 450         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 451           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 452 ctest          stop
 453         endif
 454         vbld(i+1)=dist(i,i+1)
 455         vbld_inv(i+1)=1.0d0/vbld(i+1)
 456         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 457         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 458       enddo
 459 c      if (unres_pdb) then
 460 c        if (itype(1).eq.21) then
 461 c          theta(3)=90.0d0*deg2rad
 462 c          phi(4)=180.0d0*deg2rad
 463 c          vbld(2)=3.8d0
 464 c          vbld_inv(2)=1.0d0/vbld(2)
 465 c        endif
 466 c        if (itype(nres).eq.21) then
 467 c          theta(nres)=90.0d0*deg2rad
 468 c          phi(nres)=180.0d0*deg2rad
 469 c          vbld(nres)=3.8d0
 470 c          vbld_inv(nres)=1.0d0/vbld(2)
 471 c        endif
 472 c      endif
 473       if (lside) then
 474         do i=2,nres-1
 475           do j=1,3
 476             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 477      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 478           enddo
 479           iti=itype(i)
 480           di=dist(i,nres+i)
 481 C 10/03/12 Adam: Correction for zero SC-SC bond length
 482           if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
 483      &     di=dsc(itype(i))
 484           vbld(i+nres)=di
 485           if (itype(i).ne.10) then
 486             vbld_inv(i+nres)=1.0d0/di
 487           else
 488             vbld_inv(i+nres)=0.0d0
 489           endif
 490           if (iti.ne.10) then
 491             alph(i)=alpha(nres+i,i,maxres2)
 492             omeg(i)=beta(nres+i,i,maxres2,i+1)
 493           endif
 494           if(me.eq.king.or..not.out1file)then
 495            if (lprn)
 496      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 497      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 498      &     rad2deg*alph(i),rad2deg*omeg(i)
 499           endif
 500         enddo
 501       else if (lprn) then
 502         do i=2,nres
 503           iti=itype(i)
 504           if(me.eq.king.or..not.out1file)
 505      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 506      &     rad2deg*theta(i),rad2deg*phi(i)
 507         enddo
 508       endif
 509       return
 510       end
 511 c-------------------------------------------------------------------------------
 512       subroutine sc_loc_geom(lprn)
 513       implicit real*8 (a-h,o-z)
 514       include 'DIMENSIONS'
 515 #ifdef MPI
 516       include "mpif.h"
 517 #endif
 518       include 'COMMON.LOCAL'
 519       include 'COMMON.VAR'
 520       include 'COMMON.CHAIN'
 521       include 'COMMON.INTERACT'
 522       include 'COMMON.IOUNITS'
 523       include 'COMMON.GEO'
 524       include 'COMMON.NAMES'
 525       include 'COMMON.CONTROL'
 526       include 'COMMON.SETUP'
 527       double precision x_prime(3),y_prime(3),z_prime(3)
 528       logical lprn
 529       do i=1,nres-1
 530         do j=1,3
 531           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 532         enddo
 533       enddo
 534       do i=2,nres-1
 535         if (itype(i).ne.10) then
 536           do j=1,3
 537             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 538           enddo
 539         else
 540           do j=1,3
 541             dc_norm(j,i+nres)=0.0d0
 542           enddo
 543         endif
 544       enddo
 545       do i=2,nres-1
 546         costtab(i+1) =dcos(theta(i+1))
 547         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 548         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 549         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 550         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 551         cosfac=dsqrt(cosfac2)
 552         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 553         sinfac=dsqrt(sinfac2)
 554         it=itype(i)
 555         if (it.ne.10) then
 556 c
 557 C  Compute the axes of tghe local cartesian coordinates system; store in
 558 c   x_prime, y_prime and z_prime
 559 c
 560         do j=1,3
 561           x_prime(j) = 0.00
 562           y_prime(j) = 0.00
 563           z_prime(j) = 0.00
 564         enddo
 565         do j = 1,3
 566           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 567           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 568         enddo
 569         call vecpr(x_prime,y_prime,z_prime)
 570 c
 571 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 572 C to local coordinate system. Store in xx, yy, zz.
 573 c
 574         xx=0.0d0
 575         yy=0.0d0
 576         zz=0.0d0
 577         do j = 1,3
 578           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 579           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 580           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 581         enddo
 582
 583         xxref(i)=xx
 584         yyref(i)=yy
 585         zzref(i)=zz
 586         else
 587         xxref(i)=0.0d0
 588         yyref(i)=0.0d0
 589         zzref(i)=0.0d0
 590         endif
 591       enddo
 592       if (lprn) then
 593         do i=2,nres
 594           iti=itype(i)
 595           if(me.eq.king.or..not.out1file)
 596      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 597      &      yyref(i),zzref(i)
 598         enddo
 599       endif
 600       return
 601       end
 602 c---------------------------------------------------------------------------
 603       subroutine sccenter(ires,nscat,sccor)
 604       implicit real*8 (a-h,o-z)
 605       include 'DIMENSIONS'
 606       include 'COMMON.CHAIN'
 607       dimension sccor(3,20)
 608       do j=1,3
 609         sccmj=0.0D0
 610         do i=1,nscat
 611           sccmj=sccmj+sccor(j,i)
 612         enddo
 613         dc(j,ires)=sccmj/nscat
 614       enddo
 615       return
 616       end
 617 c---------------------------------------------------------------------------
 618       subroutine bond_regular
 619       implicit real*8 (a-h,o-z)
 620       include 'DIMENSIONS'
 621       include 'COMMON.VAR'
 622       include 'COMMON.LOCAL'
 623       include 'COMMON.CALC'
 624       include 'COMMON.INTERACT'
 625       include 'COMMON.CHAIN'
 626       do i=1,nres-1
 627        vbld(i+1)=vbl
 628        vbld_inv(i+1)=1.0d0/vbld(i+1)
 629        vbld(i+1+nres)=dsc(itype(i+1))
 630        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 631 c       print *,vbld(i+1),vbld(i+1+nres)
 632       enddo
 633       return
 634       end