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:


Previous
Next Post »