8
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ Maybe a silly question, but have you tried running a small calculation with your scan settings and the keyword Polar? I would think this would be enough to get it print polarizability for each scan step (actually I had though the polarizability and dipole were printed by default even without Polar). $\endgroup$
    – Tyberius
    May 25, 2023 at 13:25
  • $\begingroup$ I tried, unfortunately the polarizability was printed just for the last point $\endgroup$
    – capitn96
    May 25, 2023 at 13:28
  • 1
    $\begingroup$ @captin96, Did the provided answer solve your question? If yes, please indicate it as accepted. If not, please let me know, and I can expand the answer. $\endgroup$ May 26, 2023 at 17:59

1 Answer 1

10
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ +1 and neat work with the code formatting man. I kind of involuntarily guessed that it was mostly you who is answering as soon as I saw the answer. $\endgroup$ May 25, 2023 at 16:24
  • $\begingroup$ Thanks! This is extremely helpful. I agree with @VandanRevanur, the code looks great $\endgroup$
    – capitn96
    May 26, 2023 at 21:57

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .