source/unres/src-HCD-5D/cart2intgrad.F

   1       subroutine cart2intgrad(n,g)
   2 ***********************************************************************
   3 * This subroutine thransforms the gradient in virtual-bond vectors to
   4 * that in the backbone and side-chain angular variables.
   5 * Adapted from the cartder subroutine.
   6 *
   7 * 03/11/20 Adam. Array fromto eliminated, computed on the fly
   8 *     Fixed the problem with vbld indices, which caused errors in
   9 *     derivatives when the backbone virtual bond lengths were not equal.
  10 ***********************************************************************
  11       implicit none
  12       include 'DIMENSIONS'
  13       include 'COMMON.IOUNITS'
  14       include 'COMMON.VAR'
  15       include 'COMMON.CHAIN'
  16       include 'COMMON.DERIV'
  17       include 'COMMON.GEO'
  18       include 'COMMON.LOCAL'
  19       include 'COMMON.INTERACT'
  20       integer n
  21       double precision g(n)
  22       double precision drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),
  23      &temp(3,3),prordt(3,3,maxres),prodrt(3,3,maxres)
  24       double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp
  25       double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,
  26      & cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl
  27       double precision fromto(3,3),aux(6)
  28       integer i,ii,j,jjj,k,l,m,indi,ind,ind1
  29 * get the position of the jth ijth fragment of the chain coordinate system
  30 * in the fromto array.
  31 c      integer indmat
  32 c      indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
  33 c      call chainbuild_extconf
  34 c      call cartprint
  35 c      call intout
  36       g=0.0d0
  37 * 3/13/20 Adam: Skip calculating backbone derivatives if SC only
  38 * requested.
  39       if (sideonly) goto 10
  40 *
  41 * calculate the derivatives of transformation matrix elements in theta
  42 *
  43       do i=1,nres-2
  44         rdt(1,1,i)=-rt(1,2,i)
  45         rdt(1,2,i)= rt(1,1,i)
  46         rdt(1,3,i)= 0.0d0
  47         rdt(2,1,i)=-rt(2,2,i)
  48         rdt(2,2,i)= rt(2,1,i)
  49         rdt(2,3,i)= 0.0d0
  50         rdt(3,1,i)=-rt(3,2,i)
  51         rdt(3,2,i)= rt(3,1,i)
  52         rdt(3,3,i)= 0.0d0
  53       enddo
  54 *
  55 * derivatives in phi
  56 *
  57       do i=2,nres-2
  58         drt(1,1,i)= 0.0d0
  59         drt(1,2,i)= 0.0d0
  60         drt(1,3,i)= 0.0d0
  61         drt(2,1,i)= rt(3,1,i)
  62         drt(2,2,i)= rt(3,2,i)
  63         drt(2,3,i)= rt(3,3,i)
  64         drt(3,1,i)=-rt(2,1,i)
  65         drt(3,2,i)=-rt(2,2,i)
  66         drt(3,3,i)=-rt(2,3,i)
  67       enddo
  68 *
  69 * Calculate backbone derivatives.
  70 * This code invlves N^2 effort and should be parallelized, to be done
  71 * later.
  72       ind1=0
  73       do i=1,nres-2
  74         ind1=ind1+1
  75 *
  76 * Derivatives of DC(i+1) in theta(i+2)
  77 *
  78 c        write (iout,*) "theta i",i
  79 c        write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
  80 c        write(iout,'(7hrdt    9f10.5)')((rdt(k,l,i),l=1,3),k=1,3)
  81 c        write(iout,*) "vbld",vbld(i+2)
  82         if (n.gt.nphi) then
  83
  84         do j=1,3
  85           do k=1,2
  86             dpjk=0.0D0
  87             do l=1,3
  88               dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
  89             enddo
  90             dp(j,k)=dpjk
  91             prordt(j,k,i)=dp(j,k)
  92           enddo
  93           dp(j,3)=0.0D0
  94 c          dcdv(j,ind1)=vbld(i+2)*dp(j,1)
  95           g(nphi+i)=g(nphi+i)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
  96         enddo
  97 c        write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
  98 *
  99 * Derivatives of SC(i+1) in theta(i+2)
 100 *
 101         xx1(1)=-0.5D0*xloc(2,i+1)
 102         xx1(2)= 0.5D0*xloc(1,i+1)
 103         do j=1,3
 104           xj=0.0D0
 105           do k=1,2
 106             xj=xj+r(j,k,i)*xx1(k)
 107           enddo
 108           xx(j)=xj
 109         enddo
 110         do j=1,3
 111           rj=0.0D0
 112           do k=1,3
 113             rj=rj+prod(j,k,i)*xx(k)
 114           enddo
 115 c          dxdv(j,ind1)=rj
 116 c          write (iout,*) "1:i",i," j",i+1,"ind1",ind1," dxdthet",rj,
 117 c     &      " gradx",gradx(j,i+1,icg)
 118           g(nphi+i)=g(nphi+i)+rj*gradx(j,i+1,icg)
 119         enddo
 120 c        write (iout,*) "dxdv",(dxdv(j,ind1),j=1,3)
 121 *
 122 * Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
 123 * than the other off-diagonal derivatives.
 124 *
 125         if (i.lt.nres-2) then
 126         do j=1,3
 127           dxoiij=0.0D0
 128           do k=1,3
 129             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 130           enddo
 131 c          dxdv(j,ind1+1)=dxoiij
 132 c          write (iout,*) "2:i",i," j",i+1,"ind1",ind1+1,
 133 c     &      " dxdthet",dxoiij," gradx",gradx(j,i+2,icg)
 134           g(nphi+i)=g(nphi+i)+dxoiij*gradx(j,i+2,icg)
 135         enddo
 136         endif
 137 c        write(iout,*)ind1+1,(dxdv(j,ind1+1),j=1,3)
 138
 139         endif
 140 *
 141 * Derivatives of DC(i+1) in phi(i+2)
 142 *
 143         if (i.gt.1) then
 144         do j=1,3
 145           do k=1,3
 146             dpjk=0.0
 147             do l=2,3
 148               dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
 149             enddo
 150             dp(j,k)=dpjk
 151             prodrt(j,k,i)=dp(j,k)
 152           enddo
 153 c          dcdv(j+3,ind1)=vbld(i+2)*dp(j,1)
 154           g(i-1)=g(i-1)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
 155         enddo
 156         endif
 157 *
 158 * Derivatives of SC(i+1) in phi(i+2)
 159 *
 160         xx(1)= 0.0D0
 161         xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
 162         xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
 163         if (i.gt.1) then
 164         do j=1,3
 165           rj=0.0D0
 166           do k=2,3
 167             rj=rj+prod(j,k,i)*xx(k)
 168           enddo
 169 c          dxdv(j+3,ind1)=-rj
 170 c          write (iout,*) "1:i",i," j",i+1,"ind1",ind1," dxdphi",-rj,
 171 c     &      " gradx",gradx(j,i+1,icg)
 172           g(i-1)=g(i-1)-rj*gradx(j,i+1,icg)
 173         enddo
 174         endif
 175 *
 176 * Derivatives of SC(i+1) in phi(i+3).
 177 *
 178         if (i.gt.1) then
 179         do j=1,3
 180           dxoiij=0.0D0
 181           do k=1,3
 182             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
 183           enddo
 184 c          dxdv(j+3,ind1+1)=dxoiij
 185           g(i-1)=g(i-1)+dxoiij*gradx(j,i+2,icg)
 186 c          write (iout,*) "2:i",i," j",i+2," ind1",ind1+1,
 187 c     &      " dxdphi",dxoiij," gradx",gradx(j,i+2,icg)
 188         enddo
 189         endif
 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 c          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           if (n.gt.nphi) then
 213 * Derivatives of virtual-bond vectors in theta
 214           do k=1,3
 215 c            dcdv(k,ind1)=vbld(j+2)*temp(k,1)
 216             g(nphi+i)=g(nphi+i)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
 217           enddo
 218 c          write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
 219 * Derivatives of SC vectors in theta
 220           do k=1,3
 221             dxoijk=0.0D0
 222             do l=1,3
 223               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 224             enddo
 225 c            dxdv(k,ind1+1)=dxoijk
 226 c            write (iout,*) "3:i",i+1," j",j+2,"ind1",ind1+1,
 227 c     &      " dxdthet",dxoijk," gradx",gradx(k,j+2,icg)
 228             g(nphi+i)=g(nphi+i)+dxoijk*gradx(k,j+2,icg)
 229           enddo
 230 c          write(iout,'(7htheta  3f10.5)')(dxdv(k,ind1),k=1,3)
 231           endif
 232 *
 233 *--- Calculate the derivatives in phi
 234 *
 235           do k=1,3
 236             do l=1,3
 237               tempkl=0.0D0
 238               do m=1,3
 239                 tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
 240               enddo
 241               temp(k,l)=tempkl
 242             enddo
 243           enddo
 244           if (i.gt.1) then
 245           do k=1,3
 246 c            dcdv(k+3,ind1)=vbld(j+2)*temp(k,1)
 247             g(i-1)=g(i-1)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
 248           enddo
 249           do k=1,3
 250             dxoijk=0.0D0
 251             do l=1,3
 252               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
 253             enddo
 254 c            dxdv(k+3,ind1+1)=dxoijk
 255             g(i-1)=g(i-1)+dxoijk*gradx(k,j+2,icg)
 256 c            write (iout,*) "3:i",i," j",j+2," ind1",ind1+1,
 257 c     &      " dxdphi",dxoijk," gradx",gradx(k,j+2,icg)
 258           enddo
 259           endif
 260         enddo
 261       enddo
 262
 263       if (nvar.le.nphi+ntheta) return
 264
 265    10 continue
 266 *
 267 * Derivatives in alpha and omega:
 268 *
 269       do i=2,nres-1
 270         if (iabs(itype(i)).eq.10 .or. itype(i).eq.ntyp1!) cycle
 271      &    .or. mask_side(i).eq.0 ) cycle
 272         ii=ialph(i,1)
 273         dsci=vbld(i+nres)
 274 #ifdef OSF
 275         alphi=alph(i)
 276         omegi=omeg(i)
 277         if(alphi.ne.alphi) alphi=100.0
 278         if(omegi.ne.omegi) omegi=-100.0
 279 #else
 280         alphi=alph(i)
 281         omegi=omeg(i)
 282 #endif
 283 cd      print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
 284         cosalphi=dcos(alphi)
 285         sinalphi=dsin(alphi)
 286         cosomegi=dcos(omegi)
 287         sinomegi=dsin(omegi)
 288         temp(1,1)=-dsci*sinalphi
 289         temp(2,1)= dsci*cosalphi*cosomegi
 290         temp(3,1)=-dsci*cosalphi*sinomegi
 291         temp(1,2)=0.0D0
 292         temp(2,2)=-dsci*sinalphi*sinomegi
 293         temp(3,2)=-dsci*sinalphi*cosomegi
 294         theta2=pi-0.5D0*theta(i+1)
 295         cost2=dcos(theta2)
 296         sint2=dsin(theta2)
 297         jjj=0
 298 cd      print *,((temp(l,k),l=1,3),k=1,2)
 299         do j=1,2
 300           xp=temp(1,j)
 301           yp=temp(2,j)
 302           xxp= xp*cost2+yp*sint2
 303           yyp=-xp*sint2+yp*cost2
 304           zzp=temp(3,j)
 305           xx(1)=xxp
 306           xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
 307           xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
 308           do k=1,3
 309             dj=0.0D0
 310             do l=1,3
 311               dj=dj+prod(k,l,i-1)*xx(l)
 312             enddo
 313 c            dxds(jjj+k,i)=dj
 314             aux(jjj+k)=dj
 315           enddo
 316           jjj=jjj+3
 317         enddo
 318         do k=1,3
 319           g(ii)=g(ii)+aux(k)*gradx(k,i,icg)
 320           g(ii+nside)=g(ii+nside)+aux(k+3)*gradx(k,i,icg)
 321         enddo
 322       enddo
 323       return
 324       end
 325 c-----------------------------------------------------------------------------
 326       subroutine build_fromto(i,j,fromto)
 327       implicit none
 328       include 'DIMENSIONS'
 329       include 'COMMON.CHAIN'
 330       include 'COMMON.IOUNITS'
 331       integer i,j,jj,k,l,m
 332       double precision fromto(3,3),temp(3,3),dp(3,3)
 333       double precision dpkl
 334       save temp
 335 *
 336 * generate the matrix products of type r(i)t(i)...r(j)t(j) on the fly
 337 *
 338 c      write (iout,*) "temp on entry"
 339 c      write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
 340 c      do i=2,nres-2
 341 c        ind=indmat(i,i+1)
 342       if (j.eq.i+1) then
 343         do k=1,3
 344           do l=1,3
 345             temp(k,l)=rt(k,l,i)
 346           enddo
 347         enddo
 348         do k=1,3
 349           do l=1,3
 350             fromto(k,l)=temp(k,l)
 351           enddo
 352         enddo
 353       else
 354 c        do j=i+1,nres-2
 355 c          ind=indmat(i,j+1)
 356           do k=1,3
 357             do l=1,3
 358               dpkl=0.0d0
 359               do m=1,3
 360                 dpkl=dpkl+temp(k,m)*rt(m,l,j-1)
 361               enddo
 362               dp(k,l)=dpkl
 363               fromto(k,l)=dpkl
 364             enddo
 365           enddo
 366           do k=1,3
 367             do l=1,3
 368               temp(k,l)=dp(k,l)
 369             enddo
 370           enddo
 371       endif
 372 c      write (iout,*) "temp upon exit"
 373 c      write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
 374 c        enddo
 375 c      enddo
 376       return
 377       end