液體核磁共振中多核體系的量子模擬

黃珊珊，李鵬，何培忠

1.上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海200093；2.上海健康醫(yī)學(xué)院醫(yī)學(xué)影像學(xué)院，上海201318

前言

核磁共振（NMR）技術(shù)是用射頻脈沖激發(fā)處在外磁場(chǎng)中具有自旋的原子核，使之在相鄰能級(jí)間發(fā)生共振躍遷，采樣得到一個(gè)自由感應(yīng)衰減（Free Induction Decay,FID）信號(hào)，經(jīng)過(guò)傅里葉變換得到頻譜圖，用來(lái)測(cè)定物質(zhì)的化學(xué)結(jié)構(gòu)［1］。NMR 譜圖中主要包含3個(gè)信息：化學(xué)位移、耦合常數(shù)、譜線的積分面積，其中耦合常數(shù)是自旋核之間磁相互作用大小的量度，在圖譜中表現(xiàn)為特征譜線裂分形成的多重峰。多核自旋體系中相鄰兩核之間的核磁矩通過(guò)化學(xué)鍵傳遞產(chǎn)生的間接相互作用稱為J 耦合［2］，耦合常數(shù)包含自旋核之間化學(xué)鍵的相關(guān)信息，是NMR 波譜技術(shù)中解析分子幾何結(jié)構(gòu)的重要依據(jù)［3］。液體NMR中分子的隨機(jī)運(yùn)動(dòng)平均了直接耦合作用，J 耦合作用是液態(tài)物質(zhì)特征譜線裂分的主要原因［4］。

量子模擬利用量子力學(xué)原理構(gòu)建量子系統(tǒng)來(lái)模擬實(shí)際實(shí)驗(yàn)，可從模擬數(shù)據(jù)中提取相關(guān)的參數(shù)信息，用于設(shè)計(jì)新的實(shí)驗(yàn)和促進(jìn)理論研究的發(fā)展［5］。核自旋系統(tǒng)是一個(gè)系綜［6］，其量子初態(tài)和整個(gè)演化過(guò)程遵循量子力學(xué)原理，使用數(shù)值計(jì)算方法對(duì)一個(gè)NMR 實(shí)驗(yàn)進(jìn)行量子力學(xué)建模，理論上可得到與實(shí)際實(shí)驗(yàn)一樣的結(jié)果。目前比較成功的NMR 模擬程序有GAMMA［7］、SIMPSON［8］、SPINEVOLUTION［9］，它們的共同點(diǎn)是使用哈密頓算符來(lái)表示NMR 中任一時(shí)刻的全部相互作用，從量子初態(tài)出發(fā)到最后采樣，得到一個(gè)FID信號(hào)，再做傅里葉變換獲得與實(shí)際實(shí)驗(yàn)結(jié)果十分相近的譜。在物理學(xué)原理上，磁共振成像（MRI）和NMR 波譜是一致的，區(qū)別在于MRI 處于低場(chǎng)，而且臨床上只對(duì)質(zhì)子信號(hào)激發(fā)和采樣，因此MRI也可用Bloch 矢量模型［10］進(jìn)行解釋。為了統(tǒng)一這兩個(gè)場(chǎng)景的模擬，Chang等［11］發(fā)布了一個(gè)新的模擬軟件Spin-Scenario，對(duì)NMR 和MRI 進(jìn)行統(tǒng)一的量子力學(xué)模擬。上述模擬程序都能夠模擬許多傳統(tǒng)的NMR實(shí)驗(yàn)，且與實(shí)驗(yàn)結(jié)果高度吻合，在NMR 研究中得到廣泛的應(yīng)用［12］。但目前NMR 模擬最大的不足在于能夠模擬的分子體系都非常小，而且尚無(wú)統(tǒng)一和有效的方法來(lái)解決分子快速運(yùn)動(dòng)問(wèn)題。

自旋量子數(shù)為12 的原子核最適合且最常用于NMR 鑒定研究［13］，本文用量子力學(xué)觀點(diǎn)描述液體NMR 中自旋量子數(shù)為12 的多核體系間的J 耦合作用，并將自旋系統(tǒng)的化學(xué)位移、耦合作用寫入密度矩陣的演化方程，在MATLAB 2015b 軟件平臺(tái)模擬出方程演化結(jié)果以及所對(duì)應(yīng)的NMR 譜圖，整個(gè)計(jì)算過(guò)程完全遵循量子力學(xué)規(guī)律，模擬的結(jié)果是客觀且精確的。

1 J耦合作用的量子描述

具有自旋的原子核在無(wú)外加磁場(chǎng)的環(huán)境中的自旋取向是任意的，磁性相互抵消，不具備宏觀磁性。在外磁場(chǎng)B0中，自旋量子數(shù)為12 的核有了兩種取向［14］，即自旋向上的低能態(tài)和自旋向下的高能態(tài)，兩種狀態(tài)的能量差為ΔE，如圖1所示。

圖1 自旋核能量狀態(tài)Fig.1 Spin nuclear energy state

在NMR 中，與被觀測(cè)核相鄰的原子核的自旋擾亂了附近電子的自旋分布，進(jìn)而擾動(dòng)觀測(cè)核的能級(jí)，由于擾動(dòng)核的自旋具有兩種狀態(tài)，當(dāng)自旋狀態(tài)與觀測(cè)核相同時(shí)觀測(cè)核的能量對(duì)應(yīng)地升高，當(dāng)自旋狀態(tài)與觀測(cè)核相反時(shí)觀測(cè)核能量降低［14］，最終在譜圖上表現(xiàn)為被觀測(cè)核譜線的裂分。

自旋量子數(shù)為12 的一對(duì)耦合原子核組成IS 系統(tǒng)，當(dāng)外磁場(chǎng)B0平行于Z 方向，在垂直于B0的方向上施加角頻率為ω的射頻場(chǎng)B1，系統(tǒng)的哈密頓量為［14］：

其中，ωI、ωS分別為I核和S核的共振頻率；Iz、Sz分別代表I核與S核在Z 軸方向的角動(dòng)量算符；JIS是兩核之間的耦合強(qiáng)度。

系數(shù)c1、c2、c3、c4是疊加系數(shù)，需滿足歸一化關(guān)系：

NMR 是自旋核相鄰能級(jí)之間的能量躍遷，能級(jí)圖如圖2所示［2］。

圖2 雙核自旋體系能級(jí)圖Fig.2 Energy level diagram of a double spin system

對(duì)應(yīng)的能量為：

相鄰能級(jí)之間躍遷的能量差為：

J耦合作用在NMR圖譜中的表現(xiàn)如圖3所示。

圖3 NH4Cl的14N-NMR譜Fig.3 14N-NMR spectrum of NH4Cl

這是NH4Cl 的14N 核磁共振一維譜［15］，耦合常數(shù)在圖中標(biāo)記為1J14N-1H，表示14N與1H之間存在耦合，4個(gè)氫質(zhì)子的耦合作用使得14N 的特征譜線分裂成五重峰。自旋量子數(shù)為12的原子核一維NMR譜的耦合規(guī)律如下［3］：觀測(cè)核與n個(gè)磁等價(jià)的自旋核耦合時(shí)，產(chǎn)生n+ 1條譜線，譜線中心為觀測(cè)核化學(xué)位移特征峰的位置，分裂譜線之間的距離記作耦合常數(shù)，單位為Hz，多重峰之間相對(duì)強(qiáng)度對(duì)應(yīng)于(a+b)n展開(kāi)式的各項(xiàng)系數(shù)。

2 IS體系的密度矩陣

NMR的自旋系綜遵從量子力學(xué)中的劉維爾馮紐曼（Liouville-von Neumann）方程［16］，在模擬前需要將所有算符的矩陣形式都寫出來(lái)。

對(duì)于自旋量子數(shù)為12的單核自旋，用泡利矩陣表示它在各坐標(biāo)方向的角動(dòng)量算符［17］：

對(duì)于雙核自旋系統(tǒng)，角動(dòng)量算符的乘積算子可由單自旋矩陣與單位矩陣的張量積得到［18］，計(jì)算求得IS體系的角動(dòng)量算符如下：

計(jì)算驗(yàn)證上述角動(dòng)量算符滿足歸一化關(guān)系：

其中，E是單位矩陣，與自旋角動(dòng)量算符共同描述密度矩陣的演化。自旋對(duì)的狀態(tài)函數(shù)與系數(shù)組成4×4密度矩陣，矩陣元素提供了關(guān)于系統(tǒng)各種狀態(tài)的信息，在哈密頓的作用下，每個(gè)元素所包含的信息隨時(shí)間演化。

IS體系在Z方向初始狀態(tài)的密度矩陣表示為［18］：

同理可證各矩陣其余元素的狀態(tài)關(guān)系。

密度矩陣隨時(shí)間演化的過(guò)程遵循如下演化過(guò)程［19］：

其中，σ(t)代表t時(shí)刻密度矩陣的狀態(tài)，R(τ)是時(shí)間演化超算符，在式中表示哈密頓在τ時(shí)間段內(nèi)對(duì)密度矩陣的作用。

密度矩陣包含自旋體系演化狀態(tài)的全部信息，但是能被儀器檢測(cè)到的只有單量子躍遷過(guò)程，I核與S核的信號(hào)分別用觀測(cè)算符I+、S+來(lái)檢測(cè)［14］：

在哈密頓的作用下密度矩陣的對(duì)角元代表自旋系統(tǒng)在對(duì)應(yīng)本征態(tài)處的被發(fā)現(xiàn)的概率，通過(guò)計(jì)算密度矩陣與觀測(cè)算符乘積的跡得到矩陣中特定的乘積算子，作為磁共振的信號(hào)。

3 多核體系的波譜模擬

3.1 雙核體系波譜模擬

NMR 中，IS 體系在熱平衡態(tài)時(shí)的密度矩陣記為σ0：

因IS 為異核體系，沿X 軸正方向?qū)υ撓到y(tǒng)施加90°脈沖激發(fā)I核共振，I核的激勵(lì)脈沖對(duì)S核的密度矩陣無(wú)影響，此時(shí)IS體系的密度矩陣演化記作σ1：

表示I核的磁化矢量在脈沖的作用下發(fā)生90°翻轉(zhuǎn)，倒向-Y 軸，同時(shí)IS體系發(fā)生耦合作用，此時(shí)系統(tǒng)的哈密頓量為：

其中，ωI、ωS分別為I核與S核的共振頻率；JIS為兩核之間的耦合常數(shù)。在化學(xué)位移和脈沖的作用下I核開(kāi)始自由進(jìn)動(dòng)和弛豫，角動(dòng)量算符在空間中發(fā)生旋轉(zhuǎn)，密度矩陣在哈密頓作用下演化的方程ρ()t表示為：

密度矩陣用I核的觀測(cè)算符檢測(cè)，此時(shí)得到的共振信號(hào)表示為：

式中tr表示求密度矩陣的跡。

考慮到液體中FID信號(hào)的衰減為洛侖茲線型，因此需對(duì)共振信號(hào)加上e指數(shù)窗函數(shù)模擬宏觀信號(hào)［20］，最終的時(shí)域信號(hào)為：

以上就是IS 體系在單脈沖核磁共振模擬實(shí)驗(yàn)中密度矩陣演化的全過(guò)程。最終的FID 信號(hào)可通過(guò)設(shè)置采樣參數(shù)進(jìn)行信號(hào)采樣，經(jīng)過(guò)傅里葉變換得到模擬譜圖。

以IS 異核體系的I核譜圖為例，I核化學(xué)位移為400 Hz，耦合常數(shù)JIS=200 Hz，采樣間隔為1 ms，采樣時(shí)間t=1 s，窗函數(shù)的時(shí)間T=0.15 s。在MATLAB 2015b 軟件平臺(tái)進(jìn)行數(shù)值模擬，得到的時(shí)域信號(hào)如圖4所示。

信號(hào)進(jìn)行傅里葉變換之后得到的模擬譜如圖5所示。

從譜圖中可知I核的單峰裂分為強(qiáng)度相等的二重峰，測(cè)得分裂峰之間的距離為200 Hz，波峰的中心為I核的化學(xué)位移400 Hz，模擬結(jié)果與耦合關(guān)系一致。

如果IS 是同核體系，則90°激勵(lì)脈沖對(duì)兩核均起作用，密度矩陣表示為：

圖4 I核的FID信號(hào)Fig.4 Free induction decay(FID)signal of nucleus I

