李健++馮劍

摘要：將分子模擬數(shù)據(jù)可視化是重要的形象化描述方法。該文利用MATLAB強大的圖形處理能力給出常用的球體和等值面繪制編程方法，并對Gemini表面活性劑自組裝的耗散粒子動力學(xué)模擬給出了具體的結(jié)果。研究表明MATLAB不僅可以很好地繪制了模擬結(jié)果，也大大提高了編程效率。

關(guān)鍵詞：MATLAB；分子模擬；可視化

中圖分類號：TP391 文獻(xiàn)標(biāo)識碼：A 文章編號：1009-3044（2017）30-0181-03

MATLAB Implementation of Visualization of Molecular Simulation Data

LI Jian， FENG Jian

（School of Material Science and Chemical Engineering， Chuzhou University， Chuzhou 239000， China）

Abstract： Visualization of the molecular simulation data is an important method to visualization description. In this paper， MATLAB is been used as a common programming method in the sphere and the isosurface rendering for the powerful graphics capabilities， and self-assembly results of the Gemini surfactant by dissipative particle dynamics simulation are given. The results show that MATLAB are not only good to draw the simulation results， also greatly improve the programming efficiency.

Key words： MATLAB； Molecular Simulation； Visualization

1 概述

計算化學(xué)或分子模擬研究中需要對體系的結(jié)構(gòu)性質(zhì)進(jìn)行繪圖，以便更形象地表示體系的演化或結(jié)構(gòu)性質(zhì)變化。很多計算化學(xué)軟件帶有自己的作圖模塊，如Material Studio。還有一些，較通用的可視化軟件可供選擇，如Rasmol、VMD等。盡管沒有自帶繪圖的計算或模擬程序，如LAMMPS、GROMACS等，一般都有輸出為常用可視化軟件的數(shù)據(jù)格式接口。對于那些自行編制計算或模擬軟件的化學(xué)工作者，一種選擇這些通用繪圖軟件，但在也會有很多有不方便之處，尤其是對于粗粒化的模型或特殊的結(jié)構(gòu)體系。這些化學(xué)工作者也需要編制與研究相關(guān)的可視化軟件。

可視化軟件編制既涉及專門的繪圖語言也涉及具體工作的特定算法，目前主要有兩種方法可供選擇。一是DirextX，它主要面向Windows平臺；另一種是OpenGL，它是跨平臺的編程接口。這些方法的選擇不可避免要涉及很多底層操作，編程人員也需要花費很多精力來學(xué)習(xí)。另外，一些特定的繪圖算法編程難度也很大，如繪制等值面的Marching Cubes算法[1]。

MATLAB是一種集科學(xué)數(shù)據(jù)運算和圖形處理的程序語言，它簡單易用的特性給化學(xué)工作者提供了極大的方便。MATLAB有著強大的繪圖能力，它通過對圖形的線型、立面、色彩、光線以及視角等屬性的處理，將計算數(shù)據(jù)的特性表現(xiàn)得淋漓盡致[2]，在圖像處理方面應(yīng)用廣泛[3]。

本文主要討論利用MATLAB的繪圖接口編制程序繪制表面活性劑分子在水溶液的聚集的膠束形態(tài)。其數(shù)據(jù)來自于耗散粒子動力學(xué)（DPD）模擬[4]。其中主要涉及粒子和等值面繪制。

2 方法

2.1 粒子的繪制

粒子可以是具體的原子，也可以是一個粗?；牧Ｗ?，如DPD模擬中那樣。在分子模擬中繪制粒子是比較基本的操作。在MATLAB首先要創(chuàng)建球模型，調(diào)用語句為：

[x y z] = sphere（num）；

其中，x， y， z是系統(tǒng)創(chuàng)建球體自動分配的值，num是構(gòu)造球體的“細(xì)膩度”。此參數(shù)可以不填，默認(rèn)的是20，意思就是這個球面是20×20個小網(wǎng)格組成的，值越大，球面看上去越真實，但計算量也就越大，繪制速度也越慢。在實際應(yīng)用過程中，可以首先設(shè)置一個較小的值，在匯報工作或撰寫論文時再使用一個更大的值，繪制出精美的圖像。

創(chuàng)建了球模型接下來就需要在坐標(biāo)系中繪制它，調(diào)用的函數(shù)為：

