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