28 February 2014

559. Briefly: Nudge Elastic Band in NWChem

I'm currently exploring a method called Nudge Elastic Band to find the minimum energy pathway in a reaction involving a large metal cluster. While the NWChem documentation isn't bad, it could be clearer for happy amateurs like myself, which is the impetus for this post.

As usual, this is how I understand it -- which may be wrong.

The NEB method is described here.


There are a number of ways of modelling reaction pathways computationally.

Brute force PES scan (g09 or nwchem)
The fastest, cheapest, least accurate one would be to make a simple Potential Energy Surface (PES) scan by keeping all coordinates frozen with the exception of the ones directly involved in a reaction
e.g. for Me-OH -> Me + OH you'd simply vary the C-O distance.

If you don't let the structures relax for each change in bond length you will overestimate the activation energy. If you do let the structures relax you run the risk of every single snapshot relaxing back to either pure methanol or pure Me and OH, rather than giving you a series of transitional energies.

Transition state search e.g. QST2 (g09)
A somewhat more sophisticated method is to try to identify a proper transition state (TS). You can do this by making a lucky guess and try to do a transition state optimisation. You can improve your chances by using e.g. QST2 (or QST3, which also uses a specific TS guess), which takes the product and the starting material and interpolates the coordinates to generate a transition state guess. This has worked pretty well for me for cases where a simple transition state is expected (i.e. a hydride transfer). You can then model the reaction path by doing an IRC calc.

Minimum energy pathway methods (chain-of-states) e.g. Nudge Elastic Band (nwchem)

A more sophisticated method is to use a chain-of-state method such as Nudge Elastic Band (NEB). There's a nice page about it here: http://theory.cm.utexas.edu/henkelman/research/saddle/

NEB generates a series of structures based on the starting point and the product. These are initially generated a simple interpolations between the coordinates of the starting point and the end point -- this is very similar to a brute for PES scan and gives an unnaturally high transition energy. However, the individual structures are then optimised in passes.

10 passes. Structural convergence is not achieved.

Not clear? Well, here's a rough algorithm:

1. Generate a series of linearly interpolated structures (beads) based on the Starting Point (SP) and End Point (EP) structures. In the simplest case, for N structures, each structure n is SP+(n/N)*(EP-SP).

Now, typically you're more interested in structures in the middle than at the ends (i.e. to get better resolution for the transition) and that is manipulated using the spring constant, kbeads. See figure 3 in http://theory.cm.utexas.edu/henkelman/pubs/jonsson98_385.pdf for an example of varying the spring constant. 1.0 seems to be reasonably safe.

The number of structures is set by nbeads. Make sure to use a reasonable number so that you get enough resolution.

2.  Each structure thus generated is submitted to a single geometry optimisation step i.e. it's not optimised until convergence, but only takes a single step along the energy gradient. The speed at which it is decending along the energy gradient is set using the stepsize.

The largest possible stepsize will give the fastest descent (which is good), but too large a stepsize may cause issues due to the atoms being moved too far (which is bad) and in turn cause SCF convergence failures due to the structures becoming unreasonable. So try 1.0 and hope for the best.

Rule of thumb: if you decrease the stepsize by a factor of 10 from 1.0 to 0.1, you should increase maxiter by a factor of 10.

3. The energies of the new structures are calculated and a reaction profile is generated.

4. Repeat step 2 and 3 until structural convergence or maxiter steps.

5. If you're lucky your calculation will end due to structural convergence. If not, and if you want to restart your calculation you can read in the last set of structures using xyz_path.

But make sure that you do print the intermediate structures in the first place by using print_shift 1.

558. More options for PDF annotation: Master PDF Editor and I, Librarian

This is an update to my previous posts (here and here) about annotating PDFs on Linux.

Master PDF Editor
The linux version of this software is free and closed source. It will not run on Debian Wheezy, as it requires glibc >= 2.14, whereas Wheezy has version 2.13. You can either pull in a newer libc from testing, or simply use testing (Jessie). See here for how to use glibc > 2.13 on wheezy: http://verahill.blogspot.com.au/2014/03/562-pulling-in-glibc-214-from-testing.html

So I tried it on Jessie.

While it almost works, I can't prevent the program from screwing up the fonts in the process. Also, the annotations don't show up in acroread for some reason. Saving as PDF/A didn't make any positive difference.

The annotation shows up in evince -- but the font's weird now

Same goes for Xournal

