Use ORBKIT for rendering MO #orbkit #rdkit #psikit #quantum_chemistry

As you know, there many packages for quantum chemistry not only commercial software but also free tools. And each soft has own output format. So user need to understand how to use it. But it is very time consuming step I think.

ORBIKIT is a modular python toolbox for cross-platform post processing of quantum chemical wavefunction data.

The reference URL is below.
https://arxiv.org/abs/1601.03069

It seems nice isn’t it? Fortunately the source code is disclosed on github.

So I tried to use orbkit with psikit. At first, I installed orbkit while referring the original site.

At first, I installed cython, numpy, scipy, h5py with conda and mayavi installed with pip (mayavi couldn’t install with conda due to package conflicts).

Then install orbkit.

$ cd $HOME
$ git clone https://github.com/orbkit/orbkit.git
$ export ORBKITPATH=$HOME/orbkit
$ cd $ORBKITPATH
$ python setup.py build_ext --inplace clean
$ export PYTHONPATH=$PYHONPATH:$ORBKITPATH

And also I added environment variables to my .bashrc.

Now ready, let’s use orbkit. I used psikit for making fchk file for orbkit. orbkit supports many kinds of input format. If reader has interest the package, please check the original document.

Following code shows how to use psikit and visualize MO with orbkit.

After running the code, mayavi viewer launched and could get MO view, the image is shown below.

Of course current psikit has some functions for communicate pymol for rendering MO but orbkit is also useful package for quantum chemistry.

I uploaded today’s code on my github repo.

https://github.com/iwatobipen/playground/blob/master/orbkit_psikit_test.ipynb

Rendering molecular orbital on Jupyter notebook #psikit #py3dmol #rdkit #memo

@fmkz___ and I( @iwatobipen ) are developing psikit which is a thin wrapper of psi4 and rdkit. I hope the package integrates quantum chemistry (Psi4) and chemoinformatics (RDKit).

By using psikit, user can make molecular orbital data very convinienlry.

Rendering MO is useful for understanding molecular electrostatic shape and nature, but sometime it is difficult to medicinal chemist because it requires CADD skills for rendering it.

Psikit has useful method named ‘create_cube_files’, the method make cube files of HOMO, LUMO and ESP. And also psikit can call pymol to rendering them. However it can’t render MO directory on jupyter notebook so I would like to show you how to render MO on jupyter notebook today.

Fortunately, py3Dmol is easy tool for rendering CUBE file!

Sample code is below. I used psikit for quantum calculation and used py3dmol for rendering. First, calculate energy and create cube files. Default gridspace value is set 0.3 but I recommend it should be changed more large r value. Because cube file with small gridspace value will become huge size file.

import py3Dmol
from psikit import Psikit
from rdkit import Chem
pk = Psikit()
pk.read_from_smiles('O')
pk.energy()
pk.create_cube_files(gridspace=0.5)

After running the method, cube files are generated. So read the file.

homo_voldata = open("Psi_a_5_1-A\"_HOMO.cube", "r").read()
lumo_voldata = open("Psi_a_6_5-A\'_LUMO.cube", "r").read()

Now ready to render. Make view object and attach some dataset. addVolumetricData method loads cube data with options. Then call addModel with molblock data. RDKit can convert molobject to molblock by using MolToMolBlock method.

v = py3Dmol.view()
v.addVolumetricData(homo_voldata, "cube", {'isoval': -0.03, 'color': "red", 'opacity': 0.75})
v.addVolumetricData(homo_voldata, "cube", {'isoval': 0.03, 'color': "blue", 'opacity': 0.75})
v.addModel(Chem.MolToMolBlock(pk.mol), 'mol')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

Works fine. The molecule object is rotatable. LUMO can be rendered with same manner.

In summary, psi4, rdkit and py3dmol integration can render molecule with MO at the touch of a button. I uploaded today’s code on my repository.

https://github.com/iwatobipen/playground/blob/master/mo/RenderMO.ipynb

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>*      ))))),
#%%%%##%  %%%#      .%%%%%%/      (%%%%%%,  ///////.     /</em></strong></strong>      ))))).
##     %%%      .%%%%%%/      <em>%%%%%%,  ////////.      <strong></em></strong>/     ,)))))
#%%%%#      /%%%%%%/      (%%%%%%      /)/)//       <strong><em></strong></em>      ))))))
##      .%%%%%%/      (%%%%%%,                  <strong><em>*</em></strong>      ))))))
.%%%%%%/      *%%%%%%,  <strong>.             /<strong><em>*</em></strong>      .))))))
*%%%%%%/      (%%%%%%   <strong><em><em></em></strong></em>/<em>..,</em>/<strong><em>*</em></strong></strong>       *))))))
#%%/      (%%%%%%,    <strong><em><em></em></strong><strong><em></em></em></strong><strong><em>*</em></strong>/        )))))))
*%%%%%%,         ,<strong><em><em></em></strong><strong><em></em></em></strong>/         ,))))))/
(%%%%%%   ()                              ))))))))
#%%%%,  ())))))                        ,)))))))),
#,    ())))))))))                ,)))))))))).
()))))))))))))))))))))))))))))))/
())))))))))))))))))))))))).
())))))))))))))),

-=#  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

Example usage of psi4-openmm-interface #Psi4 #OpenMM #RDKit

Molecular dynamics and Quantum Chemistry are important tools for CADD. I have interested in these topics and OpenMM and Psi4 are nice tool to handing MD and QM.

Today I tried to use psi4-openmm-interface which allows passing of molecular systems between each program.

I reviewed test script and found that the package pass the molecule which optimized with openmm to psi4 as xyz format. And also the package can conduct MD calculation very simply.

I used it, following code is almost same as original repo.

At first, I install the package you can find how to install it in following URL.
https://github.com/mzott/Psi4-OpenMM-Interface

Install psi4, openmm and set PYTHONPASS. It was simple. And I need to modify some code to import the package correctly.

  1. It uses PeriodicTable class of qcelemental, but the module can’t import _temp_symbol attribute from periodictable.py, so I added _temp_symbol line to qcelemental/periodic_table.py
