custom simmulated annealing program doesn't produce desired result

I need to perform global optimization using simulated annealing algorithm. To do that I prepared a custom simulated annealing program which is basically optimized a distance function by varying few set of parameters. The program runs without any problem but it doesn't optimize the parameters. The code is attached in the following section.

if true
  atype=a{1};
n=size(atype,1);
x= a{2};
y= a{3};
z= a{4};
qm_atom_charge=aa{3};
fclose(fileID1);
ah0=aaa{1};
bh0=aaa{2};
en0=aaa{3};
na=size(ah0,1);
nb=size(bh0,1);
bound=3.0;
vstep=1.0;
  bound1(1:na)=bound;
  bound2(1:nb)=bound;
  ahu=ah0+bound1';
  ahl=ah0-bound1';
  bhu=bh0+bound2';
  bhl=bh0-bound2';
  enu=en0+bound1';
  enl=en0-bound1';
  Temp=100;            
  iter1=0;
  acc=0.0;
  rejc=0.0;
      tmpstep=5;
      mcstep=5;
      maxit=tmpstep.*mcstep.*na.*na.*nb;
      cost_function(1:maxit)=0.0;
              ah=ah0;
              bh=bh0;
              en=en0;
              ah1=ah0;
              bh1=bh0;
              en1=en0;
              atom_charge=qeq_charge(n,x,y,z,atype,ah,bh,en);
              distance=atom_charge-qm_atom_charge;
              cost_current=norm(distance);
              cost_tol=9999;
for itemp=1:tmpstep
      Temp=Temp.*0.01;
  for iter=1:mcstep
      rv1=2.*bound.*rand(1:na)-bound;
      rv2=2.*bound.*rand(1:nb)-bound;
      rv3=2.*bound.*rand(1:na)-bound;
      rv1=vstep.*rv1;
      rv2=vstep.*rv2;
      rv3=vstep.*rv3;
      for m1=1:na
                 ah(m1)=ah1(m1)+rv1(m1);
                 if ah(m1)<=0.0
                     ah(m1)=ah(m1)-rv1(m1);
                 else
                 ind1=(ah(m1)>=ahl(m1)&ah(m1)<=ahu(m1));
                 if (ind1==0)
                      ah(m1)=ah(m1)-rv1(m1);
                 end
                 end
          for m2=1:nb  
                 bh(m2)=bh1(m2)+rv2(m2);
                 if bh(m2)<=0.0
                     bh(m2)=bh(m2)-rv2(m2);
                 else
                 ind2=(bh(m2)>=bhl(m2)&bh(m2)<=bhu(m2));
                 if (ind2==0)
                      bh(m2)=bh(m2)-rv2(m2);
                 end
                 end
              for m3=1:na 
                 en(m3)=en1(m3)+rv3(m3);
                 if en(m3)<=0.0
                     en(m3)=en(m3)-rv3(m3);
                 else
                 ind3=(en(m3)>=enl(m3)&en(m3)<=enu(m3));
                 if (ind3==0)
                      en(m3)=en(m3)-rv3(m3);
                 end
                 end

           iter1=iter1+1;
                chk(1)=isequal(ah,ah1);
                chk(2)=isequal(bh,bh1);
                chk(3)=isequal(en,en1);
            if any(chk)==0                   
            atom_charge=qeq_charge(n,x,y,z,atype,ah,bh,en);
            distance=atom_charge-qm_atom_charge;
            cost_old=cost_current;
            cost_current=norm(distance);
            if cost_current<=cost_old
            ah1=ah;
            bh1=bh;
            en1=en;
            acc=acc+1;
            cost_function(acc)=cost_current;
            if cost_current<=0.001
            cost_tol=1;   
            break
            else
            end
            else
            p1=rand(1);
            p=exp((cost_old-cost_current)/Temp);
            if p1<p
            acc=acc+1;
            ah1=ah;
            bh1=bh;
            en1=en;
            cost_function(acc)=cost_current;
            else
            rejc=rejc+1;
            end
                    if cost_tol==1
                    break
                    end 
            end
                    if cost_tol==1
                    break
                    end 

            end
                    if cost_tol==1
                    break
                    end         

               end
                    if cost_tol==1
                    break
                    end 
             end
                    if cost_tol==1
                    break
                    end 
       end
                    if cost_tol==1
                    break
                    end 
    end
                    if cost_tol==1
                    break
                    end 
end
  end

Here, the function qeq_charge(n,x,y,z,atype,ah,bh,en) outputs an array, atom_charge and then I calculate the distance function. My intention is to minimize the distance function by varying ah,bh,en parameters. However, it doesn't work properly. I have been really puzzled in this issue since last week. Please help me out in this regard.

Thanks in advance sudipta