In acrobat reader before editing with Master PDF Editor
in acrobat reader after adding an annotation. The annotation doesn't show up, but the fonts are screwed up.

So we're almost there -- Master PDF Editor does everything I want in terms of PDF annotation, but at least on linux there are issues with the fonts in the resulting PDFs.

I, Librarian
I, Librarian is both free and open source. It claims to be web-based -- which is true but can be misunderstood. What it means is that it is browser-based and runs a server on your computer. My first thought when I see 'web-based' is that I'm handing over my data to someone else, but luckily that's not the case here.

I installed the .deb file meant for ubuntu.

Using I, Librarian is quite straight forward, but I could not see the annotations in any other program that I tried, which makes it of little use for me -- I make annotations in galley drafts for journal editors, and for students to give them feedback on latex documents.

Either way, to start, just go to

Make sure to edit the settings to make I, Librarian use the internal pdf viewer in order for editing to work.

No matter what I did, I could not export my pinned annotations though. They did not show up in either evince or acroread.

19 February 2014

557. Briefly: Drawing Molecules on Linux: ACD/ChemSketch in Wine

This is another post where I'm cheating a bit -- this is a windows program running in wine.

Let me explain. Sometimes you absolutely need to get a task done -- I use linux at work and I can't not do something just because I can't write/find a piece of software that does it for me. Science is the goal, and the route there is immaterial. So sometimes you need to be pragmatic about your software choices.

I much prefer to use free and open source software since it allows me to keep a copy of the source code which I can recompile in case the prebuilt binary suddenly won't work with a future version of debian. Having the source makes me feel more secure in the assertion that it's worth my time learning how to use that piece of software.

Second to that, I prefer native linux software -- although unfortunately 'native' here often means 'written in java', and -- while not knowing too much about java -- from a user point of view java software tends to be comparatively slow to load and run. Certainly it appears slower than a comparable C/C++ program.

As a last resort, I can accept having to run a windows binary in wine. It doesn't make me happy, and often there are small, niggling issues associated with it -- but if it can get the job done, so be it.

Note that I'm unwilling to actually run a native windows program on a native windows installation -- one has to have some standards...

Either way, this is why I look at windows programs as well for drawing structures and processing NMR spectra.

So I recently tested marvinsketch (linux), bkchem (linux), easychem (linux), chemtool (linux), gchempaint (linux), ISISDraw (wine) and ChemSketch (wine). Out of those I prefer MarvinSketch. However, I'm still exploring the alternatives.

ACD/ChemSketch for Linux is really just an older (2006) version of ACD/LABS free ChemSketch (the windows version is from 2013). To get it, go to http://www.acdlabs.com/resources/freeware/index.php and click on Download. I tried the version 'for linux', and not the windows one. I'm making the presumption that the newer binary won't play well with wine.

You don't need to be an academic to qualify for the free version.

wine ~/Downloads/chemsk50.exe

Once it's installed you can explore ChemSketch. Overall it's fairly similar to e.g. ISISDraw.

18 February 2014

556. Briefly: Drawing molecules on linux: ghempaint, chemtool and easychem

chemtool, easychem and gchempaint are in the Wheezy repos.
See also the following posts for ISISDraw, MarvinSketch and BKChem.

I didn't get along with this at all. I managed to import a few structures from the templates, but e.g. the clean molecule function didn't actually work, although it's there. Selecting, deleting etc. were all a bit non-intuitive (to me -- this is obviously very subjective)

I can't find any templates. Drawing becomes a bit tedious since I didn't find an easy way of converting an existing bond from e.g. single to double.

It's pretty straightforward to use and I like it the best of the three. However, it's still very sparse even compared to e.g. bkchem. It certainly works if you quickly need to make a simple structure. There's no clean function.

555. Very Briefly: Drawing molecules on linux: bkchem

I can't evaluate drawing programs on linux without giving bkchem a look. In fact, up till now this has been my preferred drawing application on linux, in particular since it uses .svg as it's native file format i.e. it's very easy to post-process drawing and turn them into proper .eps files for inclusion in latex documents.

bkchem is in the debian repos:
sudo apt-get install bkchem

It's a very basic program, although looking through the options it isn't too feature sparse either.

However, there's one important feature that ISISDraw and MarvinSketch have that I can't find in bkchem: it's lacking a function to clean molecules i.e. to make them look nice.

554. Briefly: Drawing molecules on linux: Marvinsketch

I don't remember how I installed this, but marvinsketch is by far the fanciest molecular drawing application on linux. It's written in java, which means that it runs natively (well, as natively as something like java can run).

Note: it's free for academia only and you will need to register.

To get MarvinBeans, go to https://www.chemaxon.com/ and register as an academic. Then go to http://www.chemaxon.com/download/marvin-suite/ and download the 64 bit linux installer with jre.

Presumable you can install it by doing
sh marvinbeans-6.1.0-linux_with_jre64.sh

in the download directory.
Once it's installed create a file called ~/.local/share/applications/marvinsketch.desktop:
[Desktop Entry] Name=MarvinSketch GenericName=MarvinSketch Comment=Software for drawing molecules Exec=sh /home/me/ChemAxon/MarvinBeans/MarvinSketch Terminal=false Type=Application Categories=Science Version=6.1

You should now be able to find it in GNOME 3.

There are a lot of similarities between ISISDraw and MarvinSketch -- and that's a good thing since it allows people like myself to draw molecules quickly without thinking too much about it.
A Blank slate

Running Chem Inspector

It didn't like my wedges

I turned the structure into 3D and rotated it

Elemental analysis

553. Briefly: Drawing molecules on Debian Wheezy: ISISDraw under wine (very briefly)

I'm cheating here: this isn't a program for linux. Instead we're running a windows program under wine, which is not ideal. On the other hand, it's so easy to do that it's worth exploring as an option to 'pure' linux offerings.

I'll be posting about native linux offerings later.

Either way, I've kept a copy of ISISDraw 2.5 around since the early 2000s. Luckily, other people are mirroring the installation file for it as well. See e.g. http://fc.smcdsb.on.ca/~rszerminski/addl_files_for_web/DRAW24.EXE for ISISDraw 2.4

 You'll need to install wine:i386 and for that you'll need to enable multiarch if you're on amd64:
sudo dpkg --add-architecture
sudo apt-get update
sudo apt-get install wine:i386
wget http://fc.smcdsb.on.ca/~rszerminski/addl_files_for_web/DRAW24.EXE
wine DRAW24.EXE

As always, if you have problems with e.g. weird fonts, use winetricks to install at least corefonts, or ideally allfonts:
wget http://winetricks.org/winetricks
sudo mv winetricks /usr/bin/winetricks
chmod +x /usr/bin/winetricks 
winetricks allfonts

Then run ISISDraw:
A blank slate
 ISISDraw has a lot of options, including templates and journal settings.
Chem Inspector
 There's a bug under wine: to get the menus for the buttons on the left side to show up, click on the button on the left panel, and then move the mouse left while still holding the button down.
To get the menus to work, click and while still holding the left button down move right
 ISISDraw is not going to blow you away, but it's free and it works under wine.

15 February 2014

552. Very briefly: Enthalpy of correction and different number of reactants and products in nwchem and gaussian

Those who are well-versed in the computational arts won't get much out of this post. On the other hand, happy amateurs like myself might find it a bit more useful as it clarifies something that's been bothering  me.

Short version: Hcorr in Nwchem and Gaussian include PV.

As usual, I am not an expert in either computational or theoretical chemistry. I try to use it as a tool, and I try to use it as well as I can. But I am not an authority. Also, if you consult older posts on the blog you'll find plenty of examples of me misunderstanding basic computational concepts (with 550+ posts it's difficult to go back and erase all the embarrassing little gems).

The background
I had a bit of a fright the other day. I'm currently working on computationally characterising a system which undergoes a reaction that can lead to a large number of isomers. Only one of them is experimentally observed, and so it was of interest to see whether this is the thermodynamically favoured product as predict by reasonably cheap computational methods.

Because DFT calculations like these are based on gas phase reactions (even if you use a solvation model) the free energy that you get is based on the standard conditions for gas phase corrections i.e. 1 bar of partial pressure of each reactant.

If you want to calculate the free energy of reaction in solution you need to use concentrations instead. As you will see, you'll only have to worry about this for reactions that involve an unequal number of reactants and products. Normally your best results will be obtained by using isodesmic reaction schemes anyway, which is a great way of avoiding this. However, if you do have unequal numbers of reactants and products you /must/ correct for it when making solution phase predictions.

A gas at 1 bar of pressure is ca
V=n*R*T/P= 1*8.314*298.15/101350=0.024458 m^3= 24.46 litres

So for an example reaction like this:
A+B -> C


