1 C----------------------------------------------------------------------------
3 C Check the gradient of the energy in Cartesian coordinates.
6 include 'COMMON.CONTROL'
9 include 'COMMON.IOUNITS'
11 include 'COMMON.CONTACTS'
15 double precision ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
16 & g(maxvar),grad_s(6,maxres)
17 double precision energia(0:n_ene),energia1(0:n_ene)
18 double precision aincr2,etot,etot1,etot2
19 double precision dist,alpha,beta
20 double precision funcgrad,ff
24 double precision urparm(1)
31 print '("Calling CHECK_ECART",1pd12.3)',aincr
32 write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr
35 call geom_to_var(nvar,x)
36 call etotal(energia(0))
38 call enerprint(energia(0))
42 call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
46 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
50 grad_s(j,i)=gradc(j,i,icg)
51 grad_s(j+3,i)=gradx(j,i,icg)
55 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
66 c(j,k+nres)=c(j,k+nres)+aincr
68 call etotal(energia1(0))
70 ggg(j)=(etot1-etot)/aincr
74 c(j,k+nres)=c(j,k+nres)-aincr
78 c(j,i+nres)=c(j,i+nres)+aincr
79 dc(j,i+nres)=dc(j,i+nres)+aincr
80 call etotal(energia1(0))
82 ggg(j+3)=(etot1-etot)/aincr
86 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
87 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
91 c----------------------------------------------------------------------------
92 subroutine check_ecartint
93 C Check the gradient of the energy in Cartesian coordinates.
96 include 'COMMON.CONTROL'
97 include 'COMMON.CHAIN'
98 include 'COMMON.DERIV'
99 include 'COMMON.IOUNITS'
101 include 'COMMON.CONTACTS'
103 include 'COMMON.LOCAL'
104 include 'COMMON.SPLITELE'
107 double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
108 & x(maxvar),g(maxvar)
109 double precision dcnorm_safe(3),dxnorm_safe(3)
110 double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
111 double precision phi_temp(maxres),theta_temp(maxres),
112 & alph_temp(maxres),omeg_temp(maxres)
113 double precision energia(0:n_ene),energia1(0:n_ene)
115 double precision urparm(1)
118 double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
119 double precision dist,alpha,beta
126 call int_from_cart1(.false.)
129 c call checkintcartgrad
133 print '("Calling CHECK_ECARTINT",1pd12.3)',aincr
134 write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr
138 call geom_to_var(nvar,x)
139 if (.not.split_ene) then
140 call etotal(energia(0))
142 call enerprint(energia(0))
143 c write (iout,*) "enter cartgrad"
146 c write (iout,*) "exit cartgrad"
149 write (iout,'(//27(1h*)," Checking energy gradient ",27(1h*))')
150 write (iout,'(//4x,3a12,3x,3a12)')"gcart_x","gcart_y","gcart_z",
151 & "gxcart_x","gxcart_y","gxcart_z"
153 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
154 & (gxcart(j,i),j=1,3)
157 grad_s(j,0)=gcart(j,0)
161 grad_s(j,i)=gcart(j,i)
162 grad_s(j+3,i)=gxcart(j,i)
166 !- split gradient check
168 call etotal_long(energia(0))
169 call enerprint(energia(0))
171 c write (iout,*) "enter cartgrad"
174 c write (iout,*) "exit cartgrad"
177 write (iout,*) "longrange grad"
179 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
180 & (gxcart(j,i),j=1,3)
183 grad_s(j,0)=gcart(j,0)
187 grad_s(j,i)=gcart(j,i)
188 grad_s(j+3,i)=gxcart(j,i)
192 call etotal_short(energia(0))
193 call enerprint(energia(0))
195 c write (iout,*) "enter cartgrad"
198 c write (iout,*) "exit cartgrad"
201 write (iout,*) "shortrange grad"
203 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
204 & (gxcart(j,i),j=1,3)
207 grad_s1(j,0)=gcart(j,0)
211 grad_s1(j,i)=gcart(j,i)
212 grad_s1(j+3,i)=gxcart(j,i)
216 c write (iout,*) "Vector dc"
218 c write (iout,'(i5,2(3f10.5,5x))')
219 c & i,(dc(j,i),j=1,3),(dc(j,i+nres),j=1,3)
221 c write (iout,*) "Coordinates after chainbuild_cart"
222 c call chainbuild_cart
224 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
232 dcnorm_safe(k)=dc_norm(k,i)
233 dxnorm_safe(k)=dc_norm(k,i+nres)
240 c Broadcast the order to compute internal coordinates to the slaves.
242 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
244 c call int_from_cart1(.false.)
245 if (.not.split_ene) then
246 call etotal(energia1(0))
250 call etotal_long(energia1(0))
252 call etotal_short(energia1(0))
254 c write (iout,*) "etot11",etot11," etot12",etot12
256 !- end split gradient
257 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
261 c call int_from_cart1(.false.)
262 if (.not.split_ene) then
263 call etotal(energia1(0))
265 ggg(j)=(etot1-etot2)/(2*aincr)
268 call etotal_long(energia1(0))
270 ggg(j)=(etot11-etot21)/(2*aincr)
271 call etotal_short(energia1(0))
273 ggg1(j)=(etot12-etot22)/(2*aincr)
274 !- end split gradient
275 c write (iout,*) "etot21",etot21," etot22",etot22
277 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
282 dc(j,i+nres)=ddx(j)+aincr
284 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
285 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
286 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
287 c write (iout,*) "dxnormnorm",dsqrt(
288 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
289 c write (iout,*) "dxnormnormsafe",dsqrt(
290 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
292 if (.not.split_ene) then
293 call etotal(energia1(0))
297 call etotal_long(energia1(0))
299 call etotal_short(energia1(0))
302 !- end split gradient
303 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
304 dc(j,i+nres)=ddx(j)-aincr
306 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
307 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
308 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
310 c write (iout,*) "dxnormnorm",dsqrt(
311 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
312 c write (iout,*) "dxnormnormsafe",dsqrt(
313 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
314 if (.not.split_ene) then
315 call etotal(energia1(0))
317 ggg(j+3)=(etot1-etot2)/(2*aincr)
320 call etotal_long(energia1(0))
322 ggg(j+3)=(etot11-etot21)/(2*aincr)
323 call etotal_short(energia1(0))
325 ggg1(j+3)=(etot12-etot22)/(2*aincr)
326 !- end split gradient
328 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
332 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
333 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
335 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
336 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
338 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
339 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
340 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
345 c-------------------------------------------------------------------------
346 subroutine int_from_cart1(lprn)
353 include 'COMMON.IOUNITS'
355 include 'COMMON.CHAIN'
357 include 'COMMON.INTERACT'
358 include 'COMMON.LOCAL'
359 include 'COMMON.NAMES'
360 include 'COMMON.SETUP'
361 include 'COMMON.TIME1'
364 double precision dnorm1,dnorm2,be
365 double precision time00
366 double precision dist,alpha,beta
367 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
371 #if defined(PARINT) && defined(MPI)
372 do i=iint_start,iint_end
379 C print *,i,dnorm1,dnorm2
381 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
382 & +(c(j,i+1)-c(j,i))/dnorm2)
386 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
387 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
388 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
390 if (itype(i-1).ne.10) then
391 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
392 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
393 omicron(2,i)=alpha(i-1+nres,i-1,i)
395 if (itype(i).ne.10) then
396 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
399 omeg(i)=beta(nres+i,i,maxres2,i+1)
401 alph(i)=alpha(nres+i,i,maxres2)
403 theta(i+1)=alpha(i-1,i,i+1)
406 vbld_inv(i)=1.0d0/vbld(i)
407 vbld(nres+i)=dist(nres+i,i)
408 C print *,vbld(i+nres)
410 if (itype(i).ne.10) then
411 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
413 vbld_inv(nres+i)=0.0d0
416 #if defined(PARINT) && defined(MPI)
417 if (nfgtasks1.gt.1) then
418 cd write(iout,*) "iint_start",iint_start," iint_count",
419 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
420 cd & (iint_displ(i),i=0,nfgtasks-1)
421 cd write (iout,*) "Gather vbld backbone"
424 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
425 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
426 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
427 cd write (iout,*) "Gather vbld_inv"
429 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
430 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
431 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
432 cd write (iout,*) "Gather vbld side chain"
434 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
435 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
436 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
437 cd write (iout,*) "Gather vbld_inv side chain"
439 call MPI_Allgatherv(vbld_inv(iint_start+nres),
440 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
441 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
442 cd write (iout,*) "Gather theta"
444 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
445 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
446 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
447 cd write (iout,*) "Gather phi"
449 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
450 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
451 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
453 cd write (iout,*) "Gather alph"
455 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
456 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
457 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
458 cd write (iout,*) "Gather omeg"
460 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
461 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
462 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
464 time_gather=time_gather+MPI_Wtime()-time00
469 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
474 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
479 write (iout,1212) restyp(itype(i)),i,vbld(i),
480 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
481 &rad2deg*alph(i),rad2deg*omeg(i)
484 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
486 time_intfcart=time_intfcart+MPI_Wtime()-time01
490 c----------------------------------------------------------------------------
491 subroutine check_eint
492 C Check the gradient of energy in internal coordinates.
495 include 'COMMON.CONTROL'
496 include 'COMMON.CHAIN'
497 include 'COMMON.DERIV'
498 include 'COMMON.IOUNITS'
503 double precision x(maxvar),gana(maxvar),gg(maxvar)
505 double precision urparm(1)
506 double precision energia(0:n_ene),energia1(0:n_ene),
509 double precision fdum
511 double precision funcgrad,ff
514 double precision xi,etot,etot1,etot2
517 print '("Calling CHECK_INT",1pd12.3)',aincr
518 write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
522 call geom_to_var(nvar,x)
523 call var_to_geom(nvar,x)
527 call etotal(energia(0))
529 call enerprint(energia(0))
532 if (MyID.ne.BossID) then
533 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
541 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
542 c write (iout,*) "Before gradient"
547 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
549 c write (iout,*) "After gradient"
551 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
556 call var_to_geom(nvar,x)
557 call chainbuild_extconf
558 call etotal(energia1(0))
561 call var_to_geom(nvar,x)
562 call chainbuild_extconf
563 call etotal(energia2(0))
565 gg(i)=(etot2-etot1)/aincr
566 c write (iout,*) i,etot1,etot2
569 write (iout,'(/2a)')' Variable Numerical Analytical',
575 else if (i.le.nphi+ntheta) then
578 else if (i.le.nphi+ntheta+nside) then
582 ii=i-(nphi+ntheta+nside)
585 write (iout,'(i3,a,i3,3(1pd16.6))')
586 & i,key,ii,gg(i),gana(i),
587 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)