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

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