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