*
* Version of March '95, based on an early version of November '91.
*
+* 03/11/20 Adam. Array fromto eliminated, computed on the fly
+* Fixed the problem with vbld indices, which caused errors in
+* derivatives when the backbone virtual bond lengths were not equal.
***********************************************************************
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.GEO'
include 'COMMON.LOCAL'
include 'COMMON.INTERACT'
- dimension drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),temp(3,3),
- & fromto(3,3,maxdim),prordt(3,3,maxres),prodrt(3,3,maxres)
- dimension xx(3),xx1(3)
- common /przechowalnia/ fromto
+ double precision drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),
+ &temp(3,3),prordt(3,3,maxres),prodrt(3,3,maxres)
+ double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp
+ double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,
+ & cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl
+ double precision fromto(3,3)
+ integer i,ii,j,jjj,k,l,m,indi,ind,ind1
* get the position of the jth ijth fragment of the chain coordinate system
* in the fromto array.
+ integer indmat
indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
+ call chainbuild_extconf
+ call cartprint
+ call intout
*
* calculate the derivatives of transformation matrix elements in theta
*
drt(3,3,i)=-rt(2,3,i)
enddo
*
-* generate the matrix products of type r(i)t(i)...r(j)t(j)
-*
- do i=2,nres-2
- ind=indmat(i,i+1)
- do k=1,3
- do l=1,3
- temp(k,l)=rt(k,l,i)
- enddo
- enddo
- do k=1,3
- do l=1,3
- fromto(k,l,ind)=temp(k,l)
- enddo
- enddo
- do j=i+1,nres-2
- ind=indmat(i,j+1)
- do k=1,3
- do l=1,3
- dpkl=0.0d0
- do m=1,3
- dpkl=dpkl+temp(k,m)*rt(m,l,j)
- enddo
- dp(k,l)=dpkl
- fromto(k,l,ind)=dpkl
- enddo
- enddo
- do k=1,3
- do l=1,3
- temp(k,l)=dp(k,l)
- enddo
- enddo
- enddo
- enddo
-*
* Calculate derivatives.
*
ind1=0
do i=1,nres-2
- ind1=ind1+1
+ ind1=ind1+1
*
* Derivatives of DC(i+1) in theta(i+2)
*
+c write (iout,*) "theta i",i
+c write(iout,'(7hprod 9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
+c write(iout,'(7hrdt 9f10.5)')((rdt(k,l,i),l=1,3),k=1,3)
+c write(iout,*) "vbld",vbld(i+2)
do j=1,3
do k=1,2
dpjk=0.0D0
prordt(j,k,i)=dp(j,k)
enddo
dp(j,3)=0.0D0
- dcdv(j,ind1)=vbld(i+1)*dp(j,1)
+c dcdv(j,ind1)=vbld(i+1)*dp(j,1)
+ dcdv(j,ind1)=vbld(i+2)*dp(j,1)
enddo
+c write(iout,'(7hdcdv 3f10.5)')(dcdv(k,ind1),k=1,3)
*
* Derivatives of SC(i+1) in theta(i+2)
*
enddo
dxdv(j,ind1)=rj
enddo
+c write (iout,*) "dxdv",(dxdv(j,ind1),j=1,3)
*
* Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
* than the other off-diagonal derivatives.
enddo
dxdv(j,ind1+1)=dxoiij
enddo
-cd print *,ind1+1,(dxdv(j,ind1+1),j=1,3)
+c write(iout,*)ind1+1,(dxdv(j,ind1+1),j=1,3)
*
* Derivatives of DC(i+1) in phi(i+2)
*
dp(j,k)=dpjk
prodrt(j,k,i)=dp(j,k)
enddo
- dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
+c dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
+ dcdv(j+3,ind1)=vbld(i+2)*dp(j,1)
enddo
*
* Derivatives of SC(i+1) in phi(i+2)
* theta(nres) and phi(i+3) thru phi(nres).
*
do j=i+1,nres-2
- ind1=ind1+1
- ind=indmat(i+1,j+1)
-cd print *,'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
+ ind1=ind1+1
+ ind=indmat(i+1,j+1)
+c write(iout,*)'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
+ call build_fromto(i+1,j+1,fromto)
+c write(iout,'(7hfromto 9f10.5)')((fromto(k,l),l=1,3),k=1,3)
do k=1,3
do l=1,3
tempkl=0.0D0
do m=1,2
- tempkl=tempkl+prordt(k,m,i)*fromto(m,l,ind)
+ tempkl=tempkl+prordt(k,m,i)*fromto(m,l)
enddo
temp(k,l)=tempkl
enddo
enddo
-cd print '(9f8.3)',((fromto(k,l,ind),l=1,3),k=1,3)
-cd print '(9f8.3)',((prod(k,l,i),l=1,3),k=1,3)
-cd print '(9f8.3)',((temp(k,l),l=1,3),k=1,3)
+c write(iout,'(7hfromto 9f10.5)')((fromto(k,l,ind),l=1,3),k=1,3)
+c write(iout,'(7hprod 9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
+c write(iout,'(7htemp 9f10.5)')((temp(k,l),l=1,3),k=1,3)
* Derivatives of virtual-bond vectors in theta
do k=1,3
- dcdv(k,ind1)=vbld(i+1)*temp(k,1)
+c dcdv(k,ind1)=vbld(i+1)*temp(k,1)
+ dcdv(k,ind1)=vbld(j+2)*temp(k,1)
enddo
-cd print '(3f8.3)',(dcdv(k,ind1),k=1,3)
+c write(iout,'(7hdcdv 3f10.5)')(dcdv(k,ind1),k=1,3)
* Derivatives of SC vectors in theta
do k=1,3
dxoijk=0.0D0
enddo
dxdv(k,ind1+1)=dxoijk
enddo
+c write(iout,'(7htheta 3f10.5)')(dxdv(k,ind1),k=1,3)
*
*--- Calculate the derivatives in phi
*
+#ifdef FIVEDIAG
+ do k=1,3
+ do l=1,3
+ tempkl=0.0D0
+ do m=1,3
+ tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
+ enddo
+ temp(k,l)=tempkl
+ enddo
+ enddo
+#else
do k=1,3
do l=1,3
tempkl=0.0D0
temp(k,l)=tempkl
enddo
enddo
+#endif
do k=1,3
- dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
- enddo
+c dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
+ dcdv(k+3,ind1)=vbld(j+2)*temp(k,1)
+ enddo
do k=1,3
dxoijk=0.0D0
do l=1,3
enddo
enddo
enddo
+#ifdef DEBUG
+ write (iout,*)
+ write (iout,'(a)') '****************** ddc/dtheta'
+ write (iout,*)
+ do i=1,nres-2
+ do j=i+1,nres-1
+ ii = indmat(i,j)
+ write (iout,'(2i4,3e14.6)') i,j,(dcdv(k,ii),k=1,3)
+ enddo
+ enddo
+ write (iout,*)
+ write (iout,'(a)') '******************* ddc/dphi'
+ write (iout,*)
+ do i=1,nres-3
+ do j=i+2,nres-1
+ ii = indmat(i+1,j)
+ write (iout,'(2i4,3e14.6)') i,j,(dcdv(k+3,ii),k=1,3)
+ write (iout,'(a)')
+ enddo
+ enddo
+ write (iout,'(a)')
+ write (iout,'(a)') '**************** dx/dtheta'
+ write (iout,'(a)')
+ do i=3,nres
+ do j=i-1,nres-1
+ ii = indmat(i-2,j)
+ write (iout,'(2i4,3e14.6)') i,j,(dxdv(k,ii),k=1,3)
+ enddo
+ enddo
+ write (iout,'(a)')
+ write (iout,'(a)') '***************** dx/dphi'
+ write (iout,'(a)')
+ do i=4,nres
+ do j=i-1,nres-1
+ ii = indmat(i-2,j)
+ write (iout,'(2i4,3e14.6)') i,j,(dxdv(k+3,ii),k=1,3)
+ write(iout,'(a)')
+ enddo
+ enddo
+#endif
*
* Derivatives in alpha and omega:
*
if(alphi.ne.alphi) alphi=100.0
if(omegi.ne.omegi) omegi=-100.0
#else
- alphi=alph(i)
- omegi=omeg(i)
+ alphi=alph(i)
+ omegi=omeg(i)
#endif
cd print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
- cosalphi=dcos(alphi)
- sinalphi=dsin(alphi)
- cosomegi=dcos(omegi)
- sinomegi=dsin(omegi)
- temp(1,1)=-dsci*sinalphi
- temp(2,1)= dsci*cosalphi*cosomegi
- temp(3,1)=-dsci*cosalphi*sinomegi
- temp(1,2)=0.0D0
- temp(2,2)=-dsci*sinalphi*sinomegi
- temp(3,2)=-dsci*sinalphi*cosomegi
- theta2=pi-0.5D0*theta(i+1)
- cost2=dcos(theta2)
- sint2=dsin(theta2)
- jjj=0
+ cosalphi=dcos(alphi)
+ sinalphi=dsin(alphi)
+ cosomegi=dcos(omegi)
+ sinomegi=dsin(omegi)
+ temp(1,1)=-dsci*sinalphi
+ temp(2,1)= dsci*cosalphi*cosomegi
+ temp(3,1)=-dsci*cosalphi*sinomegi
+ temp(1,2)=0.0D0
+ temp(2,2)=-dsci*sinalphi*sinomegi
+ temp(3,2)=-dsci*sinalphi*cosomegi
+ theta2=pi-0.5D0*theta(i+1)
+ cost2=dcos(theta2)
+ sint2=dsin(theta2)
+ jjj=0
cd print *,((temp(l,k),l=1,3),k=1,2)
do j=1,2
- xp=temp(1,j)
- yp=temp(2,j)
- xxp= xp*cost2+yp*sint2
- yyp=-xp*sint2+yp*cost2
- zzp=temp(3,j)
- xx(1)=xxp
- xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
- xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
- do k=1,3
- dj=0.0D0
- do l=1,3
- dj=dj+prod(k,l,i-1)*xx(l)
+ xp=temp(1,j)
+ yp=temp(2,j)
+ xxp= xp*cost2+yp*sint2
+ yyp=-xp*sint2+yp*cost2
+ zzp=temp(3,j)
+ xx(1)=xxp
+ xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
+ xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
+ do k=1,3
+ dj=0.0D0
+ do l=1,3
+ dj=dj+prod(k,l,i-1)*xx(l)
enddo
- dxds(jjj+k,i)=dj
+ dxds(jjj+k,i)=dj
enddo
- jjj=jjj+3
- enddo
+ jjj=jjj+3
+ enddo
enddo
return
end
-
projects / unres.git / blobdiff