Contrainte (chimie computationnelle)

En chimie computationnelle , un algorithme de contrainte est une méthode permettant de satisfaire le mouvement newtonien d'un corps rigide constitué de points de masse. Un algorithme de contrainte est utilisé pour garantir que la distance entre les points de masse est maintenue. Les étapes générales impliquées sont les suivantes : (i) choisir de nouvelles coordonnées non contraintes (coordonnées internes), (ii) introduire des forces de contrainte explicites, (iii) minimiser les forces de contrainte implicitement par la technique des multiplicateurs de Lagrange ou des méthodes de projection.

Les algorithmes de contraintes sont souvent appliqués aux simulations de dynamique moléculaire . Bien que ces simulations soient parfois effectuées à l'aide de coordonnées internes qui satisfont automatiquement aux contraintes de longueur de liaison, d'angle de liaison et d'angle de torsion, les simulations peuvent également être effectuées à l'aide de forces de contrainte explicites ou implicites pour ces trois contraintes. Cependant, les forces de contrainte explicites donnent lieu à une inefficacité ; une puissance de calcul plus importante est nécessaire pour obtenir une trajectoire d'une longueur donnée. Par conséquent, les coordonnées internes et les solveurs de contraintes à force implicite sont généralement préférés.

Les algorithmes de contraintes permettent d'obtenir une efficacité de calcul en négligeant le mouvement le long de certains degrés de liberté. Par exemple, dans la dynamique moléculaire atomistique, la longueur des liaisons covalentes avec l'hydrogène est généralement limitée ; cependant, les algorithmes de contraintes ne doivent pas être utilisés si les vibrations le long de ces degrés de liberté sont importantes pour le phénomène étudié.

Le mouvement d'un ensemble de N particules peut être décrit par un ensemble d'équations différentielles ordinaires du second ordre, la deuxième loi de Newton, qui peut être écrite sous forme matricielle

où M est une matrice de masse et q est le vecteur de coordonnées généralisées qui décrivent les positions des particules. Par exemple, le vecteur q peut être une matrice cartésienne 3N des positions des particules r _k , où k va de 1 à N ; en l'absence de contraintes, M serait la matrice carrée diagonale 3N x 3N des masses des particules. Le vecteur f représente les forces généralisées et le scalaire V ( q ) représente l'énergie potentielle, tous deux étant des fonctions des coordonnées généralisées q .

Si M contraintes sont présentes, les coordonnées doivent également satisfaire M équations algébriques indépendantes du temps

où l'indice j va de 1 à M. Par souci de concision, ces fonctions g _i sont regroupées dans un vecteur g de dimension M ci-dessous. La tâche consiste à résoudre l'ensemble combiné d'équations algébriques différentielles (DAE), au lieu des simples équations différentielles ordinaires (ODE) de la deuxième loi de Newton.

Ce problème a été étudié en détail par Joseph Louis Lagrange , qui a exposé la plupart des méthodes pour le résoudre. L'approche la plus simple consiste à définir de nouvelles coordonnées généralisées qui ne sont pas contraintes ; cette approche élimine les équations algébriques et réduit le problème à la résolution d'une équation différentielle ordinaire. Une telle approche est utilisée, par exemple, pour décrire le mouvement d'un corps rigide ; la position et l'orientation d'un corps rigide peuvent être décrites par six coordonnées indépendantes et non contraintes, plutôt que de décrire les positions des particules qui le composent et les contraintes entre elles qui maintiennent leurs distances relatives. L'inconvénient de cette approche est que les équations peuvent devenir difficiles à manier et complexes ; par exemple, la matrice de masse M peut devenir non diagonale et dépendre des coordonnées généralisées.

Une deuxième approche consiste à introduire des forces explicites qui agissent pour maintenir la contrainte ; par exemple, on pourrait introduire des forces de ressort fortes qui imposent des distances entre les points de masse au sein d'un corps « rigide ». Les deux difficultés de cette approche sont que les contraintes ne sont pas satisfaites exactement et que les forces fortes peuvent nécessiter des pas de temps très courts, ce qui rend les simulations inefficaces sur le plan informatique.

