Ab initio multiple spawning dynamics

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

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.

The AIMS program is controlled by the Molpro input deck and a number of external text files:

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

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

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 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 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 Hessian.dat and Frequencies.dat respectively) must be present. The format of FrequenciesMP.dat follows the Molpro output of a Frequencies calculation:

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

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.

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.

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.

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.

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 examples). Control.dat should then be modified for your system and desired simulation parameters (section Control.dat).

2) Next, the initial conditions should be defined. This involves providing an initial geometry in Geometry.dat (section Geometry.dat) and normal-mode frequencies and displacements in FrequenciesMP.dat (section 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.

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

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

examples/aims_ethyl_HF.inp
! FMS-Molpro does not support parallel execution.
if(NPROC_MPPX.gt.1) then
skipped
endif
if(.not.modul_aims) then
skipped
endif

SYMMETRY,nosym
ORIENT,noorient
GEOMTYP=xyz
GEOMETRY={
6
ethylene ground-state geom
 C    0.0000000031    0.0000000000    0.6768496105
 C   -0.0000000013    0.0000000000   -0.6768496105
 H   -0.0000000017    0.9131426983   -1.2402810929
 H   -0.0000000095    0.9131426983    1.2402810929
 H   -0.0000000159   -0.9131426983    1.2402810929
 H    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'
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
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

Ethylene CASSCF Spawning Dynamics

See examples/aims_ethyl_CASSCF

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
 C    0.0000000031    0.0000000000    0.6768496105
 C   -0.0000000013    0.0000000000   -0.6768496105
 H   -0.0000000017    0.9131426983   -1.2402810929
 H   -0.0000000095    0.9131426983    1.2402810929
 H   -0.0000000159   -0.9131426983    1.2402810929
 H    0.0000000047   -0.9131426983   -1.2402810929
}
basis=6-31g*


!Aims variables - do not change!
! FMS-Molpro does not support parallel execution.
if(NPROC_MPPX.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,ci=cirec
  if (status.lt.0) then
    text, MCSCF failed in MPINIT.
    text, Attempting uncoupled iterations.
    {multi
    pspace,10.0
    wf;state,nstate
    orbital,orbrec
    ciguess,cirec
    save,ci=cirec
    {iterations
      do,uncouple,1,to,15}
    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,ci=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,ci=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,ci=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,ci=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,ci=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,ci=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'
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
/
examples/aims_ethyl_CASSCF/Geometry.dat
       UNITS=BOHR
        6
       C        .0000000059         .0000000000        1.2790603970
       C       -.0000000025         .0000000000       -1.2790603970
       H       -.0000000032        1.7255896200       -2.3437915934
       H       -.0000000180        1.7255896200        2.3437915934
       H       -.0000000300       -1.7255896200        2.3437915934
       H        .0000000089       -1.7255896200       -2.3437915934
examples/aims_ethyl_CASSCF/Frequencies.dat
         6        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