X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Fcheckder_p.F;fp=source%2Funres%2Fsrc_CSA_DiL%2Fcheckder_p.F;h=0000000000000000000000000000000000000000;hp=3fc2e62d5b91c3fa55b2911a04c234ba5a38e474;hb=de4bc5453ea46e111d936cb85e1758ed21c08fcd;hpb=b75425747e3e2b448ca5e0ef8367712e1f339124 diff --git a/source/unres/src_CSA_DiL/checkder_p.F b/source/unres/src_CSA_DiL/checkder_p.F deleted file mode 100644 index 3fc2e62..0000000 --- a/source/unres/src_CSA_DiL/checkder_p.F +++ /dev/null @@ -1,694 +0,0 @@ - subroutine check_cartgrad -C Check the gradient of Cartesian coordinates in internal coordinates. - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.DERIV' - dimension temp(6,maxres),xx(3),gg(3) - indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1 -* -* Check the gradient of the virtual-bond and SC vectors in the internal -* coordinates. -* - aincr=1.0d-7 - aincr2=5.0d-8 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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. - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CHAIN' - include 'COMMON.DERIV' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CONTACTS' -#ifdef MOMENT - include 'COMMON.CONTACTS.MOMENT' -#endif - common /srutu/ icall - dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar) - dimension grad_s(6,maxres) - double precision energia(0:n_ene),energia1(0:n_ene) - integer uiparm(1) - double precision urparm(1) - external fdum - icg=1 - nf=0 - nfl=0 - call zerograd - aincr=1.0D-7 - print '(a)','CG processor',me,' calling CHECK_CART.' - nf=0 - icall=0 - call geom_to_var(nvar,x) - call etotal(energia(0)) - etot=energia(0) - call enerprint(energia(0)) - call gradient(nvar,x,nf,g,uiparm,urparm,fdum) - icall =1 - do i=1,nres - write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) - enddo - do i=1,nres - do j=1,3 - grad_s(j,i)=gradc(j,i,icg) - grad_s(j+3,i)=gradx(j,i,icg) - enddo - enddo - call flush(iout) - write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' - do i=1,nres - do j=1,3 - xx(j)=c(j,i+nres) - ddc(j)=dc(j,i) - ddx(j)=dc(j,i+nres) - enddo - do j=1,3 - dc(j,i)=dc(j,i)+aincr - do k=i+1,nres - c(j,k)=c(j,k)+aincr - c(j,k+nres)=c(j,k+nres)+aincr - enddo - call etotal(energia1(0)) - etot1=energia1(0) - ggg(j)=(etot1-etot)/aincr - dc(j,i)=ddc(j) - do k=i+1,nres - c(j,k)=c(j,k)-aincr - c(j,k+nres)=c(j,k+nres)-aincr - enddo - enddo - do j=1,3 - c(j,i+nres)=c(j,i+nres)+aincr - dc(j,i+nres)=dc(j,i+nres)+aincr - call etotal(energia1(0)) - etot1=energia1(0) - ggg(j+3)=(etot1-etot)/aincr - c(j,i+nres)=xx(j) - dc(j,i+nres)=ddx(j) - enddo - write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)') - & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6) - enddo - return - end -c---------------------------------------------------------------------------- - subroutine check_ecartint -C Check the gradient of the energy in Cartesian coordinates. - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CONTROL' - include 'COMMON.CHAIN' - include 'COMMON.DERIV' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CONTACTS' -#ifdef MOMENT - include 'COMMON.CONTACTS.MOMENT' -#endif - include 'COMMON.MD_' - include 'COMMON.LOCAL' - include 'COMMON.SPLITELE' - common /srutu/ icall - dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar), - & g(maxvar) - dimension dcnorm_safe(3),dxnorm_safe(3) - dimension 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 energia(0:n_ene),energia1(0:n_ene) - integer uiparm(1) - double precision urparm(1) - external fdum - r_cut=2.0d0 - rlambd=0.3d0 - icg=1 - nf=0 - nfl=0 - call intout -c call intcartderiv -c call checkintcartgrad - call zerograd - aincr=1.0D-5 - write(iout,*) 'Calling CHECK_ECARTINT.' - nf=0 - icall=0 - call geom_to_var(nvar,x) - if (.not.split_ene) then - call etotal(energia(0)) - etot=energia(0) - 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 - 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 - 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 -!- split gradient check - write (iout,*) "split_ene not supported" -c call zerograd -c call etotal_long(energia(0)) -c call enerprint(energia(0)) -c call flush(iout) -c write (iout,*) "enter cartgrad" -c call flush(iout) -c call cartgrad -c write (iout,*) "exit cartgrad" -c call flush(iout) -c icall =1 -c write (iout,*) "longrange grad" -c do i=1,nres -c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), -c & (gxcart(j,i),j=1,3) -c enddo -c do j=1,3 -c grad_s(j,0)=gcart(j,0) -c enddo -c do i=1,nres -c do j=1,3 -c grad_s(j,i)=gcart(j,i) -c grad_s(j+3,i)=gxcart(j,i) -c enddo -c enddo -c call zerograd -c call etotal_short(energia(0)) -c call enerprint(energia(0)) -c call flush(iout) -c write (iout,*) "enter cartgrad" -c call flush(iout) -c call cartgrad -c write (iout,*) "exit cartgrad" -c call flush(iout) -c icall =1 -c write (iout,*) "shortrange grad" -c do i=1,nres -c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), -c & (gxcart(j,i),j=1,3) -c enddo -c do j=1,3 -c grad_s1(j,0)=gcart(j,0) -c enddo -c do i=1,nres -c do j=1,3 -c grad_s1(j,i)=gcart(j,i) -c grad_s1(j+3,i)=gxcart(j,i) -c enddo -c enddo - endif - write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' - do i=0,nres - do j=1,3 - xx(j)=c(j,i+nres) - ddc(j)=dc(j,i) - ddx(j)=dc(j,i+nres) - do k=1,3 - dcnorm_safe(k)=dc_norm(k,i) - dxnorm_safe(k)=dc_norm(k,i+nres) - enddo - enddo - do j=1,3 - dc(j,i)=ddc(j)+aincr - call chainbuild_cart -#ifdef MPI -c Broadcast the order to compute internal coordinates to the slaves. -c if (nfgtasks.gt.1) -c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR) -#endif -c call int_from_cart1(.false.) - if (.not.split_ene) then - call etotal(energia1(0)) - etot1=energia1(0) - else -!- split gradient -c call etotal_long(energia1(0)) -c etot11=energia1(0) -c call etotal_short(energia1(0)) -c etot12=energia1(0) -c write (iout,*) "etot11",etot11," etot12",etot12 - endif -!- end split gradient -c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 - dc(j,i)=ddc(j)-aincr - call chainbuild_cart -c call int_from_cart1(.false.) - if (.not.split_ene) then - call etotal(energia1(0)) - etot2=energia1(0) - ggg(j)=(etot1-etot2)/(2*aincr) - else -!- split gradient -c call etotal_long(energia1(0)) -c etot21=energia1(0) -c ggg(j)=(etot11-etot21)/(2*aincr) -c call etotal_short(energia1(0)) -c etot22=energia1(0) -c ggg1(j)=(etot12-etot22)/(2*aincr) -!- end split gradient -c write (iout,*) "etot21",etot21," etot22",etot22 - endif -c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 - dc(j,i)=ddc(j) - call chainbuild_cart - enddo - do j=1,3 - dc(j,i+nres)=ddx(j)+aincr - call chainbuild_cart -c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm" -c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3) -c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3) -c write (iout,*) "dxnormnorm",dsqrt( -c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2) -c write (iout,*) "dxnormnormsafe",dsqrt( -c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) -c write (iout,*) - if (.not.split_ene) then - call etotal(energia1(0)) - etot1=energia1(0) - else -!- split gradient -c call etotal_long(energia1(0)) -c etot11=energia1(0) -c call etotal_short(energia1(0)) -c etot12=energia1(0) - endif -!- end split gradient -c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 - dc(j,i+nres)=ddx(j)-aincr - call chainbuild_cart -c write (iout,*) "i",i," j",j," dxnorm- and dxnorm" -c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3) -c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3) -c write (iout,*) -c write (iout,*) "dxnormnorm",dsqrt( -c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2) -c write (iout,*) "dxnormnormsafe",dsqrt( -c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) - if (.not.split_ene) then - call etotal(energia1(0)) - etot2=energia1(0) - ggg(j+3)=(etot1-etot2)/(2*aincr) - else -!- split gradient -c call etotal_long(energia1(0)) -c etot21=energia1(0) -c ggg(j+3)=(etot11-etot21)/(2*aincr) -c call etotal_short(energia1(0)) -c etot22=energia1(0) -c ggg1(j+3)=(etot12-etot22)/(2*aincr) -!- end split gradient - endif -c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 - dc(j,i+nres)=ddx(j) - call chainbuild_cart - 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 -c------------------------------------------------------------------------- - subroutine int_from_cart1(lprn) - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' -#ifdef MPI - include 'mpif.h' - integer ierror -#endif - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.CHAIN' - include 'COMMON.GEO' - include 'COMMON.INTERACT' - include 'COMMON.LOCAL' - include 'COMMON.NAMES' - include 'COMMON.SETUP' - include 'COMMON.TIME1' - logical lprn - if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates' -#ifdef TIMING - time01=MPI_Wtime() -#endif -#if defined(PARINT) && defined(MPI) - do i=iint_start,iint_end -#else - do i=2,nres -#endif - dnorm1=dist(i-1,i) - dnorm2=dist(i,i+1) - do j=1,3 - c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 - & +(c(j,i+1)-c(j,i))/dnorm2) - enddo - be=0.0D0 - if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1) - omeg(i)=beta(nres+i,i,maxres2,i+1) - alph(i)=alpha(nres+i,i,maxres2) - theta(i+1)=alpha(i-1,i,i+1) - vbld(i)=dist(i-1,i) - vbld_inv(i)=1.0d0/vbld(i) - vbld(nres+i)=dist(nres+i,i) - if (itype(i).ne.10) then - vbld_inv(nres+i)=1.0d0/vbld(nres+i) - else - vbld_inv(nres+i)=0.0d0 - endif - enddo -#if defined(PARINT) && defined(MPI) - if (nfgtasks1.gt.1) then -cd write(iout,*) "iint_start",iint_start," iint_count", -cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ", -cd & (iint_displ(i),i=0,nfgtasks-1) -cd write (iout,*) "Gather vbld backbone" -cd call flush(iout) - time00=MPI_Wtime() - call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -cd write (iout,*) "Gather vbld_inv" -cd call flush(iout) - call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -cd write (iout,*) "Gather vbld side chain" -cd call flush(iout) - call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -cd write (iout,*) "Gather vbld_inv side chain" -cd call flush(iout) - call MPI_Allgatherv(vbld_inv(iint_start+nres), - & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1), - & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -cd write (iout,*) "Gather theta" -cd call flush(iout) - call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -cd write (iout,*) "Gather phi" -cd call flush(iout) - call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -#ifdef CRYST_SC -cd write (iout,*) "Gather alph" -cd call flush(iout) - call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -cd write (iout,*) "Gather omeg" -cd call flush(iout) - call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1), - & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0), - & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) -#endif - time_gather=time_gather+MPI_Wtime()-time00 - endif -#endif - do i=1,nres-1 - do j=1,3 - dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) - enddo - enddo - do i=2,nres-1 - do j=1,3 - dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) - enddo - enddo - if (lprn) then - do i=2,nres - write (iout,1212) restyp(itype(i)),i,vbld(i), - &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i), - &rad2deg*alph(i),rad2deg*omeg(i) - enddo - endif - 1212 format (a3,'(',i3,')',2(f15.10,2f10.2)) -#ifdef TIMING - time_intfcart=time_intfcart+MPI_Wtime()-time01 -#endif - return - end -c---------------------------------------------------------------------------- - subroutine check_eint -C Check the gradient of energy in internal coordinates. - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CHAIN' - include 'COMMON.DERIV' - include 'COMMON.IOUNITS' - include 'COMMON.VAR' - include 'COMMON.GEO' - common /srutu/ icall - dimension x(maxvar),gana(maxvar),gg(maxvar) - integer uiparm(1) - double precision urparm(1) - double precision energia(0:n_ene),energia1(0:n_ene), - & energia2(0:n_ene) - character*6 key - external fdum - call zerograd - aincr=1.0D-7 - print '(a)','Calling CHECK_INT.' - nf=0 - nfl=0 - icg=1 - call geom_to_var(nvar,x) - call var_to_geom(nvar,x) - call chainbuild - icall=1 - print *,'ICG=',ICG - call etotal(energia(0)) - etot = energia(0) - call enerprint(energia(0)) - print *,'ICG=',ICG -#ifdef MPL - if (MyID.ne.BossID) then - call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID) - nf=x(nvar+1) - nfl=x(nvar+2) - icg=x(nvar+3) - endif -#endif - nf=1 - nfl=3 -cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) - call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) -cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar) - icall=1 - do i=1,nvar - xi=x(i) - x(i)=xi-0.5D0*aincr - call var_to_geom(nvar,x) - call chainbuild - call etotal(energia1(0)) - etot1=energia1(0) - x(i)=xi+0.5D0*aincr - call var_to_geom(nvar,x) - call chainbuild - call etotal(energia2(0)) - etot2=energia2(0) - gg(i)=(etot2-etot1)/aincr - write (iout,*) i,etot1,etot2 - x(i)=xi - enddo - write (iout,'(/2a)')' Variable Numerical Analytical', - & ' RelDiff*100% ' - do i=1,nvar - if (i.le.nphi) then - ii=i - key = ' phi' - else if (i.le.nphi+ntheta) then - ii=i-nphi - key=' theta' - else if (i.le.nphi+ntheta+nside) then - ii=i-(nphi+ntheta) - key=' alpha' - else - ii=i-(nphi+ntheta+nside) - key=' omega' - endif - write (iout,'(i3,a,i3,3(1pd16.6))') - & i,key,ii,gg(i),gana(i), - & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr) - enddo - return - end