annually , at least 17 , 000 metric tons of fluorine are produced . while it is cheap ($ 5 - 8 per kilogram ) in the form of solid compounds ( uf 6 and sf 6 ), fluorine gas is much more expensive because of difficulties associated with its storage and handling . elemental fluorine is widely used , for example : as uf 6 , used for producing nuclear fuels ; as sf 6 , a dielectric medium in high voltage transformers ; hexafluorides of re and w that are used in the electronics industry ; and various organic fluorides ( e . g . in pharmaceuticals ). elemental chlorine has a similar volume of production and also a wide range of applications , mostly in the production of organic compounds ( e . g . polyvinylchloride ). production of these compounds would be greatly facilitated if storage of f and cl could become safer and more efficient . recently , the possibility of existence of totally new and unexpected compounds at high pressures has been suggested and experimentally proven , where such compounds as nacl 7 , nacl 3 , na 3 cl 2 , na 2 cl , na 3 cl have been found . methods for predicting the feasibility of these and other compounds , as well as optimized crystal structures for these compounds , include those disclosed in u . s . patent application publication no . 2012 / 0330632 , the entire disclosure of which is incorporated by reference herein . in accordance with the present disclosure heavier alkali metals are used to form compounds with fluorine or chlorine that are stable at zero pressure . using the methods disclosed in u . s . patent application publication no . 2012 / 0330632 , kcl 3 , is identified as thermodynamically stable at zero pressure and zero kelvin temperature , with even greater stability for chlorine - rich compounds of rb and cs , and similarly interesting results obtained for fluorides of cesium ( fig1 ). briefly , the structures of the compounds of the present disclosure , utilizing the methods of u . s . patent application publication no . 2012 / 0330632 , were as follows . as disclosed therein , a variable - composition evolutionary algorithm , known as universal structure predictor : evolutionary xtallography ( uspex ), which initializes the population using randomly selected chemical compositions , with atoms placed at random within randomly selected space groups , is used . following relaxations , the structures developed are ranked using their fitness function , and the fittest ones are allowed to produce daughter structures using crossover ( heredity ) and mutation operators . compositions that are on the thermodynamic convex hull , are by definition thermodynamically stable . more specifically , to arrive at the compounds of the present disclosure following the methods of u . s . patent application publication no . 2012 / 0330632 , characterization information of the desired chemical structure was input . the characterization information includes chemical formulae , an atomic number of an atom , a number of atoms in a unit cell , valences of atoms , and unit cell volume . the characterization information may also include lattice parameters known from experimental data and other information . following the inputting of characterization information , a first generation of crystal structures is obtained , preferably by a random input , seed input or symmetrical initialization . the first generation of crystal structures is generated using a symmetrical initialization algorithm based on the input characterization information . a randomly produced crystal structure can be utilized to commence generation of crystal structures , with the symmetrical initialization algorithm randomly selecting a space group compatible with the number of atoms in a unit cell of the chemical structure . for example , the space group can be randomly selected from space group numbers 1 to 230 using a mersenne twister algorithm . see matsumoto , m ., et al ., mersenne twister : a 623 - dimensionally equidistributed uniform pseudo - random number generator , acm transactions on modelling and computer simulation , 8 ( 1 ) ( 1998 ), pp . 3 - 30 . the space group can also be randomly selected by a random number generator algorithm . the symmetrical initialization selects from either an entire list of two hundred and thirty space groups or selects from a user - provided list of space groups . a symmetrical initialization unit cell is then configured based on the selected space group and a predetermined initial volume . the symmetrical initialization preferably includes adding symmetrical random structures in each generation to improve population diversity and provide new crystal structures to be explored by the evolutionary algorithm . the symmetrical initialization includes unit cell splitting to provide a reduced number of degrees of freedom of the crystal structures . the symmetrical initialization introduces additional translational symmetry , i . e . pseudosymmetry , to the crystal structures . the symmetrical initialization uses up to two hundred and thirty space groups , allowing the user to produce a diverse first generation of crystal structures . using the symmetrical initialization algorithm , a plurality of atoms are randomly placed on separate first wyckoff positions , herein a first position , and coordinates of each of the atoms are multiplied by space group symmetry operations . the first position of each of the atoms includes all points of a space group for which site - symmetry groups are conjugate subgroups of the space group . if an inter - atomic distance of two or more symmetry - related atoms is within a user - defined threshold , the symmetry - related atoms are merged into one atom on a second position defined by averaging of coordinates of symmetry - related atoms . a sum of the merged atoms and remaining atoms outside of the user - defined threshold for the inter - atomic distance is compatible with the space group symmetry operation . however , discarding the merged atoms from a crystal structure of the first generation of crystal structures may be necessary if the sum is incompatible with the space group symmetry operation . new atoms are added one at a time to ensure that no atoms , including symmetrically unrelated atoms , are too closely positioned . to obtain an exact number of atoms necessary in a unit cell , certain atoms are placed on special wyckoff positions having a correct multiplicity . the multiplicity of the wyckoff position refers to a number of atoms obtained after multiplying the atom in the wyckoff position using a group symmetry element . for example , if six atoms are needed to reach the number of atoms in the unit cell , atoms are input into positions with multiplicities of 1 - 6 , but not with a multiplicity larger than six . otherwise , after multiplication , more atoms than are desired for the unit cell are obtained . this process of symmetrical initialization provides an unbiased procedure for sampling of each allowed space group and placing of atoms in positions that ensure minimum distances between all atoms . if symmetrically related atoms are closer than the user - defined threshold , they are merged into the second position . if symmetrically non - equivalent atoms are within a user - defined threshold for distance between atoms , the crystal structure is discarded without relaxation . the user - defined threshold is one of the user input parameters , which , in embodiments , is set to a half sum of the covalent radii of the corresponding atoms . local optimization of the chemical structure of the first generation of crystal structures according to the input parameters then occurs . local optimization , i . e . relaxation , does not break individual cell symmetry , but increases symmetry of the entire supergroup . preferably , random displacement of atomic positions provides a possibility of breaking symmetry . variation operators are also provided to break symmetry , and to remove bias that may have been caused by the selected space group . after optimization , the optimized first generation of crystal structures is input into a niching algorithm to select two separate groups of crystal structures from the optimized first generation of crystal structures . the niching algorithm calculates a set of conventional fitness functions and a set of fingerprint functions for each crystal structure of the first generation of crystal structures . the set of fitness functions allow the user to assess the quality of the crystal structures in the first generation of crystal structures according to a chemical property of the chemical structure that the user wants to optimize . chemical properties include , for example , free energy , density , hardness , dielectric constant and / or bandgap . the set of fingerprint functions allows the user to assess the similarity between each of the crystal structures of optimal group of crystal structures and the parent group of crystal structures . a fingerprint function of a crystal structure is calculated according to equation ( 1 ) found in u . s . patent application publication no . 2012 / 0330632 , reproduced below : where v is the unit cell volume , n b is a number of type b atoms in the unit cell , r ij is inter - atomic distance between atoms a i and b j , δ ( r − r ij ) is a gaussian - smeared delta function , which absorbs numerical errors and makes f ( r ) a smooth function . gaussian smearing is controlled by an input parameter sigmafing . f ( r ) is discretized over bins of width δ to obtain a fingerprint vector fp a i b . bins of width δ are preferably specified by a user as an input parameter deltafing . the fingerprint function input parameters are preferably set as 0 . 03 for sigmafing , 0 . 08 for deltafing , and 8 . 0 for rmaxfing . the fingerprint function of equation ( 1 ) is provided for an individual type a i atoms relative to all b j type atoms , referred to herein as atoms a i and b j , respectively , surrounding atom a i with a sum running over all b j atoms within a distance threshold r max , which preferably is user specified as an input parameter rmaxfing , with r max referring to a sum over all atoms b j for which r ij & lt ; r max . the niching algorithm includes determining a degree of similarity between a plurality of fingerprint functions from the set of fingerprint functions , utilizing equation ( 2 ) as described in u . s . patent application publication no . 2012 / 0330632 , reproduced below : where fp i is a fingerprint vector of crystal structure of a i , ∥ fp i ∥ is a norm of the fp i fingerprint vector , fp j is a fingerprint vector of crystal structure of b j , and ∥ fp j ∥ is a norm of a fp j fingerprint vector . as discussed above , the fingerprint vector is a result of discretization of the fingerprint function calculated according to equation ( 1 ). utilizing these equations , where identical structures are determined , only those that have a best fitness function , i . e . structures having a lowest free energy , are selected to participate in creation of child crystal structures . since pairs of parent crystal structures cannot be fully identical , the diversity of the next population improves by employing the niching algorithm . an optimal group of crystal structures is selected . to improve diversity when selecting the optimal group of crystal structures , a first optimal crystal structure from the parent group of crystal structures is discarded when a distance between the fingerprint function of the first optimal crystal structure and the fingerprint function of a second optimal crystal structure is smaller than a user - defined distance threshold . the user - defined threshold for the distance between the fingerprint functions of the optimal group of crystal structures is input as input parameter tolerancebesthm . however , the user - defined thresholds for the optimal group of crystal structures and the parent group of crystal structures may overlap . output of the niching algorithm provides an optimal group of crystal structures that survive unchanged into a next generation . the optimal group of crystal structures is selected as the first selected group , which represents crystal structures having most desirable crystal structures , i . e . structures with lowest free energy , from the optimized first generation of crystal structures . output of the niching algorithm also provides a parent group of crystal structures having crystal structures with a best fitness function , i . e . the lowest free energy . crystal structures with a poor fitness function are discarded . crystal structures that are duplicates of crystal structures within the optimal group and the parent group of crystal structures are preferably also discarded . the evolutionary algorithm may be tuned using input parameters to allow a user to optimize desired properties , including free energy , density , hardness , dielectric constant and bandgap , of the chemical structure . these input parameters include sigmafing , deltafing , rmaxfing , tolerancefing , tolerancebesthm , and goodbonds . approximately sixty to seventy optional input parameters may be employed depending on the properties of the chemical structure that the user wishes to tune . the parent group of crystal structures thus selected is output from the optimized first generation of crystal structures according to the fitness function . the parent group of crystal structures selected from the optimized first generation of crystal structures are chosen as a user defined function of the optimal group of crystal structures . thereafter , a child group of crystal structures are produced from the parent group of crystal structures using a variation operator that includes one of a lattice mutation operator , a coordinate mutation operator , a spatial heredity operator and an atomic permutation operator , to produce a child group of crystal structures from the parent group of crystal structures according to user - defined proportions . for example , 25 % of the parent group of crystal structures will be input into each of the lattice mutation operator , the coordinate mutation operator , the spatial heredity operator , and the atomic permutation operator to produce the child group of crystal structures . random variation operators provide a diverse population by exploring an entire search space . randomness ensures an unbiased search and limits genetic drift . however , completely random variation operators are typically inefficient when applied to crystal structures with a unit cell size of over forty atoms . to improve the variation operators , a local order element is introduced according to equation ( 3 ) in u . s . patent application publication no . 2012 / 0330632 , reproduced below : where i =[ 1 , . . . , n ] denotes individual atoms in a unit cell , n is a total number of atoms in the unit cell , v is a unit cell volume , fp a , b is a fingerprint vector of a i and b represents all atom types present in the unit cell . π i is the local order element for atom i and π i characterizes environment quality and atomic position symmetry , with a smaller value of π i for a given atom indicating a larger probability that the given atom will mutate . the local order element quantifies a degree of symmetry of the environment of a given atom and also quantifies the average degree of order of a fragment of a crystal structure . thus , the less distorted and less defective regions of the crystal structure may be input into the spatial heredity operator and more defective regions of the crystal structures may be input into the lattice mutation operator and the atomic permutation operator . when the local order element is low , a number of defects and distortion are high . where the local order element is high , the number of defects and distortion are low . however , a value of the local order element depends on factors that include atoms involved and atomic positions . if distortions or defects are favorable , the user may select those regions of the crystal structures accordingly . also a plurality of slices are selected from a parent group of crystal structures in a spatial heredity operator as a function of an order - fitness function correlation coefficient . each slice is a portion of a crystal structure formed when the crystal structure is cut along two parallel planes . all atoms located between these planes form the slice . specifically , two parallel cuts are made through the crystal structures and the portion of crystal structure between the parallel cuts is the slice . the spatial heredity operator preferably includes generating an amount of slices of a unit cell and selects a slice from slices having a highest average order . specifically , in the spatial heredity operator , the slices of the unit cell are selected from a plurality of crystal structures from the parent group of crystal structures and the slices are then combined to produce the child group of crystal structures . the slices selected from the parent group of crystal structures in the spatial heredity operator function as an order - fitness correlation coefficient . to improve diversity of crystal structures , a user may incorporate portions of both the less distorted and defective regions , as well as the more distorted and defective regions . to determine which portions to incorporate , correlation between the average order of the crystal structures and the fitness value of the crystal structures is calculated . depending on the degree of correlation , the degree to which the local order element influences the variation operators is determined . when there is a positive correlation , n , attempts to cut the slices are made from each parent and the most ordered slices are selected . that is , the number of slices of the unit cell from which the slice with the highest average order is selected is determined according to equation ( 4 ) in u . s . patent application publication no . 2012 / 0330632 , reproduced below : where n s is the number of slices , r ∈[− 1 ; 1 ] is a pearson product - moment correlation coefficient , l is a thickness of a unit cell in a direction of the slice for the spatial heredity operator , and l char is a length of half of a cubic root of an average atomic volume , with n s rounded to the nearest integer . thus , n , changes from a value of one in the absence of correlation ( r = 0 ) to l / l char for perfect correlation since the user does not need an exact and computationally expensive ab initio dynamical matrix to calculate phonon modes , the dynamical matrix may be calculated from bond hardness coefficients while setting atomic masses to unity . the soft - mode mutation operator includes calculating the dynamical matrix according to equation ( 5 ) in u . s . patent application publication no . 2012 / 0330632 , reproduced below : with coefficients α , β denoting coordinates ( x , y , z ); coefficients a , b , i , j describing an atom in a unit cell ; and coefficients l , m , n describing a unit cell number . r i , j l , n is a distance between atom i in the unit cell l and atom j in the unit cell n , while r 0i , j l , n the corresponding bond distance . h i , j l , n is a bond hardness coefficient between the atom i in the unit cell / and atom j in the unit cell n . for example , atoms are preferably treated as point particles and bonds as springs that connect the atoms with a stiffness determined by h i , j l , n . bonds with a valence strength greater than input parameter goodbonds are included in h i , j l , n . eigenvectors of d αβ ( a , b ) correspond to lowest non - zero eigenvalues and determine a direction of the soft mode mutation operator . valence strength is a chemical property of the bond between two atoms . a higher valence strength indicates a stronger bond between atoms . the input parameter goodbonds sets a threshold that determines which bonds are considered sufficiently strong to participate in the dynamical matrix . the dynamical matrix , as explained above , is used to calculate the eigenvectors and eigenvalues . the eigenvalue together with corresponding eigenvector is referred to as a mode . a mode having a relatively small eigenvalue is considered a soft mode in which mutation consists of movement of atoms along the eigenvector . the child group of crystal structures that are obtained are added to the optimal group of crystal structures and a next generation is output . following output of the next generation of crystal structures , the above - described steps are repeated a predetermined number of times for a predetermined number of generations , to output an optimized crystal structure , preferably in vienna ab initio simulation package ( vasp ) format describing the crystal structure in terms of unit cell distances to provide dimensions for each atom in the unit cell in each direction and providing coordinates of each atom within the unit cell . in embodiments , a method for predicting an optimized surface structure following the teachings of u . s . patent application publication no . 2012 / 0330632 includes inputting characterization information of a chemical structure and input parameters , with the characterization information including a number of atoms , a surface energy and a volume of a unit cell . a first convex hull of a first generation of surface structures is generated from the characterization information using a docking algorithm . surface structures of the first generation of surface structures include a bulk substrate layer , an adatom layer and a vacuum layer . the first generation of surface structures is restricted based on a convex hull algorithm according to a user - defined chemical potential constraint . the docking algorithm includes scanning the adatom layer for vacant sites and adding atoms to the vacant sites within the adatom layer based on predetermined conditions that include a coordination number and a distance from adjacent atoms . the docking algorithm includes outputting a configuration of the atoms added to the vacant sites . the method includes outputting a second convex hull of the first generation of surface structures and optimizing the second convex hull of the first generation of surface structures according to the input parameters . the optimized second convex hull is input into a niching algorithm to select an optimal group of surface structures and a parent group of surface structures based on a fitness function . a child group of surface structures is produced from the parent group of surface structures by applying a variation operator to the adatom layer of the surface structures of the first generation of surface structures . the child group of surface structures is added to the optimal group of surface structures to form a second generation of surface structures . the optimizing through adding steps are repeated for a predetermined number of generations , surface structures having a chemical potential above a user - defined threshold are discarded . in embodiments of the present disclosure , utilizing the above methods , a series of compounds of the general formula csf n ( n = 2 , 3 , 5 ) are provided which are stable at zero temperature . these compounds are also stable at zero pressure . for example , using csf 5 as a model , one mole of such compound ( 228 grams , occupying the volume of about 0 . 1 liter ) contains 2 moles of f 2 gas which , in the free state , would occupy a volume of about 44 . 8 liters . thus , storage in the form of csf 5 is three orders of magnitude more efficient , and much safer , than in the form of pure fluorine gas . fig1 is a convex hull diagram for the cs - f system at zero pressure and temperature , showing the formation of numerous ( meta ) stable fluorine - rich compounds , including csf 2 , csf 3 , csf 5 . the predicted materials are thermodynamically stable at ambient pressure and temperatures up to their respective decomposition temperatures , from about 250 ° k for the most fluorine - rich compound , csf 5 , to about 400 ° k for the least fluorine - rich compound , csf 2 , and their stability increases under pressure . one may be able to synthesize the compounds from the gas phase at normal pressure , but a more desirable method of synthesis is by the chemical reaction of csf with excess fluorine at normal or elevated pressure . is thermodynamically unfavorable at zero temperature ( the formation energy is approximately 0 . 18 ev / atom for this chemical reaction ), but will be favorable at elevated temperatures , due to the higher entropy over the f 2 molecules . for example , a procedure for the synthesis of csf 2 may include the following . csf 2 is stable at ambient pressure at temperatures below 400 ° k , and the maximum temperature of its stability increases with pressure . 1 mole ( 6 . 023 * 1023 molecules ) of f 2 weighs 38 grams . 1 mole of csf weighs 151 . 9 grams . to obtain 1 mole of csf 2 ( 170 . 9 grams ), one should mix 1 mole of csf ( 151 . 9 grams ) and 0 . 5 moles of f 2 ( 19 grams ). of course , one can proportionally increase the amounts of reactants to obtain a correspondingly larger amount of product . in reality , yield y is always less than 100 %, and can be modelled as : where a and b are constants ( of dimensionality k and s , respectively ), t is time and t is temperature . in a realistic case , at t = 400 ° k and t = 1000 seconds , y = 0 . 5 . increasing time and / or temperature will increase the yield . in embodiments , the synthesis is conducted at elevated pressures , in embodiments about 2 gigapascal ( gpa ), to enable higher temperatures for stable production , in embodiments about 600 ° k , and for longer periods of time , in embodiments about 3000 seconds . this way , one can expect the yield to increase to from about 80 % to about 99 %. analogous considerations apply also to kcl 3 . kcl 3 structures may be developed using the methods of u . s . patent application publication no . 2012 / 0330632 . optimal pressure - temperature conditions of synthesis of kcl 3 are similar to the formation of csf 2 described above , e . g . at a pressure of about 2 gpa and a temperature of about 600 ° k . to synthesize 1 mole of kcl 3 ( about 145 . 44 grams ), one can react 1 mole of kcl ( about 74 . 54 grams ) and 1 mole of cl 2 ( about 70 . 90 grams ). here also , one can proportionally increase the amounts of reactants to obtain a correspondingly larger amount of product . the general reaction scheme is as follows : the volumetric density of fluorine or chlorine in the compounds of the present disclosure is extremely high , 3 orders of magnitude higher than for a pure fluorine or chlorine gas at atmospheric pressure , with the compounds possessing reversible absorbed fluorine or chlorine , thus providing favorable advantages for the efficient and safe storage of fluorine and / or chlorine . compounds of the present disclosure thus prepared may be stored using any conventional means , including , for example , simple refrigeration . the present disclosure provides compounds that are very stable and effective storage materials for fluorine and / or chlorine . compounds of the present disclosure , including csf 2 , csf 3 , and csf 5 , contain polyfluoride anions , and fluorine can be easily recovered from these compounds in the form of f 2 gas by moderate heating . in embodiments , fluorine and / or chlorine may be recovered from the compounds of the present disclosure by moderate heating to temperatures from about 250 ° k to about 400 ° k , in embodiments from about 275 ° k to about 375 ° k , and ambient pressure . the cesium ( for the csf n compounds ), and the potassium ( for the kcl compounds ) may be recovered ( as csf for the csf n compound ) and then re - utilized by combining with fluorine and / or chlorine to form additional compounds of the present disclosure . the recyclability of these materials make the compounds of the present disclosure both very efficient for their use as storage materials as well as environmentally friendly . being very active , fluorine and chlorine easily form stable chemical compounds , from which they can be safely , easily and reversibly extracted . as f and cl have a huge range of industrial applications , the present disclosure provides great benefits for the storage of these materials . the compounds herein can reach high storage capacity , stability and reversibility . while the above description contains many specific details of methods in accordance with this disclosure , these specific details should not be construed as limitations on the scope of the disclosure , but merely as exemplifications of preferred embodiments thereof . those skilled in the art will envision many other possible variations that are all within the scope and spirit of the disclosure .