Runners

Runners are the drivers which handle the results of the interface calculations.
A basic and most frequently used runner is optimizer which reads the energy, gradients (and optionally Hessian) from the external molecular modeling program and drives the structure optimization.

Types

Single point

It is the simplest type of calculations, which returns only the energy of initial molecular structure.
This runner is invoked by:

run{
    single
}

Keywords

No keywords available yet.

Energy and Gradient

It is and extension of the Single point calculation. This runner returns energy and gradient vector of the initial molecular structure. The coordinates are not affected by this job type.
This runner is invoked by:

run{
    singlegrad
}

Keywords

No keywords available yet.

Optimization

It is the most frequently used job type. The program minimize the energy of the system by changing the coordinates of the molecule. There are few modes of optimization but the default one is the optimization in the normal mode coordinates using BFGS methodology (Diagonal matrix as an initial guess for the Hessian matrix and BFGS formula for the Hessian update). It is invoked by:

run{
    optimize
}

Keywords

The keywords of optimizer are defined within the optimizer{ } datagroup:

optimizer{
    key1 = val1
    key2 = val2
    ...
}

The valid keywords has to be defined within this datagroup. The optional keywords are:

max_cycles = <int> - Maximal number of optimization steps (default = 200)
en_change = <double> - The energy convergence criteria (default = 2.7E-4 [eV]
max_grad = <double> - Maximal gradient component per degree of freedom (default = 1.0E-2 [eV/A])

Transition Structure Optimization

MonaLisa may optimize transition structures using quasi-Newton approach starting from reasonable guess of the transition state. Reasonable guess means the structure which possess the relevant imaginary frequency in harmonic approximation. The MonaLisa optimizer may follow this transition mode "uphill" the PES, whereas the other modes are followed "downhill". In this way the saddle point of the first order may be find. To use this optimization one has to start from calculations of the hessian matrix. The initial hessian matrix may be updated during the optimization using gradients and the Powell formula (by default). One may also start the optimization reading the old hessian matrix. This option enables the "hybrid" approach in which the initial guess of the hessian matrix is taken from the lower level calculation, whereas it is updated using the high level gradients. This approach works, however, not always, but in the case of the methods where the hessian calculations are too expensive it represents an alternative worth to consider. To run this calculation few additional options have to be specified:

optimizer{

   trans_vector =  1 # Vector which will be followed during TS optimization (normally the lowest frequency, i.e. nr. 1)

   keepTmode = false # Default true - keep always the mode number "trans_vector" (the vectors are numbered from the lowest to the highest, the lowest is no 1), if false - it keeps the initial mode (which may change the number during optimization)
   init_hess  = calc   # calculate initial hessian, default value is "unit" which uses the unit  matrix, but then it tries to read the monalisa.hessian, so if the option is not specified it will read the old monalisa.hessian (if present) which gives opportunity to use precalculated at the lower level hessian matrix

   hess_calc = true # calculation of the hessian, if true then hessian will be calculated every hess_calc_freq steps (by default at every step)

   hess_calc_freq = 30 # how often hessian matrix will be recalculated, most safe, but also most expensive is recalculation at each step, i.e. =1
  ...
}

Partial Optimization

MonaLisa by default optimize all coordinates of the system, but it may also perform the partial optimization, where only selected atoms are allowed to move. Specification of the partially optimized structure may be done using two complementary options:

  • carfreeze - specification of the atoms which should be frozen (the rest will be optimized). Example:
optimizer{
    $carfreeze = 1,3-5,20  #atoms 1,3,4,5 and 20 will be fixed during the optimization 
    ...
}
  • caropt - specification of the atoms which will be optimized (the rest will be frozen). Example:
optimizer{
    caropt = 1,3-5,20 #atoms 1,3,4,5 and 20 will be optimized and the rest will be frozen during the optimization 
    ...
}

Beware! the caropt option is superior to the carfreeze, so if the both keywords are present in the input file, only the former one will be taken into account

Notice In some leaf interfaces (e.g. vasp and turbomole) the partial optimization option force the programs to calculate the partial hessian

Constrained Optimization

MonaLisa can constrain bond distance, bond angles and dihedral angles. The initial values of these parameters are constrained during the optimization, so if you want to change e.g. selected bond distance to some specified value you should change this distance to desired value in some external builder and after this run constrained optimization. If you want to skip this step you may use the restrain option described later. To run the constrained optimization you should specify:

optimizer{
    ...
    $freeze
        bond   1 2
        angle  1 2 3
        tangle 1 2 3 4
    $end
}

In this example we constrain: - the bond distance between atoms 1 and 2 - the bond angle between bonds 1-2 and 2-3 - the dihedral angle between planes 1-2-3 and 2-3-4

To add another constrains one have to repeat the whole line (e.g. another bond line)

Restrained Optimization

Another option to constrain the bond distance, bond angles and dihedral angles is using harmonic potentials to force the predefined structure parameters. As opposed to constrained optimization described above, there is no need to modify the structure before running the actual calculation. The structure restrain is invoked by adding the harmonic potential with predefined structure parameter and the spring strength. Until now only the bond angle restrain is implemented, but the implementation of other restraints is straightforward and will be done soon. To run the restrained optimization you should specify:

optimizer{
    ...
    springs{
        bond 1 2 1.5 1000.0
        bond 1 3 2.5 500.0
        collective 2  1.0 1 2  -1.0  2  3  2.5  1000.0
    }
}

In this example we restrain:

  • the bond distance between atoms 1 and 2 to 1.5 Angstroms using the spring constant of 1000 eV/A
  • the bond distance between atoms 1 and 3 to 2.5 Angstroms using the spring constant of 500 eV/A
  • the collective variable (CV) which consist of two bond distances (first number). The first bond (between atoms 1 and 2) is taken with the weight 1.0, the second bond (between 2 and 3) with weight -1.0, the CV value is equal 2.5 A, and the force constant is 1000.0 1/eV. The here defined collective variable is of the form:

Beware: One have to be careful with selecting correct spring constant. If it is too small, the bond (or other property) will be not equal to predefined one, if it will be too big, some numerical problems with optimization may occur.

Molecular Dynamics

In this job type the trajectories (position and velocity as a function of time) of atoms in the system are determined by numerically solving Newton's equations of motion, where forces acting on atoms are provided by Interfaces.

This job is invoked by:

run{
    md
}

keywords

All keywords are specified within md datagroup:

md{
    key1 = val1
    key2 = val2
    ...
}

The default values specified after =
ensemble = nve statistical ensemble (options: nve (microcanonical), nvt (canonical))
time_step = 1.0 time step [fs]
steps = 100 number of integration steps
restart = off restart flag (options: on, off)
thermost = off thermostat (options: off, nose, gaussian, andersen)
temp = 298.0 temperature in K
dump_freq = 0 dump each n'th step; default (0 - no dump)
mass = off atomic masses; default:off (options: off, physical, equal, user (user defined))
constraints = off impose constraint (options: off, on (read them from .mdconstr file()
seed = 12345 random number generator seed
ramps = 1 number of panels (sections)
verb = 1 output verbosity (options: 0-10)
verb_therm = 1 output verbosity for thermostat (options: 0-10)
verb_constr = 1 output verbosity for constrained md (options: 0-10)
proj_freq = 1 project out translations and rotations each n'th step
shake_tol = 1.0E-5 SHAKE (and RATTLE)) algorithm deviation tolerance
time_zero = 0.0 time of performing simulation [fs]
nh_q = 0.0 Nose-Hoover thermal bath "mass" (for both chains)
nh_q1 = 5.0 Nose-Hoover thermal bath "mass" (chain 1)
nh_q2 = 5.0 Nose-Hoover thermal bath "mass" (chain 2)
ander_coupl = 0.1 Andersen bath coupling (collision probability)
trans_tol = 1.0E-5 Tolerance for projecting out translations
rot_tol 1.0E-5 Tolerance for projecting out rotations
proj_out_t = true project out translations
proj_out_r = true project out rotations
proj_out_tr = true project out translations and rotations
write_math = false write mathematica(R) file
printcif = true print cif file with trajectory

Nudged Elastic Band

The nudged elastic band (NEB) is a method for finding saddle points and minimum energy paths between known reactants and products. The method works by optimizing a number of intermediate images along the reaction path. Each image finds the lowest energy possible while maintaining equal spacing to neighboring images. This constrained optimization is done by adding spring forces along the band between images and by projecting out the component of the force due to the potential perpendicular to the band.

This runner is invoked by:

run{
    neb
}

Keywords

NEB runner is inherit class of optimizer, therefore all keywords are specified within optimizer datagroup:

optimizer{
    key1 = val1
    ...
}

The default values specified after =
number_configurations = 10 number of images
neb_begin = monalisa_init.car initial structure (image 0)
neb_end = monalisa_final final structure (image number_configurations + 1)
neb_transition_guess = initial guess for transition state
kmax = 30.0 maximal force constant
kmin = 5.0 minimal force constant
print_amount = 2 print level neb_type = neb_ci NEB algorithm (options: neb, neb_ci, neb_ci2)
max_cycles = 200 Maximal number of optimization steps
en_change = 2.7E-4 Energy convergence criteria [eV]
max_grad = 1.0E-2 Maximal gradient component per degree of freedom [eV/A]

Thermodynamic calculations

This type of job, calculates thermodynamic potentials of the system (Gibbs free energy, enthalpy, entropy ect.). Thermodynamic calculations within harmonic approximation rely on the frequencies of the normal modes of the system, calculated by diagonalizing of the mass weighted hessian matrix. Therefore the hessian matrix has to be provided by interfaces.
The program generates few addition output files which are placed in the folder Thermo_data/. The most important of them is "Tscan.out" which consists thermodynamic potentials for different temperatures and pressures. To start thermo calculations one has to specify in the input file

run{
    thermo
}

keywords

All keywords are specified within thermo datagroup:

thermo{
    key1 = val1
    key2 = val2
    ...
}

The default values specified after =
temp_init = 298 Initial temperature (in K)
temp_final = 298 Final temperature (in K)
temp_step = 10 Temperature interval of scan between initial and final temperature (in K)
press_init = 0.100 Initial pressure (in MPa)
press_final = 0.100 Final pressure (in MPa)
press_step = 0.100 Pressure interval of scan between initial and final values (in MPa)
phase = gas Phase of the system (options: "gas" or "solid")
symmetry = c1 Symmetry of the molecule (options: c1, cs, c2, c2v, c3v, c2h, coov, d2h, d3h, d5h, dooh, td, oh)
linear = false Is the molecule linear? (options: true or false)
sub_molecules = String containing atoms consisting a partial hessian (default: ""). One may define few partial hessians by separating them with semicolon ";" e.g. "1,3,5-7; 1-20; 7-45,3,5,6"
nfreq = 0 Use N first frequencies to perform calculations (default: "0", i.e. it does not use this option)
freq_cutoff = 0.0 All frequencies below the absolute value of this threshold will be replaced by its value (for example replace all frequencies below 30 cm-1 by 30 cm-1, also imaginary frequencies higher than this value will be replaced)
numerical_hessian = false Calculate hessian numerically using central difference method (two displacements per degree of freedom). By default a differentiation step of 0.015 A is used.
numerical_displ = 0.015 Differentiation step for numerical hessian calculations (in A)