source/unres/src-HCD-5D/gradient_p.F.org.debug

   1       subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'COMMON.CONTROL'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.DERIV'
   7       include 'COMMON.VAR'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.FFIELD'
  10       include 'COMMON.MD'
  11       include 'COMMON.QRESTR'
  12       include 'COMMON.IOUNITS'
  13       integer n,nf
  14       double precision ufparm
  15       external ufparm
  16       integer uiparm(1)
  17       double precision urparm(1)
  18       double precision x(n),g(n)
  19       integer i,j,k,ind,ind1
  20       double precision f,gthetai,gphii,galphai,gomegai
  21 c
  22 c This subroutine calculates total internal coordinate gradient.
  23 c Depending on the number of function evaluations, either whole energy
  24 c is evaluated beforehand, Cartesian coordinates and their derivatives in
  25 c internal coordinates are reevaluated or only the cartesian-in-internal
  26 c coordinate derivatives are evaluated. The subroutine was designed to work
  27 c with SUMSL.
  28 c
  29 c
  30       icg=mod(nf,2)+1
  31
  32 cd      print *,'grad',nf,icg
  33       if (nf-nfl+1) 20,30,40
  34    20 call func(n,x,nf,f,uiparm,urparm,ufparm)
  35 c     write (iout,*) 'grad 20'
  36       if (nf.eq.0) return
  37       goto 40
  38    30 call var_to_geom(n,x)
  39       call chainbuild
  40 c     write (iout,*) 'grad 30'
  41 C
  42 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
  43 C
  44    40 call cartder
  45 c     write (iout,*) 'grad 40'
  46 c     print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon
  47 C
  48 C Convert the Cartesian gradient into internal-coordinate gradient.
  49 C
  50       ind=0
  51       ind1=0
  52       do i=1,nres-2
  53         gthetai=0.0D0
  54         gphii=0.0D0
  55         do j=i+1,nres-1
  56           ind=ind+1
  57 c         ind=indmat(i,j)
  58 c         print *,'GRAD: i=',i,' jc=',j,' ind=',ind
  59           do k=1,3
  60             gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
  61           enddo
  62           do k=1,3
  63             gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
  64           enddo
  65         enddo
  66         do j=i+1,nres-1
  67           ind1=ind1+1
  68 c         ind1=indmat(i,j)
  69 c         print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1
  70           write (iout,*) "i",i," j",j," ind1",ind1
  71           write (iout,*) "dxdv",(dxdv(k,ind1),k=1,6)
  72           write (iout,*) "gradx",(gradx(k,j,icg),k=1,3)
  73           do k=1,3
  74             gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
  75             gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
  76           enddo
  77         enddo
  78         if (i.gt.1) g(i-1)=gphii
  79         if (n.gt.nphi) g(nphi+i)=gthetai
  80       enddo
  81       if (n.le.nphi+ntheta) goto 10
  82       do i=2,nres-1
  83         if (itype(i).ne.10) then
  84           galphai=0.0D0
  85           gomegai=0.0D0
  86           do k=1,3
  87             galphai=galphai+dxds(k,i)*gradx(k,i,icg)
  88           enddo
  89           do k=1,3
  90             gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
  91           enddo
  92           g(ialph(i,1))=galphai
  93           g(ialph(i,1)+nside)=gomegai
  94         endif
  95       enddo
  96 C
  97 C Add the components corresponding to local energy terms.
  98 C
  99    10 continue
 100 c Add the usampl contributions
 101       if (usampl) then
 102          do i=1,nres-3
 103            gloc(i,icg)=gloc(i,icg)+dugamma(i)
 104          enddo
 105          do i=1,nres-2
 106            gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
 107          enddo
 108       endif
 109       do i=1,nvar
 110 cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
 111         g(i)=g(i)+gloc(i,icg)
 112       enddo
 113 C Uncomment following three lines for diagnostics.
 114 cd    call intout
 115 cd    call briefout(0,0.0d0)
 116 cd    write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
 117       return
 118       end
 119 C-------------------------------------------------------------------------
 120       subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
 121       implicit none
 122       include 'DIMENSIONS'
 123       include 'COMMON.CHAIN'
 124       include 'COMMON.DERIV'
 125       include 'COMMON.VAR'
 126       include 'COMMON.INTERACT'
 127       include 'COMMON.FFIELD'
 128       include 'COMMON.IOUNITS'
 129       integer n,nf
 130       double precision ufparm
 131       external ufparm
 132       integer uiparm(1)
 133       double precision urparm(1)
 134       double precision x(maxvar),g(maxvar)
 135       integer i,j,k,ig,ind,ij,igall
 136       double precision f,gthetai,gphii,galphai,gomegai
 137
 138       icg=mod(nf,2)+1
 139       if (nf-nfl+1) 20,30,40
 140    20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
 141 c     write (iout,*) 'grad 20'
 142       if (nf.eq.0) return
 143       goto 40
 144    30 continue
 145 #ifdef OSF
 146 c     Intercept NaNs in the coordinates
 147 c      write(iout,*) (var(i),i=1,nvar)
 148       x_sum=0.D0
 149       do i=1,n
 150         x_sum=x_sum+x(i)
 151       enddo
 152       if (x_sum.ne.x_sum) then
 153         write(iout,*)" *** grad_restr : Found NaN in coordinates"
 154         call flush(iout)
 155         print *," *** grad_restr : Found NaN in coordinates"
 156         return
 157       endif
 158 #endif
 159       call var_to_geom_restr(n,x)
 160       call chainbuild
 161 C
 162 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
 163 C
 164    40 call cartder
 165 C
 166 C Convert the Cartesian gradient into internal-coordinate gradient.
 167 C
 168
 169       ig=0
 170       ind=nres-2
 171       do i=2,nres-2
 172        IF (mask_phi(i+2).eq.1) THEN
 173         gphii=0.0D0
 174         do j=i+1,nres-1
 175           ind=ind+1
 176           do k=1,3
 177             gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
 178             gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
 179           enddo
 180         enddo
 181         ig=ig+1
 182         g(ig)=gphii
 183        ELSE
 184         ind=ind+nres-1-i
 185        ENDIF
 186       enddo
 187
 188
 189       ind=0
 190       do i=1,nres-2
 191        IF (mask_theta(i+2).eq.1) THEN
 192         ig=ig+1
 193         gthetai=0.0D0
 194         do j=i+1,nres-1
 195           ind=ind+1
 196           do k=1,3
 197             gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
 198             gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
 199           enddo
 200         enddo
 201         g(ig)=gthetai
 202        ELSE
 203         ind=ind+nres-1-i
 204        ENDIF
 205       enddo
 206
 207       do i=2,nres-1
 208         if (itype(i).ne.10) then
 209          IF (mask_side(i).eq.1) THEN
 210           ig=ig+1
 211           galphai=0.0D0
 212           do k=1,3
 213             galphai=galphai+dxds(k,i)*gradx(k,i,icg)
 214           enddo
 215           g(ig)=galphai
 216          ENDIF
 217         endif
 218       enddo
 219
 220
 221       do i=2,nres-1
 222         if (itype(i).ne.10) then
 223          IF (mask_side(i).eq.1) THEN
 224           ig=ig+1
 225           gomegai=0.0D0
 226           do k=1,3
 227             gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
 228           enddo
 229           g(ig)=gomegai
 230          ENDIF
 231         endif
 232       enddo
 233
 234 C
 235 C Add the components corresponding to local energy terms.
 236 C
 237
 238       ig=0
 239       igall=0
 240       do i=4,nres
 241         igall=igall+1
 242         if (mask_phi(i).eq.1) then
 243           ig=ig+1
 244           g(ig)=g(ig)+gloc(igall,icg)
 245         endif
 246       enddo
 247
 248       do i=3,nres
 249         igall=igall+1
 250         if (mask_theta(i).eq.1) then
 251           ig=ig+1
 252           g(ig)=g(ig)+gloc(igall,icg)
 253         endif
 254       enddo
 255
 256       do ij=1,2
 257       do i=2,nres-1
 258         if (itype(i).ne.10) then
 259           igall=igall+1
 260           if (mask_side(i).eq.1) then
 261             ig=ig+1
 262             g(ig)=g(ig)+gloc(igall,icg)
 263           endif
 264         endif
 265       enddo
 266       enddo
 267
 268 cd      do i=1,ig
 269 cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
 270 cd      enddo
 271       return
 272       end
 273 C-------------------------------------------------------------------------
 274       subroutine cartgrad
 275       implicit none
 276       include 'DIMENSIONS'
 277 #ifdef MPI
 278       include 'mpif.h'
 279 #endif
 280       include 'COMMON.CONTROL'
 281       include 'COMMON.CHAIN'
 282       include 'COMMON.DERIV'
 283       include 'COMMON.VAR'
 284       include 'COMMON.INTERACT'
 285       include 'COMMON.FFIELD'
 286       include 'COMMON.MD'
 287       include 'COMMON.QRESTR'
 288       include 'COMMON.IOUNITS'
 289       include 'COMMON.TIME1'
 290       integer i,j,kk
 291 c
 292 c This subrouting calculates total Cartesian coordinate gradient.
 293 c The subroutine chainbuild_cart and energy MUST be called beforehand.
 294 c
 295 #ifdef TIMING
 296       time00=MPI_Wtime()
 297 #endif
 298       icg=1
 299 #ifdef DEBUG
 300       write (iout,*) "Before sum_gradient"
 301       do i=1,nres-1
 302         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
 303         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
 304       enddo
 305       write (iout,*) "gsaxsc, gsaxcx"
 306       do i=1,nres-1
 307         write (iout,*) i," gsaxsc  ",(gsaxsc(j,i),j=1,3)
 308         write (iout,*) i," gsaxsx  ",(gsaxsx(j,i),j=1,3)
 309       enddo
 310 #endif
 311       call sum_gradient
 312 #ifdef TIMING
 313 #endif
 314 #ifdef DEBUG
 315       write (iout,*) "After sum_gradient"
 316       do i=1,nres-1
 317         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
 318         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
 319       enddo
 320 #endif
 321 c If performing constraint dynamics, add the gradients of the constraint energy
 322       if(usampl.and.totT.gt.eq_time) then
 323 #ifdef DEBUG
 324          write (iout,*) "dudconst, duscdiff, dugamma,dutheta"
 325          write (iout,*) "wumb",wumb
 326          do i=1,nct
 327            write (iout,'(i5,3f10.5,5x,3f10.5,5x,2f10.5)')
 328      &      i,(dudconst(j,i),j=1,3),(duscdiff(j,i),j=1,3),
 329      &      dugamma(i),dutheta(i)
 330          enddo
 331 #endif
 332          do i=1,nct
 333            do j=1,3
 334              gradc(j,i,icg)=gradc(j,i,icg)+
 335      &          wumb*(dudconst(j,i)+duscdiff(j,i))
 336              gradx(j,i,icg)=gradx(j,i,icg)+
 337      &          wumb*(dudxconst(j,i)+duscdiffx(j,i))
 338            enddo
 339          enddo
 340          do i=1,nres-3
 341            gloc(i,icg)=gloc(i,icg)+wumb*dugamma(i)
 342          enddo
 343          do i=1,nres-2
 344            gloc(nphi+i,icg)=gloc(nphi+i,icg)+wumb*dutheta(i)
 345          enddo
 346       endif
 347 #ifdef TIMING
 348       time01=MPI_Wtime()
 349 #endif
 350       call intcartderiv
 351 #ifdef TIMING
 352       time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
 353 #endif
 354 cd      call checkintcartgrad
 355 cd      write(iout,*) 'calling int_to_cart'
 356 #ifdef DEBUG
 357       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 358 #endif
 359       do i=1,nct
 360         do j=1,3
 361           gcart(j,i)=gradc(j,i,icg)
 362           gxcart(j,i)=gradx(j,i,icg)
 363         enddo
 364 #ifdef DEBUG
 365         if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 366         write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
 367      &    (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg),
 368      &    gloc(ialph(i,1),icg),gloc(ialph(i,1)+nside,icg)
 369         else
 370         write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
 371      &    (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg)
 372         endif
 373         call flush(iout)
 374 #endif
 375       enddo
 376 #ifdef TIMING
 377       time01=MPI_Wtime()
 378 #endif
 379       call int_to_cart
 380 #ifdef TIMING
 381       time_inttocart=time_inttocart+MPI_Wtime()-time01
 382 #endif
 383 #ifdef DEBUG
 384       write (iout,*) "gcart and gxcart after int_to_cart"
 385       do i=0,nres-1
 386         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 387      &      (gxcart(j,i),j=1,3)
 388       enddo
 389 #endif
 390 #ifdef TIMING
 391       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 392 #endif
 393       return
 394       end
 395 c---------------------------------------------------------------------------
 396 #ifdef FIVEDIAG
 397       subroutine grad_transform
 398       implicit none
 399       include 'DIMENSIONS'
 400 #ifdef MPI
 401       include 'mpif.h'
 402 #endif
 403       include 'COMMON.CONTROL'
 404       include 'COMMON.CHAIN'
 405       include 'COMMON.DERIV'
 406       include 'COMMON.VAR'
 407       include 'COMMON.INTERACT'
 408       include 'COMMON.FFIELD'
 409       include 'COMMON.MD'
 410       include 'COMMON.QRESTR'
 411       include 'COMMON.IOUNITS'
 412       include 'COMMON.TIME1'
 413       integer i,j,kk
 414 #ifdef DEBUG
 415       write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
 416 #endif
 417       do i=nres,1,-1
 418         do j=1,3
 419           gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 420 !          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 421         enddo
 422 !        write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
 423 !            (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
 424       enddo
 425 ! Correction: dummy residues
 426       if (nnt.gt.1) then
 427         do j=1,3
 428           gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
 429         enddo
 430       endif
 431       if (nct.lt.nres) then
 432         do j=1,3
 433 !          gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
 434           gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
 435         enddo
 436       endif
 437 #ifdef DEBUG
 438       write (iout,*) "CA/SC gradient"
 439       do i=1,nres
 440         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 441      &      (gxcart(j,i),j=1,3)
 442       enddo
 443 #endif
 444       return
 445       end
 446 #endif
 447 C-------------------------------------------------------------------------
 448       subroutine zerograd
 449       implicit none
 450       include 'DIMENSIONS'
 451       include 'COMMON.DERIV'
 452       include 'COMMON.CHAIN'
 453       include 'COMMON.VAR'
 454       include 'COMMON.MD'
 455       include 'COMMON.SCCOR'
 456       include 'COMMON.SHIELD'
 457       integer i,j,kk,intertyp,maxshieldlist
 458       maxshieldlist=0
 459 C
 460 C Initialize Cartesian-coordinate gradient
 461 C
 462       do i=-1,nres
 463         do j=1,3
 464           gvdwx(j,i)=0.0D0
 465           gradx_scp(j,i)=0.0D0
 466           gvdwc(j,i)=0.0D0
 467           gvdwc_scp(j,i)=0.0D0
 468           gvdwc_scpp(j,i)=0.0d0
 469           gelc (j,i)=0.0D0
 470 C below is zero grad for shielding in order: ees (p-p)
 471 C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
 472           gshieldx(j,i)=0.0d0
 473           gshieldc(j,i)=0.0d0
 474           gshieldc_loc(j,i)=0.0d0
 475           gshieldx_ec(j,i)=0.0d0
 476           gshieldc_ec(j,i)=0.0d0
 477           gshieldc_loc_ec(j,i)=0.0d0
 478           gshieldx_t3(j,i)=0.0d0
 479           gshieldc_t3(j,i)=0.0d0
 480           gshieldc_loc_t3(j,i)=0.0d0
 481           gshieldx_t4(j,i)=0.0d0
 482           gshieldc_t4(j,i)=0.0d0
 483           gshieldc_loc_t4(j,i)=0.0d0
 484           gshieldx_ll(j,i)=0.0d0
 485           gshieldc_ll(j,i)=0.0d0
 486           gshieldc_loc_ll(j,i)=0.0d0
 487 C end of zero grad for shielding
 488           gelc_long(j,i)=0.0D0
 489           gradb(j,i)=0.0d0
 490           gradbx(j,i)=0.0d0
 491           gvdwpp(j,i)=0.0d0
 492           gel_loc(j,i)=0.0d0
 493           gel_loc_long(j,i)=0.0d0
 494           ghpbc(j,i)=0.0D0
 495           ghpbx(j,i)=0.0D0
 496           gsaxsc(j,i)=0.0D0
 497           gsaxsx(j,i)=0.0D0
 498           gcorr3_turn(j,i)=0.0d0
 499           gcorr4_turn(j,i)=0.0d0
 500           gradcorr(j,i)=0.0d0
 501           gradcorr_long(j,i)=0.0d0
 502           gradcorr5_long(j,i)=0.0d0
 503           gradcorr6_long(j,i)=0.0d0
 504           gcorr6_turn_long(j,i)=0.0d0
 505           gradcorr5(j,i)=0.0d0
 506           gradcorr6(j,i)=0.0d0
 507           gcorr6_turn(j,i)=0.0d0
 508           gsccorc(j,i)=0.0d0
 509           gsccorx(j,i)=0.0d0
 510           gradc(j,i,icg)=0.0d0
 511           gradx(j,i,icg)=0.0d0
 512           gscloc(j,i)=0.0d0
 513           gsclocx(j,i)=0.0d0
 514           gliptranc(j,i)=0.0d0
 515           gliptranx(j,i)=0.0d0
 516           gradafm(j,i)=0.0d0
 517           grad_shield(j,i)=0.0d0
 518           gg_tube(j,i)=0.0d0
 519           gg_tube_sc(j,i)=0.0d0
 520 C grad_shield_side is Cbeta sidechain gradient
 521           do kk=1,maxshieldlist
 522            grad_shield_side(j,kk,i)=0.0d0
 523            grad_shield_loc(j,kk,i)=0.0d0
 524
 525 C grad_shield_side_ca is Calfa sidechain gradient
 526
 527
 528 C           grad_shield_side_ca(j,kk,i)=0.0d0
 529           enddo
 530           do intertyp=1,3
 531            gloc_sc(intertyp,i,icg)=0.0d0
 532           enddo
 533         enddo
 534       enddo
 535 #ifndef DFA
 536       do i=1,nres
 537         do j=1,3
 538           gdfad(j,i)=0.0d0
 539           gdfat(j,i)=0.0d0
 540           gdfan(j,i)=0.0d0
 541           gdfab(j,i)=0.0d0
 542         enddo
 543       enddo
 544 #endif
 545 C
 546 C Initialize the gradient of local energy terms.
 547 C
 548       do i=1,4*nres
 549         gloc(i,icg)=0.0D0
 550       enddo
 551       do i=1,nres
 552         gel_loc_loc(i)=0.0d0
 553         gcorr_loc(i)=0.0d0
 554         g_corr5_loc(i)=0.0d0
 555         g_corr6_loc(i)=0.0d0
 556         gel_loc_turn3(i)=0.0d0
 557         gel_loc_turn4(i)=0.0d0
 558         gel_loc_turn6(i)=0.0d0
 559         gsccor_loc(i)=0.0d0
 560       enddo
 561 c initialize gcart and gxcart
 562       do i=0,nres
 563         do j=1,3
 564           gcart(j,i)=0.0d0
 565           gxcart(j,i)=0.0d0
 566         enddo
 567       enddo
 568       return
 569       end
 570 c-------------------------------------------------------------------------
 571       double precision function fdum()
 572       fdum=0.0D0
 573       return
 574       end