Make input file for geomeTRIC #openbabel #psi4 #geomeTRIC

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) -&gt; target: 0.0003 -&gt; 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.

This image has an empty alt attribute; its file name is fig.png

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.

https://github.com/iwatobipen/qchem

Advertisement

Published by iwatobipen

I'm medicinal chemist in mid size of pharmaceutical company. I love chemoinfo, cording, organic synthesis, my family.

2 thoughts on “Make input file for geomeTRIC #openbabel #psi4 #geomeTRIC

  1. 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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d bloggers like this: