I am creating a simple potential based on geometry optimised structures from GULP.

For each atom in the structure, I employ harmonic potentials for their nearest neighbours as well as three-body potentials. These potentials are using the equilibrium distances and angles of the structure such that the total lattice energy is 0.00000000 eV and there are no negative phonon frequencies.

The problem is that for certain structures, the bulk modulus in the ouptut file is 0 even if the structure is dynamically stable. Would you have advice on why this might be occurring? Most structures tested behave as expected with positive bulk moduli.

Thank you for your help.
Best wishes,
Quentin

P.S. Apologies - I am not able to upload files as a new user. I hope the screenshots help a bit.

Hi Quentin,
This is a fun one. You’ve actually found a hiccup in the code due to the unusual nature of the model. The problem is that you have only bonded potentials with manually specified connectivity. For all potentials GULP assigns a cutoff & in the case of the bonded ones it’s usually based on the sum of the covalent radii times a tolerance. Because your bonds are unphysically long, they exceed the cutoff and so the harmonic terms are neglected. There is a solution to this, which is to increase the tolerance on bonds. If you add:

rtol 2.0

then you should see that your interactions at ~5.9 Ang are now within the sum of the covalent radii plus tolerance. You could have also saved yourself lots of work by just using a distance-based harmonic with a cut-off rather than bonding everything (i.e. 1 potential to input and no bonds & so a lot quicker and simpler!).
Hope that clarifies things.
Regards,
Julian