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       double precision time00,time01
 214 c
 215 c This subrouting calculates total Cartesian coordinate gradient.
 216 c The subroutine chainbuild_cart and energy MUST be called beforehand.
 217 c
 218 #ifdef TIMING
 219       time00=MPI_Wtime()
 220 #endif
 221       icg=1
 222 #ifdef DEBUG
 223       write (iout,*) "Before sum_gradient"
 224       do i=1,nres-1
 225         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
 226         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
 227       enddo
 228       write (iout,*) "gsaxsc, gsaxcx"
 229       do i=1,nres-1
 230         write (iout,*) i," gsaxsc  ",(gsaxsc(j,i),j=1,3)
 231         write (iout,*) i," gsaxsx  ",(gsaxsx(j,i),j=1,3)
 232       enddo
 233 #endif
 234       call sum_gradient
 235 #ifdef TIMING
 236 #endif
 237 #ifdef DEBUG
 238       write (iout,*) "After sum_gradient"
 239       do i=1,nres-1
 240         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
 241         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
 242       enddo
 243 #endif
 244 c If performing constraint dynamics, add the gradients of the constraint energy
 245       if(usampl.and.totT.gt.eq_time) then
 246 #ifdef DEBUG
 247          write (iout,*) "dudconst, duscdiff, dugamma,dutheta"
 248          write (iout,*) "wumb",wumb
 249          do i=1,nct
 250            write (iout,'(i5,3f10.5,5x,3f10.5,5x,2f10.5)')
 251      &      i,(dudconst(j,i),j=1,3),(duscdiff(j,i),j=1,3),
 252      &      dugamma(i),dutheta(i)
 253          enddo
 254 #endif
 255          do i=1,nct
 256            do j=1,3
 257              gradc(j,i,icg)=gradc(j,i,icg)+
 258      &          wumb*(dudconst(j,i)+duscdiff(j,i))
 259              gradx(j,i,icg)=gradx(j,i,icg)+
 260      &          wumb*(dudxconst(j,i)+duscdiffx(j,i))
 261            enddo
 262          enddo
 263          do i=1,nres-3
 264            gloc(i,icg)=gloc(i,icg)+wumb*dugamma(i)
 265          enddo
 266          do i=1,nres-2
 267            gloc(nphi+i,icg)=gloc(nphi+i,icg)+wumb*dutheta(i)
 268          enddo
 269       endif
 270 #ifdef TIMING
 271       time01=MPI_Wtime()
 272 #endif
 273       call intcartderiv
 274 #ifdef TIMING
 275       time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
 276 #endif
 277 cd      call checkintcartgrad
 278 cd      write(iout,*) 'calling int_to_cart'
 279 #ifdef DEBUG
 280       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 281 #endif
 282       do i=0,nct
 283         do j=1,3
 284           gcart(j,i)=gradc(j,i,icg)
 285           gxcart(j,i)=gradx(j,i,icg)
 286         enddo
 287 #ifdef DEBUG
 288         if (i.eq.0) then
 289         write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
 290      &    (gxcart(j,i),j=1,3)
 291         else if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
 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      &    gloc(ialph(i,1),icg),gloc(ialph(i,1)+nside,icg)
 295         else
 296         write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
 297      &    (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg)
 298         endif
 299         call flush(iout)
 300 #endif
 301       enddo
 302 #ifdef TIMING
 303       time01=MPI_Wtime()
 304 #endif
 305       call int_to_cart
 306 #ifdef TIMING
 307       time_inttocart=time_inttocart+MPI_Wtime()-time01
 308 #endif
 309 #ifdef DEBUG
 310       write (iout,*) "gcart and gxcart after int_to_cart"
 311       do i=0,nres-1
 312         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 313      &      (gxcart(j,i),j=1,3)
 314       enddo
 315 #endif
 316 #ifdef TIMING
 317       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 318 #endif
 319       return
 320       end
 321 c---------------------------------------------------------------------------
 322 #ifdef FIVEDIAG
 323       subroutine grad_transform
 324       implicit none
 325       include 'DIMENSIONS'
 326 #ifdef MPI
 327       include 'mpif.h'
 328 #endif
 329       include 'COMMON.CONTROL'
 330       include 'COMMON.CHAIN'
 331       include 'COMMON.DERIV'
 332       include 'COMMON.VAR'
 333       include 'COMMON.INTERACT'
 334       include 'COMMON.FFIELD'
 335       include 'COMMON.MD'
 336       include 'COMMON.QRESTR'
 337       include 'COMMON.IOUNITS'
 338       include 'COMMON.TIME1'
 339       integer i,j,kk
 340 #ifdef DEBUG
 341       write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
 342       write (iout,*) "dC/dX gradient"
 343       do i=0,nres
 344         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 345      &      (gxcart(j,i),j=1,3)
 346       enddo
 347 #endif
 348       do i=nres,1,-1
 349         do j=1,3
 350           gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 351 !          gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
 352         enddo
 353 !        write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
 354 !            (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
 355       enddo
 356 ! Correction: dummy residues
 357       do i=2,nres
 358         if (itype(i-1).eq.ntyp1 .and. itype(i).ne.ntyp1) then
 359           gcart(:,i)=gcart(:,i)+gcart(:,i-1)
 360         else if (itype(i-1).ne.ntyp1 .and. itype(i).eq.ntyp1) then
 361           gcart(:,i-1)=gcart(:,i-1)+gcart(:,i)
 362         endif
 363       enddo
 364 c      if (nnt.gt.1) then
 365 c        do j=1,3
 366 c          gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
 367 c        enddo
 368 c      endif
 369 c      if (nct.lt.nres) then
 370 c        do j=1,3
 371 c!          gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
 372 c          gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
 373 c        enddo
 374 c      endif
 375 #ifdef DEBUG
 376       write (iout,*) "CA/SC gradient"
 377       do i=1,nres
 378         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 379      &      (gxcart(j,i),j=1,3)
 380       enddo
 381 #endif
 382       return
 383       end
 384 #endif
 385 C-------------------------------------------------------------------------
 386       subroutine zerograd
 387       implicit none
 388       include 'DIMENSIONS'
 389       include 'COMMON.DERIV'
 390       include 'COMMON.CHAIN'
 391       include 'COMMON.VAR'
 392       include 'COMMON.MD'
 393       include 'COMMON.SCCOR'
 394       include 'COMMON.SHIELD'
 395       integer i,j,kk,intertyp,maxshieldlist
 396       maxshieldlist=0
 397 C
 398 C Initialize Cartesian-coordinate gradient
 399 C
 400       do i=-1,nres
 401         do j=1,3
 402           gvdwx(j,i)=0.0D0
 403           gradx_scp(j,i)=0.0D0
 404           gvdwc(j,i)=0.0D0
 405           gvdwc_scp(j,i)=0.0D0
 406           gvdwc_scpp(j,i)=0.0d0
 407           gelc (j,i)=0.0D0
 408 C below is zero grad for shielding in order: ees (p-p)
 409 C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
 410           gshieldx(j,i)=0.0d0
 411           gshieldc(j,i)=0.0d0
 412           gshieldc_loc(j,i)=0.0d0
 413           gshieldx_ec(j,i)=0.0d0
 414           gshieldc_ec(j,i)=0.0d0
 415           gshieldc_loc_ec(j,i)=0.0d0
 416           gshieldx_t3(j,i)=0.0d0
 417           gshieldc_t3(j,i)=0.0d0
 418           gshieldc_loc_t3(j,i)=0.0d0
 419           gshieldx_t4(j,i)=0.0d0
 420           gshieldc_t4(j,i)=0.0d0
 421           gshieldc_loc_t4(j,i)=0.0d0
 422           gshieldx_ll(j,i)=0.0d0
 423           gshieldc_ll(j,i)=0.0d0
 424           gshieldc_loc_ll(j,i)=0.0d0
 425 C end of zero grad for shielding
 426           gelc_long(j,i)=0.0D0
 427           gradb(j,i)=0.0d0
 428           gradbx(j,i)=0.0d0
 429           gvdwpp(j,i)=0.0d0
 430           gel_loc(j,i)=0.0d0
 431           gel_loc_long(j,i)=0.0d0
 432           ghpbc(j,i)=0.0D0
 433           ghpbx(j,i)=0.0D0
 434           gsaxsc(j,i)=0.0D0
 435           gsaxsx(j,i)=0.0D0
 436           gcorr3_turn(j,i)=0.0d0
 437           gcorr4_turn(j,i)=0.0d0
 438           gradcorr(j,i)=0.0d0
 439           gradcorr_long(j,i)=0.0d0
 440           gradcorr5_long(j,i)=0.0d0
 441           gradcorr6_long(j,i)=0.0d0
 442           gcorr6_turn_long(j,i)=0.0d0
 443           gradcorr5(j,i)=0.0d0
 444           gradcorr6(j,i)=0.0d0
 445           gcorr6_turn(j,i)=0.0d0
 446           gsccorc(j,i)=0.0d0
 447           gsccorx(j,i)=0.0d0
 448           gradc(j,i,icg)=0.0d0
 449           gradx(j,i,icg)=0.0d0
 450           gscloc(j,i)=0.0d0
 451           gsclocx(j,i)=0.0d0
 452           gliptranc(j,i)=0.0d0
 453           gliptranx(j,i)=0.0d0
 454           gradafm(j,i)=0.0d0
 455           grad_shield(j,i)=0.0d0
 456           gg_tube(j,i)=0.0d0
 457           gg_tube_sc(j,i)=0.0d0
 458 C grad_shield_side is Cbeta sidechain gradient
 459           do kk=1,maxshieldlist
 460            grad_shield_side(j,kk,i)=0.0d0
 461            grad_shield_loc(j,kk,i)=0.0d0
 462
 463 C grad_shield_side_ca is Calfa sidechain gradient
 464
 465
 466 C           grad_shield_side_ca(j,kk,i)=0.0d0
 467           enddo
 468           do intertyp=1,3
 469            gloc_sc(intertyp,i,icg)=0.0d0
 470           enddo
 471         enddo
 472       enddo
 473 #ifndef DFA
 474       do i=1,nres
 475         do j=1,3
 476           gdfad(j,i)=0.0d0
 477           gdfat(j,i)=0.0d0
 478           gdfan(j,i)=0.0d0
 479           gdfab(j,i)=0.0d0
 480         enddo
 481       enddo
 482 #endif
 483 C
 484 C Initialize the gradient of local energy terms.
 485 C
 486       do i=1,4*nres
 487         gloc(i,icg)=0.0D0
 488       enddo
 489       do i=1,nres
 490         gel_loc_loc(i)=0.0d0
 491         gcorr_loc(i)=0.0d0
 492         g_corr5_loc(i)=0.0d0
 493         g_corr6_loc(i)=0.0d0
 494         gel_loc_turn3(i)=0.0d0
 495         gel_loc_turn4(i)=0.0d0
 496         gel_loc_turn6(i)=0.0d0
 497         gsccor_loc(i)=0.0d0
 498       enddo
 499 c initialize gcart and gxcart
 500       do i=0,nres
 501         do j=1,3
 502           gcart(j,i)=0.0d0
 503           gxcart(j,i)=0.0d0
 504         enddo
 505       enddo
 506       return
 507       end
 508 c-------------------------------------------------------------------------
 509       double precision function fdum()
 510       fdum=0.0D0
 511       return
 512       end