source/wham/src-M/readpdb.unr

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