очевидно, что водород с кислородом внутри одной молекулы взаимодействуют совсем по-другому, нежели тот же кислород с водородом соседней молекулы. Таким образом, разделяют ВНУТРИмолекулярное (intramolecular) и МЕЖмолекулярное (intermolecular) взаимодействие. Межмолекулярное взаимодействие можно задать короткодействующими и Кулоновскими парными потенциалами, о которых речь шла в предыдущих статьях. Здесь же сконцентрируемся на внутримолекулярном.
Наиболее распространённый тип внутримолекулярного взаимодействия – это химические (валентные) связи. Химические связи задаются функциональной зависимостью потенциальной энергии от расстояния между связанными атомами, т.е., по сути, тем же парным потенциалом. Но в отличие от обычных парных потенциалов, данное взаимодействие задается не для определенных типов частиц, а для конкретной пары частиц (по их индексам). Наиболее распространённые функциональные формы для потенциалов химической связи — это гармонический потенциал:
где r – расстояние между частицами, k – константа жесткости связи, а r0 – равновесная длина связи; и потенциал Морзе:
где D – глубина потенциальной ямы, параметр α характеризует ширину потенциальной ямы.
Следующий тип внутримолекулярных взаимодействий – это валентные углы. Обратимся снова к заглавной картинке. Почему молекула изображена уголком, ведь силы электростатики должны были обеспечить максимальное расстояние между ионами водорода, которое соответствует углу H-O-H равному 180°? Дело в том, что на рисунке нарисовано не все. Из школьного курса химии, можно вспомнить, что у кислорода есть ещё 2 неподелённые электронные пары, взаимодействие с которыми и искажает угол:
В классической молекулярной динамике обычно не вводят такие объекты как электроны или электронные облака, поэтому для имитации «правильных» углов используют потенциал валентного угла, т.е. функциональную зависимость потенциальной энергии от координат 3х частиц. Одним из наиболее удобных таких потенциалов является гармонический косинус:
где θ – угол образованной тройкой частиц, k – константа жесткости, а θ0 – равновесный угол.
Существуют внутримолекулярные потенциалы и более высокого порядка, например, торсионные углы, но они ещё более искусственные, чем валентные углы.
Добавить взаимодействия между частицами с заранее известными индексами — дело тривиальное. Для связей храним массив, содержащий индексы связанных частиц и тип взаимодействия. Каждому потоку выдаем свой диапазон связей на обработку и готово. Аналогично с валентными углами. Поэтому мы сразу себе усложним жизнь задачу: добавим возможность создавать/удалять химические связи и валентные углы runtime. Это сразу выводит нас из плоскости классической молекулярной динамики и открывает новый горизонт возможностей. Иначе можно было просто скачать что-нибудь из существующих пакетов, например LAMMPS, DL_POLY или GROMACS , тем более что они распространяются бесплатно.
Ну а теперь немного кода. Добавим соответствующие поля в основную структуру:
//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
Количество связей и углов переменное, но практически всегда можно оценить максимально возможное и выделить память сразу под максимум, чтобы не перевыделять память, поля nBond и mxBond, соответственно, означают текущее число связей и максимальное. Массив bonds будет содержать индексы связываемых атомов, тип связи и время создания связи (если нас вдруг заинтересует такая статистика, как среднее время жизни связи). Массив angles будет хранить индексы тройки атомов, образующих валентный угол, и тип валентного угла. Массивы bondTypes и angleTypes будут содержать характеристики возможных потенциалов валентных связей и углов. Вот их структуры:
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);
};
Поле type определяет функциональный форму потенциала, new_type, new_spec1 и new_spec – это индексы типа связи и типов связываемых атомов, после того как связь изменится (оборвется или превратится в связь другого типа). Эти поля представлены в виде массивов с двумя элементами. Первый соответствует ситуации, когда длина станет короче r2min1/2, второй – когда превысит r2max1/2. Наиболее сложная часть алгоритма – это применение свойств всех связей, с учетом возможности их обрыва и превращения, а также того, что другие потоки могли оборвать соседние связи, что привело к изменению типа связанных атомов. Поясню на примере той же воды. Изначально молекула электронейтральна, химические связи образованы электронами общими для водорода и кислорода. Грубо говоря можно сказать, что заряды на атомах водорода и кислорода нулевые (на самом деле электронная плотность смещена к кислороду, поэтому на водороде небольшой плюс, δ+, а на кислороде – 2δ-). Если мы будем рвать связь, кислород окончательно заберет себе электрон, а водород – его отдаст. Получатся частицы H+ и O—. Итого получается 5 типов частиц, условно обозначим их: H, H+, O, O—, O2-. Последняя образуется, если мы оторвем оба водорода от молекулы воды. Соответственно, реакции:
H2O -> H+ + OH—
и
OH— -> H+ + O2-.
Знатоки химии меня поправят, что при стандартных условиях для воды и первая то стадия распада практически не реализуется (в равновесии, только одна молекула из 107 диссоциирована на ионы, да и то не совсем такие, как написано). Но для описания алгоритмов, такие схемы будут наглядными. Допустим, один поток обрабатывает одну связь в молекуле воды, а другой поток – вторую связь этой же молекулы. И так получилось, что обе связи нужно разорвать. Тогда один поток должен превратить атомы в H+ и O—, а второй в H+ и O2-. Но если потоки это делают одновременно, на момент начала процедуры кислород находится в состоянии O и оба потока его превращают в O—, что неверно. Вот такие ситуации нам необходимо как-то предотвращать. Блок-схема функции, обрабатывающей химическую связь:
Проверяем соответствуют ли текущие типы атомов типу связи, если нет, то берем из таблицы дефолтных типов (она должна быть заранее составлена), далее определяем квадрат расстояния между атомами (r2) и, если связь подразумевает максимальную или минимальную длину, проверяем, не вышли ли мы за эти границы. Если вышли, то нам необходимо поменять тип связи или удалить её и в обоих случаях поменять типы атомов. Для этого будет использована функция atomicCAS – сравниваем текущий тип атома с тем, который должен быть и в этом случае заменяем на новый. Если тип атома уже был изменен другим потоком, возвращаемся в начало, к переопределению типа связи. Наихудший вариант, если тип первого атома мы успели поменять, а второго – нет. Откатываться назад уже поздно, ведь после того, как мы поменяли первый атом, другие потоки уже могли что-то с ним сделать. Какой же выход? Предлагаю притвориться, что мы рвём/изменяем связь другого типа, а не того, за которую взялись вначале. Находим, какой тип связи должен был быть между начальным первым атомом и измененным вторым и обрабатываем её по тем же правилам, что и изначально ожидаемую. Если и в этом случае тип атома снова успел поменяться – снова используем эту же схему. Тут неявно подразумевается, что новый тип связи обладает теми же свойствами – рвется при той же длине и т.д., а частицы, образуемые в ходе разрыва такие, какие нужно. Поскольку эту информацию задет пользователь, мы как бы перекладываем ответственность с нашей программы на него, он должен задать все корректно. Код:
__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 wrongn");
#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 = %dn", 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);
}
В коде я использовал директивы препроцессора, чтобы включить проверки на те ситуации, которые могут возникнуть по недосмотру пользователя. Для ускорения быстродействия их можно отключить. Функция реализует указанную выше схему, но завернута ещё в один цикл, который пробегает по диапазону связей, за который отвечает данный поток. Здесь и далее идентификатор типа связи может быть отрицательный, это означает, что порядок атомов в связи нужно поменять местами (например, связь O-H и H-O – это одна и та же связь, но в алгоритме порядок важен, чтобы это указать я использую индексы с противоположным знаком), делает это функция invert_bond, слишком тривиальная, чтобы её описывать. Функция delta_periodic применяет периодические граничные условия к разнице координат. Если нам нужно менять не только связи, но и валентные углы (директива USE_NEWANG), то необходимо помечать атомы, у которых мы поменяли тип (об этом далее). Для исключения повторного связывания одних и тех же атомов связью, массив parents хранит индекс одного из атомов связанных с данным (эта подстраховка работает не во всех случаях, но для моих этого достаточно). Если мы рвем какую-то связь, то нужно убрать из массива parents соответствующие индексы атомов, это делает функция 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;
}
Функция, к сожалению, пробегает по всему массиву связей. Обрабатывать и удалять связи мы научились, теперь надо научится их создавать. Следующая функция помечает атомы, пригодные для образования химической связи:
__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;
}
}
}
}
В матрице bindBonds хранится информация о том, могут ли данные типы атомов образовывать связь и если могут – то какую. В матрице bindR2 хранится максимальное расстояние между атомами необходимое для связывания. Если все условия благоприятные, то проверяем, нет ли у атомов других соседей пригодных для связывания, но поближе. Информация об ближайшем расстоянии до соседа хранится в массиве r2Min (для удобства массив имеет тип int и значения приводятся к нему с умножением на константу, 100). Если обнаруженный сосед – наиближайший, то запоминаем его индекс в массиве neighToBind и ставим флаг canBind. Есть правда опасность, что пока мы перешли к обновлению индекса, другой поток перезаписал значение минимума, но это не так критично. Данную функцию целесообразно вызывать в функциях, выполняющих обход по парам атомов, например cell_list или all_pair, описанных в первой части. Само связывание:
__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, %dn", 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
}
Функция блокирует один атом, затем второй и, если удалось заблокировать оба, создает связь между ними. В самом начале функции идёт упорядочивание индексов атомов по порядку, чтобы исключить ситуацию, когда один поток блокирует первый атом в паре, а другой – второй атом в той же паре, оба потока успешно проходят первую проверку и заваливаются на второй, в результате ни один из них не создает связь. Ну и, наконец, нам нужно удалять те связи, которые мы пометили на удаление в функции 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;
}
Просто перемещаем «зануленные» связи в конец массива и уменьшаем фактическое количество связей. К сожалению, код серийный, но я не уверен, что его параллелизация принесет ощутимый эффект. Ещё за бортом остались функции, вычисляющие собственно энергию связи и силы на атомах, те на которые указывают поля force_eng структуры cudaBond, но они полностью аналогичны функциям парных потенциалов, описанных в первом разделе.
Валентные углы
С валентными углами я введу несколько допущений, чтобы облегчить алгоритмы и функции, в результате они будут даже проще, чем для валентных связей. Во-первых, параметры валентных углов должны зависеть от всех трех атомов, но здесь мы будем считать, что тип валентного угла определяет исключительно атом в его вершине. Алгоритм образования/удаления углов предлагаю простейший: всегда, когда изменяли тип атома, запоминаем этот факт в соответствующем массиве oldTypes[]. Размер массива равен числу атомов, изначально он заполнен -1. Если какая-то функция меняет тип атома, она заменяет -1 на индекс первоначального типа. У всех атомов, которые поменяли тип, удаляем все валентные углы и пробегаем по всем связям этого атома, чтобы добавить соответствующие углы:
__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;
}
}
Массив specAngles содержит идентификаторы валентных углов, соответствующих данному типу атома. Следующая функция вызывает вычисление энергии и сил для всех углов:
__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);
}Ну и для примера потенциала таких углов, даю функцию гармонического косинуса, на которую может указывать поле force_eng структуры 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;
}
Функцию для удаления «зануленных» углов приводить не буду, она ничем принципиально не отличается от clear_bonds.
Примеры
Не претендуя на точность, я попробовал изобразить сборку молекул воды из отдельных ионов. Парные потенциалы задал произвольно в форме потенциала Букингема, а затем добавил возможность создавать связи в виде гармонического потенциала, с равновесным расстоянием равным длине связи H-O в воде, 0.96 Å. Кроме того, при связывании второго протона с кислородом добавлялся валентный угол с вершиной в кислороде. После 100’000 шагов из случайно разбросанных ионов появились первые молекулы. На рисунке показаны начальная (слева) и конечная (справа) конфигурации:
Можно поставить вот такой эксперимент: пусть вначале все атомы одинаковые, однако оказавшись рядом с друг другом они образуют связь. Позволим связанным атомам образовывать ещё одну связь либо со свободным атомом, либо с ещё одной такой же связанной молекулой. В результате у нас получится вот какая-то такая самоорганизация, где атомы выстраиваются в цепочки:
Заключительные комментарии
1. Здесь для связывания мы использовали лишь один критерий – расстояние, хотя в принципе, могут быть и другие, например энергия системы. В реальности при образовании химической связи, как правило, энергия выделяется в виде теплоты. Здесь это никак не учтено, но можно попробовать это реализовать, например изменять скорости частиц.
2. Взаимодействия между частицами через потенциал химической связи не отменяют того, что частицы могут ещё взаимодействовать через межмолекулярные парные потенциалы и Кулоновское взаимодействие. Можно было бы, конечно, для связанных атомов не вычислять межмолекулярные взаимодействия, но это, в общем случае, требует долгих проверок. Поэтому проще задать потенциал химической связи таким образом, чтобы его сумма с другими потенциалами давала нужную функцию.
3. Параллельная реализация связывания частиц не только дает выигрыш в скорости, но даже выглядит более реалистично, поскольку частицы соревнуются между собой.
Ну и на хабре есть несколько проектов, очень близких к данному:
