Documentation
Notes: here are described mainly technical details about the program. Although I used the premise that the Reader has the specialized knowledge, I have nevertheless tried to formulate things in a more comprehensible form in several places.
Table of contents
From the very first minute of development, the most fundamental aspect was that the application should be useable by anyone, for free, without any registration.
Calculations of phase equilibria are typically complicated, and need a variety of data. So the second aspect became that the program should ask only the absolutely necessary data from the user. It should be as easytouse as possible, so the user can focus on the essence of the task. There should be no need for bothering about data that are necessary for the calculations but not interesting for the user.
Therefore, the program was maked very selfpropelled. Lots of data were collected, and the app selects what is necessary for the actual task  be it drawing, process modeling, or whatever. A wide variety of materials and their combinations can be applied for a range of issues  whether in education or in industrial practice.
Another aspect of development was reliability. The data of the substances were collected from several sources, as far as possible, in order to eliminate erroneous data. Methods of calculations provide results with the accuracy expected by current professional needs. The application is useful for demonstrational or industrial uses, even though it is not perfect, and there is still room for further improvement.
This application deals basically only with the equilibrium and transition between the vapor and the liquid state, not the solid state. The reason of this is that the nonideality models used here are not (fully) suitable for handling the solid state. Additionally, solid state addition compounds give much more complicated phase diagrams, than liquid miscibility diagrams or the vaporliquid equilibrium curves. (In fact, the solid phase adducts cannot be predicted by the methods used by VLE  Calc .com.) Therefore solidliquid equilibrium (SLE as solubility) became an another project. (Its development is currently suspended).
At the same time, liquidliquid equilibrium (LLE) and vaporliquid equilibrium (VLE) are handled together, which has also an important reason, that the model works well only used them together (VLLE). In case of multicomponent systems, before doing any VLE calculation, the program tests the system for liquid stability (is there phase split or not), and if the liquid splits to multiple phases, the program determines the accurate liquid compositions, and only then begins to calculate VLE.
These VLLE calculations are currently used for two purposes: drawing various phase diagrams and modeling distillation.
Approximately 300 compounds can be found in the database, with their identifyers, synonymous names, formulas, basic physical properties, and UNIFAC functional groups describing the structures. The data were taken from several reference sources, in order to eliminate erroneous data.
The physical properties:
 molecular weight (g/mol)
 liquid density (g/cm^{3})
 melting point (°C)
 boiling point (°C) at 1 atm.
 critical temperature (K)
 critical pressure (bar)
 critical volume (cm^{3}/mol)
 heat capacity of liquid (J/mol*K)
 heat capacity of solid (J/mol*K)
 heat of vaporization (kJ/mol) *
 heat of fusion (kJ/mol)
 heat of sublimation (kJ/mol)
 Antoine constants
 Acentric factor
 MatiasCopeman parameters
*
For some substances either Antoine constans, nor heat of vaporization data could be found in the literature.
But when the normal boiling point and a vapor pressure data at another temperature
(generally room temperature) were available, the heat of vaporization were
calculated from these data using the ClausiusClapeyron equation (see
chapter 4.1).
Vapor pressure of a pure substance is calculated usually with the Antoineformula, which is an empirical equation with three materialspecific constants.
where:
 p_{sat.,i} is the vapor pressure (bar)
 T is the temperature
 A,B,C are the Antoine constants
For the majority of substances in the database of VLE  Calc .com these constants are stored. Compounds for which the Antoine constans are missing, vapor pressure is calculated by the ClausiusClapeyron equation:
where:
 p_{sat.,i} is the vapor pressure (bar)
 ΔH_{V,i} is the heat of vaporization (ethalpy of vaporization) (J/mol)
 R is the universal gas constant (8.314 J/mol⋅K)
 T is the temperature (Kelvin)
 T_{BP} is the normal boiling point (at 1 atmosphere) (Kelvin)
In case of a multicomponent system, and when the temperature and pressure are not too close to the critical state (the vapor phase can be considered as ideal gas), one can use Raoult's law, modified by the activity coefficient:
where:
 p_{mixture} is the vapor pressure of the mixture (bar)
 γ_{i}(gamma) is the activity coefficient of component i (see in chapter 4.4)
 p_{sat.,i} is the vapor pressure of pure component i (bar)
 x_{i} is the liquid phase mole fraction of component i
 y_{i} is the vapor phase mole fraction of component i
