The Born Oppenheimer Approximation In Quantum Chemistry Engineering Essay

Published: November 21, 2015 Words: 4205

In the assumption, the coupling between the nuclei and electronic motion is neglected. This allows the electronic part be solved with the nuclear positions as parameters, leading to potential energy surface (PES). It forms the basis for solving the nuclear motion.

Why? Because electrons' masses are small and they move fast, they cannot be described correctly even qualitatively by classical mechanics while nuclei are sufficiently heavy and move slow. Quantum effects like excitation energy are obvious on electrons and negligible on nuclei. Nuclei are stationary in the viewpoint of electrons, and electrons adjust instantly to the change of nuclear geometry.

This approximation leads to the following descriptions of molecular systems and chemical reactions:

Stable molecules correspond to minima energy on the potential energy surface within the Born-Oppenheimer approximation and a chemical reaction can be described as nuclei moving from one minimum to another. In the lowest degree of approximation, the motion is assumed to occur along the path of least energy, and this path forms the basis for transition state theory.

Write an expression for the energy of a Slater determinant (in terms of integrals over different parts of the Hamiltonian operator and using the spin-orbitals of the Slater determinant). Explain the different terms in this energy expression. (2p)

Eel = iIi + 1/2 i,j'(Jij - Kij)

Here, iIi is the kinetic energy of the electron plus the Coulomb attraction between all electrons and nuclei.

1/2 i,j' Jij is the Coulomb repulsion between the i and j electrons,

1/2 i,j' Kij is the exchange interaction energy, the i- and j-th electrons with no spin.

Rewrite the integrals as Eel = iIi + 1/2 i,j'(<Ψi(x1)Ψj(x2)|1/rij|Ψi(x1)Ψj(x2)> - <Ψi(x1)Ψj(x2)|1/rij|Ψj(x1)Ψi(x2)>).

Which are the approximations made in the Hartree-Fock method? (2p)

Born-Oppenheimer approximation;

one-particle functions are spin-orbitals;

Wave-function is a Slater determinant, i.e. Pauli Exclusion Principle;

Hamiltonian contains only electrostatic interaction, i.e. ignoring the correlation between electrons.

Derive the Hartree-Fock equations in canonical form. Formulas and explanations must be given for all non-trivial steps. (4p)

Assume that Φ0 is a Slater determinant built from the optimal orbitals,

E= < Φ0|H|Φ0>, < Φ0|Φ0> = 1.

A variation of the wave function

Φ0 → Φ0+δΦ0, the energy of Φ0 is at a minimum: δE = <Φ0+δΦ0|H|Φ0+δΦ0> - <Φ0|H|Φ0>=0.

= <Φ0|H|Φ0> + <Φ0|H|δΦ0> + <δΦ0|H|Φ0> + <δΦ0|H|δΦ0> - <Φ0|H|Φ0> ≈ 2 Re {<δΦ0|H|Φ0>} = 0

(Neglecting the second order term, < δΦ0|H|δΦ0>)

< Φ0+δΦ0|Φ0+δΦ0>-{<Φ0|Φ0>} = 0.

we can write the two criteria that the optimal orbitals must fulfill:

Re {<δΦ0|H|Φ0>}=0

Re {<δΦ0|Φ0>}=0.ï€½ï€ ï€°

Φ0 is made up of spin-orbitals: Ψ1, Ψ2, Ψ3, … Ψn (n electrons), taken from a complete set of orthogonal spin-orbitals:

Ψ1, Ψ2, Ψ3, … Ψn Ψn+1, Ψn+2, …

occupied virtual (non-occupied)

To form δΦ0, vary the spin-orbital: Ψi'= Ψi + δΨi (δΦ0: Ψi → δΨi)

Spin-orbitals are kept orthogonal => <δΨi|Ψj> = 0 (j=1, 2,…n, j≠i)

<δΦ0|Φ0> = <δΨi|Ψi> = 0 and the criteria 2 is fulfilled.

Expand δΨi in the complete set of orbitals: δΨi = Σl∞αilΨl

Orthogonality restriction: αil = 0 for l ≤ n

For simplicity, δΨi = αΨk k > n, i.e. Ψk is a virtual orbital.

Φ0+δΦ0 = | Ψ1,…Ψi-1,(Ψi+ αΨk),…Ψn| = | Ψ1,…Ψi-1, Ψi,…Ψn| + α| Ψ1,…Ψi-1,Ψk,…Ψn|

δΦ0 = αΦi→k singel excitation

Re{α<Φi→k|H|Φ0>} = 0 (the criteria 1)

Choose α so that {α<Φi→k|H|Φ0>} is real, <Φi→k|H|Φ0> = 0 (Brillion's theorem)

Single excitations do not interact with the HF determinant (orbitals are optimal).

H = Σin f(xi) + 1/2 Σi,j'1/rij

<Φi→k|H|Φ0>

= <Ψk| f |Ψi> + Σjn (<ΨkΨj|1/ rij|ΨiΨj> - <ΨkΨj|1/rij|ΨjΨi>) = 0

Define Fock operator: F = f + Σjn (Jj - Kj)

Jj and Kj are the Coulomb and Exchange operators:

<ΨkΨj|1/rij|ΨiΨj> = <Ψk|Jj|Ψi> <ΨkΨj|1/rij|ΨjΨi> = <Ψk|Kj|Ψi>

Then, <Φi→k|H|Φ0> = 0 can be written: <Ψk|F|Ψi> = 0 (i ≤ n, k > n)

Correspond to a block-diagonal Fock matrix and define orbitals with the lowest energy.

To find a set of unique orbitals, expand FΨi in the complete set of orbitals, FΨi = Σj∞εjiΨj

<Ψk|F|Ψi> = Σjεji <Ψk| Ψj> = εki <Ψk| Ψk> = εki

εji = <Ψj|F|Ψi> = 0 for j > n, i ≤ n, i.e. the ε matrix is block-diagonal

FΨi = ΣjnεjiΨj for i ≤ n

Make a unitary transformation of the orbitals, such that the ε matrix becomes diagonal.

FΨi = εiΨi for i = 1, 2,…n Hartree-Fock equations for the occupied orbitals

e. What is the difference between restricted and unrestricted Hartree-Fock, answer in terms of definitions of the methods and consequences for the calculations. (2p)

The Slater determinant can be written in terms of spin-orbitals, being combination of a spatial orbital and a spin function. If there are no restrictions on the form of the spatial orbitals, the function is an Unrestricted Hartree-Fock (UHF) wave function.

The term different orbitals for different spins (DODS) is also sometimes used. If we are interested in systems with an even number of electrons and a singlet type of wave function, we will make the restriction that each spatial orbital should have two electrons, with two kinds of spins. Such wave functions are Restricted Hartree-Fock (RHF). Restricted HF method forces the α-spin and β-spin cases have the same space orbitals, while Unrestricted Hartree-Fock method (UHF) allows different space orbitals for them.

In closed shell case (mα = mβ), space orbitals are the same, so, RHF and UHF give the same result. In open shell case (mα ≠mβ), the operators are different for α-spin and β-spin. RHF can get eigenfunctions to the spin operaters, only the open shells (unpaired) electrons contribute to the total spin. However, UHF calculated results are not eigenfunctions of the S2-operator.

Some aspects of quantum chemical calculations

What is the difference between STO and GTO basis sets? Which type is mostly used and why? What is meant by a contracted basis set? (2p)

The difference between the STO and GTO is in the "r" of their expressions. The GTO squares the "r" so that r2 dependence in the exponential makes the GTOs inferior to the STOs in two respects. At the nucleus, a GTO has a zero slope in contrast to a STO which has a discontinuous derivative, and GTOs consequently have problems representing the behaviour near the nucleus. The other problem is that the GTO falls off rapidly far from the nucleus compared with an STO, and the "tail" of the wave function is represented poorly. In Slater type orbitals, there are no radial nodes; and the nodes in the radial part are introduced by making linear combinations of STOs. STOs are usually used for atomic and diatomic systems where high accuracy is needed and in semi-empirical methods where all three- and four-centre integrals are neglected.

Both STOs and GTOs can be chosen to form a complete basis, but more GTOs are necessary for getting an accuracy compared with STOs. In terms of computational efficiency, GTOs are preferred and are used almost universally as basis functions in electronic structure calculations.

Contracted basis sets means the atomic orbitals are made up of a linear combination of primitive STOs or GTOs.

Describe in general terms Mulliken population analysis. What is meant by gross population and overlap population? What chemical concepts can these quantities be considered to describe? (2p)

If the coefficients of the basis functions in the molecular orbital are Cμi for the μ'th basis function in the i'th molecular orbital, the density matrix terms are Dμv=2iCμvC*vi for a closed shell system, each molecular orbital is double occupied. The population matrix P then has terms Pμv=(DS) μv, , S is the overlap matrix of the basis functions. The sum of all terms of Pμv is N -- the total number of electrons. Overlap populations Pμv between the corresponding orbital populations Pvv and Pμμ. The Mulliken population analysis aims first to divide N among all the basis functions. This is done by taking the diagonal element of and then dividing the off-diagonal elements equally between the two appropriate basis functions. Since the off-diagonal terms include Pμv and Pvμ, this simplifies to just the sum of a row. This defines the gross orbital population (GOP) as GOPμ=vPμv..

Mulliken charges come from the Mulliken population analysis and provide a means of estimating partial atomic charges from calculations carried out by the methods of computational chemistry.

Explain the problem with doing calculations on low-spin coupled open-shell states (antiferromagnetic coupling) using single determinant methods (DFT and HF). Give an example of how this can be reasonably well handled. Use a simple example with two unpaired electrons. (2p)

The low-spin coupled open-shell state, i.e. anti-ferromagnetic coupling, cannot be described by a single determinant. The determinants of these states are not eigenfunctions of the S2 spin operator. To correct the wrong spin coupling, we can break the symmetry in calculation. For example, two unpaired electrons are located at the two centers, Heisenberg spin coupling of the two electrons: H = 2JSASB, J>0 for anti-ferromagnetic coupling.

The high spin state is a pure spin state Smax=SA+SB; For the broken symmetry state, it is a mixture of the spin states,

MS=| SA-SB |.

From EHS - EBS = - 4JSASB, weak interaction limit J is determined.

Using Lande interval rule ΔES, (S-1) = - 2JS, we get the energy of the true low spin state ELS.

Explain how symmetry can be used in quantum chemistry. (2p)

Symmetry operation can transform the molecule into a final position that is indistinguishable from the initial position. Symmetry element is geometrical entity corresponding to which a symmetry operation is carried out. First, we determine the point group of the molecule (symmetry elements), then classify the molecular orbitals in irreducible representations of the point group, third, we should choose electronic configuration in terms of these orbitals, and finally determine the total electronic state (from open shells only) in terms of irreducible representation and spin.

Potential energy surfaces and transition state theory

a. A chemical reaction occurs in solution: A->B->C. The relative energies are the following (all related to the reactant A): A: 0 kcal/mol, B: -5 kcal/mol, C: -8 kcal/mol, TSAB : 12 kcal/mol, TSBC: 10 kcal/mol (where TSAB means the transition state between A and B). Which is the rate limiting transition state and what is the activation energy? Motivate your answer. (l p)

The rate limiting transition state is the TSBC, because TSBC-EB>TSAB. and the activation energy is 10-(-5)=15 kcal/mol, it is the energy that must be overcome in order for a chemical reaction to occur.

b. A chemical reaction occurs in solution: A->B->C->D. The relative energies are the following (all related to the reactant A): A: 0 kcal/mol, B: -8 kcal/mol, C: -5 kcal/mol, D: -10 kcal/mol, TSAB: 14 kcal/mol, TSBC: 4 kcal/mol, TSCD: 8 kcal/mol (where TSAB means the transition state between A and B). Which is the rate limiting transition state and what is lhe activation energy? Motivate your answer. (l p)

The rate limiting transition state is the TSAB, and the activation energy is 14-0=14 kcal/mol, it is the energy that must be overcome in order for a chemical reaction to occur.

c. Which are the assumptions made in transition state theory? (1 p)

There is an equilibrium energy distribution in all the possible quantum states at all points along the reaction coordinates, and the reacting molecule is in thermal equilibrium with the environment; atomic nuclei move according to classic mechanics, Stable molecules reach minima on the PES and chemical reactions are nuclear motions from one minimum to another. Transition State is an intermediate first-order saddle point. All molecules passing through TS will go on to form product, with no reflection to form the initial reactant.

d. Describe which type of information is needed to calculate macroscopic properties like Gibbs free energy, and how this information is obtained from quantum chemical calculations. In particular, how is the zero point correction to the enthalpy calculated? (3p)

Gibbs free energy can be got by G = H -TS. For a macroscopic ensemble of particles, the H and S may be calculated by statistical mechanics from properties of a relatively few molecular.

Then we get the partition function Q=exp(-E(r,p)/kT)dr dp, according to the thermodynamics, we get many functions: G = H -TS=kTV(lnQ/V) T -kTlnQ

From QC results, we can get the zero-point corrected results by Hessian matrix of the optimized structure.

The zero-point corrected energy is the vibration energy of the system in 0 K and in general, the molecular are not at this temperature, so correction is needed.

e. The experimental time-constant at room temperature for a species is 500μs.How large is the free activation energy? (Show your calculations.) (l p)

k = (kBT/h)* exp(-ΔG#/RT) => ΔG# = -RT ln[k/( kBT/h)], k = 1/500 [μs-1] = 2000 [s-1],

For 300K, RT = 0.592 kcal/mol, kBT/h = 6.2Ã-1012 [s-1],

So, the free activation energy ΔG# = -0.592 * ln[2000/(6.2*1012)] = 12.94 kcal/mol.

f. One reaction is 10000 times faster than another one. What is the difference in barrier heights (in kcal/mol) for the two reactions? (Show your calculations.) (1 p)

k = kBT/h * exp(-ÄG#/RT) so, k1/ k2 = 10000, ΔG2#-ΔG1# = RT*ln(k1/ k2) = 9.2*0.592=5.45kcal/mol

The slow reaction barrier height is 5.45kcal/mol bigger than the fast reaction barrier height.

Your quantum chemical calculations give the following results:

Reactant TS

Electronic energy 0 15.0 kcal/mol

Zero point energy 28.0 22.0 kcal/mol

Entropy (S) 78.0 81.0 cal/mol/Kelvin

Calculate the rate-constant for the reaction at room temperature. (Neglect the temperature dependence of the enthalpy.) Show your calculations. (1p)

Gibbs free energy corrections: G# = H# -TS# = 15.0+22.0-300*81.0*0.001 = 12.7 kcal/mol

G0 = H0 -TS0 = 0+28.0-300*78.0*0.001 = 4.6 kcal/mol

The reaction rate at 300K: k = kBT/h * exp(-ÄG#/RT) = 6.2Ã-1012 * exp[-(12.7-4.6)/0.592] = 7.1Ã-106 [s-1]

4. Molecular mechanics

a. Force fields can be used as a complement or an augmentation to quantum chemistry. Give two different examples of this and describe what is gained. (2p)

The major problems in QC is calculating the electronic energy for a given nuclear configuration to give a potential energy surface. This step is bypassed in force field (FF) by writing the electronic energy as a parametric function of the nuclear coordinates, and fitting the parameters to experimental or higher level computational data.

The quantum aspects of the nuclear motion are also neglected in Force field, which means that the dynamics of the atoms is treated by classical mechanics, i.e. Newton's second law.

So, it is relatively fast, and can calculate large systems combine with quantum mechanics, and study dynamical information.

b. Write the energy expression that is used in force field (FF) methods, and explain the different terms. Which are the basic assumptions that these methods rely on? (4p)

EFF = Estr + Ebend + Etors + Evdw + Eel + Ecross

Estr: bond stretching energy

Ebend: angle bending energy

Etors: torsion (bond rotation) energy

Evdw: van der Waals interaction energy

Eel: electrostatic energy

Ecross: couplings between Estr, Ebend, Etors

Assumption: For a particular pair of atom types i) the optimal bond distance is the same, and ii) the energy change upon distortion is the same (force constant).

c. How would an ab initio FF and a semi-empirical FF be derived for a mixture of water and methanol? (2p)

We must decide how to describe them in terms of the selected force field. Then decide three sets of information: atom types, how are they connected, and start guess of the geometry. i) empirical FF comes from different types of experimental data and ii) ab initio data is from electronic structure calculations.

