1 subroutine intcartderiv
2 implicit real*8 (a-h,o-z)
11 include 'COMMON.INTERACT'
12 include 'COMMON.DERIV'
13 include 'COMMON.IOUNITS'
14 include 'COMMON.LOCAL'
15 include 'COMMON.SCCOR'
16 double precision dcostheta(3,2,maxres),
17 & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
18 & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
19 & dsinomega(3,3,maxres),vo1(3),vo2(3),vo3(3),
20 & dummy(3),vp1(3),vp2(3),vp3(3),vpp1(3),n(3)
22 #if defined(MPI) && defined(PARINTDER)
23 if (nfgtasks.gt.1 .and. me.eq.king)
24 & call MPI_Bcast(8,1,MPI_INTEGER,king,FG_COMM,IERROR)
29 c write (iout,*) "iphi1_start",iphi1_start," iphi1_end",iphi1_end
39 c Derivatives of theta's
40 #if defined(MPI) && defined(PARINTDER)
41 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
42 do i=max0(ithet_start-1,3),ithet_end
47 sint=sqrt(1-cost*cost)
49 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
51 if (itype(i-1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
52 dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
54 if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
57 #if defined(MPI) && defined(PARINTDER)
58 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
59 do i=max0(ithet_start-1,3),ithet_end
63 if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
64 cost1=dcos(omicron(1,i))
65 sint1=sqrt(1-cost1*cost1)
66 cost2=dcos(omicron(2,i))
67 sint2=sqrt(1-cost2*cost2)
69 CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1)
70 dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
71 & cost1*dc_norm(j,i-2))/
73 domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)
74 dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
75 & +cost1*(dc_norm(j,i-1+nres)))/
77 domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)
78 CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
79 CC Looks messy but better than if in loop
80 dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
81 & +cost2*dc_norm(j,i-1))/
83 domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
84 dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
85 & +cost2*(-dc_norm(j,i-1+nres)))/
87 c write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
88 domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
94 c If phi is 0 or 180 degrees, then the formulas
95 c have to be derived by power series expansion of the
96 c conventional formulas around 0 and 180.
98 do i=iphi1_start,iphi1_end
102 c if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
103 c & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
104 c the conventional case
106 sint1=dsin(theta(i-1))
109 cost1=dcos(theta(i-1))
111 scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
112 fac0=1.0d0/(sint1*sint)
115 fac3=cosg*cost1/(sint1*sint1)
116 fac4=cosg*cost/(sint*sint)
117 c Obtaining the gamma derivatives from sine derivative
118 if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
119 & phi(i).gt.pi34.and.phi(i).le.pi.or.
120 & phi(i).gt.-pi.and.phi(i).le.-pi34) then
121 call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
122 call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
123 call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
128 c if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
129 dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
130 & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
131 dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
133 & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
134 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
135 dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
136 dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
137 & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
138 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
139 dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
141 c Bug fixed 3/24/05 (AL)
143 c Obtaining the gamma derivatives from cosine derivative
146 c if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
147 dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
148 & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
149 & dc_norm(j,i-3))/vbld(i-2)
150 dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
151 dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
152 & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
154 dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
155 dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
156 & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
157 & dc_norm(j,i-1))/vbld(i)
158 dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
163 Calculate derivative of Tauangle
165 do i=itau_start,itau_end
169 if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
170 c if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
171 c & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
172 cc dtauangle(j,intertyp,dervityp,residue number)
173 cc INTERTYP=1 SC...Ca...Ca..Ca
174 c the conventional case
176 sint1=dsin(omicron(2,i-1))
177 sing=dsin(tauangle(1,i))
179 cost1=dcos(omicron(2,i-1))
180 cosg=dcos(tauangle(1,i))
182 dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
183 cc write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
185 scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
186 fac0=1.0d0/(sint1*sint)
189 fac3=cosg*cost1/(sint1*sint1)
190 fac4=cosg*cost/(sint*sint)
191 cc write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
192 c Obtaining the gamma derivatives from sine derivative
193 if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
194 & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
195 & tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
196 call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
197 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
198 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
203 dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
204 &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
205 & *vbld_inv(i-2+nres)
206 dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
208 & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
209 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
210 c write(iout,*) "dsintau", dsintau(j,1,2,i)
211 dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
212 c Bug fixed 3/24/05 (AL)
213 dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
214 & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
215 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
216 dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
218 c Obtaining the gamma derivatives from cosine derivative
221 dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
222 & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
223 & (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
224 dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
225 dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
226 & dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
228 dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
229 dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
230 & dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
231 & dc_norm(j,i-1))/vbld(i)
232 dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
233 c write (iout,*) "else",i
237 c write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)
240 CC Second case Ca...Ca...Ca...SC
242 do i=itau_start,itau_end
246 if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
247 & (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
248 c the conventional case
249 sint=dsin(omicron(1,i))
250 sint1=dsin(theta(i-1))
251 sing=dsin(tauangle(2,i))
252 cost=dcos(omicron(1,i))
253 cost1=dcos(theta(i-1))
254 cosg=dcos(tauangle(2,i))
256 c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
258 scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
259 fac0=1.0d0/(sint1*sint)
262 fac3=cosg*cost1/(sint1*sint1)
263 fac4=cosg*cost/(sint*sint)
264 c Obtaining the gamma derivatives from sine derivative
265 if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
266 & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
267 & tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
268 call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
269 call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
270 call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
275 dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
276 & +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
277 c write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
278 c &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
279 dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
281 & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
282 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
283 c write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
284 c & sing*ctgt*domicron(j,1,2,i),
285 c & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
286 dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
287 c Bug fixed 3/24/05 (AL)
288 dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
289 & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
290 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
291 dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
293 c Obtaining the gamma derivatives from cosine derivative
296 dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
297 & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
298 & dc_norm(j,i-3))/vbld(i-2)
299 dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
300 dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
301 & dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
302 & dcosomicron(j,1,1,i)
303 dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
304 dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
305 & dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
306 & dc_norm(j,i-1+nres))/vbld(i-1+nres)
307 dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
308 c write(iout,*) i,j,"else", dtauangle(j,2,3,i)
313 CCC third case SC...Ca...Ca...SC
316 do i=itau_start,itau_end
320 c the conventional case
321 if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
322 &(itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
323 sint=dsin(omicron(1,i))
324 sint1=dsin(omicron(2,i-1))
325 sing=dsin(tauangle(3,i))
326 cost=dcos(omicron(1,i))
327 cost1=dcos(omicron(2,i-1))
328 cosg=dcos(tauangle(3,i))
330 dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
331 c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
333 scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
334 fac0=1.0d0/(sint1*sint)
337 fac3=cosg*cost1/(sint1*sint1)
338 fac4=cosg*cost/(sint*sint)
339 c Obtaining the gamma derivatives from sine derivative
340 if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
341 & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
342 & tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
343 call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
344 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
345 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
350 dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
351 & -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
352 & *vbld_inv(i-2+nres)
353 dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
355 & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
356 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
357 dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
358 c Bug fixed 3/24/05 (AL)
359 dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
360 & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
361 & *vbld_inv(i-1+nres)
362 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
363 dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
365 c Obtaining the gamma derivatives from cosine derivative
368 dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
369 & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
370 & dc_norm2(j,i-2+nres))/vbld(i-2+nres)
371 dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
372 dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
373 & dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
374 & dcosomicron(j,1,1,i)
375 dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
376 dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
377 & dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
378 & dc_norm(j,i-1+nres))/vbld(i-1+nres)
379 dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
380 c write(iout,*) "else",i
386 c Derivatives of side-chain angles alpha and omega
387 #if defined(MPI) && defined(PARINTDER)
388 do i=ibond_start,ibond_end
392 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
393 fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
397 fac9=fac5/vbld(i+nres)
398 scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
399 scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
400 cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
401 & scalar(dC_norm(1,i),dC_norm(1,i+nres))
402 & -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
403 sina=sqrt(1-cosa*cosa)
405 c write (iout,*) "i",i," cosa",cosa," sina",sina," sino",sino
407 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
408 & dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
409 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
410 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
411 & scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
412 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
413 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
414 & dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
416 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
418 c obtaining the derivatives of omega from sines
419 if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
420 & omeg(i).gt.pi34.and.omeg(i).le.pi.or.
421 & omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
422 fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
424 fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
425 fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
426 call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
427 call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
428 call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
429 coso_inv=1.0d0/dcos(omeg(i))
431 dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
432 & +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
433 & sino*dc_norm(j,i-1))/vbld(i)
434 domega(j,1,i)=coso_inv*dsinomega(j,1,i)
435 dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
436 & +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
437 & -sino*dc_norm(j,i)/vbld(i+1)
438 domega(j,2,i)=coso_inv*dsinomega(j,2,i)
439 dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
440 & fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
442 domega(j,3,i)=coso_inv*dsinomega(j,3,i)
445 c obtaining the derivatives of omega from cosines
446 fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
447 fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
452 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
453 & dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
454 & (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
455 & fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
456 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
457 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
458 & dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
459 & dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
460 & (scala2-fac11*cosa)*(0.25d0*sina/fac10*
461 & dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
463 domega(j,2,i)=-1/sino*dcosomega(j,2,i)
464 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
465 & scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
466 & (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
467 domega(j,3,i)=-1/sino*dcosomega(j,3,i)
480 #if defined(MPI) && defined(PARINTDER)
481 if (nfgtasks.gt.1) then
483 cd write (iout,*) "Gather dtheta"
485 write (iout,*) "dtheta before gather"
487 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
490 call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
491 & MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
492 & king,FG_COMM,IERROR)
494 cd write (iout,*) "Gather dphi"
496 write (iout,*) "dphi before gather"
498 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
501 call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
502 & MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
503 & king,FG_COMM,IERROR)
504 cd write (iout,*) "Gather dalpha"
507 call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
508 & MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
509 & king,FG_COMM,IERROR)
510 cd write (iout,*) "Gather domega"
512 call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
513 & MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
514 & king,FG_COMM,IERROR)
519 write (iout,*) "dtheta after gather"
521 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),k=1,2)
523 write (iout,*) "dphi after gather"
525 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
527 write (iout,*) "dalpha after gather"
529 write (iout,'(i3,3(3f8.5,3x))') i,((dalpha(j,k,i),j=1,3),k=1,3)
531 write (iout,*) "domega after gather"
533 write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
539 subroutine checkintcartgrad
540 implicit real*8 (a-h,o-z)
545 include 'COMMON.CHAIN'
548 include 'COMMON.INTERACT'
549 include 'COMMON.DERIV'
550 include 'COMMON.IOUNITS'
551 include 'COMMON.SETUP'
552 double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
553 & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
554 double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
555 & omeg_s(maxres),dc_norm_s(3)
556 double precision aincr /1.0d-5/
564 c Check theta gradient
566 & "Analytical (upper) and numerical (lower) gradient of theta"
573 call int_from_cart1(.false.)
574 dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
577 dc(j,i-1)=dc(j,i-1)+aincr
579 dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
582 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
583 & (dtheta(j,2,i),j=1,3)
584 write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
585 & (dthetanum(j,2,i),j=1,3)
586 write (iout,'(5x,3f10.5,5x,3f10.5)')
587 & (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
588 & (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
591 c Check gamma gradient
593 & "Analytical (upper) and numerical (lower) gradient of gamma"
599 dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
604 dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
607 dc(j,i-1)=dc(j,i-1)+aincr
609 dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
612 write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
613 & (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
614 write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
615 & (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
616 write (iout,'(5x,3(3f10.5,5x))')
617 & (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
618 & (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
619 & (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
622 c Check alpha gradient
624 & "Analytical (upper) and numerical (lower) gradient of alpha"
626 if(itype(i).ne.10) then
631 dalphanum(j,1,i)=(alph(i)-alph_s(i))
637 dalphanum(j,2,i)=(alph(i)-alph_s(i))
641 dc(j,i+nres)=dc(j,i+nres)+aincr
643 dalphanum(j,3,i)=(alph(i)-alph_s(i))
648 write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
649 & (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
650 write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
651 & (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
652 write (iout,'(5x,3(3f10.5,5x))')
653 & (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
654 & (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
655 & (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
658 c Check omega gradient
660 & "Analytical (upper) and numerical (lower) gradient of omega"
662 if(itype(i).ne.10) then
667 domeganum(j,1,i)=(omeg(i)-omeg_s(i))
673 domeganum(j,2,i)=(omeg(i)-omeg_s(i))
677 dc(j,i+nres)=dc(j,i+nres)+aincr
679 domeganum(j,3,i)=(omeg(i)-omeg_s(i))
684 write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
685 & (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
686 write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
687 & (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
688 write (iout,'(5x,3(3f10.5,5x))')
689 & (domeganum(j,1,i)/domega(j,1,i),j=1,3),
690 & (domeganum(j,2,i)/domega(j,2,i),j=1,3),
691 & (domeganum(j,3,i)/domega(j,3,i),j=1,3)
696 c------------------------------------------------------------
697 subroutine chainbuild_cart
698 implicit real*8 (a-h,o-z)
703 include 'COMMON.SETUP'
704 include 'COMMON.CHAIN'
705 include 'COMMON.LOCAL'
706 include 'COMMON.TIME1'
707 include 'COMMON.IOUNITS'
710 if (nfgtasks.gt.1) then
711 c write (iout,*) "BCAST in chainbuild_cart"
713 c Broadcast the order to build the chain and compute internal coordinates
714 c to the slaves. The slaves receive the order in ERGASTULUM.
716 c write (iout,*) "CHAINBUILD_CART: DC before BCAST"
718 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
719 c & (dc(j,i+nres),j=1,3)
722 & call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
723 time_bcast7=time_bcast7+MPI_Wtime()-time00
725 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
727 c write (iout,*) "CHAINBUILD_CART: DC after BCAST"
729 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
730 c & (dc(j,i+nres),j=1,3)
732 c write (iout,*) "End BCAST in chainbuild_cart"
734 time_bcast=time_bcast+MPI_Wtime()-time00
735 time_bcastc=time_bcastc+MPI_Wtime()-time01
744 c(j,i)=c(j,i-1)+dc(j,i-1)
749 c(j,i+nres)=c(j,i)+dc(j,i+nres)
752 c write (iout,*) "CHAINBUILD_CART"
754 call int_from_cart1(.false.)