Une troisième approche consiste à utiliser une méthode telle que les multiplicateurs de Lagrange ou la projection sur la variété de contraintes pour déterminer les ajustements de coordonnées nécessaires pour satisfaire les contraintes.

Enfin, il existe diverses approches hybrides dans lesquelles différents ensembles de contraintes sont satisfaits par différentes méthodes, par exemple, des coordonnées internes, des forces explicites et des solutions de force implicite.

L'approche la plus simple pour satisfaire aux contraintes de minimisation d'énergie et de dynamique moléculaire consiste à représenter le système mécanique dans des coordonnées internes correspondant à des degrés de liberté indépendants et non contraints du système. Par exemple, les angles dièdres d'une protéine sont un ensemble indépendant de coordonnées qui spécifient les positions de tous les atomes sans nécessiter de contraintes. La difficulté de telles approches à coordonnées internes est double : les équations newtoniennes du mouvement deviennent beaucoup plus complexes et les coordonnées internes peuvent être difficiles à définir pour des systèmes cycliques de contraintes, par exemple dans le cas du plissement d'anneaux ou lorsqu'une protéine a une liaison disulfure.

Les méthodes originales de minimisation efficace de l'énergie récursive dans les coordonnées internes ont été développées par Gō et ses collègues.

Des solveurs efficaces de contraintes récursives à coordonnées internes ont été étendus à la dynamique moléculaire. Des méthodes analogues ont été appliquées ultérieurement à d'autres systèmes.

Dans la plupart des simulations de dynamique moléculaire qui utilisent des algorithmes de contraintes, les contraintes sont appliquées à l'aide de la méthode des multiplicateurs de Lagrange. Étant donné un ensemble de n contraintes linéaires ( holonomes ) à l'instant t ,

où et sont les positions des deux particules impliquées dans la k -ième contrainte à l'instant t et est la distance inter-particules prescrite.

Les forces dues à ces contraintes s'ajoutent dans les équations du mouvement, ce qui donne, pour chacune des N particules du système

L'ajout des forces de contrainte ne modifie pas l'énergie totale, car le travail net effectué par les forces de contrainte (pris sur l'ensemble des particules sur lesquelles agissent les contraintes) est nul. Notez que le signe est arbitraire et que certaines références ont un signe opposé.

En intégrant les deux côtés de l'équation par rapport au temps, les coordonnées contraintes des particules à l'instant, , sont données,

où est la position non contrainte (ou non corrigée) de la i -ème particule après intégration des équations de mouvement non contraintes.

Pour satisfaire les contraintes du pas de temps suivant, les multiplicateurs de Lagrange doivent être déterminés selon l'équation suivante,

Cela implique de résoudre un système d' équations non linéaires

simultanément pour les multiplicateurs de Lagrange inconnus .

Ce système d' équations non linéaires à inconnues est généralement résolu à l'aide de la méthode de Newton-Raphson où le vecteur de solution est mis à jour à l'aide de

où est le jacobien des équations σ _k :

Étant donné que toutes les particules ne contribuent pas à toutes les contraintes, il s'agit d'une matrice en blocs et elle peut être résolue individuellement pour chaque bloc-unité de la matrice. En d'autres termes, elle peut être résolue individuellement pour chaque molécule.

Au lieu de mettre à jour constamment le vecteur , l'itération peut être démarrée avec , ce qui donne des expressions plus simples pour et . Dans ce cas

puis est mis à jour vers

Après chaque itération, les positions des particules non contraintes sont mises à jour à l'aide

Le vecteur est ensuite réinitialisé à

La procédure ci-dessus est répétée jusqu'à ce que la solution des équations de contrainte, , converge vers une tolérance prescrite d'une erreur numérique.

Bien qu'il existe un certain nombre d'algorithmes pour calculer les multiplicateurs de Lagrange, ces différences ne reposent que sur les méthodes de résolution du système d'équations. Pour ces méthodes, les méthodes quasi-Newton sont couramment utilisées.

L'algorithme SETTLE résout analytiquement le système d'équations non linéaires pour des contraintes en temps constant. Bien qu'il ne soit pas évolutif vers un plus grand nombre de contraintes, il est très souvent utilisé pour contraindre les molécules d'eau rigides, qui sont présentes dans presque toutes les simulations biologiques et sont généralement modélisées à l'aide de trois contraintes (par exemple les modèles d'eau SPC/E et TIP3P ).

L'algorithme SHAKE a été développé à l'origine pour satisfaire une contrainte de géométrie de liaison lors de simulations de dynamique moléculaire. La méthode a ensuite été généralisée pour gérer toute contrainte holonome, telle que celles requises pour maintenir des angles de liaison constants ou une rigidité moléculaire.

Dans l'algorithme SHAKE, le système d'équations de contraintes non linéaires est résolu à l'aide de la méthode de Gauss-Seidel qui approxime la solution du système d'équations linéaires à l'aide de la méthode de Newton-Raphson ;

Cela revient à supposer que est diagonalement dominant et à résoudre la ième équation uniquement pour l' inconnue. En pratique, on calcule

pour tous de manière itérative jusqu'à ce que les équations de contrainte soient résolues avec une tolérance donnée.

Le coût de calcul de chaque itération est de , et les itérations elles-mêmes convergent linéairement.

Une forme non itérative de SHAKE a été développée plus tard.

Il existe plusieurs variantes de l'algorithme SHAKE. Bien qu'elles diffèrent dans la manière dont elles calculent ou appliquent les contraintes elles-mêmes, les contraintes sont toujours modélisées à l'aide de multiplicateurs de Lagrange qui sont calculés à l'aide de la méthode de Gauss-Seidel.

L'algorithme SHAKE original est capable de contraindre à la fois les molécules rigides et flexibles (par exemple l'eau, le benzène et le biphényle ) et introduit une erreur négligeable ou une dérive énergétique dans une simulation de dynamique moléculaire. Un problème avec SHAKE est que le nombre d'itérations nécessaires pour atteindre un certain niveau de convergence augmente à mesure que la géométrie moléculaire devient plus complexe. Pour atteindre une précision informatique de 64 bits (une tolérance relative de ) dans une simulation de dynamique moléculaire typique à une température de 310 K, un modèle d'eau à 3 sites ayant 3 contraintes pour maintenir la géométrie moléculaire nécessite en moyenne 9 itérations (soit 3 par site par pas de temps). Un modèle de butane à 4 sites avec 5 contraintes nécessite 17 itérations (22 par site), un modèle de benzène à 6 sites avec 12 contraintes nécessite 36 itérations (72 par site), tandis qu'un modèle de biphényle à 12 sites avec 29 contraintes nécessite 92 itérations (229 par site par pas de temps). Par conséquent, les exigences en termes de CPU de l’algorithme SHAKE peuvent devenir importantes, en particulier si un modèle moléculaire présente un degré élevé de rigidité.

Une extension ultérieure de la méthode, QSHAKE ( Quaternion SHAKE), a été développée comme une alternative plus rapide pour les molécules composées d'unités rigides, mais elle n'est pas aussi polyvalente. Elle fonctionne de manière satisfaisante pour les boucles rigides telles que les systèmes cycliques aromatiques , mais QSHAKE échoue pour les boucles flexibles, comme lorsqu'une protéine a une liaison disulfure.

D'autres extensions incluent RATTLE, WIGGLE, et MSHAKE.

Alors que RATTLE fonctionne de la même manière que SHAKE, mais en utilisant le schéma d'intégration temporelle Velocity Verlet , WIGGLE étend SHAKE et RATTLE en utilisant une estimation initiale des multiplicateurs de Lagrange basée sur les vitesses des particules. Il convient de mentionner que MSHAKE calcule les corrections sur les forces de contrainte , obtenant ainsi une meilleure convergence.

