# Constrained ESP fitting

## using Horton's cost function

According to https://theochem.github.io/horton/2.1.0/user_postproc_espfit.html, cost function is constructed as

$$ c(\mathbf{q}, \Delta V_\text{ref}) = \int_V d\mathbf{r} \omega(\mathbf{r}) \cdot \left( V_\text{ab initio}(\mathbf{r}) - \sum_{i=1}^N \frac{q_i}{\mathbf{r} - \mathbf{R}_i} - \Delta V_\text{ref} \right)^2  $$

> $\Delta V_\text{ref}$: constant that may account for differences in reference between the ab initio ESP and the point charge ESP. The need for such a term was established in the REPEAT paper for ESP fitting in 3D periodic systems.

We look at an aperiodic system, neglect this term for now. The (unconstrained) cost function takes the general quadratic form

$$ c_u(\mathbf{x}) = \mathbf{x}^T A\ \mathbf{x} - 2 \mathbf{b}^T \mathbf{x} - C $$

with 

$$ C = - \int_V d\mathbf{r} \omega(\mathbf{r}) \cdot \left[ V_\text{ab initio}(\mathbf{r}) \right]^2$$

entry i,j of matrix $A^{N\times N}$

$$  A_{ij} = \int_V d\mathbf{r}\ \omega(\mathbf{r}) \cdot \left( \frac{1}{|\mathbf{r} - \mathbf{R}_i|} \cdot \frac{1}{|\mathbf{r} - \mathbf{R}_j|} \right)  $$

and entry i of vector $\mathbf{b}$

$$  b_i = \int_V d\mathbf{r}\ \omega(\mathbf{r}) \cdot \frac{V_\text{ab initio}(\mathbf{r})}{|\mathbf{r} - \mathbf{R}_i|} $$

In the code below, first the miniumum of an unconstrained system is found by solving

$$ \frac{1}{2} \cdot \frac{\mathrm{d} c(\mathbf{x})}{\mathrm d  \mathbf{x}}  = \frac{1}{2} \cdot \nabla_\mathbf{x} c = A \mathbf{x} - \mathbf{b} = 0$$

with $\nabla (\mathbf{b}^T \mathbf{x}) = \mathbf{b}$ and 
$\nabla \cdot (\mathbf{x}^T A \mathbf{x}) = (A + A^T) \mathbf{x} 
= 2 A \mathbf{x} $ for symmetric A, as in our case. 
The (unconstrained) solution

$$ \mathbf{x}_u = A^{-1} \mathbf{b} $$

is corrected for *one* total charge constraint of the form 

$$ g(\mathbf{x}) = \mathbf{d} \cdot \mathbf{x} - q_\text{tot} = 0 $$

with all entries of $\mathbf{d}$ unity. Notice that in the code below, the whole system is normalized in order to have unity diagonal entries $A_{jj}$. We neglect this normalization here.

A Lagrange multiplier $\lambda$ is introduced into the *constrained* cost function 

$$ c_c(\mathbf{x},\lambda) = \mathbf{x}^T A\ \mathbf{x} - 2 \mathbf{b}^T \mathbf{x} - C + \lambda \cdot  g(\mathbf{x}) $$

and the system

$$ \frac{1}{2} \cdot \nabla_\mathbf{x} c = A \mathbf{x} - \mathbf{b} + \lambda \cdot \nabla_\mathbf{x} g(\mathbf{x}) = A \mathbf{x} - \mathbf{b} +\lambda \cdot \mathbf{d} = 0 $$

$$  \nabla_\mathbf{\lambda} c = g(\mathbf{x}) = \mathbf{d} \cdot \mathbf{x} - q_\text{tot} = 0 $$

is solved by finding a correction for the unconstrained solution

$$ \mathbf{x} = \mathbf{x}_u - \lambda \cdot\delta \mathbf{x} = A^{-1} \mathbf{b} - \lambda \cdot \delta \mathbf{x} $$

as 

$$ - \lambda A \delta \mathbf{x} + \lambda \cdot \mathbf{d} = 0 
\Leftrightarrow \delta \mathbf{x} = A^{-1} \mathbf{d}$$

and the Lagrange multiplier by

$$ g(\mathbf{x}) 
    = \mathbf{d} \cdot \left( \mathbf{x}_u 
        - \lambda \ \delta \mathbf{x} \right)- q_\text{tot} 
    = \mathbf{d} \cdot A^{-1} \mathbf{b} 
        - \lambda \ \mathbf{d} \cdot \delta \mathbf{x} - q_\text{tot} 
    = 0 $$

$$ \lambda = \frac{\mathbf{d} \cdot A^{-1} \mathbf{b} - q_\text{tot}}
    { \mathbf{d} \cdot \delta \mathbf{x} }
    = \frac{\mathbf{b} \cdot \delta \mathbf{x} - q_\text{tot}}
    { \mathbf{d} \cdot \delta \mathbf{x} }$$
    
    
and thereby the constrained minimum at
$$ \mathbf{x} 
    = \mathbf{x}_u 
        - \frac{\mathbf{b} \cdot \delta \mathbf{x} - q_\text{tot}}
        { \mathbf{d} \cdot \delta \mathbf{x} } 
        \cdot \delta \mathbf{x}$$
        
as implemented in HORTON


```python
# modified excerpt from horton/espfit/cost.py

import numpy as np
import h5py

def solve(_A, _B, _C, _natom, qtot=None, ridge=0.0):
    # apply regularization to atomic degrees of freedom
    A = _A.copy()
    A.ravel()[::len(A)+1][:_natom] += ridge*np.diag(A)[:_natom].mean()
    # construct preconditioned matrices
    norms = np.diag(A)**0.5
    print("Diagonal of A: \n {}".format(np.diag(A)))
    print("Vector norms: \n {}".format(norms))
    print("Matrix norms: \n {}".format(norms*norms.reshape(-1,1)))

    A = A/norms/norms.reshape(-1,1)
    print("Normed diagonal of A: \n {}".format(np.diag(A)))

    B = _B/norms
    print("Normed B: \n {}".format(B))

    x = np.linalg.solve(A, B)
    if qtot is not None:
        # Fix the total charge with a lagrange multiplier
        d = np.zeros(len(A))
        d[:_natom] = 1/norms[:_natom]
        d[_natom:] = 0.0
        aid = np.linalg.solve(A, d)
        lagrange = (np.dot(aid, B) - qtot)/np.dot(aid, d)
        x -= aid*lagrange
    x /= norms
    return x
```


```python
np.set_printoptions(precision=4)
```


```python
# Example: water
# load the cost function
cost_function = h5py.File('h2o.cost.h5')
cost = cost_function['cost']
A = cost['A'][:]
B = cost['B'][:]
C = cost['C'].value
N = cost['natom'].value
print("A: \n {}".format(A))
print("B: {}".format(B))
print("C: {}".format(C))
print("N: {}".format(N))
```

    A: 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B: [-0.0105 -0.0797 -0.0797]
    C: 0.0266115260155
    N: 3



```python
q1 = solve(A,B,C,N)
```

    Diagonal of A: 
     [ 8.5387  9.0212  9.0212]
    Vector norms: 
     [ 2.9221  3.0035  3.0035]
    Matrix norms: 
     [[ 8.5387  8.7766  8.7766]
     [ 8.7766  9.0212  9.0212]
     [ 8.7766  9.0212  9.0212]]
    Normed diagonal of A: 
     [ 1.  1.  1.]
    Normed B: 
     [-0.0036 -0.0265 -0.0265]



```python
q1
```




    array([ 0.3378, -0.1699, -0.1699])




```python
sum(q1)
```




    -0.0019095241274412478




