source/unres/src-5hdiag-tmp/checkder_p.F

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