[gmx-users] Potential energy methanol in a box

[Justin A. Lemkul]

Fabian Casteblanco wrote:
> Hello,
> 
> I'm trying to fill a box with methanol using CHARMM FF Parameters.  I
> also need to do this for ethanol and 1-propanol.
> 
> For ethanol, each individual molecule had approximately -19 kJ/mol of
> potential energy, then I placed 1000 in a box, performed nvt, npt,
> etc.  End Result:  A negative potential energy (it means that its
> stable?)  and a density of 0.785 g/cm3 compared to actual 0.789 g/cm3
> at 1 atm, 298 K.  This result seemed ok.
> 

Negative potential implies net attraction.  Positive potential implies net 
repulsion.  Systems in the condensed phase typically have negative potential 
energies, but you should not necessarily equate negative energy with stability. 
  Unstable systems can have negative potential, too.  But if your trajectory 
itself is stable, then likely so too is the system you're looking at.

> For methanol and  1-propanol, each individual molecule had a larger
> 24 - 64 kJ/mol respectively, due to Coulomb 1,4 pair interactions (for
> example, for methanol between the HA 0.09 and the H 0.43 gives a large
> 64 kJ/mol).  Again placed 1000 each in 2 separate boxes, performed
> nvt,npt,etc on each one.  End Result:  A still positive potential
> energy for each box (shouldn't this be negative after running npt? on
> both) and a density that was 0.789 g/cm3 (for 1-propanol)compared to
> actual 0.803 g/cm3 1-propanol at 1 atm, 298 K.  Methanol density a bit
> closer 0.787 compared to 0.791 g/cm3.

Negative potential does not simply arise from NPT simulation.  The potential is 
determined by the topology and functional form used.  Something about these 
systems is causing a net repulsion, and it may very well be the intramolecular 
interactions.  That seems likely, since the densities you're getting are 
slightly low, indicating that the molecules are insufficiently condensed.

> 
> My questions are:
> -There should definetly be a positive potential energy after running
> npt, correct? Otherwise it wouldn't be stable?  I understand that the

No.  See above.  Ensemble does not dictate whether or not the potential is 
positive or negative.  The topology does.

> 1,4 coulomb interactions should cancel out with nearby molecules.

No.  1-4 interactions are intramolecular.  They are not dependent upon the 
presence or absence of nearby molecules.  They are dependent upon the 
configurations of the individual molecules themselves.

> *For CHARMM, they state that unless otherwise 1,4 parameters exist
> already, they use a full scale 1.0 of the regular LJ, coulomb
> potential.
> 

Correct.

> -As for the density, I know it's somewhat close but I'm trying to be
> as accurate as possible.  I read somewhere that the average cut-off
> (rlist) is usually 1.4 for CHARMM and I had mine set to 1.  I'm
> guessing that having a larger cutoff would contribute more molecules
> to the LJ, coulomb which might rise the density a bit to the correct
> value?
> 

The value of the cutoff should be no less than 1.2-1.4 nm, according to the 
CHARMM developers (http://dx.doi.org/10.1002/jcc.21287).  It's hard to predict 
the effects of such changes, since they will affect neighbor searching, balance 
of PME terms, and others.  You should certainly be using settings recommended 
(and used) by the force field developers, unless you have a substantive reason 
not to (i.e., someone has developed better settings).

-Justin

> I'm still new to Gromacs so I appreciate all of anybodys help.  Thanks a bunch.
> 
> --
> Best regards,
> 
> Fabian F. Casteblanco
> Rutgers University --
> Chemical Engineering PhD Student
> C: +908 917 0723
> E:  fabian.casteblanco at gmail.com
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================