ChemWiki - User contributions [en]

Mod:Hunt Research Group: ES alpha

2017-09-05T14:25:12Z

Ech17:

==introduction==
*we are using g09
*apparently Dalton can also calculate excited state polarizabilities http://www.daltonprogram.org

* useful [https://joaquinbarroso.com/2012/07/27/polarizability-and-hyperpolarizability-in-gaussian/ blog]
*data below relates to information in this paper [https://pubs.acs.org/doi/abs/10.1021/jp5057623 link] we have used molecule VI from Scheme1

==process==
*optimise the ground state use int=ultrafine and scf=(conver=9)
*make sure to "turn on" relevant symmetry
*confirm with a frequency analysis (using the optimised geometry)
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) guess=read int=ultrafine freq
</pre>
:we turn on extra printing with #p
SCF Done: E(RB3LYP) = -391.474123816
:you are looking for something like this
<pre>
Electronic spatial extent (au): <R**2>= 593.6542
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= 0.0000 Y= 0.0000 Z= 2.7398 Tot= 2.7398
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -35.4903 YY= -50.1334 ZZ= -37.0641
XY= 0.0000 XZ= 0.0000 YZ= 0.0000
</pre>
:the dipole moment is in the standard orientation, and the zero is therefore the CoM
:then
<pre>
Exact polarizability: 19.474 0.000 65.813 0.000 0.000 36.738
Approx polarizability: 29.762 0.000 136.916 0.000 0.000 55.771
</pre>
:the polarizability is in lower-triangular format with xx, xy, yy, xz, yz and zz values
:the approximate value is a cruder estimate evaluated using sum-over-states perturbation theory
:the polarizability is the mean(traceX)=1/3(xx+yy+zz) ie one third of the sum of the diagonal terms
:if you do this 1/3(19.474+65.813+36.738)=40.67
:note that this tensor is given relative to the standard orientation

*then you can use the checkpoint file of the freq job and run a job with polar
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar guess=read int=ultrafine
</pre>
:you will see something like this
:were clearly the tensor is symmetric and 1=y, 2=z and 3=x
<pre>
SCF Polarizability for W= 0.000000:
1 2 3
1 0.658125D+02
2 -0.119999D-03 0.367380D+02
3 0.227376D-02 0.472986D-03 0.194741D+02
Isotropic polarizability for W= 0.000000 40.67 Bohr**3.
SCF Static Hyperpolarizability:
K= 1 block:
1
1 0.000000D+00
K= 2 block:
1 2
1 0.000000D+00
2 0.000000D+00 0.000000D+00
K= 3 block:
1 2 3
1 -0.304281D+01
2 0.000000D+00 -0.718146D+02
3 0.000000D+00 0.000000D+00 0.451401D+02
</pre>
:and a bit later this
<pre>
Electronic spatial extent (au): <R**2>= 593.6543
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= 0.0000 Y= -2.7398 Z= -0.0010 Tot= 2.7398
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -50.1333 YY= -37.0642 ZZ= -35.4903
XY= 0.0000 XZ= 0.0010 YZ= 0.0013
</pre>
:and a bit later this
<pre>
Exact polarizability: 19.474 0.000 65.813 0.000 0.000 36.738
Approx polarizability: 29.762 0.000 136.916 0.000 0.000 55.771
</pre>

*or you can do both!
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) guess=read int=ultrafine polar freq
</pre>
:you will see this
<pre>
Isotropic polarizability for W= 0.000000 40.67 Bohr**3.
</pre>
:and later this
<pre>
Exact polarizability: 19.474 0.000 65.813 0.000 0.000 36.738
Approx polarizability: 29.762 0.000 136.916 0.000 0.000 55.771
</pre>
: and even later you will see this, and the dipole at least doesn't match that given above!
:this happens just before the frequencies are given
<pre>
(Enter /Applications/g09/l716.exe)
Dipole = 3.54385899D-17-1.22258974D-14 1.07791236D+00
Polarizability= 1.94740869D+01 2.40399932D-13 6.58126199D+01
2.26775166D-13 4.95235677D-08 3.67380378D+01
</pre>

*now optimise on the first excited state, use the time-dependent options, the root=1 is the first excited state
<pre>
# opt td(root=1) b3lyp/6-31g(d,p) geom=connectivity scf=(conver=9) guess=read int=ultrafine
</pre>
:you should see something like this in the output
<pre>
Excitation energies and oscillator strengths:

Excited State 1: Singlet-A 1.6713 eV 741.84 nm f=0.0003 <S**2>=0.000
25 -> 26 0.70603
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -391.233828069
</pre>

* carry out a frequency analysis on the optimised excited state geometry
<pre>
# freq b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) td(read,root=1) guess=read int=ultrafine
</pre>
:this will print out excited state dipole moments and so on ...
<pre>
Electronic spatial extent (au): <R**2>= 592.2331
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= -0.0001 Y= -2.0915 Z= 0.0004 Tot= 2.0915
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -49.3534 YY= -35.5628 ZZ= -36.5632
XY= 0.0002 XZ= 0.0008 YZ= 0.0009
</pre>

*use the polar keyword on the excited state molecule
: use just polar to get the default analytic derivatives, in this case you might want to specify freq keyword as well
: you want to use the excited state density
: in this case polar=numerical computes the numerical derivative of the dipole moment
: you can change the step size of the electric field 0.0001N au with step=N
: you can restart a numerical calculation used polar=(numerical,restart)
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar=numerical td(read,root=1)
density=current guess=read int=ultrafine
</pre>
:you are looking for something like this
<pre>
Excited State 1: Singlet-A 1.6129 eV 768.69 nm f=0.0007 <S**2>=0.000
25 -> 26 -0.70872
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -391.407446678
</pre>
:and
<pre>
Electronic spatial extent (au): <R**2>= 592.2331
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= -0.0001 Y= -2.0915 Z= 0.0158 Tot= 2.0915
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -49.3534 YY= -35.5628 ZZ= -36.5632
XY= 0.0002 XZ= 0.0002 YZ= 0.0008
</pre>
:and
<pre>
Isotropic polarizability= 49.78 Bohr**3.
1 2 3
1 0.796564D+02
2 -0.303870D-03 0.508599D+02
3 0.349968D-03 0.169881D-03 0.188218D+02
</pre>

*of course we can also do it analytically, the default
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar=numerical td(read,root=1)
density=current guess=read int=ultrafine
</pre>
:and you should see
<pre>
Isotropic polarizability= 49.78 Bohr**3.
1 2 3
1 0.796564D+02
2 -0.292909D-03 0.508599D+02
3 0.360142D-03 0.138905D-03 0.188218D+02
</pre>
<pre>
(Enter /Applications/g09/l716.exe)
Dipole =-4.88131552D-05-8.22840998D-01-4.04220282D-05
Polarizability= 7.96563849D+01-2.92908823D-04 5.08598561D+01
3.60141967D-04 1.38905076D-04 1.88217717D+01
</pre>

