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       call cartprint
 136       nf=0
 137       icall=0
 138       call geom_to_var(nvar,x)
 139       if (.not.split_ene) then
 140         call etotal(energia(0))
 141         etot=energia(0)
 142         call enerprint(energia(0))
 143 c        write (iout,*) "enter cartgrad"
 144 c        call flush(iout)
 145         call cartgrad
 146 c        write (iout,*) "exit cartgrad"
 147 c        call flush(iout)
 148         icall =1
 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"
 152         do i=1,nres
 153           write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
 154      &      (gxcart(j,i),j=1,3)
 155         enddo
 156         do j=1,3
 157           grad_s(j,0)=gcart(j,0)
 158         enddo
 159         do i=1,nres
 160           do j=1,3
 161             grad_s(j,i)=gcart(j,i)
 162             grad_s(j+3,i)=gxcart(j,i)
 163           enddo
 164         enddo
 165       else
 166 !- split gradient check
 167         call zerograd
 168         call etotal_long(energia(0))
 169         call enerprint(energia(0))
 170         call flush(iout)
 171 c        write (iout,*) "enter cartgrad"
 172 c        call flush(iout)
 173         call cartgrad
 174 c        write (iout,*) "exit cartgrad"
 175 c        call flush(iout)
 176         icall =1
 177         write (iout,*) "longrange grad"
 178         do i=1,nres
 179           write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
 180      &    (gxcart(j,i),j=1,3)
 181         enddo
 182         do j=1,3
 183           grad_s(j,0)=gcart(j,0)
 184         enddo
 185         do i=1,nres
 186           do j=1,3
 187             grad_s(j,i)=gcart(j,i)
 188             grad_s(j+3,i)=gxcart(j,i)
 189           enddo
 190         enddo
 191         call zerograd
 192         call etotal_short(energia(0))
 193         call enerprint(energia(0))
 194         call flush(iout)
 195 c        write (iout,*) "enter cartgrad"
 196 c        call flush(iout)
 197         call cartgrad
 198 c        write (iout,*) "exit cartgrad"
 199 c        call flush(iout)
 200         icall =1
 201         write (iout,*) "shortrange grad"
 202         do i=1,nres
 203           write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
 204      &    (gxcart(j,i),j=1,3)
 205         enddo
 206         do j=1,3
 207           grad_s1(j,0)=gcart(j,0)
 208         enddo
 209         do i=1,nres
 210           do j=1,3
 211             grad_s1(j,i)=gcart(j,i)
 212             grad_s1(j+3,i)=gxcart(j,i)
 213           enddo
 214         enddo
 215       endif
 216 c      write (iout,*) "Vector dc"
 217 c      do i=0,nres
 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)
 220 c      enddo
 221 c      write (iout,*) "Coordinates after chainbuild_cart"
 222 c      call chainbuild_cart
 223 c      call cartprint
 224       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
 225       do i=0,nres
 226         print *,i
 227         do j=1,3
 228           xx(j)=c(j,i+nres)
 229           ddc(j)=dc(j,i)
 230           ddx(j)=dc(j,i+nres)
 231           do k=1,3
 232             dcnorm_safe(k)=dc_norm(k,i)
 233             dxnorm_safe(k)=dc_norm(k,i+nres)
 234           enddo
 235         enddo
 236         do j=1,3
 237           dc(j,i)=ddc(j)+aincr
 238           call chainbuild_cart
 239 #ifdef MPI
 240 c Broadcast the order to compute internal coordinates to the slaves.
 241 c          if (nfgtasks.gt.1)
 242 c     &      call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
 243 #endif
 244 c          call int_from_cart1(.false.)
 245           if (.not.split_ene) then
 246             call etotal(energia1(0))
 247             etot1=energia1(0)
 248           else
 249 !- split gradient
 250             call etotal_long(energia1(0))
 251             etot11=energia1(0)
 252             call etotal_short(energia1(0))
 253             etot12=energia1(0)
 254 c            write (iout,*) "etot11",etot11," etot12",etot12
 255           endif
 256 !- end split gradient
 257 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
 258           dc(j,i)=ddc(j)-aincr
 259           call chainbuild_cart
 260 C          print *,c(j,i)
 261 c          call int_from_cart1(.false.)
 262           if (.not.split_ene) then
 263             call etotal(energia1(0))
 264             etot2=energia1(0)
 265             ggg(j)=(etot1-etot2)/(2*aincr)
 266           else
 267 !- split gradient
 268             call etotal_long(energia1(0))
 269             etot21=energia1(0)
 270             ggg(j)=(etot11-etot21)/(2*aincr)
 271             call etotal_short(energia1(0))
 272             etot22=energia1(0)
 273             ggg1(j)=(etot12-etot22)/(2*aincr)
 274 !- end split gradient
 275 c            write (iout,*) "etot21",etot21," etot22",etot22
 276           endif
 277 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 278           dc(j,i)=ddc(j)
 279           call chainbuild_cart
 280         enddo
 281         do j=1,3
 282           dc(j,i+nres)=ddx(j)+aincr
 283           call chainbuild_cart
 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)
 291 c          write (iout,*)
 292           if (.not.split_ene) then
 293             call etotal(energia1(0))
 294             etot1=energia1(0)
 295           else
 296 !- split gradient
 297             call etotal_long(energia1(0))
 298             etot11=energia1(0)
 299             call etotal_short(energia1(0))
 300             etot12=energia1(0)
 301           endif
 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
 305           call chainbuild_cart
 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)
 309 c          write (iout,*)
 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))
 316             etot2=energia1(0)
 317             ggg(j+3)=(etot1-etot2)/(2*aincr)
 318           else
 319 !- split gradient
 320             call etotal_long(energia1(0))
 321             etot21=energia1(0)
 322             ggg(j+3)=(etot11-etot21)/(2*aincr)
 323             call etotal_short(energia1(0))
 324             etot22=energia1(0)
 325             ggg1(j+3)=(etot12-etot22)/(2*aincr)
 326 !- end split gradient
 327           endif
 328 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 329           dc(j,i+nres)=ddx(j)
 330           call chainbuild_cart
 331         enddo
 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)
 334         if (split_ene) then
 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),
 337      &   k=1,6)
 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)
 341         endif
 342       enddo
 343       return
 344       end
 345 c-------------------------------------------------------------------------
 346       subroutine int_from_cart1(lprn)
 347       implicit none
 348       include 'DIMENSIONS'
 349 #ifdef MPI
 350       include 'mpif.h'
 351       integer ierror
 352 #endif
 353       include 'COMMON.IOUNITS'
 354       include 'COMMON.VAR'
 355       include 'COMMON.CHAIN'
 356       include 'COMMON.GEO'
 357       include 'COMMON.INTERACT'
 358       include 'COMMON.LOCAL'
 359       include 'COMMON.NAMES'
 360       include 'COMMON.SETUP'
 361       include 'COMMON.TIME1'
 362       logical lprn
 363       integer i,j
 364       double precision dnorm1,dnorm2,be
 365       double precision time00
 366       double precision dist,alpha,beta
 367       double precision time01
 368       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 369 #ifdef TIMING
 370       time01=MPI_Wtime()
 371 #endif
 372 #if defined(PARINT) && defined(MPI)
 373       do i=iint_start,iint_end
 374 #else
 375       do i=2,nres
 376 #endif
 377 C        print *,i
 378         dnorm1=dist(i-1,i)
 379         dnorm2=dist(i,i+1)
 380 C        print *,i,dnorm1,dnorm2
 381         do j=1,3
 382           c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
 383      &     +(c(j,i+1)-c(j,i))/dnorm2)
 384         enddo
 385         be=0.0D0
 386         if (i.gt.2) then
 387         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
 388         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
 389          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
 390         endif
 391         if (itype(i-1).ne.10) then
 392          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
 393          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
 394          omicron(2,i)=alpha(i-1+nres,i-1,i)
 395         endif
 396         if (itype(i).ne.10) then
 397          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
 398         endif
 399         endif
 400         omeg(i)=beta(nres+i,i,maxres2,i+1)
 401 C        print *,omeg(i)
 402         alph(i)=alpha(nres+i,i,maxres2)
 403 C        print *,alph(i)
 404         theta(i+1)=alpha(i-1,i,i+1)
 405         vbld(i)=dist(i-1,i)
 406 C        print *,vbld(i)
 407         vbld_inv(i)=1.0d0/vbld(i)
 408         vbld(nres+i)=dist(nres+i,i)
 409 C            print *,vbld(i+nres)
 410
 411         if (itype(i).ne.10) then
 412           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
 413         else
 414           vbld_inv(nres+i)=0.0d0
 415         endif
 416       enddo
 417 #if defined(PARINT) && defined(MPI)
 418        if (nfgtasks1.gt.1) then
 419 cd       write(iout,*) "iint_start",iint_start," iint_count",
 420 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
 421 cd     &   (iint_displ(i),i=0,nfgtasks-1)
 422 cd       write (iout,*) "Gather vbld backbone"
 423 cd       call flush(iout)
 424        time00=MPI_Wtime()
 425        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
 426      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
 427      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 428 cd       write (iout,*) "Gather vbld_inv"
 429 cd       call flush(iout)
 430        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
 431      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
 432      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 433 cd       write (iout,*) "Gather vbld side chain"
 434 cd       call flush(iout)
 435        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
 436      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
 437      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 438 cd       write (iout,*) "Gather vbld_inv side chain"
 439 cd       call flush(iout)
 440        call MPI_Allgatherv(vbld_inv(iint_start+nres),
 441      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
 442      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 443 cd       write (iout,*) "Gather theta"
 444 cd       call flush(iout)
 445        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
 446      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
 447      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 448 cd       write (iout,*) "Gather phi"
 449 cd       call flush(iout)
 450        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
 451      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
 452      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 453 #ifdef CRYST_SC
 454 cd       write (iout,*) "Gather alph"
 455 cd       call flush(iout)
 456        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
 457      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
 458      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 459 cd       write (iout,*) "Gather omeg"
 460 cd       call flush(iout)
 461        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
 462      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
 463      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 464 #endif
 465        time_gather=time_gather+MPI_Wtime()-time00
 466       endif
 467 #endif
 468       do i=1,nres-1
 469         do j=1,3
 470           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 471         enddo
 472       enddo
 473       do i=2,nres-1
 474         do j=1,3
 475           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 476         enddo
 477       enddo
 478       if (lprn) then
 479       do i=2,nres
 480        write (iout,1212) restyp(itype(i)),i,vbld(i),
 481      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
 482      &rad2deg*alph(i),rad2deg*omeg(i)
 483       enddo
 484       endif
 485  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
 486 #ifdef TIMING
 487       time_intfcart=time_intfcart+MPI_Wtime()-time01
 488 #endif
 489       return
 490       end
 491 c----------------------------------------------------------------------------
 492       subroutine check_eint
 493 C Check the gradient of energy in internal coordinates.
 494       implicit none
 495       include 'DIMENSIONS'
 496       include 'COMMON.CONTROL'
 497       include 'COMMON.CHAIN'
 498       include 'COMMON.DERIV'
 499       include 'COMMON.IOUNITS'
 500       include 'COMMON.VAR'
 501       include 'COMMON.GEO'
 502       integer icall
 503       common /srutu/ icall
 504       double precision x(maxvar),gana(maxvar),gg(maxvar)
 505       integer uiparm(1)
 506       double precision urparm(1)
 507       double precision energia(0:n_ene),energia1(0:n_ene),
 508      &  energia2(0:n_ene)
 509       character*6 key
 510       double precision fdum
 511       external fdum
 512       double precision funcgrad,ff
 513       external funcgrad
 514       integer i,ii,nf
 515       double precision xi,etot,etot1,etot2
 516       call zerograd
 517 c      aincr=1.0D-7
 518       print '("Calling CHECK_INT",1pd12.3)',aincr
 519       write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
 520       nf=0
 521       nfl=0
 522       icg=1
 523       call geom_to_var(nvar,x)
 524       call var_to_geom(nvar,x)
 525       call chainbuild
 526       icall=1
 527       print *,'ICG=',ICG
 528       call etotal(energia(0))
 529       etot = energia(0)
 530       call enerprint(energia(0))
 531       print *,'ICG=',ICG
 532 #ifdef MPL
 533       if (MyID.ne.BossID) then
 534         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
 535         nf=x(nvar+1)
 536         nfl=x(nvar+2)
 537         icg=x(nvar+3)
 538       endif
 539 #endif
 540       nf=1
 541       nfl=3
 542 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
 543 c      write (iout,*) "Before gradient"
 544 c      call flush(iout)
 545 #ifdef LBFGS
 546       ff=funcgrad(x,gana)
 547 #else
 548       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
 549 #endif
 550 c      write (iout,*) "After gradient"
 551 c      call flush(iout)
 552 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
 553       icall=1
 554       do i=1,nvar
 555         xi=x(i)
 556         x(i)=xi-0.5D0*aincr
 557         call var_to_geom(nvar,x)
 558         call chainbuild_extconf
 559         call etotal(energia1(0))
 560         etot1=energia1(0)
 561         x(i)=xi+0.5D0*aincr
 562         call var_to_geom(nvar,x)
 563         call chainbuild_extconf
 564         call etotal(energia2(0))
 565         etot2=energia2(0)
 566         gg(i)=(etot2-etot1)/aincr
 567 c        write (iout,*) i,etot1,etot2
 568         x(i)=xi
 569       enddo
 570       write (iout,'(/2a)')' Variable        Numerical       Analytical',
 571      &    '     RelDiff*100% '
 572       do i=1,nvar
 573         if (i.le.nphi) then
 574           ii=i
 575           key = ' phi'
 576         else if (i.le.nphi+ntheta) then
 577           ii=i-nphi
 578           key=' theta'
 579         else if (i.le.nphi+ntheta+nside) then
 580            ii=i-(nphi+ntheta)
 581            key=' alpha'
 582         else
 583            ii=i-(nphi+ntheta+nside)
 584            key=' omega'
 585         endif
 586         write (iout,'(i3,a,i3,3(1pd16.6))')
 587      & i,key,ii,gg(i),gana(i),
 588      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
 589       enddo
 590       return
 591       end