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 c
 256 c     integer iistart,iiend
 257 c
 258       include 'DIMENSIONS'
 259       include 'COMMON.LOCAL'
 260 #ifdef MPI
 261       include 'mpif.h'
 262 #endif
 263       include 'COMMON.CONTROL'
 264       include 'COMMON.CHAIN'
 265       include 'COMMON.DERIV'
 266       include 'COMMON.VAR'
 267       include 'COMMON.INTERACT'
 268       include 'COMMON.FFIELD'
 269       include 'COMMON.MD'
 270       include 'COMMON.IOUNITS'
 271       include 'COMMON.TIME1'
 272       include 'COMMON.SCCOR'
 273 c
 274 c This subrouting calculates total Cartesian coordinate gradient.
 275 c The subroutine chainbuild_cart and energy MUST be called beforehand.
 276 c
 277 c        do i=1,nres
 278 c        write (iout,*) "przed sum_grad", gloc_sc(1,i,icg),gloc(i,icg)
 279 c        enddo
 280
 281 #ifdef TIMING
 282       time00=MPI_Wtime()
 283 #endif
 284       icg=1
 285 #ifdef DEBUG
 286 c     write (iout,*) "in cartgrad before sum_: duscdiff and duscdiffx"
 287       write (iout,*) "------ Before call sum_ in cargrad ------"
 288       do i=1,nres
 289          write (iout,*) i,(duscdiff(j,i),j=1,3)
 290          write (iout,*) i,(duscdiffx(j,i),j=1,3)
 291 c        write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg)
 292 c        write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
 293 c        write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
 294       enddo
 295 #endif
 296       call sum_gradient
 297 c
 298 #ifdef TIMING
 299 #endif
 300 c        do i=1,nres
 301 c        write (iout,*) "checkgrad", gloc_sc(1,i,icg),gloc(i,icg)
 302 c        enddo
 303 #ifdef DEBUG
 304         write (iout,*) "------ After sum_gradient in cartgrad ------"
 305 c       do i=1,nres-1
 306         do i=1,nres
 307          write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg)
 308          write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
 309          write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
 310         enddo
 311 #endif
 312 c If performing constraint dynamics, add the gradients of the constraint energy
 313 #ifdef DEBUG
 314       write (iout,*) "in cartgrad: dutheta, duscdiff and duscdiffx"
 315       do i=1,nres
 316         write (iout,*) i,dutheta(i)
 317         write (iout,*) i,(duscdiff(j,i),j=1,3)
 318         write (iout,*) i,(duscdiffx(j,i),j=1,3)
 319       enddo
 320 #endif
 321       if(usampl.and.totT.gt.eq_time) then
 322 c      if(usampl.and.totT.gt.eq_time .or. constr_homology.gt.0) then
 323 c
 324 c     Setting suited bounds for HM restrs
 325 c
 326          do i=1,nct
 327            do j=1,3
 328              gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)
 329              gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)
 330            enddo
 331 #ifdef DEBUG
 332          write (iout,*) "i",i," gradc",(gradc(j,i,icg),j=1,3)
 333          write (iout,*) "i",i," gradx",(gradx(j,i,icg),j=1,3)
 334 #endif
 335          enddo
 336          do i=1,nres-3
 337            gloc(i,icg)=gloc(i,icg)+dugamma(i)
 338          enddo
 339 c
 340          do i=1,nres-2
 341            gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
 342 #ifdef DEBUG
 343          write (iout,*) "nphi+i",nphi+i," gloc",gloc(nphi+i,icg)
 344 #endif
 345          enddo
 346       endif
 347 #ifdef TIMING
 348       time01=MPI_Wtime()
 349 #endif
 350       call intcartderiv
 351 #ifdef TIMING
 352       time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
 353 #endif
 354 cd      call checkintcartgrad
 355 cd      write(iout,*) 'calling int_to_cart'
 356 cd      write (iout,*) "gcart, gxcart, gloc before int_to_cart"
 357       do i=1,nct
 358         do j=1,3
 359           gcart(j,i)=gradc(j,i,icg)
 360           gxcart(j,i)=gradx(j,i,icg)
 361         enddo
 362 cd        write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),
 363 cd     &    (gxcart(j,i),j=1,3),gloc(i,icg)
 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 cd      write (iout,*) "gcart and gxcart after int_to_cart"
 373 cd      do i=0,nres-1
 374 cd        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
 375 cd     &      (gxcart(j,i),j=1,3)
 376 cd      enddo
 377 #ifdef TIMING
 378       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
 379 #endif
 380       return
 381       end
 382 C-------------------------------------------------------------------------
 383       subroutine zerograd
 384       implicit real*8 (a-h,o-z)
 385       include 'DIMENSIONS'
 386       include 'COMMON.DERIV'
 387       include 'COMMON.CHAIN'
 388       include 'COMMON.VAR'
 389       include 'COMMON.MD'
 390       include 'COMMON.SCCOR'
 391 C
 392 C Initialize Cartesian-coordinate gradient
 393 C
 394       do i=1,nres
 395         do j=1,3
 396           gvdwx(j,i)=0.0D0
 397           gvdwxT(j,i)=0.0D0
 398           gradx_scp(j,i)=0.0D0
 399           gvdwc(j,i)=0.0D0
 400           gvdwcT(j,i)=0.0D0
 401           gvdwc_scp(j,i)=0.0D0
 402           gvdwc_scpp(j,i)=0.0d0
 403           gelc (j,i)=0.0D0
 404           gelc_long(j,i)=0.0D0
 405           gradb(j,i)=0.0d0
 406           gradbx(j,i)=0.0d0
 407           gvdwpp(j,i)=0.0d0
 408           gel_loc(j,i)=0.0d0
 409           gel_loc_long(j,i)=0.0d0
 410           ghpbc(j,i)=0.0D0
 411           ghpbx(j,i)=0.0D0
 412           gcorr3_turn(j,i)=0.0d0
 413           gcorr4_turn(j,i)=0.0d0
 414           gradcorr(j,i)=0.0d0
 415           gradcorr_long(j,i)=0.0d0
 416           gradcorr5_long(j,i)=0.0d0
 417           gradcorr6_long(j,i)=0.0d0
 418           gcorr6_turn_long(j,i)=0.0d0
 419           gradcorr5(j,i)=0.0d0
 420           gradcorr6(j,i)=0.0d0
 421           gcorr6_turn(j,i)=0.0d0
 422           gsccorc(j,i)=0.0d0
 423           gsccorx(j,i)=0.0d0
 424           gradc(j,i,icg)=0.0d0
 425           gradx(j,i,icg)=0.0d0
 426           gscloc(j,i)=0.0d0
 427           gsclocx(j,i)=0.0d0
 428           do intertyp=1,3
 429            gloc_sc(intertyp,i,icg)=0.0d0
 430           enddo
 431 #ifndef DFA
 432           gdfad(j,i)=0.0d0
 433           gdfat(j,i)=0.0d0
 434           gdfan(j,i)=0.0d0
 435           gdfab(j,i)=0.0d0
 436 #endif
 437         enddo
 438       enddo
 439 c
 440 c Initialize the gradients of local restraints
 441 c
 442       do i=1,nres
 443         dutheta(i)=0.0d0
 444         dugamma(i)=0.0d0
 445         do j=1,3
 446           duscdiff(j,i)=0.0d0
 447           duscdiffx(j,i)=0.0d0
 448         enddo
 449       enddo
 450 C
 451 C Initialize the gradient of local energy terms.
 452 C
 453       do i=1,4*nres
 454         gloc(i,icg)=0.0D0
 455       enddo
 456       do i=1,nres
 457         gel_loc_loc(i)=0.0d0
 458         gcorr_loc(i)=0.0d0
 459         g_corr5_loc(i)=0.0d0
 460         g_corr6_loc(i)=0.0d0
 461         gel_loc_turn3(i)=0.0d0
 462         gel_loc_turn4(i)=0.0d0
 463         gel_loc_turn6(i)=0.0d0
 464         gsccor_loc(i)=0.0d0
 465       enddo
 466 c initialize gcart and gxcart
 467       do i=0,nres
 468         do j=1,3
 469           gcart(j,i)=0.0d0
 470           gxcart(j,i)=0.0d0
 471         enddo
 472       enddo
 473       return
 474       end
 475 c-------------------------------------------------------------------------
 476       double precision function fdum()
 477       fdum=0.0D0
 478       return
 479       end