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

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