```python
q2 = solve(A,B,C,N,qtot=0)
```

    Diagonal of A: 
     [ 8.5387  9.0212  9.0212]
    Vector norms: 
     [ 2.9221  3.0035  3.0035]
    Matrix norms: 
     [[ 8.5387  8.7766  8.7766]
     [ 8.7766  9.0212  9.0212]
     [ 8.7766  9.0212  9.0212]]
    Normed diagonal of A: 
     [ 1.  1.  1.]
    Normed B: 
     [-0.0036 -0.0265 -0.0265]



```python
q2
```




    array([ 0.3396, -0.1698, -0.1698])




```python
sum(q2)
```




    0.0




```python
q3 = solve(A,B,C,N,qtot=1)
```

    Diagonal of A: 
     [ 8.5387  9.0212  9.0212]
    Vector norms: 
     [ 2.9221  3.0035  3.0035]
    Matrix norms: 
     [[ 8.5387  8.7766  8.7766]
     [ 8.7766  9.0212  9.0212]
     [ 8.7766  9.0212  9.0212]]
    Normed diagonal of A: 
     [ 1.  1.  1.]
    Normed B: 
     [-0.0036 -0.0265 -0.0265]



```python
q3
```




    array([ 1.2558, -0.1279, -0.1279])




```python
sum(q3)
```




    0.99999999999999989



# Arbitrary number of equality constraints

We modifiy the optimization in order to allow for an arbitrary amount of constraints.

M lagrange multipliers $\lambda_k$ are introduced into the *constrained* cost function 

$$ c_c(\mathbf{x},\mathbf{\lambda}) = \mathbf{x}^T A\ \mathbf{x} - 2 \mathbf{b}^T \mathbf{x} - C + \sum_{k=1}^M \lambda_k \cdot  g_k(\mathbf{x}) $$

All our constraints (charge groups and symmetries) will be of linear form 

$$ g_k(\mathbf{x}) = \mathbf{d}_k \cdot \mathbf{x} - q_k = 0 $$

and thus can be compacted into matrix form

$$ D^{(M \times N)} \mathbf{x} - \mathbf{q}^{(M\times 1)} = 0 $$

with 

$$ D^T = [\mathbf{d}_1, \mathbf{d}_2, \dots , \mathbf{d}_M] $$

and hence 
$$ c_c(\mathbf{x}^{(N\times1)},\mathbf{\lambda}^{(M\times1)}) = \mathbf{x}^T A\ \mathbf{x} 
    - 2 \mathbf{b}^T \mathbf{x} 
    - C + \mathbf{\lambda}^T \cdot \left( D \mathbf{x} - \mathbf{q} \right) $$

Derivative 

$$ 
\begin{align}
    \nabla_\mathbf{x} \cdot c_c & = 2\ A\ \mathbf{x} 
    + \sum_{k=1}^M \lambda_k \mathbf{d}_k - 2 \mathbf{b} & = 0\\
    \nabla_\mathbf{\lambda} \cdot c_c & = D\ \mathbf{x} - \mathbf{q} & = 0
\end{align}
$$

Identify 

$$ D^T \mathbf{\lambda}= \sum_{k=1}^M \lambda_k \mathbf{d}_k  $$

and solve

$$ \tilde{A} \mathbf{y} - \tilde{\mathbf{b}} = 0$$ 

with generalized $\mathbf{y}^T = [\mathbf{x}^T, \mathbf{\lambda}^T ]$ 
as well as $(N+M)\times(N+M)$ matrix

$$ \tilde{A} = 
    \begin{bmatrix}
        2 A & D^T \\
        D & 0
    \end{bmatrix} $$
    
and $(N+M)$ vector

$$ \tilde{\mathbf{b}} =  
   \begin{bmatrix}
        2 \mathbf{b} \\
        \mathbf{q}
   \end{bmatrix} $$


```python
import numpy as np

def constrainedMinimize(A_matrix, b_vector, C_scalar, D_matrix = None, q_vector = np.array([0]), debug=False):
    N = b_vector.shape[0]
    M = q_vector.shape[0]
    
    if not isinstance(D_matrix,np.ndarray):
        D_matrix = np.atleast_2d( np.ones(N) )
    
    if debug:
        print("{:d} unknowns, {:d} equality constraints".format(N,M))
        print("A {}: \n {}".format(A_matrix.shape,A_matrix))
        print("B {}: \n {}".format(b_vector.shape,b_vector))
        print("C {}: \n {}".format(C_scalar.shape,C_scalar))
        print("D {}: \n {}".format(D_matrix.shape,D_matrix))
        print("q {}: \n {}".format(q_vector.shape,q_vector))

    A_upper = np.block([ 2*np.atleast_2d(A_matrix), np.atleast_2d(D_matrix).T ])
    A_lower = np.block([ np.atleast_2d(D_matrix), np.atleast_2d(np.zeros((M,M)))])
    if debug:
        print("upper A ({}): \n {}".format(A_upper.shape,A_upper))
        print("upper A ({}): \n {}".format(A_lower.shape,A_lower))

    A = np.block( [ [ A_upper ], [ A_lower ] ] )
    
    if debug:
        print("block A ({}): \n {}".format(A.shape,A))
    
    B = np.block( [2*np.atleast_1d(b_vector), np.atleast_1d(q_vector)] )
    
    if debug:
        print("block B ({}): \n {}".format(B.shape,B))

    C = C_scalar
    
    x = np.linalg.solve(A, B)
    
    return x
```


```python
# charges: list of charges to meet symmetry constraint, 
#   indices beginning from 0
#   can also be list of list of group of charges to be symmetric, 
#   indices beginning from 0, i.e. charges = [[0,1],[2,3,4]] requires
#   the cummulative charge of atom 0 and 1 to equal 
#   the sum over charges 2 to 4
# N: total number of charges
# symmetry: scalar or list of symmetry requirements
#   default 1 causes equality, -1 antisymmetry, 
#   values different from +/- unity result in relative scaling
def constructPairwiseSymmetryConstraints(charges, N, symmetry=1.0, debug=False):
    M = len(charges)-1

    symmetry = symmetry*np.ones(M)
    
    if debug:
        print("{:d} unknowns, {:d} pairwise equality constraints".format(N,M))
        print("symmetry list ({}):\n{}".format(symmetry.shape,symmetry))
        
    D = np.atleast_2d(np.zeros((M,N)))
    q = np.atleast_1d(np.zeros(M))
    D[:,charges[0]] = 1
    
    for j in range(M):
        D[j,charges[j+1]] = -1.0*symmetry[j]

    if debug:
        print("D ({}):\n{}".format(D.shape,D))
        print("q ({}):\n{}".format(q.shape,q))
    
    return D, q
```


```python
# chargeGroups: list of list of charges in charge groups, 
#    indices beginning from 0, i.e. chargeGroups = [[0,1],[2,3,4]]
# N: total number of charges
# q: required charges per charge group

def constructChargegroupConstraints(chargeGroups, N, q=0, debug=False):
    M = len(chargeGroups)

    q = np.atleast_1d(q*np.ones(M))

    if debug:
        print("{:d} unknowns, {:d} pairwise equality constraints".format(N,M))

    D = np.atleast_2d(np.zeros((M,N)))
    #q = np.atleast_2d(np.zeros(M))    
    
    for j in range(M):
        D[j,chargeGroups[j]] = 1.0

    if debug:
        print("D ({}):\n{}".format(D.shape,D))
        print("q ({}):\n{}".format(q.shape,q))
    
    return D, q
```


```python
q1u = solve(A,B,C,N)
```

    Diagonal of A: 
     [ 8.5387  9.0212  9.0212]
    Vector norms: 
     [ 2.9221  3.0035  3.0035]
    Matrix norms: 
     [[ 8.5387  8.7766  8.7766]
     [ 8.7766  9.0212  9.0212]
     [ 8.7766  9.0212  9.0212]]
    Normed diagonal of A: 
     [ 1.  1.  1.]
    Normed B: 
     [-0.0036 -0.0265 -0.0265]



