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 c if (itype(i-1).ne.ntyp1)
52 dtheta(j,1,i)=-dcostheta(j,1,i)/sint
53 dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
55 c if (itype(i-1).ne.ntyp1)
56 dtheta(j,2,i)=-dcostheta(j,2,i)/sint
59 #if defined(MPI) && defined(PARINTDER)
60 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
61 do i=max0(ithet_start-1,3),ithet_end
65 if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then
66 cost1=dcos(omicron(1,i))
67 sint1=sqrt(1-cost1*cost1)
68 cost2=dcos(omicron(2,i))
69 sint2=sqrt(1-cost2*cost2)
71 CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1)
72 dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+
73 & cost1*dc_norm(j,i-2))/
75 domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)
76 dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
77 & +cost1*(dc_norm(j,i-1+nres)))/
79 domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)
80 CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
81 CC Looks messy but better than if in loop
82 dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
83 & +cost2*dc_norm(j,i-1))/
85 domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
86 dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
87 & +cost2*(-dc_norm(j,i-1+nres)))/
89 c write(iout,*) "vbld", i,itype(i),vbld(i-1+nres)
90 domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
96 c If phi is 0 or 180 degrees, then the formulas
97 c have to be derived by power series expansion of the
98 c conventional formulas around 0 and 180.
100 do i=iphi1_start,iphi1_end
104 c if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
105 c & .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
106 c the conventional case
108 sint1=dsin(theta(i-1))
111 cost1=dcos(theta(i-1))
113 scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
114 fac0=1.0d0/(sint1*sint)
117 fac3=cosg*cost1/(sint1*sint1)
118 fac4=cosg*cost/(sint*sint)
119 c Obtaining the gamma derivatives from sine derivative
120 if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
121 & phi(i).gt.pi34.and.phi(i).le.pi.or.
122 & phi(i).gt.-pi.and.phi(i).le.-pi34) then
123 call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
124 call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
125 call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
130 c if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
131 dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
132 & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
133 dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
135 & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
136 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
137 dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
138 dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
139 & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
140 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
141 dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
143 c Bug fixed 3/24/05 (AL)
145 c Obtaining the gamma derivatives from cosine derivative
148 c if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then
149 dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
150 & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
151 & dc_norm(j,i-3))/vbld(i-2)
152 dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
153 dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
154 & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
156 dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
157 dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
158 & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
159 & dc_norm(j,i-1))/vbld(i)
160 dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
165 Calculate derivative of Tauangle
167 do i=itau_start,itau_end
171 if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
172 c if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or.
173 c & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
174 cc dtauangle(j,intertyp,dervityp,residue number)
175 cc INTERTYP=1 SC...Ca...Ca..Ca
176 c the conventional case
178 sint1=dsin(omicron(2,i-1))
179 sing=dsin(tauangle(1,i))
181 cost1=dcos(omicron(2,i-1))
182 cosg=dcos(tauangle(1,i))
184 dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
185 cc write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
187 scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
188 fac0=1.0d0/(sint1*sint)
191 fac3=cosg*cost1/(sint1*sint1)
192 fac4=cosg*cost/(sint*sint)
193 cc write(iout,*) "faki",fac0,fac1,fac2,fac3,fac4
194 c Obtaining the gamma derivatives from sine derivative
195 if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
196 & tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
197 & tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
198 call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
199 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
200 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
205 dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
206 &-(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres)))
207 & *vbld_inv(i-2+nres)
208 dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
210 & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i))
211 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
212 c write(iout,*) "dsintau", dsintau(j,1,2,i)
213 dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
214 c Bug fixed 3/24/05 (AL)
215 dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
216 & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
217 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
218 dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
220 c Obtaining the gamma derivatives from cosine derivative
223 dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
224 & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp*
225 & (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
226 dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
227 dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
228 & dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
230 dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
231 dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
232 & dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
233 & dc_norm(j,i-1))/vbld(i)
234 dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
235 c write (iout,*) "else",i
239 c write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)
242 CC Second case Ca...Ca...Ca...SC
244 do i=itau_start,itau_end
248 if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
249 & (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle
250 c the conventional case
251 sint=dsin(omicron(1,i))
252 sint1=dsin(theta(i-1))
253 sing=dsin(tauangle(2,i))
254 cost=dcos(omicron(1,i))
255 cost1=dcos(theta(i-1))
256 cosg=dcos(tauangle(2,i))
258 c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
260 scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
261 fac0=1.0d0/(sint1*sint)
264 fac3=cosg*cost1/(sint1*sint1)
265 fac4=cosg*cost/(sint*sint)
266 c Obtaining the gamma derivatives from sine derivative
267 if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
268 & tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
269 & tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
270 call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
271 call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
272 call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
277 dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
278 & +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
279 c write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
280 c &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
281 dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
283 & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
284 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
285 c write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
286 c & sing*ctgt*domicron(j,1,2,i),
287 c & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
288 dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
289 c Bug fixed 3/24/05 (AL)
290 dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
291 & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
292 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
293 dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
295 c Obtaining the gamma derivatives from cosine derivative
298 dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
299 & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
300 & dc_norm(j,i-3))/vbld(i-2)
301 dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
302 dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
303 & dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
304 & dcosomicron(j,1,1,i)
305 dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
306 dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
307 & dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
308 & dc_norm(j,i-1+nres))/vbld(i-1+nres)
309 dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
310 c write(iout,*) i,j,"else", dtauangle(j,2,3,i)
315 CCC third case SC...Ca...Ca...SC
318 do i=itau_start,itau_end
322 c the conventional case
323 if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or.
324 &(itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle
325 sint=dsin(omicron(1,i))
326 sint1=dsin(omicron(2,i-1))
327 sing=dsin(tauangle(3,i))
328 cost=dcos(omicron(1,i))
329 cost1=dcos(omicron(2,i-1))
330 cosg=dcos(tauangle(3,i))
332 dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
333 c dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
335 scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
336 fac0=1.0d0/(sint1*sint)
339 fac3=cosg*cost1/(sint1*sint1)
340 fac4=cosg*cost/(sint*sint)
341 c Obtaining the gamma derivatives from sine derivative
342 if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
343 & tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
344 & tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
345 call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
346 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
347 call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
352 dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1)
353 & -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres))
354 & *vbld_inv(i-2+nres)
355 dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
357 & -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i))
358 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
359 dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
360 c Bug fixed 3/24/05 (AL)
361 dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
362 & +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))
363 & *vbld_inv(i-1+nres)
364 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
365 dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
367 c Obtaining the gamma derivatives from cosine derivative
370 dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3*
371 & dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp*
372 & dc_norm2(j,i-2+nres))/vbld(i-2+nres)
373 dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
374 dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2*
375 & dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4*
376 & dcosomicron(j,1,1,i)
377 dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
378 dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
379 & dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
380 & dc_norm(j,i-1+nres))/vbld(i-1+nres)
381 dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
382 c write(iout,*) "else",i
388 c Derivatives of side-chain angles alpha and omega
389 #if defined(MPI) && defined(PARINTDER)
390 do i=ibond_start,ibond_end
394 if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
395 fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
399 fac9=fac5/vbld(i+nres)
400 scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
401 scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
402 cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
403 & scalar(dC_norm(1,i),dC_norm(1,i+nres))
404 & -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
405 sina=sqrt(1-cosa*cosa)
407 c write (iout,*) "i",i," cosa",cosa," sina",sina," sino",sino
409 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
410 & dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
411 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
412 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
413 & scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
414 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
415 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
416 & dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
418 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
420 c obtaining the derivatives of omega from sines
421 if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
422 & omeg(i).gt.pi34.and.omeg(i).le.pi.or.
423 & omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
424 fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
426 fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
427 fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
428 call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
429 call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
430 call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
431 coso_inv=1.0d0/dcos(omeg(i))
433 dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
434 & +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
435 & sino*dc_norm(j,i-1))/vbld(i)
436 domega(j,1,i)=coso_inv*dsinomega(j,1,i)
437 dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
438 & +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
439 & -sino*dc_norm(j,i)/vbld(i+1)
440 domega(j,2,i)=coso_inv*dsinomega(j,2,i)
441 dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
442 & fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
444 domega(j,3,i)=coso_inv*dsinomega(j,3,i)
447 c obtaining the derivatives of omega from cosines
448 fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
449 fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
454 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
455 & dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
456 & (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
457 & fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
458 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
459 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
460 & dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
461 & dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
462 & (scala2-fac11*cosa)*(0.25d0*sina/fac10*
463 & dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
465 domega(j,2,i)=-1/sino*dcosomega(j,2,i)
466 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
467 & scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
468 & (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
469 domega(j,3,i)=-1/sino*dcosomega(j,3,i)
482 #if defined(MPI) && defined(PARINTDER)
483 if (nfgtasks.gt.1) then
485 cd write (iout,*) "Gather dtheta"
487 write (iout,*) "dtheta before gather"
489 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
492 call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
493 & MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
494 & king,FG_COMM,IERROR)
496 cd write (iout,*) "Gather dphi"
498 write (iout,*) "dphi before gather"
500 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
503 call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
504 & MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
505 & king,FG_COMM,IERROR)
506 cd write (iout,*) "Gather dalpha"
509 call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
510 & MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
511 & king,FG_COMM,IERROR)
512 cd write (iout,*) "Gather domega"
514 call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
515 & MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
516 & king,FG_COMM,IERROR)
521 write (iout,*) "dtheta after gather"
523 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),k=1,2)
525 write (iout,*) "dphi after gather"
527 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
529 write (iout,*) "dalpha after gather"
531 write (iout,'(i3,3(3f8.5,3x))') i,((dalpha(j,k,i),j=1,3),k=1,3)
533 write (iout,*) "domega after gather"
535 write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
541 subroutine checkintcartgrad
542 implicit real*8 (a-h,o-z)
547 include 'COMMON.CHAIN'
550 include 'COMMON.INTERACT'
551 include 'COMMON.DERIV'
552 include 'COMMON.IOUNITS'
553 include 'COMMON.SETUP'
554 double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
555 & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
556 double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
557 & omeg_s(maxres),dc_norm_s(3)
558 double precision aincr /1.0d-5/
566 c Check theta gradient
568 & "Analytical (upper) and numerical (lower) gradient of theta"
575 call int_from_cart1(.false.)
576 dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
579 dc(j,i-1)=dc(j,i-1)+aincr
581 dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
584 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
585 & (dtheta(j,2,i),j=1,3)
586 write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
587 & (dthetanum(j,2,i),j=1,3)
588 write (iout,'(5x,3f10.5,5x,3f10.5)')
589 & (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
590 & (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
593 c Check gamma gradient
595 & "Analytical (upper) and numerical (lower) gradient of gamma"
601 dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
606 dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
609 dc(j,i-1)=dc(j,i-1)+aincr
611 dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
614 write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
615 & (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
616 write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
617 & (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
618 write (iout,'(5x,3(3f10.5,5x))')
619 & (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
620 & (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
621 & (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
624 c Check alpha gradient
626 & "Analytical (upper) and numerical (lower) gradient of alpha"
628 if(itype(i).ne.10) then
633 dalphanum(j,1,i)=(alph(i)-alph_s(i))
639 dalphanum(j,2,i)=(alph(i)-alph_s(i))
643 dc(j,i+nres)=dc(j,i+nres)+aincr
645 dalphanum(j,3,i)=(alph(i)-alph_s(i))
650 write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
651 & (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
652 write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
653 & (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
654 write (iout,'(5x,3(3f10.5,5x))')
655 & (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
656 & (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
657 & (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
660 c Check omega gradient
662 & "Analytical (upper) and numerical (lower) gradient of omega"
664 if(itype(i).ne.10) then
669 domeganum(j,1,i)=(omeg(i)-omeg_s(i))
675 domeganum(j,2,i)=(omeg(i)-omeg_s(i))
679 dc(j,i+nres)=dc(j,i+nres)+aincr
681 domeganum(j,3,i)=(omeg(i)-omeg_s(i))
686 write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
687 & (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
688 write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
689 & (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
690 write (iout,'(5x,3(3f10.5,5x))')
691 & (domeganum(j,1,i)/domega(j,1,i),j=1,3),
692 & (domeganum(j,2,i)/domega(j,2,i),j=1,3),
693 & (domeganum(j,3,i)/domega(j,3,i),j=1,3)
698 c------------------------------------------------------------
699 subroutine chainbuild_cart
700 implicit real*8 (a-h,o-z)
705 include 'COMMON.SETUP'
706 include 'COMMON.CHAIN'
707 include 'COMMON.LOCAL'
708 include 'COMMON.TIME1'
709 include 'COMMON.IOUNITS'
712 if (nfgtasks.gt.1) then
713 c write (iout,*) "BCAST in chainbuild_cart"
715 c Broadcast the order to build the chain and compute internal coordinates
716 c to the slaves. The slaves receive the order in ERGASTULUM.
718 c write (iout,*) "CHAINBUILD_CART: DC before BCAST"
720 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
721 c & (dc(j,i+nres),j=1,3)
724 & call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
725 time_bcast7=time_bcast7+MPI_Wtime()-time00
727 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
729 c write (iout,*) "CHAINBUILD_CART: DC after BCAST"
731 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
732 c & (dc(j,i+nres),j=1,3)
734 c write (iout,*) "End BCAST in chainbuild_cart"
736 time_bcast=time_bcast+MPI_Wtime()-time00
737 time_bcastc=time_bcastc+MPI_Wtime()-time01
746 c(j,i)=c(j,i-1)+dc(j,i-1)
751 c(j,i+nres)=c(j,i)+dc(j,i+nres)
755 c write (iout,*) "CHAINBUILD_CART"
757 call int_from_cart1(.false.)