# QTPIE and water (Part 2)

33 %
67 %
Information about QTPIE and water (Part 2)

Published on June 16, 2008

Source: slideshare.net

## Description

Slides for group meeting presentation, Spring 2008.

QTPIE and water Jiahao Chen 2008-05-20 quot;We can't solve problems by using the same kind of thinking we used when we created them.quot; - attributed to Albert Einstein

Since last time... • QTPIE has become much faster • We now know why dipole moments and polarizabilities previously weren’t translationally invariant, and why they aren’t size extensive. • We have (some) parameters for a new water model • We’ve shown that QTPIE gets the correct direction of intermolecular charge transfer

What is QTPIE: a scientiﬁc POV • A way to model polarization and intermolecular charge transfer in molecular mechanics • One of the simplest electronic structure methods, except without electrons • Give me a geometry, and I will give you a charge distribution

What is QTPIE? A numerical POV • Give me a geometry, and I will give you a charge distribution f : R = R1 , . . . , RN → (q1 , . . . , qN ) • Minimize the quadratic form1 E(q1 , . . . , qN ; R) = qi vi (R) + qi qj Jij (R) i 2 ij • subject to the constraint qi = 0 i

How to make QTPIE faster • Using GTOs in place of STOs • Integral prescreening • Sparse matrix data structure for overlap matrix • Conjugate gradients to solve linear problem (vs GMRES) • Initial guess from previous solution • Fast multipole methods (future work)

Optimizing GTOs for STOs

n n!e−α αν An (α) = αn+1 ν=0 ν! n ν n!e−α (−α) − αν Bn (α) = αn+1 ν=0 ν! k m n mn Dp = (−1) p−k k k 2m 1 (2α) R2m−1 K2 (α, β, m, n, R) = + (2αR − 2m) A2m−1 (2αR) − e−2αR R (2m)! 2n 2m−1+ν α2m+1 R2m 2n − ν ν − (βR) 2m−1,ν Dp Bp (R (α − β)) A2m+ν−1−p (R (α + β)) 2n (2m)! ν=0 ν! p=0 2m dK2 2m + 1 2m (2αR) = − + K2 + 2m (1 + 2m − 2αR) A2m−1 (2αR) + (1 + 2m) e−2αR dR R 2 R (2m)! 2n−1 2m+ν α2m+1 R2m 2n − ν − 1 ν −β (βR) 2m−1,ν+1 Dp Bp (R (α − β)) A2m+ν−p (R (α + β)) 2n (2m)! ν=0 ν! p=0 2m 2n 2m−1+ν α2m+1 R 2n − ν ν + (βR) 2m−1,ν Dp × 2n (2m)! ν=0 ν! p=0 [(α − β) Bp+1 (R (α − β)) A2m+ν−1−p (R (α + β)) + (α + β) Bp (R (α − β)) A2m+ν−p (R (α + β))] Unlike GTOs, STOs are really, really nasty to work with. ab p = a+b dK2 2p −p2 R2 erf (pR) = √ e − dR R π R2

STO-nG orbitals were deﬁned by ﬁtting an orbital of n contracted Gaussian primitives to an STO to reproduce the orbital in a least-squares sense. The conventional wisdom is n > 2 for some kind of useful, reasonable ﬁt. STO-1Gs suck when used in QEq/QTPIE. However, there is a way to get accurate results with just primitive Gaussians...

Instead of maximizing overlap in a least-squares sense, minimize the deviation in the Coulomb integrals. min J ST O − J GT O (α) α This is equivalent to minimizing J GT O − J ST O 2 2 = J GT O , J GT O − 2 J GT O , J ST O + J ST O , J ST O = J GT O , J GT O − 2J ST O + J ST O , J ST O which is minimized by the exponent such that ∂ 0 = J GT O (α) , J GT O (α) − 2J ST O ∂α ∂ GT O ∂ GT O = J (α) , J GT O (α) − 2J ST O + J GT O (α) , J (α) ∂α ∂α ∂ GT O = 2J GT O (α) − 2J ST O , J (α) ∂α i.e. ∞ 1 −αR2 J GT O (R; α) − J ST O (R) √ e =0 0 απ

Element Best Coulomb Least-squares H 0.5343 0.3101 Li 0.1668 0.0440 C 0.2069 0.1853 N 0.2214 0.2088 O 0.2240 0.2400 F 0.2312 0.2142 Na 0.0959 0.0399 Si 0.1032 0.1256 P 0.1085 0.1430 S 0.1156 0.1584 Cl 0.1137 0.1758 K 0.0602 0.0361 Br 0.0701 0.1850 Rb 0.0420 0.0402 I 0.0469 0.1735 Cs 0.0307 0.0374

