[MMTK] symeig change to eigh

ra.biehl ra.biehl at fz-juelich.de
Mon Apr 23 10:43:41 UTC 2012


Dear Konrad

scipy.linalg.eigh has the option to calculate only the set of n to m eigenmodes.
I attach a small diff snippet for Core.py and VibrationalModes.py to allow this, if
scipy.linalg.eigh is there. The others XXXXModes.py files need the same change.

On a system with 6350 atoms the full normal mode calculation time (Amber FF: about 14000s) is reduced by a factor of 3 if only the first 100 normal modes (Amber FF: about 5000s) are calculated.
Just as an Idea if you want to integrate it.

I tried to understand the change for SparseMatrixNormalModes, but for that the mentioned interface to ARPACK is necessary. Due to that i found that there is an alternative interface, again in scipy, that uses ARPACK.
It is found in scipy.sparse.linalg.eigh. Perhaps this solves the somewhere mentioned licensing problem and will allow the usage of SparseMatrixNormalModes
.

Anyway the subspace methods seems to give best result in terms of small ram usage and computing time.



Viele Grüsse
and best wishes

Ralf Biehl





  *** temp/NormalModes/Core.py  2011-09-21 17:52:26.000000000 +0200
--- download/MMTK-2.7.5/MMTK/NormalModes/Core.py        2012-04-19 13:24:52.339737448 +0200
***************
*** 19,21 ****

! symeig = None
   if array_package == "NumPy":
--- 19,21 ----

! eigh = None
   if array_package == "NumPy":
***************
*** 27,29 ****
       try:
!         from symeig import symeig
       except ImportError:
--- 27,29 ----
       try:
!         from scipy.linalg import eigh
       except ImportError:
***************
*** 259,263 ****
           # Calculate eigenvalues and eigenvectors of self.array
!         if symeig is not None:
               _symmetrize(self.array)
!             ev, modes = symeig(self.array, overwrite=True)
               self.array = N.transpose(modes)
--- 259,268 ----
           # Calculate eigenvalues and eigenvectors of self.array
!         if eigh is not None:
               _symmetrize(self.array)
!             if self.eigvals is not None:
!                 lo,hi=self.eigvals
!             else:
!                 lo,hi=0,self.array.shape[0]-1
!             print lo,hi
!             ev, modes = eigh(self.array, eigvals =(lo,hi),overwrite_a=True,overwrite_b=True)
               self.array = N.transpose(modes)
***************


  *** temp/NormalModes/VibrationalModes.py      2011-05-11 17:36:15.000000000 +0200
--- download/MMTK-2.7.5/MMTK/NormalModes/VibrationalModes.py    2012-04-18 16:24:40.603793490 +0200
***************
*** 71,73 ****
       def __init__(self, universe=None, temperature = 300.*Units.K,
!                  subspace = None, delta = None, sparse = False):
           """
--- 71,73 ----
       def __init__(self, universe=None, temperature = 300.*Units.K,
!                  subspace = None, delta = None, sparse = False,eigvals=None):
           """
***************
*** 125,127 ****
           self.temperature = temperature
!
           self.weights = N.sqrt(self.universe.masses().array)
--- 125,127 ----
           self.temperature = temperature
!         self.eigvals=eigvals
           self.weights = N.sqrt(self.universe.masses().array)
***************
*** 136,138 ****
           self.sort_index = N.argsort(self.frequencies)
!         self.array.shape = (self.nmodes, self.natoms, 3)

--- 136,138 ----
           self.sort_index = N.argsort(self.frequencies)
!         self.array.shape = (-1, self.natoms, 3)



On 02.04.2012 08:44, Konrad Hinsen wrote:
> ra.biehl writes:
>
>   >  I was trying to reduce a bit the memory usage of the normal mode
>   >  calculation of a large system.
>
> That's not easy if you want to remain general. The only approaches I
> know that can reduce memory use significantly are
>
> - reduce the system size (coarse-grained models)
>
> - reduce the number of modes (e.g. by using MMTK.Subspace)
>
> - use iterative methods with a sparse representation of the force constant
>    matrix to calculate selected modes
>
>   >  This project or better the source is not available anymore, but the wrapper for the
>   >  underlying lapack functions has moved to scipy and is included as scipy.linalg.eig
>   >  Changing the small pieces attached at the end allow the usage of scipy.linalg.eigh
>   >  instead of symeig in NormalModes/Core.py
>
> I'll apply those modifications to MMTK, thanks!
>
>   >  Testing with AmberFF and insulin.pdb in mode='all' the Normal mode
>   >  calculation is a bit faster (232 s instead of 291s ) but
>   >  unfortunately the memory usage is not changed.
>
> That's what I would expect.
>
> If you want to try iterative approaches, I should still have some old
> interface to ARPACK that does the job, assuming that ARPACK didn't
> change too much over the last 10 years.
>
> Konrad.


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



More information about the mmtk mailing list