Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision
Previous revision
ab_initio_multiple_spawning_dynamics [2020/06/11 18:17] – external edit 127.0.0.1ab_initio_multiple_spawning_dynamics [2024/02/14 10:54] (current) – external edit 127.0.0.1
Line 1: Line 1:
 +====== Ab initio multiple spawning dynamics (AIMS) ======
 +
 +//AIMS// is a multi-state first-principles dynamics program written by\\
 +Benjamin G. Levine, Joshua D. Coe, Aaron M. Virshup, Hongli Tao, Christian R. Evenhuis, William J. Glover, Toshifumi Mori, Michal Ben-Nun and Todd J. Martinez.\\
 +
 +
 +Bibliography:
 +
 +B. G. Levine, J. D. Coe, A. M. Virshup and T. J. Martinez, Chem. Phys. 347, 3 (2008).\\
 +M. Ben-Nun and T. J. Martinez, [[https://dx.doi.org/10.1016/S0009-2614(98)01115-4|Chem. Phys. Lett.]] **298**, 57 (1998).
 +
 +All publications resulting from use of this program must acknowledge the above. See also:
 +
 +B. G. Levine and T. J. Martinez, J. Phys. Chem. A 113, 12815 (2009).\\
 +T. J. Martinez, Acc. Chem. Res. 39, 119 (2006).
 +
 +The AIMS module implements the Ab Initio Multiple Spawning method to perform dynamics calculations on multiple electronic states. Although the program was designed for non-adiabatic dynamics with CASSCF wavefunctions, it can be used quite generally for first principles molecular dynamics, provided that nuclear gradients are available.
 +
 +AIMS provides a description of a time-evolving molecular wavepacket with a multi-configurational, product gaussian function basis. To allow the size of the basis set to remain small while describing the wavepacket dynamics, the centers of the full-dimensional gaussian functions follow classical trajectories while their widths are kept fixed (the frozen Gaussian approximation). When these classical trajectories encounter regions of large non-adiabatic coupling to another electronic state, new basis functions are //spawned// on the coupled state and the Time Dependent Schrödinger Equation is solved in order to describe population transfer between the electronic states. In this way, the description of the wavepacket can be systematically improved by increasing the number of initial trajectory basis functions, while reducing the coupling threshold at which trajectories spawned. Alternatively, spawning can be disabled, and the method will reduce to traditional first-principles classical molecular dynamics when a single trajectory basis function is used.
 +
 +===== Compilation =====
 +
 +AIMS is compiled as a component of Molpro; however, it requires a fortran 2003 compiler and is disabled unless the preprocessor variable MOLPRO_f2003 is defined, which depends on the version of your compiler.
 +
 +===== Structure of the input =====
 +
 +The AIMS program is controlled by the Molpro input deck and a number of external text files:
 +
 +==== Molpro input deck ====
 +
 +The AIMS program is invoked by the command AIMS. The details of the electronic structure are defined in specifically named procedures. At a minimum, for first principles dynamics, AIMS requires the following structure to the input deck:
 +
 +''%%SYMMETRY, nosym%%''\\
 +''%%ORIENT, noorient%%''\\
 +''GEOMTYP=xyz''\\
 +''%%GEOMETRY={ Initial atom specifications. Must match ordering of Geometry.dat, but need not be the same geometry }%%''
 +
 +''BASIS=''{//Basis specifications.//''%%}%%''
 +
 +''%%MPINIT={ Electronic structure routines for initial wavefunction }%%''
 +
 +''%%MPCALC={ Electronic structure routines for dynamics calculations. At a minimum, must include: ener(1)=energy(1) ener(2)=energy(2)etc force }%%''
 +
 +''%%AIMS,%%''//dir//
 +
 +where
 +
 +  * **''dir''** (string) Relative directory to AIMS input files. (Default is working directory).
 +
 +Note: the above is appropriate only for first-principles dynamics (without spawning). CASSCF spawning calculations require additional procedures and a specific structure to the MPINIT and MPCALC procedures. See the examples in section [[ab initio multiple spawning dynamics#examples|examples]].
 +
 +==== Control.dat ====
 +
 +Control.dat is a fortran NAMELIST file storing the main parameters (in au) of an AIMS simulation. At the end of the file, one can optionally specify user-defined atom types (see below). Also time steps are in au, 1 au = 0.0241888 femtoseconds.
 +
 +^  Initialization  ^^^
 +^Parameter             ^Options (Default)          ^Description                                      ^
 +|Required parameters:  |||
 +|''NumStates''     |''Integer''            |Number of electronic states.                     |
 +|''InitState''     |''Integer''            |Initial electronic state.                        |
 +|''NumParticles''  |''Integer''            |Number of atoms.                                 |
 +|Common parameters:    |||
 +|''IMethod''       |''0''                  |First principles dynamics (No NACV).             |
 +|                      |''%%(1)%%''                |CASSCF dynamics, with NACV.                      |
 +|''IRestart''      |''%%(0)/1%%''              |Restart the calculation when set to 1.           |
 +|''RestartTime''   |''%%Real (-100.0)%%''      |Time from which to restart in Checkpoint.txt. If disabled, use Last_Bundle.txt.                |
 +|''IRndSeed''      |''%%Integer (0)%%''        |Initial condition random number seed.            |
 +|''InitialCond''   |''%%('NoSample')%%''       |Read initial conditions from Geometry.dat.       |
 +|                      |''%%'Wigner'%%''           |Sample from the Wigner distribution of harmonic normal modes in Frequencies.dat.                 |
 +|                      |''%%'Quasiclassical'%%''   |As above, for the Quasiclassical distribution.   |
 +|                      |''%%'Boltzmann'%%''        |As above, for the Boltzmann distribution.        |
 +|''Temperature''   |''%%Real (0.0)%%''         |Temperature (Kelvin) of the initial conditions.  |
 +|Expert parameters:    |||
 +|''InitBright''    |''%%Logical (.false.)%%''  |Select optically brightest of ''InitState'' and ''IDark'' as initial state.              |
 +|''InitDark''      |''%%Logical (.false.)%%''  |As above, but for darkest state.                 |
 +|''IDark''         |''%%Integer (-1)%%''       |See above.                                       |
 +|''InitGap''       |''%%Real (0.0)%%''         |Select initial ground-excited energy gap to within $\pm 0.5\times$''InitGapWidth''  |
 +|''InitGapWidth''  |''%%Real (0.0)%%''         |See above.                                       |
 +|''IRestartTraj''  |''%%Integer array (0)%%''  |Subset of trajectories to restart from. Can restart a single dead trajectory.                |
 +|''NumInitBasis''  |''%%Integer (1)%%''        |Number of initial coupled basis functions.       |
 +|''FirstGauss''    |''%%Logical (.false.)%%''  |Place the first basis function at the molecular origin (that of Geometry.dat), with zero momentum.                              |
 +|''MirrorBasis''   |''%%Logical (.false.)%%''  |Mirror the initial basis on another electronic state, ''MirrorState''          |
 +|''MirrorState''   |''%%Integer (1)%%''        |See above.                                       |
 +|''EnergyAdjust''  |''%%Logical (.false.)%%''  |Adjust momenta of mirrored basis (along NACV) to match original basis energy.            |
 +|''SharpEnergy''   |''%%Logical (.false.)%%''  |Adjust momenta of all initial basis functions to match first basis function energy.            |
 +
 +
 +^Integration and numerics  ^^^
 +^Parameter                   ^Options (Default)              ^Description                                 ^
 +|Required parameters:        |||
 +|''TimeStep''            |''Real''                   |Simulation timestep.                        |
 +|''SimulationTime''      |''Real''                   |Simulation length.                          |
 +|Common parameters:          |||
 +|''ExShift''             |''%%Real (0.0)%%''             |Apply a shift of ''%%+ExShift%%'' to diagonals of Hamiltonian to improve amplitude/phase propagation. Good choice is average potential energy of ground and excited states.               |
 +|''CoupTimeStep''        |''%%Real (TimeStep/4)%%''      |Timestep to use in regions of large nonadiabatic coupling.                      |
 +|Expert parameters:          |||
 +|''TimeStepRejection''   |''%%Logical (.true.)%%''       |Enables adaptive time-stepping when encountering integration problems.          |
 +|''MinTimeStep''         |''%%Real (CoupTimeStep/4)%%''  |Minimum timestep to use when adaptive time-stepping.                     |
 +|''DieOnMinStep''        |''%%Logical (.true.)%%''       |Stop calculation if minimum timestep is reached.                                 |
 +|''EnergyStepCons''      |''%%Real (0.005)%%''           |Allowed violation of classical energy conservation before timestep is rejected.                                   |
 +|''NormCons''            |''%%Real (0.25)%%''            |Permitted deviation in Norm from beginning of simulation.                    |
 +|''RejectAllStateFlip''  |''%%Logical (.true.)%%''       |Reject a timestep if any two electronic states flip. If ''.false.'' then only reject if state is coupled to occupied state.                          |
 +|''OlapThresh''          |''%%Real (0.001)%%''           |Threshold overlap between two trajectories below which matrix elements not calculated.                    |
 +|''RegThresh''           |''%%Real (0.0001)%%''          |Threshold to regularize the overlap matrix before inversion.                    |
 +|''Constrain''           |''%%Logical (.false.)%%''      |Apply geometric constraints? If yes, supply Constraints.dat                      |
 +
 +^Spawning  ^^^
 +^Parameter                ^Options (Default)          ^Description                                           ^
 +|Common parameters:       |||
 +|''CSThresh''         |''%%Real (0.1)%%''         |Threshold Coupling, above which to spawn.             |
 +|''PopToSpawn''       |''%%Real (0.1)%%''         |Minimum population a trajectory must have before it can spawn.                                  |
 +|''OMax''             |''%%Real (0.8)%%''         |Prevent spawning if overlap of spawned trajectory with existing basis greater than ''Omax''                                   |
 +|''MaxTraj''          |''%%Real (100)%%''         |Maximum number of allowed trajectories. Above which, spawning is disabled.                    |
 +|''IgnoreState''      |''%%Integer (0)%%''        |Stop dynamics of trajectories on ''IgnoreState''. if not coupled for at least ''DecoherenceTime''. Choose 1 if only interested in excited states.        |
 +|''MaxEDiff''         |''%%Real (0.03)%%''        |Energy gap threshold for calculating nonadiabatic couplings, above which, coupling taken to be zero.|
 +|Expert parameters:       |||
 +|''DecoherenceTime''  |''%%Real (200.0)%%''       |See above.                                            |
 +|''SpawnCoupV''       |''%%Logical (.false.)%%''  |Use Coup.velocity as criterion for spawning rather than Coupling magnitude alone.                 |
 +|''EShellOnly''       |''%%Logical (.true.)%%''   |Require spawned trajectories to have same classical energy as parent by adjusting momentum along NACV.                                  |
 +
 +^  Output  ^^^
 +^Parameter                        ^Options (Default)          ^Description                                           ^
 +|Common parameters:               |||
 +|''zAllText''                 |''%%Logical (.false.)%%''  |Write all output files?                               |
 +|''zEDatFile''                |''%%Logical (.true.)%%''   |Write Energy file?                                    |
 +|''zNDatFile''                |''%%Logical (.true.)%%''   |Write Population file?                                |
 +|''zTrajFile''                |''%%Logical (.false.)%%''  |Write Trajectory files?                               |
 +|''zAmpFile''                 |''%%Logical (.false.)%%''  |Write Amplitude files?                                |
 +|''zPotEnFile''               |''%%Logical (.false.)%%''  |Write Potential Energy files?                         |
 +|''zXYZ''                     |''%%Logical (.false.)%%''  |Write XYZ Geometry files?                             |
 +|''zCoupFile''                |''%%Logical (.false.)%%''  |Write Coupling files?                                 |
 +|''zChargeFile''              |''%%Logical (.false.)%%''  |Write Charge files?                                   |
 +|''zDipoleFile''              |''%%Logical (.false.)%%''  |Write Dipole files?                                   |
 +|''zTDipoleFile''             |''%%Logical (.false.)%%''  |Write Transition Dipole files?                        |
 +|''zQMRRFile''                |''%%Logical (.false.)%%''  |Write $\langle r^2 \rangle$ quadrupole moment files?  |
 +|''zCIVecFile''               |''%%Logical (.false.)%%''  |Write CI Vector files? (Not recommended for large CI)                             |
 +|''zCorrfile''                |''%%Logical (.false.)%%''  |Write Correlation function file?                      |
 +|''zBundMatFile''             |''%%Logical (.false.)%%''  |Write Matrices (H, S, etc)?                           |
 +|''WriteMolden''              |''%%Logical (.false.)%%''  |Write Molden files?                                   |
 +|''MoldenStep''               |''%%Integer (200)%%''      |How often (every number of timesteps) Molden files are written.                  |
 +|''RestartStep''              |''%%Integer (4)%%''        |How often (every number of timesteps) Restart files are written.                 |
 +|Expert parameters:               |||
 +|''NStepToPrint''             |''%%Integer (1)%%''        |How often (every number of timesteps) output files are written.                  |
 +|''WriteEveryStep''           |''%%Logical (.true.)%%''   |Write output every timestep, including subtimesteps of adaptive propagation. Overrides ''NStepToPrint''. |
 +|''NBonds''                   |''%%Integer (0)%%''        |Number of bond lengths to print. Define IBond below.                            |
 +|''%%IBond(2,NBonds)%%''          |''%%Integer (0)%%''        |List of atom indices that are bonded. See NBonds above.                         |
 +|''NAngles''                  |''%%Integer (0)%%''        |Number of bond angles to print. Define IAngle below.                           |
 +|''%%IAngle(3,NAngles)%%''        |''%%Integer (0)%%''        |List of atom indices that define bond angles.                                   |
 +|''NDihedrals''               |''%%Integer (0)%%''        |Number of dihedral angles to print. Define IDihedral below.                        |
 +|''%%IDihedral(4,NDihedrals)%%''  |''%%Integer (0)%%''        |List of atom indices that define dihedral angles.                               |
 +|''NPyrams''                  |''%%Integer (0)%%''        |Number of pyramidalization angles to print. Define IPyram below.                        |
 +|''%%IPyram(4,NPyrams)%%''        |''%%Integer (0)%%''        |List of atom indices that define pyramidalization angles.                       |
 +
 +**User-defined atom types**
 +
 +AIMS recognizes the atoms H, C, N, O, S, F, Cl and assigns appropriate masses and gaussian width parameters for these atoms. If you wish to modify these defaults, or if your system contains other atoms, you can define custom atom types at the end of Control.dat. The format is (in atomic units):
 +
 +''n''\\
 +''%%AtomName(1) AtomicNumber(1) AtomicMass(1) GaussianWidth(1)%%''\\
 +$\vdots$\\
 +''%%AtomName(n) AtomicNumber(n) AtomicMass(n) GaussianWidth(n)%%''
 +
 +==== Geometry.dat ====
 +
 +Geometry.dat contains the molecular geometry that will be used in the generation of initial conditions. Optionally, the momenta can also be specified, which is useful if\\
 +''%%InitialCond='NoSample'%%'':
 +
 +''%%UNITS=ANG or UNITS=BOHR %%''\\
 +''natom''\\
 +''%%AtomName(1)     x(1)      y(1)      z(1)%%''\\
 +$\vdots$\\
 +''%%AtomName(natom) x(natom)  y(natom)  z(natom)%%''\\
 +''%%Momenta:%%''\\
 +''%%AtomName(1)     px(1)     py(1)     pz(1)%%''\\
 +$\vdots$\\
 +''%%AtomName(natom) px(natom) py(natom) pz(natom)%%''
 +
 +==== FrequenciesMP.dat ====
 +
 +FrequenciesMP.dat contains the normal-mode frequencies and displacements of the molecule in its initial geometry and electronic state (usually the ground state). This information is used in the generation of initial conditions. If ''InitialCond''$\neq$''%%'NoSample'%%'' then either this file, Hessian.dat or Frequencies.dat (Sections [[ab initio multiple spawning dynamics#Hessian.dat|Hessian.dat]] and [[ab initio multiple spawning dynamics#Frequencies.dat|Frequencies.dat]] respectively) must be present. The format of FrequenciesMP.dat follows the Molpro output of a Frequencies calculation (note: in the first line, the number of normal modes has to be given; the rest is as in the Molpro output):
 +
 +<code - FrequenciesMP.dat>
 +6
 +                              1 A         2 A         3 A         4 A         5 A
 + Wavenumbers [cm-1]         1116.71     1725.49     1728.47     3435.76     3566.01
 + Intensities [km/mol]        157.34       17.93       16.24        1.67        0.53
 + Intensities [relative]      100.00       11.40       10.32        1.06        0.33
 +           GX1             -0.10597    -0.00088    -0.00021    -0.03801    -0.00009
 +           GY1              0.00067    -0.05838    -0.02628    -0.00017    -0.07554
 +           GZ1              0.00004     0.02639    -0.05807     0.00078     0.00226
 +           GX2              0.49134    -0.14946     0.19728     0.17300    -0.07729
 +           GY2              0.03972     0.54519     0.46553    -0.12339     0.02439
 +           GZ2             -0.19161     0.21763    -0.00782     0.51879    -0.19043
 +           GX3              0.49290    -0.09010    -0.22955     0.17771    -0.20720
 +           GY3              0.13665     0.43154    -0.28739    -0.39612     0.35792
 +           GZ3              0.13722    -0.52910     0.11024    -0.37702     0.36169
 +           GX4              0.48838     0.25184     0.03520     0.17755     0.28570
 +           GY4             -0.18573    -0.16545     0.18702     0.52182     0.66736
 +           GZ4              0.05376    -0.05521     0.70450    -0.15258    -0.20268
 +
 +                              6 A
 + Wavenumbers [cm-1]         3570.66
 + Intensities [km/mol]          0.61
 + Intensities [relative]        0.39
 +           GX1             -0.00060
 +           GY1             -0.00245
 +           GZ1             -0.07574
 +           GX2              0.28798
 +           GY2             -0.16745
 +           GZ2              0.68133
 +           GX3             -0.20790
 +           GY3              0.37733
 +           GZ3              0.33817
 +           GX4             -0.07175
 +           GY4             -0.17584
 +           GZ4              0.03308
 +</code>
 +
 +Note: Only the real, non-zero frequencies should be included in this file (i.e. those corresponding to vibrations). The first line of the file must be the number of frequencies. See section [[harmonic vibrational frequencies (FREQUENCIES)]] for details on how to perform a normal-mode frequencies calculation in Molpro.
 +
 +==== Hessian.dat ====
 +
 +As an alternative to FrequenciesMP.dat, AIMS can generate the normal mode frequencies after reading in the Hessian from Hessian.dat. The format of this file is:
 +
 +''natom''\\
 +''%%Hess(1,1) Hess(2,1) \ldots Hess(3*natom,3*natom)%%''
 +
 +Note: the Hessian should not be mass weighted. It is permissible to have carriage returns in the list of Hessian components.
 +
 +==== Frequencies.dat ====
 +
 +As an alternative to FrequenciesMP.dat and Hessian.dat, AIMS can read in normal-mode frequencies and vectors in AIMS format in Frequencies.dat. The format of this file is:
 +
 +''%%natom   nmodes%%''\\
 +''%%freq(1) freq(2) \ldots freq(nmodes)%%''\\
 +''%%modevec(1,1) modevec(2,1) \ldots modevec(3*natom,nmodes)%%''
 +
 +Note: Frequencies must be in cm$^{-1}$, but normal mode vectors can be mass scaled or not- AIMS will detect. Also, it is permissible to have carriage returns in the list of frequencies and mode vectors. This file is generated by AIMS if either FrequenciesMP.dat or Hessian.dat is used initially.
 +
 +==== Constraints.dat ====
 +
 +If ''Constraint=.true.'' then Constraints.dat must be supplied. This is a fortran NAMELIST file that defines the geometric constraints in the system. These constraints are enforced with the RATTLE algorithm. The parameters contained in this file are:
 +
 +^Constraint parameters  ^^^
 +^Parameter                          ^Options (Default)        ^Description                           ^
 +|''NumBonds''                   |''%%Integer (0)%%''      |Number of bond lengths to constrain. Define IBondPtcles below.             |
 +|''%%IBondPtcles(2,NumBonds)%%''    |''%%Integer (0)%%''      |List of atom indices that are bonded. See NumBonds above.           |
 +|''NumAngles''                  |''%%Integer (0)%%''      |Number of bond angles to constrain. Define IAnglePtcles below.            |
 +|''%%IAnglePtcles(3,NumAngles)%%''  |''%%Integer (0)%%''      |List of atom indices that define bond angles.                   |
 +|''NComAtoms''                  |''%%Integer (0)%%''      |Number of atoms that define the center-of-mass constraint. Define ICOMPtcles below.                     |
 +|''%%ICOMPtcles(NComAtoms)%%''      |''%%Integer (0)%%''      |List of atom indices that define the center of mass.                       |
 +|''Toler''                      |''%%Real (1E-4)%%''      |Rattle tolerance                      |
 +|''MaxIts''                     |''%%Integer (10000)%%''  |Maximum number of Rattle iterations   |
 +
 +Note, a maximum of 1000 constraints is allowed.
 +
 +===== Running a CASSCF spawning calculation =====
 +
 +This section gives an overview of how to set up and run a CASSCF ab initio multiple spawning calculation.
 +
 +1) The first step is to create a molpro input deck that controls how the CASSCF electronic structure calculations are performed. AIMS requires a fairly complicated structure to the molpro deck and it is highly recommended that you use the ethylene CASSCF example as a starting point (section [[ab initio multiple spawning dynamics#examples|examples]]). Control.dat should then be modified for your system and desired simulation parameters (section [[ab initio multiple spawning dynamics#Control.dat|Control.dat]]).
 +
 +2) Next, the initial conditions should be defined. This involves providing an initial geometry in Geometry.dat (section [[ab initio multiple spawning dynamics#Geometry.dat|Geometry.dat]]) and normal-mode frequencies and displacements in FrequenciesMP.dat (section [[ab initio multiple spawning dynamics#FrequenciesMP.dat|FrequenciesMP.dat]]), which will be used to generate a phase-space distribution that the initial conditions sample. Alternatively, one can also provide momenta in Geometry.dat and set ''%%InitialCond='NoSample'%%'' in Control.dat, and directly specify the initial conditions.
 +
 +3) The AIMS calculation is then initiated by running molpro on the input deck, in the usual fashion. In order to sample different initial conditions, multiple AIMS simulations should be run with different values of the random number seed ''IRndSeed'' (or different Geometry.dat files) - each of these should be run in a unique directory.
 +
 +===== Output =====
 +
 +Beyond the regular molpro output, which contains the electronic structure output, AIMS writes to a number of text files, depending on the output parameters in Control.dat (section [[ab initio multiple spawning dynamics#Control.dat|Control.dat]]). All values are given in atomic units unless otherwise specified. Electronic properties of a trajectory are evaluated at the center geometry of the trajectory basis function.
 +
 +  * **FMS.out** Current progress of simulation and any warning/error messages.
 +  * **E.dat** Wavepacket energies. QM energies are expectation values of the quantum operators, including coherences between the trajectory basis functions. CL energies are averages of the trajectories’ classical energies, weighted by the amplitude norm of each trajectory.
 +  * **N.dat** Electronic populations, the sum of which should be conserved.
 +  * **TrajDump.i** Dynamic variables for trajectory i.
 +  * **Amp.i** Amplitudes of trajectory i. Note: since the product Gaussian basis is non-orthogonal, the amplitude norms do not necessarily sum to 1.
 +  * **PotEn.i** Electronic energies for trajectory i. The last column displays the total classical energy of the trajectory, which should be conserved.
 +  * **positions.i.xyz** Snapshots of trajectory i’s center geometry in XYZ format (Angstrom).
 +  * **Coup.i** Nonadiabatic couplings of each state to the occupied state for trajectory i. The first nstate columns contain the magnitude of coupling vectors and the last nstate columns contain the projection of the coupling onto the velocity vector.
 +  * **Charge.i** Partial atomic charges corresponding to the occupied electronic state of trajectory i. The molpro input deck must contain a Mulliken population analysis directive in the ''MPCALC'' procedure.
 +  * **Dipole.i** Dipole moments of each electronic state for trajectory i.
 +  * **TDip.i** Transition dipoles from ground to excited states for trajectory i.
 +  * **QMRR.i** Isotropic quadrupole moment $\langle r^2 \rangle$ of each electronic state for trajectory i. The molpro input deck must contain ''%%EXPEC,QMRR%%''.
 +  * **CIVec.i** CI Vectors for each state of trajectory i.
 +  * **CFxn.dat** The overlap correlation function $\langle \Psi(t) | \Psi(0) \rangle$.
 +  * **H.dat** Hamiltonian matrix.
 +  * **S.dat** Overlap matrix. Note: the overlap between trajectories on different states is zero, but is displayed here as if they were on the same state.
 +  * **SDot.dat** The time derivative of the Overlap matrix.
 +  * **i.j.molf** Molden file for trajectory i at time j.
 +  * **Checkpoint.txt** Restart file with history.
 +  * **Last_Bundle.txt** Restart file for the last written timestep.
 +  * **Spawn.log** Details of successful spawns.
 +  * **FailSpawn.log** Details of unsuccessful spawns.
 +  * **Spawn.i** XYZ file of geometry where trajectory i was spawned.
 +
 +===== Examples =====
 +
 +The following two examples show how AIMS can be run in either single-state first-principle molecular dynamics mode, or multi-state spawning mode. The parameters found in these input files were chosen for illustrative purposes and are not necessarily appropriate for your system of interest. It is recommended the literature be consulted before performing AIMS simulations.
 +
 +Ethylene Hartree-Fock First-Principles Molecular Dynamics
 +
 +See examples/aims_ethyl_HF
 +
 +<code - examples/aims_ethyl_HF.inp>
 +! FMS-Molpro does not support parallel execution.
 +if(NPROC_MPP.gt.1) then
 +skipped
 +endif
 +if(.not.modul_aims) then
 +skipped
 +endif
 +
 +SYMMETRY,nosym
 +ORIENT,noorient
 +GEOMTYP=xyz
 +GEOMETRY={
 +6
 +ethylene ground-state geom
 +    0.0000000031    0.0000000000    0.6768496105
 +   -0.0000000013    0.0000000000   -0.6768496105
 +   -0.0000000017    0.9131426983   -1.2402810929
 +   -0.0000000095    0.9131426983    1.2402810929
 +   -0.0000000159   -0.9131426983    1.2402810929
 +    0.0000000047   -0.9131426983   -1.2402810929
 +}
 +BASIS=6-31g**
 +
 +MPINIT={
 +hf
 +}
 +
 +MPCALC={
 + ! Calculate Energy and Forces
 + show,itimestep,itrj;
 + hf
 + ener( 1 )=energy( 1 )
 + if(status.lt.0) then
 +    text,ERROR: could not converge SCF.
 +    STOP
 + endif
 + pop;density
 + force
 +}
 +
 +aims,dir='./aims_ethyl_HF'
 +</code>
 +
 +<code - examples/aims_ethyl_HF/Control.dat>
 + &control
 +
 +!------ Initialization -------!
 + InitState=1
 + iRndSeed=1
 +
 +!------ Simulation Properties -------!
 + IMethod=0              !Single- state dynamics - no couplings.
 + NumParticles=6         !Number of atoms
 + NumStates=1            !number of electronic states
 +
 +!-------- Integration and numerics --------!
 + TimeStep=20.0          !timestep for integration
 + SimulationTime=2000.0  !integration time
 + MinTimeStep=1.0        !minimum adaptive timestep
 +!------- Ouput options ---------!
 +zAllText=.false.
 +zEDatFile=.true.       !Energy file
 +zNDatFile=.false.       !Population file
 +zTrajFile=.true.      !Trajectory data file
 +zAmpFile=.false.        !Amplitude files
 +zPotEnFile=.true.      !Potential energy files
 +zXYZ=.true.            !XYZ files
 +zCoupFile=.false.       !Coupling files
 +zChargeFile=.false.     !Charge files
 +zDipoleFile=.false.     !Dipole files
 +zTDipoleFile=.false.    !Transition dipole files
 +
 +
 +/ !using default atom types, so user-defined types are not necessary
 + 
 +</code>
 +
 +<code - examples/aims_ethyl_HF/Geometry.dat>
 +   UNITS=BOHR
 +        6
 +C  0.000000  0.000000  1.279060
 +C  0.000000  0.000000 -1.279060
 +H  0.000000  1.725590 -2.343792
 +H  0.000000  1.725590  2.343792
 +H  0.000000 -1.725590  2.343792
 +H  0.000000 -1.725590 -2.343792
 +#Momenta
 +   0.000     0.000     0.000
 +   0.000     0.000     0.000
 +   0.000     0.000     0.000
 +   0.000     0.000     0.000
 +   0.000     0.000     0.000
 +   0.000     0.000     0.000
 +</code>
 +
 +Ethylene CASSCF Spawning Dynamics
 +
 +See examples/aims_ethyl_CASSCF
 +
 +<code - examples/aims_ethyl_CASSCF.inp>
 +!Ethylene (2/2)-CASSCF AIMS deck.
 +!This can be used as a template for other simulations by changing
 +!User options below. The deck calculates couplings for up to three
 +!states simultaneously coupled to the occupied state (NUMCOUP=3.0),
 +!which should suffice for most applications. If this needs to be increased,
 +!modify MPCALC accordingly.
 +
 +!User options
 +symmetry,nosym  !Do not change!
 +occ,9,context=MCSCF
 +closed,7,context=MCSCF
 +wf,16,1,0
 +nstate=3
 +orient,noorient !Do not change!
 +geomtyp=xyz     !Do not change!
 +geometry={
 +6
 +ethylene ground-state geom
 +    0.0000000031    0.0000000000    0.6768496105
 +   -0.0000000013    0.0000000000   -0.6768496105
 +   -0.0000000017    0.9131426983   -1.2402810929
 +   -0.0000000095    0.9131426983    1.2402810929
 +   -0.0000000159   -0.9131426983    1.2402810929
 +    0.0000000047   -0.9131426983   -1.2402810929
 +}
 +basis=6-31g*
 +
 +
 +!Aims variables - do not change!
 +! FMS-Molpro does not support parallel execution.
 +if(NPROC_MPP.gt.1) then
 +skipped
 +endif
 +if(.not.modul_aims) then
 +skipped
 +endif
 +
 +! tighten thresholds:
 +gthresh,twoint=1.0d-13
 +gthresh,energy=1.0d-7,gradient=1.0d-2
 +! records used by AIMS
 +orbrec=2101.2
 +cirec=2501.2
 +diabrec=3000.2
 +gradrec0=5101.1
 +do icoup=1,nstate-1
 +  gradrec(icoup)=gradrec0+icoup
 +enddo
 +
 +MPINIT={ !Find the initial CASSCF wavefunction
 +if(orbital.eq.0) then
 +  hf
 +  multi
 +  pspace,10.0
 +  wf;state,nstate
 +  maxiter,20
 +  orbital,orbrec
 +  ciguess,cirec
 +  save,det=cirec
 +  if (status.lt.0) then
 +    text, MCSCF failed in MPINIT.
 +    text, Attempting uncoupled iterations.
 +    {multi,so
 +    pspace,10.0
 +    wf;state,nstate
 +    orbital,orbrec
 +    ciguess,cirec
 +    save,det=cirec
 +    maxiter,50}
 +  endif
 +  if(status.lt.0) then
 +    text,ERROR: could not converge MCSCF.
 +    text,FMS should now stop
 +    STOP
 +  endif
 +endif
 +}
 +
 +MPCALC={ !Calculate Energy and Forces
 +if (jtrj.eq.0) then
 +  show,itimestep,itrj;
 +else
 +  show,itimestep,itrj,jtraj
 +endif
 +data,truncate,diabrec
 +data,truncate,gradrec0
 +do icoup=1,nstate-1
 +  data,truncate,gradrec(icoup)
 +enddo
 +data,copy, orbrec, diabrec
 +{multi
 +pspace,10.0
 +wf;state,nstate
 +maxiter,20
 +orbital,orbrec
 +ciguess,cirec
 +dm,orbrec
 +save,det=cirec
 +diab,diabrec,save=orbrec
 +}
 +if (status.lt.0) then
 +  text, MCSCF failed in MPCALC.
 +  text, Attempting uncoupled iterations.
 +  {multi
 +  pspace,10.0
 +  state,nstate
 +  orbital,orbrec
 +  ciguess,cirec
 +  dm,orbrec
 +  save,det=cirec
 +  {iterations
 +    do,uncouple,1,to,10}
 +  maxiter,30
 +  diab,diabrec,save=orbrec
 +  }
 +endif
 +if(status.lt.0) then
 +  text,ERROR: could not converge MCSCF.
 +  text,FMS should now stop
 +  STOP
 +endif
 +do istate=1,nstate
 +  ener( istate )=energy( istate )
 +enddo
 +pop;density,orbrec,state=stateind
 +if (numcoup.eq.0.0) then
 +  {multi
 +  pspace,10.0
 +  wf;state,nstate
 +  maxiter,20
 +  orbital,orbrec
 +  ciguess,cirec
 +  save,det=cirec
 +  diab,diabrec,save=orbrec
 +  cpmcscf,grad,stateind,save=gradrec0;
 +  }
 +elseif (numcoup.eq.1.0) then
 +  {multi
 +  pspace,10.0
 +  wf;state,nstate
 +  maxiter,20
 +  orbital,orbrec
 +  ciguess,cirec
 +  save,det=cirec
 +  diab,diabrec,save=orbrec
 +  cpmcscf,grad,stateind,save=gradrec0;
 +  cpmcscf,nacm,stateind,stateind2(1),save=gradrec(1);
 +  }
 +elseif (numcoup.eq.2.0) then
 +  {multi
 +  pspace,10.0
 +  wf;state,nstate
 +  maxiter,20
 +  orbital,orbrec
 +  ciguess,cirec
 +  save,det=cirec
 +  diab,diabrec,save=orbrec
 +  cpmcscf,grad,stateind,save=gradrec0;
 +  cpmcscf,nacm,stateind,stateind2(1),save=gradrec(1);
 +  cpmcscf,nacm,stateind,stateind2(2),save=gradrec(2);
 +  }
 +elseif (numcoup.eq.3.0) then
 +  {multi
 +  pspace,10.0
 +  wf;state,nstate
 +  maxiter,20
 +  orbital,orbrec
 +  ciguess,cirec
 +  save,det=cirec
 +  diab,diabrec,save=orbrec
 +  cpmcscf,grad,stateind,save=gradrec0;
 +  cpmcscf,nacm,stateind,stateind2(1),save=gradrec(1);
 +  cpmcscf,nacm,stateind,stateind2(2),save=gradrec(2);
 +  cpmcscf,nacm,stateind,stateind2(3),save=gradrec(3);
 +  }
 +elseif (numcoup.ge.4.0) then
 +  text, Number of couplings requested: $numcoup beyond limits of input deck.
 +  text, Please modify input deck accordingly.
 +  STOP
 +endif
 +forces;samc,gradrec0
 +}
 +
 +MGCOUP={
 +forces;nacm,gradrec(ncoup)
 +}
 +
 +! writes a molden file
 +WRITEMOLDEN={
 +moldenfile='$itrj.$itimestep.molf'
 +put,molden,$moldenfile;
 +}
 +
 +AIMS,dir='./aims_ethyl_CASSCF'
 +</code>
 +
 +<code - examples/aims_ethyl_CASSCF/Control.dat>
 + &control
 +
 +!------ Initialization -------!
 + NumStates=3            !Number of electronic states
 + InitState=2            !Initial trajectory's state
 + NumParticles=6         !Number of atoms
 + IMethod=1              !CASSCF dynamics with couplings
 + IRestart=0             !Do not restart
 + IRndSeed=1             !Random number seed
 + InitialCond='WIGNER'   !Sample from Wigner Distribution
 +
 +!-------- Integration and numerics --------!
 + TimeStep=20.0          !Simulation timestep
 + SimulationTime=800.0   !Simulation length - note this is short for illustrative purposes
 + ExShift=77.8381            !Diagonal shift of Hamiltonian
 + CoupTimeStep=5.0       !Timestep in coupling regions
 + MinTimeStep=1.0        !Minimum adaptive timestep
 +
 +!-------- Spawning -------!
 + CSThresh=3.0           !Coupling threshold to enter spawning mode
 + PopToSpawn=0.01        !Minimum population before spawning
 + OMax=0.6               !Prevent spawning if spawned trajectory has at least
 +                        !this much overlap with existing basis
 + MaxTraj=50             !Max number of trajectories
 + IgnoreState=1          !Remove uncoupled trajectories on S0
 + MaxEDiff=0.04          !Energy gap threshold for calculating couplings
 +
 +!------ Output Options -----!
 + zEDatFile=.true.       !Energy file
 + zNDatFile=.true.       !Population file
 + zTrajFile=.false.      !Trajectory data file
 + zAmpFile=.true.        !Amplitude files
 + zPotEnFile=.true.      !Potential energy files
 + zXYZ=.true.            !XYZ files
 + zCoupFile=.true.       !Coupling files
 + zChargeFile=.true.     !Charge files
 + zDipoleFile=.true.     !Dipole files
 + zTDipoleFile=.true.    !Transition dipole files
 + WriteMolden=.true.     !Write molden files?
 + MoldenStep=100         !time between molden files
 +/
 + 
 +</code>
 +
 +<code - examples/aims_ethyl_CASSCF/Geometry.dat>
 +       UNITS=BOHR
 +        6
 +              .0000000059         .0000000000        1.2790603970
 +             -.0000000025         .0000000000       -1.2790603970
 +             -.0000000032        1.7255896200       -2.3437915934
 +             -.0000000180        1.7255896200        2.3437915934
 +             -.0000000300       -1.7255896200        2.3437915934
 +              .0000000089       -1.7255896200       -2.3437915934
 +</code>
 +
 +<code - examples/aims_ethyl_CASSCF/Frequencies.dat>
 +                12
 +      903.7528062134      956.3515160547      973.2169016266     1081.4111041520
 +     1347.2742737477     1428.7725115327     1619.2543810461     1733.9363838709
 +     3311.5802800265     3331.2669071062     3396.3102184138     3424.0982861930
 +         .0000000013        -.0390884910         .0000000000        -.0000000013
 +        -.0390884910         .0000000000        -.0000000011         .2328927096
 +         .4297450592        -.0000000022         .2328927096        -.4297450592
 +        -.0000000083         .2328927096         .4297450592         .0000000116
 +         .2328927096        -.4297450592
 +        -.1224921629        -.0000000003        -.0000000019         .1224921629
 +        -.0000000003         .0000000009        -.3982789485         .0000000028
 +         .0000000065         .3982789485         .0000000046        -.0000000004
 +         .3982789485        -.0000000010         .0000000115        -.3982789485
 +         .0000000009        -.0000000053
 +         .0773477556        -.0000000004         .0000000014         .0773477556
 +         .0000000004         .0000000002        -.4608448140        -.0000000009
 +        -.0000000015        -.4608448140        -.0000000010        -.0000000078
 +        -.4608448140         .0000000036        -.0000000077        -.4608448140
 +        -.0000000017        -.0000000014
 +         .0000000000        -.0000000011         .0000000006         .0000000000
 +         .0000000008        -.0000000006        -.4980193316        -.0000000032
 +        -.0000000070         .4980193316         .0000000021         .0000000059
 +        -.4980193316         .0000000047        -.0000000049         .4980193316
 +        -.0000000005         .0000000059
 +        -.0000000001        -.1223605066         .0000000000        -.0000000001
 +         .1223605066         .0000000000         .0000000006        -.1174845526
 +        -.3808090966        -.0000000042         .1174845526        -.3808090966
 +         .0000000059         .1174845526         .3808090966         .0000000012
 +        -.1174845526         .3808090966
 +        -.0000000002         .0000000000         .1383677568         .0000000005
 +         .0000000000        -.1383677568         .0000000002        -.1068007050
 +        -.3500649362        -.0000000021        -.1068007050         .3500649362
 +         .0000000013         .1068007050         .3500649362        -.0000000032
 +         .1068007050        -.3500649362
 +         .0000000007         .0000000000        -.0652027245         .0000000004
 +         .0000000000        -.0652027245        -.0000000016         .2679058162
 +         .3884836376        -.0000000067        -.2679058162         .3884836376
 +        -.0000000048         .2679058162         .3884836376         .0000000003
 +        -.2679058162         .3884836376
 +        -.0000000004         .0000000000         .1409332482        -.0000000003
 +         .0000000000        -.1409332482        -.0000000011         .2552526222
 +         .2540250265         .0000000050         .2552526222        -.2540250265
 +         .0000000053        -.2552526222        -.2540250265        -.0000000014
 +        -.2552526222         .2540250265
 +        -.0000000010         .0000000000         .0416086531        -.0000000002
 +         .0000000000         .0416086531        -.0000000005         .4198210670
 +        -.2479080594         .0000000046        -.4198210670        -.2479080594
 +         .0000000075         .4198210670        -.2479080594         .0000000025
 +        -.4198210670        -.2479080594
 +         .0000000009         .0000000000        -.0511876036        -.0000000001
 +         .0000000000         .0511876036        -.0000000002         .4140808652
 +        -.2468786768        -.0000000049         .4140808652         .2468786768
 +        -.0000000074        -.4140808652         .2468786768         .0000000023
 +        -.4140808652        -.2468786768
 +         .0000000002         .0662092557         .0000000000         .0000000002
 +        -.0662092557         .0000000000         .0000000002         .3992928311
 +        -.2499404332         .0000000054        -.3992928311        -.2499404332
 +        -.0000000079        -.3992928311         .2499404332        -.0000000027
 +         .3992928311         .2499404332
 +         .0000000002         .0667440272         .0000000000        -.0000000002
 +         .0667440272         .0000000000         .0000000004        -.3976668561
 +         .2516792378         .0000000053        -.3976668561        -.2516792378
 +        -.0000000078        -.3976668561         .2516792378         .0000000021
 +        -.3976668561        -.2516792378
 +</code>