Basis set optimization with chemtools

In this tutorial we will go over a few examples illustrating how to use the chemtools package to optimize the exponents of orbital basis sets using various scenarios.

First some important imports:

In [1]:
from chemtools.basisset import BasisSet
from chemtools.basisopt import BSOptimizer
from chemtools.molecule import Molecule
from chemtools.calculators.molpro import Molpro

Define the program that will be used to perform the energy calculations, in this particular case it will be Molpro.

In [2]:
mp = Molpro(exevar="MOLPRO_EXE", runopts=["-s", "-n", "1", "-d", "."])

Optimization of even tempered parameters at the Hartree-Fock level for Be atom

In the first example we will optimize the s exponents for the HF calculations.

define the system for which the optimization will be performed

In [3]:
be = Molecule('Be', atoms=[('Be',)])

template for the input file to be used in sigle point calculations

In [4]:
templ =  '''***,be
memory,100,m

%geometry

%basis

gthresh,energy=1.0e-9
{rhf; wf,4,1,0}

'''
In [5]:
bso = BSOptimizer(objective='hf total energy', template=templ, code=mp, mol=be,
                  fsopt={'Be' : [('s', 'et', 8, (0.1, 2.0)),],}, verbose=False)
In [6]:
bso.run()
Script name : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages/ipykernel_launcher.py
Workdir     : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages
Start time  : 2017-11-27 14:56:58.315938
================================================================================
=============================STARTING OPTIMIZATION==============================
================================================================================

======================================CODE======================================
<Molpro(
        name=Molpro,
        molpropath=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin,
        executable=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin/molpro,
        scratch=/home/lmentel/scratch,
        runopts=['-s', '-n', '1', '-d', '.'],
)>
======================================MOL=======================================
Name: Be         Charge: 0          Multiplicty: 1          Electrons: 4
Atoms:
Element    Nuclear Charge              x              y              z
Be                   4.00               0.00000        0.00000        0.00000
=====================================OPTALG=====================================
{'jacob': None,
 'method': 'Nelder-Mead',
 'options': {'disp': True, 'maxiter': 100},
 'tol': 0.0001}Optimization terminated successfully.
         Current function value: -14.566522
         Iterations: 39
         Function evaluations: 73
 final_simplex: (array([[ 0.07025538,  3.55536302],
       [ 0.07025056,  3.55544753],
       [ 0.07025069,  3.55534628]]), array([-14.56652238, -14.56652238, -14.56652238]))
           fun: -14.566522375296
       message: 'Optimization terminated successfully.'
          nfev: 73
           nit: 39
        status: 0
       success: True
             x: array([ 0.07025538,  3.55536302])Elapsed time :               18.503 sec

We can retrieve the optimized basis directly from the optimized through:

In [8]:
bso.get_basis()
Out[8]:
{'Be': <BasisSet(
    name                 = None
    element              = None
    family               = None
    kind                 = None

 ================s shell=================
 Uncontracted:
     1      504.5070405637     1.00000000
     2      141.9002891736     1.00000000
     3       39.9116175763     1.00000000
     4       11.2257503268     1.00000000
     5        3.1574132559     1.00000000
     6        0.8880705680     1.00000000
     7        0.2497833732     1.00000000
     8        0.0702553781     1.00000000
 )>}

The raw optimization results are also available

In [9]:
bso.result
Out[9]:
final_simplex: (array([[ 0.07025538,  3.55536302],
       [ 0.07025056,  3.55544753],
       [ 0.07025069,  3.55534628]]), array([-14.56652238, -14.56652238, -14.56652238]))
           fun: -14.566522375296
       message: 'Optimization terminated successfully.'
          nfev: 73
           nit: 39
        status: 0
       success: True
             x: array([ 0.07025538,  3.55536302])

Optimization of diffuse functions for cc-pvdz for Be

Since Be atom doesn’t form a stable negative ion diffuse functions for Be are optimized for BeH-

In [10]:
beh = Molecule(name="BeH-", atoms=[('Be',), ('H', (0.0, 0.0, 2.724985))], sym="cnv 2", charge=-1, multiplicity=1)
In [11]:
bsstr = '''basis={
s,Be,2.940000E+03,4.412000E+02,1.005000E+02,2.843000E+01,9.169000E+00,3.196000E+00,1.159000E+00,1.811000E-01,5.890000E-02;
c,1.9,6.800000E-04,5.236000E-03,2.660600E-02,9.999300E-02,2.697020E-01,4.514690E-01,2.950740E-01,1.258700E-02,-3.756000E-03;
c,1.9,-1.230000E-04,-9.660000E-04,-4.831000E-03,-1.931400E-02,-5.328000E-02,-1.207230E-01,-1.334350E-01,5.307670E-01,5.801170E-01;
c,9.9,1.000000E+00;

p,Be,3.619000E+00,7.110000E-01,1.951000E-01,6.018000E-02;
c,1.4,2.911100E-02,1.693650E-01,5.134580E-01,4.793380E-01;
c,4.4,1.000000E+00;

d,Be,2.354000E-01;
c,1.1,1.000000E+00;

s, H , 13.0100000, 1.9620000, 0.4446000, 0.1220000, 0.0297400
c, 1.3, 0.0196850, 0.1379770, 0.4781480
c, 4.4, 1
c, 5.5, 1
p, H , 0.7270000, 0.1410000
c, 1.1, 1
c, 2.2, 1
}
'''
diffusetmp =  '''***,be
memory,100,m

%geometry

%basis

%core

gthresh,energy=1.0e-9
{rhf; wf,6,1,0}
cisd

'''
In [12]:
bsd = BasisSet.from_str(bsstr, fmt='molpro', name='cc-pvdz')
In [13]:
difffs = {'Be' : [('s', 'exp', 1, (0.02,)), ('p', 'exp', 1, (0.01,)), ('d', 'exp', 1, (0.07,))]}
In [14]:
bso = BSOptimizer(objective='cisd total energy', template=diffusetmp, code=mp, mol=beh,
                  fsopt=difffs, staticbs=bsd, core=[1,0,0,0,0,0,0,0], verbose=False)
In [15]:
bso.run()
Script name : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages/ipykernel_launcher.py
Workdir     : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages
Start time  : 2017-11-27 15:00:26.887764
================================================================================
=============================STARTING OPTIMIZATION==============================
================================================================================

======================================CODE======================================
<Molpro(
        name=Molpro,
        molpropath=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin,
        executable=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin/molpro,
        scratch=/home/lmentel/scratch,
        runopts=['-s', '-n', '1', '-d', '.'],
)>
======================================MOL=======================================
Name: BeH-       Charge: -1         Multiplicty: 1          Electrons: 6
Atoms:
Element    Nuclear Charge              x              y              z
Be                   4.00               0.00000        0.00000        0.00000
H                    1.00               0.00000        0.00000        2.72499
=====================================OPTALG=====================================
{'jacob': None,
 'method': 'Nelder-Mead',
 'options': {'disp': True, 'maxiter': 100},
 'tol': 0.0001}Optimization terminated successfully.
         Current function value: -15.158347
         Iterations: 71
         Function evaluations: 130
 final_simplex: (array([[-3.87366099, -1.66449926, -2.56251311],
       [-3.87366099, -1.66456835, -2.56251583],
       [-3.8736129 , -1.66454853, -2.56246558],
       [-3.87374307, -1.66455336, -2.56249475]]), array([-15.15834699, -15.15834699, -15.15834699, -15.15834699]))
           fun: -15.158346988389001
       message: 'Optimization terminated successfully.'
          nfev: 130
           nit: 71
        status: 0
       success: True
             x: array([-3.87366099, -1.66449926, -2.56251311])Elapsed time :               47.167 sec
In [16]:
bso.get_basis()
Out[16]:
{'Be': <BasisSet(
    name                 = None
    element              = None
    family               = None
    kind                 = None

 ================s shell=================
 Contracted:
     1     2940.0000000000     0.00068000    -0.00012300
     2      441.2000000000     0.00523600    -0.00096600
     3      100.5000000000     0.02660600    -0.00483100
     4       28.4300000000     0.09999300    -0.01931400
     5        9.1690000000     0.26970200    -0.05328000
     6        3.1960000000     0.45146900    -0.12072300
     7        1.1590000000     0.29507400    -0.13343500
     8        0.1811000000     0.01258700     0.53076700
     9        0.0589000000    -0.00375600     0.58011700
 Uncontracted:
    10        0.0207821468     1.00000000
    11        0.0589000000     1.00000000

 ================p shell=================
 Contracted:
     1        3.6190000000     0.02911100
     2        0.7110000000     0.16936500
     3        0.1951000000     0.51345800
     4        0.0601800000     0.47933800
 Uncontracted:
     5        0.1892854161     1.00000000
     6        0.0601800000     1.00000000

 ================d shell=================
 Uncontracted:
     1        0.0771107088     1.00000000
     2        0.2354000000     1.00000000
 )>, 'H': <BasisSet(
    name                 = cc-pvdz
    element              = H
    family               = None
    kind                 = None

 ================s shell=================
 Contracted:
     1       13.0100000000     0.01968500
     2        1.9620000000     0.13797700
     3        0.4446000000     0.47814800
 Uncontracted:
     4        0.1220000000     1.00000000
     5        0.0297400000     1.00000000

 ================p shell=================
 Uncontracted:
     1        0.7270000000     1.00000000
     2        0.1410000000     1.00000000
 )>}
In [17]:
bso.result
Out[17]:
final_simplex: (array([[-3.87366099, -1.66449926, -2.56251311],
       [-3.87366099, -1.66456835, -2.56251583],
       [-3.8736129 , -1.66454853, -2.56246558],
       [-3.87374307, -1.66455336, -2.56249475]]), array([-15.15834699, -15.15834699, -15.15834699, -15.15834699]))
           fun: -15.158346988389001
       message: 'Optimization terminated successfully.'
          nfev: 130
           nit: 71
        status: 0
       success: True
             x: array([-3.87366099, -1.66449926, -2.56251311])

Optimization of tight functions for cc-pvdz Be

In [18]:
bsd = BasisSet.from_str(bsstr, fmt='molpro', name='cc-pvdz')
pvdzbe = {k:v for k, v in bsd.items() if k=='Be'}
In [19]:
tightfs = {'Be' : [('s', 'exp', 1, (1.8,)), ('p', 'exp', 1, (4.2,))]}
In [20]:
tighttmp = '''***,be-core
memory,100,m                            !allocate 500 MW dynamic memory

%geometry

%basis

%core

{rhf; wf,4,1,0}
cisd
'''
In [21]:
bso = BSOptimizer(objective='cisd total energy', template=tighttmp, code=mp, mol=be,
                  fsopt=tightfs, staticbs=pvdzbe, runcore=True,
                  core=[[1,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0]], verbose=False)
In [22]:
bso.run()
Script name : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages/ipykernel_launcher.py
Workdir     : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages
Start time  : 2017-11-27 15:02:52.166124
================================================================================
=============================STARTING OPTIMIZATION==============================
================================================================================

======================================CODE======================================
<Molpro(
        name=Molpro,
        molpropath=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin,
        executable=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin/molpro,
        scratch=/home/lmentel/scratch,
        runopts=['-s', '-n', '1', '-d', '.'],
)>
======================================MOL=======================================
Name: Be         Charge: 0          Multiplicty: 1          Electrons: 4
Atoms:
Element    Nuclear Charge              x              y              z
Be                   4.00               0.00000        0.00000        0.00000
=====================================OPTALG=====================================
{'jacob': None,
 'method': 'Nelder-Mead',
 'options': {'disp': True, 'maxiter': 100},
 'tol': 0.0001}Optimization terminated successfully.
         Current function value: -0.031657
         Iterations: 29
         Function evaluations: 57
 final_simplex: (array([[ 0.62047435,  1.81858876],
       [ 0.62049921,  1.81859912],
       [ 0.62051376,  1.81852429]]), array([-0.0316572, -0.0316572, -0.0316572]))
           fun: -0.031657195978000985
       message: 'Optimization terminated successfully.'
          nfev: 57
           nit: 29
        status: 0
       success: True
             x: array([ 0.62047435,  1.81858876])Elapsed time :               20.818 sec
In [23]:
bso.result
Out[23]:
final_simplex: (array([[ 0.62047435,  1.81858876],
       [ 0.62049921,  1.81859912],
       [ 0.62051376,  1.81852429]]), array([-0.0316572, -0.0316572, -0.0316572]))
           fun: -0.031657195978000985
       message: 'Optimization terminated successfully.'
          nfev: 57
           nit: 29
        status: 0
       success: True
             x: array([ 0.62047435,  1.81858876])