*so what are the results?
:for our system the
::ground state dipole 2.74 D and isotropic polarizability 41 Bohr**3
::excited state dipole 2.09 D and isotropic polarizability 50 Bohr**3.
:according to the paper (M06-2X with Sadlej's basis set
::ground state dipole 1.02 D and polarizability 50 Bohr**3
::excited state dipole 0.66 D and polarizability 57 Bohr**3.
:so these differences are quite large!
:they could be due to
::basis set
::functional
::larger finite field
::convergence on energy
::convergence on the geometry

==Reproduction of the exact results==
Let us know try to reproduce the exact same results as described by Jacquemin and coworkers by applying the same procedure. First, optimise on the M06-2X/6-31G(d) level of theory using an improved energy threshold, a tight geometry optimisation criterion and a high level DFT integration grid.
<pre>
#p m062x/6-31g(d) opt=tight int=ultrafine scf=(conver=10)
</pre>
The optimised geometry is
<pre>
N -0.142626 0.000000 0.028260
N -0.072667 0.000000 1.258841
C 1.363634 0.000000 1.690124
N 2.090317 0.000000 0.518596
C 1.235531 0.000000 -0.563026
O 1.468028 0.000000 -1.729522
H 3.098469 0.000000 0.461269
O 1.726826 0.000000 2.822758
</pre>
with an energy of -391.3255482 Ha. Check that the geometry is a true minimum on the potential energy surface by conducting a frequency calculation
<pre>
#p m062x/6-31g(d) geom=checkpoint int=ultrafine scf=(conver=10)
guess=read freq
</pre>
(check for NImag=0).
The dipole moment and polarizability was calculated using the M06-2X functional and Sadlej's pVTZ basis set obtained from the EMSL Bais Set Exchange on the optimized geometry via
<pre>
#p m062x/gen int=ultrafine scf=(conver=10) guess=read polar
geom=checkpoint

H 0
S 4 1.00
33.8650140 0.0060680
5.0947880 0.0453160
1.1587860 0.2028460
0.3258400 0.5037090
S 1 1.00
0.1027410 1.0000000
S 1 1.00
0.0324000 1.0000000
P 2 1.00
1.1588000 0.1884400
0.3258000 0.8824200
P 2 1.00
0.1027000 0.1178000
0.0324000 0.0042000
****
C 0
S 5 1.00
5240.6353000 0.0009370
782.2048000 0.0072280
178.3508300 0.0363440
50.8159420 0.1306000
16.8235620 0.3189310
S 2 1.00
6.1757760 0.4387420
2.4180490 0.2149740
S 1 1.00
0.5119000 1.0000000
S 1 1.00
0.1565900 1.0000000
S 1 1.00
0.0479000 1.0000000
P 4 1.00
18.8418000 0.0138870
4.1592400 0.0862790
1.2067100 0.2887440
0.3855400 0.4994110
P 1 1.00
0.1219400 1.0000000
P 1 1.00
0.0385680 1.0000000
D 2 1.00
1.2067000 0.2628500
0.3855000 0.8043000
D 2 1.00
0.1219000 0.6535000
0.0386000 0.8636000
****
N 0
S 5 1.00
8104.0716000 0.0008020
1216.0215000 0.0061740
277.2342800 0.0312330
76.9040230 0.1151980
25.8744190 0.2969510
S 2 1.00
9.3467670 0.4473490
3.5797940 0.2450030
S 1 1.00
0.7396100 1.0000000
S 1 1.00
0.2226170 1.0000000
S 1 1.00
0.0670060 1.0000000
P 4 1.00
26.8689870 0.0144780
5.9912270 0.0911560
1.7508420 0.2974200
0.5605110 0.4937960
P 1 1.00
0.1759480 1.0000000
P 1 1.00
0.0552310 1.0000000
D 2 1.00
1.7508000 0.2247700
0.5605000 0.6595600
D 2 1.00
0.1795900 0.8713600
0.0552000 0.7042200
****
O 0
S 5 1.00
10662.2850000 0.0007990
1599.7097000 0.0061530
364.7252600 0.0311570
103.6517900 0.1155960
33.9058050 0.3015520
S 2 1.00
12.2874690 0.4448700
4.7568050 0.2431720
S 1 1.00
1.0042710 1.0000000
S 1 1.00
0.3006860 1.0000000
S 1 1.00
0.0900300 1.0000000
P 4 1.00
34.8564630 0.0156480
7.8431310 0.0981970
2.3062490 0.3077680
0.7231640 0.4924700
P 1 1.00
0.2148820 1.0000000
P 1 1.00
0.0638500 1.0000000
D 2 1.00
2.3062000 0.2027000
0.7232000 0.5791000
D 2 1.00
0.2149000 0.7854500
0.0639000 0.5338700
****
</pre>
yielding an energy of -391.4434344 Ha, a dipole moment of 1.02 au. and a polarizability of 50 au. (trace of the polarizability tensor)
<pre>
Dipole=1.0190643,0.,-0.0579513
Polar=45.7964203,0.,30.0056759,1.6195023,0.,74.1929892
</pre>
which corresponds exactly to the ground state gas phase results listed in the reference. Now optimise the excited state and check for imaginary frequencies via
<pre>
#p m062x/6-31g(d) opt=tight td(root=1) int=ultrafine scf=(conver=10)
guess=read geom=connectivity freq
</pre>
where the optimised geometry is
<pre>
N -0.099617 0.000000 0.023379
N -0.029402 0.000000 1.258821
C 1.333978 0.000000 1.722300
N 2.044096 0.000000 0.521225
C 1.202401 0.000000 -0.591640
O 1.500232 0.000000 -1.747732
H 3.054940 0.000000 0.463750
O 1.760884 0.000000 2.837197
</pre>
Dipole moment and polarizabilities can then again be calculated at the M06-2X/Sadlej level of theory
<pre>
#p m062x/gen td(root=1) int=ultrafine scf=(conver=10) guess=read
polar=numerical density=current geom=checkpoint
</pre>
yielding
<pre>
Dipole=0.6557277,0.,-0.03725
Polar=54.0435904,0.,29.7590519,1.8414148,0.,86.3255147
</pre>
which is a total dipole moment of 0.66 au. and a polarizability of 57 au., which corresponds perfectly to the reference. An alternative approach to calculate the excited state polarizability is via the use of Stark's relation, which links the energy of a state to the strength of an applied electric field.
<math display="block"> E(F)=E^0-\sum_a \mu_a F_a - \frac{1}{2}\sum_{a,b} \alpha_{ab} F_a F_b </math>
The dipole moment is then the first derivative of the energy with respect to the field at zero field strength. Make sure to take into account all directions, as <math display="inline">\mu_x=\big( \frac{\mathrm{d}E}{\mathrm{d}F_x} \big) _{F_x=0} </math>,
<math display="inline">\mu_y=\big( \frac{\mathrm{d}E}{\mathrm{d}F_y} \big) _{F_y=0} </math> and <math display="inline">\mu_z=\big( \frac{\mathrm{d}E}{\mathrm{d}F_z} \big) _{F_z=0} </math> and the total dipole moment is then <math display="inline">\mu=\sqrt{\mu_x^2+\mu_y^2+\mu_z^2}</math>. Analogously, the polarizability is obtained via <math display="inline">\alpha_{xx}=\big( \frac{\mathrm{d}^2E}{\mathrm{d}F_x^2} \big) _{F_x=0} </math>, <math display="inline">\alpha_{yy}=\big( \frac{\mathrm{d}^2E}{\mathrm{d}F_y^2} \big) _{F_y=0} </math>, <math display="inline">\alpha_{zz}=\big( \frac{\mathrm{d}^2E}{\mathrm{d}F_z^2} \big) _{F_z=0} </math> and <math display="inline">\alpha=\frac{\alpha_{xx}+\alpha_{yy}+\alpha_{zz}}{3}</math>. The electric fields are chosen such that the Romberg differentiation procedure [http://dx.doi.org/10.1016/j.theochem.2007.08.028 link] can be applied.
<math display="block"> F_a=2^k \cdot 0.0004 \; \mathrm{au.} \;\;\;\;\;\; \mathrm{with} \;k=1,2,...5 \;\;\; \mathrm{for} \; a=x,y,z</math>
To apply for example a field of 0.0008 au. in the x-direction, use
<pre>
#p m062x/gen td(root=1) int=ultrafine scf=(conver=10)
guess=read density=current field=x+8 geometry=checkpoint
</pre>
Search the output for the energy of the first excited state
<pre>
Excited State 1: Singlet-A'' 1.6794 eV 738.27 nm f=0.0008
<S**2>=0.000
25 -> 26 0.70425
25 <- 26 -0.10327
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -391.373283614
</pre>
Here, we get <math display="inline">E(F_x=0.0008)=</math>-391.373283614 Ha. Repeat the calculation for all field strengths, which finally yields
<pre>
F [au.] E(Fx=F) [au.] E(Fy=F) [au.] E(Fz=F) [au.]
-0.0128 -391.380340688 -391.386106671 -391.375696126
-0.0064 -391.375024593 -391.378566574 -391.373865503
-0.0032 -391.373697934 -391.375634101 -391.373408310
-0.0016 -391.373366469 -391.374375865 -391.373294027
-0.0008 -391.373283571 -391.373798627 -391.373265458
0.0000 -391.373255935 -391.373255935 -391.373255935
0.0008 -391.373283614 -391.372747763 -391.373265458
0.0016 -391.373366555 -391.372274081 -391.373294027
0.0032 -391.373698106 -391.371430018 -391.373408310
0.0064 -391.375024938 -391.370154501 -391.373865503
0.0128 -391.380341380 -391.369251683 -391.375696126
</pre>
The Romberg differentiation then uses
<math display="block"> \big( \frac{\mathrm{d}^2f}{\mathrm{d}x^2} \big)_{x=0}= \frac{f(-2^kh)+f(2^kh)-2f(0)}{(2^kh)^2} </math>
where h is the minimum step size (here 0.0004 au.) to evaluate the second derivative and then uses the iterative formula
<math display="block"> P^{p,k} = \frac{4^pP^{p-1,k}-P^{p-1,k+1}}{4^p-1} </math>
where p is the number of Romberg iteration and P is the required differentiation, to improve the obtained value. The value at p=0 corresponds to the initially obtained uncorrected differentiation value. The iterative procedure
<pre>
k 2^k h d^2E/dFx^2 p=1 p=2 p=3 p=4
1 0.0008 -86.4 -86.4 -86.4 -86.4 -86.4
2 0.0016 -86.4 -86.4 -86.4 -86.4
3 0.0032 -86.3 -86.3 -86.3
4 0.0064 -86.4 -86.3
5 0.0128 -86.5
</pre>
then yields <math display="inline">\alpha_{xx}</math>= 86 au. Analogue evaluation in y and z direction yields <math display="inline">\alpha_{yy}</math>= 54 au and <math display="inline">\alpha_{zz}</math>= 30 au. The polarizability is then
57 au, which is the same as the value in Reference, as well as the value obtained using the default routine for the
calculation of polarizability in Gaussian09.

Finally, try to reproduce the correct dipole moment and polarisability in implicit solvent (dichloromethane). Optimise the ground state via
<pre>
#p m062x/6-31g(d) opt=tight int=ultrafine scf=(conver=10)
SCRF=(IEFPCM,SOLVENT=dichloromethane) guess=read
</pre>
yielding the geometry
<pre>
N -0.141332 0.000000 0.028556
N -0.071415 0.000000 1.258401
C 1.367139 0.000000 1.684979
N 2.090981 0.000000 0.518558
C 1.239597 0.000000 -0.558312
O 1.461088 0.000000 -1.729021
H 3.101578 0.000000 0.461093
O 1.719875 0.000000 2.823046
</pre>
at an energy of -391.335724 Ha. Dipole moments and polarizabilities at the optimized geometry at the M06-2X/Sadlej level of theory were then obtained via
<pre>
#p m062x/gen int=ultrafine scf=(conver=10) guess=read polar
SCRF=(IEFPCM,SOLVENT=dichloromethane) geom=checkpoint
</pre>
to yield the following output:
<pre>
Dipole=1.2615935,0.,-0.0717451
Polar=55.4963369,0.,35.741816,2.3407747,0.,96.5379704
</pre>
The dipole moment in groundstate is therefore 1.26 au. and the
polarizability 63 au.

The excited state in implicit solvent was also optimized:
<pre>
#p m062x/6-31g(d) opt=tight td(root=1) int=ultrafine
scf=(conver=10) guess=read geom=(connectivity,checkpoint)
SCRF=(IEFPCM,SOLVENT=dichloromethane)
</pre>
using the optimized excited state geometry in gas phase as input. The optimized geometry is then
<pre>
N -0.101348 0.000000 0.023438
N -0.031138 0.000000 1.258959
C 1.337822 0.000000 1.716250
N 2.046763 0.000000 0.521074
C 1.206893 0.000000 -0.586065
O 1.494126 0.000000 -1.746788
H 3.059686 0.000000 0.463479
O 1.754707 0.000000 2.836953
</pre>
Calculation of the dipole moment and the polarizability is then done on the M06-2X/Sadlej level of theory using the default calculation routine.
<pre>
#p m062x/gen td(root=1) int=ultrafine scf=(conver=10)
guess=read polar=numerical density=current geom=checkpoint
SCRF=(IEFPCM,SOLVENT=dichloromethane)
</pre>
The dipole moment is then 0.84 au. and the polarizability 73 au. which is for both ground and excited state in perfect agreement with the published results in the paper by Jacquemin and coworkers.
<pre>
Dipole=0.8376841,0.,-0.047563
Polar=67.8846862,0.,35.5510398,2.7715247,0.,116.4748693
</pre>

==Atomic polarizabilities==
We now want to calculate the atomic contributions to the molecular polarizability. Again, we use Stark's relation, namely that the change in dipole moment with respect to an external field is the polarizability. Changes in dipoles arise from atomic dipoles (asymmetric distribution of charge around the nucleus) that change upon applying an external field and from charge transfer between the atoms. The GDMA program written by Anthony Stone can calculate atomic dipoles
and partial charges at specific sites in the molecule, here the nuclei. We reuse the calculations at external fields of 0.0008au in negative and positive x, y and z direction from above. The respective checkpoint files need to be converted to formatted checkpoint files. From these the atomic dipoles can be calculated using the following input script for the ground state:
<pre>
Title "title"
File file.fchk

Angstrom
Multipoles
Limit 1
Start

Finish
</pre>
For the excited state, the electron density needs to be specified by putting the keyword "DENSITY "CI" at the end of the "File" line. This input script tells GDMA to calculate multipoles up to the rank of dipoles at all nuclear sites. The output (here for Fx=0.0008au) looks like:
<pre>
[...]
Positions and radii in angstrom
Multipole moments in atomic units, ea_0^k for rank k

N x = -0.616284 y = -1.257223 z = 0.000000 angstrom
Maximum rank = 1 Radius = 0.650 angstrom
Q00 = -0.056538
|Q1| = 0.511629 Q11c = 0.331633 Q11s = 0.389594

N x = 0.616284 y = -1.257222 z = 0.000000 angstrom
Maximum rank = 1 Radius = 0.650 angstrom
Q00 = -0.060705
|Q1| = 0.516025 Q11c = -0.335883 Q11s = 0.391746
[...]
x = 2.279816 y = 0.450604 z = 0.000000 angstrom
Maximum rank = 1 Radius = 0.650 angstrom
Q00 = -0.280276
|Q1| = 0.425044 Q11c = -0.405995 Q11s = -0.125819

Total multipoles referred to origin at
x = 0.000000, y = 0.000000, z = 0.000000 angstrom
Q00 = -0.000002
|Q1| = 1.022461 Q11c = -0.059386 Q11s = 1.020735
</pre>
where [...] indicates that not all data is shown. For the first nitrogen atom, a atomic dipole (0.331633,0.389594,0) is observed (Q11c is the x direction, Q11s the y direction and Q10 the z direction). For the calculation of atomic polarizability from atomic dipoles we need to monitor the x-component of the dipole when applying a field in x-direction, the y component for fields in y-direction and so on. For the nitrogen atom we get:
<pre>
Fx= 0.0008au , mu_x= 0.331633
Fx=-0.0008au , mu_x= 0.335880
Fy= 0.0008au , mu_y= 0.387470
Fy=-0.0008au , mu_y= 0.393892
Fz= 0.0008au , mu_z=-0.003752
Fz=-0.0008au , mu_z= 0.003752
</pre>
Numerical differentiation yields <math display="inline">\alpha_{xx}=-\big( \frac{\mathrm{d}\mu_x}{\mathrm{d}F_x} \big) _{F_x=0} = -\frac{\mu_x(F_x)-\mu_x(-F_x)}{2F_x} = 2.65 au^3</math>, and analogously
<math display="inline">\alpha_{yy}=4.01 au^3</math> and <math display="inline">\alpha_{zz}=4.69 au^3</math> yielding an average polarizability of 3.8au^3. The same is done for all other atoms and then repeated for the excited state. One then ends up with the atomic polarizabilities
<pre>
Atomic alpha from polarization:
GS ES
N 3.8 4.1
N 3.8 4.1
C 1.2 1.4
N 2.5 2.3
C 1.2 1.4
O 4.6 4.7
H 0.8 0.8
O 4.6 4.7
</pre>
Summed up, we end up with a polarizability from atomic dipoles of 22.4 au^3 in the ground state and 23.4 au^3 in the excited state. The missing 27.5 au^3 in the ground and 33.3 au^3 in the excited state stem from charge transfer. The overall contribution of charge transfer over all atoms is easy to calculate via <math display="inline">\mu_{charge transfer}=\sum_ir_iq_i</math> where <math display="inline">r_i</math> are the coordinates of each atom i and <math display="inline">q_i</math> the charge at that site. The charge at a specific site can be found in the GDMA output script and is labelled "Q00". The coordinates are also printed in the output script. Summing up these values gives the desired polarizability of 27.5 au^3 in the ground and 33.3 au^3 in the excited state. However, the separation of this contribution to different atomic sites is difficult, as the sum <math display="inline">\mu_ct=\sum_ir_iq_i</math> is origin independent (for a neutral molecule) but the respective contributions <math display="inline">r_iq_i</math> aren't. A possible solution to this problem is the definition of bond charges <math display="inline">b_{i,j}</math> between the bonded atoms i and j. Each partial charge is assumed to be the sum of all bond charges, the sum of all bond charges in a ring is zero and same bond charges in different directions add up to 0 (<math display="inline">b_{i,j}+b_{j,i}=0</math>). These relations can be used to set up linear independent equations, which are solved for <math display="inline">b_{i,j}</math>.

For our molecule, we need to find a set of non-unique bond charges, for example:

N(1)-N(2) => <math display="inline">b_{1,2}</math> 
N(1)-C(5) => <math display="inline">b_{1,5}</math> 
N(2)-C(3) => <math display="inline">b_{2,3}</math> 
C(3)-N(4) => <math display="inline">b_{3,4}</math> 
C(3)-O(8) => <math display="inline">b_{3,8}</math> 
N(4)-C(5) => <math display="inline">b_{4,5}</math> 
N(4)-H(7) => <math display="inline">b_{4,7}</math> 
C(5)-O(6) => <math display="inline">b_{5,6}</math> 

Thus, the following linear equations evolve:

q(N1)= <math display="inline">b_{1,2}-b_{1,5}</math> 
q(N2)= -<math display="inline">b_{1,2}+b_{2,3}</math> 
q(C3)= -<math display="inline">b_{2,3}+b_{3,4}+b_{3,8}</math> 
q(N4)= -<math display="inline">b_{3,4}+b_{4,5}+b_{4,7}</math> 
q(C5)= -<math display="inline">b_{4,5}-b_{1,5}+b_{5,6}</math> 
q(O6)= -<math display="inline">b_{5,6}</math> 
q(H7)= -<math display="inline">b_{4,7}</math> 
q(O8)= -<math display="inline">b_{3,8}</math> 
<math display="inline">0 = b_{1,2}+b_{2,3}+b_{3,4}+b_{4,5}-b_{1,5}</math> (ring condition) 

This set of equations (9 equations for 8 variables) is solved using numpy (the charges are known, see GDMA output), to yield all <math display="inline">b_{i,j}</math>. We get for the field in positive x-direction

<math display="inline">b_{1,2}</math>= 0.0036 
<math display="inline">b_{1,5}</math>= -0.0602 
<math display="inline">b_{2,3}</math>= -0.0571 
<math display="inline">b_{3,4}</math>= -0.0538 
<math display="inline">b_{3,8}</math>= 0.2803 
<math display="inline">b_{4,5}</math>= 0.0471 
<math display="inline">b_{4,7}</math>= -0.2531 
<math display="inline">b_{5,6}</math>= 0.2709 

The atomic contribution to the dipole moment arising from charge transfer is then <math display="inline">\mu_{charge transfer, i}=\sum_j(r_i-\frac{r_i-r_{i,j}}{2})*b_{i,j}</math> where <math display="inline">r_{i,j}</math> was simply assumed to be <math display="inline">\frac{r_i+r_j}{2}</math>, thus we get

<math display="inline">\mu_{charge transfer, 1}</math> = [ -0.033, 0.080, 0.0000] in au 
<math display="inline">\mu_{charge transfer, 2}</math> = [ 0.023, 0.076, 0.0000] in au 
... 

for the dipole moments when applying a field of 0.0008au in x-direction.

The same procedure as explained above for the atomic polarization contribution is applied here: the atomic charge transfer dipole moments are calculated at different fields and differentiated numerically to yield the average polarizabilities from charge transfer
<pre>
Atomic alpha from charge transfer
GS ES
N 3.3 4.1
N 3.3 4.1
C 5.7 6.9
N 4.5 5.1
C 5.7 6.9
O 2.2 2.7
H 0.7 0.6
O 2.2 2.7
</pre>
To obtain atomic polarizabilities we now sum up the contributions from polarization and charge transfer
<pre>
Overall atomic alpha
GS ES
N 7.1 8.2
N 7.1 8.2
C 6.8 8.4
N 7.0 7.5
C 6.8 8.4
O 6.9 7.4
H 1.5 1.4
O 6.9 7.4
</pre>
In this molecule, the atomic polarization increases slightly upon excitation (in average about 6%), and charge transfer polarizability to a larger extent (about 19%). The calculation of atomic polarizabilities (polarization and charge transfer term) from the GDMA output files was automatized during this project and the script will be available soon at http://www.mdy.univie.ac.at

Mod:Hunt Research Group: ES alpha

2017-03-10T09:51:52Z

Ech17:

==introduction==
*we are using g09
*apparently Dalton can also calculate excited state polarizabilities http://www.daltonprogram.org

* useful [https://joaquinbarroso.com/2012/07/27/polarizability-and-hyperpolarizability-in-gaussian/ blog]
*data below relates to information in this paper [https://pubs.acs.org/doi/abs/10.1021/jp5057623 link] we have used molecule VI from Scheme1

==process==
*optimise the ground state use int=ultrafine and scf=(conver=9)
*make sure to "turn on" relevant symmetry
*confirm with a frequency analysis (using the optimised geometry)
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) guess=read int=ultrafine freq
</pre>
:we turn on extra printing with #p
SCF Done: E(RB3LYP) = -391.474123816
:you are looking for something like this
<pre>
Electronic spatial extent (au): <R**2>= 593.6542
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= 0.0000 Y= 0.0000 Z= 2.7398 Tot= 2.7398
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -35.4903 YY= -50.1334 ZZ= -37.0641
XY= 0.0000 XZ= 0.0000 YZ= 0.0000
</pre>
:the dipole moment is in the standard orientation, and the zero is therefore the CoM
:then
<pre>
Exact polarizability: 19.474 0.000 65.813 0.000 0.000 36.738
Approx polarizability: 29.762 0.000 136.916 0.000 0.000 55.771
</pre>
:the polarizability is in lower-triangular format with xx, xy, yy, xz, yz and zz values
:the approximate value is a cruder estimate evaluated using sum-over-states perturbation theory
:the polarizability is the mean(traceX)=1/3(xx+yy+zz) ie one third of the sum of the diagonal terms
:if you do this 1/3(19.474+65.813+36.738)=40.67
:note that this tensor is given relative to the standard orientation

*then you can use the checkpoint file of the freq job and run a job with polar
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar guess=read int=ultrafine
</pre>
:you will see something like this
:were clearly the tensor is symmetric and 1=y, 2=z and 3=x
<pre>
SCF Polarizability for W= 0.000000:
1 2 3
1 0.658125D+02
2 -0.119999D-03 0.367380D+02
3 0.227376D-02 0.472986D-03 0.194741D+02
Isotropic polarizability for W= 0.000000 40.67 Bohr**3.
SCF Static Hyperpolarizability:
K= 1 block:
1
1 0.000000D+00
K= 2 block:
1 2
1 0.000000D+00
2 0.000000D+00 0.000000D+00
K= 3 block:
1 2 3
1 -0.304281D+01
2 0.000000D+00 -0.718146D+02
3 0.000000D+00 0.000000D+00 0.451401D+02
</pre>
:and a bit later this
<pre>
Electronic spatial extent (au): <R**2>= 593.6543
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= 0.0000 Y= -2.7398 Z= -0.0010 Tot= 2.7398
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -50.1333 YY= -37.0642 ZZ= -35.4903
XY= 0.0000 XZ= 0.0010 YZ= 0.0013
</pre>
:and a bit later this
<pre>
Exact polarizability: 19.474 0.000 65.813 0.000 0.000 36.738
Approx polarizability: 29.762 0.000 136.916 0.000 0.000 55.771
</pre>

*or you can do both!
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) guess=read int=ultrafine polar freq
</pre>
:you will see this
<pre>
Isotropic polarizability for W= 0.000000 40.67 Bohr**3.
</pre>
:and later this
<pre>
Exact polarizability: 19.474 0.000 65.813 0.000 0.000 36.738
Approx polarizability: 29.762 0.000 136.916 0.000 0.000 55.771
</pre>
: and even later you will see this, and the dipole at least doesn't match that given above!
:this happens just before the frequencies are given
<pre>
(Enter /Applications/g09/l716.exe)
Dipole = 3.54385899D-17-1.22258974D-14 1.07791236D+00
Polarizability= 1.94740869D+01 2.40399932D-13 6.58126199D+01
2.26775166D-13 4.95235677D-08 3.67380378D+01
</pre>

*now optimise on the first excited state, use the time-dependent options, the root=1 is the first excited state
<pre>
# opt td(root=1) b3lyp/6-31g(d,p) geom=connectivity scf=(conver=9) guess=read int=ultrafine
</pre>
:you should see something like this in the output
<pre>
Excitation energies and oscillator strengths:

Excited State 1: Singlet-A 1.6713 eV 741.84 nm f=0.0003 <S**2>=0.000
25 -> 26 0.70603
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -391.233828069
</pre>

* carry out a frequency analysis on the optimised excited state geometry
<pre>
# freq b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) td(read,root=1) guess=read int=ultrafine
</pre>
:this will print out excited state dipole moments and so on ...
<pre>
Electronic spatial extent (au): <R**2>= 592.2331
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= -0.0001 Y= -2.0915 Z= 0.0004 Tot= 2.0915
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -49.3534 YY= -35.5628 ZZ= -36.5632
XY= 0.0002 XZ= 0.0008 YZ= 0.0009
</pre>

*use the polar keyword on the excited state molecule
: use just polar to get the default analytic derivatives, in this case you might want to specify freq keyword as well
: you want to use the excited state density
: in this case polar=numerical computes the numerical derivative of the dipole moment
: you can change the step size of the electric field 0.0001N au with step=N
: you can restart a numerical calculation used polar=(numerical,restart)
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar=numerical td(read,root=1)
density=current guess=read int=ultrafine
</pre>
:you are looking for something like this
<pre>
Excited State 1: Singlet-A 1.6129 eV 768.69 nm f=0.0007 <S**2>=0.000
25 -> 26 -0.70872
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -391.407446678
</pre>
:and
<pre>
Electronic spatial extent (au): <R**2>= 592.2331
Charge= 0.0000 electrons
Dipole moment (field-independent basis, Debye):
X= -0.0001 Y= -2.0915 Z= 0.0158 Tot= 2.0915
Quadrupole moment (field-independent basis, Debye-Ang):
XX= -49.3534 YY= -35.5628 ZZ= -36.5632
XY= 0.0002 XZ= 0.0002 YZ= 0.0008
</pre>
:and
<pre>
Isotropic polarizability= 49.78 Bohr**3.
1 2 3
1 0.796564D+02
2 -0.303870D-03 0.508599D+02
3 0.349968D-03 0.169881D-03 0.188218D+02
</pre>

*of course we can also do it analytically, the default
<pre>
#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar=numerical td(read,root=1)
density=current guess=read int=ultrafine
</pre>
:and you should see
<pre>
Isotropic polarizability= 49.78 Bohr**3.
1 2 3
1 0.796564D+02
2 -0.292909D-03 0.508599D+02
3 0.360142D-03 0.138905D-03 0.188218D+02
</pre>
<pre>
(Enter /Applications/g09/l716.exe)
Dipole =-4.88131552D-05-8.22840998D-01-4.04220282D-05
Polarizability= 7.96563849D+01-2.92908823D-04 5.08598561D+01
3.60141967D-04 1.38905076D-04 1.88217717D+01
</pre>

*so what are the results?
:for our system the
::ground state dipole 2.74 D and isotropic polarizability 41 Bohr**3
::excited state dipole 2.09 D and isotropic polarizability 50 Bohr**3.
:according to the paper (M06-2X with Sadlej's basis set
::ground state dipole 1.02 D and polarizability 50 Bohr**3
::excited state dipole 0.66 D and polarizability 57 Bohr**3.
:so these differences are quite large!
:they could be due to
::basis set
::functional
::larger finite field
::convergence on energy
::convergence on the geometry

==Reproduction of the exact results==
Let us know try to reproduce the exact same results as described by Jacquemin and coworkers by applying the same procedure. First, optimise on the M06-2X/6-31G(d) level of theory using an improved energy threshold, a tight geometry optimisation criterion and a high level DFT integration grid.
<pre>
#p m062x/6-31g(d) opt=tight int=ultrafine scf=(conver=10)
</pre>
The optimised geometry is
<pre>
N -0.142626 0.000000 0.028260
N -0.072667 0.000000 1.258841
C 1.363634 0.000000 1.690124
N 2.090317 0.000000 0.518596
C 1.235531 0.000000 -0.563026
O 1.468028 0.000000 -1.729522
H 3.098469 0.000000 0.461269
O 1.726826 0.000000 2.822758
</pre>
with an energy of -391.3255482 Ha. Check that the geometry is a true minimum on the potential energy surface by conducting a frequency calculation
<pre>
#p m062x/6-31g(d) geom=checkpoint int=ultrafine scf=(conver=10)
guess=read freq
</pre>
(check for NImag=0).
The dipole moment and polarizability was calculated using the M06-2X functional and Sadlej's pVTZ basis set obtained from the EMSL Bais Set Exchange on the optimized geometry via
<pre>
#p m062x/gen int=ultrafine scf=(conver=10) guess=read polar
geom=checkpoint

H 0
S 4 1.00
33.8650140 0.0060680
5.0947880 0.0453160
1.1587860 0.2028460
0.3258400 0.5037090
S 1 1.00
0.1027410 1.0000000
S 1 1.00
0.0324000 1.0000000
P 2 1.00
1.1588000 0.1884400
0.3258000 0.8824200
P 2 1.00
0.1027000 0.1178000
0.0324000 0.0042000
****
C 0
S 5 1.00
5240.6353000 0.0009370
782.2048000 0.0072280
178.3508300 0.0363440
50.8159420 0.1306000
16.8235620 0.3189310
S 2 1.00
6.1757760 0.4387420
2.4180490 0.2149740
S 1 1.00
0.5119000 1.0000000
S 1 1.00
0.1565900 1.0000000
S 1 1.00
0.0479000 1.0000000
P 4 1.00
18.8418000 0.0138870
4.1592400 0.0862790
1.2067100 0.2887440
0.3855400 0.4994110
P 1 1.00
0.1219400 1.0000000
P 1 1.00
0.0385680 1.0000000
D 2 1.00
1.2067000 0.2628500
0.3855000 0.8043000
D 2 1.00
0.1219000 0.6535000
0.0386000 0.8636000
****
N 0
S 5 1.00
8104.0716000 0.0008020
1216.0215000 0.0061740
277.2342800 0.0312330
76.9040230 0.1151980
25.8744190 0.2969510
S 2 1.00
9.3467670 0.4473490
3.5797940 0.2450030
S 1 1.00
0.7396100 1.0000000
S 1 1.00
0.2226170 1.0000000
S 1 1.00
0.0670060 1.0000000
P 4 1.00
26.8689870 0.0144780
5.9912270 0.0911560
1.7508420 0.2974200
0.5605110 0.4937960
P 1 1.00
0.1759480 1.0000000
P 1 1.00
0.0552310 1.0000000
D 2 1.00
1.7508000 0.2247700
0.5605000 0.6595600
D 2 1.00
0.1795900 0.8713600
0.0552000 0.7042200
****
O 0
S 5 1.00
10662.2850000 0.0007990
1599.7097000 0.0061530
364.7252600 0.0311570
103.6517900 0.1155960
33.9058050 0.3015520
S 2 1.00
12.2874690 0.4448700
4.7568050 0.2431720
S 1 1.00
1.0042710 1.0000000
S 1 1.00
0.3006860 1.0000000
S 1 1.00
0.0900300 1.0000000
P 4 1.00
34.8564630 0.0156480
7.8431310 0.0981970
2.3062490 0.3077680
0.7231640 0.4924700
P 1 1.00
0.2148820 1.0000000
P 1 1.00
0.0638500 1.0000000
D 2 1.00
2.3062000 0.2027000
0.7232000 0.5791000
D 2 1.00
0.2149000 0.7854500
0.0639000 0.5338700
****
</pre>
yielding an energy of -391.4434344 Ha, a dipole moment of 1.02 au. and a polarizability of 50 au. (trace of the polarizability tensor)
<pre>
Dipole=1.0190643,0.,-0.0579513
Polar=45.7964203,0.,30.0056759,1.6195023,0.,74.1929892
</pre>
which corresponds exactly to the ground state gas phase results listed in the reference. Now optimise the excited state and check for imaginary frequencies via
<pre>
#p m062x/6-31g(d) opt=tight td(root=1) int=ultrafine scf=(conver=10)
guess=read geom=connectivity freq
</pre>
where the optimised geometry is
<pre>
N -0.099617 0.000000 0.023379
N -0.029402 0.000000 1.258821
C 1.333978 0.000000 1.722300
N 2.044096 0.000000 0.521225
C 1.202401 0.000000 -0.591640
O 1.500232 0.000000 -1.747732
H 3.054940 0.000000 0.463750
O 1.760884 0.000000 2.837197
</pre>
Dipole moment and polarizabilities can then again be calculated at the M06-2X/Sadlej level of theory
<pre>
#p m062x/gen td(root=1) int=ultrafine scf=(conver=10) guess=read
polar=numerical density=current geom=checkpoint
</pre>
yielding
<pre>
Dipole=0.6557277,0.,-0.03725
Polar=54.0435904,0.,29.7590519,1.8414148,0.,86.3255147
</pre>
which is a total dipole moment of 0.66 au. and a polarizability of 57 au., which corresponds perfectly to the reference. An alternative approach to calculate the excited state polarizability is via the use of Stark's relation, which links the energy of a state to the strength of an applied electric field.
<math display="block"> E(F)=E^0-\sum_a \mu_a F_a - \frac{1}{2}\sum_{a,b} \alpha_{ab} F_a F_b </math>
The dipole moment is then the first derivative of the energy with respect to the field at zero field strength. Make sure to take into account all directions, as <math display="inline">\mu_x=\big( \frac{\mathrm{d}E}{\mathrm{d}F_x} \big) _{F_x=0} </math>,
<math display="inline">\mu_y=\big( \frac{\mathrm{d}E}{\mathrm{d}F_y} \big) _{F_y=0} </math> and <math display="inline">\mu_z=\big( \frac{\mathrm{d}E}{\mathrm{d}F_z} \big) _{F_z=0} </math> and the total dipole moment is then <math display="inline">\mu=\sqrt{\mu_x^2+\mu_y^2+\mu_z^2}</math>. Analogously, the polarizability is obtained via <math display="inline">\alpha_{xx}=\big( \frac{\mathrm{d}^2E}{\mathrm{d}F_x^2} \big) _{F_x=0} </math>, <math display="inline">\alpha_{yy}=\big( \frac{\mathrm{d}^2E}{\mathrm{d}F_y^2} \big) _{F_y=0} </math>, <math display="inline">\alpha_{zz}=\big( \frac{\mathrm{d}^2E}{\mathrm{d}F_z^2} \big) _{F_z=0} </math> and <math display="inline">\alpha=\frac{\alpha_{xx}+\alpha_{yy}+\alpha_{zz}}{3}</math>. The electric fields are chosen such that the Romberg differentiation procedure [http://dx.doi.org/10.1016/j.theochem.2007.08.028 link] can be applied.
<math display="block"> F_a=2^k \cdot 0.0004 \; \mathrm{au.} \;\;\;\;\;\; \mathrm{with} \;k=1,2,...5 \;\;\; \mathrm{for} \; a=x,y,z</math>
To apply for example a field of 0.0008 au. in the x-direction, use
<pre>
#p m062x/gen td(root=1) int=ultrafine scf=(conver=10)
guess=read density=current field=x+8 geometry=checkpoint
</pre>
Search the output for the energy of the first excited state
<pre>
Excited State 1: Singlet-A'' 1.6794 eV 738.27 nm f=0.0008
<S**2>=0.000
25 -> 26 0.70425
25 <- 26 -0.10327
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -391.373283614
</pre>
Here, we get <math display="inline">E(F_x=0.0008)=</math>-391.373283614 Ha. Repeat the calculation for all field strengths, which finally yields
<pre>
F [au.] E(Fx=F) [au.] E(Fy=F) [au.] E(Fz=F) [au.]
-0.0128 -391.380340688 -391.386106671 -391.375696126
-0.0064 -391.375024593 -391.378566574 -391.373865503
-0.0032 -391.373697934 -391.375634101 -391.373408310
-0.0016 -391.373366469 -391.374375865 -391.373294027
-0.0008 -391.373283571 -391.373798627 -391.373265458
0.0000 -391.373255935 -391.373255935 -391.373255935
0.0008 -391.373283614 -391.372747763 -391.373265458
0.0016 -391.373366555 -391.372274081 -391.373294027
0.0032 -391.373698106 -391.371430018 -391.373408310
0.0064 -391.375024938 -391.370154501 -391.373865503
0.0128 -391.380341380 -391.369251683 -391.375696126
</pre>
The Romberg differentiation then uses
<math display="block"> \big( \frac{\mathrm{d}^2f}{\mathrm{d}x^2} \big)_{x=0}= \frac{f(-2^kh)+f(2^kh)-2f(0)}{(2^kh)^2} </math>
where h is the minimum step size (here 0.0004 au.) to evaluate the second derivative and then uses the iterative formula
<math display="block"> P^{p,k} = \frac{4^pP^{p-1,k}-P^{p-1,k+1}}{4^p-1} </math>
where p is the number of Romberg iteration and P is the required differentiation, to improve the obtained value. The value at p=0 corresponds to the initially obtained uncorrected differentiation value. The iterative procedure
<pre>
k 2^k h d^2E/dFx^2 p=1 p=2 p=3 p=4
1 0.0008 -86.4 -86.4 -86.4 -86.4 -86.4
2 0.0016 -86.4 -86.4 -86.4 -86.4
3 0.0032 -86.3 -86.3 -86.3
4 0.0064 -86.4 -86.3
5 0.0128 -86.5
</pre>
then yields <math display="inline">\alpha_{xx}</math>= 86 au. Analogue evaluation in y and z direction yields <math display="inline">\alpha_{yy}</math>= 54 au and <math display="inline">\alpha_{zz}</math>= 30 au. The polarizability is then
57 au, which is the same as the value in Reference, as well as the value obtained using the default routine for the
calculation of polarizability in Gaussian09.

Finally, try to reproduce the correct dipole moment and polarisability in implicit solvent (dichloromethane). Optimise the ground state via
<pre>
#p m062x/6-31g(d) opt=tight int=ultrafine scf=(conver=10)
SCRF=(IEFPCM,SOLVENT=dichloromethane) guess=read
</pre>
yielding the geometry
<pre>
N -0.141332 0.000000 0.028556
N -0.071415 0.000000 1.258401
C 1.367139 0.000000 1.684979
N 2.090981 0.000000 0.518558
C 1.239597 0.000000 -0.558312
O 1.461088 0.000000 -1.729021
H 3.101578 0.000000 0.461093
O 1.719875 0.000000 2.823046
</pre>
at an energy of -391.335724 Ha. Dipole moments and polarizabilities at the optimized geometry at the M06-2X/Sadlej level of theory were then obtained via
<pre>
#p m062x/gen int=ultrafine scf=(conver=10) guess=read polar
SCRF=(IEFPCM,SOLVENT=dichloromethane) geom=checkpoint
</pre>
to yield the following output:
<pre>
Dipole=1.2615935,0.,-0.0717451
Polar=55.4963369,0.,35.741816,2.3407747,0.,96.5379704
</pre>
The dipole moment in groundstate is therefore 1.26 au. and the
polarizability 63 au.

The excited state in implicit solvent was also optimized:
<pre>
#p m062x/6-31g(d) opt=tight td(root=1) int=ultrafine
scf=(conver=10) guess=read geom=(connectivity,checkpoint)
SCRF=(IEFPCM,SOLVENT=dichloromethane)
</pre>
using the optimized excited state geometry in gas phase as input. The optimized geometry is then
<pre>
N -0.101348 0.000000 0.023438
N -0.031138 0.000000 1.258959
C 1.337822 0.000000 1.716250
N 2.046763 0.000000 0.521074
C 1.206893 0.000000 -0.586065
O 1.494126 0.000000 -1.746788
H 3.059686 0.000000 0.463479
O 1.754707 0.000000 2.836953
</pre>
Calculation of the dipole moment and the polarizability is then done on the M06-2X/Sadlej level of theory using the default calculation routine.
<pre>
#p m062x/gen td(root=1) int=ultrafine scf=(conver=10)
guess=read polar=numerical density=current geom=checkpoint
SCRF=(IEFPCM,SOLVENT=dichloromethane)
</pre>
The dipole moment is then 0.84 au. and the polarizability 73 au. which is for both ground and excited state in perfect agreement with the published results in the paper by Jacquemin and coworkers.
<pre>
Dipole=0.8376841,0.,-0.047563
Polar=67.8846862,0.,35.5510398,2.7715247,0.,116.4748693
</pre>

==Atomic polarizabilities==
We now want to calculate the atomic contributions to the molecular polarizability. Again, we use Stark's relation, namely that the change in dipole moment with respect to an external field is the polarizability. Changes in dipoles arise from atomic dipoles (asymmetric distribution of charge around the nucleus) that change upon applying an external field and from charge transfer between the atoms. The GDMA program written by Anthony Stone can calculate atomic dipoles
and partial charges at specific sites in the molecule, here the nuclei. We reuse the calculations at external fields of 0.0008au in negative and positive x, y and z direction from above. The respective checkpoint files need to be converted to formatted checkpoint files. From these the atomic dipoles can be calculated using the following input script for the ground state:
<pre>
Title "title"
File file.fchk

Angstrom
Multipoles
Limit 1
Start

Finish
</pre>
For the excited state, the electron density needs to be specified by putting the keyword "DENSITY "CI" at the end of the "File" line. This input script tells GDMA to calculate multipoles up to the rank of dipoles at all nuclear sites. The output (here for Fx=0.0008au) looks like:
<pre>
[...]
Positions and radii in angstrom
Multipole moments in atomic units, ea_0^k for rank k

N x = -0.616284 y = -1.257223 z = 0.000000 angstrom
Maximum rank = 1 Radius = 0.650 angstrom
Q00 = -0.056538
|Q1| = 0.511629 Q11c = 0.331633 Q11s = 0.389594

N x = 0.616284 y = -1.257222 z = 0.000000 angstrom
Maximum rank = 1 Radius = 0.650 angstrom
Q00 = -0.060705
|Q1| = 0.516025 Q11c = -0.335883 Q11s = 0.391746
[...]
x = 2.279816 y = 0.450604 z = 0.000000 angstrom
Maximum rank = 1 Radius = 0.650 angstrom
Q00 = -0.280276
|Q1| = 0.425044 Q11c = -0.405995 Q11s = -0.125819

Total multipoles referred to origin at
x = 0.000000, y = 0.000000, z = 0.000000 angstrom
Q00 = -0.000002
|Q1| = 1.022461 Q11c = -0.059386 Q11s = 1.020735
</pre>
where [...] indicates that not all data is shown. For the first nitrogen atom, a atomic dipole (0.331633,0.389594,0) is observed (Q11c is the x direction, Q11s the y direction and Q10 the z direction). For the calculation of atomic polarizability from atomic dipoles we need to monitor the x-component of the dipole when applying a field in x-direction, the y component for fields in y-direction and so on. For the nitrogen atom we get:
<pre>
Fx= 0.0008au , mu_x= 0.331633
Fx=-0.0008au , mu_x= 0.335880
Fy= 0.0008au , mu_y= 0.387470
Fy=-0.0008au , mu_y= 0.393892
Fz= 0.0008au , mu_z=-0.003752
Fz=-0.0008au , mu_z= 0.003752
</pre>
Numerical differentiation yields <math display="inline">\alpha_{xx}=-\big( \frac{\mathrm{d}\mu_x}{\mathrm{d}F_x} \big) _{F_x=0} = -\frac{\mu_x(F_x)-\mu_x(-F_x)}{2F_x} = 2.65 au^3</math>, and analogously
<math display="inline">\alpha_{yy}=4.01 au^3</math> and <math display="inline">\alpha_{zz}=4.69 au^3</math> yielding an average polarizability of 3.8au^3. The same is done for all other atoms and then repeated for the excited state. One then ends up with the atomic polarizabilities
<pre>
Atomic alpha from polarization:
GS ES
N 3.8 4.1
N 3.8 4.1
C 1.2 1.4
N 2.5 2.3
N 1.2 1.4
O 4.6 4.7
H 0.8 0.8
O 4.6 4.7
</pre>
Summed up, we end up with a polarizability from atomic dipoles of 22.4 au^3 in the ground state and 23.4 au^3 in the excited state. The missing 27.5 au^3 in the ground and 33.3 au^3 in the excited state stem from charge transfer. The overall contribution of charge transfer over all atoms is easy to calculate via <math display="inline">\mu_{charge transfer}=\sum_ir_iq_i</math> where <math display="inline">r_i</math> are the coordinates of each atom i and <math display="inline">q_i</math> the charge at that site. The charge at a specific site can be found in the GDMA output script and is labelled "Q00". The coordinates are also printed in the output script. Summing up these values gives the desired polarizability of 27.5 au^3 in the ground and 33.3 au^3 in the excited state. However, the separation of this contribution to different atomic sites is difficult, as the sum <math display="inline">\mu_ct=\sum_ir_iq_i</math> is origin independent (for a neutral molecule) but the respective contributions <math display="inline">r_iq_i</math> aren't. A possible solution to this problem is the definition of bond charges <math display="inline">b_{i,j}</math> between the bonded atoms i and j. Each partial charge is assumed to be the sum of all bond charges, the sum of all bond charges in a ring is zero and same bond charges in different directions add up to 0 (<math display="inline">b_{i,j}+b_{j,i}=0</math>). These relations can be used to set up linear independent equations, which are solved for <math display="inline">b_{i,j}</math>. The atomic contribution to the dipole moment arising from charge transfer is then <math display="inline">\mu_{charge transfer, i}=\sum_j(r_i-\frac{r_i-r_{i,j}}{2})*b_{i,j}</math> where <math display="inline">r_{i,j}</math> was simply assumed to be <math display="inline">\frac{r_i+r_j}{2}</math>. The same procedure as explained above for the atomic polarization contribution is applied here: the atomic charge transfer dipole moments are calculated at different fields and differentiated numerically to yield the average polarizabilities from charge transfer
<pre>
Atomic alpha from charge transfer
GS ES
N 3.3 4.1
N 3.3 4.1
C 5.7 6.9
N 4.5 5.1
N 5.7 6.9
O 2.2 2.7
H 0.7 0.6
O 2.2 2.7
</pre>
To obtain atomic polarizabilities we now sum up the contributions from polarization and charge transfer
<pre>
Overall atomic alpha
GS ES
N 7.1 8.2
N 7.1 8.2
C 6.8 8.4
N 7.0 7.5
N 6.8 8.4
O 6.9 7.4
H 1.5 1.4
O 6.9 7.4
</pre>
In this molecule, the atomic polarization increases slightly upon excitation (in average about 6%), and charge transfer polarizability to a larger extent (about 19%).

Talk:Mod:Hunt Research Group/VMDTips

2017-02-08T17:00:24Z

Ech17: /* Color by random property */ new section

=Tricks and Tips=
==quick way to identify commands==
*Turn on the TkConsole and log commands
:go to Extensions and turn on the TkConsole, then goto File select log Tcl commands to Console, this logs the command of anything you do on the GUI and thus allows you to quickly determine a number of commands for scripts

==hide a molecule==
*double click on the "D" under molecule in the main VMD window

==change one of the predefined colours==
*this is the color command [http://www.ks.uiuc.edu/Research/vmd/current/ug/node123.html link]
*must be applied to an existing colour name
*format is where color=color_name and r g b are numerical values (can be read off the slider)
:change rgb color r g b
*example which changes the colour blue2 to a light blue is
:color change rgb blue2 0.500000 0.700000 0.970000

== Color by random property ==

Sometimes it is convenient to color the atoms of a molecule according to a specific property. VMD provides the possibility to do so, namely via the implemented color option according to the B-factor of an atom. To do so, make a PDB file of the molecule of interest, e.g. water
<pre>
ATOM 1 OH OSP3 1 4.013 0.831 -9.083 1.00 0.00
ATOM 2 1HH OSP3 1 4.941 0.844 -8.837 1.00 0.00
ATOM 3 2HH OSP3 1 3.750 -0.068 -9.293 1.00 0.00
TER
</pre>
and change the B-factor array (which is the last column, here all 0.00) to the desired property, for example
<pre>
ATOM 1 OH OSP3 1 4.013 0.831 -9.083 1.00 2.00
ATOM 2 1HH OSP3 1 4.941 0.844 -8.837 1.00 1.00
ATOM 3 2HH OSP3 1 3.750 -0.068 -9.293 1.00 3.00
TER
</pre>
Save the changes, then open the PDB file in VMD.
Click on <tt>Graphics</tt> => <tt>Representations</tt>, then select "Beta" in the <tt>Coloring Method</tt> drop-down menu. The color scale data range can be changed when clicking on <tt>Trajectory</tt> in the same window and changing the values in <tt>Color Scale Data Range</tt>.
[[File:WIKI_1.png|center|Graphical Representations Window]]
Thus, every property, such as partial charge change, rmsd, how much you like an atom :) ... can be visualised easily.
[[File:WIKI_2.png|center|Colored water molecule according to the property saved to the B-factor array.]]

File:WIKI 2.png

2017-02-08T16:57:02Z

Ech17:

File:WIKI 1.png

2017-02-08T16:53:25Z

Ech17:

Mod:Hunt Research Group: ES alpha

2017-02-03T15:45:13Z

Ech17: /* Reproduction of the exact results */

Mod:Hunt Research Group: ES alpha

2017-02-02T16:34:26Z

Ech17: