على سبيل المثال ، الماء:
من الواضح أن الهيدروجين والأكسجين داخل جزيء واحد يتفاعلان بطريقة مختلفة تمامًا عن نفس الأكسجين مع الهيدروجين للجزيء المجاور. وهكذا ، يتم تمييز التفاعلات داخل الجزيئية وبين الجزيئية. يمكن تحديد التفاعلات بين الجزيئات بإمكانيات زوجية قصيرة المدى وكولومب ، والتي تمت مناقشتها في المقالات السابقة. هنا سوف نركز على الجزيئية.
النوع الأكثر شيوعًا للتفاعل داخل الجزيء هو الروابط الكيميائية (التكافؤ). يتم تعيين الروابط الكيميائية من خلال الاعتماد الوظيفي للطاقة الكامنة على المسافة بين الذرات المقيدة ، أي ، في الواقع ، من خلال نفس الزوج المحتمل. ولكن ، على عكس إمكانات الزوج العادية ، لا يتم تحديد هذا التفاعل لأنواع معينة من الجسيمات ، ولكن لزوج معين من الجسيمات (بمؤشراتها). الأشكال الوظيفية الأكثر شيوعًا لإمكانيات الروابط الكيميائية هي الجهود التوافقية:
حيث r هي المسافة بين الجسيمات ، k هو ثابت صلابة الرابطة ، و r 0 هو طول رابطة التوازن ؛ و مورس إمكانات :
حيث D هو عمق البئر المحتمل ، فإن المعلمة α تحدد عرض البئر المحتمل.
النوع التالي من التفاعلات داخل الجزيئية هو زوايا الرابطة. لنلق نظرة على صورة العنوان مرة أخرى. لماذا يتم تصوير الجزيء بزاوية ، لأن القوى الكهروستاتيكية كان من المفترض أن توفر أقصى مسافة بين أيونات الهيدروجين ، والتي تقابل زاوية HOH تساوي 180 درجة؟ الحقيقة هي أنه لا يتم رسم كل شيء في الشكل. من دورة الكيمياء المدرسية ، يمكنك أن تتذكر أن الأكسجين يحتوي على زوجين آخرين من الإلكترونات المنفردة ، والتفاعل معها يشوه الزاوية:
في الديناميات الجزيئية الكلاسيكية ، لا يتم عادةً إدخال كائنات مثل الإلكترونات أو سحب الإلكترون ، لذلك ، لمحاكاة الزوايا "الصحيحة" ، يتم استخدام إمكانات زاوية الرابطة ، أي الاعتماد الوظيفي للطاقة الكامنة على إحداثيات 3 جسيمات. واحدة من أكثر هذه الإمكانات ملاءمة هي جيب التمام التوافقي:
حيث θ هي الزاوية التي شكلها ثلاثي الجسيمات ، و k هي ثابت الصلابة ، و θ 0 هي زاوية التوازن.
هناك إمكانات داخل الجزيئية ذات ترتيب أعلى ، على سبيل المثال ، زوايا الالتواء ، لكنها أكثر اصطناعية من زوايا الرابطة.
إن إضافة تفاعلات بين الجسيمات بمؤشرات محددة مسبقًا أمر تافه. بالنسبة للروابط ، نقوم بتخزين مصفوفة تحتوي على مؤشرات الجسيمات المرتبطة ونوع التفاعل. نعطي كل خيط نطاقه الخاص من الروابط للمعالجة وقد انتهيت. وبالمثل مع زوايا الرابطة. لذلك ، سنقوم على الفور بتعقيد المهمة لأنفسنا: سنضيف القدرة على إنشاء / إزالة روابط كيميائية وزوايا رابطة وقت التشغيل. يأخذنا هذا على الفور من مستوى الديناميكيات الجزيئية الكلاسيكية ويفتح أفقًا جديدًا من الاحتمالات. خلاف ذلك، هل يمكن ببساطة تحميل شيء من الحزم الحالية، على سبيل المثال 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 ، على التوالي ، يعنيان العدد الحالي للروابط والحد الأقصى. ستحتوي مصفوفة السندات على مؤشرات الذرات المراد ربطها ، ونوع الرابطة ووقت تكوين الرابطة (إذا كنا مهتمين فجأة بإحصائيات مثل متوسط عمر الرابطة). ستحمل صفيف الزوايا مؤشرات ثلاثي الذرات التي تشكل زاوية الرابطة ونوع زاوية الرابطة. و 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);
};
يعرّف حقل النوع الشكل الوظيفي للنوع المحتمل ، والنوع الجديد ، والنوع الجديد ، والنوع الجديد ، والنوع الجديد ، وهي مؤشرات لنوع الرابطة وأنواع الذرات التي يجب ربطها بعد تغير الرابطة (تنكسر أو تتحول إلى نوع مختلف من الرابطة). يتم تمثيل هذه الحقول كمصفوفات مع عنصرين. الأول يتوافق مع الموقف عندما يصبح الطول أقصر من r2min 1/2 ، والثاني - عندما يتجاوز r2max 1/2... أصعب جزء في الخوارزمية هو تطبيق خصائص جميع الروابط ، مع الأخذ في الاعتبار إمكانية كسرها وتحويلها ، وكذلك حقيقة أن التدفقات الأخرى يمكن أن تقطع الروابط المجاورة ، مما أدى إلى تغيير في نوع الذرات المقيدة. اسمحوا لي أن أشرح باستخدام مثال نفس الماء. في البداية ، يكون الجزيء متعادل كهربائيًا ، وتتكون الروابط الكيميائية بواسطة إلكترونات مشتركة بين الهيدروجين والأكسجين. بشكل تقريبي ، يمكننا القول أن الشحنات على ذرات الهيدروجين والأكسجين صفرية (في الواقع ، يتم تحويل كثافة الإلكترون إلى الأكسجين ، لذلك ، هناك إضافة صغيرة للهيدروجين ، δ + ، وعلى الأكسجين - 2δ-). إذا كسرنا الرابطة ، سيأخذ الأكسجين أخيرًا إلكترونًا لنفسه ، وسيعطيه الهيدروجين بعيدًا. الجسيمات الناتجة هي H + و O - . في المجموع ، نحصل على 5 أنواع من الجسيمات ، دعنا نسميها تقليديًا: H ، H + ، O ، O- ، يا 2- . يتشكل الأخير إذا فصلنا كلا الهيدروجين من جزيء الماء. وفقًا لذلك ، التفاعلات:
H 2 O -> H + + OH -
and
OH - -> H + + O 2- .
سيصححني خبراء الكيمياء أنه في ظل الظروف القياسية للمياه ، لا يتم تنفيذ المرحلة الأولى من التحلل عمليًا (في حالة التوازن ، جزيء واحد فقط من 10 7تنفصل إلى أيونات ، وحتى ذلك الحين ليس كما هو مكتوب تمامًا). لكن بالنسبة لوصف الخوارزميات ، ستكون هذه المخططات توضيحية. لنفترض أن تيارًا ما يعالج رابطة واحدة في جزيء ماء ، ويعالج تيار آخر الرابطة الثانية لنفس الجزيء. وقد حدث أن كلا الارتباطين بحاجة إلى قطع. ثم يجب أن يحول تيار واحد الذرات إلى H + و O - ، والثاني إلى H + و O 2- . ولكن إذا كانت التدفقات تقوم بذلك في وقت واحد ، في وقت بداية الإجراء ، يكون الأكسجين في حالة O ويقوم كلا التدفقات بتحويله إلى O - ، وهذا غير صحيح. نحن بحاجة لمنع مثل هذه المواقف بطريقة ما. رسم تخطيطي لوظيفة تتعامل مع رابطة كيميائية:
نتحقق مما إذا كانت الأنواع الحالية من الذرات تتوافق مع نوع الاتصال ، وإذا لم يكن الأمر كذلك ، فإننا نأخذ من جدول الأنواع الافتراضية (يجب تجميعها مسبقًا) ، ثم نحدد مربع المسافة بين الذرات (r 2 ) ، وإذا كان الاتصال يعني حدًا أقصى أو أدنى للطول ، فإننا نتحقق مما إذا لم يخرج سواء كنا خارج هذه الحدود. إذا فعلنا ذلك ، فنحن بحاجة إلى تغيير نوع الاتصال أو حذفه وفي كلتا الحالتين تغيير أنواع الذرات. لهذا ، سيتم استخدام الدالة 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 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);
}
في الكود ، استخدمت توجيهات ما قبل المعالج لتمكين عمليات التحقق من المواقف التي قد تنشأ بسبب إشراف المستخدم. يمكنك إيقاف تشغيلها لتسريع الأداء. تقوم الوظيفة بتنفيذ المخطط أعلاه ، ولكنها ملفوفة في حلقة واحدة أخرى تمر عبر نطاق الروابط التي يكون هذا الخيط مسؤولاً عنها. فيما يلي ، يمكن أن يكون معرف نوع الرابطة سالبًا ، وهذا يعني أنه يجب عكس ترتيب الذرات في الرابطة (على سبيل المثال ، رابطة OH و H O هي نفس الرابطة ، ولكن في الخوارزمية الترتيب مهم ، للإشارة إلى ذلك ، أستخدم المؤشرات مع العكس علامة) ، فإن وظيفة invert_bond تجعلها تافهة للغاية بحيث لا يمكن وصفها. دالة Delta_periodicتطبق شروط الحدود الدورية لتنسيق الاختلافات. إذا كنا بحاجة إلى تغيير ليس فقط الروابط ، ولكن أيضًا زوايا الرابطة (توجيه USE_NEWANG ) ، فنحن بحاجة إلى تحديد الذرات التي قمنا بتغيير النوع من أجلها (المزيد حول ذلك لاحقًا). لاستبعاد إعادة ربط نفس الذرات برابطة ، تخزن مجموعة الوالدين فهرس إحدى الذرات المرتبطة بالبيانات (شبكة الأمان هذه لا تعمل في جميع الحالات ، ولكنها كافية بالنسبة لي). إذا قطعنا نوعًا من الاتصال ، فسنحتاج إلى إزالة المؤشرات الذرية المقابلة من مصفوفة الوالدين ، ويتم ذلك عن طريق وظيفة الاستبعاد :
__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). إذا كان الجار المكتشف هو الأقرب ، فإننا نتذكر فهرسه في مصفوفة الجوار وقم بتعيين علامة 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, %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
}
تقوم الوظيفة بحظر ذرة واحدة ، ثم الثانية ، وإذا نجحت في منع كليهما ، فإنها تنشئ اتصالًا بينهما. في بداية الوظيفة ، يتم فرز مؤشرات الذرات من أجل استبعاد الموقف عندما يقوم أحد الخيوط بسد الذرة الأولى في زوج ، بينما يقوم الخيط الآخر بكتلة الذرة الثانية في نفس الزوج ، حيث يجتاز كلا الخيطين الاختبار الأول بنجاح ويفشلان في الثاني ، ونتيجة لذلك ، الاتصال لا يخلقها. وأخيرًا ، نحتاج إلى إزالة تلك الروابط التي حددناها للحذف في وظيفة application_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);
}حسنًا ، على سبيل المثال ، إمكانات مثل هذه الزوايا ، مما يعطي دالة جيب التمام التوافقي ، والتي قد تشير إلى بنية قوة المجال : زاوية الزاوية :
__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 Å. بالإضافة إلى ذلك ، عندما يرتبط البروتون الثاني بالأكسجين ، تمت إضافة زاوية رابطة مع قمة الأكسجين. بعد 100000 خطوة ، ظهرت الجزيئات الأولى من الأيونات المتناثرة بشكل عشوائي. يوضح الشكل التكوينات الأولية (اليسرى) والنهائية (اليمنى):
يمكنك إعداد تجربة مثل هذه: دع الذرات في البداية متماثلة ، ولكن عندما تكون بجانب بعضها البعض ، فإنها تشكل رابطة. دع الذرات المربوطة تشكل رابطة أخرى إما مع ذرة حرة أو مع جزيء مرتبط آخر مشابه. نتيجة لذلك ، نحصل على نوع من التنظيم الذاتي ، حيث تصطف الذرات في سلاسل:
التعليقات النهائية
- استخدمنا هنا معيارًا واحدًا فقط للربط - المسافة ، على الرغم من أنه قد يكون هناك معايير أخرى ، على سبيل المثال ، طاقة النظام. في الواقع ، عندما تتشكل رابطة كيميائية ، كقاعدة عامة ، يتم إطلاق الطاقة في شكل حرارة. لا يؤخذ هذا في الاعتبار هنا ، ولكن يمكنك محاولة تنفيذه ، على سبيل المثال ، تغيير سرعة الجسيمات.
- لا تلغي التفاعلات بين الجسيمات من خلال إمكانات الرابطة الكيميائية حقيقة أن الجسيمات لا يزال بإمكانها التفاعل من خلال إمكانات الزوج بين الجزيئات وتفاعل كولوم. سيكون من الممكن ، بالطبع ، عدم حساب التفاعلات بين الجزيئات للذرات المقيدة ، ولكن هذا ، في الحالة العامة ، يتطلب فحوصات طويلة. لذلك ، من الأسهل ضبط إمكانات الرابطة الكيميائية بحيث يعطي مجموعها مع الإمكانات الأخرى الوظيفة المطلوبة.
- لا يؤدي التنفيذ الموازي لربط الجسيمات إلى زيادة السرعة فحسب ، بل يبدو أيضًا أكثر واقعية ، حيث تتنافس الجسيمات مع بعضها البعض.
حسنًا ، هناك العديد من المشاريع في حبري قريبة جدًا من هذا المشروع: