source/unres/src-5hdiag-tmp/gradient_p.F

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