NEWTON-X a package for Newtonian dynamics close to the crossing seam

2 About Newton-X

Newton-X ¹^,² is a general purpose program package for molecular dynamics in the electronic excited states, including mixed quantum-classical methods (surface hopping).

Newton-X modular development allows it to be easily linked to any quantum chemistry package that can provide energy gradients and (optionally) non-adiabatic coupling vectors.

In the current version, Newton-X can perform dynamics using Columbus ³^,⁴, Turbomole ⁵, Dftb ⁶, DFTB+, and Gaussian 03 ⁷ program packages. Nonadiabatic dynamics using hybrid gradients (including QM/MM approach) is available as well for combinations between Tinker ⁸ and one of the following programs Columbus, Turbomole, DFTB+, and analytical models. Other third-party programs providing energies and nonadiabatic couplings can be easily integrated to work with Newton-X as well.

Newton-X code is distributed free of charge for non-commercial and non-profit uses. You may use the program freely and adapt the code to your needs. Please, if you have any enhancements, share them with us. You are, however, not allowed to re-distribute code or binaries in parts or total. Anyone intending to use Newton-X must contact us.

Newton-X is shipped to the user "as is". There is no guarantee that it will work as described. Usage is at your own risk, there is no liability taken for any damage or loss of data and/or money. If you experience problems please tell the developers supplying the version-number – they are always grateful for any hint, but bear in mind that there is no kind of official support for Newton-X.

Contact:

Mario Barbatti

Max-Planck-Institut für Kohlenforschung

Kaiser-Wilhelm-Platz 1

D-45470 Mülheim an der Ruhr, Germany

barbatti@kofo.mpg.de

www.mpi-muelheim.mpg.de/private/barbatti

3 Mixed quantum-classical dynamics simulations

Mixed quantum-classical approaches⁹ are probably the most employed class of methods to perform excited-state molecular dynamics simulations including non-adiabatic effects. In these approaches, which include the surface hopping and the Ehrenfest methods, the nuclear time evolution is treated classically by means of the Newton’s equations, while the time evolution of the population of each electronic state is treated separately. In the surface hopping method, the time evolution of the population is obtained in two steps: first, non-adiabatic transition probability between each pair of states is computed and a stochastic algorithm is applied to decide in which state the classical trajectory is propagated in the next time step. Second, statistics over a large set of independently computed trajectories allows getting the fraction of trajectories (occupation) in each state as a function of time. The main hypothesis behind the surface hopping approach is that the occupation and the quantum population of each electronic state are the same if an infinite number of trajectories are computed.¹⁰

There are several proposed ways to evaluate the non-adiabatic transition probabilities, since the most simple methods which just assume that the probability is the unity if the energy gap between the states is smaller than some threshold, to more sophisticated approaches which take into account the variation of wavefunction coefficients¹¹ or compute the Landau-Zener transition probability.¹² The most reliable procedure to compute the non-adiabatic transition probability for surface hopping simulations is the Tully’s fewest switches algorithm.¹⁰ In this approach, the time-dependent Schroedinger equation (TDSE) is integrated simultaneously to the classical trajectory.¹³ To cope with the lack of non-local information introduced by the on-the-fly approach, non-local terms in the TDSE are neglect and the nuclear wavefunction is supposed to be entirely localized at the classical position determined by the Newton’s equation. The integration of this semi-classical version of the TDSE gives the adiabatic populations of the electronic states, which are then used to compute the probability using the fewest-switches formula.

The integration of the TDSE depends on nonadiabatic coupling terms connecting different states. If adiabatic representation is used to expand the molecular wavefunction, non-adiabatic coupling vectors should be computed. Alternatively, if diabatic representation is used, non-diagonal Hamiltonian matrix elements should be computed. Either way, the computation of the non-adiabatic coupling terms are the bottleneck for non-adiabatic dynamics approaches. These terms are not usually available for most of quantum chemical methods and when they are, their computational cost increases with the square of the number of electronic states.¹⁴ These difficulties have motivated, on one hand, the search for approximated hopping algorithms as those mentioned above, and on the other hand the computation of the coupling terms based on wavefunction overlaps.¹¹^,^14-18 Besides that, it has been an important achievement the development of procedures for analytical computation of them at the multireference configuration interaction (MRCI) and at the multiconfigurational self-consistent field (MCSCF) methods.¹⁹^,²⁰

One consequence of the hyperlocalization of the nuclear wavefunction in mixed quantum-classical approaches is that non-diagonal terms in the density matrix do not vanish with time as they should do.²¹^,²² In surface-hopping, this results in an excessive number of hopping events from lower to upper states, which disturbs the accomplishment of the occupation/population hypothesis mentioned above. Decoherence can be imposed by applying an ad hoc correction to the adiabatic population every time step, which forces the non-diagonal terms in the density matrix to damp to zero within a certain time constant.²²

When a hopping between two states takes place, it usually does through finite energy gap. In order to keep the total energy constant in the subsequent trajectory, it is necessary to correct the kinetic energy, for example, by rescaling the momentum or by adding more momentum at the direction of the non-adiabatic coupling vector.¹³^,²³ It may also happen that the stochastic algorithm attempts to make a hopping from a lower to an upper state in a region where there is not enough energy to do so. Such cases have been usually treated by forbidding the hopping occurrence.¹³ The momentum can be kept or reversed afterwards. Another possibility is to take the time uncertainty principle to search for a geometry nearby where the hopping is allowed.²⁴

Because of the stochastic nature of the fewest-switches surface-hopping approach, trajectories starting with the same initial conditions will give rise to different time development. Moreover, the initial conditions should reflect the initial phase space distribution. Therefore the averages that define the state occupation should in principle be performed over this double ensemble of trajectories starting in different points of the phase space, several times in each one. Because of computational limitations, this procedure is usually reduced to a single ensemble of trajectories starting in different points of the phase space only once in each one.

The initial condition ensemble can be generated in a diversity of ways. For instance, the simulation of an instantaneously excited wave packet into the Franck-Condon region may be done by selecting geometries and velocities from a dynamics in the grounds state, regarding this dynamics run for long enough period as to allow an adequate sampling of the phase space. Alternatively, each nuclear degree of freedom can be treated within the harmonic approximation and a Wigner distribution can be build.

Most of the methods and algorithms mentioned in this short introduction are implemented in the Newton-X program.

4 Developers, Collaborators and Contributors

The Newton-X program has been developed in a multi-institutional collaboration, involving researchers from several countries. The project in headed at the Institute for Theoretical Chemistry of the University of Vienna.

People working in the general Newton-X development are:

Mario Barbatti MPI Mülheim Germany

Hans Lischka University of Vienna Austria

Matthias Ruckenbauer University of Vienna Austria

Felix Plasser University of Vienna Austria

Many other people have contributed in the past or are still contributing to development of specific algorithms in Newton-X. They are:

Surface hopping routines and initial condition distributions

Giovanni Granucci University of Pisa Italy

Maurizio Persico University of Pisa Italy

Time-dependent non-adiabatic couplings (MRCI and TD-DFT)

Jiri Pittner J. Heyrovsky Institute Czech Republic

Random initial conditions

Bernhard Sellner University of Vienna Austria

Andersen thermostat

Vladimir Lukeš Slovak University of Technology Slovak Republic

5 How to reference Newton-X

Please, cite Newton-X as:

M. Barbatti, G. Granucci, M. Persico, M. Ruckenbauer, M. Vazdar, M. Eckert-Maksić and H. Lischka, J. Photochem. Photobio. A 190, 228 (2007).

M. Barbatti, G. Granucci, M. Ruckenbauer, F. Plasser, J. Pittner, M. Persico, H. Lischka, Newton-X: a package for Newtonian dynamics close to the crossing seam, version 1.2, www.newtonx.org (2011).

References to specific methods, algorithms and third-party programs used in Newton-X are given along this documentation. For a summary, see Chapter 7.

6 Quick start

6.1 Newton-X Tutorial

A Newton-X tutorial with step-by-step procedures for several examples is available at www.univie.ac.at/newtonx/docs/tutorial.pdf

6.2 The NXINP tool

Most of Newton-X inputs are prepared with nxinp program. To execute nxinp, type:

$NX/nxinp

The first screen should looks like:

============================================================

NEWTON-X

Newtonian dynamics close to the crossing seam

www.newtonx.org

============================================================

MAIN MENU

1. GENERATE INITIAL CONDITIONS

2. SET BASIC INPUT

3. SET GENERAL OPTIONS

4. SET NON-ADIABATIC DYNAMICS

5. GENERATE TRAJECTORIES AND SPECTRUM

6. SET STATISTICAL ANALYSIS

7. EXIT

Select one option (1-7):_

The next sections will guide you through each one of these options. By now, it is enough to note that nxinp is self-explaining. When you select one of the options, say, option 2 (Set basic input), you are asked about a series of parameters. A sequence of input may be, for example:

============================================================

NEWTON-X

Newtonian dynamics close to the crossing seam

www.newtonx.org

============================================================

SET BASIC OPTIONS

nat: Number of atoms.

There is no value attributed to nat

Enter the value of nat : 12

Setting nat = 12

nstat: Number of states.

The current value of nstat is: 2

Enter the new value of nstat : 3

Setting nstat = 3

nstatdyn: Initial state (1 - ground state).

The current value of nstatdyn is: 2

Enter the new value of nstatdyn :

Setting nstatdyn = 2

dt: Time step for the classical equations.

The current value of dt (fs) is: 0.5

Enter the new value of dt (fs) : 0.1

Setting dt = _

Each parameter contains a short description and, most of time, an attributed default value. To use the default values, just press <ENTER>. More information about each parameter can be found in this documentation.

6.3 Initial conditions generation

Input

1- Prepare geom file with equilibrium geometry. Use xyz2nx to create geom from xyz files.

2- Prepare force.out file with the harmonic frequencies using, for example, Turbolome.

3- Run nxinp program. Select option 1. GENERATE INITIAL CONDITIONS. nxinp will help you to select the input parameters to control de initial condition generation. Exit nxinp.

Run

$NX/initcond.pl > initcond.log

Output

final_output: Initial conditions.

initcond.log: Log file.

Further options

Optionally, you can control the number of excitation quanta in each vibrational modes by providing a file qvector (see Chapter 10) together with the other input files.

