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