How come I can optimize 'Indole' but not '1-Methyl-1H-indole-3-carboxamide'?

[lewiso1]

Hi PSI4, loving the library so far, thankyou.

Im using psi4 to perform SCF calculations on two drug-like molecules to rationalize some SAR observed by a colleague. For example, I have no trouble optimizing and obtaining a wavefunction of the ‘indole’ fragment ( Indole | C8H7N - PubChem ):

%%time

###Indole:
xyz = ['',
 '0 1',
 'C -2.012882216046254 -0.8640647065717614 0.18471956939653567',
 'C -2.113550305859065 0.544030066073846 0.11946599992032274',
 'C -0.9569053634132676 1.3365674981928146 -0.0018746843816307176',
 'C 0.2748791028187773 0.6853861307938023 -0.054357714240697914',
 'C 0.3709789697554246 -0.6796501235447245 0.00902058514200702',
 'C -0.754684697546061 -1.4920257568928068 0.12922262748886626',
 'C 1.708150859298044 -1.0025591528463575 -0.0676433930183273',
 'C 2.3704810170844954 0.21647847914010918 -0.17641565245003507',
 'N 1.4906850171304467 1.2443274977882215 -0.16783569447486663',
 'H -2.9089772127476583 -1.4636367379156243 0.27794436854445626',
 'H -3.086300847912362 1.0162855927526462 0.1630539704840909',
 'H -1.0230192953236514 2.41567452818668 -0.052660479343215216',
 'H -0.6662665264298884 -2.5697144490618453 0.17836028699074521',
 'H 2.1500181910090075 -1.9899991538037467 -0.047790067676360695',
 'H 3.4422125462870814 0.3409808476235725 -0.2575187445878446',
 'H 1.7151807618949582 2.2619194400851836 -0.23569097779405335',
 '']

#prep:
xyz = '\n'.join(xyz)
geom = psi4.geometry(xyz)

#optimize:
basis_sets= "scf/6-31g**"
scf_energy, wfn = psi4.optimize(basis_sets, return_wfn=True, molecule=geom)

which runs in ~30s.

We thought it might be important to represent the whole molecule more completely, so I’ve tried the same settings (plus higher maxiter) with 1-Methyl-1H-indole-3-carboxamide ( 1-Methyl-1H-indole-3-carboxamide | C10H10N2O - PubChem ):

%%time

#altered
xyz = ['',
 '0 1',
 'C -1.2066532557975438 -2.727127027225506 -0.32290182938062345',
 'C -2.355813724113892 -1.9403237111532872 -0.10049798444915797',
 'C -2.24071801484374 -0.5552053998214622 0.10490019311506109',
 'C -0.96185743804853 0.007811167472659488 0.08083063343228829',
 'C 0.1603241320947546 -0.7564922003677576 -0.13713673937758478',
 'C 0.06958623925785978 -2.1372867220949088 -0.3428786512776029',
 'C 1.252122086331245 0.10397622674745564 -0.09860857314657753',
 'C 2.6830329470495466 -0.23077000092910427 -0.28137010859240147',
 'O 3.023717229825934 -1.4134617646431078 -0.5467415133094116',
 'N 3.6982007211206085 0.768791828584701 -0.16320690730501095',
 'C 0.7180906250157293 1.3700913715404193 0.14549876394203531',
 'N -0.6286121548990626 1.3008645440452558 0.2536725140327987',
 'C -1.533870995674511 2.418679093411311 0.5033680487700258',
 'H -1.3079060065158599 -3.793264938102847 -0.47911282529135396',
 'H -3.3325115959251628 -2.406279166926703 -0.08749569431028571',
 'H -3.1205236261193305 0.051116394568693234 0.2762611171896728',
 'H 0.9431371792508623 -2.7496912754980243 -0.5126554486816562',
 'H 3.473700651593924 1.7595396810516233 0.07119750647259149',
 'H 4.700448471191456 0.5114689557053794 -0.30348883953781575',
 'H 1.2766678250610202 2.291853377081106 0.23541362967350726',
 'H -2.1040818013311804 2.236573236451151 1.438458862448881',
 'H -2.2421654072390287 2.522812055703182 -0.3452442062355543',
 'H -0.9643140872898155 3.3663242743914235 0.6117380518186651',
 '']

#prep:
xyz = '\n'.join(xyz)
geom = psi4.geometry(xyz)

#optimize:
basis_sets= "scf/6-31g**"
scf_energy, wfn = psi4.optimize(basis_sets, return_wfn=True, molecule=geom, options={'MAXITER':20})

which exits with OptimizationConvergenceError: Could not converge geometry optimization in 4 iterations.

My guess is it’s just a scale issue- the molecule is four atoms larger, thus takes more iterations to solve. But, a few questions please:

• How come MAXITER is apparently being ignored?
• Is there a functional or basis set better suited to larger molecules? Currently using scf/6-31g** based on rudimentary knowledge
• Is it OK to use the wfn coming from an un-converged optimization? The use case is simply visualizing electrostatic potential surfaces, so I won’t be doing quantitative analysis. To be honest, I’m fine with the RDKit-optimized geometry and a single call to psi4.energy1, but I’d like to know if I’ve gone wrong somewhere first.

Thanks for your time!

-Lewis

edit: psi4.version prints ‘1.3.2’, installed from conda.

[jmisiewicz]

For the first question, you’re probably hitting some other error that causes the optimization to end early. I can’t identify the error without an output file.

I’ll leave the other questions to people who use DFT and electrostatic potentials on a regular basis.

[lewiso1]

Well, turns out it was just me inputting the wrong thing. psi4.optimize does not have an options kwarg, and the options dict should instead be given to psi4.set_options(options) prior to running optimize.

So now it can optimize fine. Thanks for your time