ISO: Modeling top views for structures of atoms in surface planes for liquids

Share icon

I would like to create a pattern of circles on a grid to model the structure of atoms or molecules looking at a top view in a surface layer. I'd like to do this for 2-D random amorphous/glassy structures (i.e. liquid surfaces).

I just want the top view, not a 3-D view.

My idea is to have a control panel with an input for relative separation distance factor <R> = <r>/r_o where r_o is the zero-point equilibrium separation and <r> is to be the mean separation of the circles averaged over all circles. Calculations would place N circles with radius r_o randomly on a grid to meet a constraint that the separation factor averages to <R>. Perhaps this would also need a sigma distribution factor input to drive the random placements. Perhaps this should also have an input for N (I am thinking to keep this fixed initially at some manageable number for computations).

Does anyone have a suggestion for a starting point that I can use to build this type of visualization?

That's an interesting problem. It sounds easy, but I don't think it is. You may have to start with a close-packed hexagonal structure and then do a Monte Carlo simulation where energy is added to the system for x timesteps.

March 18, 2022 at 05:46 am - Permalink

Maybe you could start with drawing one circle in the center of your image. You then assign probabilities to every pixel of your image for being the position of the next circle. The probability would be zero if the new circle would overlap with an existing circle, but increasingly positive up until the equilibrium distance you want and then decreases after that.

Based on all those probabilities you select the next pixel to be the center of a circle. You now recalculate all probabilities and based on those add a third circle, etc... You would then stepwise add circles to your image until your image is full. The borders could end up being a little funny, depending on when you stop adding circles, but you could fix that by cropping the borders away before displaying the image.

March 18, 2022 at 06:11 am - Permalink

I am not wanting to step into doing a time evolution problem. I thought in essence to make an NxN grid, randomly choose locations on the grid while disallowing for overlap of circles, and set circles on the grid at the points chosen. If I break this down, I might ask in sequence

* What is the best way to make an N x N image/graph with a circle at a random location on the image/graph? Should I make an NxN image matrix ... but then how to I add a circle? Should I make a y(x) array to display "points as circles" ... but then should the array be M = N^2 in size?

* What is the best way to sequence random choices in (x, y) locations to remove points in regions that have already been selected? Should I just overfill the matrix / array and then post-process? Should I be able to pre-eliminate regions after each point is selected so that the next selection ignores the bounds?

This is my first step to model structure at liquid surfaces ... model a perfectly flat plane with a disordered arrangement of molecules (circles). A next might be to increase temperature ... move the molecules apart. Another next might be to add forces below the plane that pull the molecules inward ... make the circles smaller. Another next step might be to repeat the plane to show one layer below.

I was driven to this as I make a summary report about the structure on liquid surfaces not being any different laterally in the plane than what we find in planes sliced in the bulk liquid. I wanted to create this image without have to do the work to manually move each circle.

It may be my eventual next, next, next ... step (this summer to fall) to apply this practice to show ordered patterns, e.g. the nine (3 x 3) top views of the three low index plane families in three cubic metal single crystal. But slicing planes already exist in software packages, and I can draw well-structured arrangements in drawing apps.
liquid surface planes (704.12 KB)

March 18, 2022 at 06:13 am - Permalink

I think time evolution is the only way to get a scientifically accurate representation, but the simpler approaches might work well enough for a crude visual representation of an amorphous surface.

About liquid surfaces not being any different laterally than a sliced bulk plane. I think that is an approximation. I would expect the outer-most surface plane of a liquid to have a slighter denser packaging than the bulk, to compensate for the missing bonds at the surface. For a liquid consisting of molecules you may also have a preferred orientation towards the air/vacuum, and the surface itself may even have a different melting point than the bulk.

March 18, 2022 at 09:12 am - Permalink

