split matrices working reading pdb

diff --git a/source/unres/control.F90 b/source/unres/control.F90
index 87a0793..9b5f441 100644 (file)
--- a/source/unres/control.F90
+++ b/source/unres/control.F90
@@ -490,7 +490,7 @@
             ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,&
             ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,&
             ichunk,int_index_old
-
+      integer,dimension(5) :: nct_molec
 !el      allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
 !el      allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
 
@@ -632,7 +632,11 @@
       ispp=4 !?? wham ispp=2
 #ifdef MPI
 ! Now partition the electrostatic-interaction array
-      npept=nct-nnt
+      if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then
+      npept=nres_molec(1)-nnt-1
+      else
+      npept=nres_molec(1)-nnt
+      endif
       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
       if (lprint) &
@@ -644,7 +648,13 @@
       iatel_e=0
       ind_eleint=0
       ind_eleint_old=0
-      do i=nnt,nct-3
+      if (itype(nres_molec(1),1).eq.ntyp_molec(1)) then
+      nct_molec(1)=nres_molec(1)-1
+      else
+      nct_molec(1)=nres_molec(1)
+      endif
+
+      do i=nnt,nct_molec(1)-3
         ijunk=0
         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,&
           iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
@@ -660,7 +670,7 @@
       ind_eleint_vdw_old=0
       iatel_s_vdw=0
       iatel_e_vdw=0
-      do i=nnt,nct-3
+      do i=nnt,nct_molec(1)-3
         ijunk=0
         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,&
           my_ele_inde_vdw,i,&
@@ -673,16 +683,16 @@
    15 continue
 #else
       iatel_s=nnt
-      iatel_e=nct-5 ! ?? wham iatel_e=nct-3
+      iatel_e=nct_molec(1)-5 ! ?? wham iatel_e=nct-3
       do i=iatel_s,iatel_e
         ielstart(i)=i+4 ! ?? wham +2
-        ielend(i)=nct-1
+        ielend(i)=nct_molec(1)-1
       enddo
       iatel_s_vdw=nnt
-      iatel_e_vdw=nct-3
+      iatel_e_vdw=nct_molec(1)-3
       do i=iatel_s_vdw,iatel_e_vdw
         ielstart_vdw(i)=i+2
-        ielend_vdw(i)=nct-1
+        ielend_vdw(i)=nct_molec(1)-1
       enddo
 #endif
       if (lprint) then
@@ -707,7 +717,7 @@
       iatscp_e=0
       ind_scpint=0
       ind_scpint_old=0
-      do i=nnt,nct-1
+      do i=nnt,nct_molec(1)-1
         if (i.lt.nnt+iscp) then
 !d        write (iout,*) 'i.le.nnt+iscp'
           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,&
@@ -731,12 +741,12 @@
    14 continue
 #else
       iatscp_s=nnt
-      iatscp_e=nct-1
-      do i=nnt,nct-1
+      iatscp_e=nct_molec(1)-1
+      do i=nnt,nct_molec(1)-1
         if (i.lt.nnt+iscp) then
           nscp_gr(i)=1
           iscpstart(i,1)=i+iscp
-          iscpend(i,1)=nct
+          iscpend(i,1)=nct_molec(1)
         elseif (i.gt.nct-iscp) then
           nscp_gr(i)=1
           iscpstart(i,1)=nnt
@@ -746,7 +756,7 @@
           iscpstart(i,1)=nnt
           iscpend(i,1)=i-iscp
           iscpstart(i,2)=i+iscp
-          iscpend(i,2)=nct
+          iscpend(i,2)=nct_molec(1)
         endif 
       enddo ! i
 #endif
@@ -760,47 +770,48 @@
       endif ! lprint
 ! Partition local interactions
 #ifdef MPI
-      call int_bounds(nres-2,loc_start,loc_end)
+      call int_bounds(nres_molec(1)-2,loc_start,loc_end)
       loc_start=loc_start+1
       loc_end=loc_end+1
-      call int_bounds(nres-2,ithet_start,ithet_end)
+      call int_bounds(nres_molec(1)-2,ithet_start,ithet_end)
       ithet_start=ithet_start+2
       ithet_end=ithet_end+2
