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