include "COMMON.CHAIN"
include "COMMON.INTERACT"
include "COMMON.IOUNITS"
- integer i,ii,ipi,ipj,ipmin,j,jmin,k,ix,iy,iz,
- & ixmin,iymin,izmin,ir_start,ir_end
- integer iper(maxchain),iaux
- double precision dchain,dchainmin,cmchain(3,20)
- cmchain=0.0d0
- do i=1,nchain
- ii=0
- do j=chain_border(1,i),chain_border(2,i)
- if (itype(j).eq.ntyp1) cycle
- ii=ii+1
- do k=1,3
- cmchain(k,i)=cmchain(k,i)+c(k,j)
+ integer i,j,k,l,ichain,jchain,ii,lshift(maxchain),lmin(maxchain),
+ & nrestot
+ double precision sumodl,sumodl_min,shift_coord(3,maxchain),
+ & boxsize(3),cm(3,maxchain),ccm(3),ccm_prev(3)
+ logical previous,change
+ boxsize(1)=boxxsize
+ boxsize(2)=boxysize
+ boxsize(3)=boxzsize
+ nrestot=0
+ do ichain=1,nchain
+c print *,"ichain",ichain
+ nrestot=nrestot+chain_length(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 k=1,3
- cmchain(k,i)=cmchain(k,i)/ii
+ do j=1,3
+ cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
enddo
enddo
+#ifdef DEBUG
+ write (iout,*) "CM before shift"
do i=1,nchain
- iper(i)=i
+ write (iout,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
enddo
- do i=1,nchain-1
- dchainmin=1.0d10
- do j=i+1,nchain
- ipi=iper(i)
- ipj=iper(j)
- do ix=-1,1
- do iy=-1,1
- do iz=-1,1
- dchain=(cmchain(1,ipj)-cmchain(1,ipi)+ix*boxxsize)**2+
- & (cmchain(2,ipj)-cmchain(2,ipi)+iy*boxysize)**2+
- & (cmchain(3,ipj)-cmchain(3,ipi)+iz*boxzsize)**2
-c write (iout,*) "i",i," ipi",ipi," j",j," ipj",ipj," d",
-c & dsqrt(dchain)," dmin",dsqrt(dchainmin)," jmin",jmin
- if (dchain.lt.dchainmin) then
- dchainmin=dchain
- ixmin=ix
- iymin=iy
- izmin=iz
- jmin=j
- endif
- enddo
+#endif DEBUG
+ 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
+#ifdef DEBUG
+ write (iout,'(10i3)') (lshift(ichain),ichain=1,nchain)
+#endif
+ 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
+#ifdef DEBUG
+ write(iout,*) "sumodl",sumodl," sumodl_min",sumodl_min
+#endif
+ if (sumodl.lt.sumodl_min) then
+ sumodl_min=sumodl
+ lmin=lshift
+ endif
enddo
- if (ixmin.eq.0 .and. iymin.eq.0 .and. izmin.eq.0) cycle
- ipj=iper(jmin)
- cmchain(1,ipj)=cmchain(1,ipj)+ixmin*boxxsize
- cmchain(2,ipj)=cmchain(2,ipj)+iymin*boxysize
- cmchain(3,ipj)=cmchain(3,ipj)+izmin*boxzsize
- ir_start=chain_border(1,ipj)
- if (ir_start.gt.1) ir_start=ir_start-1
- ir_end=chain_border(2,ipj)
- if (ir_end.lt.nres) ir_end=ir_end+1
- do k=ir_start,ir_end
- c(1,k)=c(1,k)+ixmin*boxxsize
- c(2,k)=c(2,k)+iymin*boxysize
- c(3,k)=c(3,k)+izmin*boxzsize
- c(1,k+nres)=c(1,k+nres)+ixmin*boxxsize
- c(2,k+nres)=c(2,k+nres)+iymin*boxysize
- c(3,k+nres)=c(3,k+nres)+izmin*boxzsize
+ 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
-c write (iout,*) "jmin",jmin," ipj",ipj,
-c & " ixmin",ixmin," iymin",iymin," izmin",izmin
- iaux=iper(i+1)
- iper(i+1)=iper(jmin)
- iper(jmin)=iaux
enddo
+ do j=1,3
+ ccm(j)=ccm(j)/nrestot
+ enddo
+#ifdef DEBUG
+ write (iout,'(a,3f10.5)') "ccm ",(ccm(j),j=1,3)
+ write (iout,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
+ write (iout,'(a,3f10.5)') "diff ",(ccm(j)-ccm_prev(j),j=1,3)
+ write (iout,*) "shift_coord"
+ do i=1,nchain
+ write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+ enddo
+#endif;
+#ifdef PREVIOUS
+c
+c Shift the system by box dimensions so that its CM be closest to the
+c CM of the pevious snapshot
+c
+ 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
+c Simply shift the CM of the whole oligomer to the origin of the
+c coordinate system
+ do ichain=1,nchain
+ do j=1,3
+ shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
+ enddo
+ enddo
+#endif
+#ifdef DEBUG
+ write (iout,*) "shift_coord"
+ do i=1,nchain
+ write (iout,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+ enddo
+#endif
+ do ichain=1,nchain
+ do i=chain_border1(1,ichain),chain_border1(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
+
return
end
projects / unres.git / blobdiff