NEWCORR5D working with 15k, work and iwork in random_vel might need testing
[unres4.git] / source / unres / geometry.F90
index 1d15584..866f5e2 100644 (file)
@@ -83,7 +83,7 @@
       nres2=2*nres
 ! Set lprn=.true. for debugging
       lprn = .false.
-      print *,"I ENTER CHAINBUILD"
+!      print *,"I ENTER CHAINBUILD"
 !
 ! Define the origin and orientation of the coordinate system and locate the
 ! first three CA's and SC(2).
         write (iout,1212) restyp(itype(i,1),1),i,dist(i-1,i),&
         alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
       enddo   
- 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
+ 1212 format (a3,'(',i6,')',2(f10.5,2f10.2))
 
       endif
 
 !        print *,i,vbld(i),"vbld(i)"
         vbld_inv(i)=1.0d0/vbld(i)
         vbld(nres+i)=dist(nres+i,i)
-        if (itype(i,1).ne.10) then
+        if ((itype(i,1).ne.10).and.(molnum(i).lt.4)) then
           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
         else
           vbld_inv(nres+i)=0.0d0
        rad2deg*alph(i),rad2deg*omeg(i)
       enddo
       endif
- 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
+ 1212 format (a3,'(',i6,')',2(f15.10,2f10.2))
 #ifdef TIMING
       time_intfcart=time_intfcart+MPI_Wtime()-time01
 #endif
       dV=2.0D0*5.0D0*deg2rad*deg2rad
       print *,'dv=',dv
       do 10 it=1,1 
-        if (it.eq.10) goto 10 
+        if ((it.eq.10).or.(it.eq.ntyp1)) goto 10 
         open (20,file=restyp(it,1)//'_distr.sdc',status='unknown')
-        call gen_side(it,90.0D0 * deg2rad,al,om,fail)
+        call gen_side(it,90.0D0 * deg2rad,al,om,fail,1)
         close (20)
         goto 10
         open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown')
           enddo
         enddo
         do isample=1,MaxSample
-          call gen_side(it,90.0D0 * deg2rad,al,om,fail)
+          call gen_side(it,90.0D0 * deg2rad,al,om,fail,1)
           indal=rad2deg*al/2
           indom=(rad2deg*om+180.0D0)/5
           prob(indom,indal)=prob(indom,indal)+delt
       subroutine gen_rand_conf(nstart,*)
 ! Generate random conformation or chain cut and regrowth.
       use mcm_data
+      use random, only: iran_num,ran_number
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.CHAIN'
       integer :: it1,it2,it,j
 !d    print *,' CG Processor',me,' maxgen=',maxgen
       maxsi=1000
-!      write (iout,*) 'Gen_Rand_conf: nstart=',nstart
+      write (iout,*) 'Gen_Rand_conf: nstart=',nstart,nres
       if (nstart.lt.5) then
         it1=iabs(itype(2,1))
         phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1)))
 !       write(iout,*)'phi(4)=',rad2deg*phi(4)
-        if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4))
+        if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4),molnum(2))
 !       write(iout,*)'theta(3)=',rad2deg*theta(3) 
         if ((it1.ne.10).and.(it1.ne.ntyp1)) then
           nsi=0
           fail=.true.
           do while (fail.and.nsi.le.maxsi)
-            call gen_side(it1,theta(3),alph(2),omeg(2),fail)
+            call gen_side(it1,theta(3),alph(2),omeg(2),fail,molnum(2))
+            write (iout,*) 'nsi=',nsi,maxsi
             nsi=nsi+1
           enddo
           if (nsi.gt.maxsi) return 1
         endif ! it1.ne.10
+        write(iout,*) "before origin_frame"
         call orig_frame
+        write(iout,*) "after origin_frame"
         i=4
         nstart=4
       else
       niter=0
       back=.false.
       do while (i.le.nres .and. niter.lt.maxgen)
+        write(iout,*) 'i=',i,'back=',back
         if (i.lt.nstart) then
           if(iprint.gt.1) then
           write (iout,'(/80(1h*)/2a/80(1h*))') &
        it1=iabs(itype(i-1,molnum(i-1)))
        it2=iabs(itype(i-2,molnum(i-2)))
        it=iabs(itype(i,molnum(i)))
+        if ((it.eq.ntyp1).and.(it1.eq.ntyp1)) &
+          vbld(i)=ran_number(30.0D0,40.0D0)
 !       print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,&
 !        ' nit=',nit,' niter=',niter,' maxgen=',maxgen
        phi(i+1)=gen_phi(i+1,it1,it)
        if (back) then
           phi(i)=gen_phi(i+1,it2,it1)
 !         print *,'phi(',i,')=',phi(i)
-         theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
+         theta(i-1)=gen_theta(it2,phi(i-1),phi(i),molnum(i))
 !          print *,"theta",theta(i-1),phi(i)
          if ((it2.ne.10).and.(it2.ne.ntyp1)) then
             nsi=0
             fail=.true.
             do while (fail.and.nsi.le.maxsi)
-              call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
+              call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail,molnum(i-2))
               nsi=nsi+1
             enddo
             if (nsi.gt.maxsi) return 1
           endif
          call locate_next_res(i-1)
        endif
-       theta(i)=gen_theta(it1,phi(i),phi(i+1))
-        write(iout,*) "theta(i),",theta(i)
+       theta(i)=gen_theta(it1,phi(i),phi(i+1),molnum(i))
+!        write(iout,*) "theta(i),",theta(i)
         if ((it1.ne.10).and.(it1.ne.ntyp1)) then 
         nsi=0
         fail=.true.
         do while (fail.and.nsi.le.maxsi)
-          call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
+          call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail,molnum(i))
 !                  write(iout,*)"alpha,omeg(i-1)",alph(i-1),omeg(i-1),i,nsi,maxsi
           nsi=nsi+1
         enddo
            endif
           endif
         else
-          write(iout,*) "tu dochodze"
+!          write(iout,*) "tu dochodze"
          back=.false.
          nit=0 
          i=i+1
       nres2=2*nres
       data redfac /0.5D0/
       overlap=.false.
-      iti=iabs(itype(i,1))
+      iti=iabs(itype(i,molnum(i)))
       if (iti.gt.ntyp) return
 ! Check for SC-SC overlaps.
 !d    print *,'nnt=',nnt,' nct=',nct
       do j=nnt,i-1
+!        print *, "molnum(j)",j,molnum(j)
+        if (molnum(j).eq.1) then
         itj=iabs(itype(j,1))
+        if (itj.eq.ntyp1) cycle
         if (j.lt.i-1 .or. ipot.ne.4) then
           rcomp=sigmaii(iti,itj)
         else 
 !d      print *,'j=',j
        if (dist(nres+i,nres+j).lt.redfac*rcomp) then
           overlap=.true.
+        
 !        print *,'overlap, SC-SC: i=',i,' j=',j,
 !     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
 !     &     rcomp
          return
         endif
+        else if (molnum(j).eq.2) then
+        itj=iabs(itype(j,2))
+        if (dist(nres+i,nres+j).lt.redfac*sigma_nucl(iti,itj)) then
+          overlap=.true.
+
+!        print *,'overlap, SC-SC: i=',i,' j=',j,
+!     &     ' dist=',dist(nres+i,nres+j),' rcomp=',
+!     &     rcomp
+          return
+        endif
+        
+      endif
       enddo
 ! Check for overlaps between the added peptide group and the preceding
 ! SCs.
        c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
       enddo
       do j=nnt,i-2
+        if (molnum(j).ne.1) cycle
        itj=iabs(itype(j,1))
 !d      print *,'overlap, p-Sc: i=',i,' j=',j,
 !d   &         ' dist=',dist(nres+j,maxres2+1)
 ! Check for overlaps between the added side chain and the preceding peptide
 ! groups.
       do j=1,nnt-2
+        if (molnum(j).ne.1) cycle
        do k=1,3
          c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1))
         enddo
       enddo
 ! Check for p-p overlaps
       do j=1,3
-       c(j,nres2+4)=0.5D0*(c(j,i)+c(j,i+1))
+       c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
       enddo
       do j=nnt,i-2
+!        if (molnum(j).eq.1) then
         itelj=itel(j)
        do k=1,3
          c(k,nres2+4)=0.5D0*(c(k,j)+c(k,j+1))
         enddo
 !d      print *,'overlap, p-p: i=',i,' j=',j,
 !d   &         ' dist=',dist(maxres2+1,maxres2+2)
+        if (molnum(j).eq.1) then
         if(iteli.ne.0.and.itelj.ne.0)then
         if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then
           overlap=.true.
           return
         endif
         endif
+        else if (molnum(j).eq.2) then
+        if (dist(nres2+3,nres2+4).lt.3.0) then
+          overlap=.true.
+          return
+        endif
+      endif
       enddo
       return
       end function overlap
       return
       end function gen_phi
 !-----------------------------------------------------------------------------