surf（x0+x*r， y0+y*r， z0+z*r， color）；

該函數(shù)表示的是在（x0，y0，z0）處繪制半徑為r，顏色為color的球體。其中的x， y， z就是sphere函數(shù)中的x， y， z值。對于這個函數(shù)的最后一個參數(shù)color，是指定球體顏色的參數(shù)。為了區(qū)分粒子，需要給每種不同粒子類型賦予不同顏色。以下是從一個數(shù)據(jù)文件，讀取數(shù)據(jù)繪制球體的主要代碼，其中數(shù)據(jù)文件由四列組成前三列為粒子的xyz坐標(biāo)，最后一列為粒子類型。另外還需要從一個配置文件中讀取粒子總數(shù)和不同粒子類型的直徑大小。

fidData = fopen（fileData， 'r'）；

%將數(shù)據(jù)文件（n行×4列）讀取到data矩陣中

data = textscan（fidData，'%f %f %f %f'）；

%從data各列依次讀取坐標(biāo)和粒子類型

x = data{1，1}； y = data{1，2}； z = data{1，3}；endprint

type = data{1，4}；

[rX rY rZ] = sphere（）； %創(chuàng)建球體模型

%將各粒子繪制在三維坐標(biāo)軸中

for k=1：total

%type（k）是粒子類型

color = ones（size（X））；

color（：， ：， 1） = type（k）；

if target==type（k）

%在[x（k）， y（k）， z（k）]點作半徑radius（type（k））的球體

surf（radius（type（k）） * rX + x（k）， radius（type（k）） * rY + y（k）， radius（type（k）） * rZ + z（k）， color）；

hold on； %保留當(dāng)前繪圖和坐標(biāo)軸屬性等

end

end

2.2 等值面的繪制

表面活性劑水溶液會形成不同的膠束結(jié)構(gòu)，它有很多表面活性劑分子聚集而成。使用粒子圖，由于膠束外部為各種粒子，不能很清楚地看出聚集體的輪廓，更難以反映不同聚集體的結(jié)構(gòu)差別。為表現(xiàn)聚集體的輪廓圖，常需要繪制其等值面。將粒子圖轉(zhuǎn)化為等值面圖，需要在空間構(gòu)建網(wǎng)格，并將粒子數(shù)據(jù)轉(zhuǎn)換（分散）為網(wǎng)格點的數(shù)值，再通過這些離散的數(shù)值繪制等值面。構(gòu)建網(wǎng)格數(shù)據(jù)需要根據(jù)具體的研究確定特定算法。

本文使用的算法是將模擬盒子分解成等分的小正方體格子，每個小正方體的頂點（即格點）的值初始化為0。當(dāng)具有一定體積的某種類型的粒子占據(jù)該格點時，該格點的值賦值為1。然后將取值大于0.5的格點連接起來作等值面[5]。

在MATLAB中繪制等值面的函數(shù)為isosurface函數(shù)，調(diào)用格式為：

[fc vt] = isosurface（x， y， z， v， value， color）；

x， y， z指的是格點的坐標(biāo)。V為格點的數(shù)據(jù)值，是一個三維矩陣，包含了所有點（x， y， z）的值，且值的多少與配置文件中規(guī)定的坐標(biāo)軸范圍有關(guān)。value指的是將哪個值連接成面，如果指定是0.5，那么isosurface函數(shù)會進(jìn)一步對格點進(jìn)行插值運行獲得所有值為0.5的點，這些點可能處于這些元胞的定點、楞或面上，再將這些值連結(jié)起來繪制成值為0.5的等值面。一般說來，給定具體值的點在元胞上連結(jié)形成三角形，該值的等值面就是有很多個三角形組成，空間元胞分得越多，三角形就越多，等值面就越細(xì)膩，繪制工作量越大，計算開銷越大。color指的是繪制面的顏色，可以用顏色代碼來指定。

如果繪制一個基本的等值面，直接用上述的函數(shù)就足夠了，為了對等值面進(jìn)一步操作，而且光照及顏色效果更佳，則需要使用patch函數(shù)，調(diào)用格式如：

P = patch（‘Faces， fc， ‘Vertices， vt）；