# qcelemental/periodic_table.py 
class PeriodicTable:
    """Nuclear and mass data about chemical elements from NIST.

    Parameters
    ----------
    None

    Attributes
    ----------
    A : list of int
        Mass number, number of protons and neutrons, starting with 0 for dummies.
    Z : list of int
        Atomic number, number of protons, starting with 0 for dummies.
    E : list of str
        Element symbol from periodic table, starting with "X" for dummies. "Fe" capitalization.
    EA : list of str
        Nuclide symbol in E + A form, e.g., "Li6".
        List `EA` is a superset of `E`; that is, both "Li6" and "Li" present.
        For hydrogen, "D" and "T" also included.
    mass : list of :py:class:`decimal.Decimal`
        Atomic mass [u].
        For nuclides (e.g., "Li6"), the reported mass.
        For stable elements (e.g., "Li"), the mass of the most abundant isotope ("Li7").
        For unstable elements (e.g., "Pu"), the mass of the longest-lived isotope ("Pu244").
    name : list of str
        Element name from periodic table, starting with "Dummy". "Iron" capitalization.

    """

    def __init__(self):

        from . import data

        # Of length number of elements
        self.Z = data.nist_2011_atomic_weights["Z"]
        self.E = data.nist_2011_atomic_weights["E"]
        self.name = data.nist_2011_atomic_weights["name"]

        self._el2z = dict(zip(self.E, self.Z))
        self._z2el = collections.OrderedDict(zip(self.Z, self.E))
        self._element2el = dict(zip(self.name, self.E))
        self._el2element = dict(zip(self.E, self.name))

        # Of length number of isotopes
        self._EE = data.nist_2011_atomic_weights["_EE"]
        self.EA = data.nist_2011_atomic_weights["EA"]
        self.A = data.nist_2011_atomic_weights["A"]
        self.mass = data.nist_2011_atomic_weights["mass"]

        self._eliso2mass = dict(zip(self.EA, self.mass))
        self._eliso2el = dict(zip(self.EA, self._EE))
        self._eliso2a = dict(zip(self.EA, self.A))
        self._el2a2mass = collections.defaultdict(dict)
        for EE, m, A in zip(self._EE, self.mass, self.A):
            self._el2a2mass[EE][A] = float(m)
       # following code was added
        self._temp_symbol =  [
    "X", "H", "HE", "LI", "BE", "B", "C", "N", "O", "F", "NE", "NA", "MG", "AL", "SI", "P", "S", "CL", "AR", "K", "CA",
    "SC", "TI", "V", "CR", "MN", "FE", "CO", "NI", "CU", "ZN", "GA", "GE", "AS", "SE", "BR", "KR", "RB", "SR", "Y",
    "ZR", "NB", "MO", "TC", "RU", "RH", "PD", "AG", "CD", "IN", "SN", "SB", "TE", "I", "XE", "CS", "BA", "LA", "CE",
    "PR", "ND", "PM", "SM", "EU", "GD", "TB", "DY", "HO", "ER", "TM", "YB", "LU", "HF", "TA", "W", "RE", "OS", "IR",
    "PT", "AU", "HG", "TL", "PB", "BI", "PO", "AT", "RN", "FR", "RA", "AC", "TH", "PA", "U", "NP", "PU", "AM", "CM",
    "BK", "CF", "ES", "FM", "MD", "NO", "LR", "RF", "DB", "SG", "BH", "HS", "MT", "DS", "RG", "UUB", "UUT", "UUQ",
    "UUP", "UUH", "UUS", "UUO"
]
.....snip...

2. Line 81 of psiomm/molecule.py should be change below.

            # from
            #z_vals.append(periodictable.el2z[Zxyz_list[0].upper()])
            #To 
            z_vals.append(periodictable._el2z[Zxyz_list[0].upper()])

3. gaff2.xml which is forcefield parameters is placed in openmm installed path openmm/app/data/. Because for my env, gaff2.xml file couldn’t find in the folder.

4. Got cov_radii.py from qcdb/qcdb and placed in psi4/driver/qcdb folder.

Now ready. Let’s write code.

Import some packages at first I used psikit for molecular geometry generation.

import numpy as np

from psikit import Psikit
from psiomm import molecule
from psiomm import psi_omm as po

from simtk.openmm.app import Simulation
from simtk.openmm.app import ForceField
from simtk.openmm.app import Modeller
from simtk.openmm.app import PDBReporter

from simtk.openmm import Platform
from simtk.openmm import LangevinIntegrator

from simtk.openmm import Vec3

from simtk.unit import picosecond
from simtk.unit import kelvin
from simtk.unit import kilocalorie_per_mole
from simtk.unit import angstrom

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole

Work flow

  • make molecule from xyz string => solute

  • generate bonds

  • generate atom types

  • generate charges

  • make topology

  • define forcefield

  • generate template with FF, topology and solute

  • make modeller

  • perform simulation to get trajectory

pk = Psikit()
pk.read_from_smiles('CCOc1cccnc1')
pk.rdkit_optimize()

Psikit generate xyz strings with charge and multiplicity information but psiomm can’t read the line so I remove them and passed the string to psiomm. Then generate some molecular properties.

# remove mutiplicity, formal charge
xyzstring = '\n'.join(pk.mol2xyz().split('\n')[2:])
solute = molecule.Molecule.from_xyz_string(xyzstring)
solute.generate_bonds()
solute.generate_atom_types()
solute.generate_charges()
# check generated atoms and charges
for i, at in enumerate(solute.atom_types):
    print(i, at, solute.charges[i])

