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 #endif
 339       do i=nres,1,-1
 340         do j=1,3
 341           gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 342 !          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 343         enddo
 344 !        write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
 345 !            (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
 346       enddo
 347 ! Correction: dummy residues
 348       if (nnt.gt.1) then
 349         do j=1,3
 350           gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
 351         enddo
 352       endif
 353       if (nct.lt.nres) then
 354         do j=1,3
 355 !          gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
 356           gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
 357         enddo
 358       endif
 359 #ifdef DEBUG
 360       write (iout,*) "CA/SC gradient"
 361       do i=1,nres
 362         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 363      &      (gxcart(j,i),j=1,3)
 364       enddo
 365 #endif
 366       return
 367       end
 368 #endif
 369 C-------------------------------------------------------------------------
 370       subroutine zerograd
 371       implicit none
 372       include 'DIMENSIONS'
 373       include 'COMMON.DERIV'
 374       include 'COMMON.CHAIN'
 375       include 'COMMON.VAR'
 376       include 'COMMON.MD'
 377       include 'COMMON.SCCOR'
 378       include 'COMMON.SHIELD'
 379       integer i,j,kk,intertyp,maxshieldlist
 380       maxshieldlist=0
 381 C
 382 C Initialize Cartesian-coordinate gradient
 383 C
 384       do i=-1,nres
 385         do j=1,3
 386           gvdwx(j,i)=0.0D0
 387           gradx_scp(j,i)=0.0D0
 388           gvdwc(j,i)=0.0D0
 389           gvdwc_scp(j,i)=0.0D0
 390           gvdwc_scpp(j,i)=0.0d0
 391           gelc (j,i)=0.0D0
 392 C below is zero grad for shielding in order: ees (p-p)
 393 C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
 394           gshieldx(j,i)=0.0d0
 395           gshieldc(j,i)=0.0d0
 396           gshieldc_loc(j,i)=0.0d0
 397           gshieldx_ec(j,i)=0.0d0
 398           gshieldc_ec(j,i)=0.0d0
 399           gshieldc_loc_ec(j,i)=0.0d0
 400           gshieldx_t3(j,i)=0.0d0
 401           gshieldc_t3(j,i)=0.0d0
 402           gshieldc_loc_t3(j,i)=0.0d0
 403           gshieldx_t4(j,i)=0.0d0
 404           gshieldc_t4(j,i)=0.0d0
 405           gshieldc_loc_t4(j,i)=0.0d0
 406           gshieldx_ll(j,i)=0.0d0
 407           gshieldc_ll(j,i)=0.0d0
 408           gshieldc_loc_ll(j,i)=0.0d0
 409 C end of zero grad for shielding
 410           gelc_long(j,i)=0.0D0
 411           gradb(j,i)=0.0d0
 412           gradbx(j,i)=0.0d0
 413           gvdwpp(j,i)=0.0d0
 414           gel_loc(j,i)=0.0d0
 415           gel_loc_long(j,i)=0.0d0
 416           ghpbc(j,i)=0.0D0
 417           ghpbx(j,i)=0.0D0
 418           gsaxsc(j,i)=0.0D0
 419           gsaxsx(j,i)=0.0D0
 420           gcorr3_turn(j,i)=0.0d0
 421           gcorr4_turn(j,i)=0.0d0
 422           gradcorr(j,i)=0.0d0
 423           gradcorr_long(j,i)=0.0d0
 424           gradcorr5_long(j,i)=0.0d0
 425           gradcorr6_long(j,i)=0.0d0
 426           gcorr6_turn_long(j,i)=0.0d0
 427           gradcorr5(j,i)=0.0d0
 428           gradcorr6(j,i)=0.0d0
 429           gcorr6_turn(j,i)=0.0d0
 430           gsccorc(j,i)=0.0d0
 431           gsccorx(j,i)=0.0d0
 432           gradc(j,i,icg)=0.0d0
 433           gradx(j,i,icg)=0.0d0
 434           gscloc(j,i)=0.0d0
 435           gsclocx(j,i)=0.0d0
 436           gliptranc(j,i)=0.0d0
 437           gliptranx(j,i)=0.0d0
 438           gradafm(j,i)=0.0d0
 439           grad_shield(j,i)=0.0d0
 440           gg_tube(j,i)=0.0d0
 441           gg_tube_sc(j,i)=0.0d0
 442 C grad_shield_side is Cbeta sidechain gradient
 443           do kk=1,maxshieldlist
 444            grad_shield_side(j,kk,i)=0.0d0
 445            grad_shield_loc(j,kk,i)=0.0d0
 446
 447 C grad_shield_side_ca is Calfa sidechain gradient
 448
 449
 450 C           grad_shield_side_ca(j,kk,i)=0.0d0
 451           enddo
 452           do intertyp=1,3
 453            gloc_sc(intertyp,i,icg)=0.0d0
 454           enddo
 455         enddo
 456       enddo
 457 #ifndef DFA
 458       do i=1,nres
 459         do j=1,3
 460           gdfad(j,i)=0.0d0
 461           gdfat(j,i)=0.0d0
 462           gdfan(j,i)=0.0d0
 463           gdfab(j,i)=0.0d0
 464         enddo
 465       enddo
 466 #endif
 467 C
 468 C Initialize the gradient of local energy terms.
 469 C
 470       do i=1,4*nres
 471         gloc(i,icg)=0.0D0
 472       enddo
 473       do i=1,nres
 474         gel_loc_loc(i)=0.0d0
 475         gcorr_loc(i)=0.0d0
 476         g_corr5_loc(i)=0.0d0
 477         g_corr6_loc(i)=0.0d0
 478         gel_loc_turn3(i)=0.0d0
 479         gel_loc_turn4(i)=0.0d0
 480         gel_loc_turn6(i)=0.0d0
 481         gsccor_loc(i)=0.0d0
 482       enddo
 483 c initialize gcart and gxcart
 484       do i=0,nres
 485         do j=1,3
 486           gcart(j,i)=0.0d0
 487           gxcart(j,i)=0.0d0
 488         enddo
 489       enddo
 490       return
 491       end
 492 c-------------------------------------------------------------------------
 493       double precision function fdum()
 494       fdum=0.0D0
 495       return
 496       end