-      real(kind=8) function gen_theta(it,gama,gama1)
-      use random,only:binorm
+      real(kind=8) function gen_theta(it,gama,gama1,mnum)
+      use random,only:binorm,ran_number
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'COMMON.LOCAL'
       real(kind=8),dimension(2) :: y,z
       real(kind=8) :: theta_max,theta_min,sig,ak
 !el local variables
-      integer :: j,it,k
+      integer :: j,it,k,mnum
       real(kind=8) :: gama,gama1,thet_pred_mean,theta_temp
 !     print *,'gen_theta: it=',it
       theta_min=0.05D0*pi
       else
        z(1)=0.0D0
        z(2)=0.0D0
-      endif  
+      endif 
+      if (it.eq.ntyp1) then
+      gen_theta=ran_number(theta_max/2.0,theta_max)
+      else if (mnum.eq.1) then
+             
       thet_pred_mean=a0thet(it)
-      write(iout,*),it,thet_pred_mean,"gen_thet"
+!      write(iout,*),it,thet_pred_mean,"gen_thet"
       do k=1,2
         thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) &
            +bthet(k,it,1,1)*z(k)
       if (theta_temp.gt.theta_max) theta_temp=theta_max
       gen_theta=theta_temp
 !     print '(a)','Exiting GENTHETA.'
+      else if (mnum.eq.2) then
+       gen_theta=2.0d0 + ran_number(0.0d0,0.34d0)
+      else
+              gen_theta=ran_number(theta_max/2.0,theta_max)
+       endif
       return
       end function gen_theta
 !-----------------------------------------------------------------------------
-      subroutine gen_side(it,the,al,om,fail)
+      subroutine gen_side(it,the,al,om,fail,mnum)
       use random, only:ran_number,mult_norm1
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
       real(kind=8) :: Big=10.0D0
       logical :: lprint,fail,lcheck
 !el local variables
-      integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob
+      integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob,mnum
       real(kind=8) :: the,al,om,detApi,wart,y2,wykl,radmax
       real(kind=8) :: tant,zz1,W1i,radius,zk,fac,dV,sum,sum1
       real(kind=8) :: which_lobe
       lcheck=.false.
       lprint=.false.
       fail=.false.
+      if (mnum.eq.1) then
       if (the.eq.0.0D0 .or. the.eq.pi) then
 #ifdef MPI
         write (*,'(a,i4,a,i3,a,1pe14.5)') &
         fail=.true.
         return
       endif
+      if (nlobit.eq.0) then
+         al=ran_number(0.05d0,pi/2)
+         om=ran_number(-pi,pi)
+         return
+      endif
       tant=dtan(the-pipol)
       nlobit=nlob(it)
       allocate(z(3,nlobit))
         write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
         write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
       endif
-      if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
-#ifdef MPI
-        write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax
-        write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
-#else
+!      if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
+!#ifdef MPI
+!        write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax
+!        write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
+!#else
 !        write (iout,'(a)') 'Bad sampling box.'
-#endif
-        fail=.true.
-        return
-      endif
+!#endif
+!        fail=.true.
+!        return
+!      endif
       which_lobe=ran_number(0.0D0,sumW(nlobit))
 !     print '(a,1pe14.4)','which_lobe=',which_lobe
       do i=1,nlobit
       if (fail) return
       al=y(1)
       om=pinorm(y(2))
+      else if (mnum.eq.2) then
+       al=0.7+ran_number(0.0d0,0.2d0)
+       om=ran_number(0.0d0,3.14d0)
+      endif
+      
 !d    print *,'al=',al,' om=',om
 !d    stop
       return
         do ires=1,ioverlap_last 
           i=ioverlap(ires)
           iti=iabs(itype(i,1))
-          if ((iti.ne.10).and.(molnum(i).ne.5).and.(iti.ne.ntyp1)) then
+          if ((iti.ne.10).and.(molnum(i).lt.3).and.(iti.ne.ntyp1)) then
             nsi=0
             fail=.true.
             do while (fail.and.nsi.le.maxsi)
-              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
+              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
               nsi=nsi+1
             enddo
             if(fail) goto 999
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
+        print *,i,itypi,"sc_move"
         dsci_inv=dsc_inv(itypi)
 !
        do iint=1,nint_gr(i)
          if (itype(j,molnum(j)).eq.ntyp1_molec(molnum(j))) cycle
             ind=ind+1
             itypj=iabs(itype(j,molnum(j)))
+        print *,j,itypj,"sc_move"
+
             dscj_inv=dsc_inv(itypj)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)
       chiom1=chi1*om1
       chiom2=chi2*om2
       facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
