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