QM 004: Two-electron integral manipulation

In the independent-particle model, the electronic energy of an \(N\)-electron system is frequently written as
\begin{equation}E = \sum_{i}^{N} h_i + 1/2 \sum_{i, j} J_{ij} - K_{ij}.\label{eq:energy}\end{equation}
When inserting the expressions for the Coulomb and exchange operators, \(J\) and \(K\), and calculating the variation in energy, \(\delta E \), given by (in physicists notation)
\begin{equation}\delta E = \sum_i \langle \delta \phi_i | h_i | \phi_i \rangle + \langle \phi_i | h_i | \delta \phi_i \rangle + 1/2 \sum_{i, j} \left[ \langle \delta \phi_i \phi_j | \phi_i \phi_j \rangle   + \langle \phi_i \delta \phi_j | \phi_i \phi_j  \rangle  + \phi_i \phi_j | \delta \phi_i \phi_j \rangle + \langle \phi_i \phi_j | \phi_i \delta \phi_j \rangle \\
- \langle\delta\phi_i\phi_j | \phi_j\phi_i\rangle - \langle\phi_i\delta\phi_j | \phi_j\phi_i\rangle -  \langle\phi_i\phi_j | \delta\phi_j\phi_i\rangle - \langle\phi_i\phi_j | \phi_j\delta\phi_i\rangle\right], \label{eq:variation}\end{equation}
it can be found that four pairs of terms in the second sum-expression are equal and thus the factor of \(1/2\) can be cancelled.
For the Coulomb integrals this can be done by considering the definition of the bracket notation
$$ \langle \delta\phi_i(1)\phi_j(2) | \phi_i(1) \phi_j(2)\rangle = \int d{\bf x}_1 d{\bf x}_2 \delta\chi_i^*(1)\chi_j^*(2){\bf r}_{12}^{-1} \chi_i(1)\chi_j(2)$$
and noting that when the orbitals are understood to be real, i.e. \(\phi^* = \phi\), and the ordering of electron labels is considered, the orbitals of an electron \(m \in \{1, 2 \}\) can be swapped such that
$$ \int d{\bf x}_1 d{\bf x}_2 \delta\chi_i(1)\chi_j(2){\bf r}_{12}^{-1} \chi_i(1)\chi_j(2) = \int d{\bf x}_1 d{\bf x}_2 \chi_i(1)\chi_j(2){\bf r}_{12}^{-1} \delta\chi_i(1)\chi_j(2) $$
and thus, in bracket notation,
\begin{equation} \langle \delta\phi_i\phi_j | \phi_i \phi_j \rangle = \langle\phi_i\phi_j | \delta\phi_i \phi_j \rangle \label{eq:manipulation_coulomb}\end{equation}
which we recognize as the the first and third term in the second sum of Eq. \ref{eq:variation}.
For the exchange integrals, a similar expression exists (Szabo, Ostlund, Eq. 2.94)
$$\langle ij | kl \rangle = \langle ji | lk \rangle .$$
To show this, we write the integral explicitly and after first exchanging the dummy variables \(1\), \(2\) and then reordering the orbitals in order to restore the conventional \(1, 2\) order of electrons (keeping "the orbital on the electron"), we obtain
\begin{equation} \int d{\bf x}_1 d{\bf x}_2 \chi_i^*(1)\chi_j^*(2){\bf r}_{12}^{-1} \chi_k(1)\chi_l(2) = \int d{\bf x}_1 d{\bf x}_2 \chi_j^*(1)\chi_i^*(2){\bf r}_{12}^{-1} \chi_l(1)\chi_k(2). \label{eq:manipulation_exchange}\end{equation}
Using this approach
\begin{equation} \langle \delta\phi_i(1)\phi_j(2) | \phi_j(1)\phi_i(2) \rangle = \langle \phi_j(1)\delta\phi_i(2)|\phi_i(1)\phi_j(2) \rangle = \langle \phi_i(1)\phi_j(2)|\phi_j(1)\delta\phi_i(2) \rangle \label{eq:exchange}\end{equation}
where Eq. \ref{eq:manipulation_exchange} was used to obtain the first equality and Eq. \ref{eq:manipulation_coulomb} was used to swap \(\phi_j(1)\) with \(\phi_i(1)\) and \(\delta\phi_i(2)\) with \(\phi_j(2)\) to obtain the second equality.
The first and last integral are recognized as the first and fourth exchange integral of Eq. \ref{eq:variation}.
Applying these operations on the remaining integrals allows to factor out a factor of \(2\) which then cancels with the \(1/2\).

QM 003: Activation free energy in PCM of SN2 reaction

Consider the reaction of chloride ion with methyl iso-cyanide to form methyl chloride and a cyanide ion. What is the activation free energy for this reaction in water?
$$Cl^{-} + H_3C-NC \rightarrow Cl-H_3C + NC^{-}$$
The activation free energy is given by the difference between the free energy of the transition state and the separated species
$$\Delta G = G^{\ddagger} - G^{Reac}.$$
The free energy in solvent is obtained by correcting the energy of the solute, \(E_{solute}\), with the solvation free energy, i.e. the energy required to place the solute in a cavity of solvent
$$G_{solution} = E_{solute} + 1/2 \int_S \sigma(\vec{r}) V(\vec{r}) d\vec{a}.$$
\(G_{solution}\) is output from a PCM calculation as Free Energy in Solvent, meaning the program computes first \(E_{solute}\) and then the solvation free energy (obtained from integrating the electrostatic potential times surface charge density over the surface of the solute).
This value is corrected for free energy contributions from translation, rotation and vibration and ZPE. The translational contribution is adjusted to report the correction for a 1mol/L solution instead of the ideal gas (which is 1mol/24.5L).

The reaction is bimolecular, meaning two particles collide to form one particle (the transition state) which then decays to products. For both reacting particles, the translational free energy correction has to be evaluated (since in solution, other than in gas phase, it is not expected that the particles form a coordinated complex [Vayner et al.]).

From HF/3-21G//HF/3-21G (the only method from which genuine transition state structures could be located), the following values are obtained.

			React 1	React 2	Sum	TS	Products	Activation Energy
PCM			Cl-	H3C-NC	Cl- + H3C-NC	TS	Cl-H3C + NC-	TS - (Cl- + H3C-NC)
HF/3-21G//HF/3-21G	E_elec	Ha	-457.5	-131.2	-588.6	-588.6	-588.6
ZPE	ZPE	Ha	0.0	0.0	0.0	0.0	0.0
E(0K)	E_elec + ZPE	Ha	-457.5	-131.1	-588.6	-588.5	-588.5
E(0K)		kcal/mol	-287068.1	-82277.2	-369345.2	-369298.3	-369298.9	46.9
G_elec		kcal/mol	0.0	0.0	0.0	0.0	0.0
G_trans	24.5l	kcal/mol	-9.4	-9.6	-19.0	-10.1	-10.1	8.9
G_t(1l) = G_t(24.5l) - RTln(24.5)	RT*ln(24.5)	kcal/mol	1.9	1.9	1.9	1.9	1.9
G_trans	1l	kcal/mol	-11.3	-11.5	-22.8	-12.0	-12.0	10.8
G_rot		kcal/mol	0.0	-5.2	-5.2	-6.5	-6.5	-1.2
G_vib		kcal/mol	0.0	30.5	30.5	28.2	28.5	-2.2
G_trv(298)	1mol/24.5l	kcal/mol	-9.4	15.7	6.3	11.7	11.9	5.4
G_trv(298)	1mol/l	kcal/mol	-11.3	13.8	2.5	9.8	10.0	7.3
S_trans	24.5l	cal/(mol*K)	36.6	37.1	73.6	38.9		-34.7
S_trans	1l	cal/(mol*K)	34.7	35.2	69.8	37.0		-32.8
T*S_trans	298K, 24.5l	kcal/(mol)	10.9	11.0	21.9	11.6		-10.4
T*S_trans	298K, 1l	kcal/(mol)	10.3	10.5	20.8	11.0		-9.8
E(0K) + G_trv(298)	24.5l	kcal/mol	-287077.5	-82261.5	-369339.0	-369286.6	-369287.0	52.346
E(0K) + G_trv(298)	1l	kcal/mol	-287079.4	-82263.4	-369342.8	-369288.5	-369288.9	54.250
Free energy in solvent (PCM)		kcal/mol	-287067.9	-82307.8	-369375.7	-369327.1	-369328.7	48.530
Free energy in solvent + ZPE + G_trv	24.5l	kcal/mol	-287077.3	-82261.4	-369338.7	-369286.4	-369286.7	52.348
Free energy in solvent + ZPE + G_trv	1l	kcal/mol	-287079.2	-82263.3	-369342.5	-369288.3	-369288.6	54.252

