CNDO

Miscellaneous

A simple molecular dynamics program and some thinkings about free energy calculations

This program was designed to investigate the state-of-the-art approaches for free energy calculations, based on the Lennard-Jones interatomic interaction model. The program is supposed to do the following things:

Implement thermodynamical integration (TI) free energy calculation for a single atom mutation in a Lennard-Jones (LJ) fluid;
Implement the adiabatic switching (AS) and reversible scaling (RS) technique;
Implement the competitive lambda dynamics approach for screening the best bound atom to certain thermal environment;
Lattice dynamics and phonon spectrum calculation for obtaining the reference free energy near absolute zero temperature;

By playing with this program, one can understand some basic concepts in free energy calculations using the MD method which are otherwise buried in the complexity of protein structures and the force fields. Special effort is made on checking the effectiveness of the competitive lambda dynamics for such a simple model. Due to time limitation, I only did free energy calculations for the formation of an impurity in a pure system. My conclusions based on these benchmarking calculations can be summarized as follows (Of course I could be wrong):

With an appropriate switching function and switch time, the results given by the RS method are in good agreement with those by the TI method, however, in the high temperature regime above the melting point, disagreement becomes more and more obvious. This could be due to the fact that the system enters the superheated regime and the free energy change through a superheating route is different from the free energy difference between thermodynamically equilibrated states. Result: @See

Thermodynamical properties of a LJ fluid
Different switching functions
Comparison of free energy results using different methods
I believe there are some theoretical issues with the lambda dynamics approach implemented in Charmm V.27 that need to be cleared. I remain skeptical that that version of the lambda dynamics (LD) had been tested enough. The major question I came up with the simulations for a simple LJ fluid is that the introduction of the Langrange multipler would eventually kill all the other competitors and lead to the result that only the champion can survive after a very long time simulation. Consider the fact that nothing in the equations of motion could stop the (two) finalists from competing with each other if we run the simulation for longer time. No matter how weak the energy difference could be, it will sort the rivals out eventually, given adequately long time, just like fluids will flow as long as there exists a net chemical gradient. Therefore, theoretically speaking, the equation of motion in the LD might not be suitable for the computation of free energy differences---remember the Gibbs postulate about the ergodicity hypothesis, if the equation of motion leads clearly to a final state, it is non-ergodic and therefore any free energy calculation sampling based on this nonequilibrium time ensemble becomes problematical. Conservatively speaking, the LD method serves well in finding the best ligand, but it is doubtful that it can be used to figure out the free energy differences in its present form. Removal of the multipler (i.e. the second-order constraint) seems to put them back, however, it is not preferrable as the holonomic constraint usually requires the Langrange factor. This question seems to have put me on hold. As the energy landscape of a LJ system is much more ruggedless than that of a complicated protein-ligand setting, the system should in principle be more mobile in the phase space. This could be the reason that we see all the other components die off after a long time. In a protein-ligand simulation setting, as the energy landscape is far more rugged, it might not be easy for the ligands to die off in the lambda dynamics. This could have two completely different implications: One is that the lambda dynamics is fine as it might be the case in reality; The other is that the lambda dynamics is misleading before the system has explored enough of the phase space (as the method claims that it is able to rapidly screen the most favorable ligand out of a bunch, this is of particular importance).
The LD method, with appropriate modifications, can be applied to a number of problems in materials simulations, for instance, surface chemisorption, solid solutions, and phase stabilities of alloys. In surface physics, what kind of surface can bind what kind of gas molecules is a hot topic, as the problem is believed important for example to car depollution research. By constructing hybrid gas molecules and applying the LD method, we can rapidly screen which type of molecule will bind to a given (relaxed and reconstructed) metallic surface the most easily. In composite materials research, it is also important to investigate the phase stabilities of alloys. For example, suppose we want to have a B2 alloy which must be composed of metal Cu, and we then need to know which metal can form a B2 structure with a given lattice constant with Cu. The LD approach stands for a novel method of doing this---We build a hybrid atom which consists of different components of Au, Bi, Fe and so on, put this hybrid atom on one of the two vertices in the B2 structure, and the LD will tell us which metal is prefered. This work has been partly done: With a Lennard-Jones potential model, the LD method is clearly capable of picking up the most stable phase.

HP-UX version and Linux version are both available. A short (but early and incomplete) description of the study on the lambda dynamics method can be found in a postscript file.

A tight-binding electron transfer program

A simple playground program for electron transfer dynamics on the tight-binding model can be found here. With this program, the user can play simple tight-binding model, do two-state and three-state reductions, and check the two-state and three-state approximations for simple models, and perhaps get some intuition out of them. Here are some postscript figures: (1) Electron dynamics with different donor-accept distances; (2) Electron dynamics when the bridge frequency is different.

I would like to point out here that the multi-state reduction may be an interesting problem. But we have not been able to explore more with it. It is mentionworthy that the multi-state resonance might be nontrivial in some circumstance. Multiple state resonances are usually more difficult to find, dynamical multiple state resonances are even harder to detect. But can we rule out the possibility of multiple state resonances and their contributions to extraordinary electron transfer? For details, read a short note.

Demo applets for electron transfer dynamics in a tight-binding model: Static bridge and Dynamic bridge may be interesting for a beginner in electron transfer theory.

The Floquet program

We have spent efforts on exploring the applicability of the Floquet formulism to the electron transfer problem, and eventually came up with the conclusion after a thorough discussion with Professor David Beratan that it might not be relevant to the biological electron transfer problem. However, the Floquet picture is still helpful in understanding electron propogation behavior through a dynamical bridge. Hence it might be useful to present the program we have written here.

We present three interesting results which were produced using this program. The results are about the electron dynamics on a three-state system with the middle state vibrating. The first figure shows the dependence of the propogation on the bridge frequency. The second figure shows the dependence of the propogation on the barrier height. The third figure shows quasienergy spectrum of the dynamical system.

Trivially, I should add that the time evolution produced by the Floquet propogator is in exact agreement with that given by the predictor-corrector time integrator. This serves as a double check for our code.

Subroutine floquet.f constructs the Floquet matrix. HP-UX version and Linux version are both available.

Some trivial programs

Statistics for T_DA at Landau-Zener crossing points

You can also find a Hamiltonian series stored on tape (for Met121-Gln107, the whole beta-sheet).

Back to the Index Page