Software for Dynamics of Electrons and Nuclei in Molecules.

Author: Erik Deumens.

ENDyne is an application that implements the Electron Nuclear Dynamics (END) theory for studying the interaction between molecular geometry and electronic structure in a time-dependent and self-consistent way.

The theory is somewhat unfamiliar to most people and the software is not very user friendly, for that reason we do not make the code available to a general audience at this time. However, if someone is interested in the code, they can come and study with us for about a month and they get to take the code with them at the end. Please send e-mail. for more information.

Program structure

         ENDyne              Python
         _______      _______________________
         endmain      "import endyne"  PyRun_
            |            |               ^
            v            v               |             
        .   \ \          |               |	.
        .    \ v         v               |	.
        .     \ pythonendcmds      pythonembed	.
        .      \      |                ^	.
                 \    |              /
                  \   |             /
                   \  |            /
                    v v           /


Feaures under development for version 5

  • Multi configuration wave functions for electrons and nuclei
  • Frozen electron cores
  • Semi-empirical integrals
  • Density Functional implementation
  • Python programming interface
  • Source in Fortran 95, ANSI C and python.
  • MPI massively parallel
  • POSIX threads shared memory parallel
  • Use of QTIP integrals

New features in version 5


List of developers

Author Project
Trygve Helgaker McMurchie-Davisson gaussian orbital integrals
Augie Diz SCF and state projection for diatomic
Hugh Taylor vibration analysis for diatomic
Benny Mogensen complex gaussian orbital integrals
Jorge Morales differential cross sections with semiclassical corrections (Airy and uniform approximation)
Mauricio Coutinho general state projection
Anatol Blass general rotation-vibration analysis
Remigio Trujillo differential cross sections with semiclassical corrections (Schiff approcimation); make movies with sphere representation of the Mulliken population.
Denis Jacquemin Semi-empirical integrals (AM1) and Obara-Saika (PRISM) integrals.
Oliver Quinet Python scripts for managing differential cross-section computations; make movies with true electron densities.
Ben Hall Semi-classical nuclear wave function.
Python interface The program evolve no longer exists in version 5. It is replaced by Python. Python is an Opensource project featuring a rapid prototying laguage that is both a simple scripting language and a fully featured object oriented language.ENDyne is defined as an extension of Python, which means that you can start the python interpreter and load ENDyne to make the ENDyne commands and data available in the interactive session or Python script. A session looks like:

