Hi all,

            I performed a MD about the interaction analysis between ionic liquid and metal oxide solid Al2O3, ionic liquid is a confined layer, metal oxide is the other layer. when I finished the geometry optimization to relax the cell, I continue to equilibrate the cell using NPT, then did NVT to do the production, finally, I want to get the concentration profile for ionic liquid as the distance goes away from the interface. Now when I ran the concentration profile, I can only get the relative concentration which has no units for Y-axis, I want to convert it to atomic density profile which has the unit g/cc, that will be more easy to interpret, I just get a script for the accelrys's technician supporter, Dr. Jason, but I am really a new one, I don't know how to write or revise a script, so is there anyone who knows that can tell me how to revise this script?

here is the script , also you can get it from the attachment.

what I need to do is that I need to revise the user input part? and secondly, how can I get the number of bins when I did'nt specify it.

Hoping for anyone can help me~~~

#!perl

use strict;
use Getopt::Long;
use MaterialsScript qw(:all);

# Calculate the density profile of a given species across one spatial direction of the cell
# Author:  Jason DeJoannis
# Date:  Feb. 15, 2013

# The Forcite | Analaysis | ConcentrationProfile method is used to generate the relative
# number concentration of atoms of each type. These are converted to a mass density and added up.
# For trajectories with variable volume, the average bin volume is used.

# User inputs
my \\\$doc  = \\\$Documents{"alumina_surf_300K.xtd"};
my \\\$direction = "C";
my \\\$nbins = 20;
my \\\$FrameRange = "ALL";  # n-m or ALL
my @sets = ("BMIM+", "PF6-", "Alumina");

# Forcite settings
my \\\$H = 0; my \\\$K = 0; my \\\$L = 0;
die "Direction can only be one of A, B, C" unless \\\$direction =~ /^[ABC]\\\$/;
\\\$H = 1 if (\\\$direction eq "A");
\\\$K = 1 if (\\\$direction eq "B");
\\\$L = 1 if (\\\$direction eq "C");
my \\\$Forcite = Modules->Forcite;
\\\$Forcite->ChangeSettings([
ConcentrationProfileUserNumberOfBins => "Yes",
ConcentrationProfileNumberOfBins => \\\$nbins,
ConcentrationProfileSetA  => "Concentration analysis",
ConcentrationProfileSpecifiedDirection => "Yes",
ConcentrationProfileUserH  => \\\$H,
ConcentrationProfileUserK  => \\\$K,
ConcentrationProfileUserL  => \\\$L,
ActiveDocumentFrameRange  => \\\$FrameRange
]);

# Find average bin volume of trajectory
my \\\$start = 1;
my \\\$end = \\\$doc->Trajectory->NumFrames;
if (\\\$FrameRange =~ m/^(\\d+)-(\\d+)\\\$/)
{
\\\$start = \\\$1;
\\\$end = \\\$2;
}
my \\\$volume = \\\$doc->SymmetrySystem->Volume;
if (\\\$doc->Type eq "3DAtomisticTrajectory")
{
\\\$volume = 0;
for (my \\\$frame=\\\$start; \\\$frame<=\\\$end; \\\$frame++)
{
  \\\$doc->Trajectory->CurrentFrame = \\\$frame;
  \\\$volume += \\\$doc->SymmetrySystem->Volume;
}
\\\$volume /= \\\$doc->Trajectory->NumFrames;
}
my \\\$binvol = \\\$volume / \\\$nbins;

# Main loop over each set for analysis

my \\\$setnr=1;
my \\\$output;

foreach my \\\$setname (@sets)
{
my \\\$set;
eval { \\\$set = \\\$doc->UnitCell->Sets(\\\$setname); };
if (\\\$@)
{
  warn "WARNING: Unable to open set \\\$setname\\n";
  next;
}
print "Found set \\\$setname\\n";

# Find elements in the set
my %elements;
foreach my \\\$atom (@{\\\$set->Atoms})
{
  \\\$elements{\\\$atom->ElementSymbol} = 1;
}
my @elem = keys %elements;

# Divine element properties
my %count;
my %mass;
foreach my \\\$element (@elem)
{
  \\\$count{\\\$element}=0;
  foreach my \\\$atom (@{\\\$set->Atoms})
  {
   if (\\\$atom->ElementSymbol eq \\\$element)
   {
    \\\$count{\\\$element}++;
    if (\\\$count{\\\$element} == 1)
    {
     \\\$mass{\\\$element} = \\\$atom->Mass;
    }
   }
  }
}

# Print for validation
foreach my \\\$element (keys %count)
{
  printf "Element %s atoms %d mass %f\\n", \\\$element, \\\$count{\\\$element}, \\\$mass{\\\$element};
}

# Initialize density array
my @density; my @distance;
for (my \\\$i=0; \\\$i<\\\$nbins; \\\$i++)
{
  \\\$density[\\\$i]=0;
}

# Analyze concentration for each element
foreach my \\\$element (@elem)
{
  my @atoms;
  foreach my \\\$atom (@{\\\$set->Atoms})
  {
   if (\\\$atom->ElementSymbol eq \\\$element)
   {
    push(@atoms, \\\$atom);
   }
  }
  my \\\$tmpSet = \\\$doc->CreateSet("Concentration analysis", \\@atoms);
  my \\\$results = \\\$Forcite->Analysis->ConcentrationProfile(\\\$doc);
  \\\$results->ConcentrationProfileChart->Delete;
  my \\\$std = \\\$results->ConcentrationProfileChartAsStudyTable;
  \\\$std->Name = \\\$element;
  \\\$std->ColumnHeading(2) = "Atoms";
  \\\$std->ColumnHeading(3) = "Density (g/cc)";
  my \\\$navg = \\\$count{\\\$element}/\\\$nbins;
  my \\\$densfac = 1.66053886 * \\\$mass{\\\$element} / \\\$binvol;
  for (my \\\$row=0; \\\$row<\\\$std->RowCount; \\\$row++)
  {
   \\\$std->Cell(\\\$row,2) = \\\$std->Cell(\\\$row,1) * \\\$navg;
   \\\$std->Cell(\\\$row,3) = \\\$std->Cell(\\\$row,2) * \\\$densfac;
   \\\$density[\\\$row] += \\\$std->Cell(\\\$row,3);
   push @distance, \\\$std->Cell(\\\$row,0);
  }
  \\\$tmpSet->Delete;
  \\\$std->Delete;
}

# Put results in a new study table
if (\\\$setnr == 1)
{
  \\\$output = Documents->New(\\\$doc->Name." Density Profile.std");
  \\\$output->ColumnHeading(0) = "Distance (A)";
}
\\\$output->ColumnHeading(\\\$setnr) = "\\\$setname Density (g/cc)";
for (my \\\$row=0; \\\$row<@density; \\\$row++)
{
  \\\$output->Cell(\\\$row,0) = \\\$distance[\\\$row] if (\\\$setnr == 1);
  \\\$output->Cell(\\\$row,\\\$setnr) = \\\$density[\\\$row];
}

\\\$setnr++;
}