source/unres/src_MD/cartder.F

   1       subroutine cartder
   2 ***********************************************************************
   3 * This subroutine calculates the derivatives of the consecutive virtual
   4 * bond vectors and the SC vectors in the virtual-bond angles theta and
   5 * virtual-torsional angles phi, as well as the derivatives of SC vectors
   6 * in the angles alpha and omega, describing the location of a side chain
   7 * in its local coordinate system.
   8 *
   9 * The derivatives are stored in the following arrays:
  10 *
  11 * DDCDV - the derivatives of virtual-bond vectors DC in theta and phi.
  12 * The structure is as follows:
  13 *
  14 * dDC(x,2)/dT(3),...,dDC(z,2)/dT(3),0,             0,             0
  15 * 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)
  16 *         . . . . . . . . . . . .  . . . . . .
  17 * 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)
  18 *                          .
  19 *                          .
  20 *                          .
  21 * 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)
  22 *
  23 * DXDV - the derivatives of the side-chain vectors in theta and phi.
  24 * The structure is same as above.
  25 *
  26 * DCDS - the derivatives of the side chain vectors in the local spherical
  27 * andgles alph and omega:
  28 *
  29 * 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)
  30 * 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)
  31 *                          .
  32 *                          .
  33 *                          .
  34 * 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)
  35 *
  36 * Version of March '95, based on an early version of November '91.
  37 *
  38 ***********************************************************************
  39       implicit real*8 (a-h,o-z)
  40       include 'DIMENSIONS'
  41       include 'COMMON.VAR'
  42       include 'COMMON.CHAIN'
  43       include 'COMMON.DERIV'
  44       include 'COMMON.GEO'
  45       include 'COMMON.LOCAL'
  46       include 'COMMON.INTERACT'
  47       dimension drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),temp(3,3),
  48      &     fromto(3,3,maxdim),prordt(3,3,maxres),prodrt(3,3,maxres)
  49       dimension xx(3),xx1(3)
  50 c      common /przechowalnia/ fromto
  51 * get the position of the jth ijth fragment of the chain coordinate system
  52 * in the fromto array.
  53       indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
  54 *
  55 * calculate the derivatives of transformation matrix elements in theta
  56 *
  57       do i=1,nres-2
  58         rdt(1,1,i)=-rt(1,2,i)
  59         rdt(1,2,i)= rt(1,1,i)
  60         rdt(1,3,i)= 0.0d0
  61         rdt(2,1,i)=-rt(2,2,i)
  62         rdt(2,2,i)= rt(2,1,i)
  63         rdt(2,3,i)= 0.0d0
  64         rdt(3,1,i)=-rt(3,2,i)
  65         rdt(3,2,i)= rt(3,1,i)
  66         rdt(3,3,i)= 0.0d0
  67       enddo
  68 *
  69 * derivatives in phi
  70 *
  71       do i=2,nres-2
  72         drt(1,1,i)= 0.0d0
  73         drt(1,2,i)= 0.0d0
  74         drt(1,3,i)= 0.0d0
  75         drt(2,1,i)= rt(3,1,i)
  76         drt(2,2,i)= rt(3,2,i)
  77         drt(2,3,i)= rt(3,3,i)
  78         drt(3,1,i)=-rt(2,1,i)
  79         drt(3,2,i)=-rt(2,2,i)
  80         drt(3,3,i)=-rt(2,3,i)
  81       enddo
  82 *
  83 * generate the matrix products of type r(i)t(i)...r(j)t(j)
  84 *
  85       do i=2,nres-2
  86         ind=indmat(i,i+1)
  87         do k=1,3
  88           do l=1,3
  89             temp(k,l)=rt(k,l,i)
  90           enddo
  91         enddo
  92         do k=1,3
  93           do l=1,3
  94             fromto(k,l,ind)=temp(k,l)
  95           enddo
  96         enddo
  97         do j=i+1,nres-2
  98           ind=indmat(i,j+1)
  99           do k=1,3
 100             do l=1,3
 101               dpkl=0.0d0
 102               do m=1,3
 103                 dpkl=dpkl+temp(k,m)*rt(m,l,j)
 104               enddo
 105               dp(k,l)=dpkl
 106               fromto(k,l,ind)=dpkl
 107             enddo
 108           enddo
 109           do k=1,3
 110             do l=1,3
 111               temp(k,l)=dp(k,l)
 112             enddo
 113           enddo
 114         enddo
 115       enddo
 116 *
 117 * Calculate derivatives.
 118 *
 119       ind1=0
 120       do i=1,nres-2
 121         ind1=ind1+1
 122 *
 123 * Derivatives of DC(i+1) in theta(i+2)
 124 *
 125         do j=1,3
 126           do k=1,2
 127             dpjk=0.0D0
 128             do l=1,3
 129               dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
 130             enddo
 131             dp(j,k)=dpjk
 132             prordt(j,k,i)=dp(j,k)
 133           enddo
 134           dp(j,3)=0.0D0
 135           dcdv(j,ind1)=vbld(i+1)*dp(j,1)
 136         enddo
 137 *
 138 * Derivatives of SC(i+1) in theta(i+2)
 139 *
 140         xx1(1)=-0.5D0*xloc(2,i+1)
 141         xx1(2)= 0.5D0*xloc(1,i+1)
 142         do j=1,3
 143           xj=0.0D0
 144           do k=1,2
 145             xj=xj+r(j,k,i)*xx1(k)
 146           enddo
 147           xx(j)=xj
 148         enddo
 149         do j=1,3
 150           rj=0.0D0
 151           do k=1,3
 152             rj=rj+prod(j,k,i)*xx(k)
 153           enddo
 154           dxdv(j,ind1)=rj
 155         enddo
 156 *
 157 * Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
 158 * than the other off-diagonal derivatives.
 159 *
 160         do j=1,3
 161           dxoiij=0.0D0
 162           do k=1,3
 163             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 164           enddo
 165           dxdv(j,ind1+1)=dxoiij
 166         enddo
 167 cd      print *,ind1+1,(dxdv(j,ind1+1),j=1,3)
 168 *
 169 * Derivatives of DC(i+1) in phi(i+2)
 170 *
 171         do j=1,3
 172           do k=1,3
 173             dpjk=0.0
 174             do l=2,3
 175               dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
 176             enddo
 177             dp(j,k)=dpjk
 178             prodrt(j,k,i)=dp(j,k)
 179           enddo
 180           dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
 181         enddo
 182 *
 183 * Derivatives of SC(i+1) in phi(i+2)
 184 *
 185         xx(1)= 0.0D0
 186         xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
 187         xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
 188         do j=1,3
 189           rj=0.0D0
 190           do k=2,3
 191             rj=rj+prod(j,k,i)*xx(k)
 192           enddo
 193           dxdv(j+3,ind1)=-rj
 194         enddo
 195 *
 196 * Derivatives of SC(i+1) in phi(i+3).
 197 *
 198         do j=1,3
 199           dxoiij=0.0D0
 200           do k=1,3
 201             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 202           enddo
 203           dxdv(j+3,ind1+1)=dxoiij
 204         enddo
 205 *
 206 * Calculate the derivatives of DC(i+1) and SC(i+1) in theta(i+3) thru
 207 * theta(nres) and phi(i+3) thru phi(nres).
 208 *
 209         do j=i+1,nres-2
 210           ind1=ind1+1
 211           ind=indmat(i+1,j+1)
 212 cd        print *,'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
 213           do k=1,3
 214             do l=1,3
 215               tempkl=0.0D0
 216               do m=1,2
 217                 tempkl=tempkl+prordt(k,m,i)*fromto(m,l,ind)
 218               enddo
 219               temp(k,l)=tempkl
 220             enddo
 221           enddo
 222 cd        print '(9f8.3)',((fromto(k,l,ind),l=1,3),k=1,3)
 223 cd        print '(9f8.3)',((prod(k,l,i),l=1,3),k=1,3)
 224 cd        print '(9f8.3)',((temp(k,l),l=1,3),k=1,3)
 225 * Derivatives of virtual-bond vectors in theta
 226           do k=1,3
 227             dcdv(k,ind1)=vbld(i+1)*temp(k,1)
 228           enddo
 229 cd        print '(3f8.3)',(dcdv(k,ind1),k=1,3)
 230 * Derivatives of SC vectors in theta
 231           do k=1,3
 232             dxoijk=0.0D0
 233             do l=1,3
 234               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 235             enddo
 236             dxdv(k,ind1+1)=dxoijk
 237           enddo
 238 *
 239 *--- Calculate the derivatives in phi
 240 *
 241           do k=1,3
 242             do l=1,3
 243               tempkl=0.0D0
 244               do m=1,3
 245                 tempkl=tempkl+prodrt(k,m,i)*fromto(m,l,ind)
 246               enddo
 247               temp(k,l)=tempkl
 248             enddo
 249           enddo
 250           do k=1,3
 251             dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
 252           enddo
 253           do k=1,3
 254             dxoijk=0.0D0
 255             do l=1,3
 256               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 257             enddo
 258             dxdv(k+3,ind1+1)=dxoijk
 259           enddo
 260         enddo
 261       enddo
 262 *
 263 * Derivatives in alpha and omega:
 264 *
 265       do i=2,nres-1
 266 c       dsci=dsc(itype(i))
 267         dsci=vbld(i+nres)
 268 #ifdef OSF
 269         alphi=alph(i)
 270         omegi=omeg(i)
 271         if(alphi.ne.alphi) alphi=100.0
 272         if(omegi.ne.omegi) omegi=-100.0
 273 #else
 274         alphi=alph(i)
 275         omegi=omeg(i)
 276 #endif
 277 cd      print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
 278         cosalphi=dcos(alphi)
 279         sinalphi=dsin(alphi)
 280         cosomegi=dcos(omegi)
 281         sinomegi=dsin(omegi)
 282         temp(1,1)=-dsci*sinalphi
 283         temp(2,1)= dsci*cosalphi*cosomegi
 284         temp(3,1)=-dsci*cosalphi*sinomegi
 285         temp(1,2)=0.0D0
 286         temp(2,2)=-dsci*sinalphi*sinomegi
 287         temp(3,2)=-dsci*sinalphi*cosomegi
 288         theta2=pi-0.5D0*theta(i+1)
 289         cost2=dcos(theta2)
 290         sint2=dsin(theta2)
 291         jjj=0
 292 cd      print *,((temp(l,k),l=1,3),k=1,2)
 293         do j=1,2
 294           xp=temp(1,j)
 295           yp=temp(2,j)
 296           xxp= xp*cost2+yp*sint2
 297           yyp=-xp*sint2+yp*cost2
 298           zzp=temp(3,j)
 299           xx(1)=xxp
 300           xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
 301           xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
 302           do k=1,3
 303             dj=0.0D0
 304             do l=1,3
 305               dj=dj+prod(k,l,i-1)*xx(l)
 306             enddo
 307             dxds(jjj+k,i)=dj
 308           enddo
 309           jjj=jjj+3
 310         enddo
 311       enddo
 312       return
 313       end
 314