Microscopic Equations of State , Thermodynamical Properties , and System Collapse : Application of the Dynamical Equations of Period Vectors in Crystal Structure Prediction under External Temperature and Stress / Pressure ( Version 3 )

In crystal structure prediction, Newton’s second law can always be applied to particles (atoms or ions) in a cell to determine their positions. However the values of the period vectors (crystal cell edge vectors) should be determined as well. Here we applied the dynamical equations of period vectors derived recently based on Newton’s laws (doi:/10.1139/cjp-2014-0518) for that purpose, where the period vectors are driven by the imbalance between the internal and external pressures/stresses. For equilibrium states, they became equations of state, which essentially turn out the equilibrium conditions of crystals from mechanical point of view. Additionally for external pressure, they became Mie-Gruneisen equation supposing the Gruneisen constant is 1/3, which means phonon frequency is inversely proportional to the cell length. Since the internal stress has both a full kinetic energy term and a full interaction term, the influences of both external temperature and stress/pressure on crystal structures can be calculated, then thermodynamical properties and processes were presented. Contrary to usual ideas, the equations show that pure harmonic vibrating phonons can result in thermal expansion in crystals when the external temperature is changed. Finally, crystal system collapse due to temperature and/or stress/pressure change was discussed.


I. INTRODUCTION
Crystal structure prediction from being questioned becomes more and more feasible and important these years [1][2][3][4][5].Actually minimizing the enthalpy H = E + P V is normally employed to determine crystal structures for constant external pressure P , where E and V are the total potential energy and volume of a crystal cell respectively.Then, one has where n is the total number of particles (atoms or ions) in a cell, r i is the position vector of particle i, and h is any one of the period vectors (cell edge vectors) a, b, and c.Eq.( 1) means where F i is the net force acting on particle i. Eq.( 2) means where σ h = ∂V /∂h is the cell surface vector with respect to h.Actually the way in our previous paper [6] is equivalent to this.However the temperature effect is not shown there.
In this paper, we will apply the dynamical equation of period vectors derived recently [18] in crystal structure prediction, with temperature including phonon effect applied, and external stress considered.

II. DYNAMICAL EQUATION OF PERIOD VECTORS IN CRYSTAL STRUCTURE PREDICTION
In our recent paper [18], particles are always obeying Newton's second law (5) When the system reaches an equilibrium state, Eq.( 5) becomes Eq.(3).For determing the period vectors, the dynamical equation of them were derived based on Newton's laws as [18] where α h,h is some mass, tensor π is the internal stress, and tensor Υ is the external stress.As a matter of fact, crystal structure normally means the values of particle positon vectors and the period vectors in an equilibrium state, then let us regard all the so-called degrees of freedom (DOF) of the system in the paper [18], which are r 1 , r 2 , • • • , r n , a, b, and c, as those in the equilibrium state.Since the purpose of pursuing a crystal structure is to determine their final values from some guessed ones, the velocities of them can be completely forgotten.Then let us ignore the left side of Eq.( 6), but only change the period vectors in the directions of the calculated result of the right side of Eq.( 6).

III. EQUATIONS OF CRYSTAL STRUCTURE
The internal stress in Eq.( 6) has two terms The first term where E k,M D is the total kinetic-energy of the particles in the center cell and I is an identity tensor, then this term should be temperature dependent.Although no kinetic-energy for the DOF variables, which are the crystal framework in equilibrium states, other kinetic-energy do exist in the system separately.The kinetic-energy of the particles can be separated into two additive terms from ordered motion and disordered motion respectively.The kinetic-energy of the particles in the ordered motion can be further separated into two additive terms from the DOF variables (but zero here) and from vibrations around the DOF variables.The disordered one means thermal motion.Additionally kinetic-energy of free electrons should also be considered, while potential energy associated with electrons is always included in the total potential energy in any form.Then Eq.( 8) can be further written as where τ vb , τ th , and τ el are contribued from the kineticenergies of particle's vibrations (ordered), particle's thermal motion (disordered), and free electron's motion respectively.Specifically, where E k,M D,vb and E k,M D,el are the kinetic energies of the particle's vibrations and free electron's motion per cell.For the particle's thermal motion, it could be written as where k is Boltzmann constant and T is the temperature.However any of these tensor items may be ignored if the corresponding actual physics process is negligible.
The second term of the internal stress is where E p,M D , a function of the DOF, is the total potential energy of a cell, which should be the same as E in Eq.( 1).Explicitly, it is Considering Eq.( 3) (equilibrium states), it becomes Since in equilibrium states ḧ = 0 and using h • σ h = δ h,h V , where h = a or b or c, Eq.( 6) is Recalling the derivation of Eq.( 6) in our paper [18], Eq.( 16) actually means that if the crystal is cut into two halves, the net of all forces on each half is zero.Then the combination of Eq.( 3) and Eq.( 16) determine the crystal structure for given external temperature and stress.
For equilibrium states, the internal temperature should be the same as the external temperature.As the period vectors determine the cell volume, they are actually also equations of state of crystals, but in a microscopic form.
For the special case of external pressure P , Υ = −P I (real number P shoule be positive for compressing and negative for stretching), Eq.( 16) reduces to Comparing Eq.( 4) and Eq.( 17), we have an extra (kinetic energy) term τ • σ h which reflects the influence of temperature on crystal structures.

IV. THERMODYNAMICS
Since crystal structures can be determined based on the microscopic equations of state (Eq.(3) and Eq.( 16)), let us explore their thermodynamical properties and processes in the following three cases.

