Quantum chemistry calculation with RDKit.

I enjoyed JCUP VII. All presentations were very impressive for me and I interested in psi4 that is open source tool for quantum chemistry.
http://www.psicode.org/
Fortunately, I could install psi4 using conda command. It supports only python2.x.
If user want to install psi4, just type “conda install -c psi4 psi4 “. That’s all !
Ready to run psi4. Next, I need to prepare inputfile.
Basic structure of input file of psi4 is following. (DFT energy calculation)

#! Sample UHF/6-31G** CH2 computation <==comment here

memory 250 mb <== option

molecule ch2 {
  0 1 <= 0 means 0 the molecule is neutral, 1 means the spin mautiplicity
  C position.x position.y position.z <= atom symbole, X, Y, Z
  .....
  
}

set basis 6-31G** <==basis function
set reference uhf <==unrestricted Hartree-Fock 
energy ('scf')

molecule was represented as atom coordination of x, y, z. It is not common structure for me.
So, I wrote simple convert code that transform sdf to input file.
My code is following.(This code translate only first molecule of SDF.)

from __future__ import print_function
import sys, os, argparse
from rdkit import Chem
from rdkit.Chem import AllChem

inf = sys.argv[ 1 ]
sdf = Chem.SDMolSupplier( inf )
mol = sdf.next()
print(mol.GetNumConformers())
if mol.GetNumConformers() <= 1:
    hmol = Chem.AddHs( mol )
    AllChem.EmbedMolecule(  hmol,
                            useExpTorsionAnglePrefs=True,
                            useBasicKnowledge=True )
    #AllChem.EmbedMolecule( hmol )
else:
    hmol = Chem.AddHs( mol )

atoms = [ atom for atom in hmol.GetAtoms() ]


def atomposition2string( atom ):
    line = "{} {} {} {}"
    conf = atom.GetOwningMol().GetConformer()
    posi = conf.GetAtomPosition( atom.GetIdx() )
    line = line.format( atom.GetSymbol(), posi.x, posi.y, posi.z )
    line +="\n\n"
    return line


outf = open('in.dat', 'w')

header="#! psi4input\n\n"
molstring ="molecule inputmol "

setstring = """
set basis 6-31G**\n
set reference uhf\n
energy( "scf" )\n
"""

outf = open('in.dat', 'w')
outf.write(header)

molstring += "{\n0 1\n\n"
for atom in atoms:
    l=atomposition2string(atom)
    molstring += l
molstring += "}\n"

outf.write( molstring )
outf.write(setstring)
outf.close()

All option was hardcoding.

Next I made sample molecule using RDKIT.

from rdkit import Chem
mol = Chem.MolFromSmiles( "Cc1ccccc1" )
w = Chem.SDWriter( "tol.sdf" )
w.write( mol )
w.close()

OK ready, then convert sdf to inputfile.

iwatobipen$ python sdf2psi.py toluene.sdf
iwatobipen$ cat in.dat 
#! psi4input

molecule inputmol {
0 1

C 2.12643946569 -0.148076377601 0.280172726109

C 0.737112657203 -0.0563266799793 -0.220757232163

C 0.0894292148607 1.17309385348 -0.272662802846

C -1.26388345191 1.29962103669 0.0220149512856

C -1.99460856032 0.133837305932 0.199621363415

C -1.42391701302 -1.05115760794 -0.254293865166

C -0.0465313107692 -1.17873499271 -0.410958831529

H 2.06321195027 -0.624503035816 1.29295672814

H 2.54504673868 0.854508459836 0.429328334067

H 2.763215791 -0.825994020449 -0.31296772677

H 0.680672612474 2.01686875154 -0.629481432558

H -1.64020718362 2.2772793939 0.325940161739

H -3.02108403858 0.17949008105 0.519270208312

H -2.0538300122 -1.92491235886 -0.416048311006

H 0.438933140245 -2.12499380908 -0.552134271027

}

set basis 6-31G**

set reference uhf

energy( "scf" )

It seems to work fine! Then run calculation.
Just type, psi4 input.dat output.dat -n

iwatobipen$ time psi4 in.dat out.dat -n 6

real	0m4.868s
user	0m14.404s
sys	0m0.537s

I got out.dat file.
I’m reading document of psi4.
I think psi4 is one of nice open source tool of computational chemistry.
I’ll study quantum chemistry.

