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