-      call int_bounds(nct-nnt-2,iturn3_start,iturn3_end) 
+      call int_bounds(nct_molec(1)-nnt-2,iturn3_start,iturn3_end) 
       iturn3_start=iturn3_start+nnt
       iphi_start=iturn3_start+2
       iturn3_end=iturn3_end+nnt
       iphi_end=iturn3_end+2
       iturn3_start=iturn3_start-1
       iturn3_end=iturn3_end-1
-      call int_bounds(nres-3,itau_start,itau_end)
+      call int_bounds(nres_molec(1)-3,itau_start,itau_end)
       itau_start=itau_start+3
       itau_end=itau_end+3
-      call int_bounds(nres-3,iphi1_start,iphi1_end)
+      call int_bounds(nres_molec(1)-3,iphi1_start,iphi1_end)
       iphi1_start=iphi1_start+3
       iphi1_end=iphi1_end+3
-      call int_bounds(nct-nnt-3,iturn4_start,iturn4_end) 
+      call int_bounds(nct_molec(1)-nnt-3,iturn4_start,iturn4_end) 
       iturn4_start=iturn4_start+nnt
       iphid_start=iturn4_start+2
       iturn4_end=iturn4_end+nnt
       iphid_end=iturn4_end+2
       iturn4_start=iturn4_start-1
       iturn4_end=iturn4_end-1
+      print *,"TUTUTU",nres_molec(1),nres
       call int_bounds(nres-2,ibond_start,ibond_end) 
       ibond_start=ibond_start+1
       ibond_end=ibond_end+1
-      call int_bounds(nct-nnt,ibondp_start,ibondp_end) 
+      call int_bounds(nct_molec(1)-nnt,ibondp_start,ibondp_end) 
       ibondp_start=ibondp_start+nnt
       ibondp_end=ibondp_end+nnt
-      call int_bounds1(nres-1,ivec_start,ivec_end) 
+      call int_bounds1(nres_molec(1)-1,ivec_start,ivec_end) 
 !      print *,"Processor",myrank,fg_rank,fg_rank1,
 !     &  " ivec_start",ivec_start," ivec_end",ivec_end
       iset_start=loc_start+2
       iset_end=loc_end+2
-      call int_bounds(nres,ilip_start,ilip_end)
+      call int_bounds(nres_molec(1),ilip_start,ilip_end)
       ilip_start=ilip_start
       ilip_end=ilip_end
-      call int_bounds(nres-1,itube_start,itube_end)
+      call int_bounds(nres_molec(1)-1,itube_start,itube_end)
       itube_start=itube_start
       itube_end=itube_end
       if (ndih_constr.eq.0) then
@@ -1309,17 +1320,17 @@
       endif
 #else
       loc_start=2
-      loc_end=nres-1
+      loc_end=nres_molec(1)-1
       ithet_start=3 
-      ithet_end=nres
+      ithet_end=nres_molec(1)
       iturn3_start=nnt
-      iturn3_end=nct-3
+      iturn3_end=nct_molec(1)-3
       iturn4_start=nnt
-      iturn4_end=nct-4
+      iturn4_end=nct_molec(1)-4
       iphi_start=nnt+3
-      iphi_end=nct
+      iphi_end=nct_molec(1)
       iphi1_start=4
-      iphi1_end=nres
+      iphi1_end=nres_molec(1)
       idihconstr_start=1
       idihconstr_end=ndih_constr
       ithetaconstr_start=1
@@ -1327,21 +1338,21 @@
       iphid_start=iphi_start
       iphid_end=iphi_end-1
       itau_start=4
-      itau_end=nres
+      itau_end=nres_molec(1)
       ibond_start=2
-      ibond_end=nres-1
+      ibond_end=nres_molec(1)-1
       ibondp_start=nnt
-      ibondp_end=nct-1
+      ibondp_end=nct_molec(1)-1
       ivec_start=1
-      ivec_end=nres-1
+      ivec_end=nres_molec(1)-1
       iset_start=3
-      iset_end=nres+1
+      iset_end=nres_molec(1)+1
       iint_start=2
-      iint_end=nres-1
+      iint_end=nres_molec(1)-1
       ilip_start=1
-      ilip_end=nres
+      ilip_end=nres_molec(1)
       itube_start=1
