subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm) implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.MD' include 'COMMON.QRESTR' include 'COMMON.IOUNITS' integer n,nf double precision ufparm external ufparm integer uiparm(1) double precision urparm(1) double precision x(n),g(n) integer i,j,k,ind,ind1 double precision f,gthetai,gphii,galphai,gomegai 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 c Add the usampl contributions if (usampl) then 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 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 none include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.IOUNITS' integer n,nf double precision ufparm external ufparm integer uiparm(1) double precision urparm(1) double precision x(maxvar),g(maxvar) integer i,j,k,ig,ind,ij,igall double precision f,gthetai,gphii,galphai,gomegai 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 none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.MD' include 'COMMON.QRESTR' include 'COMMON.IOUNITS' include 'COMMON.TIME1' integer i,j,kk 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 #ifdef DEBUG write (iout,*) "Before sum_gradient" do i=1,nres-1 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) enddo write (iout,*) "gsaxsc, gsaxcx" do i=1,nres-1 write (iout,*) i," gsaxsc ",(gsaxsc(j,i),j=1,3) write (iout,*) i," gsaxsx ",(gsaxsx(j,i),j=1,3) enddo #endif call sum_gradient #ifdef TIMING #endif #ifdef DEBUG write (iout,*) "After sum_gradient" do i=1,nres-1 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3) write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3) enddo #endif c If performing constraint dynamics, add the gradients of the constraint energy if(usampl.and.totT.gt.eq_time) then #ifdef DEBUG write (iout,*) "dudconst, duscdiff, dugamma,dutheta" write (iout,*) "wumb",wumb do i=1,nct write (iout,'(i5,3f10.5,5x,3f10.5,5x,2f10.5)') & i,(dudconst(j,i),j=1,3),(duscdiff(j,i),j=1,3), & dugamma(i),dutheta(i) enddo #endif do i=1,nct do j=1,3 gradc(j,i,icg)=gradc(j,i,icg)+ & wumb*(dudconst(j,i)+duscdiff(j,i)) gradx(j,i,icg)=gradx(j,i,icg)+ & wumb*(dudxconst(j,i)+duscdiffx(j,i)) enddo enddo do i=1,nres-3 gloc(i,icg)=gloc(i,icg)+wumb*dugamma(i) enddo do i=1,nres-2 gloc(nphi+i,icg)=gloc(nphi+i,icg)+wumb*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' #ifdef DEBUG write (iout,*) "gcart, gxcart, gloc before int_to_cart" #endif do i=1,nct do j=1,3 gcart(j,i)=gradc(j,i,icg) gxcart(j,i)=gradx(j,i,icg) enddo #ifdef DEBUG if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg), & gloc(ialph(i,1),icg),gloc(ialph(i,1)+nside,icg) else write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg) endif call flush(iout) #endif enddo #ifdef TIMING time01=MPI_Wtime() #endif call int_to_cart #ifdef TIMING time_inttocart=time_inttocart+MPI_Wtime()-time01 #endif #ifdef DEBUG write (iout,*) "gcart and gxcart after int_to_cart" do i=0,nres-1 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo #endif #ifdef TIMING time_cartgrad=time_cartgrad+MPI_Wtime()-time00 #endif return end c--------------------------------------------------------------------------- #ifdef FIVEDIAG subroutine grad_transform implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' #endif include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.FFIELD' include 'COMMON.MD' include 'COMMON.QRESTR' include 'COMMON.IOUNITS' include 'COMMON.TIME1' integer i,j,kk #ifdef DEBUG write (iout,*)"Converting virtual-bond gradient to CA/SC gradient" #endif do i=nres,1,-1 do j=1,3 gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i) ! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i) enddo ! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), & ! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3) enddo ! Correction: dummy residues if (nnt.gt.1) then do j=1,3 gcart(j,nnt)=gcart(j,nnt)+gcart(j,1) enddo endif if (nct.lt.nres) then do j=1,3 ! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres) gcart(j,nct)=gcart(j,nct)+gcart(j,nres) enddo endif #ifdef DEBUG write (iout,*) "CA/SC gradient" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo #endif return end #endif C------------------------------------------------------------------------- subroutine zerograd implicit none include 'DIMENSIONS' include 'COMMON.DERIV' include 'COMMON.CHAIN' include 'COMMON.VAR' include 'COMMON.MD' include 'COMMON.SCCOR' include 'COMMON.SHIELD' integer i,j,kk,intertyp,maxshieldlist maxshieldlist=0 C C Initialize Cartesian-coordinate gradient C do i=-1,nres do j=1,3 gvdwx(j,i)=0.0D0 gradx_scp(j,i)=0.0D0 gvdwc(j,i)=0.0D0 gvdwc_scp(j,i)=0.0D0 gvdwc_scpp(j,i)=0.0d0 gelc (j,i)=0.0D0 C below is zero grad for shielding in order: ees (p-p) C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain gshieldx(j,i)=0.0d0 gshieldc(j,i)=0.0d0 gshieldc_loc(j,i)=0.0d0 gshieldx_ec(j,i)=0.0d0 gshieldc_ec(j,i)=0.0d0 gshieldc_loc_ec(j,i)=0.0d0 gshieldx_t3(j,i)=0.0d0 gshieldc_t3(j,i)=0.0d0 gshieldc_loc_t3(j,i)=0.0d0 gshieldx_t4(j,i)=0.0d0 gshieldc_t4(j,i)=0.0d0 gshieldc_loc_t4(j,i)=0.0d0 gshieldx_ll(j,i)=0.0d0 gshieldc_ll(j,i)=0.0d0 gshieldc_loc_ll(j,i)=0.0d0 C end of zero grad for shielding 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 gsaxsc(j,i)=0.0D0 gsaxsx(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 gliptranc(j,i)=0.0d0 gliptranx(j,i)=0.0d0 gradafm(j,i)=0.0d0 grad_shield(j,i)=0.0d0 gg_tube(j,i)=0.0d0 gg_tube_sc(j,i)=0.0d0 C grad_shield_side is Cbeta sidechain gradient do kk=1,maxshieldlist grad_shield_side(j,kk,i)=0.0d0 grad_shield_loc(j,kk,i)=0.0d0 C grad_shield_side_ca is Calfa sidechain gradient C grad_shield_side_ca(j,kk,i)=0.0d0 enddo do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo enddo enddo #ifndef DFA do i=1,nres do j=1,3 gdfad(j,i)=0.0d0 gdfat(j,i)=0.0d0 gdfan(j,i)=0.0d0 gdfab(j,i)=0.0d0 enddo enddo #endif 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