View source for Mod:Hunt Research Group: ES alpha

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

Thus, the following linear equations evolve:

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

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

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 <br />
<math display="inline">\mu_{charge transfer, 2}</math> = [  0.023, 0.076, 0.0000] in au <br />
...<br />

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