iwatobipen$ cat out.dat 
    -----------------------------------------------------------------------
          Psi4: An Open-Source Ab Initio Electronic Structure Package
                              Psi4 1.0rc0 Driver

                          Git: Rev {_conda_cache_origin_head} 15fc63c 

    J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein,
    F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke,
    M. L. Abrams, N. J. Russ, M. L. Leininger, C. L. Janssen, E. T. Seidl,
    W. D. Allen, H. F. Schaefer, R. A. King, E. F. Valeev, C. D. Sherrill,
    and T. D. Crawford, WIREs Comput. Mol. Sci. 2, 556-565 (2012)
    (doi: 10.1002/wcms.93)

                         Additional Contributions by
    A. E. DePrince, U. Bozkaya, A. Yu. Sokolov, D. G. A. Smith, R. Di Remigio,
    R. M. Richard, J. F. Gonthier, H. R. McAlexander, M. Saitow, and
    B. P. Pritchard
    -----------------------------------------------------------------------


    Psi4 started on: Sat May 21 14:22:12 2016

    Process ID:   5371
    PSI4DATADIR: /Users/iwatobipen/.pyenv/versions/anaconda-2.4.0/share/psi4
    Memory level set to 256.000 MB

  ==> Input File <==

--------------------------------------------------------------------------
#! psi4input

molecule inputmol {
0 1

C 2.12643946569 -0.148076377601 0.280172726109

C 0.737112657203 -0.0563266799793 -0.220757232163

C 0.0894292148607 1.17309385348 -0.272662802846

C -1.26388345191 1.29962103669 0.0220149512856

C -1.99460856032 0.133837305932 0.199621363415

C -1.42391701302 -1.05115760794 -0.254293865166

C -0.0465313107692 -1.17873499271 -0.410958831529

H 2.06321195027 -0.624503035816 1.29295672814

H 2.54504673868 0.854508459836 0.429328334067

H 2.763215791 -0.825994020449 -0.31296772677

H 0.680672612474 2.01686875154 -0.629481432558

H -1.64020718362 2.2772793939 0.325940161739

H -3.02108403858 0.17949008105 0.519270208312

H -2.0538300122 -1.92491235886 -0.416048311006

H 0.438933140245 -2.12499380908 -0.552134271027

}

set basis 6-31G**

set reference uhf

energy( "scf" )

--------------------------------------------------------------------------

*** tstart() called on takayukis-MacBook-Pro.local
*** at Sat May 21 14:22:13 2016


         ---------------------------------------------------------
                                   SCF
            by Justin Turney, Rob Parrish, and Andy Simmonett
                              UHF Reference
                        6 Threads,    256 MiB Core
         ---------------------------------------------------------

  ==> Geometry <==

    Molecular point group: c1
    Full point group: C1

    Geometry (in Angstrom), charge = 0, multiplicity = 1:

       Center              X                  Y                   Z               Mass       
    ------------   -----------------  -----------------  -----------------  -----------------
           C          2.338487050599    -0.168643622045     0.358601540853    12.000000000000
           C          0.949160242112    -0.076893924424    -0.142328417419    12.000000000000
           C          0.301476799769     1.152526609036    -0.194233988102    12.000000000000
           C         -1.051835867001     1.279053792246     0.100443766030    12.000000000000
           C         -1.782560975411     0.113270061488     0.278050178159    12.000000000000
           C         -1.211869428111    -1.071724852384    -0.175865050422    12.000000000000
           C          0.165516274139    -1.199302237154    -0.332530016785    12.000000000000
           H          2.275259535179    -0.645070280260     1.371385542884     1.007825032070
           H          2.757094323589     0.833941215392     0.507757148811     1.007825032070
           H          2.975263375909    -0.846561264893    -0.234538912026     1.007825032070
           H          0.892720197383     1.996301507096    -0.551052617814     1.007825032070
           H         -1.428159598711     2.256712149456     0.404368976483     1.007825032070
           H         -2.809036453671     0.158922836606     0.597699023056     1.007825032070
           H         -1.841782427291    -1.945479603304    -0.337619496262     1.007825032070
           H          0.650980725154    -2.145561053524    -0.473705456283     1.007825032070

  Running in c1 symmetry.

  Rotational constants: A =      0.17882  B =      0.08793  C =      0.06227 [cm^-1]
  Rotational constants: A =   5360.95290  B =   2636.07055  C =   1866.71219 [MHz]
  Nuclear repulsion =  271.879087177987628

  Charge       = 0
  Multiplicity = 1
  Electrons    = 50
  Nalpha       = 25
  Nbeta        = 25

  ==> Algorithm <==

  SCF Algorithm Type is DF.
  DIIS enabled.
  MOM disabled.
  Fractional occupation disabled.
  Guess Type is CORE.
  Energy threshold   = 1.00e-06
  Density threshold  = 1.00e-06
  Integral threshold = 0.00e+00

  ==> Primary Basis <==

  Basis Set: 6-31G**
    Number of shells: 66
    Number of basis function: 145
    Number of Cartesian functions: 145
    Spherical Harmonics?: false
    Max angular momentum: 2

  ==> Pre-Iterations <==

   -------------------------------------------------------
    Irrep   Nso     Nmo     Nalpha   Nbeta   Ndocc  Nsocc
   -------------------------------------------------------
     A        145     145       0       0       0       0
   -------------------------------------------------------
    Total     145     145      25      25      25       0
   -------------------------------------------------------

  ==> Integral Setup <==

  ==> DFJK: Density-Fitted J/K Matrices <==

    J tasked:                  Yes
    K tasked:                  Yes
    wK tasked:                  No
    OpenMP threads:              6
    Integrals threads:           6
    Memory (MB):               183
    Algorithm:                Core
    Integral Cache:           NONE
    Schwarz Cutoff:          1E-12
    Fitting Condition:       1E-12

   => Auxiliary Basis Set <=

  Basis Set: 
    Number of shells: 240
    Number of basis function: 767
    Number of Cartesian functions: 767
    Spherical Harmonics?: false
    Max angular momentum: 3

  Minimum eigenvalue in the overlap matrix is 6.6381480018E-04.
  Using Symmetric Orthogonalization.
  SCF Guess: Core (One-Electron) Hamiltonian.

  ==> Iterations <==

                           Total Energy        Delta E     RMS |[F,P]|

   @DF-UHF iter   1:  -205.23304615372797   -2.05233e+02   6.04480e-02 
   @DF-UHF iter   2:  -170.13079066988712    3.51023e+01   5.20196e-02 DIIS
   @DF-UHF iter   3:  -226.19365660788827   -5.60629e+01   4.49774e-02 DIIS
   @DF-UHF iter   4:  -254.10885985388944   -2.79152e+01   3.33189e-02 DIIS
   @DF-UHF iter   5:  -268.80837968087422   -1.46995e+01   8.19862e-03 DIIS
   @DF-UHF iter   6:  -269.58619393756834   -7.77814e-01   2.85186e-03 DIIS
   @DF-UHF iter   7:  -269.69754836796079   -1.11354e-01   8.65443e-04 DIIS
   @DF-UHF iter   8:  -269.71152278937490   -1.39744e-02   3.29262e-04 DIIS
   @DF-UHF iter   9:  -269.71354390608400   -2.02112e-03   1.57613e-04 DIIS
   @DF-UHF iter  10:  -269.71419440323365   -6.50497e-04   5.79147e-05 DIIS
   @DF-UHF iter  11:  -269.71433938905835   -1.44986e-04   2.62495e-05 DIIS
   @DF-UHF iter  12:  -269.71437022304207   -3.08340e-05   1.44521e-05 DIIS
   @DF-UHF iter  13:  -269.71438177501813   -1.15520e-05   3.78687e-06 DIIS
   @DF-UHF iter  14:  -269.71438263614840   -8.61130e-07   1.01334e-06 DIIS
   @DF-UHF iter  15:  -269.71438267844809   -4.22997e-08   3.50502e-07 DIIS

  ==> Post-Iterations <==

   @Spin Contamination Metric:  -8.881784197E-14
   @S^2 Expected:                0.000000000E+00
   @S^2 Observed:               -8.881784197E-14
   @S   Expected:                0.000000000E+00
   @S   Observed:                0.000000000E+00

    Orbital Energies (a.u.)
    -----------------------

    Alpha Occupied:                                                       

       1A    -11.232193     2A    -11.228721     3A    -11.228561  
       4A    -11.228179     5A    -11.223100     6A    -11.222292  
       7A    -11.220477     8A     -1.153293     9A     -1.043927  
      10A     -1.003502    11A     -0.928172    12A     -0.818142  
      13A     -0.787929    14A     -0.687707    15A     -0.627932  
      16A     -0.626262    17A     -0.576379    18A     -0.569458  
      19A     -0.561284    20A     -0.539412    21A     -0.504515  
      22A     -0.483283    23A     -0.452811    24A     -0.321967  
      25A     -0.304862  

    Alpha Virtual:                                                        

      26A      0.137085    27A      0.148732    28A      0.243181  
      29A      0.263514    30A      0.288079    31A      0.301195  
      32A      0.316476    33A      0.318434    34A      0.334419  
      35A      0.343650    36A      0.371991    37A      0.417017  
      38A      0.479660    39A      0.491313    40A      0.525890  
      41A      0.529168    42A      0.652208    43A      0.734546  
      44A      0.741592    45A      0.750655    46A      0.768670  
      47A      0.783180    48A      0.790292    49A      0.822630  
      50A      0.834993    51A      0.838154    52A      0.845380  
      53A      0.868186    54A      0.897361    55A      0.909334  
      56A      0.931892    57A      0.969084    58A      0.998175  
      59A      1.033428    60A      1.070228    61A      1.077825  
      62A      1.088957    63A      1.095137    64A      1.116538  
      65A      1.151673    66A      1.160011    67A      1.187652  
      68A      1.205487    69A      1.224702    70A      1.237064  
      71A      1.249958    72A      1.297255    73A      1.339902  
      74A      1.360023    75A      1.378131    76A      1.443170  
      77A      1.539029    78A      1.559924    79A      1.594993  
      80A      1.646101    81A      1.722888    82A      1.724744  
      83A      1.837813    84A      1.849758    85A      1.955008  
      86A      1.987361    87A      2.078420    88A      2.090871  
      89A      2.133700    90A      2.135491    91A      2.186987  
      92A      2.226221    93A      2.267664    94A      2.282939  
      95A      2.298196    96A      2.324082    97A      2.340714  
      98A      2.372543    99A      2.374890   100A      2.391041  
     101A      2.415743   102A      2.447705   103A      2.527864  
     104A      2.594864   105A      2.616772   106A      2.664942  
     107A      2.694929   108A      2.716371   109A      2.731899  
     110A      2.741392   111A      2.769221   112A      2.802090  
     113A      2.834865   114A      2.837176   115A      2.867525  
     116A      2.889260   117A      2.909152   118A      2.932169  
     119A      2.967933   120A      3.019370   121A      3.037355  
     122A      3.087291   123A      3.128255   124A      3.178227  
     125A      3.188039   126A      3.218667   127A      3.421825  
     128A      3.429738   129A      3.475633   130A      3.551688  
     131A      3.624122   132A      3.641477   133A      3.708909  
     134A      3.778148   135A      3.819525   136A      3.828444  
     137A      3.905211   138A      4.280298   139A      4.580291  
     140A      4.598241   141A      4.617029   142A      4.801975  
     143A      4.885618   144A      4.963210   145A      5.228419  

    Beta Occupied:                                                        

       1A    -11.232193     2A    -11.228721     3A    -11.228561  
       4A    -11.228179     5A    -11.223100     6A    -11.222292  
       7A    -11.220477     8A     -1.153293     9A     -1.043927  
      10A     -1.003502    11A     -0.928172    12A     -0.818142  
      13A     -0.787929    14A     -0.687707    15A     -0.627932  
      16A     -0.626262    17A     -0.576379    18A     -0.569458  
      19A     -0.561284    20A     -0.539412    21A     -0.504515  
      22A     -0.483283    23A     -0.452811    24A     -0.321967  
      25A     -0.304862  

    Beta Virtual:                                                         

      26A      0.137085    27A      0.148732    28A      0.243181  
      29A      0.263514    30A      0.288079    31A      0.301195  
      32A      0.316476    33A      0.318434    34A      0.334419  
      35A      0.343650    36A      0.371991    37A      0.417017  
      38A      0.479660    39A      0.491313    40A      0.525890  
      41A      0.529168    42A      0.652208    43A      0.734546  
      44A      0.741592    45A      0.750655    46A      0.768670  
      47A      0.783180    48A      0.790292    49A      0.822630  
      50A      0.834993    51A      0.838154    52A      0.845380  
      53A      0.868186    54A      0.897361    55A      0.909334  
      56A      0.931892    57A      0.969084    58A      0.998175  
      59A      1.033428    60A      1.070228    61A      1.077825  
      62A      1.088957    63A      1.095137    64A      1.116538  
      65A      1.151673    66A      1.160011    67A      1.187652  
      68A      1.205487    69A      1.224702    70A      1.237064  
      71A      1.249958    72A      1.297255    73A      1.339902  
      74A      1.360023    75A      1.378131    76A      1.443170  
      77A      1.539029    78A      1.559924    79A      1.594993  
      80A      1.646101    81A      1.722888    82A      1.724744  
      83A      1.837813    84A      1.849758    85A      1.955008  
      86A      1.987361    87A      2.078420    88A      2.090871  
      89A      2.133700    90A      2.135491    91A      2.186987  
      92A      2.226221    93A      2.267664    94A      2.282939  
      95A      2.298196    96A      2.324082    97A      2.340714  
      98A      2.372543    99A      2.374890   100A      2.391041  
     101A      2.415743   102A      2.447705   103A      2.527864  
     104A      2.594864   105A      2.616772   106A      2.664942  
     107A      2.694929   108A      2.716371   109A      2.731899  
     110A      2.741392   111A      2.769221   112A      2.802090  
     113A      2.834865   114A      2.837176   115A      2.867525  
     116A      2.889260   117A      2.909152   118A      2.932169  
     119A      2.967933   120A      3.019370   121A      3.037355  
     122A      3.087291   123A      3.128255   124A      3.178227  
     125A      3.188039   126A      3.218667   127A      3.421825  
     128A      3.429738   129A      3.475633   130A      3.551688  
     131A      3.624122   132A      3.641477   133A      3.708909  
     134A      3.778148   135A      3.819525   136A      3.828444  
     137A      3.905211   138A      4.280298   139A      4.580291  
     140A      4.598241   141A      4.617029   142A      4.801975  
     143A      4.885618   144A      4.963210   145A      5.228419  

    Final Occupation by Irrep:
              A 
    DOCC [    25 ]
    SOCC [     0 ]

  Energy converged.

  @DF-UHF Final Energy:  -269.71438267844809

   => Energetics <=

    Nuclear Repulsion Energy =            271.8790871779876284
    One-Electron Energy =                -902.0381267071792308
    Two-Electron Energy =                 360.4446568507435131
    DFT Exchange-Correlation Energy =       0.0000000000000000
    Empirical Dispersion Energy =           0.0000000000000000
    PCM Polarization Energy =               0.0000000000000000
    EFP Energy =                            0.0000000000000000
    Total Energy =                       -269.7143826784481462


  Saving occupied orbitals to File 180.

  UHF NO Occupations:
  HONO-2 :   23  A 2.0000000
  HONO-1 :   24  A 2.0000000
  HONO-0 :   25  A 2.0000000
  LUNO+0 :   26  A 0.0000000
  LUNO+1 :   27  A 0.0000000
  LUNO+2 :   28  A 0.0000000
  LUNO+3 :   29  A 0.0000000


*** tstop() called on takayukis-MacBook-Pro.local at Sat May 21 14:22:17 2016
Module time:
	user time   =      13.83 seconds =       0.23 minutes
	system time =       0.40 seconds =       0.01 minutes
	total time  =          4 seconds =       0.07 minutes
Total time:
	user time   =      13.83 seconds =       0.23 minutes
	system time =       0.40 seconds =       0.01 minutes
	total time  =          4 seconds =       0.07 minutes


Properties will be evaluated at   0.000000,   0.000000,   0.000000 Bohr

Properties computed using the SCF density matrix

  Nuclear Dipole Moment: (a.u.)
     X:     3.2552      Y:    -0.3157      Z:     1.2040

  Electronic Dipole Moment: (a.u.)
     X:    -3.0898      Y:     0.2837      Z:    -1.0784

  Dipole Moment: (a.u.)
     X:     0.1654      Y:    -0.0320      Z:     0.1256     Total:     0.2101

  Dipole Moment: (Debye)
     X:     0.4204      Y:    -0.0814      Z:     0.3192     Total:     0.5341


*** Psi4 exiting successfully. Buy a developer a beer!
広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中