Document code for reactive monte carlo. ping @gitesei

The AtomData structure wires statically all possible properties of the atom type, though they may be used only at a single specific point, e.g., in a coefficient matrix of a pair potential. Moreover addition of a new pair potential often means extension of the AtomType structure.

The atomic properties which are (solely) pair potential coefficients should have a free structure, e.g., the keyword-float pairs. The potentials will then have additional arguments with reasonable defaults, mapping the potential coefficient to the atomic properties, e.g., sigma: "sigma", epsilon: "epsilon". The only danger I can see so far is a conflict, if we decide to use a name chosen by the user, as a keyword with a different purpose in a later version of Faunus.

This approach reduces the AtomData size as well as makes the potentials decoupled from the hardcoded atomic properties. An easy representation could be a dynamically allocated Eigen matrix with cols/rows identified by atom types and rows/cols by used coefficient types. The mapping of the coefficient names to the matrix indices could be stored as a vector or similarly. Refactorization of AtomData and potentials shall be quite straight forward. However, it shall be remembered that AtomData can be access from any place in the code.

Critical discussion welcome.

I tried to simulate a concentrated solution of ions using modified seawater example. To do it I divided the simulation into few steps: 1) equilibration using Weeks-Chandler-Andersen potential (to avoid overlaps of particles) 2) equilibration using Coulomb hard sphere potential 3) production run using Coulomb hard sphere potential. When I try to perform such simulations for small system (cuboid 10 nm, around 600 particles), it works properly. It doesn't work for similar system but 8 times bigger (cuboid 20 nm, around 4800 particles). I have no idea what may cause this issue. Please find attached the bash script used to run simulation and source code of programs to run equilibration and production runs.
seawater_example.zip

Document the use of the following potentials (in topology.md, energy.md):

• squarewell (a8f7ae5). sigma_sw is not a good as it would suggest that it's a diameter but is interpreted as a radius. Perhaps threshold_sw?
• hertz (a00834b)
• dipoledipole (782c284)
• stockmayer (782c284)

Ping @bjornstenqvist

1. The coefficient matrices are implemented like PairMatrix inside the respective potential.
2. PairMatrix is constructed using a Mixer class which takes a vector of AtomType, a lambda function to extract the property (e.g., sigma, epsilon) from the AtomType and an actual combining function (e.g., arithmetic or geometric average). We do not assume the mixing depends on other parameters, i.e., the sigma value can be computed solely from a pair of sigma values.
3. The Potential class can override any default pair parameters. That will be used for custom pair values. The custom pair values are internally treated as AtomType.
4. Serialization to JSON shall check if the parameters in the matrices correspond to the mixing rules; otherwise they shall be output as custom pairs.
5. Serialization from JSON shall be greatly simplified.

Faunus do – beside running MC simulations – also analysis of data. It would be beneficial to perform the same analysis on an existing trajectory without actually running the simulation. It shall be considered that several properties may change during the simulation, e.g., an atomic charge, and not all formats can handle them, e.g., xtc format stores only atomic coordinates, pqr format stores coordinates, charges and radii.

The analysis is performed in parallel to the simulation. However only summary values are usually written at the end. It can be reasonable to write out the analysis result progressively. Either values for every (n-th) step if desirable or as an average over a block of n steps. The output can be written to separate files or in various formats, e.g., csv, hdf5, depending on the nature of data.

The nonbonded interactions have to be efficiently computed when a cut off distance is introduced. A good strategy can be to limit the particle candidates to a content of cells in a cell list that are within the cut off distance.

1. The system can have many cell lists, grouped based on the cut off distances. The cell lists do not interact with each other if the cut off distances differ, i.e., the cell list universe is determined by the cut off distance. The location of a cell in each universe is determined by a single number computed from its k, l, m coordinates.

2. Different strategies can be used to computed nonbonded energies.

3. In the most simple way, a single cell list containing all atoms is constructed for the given cut off distance. The particle candidates are then selected from corresponding neighboured cells.

4. Another approach is preferred in systems of rigid proteins in an implicit solvent / explicit salt environment. Every molecule/atom group holds its own cell list containing only its atoms. To compute the energy, the outer cycle loops over the groups, skipping the unwanted ones (e.g., the same molecule), and the inner loop determines the cells within the cut off distance and compute the energy for the atom pairs.

5. The cell list in the previous system can be implemented as sparse, i.e., a map of cell indices → sets of atoms, or as dense, i.e., a vector over all indices → sets of atoms. The sparse list is typically suitable for big molecules like proteins, while the dense list suits to disperse particles like ions.

6. The cell lists are updated upon the move and are completely rebuild when the volume of the box changes.

7. The current g2g function and its reimplementation according to 4 should be benchmarked.

Depends on #162.

The test fails showing an error message indicating that eps value is not set. The error occurrence depends on the build environment. I am not able to reproduce it in my Linux system using either gcc 8.3 or clang 7.1 and also travis passes without issues.

The error has appeared in Giulio's @gitesei system (MacOS with clang) so far. It stems from the copy assignment json j_remaining_properties = val; in atomdata.cpp. The new json object is not iterated as a json-object but as a json-array instead, i.e., the keys become numeric indices instead of attribute names. I did not find any indication for such behaviour in the nlohmann json library documentation.

Nest step: Learn more about Giulio's system, e.g., clang version.

Travis x64/clang fails in what appears to be related to pybind/pybind11#1604. Hopefully this is fixed in the next pybind release; a patch has been merged with the master branch (pybind/pybind11#1957)

In Nonbonded::g_internal() this uncommented code, which should be equivalent to the active code, causes the speciation test to fail. Why?

Tensors (see particle.cpp) are currently rotated using a rotation matrix. It would be more efficient to change this to quaternions as is done for vectors (now particle rotation functions required both matrix and quaternion). Eigen transformations.

ping @bjornstenqvist

Energy::Hamiltonian is now exposed to python, see

https://github.com/mlund/faunus/blob/master/examples/pythontest.py

At the end of the script, SASA energies are calculated using a pair-potential and PowerSASA, respectively. As it’s just two particles, the two methods should return identical energies, which they do, except at long distances, where PowerSASA misbehaves. The pair-potential uses mixing rules for atomic tfe and tension which will need to be discussed further.
The bug seems to be sporatic; running the test several times, the correct value occasionally appears.

Interesting links about SASA calc.:

• FreeSASA, C library
• C++ code: https://github.com/nioroso-x3/dr_sasa_n
• more code: https://www.cgl.ucsf.edu/chimera/data/sasa-nov2013/sasa.html
• c++ code: https://github.com/atfrank/PCASA
• review of methods

Would it be possible to get a Jupyter notebook (& the environment it was set up in), to calculate charge and second virial coefficient over salt & pH space? Similar to Figure 2 & 3 in your Faraday Discussions paper?

Kurut, Anıl, and Mikael Lund. "Solution electrostatics beyond pH: a coarse grained approach to ion specific interactions between macromolecules." Faraday discussions 160 (2013): 271-278.
 

Currently Gouy-Chapman (externalpotential) has a number of limitations. From the manual:

• this is experimental and subject to change
• does not work with volume fluctuations in the z-direction
• the salt concentration is assumed equal to the ionic strength, i.e. 1:1 salt only
• temperature is hardcoded to 300 K.

I think the way default values are provided may be incorrect. Consider this code:

  inline UnitTest::UnitTest(Tmjson &j) {
    cnt=0;
    auto _j = j["system"]["unittest"];
    file    = _j["testfile"] | string();
    stable  = _j["stable"] | true;  // <-- This always evaluates to true
    include( file );
  }

A better way would be:

    stable  = _j.value("stable", true);

Which works as expected: look up the value, or return the default when it can not be found.

I guess this replacement is only necessary if:

• the value provided by the config could be falsey (0, nullptr, false or "")
• the default value is truthy (ie not falsey)

• stockmayer fluid example
• unittests
• documentation

ping @bjornstenqvist

Pushing to the master branch triggers a documentation build at http://faunus.readthedocs.io.
Some equations do not render since LaTeX _ is interpreted as markdown italic in some equations. This can be fixed using \_ but this would break the pdf build using make manual

• fix equations
• ensure that readthedocs, make manual and make manual_html works
• skip documentation upload to mlund.github.io/faunus? If so, redirect the Docs tab to readthedocs.

This thread is for discussing how energy evaluation may be enhanced using GPU's.

Some ideas:

• STL for GPU: https://developer.nvidia.com/thrust
  • see i.e. their "fancy iterators" which can be used to zip up arrays of data
  • also runs on CPUs using OpenMP or TBB: https://github.com/thrust/thrust/wiki/Device-Backends. Use e.g. -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_X where X=CPP, OMP, CUDA, TBB
  • header only
  • thrust is also used in the Aboria particle library, see also #162
• force biased MC where forces are calculated using GPU code, either internally or via OpenMM

The program produces the following errors while trying to simulate a single ion insertion using the Widom sampling method (input file in .txt attached - to be converted to .json):
in_widom_single.txt

[json.exception.type_error.316] invalid UTF-8 byte at index 25: 0x91

or

faunus(52245,0x11173b5c0) malloc: Incorrect checksum for freed object 0x7fc6d8600048: probably modified after being freed.
Corrupt value: 0x40500db6597a4841
faunus(52245,0x11173b5c0) malloc: *** set a breakpoint in malloc_error_break to debug
Abort trap: 6

or, if it happens to execute to 100%

segmentation fault: 11

Adding a coulomb term to the WCA pair-potential, i.e. - coulomb: {type: plain, epsr: 80}
the sodiumphosphate/ example encounters a particle-particle distance of exactly zero, which triggers an assertion error in CoulombGalore. The problem doesn't happen if starting from the test restart file. I have confirmed that this happens also on a simple GCMC salt system.

ping @gitesei @rc83

The virtual volume for pressure calculation can be unstable. Isabel has a setup that triggers the bug.

To speed up the lookups when energies of non-bonded potentials are computed.

The following assertion fails on clang8/gcc8/gcc7 on xubuntu 19.04:

filter is set in the constructor and an identical assertion passes there. I cannot reproduce this on macOS.
Ping @rc83

Faunus heavily rely on configuration provided by a user in JSON file. However the JSON configuration is not validated at all, hence unknown properties may be silently ignored. This can happen for example by typo in optional parameter or using newer configuration with older version of Faunus.

To mitigate this issue, in some parts of the code SingleUseJSON class used which removes its content, a key by key when it is asked for the given key. This is less than optimal as it breaks const json& semantic when construction of compounded objects is involved.

Using the JSON schema validation can resolve this issue. There is a C++ library, based on JSON for Modern C++ library already used by Faunus, which can be employed. Writing a schema is pretty straightforward. Hence this can be a way to go.

May be I missed something about shift potential or doing something wrong, but when I calculate energy by p2p function and using DebyeHuckelShift potential I get -nan. When I change that to DebyeHuckel potential I get real number. I tested this on 2 residue system as well as large one.
Details:
//typedef CombinedPairPotential<DebyeHuckel,LennardJones> Tpairpot;
typedef CombinedPairPotential<DebyeHuckelShift,LennardJones> Tpairpot;
energy=pot.p2p(spc.p[i],spc.p[j]);
cout << energy;
If needed I can submit all the input files.

Since it's more stable and works for all pair-potentials (incl. anisotropic), the default nonbonded scheme should be FunctorPotential. To keep backwards compatibility, we could define the following:

OK?

ping @gitesei @rc83

The purpose of Space::Change is to stored information about change to be applied to a system. A move can for example specify such a change which makes it trivial to calculate i.e. energies. This issue is for discussing future changes to energy evaluation in faunus.

Vidar would like to have a system where Lennard-Jones parameters of a given element can change as a result of move.

As a quick implementation option could be to expose the init() method in an API of the MixerPairPotentialBase, possibly allowing recalculation only of affected pairs.

A better but more complex solution would be perhaps an observer pattern.

Neither gcc nor clang on macos seem to auto-vectorize code in scatter_test.h despite enabling e.g. -march=native. In contrast, running in a ubuntu docker image on the same hardware, a significant speed boost is observed, presumably due to the availability of libmvec. How can this be improved?