Summary • STOs are really, really nasty to work with • GTOs ﬁt to STOs by reproducing the Coulomb self-repulsion integrals yield excellent approximations • Use in QEq and QTPIE result in very little error (< 0.00001e) • There is very little basis to the claim that STOs can be used “for extra accuracy”

Sparse matrices in QTPIE

Matrix-vector multiplication in computer memory   1.0 0.4 0.0 do i = 1, N do j = 1, N  0.4 1.0 0.0  u(i) = M(i, j) * v(j) end do 0.0 0.0 1.0 end do In (linear) memory, matrix data structure looks like this: M(:) 1.0 0.4 0.0 0.4 1.0 0.0 0.0 0.0 1.0 k = 0 do i = 1, N do j = 1, N k = k + 1 u(i) = M(k) * v(j) end do end do

Conventional matrix data structure 1.0 0.4 0.0 0.4 1.0 0.0 0.0 0.0 1.0   k = 0 1.0 0.4 0.0 do i = 1, N do j = 1, N  0.4 1.0 0.0  k = k + 1 u(i) = M(k) * v(j) 0.0 0.0 1.0 end do end do Compressed sparse row (CSR) data structure row start 1 3 5 6 do i = 1, N do k=M%r(i),M%r(i+1)-1 j = M%c(k) column 1 2 1 2 3 u(i) = M%d(k) * v(j) end do data 1.0 0.4 0.4 1.0 1.0 end do Fewer operations and lower memory latency, so faster!

Calculating the linear coefﬁcients in QTPIE • In QTPIE, minimize 1 E(q1 , . . . , qN ; R) = qi vi (R) + qi qj Jij (R) i 2 ij • The linear coefﬁcients are given by j (χi − χj )Sij vi = Sij • The main costs are matrix-vector j multiplication and memory latency

Using conjugate gradients for constrained problems

QTPIE is a saddle-point problem • In matrix notation, 1 T min q v + q Jq T q: q·1=0 2 • Solving QTPIE with a Lagrange multiplier, J 1 q −v = 1T 0 µ 0 • J is positive deﬁnite but this constrained problem is not (the constraint introduces a negative eigenvalue). Conjugate gradients can (and does) fail by going uphill, thinking it’s going downhill.

However...

Solving the saddle-point problem • There exists a block inversion formula1 for 2 x 2 structured matrices −1 J 1 J −1 + J −1 1S1T J −1 −J −1 1S = 1T 0 −S1T J −1 S 1 S = − T −1 1 J 1 • The analytic solution is given by −1 q J 1 −v −J −1 (v − µ1) = = µ 1T 0 0 1T J −1 v S 1. e.g. P.-O. Löwdin, Linear algebra for quantum mechanics

Solving the saddle-point problem with CG • Analytic solution can be solved numerically with two symmetric, positive deﬁnite problems (CG is guaranteed to work) q −J −1 (v − µ1) Jw = 1 = w·v µ 1T J −1 v S µ=− w·1 1 Jy = −v S = − T −1 1 J 1 q = y − µw

II. Fixing the translational invariance of dipole moments and polarizabilities

Dipoles and polarizabilities in QEq • The regular story: the charge model is solved by q = −J −1 v • The energy is therefore minimized by 1 E = − v T J −1 v 2 • Want to calculate dipole moments and polarizabilities ∂E ∂dν dν = ανλ = ∂Riν j ∂Rjλ i

How to calculate dipoles and polarizabilities • Dipole coupling prescription 1 E(q1 , . . . , qN ; R, ) = qi vi (R) + qi qj Jij (R) − q i Ri · i 2 ij i • The external ﬁeld shifts the voltages on each atom by an external potential 1 E(q1 , . . . , qN ; R, ) = qi vi (R) − Ri · + qi qj Jij (R) i 2 ij • Now calculate dipole moments and polarizabilities ∂E ∂2E dν = ανλ = ∂ ν ∂ ν∂ λ

Physically, the universe is translationally invariant. Therefore, electrostatic properties (mostly) don’t depend on the choice of origin. x→x+δ ⇒ d → d + δQ ανλ → ανλ

This is not the case in QEq! d→d+δ 1 J T −1 v ανλ → ανλ − 1 JT −1 · (δν Rλ + Rν δλ ) − δν δλ 1 J T −1 1 The solution in the literature is to ﬁx an origin arbitrarily, even though nobody has a good physical reason why that should be the “correct” origin

