[MMTK] Orthogonality of normal modes

ra.biehl ra.biehl at fz-juelich.de
Tue Apr 17 13:51:20 UTC 2012


Hi

sorry but this is not ok
The mistake in your writing are these two:
>e'=1/sqrt(M)e   this is wrong as shown below
in
>m1cart = np.ravel((m1.array) * np.sqrt(m[:, np.newaxis]))
m1 is already in cartesian coordinates (unweighted)
for Energetic Modes you are right (weight =1)
for Vibrational modes not


I try to explain:
I use ' for time derivative  (see text books or one of Konrads papers)

M*r''=-K*r     (1)       as the equation of motion with force and mass
                         matrix K and M

r= A*exp(-i*w*t)      is assumed as a solution to transform the double
                      time derivative to
r''=-w**2*A*exp(-i*w*t)=-w**2*r

(1) ->  M*w^2*r=K*r

This look like an eigenequation    K^*e^=lam*e^   with eigenvalue lam
and eigenvector e^
We get this if we divide by M, but this is wrong because we deal with
vectors and matrices and not with scalars.

So we use
a :   M=sqrt(M)*sqrt(M) and
b:     1=sqrt(M)**-1 * sqrt(M)
in (1)
->     a                              b
  sqrt(M)*sqrt(M) * w^2*r=K* sqrt(M)**-1*sqrt(M)  *r    (2)

now we multiply both sides with sqrt(M)**-1 from left side
and take w**2 in front

  w**2 * sqrt(M) *r=sqrt(M)**-1*K* sqrt(M)**-1 * sqrt(M)  *r

with the new mass weighted coordinates
  r^ = sqrt(M) *r
and mass weighted  force matrix
  K^ = sqrt(M)**-1*K* sqrt(M)**-1

  w**2 * r^=K^ * r^        (3)    from (2) is reached

This is solved and gives the mass weighted eigenvectors r^_i for
eigenvalue w**2_i
which are orthogonal
  r^_i * r^_j=delta_ij = 1 for i=j ; else 0
in the normal dotProduct *

r^_i * r^_j=delta_ij     transformed to unweighted coordinates is

        ->   a
r_i * sqrt(M)*sqrt(M) * r_j=delta_ij

r_i _M*_ r_j=delta_ij      (4)  with _M*_=*M* as massWeightedDotProduct
                                for unweighted coordinates

this does not conclude  r_i * r_j=delta_ij

Now the modes of MMTK are the unweighted space coordinates r_i of a mode
in amplitude determined by the temperature. For T=0 normalized to 1 as
Konrad pointed out.
This is also true for Frictional modes with a different weight or
Energetic with weight=1.
Orthogonal eigenvectors can be accessed as modes.rawMode(i) and tested
for orthogonality

To test for orthogonality you can use also as in (4)
r_i.massWeightedDotProduct(r_j)
or as Konrad wrote it modes[i].massWeightedDotProduct(modes[j])

Perhaps this helps

Best Ralf



Viele Grüsse
and best wishes

Ralf Biehl


--

Dr. Ralf Biehl
Jülich Centre for Neutron Science (JCNS-1) &
Institute for Complex Systems (ICS-1)
Forschungszentrum Jülich
52425 Jülich
Germany

Tel.:   +49-(0)2461-61-4685
Fax:    +49-(0)2461-61-2610
Email:  ra.biehl at fz-juelich.de<mailto:ra.biehl at fz-juelich.de>
Web:    http://www.fz-juelich.de/iff/d_ins

On 17.04.2012 11:54, Ramon Crehuet wrote:
Thanks Ralf and Konrad for the explanation,
Let me see if I got it right:

* EnergeticNormalModes diagonalizes the hessian matrix (K) , in cartesian coordinates. If temperature = None, eigenvectors are normalized, otherwise the norm is temperature-dependent. Therefore dotProduct directly tests for orthogonalisation.

* VibrationalNormalModes (or default NormalModes) diagonalizes the mass-weighted hessian matrix : K' = 1/sqrt(M) K 1/sqrt(M))  where M are the masses of the atoms. The resulting eigenvectors are expressed also in mass-weighted cartesian coordinates. Therefore one needs to use the massWeighted versions of dotProduct and norm to work with them. In particular, a vector e' in mass-weighted coords is related to the cartesian vector e as:
e'=1/sqrt(M)e

So in order to have the eigenvectors of K' in cartesian corrdinates one should do something like (given m1 and m2 as vibrational modes):
m = m1.universe.masses().array
m1cart = np.ravel((m1.array) * np.sqrt(m[:, np.newaxis]))
m2cart = np.ravel((m2.array) * np.sqrt(m[:, np.newaxis]))
It is then true that m1cart.dot(m2cart)=0. and m1cart.dot(m1cart)=1.
(dotProduct no longer applies because these are numpy arrays, not VibrationalMode objects)
I hope this is correct and I think it will help future readers of the list.
Cheers,
Ramon



2012/4/12 Konrad Hinsen <research at khinsen.fastmail.net<mailto:research at khinsen.fastmail.net>>
--On 12 avril 2012 13:13:04 +0200 Ramon Crehuet <rcrehuet at gmail.com<mailto:rcrehuet at gmail.com>> wrote:

Hi,
I have performed a normal mode calculation (Calpha) and I am not sure I
understand how to scale the modes so that they are orthonormalized.

Which kind of mode do you use?

EnergeticModes: the mode vectors are orthogonal to each other, i.e.

 modes[i].dotProduct(modes[j])

is zero for i != j. For i==j, you get the thermal amplitude of the mode.

VibrationalModes: the mode vectors are in unweighted Cartesian coordinates, meaning they are not orthogonal. However, you have

 modes[i].massWeightedDotProduct(modes[j])

which should be zero for i == j.

For both kinds of modes, you can specify temperature=None to skip scaling by the thermal fluctuation amplitudes. Then you should see the norms being equal to one.

Konrad.





_______________________________________________
mmtk maillist  -  mmtk at starship.python.net<mailto:mmtk at starship.python.net>
http://starship.python.net/mailman/listinfo/mmtk



------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Forschungszentrum Juelich GmbH
52425 Juelich
Sitz der Gesellschaft: Juelich
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
Vorsitzender des Aufsichtsrats: MinDir Dr. Karl Eugen Huthmacher
Geschaeftsfuehrung: Prof. Dr. Achim Bachem (Vorsitzender),
Karsten Beneke (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt,
Prof. Dr. Sebastian M. Schmidt
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------

Kennen Sie schon unsere app? http://www.fz-juelich.de/app
-------------- n?chster Teil --------------
Ein Dateianhang mit HTML-Daten wurde abgetrennt...
URL: <http://starship.python.net/pipermail/mmtk/attachments/20120417/32a63054/attachment.html>


More information about the mmtk mailing list