-      itube_end=nres
+      itube_end=nres_molec(1)
 #endif
 !el       common /przechowalnia/
 !      deallocate(iturn3_start_all)
@@ -1734,7 +1745,8 @@
       endif
       else if (molecule.eq.2) then
       do i=1,ntyp1_molec(molecule)
-        if (nam(1:1).eq.restyp(i,molecule)) then
+         print *,nam(1:1),restyp(i,molecule)(1:1) 
+        if (nam(1:1).eq.restyp(i,molecule)(1:1)) then
           rescode=i
           return
         endif

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index f7020ff..c21d21d 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -5184,6 +5184,9 @@
               usumsqder=usumsqder+ud(j)*uprod2   
             enddo
             estr=estr+uprod/usum
+            if (energy_dec) write (iout,*) &
+            "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,&
+            AKSC(1,iti),AKSC(1,iti)*diff*diff
             do j=1,3
              gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
             enddo

diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90
index 038dede..7b6c961 100644 (file)
--- a/source/unres/geometry.f90
+++ b/source/unres/geometry.f90
@@ -1,4 +1,4 @@
-      module geometry
+             module geometry
 !-----------------------------------------------------------------------------
       use io_units
       use names
@@ -3039,6 +3039,7 @@
       integer :: i,j,ires,nscat
       real(kind=8),dimension(3,20) :: sccor
       real(kind=8) :: sccmj
+        print *,"I am in sccenter",ires,nscat
       do j=1,3
         sccmj=0.0D0
         do i=1,nscat

diff --git a/source/unres/io.f90 b/source/unres/io.f90
index 8505d11..24549e8 100644 (file)
--- a/source/unres/io.f90
+++ b/source/unres/io.f90
@@ -742,6 +742,7 @@
       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
       allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
       allocate(itype(maxres,5)) !(maxres)
+      allocate(istype(maxres))
 !
 ! Zero out tables.
 !
@@ -996,7 +997,7 @@
 
         do i=1,nres_molec(molec)
           istype(i)=sugarcode(sequence(i,molec)(1:1),i)
-          itype(i,1)=rescode(i,sequence(i,molec)(2:4),iscode,molec)
+          itype(i,1)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
         enddo
        endif
 

diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index 277b6ba..cf23282 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -2314,7 +2314,7 @@
       use control_data
       use compare_data
       use MPI_data
-      use control, only: rescode
+      use control, only: rescode,sugarcode
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
@@ -2327,7 +2327,7 @@
 !      include 'COMMON.CONTROL'
 !      include 'COMMON.DISTFIT'
 !      include 'COMMON.SETUP'
-      integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
+      integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
 !        ishift_pdb
       logical :: lprn=.true.,fail
       real(kind=8),dimension(3) :: e1,e2,e3
@@ -2336,17 +2336,23 @@
       character(len=5) :: atom
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor
-      integer :: kkk,lll,icha,kupa     !rescode,
+      integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin      !rescode,
+      integer :: isugar
+      character*1 :: sugar
       real(kind=8) :: cou
+      real(kind=8),dimension(3) :: ccc
 !el local varables
       integer,dimension(2,maxres/3) :: hfrag_alloc
       integer,dimension(4,maxres/3) :: bfrag_alloc
       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
-
+      real(kind=8),dimension(:,:), allocatable  :: c_temporary
+      integer,dimension(:,:) , allocatable  :: itype_temporary
       efree_temp=0.0d0
       ibeg=1
       ishift1=0
       ishift=0
+      molecule=1
+      counter=0
 !      write (2,*) "UNRES_PDB",unres_pdb
       ires=0
       ires_old=0
@@ -2384,9 +2390,11 @@
         else if (card(:3).eq.'TER') then
 ! End current chain
           ires_old=ires+2
+          ishift=ishift+1
           ishift1=ishift1+1
-          itype(ires_old,1)=ntyp1
-          itype(ires_old-1,1)=ntyp1
+          itype(ires_old,molecule)=ntyp1_molec(molecule)
+          itype(ires_old-1,molecule)=ntyp1_molec(molecule)
+          nres_molec(molecule)=nres_molec(molecule)+2
           ibeg=2
 !          write (iout,*) "Chain ended",ires,ishift,ires_old
           if (unres_pdb) then
