source/unres/src-HCD-5D/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 * 03/11/20 Adam. Array fromto eliminated, computed on the fly
  39 *     Fixed the problem with vbld indices, which caused errors in
  40 *     derivatives when the backbone virtual bond lengths were not equal.
  41 ***********************************************************************
  42       implicit none
  43       include 'DIMENSIONS'
  44       include 'COMMON.IOUNITS'
  45       include 'COMMON.VAR'
  46       include 'COMMON.CHAIN'
  47       include 'COMMON.DERIV'
  48       include 'COMMON.GEO'
  49       include 'COMMON.LOCAL'
  50       include 'COMMON.INTERACT'
  51       double precision drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),
  52      &temp(3,3),prordt(3,3,maxres),prodrt(3,3,maxres)
  53       double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp
  54       double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,
  55      & cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl
  56       double precision fromto(3,3)
  57       integer i,ii,j,jjj,k,l,m,indi,ind,ind1
  58 * get the position of the jth ijth fragment of the chain coordinate system
  59 * in the fromto array.
  60       integer indmat
  61       indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
  62       call chainbuild_extconf
  63       call cartprint
  64       call intout
  65 *
  66 * calculate the derivatives of transformation matrix elements in theta
  67 *
  68       do i=1,nres-2
  69         rdt(1,1,i)=-rt(1,2,i)
  70         rdt(1,2,i)= rt(1,1,i)
  71         rdt(1,3,i)= 0.0d0
  72         rdt(2,1,i)=-rt(2,2,i)
  73         rdt(2,2,i)= rt(2,1,i)
  74         rdt(2,3,i)= 0.0d0
  75         rdt(3,1,i)=-rt(3,2,i)
  76         rdt(3,2,i)= rt(3,1,i)
  77         rdt(3,3,i)= 0.0d0
  78       enddo
  79 *
  80 * derivatives in phi
  81 *
  82       do i=2,nres-2
  83         drt(1,1,i)= 0.0d0
  84         drt(1,2,i)= 0.0d0
  85         drt(1,3,i)= 0.0d0
  86         drt(2,1,i)= rt(3,1,i)
  87         drt(2,2,i)= rt(3,2,i)
  88         drt(2,3,i)= rt(3,3,i)
  89         drt(3,1,i)=-rt(2,1,i)
  90         drt(3,2,i)=-rt(2,2,i)
  91         drt(3,3,i)=-rt(2,3,i)
  92       enddo
  93 *
  94 * Calculate derivatives.
  95 *
  96       ind1=0
  97       do i=1,nres-2
  98         ind1=ind1+1
  99 *
 100 * Derivatives of DC(i+1) in theta(i+2)
 101 *
 102 c        write (iout,*) "theta i",i
 103 c        write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
 104 c        write(iout,'(7hrdt    9f10.5)')((rdt(k,l,i),l=1,3),k=1,3)
 105 c        write(iout,*) "vbld",vbld(i+2)
 106         do j=1,3
 107           do k=1,2
 108             dpjk=0.0D0
 109             do l=1,3
 110               dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
 111             enddo
 112             dp(j,k)=dpjk
 113             prordt(j,k,i)=dp(j,k)
 114           enddo
 115           dp(j,3)=0.0D0
 116 c          dcdv(j,ind1)=vbld(i+1)*dp(j,1)
 117           dcdv(j,ind1)=vbld(i+2)*dp(j,1)
 118         enddo
 119 c        write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
 120 *
 121 * Derivatives of SC(i+1) in theta(i+2)
 122 *
 123         xx1(1)=-0.5D0*xloc(2,i+1)
 124         xx1(2)= 0.5D0*xloc(1,i+1)
 125         do j=1,3
 126           xj=0.0D0
 127           do k=1,2
 128             xj=xj+r(j,k,i)*xx1(k)
 129           enddo
 130           xx(j)=xj
 131         enddo
 132         do j=1,3
 133           rj=0.0D0
 134           do k=1,3
 135             rj=rj+prod(j,k,i)*xx(k)
 136           enddo
 137           dxdv(j,ind1)=rj
 138         enddo
 139 c        write (iout,*) "dxdv",(dxdv(j,ind1),j=1,3)
 140 *
 141 * Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
 142 * than the other off-diagonal derivatives.
 143 *
 144         do j=1,3
 145           dxoiij=0.0D0
 146           do k=1,3
 147             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 148           enddo
 149           dxdv(j,ind1+1)=dxoiij
 150         enddo
 151 c        write(iout,*)ind1+1,(dxdv(j,ind1+1),j=1,3)
 152 *
 153 * Derivatives of DC(i+1) in phi(i+2)
 154 *
 155         do j=1,3
 156           do k=1,3
 157             dpjk=0.0
 158             do l=2,3
 159               dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
 160             enddo
 161             dp(j,k)=dpjk
 162             prodrt(j,k,i)=dp(j,k)
 163           enddo
 164 c          dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
 165           dcdv(j+3,ind1)=vbld(i+2)*dp(j,1)
 166         enddo
 167 *
 168 * Derivatives of SC(i+1) in phi(i+2)
 169 *
 170         xx(1)= 0.0D0
 171         xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
 172         xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
 173         do j=1,3
 174           rj=0.0D0
 175           do k=2,3
 176             rj=rj+prod(j,k,i)*xx(k)
 177           enddo
 178           dxdv(j+3,ind1)=-rj
 179         enddo
 180 *
 181 * Derivatives of SC(i+1) in phi(i+3).
 182 *
 183         do j=1,3
 184           dxoiij=0.0D0
 185           do k=1,3
 186             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 187           enddo
 188           dxdv(j+3,ind1+1)=dxoiij
 189         enddo
 190 *
 191 * Calculate the derivatives of DC(i+1) and SC(i+1) in theta(i+3) thru
 192 * theta(nres) and phi(i+3) thru phi(nres).
 193 *
 194         do j=i+1,nres-2
 195           ind1=ind1+1
 196           ind=indmat(i+1,j+1)
 197 c          write(iout,*)'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
 198           call build_fromto(i+1,j+1,fromto)
 199 c          write(iout,'(7hfromto 9f10.5)')((fromto(k,l),l=1,3),k=1,3)
 200           do k=1,3
 201             do l=1,3
 202               tempkl=0.0D0
 203               do m=1,2
 204                 tempkl=tempkl+prordt(k,m,i)*fromto(m,l)
 205               enddo
 206               temp(k,l)=tempkl
 207             enddo
 208           enddo
 209 c          write(iout,'(7hfromto 9f10.5)')((fromto(k,l,ind),l=1,3),k=1,3)
 210 c          write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
 211 c          write(iout,'(7htemp   9f10.5)')((temp(k,l),l=1,3),k=1,3)
 212 * Derivatives of virtual-bond vectors in theta
 213           do k=1,3
 214 c            dcdv(k,ind1)=vbld(i+1)*temp(k,1)
 215             dcdv(k,ind1)=vbld(j+2)*temp(k,1)
 216           enddo
 217 c          write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
 218 * Derivatives of SC vectors in theta
 219           do k=1,3
 220             dxoijk=0.0D0
 221             do l=1,3
 222               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 223             enddo
 224             dxdv(k,ind1+1)=dxoijk
 225           enddo
 226 c          write(iout,'(7htheta  3f10.5)')(dxdv(k,ind1),k=1,3)
 227 *
 228 *--- Calculate the derivatives in phi
 229 *
 230 #ifdef FIVEDIAG
 231           do k=1,3
 232             do l=1,3
 233               tempkl=0.0D0
 234               do m=1,3
 235                 tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
 236               enddo
 237               temp(k,l)=tempkl
 238             enddo
 239           enddo
 240 #else
 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 #endif
 251           do k=1,3
 252 c            dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
 253             dcdv(k+3,ind1)=vbld(j+2)*temp(k,1)
 254           enddo
 255           do k=1,3
 256             dxoijk=0.0D0
 257             do l=1,3
 258               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 259             enddo
 260             dxdv(k+3,ind1+1)=dxoijk
 261           enddo
 262         enddo
 263       enddo
 264 #ifdef DEBUG
 265       write (iout,*)
 266       write (iout,'(a)') '****************** ddc/dtheta'
 267       write (iout,*)
 268       do i=1,nres-2
 269         do j=i+1,nres-1
 270           ii = indmat(i,j)
 271           write (iout,'(2i4,3e14.6)') i,j,(dcdv(k,ii),k=1,3)
 272         enddo
 273       enddo
 274       write (iout,*)
 275       write (iout,'(a)') '******************* ddc/dphi'
 276       write (iout,*)
 277       do i=1,nres-3
 278         do j=i+2,nres-1
 279           ii = indmat(i+1,j)
 280           write (iout,'(2i4,3e14.6)') i,j,(dcdv(k+3,ii),k=1,3)
 281           write (iout,'(a)')
 282         enddo
 283       enddo
 284       write (iout,'(a)')
 285       write (iout,'(a)') '**************** dx/dtheta'
 286       write (iout,'(a)')
 287       do i=3,nres
 288         do j=i-1,nres-1
 289           ii = indmat(i-2,j)
 290           write (iout,'(2i4,3e14.6)') i,j,(dxdv(k,ii),k=1,3)
 291         enddo
 292       enddo
 293       write (iout,'(a)')
 294       write (iout,'(a)') '***************** dx/dphi'
 295       write (iout,'(a)')
 296       do i=4,nres
 297         do j=i-1,nres-1
 298           ii = indmat(i-2,j)
 299           write (iout,'(2i4,3e14.6)') i,j,(dxdv(k+3,ii),k=1,3)
 300           write(iout,'(a)')
 301         enddo
 302       enddo
 303 #endif
 304 *
 305 * Derivatives in alpha and omega:
 306 *
 307       do i=2,nres-1
 308 c       dsci=dsc(itype(i))
 309         dsci=vbld(i+nres)
 310 #ifdef OSF
 311         alphi=alph(i)
 312         omegi=omeg(i)
 313         if(alphi.ne.alphi) alphi=100.0
 314         if(omegi.ne.omegi) omegi=-100.0
 315 #else
 316         alphi=alph(i)
 317         omegi=omeg(i)
 318 #endif
 319 cd      print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
 320         cosalphi=dcos(alphi)
 321         sinalphi=dsin(alphi)
 322         cosomegi=dcos(omegi)
 323         sinomegi=dsin(omegi)
 324         temp(1,1)=-dsci*sinalphi
 325         temp(2,1)= dsci*cosalphi*cosomegi
 326         temp(3,1)=-dsci*cosalphi*sinomegi
 327         temp(1,2)=0.0D0
 328         temp(2,2)=-dsci*sinalphi*sinomegi
 329         temp(3,2)=-dsci*sinalphi*cosomegi
 330         theta2=pi-0.5D0*theta(i+1)
 331         cost2=dcos(theta2)
 332         sint2=dsin(theta2)
 333         jjj=0
 334 cd      print *,((temp(l,k),l=1,3),k=1,2)
 335         do j=1,2
 336           xp=temp(1,j)
 337           yp=temp(2,j)
 338           xxp= xp*cost2+yp*sint2
 339           yyp=-xp*sint2+yp*cost2
 340           zzp=temp(3,j)
 341           xx(1)=xxp
 342           xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
 343           xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
 344           do k=1,3
 345             dj=0.0D0
 346             do l=1,3
 347               dj=dj+prod(k,l,i-1)*xx(l)
 348             enddo
 349             dxds(jjj+k,i)=dj
 350           enddo
 351           jjj=jjj+3
 352         enddo
 353       enddo
 354       return
 355       end
 356 c-----------------------------------------------------------------------------
 357       subroutine build_fromto(i,j,fromto)
 358       implicit none
 359       include 'DIMENSIONS'
 360       include 'COMMON.CHAIN'
 361       include 'COMMON.IOUNITS'
 362       integer i,j,jj,k,l,m
 363       double precision fromto(3,3),temp(3,3),dp(3,3)
 364       double precision dpkl
 365       save temp
 366 *
 367 * generate the matrix products of type r(i)t(i)...r(j)t(j) on the fly
 368 *
 369 c      write (iout,*) "temp on entry"
 370 c      write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
 371 c      do i=2,nres-2
 372 c        ind=indmat(i,i+1)
 373       if (j.eq.i+1) then
 374         do k=1,3
 375           do l=1,3
 376             temp(k,l)=rt(k,l,i)
 377           enddo
 378         enddo
 379         do k=1,3
 380           do l=1,3
 381             fromto(k,l)=temp(k,l)
 382           enddo
 383         enddo
 384       else
 385 c        do j=i+1,nres-2
 386 c          ind=indmat(i,j+1)
 387           do k=1,3
 388             do l=1,3
 389               dpkl=0.0d0
 390               do m=1,3
 391                 dpkl=dpkl+temp(k,m)*rt(m,l,j-1)
 392               enddo
 393               dp(k,l)=dpkl
 394               fromto(k,l)=dpkl
 395             enddo
 396           enddo
 397           do k=1,3
 398             do l=1,3
 399               temp(k,l)=dp(k,l)
 400             enddo
 401           enddo
 402       endif
 403 c      write (iout,*) "temp upon exit"
 404 c      write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
 405 c        enddo
 406 c      enddo
 407       return
 408       end