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)
92 !-----------------------------------------------------------------------------
93 subroutine check_ecartint
94 ! Check the gradient of the energy in Cartesian coordinates.
97 include 'COMMON.CONTROL'
98 include 'COMMON.CHAIN'
99 include 'COMMON.INTERACT'
100 include 'COMMON.DERIV'
101 include 'COMMON.IOUNITS'
103 include 'COMMON.CONTACTS'
105 include 'COMMON.LOCAL'
106 include 'COMMON.SPLITELE'
109 double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
110 & x(maxvar),g(maxvar)
111 double precision dcnorm_safe(3),dxnorm_safe(3)
112 double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
113 double precision phi_temp(maxres),theta_temp(maxres),
114 & alph_temp(maxres),omeg_temp(maxres)
115 double precision ddc1(3),ddcn(3),dcnorm_safe1(3),dcnorm_safe2(3)
116 double precision energia(0:n_ene),energia1(0:n_ene)
118 double precision urparm(1)
119 double precision fdum
122 double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
123 double precision dist,alpha,beta
129 ! call checkintcartgrad
132 write(iout,*) 'Calling CHECK_ECARTINT.'
135 write (iout,*) "Before geom_to_var"
136 call geom_to_var(nvar,x)
137 write (iout,*) "after geom_to_var"
138 write (iout,*) "split_ene ",split_ene
140 if (.not.split_ene) then
141 write(iout,*) 'Calling CHECK_ECARTINT if'
143 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
145 write (iout,*) "etot",etot
146 call enerprint(energia(0))
148 !el call enerprint(energia)
149 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
151 c write (iout,*) "enter cartgrad"
154 c Transform the gradient to the CA-SC basis
156 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
157 c write (iout,*) "exit cartgrad"
161 c write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
164 grad_s(j,0)=gcart(j,0)
166 !elwrite(iout,*) 'Calling CHECK_ECARTINT if'
169 grad_s(j,i)=gcart(j,i)
170 grad_s(j+3,i)=gxcart(j,i)
174 c write(iout,*) 'Calling CHECK_ECARTIN else.'
175 !- split gradient check
177 call etotal_long(energia)
178 call enerprint(energia(0))
179 !el call enerprint(energia)
181 c write (iout,*) "enter cartgrad"
184 c Transform the gradient to CA-SC coordinates
186 c write (iout,*) "exit cartgrad"
189 write (iout,*) "longrange grad"
191 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
192 & (gxcart(j,i),j=1,3)
195 grad_s(j,0)=gcart(j,0)
199 grad_s(j,i)=gcart(j,i)
200 grad_s(j+3,i)=gxcart(j,i)
204 call etotal_short(energia)
205 call enerprint(energia(0))
207 c write (iout,*) "enter cartgrad"
210 c write (iout,*) "exit cartgrad"
212 c Transform the gradient to CA-SC basis
215 write (iout,*) "shortrange grad"
217 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
218 & (gxcart(j,i),j=1,3)
221 grad_s1(j,0)=gcart(j,0)
225 grad_s1(j,i)=gcart(j,i)
226 grad_s1(j+3,i)=gxcart(j,i)
230 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
235 if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
236 if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
239 dcnorm_safe1(j)=dc_norm(j,i-1)
240 dcnorm_safe2(j)=dc_norm(j,i)
241 dxnorm_safe(j)=dc_norm(j,i+nres)
245 if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
246 if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
247 if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
248 dc(j,i)=c(j,i+1)-c(j,i)
249 dc(j,i+nres)=c(j,i+nres)-c(j,i)
250 call int_from_cart1(.false.)
251 if (.not.split_ene) then
252 call etotal(energia1)
254 c write (iout,*) "ij",i,j," etot1",etot1
257 call etotal_long(energia1)
259 call etotal_short(energia1)
262 !- end split gradient
263 ! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
265 if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
266 if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
267 if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
268 dc(j,i)=c(j,i+1)-c(j,i)
269 dc(j,i+nres)=c(j,i+nres)-c(j,i)
270 call int_from_cart1(.false.)
271 if (.not.split_ene) then
272 call etotal(energia1)
274 c write (iout,*) "ij",i,j," etot2",etot2
275 ggg(j)=(etot1-etot2)/(2*aincr)
278 call etotal_long(energia1)
280 ggg(j)=(etot11-etot21)/(2*aincr)
281 call etotal_short(energia1)
283 ggg1(j)=(etot12-etot22)/(2*aincr)
284 !- end split gradient
285 ! write (iout,*) "etot21",etot21," etot22",etot22
287 ! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
289 if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
290 if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
291 if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
292 dc(j,i)=c(j,i+1)-c(j,i)
293 dc(j,i+nres)=c(j,i+nres)-c(j,i)
294 dc_norm(j,i-1)=dcnorm_safe1(j)
295 dc_norm(j,i)=dcnorm_safe2(j)
296 dc_norm(j,i+nres)=dxnorm_safe(j)
299 c(j,i+nres)=ddx(j)+aincr
300 dc(j,i+nres)=c(j,i+nres)-c(j,i)
301 call int_from_cart1(.false.)
302 if (.not.split_ene) then
303 call etotal(energia1)
307 call etotal_long(energia1)
309 call etotal_short(energia1)
312 !- end split gradient
313 c(j,i+nres)=ddx(j)-aincr
314 dc(j,i+nres)=c(j,i+nres)-c(j,i)
315 call int_from_cart1(.false.)
316 if (.not.split_ene) then
317 call etotal(energia1)
319 ggg(j+3)=(etot1-etot2)/(2*aincr)
322 call etotal_long(energia1)
324 ggg(j+3)=(etot11-etot21)/(2*aincr)
325 call etotal_short(energia1)
327 ggg1(j+3)=(etot12-etot22)/(2*aincr)
328 !- end split gradient
330 ! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
332 dc(j,i+nres)=c(j,i+nres)-c(j,i)
333 dc_norm(j,i+nres)=dxnorm_safe(j)
334 call int_from_cart1(.false.)
336 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
337 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
339 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
340 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
342 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
343 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
344 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
350 c----------------------------------------------------------------------------
351 subroutine check_ecartint
352 C Check the gradient of the energy in Cartesian coordinates.
355 include 'COMMON.CONTROL'
356 include 'COMMON.CHAIN'
357 include 'COMMON.DERIV'
358 include 'COMMON.IOUNITS'
360 include 'COMMON.CONTACTS'
362 include 'COMMON.LOCAL'
363 include 'COMMON.SPLITELE'
366 double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
367 & x(maxvar),g(maxvar)
368 double precision dcnorm_safe(3),dxnorm_safe(3)
369 double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
370 double precision phi_temp(maxres),theta_temp(maxres),
371 & alph_temp(maxres),omeg_temp(maxres)
372 double precision energia(0:n_ene),energia1(0:n_ene)
374 double precision urparm(1)
377 double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
378 double precision dist,alpha,beta
385 call int_from_cart1(.false.)
388 c call checkintcartgrad
392 print '("Calling CHECK_ECARTINT",1pd12.3)',aincr
393 write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr
396 call geom_to_var(nvar,x)
397 if (.not.split_ene) then
398 call etotal(energia(0))
400 call enerprint(energia(0))
401 c write (iout,*) "enter cartgrad"
404 c write (iout,*) "exit cartgrad"
407 write (iout,'(//27(1h*)," Checking energy gradient ",27(1h*))')
408 write (iout,'(//4x,3a12,3x,3a12)')"gcart_x","gcart_y","gcart_z",
409 & "gxcart_x","gxcart_y","gxcart_z"
411 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
412 & (gxcart(j,i),j=1,3)
415 grad_s(j,0)=gcart(j,0)
419 grad_s(j,i)=gcart(j,i)
420 grad_s(j+3,i)=gxcart(j,i)
424 !- split gradient check
426 call etotal_long(energia(0))
427 call enerprint(energia(0))
429 write (iout,*) "enter cartgrad"
432 write (iout,*) "exit cartgrad"
435 write (iout,*) "longrange grad"
437 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
438 & (gxcart(j,i),j=1,3)
441 grad_s(j,0)=gcart(j,0)
445 grad_s(j,i)=gcart(j,i)
446 grad_s(j+3,i)=gxcart(j,i)
450 call etotal_short(energia(0))
451 call enerprint(energia(0))
453 write (iout,*) "enter cartgrad"
456 write (iout,*) "exit cartgrad"
459 write (iout,*) "shortrange grad"
461 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
462 & (gxcart(j,i),j=1,3)
465 grad_s1(j,0)=gcart(j,0)
469 grad_s1(j,i)=gcart(j,i)
470 grad_s1(j+3,i)=gxcart(j,i)
474 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
482 dcnorm_safe(k)=dc_norm(k,i)
483 dxnorm_safe(k)=dc_norm(k,i+nres)
490 c Broadcast the order to compute internal coordinates to the slaves.
492 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
494 c call int_from_cart1(.false.)
495 if (.not.split_ene) then
496 call etotal(energia1(0))
500 call etotal_long(energia1(0))
502 call etotal_short(energia1(0))
504 c write (iout,*) "etot11",etot11," etot12",etot12
506 !- end split gradient
507 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
511 c call int_from_cart1(.false.)
512 if (.not.split_ene) then
513 call etotal(energia1(0))
515 ggg(j)=(etot1-etot2)/(2*aincr)
518 call etotal_long(energia1(0))
520 ggg(j)=(etot11-etot21)/(2*aincr)
521 call etotal_short(energia1(0))
523 ggg1(j)=(etot12-etot22)/(2*aincr)
524 !- end split gradient
525 c write (iout,*) "etot21",etot21," etot22",etot22
527 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
532 dc(j,i+nres)=ddx(j)+aincr
534 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
535 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
536 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
537 c write (iout,*) "dxnormnorm",dsqrt(
538 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
539 c write (iout,*) "dxnormnormsafe",dsqrt(
540 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
542 if (.not.split_ene) then
543 call etotal(energia1(0))
547 call etotal_long(energia1(0))
549 call etotal_short(energia1(0))
552 !- end split gradient
553 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
554 dc(j,i+nres)=ddx(j)-aincr
556 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
557 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
558 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
560 c write (iout,*) "dxnormnorm",dsqrt(
561 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
562 c write (iout,*) "dxnormnormsafe",dsqrt(
563 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
564 if (.not.split_ene) then
565 call etotal(energia1(0))
567 ggg(j+3)=(etot1-etot2)/(2*aincr)
570 call etotal_long(energia1(0))
572 ggg(j+3)=(etot11-etot21)/(2*aincr)
573 call etotal_short(energia1(0))
575 ggg1(j+3)=(etot12-etot22)/(2*aincr)
576 !- end split gradient
578 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
582 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
583 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
585 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
586 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
588 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
589 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
590 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
596 c-------------------------------------------------------------------------
597 subroutine int_from_cart1(lprn)
604 include 'COMMON.IOUNITS'
606 include 'COMMON.CHAIN'
608 include 'COMMON.INTERACT'
609 include 'COMMON.LOCAL'
610 include 'COMMON.NAMES'
611 include 'COMMON.SETUP'
612 include 'COMMON.TIME1'
615 double precision dnorm1,dnorm2,be
616 double precision time00
617 double precision dist,alpha,beta
618 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
622 #if defined(PARINT) && defined(MPI)
623 do i=iint_start,iint_end
630 C print *,i,dnorm1,dnorm2
632 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
633 & +(c(j,i+1)-c(j,i))/dnorm2)
637 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
638 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
639 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
641 if (itype(i-1).ne.10) then
642 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
643 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
644 omicron(2,i)=alpha(i-1+nres,i-1,i)
646 if (itype(i).ne.10) then
647 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
650 omeg(i)=beta(nres+i,i,maxres2,i+1)
652 alph(i)=alpha(nres+i,i,maxres2)
654 theta(i+1)=alpha(i-1,i,i+1)
657 vbld_inv(i)=1.0d0/vbld(i)
658 vbld(nres+i)=dist(nres+i,i)
659 C print *,vbld(i+nres)
661 if (itype(i).ne.10) then
662 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
664 vbld_inv(nres+i)=0.0d0
667 #if defined(PARINT) && defined(MPI)
668 if (nfgtasks1.gt.1) then
669 cd write(iout,*) "iint_start",iint_start," iint_count",
670 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
671 cd & (iint_displ(i),i=0,nfgtasks-1)
672 cd write (iout,*) "Gather vbld backbone"
675 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
676 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
677 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
678 cd write (iout,*) "Gather vbld_inv"
680 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
681 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
682 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
683 cd write (iout,*) "Gather vbld side chain"
685 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
686 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
687 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
688 cd write (iout,*) "Gather vbld_inv side chain"
690 call MPI_Allgatherv(vbld_inv(iint_start+nres),
691 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
692 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
693 cd write (iout,*) "Gather theta"
695 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
696 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
697 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
698 cd write (iout,*) "Gather phi"
700 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
701 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
702 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
704 cd write (iout,*) "Gather alph"
706 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
707 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
708 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
709 cd write (iout,*) "Gather omeg"
711 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
712 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
713 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
715 time_gather=time_gather+MPI_Wtime()-time00
720 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
725 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
730 write (iout,1212) restyp(itype(i)),i,vbld(i),
731 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
732 &rad2deg*alph(i),rad2deg*omeg(i)
735 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
737 time_intfcart=time_intfcart+MPI_Wtime()-time01
741 c----------------------------------------------------------------------------
742 subroutine check_eint
743 C Check the gradient of energy in internal coordinates.
746 include 'COMMON.CONTROL'
747 include 'COMMON.CHAIN'
748 include 'COMMON.DERIV'
749 include 'COMMON.IOUNITS'
754 double precision x(maxvar),gana(maxvar),gg(maxvar)
756 double precision urparm(1)
757 double precision energia(0:n_ene),energia1(0:n_ene),
760 double precision fdum
762 double precision funcgrad,ff
765 double precision xi,etot,etot1,etot2
768 print '("Calling CHECK_INT",1pd12.3)',aincr
769 write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
773 call geom_to_var(nvar,x)
774 call var_to_geom(nvar,x)
778 call etotal(energia(0))
780 call enerprint(energia(0))
783 if (MyID.ne.BossID) then
784 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
792 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
793 c write (iout,*) "Before gradient"
798 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
800 c write (iout,*) "After gradient"
802 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
807 call var_to_geom(nvar,x)
808 call chainbuild_extconf
809 call etotal(energia1(0))
812 call var_to_geom(nvar,x)
813 call chainbuild_extconf
814 call etotal(energia2(0))
816 gg(i)=(etot2-etot1)/aincr
817 c write (iout,*) i,etot1,etot2
820 write (iout,'(/2a)')' Variable Numerical Analytical',
826 else if (i.le.nphi+ntheta) then
829 else if (i.le.nphi+ntheta+nside) then
833 ii=i-(nphi+ntheta+nside)
836 write (iout,'(i3,a,i3,3(1pd16.6))')
837 & i,key,ii,gg(i),gana(i),
838 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)