As it turns out, QEq and other ﬂuctuating charge models do obey the correct translational properties, as long as one works with the correct solution of the constrained minimization q = −J −1 v J −1 1 q = −J −1 (v − µ1) = −J −1 v − T −1 1 J 1 The second term actually generates counterterms under translation that kill all the pathological terms.

The analytic solutions for the dipole moment and polarizability in QTPIE are given by µ T −1 1T J−1 Rµ dµ = − (R ) J v − 1T J−1 v + Q 1T J−1 1 µ T −1 ν 1T J−1 Rµ 1T J−1 Rν αµν = − (R ) J R − 1T J−1 1

T dµ → − (R + δ 1) J−1 v µ µ 1T J−1 − 1T J−1 v + Q T −1 (Rµ + δ µ 1) 1 J 1 T −1 1T J−1 Rµ = − (Rµ ) J v + 1T J−1 v + Q 1T J−1 1 1T J−1 1 −δ µ 1T J−1 v + δ µ 1T J−1 v + Q T −1 1 J 1 = dµ − δ µ 1T J−1 v + δ µ 1T J−1 v + Q = dµ + δ µ Q

T αµν → − (Rµ + δ µ 1) J−1 (Rν + δ ν 1) 1T J−1 (Rµ + δ µ 1) 1T J−1 (Rν + δ ν 1) − 1T J−1 1 T = − (Rµ ) J−1 Rν − δ µ δ ν 1T J−1 1 T −δ µ 1T J−1 Rν − δ ν (Rµ ) J−1 1 1T J−1 Rµ 1T J−1 Rν − 1T J−1 1 1T J−1 1 1T J−1 1 −δ µ δ ν 1T J−1 1 1T J−1 1 1T J−1 Rν −δ µ 1T J−1 1 1T J−1 Rµ 1T J−1 1 −δ ν 1T J−1 1 = αµν

III. Fixing the size extensivity of dipole moments and polarizabilities

Not size extensive! 2n lim = 2 n→∞ 4 n Why?

Work it out analytically for N identical, noninteracting subsystems...       ¯ Rµ ¯ Rµ 0  ¯ Rµ + ∆µ 1   Rµ ¯   1      µ  R = µ . = . +∆  .   . .   . .   . .  ¯ Rµ + (N − 1) ∆µ 1 ¯ Rµ (N − 1) 1  ¯ 0 ··· 0    J ¯ v  .. .   0 J ¯ . .  .    ¯ v   J=  . .. .  v= . .  . ..    . . . .  . . 0 ··· ··· J ¯ ¯ v

T 1 dµ = − (Rµ ) J−1 v − T −1 1T J−1 v + Q 1 J 1    −1    ¯µ R ¯ v J ¯ ¯ J−1 1  ¯ Rµ + ∆µ 1 ¯ ¯   J−1 v  N 1T J−1 v + Q  ¯ ¯ ¯ J−1 1        = − .  ·  . − ¯  .   . .   . .  N 1T J−1 1  . .  ¯ Rµ + (N − 1) ∆µ 1 ¯ ¯ J−1 v ¯ J−1 1 T −1 ¯ = −N ¯ Rµ T ¯ −1 v − 1 J v + Q/N Rµ J ¯ ¯ T ¯ J−1 1 1T J−1 1 N −1 T ¯ −1 ¯ 1T J−1 v + Q/N T ¯ −1 −∆ µ n 1 J ¯ v− 1 J 1 n=0 1T J−1 1 ¯ (N − 1) (N − 2) µ T ¯ −1 = N dµ − ¯ ¯ ∆ 1 J v − 1T J−1 v − Q/N 2 ¯ (N − 1) (N − 2) Q µ = N dµ − ∆ 2N

