- module geometry
+ module geometry
!-----------------------------------------------------------------------------
use io_units
use names
use energy_data
implicit none
!-----------------------------------------------------------------------------
+! commom.bounds
+! common /bounds/
+!-----------------------------------------------------------------------------
! commom.chain
! common /chain/
! common /rotmat/
nres2=2*nres
! Set lprn=.true. for debugging
lprn = .false.
+ print *,"I ENTER CHAINBUILD"
!
! Define the origin and orientation of the coordinate system and locate the
! first three CA's and SC(2).
!
+!elwrite(iout,*)"in chainbuild"
call orig_frame
+!elwrite(iout,*)"after orig_frame"
!
! Build the alpha-carbon chain.
!
do i=4,nres
call locate_next_res(i)
enddo
+!elwrite(iout,*)"after locate_next_res"
!
! First and last SC must coincide with the corresponding CA.
!
write (iout,'(/a)') 'Recalculated internal coordinates'
do i=2,nres-1
do j=1,3
- c(j,nres2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres
+ c(j,nres2+2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres
enddo
be=0.0D0
if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
- be1=rad2deg*beta(nres+i,i,nres2,i+1)
+ be1=rad2deg*beta(nres+i,i,nres2+2,i+1)
alfai=0.0D0
if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
- write (iout,1212) restyp(itype(i)),i,dist(i-1,i),&
- alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2),be1
+ 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))
real(kind=8) :: alphi,omegi,theta2
real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
real(kind=8) :: xp,yp,zp,cost2,sint2,rj
-! dsci=dsc(itype(i))
-! dsci_inv=dsc_inv(itype(i))
+! dsci=dsc(itype(i,1))
+! dsci_inv=dsc_inv(itype(i,1))
dsci=vbld(i+nres)
dsci_inv=vbld_inv(i+nres)
#ifdef OSF
xx(1)= xp*cost2+yp*sint2
xx(2)=-xp*sint2+yp*cost2
xx(3)= zp
-!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
+!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i,1)),i,
!d & xp,yp,zp,(xx(k),k=1,3)
do j=1,3
xloc(j,i)=xx(j)
#ifdef WHAM_RUN
vbld(nres+1)=0.0d0
+!write(iout,*)"geometry warring, vbld=",(vbld(i),i=1,nres+1)
vbld(2*nres)=0.0d0
vbld_inv(nres+1)=0.0d0
vbld_inv(2*nres)=0.0d0
dnorm1=dist(i-1,i)
dnorm2=dist(i,i+1)
do j=1,3
- c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 &
+ c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 &
+(c(j,i+1)-c(j,i))/dnorm2)
enddo
be=0.0D0
if (i.gt.2) then
if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
- if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
+ if ((itype(i,1).ne.10).and.(itype(i-1,1).ne.10)) then
tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
endif
- if (itype(i-1).ne.10) then
+ if (itype(i-1,1).ne.10) then
tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
omicron(1,i)=alpha(i-2,i-1,i-1+nres)
omicron(2,i)=alpha(i-1+nres,i-1,i)
endif
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
endif
endif
- omeg(i)=beta(nres+i,i,nres2,i+1)
- alph(i)=alpha(nres+i,i,nres2)
+ omeg(i)=beta(nres+i,i,nres2+2,i+1)
+ alph(i)=alpha(nres+i,i,nres2+2)
theta(i+1)=alpha(i-1,i,i+1)
vbld(i)=dist(i-1,i)
+! print *,i,vbld(i),"vbld(i)"
vbld_inv(i)=1.0d0/vbld(i)
vbld(nres+i)=dist(nres+i,i)
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
vbld_inv(nres+i)=1.0d0/vbld(nres+i)
else
vbld_inv(nres+i)=0.0d0
#endif
do i=1,nres-1
do j=1,3
+!#ifdef WHAM_RUN
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ dc(j,i)=c(j,i+1)-c(j,i)
+#endif
dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
enddo
enddo
do i=2,nres-1
do j=1,3
+!#ifdef WHAM_RUN
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+#endif
dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
enddo
enddo
if (lprn) then
do i=2,nres
- write (iout,1212) restyp(itype(i)),i,vbld(i),&
+ write (iout,1212) restyp(itype(i,1),1),i,vbld(i),&
rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),&
rad2deg*alph(i),rad2deg*omeg(i)
enddo
#endif
return
end subroutine int_from_cart1
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! check_sc_distr.f
!-----------------------------------------------------------------------------
print *,'dv=',dv
do 10 it=1,1
if (it.eq.10) goto 10
- open (20,file=restyp(it)//'_distr.sdc',status='unknown')
+ open (20,file=restyp(it,1)//'_distr.sdc',status='unknown')
call gen_side(it,90.0D0 * deg2rad,al,om,fail)
close (20)
goto 10
- open (20,file=restyp(it)//'_distr1.sdc',status='unknown')
+ open (20,file=restyp(it,1)//'_distr1.sdc',status='unknown')
do i=0,90
do j=0,72
prob(j,i)=0.0D0
ii=ialph(i,2)
alph(ii)=x(nphi+ntheta+i)
omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
+!elwrite(iout,*) "alph",ii,alph
+!elwrite(iout,*) "omeg",ii,omeg
enddo
endif
do i=4,nres
phi(i)=x(i-3)
+!elwrite(iout,*) "phi",i,phi
enddo
if (n.eq.nphi) return
do i=3,nres
theta(i)=x(i-2+nphi)
+!elwrite(iout,*) "theta",i,theta
if (theta(i).eq.pi) theta(i)=0.99d0*pi
x(i-2+nphi)=theta(i)
enddo
thetnorm=xx
return
end function thetnorm
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
subroutine var_to_geom_restr(n,xx)
!
maxsi=100
!d write (iout,*) 'Gen_Rand_conf: nstart=',nstart
if (nstart.lt.5) then
- it1=iabs(itype(2))
- phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
+ it1=iabs(itype(2,1))
+ phi(4)=gen_phi(4,iabs(itype(2,1)),iabs(itype(3,1)))
! write(iout,*)'phi(4)=',rad2deg*phi(4)
- if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
+ if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2,1)),pi,phi(4))
! write(iout,*)'theta(3)=',rad2deg*theta(3)
if (it1.ne.10) then
nsi=0
endif
return 1
endif
- it1=iabs(itype(i-1))
- it2=iabs(itype(i-2))
- it=iabs(itype(i))
+ 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)
nres2=2*nres
data redfac /0.5D0/
overlap=.false.
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
if (iti.gt.ntyp) return
! Check for SC-SC overlaps.
!d print *,'nnt=',nnt,' nct=',nct
do j=nnt,i-1
- itj=iabs(itype(j))
+ itj=iabs(itype(j,1))
if (j.lt.i-1 .or. ipot.ne.4) then
rcomp=sigmaii(iti,itj)
else
! SCs.
iteli=itel(i)
do j=1,3
- c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+! c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+ c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
enddo
do j=nnt,i-2
- itj=iabs(itype(j))
+ itj=iabs(itype(j,1))
!d print *,'overlap, p-Sc: i=',i,' j=',j,
!d & ' dist=',dist(nres+j,maxres2+1)
- if (dist(nres+j,nres2+1).lt.4.0D0*redfac) then
+ if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
overlap=.true.
return
endif
! groups.
do j=1,nnt-2
do k=1,3
- c(k,nres2+1)=0.5D0*(c(k,j)+c(k,j+1))
+ c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1))
enddo
!d print *,'overlap, SC-p: i=',i,' j=',j,
!d & ' dist=',dist(nres+i,maxres2+1)
- if (dist(nres+i,nres2+1).lt.4.0D0*redfac) then
+ if (dist(nres+i,nres2+3).lt.4.0D0*redfac) then
overlap=.true.
return
endif
enddo
! Check for p-p overlaps
do j=1,3
- c(j,nres2+2)=0.5D0*(c(j,i)+c(j,i+1))
+ c(j,nres2+4)=0.5D0*(c(j,i)+c(j,i+1))
enddo
do j=nnt,i-2
itelj=itel(j)
do k=1,3
- c(k,nres2+2)=0.5D0*(c(k,j)+c(k,j+1))
+ c(k,nres2+4)=0.5D0*(c(k,j)+c(k,j+1))
enddo
!d print *,'overlap, p-p: i=',i,' j=',j,
!d & ' dist=',dist(maxres2+1,maxres2+2)
if(iteli.ne.0.and.itelj.ne.0)then
- if (dist(nres2+1,nres2+2).lt.rpp(iteli,itelj)*redfac) then
+ if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then
overlap=.true.
return
endif
endif
if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
#ifdef MPI
- write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
+ write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax
write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
#else
! write (iout,'(a)') 'Bad sampling box.'
do ires=1,ioverlap_last
i=ioverlap(ires)
- iti=iabs(itype(i))
+ iti=iabs(itype(i,1))
if (iti.ne.10) then
nsi=0
fail=.true.
! print *,'>>overlap_sc nnt=',nnt,' nct=',nct
ind=0
do i=iatsc_s,iatsc_e
- itypi=iabs(itype(i))
- itypi1=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)
zi=c(3,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))
+ itypj=iabs(itype(j,molnum(j)))
dscj_inv=dsc_inv(itypj)
sig0ij=sigma(itypi,itypj)
chi1=chi(itypi,itypj)
end subroutine sc_angular
!-----------------------------------------------------------------------------
! initialize_p.F
+ subroutine sc_angular_nucl
+! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
+! om12. Called by ebp, egb, and egbv.
+! use calc_data
+! implicit none
+! include 'COMMON.CALC'
+! include 'COMMON.IOUNITS'
+ use comm_locel
+ use calc_data_nucl
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ chiom12=chi12*om12
+! Calculate eps1(om12) and its derivative in om12
+ faceps1=1.0D0-om12*chiom12
+ faceps1_inv=1.0D0/faceps1
+ eps1=dsqrt(faceps1_inv)
+! Following variable is eps1*deps1/dom12
+ eps1_om12=faceps1_inv*chiom12
+! diagnostics only
+! faceps1_inv=om12
+! eps1=om12
+! eps1_om12=1.0d0
+! write (iout,*) "om12",om12," eps1",eps1
+! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
+! and om12.
+ om1om2=om1*om2
+ chiom1=chi1*om1
+ chiom2=chi2*om2
+ facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
+ sigsq=1.0D0-facsig*faceps1_inv
+ sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
+ sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
+ sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
+ chipom1=chip1*om1
+ chipom2=chip2*om2
+ chipom12=chip12*om12
+ facp=1.0D0-om12*chipom12
+ facp_inv=1.0D0/facp
+ facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
+! write (iout,*) "chipom1",chipom1," chipom2",chipom2,
+! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv
+! Following variable is the square root of eps2
+ eps2rt=1.0D0-facp1*facp_inv
+! Following three variables are the derivatives of the square root of eps
+! in om1, om2, and om12.
+ eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
+ eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
+ eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
+! Evaluate the "asymmetric" factor in the VDW constant, eps3
+ eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
+! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
+! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
+! & " eps2rt_om12",eps2rt_om12
+! Calculate whole angle-dependent part of epsilon and contributions
+! to its derivatives
+ return
+ end subroutine sc_angular_nucl
+
!-----------------------------------------------------------------------------
subroutine int_bounds(total_ints,lower_bound,upper_bound)
! implicit real*8 (a-h,o-z)
dist=dsqrt(x12*x12+y12*y12+z12*z12)
return
end function dist
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! local_move.f
!-----------------------------------------------------------------------------
endif
endif
do i=1,nres-1
- iti=itype(i)
- if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+! if (molnum(i).ne.1) cycle
+!in wham do i=1,nres
+ iti=itype(i,1)
+ if (((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
+ (iti.ne.ntyp1 .and. itype(i+1,1).ne.ntyp1)).and.molnum(i).eq.1) then
write (iout,'(a,i4)') 'Bad Cartesians for residue',i
!test stop
endif
+!#ifndef WHAM_RUN
vbld(i+1)=dist(i,i+1)
vbld_inv(i+1)=1.0d0/vbld(i+1)
+!#endif
if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
enddo
+!el -----
+!#ifdef WHAM_RUN
+! if (itype(1,1).eq.ntyp1) then
+! do j=1,3
+! c(j,1)=c(j,2)+(c(j,3)-c(j,4))
+! enddo
+! endif
+! if (itype(nres,1).eq.ntyp1) then
+! do j=1,3
+! c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
+! enddo
+! endif
+!#endif
! if (unres_pdb) then
-! if (itype(1).eq.21) then
+! if (itype(1,1).eq.21) then
! theta(3)=90.0d0*deg2rad
! phi(4)=180.0d0*deg2rad
! vbld(2)=3.8d0
! vbld_inv(2)=1.0d0/vbld(2)
! endif
-! if (itype(nres).eq.21) then
+! if (itype(nres,1).eq.21) then
! theta(nres)=90.0d0*deg2rad
! phi(nres)=180.0d0*deg2rad
! vbld(nres)=3.8d0
if (lside) then
do i=2,nres-1
do j=1,3
- c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
+ c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
+(c(j,i+1)-c(j,i))*vbld_inv(i+1))
+! in wham c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)
enddo
- iti=itype(i)
+ iti=itype(i,1)
di=dist(i,nres+i)
+!#ifndef WHAM_RUN
! 10/03/12 Adam: Correction for zero SC-SC bond length
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1 .and. di.eq.0.0d0) &
- di=dsc(itype(i))
+
+ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
+ di=dsc(itype(i,molnum(i)))
vbld(i+nres)=di
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
vbld_inv(i+nres)=1.0d0/di
else
vbld_inv(i+nres)=0.0d0
endif
+!#endif
if (iti.ne.10) then
- alph(i)=alpha(nres+i,i,nres2)
- omeg(i)=beta(nres+i,i,nres2,i+1)
+ alph(i)=alpha(nres+i,i,nres2+2)
+ omeg(i)=beta(nres+i,i,nres2+2,i+1)
endif
+ if (iti.ne.0) then
if(me.eq.king.or..not.out1file)then
if (lprn) &
- write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),&
+ write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
rad2deg*alph(i),rad2deg*omeg(i)
endif
+ else
+ if(me.eq.king.or..not.out1file)then
+ if (lprn) &
+ write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
+ rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),&
+ rad2deg*alph(i),rad2deg*omeg(i)
+ endif
+ endif
enddo
else if (lprn) then
do i=2,nres
- iti=itype(i)
+ iti=itype(i,1)
if(me.eq.king.or..not.out1file) &
- write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),&
+ write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,dist(i,i-1),&
rad2deg*theta(i),rad2deg*phi(i)
enddo
endif
enddo
enddo
do i=2,nres-1
- if (itype(i).ne.10) then
+ if (itype(i,1).ne.10) then
do j=1,3
dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
enddo
cosfac=dsqrt(cosfac2)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
- it=itype(i)
- if (it.ne.10) then
+ it=itype(i,1)
+
+ if ((it.ne.10).and.(it.ne.ntyp1)) then
+!el if (it.ne.10) then
!
! Compute the axes of tghe local cartesian coordinates system; store in
! x_prime, y_prime and z_prime
enddo
if (lprn) then
do i=2,nres
- iti=itype(i)
+ iti=itype(i,1)
if(me.eq.king.or..not.out1file) &
- write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),&
+ write (iout,'(a3,i4,3f10.5)') restyp(iti,1),i,xxref(i),&
yyref(i),zzref(i)
enddo
endif
+
return
end subroutine sc_loc_geom
!-----------------------------------------------------------------------------
integer :: i,j,ires,nscat
real(kind=8),dimension(3,20) :: sccor
real(kind=8) :: sccmj
+! print *,"I am in sccenter",ires,nscat
do j=1,3
sccmj=0.0D0
do i=1,nscat
- sccmj=sccmj+sccor(j,i)
+ sccmj=sccmj+sccor(j,i)
+!C print *,"insccent", ires,sccor(j,i)
enddo
dc(j,ires)=sccmj/nscat
enddo
return
end subroutine sccenter
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
subroutine bond_regular
use calc_data
! 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))
- vbld_inv(i+1+nres)=dsc_inv(itype(i+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
! calculating dE/ddc1
!el local variables
integer :: j,i
+! print *,"gloc",gloc(:,:)
+! print *, "gcart",gcart(:,:)
if (nres.lt.3) go to 18
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)
- if(itype(2).ne.10) then
+ if ((itype(2,1).ne.10).and.&
+ (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) 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).ne.10) then
+ if(itype(2,1).ne.10) 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
- if(itype(3).ne.10) then
+ if(itype(3,1).ne.10) then
gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
endif
gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
dtheta(j,1,5)
- if(itype(3).ne.10) then
+! if(itype(3,1).ne.10) then
+ if ((itype(3,1).ne.10).and.&
+ (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) then
gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
endif
- if(itype(4).ne.10) then
+! if(itype(4,1).ne.10) then
+ if ((itype(4,1).ne.10).and.&
+ (itype(4,molnum(4)).ne.ntyp1_molec(molnum(4)))) then
gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
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).ne.10) then
+ if(itype(i,1).ne.10) 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
- if(itype(i+1).ne.10) then
+ if(itype(i+1,1).ne.10) then
gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
+gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
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).ne.10) then
+ if(itype(nres-2,1).ne.10) 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).ne.10) then
+ if(itype(nres-1,1).ne.10) 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)
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)
- if(itype(nres-1).ne.10) then
+ if(itype(nres-1,1).ne.10) 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)
enddo
! The side-chain vector derivatives
do i=2,nres-1
- if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if(itype(i,1).ne.10 .and. &
+ itype(i,molnum(i)).ne.ntyp1_molec(molnum(i))) 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)
! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
! enddo
if (nres.lt.2) return
- if ((nres.lt.3).and.(itype(1).eq.10)) return
- if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+ 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
do j=1,3
!c Derviative was calculated for oposite vector of side chain therefore
! there is "-" sign before gloc_sc
dtauangle(j,1,1,3)
gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
dtauangle(j,1,2,3)
- if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+ if ((itype(2,1).ne.10).and. &
+ (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) 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)* &
endif
enddo
endif
- if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1)) &
+ if ((nres.ge.3).and.(itype(3,molnum(3)).ne.10).and.&
+ (itype(3,molnum(3)).ne.ntyp1_molec(molnum(3)))) &
then
do j=1,3
gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
! Calculating the remainder of dE/ddc2
do j=1,3
- if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
- if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
+ if((itype(2,1).ne.10).and. &
+ (itype(2,molnum(2)).ne.ntyp1_molec(molnum(2)))) then
+ if ((itype(1,1).ne.10).and.&
+ ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))))&
+ gxcart(j,2)=gxcart(j,2)+ &
gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
- if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1)) &
+ if ((itype(3,1).ne.10).and.(nres.ge.3).and.(itype(3,molnum(3)).ne.ntyp1_molec(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
gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
endif
if (nres.gt.3) then
+! if ((itype(1,1).ne.10).and.&
+! ((itype(1,molnum(1)).ne.ntyp1_molec(molnum(1))))) &
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,1,4),"gx"
endif
endif
- if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+ if ((itype(1,1).ne.10).and.&
+ (itype(1,molnum(1)).ne.ntyp1_molec(molnum(1)))) 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).ne.10).and.(nres.ge.3)) then
+ if ((itype(3,1).ne.10).and.(nres.ge.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).ne.10).and.(nres.ge.4)) then
+ if ((itype(4,1).ne.10).and.(nres.ge.4)) 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),itype(1),itype(3), "gcart2"
+! write(iout,*) gcart(j,2),itype(2,1),itype(1,1),itype(3,1), "gcart2"
enddo
! 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).ne.10).and.(itype(i).ne.ntyp1)) then
+ if ((itype(i,1).ne.10).and.&
+ (itype(i,molnum(i)).ne.ntyp1_molec(molnum(i)))) 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)
*dtauangle(j,1,2,i+2)
! write(iout,*) "new",j,i,
! & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
- if (itype(i-1).ne.10) then
+! 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
+
gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
*dtauangle(j,3,3,i+1)
endif
- if (itype(i+1).ne.10) then
- gxcart(j,i)=gxcart(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)))) 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) &
*dtauangle(j,3,2,i+2)
endif
endif
- if (itype(i-1).ne.10) then
+! 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
gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
dtauangle(j,1,3,i+1)
endif
- if (itype(i+1).ne.10) then
+! 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
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),
! & dtauangle(j,2,2,i+2)
endif
- if (itype(i+2).ne.10) then
+! 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
gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
dtauangle(j,2,1,i+3)
endif
! Setting dE/ddnres-1
if(nres.ge.4) then
do j=1,3
- if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
+ if ((itype(nres-1,1).ne.10).and.&
+ (itype(nres-1,molnum(nres-1)).ne.ntyp1_molec(molnum(nres-1)))) 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).ne.10) then
- gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
+! 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
+ gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
*dtauangle(j,3,3,nres)
endif
- if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+ if ((itype(nres,1).ne.10).and.&
+ (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) 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) &
*dtauangle(j,3,2,nres+1)
endif
endif
- if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
+ if ((itype(nres-2,1).ne.10).and.&
+ (itype(nres-2,molnum(nres-2)).ne.ntyp1_molec(molnum(nres-2)))) then
gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
dtauangle(j,1,3,nres)
endif
- if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+ if ((itype(nres,1).ne.10).and.(itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres)))) 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),
-! & dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres)
+! & dtauangle(j,2,2,nres+1), itype(nres-1,1),itype(nres,1)
endif
enddo
endif
! Settind dE/ddnres
- if ((nres.ge.3).and.(itype(nres).ne.10).and. &
- (itype(nres).ne.ntyp1))then
+ if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
+ (itype(nres,molnum(nres)).ne.ntyp1_molec(molnum(nres))))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) &
enddo
endif
! The side-chain vector derivatives
+! print *,"gcart",gcart(:,:)
return
end subroutine int_to_cart
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! readrtns_CSA.F
!-----------------------------------------------------------------------------
!d & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
!d & dhpb(i),forcon(i)
!d enddo
- deallocate(itype_pdb)
+! deallocate(itype_pdb)
return
end subroutine gen_dist_constr
write (iout,100)
do i=1,nres
- write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+ write (iout,110) restyp(itype(i,1),1),i,c(1,i),c(2,i),&
c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
enddo
100 format (//' alpha-carbon coordinates ',&
! real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2)
! real(kind=8),dimension(:,:),allocatable :: dc
allocate(dc_old(3,0:nres2))
- if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2)) !(3,0:maxres2)
+! if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
+ if(.not.allocated(dc_norm2)) then
+ allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
+ dc_norm2(:,:)=0.d0
+ endif
+!
!el if(.not.allocated(dc_norm))
!elwrite(iout,*) "jestem w alloc geo 1"
- if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:nres2)) !(3,0:maxres2)
+ if(.not.allocated(dc_norm)) then
+ allocate(dc_norm(3,0:nres2+2)) !(3,0:maxres2)
+ dc_norm(:,:)=0.d0
+ endif
!elwrite(iout,*) "jestem w alloc geo 1"
allocate(xloc(3,nres),xrot(3,nres))
!elwrite(iout,*) "jestem w alloc geo 1"
- do i=1,nres
- do j=1,3
- xloc(j,i)=0.0D0
- enddo
- enddo
+ xloc(:,:)=0.0D0
!elwrite(iout,*) "jestem w alloc geo 1"
allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
! common /rotmat/
if(.not.allocated(yyref)) allocate(yyref(nres))
if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
allocate(ialph(nres,2)) !(maxres,2)
+ ialph(:,1)=0
+ ialph(:,2)=0
allocate(ivar(4*nres2)) !(4*maxres2)
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ allocate(vbld(2*nres))
+ vbld(:)=0.d0
+ allocate(vbld_inv(2*nres))
+ vbld_inv(:)=0.d0
+#endif
+
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
projects / unres4.git / blobdiff