I imagine that we cannot create a spontaneous (self-directed) change in lateral surface area number density by invoking bond forces that are asymmetric only in normal directions to the surface plane. The lack of bond forces above the liquid surface plane can however be considered to cause the top surface layer to move (relax) inward toward the bulk liquid. Otherwise, going from an arbitrary initial state arrangement to a denser area number density in a lateral surface plane is also entropically unfavored (think going to fewer partition states / choices or going toward higher order). I cannot imagine therefore that laterally denser layers will spontaneously form unless external work is done on them (e.g. to compress them). Differences in melting points are solid surface considerations.

In any case, I'm just after a crude approximation as a possible hobby. I leave the time evolution, sophisticated modeling to folks who know how to do it far better.

https://pubs.acs.org/doi/10.1021/la403421b

March 18, 2022 at 09:39 am - Permalink

Generally speaking, when you remove one bond from an atom the remaining bonds are strengthened. For surfaces of crystalline solids this means the bond lengths between the first and second layer gets shorter, as you describe. The bonds laterally within the surface plane do not shorten because that would make the surface structure incommensurate with the bulk.

However, the surface of an amorphous liquid does not have to be commensurate with the bulk and hence there is no reason why the bonds in the surface would not contract and make the density in the surface couldn't be slighter higher than in bulk.

However, I freely admit that I know a lot more about crystalline surfaces than amorphous surfaces.

March 18, 2022 at 10:26 am - Permalink

Chemistry. Especially considering secondary bonds (that are non-directional and are not "removable" as such) instead of primary bonds (that can be "removed"). Then put this all in a liquid and try to argue whether the surface is more or less dense laterally as well as perpendicularly to the surface. Add a dynamic event such as suspending a water drop on a flat water surface. Try to teach this all to a class of graduate students with backgrounds ranging from chemistry to mechanical engineering. Do so without embarrassing yourself at your own ignorance. Fun stuff. I'd certainly prefer to be teaching about crystalline solid surfaces that are held in UHV. :-)

March 18, 2022 at 10:35 am - Permalink

Or the stuff we saw in TEM samples of intensely sheared powdered rocks (think ceramics) that looked to the TEM to be amorphous.

That was 26 years ago. I don't know what the state of the art is now!

March 18, 2022 at 11:41 am - Permalink

This is my solution. It's a big pile of pseudo science that would never pass peer review, but it might be good enough for what you had in mind.

I build up the amorphous structure one atom at a time. For each atom I calculate a Lennard-Jones-like potential around the atom, and I assume the combined potential when multiple atoms exist is just the sum of the individual potentials.

Based on the potential energy landscape I then add additional atoms one at a time using a Boltzmann-like probability distribution.

 

Function PseudoScience()
//  Creates a pseudo-scientific representation of 2D amorphous surface
 
    //  The size of the surface and the number of atoms to add to it
    Variable XSize=500, YSize=500
    Variable NumberOfAtoms=400
    Variable AtomicRadius=8         //  This is for the visual representation only
    Variable LennardJonesRadius=20      //  This is for calculating the potential
    Variable LennardJonesWell=20e3      //  20 kJ/mol potential-well depth
    Variable GasConstant=8.31
    Variable Temperature=500
    
    //  Creates four waves describing the area immediately around a single isolated atom
    //  The atom is considered to have no effect beyond this range ( > 2*LennardJonesRadius away from the center of the atom )
    Make/O/N=(4*LennardJonesRadius+1, 4*LennardJonesRadius+1) RadialDistanceMap, SinglePotentialEnergyMap, SingleProbabilityMap, SingleVisualRepresentation
 
    //  A radial distance map used below to create the potential and visual representation of a single isolated atom
    MultiThread RadialDistanceMap=Sqrt((p-2*LennardJonesRadius)^2+(q-2*LennardJonesRadius)^2)   
    
    //  The 2D potential energy surface for a single isolated atom, assuming a Lennard-Jones potential straight out of Wikipedia
    MultiThread SinglePotentialEnergyMap[][]=4*LennardJonesWell*((LennardJonesRadius/RadialDistanceMap[p][q])^12-(LennardJonesRadius/RadialDistanceMap[p][q])^6)
    SinglePotentialEnergyMap[2*LennardJonesRadius][2*LennardJonesRadius]=inf        //  Changes the center point from NaN to +inf
    
    //  Calculates a Boltzmann distribution around the atom, based on the Lennard-Jones potential
    FastOP SingleProbabilityMap=(-1/(GasConstant*Temperature))*SinglePotentialEnergyMap
    MultiThread SingleProbabilityMap=exp(SingleProbabilityMap)
 
    //  Creates a visual representation of an isolated atom as a circle
    MultiThread SingleVisualRepresentation[][]=((RadialDistanceMap[p][q]>AtomicRadius) ? 0 : 1)
 
    //  The combined 2D potential energy surface for every atom added, the combined Boltzmann probability map and the combined visual representation
    Make/O/N=(XSize, YSize) CombinedPotentialEnergyMap, CombinedProbabilityMap, CombinedVisualRepresentation
    FastOP CombinedPotentialEnergyMap=0
    FastOP CombinedVisualRepresentation=0
    FastOP CombinedProbabilityMap=1
    
    //  The wave containing the X and Y positions of every atom added
    Make/O/N=(NumberOfAtoms, 2) AtomicPositions
    FastOP AtomicPositions=(NaN)
    
    //  Needed in the loop to select a random position
    Make/O/N=(XSize*YSize+1) IntegralWave
    
    //  Adds one atoms at a time
    Variable i=0, SelectedPosition=0, XOffset=0, YOffset=0, X1=0, X2=0, Y1=0, Y2=0
    for (i=0; i<NumberOfAtoms; i+=1)
    
        //  Adds an atom at a random position. For the first step in the loop all positions have equal probabilities
        //  The probability map is first converted to a one-dimensional wave to make it compatible with the Integrate command
        Redimension/N=(NumPnts(CombinedProbabilityMap)) CombinedProbabilityMap
        SelectedPosition=SelectRandomPosition(CombinedProbabilityMap, IntegralWave)
        Redimension/N=(XSize, YSize) CombinedProbabilityMap
        
        //  Converts the 1D index returned by SelectRandomPosition to X and Y positions in the image
        AtomicPositions[i][0]=mod(SelectedPosition, XSize)      //  X position
        AtomicPositions[i][1]=Floor(SelectedPosition/XSize)     //  Y Position
        
        //  Calculates the X and Y offsets required to get SingleProbabilityMap into the right positions
        XOffset=AtomicPositions[i][0]-(DimSize(SingleProbabilityMap, 0)-1)/2
        YOffset=AtomicPositions[i][1]-(DimSize(SingleProbabilityMap, 1)-1)/2
 
        //  Calculates the range of data points to iterate over to prevent going out of range
        X1=max(0, XOffset)
        X2=min(DimSize(CombinedProbabilityMap, 0)-1, AtomicPositions[i][0]+(DimSize(SingleProbabilityMap, 0)-1)/2)
        
        Y1=max(0, YOffset)
        Y2=min(DimSize(CombinedProbabilityMap, 1)-1, AtomicPositions[i][1]+(DimSize(SingleProbabilityMap, 1)-1)/2)
        
        //  Calculates the new Boltzmann probability distribution, taking advantage of exp(A+B) = exp(A) * exp(B)
        //  Meaning the if the potential energy landscape is a sum of two contributions, the Boltzmann distribution becomes the product of the individual distributions
        //  This means I don't have to calculate exp( ) inside the loop
        MultiThread CombinedProbabilityMap[X1, X2][Y1, Y2]=CombinedProbabilityMap[p][q]*SingleProbabilityMap[p-XOffset][q-YOffset]
 
        //  Adds the new atom to the visual representation
        MultiThread CombinedVisualRepresentation[X1, X2][Y1, Y2]=CombinedVisualRepresentation[p][q]+SingleVisualRepresentation[p-XOffset][q-YOffset]
 
    endfor
 
    //  Displays the visual representation
    DoWindow /K VisualRepresentationWin
    Display /W=(200,200,800,700) /K=1 /N=VisualRepresentationWin
    AppendImage /W=VisualRepresentationWin CombinedVisualRepresentation
