Error when computing energy grid for rotatable molecules in a charged crystal

[qlx17]

I would like to report two potential bugs in grid.jl, specifically on lines 342 and 347. These issues arise when computing the energy grid for molecules that can rotate within a charged crystal.

While computing the energy grid for CO₂ in a charged crystal using the following code:

include("../src/PorousMaterials.jl")
using .PorousMaterials
# using PorousMaterials

xtal = Crystal("zif71_bogus_charges.cif")
strip_numbers_from_atom_labels!(xtal)
molecule = Molecule("CO2")
ljforcefield = LJForceField("UFF", r_cutoff=14.0, mixing_rules="Lorentz-Berthelot")
grid = energy_grid(xtal, molecule, ljforcefield, 
                   resolution=1.0, units=:kJ_mol, temperature=300.0, n_rotations=2)

I encountered the following error:

WARNING: replacing module PorousMaterials.
WARNING: using PorousMaterials.LJForceField in module Main conflicts with an existing identifier.
WARNING: using PorousMaterials.energy_grid in module Main conflicts with an existing identifier.
WARNING: using PorousMaterials.Molecule in module Main conflicts with an existing identifier.
        Ewald summation parameters:
                k-replication factors: 7 7 7
                Ewald convergence param. α = 0.219398
                short-range cutoff radius (Å): 14.000000
                823 kvectors
                (calculated with specified precision 1.000000e-06)
Computing energy grid of CO2 in zif71_bogus_charges.cif
        Regular grid (in fractional space) of 30 by 30 by 30 points superimposed over the unit cell.
        2 molecule rotations per grid point with temperature 300.000000 K.
ERROR: LoadError: MethodError: no method matching random_rotation!(::Molecule{Frac})

Closest candidates are:
  random_rotation!(::Molecule{Frac}, ::Box)
   @ Main.PorousMaterials ~/MyGitDownload/PorousMaterials_qlx.jl/src/molecule.jl:269
  random_rotation!(::Molecule{Cart})
   @ Main.PorousMaterials ~/MyGitDownload/PorousMaterials_qlx.jl/src/molecule.jl:255

Stacktrace:
 [1] energy_grid(crystal::Crystal, molecule::Molecule{Cart}, ljforcefield::LJForceField; resolution::Float64, n_rotations::Int64, temperature::Float64, units::Symbol, center::Bool, verbose::Bool)
   @ Main.PorousMaterials ~/MyGitDownload/PorousMaterials_qlx.jl/src/grid.jl:342
 [2] top-level scope
   @ ~/MyGitDownload/PorousMaterials_qlx.jl/test/PES_MWE.jl:9
 [3] include(fname::String)
   @ Base.MainInclude ./client.jl:478
 [4] top-level scope
   @ REPL[2]:1

This error seems to originate from a bug on line 342 of grid.jl. To troubleshoot, I modified this line to random_rotation!(molecule, crystal.box), and ran the code again.

However, I then encountered a different error:

Ewald summation parameters:
                k-replication factors: 7 7 7
                Ewald convergence param. α = 0.219398
                short-range cutoff radius (Å): 14.000000
                823 kvectors
                (calculated with specified precision 1.000000e-06)
Computing energy grid of CO2 in zif71_bogus_charges.cif
        Regular grid (in fractional space) of 30 by 30 by 30 points superimposed over the unit cell.
        2 molecule rotations per grid point with temperature 300.000000 K.
ERROR: LoadError: type PotentialEnergy has no field coulomb
Stacktrace:
 [1] setproperty!(x::PotentialEnergy, f::Symbol, v::Float64)
   @ Base ./Base.jl:38
 [2] energy_grid(crystal::Crystal, molecule::Molecule{Cart}, ljforcefield::LJForceField; resolution::Float64, n_rotations::Int64, temperature::Float64, units::Symbol, center::Bool, verbose::Bool)
   @ Main.PorousMaterials ~/MyGitDownload/PorousMaterials_qlx.jl/src/grid.jl:348
 [3] top-level scope
   @ ~/MyGitDownload/PorousMaterials_qlx.jl/test/PES_MWE.jl:9
 [4] include(fname::String)
   @ Base.MainInclude ./client.jl:478
 [5] top-level scope
   @ REPL[2]:1

This error appears to be due to line 347 of grid.jl, where a non-existent property coulomb of type PotentialEnergy is referenced. I believe this line should be modified to:

energy.es = total(electrostatic_potential_energy(crystal, molecule,
                                                     eparams, eikr))

Could you please look into these issues?

[eahenle]

Thanks @qlx17 for the report! I will take a look, but @SimonEnsemble (currently travelling) may need to be involved in the fix.