Native, GCC9

Ubuntu docker image, GCC9

Notes:

1. Testing also Clang pragma's and diagnostics described here
2. Manually unrolling loops in the SIMD case improve speed slightly significantly.

Re-starting with state file generated by the following input fails:

./input.yml | faunus -s state.json

temperature: 298
random: {seed: fixed}
geometry: {type: hexagonal, length: 100, radius: 100}
mcloop: {macro: 10, micro: 10}

atomlist:
    - HD:   {sigma: 4.64, eps: 2.252473636363636, dp: 2, alphax: -11}
    - TL1:  {sigma: 6.56, eps: 2.252473636363636, dp: 2, alphax: -20}
    - TL2:  {sigma: 8, eps: 2.252473636363636, dp: 2, alphax: -20}
    - AM:   {sigma: 2.4, eps: 2.252473636363636, dp: 2}
    - HAM:  {sigma: 2.4, eps: 2.252473636363636, dp: 2, q: 1}
    - NA:   {sigma: 4.6, eps: 0.01, dp: 30, q: 1}
    - CL:   {sigma: 4.6, eps: 0.01, dp: 30, q: -1}
    - H:    {pactivity: 4, implicit: True}

moleculelist:
    - Na: {atoms: [NA], atomic: True, activity: 0.112}
    - Cl: {atoms: [CL], atomic: True, activity: 0.112}
    - lipid:
        structure:
            - AM:  [0,0,0]
            - HD:  [0,0,3.5199999999999996]
            - TL1: [0,0,9.12]
            - TL2: [0,0,16.4]
            
        bondlist:
            - fene: {index: [0,1], k: 1.0558470170454544, rmax: 5.279999999999999,
                    sigma: 3.5199999999999996, eps: 2.252473636363636}
            - fene: {index: [1,2], k: 1.0558470170454544, rmax: 8.399999999999999,
                    sigma: 5.6, eps: 2.252473636363636}
            - fene: {index: [2,3], k: 1.0558470170454544, rmax: 10.919999999999998,
                    sigma: 7.279999999999999, eps: 2.252473636363636}

            - harmonic: {index: [0,2], k: 0.35194900568181814, req: 18.159999999999997}
            - harmonic: {index: [1,3], k: 0.35194900568181814, req: 25.64}
            
insertmolecules:
    - Na: {N: 200, inactive: True}
    - Cl: {N: 200, inactive: True}
    - lipid: {N: 160}

energy:
    - customexternal: 
        molecules: [lipid]
        com: false
        constants: {radius: 12}
        function:
            var r2 := x^2 + y^2;
            if ( r2 < radius^2 )
               1e10;
            else
               0;
    - bonded: {}
    - nonbonded:
        cutoff_g2g: 50
        default:
            - wca: {mixing: LB}
            
        TL1 TL1:
            - wca:
                mixing: LB
            - cos2: {rc: 7.36335072, eps: 2.252473636363636, wc: 11.152}
        TL1 TL2:
            - wca:
                mixing: LB
            - cos2: {rc: 8.17152336, eps: 2.252473636363636, wc: 12.376}
        TL2 TL2:
            - wca:
                mixing: LB
            - cos2: {rc: 8.979696, eps: 2.252473636363636, wc: 13.6}
        NA HD:
            - wca:
                mixing: LB
            - polar: {epsr: 78.47995215311005}
        NA TL1:
            - wca:
                mixing: LB
            - polar: {epsr: 78.47995215311005}
        NA TL2:
            - wca:
                mixing: LB
            - polar: {epsr: 78.47995215311005}
        CL HD:
            - wca:
                mixing: LB
            - polar: {epsr: 78.47995215311005}
        CL TL1:
            - wca:
                mixing: LB
            - polar: {epsr: 78.47995215311005}
        CL TL2:
            - wca:
                mixing: LB
            - polar: {epsr: 78.47995215311005}

moves:
    - moltransrot: {molecule: lipid, dp: 0.5, dprot: 0.5, repeat: 10}
        
analysis:
    - systemenergy: {file: energy_0.dat, nstep: 10}
    - savestate: {file: state.json}
    - savestate: {file: confout.pqr}
    - sanity: {nstep: 1}

