The SetAll dataset consists of the server models with a GDT-TS value published on the CASP7 website, a total of 23,423 models. To calculate a score on a protein model SELECTpro requires the backbone coordinates (N, C_α, C) for all model residues as input. A total of 8,812 models in SetAll have only a C_α trace or have no coordinates for one or more residues. Modeller8v1 [28-30] was used to generate complete models from the incomplete ones, and then the complete models were scored by SELECTpro. In most cases the complete models were built appropriately from the incomplete models; however, in some cases the final model was a fully extended chain due to an error in our application of Modeller. We failed to identify this problem until after the completion of the CASP7 competition. The SELECTpro scores versus GDT-TS scores for all models of target T0305 are displayed in plot A of Figure 4. The circled outliers with very low confidence scores and high GDT-TS scores are models that were incomplete and the complete models generated by Modeller were fully extended chains. The Pearson correlation on the set of all models for T0305 is .641. The SELECTpro scores versus GDT-TS scores for complete models only are displayed in plot B of Figure 4, and the Pearson correlation is .966.

The scores produced by SELECTpro are comparable on complete models of the same sequence. There is no standard for the handling of incomplete models and we assume that participating groups took a variety of approaches. Using only complete models ensures that the MQAP scores are calculated from the same coordinates. Thus, the models retained in SetComplete are screened first for completeness. Models missing backbone coordinates for one or more residues are removed. This leaves 14,611 models.

Structure-based MQAPs are susceptible to local structural irregularities in models, and will tend to score such models poorly. This is why methods developed to select near-native models from sets of decoys remove such models from consideration [31]. We apply additional filters (described below) for C_α-C_α clashes, C_α-C_α chain breaks, and expanded termini to remove an additional 1,217 models leaving 13,494 more plausible models in SetComplete.

The C_α-C_α clash model filter enforces a squared difference penalty for C_α-C_α distances less than 3.6 Ǻ. The distance between the C_α atoms of residue i and j is denoted by r_i,Cα,j,Cα and N is the protein length. The constant 13.52 in the threshold below corresponds to two severe clashes where r_i,Cα,j,Cα = 1.0 Ǻ. Models with a sum of squared differences greater than 13.52 per 100 residues are filtered out.

∑ i > j max ⁡ { 0 , 3.6 − r i , C α , j , C α } 2 > 13.52 ( N / 100 )

The C_α-C_α chain break model filter enforces a squared difference penalty for r_{i,Cα,i+1,Cα} distances greater than 4.0 Ǻ. The constant 16.0 in the threshold below corresponds to a single chain break where r_{i,Cα,i+1,Cα} = 8.0 Ǻ. Models with a sum of squared differences greater than 16.0 per 100 residues are filtered out.

∑ i max ⁡ { 0 , r i , C α , i + 1 , C α − 4.0 } 2 > 16.0 ( N / 100 )

The expanded termini filter removes models where a large portion of the structure is modeled as expanded chain with no non-local interactions. The screening procedure is: scan from the N-terminus until three consecutive residues have a contact number of at least 10, and repeat from the C-terminus. The contact number of a residue is defined here as the number of other C_β atoms within 10 Ǻ of the residue's C_β [3]. If the sum of low contact number termini residues is at least 20% of N, the model is filtered out.

In the reduced representation the heavy backbone atoms, carbonyl oxygen, amide hydrogen (N, C_α, C, O, H), and C_β are represented explicitly. For glycine residues a pseudo C_β is calculated. The side-chain atoms are represented by a single united point (centroid) [32,33]. The centroid is calculated as the mean of the position of the heavy side-chain atoms. For glycine and alanine the centroid (CT) is set to the C_β atom. Only the heavy backbone atoms (N, C_α, C) are used as input to SELECTpro and the positions of additional atoms and centroids are calculated from these.

In the all heavy-atom representation the centroid is removed and the heavy side chain atoms are represented explicitly. The side-chains are initially placed onto the backbone of the reduced representation in their most likely conformation according to the SCWRL backbone-dependent rotamer library [34]. The side-chain placements are then optimized using the SELECTpro all-atom energy terms (described below) in conjunction with the rotamer library.

The parameter weights were determined by repeatedly varying individual weights and maximizing the sum of the GDT-TS of the lowest E_FINAL models on a training set built from CASP6 protein domains. For each CASP6 protein domain a set of 500 decoy models was generated using fragment assembly with the RMSD to native as the dominant term in the objective function [3].

E_PRED-SS: predicted secondary structure

E_PRED-ACC: predicted solvent accessibility

E_PRED-CM: predicted contact map

E_BETA: sheet formation

E_BB-REP: backbone repulsion

E_CT-REP: centroid repulsion

E_STAT-ENV: residue environment potential [3]

E_STAT-PW-CI: context independent pair-wise potential [3,16]

E_STAT-PW-CD: context dependent pair-wise potential [6]

E_ROG: compactness

E_SC-HB: side-chain hydrogen bonding

E_LEN-JONES: van der Waals forces [10]

E_SOLVATION: solvation effects [35]

E_ELECTRO: electrostatic interactions

Throughout this work the convention of all capital letters referring to global energy and all lower case referring to local energy is used. For instance, E_PRED-CM refers to the global contact map energy and E_pred-cm(i,j) refers to the contact map energy between residues i and j.

r_i,x,j,y: distance between atom x of residue i and atom y of residue j

r_x,y: distance between atom x and atom y

v_i,x,j,y: vector from atom x of residue i to atom y of residue j

u_i,x,j,y: unit vector calculated from v_i,x,j,y

N_i: number of residues in contact with residue i, with contact defined as r_i,Cβ,j,Cβ < 10 Ǻ

phi_i: Phi angle of residue i

psi_i: Psi angle of residue i

aa_i: amino acid type of residue i

ss_i: predicted secondary structure of residue i (H,E,C)

acc_i: predicted solvent accessibility of residue i ('e': exposed, '-', buried)

cmap_i,j: predicted contact/non-contact between residues i and j, with contact defined as r_i,Cα,j,Cα < 12 Ǻ

I_value: ideal parameter value for a given calculation

σ_value: standard deviation value for a given calculation

vdw_x: van der Waals radius of atom x

vdw_x+y: vdw_x + vdw_y

Ω_stat-env: pre-calculated statistics for use in E_STAT-ENV

Ω_stat-pw-oi: pre-calculated statistics for use in E_STAT-PW-CI

Ω_stat-pw-od: pre-calculated statistics for use in E_STAT-PW-CD

D_min,pw-od: minimum interaction distance for centroid pairs used in E_STAT-PW-CD

D_max,pw-od: maximum interaction distance for centroid pairs used in E_STAT-PW-CD

D_min-CT: minimum distances between centroids of amino acid pairs observed in pdb_select25 [36].

The predicted structural feature predictions used in E_PRED-SS, E_PRED-ACC, and E_PRED-CM come from the SCRATCH suite of predictors [25]. Each predictor is trained in a supervised fashion using curated non-redundant datasets extracted from the PDB [37]. The secondary structure (SSpro [38]) and solvent accessibility (ACCpro [39]) predictors use ensembles of 1D-RNN (one dimensional-recursive neural network) architectures [40]. The contact map predictor (CMAPpro [41]) uses ensembles of 2D-RNN architectures [40].

The predicted secondary structure term E_PRED-SS penalizes deviation of the torsion angles from the torsion angle parameters for helices and strands predicted by SSpro. There is no penalty for predicted coils. The parameter values for helix residues are: I_Hφ = -65.3, σ_Hφ = 11.9, I_Hψ = -39.4, σ_Hψ = 11.3. The parameter values for strand residues are: I_Eφ = -135.0, σ_Eφ = 15.0, I_Eψ = 135.0, σ_Eψ = 15.0. Only torsion angles that are more than two standard deviations from the ideal are penalized, with the penalty defined as follows:

E P R E D − S S = ∑ s s i = H E p r e d − h e l i x ( i ) + ∑ s s j = E E p r e d − s t r a n d ( j ) E p r e d − h e l i x ( i ) = Θ H Φ ( i ) ( | p h i i − I H Φ | − 2 σ H Φ ) 2 + Θ H Ψ ( i ) ( | p s i i − I H Ψ | − 2 σ H Ψ ) 2 Θ H Φ ( i ) = { 1 ,     if | p h i i − I H Φ | > 2 σ H Φ 0,    otherwise Θ H Ψ ( i ) = { 1 ,     if | p s i i − I H Ψ | > 2 σ H Ψ 0,    otherwise

The definition of E_pred-strand(j) is equivalent to E_pred-helix(i), but with I_Eφ, σ_Eφ, I_Eψ and σ_Eψ in place of the corresponding helical values.

The solvent accessibility predictor ACCpro predicts the percent of solvent accessibility in 5% increments for each residue. Using 25% exposure as a binary threshold the accuracy of the predictor is ~77% [39]. The binary exposure ('e')/burial ('-') prediction is used as the predicted solvent accessibility for E_PRED-ACC. In the reduced representation the solvent accessibility of residue i is estimated by its contact number (N_i), where N_i > 16 is considered buried [3]. If the predicted status of a residue is not realized in the model, the penalty is calculated as:

E P R E D − A C C = ∑ i E p r e d − a c c ( i ) E p r e d − a c c ( i ) = { ( 17 − N i ) 2 ,  if  a c c i = ' − ' and  N i ≤ 16 ( N i − 16 ) 2 , if  a c c i = ' e ' and  N i > 16 0 , otherwise

The contact map predictor CMAPpro predicts the probability of contact or non-contact between C_α atoms, with a contact threshold of 12 Å. The strategy utilized to infer predicted contacts from the probability matrix [41] results in maps that are sparse when compared to those of real proteins; thus, unrealized contacts are penalized while non-contacts are not. The constant 1.0 is added to the penalty to ensure that all unrealized contacts make a significant contribution to E_PRED-CM.

E P R E D − C M = ∑ j > i Θ   ( i , j ) E p r e d − c m ( i , j ) Θ   ( i , j ) = { 1 , if  r i , C α , j , C α > I c m − t h r e s h 0 , otherwise E p r e d − c m ( i , j ) = c m a p i , j ( 1.0 + r i , C α , j , C α 2 I c m − t h r e s h 2 )

The predicted contact map can help identify the highest GDT-TS models in the set, even when they are not highly similar to native. A good example of this is CASP7 target T0304 is a 122 residue α/β protein where the highest GDT-TS model in the set is Zhang-Server_TS1 (GDT-TS = 45.55). Most secondary structure predictors (including SSpro) failed to predict the first two strands making this target especially difficult. No QA method ranked the highest GDT-TS model first; however, SELECTpro ranked it second and the model ranked first by SELECTpro (T0304.Zhang-Server_TS4) has the second highest GDT-TS. These models have the lowest E_PRED-CM of any models in the set, but the native structure has an even lower E_PRED-CM. Figure 5 compares the native and predicted contact maps for target T0304.

The formation of hydrogen bonds between the residues of β-strand partners is a major determinant of the tertiary structure of β and α/β proteins. The β hydrogen bonding treatment described here favors realistic strand pairing and sheet formation. The treatment also efficiently accommodates bulges in strands because it does not force the register between two paired strands. E_BETA is the global strand pairing energy that penalizes the hydrogen bonding of β residues between strand pairs. E_beta-sp(β_k→β_w) is the strand pairing energy of strand β_k to strand β_w. E_beta-sp is only commutative if the two strands have the same length. E_beta-hb(i,j) is the hydrogen bonding penalty between residues i and j.

E_beta-sp is calculated for all possible strand pairings, but only the two lowest energies from each strand are used in E_BETA. Other strand-strand interactions are ignored. In the equations below S is the set of all strands in the protein, β_m1 is the strand with the minimum pairing energy from β_k, and β_m2 is the strand with the next lowest pairing energy from β_k. If the strand count is less than six at least two of the strands must be edge strands. This is accounted for by only considering the single best strand partner for two strands.

E B E T A = ∑ β k ∈ S E b e t a − s p ( β k → β m 1 ) + E b e t a − s p ( β k → β m 2 ) β m 1 = { β x : min ⁡ β x ∈ S / { β k } E b e t a − s p ( β k → β x ) } β m 2 = { β y : min ⁡ β y ∈ S / { β k , β m 1 } E b e t a − s p ( β k → β y ) }

In the equations for E_beta-sp below, S_k is the set of all residues in strand β_k. Each time E_beta-hb is calculated the pair (i,j) is chosen with i from S_k and j from S_w, such that E_beta-hb is minimized. Then residue i is removed from S_k, and residue j is removed from S_w. E_beta-hb is calculated once for each residue in S_k. If S_k has more residues than S_w each unpaired residue is given maximum penalty of E_beta-hb.

E b e t a − s p ( β k → β w ) = ∑ S k ≠ ∅ , E b e t a − h b ( i , j ) ( i , j ) = { ( x , y ) : min ⁡ x ∈ S k , y ∈ S w E b e t a − h b ( x , y ) } S k = S k / { i } ,   S w = S w / { j }

Between two anti-parallel strand partners, only every other pair of residues is hydrogen bonded. For the pairs that are not hydrogen bonded, a pseudo-bonding calculation is used. The hydrogen bonding energy and pseudo-bonding energy are both calculated and the minimum of the two is used in E_beta-hb(i,j).

If residues i and j are paired in parallel strands, either i forms hydrogen bonds with j-1 and j+1, or j forms hydrogen bonds with i-1 and i+1. No hydrogen bonds are formed between the atoms of residues i and j. The hydrogen bonding energy is calculated for both possible conformations and only the minimum of the two is used in E_beta-hb(i,j).

E b e t a − h b ( i , j ) = { min{ Φ ( i → j ) + Φ ( j → i ) , Φ p s ( i → j ) + Φ p s ( j → i ) } , if strands are anti-parallel min{ Φ ( i → j + 1 ) + Φ ( j − 1 → i ) , Φ ( i − 1 → j ) + Φ ( j → i + 1 ) } , if strands are parallel

Φ(a→d) is the directional energy calculation for a single hydrogen bond where a is the index of the acceptor residue and d is the index of the donor residue. Three geometrical measures are used to estimate the strength of hydrogen bonds: the distance between the acceptor and the hydrogen atoms (r_a,O,d,H), the angle at the acceptor atom (u_a,C,a,O · u_a,O,d,H), and the angle between the acceptor and donor atom vectors (u_a,C,a,O · u_d,N,d,H). The distance and acceptor atom angle parameters are motivated by the orientation-dependent hydrogen bonding potential described in [42]. The following parameters were set based on idealized hydrogen bonding between β residues, with standard deviation values set such that two standard deviations approximate the cut-off in true hydrogen bonds. The ideal distance from hydrogen atom to accepting oxygen is I_hb-dist = 1.9 Ǻ, with standard deviation σ_hb-dist = 0.5 Ǻ. The ideal angle at the acceptor atom is 0°, so the ideal (u_a,C,a,O · u_a,O,d,H) is I_acc-dp = 1.0, with standard deviation σ_acc-dp = 0.11. The ideal angle between the acceptor and donor atom vectors is 180°, so the ideal (u_a,C,a,O · u_d,N,d,H) is I_acc-don-dp = -1.0, with standard deviation σ_acc-dp = 0.15. The parameters for pseudo-bonded residues are as follows: the ideal distance for r_a,O,d,H is I_ps-hb-dist = 7.9 Ǻ, I_ps-acc-dp = -1.0, and I_{ps-acc-don-dp} = -1.0. The standard deviations from the corresponding hydrogen bonding parameters above are used in Φ_ps(a→d).

Φ ( a → d ) = { 54.0 , if  r a , O , d , H > 7.0 Ψ ( r a , O , d , H , I h b − d i s t , σ h b − d i s t ) + Ψ (u a , C , a , O • u d , N , d , H , I a c c − d o n − d p , σ a c c − d o n − d p ) + Ψ (u a , C , a , O • u a , O , d , H , I a c c − d p , σ a c c − d p ) , otherwise Φ p s ( a → d ) = { 54.0 , if  r a , O , d , H > 10.0 Ψ ( r a , O , d , H , I p s − h b − d i s t , σ p s − h b − d i s t ) + Ψ (u a , C , a , O • u d , N , d , H , I a c c − d o n − d p , σ a c c − d o n − d p ) + Ψ (u a , C , a , O • u a , O , d , H , I p s − a c c − d p , σ a c c − d p ) , otherwise

The penalty for the observed value (x) increases up to 6 standard deviations from the ideal value (μ).

Ψ ( x , μ , σ ) = { ( x − μ ) 2 2 σ 2 , if  | x − μ | < 6 σ ( 6 σ ) 2 2 σ 2 = 18.0 , otherwise

E_SC-HB penalizes unsatisfied hydrogen bond donor and acceptor atoms that are at least partially buried. There is no penalty for fully exposed donor or acceptor atoms. Exposure percent (ΔGxslv %) is calculated as ΔGxslv/ΔGxref. The definitions of ΔGxslv and ΔGxref are provided in the description of E_SOLVATION in the Appendix. Atoms at least 75% exposed are considered fully exposed and atoms less than 25% exposed are considered fully buried. For 25% <ΔGxslv% < 75% the penalty weight is reduced linearly from 1.0 at 25% to 0 at 75%. The ideal distance from the acceptor atom to donor atom is I_hb-da-dist = 2.9 Ǻ. In the equations below donors is the set of all side-chain hydrogen donor atoms and acceptors is the set of all side-chain hydrogen acceptor atoms.

E S C − H B = ∑ x ∈ a c c e p t o r s E h b − a c c ( x ) + ∑ x ∈ d o n o r s E h b − d o n ( x ) E h b − a c c ( x ) = λ ( x ) min ⁡ y ∈ d o n o r s | r x , y − I h b − d a − d i s t | 2 E h b − d o n ( x ) = λ ( x ) min ⁡ y ∈ a c c e p t o r s | r x , y − I h b − d a − d i s t | 2 λ ( x ) = { 1 if  Δ G x s l v % < .25 0 if  Δ G x s l v % > .75 2 ( .75 − Δ G x s l v % ) otherwise

This term penalizes steric clashes between non-bonded atoms explicitly represented in the reduced representation. The penalty for overlapping atoms is the overlap distance squared as defined here:

E B B − R E P = ∑ j > i E b b − r e p ( i , j ) E b b − r e p ( i , j ) = ∑ x ∑ y Θ ( i , x , j , y ) ( v d w x + y − r i , x , j , y ) 2 Θ   ( i , x , j , y ) = { 1 , if   r i , x , j , y < v d w x + y 0 , otherwise

A centroid-centroid repulsive term is used to reduce the overcrowding of side-chains in the reduced representation. The minimum distance between two centroids in the calculation is the minimum observed for each pair of residue types – D_CT-min(aa_i,aa_j) – in pdb_select25. The penalty for centroid-centroid overlaps is defined as the overlap distance squared:

E C T − R E P = ∑ j > i Θ ( i , j ) [ D C T − min ⁡ ( a a i , a a j ) − r i , C T , j , C T ] 2 Θ ( i , j ) = { 1 , if  r i , C T , j , C T < D C T − min ⁡ ( a a i , a a j ) 0 , otherwise

The motivation for this term is to model the hydrophobic effect. The level of burial for each residue in the model is estimated by the number of other C_β atoms within 10 Ǻ (the contact number N_i) [3]. The values in the table Ω_stat-env reflect the likelihood of observing a particular N_i for each residue type. For model residues near both termini the contact number is artificially increased to account for the missing neighbors along the chain.

E S T A T − E N V = ∑ i Ω s t a t − e n v ( a a i , N i * ) N i * = { N i + 4 − i , if  i < 4 N i + 4 − | i − N | , if  | i − N | < 4 N i , otherwise

This context independent pair-wise potential comes from Equation 6 of [3]. The potential considers the likelihood of observing the pair of centroids in a given distance bin relative to the background, with distance bins of < 5, 5–7, 7–10, 10–12, and > 12 Å. The advantage of a context independent pair-wise potential is that it is less vulnerable to over-fitting by a conformational search because of its generality.

E S T A T − P W − C I = ∑ j > i E s t a t − p w − c i ( i , j ) E s t a t − p w − c i ( i , j ) = Ω s t a t − p w − c i [ a a i , a a j , r b i n ( i , j ) ] r b i n ( i , j ) = { 0 − 5 if  r i , C T , j , C T ≤ 5.0 5 − 7 if 5.0 < r i , C T , j , C T ≤ 7.0 7 − 10 if  7.0 < r i , C T , j , C T ≤ 10.0 10 − 12 if  10.0 < r i , C T , j , C T ≤ 12.0 12 + if  r i , C T , j C T > 12.0

This context specific pair-wise potential is from [6]. This pair-wise potential depends on the local structure and relative orientation of both amino acids in the interaction. The statistics are calculated independently for each combination of local structures and relative orientations. At each position the local structure is considered either compact or open and the relative orientation is determined by the dot product of the C_α to C_β unit vectors of each residue and divided into three classes: parallel, anti-parallel, and intermediate.

E S T A T − P W − C D = ∑ j > i E s t a t − p w − c d ( i , j ) E s t a t − p w − o d ( i , j ) = Θ ( i , j ) Ω s t a t − p w − o d [ a a i , a a j , λ ( i ) , λ ( j ) , Φ ( i , j ) ] Θ ( i , j ) = { 1 , if  r i , C T , j , C T > D min ⁡ , p w − o d [ a a i , a a j , λ ( i ) , λ ( j ) , Φ ( i , j ) ] and r i , C T , j , C T < D max ⁡ , p w − o d [ a a i , a a j , λ ( i ) , λ ( j ) , Φ ( i , j ) ] 0 , otherwise λ ( i ) = { compact , if  r i − 1 , C α , i + 1 , C α < 6.0 open, otherwise Φ ( i , j ) = { parallel , if u i , C α , i , C β • u j , C α , j , C β > .5 antiparallel , if u i , C α , i , C β • u j , C α , j , C β < − .5 intermediate , otherwise

The radius of gyration is a simple measure of the global compactness of a domain. E_ROG penalizes models that are less compact than expected according to [44]. If the radius of gyration of the model (λ) is less than the expected value (2.2N^.38), there is no penalty. If it is greater, then the penalty is the squared difference between observed and expected. In the equation below r_i,mean is the distance between the C_α of residue i and the mean of all C_αs in the model.

E R O G = Θ ( λ − 2.2 N .38 ) 2 λ = ∑ r i , m e a n 2 N Θ = { 1 , if  λ > 2.2 N .38 0 , otherwise

A fundamental characteristic of native globular protein structures is their efficient steric packing of atoms in the protein core. A Lennard-Jones 12-6 potential with damped repulsion (E_LEN-JONES) is used to measure the quality of steric packing. E_LEN-JONES is the sum of local energy calculations E_len-jones(x,y) performed on all pairs of non-bonded atoms. Since the repulsive portion of the standard Lennard-Jones 12-6 potential will overwhelm the entire energy function with a single significant atom-atom clash – repulsion is handled by a linear ramp from 0 to 10 as shown in the equation below [10]. Since E_len-jones = 0 when (vdw_x,y/r_x,y) = 26 independent of atom types, the switch to a linear ramp occurs when (vdw_x,y/r_x,y) > 26.

E L E N − J O N E S = ∑ y > x E l e n − j o n e s ( x , y ) E l e n − j o n e s ( x , y ) = { 10.0 ( 1 − 2 6 v d w x , y / r x , y ) , if  ( v d w x , y / r x , y )   > 2 6 ε x ε y [ ( v d w x , y r x , y ) 12 − 2 ( v d w x , y r x , y ) 6 ] , otherwise

Solvation energy is calculated using the implicit solvation model described in [35] with the following adjustment: for overlapping atoms, the sum of their van der Waals radii is used in the calculation in place of the observed atom-atom distance in the model. This restricts the amount a single atom can contribute to the burial of another atom. Without this adjustment overlapping atoms will bias the calculation to indicate an atom is more buried than it would be otherwise. In the solvation model ΔGxslv is the observed solvation free energy of atom x in the model, calculated as the free energy of the fully exposed atom (ΔGxref) minus the reduction in solvation caused by the surrounding atoms. ΔGxfree was determined empirically by setting it equal to ΔGxref and increasing its magnitude until ΔGxslv of deeply buried atoms became zero. λ_x is the correlation length of atom x. V_y is the volume neighboring atom y. The values of these parameters come from [3535], with the exception of ΔGxref[45]. The equation for ΔGxslv below is the combination of Equations 5,6, and 7 of [35], with the atom overlap adjustment.

E S O L V A T I O N = ∑ x Δ G x s l v Δ G x s l v = Δ G x r e f − Δ G x f r e e 2 λ x π π ∑ x ≠ y e − [ ( r x , y * − v d w x λ x ) 2 ] V y r x , y * 2 r x , y * = { v d w x + y , if   r x , y < v d w x + y r x , y , otherwise

Electrostatic interactions between charged atoms are treated by simple repulsion and attraction according to inverse distance squared. The use of distance squared rather than linear distance encourages the formation of salt bridges in the models. There is a correction for atom-atom distance below the minimum realistic value. The ideal distance between oppositely charged atoms is I_hb-da-dist = 2.75 Ǻ. In the equations below pos is the set of all positively charged atoms and neg is the set of all negatively charged atoms.

E E L E C T R O = ∑ y > x ∈ p o s ∪ n e g Θ ( x ) Θ ( y ) / r x , y * 2 Θ ( x ) = { 1 , if  x ∈ p o s − , 1 if  x ∈ n e g r x , y * = { I e _ d i s t , if   r x , y < I e _ d i s t r x , y , otherwise