@@ -2402,6 +2410,7 @@
 ! Read free energy
         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 ! Fish out the ATOM cards.
+!        write(iout,*) 'card',card(1:20)
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
 !          write (iout,*) "! ",atom," !",ires
@@ -2436,6 +2445,7 @@
               if (res.ne.'GLY' .and. res.ne. 'ACE') then
                 ishift=ishift-1
                 itype(1,1)=ntyp1
+                nres_molec(molecule)=nres_molec(molecule)+1
               endif
               ires=ires-ishift+ishift1
               ires_old=ires
@@ -2448,17 +2458,27 @@
               ishift1=ishift1-1    !!!!!
 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
               ires=ires-ishift+ishift1
+              print *,ires,ishift,ishift1
               ires_old=ires
               ibeg=0
             else
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
               ires=ires-ishift+ishift1
               ires_old=ires
-            endif
+            endif 
+            print *,'atom',ires,atom
             if (res.eq.'ACE' .or. res.eq.'NHE') then
               itype(ires,1)=10
             else
-              itype(ires,1)=rescode(ires,res,0,1)
+             if (atom.eq.'CA  '.or.atom.eq.'N   ') then
+             molecule=1
+              itype(ires,molecule)=rescode(ires,res,0,molecule)
+!              nres_molec(molecule)=nres_molec(molecule)+1
+            else
+             molecule=2
+              itype(ires,molecule)=rescode(ires,res(3:4),0,molecule)
+!              nres_molec(molecule)=nres_molec(molecule)+1
+            endif
             endif
           else
             ires=ires-ishift+ishift1
@@ -2471,23 +2491,69 @@
           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)
+!              print *,ires,ishift,ishift1
 !            write (iout,*) "backbone ",atom
 #ifdef DEBUG
             write (iout,'(2i3,2x,a,3f8.3)') &
             ires,itype(ires,1),res,(c(j,ires),j=1,3)
 #endif
             iii=iii+1
+              nres_molec(molecule)=nres_molec(molecule)+1
             do j=1,3
               sccor(j,iii)=c(j,ires)
             enddo
+          else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
+               atom.eq."C2'" .or. atom.eq."C3'" &
+               .or. atom.eq."C4'" .or. atom.eq."O4'")) then
+            read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
+!c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
+              print *,ires,ishift,ishift1
+            counter=counter+1
+!            iii=iii+1
+!            do j=1,3
+!              sccor(j,iii)=c(j,ires)
+!            enddo
+            do j=1,3
+              c(j,ires)=c(j,ires)+ccc(j)/5.0
+            enddo
+             if (counter.eq.5) then
+!            iii=iii+1
+              nres_molec(molecule)=nres_molec(molecule)+1
+!            do j=1,3
+!              sccor(j,iii)=c(j,ires)
+!            enddo
+             counter=0
+           endif
+            print *, "ATOM",atom(1:3)
+          else if (atom(1:3).eq."C5'") then
+             read (card(19:19),'(a1)') sugar
+             isugar=sugarcode(sugar,ires)
+            if (ibeg.eq.1) then
+              istype(1)=isugar
+            else
+              istype(ires)=isugar
+            endif
+            if (unres_pdb) then
+              read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+            else
+              iii=iii+1
+              read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+            endif
 !            write (*,*) card(23:27),ires,itype(ires,1)
           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
+                   atom.ne.'OXT' .and. atom(:2).ne.'3H' &
+                   .and. atom.ne.'P  '.and. &
+                  atom(1:1).ne.'H' .and. &
+                  atom.ne.'OP1' .and. atom.ne.'OP2 ') then
 !            write (iout,*) "sidechain ",atom
+!            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
+                 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
+!                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
             iii=iii+1
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+              endif
           endif
         endif
       enddo
@@ -2496,6 +2562,11 @@
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires
+      if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
+!      print *,'I have', nres_molec(:)
+      
+      do k=1,5 
+       if (nres_molec(k).eq.0) cycle
       do i=2,nres-1
 !        write (iout,*) i,itype(i,1)
 !        if (itype(i,1).eq.ntyp1) then
@@ -2506,8 +2577,21 @@
 !            dc(j,i)=c(j,i)
 !          enddo
 !        endif
