Minimizing energy of a molecule

Question

I am minimizing energy of butane on the CHARMM Force Field and I am confused on how to go about getting a uniform unit for each different PE as in the potential energy for bond, angle, torsion, VDW, ES of butane. For instance if you compute the angle PE you will be using angles vs bond length PE which will be in Ångstrom. I am confused on how to set everything within the same units or parameter.

Also I am programming this in python using numerical methods and will not be using a simulator.

Answer 1

First I only wanted to comment, but then it got too long... But I wanted to pass on this info... this isn't a real answer...

Every simulator I know uses atomic units. I would check out the NAMD and CHARMM manual for their units, I imagine they use nanometer for distance and ultimately kJ/mol for energy. You will need a conversion for the electrostatic since you always input +/- charge, where charge is usually between 0 and 1, but units of charge are actually Coulomb. The conversion factor for this is ~138.9

This is from my own code... qq_convert is the electrostatic conversion I mentioned above. I use nm for all distances. I try to follow the same procedure as Gromacs.

const na = 6.02214129e23      # mol^-1
const R = 8.3144598e-3        # kJ mol^-1 K^-1
const kb = R / na             # kJ K^-1
const e = 1.6021765653e-19     # C
ϵ₀ = 8.854187817e-12 # C²/(J ⋅ m)
ϵ₀ *= 1e-9  # C²/(J ⋅ nm)
ϵ₀ *= 1000   # C²/(kJ ⋅ nm)
const qq_convert = e^2 / ϵ₀ * na / 4.0 / pi # 138.935458  # kJ mol^-1 nm e^-2

after each electrostatic interaction I multiply the energy by qq_convert

For instance the real part contribution of the EWALD summation is

@inline function ERealSingle(diff::MVector{3,Float64},qqCut²::Float64, coord1::SVector{3,Float64},
        coord2::SVector{3,Float64},q1::Float64, q2::Float64, kappa::Float64,
        overflow²::Float64, L::SVector{3,Float64})

    @inbounds for i=1:3 
        diff[i] = vector1D(coord1[i],coord2[i],L[i]) # mirror image separation
    end

    rij_sq = diff[1]*diff[1] + diff[2]*diff[2] + diff[3]*diff[3]
    rij = sqrt(rij_sq)

    if rij_sq < qqCut²
        return q1 * q2 * erfc( kappa * rij ) / rij
    else
        return 0.0
    end
end

This is for a single charge interacting with another charge... I loop through all charges etc and then I multiply the summed energy by qq_convert.

I highly recommend checking your code against these SPC/E water benchmarks offered by NIST... really saved my bacon having a test set to duplicate... NIST Test Set. This set allows you to test your Lennard-Jones, all parts of the EWALD Summation and also tail corrections.

Lastly, I recommend using the Julia programming language... it is much faster than Python... Python can almost keep up if you use the Numba JIT decorators, but Julia has other perks.

Edit:

It looks to me like charmm uses Angstrom as you have said, so that is the only units of length you should be using. It would be dangerous for them to deviate.

In my experience with Gromacs, the parameters of each potential energy function will have units such that everything cancels out leaving you with $\ce{kJ/mol}$ (or $\ce{kcal/mol}$ in CHARMM I assume). In Gromacs all I need to do is make sure all units of length are nm. I imagine in CHARMM it is the same, but for Angstrom. Look at the angle parameter units, they do have units and it will work out that the energy is $\ce{kcal/mol}$ or else I say leave CHARMM.