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.
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 ***********************************************************************
13 include 'COMMON.IOUNITS'
15 include 'COMMON.CHAIN'
16 include 'COMMON.DERIV'
18 include 'COMMON.LOCAL'
19 include 'COMMON.INTERACT'
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.
32 c indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
33 c call chainbuild_extconf
37 * 3/13/20 Adam: Skip calculating backbone derivatives if SC only
41 * calculate the derivatives of transformation matrix elements in theta
69 * Calculate backbone derivatives.
70 * This code invlves N^2 effort and should be parallelized, to be done
76 * Derivatives of DC(i+1) in theta(i+2)
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)
88 dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
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)
97 c write(iout,'(7hdcdv 3f10.5)')(dcdv(k,ind1),k=1,3)
99 * Derivatives of SC(i+1) in theta(i+2)
101 xx1(1)=-0.5D0*xloc(2,i+1)
102 xx1(2)= 0.5D0*xloc(1,i+1)
106 xj=xj+r(j,k,i)*xx1(k)
113 rj=rj+prod(j,k,i)*xx(k)
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)
120 c write (iout,*) "dxdv",(dxdv(j,ind1),j=1,3)
122 * Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
123 * than the other off-diagonal derivatives.
125 if (i.lt.nres-2) then
129 dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
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)
137 c write(iout,*)ind1+1,(dxdv(j,ind1+1),j=1,3)
141 * Derivatives of DC(i+1) in phi(i+2)
148 dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
151 prodrt(j,k,i)=dp(j,k)
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)
158 * Derivatives of SC(i+1) in phi(i+2)
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)
167 rj=rj+prod(j,k,i)*xx(k)
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)
176 * Derivatives of SC(i+1) in phi(i+3).
182 dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
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)
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).
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)
204 tempkl=tempkl+prordt(k,m,i)*fromto(m,l)
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)
213 * Derivatives of virtual-bond vectors in theta
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)
218 c write(iout,'(7hdcdv 3f10.5)')(dcdv(k,ind1),k=1,3)
219 * Derivatives of SC vectors in theta
223 dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
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)
230 c write(iout,'(7htheta 3f10.5)')(dxdv(k,ind1),k=1,3)
233 *--- Calculate the derivatives in phi
239 tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
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)
252 dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
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)
263 if (nvar.le.nphi+ntheta) return
267 * Derivatives in alpha and omega:
270 if (iabs(itype(i)).eq.10 .or. itype(i).eq.ntyp1!) cycle
271 & .or. mask_side(i).eq.0 ) cycle
277 if(alphi.ne.alphi) alphi=100.0
278 if(omegi.ne.omegi) omegi=-100.0
283 cd print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
288 temp(1,1)=-dsci*sinalphi
289 temp(2,1)= dsci*cosalphi*cosomegi
290 temp(3,1)=-dsci*cosalphi*sinomegi
292 temp(2,2)=-dsci*sinalphi*sinomegi
293 temp(3,2)=-dsci*sinalphi*cosomegi
294 theta2=pi-0.5D0*theta(i+1)
298 cd print *,((temp(l,k),l=1,3),k=1,2)
302 xxp= xp*cost2+yp*sint2
303 yyp=-xp*sint2+yp*cost2
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)
311 dj=dj+prod(k,l,i-1)*xx(l)
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)
325 c-----------------------------------------------------------------------------
326 subroutine build_fromto(i,j,fromto)
329 include 'COMMON.CHAIN'
330 include 'COMMON.IOUNITS'
332 double precision fromto(3,3),temp(3,3),dp(3,3)
333 double precision dpkl
336 * generate the matrix products of type r(i)t(i)...r(j)t(j) on the fly
338 c write (iout,*) "temp on entry"
339 c write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
350 fromto(k,l)=temp(k,l)
360 dpkl=dpkl+temp(k,m)*rt(m,l,j-1)
372 c write (iout,*) "temp upon exit"
373 c write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)