source/cluster/wham/src-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       character*3 seq,atom,res
  14       character*80 card
  15       dimension sccor(3,20)
  16       integer rescode
  17 c      call permut(symetr)
  18       ibeg=1
  19       write(iout,*) 'pdbread'
  20       do
  21         read (ipdbin,'(a80)',end=10) card
  22         if (card(:3).eq.'END') then
  23           goto 10
  24         else if (card(:3).eq.'TER') then
  25 C End current chain
  26           ires_old=ires+2
  27           itype(ires_old-1)=ntyp1
  28           itype(ires_old)=ntyp1
  29           ibeg=2
  30 c          write (iout,*) "Chain ended",ires,ishift,ires_old
  31           call sccenter(ires,iii,sccor)
  32         endif
  33 C Fish out the ATOM cards.
  34         if (index(card(1:4),'ATOM').gt.0) then
  35           read (card(14:16),'(a3)') atom
  36           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  37 C Calculate the CM of the preceding residue.
  38             if (ibeg.eq.0) then
  39               call sccenter(ires,iii,sccor)
  40             endif
  41 C Start new residue.
  42 c            write (iout,'(a80)') card
  43             read (card(24:26),*) ires
  44             read (card(18:20),'(a3)') res
  45             if (ibeg.eq.1) then
  46               ishift=ires-1
  47               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  48                 ishift=ishift-1
  49                 itype(1)=ntyp1
  50               endif
  51 c              write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
  52               ibeg=0
  53             else if (ibeg.eq.2) then
  54 c Start a new chain
  55               ishift=-ires_old+ires-1
  56 c              write (iout,*) "New chain started",ires,ishift
  57               ibeg=0
  58             endif
  59             ires=ires-ishift
  60 c            write (2,*) "ires",ires," ishift",ishift
  61             if (res.eq.'ACE') then
  62               ity=10
  63             else
  64               itype(ires)=rescode(ires,res,0)
  65             endif
  66             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  67             write (iout,'(2i3,2x,a,3f8.3)')
  68      &       ires,itype(ires),res,(c(j,ires),j=1,3)
  69             iii=1
  70             do j=1,3
  71               sccor(j,iii)=c(j,ires)
  72             enddo
  73           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  74      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  75             iii=iii+1
  76             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  77           endif
  78         endif
  79       enddo
  80    10 write (iout,'(a,i5)') ' Nres: ',ires
  81 C Calculate dummy residue coordinates inside the "chain" of a multichain
  82 C system
  83       nres=ires
  84       do i=2,nres-1
  85 c        write (iout,*) i,itype(i)
  86         if (itype(i).eq.ntyp1) then
  87          if (itype(i+1).eq.ntyp1) then
  88
  89 c          write (iout,*) "dummy",i,itype(i)
  90 C          do j=1,3
  91 C            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
  92 c            c(j,i)=(c(j,i-1)+c(j,i+1))/2
  93 C            dc(j,i)=c(j,i)
  94 C          enddo
  95            do j=1,3
  96              dcj=(c(j,i-2)-c(j,i-3))/2.0
  97              c(j,i)=c(j,i-1)+dcj
  98              c(j,nres+i)=c(j,i)
  99            enddo
 100 C          endif   !unres_pdb
 101          else     !itype(i+1).eq.ntyp1
 102            do j=1,3
 103             dcj=(c(j,i+3)-c(j,i+2))/2.0
 104             c(j,i)=c(j,i+1)-dcj
 105             c(j,nres+i)=c(j,i)
 106            enddo
 107 C          endif !unres_pdb
 108          endif !itype(i+1).eq.ntyp1
 109         endif  !itype.eq.ntyp1
 110       enddo
 111 C Calculate the CM of the last side chain.
 112       call sccenter(ires,iii,sccor)
 113       nsup=nres
 114       nstart_sup=1
 115       if (itype(nres).ne.10) then
 116         nres=nres+1
 117         itype(nres)=ntyp1
 118         do j=1,3
 119           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
 120           c(j,nres)=c(j,nres-1)+dcj
 121           c(j,2*nres)=c(j,nres)
 122         enddo
 123       endif
 124       do i=2,nres-1
 125         do j=1,3
 126           c(j,i+nres)=dc(j,i)
 127         enddo
 128       enddo
 129       do j=1,3
 130         c(j,nres+1)=c(j,1)
 131         c(j,2*nres)=c(j,nres)
 132       enddo
 133       if (itype(1).eq.ntyp1) then
 134         nsup=nsup-1
 135         nstart_sup=2
 136         do j=1,3
 137           dcj=(c(j,4)-c(j,3))/2.0
 138           c(j,1)=c(j,2)-dcj
 139           c(j,nres+1)=c(j,1)
 140         enddo
 141       endif
 142 C Calculate internal coordinates.
 143       do ires=1,nres
 144         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 145      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 146      &    (c(j,nres+ires),j=1,3)
 147       enddo
 148       call int_from_cart1(.false.)
 149       call int_from_cart(.true.,.false.)
 150       call sc_loc_geom(.true.)
 151       write (iout,*) "After int_from_cart"
 152       call flush(iout)
 153       do i=1,nres-1
 154         do j=1,3
 155           dc(j,i)=c(j,i+1)-c(j,i)
 156           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 157         enddo
 158       enddo
 159       do i=2,nres-1
 160         do j=1,3
 161           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 162           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 163         enddo
 164 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 165 c     &   vbld_inv(i+nres)
 166       enddo
 167       do i=1,nres
 168         thetaref(i)=theta(i)
 169         phiref(i)=phi(i)
 170 c
 171         phi_ref(i)=phi(i)
 172         theta_ref(i)=theta(i)
 173         alph_ref(i)=alph(i)
 174         omeg_ref(i)=omeg(i)
 175       enddo
 176
 177 c      call chainbuild
 178 C Copy the coordinates to reference coordinates
 179 c      do i=1,2*nres
 180 c        do j=1,3
 181 c          cref(j,i)=c(j,i)
 182 c        enddo
 183 c      enddo
 184
 185       kkk=1
 186       lll=0
 187       cou=1
 188       do i=1,nres
 189       lll=lll+1
 190 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 191       if (i.gt.1) then
 192       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
 193       chain_length=lll-1
 194       kkk=kkk+1
 195 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 196       lll=1
 197       endif
 198       endif
 199         do j=1,3
 200           cref(j,i,cou)=c(j,i)
 201           cref(j,i+nres,cou)=c(j,i+nres)
 202           if (i.le.nres) then
 203           chain_rep(j,lll,kkk)=c(j,i)
 204           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
 205           endif
 206          enddo
 207       enddo
 208       do j=1,3
 209       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 210       chain_rep(j,chain_length+nres,symetr)
 211      &=chain_rep(j,chain_length+nres,1)
 212       enddo
 213
 214       if (symetr.gt.1) then
 215        call permut(symetr)
 216        nperm=1
 217        do i=1,symetr
 218        nperm=nperm*i
 219        enddo
 220 c       do i=1,nperm
 221 c       write(iout,*) "tabperm", (tabperm(i,kkk),kkk=1,4)
 222 c       enddo
 223        do i=1,nperm
 224         cou=0
 225         do kkk=1,symetr
 226          icha=tabperm(i,kkk)
 227 c         write (iout,*) i,icha
 228          do lll=1,chain_length
 229           cou=cou+1
 230            if (cou.le.nres) then
 231            do j=1,3
 232             kupa=mod(lll,chain_length)
 233             iprzes=(kkk-1)*chain_length+lll
 234             if (kupa.eq.0) kupa=chain_length
 235 c            write (iout,*) "kupa", kupa
 236             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 237             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 238           enddo
 239           endif
 240          enddo
 241         enddo
 242        enddo
 243        endif
 244
 245 C-koniec robienia kopidm
 246       nperm=0
 247       do kkk=1,nperm
 248       write (iout,*) "nowa struktura", nperm
 249       do i=1,nres
 250         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 251      &cref(2,i,kkk),
 252      &cref(3,i,kkk),cref(1,nres+i,kkk),
 253      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 254       enddo
 255   100 format (//'              alpha-carbon coordinates       ',
 256      &          '     centroid coordinates'/
 257      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 258      &                          10X,'X',11X,'Y',11X,'Z')
 259   110 format (a,'(',i3,')',6f12.5)
 260        enddo
 261
 262
 263       ishift_pdb=ishift
 264       return
 265       end
 266 c---------------------------------------------------------------------------
 267       subroutine int_from_cart(lside,lprn)
 268       implicit real*8 (a-h,o-z)
 269       include 'DIMENSIONS'
 270       include 'COMMON.LOCAL'
 271       include 'COMMON.VAR'
 272       include 'COMMON.CHAIN'
 273       include 'COMMON.INTERACT'
 274       include 'COMMON.IOUNITS'
 275       include 'COMMON.GEO'
 276       include 'COMMON.NAMES'
 277       character*3 seq,atom,res
 278       character*80 card
 279       dimension sccor(3,20)
 280       integer rescode
 281       logical lside,lprn
 282       if (lprn) then
 283         write (iout,'(/a)')
 284      &  'Internal coordinates calculated from crystal structure.'
 285         if (lside) then
 286           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 287      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 288      & '     Omega'
 289         else
 290           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 291      & '       Phi'
 292         endif
 293       endif
 294       call flush(iout)
 295       do i=nnt+1,nct
 296         iti=itype(i)
 297 c        write (iout,*) i,dist(i,i-1)
 298         if ((dist(i,i-1).lt.2.0D0 .or. dist(i,i-1).gt.5.0D0)
 299      &.and.(iti.ne.ntyp1).and.(itype(i-1).ne.ntyp1)) then
 300           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 301           stop
 302         endif
 303         theta(i+1)=alpha(i-1,i,i+1)
 304         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 305       enddo
 306       if (lside) then
 307         do i=2,nres-1
 308           do j=1,3
 309             c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
 310           enddo
 311           iti=itype(i)
 312           di=dist(i,nres+i)
 313           if (iti.ne.10) then
 314             alph(i)=alpha(nres+i,i,maxres2)
 315             omeg(i)=beta(nres+i,i,maxres2,i+1)
 316           endif
 317           if (lprn)
 318      &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 319      &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,rad2deg*alph(i),
 320      &    rad2deg*omeg(i)
 321         enddo
 322       else if (lprn) then
 323         do i=2,nres
 324           iti=itype(i)
 325           write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 326      &    rad2deg*theta(i),rad2deg*phi(i)
 327         enddo
 328       endif
 329       return
 330       end
 331
 332 c-------------------------------------------------------------------------------
 333       subroutine sc_loc_geom(lprn)
 334       implicit real*8 (a-h,o-z)
 335       include 'DIMENSIONS'
 336       include 'COMMON.LOCAL'
 337       include 'COMMON.VAR'
 338       include 'COMMON.CHAIN'
 339       include 'COMMON.INTERACT'
 340       include 'COMMON.IOUNITS'
 341       include 'COMMON.GEO'
 342       include 'COMMON.NAMES'
 343       include 'COMMON.CONTROL'
 344       include 'COMMON.SETUP'
 345       double precision x_prime(3),y_prime(3),z_prime(3)
 346       logical lprn
 347       do i=1,nres-1
 348         do j=1,3
 349           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 350         enddo
 351       enddo
 352       do i=2,nres-1
 353         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 354           do j=1,3
 355             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 356           enddo
 357         else
 358           do j=1,3
 359             dc_norm(j,i+nres)=0.0d0
 360           enddo
 361         endif
 362       enddo
 363       do i=2,nres-1
 364         costtab(i+1) =dcos(theta(i+1))
 365         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 366         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 367         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 368         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 369         cosfac=dsqrt(cosfac2)
 370         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 371         sinfac=dsqrt(sinfac2)
 372         it=itype(i)
 373         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 374 c
 375 C  Compute the axes of tghe local cartesian coordinates system; store in
 376 c   x_prime, y_prime and z_prime
 377 c
 378         do j=1,3
 379           x_prime(j) = 0.00
 380           y_prime(j) = 0.00
 381           z_prime(j) = 0.00
 382         enddo
 383         do j = 1,3
 384           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 385           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 386         enddo
 387         call vecpr(x_prime,y_prime,z_prime)
 388 c
 389 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 390 C to local coordinate system. Store in xx, yy, zz.
 391 c
 392         xx=0.0d0
 393         yy=0.0d0
 394         zz=0.0d0
 395         do j = 1,3
 396           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 397           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 398           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 399         enddo
 400
 401         xxref(i)=xx
 402         yyref(i)=yy
 403         zzref(i)=zz
 404         else
 405         xxref(i)=0.0d0
 406         yyref(i)=0.0d0
 407         zzref(i)=0.0d0
 408         endif
 409       enddo
 410       if (lprn) then
 411         do i=2,nres
 412           iti=itype(i)
 413           if(me.eq.king.or..not.out1file)
 414      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 415      &      yyref(i),zzref(i)
 416         enddo
 417       endif
 418       return
 419       end
 420 c---------------------------------------------------------------------------
 421       subroutine sccenter(ires,nscat,sccor)
 422       implicit real*8 (a-h,o-z)
 423       include 'DIMENSIONS'
 424       include 'COMMON.CHAIN'
 425       dimension sccor(3,20)
 426       do j=1,3
 427         sccmj=0.0D0
 428         do i=1,nscat
 429           sccmj=sccmj+sccor(j,i)
 430         enddo
 431         dc(j,ires)=sccmj/nscat
 432       enddo
 433       return
 434       end