nqcd / nqcdynamics.jl Goto Github PK

We discussed that it would be nice to also showcase that the code in principle can perform on the fly dynamics on realistic systems.
One way to do that is to take one of the H2 on Ag(111) trajectories that @oadoug has simulated with the code and run on the fly dynamics with NQCD and FHIaims via the ASE interface.
This means the NQCD code drives the dynamics using an ASE calculator for FHIaims.
That calculator has a get_forces() and get_friction() function.
It's enough to do this for one set of initial conditions to show that it is principally possible.

The comparison of the evaluated probabilities to the random number is not quite right.
Should be doing a < comparison rather than a >.

See here for correct version:
https://github.com/smparker/mudslide/blob/9fffecf749213936f913237b2734046462c16d7c/mudslide/trajectory_sh.py#L421

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/8dcd94f85b7dc6ee5e3530c90fb082cc4de258fc/docs/src/nonadiabaticmodels/overview.md#L35

Are nstates and the degrees of freedom hard-coded into the model? If not, I would specify them as arguments in the example, so it's clearer how to change them. If they are, it should be mentioned that they are hard-coded.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/6039a6c2e0c9be83868f8f0fc9a04d6fd7e978f6/docs/src/dynamicssimulations/dynamicsmethods/classical.md#L23

In the 'Getting Started' section, it's written that simulation is set up as Simulation{...}, however, this does not seem to be the case here. Is that a mistake?

If the models are publicly hosted we can use BinaryBuilder.jl for compiled models and Artifacts for data files. This will make it very easy to use ML models.

Could also consider moving the EANN type models to a separate repository. This can be an add-on to the base package that provides extra models that we do not want to release publicly.

I think a possible option is to define a type that contains a vector of DynamicalVariables instances, where each instance is a separate bead of the system. This type can also contain various ring polymer quantities like the normal mode transformation and the harmonic spring matrix.

Functions can then be defined to dispatch on either DynamicalVariables or Replicas{DynamicalVariables}, which would usually involve repeating the propagation for each copy of the system.

This allows for a separation of the typical usage with classical nuclei and the ring polymer alternative that should make it simpler for people who are not interested in ring polymer methods.

E.g

User:
cutoff_constraint = [H2 bond length > 4 + H2 molecule height > 3 or time > 400]

dynamics.attach(cutoff_constraint)

maybe?

This should be implemented for MDEF and for thermal Langevin dynamics. Also worth considering if combining the ring polymer algorithm with the classical algorithm is a good idea.

@oadoug has encountered an issue where the model does not work with the parallel ensembles. Both the EnsembleThreads and EnsembleDistributed algorithms are broken.

The EnsembleThreads issue makes sense since this is shared memory and variables are allocated many times. This is likely not possible to fix.

EnsembleDistributed should work, since it is distributed memory but gets an access to undefined reference error on worker 2. This happens on line 62 of EANN_H2Ag.jl where the ccall derivative is performed. Not sure how to solve this?

This will largely involve porting over existing code from NonAdiabaticRingPolymers.jl.

It would be good to use unitful units so it is clear what is going on but is a bit annoying to have to deal with the units all the time. What is the best way to do this?

I forgot to do this.

In the documentation there should be a section Dynamics that will have a separate page for each of the types of dynamics included so far.

Each of these dynamics types should have a description of the theory behind the method with references included. Then an example of how to run a single trajectory and plot some observable from that.

Any quirks or specifics of the implementation should also be mentioned here.

Methods to include

• NRPMD
• Ehrenfest (also ring polymer version) @wgst
• FSSH/RPSH (should these be separate? probably not)
• Classical dynamics
• Thermal Langevin dynamics
• MDEF @oadoug

Now we have two versions of surface hopping (three including RPSH), it would be good to move these into a sub directory and ensure code repetition is kept to a minimum. It should be possible to have them all follow a basic structure, then have specialised methods for each dynamics type. This will also make it easier to add more surface hopping methods in the future.

Anything accessible by the user and possibly any other non-trivial functions should have docstrings added, and these should be setup to show in a section of the documentation.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/9f128283e060b6d7e1638df1c6c0d441266e0fea/docs/src/nonadiabaticmodels/overview.md#L61

