readpdb & fort.2 from unres
[unres.git] / source / unres / src-HCD-5D / readpdb-mult.F
index 8fd1da8..9b99e64 100644 (file)
@@ -21,7 +21,7 @@ C geometry.
       integer rescode,iterter(maxres),cou
       logical fail
       integer i,j,iii,ires,ires_old,ishift,ishift1,ibeg
-      double precision dcj
+      double precision dcj,efree_temp
       bfac=0.0d0
       do i=1,maxres
          iterter(i)=0
@@ -33,103 +33,141 @@ C geometry.
       nbfrag=0
       do
         read (ipdbin,'(a80)',end=10) card
+!       write (iout,'(a)') 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)
+          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----------------------------------------
+          nbfrag=nbfrag+1
+          lsecondary=.true.
+          read(card(24:26),*) bfrag(1,nbfrag)
+          read(card(35:37),*) bfrag(2,nbfrag)
+!rc----------------------------------------
+!rc  to be corrected !!!
+          bfrag(3,nbfrag)=bfrag(1,nbfrag)
+          bfrag(4,nbfrag)=bfrag(2,nbfrag)
+!rc----------------------------------------
         endif
         if (card(:3).eq.'END') then
           goto 10
         else if (card(:3).eq.'TER') then
-C End current chain
+! End current chain
           ires_old=ires+2
-          itype(ires_old-1)=ntyp1 
+          itype(ires_old-1)=ntyp1
           iterter(ires_old-1)=1
           itype(ires_old)=ntyp1
-          iterter(ires_old)=1
+          ishift1=ishift1+1
           ibeg=2
-          write (iout,*) "Chain ended",ires,ishift,ires_old
+!          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
             do j=1,3
               dc(j,ires)=sccor(j,iii)
             enddo
-          else 
+          else
             call sccenter(ires,iii,sccor)
           endif
+c          iii=0
         endif
-C Fish out the ATOM cards.
-c        if (index(card(1:4),'ATOM').gt.0) then  
-c          read (card(14:16),'(a3)') atom
-c          if (atom.eq.'CA' .or. atom.eq.'CH3') then
-C Calculate the CM of the preceding residue.
-        read (card(23:26),*) ires
-        read (card(18:20),'(a3)') res
-        if (ires-ishift+ishift1.ne.ires_old) then
-          if (ibeg.eq.0) then
-            if (unres_pdb) then
-              do j=1,3
-                dc(j,ires+nres)=sccor(j,iii)
-              enddo
+! Read free energy
+        if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
+! Fish out the ATOM cards.
+        if (index(card(1:4),'ATOM').gt.0) then  
+          read (card(12:16),*) atom
+c          write (2,'(a)') card
+!          write (iout,*) "! ",atom," !",ires
+!          if (atom.eq.'CA' .or. atom.eq.'CH3') then
+          read (card(23:26),*) ires
+          read (card(18:20),'(a3)') res
+!          write (iout,*) "ires",ires,ires-ishift+ishift1,
+!     &      " ires_old",ires_old
+!          write (iout,*) "ishift",ishift," ishift1",ishift1
+!          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+          if (ires-ishift+ishift1.ne.ires_old) then
+! Calculate the CM of the preceding residue.
+!            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
+            if (ibeg.eq.0) then
+!              write (iout,*) "Calculating sidechain center iii",iii
+              if (unres_pdb) then
+                do j=1,3
+                  dc(j,ires+nres)=sccor(j,iii)
+                enddo
+              else
+                call sccenter(ires_old,iii,sccor)
+              endif
+              iii=0
+            endif
+! Start new residue.
+            if (res.eq.'Cl-' .or. res.eq.'Na+') then
+              ires=ires_old
+              cycle
+            else if (ibeg.eq.1) then
+c              write (iout,*) "BEG ires",ires
+              ishift=ires-1
+              if (res.ne.'GLY' .and. res.ne. 'ACE') then
+                ishift=ishift-1
+                itype(1)=ntyp1
+              endif
+              ires=ires-ishift+ishift1
+              ires_old=ires
+!              write (iout,*) "ishift",ishift," ires",ires,&
+!               " ires_old",ires_old
+              ibeg=0 
+            else if (ibeg.eq.2) then
+! Start a new chain
+              ishift=-ires_old+ires-1 !!!!!
+              ishift1=ishift1-1    !!!!!
+!              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
+              ires=ires-ishift+ishift1
+              ires_old=ires
+              ibeg=0
             else
-              call sccenter(ires,iii,sccor)
+              ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+              ires=ires-ishift+ishift1
+              ires_old=ires
             endif
-          endif
-C Start new residue.
-c            write (iout,'(a80)') card
-          if (ibeg.eq.1) then
-            ishift=ires-1
-            if (res.ne.'GLY' .and. res.ne. 'ACE') then
-              ishift=ishift-1
-              itype(1)=ntyp1
+            if (res.eq.'ACE' .or. res.eq.'NHE') then
+              itype(ires)=10
+            else
+              itype(ires)=rescode(ires,res,0)
             endif
-c             write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
-            ibeg=0          
-          else if (ibeg.eq.2) then
-c Start a new chain
-            ishift=-ires_old+ires-1
-c              write (iout,*) "New chain started",ires,ishift
-            ibeg=0
-          endif
-        else
-          ires=ires-ishift
-c            write (2,*) "ires",ires," ishift",ish
-        endif
-        if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
-     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
-          if (res.eq.'ACE') then
-            itype(ires)=10
           else
-            itype(ires)=rescode(ires,res,0)
+            ires=ires-ishift+ishift1
           endif
-          read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-          read(card(61:66),*) bfac(ires)
-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
+!          write (iout,*) "ires_old",ires_old," ires",ires
+          if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+!            ishift1=ishift1+1
+          endif
+!          write (2,*) "ires",ires," res ",res!," ity"!,ity 
+          if (atom.eq.'CA' .or. atom.eq.'CH3' .or. 
+     &       res.eq.'NHE'.and.atom(:2).eq.'HN') then
+            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+!            write (iout,*) "backbone ",atom
+#ifdef DEBUG
+            write (iout,'(2i3,2x,a,3f8.3)') 
+     &      ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
             iii=iii+1
-          read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+            do j=1,3
+              sccor(j,iii)=c(j,ires)
+            enddo
+c            write (2,*) card(23:27),ires,itype(ires),iii
+          else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. 
+     &             atom.ne.'N' .and. atom.ne.'C' .and. 
+     &             atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. 
+     &             atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+!            write (iout,*) "sidechain ",atom
+            iii=iii+1
+            read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+c            write (2,*) "iii",iii
+          endif
         endif
       enddo
    10 if(me.eq.king.or..not.out1file) 
      & write (iout,'(a,i5)') ' Nres: ',ires
+c      write (iout,*) "iii",iii
 C Calculate dummy residue coordinates inside the "chain" of a multichain
 C system
       nres=ires
@@ -191,6 +229,7 @@ C Calculate the CM of the last side chain.
           dc(j,ires)=sccor(j,iii)
         enddo
       else 
+c        write (iout,*) "Calling sccenter iii",iii
         call sccenter(ires,iii,sccor)
       endif
       nsup=nres
@@ -251,11 +290,15 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue
       endif
 C Calculate internal coordinates.
       if(me.eq.king.or..not.out1file)then
-       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
+      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)"
+      do ires=1,nres
+        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+     &    (c(j,ires+nres),j=1,3)
+      enddo
       endif
       call flush(iout)
 c      write(iout,*)"before int_from_cart nres",nres