source/unres/src_Eshel/checkder_p.F

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