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
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
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)
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
!-------------------------------------------------------------------------------------------------------
projects / unres4.git / blobdiff