Hello, everyone,

I'm a new user of Material Studio 5.5, and I have read many posts or scripts here on different systems by different modules which are very useful for me.

Recently I worked on the micellization of surfactant(SDS) in solution, and want to find out the critical micellar concentration of the surfactant. Note that there are different size of aggregation of SDS, we must calculate the number of free surfactants( which are defined as those surfactants existing alone in the solution and do not belong to any aggregate) and find out how to define a micelle. And from literature, most researchers distinguished aggregates using a criterion to decide whether two surfactants are considered to belong to a same cluster. Now I tried to define a Radius Rc, if the distance between two particles is less than Rc, then these two particles belong to the same cluster.

Though I found a similar script to calculate the number of molecules in a layer, but I don't know how to edit it for my specific system.Could you be kind enough to help me find out if there are any errors in the script below? I got a trajectory file "61 beads of SDS.xtd" after Mesocite DPD simulation.The molecule SDS is modeled as H1T4, and water W1, where H, T,W represent head bead, tail bead of SDS and water bead, respectively.

#!perl

# This script calculates the number of molecules in the micelle. The molecule is counted in only if
#the distance between the molecular center and any bead of the targeted moelcule is less than
#the cutoff distance

use strict;
use MaterialsScript qw(:all);

##################################################################
#Begin user editable settings

#Starting file
my \\\$docstart="61beads of SDS.xtd";

#Study Table file
my \\\$studytable="number.std";

#molecule name to check distance with
my \\\$moleculename = "SDS";

#Cutoff distance
my \\\$cutoff=10.0;

#End user editable settings
##################################################################

#The trajectory file
my \\\$doc = \\\$Documents{\\\$docstart};

my \\\$lengthX = \\\$doc->Lattice3D->LengthA;
my \\\$lengthY = \\\$doc->Lattice3D->LengthB;
my \\\$lengthZ = \\\$doc->Lattice3D->LengthC;

#The study table file
my \\\$newStudyTable = Documents->New(\\\$studytable);
my \\\$calcSheet = \\\$newStudyTable->ActiveSheet;

#Specify the heads of the colums in the study table
\\\$calcSheet->ColumnHeading(0) = "time";
\\\$calcSheet->ColumnHeading(1) = "number of molecules";

#Find the total number of frames
my \\\$numFrames = \\\$doc->Trajectory->NumFrames;
print "Number of frames being analyzed is \\\$numFrames\\n";

#This option is not working
#my \\\$SDSMolecules = \\\$doc->UnitCell->Sets("SDS")->Molecules;

#This is the loop over all the frames
for (my \\\$counter = 1; \\\$counter <= \\\$numFrames; ++\\\$counter) {

#Find the current frame and write the current time value into the study table
\\\$doc->Trajectory->CurrentFrame = \\\$counter;
printf "time=%f\\n", \\\$doc->Trajectory->FrameTime;
\\\$calcSheet->Cell(\\\$counter-1, 0) = \\\$doc->Trajectory->FrameTime;

#initializing the number of molecules count
my \\\$numberofmolecules = 0;

my \\\$distance1X = 0; # Define some variables to store the distances in
my \\\$distance1Y = 0;
my \\\$distance1Z = 0;

#handle to all the molecules and atoms in the check distance molecule
my \\\$molecules = \\\$doc->UnitCell->Molecules;
my \\\$moldist = \\\$molecules->Item(\\\$moleculename);

#Loop of all the molecules
foreach my \\\$molecule (@\\\$molecules){

next unless (\\\$molecule->Name ne \\\$moleculename); #exlude the molecule to check distance with

# Get Cartesian coordinates of the center of the molecule
my \\\$molcenterX = \\\$molecule->Center->X;
my \\\$molcenterY = \\\$molecule->Center->Y;
my \\\$molcenterZ = \\\$molecule->Center->Z;

#Loop of the atoms in the molecule to check distance with
foreach my \\\$atom (@{\\\$moldist->Atoms}) {

#Calculate the distance between the molecular center and the atom of the target molecule

\\\$distance1X = distanceCalc (\\\$molcenterX,\\\$atom->X, \\\$lengthX);
\\\$distance1Y = distanceCalc (\\\$molcenterY,\\\$atom->Y, \\\$lengthY);
\\\$distance1Z = distanceCalc (\\\$molcenterZ,\\\$atom->Z, \\\$lengthZ);

# Using the different calculations x y, and z distances, calculate the actual distance

my \\\$dist = sqrt (\\\$distance1X*\\\$distance1X + \\\$distance1Y*\\\$distance1Y + \\\$distance1Z*\\\$distance1Z);

#To count molecule in the distance should be less than cutoff
if(\\\$dist < \\\$cutoff) {
\\\$numberofmolecules++;
last;
}

}

}

#Fill in the study table
printf "%s :%f\\n", "Number of molecules", \\\$numberofmolecules;
\\\$calcSheet->Cell(\\\$counter-1, 1) = \\\$numberofmolecules;

}

Thank you very much!

Best wishes