X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_Eshel%2Fgradient_p.F;fp=source%2Funres%2Fsrc_Eshel%2Fgradient_p.F;h=9a2549709805450baf542c63695821239eabf687;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_Eshel/gradient_p.F b/source/unres/src_Eshel/gradient_p.F new file mode 100644 index 0000000..9a25497 --- /dev/null +++ b/source/unres/src_Eshel/gradient_p.F @@ -0,0 +1,421 @@ + 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