```python
q1c = constrainedMinimize(A,B,C,debug=True)
```

    3 unknowns, 1 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (1, 3): 
     [[ 1.  1.  1.]]
    q (1,): 
     [0]
    upper A ((3, 4)): 
     [[ 17.0774  17.0433  17.0433   1.    ]
     [ 17.0433  18.0423  16.7912   1.    ]
     [ 17.0433  16.7912  18.0423   1.    ]]
    upper A ((1, 4)): 
     [[ 1.  1.  1.  0.]]
    block A ((4, 4)): 
     [[ 17.0774  17.0433  17.0433   1.    ]
     [ 17.0433  18.0423  16.7912   1.    ]
     [ 17.0433  16.7912  18.0423   1.    ]
     [  1.       1.       1.       0.    ]]
    block B ((4,)): 
     [-0.021  -0.1594 -0.1594  0.    ]



```python
q1u
```




    array([ 0.3378, -0.1699, -0.1699])




```python
q1c
```




    array([ 0.3396, -0.1698, -0.1698, -0.0326])




```python
q2_horton = solve(A,B,C,N,qtot=1)
```

    Diagonal of A: 
     [ 8.5387  9.0212  9.0212]
    Vector norms: 
     [ 2.9221  3.0035  3.0035]
    Matrix norms: 
     [[ 8.5387  8.7766  8.7766]
     [ 8.7766  9.0212  9.0212]
     [ 8.7766  9.0212  9.0212]]
    Normed diagonal of A: 
     [ 1.  1.  1.]
    Normed B: 
     [-0.0036 -0.0265 -0.0265]



```python
q2c = constrainedMinimize(A,B,C,q_vector=np.array([1]),debug=True)
```

    3 unknowns, 1 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (1, 3): 
     [[ 1.  1.  1.]]
    q (1,): 
     [1]
    upper A ((3, 4)): 
     [[ 17.0774  17.0433  17.0433   1.    ]
     [ 17.0433  18.0423  16.7912   1.    ]
     [ 17.0433  16.7912  18.0423   1.    ]]
    upper A ((1, 4)): 
     [[ 1.  1.  1.  0.]]
    block A ((4, 4)): 
     [[ 17.0774  17.0433  17.0433   1.    ]
     [ 17.0433  18.0423  16.7912   1.    ]
     [ 17.0433  16.7912  18.0423   1.    ]
     [  1.       1.       1.       0.    ]]
    block B ((4,)): 
     [-0.021  -0.1594 -0.1594  1.    ]



```python
q2_horton
```




    array([ 1.2558, -0.1279, -0.1279])




```python
q2c
```




    array([  1.2558,  -0.1279,  -0.1279, -17.1071])




```python
d1 = np.array([0,1,-2])
```


```python
q3c = constrainedMinimize(A,B,C,D_matrix=d1,q_vector=np.array([0]),debug=True)
```

    3 unknowns, 1 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (3,): 
     [ 0  1 -2]
    q (1,): 
     [0]
    upper A ((3, 4)): 
     [[ 17.0774  17.0433  17.0433   0.    ]
     [ 17.0433  18.0423  16.7912   1.    ]
     [ 17.0433  16.7912  18.0423  -2.    ]]
    upper A ((1, 4)): 
     [[ 0.  1. -2.  0.]]
    block A ((4, 4)): 
     [[ 17.0774  17.0433  17.0433   0.    ]
     [ 17.0433  18.0423  16.7912   1.    ]
     [ 17.0433  16.7912  18.0423  -2.    ]
     [  0.       1.      -2.       0.    ]]
    block B ((4,)): 
     [-0.021  -0.1594 -0.1594  0.    ]



```python
q3c
```




    array([ 0.2884, -0.1935, -0.0967,  0.0403])




```python
d1 = np.array([0,1,-2])
```


```python
d2 = np.array([1,1,1])
```


```python
D = np.block([[d1],[d2]])
```


```python
D
```




    array([[ 0,  1, -2],
           [ 1,  1,  1]])




```python
q = np.array([0,3])
```


```python
q4c = constrainedMinimize(A,B,C,D_matrix=D,q_vector=q,debug=True)
```

    3 unknowns, 2 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (2, 3): 
     [[ 0  1 -2]
     [ 1  1  1]]
    q (2,): 
     [0 3]
    upper A ((3, 5)): 
     [[ 17.0774  17.0433  17.0433   0.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.    ]
     [ 17.0433  16.7912  18.0423  -2.       1.    ]]
    upper A ((2, 5)): 
     [[ 0.  1. -2.  0.  0.]
     [ 1.  1.  1.  0.  0.]]
    block A ((5, 5)): 
     [[ 17.0774  17.0433  17.0433   0.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.    ]
     [ 17.0433  16.7912  18.0423  -2.       1.    ]
     [  0.       1.      -2.       0.       0.    ]
     [  1.       1.       1.       0.       0.    ]]
    block B ((5,)): 
     [-0.021  -0.1594 -0.1594  0.      3.    ]



```python
q4c
```




    array([  3.0755e+00,  -5.0328e-02,  -2.5164e-02,   1.0487e-02,  -5.1256e+01])




```python
sum(q4c)
```




    -48.245273601712007




```python
D2,q2 = constructPairwiseSymmetryConstraints([1,3,8],10,debug=True)
```

    10 unknowns, 2 pairwise equality constraints
    symmetry list ((2,)):
    [ 1.  1.]
    D ((2, 10)):
    [[ 0.  1.  0. -1.  0.  0.  0.  0.  0.  0.]
     [ 0.  1.  0.  0.  0.  0.  0.  0. -1.  0.]]
    q ((2,)):
    [ 0.  0.]



```python
D3,q3 = constructPairwiseSymmetryConstraints([1,3,8],10,symmetry=-1,debug=True)
```

    10 unknowns, 2 pairwise equality constraints
    symmetry list ((2,)):
    [-1. -1.]
    D ((2, 10)):
    [[ 0.  1.  0.  1.  0.  0.  0.  0.  0.  0.]
     [ 0.  1.  0.  0.  0.  0.  0.  0.  1.  0.]]
    q ((2,)):
    [ 0.  0.]



```python
D4,q4 = constructPairwiseSymmetryConstraints([1,3,8],10,symmetry=np.array([1,-1]),debug=True)
```

    10 unknowns, 2 pairwise equality constraints
    symmetry list ((2,)):
    [ 1. -1.]
    D ((2, 10)):
    [[ 0.  1.  0. -1.  0.  0.  0.  0.  0.  0.]
     [ 0.  1.  0.  0.  0.  0.  0.  0.  1.  0.]]
    q ((2,)):
    [ 0.  0.]



```python
D5,q5 = constructPairwiseSymmetryConstraints([[1,2],[3,4,5],[8,9]],10,symmetry=np.array([1,-1]),debug=True)
```

    10 unknowns, 2 pairwise equality constraints
    symmetry list ((2,)):
    [ 1. -1.]
    D ((2, 10)):
    [[ 0.  1.  1. -1. -1. -1.  0.  0.  0.  0.]
     [ 0.  1.  1.  0.  0.  0.  0.  0.  1.  1.]]
    q ((2,)):
    [ 0.  0.]


### full example
Water with antisymmetric hydrogens and +1 integer charge group of oxygen and 1 hydrogen