µ T 1T J−1 Rµ 1T J−1 Rν αµν = − (R ) J R + −1 ν 1T J−1 1  T  J−1 ¯ 0 0   R¯µ ··· ¯ Rν  ¯ Rµ + ∆µ 1   .. .  .  ¯ Rν + ∆ν 1     0 ¯ J−1 . .   = − .    . .  . .   . .. .. .  .  . .   . . . . R¯ µ + (N − 1) ∆µ 1 ¯ ¯ ν + (N − 1) ∆ν 1 R 0 · · · · · · J−1    ¯ ¯ 1T J−1 Rµ ¯ ¯ 1T J−1 Rν  ¯ ¯ 1T J−1 Rµ + ∆µ 1  ¯ ¯ 1T J−1 Rν + ∆ν 1      . .  . .   .  .  ¯ ¯ 1T J−1 Rµ + (N − 1) ∆µ 1 ¯ ¯ 1T J−1 Rν + (N − 1) ∆ν 1 + ¯ N 1T J−1 1 N −1 = − ¯ Rµ T ¯ ¯ ¯ J−1 Rν + n∆ν Rµ T ¯ ¯ ¯ ¯ J−1 1 + n∆µ 1J−1 Rν + n2 ∆µ ∆ν 1T J−1 1 n=0 N −1 ¯ Rµ T ¯ ¯ J−1 1 + n∆µ 1T J−1 1 N −1 ¯ Rν T ¯ ¯ J−1 1 + n ∆ν 1T J−1 1 n=0 n =0 + ¯ N 1T J−1 1 ¯ T¯ ¯ (N − 1) (N − 2) ν ¯ µ T ¯ −1 = −N Rµ J−1 Rν − ∆ R J 1 2 (N − 1) (N − 2) µ ¯ −1 ¯ ν (N − 1) (N − 2) (2N − 3) µ ν T ¯ −1 − ∆ 1J R − ∆ ∆ 1 J 1 2 6 ¯ T¯ Rµ J−1 1 ¯ T¯ Rν J−1 1 2 2 (N − 1) (N − 2) µ ν T ¯ −1 +N ¯ + ∆ ∆ 1 J 1 1T J−1 1 N (N − 1) (N − 2) ¯ T¯ ¯ T¯ + Rµ J−1 1∆ν + Rν J−1 1∆µ 2 (N − 1) (N − 2) ¯ = N αµν − ¯ N 2 − 3N − 6 ∆µ ∆ν 1T J−1 1 6N

Dipole moments have the correct translational properties because the terms and counterterms cancel perfectly. The terms in the polarizabilities cancel only to second order; the cubic terms do not cancel perfectly, giving rise to anomalous cubic scaling.

Modiﬁed dipole coupling In QEq, the external ﬁeld shifts the electronegativities on each atom by an external potential 1 E(q1 , . . . , qN ; R) = qi vi (R) + qi qj Jij (R) i 2 ij vi (R) = χi → χi − Ri · 1 E(q1 , . . . , qN ; R, ) = qi vi (R) − Ri · + qi qj Jij (R) i 2 ij

Modiﬁed dipole coupling We propose to apply the same coupling in QTPIE, which shifts the atomic voltages in a less trivial manner 1 E(q1 , . . . , qN ; R) = qi vi (R) + qi qj Jij (R) i 2 ij j (χi − χj )Sij vi = j Sij χi → χi − Ri · (χi − χj ) Sij j Ri − Rj Sij j vi (R) → vi (R, ) = − · j Sij j Sij

With this coupling, the dipole moments and polarizabilities become µ µ Sik (Ri − Rk ) −1 dµ = k J v i i l Sil P “ ” µ µ Si k R −R 1 J T −1 v i 1 J T −1 i k P Si i k l l − 1T J−1 1 µ µ k Sik (Ri − Rk ) J−1 ij l ν ν Sjl Rj − Rl αµν = − ij k Sik l Sjl P P k Sik (Ri −Rk ) (Ri −Rk ) µ µ ν ν Si k i 1 J T −1 P i 1 J T −1 k P i l Sil i l Si l + 1T J−1 1

More importantly, these expressions are still translationally invariant, but are now size-extensive dµ → dµ α→α dµ = N dµ α = Nα ¯ The key difference is that the overlap matrix attenuates any contributions across subsystems to terms of this form  ¯ −1   ¯  J 0 ··· 0 S 0 ··· 0  .. . .   .. . .   0 ¯ J−1 . .   0 S¯ . .  J−1 ν ν Sjl Rj − Rl =  .  .     . . .   ij  . .. .. .  . .. .. . l . . . .  l . . .  0 ··· ··· ¯ J −1 0 ··· ··· ¯ S ij jl     ¯ Rν ¯ Rν  ¯ Rν + ∆ν 1   ¯ Rν + ∆ν 1        ×  . .  − . .    .   .   ¯ Rν + (N − 1) ∆ν 1 ¯ Rν + (N − 1) ∆ν 1 j l = ¯ J−1 ¯ ¯ν ¯ν Sjl Rj − Rl ij l

