Draw similarity network #RDKit #Cyjupyter

Recently Kei Ono who is developer of cytoscape developed cyjupyter.
https://pypi.org/project/cyjupyter/0.2.0/
It seems attractive for me because the library can draw network diagram on jupyter notebook.
There are many network structured data in chemoinformatics. For example molecule, molecular similarity map and MMP etc… I used the library to draw similarity map of molecules today.
I am newbie of the library, so following code is very simple but there are several useful examples are provided in official repository.
At first, load modules.

import os
import numpy as np
import igraph
from py2cytoscape import util
from cyjupyter import Cytoscape
from rdkit import Chem
from rdkit.Chem import DataStructs
from rdkit.Chem import AllChem
from rdkit import RDConfig
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole

I used CDK2.sdf as a sample dataset

filedir = os.path.join(RDConfig.RDDocsDir,'Book/data/cdk2.sdf')
mols = [mol for mol in Chem.SDMolSupplier(filedir) if mol != None]
for mol in mols:
    AllChem.Compute2DCoords(mol)
fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2) for mol in mols]
smiles_list = [Chem.MolToSmiles(mol) for mol in mols]

Then make graph object, node as an each molecule and make edge if tanimoto similarity more 0.5.

g = igraph.Graph()
for smiles in smiles_list:
    g.add_vertex(name=smiles)
for i in range(len(mols)):
    for j in range(i):
        tc = DataStructs.TanimotoSimilarity(fps[i], fps[j])
        if tc >= 0.5:
            g.add_edge(smiles_list[i], smiles_list[j])

Finally convert graph object to json by using py2cytoscape and Draw graph with default settings.

graph_data = util.from_igraph(g)
Cytoscape(data=graph_data)

I could get following image.

And check network


Cyjupyter is useful network drawing tool for jupyter notebook user. I would like to check the way to control visualization.
My example code is uploaded my repository. URL is below.
https://nbviewer.jupyter.org/github/iwatobipen/chemo_info/blob/master/rdkit_notebook/Cyjupyter.ipynb

Advertisements

Calculate HOMO and LUMO with Psi4 reviced #RDKit #Psi4

Yesterday, I got comments from reader.
Regarding the comment, to calculate HOMO LUMO with psi4 correct way is below.

import psi4
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
psi4.core.set_output_file("output1.dat", True)
def mol2xyz(mol):
    mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol, useExpTorsionAnglePrefs=True,useBasicKnowledge=True)
    AllChem.UFFOptimizeMolecule(mol)
    atoms = mol.GetAtoms()
    string = string = "\n"
    for i, atom in enumerate(atoms):
        pos = mol.GetConformer().GetAtomPosition(atom.GetIdx())
        string += "{} {} {} {}\n".format(atom.GetSymbol(), pos.x, pos.y, pos.z)
    string += "units angstrom\n"
    return string, mol

Next, calculate HOMO-LUMO of benzene with the function and psi4.

mol = Chem.MolFromSmiles("c1ccccc1")
xyz, mol=mol2xyz(mol)
psi4.set_memory('4 GB')
psi4.set_num_threads(4)
benz = psi4.geometry(xyz)
%time scf_e, scf_wfn = psi4.energy("B3LYP/cc-pVDZ", return_wfn=True)
> CPU times: user 13.8 s, sys: 226 ms, total: 14 s
> Wall time: 4.51 s

After the calculation, I could access HOMO-LUMO, the code is below.

# HOMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha()]
# LUMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha() + 1]
HOMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha()-1]
LUMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha()]
print(HOMO, LUMO, scf_e)
> -0.2529021800443842 -0.006506519238935586 -232.26252757817264

Check log file

!cat out2.dat

  Memory set to   3.725 GiB by Python driver.
  Threads set to 4 by Python driver.

*** tstart() called on takayukis-MacBook-Pro.local
*** at Mon Aug 27 22:26:56 2018

   => Loading Basis Set <=

    Name: CC-PVDZ
    Role: ORBITAL
    Keyword: BASIS
    atoms 1-6  entry C          line   138 file /Users/iwatobipen/.pyenv/versions/anaconda3-4.2.0/share/psi4/basis/cc-pvdz.gbs 
    atoms 7-12 entry H          line    22 file /Users/iwatobipen/.pyenv/versions/anaconda3-4.2.0/share/psi4/basis/cc-pvdz.gbs 

    There are an even number of electrons - assuming singlet.
    Specify the multiplicity in the molecule input block.


         ---------------------------------------------------------
                                   SCF
            by Justin Turney, Rob Parrish, Andy Simmonett
                             and Daniel Smith
                              RKS Reference
                        4 Threads,   3814 MiB Core
         ---------------------------------------------------------

  ==> Geometry <==

    Molecular point group: c1
    Full point group: C1

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

       Center              X                  Y                   Z               Mass       
    ------------   -----------------  -----------------  -----------------  -----------------
         C            1.160791903590    -0.776928934439    -0.079150803890    12.000000000000
         C            1.255739388827     0.616792238521    -0.002681541758    12.000000000000
         C            0.094947479782     1.393721178119     0.076469246393    12.000000000000
         C           -1.160791914824     0.776928940041     0.079150771790    12.000000000000
         C           -1.255739387065    -0.616792248432     0.002681608472    12.000000000000
         C           -0.094947472815    -1.393721173263    -0.076469222838    12.000000000000
         H            2.058813166717    -1.377983059255    -0.140384154914     1.007825032070
         H            2.227214722117     1.093960079368    -0.004756245235     1.007825032070
         H            0.168401517636     2.471943131887     0.135627614553     1.007825032070
         H           -2.058813208566     1.377983069413     0.140383740077     1.007825032070
         H           -2.227214704360    -1.093960120700     0.004756292428     1.007825032070
         H           -0.168401463717    -2.471943107238    -0.135627939511     1.007825032070

  Running in c1 symmetry.

  Rotational constants: A =      0.18924  B =      0.18924  C =      0.09462 [cm^-1]
  Rotational constants: A =   5673.32397  B =   5673.32388  C =   2836.66196 [MHz]
  Nuclear repulsion =  203.019334728971273

  Charge       = 0
  Multiplicity = 1
  Electrons    = 42
  Nalpha       = 21
  Nbeta        = 21

  ==> Algorithm <==

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

  ==> Primary Basis <==

  Basis Set: CC-PVDZ
    Blend: CC-PVDZ
    Number of shells: 54
    Number of basis function: 114
    Number of Cartesian functions: 120
    Spherical Harmonics?: true
    Max angular momentum: 2

  ==> DFT Potential <==

   => Composite Functional: B3LYP <= 

    B3LYP Hyb-GGA Exchange-Correlation Functional

    P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Phys. Chem. 98, 11623 (1994)

    Deriv               =              1
    GGA                 =           TRUE
    Meta                =          FALSE

    Exchange Hybrid     =           TRUE
    MP2 Hybrid          =          FALSE

   => Exchange Functionals <=

    0.0800   Slater exchange
    0.7200         Becke 88

   => Exact (HF) Exchange <=

    0.2000               HF 

   => Correlation Functionals <=

    0.1900   Vosko, Wilk & Nusair (VWN5_RPA)
    0.8100   Lee, Yang & Parr

   => Molecular Quadrature <=

    Radial Scheme       =       TREUTLER
    Pruning Scheme      =           FLAT
    Nuclear Scheme      =       TREUTLER

    BS radius alpha     =              1
    Pruning alpha       =              1
    Radial Points       =             75
    Spherical Points    =            302
    Total Points        =         266172
    Total Blocks        =           2084
    Max Points          =            255
    Max Functions       =            114

   => Loading Basis Set <=

    Name: (CC-PVDZ AUX)
    Role: JKFIT
    Keyword: DF_BASIS_SCF
    atoms 1-6  entry C          line   121 file /Users/iwatobipen/.pyenv/versions/anaconda3-4.2.0/share/psi4/basis/cc-pvdz-jkfit.gbs 
    atoms 7-12 entry H          line    51 file /Users/iwatobipen/.pyenv/versions/anaconda3-4.2.0/share/psi4/basis/cc-pvdz-jkfit.gbs 

  ==> Pre-Iterations <==

   -------------------------------------------------------
    Irrep   Nso     Nmo     Nalpha   Nbeta   Ndocc  Nsocc
   -------------------------------------------------------
     A        114     114       0       0       0       0
   -------------------------------------------------------
    Total     114     114      21      21      21       0
   -------------------------------------------------------

  ==> Integral Setup <==

  DFHelper Memory: AOs need 0.070 [GiB]; user supplied 2.794 [GiB]. Using in-core AOs.

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

    J tasked:                   Yes
    K tasked:                   Yes
    wK tasked:                   No
    OpenMP threads:               4
    Memory (MB):               2861
    Algorithm:                 Core
    Schwarz Cutoff:           1E-12
    Mask sparsity (%):       0.3693
    Fitting Condition:        1E-12

   => Auxiliary Basis Set <=

  Basis Set: (CC-PVDZ AUX)
    Blend: CC-PVDZ-JKFIT
    Number of shells: 198
    Number of basis function: 558
    Number of Cartesian functions: 636
    Spherical Harmonics?: true
    Max angular momentum: 3

  Minimum eigenvalue in the overlap matrix is 3.7184237071E-04.
  Using Symmetric Orthogonalization.

  SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF.

  ==> Iterations <==

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

   @DF-RKS iter   0:  -232.99504225236745   -2.32995e+02   7.10044e-02 
   @DF-RKS iter   1:  -232.10322360969730    8.91819e-01   8.92096e-03 
   @DF-RKS iter   2:  -232.08890690303650    1.43167e-02   9.91157e-03 DIIS
   @DF-RKS iter   3:  -232.26148696855563   -1.72580e-01   7.72966e-04 DIIS
   @DF-RKS iter   4:  -232.26244386501853   -9.56896e-04   1.55832e-04 DIIS
   @DF-RKS iter   5:  -232.26249567324737   -5.18082e-05   1.40053e-04 DIIS
   @DF-RKS iter   6:  -232.26252736183829   -3.16886e-05   1.12370e-05 DIIS
   @DF-RKS iter   7:  -232.26252757817264   -2.16334e-07   5.75991e-07 DIIS

  ==> Post-Iterations <==

    Orbital Energies [Eh]
    ---------------------

    Doubly Occupied:                                                      

       1A    -10.190567     2A    -10.190357     3A    -10.190356  
       4A    -10.189878     5A    -10.189877     6A    -10.189659  
       7A     -0.851906     8A     -0.745429     9A     -0.745428  
      10A     -0.602584    11A     -0.602584    12A     -0.521871  
      13A     -0.463024    14A     -0.444273    15A     -0.421160  
      16A     -0.421159    17A     -0.365359    18A     -0.344262  
      19A     -0.344262    20A     -0.252905    21A     -0.252902  

    Virtual:                                                              

      22A     -0.006507    23A     -0.006506    24A      0.059666  
      25A      0.099451    26A      0.099452    27A      0.133761  
      28A      0.133764    29A      0.150228    30A      0.153828  
      31A      0.277077    32A      0.277080    33A      0.295653  
      34A      0.295656    35A      0.405118    36A      0.408148  
      37A      0.464413    38A      0.465367    39A      0.512461  
      40A      0.512464    41A      0.517424    42A      0.517430  
      43A      0.519856    44A      0.519868    45A      0.526820  
      46A      0.540745    47A      0.592924    48A      0.592924  
      49A      0.647786    50A      0.647787    51A      0.650076  
      52A      0.650083    53A      0.667968    54A      0.720592  
      55A      0.788201    56A      0.828354    57A      0.857690  
      58A      0.883159    59A      0.883161    60A      0.924558  
      61A      0.992120    62A      0.992122    63A      1.005964  
      64A      1.005973    65A      1.011555    66A      1.011587  
      67A      1.059132    68A      1.071624    69A      1.071625  
      70A      1.216812    71A      1.280232    72A      1.280245  
      73A      1.456675    74A      1.476446    75A      1.476447  
      76A      1.510782    77A      1.538630    78A      1.598096  
      79A      1.598099    80A      1.606414    81A      1.606425  
      82A      1.633651    83A      1.682795    84A      1.682801  
      85A      1.692381    86A      1.692388    87A      1.725702  
      88A      1.795113    89A      1.795115    90A      1.839139  
      91A      1.839159    92A      1.844378    93A      1.848385  
      94A      1.848386    95A      1.883194    96A      1.981556  
      97A      1.981559    98A      1.989993    99A      1.990001  
     100A      2.021139   101A      2.191006   102A      2.261460  
     103A      2.371980   104A      2.426824   105A      2.426829  
     106A      2.435274   107A      2.435279   108A      2.641049  
     109A      2.641053   110A      2.695891   111A      2.795913  
     112A      2.914328   113A      2.914351   114A      3.646462  

    Final Occupation by Irrep:
              A 
    DOCC [    21 ]

  Energy converged.

  @DF-RKS Final Energy:  -232.26252757817264

   => Energetics <=

    Nuclear Repulsion Energy =            203.0193347289712733
    One-Electron Energy =                -713.2777384427657807
    Two-Electron Energy =                 306.1488686089047064
    DFT Exchange-Correlation Energy =     -28.1529924732828292
    Empirical Dispersion Energy =           0.0000000000000000
    VV10 Nonlocal Energy =                  0.0000000000000000
    Total Energy =                       -232.2625275781726373



Properties will be evaluated at   0.000000,   0.000000,   0.000000 [a0]

Properties computed using the SCF density matrix

  Nuclear Dipole Moment: [e a0]
     X:     0.0000      Y:    -0.0000      Z:    -0.0000

  Electronic Dipole Moment: [e a0]
     X:     0.0000      Y:    -0.0000      Z:     0.0000

  Dipole Moment: [e a0]
     X:     0.0000      Y:    -0.0000      Z:    -0.0000     Total:     0.0000

  Dipole Moment: [D]
     X:     0.0000      Y:    -0.0000      Z:    -0.0000     Total:     0.0000


*** tstop() called on takayukis-MacBook-Pro.local at Mon Aug 27 22:27:00 2018
Module time:
	user time   =      13.71 seconds =       0.23 minutes
	system time =       0.22 seconds =       0.00 minutes
	total time  =          4 seconds =       0.07 minutes
Total time:
	user time   =      27.58 seconds =       0.46 minutes
	system time =       0.47 seconds =       0.01 minutes
	total time  =         32 seconds =       0.53 minutes

That’s all. I am happy because I can get many response through the my blog post and I can have many opportunity to learn many things.

Calculate HOMO and LUMO with Psi4 #RDKit #Psi4

You know Psi4 is an open-source suite of ab initio quantum chemistry programs designed for efficient, high-accuracy simulations of a variety of molecular properties. It is very easy to use and has an optional Python interface.
It is useful for us I think. Because Psi4 can use in python, it means we can integrate many libraries in python!  And it is worth to know that, to communicate numpy and psi4 is very easy.