```python
# Example: water
import h5py
#load the cost function
cost_function = h5py.File('h2o.cost.h5')
cost = cost_function['cost']
A = cost['A'][:]
B = cost['B'][:]
C = cost['C'].value
N = cost['natom'].value
print("A: {}".format(A))
print("B: {}".format(B))
print("C: {}".format(C))
print("N: {}".format(N))


```

    A: [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B: [-0.0105 -0.0797 -0.0797]
    C: 0.0266115260155
    N: 3



```python
D1,q1 = constructPairwiseSymmetryConstraints(
    charges=[1,2],N=3,symmetry=-1,debug=True)
```

    3 unknowns, 1 pairwise equality constraints
    symmetry list ((1,)):
    [-1.]
    D ((1, 3)):
    [[ 0.  1.  1.]]
    q ((1,)):
    [ 0.]



```python
D1, q1
```




    (array([[ 0.,  1.,  1.]]), array([ 0.]))




```python
D2,q2 = constructChargegroupConstraints(
    chargeGroups=[[0,1]],N=3,q=1,debug=True)
```

    3 unknowns, 1 pairwise equality constraints
    D ((1, 3)):
    [[ 1.  1.  0.]]
    q ((1,)):
    [ 1.]



```python
D2, q2
```




    (array([[ 1.,  1.,  0.]]), array([ 1.]))




```python
D = np.vstack([D1,D2])
```


```python
D
```




    array([[ 0.,  1.,  1.],
           [ 1.,  1.,  0.]])




```python
q = np.hstack([q1,q2])
```


```python
q
```




    array([ 0.,  1.])




```python
X = constrainedMinimize(A,B,C,D,q,debug=True)
```

    3 unknowns, 2 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (2, 3): 
     [[ 0.  1.  1.]
     [ 1.  1.  0.]]
    q (2,): 
     [ 0.  1.]
    upper A ((3, 5)): 
     [[ 17.0774  17.0433  17.0433   0.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.    ]
     [ 17.0433  16.7912  18.0423   1.       0.    ]]
    upper A ((2, 5)): 
     [[ 0.  1.  1.  0.  0.]
     [ 1.  1.  0.  0.  0.]]
    block A ((5, 5)): 
     [[ 17.0774  17.0433  17.0433   0.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.    ]
     [ 17.0433  16.7912  18.0423   1.       0.    ]
     [  0.       1.       1.       0.       0.    ]
     [  1.       1.       0.       0.       0.    ]]
    block B ((5,)): 
     [-0.021  -0.1594 -0.1594  0.      1.    ]



```python
X
```




    array([ 0.1267,  0.8733, -0.8733, -1.2267, -2.1852])




```python
Q = X[:N]
```


```python
Lagrange = X[N:]
```


```python
Q # all constraints fulfilled
```




    array([ 0.1267,  0.8733, -0.8733])




```python
Q[0]+Q[1]
```




    1.0




```python
Q[1]+Q[2]
```




    1.1102230246251565e-16




```python
Lagrange
```




    array([-1.2267, -2.1852])




```python

```


```python
D3,q3 = constructPairwiseSymmetryConstraints(
    charges=[1,2],N=3,symmetry=1,debug=True)
```

    3 unknowns, 1 pairwise equality constraints
    symmetry list ((1,)):
    [ 1.]
    D ((1, 3)):
    [[ 0.  1. -1.]]
    q ((1,)):
    [ 0.]



```python
D4,q4 = constructPairwiseSymmetryConstraints(
    charges=[1,2],N=3,symmetry=-1,debug=True)
```

    3 unknowns, 1 pairwise equality constraints
    symmetry list ((1,)):
    [-1.]
    D ((1, 3)):
    [[ 0.  1.  1.]]
    q ((1,)):
    [ 0.]



```python
D5,q5 = constructChargegroupConstraints(chargeGroups=[[0,1,2],[0,1,2]],N=3,q=[0,1],debug=True)
```

    3 unknowns, 2 pairwise equality constraints
    D ((2, 3)):
    [[ 1.  1.  1.]
     [ 1.  1.  1.]]
    q ((2,)):
    [ 0.  1.]



```python
D = np.vstack([D3,D4,D5])
```


```python
D
```




    array([[ 0.,  1., -1.],
           [ 0.,  1.,  1.],
           [ 1.,  1.,  1.],
           [ 1.,  1.,  1.]])




```python
q = np.hstack([q3,q4,q5])
```


```python
q
```




    array([ 0.,  0.,  0.,  1.])




```python
X = constrainedMinimize(A,B,C,D,q,debug=True)
```

    3 unknowns, 4 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (4, 3): 
     [[ 0.  1. -1.]
     [ 0.  1.  1.]
     [ 1.  1.  1.]
     [ 1.  1.  1.]]
    q (4,): 
     [ 0.  0.  0.  1.]
    upper A ((3, 7)): 
     [[ 17.0774  17.0433  17.0433   0.       0.       1.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.       1.       1.    ]
     [ 17.0433  16.7912  18.0423  -1.       1.       1.       1.    ]]
    upper A ((4, 7)): 
     [[ 0.  1. -1.  0.  0.  0.  0.]
     [ 0.  1.  1.  0.  0.  0.  0.]
     [ 1.  1.  1.  0.  0.  0.  0.]
     [ 1.  1.  1.  0.  0.  0.  0.]]
    block A ((7, 7)): 
     [[ 17.0774  17.0433  17.0433   0.       0.       1.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.       1.       1.    ]
     [ 17.0433  16.7912  18.0423  -1.       1.       1.       1.    ]
     [  0.       1.      -1.       0.       0.       0.       0.    ]
     [  0.       1.       1.       0.       0.       0.       0.    ]
     [  1.       1.       1.       0.       0.       0.       0.    ]
     [  1.       1.       1.       0.       0.       0.       0.    ]]
    block B ((7,)): 
     [-0.021  -0.1594 -0.1594  0.      0.      0.      1.    ]



```python
X[:N]
```




    array([-0.004 ,  0.0258, -0.0218])




```python
X[N:]
```




    array([ -3.1395e-02,  -1.2500e-01,   1.4412e+17,  -1.4412e+17])




```python
X = constrainedMinimize(A,B,C,D5,q5,debug=True)
```

    3 unknowns, 2 equality constraints
    A (3, 3): 
     [[ 8.5387  8.5216  8.5216]
     [ 8.5216  9.0212  8.3956]
     [ 8.5216  8.3956  9.0212]]
    B (3,): 
     [-0.0105 -0.0797 -0.0797]
    C (): 
     0.0266115260155
    D (2, 3): 
     [[ 1.  1.  1.]
     [ 1.  1.  1.]]
    q (2,): 
     [ 0.  1.]
    upper A ((3, 5)): 
     [[ 17.0774  17.0433  17.0433   1.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.    ]
     [ 17.0433  16.7912  18.0423   1.       1.    ]]
    upper A ((2, 5)): 
     [[ 1.  1.  1.  0.  0.]
     [ 1.  1.  1.  0.  0.]]
    block A ((5, 5)): 
     [[ 17.0774  17.0433  17.0433   1.       1.    ]
     [ 17.0433  18.0423  16.7912   1.       1.    ]
     [ 17.0433  16.7912  18.0423   1.       1.    ]
     [  1.       1.       1.       0.       0.    ]
     [  1.       1.       1.       0.       0.    ]]
    block B ((5,)): 
     [-0.021  -0.1594 -0.1594  0.      1.    ]



```python
X
```




    array([  3.5032e-01,  -1.6107e-01,  -1.8995e-01,   1.4412e+17,  -1.4412e+17])




```python

```

c741ea3893b3f874d783dca54f5e73915e0629c4
ipynb
Jupyter Notebook
examples/horton-esp-fit-constrained.ipynb
jotelha/smampppp
729e4733b436e68adfe07bcaa39a47727d0c8dd8
MIT
2022-03-15T17:23:52.000Z
2022-03-15T17:23:52.000Z
examples/horton-esp-fit-constrained.ipynb
jotelha/smampppp
729e4733b436e68adfe07bcaa39a47727d0c8dd8
MIT
examples/horton-esp-fit-constrained.ipynb
jotelha/smampppp
729e4733b436e68adfe07bcaa39a47727d0c8dd8
MIT
2018-04-06T11:29:44.000Z
2018-04-06T11:29:44.000Z
Qwen/Qwen-72B
1. YES
2. YES