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