PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Phys Chem A. Author manuscript; available in PMC Dec 3, 2008.
Published in final edited form as:
PMCID: PMC2593116
NIHMSID: NIHMS63539
First-Principles Simulation of Amide and Aromatic Side Chain Ultraviolet Spectroscopy of a Cyclic Dipeptide
Zhenyu Li and Shaul Mukamel
Department of Chemistry, University of California, Irvine, CA 92697 USA
By combining time-dependent density functional theory (TDDFT) and molecular dynamics (MD) simulations, we calculate the ultraviolet absorption and circular dichroism (CD) of a cyclic dipeptide, cyclo(L-Pro-D-Tyr), in the 185−300 nm region. The absorption is dominated by the phenol chromophore of tyrosine. The CD spectrum shows both phenol and amide units transitions. A crude coherent two dimensional ultraviolet spectrum (2DUV) calculated by neglecting the two-excitation states shows a cross peak between two transitions of the phenol in the tyrosine side chain. Additional cross peaks between the side chain and the backbone are observed when using a chirality-induced pulse polarization configuration.
Optical techniques are widely used for the structural characterization and probing ultrafast dynamics of proteins.1 Infrared spectra monitor the vibrations, whereas visible spectra probe specific chromophores such as porphyrins or retinal. There are two primary types of UV transitions in proteins: the first is the peptide backbone nπ* and ππ* transitions in the far-UV range, and the second one is the aromatic side chains. There are only three amino acids with aromatic side chains, phenylalanine (Phe), tyrosine (Tyr), and tryptophan (Trp). According to Platt's notation,2 there are four valence electronic excitations for each aromatic side chain, namely La, Lb, Ba, and Bb. Aromatic side chains are relatively rare, making them easier to characterize, and better local probes for structure. UV resonant Raman excitation profiles are very sensitive to protein structure.3-7
Theoretical simulations are necessary for correlating spectra with molecular structure. The exciton model has been widely used for the interpretation of UV spectra of protein.8,9 The necessary parameters are typically obtained from detailed studies on small model molecules.10 Thanks to their simplicity and structural rigidity, cyclic dipeptides, also known as diketopiperazines (DKP) or 2,5-piperazinediones, have been an attractive model system for characterizing the interactions between two amide groups.11-24 The simplest cyclic dipeptide is the unsubstituted DKP, i.e. cyclo(Gly-Gly). Its vibrational and electronic spectra have been measured25-28 and calculated using correlated ab initio techniques.29 By examining the splitting between the two ππ* transitions, Hirst et al.29 concluded that the amide-amide interactions are primary electrostatics.
CD is widely used for probing the secondary structure and folding dynamics of proteins.30,31 It will be interesting to study chiral substituted DKP, which are optically active (unsubstituted DKP is non-chiral). In this article, we consider the cyclic dipeptide PTDKP composed by proline and tyrosine, i.e. cyclo(L-Pro-D-Tyr) [see Fig. 1]. Its circular dichroism (CD) spectrum has been measured.19
FIG. 1
FIG. 1
Optimized geometry of PTDKP. Nitrogen (blue), oxygen (red), carbon (green), and hydrogen (white).
Proline is a unique amino acid, since its side chain is covalently bonded to the nitrogen atom of the peptide backbone. As suggested by CD, the cyclic five-membered ring of proline imposes rigid constraints on the structure of the DKP ring.19 For unsubstituted DKP, two boat enantiomers lie in a shallow well, separated by a very small barrier corresponding to the planar structure.29 However, cyclic dipeptides containing proline have a stable boat configuration. On the other hand, the poly-peptide of proline can form the PII conformation,32 which shows a very different CD pattern than α-helix and β-sheet. It is hard to simulate the CD spectra of the PII structure using an exciton model, even with several polarizability parameters.33 Cyclic dipeptides which contain a proline show similar CD pattern to the PII structure,19 making them good candidates for theoretical modeling.
The other amino acid in PTDKP, tyrosine, has an aromatic side chain with a phenol chromophore. The phenol La excitation is close in energy to the backbone nπ* transition, and its Ba and Bb transitions are close to the amide ππ* transition. The aromatic side chain spectra are very sensitive to their immediate environment, making them excellent probes of protein structure. PTDKP is therefore also an ideal model system for studying backbone-side chain interactions.
Except for the simplest cyclo(Gly-Gly), previous theoretical modelings of cyclic dipeptide, were carried out using model Hamiltonians. First principles simulations of spectra of molecule of this size are now feasible.34 In this article, we simulate the spectra of PTDKP by combining time-dependent density functional theory (TDDFT) with first-principles molecular dynamics (MD). Absorption and CD are studied first. Beyond these linear spectroscopic techniques, two-dimensional (2D) coherent third-order optical techniques may reveal additional details.35-38 Such techniques are becoming available due to the recent development of ultrafast UV laser sources. The cross peaks in 2DUV spectra may provide a good probe for the backbone-side chain interaction, which may not be probed by CD directly. Previously, we have used a model hamiltonian to simulate 2DUV of protein. We demonstrated that the technique can distinguish between two different parameterized Hamiltonians both providing a good fit for CD.39 Here, simulated 2DUV of PTDKP show clear cross peaks between side chain and backbone. We describe the computational details in section II. The electronic structure, the MD fluctuations, and the simulated UV spectra are presented in section III. Finally, we conclude in section IV.
The simulations were performed in three steps. First, a MD simulation is used to generate a trajectory of PTDKP on the ground state potential energy surface. Repeated TDDFT calculations then provide the excited states information at different geometries. Finally, the spectra are simulated using the sum over states expressions for the response functions.35
Molecular dynamics simulations
We solve the classical equations of motion of the nuclei in 1 fs time steps by means of the Verlet leapfrog algorithm,40 as available in the FROG module41 of TURBOMOLE program suite.42 SV(P) basis set43 were used. Energies and forces were evaluated with the BP86 functional44,45 in combination with the RI-J approximation46,47 and a small quadrature grid (size 2).48 Implicit solvent effects were taken into account by the COSMO model49 using the dielectric constant of water. The Nose-Hoover thermostat50 is included to run constant temperature MD. First, a 2 ps trajectory is generated at 800 K, and 20 snapshots are taken from this trajectory. Starting from these initial configurations, 20 trajectories were launched at 300 K. Each trajectory comprised 2000 1 fs time steps, corresponding to a total simulation time of 40 ps. After an initial 200 step equilibration, we sampled structural parameters and computed the UV spectra every twenty steps, for a total of N=1800 snapshots.
Time-dependent density functional calculations
Using TDDFT calculations for the 1800 snapshots, we extracted 30 excitations to cover the 185 − 300 nm range. The hybrid density functional of Perdew, Burke, and Ernzerhof (PBE0)51 was used. It contains 25% of Hartree-Fock exchange, is more robust for higher excitation energies52 and is less susceptible to charge-transfer error than other functionals. SV(P) basis set were used. Our test on one snapshot indicated that SV(P) gave similar result with larger TVZP basis set in the frequency range considered. The TDDFT response calculations were based on ground state DFT wavefunctions with COSMO continuum solvent model. All calculations were carried out with TURBOMOLE.
Spectra simulation
Spectra were calculated from the response functions,35 by using the SPECTRON package.53 Inhomogeneous broadening due to slow fluctuations is assumed. Fluctuations are taken from fully quantum-mechanical Born-Oppenheimer MD simulations, where ensemble averages are replaced by time averages.54,55 For each snapshot, the response function is calculated by sum over the lowest 30 excited states. Transition energies, transition dipoles, and magnetic dipoles of these excited states to the ground state are obtained from the TDDFT. There is a minor difference between the velocity and length representations of a transition dipole, the mean value was used in our calculations. Two-exciton states are not available from linear response TDDFT within the adiabatic approximation,56,57 and were neglected. This does not affect absorption and CD spectra, but should affect the 2DUV signals.
A. Geometry and electronic structure
The optimized PTDKP geometry with solvent effects included by the COSMO continuum model is shown in Fig. 1. As suggested by previous force field calculations,19 the DKP ring in the optimized structure has a boat configuration. Its mean deviation from a plane is 0.096 Å. The phenol ring is almost planar, with mean deviation 0.002 Å. For such systems, it has been suggested that, due to interactions between aromatic and peptide π systems, the folded conformer predominates for several solvents and temperatures.20 Instead of using the two torsion angles γ1 and γ258 to describe the folding configuration, we used a single parameter, the folding angle θ between the two least-squares planes of the DKP ring and the phenol ring. The optimized geometry gives θ = 59.0 .
Based on the optimized geometry, the electronic structure of PTDKP is calculated with PBE0 hybrid density functional. A 6.07 eV gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is obtained. As shown in Fig. 2, the HOMO (orbital 69) is a π orbital within the phenol ring of tyrosine. LUMO (orbital 70) and LUMO+1 are both antibonding π orbitals of the phenol ring. Lb and La represent transitions from HOMO to these two orbitals. We shall label the Pro-Tyr peptide unit and the Tyr-Pro peptide unit as 1 and 2 respectively. Orbital 68 is then a combination π2 and n1. Orbital 67 is a combination of n2 and π1. Orbital 65 is mainly the phenol π orbital, with small contributions from n1 and n2 electrons. Orbital 64 is more complicated, it is a combination of n2, n1, π1, and πphenol. For unoccupied orbitals, the mixing between the two peptide units is smaller than in the lower occupied states. Orbital 72 is equation M1, and orbital 73 is equation M2.
FIG. 2
FIG. 2
Molecular orbitals and their energy levels. The highest occupied molecular orbital (HOMO, orbital 69) is at −6.48 eV, and the lowest unoccupied molecular orbital (LUMO, orbital 70) is at −0.41 eV.
B. Structural fluctuations from ab initio molecular dynamics
Geometry optimization may depend on the starting configuration. To generate a good ensemble of structures, we carried out an MD simulation. CD measurements suggest that the PTDKP conformation does not vary strongly for different solvents.19 The Continuum solvent model should thus be adequate.
Since the DKP and phenol rings are the chromophores, we shall mainly consider fluctuations of the rigidity of these two rings and their folding configuration. The ring rigidity is measured by the mean deviation of the six atoms on the ring from their least-squares plane. These will be denoted Δ1 (DKP) and Δ2 (phenol). The distributions of the mean deviations of these two rings are displayed in Fig. 3. As expected, the phenol ring is almost planar, with mean Δ2 about 0.03 Å. It is also very rigid, and the variance of Δ2 is only 0.014 Å. The dipeptide ring is more flexible, with the mean value and variance of Δ1 0.12 and 0.032 Å, respectively. The folding angle θs centered at 50.6 degree (variance 7.48 degree) indicates a folded configuration. No jump from folded to unfolded configurations is observed in our MD simulations.
FIG. 3
FIG. 3
(a) Distributions of the distortions of the dipeptide ring Δ1 (red) and the phenol ring Δ2 (blue). (b) Distribution of the folding angle θ.
C. Absorption and CD spectra
The simulated absorption spectrum is depicted in Fig. 4A. A 150 cm1 Lorentzian broadening (FWHM) is added to each snapshot. There are mainly three peaks in the spectrum. We can assign these peaks from the stick spectrum by examining the excitation configurations in the optimized geometry. The absorption is dominated by the phenol side chain of tyrosine. The strongest ~184 nm peak at the blue side is the Ba transition. The second ~221 nm peak is the La state, and the weakest ~252 nm peak corresponds to the Lb transition. The peptide nπ* transitions are forbidden. One of the two peptide ππ* transitions is covered by the phenol Ba peak, the other has higher energy, outside the frequency range of our plot.
FIG. 4
FIG. 4
Simulated (a) absorption and (b) circular dichroism of PTDKP. Black line is the experimental CD.19 The intensity of the experimental CD in the 250−300 nm range is multiplied by 40. Stick spectra are calculated from optimized PTDKP geometry.
The CD spectrum simulated with the same protocol is displayed in Fig. 4B. Its main features agree with experiment.19 At the blue side, ~192 nm, there is a strong negative peak, which mainly comes from equation M3 transition within the Tyr-Pro peptide unit. This is the main feature of all proline containing cyclic dipeptides. It is interesting to note that the oscillator strength of Ba is much larger than equation M4, but the rotatory strength is weaker. Therefore, by combining the absorption and CD spectra, we can determine the transition energies of both states. The red side is dominated by the negative Lb peak, which is stronger and blue shifted compared to experiment. The measured intensity of Lb strongly depends on the type of solvent,19 and the difference in intensity between our calculation and experiment may be caused by the continuum model used for the solvent. Another possible reason is the under-estimated frequency of the TDDFT charge transfer state, as will be discussed later.
Between equation M5 and Lb, there are three additional states, namely the two peptide nπ* states and the phenol La state. The two nπ* states are almost degenerate. In experiment,19 only one positive peak is observed between the two negative peaks, and there is a small frequency range around 250 nm where no experimental data are available. Madison et al.19 had assigned this positive peak to both nπ* and La. Our simulated CD show two positive and one negative peak between equation M6 and Lb. The stick spectrum analysis indicates that the positive 234 nm CD peak is mainly nπ*, the other positive peak is mainly a charge transfer state between the dipeptide and the phenol rings, and the negative peak is La. Our TDDFT calculations give an incorrect order for La and nπ*, and shift an high energy charge transfer state to the experimental frequency range. TDDFT is well known to underestimate charge transfer state energies, compared with higher level of theory.59,60 In our TDDFT calculations, the charge transfer states mix with the nπ* and La states. This may be the reason why their order is reversed.
The negative La CD peak is slightly shifted with respect to the corresponding absorption peak. This is because the negative La CD peak is close to the positive nπ* peak, and they interfere. Therefore, the nπ* CD peak should be slightly blue shifted compared to the transition energy (a reasonable estimate is 231 nm). By combining the absorption and the CD, we obtained ensemble averaged TDDFT transition energies for all transitions in the experimental frequency range, and it is shown in Table I. These will be used in the next section to assign cross peaks in 2D signals.
Table 1
Table 1
TDDFT ensemble averaged transition frequencies and experimental CD peak positions. Both are in nm.
D. Two dimensional electronic correlation spectroscopy
A 2D photon echo (2DPE) experiment uses three pulses with wavevectors k1, k2, and k3, which interact with the molecule to generate a coherent signal in the direction ks = −k1 + k2+k3.35-38 There are three time delays (t1, t2, and t3) between the three incoming pulses and the signal pulse. 2D correlation plots of the signals are plotted as Fourier transforms with respect to t1 and t3, and we set t2 = 0. Different polarization and wavevector configurations are denoted by μ4μ3μ2μ1(α4α3α2α1), where μi is the polarization direction of field i and αi is its wavevector direction. The αi indices will be omitted when the signal does not depend on wavevectors of the field.
We have simulated the two dimensional photon echo spectrum of PTDKP by calculating third order response function (Eq. 3.6 of Ref.61). The imaginary part of the signal for the all parallel xxxx pulse polarization configuration is shown in the top panel in Fig. 5. Along the diagonal line, as expected, the 2D signal is very similar to linear absorption (see Fig. 6a), and is dominated by two peaks from phenol La and Ba transitions. Most interesting is the cross peak observed between these two transitions. The diagonal 2DPE peaks are broad along the diagonal direction, but narrower in the antidiagonal direction. This characteristic selective elimination of inhomogeneous broadening along the antidiagonal direction in 2DPE allows us to extract information about spectral fluctuations and dephasing.
FIG. 5
FIG. 5
Top: Simulated xxxx 2DUV signal. Diagonal line is marked by dotted line. The three solid grid lines indicate energies of phenol side chain transitions Ba (184 nm), La (221 nm), and Lb (252 nm), respectively. The two dashed grid lines correspond to energies (more ...)
FIG. 6
FIG. 6
(a) Absorption (red) and the minus diagonal section (blue) of the xxxx 2DUV signal. (b) CD spectrum (red) and the diagonal section (blue) of the chiral xxxy(zzzz) signal.
Chirality induced 2DUV signals are extension of CD, which have been predicted to be a powerful technique.62,63 These signals also depend on the wavevectors of the fields. We adopted the configurations where all fields propagate along z, and simulated the real part of the chiral 2DPE xxxy(zzzz) signal. The chiral signal is weaker than the nonchiral one, but it contains more information. As shown in Fig. 6b, the diagonal section of the chiral 2DPE is similar to CD spectrum. It is interesting to note that the charge transfer peak in the simulated CD at 205 nm does not appear in the diagonal spectrum of the xxxy(zzzz) signal. In the nondiagonal region, we see many cross peaks. The upper triangle region shows cross peaks between nπ* and La, nπ* and ππ*, and nπ* and Ba. In the lower triangle region, we have cross peaks between nπ* and ππ*, La and ππ*, and ππ* and Ba.
An extremely interesting feature of the chirality induced 2DUV signals is that they show cross peaks between the aromatic side chain and the backbone, which are absent in the nonchiral xxxx signal. Cross peaks reflect the interactions between transitions. Therefore, this UV technique can directly probe the interaction between side chains and backbone. Amino acids with aromatic chromophores are relatively rare in proteins, and can thus provide a site specific information, as is done by isotopic labeling in the infrared spectroscopy.
In our simulations of 2DUV, the two-excitation states are not included, since they are not available from TDDFT. This may be justified if two-exciton couplings are very large, shifting these states outside the spectral range of interest. At this point, we can not tell how good this approximation is for PTDKP. It is a formidable challenge, however, to calculate the two-excitation states accurately with an affordable computational cost. Higher level electronic structure techniques which can describe the two-excitation states, such as full configuration interaction (CI), multi-reference CI (MRCI), second order multiconfigurational perturbation theory (CASPT2), and couple-cluster singles and doubles (CCSD), are only feasible for very small systems. Some attempts to go beyond the adiabatic approximation in TDDFT to include double excitations have been made recently.56,64 2DUV could provide a powerful probe for two excitation states.65
The electronic structure and spectra of PTDKP are studied at the ab initio TDDFT level. This cyclic dipeptide mainly assumes the folded configuration in water. The absorption is dominated by the phenol side chain excitations. The simulated CD reproduces the two negative experimental peaks, while an extra charge transfer state appears between them. TDDFT gives the wrong order for nπ* and La transitions. The simulated photon echo spectrum shows peaks from phenol La and Ba states, and their cross peak. Chiral 2DPE further enhances the resolution and reveals numerous additional features. Several cross peaks are observed. We expect the possibility to observe 2DUV cross peaks between backbone and side chains to be very useful for probing protein structure and folding dynamics.
acknowledgements
We wish to thank Dr. Darius Abramavicius and Dr. Filipp Furche for many helpful discussions. This work was supported by the National Institutes of Health Grant GM59230 and the National Science Foundation Grant CHE-0446555.
1. Hammes GG. Spectroscopy for the Biological Sciences. Wiley; New York: 2005.
2. Platt JR. J. Chem. Phys. 1949;17:484.
3. Lednev IK, Karnoup AS, Sparrow MC, Asher SA. J. Am. Chem. Soc. 1999;121:8074.
4. Lednev IK, Karnoup AS, Sparrow MC, Asher SA. J. Am. Chem. Soc. 1001;123:2388. [PubMed]
5. Asher SA, Mikhonin A, Bykov S. J. Am. Chem. Soc. 2004;126:8433. [PubMed]
6. Ahmed Z, Beta IA, Mikhonin AV, Asher SA. J. Am. Chem. Soc. 2004;127:10943. [PubMed]
7. Mikhonin AV, Asher SA. J. Am. Chem. Soc. 2006;128:13879. [PubMed]
8. Bayley PM, Nielsen EB, Schellman JA. J. Phys. Chem. 1969;73:228. [PubMed]
9. Woody RW. J. Chem. Phys. 1968;49:4797. [PubMed]
10. Besley NA, Hirst JD. J. Am. Chem. Soc. 1999;121:9636.
11. Balasubramanian D, Wetlaufer WB. J. Am. Chem. Soc. 1966;88:3449.
12. Urry DW. Annu. Rev. Phys. Chem. 1968;19:477.
13. Kopple KD, Ohnishi M. J. Am. Chem. Soc. 1969;91:962. [PubMed]
14. Sletten E. J. Am. Chem. Soc. 1970;92:172.
15. Strickland EH, Wilchek M, Horwitz J, Billups C. J. Bio. Chem. 1970;245:4168. [PubMed]
16. Karle IL. J. Am. Chem. Soc. 1972;94:81. [PubMed]
17. Karle IL, Ottenheym HC. J. Am. Chem. Soc. 1974;96:539.
18. Snow JW, Hooker TM., Jr. J. Am. Chem. Soc. 1975;97:3506. [PubMed]
19. Madison V, Young PE, Blout ER. J. Am. Chem. Soc. 1976;98:5358. [PubMed]
20. Young PE, Madison V, Blout ER. J. Am. Chem. Soc. 1976;98:5365. [PubMed]
21. Sammes PG, Weedon AC. J. Chem. Soc., Perkin Trans. 1. 1979:3048.
22. Bowman RL, Kellerman M, Hohnson WC., Jr. Biopolymers. 1983;22:1045.
23. Fleischhauer J, Grotzinger J, Kramer B, Kruger P, Wollmer A, Woody RW, Zobel E. Biophys. Chem. 1994;49:141. [PubMed]
24. Besley NA, Brienne M-J, Hirst JD. J. Phys. Chem. B. 2000;104:12371.
25. Kaya K, Nagakura S. Theor. Chim. Acta. 1967;7:124.
26. Nielsen EB, Schellman JA. J. Phys. Chem. 1967;71:2297. [PubMed]
27. Kaya K, Nagakuru S. J. Mol. Spectrosc. 1972;44:279.
28. Song S, Asher S, Krimm S, Shaw KD. J. Am. Chem. Soc. 1991;113:1155.
29. Hirst JD, Persson BJ. J. Phys. Chem. A. 1998;102:7519.
30. Nakanishi K, Berova N, Woody RW. Circular Dichroism Principles and Applications. VCH; New York: 1994.
31. Fasman GD. Circular Dichroism and the Conformational Analysis of Biomolecules. Plenum; New York: 1996.
32. Creighton TE. Proteins, structures and Molecular Properties. 2rd Ed. W. H. Freeman; New York: 1993.
33. Woody RW. Pravite Communication.
34. Frelek J, Kowalska P, Masnyk M, Kazimierski A, Korda A, Woznica M, Chmielewski M, Furche F. Chem. Eur. J. 2007;13:6732. [PubMed]
35. Mukamel S. Principles of Nonlinear Optical Spectroscopy. Oxford University Press; New York: 1995.
36. Mukamel S. Ann. Rev. Phys. Chem. 2000;51:691. [PubMed]
37. Asplund MC, Zanni MT, Hochstrasser RM. Proc. Natl. Acad. Sci. USA. 2000;97:8219. [PubMed]
38. Brixner T, Stenger J, Vaswani HM, Cho M, Blankenship RE, Fleming GR. Nature. 2005;434:625. [PubMed]
39. Li Z, Abramavicius D, Zhuang W, Mukamel S. Chem. Phys. 2007 in press. DOI: 10.1016/j.chemphys.2007.03.029.
40. Verlet L. Phys. Rev. 1967;159:98.
41. Elliott SD, Ahlrichs R, Hampe O, Kappes MM. Phys. Chem. Chem. Phys. 2000;2:3415.
42. Ahlrichs R, Bar M, Haser M, Horn H, Kolmel C. Chem. Phys. Lett. 1989;162:165.
43. Schafer A, Horn H, Ahlrichs R. J. Chem. Phys. 1992;97:2571.
44. Becke AD. Phys. Rev. A. 1988;38:3098. [PubMed]
45. Perdew JP. Phys. Rev. B. 1986;33:8822. [PubMed]
46. Eichkorn K, Treutler O, Ohm H, Haser M, Ahlrichs R. Chem. Phys. Lett. 1995;242:652.
47. Eichkorn K, Weigend F, Treutler O, Ahlrichs R. Theor. Chem. Acc. 1997;97:119.
48. Treutler O, Ahlrichs R. J. Chem. Phys. 1995;102:346.
49. Klamt A, Schurmann G. J. Chem. Soc. Perkin Trans. 1993;5:799.
50. Hoover WG. Phys. Rev. A. 1985;31:1695. [PubMed]
51. Perdew JP, Ernzerhof M, Burke K. J. Chem. Phys. 1996;105:9982.
52. Adamo C, Scuseria GE, Barone V. J. Chem. Phys. 1999;111:2889.
53. Zhuang W, Abramavicius D, Hayashi T, Mukamel S. J. Phys. Chem. B. 2006;110:3362. [PMC free article] [PubMed]
54. Allen MP, Tildesley DJ. Computer Simulation of Liquids. Oxford University Press; Oxford: 1987.
55. Frenkel D, Smit B. Understanding Molecular Simulation: From Algorithms to applications. 2nd Ed. Academic Press; San Diego: 2001.
56. Maitra NT, Zhang F, Cave RJ, Burke K. J. Chem. Phys. 2004;120:5932. [PubMed]
57. Casida M. J. Chem. Phys. 2005;122:54111. [PubMed]
58. IUPAC Biochemistry. 1970;9:3471. [PubMed]
59. Dreuw A, Head-Gordon M. Chem. Rev. 2005;105:4009. [PubMed]
60. Oakley MT, Hirst JD. J. Am. Chem. Soc. 2006;128:12414. [PubMed]
61. Abramavicius D, Mukamel S. Chem. Rev. 2004;104:2073. [PubMed]
62. Abramavicius D, Mukamel S. J. Chem. Phys. 2006;124:34113. [PubMed]
63. Abramavicius D, Zhuang W, Mukamel S. J. Phys. B: At. Mol. Opt. Phys. 2006;39:5051.
64. Cave RJ, Zhang R, Maitra NT, Burke K. Chem. Phys. Lett. 2004;389:39.
65. Yang L, Mukamel S. Phys. Rev. Lett. 2007 submitted.