text:
# Cooperativity in Two-State Protein Folding Kinetics

## I Introduction

Protein folding kinetics is usually modeled in either of three ways. First, there are mass-action models that capture the amplitudes and decay rates of the exponentials in the folding or unfolding relaxation process Ikai and Tanford (1971); Tsong et al. (1971); Dill and Chan (1997); Englander (2000). Mass-action models are useful for cataloging the different types of kinetic behavior, but give no insight into molecular structures or mechanisms. Such models do not predict other experimental properties, such as $`\mathrm{\Phi }`$-values. Second, there are all-atom or lattice model simulations that can explore sequence-structure relationships (see, e.g., Duan and Kollman (1998); Shea and Brooks (2001); Daggett (2002)). They are usually limited by computational power to short time scales and to studying restricted conformational ensembles. Third, between these macroscopic and microscopic extremes, another type of model has recently emerged. This class of models uses knowledge of the native structure to infer the sequences of folding events Munoz and Eaton (1999); Alm and Baker (1999); Galzitskaya and Finkelstein (1999); Debe and Goddard (1999); Shoemaker et al. (1999); Clementi et al. (2000); Hoang and Cieplak (2000); Li and Shakhnovich (2001); Portman et al. (2001); Ivankov and Finkelstein (2001); Alm et al. (2002); Klimov and Thirumalai (2002); Flammini et al. (2002); Bruscolini and Pelizzola (2002); Bruscolini and Cecconi (2002); Micheelsen et al. (2003). Some of these models define partially folded states with one or two contiguous sequences of native-like ordered residues Munoz and Eaton (1999); Alm and Baker (1999); Alm et al. (2002); Galzitskaya and Finkelstein (1999). Others are based on a Go-model energy function that enforces the global stability of the native state Clementi et al. (2000); Hoang and Cieplak (2000); Li and Shakhnovich (2001).

We describe here a folding model of the third type. Our model uses knowledge of the native structure to predict the kinetics. However, it differs from previous models in several respects. First, our model focuses on chain entropies and estimates loop lengths from the graph-theoretical concept of effective contact order ECO (see below). We follow time sequences of loop-closure events because we expect that these events reveal how the kinetics is encoded in the native structure. We assume that folding proceeds mostly through closures of small loops, and that large-loop closures are much slower and less important processes. Second, our model focuses on contacts within the chain, not on whether residues are native-like or not Munoz and Eaton (1999); Alm and Baker (1999); Galzitskaya and Finkelstein (1999), because we think the formation of contacts is a more physical description of the folding process. Therefore, in our model partially folded states are characterized by formed contacts, not by contiguous stretches of native-like ordered residues as in other simple models. Third, the folding kinetics is described by a master equation that can be solved directly for the macrostates considered here, without stochastic simulations such as molecular dynamics or Monte Carlo. Hence the present treatment can handle the full spectrum of temporal events.

The present work is related to a recent model of protein zipping Weikl et al. 2003a ; Weikl et al. 2003b . Our fundamental units of protein structure are contact clusters. A contact cluster is a collection of contacts that is localized on a contact map, corresponding roughly to the main structural elements of the native structure. Examples of contact clusters are turns, $`\alpha `$-helices, $`\beta `$-strand pairings, and tertiary pairings of helices. A central quantity in our models is the effective contact order (ECO)Fiebig and Dill (1993); Dill et al. (1993). The ECO is the length of the loop that has to be closed in order to form a contact, given a set of previously formed contacts or contact clusters. The premise is that the formation of the nonlocal contact clusters requires the prior formation of other, more local, clusters.

Our model predicts average $`\mathrm{\Phi }`$-values for secondary structural elements that are in good agreement with the experimentally observed values for several two-state proteins. It shows that $`\mathrm{\Phi }`$-value distributions can be understood from loop-closure events that are defined by the native topology of a protein. The importance of topology for routes and $`\mathrm{\Phi }`$-values has also been previously noted by other groups Munoz and Eaton (1999); Alm and Baker (1999); Alm et al. (2002); Clementi et al. (2000); Vendruscolo et al. (2001).

To compute the dynamics, we use a master equation. Several previous studies of the folding kinetics of lattice heteropolymer models have also been based on master equation methods Leopold et al. (1992); Chan and Dill (1993); Cieplak et al. (1998); Ozkan et al. (2001, 2002, 2003); Schonbrun and Dill . These methods have the advantage that they require no ad hoc assumptions about what the transition state is. The transition state emerges in a direct physical way from the solution to the master equation. However, the lattice models are too simplified to treat specific amino acid sequences or specific protein structures. Lattice models focus on transitions between microstates, the individual chain conformations, since these are the fundamental units of structure in such models. Our present master equation describes transitions between macrostates, defined by the contact clusters of a given protein structure. In this way, the present model aims to make closer contact with experiments.

## II The Model

### Contact clusters

To compute the folding kinetics, we start with the native contact map, the matrix in which element $`(i,j)`$ equals 1 if the residues $`i`$ and $`j`$ are in contact, and equals 0 otherwise. Two residues are defined as being in contact if the distance between their $`C_\alpha `$ or $`C_\beta `$ atoms is less than 6 Angstroms.

Next, we divide the native contact map into contact clusters. Each contact cluster corresponds to a structural element of the protein. Two contacts $`(i,j)`$ and $`(k,l)`$ are defined as being in the same cluster if they are close together on the contact map, according to the distance criterion that $`|ik|+|jl|4`$. We define two types of clusters: local and nonlocal. Clusters are local if they contain at least one local contact $`(i,j)`$ having contact order $`CO=|ij|6`$. Local clusters include helices, turns, or $`\beta `$-hairpins, for example. A cluster is nonlocal if it has no local contacts; examples include $`\beta `$-strand pairings other than hairpins, and the tertiary interactions of helices. To qualify as nonlocal, a cluster must also have more than two contacts; isolated nonlocal contacts are not considered to be clusters. Similarly, we do not consider as contributing to clusters any ‘peripheral’ contacts $`(i,j)`$ with a minimum distance $`|ik|+|jl|=4`$ to the other contacts in the cluster. In general, typical contact maps have only a few isolated nonlocal or peripheral contacts. Fig. 1 shows examples of clusters, specifically for chymotrypsin inhibitor 2 (CI2) and the src SH3 domain. By our criteria, CI2 has 5 local clusters and 2 nonlocal clusters ($`\beta _2\beta _3`$ and $`\beta _1\beta _4`$), and the src SH3 domain has 6 local clusters and 2 nonlocal clusters (RT-$`\beta _4`$ and $`\beta _1\beta _5`$).

### States and free energies

We assume that each cluster is either formed or not; we neglect partial degrees of formation. Thus, for a protein with $`M`$ clusters, there are $`2^M`$ possible states that describe the progression to the native state. Each of these macrostates is characterized by a vector $`n=\{n_1,n_2,\mathrm{},n_M\}`$, where $`n_i=1`$ indicates that cluster $`i`$ is formed and $`n_i=0`$ indicates that cluster $`i`$ is not formed.

In our model, the free energy of the protein as a function of the state $`n`$ of cluster formation is given by:

$$F_n=\underset{i=1}{\overset{M}{}}n_i\left[c\mathrm{}_i(n)+f_i\right]$$
(1)
Each cluster $`i`$ that is formed ($`n_i=1`$) contributes to the free energy $`F_n`$ of the state $`n`$ with two terms: A state-dependent free energy of loop closure $`c\mathrm{}_i(n)`$ (‘initiation’ free energy), and a free energy $`f_i`$ for forming the cluster contacts (‘propagation’ free energy). Here, $`c`$ is a loop-closure parameter. The quantity $`\mathrm{}_i=\mathrm{}_i(n)`$ is the initiation ECO Weikl et al. 2003a  for cluster $`i`$. The initiation ECO of a cluster is the length of the smallest loop that must be closed in order to form that cluster from the other existing clusters. For a local cluster, the initiation ECO is the smallest CO among the contacts. For a nonlocal cluster, the initiation ECO depends on the presence of other clusters in the state $`n`$.

In general, the initiation ECO also depends on the sequence through which those clusters are formed. However, in order to apply the master equation formalism, we need a free energy and thus we require a definition of initiation ECO that is only a function of state. For this purpose, we use the following scheme. If only one nonlocal cluster is formed in a certain state, the initiation ECO of that cluster is the smallest ECO among the cluster contacts, given all the local clusters formed in that state. If multiple nonlocal clusters are present in a state, we consider all the possible sequences along which these clusters can form, and determine the one having the smallest sum of ECOs. For instance, for a state with two nonlocal clusters $`C_i`$ and $`C_j`$, there are two sequences: (1) $`C_i`$ $``$ $`C_j`$, and (2) $`C_j`$ $``$ $`C_i`$. The minimum ECOs for the clusters are determined sequentially: $`\mathrm{}_i^{(1)}`$ and $`\mathrm{}_j^{(1)}`$ along sequence (1), and $`\mathrm{}_i^{(2)}`$ and $`\mathrm{}_j^{(2)}`$ along sequence (2). If $`\mathrm{}_i^{(1)}+\mathrm{}_j^{(1)}`$ is smaller than $`\mathrm{}_i^{(2)}+\mathrm{}_j^{(2)}`$, the initiation ECOs $`\mathrm{}_i`$ and $`\mathrm{}_j`$ of the clusters $`i`$ and $`j`$ in the given state are taken to be $`\mathrm{}_i^{(1)}`$ and $`\mathrm{}_j^{(1)}`$. The initiation ECOs $`\mathrm{}_i`$ and $`\mathrm{}_j`$ are an estimate for the smallest loop lengths required to form the two clusters in the state.

In eq. (1), the free energy cost of the loops is estimated by a simple linear approximation in the loop length. This is not unreasonable since the range of relevant ECOs only spans roughly one order of magnitude, from about $`\mathrm{}3`$ to $`\mathrm{}30`$ or 40. In general, determining the free energy of a chain molecule with multiple constraints or contacts is a complicated and unsolved problem. For the simpler problem of hairpin-like loop closures, several estimates have been given in the literature (see, e.g., Chan and Dill (1990); Galzitskaya and Finkelstein (1999); Ivankov and Finkelstein (2001)).

In principle, this model could treat the detailed energetics of each folding route, if each of the $`M`$ clusters were characterized by its own free energy $`f_i`$. But here we consider a simpler version of the model. We assume that there are only two parameters for the free energy of formation: $`f_i=f_l`$ for propagating any local cluster, and $`f_i=f_{nl}`$ for propagating any nonlocal cluster. To obtain two-state folding and agreement with experimental $`\mathrm{\Phi }`$-values, we find that $`f_l`$ must be nonnegative and $`f_{nl}`$ must be negative. This is consistent with the experimental observation that local structures, such as helices or $`\beta `$-hairpins, are generally unstable in isolation. Similar in spirit, the diffusion-collision model of Karplus and Weaver assumes that microdomains, e.g. helices, are individually unstable Karplus and Weaver (1976, 1994). Thus, the rate-limiting barrier to folding in our model turns out to be the formation of mostly local structures needed to reduce the ECOs of nonlocal clusters. The driving force for overcoming this barrier is the favorable free energy $`f_{nl}`$ of assembling the nonlocal clusters.

The predicted free energy landscape of the src SH3 domain is shown in Fig. 2, using the parameters $`f_l=0`$ and $`c=0.5`$ $`k_BT`$, where $`k_BT`$ is Boltzmann’s constant $`\times `$ temperature. The value of $`f_{nl}`$ is chosen so that the equilibrium probability that the two nonlocal clusters RT-$`\beta _4`$ and $`\beta _1\beta _5`$ are both folded (‘native state’) is 0.9, which gives $`f_{nl}=6.6`$ $`k_BT`$ for src SH3. With these parameter settings, we obtain a good agreement with average experimental $`\mathrm{\Phi }`$-values for the src SH3 domain and other two-state folders (see below). For clarity, we show in the figure only a reduced set of states based on the 5 major clusters $`RT`$, $`\beta _2\beta _3`$, $`\beta _3\beta _4`$, RT-$`\beta _4`$, and $`\beta _1\beta _5`$. The three small clusters T, DT, and H have negligible effects on the folding kinetics and on the $`\mathrm{\Phi }`$-values. Only states differing by the formation of a single cluster are kinetically connected. The uphill steps in this model either are steps in which a local cluster is formed, or steps involving high ECOs. The downhill steps are steps in which a nonlocal cluster is formed with a low ECO, or steps in which a local cluster significantly reduces the ECOs of previously formed nonlocal clusters. The model predicts two main folding routes. Along the upper route (E) $`\beta _1\beta _5`$ folds after (D) RT-$`\beta _4`$; along the lower route, they form in the opposite order. Along these routes, the barriers (highest free energies states) are the states in which two clusters are formed: BD and BC for the upper route, and AC for the lower route.

### Master equation

In this section, we describe the folding dynamics. We use the master equation,

$$\frac{\text{d}P_n(t)}{\text{d}t}=\underset{mn}{}\left[w_{nm}P_m(t)w_{mn}P_n(t)\right],$$
(2)
which gives the time evolution of the probability $`P_n(t)`$ that the protein is in state $`n`$ at time $`t`$. Here, $`w_{nm}`$ is the transition rate from state $`m`$ to $`n`$. The master equation can be written in matrix form

$$\frac{\text{d}𝑷(t)}{\text{d}t}=𝑾𝑷(t)$$
(3)
where $`𝑷(t)`$ is the vector with elements $`P_n(t)`$, and the matrix elements of $`𝑾`$ are given by

$$W_{nm}=w_{nm}\text{for}nm;W_{nn}=\underset{mn}{}w_{mn}.$$
(4)
The transition rates are given in terms of the free energies by

$$w_{nm}=\frac{\delta _{|nm|,1}}{t_o}\left[1+\mathrm{exp}\left(\frac{F_nF_m}{k_BT}\right)\right]^1$$
(5)
where $`t_o`$ is a reference time scale. The only transitions that are assigned to have nonzero rates $`w_{nm}`$ are those incremental steps that change the state $`n`$ by a single cluster unit. This is enforced by the term $`\delta _{|nm|,1}`$ in eq. (5) where the Kronecker $`\delta _{i,j}`$ is one for $`i=j`$ and zero otherwise. The condition $`|nm|=1`$ is only satisfied by pairs of states $`n=\{n_1,\mathrm{},n_M\}`$ and $`m=\{m_1,\mathrm{},m_M\}`$ with $`n_km_k`$ for a single cluster $`k`$, and with $`n_k=m_k`$ for all other clusters. The transition rates (5) satisfy detailed balance, $`w_{nm}P_m^e=w_{mn}P_n^e`$ where $`P_n^e\mathrm{exp}[F_n/(k_BT)]`$ is the equilibrium weight for the state $`n`$. We have chosen here the ‘Glauber dynamics’ with $`w_{nm}(1+\mathrm{exp}[(F_nF_m)/(k_BT)])^1`$. Another standard choice satisfying detailed balance is the Metropolis dynamics, which should lead to equivalent results.

The detailed balance property of the transition rates implies that the eigenvalues of the matrix $`𝑾`$ are real. One of the eigenvalues is zero, corresponding to the equilibrium distribution, while all other eigenvalues are positive van Kampen (1992). The solution to the master equation is given by

$$𝑷(t)=\underset{\lambda }{}c_\lambda 𝒀_\lambda \mathrm{exp}[\lambda t]$$
(6)
where $`𝒀_\lambda `$ is the eigenvector corresponding to the eigenvalue $`\lambda `$, and the coefficients $`c_\lambda `$ are determined by the initial condition $`𝑷(t=0)`$. For $`t\mathrm{}`$, the probability distribution $`𝑷(t)`$ tends towards the equilibrium distribution $`𝑷^e𝒀_0`$ where $`𝒀_0`$ is the eigenvector with eigenvalue $`\lambda =0`$.

Solving the master equation gives a set of $`2^M`$ eigenvalues, each with its associated eigenvector. Each eigenvalue represents a relaxation rate. As initial conditions at $`t=0`$, we start from the state in which no clusters are formed. This corresponds to folding from high temperatures or high denaturant concentrations.

## III Results

### The cooperativity in two-state kinetics

The signature of two-state kinetics is the existence of one slow relaxation process (described by a single exponential), separated in time from $`2^M1`$ fast relaxations (a ‘burst’ phase). Fig. 3 shows the eigenvalue spectra for CI2 and the src SH3 domain, based on using the parameters $`c=0.5`$ $`k_BT`$, local cluster free energy $`f_l=0`$, and a nonlocal cluster free energy chosen so that the equilibrium ‘native’ population with all nonlocal clusters formed has probability 0.9. The latter condition leads to $`f_{nl}=7.9`$ $`k_BT`$ for CI2, and $`f_{nl}=6.6`$ $`k_BT`$ for src SH3. Fig. 4 shows the predicted folding dynamics for the src SH3 domain.

The spectra in Fig. 3 show that for these proteins, the eigenvalues do indeed separate into a slow single-exponential step and a burst phase, consistent with the experimental observation of two-state behavior. The slowest relaxation rate $`\lambda _1`$ is about one order of magnitude smaller than the other nonzero eigenvalues (see Fig. 3). At times $`t1/\lambda _1`$, the probability distribution (6) is well approximated by $`𝑷(t)c_0𝒀_0+c_1𝒀_1\mathrm{exp}[\lambda _1t]`$ where $`𝒀_0`$ is the eigenvector with eigenvalue 0, which characterizes the equilibrium state, and $`𝒀_1`$ is the eigenvector with eigenvalue $`\lambda _1`$.

The typical time evolution of the folding process predicted by the model is as follows. We have two time scales, $`tt_o`$ and $`t_F1/\lambda _1`$. Time $`t_o`$ is a characteristic of the burst phase in the model and $`t_F`$ is the single-exponential folding time. At the earliest times, $`t<t_o`$, single local clusters start to form: examples are the clusters A, B, and C of the src SH3 domain, see Fig. 2. As shown in Fig. 4, on this time scale, each cluster is only weakly populated, with a probability less than 10$`\%`$. Any structures having larger-scale organization – cluster pairs, triplets, etc. – have negligible populations. At intermediate times, $`t_Ftt_o`$, there is a crossover from the burst phase to the single-exponential folding process. During these intermediate times, cluster pairs (AC, BC, BD) begin to form. Fig. 2 shows that these pairwise clusters are the barrier events, i.e., they represent the conformational states of maximum free energy obtained during folding. Finally, on the longest time scale, $`tt_F`$, the pairwise and triplet clusters reach sufficiently high populations to assemble into multi-cluster complexes, proceeding downhill in free energy to the native structure.

What is the basis for the cooperativity of folding in our model, i.e. for the separation of time scales? First, the formation of local structures in our model reduces the loop-closure entropies for the formation of the nonlocal structures. Second, only the nonlocal structures have favorable propagation free energies $`f_i=f_{nl}<0`$. Hence, the formation of the nonlocal structures stabilizes the overall fold, and thus also the local structures. The barrier arises from the positive free energies in eq. (1) due to the formation of local structures and loops (see Fig. 2). Interestingly, if we set the free energies for local structure formation to be negative by several $`k_BT`$, we obtain fast multi-exponential downhill folding, without a barrier. Based on experiments and theory, such downhill folding has been recently postulated for the protein BBL Garcia-Mira et al. (2002).

To understand the cooperative folding in the model, it is instructive to turn off the loop-closure term in eq. (1) by setting $`c=0`$. Then all $`M`$ clusters are independent of each other. In that case, there is no cooperativity. It can be shown that the matrix $`𝑾`$ then has the eigenvalues $`\lambda =j/t_o`$ where $`j`$ is an integer between 0 and M, the number of clusters. Each of these eigenvalues has a population that is given by the binomial coefficient $`j!/[j!(Mj)!]`$. This gives a broad non-two-state spectrum. Hence, the separation of time scales – and the two-state cooperativity – arise in this model from the coupling of the clusters via the loop-closure term in eq. (1).

To see the magnitude of the barrier, note that the folding rate $`\lambda _1`$ is related to the height of the energy barrier on the folding landscape. For comparison, consider a mass-action model with three states D $``$ T $``$ N (denatured state, ‘transition’ state, native state) and transition rates as in eq. (5). The folding rate is given, to a very good approximation, by $`(1/2)t_o^1\mathrm{exp}[F_b^{}/(k_BT)]`$ for barrier energies $`F_b^{}=F_T^{}F_D^{}k_BT`$. The factor of 1/2 comes from the fact that a molecule in state T can jump both to D and N, with almost equal probability, since both sides of high-barrier transition states are steep downhills. Now, for the energy landscape of the src SH3 domain shown in Fig. 2, the minimum barrier has free energy 2.4 $`k_BT`$ for state BD. The corresponding barrier crossing rate of $`(1/2)t_o^1\mathrm{exp}[2.4]`$ is in good agreement with the folding rate $`\lambda _10.05/t_o`$ (see Fig. 3).

Experiments have been interpreted either as indicating that burst phases involve structure formation or that burst phases are processes of non-structured polymer collapse, depending on the protein and the experimental method Englander (2000); Callender et al. (1998); Gruebele et al. (1998); Eaton et al. (1998); Parker and Marqusee (2000); Ferguson and Fersht (2003). In our model, the burst phase is a process of structure formation. Non-structured collapse is beyond the scope, or resolution, of our model, because the model has only a single fully unstructured state – the state in which none of the clusters is formed. The burst phase in our model captures fast preequilibration events within the denatured state in response to initiating the folding conditions at $`t=0`$. In the model, this denatured state is an ensemble of macrostates on one side of the barrier in the energy landscape (see Fig. 2). It is reasonable to assume that such preequilibration events within the denatured state exist also for real proteins. However, whether these events can be detected as burst phases in experiments should depend on the initial conditions, experimental probes, etc.

During folding or unfolding, certain conformations will be populated transiently. If the populations of those conformations are always small, we call them ‘hidden intermediates’ Ozkan et al. (2002). The population of a hidden intermediate conformation rises to a maximum, $`P_{max}`$, then falls as the protein ultimately becomes fully folded. The term ‘hidden’ means that $`P_{max}`$ is always small enough that it does not contribute an additional kinetic phase; i.e., the folding kinetics is two-state. Here, we consider two quantities. (1) We compute $`P_{max}`$ for the transient states. For simplicity, we consider only the 5 major clusters $`RT`$, $`\beta _2\beta _3`$, $`\beta _3\beta _4`$, RT-$`\beta _4`$, and $`\beta _1\beta _5`$. (2) We look at the elements of the eigenvector $`𝐘_1`$, the eigenvector corresponding to the smallest eigenvalue $`\lambda _1`$. These elements show how the various conformations grow and decay with rate $`\lambda _1`$ as folding proceeds. Table 1 shows that the maximum population $`P_{max}`$ correlates well with the elements of $`𝐘_1`$. For a typical route of src SH3, Fig. 4 (bottom) illustrates the decay of the denatured state and hidden intermediates and the growth of the native state, all with rate $`\lambda _1`$.

### Average $`\mathrm{\Phi }`$-values for secondary structural elements

The effects of a mutation on the folding kinetics are often explored through experimental measurements of a $`\mathrm{\Phi }`$-value, which is defined as

$$\mathrm{\Phi }=\frac{k_BT\mathrm{ln}(k_f^{}/k_f)}{\mathrm{\Delta }G^{}\mathrm{\Delta }G}$$
(7)
where $`k_f`$ is the folding rate of the native protein and $`\mathrm{\Delta }G`$ is its stability, and $`k_f^{}`$ and $`\mathrm{\Delta }G^{}`$ are the corresponding quantities for the mutant protein.

Since the minimal structural units in our model are clusters of contacts, we do not calculate $`\mathrm{\Phi }`$-values for single-residue mutations. Rather, we consider whole helices and strands as units. To compare with experiments, we average the experimental $`\mathrm{\Phi }`$-values over all the residues composing a given secondary structural element.

To calculate average $`\mathrm{\Phi }`$-values for secondary structures, we consider ‘mutations’ that change the free energy $`f_i`$ of a contact cluster according to

$$\mathrm{\Delta }f_i(j)=x_{ji}ϵ$$
(8)
where $`x_{ji}`$ is the fraction of residues of the secondary structural element $`j`$ that are involved in contacts of the cluster $`i`$, and $`ϵ`$ is a small energy. For example, if the secondary structural element $`j`$ contains $`m_1`$ residues, and $`m_2m_1`$ of these residues appear in contacts of the cluster $`i`$, we have $`x_{ji}=m_2/m_1`$. Note that $`0x_{ji}1`$, where the value $`x_{ji}=1`$ is obtained if the whole secondary structural element $`j`$ has contacts in cluster $`i`$. Thus the $`\mathrm{\Phi }`$-value for the secondary structural element $`j`$ is given by eq. (7) with

$$\mathrm{ln}(k_f^{}/k_f)=\mathrm{ln}(\lambda _1^{}/\lambda _1)$$
(9)
where $`\lambda _1^{}`$ is the smallest nonzero eigenvalue of the mutant with cluster free energies $`f_if_i+\mathrm{\Delta }f_i(j)`$, and

$$\mathrm{\Delta }G^{}\mathrm{\Delta }G=\underset{i}{}\mathrm{\Delta }f_i(j)$$
(10)
For $`ϵk_BT`$, we find that the calculated $`\mathrm{\Phi }`$-values are nearly independent of $`ϵ`$. We choose here $`ϵ=0.01`$ $`k_BT`$.

Predicted $`\mathrm{\Phi }`$-values are compared with experiments in Fig. 5. The theoretical $`\mathrm{\Phi }`$-values were calculated with the same parameters for all four proteins (see figure caption). The predicted values agree well with the experimental values. This comparison indicates that the folding kinetics of these proteins is dominated by generic features of the fold topology, rather than by the specific energetic details – i.e., which residues form contacts, how much hydrogen bonds or hydrophobic interactions are worth, the details of sidechain packing, etc. In the case of protein G (see Fig. 6), the experimental $`\mathrm{\Phi }`$-value distribution is largely reproduced by making the additional assumption that the $`\alpha `$-helix cluster has a free energy $`f_i=2.0`$ $`k_BT`$, rather than the value $`f_i=0`$ $`k_BT`$ that we have otherwise used for local clusters (see Fig. 5). However, even without changing this parameter, the $`\mathrm{\Phi }`$-value distribution reflects the features of the experimental distribution that the $`\mathrm{\Phi }`$-values for the strands $`\beta _3`$ and $`\beta _4`$ are larger than those for $`\beta _1`$ and $`\beta _2`$.

## IV Conclusions

We have developed a simple model of the folding kinetics of two-state proteins. The model aims to predict the folding rates of the fast and slow processes, the folding routes, and $`\mathrm{\Phi }`$-values for a protein, if the native structure is given. The dominant folding routes are found to be those having small ECOs, i.e., steps that involve only small ‘loop closures’. The model parameters include: $`c`$, an intrinsic free energy for loop closure; $`f_l`$, the free energy for propagating contacts in local structures; and $`f_{nl}`$, the free energy for propagating nonlocal contacts. The model predicts that the barrier to two-state folding is the formation of local structural elements like helices and hairpins, and that the steps involving their assembly into larger and more native-like structure are downhill in free energy.


format: text

text:
# First-principle molecular dynamics with ultrasoft pseudopotentials: parallel implementation and application to extended bio-inorganic systems

## I Introduction

There is increasing interest in studying the electronic structure of complex biological molecules. This is an essential step to understand e.g. enzymatic and/or biomimetic catalysis. Modeling bio-catalytic systems is however very challenging, because a proper description of the active site needs the inclusion of a large number of atoms (from several tens to a few hundreds) treated at a high level of quantum chemical theory. In this respect a good compromise in terms of accuracy and computational cost is provided by Density-Functional Theory (DFT)DFT , whose use to model the electronic structure of protein active sites is becoming increasingly popular. In combination with Car-Parrinello (CP)CP  first-principle molecular dynamics (MD), DFT allows us to optimize molecular structures, study dynamical and finite-temperature properties, and model reaction paths.

In most standard implementations, the CP method employs a plane-wave (PW) basis set. An advantage of PWs is that they do not depend on atomic positions and are free of basis-set superposition errors. Total energies and forces on the atoms can be calculated using computationally efficient Fast Fourier transform (FFT) techniques. Finally, the convergence of a calculation can be controlled in a simple way, since it depends only upon the number of PWs included in the expansion of the electron density. The dimension of a PW basis set is controlled by a cutoff in the kinetic energy of the PWs, which is usually measured in Ry units. A disadvantage of PWs is their extremely slow convergence in describing core states. To deal with this difficulty, one usually employs norm-conserving (NC) pseudopotentials (PPs) NCPP  to model the interaction of the valence electrons with the ionic core (nucleus + core electrons). Parallel implementations of PW calculations based on NC PPs are well documented in the literature (see e.g. ParaFFT ; MarxHutter ; Cavazzoni ).

When using NC PPs, very large PW basis sets are needed to accurately represent the contracted $`p`$ orbitals of the first-row elements O, N, F, and the $`3d`$ orbitals of the transition metal block. These orbitals belong to elements for which the NC PPs are “hard”, typically requiring cutoffs of more than 70 Ry in order to yield sufficiently converged results. As a comparison, calculations on elements like Al, Si, P, for which the corresponding NC PPs are “soft”, are usually well converged with a cutoff of 20 Ry or less. A consequence of the delocalized nature of the PWs is that the presence of a single hard PP in a system requires the use of a correspondingly high cutoff for all the other PPs. This difficulty is particularly serious for metallorganic systems containing one or more transition metal centers. The high cutoff required for such atoms translates into a very large number of PWs, which in turn implies long execution times and large memory requirements.

An approach that drastically reduces the PW cutoff was proposed by Vanderbilt Van , who introduced “ultrasoft” (US) PPs. In this approach the normalized charge density is the sum of two terms, a soft part represented in terms of smooth orbitals, and a hard part which is treated as an augmented charge. A closely related approach, the “Projector-augmented wave” (PAW) method introduced by Blöchl Blo , is an all-electron rather than a PP electronic structure method. The PAW approach provides a simple and effective algorithm for reconstructing all-electron orbitals from pseudo-orbitalsReconstruct . An efficient serial implementation of the CP scheme with US PPs is described in Ref. CarPasquarello .

In this paper we present in detail a parallel implementation of the CP scheme using US PPs. We also provide further details on the procedures used in serial implementation. We compare the relative efficiency of US and NC PPs in realistic calculations for large molecules, performed on parallel machines. Our test molecules are a reduced and an extended model of the active site of myoglobin, containing the iron-porphyrin motif. We focus on metalloporphyrin systems because they are representative of a large class of biomolecules which can be modeled efficiently with US PPs. Our benchmarks indicate that US calculations are at least 2-3 times less expensive than NC calculations of comparable accuracy.

The use of a PW basis set implies that Periodic Boundary Conditions (PBC) are imposed, i.e. an isolated molecule has to be placed into a periodically repeated box (a “supercell”). The supercell must be large enough to ensure that the total potential is vanishingly small at the box boundary, thus minimizing spurious interactions between periodic replicas. For neutral systems with small dipole/quadrupole moments, supercells of reasonable size can be safely used. For charged molecules, or molecules with large dipole/quadrupole moments, however, the error induced by PBC may be rather large unless exceedingly large supercells are used.

We follow here a technique introduced by Makov and Payne (MP)MakovPayne  to eliminate the spurious electrostatic interactions in the latter case. The MP technique is approximate because it is not self-consistent and takes into account only moments up to quadrupole. To check the accuracy of the MP technique, we compare CP calculations on highly charged manganese porphyrins, performed using PWs and US PPs, with calculations on the same systems using localized basis sets which do not require PBC. The comparison shows that the MP correction yields results that to all practical effects are indistinguishable from results obtained without PBC.

The paper is organized as follows. In Sec.2, we recall the main aspects of US-PP implementation in the serial case. In Sec.3, we describe our parallel implementation. In Sec.4, we compare the computer performances of US and NC PPs for a reduced and an extended model of the myoglobin active site. In Sec.5, we compare CP calculations with localized basis-set calculations not requiring PBC. The test systems are highly-charged isomeric meso-substituted manganese porphyrins. Sec.6 contains our conclusions.

## II Plane-Wave Ultrasoft Pseudopotential Implementation

### II.1 Kohn-Sham equations with Ultrasoft Pseudopotentials

The implementation of CP molecular dynamics within a US-PP framework is described in Ref.CarPasquarello . Here we briefly remind the main formulas, using the same notation of Ref.CarPasquarello .

The total energy of a system of $`N_v`$ valence electrons, having one-electron Kohn-Sham (KS) orbitals $`\varphi _i`$, is given by

$`E_{\mathrm{tot}}[\{\varphi _i\},\{𝐑_I\}]`$ $`=`$ $`{\displaystyle \underset{i}{}}\varphi _i|{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_{\mathrm{NL}}|\varphi _i+E_\mathrm{H}[n]+E_{\mathrm{xc}}[n]`$ (1)
$`+{\displaystyle 𝑑𝐫V_{\mathrm{loc}}^{\mathrm{ion}}(𝐫)n(𝐫)}+U(\{𝐑_I\}),`$
where $`n(𝐫)`$ is the electron density, $`E_\mathrm{H}[n]`$ is the Hartree energy:

$`E_\mathrm{H}[n]={\displaystyle \frac{e^2}{2}}{\displaystyle 𝑑𝐫𝑑𝐫^{}\frac{n(𝐫)n(𝐫^{})}{|𝐫𝐫^{}|}},`$ (2)
$`E_{\mathrm{xc}}[n]`$ is the exchange and correlation energy, $`U(\{𝐑_I\})`$ is the ion-ion interaction energy, and $`𝐑_I`$ indicate atomic positions. In the following, potentials have energy dimensions. The PP is composed of a local part $`V_{\mathrm{loc}}^{\mathrm{ion}}`$, given by a sum of atom-centered radial potentials:

$`V_{\mathrm{loc}}^{\mathrm{ion}}(𝐫)={\displaystyle \underset{I}{}}V_{\mathrm{loc}}^I(|𝐫𝐑_I|)`$ (3)
and a nonlocal part $`V_{\mathrm{NL}}`$:

$`V_{\mathrm{NL}}={\displaystyle \underset{nm,I}{}}D_{nm}^{(0)}|\beta _n^I\beta _m^I|,`$ (4)
where the functions $`\beta _n^I`$ and the coefficients $`D_{nm}^{(0)}`$ characterize the PP and are specific for each atomic species. For simplicity, we consider a single atomic species only in what follows. The $`\beta _n^I`$ functions, centered at site $`𝐑_I`$, depend on the ionic positions through

$`\beta _n^I(𝐫)=\beta _n(𝐫𝐑_I).`$ (5)
$`\beta _n`$ is an angular momentum eigenfunction in the angular variables, times a radial function which vanishes outside the core region; the indices $`n`$ and $`m`$ in Eq. (4) run over the total number $`N_\beta `$ of these functions.

The electron density in Eq. (1) is given by

$`n(𝐫)={\displaystyle \underset{i}{}}\left[|\varphi _i(𝐫)|^2+{\displaystyle \underset{nm,I}{}}Q_{nm}^I(𝐫)\varphi _i|\beta _n^I\beta _m^I|\varphi _i\right],`$ (6)
where the sum runs over occupied KS orbitals. The augmentation functions $`Q_{nm}^I(𝐫)=Q_{nm}(𝐫𝐑_I)`$ are localized in the core. The ultrasoft PP is fully determined by the quantities $`V_{\mathrm{loc}}^I(r),D_{nm}^{(0)},Q_{nm}(𝐫),`$ and $`\beta _n(𝐫)`$. The functions $`Q_{nm}(𝐫)`$ are defined in terms of atomic orbitals as: $`Q_{nm}(𝐫)=\psi _n^{ae}(𝐫)\psi _m^{ae}(𝐫)\psi _n^{ps}(𝐫)\psi _m^{ps}(𝐫)`$ where $`\psi ^{ae}`$ are atomic one-electron orbitals (not necessarily bound), and $`\psi ^{ps}`$ are the corresponding pseudo-orbitals. The $`Q_{nm}(𝐫)`$ are pseudized as described in Ref.CarPasquarello . This enables us to treat the $`Q_{nm}(𝐫)`$ with Fourier transform techniques.

The KS orbitals obey generalized orthonormality conditions

$`\varphi _i|S(\{𝐑_I\})|\varphi _j=\delta _{ij},`$ (7)
where $`S`$ is a Hermitian overlap operator given by

$`S=1+{\displaystyle \underset{nm,I}{}}q_{nm}|\beta _n^I\beta _m^I|,`$ (8)
and

$`q_{nm}={\displaystyle 𝑑𝐫Q_{nm}(𝐫)}.`$ (9)
The orthonormality condition (7) is consistent with the conservation of the charge $`𝑑𝐫n(𝐫)=N_v`$. Note that the overlap operator $`S`$ depends on ionic positions through the $`|\beta _n^I`$.

The ground-state orbitals $`\varphi _i`$ minimize the total energy (1) subject to the constraints (7),

$`{\displaystyle \frac{\delta E_{\mathrm{tot}}}{\delta \varphi _i^{}(𝐫)}}=ϵ_iS\varphi _i(𝐫),`$ (10)
where the $`ϵ_i`$ are Lagrange multipliers. Eq.(11) is the KS equation:

$`H|\varphi _i=ϵ_iS|\varphi _i`$ (11)
where

$`H={\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_{\mathrm{eff}}+{\displaystyle \underset{nm,I}{}}D_{nm}^I|\beta _n^I\beta _m^I|.`$ (12)
$`V_{\mathrm{eff}}`$ is a screened effective local potential,

$`V_{\mathrm{eff}}(𝐫)=V_{\mathrm{loc}}^{\mathrm{ion}}(𝐫)+V_\mathrm{H}(𝐫)+\mu _{\mathrm{xc}}(𝐫).`$ (13)
$`\mu _{\mathrm{xc}}(𝐫)`$ is the exchange-correlation potential:

$`\mu _{\mathrm{xc}}(𝐫)={\displaystyle \frac{\delta E_{\mathrm{xc}}[n]}{\delta n(𝐫)}},`$ (14)
and $`V_\mathrm{H}(𝐫)`$ is the Hartree potential:

$`V_\mathrm{H}(𝐫)=e^2{\displaystyle 𝑑𝐫^{}\frac{n(𝐫^{})}{|𝐫𝐫^{}|}}.`$ (15)
The “screened” coefficients $`D_{nm}^I`$ appearing in Eq. (12) are defined as:

$`D_{nm}^I=D_{nm}^{(0)}+{\displaystyle 𝑑𝐫V_{\mathrm{eff}}(𝐫)Q_{nm}^I(𝐫)}.`$ (16)
They depend on the KS orbitals through $`V_{\mathrm{eff}}`$, Eq. (13), and the charge density, Eq. (6).

### II.2 Molecular dynamics with Ultrasoft Pseudopotentials

In the CP approach,CP  the electronic orbitals and the ionic coordinates evolve according to a classical Lagrangian

$`=\mu {\displaystyle \underset{i}{}}{\displaystyle 𝑑𝐫|\dot{\varphi }_i(𝐫)|^2}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{I}{}}M_I\dot{𝐑}_I^2E_{\mathrm{tot}}(\{\varphi _i\},\{𝐑_I\}),`$ (17)
subject to a set of constraints

$`𝒩_{ij}(\{\varphi _i\},\{𝐑_I\})=\varphi _i|S|\varphi _j\delta _{ij}=0.`$ (18)
Here $`\mu `$ is a fictitious mass parameter for the electronic degrees of freedom, $`M_I`$ are the masses of the atoms, and $`E_{\mathrm{tot}}`$ and $`S`$ are given in Eqs. (1) and (8), respectively. The holonomic orthonormality constraints (18) do not cause energy dissipation in an MD run. They may be incorporated in the Euler equations of motion by introducing Lagrange multipliers $`\mathrm{\Lambda }_{ij}`$:

$`\mu \ddot{\varphi }_i`$ $`=`$ $`{\displaystyle \frac{\delta E_{\mathrm{tot}}}{\delta \varphi _i^{}}}+{\displaystyle \underset{j}{}}\mathrm{\Lambda }_{ij}S\varphi _j,`$ (19)
$`𝐅_I=M_I\ddot{𝐑}_I`$ $`=`$ $`{\displaystyle \frac{E_{\mathrm{tot}}}{𝐑_I}}+{\displaystyle \underset{ij}{}}\mathrm{\Lambda }_{ij}\varphi _i|{\displaystyle \frac{S}{𝐑_I}}|\varphi _j.`$ (20)
At equilibrium, Eq. (19) reduces to the electronic KS equations (10) or (11). A unitary rotation brings the $`\mathrm{\Lambda }`$ matrix into diagonal form: $`\mathrm{\Lambda }_{ij}=ϵ_i\delta _{ij}`$. The equilibrium for the ions is achieved when the ionic forces $`𝐅_I`$ in Eq. (20) vanish.

In deriving explicit expressions for the forces, Eq. (20), one should keep in mind that the electron density also depends on $`𝐑_I`$ through $`Q_{nm}^I`$ and $`\beta _n^I`$. Introducing the quantities

$`\rho _{nm}^I={\displaystyle \underset{i}{}}\varphi _i|\beta _n^I\beta _m^I|\varphi _i,`$ (21)
and

$`\omega _{nm}^I={\displaystyle \underset{ij}{}}\mathrm{\Lambda }_{ij}\varphi _j|\beta _n^I\beta _m^I|\varphi _i,`$ (22)
we arrive at the expression

$`𝐅_I`$ $`=`$ $`{\displaystyle \frac{U}{𝐑_I}}{\displaystyle 𝑑𝐫\frac{V_{\mathrm{loc}}^{\mathrm{ion}}}{𝐑_I}n(𝐫)}{\displaystyle 𝑑𝐫V_{\mathrm{eff}}(𝐫)\underset{nm}{}\frac{Q_{nm}^I(𝐫)}{𝐑_I}\rho _{nm}^I}`$ (23)
$`{\displaystyle \underset{nm}{}}D_{nm}^I{\displaystyle \frac{\rho _{nm}^I}{𝐑_I}}+{\displaystyle \underset{nm}{}}q_{nm}{\displaystyle \frac{\omega _{nm}^I}{𝐑_I}},`$
where $`D_{nm}^I`$ and $`V_{\mathrm{eff}}`$ have been defined in Eqs. (16) and (13), respectively. The last term of Eq. (23) gives the constraint contribution to the forces. Since the PW basis set does not depend on atomic positions, Pulay-type corrections Pulay  do not appear in the expression for the forces.

### II.3 Discretization of the equation of motion and orthonormality constraints

The equations of motion (19) and (20) are usually discretized using the Verlet or the velocity-Verlet algorithms. The following discussion, including the treatment of the $`𝐑_I`$-dependence of the orthonormality constraints, applies to the Verlet algorithm when using the Fourier acceleration scheme of Ref.prec . In this framework the fictitious electron mass is represented by an operator $`\mathrm{\Theta }`$, whose matrix elements between PWs are given by

$$\mathrm{\Theta }_{𝐆,𝐆^{}}=\text{max}(\mu ,\mu \frac{\mathrm{}^2G^2}{2mE_c})\delta _{𝐆,𝐆^{}},$$
(24)
where $`𝐆,𝐆^{}`$ are the wave vector of PWs, $`E_c`$ is a cutoff (typically a few Ry) which defines the threshold for Fourier acceleration. The fictitious electron mass depends on $`G`$ as the kinetic energy for large $`G`$, it is constant for small $`G`$. This scheme allows us to use larger time steps with negligible computational overhead.

The electronic orbitals at time $`t+\mathrm{\Delta }t`$ are given by:

$`\varphi _i(t+\mathrm{\Delta }t)=2\varphi _i(t)\varphi _i(t\mathrm{\Delta }t)(\mathrm{\Delta }t)^2\mathrm{\Theta }^1\left[{\displaystyle \frac{\delta E_{\mathrm{tot}}}{\delta \varphi _i^{}}}{\displaystyle \underset{j}{}}\mathrm{\Lambda }_{ij}(t+\mathrm{\Delta }t)S(t)\varphi _j(t)\right],`$ (25)
where $`\mathrm{\Delta }t`$ is the time step, and $`S(t)`$ indicates the operator $`S`$ evaluated for ionic positions $`R_I(t)`$. Similarly the ionic coordinates at time $`t+\mathrm{\Delta }t`$ are given by:

$`𝐑_I(t+\mathrm{\Delta }t)`$ $`=`$ $`2𝐑_I(t)𝐑_I(t\mathrm{\Delta }t)`$ (26)
$`{\displaystyle \frac{(\mathrm{\Delta }t)^2}{M_I}}\left[{\displaystyle \frac{E_{\mathrm{tot}}}{𝐑_I}}{\displaystyle \underset{ij}{}}\mathrm{\Lambda }_{ij}(t+\mathrm{\Delta }t)\varphi _i(t)|{\displaystyle \frac{S(t)}{𝐑_I}}|\varphi _j(t)\right].`$
The orthonormality conditions must be imposed at each time-step:

$$\varphi _i(t+\mathrm{\Delta }t)|S(t+\mathrm{\Delta }t)|\varphi _j(t+\mathrm{\Delta }t)=\delta _{ij},$$
(27)
leading to the following matrix equation:

$`A+\lambda B+B^{}\lambda ^{}+\lambda C\lambda ^{}=1`$ (28)
where the unknown matrix $`\lambda `$ is related to the matrix of Lagrange multipliers $`\mathrm{\Lambda }`$ at time $`t+\mathrm{\Delta }t`$ via $`\lambda =(\mathrm{\Delta }t)^2\mathrm{\Lambda }^{}(t+\mathrm{\Delta }t)`$. In Eq.(28) the dagger indicates Hermitian conjugate ($`\lambda =\lambda ^{}`$). The matrices $`A`$, $`B`$, and $`C`$ are given by:

$`A_{ij}`$ $`=`$ $`\overline{\varphi }_i|S(t+\mathrm{\Delta }t)|\overline{\varphi }_j,`$
$`B_{ij}`$ $`=`$ $`\mathrm{\Theta }^1S(t)\varphi _i(t)|S(t+\mathrm{\Delta }t)|\overline{\varphi }_j,`$
$`C_{ij}`$ $`=`$ $`\mathrm{\Theta }^1S(t)\varphi _i(t)|S(t+\mathrm{\Delta }t)|\mathrm{\Theta }^1S(t)\varphi _j(t),`$ (29)
with

$`\overline{\varphi }_i`$ $`=`$ $`2\varphi _i(t)\varphi _i(t\mathrm{\Delta }t)(\mathrm{\Delta }t)^2\mathrm{\Theta }^1{\displaystyle \frac{\delta E_{\mathrm{tot}}(t)}{\delta \varphi _i^{}}}.`$ (30)
The solution of Eq. (28) in the ultrasoft PP case is not obvious, because Eq. (26) is not a closed expression for $`𝐑_I(t+\mathrm{\Delta }t)`$. The problem is that $`\mathrm{\Lambda }(t+\mathrm{\Delta }t)`$ appearing in Eq. (26) depends implicitly on $`𝐑_I(t+\mathrm{\Delta }t)`$ through $`S(t+\mathrm{\Delta }t)`$. Consequently, it is in principle necessary to solve iteratively for $`𝐑_I(t+\mathrm{\Delta }t)`$ in Eq. (26).

A simple solution to this problem is given in Ref. CarPasquarello . $`\mathrm{\Lambda }(t+\mathrm{\Delta }t)`$ is extrapolated using two previous values:

$$\mathrm{\Lambda }_{ij}^{(0)}(t+\mathrm{\Delta }t)=2\mathrm{\Lambda }_{ij}(t)\mathrm{\Lambda }_{ij}(t\mathrm{\Delta }t).$$
(31)
Eq. (26) is used to find $`𝐑_I^{(0)}(t+\mathrm{\Delta }t)`$, which is correct to $`O(\mathrm{\Delta }t^4)`$. From $`𝐑_I^{(0)}(t+\mathrm{\Delta }t)`$ we can obtain a new set $`\mathrm{\Lambda }_{ij}^{(1)}(t+\mathrm{\Delta }t)`$ and repeat the procedure until convergence is achieved. It turns out that in most practical applications the procedure converges at the very first iteration. Thus, the operations described above are generally executed only once per time step.

The solution of Eq. (28) is found using a modified version CarPasquarello ; Cavazzoni  of the iterative procedure of Ref. cp89 . The matrix $`B`$ is decomposed into hermitian $`(B_h)`$ and antihermitian $`(B_a)`$ parts,

$`B=B_h+B_a,`$ (32)
and the solution is obtained by iteration:

$`\lambda ^{(n+1)}B_h+B_h\lambda ^{(n+1)}=1A\lambda ^{(n)}B_aB_a^{}\lambda ^{(n)}\lambda ^{(n)}C\lambda ^{(n)}.`$ (33)
The initial guess $`\lambda ^{(0)}`$ can be obtained from

$`\lambda ^{(0)}B_h+B_h\lambda ^{(0)}=1A.`$ (34)
Here the $`B_a`$\- and $`C`$-dependent terms are neglected because they are of higher order in $`\mathrm{\Delta }t`$ ($`B_a`$ vanishes for vanishing $`\mathrm{\Delta }t`$). Eqs. (34) and (33) have the same structure:

$`\lambda B_h+B_h\lambda =X`$ (35)
where $`X`$ a Hermitian matrix. Eq. (35) can be solved exactly by finding the unitary matrix $`U`$ that diagonalizes $`B_h`$: $`U^{}B_hU=D`$, where $`D_{ij}=d_i\delta _{ij}`$. The solution is obtained from

$`(U^{}\lambda U)_{ij}=(U^{}XU)_{ij}/(d_i+d_j).`$ (36)
When $`X=1A`$ Eq.(36) yields the starting $`\lambda ^{(0)}`$, while $`\lambda ^{(n+1)}`$ is obtained from $`\lambda ^{(n)}`$ by solving Eq. (36) with $`X`$ given by Eq. (33). This iterative procedure usually converges in ten steps or less.

### II.4 Ultrasoft pseudopotential implementation in the serial case

#### II.4.1 Plane-wave expansion

Let $`\{𝐑\}`$ be the translation vectors of the periodically repeated supercell. The corresponding reciprocal lattice vectors $`\{𝐆\}`$ obey the conditions $`𝐑_i𝐆_j=2\pi n`$, with $`n`$ an integer number.

The KS orbitals can be expanded in a PW basis up to a kinetic energy cutoff $`E_\mathrm{c}^{\mathrm{wf}}`$:

$`\varphi _{j,𝐤}(𝐫)={\displaystyle \frac{1}{\sqrt{\mathrm{\Omega }}}}{\displaystyle \underset{𝐆\{𝐆_\mathrm{c}^{\mathrm{wf}}\}}{}}\varphi _{j,𝐤}(𝐆)e^{i(𝐤+𝐆)𝐫},`$ (37)
where $`\mathrm{\Omega }`$ is the volume of the cell, $`\{𝐆_\mathrm{c}^{\mathrm{wf}}\}`$ is the set of $`𝐆`$ vectors satisfying the condition

$`{\displaystyle \frac{\mathrm{}^2}{2m}}|𝐤+𝐆|^2<E_\mathrm{c}^{\mathrm{wf}},`$ (38)
and k is the Bloch vector of the electronic states. In crystals, one must use a grid of k-points dense enough to sample the Brillouin zone (the unit cell of the reciprocal lattice). In molecules, liquids and in general if the simulation cell is large enough, the Brillouin zone can be sampled using only the $`𝐤=0`$ ($`\mathrm{\Gamma }`$) point. An advantage of this choice is that the orbitals can be taken to be real in r-space. In the following we will drop the k vector index. Functions in real space and their Fourier transforms will be denoted by the symbols, if this does not originate ambiguity.

The $`\varphi _j(𝐆)`$’s are the electronic variables. The calculation of $`H\varphi _j`$ and of the forces acting on the ions are the basic ingredients of the computation. Scalar products $`\varphi _j|\beta _n^I`$ and their spatial derivatives are typically evaluated in G-space. An important advantage of working in G-space is that atom-centered functions like $`\beta _n^I`$ and $`Q_{nm}^I`$ are easily evaluated at any atomic position, for example:

$`\beta _n^I(𝐆)=\beta _n(𝐆)e^{i𝐆𝐑_I}.`$ (39)
Thus:

$`\varphi _j|\beta _n^I={\displaystyle \underset{𝐆\{𝐆_\mathrm{c}^{\mathrm{wf}}\}}{}}\varphi _j^{}(𝐆)\beta _n(𝐆)e^{i𝐆𝐑_I},`$ (40)
and

$`\varphi _j|{\displaystyle \frac{\beta _n^I}{𝐑_I}}=i{\displaystyle \underset{𝐆\{𝐆_\mathrm{c}^{\mathrm{wf}}\}}{}}𝐆\varphi _j^{}(𝐆)\beta _n(𝐆)e^{i𝐆𝐑_I}.`$ (41)
The kinetic energy term is diagonal in G-space and is easily calculated:

$`\left(^2\varphi _j\right)(𝐆)=G^2\varphi _j(𝐆).`$ (42)
In summary, the kinetic and nonlocal PP terms in $`H\varphi _j`$ are calculated in G-space.

#### II.4.2 Dual space technique

The local potential term $`V_{\mathrm{eff}}\varphi _j`$ could be calculated in G-space, but it is more convenient to use a different (“dual space”) technique. The idea is to switch from G\- to r-space, back and forth, using FFT, and to perform the calculation in the space where it is more convenient. The KS orbitals are first Fourier-transformed to r-space; then, $`(V_{\mathrm{eff}}\varphi _j)(𝐫)=V_{\mathrm{eff}}(𝐫)\varphi _j(𝐫)`$ is calculated in r-space, where $`V_{\mathrm{eff}}`$ is diagonal; finally $`(V_{\mathrm{eff}}\varphi _j)(𝐫)`$ is Fourier-transformed back to $`(V_{\mathrm{eff}}\varphi _j)(𝐆)`$.

In order to use FFT, one discretizes the r-space by a uniform grid spanning the unit cell:

$$f(m_1,m_2,m_3)f(𝐫_{m_1,m_2,m_3}),𝐫_{m_1,m_2,m_3}=m_1\frac{𝐚_1}{N_1}+m_2\frac{𝐚_2}{N_2}+m_3\frac{𝐚_3}{N_3},$$
(43)
where $`𝐚_1,𝐚_2,𝐚_3`$ are lattice basis vectors, the integer index $`m_1`$ runs from $`0`$ to $`N_11`$, and similarly for $`m_2`$ and $`m_3`$. In the following we will assume for simplicity that $`N_1,N_2,N_3`$ are even numbers. The FFT maps a discrete periodic function in real space $`f(m_1,m_2,m_3)`$ into a discrete periodic function in reciprocal space $`\stackrel{~}{f}(n_1,n_2,n_3)`$ (where $`n_1`$ runs from $`0`$ to $`N_11`$, and similarly for $`n_2`$ and $`n_3`$), and vice versa.

The link between G-space components and FFT indices is:

$$\stackrel{~}{f}(n_1,n_2,n_3)f(𝐆_{n_1^{},n_2^{},n_3^{}}),𝐆_{n_1^{},n_2^{},n_3^{}}=n_1^{}𝐛_1+n_2^{}𝐛_2+n_3^{}𝐛_3$$
(44)
where $`n_1=n_1^{}`$ if $`n_1^{}0`$, $`n_1=n_1^{}+N_1`$ if $`n_1^{}<0`$, and similarly for $`n_2`$ and $`n_3`$. The FFT dimensions $`N_1,N_2,N_3`$ must be big enough to include all non negligible Fourier components of the function to be transformed: ideally the Fourier component corresponding to $`n_1^{}=N_1/2`$, and similar for $`n_2^{}`$ and $`n_3^{}`$, should vanish. In the following, we will refer to the set of indices $`n_1,n_2,n_3`$ and to the corresponding Fourier components as the “FFT grid”.

The soft part of the charge density: $`n_{\mathrm{soft}}(𝐫)=_j|\varphi _j(𝐫)|^2`$, contains Fourier components up to a kinetic energy cutoff $`E_\mathrm{c}^{\mathrm{soft}}=4E_\mathrm{c}^{\mathrm{wf}}`$. This is evident from the formula:

$`n_{\mathrm{soft}}(𝐆)={\displaystyle \underset{𝐆^{}\{𝐆_\mathrm{c}^{\mathrm{wf}}\}}{}}{\displaystyle \underset{j}{}}\varphi _j^{}(𝐆𝐆^{})\varphi _j(𝐆^{}).`$ (45)
In the case of NC PPs, the entire charge density is given by $`n_{\mathrm{soft}}(𝐫)`$.

$`V_{\mathrm{eff}}`$ should be expanded up to the same $`E_\mathrm{c}^{\mathrm{soft}}`$ cutoff since all the Fourier components of $`V_{\mathrm{eff}}\varphi _j`$ up to $`E_\mathrm{c}^{\mathrm{wf}}`$ are required. Let us call $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$ the set of G-vectors such that

$`{\displaystyle \frac{\mathrm{}}{2m}}G^2<E_\mathrm{c}^{\mathrm{soft}}.`$ (46)
The soft part of the charge density is conveniently calculated in r-space, by Fourier-transforming $`\varphi _j(𝐆)`$ into $`\varphi _j(𝐫)`$ and summing over the occupied states.

The exchange-correlation potential $`\mu _{\mathrm{xc}}(𝐫)`$, Eq. (14), is a function of the local charge density and – for gradient-corrected functionals – of its gradient at point r:

$`\mu _{\mathrm{xc}}(𝐫)=V_{\mathrm{xc}}(n(𝐫),|n(𝐫)|).`$ (47)
The gradient $`n(𝐫)`$ is conveniently calculated from the charge density in G-space, using $`(n)(𝐆)=i𝐆n(𝐆)`$. The Hartree potential $`V_\mathrm{H}(𝐫)`$, Eq. (15), is also conveniently calculated in G-space:

$`V_\mathrm{H}(𝐆)={\displaystyle \frac{4\pi }{\mathrm{\Omega }}}{\displaystyle \frac{n(𝐆)^{}}{G^2}}.`$ (48)
Thus in the NC-PP case, a single FFT grid, large enough to accommodate the $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$ set, can be used for orbitals, charge density, and potential.

The use of FFT is mathematically equivalent to a pure G-space description (we neglect here a small inconsistency in exchange-correlation potential and energy density, due to the presence of a small amount of components beyond the $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$ set). This has important consequences: working in G-space means that translational invariance is exactly conserved and that forces are analytical derivatives of the energy (apart from the effect of the small inconsistency mentioned above). Forces that are analytical derivatives of the energy ensure that the constant of motion (i.e., the sum of kinetic and potential energy of the ions in Newtonian dynamics) is conserved during the evolution.

#### II.4.3 Double-grid technique

Let us focus on US PPs. In G-space the charge density is:

$`n(𝐆)=n_{\mathrm{soft}}(𝐆)+{\displaystyle \underset{i,nm,I}{}}Q_{mn}^I(𝐆)\varphi _i|\beta _n^I\beta _m^I|\varphi _i.`$ (49)
If $`E_\mathrm{c}^{\mathrm{wf}}`$ is the cutoff for the KS orbitals, the cutoff for the soft part of the charge density is $`E_\mathrm{c}^{\mathrm{soft}}=4E_\mathrm{c}^{\mathrm{wf}}`$. The augmentation term often requires a cutoff higher than $`E_\mathrm{c}^{\mathrm{soft}}`$, and as a consequence a larger set of G-vectors. Let us call $`\{𝐆_\mathrm{c}^{\mathrm{dens}}\}`$ the set of G-vectors that are needed for the augmented part:

$`{\displaystyle \frac{\mathrm{}^2}{2m}}G^2<E_\mathrm{c}^{\mathrm{dens}}.`$ (50)
In typical situations, using pseudized augmented charges, $`E_\mathrm{c}^{\mathrm{dens}}`$ ranges from $`E_\mathrm{c}^{\mathrm{soft}}`$ to $`23E_\mathrm{c}^{\mathrm{soft}}`$.

The same FFT grid could be used both for the augmented charge density and for KS orbitals. This however would imply using an oversized FFT grid in the most expensive part of the calculation, dramatically increasing computer time. A better solution is to introduce two FFT grids:

* a coarser grid (in r-space) for the KS orbitals and the soft part of the charge density. The FFT dimensions $`N_1,N_2,N_3`$ of this grid are big enough to accommodate all G-vectors in $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$
* a denser grid (in r-space) for the total charge density and the exchange-correlation and Hartree potentials. The FFT dimensions $`M_1N_1,M_2N_2,M_3N_3`$ of this grid are big enough to accommodate all G-vectors in $`\{𝐆_\mathrm{c}^{\mathrm{dens}}\}`$.

In this framework, the soft part of the electron density $`n_{\mathrm{soft}}`$, is calculated in r-space using FFTs on the coarse grid and transformed in G-space using a coarse-grid FFT on the $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$ grid. The augmented charge density is calculated in G-space on the $`\{𝐆_\mathrm{c}^{\mathrm{dens}}\}`$ grid, using Eq. (49) as described in the next section. $`n(𝐆)`$ is used to evaluate the Hartree potential, Eq.(48). Then $`n(𝐆)`$ is Fourier-transformed in r-space on the dense grid, where the exchange-correlation potential, Eq.(47), is evaluated.

In real space, the two grids are not necessarily commensurate. Whenever the need arises to go from the coarse to the dense grid, or vice versa, this is done in G-space. For instance, the potential $`V_{\mathrm{eff}}`$, Eq. (13), is needed both on the dense grid to calculate quantities such as the $`D_{nm}^I`$, Eq. (16), and on the coarse grid to calculate $`V_{\mathrm{eff}}\varphi _j`$, Eq. (11). The connection between the two grids occurs in G-space, where Fourier filtering is performed: $`V_{\mathrm{eff}}`$ is first transformed in G-space on the dense grid, then transferred to the coarse G-space grid by eliminating components incompatible with $`E_\mathrm{c}^{\mathrm{soft}}`$, and then back-transformed in r-space using a coarse-grid FFT.

We remark that for each time-step only a few dense-grid FFT are performed, while the number of necessary coarse-grid FFTs is much larger, proportional to the number of KS states $`N_{\mathrm{ks}}`$.

#### II.4.4 Augmentation boxes

Let us consider the augmentation functions $`Q_{nm}`$, which appear in the calculation of the electron density, Eq. (49), in the calculation of $`D_{nm}^I`$, Eq. (16), and in the integrals involving $`Q_{nm}^I/𝐑_I`$ needed to compute the ionic forces, Eq. (23). The calculation of the $`Q_{nm}`$ in G-space has a large computational cost because the cutoff for the $`Q_{nm}`$ is the large cutoff $`E_\mathrm{c}^{\mathrm{dens}}`$. The computational cost can be significantly reduced if we take advantage of the localization of the $`Q_{nm}`$ in the core region.

We call “augmentation box” a fraction of the supercell, containing a small portion of the dense grid in real space. An augmentation box is defined only for atoms described by US PPs. The augmentation box for atom $`I`$ is centered at the point of the dense grid that is closer to the position $`𝐑_I`$. During a MD run, the center of the $`I`$th augmentation box makes discontinuous jumps to one of the neighboring grid points whenever the position vector $`𝐑_I`$ gets closer to such grid point. In a MD run, the augmentation box must always contain completely the augmented charge belonging to the $`I`$th atom; otherwise, the augmentation box must be as small as possible. Augmentation boxes of different sizes for different atoms could in principle be used, but in our implementation the same box size is chosen for all the atoms. Thus the atomic species having the less localized augmented charge determines the size of all the augmentation boxes.

The volume of the augmentation box is much smaller than the volume of the supercell. The number of G-vectors in the reciprocal space of the augmentation box is smaller than the number of G-vectors in the dense grid by the ratio of the volumes of the augmentation box and of the supercell. As a consequence, the cost of calculations on the augmentation boxes increases linearly with the number of atoms described by US PPs.

Augmentation boxes are used twice in the calculation:

* to construct the augmented charge density, Eq. (6),
* to calculate the self-consistent contribution to the coefficients of the nonlocal PP, Eq. (16).

In case (i), the augmented charge is conveniently calculated in G-space, following Ref.CarPasquarello , and Fourier-transformed in r-space. All these calculations are done on the augmentation box grid. Then the calculated contribution at each r-point of the augmentation box grid is added to the charge density at the same point in the dense grid. In case (ii), it is convenient to calculate $`D_{nm}^I`$ as follows: for every US atom, take the Fourier transform of $`V_{\mathrm{eff}}(𝐫)`$ on the corresponding augmentation box grid and evaluate the integral of Eq. (16) in G-space.

## III Parallel Ultrasoft Pseudopotential Implementation

Various parallelization strategies for PW-PP calculations have been described in the literature. A strategy that ensures excellent scalability in terms of both computer time and memory consists in distributing the PW basis set and the FFT grid points in real and reciprocal space across processors. A crucial issue for the success of this approach is the FFT algorithm, which must be capable of performing three-dimensional FFT on data shared across different processors with good load balancing ParaFFT . This algorithm can be generalized to the US case as described in the following subsection.

### III.1 Parallel FFT in the US case

Partitioning a real-space FFT grid across processors is straightforward. The FFT grid, Eq.(43), is subdivided in a number of slices equal to the number of processors, so that each processor can take care of a different slice. The slices are cut along planes orthogonal to a chosen crystallographic direction. We label the crystallographic directions by 1,2,3. For instance, let us consider a FFT grid with $`N_3`$ planes along direction 3, which is distributed across $`N_p`$ processors. If $`N_p`$ is a divisor of $`N_3`$, good load balancing is achieved if each slice contains the same number $`(N_3/N_p)`$ of planes. Processor $`p`$ contains planes with $`m_3`$ values such that: $`(p1)(N_3/N_p)m_3p(N_3/N_p)1`$. If $`N_p`$ is not a divisor of $`N_3`$, all the slices cannot be equal. In this case their dimension is chosen in such a way as to minimize load imbalance. If $`N_p`$ exceeds the number of planes $`N_3`$, this strategy has to be refined.

The partition of the G-space grid is more involved. The Fourier components of the quantities of interest (e.g. the orbitals, the charge density, etc.) are stored as vectors (one-dimensional arrays): $`f(i)f(𝐆_i)`$, where the index $`i`$ spans one of the three sets of G-vectors defined above, namely the set $`\{𝐆_\mathrm{c}^{\mathrm{wf}}\}`$, the set $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$, and the set $`\{𝐆_\mathrm{c}^{\mathrm{dens}}\}`$. When a FFT is needed, the Fourier components have to be transferred to one of the two grids (three-dimensional arrays), defined by Eq.(44). The two grids are either the coarse grid, with dimensions $`N_1,N_2,N_3`$, or the dense grid, with dimensions $`M_1,M_2,M_3`$. The Fourier components must be evenly distributed across processors in order to achieve optimal load balancing for operations like scalar products. At the same time, their distribution across processors should achieve good load balancing in the FFTs and minimize the amount of data communication needed to perform the FFTs.

For each pair $`n_1^{},n_2^{}`$ in Eq.(44) we define a “column” in G-space, including all $`𝐆_{n_1^{},n_2^{},n_3^{}}`$ with $`M_3/2n_3^{}M_3/2`$. Since the KS orbitals have nonzero Fourier component only for G-vectors belonging to the set $`\{𝐆_\mathrm{c}^{\mathrm{wf}}\}`$, only a subset of all the columns contributes to a one-dimensional FFT of a KS orbital in the direction $`3`$. We call these columns “active columns” for the set $`\{𝐆_\mathrm{c}^{\mathrm{wf}}\}`$. In general, the number of nonzero Fourier components is different for each active column. Ideally, we would like to distribute the active columns across the processors, so that each processor receives the same number of active columns and the same number of Fourier components. Although not possible in general, this can be achieved to a good extent with a simple algorithmSdG : 1) create a list of columns, ordered by decreasing number of nonzero Fourier components; 2) assign the column to the processors, following the order in the list; 3) when all the processors contain at least one column, assign the following column in the list to the processor with the smallest number of nonzero Fourier components. This algorithm works nicely when the number of columns per processor is large enough.

After assigning to the processors all the columns that are active for the set $`\{𝐆_\mathrm{c}^{\mathrm{wf}}\}`$, we distribute across the processors the remaining columns, that are active for the set $`\{𝐆_\mathrm{c}^{\mathrm{soft}}\}`$, using the same algorithm. Finally we distribute across the processors the remaining columns, that are active for the set $`\{𝐆_\mathrm{c}^{\mathrm{dens}}\}`$, again using the same algorithm. The remaining columns are not active for any set of G-vectors and play no role.

After distributing all the columns across the processors, a one-dimensional FFT along direction 3 is done on local data (on a single processor). However, the data on the planes orthogonal to direction 3 are distributed across the processors. In order to perform FFTs in each of these planes, the corresponding data must be made local to a processor. This is achieved by a parallel transpose operation, performed with a single call to the appropriate MPIMPI  library routine. Two-dimensional FFTs can then be performed on the planes, with each processor operating on local data. Nonzero contributions are present only for $`(n_1^{},n_2^{})`$ pairs corresponding to active columns. This fact can be exploited to reduce the number of FFT operations, by performing only the FFTs along direction 1 (or 2) that include nonzero contributions. The strategy for parallel three-dimensional FFT that we have presented requires the number of processors $`N_p`$ be smaller than or equal to the number of planes $`N_3`$. The FFT from r\- to G-space uses the same algorithm in reversed sequence.

In calculations using only the $`\mathrm{\Gamma }`$ point ($`𝐤=0`$), the KS orbitals can be chosen to be real functions in r-space, so that $`\varphi (𝐆)=\varphi ^{}(𝐆)`$. This allows us to store only half of the Fourier components. Moreover, two real FFTs can be performed as a single complex FFT. To this end the auxiliary function $`\mathrm{\Phi }`$ is introduced:

$`\mathrm{\Phi }(𝐫)=\varphi _j(𝐫)+i\varphi _{j+1}(𝐫)`$ (51)
whose Fourier transform $`\mathrm{\Phi }(𝐆)`$ yields

$`\varphi _j(𝐆)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }(𝐆)+\mathrm{\Phi }^{}(𝐆)}{2}}`$ (52)
$`\varphi _{j+1}(𝐆)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }(𝐆)\mathrm{\Phi }^{}(𝐆)}{2i}}.`$ (53)
A side effect on parallelization is that $`𝐆`$ and $`𝐆`$ must reside on the same processor. As a consequence, pairs of columns with $`𝐆_{n_1^{},n_2^{},n_3^{}}`$ and $`𝐆_{n_1^{},n_2^{},n_3^{}}`$ (with the exception of the case $`n_1^{}=n_2^{}=0`$), must be assigned to the same processor.

### III.2 Scalar products

All scalar products $`f|g=_if_i^{}g_i,i=1,n`$ where $`i`$ runs on a distributed grid, can be calculated by calling standard optimized library routines (like BLAS from NetLibNetLib ) on each processor, and subsequently by summing the partial results of all processors, using a call to standard MPIMPI  libraries. Scalar products between vectors for which only half of the Fourier components are stored require a special treatment. Let $`n_p`$ be the number of Fourier components stored on processor $`p`$. The contribution of this processor to the scalar product is $`f|g_p=2_{i=1,n_p}f_i^{}g_i`$ if the $`𝐆=0`$ components are not within the set of $`n_p`$ components. If instead the $`𝐆=0`$ components, identified by $`i=1`$, are stored on processor $`p`$, the contribution of processor $`p`$ to the scalar product is $`f|g_p=f_1g_1+2_{i=2,n_p}f_i^{}g_i`$.

### III.3 Iterative Orthonormalization

The scalar products in the matrix elements, Eq. (29), needed to compute the Lagrange multipliers are calculated in parallel, following the procedure of the previous subsection. In the present implementation, the solution of the matrix equation (35), involving square matrices of dimension equal to the number $`N_{\mathrm{ks}}`$ of KS orbitals, is not parallelized but replicated on all the processors. Usually the time spent in the non parallelized part of the iterative orthonormalization is only a small fraction of the total time of the calculation.