The increase of the barrier due to the loss in translational entropy is found in the third and fourth column, where '+' indicates the summation of the values from the first and second column (i.e. S\(_{trans, 1M}\) = 69.8 cal/(mol*K) = (34.69 + 35.16) cal/(mol*K), which are obtained from separate calculations of the two species (Cl\(^-\) and H\(_3\)C-NC). This is not the same as calculating S\(_{trans, 1M}\) from a coordinated complex, which consists only from one particle.
The translational entropy (for a 1M state) of the TS however is 37.0 cal/(mol*K). At \(T=298K\), the total entropic change evaluates to (11.59 - 21.95) kcal/mol = -10.35 kcal/mol (for 1mol/24.5L) and from
$$G = H - TS$$
it is seen that this loss of entropy will lead to an increase of the barrier. The correction in translational free energy when going from the gas phase standard state to the solution standard state is moderate (around 2 kcal/mol).

It is seen that given such an approach, the barrier would be estimated to be around 54kcal/mol.
The following points are unclear to me:

1. Should the free energy in solution be corrected for ZPE and free energy contributions from translation, rotation and vibration?
2. Are the free energy corrections of the reactants correct when calculated from a coordinated structure (as opposed to two individual structures in separate calculations)?
3. Which effects not included (apart from electronic structure properties) would be most relevant when extending the model and how could they be obtained? Conformational free energy correction from MD snap shots?

Surprisingly, the activation energy in E(0K) + G\(_{trv, 298, 1L}\) without considering solvation (54.25 kcal/mol) is the same as the activation free energy in solvent. Is this a coincidence or does this make sense or not?

PyMOL 006: Removing a bond between two atoms

PyMOL sometimes displays bonds which are not correct and can be confusing. To remove a bond between to atoms, first set the selection argument to "Atoms", then select the first atom involved in the bond and rename its selection identifier to, e.g., "atm1". Unselect the atom, select the second atom of the bond, rename the selection identifier to "atm2".

Then enter the following at the PyMOL prompt:

cmd.unbond("(atm1)", "(atm2)")

Be sure to include all parentheses and quotation marks.

Python 005: The so called "sliding window"

Say you find yourself required to locate the line number of the last atom of residue '2' in this PDB (following the PDB internal numbering scheme, of course we assume you don't know a priori its called H01):

ATOM      1  N   ALA     2      -0.677  -1.230  -0.491  1.00  0.00           N
ATOM      2  CA  ALA     2      -0.001   0.064  -0.491  1.00  0.00           C
ATOM      3  C   ALA     2       1.499  -0.110  -0.491  1.00  0.00           C
ATOM      4  O   ALA     2       2.065  -0.922   0.251  1.00  0.00           O
ATOM      5  CB  ALA     2      -0.509   0.856   0.727  1.00  0.00           C
ATOM      6  H   ALA     2      -0.131  -2.162  -0.491  1.00  0.00           H
ATOM      7  HA  ALA     2      -0.269   0.603  -1.418  1.00  0.00           H
ATOM      8 1HB  ALA     2      -1.605   1.006   0.691  1.00  0.00           H
ATOM      9 2HB  ALA     2      -0.285   0.342   1.681  1.00  0.00           H
ATOM     10 3HB  ALA     2      -0.053   1.861   0.784  1.00  0.00           H
ATOM     11  H01 ALA     2      -1.261  -1.244  -1.314  1.00  0.00           H
ATOM     12  N   ASN     3       2.311   0.711  -1.400  1.00  0.00           N
ATOM     13  CA  ASN     3       3.700   0.321  -1.173  1.00  0.00           C
ATOM     14  C   ASN     3       4.606   1.530  -1.169  1.00  0.00           C
ATOM     15  O   ASN     3       4.515   2.417  -2.017  1.00  0.00           O
ATOM     16  CB  ASN     3       4.140  -0.693  -2.267  1.00  0.00           C
ATOM     17  CG  ASN     3       3.309  -1.974  -2.400  1.00  0.00           C
ATOM     18  ND2 ASN     3       3.497  -2.939  -1.540  1.00  0.00           N
ATOM     19  OD1 ASN     3       2.462  -2.114  -3.270  1.00  0.00           O
ATOM     20  H   ASN     3       1.953   1.454  -2.096  1.00  0.00           H
ATOM     21  HA  ASN     3       3.774  -0.146  -0.174  1.00  0.00           H
ATOM     22 2HB  ASN     3       5.198  -0.981  -2.124  1.00  0.00           H
ATOM     23 3HB  ASN     3       4.123  -0.199  -3.258  1.00  0.00           H
ATOM     24 1HD2 ASN     3       2.793  -3.677  -1.623  1.00  0.00           H
ATOM     25 2HD2 ASN     3       4.145  -2.756  -0.775  1.00  0.00           H
ATOM     26  N   PRO     4       5.629   1.649  -0.120  1.00  0.00           N
ATOM     27  CA  PRO     4       6.360   2.889  -0.340  1.00  0.00           C
ATOM     28  C   PRO     4       6.031   3.505  -1.710  1.00  0.00           C
ATOM     29  O   PRO     4       5.228   2.956  -2.465  1.00  0.00           O
ATOM     30  CB  PRO     4       7.819   2.420  -0.227  1.00  0.00           C
ATOM     31  CG  PRO     4       7.769   0.962  -0.705  1.00  0.00           C
ATOM     32  CD  PRO     4       6.440   0.451  -0.143  1.00  0.00           C
ATOM     33  HA  PRO     4       6.122   3.617   0.439  1.00  0.00           H
ATOM     34 2HB  PRO     4       8.119   2.445   0.823  1.00  0.00           H
ATOM     35 3HB  PRO     4       8.509   3.035  -0.810  1.00  0.00           H
ATOM     36 2HG  PRO     4       8.619   0.375  -0.350  1.00  0.00           H
ATOM     37 3HG  PRO     4       7.737   0.935  -1.796  1.00  0.00           H
ATOM     38 2HD  PRO     4       6.564   0.075   0.875  1.00  0.00           H
ATOM     39 3HD  PRO     4       5.991  -0.316  -0.778  1.00  0.00           H
ATOM     40  H01 PRO     4       6.512   4.435  -2.015  1.00  0.00           H

As you can see, it's the Atom 11 (H01) of an Alanine. How would you do that? Choose a for-loop, choose a nicely indented block of if and else statements. Choose a number of different counters, each with their very own meaningful name. Choose endless meetings with IndexErrors, comparisons, type castings and list slicings. I say no, I choose a generator. In fact, a pair wise sliding window.

def make_sliding_window(val):
    ln_cnt = 0
    while ln_cnt < len(val):
        yield ln_cnt, ln_cnt-1
        ln_cnt += 1 

The 'yield' statement characterizes this as a generator. What does that mean? First, a generator is *not* a list, or array or something like that. It's a function. With special properties. When called through 'next(generator_function)', it executes until it reaches a 'yield' statement, from where it returns. When called again (with a new 'next(generator_function)') it resumes from where it left of. Yes, and it returns a generator object, when you call it as such (no 'next'), but then it does not execute.
So lets instantiate the object.

res=make_sliding_window(val)

'res' is now the generator. Only when you call 'next(res)', the generator executes, up until the 'yield' statement, from where it returns. In addition, it can be used directly in a for-loop. The for-loop takes care of implicitly calling the 'next' method. Here's an example:

def get_line_num_last_atom_res_2(val):
    # Form a generator 'res'.
    res = make_sliding_window(val)
    for a, b in res:
        if val[a][25] == '3':
            return val[b], a

print get_line_num_last_atom_res_2(open('randompdbfile.pdb', 'r').readlines())

As you can see, we're now using it as-if it were a list. But notice, if you think it's a list, you'll not understand what is happening. It has a 'next' method, that's the important point and it's being called, even if we dont see it explicitly. So anyway, upon reaching 'yield', we get apparently two line numbers back, where it's two times the same one, but just offset by one. Now we can easily check if the upcoming atom already is part of the next residue, and if not, return the current atom. And its line number.

Python 004: My first closure

Depending on the presence of a number of files, a couple of functions had to be called. But depending on the kind of file, a different function was required. However, for each file, the logic to determine if the corresponding function should be called or not is the same. For technical reasons, the respective functions are in different modules (executed through some Ajax).

So it seemed to me I would have to write the logic part specifically for each case, i.e. for each case where a file is available or not.

Something as below:

if file_exists('-ini'):
    function_one()
if not file_exists('-ini'):
    generate_ini_file()
    function_one()
if file_exits('-fin'):
    function_two()
if not file_exists('-fin'):
    generate_fin_file()
    function_two()

and so on. However, I remembered, in Python everything is an object. Even functions. So, why not prepare a function, which takes care of the logic and the execution of the function which I want to apply to each of the corresponding files. But then, where should the applied function come from? Exactly, from a closure. This basically is a function which generates functions at runtime. What kind of arguments does it take? Precicely, the function you actually want to call. Huh? Yes, of course the applied function needs to be defined somewhere, but thats another story, lets assume we have it. It could look something like this:

def translate_com(target):
    """Calling VMD and returning the output directly."""
    vmd_src = 'mol load pdb %s-fix.pdb\n' % (bio_lib.pdb_base_path + target)\
            + 'set HOME ' + bio_lib.vmd_base_path + '\n'\
            + 'molinfo 0 get center\n'\
            + 'env\n'
    comp = Popen([bio_lib.vmd_cmd_path, '-dispdev', 'text', '-eofexit'],
                  stdout=PIPE, stdin=PIPE, stderr=PIPE, shell=False)
    comp.stdin.write(vmd_src)
    vmd_com_result = comp.communicate()[0]

    for line in vmd_com_result.split('\n'):
        if len(line) > 0 and line[0] == '{':
            com_xyz = [float(i) for i in line[1:-1].split()]
    com_xyz = [i*(-1) for i in com_xyz] 
    # Calling the VMD reorientation script.
    reo_src = bio_lib.get_reo_src(com_xyz, target)
    vmdp = Popen([bio_lib.vmd_cmd_path, '-dispdev', 'text', '-eofexit'],
                  stdin=PIPE, stdout=PIPE, stderr=PIPE, shell=False)
    vmdp.stdin.write(reo_src)
    print vmdp.communicate()[0]

It does something, other functions are doing similar things. This is, e.g., the 'function_two' above. So for each of these, there would be a logic block as this:

def check_availability_and_translate(target):
    # Use upload.
    if uploaded:
        # Force overwrite.
        if overwrite:
            translate_com(target)
        # Does not exist.
        if not os.path.exists(bio_lib.pdb_base_path + target + '-reo.pdb'):
            translate_com(target)
    # Download.
    else:
        # Use existing.
        if os.path.exists(bio_lib.pdb_base_path + target + '-reo.pdb'):
            print "Structure available." 
        # Fixed structure not available.
        else:
            reo.translate_com(target)

In the above, there's also a bit of user interaction coming in, but that's not the point. The point is, all this is the same for each of my functions. So, lets generate a function that returns this logic, but takes the 'translate_com' function as an argument and calls it whenever *a* function is called:

def closure(state):
    """
        'state': '-ini', '-fin', '-reo'
    """
    def check_availability_and_apply(target, uploaded, overwrite, func):
        # Use upload.
        if uploaded:
            print "Uploaded"
            # Force overwrite.
            if overwrite:
                print "Overwriting"
                func(target)
            # Does not exist.
            if not os.path.exists(pdb_base_path + target + state):
                print "Not found, fixing."
                func(target)
            # Use existing.
        # Download.
        else:
            # Use existing.
            if os.path.exists(pdb_base_path + target + state):
                print "Structure available."
            # Structure not available.
            else:
                func(target)
    return check_availability_and_apply

The fun is in the last line. As one can see, a function is returned which was built within the closure. So how can I now call (or maybe better "apply") the 'translate_com' function?

closure('-ini.pdb')(target, uploaded, overwrite, translate_com)

And how about another function?

closure('-fix.pdb')(target, uploaded, overwrite, fix_pdb)

This is really cool, because now file checks are in one place and application is in a separate place. So only one place to edit.
I would really like to hear my boss asking me now: "So how many lines of code did you write yesterday?" Then I could answer: "Well, actually like -200."

Vim 006: Transform Block to Lines

A block of text can be split up into individual lines by hitting "gww" within the respective block. The reverse is possible when highlighting a number of lines in visual mode and hitting "J", all lines are combined into a single, continuous block.

QM 002: Comparing PCM and APBS

I compare the solvation energy of the Born ion using "quantum mechanics" and electrostatics. The system of interest is an $H^+$ ion, so the reason for the quotation marks around quantum mechanics is that there is no electron, so how can there be any quantum mechanics? Nevermind, I use the term because the PCM model implemented in the GAMESS program is compared with the solution to the Poisson-Boltzmann equation as can be found using the Adaptive Poisson-Boltzmann Solver (APBS) program.
The free energy of solution can be obtained from the equation

$$G_{solution} = E_{solute} + \frac{1}{2} \int_s d\vec{r} \sigma (\vec{r}) V(\vec{r})$$

First, the solution obtained from GAMESS/PCM is presented. The first term on the right side in the equation for the free energy in solution turns out to be exacltly zero because there is no electron in the system. Next, the second term on the right side is evaluated using the below PCM input for GAMESS.

1  $BASIS gbasis=STO ngauss=2 $end
2  $contrl mult=1 scftyp=rhf icharg=1 $end
3  $pcm solvnt=water $end
4  $DATA
5
6 C1
7  H 1.0 0.0 0.0 0.0
8  $END

At the end of the calculation, the following is printed to the log file:

 ----------------------------------------------
-------   RESULTS OF PCM CALCULATION   -------
----------------------------------------------

FREE ENERGY IN SOLVENT =                            =       -0.1813647781 A.U.
INTERNAL ENERGY IN SOLVENT =                        =        0.0000000000 A.U.
DELTA INTERNAL ENERGY =                             =        0.0000000000 A.U.
ELECTROSTATIC INTERACTION                           =       -0.1813647781 A.U.
PIEROTTI CAVITATION ENERGY                          =        0.0000000000 A.U.
DISPERSION FREE ENERGY                              =        0.0000000000 A.U.
REPULSION FREE ENERGY                               =        0.0000000000 A.U.
TOTAL INTERACTION (DELTA + ES + CAV + DISP + REP)   =       -0.1813647781 A.U.
TOTAL FREE ENERGY IN SOLVENT                        =       -0.1813647781 A.U.

The solvation energy (the energy released when transfering the ion from vacuum into the solvent) is now calculated as

$$(-0.1813 - 0.0)Ha*2526.5kJ/(mol*Ha)*\frac{1}{2} = -229.03kJ/mol.$$

Next, we us APBS to solve the problem. This requires a similar approach, involving the calculation of the gas phase and solution phase electrostatic energy. The below input for APBS carries out both calculations in the same run.

# READ IN MOLECULES
# Input copied from poissonboltzmann.org
read
mol pqr ion.pqr
end

# ...................................................................
# Electrostatics calculation on the solvated state
elec name solv        
mg-manual             # Specify the mode for APBS to run
dime 97 97 97         # The grid dimensions
#nlev 4               # Multigrid level parameter, not available in APBS 1.3
grid 0.33 0.33 0.33   # Grid spacing
gcent mol 1           # Center the grid on molecule 1
mol 1                 # Perform the calculation on molecule 1
lpbe                  # Solve the linearized Poisson-Boltzmann equation
bcfl mdh              # Use all multipole moments when calculating the potential
pdie 1.0              # Solute dielectric
sdie 78.54            # Solvent dielectric
chgm spl2             # Spline-based discretization of the delta functions
srfm mol              # Molecular surface definition
srad 1.4              # Solvent probe radius (for molecular surface)
swin 0.3              # Solvent surface spline window (not used here)
sdens 10.0            # Sphere density for accessibility object
temp 298.15           # Temperature
calcenergy total      # Calculate energies
calcforce no          # Do not calculate forces
end

# ...................................................................
# Calculate potential for reference (vacuum) state
elec name ref
mg-manual
dime 97 97 97
#nlev 4
grid 0.33 0.33 0.33
gcent mol 1
mol 1
lpbe
bcfl mdh
pdie 1.0
sdie 1.0
chgm spl2
srfm mol
srad 1.4
swin 0.3
sdens 10.0
temp 298.15
calcenergy total
calcforce no
end

# ...................................................................
# Calculate solvation energy
print elecEnergy solv - ref end
quit

The ion is defined by a PDB formatted input file:

ATOM      1  I   ION     1       0.000   0.000  0.000  1.00  3.0

The charge of the ion is defined in the second to last column, and the radius in the last one. The calculated energy of solvation is the energy of the solvated ion minus the energy of the gas phase ion:

print energy 1 (solv) - 2 (ref) end
Local net energy (PE 0) = -2.295878286327E+02 kJ/mol
Global net ELEC energy = -2.295878286327E+02 kJ/mol

As can be seen, the values from PCM and Poisson-Boltzmann equation are very close.

The problem can also be solved analytically. It involves the calculation of an integral over the surface of a sphere, which can be done as follows. First, the expression for the surface charge densitity $\sigma (r)$ needs to be evaluated. From the Poisson equation, it follows that $\sigma (\vec{r}) = -\frac{1}{4\pi}\left(1- \frac{1}{\epsilon}\right)Q\frac{1}{r^2}$ and using a potential $V(r) = Q\frac{1}{r}$, the integral to solve is thus

$$-\frac{1}{2}\int_s d\vec{r} \sigma (\vec{r}) V(\vec{r}) = -\frac{1}{8\pi} \left(1- \frac{1}{\epsilon}\right)Q^2\int_s d\vec{r} \frac{1}{r^3}$$

After realizing that the integration over the surface of the sphere requires the insertion of the Jacobi-Determinant $r^2 \sin \theta$, the new integration limits can be applied, the integral becomes independent of radius (setting $r \rightarrow R$) and one obtaines

$$-\frac{1}{8\pi}\left(1-\frac{1}{\epsilon}\right)Q^2\frac{1}{R^3}R^2\int_0^{2\pi} d\phi \int_0^{\pi} d\theta \sin \theta$$

which can be evaluated to

$$-\frac{1}{2}\left(1-\frac{1}{\epsilon}\right)\frac{Q^2}{R}.$$

From this it is seen that the solvation energy is depending more on the charge of the ion, than on its radius or the dielectric constant.

QM 001: Getting the barrier for rotation in H2C=CH2, correctly

Upcoming post on the calculation of the energy barrier around the sigma bond in H2C=CH2.

Bash if statement

The Bash if statement is a bit funny. Say for every files-i.txt, there should be a files-i.dat.
This example can be used to filter out files that dont have the corresponding .dat file present.

for i in files-*txt
do
    if [ ! -e ${i/txt/dat} ]
    then
         echo $i
    fi
done

Note the spaces after and before the opening and closing hard bracket. Also note that "then" is on its own line.

PyMOL 005: Partial structure optimization

PyMOL is amazing. Say you have a structure where, eg., water molecules are too close to your protein atoms. Maybe you want to optimize the water molecules (only the water molecules) by some crude approach before submitting to more sophisticated theory.
Here's how you can do it at the PyMOL prompt. You loaded the structure as "enz". The waters in "enz" all have the same id 500. So typing

select wat, enz/PROT/A/500/

should select your waters. "PROT" is some sort of segment and "A" is the chain of the structure, these depend and maybe need some adjustment in your system. Alternatively, try

select wat, enz///500/

Having selected the waters (if nothing else then manually), enter the following commands

cmd.protect('(not wat)')
cmd.sculpt_activate('enz')
cmd.sculpt_iterate('enz', cycles=5000)
cmd.sculpt_deactivate('enz')
cmd.deprotect()

The number of cycles can be adjusted. Also, it is possible to optimize locally only a side chain of one residue in the protein. Then you would enter

cmd.protect('(not enz///105/) or name C+CA+O+N+OXT')
cmd.sculpt_activate('enz')
[same commands like above]

Dont forget the trailing forward slash in the selection commands.

Displaying the Path in Mac OS X

A great little trick for Mac OS X

$ defaults write com.apple.finder _FXShowPosixPathInTitle -bool YES
$ osascript -e 'tell app "Finder" to quit'