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,1).ne.10) then
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)
endif
endif
do i=1,nres-1
+! 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)) 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,1).ne.ntyp1)).and.molnum(i).eq.1) then
write (iout,'(a,i4)') 'Bad Cartesians for residue',i
!test stop
endif
di=dist(i,nres+i)
!#ifndef WHAM_RUN
! 10/03/12 Adam: Correction for zero SC-SC bond length
+
if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
- di=dsc(itype(i,1))
+ di=dsc(itype(i,molnum(i)))
vbld(i+nres)=di
if (itype(i,1).ne.10) then
vbld_inv(i+nres)=1.0d0/di
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,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
integer :: i,j,ires,nscat
real(kind=8),dimension(3,20) :: sccor
real(kind=8) :: sccmj
- print *,"I am in sccenter",ires,nscat
+! 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
! 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,1).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
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,1).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,1).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
enddo
! The side-chain vector derivatives
do i=2,nres-1
- if(itype(i,1).ne.10 .and. itype(i,1).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)
! enddo
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,1).ne.ntyp1)) then
+ 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,1).ne.10).and.(itype(2,1).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,1).ne.10).and.(itype(3,1).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,1).ne.10).and.(itype(2,1).ne.ntyp1)) then
- if (itype(1,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,1).ne.10).and.(nres.ge.3).and.(itype(3,1).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,1).ne.10).and.(itype(1,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
do i=3,nres-2
do j=1,3
! write(iout,*) "before", gcart(j,i)
- if ((itype(i,1).ne.10).and.(itype(i,1).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,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,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,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,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,1).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,1).ne.10).and.(itype(nres-1,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,1).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,1).ne.10).and.(itype(nres,1).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,1).ne.10).and.(itype(nres-2,1).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,1).ne.10).and.(itype(nres,1).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),
endif
! Settind dE/ddnres
if ((nres.ge.3).and.(itype(nres,1).ne.10).and. &
- (itype(nres,1).ne.ntyp1))then
+ (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
#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
projects / unres4.git / blobdiff