另外，使用函數(shù)set（P， 'FaceColor'， 'green'， 'EdgeColor'， 'none'）可以改變表面的顏色和邊沿的顏色，該語句使用的是表面為綠色（或RGB值[0 1 0]表示），無邊沿。

%構(gòu)建（scaleX scaleY scaleZ）網(wǎng)格矩陣

[X Y Z] = meshgrid（0：1：scaleX， 0：1：scaleY， 0：1：scaleZ）；

[rows cols heis] = size（X）； %獲得剛定義矩陣的長， 寬， 高

V = zeros（rows， cols， heis）； %定義一個初始化為0的數(shù)據(jù)矩陣

%把粒子的坐標(biāo)值轉(zhuǎn)化成正數(shù)的新坐標(biāo)

x = x - rangeMinX； y = y - rangeMinY； z = z - rangeMinZ；

%讀取每個粒子

for i=1：total

%僅對type（i）粒子處理等值面，構(gòu)建頂點值

if targetSur==type（i）

%坐標(biāo)點直接轉(zhuǎn)成格點

x_t = round（x（i）*scaleX/（rangeMaxX-rangeMinX）） + 1；

y_t = round（y（i）*scaleY/（rangeMaxY-rangeMinY）） + 1；

z_t = round（z（i）*scaleZ/（rangeMaxZ-rangeMinZ）） + 1；

V（x_t， y_t， z_t） = 1；

end

end

%繪制等值面

[face vertice] = isosurface（X， Y， Z， V， 0.5）；

P = patch（'Faces'， face， 'Vertices'， vertice）；

set（P， 'FaceColor'， [0 1 0]， 'EdgeColor'， 'none'）； %顏色設(shè)置

axis equal； %坐標(biāo)軸單位相等

view（3）； %設(shè)定三維視點的默認(rèn)值

light（'style'，'infinite'）； %光照設(shè)置

camlight right；

lighting phong；

axis （[0 scaleX 0 scaleY 0 scaleZ]）； %指定坐標(biāo)軸大小

2.3 光照處理

盡管指定了體系中的每種粒子或等值面的顏色，但未必能看到它。就像黑夜里看不到各種物體，更別談物體的顏色了。所以為了顯示繪制的球體等對象，需要設(shè)置光照。它與真實世界一樣，不同顏色光射在物體上，物體的顏色也有差別。只有白色光照反映物體的真實顏色，其它顏色光照，物體反映出來的顏色是由物體自身顏色和投射光顏色的疊加。設(shè)置光源的調(diào)用格式為：endprint

light（‘顏色， 顏色選項， ‘樣式， 樣式選項， ‘位置， 具體位置）；

光源的完整設(shè)置包括三個性質(zhì)（選項），分別是顏色、樣式和位置，次序任意。第一個性質(zhì)是顏色字符，如r表示紅色，g表示綠色。第二個參數(shù)設(shè)置光的類型，分為近光（‘local）和遠(yuǎn)光（‘infinite）兩種類型。當(dāng)選擇近光時，與觀察著相同距離平面上的粒子反光點不一樣選擇遠(yuǎn)光時，這些反光點是一致的。選擇不同的樣式類型，使得繪制的系統(tǒng)更真實。第三個選項，光源所在的位置，填寫三維坐標(biāo)。對于光照的角度也可以使用camlight right這樣的指令，它表示光的照射方向在右邊。

不同光源照射到同一物體有不同效果，如一個很大的平面光源和一個點光源差別就很大。距離平面光源不同位置的粒子其投影大小相同并等于粒子大小，而距離點光源的近的粒子投射陰影則更大。照明效果（算法）設(shè)置為：

lighting <選項>

其中的選項有，flat（平行光，對每個面使用均勻光）、gouraud（計算頂點法向量，再對穿越的每個面進(jìn)行線性插值，用于曲面）和none（關(guān)閉光源）。

另外，物體的材質(zhì)不一樣，反射光不一樣，看到的效果也不一樣。如相同顏色的木制球體和金屬球體，光在木制球體上產(chǎn)生的投射光點具有啞光性質(zhì)、甚至沒有反光點，而在金屬球體上顯得比較明亮。調(diào)用格式為：

material <選項>

其中選項有，shiny（較亮，反射光顏色僅取決于光源顏色）、dull（較暗，沒有鏡面亮點，反射光顏色僅取決于光源顏色）、metal（金屬光澤，反射光顏色取決于光源和對象表面兩者的顏色）和default（默認(rèn)設(shè)置）。

如繪制完粒子代碼后面追加如下光照語句：

light（'style'，'infinite'）；

camlight right；

shading interp； % lighting phong； %等值面使用該語句

3 結(jié)果

3.1 模擬細(xì)節(jié)

本文模擬的為Gemini表面活性劑在水溶液中的自組裝行為。所謂的Gemini表面活性劑就是兩種傳統(tǒng)的表面活性劑的頭基通過橋接基團連結(jié)在一起。本文粗?；腉emini表面活性劑結(jié)構(gòu)為（B3A）3C，其中A是頭基，B3為尾鏈，C為橋接基團。水分子使用W表示。表面活性劑的密度為0.05。模擬盒子為正方體盒子，長寬高均為20。所有相同粒子間的相互作用力參數(shù)為25，aAC=aAW=aWC=25，aAB=aCB=aWB=40。

3.2 粒子與等值面圖

在本文給定的參數(shù)下，Gemini表面活性劑在水中形成了膠束。圖1是僅對表面活性劑聚集體繪圖，表面活性劑的所有基團都使用粒子呈現(xiàn)，為了顯示更為圓滑的球體，sphere中的參數(shù)設(shè)為40。MATLAB繪制粒子圖需要很大的運算量，如果對包含水的所有24000個粒子繪圖，計算機顯卡較弱時是無法繪制或響應(yīng)時間很長。

要想從圖1中看出膠束的輪廓圖是非常困難的。圖2是選取表面活性劑尾鏈B值為0.5時的等值面圖，從該圖可以更清楚看出形成膠束輪廓結(jié)構(gòu)。當(dāng)然通過該圖可以看出膠束的輪廓圖過渡還是比較生硬，很多地方不夠圓滑。則需要對由粒子坐標(biāo)構(gòu)建網(wǎng)格數(shù)據(jù)的算法進(jìn)行調(diào)整。

單純繪制粒子無法看清楚膠束的形態(tài)，僅僅繪制膠束的等值面圖，往往只能就一種類型粒子進(jìn)行。為了既能看清楚膠束輪廓，又能看到某些重要粒子的分布情況需要同時繪制等值面圖和粒子圖。在算法方面只需要將前面介紹的兩種繪制方法的坐標(biāo)軸設(shè)置為相同，并進(jìn)行疊加即可。圖3給出了體系的尾鏈的等值面和橋接基團的粒子圖。本文選擇的力常數(shù)表明橋接基團是親水的，因而它在尾鏈聚集體的外側(cè)，與水交接的界面層中。

4 結(jié)論

本文僅討論了粒子和等值面的繪制方法，其他如化學(xué)鍵的棒狀顯示，以及某官能團的特定顯示方法并沒有提到，當(dāng)這些特定的顯示方法MATLAB均能做到。MATLAB將這些作圖函數(shù)進(jìn)行了高級別封裝，極大降低了科學(xué)可視化的編程工作量?；瘜W(xué)工作者進(jìn)行相應(yīng)研究工作時可以將更多的時間和精力放在特定問題的特殊處理上，從而更清楚地顯示研究結(jié)果。

參考文獻(xiàn)：

[1] Thomas Lewiner ， Hélio Lopes ， Ant?nio Wilson Vieira， Efficient implementation of Marching Cubes cases with topological guarantees[J]， Journal of Graphics Tools， 2003， 8（2）：1-15.

[2] 李麗， 王振領(lǐng). MATLAB工程計算及應(yīng)用[M].北京： 人民郵電出版社， 2001：47-106.

[3] 李卓，李益民. Matlab 與VC++混合編程技術(shù)在圖像處理中的應(yīng)用研究[M].電腦知識與技術(shù)， 2011，7（22）：5450-5452.

[4] Robert D. Groot， Patrick B. Warren. Dissipative particle dynamics： Bridging the gap between atomistic and mesoscopic simulation[J]. J. Chem. Phys. 1997， 107（11）：4423-4435.

[5] Robert D. Groot， Patrick B. Warren. Dynamic simulation of diblock copolymer microphase separation[J]. J. Chem. Phys. 1998， 108（20）：8713-8724.endprint