MS計算軌迹的分子間作用能

之前寫過一篇關于計算分子間作用能的文章，那如何計算軌迹的分子間互相作用能呢？此時，需要利用perl腳本進行計算。

MS的help已經給出了計算單個xsd檔案互相作用能的腳本，我們需要做的就是将它進一步優化，用于軌迹檔案分子間作用能的計算。主要思想包括以下幾個部分：（1）提取各個Layer,用于後面的計算（2）計算各個Layer能量及總能量（3）采用循環語句，循環（1）（2）步

1 建立清單 這是所有工作之前的準備工作，需要建立一個清單存儲layer和計算資料

my $doc = $Documents{"Layer.xtd"};
my $newStudyTable = Documents->New("InteractionEnergy.std");
my $calcSheet = $newStudyTable->ActiveSheet;
# 建立清單擡頭
$calcSheet->ColumnHeading(0) = "Cell";
$calcSheet->ColumnHeading(1) = "Total Energy of Cell";
$calcSheet->ColumnHeading(2) = "layer1";
$calcSheet->ColumnHeading(3) = "Energy of layer1";
$calcSheet->ColumnHeading(4) = "layer2";
$calcSheet->ColumnHeading(5) = "Energy of layer2";
           

2 提取各個Layer

my $layer1Doc = Documents->New("bottom.xsd");
my $layer2Doc = Documents->New("top.xsd");
	
$allDoc->CopyFrom($doc);
$layer1Doc->CopyFrom($doc);
$layer1Doc->UnitCell->Sets("Layer 2")->Atoms->Delete; 
$layer2Doc->CopyFrom($doc);
$layer2Doc->UnitCell->Sets("Layer 1")->Atoms->Delete; 
		
#取消所有原子的固定
my $atoms = $allDoc->UnitCell->Atoms;
$atoms->Unfix("XYZ");
my $atoms = $layer1Doc->UnitCell->Atoms;
$atoms->Unfix("XYZ");	
	
#将結構放到表格對應的位置
$calcSheet->Cell($counter-1,0) = $allDoc;
$calcSheet->Cell($counter-1,2) = $layer1Doc;
$calcSheet->Cell($counter-1,4) = $layer2Doc;
           

3 能量計算 首先是基本參數設定，然後進行計算。其中，由于有一步是删除計算後的.xsd檔案，這是因為計算過程中會産生大量的過程檔案，這樣可以防止記憶體不夠。

#首先是forcite的設定，這一步你可以直接打開forcite，設定完所有選項後粘貼過來
my $forcite = Modules->Forcite;
$forcite->ChangeSettings(Settings(Quality => "Fine",
	CurrentForcefield => "COMPASS"));

#計算每個layer的能量及總能量
$forcite->Energy->Run($allDoc);
$calcSheet->Cell($counter-1, 1) = $allDoc->PotentialEnergy; 
$forcite->Energy->Run($layer1Doc);
$calcSheet->Cell($counter-1, 3) = $layer1Doc->PotentialEnergy; 
$forcite->Energy->Run($layer2Doc);
$calcSheet->Cell($counter-1, 5) = $layer2Doc->PotentialEnergy; 

$calcSheet->Cell($counter-1,0) = $allDoc-> UnitCell->Atoms->Delete; 
$calcSheet->Cell($counter-1,2) = $layer1Doc-> UnitCell->Atoms->Delete;
$calcSheet->Cell($counter-1,4) = $layer2Doc-> UnitCell->Atoms->Delete;

my $interactionEnergy = $totalEnergy - ($layer1Energy + $layer2Energy);
$calcSheet->Cell($counter-1, 6) = $interactionEnergy;

#删除過程檔案
my $totalEnergy = $calcSheet->Cell($counter-1, 1);
my $layer1Energy = $calcSheet->Cell($counter-1, 3);
my $layer2Energy = $calcSheet->Cell($counter-1, 5);	
           

4循環計算 以軌迹的幀數為循環次數，設定自己的循環步長，然後将計算主體塞入循環即可

my $numFrames = $doc->Trajectory->NumFrames;
for (my $counter = 1; $counter <= $numFrames ; ++$counter) {
	#sets the current frame to the value of counter
	$doc->Trajectory->CurrentFrame = $counter;
	...循環主體...
}
           

最後提醒大家，help檔案是學習每個軟體最好的教程。對于這類程式的編寫，基本上都可以通過help檔案進行改編，整合為自己所用。也希望大家多多貢獻自己的代碼，互相學習！