projects
/
unres4.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
working gradient for lipid and shielding
[unres4.git]
/
source
/
unres
/
geometry.f90
diff --git
a/source/unres/geometry.f90
b/source/unres/geometry.f90
index
780684e
..
951e8eb
100644
(file)
--- a/
source/unres/geometry.f90
+++ b/
source/unres/geometry.f90
@@
-9,6
+9,9
@@
use energy_data
implicit none
!-----------------------------------------------------------------------------
use energy_data
implicit none
!-----------------------------------------------------------------------------
+! commom.bounds
+! common /bounds/
+!-----------------------------------------------------------------------------
! commom.chain
! common /chain/
! common /rotmat/
! commom.chain
! common /chain/
! common /rotmat/
@@
-80,17
+83,21
@@
nres2=2*nres
! Set lprn=.true. for debugging
lprn = .false.
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).
!
!
! 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
call orig_frame
+!elwrite(iout,*)"after orig_frame"
!
! Build the alpha-carbon chain.
!
do i=4,nres
call locate_next_res(i)
enddo
!
! 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.
!
!
! First and last SC must coincide with the corresponding CA.
!
@@
-111,15
+118,15
@@
write (iout,'(/a)') 'Recalculated internal coordinates'
do i=2,nres-1
do j=1,3
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)
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=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
+ alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
enddo
1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
enddo
1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
@@
-407,6
+414,7
@@
#ifdef WHAM_RUN
vbld(nres+1)=0.0d0
#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
vbld(2*nres)=0.0d0
vbld_inv(nres+1)=0.0d0
vbld_inv(2*nres)=0.0d0
@@
-420,7
+428,7
@@
dnorm1=dist(i-1,i)
dnorm2=dist(i,i+1)
do j=1,3
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
+(c(j,i+1)-c(j,i))/dnorm2)
enddo
be=0.0D0
@@
-438,8
+446,8
@@
tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
endif
endif
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)
vbld_inv(i)=1.0d0/vbld(i)
theta(i+1)=alpha(i-1,i,i+1)
vbld(i)=dist(i-1,i)
vbld_inv(i)=1.0d0/vbld(i)
@@
-503,11
+511,19
@@
#endif
do i=1,nres-1
do j=1,3
#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
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
dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
enddo
enddo
@@
-524,7
+540,7
@@
#endif
return
end subroutine int_from_cart1
#endif
return
end subroutine int_from_cart1
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! check_sc_distr.f
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! check_sc_distr.f
!-----------------------------------------------------------------------------
@@
-642,14
+658,18
@@
ii=ialph(i,2)
alph(ii)=x(nphi+ntheta+i)
omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
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)
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)
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
if (theta(i).eq.pi) theta(i)=0.99d0*pi
x(i-2+nphi)=theta(i)
enddo
@@
-751,7
+771,7
@@
thetnorm=xx
return
end function thetnorm
thetnorm=xx
return
end function thetnorm
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
subroutine var_to_geom_restr(n,xx)
!
!-----------------------------------------------------------------------------
subroutine var_to_geom_restr(n,xx)
!
@@
-953,13
+973,14
@@
! SCs.
iteli=itel(i)
do j=1,3
! 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))
!d print *,'overlap, p-Sc: i=',i,' j=',j,
!d & ' dist=',dist(nres+j,maxres2+1)
enddo
do j=nnt,i-2
itj=iabs(itype(j))
!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
overlap=.true.
return
endif
@@
-968,28
+989,28
@@
! groups.
do j=1,nnt-2
do k=1,3
! 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)
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
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
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
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
overlap=.true.
return
endif
@@
-1321,7
+1342,7
@@
endif
if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
#ifdef MPI
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.'
write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
#else
! write (iout,'(a)') 'Bad sampling box.'
@@
-1807,7
+1828,7
@@
dist=dsqrt(x12*x12+y12*y12+z12*z12)
return
end function dist
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
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! local_move.f
!-----------------------------------------------------------------------------
@@
-2830,16
+2851,33
@@
endif
endif
do i=1,nres-1
endif
endif
do i=1,nres-1
+!in wham do i=1,nres
iti=itype(i)
iti=itype(i)
- if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+ if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
+ (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1)) then
write (iout,'(a,i4)') 'Bad Cartesians for residue',i
!test stop
endif
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)
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
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).eq.ntyp1) then
+! do j=1,3
+! c(j,1)=c(j,2)+(c(j,3)-c(j,4))
+! enddo
+! endif
+! if (itype(nres).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
! theta(3)=90.0d0*deg2rad
! if (unres_pdb) then
! if (itype(1).eq.21) then
! theta(3)=90.0d0*deg2rad
@@
-2857,11
+2895,13
@@
if (lside) then
do i=2,nres-1
do j=1,3
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))
+(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)
di=dist(i,nres+i)
enddo
iti=itype(i)
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))
! 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))
@@
-2871,9
+2911,10
@@
else
vbld_inv(i+nres)=0.0d0
endif
else
vbld_inv(i+nres)=0.0d0
endif
+!#endif
if (iti.ne.10) then
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(me.eq.king.or..not.out1file)then
if (lprn) &
endif
if(me.eq.king.or..not.out1file)then
if (lprn) &
@@
-2940,7
+2981,9
@@
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
it=itype(i)
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
it=itype(i)
- if (it.ne.10) then
+
+ 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
!
! Compute the axes of tghe local cartesian coordinates system; store in
! x_prime, y_prime and z_prime
@@
-2985,6
+3028,7
@@
yyref(i),zzref(i)
enddo
endif
yyref(i),zzref(i)
enddo
endif
+
return
end subroutine sc_loc_geom
!-----------------------------------------------------------------------------
return
end subroutine sc_loc_geom
!-----------------------------------------------------------------------------
@@
-3004,7
+3048,7
@@
enddo
return
end subroutine sccenter
enddo
return
end subroutine sccenter
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
subroutine bond_regular
use calc_data
!-----------------------------------------------------------------------------
subroutine bond_regular
use calc_data
@@
-3375,7
+3419,7
@@
! The side-chain vector derivatives
return
end subroutine int_to_cart
! The side-chain vector derivatives
return
end subroutine int_to_cart
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! readrtns_CSA.F
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! readrtns_CSA.F
!-----------------------------------------------------------------------------
@@
-3431,7
+3475,7
@@
!d & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
!d & dhpb(i),forcon(i)
!d enddo
!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
return
end subroutine gen_dist_constr
@@
-3479,18
+3523,22
@@
! real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2)
! real(kind=8),dimension(:,:),allocatable :: dc
allocate(dc_old(3,0:nres2))
! 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"
!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"
!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/
!elwrite(iout,*) "jestem w alloc geo 1"
allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
! common /rotmat/
@@
-3533,8
+3581,17
@@
if(.not.allocated(yyref)) allocate(yyref(nres))
if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
allocate(ialph(nres,2)) !(maxres,2)
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)
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
!-----------------------------------------------------------------------------
return
end subroutine alloc_geo_arrays
!-----------------------------------------------------------------------------