Currently the only way to skip interactions is to use a group-to-group based cutoff for molecular groups. This is inefficient for small molecules, of if atomic groups are in use. The idea is to implement a cell list where particles and particle mass centra are assigned to a grid point. A basic structure for this already exists.

Things to consider:

• who holds the cell list (space, nonbonded?)
• periodic boundaries or not
• updating (particle movement, volume scaling)

The Aboria library is a particle container, i.e. corresponding to our std::vector<Particle> and implements for neighbour list routines. It requires boost, though.

This issue is for discussing ways to ease input generation for Faunus using i.e. IPython with tab completion. For example:

from faunus import config

c = config(indent=4, type='json')
c.clear()

c.temperature=298.15;
c.geometry.cuboid( length=100  )

c.analysis.widommolecule(nstep=20, ninsert=100, group="water")
c.analysis.widommolecule.clear()

c.moves.atomtranslate.add(group="salt", peratom=True)

c.atom.add(name="Na", charge=1.0, Ninit=10)
c.molecule.add(name="water", Ninit=10, structure="spce.pqr")

print(c) # --> dict
c.save('input.json')

Possible ideas:

• http://stackoverflow.com/questions/13870241/ipython-tab-completion-for-custom-dict-class
• https://pypi.python.org/pypi/schema
• https://pypi.python.org/pypi/voluptuous
• https://www.blog.pythonlibrary.org/2014/02/14/python-101-how-to-change-a-dict-into-a-class/

The Hamiltonian of many force fields does not include a non-bonded energy between particles which are connected by a bond (1-2 interaction). Furthermore, this can be also extended to 1-3 and 1-4 interactions where the mutual interaction energy is already fully determined by the angular and dihedral potentials, respectively.

Therefore the specification of particle pairs with excluded non-bonded energy should be permitted in the topology. Only particle pairs within a single molecule are allowed. Beside an explicit configuration option, convenient shortcuts can be also provided, e.g., an extra attribute of the bonded interaction.

Subtraction of the non-bonded energy is perhaps a good implementation strategy, accounting for the performance. However a special attention has to be paid for high or even infinite energy values, e.g., overlapping hard spheres (the sphere diameter is greater than the bond length). Van der Waals interaction and electrostatics should be handled independently internally.

The implementation can also account for a 1-4 pair interaction, i.e., a bonded potential in a form of the scaled non-bonded potential.

Create a framework for scattering form factors, F(q), for particles. This could be atoms, coarse grained amino acids etc. The idea is to inject into the existing scattering routines.

when pqr and xtc files are printed after move they do not correspond to each other. When first particle moves, in xtc it is movement of second particle.

Todo list

• Get swap move working
• Merge PR #102 from @bearnix or new one from @gitesei
• Decide on the input format
• Create tests
• Write documentation
• Does the atom / molecule tracker give a speed gain? If not, it is simpler to call existing Space
• rename pactivity --> pX? In to_json(), write only activity.
  functions

Here is my initial suggestion for how the user input could look like. Currently activities for atoms are given as molecule (with atomic==True) properties, but I think it's more intuitive to specify it for the atomic species. Any thoughts on this?

implicit protons

atomlist:
    H+: { q: 1.0, activity: 1e-7, implicit: True}
    HA: { q: 0.0, ... }
    A-: { q: -1.0, ... }
moleculelist:
    protein: {atoms: [HA, HA, ...]}
reactionlist:
    "A- + H+ = HA": {pK: 4.75}

Analysing the reaction, we automatically detect that H+ is implicit and identify the reaction as a swap. pH is taken from the activity in the atomlist. For implicit titration, there's in this way no need to specify a fictitious molecule for protons.

implicit protons and GC salt bath

atomlist:
    H+: { q: 1.0, activity: 1e-7, implicit: True} # alternatively `pactivity: 7.0`
    HA: { q: 0.0, ... }
    A-: { q: -1.0, ... }
    Na+: {q: 1.0, ..., activity: 0.01 }
    Cl-: {q: -1.0, ..., activity: 0.01 }
