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", "."])
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])
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])
[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])
[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]:
Software | Version |
---|---|
Python | 3.6.3 64bit [GCC 7.2.0] |
IPython | 6.2.1 |
OS | Linux 4.9.0 4 amd64 x86_64 with debian 9.1 |
chemtools | 0.8.4 |
Mon Nov 27 15:04:18 2017 CET |