Le logiciel de modélistion moléculaire dont nous disposons, CAChe récement renommé en SciGress par Fujitsu, propose un module d'analyse conformationnel appelé ConFlex qui permet l'exploration rapide et exaustif de l'espace conformationnel, malheureusement ConFlex est limité au vieillisant champ de force MM3. Le problème posé par Bertrand Castro est relatif à la récupération des conformations issues de ConFlex et de leur retraitement en utilisant d'autres méthodes par exemple MOPAC ou Gaussian qui sont deux programmes de chimie quantiques. Dans la mesure où Gaussian et MOPAC sont capables de déterminer de parametres thermodynamique, dont le plus important : l'enthalpie libre de Gibbs (G), serait il possible d'utiliser les conformères trouvés pour déterminer les propriétés thermodynamiques de l'espece chimique et non d'un conformère unique (dans le vide).
Le module ConFlex[ref] inclu dans SciGress propose une analyse conformationnelle en effectuant des inversion chaises/bateaux et en tournant les angles dihèdres des chaines non cycliques. L'algorythme est très performant et permet de générer en quelques heures toutes les conformation probables pour une petite moléule «druglike» (il faut compter une semaine pour un tétrasaccharide). Malheureusement SciGress est limité à MM3 [Allinger] qui est paramétré pour les hydrocarbures et n'est pas très précis. Après quelques échanges avec l'editeur du logicielSciGress, ce dernier nous a donné un moyen de faire afficher les coordonnées des atomes de chaque conformère (initialement dans le fichier project.io/conflex.ls3). ConFlex utilise un moteur basé sur les angles dièdres pour déterminer si deux conformères sont identiques.
Chaque conformère ainsi extrait par le programme ExtractConflex.py (Voir le code source, en annexe) peut être optimisé en géométrie à l'aide d'une méthode plus précise que MM3, par exemple un autre champ de force (MMFF94s), MOPAC ou Gaussian.
Le Project Leader, un des modules de SciGress, est tout a fait adapté au retraitement par lot des structures extraites par ExtractConflex.py. Le défaut principal de cette méthode vient de l'absence de paralèlisation des calculs ce qui prends des semaines de calculs.
La méthode de Kabsch a été publiée en 1976 (et corrigée en 1978) et propose un moyen de superposer un jeu de vecteurs (les positions atomiques) et détermine l'écart type en distance (RMSD). La première étape consiste en un changement de repère pour se placer au niveau du barycentre moléculaire de chacune des molécules (non pondéré par la masse). La seconde consiste à trouver la matrice de rotation permettant la minimisation du RMSD. Si le RMSD calculé est inférieur à une valeur limite (choisi à 0.1Å), les deux conformères sont considérés comme identiques est celui qui a l'énergie la plus faible est conservé. Cette méthode présente une problème lorsque un fragment moléculaire présente une symétrie interne. Par exemple un groupe phenyl qui est inva
cacheparser.py (à intégrer dans cclib0.6.1, disponible surhttp:cclib.sourceforge.net)
ExtractConflex.py
GridPostProcess.py
===== De l'énergie électronique aux propriétés thermodynamiques : Utilisation de la thermodynamique Statisitque =====
Thermodynamic quantities, such as heat capacity, entropy, and internal energy, can be calculated using the vibrational frequencies (energies), moments of inertia of the molecule, its symmetry number, and a knowledge of the temperature. In this section, the relationships of these quantities is described.Some derived quantities, which will be used in this section only are: Moment of inertia:I(1 amu Ångstrom^2^= 1.660540×10^-40^g cm^2^).Rotational constants: A, B, and C (e.g. ) A = h/(8p^2^I) With I in amu Ångstroms^2^ then: A (in MHz) = 5.053791×10^5^/I A (in cm^-1^) = 5.053791×10^5^/cI = 16.85763/I
====== Thermochemistry from ab initio MO methods ======
Ab initio MO methods provide total energies, E,,eq,,, as the sum of electronic and nuclear-nuclear repulsion energies for molecules, isolated in vacuum, without vibration at 0 K.
From the 0 K potential surface and using the harmonic oscillator approximation, we can calculate the vibrational frequencies, n,,i,,, of the normal modes of vibration. Using these, we can calculate vibrational, rotational and translational contributions to the thermodynamic quantities such as the partition function and heat capacity which arise from heating the system from 0 to T K.
Q: partition function, E: energy, S: entropy, and C: Heat capacity at constant pressure = C,,p,,. In ab initio calculations, the heat capacity calculated is C,,v,,. The relationship between C,,p,, and C,,v,, (in cal.degree^-1^.mol^-1^) is:
===== Vibrational terms =====
The vibrational contribution to the internal energy arises from population of the vibrational energy levels. The vibrational partition coefficient, Q,,vib,,, is given by:
E,,vib,,, for a molecule at the temperature T as:
where h is Planck's constant, n,,i ,,the i-th normal vibration frequency, and k the Boltzmann constant. For 1 mole of molecules, E,,vib ,,should be multiplied by the Avogadro number N,,a,, = R/k. Thus:
||<tablewidth=“100%”>
||Equation 1 ||
Note that the first term in the above equation is the zero-point vibration energy, E,,zpe,,. Hence, the second term is the additional vibrational contribution due to the temperature increase from 0 K to T K. Namely,
||<tablewidth=“100%”>E,,vib,, ||<5 %>= || E,,zpe,,+,, ,,E,,vib ,,(T) ||<style=“vertical-align: top;”> ||
||E,,zpe,, ||<5 %>= ||
||<style=“vertical-align: top;”> ||
||E,,vib,,(T) ||<5 %>= ||
||Equation 2 ||
The value of E,,vib ,,from GAUSSIAN 82 and 86 includes E,,zpe ,,as defined by Equation 1 and Equation 2.
||<tablewidth=“100%“8 %>S,,vib,, ||<5 %>= ||
||
||<8 %>C,,vib,, ||<5 %>= ||
||
At temperature T>0 K, a molecule rotates about the x, y, and z-axes and translates in x, y, and z-directions. By assuming the equipartition of energy, energies for rotation and translation, E,,rot ,,and E,,tr,,, are calculated.
===== Rotational terms =====
s is the symmetry number (Examples of symmetry numbers are shown in the Table). I is moment of inertia. I,,A,,, I,,B,,, and I,,C,, are moments of inertia about A, B, and C axes.
Table: Table of Symmetry Numbers
||<tablewidth=“100%“7 %>C,,1,, ||<9 %>C,,I,, ||C,,S,,: ||<5 %>1 ||<5 %> ||<7 %>D,,2,, ||<9 %>D,,2,,,,d,, ||D,,2,,,,h,,: ||<7 %>4 ||<5 %> ||C
: ||<8 %>1 ||
||<7 %>C,,2,, ||<9 %>C,,2,,,,v,, ||C,,2,,,,h,,: ||<5 %>2 ||<5 %> ||<7 %>D,,3,, ||<9 %>D,,3,,,,d,, ||D,,3,,,,h,,: ||<7 %>6 ||<5 %> ||D
: ||<8 %>2 ||
||<7 %>C,,3,, ||<9 %>C,,3,,,,v,, ||C,,3,,,,h,,: ||<5 %>3 ||<5 %> ||<7 %>D,,4,, ||<9 %>D,,4,,,,d,, ||D,,4,,,,h,,: ||<7 %>8 ||<5 %> ||T, T,,h,, T,,d,,: ||<8 %>12 ||
||<7 %>C,,4,, ||<9 %>C,,4,,,,v,, ||C,,4,,,,h,,: ||<5 %>4 ||<5 %> ||<7 %>D,,5,, ||<9 %>D,,5,,,,d,, ||D,,5,,,,h,,: ||<7 %>10 ||<5 %> ||O, O,,h,,: ||<8 %>24 ||
||<7 %>C,,5,, ||<9 %>C,,5,,,,v,, ||C,,5,,,,h,,: ||<5 %>5 ||<5 %> ||<7 %>D,,6,, ||<9 %>D,,6,,,,d,, ||D,,6,,,,h,,: ||<7 %>12 ||<5 %> ||I, I,,h,,: ||<8 %>60 ||
||<7 %>C,,6,, ||<9 %>C,,6,,,,v,, ||C,,6,,,,h,,: ||<5 %>6 ||<5 %> ||<7 %>D,,7,, ||<9 %>D,,7,,,,d,, ||D,,7,,,,h,,: ||<7 %>14 ||<5 %> ||S,,4,,: ||<8 %>2 ||
||<7 %>C,,7,, ||<9 %>C,,7,,,,v,, ||C,,7,,,,h,,: ||<5 %>7 ||<5 %> ||<7 %>D,,8,, ||<9 %>D,,8,,,,d,, ||D,,8,,,,h,,: ||<7 %>16 ||<5 %> ||S,,6,,: ||<8 %>3 ||
||<7 %>C,,8,, ||<9 %>C,,8,,,,v,, ||C,,8,,,,h,,: ||<5 %>8 ||<5 %> ||<7 %> ||<9 %> || ||<7 %> ||<5 %> ||S,,8,,: ||<8 %>4 ||
==== Linear molecule ====
Values for Q,,rot,,, E,,rot,,, and S,,rot,, for a linear molecule are defined below:
||<tablewidth=“100%”>Q,,rot,, ||<6 %>= ||
||
||E,,rot,, ||<6 %>= ||(2/2)RT ||
||S,,rot,, ||<6 %>= ||
||
.
||<6 %>= ||
||
where
.
==== Non-linear molecule ====
Values for Q,,rot,,, E,,rot,,, and S,,rot,, for a non-linear molecule are defined below:
||<tablewidth=“100%“8 %>Q,,rot,, ||<3 %>= ||
||
.
||<3 %>= ||
||
||<8 %>E,,rot,, ||<3 %>= ||(3/2)RT ||
||<8 %>S,,rot,, ||<3 %>= ||
||
.
||<3 %>= ||
||
Here, -5.386 3921 is calculated as:
===== Translational terms =====
Given that M is the molecular weight, then the values for Q,,tra,,, E,,tra,,, S,,tra,,, and C,,tra ,,for a molecule are as defined below:
||<tablewidth=“457px”width=“37px”>Q,,tra,, ||<width=“19px”>= ||<width=“389px”>
||
||<width=“37px”>E,,tra,, ||<width=“19px”>= ||<width=“389px”>(3/2)RT ||
||<width=“37px”>S,,tra,, ||<width=“19px”>= ||<width=“389px”>
||
.
||<width=“19px”>= ||<width=“389px”>
||
||<width=“37px”>C,,tra,, ||<width=“19px”>= ||<width=“389px”>(5/2)R ||
or H,,tra ,,= (5/2)RT due to the PV term (cf. H = U + PV). The internal energy U at T is:
U = E,,eq,, + [E,,vib ,,+ E,,rot,, + E,,tra,,]
or, in terms of the zero point energy, E,,zpe,,, and real vibrations, E,,vib ,,(T),
U = E,,eq,, + [(E,,zpe ,,+,, ,,E,,vib ,,(T))+ E,,rot,, + E,,tra,,].
Enthalpy H for one mole of gas is defined as
H = U + PV
Assumption of an ideal gas (i.e., PV = RT) leads to
H = U + PV = U + RT
Thus, Gibbs free energy G can be calculated as:
G = H - T S
====== Thermochemistry in MOPAC ======
It should be noted that M.O. parameters for MNDO, AM1, etc., are optimized so as to reproduce the experimental heat of formation (i.e., standard enthalpy of formation or the enthalpy change to form a mole of compound at 25^o^C from its elements in their standard state) as well as observed geometries (mostly at 25^o^C), and not to reproduce the E,,eq,, and equilibrium geometry at 0 K.
In this sense, E,,SCF,, (defined as Heat of formation, DH,,f,,), force constants, normal vibration frequencies, etc. are all related to the values at 25^o^C, not to 0 K. Therefore, the E,,0 ,,calculated in FORCE is not the true E,,0,,. Its use as E,,0 ,,should be made at your own risk, bearing in mind the situation discussed above.
Since E,,SCF,, is standard enthalpy of formation (at 25
C):
E,,SCF,, = E,,eq,, + E,,zpe,, + E,,vib,,(298.15) + E,,rot,, + E,,tra,, + PV + S[-E,,elec,,(atom) + DH,,f,,(atom)].
To avoid the complication arising from the definition of E,,SCF,,, within the thermodynamics calculation the Standard Enthalpy of Formation, DH, is calculated by
DH = E,,SCF,, + (H,,T,, - H,,298,,).
Here, E,,SCF,, is the heat of formation (at 25^o^C) given in the output list, and H,,T,, and H,,298,, are the enthalpy contributions for the increase of the temperature from 0 K to T and 298.15, respectively. In other words, the enthalpy of formation is corrected for the difference in temperature from 298.15 K to T.
There is a problem in that H,,T,, is the heat of formation at T relative to the heat of formation of the elements in their standard state at 298K. This involves mixing standard and not standard terms. There is no easy way to get the correct value for H,,T,,, but for rough work H,,T,, is useful. For more correct work, calculate DH for the elements in their standard state at T, and use these DH's to get the DH for the compound you're working with (or use tables from the literature).
This problem is, however, not normally important, because the most common use of H,,T,, is for calculating the thermodynamics of reactions at temperatures other than 298K. For all reactions, the types and number of atoms must be the same in reactants and products, therefore the fact that the H,,T,, are relative to the elements in their standard state at 298K is irrelevant. Consider the simple Diels-Alder reaction:
C,,2,,H,,4,, + C,,4,,H,,6,, = C,,6,,H,,10,,
The heat of this reaction at 298K is H,,298,,( C,,6,,H,,10,,) - H,,298,,(C,,2,,H,,4,,) - H,,298,,( C,,4,,H,,6,,). At any other temperature, the heat of reaction would be:
H,,R ,,=,, ,,H,,T,,( C,,6,,H,,10,,) - H,,T,,(C,,2,,H,,4,,) - H,,T,,( C,,4,,H,,6,,).
Care must, however, be taken to account for changes in volume - if any of the reactants or products are gaseous, then appropriate corrections must be made to H,,R,, . Complications arise only if absolute heats of formation are needed. Thus, if the heat of formation of benzoic acid (C,,7,,H,,6,,O,,2,,) at 398K (100C) is needed, the H,,398,,(C,,7,,H,,6,,O,,2,,) generated by MOPAC would be for the reaction:
7H,,298,,(graphite) + 3H,,298,,(H,,2,,) + H,,298,,(O,,2,,) ⇒ H,,398,,(C,,7,,H,,6,,O,,2,,)
. 0 + 0 + 0 ⇒H,,398,,(C,,7,,H,,6,,O,,2,,)
Note that on the left side, the temperatures are 298K. For H,,2,, and O,,2,,, the heats of formation at 398K can readily be calculated, but for graphite the calculation is more complicated. The easiest way to generate a balanced equation would be to use tables of heats of formation of the elements at non-standard temperatures.
Finally, as mentioned above, changes in volume must also be taken into account: if the reaction volume changes, then RDN(T-298) must be added or subtracted, where R is the gas constant (~ 2cals/degree/mol), and DN is the change in volume. Thus for the formation of methane from graphite and hydrogen, 2 volumes of reactant (H,,2,, + graphite) yield 1 volume of methane, therefore DN = 1.
The method of calculation for T and H,,298,, will be given below.
In MOPAC, the variables defined below are used: