> On Mar 9, 2015, at 5:11 PM, Sanmeet Chahal <schah063.uottawa.ca> wrote:
>
> Hello all,
>
> Sorry to bring up an old issue, but I was reading the amber14 manual and on
> page 245 found the following example:
>
> # Make atoms 40, 41, 42, and 57 a new LJ type with 0s for
>
> # their parameters to remove all of their LJ interactions
>
> # with every other atom
>
> addLJType .40-42,57 0.0 0.0
>
>
> I'm confused because it seems that having the epsilon as 0.0 was what was
> causing the segfault when I ran the run minimization step of Tutorial B0
> because after changing the epsilon to be a non-zero value the error stopped
> occurring. Could some one please explain when it is ok to set LJ parameters
> to 0 and when they could cause an error?
This was simply meant to illustrate an example of what ParmEd can do with addLJType. It was not meant as an example that shows what you *should* do. The most common example I can think of where you would *want* to zero-out Lennard-Jones interactions is if you were doing some kind of alchemical transformation like thermodynamic integration or free energy perturbation (TI or FEP). Imagine trying to find the binding free energy difference between phenol and benzene. In that case, the phenol has a hydroxyl group where the benzene has a hydrogen. One transformation you can do is to turn the phenol O into an H and the hydroxyl H into a dummy atom with no charge or Lennard-Jones interactions (i.e., so it doesn’t interact with *anything*).
In this way, ParmEd lets you create a series of topology files that interpolate between two end-states of an alchemical transformation much more easily than if you tried to do this whole thing through tleap. While this process is not needed for TI as it is implemented in Amber, if you wanted to use the replica exchange free energy perturbation method (through Hamiltonian replica exchange), you *do* need to use this approach.
The reason that your system was seeing problems when you set epsilon = 0 is that your hydrogen had no Lennard-Jones interactions *AND* it had a charge. So there are pathological cases where an atom could collapse onto your hydrogen and the force due to the charges becomes VERY large (since the force varies as 1/r^2). If the distance between those atoms are 0, the force becomes infinite, and your system blows up. Or if it’s very small, it’s not *infinite*, but is still so large that it behaves effectively the same. [1]
In the example I gave above with alchemical transformations, you would change both the Lennard-Jones *and* the charge of the dummy atom to 0, so no matter how close that atom got to any other atom, their electrostatic and van der Waals forces would both be zero, thereby avoiding the force catastrophe that destroyed your previous simulation. The likelihood of two noninteracting particles becoming *completely* coincident (so that r between them is exactly equal to 0) is so small as to be insignificant.
HTH,
Jason
[1] Normally the Lennard-Jones interactions of the *oxygen* is enough to keep atoms from collapsing onto the hydroxyl hydrogen atom. However, the atoms that the oxygen forms a bond *or* an angle with are excluded from the nonbonded interaction with that oxygen atom. So an atom that forms an angle with the oxygen -- which forms a dihedral with the hydrogen -- may collapse into the hydrogen since the oxygen’s Lennard-Jones interaction is not computed with that atom and the hydrogen *has* no Lennard-Jones radius. The angle force constant may not strong enough to prevent this collapse
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Mar 09 2015 - 19:30:02 PDT
