มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 3

ต่อจากตอนที่แล้ว http://www.sakngoi.com/?p=238

พอเห็นภาพคร่าวๆของสมการการเคลื่อนที่ที่เราใช้แล้วนะครับ ตอนนี้เราจะมาลองดูเรื่องคุณสมบัติทางเทอโมไดนามิค (อุณหภูมิ, ความดัน) ของระบบที่เรากำลังสนใจครับ เพราะในการคำนวณจริงๆ เราต้องคิดถึงว่าระบบของเรานั้นอยู่ในสภาวะแบบไหน อุณหภูมิคงที่หรือเปล่า ความดันล่ะคงที่ด้วยหรือเปล่า จำนวนของอะตอมคงที่หรือมีการหายไปด้วยอีกหรือเปล่า เพื่อที่ว่าจะได้ดูสมจริงมากที่สุดเมื่อเอาไปเปรียบเทียบกับของจริงในธรรมชาติ

ในโปรแกรม NAMD ที่ผมใช้นี้ การควบคุมในอุณหภูมิคงของระบบ ก็มีอยู่สามวิธีครับ คือ
1) rescaling of velocities
2) weak coupling with a heat bath
3) Langevin dynamics

ในวิธีที่ 1
ก็คือค่อยๆเปลี่ยนอุณหภูมิไปเรื่อยๆที่ล่ะนิด ด้วยการคูณอุณหภูมิของระบบด้วยค่าคงที่ (จริงๆแล้วคูณกับความเร็ว)ครับ
อุณภูมิมันสัมพันธ์กับพลังงานจลน์ของระบบ ดูสมการด้านล่างครับ

\Delta T = \frac{1}{3}\sum_{i=1}^{N}\frac{m_i (\lambda v_i)^2}{Nk_b}

-\frac{1}{3}\sum_{i=1}^{N}\frac{m_i\lambda  v_i^2}{Nk_b}

\Delta T = \left( \lambda^2 - 1\right) T(t)

\lambda  = \sqrt{T_{req}/T(t)}

ใน NAMD จะมีพารามีเตอร์คือ rescaleTemp – อุณหภูมิที่ต้องการ กับ rescaleFreq -ความถี่ที่จะคูณเจ้าค่าคงที่นี้ประมาณว่าจะคูณมันกี่step

วิธีที่ 2

วิธีนี่ก็คล้ายกับวิธีแรก แต่จะมองเหมือนว่าระบบของเราถูกแช่อยู่ในอ่างของความร้อนที่มีอุณหภูมิที่เราต้องการ ระบบของเรากับอ่างนี้ก็จะมีการแลกเปลี่ยนความร้อนกัน โดยที่อัตราการแลกเปลี่ยนความร้อนหรืออุณหภูมิสามารถเขียนได้ว่า

\frac{dT(t)}{dt}=\frac{1}{\tau}\left(T_{bath}-T(t)\right)

หรือในแต่ล่ะstep อุณหภูมิจะต่างกัน

\Delta T  = \frac{\delta t}{\tau}\left(T_{bath}-T(t)\right)

\lambda^2 = 1+ \frac{\delta t}{\tau} \left( \frac{T_{bath}}{T(t)}-1 \right)

วิธีนี้ NAMD มีพารามีเตอร์  tCoupleTemp- ค่าของอุณหภูมิที่ต้องการ (T_{bath}) แต่ต้องมีการเรียก tCouple = on ก่อนนะครับถึงจะใช้ได้

วิธีที่ 3

การควบคุมอุณหภูมิด้วยวิธีที่เรียกว่า Langevin dynamics เป็นอีกวิธีที่นิยมใช้กัน
วิธีนี้จะไม่ได้แก้ไขความเร็วของอะตอมหรือโมเลกุลโดยตรงเหมือนสองวิธีแรก แต่จะมีการเพิ่มแรงนอกเหนือจากแรงที่คำนวณจาก
force field ที่กระทำกับอะตอมหรือโมเลกุลที่สนใจแบบสุ่มเข้าไป (R_i(t))
แล้วลบออกด้วยแรงเสียดทานที่มีค่าสัมประสิทธิ์เสียดทาน หรือ damp constant \eta แรงที่เพิ่มเข้านี้จะเป็นเหมือนการเพิ่มพลังงานหรืออุณหภูมิให้ระบบในขณะที่แรงเสียดทานจะเหมือนการเอาพลังงานที่เพิ่มเข้าไปนี้ออก

m_i\frac{dr_i^2}{dt} = F_{intra}+R_i(t)-\eta m_i\frac{dr_i}{dt}

โดยที่แรงสุ่ม (R_i(t)) จะมาจาก Gaussian distribution ครับ ซึ่งจะมีคุณสมบัติตามนี้

R_i(t)= 0

R_i(t)R_i(t) = 2k_b \eta T_o \deta(t)

โดย T_0 คืออุณหภูมิที่ต้องการครับ 🙂

ถ้าใช้ใน NAMD จะต้องตั้งค่า langevin เป็น on และ
langevinTemp เป็นอุณหภูมิที่ต้องการ ส่วนค่าของ \eta ก็คือ langevinDamping ครับ ลองเล่นดู

ต่อตอนที่ 4