To efficiently perform calculations on very large systems, using a large number of processors, the solution of Eq. (35) should also be parallelized. The time consuming steps are matrix-matrix multiplication and the diagonalization of the $`B`$ matrix. Both calculations require $`𝒪(N_{\mathrm{ks}}^3)`$ floating-point operations. A convenient parallelization approach is described in Ref.Cavazzoni .

### III.4 Augmentation boxes

The parallelization of the calculations performed on the augmentation boxes is not obvious for two reasons: 1) each augmentation box has a grid which is a portion of the dense grid and is distributed across processors; 2) the boxes follow the atoms in the MD evolution, causing the portion of the dense grid to change with time. In the present implementation, we deal with these difficulties as follows. We keep on all processors a copy of all the quantities defined on the augmentation boxes. Calculations on the grid of a given augmentation-box are performed only in the processors that contain at least a fraction of the given augmentation-box grid. This causes some replication of the calculations. FFTs on the augmentation box grid are performed locally on each processor. In order to reduce the amount of replication, in the FFTs from G\- to r-space, the two-dimensional FFTs along planes orthogonal to direction 3 are performed only in the planes belonging to the slice of the dense grid that is local to a given processor. No communication is needed to copy the augmented charge in r-space from the augmentation-box grid to the dense grid (see Sec. II.4.4). In the calculation of $`D_{nm}^I`$, we evaluate $`Q_{nm}^I`$ in G-space, transform it in r-space using augmented-box FFT, evaluate the integral of Eq. (16) in r-space,and sum the final result over all processors. This approach keeps communications to a minimum, at the expense of a number of augmentation-box FFTs larger than in the serial case.

Augmentation-box grid related calculations constitute a very small part of the overall computational cost, both in computer time and in memory. Therefore the simple approach that we have presented is convenient even if some calculations are replicated on few processors and the load balance is not optimal.

## IV Test case: iron porphyrins

We report here a comparison of computer performances for US and standard PPs in CP calculations. Our test systems–prototypes of systems containing the iron-porphyrin motif–are a reduced and an extended model of the active site of myoglobin.

### IV.1 Models and computational details

The reduced model is composed of an iron-porphyrin-imidazole complex, already investigated using the CP method by Rovira et al. Rovira ; the metallic pentacoordinated center is bound to the four planar porphyrin nitrogens, with the imidazole nitrogen occupying one of the axial sites, binding approximately orthogonal to the porphyrin plane, see Figure 1. The chemical formula is \[FeN<sub>6</sub>C<sub>23</sub>H<sub>16</sub>\]. A simple cubic cell of size 15.875 Å, containing a total of 46 atoms and 154 electrons, is used. For the reduced model, we compare both spin-restricted and unrestricted (S=2) calculations.

The extended model is composed by a large portion of the myoglobin active site, defined by the full heme group (same coordination as for the reduced model) plus the 13 surrounding residues which were comprised within a sphere of 8 Å radius centered on the iron atom, see Figure 2. The initial geometry has been taken from the X-ray experimental structure of the O<sub>2</sub>-myoglobin complex X-ray  and the included residues have been terminated by NH<sub>2</sub> groups, resulting in a total of 332 atoms and 902 electrons. The chemical formula is \[FeO<sub>19</sub>N<sub>35</sub>C<sub>106</sub>H<sub>173</sub>\]. A simple cubic cell of size 25.4 Å has been used, ensuring a minimum separation of 5 Å between periodic replicas. For the extended model we discuss only the performances of the more computationally demanding spin-unrestricted (S=2) calculations.

For a correct comparison of performances we need to compare data of similar quality, in terms of accuracy of chemical properties, obtained with algorithms of comparable quality in terms of serial speed and parallel speedup. In order to satisfy the first requirement, we need to determine a set of cutoffs for the US and standard calculations yielding comparable structural properties. To this end we compared the optimized geometry of the triplet ground state for the reduced model from our US calculations and from published NC results, performed at $`E_\mathrm{c}^{\mathrm{wf}}=70`$ Ry Rovira . Our geometrical parameters are converged to the same extent as those in Ref. Rovira  at $`E_\mathrm{c}^{\mathrm{wf}}=25`$ Ry and $`E_\mathrm{c}^{\mathrm{dens}}=200`$ Ry; in particular the critical Fe-N distances in the porphyrin and in the imidazole are found to be 2.00 and 2.12 Å respectively, vs. 2.00 and 2.14 Å of Ref. Rovira , with the out-of-plane displacement of the iron atom computed to be 0.14 vs. 0.15 Å. The overall agreement between the two sets of results is excellent, and the residual difference can be attributed to the different functionals used: BP86 BP86  in Ref.Rovira , and PW91 PW91  in our calculations. We estimate that we can safely compare US-PP calculations performed at $`E_\mathrm{c}^{\mathrm{wf}}=25`$ Ry to standard calculations at $`E_\mathrm{c}^{\mathrm{wf}}=70`$ Ry. For a fair comparison we use the same cutoff for the charge density in the US case as in the standard case ($`E_\mathrm{c}^{\mathrm{dens}}=280`$ Ry). We also performed calculations at $`E_\mathrm{c}^{\mathrm{wf}}=35`$ Ry for the US case, at $`E_\mathrm{c}^{\mathrm{wf}}=100`$ Ry for the standard case ($`E_\mathrm{c}^{\mathrm{dens}}=400`$ Ry in both cases).

In order to compare algorithms of similar quality, all calculations were performed using the same code CPV  (standard PPs are just a special case of US PPs). The PPs used in the standard case tests were generated using the technique of Troullier and Martinstm . In the US case, we use US PPs for all atoms, including H. The PW91 functionalPW91  is used in all calculations.

### IV.2 Results

The results for the reduced models were obtained on a 32-node IBM SP3 ($`4\times 375`$MHz power3 processors per node), while the extended model calculations were performed on a 64-processor SGI Origin ($`64\times 300`$MHz RISC 12000 processors), both at the Keck Materials Science Laboratory, Princeton University.

The reported execution times are an average over 20 time steps of the measured wall time (the sum of CPU and system time, differing by only a few percent from pure CPU time).

The parallelization performances of the US vs. NC-PPs implementation for the spin-restricted case of the reduced model are shown in Table 1. The execution times for a CP step (including calculation of forces and time evolution of atomic positions) are about 15 % larger than for a purely electronic step (orbital time evolution only), as expected. The small superlinear speedup observed both for US and NC PPs is a consequence of caching: since the memory per processor decreases almost linearly with the number of processors, better caching can be achieved with an increased number of processors, thus increasing the serial speed of the code. It is worth noting that US calculations are faster by a factor $`2.5`$ with respect to the NC case and require half RAM memory and $`1/4`$ disk space with respect to standard calculations.

The performances of US vs. NC-PPs calculations at higher cutoff are shown in Table 2. The number of PWs is approximately $`(35/25)^{3/2}1.65`$ times larger than in the preceding case. Execution times should theoretically be proportional to the same factor. The factor is actually somewhat larger ($`1.9`$), but the cache effects mentioned above and the effect of the discreteness of the FFT grid explain the difference. Again, US calculations are faster by a factor $`2.5`$ with respect to the NC case and require half RAM memory with respect to standard calculations.

In Table 3 we report spin-unrestricted results, showing an approximate doubling of execution time and of memory requirements, in line with expectations.

Table 4 contains the performances of US calculations on the quintet state ($`S=2`$), which is experimentally known to be the ground spin state in myoglobin, of our extended model of myoglobin at $`E_\mathrm{c}^{\mathrm{wf}}=25`$ Ry and $`E_\mathrm{c}^{\mathrm{dens}}=200`$ Ry. The scaling with the number of processors is excellent in this case too, with a slight superlinear scaling up to 32 processors. The ratio between execution times for electronic and CP time-steps is almost the same in the extended and in the reduced models. A geometry optimization time step requires less than 20 minutes in the spin-unrestricted case on 48 processors, with a total RAM usage of $`21`$ Gb. The total size of the file containing the KS orbitals is 6.7 Gb. In this case no standard calculation was attempted; indeed, an extrapolation from the results of Tables 1, 3 and 4 points to a total memory requirement of more than 40 Gb.

A typical local geometry optimization requires $`250`$ time-steps. An execution time of 20 minutes per time step thus translates into less than 4 days for the optimization to complete. This is a perfectly feasible calculation, not even requiring a state-of-the art massive parallel computer. On the other hand, a typical MD run requires no less than 10000 time-steps, corresponding in the present case to a few ps of simulation time. A true dynamical simulation would therefore become accessible on a state-of-the art massive parallel computer.

The relevance of simulating extended portions of myoglobin active site can be understood by considering the spin density distribution of the quintet spin state of the extended model, reported in Figure 3 together with selected integrated spin density values. The spin density is mainly localized on the iron atom (ca. 88 $`\%`$), even though a sizable contribution (ca. 12 $`\%`$) is computed to be delocalized over the rest of the system, with largest contributions arising from the propionate groups bound to the porphyrin ring. This finding is of particular interest, considering that CO rebinding in myoglobin has been recently related to a spin crossover from the quintet spin state, corresponding to unbound CO + myoglobin, to the singlet spin state characterizing the bound configuration of the CO-myoglobin complex. Harvey  Since explicit inclusion of the protein environment alters the spin distribution of the quintet state in our extended model, an effect on the relative energy of the different spin states can be expected.

## V Accuracy of Periodic Boundary Conditions

The use of PBC to describe molecules is perfectly appropriate for neutral molecules with small dipole/quadrupole moments, provided that the chosen supercell is large enough to minimize the spurious interactions between periodic replicas. This goal can usually be reached with supercells that leave a few Å of empty space between periodic replicas.

Charged molecules should be described by charged supercells, but these have infinite electrostatic energy. A finite energy is obtained by setting to zero the divergent $`G=0`$ contributions to the energy, as if the system were neutral. This is equivalent to adding a neutralizing background. Energies obtained in this way will be referred to as “uncorrected”. The direct comparison of uncorrected energies between different charge states is usually meaningless, because the error induced by PBC is large in this case. The long-range character of Coulomb interactions would require unpractically large supercells. Uncorrected energies may be affected by a large error also in molecules with large dipole/quadrupole moments.

Several techniques have been devised to overcome such limitation. The Hockney techniqueHockney  yields an exact treatment of charged species using PWs without imposing PBC. This is achieved by cutting the Coulomb potential in real space beyond a suitably chosen cutoff that excludes all spurious interactions between periodic replica, still taking into account intramolecular interactions. This technique is rather expensive, since it requires the definition of an enlarged FFT grid for the Coulomb potential. A similar techniqueMartyna  where the cutoff acts in reciprocal space allows for faster execution with minor loss of accuracy.

A much simpler and approximate technique, due to Makov and Payne (MP) MakovPayne , consists in calculating the leading electrostatic correction terms and removing them from the uncorrected total energy. The MP correction is performed a posteriori on the energy only: the effect of the charge on the potential and on atomic forces is therefore neglected.

The Makov-Payne corrected energy $`E_{MP}`$ for a cubic supercell has the form:

$$E_{MP}=E+\frac{q^2\alpha }{2L}\frac{2\pi qQ}{3L^3}$$
(54)
where $`E`$ is the uncorrected energy, $`q`$ is the net charge, $`Q`$ is the quadrupole moment, $`\alpha `$ is the Madelung constant, $`L`$ is the supercell side. This is Eq.(15) of Ref.MakovPayne  with the correct sign. When applying Eq. (54), the origin has to be translated so that the dipole moment vanishes MakovPayne . Therefore, the calculation of the Makov-Payne correction is straightforward since it requires only calculation of the dipole and quadrupole moments.

### V.1 Models and computational details

We want to verify the ability of PW-PP calculations with the MP correction to reproduce the electronic and structural properties of highly charged species. We compute the energy difference between two meso-substituted Mn(V) porphyrins: the oxo-aquo-Mn(V)TM-2-Pyridyl (2-Pyp) and -Mn(V)TM-4-Pyridyl (4-Pyp) porphyrins, I (see Figure 4), and between the corresponding oxo-hydroxo species, II, in which the axial water molecule has been replaced by a OH<sup>-</sup> ligand. We compare our results to calculations employing a localized basis set of Slater type orbitals (STO).

The oxo-aquo and oxo-hydroxo porphyrins, recently experimentally characterized as mimics of the halide oxidation reaction catalyzed by haloperoxidasesGroves1 ; Groves2 , have a charge $`+5`$ and $`+4`$, respectively, with 4 positive charges approximately localized on the aryl nitrogens and, in the case of the oxo-aquo species, the residual positive charge located at the metal center; the two isomeric porphyrins differ, both in the oxo-aquo and oxo-hydroxo form, for the position of the methyl-substituted nitrogen in the aromatic ring attached to the meso porphyrin carbons, which should lead to a considerable increase of the quadrupole moment from the 2-Pyp to the 4-Pyp isomer. Due to charge differences between oxo-aquo and oxo-hydroxo species, to the high total charge and to the large expected difference in the quadrupole moment between the 2-Pyp and 4-Pyp porphyrins, we believe that the calculation of the relative energies of the two isomeric porphyrins represents a severe test for PW calculations within PBC.

We consider a reduced model of the real systems, in which the methyl groups bound the the aryl nitrogens are replaced by hydrogens. The US-PP results are compared to those obtained by using STOs in the frozen core approximation. US-PP calculations were performed by using a cutoff of 25-200 Ry for the KS orbitals and density, respectively, and a cubic cell of side 19.05 Å, without any symmetry constraints.

STOs results were obtained using the ADF program ADF1 ; ADF2 ; the frozen cores include 1$`s`$-2$`p`$ states for Mn, 1$`s`$ states for O, N and C. The KS orbitals were expanded in an uncontracted DZ STO, standard basis set II basis  for all atoms with the exception of the transition metal for which we used a standard basis set IV basis , of TZP quality. STOs calculations were performed within C<sub>2v</sub> and C<sub>S</sub> symmetry constraints for species I and II, respectively.

### V.2 Results

In Table 5 we compare the relative energy of the singlet ground states Groves1  of the two isomeric 2-Pyp and 4-Pyp porphyrins, for both the oxo-aquo and oxo-hydroxo species. For the PP calculations we report both uncorrected and Makov-Payne corrected energy differences. Table 6 contains the calculated values of the quadrupole moment used in Eq. (54).

Results obtained with STOs localized basis sets compute the oxo-aquo 4-Pyp system to be ca. 26 kcal/mol more stable than the 2-Pyp one. The Makov-Payne corrected US-PP results are in excellent agreement with STOs results, while uncorrected PP results indicate instead that the 4-Pyp isomer is more stable than the 2-Pyp one by only 9.5 kcal/mol. Moreover, in the case of the oxo-hydroxo species, the uncorrected results yield an incorrect energy ordering, with the 2-Pyp isomer computed to be more stable than the 4-Pyp one by 7.0 kcal/mol. The discrepancy is resolved upon correcting the total energies with the Makov-Payne term, resulting again in an excellent agreement with the STOs energy differences.

Interestingly, the geometrical structures of the investigated charged systems calculated using the US-PP approach with PBC, turn out to quantitatively compare with results obtained using localized basis sets, see Table 7 for a comparison of main optimized geometrical parameters of species I. The main discrepancy (0.03 Å) is computed for the formally triple Mn$``$O bond, probably because of the lack of polarization functions in the O STO basis set which in turn leads to an overestimate of such parameter. The agreement between US PPs and STOs results suggests that the error introduced on the electrostatic potential by the presence of a charge in PBC does not significantly affect the structural properties. On the other hand, the effect on the total energy is sizable but mostly corrected by the use of Eq. (54).

### V.3 Conclusions

We believe that our results demonstrate that the Car-Parrinello approach in conjunction with ultrasoft pseudopotentials represents a valuable and relatively cheap tool to describe the electronic and geometrical properties of complex bio-inorganic systems, including highly charged and open-shell species. The study of the electronic and geometrical properties of such systems can now be achieved at a reasonable computational cost on conventional parallel machines with a limited number of processors. Moreover, a simple correction allows us to calculate energy differences in charged systems with an accuracy that is comparable to that of localized basis-set calculations not using Periodic Boundary Conditions. First-principle molecular dynamics simulations of extended bio-inorganic systems may still be feasible only on state-of-the art massive parallel computer. The key to reach such a goal is the availability of parallel machines with increased performances and number of processors and highly optimized scalable algorithms. We believe that the present parallel implementation of the Car-Parrinello method using ultrasoft pseudopotentials provides such an algorithm that will allow in the near future the simulation of the dynamical properties of such complex systems.

Acknowledgments Calculations were performed at Keck computational facility of the Princeton Materials Institute and at ISTM. Support from the NSF (CHE-0121432) is acknowledged. We wish to thank Prof. J. T. Groves, Prof. T. G. Spiro and Dr. A. Jarzecki for helpful discussions. FDA thanks CNR (Progetto Finalizzato “Materiali Speciali per Tecnologie Avanzate II”) for financial support. PG thanks MIUR grant PRIN 2001-028432 for partial support.


text:
# Jet Production at HERA

## 1 Introduction

Jet production allows precise tests of perturbative QCD. At the electron-proton collider HERA the production of jets is investigated in inelastic electron (positron) proton scattering over a wide range of photon virtualities $`Q^2`$, from photoproduction ($`Q^20\text{GeV}^2`$) to deep inelastic scattering (DIS). This article presents recent measurements of jet cross sections and jet shapes and their comparison to perturbative QCD (pQCD) calculations. These calculations are performed up to next-to-leading order (NLO) for all measurements presented. Jet cross sections at HERA are calculated by convoluting the short distance subprocess cross sections with parton density functions (pdfs) of the proton. At very small $`Q^2`$ (photoproduction) the photon may fluctuate into partons before the hard scatter. This is parametrised by the pdfs of the photon.

The standard choice since about 5 years at HERA is to define jets according to the theoretically preferred invariant $`k_{}`$ algorithm as proposed by Ellis and Soper. ZEUS and H1 have reduced their energy scale uncertainty for high transverse energy ($`E_T`$) jets to $`1\%`$, respectively $`2\%`$, which significantly improves the data precision.

The data at high $`Q^2`$ and high $`E_T`$ allows quantitative tests of pQCD and extractions of $`\alpha _s`$ and the proton pdfs. Photoproduction data at high $`E_T`$ are, in addition, sensitive to the photon pdfs. At moderate $`Q^2`$ and $`E_T`$ the limits of NLO pQCD are exploited. These measurements provide additional and complementary information to $`p\overline{p}`$ and $`e^+e^{}`$ data.

## 2 Jets at moderate and high $`Q^2`$ DIS

Jet production in deep inelastic scattering events has been investigated over a wide range in $`Q^2`$. A recent H1 measurement  of inclusive jet production is performed for the medium $`Q^2`$ range between 5 and 100 $`\text{GeV}^2`$. The data are compared to LO and NLO pQCD calculations using $`E_T`$ as renormalization scale. The agreement is good, although NLO pQCD fails to describe the data in the forward region ($`\eta _{lab}>1.5`$), when both $`E_T`$ and $`Q^2`$ are small. This region corresponds to the largest NLO/LO corrections. Figure 1 shows the inclusive jet cross section $`d\sigma /dE_T`$ for $`1.5<\eta _{lab}<2.8`$ as a function of the transverse energy in the Breit frame $`E_T`$, in different regions of $`Q^2`$.

Especially at $`5<E_T<20`$ GeV and $`5<Q^2<10\text{GeV}^2`$ the renormalization scale uncertainties do not cover the uncertainties of the data. These theory uncertainties are estimated by a variation of the scale by a factor of 2 (0.5).

Both the large differences between NLO and LO calculations and the poor agreement between data and NLO shows the necessity of further theoretical progress.

As an example of the many jet measurements at high $`Q^2`$ a new ZEUS measurement of inclusive jet cross sections  for photon virtualities $`Q^2>125\text{GeV}^2`$ is discussed. A good agreement between data and NLO calculations is found at high $`Q^2`$. Figure 2 shows the differential cross section $`d\sigma /dQ^2`$ as a function of $`Q^2`$ for inclusive jet production requiring $`E_T>8`$ GeV in the Breit frame. Since there is good agreement between data and the NLO calculation the measured cross sections were used to determine $`\alpha _s(M_Z)`$. The QCD fit to H1 and ZEUS high $`Q^2`$ jet data yield $`\alpha _s`$ values in good agreement with the world average and of similar precision to that obtained by other experiments.

## 3 Internal jet structure

Observables which give insight into the internal jet structure are the mean subjet multiplicity and the integrated jet shape. At high $`E_T`$, where fragmentation effects are small, these observables are calculable in pQCD. The subjet multiplicity $`<n_{subjet}>`$, for instance, is derived by repeating the $`k_T`$ jet algorithm with a smaller resolution scale $`y_{cut}`$ considering particles lying within a jet and counting the number of subjets found inside it. ZEUS has compared new measurements  of these observables to QCD calculations at NLO. The jets are found in the laboratory frame. Figure 3 shows the mean subjet multiplicity as a function of $`y_{cut}`$ for $`E_T>15\text{GeV}`$ and $`Q^2>125\text{GeV}^2`$. The NLO pQCD prediction reproduces the data well. This shows, for the first time, the ability of NLO pQCD in describing the internal jet structure. A QCD fit of these data were used to derive $`\alpha _s`$ values, consistent with the world average.

## 4 Jets in photoproduction

Dijet photoproduction has been investigated by the H1 and ZEUS experiments . Both analyses present dijet cross sections as functions of various jet observables. Asymmetric cuts on the $`E_T`$ of the two jets are applied in order to avoid regions of phase space which are infrared sensitive in NLO calculations. Thus the H1 (ZEUS) analysis requires $`E_T`$ for the first and second jet of $`25`$ and $`15`$ GeV ($`14`$ and $`11`$ GeV). In both analyses the ratio of the measured cross sections to the theoretical prediction varies at small $`x_\gamma `$ by $`\pm `$10-15% when the cut on $`E_T`$ of the second jet is varied by $`\pm 5`$ GeV. This effect is covered by the NLO scale uncertainty at low $`x_\gamma `$ of $`1520`$%.

Figure 4 shows the dijet cross section, as measured by ZEUS, as a function of $`x_\gamma `$, which is an estimator of the momentum fraction carried by partons out of the photon.

Although the shape of the data for different $`E_T`$ is not perfectly modelled by the NLO prediction, strong conclusions can’t be drawn, because of the size of the theoretical uncertainties. These data are nevertheless sensitive enough to constrain the photon pdfs (quarks and gluon). Once more further theoretical understanding would be helpful in the interpretation of the data.

In Figure 5 the H1 dijet cross section as a function of the mean transverse energy of the two jets $`d\sigma /dE_{T,mean}`$ and of the transverse energy of the leading jet $`d\sigma /dE_{T,max}`$ is shown. The data are well described up to $`E_T70\text{GeV}`$ by the NLO calculation. The NLO scale uncertainty is not reduced significantly with increasing $`E_{T,max}`$, but decreases from $`\pm 20`$% for the first bins to less than $`\pm 5`$% for increasing $`E_{T,mean}`$. This shows the sensitivity to the choice of dijet $`E_T`$ cuts.

The momentum fraction carried by partons out of the proton is termed $`x_P`$. Figure 6 shows the cross section $`d\sigma /dx_p`$ as a function of $`x_p`$ for two different $`x_\gamma `$ regions.

Even at the highest $`x_p`$ the measured cross sections agree well with the QCD predictions. In this part of the phase space about 40% of the cross section is attributed to processes induced by gluons in the proton. The recent data could be used to further constrain the existing proton pdfs, particularly at medium $`x_P`$.

## 5 Summary

Recent HERA jet data show that NLO QCD calculations and current pdf parameterisations are mostly compatible with high $`E_T`$ and high $`Q^2`$ jet data, which allows precise determinations of $`\alpha _s(M_Z)`$ and pdfs. The increasing precision of the data shows that a deeper understanding of jet physics will require theoretical progress.