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       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 368 #ifdef TIMING
 369       time01=MPI_Wtime()
 370 #endif
 371 #if defined(PARINT) && defined(MPI)
 372       do i=iint_start,iint_end
 373 #else
 374       do i=2,nres
 375 #endif
 376 C        print *,i
 377         dnorm1=dist(i-1,i)
 378         dnorm2=dist(i,i+1)
 379 C        print *,i,dnorm1,dnorm2
 380         do j=1,3
 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)
 383         enddo
 384         be=0.0D0
 385         if (i.gt.2) then
 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)
 389         endif
 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)
 394         endif
 395         if (itype(i).ne.10) then
 396          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
 397         endif
 398         endif
 399         omeg(i)=beta(nres+i,i,maxres2,i+1)
 400 C        print *,omeg(i)
 401         alph(i)=alpha(nres+i,i,maxres2)
 402 C        print *,alph(i)
 403         theta(i+1)=alpha(i-1,i,i+1)
 404         vbld(i)=dist(i-1,i)
 405 C        print *,vbld(i)
 406         vbld_inv(i)=1.0d0/vbld(i)
 407         vbld(nres+i)=dist(nres+i,i)
 408 C            print *,vbld(i+nres)
 409
 410         if (itype(i).ne.10) then
 411           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
 412         else
 413           vbld_inv(nres+i)=0.0d0
 414         endif
 415       enddo
 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"
 422 cd       call flush(iout)
 423        time00=MPI_Wtime()
 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"
 428 cd       call flush(iout)
 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"
 433 cd       call flush(iout)
 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"
 438 cd       call flush(iout)
 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"
 443 cd       call flush(iout)
 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"
 448 cd       call flush(iout)
 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)
 452 #ifdef CRYST_SC
 453 cd       write (iout,*) "Gather alph"
 454 cd       call flush(iout)
 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"
 459 cd       call flush(iout)
 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)
 463 #endif
 464        time_gather=time_gather+MPI_Wtime()-time00
 465       endif
 466 #endif
 467       do i=1,nres-1
 468         do j=1,3
 469           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 470         enddo
 471       enddo
 472       do i=2,nres-1
 473         do j=1,3
 474           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 475         enddo
 476       enddo
 477       if (lprn) then
 478       do i=2,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)
 482       enddo
 483       endif
 484  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
 485 #ifdef TIMING
 486       time_intfcart=time_intfcart+MPI_Wtime()-time01
 487 #endif
 488       return
 489       end
 490 c----------------------------------------------------------------------------
 491       subroutine check_eint
 492 C Check the gradient of energy in internal coordinates.
 493       implicit none
 494       include 'DIMENSIONS'
 495       include 'COMMON.CONTROL'
 496       include 'COMMON.CHAIN'
 497       include 'COMMON.DERIV'
 498       include 'COMMON.IOUNITS'
 499       include 'COMMON.VAR'
 500       include 'COMMON.GEO'
 501       integer icall
 502       common /srutu/ icall
 503       double precision x(maxvar),gana(maxvar),gg(maxvar)
 504       integer uiparm(1)
 505       double precision urparm(1)
 506       double precision energia(0:n_ene),energia1(0:n_ene),
 507      &  energia2(0:n_ene)
 508       character*6 key
 509       double precision fdum
 510       external fdum
 511       double precision funcgrad,ff
 512       external funcgrad
 513       integer i,ii,nf
 514       double precision xi,etot,etot1,etot2
 515       call zerograd
 516 c      aincr=1.0D-7
 517       print '("Calling CHECK_INT",1pd12.3)',aincr
 518       write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
 519       nf=0
 520       nfl=0
 521       icg=1
 522       call geom_to_var(nvar,x)
 523       call var_to_geom(nvar,x)
 524       call chainbuild
 525       icall=1
 526       print *,'ICG=',ICG
 527       call etotal(energia(0))
 528       etot = energia(0)
 529       call enerprint(energia(0))
 530       print *,'ICG=',ICG
 531 #ifdef MPL
 532       if (MyID.ne.BossID) then
 533         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
 534         nf=x(nvar+1)
 535         nfl=x(nvar+2)
 536         icg=x(nvar+3)
 537       endif
 538 #endif
 539       nf=1
 540       nfl=3
 541 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
 542 c      write (iout,*) "Before gradient"
 543 c      call flush(iout)
 544 #ifdef LBFGS
 545       ff=funcgrad(x,gana)
 546 #else
 547       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
 548 #endif
 549 c      write (iout,*) "After gradient"
 550 c      call flush(iout)
 551 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
 552       icall=1
 553       do i=1,nvar
 554         xi=x(i)
 555         x(i)=xi-0.5D0*aincr
 556         call var_to_geom(nvar,x)
 557         call chainbuild_extconf
 558         call etotal(energia1(0))
 559         etot1=energia1(0)
 560         x(i)=xi+0.5D0*aincr
 561         call var_to_geom(nvar,x)
 562         call chainbuild_extconf
 563         call etotal(energia2(0))
 564         etot2=energia2(0)
 565         gg(i)=(etot2-etot1)/aincr
 566 c        write (iout,*) i,etot1,etot2
 567         x(i)=xi
 568       enddo
 569       write (iout,'(/2a)')' Variable        Numerical       Analytical',
 570      &    '     RelDiff*100% '
 571       do i=1,nvar
 572         if (i.le.nphi) then
 573           ii=i
 574           key = ' phi'
 575         else if (i.le.nphi+ntheta) then
 576           ii=i-nphi
 577           key=' theta'
 578         else if (i.le.nphi+ntheta+nside) then
 579            ii=i-(nphi+ntheta)
 580            key=' alpha'
 581         else
 582            ii=i-(nphi+ntheta+nside)
 583            key=' omega'
 584         endif
 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)
 588       enddo
 589       return
 590       end