Method Article
Les concentrations atmosphériques de grappes moléculaires faiblement liées peuvent être calculées à partir des propriétés thermochimiques des structures à faible énergie trouvées grâce à une méthodologie d’échantillonnage configuration en plusieurs étapes utilisant un algorithme génétique et une chimie quantique semi-empirique et ab initio.
L’étude informatique de la formation et de la croissance des aérosols atmosphériques nécessite une surface d’énergie sans Gibbs précise, qui peut être obtenue à partir de la structure électronique de phase gazeuse et des calculs de fréquence vibratoire. Ces quantités sont valables pour les amas atmosphériques dont les géométries correspondent au minimum sur leurs surfaces énergétiques potentielles. L’énergie libre Gibbs de la structure énergétique minimale peut être utilisée pour prédire les concentrations atmosphériques de l’amas dans une variété de conditions telles que la température et la pression. Nous présentons une procédure computationnellement peu coûteuse construite sur un échantillonnage configurationnel basé sur un algorithme génétique suivi d’une série de calculs de dépistage de plus en plus précis. La procédure commence par générer et faire évoluer les géométries d’un grand ensemble de configurations à l’aide de modèles semi-empiriques puis affine les structures uniques résultantes à une série de niveaux ab initio de haut niveau de la théorie. Enfin, les corrections thermodynamiques sont calculées pour l’ensemble résultant de structures d’énergie minimale et utilisées pour calculer les énergies libres Gibbs de formation, les constantes d’équilibre, et les concentrations atmosphériques. Nous présentons l’application de cette procédure à l’étude des grappes de glycine hydratées dans des conditions ambiantes.
Le paramètre le plus incertain dans les études atmosphériques sur le changement climatique est la mesure exacte dans laquelle les particules de nuages reflètent le rayonnement solaire entrant. Les aérosols, qui sont des particules suspendues dans un gaz, forment des particules de nuages appelées noyaux de condensation des nuages (CCN) qui dispersent le rayonnement entrant, empêchant ainsi son absorption et le chauffage subséquent de l’atmosphère1. Une compréhension détaillée de cet effet de refroidissement net nécessite une compréhension de la croissance des aérosols en CCN, ce qui nécessite une compréhension de la croissance des petits amas moléculaires en particules d’aérosols. Des travaux récents ont suggéré que la formation d’aérosols est initiée par des grappes moléculaires de 3 nm de diamètre ou moins2; cependant, ce régime de taille est difficile d’accès en utilisant des techniques expérimentales3,4. Par conséquent, une approche de modélisation computationnelle est souhaitée afin de surmonter cette limitation expérimentale.
En utilisant notre approche de modélisation décrite ci-dessous, nous pouvons analyser la croissance de n’importe quel cluster hydraté. Parce que nous nous intéressons au rôle de l’eau dans la formation de grandes molécules biologiques à partir de petits constituants dans des environnements prébiotiques, nous illustrerons notre approche avec la glycine. Les défis rencontrés et les outils nécessaires pour répondre à ces questions de recherche sont très similaires à ceux impliqués dans l’étude des aérosols atmosphériques et des grappes de prénucléation5,6,7,8,9,10,11,12,13,14,15. Ici, nous examinons les amas de glycine hydratés à partir d’une molécule de glycine isolée suivie d’une série d’ajouts étape de jusqu’à cinq molécules d’eau. L’objectif final est de calculer les concentrations d’équilibre des grappes Gly (H2O)n-0-5 dans l’atmosphère à température ambiante au niveau de la mer et une humidité relative (RH) de 100 %.
Un petit nombre de ces amas moléculaires sous-nanomètres se développent en un amas critique métastable (1-3 nm de diamètre) soit en ajoutant d’autres molécules de vapeur ou en coagulant sur des amas existants. Ces grappes critiques ont un profil de croissance favorable menant à la formation de noyaux de condensation de nuages beaucoup plus grands (jusqu’à 50-100 nm) (CCN), qui affectent directement l’efficacité des précipitations des nuages ainsi que leur capacité à refléter la lumière incidente. Par conséquent, avoir une bonne compréhension de la thermodynamique des amas moléculaires et de leurs distributions d’équilibre devrait conduire à des prédictions plus précises de l’impact des aérosols sur le climat mondial.
Un modèle descriptif de formation d’aérosols nécessite une thermodynamique précise de la formation de grappes moléculaires. Le calcul de la thermodynamique précise de la formation des grappes moléculaires nécessite l’identification des configurations les plus stables, ce qui implique de trouver les minima globaux et locaux sur la surface énergétique potentielle de l’amas (PES)16. Ce processus est appelé échantillonnage de configuration et peut être réalisé grâce à une variété de techniques, y compris celles basées sur la dynamique moléculaire (MD)17,18,19,20, Monte Carlo (MC)21,22, et les algorithmes génétiques (GA)23,24,25.
Différents protocoles ont été développés au fil des ans pour obtenir la structure et la thermodynamique des hydrates atmosphériques à un niveau élevé de théorie. Ces protocoles différaient dans le choix de (i) méthode d’échantillonnage configurationnel, (ii) la nature de la méthode de bas niveau utilisée dans l’échantillonnage de configuration, et (iii) la hiérarchie des méthodes de niveau supérieur utilisées pour affiner les résultats dans les étapes suivantes.
Les méthodes d’échantillonnage configurationnels comprenaient l’intuition chimique26, l’échantillonnage aléatoire27,28, la dynamique moléculaire (MD)29,30, le saut de bassin (BH)31, et l’algorithme génétique (GA)24,25,32. Les méthodes de faible niveau les plus courantes utilisées avec ces méthodes d’échantillonnage sont les champs de force ou les modèles semi-empiriques tels que pm6, PM7 et SCC-DFTB. Ceux-ci sont souvent suivis par des calculs DFT avec des ensembles de base de plus en plus grands et des fonctions plus fiables des échelons supérieurs de l’échelle de Jacob33. Dans certains cas, ceux-ci sont suivis par des méthodes de fonction d’onde de niveau supérieur telles que MP2, CCSD (T), et le DLPNO-CCSD(T)34,35.
Kildgaard et coll.36 ont développé une méthode systématique où des molécules d’eau sont ajoutées à des points sur les sphères Fibonacci37 autour de plus petits amas hydratés ou non hydratés pour générer des candidats pour de plus grands groupes. Les candidats non physiques et redondants sont supprimés en fonction des seuils de contact étroits et de la distance carrée par les racines entre les différents conformistes. Les optimisations ultérieures utilisant la méthode semi-empirique PM6 et une hiérarchie des méthodes DFT et de fonction d’onde sont utilisées pour obtenir un ensemble de conformations à faible énergie à un niveau élevé de théorie.
L’algorithme38 de la colonie d’abeilles artificielles (ABC) est une nouvelle approche d’échantillonnage configurationnel qui a récemment été mise en œuvre par Zhang et coll. pour étudier les grappes moléculaires dans le cadre d’un programme appelé ABCluster39. Kubecka et coll.40 ont utilisé ABCluster pour l’échantillonnage configuration, suivis de réoptimisations de bas niveau à l’aide de la méthode semi-empirique GFN-xTB41. Ils ont en outre affiné les structures et les énergies à l’aide de méthodes DFT suivies d’énergies finales à l’aide de DLPNO-CCSD (T).
Quelle que soit la méthode, l’échantillonnage de configuration commence par une distribution aléatoire ou non-générée de points sur le PSE. Chaque point correspond à une géométrie spécifique du cluster moléculaire en question et est généré par la méthode d’échantillonnage. Ensuite, le minimum local le plus proche se trouve pour chaque point en suivant la direction «descente» sur le PSE. L’ensemble de minima ainsi trouvé correspondent à ces géométries de l’amas moléculaire qui sont stables, au moins pendant un certain temps. Ici, la forme du PSE et l’évaluation de l’énergie à chaque point de la surface seront sensibles à la description physique du système où une description physique plus précise entraîne un calcul de l’énergie plus coûteux. Nous utiliserons spécifiquement la méthode GA mise en œuvre dans le programme OGOLEM25, qui a été appliquée avec succès à une variété de problèmes d’optimisation globale et d’échantillonnage de configuration42,43,44,45, pour générer l’ensemble initial de points d’échantillonnage. Le PSE sera décrit par le modèlePM7 46 mis en œuvre dans le programme MOPAC201647. Cette combinaison est utilisée parce qu’elle génère une plus grande variété de points par rapport aux méthodes MD et MC et trouve les minima locaux plus rapidement que les descriptions plus détaillées du PSE.
L’ensemble des minima locaux optimisés par les GA sont pris comme les géométries de départ pour une série d’étapes de dépistage, qui conduisent à un ensemble d’énergie minimale de faible altitude. Cette partie du protocole commence par l’optimisation de l’ensemble de structures uniques optimisées par les GA à l’aide de la théorie fonctionnelle de la densité (DFT) avec un petit ensemble de base. Cet ensemble d’optimisations donnera généralement un plus petit ensemble de structures minimales locales uniques qui sont modélisés plus en détail par rapport aux structures semi-empiriques optimisées par les GA. Ensuite, une autre série d’optimisations DFT sont effectuées sur ce plus petit ensemble de structures en utilisant un ensemble de base plus large. Encore une fois, cette étape donnera généralement un plus petit ensemble de structures uniques qui sont modélisés plus en détail par rapport à la petite base DFT étape. L’ensemble final de structures uniques est ensuite optimisé pour une convergence plus serrée et les fréquences vibratoires harmoniques sont calculées. Après cette étape, nous avons tout ce dont nous avons besoin pour calculer les concentrations d’équilibre des grappes dans l’atmosphère. L’approche globale est résumée schématiquement à la figure 1. Nous utiliserons la corrélation d’échange PW9148 de gradient généralisé (GGA) fonctionnelle dans la mise en œuvre de la DFT gaussian0949 ainsi que deux variantes de l’ensemble de base Pople50 (6-31 GD POUR la petite étape de base et 6-311 GD pour la grande étape de base). Cette combinaison particulière de l’ensemble fonctionnel et de base de corrélation d’échange a été choisie en raison de son succès précédent dans l’informatique précis Gibbs énergies libres de formation pour les amas atmosphériques51,52.
Ce protocole suppose que l’utilisateur a accès à un cluster informatique haute performance (HPC) avec le système de lots portable53 (PBS), MOPAC2016 (http://openmopac.net/MOPAC2016.html)47, OGOLEM (https://www.ogolem.org)25, Gaussian 09 (https://gaussian.com)49, et OpenBabel54 (http://openbabel.org/wiki/Main_Page) logiciel installé suivant leurs instructions d’installation spécifiques. Chaque étape de ce protocole utilise également un ensemble de scripts en coquille interne et Python 2.7 qui doivent être enregistrés à un répertoire qui est inclus dans la variable environnementale $PATH de l’utilisateur. Tous les modules environnementaux nécessaires et les autorisations d’exécution pour exécuter tous les programmes ci-dessus doivent également être chargés dans la session de l’utilisateur. L’utilisation du disque et de la mémoire par le code GA (OGOLEM) et les codes semi-empiriques (MOPAC) sont très faibles par les normes modernes de ressources informatiques. La mémoire globale et l’utilisation du disque pour OGOLEM/MOPAC dépend du nombre de threads que l’on veut utiliser, et même alors, l’utilisation des ressources sera faible par rapport aux capacités de la plupart des systèmes HPC. Les besoins en ressources des méthodes de QM dépendent de la taille des grappes et du niveau de théorie utilisé. L’avantage de l’utilisation de ce protocole est que l’on peut varier le niveau de théorie pour être en mesure de calculer l’ensemble final de structures à faible énergie, en gardant à l’esprit que les calculs généralement plus rapides conduisent à plus d’incertitude dans l’exactitude des résultats.
Par souci de clarté, l’ordinateur local de l’utilisateur sera appelé «ordinateur local» tandis que le cluster HPC auquel il a accès sera appelé «cluster à distance».
1. Trouver la structure énergétique minimale de la glycine et de l’eau isolées
REMARQUE : L’objectif ici est double : (i) d’obtenir des structures énergétiques minimales de molécules d’eau et de glycine isolées pour une utilisation dans l’échantillonnage de configuration de l’algorithme génétique, (ii) et de calculer les corrections thermodynamiques aux énergies de phase gazeuse de ces molécules pour une utilisation dans le calcul des concentrations atmosphériques.
2. Échantillonnage de configuration basé sur l’algorithme génétique des clusters Gly (H2O)n-1-5
REMARQUE: L’objectif ici est d’obtenir un ensemble de structures à faible consommation d’énergie pour Gly (H2O)n '1-5 au niveau semi-empirique peu coûteux de la théorie, en utilisant le modèle PM746 mis en œuvre dans MOPAC47. Il est impératif que le répertoire de travail ait l’organisation et la structure exactes, comme le montre la figure 2. Il s’agit de s’assurer que la coque personnalisée et les scripts Python fonctionnent sans échecs.
3. Raffinement à l’aide de la méthode QM avec un petit ensemble de base
REMARQUE : L’objectif ici est d’affiner l’échantillonnage configuré des clusters Gly(H2O)n'1-5 en utilisant une meilleure description quantique-mécanique pour obtenir un ensemble plus petit mais plus précis de structures de cluster Gly(H2O)n-1-5. Les structures de départ de cette étape sont les sorties de l’étape 2.
4. Amélioration supplémentaire de la méthode QM avec un ensemble de base
REMARQUE : L’objectif ici est d’affiner davantage l’échantillonnage de configuration des clusters Gly(H2O)n'1-5 à l’aide d’une meilleure description quantique-mécanique. Les structures de départ de cette étape sont les sorties de l’étape 3.
5. Calculs de correction finales de l’énergie et de la thermodynamique
REMARQUE : L’objectif ici est d’obtenir la structure vibratoire et les énergies des clusters Gly(H2O)n'1-5 à l’aide d’un grand ensemble de base et d’une grille d’intégration ultrafine afin de calculer les corrections thermochimiques souhaitées.
6. Calcul des concentrations atmosphériques des amas Gly (H2O)n-0-5 à température ambiante au niveau de la mer
REMARQUE : Ceci est accompli en copiant d’abord les données thermodynamiques générées dans l’étape précédente dans une feuille de calcul et en calculant les énergies libres Gibbs de l’hydratation séquentielle. Ensuite, les énergies libres Gibbs sont utilisées pour calculer les constantes d’équilibre pour chaque hydratation séquentielle. Enfin, un ensemble d’équations linéaires sont résolus pour obtenir la concentration d’équilibre des hydrates pour une concentration donnée de monomères, de température et de pression.
La première série de résultats de ce protocole devrait être un ensemble de structures à faible consommation d’énergie de Gly (H2O)n'1-5 trouvé par la procédure d’échantillonnage configurationnel. Ces structures ont été optimisées au niveau de théorie PW91/6-311 et sont supposées exactes aux fins du présent document. Rien n’indique que PW91/6-311 GMD sous-estime ou surestime constamment l’énergie contraignante de ces grappes. Sa capacité à prédire les énergies contraignantes par rapport à MP2/CBS32 et [DLPNO-]CCSD(T)/CBS60,61 estimations et l’expérience52 montre beaucoup de fluctuations. Il en va de même pour la plupart des autres fonctions de densité. En général, chaque valeur de n 1 à 5 devrait donner une poignée de structures à faible énergie dans un rayon d’environ 5 kcal mol-1 de la structure la plus faible en énergie. Ici, nous nous concentrons sur la première structure produite par le script run-thermo-pw91.csh pour la brièveté. La figure 3 montre les plus faibles isomères d’énergie électronique des grappes Gly (H2O)n-0-5. On peut voir que le réseau de liaisons d’hydrogène augmente en complexité à mesure que le nombre de molécules d’eau augmente, et passe même d’un réseau principalement planaire à une structure tridimensionnelle en forme de cage à n 5. Le reste de ce texte utilise les énergies et les quantités thermodynamiques correspondant à ces cinq clusters spécifiques.
Le tableau 1 contient les quantités thermodynamiques nécessaires à l’application du protocole. Le tableau 2 montre un exemple de la sortie du script run-thermo-pw91.csh où les énergies électroniques, les corrections vibratoires à zéro point et les corrections thermodynamiques à trois températures différentes sont imprimées. Pour chaque cluster (ligne), le niveau de théorie E[PW91/6-311-G] correspond aux énergies électroniques de phase gazeuse du niveau de théorie PW91/6-311-GMD calculé sur les réseaux d’intégration ultrafin dans les unités de Hartree, ainsi qu’à l’énergie vibratoire zéro point(ZPVE) dans les unités de m mol kcal-1. À chaque température, 216,65 K, 273,15 K, et 298,15 K, les corrections thermodynamiques sont énumérées, H l’enthalpie de formation dans les unités de kcal mol-1, S l’entropie de formation dans les unités de cal mol-1, et G l’énergie sans Gibbs de formation dans les unités de m mol kcal-1. Le tableau 3 montre un exemple de calcul du changement total d’énergie libre Gibbs de l’hydratation, ainsi que pour l’hydratation séquentielle. Un exemple de calcul du changement total d’énergie libre Gibbs de l’hydratation pour la réaction
commence par le calcul de l’énergie électronique EPW91 comme
où EPW91[Gly(H2O)] est tiré de la colonne C du tableau 2, et EPW91[Gly] et EPW91[H2O] sont tirés de la colonne B du tableau 1. Ensuite, nous calculons le changement énergétique de phase totale du gaz'E(0) en incluant le changement dans l’énergie vibratoire zéro point de la réaction comme
pour obtenir la colonne D. Ici, leTPS91/6-311-GMD est tiré de la colonne C du tableau 3, de l’EZPVE[Gly (H2O)] de la colonne D du tableau 2, et de l’EZPVE[Gly] et de l’EZPVE[H2O] de la colonne C du tableau 1. Pour des raisons de brièveté, nous passerons à des grappes de température ambiante, donc nous passons sur les données de 216,65 K et 273,15 K. À température ambiante, nous calculons ensuite le changement d’enthalpie de la réactionH en corrigeant le changement d’énergie de phase de gaz
où l’onpuise dans la colonne du tableau 3 D,le H[GlyMD(H2O)] est tiré de la colonne K du tableau 2, et les articlesde H[Gly] etH[H2O] sont tirés de la colonne J du tableau 1. Enfin, nous calculons le changement d’énergie libre Gibbs de la réactionG comme
où le tableau 2 est tiré de la colonne du tableau 3 I, S[GlyMD(H2O)] est tiré de la colonne L du tableau 2, et S[Gly] et S[H2O] sont tirés de la colonne K. Note ici que les valeurs d’entropie doivent être converties en unités de mol kcal-1 K-1 pendant cette étape.
Nous avons maintenant les quantités nécessaires pour calculer les concentrations atmosphériques de glycine hydratée comme le montre l’étape 6. Les résultats devraient ressembler aux données présentées au tableau 4,mais il faut s’attendre à de petites différences numériques. Le tableau 4 montre les concentrations d’hydrate d’équilibre trouvées à partir de la formulation du système de six équations de l’étape 6.2 en une seule équation de matrice et sa solution subséquente. Nous commençons par reconnaître le fait que le système d’équations peut être
où Kn est la constanted’équilibre pour l’hydratation séquentielle n de glycine, w est la concentration d’eau dans l’atmosphère, g est la concentration initiale de glycine isolée dans l’atmosphère, et gn est la concentration d’équilibre de Gly (H2O)n. Si nous réécrivons l’équation ci-dessus comme Ax b, nous obtenons x 'A'1b où A'1 est l’inverse de la matrice A. Cet inverse peut être facilement calculé à l’aide de fonctions intégrées de feuilles de calcul comme indiqué dans le tableau 4 pour obtenir les résultats finaux.
La figure 4 montre la concentration d’équilibre de la glycine hydratée calculée dans le tableau 4 en fonction de la température à 100 % d’humidité relative et à 1 pression de l’atmosphère. Il montre que, comme la température diminue de 298,15K à 216.65K, la concentration de glycine non hydratée (n '0) diminue et ceux de la glycine hydratée augmente. Le dihydrate de glycine (n-2) augmente considérablement avec la température décroissante tandis que le changement dans la concentration d’autres hydrates est moins perceptible. Ces corrélations inverses entre la température et la concentration d’hydrates sont compatibles avec l’attente que les énergies libres inférieures De Gibbs d’hydratations à des températures plus basses favorisent la formation d’hydrates.
La figure 5 illustre la dépendance relative à l’humidité de la concentration d’équilibre des hydrates de glycine à 298,15 K et 1 pression de l’atmosphère. Il démontre clairement qu’à mesure que les RH augmentent de 20 % à 100 %, la concentration d’hydrates (n-gt;0) augmente au détriment de la glycine non hydratée (n-0). Une fois de plus, la corrélation directe entre l’humidité relative et la concentration des hydrates est compatible avec l’idée que la présence d’un plus grand nombre de molécules d’eau à un RH plus élevé favorise la formation d’hydrates.
Tel que présenté, ce protocole donne une compréhension qualitative des populations de glycine hydratées dans l’atmosphère. En supposant une concentration initiale de glycine isolée de 2,9 millions de molécules par centimètre cube, nous constatons que la glycine non hydratée (n-0) est l’espèce la plus abondante dans la plupart des conditions, à l’exception de T-216.65K et RH-100%. Le dihydrate (n'2), qui a l’énergie sans Gibbs séquentielle la plus faible de l’hydratation à toutes les trois températures, est l’hydrate le plus abondant dans les conditions considérées ici. On prévoit que les hydrates monohydrat (n-1) et les plus grands hydrates (n'3) se trouveraient en quantités négligeables. Lors de l’inspection de la figure 3, l’abondance des grappes n 1-4 peut être liée à la stabilité et à la souche du réseau de liaisons hydrogène des grappes. Ces grappes ont les molécules d’eau de l’hydrogène collé à la mésety d’acide carboxylique de la glycine dans une géométrie ressemblant étroitement à celles de diverses structures d’anneau à hydrogène, ce qui les rend particulièrement stables.
Figure 1 : Description schématique de la procédure actuelle. Un grand bassin de structures de devinettes générées par l’algorithme génétique (GA) est affiné par une série d’optimisations de géométrie PW91 jusqu’à ce qu’un ensemble de structures convergentes soient obtenus. Les fréquences vibratoires de ces structures sont calculées et utilisées pour calculer l’énergie libre Gibbs de formation, qui est à son tour utilisée pour calculer les concentrations d’équilibre des amas dans des conditions ambiantes. S’il vous plaît cliquez ici pour voir une version plus grande de ce chiffre.
Figure 2 : Structure d’annuaire représentative pour chaque cluster. Les scripts internes inclus dans ce protocole nécessitent la structure du répertoire indiqué ci-dessus, où n est le nombre de molécules d’eau. Pour chaque n en gly-h2o-n, il y a les sous-directeurs suivants : GA pour algorithme génétique avec un répertoire GA/pm7, QM pour la mécanique quantique avec QM/pw91-sb pour PW91/6-31-GMD, QM/pw91-lb pour PW91/6-311 GMD, et QM/pw91-lb/ultrafine pour les optimisations et les calculs vibratoires finaux sur les grilles d’intégration ultrafines. S’il vous plaît cliquez ici pour voir une version plus grande de ce chiffre.
Figure 3 : Structures représentatives à faible consommation d’énergie de Gly(H2O)n'0-5. Ces grappes étaient les minima mondiaux d’énergie électronique optimisés au niveau de théorie PW91/6-311-G. S’il vous plaît cliquez ici pour voir une version plus grande de ce chiffre.
Figure 4 : Dépendance à la température de Gly (H2O)n-0-5 comme humidité relative de 100 % et 1 pression au guichet. La concentration des hydrates est donnée dans des unités de molécules cm-3. S’il vous plaît cliquez ici pour voir une version plus grande de ce chiffre.
Figure 5 : Dépendance relative à l’humidité de Gly (H2O)n-0-5 comme 298,15 K et 1 pression au guichet. La concentration des hydrates est donnée dans des unités de molécules cm-3. S’il vous plaît cliquez ici pour voir une version plus grande de ce chiffre.
E[PW91/6-311-G] | 216,65 K | 273,15 K | 298,15 K | ||||||||
LB-UF | ZPVE (ZPVE) | H | S | G G | H | S | G G | H | S | G G | |
Eau | -76.430500 | 13.04 | 1.72 | 42.59 | 5.54 | 2.17 | 44.44 | 3.08 | 2.37 | 45.14 | 1.96 |
Glycine | -284.434838 | 48.55 | 2.65 | 69.53 | 36.14 | 3.70 | 73.81 | 32.09 | 4.22 | 75.61 | 30.22 |
Tableau 1 : Énergies monomer. Les énergies électroniques sont dans les unités de Hartree tandis que toutes les autres quantités sont dans des unités de mol kcal-1. L’eau et la glycine ont été optimisées au niveau de théorie et de fréquences vibratoires PW91/6-311-GMD. Les corrections thermodynamiques pour une pression de 1 atm et la température de 298,15 K ont été calculées à l’aide du script thermo.pl.
E[PW91/6-311-G] | 0 K | 216,65 K | 273,15 K | 298,15 K | ||||||||
¡n | Nom | LB-UF | ZPVE (ZPVE) | H | S | G G | H | S | G G | H | S | G G |
1 | gly-h2o-1 | -360.88481 | 63.96 | 3.61 | 80.12 | 50.22 | 5.12 | 86.27 | 45.52 | 5.85 | 88.83 | 43.33 |
2 | gly-h2o-2 | -437.33763 | 79.33 | 4.53 | 90.86 | 64.17 | 6.46 | 98.78 | 58.81 | 7.40 | 102.06 | 56.30 |
3 | gly-h2o-3 | -513.78620 | 94.52 | 5.67 | 105.08 | 77.42 | 8.08 | 114.94 | 71.19 | 9.23 | 119.00 | 68.27 |
4 | gly-h2o-4 | -590.23667 | 109.80 | 6.03 | 104.98 | 91.30 | 8.78 | 116.21 | 84.40 | 10.11 | 120.87 | 81.14 |
5 | gly-h2o-5 | -666.68845 | 125.80 | 7.26 | 121.70 | 106.69 | 10.47 | 134.83 | 99.44 | 12.01 | 140.24 | 96.00 |
Tableau 2 : Énergies de cluster. Les énergies des structures Gly (H2O)n-1-5 à faible énergie trouvées à l’aide de notre procédure décrites dans la figure 1. Les énergies électroniques sont dans les unités de Hartree tandis que toutes les autres quantités sont dans des unités de mol kcal-1.
Hydratation totale: Gly nH2O 'lt;-gt; Gly (H2O)n | Hydratation séquentielle: Gly (H2O)n-1 - H2O 'lt;-gt; Gly (H2O)n | ||||||||||||||||
E[PW91/6-311-G] | 216.65 | 273.15 | 298.15 | 216.65 | 273.15 | 298.15 | |||||||||||
¡n | nom du système | LB-UF | E (0) | H (T) | G (T) | H (T) | G (T) | H (T) | G (T) | LB-UF | E (0) | H (T) | G (T) | H(T) | G (T) | H (T) | G (T) |
1 | gly-h2o-1 | -12.22 | -9.85 | -10.61 | -3.68 | -10.61 | -1.87 | -10.59 | -1.07 | -12.22 | -9.85 | -10.61 | -3.68 | -10.61 | -1.87 | -10.59 | -1.07 |
2 | gly-h2o-2 | -26.22 | -21.53 | -23.10 | -9.27 | -23.11 | -5.66 | -23.09 | -4.06 | -14.00 | -11.68 | -12.49 | -5.59 | -12.50 | -3.79 | -12.50 | -2.99 |
3 | gly-h2o-3 | -37.56 | -30.72 | -32.88 | -12.90 | -32.87 | -7.69 | -32.82 | -5.38 | -11.34 | -9.19 | -9.78 | -3.63 | -9.76 | -2.03 | -9.73 | -1.32 |
4 | gly-h2o-4 | -50.10 | -40.34 | -43.48 | -15.87 | -43.54 | -8.71 | -43.51 | -5.55 | -12.54 | -9.62 | -10.60 | -2.97 | -10.67 | -1.02 | -10.69 | -0.17 |
5 | gly-h2o-5 | -63.45 | -51.41 | -55.42 | -20.58 | -55.51 | -11.48 | -55.48 | -7.45 | -13.35 | -11.07 | -11.94 | -4.71 | -11.97 | -2.77 | -11.97 | -1.90 |
Tableau 3 : Énergies d’hydratation. L’énergie totale de l’hydratation et de l’énergie de l’hydratation séquentielle pour Gly (H2O)n 1-5 dans les unités de kcal mol-1. Ici, E[PW91/6-311-G] est le changement dans l’énergie électronique, E(0) est l’énergie vibratoire zéro point (ZPVE) changement corrigé de l’énergie, H(T) est le changement enthalpy à la température T, et G(T) est le changement d’énergie libre Gibbs de l’hydratation de chaque Gly (H2O)n '1-5 cluster.
Répartition des hydrates d’équilibre en fonction de la température et de l’humidité relative | |||||||||
T-298.15K | T-273.15K | T-216.65K | |||||||
Gly (H2O)n | RH 100% | RH-50% | RH-20% | RH 100% | RH-50% | RH-20% | RH 100% | RH-50% | RH-20% |
0 | 1.3E 06 | 2.2E-06 | 2.7E 06 | 1.1E-06 | 2.0E 06 | 2.7E 06 | 6.1E-05 | 1.5E 06 | 2.5E 06 |
1 | 2.3E 05 | 1.9E 05 | 9.5E-04 | 2.0E 05 | 1.9E 05 | 9.9E-04 | 1.2E 05 | 1.5E 05 | 9.5E-04 |
2 | 1.0E 06 | 4.3E 05 | 8.4E 04 | 1.3E 06 | 6.1E-05 | 1.3E 05 | 1.8E 06 | 1.1E-06 | 3.0E 05 |
3 | 2.8E 05 | 5.8E 04 | 4.5E 03 | 3.2E 05 | 7.4E 04 | 6.3E 03 | 3.1E-05 | 9.6E-04 | 1.0E 04 |
4 | 1.1E-04 | 1.1E 03 | 3.4E 01 | 1.3E 04 | 1.5E 03 | 5.0E 01 | 1.1E-04 | 1.8E 03 | 7.5E 01 |
5 | 7.5E 03 | 3.9E 02 | 4.9E 00 | 1.2E 04 | 7.2E 02 | 9.7E-00 | 2.4E 04 | 1.9E 03 | 3.1E-01 |
Tableau 4 : Les concentrations d’hydrates d’équilibre de Gly (H2O)n-0-5 en fonction de la température de la fonction (T-298.15K, 273.15K, 216.65K) et l’humidité relative (RH-100%, 50%, 20%). La concentration des hydrates est donnée dans des unités de molécules cm-3 en supposant des valeurs expérimentales56,57,58, de [Gly]0 - 2,9 x 106 cm-3 et [H2O] - 7,7 x 101017 cm-3, 1,6 x 1017 cm-3 et 9,9 x 1014 cm-3 à 100% humidité relative et T 298,15 K, 273,15 K, et 216,65 K, respectivement59.
Fichiers supplémentaires. S’il vous plaît cliquez ici pour télécharger ces fichiers.
L’exactitude des données générées par ce protocole dépend principalement de trois choses : (i) la variété des configurations échantillonnées par l’étape 2, (ii) la précision de la structure électronique du système, (iii) et la précision des corrections thermodynamiques. Chacun de ces facteurs peut être abordé en modifiant la méthode en éditant les scripts inclus. Le premier facteur est facilement surmonté avec l’utilisation d’un plus grand bassin initial de structures générées au hasard, des itérations plus nombreuses de l’AG, et une définition plus souple des critères impliqués dans l’AG. En outre, on peut utiliser une méthode semi-empirique différente comme le modèle de fixation serrée-fonctionnelle de densité de charge auto-cohérente (SCC-DFTB)62 modèle et le potentiel de fragment efficace (EFP)63 modèle afin d’explorer les effets de différentes descriptions physiques. La principale limitation ici est l’incapacité de la méthode à former ou briser les liaisons covalentes, ce qui signifie que les monomères sont gelés. La procédure d’AG ne trouve que les positions relatives les plus stables de ces monomères congelés selon la description semi-empirique.
L’exactitude de la structure électronique du système peut être améliorée de diverses façons, chacune avec son coût de calcul. On peut choisir une meilleure densité fonctionnelle, comme M06-2X64 et wB97X-V65, ou mécanique quantique (QM) méthode comme le M-ller-Plesset66,67,68 (MPn) théories de perturbation et couplé-cluster69 (CC) méthodes afin d’améliorer la description physique du système. Dans la hiérarchie des fonctions, les performances s’améliorent généralement en passant de fonctions d’approximation généralisée-gradient (GGA) comme PW91 aux fonctions hybrides séparées à la portée comme wB97X-D et méta-GGA fonctions hybrides comme M06-2X.
L’inconvénient des méthodes DFT est qu’une convergence systématique vers une valeur précise n’est pas possible; cependant, les méthodes DFT sont computationnellement peu coûteuses et il y a une grande variété de fonctions pour une grande variété d’applications.
Energies calculées à l’aide de méthodes de fonction d’onde comme MP2 et CCSD(T) en conjonction avec des ensembles de base cohérents de corrélation de nombre cardinal croissant ([aug-]cc-pV[D,T,Q,...] Z) convergent vers leur limite complète de base définie systématiquement, mais le coût de calcul de chaque calcul devient prohibitif à mesure que la taille du système augmente. Un plus grand raffinement de la structure électronique peut être réalisé en utilisant des ensembles de base explicitement corrélés70 et en extrapolant à la limite de base complète (CBS)71. Nos travaux récents suggèrent qu’une approche perturbative de second ordre explicitement corrélée et en corrélation explicite de densité M’ller-Plesset (DF-MP2-F12) donne des énergies approchant celle des calculs MP2/CBS32. La modification du protocole actuel pour utiliser différentes méthodes de structure électronique comporte deux étapes : (i) préparer un fichier d’entrée de modèle suivant la syntaxe donnée par le logiciel, (ii) et modifier les scripts run-pw91-sb.csh, run-pw91-lb.csh, et exécuter-pw91-lb-ultrafine.csh scripts pour générer la syntaxe fichier d’entrée correcte ainsi que le script de soumettre correctement pour le logiciel.
Enfin, l’exactitude des corrections thermodynamiques dépend de la méthode de la structure électronique ainsi que de la description du PSE autour du minimum global. Une description précise du PSE exige le calcul des dérivés de troisième et plus haut débit du PSE en ce qui concerne les déplacements dans les degrés nucléaires de liberté, tels que le champ de la force quartique72,73 (QFF), ce qui est une tâche exceptionnellement coûteuse. Le protocole actuel utilise l’approximation harmonique des oscillateurs aux fréquences vibratoires, ce qui entraîne la nécessité de calculer seulement jusqu’à la deuxième dérivé du PSE. Cette approche devient problématique dans les systèmes à forte harmonie, tels que les molécules très disquettes et les potentiels symétriques à double puits en raison de la grande différence dans le vrai PSE et le PSE harmonique. En outre, le coût d’avoir un PSE de haute qualité à partir d’une méthode de structure électronique exigeante en calcul ne fait qu’aggraver le problème du coût des calculs de fréquence vibratoire. Une approche pour surmonter cela est d’utiliser les énergies électroniques à partir d’un calcul de structure électronique de haute qualité avec des fréquences vibratoires calculées sur un PSE de qualité inférieure, résultant en un équilibre entre le coût et la précision. Le protocole actuel peut être modifié pour utiliser différentes descriptions de PSE telles que décrites dans le paragraphe précédent; cependant, on peut également modifier les mots-clés de fréquence vibratoire dans les scripts et les modèles pour calculer les fréquences vibratoires anharmoniques.
Deux questions cruciales pour tout protocole d’échantillonnage de configuration sont la méthode initiale pour l’échantillonnage de la surface énergétique potentielle et les critères utilisés pour identifier chaque cluster. Nous avons fait un usage intensif d’une variété de méthodes dans nos travaux précédents. Pour le premier numéro, la méthode initiale pour l’échantillonnage de la surface énergétique potentielle, nous avons fait le choix d’utiliser l’AG avec des méthodes semi-empiriques basées sur ces facteurs. Échantillonnage configurationnel à l’aide de l’intuition chimique26, échantillonnage aléatoire, et la dynamique moléculaire (MD)29,30, ne parviennent pas à trouver des minima mondiaux putatifs régulièrement pour les grappes de plus de 10 monomers, comme nous l’avons observé dans nos études sur les grappes d’eau18. Nous avons utilisé avec succès le saut de bassin (BH) pour étudier le PES complexe de (H2O)1174, mais il a exigé l’inclusion manuelle de certains isomers potentiels à faible énergie l’algorithme BH n’a pas trouvé. Une comparaison de la performance de BH et GA dans la recherche du minimum mondial de grappes d’eau, (H2O)n '10-20 a démontré que l’AG a toujours trouvé le minimum mondial plus rapide que BH75. GA tel qu’il est mis en œuvre dans OGOLEM et CLUSTER est très polyvalent parce qu’il peut être appliqué à n’importe quel cluster moléculaire et il peut l’interface avec un grand nombre de paquets avec champ de force classique, semi-empirique, densité fonctionnelle, et ab initio capacités. Le choix du PM7 est motivé par sa vitesse et sa précision raisonnable. Pratiquement n’importe quelle autre méthode semi-empirique aurait un coût de calcul significativement plus élevé.
En ce qui concerne le deuxième numéro, nous avons exploré en utilisant différents critères pour identifier des structures uniques allant des énergies électroniques, des moments dipole, des RMSD qui se chevauchent et des constantes de rotation. L’utilisation de moments dipole s’est avérée difficile parce que les composants du moment dipole dépendaient de l’orientation de la molécule et le moment total de dipole était très sensible aux différences de géométrie de telle sorte qu’il était difficile de fixer des seuils déterminant est que les structures sont les mêmes ou uniques. Une combinaison d’énergies électroniques et de constantes de rotation s’est avérée la plus utile.
Les critères actuels pour juger deux structures uniques sont basés sur un seuil de différence d’énergie de 0,10 kcal mol-1 et une différence constante de rotation de 1%. Par conséquent, deux structures sont considérées comme différentes si leurs énergies diffèrent de plus de 0,10 kcal mol-1 (0,00015 a.u.) ET l’une de leurs trois constantes de rotation (A, B, C) diffèrent de plus de 1%. D’importants repères internes au fil des ans ont conclu que ces seuils étaient des choix raisonnables. Notre approche d’échantillonnage configurationnel et notre méthodologie de dépistage ont été appliquées à des grappes très faiblement liées telles que les hydrocarbures polyaromatiques complexes avec de l’eau76,77 ainsi que des hydrates de sulfate sternary fortement liés contenant de l’ammoniac et des amines32. Pour les grappes où il ya différents états de protonation à considérer, la meilleure approche est d’exécuter divers calculs d’AG, chacun commençant par des monomères dans différents états de protonation. Cela garantit que les structures avec différents états de protonation sont soigneusement examinées. Cependant, les calculs DFT de bas niveau permettent souvent aux états de protonation de changer au cours de l’optimisation de la géométrie, produisant ainsi l’état de protonation le plus stable indépendamment de la géométrie de départ.
Nos méthodes d’échantillonnage configurationnelle GA devraient bien fonctionner même pour les molécules disquettes tant que les codes GA sont interatronisés avec des méthodes générales non paramétrisées qui permettent aux monomères d’adopter différentes configurations au cours de la course GA. Par exemple, l’interaction de l’AG avec le PM7 permettrait aux structures des monomères de changer, mais si leurs liens se brisent comme cela se produirait lorsque les États protonation changent, les structures peuvent être rejetées en tant que candidats inacceptables.
Nous avons examiné différentes façons de corriger les lacunes de l’approximation harmonique, en particulier celles résultant de basses fréquences vibratoires. Il n’est pas difficile d’intégrer l’approximation quasi harmonique dans la méthodologie actuelle. Cependant, il ya encore des questions sur la méthode quasi-harmonique, en particulier quand il s’agit de la fréquence de coupure ci-dessous laquelle il sera appliqué. En outre, il n’y a pas d’ouvrages de benchmarking rigoureux examinant la fiabilité de l’approximation quasi-RRHO, même si la sagesse conventionnelle suggère qu’il devrait être une amélioration par rapport à l’approximation RRHO.
Le protocole ainsi présenté peut être généralisé à n’importe quel système de clusters moléculaires de phase de gaz non liés. Il peut également être généralisé d’utiliser n’importe quelle méthode semi-empirique, méthode de structure électronique et logiciel, et méthode d’analyse vibratoire et logiciel en éditant les scripts et les modèles. Cela suppose que l’utilisateur est à l’aise avec l’interface Linux de ligne de commande, le script Python, et l’informatique haute performance. La syntaxe et l’apparence inconnues du système d’exploitation Linux et le manque d’expérience de script est le plus grand piège dans ce protocole et est l’endroit où les nouveaux étudiants luttent le plus. Ce protocole a été utilisé avec succès dans une variété d’implémentations pendant des années dans notre groupe, principalement en se concentrant sur les effets de l’acide sulfurique et l’ammoniac sur la formation d’aérosols. D’autres améliorations à ce protocole impliqueront une interface plus robuste à plus de logiciels de structure électronique, des implémentations alternatives de l’algorithme génétique, et peut-être l’utilisation de nouvelles méthodes pour des calculs plus rapides des énergies électroniques et vibratoires. Nos applications actuelles de ce protocole explorent l’importance des acides aminés dans les premiers stades de la formation d’aérosols dans l’atmosphère actuelle et dans la formation de molécules biologiques plus grandes dans des environnements prébiotiques.
Aucun.
Ce projet a été soutenu par des subventions CHE-1229354, CHE-1662030, CHE-1721511 et CHE-1903871 de la National Science Foundation (GCS), l’Arnold and Mabel Beckman Foundation Beckman Scholar Award (AGG), et la bourse Barry M. Goldwater (AGG). Des ressources informatiques haute performance du Consortium MERCURY (http://www.mercuryconsortium.org) ont été utilisées.
Name | Company | Catalog Number | Comments |
Avogadro | https://avogadro.cc | Open-source molecular visualization program | |
Gaussian [09/16] Software | http://www.gaussian.com/ | Commercial ab initio electronic structure program | |
MOPAC 2016 | http://openmopac.net/MOPAC2016.html | Open-source semi-empirical program | |
OGOLEM Software | https://www.ogolem.org | Genetic algorithm-based global optimization program | |
OpenBabel | http://openbabel.org/wiki/Main_Page | Open-source cheminformatics library | |
calcRotConsts.py | Shields Group, Department of Chemistry, Furman University | Python script to compute rotational constants | |
calcSymmetry.csh | Shields Group, Department of Chemistry, Furman University | Shell script to calculate symmetry number of a molecule given Cartesian coordinates | |
combine-GA.csh | Shields Group, Department of Chemistry, Furman University | Shell script to combine energy and rotational constants from different GA directories | |
combine-QM.csh | Shields Group, Department of Chemistry, Furman University | Shell script to combine energy and rotational constants from different QM directories | |
gaussianE.csh | Shields Group, Department of Chemistry, Furman University | Shell script to extract Gaussian 09 energies | |
gaussianFreqs.csh | Shields Group, Department of Chemistry, Furman University | Shell script to extract Gaussian 09 vibrational frequencies | |
getrotconsts | Shields Group, Department of Chemistry, Furman University | Executable to calculate rotational constants given a molecule's Cartesian coordinates | |
getRotConsts-dft-lb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of large basis DFT optimized structures | |
getRotConsts-dft-lb-ultrafine.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of ultrafine DFT optimized structures | |
getRotConsts-dft-sb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of small basis DFT optimized structures | |
getRotConsts-GA.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute rotational constants for a batch of genetic algorithm optimized structures | |
global-minimum-coords.xyz | Shields Group, Department of Chemistry, Furman University | Cartesian coordinates of global minimum structures of gly-(h2o)n, where n=0-5 | |
make-thermo-gaussian.csh | Shields Group, Department of Chemistry, Furman University | Shell script to extract data from Gaussian output files and make input files for the thermo.pl script | |
ogolem-input-file.ogo | Shields Group, Department of Chemistry, Furman University | Ogolem sample input file | |
ogolem-submit-script.pbs | Shields Group, Department of Chemistry, Furman University | PBS batch submission file for Ogolem calculations | |
README.docx | Shields Group, Department of Chemistry, Furman University | Clarifications to help readers use the scripts effectively | |
runogolem.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run OGOLEM | |
run-pw91-lb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run a batch of large basis DFT optimization calculations | |
run-pw91-lb-ultrafine.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run a batch of ultrafine DFT optimization calculations | |
run-pw91-sb.csh | Shields Group, Department of Chemistry, Furman University | Shell script to run a batch of small basis DFT optimization calculations | |
run-thermo-pw91.csh | Shields Group, Department of Chemistry, Furman University | Shell script to compute the thermodynamic corrections for a batch of DFT optimized structures | |
similarityAnalysis.py | Shields Group, Department of Chemistry, Furman University | Python script to determine unique structures based on rotational constants and energies | |
symmetry | Shields Group, Department of Chemistry, Furman University | Executable to calculate molecular symmetry given Cartesian coordinates | |
symmetry.c | (C) 1996, 2003 S. Patchkovskii, Serguei.Patchkovskii@sympatico.ca | C code to determine the molecular symmstry of a molecule given Cartesian coordinates | |
template-marcy.pbs | Shields Group, Department of Chemistry, Furman University | Template for a PBS submit script which uses OGOLEM | |
template-pw91.com | Shields Group, Department of Chemistry, Furman University | Template Gaussian 09 input | |
template-pw91-HL.com | Shields Group, Department of Chemistry, Furman University | Template Gaussian 09 input for ultrafine DFT optimization | |
thermo.pl | https://www.nist.gov/mml/csd/chemical-informatics-research-group/products-and-services/program-computing-ideal-gas | Perl open-source script to compute ideal gas thermodynamic corrections | |
gly-h2o-n.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet for the complete protocol | |
table-1.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
table-2.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
table-3.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
table-4.xlsx | Shields Group, Department of Chemistry, Furman University | Excel spreadsheet | |
water.xyz | Shields Group, Department of Chemistry, Furman University | Cartesian coordinates of water | |
glycine.xyz | Shields Group, Department of Chemistry, Furman University | Cartesian coordinates of glycine |
Demande d’autorisation pour utiliser le texte ou les figures de cet article JoVE
Demande d’autorisationThis article has been published
Video Coming Soon