It is also possible to pick points from previous dynamics calculations or to generate random velocities to be used as initial conditions.

6.4 Dynamics Input

6.4.1 Adiabatic dynamics (dynamics on one surface)

Input

1- Create a directory called JOB_AD containing a set of input files for geometry optimization (1 step) with the chosen electronic structure program.

2- Run nxinp program. First, select option 2. SET BASIC INPUT. nxinp will help you to select the input parameters to control the dynamics. After that, select option 3. SET GENERAL OPTIONS if you want to change some more technical parameter. Exit nxinp.

Hint

Check control.dyn file. If the keyword “thres” is defined there, be sure that its value is 0 (thres = 0).

6.4.2 Non-adiabatic dynamics (Surface hopping)

Input

1- Create a directory called JOB_NAD containing an input for non-adiabatic coupling (single point) with the chosen electronic structure program..

2- Run nxinp program ($NX/nxinp). First, select option 2. SET BASIC INPUT. nxinp will help you to select the input parameters to control the dynamics. After that, select option 3. SET GENERAL OPTIONS if you want to change some more technical parameter. Optionally, select option 4. SET NON-ADIABATIC DYNAMICS to change the non-adiabatic-dynamics options. Exit nxinp.

Hint

Check control.dyn file. If the keyword “thres” is defined there, be sure that its value is 100 (thres = 100).

6.4.3 Mixed adiabatic and non-adiabatic dynamics

Input

1- Create a directory called JOB_AD containing a set of input files for geometry optimization (1 step) with the chosen electronic structure program.

2- Create a directory called JOB_NAD containing an input for non-adiabatic coupling (single point).

3- Run nxinp program ($NX/nxinp). First, select option 2. SET BASIC INPUT. nxinp will help you to select the input parameters to control the dynamics. Optionally, select option 3. SET GENERAL OPTIONS if you want to change some more technical parameter. Optionally, select option 4. SET NON-ADIABATIC DYNAMICS to change the non-adiabatic-dynamics options. Exit nxinp.

Hint

Check control.dyn file. The keyword “thres” must be defined there. Its value is the energy difference threshold (eV) bellow to which the non-adiabatic dynamics starts. This option is not fully tested.

6.5 Creating the trajectory inputs

Input

1- Copy final_output file created after the initial conditions generation, section 6.3, into the same directory containing all files and directories created in the dynamics input, section 6.4.

2- If the jobs will run in a batch system, also copy the submission script file to that directory. In $NX/../batch you may find several examples of submission scripts.

3- Run nxinp program ($NX/nxinp). Select option 5. GENERATE TRAJECTORIES AND SPECTRUM. nxinp will help you to select the input parameters to control the trajectory inputs generation. At the end of the input selection, Newton-X will automatically generate the trajectory directories. This process can take some few minutes.

Output

The trajectories inputs were written to TRAJECTORIES/TRAJn, where n is the trajectory number.

6.6 Running the dynamics

1- Go to each TRAJECTORIES/TRAJn (see section 6.5) and run

$NX/moldyn.pl > moldyn.log &

or, if is this the case, submit the job to the batch system.

2- Alternativelly, go to TRAJECTORIES and run

$NX/submit.pl

It will allow you to automatically submit several sequential jobs to the batch system.

6.7 Where are the results?

1- During or after the dynamics, go to directory TRAJECTORIES/TRAJn/RESULTS and run:

$NX/plot to generate "energy x time" graph with Gnuplot.

molden dyn.mld to see motion with Molden or any visualization package (xyz format).

$NX/arrow to generate a Molden file for a specific time step, containing the velocity and, in the case of dynamics with Columbus, also molecular orbitals and non-adiabatic coupling vectors.

2- dyn.out file contains details about the geometry, velocity, energy and wave function (adiabatic coefficients) along the trajectory. tprob contains information about the hopping probability at each time step. en.dat contains information about the energy of each state. sh.out contains further information about the TDSE integration.

3- In TRAJECTORIES/TRAJn, the standard output (moldyn.log) contains the log information of the job and information about states, gradients and non-adiabatic couplings. The standard output is also written to RESULTS/nx.log.

4- In TRAJECTORIES/TRAJn/DEBUG, runnx.error contains occasional error messages. log.conv contais the information about convergence of ab initio calculations.

5- TRAJECTORIES/TRAJn/INFO_RESTART contains a complete set of files to restart the dynamics from the last time step that run.

6.8 Overview of the file structure

6.8.1 Dynamics simulations

The basic structure of directories and files during the dynamics simulations is shown in Figure 1.

In the case of hybrid calculations (QM/MM) the structure is very similar. The main difference is that in this case the JOB_AD and JOB_NAD directories have a substructure specifying the several jobs. The inset in Figure 1 shows an example of this substructure for a QM/MM job with Columbus and Tinker.

tree_3.tif

Figure 1 Basic structure of files and directories for a dynamics simulation. The inset shows an example for QM/MM simulations.

7 Capabilities

7.1 General features

Features	Options
Dynamics	On-the-fly dynamics with the velocity-Verlet algorithm ²⁵^,²⁶
	Mixed quantum classical (non-adiabatic) dynamics ⁹ Surface hopping fewest switches algorithms ²⁵^,²⁶ Dynamics with non-adiabatic coupling vectors and time-derivative couplings ¹⁴^,¹⁵ Dynamics with Local Diabatization method ¹⁶ Decoherence corrections ²² Lagrangean extrapolation of orbitals ²⁷ Damped dynamics (kinetic energy always null) Andersen thermostat ²⁸^,²⁹ Hybrid dynamics (QM/MM) ³⁰
Interfaces	Columbus (MCSCF,MRCI) ³^,⁴ Turbomole (TD-DFT, RI-CC2, ADC(2))⁵ Dftb (TD-DFTB) ⁶ Dftb+ (DFTB) Tinker ⁸ Gaussian (CASSCF) ⁷ Dft-mrci ³¹
Initial conditions	Wigner distribution
	Quantum and classical harmonic oscillator distributions
	Pick points from previous dynamics
	Random velocity generation ³²
Spectrum simulation	UV photoabsorption cross section^33-35 UV absorption and emission spectra.
File management	Friendly input via nxinp facility
	Automatic management of files and directories in multiple trajectories
Output and analysis	Statistical analysis of results (internal coordinates and forces, energies, wave function properties)
	Normal Mode Analysis, Essential Dynamics Analysis ³⁶^,³⁷
	Recomputation of internal coordinates
	On-the-fly graphical outputs (Molden, Gnuplot)

7.2 Overview of the available interfaces

			Initial conditions / Spectrum		Dynamics
Program	Version	Method	Frequency reading	Spectrum simulation	Adiabatic	Nonadiabatic (surface hopping)
						NAC	CIO	LD
Columbus	5.9	MRCI / MCSCF
	7	MRCI / MCSCF
Turbomole		TD-DFT
		CC2 / ADC(2)
Dftb		TD-DFTB
Dftb+		DFTB
Gaussian	03	CASSCF
Tinker		MM
Dft-mrci		DFT-MRCI
Molden		-
Analytical		User defined

Available features

NAC – Based on nonadiabatic coupling vectors.

CIO – Based on wavefunction overlaps.

LD – Based on local diabatization.

8 How to get Newton-X

Newton-X code is distributed free of charge for non-commercial and non-profit uses. To request a copy, send an e-mail to barbatti@kofo.mpg.de containing your name and affiliation.

The third-party quantum chemistry and visualization programs are not distributed with Newton-X. They should be directly obtained from the respective distributors and owners.

9 How to install Newton-X

9.1 To install and run Newton-X you need

· The Newton-X source code. See Chapter 8 to see how to get the Newton-X source code.

· A FORTRAN 90 compiler installed, preferentially the Intel Fortran compiler.

· PERL 5.8.0 or higher installed (www.perl.com).

· BLAS, LAPACK, and GSL libraries installed. More information on how to set the paths to these libraries can be found in the README.install file, which can be found in the Newton-X distribution.

· At least one of the interfaced third-party programs (Columbus, Turbomole, Dftb, Tinker, Gaussian). See Chapter 7 of this documentation for an updated list of available interfaces.

For the full functionality, you also need:

· Gnuplot (www.gnuplot.info).

· Some molecular visualization program able to read xyz files (e.g., Molden, Vmd or Molekel).

9.2 Installation procedure

Untar and uncompress the source code file

tar –zxf nx-<version>.tgz

or, alternatively,

gunzip nx-<version>.tgz

tar –xf nx-<version>.tar

The files will be written to the NX-pack directory. Go to the NX-pack/install subdirectory and run the script

./installnx.pl

It will ask you for an installation address and for which packages you want to install. If the installation proceeds without problem, you will get in the screen something like:

================================================

Welcome to the NEWTON-X installation routine

www.newtonx.org

================================================

Please enter directory for the NEWTON-X installation

/home/users/NX ! Give the path in your system

This directory seems not to exist - create it? (y/n) y

Created directory /home/users/NX

Select packages to be installed.

[0] Complete installation [Default option]

[1] Copy files and install tools (TOOLS)

[2] Initial conditions (TOOLS+INITCOND)

[3] Molecular dynamics (TOOLS+MOLDYN+MODEL)

[4] Time-dependent overlaps (TOOLS+CIOVERLAP)

Enter comma separated list or dash separated range: 0 ! Complete installation

Starting NEWTON-X installation at Mon May 17 09:27:16 CEST 2010

Compiler keyword: ifort

Flag keyword: -static -vec-report0

Program version: 1.2.0a (01/Jun/2010)

Fortran version: /opt/intel/Compiler/11.1/069/bin/intel64/ifort

Perl version: /usr/bin/perl

Original directory: /home/users/NX-pack/install

Installation directory: /home/users/NX

BLAS: /opt/intel/mkl/10.1.3.027/lib/em64t/

GCC for LA libraries: /usr/bin/g++

GCC for cioverlap: /usr/bin/g++

BLAS for cioverlap: /usr/lib64

LAPACK for cioverlap: /usr/lib64

GSL for cioverlap: /usr/lib64

GSLCBLAS for cioverlap: /usr/lib64

Packages selected:

1 TOOLS

2 INITCOND

3 MOLDYN

4 CIOVERLAP

Checking compilers, libraries and dependencies...

No problem detected for the selected packages.

Copying file .......................................................................... ..........................DONE

Compiling INITCOND ........................DONE

Compiling MOLDYN ................................................................DONE

Compiling TOOLS ................................DONE

Compiling MODEL ............DONE

Compiling CIOVERLAP ................................DONE

Setting permissions ........DONE

Installation completed!

Don't forget to set the environment variable $NX to /home/users/NX/bin

Do not forget to set the shell-variable $NX to <install-path>/bin after installation.

For c-shell users, this is done by:

setenv NX <install-path>/bin

For bash-shell users, this is done by:

export NX=<install-path>/bin

A complete installation log is written to installnx.log file.

As default, instalnx.pl is set to use the Intel Fortran Compiler for Linux. To change the compiler and the compiler options follow the instructions contained in the installnx.pl body. More information on compilers can be found in the README.install file, which can be found in the Newton-X distribution. Newton-X has been compiled and tested in 32-bit, EM64T, and Opteron machines.

Newton-X works using third-party programs, such as Columbus, Turbomole, Gaussian, and others. These programs should be separately obtained from their distributors and installed according to their specific instructions.

Newton-X assumes that the variables $COLUMBUS, $DFTB, $DFTPLUS, $g03root, $DFTCI pointing to, respectively, Columbus, Dftb, Dftb+, Gaussian 03, and Dft-mrci directories are defined in the system. Turbomole and Tinker are assumed to be in the default $PATH. For Dft-mrci, the variable $KSCIHOME pointing to the directory containing the parameters files must also be defined.

9.3 Problems with the Installation

The installation routine will check the libraries and dependencies related to the chosen packages. If something is wrong or missing, the installation output will look like:

================================================

Welcome to the NEWTON-X installation routine

www.newtonx.org

================================================

Please enter directory for the NEWTON-X installation

/home/users/NX ! Give the path in your system

This directory seems not to exist - create it? (y/n) y

Created directory /home/users/NX

Select packages to be installed.

[0] Complete installation [Default option]

[1] Copy files and install tools (TOOLS)

[2] Initial conditions (TOOLS+INITCOND)

[3] Molecular dynamics (TOOLS+MOLDYN+MODEL)

[4] Time-dependent overlaps (TOOLS+CIOVERLAP)

Enter comma separated list or dash separated range: 0 ! Complete installation

Starting NEWTON-X installation at Mon May 17 09:27:16 CEST 2010

Compiler keyword: ifort

Flag keyword: -static -vec-report0

Program version: 1.2.0a (01/Jun/2010)

Fortran version: /usr/users/barbatti/intel/Compiler/11.1/069/bin/intel64/ifort

Perl version: /usr/bin/perl

Original directory: /ns80th/nas/users/barbatti/NX-DEVELOPMENT/NX-1.2.0a/install

Installation directory: /ns80th/nas/users/barbatti/TEST_NX/NX-1.2.1a

BLAS: /opt/intel/mkl/10.1.3.027/lib/em64t/ <<NOT FOUND>>

GCC for LA libraries: /usr/bin/g++

GCC for cioverlap: /usr/bin/g++

BLAS for cioverlap: <<NOT FOUND>>

LAPACK for cioverlap: <<NOT FOUND>>

GSL for cioverlap: <<NOT FOUND>>

GSLCBLAS for cioverlap: <<NOT FOUND>>

SIFS integer format: 64 bits

Packages selected:

1 TOOLS

2 INITCOND

3 MOLDYN

4 CIOVERLAP

Checking compilers, libraries and dependencies...

Possible problems were detected.

Detailed information was written to installnx.log.

If the installation proceeds, it will probably fail.

Do you want to continue the installation? (y/n) [Default=y] (default=y) n

Installation aborted.

In this specific example, the following problems occurred:

1) The path to BLAS library defined in installnx.pl was not found.

2) The BLAS, LAPACK, and GSL libraries used in the CIOVERLAP package were not found.

The installnx.log file will contain more detailed information about each of these problems. For instance, for the first problem, it will issue the following warning message:

...

===================================================================

Checking compilers, libraries and dependencies...

*** WARNING ***

Cannot find $BLAS libraries as defined in the beginning of installnx.pl.

Check whether the $BLAS value is appropriated for your system.

If $BLAS=-lblas then it is possible that Blas libraries are not installed

in this system. If $BLAS has an absolute path ($BLAS=-L/...) then, check

whether the path is valid. When using MKL libraries (INTEL), the absolute path

will change from version to version. You may have to change the $BLAS value

in the beginning of installnx.pl.

...

9.4 Verification of the installation

To test the installation, create a directory to execute the tests. Move into it and select the tests that you want to perform by running

$NX/inp-testnx.pl

It will return a list of all available tests. You can select particular tests or all of them. It will create a file called test.inp containing a space-separated list with the number of the tests that should run.

Then, run the program

$NX/test-nx.pl > test-nx.log &

This program will run the selected tests. If test.inp is not present, all tests are selected by default. The tests are a series of quick pre-build examples. It will check whether they have normal termination or not, and it will compare the results to those of standard files. Check the results in test-nx.log.

9.5 Compatibility between DALTON and NEWTON-X installations

Even after a successful installation, NEWTON-X jobs using COLUMBUS may return the following error message:

readsifs returned with error! Are the integer formats between Columbus and NX matching?

In such cases, the probable cause is an incompatibility issue between DALTON program (used by COLUMBUS) and NEWTON-X. It can be solved by recompiling CIOVERLAP program, by the following these steps:

1. Go to install directory.

2. Edit installnx.pl file and change the value of variable $status_record32. The default value is “on” (SIFS integer format = 64 bits), which means that DALTON is assumed to have been compiled for 64 bits. “off” (SIFS integer format = 32 bits) means that DALTON is assumed to have been compiled for 32 bits.

3. Run the installation program again, but this time select option

[4] Time-dependent overlaps (TOOLS+CIOVERLAP).

4. It is not necessary to recompile LA library.

A previous installation of the LA library already exits.

Do you want to replace it? (y/n) n

10 Initial conditions generation

INITCOND is a set of programs developed to generate initial conditions to the molecular dynamics and spectrum simulation.

It is possible to generate initial conditions by sampling the data according to harmonic oscillator distributions (classical and quantum). It is also possible to pick at random points from a previous dynamics calculation performed with Newton-X or to generate random velocities for some specific geometry.

For NACT = 1, 2 or 3 the initial conditions are based on normal modes. The normal modes themselves should be generated in a separated run and given as input to Newton-X. Newton-X can read the normal modes directly from the output files of several third-party programs (see IPROG below). Despite the initial format, the normal modes are internally transformed into mass-weighted normal modes L. For NACT = 1 or 2, for each normal mode, an amplitude and a momentum are randomly selected from a harmonic oscillator distribution. The Cartesian coordinates and velocities of all atoms are then determined from the normal coordinates by the inverse transformation³⁷. If NACT = 3, only is randomly selected, while is scaled as to give the harmonic oscillator energy. For NACT = 2, if the vibrational numbers are zero, the distributions matches the Wigner distribution for the quantum harmonic oscillator ³⁸^,³⁹. For higher vibrational quantum numbers, the distribution for each normal mode is , where and are respectively the harmonic oscillator wavefunctions in the coordinate and momentum spaces.

If NACT = 4, random points are picked from a previously performed dynamics. If NACT = 5, random velocities are generated for a fixed velocity according to the algorithm described in Ref. ³². If NACT = 6, single point electronic structure calculations are performed for a previously computed set of initial conditions.

10.1 How to execute INITICOND

Run

$NX/initcond.pl > initcond.log

10.2 Input parameters

File: initqp_input

Namelist: DAT

