Adam's changes
[unres.git] / source / wham / src-M-SAXS.safe / cartder.f
1       subroutine cartder 
2       implicit real*8 (a-h,o-z)
3 ***********************************************************************
4 * This subroutine calculates the derivatives of the consecutive virtual
5 * bond vectors and the SC vectors in the virtual-bond angles theta and
6 * virtual-torsional angles phi, as well as the derivatives of SC vectors
7 * in the angles alpha and omega, describing the location of a side chain
8 * in its local coordinate system.
9 *
10 * The derivatives are stored in the following arrays:
11 *
12 * DDCDV - the derivatives of virtual-bond vectors DC in theta and phi.
13 * The structure is as follows:
14 *
15 * dDC(x,2)/dT(3),...,dDC(z,2)/dT(3),0,             0,             0
16 * 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)
17 *         . . . . . . . . . . . .  . . . . . .
18 * 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)
19 *                          .
20 *                          .
21 *                          .
22 * 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)
23 *
24 * DXDV - the derivatives of the side-chain vectors in theta and phi. 
25 * The structure is same as above.
26 *
27 * DCDS - the derivatives of the side chain vectors in the local spherical
28 * andgles alph and omega:
29 *
30 * 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)
31 * 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)
32 *                          .
33 *                          .
34 *                          .
35 * 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)
36 *
37 * Version of March '95, based on an early version of November '91.
38 *
39 *********************************************************************** 
40       include 'DIMENSIONS'
41       include 'DIMENSIONS.ZSCOPT'
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       dimension drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),temp(3,3),
49      &     fromto(3,3,maxdim),prordt(3,3,maxres),prodrt(3,3,maxres)
50       dimension xx(3),xx1(3)
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         do j=1,3
126           do k=1,2
127             dpjk=0.0D0
128             do l=1,3
129               dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
130             enddo
131             dp(j,k)=dpjk
132             prordt(j,k,i)=dp(j,k)
133           enddo
134           dp(j,3)=0.0D0
135           dcdv(j,ind1)=vbl*dp(j,1)       
136         enddo
137 *
138 * Derivatives of SC(i+1) in theta(i+2)
139
140         xx1(1)=-0.5D0*xloc(2,i+1)
141         xx1(2)= 0.5D0*xloc(1,i+1)
142         do j=1,3
143           xj=0.0D0
144           do k=1,2
145             xj=xj+r(j,k,i)*xx1(k)
146           enddo
147           xx(j)=xj
148         enddo
149         do j=1,3
150           rj=0.0D0
151           do k=1,3
152             rj=rj+prod(j,k,i)*xx(k)
153           enddo
154           dxdv(j,ind1)=rj
155         enddo
156 *
157 * Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
158 * than the other off-diagonal derivatives.
159 *
160         do j=1,3
161           dxoiij=0.0D0
162           do k=1,3
163             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
164           enddo
165           dxdv(j,ind1+1)=dxoiij
166         enddo
167 cd      print *,ind1+1,(dxdv(j,ind1+1),j=1,3)
168 *
169 * Derivatives of DC(i+1) in phi(i+2)
170 *
171         do j=1,3
172           do k=1,3
173             dpjk=0.0
174             do l=2,3
175               dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
176             enddo
177             dp(j,k)=dpjk
178             prodrt(j,k,i)=dp(j,k)
179           enddo 
180           dcdv(j+3,ind1)=vbl*dp(j,1)
181         enddo
182 *
183 * Derivatives of SC(i+1) in phi(i+2)
184 *
185         xx(1)= 0.0D0 
186         xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
187         xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
188         do j=1,3
189           rj=0.0D0
190           do k=2,3
191             rj=rj+prod(j,k,i)*xx(k)
192           enddo
193           dxdv(j+3,ind1)=-rj
194         enddo
195 *
196 * Derivatives of SC(i+1) in phi(i+3).
197 *
198         do j=1,3
199           dxoiij=0.0D0
200           do k=1,3
201             dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
202           enddo
203           dxdv(j+3,ind1+1)=dxoiij
204         enddo
205 *
206 * Calculate the derivatives of DC(i+1) and SC(i+1) in theta(i+3) thru 
207 * theta(nres) and phi(i+3) thru phi(nres).
208 *
209         do j=i+1,nres-2
210           ind1=ind1+1
211           ind=indmat(i+1,j+1)
212 cd        print *,'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
213           do k=1,3
214             do l=1,3
215               tempkl=0.0D0
216               do m=1,2
217                 tempkl=tempkl+prordt(k,m,i)*fromto(m,l,ind)
218               enddo
219               temp(k,l)=tempkl
220             enddo
221           enddo  
222 cd        print '(9f8.3)',((fromto(k,l,ind),l=1,3),k=1,3)
223 cd        print '(9f8.3)',((prod(k,l,i),l=1,3),k=1,3)
224 cd        print '(9f8.3)',((temp(k,l),l=1,3),k=1,3)
225 * Derivatives of virtual-bond vectors in theta
226           do k=1,3
227             dcdv(k,ind1)=vbl*temp(k,1)
228           enddo
229 cd        print '(3f8.3)',(dcdv(k,ind1),k=1,3)
230 * Derivatives of SC vectors in theta
231           do k=1,3
232             dxoijk=0.0D0
233             do l=1,3
234               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
235             enddo
236             dxdv(k,ind1+1)=dxoijk
237           enddo
238 *
239 *--- Calculate the derivatives in phi
240 *
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           do k=1,3
251             dcdv(k+3,ind1)=vbl*temp(k,1)
252           enddo
253           do k=1,3
254             dxoijk=0.0D0
255             do l=1,3
256               dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
257             enddo
258             dxdv(k+3,ind1+1)=dxoijk
259           enddo
260         enddo
261       enddo
262 *
263 * Derivatives in alpha and omega:
264 *
265       do i=2,nres-1
266         dsci=dsc(iabs(itype(i)))
267         alphi=alph(i)
268         omegi=omeg(i)
269 cd      print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
270         cosalphi=dcos(alphi)
271         sinalphi=dsin(alphi)
272         cosomegi=dcos(omegi)
273         sinomegi=dsin(omegi)
274         temp(1,1)=-dsci*sinalphi
275         temp(2,1)= dsci*cosalphi*cosomegi
276         temp(3,1)=-dsci*cosalphi*sinomegi
277         temp(1,2)=0.0D0
278         temp(2,2)=-dsci*sinalphi*sinomegi
279         temp(3,2)=-dsci*sinalphi*cosomegi
280         theta2=pi-0.5D0*theta(i+1)
281         cost2=dcos(theta2)
282         sint2=dsin(theta2)
283         jjj=0
284 cd      print *,((temp(l,k),l=1,3),k=1,2)
285         do j=1,2
286           xp=temp(1,j)
287           yp=temp(2,j)
288           xxp= xp*cost2+yp*sint2
289           yyp=-xp*sint2+yp*cost2
290           zzp=temp(3,j)
291           xx(1)=xxp
292           xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
293           xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
294           do k=1,3
295             dj=0.0D0
296             do l=1,3
297               dj=dj+prod(k,l,i-1)*xx(l)
298             enddo
299             dxds(jjj+k,i)=dj
300           enddo
301           jjj=jjj+3
302         enddo
303       enddo
304       return
305       end
306