X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc-HCD-5D%2Fcartder.F;h=36b4e6383e9e70317bd5bafd0220fe553c3387cb;hb=75c81b9dbe2f7e18e73a9d61ee02980790335164;hp=dd2b3f1bb21d54662bda9cf7e07a5a518e87c9c1;hpb=020e579626d686ec20ecd9f0cc4c8313f474e152;p=unres.git diff --git a/source/unres/src-HCD-5D/cartder.F b/source/unres/src-HCD-5D/cartder.F index dd2b3f1..36b4e63 100644 --- a/source/unres/src-HCD-5D/cartder.F +++ b/source/unres/src-HCD-5D/cartder.F @@ -35,22 +35,33 @@ * * 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 * @@ -80,48 +91,18 @@ 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 @@ -132,8 +113,10 @@ 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) * @@ -153,6 +136,7 @@ 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. @@ -164,7 +148,7 @@ 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) * @@ -177,7 +161,8 @@ cd print *,ind1+1,(dxdv(j,ind1+1),j=1,3) 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) @@ -207,26 +192,29 @@ cd print *,ind1+1,(dxdv(j,ind1+1),j=1,3) * 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 @@ -235,9 +223,21 @@ cd print '(3f8.3)',(dcdv(k,ind1),k=1,3) 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 @@ -247,9 +247,11 @@ cd print '(3f8.3)',(dcdv(k,ind1),k=1,3) 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 @@ -259,6 +261,46 @@ cd print '(3f8.3)',(dcdv(k,ind1),k=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: * @@ -271,44 +313,43 @@ c dsci=dsc(itype(i)) 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 -