Parameter	Default	Description
NUMAT	= [3]	number of atoms
NPOINTS	= [1]	number of initial conditions generated
NACT	= [2]	1 - Each normal mode contains a given energy in the initial state, according to the vibrational quantum number, and the distribution of coordinates and momenta is that of a classical harmonic oscillator. 2 - The distribution of coordinates and momenta is that of a quantum mechanical harmonic oscillator in the specified vibrational state. The sample of the coordinates and momenta is uncorrelated. In the ground vibrational state, this distribution matches the Wigner distribution for the harmonic oscillator. 3 - The distribution of coordinates and momenta is that of a quantum mechanical harmonic oscillator in the specified vibrational state. Only coordinates are sampled. For each normal mode, the momenta are scaled to the coordinates up to the have the correct harmonic vibrational energy. 4 - NPOINTS points are picked up from some previously run dynamics simulations. 5 - Gaussian-distributed random velocities are generated for some specific geometry. 6 - Single point electronic structure calculations for a previously computed set of initial conditions.
ISEED	= [1]	Random number seed 0 Use standard value for the random number seed. -1 Use random seed. Any interger > 0 is used as the seed itself. Use ISEED = -1 if the job will be split (see section 15.7).
LVPRT	= [1]	1 - Standard print level. 2 - Debug print level.
If NACT ≤ 4
N_PICK	= [-1]	Method to pick up points: -1 - pick random points from trajectories. n (integer) - pick time step n from trajectories.
NIS	= [1]	Intial state (Ground state = 1).
NFS	= [2]	Final state. (Information for all states between NIS and NFS are stored.)
CHK_E	= [0]	0 - Do not check the energies. 1 - Check the energies between states NIS and NFS. For ground state dynamics, set CHK_E = 0.
CMP_E	= [0]	(Only if NACT = 4 and CHK_E = 1) 0 - Read energies from dynamics output 1 - Compute energies.
PROG	= [1.0]	(Only if CHK_E = 1 and NACT = 3 or CMP_E = 1 and NACT = 4) Program to compute vertical excitation energies: 1.0 - Columbus 2.0 - Turbomole RI-CC2 / ADC(2) 2.1 - Turbomole TD-DFT 5.0 - Dftb 9.0 - Dft-mrci
EVERT	= [5.0]	Required vertical energy (eV).
DE	= [0.5]	Allowed variation EVERT (eV).
If NACT ≤ 3
FILE_GEOM	= [geom]	file containing the equilibrium geometry of the molecule in Columbus and Newton-X format. (See section 15.6 for conversion tools)
FILE_NMODES	= [nmodes]	file containing the normal coordinates and the vibrational frequencies (in cm^-1). Normally it is an output file from some quantum chemistry program (see option IPROG).
FILE_OUT	= [ini_qv]	output file from INIQP, containing NPOINTS set of generated Q and V.
FILE_VIB	= [qvector]	file with the vibrational quantum numbers. See description bellow.
ANH_F	= [1.0]	Multiply all frequencies by ANH_F. Useful to add anharmonic effects.
KVERT	= [1]	0 - Use the provided EVERT. 1 - Use Evert of the equilibrium geometry that is calculated in the first step of the program.
IPROG	= [1]	Read vibrational modes from: 1 - Gamess 2 - Turbomole 3 - Columbus 4 - Gaussian 5 - Molden 6 - Dftb 7 - Aces2
NM_FLAG	= [0]	(Only if IPROG = 5) If normal modes should be read from a Molden file, then the type of normal mode should be given. 1 - Cartesian normal modes (amu^-1/2) 2 - Normalized Cartesian normal modes 3 - Mass weighted normal modes The default 0 will terminate the program.
IF NACT = 4
ADDRESS	= []	Path to the set of previous trajectories from which the initial conditions must be picked. The default address [/home/old_dyn/TRAJECTORIES] certainly should be changed.
TRAJI	= [1]	Initial trajectory to be read.
TRAJF	= [10]	Final trajectory to be read.
TI	= [0]	(Only if N_PICK = -1) Start to pick points only after time TI (fs) of each trajectory.
TF	= [-1]	(Only if N_PICK = -1) Do not pick points after time TF (fs) of each trajectory. Diagnostic.pl (see section 15.5) is always called in this kind of run. The time suggested in its output (diag.log) is used if it is smaller than TF. If TF = -1, always use the time suggested in diag.log file.
REORDER	= [1]	(Only if N_PICK = -1) Sort the points according to the trajectory number and time step. 0 - Do not sort. 1 -Sort.
ETOT_DEV	= [0.5]	Allowed variation in the total energy (eV).
POP_DEV	= [0.1]	Allowed variation in the norm of the adiabatic population.
IF NACT = 5
FILE_GEOM	= [geom]	file containing the initial geometry of the molecule in the Columbus and Newton-X format. (See section 15.6 for conversion tools)
EKIN	= [0]	Kinetic energy (eV).
TEMP	= [0]	Temperature (K). If TEMP = 0 (default), this option generates random velocities corresponding to kinetic energy Ekin (microcanonical ensemble). If TEMP > 0, canonical ensemble of random velocities is generated with mean kinetic energy EKIN and standard deviation s = k_BT(3N/2)^1/2. In any case, translational and rotational velocities are zero.
IF NACT = 6
NIS	= [1]	Intial state (Ground state = 1).
NFS	= [2]	Final state. (Information for all states between NIS and NFS are stored.)
PROG	= [1.0]	Program to compute vertical excitation energies: 1.0 - Columbus 2.0 - Turbomole RI-CC2 / ADC(2) 2.1 - Turbomole TD-DFT 5.0 - Dftb 9.0 - Dft-mrci 20 – Hybrid energies (QM/MM)

Example of initqp_input file:

&DAT

nact= 2,

numat= 2,

npoints= 20,

file_geom= hf.geom,

file_nmodes=hf.nmodes,

file_out= qvector,

evert= 5.0,

de= 0.25,

kvert= 1,

chk_e= 1,

prog = 1.0,

iprog = 2.0,

nis = 1,

nfs = 2,

&END

File : columbus.par

Optional file containing information for Columbus jobs.

Parameter

Default

Description

MEM

[100]

Columbus core memory (1 GB = 134 Mwords).

File : turbomole.par

Optional file containing information for Turbomole jobs.

Parameter

Default

Description

PARALLEL

[1]

Number of cores to use for parallel Turbomole (SMP only, no MPI!)

File : dftb.par

Optional file containing information for Dftb jobs.

Parameter

Default

Description

DFTB_EXEC

= [dftb]

Name of the Dftb executable file. $DFTB variable must be defined in the system. Newton-X will run $DFTB/<dftb_exec>.

File : dftb+.par

Optional file containing information for DFTB+ jobs.

Parameter

Default

Description

DFTBP_EXEC

= [dftb+]

Name of the Dftb+ executable file. $DFTBPLUS variable must be defined in the system. Newton-X will run $DFTBPLUS/<dftb_exec>.

File : dftci.par

Optional file containing information for Dft-MRCI jobs.

Parameter

Default

Description

maxiter

= [6]

Maximum number of MRCI cycles.

10.3 What you need to execute

· initqp_input with the set of parameters. (It can be generated with nxinp tool.)

10.3.1 Additional files

If NACT ≤ 3:

· file with the equilibrium geometry (see FILE_GEOM keyword). Hint: the program tm2nx converts the coord file (Turbomole format) to the Newton-X format.

· file with the vibrational modes (see FILE_NMODES keyword).

· file with the vibrational quantum numbers (see FILE_VIB keyword). The default is the ground state (0 quantum in each mode). If you want to change this default values, write to FILE_VIB a list of quanta in each mode. Example:

0 0 0 0 0

0 1

This example puts one quantum at the highest frequency mode of a system with 12 modes (6 atoms). Keep the format with 5 columns of integers. The order of the modes is the same as in FILE_NMODE. If FILE_VIB des not exist, the program assumes 0 for all modes. To set -0.5 for a mode makes this mode contribute with 0 to the initial energy.

· JOB_AD directory containing input files for excited-state single point calculation with a third-party quantum chemistry program (only if CHK_E=1).

· save_file, optional file (see section 10.3.3).

If NACT = 4:

· The TRAJn directories of some dynamics calculations previously performed.

If NACT = 5:

· Geometry file with the initial geometry in Newton-X format.

If NACT = 6:

· A previously computed set of initial conditions in the Newton-X format renamed final_output.old.

· JOB_AD directory containing input files for excited-state single point calculation with a third-party quantum chemistry program.

· save_file, optional file (see section 10.3.3).

10.3.2 The JOB_AD directory

The JOB_AD directory contains input files for excited-state single point calculation with a third-party quantum chemistry program.

For the cases that this directory is needed, set the oscillator strength to be calculated by the quantum chemistry program. This will give you the possibility of generating UV spectra as well as to use the transition probability to select the initial conditions.

Specifically for Columbus calculations (PROG = 1.0): Input files must be prepared for a single point MCSCF or CI calculation in C₁ point group (no symmetry).

Specifically for Dftb calculations (PROG = 5.0): Input files must be prepared for excited-state spectrum calculation using code 10. The Dftb input parameter and geometry files must be named dftb.in and in.gen, respectively.

For Dftb+ calculations (PROG = 8.0) Input files must be prepared for a DFTB calculations yielding energy and gradient. The Dftb input parameter and geometry files must be named dftb.in and in.gen, respectively. At the present, only ground state dynamics (no TD) is available with this program.

For Dft-mrci calculations (PROG = 9.0): JOB_AD should contain a complete set of input files for Turbomole DFT calculation, including the auxiliary basis set, and the specific input file for the Dft-mrci program. The Dft-mrci input file must be named mrci.inp and should be set to the total number of excited states.

10.3.3 Save files option

When electronic structure calculations are performed during the initial condition generation (NACT ≤ 3 with CHK_E = 1 or NACT = 6), it is possible to give a list of files or directories that should be saved for every initial condition. This is simply done by creating a file name save_file and including a list of file names to be saved. save_file must contain only one file name per line. For instance, you may want to save the molecular orbitals (mos file) and the scf log file (dscf.out) resulting from the Turbomole calculations for each initial condition. In this case the following save_file should be given together with the other input files:

Example of save_file file:

dscf.out

mos

The required files are tar-compressed and written to the DEBUG directory renamed as ic<card>-filename.tgz, where <card> is the number of the initial condition as given in the final_output file. For the example above, DEBUG will contain:

ic0-dscf.out.tgz

ic0-mos.tgz

ic1-dscf.out.tgz

ic1-mos.tgz

ic2-dscf.out.tgz

ic2-mos.tgz

…

ic0-<NPOINTS>.out.tgz

ic0-<NPOINTS>.tgz

The routine that saves the files runs right after the execution of the third-party program, before the cleanup of execution directory.

If the file that should be saved is inside a directory, the relative path must be provided. For instance, if we want to save the MCSCF log file of Columbus jobs (mcscfsm.sp), save_file should contain “LISTINGS/mcscfsm.sp”. The full LISTINGS directory will be saved is save_file contains “LISTINGS”.

10.4 Output

The initial conditions are written to final_output file. The data from spectrum simulation are written to the cross-section.dat, spectrum.dat and spectrum-hist.dat files.

Initial condition for multiple states are written to final_output.[initial-state].[final-state]. Thus, for example, final_output.1.3 contains initial conditions for transitions from state 1 (ground state) to state 3. The final_output file is the same as final_output.[NIS].[NFS].

DEBUG directory contains the error log information (runnx.error) and files that should be saved when save_file file is provided (see section 10.3.3).

10.5 Running in several computers

It is possible to split the spectrum and initial condition generation job to run in different computers and merge them again afterwards. See section 15.7 for general instructions.

Do not forget to set ISEED = -1 when the jobs are split.

11 Trajectory-input generation, initial conditions for multiple states, and spectra

11.1 nxinp and options in the makedir.pl program

The section “5. GENERATE TRAJECTORIES AND SPECTRUM” of nxinp input program allows multiple different tasks, which can be selected in the sub menu:

============================================================

NEWTON-X

Newtonian dynamics close to the crossing seam

www.newtonx.org

============================================================

GENERATE TRAJECTORIES AND SPECTRUM

type: What do you want to do?

1 - Generate spectrum

2 - Select initial conditions for multiple initial states

3 - Generate trajectories

4 - Return to main menu

The current value of type is: 3

Enter the new value of type :

After selecting one of these options, you will be asked a series of question about the specificities of your job. At the end, the file mkd.inp will be generated and the program makedir.pl will be automatically running when you see the message:

Processing data: This may take some minutes. Please, wait...

If you want, after leaving nxinp, makedir.pl can be executed again by running:

$NX/makedir.pl > makedir.log

The options in makedir.pl may be changed by setting the following keywords in mkd.inp.

Name: mkd.inp

Parameter	Default	Description
TYPE	= [3]	1 - Generate spectrum 2 - Select initial conditions for multiple initial states 3 - Generate trajectories
nis	= [1]	Initial state. 1 is the ground state. (Spectrum simulation.)
nfs	= [2]	Array of final states (space separated, e.g., 2 3 4). If not only spectrum, but trajectories input is required as well, only one final state is allowed.
SCREEN	= [0]	Energy restriction: 0 - don't apply any restriction 1 - use the original energy restriction written in the final_output files 2 - apply new energy restriction
E_CENTER	= [0.0]	Center of the energy restriction (Only if SCREEN = 2) x - value of the center of restriction (eV) ref n - use the vertical excitation of final_output.nis.n file
E_VAR	= [0.5]	Width of the energy restriction. (eV) (Only if SCREEN = 2)
OS_CONDON	= [-1]	Oscillator strength: -1 - try to read from final_output file x - oscillator strength is always x (Condon approximation)
prob_kind	= [F]	Probabilities in the spectrum generation will be computed according to A - Einstein coefficients A (spontaneous emission) B - Einstein coefficients B (induced absorption or induced emission) F – Oscillator strength (photoabsorption cross section)
NORM	= [local]	Normalization of the Eintein's coefficients: local - Use energy-restricted data set global - Use complete data set
seed	= [0]	Seed for the random number generation 0 - a default random number seed is used 1 - a randomized seed is used Any other positive integer is used itself as the random number seed.
if TYPE = 1.The next three keywords define the parameters for the Gaussian broadening method.
L_SHAPE	= [lorentz]	Line shape: gauss – Normalized gaussian function lorentz – Normalized Lorentzian function
delta	= [0.05]	Phenomenological broadening of the spectrum using the Gaussian method ¹ (in eV). The default D corresponds to a small broadening, equivalent to 0.5 nm.
Temp	= [0]	Temperature correction for the spectrum (in K).
NREF	= [1]	Refraction index
eps	= [0.005]	Distance between consecutive points in the spectrum using the Gaussian method. (in eV)
kappa	= [3]	k (integer) is used to define the range of the spectrum between and .
TYPE = 3. The next three keywords define the batch system behaviour.
title	= [nx]	Title. Useful if the job runs in batch.
subfile	= [pmold]	If the job run in batch, subfile is the name of the submission script.
batchdef	= [-N]	Definition of the batch system. Default is SGE but Newton-X attempts to read it from SUBFILE and change it when some suitable value is found.

11.2 Transition spectra

Newton-X computes the spectrum by means of the line broadening method ¹^,³⁵, by assigning to each initial condition a line shape function g with the height P and a width representing some phenomenological broadening (D) and plotting the sum (S) of these line shape functions as a function of the transition energy E, i.e.,

E (eV) and S(E) (arbitrary unities) are written to spectrum.dat file. The line shape options are the Gaussian function

or the Lorentzian function

The intensity of each line P can be taken as the Einstein coefficient A and B or the oscillator strength f. Absolute values for the Einstein’s coefficients ⁴⁰ can be obtained from the log information written in makedir.log. First, note that the Einstein’s coefficients are given by

with and . As usual, e, h, m, c and e₀ are, respectively, fundamental charge, Planck’s constant, electron mass, speedy of light, and free-space permittivity. Therefore, to compute the absolute values of the Einstein’s coefficients, first check the units of vertical excitation energy in final_output file. If they are in eV, then multiply coefficients A by 0.43392×10⁸ or B by 0.11445×10¹⁵. If they are in hartree, then multiply coefficients A by 0.32130×10¹¹ or B by 0.56800×10¹⁰. The final units of A will be s^-1 and of B will be m³s^-2J^-1.

Note that the algorithm assumes the same degeneracy factor for both states (g_i = g_f) and refraction index n = 1. When spectrum involving multiple states is selected by attributing several values to the nfs array, final_output.nis.isf files (isf is each value in nfs array) for each transition must be present.

The file spectrum.dat is composed of two columns:

E (eV) S (arbitrary units)

When the intensity is selected to be f (PROB_KIND = F), in addition to the spectrum computed as above, the photoabsorption cross section is computed as well and written to the file cross-section.dat. The cross section is given by³⁵

where the internal sum runs over each of the NPOINTS initial coordinate R_k and the external sum runs over the several states (NFS). The factor g is given by³³

In order to compare the simulations to the experiments it useful to remind that photoabsorption cross sections (cm²) and molar extinction coefficients (e in M^-1cm^-1) can be interconverted by the relation⁴¹

where N_A is the Avogadro’s number.

The file cross-section.dat is composed of three columns:

E (eV) l (nm) s (Å².molecule^-1)

During spectrum and initial condition generation, the DE energies (eV) of all accepted initial conditions are written to spectrum-hist.dat file. A histogram produced with this data will also give the absorption spectrum. The graphical tool is available to produce simple histograms. From the directory containing spectrum-hist.dat, call $NX/spectrum.pl arg, where arg is an optional argument with the number of bins in the histogram. If arg is not given, the bin width is assumed to be 0.05 eV. The program produces a simple Gnuplot histogram that represents the spectrum. (Counts x eV) Average and standard deviation values are written to hist.out. Frequency table is written in hist.dat.

11.3 Initial conditions for multiple states

In Chapter 10 it was explained how to generate initial conditions for starting the dynamics in a single excited state. That procedure may be generalized to create initial conditions for multiple adiabatic states, weighting each one according to their oscillator strength.

Having the final_output.[nis].[nfs] files in the same directory, run nxinp and select option

2 - Select initial conditions for multiple initial states

The energy restriction may be changed or not. After running makedir.pl, a new directory called SELECTED_INITIAL_CONDITIONS is created. Inside this directory, there are new final_output.[nis].[nfs] files containing initial conditions for each nfs state. These files can be indivisually used to start the dynamics in each nfs state.

Suppose you want to generate initial conditions to start the dynamics simultaneously from S₂ and S₁ states. SELECTED_INITIAL_CONDITIONS directory should contain the files final_output.1.2 and final_output.1.3. The number of initial conditions in each one reflects the dipole transition probability from S₀ into these two states. If, for example, if S₁ state has pp* character and S₂ state has np* character, the number of points in final_output.1.2 will be substantially larger than in final_output.1.3. When the dynamics is performed the number of trajectories starting in S₁ and S₂ states should reflect this proportion.

11.4 Managing several trajectories

For managing several trajectories, it is useful to use the script makedir.pl. makedir.pl reads final_output file generated by the initial condition generation procedure, and writes a new Newton-X input directory for each initial condition accepted.

After preparing the inputs using nxinp tool, be sure that you have in the same directory:

· final_output (always)

· control.dyn (always)

· JOB_AD (if necessary for the specific job)

· JOB_NAD (if necessary for the specific job)

· jiri.inp (if necessary for the specific job)

· sh.inp (if necessary for the specific job)

· wf.inp (if necessary for the specific job)

· <third-party program>.par (if necessary for the specific job)

· therm.inp (if necessary for the specific job)

· freezing.inp (if necessary for the specific job)

· therm.freeze (if necessary for the specific job)

· boundaries.inp (if necessary for the specific job)

· mkd.inp (optional)

· pmold (if available, see below)

Run

$NX/makedir.pl > makedir.log

If TYPE = 3, makedir.pl will create the directories TRAJECTORIES/TRAJi, where i is the number of the accepted condition.

11.5 About energy restrictions

If energy restrictions are applied (SCREEN = 2) or spectrum generation is chosen (TYPE = 1), vertical excitation energy (DE) and the oscillator strength (f) are collected from each initial condition in final_output.nis.isf or final_output files (where isf is each value in nfs array). For each initial condition within E_CENTER ± E_VAR, the quantity f, or proportional to the oscillator strength or Einstein’s coefficients is computed according to prob_kind keyword. The normalization may be done by the maximum A or B values within the window restriction (NORM = local) or within the complete set of values (NORM = global). The relative "transition probability" P = f/f_max, A/A_max or P = B/B_max is compared to a random number in order to select whether the initial condition will be accepted or not.

12 Dynamics inputs

12.1 What is necessary to run

For the input you need the following files:

Always:

control.dyn - with parameters to control de dynamics.

geom - with initial geometry.

veloc - with initial velocity.

Depending on the settings in control.dyn:

jiri.inp - with parameters to control non-adiabatic dynamics.

sh.inp - with parameters to control non-adiabatic dynamics.

<third-party>.par - with third-party program options.

therm.inp - with thermostat options.

freezing.inp - with frozen atoms information.

therm.freeze - with atoms not affected by thermostat.

boundaries.inp - with boundary condition options.

NX_analysis - with instructions for customized analysis of the results.

You also need the input files for the program that will calculate the energies, gradients and non-adiabatic couplings.

JOB_AD - directory containing input files for adiabatic calculation. Energy and gradient for one state.

JOB_NAD - directory containing input files for non-adiabatic calculations. Energy, gradient and non-adiabatic coupling.

For adiabatic dynamics or non-adiabatic dynamics using time-derivative couplings¹⁴, only JOB_AD is necessary.

For mixed, non- and adiabatic dynamics (using the THRES keyword), both JOB_AD and JOB_NAD are necessary.

For non-adiabatic dynamics using non-adiabatic coupling vectors, only JOB_NAD is necessary.

12.1.1 control.dyn

File control.dyn should contain the main parameters for the dynamics:

File: control.dyn

Namelist: input

Parameter	Default	Description
Nat	= [3]	Number of atoms.
nstat	= [2]	Number of states (dynamics will be performed on the highest one). nstat ≥ nstatdyn.
nstatdyn	= [nstat]	Initial state.
ndamp	= [0]	0 : normal dynamics. 1 : velocity is dumped to zero at each time step.
kt	= [1]	Print output at each kt steps.
dt	= [0.5]	Time step (fs).
t	= [0.0]	Initial time (fs).
tmax	= [10.0]	Maximum time (fs).
nintc	= [3*Nat-6]	Number of internal coordinates to read from the ab intio program. Change default if the system is linear or has redundant coordinates.
mem	= [50]	Core memory for the ab initio program (Mwords). (Relevant if PROG = 1)
killstat	= [1]	Finish dynamics if after a hopping the system remains more than timekill fs on state killstat.
timekill	= [0]	See killstat. 0 : deactivate killstat and timekill
prog	= [1]	Program to compute energies, gradients and non-adiabatic coupling vectors (if available): 0 : Analytical model (see section "Using analytical models") 1 : Columbus 2.0 : Turbomole RI-CC2 / ADC(2) 2.1 : Turbomole TD-DFT 5 : Dftb 6 : Gaussian 7 : Tinker 8 : Dftb+ 20 : Hybrid jobs
thres	= [0 for PROG = 2.x and 5] = [100 for PROG = 0, 1 and 6]	Energy difference threshold to initiate non-adiabatic dynamics (eV).
lvprt	= [1]	Amount of output to print and output files to keep: 0 : Minimal level 1 : Normal level 2 : Debug level (huge amount of data)
Etot_jump	= [0.5]	Kill trajectory if total energy deviate more than Etot_jump eV in one time step.
Etot_drift	= [0.5]	Kill trajectory if total energy deviate more than Etot_drift (eV) in comparison to the value in t = 0.
NXRESTART	= [0]	Restart options: 0: new job. 1: restart job using the content of INFO_RESTART directory. See details in Section 12.4.

[default values] written in inp.f90.

Example 1:

&input

Nat = 6 ! Number of atoms

nstat = 3 ! Number of states

nstatdyn = 3 ! Initial state

dt = 0.5 ! Time step (fs)

tmax = 200.0 ! Total time (fs)

prog = 2.1 ! Dynamics with Turbomole (TD-DFT).

/&end ! Do not forget &input and /&end.

Example 2: Ethylene S₁-state dynamics, calculating S₀, S₁ and S₂ energies. Time step of 0.1 fs during 200 fs and output printed at each 2 fs. 100,000,000 words of memory. If potential energy difference between two consecutive states drops below 2.0 eV, start non-adiabatic dynamics. Finish dynamics if after a hopping to S₀, the system remains more than 10 fs on this state. Get energies, gradients and non-adiabatic couplings with Columbus.

&input

Nat = 6

nstat = 3

nstatdyn = 2

kt = 20

dt = 0.10

tmax = 200.0

mem = 100

thres = 2.0

killstat = 1

timekill = 10.0

/&end

Do not forget "&input" and "/&end".

12.1.2 Geometry

The geometry input is (free format, au):

Name: geom

Symbol_1 Z_1 x_1 y_1 z_1 M_1

Symbol_2 Z_2 x_2 y_2 z_2 M_2

Symbol_nat Z_nat x_nat y_nat z_nat M_nat

where Z is the nuclear charge, x,y,z are the Cartesian coordinates, and M, the atomic masses.

Example:

C 6. -1.27572383 0.00000000 0.00000000 12.00000000

C 6. 1.27572383 0.00000000 0.00000000 12.00000000

H 1. -2.34867651 0.00000000 -1.75798067 1.00782504

H 1. -2.34867651 0.00000000 1.75798067 1.00782504

H 1. 2.34867651 0.00000000 -1.75798067 1.00782504

H 1. 2.34867651 0.00000000 1.75798067 1.00782504

12.1.3 Velocity

The velocity input is (free format, au):

Name: veloc

vx_1 vy_1 vz_1

vx_2 vy_2 vz_2

vx_nat vy_nat vz_nat

Example:

3.226196727134333E-004 -7.803939823701649E-004 3.501212063660452E-004

8.831202616257374E-005 1.103339279899924E-003 -7.468292758584672E-004

-1.994896289256706E-004 4.345679802152278E-004 -6.123248174920957E-004

-3.735660908785368E-003 2.225120145326573E-003 1.414904777249870E-003

-2.272632671568894E-003 -2.177375422081543E-003 3.098061686476041E-003

-6.469986389583130E-004 1.402416374342127E-004 2.362055042392439E-003

12.1.4 Freezing atoms

Cartesian coordinates of specific atoms can be kept frozen along the dynamics. For that, a list of atoms should be given in a file called freezing.inp. This file should be given together with the other input files. The atoms in the list are identified by its positions in the geom file and the list is blank separated.

Example of freezing.inp:

1 3

(Cartesian coordinate of atoms one and three in geom file will not change during the dynamics.)

If a file freezing.inp is present in the input, the atoms defined there will also not be affected by the thermostat.

The velocity of the atoms to be kept frozen should be set to zero in the veloc file (section 12.1.3). If the velocity of the frozen atoms is not initially set to zero, the program will stop with an error message.

Be aware that this algorithm may introduce spurious rotations and translation in the molecule depending on how the initial velocities of the remaining atoms were generated.

12.1.5 Specific input for quantum-chemistry electronic-structure calculations

Please, refer to the documentation of each particular program to see details on their inputs. Here, we present a brief summary of the main option that should be selected for some frequent jobs.

12.1.5.1 Using analytical models

Analytical models for potential energies, gradients and non-adiabatic coupling vectors can be used in the molecular dynamics.

The analytical model is supposed to be a program in any language that reads molecular geometry from geom and velocities (if needed) from veloc, and writes the potential energies to epot, gradients to grad, and non-adiabatic coupling vectors to nad_vectors. The format of each one of these files is described in the section 16.3. If the analytical model requires additional files besides the executable file (for example, parameter inputs) they must be put in JOB_NAD directory. During the dynamics, the content of this directory is copied to the same location as the other Newton-X input files.

The specific location of the analytical model can be set via analyt.par. The options are given below.

Name: analyt.par

Parameter

Default

Description

anmod

[analytical.model]

File name of the executable containing the analytical model. The default is the two-dimension conical-intersection analytical model proposed by Ferretti et al. ¹³.

path

[$NX]

Absolute path to anmod file.

Example of analyt.par file:

anmod = my_model.x

path = /home/model/

The parameters of the two-dimension conical intersection model (default option) can be set via file con_int.dat in JOB_NAD directory. The options are given below.

Name: JOB_NAD/con_int.dat

namelist DAT

Parameter	Default	Description
a	= [3.0]	See Ref. ¹³.
b	= [1.5]
K_x	= [0.02]
K_y	= [0.10]
D	= [0.01]
X₁	= [4.0]
X₂	= [3.0]
X₃	= [3.0]
g	= [0.04]

Example of con_int.dat

&DAT

alpha=3.0

beta=1.5

kx=0.02

ky=0.10,

delta=0.01

x1=4.0

x2=3.0

x3=3.0

gamma=0.04

&END

The atomic masses must be set directly in the geom file.

12.1.5.2 Columbus

Columbus 5.9

Adiabatic dynamics: JOB_AD directory

Prepare a set of input files for a Columbus job (geometry optimization, one iteration, NROOT option). For MCSCF jobs, prepare input for CI gradient, but set "Maximum excitation level" to 0 in the CIDRT input. It is also possible to use MCSCF gradients, but in this case only single state dynamics are allowed (NSTAT = 1).

Non-adiabatic dynamics using non-adiabatic coupling vectors: JOB_NAD directory

Prepare input files for a single-point non-adiabatic coupling calculation.

Non-Adiabatic dynamics using time-derivative couplings: JOB_AD directory

Columbus 7.0

Adiabatic and non-adiabatic (using coupling vectors or time-derivative couplings) dynamics are directly available at the SA-MCSCF and MR-CI levels. All inputs can be performed using the “non-adiabatic coupling” option in Columbus and only the necessary terms are computed for each case. Optionally, adiabatic and time-derivative MR-CI dynamics may still be performed through “geometry optimization”. Please set version=7.0 when using Columbus 7.0.

Parameters to control Columbus jobs

Name: columbus.par

Parameter	Default	Description
VERSION	= [5.9]	Columbus version
mem	= [200]	Core memory for Columbus (Mwords). The same as mem in control.dyn, but with priority over it. (1 GB = 134 Mwords).
cirestart	= [0]	0: do not use previous CI vector 1: use previous CI vector
mocoef	= [4]	0: use the same mocoef file in all time steps 1: use the mocoef file from the previous time step k: Lagrangean extrapolation of molecular orbitals at order k-1 (k ≤ 11) (see Section 12.1.10)
prt_mo	= [20]	Save mocoef file to DEBUG file every prt_mo timesteps. (Only if mocoef = 1.)
ivmode	= [8]	Initial CI-vector generation mode. See Columbus docs (ciudg program) for a list of options. Cirestart keyword has priority over ivmode keyword.
CITOL	=[“1E-4”]	Tolerance for MR-CI calculations
reduce_tol	= [1]	0: keep rtolci and rtolbk in ciudgin as CITOL for all states. 1: use CITOL only to nstatdyn. For all other states use 1E-3. 2: use CITOL only to nstatdyn. Use 1E-3 for states coupled to nstatdyn according to transmomin Columbus file. For all other states use 1E-1. Note that energies computed with 1E-1 are not reliable estimates. Although they are written to the output files, they do not have significance.
mc_conv	= [0]	If mcscf calculation do not converge: 0: warn and continue 1: kill trajectory
ci_conv	= [0]	If CI calculation do not converge: 0: warn and continue 1: kill trajectory
quad_conv	= [60]	Set value of NCOUPL in mcscfin
MULL_pop	= [0]	1: print out the Mulliken populations from the Columbus calculation to nx.log

12.1.5.3 Turbomole

Adiabatic dynamics: JOB_AD directory

Prepare a set of input files for a Turbomole job at TD-DFT or RI-CC2 level without symmetry and copy them to JOB_AD directory.

The initial geometry (Turbomole coord file), the number of states (Turbomole control file), and the state for which the gradient should be computed (Turbomole control file) are automatically set by Newton-X during the program execution, overwriting the options in JOB_AD.

For TD-DFT dynamics, Turbomole control file is automatically changed to have:

$soes all NSTAT-1

$exopt NSTATDYN-1

For RI-CC2 or ADC(2) dynamics, Turbomole control file is automatically changed to have:

$excitations

irrep=a nexc=NSTAT-1

exgrad states=(a NSTATDYN-1)

If Turbomole auxbasis file is provided in JOB_AD, Newton-X assumes that the resolution-of-identity (RI) should be used. For CC2 or ADC(2) dynamics, auxbasis must be always provided and the ricc2 program is invoked every time step. For DFT dynamics, the presence of auxbasis file is optional. If it is provided, ridft program is invoked, otherwise dscf program is invoked.

The Newton-X input for CC2 and ADC(2) dynamics are exactly the same. The only difference is in the Turbomole control file within JOB_AD directory, which should be adequate to the desired method. During the execution, Newton-X search for ADC(2) keyword in the Turbomole control file. If it is not found, the CC2 ground state energy is used. If it is found, MP2 ground state energy is used.

If you want to perform only adiabatic dynamics, set THRES = 0 (default in control.dyn file) to avoid that Newton-X starts the non-adiabatic dynamics calculations.

The Turbomole mos file is not updated during the dynamics.

Non-Adiabatic dynamics using TD-DFT time-derivative couplings: JOB_AD directory

Prepare the Turbomole input files as explained in the previous item. The non-adiabatic dynamics options are described in section 12.1.7.

Set THRES = 100 in control.dyn file.

Parallel Turbomole (SMP)

If you have the SMP parallelized version of Turbomole executables available (dscf_smp, grad_smp, egrad_smp and ricc2_smp) you can use them by specifying

parallel = <ncores>

in the file turbomole.par with <ncores> being the number of cores to use. The environment variable $OMP_NUM_THREADS is set automatically according to this number and needs not to be set by the user.

Parameters to control Turbomole jobs

Name: turbomole.par

Parameter

Default

Description

PARALLEL

= [1]

Number of cores to use for parallel Turbomole (smp, no mpi!)

Inputs for old versions of Newton-X

For Newton-X versions prior the 1.0.8, additional input files should be provided for TD-DFT simulations if NSTATDYN = 1 and nstat > nstatdyn. In this case, the following files are required in JOB_AD directory:

control.opt - Turbomole control file for ground state gradient (GRAD) calculation.

control.sp - Turbomole control file for excited state single point (ESCF) calculation.

For excited state gradient calculations (nstatdyn > 1), no additional files are needed.

12.1.5.4 Dftb

Adiabatic dynamics in the ground and excited states can be performed using the time-dependent density functional theory tight binding (TD-DFTB) method. Information about the program can be obtained at www.dftb.org.

Input for ground state dynamics: JOB_AD directory

Prepare Dftb input for conjugated-gradient ground-state calculation (code 4 in the Dftb program). The Dftb input and geometry information must be called dftb.in and in.gen, respectively.

Input for excited state dynamics: JOB_AD directory

Prepare TD-Dftb input for conjugated-gradient excited-state calculation (code 4 in the Dftb program). The Dftb input and geometry information must be called dftb.in and in.gen, respectively.

Parameters to control the Dftb jobs

Name: dftb.par

Parameter	Default	Description
dftb_exec	[dftb]	Name of the Dftb executable file. $DFTB variable must be defined in the system. Newton-X will run $DFTB/<dftb_exec>.
other_state	[1]	It is possible to perform dynamics on one surface, and at the same time to keep track of the properties of other states (energies and oscillator strengths). In this case, however, Dftb must run twice per time step, which makes the calculation more expensive. 0 - Do not monitor other states 1 - Monitor other states
mult	[S]	Multiplicity of the states. Only states with the same multiplicity are allowed. S - singlet T - triplet

12.1.5.5 Gaussian 03

Surface-hopping non-adiabatic dynamics between the ground and the first excited state can be performed at CASSCF level with Gaussian 03.

Content of JOB_NAD directory

Prepare input files for a single point conical intersection calculation at CASSCF level. The input file must be called gaussian.com and the geometry must be given in Cartesian coordinates. The check point file named gaussian.chk containing the initial molecular orbitals must be provided as well. A suitable example of gaussian.com content is:

%chk=gaussian

%mem=20000000

#P OPT=(Conical,MaxCycle=1) CAS(4,3) IOp(5/7=200) Guess=read 3-21G Nosymm UNIT=AU

methaniminium

1 1

N 0.00000059 0.00000062 1.20440093

C 0.00000002 0.00000002 -1.32963740

H -1.22595386 1.08158260 2.22587621

H 1.22595117 -1.08158568 2.22587634

H 1.16944597 1.32484957 -2.41516400

H -1.16944649 -1.32484997 -2.41516296

Note that OPT=(Conical,MaxCycle=1), Nosymm, Unit=AU, and Guess=read keywords and options should be given exactly as in this example. IOp(5/7=MaxIt) controls the maximum number of iterations in the CASSCF calculation.

Parameters to control Gaussian jobs

Name: g03.par

Parameter

Default

Description

mocoef

= [1]

0: use the initial molecular orbitals in all time steps

1: use the check point file from the previous time step as the source of molecular orbitals

prt_mo

= [20]

Save gaussian.chk file to DEBUG file every prt_mo timesteps. (Only if mocoef = 1.)

Newton-X executes Gaussian 03 by invoking the command:

. $g03root/g03/bsd/g03.profile;$g03root/g03/g03 gaussian.com

If this path is not adequate for your system, you can change it in $NX/run-g03.pl program.

12.1.5.6 Tinker

Adiabatic ground-state dynamics with forcefield-gradients and –energies can be performed with Tinker.

Content of the JOB_AD-directory

An Tinker-XYZ file called tinkin.xyz and the corresponding keyword-file tinkin.key have to be present. The file tinkin.key has to state at least the name of the parameter file ('parameters <filename>') and 'digits 12'.

The command executed by Newton-X is

testgrad tinkin.xyz –k tinkin.key y n n

If the Tinker-executables are not in the $PATH, testgrad.x can also be provided as binary in the JOB_AD directory.

12.1.5.7 Dftb+

Adiabatic dynamics in the ground can be performed using the density functional theory tight binding (DFTB) method.

Input for ground state dynamics: JOB_AD directory

Prepare Dftb+ input for conjugated-gradient ground-state calculation (code 4 in the Dftb program). The Dftb+ input and geometry information must be called dftb.in and in.gen, respectively.

Parameters to control the Dftb+ jobs

Name: dftb+.par

Parameter

Default

Description

dftbP_exec

[dftb+]

Name of the Dftb+ executable file. $DFTBPLUS variable must be defined in the system. Newton-X will run $DFTBPLUS/<dftbp_exec>.

12.1.5.8 Hybrid Gradients (QM/MM)

The hybrid gradients module allows adiabatic and non-adiabatic dynamics with combinations of programs (Columbus, Tinker, Turbomole and Analytical Model). The current implementation is described in Ref.³⁰.

General explanations of hybrid calculations

Energies and gradients for subsets of atoms are treated with different programs and the partial gradients are then joint into a resulting total (hence 'hybrid') energy and gradient. For that purpose the set of atoms of the whole system is split in disjoint regions. These regions need not follow physical reasoning (they often will, but the can also pick e.g. single atoms out of molecules), but are logical entities for the definition of the single partial calculations.

The whole calculation is split into jobs. Each of these jobs can treat one or more regions of atoms and the partial result are multiplied by a user-defined factor before being added to the total. One region can be treated by multiple jobs and to care about 'double counting' of atoms is left to the user completely.

For Columbus and Turbomole there is the possibility to include regions only as point charges and not as atoms with basis-set. A hybrid setup (which is done in the JOB_AD or JOB_NAD directory) consists of some general information about the hybrid setup ('control' parameters), the definition of the atoms, some of their properties and membership to the regions and the definition of the partial jobs, which regions they are concerned with and how they shall be put together to the overall result. The main input file is named hybrid_gradients.inp.

The nonadiabatic couplings, oscillator strengths and other nonadditive properties are given by only one job. It is possible to restrict the hybrid NAD-vectors to some regions. The NAD-vector components on all atoms not belonging to these regions are replaced by Zeros in this case. Also for back hoppings only the kinetic energy of these nad-regions is regarded as available energy.

For the treatment of bonded interaction between two regions link atoms can be inserted in bond connecting them. The gradient- and nonadiabatic coupling vector elements of these link atoms will be distributed to the two atoms between which it has been inserted. To avoid over-polarization effects, the point charges near the link atom can be set to zero and, so as to retain the overall charge, re-distributed to a set of other atoms. The job, where the link atom is inserted is (suggestively) called ‘QM_JOB’ regardless of the method really used in this job. Similarly is the naming of ‘QM_ATOM’ and ‘MM_ATOM’.

Format of the hybrid_gradients.inp file

Each section begins with $<name> and ends with $end, where <name> can be 'control', 'job' or 'atoms'. Parameters for sections '$control' and '$job' are set in key=value pairs. All key=value pairs have to be separated with blanks.

Name: hybrid_gradient.inp

Section	Keyword	Description
$control	natoms	number of atoms (optional)
	properties	ID of the job, that gives the properties (oscillator strength, nonadiabatic couplings, optional)
	NADREGIONS	array of regions to which the hybrid nad-vector shall be restricted
$job	id	unique numeric (integer) identification for the job
	program	name of the program to be used for this job (columbus, tinker, turbomole)
	regions	array of regions, that shall be treated by this job (comma seperated)
	factor	pre-factor with that the result of this job enters the total result (energy, gradients)
	pointcharges	array of regions, that shall be treated as point-charges (for columbus and turbomole, comma seperated, optional)
$atoms	(For each atom give in one single line and in this order:)
	<name>	Label for an atom (compulsory). Can be anything without blanks, comma and similar. The type of atom is recognized from the nuclear charge.
	<nuclear charge>	Nuclear charge.
	<x,y,z coords>	Cartesian coordinate (Bohr).
	<pointcharge>	Effective point charge for this atom. If this atom is not treated as point charge this can be left to 0.0000.
	<mass>	Atomic mass for this atom (amu).
	<region>	The region in which this atom is.

link	QM_JOB	ID of the job where the link atom shall be inserted.
	QM_ATOM	number (=position in input) of the atom connected to the link atom at the QM-side.
	MM_ATOM	number (=position in input) of the atom connected to the link atom at the MM side.
	RATIO	ratio of distances QMat→LINKat/QMat→MMat.
	ZERO	comma seperated list of atoms for which the point charge shall be set to 0.000 in the QM_JOB.
	SCATTER	comma seperated list of atoms to which the accumulated charge of the ZERO-atoms shall be distributed.

The section $atoms is an extension of the Newton-X geom. with the two additional columns, <pointcharge> and <region>.

Each atom has to be in exactly one region.

All values in section $atoms have to be separated with blanks.

Every section has to be ended with $end. No section (except '$jobs') is allowed to appear twice.

There are no further formatting rules to the file.

Example of hybrid_gradient.inp file:

$job ID = 1 regions = 1,2 program = columbus pointcharges = 2 factor = 1 $end

$job ID = 2 regions = 1,2 program = tinker factor = 1 $end

$job ID = 3 regions = 1 program = tinker factor = -1 $end

$control properties = 1 nadregions = 1 natoms = 9 $end

$atoms

C 6.0 -0.58698100 -0.10086826 0.12744020 12.00000000 0.0000 1

O 8.0 1.68492145 0.00195156 0.55748655 15.99491464 0.0000 1

N 7.0 -1.73272850 -2.27661770 -0.35180374 14.00307401 0.0000 1

H 1.0 -1.76941972 1.58570542 0.11857328 1.00782504 0.0000 1

H 1.0 -3.64014444 -2.31932440 -0.70744017 1.00782504 0.0000 1

H 1.0 -0.69532029 -3.86742067 -0.35771869 1.00782504 0.0000 1

OW 8.0 -6.90020375 -2.40598424 -1.31691242 15.99491464 -0.8340 2

HW 1.0 -7.75621060 -2.20340156 -2.90625351 1.00782504 0.4170 2

HW 1.0 -8.26652850 -2.54571452 -0.12804111 1.00782504 0.4170 2

$end

This example is explained in details in the Newton-X tutorial.

Example including link atoms (some lines not displayed):

# 2-butene with cyclohexane rings attached to the sides

# Z2-scheme is used for point charges,i.e. first and second

# neighbour to the link atom have zero charges and the

# their charge is distributed to their next-neighbours on the MM side

$job ID=1 regions=1,2,3 program=turbomole turbo_method=td-dft pointcharges=2,3 factor=1 $end

$job ID = 2 regions = 1,2,3 program = tinker factor = 1 $end

$job ID = 3 regions = 1 program = tinker factor = -1 $end

$link qm_job=1 qm_atom=1 mm_atom=11 ratio=0.713 zero=11,12,13,14 scatter=15,16,17 $end

$link qm_job=1 qm_atom=4 mm_atom=34 ratio=0.713 zero=34,35,36,37 scatter=38,39,40 $end

$control properties = 1 nadregions=1 natoms = 56 $end

$atoms

C 6.0 1.44087838 0.85091533 -3.47469613 12.00000000 0.0000 1 #atom 1

# linked to atom 11

C 6.0 1.23567680 -0.48733580 -6.00452731 12.00000000 0.0000 1

C 6.0 2.80448775 -2.25355699 -6.92713507 12.00000000 0.0000 1

C 6.0 5.12903987 -3.34336584 -5.64481794 12.00000000 0.0000 1 #atom 4

# linked to atom 34

H 1.0 1.17513186 2.91664300 -3.78035178 1.00782504 0.0000 1

H 1.0 3.36724819 0.62483228 -2.66873550 1.00782504 0.0000 1

H 1.0 -0.40125500 0.07262218 -7.17351567 1.00782504 0.0000 1

H 1.0 2.39776578 -3.02448967 -8.82452971 1.00782504 0.0000 1

H 1.0 4.98819669 -5.44423296 -5.69313068 1.00782504 0.0000 1

H 1.0 5.18046121 -2.82088113 -3.61093594 1.00782504 0.0000 1

CT 6.0 -0.52739233 -0.13060275 -1.57905516 12.00000000 -0.1200 2 #atom 11

HC 1.0 -0.25927609 -2.20323737 -1.33315643 1.00782504 0.0600 2

HC 1.0 -2.45999068 0.14865342 -2.36823880 1.00782504 0.0600 2

CT 6.0 -0.33457412 1.20138205 0.99721226 12.00000000 -0.1200 2

HC 1.0 1.56440222 0.83732631 1.83045298 1.00782504 0.0600 2

HC 1.0 -0.45614398 3.28079700 0.68304340 1.00782504 0.0600 2

: : : : : : : : : : : :# lines omitted

CT 6.0 -2.23690095 2.08870485 5.28849685 12.00000000 -0.1200 2

HC 1.0 -0.32070164 1.91506358 6.14394752 1.00782504 0.0600 2

HC 1.0 -2.48154867 4.11306830 4.76133507 1.00782504 0.0600 2

CT 6.0 7.60087913 -2.51304285 -6.92702736 12.00000000 -0.1200 3 #atom34

HC 1.0 9.21362621 -3.51540650 -6.01449183 1.00782504 0.0600 3

HC 1.0 7.56925645 -3.13134234 -8.93790950 1.00782504 0.0600 3

CT 6.0 8.05428868 0.35505686 -6.77227957 12.00000000 -0.1200 3

HC 1.0 7.92399962 0.94155983 -4.75349648 1.00782504 0.0600 3

HC 1.0 6.51305828 1.36826754 -7.78612331 1.00782504 0.0600 3

CT 6.0 10.63311478 1.21580066 -7.83079266 12.00000000 -0.1200 3

HC 1.0 12.15373961 0.13167423 -6.84744154 1.00782504 0.0600 3

: : : : : : : : : : : :# lines omitted

$end

Generating and updating the remaining input files

After having written hybrid_gradient.inp file, the complete set of input files can be created by run

$NX/hybrid_read_onefile.pl

Subdirectories for the different jobs will be created as well. Appropriate geom files for each job are provided in these subdirectories. Set up the third-party single jobs inputs in the subdirectories in the same way as for normal Newton‑X dynamics. If you chose to include some atoms as point-charges to a Columbus- or Turbomole-job you have to set the jobs up appropriately. Refer to the documentations of the programs for information on how to do that (hybrid_read_onefile.pl will provide some useful files containing most of things, that have to be done extra to a normal setup).

In the subdirectories for partial hybrid jobs to be computed with Columbus including point-charges two files 'potential.xyz' and 'elpotin' are provided. Do not delete these files!

In the subdirectories for partial hybrid jobs to be computed with Turbomole including point-charges two files ‘pointcharges’ and ‘control-additions’ are provided. Do not delete the pointcharges-file. The contents of the file ‘control-additions’ have to be inserted in the control-file at an appropriate position after setting up the Turbomole calculation.

If any change is done to hybrid_gradient.inp, hybrid_read_onefile.pl should be run again to update the remaining files.

12.1.6 Thermostat control

Thermal-equilibration may be obtained during the dynamics by using a thermostat ²⁹. In the current Newton-X version, the Andersen thermostat is available ²⁸. The Andersen-Lowe thermostat ⁴² is being implemented and it will available soon. The parameters of the thermostat may be adjusted by using nxinp tool, in the input section “Set General Options”. The options are given below.

Name: therm.inp

namelist therm

Parameter	Default	Description
ktherm	= [1]	Thermostat type: 0 - No thermostat. 1 - Andersen ²⁸. 2 - Andersen-Lowe ⁴² (in development).
temp	= [300]	Temperature (K).
kts	= [1]	Turn the thermostat on at time step kts. (kts ³ 1, integer)
lts	= [-1]	Turn the thermostat off at time step lts. (lts > kts, integer) -1 - Do not turn it off.
nstherm	= [1]	If the system is in the excited state: 0 - turn the thermostat off. 1 - apply the thermostat.
gamma	= [0.2]	Collision frequency (fs^-1).
radius	= [10.0]	Collision radius (bohr). (Only in Andersen-Lowe thermostat.)
iseed	= [1]	Random number seed for the thermostat. 0 - default seed value. 1 - generate random seed. > 1 - use this value (integer) as the seed.
lvp	= [1]	Print level. (In any case, relevant information is written to Newton-X log.) 1 - Do not print thermostat log file. 2 - Print minimum thermostat log file. 3 - Print debug-level thermostat log file.

Example of therm.inp file:

&therm

ktherm = 1

kts = 1

lts = -1

nstherm = 1

temp = 300

gamma = 0.6

iseed = 0

lvp = 3

&end

You can choose a subset of atoms that will not be affected by the thermostat. For that, a list of atoms should be given in a file called therm.freeze. This file should be given together with the other input files. The atoms in the list are identified by its positions in the geom file and the list is blank separated.

Example of therm.freeze:

1 3

(Atoms one and three in geom file will not be affected by the thermostat.)

If a file freezing.inp is present in the input, the atoms defined there will also not be affected by the thermostat (see section 12.1.4).

12.1.7 Non-adiabatic dynamics control

The non-adiabatic dynamics is controlled by two input files, sh.inp and jiri.inp.

SH is a stand-alone module for surface-hopping. It reads nuclear velocities and non-adiabatic couplings and gives as output the electronic wavefunction expansion coefficients on the adiabatic set of states. Switch from one adiabatic PES to another is ruled by the Tully's fewest switches algorithm ¹⁰.

The parameters for non-adiabatic dynamics can be set via nxinp tool, at the input option “Set non-adiabatic dynamics”. The options are given below.

Name: sh.inp

namelist shinp

Parameter	Default	Description
vDOTH	= [0]	Dynamics with non-adiabatic coupling vectors or with time-derivative couplings (HST model ¹⁴^,¹⁵). -1 – Local-diabatization method.¹⁶ 0 - Compute non-adiabatic coupling vectors. 1 - Compute time-derivative couplings (Columbus and Turbomole TD-DFT).¹⁴ The default is reset to 1 if PROG = 2.1.
seed	= [1]	0 a default random number seed is used. 1 a randomized seed is used >1 random number seed
integrator	= [5]	selects the integrator of the TDSE: 0 - a "home-made" integrator, see Ref. ¹³. 1 - standard 4th order Runge-Kutta. 2 - Adams Moulton predictor-corrector, 5th order. 3 - Adams Moulton predictor-corrector, 6th order. 4 - Unitary propagator. 5 - Butcher, 5th order.⁴³ 6 – Unitary propagator for Local-Diabatization method. Good choices are usually 3, 4, and 5. For vdoth = -1 (local diabatization) integrator must be 6.
phase	= [1]	Integrate the phase along the trajectory (only for debug purposes): 0 - No. Phase is always zero. 1 - Yes.
nohop	= [0]	Force the hopping at the certain time step. 0 - The trans. prob. are computed and hopping is allowed (normal surface hopping). -1 - Hopping is not allowed at any time. n - Hopping is forced at (and only at) timestep n (n = positive integer).
forcesurf	= [1]	Force the hopping to the surface forcesurf (ground state = 1).
nrelax