+!      print *,"TUT?",om1*chiom1,facsig,om1,om2,om12
       sigsq=1.0D0-facsig*faceps1_inv
       sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
       sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
       integer :: total_ints,lower_bound,upper_bound,nint
       integer,dimension(0:nfgtasks) :: int4proc,sint4proc      !(0:max_fg_procs)
       integer :: i,nexcess
+      if (total_ints.le.0) then
+      lower_bound=1
+      upper_bound=0
+      return
+      endif
       nint=total_ints/nfgtasks
       do i=1,nfgtasks
         int4proc(i-1)=nint
       character(len=3) :: seq,res
 !      character*5 atom
       character(len=80) :: card
-      real(kind=8),dimension(3,20) :: sccor
+      real(kind=8),dimension(3,40) :: sccor
       integer :: i,j,iti !el  rescode,
       logical :: lside,lprn
       real(kind=8) :: di,cosfac,sinfac
 !      include 'DIMENSIONS'
 !      include 'COMMON.CHAIN'
       integer :: i,j,ires,nscat
-      real(kind=8),dimension(3,20) :: sccor
+      real(kind=8),dimension(3,40) :: sccor
       real(kind=8) :: sccmj
 !        print *,"I am in sccenter",ires,nscat
       do j=1,3
        do j=1,3
          gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
            +gloc(nres-2,icg)*dtheta(j,1,3)      
+!         write(iout,*) "pierwszy gcart", gcart(j,2)
           if ((itype(2,1).ne.10).and.&
-          (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+          (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).lt.3))) then
           gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
           gloc(ialph(2,1)+nside,icg)*domega(j,1,2)             
         endif
        do j=1,3
          gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
         gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
-        if(itype(2,1).ne.10) then
+!         write(iout,*) "drugi gcart", gcart(j,2)
+        if((itype(2,1).ne.10).and.(molnum(2).lt.3)) then
           gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
           gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
         endif
           +gloc(i-1,icg)*dphi(j,2,i+2)+ &
           gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
           gloc(nres+i-3,icg)*dtheta(j,1,i+2)
-          if(itype(i,1).ne.10) then
+          if((itype(i,1).ne.10).and.(molnum(nres-1).lt.3)) then
            gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
            gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
           endif
           dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
            +gloc(2*nres-6,icg)* &
            dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
-          if(itype(nres-2,1).ne.10) then
+          if((itype(nres-2,1).ne.10).and.(molnum(nres-1).lt.3)) then
               gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
              dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
               domega(j,2,nres-2)
           endif
-          if(itype(nres-1,1).ne.10) then
+          if((itype(nres-1,1).ne.10).and.(molnum(nres-1).lt.3)) then
              gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
             dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
              domega(j,1,nres-1)
       endif 
 !  Settind dE/ddnres-1       
 !#define DEBUG
-#ifdef DEBUG
-          j=1
-              write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
-        gloc(2*nres-5,icg),dtheta(j,2,nres)
+!#ifdef DEBUG
+!          j=1
+!              write(iout,*)"in int to carta",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
+!        gloc(2*nres-5,icg),dtheta(j,2,nres)
 
-#endif
+!#endif
 !#undef DEBUG
 
        do j=1,3
         gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
        gloc(2*nres-5,icg)*dtheta(j,2,nres)
 !#define DEBUG
-#ifdef DEBUG
-              write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
-        gloc(2*nres-5,icg),dtheta(j,2,nres)
-
-#endif
+!#ifdef DEBUG
+!              write(iout,*)"in int to cartb",nres-1,gcart(j,nres-1),gloc(nres-3,icg),dphi(j,3,nres), &
+!        gloc(2*nres-5,icg),dtheta(j,2,nres)
+!
+!#endif
 !#undef DEBUG
-        if(itype(nres-1,1).ne.10) then
+        if((itype(nres-1,1).ne.10).and.(molnum(nres-1).lt.3)) then
           gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
          dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
           domega(j,2,nres-1)
 !#define DEBUG
-#ifdef DEBUG
-              write(iout,*)"in int to cart2",i,gcart(j,nres-1),gloc(ialph(nres-1,1),icg)* &
-          dalpha(j,2,nres-1),gloc(ialph(nres-1,1)+nside,icg), &
-          domega(j,2,nres-1)
+!#ifdef DEBUG
+!              write(iout,*)"in int to cart2",i,gcart(j,nres-1),gloc(ialph(nres-1,1),icg)* &
+!          dalpha(j,2,nres-1),gloc(ialph(nres-1,1)+nside,icg), &
+!          domega(j,2,nres-1)
 