There are substances whose molecules form associates in the vapor phase, and therefore erroneous vapor composition results using the modified Raoult's law. A typical example is the pair of acetic acid and water: the modified Raoult's law alone results in a minimal boiling azeotrope, although in reality this system is zeotropic. In this case, the association must be taken into account even at low pressures.
The method of J. Gmehling (see chapter 7) was adopted, whose essence is considering the association as an equilibrium chemical reaction.
To this equilibrium reaction we assign an equilibrium constant that depends only on the temperature:
K_{D} = p_{i,D} / p^{2}_{i,M}
where p_{i,M} is the partial pressure of the monomers, p_{i,D} is the partial pressure of the dimers.
We apply the following simplifications:
 we only deal with the association of carboxylic acids
 only dimerization is assumed
 crossdimerization (association of different compounds) is not taken into account
 the same equilibrium constant is used for all substances containing carboxylic acid groups (no substancespecific equilibrium constants are used)
 (The liquid phase association is included in the activity coefficient.)
The calculation of vapor phase dimerization goes according to the following:
 At first the equilibrium constant of the dimerization is calculated at the given temperature T:
ln K_{D} = A_{D} + (B_{D} / T)
(A_{D} and B_{D} are constants.)

Then the partial pressures of the monomers and dimers have to be calculated:
p_{sat.,i,M} = (−1 + √1 + 4⋅K_{D}⋅p_{sat.,i}) / 2⋅K_{D}
p_{i,M} = x_{i} ⋅ γ_{i} ⋅ p_{sat.,i,M}
p_{i,D} = K_{D} ⋅ p^{2}_{i,M}
 Finally it should be considered that the dimer species counts as 1 molecule in the vapor pressure,
