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' include 'COMMON.SCCOR' 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' include 'COMMON.SCCOR' 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' include 'COMMON.SCCOR' 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)) c do i=1,nres c write (iout,*) "atu?", gloc_sc(1,i,icg),gloc(i,icg) c enddo etot=energia(0) call enerprint(energia(0)) call flush(iout) write (iout,*) "enter cartgrad" c do i=1,nres c write (iout,*) gloc_sc(1,i,icg) c enddo 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)) c do i=1,nres c write (iout,*) gloc_sc(1,i,icg) c enddo 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+1 #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