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

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