- subroutine check_cartgrad
-C Check the gradient of Cartesian coordinates in internal coordinates.
- implicit none
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.IOUNITS'
- include 'COMMON.VAR'
- include 'COMMON.CHAIN'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.DERIV'
- double precision temp(6,maxres),xx(3),gg(3),thet,theti,phii,alphi,
- & omegi,aincr2
- integer indmat
- integer i,ii,j,k
- indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
- integer nf
-*
-* Check the gradient of the virtual-bond and SC vectors in the internal
-* coordinates.
-*
- print '("Calling CHECK_ECART",1pd12.3)',aincr
- write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr
- aincr2=0.5d0*aincr
- call chainbuild_extconf
- call cartder
- write (iout,'(a)') '**************** dx/dalpha'
- write (iout,'(a)')
- do i=2,nres-1
- alphi=alph(i)
- alph(i)=alph(i)+aincr
- do k=1,3
- temp(k,i)=dc(k,nres+i)
- enddo
- call chainbuild_extconf
- do k=1,3
- gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
- xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr))
- enddo
- write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)')
- & i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3)
- write (iout,'(a)')
- alph(i)=alphi
- call chainbuild_extconf
- enddo
- write (iout,'(a)')
- write (iout,'(a)') '**************** dx/domega'
- write (iout,'(a)')
- do i=2,nres-1
- omegi=omeg(i)
- omeg(i)=omeg(i)+aincr
- do k=1,3
- temp(k,i)=dc(k,nres+i)
- enddo
- call chainbuild_extconf
- do k=1,3
- gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
- xx(k)=dabs((gg(k)-dxds(k+3,i))/
- & (aincr*dabs(dxds(k+3,i))+aincr))
- enddo
- write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)')
- & i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3)
- write (iout,'(a)')
- omeg(i)=omegi
- call chainbuild_extconf
- enddo
- write (iout,'(a)')
- write (iout,'(a)') '**************** dx/dtheta'
- write (iout,'(a)')
- do i=3,nres
- theti=theta(i)
- theta(i)=theta(i)+aincr
- do j=i-1,nres-1
- do k=1,3
- temp(k,j)=dc(k,nres+j)
- enddo
- enddo
- call chainbuild_extconf
- do j=i-1,nres-1
- ii = indmat(i-2,j)
-c print *,'i=',i-2,' j=',j-1,' ii=',ii
- do k=1,3
- gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
- xx(k)=dabs((gg(k)-dxdv(k,ii))/
- & (aincr*dabs(dxdv(k,ii))+aincr))
- enddo
- write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
- & i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3)
- write(iout,'(a)')
- enddo
- write (iout,'(a)')
- theta(i)=theti
- call chainbuild_extconf
- enddo
- write (iout,'(a)') '***************** dx/dphi'
- write (iout,'(a)')
- do i=4,nres
- phi(i)=phi(i)+aincr
- do j=i-1,nres-1
- do k=1,3
- temp(k,j)=dc(k,nres+j)
- enddo
- enddo
- call chainbuild_extconf
- do j=i-1,nres-1
- ii = indmat(i-2,j)
-c print *,'ii=',ii
- do k=1,3
- gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
- xx(k)=dabs((gg(k)-dxdv(k+3,ii))/
- & (aincr*dabs(dxdv(k+3,ii))+aincr))
- enddo
- write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
- & i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3)
- write(iout,'(a)')
- enddo
- phi(i)=phi(i)-aincr
- call chainbuild_extconf
- enddo
- write (iout,'(a)') '****************** ddc/dtheta'
- do i=1,nres-2
- thet=theta(i+2)
- theta(i+2)=thet+aincr
- do j=i,nres
- do k=1,3
- temp(k,j)=dc(k,j)
- enddo
- enddo
- call chainbuild_extconf
- do j=i+1,nres-1
- ii = indmat(i,j)
-c print *,'ii=',ii
- do k=1,3
- gg(k)=(dc(k,j)-temp(k,j))/aincr
- xx(k)=dabs((gg(k)-dcdv(k,ii))/
- & (aincr*dabs(dcdv(k,ii))+aincr))
- enddo
- write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
- & i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3)
- write (iout,'(a)')
- enddo
- do j=1,nres
- do k=1,3
- dc(k,j)=temp(k,j)
- enddo
- enddo
- theta(i+2)=thet
- enddo
- write (iout,'(a)') '******************* ddc/dphi'
- do i=1,nres-3
- phii=phi(i+3)
- phi(i+3)=phii+aincr
- do j=1,nres
- do k=1,3
- temp(k,j)=dc(k,j)
- enddo
- enddo
- call chainbuild_extconf
- do j=i+2,nres-1
- ii = indmat(i+1,j)
-c print *,'ii=',ii
- do k=1,3
- gg(k)=(dc(k,j)-temp(k,j))/aincr
- xx(k)=dabs((gg(k)-dcdv(k+3,ii))/
- & (aincr*dabs(dcdv(k+3,ii))+aincr))
- enddo
- write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
- & i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3)
- write (iout,'(a)')
- enddo
- do j=1,nres
- do k=1,3
- dc(k,j)=temp(k,j)
- enddo
- enddo
- phi(i+3)=phii
- enddo
- return
- end
C----------------------------------------------------------------------------
subroutine check_ecart
C Check the gradient of the energy in Cartesian coordinates.
double precision energia(0:n_ene),energia1(0:n_ene)
double precision aincr2,etot,etot1,etot2
double precision dist,alpha,beta
+ double precision funcgrad,ff
+ external funcgrad
integer nf
integer uiparm(1)
double precision urparm(1)
call etotal(energia(0))
etot=energia(0)
call enerprint(energia(0))
+#ifdef LBFGS
+ ff=funcgrad(x,g)
+#else
call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
+#endif
icall =1
do i=1,nres
write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
enddo
return
end
-#ifdef FIVEDIAG
-!-----------------------------------------------------------------------------
- subroutine check_ecartint
-! Check the gradient of the energy in Cartesian coordinates.
- implicit none
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.CHAIN'
- include 'COMMON.INTERACT'
- include 'COMMON.DERIV'
- include 'COMMON.IOUNITS'
- include 'COMMON.VAR'
- include 'COMMON.CONTACTS'
- include 'COMMON.MD'
- include 'COMMON.LOCAL'
- include 'COMMON.SPLITELE'
- integer icall
- common /srutu/ icall
- double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
- & x(maxvar),g(maxvar)
- double precision dcnorm_safe(3),dxnorm_safe(3)
- double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
- double precision phi_temp(maxres),theta_temp(maxres),
- & alph_temp(maxres),omeg_temp(maxres)
- double precision ddc1(3),ddcn(3),dcnorm_safe1(3),dcnorm_safe2(3)
- double precision energia(0:n_ene),energia1(0:n_ene)
- integer uiparm(1)
- double precision urparm(1)
- double precision fdum
- external fdum
- integer i,j,k,nf
- double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
- double precision dist,alpha,beta
- icg=1
- nf=0
- nfl=0
- call intout
-! call intcartderiv
-! call checkintcartgrad
- call zerograd
- aincr=1.0D-5
- write(iout,*) 'Calling CHECK_ECARTINT.'
- nf=0
- icall=0
- write (iout,*) "Before geom_to_var"
- call geom_to_var(nvar,x)
- write (iout,*) "after geom_to_var"
- write (iout,*) "split_ene ",split_ene
- call flush(iout)
- if (.not.split_ene) then
- write(iout,*) 'Calling CHECK_ECARTINT if'
- call etotal(energia)
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
- etot=energia(0)
- write (iout,*) "etot",etot
- call enerprint(energia(0))
- call flush(iout)
-!el call enerprint(energia)
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
- call flush(iout)
- write (iout,*) "enter cartgrad"
- call flush(iout)
- call cartgrad
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
- write (iout,*) "exit cartgrad"
- call flush(iout)
- icall =1
- do i=1,nres
- write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
- enddo
- do j=1,3
- grad_s(j,0)=gcart(j,0)
- enddo
-!elwrite(iout,*) 'Calling CHECK_ECARTINT if'
- do i=1,nres
- do j=1,3
- grad_s(j,i)=gcart(j,i)
- grad_s(j+3,i)=gxcart(j,i)
- enddo
- enddo
- else
- write(iout,*) 'Calling CHECK_ECARTIN else.'
-!- split gradient check
- call zerograd
- call etotal_long(energia)
- call enerprint(energia(0))
-!el call enerprint(energia)
- call flush(iout)
- write (iout,*) "enter cartgrad"
- call flush(iout)
- call cartgrad
- write (iout,*) "exit cartgrad"
- call flush(iout)
- icall =1
- write (iout,*) "longrange grad"
- 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
- do j=1,3
- grad_s(j,0)=gcart(j,0)
- enddo
- do i=1,nres
- do j=1,3
- grad_s(j,i)=gcart(j,i)
- grad_s(j+3,i)=gxcart(j,i)
- enddo
- enddo
- call zerograd
- call etotal_short(energia)
- call enerprint(energia(0))
- call flush(iout)
- write (iout,*) "enter cartgrad"
- call flush(iout)
- call cartgrad
- write (iout,*) "exit cartgrad"
- call flush(iout)
- icall =1
- write (iout,*) "shortrange grad"
- 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
- do j=1,3
- grad_s1(j,0)=gcart(j,0)
- enddo
- do i=1,nres
- do j=1,3
- grad_s1(j,i)=gcart(j,i)
- grad_s1(j+3,i)=gxcart(j,i)
- enddo
- enddo
- endif
- write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
-! do i=1,nres
-c do i=nnt,nct
- do i=1,nres
- do j=1,3
- if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1)
- if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres)
- ddc(j)=c(j,i)
- ddx(j)=c(j,i+nres)
- dcnorm_safe1(j)=dc_norm(j,i-1)
- dcnorm_safe2(j)=dc_norm(j,i)
- dxnorm_safe(j)=dc_norm(j,i+nres)
- enddo
- do j=1,3
- c(j,i)=ddc(j)+aincr
- if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr
- if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr
- if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
- dc(j,i)=c(j,i+1)-c(j,i)
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- call int_from_cart1(.false.)
- if (.not.split_ene) then
- call etotal(energia1)
- etot1=energia1(0)
-c write (iout,*) "ij",i,j," etot1",etot1
- else
-!- split gradient
- call etotal_long(energia1)
- etot11=energia1(0)
- call etotal_short(energia1)
- etot12=energia1(0)
- endif
-!- end split gradient
-! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
- c(j,i)=ddc(j)-aincr
- if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr
- if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr
- if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
- dc(j,i)=c(j,i+1)-c(j,i)
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- call int_from_cart1(.false.)
- if (.not.split_ene) then
- call etotal(energia1)
- etot2=energia1(0)
-c write (iout,*) "ij",i,j," etot2",etot2
- ggg(j)=(etot1-etot2)/(2*aincr)
- else
-!- split gradient
- call etotal_long(energia1)
- etot21=energia1(0)
- ggg(j)=(etot11-etot21)/(2*aincr)
- call etotal_short(energia1)
- etot22=energia1(0)
- ggg1(j)=(etot12-etot22)/(2*aincr)
-!- end split gradient
-! write (iout,*) "etot21",etot21," etot22",etot22
- endif
-! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
- c(j,i)=ddc(j)
- if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)
- if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)
- if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1)
- dc(j,i)=c(j,i+1)-c(j,i)
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- dc_norm(j,i-1)=dcnorm_safe1(j)
- dc_norm(j,i)=dcnorm_safe2(j)
- dc_norm(j,i+nres)=dxnorm_safe(j)
- enddo
- do j=1,3
- c(j,i+nres)=ddx(j)+aincr
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- call int_from_cart1(.false.)
- if (.not.split_ene) then
- call etotal(energia1)
- etot1=energia1(0)
- else
-!- split gradient
- call etotal_long(energia1)
- etot11=energia1(0)
- call etotal_short(energia1)
- etot12=energia1(0)
- endif
-!- end split gradient
- c(j,i+nres)=ddx(j)-aincr
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- call int_from_cart1(.false.)
- if (.not.split_ene) then
- call etotal(energia1)
- etot2=energia1(0)
- ggg(j+3)=(etot1-etot2)/(2*aincr)
- else
-!- split gradient
- call etotal_long(energia1)
- etot21=energia1(0)
- ggg(j+3)=(etot11-etot21)/(2*aincr)
- call etotal_short(energia1)
- etot22=energia1(0)
- ggg1(j+3)=(etot12-etot22)/(2*aincr)
-!- end split gradient
- endif
-! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
- c(j,i+nres)=ddx(j)
- dc(j,i+nres)=c(j,i+nres)-c(j,i)
- dc_norm(j,i+nres)=dxnorm_safe(j)
- call int_from_cart1(.false.)
- enddo
- write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
- & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
- if (split_ene) then
- write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
- & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
- & k=1,6)
- write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
- & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
- & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
- endif
- enddo
- return
- end
-#else
c----------------------------------------------------------------------------
subroutine check_ecartint
C Check the gradient of the energy in Cartesian coordinates.
c aincr=1.0D-7
print '("Calling CHECK_ECARTINT",1pd12.3)',aincr
write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr
+ call cartprint
nf=0
icall=0
call geom_to_var(nvar,x)
call etotal_long(energia(0))
call enerprint(energia(0))
call flush(iout)
- write (iout,*) "enter cartgrad"
- call flush(iout)
+c write (iout,*) "enter cartgrad"
+c call flush(iout)
call cartgrad
- write (iout,*) "exit cartgrad"
- call flush(iout)
+c write (iout,*) "exit cartgrad"
+c call flush(iout)
icall =1
write (iout,*) "longrange grad"
do i=1,nres
- write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+ write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
& (gxcart(j,i),j=1,3)
enddo
do j=1,3
call etotal_short(energia(0))
call enerprint(energia(0))
call flush(iout)
- write (iout,*) "enter cartgrad"
- call flush(iout)
+c write (iout,*) "enter cartgrad"
+c call flush(iout)
call cartgrad
- write (iout,*) "exit cartgrad"
- call flush(iout)
+c write (iout,*) "exit cartgrad"
+c call flush(iout)
icall =1
write (iout,*) "shortrange grad"
do i=1,nres
- write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+ write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
& (gxcart(j,i),j=1,3)
enddo
do j=1,3
enddo
enddo
endif
+c write (iout,*) "Vector dc"
+c do i=0,nres
+c write (iout,'(i5,2(3f10.5,5x))')
+c & i,(dc(j,i),j=1,3),(dc(j,i+nres),j=1,3)
+c enddo
+c write (iout,*) "Coordinates after chainbuild_cart"
+c call chainbuild_cart
+c call cartprint
write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
do i=0,nres
print *,i
enddo
return
end
-#endif
c-------------------------------------------------------------------------
subroutine int_from_cart1(lprn)
implicit none
character*6 key
double precision fdum
external fdum
+ double precision funcgrad,ff
+ external funcgrad
integer i,ii,nf
double precision xi,etot,etot1,etot2
call zerograd
nf=1
nfl=3
cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
+c write (iout,*) "Before gradient"
+c call flush(iout)
+#ifdef LBFGS
+ ff=funcgrad(x,gana)
+#else
call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
+#endif
+c write (iout,*) "After gradient"
+c call flush(iout)
cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
icall=1
do i=1,nvar
call etotal(energia2(0))
etot2=energia2(0)
gg(i)=(etot2-etot1)/aincr
- write (iout,*) i,etot1,etot2
+c write (iout,*) i,etot1,etot2
x(i)=xi
enddo
write (iout,'(/2a)')' Variable Numerical Analytical',
projects / unres.git / blobdiff