Take snapshots of molecule positions during GCMC and record adsorbate density grid

[SimonEnsemble]

@ahyork

Looks great! Before merging this PR:

• add to documentation the new optional arguments to gcmc_simulation()
• snapshots to more informative write_adsorbate_snapshots since now it is only writing them to file.
• xyz_snapshot_file = "" is initialized as a string. Doesn't that introduce type instability? Does xyz_snapshot_file = IOStream work better, initializing an empty IOStream?
• comment=filename_comment * "adsorbate_positions" switch order so comment is always at end.
• " with each point being 0.1 units apart." not true anymore
• if (outer_cycle % snapshot_frequency == 0) && (outer_cycle > n_burn_cycles) for speed is it better to reverse these? I imagine > is faster than a %?
• write_xyz_to_file can you change to write_xyz to be consistent with the other write_xyz? i.e. overload it so one less function name to remember.
• Make snapshot_frequency=1 by default.
• "- snapshot::Bool: Whether this is the name for a snapshot file" is a vestige of an old version of gcmc_result_savename?
• "Takes the fractional coordinates of a point in a unit cell and returns the indices for storing that data inside a Grid object" to be more specific, change to: "returns the indices of the voxel in which it falls when a unit cube is partitioned into a regular grid of n_pts[1] by n_pts[2] by n_pts[3] voxels.
• write_xyz_to_file(molecules, xyz_file) doc string does not match function. specify in doc string it is writing cartesian not fractional
• atoms_in_molecules is redundant given type assertion, how about n_atoms(molecules::Array{Molecules})?
• in update_density!, the problem here is that for CO2, it would increment the voxel for both C and O then you'd get a density grid that is confusingly keeping track of C and O. I think we should put the species granularity on the atom within the molecule (not the molecule itself), so e.g. the user would pass :C to keep track of :C within CO2. so if molecule.atoms.species[i] == species I think should be the condition for writing

[SimonEnsemble]

one more thing: the GCMC prints off a bunch of stuff in the beginning.

if snapshots are being taken, can you printf "writing snapshots of adsorption positions every %d cycles (after burn cycles)\n"

if density grid is being written, can you printf "adsorption spatial probability density grid written with spacing %.3f Angstrom every %d cycles (after burn cycles)"

[ahyork]

I went through all suggestions and implemented them. The branch passes all tests as well.

This should be ready to merge, I am going to update the docs site so the new function appear there as well.

[SimonEnsemble]

• I just realized snapshot_frequency::Int: is a misleading name for the variable since high frequency implies often. how about snapshot_period::Int, and it is the period of the snapshots expressed in number of cycles (not samples, which is different since there are many samples in one cycle). would be nice to change sample_frequency and checkpoint_frequency to *_period as well.
• I discovered a bug! So currently you make the grid before the repfactors are applied to the framework. This should be done after framework = replicate(framework, repfactors) since the fractional coordinates will change. This only affects frameworks that are replicated, but still will be a bug for these. so move the snapshots stuff after that line. also plz done with a comment # snapshot / density grid set up for clarity.
• close(xyz_snapshot_file) with new intialization, do u need to close it every time now?
• num_atoms = 0; semi-colon not required
• @test all(xf_to_id(n_pts, [0.0001, 0.24, 0.26]) .== [1, 2, 3]) for another test
• any test for update_density!? how about:

  # test update_density! by placing CO2 strategically to be in voxel (1, 2, 4)
  box = UnitCube()
  n_pts = (4, 4, 4)
  density_grid = Grid(UnitCube(), n_pts, zeros(n_pts...), :inverse_A3, [0.0, 0.0, 0.0])
  molecule = Molecule("CO2")
  translate_to!(molecule, [0.24, 0.26, 0.99])
  update_density!(density_grid, [molecule], :C_CO2)
  @test isapprox(sum(density_grid.data), 1.0)
  @test isapprox(density_grid.data[1, 2, 4], 1.0)

[SimonEnsemble]

@printf("\tAdsorption spatial probability density grid written with spacing %.3f Angstroms every %d cycles (after burn cycles)\n", density_grid_dx, snapshot_frequency) also can you print out density_grid.n_pts to alert the user if the grid is so fine that it is going to be a huge file?

[ahyork]

I am wrapping up with these edits soon, but I was just thinking about changing the *_frequency to *_period. This is going to break a lot of backwards compatibility because these are optional arguments. People will need to change their scripts in order to keep this working.

I think it is a good idea to do this, but we should release the next version as 1.0.0 since PorousMaterials.jl is no longer going through rapid development and the API won't be compatible with scripts written for older versions.

[ahyork]

I found an issue when writing a test for the O_CO2 in CO2. One of the atoms is outside the box, and so it cannot be placed at a negative index. I am going to see if I can use NearestImage to get the indices to wrap.

[SimonEnsemble]

:fireworks: hooray for tests! Didn't think of that...

[ahyork]

I think we should leave renaming *_frequency to *_period for a later time.

I agree that informative parameter names are important, but that small change breaks some backwards compatibility.

We should leave it out for now, but we should go through and decide what 'stable' PorousMaterials.jl looks like. Then we can fix small errors like this without creating a chain of incompatible releases

[SimonEnsemble]

sure, that is fine with me!

[ahyork]

This should be good to go now! Take a look at the new xf_to_id. It checks to see if it is above and below because it has to treat them with different cases.

[SimonEnsemble]

thanks @ahyork

• I allowed user to pass in which species to keep track of in the case that there are multiple atoms in a molecule. (before it was keeping track of molecule.species which might not even be a species atom; e.g. the grid would be empty for Molecule("CO2"))
• the xf_to_id! currently looks like it would modify the coordinates of the molecule? to avoid changing the coordinates of the molecule since arrays are passed by reference, I instead changed it so that periodic boundary conditions are applied to voxel_id

do these changes look ok to you?

[ahyork]

Those changes look good to me, I am going to add these new functions to the docs and then merge unless there is anything else you think we should add?