--- /dev/null
+ 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'
+ include 'COMMON.SCCOR'
+ 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'
+ include 'COMMON.SCCOR'
+c
+c This subrouting calculates total Cartesian coordinate gradient.
+c The subroutine chainbuild_cart and energy MUST be called beforehand.
+c
+c do i=1,nres
+c write (iout,*) "przed sum_grad", gloc_sc(1,i,icg),gloc(i,icg)
+c enddo
+
+#ifdef TIMING
+ time00=MPI_Wtime()
+#endif
+ icg=1
+ call sum_gradient
+#ifdef TIMING
+#endif
+c do i=1,nres
+c write (iout,*) "checkgrad", gloc_sc(1,i,icg),gloc(i,icg)
+c enddo
+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
+c if(usampl.and.totT.gt.eq_time) then
+c do i=1,nct
+c do j=1,3
+c gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i)
+c gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i)
+c enddo
+c enddo
+c do i=1,nres-3
+c gloc(i,icg)=gloc(i,icg)+dugamma(i)
+c enddo
+c do i=1,nres-2
+c gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+c enddo
+c 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'
+ include 'COMMON.SCCOR'
+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
+ do intertyp=1,3
+ gloc_sc(intertyp,i,icg)=0.0d0
+ enddo
+ 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