0 c3 -0.16240415191041357
1 c3 0.03679384783580364
2 os -0.263056064493707
   ---snip---
17 h4 0.07500144690962385

Following workflow is based on openmm. GAFF2 is implemented in AMBER and TIP3P is water model which is implemented in CHARMM.

topology = po.make_topology(solute)
forcefield = ForceField('gaff2.xml', 'tip3p.xml')
po.generateTemplate(forcefield, topology, solute)
coords = []
for i in range(len(solute.xyz)):
    print(solute.xyz[i])
    coords.append(Vec3(solute.xyz[i][0], solute.xyz[i][1], solute.xyz[i][2]))
positions = coords * angstrom

Make modeller with topology and positions. I added 50 water molecules to the model. And defined the simulation system.

modeller = Modeller(topology, positions)
modeller.addSolvent(forcefield, numAdded=50)
omm_sys = forcefield.createSystem(modeller.topology)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picosecond)

simulation = Simulation(modeller.topology,
                       omm_sys,
                       integrator,
                       platform=Platform.getPlatformByName('Reference'))

I used PDBReporter module to store the simulation. Simulation step set 400 and reporter object save data each 10 step. So I could get output pdb with 40 (400 /10 ) states.

simulation.context.setPositions(modeller.positions)
simulation.reporters.append(PDBReporter('output.pdb', 10))
simulation.step(400)
simulation.minimizeEnergy(tolerance=2*kilocalorie_per_mole)
MD movie

It worked fine. following code didn’t use psi4 based calculation, to do it I need to modify input file because last data has not only target molecule but also 50 water molecules. I’ll check the package and psi4 more.

Today’s code can find from my gist/github repo.

https://github.com/iwatobipen/playground/blob/master/openmm-psikit.ipynb

Call psi4 from web app! #Psi4 #QuantumChemistry

Psi4 is a python package for quantum chemistry. I like it. Native psi4 has some difficulties for preparing input file. So psikit is useful I think. And I think webmo is useful package for calling psi4. I’ve always been interested webmo. So I tried to use it today.
Basic function of webmo is provided free licence.

At first I registered webmo for getting licence code. https://www.webmo.net/index.html

After registration, I could get webmo download link, licence number, and password.

Second, I installed several packages with apt command.

apt install apache2 apache2-suexec-custom libcgi-pm-perl

Then create directory structure.

$ cd ~
$ mkdir public_html
$ chmod 755 public_html
$ cd public_html
$ mkdir cgi-bin
$ chmod 755 cgi-bin
$ exit
Start apache
# a2enmod userdir
# a2enmod suexec
# a2enmod cgid
# service apache2 start

Then eddit /etc/apache2/mods-available/userdir.conf, to enable cgi script.


        UserDir public_html
        UserDir disabled root

        
                AllowOverride FileInfo AuthConfig Limit Indexes
                Options MultiViews Indexes SymLinksIfOwnerMatch IncludesNoExec ExecCGI #ExecCGI added
                AddHAndler cgi-script .cgi #Added
                Require method GET POST OPTIONS
        


Then restart apache ;)

# service apache2 restart

Then install webmo with perl script.

$ tar xzf WebMO.9.x.xxx.tar.gz
$ cd WebMO.install
$ perl setup.pl

Now ready, after restarting apache, I could access local webmo HP via web browser. Default setting is username is admin, password is blank.

After login as admin, system asked to change admin password. So I changed and the create user for calling Psi4.

sample URL is below
http://localhost/~yourusername/cgi-bin/webmo/login.cgi

You can get more details about installation is following URL.
https://www.webmo.net/support/installation.html

Then from admin top page, I set psi4 as quantum chemistry engine. Click interface tab, change psi4 from disable to enable. Then click edit button and set path for psi4 execution.

Tips, If reader installed psi4 via conda, the path schould be ‘/home/***/anaconda/envs/yourenv’. In the case, psi4 path is /home/***/anaconda/envs/yourenv/bin/psi4. bin/psi4 should not be included setting path!!!

After completing the configuration, I made user for QM calculation and login the system as a user.

Making input file is very easy, import file or draw molecule with implemented molecular editor. I used editor for the test.

After drawing the molecule, minimize energy and symmetrize(if possible) is done from selecting upper menu bar. It’s very easy isn’t it?

Then set job options which defines calculation method, basis function, reference and several QM related options. Also the setting can do from GUI.

After setting the options, input data can check from preview tab.

All setting is done, click submit job to run the calculation.

Following screenshot shows water as an example because benzoic acid required looooooooooooooong calculation time in my hobby PC.

Result can see jobmanager. I calculated water molecule vibration. So I could see calculated IR values and spectra.

The spectrum is not so accurate, the reason is that I used low level basis set I think.

Of course psi4 raw output file can download from webserver.

With free charge licence, it is not able to draw MO. But with commercial licence, rendering MO of molecule is available.

WebMO is very useful tool for education, and wet researchers. Because it is easy to use and rich GUI.

If reader who has interested in WebMO, try it! ;)

Calculate solvent effect in Psi4 #psi4 #quantumchemistry

Recently I use not only chemoinformatics tools but also quantum chemistry tool, my favorite is Psi4. Psi4 has many options and plug-ins for quantum calculation. Most setting of calculation is vacuum, but it actually true. So considering the solvent around the molecules is important.

Can psi4 perform calculation with solvent effect?

Yes!

PCMSolver is plugin for PCM calculation. It can be installed via conda.

$ conda install -c psi4 pcmsolver

PCM means ‘polarizable continuum model’. PCM describes solute-solvent interactions only by means of electrostatics and polarization between the solute molecule and solvent. The solvent is modeled as a dielectric continuum with a given permittivity ε.

To use PCMsolver in psi4, user need to make specific input file for the calculation which should has kind of solvent, basis function, CPCM or IEFPCM, etc.

I made to input file from toulene.mol file which is retrieved from ChEMBL.

input files are below. PCMsolver can not call from python interpreter so I need to make input.dat.

# toluene in water
#! pcm

molecule mol {

0 1
C -2.0705902285428452 0.09342117730391691 0.21438842743524625
C 0.06881622580040792 1.2107631870687914 0.004377262725162135
C -0.043709196272760036 -1.211640074147377 -0.026831077837745607
C -1.3202142335884217 1.271231179761799 0.15648583022799556
C 0.714364020910022 -0.03184758851751915 -0.10370200526939809
C 2.2076605753403253 -0.0996279156897954 -0.22684462930707974
C -1.4325455046382927 -1.1469904802166477 0.1253309895373383
H -3.1452716192643337 0.14180771794289201 0.3335057008466673
H 0.6396267675023236 2.1305648085368674 -0.0295934787763336
H 0.4394292831589224 -2.1792111515476216 -0.08511724686650844
H -1.8144108252945523 2.2311166905139106 0.23329105314598383
H 2.59680405039162 0.7889905119305699 -0.7681174009541791
H 2.5134086150024317 -1.0062726983460601 -0.7912601501934461
H 2.660308829389096 -0.1337122631434785 0.78606242092797
H -2.013676759894258 -2.058593101450173 0.17802430435839692
no_reorient
no_com
units angstrom

}

set {
basis 6-31G**
scf_type pk
pcm true
pcm_scf_type total
}

pcm = {
Units = Angstrom
Medium {
SolverType = IEFPCM
Solvent = Water
}

Cavity {
RadiiSet = UFF
Type = GePol
Scaling = false
Area = 0.3
Mode = Implicit
}
}

energy_pte, wfn = energy(‘scf’, return_wfn=True)