-        if (itype(i,1).eq.ntyp1) then
-         if (itype(i+1,1).eq.ntyp1) then
+        if (itype(i,k).eq.ntyp1_molec(k)) then
+         if (itype(i+1,k).eq.ntyp1_molec(k)) then
+          if (itype(i+2,k).eq.0) then 
+           print *,"masz sieczke"
+           do j=1,5
+            if (itype(i+2,j).ne.0) then
+            itype(i+1,k)=0
+            itype(i+1,j)=ntyp1_molec(j)
+            nres_molec(k)=nres_molec(k)-1
+            nres_molec(j)=nres_molec(j)+1
+            go to 3331
+            endif
+           enddo
+ 3331      continue
+          endif
 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
@@ -2556,6 +2640,7 @@
         endif  !itype.eq.ntyp1
 
       enddo
+      enddo
 ! Calculate the CM of the last side chain.
       if (iii.gt.0)  then
       if (unres_pdb) then
@@ -2569,9 +2654,11 @@
 !      nres=ires
       nsup=nres
       nstart_sup=1
+      print *,"molecule",molecule
       if (itype(nres,1).ne.10) then
         nres=nres+1
-        itype(nres,1)=ntyp1
+        itype(nres,molecule)=ntyp1_molec(molecule)
+        nres_molec(molecule)=nres_molec(molecule)+1
         if (unres_pdb) then
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
@@ -2591,6 +2678,8 @@
         enddo
         endif
       endif
+     print *,'I have',nres, nres_molec(:)
+
 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
 #ifdef WHAM_RUN
       if (nres.ne.nres0) then
@@ -2638,6 +2727,7 @@
         c(j,nres+1)=c(j,1)
         c(j,2*nres)=c(j,nres)
       enddo
+      
       if (itype(1,1).eq.ntyp1) then
         nsup=nsup-1
         nstart_sup=2
@@ -2660,6 +2750,24 @@
         enddo
         endif
       endif
+! First lets assign correct dummy to correct type of chain
+! 1) First residue
+      if (itype(1,1).eq.ntyp1) then
+        if (itype(2,1).eq.0) then
+         do j=2,5
+           if (itype(2,j).ne.0) then
+           itype(1,1)=0
+           itype(1,j)=ntyp1_molec(j)
+           nres_molec(1)=nres_molec(1)-1
+           nres_molec(j)=nres_molec(j)+1
+           go to 3231
+           endif
+         enddo
+3231    continue
+        endif
+       endif
+       print *,'I have',nres, nres_molec(:)
+
 ! Copy the coordinates to reference coordinates
 !      do i=1,2*nres
 !        do j=1,3
@@ -2673,8 +2781,8 @@
       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,1),1),ires,(c(j,ires),j=1,3),&
+        write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
+          (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
           (c(j,ires+nres),j=1,3)
       enddo
       endif
@@ -2689,7 +2797,47 @@
           (c(j,nres+ires),j=1,3)
        enddo
       endif
+! NOW LETS ROCK! SORTING
+      allocate(c_temporary(3,2*nres))
+      allocate(itype_temporary(nres,5))
+       itype_temporary(:,:)=0
+      seqalingbegin=1
+      do k=1,5
+        do i=1,nres
+         if (itype(i,k).ne.0) then
+          do j=1,3
+          c_temporary(j,seqalingbegin)=c(j,i)
+          c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
+
+          enddo
+          itype_temporary(seqalingbegin,k)=itype(i,k)
+          seqalingbegin=seqalingbegin+1
+         endif
+        enddo
+       enddo
+       do i=1,2*nres
+        do j=1,3
+        c(j,i)=c_temporary(j,i)
+        enddo
+       enddo
+       do k=1,5
+        do i=1,nres
+         itype(i,k)=itype_temporary(i,k)
+        enddo
+       enddo
+      if (lprn) then
+      write (iout,'(/a)') &
+        "Cartesian coordinates of the reference structure after sorting"
+      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,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
+          (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
+          (c(j,ires+nres),j=1,3)
+      enddo
+      endif
 
+       print *,seqalingbegin,nres
       if(.not.allocated(vbld)) then
        allocate(vbld(2*nres))
        do i=1,2*nres