Zum Beispiel Wasser:
Es ist offensichtlich, dass Wasserstoff und Sauerstoff in einem Molekül auf völlig andere Weise als der gleiche Sauerstoff mit dem Wasserstoff eines benachbarten Moleküls interagieren. Somit werden die intramolekularen und intermolekularen Wechselwirkungen unterschieden. Intermolekulare Wechselwirkungen können durch Nahbereichs- und Coulomb-Paarpotentiale spezifiziert werden, die in früheren Artikeln diskutiert wurden. Hier konzentrieren wir uns auf das Intramolekulare.
Die häufigste Art der intramolekularen Wechselwirkung sind chemische (Valenz-) Bindungen. Chemische Bindungen werden durch die funktionelle Abhängigkeit der potentiellen Energie vom Abstand zwischen den gebundenen Atomen festgelegt, d. H. Tatsächlich durch das gleiche Paarpotential. Im Gegensatz zu gewöhnlichen Paarpotentialen wird diese Wechselwirkung jedoch nicht für bestimmte Partikeltypen, sondern für ein bestimmtes Partikelpaar (anhand ihrer Indizes) spezifiziert. Die häufigsten funktionellen Formen für chemische Bindungspotentiale sind harmonische Potentiale:
wobei r der Abstand zwischen Partikeln ist, k die Bindungssteifigkeitskonstante ist und r 0 die Gleichgewichtsbindungslänge ist; und Morsepotential :
wobei D die Tiefe der Potentialwanne ist, charakterisiert der Parameter α die Breite der Potentialwanne.
Die nächste Art von intramolekularen Wechselwirkungen sind Bindungswinkel. Schauen wir uns noch einmal das Titelbild an. Warum ist das Molekül durch eine Ecke dargestellt, weil elektrostatische Kräfte den maximalen Abstand zwischen Wasserstoffionen liefern sollten, der einem Winkel HOH von 180 ° entspricht? Tatsache ist, dass nicht alles in der Figur gezeichnet ist. Aus dem Chemiekurs der Schule können Sie sich daran erinnern, dass Sauerstoff zwei weitere Einzelelektronenpaare hat, deren Wechselwirkung den Winkel verzerrt:
In der klassischen Molekulardynamik werden Objekte wie Elektronen oder Elektronenwolken normalerweise nicht eingeführt, daher wird das Valenzwinkelpotential verwendet, um "korrekte" Winkel zu simulieren; funktionelle Abhängigkeit der potentiellen Energie von den Koordinaten von 3 Teilchen. Eines der bequemsten derartigen Potentiale ist der harmonische Kosinus:
wobei θ der Winkel ist, der durch das Triplett von Partikeln gebildet wird, k die Steifheitskonstante ist und θ 0 der Gleichgewichtswinkel ist.
Es gibt intramolekulare Potentiale höherer Ordnung, zum Beispiel Torsionswinkel , aber sie sind noch künstlicher als Bindungswinkel.
Das Hinzufügen von Wechselwirkungen zwischen Partikeln mit vorgegebenen Indizes ist trivial. Für Links speichern wir ein Array, das die Indizes der verknüpften Partikel und die Art der Interaktion enthält. Wir geben jedem Thread seine eigenen Links zur Verarbeitung und Sie sind fertig. Ebenso bei Bindungswinkeln. Daher werden wir die Aufgabe sofort für uns selbst komplizieren: Wir werden die Möglichkeit hinzufügen, chemische Bindungen und Laufzeitbindungswinkel zu erzeugen / entfernen. Dies führt uns sofort aus der Ebene der klassischen Molekulardynamik heraus und eröffnet einen neuen Horizont von Möglichkeiten. Andernfalls können Sie einfach etwas von den vorhandenen Paketen herunterladen, z. B. LAMMPS , DL_POLY oder GROMACS , zumal diese kostenlos verteilt werden.
Nun zu etwas Code. Fügen wir der Hauptstruktur die entsprechenden Felder hinzu:
//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
Die Anzahl der Verknüpfungen und Winkel ist variabel, aber fast immer können Sie das maximal mögliche schätzen und den Speicher sofort dem Maximum zuordnen, um den Speicher nicht zu stark zuzuweisen. Die Felder nBond und mxBond bedeuten die aktuelle Anzahl der Verknüpfungen und das Maximum. Das Bindungsarray enthält die Indizes der zu bindenden Atome, die Art der Bindung und den Zeitpunkt der Bindungsbildung (wenn wir plötzlich an Statistiken wie der durchschnittlichen Bindungslebensdauer interessiert sind). Das Winkelarray enthält die Indizes des Tripletts von Atomen, die den Bindungswinkel und die Art des Bindungswinkels bilden. Die Arrays bondTypes und angleTypes enthalten die Eigenschaften der möglichen Bindungspotentiale und -winkel. Hier sind ihre Strukturen:
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);
};
Das Typ - Feld definiert die Funktionsform des Potentials, new_type , new_spec1 und new_spec sind die Indizes der Bindungstyp und die Art der Atome werden gebunden , nachdem die Bindung der Änderungen (Brüche oder verwandelt sich in eine andere Art der Bindung). Diese Felder werden als Arrays mit zwei Elementen dargestellt. Die erste entspricht der Situation, wenn die Länge kürzer als r2min 1/2 wird , die zweite - wenn sie r2max 1/2 überschreitet... Der schwierigste Teil des Algorithmus ist die Anwendung der Eigenschaften aller Bindungen unter Berücksichtigung der Möglichkeit ihres Bruchs und ihrer Umwandlung sowie der Tatsache, dass andere Strömungen benachbarte Bindungen abbrechen könnten, was zu einer Änderung der Art der gebundenen Atome führte. Lassen Sie mich am Beispiel des gleichen Wassers erklären. Anfangs ist das Molekül elektrisch neutral, chemische Bindungen werden durch Elektronen gebildet, die Wasserstoff und Sauerstoff gemeinsam haben. Grob gesagt können wir sagen, dass die Ladungen an den Wasserstoff- und Sauerstoffatomen Null sind (tatsächlich ist die Elektronendichte zu Sauerstoff verschoben, daher gibt es ein kleines Plus für Wasserstoff, δ + und für Sauerstoff - 2δ-). Wenn wir die Bindung aufbrechen, nimmt Sauerstoff schließlich ein Elektron für sich auf und Wasserstoff gibt es ab. Die resultierenden Partikel sind H + und O - . Insgesamt erhalten wir 5 Arten von Partikeln, die wir herkömmlicherweise bezeichnen: H, H + , O, O.- , O 2- . Letzteres entsteht, wenn wir beide Wasserstoffatome vom Wassermolekül ablösen. Dementsprechend sind die Reaktionen:
H 2 O -> H + + OH -
und
OH - -> H + + O 2- .
Chemieexperten werden mich korrigieren, dass unter Standardbedingungen für Wasser die erste Zersetzungsstufe praktisch nicht durchgeführt wird (im Gleichgewicht nur ein Molekül von 10 7dissoziiert in Ionen und selbst dann nicht ganz so wie geschrieben). Für die Beschreibung von Algorithmen werden solche Schemata jedoch veranschaulichend sein. Angenommen, ein Strom verarbeitet eine Bindung in einem Wassermolekül und ein anderer Strom verarbeitet die zweite Bindung desselben Moleküls. Und so kam es, dass beide Bindungen gebrochen werden müssen. Dann sollte ein Strom Atome in H + und O - und der zweite in H + und O 2- umwandeln . Wenn die Ströme dies jedoch gleichzeitig tun, befindet sich Sauerstoff zu Beginn des Verfahrens im O-Zustand und beide Ströme wandeln ihn in O - um , was falsch ist. Wir müssen solche Situationen irgendwie verhindern. Blockdiagramm einer Funktion, die eine chemische Bindung handhabt:
Wir prüfen, ob die aktuellen Atomtypen der Art der Verbindung entsprechen. Wenn nicht, entnehmen wir der Tabelle der Standardtypen (sie muss im Voraus kompiliert werden), bestimmen das Quadrat des Abstands zwischen Atomen (r 2 ) und prüfen, ob die Verbindung eine maximale oder minimale Länge impliziert, ob sie nicht herausgekommen ist ob wir jenseits dieser Grenzen sind. Wenn ja, müssen wir die Art der Verbindung ändern oder löschen und in beiden Fällen die Art der Atome ändern. Hierzu wird die atomicCAS- Funktion verwendet- Wir vergleichen den aktuellen Atomtyp mit dem, der sein sollte, und ersetzen ihn in diesem Fall durch einen neuen. Wenn der Atomtyp bereits von einem anderen Thread geändert wurde, kehren wir zum Anfang zurück, um den Verknüpfungstyp zu überschreiben. Das schlimmste Szenario ist, wenn es uns gelungen ist, den Typ des ersten Atoms zu ändern, nicht jedoch des zweiten. Es ist schon zu spät, um zurück zu gehen, denn nachdem wir das erste Atom gewechselt haben, könnten andere Threads bereits etwas damit anfangen. Was ist der Ausweg? Ich schlage vor, dass wir so tun, als würden wir eine Verbindung eines anderen Typs unterbrechen / ändern, und nicht die, die wir am Anfang angegangen sind. Wir finden heraus, welche Art von Verbindung zwischen dem ersten und dem geänderten zweiten Atom hätte bestehen sollen, und verarbeiten sie nach den gleichen Regeln wie ursprünglich erwartet. Wenn sich in diesem Fall der Atomtyp erneut geändert hat, verwenden wir das gleiche Schema erneut. Es ist hier implizit impliziert,dass eine neue Art von Bindung die gleichen Eigenschaften hat - sie bricht bei der gleichen Länge usw. und die während des Bruchs gebildeten Partikel sind nach Bedarf. Da diese Informationen vom Benutzer berührt werden, verlagern wir die Verantwortung von unserem Programm auf ihn, er muss alles richtig einstellen. Der 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);
}
Im Code habe ich Präprozessoranweisungen verwendet, um Überprüfungen für Situationen zu ermöglichen, die aufgrund von Benutzerüberwachung auftreten können. Sie können sie ausschalten, um die Leistung zu beschleunigen. Die Funktion implementiert das obige Schema, ist jedoch in eine weitere Schleife eingeschlossen, die den Linkbereich durchläuft, für den dieser Thread verantwortlich ist. Im Folgenden kann der Bezeichner des Bindungstyps negativ sein. Dies bedeutet, dass die Reihenfolge der Atome in der Bindung umgekehrt werden muss (zum Beispiel ist die OH- und HO-Bindung dieselbe Bindung, aber im Algorithmus ist die Reihenfolge wichtig, um dies anzuzeigen, verwende ich Indizes mit dem Gegenteil sign), die Funktion invert_bond macht es zu trivial, um es zu beschreiben. Delta_periodische Funktionwendet periodische Randbedingungen an, um Differenzen zu koordinieren. Wenn wir nicht nur Bindungen, sondern auch Bindungswinkel ändern müssen ( USE_NEWANG- Direktive ), müssen wir die Atome markieren, für die wir den Typ geändert haben (dazu später mehr). Um die erneute Bindung derselben Atome mit einer Bindung auszuschließen, speichert das übergeordnete Array den Index eines der mit den Daten verknüpften Atome (dieses Sicherheitsnetz funktioniert nicht in allen Fällen, aber für meine ist es ausreichend). Wenn wir eine Verbindung trennen , müssen wir die entsprechenden Atomindizes aus dem übergeordneten Array entfernen. Dies geschieht durch die Funktion 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;
}
Die Funktion läuft leider über das gesamte Array von Links. Wir haben gelernt, wie man Links verarbeitet und löscht. Jetzt müssen wir lernen, wie man sie erstellt. Die folgende Funktion markiert die Atome, die zur Bildung einer chemischen Bindung geeignet sind:
__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;
}
}
}
}
Die bindBonds- Matrix speichert Informationen darüber, ob und welche Arten von Atomen eine Bindung eingehen können. Die bindR2- Matrix speichert den maximalen Abstand zwischen Atomen, der für die Bindung erforderlich ist. Wenn alle Bedingungen günstig sind, prüfen wir, ob die Atome anderer Nachbarn für die Bindung geeignet sind, aber näher.
Informationen über die nächstgelegene Entfernung zum Nachbarn werden im r2Min- Array gespeichert (der Einfachheit halber ist das Array vom Typ int und die Werte werden mit Multiplikation mit einer Konstanten 100 in dieses Array konvertiert ). Wenn der erkannte Nachbar der nächstgelegene ist, merken wir uns seinen Index im Array neighToBind und setzen das canBind- Flag... Es besteht die reale Gefahr, dass ein anderer Thread den Mindestwert überschrieb, während wir mit der Aktualisierung des Index fortfuhren. Dies ist jedoch nicht so kritisch. Es ist ratsam, diese Funktion in Funktionen aufzurufen , die Atompaare durchlaufen, z. B. cell_list oder all_pair , wie im ersten Teil beschrieben . Die Bindung selbst:
__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
}
Die Funktion blockiert ein Atom, dann das zweite und stellt, wenn es gelingt, beide zu blockieren, eine Verbindung zwischen ihnen her. Ganz am Anfang der Funktion werden die Atomindizes sortiert, um die Situation auszuschließen, dass ein Thread das erste Atom in einem Paar blockiert und der andere das zweite Atom in demselben Paar blockiert. Beide Threads bestehen die erste Prüfung erfolgreich und schlagen beim zweiten fehl, was dazu führt, dass keiner der beiden Die Verbindung erstellt sie nicht. Und schließlich müssen wir die Links entfernen, die wir zum Löschen in der Funktion apply_bonds markiert haben :
__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;
}
Wir verschieben einfach die "annullierten" Links an das Ende des Arrays und verringern die tatsächliche Anzahl der Links. Leider ist der Code seriell, aber ich bin mir nicht sicher, ob die Parallelisierung einen spürbaren Effekt hat. Die Funktionen, die die tatsächliche Bindungsenergie und die Kräfte auf Atome berechnen, die durch die force_eng- Felder der cudaBond- Struktur angezeigt werden, bleiben weiterhin aus , sind jedoch völlig analog zu den im ersten Abschnitt beschriebenen Funktionen der Paarpotentiale.
Valenzwinkel
Bei Valenzwinkeln werde ich einige Annahmen einführen, um die Algorithmen und Funktionen zu vereinfachen. Infolgedessen sind sie noch einfacher als bei Valenzbindungen. Erstens sollten die Parameter der Bindungswinkel von allen drei Atomen abhängen, aber hier nehmen wir an, dass die Art des Bindungswinkels ausschließlich das Atom an seinem Scheitelpunkt bestimmt. Ich schlage den einfachsten Algorithmus zum Bilden / Entfernen von Ecken vor: Wenn wir den Typ eines Atoms ändern, erinnern wir uns an diese Tatsache im entsprechenden Array oldTypes [] . Die Größe des Arrays entspricht der Anzahl der Atome, zunächst ist es mit -1 gefüllt. Wenn eine Funktion den Typ eines Atoms ändert, ersetzt sie -1 durch den Index des ursprünglichen Typs. Entfernen Sie für alle Atome, die den Typ geändert haben, alle Bindungswinkel und fahren Sie über alle Bindungen dieses Atoms, um die entsprechenden Winkel hinzuzufügen:
__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;
}
}
Das specAngles-Array enthält die Bindungswinkelkennungen, die dem angegebenen Atomtyp entsprechen. Die folgende Funktion ruft die Berechnung von Energie und Kräften für alle Winkel auf:
__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);
}Nun, zum Beispiel das Potential solcher Winkel, das eine harmonische Kosinusfunktion ergibt, die das Feld force_eng Struktur cudaAngle anzeigen kann :
__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;
}
Ich werde keine Funktion zum Entfernen von "nullifizierten" Ecken angeben , sie unterscheidet sich nicht grundlegend von clear_bonds .
Beispiele von
Ohne vorzugeben, genau zu sein, versuchte ich, die Anordnung von Wassermolekülen aus einzelnen Ionen darzustellen. Gepaarte Potentiale wurden willkürlich in Form des Buckinghamschen Potentials eingestellt und fügten dann die Fähigkeit hinzu, Bindungen in Form eines harmonischen Potentials mit einem Gleichgewichtsabstand zu erzeugen, der der Länge der HO-Bindung in Wasser von 0,96 Å entspricht. Wenn das zweite Proton mit Sauerstoff gebunden wurde, wurde zusätzlich ein Bindungswinkel mit der Spitze in Sauerstoff hinzugefügt. Nach 100'000 Schritten erschienen die ersten Moleküle aus zufällig gestreuten Ionen. Die Abbildung zeigt die anfängliche (links) und endgültige (rechts) Konfiguration:
Sie können ein Experiment wie dieses einrichten: Lassen Sie zunächst alle Atome gleich sein, aber wenn sie nebeneinander liegen, bilden sie eine Bindung. Lassen Sie die gebundenen Atome entweder mit einem freien Atom oder mit einem anderen ähnlich gebundenen Molekül eine weitere Bindung eingehen. Als Ergebnis erhalten wir eine Art Selbstorganisation, bei der die Atome in Ketten angeordnet sind:
Letzte Kommentare
- Hier haben wir nur ein Kriterium für die Bindungsentfernung verwendet, obwohl es im Prinzip andere geben kann, zum Beispiel die Energie des Systems. In der Realität wird bei der Bildung einer chemischen Bindung in der Regel Energie in Form von Wärme freigesetzt. Dies wird hier nicht berücksichtigt, aber Sie können versuchen, es zu implementieren, indem Sie beispielsweise die Partikelgeschwindigkeit ändern.
- Wechselwirkungen zwischen Partikeln durch das chemische Bindungspotential negieren nicht die Tatsache, dass Partikel immer noch durch intermolekulare Paarpotentiale und die Coulomb-Wechselwirkung interagieren können. Es wäre natürlich möglich, intermolekulare Wechselwirkungen für gebundene Atome nicht zu berechnen, aber dies erfordert im allgemeinen Fall lange Überprüfungen. Daher ist es einfacher, das chemische Bindungspotential so einzustellen, dass seine Summe mit anderen Potentialen die gewünschte Funktion ergibt.
- Die parallele Implementierung der Partikelbindung führt nicht nur zu einem Geschwindigkeitsgewinn, sondern sieht sogar realistischer aus, da die Partikel miteinander konkurrieren.
Nun, es gibt mehrere Projekte auf Habré, die diesem sehr nahe stehen: