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