ขั้นตอนวิธีการของหวังและลันเดา

ขั้นตอนวิธีของหวังและลันเดา (Wang and Landau algorithm) เป็นขั้นตอนวิธีที่ถูกเพิ่มตามขั้นตอนวิธีการแบบ  Metropolis Monte Carlo ซึ่งถูกเสนอขึ้นโดยฟุเกา หวัง (Fugao Wang) และ เดวิด พี.ลันเดา (David P. Landau) ซึ่งขั้นตอนวิธีการนี้ถูกออกแบบเพื่อให้สามารถคำนวณความหนาแน่นของสภาวะของสารต่าง ๆ ที่ถูกจำลองโดยระบบคณิตกร (คอมพิวเตอร์) เช่น แบบจำลองไอซิ่งของแก้วที่ถูกหมุน หรือ แบบจำลองอะตอมในสนามของแรงในโมเลกุลต่าง ๆ แนวคิดหลักของขั้นตอนวิธีการนี้คือการสร้างความหนาแน่นของสภาวะของสารโดยใช้ความจริงที่ว่า ระบบนั้นเมื่อถูกจำลองแล้ว จะสร้างความหนาแน่นของสภาวะของสารขึ้น ขั้นตอนวิธีการนี้เป็นวิธีการทั่วไปที่สามารถทำให้รู้ได้ถึงความหนาแน่นของสภาวะของสาร ซึ่งเป็นที่ยอมรับโดยทั่วไป
วิธีการแบบ Metropolis Monte Carlo ถูกตั้งขึ้นตามความคิดของโบลท์แมนที่ว่า เมื่อเพิ่มอุณหภูมิให้กับโมเลกุลที่กระจัดกระจายในระหว่างช่วงพลังงานที่สูง หรือเรียกอีกอย่างว่าในสถานะที่ไม่น่าพึงพอใจ กับ ช่วงพลังงานต่ำ หรือเรียกอีกอย่างว่าในสถานะที่น่าพึงพอใจ ซึ่งมีความน่าจะเป็นที่พลังงานนั้นจะมีผลต่อความหนาแน่นของสภาวะของสาร ดังนั้นถ้าหากว่ามีสารใดสารหนึ่งอยู่ในภาวะที่พลังงานต่ำ แล้วได้รับพลังงานมากเพียงพอที่จะทำให้สถานะของสารเปลี่ยนแปลงไป ซึ่งทำให้สภาวะความหนาแน่นของสารที่มีพลังงานมากถูกทำลายไป ยกตัวอย่างเข่นการที่ขี้ผึ้งละลาย ในสภาวะอุณหภูมิต่ำ เฉพาะสถานะของแข็งเท่านั้นที่สามารถจับต้องได้ เมื่ออุณหภูมิเพิ่มขึ้น สถานนะที่ไม่น่าพึงพอใจจะเพิ่มมากขึ้น จนในที่สุดทำให้ของแข็งนั้นหลอมเหลวกลายเป็นของเหลว
ในทฤษฎีของขั้นตอนวิธีการของหวังและลันเดาสามารถนำไปใช้กับระบบใด ๆ ที่ถูกกำหนดสถานะโดยพลังงาน (หรือค่าอื่น ๆ) ในปัจจุบัน ขั้นตอนวิธีการนี้ถูกนำไปใช้กับการแก้ปัญหาปริพันธ์ต่าง ๆ และการหดตัวของโปรตีน
== ขั้นตอนวิธี ==
ขั้นแรก ถ้าเกิดว่าการจัดสรรแบบพลวัตไม่ถูกนำมาใช้ จำนวนน้อยสุดกับมากสุดที่เป็นไปได้ของสภาวะของระบบจะต้องถูกคำนวณก่อน ในขณะเดียวกัน ถ้าหากเราพิจารณาแบบจำลองไอซิ่งในขณะที่อยู่ในสถานะที่น่าพอใจเป็นอย่างยิ่งและอยู่รอบจุดหมุนที่มีทิศทางเดียวกัน และเมื่ออยู่ในสถานะที่น่าพึงพอใจน้อยลงมา จะสมมติการใช้แบบแผนแบบตารางหมากรุก ช่วงของค่าจะถูกนำมาหารเป็นค่าที่แยกจากกันโดยสิ้นเชิงในกราฟแสดงความถี่เชิงสถิติ (Discrete Histogram) เราจะทำการกำหนดค่ามากสุดและน้อยสุดของกระบวนการคำนวณนี้เป็น Emax และ Emin ซึ่งจำเป็นต้องใช้ค่าความต่างระหว่างพลังงานมากสุดและน้อยสุดที่แม่นยำ (Δ) จะได้สมการว่า
N = \frac{E_{\rm max}-E_{\rm min}}{\Delta}
หลังจากนั้นอาเรย์ซึ่งมีจำนวน N ช่องจะต้องเก็บค่าความหนาแน่นของสภาวะของสาร ในขั้นตอนนี้จะต้องใส่ใจเรื่องความถูกต้องของ Δ เพราะว่ายิ่งค่าNเพิ่มขึ้นเวลาที่ต้องใช้ในการจำลองรูปแบบก็จะเพิ่มอย่างรวดเร็ว
ในขั้นต้นค่าความหนาแน่นของสภาวะของสาร g (E) นั้นยังไม่สามารถทราบได้ ดังนั้นค่าทุกช่องในอาเรย์จะต้องถูกตั้งเป็นค่า ๆ หนึ่งซึ่งเป็นค่าปริยาย เนื่องจากค่าความหนาแน่นโดยทั่วไปนั้นจะมีค่าระหว่างพหุคูณของมิติ ดังนั้นเป็นการเก็บข้อมูลควรเก็บเป็นค่าลอการิททึ่มของค่าความหนาแน่นของสภาวะของสาร คือ ln [g (E)] ยิ่งไปกว่านั้น ค่าของกราฟแสดงความถี่เชิงสถิติ H (E) จะต้องถูกทำให้คงที่ไว้
ค่าเริ่มต้นของจุดที่ยังไม่ได้ถูกพบของ ln[g (E)] และ H (E) จะต้องถูกตั้งเป็น0 ในจุดแต่ละจุดถูกแต้มสีเหมือนกับการจำลองแบบMonte-Carlo ซึ่งอาจจะถูกยอมรับโดยพิจารณาจากค่าความหนาแน่นของสภาวะของสารแทนที่จะเป็นการพิจารณาจากเกณฑ์ของMetropolis ซึ่งการเลือกจุดจะถูกยอมรับถ้า
p 
โดยpเป็นค่า ๆ หนึ่งซึ่งมากกว่าเท่ากับ0แต่น้อยกว่า1 Eคือค่าพลังงานของสถานะในปัจจุบัน และE' เป็นค่าพลังงานที่จะต้องเปลี่ยนแปลงไป หลังจากที่การเคลื่อนย้ายของขุดถูกยอมรับหรือถูกปฏิเสธ ค่าการถูกพบในกราฟแสดงความถี่เชิงสถิติ H (E) จะถูกเพิ่มขึ้น1 และค่าความหนาแน่นของสภาวะของสารในกราฟแสดงความถี่เชิงสถิติ H (E) จะถูกคูณโดยค่าคงที่ค่าหนึ่ง
g (E) \rightarrow g (E) \times f
ซึ่งค่าคงที่มักถูกแทนที่ด้วย f=c\simeq 2.72 (ตัวเลขของออยเลอร์) เมื่อกราฟแสดงความถี่เชิงสถิติ H (E) ถูกค้นพบจนทั่วแล้ว g (E) จะมีค่าประมาณที่แม่นยำของค่าความหนาแน่นของสภาวะของสาร ซึ่งอยู่ในขอบเขตของการเปลี่ยนแปลงของค่าคงที่
ในขั้นตอนนี้กราฟแสดงความถี่เชิงสถิติ H (E) จะถูกตั้งค่าเป็น0และค่าการเปลี่ยนแปลงของค่าคงที่จะถูกลดดังนี้
f_{i+1} \rightarrow \sqrt{f_{i}}
การทำขั้นตอนข้างบนซ้ำ ๆ จะทำให้เกิดปัญหาที่ว่า g (E) จะต้องเก็บข้อมูลจำนวนมหาศาลตามที่ได้กล่าวไว้ในเบื้องต้น ถ้าต้องการหลีกเลี่ยงปัญหานี้ ค่าของ
\ln [ g (E) ] \rightarrow \ln [ g (E) ] + \ln [ f ] 
จะต้องถูกนำมาใช้ ซึ่งขั้นตอนนี้ค่า f = e ที่ได้ตั้งไว้ก่อนหน้านี้จะมีประโยชน์เพราะ ln (e) = 1 และค่าคงที่จะถูกเปลี่ยนไปดังนี้
\ln (f_{i+1}) \rightarrow (1/2) \ln (f_{i}) 
ขั้นตอนวิธีการนี้สามารถใช้งานได้เพราะการเปลี่ยนแปลงค่าในขั้นตอนต่าง ๆ นั้นเป็นเหมือนการจำลองค่าความหนาแน่นของสภาวะของสารจริง ๆ ดังนั้นถ้าค่าความหนาแน่นของสภาวะของสารถูกสร้างโดยกราฟแสดงความถี่เชิงสถิติถูกยอมรับได้เมื่อเป็นค่าจำนวนเลขที่กลับกัน ค่าที่ได้จะเป็นค่าประมาณที่มีความแม่นยำ ซึ่งในความเป็นจริงแล้วจริงแล้ว H (E) จะไม่สามารถไปสู่จุดที่ราบอย่างสมบูรณ์แบบได้ ดังนั้นบางเกณฑ์อาจกล่าวว่า ค่าที่บางที่สุดที่เป็นไปได้จะอยู่ในช่วงระหว่าง10%ของจุดสูงสุดและจุดต่ำสุด
การเปลี่ยนแปลงของค่าคงที่โดยใช้รากที่สองอาจทำให้เกิดความผิดพลาดเพราะการอิ่มตัวได้ ดังนั้นวิธีที่จะสามารถหลีกเลี่ยงปัญหานี้ได้คือใช้ค่าคงที่เป็น1 / t โดยที่ tคือเวลาในการจำลอง
ในแบบจำลองง่าย ๆ เช่นแบบจำลองไอซิ่ง ปัญหานี้อาจไม่ได้มีปัญหามากนัก แต่หากเป็นระบบที่ซับซ้อนเช่นโปรตีน เยื่อหุ้มเซล และการหาปริพันธ์หลายมิติจะเป็นปัญหาใหญ่มาก เพราะต้องใช้เวลาคำนวณทีเยอะมาก
การหาค่าความหนาแน่นของสภาวะของสารนั้นไม่ขึ้นต่ออุณหภูมิ แต่ว่าจะสามารถพิจารณาผลลัพธ์เมื่อมีค่าของอุณหภูมิมาเกี่ยวข้องได้
== การทดสอบระบบ ==
เราต้องการหาค่าพลังงานศักย์ของการแกว่งแบบฮาร์โมนิคซึ่งมีสูตรดังนี้
E (x) = x2
จากการพิจารณาค่าความหนาแน่นของสาร ได้ว่า
g (E) =\int\delta (E (x) -E_0) dx=\int\delta (x^2-E_0) dx
เมื่อทำการแก้ปริพันธ์แล้วจะได้
g (E) \propto E^{-1/2}
โดยทั่วไปนั้นการหาค่าความหนาแน่นของสภาวะของสารสำหรับการแกว่งแบบฮาร์โมนิคหลายมิติ จะถูกกำหนดโดยพลังงานE ซึ่งเลขชี้กำลังจะถูกกำหนดตามมิติของระบบ
=== รหัสเทียม ===
รหัสเทียมนี้เป็นรหัสแสดงขั้นตอนวิธีการของหวังและลันเดาในภาษาการเขียนโปรแกรม D
real[real] log_g; // อาเรย์นี้เก็บค่าล็อกการิททึ่มของค่าความหนาแน่นของสภาวะของสาร
uint[real] H; // อาเรย์เก็บค่ากราฟแสดงความถี่ทางสถิติ
real E_old = the_system.energy ; // เก็บค่าพลังงานเก่า
real F = 1; //เก็บค่าคงที่
while (F > epsilon) {
the_system.evolve ; // สร้างการพัฒนาการของระบบ
real E_new = the_system.energy ; // คำนวณค่าพลังงานของระบบปัจจุบัน
if (log_g[E_new] 
รหัสเทียมเป็นตัวอักษร
เริ่มต้นจาก ตั้งค่า f=e
do{
do{
Metropolis แก้ไขค่าซึ่งขึ้นกับความน่าจะเป็น W[i->j] = min[1, r (Ei) / r (Ej)]
ปรับตัว r (E) ในแต่ละขั้น : r (E)  g (i) +f โดย f เรียกว่า modification factor
เมื่อแต่ละ bin ใน histogram ถูก update เป็นจำนวน 10-15 รอบ ค่า f จะถูกลดทอนลงด้วย f (n+1) =1/2f (n) โดยที่ n คือรอบที่ n ของการ update
การ simulate นั้นจะทำไปเรื่อย ๆ จนกว่า modification factor จะน้อยกว่าค่า threshold ที่ได้กำหนดไว้ เช่น fn 0.05จูล
ได้ผลดังนี้
== อื่น ๆ ==
=== การเพิ่มสมรรถภาพของขั้นตอนวิธีของหวังและลันเดาโดยใช้ขอบเขตที่สามารถเปลี่ยนแปลงได้ ===
ขั้นตอนวิธีของหวังและลันเดาในการจำลองระบบขนาดใหญ่ จำเป็นต้องทำการหารช่วงของพลังงานลงเป็น ขอบเขต หลาย ๆ ขอบเขต แล้วนำมารวมกันเพื่อหาคำตอบ ผลของความหนาแน่นของสถานะของสารที่มีระบบขนาดใหญ่มากเกินไป ทำให้เกิดปัญหาการจำลองพลังงานแลททิซของโพลีเมอร์และแบบจำลอง 5ขั้นตอนของ Potts เพราะการจำลองของทั้งสองอย่างนั้นต้องใช้การจำลองที่มีขนาดใหญ่มาก
วิธีการแก้ปัญหาที่เกิดจากระบบมีขนาดใหญ่มากเกินไปทำได้โดยใช้ ขอบเขตที่สามารถเปลี่ยนแปลงได้แทนที่จะใช้ขอบเขตที่มีขนาดคงที่ ค่าของขอบเขตนั้นขึ้นกับค่าของพลังงานในกลุ่มในกราฟแสดงความถี่เชิงสถิติ (Discrete Histogram) มีลักษณะแบนราบหรือไม่เป็นลำดับขั้นของการจำลอง
การเปลี่ยนแปลงขนาดของขอบเขตในแต่ละครั้งของการเปลี่ยนค่าf เราจะทำการตัดขอบของค่า ที่มีผลต่อการจำลองซึ่งใช้ในการจำลองโครงสร้างของขอบเขตปกติ ขอบเขตแบบเปลี่ยนแปลงได้สามารถเพิ่มช่วงของระบบที่จะนำมาจำลองสภานะสภาวะของสารได้อย่างมีประสิทธิภาพ
== ความคิดเห็นของผู้เขียน ==
ขั้นตอนวิธีนี้เป็นขั้นตอนวิธีที่สามารถใช้ประโยชน์ได้มาก นอกจากที่จะใช้ในเชิงฟิสิกส์แล้ว ยังสามารถใช้ในเชิงพาณิชย์ได้ เช่นการทำซอฟต์แวร์ การติดตามการเคลื่อนที่แบบทันทีทันใด สามารถใช้ในการตรวจสอบการเคลื่อนไหวของนักเทนนิส อาจดัดแปลงไปเป็นการทำงานอื่น ๆ ที่มีประโยชน์ต่อมวลมนุษยชาติก็เป็นได้
== อ้างอิง ==
Wang, Fugao and Landau, D. P. (Mar 2001). "Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States". Phys. Rev. Lett. (American Physical Society) 86 (10) : 2050–2053. Bibcode 2001PhRvL..86.2050W. doi:10.1103/PhysRevLett.86.2050. PMID 11289852.
R. E. Belardinelli and S. Manzi and V. D. Pereyra (Dec 2008). "Analysis of the convergence of the 1∕t and Wang-Landau algorithms in the calculation of multidimensional integrals". Phys. Rev. E (American Physical Society) 78: 067701. doi:10.1103/PhysRevE.78.067701.
P. Ojeda and M. Garcia and A. Londono and N.Y. Chen (Feb 2009). "Monte Carlo Simulations of Proteins in Cages: Influence of Confinement on the Stability of Intermediate States". Biophys. Jour. (Biophysical Society) 96 (3) : 1076–1082. Bibcode 2009BpJ....96.1076O. doi:10.1529/biophysj.107.125369.
Belardinelli, R. E. and Pereyra, V. D. (2007). "Wang-Landau algorithm: A theoretical analysis of the saturation of the error". Jour. Chem. Phys. 127 (18) : 184105. doi:10.1063/1.2803061.
http://cv.snu.ac.kr/research/~wlmctracker/index.html
http://www.mendeley.com/research/improving-wanglandau-sampling-with-adaptive-windows/
http://www.comp-phys.org/lugano04/Talks/QWL.pdf 
ฟิสิกส์เชิงคำนวณ
ขั้นตอนวิธี