Much more accessible now! However, how does it matter to people reading this particular page that the models have different abstract types? Does each abstract type come with different ways of passing in the variables? If so, this should be written here (just one sentence saying how the abstract types differ, just repeat/extend the question in case that's right). What is the API-section and can we give a link to there?
I like the link to How-to-implement the model, that is very useful here.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/6039a6c2e0c9be83868f8f0fc9a04d6fd7e978f6/docs/src/getting_started.md#L176

This is unclear. How do you access this table and this u?

Ideally, the user would provide their system, select the type of dynamics to perform and press go. We are able to hide a lot of the specifics so this is as simple as possible.

Previously, in NonAdiabaticRingPolymers.jl, I did this using the DifferentialEquations.EnsembleProblem interface. I had a function that ran population dynamics and it automatically picked a good solver algorithm and provided the averaged population output. We have the option to do something similar where we can write functions for specific outputs. Or keep things more general.

I never provided a simple interface for running single trajectories, only ever ensembles, as I didn't ever need to run just a single trajectory. However, we could also include simplified functions to do this if this would be a desired use case.

This will make it easier to see the available outputs and removes ambiguity surrounding the units of the results.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/685eb31fd86a54c6f21ef75f4e61bdd389960dd4/docs/src/nonadiabaticmodels/ldfa.md#L17

For which systems is the friction coefficient evaluated and under what circumstances is it evaluated (i.e., where does the electron density come from)? Are these numbers you just read in that are tabulated in the paper or is this a routine that you need to call first? If you need to call it first, please provide an example of how to do that.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/9f128283e060b6d7e1638df1c6c0d441266e0fea/docs/src/nonadiabaticmodels/overview.md#L58

Good example and explanation!
It would be good to provide another example for higher dimensions here, however brief (just because going from the 1D to 3D case is usually not trivial).

Here a proposal for a feature enhancement:

My understanding is that @oadoug and @jamesgardner1421 have already implemented a 2TM-MDEF dynamics approach, so given a T_el(t) profile, MDEF dynamics can be performed where the temperature of the electronic bath is coupled to the 2TM temperature profile. This takes a precalculated T_el profile.

It would be interesting to extend this approach to also couple to the T_phonon profile to the MD that is calculated in the 2 temperature model. See approaches by Peter Saalfrank and Maite Alducin:
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.95.125439

The way how this works is that the atoms of the surface are coupled to a standard thermostat (Nose or Langevin or whatever) and the temperature of that thermostat is coupled to the time-dependent T_phonon profile from the 2TM calculation.

Side note:
@oadoug mentioned at some point that he was interested to port his (currently Python-based) 2TM equation implementation to JL. When using precalculated T_el profiles, that isn't really a high priority right now, but it could be interesting if, in the future, we might be interested to look into some self-consistent approach where the atomistic dynamics are propagated simultaneously with the 2TM and they interact with each other.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/8dcd94f85b7dc6ee5e3530c90fb082cc4de258fc/docs/src/api/nonadiabaticmodels/diabaticmodels.md#L5

Where did the text go for the Diabatic models? I can't seem to get to it when clicking on 'edit in GitHub' on the documentation page.

In any case, long-term, we need to add all the equations and parameters we implement in this model in the description. The papers can be behind a paywall, so not everyone might be able to see the equations.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/685eb31fd86a54c6f21ef75f4e61bdd389960dd4/docs/src/nonadiabaticmodels/analyticmodels.md#L21

This plot and a few others do not have access labels. Could you add them, please? (I would also provide a nice example of how to do that.)

For the Tully-Models and a few others, plot legends are missing and would be helpful, both in terms of giving them as an example and understanding what the models show.

This should be a direct copy of the old code.

The position vector must be converted to the correct units (bohr->Angstrom)

The forces and energies must be converted back into atomic units.

These can be done using Unitful and UnitfulAtomic.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/3132608e9ee40331f6e72ea578abc4ef8887cd0e/docs/src/getting_started.md#L55

Cool box!
What are 'keyword constructors' and what is 'nice printing'? (i.e., is nice printing just formating or is it some kind of explanation?
Perhaps write what keyword constructors do and than put (keyword constructors) behind that. For nice printing, is there a way to describe more clearly what that means?

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/685eb31fd86a54c6f21ef75f4e61bdd389960dd4/docs/src/nonadiabaticmodels/ldfa.md#L7

The description of MDEF and the propagation by Langevin equations needs a citation. Possibly Allen and Tildesley (Computer Simulation of Liquids, Oxford), but you may want to ask Connor if he has any other suggestions.

The Langevin equation should be written out for the general case, so that we can then talk about LDFA and how it relates to the general case.

In NonAdiabaticRingPolymers.jl I've created the NuclearDistribution type which just stores a vector of all the positions with some functions defined to pick a configuration at random and perform simple statistical analysis (mean, var).

Is this still the way to go or is there a better alternative?

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/685eb31fd86a54c6f21ef75f4e61bdd389960dd4/docs/src/nonadiabaticmodels/ldfa.md#L24

An example is missing that shows how you set up this type of calculation.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/6039a6c2e0c9be83868f8f0fc9a04d6fd7e978f6/docs/src/getting_started.md#L153

Eventually, a link to a page about explaining how to run the dynamics with how to set different algorithms would be very helpful here.

The name System for the overall parameter type is not very appropriate. The type doesn't really represent the system, that name is perhaps better reserved for the atomic parameters type which currently has a pretty rubbish name. A better name could be Simulation.

Then we would have Simulation{MDEF} or Simulation{MonteCarlo} and so on.

Then we can have general functions that perform the simulations. Such as: run_trajectory(), run_trajectory_ensemble().

There are a couple of things that could be done to improve performance of MDEF with tensorial and diagonal friction:

• Ensure friction tensors are only 3*n_friction_atoms rather than 3*n_atoms in dimension. Currently there are lots of zeros which slows the eigen decomposition and matrix multiplications.
• Make friction tensors Symmetric type, this will simplify representation of data.
• Switch LDFA algorithm to expect only a vector of friction, rather than a diagonal matrix.

By default each atom will have three degrees of freedom. However, for the model systems, these typically have only one per atom.
Some of the containers and force evaluations, especially in the ring polymer case do not properly account for this. But this will be most easily solved after completing the schnetpack interface.

Similar to #61 we want to make setup as automated as possible. Currently PyCall is required for using ase.io. But loading PyCall also loads extra code like the LDFA model. The problem here is that the LDFA model has extra dependencies (numpy for example). Ideally we don't want to have to install these things if we're not using them.

One solution could be to split weird models into separate packages that can be optionally loaded. Alternatively we can look at automated setup of python.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/685eb31fd86a54c6f21ef75f4e61bdd389960dd4/docs/src/nonadiabaticmodels/analyticmodels.md#L5

I really like this page.
What does this tip mean for anyone using the code? Do they have to load them in with the `using' statement? If so, can you put that into the examples, so it's really blindingly obvious?

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/6039a6c2e0c9be83868f8f0fc9a04d6fd7e978f6/docs/src/getting_started.md#L180

This section needs examples.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/9f128283e060b6d7e1638df1c6c0d441266e0fea/docs/src/nonadiabaticmodels/overview.md#L75

I really like the link to the different models here, that's super useful. The figure is also super helpful.
What is the sidebar you're talking about in the sentence and where exactly are the Extra models (do they fit somehow into the figure?)

The biggest challenge is likely reading cube files. Could use the existing python code or rewrite in julia for possible speed increase. Shouldn't take too long.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/6039a6c2e0c9be83868f8f0fc9a04d6fd7e978f6/docs/src/getting_started.md#L66

The hcat box is really useful!
Before explaining about hcat, please explain how to set up R in general, so by an example for a 1D and 3D case. At the moment, I would not know how to define positions for my system. Could also be a link to another page.

https://github.com/NQCD/NonadiabaticMolecularDynamics.jl/blob/3132608e9ee40331f6e72ea578abc4ef8887cd0e/docs/src/getting_started.md#L11

Very good basic description!
It still needs a bit extension, mainly:

• how can you define a large array of atomic positions, e.g. for 56 Au atoms? I remember from working with the code that it was a bit more complicated than just passing the atom-label into an array and writing it all by hand would be very cumbersome.
• how could you read in a file with already existing geometries? (This could also reference to another page)

It's possible that performance improvements could be gained by using static arrays instead in the calculator/model part. Will have to see if models providing static arrays is faster than the current setup.

Create Models for the potentials in the ML-model-repository.

These should look similar to the existing schnet model.

Requires.jl allows you to load parts of the code only when certain packages are loaded first. This means we can remove PyCall, Conda, JuLIP... Lots and lots. There's a few others too.

It's probably a good idea to port a lot of the initialisation that currently sits in RingPolymers over to this package together with use cases, examples and tests

Provide a flexible, composable interface that allows the user the provide or select from a set of problem, reduction and output functions.
This will be heavily interfaced with the DifferentialEquations.jl ensembles, acting as a light wrapper that streamlines the process.