Job types
MonaLisa offers few job types:
- single point calculation
- structure optimization
- nudge elastic band calculations
- thermodynamic calculations
- molecular dynamics
Job type is specified in the run datagroup which is define in the MonaLisa input file as follows:
run{
<jobtype>
}
Single point calculation
It is the simplest calculation type which returns the energy of investigated system. This type of calculation is invoked by:
run{
single
}
Structure 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
}
To modify the options of optimizer the optimizer datagroup has to be defined:
optimizer{
keyword1
keyword2
...
}
To specify the options of optimizer the following keywords may be used:
keywords:
max_cycles = <integer>// Maximal number of optimization steps (default: "200")
en_change = <double>// The energy convergence criteria (default: "2.7e-4" in eV)
max_grad = <double>// Maximal gradient component convergence criteria (default: "1.0e-2" in eV/A)
rms_grad = <double>// Root mean square of gradient convergence criteria (default: "5.0e-3" in eV/A)
max_step = <double>// Maximal displacement convergence criteria (default: "5.0e-3" in A)
rms_step = <double>// Root mean square of displacement convergence criteria (default: "1.0e-3" in 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 posses the relevant imaginary frequency in harmonic approximation. The MonaLisa optimizer may follow this transition mode "uphill" the PES, whereas the oder 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 nr. 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 matrixhess_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 expenive is recalculation at each step, i.e. =1
...
}
Partial optimization:
MonaLisa may optimize all coordinates of the system, but it can also optimize only selected atom which is called partial optimization. Specification of the partially optimized structure may be done using two complementary options:
(1) 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
...
}
(2) 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, only the former one will be taken into account
In some leaf interfaces (e.g. vasp and turbomole) the partial optimization option force the programs to calculate the partial hessian
Constrained optimization:
MonaLisa may 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{
keyword1
keyword2$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:
run{
optimize}
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 (firs 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 eV^-1. The here defined collective variable is of the form: CV = 1.0*d(1,2)-1.0*d(2,3).
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.
Thermodynamic calculations
This type of job, shortly called thermo, calculates thermodynamic functions of the system (Gibbs free energy, enthalpy, entropy ect.) on the basis of harmonic oscillator approximation. In the near future also anharmonic diagonal oscillator approach of Piccini and Sauer will be available.
Thermodynamic calculations within harmonic approximation rely on the frequencies of the normal modes of the system, calculated by diagonalize of the mass weighted hessian matrix. Therefore the hessian matrix has to be provided by interfaces.
The program generates few addition output files. The most important is "Tscan.out" which consists the most important thermodynamic functions for different temperatures and pressures.
To start thermo calculations one has to specify:
run{
thermo
}
keywords:
temp_init
= <double>// Initial temperature (default: "298" in K)temp_final
= <double>// Final temperature (default: "298" in K)temp_step
= <double>// Temperature interval of scan between initial and final temperature (default: "10" in K)press_init
= <double>// Initial pressure (default: "0.100" in MPa)press_final
= <double>// Final pressure (default: "0.100" in MPa)press_step
= <double>// Pressure interval of scan between initial and final values (default: "0.1" in MPa)
phase = <string>// Phase of the system (default: "gas", options: "gas" or "solid")
symmetry = <string>// Symmetry of the molecule (default: "c1", options: c1, cs, c2, c2v, c3v, c2h, coov, d2h, d3h, d5h, dooh, td, oh)
linear= <bool> // Is the molecule linear? (default: false, options: true or false)
sub_molecules = <string>// 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 = <integer>// Use N first frequencies to perform calculations (default: "0", i.e. it does not use this option)
freq_cutoff = <double>// 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 = <bool>// Calculate hessian numerically using central difference method (two displacements per degree of freedom). (default: "false", i.e. it does not use this option). By default a differentiation step of 0.015 A is used.
numerical_displ = <double>// Differentiation step for numerical hessian calculations (default: "0.015" n A).
Example:
Here you can find a useful example of the thermal analysis performed on the basis of precalculated VASP hessian placed in the OUTCAR file one level above the current directory:
# Interface which provides the hessian matrix (here VASP)
toplevel{
interface = vasp{
program = aprun -B /sw/chem/vasp/5.3.5/mpp1/INTEL/vasp_gamma_cd
read_hessian = ../../OUTCAR # Reads the hessian from the OUTCAR file placed in the folder above (two folders up from the toplevel/)
read_energy = ../../OUTCAR # Reads the energy from the OUTCAR file placed in the folder above (two folders up from the toplevel/)
}
}
# Structure of the molecule
structure POSCAR
# Call the thermo job
run{
thermo
}
# Thermo calculation options
thermo{
temp_init = 323
temp_final = 1023
temp_step = 10
press_init = 0.1
press_final = 0.1
symmetry = c1
linear = false
phase = solid
}
Don't forget to run setup before thermo calculations