Skip to content

Commit

Permalink
Merge pull request #4 from m3g/release-1.3.2
Browse files Browse the repository at this point in the history 
Release 1.3.2
  • Loading branch information
@lmiq
lmiq authored Dec 21, 2023
2 parents b1568b8 + 600632d commit 4a1d006
Show file tree
Hide file tree
Showing 4 changed files with 132 additions and 65 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "MolSimToolkit"
uuid = "054db54f-6694-444d-9bbb-e9ecdbfe77be"
authors = ["Leandro Martinez <lmartine@unicamp.br> and contributors"]
version = "1.3.2-DEV"
version = "1.3.2"

[deps]
AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
Expand Down
6 changes: 3 additions & 3 deletions docs/src/system_setup.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ and sizes necessary to build different systems. Currently (as of version 1.2.0)
supports the construction of systems of a solute solvated by a mixture of two solvents.

!!! compat
The functionality described here is available in MolSimToolkit version 1.3.0.
The functionality described here is available in MolSimToolkit version 1.3.2.

### How to use it

Expand Down Expand Up @@ -60,10 +60,10 @@ The concentration units can be one of `"mol/L"` (molarity), `"x"` (molar fractio
to be in `g/mL`.

!!! tip
The density table can be converted among different units with the function `convert_density_table`,
The density table can be converted among different units with the function `convert_density_table!`,
which acts on the `SystemBox` object. For example:
```julia-repl
julia> convert_density_table(system, "mol/L")
julia> convert_density_table!(system, "mol/L")
```

Finally, we can generate an input file for `Packmol` with:
Expand Down
120 changes: 81 additions & 39 deletions src/PackmolInputCreator/PackmolInputCreator.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ using ..MolSimToolkit: center_of_mass
using PDBTools
using TestItems

export convert_concentration, convert_density_table
export convert_concentration, convert_density_table!
export density_pure_solvent, density_pure_cossolvent
export write_packmol_input
export SolutionBoxUSC
Expand All @@ -15,15 +15,15 @@ density_pure_cossolvent(system::SolutionBox) = system.density_table[end, 2]

const PackmolInputCreatorDirectory = @__DIR__

@kwdef mutable struct SolutionBoxUSC <: SolutionBox
concentration_units::String = "x"
mutable struct SolutionBoxUSC <: SolutionBox
concentration_units::String
solute_pdbfile::String
solvent_pdbfile::String
cossolvent_pdbfile::String
density_table::Matrix{Float64}
solute_molar_mass::Float64 = mass(readPDB(solute_pdbfile))
solvent_molar_mass::Float64 = mass(readPDB(solvent_pdbfile))
cossolvent_molar_mass::Float64 = mass(readPDB(cossolvent_pdbfile))
solute_molar_mass::Float64
solvent_molar_mass::Float64
cossolvent_molar_mass::Float64
end

"""
Expand All @@ -41,7 +41,49 @@ The concentration units of the density table can be provided explicitly and
are assumed by default to be the molar fraction, `x`, of the cossolvent.
"""
function SolutionBoxUSC end
function SolutionBoxUSC(;
solute_pdbfile::String,
solvent_pdbfile::String,
cossolvent_pdbfile::String,
density_table::Matrix{Float64},
concentration_units::Union{Nothing,String} = nothing,
)
if isnothing(concentration_units)
concentration_units = "x"
@warn "Concentration units not provided, assuming molar fraction." _file=nothing _line=nothing
end
if !(density_table[begin, 1] == 0.0)
throw(ArgumentError("First line of density table be the density of pure solvent, with cossolvent concentration equal to 0.0"))
end
if concentration_units in ("x", "vv", "mm") && !(density_table[end, 1] == 1.0)
throw(ArgumentError("Last line of density table be the density of pure cossolvent, with cossolvent concentration equal to 1.0"))
end
solute_molar_mass = mass(readPDB(solute_pdbfile))
solvent_molar_mass = mass(readPDB(solvent_pdbfile))
cossolvent_molar_mass = mass(readPDB(cossolvent_pdbfile))
system = SolutionBoxUSC(
concentration_units,
solute_pdbfile,
solvent_pdbfile,
cossolvent_pdbfile,
density_table,
solute_molar_mass,
solvent_molar_mass,
cossolvent_molar_mass,
)
if concentration_units == "mol/L"
x = convert_concentration(system, density_table[end, 1], "mol/L" => "x")
if !(x ≈ 1)
throw(ArgumentError(chomp("""
Conversion of concentration of last line into molar fraction gives: $x which is different from 1.0.
The last line of the density_table must contain the density of pure cossolvent.
""")))
end
end
return system
end

function unit_name(u::String)
u == "mol/L" && return "molarity"
Expand Down Expand Up @@ -101,8 +143,7 @@ fixrange(x) = x < 0 ? 0 : (x > 1 ? 1 : x)
convert_concentration(
system::SolutionBoxUSC,
input_concentration,
units;
density, # of the solution, required when converting from/to mol/L
units
)
Convert concentration from one unit to another. The input
Expand Down Expand Up @@ -133,33 +174,29 @@ where `system` is a `SolutionBoxUSC` object, and `55.5` is the molarity.
function convert_concentration(
system::SolutionBoxUSC,
input_concentration::Real,
units::Pair{String,String};
density::Union{Real,Nothing} = nothing, # of the solution
units::Pair{String,String}
)

(; solvent_molar_mass, cossolvent_molar_mass) = system

# If the units didn't change, just return the input concentrations
first(units) == last(units) && return input_concentration

ρ = density # density of the solution
# Obtain density of the solution by interpolation
if first(units) != system.concentration_units
input_concentration_units = system.concentration_units
convert_density_table!(system, first(units))
ρ = interpolate_concentration(system, input_concentration)
convert_density_table!(system, input_concentration_units)
else
ρ = interpolate_concentration(system, input_concentration)
end

ρc = density_pure_cossolvent(system) # density of the pure cossolvent
ρw = density_pure_solvent(system) # density of the pure solvent
Mw = solvent_molar_mass # molar mass of solvent
Mc = cossolvent_molar_mass # molar mass of the cossolvent

# Check if density of solution was provided, in the case of molar
# fraction conversions. Also, to converto to other units from molarity, we need ρ
if isnothing(ρ)
if last(units) == "mol/L" || last(units) == "mol/L"
throw(ArgumentError(
"""
Density of solution is required to convert from/to molarity.
Use the optional keyword argument `density` to provide it.
"""))
end
end

# nc and nw are the molar concentrations
if first(units) == "vv"
if !(0 <= input_concentration <= 1)
Expand Down Expand Up @@ -247,7 +284,7 @@ function convert_concentration(
end

"""
convert_density_table(system::SolutionBoxUSC, target_units)
convert_density_table!(system::SolutionBoxUSC, target_units)
Converts the density table of the system from one unit to another. Returns the
input `system` with the density table converted to the new units.
Expand All @@ -257,22 +294,23 @@ The target units may be one of: `"mol/L"`, `"x"`, `"vv"`, `"mm"`.
## Example
```julia
convert_density_table(system, "mol/L")
convert_density_table!(system, "mol/L")
```
"""
function convert_density_table(
function convert_density_table!(
system::SolutionBoxUSC,
target_units::String;
)
current_units = system.concentration_units
density_table = system.density_table
for irow in eachindex(eachrow(density_table))
cin = density_table[irow, 1]
ρ = density_table[irow, 2]
cout = convert_concentration(system, cin, current_units => target_units; density = ρ)
density_table[irow, 1] = cout
new_density_table = copy(system.density_table)
for irow in eachindex(eachrow(new_density_table))
cin = new_density_table[irow, 1]
ρ = new_density_table[irow, 2]
cout = convert_concentration(system, cin, current_units => target_units)
new_density_table[irow, 1] = cout
end
system.density_table .= new_density_table
system.concentration_units = target_units
return system
end
Expand Down Expand Up @@ -333,10 +371,10 @@ function write_packmol_input(
ρ = interpolate_concentration(system, concentration)

# Obtain the concentration in all units, for testing
c_x = convert_concentration(system, concentration, cunit => "x"; density = ρ)
c_vv = convert_concentration(system, concentration, cunit => "vv"; density = ρ)
cc_mol = convert_concentration(system, concentration, cunit => "mol/L"; density = ρ)
cc_mm = convert_concentration(system, concentration, cunit => "mm"; density = ρ)
c_x = convert_concentration(system, concentration, cunit => "x")
c_vv = convert_concentration(system, concentration, cunit => "vv")
cc_mol = convert_concentration(system, concentration, cunit => "mol/L")
cc_mm = convert_concentration(system, concentration, cunit => "mm")

# aliases for clearer formulas
ρc = density_pure_cossolvent(system)
Expand Down Expand Up @@ -368,7 +406,11 @@ function write_packmol_input(
nc = round(Int, cc * vs)

# Number of solvent molecules
nw = round(Int, nc * (1 - c_x) / c_x)
if nc != 0
nw = round(Int, nc * (1 - c_x) / c_x)
else
nw = round(Int, vs * (ρw/CMV) / Mw)
end

# Final density of the solution (not inclusing solute volume)
ρ = CMV * (Mc * nc + Mw * nw) / vs
Expand Down Expand Up @@ -469,7 +511,7 @@ function write_packmol_input(
"""))

if debug
return nw, nc, l
return nw, nc, 2*l
else
return nothing
end
Expand Down
69 changes: 47 additions & 22 deletions src/PackmolInputCreator/test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,16 +13,16 @@
mw = 55.508250191225926
for x in (0.0, 0.2, 0.5, 0.7, 1.0)
@test convert_concentration(system, x, "x" => "vv") ≈ x atol = 1e-3
@test convert_concentration(system, x, "x" => "mol/L"; density = 1.0) ≈ x * mw atol = 1e-3
@test convert_concentration(system, x, "x" => "mol/L") ≈ x * mw atol = 1e-3
@test convert_concentration(system, x, "x" => "mm") ≈ x atol = 1e-3

@test convert_concentration(system, x, "vv" => "x") ≈ x atol = 1e-3
@test convert_concentration(system, x, "vv" => "mol/L"; density = 1.0) ≈ x * mw atol = 1e-3
@test convert_concentration(system, x, "vv" => "mol/L") ≈ x * mw atol = 1e-3
@test convert_concentration(system, x, "vv" => "mm") ≈ x atol = 1e-3

@test convert_concentration(system, x * mw, "mol/L" => "x"; density = 1.0) ≈ x atol = 1e-3
@test convert_concentration(system, x * mw, "mol/L" => "vv"; density = 1.0) ≈ x atol = 1e-3
@test convert_concentration(system, x * mw, "mol/L" => "mm"; density = 1.0) ≈ x atol = 1e-3
@test convert_concentration(system, x * mw, "mol/L" => "x") ≈ x atol = 1e-3
@test convert_concentration(system, x * mw, "mol/L" => "vv") ≈ x atol = 1e-3
@test convert_concentration(system, x * mw, "mol/L" => "mm") ≈ x atol = 1e-3
end

# system with ideal solution
Expand Down Expand Up @@ -50,18 +50,40 @@
# Concentration conversions
for x in (0.0, 0.2, 0.5, 0.7, 1.0)
@test convert_concentration(system, x, "x" => "vv") ≈ vv(x)
@test convert_concentration(system, x, "x" => "mol/L"; density = ρ(x)) ≈ mx(x)
@test convert_concentration(system, x, "x" => "mol/L") ≈ mx(x)
@test convert_concentration(system, x, "x" => "mm") ≈ x

@test convert_concentration(system, vv(x), "vv" => "x") ≈ x
@test convert_concentration(system, vv(x), "vv" => "mol/L"; density = ρ(x)) ≈ mx(x)
@test convert_concentration(system, vv(x), "vv" => "mol/L") ≈ mx(x)
@test convert_concentration(system, vv(x), "vv" => "mm") ≈ x

@test convert_concentration(system, mx(x), "mol/L" => "x"; density = ρ(x)) ≈ x
@test convert_concentration(system, mx(x), "mol/L" => "vv"; density = ρ(x)) ≈ vv(x)
@test convert_concentration(system, mx(x), "mol/L" => "mm"; density = ρ(x)) ≈ x
@test convert_concentration(system, mx(x), "mol/L" => "x") ≈ x
@test convert_concentration(system, mx(x), "mol/L" => "vv") ≈ vv(x)
@test convert_concentration(system, mx(x), "mol/L" => "mm") ≈ x
end

# Water as cossolvent
dw = [
0.0 0.7906
0.2214 0.8195
0.3902 0.845
0.5231 0.8685
0.6305 0.8923
0.7191 0.9151
0.7934 0.9369
0.8566 0.9537
0.911 0.9685
0.9584 0.982
1.0 0.9981
]
system = SolutionBoxUSC(
solute_pdbfile = "$test_dir/data/poly_h.pdb",
solvent_pdbfile = "$test_dir/data/ethanol.pdb",
cossolvent_pdbfile = "$test_dir/data/water.pdb",
density_table = dw
)
@test convert_concentration(system, 1.0, "x" => "mol/L") ≈ 55.40278451586260

end

@testitem "Write packmol input" begin
Expand All @@ -82,7 +104,7 @@ end
Mc = 1000 * density_pure_cossolvent(system) / system.cossolvent_molar_mass # mol / L pure ethanol

# Test concentration conversions for real data
@test convert_concentration(system, 1.0, "x" => "mol/L"; density = density_pure_cossolvent(system)) ≈ Mc
@test convert_concentration(system, 1.0, "x" => "mol/L") ≈ Mc

mm = 0.3997224931406948
vv = 0.45671854335897716
Expand All @@ -91,40 +113,43 @@ end
ρ = 0.9369

@test convert_concentration(system, x, "x" => "vv") ≈ vv
@test convert_concentration(system, x, "x" => "mol/L"; density = ρ) ≈ M
@test convert_concentration(system, x, "x" => "mol/L") ≈ M
@test convert_concentration(system, x, "x" => "mm") ≈ mm

@test convert_concentration(system, vv, "vv" => "x") ≈ x
@test convert_concentration(system, vv, "vv" => "mol/L"; density = ρ) ≈ M
@test convert_concentration(system, vv, "vv" => "mol/L") ≈ M
@test convert_concentration(system, vv, "vv" => "mm") ≈ mm

@test convert_concentration(system, M, "mol/L" => "x"; density = ρ) ≈ x
@test convert_concentration(system, M, "mol/L" => "mm"; density = ρ) ≈ mm
@test convert_concentration(system, M, "mol/L" => "vv"; density = ρ) ≈ vv
@test convert_concentration(system, M, "mol/L" => "x") ≈ x
@test convert_concentration(system, M, "mol/L" => "mm") ≈ mm
@test convert_concentration(system, M, "mol/L" => "vv") ≈ vv

@test convert_concentration(system, mm, "mm" => "x") ≈ x
@test convert_concentration(system, mm, "mm" => "mol/L"; density = ρ) ≈ M
@test convert_concentration(system, mm, "mm" => "mol/L") ≈ M
@test convert_concentration(system, mm, "mm" => "vv") ≈ vv

tmp_input_file = tempname()

convert_density_table(system, "x")
convert_density_table!(system, "x")
r1 = write_packmol_input(system; concentration = 0.0, box_sides=[120,120,120], input = tmp_input_file, debug = true)
@test r1[1] == 57322

r1 = write_packmol_input(system; concentration = 0.5, margin = 20.0, input = tmp_input_file, debug = true)
@test isfile(tmp_input_file)

convert_density_table(system, "mol/L")
convert_density_table!(system, "mol/L")
r2 = write_packmol_input(system; concentration = 13.488667939471432, margin = 20.0, input = tmp_input_file, debug = true)
@test all(isapprox.(r2,r1,rtol=0.005))

convert_density_table(system, "vv")
convert_density_table!(system, "vv")
r3 = write_packmol_input(system; concentration = 0.7635032204047275, margin = 20.0, input = tmp_input_file, debug = true)
@test all(isapprox.(r3,r1,rtol=0.005))

convert_density_table(system, "mm")
convert_density_table!(system, "mm")
r4 = write_packmol_input(system; concentration = 0.7188817400010237, margin = 20.0, input = tmp_input_file, debug = true)
@test all(isapprox.(r4,r1,rtol=0.005))

convert_density_table(system, "mol/L")
convert_density_table!(system, "mol/L")
system.concentration_units = "x"
@test_throws ArgumentError write_packmol_input(system; concentration = 0.5, margin = 20.0, input = tmp_input_file, debug = true)
system.concentration_units = "mol/L"
Expand Down

2 comments on commit 4a1d006

@lmiq
Copy link
Member Author

@lmiq lmiq commented on 4a1d006 Dec 21, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

Release notes:

  • fixed bugs in input creator

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/97573

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.3.2 -m "<description of version>" 4a1d00636dfa706465a23c8f707bc934c7422566
git push origin v1.3.2

Please sign in to comment.