Mod:Hunt Research Group: ES alpha

introduction

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

useful blog
data below relates to information in this paper 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)

#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) guess=read int=ultrafine freq

we turn on extra printing with #p

SCF Done: E(RB3LYP) = -391.474123816

you are looking for something like this

 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

the dipole moment is in the standard orientation, and the zero is therefore the CoM

then

 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

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

#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar guess=read int=ultrafine

you will see something like this

were clearly the tensor is symmetric and 1=y, 2=z and 3=x

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

and a bit later this

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

and a bit later this

 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

or you can do both!

#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) guess=read int=ultrafine polar freq

you will see this

 Isotropic polarizability for W=    0.000000       40.67 Bohr**3.

and later this

  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

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

 (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

now optimise on the first excited state, use the time-dependent options, the root=1 is the first excited state

# opt td(root=1) b3lyp/6-31g(d,p) geom=connectivity scf=(conver=9) guess=read int=ultrafine

you should see something like this in the output

 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

carry out a frequency analysis on the optimised excited state geometry

# freq b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) td(read,root=1) guess=read int=ultrafine

this will print out excited state dipole moments and so on ...

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

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)

#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar=numerical td(read,root=1) 
  density=current guess=read int=ultrafine

you are looking for something like this

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

and

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

and

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

of course we can also do it analytically, the default

#p b3lyp/6-31g(d,p) geom=checkpoint scf=(conver=9) polar=numerical td(read,root=1) 
  density=current guess=read int=ultrafine

