界面张力计算脚本（配合Forcite使用）

本期推送我在MS官方售后论坛里面淘来的界面张力计算脚本

本文的脚本，来自于https://bioviacommunity.force.com/在不登陆账号的情况下，可以自由下载查看的内容

计算界面张力的脚本，脚本作者是Jason DeJoannis，就是写了xlink脚本的那位大佬。对xlink交联脚本感兴趣的同学可以点击Forcite交联脚本

该脚本需要用户自定义的数据在脚本中有着明显的提示，这里就不多赘述了。针对脚本使用前的体系该如何模拟，Jason在论坛中给出了自己的回答，这里很有必要跟大家说一下。

这里可以计算液-液界面和液-真空界面的界面张力，使用的公式在前面的注释部分给出了，如果需要引用文献可以引用J Chem Phys B, 2007, 11, 7812; Molecular Simulation, 2007, 33, 27

*如果计算液-真空界面，就像下图一样。一定要使用NVT系综进行计算。为啥不能用NPT？这个大家建议建模如下图然后尝试一样NPT系综模拟，这样你得到的宝贵经验会指导你以后的模拟计算。不想亲自模拟的，可以在本文的角落处找到我给出的问题答案。

*如果是研究液-液界面，那么需要先进行NPT作为平衡过程然后NVT生产过程收集轨迹文件。

无论哪种界面，都要注意，要保持层层分明，一不小心搞成乳液了就没办法算了。

#!perl

use strict;
use MaterialsScript qw(:all);

###################################################################
#
# INTERFACIAL TENSION
#
# Author:   Jason DeJoannis
# Date:   Sept. 23, 2012
#
# Calculates the average interfacial tension from an atomistic
# or Mesocite trajectory using the stress tensor stored in
# each frame (Stress->Eij). The assumption is that there
# are two interfaces and the box size is fixed.
#
#    gamma    =  0.5 L < P_normal - P_lateral >
#    L     = box length in direction normal to surfaces
#    P_normal    = pressure in normal direction
#    P_lateral   = pressure in lateral direction
#
########## USER DEFINED ###########################################

# Surface normal direction X, Y or Z
my $surfnormal = "Z";

# Trajectory with stress tensor
my $doc = $Documents{"water slab.xtd"}; 

# Begin with this frame number
my $startFrame = 1;

###################################################################

# Results study table
my $std = Documents->New("SurfTens.std");
$std->ColumnHeading(0) = "Time (ps)";
$std->ColumnHeading(1) = "P_norm (GPa)";
$std->ColumnHeading(2) = "P_lat (GPa)";
$std->ColumnHeading(3) = "Surface Tension (dyne/cm)";

# Initialize
my $Pn = 0; my $Pl = 0; my $surftens = 0;
my $Lnorm; my $axisN; my $axisL1; my $axisL2;
if    ($surfnormal eq "X") 
{ 
  $Lnorm = $doc->Lattice3D->LengthA; 
  $axisN = 1; $axisL1 = 2; $axisL2 = 3;
} 
elsif ($surfnormal eq "Y") { 
  $Lnorm = $doc->Lattice3D->LengthB; 
  $axisN = 2; $axisL1 = 1; $axisL2 = 3;
} 
elsif ($surfnormal eq "Z") { 
  $Lnorm = $doc->Lattice3D->LengthC; 
  $axisN = 3; $axisL1 = 1; $axisL2 = 2;
}
else { die "Surface normal must be one of X, Y or Z\n"; }
my $trj = $doc->Trajectory;
my $nframes = $trj->NumFrames;
my $frameCount = 0;

# Loop over frames, and use the stress tensor in each
for (my $i=$startFrame; $i<=$nframes; $i++) 
{
  $trj->CurrentFrame = $i;

  # Instantaneous normal (Pn_i) and lateral (Pl_i) pressure components
  my $Pn_i = -$doc->SymmetrySystem->Stress->Eij($axisN,$axisN);
  my $Pl_i = -0.5 * ( $doc->SymmetrySystem->Stress->Eij($axisL1,$axisL1)
               +$doc->SymmetrySystem->Stress->Eij($axisL2,$axisL2) );

  # Surface tension assuming two interfaces and converting to dyne/cm
  my $surftens_i = 100.0*$Lnorm*0.5*($Pn_i-$Pl_i);
  
  # Report raw data to study table
  $std->Cell($frameCount,0) = $trj->FrameTime;
  $std->Cell($frameCount,1) = $Pn_i;
  $std->Cell($frameCount,2) = $Pl_i;
  $std->Cell($frameCount,3) = $surftens_i;

  # Averaging
  $Pn += $Pn_i;
  $Pl += $Pl_i;
  $surftens += $surftens_i;
  $frameCount++;

}

# Normalize and report averages
$Pn /= $frameCount;
$Pl /= $frameCount;
$surftens /= $frameCount;
$std->Title = "Raw Data";
$std->InsertSheet(0, "Averages");
$std->ColumnHeading(1) = "P_norm (GPa)";
$std->ColumnHeading(2) = "P_lat (GPa)";
$std->ColumnHeading(3) = "Surface Tension (dyne/cm)";
$std->Cell(0,0) = "Avg";
$std->Cell(0,1) = $Pn;
$std->Cell(0,2) = $Pl;
$std->Cell(0,3) = $surftens;

END

往期推荐

液-真空体系做NPT平衡，最终真空层会消失。如果忘记设置压力数值，会出现跑飞了的情况。

感谢与热爱计算的你相遇↓↓↓