10 June 2014

580. Python script: Interpolate between structures in a multi-xyz file

I'm doing a lot of NEB (nudged elastic band) calculations using nwchem at the moment, and while getting 'neb convergence' is simple enough, I often get an error along the lines of
3965547 @neb NEB calculation converged
3965548 @neb However NEB Gmax not converged...Try increasing the number of beads.

While that sounds simple enough it's nicer if you don't have to go back to the beginning and e.g. run a more fine-grained PES job to generate a reasonable trajectory (straight linear interpolation often doesn't work), then keep on running neb iterations. One way to cut down on time (presumably) is to simply pad the neb path xyz with intermediate structures, and that is what this python (2.7) script does.

Oh, and I really wish blogspot would support code inclusion better...

How to use:
python nebinterpolate.py -i neb_A_F.neb_final.xyz  -o test.xyz

Example:
Say we have the structure of methanol, and methanol in which the oxygen-carbon distance is 3.0 Ångström:

Here's the corresponding xyz file, which we'll call first.xyz:
6 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 0.83960 0.81877 0.83960 H 1.38912 0.20156 1.38912 6 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.76087 1.75017 1.76087 H 2.31039 1.13296 2.31039

Run nebinterpolate -i first.xyz -o second.xyz and you'll get a new xyz file with three structures -- the first one plus an intermediate structure:

Here's the new file, second.xyz:
6 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 0.83960 0.81877 0.83960 H 1.38912 0.20156 1.38912 6 structure 2 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.30024 1.28447 1.30024 H 1.84976 0.66726 1.84976 6 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.76087 1.75017 1.76087 H 2.31039 1.13296 2.31039

 Run it again, nebinterpolate -i second.xyz -o third.xyz, and you'll get:

Here's the new file, third.xyz:
6 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 0.83960 0.81877 0.83960 H 1.38912 0.20156 1.38912 6 structure 2 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.06992 1.05162 1.06992 H 1.61944 0.43441 1.61944 6 structure 2 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.30024 1.28447 1.30024 H 1.84976 0.66726 1.84976 6 structure 4 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.53056 1.51732 1.53056 H 2.08007 0.90011 2.08007 6 C 0.03517 0.00549 0.03517 H -0.61778 -0.63407 0.66798 H 0.66798 -0.63407 -0.61778 H -0.60514 0.64647 -0.60514 O 1.76087 1.75017 1.76087 H 2.31039 1.13296 2.31039

Now, the real use for this is when you've been optimising a string of structures using NEB and want to increase the number of images because you're not getting gmax convergence, or if you want to do a quick and rough optimisation and then get a prettier looking set of coordinates.

You can load the multi-xyz file in nwchem by using

NEB
   ...
  XYZ_PATH path.xyz
END


nebinterpolate.py 

#!/usr/bin/python

import sys

def getvars(arguments):
 switches={}

 version='0.1'
 
 try:
  if "-i" in arguments:
   switches['in_one']=arguments[arguments.index('-i')+1]
   print 'Input: %s '% (switches['in_one'])
  else:
   arguments="--help";
 except:
  arguments="--help";
  
 try:
  if "-o" in arguments:
   switches['o']=arguments[arguments.index('-o')+1].lower()
   print 'Output: %s'% switches['o']
  else:
   arguments="--help";
 except:
  arguments="--help";

 try:
  if "-w" in arguments:
   switches['w']=float(arguments[arguments.index('-w')+1])
   print 'Weighting: %i'% switches['w']
  else:
   print 'Assuming no weighting'
   switches['w']=1.0;
 except:
  switches['w']=1.0;

 doexit=0
 try:
  if ("-h" in arguments) or ("--help" in arguments):
   print '\t\t bytes2words version %s' % version
   print 'Creates interpolated structures'
   print 'from an multixyz file'
   print '--input \t-i\t multi-xyz file to morph'
   print '--output\t-o\t output file'
   print '--weight\t-w\t weight first structure vs second one (1=average; 0=start; 2=end)'
   print 'Exiting'
   doexit=1
 except:
  a=0 # do nothing
 if doexit==1:
  sys.exit(0)

 return switches

def getcmpds(switches):
 
 cmpds={}
 
 g=open(switches['in_one'],'r') 
 n=0
 xyz=[]
 atoms=[]
 structure_id=1

 for line in g:

  try:
   if len(xyz)==cmpds['atoms_'+str(structure_id)]:
    cmpds['coords_'+str(structure_id)]=xyz
    cmpds['elements_'+str(structure_id)]=atoms   
    structure_id+=2
    n=0
    atoms=[]
    xyz=[]
  except:
   pass

  n+=1
  line=line.rstrip('\n')
   
  if n==1:
   cmpds['atoms_'+str(structure_id)]=int(line)
  elif n==2:
   cmpds['title_'+str(structure_id)]=line
  else:
   line=line.split(' ')
   line=filter(None,line)
   xyz+=[[float(line[1]),float(line[2]),float(line[3])]]
   atoms+=[line[0].capitalize()]
 g.close

 cmpds['coords_'+str(structure_id)]=xyz
 cmpds['elements_'+str(structure_id)]=atoms   
 cmpds['w']=switches['w']
 cmpds['structures']=(structure_id)
 
 return cmpds

def morph(cmpds,structure_id):
 coords_one=cmpds['coords_'+str(structure_id)]
 coords_two=cmpds['coords_'+str(structure_id+2)]
 
 coords_morph=[]
 coords_diff=[]
 for n in range(0,cmpds['atoms_'+str(structure_id)]):
  morph_x=coords_one[n][0]+cmpds['w']*(coords_two[n][0]-coords_one[n][0])/2.0
  morph_y=coords_one[n][1]+cmpds['w']*(coords_two[n][1]-coords_one[n][1])/2.0
  morph_z=coords_one[n][2]+cmpds['w']*(coords_two[n][2]-coords_one[n][2])/2.0
  diff_x=coords_two[n][0]-coords_one[n][0]
  diff_y=coords_two[n][1]-coords_one[n][1]
  diff_z=coords_two[n][2]-coords_one[n][2]
  coords_morph+=[[morph_x,morph_y,morph_z]]
  coords_diff+=[[diff_x,diff_y,diff_z]]

 cmpds['atoms_'+str(structure_id+1)]=cmpds['atoms_'+str(structure_id)]
 cmpds['elements_'+str(structure_id+1)]=cmpds['elements_'+str(structure_id)]
 cmpds['title_'+str(structure_id+1)]='structure '+str(structure_id+1)
 cmpds['coords_'+str(structure_id+1)]=coords_morph

 return cmpds

def genxyzstring(coords,element):
 x_str='%10.5f'% coords[0]
 y_str='%10.5f'% coords[1]
 z_str='%10.5f'% coords[2]
 
 xyz_string=element+(3-len(element))*' '+10*' '+\
 (8-len(x_str))*' '+x_str+10*' '+(8-len(y_str))*' '+y_str+10*' '+(8-len(z_str))*' '+z_str+'\n'
 
 return xyz_string

def writemorph(cmpds,outfile):
 g=open(outfile,'w') 

 for m in range(1,cmpds['structures']+1):
  g.write(str(cmpds['atoms_'+str(m)])+'\n')
  g.write(str(cmpds['title_'+str(m)])+'\n')
  for n in range(0,cmpds['atoms_'+str(m)]):
   coords=cmpds['coords_'+str(m)][n]
   g.write(genxyzstring(coords, cmpds['elements_'+str(m)][n]))
 g.close

 return 0

if __name__=="__main__":
 arguments=sys.argv[1:len(sys.argv)]
 switches=getvars(arguments)
 cmpds=getcmpds(switches)

#check that the structures are compatible
 for n in range(1,cmpds['structures'],2):

  if cmpds['atoms_'+str(n)]!=cmpds['atoms_'+str(n+2)]:
   print 'The number of atoms differ. Exiting'
   sys.exit(1)
  elif cmpds['elements_'+str(n)]!=cmpds['elements_'+str(n+2)]:
   print 'The types of atoms differ. Exiting'
   sys.exit(1)
  cmpds=morph(cmpds,n)
 success=writemorph(cmpds,switches['o'])
 if success==0:
  print 'Conversion seems successful'

29 May 2014

579. TEMPerHum (0c45:7402 ) on debian linux wheezy

Yet another little usb device from PC sensors (I seem to have been giving them a fair amount of money recently ).

This is the device in question: http://pcsensor.com/index.php?_a=product&product_id=178 (note that I buy my devices via ebay where the prices are apparently always the campaign ones -- I paid AUD25 , not USD80)

It didn't come in a box so I have no scans or shots to show.

[See here for the regular TEMPer USB device (0c45:7401), and here for the TEMPer 1K4  (0c45:7403) USB thermocouple reader. Note that since the temper-usb code is in  a lot of flux you can't use the line numbers in those posts directly -- you'll have to read and understand the code before pasting it in. Luckily, it's quite simple -- even I managed to sort it out!]

I couldn't get temper-usb to work even when making (what I consider) the necessary edits, but instead got lots of errors (including " usb.USBError: could not detach kernel driver from interface 1: No data available"). So I finally gave up.

Instead, web searching led me to http://edorfaus.wordpress.com/2012/07/02/new-library-examples-and-features/ -- one of the replies by eg1l spelled out the solution -- I'll grant myself the liberty to repost it here, but please remember where it originally came from and link to the original article (exclusively or in addition).

mkdir ~/tmp
cd ~/tmp/
sudo apt-get install libudev-dev libusb-1.0.0-dev libfox-1.6-dev autoconf cmake

git clone git://github.com/signal11/hidapi.git
cd hidapi/
./bootstrap
./configure
cd libusb/
make
sudo make install

cd ../../

git clone git://github.com/edorfaus/TEMPered.git
cd TEMPered/
mkdir build
cd build/
cmake ..
make
cd utils/
sudo ./tempered 
cd ../
sudo make install

sudo ln -s /usr/local/lib/x86_64-linux-gnu/libtempered.so.0 /usr/local/lib/libtempered.so.0
sudo ln -s /usr/local/lib/x86_64-linux-gnu/libtempered-util.so.0 /usr/local/lib/libtempered-util.so.0
sudo ldconfig

sudo tempered
0006:000e:01 0: temperature 21.75 °C, relative humidity 49.1%, dew point 10.6 °C

Create an 80-temper.rules file in /etc/udev/rules.d:
SUBSYSTEMS=="usb", ATTRS{idVendor}=="0c45", ATTRS{idProduct}=="7401", GROUP="users", MODE="0666"
SUBSYSTEMS=="usb", ATTRS{idVendor}=="0c45", ATTRS{idProduct}=="7402", GROUP="users", MODE="0666"
SUBSYSTEMS=="usb", ATTRS{idVendor}=="0c45", ATTRS{idProduct}=="7403", GROUP="users", MODE="0666"

Then do
sudo usermod -a -G users $USER

And
sudo service udev restart

Unplug and re-plug your device, then open a new terminal and you're set (type group to make sure that users show up).

me@boron:~$ tempered
0006:000f:01 0: temperature 23.34 °C, relative humidity 46.8%, dew point 11.3 °C

20 May 2014

578. More Gnome issues -- semi-rant

I was able to get on with work in spite of transitioning to gnome 3 -- but only thanks to the frippery extensions:
http://verahill.blogspot.com.au/2011/11/debian-testing-64-bit-gnome-3gnome.html 
http://verahill.blogspot.com.au/2013/04/402-very-briefly-what-i-forgot-about.html

And then gnome-screenshot got crippled, but I managed to cope with that by patching it:
http://verahill.blogspot.com.au/2012/06/fixing-gnome-screenshot-341-in-debian.html

I'm using debian testing (jessie) on my laptop and since I normally don't do much on it other than occasionally log into work to check on jobs I have been able to ignore the issues that are so apparent in gnome 3.12, two of which are:

*  gnome-terminal doesn't have a transparent background option since ersion 3.8 -- instead of being able to read what's underneath (e.g. a blog post with a how-to), and thus making good use of the screen real estate, my laptop screen is now feeling very, very small.

[And the developer seems to have the usual gnome attitude issue: https://bugzilla.gnome.org/show_bug.cgi?id=698544 ]

gnome-terminal 3.12
 There's basically nothing that has forced me to stay with gnome-terminal, so switching terminals was a simple matter of installing xfce4-terminal instead, which looks pretty much like gnome-terminal used to. Besides, I use guake more of the time anyway.

xfce4-terminal

 * you can't resize 'native' gnome programs, such as nautilus, that are in full-screen mode. Well, you can -- by holiding super (i.e. windows button) + down arrow. But the icons in the top right corner are no longer, and I find that incredibly annoying.
Nautilus. Nautilus shows when dropbox is synced -- thunar doesn't.

Thunar -- like nautilus of days gone by
 There are two things that I need from nautilus (other than working as a filemanager, of course, and allowing me to open terminals where I want which is something thunar supports 'out of the box') -- the ability to batch resize images, and dropbox. Now, dropbox has nothing to do with nautilus and the dropbox server will run happily in the background whether you're using a file manager or not. It would be nice if thunar would show whether the dropbox folder is synced or not, but it's hardly crucial.

Image resizing is a different matter since we maintain our own little website with personal photos for overseas members of our families. Ergo, it's crucial.

Luckily, this is pretty easy:
sudo apt-get install simple-image-reducer

Go to Edit/Configure custom actions... and set command to simple-image-reducer %F, and scope to apply to images, as shown in the screenshots.





I'm sure there are lots of other small issues in gnome 3.12 that I either keep missing or wilfully ignore in the interest of maintaining a low blood pressure. The whole idea of removing features that are actually useful with no way of re-enabling them smacks of 'we know better than you', and that irks me.