I've got a problem with the implementation of the Powell method: I've created a box which contains a certain amount of circles (in a simple case all of the same diameter). The positions are randomally allocated, so it is possible that they overlap.
Furthermore I wrote an energy function which draws the overlap as a feather energy.
The powell routine shall minimize the energy to find an optimal configuration. At first sight, it does, but when I let the program show me the forces on every particle, these are not all equal zero. It turns out that the circles at the left and bottom border are those particles with net force != 0.
Does anyone have an idea? All that I do is giving the function to calculate the energy, isn't it?

double TotalEnergy(VecDoub_I &particles)
{
double totalEnergy = 0.;
double distance;
double k = federkonstante;
int N = 3*amount;
double dx,dy;
double sigma=100;

for (int i=0; i<N; i+=3) {
for (int j=0; j<N; j+=3) {

dx = particles[j]-particles[i];
dy = particles[j+1]-particles[i+1];

while (dx < -X_MAX/2) dx = X_MAX + dx;
while (dx > X_MAX/2) dx = X_MAX - dx;
while (dy > Y_MAX/2) dy = Y_MAX - dy;
while (dy < -Y_MAX/2) dy = Y_MAX + dy;

distance = sqrt(pow(dx, 2.)+pow(dy, 2.));

if (distance <= sigma) {
totalEnergy+= k*pow((sigma-distance)/(2.*a), 2.);
}

}
}
return totalEnergy;
}

The particles vector is a vector in which all the x,y positions of every particle are written one after another

ok, i posted that questions several weeks ago and still there's no answer. maybe my code isn't that clear. so i will post a more structured one, perhaps it's better understandable like that:

int main (int argc, const char * argv[])
{
char filename[50] = "MinPacking.txt";

double volume = RAD*RAD*M_PI;
length = sqrt(volume/PACKFRAC*AMOUNT);

VecDoub p(2*AMOUNT);

//random positions
for (int i=0; i<2*AMOUNT; i++) {
p[i] = length*static_cast<double>(rand())/RAND_MAX;
}


//powell minimization
Powell<Doub (VecDoub_I &)> powell(TotalEnergy);
p=powell.minimize(p);


//save optimal positions to file
SavePositions(filename, p);

//calculate force on every particle
Force(p);
return 0;
}

Doub TotalEnergy(VecDoub_I &x)
{
double totalEnergy = 0.,
distance, dx, dy;



for (int i=0; i<2*AMOUNT; i+=2) {
for (int j=0; j<2*AMOUNT; j+=2) {

dx = x[j]-x[i];
dy = x[j+1]-x[i+1];

while (dx > length/2) dx = length-dx;
while (dx < -length/2) dx = length +dx;
while (dy > length/2) dy = length-dy;
while (dy < -length/2) dy = length +dy;

distance = sqrt(dx*dx+dy*dy);

if (distance < 2*RAD) {
totalEnergy += k*pow(distance-2*RAD,2.);
}

}
}

return totalEnergy;
}

Again my problem: if i let this code run, it does find positions which are near to the minimum, but not exactly because the force on every particle is !=0 (near to 0 but not exactly). so i wonder whether this is my fault or if it is a numerical error which i have to accept. firstly i want to know the reason for that...
thanks