Osmotic Pressure

globals [
  left-side                                 ;; left side of the membrane
  right-side                                ;; right side of the membrane
  split                                     ;; location of membrane
  pre-split                                 ;; location of membrane in prior step, used to calculate the change in membrane location
  equilibrium                               ;; the difference in number of particles moving each direction across membrane
  tick-delta                                ;; how much we advance the tick counter this time through
  max-tick-delta                            ;; the largest tick-delta is allowed to be
]

breed [solutes solute]
breed [solvents solvent]
breed [membranes membrane]

turtles-own [
  speed
  mass
  energy
  last-collision
  old-pos
  new-pos
  stick-count                              ;; counter for solutes that are stuck to solvents
  linked?                                  ;; flag for solvents stating whether or not they are stuck
]

to setup
  clear-all
  set-default-shape solutes "circle"
  set-default-shape solvents "circle"
  set-default-shape membranes "square"
  set max-tick-delta 0.1073
  create-container
  spawn-particles
  set split 0
  calculate-wall
  reset-ticks
end 

to create-container
  ask patches with [pycor = max-pycor or pycor = min-pycor] [              ;; sets the walls of the container
    set pcolor blue
  ]
  ask patches with [pxcor = max-pxcor or pxcor = min-pxcor] [
    set pcolor blue
  ]
  set left-side patches with [pxcor < 0]                                   ;; defines the left and right sides of the container equally
  set right-side patches with [pxcor > 0]
  set equilibrium 0                                                        ;; sets equilibrium to starting value (the middle of the container)
  set split 0                                                              ;; sets the starting location of the membrane
  ask patches with [pxcor = 0 and pycor != max-pycor and pycor != min-pycor and abs pycor mod 2 != 0] [
    sprout-membranes 1 [                                                   ;; create the membrane
      set color red
      set size 1.3
    ]
  ]
end 

to setup-particles                                                      ;; set variable values for each new particle
  set speed 10
  set mass 1
  set energy (0.5 * mass * speed * speed)
  if breed = solutes [
    set stick-count 0
    set size 1.3
  ]
  set heading random 360
  set linked? false
  set last-collision nobody
end 

;; Osmotic pressure is a colligative property, meaning the number of particles matter more than the identity of the particles.
;; So when spawning particles, we have to take into account whether or not particles dissociate (break up into ions) in the solvent.
;; Covalent compounds, such as sugar, do not break up, while ionic compounds, such as sodium chloride do.  Therefore, when we create 10 solute
;; particles of sugar, 10 solute particles are created.  However, when we create 10 solute particles of sodium chloride, the particles break up into
;; sodium and chloride ions -- so for every one particle of sodium chloride, we get 2 particles of ions (in the model, these still retain the "solute"
;; breed).  Magnesium chloride breaks up into one magnesium ion, and two chloride ions (so 1 magnesium chloride gets you 3 ion particles),
;; and aluminum chloride breaks up into one aluminum ion, and three chloride ions (so 1 aluminum chloride gets you 4 ion particles).
;; 
;; In the code below, this is accomplished by creating solute particles based on the slider value, and then "hatching" new solute particles based on
;; the number of ions the compound forms.  So for sodium chloride, we create the number of particles indicated by the slider, and then each solute
;; spawns an extra solute particle to make a total of two ion particles for each solute particle added.

to spawn-particles
  create-solvents 1000 [                                               ;; creates 1000 solvent molecules
    set color blue + 2
    setup-particles
    move-to one-of patches with [pcolor = black]                       ;; randomly distribute them throughout the world
    while [any? other turtles-here] [
      move-to one-of patches with [pcolor = black]
    ]
  ]

  if Solute_Type = "Sugar" [                       ;; checks to see the identity of the solute based on the chooser
    output-type "C12H22O11"                        ;; write the chemical formula in the output area
    create-solutes-and-ions 1                      ;; since sugar is covalent and it does not break into ions in solvent, one molecule is formed
  ]

  if Solute_Type = "Sodium Chloride" [
    output-type "NaCl"
    create-solutes-and-ions 2                      ;; Sodium Chloride breaks up into Na+ and Cl- ions, so two ions are formed
  ]

  if Solute_Type = "Magnesium Chloride" [
    output-type "MgCl2"
    create-solutes-and-ions 3                      ;; Magnesium Chloride breaks up into Mg2+ and two Cl- ions, so three ions are formed
  ]

  if Solute_Type = "Aluminum Chloride" [
    output-type "AlCl3"
    create-solutes-and-ions 4                      ;; Aluminum Chloride breaks up into Al3+ and three Cl- ions, so four ions are formed
  ]
end 

to create-solutes-and-ions [ions]
  create-solutes solute-left [                              ;; create the number of particles on the left indicated by the slider
      set color white
      setup-particles
      move-to one-of left-side with [pcolor = black]        ;; move to an open space on the left side
      while [any? other turtles-here] [
        move-to one-of left-side with [pcolor = black] ]
      hatch-solutes (ions - 1) [                            ;; since one particle was already created above, 
        set color white                                     ;; we hatch 1 less than the total number of ions the substance has
        setup-particles
        fd 0.2                                              ;; move forward a bit so we can see the particles better
      ]
    ]
    create-solutes solute-right [                           ;; create the number of particles on the right indicated by the slider
      set color white
      setup-particles
      move-to one-of right-side with [pcolor = black]       ;; move to an open space on the right side
      while [any? other turtles-here] [
        move-to one-of right-side with [pcolor = black] ]
      hatch-solutes (ions - 1) [
        set color white
        setup-particles
        fd 0.2
      ]
    ]
end 

to-report particles
  report (turtle-set solutes solvents)
end 

to particle-jump
  ask solutes with [xcor >= pre-split and xcor <= split] [       ;; check for solutes in the way of a membrane jump
    set heading 90                                                                           ;; turn proper direction and jump
    jump (split - pre-split) + 1
  ]
  ask solutes with [xcor <= pre-split and xcor >= split] [
    set heading 270
    jump (pre-split - split) + 1
  ]
end 

to calculate-wall
  set pre-split split                                                  ;; save location of the membrane before it moves
  let nudge ((equilibrium * -1) / 10)                                  ;; calculates the amount to move the membrane
  set split split + nudge                                              ;; set split to be the new location of the membrane
  particle-jump                                                        ;; move solutes in the way of the membrane jump
  ask membranes [
    set xcor split                                                     ;; move membrane turtles according to nudge value
  ]

  set left-side patches with [pxcor < split]                           ;; redefine right and left sides
  set right-side patches with [pxcor > split]
  set equilibrium 0                                                    ;; reset equilibrium so we can keep track of what happens during the next tick
end 

to go
  if count turtles-on right-side = 0 or count turtles-on left-side = 0 [                    ;; stops model if membrane reaches either edge
    user-message "The membrane has burst! Make sure you have some solute on both sides!"
    stop
  ]
  ask particles [ set old-pos xcor ]                            ;; saves starting position of particles
  ask particles [ bounce ]                                      ;; all particles bounce off of walls and solutes bounce off of the membrane
  ask particles with [ linked? = false ] [ move ]               ;; only particles not linked should move -- these are "unstuck" particles
  ask links [ check-for-release ]                               ;; links have a chance of dying -- stuck particles have a chance of getting free
  ask particles [ check-for-collision ]                         ;; particles bounce off of each other
  ask solvents [ check-for-stick ]                              ;; solvent particles can stick to solute molecules

  ask solvents [
    if old-pos < split and xcor >= split [                      ;; if a solvent moves from the left of the membrane to the right of the membrane
      set equilibrium equilibrium + 1                           ;; add one to equilibrium
    ]
    if old-pos > split and xcor <= split [                      ;; if a solvent moves from the right of the membrane to the left of the membrane
      set equilibrium equilibrium - 1                           ;; subtract one from equilibrium
    ]
  ]
  tick-advance tick-delta
  update-plots
  calculate-tick-delta
  display
  calculate-wall                                                ;; recalculate the new location of the wall
end 

to calculate-tick-delta
  ;; tick-delta is calculated in such way that even the fastest
  ;; particle will jump at most 1 patch length when we advance the
  ;; tick counter. As particles jump (speed * tick-delta) each time, making
  ;; tick-delta the inverse of the speed of the fastest particle
  ;; (1/max speed) assures that. Having each particle advance at most
  ;; one patch-length is necessary for it not to "jump over" a wall
  ;; or another particle.
  ifelse any? particles with [speed > 0]
    [ set tick-delta min list (1 / (ceiling max [speed] of particles)) max-tick-delta ]
    [ set tick-delta max-tick-delta ]
end 

to bounce  ;; particle procedure
  ;; get the coordinates of the patch we'll be on if we go forward 1
  let new-patch patch-ahead 1
  let new-px [pxcor] of new-patch
  let new-py [pycor] of new-patch

  ;; if hitting the membrane...
  if any? membranes in-cone 3 180 and breed = solutes
    [ set heading (- heading) ]

  ;; if we're not about to hit a wall, we don't need to do any further checks
  if not shade-of? blue [pcolor] of new-patch
    [ stop ]

  ;; if hitting left or right wall, reflect heading around x axis
  if (abs new-px = max-pxcor) or (abs new-px = min-pxcor)
    [ set heading (- heading) ]
  ;; if hitting top or bottom wall, reflect heading around y axis
  if (abs new-py = max-pycor) or (abs new-py = min-pycor)
    [ set heading (180 - heading)]
end 

to move  ;; particle procedure
  if patch-ahead (speed * tick-delta) != patch-here
    [ set last-collision nobody ]
  jump (speed * tick-delta)
end 

to check-for-release
  ;; there is a 5% chance that a stuck particle will become unstuck
  if random 100 <= 5 [
    ask end1 [
      set stick-count stick-count - 1                                   ;; lower the stick counter of the solute
    ]
    ask end2 [
      set linked? false                                                 ;; set the solvent's flag to "unstuck"
      set color color + 2                                               ;; return the color to normal
    ]
    die
  ]
end 

to check-for-collision  ;; particle procedure
  ;; Here we impose a rule that collisions only take place when there
  ;; are exactly two particles per patch.  We do this because we want them to
  ;; form a uniform wavefront.
  ;;
  ;; Why do we want a uniform wavefront?  Because it is actually more
  ;; realistic.  
  ;;
  ;; Why is it realistic to assume a uniform wavefront?  Because in reality,
  ;; whether a collision takes place would depend on the actual headings
  ;; of the particles, not merely on their proximity.  Since the particles
  ;; in the wavefront have identical speeds and near-identical headings,
  ;; in reality they would not collide.  So even though the two-particles
  ;; rule is not itself realistic, it produces a realistic result.  Also,
  ;; unless the number of particles is extremely large, it is very rare
  ;; for three or more particles to land on the same patch (for example,
  ;; with 400 particles it happens less than 1% of the time).  So imposing
  ;; this additional rule should have only a negligible effect on the
  ;; aggregate behavior of the system.
  ;;
  ;; Why does this rule produce a uniform wavefront?  The particles all
  ;; start out on the same patch, which means that without the only-two
  ;; rule, they would all start colliding with each other immediately,
  ;; resulting in much random variation of speeds and headings.  With
  ;; the only-two rule, they are prevented from colliding with each other
  ;; until they have spread out a lot.  (And in fact, if you observe
  ;; the wavefront closely, you will see that it is not completely smooth,
  ;; because some collisions eventually do start occurring when it thins out while fanning.)

  if (count other turtles-here with [breed != membranes] = 1)
  [
    ;; the following conditions are imposed on collision candidates:
    ;;   1. they must have a lower who number than my own, because collision
    ;;      code is asymmetrical: it must always happen from the point of view
    ;;      of just one particle.
    ;;   2. they must not be the same particle that we last collided with on
    ;;      this patch, so that we have a chance to leave the patch after we've
    ;;      collided with someone.
    let candidate one-of other turtles-here with [breed != membranes and who < [who] of myself and myself != last-collision]
    ;; we also only collide if one of us has non-zero speed. It's useless
    ;; (and incorrect, actually) for two particles with zero speed to collide.
    if (candidate != nobody) and (speed > 0 or [speed] of candidate > 0)
    [
      collide-with candidate
      set last-collision candidate
      ask candidate [ set last-collision myself ]
    ]
  ]
end 

;; implements a collision with another particle.
;;
;; THIS IS THE HEART OF THE PARTICLE SIMULATION, AND YOU ARE STRONGLY ADVISED
;; NOT TO CHANGE IT UNLESS YOU REALLY UNDERSTAND WHAT YOU'RE DOING!
;;
;; The two particles colliding are self and other-particle, and while the
;; collision is performed from the point of view of self, both particles are
;; modified to reflect its effects. This is somewhat complicated, so I'll
;; give a general outline here:
;;   1. Do initial setup, and determine the heading between particle centers
;;      (call it theta).
;;   2. Convert the representation of the velocity of each particle from
;;      speed/heading to a theta-based vector whose first component is the
;;      particle's speed along theta, and whose second component is the speed
;;      perpendicular to theta.
;;   3. Modify the velocity vectors to reflect the effects of the collision.
;;      This involves:
;;        a. computing the velocity of the center of mass of the whole system
;;           along direction theta
;;        b. updating the along-theta components of the two velocity vectors.
;;   4. Convert from the theta-based vector representation of velocity back to
;;      the usual speed/heading representation for each particle.
;;   5. Perform final cleanup and update derived quantities.

to collide-with [ other-particle ] ;; particle procedure
  ;;; PHASE 1: initial setup

  ;; for convenience, grab some quantities from other-particle
  let mass2 [mass] of other-particle
  let speed2 [speed] of other-particle
  let heading2 [heading] of other-particle

  ;; since particles are modeled as zero-size points, theta isn't meaningfully
  ;; defined. we can assign it randomly without affecting the model's outcome.
  let theta (random-float 360)



  ;;; PHASE 2: convert velocities to theta-based vector representation

  ;; now convert my velocity from speed/heading representation to components
  ;; along theta and perpendicular to theta
  let v1t (speed * cos (theta - heading))
  let v1l (speed * sin (theta - heading))

  ;; do the same for other-particle
  let v2t (speed2 * cos (theta - heading2))
  let v2l (speed2 * sin (theta - heading2))



  ;;; PHASE 3: manipulate vectors to implement collision

  ;; compute the velocity of the system's center of mass along theta
  let vcm (((mass * v1t) + (mass2 * v2t)) / (mass + mass2) )

  ;; now compute the new velocity for each particle along direction theta.
  ;; velocity perpendicular to theta is unaffected by a collision along theta,
  ;; so the next two lines actually implement the collision itself, in the
  ;; sense that the effects of the collision are exactly the following changes
  ;; in particle velocity.
  set v1t (2 * vcm - v1t)
  set v2t (2 * vcm - v2t)



  ;;; PHASE 4: convert back to normal speed/heading

  ;; now convert my velocity vector into my new speed and heading
  set speed sqrt ((v1t ^ 2) + (v1l ^ 2))
  set energy (0.5 * mass * speed ^ 2)
  ;; if the magnitude of the velocity vector is 0, atan is undefined. but
  ;; speed will be 0, so heading is irrelevant anyway. therefore, in that
  ;; case we'll just leave it unmodified.
  if v1l != 0 or v1t != 0
    [ set heading (theta - (atan v1l v1t)) ]

  ;; and do the same for other-particle
  ask other-particle [
    set speed sqrt ((v2t ^ 2) + (v2l ^ 2))
    set energy (0.5 * mass * (speed ^ 2))
    if v2l != 0 or v2t != 0
      [ set heading (theta - (atan v2l v2t)) ]
  ]
end 

to check-for-stick ;; particle procedure
  ;; The stick code is similar to the collision code.
  ;; If there are more than one particles on the same patch,
  ;; and the solvent isn't already stuck to another particle,
  ;; the two particles have a chance to stick.
  if (count other turtles-here with [breed != membranes] = 1) and (linked? = false)
  [
    ;; the stick candidate must be a solute and not already be stuck to too many solvents
    let candidate one-of other solutes-here with [stick-count < 5]
    ;; we also only stick if one of us has non-zero speed
    ;; there is a 50% chance to stick
    if (candidate != nobody) and (speed > 0 or [speed] of candidate > 0) and (random 100 <= 50)
    [
      ask candidate [
        create-link-to myself [
          tie
          hide-link
        ]
        set stick-count stick-count + 1                                    ;; keep a running total of how many solvents are stuck
      ]
      set linked? true                                                     ;; flag the solvent as "stuck"
      set color color - 2                                                  ;; tweak the color so we can see it's stuck
    ]
  ]
end