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

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