A. Isovolumic
Now let us consider the temperature is changed from T to T but keep the volume fixed.Actually, the fixed volume means the crystal structure does not change, then Eq.( 3) is always satisfied, and the total potential energy is also fixed.However, in order to satisfy Eq.( 16), since ∂E/∂h does not change, the externally applied stress must change accordingly from Υ to Υ as the temperature dependent (kinetic-energy) part of the internal stress changes from τ to τ The heat capacity at constant volume is then When the temperature is changed from T to T , even with the external stress Υ kept fixed, the crystal structure must change accordingly in order to satisfy Eq.( 3) and Eq.( 16).As a result, the volume should change from V to V , for an example.Then the coefficient of isobaric thermal expansion can be quickly calculated as The heat capacity at constant external pressure can also be calucalated the same way.
An alternative more accurate way to calculate them is to take derivatives of Eq.( 3) and Eq.( 16) with respect to temperature, since the equations should always be satisfied even if the external conditions are changed: Please note that all these are vector equations, then they are total 3(n + 3) scalar equations.In Cardesian coordinates, for every coordinate component k = x, y, z, the component F i,k of the force F i is a function of DOF, then As all the partial derivatives ∂F i,k /∂z are supposed known, all those total 3n equations are linear ones about the derivatives of DOF with respect to temperature (dr 1 /dT , dr 2 /dT , • • •, dr n /dT , da/dT , db/dT , and dc/dT ).Very similarly, are all linear combinations of the derivatives of DOF with respect to temperature.For example, h = a, considering constant stress Υ, then also a linear combination of the derivatives of DOF with respect to temperature.Assuming d(τ • σ h )/dT can also be expanded as a linear combination of the derivatives of DOF with respect to temperature, we will get a system of 3(n + 3) linear equations about the derivatives of DOF with respect to temperature.Once solved, the coefficients of isobaric thermal expansion can be calculated as The heat capacity at constant external pressure can be calculated the same way.
If the system expansion is assumed/proven evenly to certain degrees, the number of independent equations of the linear system should be reduced significantly.
With no external force, an interesting observation is when a crystal is formed at zero absolute temperature.Then there must be some attractive forces between particles.However since under this situation Eq.( 16) becomes ∂E/∂h = 0, there must also be repulsive forces to balance the attractive ones.Actually repulsive forces are believed only occuring between close neighbour particles.Then when the temperature, which acts as repulsive effect, is increased, in order to satisfy Eq.( 16), more attractive forces are needed.If the particles would go closer, the repulsive forces between particles would be increased too much, then the particles have to go apart more.Then the repulsive and attractive forces between particles both become less, but the repulsive lose more.As a result the inside forces between particles show more attractive featrue as a total to balance the temperature increasement.So the cell volume becomes bigger with temperature.

C. Isothermal
In this case, an external pressure, rather than a stress, is applied to the system.Everything is similar to the isobaric case, but the derivatives shoule be taken of Eq.(3) and Eq.( 17) with respect to pressure: Considering dT /dP = 0 for isothermal condition and supposing the above equations also generate a system of 3(n + 3) linear equations but about the derivatives of DOF with respect to pressure, the coefficients of isothermal compressibility can be calculated as

V. THERMODYNAMICS BASED ON PHONONS
Since many studies show that phonons play a very important role in crystal thermodynamical processes [7][8][9][10][11][12][13] where h = h/2π and h is Planck's constant.Considering Cardesian coordinates and using r i ,k as k = x, y, z component of particle i position vector r i , let us define Then all G i ,k ,i,k form a 3n by 3n matrix.The squareroots of the eigenvalues of the matrix are the total 3n frequencies of the phonons [17].As the DOF is the equilibrium values, instantaneously the total kinetic energy of phonons per cell is their total energy Then Eq.( 16) becomes This equation (and Eq.( 6)) tells us that the period vectors should change if the temperature is changed even if the external stress is kept fixed, although only harmonic phonons are considered.This thermal expansion by pure harmonic vibrating phonons is contrary to some researcher's ideas.An alternative choice is to employ Debye model for phonon motion, as it was demonstrated so successful [7][8][9][10][11][12][13].This can be done by bringing the total energy of phonons, for example Eq. (32) of reference [9], as their total kinetic energy into Eq.(34).

VI. SYSTEM COLLAPSE
Careful inspection on Eq.( 16) (also based on Eq.( 6)), actually the temperature effect, the first term on the left side of Eq.( 16), is to try to "break" the crystals, i.e. playing a repulsive role inside the system, then particles inside the system should attract each other enough to overcome it in order to keep crystal stable.However, there is a limit to the overall attraction inside the crystal (the middle term of the left side of Eq.( 16)).For fixed external stress, if the temperature is increased so much that its repulsive effect is bigger than the internal attraction limit and the external stress compensation, the crystal will collapse, where Eq.( 16) can no longer be satisfied.On the other hand, for whatever temperature, if the external stress is stretching and bigger than the internal attraction limit, the system will also collapse.The higher temperature, the less stretching is needed to break the crystal system.

VII. SUMMARY
By employing Newton's second law on particles in a cell and applying the dynamical equation of period vectors, crystal structures can be predicted.Not only the action of the external stress/pressure is presented in these equations, but also that of the (external) temperature via the kinetic energy term of the internal stress.They are actually equations of state of crystals in a microscopic form, then thermodynamical properties and processes can be calculated as well.Contrary to usual ideas, the equations show that pure harmonic vibrating phonons can cause thermal expansion in crystals.Furthermore crystals are also shown to collapse when the equations can not be satisfied.