Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
vibrational_scf_programs [2021/07/28 10:25] rauhuterfortvibrational_scf_programs [2022/08/16 14:15] (current) – external edit 127.0.0.1
Line 1: Line 1:
-====== Vibrational SCF programs ====== 
- 
 ===== The VSCF programs (VSCF) ===== ===== The VSCF programs (VSCF) =====
  
Line 44: Line 42:
   * **''VAM''=//n//** The 0D terms of the vibrational angular momentum terms, i.e. $\frac{1}{2} \sum_{\alpha\beta} \hat{\pi}_\alpha\mu_{\alpha\beta} \hat{\pi}_\beta$, and the Watson correction term are by default (''VAM=2'') included.\\   * **''VAM''=//n//** The 0D terms of the vibrational angular momentum terms, i.e. $\frac{1}{2} \sum_{\alpha\beta} \hat{\pi}_\alpha\mu_{\alpha\beta} \hat{\pi}_\beta$, and the Watson correction term are by default (''VAM=2'') included.\\
 ''VAM=1'' adds the Watson correction term (see Eq. \eqref{eq:1}) as a pseudo-potential like contribution to the fine grid of the potential.\\ ''VAM=1'' adds the Watson correction term (see Eq. \eqref{eq:1}) as a pseudo-potential like contribution to the fine grid of the potential.\\
-''VAM=2'' allows for the calculation of the integrals of the ''VAM'' operator using the approximation that the $\mu$ tensor is given as the inverse of the moment of inertia tensor at equilibrium geometry.\\+''VAM=2'' adds the ''VAM'' operator using the approximation that the $\mu$ tensor is given as the inverse of the moment of inertia tensor at equilibrium geometry to the quasi Fock operator.\\ 
 +''VAM=3'' allows for the calculation of the integrals of the ''VAM'' operator using the approximation that the $\mu$ tensor is given as the inverse of the moment of inertia tensor at equilibrium geometry.\\
 When using ''VAM=4'' the expansion of the effective moment of inertia tensor will be truncated after the 1D terms (rather than the 0D term in case of ''VAM=2''). Note that values higher than 2 are only active for non-linear molecules. ''VAM=5'' truncates the series after the 2D term. In almost all cases ''VAM=2'' is fully sufficient. Vibrational angular momentum terms are accounted for in a perturbational manner and do not affect the wavefunction. When using ''VAM=4'' the expansion of the effective moment of inertia tensor will be truncated after the 1D terms (rather than the 0D term in case of ''VAM=2''). Note that values higher than 2 are only active for non-linear molecules. ''VAM=5'' truncates the series after the 2D term. In almost all cases ''VAM=2'' is fully sufficient. Vibrational angular momentum terms are accounted for in a perturbational manner and do not affect the wavefunction.
  
Line 79: Line 78:
 vscf,pot=poly                            !(5) do a VSCF calculation vscf,pot=poly                            !(5) do a VSCF calculation
 </code> </code>
