Generating the displaced structures from a Gaussian frequencies calculation

Requirements:

• A vibrational frequencies calculation on a well optimized structure

The vibrational frequencies should have the highest precision possible as the displacements made are small to remain in the harmonic approximztion. In gaussian this can be done from the normal output with the Freq=(HPModes) or with the formatted checkpoint file. The frequencies must be saved to the checkpoint file which is not a default with Freq=(SaveNormalModes).

from vibrav.util.io import uncompress_file
from vibrav.base import resource
import os

This next step is only to decompress the resource output file. In most user cases this can be skipped as the output files are not compressed binaries.

decomp = uncompress_file(resource('g16-nitromalonamide-freq.out.xz'))

Use the quantum code parser of choice and get the atom, frequency and extended frequency data frames. In this case we are using the gaussian parser from the Exatomic package.

from exa import logging
logging.disable()
from exatomic import gaussian

ed = gaussian.Output(decomp)
ed.parse_atom()
ed.parse_frequency()
ed.parse_frequency_ext()

Parsing frequency normal modes from HPModes output
Parsing frequency normal modes from HPModes output

Remove the uncompressed file for some clean up.

os.remove(decomp)

Generating the displaced structures

from vibrav.util.gen_displaced import Displace

inputs = Displace(cls=ed)

The disp class attribute holds the coordinates of all the displaced coordinates that are generated. The frequency indeces are as follows: - 0 is reserved for the equilibrium structure - From 1 up to and including the number of normal modes (39 in this example) are the ones displaced in the positive direction. - From the number of normal modes plus 1 (40 in this example) up to and including twice the number of normal modes (78 in this example) are the negative displacements.

It should be mentioned that the positive and negative displacements are completely arbitrary. Meaning, that we only multiply the normal modes by +1 or -1 for the positive and negative displacements, respectively.

inputs.disp

	x	y	z	freqdx	Z	symbol	frequency	frame
atom
0	0.213964	-4.244614	0.0	0	1	H	0.0000	0
1	-5.885437	-2.178947	0.0	0	1	H	0.0000	0
2	-4.921215	1.041846	0.0	0	1	H	0.0000	0
3	6.295407	-0.573277	0.0	0	1	H	0.0000	0
4	4.506458	2.260056	0.0	0	1	H	0.0000	0
...	...	...	...	...	...	...	...	...
1180	-1.771898	-4.006233	0.0	78	8	O	3696.3886	78
1181	-0.358161	2.748341	0.0	78	7	N	3696.3886	78
1182	1.505387	4.175983	0.0	78	8	O	3696.3886	78
1183	-2.542736	3.630069	0.0	78	8	O	3696.3886	78
1184	-4.537880	-0.833527	0.0	78	7	N	3696.3886	78

1185 rows × 8 columns

Generate the ‘xyz’ coordinate files to view later if something goes wrong.

xyz_dir = 'xyz'
if not os.path.exists(xyz_dir):
    os.mkdir(xyz_dir)
for frame in range(inputs.disp.nframes):
    filename = 'nitromal-{:03d}.xyz'.format(frame)
    xyz_file = os.path.join(xyz_dir, filename)
    with open(xyz_file, 'w') as fn:
        comments = "{:03d} displacement for nitromalonamide".format(frame)
        fn.write(inputs.disp.to_xyz(header=True, comments=comments,
                                    frame=frame))

Define the templates for the gaussian input files to be generated.

grad_template = '''\
%Mem={mem}
%Chk={chk}
%NProc={nproc}
#P Force B3LYP/6-311++G** SCF=Tight Int=UltraFine NoSymm

{comment}

{charge} {mult}
{coord}
'''
prop_template = '''\
%Mem={mem}
%Chk={chk}
%NProc={nproc}
#P PBE1PBE/cc-pVDZ NMR SCF=Tight Int=UltraFine NoSymm

{comment}

{charge} {mult}
{coord}
'''

Generate the gradient inputs with a B3LYP/6-311++G** functional and basis.

input_dir = 'input'
if not os.path.exists(input_dir):
    os.mkdir(input_dir)
for frame in range(inputs.disp.nframes):
    filename = 'nitromal-grad-{:03d}.inp'.format(frame)
    grad_file = os.path.join(input_dir, filename)
    with open(grad_file, 'w') as fn:
        comment = "Gradient calculation for the {:03d} ".format(frame) \
                  +"displacement of nitromalonamide"
        nproc = 2
        mem = "1GB"
        chk = "chk-{:03d}.chk".format(frame)
        charge = 0
        mult = 1
        coord = inputs.disp.to_xyz(header=False, frame=frame)
        fn.write(grad_template.format(mem=mem, chk=chk, nproc=nproc,
                                      comment=comment, charge=charge,
                                      mult=mult, coord=coord))

Generate the NMR property inputs with a PBE1PBE/cc-pVDZ functional and basis.

input_dir = 'input'
if not os.path.exists(input_dir):
    os.mkdir(input_dir)
for frame in range(inputs.disp.nframes):
    filename = 'nitromal-prop-{:03d}.inp'.format(frame)
    prop_file = os.path.join(input_dir, filename)
    with open(prop_file, 'w') as fn:
        comment = "Property calculation for the {:03d} ".format(frame) \
                  +"displacement of nitromalonamide"
        nproc = 2
        mem = "1GB"
        chk = "chk-{:03d}.chk".format(frame)
        charge = 0
        mult = 1
        coord = inputs.disp.to_xyz(header=False, frame=frame)
        fn.write(prop_template.format(mem=mem, chk=chk, nproc=nproc,
                                      comment=comment, charge=charge,
                                      mult=mult, coord=coord))

Make sure that the output directory is present.

output_dir = 'output'
if not os.path.exists(output_dir):
    os.mkdir(output_dir)