20 April 2014

573. PCSensors K-type USB thermocouple adapter TEMPer1K4 (0c45:7403) on debian

UPDATE 1 May 2014:
I've rewritten the code now to properly deal with the presence of both a 0c45:7401 and a 0c45:7403 simultaneously. Note that I'm not convince about the accuracy of the temperature readings (accuracy as in whether there's any systematic bias to the output. The reproducibility seems to be excellent though).

Original post:
I just received this usb thermocouple reader (TEMPer1k4): http://www.pcsensor.com/index.php?_a=product&product_id=104

About the product:
lsusb shows 0c45:7403 (Microdia Foot Switch)

dmesg shows
[13448.536120] usb 6-1: new low-speed USB device number 3 using uhci_hcd [13448.709126] usb 6-1: New USB device found, idVendor=0c45, idProduct=7403 [13448.709139] usb 6-1: New USB device strings: Mfr=1, Product=2, SerialNumber=0 [13448.709146] usb 6-1: Product: TH1000isoV1.5 [13448.709152] usb 6-1: Manufacturer: RDing

Getting it to work:
You can use the same program as in this post: http://verahill.blogspot.com.au/2013/12/532-temper-temperature-monitoring-usb.html with some minor modifications. Before running
sudo python setup.py install

make some changes in temperusb/temper.py.

Firstly, you need to allow for the detection of devices with the 0c45:7403 ID
15 VIDPIDs = [(0x0c45L,0x7401L),(0x0c45L,0x7403L)]
Secondly, you need to make sure that the program knows which device is which and treats them differently. Finally, you need to collect the right data-- the TEMPer1k4 returns a data_s with a length of 8 (caveat -- I haven't checked what the output from 0c45:7401 looks like. Maybe it too yields 8 fields), and the [2:4] data block contains the internal temperature of the USB device, while the [4:6] data block holds the thermocouple junction temperature. Of a sort -- you get a number which you need to translate into a temperature. See the bottom of this post for how I worked it all out. Please note that you should calibrate the thermocouple -- the uncorrected output seems to be at least 2-3 degrees too high.
 15 VIDPIDs = [(0x0c45L,0x7401L),(0x0c45L,0x7403L)]

 63     def getid(self):
 64         return self._device.idProduct

125             if self._device.idProduct==29697:
126                 temp_c = 125.0/32000.0*(struct.unpack('>h', data_s[2:4])[0])+0.0025
127                 if format == 'celsius':
128                     temp_c = temp_c * self._scale + self._offset
129                     return temp_c
130                 elif format == 'fahrenheit':
131                     return temp_c*1.8+32.0
132                 elif format == 'millicelsius':
133                     return int(temp_c*1000)
134                 else:
135                     raise ValueError("Unknown format")
136             elif self._device.idProduct==29699:
137                 temp_c = (struct.unpack('>h', data_s[4:6])[0])*0.25-2.00
138                 temp_internal=(125.0/32000.0*(struct.unpack('>h', data_s[2:4])[0]))+0.0025
139                 if format == 'celsius':
140                     temp_c = temp_c * self._scale + self._offset
141                     return temp_c
142                 elif format == 'fahrenheit':
143                     if self._device.idProduct==29697: #0x7401
144                         return temp_c*1.8+32.0
145                     elif self._device.idProduct==29699: #0x7403
146                         return temp_internal

 31     for i, dev in enumerate(devs):
 32         readings.append({'device': i,
 33                          'id':dev.getid(),
 34                          'temperature_c': dev.get_temperature(),
 35                          'temperature_f':
 36                          dev.get_temperature(format="fahrenheit"),
 37                          'ports': dev.get_ports(),
 38                          'bus': dev.get_bus()
 39                          })
 41     for reading in readings:
 42         if disp_ports:
 43             portinfo = " (bus %s - port %s)" % (reading['bus'],
 44                                                 reading['ports'])
 45         else:
 46             portinfo = ""
 47         if reading['id']==29697:
 48             print('Device #%i%s: %0.1f°C %0.1f°F'
 49                   % (reading['device'],
 50                      portinfo,
 51                      reading['temperature_c'],
 52                      reading['temperature_f']))
 53         elif reading['id']==29699:
 54             print('Device #%i%s: %0.1f°C %0.1f°C'
 55                   % (reading['device'],
 56                      portinfo,
 57                      reading['temperature_c'],
 58                      reading['temperature_f']))

Beyond this, follow the instructions in http://verahill.blogspot.com.au/2013/12/532-temper-temperature-monitoring-usb.html and make sure to set up a rules file with the correct USB ID for this device.


If you're curious about the details, have a look below.

Analysing USB traffic:
Since I wasn't quite sure where to start I figured the easiest approach would be to sniff the traffic between the usb device and the offically supported programme from PC Sensors. So I installed it Win XP in virtualbox, read a couple of temperatures and compared them with the captured data.

To capture data I did:
sudo apt-get install tshark
sudo modprobe usbmon
Bus 008 Device 003: ID 17ef:4815 Lenovo Integrated Webcam [R5U877] Bus 008 Device 001: ID 1d6b:0002 Linux Foundation 2.0 root hub Bus 006 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub Bus 005 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub Bus 004 Device 043: ID 0c45:7403 Microdia Foot Switch Bus 004 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub Bus 004 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub Bus 007 Device 001: ID 1d6b:0002 Linux Foundation 2.0 root hub Bus 003 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub Bus 002 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub Bus 001 Device 001: ID 1d6b:0001 Linux Foundation 1.1 root hub
sudo tshark -D
1. eth0 2. wlan0 3. nflog 4. nfqueue 5. usbmon1 6. usbmon2 7. usbmon3 8. usbmon4 9. usbmon5 10. usbmon6 11. usbmon7 12. usbmon8 13. any 14. lo (Loopback)
touch 1.pcap chmod ugo+w 1.pcap sudo tshark -i usbmon4 -w $HOME/1.pcap

I opened the 1.pcap file in wireshark and looked for fields called leftover data.

Here are two examples of such 'leftover data' fields:
Host to USB: 01 80 00 00 00 00 00 00
USB to Host: 80 06 17 f0 00 5c 0f ff

Subjecting the system to different temperatures helped me figure out what fields where changing -- I've marked them in bold above (blue for the thermocouple, red for the internal temperature). Interestingly, I got overflow when freezing the thermocouple between two ice cubes, which indicated that there was a fairly high lower limit (1.5 degrees)
The following are some of the data I harvested. The temperature is in the first column, the hex code is in the second one and the decimal value is in the third one:
23.63 17a 378
23.69 17b 379
23.94 17f 383
23.75 17c 380
21.88 15e 350
22.31 165 357
22.88 16e 366

Thermocouple: 23.75 05f 95 23.50 05e 9 38.25 099 153 24.00 066 102 87.50 15e 350 79.50 13e 318
In other words, T(thermo)=output(thermo)*0.25+1.50 and T(int.)=(output(int.)-0.0025)/0.0625. Because the internal data is in fields 2:4 it's seen as e.g. 17b0 instead of 17b, so the decimal version needs to be divided by 16 (line 126 above).

I also set up an /etc/temper.conf file and did chown $USER /etc/temper.conf in order to be able to read it (must be a better way). I tried to calibrate the thermocouple in my kitchen with iced water and boiling water. The boiling water gave a reading of 102 degrees, while the iced water gave me 7.5 degrees Celsius even after what should've been long enough to achieve equilibrium. I'll calibrate it in the lab later. So for now a simple offset might be enough to give a working temperature.
temper-poll -p
Device #0 (bus 2 - port 2)
so temper.conf became
2-2: scale = 1.00, offset = -4.0

The temper-usb code is changing too quickly!. You'll thus need to read and understand the code rather than just pasting it in. Here are the full, edited cli.py an temper.py files.

# encoding: utf-8
from __future__ import print_function
from temper import TemperHandler

def main():
    th = TemperHandler()
    devs = th.get_devices()
    readings = []
    print("Found %i devices" % len(devs))

    for i, dev in enumerate(devs):
        readings.append({'device': i,
                         'temperature_c': dev.get_temperature(),

    for reading in readings:
  if reading['id']==29697:
   print('Device #%i: %0.1f°C %0.1f°F' % (reading['device'],
  elif reading['id']==29699:
   print('Device #%i: %0.1f°C %0.1f°C' % (reading['device'],

# encoding: utf-8
# Handles devices reporting themselves as USB VID/PID 0C45:7401 (mine also says RDing TEMPerV1.2).
# Copyright 2012, 2013 Philipp Adelt 
# This code is licensed under the GNU public license (GPL). See LICENSE.md for details.

import usb
import sys
import struct

VIDPIDs = [(0x0c45L,0x7401L),(0x0c45L,0x7402L),(0x0c45L,0x7403L)]
TIMEOUT = 2000

class TemperDevice():
    def __init__(self, device):
        self._device = device
        self._handle = None

    def get_temperature(self, format='celsius'):
            if not self._handle:
                self._handle = self._device.open()
                except usb.USBError:
                except usb.USBError:
                self._handle.controlMsg(requestType=0x21, request=0x09, value=0x0201, index=0x00, buffer="\x01\x01", timeout=TIMEOUT) # ini_control_transfer

            self._control_transfer(self._handle, "\x01\x80\x33\x01\x00\x00\x00\x00") # uTemperatura
            self._control_transfer(self._handle, "\x01\x82\x77\x01\x00\x00\x00\x00") # uIni1
            self._control_transfer(self._handle, "\x01\x86\xff\x01\x00\x00\x00\x00") # uIni2
            self._control_transfer(self._handle, "\x01\x80\x33\x01\x00\x00\x00\x00") # uTemperatura
            data = self._interrupt_read(self._handle)
            data_s = "".join([chr(byte) for byte in data])

            if self._device.idProduct==29697: 
    temp_c = 125.0/32000.0*(struct.unpack('>h', data_s[2:4])[0])+0.0025
    if format == 'celsius':
     return temp_c
    elif format == 'fahrenheit':
     return temp_c*1.8+32.0
    elif format == 'millicelsius':
     return int(temp_c*1000)
     raise ValueError("Unknown format")
            elif self._device.idProduct==29699: 
    temp_c = (struct.unpack('>h', data_s[4:6])[0])*0.25-4.00
    temp_internal=(125.0/32000.0*(struct.unpack('>h', data_s[2:4])[0]))+0.0025
    if format == 'celsius':
     return temp_c
    elif format == 'fahrenheit':
     if self._device.idProduct==29697: #0x7401
      return temp_c*1.8+32.0
     elif self._device.idProduct==29699: #0x7403
      return temp_internal
    elif format == 'millicelsius':
     return int(temp_c*1000)
     raise ValueError("Unknown format")
        except usb.USBError, e:
            if "not permitted" in str(e):
                raise Exception("Permission problem accessing USB. Maybe I need to run as root?")
    def getid(self):
        #print self._device.idProduct
        return self._device.idProduct

    def close(self):
        if self._handle:
            except ValueError:
            self._handle = None

    def _control_transfer(self, handle, data):
        handle.controlMsg(requestType=0x21, request=0x09, value=0x0200, index=0x01, buffer=data, timeout=TIMEOUT)

    def _interrupt_read(self, handle):
        return handle.interruptRead(0x82, REQ_INT_LEN)

class TemperHandler():
    def __init__(self):
        busses = usb.busses()
        self._devices = []
        for bus in busses:
            self._devices.extend([TemperDevice(x) for x in bus.devices if (x.idVendor,x.idProduct) in VIDPIDs])

    def get_devices(self):
        return self._devices

16 April 2014

572. autorotate/superimpose python script

If you want to calculate reaction coordinates between two structures you need to make sure that the structures haven't been rotated or translated, something which easily happens if you allow symmetry in gaussian and (it seems) z-matrix in nwchem.

I've written a script that lets you take two structures and align and superimpose them so that only the atoms that take part in the reaction move.

It works by you defining a minimum of four atoms that aren't supposed to move /relative to each other/ (i.e. they can be translated/rotated --- just not relative to each other) between the two structures. Four atoms far from each other are ideal. You need to make sure that they also don't lie in the same plane, but form a three-dimensional space.

I've tried this on real molecules too and it works better than I'd ever dared to hope for. The more 'stationary' atoms that you can use to make up the transformation matrix, the better.

In this example atom F is in a different location in structures a and b. The structures have also been rotated relative to each other.

6 A A 0.00000 0.00000 1.00000 B 0.00000 1.00000 0.00000 C 1.00000 0.00000 0.00000 D -1.00000 0.00000 0.00000 E 0.00000 0.00000 -1.00000 F 0.00000 0.500000 0.00000
6 B A 1 0 0 B 0 1 0 C 0 0 -1 D 0 0 1 E -1 0 0 F 0 -1 0

./autorotate.py a.xyz b.xyz '1,2,3,4'
Selected atoms in molecules 1 and 2 ['A', 0.0, 0.0, 1.0] ['A', 1.0, 0.0, 0.0] ['B', 0.0, 1.0, 0.0] ['B', 0.0, 1.0, 0.0] ['C', 1.0, 0.0, 0.0] ['C', 0.0, 0.0, -1.0] ['D', -1.0, 0.0, 0.0] ['D', 0.0, 0.0, 1.0] Transformation max error: 3.33066907388e-16 Writing to a.rot.xyz

6 A A 1.00000 0.00000 0.00000 B -0.00000 1.00000 -0.00000 C -0.00000 0.00000 -1.00000 D -0.00000 0.00000 1.00000 E -1.00000 -0.00000 -0.00000 F -0.00000 0.50000 -0.00000
This is how it looks (note that the axis aren't aligned with (1,0,0; 0,1,0; 0,0,1) but seem to go through the centre of the molecule):


import sys
import numpy as np

#autorotate input_1.xyz input_2.xyz '1,2,3,4'
# need to pick at least four atoms that are not in the same plane
# input_1.xyz will be rotated to align with input_2.xyz
# you pick at least four atoms that should have the same positions
# relative to one another (i.e. distance and relative geometry). These 
# are then used to calculate an affine transform matrix which is used 
# to rotate and translate structure input_1.xyz to overlap with 
# structure 2

def formatinput(argument):

 for n in atoms:
 return infile,coord_sys
def getrawdata(infile):
 for line in f:
  if n<2: data-blogger-escaped-if="" data-blogger-escaped-line.rstrip="" data-blogger-escaped-n="" data-blogger-escaped-preamble="">1:
 return xyz, preamble

def getcoords(rawdata,preamble,atoms):
 for structure in rawdata:
  num="%03d" % (n,)
  for item in structure:
   coordx=filter(None,item.split(' '))
   if len(coordx)>len(coordy):
 return cartesian

def getstructures(rawdata,preamble):
 for structure in rawdata:
  num="%03d" % (n,)
  for item in structure:
   coordx=filter(None,item.split(' '))
   if len(coordx)>len(coordy):
 return cartesian

def affine_transform(atoms,structures):
# from http://stackoverflow.com/questions/20546182/how-to-perform-coordinates-affine-transformation-using-python-part-2
 for n in atoms:

 for n in atoms:

 primary = np.array(primaryatomcoords)
 secondary = np.array(secondaryatomcoords)
 primaryfull = np.array(structures[0])

 # Pad the data with ones, so that our transformation can do translations too
 n = primary.shape[0]
 pad = lambda x: np.hstack([x, np.ones((x.shape[0], 1))])
 unpad = lambda x: x[:,:-1]
 X = pad(primary)
 Y = pad(secondary)
 Xp= pad(primaryfull)

 # Solve the least squares problem X * A = Y
 # to find our transformation matrix A
 A, res, rank, s = np.linalg.lstsq(X, Y)

 transform = lambda x: unpad(np.dot(pad(x), A))

# print "Max error should be as small as possible if the rotation is successful"
# print "If max error is large you may have selected a bad set of atoms"
 print "Transformation max error:", np.abs(secondary - transform(primary)).max()
 return secondaryfull

def transform_xyz(tmatrix,newxyz):
 for n in newxyz:
  final_xyz+=[[ newcoord[0][0],newcoord[1][0],newcoord[2][0]]]
 return final_xyz

def genxyzstring(coords,elementnumber):
 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 write_aligned(aligned_structure,atom_coords,preamble,outfile):
 print "Writing to ",outfile
 for n in range(0,len(aligned_structure)):
 return 0
if __name__ == "__main__":
 #get raw data

 #collect structures from raw data
 print "Selected atoms in molecules 1 and 2"
 for n in atoms:
  print atom_coords[0][n],atom_coords[1][n]
 #transform structure

09 April 2014

571. Briefly: Dodgy/underpowered UPS?

I've built quite a few computers in the past, and in general I haven't had any issues beyond the odd dodgy RAM stick.

However, a while back I became careless and built a box ('Oxygen') where the motherboard didn't officially support the CPU. Swapping CPUs with another box seemed to solve the issues I had.

See e.g.

In the past couple of weeks I've begun to see some worrying signs that all isn't right. In particular I noticed the following in the dmesg output:
[693166.514897] [Hardware Error]: MC2 Error: VB Data ECC or parity error. [693166.514926] [Hardware Error]: Error Status: Corrected error, no action required. [693166.514934] [Hardware Error]: CPU:6 (15:1:2) MC2_STATUS[-|CE|MiscV|-|-|-|-|CECC]: 0x98414000010c0176 [693166.514955] [Hardware Error]: cache level: L2, tx: DATA, mem-tx: EV

A few days after that, the computer turned itself off without returning any additional error messages. It did cause me to look at the sensor output though (I've been logging it every two minutes for months), and I compared it with another computer ('Neon') which is completely stable. Note that both computers have been running the same types of jobs recently (large memory frequency jobs).

Hardware specs:
Oxygen: AMD FX8150, 32 gb ram, Corsair GS700, asrock 990 fx extreme3
Neon: AMD FX8350, 32 gb ram, Corsair GS800, gigabyte 990 fxa

Anyway, this is what I found:
On Neon the power output is very stable, while on Oxygen it jumps up and down between ca 45 W and 130 W.

Has it been a crappy UPS that has been causing the issues all along? Or do these plot mean nothing?

570. Briefly: restarting a g09 frequency job with SGE, using same queue

I've had g09 frequency jobs die on me, and in g09 analytical frequency jobs can only be restarted using the .rwf. Because the .rwf files are 160 gb, I don't want to be copying them back and forth between nodes. It's easier then to simply make sure that the restarted job is run on the same node as the original job.

A good resource for SGE related stuff is http://rous.mit.edu/index.php/SGE_Instructions_and_Tips#Submitting_jobs_to_specific_queues

Either way, first figure out what node the job ran on. Assuming that the job number was 445:
qacct -j 445|grep hostname
hostname compute-0-6.local

Next figure out the PID, as this is used to name the Gau-[PID].rwf file:
grep PID g03.g03out
Entering Link 1 = /share/apps/gaussian/g09/l1.exe PID= 24286.

You can now craft your restart file, g09_freq.restart -- you'll need to make sure that the paths are appropriate for your system:
%nprocshared=8 %Mem=900000000 %rwf=/scratch/Gau-24286.rwf %Chk=/home/me/jobs/testing/delta_631gplusstar-freq/delta_631gplusstar-freq.chk #P restart
(having empty lines at the end of the file is important) and a qsub file, g09_freq.qsub:
#$ -S /bin/sh #$ -cwd #$ -l h_rt=999:30:00 #$ -l h_vmem=8G #$ -j y #$ -pe orte 8 export GAUSS_SCRDIR=/tmp export GAUSS_EXEDIR=/share/apps/gaussian/g09/bsd:/share/apps/gaussian/g09/local:/share/apps/gaussian/g09/extras:/share/apps/gaussian/g09 /share/apps/gaussian/g09/g09 g09_freq.restart > g09_freq.out
Then submit it to the correct queue by doing
qsub -q all.q@compute-0-6.local g09_freq.qsub

The output goes to g09_freq.log. You know if the restart worked properly if it says
Skip MakeAB in pass 1 during restart.
Resume CPHF with iteration 214.

Note that restarting analytical frequency jobs in g09 can be a hit and miss affair. Jobs that run out of time are easy to restart, and some jobs that die silently have also been restarted successfully. On the other hand, a job that died because my resource allocations ran out couldn't be restarted i.e. restart started the freq job from scratch. The same happened with one a node of mine that has what seems like a dodgy PSU. Finally, I also couldn't restart jobs that died silently due to allocation all the RAM to g09 without leaving any to the OS (or at least that's the current best theory). It may thus be a good idea to back up the rwf file every now and again, in spite of the unwieldy size.

03 April 2014

569. Briefly: Dual monitor set-up on Debian Wheezy with Gnome and a single nvidia graphics card

This is very easy, but I might as well document it here anyway.

This morning another group threw out two functioning monitors and I grabbed both. While I haven't yet decided on what to do with the second one I decided to use one to make a dual monitor set-up for my work station.

My desktop has both onboard nvidia graphics and a separate pci-e nvidia (GT 430) graphics card. Using lspci only the external graphics card shows up, probably because the bios prioritises the external card and disables the onboard graphics.

The nvidia GT 430 card has three output ports: vga, hdmi and dvi. My main monitor (Dell P2411H, 1920x1080) has both vga and dvi, and my 'new' monitor (HP S1932, 1366x768) only has vga.

The first step was to physically connect both monitors to my computer. I originally thought I had to use one card per monitor, which would've necessitated me to reboot, change the bios and probably hand-craft an xorg.conf. I don't like rebooting, so I looked at the alternatives.

Apparently you can simply connect both monitors to the same card by using the different ports, so I hooked up the small screen to the vga port and the big one to the hdmi port.

After that it was a simple matter of opening 'displays' in the gnome 3 systems settings, setting both monitors to 'on', and arranging them side by side correctly by dragging them with the mouse:

I also had a look at it in nvidia-settings:

The only issue that remained was guake -- it was showing up in the 'wrong' screen (i.e. the left-most, smaller one).  This post showed how to edit: http://haifzhan.blogspot.com.au/2013/10/guake-dual-monitor-setup.html

This is how to do it on the version currently in wheezy, 0.4.3:
sudo vim `which guake`
814 def get_final_window_rect(self): 815 """Gets the final size of the main window of guake. The height 816 is the window_height property, width is window_width and the 817 horizontal alignment is given by window_alignment. 818 """ 819 screen = self.window.get_screen() 820 height = self.client.get_int(KEY('/general/window_height')) 821 width = 100 822 halignment = self.client.get_int(KEY('/general/window_halignment')) 823 824 # get the rectangle just from the first/default monitor in the 825 # future we might create a field to select which monitor you 826 # wanna use 827 828 #monitor = 0 # use the left most monitor 829 monitor = screen.get_n_monitors() - 1 # use the rightmost monitor 830 831 monitor_rect = screen.get_monitor_geometry(monitor) 832 window_rect = monitor_rect.copy() 833 total_width = window_rect.width 834 window_rect.height = window_rect.height * height / 100 835 window_rect.width = window_rect.width * width / 100 836 837 if width < monitor_rect.width: 838 if halignment == ALIGN_CENTER: 839 window_rect.x = monitor_rect.x+(monitor_rect.width - window_rect.width) / 2 840 elif halignment == ALIGN_LEFT: 841 window_rect.x = monitor_rect.x 842 elif halignment == ALIGN_RIGHT: 843 window_rect.x = monitor_rect.x+monitor_rect.width-window_rect.width 844 window_rect.y = monitor_rect.y 845 return window_rect 846

Note that the edited version will be overwritten when you upgrade guake.