d. Describe the limitations of force field methods? (2p)

The accuracy and predictability may be not so good as the lack of good parameters. Transferability and reactions are very difficult.

e. What is a hybrid quantum mechanics - molecular mechanics method? Describe also the ONIOM method available in the Gaussian package. (2p)

A hybrid QM/MM method is always applied to large systems. The active part is calculated by electronic structure methods (i.e. quantum mechanics) such as: ab initio or DFT, while the rest part-the backbone is treated in force field methods (i.e. molecular mechanics). This is the hybrid QM/MM method.

EQM/MM = EQMmodel + EMMr-m + EMM(r-m)*model, EQMmodel is the energy of the important group; EMMr-m is the energy of the backbone, and EMM(r-m)*model is the interaction between them.

ONIOM method: For a two-layer, the small system is calculated at both the low and high levels of theory, while the large system is calculated at the low level of theory. The energy in ONIOM EONIOM = EQMmodel + EMMreal - EMMmodel. The ONIOM model can be constructed from the derivative of the underlying methods, and thus possible to perform geometry optimizations and vibrational analysis by using the ONIOM function.

5. Electron correlation methods

a. Define correlation energy and describe its physical origin. (2p)

Correlation energy is the difference between the Hartree-Fock energy calculated for a system and the exact non-relativistic energy of that system. The correlation energy comes from the approximate representation of the electron-electron repulsions in the Hartree-Fock method.