圖5 I核的頻譜Fig.5 Spectrum of nucleus I

I核與S核的磁化矢量均倒向-Y 軸，在耦合作用下系統(tǒng)的哈密頓與異核IS 體系保持一致，此時(shí)密度矩陣的演化方程：

使用兩核的觀測(cè)算符，最后共振信號(hào)可表示為：

IS同核體系中，兩個(gè)核處于不同的化學(xué)環(huán)境，I核的化學(xué)位移為200 Hz，S 核的化學(xué)位移為600 Hz，耦合常數(shù)JIS=30 Hz，采樣參數(shù)設(shè)置與IS 異核體系波譜模擬一致，得到的時(shí)域信號(hào)如圖6所示。信號(hào)經(jīng)過(guò)傅里葉變換之后得到的模擬譜圖如圖7所示。

圖6 IS核的FID信號(hào)Fig.6 FID signal of nucleus IS

圖7 IS核的頻譜Fig.7 Spectrum of nucleus IS

分析譜圖可知，I核的特征峰中心在200 Hz 處，譜線裂分為強(qiáng)度相等的二重峰，波峰間隔為30 Hz；S核的特征峰中心在600 Hz 處，譜線也裂分為強(qiáng)度相等的二重峰，波峰間隔為30 Hz，模擬結(jié)果與耦合關(guān)系的描述一致。

3.2 三核體系波譜模擬

將自旋量子數(shù)為12 的3 個(gè)耦合原子核組成的體系記作I2S，其中I1、I2核是為自然豐度高的兩個(gè)同類核，S核是自然豐度低的稀核。做S核的一維NMR譜時(shí)，I1、I2核對(duì)S核均有耦合作用，耦合常數(shù)分別記作J1S、J2S；I1、I2之間也存在耦合關(guān)系，耦合常數(shù)記作J12；3個(gè)核的共振頻率分別記為ωI1、ωI2、ωS。

由IS 耦合體系的自旋算符推導(dǎo)，同理可得I2S體系各核的自旋算符矩陣，例如I1核各方向的角動(dòng)量算符的矩陣分別表示為：

對(duì)I2S體系的S核作波譜模擬，沿X軸正方向?qū)核施加90°脈沖，此時(shí)系統(tǒng)的哈密頓函數(shù)為［13］：

重復(fù)上文雙核體系NMR 模擬的演化步驟可得到模擬信號(hào)。S核的化學(xué)位移為400 Hz，耦合常數(shù)J1S=J2S=100 Hz、J12=30 Hz；采樣間隔為1 ms，采樣時(shí)間設(shè)為t=1 s，窗函數(shù)的時(shí)間T=0.2 s。得到的時(shí)域信號(hào)如圖8所示。

圖8 S核的FID信號(hào)Fig.8 FID signal of nucleus S

信號(hào)進(jìn)行傅里葉變換之后得到的模擬波譜如圖9所示。

圖9 S核的頻譜Fig.9 Spectrum of nucleus S

圖譜中，S核的特征峰中心在400 Hz 處，特征譜峰分裂為強(qiáng)度比為1:2:1 的譜線，分裂譜線與中心距離為100 Hz，模擬結(jié)果與描述的耦合關(guān)系一致，譜峰強(qiáng)度之比符合三核體系譜線耦合裂分的規(guī)律。

4 結(jié)束語(yǔ)

本文使用量子力學(xué)數(shù)值模擬算法，模擬了液體NMR 中多核體系的自旋動(dòng)力學(xué)演化。首先將液體NMR 中自旋體系所涉及的相互作用寫入哈密頓函數(shù)，并將哈密頓量作用于密度矩陣演化方程，精確記錄密度矩陣在演化過(guò)程中的狀態(tài)，采樣所得的時(shí)域信號(hào)經(jīng)過(guò)傅里葉變換得到譜圖，分別模擬了兩核和三核體系中的液體核磁共振實(shí)驗(yàn)，其結(jié)果與對(duì)應(yīng)的理論預(yù)測(cè)及預(yù)期的實(shí)驗(yàn)結(jié)果一致。這種算法可推廣到更高數(shù)量的自旋體系，可幫助研究人員理解NMR的物理學(xué)原理，也可以用于輔助分析真實(shí)的NMR實(shí)驗(yàn)。