subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.MD' include 'COMMON.IOUNITS' external ufparm integer uiparm(1) double precision urparm(1) dimension x(maxvar),g(maxvar) c c This subroutine calculates total internal coordinate gradient. c Depending on the number of function evaluations, either whole energy c is evaluated beforehand, Cartesian coordinates and their derivatives in c internal coordinates are reevaluated or only the cartesian-in-internal c coordinate derivatives are evaluated. The subroutine was designed to work c with SUMSL. c c icg=mod(nf,2)+1 cd print *,'grad',nf,icg if (nf-nfl+1) 20,30,40 20 call func(n,x,nf,f,uiparm,urparm,ufparm) c write (iout,*) 'grad 20' if (nf.eq.0) return goto 40 30 call var_to_geom(n,x) call chainbuild c write (iout,*) 'grad 30' C C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. C 40 call cartder c write (iout,*) 'grad 40' c print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon C C Convert the Cartesian gradient into internal-coordinate gradient. C ind=0 ind1=0 do i=1,nres-2 gthetai=0.0D0 gphii=0.0D0 do j=i+1,nres-1 ind=ind+1 c ind=indmat(i,j) c print *,'GRAD: i=',i,' jc=',j,' ind=',ind do k=1,3 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) enddo do k=1,3 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) enddo enddo do j=i+1,nres-1 ind1=ind1+1 c ind1=indmat(i,j) c print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1 do k=1,3 gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg) gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg) enddo enddo if (i.gt.1) g(i-1)=gphii if (n.gt.nphi) g(nphi+i)=gthetai enddo if (n.le.nphi+ntheta) goto 10 do i=2,nres-1 if (itype(i).ne.10) then galphai=0.0D0 gomegai=0.0D0 do k=1,3 galphai=galphai+dxds(k,i)*gradx(k,i,icg) enddo do k=1,3 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) enddo g(ialph(i,1))=galphai g(ialph(i,1)+nside)=gomegai endif enddo C C Add the components corresponding to local energy terms. C 10 continue do i=1,nvar cd write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg) g(i)=g(i)+gloc(i,icg) enddo C Uncomment following three lines for diagnostics. cd call intout cd call briefout(0,0.0d0) cd write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n) return end C------------------------------------------------------------------------- subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' external ufparm integer uiparm(1) double precision urparm(1) dimension x(maxvar),g(maxvar) icg=mod(nf,2)+1 if (nf-nfl+1) 20,30,40 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm) c write (iout,*) 'grad 20' if (nf.eq.0) return goto 40 30 continue #ifdef OSF c Intercept NaNs in the coordinates c write(iout,*) (var(i),i=1,nvar) x_sum=0.D0 do i=1,n x_sum=x_sum+x(i) enddo if (x_sum.ne.x_sum) then write(iout,*)" *** grad_restr : Found NaN in coordinates" call flush(iout) print *," *** grad_restr : Found NaN in coordinates" return endif #endif call var_to_geom_restr(n,x) call chainbuild C C Evaluate the derivatives of virtual bond lengths and SC vectors in variables. C 40 call cartder C C Convert the Cartesian gradient into internal-coordinate gradient. C ig=0 ind=nres-2 do i=2,nres-2 IF (mask_phi(i+2).eq.1) THEN gphii=0.0D0 do j=i+1,nres-1 ind=ind+1 do k=1,3 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg) gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg) enddo enddo ig=ig+1 g(ig)=gphii ELSE ind=ind+nres-1-i ENDIF enddo ind=0 do i=1,nres-2 IF (mask_theta(i+2).eq.1) THEN ig=ig+1 gthetai=0.0D0 do j=i+1,nres-1 ind=ind+1 do k=1,3 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg) gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg) enddo enddo g(ig)=gthetai ELSE ind=ind+nres-1-i ENDIF enddo do i=2,nres-1 if (itype(i).ne.10) then IF (mask_side(i).eq.1) THEN ig=ig+1 galphai=0.0D0 do k=1,3 galphai=galphai+dxds(k,i)*gradx(k,i,icg) enddo g(ig)=galphai ENDIF endif enddo do i=2,nres-1 if (itype(i).ne.10) then IF (mask_side(i).eq.1) THEN ig=ig+1 gomegai=0.0D0 do k=1,3 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg) enddo g(ig)=gomegai ENDIF endif enddo C C Add the components corresponding to local energy terms. C ig=0 igall=0 do i=4,nres igall=igall+1 if (mask_phi(i).eq.1) then ig=ig+1 g(ig)=g(ig)+gloc(igall,icg) endif enddo do i=3,nres igall=igall+1 if (mask_theta(i).eq.1) then ig=ig+1 g(ig)=g(ig)+gloc(igall,icg) endif enddo do ij=1,2 do i=2,nres-1 if (itype(i).ne.10) then igall=igall+1 if (mask_side(i).eq.1) then ig=ig+1 g(ig)=g(ig)+gloc(igall,icg) endif endif enddo enddo cd do i=1,ig cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i) cd enddo return end C------------------------------------------------------------------------- subroutine cartgrad implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.MD' include 'COMMON.IOUNITS' include 'COMMON.TIME1' c c This subrouting calculates total Cartesian coordinate gradient. c The subroutine chainbuild_cart and energy MUST be called beforehand. c #ifdef TIMING time00=MPI_Wtime() #endif icg=1 call sum_gradient #ifdef TIMING #endif cd write (iout,*) "After sum_gradient" cd do i=1,nres-1 cd write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) cd write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) cd enddo c If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then do i=1,nct do j=1,3 gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i) gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i) enddo enddo do i=1,nres-3 gloc(i,icg)=gloc(i,icg)+dugamma(i) enddo do i=1,nres-2 gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i) enddo endif #ifdef TIMING time01=MPI_Wtime() #endif call intcartderiv #ifdef TIMING time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01 #endif cd call checkintcartgrad cd write(iout,*) 'calling int_to_cart' cd write (iout,*) "gcart, gxcart, gloc before int_to_cart" do i=1,nct do j=1,3 gcart(j,i)=gradc(j,i,icg) gxcart(j,i)=gradx(j,i,icg) enddo cd write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3), cd & (gxcart(j,i),j=1,3),gloc(i,icg) enddo #ifdef TIMING time01=MPI_Wtime() #endif call int_to_cart #ifdef TIMING time_inttocart=time_inttocart+MPI_Wtime()-time01 #endif cd write (iout,*) "gcart and gxcart after int_to_cart" cd do i=0,nres-1 cd write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), cd & (gxcart(j,i),j=1,3) cd enddo #ifdef TIMING time_cartgrad=time_cartgrad+MPI_Wtime()-time00 #endif return end C------------------------------------------------------------------------- subroutine zerograd implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.DERIV' include 'COMMON.CHAIN' include 'COMMON.VAR' include 'COMMON.MD' C C Initialize Cartesian-coordinate gradient C do i=1,nres do j=1,3 gvdwx(j,i)=0.0D0 gvdwxT(j,i)=0.0D0 gradx_scp(j,i)=0.0D0 gvdwc(j,i)=0.0D0 gvdwcT(j,i)=0.0D0 gvdwc_scp(j,i)=0.0D0 gvdwc_scpp(j,i)=0.0d0 gelc (j,i)=0.0D0 gelc_long(j,i)=0.0D0 gradb(j,i)=0.0d0 gradbx(j,i)=0.0d0 gvdwpp(j,i)=0.0d0 gel_loc(j,i)=0.0d0 gel_loc_long(j,i)=0.0d0 ghpbc(j,i)=0.0D0 ghpbx(j,i)=0.0D0 gcorr3_turn(j,i)=0.0d0 gcorr4_turn(j,i)=0.0d0 gradcorr(j,i)=0.0d0 gradcorr_long(j,i)=0.0d0 gradcorr5_long(j,i)=0.0d0 gradcorr6_long(j,i)=0.0d0 gcorr6_turn_long(j,i)=0.0d0 gradcorr5(j,i)=0.0d0 gradcorr6(j,i)=0.0d0 gcorr6_turn(j,i)=0.0d0 gsccorc(j,i)=0.0d0 gsccorx(j,i)=0.0d0 gradc(j,i,icg)=0.0d0 gradx(j,i,icg)=0.0d0 gscloc(j,i)=0.0d0 gsclocx(j,i)=0.0d0 enddo enddo C C Initialize the gradient of local energy terms. C do i=1,4*nres gloc(i,icg)=0.0D0 enddo do i=1,nres gel_loc_loc(i)=0.0d0 gcorr_loc(i)=0.0d0 g_corr5_loc(i)=0.0d0 g_corr6_loc(i)=0.0d0 gel_loc_turn3(i)=0.0d0 gel_loc_turn4(i)=0.0d0 gel_loc_turn6(i)=0.0d0 gsccor_loc(i)=0.0d0 enddo c initialize gcart and gxcart do i=0,nres do j=1,3 gcart(j,i)=0.0d0 gxcart(j,i)=0.0d0 enddo enddo return end c------------------------------------------------------------------------- double precision function fdum() fdum=0.0D0 return end