This tutorial will describe the steps needed to calculate the classical polarizability and hyperpolarizability of two water molecules interacting, using the localized properties of each molecule obtained seperately from quantum mechanics.
The two interacting water molecules will be of the TIP3P model, separated at 5 bohr in the z-axis direction.
The DALTON source can be obtained from http://www.daltonprogram.org/. It uses an LGPL license.
To run a DALTON calculation once it is installed, execute the dalton runscript:
>>> dalton -get "AOPROPER AOONEINT" input.dal molecule.mol
The get argument saves the crucial binary files holding the one-electron integrals, and atomic properties.
Upon success, DALTON will produce two output files, one of which is input_molecule.tar.gz.
To see a list of options such as MPI core usage and temporary directory setup, execute:
>>> dalton
**DALTON INPUT
.RUN RESPONSE
.DIRECT
.PARALLELL
**WAVE FUNCTION
.HF
**RESPONSE
.PROPAV
XDIPLEN
.PROPAV
YDIPLEN
.PROPAV
ZDIPLEN
*QUADRATIC
.QLOP
.DIPLEN
**END OF DALTON INPUT
ATOMBASIS
Comment 1
Comment 2
Atomtypes=2 Charge=0 Nosymm
Charge=8.0 Atoms=1 Basis=ano-1 4 3 1 0
O 0.00000 0.00000 0.00000
Charge=1.0 Atoms=2 Basis=ano-1 2 0 0 0
H 1.43043 0.00000 1.10716
H -1.43043 0.00000 1.10716
Get the source code to vahtras/loprop.git:
>>> git clone --recursive https://github.com/vahtras/loprop.git
With the loprop/loprop.py script in your PATH, execute
>>> loprop.py -l 1 -a 22 -B 1 -f input_molecule.tar.gz
The loprop script outputs in left-to-right (ignoring the first digit) order:
1) The coordinates
2) The multipole moments. Charges are always included, i.e. "-l 0". Up to quadrupoles are supported at this time.
3) Polarizability. "-a 1" is isotropic. "-a 2" or "-a 22" is anisotropic xx, xy, xz, yy, yz, zz, components.
4) Hyperpolarizability. "-B 1" is the upper triangular form components xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz,
Running:
>>> loprop.py -l 1 -a 22 -B 1 -f input_molecule.tar.gz
should give the following output:
AU
3 1 22 1
1 0.000 0.000 0.000 -0.664 0.000 0.000 0.345 3.783 0.000 -0.000 3.966 -0.000 3.527 0.000 0.000 -2.984 -0.000 0.000 -0.000 0.000 1.267 0.000 2.167
1 1.430 0.000 1.107 0.332 -0.161 0.000 -0.113 1.552 -0.000 1.155 0.609 -0.000 1.211 -4.468 0.000 -4.559 -0.056 0.000 -3.720 0.000 0.460 0.000 -2.404
1 -1.430 0.000 1.107 0.332 0.161 0.000 -0.113 1.552 0.000 -1.155 0.609 -0.000 1.211 4.468 0.000 -4.559 0.056 -0.000 3.720 0.000 0.460 0.000 -2.404
Time used in Loprop : 0.53 (cpu) 0.14 (wall) ~
The Applequist equations are solved using the source code obtained from vahtras/pd.git:
>>> git clone https://github.com/vahtras/pd.git
The pd/particles.py contains the following PointDipoleList class, which is used to obtain the polarizability or hyperpolarizability
class PointDipoleList( list ):
...
It is used with the output from LoProp to generate interacting particles, but care has to be taken to modify each particles group number (The first digit from the LoProp output for each atom in the molecule). Particles with the same group number will thus not interact with each other.
Note that only the first two lines specifying the format, and the lines containing atoms are used for the strings.
Since we have two water molecules we modify the first second line:
3 1 22 1to
6 1 22 1We then duplicate the 3 succeeding lines, while also changing their group ID to 2.
Furthermore, for the second water molecule, we add 5.0 to the z-coordiante (4th column) of the 3 atoms in the molecule.mol.
Their properties are individually identical.
The resulting string will be defined as "water" in the runscript, and fully model 2 water molecules classically, where each water on itself can not interact within itself, but only with the other.
#!/usr/bin/env python
from particles import PointDipoleList as pdl
waters = """AU
6 1 22 1
1 0.000 0.000 0.000 -0.664 0.000 0.000 0.345 3.783 0.000 -0.000 3.966 -0.000 3.527 0.000 0.000 -2.984 -0.000 0.000 -0.000 0.000 1.267 0.000 2.167
1 1.430 0.000 1.107 0.332 -0.161 0.000 -0.113 1.552 -0.000 1.155 0.609 -0.000 1.211 -4.468 0.000 -4.559 -0.056 0.000 -3.720 0.000 0.460 0.000 -2.404
1 -1.430 0.000 1.107 0.332 0.161 0.000 -0.113 1.552 0.000 -1.155 0.609 -0.000 1.211 4.468 0.000 -4.559 0.056 -0.000 3.720 0.000 0.460 0.000 -2.404
2 0.000 0.000 5.000 -0.664 0.000 0.000 0.345 3.783 0.000 -0.000 3.966 -0.000 3.527 0.000 0.000 -2.984 -0.000 0.000 -0.000 0.000 1.267 0.000 2.167
2 1.430 0.000 6.107 0.332 -0.161 0.000 -0.113 1.552 -0.000 1.155 0.609 -0.000 1.211 -4.468 0.000 -4.559 -0.056 0.000 -3.720 0.000 0.460 0.000 -2.404
2 -1.430 0.000 6.107 0.332 0.161 0.000 -0.113 1.552 0.000 -1.155 0.609 -0.000 1.211 4.468 0.000 -4.559 0.056 -0.000 3.720 0.000 0.460 0.000 -2.404"""
pdlist = pdl.from_string( waters )
print "alpha:"
print pdlist.alpha()
print "beta:"
print pdlist.beta()
Assuming that particles is in the $PYTHONPATH variable, or we are in the same vahtras/pd directory, running:
python runscript.py
outputs:
[[ 1.28427223e+01 0.00000000e+00 -2.22044605e-16]
[ 0.00000000e+00 9.99993520e+00 0.00000000e+00]
[ -2.22044605e-16 0.00000000e+00 1.28192101e+01]]
[[[ 0.00000000e+00 0.00000000e+00 -2.35412114e+01]
[ 0.00000000e+00 6.93889390e-18 0.00000000e+00]
[ -2.35412114e+01 0.00000000e+00 8.88178420e-16]]
[[ 0.00000000e+00 -1.38777878e-17 0.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 4.37513260e+00]
[ 0.00000000e+00 4.37513260e+00 0.00000000e+00]]
[[ -2.35412114e+01 0.00000000e+00 4.44089210e-16]
[ 0.00000000e+00 4.37513260e+00 0.00000000e+00]
[ 4.44089210e-16 0.00000000e+00 -5.31056471e+00]]]
See tutorial https://github.com/fishstamp82/moltools/blob/docs/share/tutorial/moltools_tutorial.ipynb
In [ ]:
In [ ]: