506. Extracting optimized structures from a potential energy scan in nwchem

Another update:
It now dumps the energies in a file, energies.dat, as well.

Update:
some programmes, like ecce, are more picky about the xyz format than others (e.g. jmol, vmd). I've updated the code to output xyz files that ecce too can read.

Original post:
When you use scan_input() in nwchem to do a PES scan (see e.g. here: http://verahill.blogspot.com.au/2013/08/503-relaxed-pes-scanning-in-nwchem.html) you get the energies and the gradients for the optimized structures returned as the results. However, for a casual user the atomic actual coordinates is more informative.

Here's a very simple parser written in python (2.7) which extracts the optimized structures from the output file:

#!/usr/bin/python
import sys

def getrawdata(infile):
        f=open(infile,'r')
        opt=0
        geo=0
        energy=[]
        energies=[]
        struct=[]
        structure=[]
        for line in f:
                if "Total DFT" in line:
                        line=filter(None,line.rstrip('\n').split(' '))
                        energy=float(line[4])
                if 'Optimization converged' in line:
                        opt=1
                if opt==1 and 'Geometry' in line:
                        geo=1
                if      'Atomic Mass' in line and (opt==1 and geo==1):
                        opt=0
                        geo=0
                        struct+=[structure]
                        energies+=[energy]
                        structure=[]
                if opt==1 and geo==1:
                        structure+=[line.rstrip()]
        return struct,energies

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 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:
                        
                        if itson and not(item==""):
                                coords=filter(None,item.split(' '))
                                coordinates=[float(coords[3]),float(coords[4]),float(coords[5])]
                                element=coords[1]
                                cartesian+=[genxyzstring(coordinates,element)]
                                #cartesian+=[coords[1]+'\t'+coords[3]+'\t'+coords[4]+'\t'+coords[5]+'\n']
                
                        if "---" in item:
                                itson=True
                        if item=="" and itson==True:
                                itson=False
                                if not(len(cartesian)==0):
                                        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()

Presuming that you've saved it as pes_parse.py you can then generate a series of xyz files with the structures, catenate them into a trajectory file, and open it in e.g. jmol. I'm using the output from example 1 in http://verahill.blogspot.com.au/2013/08/503-relaxed-pes-scanning-in-nwchem.html as the example:

chmod +x pes_parse.py
./pes_parse.py nwch.nwout
ls

nwch.nwout    structure_001.xyz  structure_003.xyz  structure_005.xyz  structure_007.xyz  structure_009.xyz  structure_011.xyz  structure_013.xyz  structure_015.xyz  structure_017.xyz  structure_019.xyz
pes_parse.py  structure_002.xyz  structure_004.xyz  structure_006.xyz  structure_008.xyz  structure_010.xyz  structure_012.xyz  structure_014.xyz  structure_016.xyz  structure_018.xyz

cat structure_*.xyz >> trajectory.xyz 
jmol trajectory.xyz

You can go through the structures by clicking on the arrows indicated by the white arrow:

Finally, using VMD it's easy to make videos -- note that they for some reason look awful here (seems like a lot of frames are removed, in particular from the beginning of the run):

And here's the SN2 reaction from post 503: