I am a math student trying to implement a simple force field in python for determining the energy of a small molecule and trying to minimize it in order to find a stable configuration. I was wondering if there is a prefered system of units to be used am currently using the SI system of units however when i calculate Van der waals energies the numbers are quite difficult to deal with . Any reference of basic force fields would be much appreciated
1
-
1$\begingroup$ You usually want to use atomic units (Hartree are very common) for such calculations, since it makes the math easier and there are less possibilities for mistakes, bugs and stuff like that. See: muchomas.lassp.cornell.edu/P480/Notes/dft/node13.html You can always translate those into SI after the calculation. $\endgroup$– DSVAFeb 11, 2017 at 14:20
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
Atomic units are the standard in almost any ab-initio (LAMMPS, NWChem, Orca, QChem) and/or classical force field (LAMMPS, MPMC, CHARMM, AMBER, GROMACS etc.) modeling code because you end up computing numbers much closer to 1.
Things can get bad very quickly for small values in python. To prevent this, I would use atomic units, and wherever arithmetic error occurs, try using the decimal type, instead of default float, via float -> string -> decimal:
from decimal import *
energy = Decimal(str(3.4280284e-19))
Here is a sample file I've used for atomic units in python:
pi = 3.14159265359
kb = 1.38064852e-23 # boltzmann / J/K
ab = 5.2917721092e-11 # bohr radius / meters
hbar = 1.054571726e-34 # planck / Js
h = hbar*2*pi
me = 9.10938291e-31 # electron rest mass / kg
e = 1.602176565e-19 # electron charge / C
ke = 8.9875517873681e9 # coulomb constant / kgm^3 s^-2 C^-2
Eh = 4.35974417e-18 # hartree energy / J
c = 299792458 # speed of light, m/s
cs = 340.29 # speed of sound, m/s
t = hbar / Eh # time / s
v = ab*Eh/hbar # velocity m/s
F = Eh/ab #force
T = Eh/kb # temp
P = Eh/(ab**3) # pressure
Ef = Eh/(e*ab) # electric field
Ep = Eh/e # electric potential
u = e*ab # electric dipole moment
$\begingroup$
$\endgroup$
This very much depends on whether you're trying to calculate the total energy, or whether you want to know the relative stability of different structures. The most convenient choice would be akin to Lennard-Jones dimensionless units (https://en.wikipedia.org/wiki/Lennard-Jones_potential#Dimensionless_.28reduced.29_units), if you're using a relatively simple model. If you're after real values, and have some empirical potential, you may want to consider atomic units (Hartrees for energy, Bohr radii for length, atomic time units for time), or alternatively electron volts, angstroms, femtoseconds, etc.
