source/unres/src-HCD-5D/cartder.F.org

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