source/unres/src_MD-M/cartder.F.diag

   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       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         write (iout,*) "theta i",i
 126         write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
 127         write(iout,'(7hrdt    9f10.5)')((rdt(k,l,i),l=1,3),k=1,3)
 128         write(iout,*) "vbld",vbld(i+1)
 129         do j=1,3
 130           do k=1,2
 131             dpjk=0.0D0
 132             do l=1,3
 133               dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
 134             enddo
 135             dp(j,k)=dpjk
 136             prordt(j,k,i)=dp(j,k)
 137           enddo
 138           dp(j,3)=0.0D0
 139           dcdv(j,ind1)=vbld(i+1)*dp(j,1)
 140         enddo
 141         write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
 142 *
 143 * Derivatives of SC(i+1) in theta(i+2)
 144 *
 145         xx1(1)=-0.5D0*xloc(2,i+1)
 146         xx1(2)= 0.5D0*xloc(1,i+1)
 147         do j=1,3
 148           xj=0.0D0
 149           do k=1,2
 150             xj=xj+r(j,k,i)*xx1(k)
 151           enddo
 152           xx(j)=xj
 153         enddo
 154         do j=1,3
 155           rj=0.0D0
 156           do k=1,3
 157             rj=rj+prod(j,k,i)*xx(k)
 158           enddo
 159           dxdv(j,ind1)=rj
 160         enddo
 161 *
 162 * Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
 163 * than the other off-diagonal derivatives.
 164 *
 165         do j=1,3
 166           dxoiij=0.0D0
 167           do k=1,3
 168             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 169           enddo
 170           dxdv(j,ind1+1)=dxoiij
 171         enddo
 172 cd      print *,ind1+1,(dxdv(j,ind1+1),j=1,3)
 173 *
 174 * Derivatives of DC(i+1) in phi(i+2)
 175 *
 176         do j=1,3
 177           do k=1,3
 178             dpjk=0.0
 179             do l=2,3
 180               dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
 181             enddo
 182             dp(j,k)=dpjk
 183             prodrt(j,k,i)=dp(j,k)
 184           enddo
 185           dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
 186         enddo
 187 *
 188 * Derivatives of SC(i+1) in phi(i+2)
 189 *
 190         xx(1)= 0.0D0
 191         xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
 192         xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
 193         do j=1,3
 194           rj=0.0D0
 195           do k=2,3
 196             rj=rj+prod(j,k,i)*xx(k)
 197           enddo
 198           dxdv(j+3,ind1)=-rj
 199         enddo
 200 *
 201 * Derivatives of SC(i+1) in phi(i+3).
 202 *
 203         do j=1,3
 204           dxoiij=0.0D0
 205           do k=1,3
 206             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 207           enddo
 208           dxdv(j+3,ind1+1)=dxoiij
 209         enddo
 210 *
 211 * Calculate the derivatives of DC(i+1) and SC(i+1) in theta(i+3) thru
 212 * theta(nres) and phi(i+3) thru phi(nres).
 213 *
 214         do j=i+1,nres-2
 215           ind1=ind1+1
 216           ind=indmat(i+1,j+1)
 217           write(iout,*)'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
 218 cd        print *,'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
 219           do k=1,3
 220             do l=1,3
 221               tempkl=0.0D0
 222               do m=1,2
 223                 tempkl=tempkl+prordt(k,m,i)*fromto(m,l,ind)
 224               enddo
 225               temp(k,l)=tempkl
 226             enddo
 227           enddo
 228           write(iout,'(7hfromto 9f10.5)')((fromto(k,l,ind),l=1,3),k=1,3)
 229           write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
 230           write(iout,'(7htemp   9f10.5)')((temp(k,l),l=1,3),k=1,3)
 231 cd        print '(9f8.3)',((fromto(k,l,ind),l=1,3),k=1,3)
 232 cd        print '(9f8.3)',((prod(k,l,i),l=1,3),k=1,3)
 233 cd        print '(9f8.3)',((temp(k,l),l=1,3),k=1,3)
 234 * Derivatives of virtual-bond vectors in theta
 235           do k=1,3
 236             dcdv(k,ind1)=vbld(i+1)*temp(k,1)
 237           enddo
 238           write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
 239 * Derivatives of SC vectors in theta
 240           do k=1,3
 241             dxoijk=0.0D0
 242             do l=1,3
 243               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 244             enddo
 245             dxdv(k,ind1+1)=dxoijk
 246           enddo
 247 *
 248 *--- Calculate the derivatives in phi
 249 *
 250           do k=1,3
 251             do l=1,3
 252               tempkl=0.0D0
 253               do m=1,3
 254                 tempkl=tempkl+prodrt(k,m,i)*fromto(m,l,ind)
 255               enddo
 256               temp(k,l)=tempkl
 257             enddo
 258           enddo
 259           do k=1,3
 260             dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
 261           enddo
 262           do k=1,3
 263             dxoijk=0.0D0
 264             do l=1,3
 265               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 266             enddo
 267             dxdv(k+3,ind1+1)=dxoijk
 268           enddo
 269         enddo
 270       enddo
 271 *
 272 * Derivatives in alpha and omega:
 273 *
 274       do i=2,nres-1
 275 c       dsci=dsc(itype(i))
 276         dsci=vbld(i+nres)
 277 #ifdef OSF
 278         alphi=alph(i)
 279         omegi=omeg(i)
 280         if(alphi.ne.alphi) alphi=100.0
 281         if(omegi.ne.omegi) omegi=-100.0
 282 #else
 283         alphi=alph(i)
 284         omegi=omeg(i)
 285 #endif
 286 cd      print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
 287         cosalphi=dcos(alphi)
 288         sinalphi=dsin(alphi)
 289         cosomegi=dcos(omegi)
 290         sinomegi=dsin(omegi)
 291         temp(1,1)=-dsci*sinalphi
 292         temp(2,1)= dsci*cosalphi*cosomegi
 293         temp(3,1)=-dsci*cosalphi*sinomegi
 294         temp(1,2)=0.0D0
 295         temp(2,2)=-dsci*sinalphi*sinomegi
 296         temp(3,2)=-dsci*sinalphi*cosomegi
 297         theta2=pi-0.5D0*theta(i+1)
 298         cost2=dcos(theta2)
 299         sint2=dsin(theta2)
 300         jjj=0
 301 cd      print *,((temp(l,k),l=1,3),k=1,2)
 302         do j=1,2
 303           xp=temp(1,j)
 304           yp=temp(2,j)
 305           xxp= xp*cost2+yp*sint2
 306           yyp=-xp*sint2+yp*cost2
 307           zzp=temp(3,j)
 308           xx(1)=xxp
 309           xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
 310           xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
 311           do k=1,3
 312             dj=0.0D0
 313             do l=1,3
 314               dj=dj+prod(k,l,i-1)*xx(l)
 315             enddo
 316             dxds(jjj+k,i)=dj
 317           enddo
 318           jjj=jjj+3
 319         enddo
 320       enddo
 321       return
 322       end
 323