changes need for chamm-gui
[unres4.git] / source / unres / geometry.F90
index c89442f..60c888f 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).
 !        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).and.(molnum(i).ne.5)) 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
         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)
       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
            +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)).and.(molnum(2).ne.5))) 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
          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)
 !         write(iout,*) "drugi gcart", gcart(j,2)
-        if((itype(2,1).ne.10).and.(molnum(2).ne.5)) then
+        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).and.(molnum(nres-1).ne.5)) 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).and.(molnum(nres-1).ne.5)) 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).and.(molnum(nres-1).ne.5)) 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)
          enddo
       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)
+!#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)
 
-#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).and.(molnum(nres-1).ne.5)) 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
         do i=2,nres-1
          if(itype(i,1).ne.10 .and.  &
            itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)).and.&
-             (molnum(i).ne.5)) then    
+             (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
-#undef 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
+!#undef DEBUG
             enddo
          endif     
        enddo                                                                                                                                                   
        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))).and.(molnum(1).ne.5)) 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))).and.(molnum(2).ne.5)) 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)* &
 !     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)).and.(molnum(2).ne.5))) 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)))).and.(molnum(1).ne.5))&
+              ((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)).and.molnum(3).ne.5) &
+        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
           endif
          endif
          if ((itype(1,1).ne.10).and.&
-         (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))).and.(molnum(1).ne.5)) 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).and.(molnum(3).ne.5)) 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).and.(molnum(4).ne.5)) 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
          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))).and.(molnum(i).ne.5)) 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)
 !          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))).and.&
-           (molnum(i+1).ne.5)) then
+           (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) &
 !          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))).and.&
-         (molnum(i-1).ne.5)) then
+         (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)))&
-       .and. (molnum(i+1).ne.5)) then
+       .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))).and.(molnum(i+2).ne.5)) 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))).and.(molnum(nres-1).ne.5)) 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))).and.(molnum(nres-2).ne.5)) 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))).and.(molnum(nres).ne.5)) 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))).and.(molnum(nres-2).ne.5)) 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))).and.(molnum(nres).ne.5)) 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))).and.(molnum(nres).ne.5))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) &
 !      common /refstruct/
       if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
 !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)
       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)
-          if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
-          c(2,i)=dmod(c(2,i),boxysize)
-          if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize
-          c(3,i)=dmod(c(3,i),boxzsize)
-          if (c(3,i).lt.0) c(3,i)=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
 !-------------------------------------------------------------------------------------------------------