A method of determining one or more stable crystal structures of a system, each crystal structure comprising a repeated unit cell, by calculating a configurational energy of the system with periodic boundary conditions, the configurational energy resulting from non-electrostatic interactions and electrostatic interactions, the system having a plurality of particles, wherein the particles comprise atoms and molecules, said method comprising the computer-implemented steps (i) to (ix) of: i) defining a cut-off radius, said radius defining a non-electrostatic interaction potential cut-off distance between the particles in the system; ii) defining a set of cell vectors to generate a simulation cell; iii) defining a set of supercell vectors to generate a supercell, wherein said supercell comprises a plurality of replicas of the simulation cell; iv) calculating, for any one of the particles located within the supercell, non-electrostatic pair potentials between the one of the particles and first surrounding ones of the particles located from the one of the particles within a boundary defined by the cut-off radius, said non-electrostatic pair potentials resulting from the interaction between the one of the particles with the first surrounding ones of the particles, wherein the first surrounding ones of the particles include ones of particles lying outside the supercell; v) summing the calculated non-electrostatic pair potentials for the particles located within the supercell to provide a non-electrostatic configurational energy of the system; vi) accounting for electrostatic interactions, wherein each particle comprises one or more particle multipoles and the electrostatic interactions between the particles comprise multipole interactions between the corresponding one or more particle multipoles that are combined using an Ewald sum; vii) calculating a global minimum by dividing the calculated configurational energy of the system by a total number of simulation cells to define a potential energy per simulation cell; and viii) determining the global minimum of the potential energy per simulation cell by obtaining, by gradient descent, a plurality of local minima of the configurational energy, and by selecting, using a basin-hopping global optimisation algorithm, the lowest one of the local minima to be the global minimum, and wherein the unit cell is equivalent to the simulation cell; and ix) selecting the crystal structure having the global minimum to be a candidate for the stable crystal structure; the method further comprising: determining potential properties of the selected stable crystal structure and accordingly controlling the necessary parameters for manufacturing the stable crystal structure i) defining a cut-off radius, said radius defining a non-electrostatic interaction potential cut-off distance between the particles in the system; ii) defining a set of cell vectors to generate a simulation cell; iii) defining a set of supercell vectors to generate a supercell, wherein said supercell comprises a plurality of replicas of the simulation cell; iv) calculating, for any one of the particles located within the supercell, non-electrostatic pair potentials between the one of the particles and first surrounding ones of the particles located from the one of the particles within a boundary defined by the cut-off radius, said non-electrostatic pair potentials resulting from the interaction between the one of the particles with the first surrounding ones of the particles, wherein the first surrounding ones of the particles include ones of particles lying outside the supercell; v) summing the calculated non-electrostatic pair potentials for the particles located within the supercell to provide a non-electrostatic configurational energy of the system; vi) accounting for electrostatic interactions, wherein each particle comprises one or more particle multipoles and the electrostatic interactions between the particles comprise multipole interactions between the corresponding one or more particle multipoles that are combined using an Ewald sum; vii) calculating a global minimum by dividing the calculated configurational energy of the system by a total number of simulation cells to define a potential energy per simulation cell; and viii) determining the global minimum of the potential energy per simulation cell by obtaining, by gradient descent, a plurality of local minima of the configurational energy, and by selecting, using a basin-hopping global optimisation algorithm, the lowest one of the local minima to be the global minimum, and wherein the unit cell is equivalent to the simulation cell; and ix) selecting the crystal structure having the global minimum to be a candidate for the stable crystal structure;