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
29 c Derivatives of theta's
30 #if defined(MPI) && defined(PARINTDER)
31 c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
32 do i=max0(ithet_start-1,3),ithet_end
37 sint=sqrt(1-cost*cost)
39 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/
41 dtheta(j,1,i)=-1/sint*dcostheta(j,1,i)
42 dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/
44 dtheta(j,2,i)=-1/sint*dcostheta(j,2,i)
49 c If phi is 0 or 180 degrees, then the formulas
50 c have to be derived by power series expansion of the
51 c conventional formulas around 0 and 180.
53 do i=iphi1_start,iphi1_end
57 c the conventional case
59 sint1=dsin(theta(i-1))
62 cost1=dcos(theta(i-1))
64 scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
65 fac0=1.0d0/(sint1*sint)
68 fac3=cosg*cost1/(sint1*sint1)
69 fac4=cosg*cost/(sint*sint)
70 c Obtaining the gamma derivatives from sine derivative
71 if (phi(i).gt.-pi4.and.phi(i).le.pi4.or.
72 & phi(i).gt.pi34.and.phi(i).le.pi.or.
73 & phi(i).gt.-pi.and.phi(i).le.-pi34) then
74 call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
75 call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
76 call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
81 dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
82 & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
83 dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
85 & -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i))
86 & -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
87 dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
88 c Bug fixed 3/24/05 (AL)
89 dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i)
90 & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
91 c & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
92 dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
94 c Obtaining the gamma derivatives from cosine derivative
97 dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
98 & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
99 & dc_norm(j,i-3))/vbld(i-2)
100 dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
101 dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
102 & dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
104 dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
105 dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4*
106 & dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp*
107 & dc_norm(j,i-1))/vbld(i)
108 dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
113 c Derivatives of side-chain angles alpha and omega
114 #if defined(MPI) && defined(PARINTDER)
115 do i=ibond_start,ibond_end
119 if(itype(i).ne.10) then
120 fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
124 fac9=fac5/vbld(i+nres)
125 scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
126 scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
127 cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))*(
128 & scalar(dC_norm(1,i),dC_norm(1,i+nres))
129 & -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
130 sina=sqrt(1-cosa*cosa)
133 dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)-
134 & dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
135 dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
136 dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)-
137 & scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
138 dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
139 dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)-
140 & dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/
142 dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
144 c obtaining the derivatives of omega from sines
145 if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or.
146 & omeg(i).gt.pi34.and.omeg(i).le.pi.or.
147 & omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
148 fac15=dcos(theta(i+1))/(dsin(theta(i+1))*
150 fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
151 fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
152 call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
153 call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
154 call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
155 coso_inv=1.0d0/dcos(omeg(i))
157 dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1)
158 & +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)-(
159 & sino*dc_norm(j,i-1))/vbld(i)
160 domega(j,1,i)=coso_inv*dsinomega(j,1,i)
161 dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1)
162 & +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j)
163 & -sino*dc_norm(j,i)/vbld(i+1)
164 domega(j,2,i)=coso_inv*dsinomega(j,2,i)
165 dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)-
166 & fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/
168 domega(j,3,i)=coso_inv*dsinomega(j,3,i)
171 c obtaining the derivatives of omega from cosines
172 fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
173 fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
178 dcosomega(j,1,i)=(-(0.25d0*cosa/fac11*
179 & dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+
180 & (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina*
181 & fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
182 domega(j,1,i)=-1/sino*dcosomega(j,1,i)
183 dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2*
184 & dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11*
185 & dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+
186 & (scala2-fac11*cosa)*(0.25d0*sina/fac10*
187 & dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)
189 domega(j,2,i)=-1/sino*dcosomega(j,2,i)
190 dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)-
191 & scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+
192 & (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
193 domega(j,3,i)=-1/sino*dcosomega(j,3,i)
199 #if defined(MPI) && defined(PARINTDER)
200 if (nfgtasks.gt.1) then
202 cd write (iout,*) "Gather dtheta"
204 write (iout,*) "dtheta before gather"
206 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
209 call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),
210 & MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,
211 & king,FG_COMM,IERROR)
213 cd write (iout,*) "Gather dphi"
215 write (iout,*) "dphi before gather"
217 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
220 call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),
221 & MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,
222 & king,FG_COMM,IERROR)
223 cd write (iout,*) "Gather dalpha"
226 call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),
227 & MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
228 & king,FG_COMM,IERROR)
229 cd write (iout,*) "Gather domega"
231 call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),
232 & MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,
233 & king,FG_COMM,IERROR)
238 write (iout,*) "dtheta after gather"
240 write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),j=1,2)
242 write (iout,*) "dphi after gather"
244 write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
250 subroutine checkintcartgrad
251 implicit real*8 (a-h,o-z)
256 include 'COMMON.CHAIN'
259 include 'COMMON.INTERACT'
260 include 'COMMON.DERIV'
261 include 'COMMON.IOUNITS'
262 include 'COMMON.SETUP'
263 double precision dthetanum(3,2,maxres),dphinum(3,3,maxres)
264 & ,dalphanum(3,3,maxres), domeganum(3,3,maxres)
265 double precision theta_s(maxres),phi_s(maxres),alph_s(maxres),
266 & omeg_s(maxres),dc_norm_s(3)
267 double precision aincr /1.0d-5/
275 c Check theta gradient
277 & "Analytical (upper) and numerical (lower) gradient of theta"
284 call int_from_cart1(.false.)
285 dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
288 dc(j,i-1)=dc(j,i-1)+aincr
290 dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
293 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),
294 & (dtheta(j,2,i),j=1,3)
295 write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),
296 & (dthetanum(j,2,i),j=1,3)
297 write (iout,'(5x,3f10.5,5x,3f10.5)')
298 & (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),
299 & (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
302 c Check gamma gradient
304 & "Analytical (upper) and numerical (lower) gradient of gamma"
310 dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
315 dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
318 dc(j,i-1)=dc(j,i-1)+aincr
320 dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
323 write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),
324 & (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
325 write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),
326 & (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
327 write (iout,'(5x,3(3f10.5,5x))')
328 & (dphinum(j,1,i)/dphi(j,1,i),j=1,3),
329 & (dphinum(j,2,i)/dphi(j,2,i),j=1,3),
330 & (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
333 c Check alpha gradient
335 & "Analytical (upper) and numerical (lower) gradient of alpha"
337 if(itype(i).ne.10) then
342 dalphanum(j,1,i)=(alph(i)-alph_s(i))
348 dalphanum(j,2,i)=(alph(i)-alph_s(i))
352 dc(j,i+nres)=dc(j,i+nres)+aincr
354 dalphanum(j,3,i)=(alph(i)-alph_s(i))
359 write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),
360 & (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
361 write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),
362 & (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
363 write (iout,'(5x,3(3f10.5,5x))')
364 & (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),
365 & (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),
366 & (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
369 c Check omega gradient
371 & "Analytical (upper) and numerical (lower) gradient of omega"
373 if(itype(i).ne.10) then
378 domeganum(j,1,i)=(omeg(i)-omeg_s(i))
384 domeganum(j,2,i)=(omeg(i)-omeg_s(i))
388 dc(j,i+nres)=dc(j,i+nres)+aincr
390 domeganum(j,3,i)=(omeg(i)-omeg_s(i))
395 write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),
396 & (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
397 write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),
398 & (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
399 write (iout,'(5x,3(3f10.5,5x))')
400 & (domeganum(j,1,i)/domega(j,1,i),j=1,3),
401 & (domeganum(j,2,i)/domega(j,2,i),j=1,3),
402 & (domeganum(j,3,i)/domega(j,3,i),j=1,3)
408 subroutine chainbuild_cart
409 implicit real*8 (a-h,o-z)
414 include 'COMMON.SETUP'
415 include 'COMMON.CHAIN'
416 include 'COMMON.LOCAL'
417 include 'COMMON.TIME1'
418 include 'COMMON.IOUNITS'
421 if (nfgtasks.gt.1) then
422 c write (iout,*) "BCAST in chainbuild_cart"
424 c Broadcast the order to build the chain and compute internal coordinates
425 c to the slaves. The slaves receive the order in ERGASTULUM.
427 c write (iout,*) "CHAINBUILD_CART: DC before BCAST"
429 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
430 c & (dc(j,i+nres),j=1,3)
433 & call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
434 time_bcast7=time_bcast7+MPI_Wtime()-time00
436 call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,
438 c write (iout,*) "CHAINBUILD_CART: DC after BCAST"
440 c write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
441 c & (dc(j,i+nres),j=1,3)
443 c write (iout,*) "End BCAST in chainbuild_cart"
445 time_bcast=time_bcast+MPI_Wtime()-time00
446 time_bcastc=time_bcastc+MPI_Wtime()-time01
454 c(j,i)=c(j,i-1)+dc(j,i-1)
459 c(j,i+nres)=c(j,i)+dc(j,i+nres)
462 c write (iout,*) "CHAINBUILD_CART"
464 call int_from_cart1(.false.)