and you should see

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

 (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

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.

#p m062x/6-31g(d) opt=tight int=ultrafine scf=(conver=10)

The optimised geometry is

  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

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

#p m062x/6-31g(d) geom=checkpoint  int=ultrafine scf=(conver=10) 
guess=read freq

(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

#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
****

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)

  Dipole=1.0190643,0.,-0.0579513
  Polar=45.7964203,0.,30.0056759,1.6195023,0.,74.1929892

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

#p m062x/6-31g(d) opt=tight td(root=1)  int=ultrafine scf=(conver=10) 
guess=read geom=connectivity freq

where the optimised geometry is

  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

Dipole moment and polarizabilities can then again be calculated at the M06-2X/Sadlej level of theory

#p m062x/gen td(root=1)  int=ultrafine scf=(conver=10) guess=read 
 polar=numerical density=current geom=checkpoint

yielding

  Dipole=0.6557277,0.,-0.03725
  Polar=54.0435904,0.,29.7590519,1.8414148,0.,86.3255147

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. $E (F) = E^{0} - \sum_{a} μ_{a} F_{a} - \frac{1}{2} \sum_{a, b} α_{a b} F_{a} F_{b}$ 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 $μ_{x} = (\frac{d E}{d F_{x}})_{F_{x} = 0}$ , $μ_{y} = (\frac{d E}{d F_{y}})_{F_{y} = 0}$ and $μ_{z} = (\frac{d E}{d F_{z}})_{F_{z} = 0}$ and the total dipole moment is then $μ = \sqrt{μ_{x}^{2} + μ_{y}^{2} + μ_{z}^{2}}$ . Analogously, the polarizability is obtained via $α_{x x} = (\frac{d^{2} E}{d F_{x}^{2}})_{F_{x} = 0}$ , $α_{y y} = (\frac{d^{2} E}{d F_{y}^{2}})_{F_{y} = 0}$ , $α_{z z} = (\frac{d^{2} E}{d F_{z}^{2}})_{F_{z} = 0}$ and $α = \frac{α_{x x} + α_{y y} + α_{z z}}{3}$ . The electric fields are chosen such that the Romberg differentiation procedure link can be applied. $F_{a} = 2^{k} \cdot 0.0004 a u . w i t h k = 1, 2, . . . 5 f o r a = x, y, z$ To apply for example a field of 0.0008 au. in the x-direction, use

#p m062x/gen td(root=1)  int=ultrafine scf=(conver=10) 
guess=read  density=current field=x+8 geometry=checkpoint

Search the output for the energy of the first excited state

 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

Here, we get $E (F_{x} = 0.0008) =$ -391.373283614 Ha. Repeat the calculation for all field strengths, which finally yields

 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

The Romberg differentiation then uses $(\frac{d^{2} f}{d x^{2}})_{x = 0} = \frac{f (- 2^{k} h) + f (2^{k} h) - 2 f (0)}{(2^{k} h)^{2}}$ where h is the minimum step size (here 0.0004 au.) to evaluate the second derivative and then uses the iterative formula $P^{p, k} = \frac{4^{p} P^{p - 1, k} - P^{p - 1, k + 1}}{4^{p} - 1}$ 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

  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

then yields $α_{x x}$ = 86 au. Analogue evaluation in y and z direction yields $α_{y y}$ = 54 au and $α_{z z}$ = 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

#p m062x/6-31g(d) opt=tight  int=ultrafine scf=(conver=10)
  SCRF=(IEFPCM,SOLVENT=dichloromethane) guess=read

yielding the geometry

  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

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

  #p m062x/gen int=ultrafine scf=(conver=10) guess=read polar
  SCRF=(IEFPCM,SOLVENT=dichloromethane) geom=checkpoint

to yield the following output:

  Dipole=1.2615935,0.,-0.0717451
  Polar=55.4963369,0.,35.741816,2.3407747,0.,96.5379704

The dipole moment in groundstate is therefore 1.26 au. and the polarizability 63 au.

The excited state in implicit solvent was also optimized:

  #p m062x/6-31g(d) opt=tight td(root=1)  int=ultrafine
  scf=(conver=10) guess=read geom=(connectivity,checkpoint)
  SCRF=(IEFPCM,SOLVENT=dichloromethane)

using the optimized excited state geometry in gas phase as input. The optimized geometry is then

  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

Calculation of the dipole moment and the polarizability is then done on the M06-2X/Sadlej level of theory using the default calculation routine.

  #p m062x/gen td(root=1)  int=ultrafine scf=(conver=10)
  guess=read polar=numerical density=current geom=checkpoint
  SCRF=(IEFPCM,SOLVENT=dichloromethane)

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.

  Dipole=0.8376841,0.,-0.047563
  Polar=67.8846862,0.,35.5510398,2.7715247,0.,116.4748693

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:

Title "title"
File file.fchk 

Angstrom
Multipoles
  Limit 1
Start

Finish

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:

[...]
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

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:

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

Numerical differentiation yields $α_{x x} = - (\frac{d μ_{x}}{d F_{x}})_{F_{x} = 0} = - \frac{μ_{x} (F_{x}) - μ_{x} (- F_{x})}{2 F_{x}} = 2.65 a u^{3}$ , and analogously $α_{y y} = 4.01 a u^{3}$ and $α_{z z} = 4.69 a u^{3}$ 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

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

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 $μ_{c h a r g e t r a n s f e r} = \sum_{i} r_{i} q_{i}$ where $r_{i}$ are the coordinates of each atom i and $q_{i}$ 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 $μ_{c} t = \sum_{i} r_{i} q_{i}$ is origin independent (for a neutral molecule) but the respective contributions $r_{i} q_{i}$ aren't. A possible solution to this problem is the definition of bond charges $b_{i, j}$ 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 ( $b_{i, j} + b_{j, i} = 0$ ). These relations can be used to set up linear independent equations, which are solved for $b_{i, j}$ .

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

N(1)-N(2) => $b_{1, 2}$
N(1)-C(5) => $b_{1, 5}$
N(2)-C(3) => $b_{2, 3}$
C(3)-N(4) => $b_{3, 4}$
C(3)-O(8) => $b_{3, 8}$
N(4)-C(5) => $b_{4, 5}$
N(4)-H(7) => $b_{4, 7}$
C(5)-O(6) => $b_{5, 6}$

Thus, the following linear equations evolve:

q(N1)= $b_{1, 2} - b_{1, 5}$
q(N2)= - $b_{1, 2} + b_{2, 3}$
q(C3)= - $b_{2, 3} + b_{3, 4} + b_{3, 8}$
q(N4)= - $b_{3, 4} + b_{4, 5} + b_{4, 7}$
q(C5)= - $b_{4, 5} - b_{1, 5} + b_{5, 6}$
q(O6)= - $b_{5, 6}$
q(H7)= - $b_{4, 7}$
q(O8)= - $b_{3, 8}$
$0 = b_{1, 2} + b_{2, 3} + b_{3, 4} + b_{4, 5} - b_{1, 5}$ (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 $b_{i, j}$ . We get for the field in positive x-direction

$b_{1, 2}$ = 0.0036
$b_{1, 5}$ = -0.0602
$b_{2, 3}$ = -0.0571
$b_{3, 4}$ = -0.0538
$b_{3, 8}$ = 0.2803
$b_{4, 5}$ = 0.0471
$b_{4, 7}$ = -0.2531
$b_{5, 6}$ = 0.2709

The atomic contribution to the dipole moment arising from charge transfer is then $μ_{c h a r g e t r a n s f e r, i} = \sum_{j} (r_{i} - \frac{r_{i} - r_{i, j}}{2}) * b_{i, j}$ where $r_{i, j}$ was simply assumed to be $\frac{r_{i} + r_{j}}{2}$ , thus we get

$μ_{c h a r g e t r a n s f e r, 1}$ = [ -0.033, 0.080, 0.0000] in au
$μ_{c h a r g e t r a n s f e r, 2}$ = [ 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

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

To obtain atomic polarizabilities we now sum up the contributions from polarization and charge transfer

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

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