Today, I conducted HOMO LUMO calculation with Psi4 and RDKit.

HOMO-LUMO gap is often used to estimate risk of drug-induced phototoxicity.

Let's start! ;-)

At first, import libraries and define the mol2xyz function. The function generates a conformation and converts RDKit mol object to xyz format.

psi4.core.set_output_file(“output1.dat”, True) is used to logging.

import psi4
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
psi4.core.set_output_file("output1.dat", True)
def mol2xyz(mol):
    mol = Chem.AddHs(mol)
    AllChem.EmbedMolecule(mol, useExpTorsionAnglePrefs=True,useBasicKnowledge=True)
    AllChem.UFFOptimizeMolecule(mol)
    atoms = mol.GetAtoms()
    string = "\n"
    for i, atom in enumerate(atoms):
        pos = mol.GetConformer().GetAtomPosition(atom.GetIdx())
        string += "{} {} {} {}\n".format(atom.GetSymbol(), pos.x, pos.y, pos.z)
    string += "units angstrom\n"
    return string, mol

Next, calculate HOMO-LUMO of benzene with the function and psi4.

mol = Chem.MolFromSmiles("c1ccccc1")
xyz, mol=mol2xyz(mol)
psi4.set_memory('4 GB')
psi4.set_num_threads(4)
benz = psi4.geometry(xyz)
%time scf_e, scf_wfn = psi4.energy("B3LYP/cc-pVDZ", return_wfn=True)
>CPU times: user 13.5 s, sys: 207 ms, total: 13.7 s
>Wall time: 4.46 s

I set return_wfn argument is True because I want to get wave function information.
Energy calculation need many CPU time. To perform many compound calculation, I need more machine power!

After the calculation, I could access HOMO-LUMO, the code is below.

HOMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha()]
LUMO = scf_wfn.epsilon_a_subset('AO', 'ALL').np[scf_wfn.nalpha() + 1]
print(HOMO, LUMO, scf_e)
>-0.006507529999155065 -0.006506586740442874 -232.26253075556204

Yah, easy… But the value is quite different from some literatures. Hmm Am I wrong ? Maybe yes…
Any advice and suggestions will be greatly appreciated.

It is very easy to check log.
I am reading API and reference of Psi4 now!

!cat out.dat

 Memory set to   3.725 GiB by Python driver.
  Threads set to 4 by Python driver.

  Memory set to   3.725 GiB by Python driver.
  Threads set to 4 by Python driver.

*** tstart() called on ********
*** at Fri Aug 24 22:21:33 2018

   => Loading Basis Set <=

    Name: CC-PVDZ
    Role: ORBITAL
    Keyword: BASIS
---------snip------------

    There are an even number of electrons - assuming singlet.
    Specify the multiplicity in the molecule input block.


         ---------------------------------------------------------
                                   SCF
            by Justin Turney, Rob Parrish, Andy Simmonett
                             and Daniel Smith
                              RKS Reference
                        4 Threads,   3814 MiB Core
         ---------------------------------------------------------

  ==> Geometry <==

    Molecular point group: c1
    Full point group: C1

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

       Center              X                  Y                   Z               Mass       
    ------------   -----------------  -----------------  -----------------  -----------------
         C            0.806497703012    -1.143092245995     0.014915282772    12.000000000000
         C            1.393280339340     0.126835003153    -0.002136692164    12.000000000000
         C            0.586782120312     1.269928194219    -0.017051393298    12.000000000000
         C           -0.806497629817     1.143092609376    -0.014913861108    12.000000000000
         C           -1.393280087652    -0.126835526704     0.002137132311    12.000000000000
         C           -0.586782336814    -1.269927969465     0.017053614836    12.000000000000
         H            1.430428517155    -2.027419681898     0.026439293008     1.007825032070
         H            2.471161063731     0.224955697903    -0.003795091078     1.007825032070
         H            1.040731718362     2.252381167101    -0.030251698606     1.007825032070
         H           -1.430425978428     2.027420563468    -0.026458153907     1.007825032070
         H           -2.471160894666    -0.224959372018     0.003781521713     1.007825032070
         H           -1.040735716617    -2.252379143548     0.030235509138     1.007825032070

  Running in c1 symmetry.

  Rotational constants: A =      0.18924  B =      0.18924  C =      0.09462 [cm^-1]
  Rotational constants: A =   5673.32470  B =   5673.32274  C =   2836.66186 [MHz]
  Nuclear repulsion =  203.019333245138824

  Charge       = 0
  Multiplicity = 1
  Electrons    = 42
  Nalpha       = 21
  Nbeta        = 21

  ==> Algorithm <==

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

  ==> Primary Basis <==

  Basis Set: CC-PVDZ
    Blend: CC-PVDZ
    Number of shells: 54
    Number of basis function: 114
    Number of Cartesian functions: 120
    Spherical Harmonics?: true
    Max angular momentum: 2

  ==> DFT Potential <==

   => Composite Functional: B3LYP <= 

    B3LYP Hyb-GGA Exchange-Correlation Functional

    P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Phys. Chem. 98, 11623 (1994)

    Deriv               =              1
    GGA                 =           TRUE
    Meta                =          FALSE

    Exchange Hybrid     =           TRUE
    MP2 Hybrid          =          FALSE

   => Exchange Functionals <=

    0.0800   Slater exchange
    0.7200         Becke 88

   => Exact (HF) Exchange <=

    0.2000               HF 

   => Correlation Functionals <=

    0.1900   Vosko, Wilk & Nusair (VWN5_RPA)
    0.8100   Lee, Yang & Parr

   => Molecular Quadrature <=

    Radial Scheme       =       TREUTLER
    Pruning Scheme      =           FLAT
    Nuclear Scheme      =       TREUTLER

    BS radius alpha     =              1
    Pruning alpha       =              1
    Radial Points       =             75
    Spherical Points    =            302
    Total Points        =         266204
    Total Blocks        =           2062
    Max Points          =            256
    Max Functions       =            114

   => Loading Basis Set <=

    Name: (CC-PVDZ AUX)
    Role: JKFIT
    Keyword: DF_BASIS_SCF
