amber tutorials
crown tutorial

Perturbation of ion in crown ether (vacuum)

Running the Gibbs program

 
  Perturbations:

	Na+ <-> K+
	K+  <-> Cs+ 

  For example, Na+ -> K+:

	% gibbs -O \
		-i gib10.in \
		-p ../setup/crown_na_k.top \
		-c ../dynamics/md0.rst \
		-o gib_na_k.out \
		-r gib_na_k.rst

  The 'forward' and 'reverse' perturbation control files, which
  use the thermodynamic integration method, are:

		gib10.in
		gib01.in

  The perturbation output files are:

	Na+ <-> K+

		gib_na_k.out
		gib_k_na.out

	K+ <-> Cs+ 

		gib_k_cs.out 
		gib_cs_k.out

  The input file for running dynamics in the lambda=0 ('forward') state is:

	 	gib0.in

  The dynamics output files are:

		gib_k.out
		gib_cs.out

Limitations and advantages of this demo

 The main limitation of this demo is the lack of explicit waters. However, the low atom count and the classic topological simplicity of the molecule allow relatively quick runs, so that different perturbation regimes (particularly long sampling) can be tried. It is often useful to do 'prototype' runs in vacuum before investing time in slower simulations with explicit solvent.

Since the crown/ion complex is topologically simple, in vacuum it is possible that not all modes of movement are explored, leading to sampling problems. This issue is examined further in the next section. Note: In molecules with freely-rotating groups such as methyls, the danger of getting trapped in a limited sampling regime in vacuum is more pronounced, i.e. methyls can spin very rapidly, while the molecule overall moves in a much lower-temperature regime.

Solvating this system would introduce an ion parameterization problem. This is because, with current methods, one must choose ionic parameters consistent with either water or the atom types in the crown molecule: it is impossible to allow the ion to act according to different nonbonded constants with different atoms. In the case of the crown, the issue could be crucial if the waters provide much of the free energy difference between different ion complexes; consider the structural finding that waters can fit in the ring with smaller ions. The reason for concern about this particular blending of water and solute force fields is that the TIP/SPC waters are oxygen-centered spheres somewhat larger than oxygen.

The optimal solution would be to add the 'separate parameters' feature to the software; the lesser is to try both ion parameter sets in the vacuum perturbation and in a pure water perturbation, and see whether the differences warrant further worry.

Sampling/hysteresis in Na+ <-> K+

 Sampling is always an issue with free energy studies: in theory the longer the dynamics run and the more gradual the perturbation, the better the chance that conformational space will be sampled adequately. One clear sign that sampling is inadequate is a difference in the free energies of the 'forward' and 'reverse' perturbations (this difference is often referred to - rather loosely - as 'hysteresis'). However, it is not unusual to get an apparently nicely-converged perturbation that diverges with longer sampling time. Therefore it is common to perform a perturbation with three or more sampling times in order to explore this effect.

A large amount of cpu time was expended at this point to show how various parameters, notably sampling, affect hysteresis. Eventually it was noticed that temperatures fell from about 250K to 50K over two nanoseconds whenever SHAKE was used, and the work was re-done. The original work. | An investigation of the SHAKE energy drain using dynamics.

Here is an exploration of how conditions affect hysteresis and free energies in the Na+ <-> K+ perturbation. Note that, since this is in vacuum, any implications that can be discerned may not be applicable for solvated systems.

In the long constant temperature run using SHAKE, the complex drifted so far from the origin that the coordinates exceeded the format of the restrt file, so center of mass motion was removed periodically (nscm=1000) for all constant temperature runs.

Free energy / hysteresis for Na+ <-> K+ in crown ether
for various regimes and sampling times.
Equilibration steps for each window: 100. 
E == constant energy, no shake
T == constant temp, shake on H, tol=0.000005
R1 == same as T, with velocity rescaling (ntt=-100, dtemp=100)
R2 == no shake, 1000 steps equil, center & temp rescaling before equil

        ------------------ data steps per window -----------------------------
        200         1000       10000       20000       200000
        ----------------------------------------------------------------------

E/.5fs  -11.0      -10.8       -10.1      -10.0        -10.2
          0.9        0.1         0.2        0.7          0.1

E/1fs   -11.6      -12.2       -11.2       -9.1         -8.3
          2.8        0.7         2.3        1.9          0.4

T/1fs   -10.9       -9.9        -9.6      -10.2         -8.9
          0.7        2.0         2.1        0.6          0.5

T/2fs    -9.2      -10.4        -8.7      -10.1         -9.1
          0.9        0.5         0.7        0.1          0.7

R1/2fs  -10.4       -9.6        -9.8       -9.2         -9.1
          0.6        2.1         0.4        1.7          0.2


R2/.5fs            -10.8        -9.3      -10.1         -9.1
                     0.1         1.8       -0.3          1.2
Forward/reverse energies | Ensemble properties of the 200000-data-step runs

The variation in these results is perplexing: there is no telling if any conditions reliably yield the "right" result. At this point, a study of equipartition was performed using molecular dynamics, to see if periodic velocity assignment might yield better coverage of the conformational space. The upshot was the discovery that there appear to be at least two relatively stable crown conformations around the Na+ ion: a relatively stiff, flat circle; and a folded, more flexible conformation. Accordingly, more perturbations were run, using the folded conformation as a starting point, to see if this would reduce the variability of the results.

It turns out that the expanding Na+ ion forces the initial folded conformation to the flat circle by the time the K+ state is reached. Thus the initial crown conformation would be expected to mostly affect the Na+ -> K+ perturbation, and the K+ -> Na+ result would be expected to vary depending on whether the crown happened to fold as the ion shrank.

In the graph below, each row of numbers that is labeled on the left margin refers to perturbations from the planar crown configuration, while the next, unlabeled row is from the corresponding perturbations from the folded configuration.

Free energy / hysteresis for Na+ <-> K+ in crown ether
Comparing 'flat' and 'folded' starting configurations
All runs constant temperature = 300K 
        ------------------ data steps per window -----------------------------
        200         1000       10000       20000       200000
        ----------------------------------------------------------------------

---- SHAKE on H, 100 steps equil

-- step=1fs

1->0    -11.2      -10.9       -10.6       -9.9         -8.6
        -11.9      -10.9        -9.7      -10.0         -9.2

0->1    -10.5       -8.9        -8.5      -10.5         -9.1
        -10.3      -10.3        -8.4       -8.5         -8.7

mean    -10.9       -9.9        -9.6      -10.2         -8.9
        -11.1      -10.6        -9.1       -9.3         -9.0

hyst      0.7        2.0         2.1        0.6          0.5
          1.6        0.6         1.3        1.5          0.5

-- step=2fs

1->0     -9.7      -10.6        -9.0      -10.0         -9.4
        -11.8      -10.9       -10.2      -10.0         -8.4

0->1     -8.6      -10.1        -8.3      -10.1         -8.7
        -10.3      -10.3        -8.4       -8.4         -8.4

mean     -9.2      -10.4        -8.7      -10.1         -9.1
        -11.1      -10.6        -9.3       -9.2         -8.4

hyst      0.9        0.5         0.7        0.1          0.7
          1.5        0.3         1.8        1.6          0.0

---- 1000 steps equil, reassign velocities at each window

-- step=0.5fs

1->0               -10.7       -10.2       -9.9         -9.7
                   -11.6       -10.3       -9.6         -9.7

0->1               -10.8        -8.4      -10.2         -8.5
                   -10.2       -10.1       -8.3         -8.6

mean               -10.8        -9.3      -10.1         -9.1
                   -10.9       -10.2       -9.0         -9.2

hyst                -0.1         1.8       -0.3         -1.2
                     1.4         0.2        1.3          1.1
Ensemble properties of the 200000-data-step runs
From these numbers, we see that the starting conformation does not significantly affect the results for this system. However, now that several series of runs have been accumulated, we can look for trends, e.g. by comparing step sizes. In the first graph below, the red and blue plots are the results for the two 0.5fs, constant temperature runs above, while the black plot is for the contant energy runs in the first table.

The plots are colored {black, red, blue} according to their order in the raw results file.

It appears that the 0.5fs time step leads to more uniform results after a given period, although less simulated time for sampling could account for this.

Temperature in K+ <-> Cs+

For K+ <-> Cs+, the problem is that even the intermediate-sized Rb+ cannot fit in the ring, and in vacuum, when it pops out of the ring, it drifts away. The solution was to run the dynamics at 10K. It turns out that, for this system, the Na+ <-> K+ perturbation gives results consistent with those above for the shortest protocol used at 300K, 100K, and 10K:

Na+ <-> K+ free energy / hysteresis
with constant temperature
step = 2 fs 
                                 Temperature
                    -------------------------------------
                    300K             100K             10K

     mean           -9.2            -10.0            -9.8
     hysteresis      1.1              0.1             0.0
Forward/reverse energies

Therefore all the perturbations are run at 10K for this vacuum "toy problem." Note that the reduction of hysteresis with lower temperatures is due to less sampling, because at the lower temperature the system is much less apt to reach higher-energy states; possibly this also prevents adequate sampling of low energy states as well, if they would only be reached by crossing energy barriers higher than the 10K ensemble would be expected to reach in the chosen sampling time.

Demo by Bill Ross.