b. What does it mean that a calculated energy is variational? Explain and give a formula. (2p)

The HF approximation has its limitation as it uses an average interaction as a substitution for the e-e interaction and it can get 99% of the total energy in a large basis set. So, electron correlation energy which is the difference between the HF and lowest energy in a basis is introduced. Ecorr=Emin-EHF.

c. Describe how the wavefunction is constructed in MCSCF and how it is optimized. (2p)

MCSCF wave function optimizations are normally carried out by expanding the energy to second order in the variational parameters, using the Newton-Raphson-based methods to force convergence to a minimum.

The MCSCF optimization is iterative as the SCF procedure. Since the number of MCSCF iterations required for achieving convergence tends to increase with the number of configurations included, the size of MCSCF wave functions that can be treated is a bit smaller than for CI methods.

d. Write down the second order energy expression for Moeller-Plesset type of perturbation theory? Use this expression and explain what is characteristic for the configurations that give the largest contributions to the second order energy? Which type of excitations contribute to the second order energy? (4p)

EMP(2) = Ói < j, a < b |<Фn(0) |Órs1/rrs |Фijab>|2/(åi+åj-åa-åb)

Difference in orbital energies in denominator: ( åi-åa )i-a, (åj-åb)j-b gives large contributions from near degeneracy.