Une dernière modification de l'algorithme SHAKE est l'algorithme P-SHAKE qui s'applique aux molécules très rigides ou semi-rigides . P-SHAKE calcule et met à jour un préconditionneur qui est appliqué aux gradients de contrainte avant l'itération SHAKE, ce qui fait que le Jacobien devient diagonal ou fortement diagonalement dominant. Les contraintes ainsi découplées convergent beaucoup plus rapidement (de manière quadratique plutôt que linéaire) pour un coût de .

L'algorithme M-SHAKE résout le système d'équations non linéaires en utilisant directement la méthode de Newton . À chaque itération, le système d'équations linéaires

est résolu exactement en utilisant une décomposition LU . Chaque itération coûte des opérations, mais la solution converge de manière quadratique , nécessitant moins d'itérations que SHAKE.

Cette solution a été proposée pour la première fois en 1986 par Ciccotti et Ryckaert sous le titre « méthode matricielle », mais elle diffère dans la résolution du système d'équations linéaires. Ciccotti et Ryckaert proposent d'inverser directement la matrice, mais en ne le faisant qu'une seule fois, lors de la première itération. La première itération coûte alors des opérations, alors que les itérations suivantes ne coûtent que des opérations (pour la multiplication matrice-vecteur). Cette amélioration a cependant un prix, puisque le Jacobien n'est plus mis à jour, la convergence est seulement linéaire , bien qu'à un rythme beaucoup plus rapide que pour l'algorithme SHAKE.

Plusieurs variantes de cette approche basées sur des techniques de matrices clairsemées ont été étudiées par Barth et al. .

L'algorithme SHAPE est un analogue multicentrique de SHAKE pour contraindre les corps rigides de trois centres ou plus. Comme SHAKE, une étape sans contrainte est effectuée puis corrigée en calculant et en appliquant directement la matrice de rotation du corps rigide qui satisfait :

Cette approche implique une diagonalisation de matrice 3×3 unique suivie de trois ou quatre itérations rapides de Newton pour déterminer la matrice de rotation. SHAPE fournit la même trajectoire que celle fournie par SHAKE itératif entièrement convergé, mais il s'avère plus efficace et plus précis que SHAKE lorsqu'il est appliqué à des systèmes impliquant trois centres ou plus. Il étend la capacité des contraintes de type SHAKE aux systèmes linéaires avec trois atomes ou plus, aux systèmes plans avec quatre atomes ou plus et aux structures rigides nettement plus grandes où SHAKE est intraitable. Il permet également de lier des corps rigides à un ou deux centres communs (par exemple des plans peptidiques) en résolvant les contraintes de corps rigides de manière itérative de la même manière de base que SHAKE est utilisé pour les atomes impliquant plus d'une contrainte SHAKE.

Une méthode de contrainte alternative, LINCS (Linear Constraint Solver) a été développée en 1997 par Hess, Bekker, Berendsen et Fraaije, et était basée sur la méthode de 1986 d'Edberg, Evans et Morriss (EEM), et une modification de celle-ci par Baranyai et Evans (BE).

LINCS applique les multiplicateurs de Lagrange aux forces de contrainte et résout les multiplicateurs en utilisant un développement en série pour approximer l'inverse du Jacobien :

à chaque étape de l'itération de Newton. Cette approximation ne fonctionne que pour les matrices dont les valeurs propres sont inférieures à 1, ce qui rend l'algorithme LINCS adapté uniquement aux molécules à faible connectivité.

Il a été rapporté que LINCS est 3 à 4 fois plus rapide que SHAKE.

Des méthodes hybrides ont également été introduites dans lesquelles les contraintes sont divisées en deux groupes ; les contraintes du premier groupe sont résolues à l'aide de coordonnées internes tandis que celles du second groupe sont résolues à l'aide de forces de contrainte, par exemple, par un multiplicateur de Lagrange ou une méthode de projection. Cette approche a été lancée par Lagrange, et aboutit à des équations de Lagrange de type mixte .