source/unres/src_MD/src/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 i=1,10000
  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' .or. card(:3).eq.'TER') goto 10
  44 C Fish out the ATOM cards.
  45         if (index(card(1:4),'ATOM').gt.0) then
  46           read (card(14:16),'(a3)') atom
  47           if (atom.eq.'CA' .or. atom.eq.'CH3') then
  48 C Calculate the CM of the preceding residue.
  49             if (ibeg.eq.0) then
  50               if (unres_pdb) then
  51                 do j=1,3
  52                   dc(j,ires+nres)=sccor(j,iii)
  53                 enddo
  54               else
  55                 call sccenter(ires,iii,sccor)
  56               endif
  57             endif
  58 C Start new residue.
  59             read (card(24:26),*) ires
  60             read (card(18:20),'(a3)') res
  61             if (ibeg.eq.1) then
  62               ishift=ires-1
  63               if (res.ne.'GLY' .and. res.ne. 'ACE') then
  64                 ishift=ishift-1
  65                 itype(1)=21
  66               endif
  67               ibeg=0
  68             endif
  69             ires=ires-ishift
  70             if (res.eq.'ACE') then
  71               ity=10
  72             else
  73               itype(ires)=rescode(ires,res,0)
  74             endif
  75             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
  76 c            if(me.eq.king.or..not.out1file)
  77 c     &       write (iout,'(2i3,2x,a,3f8.3)')
  78 c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
  79             iii=1
  80             do j=1,3
  81               sccor(j,iii)=c(j,ires)
  82             enddo
  83           else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and.
  84      &             atom.ne.'N  ' .and. atom.ne.'C   ') then
  85             iii=iii+1
  86             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
  87           endif
  88         endif
  89       enddo
  90    10 if(me.eq.king.or..not.out1file)
  91      & write (iout,'(a,i5)') ' Nres: ',ires
  92 C Calculate the CM of the last side chain.
  93       if (unres_pdb) then
  94         do j=1,3
  95           dc(j,ires+nres)=sccor(j,iii)
  96         enddo
  97       else
  98         call sccenter(ires,iii,sccor)
  99       endif
 100       nres=ires
 101       nsup=nres
 102       nstart_sup=1
 103       if (itype(nres).ne.10) then
 104         nres=nres+1
 105         itype(nres)=21
 106         if (unres_pdb) then
 107           c(1,nres)=c(1,nres-1)+3.8d0
 108           c(2,nres)=c(2,nres-1)
 109           c(3,nres)=c(3,nres-1)
 110         else
 111         do j=1,3
 112           dcj=c(j,nres-2)-c(j,nres-3)
 113           c(j,nres)=c(j,nres-1)+dcj
 114           c(j,2*nres)=c(j,nres)
 115         enddo
 116         endif
 117       endif
 118       do i=2,nres-1
 119         do j=1,3
 120           c(j,i+nres)=dc(j,i)
 121         enddo
 122       enddo
 123       do j=1,3
 124         c(j,nres+1)=c(j,1)
 125         c(j,2*nres)=c(j,nres)
 126       enddo
 127       if (itype(1).eq.21) then
 128         nsup=nsup-1
 129         nstart_sup=2
 130         if (unres_pdb) then
 131           c(1,1)=c(1,2)-3.8d0
 132           c(2,1)=c(2,2)
 133           c(3,1)=c(3,2)
 134         else
 135         do j=1,3
 136           dcj=c(j,4)-c(j,3)
 137           c(j,1)=c(j,2)-dcj
 138           c(j,nres+1)=c(j,1)
 139         enddo
 140         endif
 141       endif
 142 C Calculate internal coordinates.
 143       if(me.eq.king.or..not.out1file)then
 144        write (iout,'(a)')
 145      &   "Backbone and SC coordinates as read from the PDB"
 146        do ires=1,nres
 147         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
 148      &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
 149      &    (c(j,nres+ires),j=1,3)
 150        enddo
 151       endif
 152       call int_from_cart(.true.,.false.)
 153       call sc_loc_geom(.false.)
 154       do i=1,nres
 155         thetaref(i)=theta(i)
 156         phiref(i)=phi(i)
 157       enddo
 158       do i=1,nres-1
 159         do j=1,3
 160           dc(j,i)=c(j,i+1)-c(j,i)
 161           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 162         enddo
 163       enddo
 164       do i=2,nres-1
 165         do j=1,3
 166           dc(j,i+nres)=c(j,i+nres)-c(j,i)
 167           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 168         enddo
 169 c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
 170 c     &   vbld_inv(i+nres)
 171       enddo
 172 c      call chainbuild
 173 C Copy the coordinates to reference coordinates
 174       do i=1,2*nres
 175         do j=1,3
 176           cref(j,i)=c(j,i)
 177         enddo
 178       enddo
 179
 180
 181       do j=1,nbfrag
 182         do i=1,4
 183          bfrag(i,j)=bfrag(i,j)-ishift
 184         enddo
 185       enddo
 186
 187       do j=1,nhfrag
 188         do i=1,2
 189          hfrag(i,j)=hfrag(i,j)-ishift
 190         enddo
 191       enddo
 192
 193       return
 194       end
 195 c---------------------------------------------------------------------------
 196       subroutine int_from_cart(lside,lprn)
 197       implicit real*8 (a-h,o-z)
 198       include 'DIMENSIONS'
 199 #ifdef MPI
 200       include "mpif.h"
 201 #endif
 202       include 'COMMON.LOCAL'
 203       include 'COMMON.VAR'
 204       include 'COMMON.CHAIN'
 205       include 'COMMON.INTERACT'
 206       include 'COMMON.IOUNITS'
 207       include 'COMMON.GEO'
 208       include 'COMMON.NAMES'
 209       include 'COMMON.CONTROL'
 210       include 'COMMON.SETUP'
 211       character*3 seq,atom,res
 212       character*80 card
 213       dimension sccor(3,20)
 214       integer rescode
 215       logical lside,lprn
 216       if(me.eq.king.or..not.out1file)then
 217        if (lprn) then
 218         write (iout,'(/a)')
 219      &  'Internal coordinates calculated from crystal structure.'
 220         if (lside) then
 221           write (iout,'(8a)') '  Res  ','       dvb','     Theta',
 222      & '     Gamma','    Dsc_id','       Dsc','     Alpha',
 223      & '     Beta '
 224         else
 225           write (iout,'(4a)') '  Res  ','       dvb','     Theta',
 226      & '     Gamma'
 227         endif
 228        endif
 229       endif
 230       do i=1,nres-1
 231         iti=itype(i)
 232         if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
 233           write (iout,'(a,i4)') 'Bad Cartesians for residue',i
 234 ctest          stop
 235         endif
 236         vbld(i+1)=dist(i,i+1)
 237         vbld_inv(i+1)=1.0d0/vbld(i+1)
 238         if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
 239         if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
 240       enddo
 241 c      if (unres_pdb) then
 242 c        if (itype(1).eq.21) then
 243 c          theta(3)=90.0d0*deg2rad
 244 c          phi(4)=180.0d0*deg2rad
 245 c          vbld(2)=3.8d0
 246 c          vbld_inv(2)=1.0d0/vbld(2)
 247 c        endif
 248 c        if (itype(nres).eq.21) then
 249 c          theta(nres)=90.0d0*deg2rad
 250 c          phi(nres)=180.0d0*deg2rad
 251 c          vbld(nres)=3.8d0
 252 c          vbld_inv(nres)=1.0d0/vbld(2)
 253 c        endif
 254 c      endif
 255       if (lside) then
 256         do i=2,nres-1
 257           do j=1,3
 258             c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
 259      &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
 260           enddo
 261           iti=itype(i)
 262           di=dist(i,nres+i)
 263           vbld(i+nres)=di
 264           if (itype(i).ne.10) then
 265             vbld_inv(i+nres)=1.0d0/di
 266           else
 267             vbld_inv(i+nres)=0.0d0
 268           endif
 269           if (iti.ne.10) then
 270             alph(i)=alpha(nres+i,i,maxres2)
 271             omeg(i)=beta(nres+i,i,maxres2,i+1)
 272           endif
 273           if(me.eq.king.or..not.out1file)then
 274            if (lprn)
 275      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 276      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 277      &     rad2deg*alph(i),rad2deg*omeg(i)
 278           endif
 279         enddo
 280       else if (lprn) then
 281         do i=2,nres
 282           iti=itype(i)
 283           if(me.eq.king.or..not.out1file)
 284      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 285      &     rad2deg*theta(i),rad2deg*phi(i)
 286         enddo
 287       endif
 288       return
 289       end
 290 c-------------------------------------------------------------------------------
 291       subroutine sc_loc_geom(lprn)
 292       implicit real*8 (a-h,o-z)
 293       include 'DIMENSIONS'
 294 #ifdef MPI
 295       include "mpif.h"
 296 #endif
 297       include 'COMMON.LOCAL'
 298       include 'COMMON.VAR'
 299       include 'COMMON.CHAIN'
 300       include 'COMMON.INTERACT'
 301       include 'COMMON.IOUNITS'
 302       include 'COMMON.GEO'
 303       include 'COMMON.NAMES'
 304       include 'COMMON.CONTROL'
 305       include 'COMMON.SETUP'
 306       double precision x_prime(3),y_prime(3),z_prime(3)
 307       logical lprn
 308       do i=1,nres-1
 309         do j=1,3
 310           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 311         enddo
 312       enddo
 313       do i=2,nres-1
 314         if (itype(i).ne.10) then
 315           do j=1,3
 316             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 317           enddo
 318         else
 319           do j=1,3
 320             dc_norm(j,i+nres)=0.0d0
 321           enddo
 322         endif
 323       enddo
 324       do i=2,nres-1
 325         costtab(i+1) =dcos(theta(i+1))
 326         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 327         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 328         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 329         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 330         cosfac=dsqrt(cosfac2)
 331         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 332         sinfac=dsqrt(sinfac2)
 333         it=itype(i)
 334         if (it.ne.10) then
 335 c
 336 C  Compute the axes of tghe local cartesian coordinates system; store in
 337 c   x_prime, y_prime and z_prime
 338 c
 339         do j=1,3
 340           x_prime(j) = 0.00
 341           y_prime(j) = 0.00
 342           z_prime(j) = 0.00
 343         enddo
 344         do j = 1,3
 345           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 346           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 347         enddo
 348         call vecpr(x_prime,y_prime,z_prime)
 349 c
 350 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 351 C to local coordinate system. Store in xx, yy, zz.
 352 c
 353         xx=0.0d0
 354         yy=0.0d0
 355         zz=0.0d0
 356         do j = 1,3
 357           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 358           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 359           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 360         enddo
 361
 362         xxref(i)=xx
 363         yyref(i)=yy
 364         zzref(i)=zz
 365         else
 366         xxref(i)=0.0d0
 367         yyref(i)=0.0d0
 368         zzref(i)=0.0d0
 369         endif
 370       enddo
 371       if (lprn) then
 372         do i=2,nres
 373           iti=itype(i)
 374           if(me.eq.king.or..not.out1file)
 375      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 376      &      yyref(i),zzref(i)
 377         enddo
 378       endif
 379       return
 380       end
 381 c---------------------------------------------------------------------------
 382       subroutine sccenter(ires,nscat,sccor)
 383       implicit real*8 (a-h,o-z)
 384       include 'DIMENSIONS'
 385       include 'COMMON.CHAIN'
 386       dimension sccor(3,20)
 387       do j=1,3
 388         sccmj=0.0D0
 389         do i=1,nscat
 390           sccmj=sccmj+sccor(j,i)
 391         enddo
 392         dc(j,ires)=sccmj/nscat
 393       enddo
 394       return
 395       end
 396 c---------------------------------------------------------------------------
 397       subroutine bond_regular
 398       implicit real*8 (a-h,o-z)
 399       include 'DIMENSIONS'
 400       include 'COMMON.VAR'
 401       include 'COMMON.LOCAL'
 402       include 'COMMON.CALC'
 403       include 'COMMON.INTERACT'
 404       include 'COMMON.CHAIN'
 405       do i=1,nres-1
 406        vbld(i+1)=vbl
 407        vbld_inv(i+1)=1.0d0/vbld(i+1)
 408        vbld(i+1+nres)=dsc(itype(i+1))
 409        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 410 c       print *,vbld(i+1),vbld(i+1+nres)
 411       enddo
 412       return
 413       end
 414