Por ejemplo, agua:
es obvio que el hidrógeno y el oxígeno dentro de una molécula interactúan de una manera completamente diferente que el mismo oxígeno con el hidrógeno de una molécula vecina. Así, se distinguen las interacciones intramoleculares e intermoleculares. Las interacciones intermoleculares se pueden especificar mediante potenciales de par de Coulomb y de corto alcance, que se discutieron en artículos anteriores. Aquí nos centraremos en lo intramolecular.
El tipo más común de interacción intramolecular son los enlaces químicos (valencia). Los enlaces químicos se establecen por la dependencia funcional de la energía potencial de la distancia entre los átomos unidos, es decir, de hecho, por el mismo par de potencial. Pero, a diferencia de los potenciales de par ordinarios, esta interacción no se especifica para ciertos tipos de partículas, sino para un par específico de partículas (por sus índices). Las formas funcionales más comunes para los potenciales de enlace químico son los potenciales armónicos:
donde r es la distancia entre partículas, k es la constante de rigidez del enlace y r 0 es la longitud del enlace de equilibrio; y potencial Morse :
donde D es la profundidad del pozo potencial, el parámetro α caracteriza el ancho del pozo potencial.
El siguiente tipo de interacciones intramoleculares son los ángulos de enlace. Veamos de nuevo la imagen del título. ¿Por qué la molécula está representada con una esquina, porque se suponía que las fuerzas electrostáticas proporcionaban la distancia máxima entre los iones de hidrógeno, que corresponde a un ángulo HOH igual a 180 °? El caso es que no todo está dibujado en la figura. Del curso de química de la escuela, puedes recordar que el oxígeno tiene 2 pares de electrones solitarios más, cuya interacción distorsiona el ángulo:
En la dinámica molecular clásica, generalmente no se introducen objetos como electrones o nubes de electrones, por lo tanto, el potencial del ángulo de valencia se usa para simular ángulos "correctos"; dependencia funcional de la energía potencial de las coordenadas de 3 partículas. Uno de los potenciales más convenientes es el coseno armónico:
donde θ es el ángulo formado por el triple de partículas, k es la constante de rigidez y θ 0 es el ángulo de equilibrio.
Hay potenciales intramoleculares de orden superior, por ejemplo, ángulos de torsión , pero son incluso más artificiales que los ángulos de enlace.
Agregar interacciones entre partículas con índices predeterminados es trivial. Para los enlaces, almacenamos una matriz que contiene los índices de las partículas enlazadas y el tipo de interacción. Le damos a cada hilo su propia gama de enlaces para su procesamiento y listo. Lo mismo ocurre con los ángulos de enlace. Por lo tanto, nos complicaremos la tarea de inmediato: agregaremos la capacidad de crear / eliminar enlaces químicos y ángulos de enlace en tiempo de ejecución. Esto nos saca inmediatamente del plano de la dinámica molecular clásica y abre un nuevo horizonte de posibilidades. De lo contrario, podría simplemente descargar algo de los paquetes existentes, por ejemplo LAMMPS , DL_POLY o GROMACS , especialmente porque se distribuyen de forma gratuita.
Ahora, un poco de código. Agreguemos los campos apropiados a la estructura principal:
//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
El número de enlaces y ángulos es variable, pero casi siempre se puede estimar el máximo posible y asignar memoria inmediatamente al máximo, para no sobreasignar memoria, los campos nBond y mxBond , respectivamente, significan el número actual de enlaces y el máximo. La matriz de enlaces contendrá los índices de los átomos que se unirán, el tipo de enlace y el tiempo de creación del enlace (si de repente nos interesan estadísticas como la vida media del enlace). La matriz de ángulos contendrá los índices del triplete de átomos que forman el ángulo de enlace y el tipo de ángulo de enlace. Los arreglos bondTypes y angleTypes contendrán las características de los posibles ángulos y potenciales de enlace. Aquí están sus estructuras:
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);
};
El campo de tipo define la forma funcional del potencial, new_type , new_spec1 y new_spec son los índices del tipo de enlace y los tipos de átomos que se unirán después de que el enlace cambie (se rompa o se convierta en un tipo de enlace diferente). Estos campos se representan como matrices con dos elementos. El primero corresponde a la situación en la que la longitud se vuelve más corta que r2min 1/2 , el segundo - cuando excede r2max 1/2... La parte más difícil del algoritmo es la aplicación de las propiedades de todos los enlaces, teniendo en cuenta la posibilidad de su ruptura y transformación, así como el hecho de que otros flujos podrían romper los enlaces vecinos, lo que provocó un cambio en el tipo de átomos ligados. Déjame explicarte usando el ejemplo de la misma agua. Inicialmente, la molécula es eléctricamente neutra, los enlaces químicos están formados por electrones comunes al hidrógeno y al oxígeno. En términos generales, podemos decir que las cargas en los átomos de hidrógeno y oxígeno son cero (de hecho, la densidad de electrones se desplaza a oxígeno, por lo tanto, hay una pequeña ventaja en el hidrógeno, δ + y en el oxígeno - 2δ-). Si rompemos el enlace, el oxígeno finalmente tomará un electrón para sí mismo y el hidrógeno lo regalará. Las partículas resultantes son H + y O - . En total se obtienen 5 tipos de partículas, designémoslas convencionalmente: H, H + , O, O- , O 2- . Este último se forma si separamos ambos hidrógenos de la molécula de agua. En consecuencia, las reacciones:
H 2 O -> H + + OH -
y
OH - -> H + + O 2- .
Los expertos en química me corregirán que en condiciones estándar para el agua, la primera etapa de descomposición prácticamente no está implementada (en equilibrio, solo una molécula de 10 7disociados en iones, e incluso entonces no como está escrito). Pero para la descripción de algoritmos, tales esquemas serán ilustrativos. Suponga que una corriente procesa un enlace en una molécula de agua y otra corriente procesa el segundo enlace de la misma molécula. Y sucedió que ambos vínculos deben romperse. Entonces una corriente debería transformar átomos en H + y O - , y la segunda en H + y O 2- . Pero si los flujos hacen esto simultáneamente, en el momento del inicio del procedimiento, el oxígeno está en el estado O y ambos flujos lo convierten en O - , lo cual es incorrecto. Necesitamos prevenir tales situaciones de alguna manera. Diagrama de bloques de una función que maneja un enlace químico:
Comprobamos si los tipos de átomos actuales corresponden al tipo de conexión, si no, luego tomamos de la tabla de tipos por defecto (hay que compilarla con anticipación), luego determinamos el cuadrado de la distancia entre átomos (r 2 ) y, si la conexión implica una longitud máxima o mínima, comprobamos si no salió si estamos más allá de estos límites. Si lo hicimos, entonces debemos cambiar el tipo de conexión o eliminarlo y en ambos casos cambiar los tipos de átomos. Para ello se utilizará la función atomicCAS- comparamos el tipo de átomo actual con el que debería ser y en este caso lo reemplazamos por uno nuevo. Si el tipo del átomo ya ha sido cambiado por otro hilo, volvemos al principio para anular el tipo de enlace. El peor de los casos es si logramos cambiar el tipo del primer átomo, pero no el del segundo. Ya es demasiado tarde para volver atrás, porque después de que cambiamos el primer átomo, otros subprocesos ya podrían hacer algo con él. ¿Dónde está la salida? Sugiero que finjamos que estamos rompiendo / cambiando una conexión de un tipo diferente, y no la que tomamos al principio. Encontramos qué tipo de conexión debería haber existido entre el primer átomo inicial y el segundo modificado y lo procesamos de acuerdo con las mismas reglas que se esperaba originalmente. Si en este caso el tipo de átomo ha cambiado nuevamente, usaremos el mismo esquema nuevamente. Está implícitamente implícito aquí,que un nuevo tipo de enlace tiene las mismas propiedades: se rompe a la misma longitud, etc., y las partículas formadas durante la ruptura son las necesarias. Dado que esta información es tocada por el usuario, de alguna manera le pasamos la responsabilidad de nuestro programa a él, él debe configurar todo correctamente. El código:
__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);
}
En el código, utilicé directivas de preprocesador para permitir verificaciones de situaciones que pueden surgir debido a la supervisión del usuario. Puede apagarlos para acelerar el rendimiento. La función implementa el esquema anterior, pero envuelto en un ciclo más que se ejecuta a través del rango de enlaces de los que este hilo es responsable. En adelante, el identificador del tipo de enlace puede ser negativo, esto significa que el orden de los átomos en el enlace debe invertirse (por ejemplo, el enlace OH y HO es el mismo enlace, pero en el algoritmo el orden es importante, para indicar esto utilizo índices con el opuesto sign), la función invert_bond hace que sea demasiado trivial de describir. Función Delta_periodicaplica condiciones de contorno periódicas para coordinar las diferencias. Si necesitamos cambiar no solo los enlaces, sino también los ángulos de enlace (directiva USE_NEWANG ), entonces necesitamos marcar los átomos para los que hemos cambiado el tipo (más sobre eso más adelante). Para excluir la re-unión de los mismos átomos con un enlace, la matriz principal almacena el índice de uno de los átomos asociados con los datos (esta red de seguridad no funciona en todos los casos, pero para la mía es suficiente). Si rompemos algún tipo de conexión, entonces necesitamos eliminar los índices atómicos correspondientes de la matriz de padres, esto se hace mediante la función 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 función, desafortunadamente, se ejecuta en toda la gama de enlaces. Aprendimos a procesar y eliminar enlaces, ahora necesitamos aprender a crearlos. La siguiente función marca los átomos que son adecuados para formar un enlace químico:
__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 matriz bindBonds almacena información sobre si estos tipos de átomos pueden formar un enlace y, de ser así, cuál. La matriz bindR2 almacena la distancia máxima entre átomos requerida para la unión. Si todas las condiciones son favorables, comprobamos si los átomos de otros vecinos son adecuados para la unión, pero más cercanos.
La información sobre la distancia más cercana al vecino se almacena en la matriz r2Min (por conveniencia, la matriz es de tipo int y los valores se convierten a ella multiplicando por una constante, 100). Si el vecino detectado es el más cercano, recordamos su índice en la matriz neighToBind y establecemos la bandera canBind... Existe un peligro real de que, mientras pasamos a actualizar el índice, otro hilo haya sobrescrito el valor mínimo, pero esto no es tan crítico. Es aconsejable llamar a esta función en funciones que atraviesan pares de átomos, por ejemplo, cell_list o all_pair , descritas en la primera parte . La encuadernación en sí:
__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 función bloquea un átomo, luego el segundo y, si logra bloquear ambos, crea una conexión entre ellos. Al comienzo de la función, los índices de los átomos se ordenan para excluir la situación en la que un hilo bloquea el primer átomo en un par y el otro bloquea el segundo átomo en el mismo par, ambos hilos pasan con éxito la primera comprobación y fallan en el segundo, como resultado de ninguno de los la conexión no los crea. Y finalmente, necesitamos eliminar esos enlaces que marcamos para su eliminación en la función 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;
}
Simplemente movemos los enlaces "anulados" al final de la matriz y disminuimos el número real de enlaces. Desafortunadamente, el código es serial, pero no estoy seguro de si paralelizarlo traerá algún efecto tangible. Las funciones que calculan la energía de enlace real y las fuerzas sobre los átomos, que están indicadas por los campos force_eng de la estructura cudaBond , aún se omiten , pero son completamente análogas a las funciones de los potenciales de par descritas en la primera sección.
Ángulos de valencia
Con los ángulos de valencia, introduciré algunas suposiciones para facilitar los algoritmos y las funciones y, como resultado, serán incluso más simples que para los enlaces de valencia. Primero, los parámetros de los ángulos de enlace deberían depender de los tres átomos, pero aquí asumiremos que el tipo de ángulo de enlace determina exclusivamente el átomo en su vértice. Propongo el algoritmo más simple para formar / eliminar esquinas: siempre que cambiamos el tipo de un átomo, recordamos este hecho en el correspondiente arreglo oldTypes [] . El tamaño de la matriz es igual al número de átomos, inicialmente se llena con -1. Si una función cambia el tipo de un átomo, reemplaza -1 con el índice del tipo original. Para todos los átomos que han cambiado de tipo, elimine todos los ángulos de enlace y recorra todos los enlaces de este átomo para agregar los ángulos correspondientes:
__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;
}
}
La matriz specAngles contiene los identificadores de ángulo de enlace correspondientes al tipo de átomo dado. La siguiente función invoca el cálculo de energía y fuerzas para todos los ángulos:
__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);
}Bueno, por ejemplo, el potencial de tales ángulos, dando una función coseno armónica, que puede indicar el campo force_eng estructura 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;
}
No daré una función para eliminar esquinas "anuladas", no difiere fundamentalmente de clear_bonds .
Ejemplos de
Sin pretender ser exacto, traté de representar el ensamblaje de moléculas de agua a partir de iones individuales. Los potenciales emparejados se establecieron arbitrariamente en forma de potencial de Buckingham, y luego se agregó la capacidad de crear enlaces en forma de potencial armónico, con una distancia de equilibrio igual a la longitud del enlace HO en el agua, 0,96 Å. Además, cuando el segundo protón se unió al oxígeno, se añadió un ángulo de enlace con el vértice del oxígeno. Después de 100.000 pasos, las primeras moléculas aparecieron a partir de iones dispersos al azar. La figura muestra las configuraciones inicial (izquierda) y final (derecha):
Puede configurar un experimento como este: deje que al principio todos los átomos sean iguales, pero cuando están uno al lado del otro forman un enlace. Deje que los átomos unidos formen otro enlace con un átomo libre o con otra molécula unida similar. Como resultado, obtenemos algún tipo de autoorganización, donde los átomos se alinean en cadenas:
Comentarios finales
- Aquí usamos solo un criterio para la unión: la distancia, aunque en principio puede haber otros, por ejemplo, la energía del sistema. En realidad, cuando se forma un enlace químico, por regla general, la energía se libera en forma de calor. Esto no se tiene en cuenta aquí, pero puede intentar implementarlo, por ejemplo, cambiar la velocidad de las partículas.
- Las interacciones entre partículas a través del potencial de enlace químico no niegan el hecho de que las partículas aún pueden interactuar a través de los potenciales de pares intermoleculares y la interacción de Coulomb. Sería posible, por supuesto, no calcular interacciones intermoleculares para átomos unidos, pero esto, en el caso general, requiere comprobaciones prolongadas. Por lo tanto, es más fácil establecer el potencial de enlace químico de tal manera que su suma con otros potenciales dé la función deseada.
- La implementación paralela de la unión de partículas no solo aumenta la velocidad, sino que incluso se ve más realista, ya que las partículas compiten entre sí.
Bueno, hay varios proyectos en Habré que se acercan mucho a este: