For example, water:
it is obvious that hydrogen and oxygen inside one molecule interact in a completely different way than the same oxygen with the hydrogen of a neighboring molecule. Thus, the intramolecular and intermolecular interactions are distinguished. The intermolecular interaction can be specified by short-range and Coulomb pair potentials, which were discussed in previous articles. Here we will focus on the intramolecular.
The most common type of intramolecular interaction is chemical (valence) bonds. Chemical bonds are set by the functional dependence of the potential energy on the distance between the bound atoms, i.e., in fact, by the same pair potential. But, unlike ordinary pair potentials, this interaction is specified not for certain types of particles, but for a specific pair of particles (by their indices). The most common functional forms for chemical bond potentials are harmonic potentials:
where r is the distance between particles, k is the bond stiffness constant, and r 0 is the equilibrium bond length; and Morse potential :
where D is the depth of the potential well, the parameter α characterizes the width of the potential well.
The next type of intramolecular interactions is bond angles. Let's look at the title picture again. Why is the molecule depicted by a corner, because the forces of electrostatics were supposed to provide the maximum distance between hydrogen ions, which corresponds to an angle HOH equal to 180 °? The fact is that not everything is drawn in the figure. From the school chemistry course, you can recall that oxygen has 2 more lone electron pairs, the interaction with which distorts the angle:
In classical molecular dynamics, objects such as electrons or electron clouds are usually not introduced, therefore, the valence angle potential is used to simulate "correct" angles; functional dependence of the potential energy on the coordinates of 3 particles. One of the most convenient such potentials is the harmonic cosine:
where θ is the angle formed by the triplet of particles, k is the stiffness constant, and θ 0 is the equilibrium angle.
There are intramolecular potentials of a higher order, for example, torsion angles , but they are even more artificial than bond angles.
Adding interactions between particles with predetermined indices is trivial. For links, we store an array containing the indices of the linked particles and the type of interaction. We give each thread its own range of links for processing and you're done. Likewise with bond angles. Therefore, we will immediately complicate the task for ourselves: we will add the ability to create / remove chemical bonds and runtime bond angles. This immediately takes us out of the plane of classical molecular dynamics and opens up a new horizon of possibilities. Otherwise, you could simply download something from the existing packages, for example LAMMPS , DL_POLY or GROMACS , especially since they are distributed free of charge.
Now for some code. Let's add the appropriate fields to the main structure:
//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
The number of links and angles is variable, but almost always you can estimate the maximum possible and allocate memory immediately to the maximum, so as not to over-allocate memory, the fields nBond and mxBond , respectively, mean the current number of links and the maximum. The bonds array will contain the indices of the atoms to be bound, the type of bond and the time of the bond creation (if we are suddenly interested in such statistics as the average bond lifetime). The angles array will store the indices of the three atoms that make up the bond angle and the type of bond angle. The bondTypes and angleTypes arrays will contain the characteristics of the possible bond potentials and angles. Here are their 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);
};
The type field defines the functional form of the potential, new_type , new_spec1 and new_spec are the indices of the bond type and the types of atoms to be bonded after the bond changes (breaks or turns into a different type of bond). These fields are represented as arrays with two elements. The first corresponds to the situation when the length becomes shorter than r2min 1/2 , the second - when it exceeds r2max 1/2... The most difficult part of the algorithm is the application of the properties of all bonds, taking into account the possibility of their breakage and transformation, as well as the fact that other flows could break off neighboring bonds, which led to a change in the type of bound atoms. Let me explain using the example of the same water. Initially, the molecule is electrically neutral, chemical bonds are formed by electrons common to hydrogen and oxygen. Roughly speaking, we can say that the charges on the hydrogen and oxygen atoms are zero (in fact, the electron density is shifted to oxygen, therefore, there is a small plus on hydrogen, δ +, and on oxygen - 2δ-). If we break the bond, oxygen will finally take an electron for itself, and hydrogen will give it away. The resulting particles are H + and O - . In total, 5 types of particles are obtained, let's conventionally designate them: H, H + , O, O- , O 2- . The latter is formed if we detach both hydrogens from the water molecule. Accordingly, the reactions:
H 2 O -> H + + OH -
and
OH - -> H + + O 2- .
Experts in chemistry will correct me that under standard conditions for water, the first stage of decomposition is practically not implemented (in equilibrium, only one molecule out of 10 7dissociated into ions, and even then not quite as written). But for the description of algorithms, such schemes will be illustrative. Suppose one stream processes one bond in a water molecule, and another stream processes the second bond of the same molecule. And it so happened that both ties need to be broken. Then one stream should transform atoms into H + and O - , and the second into H + and O 2- . But if the flows do this simultaneously, at the time of the beginning of the procedure, oxygen is in the O state and both flows convert it into O - , which is incorrect. We need to prevent such situations somehow. Block diagram of a function that handles a chemical bond:
We check whether the current types of atoms correspond to the type of connection, if not, then we take from the table of default types (it must be compiled in advance), then we determine the square of the distance between atoms (r 2 ) and, if the connection implies a maximum or minimum length, we check whether it did not come out whether we are beyond these boundaries. If we did, then we need to change the type of connection or delete it and in both cases change the types of atoms. For this, the atomicCAS function will be used- we compare the current type of atom with the one that should be and in this case we replace it with a new one. If the atom's type has already been changed by another thread, we go back to the beginning to override the link type. The worst case scenario is if we managed to change the type of the first atom, but not the second. It's already too late to go back, because after we changed the first atom, other threads could already do something with it. What is the way out? I suggest that we pretend that we are breaking / changing a connection of a different type, and not the one we took in the beginning. We find what type of connection should have been between the initial first atom and the changed second one and process it according to the same rules as originally expected. If in this case the type of the atom has changed again, we will use the same scheme again. It is implicitly implied here,that a new type of bond has the same properties - it breaks at the same length, etc., and the particles formed during the break are as needed. Since this information is touched by the user, we kind of shift the responsibility from our program onto him, he must set everything correctly. The 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);
}
In the code, I used preprocessor directives to enable checks for situations that may arise due to user oversight. You can turn them off to speed up the performance. The function implements the above scheme, but wrapped in one more loop that runs through the range of links for which this thread is responsible. Hereinafter, the identifier of the bond type can be negative, this means that the order of atoms in the bond must be reversed (for example, the OH and HO bond is the same bond, but in the algorithm the order is important, to indicate this, I use indices with the opposite sign), the invert_bond function makes it too trivial to describe. Delta_periodic functionapplies periodic boundary conditions to coordinate differences. If we need to change not only bonds, but also bond angles ( USE_NEWANG directive ), then we need to mark the atoms for which we have changed the type (more on that later). To exclude the re-binding of the same atoms with a bond, the parents array stores the index of one of the atoms associated with the data (this safety net does not work in all cases, but for mine it is enough). If we break some kind of connection, then we need to remove the corresponding atomic indices from the parents array, this is done by the exclude_parents function :
__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;
}
The function, unfortunately, runs through the entire array of links. We learned how to process and delete links, now we need to learn how to create them. The following function marks the atoms that are suitable for forming a chemical bond:
__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;
}
}
}
}
The bindBonds matrix stores information about whether these types of atoms can form a bond, and if so, which one. The bindR2 matrix stores the maximum distance between atoms required for binding. If all the conditions are favorable, then we check whether the atoms of other neighbors are suitable for bonding, but closer.
Information about the nearest distance to the neighbor is stored in the r2Min array (for convenience, the array is of type int and values are converted to it with multiplication by a constant, 100). If the detected neighbor is the closest one, then we remember its index in the neighToBind array and set the canBind flag... There is a real danger that while we moved on to updating the index, another thread has overwritten the minimum value, but this is not so critical. It is advisable to call this function in functions that traverse pairs of atoms, for example, cell_list or all_pair , described in the first part . The binding itself:
__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
}
The function blocks one atom, then the second and, if it succeeds in blocking both, creates a connection between them. At the very beginning of the function, the atomic indices are sorted in order to exclude the situation when one thread blocks the first atom in a pair, and the other blocks the second atom in the same pair, both threads successfully pass the first check and fail on the second, as a result neither of the the connection does not create them. And finally, we need to remove the links that we marked for deletion in the apply_bonds function :
__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;
}
We simply move the "nullified" links to the end of the array and decrease the actual number of links. Unfortunately, the code is serial, but I'm not sure that parallelizing it will bring any tangible effect. The functions that calculate the actual binding energy and forces on atoms, which are indicated by the force_eng fields of the cudaBond structure, are still left out , but they are completely analogous to the functions of pair potentials described in the first section.
Valence angles
With valence angles, I will introduce a few assumptions to make the algorithms and functions easier, as a result they will be even simpler than for valence bonds. First, the parameters of the bond angles should depend on all three atoms, but here we will assume that the type of bond angle determines exclusively the atom at its vertex. I propose the simplest algorithm for forming / removing corners: whenever we change the type of an atom, we remember this fact in the corresponding array oldTypes [] . The size of the array is equal to the number of atoms, initially it is filled with -1. If a function changes the type of an atom, it replaces -1 with the index of the original type. For all atoms that have changed the type, remove all bond angles and run over all bonds of this atom to add the corresponding angles:
__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;
}
}
The specAngles array contains the bond angle identifiers corresponding to the given atom type. The following function invokes the calculation of energy and forces for all 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);
}Well, for example, the potential of such angles, giving a harmonic cosine function, which may indicate the field 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;
}
I will not give a function for removing "nullified" corners, it does not fundamentally differ from clear_bonds .
Examples of
Without pretending to be accurate, I tried to depict the assembly of water molecules from individual ions. Paired potentials were set arbitrarily in the form of Buckingham's potential, and then added the ability to create bonds in the form of a harmonic potential, with an equilibrium distance equal to the length of the HO bond in water, 0.96 Å. In addition, when the second proton was bound with oxygen, a bond angle with the apex in oxygen was added. After 100'000 steps, the first molecules appeared from randomly scattered ions. The figure shows the initial (left) and final (right) configurations:
You can set up an experiment like this: let at first all the atoms are the same, but when they are next to each other they form a bond. Let the bound atoms form another bond either with a free atom or with another similar bound molecule. As a result, we will get some kind of self-organization, where the atoms line up in chains:
Final comments
- Here we used only one criterion for binding - distance, although in principle there may be others, for example, the energy of the system. In reality, when a chemical bond is formed, as a rule, energy is released in the form of heat. This is not taken into account here, but you can try to implement it, for example, change the particle speed.
- Interactions between particles through the chemical bond potential do not negate the fact that particles can still interact through intermolecular pair potentials and the Coulomb interaction. It would be possible, of course, not to calculate intermolecular interactions for bound atoms, but this, in the general case, requires long checks. Therefore, it is easier to set the chemical bond potential in such a way that its sum with other potentials gives the desired function.
- Parallel implementation of particle binding not only gives a speed gain, but even looks more realistic, since the particles compete with each other.
Well, there are several projects on Habré that are very close to this one: