X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2Freadpdb.F;fp=source%2Funres%2Fsrc_Eshel%2Freadpdb.F;h=5d6acc0c3adf1dc719bdf941c8d0c6e51ce35f9c;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_Eshel/readpdb.F b/source/unres/src_Eshel/readpdb.F
new file mode 100644
index 0000000..5d6acc0
--- /dev/null
+++ b/source/unres/src_Eshel/readpdb.F
@@ -0,0 +1,441 @@
+      subroutine readpdb
+C Read the PDB file and convert the peptide geometry into virtual-chain 
+C geometry.
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.SETUP'
+      character*3 seq,atom,res
+      character*80 card
+      dimension sccor(3,20)
+      double precision e1(3),e2(3),e3(3)
+      logical fail
+      integer rescode
+      ibeg=1
+      lsecondary=.false.
+      nhfrag=0
+      nbfrag=0
+      do i=1,maxres
+        itype(i)=21
+        do j=1,3
+          c(j,i)=0.0d0
+          c(j,i+nres)=0.0d0
+        enddo
+      enddo
+      do i=1,10000
+        read (ipdbin,'(a80)',end=10) card
+        if (card(:5).eq.'HELIX') then
+         nhfrag=nhfrag+1
+         lsecondary=.true.
+         read(card(22:25),*) hfrag(1,nhfrag)
+         read(card(34:37),*) hfrag(2,nhfrag)
+        endif
+        if (card(:5).eq.'SHEET') then
+         nbfrag=nbfrag+1
+         lsecondary=.true.
+         read(card(24:26),*) bfrag(1,nbfrag)
+         read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc  to be corrected !!!
+         bfrag(3,nbfrag)=bfrag(1,nbfrag)
+         bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
+        endif
+        if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+C Fish out the ATOM cards.
+        if (index(card(1:4),'ATOM').gt.0) then  
+          read (card(14:16),'(a3)') atom
+          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+C Calculate the CM of the preceding residue.
+            if (ibeg.eq.0) then
+              if (unres_pdb) then
+                do j=1,3
+                  dc(j,ires+nres)=sccor(j,iii)
+                enddo
+              else
+                call sccenter(ires,iii,sccor)
+              endif
+            endif
+C Start new residue.
+            read (card(24:26),*) ires
+            read (card(18:20),'(a3)') res
+            if (ibeg.eq.1) then
+              ishift=ires-1
+              if (res.ne.'GLY' .and. res.ne. 'ACE') then
+                ishift=ishift-1
+                itype(1)=21
+              endif
+              ibeg=0          
+            endif
+            ires=ires-ishift
+            if (res.eq.'ACE') then
+              ity=10
+            else
+              itype(ires)=rescode(ires,res,0)
+            endif
+            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+c            if(me.eq.king.or..not.out1file)
+c     &       write (iout,'(2i3,2x,a,3f8.3)') 
+c     &       ires,itype(ires),res,(c(j,ires),j=1,3)
+            iii=1
+            do j=1,3
+              sccor(j,iii)=c(j,ires)
+            enddo
+          else if (atom.ne.'O  '.and.atom(1:1).ne.'H' .and. 
+     &             atom.ne.'N  ' .and. atom.ne.'C   ') then
+            iii=iii+1
+            read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+          endif
+        endif
+      enddo
+   10 if(me.eq.king.or..not.out1file) 
+     & write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate the CM of the last side chain.
+      if (unres_pdb) then
+        do j=1,3
+          dc(j,ires+nres)=sccor(j,iii)
+        enddo
+      else if (.not.catrace) then
+        call sccenter(ires,iii,sccor)
+      endif
+      nres=ires
+      nsup=nres
+      nstart_sup=1
+      if (itype(nres).ne.10) then
+        nres=nres+1
+        itype(nres)=21
+        if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+          enddo
+        else if (.not.catrace) then
+        do j=1,3
+          dcj=c(j,nres-2)-c(j,nres-3)
+          c(j,nres)=c(j,nres-1)+dcj
+          c(j,2*nres)=c(j,nres)
+        enddo
+        endif
+      endif
+      do i=2,nres-1
+        do j=1,3
+          c(j,i+nres)=dc(j,i)
+        enddo
+      enddo
+      do j=1,3
+        c(j,nres+1)=c(j,1)
+        c(j,2*nres)=c(j,nres)
+      enddo
+      if (itype(1).eq.21) then
+        nsup=nsup-1
+        nstart_sup=2
+        if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+          call refsys(2,3,4,e1,e2,e3,fail)
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+          do j=1,3
+            c(j,1)=c(j,2)-3.8d0*e2(j)
+          enddo
+        else if (.not.catrace) then
+        do j=1,3
+          dcj=c(j,4)-c(j,3)
+          c(j,1)=c(j,2)-dcj
+          c(j,nres+1)=c(j,1)
+        enddo
+        endif
+      endif
+C Calculate internal coordinates.
+      if(me.eq.king.or..not.out1file)then
+       write (iout,'(a)') 
+     &   "Backbone and SC coordinates as read from the PDB"
+       do ires=1,nres
+        write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') 
+     &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+     &    (c(j,nres+ires),j=1,3)
+       enddo
+      endif
+      call int_from_cart(.not.catrace,.false.)
+      if (.not.catrace) call sc_loc_geom(.false.)
+      do i=1,nres
+        thetaref(i)=theta(i)
+        phiref(i)=phi(i)
+      enddo
+      do i=1,nres-1
+        do j=1,3
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+        enddo
+      enddo
+      do i=2,nres-1
+        do j=1,3
+          dc(j,i+nres)=c(j,i+nres)-c(j,i)
+          dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+        enddo
+c        write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+c     &   vbld_inv(i+nres)
+      enddo
+c      call chainbuild
+C Copy the coordinates to reference coordinates
+      do i=1,2*nres
+        do j=1,3
+          cref(j,i)=c(j,i)
+        enddo
+      enddo
+
+
+      do j=1,nbfrag     
+        do i=1,4                                                       
+         bfrag(i,j)=bfrag(i,j)-ishift
+        enddo
+      enddo
+
+      do j=1,nhfrag
+        do i=1,2
+         hfrag(i,j)=hfrag(i,j)-ishift
+        enddo
+      enddo
+
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine int_from_cart(lside,lprn)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+#endif
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SETUP'
+      character*3 seq,atom,res
+      character*80 card
+      dimension sccor(3,20)
+      integer rescode
+      logical lside,lprn
+      if(me.eq.king.or..not.out1file)then
+       if (lprn) then 
+        write (iout,'(/a)') 
+     &  'Internal coordinates calculated from crystal structure.'
+        if (lside) then 
+          write (iout,'(8a)') '  Res  ','       dvb','     Theta',
+     & '     Gamma','    Dsc_id','       Dsc','     Alpha',
+     & '     Beta '
+        else 
+          write (iout,'(4a)') '  Res  ','       dvb','     Theta',
+     & '     Gamma'
+        endif
+       endif
+      endif
+      do i=1,nres-1
+        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+        iti=itype(i)
+        if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+          write (iout,'(a,i4)') 'Bad Cartesians for residue',i
+ctest          stop
+        endif
+        vbld(i+1)=dist(i,i+1)
+        vbld_inv(i+1)=1.0d0/vbld(i+1)
+        if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
+        if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+      enddo
+c      if (unres_pdb) then
+c        if (itype(1).eq.21) then
+c          theta(3)=90.0d0*deg2rad
+c          phi(4)=180.0d0*deg2rad
+c          vbld(2)=3.8d0
+c          vbld_inv(2)=1.0d0/vbld(2)
+c        endif
+c        if (itype(nres).eq.21) then
+c          theta(nres)=90.0d0*deg2rad
+c          phi(nres)=180.0d0*deg2rad
+c          vbld(nres)=3.8d0
+c          vbld_inv(nres)=1.0d0/vbld(2)
+c        endif
+c      endif
+      if (lside) then
+        do i=2,nres-1
+          if (itype(i).eq.ntyp1) cycle
+          do j=1,3
+            c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+     &     +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
+          enddo
+          iti=itype(i)
+          di=dist(i,nres+i)
+C 10/03/12 Adam: Correction for zero SC-SC bond length
+          if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
+     &     di=dsc(itype(i))
+          vbld(i+nres)=di
+          if (itype(i).ne.10) then
+            vbld_inv(i+nres)=1.0d0/di
+          else
+            vbld_inv(i+nres)=0.0d0
+          endif
+          if (iti.ne.10) then
+            alph(i)=alpha(nres+i,i,maxres2)
+            omeg(i)=beta(nres+i,i,maxres2,i+1)
+          endif
+          if(me.eq.king.or..not.out1file)then
+           if (lprn)
+     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+     &     rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+     &     rad2deg*alph(i),rad2deg*omeg(i)
+          endif
+        enddo
+      else if (lprn) then
+        do i=2,nres
+          iti=itype(i)
+          if(me.eq.king.or..not.out1file)
+     &     write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
+     &     rad2deg*theta(i),rad2deg*phi(i)
+        enddo
+      endif
+      return
+      end
+c-------------------------------------------------------------------------------
+      subroutine sc_loc_geom(lprn)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include "mpif.h"
+#endif
+      include 'COMMON.LOCAL'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.NAMES'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SETUP'
+      double precision x_prime(3),y_prime(3),z_prime(3)
+      logical lprn
+      do i=1,nres-1
+        do j=1,3
+          dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
+        enddo
+      enddo
+      do i=2,nres-1
+        if (itype(i).ne.10) then
+          do j=1,3
+            dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
+          enddo
+        else
+          do j=1,3
+            dc_norm(j,i+nres)=0.0d0
+          enddo
+        endif
+      enddo
+      do i=2,nres-1
+        costtab(i+1) =dcos(theta(i+1))
+        sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+        cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+        sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+        cosfac2=0.5d0/(1.0d0+costtab(i+1))
+        cosfac=dsqrt(cosfac2)
+        sinfac2=0.5d0/(1.0d0-costtab(i+1))
+        sinfac=dsqrt(sinfac2)
+        it=itype(i)
+        if (it.ne.10) then
+c
+C  Compute the axes of tghe local cartesian coordinates system; store in
+c   x_prime, y_prime and z_prime 
+c
+        do j=1,3
+          x_prime(j) = 0.00
+          y_prime(j) = 0.00
+          z_prime(j) = 0.00
+        enddo
+        do j = 1,3
+          x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+          y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+        enddo
+        call vecpr(x_prime,y_prime,z_prime)
+c
+C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
+C to local coordinate system. Store in xx, yy, zz.
+c
+        xx=0.0d0
+        yy=0.0d0
+        zz=0.0d0
+        do j = 1,3
+          xx = xx + x_prime(j)*dc_norm(j,i+nres)
+          yy = yy + y_prime(j)*dc_norm(j,i+nres)
+          zz = zz + z_prime(j)*dc_norm(j,i+nres)
+        enddo
+
+        xxref(i)=xx
+        yyref(i)=yy
+        zzref(i)=zz
+        else
+        xxref(i)=0.0d0
+        yyref(i)=0.0d0
+        zzref(i)=0.0d0
+        endif
+      enddo
+      if (lprn) then
+        do i=2,nres
+          iti=itype(i)
+          if(me.eq.king.or..not.out1file)
+     &     write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
+     &      yyref(i),zzref(i)
+        enddo
+      endif
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine sccenter(ires,nscat,sccor)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      dimension sccor(3,20)
+      do j=1,3
+        sccmj=0.0D0
+        do i=1,nscat
+          sccmj=sccmj+sccor(j,i) 
+        enddo
+        dc(j,ires)=sccmj/nscat
+      enddo
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine bond_regular
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'   
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'      
+      include 'COMMON.CALC'
+      include 'COMMON.INTERACT'
+      include 'COMMON.CHAIN'
+      do i=1,nres-1
+       vbld(i+1)=vbl
+       vbld_inv(i+1)=1.0d0/vbld(i+1)
+       vbld(i+1+nres)=dsc(itype(i+1))
+       vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
+c       print *,vbld(i+1),vbld(i+1+nres)
+      enddo
+      return
+      end