Compute polarizability after each step of relaxed scan in Gaussian 16

Question

I am computing a relaxed 2D scan with Gaussian 16 in order to reproduce some computational data. I am interested in computing the dipole moment vector and the (static) polarizability matrix after each geometry optimization. Basically this is what Orca does automatically but I have to use Gaussian since I want to use the PM6 method. I read the manual and the IOPs but to be honest I did not understand if it is doable.

Alternatively, I would be happy if I could find a way to extract the geometries of the relaxed scan done in g16, in this way I could write a script to launch several single points to extract the data that I need.

Lastly, any suggestion about an alternative software (possibly free) that does this kind of calculations using the PM6/PM7 level of theory is welcome.

Thanks in advance for all the help.

Answer 1

This python script will read the logfile from the scan and generate the corresponding gaussian input files.

import  numpy   as  np
import  sys     as  sys
import  typing  as  typing

class Coordinates:
    def __init__(self, natoms:int, coords:list, atoms:list) -> None:
        self.natoms =   natoms
        self.coords =   coords
        self.atoms  =   atoms
        print(self.atoms)

def Gaussian_writeCOM(scanpoint:float, identifier:str, coordinates: Coordinates) -> None:
    with open("coord.{:0.0f}.com".format(scanpoint), "w") as comWriter:
        comWriter.write("%mem = 16GB\n%nprocs = 12\n%chk = {}.{:0.0f}.chk\n".format(identifier, scanpoint))
        comWriter.write("\n")
        comWriter.write("# mp2/cc-pvtz scf = qc\n")
        comWriter.write("\n")
        comWriter.write("single-point energy of {} at {:0.0f}\n".format(identifier, scanpoint))
        comWriter.write("\n")
        comWriter.write("0 1\n")
        for index, coord in enumerate(coordinates.coords):
            comWriter.write(" {:<10s}{:15.5f}{:15.5f}{:15.5f}\n".format(coordinates.atoms[index], coord[0], coord[1], coord[2]))

identifier  =   sys.argv[1]
logFile     =   sys.argv[2]
nAtoms      =   sys.argv[3]
initScan    =   sys.argv[4]
nPoints     =   sys.argv[5]
stride      =   sys.argv[6]

with open(logFile) as log:
    fullData        =   log.readlines()
    initscan        =   float(initScan)
    strideValue     =   float(stride)
    npoints         =   int(nPoints)
    natoms          =   int(nAtoms)
    scanPoints      =   np.arange(initscan, initscan + (strideValue * npoints), strideValue)
    coordObjects    =   []

    atoms   =   []
    # First find the atom names
    for lnumber, line in enumerate(fullData):
        if "Symbolic Z-matrix:" in line:
            atomData    =   fullData[lnumber+2 : lnumber+2+natoms]
            for atomline in atomData:
                atomline    =   atomline[0:2].split()[0]
                atoms.append(atomline)
                print(atomline)

    # Now read the coordinates
    for lnumber, line in enumerate(fullData):
        if "Input orientation" in line:
            coordData   =   fullData[lnumber+5 : lnumber+5+natoms]
        if "Stationary point found" in line:
            coords  =   []
            for line in coordData:
                line    =   line.split()
                coord   =   [float(x) for x in line[3:]]
                coords.append(coord)
            coordObject =   Coordinates(natoms=natoms, coords=coords, atoms=atoms)
            coordObjects.append(coordObject)

for index, coordObject in enumerate(coordObjects):       
    scanPoint   =   scanPoints[index] 
    Gaussian_writeCOM(scanpoint=scanPoint, identifier=identifier, coordinates=coordObject)

You will need to provide the total number of scan points as the input. Rest of the keywords are self-explanatory. Please change the route section to match your needs.