changes in wham and unres
[unres4.git] / source / unres / geometry.f90
index 52c3b83..2287f0a 100644 (file)
@@ -1,4 +1,4 @@
-             module geometry
+                     module geometry
 !-----------------------------------------------------------------------------
       use io_units
       use names
           endif
           return 1
         endif
-       it1=iabs(itype(i-1,1))
-       it2=iabs(itype(i-2,1))
-       it=iabs(itype(i,1))
+       it1=iabs(itype(i-1,molnum(i-1)))
+       it2=iabs(itype(i-2,molnum(i-2)))
+       it=iabs(itype(i,molnum(i)))
 !       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)
 !      print *,'>>overlap_sc nnt=',nnt,' nct=',nct
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i,1))
+        if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) cycle
+        itypi=iabs(itype(i,molnum(i)))
         itypi1=iabs(itype(i+1,1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
 !
        do iint=1,nint_gr(i)
          do j=istart(i,iint),iend(i,iint)
+         if (itype(j,molnum(j)).eq.ntyp1_molec(molnum(j))) cycle
             ind=ind+1
-            itypj=iabs(itype(j,1))
+            itypj=iabs(itype(j,molnum(j)))
             dscj_inv=dsc_inv(itypj)
             sig0ij=sigma(itypi,itypj)
             chi1=chi(itypi,itypj)
 !      include 'COMMON.INTERACT'
 !      include 'COMMON.CHAIN'
       do i=1,nres-1
+       
        vbld(i+1)=vbl
        vbld_inv(i+1)=1.0d0/vbld(i+1)
-       vbld(i+1+nres)=dsc(itype(i+1,1))
-       vbld_inv(i+1+nres)=dsc_inv(itype(i+1,1))
+       vbld(i+1+nres)=dsc(itype(i+1,molnum(i)))
+       vbld_inv(i+1+nres)=dsc_inv(itype(i+1,molnum(i)))
 !       print *,vbld(i+1),vbld(i+1+nres)
       enddo
       return
       end subroutine alloc_geo_arrays
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
+      subroutine returnbox
+      integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
+      integer :: chain_end,ireturnval
+      real*8 :: difference
+!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
+!C change suggested by Ana - begin
+        allareout=1
+!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
+         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
+         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
+         else
+          if (int(c(j,i)/boxysize).eq.0) allareout=0
+         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
+          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
+         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
+         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
+            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
+        enddo
+        return
+        end       subroutine returnbox
+!-------------------------------------------------------------------------------------------------------
       end module geometry