diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst
index 6c75ba8ca302355aa556b715e0e35e802f393e4f..206c255c40ba19855b548208ce95821485e3192e 100644
--- a/doc/source/users.guide/programs/vibrot.rst
+++ b/doc/source/users.guide/programs/vibrot.rst
@@ -386,38 +386,103 @@ input:
+:kword:`DISTunit`
+ Unit used for distances in the input potential. The default is `BOHR`. Other
+ options include `ANGSTROM` and `PICOMETER`. The short form `PM` can also be used,
+ instead of `PICOMETER`.
+
+ .. xmldoc::
+ %%Keyword: DISTunit
+
+ Specifies the unit used for distances in the input potential.
+
+
+
+:kword:`ENERunit`
+ Unit used for energies in the input potential. The default is `HARTREE`. Other
+ options include `ELECTRONVOLT`, `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MEGAHERTZ`.
+ The short form `EV` can be used instead of `ELECTRONVOLT` and likewise `MHZ`
+ can be used instead of `MEGAHERTZ`.
+
+ .. xmldoc::
+ %%Keyword: ENERunit
+
+ Specifies the unit used for energies in the input potential.
+
+
+
Input example
.............
::
- &VIBROT
+&VIBROT
RoVibrational spectrum
- Title = Vib-Rot spectrum for FeNi
- Atoms = 0 Fe 0 Ni
+ Title = H2 (^1 Pi_u)
+ Atoms = 0 H 0 H
Potential
- 1.0 -0.516768
- 1.1 -0.554562
- [...]
- Plot = 1.0 10.0 0.1
- Grid = 150
- Range = 1.0 10.0
- Vibrations = 10
- Rotations = 2 10
- Orbital = 2
+ 0.4233417991952784 -93390.8116364055
+ 0.5291772489940979 -125520.5784258792
+ 0.5820949738935077 -135202.0740308874
+ 0.6350126987929174 -142230.7885620708
+ 0.6879304236923273 -147325.2117261678
+ 0.7408481485917370 -150985.4845047687
+ 0.7937658734911469 -153567.9481018878
+ 0.8466835983905567 -155331.6637865382
+ 0.8996013232899664 -156468.2460791877
+ 0.9525190481893763 -157121.6176632051
+ 1.0054367730887860 -157401.2568735270
+ 1.0583544979881960 -157391.4024626400
+ 1.1112722228876060 -157157.4776230008
+ 1.1641899477870150 -156750.6989542662
+ 1.2700253975858350 -155571.7997582064
+ 1.4816962971834740 -152450.7563927988
+ 1.6933671967811130 -149070.0021134733
+ 1.9050380963787530 -145873.2312217305
+ 2.1167089959763920 -143043.6172437684
+ 2.6458862449704900 -137805.7761879516
+ 3.1750634939645880 -134764.6588985511
+ 5.2917724899409790 -131360.0872323780
+ DistUnit = Angstrom
+ EnerUnit = cm-1
+ Grid = 450
+ Range = 0.4 5.0
+ Vibrations = 3
+ Rotations = 1 4
+ Orbital = 1
Observable
- Dipole Moment
- 1.0 0.102354
- 1.1 0.112898
- [...]
+ Dipole Moment
+ 0.4233417991952784 0.57938359
+ 0.5291772489940979 0.62852037
+ 0.5820949738935077 0.65216622
+ 0.6350126987929174 0.67506184
+ 0.6879304236923273 0.69709869
+ 0.7408481485917370 0.71821433
+ 0.7937658734911469 0.73833904
+ 0.8466835983905567 0.75741713
+ 0.8996013232899664 0.77538706
+ 0.9525190481893763 0.79219774
+ 1.0054367730887860 0.80778988
+ 1.0583544979881960 0.82211035
+ 1.1112722228876060 0.83510594
+ 1.1641899477870150 0.84672733
+ 1.2700253975858350 0.86565481
+ 1.4816962971834740 0.88532063
+ 1.6933671967811130 0.88056207
+ 1.9050380963787530 0.85474708
+ 2.1167089959763920 0.81515210
+ 2.6458862449704900 0.70549066
+ 3.1750634939645880 0.62103112
+ 5.2917724899409790 0.46501146
Plot = 1.0 10.0 0.1
+ Scale
-**Comments**: The vibrational-rotation spectrum for :math:`\ce{FeNi}`
-will be computed using the potential curve given in input. The 10
-lowest vibrational levels will be obtained and for each level the
-rotational states in the range :math:`J`\=2 to 10. The vib-rot matrix elements
+**Comments**: The vibrational-rotation spectrum for the :math:`^1\Pi_u` state of \
+:math:`\ce{H2}` will be computed using the potential curve given in the input. The 3
+lowest vibrational levels will be obtained and for each level for the
+rotational states in the range :math:`J`\=1 to 4. The mass for
+the most abundant isotope of :math:`\ce{H}` will be used. The vib-rot matrix elements
of the dipole function will also be computed. A plot file of the
-potential and the dipole function will be generated. The masses for
-the most abundant isotopes of :math:`\ce{Fe}` and :math:`\ce{Ni}` will be selected.
+potential and the dipole function will be generated.
.. xmldoc::
diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90
index d7a9f49cb43601188c580e2f04a9fb3e7476546d..477ed9c21023286148e90856ded05cff2474241b 100644
--- a/src/vibrot/vibinp.F90
+++ b/src/vibrot/vibinp.F90
@@ -19,8 +19,8 @@ subroutine Vibinp(ncase,ngrid,nvib,Umin,Umax,Rout,PotR,E0,dE0,Redm,Teas,Req,sc,t
use Vibrot_globals, only: Atom1, Atom2, dRo, EoutO, iad12, iad13, iadvib, iallrot, IfPrWf, iobs, iplot, iscale, ispc, J1A, J1B, &
J2A, J2B, lambda, n0, n02, nop, npin, npobs, npoint, nRot_Max, nvib1, nvib21, Obsin, R0o, R1o, RinO, &
- Titobs, Vibwvs, Vibwvs1, Vibwvs2
-use Constants, only: Zero, One, Five, UTOAU
+ Titobs, Vibwvs, Vibwvs1, Vibwvs2, DistUnit, EnerUnit
+use Constants, only: Zero, One, Five, UTOAU, Angstrom, auToeV, auTokcalmol, auToHz, auTocm, cal_to_J
use Definitions, only: wp, iwp, u6
implicit none
@@ -41,9 +41,9 @@ character(len=4) :: word, Diatom, Diatomx
character(len=8) :: IntCh
character(len=80) :: Title1(10), Title2(10)
character(len=180) :: Line, l84, l84x
-integer(kind=iwp), parameter :: ntab = 19
+integer(kind=iwp), parameter :: ntab = 21
character(len=*), parameter :: tabinp(ntab) = ['TITL','ATOM','GRID','RANG','VIBR','ROTA','ORBI','NOSP','OBSE','STEP', &
- 'POTE','ROVI','TRAN','ASYM','PRWF','SCAL','TEMP','ALLR','END ']
+ 'POTE','ROVI','TRAN','ASYM','PRWF','SCAL','TEMP','ALLR','DIST','ENER','END ']
integer(kind=iwp), external :: IsFreeUnit, iNuclearChargeFromSymbol, iMostAbundantIsotope
real(kind=wp), external :: dNuclearMass
character(len=180), external :: Get_Ln, Get_Ln_EOF
@@ -55,6 +55,8 @@ call SpoolInp(LuIn)
Atom1 = ''
Atom2 = ''
+DistUnit = 'BOHR'
+EnerUnit = 'HARTREE'
ipot = 0 ! Indicator for potential input
ngrid = 199 ! Maximum number of grid points
Rmin = One
@@ -402,6 +404,144 @@ input: do
iallrot = 1
case (tabinp(19))
+ ! Distance units
+ Line = Get_Ln(LuIn)
+ call Upcase(Line)
+! DistUnit = Line
+!select case (DistUnit)
+ select case (Line)
+
+ case ('BOHR')
+ ! Distance units of Bohr radii, no need for conversion
+ write(u6,*)
+ write(u6,*) 'Distance provided in units of Bohr radii.'
+ write(u6,*) 'No conversion.'
+
+ case ('ANGSTROM')
+ ! Distance units of Angstroms, convert to Bohr radii
+ write(u6,*)
+ write(u6,*) 'Distance provided in units of Angstroms.'
+ write(u6,*) 'Converting to Bohr radii.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Rin(i) = Rin(i) / Angstrom
+ end do
+ do i = 1, iobs
+ do j = 1, npobs(i)
+ RinO(j,i) = RinO(j,i) / Angstrom
+ end do
+ end do
+ end if
+
+ case ('PICOMETER','PM')
+ ! Distance units of picometers, convert to Bohr radii
+ write(u6,*)
+ write(u6,*) 'Distance provided in units of picometers.'
+ write(u6,*) 'Converting to Bohr radii.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Rin(i) = Rin(i) * 1.0e-2_wp / Angstrom
+ end do
+ do i = 1, iobs
+ do j = 1, npobs(i)
+ RinO(j,i) = RinO(j,i) * 1.0e-2_wp / Angstrom
+ end do
+ end do
+ end if
+
+ case default
+ write(u6,*)
+ write(u6,*) '********************************************'
+ write(u6,*) ' VIBINP Error: Distance unit not recognized.'
+ write(u6,*) '********************************************'
+ call Quit_OnUserError()
+ end select
+
+ case (tabinp(20))
+ ! Energy units
+ Line = Get_Ln(LuIn)
+ call Upcase(Line)
+! EnerUnit = Line
+! select case (EnerUnit)
+ select case (Line)
+
+ case ('HARTREE')
+ ! Energy units of hartrees, no need for conversion
+ write(u6,*)
+ write(u6,*) 'Energy provided in units of hartree.'
+ write(u6,*) 'No conversion.'
+
+ case ('EV','ELECTRONVOLT')
+ ! Energy units of eV, convert to hartrees
+ write(u6,*)
+ write(u6,*) 'Energy provided in electronvolts.'
+ write(u6,*) 'Converting to hartree.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Ein(i) = Ein(i) / auToeV
+ end do
+ end if
+
+ case ('KCAL/MOL')
+ ! Energy units of kcal/mol, convert to hartrees
+ write(u6,*)
+ write(u6,*) 'Energy provided in kcal/mol.'
+ write(u6,*) 'Converting to hartree.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Ein(i) = Ein(i) / auTokcalmol
+ end do
+ end if
+
+ case ('KJ/MOL')
+ ! Energy units of kJ/mol, convert to hartrees
+ write(u6,*)
+ write(u6,*) 'Energy provided in kJ/mol.'
+ write(u6,*) 'Converting to hartree.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Ein(i) = Ein(i) / (auTokcalmol * cal_to_J)
+ end do
+ end if
+
+ case ('CM-1')
+ ! Energy units of cm^(-1), convert to hartrees
+ write(u6,*)
+ write(u6,*) 'Energy provided in cm^(-1).'
+ write(u6,*) 'Converting to hartree.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Ein(i) = Ein(i) / auTocm
+ end do
+ end if
+
+ case ('MHZ','MEGAHERTZ')
+ ! Energy units of MHz, convert to hartrees
+ write(u6,*)
+ write(u6,*) 'Energy provided in MHz.'
+ write(u6,*) 'Converting to hartree.'
+
+ if (ipot /= 0) then
+ do i=1,nop
+ Ein(i) = Ein(i) * 1.0e6_wp / auToHz
+ end do
+ end if
+
+ case default
+ write(u6,*)
+ write(u6,*) '******************************************'
+ write(u6,*) ' VIBINP Error: Energy unit not recognized.'
+ write(u6,*) '******************************************'
+ call Quit_OnUserError()
+ end select
+
+ case (tabinp(21))
exit input
case default
diff --git a/src/vibrot/vibrot_globals.F90 b/src/vibrot/vibrot_globals.F90
index 97bc081864b1bf37bf88d0d8266aa453af10b544..e6009805e4f29ba62db685fa58f48ca50e606316 100644
--- a/src/vibrot/vibrot_globals.F90
+++ b/src/vibrot/vibrot_globals.F90
@@ -22,10 +22,11 @@ integer(kind=iwp) :: J1A, J2A, lambda, ispc, iobs, nop, Vibwvs, iadvib, Vibwvs1,
integer(kind=iwp) :: iadrsp(nRot_Max), iad12(nRot_Max), iad13(nRot_Max), iplot(nobs), npobs(nobs)
real(kind=wp) :: R0o(nobs), R1o(nobs), dRo(nobs), RinO(npin,nobs), Obsin(npin,nobs), EoutO(npoint+4,nobs)
character(len=2) :: Atom1, Atom2
+character(len=9) :: DistUnit, EnerUnit
character(len=80) :: Titobs(nobs)
public :: Atom1, Atom2, dRo, EoutO, iad12, iad13, iadrsp, iadvib, iallrot, ifPrWf, iobs, iplot, iscale, ispc, J1A, J1B, J2A, J2B, &
lambda, n0, n02, nop, npin, npobs, npoint, nRot_Max, nvib1, nvib21, Obsin, R0o, R1o, RinO, Titobs, Vibwvs, Vibwvs1, &
- Vibwvs2
+ Vibwvs2, DistUnit, EnerUnit
end module Vibrot_globals
diff --git a/test/standard/016.input b/test/standard/016.input
index 2049b9e79737c4a3665e7b7896076b614f85e4de..36fbc30cd83b3aab8346fba606fabcdfc41abe3d 100644
--- a/test/standard/016.input
+++ b/test/standard/016.input
@@ -6,11 +6,49 @@
* Responsible person: P.-O. Widmark 041227
* Comments: Test of VIBROT, spectroscopic constants and trans. mom.
*-------------------------------------------------------------------------------
-* Ground state, 1Sg+
+* Ground state, ^1Sigma_g+ (spectroscopic units)
*-------------------------------------------------------------------------------
&VIBROT
RoVibrational spectrum
- Title = H2 1Sg+
+ Title = H2 1Sg+ (spectroscopic units)
+ Atoms = 0 H 0 H
+ Potential
+ 0.4233417991952784 -223560.5452917840
+ 0.5291772489940979 -246545.0303080266
+ 0.5820949738935077 -252161.7065226864
+ 0.6350126987929174 -255439.1694568950
+ 0.6879304236923273 -257075.3967184710
+ 0.7408481485917370 -257550.2147172519
+ 0.7937658734911469 -257200.7606271270
+ 0.8466835983905567 -256268.3533880202
+ 0.8996013232899664 -254928.5024720499
+ 0.9525190481893763 -253310.5640313918
+ 1.0054367730887860 -251510.9071813326
+ 1.0583544979881960 -249601.9124601000
+ 1.1112722228876060 -247638.2817244752
+ 1.1641899477870150 -245661.5110427523
+ 1.2700253975858350 -241787.2175787069
+ 1.4816962971834740 -234849.6267191532
+ 1.6933671967811130 -229417.8271538202
+ 1.9050380963787530 -225549.0973716453
+ 2.1167089959763920 -223014.0578525766
+ 2.6458862449704900 -220281.0390487701
+ 3.1750634939645880 -219644.2268255862
+ 5.2917724899409790 -219468.1708616391
+ DistUnit = ANGSTROM
+ EnerUnit = CM-1
+ Grid = 450
+ Range = 0.4 5.0
+ Vibrations = 3
+ Rotations = 0 3
+ Orbital = 0
+ Scale
+*-------------------------------------------------------------------------------
+* Ground state, 1Sg+ (default units)
+*-------------------------------------------------------------------------------
+&VIBROT
+ RoVibrational spectrum
+ Title = H2 ^1Sigma_g+
Atoms = 0 H 0 H
Potential
0.800 -1.01861680
@@ -35,6 +73,8 @@
5.000 -1.00367427
6.000 -1.00077274
10.00 -0.99997057
+ DistUnit = BOHR
+ EnerUnit = HARTREE
Grid = 450
Range = 0.4 5.0
Vibrations = 3
@@ -44,7 +84,7 @@
*---
>>COPY $Project.VibWvs $Project.VibWvs_GS
*-------------------------------------------------------------------------------
-* Check ground state
+* Excited state, ^1Pi_u (default units)
*-------------------------------------------------------------------------------
&VIBROT
RoVibrational spectrum
@@ -82,12 +122,10 @@
*---
>>COPY $Project.VibWvs $Project.VibWvs_XS
*-------------------------------------------------------------------------------
-* Check excited state
-*-------------------------------------------------------------------------------
>>COPY $Project.VibWvs_GS VIBWVS1
>>COPY $Project.VibWvs_XS VIBWVS2
*-------------------------------------------------------------------------------
-* Transition moments
+* Transition moments (default units)
*-------------------------------------------------------------------------------
&VIBROT
Transition moments
@@ -120,14 +158,16 @@
>>FILE checkfile
* This file is autogenerated:
-* Molcas version 20.10-708-gb8344319
-* Linux otis 4.15.0-1073-oem #83-Ubuntu SMP Mon Feb 17 11:21:18 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux
-* Mon Jan 25 13:28:00 2021
+* Molcas version 22.10-257-g3f64bb863
+* Linux cedar5.cedar.computecanada.ca 3.10.0-1160.53.1.el7.x86_64 #1 SMP Fri Jan 14 13:59:45 UTC 2022 x86_64 GNU/Linux
+* Sun Dec 18 21:56:09 2022
*
#>> 1
-#> VIBROT_SPECTC="23070.121186095137"/2
+#> VIBROT_SPECTC="23070.119530335731"/2
#>> 2
-#> VIBROT_SPECTC="13157.186971543833"/2
+#> VIBROT_SPECTC="23070.121186095137"/2
#>> 3
+#> VIBROT_SPECTC="13157.186971543833"/2
+#>> 4
#> VIBROT_VIBTRM="0.963221285781"/6
>>EOF