Calculate the hyperpolarizability for two water molecules the fast and simple way

The tutorial https://github.com/fishstamp82/moltools/blob/docs/share/tutorial/two_water_example.ipynb
demonstrates the step-by-step process in how the Applequist equations are solved from quantum mechanical properties, all the information is given explicitly to remove doubts in how the procedure looks.

This tutorial will perform all the above steps in a more compact, and faster way, using the moltools package.

Dalton

In this tutorial, we will be running Dalton from inside python. For this we will utilize the full path to the Dalton runscript. The location will be '/home/user/repos/dalton/build_gnu/dalton', but is of course custom to every user.

Dalton can be obtained here

Moltools

Clone the entire repository, which will come with loprop and pd automatically:

$> git clone --recursive https://github.com/fishstamp82/moltools.git

(Optional, but highly recommended) If you want optimized parallel excecution, make sure you have Cython and build the cython extensions of vahtras/pd by running the following script:

$> moltools/scripts/build_cython.sh

See the requirements.txt file for packages needed

Make sure that the path to moltools is in your PYTHONPATH env variable, and run an ipython console

$> ipython



In [1]:

    
DALTON = '/home/user/repos/dalton/build_gnu/dalton'



In [2]:

    
from moltools import Water, Molecule, Cluster

The get_standard method creates a Water instance, with the geometry of a TIP3P water model. It's oxygen is placed in origo, the molecule in the xz-plane, with it' electronic dipole pointing in the positive z-axis



In [3]:

    
w = Water.get_standard()

Take a look at it



In [4]:

    
w.plot()

the molecular property of this water molecule is accesible via Water.Property or Water.p



In [5]:

    
print w.p









    



{'beta': array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]), 'alpha': array([ 0.,  0.,  0.,  0.,  0.,  0.]), 'charge': 0.0, 'dipole': array([ 0.,  0.,  0.]), 'quadrupole': array([ 0.,  0.,  0.,  0.,  0.,  0.])}

All properties are zero. Let's obtain properties via QM using Dalton.



In [6]:

    
w.props_from_qm( method = 'hfqua', #hf is Hartree-Fock, qua is Quadratic response (LoProp Beta)
                dalpath = DALTON, 
                tmpdir = '/tmp',
               )

Now the properties are non-zero



In [7]:

    
print "Charge of the water molecule:" , w.p.q
print "Charge of oxygen:" , w.o.p.q

print "Dipole moment of the water molecule:" , w.p.d
print "Dipole moment of oxygen:" , w.o.p.d









    



Charge of the water molecule: 1.00000000058e-07
Charge of oxygen: -0.6638767
Dipole moment of the water molecule: [ 0.          0.          0.85412685]
Dipole moment of oxygen: [ 0.         0.         0.3450972]

Now it is straightforward to calculate the properties of individual molecules.

We now create an additional water molecule, translate it to (x, y, z) = (0, 0, 5 ), and calculate the Applequist beta of two water molecules interacting.



In [8]:

    
water2 = Water.get_standard().t(0, 0, 5)
water2.props_from_qm( method = 'hfqua',
                dalpath = DALTON, 
                tmpdir = '/tmp',
                )

Using Cluster, we group the molecules to produce a list of Applequist point dipoles.



In [9]:

    
c = Cluster( w, water2 )
c.populate_bonds()



In [10]:

    
#Uncomment to visualize it, you can imagine how it would look like
#c.plot()



In [11]:

    
#get_pdlist generates the input format file for PointDipoleList automatically based on the molecules
#in the cluster and their unit of coordinate

pdlist = c.get_pdlist()

print "Beta zzz no Applequist:", c.p.b[9]
print "Beta zzz after Applequist:", pdlist.beta( cython = 1, num_threads = 4 )[2,2,2]









    



Beta zzz no Applequist: -5.2818056
Beta zzz after Applequist: -5.30956890099

Example creating a molecule, can be directly from reading a .xyz file or specifying atom by atom



In [12]:

    
#By default reading in coordinates in angstrom
o2 = Molecule.from_xyz_string( """2\n\nO 0 0 0\nO 0 0 1.3""")



In [13]:

    
o2.to_AU()
o2.plot()



In [14]:

    
o2.props_from_qm( method = 'hfqua',
                 dalpath = DALTON,
                 tmpdir = '/tmp',
                )



In [15]:

    
print o2.p.a
print np.allclose( np.zeros(10), o2.p.b, atol = 1e-4 )









    



[  5.51446360e+00  -1.67396000e-02   0.00000000e+00   6.08339280e+00
   0.00000000e+00   1.73507914e+01]
True



In [ ]:



In [ ]: