geomeTRIC is code for molecular structures geometry optimization. TRIC means translation-rotation-internal coordinate (TRIC) system. As described in their publication, geomeTRIC can optimize molecular geometry rapidly and few iteration numbers.
It seems useful for Quantum Chemistry. And geomeTRIC supports many QC engines Q-Chem, TeraChem, Psi4 andMolpro. I like psi4 so I tried to use geomeTRIC with psi4. However to use the package native psi4 input file is required, so I can’t use psikit. So I make input generator for geomeTRIC.
To do it, I use pybel and jinja because pybel can calculate not only charge of molecule but also spin.
And jinja2 is used for template enjine.
My script is very simple as below.
#psi4inpgen.py
import pybel
from jinja2 import Template
from argparse import ArgumentParser
tmpl = Template("""molecule {
{{FC}} {{MP}}
{{XYZ}}
}
set basis {{BASIS}}
gradient('{{METHOD}}')
""")
def getArgs():
parser = ArgumentParser()
parser.add_argument("smi", type=str, help="SMILES strings of inputmolecule")
parser.add_argument("--basis", type=str, help='BASIS SETS', default="sto-3g")
parser.add_argument("--method", type=str, help="DFT method", default="hf")
parser.add_argument("--output", type=str, help="outputfile name", default="input")
return parser
if __name__=="__main__":
parser = getArgs()
args = parser.parse_args()
obmol = pybel.readstring("smi", args.smi)
obmol.make3D()
xyz = obmol.write('xyz').split("\n")[2:]
xyz = "\n".join(xyz)
output = tmpl.render(FC=obmol.charge,
MP=obmol.spin,
XYZ=xyz,
BASIS=args.basis,
METHOD=args.method)
with open(f"{args.output}.in", "w") as f:
for line in output:
f.write(line)
input.in file is generated just typing following command from bash.
$python psi4inpgen.py "CO" # methanol as an example
cat input.in
molecule {
0 1
C 0.93502 -0.05907 0.03690
O 2.35068 -0.04556 0.04204
H 0.57490 -0.85612 0.69201
H 0.57490 -0.21982 -0.98222
H 0.56919 0.90323 0.40263
H 2.64600 -0.91154 -0.28708
}
set basis sto-3g
gradient('hf')
Now ready to geometory-optimization. I used –psi4 option for calculation. –verbose and –qdata are optional arguments, they aren’t necessary parameters.
After running code, I could get optimized structure as XYZ format file and it can read from pymol.
$geometic-optimize --psi4 input.in --verbose --qdata
<pre><code> ())))))))))))))))/
())))))))))))))))))))))))),
*)))))))))))))))))))))))))))))))))
#, ()))))))))/ .)))))))))),
#%%%%, ()))))) .))))))))*
*%%%%%%, )) .. ,))))))).
*%%%%%%, ***************/. .)))))))
#%%/ (%%%%%%, /*********************. )))))))
.%%%%%%# *%%%%%%, *******/, **********, .))))))
.%%%%%%/ *%%%%%%, ** ******** .))))))
## .%%%%%%/ (%%%%%%, ,****** /)))))
%%%%%% .%%%%%%# *%%%%%%, ,/////. ****** ))))))
#% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,)))))
#%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).</code></pre>
#%%%%. %%%%%# /%%%%%%* #%%%%%% /////) <strong><em><strong>* ))))),
#%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /</strong></em></strong> ))))).
## %%% .%%%%%%/ <em>%%%%%%, ////////. <strong></strong></em>/ ,)))))
#%%%%# /%%%%%%/ (%%%%%% /)/)// <strong><em></em></strong> ))))))
## .%%%%%%/ (%%%%%%, <strong><em>*</em></strong> ))))))
.%%%%%%/ *%%%%%%, <strong>. /<strong><em>*</em></strong> .))))))
*%%%%%%/ (%%%%%% <strong><em><em></em></em></strong>/<em>..,</em>/<strong><em>*</em></strong></strong> *))))))
#%%/ (%%%%%%, <strong><em><em></em></em></strong><em><strong><em></em></strong></em><strong><em>*</em></strong>/ )))))))
*%%%%%%, ,<strong><em><em></em></em></strong><em><strong><em></em></strong></em>/ ,))))))/
(%%%%%% () ))))))))
#%%%%, ()))))) ,)))))))),
#, ()))))))))) ,)))))))))).
()))))))))))))))))))))))))))))))/
())))))))))))))))))))))))).
())))))))))))))),
-=# geomeTRIC started. Version: 0.9.7.2 #=-
geometric-optimize called with the following command line:
/home/takayuki/anaconda3/envs/chemo37/bin/geometric-optimize --psi4 input.in --verbose --qdata
18 internal coordinates being used (instead of 18 Cartesians)
Internal coordinate system (atoms numbered from 1):
Distance 1-2
Distance 1-3
Distance 1-4
----snip---
Iter: 1 Err-dQ = 1.93889e-07 RMSD: 2.85573e-04 Damp: 1.00000e+00
Cartesian coordinates obtained after 1 microiterations (rmsd = 2.856e-04 |dQ| = 1.939e-07)
dy(i): 0.0009 dy(c) -> target: 0.0003 -> 0.2000
Step 5 : Displace = 2.617e-04/3.489e-04 (rms/max) Trust = 2.000e-01 (+) Grad = 4.562e-05/7.685e-05 (rms/max) E (change) = -113.5493994088 (-8.227e-08) Quality = 0.730
Converged! =D
<pre><code>#==========================================================================#
#| If this code has benefited your research, please support us by citing: |#
#| |#
#| Wang, L.-P.; Song, C.C. (2016) "Geometry optimization made simple with |#
#| translation and rotation coordinates", J. Chem, Phys. 144, 214108. |#
#| http://dx.doi.org/10.1063/1.4952956 |#
#==========================================================================#</code></pre>
Time elapsed since start of run_optimizer: 4.998 seconds
$ cat input_optim.xyz
6
Iteration 0 Energy -113.54722876
C 0.9350200000 -0.0590700000 0.0369000000
O 2.3506800000 -0.0455600000 0.0420400000
H 0.5749000000 -0.8561200000 0.6920100000
H 0.5749000000 -0.2198200000 -0.9822200000
---snip---
H 0.5804872645 0.8928197117 0.3986730191
H 2.6495111809 -0.9134121048 -0.2877923937
6 iteration data is shown as 6 states in pymol.
To use the output file I need to parse the file. At first I search code for using geomeTRIC and how to parse or handle the data.
I uploaded today’s code on my repository.
Is life worth living? 
Dear Iwatobipen, I used to follow your blog and I found it very helpful. I would like to suggest a correction in your code on this page. In your python script “psi4inpgen.py” if you replace the “n” –> “\n” then everything work smoothly. I guess you are trying to split and join using new line as a separator. Hope it is helpful.
Thanks again for sharing your work.
Regards,
Piyush Agrawal
Hi,
Nice catch! It’s just typo. But my github repo used ‘\n’.
Anyway, I just fixed my typo.
Thanks.