分子模擬數(shù)據(jù)可視化的MATLAB實(shí)現(xiàn)

李健++馮劍

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

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

中圖分類(lèi)號(hào)：TP391 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：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 概述

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

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

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

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

2 方法

2.1 粒子的繪制

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

[x y z] = sphere（num）；

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

創(chuàng)建了球模型接下來(lái)就需要在坐標(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值。對(duì)于這個(gè)函數(shù)的最后一個(gè)參數(shù)color，是指定球體顏色的參數(shù)。為了區(qū)分粒子，需要給每種不同粒子類(lèi)型賦予不同顏色。以下是從一個(gè)數(shù)據(jù)文件，讀取數(shù)據(jù)繪制球體的主要代碼，其中數(shù)據(jù)文件由四列組成前三列為粒子的xyz坐標(biāo)，最后一列為粒子類(lèi)型。另外還需要從一個(gè)配置文件中讀取粒子總數(shù)和不同粒子類(lèi)型的直徑大小。

fidData = fopen（fileData， 'r'）；

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

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

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

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）是粒子類(lèi)型

color = ones（size（X））；

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

if target==type（k）

%在[x（k）， y（k）， z（k）]點(diǎn)作半徑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 等值面的繪制

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

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

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

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

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

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

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

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

%構(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）； %獲得剛定義矩陣的長(zhǎng)， 寬， 高

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

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

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

%讀取每個(gè)粒子

for i=1：total

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

if targetSur==type（i）

%坐標(biāo)點(diǎn)直接轉(zhuǎn)成格點(diǎ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è)定三維視點(diǎn)的默認(rèn)值

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

camlight right；

lighting phong；

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

2.3 光照處理

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

light（‘顏色， 顏色選項(xiàng)， ‘樣式， 樣式選項(xiàng)， ‘位置， 具體位置）；

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

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

lighting <選項(xiàng)>

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

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

material <選項(xiàng)>

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

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

light（'style'，'infinite'）；

camlight right；

shading interp； % lighting phong； %等值面使用該語(yǔ)句

3 結(jié)果

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

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

3.2 粒子與等值面圖

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

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

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

4 結(jié)論

本文僅討論了粒子和等值面的繪制方法，其他如化學(xué)鍵的棒狀顯示，以及某官能團(tuán)的特定顯示方法并沒(méi)有提到，當(dāng)這些特定的顯示方法MATLAB均能做到。MATLAB將這些作圖函數(shù)進(jìn)行了高級(jí)別封裝，極大降低了科學(xué)可視化的編程工作量?；瘜W(xué)工作者進(jìn)行相應(yīng)研究工作時(shí)可以將更多的時(shí)間和精力放在特定問(wèn)題的特殊處理上，從而更清楚地顯示研究結(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工程計(jì)算及應(yīng)用[M].北京： 人民郵電出版社， 2001：47-106.

[3] 李卓，李益民. Matlab 與VC++混合編程技術(shù)在圖像處理中的應(yīng)用研究[M].電腦知識(shí)與技術(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