e. The full CI wave function and the full coupled-cluster wave function are identical, but their truncations differ. Write down the general expressions for the CI and CC wave functions. Discuss and explain the differences between CISD and CCSD. (2p)

for CC: ΨCC = eT  0 Ti are excitation operators corresponding to single, double, … excitations.

ΨCCSD = e(T1+T2)  0 = (1+T1 +T2 + 1/2! (T1 +T2)2 + …)ï€ ï¦ 0

For CI: ΨCI = (1+T)ï€ ï¦ 0

ΨCISD = (1+T1 +T2)ï€ ï¦ 0

ΨCCSD is size consistent while ΨCISD is not.

ΨCCSD is size consistent because some effects of higher order excitations are included through 1/2! (T1 +T2)2

6. Density functional theory

a. Describe the two Hohenberg-Kohn theorems. (2p)

1) The external potential Vext is a functional of electron density ñ(r).

2) The ground state ñ0(r) of the electron density minimizes the energy.

b. Describe the Kohn-Sham formalism and the Kohn-Sham equations. What is achieved by this formalism, i.e. which is the main problem in orbital-free approaches that is solved? (2p)

The idea in the KS formalism is to split the kinetic energy functional into two parts, one which can be calculated exactly, and a small correction term.

For real system,

E[] = ∫vext(r) (r) dr + T[] + Vee[] =âˆ«ï€ vext(r) (r) dr + TS[] + J[] + ( T[] - TS[] ) + ( Vee[] - J[] ) = âˆ«ï€ vext(r) (r) dr + TS[] + J[] + Exc[]

Exc[ñ] = T[ñ] - TS[ñ] + Vee[ñ] - J[ñ] the exchange-correlation energy

vxc(r) = δExc[ρ]/δρ(r) Kohn-Sham equations: [-1/2

2 + vext(r) + vcoul(r) + vxc(r)] Фi(r) = εi Фi(r)

The task in developing orbital-free models is to derive approximations to the kinetic, exchange and correlation energy. The corresponding task in Kohn-Sham theory is to derive approximations to the exchange-correlation energy functional.

c. What is the local density approximation, and what is meant by gradient approximation in DFT? (2p)

In local density approximation, the density locally can be treated as a uniform electron gas.

In real, electron density in a molecule is non-uniform, so we introduce the gradient of the density. Thus gradient corrections are needed.

d. The exact non-local exchange is calculated in HF. Why is that expression not simply used to replace the approximate exchange functionals in DFT? What is meant by a hybrid method in DFT? (3p)

Non-local method is a bit misleading since the functionals only depend on the density (and derivative) at a given point, not on a space volume as the Hartree-Fock exchange energy.

Hybrid functionals in DFT: combine the HF exchange (the exact exchange) with different gradient correction terms.

e. What is the self-interaction-error, SIE, in DFf? You can explain using the simple example of H2+. How can you correct for the SIE? (2p)

The HF-energy is Eel = <Ã-|H|Ã-> = ÓiIi + 1/2 Ói,j(Jij -Kij). For i = j, one can see Coulomb repulsion Jii is exactly cancelled by the exchange energy,

While the DFT-energy is E[ñ] = ∫vext(r)ñ(r)dr + TS[ñ] + J[ñ] + Exc[ñ], Coulomb repulsion from itself is not exactly cancelled by Ex[ñ]; instead, delocalization leads to artificially low energy.

For H2+, we can see the self interaction error in DFT leads to low energy than the HF. We can choose part of the exchange energy as a correction for the self-interaction energy.

f. Why is the dispersion interaction missing in local and semi-local DFT functionals? How can it be included in the description (2p)

Weak interactions due to dispersion forces (part of van der Waals type interactions) come from electron correlation in wave function methods, but this is poorly described by current DFT methods.

The Morokuma energy decomposition is used for analyzing weak interactions; for strong interactions the mixing term often accounts for a significant part of the total interaction, thus obscuring the decomposition.

7. Environmental effects

a. Take a suitable quantum chemical problem and describe how it would be treated by implicit and explicit solvation models. What are the pros and cons of the different approaches in that application and in general? (4p)

In the continuum model, the solvent is considered implicitly as a polarizable medium characterized by its dielectric constant. It is often applied to calculate the free energy of solute-solvent interactions in chemical processes, such as folding or conformational transitions of proteins, DNA, RNA. So, it is successful in a large number of applications and is the most used model at present. The advantage is that its computational cost is similar to that of the isolated solute, and since the solvent is described implicitly, the results may be rationalized following the same way as for the case of an isolated molecule.

Although the implicit solvent model is useful for simulations of biomolecules, this is an approximate method with certain limitations related to parameterization and treatment of ionization effects. Its major difficulty is found when specific interactions such as hydrogen bonds are needed for a proper treatment.

In explicit model, all molecular interactions in the underlying many-body formulation of the quantum chemical model are included in the calculation. The drawback relates to the complexity of performing full QM calculations of such large samples.

Taking solvent effects on the electronic absorption spectrum of camphor as an example, the implicit is less successful in predicting solvatochromic effects for solvents of low polarity, but an explicit model with a significant number of solvent molecules is good.

b. How would you evaluate a quantum chemical model calculation using an implicit solvation model? Describe acetone in water. (2p)

A direct calculation of solvating acetone in water requires simulating the transfer of an acetone molecule from the gas phase to an aqueous phase, and then by solvent reorganization.

Implicit solvation models consider the solvent as a uniform polarizable medium with a dielectric constant of ε, and with the solute M placed in a suitably shaped hole in the medium, we must consider following aspect: size and shape of the hole, cavity/dispersion contribution, charge distribution of M, describing solute M and dielectric medium in either classical (force field) or quantum.

The dielectric medium normally has a constant value of ε, but for some dynamical phenomena, it can also be allowed to be frequency dependent. The ε is the only parameter characterizing the solvent, and solvents having the same ε value are thus treated equally. The hydrogen bonding capability of 1-propanol compared with acetone will be likely make a difference.

c. Describe a case in which explicit solvation is required and how the properties would be derived. (4p)

Molecular Dynamics Simulations is a case that using explicit salvation. Molecular systems are made of an assembly of interacting particles and from the Interaction potential we can get the forces between the particles, then using Newtons 2nd law of motion, positions and velocities can be got and then through statistical mechanics, we can get average physical characteristic of the system.