source/unres/src-HCD-5D/checkder_p.F

   1 C----------------------------------------------------------------------------
   2       subroutine check_ecart
   3 C Check the gradient of the energy in Cartesian coordinates.
   4       implicit none
   5       include 'DIMENSIONS'
   6       include 'COMMON.CONTROL'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.DERIV'
   9       include 'COMMON.IOUNITS'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CONTACTS'
  12       integer i,j,k
  13       integer icall
  14       common /srutu/ icall
  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
  21       external funcgrad
  22       integer nf
  23       integer uiparm(1)
  24       double precision urparm(1)
  25       double precision fdum
  26       external fdum
  27       icg=1
  28       nf=0
  29       nfl=0
  30       call zerograd
  31       print '("Calling CHECK_ECART",1pd12.3)',aincr
  32       write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr
  33       nf=0
  34       icall=0
  35       call geom_to_var(nvar,x)
  36       call etotal(energia(0))
  37       etot=energia(0)
  38       call enerprint(energia(0))
  39 #ifdef LBFGS
  40       ff=funcgrad(x,g)
  41 #else
  42       call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
  43 #endif
  44       icall =1
  45       do i=1,nres
  46         write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
  47       enddo
  48       do i=1,nres
  49         do j=1,3
  50           grad_s(j,i)=gradc(j,i,icg)
  51           grad_s(j+3,i)=gradx(j,i,icg)
  52         enddo
  53       enddo
  54       call flush(iout)
  55       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
  56       do i=1,nres
  57         do j=1,3
  58           xx(j)=c(j,i+nres)
  59           ddc(j)=dc(j,i)
  60           ddx(j)=dc(j,i+nres)
  61         enddo
  62         do j=1,3
  63           dc(j,i)=dc(j,i)+aincr
  64           do k=i+1,nres
  65             c(j,k)=c(j,k)+aincr
  66             c(j,k+nres)=c(j,k+nres)+aincr
  67           enddo
  68           call etotal(energia1(0))
  69           etot1=energia1(0)
  70           ggg(j)=(etot1-etot)/aincr
  71           dc(j,i)=ddc(j)
  72           do k=i+1,nres
  73             c(j,k)=c(j,k)-aincr
  74             c(j,k+nres)=c(j,k+nres)-aincr
  75           enddo
  76         enddo
  77         do j=1,3
  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))
  81           etot1=energia1(0)
  82           ggg(j+3)=(etot1-etot)/aincr
  83           c(j,i+nres)=xx(j)
  84           dc(j,i+nres)=ddx(j)
  85         enddo
  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)
  88       enddo
  89       return
  90       end
  91 c----------------------------------------------------------------------------
  92       subroutine check_ecartint
  93 C Check the gradient of the energy in Cartesian coordinates.
  94       implicit none
  95       include 'DIMENSIONS'
  96       include 'COMMON.CONTROL'
  97       include 'COMMON.CHAIN'
  98       include 'COMMON.DERIV'
  99       include 'COMMON.IOUNITS'
 100       include 'COMMON.VAR'
 101       include 'COMMON.CONTACTS'
 102       include 'COMMON.MD'
 103       include 'COMMON.LOCAL'
 104       include 'COMMON.SPLITELE'
 105       integer icall
 106       common /srutu/ icall
 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)
 114       integer uiparm(1)
 115       double precision urparm(1)
 116       external fdum
 117       integer i,j,k,nf
 118       double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
 119       double precision dist,alpha,beta
 120 c      r_cut=2.0d0
 121 c      rlambd=0.3d0
 122       icg=1
 123       nf=0
 124       nfl=0
 125 c      print *,"ATU 3"
 126       call int_from_cart1(.false.)
 127       call intout
 128 c      call intcartderiv
 129 c      call checkintcartgrad
 130       call zerograd
 131 c      aincr=8.0D-7
 132 c      aincr=1.0D-7
 133       print '("Calling CHECK_ECARTINT",1pd12.3)',aincr
 134       write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr
 135       nf=0
 136       icall=0
 137       call geom_to_var(nvar,x)
 138       if (.not.split_ene) then
 139         call etotal(energia(0))
 140         etot=energia(0)
 141         call enerprint(energia(0))
 142 c        write (iout,*) "enter cartgrad"
 143 c        call flush(iout)
 144         call cartgrad
 145 c        write (iout,*) "exit cartgrad"
 146 c        call flush(iout)
 147         icall =1
 148         write (iout,'(//27(1h*)," Checking energy gradient ",27(1h*))')
 149         write (iout,'(//4x,3a12,3x,3a12)')"gcart_x","gcart_y","gcart_z",
 150      &     "gxcart_x","gxcart_y","gxcart_z"
 151         do i=1,nres
 152           write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
 153      &      (gxcart(j,i),j=1,3)
 154         enddo
 155         do j=1,3
 156           grad_s(j,0)=gcart(j,0)
 157         enddo
 158         do i=1,nres
 159           do j=1,3
 160             grad_s(j,i)=gcart(j,i)
 161             grad_s(j+3,i)=gxcart(j,i)
 162           enddo
 163         enddo
 164       else
 165 !- split gradient check
 166         call zerograd
 167         call etotal_long(energia(0))
 168         call enerprint(energia(0))
 169         call flush(iout)
 170 c        write (iout,*) "enter cartgrad"
 171 c        call flush(iout)
 172         call cartgrad
 173 c        write (iout,*) "exit cartgrad"
 174 c        call flush(iout)
 175         icall =1
 176         write (iout,*) "longrange grad"
 177         do i=1,nres
 178           write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
 179      &    (gxcart(j,i),j=1,3)
 180         enddo
 181         do j=1,3
 182           grad_s(j,0)=gcart(j,0)
 183         enddo
 184         do i=1,nres
 185           do j=1,3
 186             grad_s(j,i)=gcart(j,i)
 187             grad_s(j+3,i)=gxcart(j,i)
 188           enddo
 189         enddo
 190         call zerograd
 191         call etotal_short(energia(0))
 192         call enerprint(energia(0))
 193         call flush(iout)
 194 c        write (iout,*) "enter cartgrad"
 195 c        call flush(iout)
 196         call cartgrad
 197 c        write (iout,*) "exit cartgrad"
 198 c        call flush(iout)
 199         icall =1
 200         write (iout,*) "shortrange grad"
 201         do i=1,nres
 202           write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
 203      &    (gxcart(j,i),j=1,3)
 204         enddo
 205         do j=1,3
 206           grad_s1(j,0)=gcart(j,0)
 207         enddo
 208         do i=1,nres
 209           do j=1,3
 210             grad_s1(j,i)=gcart(j,i)
 211             grad_s1(j+3,i)=gxcart(j,i)
 212           enddo
 213         enddo
 214       endif
 215       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
 216       do i=0,nres
 217         print *,i
 218         do j=1,3
 219           xx(j)=c(j,i+nres)
 220           ddc(j)=dc(j,i)
 221           ddx(j)=dc(j,i+nres)
 222           do k=1,3
 223             dcnorm_safe(k)=dc_norm(k,i)
 224             dxnorm_safe(k)=dc_norm(k,i+nres)
 225           enddo
 226         enddo
 227         do j=1,3
 228           dc(j,i)=ddc(j)+aincr
 229           call chainbuild_cart
 230 #ifdef MPI
 231 c Broadcast the order to compute internal coordinates to the slaves.
 232 c          if (nfgtasks.gt.1)
 233 c     &      call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
 234 #endif
 235 c          call int_from_cart1(.false.)
 236           if (.not.split_ene) then
 237             call etotal(energia1(0))
 238             etot1=energia1(0)
 239           else
 240 !- split gradient
 241             call etotal_long(energia1(0))
 242             etot11=energia1(0)
 243             call etotal_short(energia1(0))
 244             etot12=energia1(0)
 245 c            write (iout,*) "etot11",etot11," etot12",etot12
 246           endif
 247 !- end split gradient
 248 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
 249           dc(j,i)=ddc(j)-aincr
 250           call chainbuild_cart
 251 C          print *,c(j,i)
 252 c          call int_from_cart1(.false.)
 253           if (.not.split_ene) then
 254             call etotal(energia1(0))
 255             etot2=energia1(0)
 256             ggg(j)=(etot1-etot2)/(2*aincr)
 257           else
 258 !- split gradient
 259             call etotal_long(energia1(0))
 260             etot21=energia1(0)
 261             ggg(j)=(etot11-etot21)/(2*aincr)
 262             call etotal_short(energia1(0))
 263             etot22=energia1(0)
 264             ggg1(j)=(etot12-etot22)/(2*aincr)
 265 !- end split gradient
 266 c            write (iout,*) "etot21",etot21," etot22",etot22
 267           endif
 268 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 269           dc(j,i)=ddc(j)
 270           call chainbuild_cart
 271         enddo
 272         do j=1,3
 273           dc(j,i+nres)=ddx(j)+aincr
 274           call chainbuild_cart
 275 c          write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
 276 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
 277 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
 278 c          write (iout,*) "dxnormnorm",dsqrt(
 279 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
 280 c          write (iout,*) "dxnormnormsafe",dsqrt(
 281 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
 282 c          write (iout,*)
 283           if (.not.split_ene) then
 284             call etotal(energia1(0))
 285             etot1=energia1(0)
 286           else
 287 !- split gradient
 288             call etotal_long(energia1(0))
 289             etot11=energia1(0)
 290             call etotal_short(energia1(0))
 291             etot12=energia1(0)
 292           endif
 293 !- end split gradient
 294 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
 295           dc(j,i+nres)=ddx(j)-aincr
 296           call chainbuild_cart
 297 c          write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
 298 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
 299 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
 300 c          write (iout,*)
 301 c          write (iout,*) "dxnormnorm",dsqrt(
 302 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
 303 c          write (iout,*) "dxnormnormsafe",dsqrt(
 304 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
 305           if (.not.split_ene) then
 306             call etotal(energia1(0))
 307             etot2=energia1(0)
 308             ggg(j+3)=(etot1-etot2)/(2*aincr)
 309           else
 310 !- split gradient
 311             call etotal_long(energia1(0))
 312             etot21=energia1(0)
 313             ggg(j+3)=(etot11-etot21)/(2*aincr)
 314             call etotal_short(energia1(0))
 315             etot22=energia1(0)
 316             ggg1(j+3)=(etot12-etot22)/(2*aincr)
 317 !- end split gradient
 318           endif
 319 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 320           dc(j,i+nres)=ddx(j)
 321           call chainbuild_cart
 322         enddo
 323         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 324      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
 325         if (split_ene) then
 326           write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 327      &   i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
 328      &   k=1,6)
 329          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 330      &   i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
 331      &   ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
 332         endif
 333       enddo
 334       return
 335       end
 336 c-------------------------------------------------------------------------
 337       subroutine int_from_cart1(lprn)
 338       implicit none
 339       include 'DIMENSIONS'
 340 #ifdef MPI
 341       include 'mpif.h'
 342       integer ierror
 343 #endif
 344       include 'COMMON.IOUNITS'
 345       include 'COMMON.VAR'
 346       include 'COMMON.CHAIN'
 347       include 'COMMON.GEO'
 348       include 'COMMON.INTERACT'
 349       include 'COMMON.LOCAL'
 350       include 'COMMON.NAMES'
 351       include 'COMMON.SETUP'
 352       include 'COMMON.TIME1'
 353       logical lprn
 354       integer i,j
 355       double precision dnorm1,dnorm2,be
 356       double precision time00
 357       double precision dist,alpha,beta
 358       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 359 #ifdef TIMING
 360       time01=MPI_Wtime()
 361 #endif
 362 #if defined(PARINT) && defined(MPI)
 363       do i=iint_start,iint_end
 364 #else
 365       do i=2,nres
 366 #endif
 367 C        print *,i
 368         dnorm1=dist(i-1,i)
 369         dnorm2=dist(i,i+1)
 370 C        print *,i,dnorm1,dnorm2
 371         do j=1,3
 372           c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
 373      &     +(c(j,i+1)-c(j,i))/dnorm2)
 374         enddo
 375         be=0.0D0
 376         if (i.gt.2) then
 377         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
 378         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
 379          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
 380         endif
 381         if (itype(i-1).ne.10) then
 382          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
 383          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
 384          omicron(2,i)=alpha(i-1+nres,i-1,i)
 385         endif
 386         if (itype(i).ne.10) then
 387          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
 388         endif
 389         endif
 390         omeg(i)=beta(nres+i,i,maxres2,i+1)
 391 C        print *,omeg(i)
 392         alph(i)=alpha(nres+i,i,maxres2)
 393 C        print *,alph(i)
 394         theta(i+1)=alpha(i-1,i,i+1)
 395         vbld(i)=dist(i-1,i)
 396 C        print *,vbld(i)
 397         vbld_inv(i)=1.0d0/vbld(i)
 398         vbld(nres+i)=dist(nres+i,i)
 399 C            print *,vbld(i+nres)
 400
 401         if (itype(i).ne.10) then
 402           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
 403         else
 404           vbld_inv(nres+i)=0.0d0
 405         endif
 406       enddo
 407 #if defined(PARINT) && defined(MPI)
 408        if (nfgtasks1.gt.1) then
 409 cd       write(iout,*) "iint_start",iint_start," iint_count",
 410 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
 411 cd     &   (iint_displ(i),i=0,nfgtasks-1)
 412 cd       write (iout,*) "Gather vbld backbone"
 413 cd       call flush(iout)
 414        time00=MPI_Wtime()
 415        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
 416      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
 417      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 418 cd       write (iout,*) "Gather vbld_inv"
 419 cd       call flush(iout)
 420        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
 421      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
 422      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 423 cd       write (iout,*) "Gather vbld side chain"
 424 cd       call flush(iout)
 425        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
 426      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
 427      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 428 cd       write (iout,*) "Gather vbld_inv side chain"
 429 cd       call flush(iout)
 430        call MPI_Allgatherv(vbld_inv(iint_start+nres),
 431      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
 432      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 433 cd       write (iout,*) "Gather theta"
 434 cd       call flush(iout)
 435        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
 436      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
 437      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 438 cd       write (iout,*) "Gather phi"
 439 cd       call flush(iout)
 440        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
 441      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
 442      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 443 #ifdef CRYST_SC
 444 cd       write (iout,*) "Gather alph"
 445 cd       call flush(iout)
 446        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
 447      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
 448      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 449 cd       write (iout,*) "Gather omeg"
 450 cd       call flush(iout)
 451        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
 452      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
 453      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 454 #endif
 455        time_gather=time_gather+MPI_Wtime()-time00
 456       endif
 457 #endif
 458       do i=1,nres-1
 459         do j=1,3
 460           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 461         enddo
 462       enddo
 463       do i=2,nres-1
 464         do j=1,3
 465           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 466         enddo
 467       enddo
 468       if (lprn) then
 469       do i=2,nres
 470        write (iout,1212) restyp(itype(i)),i,vbld(i),
 471      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
 472      &rad2deg*alph(i),rad2deg*omeg(i)
 473       enddo
 474       endif
 475  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
 476 #ifdef TIMING
 477       time_intfcart=time_intfcart+MPI_Wtime()-time01
 478 #endif
 479       return
 480       end
 481 c----------------------------------------------------------------------------
 482       subroutine check_eint
 483 C Check the gradient of energy in internal coordinates.
 484       implicit none
 485       include 'DIMENSIONS'
 486       include 'COMMON.CONTROL'
 487       include 'COMMON.CHAIN'
 488       include 'COMMON.DERIV'
 489       include 'COMMON.IOUNITS'
 490       include 'COMMON.VAR'
 491       include 'COMMON.GEO'
 492       integer icall
 493       common /srutu/ icall
 494       double precision x(maxvar),gana(maxvar),gg(maxvar)
 495       integer uiparm(1)
 496       double precision urparm(1)
 497       double precision energia(0:n_ene),energia1(0:n_ene),
 498      &  energia2(0:n_ene)
 499       character*6 key
 500       double precision fdum
 501       external fdum
 502       double precision funcgrad,ff
 503       external funcgrad
 504       integer i,ii,nf
 505       double precision xi,etot,etot1,etot2
 506       call zerograd
 507 c      aincr=1.0D-7
 508       print '("Calling CHECK_INT",1pd12.3)',aincr
 509       write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
 510       nf=0
 511       nfl=0
 512       icg=1
 513       call geom_to_var(nvar,x)
 514       call var_to_geom(nvar,x)
 515       call chainbuild
 516       icall=1
 517       print *,'ICG=',ICG
 518       call etotal(energia(0))
 519       etot = energia(0)
 520       call enerprint(energia(0))
 521       print *,'ICG=',ICG
 522 #ifdef MPL
 523       if (MyID.ne.BossID) then
 524         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
 525         nf=x(nvar+1)
 526         nfl=x(nvar+2)
 527         icg=x(nvar+3)
 528       endif
 529 #endif
 530       nf=1
 531       nfl=3
 532 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
 533 c      write (iout,*) "Before gradient"
 534 c      call flush(iout)
 535 #ifdef LBFGS
 536       ff=funcgrad(x,gana)
 537 #else
 538       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
 539 #endif
 540 c      write (iout,*) "After gradient"
 541 c      call flush(iout)
 542 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
 543       icall=1
 544       do i=1,nvar
 545         xi=x(i)
 546         x(i)=xi-0.5D0*aincr
 547         call var_to_geom(nvar,x)
 548         call chainbuild_extconf
 549         call etotal(energia1(0))
 550         etot1=energia1(0)
 551         x(i)=xi+0.5D0*aincr
 552         call var_to_geom(nvar,x)
 553         call chainbuild_extconf
 554         call etotal(energia2(0))
 555         etot2=energia2(0)
 556         gg(i)=(etot2-etot1)/aincr
 557 c        write (iout,*) i,etot1,etot2
 558         x(i)=xi
 559       enddo
 560       write (iout,'(/2a)')' Variable        Numerical       Analytical',
 561      &    '     RelDiff*100% '
 562       do i=1,nvar
 563         if (i.le.nphi) then
 564           ii=i
 565           key = ' phi'
 566         else if (i.le.nphi+ntheta) then
 567           ii=i-nphi
 568           key=' theta'
 569         else if (i.le.nphi+ntheta+nside) then
 570            ii=i-(nphi+ntheta)
 571            key=' alpha'
 572         else
 573            ii=i-(nphi+ntheta+nside)
 574            key=' omega'
 575         endif
 576         write (iout,'(i3,a,i3,3(1pd16.6))')
 577      & i,key,ii,gg(i),gana(i),
 578      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
 579       enddo
 580       return
 581       end