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