Extract atom coordinates from Gaussian .log file and compute covergence

log2xyz.py

# -*- coding: utf-8 -*-
"""
Created on Thu Feb 28 19:56:58 2019

@author: Attila Lengyel
attila@lengyel.nl
"""

import numpy as np
import os
import sys
import matplotlib.pyplot as plt

def log2xyz(inputpath,comment):
    print('Working on it...')
    
    filename = os.path.splitext(os.path.basename(inputpath))[0]
    
    # Read logfile
    with open(inputpath, 'r') as logfile:
        log = logfile.readlines()
    
    # Get linenumbers where xyz tables start
    indices = [i for i, s in enumerate(log) if 'Standard orientation' in s]
    indices = np.asarray(indices) + 5 # offset header
    
    distances = []
    
    # Process all tables
    for i, linenumber in enumerate(indices):
        lines = []
        # Store lines in list until end of table is found
        for line in log[linenumber:]:
            if '---------------------------------------------------------------------' in line:
                break
            else:
                line = line.split()
                lines.append(line)
        
        # Convert to numpy array
        data = np.asarray(lines).astype(float)
        data = data[:,[1,3,4,5]]
        
        # Output xyz file
        np.savetxt('{}_{:03d}.xyz'.format(filename, i),data,delimiter='\t', fmt=['%d','%f','%f','%f'], comments='', header='{}\n# {}'.format(data.shape[0],comment))
        
        # Calculate Euclidean distance between atoms with number atomnum
        atomnum = 26
        ccoords = data[np.where(data[:,0]==atomnum),1:]
        distances.append(np.linalg.norm(ccoords[0,0,:]-ccoords[0,1,:]))
    
    # Output file with distances
    with open('{}_distances.txt'.format(filename), 'w') as distfile:
        for i, dist in enumerate(distances):
            distfile.write('{}\t{}\n'.format(i,dist))
    
    # Output plot with distances
    plt.plot(distances,linestyle='None', marker='.', markersize=5, color='g')
    plt.xlabel('Geometric optimization steps')
    plt.ylabel('Fe-Fe separation')
    plt.savefig('{}_distances.pdf'.format(filename),bbox_inches='tight')
    
    print('Done')

if __name__ == '__main__':
    _, inputpath, comment = sys.argv
    
    log2xyz(inputpath,comment)