Re: editing the code of the non bonded pair potentials

From: Joshua Adelman (jla65_at_pitt.edu)
Date: Fri Aug 30 2013 - 14:06:12 CDT

On Aug 30, 2013, at 2:47 PM, Kenno Vanommeslaeghe wrote:

> On 08/30/2013 10:10 AM, Axel Kohlmeyer wrote:
>> On Fri, Aug 30, 2013 at 12:26 PM, Asaf Farhi <asaf.farhi_at_weizmann.ac.il> wrote:
>>> If (E_VDW>E_c)
>>> {
>>> E_VDW=E_c;
>>> f_VDW=0;
>>> }
>>> If (E_el>E_c)
>>> {
>>> E_el=E_c;
>>> f_el=0;
>>> }
>
> Asaf, it looks like you're introducing a discontinuity in the first derivative. While this is OK for MC, it violates the underlying assumptions that make the Verlet integrators used in MD well-behaved. You'll be needing a switching function.
>
>> (...) then you probably should modify an MD
>> code that is simpler to understand and modify. there are many of those
>> around.
>
> Axel, sorry for going a bit on a tangent here, but purely out of educational interest, can you give a pointer to one (or more) from the top of your head? A quick Google search turned up toyMD, but since you appear to be more knowledgeable in these things than I am, I'd like to know your opinion.
>
> Cheers,
>
> Kenno.
>

You might want to take a look at OpenMM (https://simtk.org/home/openmm) if you want to mess around with how forces are being calculated. They have a number of nice facilities for specifying custom forces, integrators, etc and their python interface makes it simple to prototype simulation protocols. I don't know about Asaf's particular use-case, but here is an example of how to code up aMD in OpenMM for example:

https://simtk.org/svn/pyopenmm/trunk/simtk/pyopenmm/extras/amd.py

Cheers,
Josh

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:35 CST