Planar water chains with z-spacing 10.00 Angstroms 0 zz-polarizability linear fit (gradient = -0.325919, intercept = 0.007088) -100 -200 -300 Polarizability (bohr^3) -400 -500 -600 -700 -800 -900 0 500 1000 1500 2000 2500 3000 Number of atoms Planar water chains with z-spacing 10.00 Angstroms 5e+09 zz-polarizability linear fit (gradient = -5773165.409344, intercept = 1434581215.053739) 0 -5e+09 Polarizability (bohr^3) -1e+10 -1.5e+10 -2e+10 -2.5e+10 0 500 1000 1500 2000 2500 3000 Number of atoms

IV. Parameters for water

Strategy 5 2 • Aim: reproduce ab initio 1 data (LCCSD(T)/aug-cc- 4 pVTZ) with three-site 3 6 water model E = EB (R12 ) + EB (R13 ) + EB (R45 ) + EB (R46 ) OH OH OH OH • +EA (R12 , R13 , R23 ) + EA (R45 , R46 , R56 ) HOH HOH First determine QTPIE +EvdW (R14 ) OO parameters (4) by ﬁtting +EvdW (R15 ) + EvdW (R16 ) + EvdW (R42 ) + EvdW (R43 ) OH OH OH OH +EvdW (R25 ) + EvdW (R26 ) + EvdW (R35 ) + EvdW (R36 ) HH HH HH HH to dipole moments +EQT P IE (R) • 2 Then ﬁt other parameters EB (r) = k OH r − r0 OH OH 2 (10) to energies (up to HOH EA (θ) = κOH θ − θ0 HOH σAB 12 σAB 6 constant) EvdW (r) = AB AB r − r

Geometries • 1,230 water monomer geometries generated by directly varying internal coordinates • 890 water dimer geometries generated by classical MD sampling (TINKER, ﬂexible SPC water model, 1000K) • O-O distance constrained by varying vdW radius of O • 10 ps equilibration then 10 geometries at 1ps intervals

Fitting method • Weighted least-squares ﬁtting with Boltzmann weight (β = 100) f (ζ) = e−βEi (dab initio − dQTPIE (ζ))2 iν iν i ν∈{x,y,z} • Fitting done with Neader-Mead downhill simplex algorithm in scipy import scipy.optimize z0 = (8.741, 13.364, 4.528, 13.890) zopt = scipy.optimize.fmin(f, z0) print quot;Optimal parameters = quot;, zopt

Best-ﬁt parameters Par./eV QEq New χ(H) 4.5280 5.6841 η(H) 13.8904 12.4166 χ(O) 8.741 7.9173 η(O) 13.364 13.1643

QTPIE Dipoles correlation ab initio

Dipole moment per water / D Number of water molecules

yy-polarizability per water / Å³ Number of water molecules

zz-polarizability per water / Å³ Number of water molecules

That’s all, folks!

## Add a comment

 User name: Comment:

## Related pages

### Max Steel s01e07 Hard Water part 2 - YouTube

Max Steel s01e07 Hard Water part 2 Nathan Kima. Subscribe Subscribed Unsubscribe 433 433. ... +YouTube; Terms; Privacy; Policy & Safety Send ...

### Sylvan – The Waters I Traveled, Part II Lyrics | Genius Lyrics

The Waters I Traveled, Part II Lyrics. Awake on a salving shore ... The surges have ceased from me It seems that a saving stream somehow has brought me ...

### Blue Water High Deutsch Liebeskummer Part2 - YouTube

Blue Water High Deutsch Liebeskummer Part2 ... Blue Water High Deutsch Folge Haialarm Part 1/2 ... Blue Water High Gewissensbisse Part1/2 ...

### Home : Waters

ACQUITY UPC 2; ACQUITY Arc UHPLC; Alliance HPLC; ... Using Quality Parts Discover how Waters ensures ... Waters has developed innovative analytical science ...

### Australian Television: Blue Water High: episode guide ...

... nobody will ever forget their time at Blue Water High. Eric and Mike face off ... gets more screen time after appearing previously in episode 2.09 as a ...

### PHASE 2) PREPARATION OF PART FOR PAINTING & PAINT 1. SAND ...

PHASE 2) PREPARATION OF PART FOR PAINTING & PAINT 1. ... WASH paintable surfaces and the tape ledge w/ 3M #7447 Scotchbrite pad, TSP, and warm water 3.

### SoCal Water\$mart – Residential Rebates

... SoCal Water\$mart; Rebates. ... be eligible for \$2 per ... water districts that provides drinking water to nearly 19 million people in parts of ...

### ASME Code - druckgeraete-online.de

(ASME Code) Der ASME Boiler ... Er erscheint alle 2 Jahre ... Superheaters, economizers, and other pressure parts connected to the boiler without ...