515. Very briefly: boot issues and nvram (anecdote)

I had to move offices at work last Friday (which partly explains my radio silence).

One of my nodes, with a Biostar 880G+ motherboard, wouldn't boot. The issue manifested itself in a few different ways (and not always at the same time):

Most of the time the bios portion of the boot would be very slow, and I'd get a loooooong beep when the voltages were echoed during what I presume is POST. I never made it to the GRUB screen. I'd had this happen once or twice before, especially when the computer had been turned off and/or unplugged for a longer period of time.

Occasionally I'd get a message about NVRAM being checked..

After
1) having checked that nothing had vibrated loose during the transport
2) having removed all cards and drives except a single RAM stick
3) having checked that RAM stick
I still wasn't getting any closer to getting booted. So I went into the BIOS. Having seen the message about nvram and having done a bit of reading online I selected to erase the NVRAM.

And suddenly booting worked...once.

After that I kept getting messages about nvram being corrupted and boot fail, or 'NVRAM...Update failed'.

I then set the CMOS jumpers to pin 2 and 3 (small jumper next to the CMOS battery)

And then it worked! It booted just fine! But I kept getting the following error echoing on my screen:

[  120.718683] hub 3-0:1.0: unable to enumerate USB device on port 2
[  120.932378] usb 3-2: new full-speed USB device number 28 using ohci_hcd
[  121.072497] usb 3-2: device descriptor read/64, error -62
[  121.316673] usb 3-2: device descriptor read/64, error -62
[  121.556904] usb 3-2: new full-speed USB device number 29 using ohci_hcd
[  121.697019] usb 3-2: device descriptor read/64, error -62
[  121.941224] usb 3-2: device descriptor read/64, error -62
[  122.181425] usb 3-2: new full-speed USB device number 30 using ohci_hcd

Unplugged all the USB devices. Turned out to be the KVM (ritmo). Turning it off (i.e. unplugging all USB input to the KVM so that it lost power) and plugging it in again solved it.

And ever since then all has been fine.

Presumably it's time to change the CMOS battery.

514. Extracting Frequency data from a gaussian 09 calculation for gnuplot

This is another python script.

Say you've done a computation along the lines of this:

#P rBP86/GEN 5D Pseudo(Read) Opt=() Freq=() SCF=(MaxCycle=999 ) Punch=(MO) Pop=()

and want the data in a neat data file, like this:

33.237 0.0023 0.0536
39.9976 0.0043 0.8305
69.7345 0.0129 0.3348
84.7005 0.0173 0.7027
[..]
3133.0068 6.2938 0.6114
3143.8021 6.3551 0.3775
3164.9242 6.4685 0.8829
3221.8787 6.6972 4.6005

Then you can use the following python (2.x) script, g09freq:

#!/usr/bin/python
# Compatible with python 2.7 
# Reads frequency output from a g09 (gaussian) calculation
# Usage ex.: g09freq g09.log ir.dat
import sys 

def ints2float(integerlist):
    for n in range(0,len(integerlist)):
        integerlist[n]=float(integerlist[n])
    return integerlist

def parse_in(infile):
    g09output=open(infile,'r')
    captured=[]
    for line in g09output:
        if ('Frequencies' in line) or ('Frc consts' in line) or ('IR Inten' in line):
            captured+=[line.strip('\n')]
    g09output.close()
    return captured
    
def format_captured(captured):
    vibmatrix=[]
    steps=len(captured)
    for n in range(0,steps,3):
        freqs=ints2float(filter(None,captured[n].split(' '))[2:5])
        forces=ints2float(filter(None,captured[n+1].split(' '))[3:6])
        intensities=ints2float(filter(None,captured[n+2].split(' '))[3:6])
        for m in range(0,3):
            vibmatrix+=[[freqs[m],forces[m],intensities[m]]]
    return vibmatrix

def write_matrix(vibmatrix,outfile):
    f=open(outfile,'w')
    for n in range(0,len(vibmatrix)):
        item=vibmatrix[n]
        f.write(str(item[0])+'\t'+str(item[1])+'\t'+str(item[2])+'\n')
    f.close()
    return 0

if __name__ == "__main__":
    infile=sys.argv[1]
    outfile=sys.argv[2]

    captured=parse_in(infile)

    if len(captured)%3==0:
        vibmatrix=format_captured(captured)
    else:
        print 'Number of elements not divisible by 3 (freq+force+intens=3)'
        exit()
    success=write_matrix(vibmatrix,outfile)
    if success==0:
        print 'Read %s, parsed it, and wrote %s'%(infile,outfile)

Run it as e.g.

g09freq g09.log test.out

The output is compatible with gnuplot:

gnuplot> set xrange [3500:0]
gnuplot> set yrange [10:-1]
gnuplot> plot './test.out' u 1:2 w impulse

It's trivial to add gaussian broadening (see e.g. this post)

513. Extracting data from a PES scan with gaussian

There are a few reasons to like gaussian, and many reasons not to. Gaussian is fast, and their whitepapers are great resources for learning computational techniques.

Without going into discussions about the commercial behaviour of Wavefunction inc., the things I don't like about gaussian is the clunky input format (nwchem has a much more readable syntax), the inscrutable error messages, and the unreadable output. Well, it's not unreadable in a literal sense, but it could certainly be clearer. On the other hand, I've having issues with running some of my PES scans in nwchem -- and I can't find a solution (more about that in a later post)

Anyway, here's a python script for extracting optimized structures and energies from a relaxed PES scan in Gaussian 09.

First, an example of a simple scan:

%nprocshared=2
%Chk=methanol.chk
#P rB3LYP/6-31g 6D 10F Opt=(modredundant) NoSymm  Punch=(MO) Pop=() 

methanol

0 1 ! charge and multiplicity
 C     0.0351714     0.00548884     0.0351714
 H     -0.617781     -0.634073     0.667983
 H     0.667983     -0.634073     -0.617781
 H     -0.605139     0.646470     -0.605139
 O     0.839603     0.818768     0.839603
 H     1.38912     0.201564     1.38912

 1 5 S 10 0.1

And here's the script, pes_parse_g09:

#!/usr/bin/python
import sys

def getrawdata(infile):
        f=open(infile,'r')
        opt=0
        geo=0
        struct=[]
        structure=[]
        energies=[]
        energy=[]
        for line in f:
                
                if opt==1 and geo==1 and not ("---" in line):
                        structure+=[line.rstrip()]
                
                if 'Coordinates (Angstroms)' in line:
                        if opt==0:
                                opt=1
                                structure=[]
                        
                if opt==1 and "--------------------------" in line:
                        if geo==0:
                                geo=1
                        elif geo==1:
                                geo=0
                                opt=0
                if 'SCF Done' in line:
                        energy=filter(None,line.rstrip('\n').split(' '))
                if      'Optimization completed' in line and (opt==0 and geo==0):
                        energies+=[float(energy[4])]
                        opt=0
                        geo=0
                        struct+=[structure]
                        structure=[]
        
        return struct, energies

def periodictable(elementnumber):
        ptable={1:'H',2:'He',\
        3:'Li', 4:'Be',5:'B',6:'C',7:'N',8:'O',9:'F',10:'Ne',\
        11:'Na',12:'Mg',13:'Al',14:'Si',15:'P',16:'S',17:'Cl',18:'Ar',\
        19:'K',20:'Ca',\
        21:'Sc',22:'Ti',23:'V',24:'Cr',25:'Mn',26:'Fe',27:'Co',28:'Ni',29:'Cu',30:'Zn',\
        31:'Ga',32:'Ge',33:'As',34:'Se',35:'Br',36:'Kr',\
        37:'Rb',38:'Sr',\
        39:'Y',40:'Zr',41:'Nb',42:'Mo',43:'Tc',44:'Ru',45:'Rh',46:'Pd',47:'Ag',48:'Cd',\
        49:'In',50:'Sn',51:'Sb',52:'Te',53:'I',54:'Xe',\
        55:'Cs',56:'Ba',\
        57:'La',58:'Ce',59:'Pr',60:'Nd',61:'Pm',62:'Sm',63:'Eu',64:'Gd',65:'Tb',66:'Dy',67:'Ho',68:'Er',69:'Tm',70:'Yb',71:'Lu',\
        72:'Hf', 73:'Ta', 74:'W',75:'Re', 76:'Os', 77:'Ir',78:'Pt', 79:'Au', 80:'Hg',\
        81:'Tl', 82:'Pb', 83:'Bi',84:'Po',85:'At',86:'Rn',\
        87:'Fr',88:'Ra',\
        89:'Ac',90:'Th',91:'Pa',92:'U',93:'Np',94:'Pu',95:'Am',96:'Cm',97:'Bk',98:'Cf',99:'Es',100:'Fm',101:'Md',102:'No',\
        103:'Lr',104:'Rf',105:'Db',106:'Sg',107:'Bh',108:'Hs',109:'Mt',110:'Ds',111:'Rg',112:'Cn',\
        113:'Uut',114:'Fl',115:'Uup',116:'Lv',117:'Uus',118:'Uuo'}
        element=ptable[elementnumber]
        return element

def genxyzstring(coords,elementnumber):
        x_str='%10.5f'% coords[0]
        y_str='%10.5f'% coords[1]
        z_str='%10.5f'% coords[2]
        element=periodictable(int(elementnumber))
        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 getstructures(rawdata):
        
        n=0
        for structure in rawdata:
                
                n=n+1
                num="%03d" % (n,)
                g=open('structure_'+num+'.xyz','w')
                itson=False
                cartesian=[]
                        
                for item in structure:
                        
                        coords=filter(None,item.split(' '))
                        coordinates=[float(coords[3]),float(coords[4]),float(coords[5])]
                        element=coords[1]
                        cartesian+=[genxyzstring(coordinates,element)]
                g.write(str(len(cartesian))+'\n')
                g.write('Structure '+str(n)+'\n')
                for line in cartesian:
                        g.write(line)
                g.close()
                cartesian=[]
        return 0
        
if __name__ == "__main__":
        infile=sys.argv[1]
        rawdata,energies=getrawdata(infile)
        structures=getstructures(rawdata)
        g=open('energies.dat','w')
        for n in range(0,len(energies)):
                g.write(str(n)+'\t'+str(energies[n])+'\n')
        g.close()

And here's what we get from the output:

g09 methanol.in |tee methanol.out
pes_parse_g09 methanol.log
cat structure* > meoh_traj.xyz

And here's a plot of energies.dat: