Generate conformers with normal mode sampling

gen-conformers.py

#!/usr/bin/env python

from __future__ import print_function

import sys
import logging
import numpy as np
from cclib.io import ccopen
from nms import nmsgenerator
from masses import element_masses


def nms_ts(n):
    """Normal Mode Temperature and Samples

    Args:
        n (int): The number of non-hydrogen atoms.

    Returns:
        temperature (float): thermal energy to populate normal normal_modes
        samples     (int):   number of samples to consider (ignored)
    """
    stable =  {1: (2000.0, 500),
               2: (1500.0, 450),
               3: (1000.0, 425),
               4: ( 600.0, 400),
               5: ( 600.0, 200),
               6: ( 600.0,  30),
               7: ( 600.0,  20),
               8: ( 450.0,   5),
               9: ( 450.0,   3)}
    if n in stable: return stable[n]
    else: return (400, 2)


for arg in sys.argv[1:]:
    # read the file with cclib
    file = ccopen(arg)
    file.logger.setLevel(logging.ERROR)
    molecule = file.parse()

    labels = molecule.atomnos
    try:
        mass = np.asarray(molecule.atommasses)
    except AttributeError:
        mass = np.asarray([element_masses[i] for i in labels]).repeat(3)

    # get the last geometry from a geometry optimization
    xyz = np.asarray(molecule.atomcoords[-1]).reshape(-1,3)

    # vibrational frequencies (remove any zeros, e.g., translations, rotations)
    freq = np.asarray(molecule.vibfreqs)
    freq = freq[freq !=0.0]

    nmo = np.asarray(molecule.vibdisps).reshape(-1, xyz.shape[0]*3).T
    nmo = nmo[np.array([not (i == 0.0).all() for i in nmo])]

    print(freq.shape, nmo.shape, mass.shape)
    # e.g., (18,) (24, 18) (24,)

    # force constants
    fcc = freq * freq * (nmo * nmo * mass).sum(axis=1) / 16.9744 / 10e4

    t, s = nms_ts(len(labels)-labels.count(1))
    g = nmsgenerator(xyz, nmo, fcc, labels, t, 0)

    #while i < 100:
    e, rxyz = g.get_random_structure(100)