#### Previous topic

Least-Squares Fitting

#### Next topic

Minimal Surface Problem

# Hessian of Two Particle Coulomb PotentialΒΆ

This is an example from the Scipy-User mailing list. The original discussion can be found on http://www.mail-archive.com/numpy-discussion@scipy.org/msg28633.html

# import a tool to use / as a symbol for normal division

#import system data
import sys, os

import scipy as sp
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# and subpackages
from scipy import *
from scipy import linalg, optimize, constants

import numpy
import algopy
import time

#-----------------------------------------------------------------------------------------
#       Hamiltonian      H=sum_i(p_i^2/(2m)+ 1/2 * m * w^2 x_i^2)+ sum_(i!=j)(a/|x_i-x_j|)
#-----------------------------------------------------------------------------------------

class classicalHamiltonian:
def __init__(self):

self.N = 2                                                      #N is a scalar, it's the number of ions in the chain
f = 1000000                                                     #f is a scalar, it's the trap frequency
self.w = 2*pi*f                                                 #w is a scalar, it's the angular velocity corresponding to the trap frequency
self.C = (4*pi*constants.epsilon_0)**(-1)*constants.e**2        #C is a scalar, it's the Coulomb constant times the electronic charge in SI
self.m = 39.96*1.66e-27                                         #m is the mass of a single trapped ion in the chain

def potential(self, positionvector):                                    #Defines the potential that is going to be minimized

x= positionvector                                               #x is an 1-d array (vector) of lenght N that contains the positions of the N ions
w= self.w
C= self.C
m= self.m

#First we consider the potential of the harmonic osszilator
Vx = 0.5 * m * (w**2) * sum(x**2)

for i in range(len(x)):
for j in range(i+1, len(x)):
Vx += C * (abs(x[i]-x[j]))**(-1)        #then we add the coulomb interaction

return Vx

def initialposition(self):              #Defines the initial position as an estimate for the minimize process

N= self.N
x_0 = r_[-(N-1)/2:(N-1)/2:N*1j]
return x_0

def normal_modes(self, eigenvalues):    #the computed eigenvalues of the matrix Vx are of the form (normal_modes)^2*m.
m = self.m
normal_modes = sqrt(eigenvalues/m)
return normal_modes

c=classicalHamiltonian()
xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10)

x = algopy.UTPM.init_hessian(xopt)
y = c.potential(x)
hessian = algopy.UTPM.extract_hessian(2, y)

print(hessian)


As output one obtains:

\$ python hessian_of_potential_function.py
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 81
Function evaluations: 153
[[  5.23748399e-12  -2.61873843e-12]
[ -2.61873843e-12   5.23748399e-12]]