Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid
[unres.git] / source / wham / src-NEWSC / cartder.f
diff --git a/source/wham/src-NEWSC/cartder.f b/source/wham/src-NEWSC/cartder.f
new file mode 100755 (executable)
index 0000000..ed14f18
--- /dev/null
@@ -0,0 +1,306 @@
+      subroutine cartder 
+      implicit real*8 (a-h,o-z)
+***********************************************************************
+* This subroutine calculates the derivatives of the consecutive virtual
+* bond vectors and the SC vectors in the virtual-bond angles theta and
+* virtual-torsional angles phi, as well as the derivatives of SC vectors
+* in the angles alpha and omega, describing the location of a side chain
+* in its local coordinate system.
+*
+* The derivatives are stored in the following arrays:
+*
+* DDCDV - the derivatives of virtual-bond vectors DC in theta and phi.
+* The structure is as follows:
+*
+* dDC(x,2)/dT(3),...,dDC(z,2)/dT(3),0,             0,             0
+* dDC(x,3)/dT(4),...,dDC(z,3)/dT(4),dDC(x,3)/dP(4),dDC(y,4)/dP(4),dDC(z,4)/dP(4)
+*         . . . . . . . . . . . .  . . . . . .
+* dDC(x,N-1)/dT(4),...,dDC(z,N-1)/dT(4),dDC(x,N-1)/dP(4),dDC(y,N-1)/dP(4),dDC(z,N-1)/dP(4)
+*                          .
+*                          .
+*                          .
+* dDC(x,N-1)/dT(N),...,dDC(z,N-1)/dT(N),dDC(x,N-1)/dP(N),dDC(y,N-1)/dP(N),dDC(z,N-1)/dP(N)
+*
+* DXDV - the derivatives of the side-chain vectors in theta and phi. 
+* The structure is same as above.
+*
+* DCDS - the derivatives of the side chain vectors in the local spherical
+* andgles alph and omega:
+*
+* dX(x,2)/dA(2),dX(y,2)/dA(2),dX(z,2)/dA(2),dX(x,2)/dO(2),dX(y,2)/dO(2),dX(z,2)/dO(2)
+* dX(x,3)/dA(3),dX(y,3)/dA(3),dX(z,3)/dA(3),dX(x,3)/dO(3),dX(y,3)/dO(3),dX(z,3)/dO(3)
+*                          .
+*                          .
+*                          .
+* dX(x,N-1)/dA(N-1),dX(y,N-1)/dA(N-1),dX(z,N-1)/dA(N-1),dX(x,N-1)/dO(N-1),dX(y,N-1)/dO(N-1),dX(z,N-1)/dO(N-1)
+*
+* Version of March '95, based on an early version of November '91.
+*
+*********************************************************************** 
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      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)
+* get the position of the jth ijth fragment of the chain coordinate system      
+* in the fromto array.
+      indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
+*
+* calculate the derivatives of transformation matrix elements in theta
+*
+      do i=1,nres-2
+        rdt(1,1,i)=-rt(1,2,i)
+        rdt(1,2,i)= rt(1,1,i)
+        rdt(1,3,i)= 0.0d0
+        rdt(2,1,i)=-rt(2,2,i)
+        rdt(2,2,i)= rt(2,1,i)
+        rdt(2,3,i)= 0.0d0
+        rdt(3,1,i)=-rt(3,2,i)
+        rdt(3,2,i)= rt(3,1,i)
+        rdt(3,3,i)= 0.0d0
+      enddo
+*
+* derivatives in phi
+*
+      do i=2,nres-2
+        drt(1,1,i)= 0.0d0
+        drt(1,2,i)= 0.0d0
+        drt(1,3,i)= 0.0d0
+        drt(2,1,i)= rt(3,1,i)
+        drt(2,2,i)= rt(3,2,i)
+        drt(2,3,i)= rt(3,3,i)
+        drt(3,1,i)=-rt(2,1,i)
+        drt(3,2,i)=-rt(2,2,i)
+        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
+*
+* Derivatives of DC(i+1) in theta(i+2)
+*
+        do j=1,3
+          do k=1,2
+            dpjk=0.0D0
+            do l=1,3
+              dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
+            enddo
+            dp(j,k)=dpjk
+            prordt(j,k,i)=dp(j,k)
+          enddo
+          dp(j,3)=0.0D0
+          dcdv(j,ind1)=vbl*dp(j,1)       
+        enddo
+*
+* Derivatives of SC(i+1) in theta(i+2)
+* 
+        xx1(1)=-0.5D0*xloc(2,i+1)
+        xx1(2)= 0.5D0*xloc(1,i+1)
+        do j=1,3
+          xj=0.0D0
+          do k=1,2
+            xj=xj+r(j,k,i)*xx1(k)
+          enddo
+          xx(j)=xj
+        enddo
+        do j=1,3
+          rj=0.0D0
+          do k=1,3
+            rj=rj+prod(j,k,i)*xx(k)
+          enddo
+          dxdv(j,ind1)=rj
+        enddo
+*
+* Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
+* than the other off-diagonal derivatives.
+*
+        do j=1,3
+          dxoiij=0.0D0
+          do k=1,3
+            dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
+          enddo
+          dxdv(j,ind1+1)=dxoiij
+        enddo
+cd      print *,ind1+1,(dxdv(j,ind1+1),j=1,3)
+*
+* Derivatives of DC(i+1) in phi(i+2)
+*
+        do j=1,3
+          do k=1,3
+            dpjk=0.0
+            do l=2,3
+              dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
+            enddo
+            dp(j,k)=dpjk
+            prodrt(j,k,i)=dp(j,k)
+          enddo 
+          dcdv(j+3,ind1)=vbl*dp(j,1)
+        enddo
+*
+* Derivatives of SC(i+1) in phi(i+2)
+*
+        xx(1)= 0.0D0 
+        xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
+        xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
+        do j=1,3
+          rj=0.0D0
+          do k=2,3
+            rj=rj+prod(j,k,i)*xx(k)
+          enddo
+          dxdv(j+3,ind1)=-rj
+        enddo
+*
+* Derivatives of SC(i+1) in phi(i+3).
+*
+        do j=1,3
+          dxoiij=0.0D0
+          do k=1,3
+            dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
+          enddo
+          dxdv(j+3,ind1+1)=dxoiij
+        enddo
+*
+* Calculate the derivatives of DC(i+1) and SC(i+1) in theta(i+3) thru 
+* 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
+          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)
+              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)
+* Derivatives of virtual-bond vectors in theta
+          do k=1,3
+            dcdv(k,ind1)=vbl*temp(k,1)
+          enddo
+cd        print '(3f8.3)',(dcdv(k,ind1),k=1,3)
+* Derivatives of SC vectors in theta
+          do k=1,3
+            dxoijk=0.0D0
+            do l=1,3
+              dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
+            enddo
+            dxdv(k,ind1+1)=dxoijk
+          enddo
+*
+*--- Calculate the derivatives in phi
+*
+          do k=1,3
+            do l=1,3
+              tempkl=0.0D0
+              do m=1,3
+                tempkl=tempkl+prodrt(k,m,i)*fromto(m,l,ind)
+              enddo
+              temp(k,l)=tempkl
+            enddo
+          enddo
+          do k=1,3
+            dcdv(k+3,ind1)=vbl*temp(k,1)
+         enddo
+          do k=1,3
+            dxoijk=0.0D0
+            do l=1,3
+              dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
+            enddo
+            dxdv(k+3,ind1+1)=dxoijk
+          enddo
+        enddo
+      enddo
+*
+* Derivatives in alpha and omega:
+*
+      do i=2,nres-1
+       dsci=dsc(itype(i))
+       alphi=alph(i)
+       omegi=omeg(i)
+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
+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)
+            enddo
+           dxds(jjj+k,i)=dj
+          enddo
+         jjj=jjj+3
+       enddo
+      enddo
+      return
+      end
+