source/unres/src-HCD-5D/checkder_p.F

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