source/unres/src_MD_DFA/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 C 10/03/12 Adam: Correction for zero SC-SC bond length
 264           if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0)
 265      &     di=dsc(itype(i))
 266           vbld(i+nres)=di
 267           if (itype(i).ne.10) then
 268             vbld_inv(i+nres)=1.0d0/di
 269           else
 270             vbld_inv(i+nres)=0.0d0
 271           endif
 272           if (iti.ne.10) then
 273             alph(i)=alpha(nres+i,i,maxres2)
 274             omeg(i)=beta(nres+i,i,maxres2,i+1)
 275           endif
 276           if(me.eq.king.or..not.out1file)then
 277            if (lprn)
 278      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
 279      &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
 280      &     rad2deg*alph(i),rad2deg*omeg(i)
 281           endif
 282         enddo
 283       else if (lprn) then
 284         do i=2,nres
 285           iti=itype(i)
 286           if(me.eq.king.or..not.out1file)
 287      &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
 288      &     rad2deg*theta(i),rad2deg*phi(i)
 289         enddo
 290       endif
 291       return
 292       end
 293 c-------------------------------------------------------------------------------
 294       subroutine sc_loc_geom(lprn)
 295       implicit real*8 (a-h,o-z)
 296       include 'DIMENSIONS'
 297 #ifdef MPI
 298       include "mpif.h"
 299 #endif
 300       include 'COMMON.LOCAL'
 301       include 'COMMON.VAR'
 302       include 'COMMON.CHAIN'
 303       include 'COMMON.INTERACT'
 304       include 'COMMON.IOUNITS'
 305       include 'COMMON.GEO'
 306       include 'COMMON.NAMES'
 307       include 'COMMON.CONTROL'
 308       include 'COMMON.SETUP'
 309       double precision x_prime(3),y_prime(3),z_prime(3)
 310       logical lprn
 311       do i=1,nres-1
 312         do j=1,3
 313           dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
 314         enddo
 315       enddo
 316       do i=2,nres-1
 317         if (itype(i).ne.10) then
 318           do j=1,3
 319             dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
 320           enddo
 321         else
 322           do j=1,3
 323             dc_norm(j,i+nres)=0.0d0
 324           enddo
 325         endif
 326       enddo
 327       do i=2,nres-1
 328         costtab(i+1) =dcos(theta(i+1))
 329         sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
 330         cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
 331         sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
 332         cosfac2=0.5d0/(1.0d0+costtab(i+1))
 333         cosfac=dsqrt(cosfac2)
 334         sinfac2=0.5d0/(1.0d0-costtab(i+1))
 335         sinfac=dsqrt(sinfac2)
 336         it=itype(i)
 337         if (it.ne.10) then
 338 c
 339 C  Compute the axes of tghe local cartesian coordinates system; store in
 340 c   x_prime, y_prime and z_prime
 341 c
 342         do j=1,3
 343           x_prime(j) = 0.00
 344           y_prime(j) = 0.00
 345           z_prime(j) = 0.00
 346         enddo
 347         do j = 1,3
 348           x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
 349           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
 350         enddo
 351         call vecpr(x_prime,y_prime,z_prime)
 352 c
 353 C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
 354 C to local coordinate system. Store in xx, yy, zz.
 355 c
 356         xx=0.0d0
 357         yy=0.0d0
 358         zz=0.0d0
 359         do j = 1,3
 360           xx = xx + x_prime(j)*dc_norm(j,i+nres)
 361           yy = yy + y_prime(j)*dc_norm(j,i+nres)
 362           zz = zz + z_prime(j)*dc_norm(j,i+nres)
 363         enddo
 364
 365         xxref(i)=xx
 366         yyref(i)=yy
 367         zzref(i)=zz
 368         else
 369         xxref(i)=0.0d0
 370         yyref(i)=0.0d0
 371         zzref(i)=0.0d0
 372         endif
 373       enddo
 374       if (lprn) then
 375         do i=2,nres
 376           iti=itype(i)
 377           if(me.eq.king.or..not.out1file)
 378      &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
 379      &      yyref(i),zzref(i)
 380         enddo
 381       endif
 382       return
 383       end
 384 c---------------------------------------------------------------------------
 385       subroutine sccenter(ires,nscat,sccor)
 386       implicit real*8 (a-h,o-z)
 387       include 'DIMENSIONS'
 388       include 'COMMON.CHAIN'
 389       dimension sccor(3,20)
 390       do j=1,3
 391         sccmj=0.0D0
 392         do i=1,nscat
 393           sccmj=sccmj+sccor(j,i)
 394         enddo
 395         dc(j,ires)=sccmj/nscat
 396       enddo
 397       return
 398       end
 399 c---------------------------------------------------------------------------
 400       subroutine bond_regular
 401       implicit real*8 (a-h,o-z)
 402       include 'DIMENSIONS'
 403       include 'COMMON.VAR'
 404       include 'COMMON.LOCAL'
 405       include 'COMMON.CALC'
 406       include 'COMMON.INTERACT'
 407       include 'COMMON.CHAIN'
 408       do i=1,nres-1
 409        vbld(i+1)=vbl
 410        vbld_inv(i+1)=1.0d0/vbld(i+1)
 411        vbld(i+1+nres)=dsc(itype(i+1))
 412        vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
 413 c       print *,vbld(i+1),vbld(i+1+nres)
 414       enddo
 415       return
 416       end
 417