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 double precision dcostheta(3,2,maxres),
16 & dcosphi(3,3,maxres),dsinphi(3,3,maxres),
17 & dcosalpha(3,3,maxres),dcosomega(3,3,maxres),
18 & dsinomega(3,3,maxres),vo1(3),vo2(3),vo3(3),
19 & dummy(3),vp1(3),vp2(3),vp3(3),vpp1(3),n(3)
21 #if defined(MPI) && defined(PARINTDER)
22 if (nfgtasks.gt.1 .and. me.eq.king)
23 & call MPI_Bcast(8,1,MPI_INTEGER,king,FG_COMM,IERROR)
28 c write (iout,*) "iphi1_start",iphi1_start," iphi1_end",iphi1_end
38 c Derivatives of theta's
39 #if defined(MPI) && defined(PARINTDER)
40 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
41 do i=max0(ithet_start-1,3),ithet_end
46 sint=sqrt(1-cost*cost)
48 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
50 if (itype(i-1).ne.21) dtheta(j,1,i)=-dcostheta(j,1,i)/sint
51 dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
53 if (itype(i-1).ne.21) dtheta(j,2,i)=-dcostheta(j,2,i)/sint
58 c If phi is 0 or 180 degrees, then the formulas
59 c have to be derived by power series expansion of the
60 c conventional formulas around 0 and 180.
62 do i=iphi1_start,iphi1_end
66 c if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle
67 c the conventional case
69 sint1=dsin(theta(i-1))
72 cost1=dcos(theta(i-1))
74 scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
75 fac0=1.0d0/(sint1*sint)
78 fac3=cosg*cost1/(sint1*sint1)
79 fac4=cosg*cost/(sint*sint)
80 c Obtaining the gamma derivatives from sine derivative
81 if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
82 & phi(i).gt.pi34.and.phi(i).le.pi.or.
83 & phi(i).gt.-pi.and.phi(i).le.-pi34) then
84 call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
85 call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
86 call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
91 if (itype(i-1).ne.21 .and. itype(i-2).ne.21) then
92 dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
93 & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
94 dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
96 & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
97 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
98 dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
99 dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
100 & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
101 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
102 dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
104 c Bug fixed 3/24/05 (AL)
106 c Obtaining the gamma derivatives from cosine derivative
109 if (itype(i-1).ne.21 .and. itype(i-2).ne.21) then
110 dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
111 & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
112 & dc_norm(j,i-3))/vbld(i-2)
113 dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
114 dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
115 & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
117 dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
118 dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
119 & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
120 & dc_norm(j,i-1))/vbld(i)
121 dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
127 c Derivatives of side-chain angles alpha and omega
128 #if defined(MPI) && defined(PARINTDER)
129 do i=ibond_start,ibond_end
133 if(itype(i).ne.10 .and. itype(i).ne.21) then
134 fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
138 fac9=fac5/vbld(i+nres)
139 scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
140 scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
141 cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
142 & scalar(dC_norm(1,i),dC_norm(1,i+nres))
143 & -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
144 sina=sqrt(1-cosa*cosa)
146 c write (iout,*) "i",i," cosa",cosa," sina",sina," sino",sino
148 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
149 & dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
150 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
151 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
152 & scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
153 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
154 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
155 & dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
157 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
159 c obtaining the derivatives of omega from sines
160 if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
161 & omeg(i).gt.pi34.and.omeg(i).le.pi.or.
162 & omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
163 fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
165 fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
166 fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
167 call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
168 call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
169 call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
170 coso_inv=1.0d0/dcos(omeg(i))
172 dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
173 & +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
174 & sino*dc_norm(j,i-1))/vbld(i)
175 domega(j,1,i)=coso_inv*dsinomega(j,1,i)
176 dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
177 & +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
178 & -sino*dc_norm(j,i)/vbld(i+1)
179 domega(j,2,i)=coso_inv*dsinomega(j,2,i)
180 dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
181 & fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
183 domega(j,3,i)=coso_inv*dsinomega(j,3,i)
186 c obtaining the derivatives of omega from cosines
187 fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
188 fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
193 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
194 & dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
195 & (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
196 & fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
197 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
198 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
199 & dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
200 & dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
201 & (scala2-fac11*cosa)*(0.25d0*sina/fac10*
202 & dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
204 domega(j,2,i)=-1/sino*dcosomega(j,2,i)
205 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
206 & scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
207 & (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
208 domega(j,3,i)=-1/sino*dcosomega(j,3,i)
221 #if defined(MPI) && defined(PARINTDER)
222 if (nfgtasks.gt.1) then
224 cd write (iout,*) "Gather dtheta"
226 write (iout,*) "dtheta before gather"
228 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
231 call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
232 & MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
233 & king,FG_COMM,IERROR)
235 cd write (iout,*) "Gather dphi"
237 write (iout,*) "dphi before gather"
239 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
242 call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
243 & MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
244 & king,FG_COMM,IERROR)
245 cd write (iout,*) "Gather dalpha"
248 call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
249 & MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
250 & king,FG_COMM,IERROR)
251 cd write (iout,*) "Gather domega"
253 call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
254 & MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
255 & king,FG_COMM,IERROR)
260 write (iout,*) "dtheta after gather"
262 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),k=1,2)
264 write (iout,*) "dphi after gather"
266 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
268 write (iout,*) "dalpha after gather"
270 write (iout,'(i3,3(3f8.5,3x))') i,((dalpha(j,k,i),j=1,3),k=1,3)
272 write (iout,*) "domega after gather"
274 write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
280 subroutine checkintcartgrad
281 implicit real*8 (a-h,o-z)
286 include 'COMMON.CHAIN'
289 include 'COMMON.INTERACT'
290 include 'COMMON.DERIV'
291 include 'COMMON.IOUNITS'
292 include 'COMMON.SETUP'
293 double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
294 & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
295 double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
296 & omeg_s(maxres),dc_norm_s(3)
297 double precision aincr /1.0d-5/
305 c Check theta gradient
307 & "Analytical (upper) and numerical (lower) gradient of theta"
314 call int_from_cart1(.false.)
315 dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
318 dc(j,i-1)=dc(j,i-1)+aincr
320 dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
323 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
324 & (dtheta(j,2,i),j=1,3)
325 write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
326 & (dthetanum(j,2,i),j=1,3)
327 write (iout,'(5x,3f10.5,5x,3f10.5)')
328 & (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
329 & (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
332 c Check gamma gradient
334 & "Analytical (upper) and numerical (lower) gradient of gamma"
340 dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
345 dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
348 dc(j,i-1)=dc(j,i-1)+aincr
350 dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
353 write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
354 & (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
355 write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
356 & (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
357 write (iout,'(5x,3(3f10.5,5x))')
358 & (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
359 & (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
360 & (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
363 c Check alpha gradient
365 & "Analytical (upper) and numerical (lower) gradient of alpha"
367 if(itype(i).ne.10) then
372 dalphanum(j,1,i)=(alph(i)-alph_s(i))
378 dalphanum(j,2,i)=(alph(i)-alph_s(i))
382 dc(j,i+nres)=dc(j,i+nres)+aincr
384 dalphanum(j,3,i)=(alph(i)-alph_s(i))
389 write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
390 & (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
391 write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
392 & (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
393 write (iout,'(5x,3(3f10.5,5x))')
394 & (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
395 & (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
396 & (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
399 c Check omega gradient
401 & "Analytical (upper) and numerical (lower) gradient of omega"
403 if(itype(i).ne.10) then
408 domeganum(j,1,i)=(omeg(i)-omeg_s(i))
414 domeganum(j,2,i)=(omeg(i)-omeg_s(i))
418 dc(j,i+nres)=dc(j,i+nres)+aincr
420 domeganum(j,3,i)=(omeg(i)-omeg_s(i))
425 write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
426 & (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
427 write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
428 & (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
429 write (iout,'(5x,3(3f10.5,5x))')
430 & (domeganum(j,1,i)/domega(j,1,i),j=1,3),
431 & (domeganum(j,2,i)/domega(j,2,i),j=1,3),
432 & (domeganum(j,3,i)/domega(j,3,i),j=1,3)
438 subroutine chainbuild_cart
439 implicit real*8 (a-h,o-z)
444 include 'COMMON.SETUP'
445 include 'COMMON.CHAIN'
446 include 'COMMON.LOCAL'
447 include 'COMMON.TIME1'
448 include 'COMMON.IOUNITS'
451 if (nfgtasks.gt.1) then
452 c write (iout,*) "BCAST in chainbuild_cart"
454 c Broadcast the order to build the chain and compute internal coordinates
455 c to the slaves. The slaves receive the order in ERGASTULUM.
457 c write (iout,*) "CHAINBUILD_CART: DC before BCAST"
459 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
460 c & (dc(j,i+nres),j=1,3)
463 & call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
464 time_bcast7=time_bcast7+MPI_Wtime()-time00
466 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
468 c write (iout,*) "CHAINBUILD_CART: DC after BCAST"
470 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
471 c & (dc(j,i+nres),j=1,3)
473 c write (iout,*) "End BCAST in chainbuild_cart"
475 time_bcast=time_bcast+MPI_Wtime()-time00
476 time_bcastc=time_bcastc+MPI_Wtime()-time01
484 c(j,i)=c(j,i-1)+dc(j,i-1)
489 c(j,i+nres)=c(j,i)+dc(j,i+nres)
492 c write (iout,*) "CHAINBUILD_CART"
494 call int_from_cart1(.false.)