end
 
 
 
Function SelectRandomPosition(ProbabilityWave, IntegralWave)
//  Selects a random point in ProbabilityWave weighted by the values in the wave
//  There must be a built-in function for this, but I couldn't find it 
Wave ProbabilityWave, IntegralWave
 
    //  Integrates all values in the wave
    Integrate /METH=2/P ProbabilityWave /D=IntegralWave
 
    //  Selects a random number between 0 and the maximum value of the integrated wave
    Variable RandomNumber=abs(enoise(IntegralWave[NumPnts(IntegralWave)-1]))
    
    //  Finds the crossing of the random number
    FindLevel/Q/P IntegralWave, RandomNumber
 
    //  Returns the corresponding data point
    Return Floor(V_LevelX)
end

 

March 21, 2022 at 05:02 am - Permalink

It takes about one second to run on my office computer.

You could pretty easily convert this into three dimensions, but the calculation may get very, very slow very quickly

March 21, 2022 at 05:07 am - Permalink

Wow! I am extremely honored that you took the time to put this together. Thanks.

March 21, 2022 at 05:29 am - Permalink

If you create the graph before the for loop and add a DoUpdate inside the loop, you get a small movie of the atoms being added

Function PseudoScience()
//  Creates a pseudo-scientific representation of 2D amorphous surface
//  Created in response to a thread by jjweimer on the Igor forums "ISO: Modeling top views for structures of atoms in surface planes for liquids "
 
    //  The size of the surface and the number of atoms to add to it
    Variable XSize=500, YSize=500
    Variable NumberOfAtoms=400
    Variable AtomicRadius=8         //  This is for the visual representation only
    Variable LennardJonesRadius=20      //  This is for calculating the potential
    Variable LennardJonesWell=20e3      //  20 kJ/mol potential-well depth
    Variable GasConstant=8.31
    Variable Temperature=500
    
    //  Creates four waves describing the area immediately around a single isolated atom
    //  The atom is considered to have no effect beyond this range ( > 2*LennardJonesRadius away from the center of the atom )
    Make/O/N=(4*LennardJonesRadius+1, 4*LennardJonesRadius+1) RadialDistanceMap, SinglePotentialEnergyMap, SingleProbabilityMap, SingleVisualRepresentation
 
    //  A radial distance map used below to create the potential and visual representation of a single isolated atom
    MultiThread RadialDistanceMap=Sqrt((p-2*LennardJonesRadius)^2+(q-2*LennardJonesRadius)^2)   
    
    //  The 2D potential energy surface for a single isolated atom, assuming a Lennard-Jones potential straight out of Wikipedia
    MultiThread SinglePotentialEnergyMap[][]=4*LennardJonesWell*((LennardJonesRadius/RadialDistanceMap[p][q])^12-(LennardJonesRadius/RadialDistanceMap[p][q])^6)
    SinglePotentialEnergyMap[2*LennardJonesRadius][2*LennardJonesRadius]=inf        //  Changes the center point from NaN to +inf
    
    //  Calculates a Boltzmann distribution around the atom, based on the Lennard-Jones potential
    FastOP SingleProbabilityMap=(-1/(GasConstant*Temperature))*SinglePotentialEnergyMap
    MultiThread SingleProbabilityMap=exp(SingleProbabilityMap)
 
    //  Creates a visual representation of an isolated atom as a circle
    MultiThread SingleVisualRepresentation[][]=((RadialDistanceMap[p][q]>AtomicRadius) ? 0 : 1)
 
    //  The combined 2D potential energy surface for every atom added, the combined Boltzmann probability map and the combined visual representation
    Make/O/N=(XSize, YSize) CombinedPotentialEnergyMap, CombinedProbabilityMap, CombinedVisualRepresentation
    FastOP CombinedPotentialEnergyMap=0
    FastOP CombinedVisualRepresentation=0
    FastOP CombinedProbabilityMap=1
    
    //  Displays the visual representation
    //  Creating the graph before the for loop allows me to show the atoms being added as a small movie using DoUpdate inside the for loop
    DoWindow /K VisualRepresentationWin
    Display /W=(200,200,800,700) /K=1 /N=VisualRepresentationWin
    AppendImage /W=VisualRepresentationWin CombinedVisualRepresentation
    
    //  The wave containing the X and Y positions of every atom added
    Make/O/N=(NumberOfAtoms, 2) AtomicPositions
    FastOP AtomicPositions=(NaN)
    
    //  Needed in the loop to select a random position
    Make/O/N=(XSize*YSize+1) IntegralWave
    
    //  Adds one atoms at a time
    Variable i=0, SelectedPosition=0, XOffset=0, YOffset=0, X1=0, X2=0, Y1=0, Y2=0
    for (i=0; i<NumberOfAtoms; i+=1)
    
        //  Adds an atom at a random position. For the first step in the loop all positions have equal probabilities
        //  The probability map is first converted to a one-dimensional wave to make it compatible with the Integrate command
        Redimension/N=(NumPnts(CombinedProbabilityMap)) CombinedProbabilityMap
        SelectedPosition=SelectRandomPosition(CombinedProbabilityMap, IntegralWave)
        Redimension/N=(XSize, YSize) CombinedProbabilityMap
        
        //  Converts the 1D index returned by SelectRandomPosition to X and Y positions in the image
        AtomicPositions[i][0]=mod(SelectedPosition, XSize)      //  X position
        AtomicPositions[i][1]=Floor(SelectedPosition/XSize)     //  Y Position
        
        //  Calculates the X and Y offsets required to get SingleProbabilityMap into the right positions
        XOffset=AtomicPositions[i][0]-(DimSize(SingleProbabilityMap, 0)-1)/2
        YOffset=AtomicPositions[i][1]-(DimSize(SingleProbabilityMap, 1)-1)/2
 
        //  Calculates the range of data points to iterate over to prevent going out of range
        X1=max(0, XOffset)
        X2=min(DimSize(CombinedProbabilityMap, 0)-1, AtomicPositions[i][0]+(DimSize(SingleProbabilityMap, 0)-1)/2)
        
        Y1=max(0, YOffset)
        Y2=min(DimSize(CombinedProbabilityMap, 1)-1, AtomicPositions[i][1]+(DimSize(SingleProbabilityMap, 1)-1)/2)
        
        //  Calculates the new Boltzmann probability distribution, taking advantage of exp(A+B) = exp(A) * exp(B)
        //  Meaning the if the potential energy landscape is a sum of two contributions, the Boltzmann distribution becomes the product of the individual distributions
        //  This means I don't have to calculate exp( ) inside the loop
        MultiThread CombinedProbabilityMap[X1, X2][Y1, Y2]=CombinedProbabilityMap[p][q]*SingleProbabilityMap[p-XOffset][q-YOffset]
 
        //  Adds the new atom to the visual representation
        MultiThread CombinedVisualRepresentation[X1, X2][Y1, Y2]=CombinedVisualRepresentation[p][q]+SingleVisualRepresentation[p-XOffset][q-YOffset]
        
        //  Updates the graph displaying to visual representation
        //  Remove this to speed up the calculations
        DoUpdate /W=VisualRepresentationWin
    endfor
end
 
 
 
Function SelectRandomPosition(ProbabilityWave, IntegralWave)
//  Selects a random point in ProbabilityWave weighted by the values in the wave
//  There must be a built-in function for this, but I couldn't find it 
Wave ProbabilityWave, IntegralWave
 
    //  Integrates all values in the wave
    Integrate /METH=2/P ProbabilityWave /D=IntegralWave
 
    //  Selects a random number between 0 and the maximum value of the integrated wave
    Variable RandomNumber=abs(enoise(IntegralWave[NumPnts(IntegralWave)-1]))
    
    //  Finds the crossing of the random number
    FindLevel/Q/P IntegralWave, RandomNumber
 
    //  Returns the corresponding data point
    Return Floor(V_LevelX)
end

 

March 22, 2022 at 01:55 am - Permalink