but it counts as 2 molecules in the vapor composition:
p_{i} = p_{i,M} + p_{i,D}
y_{i} = (p_{i,M} + 2⋅p_{i,D}) / p_{mixture}
Model calculations of multicomponent systems usually need some method to take into account the interaction of different type materials. This program uses the socalled activity coefficient. (This is a kind of measure of interaction strongness between different molecules in a liquid solution, the "extent of nonideality".)
VLE  Calc .com currently uses the UNIQUAC and the UNIFAC methods. The former is used for specific substances and works with interaction parameters derived from specific experimental measurements. The latter is a groupcontribution method, so it can be used  theoretically  for any molecule, without having any interaction parameters based on experiments. Very simplified, UNIFAC works in a way, that molecules are defined as a set of different functional groups, and the activity coefficient is calculated from the sum of the molecular interactions between each grouppairs (using rather complex formulas).
In case of some common compound pairs the predictive power (accuracy) of UNIFAC seemed to be quite poor, so for some compound pairs UNIQUAC interaction parameters have been defined and recorded in the VLECalc.com database, derived from experimental data sets. When a user deals with a multicomponent system, in which all componentpairs' UNIQUAC parameters are available, the program automatically uses them instead of UNIFAC.
Two versions of UNIFAC is used by this app: one version for equation of state (EOS) usage, another version if no EOS is used.
In principle UNIFAC enables combining substances arbitrarily, unfortunately in practice it is not so simple. Certain combinations are allowed, while others are not. It depends on whether the interaction parameters are available for all functional group pairs occurring in the applied compounds, or not. Because if not, the program rather disallow the selection of such substances at the same time, in order to prevent giving erroneous results due to missing group interaction parameters.
This type of blocking is implemented in a way, that before displaying the list of substances (for compound selection) the program checks all the compounds which can be paired with the already selected one(s) and which not. Those which cannot be paired do not appear in the compound selector list, so they cannot be selected.
At normal conditions (meaning: common solvents, from room temperature up to about 100°C, at normal or reduced pressure) Raoult's law (in chapter 4.2) is fully suitable for vapor pressure calculations (naturally modified with the activity coefficient). (The vapor phase can be considered as ideal gas.)
But when the system's temperature or/and pressure is somewhat close to the critical point (or the mixture contains supercritical compontents), the vapor phase cannot be considered as ideal gas anymore. In such case an equation of state (EOS) should be used, which can describe the properties of the mixture realistic, even near to the critical state.
VLE  Calc .com uses currently the PSRK (Predictive SoaveRedlichKwong) equation of state. It's advantage is that it uses the UNIFAC model which could be fitted excellent to the activity coefficient part (chapter 4.4).
Several specific data are needed for EOS calculations (e.g. the critical constants of the material and others). When opening the compound selector list, the program checks all substances that the EOSrelated data are available, or not. Because if not, the program either disallow the selection of such substances (by hiding them in the list), or disables the use of EOS and only enables only "normal" conditions (i.e. temperatures/pressures far from the critical state) for the VLE calculations.
Calculating the vapor pressure for a given temperature is relatively easy according to section 4.2 (Raoult'slaw). The solution of its reverse in case of a single material is also relatively simple (see: both the Antoine and the ClausiusClapeyron equations can be inverted to express ecplicitly the T). But for a multicomponent system it is not trivial how to get the boiling temperature T of p_{mix}.
Literature recommends in general iterative methods, just like the quasiNewton secant method, which have been used by early versions of VLECalc.com. Now, the App uses a method wich is based on the fact that the vapor pressure function of ln(p) − 1/T is quasilinear (even in case of a multicomponent system), altough not in a very wide temperature range. The algorithm:
 first we calcualte the boiling temparature of all of the pure components at the give pressure
 then we calculate the vapor pressure of the mixture at the lowest and at the highest temperature
 based on the 22 pressuretemperature data we now fit a linear function on ln(p) − 1/T (slope and intersect)
 based on the fitted linear function one could calculate the searched temperature for the given temperature  this is the bubble point of the mixture at the given pressure.
 For the third temperature we got, we calculate again the vapor pressure of the mixture. (This should be very close to the given pressure. If not, something were wrong along the calculations.)
 Now we may fit a cubic function (instead of a linear one) based on the calculated 33 ln(p) and 1/T values, and we get the bubble temperature of given pressure from this.
 (When vapor composition is necessary, the one should calculate one more time the mixture vapor pressure at the final bubble temperature, and the composition obtains as a byproduct of the results.)
Note: this result is relatively accurate. Error is about 0.1°C. For less demanding applications one can stop at this point. But this result can be made much more accurate at relatively low cost, with several orders of magitude.
This method is computationally cheap (thus fast), mixture vapor pressure should be calculated only fourtimes in order to get the most accurate result, while the secant method usually needs more.
Introduction
Various liquid materials are either miscible, or not. (But that is certain.) Well, it is wellknown that hydrocarbons and lots of other organic stuffs do not mix with water, but there are plenty of cases where it is not trivial at all whether the mixture splits to multiple liquid phases, or not.
In principle, liquid phase split can be predicted from the curvature of the Gibbs free energy of mixing function (g^{E} = ΔG_{mix.}(x,T)/RT): when the function is convex from below at a composition range, i.e. ∂^{2}g^{E}(x,T)/∂x^{2} < 0, then the mixture will split to multiple liquid phases. This range is the unstable zone (see as green in the picture).
In case of two components the composition of the phases can be predicted by aligning a tangent line according to the figure (see as red line). In equilibrium the activity (a_{i}=γ_{i}⋅x_{i}) of all components are equal in all phases, so a computationally equivalent solution is to find the compositions with isoactivity (a_{i,phase1} = a_{i,phase2}).
In practice, there is some trouble with this. One can notice in the figure that the g^{E} function above the red tangent line is not everywhere convex from below (see as blue in the picture). This is the metastable zone. If the gross composition of a mixture falls into the metastable zone, testing the sign of the second derivate will not show that the mixture is biphasic.
In case of three or even more components the situation becomes more complicated, because the concepts (zones) described above should be interpreted and computationally handled on two or more dimensional surfaces. Practice of this is very not trivial, several article were published in this topic in the recent decades. Former versions of VLE  Calc .com used to use methods based on second derivate tests and isoactivity searches described above, thus the algorithms could not detect heterogenity in all cases. Another critical problem was that these algorithms produced sometimes erroneous equilibrium concentrations, the tie lines on the LLEdiagrams did not lined up consistently, even crossing tielines occurred. (The latter was not due to a bug in the program code, but due to properties of the thermodynamic model  see later in this chapter).
Liquid stability test
Due to the problems described in the introduction, the method of Wasylkievich et al were adopted to current version of VLE  Calc .com. This one is a very sophisticated approach that seemes to be reliable  though it was extremely difficult to reproduce and implement here.
The essence of the method is that the local minima (or inflexion points) of the socalled TangentPlaneDistanceFunction (TPDF) have to be located, and the equilibral compositions of the liquid phases should be searched starting from the found local minima. The method has several advantage:
 indicates heterogenity even in the metastable zone
 can be used for any number of components (in theory)
 it provides very good starting points for finding the equilibrium compositions
 indicates when the found equilibral composition are false
 indicates when the system splits into more than 2 phases (in case of more than 2 components)
At a given T temperature and z gross mixture composition:
TPDF(x) = g^{E}(x) − L(x, z)
where L(x, z) is a function tangent to g^{E} at z gross mixture composition.
Is TPDF(x, z) negative in any x, then the liquid mixture is unstable, and will split to multiple phases. (See with blue on the TPDF figure)
Twocomponent systems
In case of a binary system the local minima of the TPDF (as well as inflexions and local maximum) can be found at the zeros of the 1^{st} derivates. In current implementation these roots are searched by a quasiNewton method. Not all roots will be searched, but the two probably closest to the minima of the TPDF.
In this method we exploit some useful properties of the binary g^{E} function:
 there is no always local maximum, but if it exist, there is always only one.
 the slope of the function from the borders of the composition space into inward direction is always negative
Thus even if a local maximum exist, starting a search from the borders one would find always local minuma or inflexion point. The 1^{st} derivate of the TPDF can be expressed by an approximate analytical derivative, and this expression should be made equal to zero:
∂TPDF(x_{1,2}) / ∂x_{1} ≅ (ln x_{1} + ln γ_{1}(x) − ln z_{1} − ln γ_{1}(z)) − (ln x_{2} + ln γ_{2}(x) − ln z_{2} − ln γ_{2}(z)) = 0
The algorithm:
 Compute the activity coefficients at gross composition z_{}, and store them (they are constants during the processγ_{1,2}(z)).
 Start "guessed x"with point x_{1}=0.
 Compute the activity coefficients at current x_{} (γ_{1,2}(x)).

Using the analytical derivate described above one can express the formulae below, from which the new x can be obtained:
c_{1} = ln γ_{1}(x) − ln z_{1} − ln γ_{1}(z)
c_{2} = ln γ_{2}(x) − ln z_{2} − ln γ_{2}(z)
x_{1,NEW} = exp (−c_{1} + ln x_{2} + c_{2})
x_{2,NEW} = 1 − x_{1,NEW}
 Calculate the activity coefficients again with the new x value, and iterate the process until convergence. The resultant x is the closest minimum or inflexion to x_{1}=0.
 Repeat the process described above starting from x_{1}=1. From this process can be obtained the second local minimum or inflexion  if it exist. Namely, if there is only one minimum and no inflexion then both process gets to the same result, hence the system is homogenous.
Check the sign of the TPDF at the resultant x values. If it is negative, the system is unstable (heterogenous, biphasic). Now proceed searching for equilibral compositions of the phases, using these x values as starting points (this process is the "LLEflash"). In general the RachfordRice equation is used for that (see below).
Systems with three or more components
As mentioned in the introduction, the situation in case of a multicomponent system is much more complicated as in binary case. The method of S. Wasylkievich et al handles the situation with high reliability, in case of any number of components. However the process is computationally quite intensive. I would not go into the formulas here anymore, everything is described in detail in the cited article.
The essence of the method is that in the case of n component the a TPDF is a (n−1)dimensional surface, and the local minima of it has to be find. This surface is "crumpled cloth"like, it has local minima (full red dots on the figure), local maximum can occur on it, and there can be saddle points (empty red dots on the figure). These points are called the stationary points of TPDF. In these points all the partial derivatives of the TPDF, by mole fractions of all i components, are zero ( ∂ TPDF(x) / ∂x_{i} = 0). These stationary points are connected by "ridges" and "valleys", and according to the method, theese points are searched by tracking these ridges and valleys, along the whole (n−1)dimensional composition space.
Once the minima are found, similarly to the binary case, the sign of TPDF at the minimums are checked, whether it is negative or positive. If any of the minima is negative, the system is unstable (heterogenous). Then proceed searching for the equilibral compositions of the phases, using the x values of the local minima as starting points (this process is the "LLEflash"), and with solving the RachfordRice equation is used for that (see below).
Unfortunately, in multicomponent case, the first obtained phase compositions are not surely correct  it has to be checked. It can happen that not only one, but more isoactivity pointpairs can be found for a given z gross composition. From these solutions, the one with the lowest total Gibbs free energy of the system is correct.
So the obtained phase compositions have to be tested for stability, according to the ridge tracking technique described above (let z = x_{oneofthephases}). If x_{oneofthephases} is stable (i.e. no negative local minimum found), then the found equilibral phase compositions are correct. If the stability test results unstable (negative local minimum are found), then the found equilibral phase compositions are incorrect, and the LLEflash process has to be repeated using the mole fractions of the local minima obtained in the last stability test, as initial values.
From the number of local minima, it can be assumed how many liquid phases are, but there is no guarantee that it is actually as much. E.g. if there are three minima, it is worth assuming this, and then look for three equilibrium phases based on the minima. If there were no three phases, it would be revealed either in the RachfordRice calculations (the system automatically reduces to 2phase) or in the stability test of the found phase compositions.
The RachfordRice equation
Once it is known that the system is heterogenous, it is time to calculate the accurate composition of the coexisting luqid phases. In equilibrium the activity of the components (a_{i}=γ_{i}⋅x_{i}) is equal in all phases. At the same time it has to be also take into account that the material balance of each component must be in order, i.e. the total amount of each component in the separate phases must be equal to that of in the gross mixture (z). Rachford and Rice derived a very useful little equation for this task, solving it leads to results that meet both expectations at the same time:
 F(α) is the objective function, which has to be equal to zero
 α is a number between 0 and 1, which measures the molar amount of the 1^{st} separate phase, from 1 mole gross mixture. (The molar amount of the 2^{nd} phase is logically 1−α.)
 z_{i} is the mole fraction of component i in the gross mixture
 K_{i} is the distribution ratio of component i: K_{i} = x_{i,phase 1}/x_{i,phase 2} = γ_{i,phase 2}/γ_{i,phase 1}
The RachfordRice equation can be solved by iterative way:
 In the first step calculate the distribution ratios K_{i} using the assumed phase compositions and activity coefficients originating from the liquid stability test.
 Determine α by the NewtonRaphson method. Substitue everything into F(α) and into it's first derivate according to α, and the next approximation of α can be obtained by the formula: α_{new} = α_{previous} − F(α_{previous}) / F'(α_{previous}). Repeat with the new α and iterate until convergence.
 Calculate the phase compositions with the formulae:
x_{i,phase 1} = z_{i} / (1 + α(K_{i} − 1))
x_{i,phase 2} = x_{i,phase 1} ⋅ K_{i}  Normalize the obtained mole fractions (for all phases): x_{i} = x_{i} / ∑^{n}_{j=1} x_{j}
 Calculate the activity coefficients using the new phase compositions, and calculate again all K_{i}, and iterate the whole process until convergence.
Boiling point / vapor pressure map by contour lines
Although it is not the most useful representation, but this type is the oldest triangle diagram in this app, since it is relatively simple and can be obtained quickly. Basically it represents the boiling point (bubble point, in isobaric case) or the vapor pressure (in isothermal case) of a ternary mixture, using contour lines. It's disadvantage is that the vapor composition (as additional information) does not appear in any form in this mode of representation.
At first, the triangular composition space is divided into many small triangles with a resolution of 20x20, and a VLE calculation is performed at each corner point of each small triangle. The minimum and maximum temperature or vapor pressure values are determined from all of the VLE results, and 8 values are selected between them, for the 8 contour lines. These are the "constant values", since these values are constant along the contour lines. With each constant value, all the small triangles are checked wether this value appears in it, or not. If a constant falls between the lowest and highest temperature/pressure values at the corners of the small triangle, then the contour line of the constant value passes through the given small triangle. If so, the position of the crossing section in the triangle is calculated by linear interpolation, and stored.
Plotting all of the crossing sections the contour lines outlined.
Residue curve map
For distillation problems, the residue curve is a much more useful representation than the boiling point map. The residue curve lies along the path as the composition of the bottom product changes during a batch distillation process. Logically, it moves from the more volatile regions to the less volatile regions, and the direction of the curve at each point is determined by the equilibrium x − y vector, characteristic of that point.
The shape of residue curve diagrams can be very diverse. In fact, the azeotropic points of a given ternary system, basically determine the shape of the figure. The azeotropic points (together with the corner points of the entire triangular diagram, which represent the pure components) divide the total composition space into smaller areas. Simple batch distillation processes can operate generally in these areas, in which the starting composition falls, so with the help of the residue curve set, it can be easily determined visually whether a given separation task is possible or not.
At first the azeotropic points of the system are located, then the corner points of the whole triangular diagram are added to the set. These points are called nodes. The set of nodes is arranged in triangles by Delaunay triangulation. Within each small triangle, residue curves are generated starting from one or more points.
During generation of a residue curve, the area of the small triangle is scanned from the starting point first in the x → y direction and then in the x ← y direction until a node is reached. Then another starting point is chosen and scan again.
Liquid miscibility phasediagram
Liquidliquid equilibrium diagrams have numerous types, these are not discussed here just the method of figure making.
At first the algorithm tests for liquid stability the binary subsystems representing the edges of the triangle diagram. The phase compositions of the miscibility gaps of the heterogenous binaries are stored for drawing on the one hand, on the other hand they are used for define starting points in the middle of the gaps.
Then it starts stepping from the starting points in small steps towards the inside of the triangle and determines the ternary tielines or conodes. The direction is determined by the slope of the tieline (the direction is always perpendicular to the last tieline), and the step size depends on the length of the tieline  the shorter the tieline, the smaller the step size. (These rules are basically for aesthetic purposes.)
If the algorithm finds a threephase region, it continues to stepping from the middle of the other two edges of the triangle formed by the compositions of the three phases. Finally, plots all the found tielines.
In case of ternary vaporliquid equilibrium calculations, the liquidliquid equilibrium tielines are precalculated, stored, and arranged in small and thin triangles. Before starting the VLE calculation, the program checks if the given liquid composition falls into any of the LLE triangles. If so, it determines the composition of the separated liquid phases in a simplified procedure and performs VLE calculations with these compositions. (This precalculation and storetheLLEmap mess is necessary because the liquid stability test is computationally very intensive. Less stability testing is required for an LLE diagram than for a batch distillation.)
During a single stage continuous distillation (flash distillation) process the starting liquid mixture (feed stream) is evaporated partially. The composition of the head and bottom product depends on the evaporated proportion. In case of total evaporation (no bottom product) the composition of the distillate is the same as that of the starting mixture. In case of no evaporation (only bottom product), logically the bottom product is the same as the feed. The intermediate (real) situation is calculated using the RachfordRice equation (just like in case of liquidliquid equilibrium).
where
 F(α) is the objective function, which has to be equal to zero
 α is the evaporated proportion, a number between 0 and 1, the molar amount of the vapor phase from 1 mol feed mixture. (The molar amount of the bottom product
 z_{i} is the mole fraction of component i in the feed stream
 K_{i} is the distribution ratio of component i between the vapor and liquid phase: K_{i} = y_{i}/x_{i}
The RachfordRice equation can be solved by iterative way:
 The evaporated proportion α is given, constant.
 At first perform a bubble point calculation using the z feed composition.
 Calculate distribution ratios (as astarting approximation) using the resultant y values: K_{i} = y_{i}/z_{i}
 Calculate the composition of vapor and liquid phase:
x_{i} = z_{i} / (1 + α(K_{i} − 1))
y_{i} = x_{i} ⋅ K_{i} 
Normalize the mol fractions:
x_{i} = x_{i} / ∑^{n}_{j=1} x_{j}
y_{i} = y_{i} / ∑^{n}_{j=1} y_{j}  Perform another bubble point calculation using the resultant x values, and then recalculate the distribution ratios with the new y values (K_{i} = y_{i}/x_{i}), and iterate until convergence.
Batch distillation is a timely process, the composition of both the boiling liquid and the distillate changes over time. To model this, the process is divided into n equal parts. VLECalc.com always uses a resolution of n = 200, even if a lower resolution is given in the form, and only shows as much of the results in the table as specified. The program always performs the calculations in molar units.
The initial x composition of the mixture is given. The amount of the starting mixture is always n = 200 mol as described above (this made the code simpler), which is of course converted at the end according to the amount (and unit) specified by the user. A bubble point calculation is performed from x, and with the obtained vapor composition y one can calculate the amount of all components from the first evaporated dn = 1 mol:
n_{new} = n_{old} − dn
x_{i,new} = (x_{i,old} ⋅ n_{old} − y_{i} ⋅ dn) / n_{new}
Store the new boiling liquid composition and calculate from this again the vapor composition, and so on, iterate until the liquid runs out completely.
Physicochemical models operate usually with molar amounts and mole fractions, but in practice, man can usually think in mass and volume units.
For one type of task it is practical to measure the temperature in °C, for another task in Kelvin.
It proven to be advantageous to build in an easytouse unit converter.
In this chapter one can study the converting methods and formulas.
Basic function
The unit converter module handles currently nine types of quantities, however, some of them has  for now  only 1 unit type built in. (I.e. it cannot be converted into another unit.) Every quantity type has a default unit, in which the numeric values are stored and the calculations are made with. When showing the numeric values on the page, the program calculates always from the default unit to the user selected unit.
The convertible quantity types and the selectable units are as follows:
Quantity type  Default unit  Additional units 
Molar weigth  g/mol   
Liquid density  g/cm^{3}   
Temperature  K (Kelvin)  °C (Celsius), °F (Fahrenheit) 
Pressure  bar  kPa, atm, Hg mm, psi 
Molar volume  cm^{3}/mol  dm^{3}/mol, m^{3}/mol 
Specific heat  kJ/kg⋅K  J/mol⋅K, kCal/kg⋅K, Cal/mol⋅K 
Latent heat  kJ/mol  kJ/kg, kCal/mol, kCal/kg 
Concentration or composition  mole fraction (x)  mass%, volume%, mol/dm^{3}, g/100 cm^{3} 
Quantity of material  kg  g, t, dm^{3}, cm^{3}, m^{3}, mol, mmol, kmol 
Temperature conversion
[K]  =  [°C] + 273.16 
[K]  =  ([°F] + 459.688) / 1.8 
[°C]  =  [K] − 273.16 
[°C]  =  ([°F] − 32) / 1.8 
[°F]  =  1.8 ⋅ [K] − 459.688 
[°F]  =  1.8 ⋅ [°C] + 32 
Pressure conversion
[bar]  [kPa]  [atm]  [Hg mm]  [psi]  
[bar]  =  n.a  /100  ⋅1.01325  /750.0617  /14.5037738007 
[kPa]  =  ⋅ 100  n.a  ⋅101.325  /7.500617  /1450.37738007 
[atm]  =  /1.010325  /101.325  n.a  /760  /14.6959488 
[Hg mm]  =  ⋅ 750.0617  ⋅ 7.500617  ⋅ 760  n.a  ⋅ 51.71493367 
[psi]  =  ⋅ 14.5037738007  ⋅ 0.145037738007  ⋅ 14.6959488  /51.71493367  n.a 
Specific heat conversion
[kJ/kg⋅K]  [J/mol⋅K]  [kCal/kg⋅K]  [Cal/mol⋅K]  
[kJ/kg⋅K]  =  n.a  / M  / 4.1868  /(4.1868⋅M) 
[J/mol⋅K]  =  ⋅ M  n.a  ⋅ M /4.1868  /4.1868 
[kCal/kg⋅K]  =  ⋅ 4.1868  ⋅ 4.1868 / M  n.a  / M 
[Cal/mol⋅K]  =  ⋅ 4.1868⋅M  ⋅ 4.1868  ⋅ M  n.a 
where M is the mole weigth of the material.
Latent heat conversion
[kJ/mol]  [kJ/kg]  [kCal/mol]  [kCal/kg]  
[kJ/mol]  =  n.a  ⋅ 1000 ⋅ M  / 4.1868  ⋅ 1000 ⋅ M / 4.1868 
[kJ/kg]  =  ⋅ 1000 / M  n.a  /(4.1868⋅M)  /4.1868 
[kCal/mol]  =  ⋅ 4.1868  ⋅ 4.1868 ⋅ M  n.a  ⋅ M 
[kCal/kg]  =  ⋅ 4.1868 / (1000⋅M)  ⋅ 4.1868  / (1000⋅M)  n.a 
where M is the mole weigth of the material.
Concentration / composition conversion
Conversion between various concentration types is somewhat more difficult compared to the simpler ones above, because the calculations need some of the physical data of the components.
[mole fraction] (x)  [mass%] (w)  [vol.%] (v)  [mol/dm^{3}] (c)  [g/100 cm^{3}] (ρ)  
x_{i}  =  n.a  w_{i} M_{i} ⋅ ∑^{n}_{j=1} (w_{j} / M_{j}) 
v_{i}⋅d_{i} M_{i} ⋅ ∑^{n}_{j=1} (v_{j}⋅d_{j} / M_{j}) 
c_{i} ∑^{n}_{j=1} c_{j} 
ρ_{i} M_{i} ⋅ ∑^{n}_{j=1} (ρ_{j} / M_{j}) 
w_{i}  =  100% ⋅ x_{i} ⋅ M_{i} ∑^{n}_{j=1} (x_{j} ⋅ M_{j}) 
n.a  100% ⋅ v_{i} ⋅ d_{i} ∑^{n}_{j=1} (v_{j} ⋅ d_{j}) 
100% ⋅ c_{i} ⋅ M_{i} ∑^{n}_{j=1} (c_{j} ⋅ M_{j}) 
100% ⋅ ρ_{i} ∑^{n}_{j=1} ρ_{j} 
v_{i}  =  x_{i}⋅M_{i} d_{i} ⋅ ∑^{n}_{j=1} (x_{j}⋅M_{j} / d_{j}) 
w_{i} d_{i} ⋅ ∑^{n}_{j=1} (w_{j} / d_{j}) 
n.a  c_{i}⋅M_{i} d_{i} ⋅ ∑^{n}_{j=1} (c_{j}⋅M_{j} / d_{j}) 
ρ_{i} d_{i} ⋅ ∑^{n}_{j=1} (ρ_{j} / d_{j}) 
c_{i}  =  1000 ⋅ x_{i} ∑^{n}_{j=1} (x_{j} ⋅ M_{j} / d_{j}) 
1000 ⋅ w_{i} M_{i} ⋅ ∑^{n}_{j=1} (w_{j} / d_{j}) 
1000 ⋅ v_{i}⋅d_{i} M_{i} ⋅ 100% 
n.a  1000 ⋅ ρ_{i} M_{i} ⋅ 100 
ρ_{i}  =  100 ⋅ x_{i} ⋅ M_{i} ∑^{n}_{j=1} (x_{j} ⋅ M_{j} / d_{j}) 
100 ⋅ w_{i} ∑^{n}_{j=1} (w_{j} / d_{j}) 
100 ⋅ v_{i} ⋅ d_{i} 100% 
100 ⋅ c_{i} ⋅ M_{i} 1000 
n.a 
where M is the mole weigth, and d_{i} is the liquid density of component i.
Important notes: quantities per volume calculated using these formulas are not accurate, due to the socalled volume contraction phenomenon and to the temperature dependence of liquid density, which are not taken into account by these formulas. The order of magnitude of this inaccuracy may approx. a few rel.%.
New figures
The calculation methods of VLECalc.com allow to construct new useful diagrams, as well as to extend existing figures with additional information. Some of the plans:
 pT phase diagram of multicomponent mixtures ("phase envelope")
 pV phase diagram of multicomponent mixtures
 azeotropic composition in function of the pressure
 figure illustrating batch distillation processes (composition change of distillate or residue)
 ternary VLEdiagram extended with LLEdiagram drawn on it
New functions
Some upgrades for more convenient use and for usefulness in even more areas are also in the plans.
 mixtures with more than 3 components
 rectification (fractional distillation) calculation
 handling of distillate fractions (fore run, after run) at distillation calculation
 predicting physical properties of pure compounds using group contribution methods
 predicting physical properties of mixtures
 calculation of temperaturedependent physical properties at usergiven temperature
 putting more and more new compounds into the database
 ability to handle custom (usergiven) molecule
Known bugs and deficiencies
Onfortunatelly this new version of the app has some new bugs. Although, since the release of current version several errors have been found (thanks a lot for every feedback!), what could be fixed quickly was fixed quickly.
However, there is a stubborn, hidden bug  or malfunction  that could not be fixed yet. It occurs most often when dealing with ternary, heterogeneous systems; the program stops with an error message when attempting to do VLE calculations. The probability of this runtime error is higher, when very low concentrations are there among the LLEdata. The root cause of the malfunction has been largely determined, but no reassuring solution has been found so far.
Another annoying novelty is that the new LLEprocedure (the new liquid stability test) is computationally very intensive, thus slow even on today’s fast computers. (Although it is more reliable and knows more than the old one in return.) One of the goals is to change this as well.
Outlook
Several projects are in progress or are planned, which are more or less similar to VLECalc.com, but are to be used at different areas. An example of them is solubility.info which deals mainly with the solubility of solids in solvent mixtures. Another example may be azeotrope.info which moved recently to its own domain name.
An interesting and good news may be for those, who asked about a downloadable version of VLECalc.com (which does not exist), that an Excel version is under development. Most of the work is done already, product finalization, testing and debugging is in progress now. It is planned to publish in 2021.
Methods of calculations

A. Fredenslund, R.L. Jones, J.M. Prausnitz: Groupcontribution estimation of activity coefficients in nonideal liquid mixtures, AIChE Journal, Vol. 21, No. 6, 1975, page 10861999
(DOI: 10.1002/aic.690210607)

J. Gmehling, J. Li, M. Schiller: Modified UNIFAC model. 2. Present parameter matrix and results for different thermodynamic properties, Ind. Eng. Chem. Res. 1993, 32, 178193
(DOI: 10.1021/ie00013a024)

S. Horstmann, A. Jab³oniec, J. Krafczyk, K. Fischer, J. Gmehling: PSRK group contribution equation of state: comprehensive revision and extension IV, including critical constants and alfafunction parameters for 1000 components, Fluid Phase Equilibria 227 (2005) 157–164
(DOI: 10.1016/03783812(91)85038V)

Stanislaw K. Wasylkiewicz et al: Global Stability Analysis and calculation of LLE in multicomponent mixtures, Ind. Eng. Chem. Res. 1996, 35, 13951408
(DOI: 10.1021/ie950049r)
 J. Gmehling et al.: VaporLiquid Equilibrium Data Collection series, 1981 Dechema
Activity coefficient
Equation of state
Liquid stability and LLE calculation
Vaporphase dimerization
Data of substances
 E.W. Flick: Industrial Solvents Handbook, 5th ed., 1998 William Andrew Publishing/Noyes
(eBook ISBN: 9780815518099)
 NIST Chemistry Webbook
 Wikipedia
 R. H. Perry, Don W. Green: Perry's Chemical Engineers' Handbook, 1999 McGrawHill Inc.

Jr. W. Acree, J. S. Chickos: Phase Transition Enthalpy Measurements of Organic and Organometallic Compounds. Sublimation, Vaporization and Fusion Enthalpies From 1880 to 2010, J. Phys. Chem. Ref. Data, Vol. 39, No. 4, 2010
(DOI: 10.1063/1.3309507)
 I. Mellan: Source Book of Industrial solvents, vol.II. (Halogenated hydrocarbons), Reinhold Publishing Corp., N.Y., 1957
 ChemSpider
 SigmaAldrich catalog

B. E. Poling, J. M. Prausnitz, J. P. O'Connell: The Properties of gases and liquids, 5th ed., McGrawHill, 2001
(DOI: 10.1021/ja0048634)

P. Martins, P. Sbaite, C. Benites, M. Maciel: Thermal Characterization of Orange, Lemongrass, and Basil Essential Oils; School of Chemical Engineering, University of Campinas (UNICAMP), 13083970, Campinas –SP, Brazil
(DOI: 10.3303/CET1124078)
 PubChem

S. Horstmann, A. Jab³oniec, J. Krafczyk, K. Fischer, J. Gmehling: PSRK group contribution equation of state: comprehensive revision and extension IV, including critical constants and alfafunction parameters for 1000 components, Fluid Phase Equilibria 227 (2005) 157–164
(DOI: 10.1016/03783812(91)85038V)

N. N. Greenwood, A. Earnshaw: Chemistry of elements, Second edition, ButterworthHeinemann, © Reed Educational and Professional Publishing Ltd 1984, 1997
(eBook ISBN: 9780080501093)