Fixed the addition of dummy chain ends to structures read from UNRES PDB files

[unres.git] / source / unres / src_CSA_DiL / readpdb.F

diff --git a/source/unres/src_CSA_DiL/readpdb.F b/source/unres/src_CSA_DiL/readpdb.F
index 69ec3e1..7810b24 100644 (file)
--- a/source/unres/src_CSA_DiL/readpdb.F
+++ b/source/unres/src_CSA_DiL/readpdb.F
@@ -16,6 +16,8 @@ C geometry.
       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.
@@ -104,9 +106,16 @@ C Calculate the CM of the last side chain.
         nres=nres+1
         itype(nres)=21
         if (unres_pdb) then
-          c(1,nres)=c(1,nres-1)+3.8d0
-          c(2,nres)=c(2,nres-1)
-          c(3,nres)=c(3,nres-1)
+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
         do j=1,3
           dcj=c(j,nres-2)-c(j,nres-3)
@@ -128,9 +137,16 @@ C Calculate the CM of the last side chain.
         nsup=nsup-1
         nstart_sup=2
         if (unres_pdb) then
-          c(1,1)=c(1,2)-3.8d0
-          c(2,1)=c(2,2)
-          c(3,1)=c(3,2)
+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
         do j=1,3
           dcj=c(j,4)-c(j,3)
@@ -258,8 +274,8 @@ c      endif
           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.21. and. di.eq.0.0d0)
+C 9/29/12 Adam: Correction for zero SC-SC bond length 
+          if (itype(i).ne.10 .and. itype(i).ne.21. and. di.eq.0.0d0) 
      &     di=dsc(itype(i))
           vbld(i+nres)=di
           if (itype(i).ne.10) then
@@ -410,4 +426,64 @@ c       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end
+C-----------------------------------------------------------------------------
+      subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+c this subroutine calculates unity vectors of a local reference system
+c defined by atoms (i2), (i3), and (i4). the x axis is the axis from
+c atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
+c (i2), (i3), and (i4). z axis is directed according to the sign of the
+c vector product (i3)-(i2) and (i3)-(i4). sets fail to .true. if atoms
+c (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
+c form a linear fragment. returns vectors e1, e2, and e3.
+      logical fail
+      double precision e1(3),e2(3),e3(3)
+      double precision u(3),z(3)
+      include "COMMON.CHAIN"
+      data coinc /1.0d-13/,align /1.0d-13/
+      fail=.false.
+      s1=0.0d0
+      s2=0.0d0
+      do 1 i=1,3
+      zi=c(i,i2)-c(i,i3)
+      ui=c(i,i4)-c(i,i3)
+      s1=s1+zi*zi
+      s2=s2+ui*ui
+      z(i)=zi
+    1 u(i)=ui
+      s1=sqrt(s1)
+      s2=sqrt(s2)
+      if (s1.gt.coinc) goto 2
+      write (iout,1000) i2,i3,i1
+      fail=.true.
+      return
+    2 if (s2.gt.coinc) goto 4
+      write(iout,1000) i3,i4,i1
+      fail=.true.
+      return
+    4 s1=1.0/s1
+      s2=1.0/s2
+      v1=z(2)*u(3)-z(3)*u(2)
+      v2=z(3)*u(1)-z(1)*u(3)
+      v3=z(1)*u(2)-z(2)*u(1)
+      anorm=sqrt(v1*v1+v2*v2+v3*v3)
+      if (anorm.gt.align) goto 6
+      write (iout,1010) i2,i3,i4,i1
+      fail=.true.
+      return
+    6 anorm=1.0/anorm
+      e3(1)=v1*anorm
+      e3(2)=v2*anorm
+      e3(3)=v3*anorm
+      e1(1)=z(1)*s1
+      e1(2)=z(2)*s1
+      e1(3)=z(3)*s1
+      e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
+      e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
+      e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
+ 1000 format (/1x,' * * * error - atoms',i4,' and',i4,' coincide.')
+ 1010 format (/1x,' * * * error - atoms',2(i4,2h, ),i4,' form a linear')
+      return
+      end