Optimization of mid-bond function exponents

In [24]:
be2X = Molecule(name="Be2", atoms=[('Be', (0.0, 0.0, -1.5)),
                                  ('H', (0.0, 0.0,  0.0), True),
                                  ('Be', (0.0, 0.0,  1.5))], sym="dnh 2", charge=0, multiplicity=1)
In [25]:
mbfs = {'H' : [('s', 'et', 4, (0.05, 2.0)), ('p', 'et', 4, (0.04, 2.0))]}
In [26]:
mbtmp = '''***,h2o test
memory,100,m                            !allocate 500 MW dynamic memory

%geometry

%basis

dummy, H

%core

{rhf; wf,8,1,0}
cisd

'''
In [27]:
bso = BSOptimizer(objective='cisd total energy', template=mbtmp, code=mp, mol=be2X,
                  fsopt=mbfs, staticbs=pvdzbe, core=[2,0,0,0,0,0,0,0], verbose=False)
In [28]:
bso.run()
Script name : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages/ipykernel_launcher.py
Workdir     : /home/lmentel/anaconda3/envs/chemtools/lib/python3.6/site-packages
Start time  : 2017-11-27 15:03:13.366257
================================================================================
=============================STARTING OPTIMIZATION==============================
================================================================================

======================================CODE======================================
<Molpro(
        name=Molpro,
        molpropath=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin,
        executable=/home/lmentel/Programs/molprop_2012_1_Linux_x86_64_i8/bin/molpro,
        scratch=/home/lmentel/scratch,
        runopts=['-s', '-n', '1', '-d', '.'],
)>
======================================MOL=======================================
Name: Be2        Charge: 0          Multiplicty: 1          Electrons: 8
Atoms:
Element    Nuclear Charge              x              y              z
Be                   4.00               0.00000        0.00000       -1.50000
H                    0.00               0.00000        0.00000        0.00000
Be                   4.00               0.00000        0.00000        1.50000
=====================================OPTALG=====================================
{'jacob': None,
 'method': 'Nelder-Mead',
 'options': {'disp': True, 'maxiter': 100},
 'tol': 0.0001}Warning: Maximum number of iterations has been exceeded.
 final_simplex: (array([[ 0.05114879,  1.70176733,  0.05504829,  1.9460615 ],
       [ 0.05114531,  1.70019891,  0.05505673,  1.94637608],
       [ 0.05111266,  1.70488967,  0.05504958,  1.94627325],
       [ 0.05116421,  1.7018318 ,  0.05507583,  1.94578227],
       [ 0.05111091,  1.70251424,  0.05509779,  1.94542874]]), array([-29.16470765, -29.16470765, -29.16470765, -29.16470765, -29.16470765]))
           fun: -29.164707647050999
       message: 'Maximum number of iterations has been exceeded.'
          nfev: 168
           nit: 100
        status: 2
       success: False
             x: array([ 0.05114879,  1.70176733,  0.05504829,  1.9460615 ])Elapsed time :               65.526 sec
In [29]:
bso.result
Out[29]:
final_simplex: (array([[ 0.05114879,  1.70176733,  0.05504829,  1.9460615 ],
       [ 0.05114531,  1.70019891,  0.05505673,  1.94637608],
       [ 0.05111266,  1.70488967,  0.05504958,  1.94627325],
       [ 0.05116421,  1.7018318 ,  0.05507583,  1.94578227],
       [ 0.05111091,  1.70251424,  0.05509779,  1.94542874]]), array([-29.16470765, -29.16470765, -29.16470765, -29.16470765, -29.16470765]))
           fun: -29.164707647050999
       message: 'Maximum number of iterations has been exceeded.'
          nfev: 168
           nit: 100
        status: 2
       success: False
             x: array([ 0.05114879,  1.70176733,  0.05504829,  1.9460615 ])
In [30]:
%version_information chemtools
Out[30]:
SoftwareVersion
Python3.6.3 64bit [GCC 7.2.0]
IPython6.2.1
OSLinux 4.9.0 4 amd64 x86_64 with debian 9.1
chemtools0.8.4
Mon Nov 27 15:04:18 2017 CET