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:

[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.

[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

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

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

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

%geometry

%basis

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

'''
[5]:
bso = BSOptimizer(objective='hf total energy', template=templ, code=mp, mol=be,
                  fsopt={'Be' : [('s', 'et', 8, (0.1, 2.0)),],}, verbose=False)
[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:

[8]:
bso.get_basis()
[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

[9]:
bso.result
[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-

[10]:
beh = Molecule(name="BeH-", atoms=[('Be',), ('H', (0.0, 0.0, 2.724985))], sym="cnv 2", charge=-1, multiplicity=1)
[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

'''
[12]:
bsd = BasisSet.from_str(bsstr, fmt='molpro', name='cc-pvdz')
[13]:
difffs = {'Be' : [('s', 'exp', 1, (0.02,)), ('p', 'exp', 1, (0.01,)), ('d', 'exp', 1, (0.07,))]}
[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)
[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
[16]:
bso.get_basis()
[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
 )>}
[17]:
bso.result
[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

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

%geometry

%basis

%core

{rhf; wf,4,1,0}
cisd
'''
[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)
[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
[23]:
bso.result
[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

[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)
[25]:
mbfs = {'H' : [('s', 'et', 4, (0.05, 2.0)), ('p', 'et', 4, (0.04, 2.0))]}
[26]:
mbtmp = '''***,h2o test
memory,100,m                            !allocate 500 MW dynamic memory

%geometry

%basis

dummy, H

%core

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

'''
[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)
[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
[29]:
bso.result
[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 ])
[30]:
%version_information chemtools
[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