Each walker begins at the origin and at each time-step it takes a step to the ri
ID: 659325 • Letter: E
Question
Each walker begins at the origin and at each time-step it takes a step to the right or left with equal probability, so that pright = 0.5, pleft = 0.5. Use a lattice of spacing ?x=1 and discrete time-steps, ?t=1. Number of steps in each walk: 20 Number of walkers: 10,000. Create a c++ program to find the Mean square displacement vs. time, for time from t = 0, 1, 2, 3... 20. Plot this data and fit a straight line to the curve and find the slope of the curve. Can you define a physics concept that describes what this slope represents? How would you call it? Explain.
Explanation / Answer
void CalcMSD
(
double [NDIR],
stepdata *, stepdata *,
double *,
BOOLEAN, BOOLEAN, BOOLEAN
);
void CalcOptimumDisp
(stepdata *, stepdata *, double[NDIR], BOOLEAN);
void CalcOptMSD
(double [NVOIGHT], stepdata *, stepdata *, double[NDIR], BOOLEAN, BOOLEAN);
void WrapParticle (double[NDIR], double[NDIR], float [NDIR], BOOLEAN[NDIR]);
void CalcMinimumMSD
(
double MSD[NVOIGHT],
stepdata *RefStep,
stepdata *CurStep,
BOOLEAN UseSelectedAtoms
);
void CalcSelectBox(stepdata *, stepdata *, double [2][NDIR]);
/*
************************************************************************
Main Routine
************************************************************************
*/
int main (int argc, char *argv[])
{
CFILE *InputCoordinateFile;
FILE *fin;
stepdata *s = (stepdata *) calloc (1, sizeof(stepdata));
stepdata *r = (stepdata *) calloc (1, sizeof(stepdata));
stepdata *SwapStep = NULL;
BOOLEAN UseStep = FALSE;
BOOLEAN UseRel = FALSE;
BOOLEAN UseIlist = FALSE;
BOOLEAN UseBox = FALSE;
BOOLEAN UseDisp = FALSE;
BOOLEAN PrintDisp = FALSE;
BOOLEAN UseShuffle = FALSE;
BOOLEAN Print3MSD = FALSE;
BOOLEAN Print6MSD = FALSE;
BOOLEAN PrintHead = FALSE;
BOOLEAN IsFirstStep= FALSE;
char *ilistname;
char *refname;
char *curname;
double msd;
double OptDisp[NDIR];
double MSD[NVOIGHT];
double BoxLimit[2][NDIR];
int iatom;
int step;
int c;
int iopt;
/* Make sure Cor defined types are the expected sizes on this system. */
/* On error, this prints a message and stop */
cor_test_defined_types();
/* PARSE COMMAND LINE PARAMETERS */
/* Parse while parameter has - */
iopt = 1;
while (iopt<argc && argv[iopt][0]=='-')
{
c = argv[iopt][1];
switch (c)
{
/* Name of ilist file */
case 'i':
UseIlist = TRUE;
ilistname = argv[iopt+1];
iopt += 1;
break;
/* Select atoms in box */
case 'b':
UseBox = TRUE;
BoxLimit[0][X] = 1e-8*atof(argv[iopt+1]);
BoxLimit[0][Y] = 1e-8*atof(argv[iopt+2]);
BoxLimit[0][Z] = 1e-8*atof(argv[iopt+3]);
BoxLimit[1][X] = 1e-8*atof(argv[iopt+4]);
BoxLimit[1][Y] = 1e-8*atof(argv[iopt+5]);
BoxLimit[1][Z] = 1e-8*atof(argv[iopt+6]);
iopt += 6;
break;
/* Print 3 MSD componenents */
case '3':
Print3MSD = TRUE;
break;
/* Print 6 MSD componenents */
case '6':
Print6MSD = TRUE;
break;
/* Use minimimum displacement calculation */
case 'd':
UseDisp = TRUE;
break;
/* Print column headers */
case 'h':
PrintHead = TRUE;
break;
/* Print displacement calculation */
case 'p':
PrintDisp = TRUE;
break;
/* Minimize MSD by swaping like atoms */
case 'm':
UseShuffle = TRUE;
break;
/* Calculate MSD relative to previous step */
case 'r':
UseRel = TRUE;
break;
/* Turns off warnings */
case 'w':
PrintWarn_m = FALSE;
break;
default :
fprintf (stderr,
"Command line switch %s unknown. ", argv[iopt]);
break;
}
iopt++;
}
/* PRINT USAGE MESSAGE */
if (argc-iopt<1)
{
char **s = Usage_m;
while (*s)
{
printf ("%s ", *s);
s++;
}
return 0;
}
/* Option sanity check */
if (UseShuffle && UseDisp)
{
printf ("ERROR: Cannot use -m with -c or -d. ");
return 1;
}
if (UseStep && UseRel)
{
printf ("ERROR: Cannot use -r with step option. ");
return 1;
}
/* Read file name for reference step */
refname = argv[iopt];
/* read file name for steps to compare to refernce */
/* If asterisk instead of file name, use previous file name */
/* If no more parameters, use previos file name */
if (iopt+1 >= argc || argv[iopt+1][0] == '*')
curname = argv[iopt ];
else
curname = argv[iopt+1];
/* Read reference step number */
UseStep = (argc-iopt>2);
if (UseStep)
step = atoi (argv[iopt+2]);
/* Read Reference Step */
InputCoordinateFile = OpenCoordinateFile (refname, "r");
if (InputCoordinateFile == NULL)
{
printf ("Sorry. Can not open file %s. ", refname);
return 0;
}
/* Find requested step */
if (UseStep)
while
(
ReadCoordinateFile(InputCoordinateFile,r)!=EOF && r->step!=step
);
/* Take first step */
else
ReadCoordinateFile (InputCoordinateFile,r);
CloseCoordinateFile (InputCoordinateFile);
if (UseStep && r->step!=step)
{
printf ("Sorry. ");
printf ("Reference step %i in file %s not found. ", step, refname);
return 1;
}
/* Read selected atoms from ILIST file - store selection in *s */
if (UseIlist)
{
fin = fopen (ilistname, "rt");
if (fin==NULL)
{
fprintf (stderr, "ERROR: Cannot open ilist file %s. ", ilistname);
return 1;
}
readilist (fin, s);
/* Set flag to maintain selection */
s->selkeep = TRUE;
/* Close ILIST file */
fclose (fin);
}
/* Select atom within box - store selection in *s */
if (UseBox)
{
CalcSelectBox (r, s, BoxLimit);
}
/* Copy selection info from s to r if using relative MSD */
if (UseRel && (UseIlist || UseBox) )
{
/* Allocate select array for current step if needed */
REALLOCATE(r->sel,BYTE,s->np)
LOOP (iatom, s->np)
{
r->sel[iatom] = s->sel[iatom];
}
/*
Mark selection array a permanent, not erase when reading new step
*/
r->selkeep = TRUE;
}
/* TEST FOR PRESENCE OF CURRENT FILE */
InputCoordinateFile = OpenCoordinateFile (curname, "r");
if (InputCoordinateFile==NULL)
{
printf ("Sorry. Can not open file %s. ", curname);
return 1;
}
/* Read coordinate file and calculate msd */
IsFirstStep = TRUE;
while ( ReadCoordinateFile(InputCoordinateFile,s)!=EOF )
{
if (IsFirstStep)
{
ZERO_DOUBLE_ARRAY(MSD,NVOIGHT)
ZERO_DOUBLE_ARRAY(OptDisp,NDIR)
}
else
{
CalcMSD
(
MSD,
r,
s,
OptDisp,
UseIlist || UseBox,
UseDisp,
UseShuffle
);
}
/* Column Headers */
if (PrintHead && IsFirstStep)
{
if (Print3MSD)
{
printf ("# step MSD XX MSD YY MSD ZZ");
}
else if (Print6MSD)
{
printf ("# step MSD XX MSD YY MSD ZZ");
printf (" MSD YZ MSD ZX MSD XY");
}
else
{
printf ("# step MSD");
}
if (PrintDisp && UseDisp)
{
printf (" DISP X DISP Y DISP Z");
}
/* Print new line */
printf (" ");
}
/* Column Separators */
if (PrintHead && IsFirstStep)
{
if (Print3MSD)
{
printf ("# ------ ---------- ---------- ----------");
}
else if (Print6MSD)
{
printf ("# ------ ---------- ---------- ----------");
printf (" ---------- ---------- ----------");
}
else
{
printf ("# ------ ----------");
}
if (PrintDisp && UseDisp)
{
printf (" ---------- ---------- ----------");
}
/* Print new line */
printf (" ");
}
/* Print Means Squared Displacement in Requested Format */
if (Print3MSD)
{
printf
(
"%8li %10.4lf %10.4lf %10.4lf",
s->step,
1e16*MSD[XX], 1e16*MSD[YY], 1e16*MSD[ZZ]
);
}
else if (Print6MSD)
{
printf
(
"%8li %10.4lf %10.4lf %10.4lf %10.4lf %10.4lf %10.4lf",
s->step,
1e16*MSD[XX], 1e16*MSD[YY], 1e16*MSD[ZZ],
1e16*MSD[YZ], 1e16*MSD[ZX], 1e16*MSD[XY]
);
}
else
{
printf
(
"%8li %10.4lf",
s->step,
1e16 * ( MSD[XX] + MSD[YY] + MSD[ZZ] )
);
}
/* Print optional optimum displacement */
if (PrintDisp && UseDisp)
{
printf (" %10.4lf %10.4lf %10.4lf",
1e8*OptDisp[X], 1e8*OptDisp[Y], 1e8*OptDisp[Z]);
}
/* Print newline */
printf (" ");
/* If using relative step, save current step as reference */
if (UseRel)
{
SwapStep = r;
r = s;
s = SwapStep;
}
/* Make first step */
IsFirstStep = FALSE;
}
/* while reading steps */
/* CLOSE SHOP */
CloseCoordinateFile (InputCoordinateFile);
return 0;
}
/*
************************************************************************
Local Functions
************************************************************************
*/
/*
Calculate MSD according to one of two algorithms by calling
approriate function
*/
void CalcMSD
(
double MSD[NVOIGHT],
stepdata *RefStep,
stepdata *CurStep,
double OptDisp[NDIR],
BOOLEAN UseSelectedAtoms,
BOOLEAN UseDisplacementMin,
BOOLEAN MinimizeMSD
)
{
int idir;
int ipart;
BOOLEAN OK;
BOOLEAN RefBoundary[NDIR];
BOOLEAN CurBoundary[NDIR];
double Ref [NDIR];
double Cur [NDIR];
/* Test that both steps have same number of particles */
if (RefStep->np != CurStep->np)
{
printf ("ERROR: Number of atoms (%li) for reference step (%i) is ",
RefStep->step, RefStep->np);
printf (" different from number of atoms (%li) for current step (%i). ",
CurStep->step, CurStep->np);
exit (1);
}
/* Map repeating flag to TRUE or FALSE */
LOOP (idir, NDIR)
{
RefBoundary[idir] = RefStep->boundrep[idir];
CurBoundary[idir] = CurStep->boundrep[idir];
}
/* Test that both steps have repeating boundary conditions */
OK = TRUE;
LOOP (idir, NDIR)
if (RefBoundary[idir] != CurBoundary[idir])
OK = FALSE;
if (!OK && PrintWarn_m)
{
printf ("# WARNING: Boundary conditions (%1i,%1i,%1i) ",
RefBoundary[X], RefBoundary[Y], RefBoundary[Z]);
printf ("for Reference step (%li) ",
RefStep->step);
printf ("# are difference from boundary conditions (%1i,%1i,%1i) ",
CurBoundary[X], CurBoundary[Y], CurBoundary[Z]);
printf ("for Current step (%li) ",
CurStep->step);
}
/* Test that both steps have same box sizes */
OK = TRUE;
LOOP (idir, NDIR)
if (RefBoundary[idir] && CurBoundary[idir])
if (RefStep->box[idir] != CurStep->box[idir])
OK = FALSE;
if (!OK && PrintWarn_m)
{
/* Print out repeating box warnging */
printf ("# WARNING: Repeating box sizes for Reference step ");
printf ("(%li) and Current step (%li) ",
RefStep->step, CurStep->step);
printf ("# are not equal ");
printf ("# The repeating boundary conditions are.. ");
printf ("# DIR RefStep CurStep ");
printf ("# --- ---------- ---------- ");
/* Print out table of repeating boundary box sizes */
LOOP (idir, NDIR)
{
if (RefBoundary[idir] && CurBoundary[idir])
printf
(" %c: %10.4f %10.4f ",
'X'+idir,
RefStep->box[idir],
CurStep->box[idir]
);
}
}
/*
***************
MSD Calculation
***************
*/
if (MinimizeMSD)
CalcMinimumMSD
(MSD, RefStep, CurStep, UseSelectedAtoms);
else
CalcOptMSD
(MSD, RefStep, CurStep, OptDisp, UseSelectedAtoms, UseDisplacementMin);
}
/*
Calculate MSD while re-arranging atoms to find minimum.
For example, atom 5 in the reference step might have become
atom 2 in the current step, if so these two atoms should be
compared to get displacement. This function *tries* to find
the correct atom matches between two steps by the following method.
This function traverses the list of N reference particles N times.
The first time, it matches ref particle 1 with its closest cur particle.
The second time, it matches ref particle 2 tihe *its* closest cur
particle, etc.
*/
void CalcMinimumMSD
(
double MSD[NVOIGHT],
stepdata *RefStep,
stepdata *CurStep,
BOOLEAN UseSelectedAtoms
)
{
int MsdCnt;
int ipart;
int jpart;
int jmatch;
int ivoight;
int idir;
int itype;
int RefType;
int BestPart;
int BestMatch;
double Ref[NDIR];
double Cur[NDIR];
double Disp[NDIR];
double BestDisp[NDIR];
int RefTypeCount[256];
int CurTypeCount[256];
int *Match;
float DispSqr;
float MinDispSqr;
/* Test for same number of each atom type */
LOOP (itype, 256)
CurTypeCount[itype] = RefTypeCount[itype] = 0;
LOOP (ipart, RefStep->np)
{
CurTypeCount[CurStep->t[ipart]]++;
RefTypeCount[RefStep->t[ipart]]++;
}
LOOP (itype, 256)
{
if (CurTypeCount[itype]!=RefTypeCount[itype])
{
printf ("ERROR: ");
printf ("Steps %li and %li have different numbers of each type. ",
RefStep->step, CurStep->step);
exit (1);
}
}
/* Allocate match */
Match = (int *) calloc (RefStep->np, sizeof(int) );
if (Match==NULL)
{
printf ("ERROR: ");
printf ("Cannot allocate index array for re-arranging like atoms. ");
exit (1);
}
/* Initialize MSD */
ZERO_DOUBLE_ARRAY(MSD,NVOIGHT)
/* Initialize match (keeps track of which atoms match up) */
LOOP (jmatch, RefStep->np)
Match[jmatch] = jmatch;
/* Calculate MSD for selected particles and remove Center of mass */
MsdCnt = 0 ;
LOOP (ipart, RefStep->np)
if (!UseSelectedAtoms || CurStep->sel[ipart])
{
/* Save particle coordinates */
RefType = RefStep->t[ipart];
LOOP (idir, NDIR)
Ref[idir] = RefStep->c[idir][ipart];
/* Initialize minimum displacment */
MinDispSqr =
SQR (CurStep->box[X]) +
SQR (CurStep->box[Y]) +
SQR (CurStep->box[Z]);
/* Compare current particle to all others of same type */
for (jmatch=ipart; jmatch<RefStep->np; jmatch++)
{
/* Get particle current particle index */
jpart = Match[jmatch];
/* Test if current and reference particle are same type */
if (RefType==CurStep->t[jpart])
{
/* Use copy of current coorindate */
LOOP (idir, NDIR)
{
Cur[idir] = CurStep->c[idir][jpart];
}
/* Wrap current particle close to reference particle */
WrapParticle (Ref, Cur, RefStep->box, RefStep->boundrep);
/* Calculate current displacement squared */
Disp[X] = Cur[X] - Ref[X];
Disp[Y] = Cur[Y] - Ref[Y];
Disp[Z] = Cur[Z] - Ref[Z];
DispSqr = SQR(Disp[X]) + SQR(Disp[Y]) + SQR(Disp[Z]);
/* Compare to minimum */
if (DispSqr < MinDispSqr)
{
BestMatch = jmatch;
MinDispSqr= DispSqr;
BestDisp[X] = Disp[X];
BestDisp[Y] = Disp[Y];
BestDisp[Z] = Disp[Z];
}
}
}
/* Swap particles */
BestPart = Match[BestMatch];
Match[BestMatch] = Match[ipart];
Match[ipart ] = BestPart;
/* Sum minimum squared displacement into MSD */
MSD[XX] += BestDisp[X] * BestDisp[X];
MSD[YY] += BestDisp[Y] * BestDisp[Y];
MSD[ZZ] += BestDisp[Z] * BestDisp[Z];
MSD[YZ] += BestDisp[Y] * BestDisp[Z];
MSD[ZX] += BestDisp[Z] * BestDisp[X];
MSD[XY] += BestDisp[X] * BestDisp[Y];
/* Sum number of contributions */
MsdCnt++;
}
/* Normalize */
if (MsdCnt>1)
LOOP (ivoight, NVOIGHT)
MSD[ivoight] /= MsdCnt;
}
/*
Calc MSD while optionally compensating for uniform
displacement of selected particles. Calculate *optimum*
displacement using function CalcOptimumDisp()
*/
void CalcOptMSD
(
double MSD[NVOIGHT],
stepdata *RefStep,
stepdata *CurStep,
double OptDisp[NDIR],
BOOLEAN UseSelectedAtoms,
BOOLEAN UseDisplacementMin
)
{
int MsdCnt;
int ipart;
int idir;
int ivoight;
double Ref[NDIR];
double Cur[NDIR];
double Disp[NDIR];
/*
**************************
Optimum Displacement
**************************
*/
OptDisp[X] = OptDisp[Y] = OptDisp[Z] = 0.0;
if (UseDisplacementMin)
{
CalcOptimumDisp (RefStep, CurStep, OptDisp, UseSelectedAtoms);
}
/* Calculate MSD for selected particles and remove Center of mass */
ZERO_DOUBLE_ARRAY(MSD,NVOIGHT)
MsdCnt = 0;
LOOP (ipart, RefStep->np)
if (!UseSelectedAtoms || CurStep->sel[ipart])
{
/* Save particle coordinates */
LOOP (idir, NDIR)
{
Ref[idir] = RefStep->c[idir][ipart];
Cur[idir] = CurStep->c[idir][ipart] - OptDisp[idir];
}
/* Wrap current particle close to reference particle */
WrapParticle (Ref, Cur, RefStep->box, RefStep->boundrep);
/* Calculate displacement between reference and current step */
LOOP (idir, NDIR)
{
Disp[idir] = Cur[idir] - Ref[idir];
}
/* Sum Squared Displacement Tensor */
MSD[XX] += Disp[X] * Disp[X];
MSD[YY] += Disp[Y] * Disp[Y];
MSD[ZZ] += Disp[Z] * Disp[Z];
MSD[YZ] += Disp[Y] * Disp[Z];
MSD[ZX] += Disp[Z] * Disp[X];
MSD[XY] += Disp[X] * Disp[Y];
MsdCnt++;
}
/* Normalize MSD */
if (MsdCnt>1)
LOOP (ivoight, NVOIGHT)
MSD[ivoight] /= MsdCnt;
}
/*
Calculate the uniform displacement of selected particles
which minimizes the residual mean square displacement
For Repeating Boundaries algorithm is approximate
but should work well for 'typical' cases where there is
a near uniform translation.
For non-repeating boundaries algorithm is exact.
*/
void CalcOptimumDisp
(
stepdata *RefStep,
stepdata *CurStep,
double OptDisp[NDIR],
BOOLEAN UseSelectedAtoms
)
{
int iatom;
int idir;
int NumSel;
double xsum;
double ysum;
double Angle;
double AvgAngle;
double Disp;
double AvgDisp;
double MapFactor;
/* Initialize sum of displacements to zero */
LOOP (idir, NDIR)
OptDisp[idir] = 0.0;
/* Find number of selected particles */
NumSel = 0;
LOOP (iatom, RefStep->np)
{
if (!UseSelectedAtoms || CurStep->sel[iatom])
NumSel++;
}
/*
Find optimum displacement for those directions
WITHOUT repeating boundaries.
*/
LOOP (idir, NDIR)
if (!RefStep->boundrep[idir])
{
OptDisp[idir] = 0.0;
LOOP (iatom, RefStep->np)
if (!UseSelectedAtoms || CurStep->sel[iatom])
{
OptDisp[idir] +=
CurStep->c[idir][iatom] - RefStep->c[idir][iatom];
}
if (NumSel>1)
OptDisp[idir] /= NumSel;
}
/*
Find optimum displacement for those directions
WITH repeating boundaries
*/
LOOP (idir, NDIR)
if (RefStep->boundrep[idir])
{
/*
Factor to map displacements from 0..Box[idir] to 0..2*pi
*/
MapFactor = (2.0*M_PI) / CurStep->box[idir];
/*
Mapping repeating direction onto circle perimeter,
find average angle
*/
xsum = ysum = 0.0;
LOOP (iatom, RefStep->np)
if (!UseSelectedAtoms || CurStep->sel[iatom])
{
Angle =
MapFactor *
(CurStep->c[idir][iatom] - RefStep->c[idir][iatom]);
xsum += cos(Angle);
ysum += sin(Angle);
}
if (NumSel==0.0)
AvgAngle = 0.0;
else if (xsum==0.0 && ysum==0.0)
AvgAngle = 0.0;
else
AvgAngle = atan2 (ysum, xsum);
/* Find approximate average displacement from average angle */
AvgDisp = AvgAngle / MapFactor;
#if 0
/* Calcualate average displacements with bias */
OptDisp[idir] = 0.0;
LOOP (iatom, RefStep->np)
if (!UseSelectedAtoms || CurStep->sel[iatom])
{
Disp =
CurStep->c[idir][iatom] - RefStep->c[idir][iatom]
- AvgDisp;
if (Disp<-RefStep->box[idir]) Disp += RefStep->box[idir];
else if (Disp>+RefStep->box[idir]) Disp -= RefStep->box[idir];
OptDisp[idir] += Disp;
}
if (NumSel>1)
OptDisp[idir] /= NumSel;
OptDisp[idir] += AvgDisp;
#endif
OptDisp[idir] = AvgDisp;
}
}
/*
Selects atoms within the box
*/
void CalcSelectBox
(
stepdata *ref,
stepdata *cur,
double BoxLimit[2][NDIR]
)
{
int iatom;
int idir;
BOOLEAN IsInBox;
/* stepdata structures not present, return */
if (ref==NULL || cur==NULL)
return;
/* Need to allocate sel[] array */
REALLOCATE(cur->sel,BYTE,ref->np)
if (cur->np<ref->np)
AllocateParticles(cur,ref->np);
/* Test each atom to see if withing BoxLimits[] */
LOOP (iatom, ref->np)
{
IsInBox = TRUE;
for (idir=0; idir<NDIR && IsInBox; idir++)
{
IsInBox =
(
IsInBox &&
(BoxLimit[0][idir] <= ref->c[idir][iatom]) &&
(BoxLimit[1][idir] > ref->c[idir][iatom])
);
}
/* Select this atom if in box (atom may already be selected
* previous to this function
*
*/
if (IsInBox)
cur->sel[iatom] = TRUE;
}
/* Set select flags as permanant */
cur->selkeep = TRUE;
}
/* Wrap particles in Cur by box size so that it is closest to Ref */
void WrapParticle
(
double Ref[NDIR],
double Cur[NDIR],
float Box[NDIR],
BOOLEAN Boundary[NDIR]
)
{
int idir;
double Del;
double HalfBox;
LOOP (idir, NDIR)
{
if (Boundary[idir])
{
Del = Ref[idir] - Cur[idir];
HalfBox = 0.5 * Box[idir];
if (Del < -HalfBox) Cur[idir] -= Box[idir];
else if (Del > HalfBox) Cur[idir] += Box[idir];
}
}
}
Related Questions
drjack9650@gmail.com
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.