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