# toluene in benzene
#! pcm

molecule mol {

0 1
C -2.0705902285428452 0.09342117730391691 0.21438842743524625
C 0.06881622580040792 1.2107631870687914 0.004377262725162135
C -0.043709196272760036 -1.211640074147377 -0.026831077837745607
C -1.3202142335884217 1.271231179761799 0.15648583022799556
C 0.714364020910022 -0.03184758851751915 -0.10370200526939809
C 2.2076605753403253 -0.0996279156897954 -0.22684462930707974
C -1.4325455046382927 -1.1469904802166477 0.1253309895373383
H -3.1452716192643337 0.14180771794289201 0.3335057008466673
H 0.6396267675023236 2.1305648085368674 -0.0295934787763336
H 0.4394292831589224 -2.1792111515476216 -0.08511724686650844
H -1.8144108252945523 2.2311166905139106 0.23329105314598383
H 2.59680405039162 0.7889905119305699 -0.7681174009541791
H 2.5134086150024317 -1.0062726983460601 -0.7912601501934461
H 2.660308829389096 -0.1337122631434785 0.78606242092797
H -2.013676759894258 -2.058593101450173 0.17802430435839692
no_reorient
no_com
units angstrom

}

set {
basis 6-31G**
scf_type pk
pcm true
pcm_scf_type total
}

pcm = {
Units = Angstrom
Medium {
SolverType = IEFPCM
Solvent = Benzene
}

Cavity {
RadiiSet = UFF
Type = GePol
Scaling = false
Area = 0.3
Mode = Implicit
}
}

energy_pte, wfn = energy(‘scf’, return_wfn=True)

PCMsolver can has many kinds of solvents listed below.

    1:'Water',
    2:'Propylene Carbonate',
    3:'Dimethylsulfoxide',
    4:'Nitromethane',
    5:'Acetonitrile',
    6:'Methanol',
    7:'Ethanol',
    8:'Acetone',
    9:'1,2-Dichloroethane',
    10:'Methylenechloride',
    11:'Tetrahydrofurane',
    12:'Aniline',
    13:'Chlorobenzene',
    14:'Chloroform',
    15:'Toluene',
    16:'1,4-Dioxane',
    17:'Benzene',
    18:'Carbon Tetrachloride',
    19:'Cyclohexane',
    20:'N-heptane',
    21:'Explicit' 

After making the input files, following step is as same as default psi4 calculation.

$ psi4 input.dat -o output.dat

After waiting few seconds, I could get two output results. Picked up line 355-363.

## toulene in water
@RHF Final Energy: -269.75544369636918

=> Energetics Energetics <=

Nuclear Repulsion Energy = 268.6733324038011688
One-Electron Energy = -895.7789138809855558
Two-Electron Energy = 357.3548941675961714
PCM Polarization Energy = -0.0019460266422257
Total Energy = -269.7526333362304172

Computation Completed

There is difference in PCM Polarization Energy between in water and in toluene, in water it is -0.0052558561183017A.U.(-3.298kcal/mol) and in tolunene, it is -0.0019460266422257A.U.(-1.22kcal/mol)
*1 a.u. =627.50 kcal/mol

I would like to test some molecules. Are there freely available experimental data set? I’ll try to find them.

Added Notebook Example for psikit #psi4 #RDKit

Some days ago. I updated psikit, added new function, SAPT.

The method can calculate inter-intra molecular interaction with psi4. Psikit can call the method very easily.

And I added sample notebook in the repository. Following codes are very simple example of SATP for water dimer and FSAPT for phenol dimer.

SAPT for water dimer.



from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from psikit import Sapt

w1 = Chem.MolFromMolFile('water1.mol', removeHs=False)
w2 = Chem.MolFromMolFile('water2.mol', removeHs=False)
Draw.MolsToGridImage([w1,w2])
sapt = Sapt()
sapt.monomer1_from_molfile('water1.mol')
sapt.monomer2_from_molfile('water2.mol')
sapt.make_dimer()
res = sapt.run_sapt()
---stdout is below
Initializing SAPT object...

RHF for monomer A finished in 0.55 seconds.
RHF for monomer B finished in 0.52 seconds.
Building ERI tensor...
...built ERI tensor in 3.180 seconds.
Size of the ERI tensor is 0.36 GB, 82 basis functions.

...finished initializing SAPT object in  4.49 seconds.

Starting electrostatics...
...electrostatics took a total of  0.17 seconds.

Starting exchange...
...exchange took a total of  0.63 seconds.

Starting dispersion...
...dispersion took a total of  8.92 seconds.

Starting induction...
Ind20,r (AB)           -1.37852223 mH       -0.86503511 kcal/mol
Exch-Ind20,r (AB)       0.88580457 mH        0.55585034 kcal/mol
...induction took a total of  15.90 seconds.

SAPT0 Results
----------------------------------------------------------------------
Exch10 (S^2)             10.53844851 mH        6.61297129 kcal/mol
Elst10                  -13.02830646 mH       -8.17537956 kcal/mol
Disp20                   -3.42996225 mH       -2.15233218 kcal/mol
Exch-Disp20               0.61399531 mH        0.38528758 kcal/mol
Ind20,r                  -4.33068313 mH       -2.71754264 kcal/mol
Exch-Ind20,r              2.30125737 mH        1.44405971 kcal/mol
----------------------------------------------------------------------
Total SAPT0              -7.33525065 mH       -4.60293580 kcal/mol

F-SAPT for phenol dimer

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from psikit import Sapt

p1 = Chem.MolFromMolFile('phenol1.mol', removeHs=False)
p2 = Chem.MolFromMolFile('phenol2.mol', removeHs=False)
Draw.MolsToGridImage([p1,p2])
sapt = Sapt()
sapt.monomer1_from_molfile('phenol1.mol')
sapt.monomer2_from_molfile('phenol2.mol')
sapt.make_dimer()
res = sapt.run_fisapt()

--fA.dat and fB.dat files are generated in fsapt folder the  call fsapt.py from terminal, it generates FSAPT results which is shown below.
   => Full Analysis  Reduced Analysis  F-ISAPT: Links 50-50  Full Analysis  Reduced Analysis <=

Frag1     Frag2         Elst     Exch    IndAB    IndBA     Disp    Total 
c1ccccc1_0 c1ccccc1_0    0.686    2.197    0.007   -0.208   -2.403    0.278 
c1ccccc1_0 O_0         -2.751    0.733   -0.147   -0.227   -0.675   -3.067 
O_0       c1ccccc1_0    1.392    0.720    0.222   -0.347   -0.793    1.194 
O_0       O_0         -8.421    6.218   -0.584   -1.512   -1.250   -5.549 
c1ccccc1_0 All         -2.065    2.930   -0.140   -0.435   -3.078   -2.789 
O_0       All         -7.030    6.938   -0.362   -1.859   -2.043   -4.356 
All       c1ccccc1_0    2.078    2.917    0.229   -0.556   -3.196    1.472 
All       O_0        -11.173    6.951   -0.731   -1.739   -1.925   -8.617 
All       All         -9.095    9.868   -0.502   -2.295   -5.121   -7.145 

Psikit can use SAPT, FSAPT efficiently I think. There is room for improvement for FSAPT fragmentation method. Any comments and advice will be greatly appreciated.

URL is below.

GITHUB
https://github.com/Mishima-syk/psikit
FSAPT
https://nbviewer.jupyter.org/github/Mishima-syk/psikit/blob/master/examples/example_sapt/fsapt_ex.ipynb
SAPT
https://nbviewer.jupyter.org/github/Mishima-syk/psikit/blob/master/examples/example_sapt/sapt_ex.ipynb