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