xdrf2pdb Adam's PBC correction

[unres.git] / source / xdrfpdb / src-M / xdrf2pdb-m.F

diff --git a/source/xdrfpdb/src-M/xdrf2pdb-m.F b/source/xdrfpdb/src-M/xdrf2pdb-m.F
index d2d7d9c..e3825b0 100644 (file)
--- a/source/xdrfpdb/src-M/xdrf2pdb-m.F
+++ b/source/xdrfpdb/src-M/xdrf2pdb-m.F
@@ -16,12 +16,15 @@
       logical oneletter,iblnk
       integer rescode
       external rescode
+      integer lshift(maxchain),lmin(maxchain)
+      double precision shift_coord(3,maxchain),cm(3,maxchain),ccm(3),
+     & ccm_prev(3)
       
       ifreq=1
-      if (iargc().lt.5) then
+      if (iargc().lt.8) then
         print '(2a)',
-     &   "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",
-     &   " [pdbfile] [freq]"
+     &   "Usage: xdrf2pdb-m one/three seqfile boxx boxy boxz cxfile ",
+     &   "ntraj itraj [pdbfile] [freq]"
         stop
       endif
       call getarg(1,onethree)
@@ -62,6 +65,22 @@
         print *,nres
         print '(20(a3,1x))',(sequenc(i),i=1,nres)
       endif
+      call seq2chains(nres,itype,nchain,chain_length,chain_border,
+     &  ireschain)
+      print *,"nchain",nchain
+      print *,"chain_border, chain_length"
+      nrestot=0
+      do i=1,nchain
+        nrestot=nrestot+chain_length(i)
+        print *,chain_border(1,i),chain_border(2,i),chain_length(i)
+      enddo
+      print *,"nrestot",nrestot
+      call getarg(3,arg)
+      read(arg,*) boxsize(1)
+      call getarg(4,arg)
+      read(arg,*) boxsize(2)
+      call getarg(5,arg)
+      read(arg,*) boxsize(3)
       call getarg(3,arg)
       iext = index(arg,'.cx') - 1
       if (iext.lt.0) then
@@ -137,6 +156,117 @@ c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
              c(j,i+nres)=coord(j,ii+nres)
            enddo
          enddo
+c Bring all chains to same box and to a common origin
+         do ichain=1,nchain
+c           print *,"ichain",ichain
+           do j=1,3
+             cm(j,ichain)=0.0d0
+           enddo
+           do i=chain_border(1,ichain),chain_border(2,ichain)
+c             print *,i,c(:,i)
+             do j=1,3
+               cm(j,ichain)=cm(j,ichain)+c(j,i)
+             enddo
+           enddo
+           do j=1,3
+             cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
+           enddo
+         enddo
+c         write (*,*) "CM"
+c         do i=1,nchain
+c           write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
+c         enddo
+         do j=1,3
+           lmin=0
+           sumodl_min=1.0d10
+           do i=0,3**nchain-1
+             ii=i
+             do ichain=1,nchain
+               lshift(ichain)=mod(ii,3)-1
+               ii=ii/3
+             enddo
+c             print '(10i3)',(lshift(ichain),ichain=1,nchain)
+             sumodl=0.0d0
+             do ichain=1,nchain
+               do jchain=1,ichain-1
+               sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
+     &             -cm(j,jchain)-lshift(jchain)*boxsize(j))
+               enddo
+             enddo
+c             print *,"sumodl",sumodl," sumodl_min",sumodl_min
+             if (sumodl.lt.sumodl_min) then
+               sumodl_min=sumodl
+               lmin=lshift
+             endif
+           enddo
+           do ichain=1,nchain
+             cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
+             shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
+           enddo
+         enddo
+         do ichain=1,nchain
+           do j=1,3
+             ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
+           enddo
+         enddo
+         do j=1,3
+           ccm(j)=ccm(j)/nrestot
+         enddo
+c         write (*,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
+c         write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
+c         write (*,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
+c         write (*,*) "shift_coord"
+c         do i=1,nchain
+c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+c         enddo
+#ifdef PREVIOUS
+         if (previous) then
+           do j=1,3
+             change=.true.
+             do while (change)
+             if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
+               do ichain=1,nchain
+                 shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
+               enddo
+               ccm(j)=ccm(j)-boxsize(j)
+             else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
+               do ichain=1,nchain
+                 shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
+               enddo
+               ccm(j)=ccm(j)+boxsize(j)
+             else
+               change=.false.
+             endif
+           enddo
+           enddo
+           do j=1,3
+             ccm_prev(j)=ccm(j)
+           enddo
+         else
+           previous=.true.
+           do j=1,3
+             ccm_prev(j)=ccm(j)
+           enddo
+         endif
+#else
+         do ichain=1,nchain
+           do j=1,3
+             shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
+           enddo
+         enddo 
+#endif
+c         write (*,*) "shift_coord"
+c         do i=1,nchain
+c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+c         enddo
+         do ichain=1,nchain
+           do i=chain_border(1,ichain),chain_border(2,ichain)
+             do j=1,3
+               c(j,i)=c(j,i)+shift_coord(j,ichain)
+               c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
+             enddo
+           enddo
+         enddo 
          etot=potE
          write (tytul,'(a1,i6,a8,f6.1,a6,f10.1)') "#",kk," t_bath ",
      &             t_bath, " time ",time