### Program stops unexpectedly

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

 ``1234567891011121314151617181920212223242526272829303132333435363738394041`` ``````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:

 ``123456789101112131415161718192021222324252627282930`` ``````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.
Last edited on
Hi lmd141,

`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:

 ``1234567891011121314151617181920212223`` ``````#include #include 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_limitsepsilon; 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_limitsepsilon; //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 !!

Last edited on
 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.
Topic archived. No new replies allowed.