edsparc% python
Python 2.0 (#3, Sep  3 2001, 11:50:37) [C] on sunos5
Type "copyright", "credits" or "license" for more information.
>>> import endyne
>>> endyne.version()
(5, 'A', 3, 1, 6, 'generic', 'Tue Jul 30 12:52:28 EDT 2002')

Python is also embedded in the endyne executable, so that you can run endyne with a .inp file as argument to make a calculation as in version 2:

edsparc% endyne test.inp
ENDyne version 5.A.3.1 test.inp

                             *  E N D y n e  *

                   Executable version  5 alpha development phase  3
                           Restart file format  6

Or you can run endyne with a .rsl file as argument to get the functionality of evolve:

edsparc% endyne test.rsl
ENDyne version 5.A.3.1 test.rsl
RSL file test.rsl opened with the command:
    rslfile = RSLopen('test.rsl','r')
Recorded state log (RSL) file
 name:    test.rsl
 mode:    r
 format:   6
 version:      5.0603
 begin:    0.0E+0 1
 end:      0.1 29
 current:  0.0E+0 -1
 records:  11
 flags:   ownself       T
          seekdone      T
          readonly      T
          flushbuf      F
          isopen        T
          rstisopen     F


All the functionality of ENDyne can be accessed from Python and many of the new features are only available through the Python interface. Documentation of any command is available with the standard Python method:

endyne>>> print thing.__doc__
Schiff Approximation When a bundle of trajectories as a function of orientation angles of projectile and target and of impact parameters has been calculated, it is possible to get the differential cross section with the Python script dcs.py, which use the Python library dcslib.py. The semi-classical Schiff Approximation correction can be obtained with the Python function schiff.
Electronic subspace projection The Mulliken population analysis provides away to estimate the probability of charge transfer. A more detailed way involves constructing the projection of the evolved state on the space of all determinants with given charge and spin distribution. This is done with the project method of RSLdata objects.
Thouless ordering By setting LORDV2=T, the algorithm that determines the best reference for Thouless representation now takes into account the charge of each nucleus and builds a reference state that has a number of electrons on each nucleus in agreement with chemical intuition.To allow application to negative ions, the algorithm allows putting too many elecrons on a center, after all centers are chemically neutral.
Backstep recovery Version 2 propagation relies on the RSL (recorded state log) file to recover from a problem where further integration is not possible. This happens because the integrator has propagated into a region of phase space where various operations to determine the forces on nuclei and electrons become numerically ill-conditioned. Continuing would lead to inaccurate and incorrect results. If the RSL step RDT is too large, the backstep may not yield a numerically more stable reference and then integration is aborted. This is known as a broken trajectory.Version 5 uses an internal mechanism to store the last known good state in RAM. This method is more flexible and requires little or no user intervention. It works effectively with the Hindmarsh-Gear LUSEHG=T and Shampine-Gordon LUSESG=T methods. Because the Bulirsch-Stoer LUSEBS=T method does not step forward uniformly, the RAM backstep method is really equivalent to the old RSL method.

The use of the parameter RCNDLM is much more sensitive in the new algorithm. Therefore users who have set this parameter, should remove it from their input files to ensure the use of the default. A new parameter RCNDAC has been introduced, giving the relative accuracy of the condition number returned by LAPACK routines. The version 2 algorithm considered a new state bad if

new_condition_number > RCNDLM * NDIM * old_condition_number
independent of how close in time the previous condition number was obtained. The version 5 algorithm uses the inequality
d/dt ln condition_number > RCNDLM * NDIM
within the accuracy RCNDAC. Setting RCNDAC=1. and RCNDLM large ensures that the test is never satisfied. The resulting trajectory may then inaccurate and eventually wrong.

Project and target orientation To compute differential cross sections for reactions involving non-spherical systems, it is necessary to compute sets of collision trajectories as a function of impact parameter for a set of orientations of the projectile and the target. The Python script prep.py will set up a directory tree and populate it with .inp files for ENDyne and .job files for LoadLeveler by editing tenplates provided by the user. The routines to average cross sections accurately over all orientations are provided in orientlib.py. These routines are used by dcs.py and dcslib.py.When both target and projectile are atoms, there are three relative coordinates:
1) b: the impact parameter along the b axis,
2) z: the initial projectile distance along the beam axis,
3) phi: the angle defining the b-axis beam-axis plane around the beam axis.
For the most general case of two molecules, there are nine relative coordinates: b, z, phi as in the atom-atom case, plus:
4) alpha: azimuth for the projectile axis,
5) beta: azimuth for the target axis,
6) gamma: dihedral angle between the projectile axis and target axis,
7) chi: dihedral angle between the b-axis and beam-axis plane, and the projectile axis,
8) delta: projectile spin angle around the projectile axis,
9) epsilon: target spin angle around the target axis.There are 6 distinct types with increasing number of relative coordinates.
39 The case of an atom and an atom has just three: b, z, phi.
59 The case of an atom and a diatomic molecule has five: b, z, phi, beta, chi.
69 The case of an atom and a polyatomic has six: b, z, phi, beta, chi, epsilon.
79 The case of two diatomic molecules has seven: b, z, phi, alpha, beta, gamma, chi.
89 The case of a diatomic molecule and a polyatomic has eight: b, z, phi, alpha, beta, gamma, chi, epsilon.
99 The case of two polyatomics has all nine: b, z, phi, alpha, beta, gamma, chi, delta, epsilon.The program does not care whether the first cluster is the projectile or the target, hence to scatter a polyatomic on an atom, designate the atom as fragment 1 and target and the polyatomic as fragment 2 and projectile and use type 69.
Performance Monitor The program has been instrumented with calls to performance monitoring code. On all platform on which PAPI is supported, a package hpm.c provides a simple interface. This interface is the same as the one developed for IBM POWER3 and POWER4 processors by Luis DeRose called HPM Toolkit. On these platforms, ENDyne uses the HPM Toolkit. Monitoring is turned on with the flags LTRACE and LTRPRF.
Free electrons It is possible to construct a single determinant wave function with complex orbitals using basis functions with electron translation factors (ETF) placed on centers that can be atomic nuclei or freely moving dynamical centers. We refer to such wave functions as wave functions with free electrons. The dynamical equations with dictate the motion of the free centers and of the occupation of the basis functions on them.Wave function with free electrons can be used to describe such processes as electron scattering on atoms and molecules;ionization of atoms and molecules during colissions and reactive processes; and recombination of ions with free electrons.

The use of the much more stringent backstep algorithm implemented in version 5 is essential for wave functions with free electrons, as the numerical integrator can produce nonsense much more easily than for wave functions with only atomic centers.

The use of coordinates for the center of free moving centers of basis functions introduces a duplication with the basis functions exapansion coefficients. To obtain a singularity free coordinate system for the phase space it is essential to perform a transformation. A new way to solve the dynamical equations has been implemented which constructs the complete symplectic form. The TDVP equations were solved by substitution in version 2, to save the space for the full symplectic form. However, the transformation is only meaningful on the entire symplectic form. The solution method with the full symplectic form is used when LUSPFM=T, and the singularity removing transformation is applied when LUSPTR=T. The Hamiltonian equation is solved with a general linear eqaution solver by default; with LUSPVD=T the equation is solved by singular value decomposition.

By turning on tracing LTRACE=T and LTRSPD=T detailed information about the solution of the equations is printed to the log file and, with LUSPVD=T, the singular values and right singular vectors are written to the RSL file in TRVEC1 and TRMAT1, respectively. The components of the singular vectors are in the order of the internal representation of the symplectic equation: real part of the Thouless coefficients, center positions, imaginary part of the Thouless coefficients, center momenta. TRVEC2 and TRVEC3 contain the velocities as computed by the substitution method (version 2) and by solving the symplectic form (version 5, LUSPFM=T). TRVEC4 holds the solution with the symplectic transformation LUSPTR=T. The velocity components are in the following order: Thouless coefficients, center positions, center momenta. TRVEC5 contains the force of the TDVP equation. The force components are in the same order as the velocities.

There is an extra singularity that can develop: the free center can get zero coefficients for all its basis functions. The algorithm to avoid this singularity in the nuclear metric uses an expanded state phase space. This feature can be turned on by setting LUSEXS=T.

Rules for choosing and constructing the basis set for free centers are still being developed. So far the only guideline is that there should be only basis functions of even angular momentum: s,d,g,j,… and not p,f,h,….

Optimization The strength of ENDyne has always been in the time propagation. The capability to optimize electronic parameters at a fixed geometry, or to optimize electronic parameters and geometry simultaneously has always existed but did not perform very well. For that reason a conventional SCF code based on DIIS has been implemented to allow constructing initial SCF wave functions.The optimization capability has now been completely disentangled from the propagation capability, thus increasing the stability and maintainability of the software system. The optimization functionality appears to be much improved.
Phase normalization Routines that diagonalize matrices, do not produce thesame vectors, they are determined only up to a complex phase. Now all Fock eigenvectors are normalized so that the largest component is real and positive.Reordering can produce different looking Thouless coefficientcs. Now the Thouless coefficients printed in the output at the end of a run are normalized so that the largest component is equal to 1. The speed vectors are similarly normalized, which makes no sense!
TDA and RPA spectra and modes New flags LCIEIG and LCIVEC control that the configuration interaction singles or Tamm-Dancoff approximation energies and vectors will be computed with the current state as reference, which need not be an SCF state. Similarly, LIPEIG and LIPVEC control that the interaction picture or Random Phase approximation energies and vectors will be computed. These properties require the LMKPRP is set.
Interaction picture coordinates for propagation The flag LUSEIP turns on the use of a transformation of the electronic variables to interaction picture variables, i.e. RPA modes, and propagation in these variables. This has the potential to significantly speed up propagation of electronic variables. The variables are monitored and when the forces get too big, indicating they are no longer optimal, evolution is stopped and a new transformation is forced before evolution is resumed.
New computational kernel The computational kernel has been rewritten completely. The old HFFRC routine with its complicated spin block structure has been broken up into smaller pieces with a simpler spin orbital structure more close in design to the structure required to support multi-configuration wave functions. The logic to determine a new mapping of basis functions to occupied/unoccupied spaces needed to transform to Thouless coefficients has been modified to be closer to the fully general design needed for MC. A new input switch LMKTHC has been introduced to allow turning off generation of a new Thouless map, thus keeping the existing map, e.g. from the RSL file when LRSTRT=T.
Vector Hartree Fock The Complete Active Space (CAS) multiconfiguration (MC) or multi reference (MR) (spin-)unrestricted vector coherent state wave function is being implemented and tested. It is more a generalization of Hartree-Fock than the conventional CAS-MC wave function. For that reason we call it the Vector Hartree-Fock wave function.
Semi-classical nuclear wave function The nuclear de Broglie phase factor is integrated numerically alongthe trajectories and multiple trajectories are allowed to interfere. This is an implementation of the semi-classical nuclear wave function to generalize the Schiff approximation to the classical cross sections.
Semi empirical integrals The integral generation module has been completely restructured and now delivers semi empirical integrals and ab initio integrals in a consistent way. The code to compute contractions of integral patches with densities has been restructured as well. The new algorithm for evaluating the contractions is the same for MD model-forces, MNDO methods and ab initio methods.

Software versions

Version Release Update Date Feature
2 7 6 Jan 1, 1998 Production; general molecules; electron translation factors
Version Release Update Date Feature
2 8 1 Jan 1, 2000 Fortran 77; use QTIP integrals
2 8 1 Nov 1, 2000 Fortran 90 rewrite
Version Release Phase Date Feature
5 Alpha 2 Mar 10, 2002 python integration; state projection; schiff approximation; use QTIP 1.A.4
5 Alpha 3 Jul 31, 2002 port to MIPS, Alpha, x86; use QTIP 1.A.4
5 Alpha 4 Dec 31, 2003 new spin loops; TDA and RPA; use QTIP 1.A.4
5 Alpha 6 Q3, 2008 Vector Hartree-Fock; semi-classical nuclear wave function; use QTIP 1.A.4
5 Alpha 5 Q1, 2009 semi empirical integrals; use QTIP 1.A.5
Version Release Update Date Feature
5 1 0 TBA vector Hartree-Fock wave function; use QTIP 1.1

Hardware and software platforms tested

ENDyne version 2

Hardware OS Fortran 90
SUN SPARC Solaris 7 Forte 6.2
IBM PowerPC AIX 4.3.3 XLF 7.1
SGI MIPS IRIX 6.5 Compilers 7.3
Linux RH ES 4.5 Intel Fortran 9.1

ENDyne version 5

Hardware OS C Fortran 90 Python
SUN SPARC Solaris 7/8/9 Forte 6.2 Forte 6.2 2.1
AIX 5L 5.1 VAC 6.0 XLF 8.1 2.1
Linux RH ES 4.5 gcc 3.3 Intel Fortran 9.1 2.1

Theoretical reference papers

  • E. Deumens, Y. Öhrn, and B. Weiner,
    Coherent state formulation of multiconfiguration states
    J. Math. Phys. 32, 1166-1175 (1991).
  • E. Deumens, A. Diz, R. Longo and Y. Öhrn,
    Time-dependent theoretical treatments of the dynamics of electrons and nuclei in molecular systems,
    Rev. of Mod. Phys., 66, 917-983 (1994).
  • E. Deumens and Y. Öhrn,
    Wave function phase space: An approach to the dynamics of molecular systems
    J. Chem. Soc. Faraday Trans., 93, 919 (1997).
  • Y. Öhrn and E. Deumens,
    Towards an ab initio treatment of the time-dependent Schrödinger equation of molecular systems,
    J. Phys. Chem., A103, 9545 (1999).