/****************************************************************************
 * Test 01: Simple volume macromolecules.
 *
 *    This case is simple enough that it may be solved analytically.  When I
 *    was originally working on the macromolecules, I compared against the
 *    analytically derived solution.
 *
 *    An analyically derived solution is:
 *        d_0 ~31.1
 *        d_1 ~30.0
 *        d_2 ~28.9
 *
 *    where d_0 is the number of dimer partners in the doubly-unbound state,
 *    d_1 is the number in the singly-bound state, and d_2 is the number in the
 *    doubly-bound state.  The counts below relate in the following way:
 *
 *      00 -> 2*d_0     (number of subunits in doubly unbound dimer pairs)
 *      01 -> d_1       (number of bound subunits with unbound partners)
 *      10 -> d_1       (number of unbound subunits with bound partners)
 *      11 -> 2*d_2     (number of subunits in doubly bound dimer pairs)
 *      us -> 2*d_0+d_1 (number of unbound subunits)
 *      bs -> 2*d_2+d_1 (number of unbound subunits)
 *      cam -> 800 - 2*d_2 - d_1    (number of free Cam)
 *
 *    01 and 10 should be in agreement at all times.
 *
 * Author: Jed Wing <jed@salk.edu>
 * Date:   2008-04-04
 ****************************************************************************/

basename = "01-macro_nonmixed"
countdir = "dat/" & basename & "/"
vizdir   = "./viz_dat/" & basename & "/"
sprintf(seed, "%04g", SEED)

dt = 1e-6
iterations=500000

INCLUDE_FILE="parameters.mdl"

TIME_STEP = dt
ITERATIONS = iterations
GRID_DENSITY = 15000
ACCURATE_3D_REACTIONS = FALSE

PARTITION_X = [[-box_size/2 - 0.0001 TO box_size/2 + 0.0001 STEP partition_size]]
PARTITION_Y = [[-box_size/2 - 0.0001 TO box_size/2 + 0.0001 STEP partition_size]]
PARTITION_Z = [[-box_size/2 - 0.0001 TO box_size/2 + 0.0001 STEP partition_size]]
MEMORY_PARTITION_X = mem_partition_size
MEMORY_PARTITION_Y = mem_partition_size
MEMORY_PARTITION_Z = mem_partition_size

DEFINE_MOLECULES {
  camkii_subunit  { DIFFUSION_CONSTANT_3D = 0 }
  camkii_subunitB { DIFFUSION_CONSTANT_3D = 0 }
  cam             { DIFFUSION_CONSTANT_3D = cam_diffusion }
}

DEFINE_COMPLEX_MOLECULE camkii {
  NUMBER_OF_SUBUNITS = [2, 6]
  SUBUNIT[1:2, 1:6] = camkii_subunit
  SHAPE {
    SUBUNIT[1, 1] = [-.10,  .0000,  .05]
    SUBUNIT[1, 2] = [-.05,  .0866,  .05]
    SUBUNIT[1, 3] = [ .05,  .0866,  .05]
    SUBUNIT[1, 4] = [ .10,  .0000,  .05]
    SUBUNIT[1, 5] = [ .05, -.0866,  .05]
    SUBUNIT[1, 6] = [-.05, -.0866,  .05]
    SUBUNIT[2, 1] = [-.10,  .0000, -.05]
    SUBUNIT[2, 2] = [-.05,  .0866, -.05]
    SUBUNIT[2, 3] = [ .05,  .0866, -.05]
    SUBUNIT[2, 4] = [ .10,  .0000, -.05]
    SUBUNIT[2, 5] = [ .05, -.0866, -.05]
    SUBUNIT[2, 6] = [-.05, -.0866, -.05]
  }
  SUBUNIT_RELATIONSHIPS {
    ring_negative = [ 0, -1]
    ring_positive = [ 0, +1]
    dimer_partner = [+1,  0]
  }
  RATE_RULES {
    coop_binding_rate {
      dimer_partner != camkii_subunit : fw_rate * TMP_spec_F
      DEFAULT : fw_rate
    }
    coop_unbinding_rate {
      dimer_partner != camkii_subunit : bw_rate * TMP_spec_B
      DEFAULT : bw_rate
    }
  }
}

DEFINE_REACTIONS {
  (camkii_subunit) + cam <-> (camkii_subunitB)   [> COMPLEX_RATE camkii coop_binding_rate,
                                                  < COMPLEX_RATE camkii coop_unbinding_rate]
}

mainbox OBJECT {
  b BOX {
    CORNERS = [ -box_size / 2, -box_size / 2, -box_size / 2 ], [ box_size / 2, box_size / 2, box_size / 2 ]
  }
}

INSTANTIATE world OBJECT
{
  mainbox OBJECT mainbox {
  }

  rs1 RELEASE_SITE {
    SHAPE = world.mainbox.b[ALL]
    MOLECULE = camkii
    NUMBER_TO_RELEASE = num_complexes
  }

  rs2 RELEASE_SITE {
    SHAPE = SPHERICAL
    MOLECULE = cam
    LOCATION = pt_of_release
    SITE_RADIUS = 0
    NUMBER_TO_RELEASE = num_cam
  }
}

REACTION_DATA_OUTPUT
{
    OUTPUT_BUFFER_SIZE = 1000
    STEP = 1e-5
    HEADER = "# "
    {
      COUNT[SUBUNIT{camkii : camkii_subunit  [dimer_partner == camkii_subunit]}, WORLD] : "00",
      COUNT[SUBUNIT{camkii : camkii_subunit  [dimer_partner != camkii_subunit]}, WORLD] : "01",
      COUNT[SUBUNIT{camkii : camkii_subunitB [dimer_partner == camkii_subunit]}, WORLD] : "10",
      COUNT[SUBUNIT{camkii : camkii_subunitB [dimer_partner != camkii_subunit]}, WORLD] : "11",
      COUNT[camkii_subunit, WORLD]                                                      : "us",
      COUNT[camkii_subunitB, WORLD]                                                     : "bs",
      COUNT[cam, WORLD]                                                                 : "cam"
    }  => countdir & "counts.dat"
}