In Frank, Computational Chemistry, Eq. 3.86 is
$$ \langle \mu_A \left| \bf{h} \right| \mu_A \rangle = \left< \mu_A \left| -\frac{1}{2}\nabla^2 - \bf{V}_a \right| \mu_A \right> - \sum_{a \neq A}^{Nuclei} \langle \mu_A \left| \bf{V}_a \right| \mu_A \rangle $$.
Is there a reason for not writing it as
$$ \langle \mu_A \left| \bf{h} \right| \mu_A \rangle = \left< \mu_A \left| -\frac{1}{2}\nabla^2 \right| \mu_A \right> - \sum_{a = A}^{Nuclei} \langle \mu_A \left| \bf{V}_a \right| \mu_A \rangle $$?
Sonntag, 10. März 2013
Dienstag, 26. Februar 2013
QM 005: Comparing semi-empirical methods
In a recent paper Liao et al. calculated the reaction profile for a newly suggested mechanism of acetylene hydratase. Conveniently, the atomic coordinates of all species along the reaction coordinate have been deposited in the supporting material.
Carrying out some SPE calculations of all species using MOPAC2012, we find the following profiles.
 |
| Fig. 1: Calculated Reaction profiles. |
In the SPE calculations, a dielectric constant of \(\epsilon=4\) was applied. In the reference ("B3LYP" in Fig. 1), the optimization was done using the LANL2TZ(f) pseudo-potential on W, 6-311+G(d) for S and 6-31G(d,p) for the other elements. The reference SPE calculations were carried out using 6-311+G(2d,2p) on all elements but W.
Apparently, the semi-empirical methods greatly overestimate the interaction energy between W and the acetylene unit (going from 0 to 1) and underestimate the activation energy for the proton transfer in step 4.
However, without reoptimizing the transition state, PM6 is does see a very clear transition state for the hydroxylation at step 2.
Orbital diagram from supplementary material.
However, without reoptimizing the transition state, PM6 is does see a very clear transition state for the hydroxylation at step 2.
Orbital diagram from supplementary material.
Mittwoch, 16. Januar 2013
PyMOL 007: New Feature
I recently wrote a mail to the PyMOL user mailing-list asking if PyMOL can write vector graphics based files of a session (similar to Molscript). Vector graphics are scalable, require less disc space (good for journal uploads) and look beautiful, allowing to focus attention on details that matter.
Here's the mail I wrote:
Here's the mail I wrote:
Hi PyMOL users Can PyMOL write a vector based picture of a session? Does not require to be very "fancy", but vector based would be cool. Something like in the old days with molscript. Best regards Martin
In the following, a number of people responded suggesting SVG or EPS/PDF as file format. Within three days, Jason Vertrees of Schrödinger started a poll to check which format was most preferred by the users. Find it here:
http://pymol.org/vector_poll
It looks as if this really is a feature that could be appreciated and it would be great if it was implemented.
I guess there will be some management involved in the decision making, but in any case its another example of the great support by the PyMOL developers at Schrödinger. And sometimes this deserves mentioning.
Samstag, 8. Dezember 2012
QM 004: Two-electron integral manipulation
In the independent-particle model, the electronic energy of an \(N\)-electron system is frequently written as
\begin{equation}E = \sum_{i}^{N} h_i + 1/2 \sum_{i, j} J_{ij} - K_{ij}.\label{eq:energy}\end{equation}
When inserting the expressions for the Coulomb and exchange operators, \(J\) and \(K\), and calculating the variation in energy, \(\delta E \), given by (in physicists notation)
\begin{equation}\delta E = \sum_i \langle \delta \phi_i | h_i | \phi_i \rangle + \langle \phi_i | h_i | \delta \phi_i \rangle + 1/2 \sum_{i, j} \left[ \langle \delta \phi_i \phi_j | \phi_i \phi_j \rangle + \langle \phi_i \delta \phi_j | \phi_i \phi_j \rangle + \phi_i \phi_j | \delta \phi_i \phi_j \rangle + \langle \phi_i \phi_j | \phi_i \delta \phi_j \rangle \\
- \langle\delta\phi_i\phi_j | \phi_j\phi_i\rangle - \langle\phi_i\delta\phi_j | \phi_j\phi_i\rangle - \langle\phi_i\phi_j | \delta\phi_j\phi_i\rangle - \langle\phi_i\phi_j | \phi_j\delta\phi_i\rangle\right], \label{eq:variation}\end{equation}
it can be found that four pairs of terms in the second sum-expression are equal and thus the factor of \(1/2\) can be cancelled.
For the Coulomb integrals this can be done by considering the definition of the bracket notation
$$ \langle \delta\phi_i(1)\phi_j(2) | \phi_i(1) \phi_j(2)\rangle = \int d{\bf x}_1 d{\bf x}_2 \delta\chi_i^*(1)\chi_j^*(2){\bf r}_{12}^{-1} \chi_i(1)\chi_j(2)$$
and noting that when the orbitals are understood to be real, i.e. \(\phi^* = \phi\), and the ordering of electron labels is considered, the orbitals of an electron \(m \in \{1, 2 \}\) can be swapped such that
$$ \int d{\bf x}_1 d{\bf x}_2 \delta\chi_i(1)\chi_j(2){\bf r}_{12}^{-1} \chi_i(1)\chi_j(2) = \int d{\bf x}_1 d{\bf x}_2 \chi_i(1)\chi_j(2){\bf r}_{12}^{-1} \delta\chi_i(1)\chi_j(2) $$
and thus, in bracket notation,
\begin{equation} \langle \delta\phi_i\phi_j | \phi_i \phi_j \rangle = \langle\phi_i\phi_j | \delta\phi_i \phi_j \rangle \label{eq:manipulation_coulomb}\end{equation}
which we recognize as the the first and third term in the second sum of Eq. \ref{eq:variation}.
For the exchange integrals, a similar expression exists (Szabo, Ostlund, Eq. 2.94)
$$\langle ij | kl \rangle = \langle ji | lk \rangle .$$
To show this, we write the integral explicitly and after first exchanging the dummy variables \(1\), \(2\) and then reordering the orbitals in order to restore the conventional \(1, 2\) order of electrons (keeping "the orbital on the electron"), we obtain
\begin{equation} \int d{\bf x}_1 d{\bf x}_2 \chi_i^*(1)\chi_j^*(2){\bf r}_{12}^{-1} \chi_k(1)\chi_l(2) = \int d{\bf x}_1 d{\bf x}_2 \chi_j^*(1)\chi_i^*(2){\bf r}_{12}^{-1} \chi_l(1)\chi_k(2). \label{eq:manipulation_exchange}\end{equation}
Using this approach
\begin{equation} \langle \delta\phi_i(1)\phi_j(2) | \phi_j(1)\phi_i(2) \rangle = \langle \phi_j(1)\delta\phi_i(2)|\phi_i(1)\phi_j(2) \rangle = \langle \phi_i(1)\phi_j(2)|\phi_j(1)\delta\phi_i(2) \rangle \label{eq:exchange}\end{equation}
where Eq. \ref{eq:manipulation_exchange} was used to obtain the first equality and Eq. \ref{eq:manipulation_coulomb} was used to swap \(\phi_j(1)\) with \(\phi_i(1)\) and \(\delta\phi_i(2)\) with \(\phi_j(2)\) to obtain the second equality.
The first and last integral are recognized as the first and fourth exchange integral of Eq. \ref{eq:variation}.
Applying these operations on the remaining integrals allows to factor out a factor of \(2\) which then cancels with the \(1/2\).
\begin{equation}E = \sum_{i}^{N} h_i + 1/2 \sum_{i, j} J_{ij} - K_{ij}.\label{eq:energy}\end{equation}
When inserting the expressions for the Coulomb and exchange operators, \(J\) and \(K\), and calculating the variation in energy, \(\delta E \), given by (in physicists notation)
\begin{equation}\delta E = \sum_i \langle \delta \phi_i | h_i | \phi_i \rangle + \langle \phi_i | h_i | \delta \phi_i \rangle + 1/2 \sum_{i, j} \left[ \langle \delta \phi_i \phi_j | \phi_i \phi_j \rangle + \langle \phi_i \delta \phi_j | \phi_i \phi_j \rangle + \phi_i \phi_j | \delta \phi_i \phi_j \rangle + \langle \phi_i \phi_j | \phi_i \delta \phi_j \rangle \\
- \langle\delta\phi_i\phi_j | \phi_j\phi_i\rangle - \langle\phi_i\delta\phi_j | \phi_j\phi_i\rangle - \langle\phi_i\phi_j | \delta\phi_j\phi_i\rangle - \langle\phi_i\phi_j | \phi_j\delta\phi_i\rangle\right], \label{eq:variation}\end{equation}
it can be found that four pairs of terms in the second sum-expression are equal and thus the factor of \(1/2\) can be cancelled.
For the Coulomb integrals this can be done by considering the definition of the bracket notation
$$ \langle \delta\phi_i(1)\phi_j(2) | \phi_i(1) \phi_j(2)\rangle = \int d{\bf x}_1 d{\bf x}_2 \delta\chi_i^*(1)\chi_j^*(2){\bf r}_{12}^{-1} \chi_i(1)\chi_j(2)$$
and noting that when the orbitals are understood to be real, i.e. \(\phi^* = \phi\), and the ordering of electron labels is considered, the orbitals of an electron \(m \in \{1, 2 \}\) can be swapped such that
$$ \int d{\bf x}_1 d{\bf x}_2 \delta\chi_i(1)\chi_j(2){\bf r}_{12}^{-1} \chi_i(1)\chi_j(2) = \int d{\bf x}_1 d{\bf x}_2 \chi_i(1)\chi_j(2){\bf r}_{12}^{-1} \delta\chi_i(1)\chi_j(2) $$
and thus, in bracket notation,
\begin{equation} \langle \delta\phi_i\phi_j | \phi_i \phi_j \rangle = \langle\phi_i\phi_j | \delta\phi_i \phi_j \rangle \label{eq:manipulation_coulomb}\end{equation}
which we recognize as the the first and third term in the second sum of Eq. \ref{eq:variation}.
For the exchange integrals, a similar expression exists (Szabo, Ostlund, Eq. 2.94)
$$\langle ij | kl \rangle = \langle ji | lk \rangle .$$
To show this, we write the integral explicitly and after first exchanging the dummy variables \(1\), \(2\) and then reordering the orbitals in order to restore the conventional \(1, 2\) order of electrons (keeping "the orbital on the electron"), we obtain
\begin{equation} \int d{\bf x}_1 d{\bf x}_2 \chi_i^*(1)\chi_j^*(2){\bf r}_{12}^{-1} \chi_k(1)\chi_l(2) = \int d{\bf x}_1 d{\bf x}_2 \chi_j^*(1)\chi_i^*(2){\bf r}_{12}^{-1} \chi_l(1)\chi_k(2). \label{eq:manipulation_exchange}\end{equation}
Using this approach
\begin{equation} \langle \delta\phi_i(1)\phi_j(2) | \phi_j(1)\phi_i(2) \rangle = \langle \phi_j(1)\delta\phi_i(2)|\phi_i(1)\phi_j(2) \rangle = \langle \phi_i(1)\phi_j(2)|\phi_j(1)\delta\phi_i(2) \rangle \label{eq:exchange}\end{equation}
where Eq. \ref{eq:manipulation_exchange} was used to obtain the first equality and Eq. \ref{eq:manipulation_coulomb} was used to swap \(\phi_j(1)\) with \(\phi_i(1)\) and \(\delta\phi_i(2)\) with \(\phi_j(2)\) to obtain the second equality.
The first and last integral are recognized as the first and fourth exchange integral of Eq. \ref{eq:variation}.
Applying these operations on the remaining integrals allows to factor out a factor of \(2\) which then cancels with the \(1/2\).
Montag, 15. Oktober 2012
QM 003: Activation free energy in PCM of SN2 reaction
Consider the reaction of chloride ion with methyl iso-cyanide to form methyl chloride and a cyanide ion. What is the activation free energy for this reaction in water?
$$Cl^{-} + H_3C-NC \rightarrow Cl-H_3C + NC^{-}$$
The activation free energy is given by the difference between the free energy of the transition state and the separated species
$$\Delta G = G^{\ddagger} - G^{Reac}.$$
The free energy in solvent is obtained by correcting the energy of the solute, \(E_{solute}\), with the solvation free energy, i.e. the energy required to place the solute in a cavity of solvent
$$G_{solution} = E_{solute} + 1/2 \int_S \sigma(\vec{r}) V(\vec{r}) d\vec{a}.$$
\(G_{solution}\) is output from a PCM calculation as Free Energy in Solvent, meaning the program computes first \(E_{solute}\) and then the solvation free energy (obtained from integrating the electrostatic potential times surface charge density over the surface of the solute).
This value is corrected for free energy contributions from translation, rotation and vibration and ZPE. The translational contribution is adjusted to report the correction for a 1mol/L solution instead of the ideal gas (which is 1mol/24.5L).
The reaction is bimolecular, meaning two particles collide to form one particle (the transition state) which then decays to products. For both reacting particles, the translational free energy correction has to be evaluated (since in solution, other than in gas phase, it is not expected that the particles form a coordinated complex [Vayner et al.]).
From HF/3-21G//HF/3-21G (the only method from which genuine transition state structures could be located), the following values are obtained.
The increase of the barrier due to the loss in translational entropy is found in the third and fourth column, where '+' indicates the summation of the values from the first and second column (i.e. S\(_{trans, 1M}\) = 69.8 cal/(mol*K) = (34.69 + 35.16) cal/(mol*K), which are obtained from separate calculations of the two species (Cl\(^-\) and H\(_3\)C-NC). This is not the same as calculating S\(_{trans, 1M}\) from a coordinated complex, which consists only from one particle.
The translational entropy (for a 1M state) of the TS however is 37.0 cal/(mol*K). At \(T=298K\), the total entropic change evaluates to (11.59 - 21.95) kcal/mol = -10.35 kcal/mol (for 1mol/24.5L) and from
$$G = H - TS$$
it is seen that this loss of entropy will lead to an increase of the barrier. The correction in translational free energy when going from the gas phase standard state to the solution standard state is moderate (around 2 kcal/mol).
It is seen that given such an approach, the barrier would be estimated to be around 54kcal/mol.
The following points are unclear to me:
$$Cl^{-} + H_3C-NC \rightarrow Cl-H_3C + NC^{-}$$
The activation free energy is given by the difference between the free energy of the transition state and the separated species
$$\Delta G = G^{\ddagger} - G^{Reac}.$$
The free energy in solvent is obtained by correcting the energy of the solute, \(E_{solute}\), with the solvation free energy, i.e. the energy required to place the solute in a cavity of solvent
$$G_{solution} = E_{solute} + 1/2 \int_S \sigma(\vec{r}) V(\vec{r}) d\vec{a}.$$
\(G_{solution}\) is output from a PCM calculation as Free Energy in Solvent, meaning the program computes first \(E_{solute}\) and then the solvation free energy (obtained from integrating the electrostatic potential times surface charge density over the surface of the solute).
This value is corrected for free energy contributions from translation, rotation and vibration and ZPE. The translational contribution is adjusted to report the correction for a 1mol/L solution instead of the ideal gas (which is 1mol/24.5L).
The reaction is bimolecular, meaning two particles collide to form one particle (the transition state) which then decays to products. For both reacting particles, the translational free energy correction has to be evaluated (since in solution, other than in gas phase, it is not expected that the particles form a coordinated complex [Vayner et al.]).
From HF/3-21G//HF/3-21G (the only method from which genuine transition state structures could be located), the following values are obtained.
| React 1 | React 2 | Sum | TS | Products | Activation Energy | |||
| PCM | Cl- | H3C-NC | Cl- + H3C-NC | TS | Cl-H3C + NC- | TS - (Cl- + H3C-NC) | ||
| HF/3-21G//HF/3-21G | E_elec | Ha | -457.5 | -131.2 | -588.6 | -588.6 | -588.6 | |
| ZPE | ZPE | Ha | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | |
| E(0K) | E_elec + ZPE | Ha | -457.5 | -131.1 | -588.6 | -588.5 | -588.5 | |
| E(0K) | kcal/mol | -287068.1 | -82277.2 | -369345.2 | -369298.3 | -369298.9 | 46.9 | |
| G_elec | kcal/mol | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | ||
| G_trans | 24.5l | kcal/mol | -9.4 | -9.6 | -19.0 | -10.1 | -10.1 | 8.9 |
| G_t(1l) = G_t(24.5l) - RTln(24.5) | RT*ln(24.5) | kcal/mol | 1.9 | 1.9 | 1.9 | 1.9 | 1.9 | |
| G_trans | 1l | kcal/mol | -11.3 | -11.5 | -22.8 | -12.0 | -12.0 | 10.8 |
| G_rot | kcal/mol | 0.0 | -5.2 | -5.2 | -6.5 | -6.5 | -1.2 | |
| G_vib | kcal/mol | 0.0 | 30.5 | 30.5 | 28.2 | 28.5 | -2.2 | |
| G_trv(298) | 1mol/24.5l | kcal/mol | -9.4 | 15.7 | 6.3 | 11.7 | 11.9 | 5.4 |
| G_trv(298) | 1mol/l | kcal/mol | -11.3 | 13.8 | 2.5 | 9.8 | 10.0 | 7.3 |
| S_trans | 24.5l | cal/(mol*K) | 36.6 | 37.1 | 73.6 | 38.9 | -34.7 | |
| S_trans | 1l | cal/(mol*K) | 34.7 | 35.2 | 69.8 | 37.0 | -32.8 | |
| T*S_trans | 298K, 24.5l | kcal/(mol) | 10.9 | 11.0 | 21.9 | 11.6 | -10.4 | |
| T*S_trans | 298K, 1l | kcal/(mol) | 10.3 | 10.5 | 20.8 | 11.0 | -9.8 | |
| E(0K) + G_trv(298) | 24.5l | kcal/mol | -287077.5 | -82261.5 | -369339.0 | -369286.6 | -369287.0 | 52.346 |
| E(0K) + G_trv(298) | 1l | kcal/mol | -287079.4 | -82263.4 | -369342.8 | -369288.5 | -369288.9 | 54.250 |
| Free energy in solvent (PCM) | kcal/mol | -287067.9 | -82307.8 | -369375.7 | -369327.1 | -369328.7 | 48.530 | |
| Free energy in solvent + ZPE + G_trv | 24.5l | kcal/mol | -287077.3 | -82261.4 | -369338.7 | -369286.4 | -369286.7 | 52.348 |
| Free energy in solvent + ZPE + G_trv | 1l | kcal/mol | -287079.2 | -82263.3 | -369342.5 | -369288.3 | -369288.6 | 54.252 |
The increase of the barrier due to the loss in translational entropy is found in the third and fourth column, where '+' indicates the summation of the values from the first and second column (i.e. S\(_{trans, 1M}\) = 69.8 cal/(mol*K) = (34.69 + 35.16) cal/(mol*K), which are obtained from separate calculations of the two species (Cl\(^-\) and H\(_3\)C-NC). This is not the same as calculating S\(_{trans, 1M}\) from a coordinated complex, which consists only from one particle.
The translational entropy (for a 1M state) of the TS however is 37.0 cal/(mol*K). At \(T=298K\), the total entropic change evaluates to (11.59 - 21.95) kcal/mol = -10.35 kcal/mol (for 1mol/24.5L) and from
$$G = H - TS$$
it is seen that this loss of entropy will lead to an increase of the barrier. The correction in translational free energy when going from the gas phase standard state to the solution standard state is moderate (around 2 kcal/mol).
It is seen that given such an approach, the barrier would be estimated to be around 54kcal/mol.
The following points are unclear to me:
- Should the free energy in solution be corrected for ZPE and free energy contributions from translation, rotation and vibration?
- Are the free energy corrections of the reactants correct when calculated from a coordinated structure (as opposed to two individual structures in separate calculations)?
- Which effects not included (apart from electronic structure properties) would be most relevant when extending the model and how could they be obtained? Conformational free energy correction from MD snap shots?
Surprisingly, the activation energy in E(0K) + G\(_{trv, 298, 1L}\) without considering solvation (54.25 kcal/mol) is the same as the activation free energy in solvent. Is this a coincidence or does this make sense or not?
Dienstag, 18. September 2012
PyMOL 006: Removing a bond between two atoms
PyMOL sometimes displays bonds which are not correct and can be confusing. To remove a bond between to atoms, first set the selection argument to "Atoms", then select the first atom involved in the bond and rename its selection identifier to, e.g., "atm1". Unselect the atom, select the second atom of the bond, rename the selection identifier to "atm2".
Then enter the following at the PyMOL prompt:
cmd.unbond("(atm1)", "(atm2)")
Be sure to include all parentheses and quotation marks.
Freitag, 29. Juni 2012
Python 005: The so called "sliding window"
Say you find yourself required to locate the line number of the last atom of residue '2' in this PDB (following the PDB internal numbering scheme, of course we assume you don't know a priori its called H01):
As you can see, it's the Atom 11 (H01) of an Alanine. How would you do that? Choose a for-loop, choose a nicely indented block of if and else statements. Choose a number of different counters, each with their very own meaningful name. Choose endless meetings with IndexErrors, comparisons, type castings and list slicings. I say no, I choose a generator. In fact, a pair wise sliding window.
So lets instantiate the object.
'res' is now the generator. Only when you call 'next(res)', the generator executes, up until the 'yield' statement, from where it returns. In addition, it can be used directly in a for-loop. The for-loop takes care of implicitly calling the 'next' method. Here's an example:
ATOM 1 N ALA 2 -0.677 -1.230 -0.491 1.00 0.00 N ATOM 2 CA ALA 2 -0.001 0.064 -0.491 1.00 0.00 C ATOM 3 C ALA 2 1.499 -0.110 -0.491 1.00 0.00 C ATOM 4 O ALA 2 2.065 -0.922 0.251 1.00 0.00 O ATOM 5 CB ALA 2 -0.509 0.856 0.727 1.00 0.00 C ATOM 6 H ALA 2 -0.131 -2.162 -0.491 1.00 0.00 H ATOM 7 HA ALA 2 -0.269 0.603 -1.418 1.00 0.00 H ATOM 8 1HB ALA 2 -1.605 1.006 0.691 1.00 0.00 H ATOM 9 2HB ALA 2 -0.285 0.342 1.681 1.00 0.00 H ATOM 10 3HB ALA 2 -0.053 1.861 0.784 1.00 0.00 H ATOM 11 H01 ALA 2 -1.261 -1.244 -1.314 1.00 0.00 H ATOM 12 N ASN 3 2.311 0.711 -1.400 1.00 0.00 N ATOM 13 CA ASN 3 3.700 0.321 -1.173 1.00 0.00 C ATOM 14 C ASN 3 4.606 1.530 -1.169 1.00 0.00 C ATOM 15 O ASN 3 4.515 2.417 -2.017 1.00 0.00 O ATOM 16 CB ASN 3 4.140 -0.693 -2.267 1.00 0.00 C ATOM 17 CG ASN 3 3.309 -1.974 -2.400 1.00 0.00 C ATOM 18 ND2 ASN 3 3.497 -2.939 -1.540 1.00 0.00 N ATOM 19 OD1 ASN 3 2.462 -2.114 -3.270 1.00 0.00 O ATOM 20 H ASN 3 1.953 1.454 -2.096 1.00 0.00 H ATOM 21 HA ASN 3 3.774 -0.146 -0.174 1.00 0.00 H ATOM 22 2HB ASN 3 5.198 -0.981 -2.124 1.00 0.00 H ATOM 23 3HB ASN 3 4.123 -0.199 -3.258 1.00 0.00 H ATOM 24 1HD2 ASN 3 2.793 -3.677 -1.623 1.00 0.00 H ATOM 25 2HD2 ASN 3 4.145 -2.756 -0.775 1.00 0.00 H ATOM 26 N PRO 4 5.629 1.649 -0.120 1.00 0.00 N ATOM 27 CA PRO 4 6.360 2.889 -0.340 1.00 0.00 C ATOM 28 C PRO 4 6.031 3.505 -1.710 1.00 0.00 C ATOM 29 O PRO 4 5.228 2.956 -2.465 1.00 0.00 O ATOM 30 CB PRO 4 7.819 2.420 -0.227 1.00 0.00 C ATOM 31 CG PRO 4 7.769 0.962 -0.705 1.00 0.00 C ATOM 32 CD PRO 4 6.440 0.451 -0.143 1.00 0.00 C ATOM 33 HA PRO 4 6.122 3.617 0.439 1.00 0.00 H ATOM 34 2HB PRO 4 8.119 2.445 0.823 1.00 0.00 H ATOM 35 3HB PRO 4 8.509 3.035 -0.810 1.00 0.00 H ATOM 36 2HG PRO 4 8.619 0.375 -0.350 1.00 0.00 H ATOM 37 3HG PRO 4 7.737 0.935 -1.796 1.00 0.00 H ATOM 38 2HD PRO 4 6.564 0.075 0.875 1.00 0.00 H ATOM 39 3HD PRO 4 5.991 -0.316 -0.778 1.00 0.00 H ATOM 40 H01 PRO 4 6.512 4.435 -2.015 1.00 0.00 H
As you can see, it's the Atom 11 (H01) of an Alanine. How would you do that? Choose a for-loop, choose a nicely indented block of if and else statements. Choose a number of different counters, each with their very own meaningful name. Choose endless meetings with IndexErrors, comparisons, type castings and list slicings. I say no, I choose a generator. In fact, a pair wise sliding window.
def make_sliding_window(val):
ln_cnt = 0
while ln_cnt < len(val):
yield ln_cnt, ln_cnt-1
ln_cnt += 1
The 'yield' statement characterizes this as a generator. What does that mean? First, a generator is *not* a list, or array or something like that. It's a function. With special properties. When called through 'next(generator_function)', it executes until it reaches a 'yield' statement, from where it returns. When called again (with a new 'next(generator_function)') it resumes from where it left of. Yes, and it returns a generator object, when you call it as such (no 'next'), but then it does not execute.So lets instantiate the object.
res=make_sliding_window(val)
'res' is now the generator. Only when you call 'next(res)', the generator executes, up until the 'yield' statement, from where it returns. In addition, it can be used directly in a for-loop. The for-loop takes care of implicitly calling the 'next' method. Here's an example:
def get_line_num_last_atom_res_2(val):
# Form a generator 'res'.
res = make_sliding_window(val)
for a, b in res:
if val[a][25] == '3':
return val[b], a
print get_line_num_last_atom_res_2(open('randompdbfile.pdb', 'r').readlines())
As you can see, we're now using it as-if it were a list. But notice, if you think it's a list, you'll not understand what is happening. It has a 'next' method, that's the important point and it's being called, even if we dont see it explicitly. So anyway, upon reaching 'yield', we get apparently two line numbers back, where it's two times the same one, but just offset by one. Now we can easily check if the upcoming atom already is part of the next residue, and if not, return the current atom. And its line number.
Abonnieren
Posts (Atom)