G_SATP=G(gas phase) + R*T*ln Q
=G(gas phase)+ R*T*ln([C]/([A]*[B])
=G(gas phase)+ R*T*ln((1/24.46)/((1/24.46)*(1/24.46)))
=G(gas phase)+1.89 kcal/mol 

In other words: you can't ignore the standard state change when doing solution phase calculations. This is obviously of extra importance in pH calculations, which are notoriously tricky.

So knowing the need for standard state corrections I was a bit paranoid about how to treat the reaction enthalpy and came across this document: http://chemistry.illinoisstate.edu/standard/che38037/handouts/380.37assign3.pdf

On page, equation 2 states that for Gaussian

(where Ee is the electronic energy) which is an expression that doesn't include Δ(PV). In that case
where Δn is the difference in number of moles of products and reactants.

Consulting the excellent Gaussian thermochemistry whitepaper (http://www.gaussian.com/g_whitepap/thermo.htm) offers the following:
The Gibbs free energy includes Δ(PV) , so when it's applied to calculating ΔG for a reaction, ΔNRT=ΔPV is already included. This means that ΔG will be computed correctly when the number of moles of gas changes during the course of a reaction.
[Note that H=Ee+Hcorr=Ee+Etot+kBT]

At a first glance, kBT isn't equivalent to RT, but in fact kB=R/NA -- in the words of Wikipedia: "[The gas constant R] is equivalent to the Boltzmann constant, but expressed in units of energy".

In other words, Δ(PV) is already accounted for in Hcorr in gaussian.

Somewhat clearer: the Freq() page, http://www.gaussian.com/g_tech/g_ur/k_freq.htm, on the gaussian website now states
Sum of electronic and zero-point Energies=-527.492585  E0=Eelec+ZPE
Sum of electronic and thermal Energies= -527.489751     E= E0+ Evib+ Erot+Etrans
Sum of electronic and thermal Enthalpies=-527.488807  H=E+RT

We test that claim:
E=-527.489751 Hartree
RT=(1.987*298.15/1000)/627.503=9.4410e-04 Hartree

It holds up. 

A similar example from an nwchem calculation:
Zero-Point correction to Energy = 120.416 kcal/mol ( 0.191895 au)
Thermal correction to Energy = 127.114 kcal/mol ( 0.202569 au)
Thermal correction to Enthalpy = 127.706 kcal/mol ( 0.203513 au)
(0.203513-0.202569)=9.4400e-04 Hartree(1.987*298.15/1000)/627.503

Nwchem also includes Δ(PV) in the thermal correction to enthalpy.

05 February 2014

551. Very Briefly: Getting ECCE to work with Gaussian 09 (G09) part 3: Energy gradients (EGRADVEC)

This is a really simple fix to make ECCE return the energy gradient vectors for geometry optimisation jobs in G09:

Edit apps/scripts/parsers/gaussian-03.egradvec and comment out lines 38 (if () {... ) and 70 (}):

diff --git a/ecce-v7.0/scripts/parsers/gaussian-03.egradvec b/ecce-v7.0/scripts/parsers/gaussian-03.egradvec
index a2f5633..4a4b691 100755
--- a/ecce-v7.0/scripts/parsers/gaussian-03.egradvec
+++ b/ecce-v7.0/scripts/parsers/gaussian-03.egradvec
@@ -35,7 +35,7 @@ $| = 1;
 #Don't parse gradients for runtypes with geometry optimization
 #since we can't match gradients with correct geometry.
-if (!($runtype =~ /Geo/i)) {
+#if (!($runtype =~ /Geo/i)) {
   $natom = 0;
   while () {
     if (/-----/) { last; }
@@ -67,4 +67,4 @@ if (!($runtype =~ /Geo/i)) {
   print "units:\nHartree/Bohr\n";
   print "END\n";

550. Briefly: Getting ECCE to work with Gaussian 09 (G09) part 2: Geometry trace when using NoSymm

I've made a couple of posts about how to modify ECCE to work better in general, but in particular with Gaussian 09  (G09).

See here for some of the previous posts:
Adding new basis sets, part 1: http://verahill.blogspot.com.au/2013/06/455-adding-nwchem-basis-sets-to-ecce.html
Adding new basis sets, part 2:  http://verahill.blogspot.com.au/2013/06/456-adding-nwchem-basis-sets-to-ecce.html
Adding new exchange/correlation functionals:
Adding new options to ECCE
Getting ECCE to read frequency calc output from G09 calcs:

Either way, I've noticed that when I'm using the NoSymm option during optimisation ECCE fails to import the geometries from the geometry steps.
No geometry trace option

What is even worse is that if you use an optimisation job like that as the basis for another job by selecting Duplicate Setup with Last Geometry ECCE will in fact use the starting geometry -- not the optimised one.

Luckily the fix is simple -- edit apps/parsers/gaussian-03.desc and add the following lines:
[GEOMTRACE] Script=gaussian-03.geomtrace Begin=Z-Matrix orientation\: Skip=5 Frequency=all End=-------------------------------------- [END]

You can now 'Reconnect job monitoring' also to your old jobs if you still have access to the run directory. Either way, the view will now be a sweet one:
Et voila: A geometry trace

03 February 2014

549. Very briefly: Unexpected comment spam

I get about 2-3 spam comments on my blog per day. Typically they advertise money making schemes, email marketing schemes etc., and look something like this:
This site was... how do I say it? Relevant!! Finally I have found something that helped me. Kudos! Also visit my web site: minecraft 1.0

Today I got an odd one though and I'll post it here even though it kind of works in their favour:
 Anonymous has left a new comment on your post "341. Upgrading/installing BankID on 64 bit linux":

That was enjoyable to read, thanks for posting it.

My site: Bing (www.bing.com) 

I presume that spammers don't work for free, so someone must be paying them. Is Microsoft hiring less-than-reputable web marketing companies to do their bidding?

548. A lot has been happening with Kingsoft WPS

Seems like they have been making progress.

I wrote my original post about wps (also see this post) for version, and on this site they are currently at -- Alpha 12, patch 3. Judging from all the bugs that have been fixed it seems that the devs have been quite busy. The short interval between different alpha version also indicates that the development is very active. Good stuff.

Basically, WPS is the only reasonable (subjective) solution for reading .docx files on linux that I've encountered.

Either way,if you have a previous version of wps installed, remove it
sudo apt-get autoremove wps-office:i386

Then get on with it:
mkdir ~/tmp
cd ~/tmp
sudo dpkg --add-architecture i386
sudo apt-get update
wget http://wdl.cache.ijinshan.com/wps/download/Linux/unstable/kingsoft-office_9.1.0.4244~a12p3_i386.deb
sudo dpkg -i kingsoft-office_9.1.0.4244~a12p3_i386.deb
sudo apt-get install -f

Missing fonts. Clicking on the link takes you indirectly to bbs.wpn.cn, which is in Chinese.

To sort out the fonts issue above you can join bbs.wpn.cn as shown in this post, then download the wps_symbol_fonts.zip file from this post. Install the fonts by doing
cp ~/Downloads/wps_symbol_fonts.zip ~/.fonts
cd ~/.fonts
unzip wps_symbol_fonts.zip

And start WPS. The 'symbol' issue should be solved. Either way, the look of WPS has been updated, and now it can handle my test file which even MS offerings are struggling with:
See here for a post where I compare different office suites using that file (older post): http://verahill.blogspot.com.au/2013/01/313-which-office-for-linux-users.html

To sort out mime types/file associations do
sudo update-desktop-database /usr/local/share/applications/
sudo update-mime-database /usr/share/mime/

02 February 2014

547. NMR on Debian Jessie: gSim

Note: the gsim deb needs glibc 2.14, while wheezy has version 2.13. See here for how to deal with glibc > 2.14 on wheezy: http://verahill.blogspot.com.au/2014/03/562-pulling-in-glibc-214-from-testing.html

gSim is a free program (doesn't appear to be open source -- I can't find the source anyway) for processing 1D and 2D nmr spectra. I don't have much to say about it other than that it seems to provide all the things I require of an NMR processing program, including line fitting (deconvolution).

One neat extra is that gsim has a terminal-like window at the bottom similar to xwinnmr/topspin so that you can control the processing completely via the keyboard. 

gSim can be downloaded from http://sourceforge.net/projects/gsim/

installation via gdebi

splash screen

can read bruker files


overview of processing options. I had to select bruker to correct the fid before doing ft

zooming is a bit awkward, but can be done using the mouse

I had some issue with deconvolution

awkwardly it uses ppm for the width too if the x axis is in ppm

you can always switch to hz though
gsim also reads the acqu/s files

It can output the spectrum to eps in case you're too lazy to work in gnuplot