Use SAPT/FISAPT more easily #chemoinformatics #psikit #quantumchemistry

Now I’m developing new method for psikit. So this new function is not imlemented in original repository yet but I checked it in my repo.

The new method is named SAPT.

Symmetry-adapted perturbation theory (SAPT) provides a means of directly computing the noncovalent interaction between two molecules. And as same as fragment orbital method(FMO)’s PIEDA (pair interaction energy decomposition) SAPT provides a decomposition of interaction energy: i.e., electrostatic, exchange, induction, and dispersion terms.

And F/ISAPT, functional group or intermolecular SAPT provides an effective two-body partition of the various SAPT terms to localized chemical functional groups (F-SAPT). I think SAPT and especially F-SAPT is useful for chemist however, original psi4 implementation is little bit difficult to make input files.

So I’m trying to implement the function to psikit and make easy to use FSAPT/SAPT.

For SAPT, psi4numpy repository provides useful example for SAPT calculation so most of part borrowed from the repo.

And following code is the SAPT basic usage. For convenience, to run SAPT user need to provide two mol files for calculation which has 3D structure as molfie.

My idea of the implementation, I defined Sapt class which is separately to Psikit class. And code is below.

from psikit import Sapt
sapt=Sapt()
# load monomer 1 and 2 from molfiles
sapt.monomer1_from_molfile('water1.mol')
sapt.monomer2_from_molfile('water2.mol')
# then make dimer geometry format.
# sapt.make_dimer()
# get SAPT0 result. It is easy to do it just call run_sapt()
sapt0, Exch100, Elst10, Disp200, ExchDisp20, Ind20r, ExchInd20r = sapt.run_sapt()
print('total sapt0', sapt0)
print('Exch10 (S^2)', Exch100)
print('Elst10', Elst10)
print('Disp20', Disp200)
print('Exch-Disp20', ExchDisp20)
print('Ind20,r', Ind20r)
print('Exch-Ind20,r', ExchInd20r)

Next example is FSAPT. I had some troubles for the implementation. Because FSAPT of psi4 provides useful post-processing script named fsapt.py but user need to provide fA.dat and fB.dat files which define atom index of the functional groups. At first, I tried to use rdkit pharmacophore modules to make the file. But it was failed because multiple assignment of one atom is not allowed the calculation. So I switched from pharmacophore function to RECAP function. And also I added some helper function.

Finally I checked the function with phenol dimer.

I think the dev version of psikit can perform FSAPT calculation easily.

import psikit
sapt = psikit.Sapt()
sapt.monomer1_from_molfile('phenol1.mol')
sapt.monomer2_from_molfile('phenol2.mol')
sapt.make_dimer()
sapt.run_fisapt()

After calling the run_fisapt(), fsapt folder is generated and all results are stored the folder and fA.dat, fB.dat files are saved in fsapt folder too.

Then go to fsapt folder and call fsapt.py which is post-process script.

fsapt.dat file can see functional group base interactions like below.

=> Full Analysis <=

Frag1 Frag2 Elst Exch IndAB IndBA Disp Total
O_0 O_0 -8.730 6.201 -0.512 -1.546 -1.118 -5.704
O_0 c1ccccc1_0 2.471 0.686 0.256 -0.346 -0.661 2.406
O_0 Link-1 -0.981 0.029 -0.090 -0.122 -0.149 -1.313
c1ccccc1_0 O_0 -2.992 0.670 -0.083 -0.299 -0.558 -3.263
c1ccccc1_0 c1ccccc1_0 1.833 2.115 0.032 -0.245 -2.287 1.449
c1ccccc1_0 Link-1 -1.116 0.123 -0.074 -0.047 -0.117 -1.231
Link-1 O_0 1.261 0.002 -0.050 0.178 -0.107 1.285
Link-1 c1ccccc1_0 -1.514 0.038 0.027 0.107 -0.105 -1.447
Link-1 Link-1 0.673 0.003 -0.008 0.025 -0.019 0.674
O_0 All -7.239 6.916 -0.347 -2.014 -1.927 -4.611
c1ccccc1_0 All -2.275 2.908 -0.124 -0.590 -2.963 -3.045
Link-1 All 0.419 0.043 -0.030 0.310 -0.231 0.512
All O_0 -10.461 6.873 -0.644 -1.667 -1.783 -7.681
All c1ccccc1_0 2.789 2.839 0.315 -0.484 -3.053 2.407
All Link-1 -1.423 0.155 -0.173 -0.144 -0.285 -1.870
All All -9.095 9.868 -0.502 -2.295 -5.121 -7.145

=> Reduced Analysis <=

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

And if you have pymol, launch pymol and call @rum.pymol command from pymol command session, you can get interaction view like image below.

This is still test version not fixed method I need to add testcode and need to review. And after that I would like to Pull Request to original repo.

reference
http://www.psicode.org/psi4manual/master/sapt.html?highlight=sapt
http://www.psicode.org/psi4manual/1.2/fisapt.html