-#endif
+!#endif
 !#undef DEBUG
 
         endif
 !   The side-chain vector derivatives
         do i=2,nres-1
          if(itype(i,1).ne.10 .and.  &
-           itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) then  
+           itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)).and.&
+             (molnum(i).lt.3)) then    
             do j=1,3   
               gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
               +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
 !#define DEBUG
-#ifdef DEBUG
-              write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), &
-              gloc(ialph(i,1)+nside,icg),domega(j,3,i)
-#endif
+!#ifdef DEBUG
+!              write(iout,*)"in int to cart",i, gxcart(j,i),gloc(ialph(i,1),icg),dalpha(j,3,i), &
+!              gloc(ialph(i,1)+nside,icg),domega(j,3,i)
+!#endif
 !#undef DEBUG
             enddo
          endif     
 ! INTERTYP=3 SC...Ca...Ca...SC
 !   calculating dE/ddc1      
   18   continue
+!         write(iout,*) "przed sccor gcart", gcart(1,2)
+
 !       do i=1,nres
 !       gloc(i,icg)=0.0D0
 !          write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
        if (nres.lt.2) return
        if ((nres.lt.3).and.(itype(1,1).eq.10)) return
        if ((itype(1,1).ne.10).and. &
-        (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
+        (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).lt.3)) then
         do j=1,3
 !c Derviative was calculated for oposite vector of side chain therefore
 ! there is "-" sign before gloc_sc
          gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
            dtauangle(j,1,2,3)
           if ((itype(2,1).ne.10).and. &
-        (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+        (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2))).and.(molnum(2).lt.3)) then
          gxcart(j,1)= gxcart(j,1) &
                      -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
          gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
 !   ommited 
 !     &     +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)  
 
+!         write(iout,*) "przed dE/ddc2 gcart", gcart(1,2)
+
 !     Calculating the remainder of dE/ddc2
        do j=1,3
          if((itype(2,1).ne.10).and. &
-           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+           (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)).and.(molnum(2).lt.3))) then
            if ((itype(1,1).ne.10).and.&
-              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))))&
+              ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))).and.(molnum(1).lt.3))&
             gxcart(j,2)=gxcart(j,2)+ &
                                gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
-        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3))) &
+        if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(3)).and.molnum(3).lt.3) &
          then
            gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
 !c                  the   - above is due to different vector direction
            gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
 !c                  the   - above is due to different vector direction
            gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
-!          write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart"
+!          write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart",gcart(j,2)
 !           write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
           endif
          endif
          if ((itype(1,1).ne.10).and.&
-         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) then
+         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).lt.3)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
 !           write(iout,*)  gloc_sc(1,0,icg),dtauangle(j,1,3,3)
         endif
-         if ((itype(3,1).ne.10).and.(nres.ge.3)) then
+         if ((itype(3,1).ne.10).and.(nres.ge.3).and.(molnum(3).lt.3)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
 !           write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
          endif
-         if ((itype(4,1).ne.10).and.(nres.ge.4)) then
+         if ((itype(4,1).ne.10).and.(nres.ge.4).and.(molnum(4).lt.3)) then
           gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
 !           write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
          endif
 
 !      write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
        enddo
+!         write(iout,*) "po dE/ddc2 gcart", gcart(1,2)
+
 !    If there are more than five residues
       if(nres.ge.5) then                        
         do i=3,nres-2
          do j=1,3
 !          write(iout,*) "before", gcart(j,i)
           if ((itype(i,1).ne.10).and.&
-          (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) then
+          (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))).and.(molnum(i).lt.3)) then
           gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
           *dtauangle(j,2,3,i+1) &
           -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
 !     &  gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
 !          if (itype(i-1,1).ne.10) then
           if ((itype(i-1,1).ne.10).and.&
-          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then
+          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.(molnum(i-1).eq.5)) then
 
            gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
       *dtauangle(j,3,3,i+1)
           endif
 !          if (itype(i+1,1).ne.10) then
           if ((itype(i+1,1).ne.10).and.&
-          (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then
+          (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1))).and.&
+           (molnum(i+1).lt.3)) then
           gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
       *dtauangle(j,3,1,i+2)
            gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
           endif
 !          if (itype(i-1,1).ne.10) then
           if ((itype(i-1,1).ne.10).and.&
-          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1)))) then
+          (itype(i-1,molnum(i-1)).ne.ntyp1_molec(molnum(i-1))).and.&
+         (molnum(i-1).lt.3)) then
            gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
            dtauangle(j,1,3,i+1)
           endif
 !          if (itype(i+1,1).ne.10) then
           if ((itype(i+1,1).ne.10).and.&
-          (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))) then
+          (itype(i+1,molnum(i+1)).ne.ntyp1_molec(molnum(i+1)))&
+       .and. (molnum(i+1).lt.3)) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
            dtauangle(j,2,2,i+2)
 !          write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
           endif
 !          if (itype(i+2,1).ne.10) then
           if ((itype(i+2,1).ne.10).and.&
-          (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2)))) then
+          (itype(i+2,molnum(i+2)).ne.ntyp1_molec(molnum(i+2))).and.(molnum(i+2).lt.3)) then
            gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
            dtauangle(j,2,1,i+3)
           endif
       if(nres.ge.4) then
          do j=1,3
          if ((itype(nres-1,1).ne.10).and.&
-       (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) then
+       (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1))).and.(molnum(nres-1).lt.3)) then
          gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
           *dtauangle(j,2,3,nres)
 !          write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
 !     &     dtauangle(j,2,3,nres), gxcart(j,nres-1)
 !         if (itype(nres-2,1).ne.10) then
          if ((itype(nres-2,1).ne.10).and.&
-       (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
+       (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).lt.3)) then
        gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
           *dtauangle(j,3,3,nres)
           endif
          if ((itype(nres,1).ne.10).and.&
-         (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
+         (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).lt.3)) then
         gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
           *dtauangle(j,3,1,nres+1)
         gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
           endif
          endif
          if ((itype(nres-2,1).ne.10).and.&
-         (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
+         (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2))).and.(molnum(nres-2).lt.3)) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
          dtauangle(j,1,3,nres)
          endif
-          if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) then
+          if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).lt.3)) then
             gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
            dtauangle(j,2,2,nres+1)
 !           write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
       endif
 !  Settind dE/ddnres       
        if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
-          (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))then
+          (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))).and.(molnum(nres).lt.3))then
        do j=1,3
         gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
        *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
        *dtauangle(j,2,3,nres+1)
         enddo
        endif
+!       write(iout,*) "final gcart",gcart(1,2)
 !   The side-chain vector derivatives
 !       print *,"gcart",gcart(:,:)
       return
                 '     centroid coordinates'/ &
                 '       ', 6X,'X',11X,'Y',11X,'Z',&
                                 10X,'X',11X,'Y',11X,'Z')
-  110 format (a,'(',i3,')',6f12.5)
+  110 format (a,'(',i6,')',6f12.5)
       return
       end subroutine cartprint
 !-----------------------------------------------------------------------------
 !      common /rotmat/
       allocate(t(3,3,nres),r(3,3,nres))
       allocate(prod(3,3,nres),rt(3,3,nres)) !(3,3,maxres)
+      print *,"before permutations",maxperm,maxsym
 !      common /refstruct/
       if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
+      print *,"cref"
 !elwrite(iout,*) "jestem w alloc geo 2"
-      allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
+!      allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
+        if (.not.allocated(crefjlee)) allocate (crefjlee(3,nres2+2))
       if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym)
+      print *,"chain_rep"
       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
 !      common /from_zscore/ in module.compare
 !----------------------
 !-----------------------------------------------------------------------------
       subroutine returnbox
       integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
-      integer :: chain_end,ireturnval
-      real*8 :: difference
+      integer :: chain_end,ireturnval,idum,mnumi1
+      real*8 :: difference,xi,boxsize,x,xtemp,box2shift
+      real(kind=8),dimension(3) :: boxx
+      real(kind=8),dimension(3,10000) :: xorg
+      integer,dimension(10000) :: posdummy
+
 !C change suggested by Ana - end
         j=1
         chain_beg=1
-!C        do i=1,nres
-!C       write(*,*) 'initial', i,j,c(j,i)
-!C        enddo
+        boxx(1)=boxxsize
+        boxx(2)=boxysize
+        boxx(3)=boxzsize
+        idum=0
+!        if(me.eq.king.or..not.out1file) then
+!        do i=1,nres
+!          write(iout,'(a6,2i3,6f9.3)') 'initial', i,j,(c(j,i),j=1,3),(c(j,i+nres),j=1,3)
+!        enddo
+!        endif
 !C change suggested by Ana - begin
         allareout=1
+        chain_end=0
 !C change suggested by Ana -end
-        do i=1,nres-1
-           mnum=molnum(i)
-         if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
-            .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
-          chain_end=i
-          if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxxsize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,chain_end
-              c(j,k)=c(j,k)-ireturnval*boxxsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
-            enddo
-!C Suggested by Ana
-            if (chain_beg.eq.1) &
-            dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
-!C Suggested by Ana -end
-           endif
-           chain_beg=i+1
-           allareout=1
-         else
-          if (int(c(j,i)/boxxsize).eq.0) allareout=0
-         endif
-        enddo
-         if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxxsize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,nres
-              c(j,k)=c(j,k)-ireturnval*boxxsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
-            enddo
-          endif
-!C NO JUMP 
-!C        do i=1,nres
-!C        write(*,*) 'befor no jump', i,j,c(j,i)
-!C        enddo
-        nojumpval=0
-        do i=2,nres
-           mnum=molnum(i)
-           if (itype(i,mnum).eq.ntyp1_molec(mnum)&
-              .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
-             difference=abs(c(j,i-1)-c(j,i))
-!C             print *,'diff', difference
-             if (difference.gt.boxxsize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-              endif
-              c(j,i)=c(j,i)+nojumpval*boxxsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
+        do i=1,nres
+          mnum=molnum(i)
+         if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
+         idum=idum+1
+         posdummy(idum)=i
+         if (i.ne.nres) then
+         mnumi1=molnum(i+1)
+         if ((itype(i+1,mnum).eq.ntyp1_molec(mnum)).or.(mnum.ne.mnumi1)) then
+         do j=1,3
+         xorg(j,idum)=c(j,i)-c(j,i-1)
          enddo
-       nojumpval=0
-        do i=2,nres
-           mnum=molnum(i)
-           if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxxsize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-             endif
-              c(j,i)=c(j,i)+nojumpval*boxxsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
+         else
+         do j=1,3
+         xorg(j,idum)=c(j,i)-c(j,i+1)
          enddo
-
-!C        do i=1,nres
-!C        write(*,*) 'after no jump', i,j,c(j,i)
-!C        enddo
-
-!C NOW Y dimension
-!C suggesed by Ana begins
-        allareout=1
-        j=2
-        chain_beg=1
-        do i=1,nres-1
-           mnum=molnum(i)
-         if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
-           .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
-          chain_end=i
-          if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxysize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,chain_end
-              c(j,k)=c(j,k)-ireturnval*boxysize
-             c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
-            enddo
-!C Suggested by Ana
-            if (chain_beg.eq.1) &
-            dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
-!C Suggested by Ana -end
-           endif
-           chain_beg=i+1
-           allareout=1
+         endif
          else
-          if (int(c(j,i)/boxysize).eq.0) allareout=0
+         do j=1,3
+         xorg(j,idum)=c(j,i)-c(j,i-1)
+         enddo
+         endif
          endif
         enddo
-         if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxysize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,nres
-              c(j,k)=c(j,k)-ireturnval*boxysize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
-            enddo
+ 12     continue
+        do i=1,nres
+         mnum=molnum(i)
+         if (molnum(i).ge.1) then
+          if (i.le.3) then
+           k=2
+          else
+          if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
+          k=k+1
           endif
-        nojumpval=0
-        do i=2,nres
-           mnum=molnum(i)
-           if (itype(i,mnum).eq.ntyp1_molec(mnum)&
-              .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxysize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-             else
-              nojumpval=0
-              endif
-           endif
-              c(j,i)=c(j,i)+nojumpval*boxysize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
+          endif
+!          print *,"tu2",i,k
+          if (itype(k,mnum).eq.ntyp1_molec(mnum)) k=k+1
+          if (itype(k,mnum).eq.ntyp1_molec(mnum)) k=k+1
+!          print *,"tu2",i,k
+
+          do j=1,3
+          c(j,i)=dmod(c(j,i),boxx(j))
+!          if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
+          x=c(j,i)-c(j,k)
+!          print *,"return box,before wrap",c(1,i)
+         boxsize=boxx(j)
+           xtemp=dmod(x,boxsize)
+      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+        box2shift=xtemp-boxsize
+      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+        box2shift=xtemp+boxsize
+      else
+        box2shift=xtemp
+      endif
+          xi=box2shift
+!          print *,c(1,2),xi,xtemp,box2shift
+          c(j,i)=c(j,k)+xi
          enddo
-      nojumpval=0
-        do i=2,nres
-           mnum=molnum(i)
-           if (itype(i,mnum).eq.ntyp1_molec(mnum)&
-             .and. itype(i-1,mnum).eq.ntyp1) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxysize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-            endif
-              c(j,i)=c(j,i)+nojumpval*boxysize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
+          do j=1,3
+          c(j,i+nres)=dmod(c(j,i+nres),boxx(j))
+!          if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
+          x=c(j,i+nres)-c(j,i)
+!          print *,"return box,before wrap",c(1,i)
+         boxsize=boxx(j)
+           xtemp=dmod(x,boxsize)
+      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+        box2shift=xtemp-boxsize
+      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+        box2shift=xtemp+boxsize
+      else
+        box2shift=xtemp
+      endif
+          xi=box2shift
+!          print *,c(1,2),xi,xtemp,box2shift
+          c(j,i+nres)=c(j,i)+xi
          enddo
-!C Now Z dimension
-!C Suggested by Ana -begins
-        allareout=1
-!C Suggested by Ana -ends
-       j=3
-        chain_beg=1
-        do i=1,nres-1
-           mnum=molnum(i)
-         if ((itype(i,mnum).eq.ntyp1_molec(mnum))&
-           .and.(itype(i+1,mnum).eq.ntyp1_molec(mnum))) then
-          chain_end=i
-          if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxysize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,chain_end
-              c(j,k)=c(j,k)-ireturnval*boxzsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
-            enddo
-!C Suggested by Ana
-            if (chain_beg.eq.1) dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
-!C Suggested by Ana -end
-           endif
-           chain_beg=i+1
-           allareout=1
-         else
-          if (int(c(j,i)/boxzsize).eq.0) allareout=0
          endif
         enddo
-         if (allareout.eq.1) then
-            ireturnval=int(c(j,i)/boxzsize)
-            if (c(j,i).le.0) ireturnval=ireturnval-1
-            do k=chain_beg,nres
-              c(j,k)=c(j,k)-ireturnval*boxzsize
-              c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
+        do i=1,idum
+          k=posdummy(i)
+          mnum=molnum(k)
+        if(me.eq.king.or..not.out1file) then
+          write(iout,*),"posdummy",i,k,(xorg(j,i),j=1,3)
+        endif
+!          do j=1,3
+!           if (dabs(xorg(j,i)).gt.boxx(j))then
+!           x=xorg(j,i)
+!           boxsize=boxx(j)
+!           xtemp=dmod(x,boxsize)
+!           if (dabs(-boxsize).lt.dabs(xtemp)) then
+!             box2shift=xtemp-boxsize
+!           else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+!          box2shift=xtemp+boxsize
+!!        else
+!           box2shift=xtemp
+!          endif
+!           xorg(j,i)=box2shift
+!           endif
+!          enddo
+          if (k.eq.nres) then
+           do j=1,3
+           c(j,k)=c(j,k-1)+xorg(j,i)
+           enddo
+          else
+           mnumi1=molnum(k+1)
+           if ((itype(k+1,mnum).eq.ntyp1_molec(mnum)).or.(mnum.ne.mnumi1)) then
+           do j=1,3
+             c(j,k)=c(j,k-1)+xorg(j,i)
+           enddo
+           else
+            do j=1,3
+             c(j,k)=c(j,k+1)+xorg(j,i)
             enddo
-          endif
-        nojumpval=0
-        do i=2,nres
-           mnum=molnum(i)
-           if (itype(i,mnum).eq.ntyp1_molec(mnum) .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.(boxzsize/2.0)) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-            endif
-              c(j,i)=c(j,i)+nojumpval*boxzsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
-         enddo
-       nojumpval=0
-        do i=2,nres
-           mnum=molnum(i)
-           if (itype(i,mnum).eq.ntyp1_molec(mnum) &
-            .and. itype(i-1,mnum).eq.ntyp1_molec(mnum)) then
-             difference=abs(c(j,i-1)-c(j,i))
-             if (difference.gt.boxzsize/2.0) then
-                if (c(j,i-1).gt.c(j,i)) then
-                  nojumpval=1
-                 else
-                   nojumpval=-1
-                 endif
-              else
-              nojumpval=0
-              endif
-            endif
-             c(j,i)=c(j,i)+nojumpval*boxzsize
-              c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
-         enddo
-        do i=1,nres
-         if (molnum(i).eq.5) then
-          c(1,i)=dmod(c(1,i),boxxsize)
-          c(2,i)=dmod(c(2,i),boxysize)
-          c(3,i)=dmod(c(3,i),boxzsize)
-          c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
-          c(2,i+nres)=dmod(c(2,i+nres),boxysize)
-          c(3,i+nres)=dmod(c(3,i+nres),boxzsize)
+           endif
          endif
         enddo
+!        if(me.eq.king.or..not.out1file) then
+!        do i=1,nres
+!          write(iout,'(a6,2i3,6f9.3)') 'final', i,j,(c(j,i),j=1,3),(c(j,i+nres),j=1,3)
+!        enddo
+!        endif
         return
         end       subroutine returnbox
 !-------------------------------------------------------------------------------------------------------