Adasko's dir
[unres.git] / source / cluster / wham / src-NEWSC / readpdb.f
diff --git a/source/cluster/wham/src-NEWSC/readpdb.f b/source/cluster/wham/src-NEWSC/readpdb.f
new file mode 100644 (file)
index 0000000..de5811c
--- /dev/null
@@ -0,0 +1,183 @@
+      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'
+      character*3 seq,atom,res
+      character*80 card
+      dimension sccor(3,20)
+      integer rescode
+      ibeg=1
+      do i=1,10000
+        read (ipdbin,'(a80)',end=10) card
+        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) call sccenter(ires,iii,sccor)
+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)
+            write (iout,'(2i3,2x,a,3f8.3)') 
+     &      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 write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate the CM of the last side chain.
+      call sccenter(ires,iii,sccor)
+      nres=ires
+      nsup=nres
+      nstart_sup=1
+      if (itype(nres).ne.10) then
+        nres=nres+1
+        itype(nres)=21
+        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
+      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
+        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
+C Copy the coordinates to reference coordinates
+      do i=1,2*nres
+        do j=1,3
+          cref(j,i)=c(j,i)
+        enddo
+      enddo
+C Calculate internal coordinates.
+c     do ires=1,nres
+c       write (iout,'(2i3,2x,a,3f8.3)') 
+c    &    ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3)
+c     enddo
+      call int_from_cart(.true.,.true.)
+      return
+      end
+c---------------------------------------------------------------------------
+      subroutine int_from_cart(lside,lprn)
+      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'
+      character*3 seq,atom,res
+      character*80 card
+      dimension sccor(3,20)
+      integer rescode
+      logical lside,lprn
+      if (lprn) then 
+        write (iout,'(/a)') 
+     &  'Internal coordinates calculated from crystal structure.'
+        if (lside) then 
+          write (iout,'(8a)') '  Res  ','       dvb','     Theta',
+     & '       Phi','    Dsc_id','       Dsc','     Alpha',
+     & '     Omega'
+        else 
+          write (iout,'(4a)') '  Res  ','       dvb','     Theta',
+     & '       Phi'
+        endif
+      endif
+      call flush(iout)
+      do i=nnt+1,nct
+        iti=itype(i)
+c        write (iout,*) i,dist(i,i-1)
+        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
+          stop
+        endif
+        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
+      if (lside) then
+        do i=2,nres-1
+          do j=1,3
+            c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
+          enddo
+          iti=itype(i)
+          di=dist(i,nres+i)
+          if (iti.ne.10) then
+            alph(i)=alpha(nres+i,i,maxres2)
+            omeg(i)=beta(nres+i,i,maxres2,i+1)
+          endif
+          if (lprn)
+     &    write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
+     &    rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,rad2deg*alph(i),
+     &    rad2deg*omeg(i)
+        enddo
+      else if (lprn) then
+        do i=2,nres
+          iti=itype(i)
+          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 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