Par exemple, l'eau:
il est évident que l'hydrogène et l'oxygène à l'intérieur d'une molécule interagissent d'une manière complètement différente du même oxygène avec l'hydrogène d'une molécule voisine. Ainsi, les interactions intramoléculaires et intermoléculaires sont distinguées. L'interaction intermoléculaire peut être spécifiée par des potentiels à courte portée et des paires de Coulomb, qui ont été discutés dans les articles précédents. Ici, nous allons nous concentrer sur l'intramoléculaire.
Le type d'interaction intramoléculaire le plus courant est celui des liaisons chimiques (valence). Les liaisons chimiques sont établies par la dépendance fonctionnelle de l'énergie potentielle sur la distance entre les atomes liés, c'est-à-dire, en fait, par le même potentiel de paire. Mais, contrairement aux potentiels de paires ordinaires, cette interaction n'est pas spécifiée pour certains types de particules, mais pour une paire spécifique de particules (par leurs indices). Les formes fonctionnelles les plus courantes des potentiels de liaison chimique sont les potentiels harmoniques:
où r est la distance entre les particules, k est la constante de rigidité de la liaison et r 0 est la longueur de la liaison à l'équilibre; et potentiel Morse :
où D est la profondeur du puits de potentiel, le paramètre α caractérise la largeur du puits de potentiel.
Le prochain type d'interactions intramoléculaires est les angles de liaison. Regardons à nouveau l'image du titre. Pourquoi la molécule est-elle représentée avec un coin, parce que les forces de l'électrostatique étaient censées fournir la distance maximale entre les ions hydrogène, ce qui correspond à un angle HOH égal à 180 °? Le fait est que tout n'est pas dessiné dans la figure. Du cours de chimie à l'école, vous pouvez vous rappeler que l'oxygène a 2 paires d'électrons solitaires supplémentaires, l'interaction avec laquelle fausse l'angle:
Dans la dynamique moléculaire classique, des objets tels que des électrons ou des nuages d'électrons ne sont généralement pas introduits, par conséquent, pour simuler des angles «corrects», le potentiel de l'angle de liaison est utilisé, c'est-à-dire dépendance fonctionnelle de l'énergie potentielle sur les coordonnées de 3 particules. L'un des potentiels les plus pratiques est le cosinus harmonique:
où θ est l'angle formé par le triple des particules, k est la constante de rigidité et θ 0 est l'angle d'équilibre.
Il existe des potentiels intramoléculaires d'ordre supérieur, par exemple les angles de torsion , mais ils sont encore plus artificiels que les angles de liaison.
L'ajout d'interactions entre particules avec des indices prédéterminés est trivial. Pour les liens, nous stockons un tableau contenant les indices des particules liées et le type d'interaction. Nous donnons à chaque thread sa propre gamme de liens pour le traitement et vous avez terminé. De même avec les angles de liaison. Par conséquent, nous allons immédiatement compliquer la tâche pour nous-mêmes: nous ajouterons la possibilité de créer / supprimer des liaisons chimiques et des angles de liaison d'exécution. Cela nous fait immédiatement sortir du plan de la dynamique moléculaire classique et nous ouvre un nouvel horizon de possibilités. Sinon, vous pouvez simplement télécharger quelque chose à partir des packages existants, par exemple LAMMPS , DL_POLY ou GROMACS , d'autant plus qu'ils sont distribués gratuitement.
Maintenant pour un peu de code. Ajoutons les champs appropriés à la structure principale:
//bonds:
int nBond; // number of bonds
int mxBond; // maximal number of bonds
int4* bonds; // array of bonds
int* nbonds; // count of bond for a given atom
int* neighToBind; // a neighbor of a given atom for binding
int* canBind; // flags that atom[iat] can be bind
int* r2Min; // distances for the nearest neighbor (used for binding)
int* parents; // indexes of one of the atom bonded with a given
cudaBond* bondTypes;
int** def_bonds; // array[nSpec][nSpec] of default bond types
int** bindBonds; // array[nSpec][nSpec] bond types created by binding
float** bindR2; // square of binding distance [nSpec][nSpec]
//angles:
int nAngle; // number of angles
int mxAngle;
int4* angles; // array of angles
int* nangles; // number of angles for given atom
int* oldTypes;
cudaAngle* angleTypes;
int* specAngles; // [nSp] angle type formed by given species
Le nombre de liens et d'angles est variable, mais presque toujours vous pouvez estimer le maximum possible et allouer de la mémoire immédiatement sous le maximum, afin de ne pas surallouer de la mémoire, les champs nBond et mxBond , respectivement, signifient le nombre actuel de liens et le maximum. Le tableau des liaisons contiendra les indices des atomes à lier, le type de liaison et le moment de la création de la liaison (si nous sommes soudainement intéressés par des statistiques telles que la durée de vie moyenne des liaisons). Le tableau des angles contiendra les indices du triplet d'atomes formant l'angle de liaison et le type d'angle de liaison. Les tableaux bondTypes et angleTypes contiendront les caractéristiques des potentiels et angles de liaison possibles. Voici leurs structures:
struct cudaBond
{
int type; // potential type
int spec1, spec2; // type of atoms that connected by this bond type
int new_type[2]; // bond type after mutation
int new_spec1[2], new_spec2[2];
int mxEx, mnEx; // flags: maximum or minimum of bond length exists
float p0, p1, p2, p3, p4; // potential parameters
float r2min, r2max; // square of minimal and maximal bond length
float (*force_eng)(float r2, float r, float &eng, cudaBond *bond); // return energy
int count; // quantity of such bonds
float rSumm; // summ of lentghs (for mean length calculation)
int rCount; // number of measured lengths (for mean length calculation)
int ltSumm, ltCount; // for calculation of lifetime
};
struct cudaAngle
{
int type; // potential type
float p0, p1, p2; // potential parameters
void (*force_eng)(int4* angle, cudaAngle* type, cudaMD* md, float& eng);
};
Le champ type définit la forme fonctionnelle du potentiel, new_type , new_spec1 et new_spec sont les indices du type de liaison et les types d'atomes à lier après que la liaison change (se brise ou se transforme en un autre type de liaison). Ces champs sont représentés sous forme de tableaux avec deux éléments. Le premier correspond à la situation où la longueur devient plus courte que r2min 1/2 , le second - quand elle dépasse r2max 1/2... La partie la plus difficile de l'algorithme est l'application des propriétés de toutes les liaisons, en tenant compte de la possibilité de leur rupture et de leur transformation, ainsi que du fait que d'autres flux pourraient rompre les liaisons voisines, ce qui a conduit à un changement du type d'atomes liés. Permettez-moi de vous expliquer en utilisant l'exemple de la même eau. Au départ, la molécule est électriquement neutre, les liaisons chimiques sont formées par des électrons communs à l'hydrogène et à l'oxygène. En gros, on peut dire que les charges sur les atomes d'hydrogène et d'oxygène sont nulles (en fait, la densité électronique est décalée vers l'oxygène, donc il y a un petit plus sur l'hydrogène, δ +, et sur l'oxygène - 2δ-). Si nous rompons la liaison, l'oxygène prendra finalement un électron pour lui-même et l'hydrogène le cédera. Vous obtenez des particules H + et O - . Au total, 5 types de particules sont obtenus, désignons-les classiquement: H, H + , O, O- , O 2- . Ce dernier se forme si l'on détache les deux hydrogènes de la molécule d'eau. En conséquence, les réactions:
H 2 O -> H + + OH -
et
OH - -> H + + O 2- .
Les experts en chimie me corrigeront que dans des conditions standards pour l'eau, la première étape de décomposition n'est pratiquement pas mise en œuvre (en équilibre, une seule molécule sur 10 7dissocié en ions, et même alors pas tout à fait comme écrit). Mais pour la description des algorithmes, de tels schémas seront illustratifs. Supposons qu'un flux traite une liaison dans une molécule d'eau et qu'un autre flux traite la deuxième liaison de la même molécule. Et il se trouve que les deux liens doivent être rompus. Ensuite, un courant devrait transformer les atomes en H + et O - , et le second en H + et O 2- . Mais si les flux le font simultanément, au moment du début de la procédure, l'oxygène est à l'état O et les deux flux le convertissent en O - , ce qui est incorrect. Nous devons empêcher de telles situations d'une manière ou d'une autre. Schéma fonctionnel d'une fonction qui gère une liaison chimique:
On vérifie si les types d'atomes actuels correspondent au type de connexion, sinon, on prend alors dans le tableau des types par défaut (il faut le compiler au préalable), puis on détermine le carré de la distance entre les atomes (r 2 ) et, si la connexion implique une longueur maximale ou minimale, on vérifie si elle n'est pas sortie si nous sommes au-delà de ces limites. Si c'est le cas, nous devons changer le type de connexion ou le supprimer et dans les deux cas changer les types d'atomes. Pour cela, la fonction atomicCAS sera utilisée- nous comparons le type d'atome actuel avec celui qui devrait être et dans ce cas nous le remplaçons par un nouveau. Si le type de l'atome a déjà été modifié par un autre thread, nous retournons au début pour remplacer le type de lien. Le pire des cas est si nous avons réussi à changer le type du premier atome, mais pas le second. Il est déjà trop tard pour revenir en arrière, car après avoir changé le premier atome, d'autres threads pourraient déjà en faire quelque chose. Quelle est la sortie? Je suggère que nous prétendions que nous rompons / modifions une connexion d'un type différent, et non celle que nous avons prise au début. Nous trouvons quel type de connexion aurait dû être entre le premier atome initial et le second modifié et le traitons selon les mêmes règles que celles initialement prévues. Si dans ce cas le type d'atome a changé à nouveau, nous utiliserons à nouveau le même schéma. C'est implicitement implicite ici,qu'un nouveau type de liaison a les mêmes propriétés - il se rompt à la même longueur, etc., et les particules formées pendant la rupture sont au besoin. Puisque cette information est touchée par l'utilisateur, nous transférons en quelque sorte la responsabilité de notre programme sur lui, il doit tout régler correctement. Le code:
__global__ void apply_bonds(int iStep, int bndPerBlock, int bndPerThread, cudaMD* md)
{
int def;
int id1, id2; // atom indexes
int old, old_spec2, spec1, spec2, new_spec1, new_spec2; // atom types
int new_bond_type;
int save_lt, need_r, loop; // flags to save lifetime, to need to calculate r^2 and to be in ‘while’ loop
int mnmx; // flag minimum or maximum
int action; // flag: 0 - do nothing, 1 - delete bond, 2 - transform bond
cudaBond *old_bnd, *cur_bnd; // old bond type, current bond type
float dx, dy, dz, r2, r;
float f, eng = 0.0f;
__shared__ float shEng;
#ifdef DEBUG_MODE
int cnt; // count of change spec2 loops
#endif
if (threadIdx.x == 0)
{
shEng = 0.0f;
}
__syncthreads();
int id0 = blockIdx.x * bndPerBlock + threadIdx.x * bndPerThread;
int N = min(id0 + bndPerThread, md->nBond);
int iBnd;
for (iBnd = id0; iBnd < N; iBnd++)
if (md->bonds[iBnd].z) // the bond is not broken
{
// atom indexes
id1 = md->bonds[iBnd].x;
id2 = md->bonds[iBnd].y;
// atom types
spec1 = md->types[id1];
spec2 = md->types[id2];
old_bnd = &(md->bondTypes[md->bonds[iBnd].z]);
cur_bnd = old_bnd;
save_lt = 0;
need_r = 1;
loop = 1;
#ifdef DEBUG_MODE
cnt = 0;
#endif
if ((cur_bnd->spec1 == spec1)&&(cur_bnd->spec2 == spec2))
{
//ok
}
else
if ((cur_bnd->spec1 == spec2) && (cur_bnd->spec2 == spec1))
{
invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));
//... then ok
}
else // atom types do not correspond to bond types
{
save_lt = 1;
}
// end initial stage
while (loop)
{
if (save_lt)
{
def = md->def_bonds[spec1][spec2];
if (def == 0) // these atom types do not form a bond
{
#ifdef DEBUG_MODE
printf("probably, something goes wrong\n");
#endif
action = 1; // delete
break;
}
else
{
//! change bond type and go on
if (def < 0)
{
invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));
def = -def;
}
md->bonds[iBnd].z = def;
cur_bnd = &(md->bondTypes[def]);
}
} // end if (save_lt)
// calculate distance (only once)
if (need_r)
{
dx = md->xyz[id1].x - md->xyz[id2].x;
dy = md->xyz[id1].y - md->xyz[id2].y;
dz = md->xyz[id1].z - md->xyz[id2].z;
delta_periodic(dx, dy, dz, md);
r2 = dx * dx + dy * dy + dz * dz;
need_r = 0;
}
action = 0; // 0 - just cultivate bond 1 - delete bond 2 - transform bond
if ((cur_bnd->mxEx) && (r2 > cur_bnd->r2max))
{
mnmx = 1;
if (cur_bnd->new_type[mnmx] == 0) // delete bond
action = 1;
else
action = 2; // modify bond
}
else if ((cur_bnd->mnEx) && (r2 < cur_bnd->r2min))
{
mnmx = 0;
action = 2; // at minimum only bond modification possible
}
// end select action
// try to change atom types (if needed)
if (action)
{
save_lt = 1;
new_spec1 = cur_bnd->new_spec1[mnmx];
new_spec2 = cur_bnd->new_spec2[mnmx];
//the first atom
old = atomicCAS(&(md->types[id1]), spec1, new_spec1);
if (old != spec1)
{
spec1 = old;
spec2 = md->types[id2]; // refresh type of the 2nd atom
// return to begin of the ‘while’ loop
}
else // types[id1] have been changed
{
#ifdef USE_NEWANG // save changes in atom type
atomicCAS(&(md->oldTypes[id1]), -1, spec1);
#endif
old_spec2 = spec2;
while ((old = atomicCAS(&(md->types[id2]), old_spec2, new_spec2)) != old_spec2)
{
//! the worst variant: this thread changes atom 1, other thread changes atom 2
// imagine that we had A-old bond with the same behavior
def = md->def_bonds[spec1][old];
#ifdef DEBUG_MODE
if (def == 0)
{
printf("UBEH[001]: in apply_bonds, change atom types. There are no bond types between Species[%d] and [%d]\n", spec1, old);
break;
}
#endif
if (def < 0) // spec1 -> new_spec2 spec2 -> newSpec1
{
cur_bnd = &(md->bondTypes[-def]);
new_spec2 = cur_bnd->new_spec1[mnmx];
}
else // direct order
{
cur_bnd = &(md->bondTypes[def]);
new_spec2 = cur_bnd->new_spec2[mnmx];
}
old_spec2 = old;
#ifdef DEBUG_MODE
cnt++;
if (cnt > 10)
{
printf("UBEH[002]: too many atempst to change spec2 = %d\n", spec2);
break;
}
#endif
}
#ifdef USE_NEWANG // save changes in atom type
atomicCAS(&(md->oldTypes[id2]), -1, spec2);
#endif
loop = 0;
}
//end change types
} // end if (action)
else
loop = 0; // action == 0, out of cycle
} // end while(loop)
if (action == 2)
{
new_bond_type = cur_bnd->new_type[mnmx];
if (new_bond_type < 0)
{
md->bonds[iBnd].x = id2;
md->bonds[iBnd].y = id1;
new_bond_type = -new_bond_type;
}
md->bonds[iBnd].z = new_bond_type;
cur_bnd = &(md->bondTypes[new_bond_type]);
}
// perform calculations and save mean bond length
if (action != 1) // not delete
{
r = sqrt(r2);
f = cur_bnd->force_eng(r2, r, eng, cur_bnd);
atomicAdd(&(md->frs[id1].x), f * dx);
atomicAdd(&(md->frs[id2].x), -f * dx);
atomicAdd(&(md->frs[id1].y), f * dy);
atomicAdd(&(md->frs[id2].y), -f * dy);
atomicAdd(&(md->frs[id1].z), f * dz);
atomicAdd(&(md->frs[id2].z), -f * dz);
atomicAdd(&(cur_bnd->rSumm), r);
atomicAdd(&(cur_bnd->rCount), 1);
}
else //delete bond
{
// decrease the number of bonds for atoms
atomicSub(&(md->nbonds[id1]), 1);
atomicSub(&(md->nbonds[id2]), 1);
md->bonds[iBnd].z = 0;
// change parents
exclude_parents(id1, id2, md);
}
if (save_lt)
{
keep_bndlifetime(iStep, &(md->bonds[iBnd]), old_bnd);
if (action != 1) // not delete
atomicAdd(&(cur_bnd->count), 1);
atomicSub(&(old_bnd->count), 1);
}
} // end main loop
// split energy to shared and then to global memory
atomicAdd(&shEng, eng);
__syncthreads();
if (threadIdx.x == 0)
atomicAdd(&(md->engBond), shEng);
}
Dans le code, j'ai utilisé des directives de préprocesseur pour activer les vérifications des situations pouvant survenir en raison de la surveillance de l'utilisateur. Vous pouvez les désactiver pour accélérer les performances. La fonction implémente le schéma ci-dessus, mais encapsulée dans une autre boucle qui traverse la plage de liens dont ce thread est responsable. Ci-après, l'identifiant du type de liaison peut être négatif, cela signifie que l'ordre des atomes dans la liaison doit être inversé (par exemple, la liaison OH et HO est la même liaison, mais dans l'algorithme l'ordre est important, pour l'indiquer, j'utilise des indices avec le contraire sign), la fonction invert_bond la rend trop simple à décrire. Fonction delta_periodicapplique des conditions aux limites périodiques pour coordonner les différences. Si nous devons changer non seulement les liaisons, mais aussi les angles de liaison (directive USE_NEWANG ), nous devons marquer les atomes dont nous avons changé le type (nous en reparlerons plus tard). Pour exclure la re-liaison des mêmes atomes avec une liaison, le tableau des parents stocke l'index de l'un des atomes associés aux données (ce filet de sécurité ne fonctionne pas dans tous les cas, mais pour le mien c'est suffisant). Si nous rompons une sorte de connexion, nous devons supprimer les indices atomiques correspondants du tableau des parents, cela se fait par la fonction exclude_parents :
__device__ void exclude_parents(int id1, int id2, cudaMD* md)
// exclude id1 and id2 from parents of each other (if they are)
// and seek other parents if able
{
// flags to clear parent
int clear_1 = 0;
int clear_2 = 0;
int i, flag;
if (md->parents[id1] == id2)
clear_1 = 1;
if (md->parents[id2] == id1)
clear_2 = 1;
i = 0;
while ((i < md->nBond) && (clear_1 || clear_2))
{
if (md->bonds[i].z != 0)
{
flag = 0;
if (clear_1)
{
if (md->bonds[i].x == id1)
{
md->parents[id1] = md->bonds[i].y;
flag = 1;
}
else if (md->bonds[i].y == id1)
{
md->parents[id1] = md->bonds[i].y;
flag = 1;
}
if (flag)
{
clear_1 = 0;
i++;
continue;
}
}
if (clear_2)
{
if (md->bonds[i].x == id2)
{
md->parents[id2] = md->bonds[i].y;
flag = 1;
}
else if (md->bonds[i].y == id2)
{
md->parents[id2] = md->bonds[i].y;
flag = 1;
}
if (flag)
{
clear_2 = 0;
i++;
continue;
}
}
}
i++;
}
// be on the safe side
if (clear_1)
md->parents[id1] = -1;
if (clear_2)
md->parents[id2] = -1;
}
La fonction, malheureusement, parcourt tout le tableau de liens. Nous avons appris comment traiter et supprimer des liens, nous devons maintenant apprendre à les créer. La fonction suivante marque les atomes qui conviennent pour former une liaison chimique:
__device__ void try_to_bind(float r2, int id1, int id2, int spec1, int spec2, cudaMD *md)
{
int r2Int; // (int)r2 * const
// save parents to exclude re-linking
if (md->parents[id1] == id2)
return;
if (md->parents[id2] == id1)
return;
if (md->bindBonds[spec1][spec2] != 0)
{
if (r2 < md->bindR2[spec1][spec2])
{
r2Int = (int)(r2 * 100);
if (atomicMin(&(md->r2Min[id1]), r2Int) > r2Int) // replace was sucessfull
{
md->neighToBind[id1] = id2 + 1; // as 0 is reserved for no neighbour
md->canBind[id1] = 1;
}
// similar for the second atom
if (atomicMin(&(md->r2Min[id2]), r2Int) > r2Int) // replace was sucessfull
{
md->neighToBind[id2] = id1 + 1; // as 0 is reserved for no bind
md->canBind[id2] = 1;
}
}
}
}
La matrice bindBonds stocke des informations indiquant si ces types d'atomes peuvent former une liaison, et si oui, laquelle. La matrice bindR2 stocke la distance maximale entre les atomes requise pour la liaison. Si toutes les conditions sont favorables, alors nous vérifions si les atomes des autres voisins sont adaptés pour la liaison, mais plus proches.
Les informations sur la distance la plus proche du voisin sont stockées dans le tableau r2Min (pour plus de commodité, le tableau est de type int et les valeurs y sont converties avec multiplication par une constante, 100). Si le voisin détecté est le plus proche, nous nous souvenons de son index dans le tableau neighToBind et définissons l'indicateur canBind... Il y a un risque réel que pendant que nous passions à la mise à jour de l'index, un autre thread ait écrasé la valeur minimale, mais ce n'est pas si critique. Il est conseillé d'appeler cette fonction dans les fonctions qui traversent des paires d'atomes, par exemple, cell_list ou all_pair , décrites dans la première partie . La reliure elle-même:
__global__ void create_bonds(int iStep, int atPerBlock, int atPerThread, cudaMD* md)
// connect atoms which are selected to form bonds
{
int id1, id2, nei; // neighbour index
int btype, bind; // bond type index and bond index
cudaBond* bnd;
int spec1, spec2; // species indexes
int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;
int N = min(id0 + atPerThread, md->nAt);
int iat;
for (iat = id0; iat < N; iat++)
{
nei = md->neighToBind[iat];
if (nei) // neighbour exists
{
nei--; // (nei = spec_index + 1)
if (iat < nei)
{
id1 = iat;
id2 = nei;
}
else
{
id1 = nei;
id2 = iat;
}
// try to lock the first atom
if (atomicCAS(&(md->canBind[id1]), 1, 0) == 0) // the atom is already used
continue;
// try to lock the second atom
if (atomicCAS(&(md->canBind[id2]), 1, 0) == 0) // the atom is already used
{
// unlock the first one back
atomicExch(&(md->canBind[id1]), 1);
continue;
}
// create bond iat-nei
bind = atomicAdd(&(md->nBond), 1);
#ifdef DEBUG_MODE
if (bind >= md->mxBond)
{
printf("UBEH[003]: Exceed maximal number of bonds, %d\n", md->mxBond);
}
#endif
spec1 = md->types[id1];
spec2 = md->types[id2];
#ifdef USE_NEWANG // save that we have changed atom type
atomicCAS(&(md->oldTypes[id1]), -1, spec1);
atomicCAS(&(md->oldTypes[id2]), -1, spec2);
#endif
btype = md->bindBonds[spec1][spec2];
if (btype < 0)
{
// invert atoms order
md->bonds[bind].x = id2;
md->bonds[bind].y = id1;
md->bonds[bind].z = -btype;
bnd = &(md->bondTypes[-btype]);
// change atom types according the formed bond
md->types[id1] = bnd->spec2;
md->types[id2] = bnd->spec1;
}
else
{
md->bonds[bind].x = id1;
md->bonds[bind].y = id2;
md->bonds[bind].z = btype;
bnd = &(md->bondTypes[btype]);
// change atom types according the formed bond
md->types[id1] = bnd->spec1;
md->types[id2] = bnd->spec2;
}
atomicAdd((&bnd->count), 1);
md->bonds[bind].w = iStep; // keep time of the bond creation for lifetime calculation
atomicAdd(&(md->nbonds[id1]), 1);
atomicAdd(&(md->nbonds[id2]), 1);
// replace parents if none:
atomicCAS(&(md->parents[id1]), -1, id2);
atomicCAS(&(md->parents[id2]), -1, id1);
}
}
// end loop by atoms
}
La fonction bloque un atome, puis le second et, si elle réussit à bloquer les deux, crée une connexion entre eux. Au tout début de la fonction, les indices des atomes sont triés afin d'exclure la situation où un thread bloque le premier atome d'une paire, et l'autre bloque le deuxième atome de la même paire, les deux threads réussissent le premier contrôle et échouent sur le second, en conséquence aucun des deux la connexion ne les crée pas. Et enfin, nous devons supprimer les liens que nous avons marqués pour suppression dans la fonction apply_bonds :
__global__ void clear_bonds(cudaMD* md)
// clear bonds with .z == 0
{
int i = 0;
int j = md->nBond - 1;
while (i < j)
{
while ((md->bonds[j].z == 0) && (j > i))
j--;
while ((md->bonds[i].z != 0) && (i < j))
i++;
if (i < j)
{
md->bonds[i] = md->bonds[j];
md->bonds[j].z = 0;
i++;
j--;
}
}
if ((i == j) && (md->bonds[i].z == 0))
md->nBond = j;
else
md->nBond = j + 1;
}
Nous déplaçons simplement les liens «annulés» à la fin du tableau et diminuons le nombre réel de liens. Malheureusement, le code est en série, mais je ne suis pas sûr que sa mise en parallèle aura un effet tangible. Les fonctions qui calculent l'énergie de liaison réelle et les forces sur les atomes, qui sont indiquées par les champs force_eng de la structure cudaBond , sont toujours omises , mais elles sont complètement analogues aux fonctions de potentiels de paires décrites dans la première section.
Angles de Valence
Avec les angles de valence, je vais introduire quelques hypothèses pour rendre les algorithmes et les fonctions plus faciles, et par conséquent, ils seront encore plus simples que pour les liaisons de valence. Premièrement, les paramètres des angles de liaison doivent dépendre des trois atomes, mais ici nous supposerons que le type d'angle de liaison détermine exclusivement l'atome à son sommet. Je propose l'algorithme le plus simple pour former / supprimer des coins: chaque fois que nous changeons le type d'un atome, nous nous souvenons de ce fait dans le tableau correspondant oldTypes [] . La taille du tableau est égale au nombre d'atomes, initialement il est rempli de -1. Si une fonction change le type d'un atome, elle remplace -1 par l'index du type d'origine. Pour tous les atomes qui ont changé de type, supprimez tous les angles de liaison et parcourez toutes les liaisons de cet atome pour ajouter les angles correspondants:
__global__ void refresh_angles(int iStep, int atPerBlock, int atPerThread, cudaMD *md)
// delete old angles and create new ones for atoms which change their type
{
int i, j, n, t, ang;
int nei[8]; // bonded neighbors of given atom
int cnt;
int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;
int N = min(id0 + atPerThread, md->nAt);
int iat;
for (iat = id0; iat < N; iat++)
if (md->oldTypes[iat] != -1)
{
i = 0;
n = md->nangles[iat];
while (n && (i < md->nAngle))
{
if (md->angles[i].w)
if (md->angles[i].x == iat)
{
n--;
md->angles[i].w = 0;
}
i++;
}
// create new angles
t = md->specAngles[md->types[iat]]; // get type of angle, which formed by current atom type
if (t && (md->nbonds[iat] > 1)) // atom type supports angle creating and number of bonds is enough
{
// search of neighbors by bonds
i = 0; cnt = 0;
n = md->nbonds[iat];
while (n && (i < md->nBond))
{
if (md->bonds[i].z) // if bond isn't deleted
{
if (md->bonds[i].x == iat)
{
nei[cnt] = md->bonds[i].y;
cnt++;
n--;
}
else if (md->bonds[i].y == iat)
{
nei[cnt] = md->bonds[i].x;
cnt++;
n--;
}
}
i++;
}
// add new angles based on found neighbors:
for (i = 0; i < cnt-1; i++)
for (j = i + 1; j < cnt; j++)
{
ang = atomicAdd(&(md->nAngle), 1);
md->angles[ang].x = iat;
md->angles[ang].y = nei[i];
md->angles[ang].z = nei[j];
md->angles[ang].w = t;
}
n = (cnt * (cnt - 1)) / 2;
}
md->nangles[iat] = n;
// reset flag
md->oldTypes[iat] = -1;
}
}
Le tableau specAngles contient les identificateurs d'angle de liaison correspondant au type d'atome donné. La fonction suivante invoque le calcul de l'énergie et des forces pour tous les angles:
__global__ void apply_angles(int iStep, int angPerBlock, int angPerThread, cudaMD* md)
// apply valence angle potentials
{
cudaAngle* ang;
// energies of angle potential
float eng;
__shared__ float shEng;
if (threadIdx.x == 0)
shEng = 0.0f;
__syncthreads();
int id0 = blockIdx.x * angPerBlock + threadIdx.x * angPerThread;
int N = min(id0 + angPerThread, md->nAngle);
int i;
for (i = id0; i < N; i++)
if (md->angles[i].w)
{
ang = &(md->angleTypes[md->angles[i].w]);
ang->force_eng(&(md->angles[i]), ang, md, eng);
}
// split energy to shared and then to global memory
atomicAdd(&shEng, eng);
__syncthreads();
if (threadIdx.x == 0)
atomicAdd(&(md->engAngl), shEng);
}Eh bien, par exemple, le potentiel de tels angles, donnant une fonction cosinus harmonique, qui peut indiquer le champ force_eng structure cudaAngle :
__device__ void angle_hcos(int4* angle, cudaAngle* type, cudaMD* md, float& eng)
// harmonic cosine valent angle potential:
// U = k / 2 * (cos(th)-cos(th0))^
{
float k = type->p0;
float cos0 = type->p1;
// indexes of central atom and ligands:
int c = angle->x;
int l1 = angle->y;
int l2 = angle->z;
// vector ij
float xij = md->xyz[l1].x - md->xyz[c].x;
float yij = md->xyz[l1].y - md->xyz[c].y;
float zij = md->xyz[l1].z - md->xyz[c].z;
delta_periodic(xij, yij, zij, md);
float r2ij = xij * xij + yij * yij + zij * zij;
float rij = sqrt(r2ij);
// vector ik
float xik = md->xyz[l2].x - md->xyz[c].x;
float yik = md->xyz[l2].y - md->xyz[c].y;
float zik = md->xyz[l2].z - md->xyz[c].z;
delta_periodic(xik, yik, zik, md);
float r2ik = xik * xik + yik * yik + zik * zik;
float rik = sqrt(r2ik);
float cos_th = (xij * xik + yij * yik + zij * zik) / rij / rik;
float dCos = cos_th - cos0; // delta cosinus
float c1 = -k * dCos;
float c2 = 1.0 / rij / rik;
atomicAdd(&(md->frs[c].x), -c1 * (xik * c2 + xij * c2 - cos_th * (xij / r2ij + xik / r2ik)));
atomicAdd(&(md->frs[c].y), -c1 * (yik * c2 + yij * c2 - cos_th * (yij / r2ij + yik / r2ik)));
atomicAdd(&(md->frs[c].z), -c1 * (zik * c2 + zij * c2 - cos_th * (zij / r2ij + zik / r2ik)));
atomicAdd(&(md->frs[l1].x), c1 * (xik * c2 - cos_th * xij / r2ij));
atomicAdd(&(md->frs[l1].y), c1 * (yik * c2 - cos_th * yij / r2ij));
atomicAdd(&(md->frs[l1].z), c1 * (zik * c2 - cos_th * zij / r2ij));
atomicAdd(&(md->frs[l2].x), c1 * (xij * c2 - cos_th * xik / r2ik));
atomicAdd(&(md->frs[l2].y), c1 * (yij * c2 - cos_th * yik / r2ik));
atomicAdd(&(md->frs[l2].z), c1 * (zij * c2 - cos_th * zik / r2ik));
eng += 0.5 * k * dCos * dCos;
}
Je ne donnerai pas de fonction pour supprimer les coins "annulés", cela ne diffère pas fondamentalement de clear_bonds .
Exemples de
Sans prétendre être précis, j'ai essayé de représenter l'assemblage de molécules d'eau à partir d'ions individuels. Les potentiels appariés ont été fixés arbitrairement sous la forme du potentiel de Buckingham, puis ont ajouté la capacité de créer des liaisons sous la forme d'un potentiel harmonique, avec une distance d'équilibre égale à la longueur de la liaison HO dans l'eau, 0,96 Å. De plus, lorsque le deuxième proton était lié à l'oxygène, un angle de liaison avec le sommet dans l'oxygène a été ajouté. Après 100'000 étapes, les premières molécules sont apparues à partir d'ions dispersés aléatoirement. La figure montre les configurations initiale (gauche) et finale (droite):
Vous pouvez mettre en place une expérience comme celle-ci: laissez d'abord tous les atomes sont identiques, mais lorsqu'ils sont côte à côte, ils forment une liaison. Laissez les atomes liés former une autre liaison soit avec un atome libre, soit avec une autre molécule liée similaire. En conséquence, nous obtiendrons une sorte d'auto-organisation, où les atomes s'alignent en chaînes:
Commentaires finaux
- Ici, nous n'avons utilisé qu'un seul critère pour la liaison - la distance, bien qu'en principe il puisse y en avoir d'autres, par exemple, l'énergie du système. En réalité, lorsqu'une liaison chimique est formée, en règle générale, de l'énergie est libérée sous forme de chaleur. Ceci n'est pas pris en compte ici, mais vous pouvez essayer de l'implémenter, par exemple changer la vitesse des particules.
- Les interactions entre les particules par le biais du potentiel de liaison chimique ne nient pas le fait que les particules peuvent encore interagir par le biais de potentiels de paires intermoléculaires et de l'interaction de Coulomb. Il serait bien sûr possible de ne pas calculer les interactions intermoléculaires pour les atomes liés, mais cela, dans le cas général, nécessite de longues vérifications. Par conséquent, il est plus facile de régler le potentiel de liaison chimique de telle sorte que sa somme avec d'autres potentiels donne la fonction souhaitée.
- La mise en œuvre parallèle de la liaison de particules donne non seulement un gain de vitesse, mais semble même plus réaliste, car les particules se font concurrence.
Eh bien, il y a plusieurs projets sur Habré qui sont très proches de celui-ci: