source/unres/src_MD-M/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(maxvar),g(maxvar)
  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,*) "gradc, gradx, gloc before sum_gradient"
 285       do i=1,nres-1
 286         write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gradc(j,i,icg),j=1,3),
 287      &    (gradx(j,i,icg),j=1,3),gloc(i,icg)
 288       enddo
 289 #endif
 290       call sum_gradient
 291 #ifdef DEBUG
 292       write (iout,*) "gradc, gradx, gloc after sum_gradient"
 293       do i=1,nres-1
 294         write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gradc(j,i,icg),j=1,3),
 295      &    (gradx(j,i,icg),j=1,3),gloc(i,icg)
 296       enddo
 297 #endif
 298 c If performing constraint dynamics, add the gradients of the constraint energy
 299       if(usampl.and.totT.gt.eq_time) then
 300 #ifdef DEBUG
 301          write (iout,*) "dudconst, duscdiff, dugamma,dutheta"
 302          write (iout,*) "wumb",wumb
 303          do i=1,nct
 304            write (iout,'(i5,3f10.5,5x,3f10.5,5x,2f10.5)')
 305      &      i,(dudconst(j,i),j=1,3),(duscdiff(j,i),j=1,3),
 306      &      dugamma(i),dutheta(i)
 307          enddo
 308 #endif
 309          do i=1,nct
 310            do j=1,3
 311              gradc(j,i,icg)=gradc(j,i,icg)+
 312      &          wumb*(dudconst(j,i)+duscdiff(j,i))
 313              gradx(j,i,icg)=gradx(j,i,icg)+
 314      &          wumb*(dudxconst(j,i)+duscdiffx(j,i))
 315            enddo
 316          enddo
 317          do i=1,nres-3
 318            gloc(i,icg)=gloc(i,icg)+wumb*dugamma(i)
 319          enddo
 320          do i=1,nres-2
 321            gloc(nphi+i,icg)=gloc(nphi+i,icg)+wumb*dutheta(i)
 322          enddo
 323       endif
 324 #ifdef TIMING
 325       time01=MPI_Wtime()
 326 #endif
 327       call intcartderiv
 328 #ifdef TIMING
 329       time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
 330 #endif
 331 cd      call checkintcartgrad
 332 cd      write(iout,*) 'calling int_to_cart'
 333 #ifdef DEBUG
 334       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 335 #endif
 336       do i=1,nct
 337         do j=1,3
 338           gcart(j,i)=gradc(j,i,icg)
 339           gxcart(j,i)=gradx(j,i,icg)
 340         enddo
 341 #ifdef DEBUG
 342         write (iout,'(i5,2(3f10.5,5x),2f10.5)') i,(gcart(j,i),j=1,3),
 343      &    (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg)
 344         call flush(iout)
 345 #endif
 346       enddo
 347 #ifdef TIMING
 348       time01=MPI_Wtime()
 349 #endif
 350       call int_to_cart
 351 #ifdef TIMING
 352       time_inttocart=time_inttocart+MPI_Wtime()-time01
 353 #endif
 354 #ifdef DEBUG
 355       write (iout,*) "gcart and gxcart after int_to_cart"
 356       do i=0,nres-1
 357         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 358      &      (gxcart(j,i),j=1,3)
 359       enddo
 360 #endif
 361 #ifdef TIMING
 362       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 363 #endif
 364       return
 365       end
 366 C-------------------------------------------------------------------------
 367       subroutine zerograd
 368       implicit real*8 (a-h,o-z)
 369       include 'DIMENSIONS'
 370       include 'COMMON.DERIV'
 371       include 'COMMON.CHAIN'
 372       include 'COMMON.VAR'
 373       include 'COMMON.MD'
 374       include 'COMMON.SCCOR'
 375       include 'COMMON.SHIELD'
 376       maxshieldlist=0
 377 C
 378 C Initialize Cartesian-coordinate gradient
 379 C
 380       do i=-1,nres
 381         do j=1,3
 382           gvdwx(j,i)=0.0D0
 383           gradx_scp(j,i)=0.0D0
 384           gvdwc(j,i)=0.0D0
 385           gvdwc_scp(j,i)=0.0D0
 386           gvdwc_scpp(j,i)=0.0d0
 387           gelc (j,i)=0.0D0
 388 C below is zero grad for shielding in order: ees (p-p)
 389 C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
 390           gshieldx(j,i)=0.0d0
 391           gshieldc(j,i)=0.0d0
 392           gshieldc_loc(j,i)=0.0d0
 393           gshieldx_ec(j,i)=0.0d0
 394           gshieldc_ec(j,i)=0.0d0
 395           gshieldc_loc_ec(j,i)=0.0d0
 396           gshieldx_t3(j,i)=0.0d0
 397           gshieldc_t3(j,i)=0.0d0
 398           gshieldc_loc_t3(j,i)=0.0d0
 399           gshieldx_t4(j,i)=0.0d0
 400           gshieldc_t4(j,i)=0.0d0
 401           gshieldc_loc_t4(j,i)=0.0d0
 402           gshieldx_ll(j,i)=0.0d0
 403           gshieldc_ll(j,i)=0.0d0
 404           gshieldc_loc_ll(j,i)=0.0d0
 405 C end of zero grad for shielding
 406           gelc_long(j,i)=0.0D0
 407           gradb(j,i)=0.0d0
 408           gradbx(j,i)=0.0d0
 409           gvdwpp(j,i)=0.0d0
 410           gel_loc(j,i)=0.0d0
 411           gel_loc_long(j,i)=0.0d0
 412           ghpbc(j,i)=0.0D0
 413           ghpbx(j,i)=0.0D0
 414           gcorr3_turn(j,i)=0.0d0
 415           gcorr4_turn(j,i)=0.0d0
 416           gradcorr(j,i)=0.0d0
 417           gradcorr_long(j,i)=0.0d0
 418           gradcorr5_long(j,i)=0.0d0
 419           gradcorr6_long(j,i)=0.0d0
 420           gcorr6_turn_long(j,i)=0.0d0
 421           gradcorr5(j,i)=0.0d0
 422           gradcorr6(j,i)=0.0d0
 423           gcorr6_turn(j,i)=0.0d0
 424           gsccorc(j,i)=0.0d0
 425           gsccorx(j,i)=0.0d0
 426           gradc(j,i,icg)=0.0d0
 427           gradx(j,i,icg)=0.0d0
 428           gscloc(j,i)=0.0d0
 429           gsclocx(j,i)=0.0d0
 430           gliptranc(j,i)=0.0d0
 431           gliptranx(j,i)=0.0d0
 432           gradafm(j,i)=0.0d0
 433           grad_shield(j,i)=0.0d0
 434           gg_tube(j,i)=0.0d0
 435           gg_tube_sc(j,i)=0.0d0
 436 C grad_shield_side is Cbeta sidechain gradient
 437           do kk=1,maxshieldlist
 438            grad_shield_side(j,kk,i)=0.0d0
 439            grad_shield_loc(j,kk,i)=0.0d0
 440
 441 C grad_shield_side_ca is Calfa sidechain gradient
 442
 443
 444 C           grad_shield_side_ca(j,kk,i)=0.0d0
 445           enddo
 446           do intertyp=1,3
 447            gloc_sc(intertyp,i,icg)=0.0d0
 448           enddo
 449         enddo
 450       enddo
 451 C
 452 C Initialize the gradient of local energy terms.
 453 C
 454       do i=1,4*nres
 455         gloc(i,icg)=0.0D0
 456       enddo
 457       do i=1,nres
 458         gel_loc_loc(i)=0.0d0
 459         gcorr_loc(i)=0.0d0
 460         g_corr5_loc(i)=0.0d0
 461         g_corr6_loc(i)=0.0d0
 462         gel_loc_turn3(i)=0.0d0
 463         gel_loc_turn4(i)=0.0d0
 464         gel_loc_turn6(i)=0.0d0
 465         gsccor_loc(i)=0.0d0
 466       enddo
 467 c initialize gcart and gxcart
 468       do i=0,nres
 469         do j=1,3
 470           gcart(j,i)=0.0d0
 471           gxcart(j,i)=0.0d0
 472         enddo
 473       enddo
 474       return
 475       end
 476 c-------------------------------------------------------------------------
 477       double precision function fdum()
 478       fdum=0.0D0
 479       return
 480       end