Gamess (US) frequently asked questions Part 1: SCF convergence

In spite of the very high quality of the Gamess(US) documentation, the Gamess(US) list is very often flooded with requests from new users regarding the lack of convergence of the SCF procedure. A few words of advice:

When your SCF does not converge,  you should re-run the job including a $guess guess=moread $end line, as well as the complete $VEC group present in the output PUNCH file (usually called <jobname>.dat, and present in you scratch directory).

Addendum:

Whenever you read a $VEC group from a UHF run you must assign NORB in the $GUESS group. An additional problem is that by default the $VEC group only includes the occupied orbitals, and this means that in UHF runs the $VEC group does not include equal numbers of alpha and beta orbitals (e.g., a run with 41 electrons and MULT=2) will have 21 alpha orbitals and 20 beta orbitals. Therefore, if you include

$guess guess=moread NORB=21 $end

Gamess will crash because there are not 21 beta orbitals, and if you input

$guess guess=moread NORB=20 $end

there will be another error, since there are more than 20 alpha orbitals. In these cases, you should check the number of alpha and beta orbitals. Then , copy the coefficients of the extra alpha orbitals to the end of the beta orbitals. In my example above

$guess guess=moread NORB=21 $end

will yield no problems, since the modification of the VEC group yields equal numbers of alpha and beta orbitals. There is also an option to PUNCH every orbital (occupied+virtuals) at every step. In this case, Gamess always punches a full $VEC group, making it very easy to assign NORB as one can simply inspect the output file to learn the number of orbitals. However, this yields gigantic PUNCH files, and may therefore not be feasible.

You should also experiment with changing convergers, damping, etc. Some systems are notoriously hard to converge, and may require several re-iterations of the whole process. 

Advances in peptide chemistry

Protein synthesis is nowadays achieved through molecular biology techniques: the relevant gene is cloned in an appropriate vector, over-expressed with e.g. a poly-histidine tag, and then purified through high affinity chromatography. Peptide chemistry is therefore often forgotten by biochemists, unless we need to order a short customized peptide from a commercial source.
Danishefsky et al. have now combined solid phase peptide synthesis, native chemical ligation and metal-free dethyilation to synthesize a number of analogues of human parathormone. Their strategy afforded native parathormone with higher purity than obtained from commercial sources, as well as pure analogues not achievable by any other means. These analogues were shown to be much more stable (10% decomposition in 7 days) than parathormone ,(>90% loss in 7 days), and to be as active as parathormone when injected to mice.
This is a very interesting work, which should pave the way towards the synthesis of long-lived synthetic peptide hormones, thus potentially decreasing the number of injections needed to control hormone levels in patients suffering from impaired endocrine function.

Drawing can be torture

Drawing complex three-dimensional molecules in two-dimensions can be a real torture. I am glad I have never had to draw anything as convoluted as palhinine A. Check the 3-D structure on the left, and try to draw it in less than 10 minutes in ChemDraw or ChemSketch. Good luck!

QM/MM vs. QM-only studies of large cluster models

How large must a quantum model of an enzyme active site be to achieve optimum results? Proponents of the so-called "cluster model" argue that, most often, good results may be obtained even with small models (< 100 atoms). Fahmi Himo has repeatedly shown that fully including the first layer of aminoacids surrounding the reacting substrate (i.e. to about 150 atoms) yields results that are insensitive to the inclusion of a polarizable-continuum solvent field, and has concluded from these data that such models are sufficient to capture all the relevant enzymatic effexts on catalysis.

Walter Thiel has now published a QM/MM analysis of the reaction mechanism of acetylene hydratase (previously studied by Fahmi Himo using increasingly large QM-only models). Inclusion of the surrounding protein dramatically changed the results for the largest model studied by Himo, due to the absence (in the "cluster model") of two negatively charged phosphate groups adjacent to the active site. Although these charges are quite "shielded" from the active site because of neighbouring positively-charged amino acids, they originate local charge assymmetries that interact differently with the active site during each step of the catalytic cycle. This effect is quite similar to the major influence of the internal protein dipoles on enzyme catalysis expounded by Arieh Warshel, and should be kept in mind by all of us who tend to prefer the QM-only approach: a polarizable-continuum model assumes a homogeneous environment surrounding the QM system, and in proteins "it ain't necessarily so".

An interesting hypothesis on the selection of glucose as major fuel source in neurons

Earlier this year, I wondered why neurons preferentially use glucose as fuel. I have now found an interesting paper by Dave Speijer regarding this problem. He proposes the following reasoning to explain this observation:

reactive oxygen species are generated in large amounts by NADH dehydrogenase (complex I) when the amount of oxidized ubiquinone is limited

generation of large amounts of FADH₂ increases the rate of reduction of ubiquinone, and therefore increases indirectly the amount of harmful radical species generated by NADH dehydrogenase

glucose oxidation generates a much smaller amount of FADH₂ than fatty-acid oxidation. Therefore:

Especially vulnerable cells may be expected to have evolved a preference for glucose.

Incidentally, neurons do seem to lack large amounts of one of the enzymes involved in fatty acid oxidation: thiolase.

The limits of homology modeling

The computational prediction of three-dimensional structures of protein sequences may be performed using a wide variety of techniques, such as homology modeling or threading. In threading, the correct fold is searched for by evaluating the energy of the intended sequence when it is "forced" to adopt each of the known folding patterns. In homology modeling, one looks for a high-similarity protein sequence with experimentally-determined 3D structure, and mutates it in silico until the desired sequence is obtained. Many different programs and web-servers are now available for these tasks, differing among themselves in the forcefields used, alignment algorithms, etc. Performance is usually quite good when templates with similarity >40% are used.

Recently, two small proteins with very high homology (>95%) but widely differing structure have been designed and studied. Starting from a pair of proteins with < 20 % identity and different 3D structures, the authors gradually mutated one sequence into the other, and ended up generating two sequences differing only in one amino acid, but with different folds. Attempts to unravel the precise mechanisms governing the selection of one fold over the other have however been inconclusive, because current molecular dynamics protocols and force fields are not accurate enough to measure the small energy differences involved.

Limitations of PCM

A new paper claims to compute the pKa of nitrous acidium ion from gas phase DFT computations followed by estimation of solvation effects by a Polarizable Continuum Method (PCM). It is true that most often geometries do not change too much when going from gas phase to solution, but I strongly doubt the results are as accurate as they could be: PCM does not include the contribuiton from hydrogen bonds between the solute and the solvent, and I would expect that effect to be quite different in neutral HONO and protonated H₂ONO⁺