Adam's corrections
[unres.git] / source / unres / src-HCD-5D / readpdb.F
index 68db17c..1190bb1 100644 (file)
@@ -1,7 +1,7 @@
       subroutine readpdb
 C Read the PDB file and convert the peptide geometry into virtual-chain 
 C geometry.
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
@@ -11,15 +11,17 @@ C geometry.
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.DISTFIT'
+      include 'COMMON.FRAG'
       include 'COMMON.SETUP'
       include 'COMMON.SBRIDGE'
       character*3 seq,atom,res
       character*80 card
-      dimension sccor(3,50)
+      double precision sccor(3,50)
       double precision e1(3),e2(3),e3(3)
       integer rescode,iterter(maxres),cou
       logical fail
+      integer i,j,iii,ires,ires_old,ishift,ibeg
+      double precision dcj
       bfac=0.0d0
       do i=1,maxres
          iterter(i)=0
@@ -303,273 +305,10 @@ cc enddiag
       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,50)
-      integer rescode
-      logical lside,lprn
-#ifdef MPI
-      if(me.eq.king.or..not.out1file)then
-#endif
-       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
-#ifdef MPI
-      endif
-#endif
-      do i=1,nres-1
-        iti=itype(i)
-        if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and. 
-     &      (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)
-c        write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
-c     &      vbld_inv(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
-c      print *,"A TU2"
-      if (lside) then
-        do i=2,nres-1
-          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)
-          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
-        if (lprn) then
-           i=nres
-           iti=itype(nres)
-           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
-      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
-c        write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
-c     &    " vbld",vbld_inv(i+1)
-      enddo
-      do i=2,nres-1
-        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
-          do j=1,3
-            dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
-          enddo
-c          write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
-c     &    " vbld",vbld_inv(i+nres)
-        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)
-c        write (iout,*) "i",i," costab",costtab(i+1),
-c     &    " sintab",sinttab(i+1)
-c        write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
-c        write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
-        if (it.ne.10 .and. itype(i).ne.ntyp1) 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
-c        write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c        write (iout,*) "y_prime",(y_prime(j),j=1,3)
-        call vecpr(x_prime,y_prime,z_prime)
-c        write (iout,*) "z_prime",(z_prime(j),j=1,3)
-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
-#ifdef MPI
-        if (me.eq.king.or..not.out1file) then
-#endif
-        write (iout,*) "xxref,yyref,zzref"
-        do i=2,nres
-          write (iout,'(a3,i4,3f10.5)') 
-     &     restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
-        enddo
-#ifdef MPI
-        endif
-#endif
-      endif
-      return
-      end
-c---------------------------------------------------------------------------
-      subroutine sccenter(ires,nscat,sccor)
-      implicit real*8 (a-h,o-z)
-      include 'DIMENSIONS'
-      include 'COMMON.CHAIN'
-      dimension sccor(3,50)
-      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)=vblinv
-       vbld(i+1+nres)=dsc(iabs(itype(i+1)))
-       vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
-c       print *,vbld(i+1),vbld(i+1+nres)
-      enddo
-c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
-      do i=1,nchain
-        i1=chain_border(1,i)
-        i2=chain_border(2,i)
-        if (i1.gt.1) then
-          vbld(i1)=vbld(i1)/2
-          vbld_inv(i1)=vbld_inv(i1)*2
-        endif
-        if (i2.lt.nres) then
-          vbld(i2+1)=vbld(i2+1)/2
-          vbld_inv(i2+1)=vbld_inv(i2+1)*2
-        endif
-      enddo
-      return
-      end
-c---------------------------------------------------------------------------
       subroutine readpdb_template(k)
 C Read the PDB file for read_constr_homology with read2sigma
 C and convert the peptide geometry into virtual-chain geometry.
-      implicit real*8 (a-h,o-z)
+      implicit none
       include 'DIMENSIONS'
       include 'COMMON.LOCAL'
       include 'COMMON.VAR'
@@ -579,10 +318,10 @@ C and convert the peptide geometry into virtual-chain geometry.
       include 'COMMON.GEO'
       include 'COMMON.NAMES'
       include 'COMMON.CONTROL'
-      include 'COMMON.DISTFIT'
+      include 'COMMON.FRAG'
       include 'COMMON.SETUP'
-      integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
-     &  ishift_pdb
+      integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+     &  ishift_pdb,ires_ca
       logical lprn /.false./,fail
       double precision e1(3),e2(3),e3(3)
       double precision dcj,efree_temp
@@ -846,9 +585,9 @@ C Calculate internal coordinates.
       write (iout,'(/a)') 
      &  "Cartesian coordinates of the reference structure"
       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') 
-     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+     & "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
       do ires=1,nres
-        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') 
+        write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)') 
      &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
      &    (c(j,ires+nres),j=1,3)
       enddo