# Simulating a simple intermolecular repulsion term on a 1D adsorption model

Today I am trying to simulate some molecular adsorption behavior with a really simple repulsion term. What I have so far is a a script that generates a 1D chain, selects random sites (with selection rules dictating that they can only occupy 1 in every three sites) for populating the chain with molecules (which occupy two contiguous sites) and enforces dense-packing constraints at low initial molecular coverages. I am having difficulty implementing the intermolecular repulsion term (implemented using If logic).

At the moment the intermolecular potential is enforced if a RandomReal is greater than some threshold value, with sampling of a new random population point, a test to see if that point is present in my original random sample list and if that point if more favorable than the original point. If it is, the 1D chain is populated at the resampled point.

The 1D chains are then characterized by contiguous regions of 1 or 0 and the number of molecules to populate the chain is then increased by 1.

The code takes about 30 minutes to evaluate on my laptop (depends on the randomreal threshold) and I am sure many of my kludgy workarounds are not the most efficient ways to do things.

I would also like to better populate the 1D chain by enforcing the 1-1-0-1-1-0 densest packing rule only at the highest molecular load (what I mean is that I would like to allow for 1-1-0-0-0-1-1-0-0-1-1-0 type chains, though my random sampling right now only allows for population at every third site.

Do any of you have a better idea how to accomplish this (especially the repulsion term)?

My code follows:

ClearAll["Global*"];(*clear global variables*)n = 1; rep = 1;(*set n = number of molecules, rep = number of replicate tests*)
Do[Do[b =
SparseArray[{1, 1} -> 0, {30000, 1}, 0];(*generate zero matrix*)
a = RandomSample[Range[1, 29997, 3],
n];(*select random sites to add elements*)

Do[If[b[[a[[i]] - 2]] == 1 && b[[a[[i]] + 3]] == 1,
If[RandomReal[{-1, 1}] \[GreaterSlantEqual] 0.75,
b[[a[[i]]]] = 1; b[[a[[i]] + 1]] = 1,
Label[RepulsionLoop]; c = RandomSample[Range[1, 29997, 3], 1];
If[FreeQ[a, c[[1]]] && FreeQ[a, c[[1]] - 3] &&
FreeQ[a, c[[1]] + 3], b[[c[[1]]]] = 1; b[[c[[1]] + 1]] = 1,
Goto[RepulsionLoop]],
b[[a[[i]]]] = 1; b[[a[[i]] + 1]] = 1]], {i, n}];(*insert elements*)
e[i] = b;(*generate iterated variable to record each run*), {i,
rep}];
f = Table[
e[i], {i, 1, rep}];(*generate list from iterated variable*); sites =
Length /@ Select[Split[Flatten[f]], Total[#] > 0 &];
sitesa = BinCounts[Flatten[sites], {0, 30000, 1}];(*count number of molecules - should be only populated at 2*)
voids = Length /@ Select[Split[Flatten[f]], Total[#] < 1 &]; (*count distances between molecules*)
voidsa = BinCounts[Flatten[voids], {0, 30000, 1}];
sitesb[n] = sitesa;
voidsb[n] = voidsa;
, {n, 1, 9999}];
`

edit: fixed some brackets.

-