moleculelist:
    protein: {atoms: [HA, HA, ...]}
    Na: {atomic: True, atoms: [Na+]} # activity picked from atomlist
    Cl: {atomic: True, atoms: [Cl-]}
reactionlist:
    "A- + H+ = HA + Cl": {pK: 4.75}
    "= Na + Cl": {} # bearnix, is this the way to set GC for salt? Need policy for picking either activity OR given `pK`

we still identify H+ as implicit and insert/delete chloride. The ReactionData figures out what is atomic and what is molecular. If the molecular group has atomic==True, then activities are taken from the atomlist for the correspondingly contained atom.

Is this correct?

ping: @bjornstenqvist

Add analysis that tracks the average charge product between titratable sites. The output could look like this:

i j <z_iz_j> <z_i> <z_j>

Below is a benchmark of the examples/bulk system of melted NaCl for the same number particles, trial moves, release mode; and all analysis disabled. Hardcoded, non-splined potentials are used in both cases.

1. faunus-v1.0.2: 54 secs
2. faunus-master (f5e11ef): 105 secs

Possible places to look:

• geometry::vdist()
• passing of change
• compare with nonbonded switched off, i.e. probe only trial move code
• compare clang and gcc

ping @rc83

Sometimes the simulation crashes. It would be useful to restart the simulation just before the crash to quickly reproduce the conditions for debugging. Hence the simulation state shall be regularly output into a checkpoint file in the same way as it is written at the end of the simulation as a part of the analysis (SaveState).

Cmd line arguments: Number of macro MC steps between the checkpoint file updates.

Add non-equilibrium candidate MC move (NCMC).

• basic C++ framework, see i.e. http://github.com/choderalab/openmmtools
• bistable dimer example from paper

A lightweight logging framework shall be incorporated into faunus to enable efficient logging at various levels (e.g., error, warning, notice, debug). The logging framework shall be also capable to choose different sinks like stderr and/or file. The log level and the sink shall be command line options. Currently only messages to stderr are output.

Some examples of use cases: 1. When adding a new feature you realize that some combinations of parameters will not be currently used and hence tested thus a notice or warning shall be emitted. 2. An obsolete syntax in json configuration is encountered. 3. Very improbable state is detected.

Next step: Review logging framework candidates.

The reaction coordinate for molecule-molecule separation is implemented in several ways and documentation should be updated. src/penalty.h

• Add reciprocal Ewald energy also for Functor and splined potentials, and for dipoles.
• Is Ewald energy added twice in current scheme? The selfEnergy() in CoulombGalore is added automatically, and so is the self-energy from the IonIonPolicy
• Have tests for self-energies (grand canonical Ewald, qpot?)
• Add self energies for CoulombGalore and DipoleDipoleGalore
• Complete Ewald documentation, incl. the use of Yukawa Ewald.

ping @bjornstenqvist

A benchmarking system to compare performance between commits and versions would be useful. Catch2 already has support for this and there's an open issue at doctest:

doctest/doctest#251

Here's a simple solution for doctest: https://github.com/martinus/nanobench/blob/master/docs/tutorial.md

• document atomdipdipcorr in code and in manual
• document slicedir and thickness
• is PairAngleFunctionBase really needed? What it injects could be done directly in atomdipdip which would increase clarity.

ping @bjornstenqvist

As the Hertz potential was broken it would be useful to cover it with unit tests. Björn @bjornstenqvist, can you please create a test case in potentials_test.h? The LennardJones test case can be easily taken as a template. Many thanks!

I noticed a commit that was used in a faunus-notebook, but was "orphaned", i.e. it was no longer part of the history of any tag or branch. The result is that this commit cannot be fetched by anyone who doesn't already have it.

I think this wasn't noticed before because the offending git checkout was one line of a bash script. If a command fails, bash silently ignores the error and moves to the next line. (See my comment on the faunus-notebooks repository)

The solution would be to create a tag for a commit if it is used in a notebook. Tags are very lightweight: it is a 40 bytes file which only content is a commit hash.

Pressure calc. using "virtualvolume" and "virial" give different results in bulk.cpp example if CoulombGalore is used.