10
$\begingroup$

Is there a way to easily extract all the structures and energies from an ADMP calculation done in Gaussian 09? GaussView can plot them but not export them in a useful way and babel just gives me the last structure in there.

I tried 

babel <G09output> all -oxyz <output.xyz>

but I get just 

==============================
*** Open Babel Error  in OpenAndSetFormat
  Cannot open all
1 molecule converted
1 errors 16 audit log messages

Or is there a different easy way to do this? Otherwise I guess I have to write my own script or look into cclib, but that one is using babel too if I understand it correct?

Example output file: http://dl.dropboxusercontent.com/u/6111391/Tz_TCO_dyn6.txt

Input:

%nprocshared=16
%mem=56GB
%chk=Tz_TCO_dyn6.chk
# 6-31g(d) scrf=check geom=(connectivity,crowd)admp=(stepsize=2500,maxpoints=400) rm062x IOp(1/44=990210)

[No Title]

0 1
 C                  1.08049600   -0.58720100   -2.07918800
 C                 -1.08049600    0.58720100   -2.07918800
 C                 -0.58945000    0.34949100    0.13978700
 H                 -0.52250100    1.43457800    0.23933400
 C                  0.58945000   -0.34949100    0.13978700
 H                  0.52250100   -1.43457800    0.23933400
 N                 -1.13318600   -0.73433600   -2.39136300
 N                  0.00078600   -1.35005000   -2.39250200
 N                 -0.00078600    1.35005000   -2.39250200
 N                  1.13318600    0.73433600   -2.39136300
 C                  1.84657900    0.28440600    0.66016000
 H                  2.75382900   -0.14516700    0.22100400
 H                  1.84124500    1.35503200    0.42328200
 C                 -1.84657900   -0.28440600    0.66016000
 H                 -2.75382900    0.14516700    0.22100400
 H                 -1.84124500   -1.35503200    0.42328200
 C                  1.87373600    0.08444000    2.19307200
 H                  2.72499800    0.64701600    2.59375900
 H                  2.07520100   -0.97273600    2.41061700
 C                 -1.87373600   -0.08444000    2.19307200
 H                 -2.72499800   -0.64701600    2.59375900
 H                 -2.07520100    0.97273600    2.41061700
 C                 -0.58945000   -0.50859700    2.93023900
 H                 -0.86607100   -0.70606500    3.97132400
 H                 -0.24860000   -1.47389400    2.53270700
 C                  0.58945000    0.50859700    2.93023900
 H                  0.86607100    0.70606500    3.97132400
 H                  0.24860000    1.47389400    2.53270700
 H                  2.02862900   -1.10167200   -1.96939700
 H                 -2.02862900    1.10167200   -1.96939700

 1 8 1.5 10 1.5 29 1.0
 2 7 1.5 9 1.5 30 1.0
 3 4 1.0 5 2.0 14 1.0
 4
 5 6 1.0 11 1.0
 6
 7 8 1.5
 8
 9 10 1.5
 10
 11 12 1.0 13 1.0 17 1.0
 12
 13
 14 15 1.0 16 1.0 20 1.0
 15
 16
 17 18 1.0 19 1.0 26 1.0
 18
 19
 20 21 1.0 22 1.0 23 1.0
 21
 22
 23 24 1.0 25 1.0 26 1.0
 24
 25
 26 27 1.0 28 1.0
 27
 28
 29
 30
Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Are you constrained to a specific platform? $\endgroup$
    – Martin - マーチン ♦
    Nov 12, 2016 at 14:36
  • $\begingroup$ Do you mean in terms of doing the calculations or for data analysis? Calculations are run on a supercomputer which runs on centos7 and I can use g09 and nwchem there, but for these calculations I prefer g09. For data analysis: not really, using Windows and Ubuntu here, and I'm fine with using whatever software or programming language is useful here. $\endgroup$
    – DSVA
    Nov 12, 2016 at 14:42
  • $\begingroup$ Yes I meant for analysis. I was assuming from your input stats that you are using a supercomputer. I am currently not in the office so I don't have access to all my toys, you might want to try a recent version of ChemCraft. Iirc you can save an optimisation as a xyz file with multiple structures, so it might work for this problem, too. (It is payware though, but probably the best there is, and if you don't already have it, you might want to consider getting it anyway.) $\endgroup$
    – Martin - マーチン ♦
    Nov 12, 2016 at 14:50
  • $\begingroup$ Oh thanks, that works at least for Scan and IRC and it's possible to put all the xyz coordinates in one file. But it doesn't work for dynamics. Chemcraft can read out energies and geometries (but it also says that there's an error) and whatever export option I use it's always the same single geometry. As far as I understand it it interprets this data as optimization steps and I haven't found an option to export those. But at least it shows where it found the values, so in worst case I could write a python script which exports everything. $\endgroup$
    – DSVA
    Nov 12, 2016 at 15:13
  • $\begingroup$ I recently found this again, and noticed there is a utility newzmat that has some capabilities. It's not ideal, but it should at least be available on systems that have G09/16 installed. $\endgroup$
    – Martin - マーチン ♦
    Jun 8, 2018 at 14:35

2 Answers 2

Reset to default
9
$\begingroup$

If you have the latest copy of cclib installed, it can extract both the energies and geometries and write the geometries to an XYZ trajectory file that you can open with Avogadro or VMD.

To print the HF or DFT energies to the screen,

$ ccget scfenergies admp.out

but be warned that the energies are in eV, not hartree.

To dump geometries to XYZ files, here's a simplified version of a script I use all the time:

#!/usr/bin/env python

"""cclib_extract_geom.py: Extract geometries from a quantum chemical
output file using cclib.
"""


from __future__ import print_function

import os.path

from cclib.io import ccopen, ccwrite


def getargs():
    """Get command-line arguments."""

    import argparse

    parser = argparse.ArgumentParser()

    parser.add_argument('outputfilename', nargs='+')

    parser.add_argument('--trajectory',
                        action='store_true',
                        help="""Extract all possible geometries into a trajectory-like file?""")
    parser.add_argument('--suffix')

    args = parser.parse_args()

    return args


if __name__ == '__main__':

    args = getargs()

    for outputfilename in args.outputfilename:

        job = ccopen(outputfilename)
        data = job.parse()

        stub = os.path.splitext(outputfilename)[0]
        if args.suffix:
            xyzfilename = ''.join([stub, '.', args.suffix, '.xyz'])
        else:
            xyzfilename = ''.join([stub, '.xyz'])

        ccwrite(data, outputdest=xyzfilename, allgeom=args.trajectory)

To write the energies in to a text file in hartree, you might add something like

        from cclib.parser.utils import convertor # put me at the top of the file
        with open(stub + '.energies', 'w') as fh:
            for energy in data.scfenergies:
                fh.write(str(convertor(energy, 'eV', 'hartree')) + '\n')
Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ If you're on Ubuntu, there's a cclib package, but these things are usually horribly out of date, so I'd download the latest copy from GitHub or pip install cclib. Come to chat if you need help. $\endgroup$
    – pentavalentcarbon
    Nov 12, 2016 at 21:16
  • $\begingroup$ Oh thanks! I figured that I will write some code myself, which I don't really mind because I can optimize it for whatever I need, but having a working script where I can actually see how to use it is great. I'm also quite new to python so... And now that you told me I looked and I had only cclib 1.3 installed but got the 1.5 version now. =) $\endgroup$
    – DSVA
    Nov 12, 2016 at 21:57
5
$\begingroup$

If, like me, you just want to extract the optimized structures from a potential energy scan, you can use this messy script:

#!/usr/bin/env python
# python Extract_Optimized_From_Gaussian.py filename

from __future__ import print_function
import sys, os

def extract_all(text, target):
    linenums = []
    # Start count at 1 because files start at line 1 not 0
    count = 1
    for line in text:
        if (line.find(target)) > -1:
            linenums.append(count)
        count += 1
    return linenums

code = {"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" : "Uub","113" : "Uut","114" : "Uuq","115" : "Uup", \
"116": "Uuh","117" : "Uus","118" : "Uuo"}

logfile_fn = os.path.basename(sys.argv[1])
logfile_bn = os.path.splitext(logfile_fn)[0]
logfile_fh = open(logfile_fn, 'r')
text = logfile_fh.readlines()
logfile_fh.close()

# Find all lines that contain "!   Optimized Parameters   !"
opt_param_line_nums = extract_all(text, '!   Optimized Parameters   !')
# Then find all lines that have "Input orientation:"
input_orient_line_nums = extract_all(text, 'Standard orientation')
# And all lines that have the ---- which we will use to find the ends of the coordinate sections
end_coor_line_nums = extract_all(text, '---------')
# Add 5 because the first actual coordinate is +5 lines
input_orient_line_nums = [x+5 for x in input_orient_line_nums]
# Make sure they're sorted
opt_param_line_nums.sort()
input_orient_line_nums.sort()
end_coor_line_nums.sort()

new_input_lns = []
count = 0
for x in reversed(input_orient_line_nums):
    if count > len(opt_param_line_nums)-1:
        break
    if x < list(reversed(opt_param_line_nums))[count]:
        new_input_lns.append(x)
        count += 1

new_input_lns.reverse()

new_end_lns = []
count = 0
for x in end_coor_line_nums:
    if count > len(new_input_lns)-1:
        break
    if x > new_input_lns[count]:
        new_end_lns.append(x)
        count += 1

intervals = zip(new_input_lns, new_end_lns)

# Now we convert each interval into a .xyz
outfile = open(logfile_bn+'.xyz','w')
for intvl in intervals:
    n_atoms = intvl[1] - intvl[0]
    outfile.write(str(n_atoms)+"\n\n")
    for x in range(intvl[0]-1, intvl[1]-1):
        column = text[x].split()
        print(code[column[1]],float(column[3]),float(column[4]),float(column[5]), file=outfile)

outfile.close()

I adapted this from here. I've tested it with Gaussian 09. No idea if it works with other versions. The resulting file loads into VMD. I don't know about other visualization software.

Edit: I checked the original script I had posted here against GaussView, and my script was wrong. Specifically, it was the structures titled "Standard orientation" that came before each "! Optimized Parameters !" line that I needed to extract, not the ones that came after. I have updated the script above to do that, and the output now matches that of GaussView. I removed the original script, but I have copied it to Pastbin here in case anyone still wants access to it (they're not too different).

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks, I will look into that and see if I get what's going on here:-D Always useful to see how people do this and it seems to work fine $\endgroup$
    – DSVA
    Dec 8, 2016 at 22:31
  • $\begingroup$ FYI, I realized a bit later today that I had an error in it. See my edit at the bottom. The script here is now correct. $\endgroup$
    – Nate
    Dec 9, 2016 at 4:07

Not the answer you're looking for? Browse other questions tagged or ask your own question.