[gmx-users] Force Constants and Unit Systems

[Richard Broadbent]

On Fri, 2012-05-18 at 09:40 +0100, Lara Bunte wrote: 
> Hi
> 
> I have two questions left: 
> 
> 
> 1.)
> You wrote 
> 
> 
> >If your term in question is an angle potential, then the force constant should indeed have units of energy
> 
> Could you please explain this? Why is it here consistent to have Energy as unit for a force constant?
I believe he meant that the angular force coefficient should not have
length dependence in gromacs it should be in kJ/(mol rad^2). 

The length dependence when calculating a harmonic angular force comes
from differentiating the angle theta_{ijk} with respect to the positions
this is done internally in the code. 
> 
> 2.)
> The force constant is divided by mol. What is the reason for this? In example Energy/length^2 is the force constant, why do we consider Energy/(mol * lenght^2) ?
A typical bond stretching coefficient is around 10^-10 J/m^2. If a
simulation package was built using this unit system it would hit many
problems with the finite precision of numbers in a computer. Therefore,
we pick units which make the numbers we are likely to encounter of order
1. This is done in this unit system by dividing by avagadro's constant
(which is dimensionless). This has the added benefit that it is the unit
system which chemists often use as they tend to do experiments on around
a mole of the substance.


Hope that helps,

Richard 
> 
> Thanks for helping me
> Greetings
> Lara
> 
> 
> 
> 
> 
> ----- Ursprüngliche Message -----
> Von: Justin A. Lemkul <jalemkul at vt.edu>
> An: Lara Bunte <lara.bunte at yahoo.de>; Discussion list for GROMACS users <gmx-users at gromacs.org>
> CC: 
> Gesendet: 20:37 Donnerstag, 17.Mai 2012
> Betreff: Re: [gmx-users] Force Constants and Unit Systems
> 
> 
> 
> On 5/17/12 1:33 PM, Lara Bunte wrote:
> > Hi
> > 
> >> Therefore either they have a potential of the form 1/[Length] or they
> >> weren't using the term correctly.
> > 
> > But a 1/[lenght] potential, which is a coulomb potential makes no sense for
> > springs, that have a quadratic potential, like V(x) = 1/2 * k * x^2 of a
> > harmonic oscillator.
> > 
> 
> Perhaps you should tell us the bonded term your force constant in question applies.  I don't see how this has anything to do with a Coulombic potential, as force constants are not involved.  The 1/length dependence (in terms of proportionality, not literally that your energy is calculated as 1/r) is in the units, e.g. kJ/(mol nm^2) - energy is dependent upon the length of the bond, in other words, the displacement from the equilibrium value.
> 
> If your term in question is an angle potential, then the force constant should indeed have units of energy, per the manual (Table 5.5).  If it is anything else, there is an error somewhere.
> 
> > So that means probably, that the writers of the paper did an error. Could
> > such an error ruin my hole MD? My complete force field parametrization is out
> > of this paper.
> > 
> 
> It depends on what the error is, if it exists at all.  If it is a simple typographical mistake, then there's likely no harm.  If there is some larger calculation error, then the force constants may be flawed.  We have no way to know, as you've not said what your bonded term is or what the source of the parameters is.
> 
> The results you obtain in a simulation are only as good as the physical model itself and the assumptions it makes.  If you come to find out that there is some underlying mistake in the parameterization, I would have little or no faith in the results.  Whether you need to be concerned or not at this point is quit unclear.
> 
> -Justin
> 
> -- ========================================
> 
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> ========================================
>