-----snip--------
  ==> Pre-Iterations <==

   -------------------------------------------------------
    Irrep   Nso     Nmo     Nalpha   Nbeta   Ndocc  Nsocc
   -------------------------------------------------------
     A        114     114       0       0       0       0
   -------------------------------------------------------
    Total     114     114      21      21      21       0
   -------------------------------------------------------

  ==> Integral Setup <==

  DFHelper Memory: AOs need 0.070 [GiB]; user supplied 2.794 [GiB]. Using in-core AOs.

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

    J tasked:                   Yes
    K tasked:                   Yes
    wK tasked:                   No
    OpenMP threads:               4
    Memory (MB):               2861
    Algorithm:                 Core
    Schwarz Cutoff:           1E-12
    Mask sparsity (%):       0.3693
    Fitting Condition:        1E-12

   => Auxiliary Basis Set <=

  Basis Set: (CC-PVDZ AUX)
    Blend: CC-PVDZ-JKFIT
    Number of shells: 198
    Number of basis function: 558
    Number of Cartesian functions: 636
    Spherical Harmonics?: true
    Max angular momentum: 3

  Minimum eigenvalue in the overlap matrix is 3.7184246400E-04.
  Using Symmetric Orthogonalization.

  SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF.

  ==> Iterations <==

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

   @DF-RKS iter   0:  -232.98648910745848   -2.32986e+02   7.09616e-02 
   @DF-RKS iter   1:  -232.10265029904261    8.83839e-01   8.93885e-03 
   @DF-RKS iter   2:  -232.08812748497797    1.45228e-02   9.76151e-03 DIIS
   @DF-RKS iter   3:  -232.26148987002321   -1.73362e-01   7.73194e-04 DIIS
   @DF-RKS iter   4:  -232.26244676165251   -9.56892e-04   2.20816e-04 DIIS
   @DF-RKS iter   5:  -232.26249884952190   -5.20879e-05   1.40071e-04 DIIS
   @DF-RKS iter   6:  -232.26253054068945   -3.16912e-05   1.12245e-05 DIIS
   @DF-RKS iter   7:  -232.26253075556204   -2.14873e-07   5.66817e-07 DIIS

  ==> Post-Iterations <==

    Orbital Energies [Eh]
    ---------------------

    Doubly Occupied:                                                      

       1A    -10.190570     2A    -10.190360     3A    -10.190358  
       4A    -10.189881     5A    -10.189879     6A    -10.189662  
       7A     -0.851912     8A     -0.745431     9A     -0.745426  
      10A     -0.602585    11A     -0.602585    12A     -0.521876  
      13A     -0.463024    14A     -0.444275    15A     -0.421161  
      16A     -0.421155    17A     -0.365362    18A     -0.344264  
      19A     -0.344262    20A     -0.252907    21A     -0.252901  

    Virtual:                                                              

      22A     -0.006508    23A     -0.006507    24A      0.059658  
      25A      0.099450    26A      0.099455    27A      0.133761  
      28A      0.133763    29A      0.150228    30A      0.153828  
      31A      0.277076    32A      0.277078    33A      0.295654  
      34A      0.295658    35A      0.405112    36A      0.408146  
      37A      0.464412    38A      0.465356    39A      0.512461  
      40A      0.512461    41A      0.517423    42A      0.517433  
      43A      0.519852    44A      0.519872    45A      0.526750  
      46A      0.540745    47A      0.592922    48A      0.592924  
      49A      0.647783    50A      0.647787    51A      0.650079  
      52A      0.650080    53A      0.667968    54A      0.720591  
      55A      0.788173    56A      0.828353    57A      0.857688  
      58A      0.883160    59A      0.883161    60A      0.924558  
      61A      0.992121    62A      0.992123    63A      1.005962  
      64A      1.005979    65A      1.011551    66A      1.011617  
      67A      1.059133    68A      1.071623    69A      1.071625  
      70A      1.216802    71A      1.280227    72A      1.280249  
      73A      1.456674    74A      1.476443    75A      1.476447  
      76A      1.510781    77A      1.538598    78A      1.598094  
      79A      1.598098    80A      1.606419    81A      1.606454  
      82A      1.633651    83A      1.682793    84A      1.682805  
      85A      1.692380    86A      1.692389    87A      1.725701  
      88A      1.795112    89A      1.795115    90A      1.839139  
      91A      1.839180    92A      1.844376    93A      1.848382  
      94A      1.848384    95A      1.883194    96A      1.981552  
      97A      1.981560    98A      1.989989    99A      1.990002  
     100A      2.021138   101A      2.191006   102A      2.261459  
     103A      2.371976   104A      2.426825   105A      2.426828  
     106A      2.435277   107A      2.435291   108A      2.641052  
     109A      2.641053   110A      2.695891   111A      2.795915  
     112A      2.914318   113A      2.914347   114A      3.646453  

    Final Occupation by Irrep:
              A 
    DOCC [    21 ]

  Energy converged.

  @DF-RKS Final Energy:  -232.26253075556204

   => Energetics <=

    Nuclear Repulsion Energy =            203.0193332451388244
    One-Electron Energy =                -713.2776809379226961
    Two-Electron Energy =                 306.1488028199291875
    DFT Exchange-Correlation Energy =     -28.1529858827073518
    Empirical Dispersion Energy =           0.0000000000000000
    VV10 Nonlocal Energy =                  0.0000000000000000
    Total Energy =                       -232.2625307555620680



Properties will be evaluated at   0.000000,   0.000000,   0.000000 [a0]

Properties computed using the SCF density matrix

  Nuclear Dipole Moment: [e a0]
     X:    -0.0000      Y:    -0.0000      Z:    -0.0000

  Electronic Dipole Moment: [e a0]
     X:     0.0000      Y:    -0.0000      Z:     0.0000

  Dipole Moment: [e a0]
     X:    -0.0000      Y:    -0.0000      Z:    -0.0000     Total:     0.0000

  Dipole Moment: [D]
     X:    -0.0000      Y:    -0.0000      Z:    -0.0000     Total:     0.0000


*** tstop() ********l at Fri Aug 24 22:21:37 2018
Module time:
	user time   =      13.48 seconds =       0.22 minutes
	system time =       0.21 seconds =       0.00 minutes
	total time  =          4 seconds =       0.07 minutes
Total time:
	user time   =      13.48 seconds =       0.22 minutes
	system time =       0.21 seconds =       0.00 minutes
	total time  =          4 seconds =       0.07 minutes

​

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 tkyks-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 tkyks-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!