-==== Rovibrational calculations ==== 
- 
-''ROVIB'',//options// 
- 
-By default, the ''VSCF'' program calculates purely vibrational states only. However, the ''ROVIB'' directive allows for the calculation of rovibrational transitions for molecules with Abelian point groups. This includes also the IR intensities once dipole moment surfaces have been computed. For details see:\\ 
-S. Erfort, M. Tschoepe, G. Rauhut, //Towards a fully automated calculation of rovibrational infrared intensities for semi-rigid polyatomic molecules//, J. Chem. Phys. **153**, xxxxxx (2020).\\ 
-The following //options// are available: 
- 
-  * **''HOTB''=//n//** (=0 (off) Default) The calculation of vibrational hot bands can be switched on with ''HOTB=1''. 
-  * **''IRDAT''=//string//** File name for dumping the rovibrational infrared line list. Activates calculation of rovibrational intensities. 
-  * **''IRUNIT''=//string//** The default unit for the IR intensities is km/mol. ''IRUNIT''=’HITRAN’ provides the results in HITRAN units, i.e. cm$^{-1}$/(molecule cm$^{-2}$). This key is only active in rovibrational calculations. 
-  * **''JMAX''=//n//** By default VSCF calculations will be performed for non-rotating molecules, i.e. J=0. Rovibrational levels can be computed for arbitrary numbers of J$=n$. This will perform a purely rotational calculation (RCI). To obtain approximate rovibrational energies, vibrational energies have to be added. 
-  * **''%%JMAX\_PRINT%%''=//n//** (=min(Jmax,3) Default) This option controls the printout in rovibrational calculations, i.e. the maximum J value, up to which information shall be printed. 
-  * **''NSSW''=//string//** Sequence of nuclear spin statistical weights in the order of irreps commonly used in the character table for the current molecular point group, e.g. ’1-1-3-3’ for irreps A$_1$, A$_2$, B$_1$, B$_2$ in the case of H$_2$CO. Overwrites automatically determined values. 
-  * **''%%PARTF\_R\_THR%%''=//value//** (=$10^{-4}$ Default) Threshold for the relative deviation within the iterative determination of the rotational partition function. 
-  * **''%%PARTF\_V\_THR%%''=//value//** (=$10^{-2}$ Default) Threshold for the relative deviation within the iterative determination of the vibrational partition function. 
-  * **''RBAS''=//n//** (=1 Default) Definition of the rotational basis in rovibrational calculations.\\ 
-''RBAS=1'' refers to primitive rigid rotor states $|Jk>$.\\ 
-''RBAS=2'' a symmetrized rotational basis by employing Wang combination is used, i.e. $|J K \tau> = i^\sigma/\sqrt{2} (|JK> + (-1)^{J+K+\tau}|J-K>)$. 
-  * **''RVINFO''=//n//** (=1 Default) Additional rovibrational output. By default this will print the nuclear spin statistical weights. ''RVINFO=2'' provides additional details on the calculation and assignment of nuclear spin statstical weights. ''RVINFO=3'' enables further integrals, etc. 
-  * **''%%RVINT\_THR%%''=//value//** (=10$^{-2}$ Default) Threshold for printing rovibrational intensities. 
-  * **''RVMU''=//n//** This keyword controls the order of integrals arising from the inverse moment of inertia tensor $\mu_{\alpha\beta}$ within the calculation of the partition functions as needed in rovibrational calculations. By default a constant $\mu$-tensor is assumed. 
-  * **''RVPRINT''=//n//** This keyword controls the rovibrational line list printout. ''RVPRINT=1'' prints the transition moments, ''RVPRINT=2'' the oscillator strengths, ''RVPRINT=3'' the Einstein A coefficients, ''RVPRINT=4'' symmetry information, and ''RVPRINT=5'' vibrational hot bands. Any of these numbers can be combined, e.g. ''RVPRINT=123'' prints the transition moments, the oscillator strengths and the Einstein A coefficients. This keyword or the ''IRDAT'' and/or ''RAMANDAT'' keyword have to be set in order for rovibrational intensitites to be computed. 
-  * **''TINC''=//value//** (= 100 Default, in K) Temperature increment. 
-  * **''TLIST''=//string//** (off Default, in K) List of specific temperature values, e.g. ’300-350-400’. Combinable with other temperature-keywords. 
-  * **''TMAX''=//value//** (=0 (off) Default, in K) Maximum temperature. Setting only ''TMIN'' will set ''TMAX'' to the same value. 
-  * **''TMIN''=//value//** (=0 (off) Default, in K) Minimum temperature. Setting only ''TMAX'' will set ''TMIN'' to the same value. 
- 
-==== Visualization of results ==== 
- 
-''GRAPH'',//options// 
- 
-The ''GRAPH'' directive enables the visualization of ''VSCF'' and ''VCI'' results, mostly in the form of Gnuplot scripts with corresponding data files. Note, for the more demanding postprocessing of rovibrational line lists, the ''DAT2GR'' program should be used.\\ 
-The following //options// are available: 
- 
-  * **''IR''=//string//** (=0 (off) Default) Once dipole surfaces are provided and IR intensity calculations are performed, an GnuPlot file with filename //string// for an IR spectrum can be created with this option. The file name wlll be the same as for the input file, but with the extension .gnu. 
-  * **''PARTF''=//string//** (=0 (off) Default) Within rovibrational calculations, the convergence of the vibrational and rotational partition functions as a function of the vibrational state and rotational principal quantum number $J$ can be plotted. Creates a Gnuplot file with filename //string//. 
-  * **''PRINTWF''=//n//** This key allows for the printing of the effective potentials and of the //n// lowest modals. 
- 
-==== Record handling ==== 
- 
-''DISK'',//options// 
- 
-The ''DISK'' directive allows to specify explicitly, from where the potential information shall be taken and where it shall be stored to disk. This can also be accomplished in an automated manner. These features are only relevant for the simulation of vibronic spectra as one has to deal with several PESs in the same input. For simple VCI calculations, no information is needed here. 
- 
-The following //options// are available: 
- 
-  * **''AUTO''=//n//** Rather than using the options ''START'' and ''SAVE'' one may simply assign a label //n// to a certain PES and all the records will be set automatically. 
-  * **''SAVE''=//record//** This specifies the record, where to dump the VSCF information. Usually this is the same record as specified in the ''START'' option. Note that the VSCF information is currently stored in the same record as the polynomial information. 
-  * **''START''=//record//** Polynomial and other information shall be read from the specified record. This must be the same record, to which the polynomials have been dumped in the ''POLY'' program. 
-  * **''STARTSURF''=//record//** Surface information shall be read from the specified record. This must correspond to the record, to which the potential has been dumped in the ''SURF'' program. 
- 
-===== The VMCSCF program (VMCSCF) ===== 
- 
-''VMCSCF'',//options// [vmcscf] 
- 
-The ''VMCSCF'' program is still under development and has not yet been fully optimized with respect to speed and efficiency. It is strongly recommended to run a VSCF calculation prior to the VMCSCF calculation and thus to use VSCF modals as an initial guess. By default configuration-selective VMCSCF calculations will be performed. The VMCSCF calculation based on analytical representations of the potential is significantly faster than the grid-based version and should be used whenever possible. By default the active space will be determined automatically, but it can also be controlled by the input stream. For details see:\\ 
-S. Heislbetz, G. Rauhut, //Vibrational multiconfiguration self-consistent field theory: Implementation and test calculations//, [[https://dx.doi.org/10.1063/1.3364861|J. Chem. Phys.]] **132**, 124102 (2010).\\ 
-S. Heislbetz, F. Pfeiffer, G. Rauhut, //Configuration selection within vibrational multiconfiguration self-consistent field theory: Application to bridged lithium compounds//, [[https://dx.doi.org/10.1063/1.3593714|J. Chem. Phys.]] **134**, 204108 (2011).\\ 
-P. Meier, D. Oschetzki, F. Pfeiffer, G. Rauhut, //Towards an automated and efficient calculation of resonating vibrational states based on state-averaged multiconfigurational approaches//, [[https://dx.doi.org/10.1063/1.4938280|J. Chem. Phys.]] **143**, 244111 (2015). [2ex] 
- 
-The following //options// are available: 
- 
-  * **''AVERAGE''=//n//** By default state-specific VMCSCF calculations will be performed. This may be altered by ''AVERAGE=1'', which calls the state-averaged VMCSCF program. The states, which will automatically be chosen for averaging, are displayed in th output. 
-  * **''CALCANGLE''=//n//** The procedure how to determine the rotational angles between the modals within the Jacobi rotations can be specified by this keyword. ''CALCANGLE=1'' uses the standard quadratic procedure, which is the current default. ''CALCANGLE=2'' uses a cubic equation instead. ''CALCANGLE=3'' switches to a full numerical determination of the rotational angle. 
-  * **''CITYPE''=//n//** Defines the maximum number of simultaneous excitations within the configurations generated from modals of the active space, i.e. Singles, Doubles, Triples, ... The maximum excitation level is limited to ''CITYPE=9'', the default is set to 4 for 3D potentials and to 5 for 4D potentials. 
-  * **''GSMODALS''=//n//** By default all VMCSCF calculations are based on ground-state based VSCF modals, ''GSMODALS=1''. ''GSMODALS=0'' uses state-specific modals for all VMCSCF calculations. 
-  * **''MAXMAC''=//n//** Controls the maximum number of macroiterations of the VMCSCF program. The default is 25. 
-  * **''MAXMIC''=//n//** Controls the maximum number of microiterations of the VMCSCF program. The default is 50. 
-  * **''NACT''=//n//** Is the number of active modals for each mode. The smallest meaningful value for ''NACT'' is 3, which is the current default. Configurations will only be generated within this space. Note that there is no equivalent in the VMCSCF program to the closed and core orbitals in electronic structure theory. 
-  * **''NVIRT''=//n//** Is the number of modals on top of the number of active modals. The default is ''NVIRT=3''. The virtual modals are needed for the modal rotations. Modals above the virtual modals are entirely neglected. 
-  * **''PRINT''=//n//** Provides additional information within the VMCSCF iterations, once a value larger than 0 (default) is used. 
-  * **''REF''=//n//** ''REF=0'' specifies ground-state based VMCSCF calculations, i.e. the initial guess of the VMCSCF calculation is given by the VSCF wave function of the vibrational ground state. In contrast to that, ''REF=1'' allows for a state-specific initial guess (default). 
-  * **''ROT''=//n//** Specifies, if only active-virtual modal rotations shall be considered (''ROT=0'', default) or if active-active modal rotations shall be considered as well (''ROT=1''). Formally one would need active-active modal rotations for configuration-selective VMCSCF calculations, but the effects are usually extremely small. 
-  * **''SAVE''=//record//** Specifies the record and file, on which the VMCSCF wave function shall be stored. The default is 5950.2. 
-  * **''SAVEVSCF''=//n//** Once set to 1, this option allows to store the VMCSCF modals in the record of the VSCF modals (which will be overwritten). This allows for VCI calculations with VMCSCF modals. 
-  * **''THRMIC''=//value//** Controls the convergence threshold for the microiterations of the VMCSCF program. The default is 1.d-5. 
-  * **''USERMODE''=//n//** Once vibrational states have been defined with the ''VIBSTATE'' program (section [[vibrational SCF programs#the DAT2GR program (DAT2GR)|the DAT2GR program (DAT2GR)]]), the VMCSCF program can be forced to compute just these states by the option ''USERMODE=1''. Note that the vibrational ground state will always be computed and needs not to be specified explicitly. 
-  * **''VERSION''=//n//** By default configuration-selective VMCSCF calculations will be performed (''VERSION=3''). Note that the selection is always performed prior to the VMCSCF iterations but not within them. The selection of configurations can be switched off by (''VERSION=4''). 
- 
-Besides these VMCSCF specific keywords, a number of option can be used, which are identical with those provided for the VCI program. These keywords are ''%%NDIM, VAM, POT, NBAS, DIPOLE, NDIMDIP, COMBI, INFO, MPG, DIAG, CONT, THRSEL, THRCF, ANALYZE%%''. 
- 
-==== Explicit definition of active spaces ==== 
- 
-''NACT'',//options//\\ 
-Within the VMCSCF program the active space can be specified in a general manner (option ''NACT''), which means that the same number of active modals for each mode will be used. Alternatively, one may use the ''NACT'' directive. This allows to specify active spaces for the individual modes. Typically such definition refer to just one vibrational state and thus the program should be used within the ''USERMODE'' mode. 
- 
-  * **''AUTO''=//on / off//** By default, the active space will be determined by a VMP2 based criterion (''AUTO=on''). This can be switched off by ''AUTO=off''. The options exists only for analytical representations of the potential. 
-  * **''%%MODE(n)%%''=//m//** The number of active modals for mode //n// is set to //m//. This always includes the vibrational ground state. In order to resolve Fermi resonances, one would need to choose //m// at least to be 3, which corresponds to the modals 0, 1 and 2. 
- 
-===== The VIBSTATE program (VIBSTATE) ===== 
- 
-''VIBSTATE'',//options// [vibstate] 
- 
-The ''VIBSTATE'' program allows to specify the occupation number vectors for individual vibrational states to be calculated in the following vibrational SCF and vibration correlation programs. Within the input stream, the ''VIBSTATE'' program needs to be called prior to the first call of the VSCF program. Note that, the ''VIBSTATE'' program needs only to be called, if a limited number of states shall be computed. On the contrary, the ''COMBI'' keyword within the ''VSCF'' and ''VCI'' programs generate large lists of vibrational states, which may result in significant computational effort. The ''VIBSTATE'' program necessarily requests the ''VSTATE'' directive to be called as described below. If only those vibrational states shall be computed, which are defined by the ''VIBSTATE'' program, all subsequent programs need to call the ''USERMODE'' option. 
- 
-==== Definition of vibrational states ==== 
- 
-''VSTATE'',//options//\\ 
-This directive specifies the occupation number vector of the vibrational state to be calculated. 
- 
-  * **''%%LQUANT(n)%%''=//m//** For symmetric top and linear molecules the specification of the quantum number $l$ is supported by the keyword ''LQUANT''. With ''%%LQUANT(n)=l%%'' the $n$th mode is supposed to have the quantum number $l$ corresponding to the angular momentum. The occupation number must be set by the keyword ''%%MODE(n)%%'' before. E.g., ''%%VSTATE,MODE(3)=2,LQUANT(3)=2%%''. All $l$-quantum numbers not set will have a value of zero. Note that this option can only be used for degenerate mode pairs. 
-  * **''%%MODE(n)%%''=//m//** The $n$th mode of the molecule is supposed to have the quantum number //m//. If several modes are excited, the option can be repeated accordingly, e.g. ''%%VSTATE,MODE(1)=3,MODE(3)=2,MODE(6)=2%%''. All other modes are not excited and thus have the quantum number 0. 
- 
-===== The DAT2GR program (DAT2GR) ===== 
- 
-''DAT2GR'',//options// [vibstate] 
- 
-The ''DAT2GR'' program allows for the transformation of rovibrational line lists (IR and Raman) to a file format, which can be further processed by any graphics program. The file containing the line list must be generated before by the ''VSCF'' or ''VCI'' programs. Besides the conversion, it accounts accounts for the temperature dependence of the spectra. The ''DAT2GR'' does not require any ''XSURF'' or ''VSCF'' calculations prior to its call, but requests at least a Hartree-Fock calculation, which is always needed in Molpro.\\ 
-The following //options// are available: 
- 
-  * **''DUMP''=//filename//** This key defines the output file to be generated. 
-  * **''EINC''=//value//** (= 0.1 Default, in cm$^{-1}$) Incremental energy (resolution) to be considered. 
-  * **''EMAX''=//value//** (= 5000 Default, in cm$^{-1}$) This option specifies the maximum energy of the spectral range to be considered. 
-  * **''EMIN''=//value//** (= 0 Default, in cm$^{-1}$) This option specifies the minimum energy of the spectral range to be considered. 
-  * **''EXTERN''=//filename//** This keyword is mandatory and requests the filename of the file containing the input file (rovibrational line list) as generated from the VSCF or VCI programs. 
-  * **''GAMMA''=//value//** Half width to be used for Lorentzians. 
-  * **''%%GAMMA_LIST%%''=//value//** (off Default) Define a list of scaling parametes ''%%GAMMA_SCALE%%'' for the half-width for Lorentzians used in Voigt profile, e.g. ’1-3-10’. 
-  * **''%%GAMMA_SCALE%%''=//value//** (=1.0 (off) Default) Scale the half-width for Lorentzians used in Voigt profile as $\gamma$/''%%GAMMA_SCALE%%''. 
-  * **''PUNIT''=//string//** (=’atm’ Default) Pressure unit. 
-  * **''PINC''=//value//** (= 0.1 atm Default, in ''PUNIT'') Pressure increment to be considered. 
-  * **''PMAX''=//value//** (= 1.0 atm Default, in ''PUNIT'') Maximum pressure to be considered. 
-  * **''PMIN''=//value//** (= 1.0 atm Default, in ''PUNIT'') Minimum pressure to be considered. 
-  * **''PLIST''=//string//** (off Default) List of pressures to be considered, e.g. ’1.0-1.1-1.2’. Combinable with other pressure-keywords. 
-  * **''PRINT''=//n//** (=0 (off) Default) ''PRINT=2'' enables additional printout. 
-  * **''PROFILE''=//string//** (=’Gaussian’ Default) The option controls the line shape to be used for accounting of line-broadening. ''%%PROFILE=’Voigt’%%'' specifies a line shape from numerical convolution of Gaussians and Lorentzians, which is numerically expensive. The combined half-widths $\sigma$ for the Gaussian and $\gamma$ for the Lorentzian line shape are calculated from the temperature and pressure using a physical description. ''%%PROFILE=’Gaussian’%%'' specifies a simple Gaussian line shape with half-width ''SIGMA''. 
-  * **''SIGMA''=//value//** (= 1.0 Default, in cm$^{-1}$) Half width to be used for Gaussians if ''%%PROFILE=’GAUSSIAN’%%''. 
-  * **''%%SIGMA_LIST%%''=//value//** (off Default) Define a list of scaling parametes ''%%SIGMA_SCALE%%'' for the half-width for Gaussians used in Voigt profile, e.g. ’1-3-10’. 
-  * **''%%SIGMA_SCALE%%''=//value//** (=1.0 (off) Default) Scale the half-width for Gaussians used in Voigt profile as $\sigma$/''%%SIGMA_SCALE%%''. 
-  * **''TINC''=//value//** (= 100 Default, in K) Temperature increment to be considered. 
-  * **''TMAX''=//value//** (=0 (off) Default, in K) Maximum temperature to be considered. Setting only ''TMIN'' will set ''TMAX'' to the same value. 
-  * **''TMIN''=//value//** (=0 (off) Default, in K) Minimum temperature to be considered. Setting only ''TMAX'' will set ''TMIN'' to the same value. 
-  * **''TLIST''=//string//** (off Default, in K) List of temperatures to be considered, e.g. ’200-250-300’. Combinable with other temperature-keywords. to be considered. 
- 
-The following example generates a rovibrational IR spectrum for the spectral range between 900 and 1100 cm$^{-1}$ for the given temeprature of 300 K from the line list provided in the file ''H2CS_IR_VCI_J3.dat'' and dumps the output to ''H2CS_IR_VCI_J3_graph.dat''.\\ 
- 
- 
-<code> 
-memory,50,m 
-logfile,scratch 
-mass,iso 
- 
-geometry={ 
-   S 
-   C,S,rcs 
-   H1,C,rch,S,ahcs 
-   H2,C,rch,S,ahcs,H1,180 
-} 
- 
-rcs = 1.2 ang 
-rch = 1.0 ang 
-ahcs = 120 degree 
- 
-basis=vdz 
-hf 
-dat2gr, extern='H2CS_IR_VCI_J3.dat', Emin=900, Emax=1100, Einc=1d-1 
-dat2gr, type='IR', profile='Voigt', dump='H2CS_IR_VCI_J3_graph.dat', TList='300', print=2 
-</code> 
-