Program stops unexpectedly

[lmd141]

Hi, there's a problem I've been having for a while now that I can't figure out. I have the following function:

void MCinsert(int &naccept, int &nreject, const double *box_size, double **coords, int &n_atoms, double &energy_box, CRandomMersenne* pRanGen)
{
	
		bool accept = false;
		int new_atom=n_atoms+1;

		// Select a random position to insert particle
		double x_position = rand(0, volume);
		double y_position = rand(0, volume);
                double z_position = rand(0, volume);

		// Insert a particle
		coords[new_atom][0] = x_position;
		coords[new_atom][1] = y_position;
		coords[new_atom][2] = z_position;
	
		// Calculate initial energy
		double energy1 = 0;
		energy1 = LJ_pot(new_atom, box_size, coords, n_atoms);

		// Acceptance Criteria 
		double acceptance_insert = (volume/(pow(debroglie,3)*(n_atoms+1)))*exp(beta*(chem_pot - (energy1)));
		double random1;
		random1 = pRanGen->Random();

			if (random1 < acceptance_insert)
			{
				accept = true;
				naccept = naccept + 1;
				n_atoms = n_atoms + 1;
				energy_box = energy_box + (energy1);		
			}

			else
			{
				accept = false;
				nreject = nreject + 1;
				coords[new_atom][0] = 0;
				coords[new_atom][1] = 0;
				coords[new_atom][2] = 0;
			}

Which basically decides whether or not to accept or reject a move based on the energy that has been calculated and an acceptance criteria. As you can see there is another function "LJ_pot", which is the following:

double LJ_pot(const int id_particle, const double *box_size, double **coords, int n_atoms)
{

	double particle_energy = 0; 
	for (int i = 0; i <= n_atoms-1; i++)
    {		
			double delta_x = coords[id_particle][0] - coords[i][0];
            double delta_y = coords[id_particle][1] - coords[i][1];
            double delta_z = coords[id_particle][2] - coords[i][2];

            // Apply periodic boundaries
            delta_x = make_periodic(delta_x, box_size[0]);
            delta_y = make_periodic(delta_y, box_size[1]);
            delta_z = make_periodic(delta_z, box_size[2]);

			// Calculate the LJ potential

            double r = pow((delta_x*delta_x) + (delta_y*delta_y) +
                          (delta_z*delta_z),0.5)/sigma;

			if (r>0)
			{
			double e_lj = ((4.0/pow(r,12.0))-(4.0/pow(r,6.0)));
			particle_energy = (particle_energy + e_lj)/epsilon;
			}
            
	}
	
	return particle_energy;
}

Now the problem is that after 1 move (which is accepted) the program just stops. No errors or breakpoints. This is obviously a bit confusing.

I should point out the problem seems to be with n_atoms. If I remove the & from the function then the program runs (but obviously n_atoms is not updated).

Any help would be appreciated, and if any other information is required just let me know.

[TheIdeasMan]

Hi lmd141,

I don't see an immediate answer, so I will provide these trivial comments instead:

int new_atom=n_atoms+1;

This nearly confused me, because the variables names are so close. Try to avoid doing this, it will save you one day.

nreject = nreject + 1; can be written nreject++; because it is an int. If it were a double you can do this MyDouble += 1.0; or with another variable MyDouble += AnotherDbl;

I tend to name my variables like this ParticleEnergy - no underscores and capitalise the first letter of each word. That is my preference, just putting the idea out there.

for (int i = 0; i <= n_atoms-1; i++)

can be written:

for (int i = 0; i < n_atoms; i++)


With these parts:
if (random1 < acceptance_insert) if (r>0)

Doubles are stored as binary fractions, and cannot represent al real numbers, so comparisons like that often fail. Consider this:

#include <cmath>
#include <limits>

float a = 0.1; //a == 0.0999997 
float b = 10 * a; //b == 0.9999997
if ( b == 1.0){} //false
if ((1.0 - b) == 0.0){} //false
if((1.0 -b) < 0.0){} //false
if((1.0 - b) > 0.0){} //true but not what you want

MyEpsilon = std::numeric_limits<float>epsilon;

if(abs(1.0 - b) < MyEpsilon){}  // effective equality comparison with zero

double c = 1000.0;
//scale epsilon up to the number
MyDblEpsilon = c* std::numeric_limits<double>epsilon; //only good for values from 1.0e3 to < 1.0e4

double d = 2 * 500.0;
if (c == d){} //fail
if(abs(c-d )< MyDblEpsilon){}  // effective equality comparison betwen variables c and d

The value of epsilon needs to be scaled up because you need to deal with the 'distance' between representable numbers which varies across the range of floating point values.

Changing the type to double doesn't fix the problem.

You can do something similar for < and > comparisons - I am sure you can figure that out.

Hopefully this won't be confused with your global variable epsilon.

This whole thing could be the cause of your problem in the if statement in the MCinsert and LJ_pot functions. If it isn't, then what I have mentioned is better practice anyway.

Finally, the way to solve these problems properly is to use the debugger. If you have an IDE it should be easy to use. You can set break points, have a watch list of variable values,step through the code 1 line at a time. See how the values change, and how that affects the logic.


Hope all this helps. Good luck !!

[cire]

I should point out the problem seems to be with n_atoms. If I remove the & from the function then the program runs (but obviously n_atoms is not updated).

This makes me curious to see how the original int is bound to n_atoms in the function call of MCInsert.