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       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)=ntyp1
  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)=ntyp1
  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               itype(ires)=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.ntyp1) 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)=ntyp1
 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.ntyp1) 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,nres
 213       lll=lll+1
 214 cc      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 215       if (i.gt.1) then
 216       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
 217       chain_length=lll-1
 218       kkk=kkk+1
 219 c       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
 220       lll=1
 221       endif
 222       endif
 223         do j=1,3
 224           cref(j,i,cou)=c(j,i)
 225           cref(j,i+nres,cou)=c(j,i+nres)
 226           if (i.le.nres) then
 227           chain_rep(j,lll,kkk)=c(j,i)
 228           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
 229           endif
 230          enddo
 231       enddo
 232       write (iout,*) chain_length
 233       do j=1,3
 234       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
 235       chain_rep(j,chain_length+nres,symetr)
 236      &=chain_rep(j,chain_length+nres,1)
 237       enddo
 238 c diagnostic
 239 c       write (iout,*) "spraw lancuchy",chain_length,symetr
 240 c       do i=1,4
 241 c         do kkk=1,chain_length
 242 c           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
 243 c         enddo
 244 c        enddo
 245 c enddiagnostic
 246 C makes copy of chains
 247         write (iout,*) "symetr", symetr
 248
 249       if (symetr.gt.1) then
 250        call permut(symetr)
 251        nperm=1
 252        do i=1,symetr
 253        nperm=nperm*i
 254        enddo
 255        do i=1,nperm
 256        write(iout,*) (tabperm(i,kkk),kkk=1,4)
 257        enddo
 258        do i=1,nperm
 259         cou=0
 260         do kkk=1,symetr
 261          icha=tabperm(i,kkk)
 262 c         write (iout,*) i,icha
 263          do lll=1,chain_length
 264           cou=cou+1
 265            if (cou.le.nres) then
 266            do j=1,3
 267             kupa=mod(lll,chain_length)
 268             iprzes=(kkk-1)*chain_length+lll
 269             if (kupa.eq.0) kupa=chain_length
 270 c            write (iout,*) "kupa", kupa
 271             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
 272             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
 273           enddo
 274           endif
 275          enddo
 276         enddo
 277        enddo
 278        endif
 279 C-koniec robienia kopii
 280 c diag
 281       do kkk=1,nperm
 282       write (iout,*) "nowa struktura", nperm
 283       do i=1,nres
 284         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),
 285      &cref(2,i,kkk),
 286      &cref(3,i,kkk),cref(1,nres+i,kkk),
 287      &cref(2,nres+i,kkk),cref(3,nres+i,kkk)
 288       enddo
 289   100 format (//'              alpha-carbon coordinates       ',
 290      &          '     centroid coordinates'/
 291      1          '       ', 6X,'X',11X,'Y',11X,'Z',
 292      &                          10X,'X',11X,'Y',11X,'Z')
 293   110 format (a,'(',i3,')',6f12.5)
 294
 295       enddo
 296 cc enddiag
 297       do j=1,nbfrag
 298         do i=1,4
 299          bfrag(i,j)=bfrag(i,j)-ishift
 300         enddo
 301       enddo
 302
 303       do j=1,nhfrag
 304         do i=1,2
 305          hfrag(i,j)=hfrag(i,j)-ishift
 306         enddo
 307       enddo
 308
 309       return
 310       end
 311 c---------------------------------------------------------------------------
 312       subroutine int_from_cart(lside,lprn)
 313       implicit real*8 (a-h,o-z)
 314       include 'DIMENSIONS'
 315 #ifdef MPI
 316       include "mpif.h"
 317 #endif
 318       include 'COMMON.LOCAL'
 319       include 'COMMON.VAR'
 320       include 'COMMON.CHAIN'
 321       include 'COMMON.INTERACT'
 322       include 'COMMON.IOUNITS'
 323       include 'COMMON.GEO'
 324       include 'COMMON.NAMES'
 325       include 'COMMON.CONTROL'
 326       include 'COMMON.SETUP'
 327       character*3 seq,atom,res
 328       character*80 card
 329       dimension sccor(3,20)
 330       integer rescode
 331       logical lside,lprn
 332 #ifdef MPI
 333       if(me.eq.king.or..not.out1file)then
 334 #endif
 335        if (lprn) then
 336         write (iout,'(/a)')
 337      &  'Internal coordinates calculated from crystal structure.'
 338         if (lside) then
 339           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 340      & '       Phi','    Dsc_id','       Dsc','     Alpha',
 341      & '     Omega'
 342         else
 343           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 344      & '       Phi'
 345         endif
 346        endif
 347 #ifdef MPI
 348       endif
 349 #endif
 350       do i=1,nres-1
 351         iti=itype(i)
 352         if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
 353      &      (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
 354           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 355 ctest          stop
 356         endif
 357         vbld(i+1)=dist(i,i+1)
 358         vbld_inv(i+1)=1.0d0/vbld(i+1)
 359         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 360         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 361       enddo
 362 c      if (unres_pdb) then
 363 c        if (itype(1).eq.21) then
 364 c          theta(3)=90.0d0*deg2rad
 365 c          phi(4)=180.0d0*deg2rad
 366 c          vbld(2)=3.8d0
 367 c          vbld_inv(2)=1.0d0/vbld(2)
 368 c        endif
 369 c        if (itype(nres).eq.21) then
 370 c          theta(nres)=90.0d0*deg2rad
 371 c          phi(nres)=180.0d0*deg2rad
 372 c          vbld(nres)=3.8d0
 373 c          vbld_inv(nres)=1.0d0/vbld(2)
 374 c        endif
 375 c      endif
 376       if (lside) then
 377         do i=2,nres-1
 378           do j=1,3
 379             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 380      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 381           enddo
 382           iti=itype(i)
 383           di=dist(i,nres+i)
 384           vbld(i+nres)=di
 385           if (itype(i).ne.10) then
 386             vbld_inv(i+nres)=1.0d0/di
 387           else
 388             vbld_inv(i+nres)=0.0d0
 389           endif
 390           if (iti.ne.10) then
 391             alph(i)=alpha(nres+i,i,maxres2)
 392             omeg(i)=beta(nres+i,i,maxres2,i+1)
 393           endif
 394           if(me.eq.king.or..not.out1file)then
 395            if (lprn)
 396      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 397      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 398      &     rad2deg*alph(i),rad2deg*omeg(i)
 399           endif
 400         enddo
 401       else if (lprn) then
 402         do i=2,nres
 403           iti=itype(i)
 404           if(me.eq.king.or..not.out1file)
 405      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 406      &     rad2deg*theta(i),rad2deg*phi(i)
 407         enddo
 408       endif
 409       return
 410       end
 411 c-------------------------------------------------------------------------------
 412       subroutine sc_loc_geom(lprn)
 413       implicit real*8 (a-h,o-z)
 414       include 'DIMENSIONS'
 415 #ifdef MPI
 416       include "mpif.h"
 417 #endif
 418       include 'COMMON.LOCAL'
 419       include 'COMMON.VAR'
 420       include 'COMMON.CHAIN'
 421       include 'COMMON.INTERACT'
 422       include 'COMMON.IOUNITS'
 423       include 'COMMON.GEO'
 424       include 'COMMON.NAMES'
 425       include 'COMMON.CONTROL'
 426       include 'COMMON.SETUP'
 427       double precision x_prime(3),y_prime(3),z_prime(3)
 428       logical lprn
 429       do i=1,nres-1
 430         do j=1,3
 431           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 432         enddo
 433       enddo
 434       do i=2,nres-1
 435         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 436           do j=1,3
 437             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 438           enddo
 439         else
 440           do j=1,3
 441             dc_norm(j,i+nres)=0.0d0
 442           enddo
 443         endif
 444       enddo
 445       do i=2,nres-1
 446         costtab(i+1) =dcos(theta(i+1))
 447         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 448         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 449         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 450         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 451         cosfac=dsqrt(cosfac2)
 452         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 453         sinfac=dsqrt(sinfac2)
 454         it=itype(i)
 455         if (it.ne.10 .and. itype(i).ne.ntyp1) then
 456 c
 457 C  Compute the axes of tghe local cartesian coordinates system; store in
 458 c   x_prime, y_prime and z_prime
 459 c
 460         do j=1,3
 461           x_prime(j) = 0.00
 462           y_prime(j) = 0.00
 463           z_prime(j) = 0.00
 464         enddo
 465         do j = 1,3
 466           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 467           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 468         enddo
 469         call vecpr(x_prime,y_prime,z_prime)
 470 c
 471 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 472 C to local coordinate system. Store in xx, yy, zz.
 473 c
 474         xx=0.0d0
 475         yy=0.0d0
 476         zz=0.0d0
 477         do j = 1,3
 478           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 479           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 480           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 481         enddo
 482
 483         xxref(i)=xx
 484         yyref(i)=yy
 485         zzref(i)=zz
 486         else
 487         xxref(i)=0.0d0
 488         yyref(i)=0.0d0
 489         zzref(i)=0.0d0
 490         endif
 491       enddo
 492       if (lprn) then
 493         do i=2,nres
 494           iti=itype(i)
 495 #ifdef MPI
 496           if(me.eq.king.or..not.out1file)
 497      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 498      &      yyref(i),zzref(i)
 499 #else
 500           write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
 501      &     zzref(i)
 502 #endif
 503         enddo
 504       endif
 505       return
 506       end
 507 c---------------------------------------------------------------------------
 508       subroutine sccenter(ires,nscat,sccor)
 509       implicit real*8 (a-h,o-z)
 510       include 'DIMENSIONS'
 511       include 'COMMON.CHAIN'
 512       dimension sccor(3,20)
 513       do j=1,3
 514         sccmj=0.0D0
 515         do i=1,nscat
 516           sccmj=sccmj+sccor(j,i)
 517         enddo
 518         dc(j,ires)=sccmj/nscat
 519       enddo
 520       return
 521       end
 522 c---------------------------------------------------------------------------
 523       subroutine bond_regular
 524       implicit real*8 (a-h,o-z)
 525       include 'DIMENSIONS'
 526       include 'COMMON.VAR'
 527       include 'COMMON.LOCAL'
 528       include 'COMMON.CALC'
 529       include 'COMMON.INTERACT'
 530       include 'COMMON.CHAIN'
 531       do i=1,nres-1
 532        vbld(i+1)=vbl
 533        vbld_inv(i+1)=1.0d0/vbld(i+1)
 534        vbld(i+1+nres)=dsc(iabs(itype(i+1)))
 535        vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
 536 c       print *,vbld(i+1),vbld(i+1+nres)
 537       enddo
 538       return
 539       end
 540