X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fcheckder_p.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fcheckder_p.F;h=0539e4825f9b00e33135f03a25b022f77bbccdd1;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/checkder_p.F b/source/unres/src_MD-M-newcorr/checkder_p.F new file mode 100644 index 0000000..0539e48 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/checkder_p.F @@ -0,0 +1,700 @@ + 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' + 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' + 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 + call zerograd + call etotal_long(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 + 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(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 + 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=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 + call etotal_long(energia1(0)) + etot11=energia1(0) + call etotal_short(energia1(0)) + 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 + call etotal_long(energia1(0)) + etot21=energia1(0) + ggg(j)=(etot11-etot21)/(2*aincr) + call etotal_short(energia1(0)) + etot22=energia1(0) + 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 + call etotal_long(energia1(0)) + etot11=energia1(0) + call etotal_short(energia1(0)) + 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 + call etotal_long(energia1(0)) + etot21=energia1(0) + ggg(j+3)=(etot11-etot21)/(2*aincr) + call etotal_short(energia1(0)) + etot22=energia1(0) + 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) then + if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1) + if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then + tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres) + endif + if (itype(i-1).ne.10) then + tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1) + omicron(1,i)=alpha(i-2,i-1,i-1+nres) + omicron(2,i)=alpha(i-1+nres,i-1,i) + endif + if (itype(i).ne.10) then + tauangle(2,i+1)=beta(i-2,i-1,i,i+nres) + endif + endif + 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