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 #ifdef MOMENT
 184       include 'COMMON.CONTACTS.MOMENT'
 185 #endif
 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       aincr=1.0D-7
 198       print '(a)','CG processor',me,' calling CHECK_CART.'
 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 #ifdef MOMENT
 265       include 'COMMON.CONTACTS.MOMENT'
 266 #endif
 267       include 'COMMON.MD_'
 268       include 'COMMON.LOCAL'
 269       include 'COMMON.SPLITELE'
 270       common /srutu/ icall
 271       dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
 272      &  g(maxvar)
 273       dimension dcnorm_safe(3),dxnorm_safe(3)
 274       dimension grad_s(6,0:maxres),grad_s1(6,0:maxres)
 275       double precision phi_temp(maxres),theta_temp(maxres),
 276      &  alph_temp(maxres),omeg_temp(maxres)
 277       double precision energia(0:n_ene),energia1(0:n_ene)
 278       integer uiparm(1)
 279       double precision urparm(1)
 280       external fdum
 281       r_cut=2.0d0
 282       rlambd=0.3d0
 283       icg=1
 284       nf=0
 285       nfl=0
 286       call intout
 287 c      call intcartderiv
 288 c      call checkintcartgrad
 289       call zerograd
 290       aincr=1.0D-5
 291       write(iout,*) 'Calling CHECK_ECARTINT.'
 292       nf=0
 293       icall=0
 294       call geom_to_var(nvar,x)
 295       if (.not.split_ene) then
 296         call etotal(energia(0))
 297         etot=energia(0)
 298         call enerprint(energia(0))
 299         call flush(iout)
 300         write (iout,*) "enter cartgrad"
 301         call flush(iout)
 302         call cartgrad
 303         write (iout,*) "exit cartgrad"
 304         call flush(iout)
 305         icall =1
 306         do i=1,nres
 307           write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
 308         enddo
 309         do j=1,3
 310           grad_s(j,0)=gcart(j,0)
 311         enddo
 312         do i=1,nres
 313           do j=1,3
 314             grad_s(j,i)=gcart(j,i)
 315             grad_s(j+3,i)=gxcart(j,i)
 316           enddo
 317         enddo
 318       else
 319 !- split gradient check
 320          write (iout,*) "split_ene not supported"
 321 c        call zerograd
 322 c        call etotal_long(energia(0))
 323 c        call enerprint(energia(0))
 324 c        call flush(iout)
 325 c        write (iout,*) "enter cartgrad"
 326 c        call flush(iout)
 327 c        call cartgrad
 328 c        write (iout,*) "exit cartgrad"
 329 c        call flush(iout)
 330 c        icall =1
 331 c        write (iout,*) "longrange grad"
 332 c        do i=1,nres
 333 c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 334 c     &    (gxcart(j,i),j=1,3)
 335 c        enddo
 336 c        do j=1,3
 337 c          grad_s(j,0)=gcart(j,0)
 338 c        enddo
 339 c        do i=1,nres
 340 c          do j=1,3
 341 c            grad_s(j,i)=gcart(j,i)
 342 c            grad_s(j+3,i)=gxcart(j,i)
 343 c          enddo
 344 c        enddo
 345 c        call zerograd
 346 c        call etotal_short(energia(0))
 347 c        call enerprint(energia(0))
 348 c        call flush(iout)
 349 c        write (iout,*) "enter cartgrad"
 350 c        call flush(iout)
 351 c        call cartgrad
 352 c        write (iout,*) "exit cartgrad"
 353 c        call flush(iout)
 354 c        icall =1
 355 c        write (iout,*) "shortrange grad"
 356 c        do i=1,nres
 357 c          write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 358 c     &    (gxcart(j,i),j=1,3)
 359 c        enddo
 360 c        do j=1,3
 361 c          grad_s1(j,0)=gcart(j,0)
 362 c        enddo
 363 c        do i=1,nres
 364 c          do j=1,3
 365 c            grad_s1(j,i)=gcart(j,i)
 366 c            grad_s1(j+3,i)=gxcart(j,i)
 367 c          enddo
 368 c        enddo
 369       endif
 370       write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
 371       do i=0,nres
 372         do j=1,3
 373           xx(j)=c(j,i+nres)
 374           ddc(j)=dc(j,i)
 375           ddx(j)=dc(j,i+nres)
 376           do k=1,3
 377             dcnorm_safe(k)=dc_norm(k,i)
 378             dxnorm_safe(k)=dc_norm(k,i+nres)
 379           enddo
 380         enddo
 381         do j=1,3
 382           dc(j,i)=ddc(j)+aincr
 383           call chainbuild_cart
 384 #ifdef MPI
 385 c Broadcast the order to compute internal coordinates to the slaves.
 386 c          if (nfgtasks.gt.1)
 387 c     &      call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
 388 #endif
 389 c          call int_from_cart1(.false.)
 390           if (.not.split_ene) then
 391             call etotal(energia1(0))
 392             etot1=energia1(0)
 393           else
 394 !- split gradient
 395 c            call etotal_long(energia1(0))
 396 c            etot11=energia1(0)
 397 c            call etotal_short(energia1(0))
 398 c            etot12=energia1(0)
 399 c            write (iout,*) "etot11",etot11," etot12",etot12
 400           endif
 401 !- end split gradient
 402 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
 403           dc(j,i)=ddc(j)-aincr
 404           call chainbuild_cart
 405 c          call int_from_cart1(.false.)
 406           if (.not.split_ene) then
 407             call etotal(energia1(0))
 408             etot2=energia1(0)
 409             ggg(j)=(etot1-etot2)/(2*aincr)
 410           else
 411 !- split gradient
 412 c            call etotal_long(energia1(0))
 413 c            etot21=energia1(0)
 414 c           ggg(j)=(etot11-etot21)/(2*aincr)
 415 c            call etotal_short(energia1(0))
 416 c            etot22=energia1(0)
 417 c           ggg1(j)=(etot12-etot22)/(2*aincr)
 418 !- end split gradient
 419 c            write (iout,*) "etot21",etot21," etot22",etot22
 420           endif
 421 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 422           dc(j,i)=ddc(j)
 423           call chainbuild_cart
 424         enddo
 425         do j=1,3
 426           dc(j,i+nres)=ddx(j)+aincr
 427           call chainbuild_cart
 428 c          write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
 429 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
 430 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
 431 c          write (iout,*) "dxnormnorm",dsqrt(
 432 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
 433 c          write (iout,*) "dxnormnormsafe",dsqrt(
 434 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
 435 c          write (iout,*)
 436           if (.not.split_ene) then
 437             call etotal(energia1(0))
 438             etot1=energia1(0)
 439           else
 440 !- split gradient
 441 c            call etotal_long(energia1(0))
 442 c            etot11=energia1(0)
 443 c            call etotal_short(energia1(0))
 444 c            etot12=energia1(0)
 445           endif
 446 !- end split gradient
 447 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
 448           dc(j,i+nres)=ddx(j)-aincr
 449           call chainbuild_cart
 450 c          write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
 451 c          write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
 452 c          write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
 453 c          write (iout,*)
 454 c          write (iout,*) "dxnormnorm",dsqrt(
 455 c     &  dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
 456 c          write (iout,*) "dxnormnormsafe",dsqrt(
 457 c     &      dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
 458           if (.not.split_ene) then
 459             call etotal(energia1(0))
 460             etot2=energia1(0)
 461             ggg(j+3)=(etot1-etot2)/(2*aincr)
 462           else
 463 !- split gradient
 464 c            call etotal_long(energia1(0))
 465 c            etot21=energia1(0)
 466 c           ggg(j+3)=(etot11-etot21)/(2*aincr)
 467 c           call etotal_short(energia1(0))
 468 c            etot22=energia1(0)
 469 c           ggg1(j+3)=(etot12-etot22)/(2*aincr)
 470 !- end split gradient
 471           endif
 472 c          write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
 473           dc(j,i+nres)=ddx(j)
 474           call chainbuild_cart
 475         enddo
 476         write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 477      &   i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
 478         if (split_ene) then
 479           write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 480      &   i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
 481      &   k=1,6)
 482          write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
 483      &   i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
 484      &   ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
 485         endif
 486       enddo
 487       return
 488       end
 489 c-------------------------------------------------------------------------
 490       subroutine int_from_cart1(lprn)
 491       implicit real*8 (a-h,o-z)
 492       include 'DIMENSIONS'
 493 #ifdef MPI
 494       include 'mpif.h'
 495       integer ierror
 496 #endif
 497       include 'COMMON.IOUNITS'
 498       include 'COMMON.VAR'
 499       include 'COMMON.CHAIN'
 500       include 'COMMON.GEO'
 501       include 'COMMON.INTERACT'
 502       include 'COMMON.LOCAL'
 503       include 'COMMON.NAMES'
 504       include 'COMMON.SETUP'
 505       include 'COMMON.TIME1'
 506       logical lprn
 507       if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
 508 #ifdef TIMING
 509       time01=MPI_Wtime()
 510 #endif
 511 #if defined(PARINT) && defined(MPI)
 512       do i=iint_start,iint_end
 513 #else
 514       do i=2,nres
 515 #endif
 516         dnorm1=dist(i-1,i)
 517         dnorm2=dist(i,i+1)
 518         do j=1,3
 519           c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
 520      &     +(c(j,i+1)-c(j,i))/dnorm2)
 521         enddo
 522         be=0.0D0
 523         if (i.gt.2) then
 524         if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
 525         if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
 526          tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
 527         endif
 528         if (itype(i-1).ne.10) then
 529          tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
 530          omicron(1,i)=alpha(i-2,i-1,i-1+nres)
 531          omicron(2,i)=alpha(i-1+nres,i-1,i)
 532         endif
 533         if (itype(i).ne.10) then
 534          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
 535         endif
 536         endif
 537         omeg(i)=beta(nres+i,i,maxres2,i+1)
 538         alph(i)=alpha(nres+i,i,maxres2)
 539         theta(i+1)=alpha(i-1,i,i+1)
 540         vbld(i)=dist(i-1,i)
 541         vbld_inv(i)=1.0d0/vbld(i)
 542         vbld(nres+i)=dist(nres+i,i)
 543         if (itype(i).ne.10) then
 544           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
 545         else
 546           vbld_inv(nres+i)=0.0d0
 547         endif
 548       enddo
 549 #if defined(PARINT) && defined(MPI)
 550        if (nfgtasks1.gt.1) then
 551 cd       write(iout,*) "iint_start",iint_start," iint_count",
 552 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
 553 cd     &   (iint_displ(i),i=0,nfgtasks-1)
 554 cd       write (iout,*) "Gather vbld backbone"
 555 cd       call flush(iout)
 556        time00=MPI_Wtime()
 557        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
 558      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
 559      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 560 cd       write (iout,*) "Gather vbld_inv"
 561 cd       call flush(iout)
 562        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
 563      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
 564      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 565 cd       write (iout,*) "Gather vbld side chain"
 566 cd       call flush(iout)
 567        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
 568      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
 569      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 570 cd       write (iout,*) "Gather vbld_inv side chain"
 571 cd       call flush(iout)
 572        call MPI_Allgatherv(vbld_inv(iint_start+nres),
 573      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
 574      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 575 cd       write (iout,*) "Gather theta"
 576 cd       call flush(iout)
 577        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
 578      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
 579      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 580 cd       write (iout,*) "Gather phi"
 581 cd       call flush(iout)
 582        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
 583      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
 584      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 585 #ifdef CRYST_SC
 586 cd       write (iout,*) "Gather alph"
 587 cd       call flush(iout)
 588        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
 589      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
 590      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 591 cd       write (iout,*) "Gather omeg"
 592 cd       call flush(iout)
 593        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
 594      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
 595      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 596 #endif
 597        time_gather=time_gather+MPI_Wtime()-time00
 598       endif
 599 #endif
 600       do i=1,nres-1
 601         do j=1,3
 602           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 603         enddo
 604       enddo
 605       do i=2,nres-1
 606         do j=1,3
 607           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 608         enddo
 609       enddo
 610       if (lprn) then
 611       do i=2,nres
 612        write (iout,1212) restyp(itype(i)),i,vbld(i),
 613      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
 614      &rad2deg*alph(i),rad2deg*omeg(i)
 615       enddo
 616       endif
 617  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
 618 #ifdef TIMING
 619       time_intfcart=time_intfcart+MPI_Wtime()-time01
 620 #endif
 621       return
 622       end
 623 c----------------------------------------------------------------------------
 624       subroutine check_eint
 625 C Check the gradient of energy in internal coordinates.
 626       implicit real*8 (a-h,o-z)
 627       include 'DIMENSIONS'
 628       include 'COMMON.CHAIN'
 629       include 'COMMON.DERIV'
 630       include 'COMMON.IOUNITS'
 631       include 'COMMON.VAR'
 632       include 'COMMON.GEO'
 633       common /srutu/ icall
 634       dimension x(maxvar),gana(maxvar),gg(maxvar)
 635       integer uiparm(1)
 636       double precision urparm(1)
 637       double precision energia(0:n_ene),energia1(0:n_ene),
 638      &  energia2(0:n_ene)
 639       character*6 key
 640       external fdum
 641       call zerograd
 642       aincr=1.0D-7
 643       print '(a)','Calling CHECK_INT.'
 644       nf=0
 645       nfl=0
 646       icg=1
 647       call geom_to_var(nvar,x)
 648       call var_to_geom(nvar,x)
 649       call chainbuild
 650       icall=1
 651       print *,'ICG=',ICG
 652       call etotal(energia(0))
 653       etot = energia(0)
 654       call enerprint(energia(0))
 655       print *,'ICG=',ICG
 656 #ifdef MPL
 657       if (MyID.ne.BossID) then
 658         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
 659         nf=x(nvar+1)
 660         nfl=x(nvar+2)
 661         icg=x(nvar+3)
 662       endif
 663 #endif
 664       nf=1
 665       nfl=3
 666 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
 667       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
 668 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
 669       icall=1
 670       do i=1,nvar
 671         xi=x(i)
 672         x(i)=xi-0.5D0*aincr
 673         call var_to_geom(nvar,x)
 674         call chainbuild
 675         call etotal(energia1(0))
 676         etot1=energia1(0)
 677         x(i)=xi+0.5D0*aincr
 678         call var_to_geom(nvar,x)
 679         call chainbuild
 680         call etotal(energia2(0))
 681         etot2=energia2(0)
 682         gg(i)=(etot2-etot1)/aincr
 683         write (iout,*) i,etot1,etot2
 684         x(i)=xi
 685       enddo
 686       write (iout,'(/2a)')' Variable        Numerical       Analytical',
 687      &    '     RelDiff*100% '
 688       do i=1,nvar
 689         if (i.le.nphi) then
 690           ii=i
 691           key = ' phi'
 692         else if (i.le.nphi+ntheta) then
 693           ii=i-nphi
 694           key=' theta'
 695         else if (i.le.nphi+ntheta+nside) then
 696            ii=i-(nphi+ntheta)
 697            key=' alpha'
 698         else
 699            ii=i-(nphi+ntheta+nside)
 700            key=' omega'
 701         endif
 702         write (iout,'(i3,a,i3,3(1pd16.6))')
 703      & i,key,ii,gg(i),gana(i),
 704      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
 705       enddo
 706       return
 707       end