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) phi(i+1)=beta(i-2,i-1,i,i+1)
 524         omeg(i)=beta(nres+i,i,maxres2,i+1)
 525         alph(i)=alpha(nres+i,i,maxres2)
 526         theta(i+1)=alpha(i-1,i,i+1)
 527         vbld(i)=dist(i-1,i)
 528         vbld_inv(i)=1.0d0/vbld(i)
 529         vbld(nres+i)=dist(nres+i,i)
 530         if (itype(i).ne.10) then
 531           vbld_inv(nres+i)=1.0d0/vbld(nres+i)
 532         else
 533           vbld_inv(nres+i)=0.0d0
 534         endif
 535       enddo
 536 #if defined(PARINT) && defined(MPI)
 537        if (nfgtasks1.gt.1) then
 538 cd       write(iout,*) "iint_start",iint_start," iint_count",
 539 cd     &   (iint_count(i),i=0,nfgtasks-1)," iint_displ",
 540 cd     &   (iint_displ(i),i=0,nfgtasks-1)
 541 cd       write (iout,*) "Gather vbld backbone"
 542 cd       call flush(iout)
 543        time00=MPI_Wtime()
 544        call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
 545      &   MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
 546      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 547 cd       write (iout,*) "Gather vbld_inv"
 548 cd       call flush(iout)
 549        call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
 550      &   MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
 551      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 552 cd       write (iout,*) "Gather vbld side chain"
 553 cd       call flush(iout)
 554        call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
 555      &  MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
 556      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 557 cd       write (iout,*) "Gather vbld_inv side chain"
 558 cd       call flush(iout)
 559        call MPI_Allgatherv(vbld_inv(iint_start+nres),
 560      &   iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
 561      &   iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 562 cd       write (iout,*) "Gather theta"
 563 cd       call flush(iout)
 564        call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
 565      &   MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
 566      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 567 cd       write (iout,*) "Gather phi"
 568 cd       call flush(iout)
 569        call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
 570      &   MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
 571      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 572 #ifdef CRYST_SC
 573 cd       write (iout,*) "Gather alph"
 574 cd       call flush(iout)
 575        call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
 576      &   MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
 577      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 578 cd       write (iout,*) "Gather omeg"
 579 cd       call flush(iout)
 580        call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
 581      &   MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
 582      &   MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
 583 #endif
 584        time_gather=time_gather+MPI_Wtime()-time00
 585       endif
 586 #endif
 587       do i=1,nres-1
 588         do j=1,3
 589           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 590         enddo
 591       enddo
 592       do i=2,nres-1
 593         do j=1,3
 594           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 595         enddo
 596       enddo
 597       if (lprn) then
 598       do i=2,nres
 599        write (iout,1212) restyp(itype(i)),i,vbld(i),
 600      &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
 601      &rad2deg*alph(i),rad2deg*omeg(i)
 602       enddo
 603       endif
 604  1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
 605 #ifdef TIMING
 606       time_intfcart=time_intfcart+MPI_Wtime()-time01
 607 #endif
 608       return
 609       end
 610 c----------------------------------------------------------------------------
 611       subroutine check_eint
 612 C Check the gradient of energy in internal coordinates.
 613       implicit real*8 (a-h,o-z)
 614       include 'DIMENSIONS'
 615       include 'COMMON.CHAIN'
 616       include 'COMMON.DERIV'
 617       include 'COMMON.IOUNITS'
 618       include 'COMMON.VAR'
 619       include 'COMMON.GEO'
 620       common /srutu/ icall
 621       dimension x(maxvar),gana(maxvar),gg(maxvar)
 622       integer uiparm(1)
 623       double precision urparm(1)
 624       double precision energia(0:n_ene),energia1(0:n_ene),
 625      &  energia2(0:n_ene)
 626       character*6 key
 627       external fdum
 628       call zerograd
 629       aincr=1.0D-7
 630       print '(a)','Calling CHECK_INT.'
 631       nf=0
 632       nfl=0
 633       icg=1
 634       call geom_to_var(nvar,x)
 635       call var_to_geom(nvar,x)
 636       call chainbuild
 637       icall=1
 638       print *,'ICG=',ICG
 639       call etotal(energia(0))
 640       etot = energia(0)
 641       call enerprint(energia(0))
 642       print *,'ICG=',ICG
 643 #ifdef MPL
 644       if (MyID.ne.BossID) then
 645         call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
 646         nf=x(nvar+1)
 647         nfl=x(nvar+2)
 648         icg=x(nvar+3)
 649       endif
 650 #endif
 651       nf=1
 652       nfl=3
 653 cd    write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
 654       call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
 655 cd    write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
 656       icall=1
 657       do i=1,nvar
 658         xi=x(i)
 659         x(i)=xi-0.5D0*aincr
 660         call var_to_geom(nvar,x)
 661         call chainbuild
 662         call etotal(energia1(0))
 663         etot1=energia1(0)
 664         x(i)=xi+0.5D0*aincr
 665         call var_to_geom(nvar,x)
 666         call chainbuild
 667         call etotal(energia2(0))
 668         etot2=energia2(0)
 669         gg(i)=(etot2-etot1)/aincr
 670         write (iout,*) i,etot1,etot2
 671         x(i)=xi
 672       enddo
 673       write (iout,'(/2a)')' Variable        Numerical       Analytical',
 674      &    '     RelDiff*100% '
 675       do i=1,nvar
 676         if (i.le.nphi) then
 677           ii=i
 678           key = ' phi'
 679         else if (i.le.nphi+ntheta) then
 680           ii=i-nphi
 681           key=' theta'
 682         else if (i.le.nphi+ntheta+nside) then
 683            ii=i-(nphi+ntheta)
 684            key=' alpha'
 685         else
 686            ii=i-(nphi+ntheta+nside)
 687            key=' omega'
 688         endif
 689         write (iout,'(i3,a,i3,3(1pd16.6))')
 690      & i,key,ii,gg(i),gana(i),
 691      & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)
 692       enddo
 693       return
 694       end