From bc06d5e78a81d940cda8e189c25f169d613426b6 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Fri, 11 Nov 2022 02:40:51 +0100 Subject: [PATCH 01/42] First stab at unit conversion --- src/vibrot/vibinp.F90 | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 0ecdaeb929..a13f481888 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -84,6 +84,8 @@ call RdNLst(LuIn,'VibRot') ! Read input data from input file +write(*,*)'Deyan molcas' + ntit1 = 0 skip = .false. input: do @@ -577,6 +579,13 @@ if ((ipot == 0) .and. (ncase == 1)) then call Quit_OnUserError() end if +if (ipot /= 0) then + do i=1,nop + Rin(i) = Rin(i) * 1.8897259886 + Ein(i) = Ein(i) * 0.0000045563352812122295 + end do +end if + ! Print input potential if (ipot /= 0) then -- GitLab From e8b4167a405a66338a948f6670895df67fe13464 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Sun, 13 Nov 2022 23:01:39 +0100 Subject: [PATCH 02/42] Include distance and energy unit keywords --- src/vibrot/vibinp.F90 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index a13f481888..927d6e11cb 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -32,7 +32,7 @@ integer(kind=iwp) :: LuIn, LuIn1, ntit1, ntit2, i, j, k, ii, kk, iadvi1, iadvi2, logical(kind=iwp) :: skip, exists real(kind=wp) :: del, dRp, O0, Oeq, R0p, R1p, Redm1, Redmx, Reqx, Rmax, Rmin, Rn1, Rout0, U, Umaxx, Uminx, xMass1, xMass2, xxx, & Rin(npin), Ein(npin) -character(len=2) :: At1x, At2x +character(len=2) :: At1x, At2x, DistUnit, EnerUnit character(len=4) :: word, Diatom, Diatomx ! For storing character data using gather/scatter DAFILE operations, it ! is imperative that the strings are aligned on integers. @@ -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=r8), external :: dNuclearMass character(len=180), external :: Get_Ln, Get_Ln_EOF @@ -404,6 +404,16 @@ input: do iallrot = 1 case (tabinp(19)) + ! Distance units + Line = Get_Ln(LuIn) + call Get_S(1,DistUnit(1:2),1) + + case (tabinp(20)) + ! Energy units + Line = Get_Ln(LuIn) + call Get_S(1,EnerUnit(1:2),1) + + case (tabinp(21)) exit input case default -- GitLab From bc5745ec650c02f061b33df58c6d5df06a4a6798 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Mon, 14 Nov 2022 02:05:30 +0100 Subject: [PATCH 03/42] Distance and energy units conversion --- src/vibrot/vibinp.F90 | 137 +++++++++++++++++++++++++++++++--- src/vibrot/vibrot_globals.F90 | 4 +- 2 files changed, 127 insertions(+), 14 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 927d6e11cb..34f403ac26 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -19,7 +19,7 @@ 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 + Titobs, Vibwvs, Vibwvs1, Vibwvs2, DistUnit, EnerUnit use Constants, only: Zero, One, Five, UTOAU use Definitions, only: wp, iwp, r8, u6 @@ -32,7 +32,7 @@ integer(kind=iwp) :: LuIn, LuIn1, ntit1, ntit2, i, j, k, ii, kk, iadvi1, iadvi2, logical(kind=iwp) :: skip, exists real(kind=wp) :: del, dRp, O0, Oeq, R0p, R1p, Redm1, Redmx, Reqx, Rmax, Rmin, Rn1, Rout0, U, Umaxx, Uminx, xMass1, xMass2, xxx, & Rin(npin), Ein(npin) -character(len=2) :: At1x, At2x, DistUnit, EnerUnit +character(len=2) :: At1x, At2x character(len=4) :: word, Diatom, Diatomx ! For storing character data using gather/scatter DAFILE operations, it ! is imperative that the strings are aligned on integers. @@ -55,6 +55,8 @@ call SpoolInp(LuIn) Atom1 = '' Atom2 = '' +DistUnit = 0 +EnerUnit = 0 ipot = 0 ! Indicator for potential input ngrid = 199 ! Maximum number of grid points Rmin = One @@ -84,8 +86,6 @@ call RdNLst(LuIn,'VibRot') ! Read input data from input file -write(*,*)'Deyan molcas' - ntit1 = 0 skip = .false. input: do @@ -406,12 +406,12 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) - call Get_S(1,DistUnit(1:2),1) + call Get_I1(1,DistUnit) case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) - call Get_S(1,EnerUnit(1:2),1) + call Get_I1(1,EnerUnit) case (tabinp(21)) exit input @@ -589,12 +589,125 @@ if ((ipot == 0) .and. (ncase == 1)) then call Quit_OnUserError() end if -if (ipot /= 0) then - do i=1,nop - Rin(i) = Rin(i) * 1.8897259886 - Ein(i) = Ein(i) * 0.0000045563352812122295 - end do -end if +! Distance units + +select case (DistUnit) + + case (0) + ! 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 (1) + ! 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) * 1.8897259886 + end do + end if + + case (2) + ! 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) * 0.0188973 + 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 + +! Energy units + +select case (EnerUnit) + +case (0) + ! Energy units of hartrees, no need for conversion + write(u6,*) + write(u6,*) 'Energy provided in units of hartrees.' + write(u6,*) 'No conversion.' + +case (1) + ! Energy units of eV, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of electron volts.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.0367493 + end do + end if + +case (2) + ! Energy units of kcal/mol, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of kcal/mol.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.00159360264 + end do + end if + +case (3) + ! Energy units of kJ/mol, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of kJ/mol.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.00038087983 + end do + end if + +case (4) + ! Energy units of cm^(-1), convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of cm^(-1).' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.00000455633 + end do + end if + +case (5) + ! Energy units of MHz, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of MHz.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 1.519829E-10 + 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 ! Print input potential diff --git a/src/vibrot/vibrot_globals.F90 b/src/vibrot/vibrot_globals.F90 index 97bc081864..dbe0f5d993 100644 --- a/src/vibrot/vibrot_globals.F90 +++ b/src/vibrot/vibrot_globals.F90 @@ -18,7 +18,7 @@ private integer(kind=iwp), parameter :: nRot_Max = 200, nobs = 10, npin = 500, npoint = 5000 integer(kind=iwp) :: J1A, J2A, lambda, ispc, iobs, nop, Vibwvs, iadvib, Vibwvs1, Vibwvs2, n0, nvib1, n02, nvib21, J1B, J2B, & - IfPrWf, iscale, iallrot + IfPrWf, iscale, iallrot, DistUnit, EnerUnit 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 @@ -26,6 +26,6 @@ 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 -- GitLab From b73ffaf9e61aed3c8641727d3c64f03c5560a3b8 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 12 Dec 2022 17:03:18 -0500 Subject: [PATCH 04/42] Updated the example input in the vibrot docs because: - The previous input file did not have a complete potential, which made it impossible for anyone to use, without generating their own FeNi potential (a super inconvenience). This means people can't do a quick test to see what the output of the program is supposed to look like, without doing substantial work. - The new input file contains the new feature which allows people to use units that are more common in spectroscopy (which is where vibrot is most likely to be useful). - The new input also matches test/standard/016.input which is convenient. --- doc/source/users.guide/programs/vibrot.rst | 53 ++++++++++++++-------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 6c75ba8ca3..8d98d4031c 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -392,25 +392,40 @@ Input example :: &VIBROT - RoVibrational spectrum - Title = Vib-Rot spectrum for FeNi - Atoms = 0 Fe 0 Ni - 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 - Observable - Dipole Moment - 1.0 0.102354 - 1.1 0.112898 - [...] - Plot = 1.0 10.0 0.1 + RoVibrational spectrum + Title = H2 1Sg+ + 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 **Comments**: The vibrational-rotation spectrum for :math:`\ce{FeNi}` will be computed using the potential curve given in input. The 10 -- GitLab From 7be6cc7607f78d575c5a2889721c512f151fd447 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Mon, 19 Dec 2022 00:02:46 +0100 Subject: [PATCH 05/42] Use constants from constants.F90 --- src/vibrot/vibinp.F90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 34f403ac26..021f0d6a21 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -20,7 +20,7 @@ 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, DistUnit, EnerUnit -use Constants, only: Zero, One, Five, UTOAU +use Constants, only: Zero, One, Five, UTOAU, Angstrom, auToeV, auTokcalmol, auToHz, auTocm, cal_to_J use Definitions, only: wp, iwp, r8, u6 implicit none @@ -607,7 +607,7 @@ select case (DistUnit) if (ipot /= 0) then do i=1,nop - Rin(i) = Rin(i) * 1.8897259886 + Rin(i) = Rin(i) / Angstrom end do end if @@ -619,7 +619,7 @@ select case (DistUnit) if (ipot /= 0) then do i=1,nop - Rin(i) = Rin(i) * 0.0188973 + Rin(i) = Rin(i) * 1.0d-2 / Angstrom end do end if @@ -644,60 +644,60 @@ case (0) case (1) ! Energy units of eV, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of electron volts.' + write(u6,*) 'Energy provided in units of electron volts.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.0367493 + Ein(i) = Ein(i) / auToeV end do end if case (2) ! Energy units of kcal/mol, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of kcal/mol.' + write(u6,*) 'Energy provided in units of kcal/mol.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.00159360264 + Ein(i) = Ein(i) / auTokcalmol end do end if case (3) ! Energy units of kJ/mol, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of kJ/mol.' + write(u6,*) 'Energy provided in units of kJ/mol.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.00038087983 + Ein(i) = Ein(i) / (auTokcalmol * cal_to_J) end do end if case (4) ! Energy units of cm^(-1), convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of cm^(-1).' + write(u6,*) 'Energy provided in units of cm^(-1).' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.00000455633 + Ein(i) = Ein(i) / auTocm end do end if case (5) ! Energy units of MHz, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of MHz.' + write(u6,*) 'Energy provided in units of MHz.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 1.519829E-10 + Ein(i) = Ein(i) * 1.0d6 / auToHz end do end if -- GitLab From cec0a197b077af4d825e5fad467a8dd6241deb1c Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Mon, 19 Dec 2022 00:36:35 +0100 Subject: [PATCH 06/42] Add support for units provided by strings --- src/vibrot/vibinp.F90 | 26 +++++++++++++------------- src/vibrot/vibrot_globals.F90 | 3 ++- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 021f0d6a21..b1c1beeac7 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -55,8 +55,8 @@ call SpoolInp(LuIn) Atom1 = '' Atom2 = '' -DistUnit = 0 -EnerUnit = 0 +DistUnit = '' +EnerUnit = '' ipot = 0 ! Indicator for potential input ngrid = 199 ! Maximum number of grid points Rmin = One @@ -406,12 +406,12 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) - call Get_I1(1,DistUnit) + DistUnit = Line case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) - call Get_I1(1,EnerUnit) + EnerUnit = Line case (tabinp(21)) exit input @@ -593,13 +593,13 @@ end if select case (DistUnit) - case (0) + 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 (1) + case ('ANGSTROM') ! Distance units of Angstroms, convert to Bohr radii write(u6,*) write(u6,*) 'Distance provided in units of Angstroms.' @@ -611,7 +611,7 @@ select case (DistUnit) end do end if - case (2) + case ('PICOMETER') ! Distance units of picometers, convert to Bohr radii write(u6,*) write(u6,*) 'Distance provided in units of picometers.' @@ -635,13 +635,13 @@ end select select case (EnerUnit) -case (0) +case ('', 'HARTREE') ! Energy units of hartrees, no need for conversion write(u6,*) write(u6,*) 'Energy provided in units of hartrees.' write(u6,*) 'No conversion.' -case (1) +case ('EV') ! Energy units of eV, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of electron volts.' @@ -653,7 +653,7 @@ case (1) end do end if -case (2) +case ('KCAL/MOL') ! Energy units of kcal/mol, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of kcal/mol.' @@ -665,7 +665,7 @@ case (2) end do end if -case (3) +case ('KJ/MOL') ! Energy units of kJ/mol, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of kJ/mol.' @@ -677,7 +677,7 @@ case (3) end do end if -case (4) +case ('CM-1') ! Energy units of cm^(-1), convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of cm^(-1).' @@ -689,7 +689,7 @@ case (4) end do end if -case (5) +case ('MHZ') ! Energy units of MHz, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of MHz.' diff --git a/src/vibrot/vibrot_globals.F90 b/src/vibrot/vibrot_globals.F90 index dbe0f5d993..e6009805e4 100644 --- a/src/vibrot/vibrot_globals.F90 +++ b/src/vibrot/vibrot_globals.F90 @@ -18,10 +18,11 @@ private integer(kind=iwp), parameter :: nRot_Max = 200, nobs = 10, npin = 500, npoint = 5000 integer(kind=iwp) :: J1A, J2A, lambda, ispc, iobs, nop, Vibwvs, iadvib, Vibwvs1, Vibwvs2, n0, nvib1, n02, nvib21, J1B, J2B, & - IfPrWf, iscale, iallrot, DistUnit, EnerUnit + IfPrWf, iscale, iallrot 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, & -- GitLab From e89e88a6bb068a5a319d169f881513f6f6901c02 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Sun, 18 Dec 2022 23:55:04 -0500 Subject: [PATCH 07/42] These two things: - added descriptions for the DistUnit and EnerUnit keywords. - Replaced comment about the example input file, with a new one which reflects the new example input file being used. --- doc/source/users.guide/programs/vibrot.rst | 34 +++++++++++++++++----- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 8d98d4031c..23c027e952 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -386,6 +386,28 @@ input: +:kword:`DISTunit` + Unit used for distances in the input potential. The default is `BOHR`. Other + options include `ANGSTOM` and `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 `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. + + .. xmldoc:: + %%Keyword: ENERunit + + Specifies the unit used for energies in the input potential. + + + Input example ............. @@ -427,12 +449,10 @@ Input example Orbital = 0 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 -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. +**Comments**: The vibrational-rotation spectrum for :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`\=0 to 3. The mass for +the most abundant isotope of :math:`\ce{H}` will be used. .. xmldoc:: -- GitLab From 3f64bb86379797e80e74c36f1b3d39170ffc6af7 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 19 Dec 2022 00:09:34 -0500 Subject: [PATCH 08/42] These things: - Corrected typos in test/standard/016.input (e.g. ground state -> excited state, and aomic term symbols -> molecular term symbols) - For the ground state case alraedy in 016.input, I added the default unit specifications - For the ground state case alraedy in 016.input, I made a copy of it and used spectroscopic units instead of atomic units. --- test/standard/016.input | 50 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/test/standard/016.input b/test/standard/016.input index 2049b9e797..c7704c1b0f 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 -- GitLab From 69850b900d61c88726ee4c6169c7263942a0ab74 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 19 Dec 2022 01:09:30 -0500 Subject: [PATCH 09/42] Updated test/standard: - 016.input now has a check file appended to the end of it, after running "pymolcas verify 016 --generate" on the updated 016.input --- test/standard/016.input | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/standard/016.input b/test/standard/016.input index c7704c1b0f..36fbc30cd8 100644 --- a/test/standard/016.input +++ b/test/standard/016.input @@ -158,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 -- GitLab From 7ddf8308e96d22a11919c9313c0683bde4287e52 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 09:56:29 -0500 Subject: [PATCH 10/42] Resolved the following vibrot documentaiton issues: - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890430 (typo in the word ANGSTROM in documentation) - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890430 (use of "SINGLE" instead of "CHOICE" in documentation) - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890433 (Added back the Observable and Plot keywords to the example input) --- doc/source/users.guide/programs/vibrot.rst | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 23c027e952..f7ce7d5fe8 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -388,9 +388,10 @@ input: :kword:`DISTunit` Unit used for distances in the input potential. The default is `BOHR`. Other - options include `ANGSTOM` and `PICOMETER`. + options include `ANGSTROM` and `PICOMETER`. The short form `PM` can also be used, + instead of `PICOMETER`. - .. xmldoc:: + .. xmldoc:: %%Keyword: DISTunit Specifies the unit used for distances in the input potential. @@ -401,7 +402,7 @@ input: Unit used for energies in the input potential. The default is `HARTREE`. Other options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. - .. xmldoc:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. @@ -447,12 +448,20 @@ Input example Vibrations = 3 Rotations = 0 3 Orbital = 0 + Observable + Dipole Moment + 1.0 0.102354 + 1.1 0.112898 + [...] + Plot = 1.0 10.0 0.1 Scale **Comments**: The vibrational-rotation spectrum for :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`\=0 to 3. The mass for -the most abundant isotope of :math:`\ce{H}` will be used. +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. .. xmldoc:: -- GitLab From db5ceb477abce14ebcfefb70b322a45352e06846 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 10:25:43 -0500 Subject: [PATCH 11/42] Fixed the following issues in vibinp (withiin vibrot): - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890442 (defaults are now set to BOHR and HARTREE instead of '' and '') - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890451 call UpCase(Line) was added for DistUnit and EnerUnit - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890454 (no longer explicitly accepting '' as input for DistUnit and EnerUnit) - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890458 (added ability to allow 'PM" for picometers and 'MEGAHERZ' for MHz) --- src/vibrot/vibinp.F90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index b1c1beeac7..8a98498785 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -55,8 +55,8 @@ call SpoolInp(LuIn) Atom1 = '' Atom2 = '' -DistUnit = '' -EnerUnit = '' +DistUnit = 'BOHR' +EnerUnit = 'HARTREE' ipot = 0 ! Indicator for potential input ngrid = 199 ! Maximum number of grid points Rmin = One @@ -406,11 +406,13 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) + call Upcase(line) DistUnit = Line case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) + call Upcase(line) EnerUnit = Line case (tabinp(21)) @@ -593,7 +595,7 @@ end if select case (DistUnit) - case ('', 'BOHR') + case ('BOHR') ! Distance units of Bohr radii, no need for conversion write(u6,*) write(u6,*) 'Distance provided in units of Bohr radii.' @@ -611,7 +613,7 @@ select case (DistUnit) end do end if - case ('PICOMETER') + case ('PICOMETER','PM') ! Distance units of picometers, convert to Bohr radii write(u6,*) write(u6,*) 'Distance provided in units of picometers.' @@ -635,7 +637,7 @@ end select select case (EnerUnit) -case ('', 'HARTREE') +case ('HARTREE') ! Energy units of hartrees, no need for conversion write(u6,*) write(u6,*) 'Energy provided in units of hartrees.' @@ -689,7 +691,7 @@ case ('CM-1') end do end if -case ('MHZ') +case ('MHZ','MEGAHERZ') ! Energy units of MHz, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of MHz.' -- GitLab From b7d5f5a03fb04ecee7161b96a3e0e531e41e990a Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 10:29:28 -0500 Subject: [PATCH 12/42] Fixed typo for Megahertz in vibinp.F90 within vibrot --- src/vibrot/vibinp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 8a98498785..73efe2d51f 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -691,7 +691,7 @@ case ('CM-1') end do end if -case ('MHZ','MEGAHERZ') +case ('MHZ','MEGAHERTZ') ! Energy units of MHz, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of MHz.' -- GitLab From 4fe9ba24744dea3fd1d13ed1806ea40517f41ca1 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 10:31:52 -0500 Subject: [PATCH 13/42] Just made capitalization of L in Line consistent for: DistUnit and EnerUnit parts of vibinp.F90 within vibrot --- src/vibrot/vibinp.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 73efe2d51f..e8b2e71568 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -406,13 +406,13 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) - call Upcase(line) + call Upcase(Line) DistUnit = Line case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) - call Upcase(line) + call Upcase(Line) EnerUnit = Line case (tabinp(21)) -- GitLab From 2011d95dce25b1fffd2766f341782f8cbfeee16a Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 11:04:05 -0500 Subject: [PATCH 14/42] Added to vibrot documentation: - full dipole moment function for H2 --- doc/source/users.guide/programs/vibrot.rst | 27 ++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index f7ce7d5fe8..ea6834f312 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -449,10 +449,29 @@ Input example Rotations = 0 3 Orbital = 0 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 -- GitLab From 139d0c69eb37c093894534590218058ba9708290 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Fri, 23 Dec 2022 13:19:13 +0100 Subject: [PATCH 15/42] First pass at converting distance units for observables --- src/vibrot/vibinp.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index e8b2e71568..6deec9a184 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -611,6 +611,9 @@ select case (DistUnit) do i=1,nop Rin(i) = Rin(i) / Angstrom end do + do i=1,iobs + RinO(1,i) = RinO(1,i) / Angstrom + end do end if case ('PICOMETER','PM') @@ -623,6 +626,9 @@ select case (DistUnit) do i=1,nop Rin(i) = Rin(i) * 1.0d-2 / Angstrom end do + do i=1,iobs + RinO(1,i) = RinO(1,i) * 1.0d-2 / Angstrom + end do end if case default -- GitLab From c19706c8ec9b42e623255e70108a1ed9abcdd8f1 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Fri, 23 Dec 2022 14:36:19 +0100 Subject: [PATCH 16/42] Fix decimal notation --- src/vibrot/vibinp.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 6deec9a184..a339845d1d 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -624,10 +624,10 @@ select case (DistUnit) if (ipot /= 0) then do i=1,nop - Rin(i) = Rin(i) * 1.0d-2 / Angstrom + Rin(i) = Rin(i) * 1.0e-2_wp / Angstrom end do do i=1,iobs - RinO(1,i) = RinO(1,i) * 1.0d-2 / Angstrom + RinO(1,i) = RinO(1,i) * 1.0e-2_wp / Angstrom end do end if @@ -705,7 +705,7 @@ case ('MHZ','MEGAHERTZ') if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 1.0d6 / auToHz + Ein(i) = Ein(i) * 1.0e6_wp / auToHz end do end if -- GitLab From 9aeee388e58a6bb229720b924a0e2101fa1f9f09 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Fri, 23 Dec 2022 15:03:36 -0500 Subject: [PATCH 17/42] Added LIST to vibrot.rst for: DistUnit and EnerUnit keywords as explained here: https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1220348740 --- doc/source/users.guide/programs/vibrot.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index ea6834f312..cd34a604fb 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -391,7 +391,7 @@ input: options include `ANGSTROM` and `PICOMETER`. The short form `PM` can also be used, instead of `PICOMETER`. - .. xmldoc:: + .. xmldoc:: %%Keyword: DISTunit Specifies the unit used for distances in the input potential. @@ -402,7 +402,7 @@ input: Unit used for energies in the input potential. The default is `HARTREE`. Other options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. - .. xmldoc:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. -- GitLab From 32fd63395393c723850e8be1bb0bf68a35f68496 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Fri, 11 Nov 2022 02:40:51 +0100 Subject: [PATCH 18/42] First stab at unit conversion --- src/vibrot/vibinp.F90 | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index d7a9f49cb4..a04f6c16ff 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -84,6 +84,8 @@ call RdNLst(LuIn,'VibRot') ! Read input data from input file +write(*,*)'Deyan molcas' + ntit1 = 0 skip = .false. input: do @@ -577,6 +579,13 @@ if ((ipot == 0) .and. (ncase == 1)) then call Quit_OnUserError() end if +if (ipot /= 0) then + do i=1,nop + Rin(i) = Rin(i) * 1.8897259886 + Ein(i) = Ein(i) * 0.0000045563352812122295 + end do +end if + ! Print input potential if (ipot /= 0) then -- GitLab From 41ef470fe80d4a23d7158ca73800a9ecd46fc07d Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Sun, 13 Nov 2022 23:01:39 +0100 Subject: [PATCH 19/42] Include distance and energy unit keywords --- src/vibrot/vibinp.F90 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index a04f6c16ff..12d6420047 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -32,7 +32,7 @@ integer(kind=iwp) :: LuIn, LuIn1, ntit1, ntit2, i, j, k, ii, kk, iadvi1, iadvi2, logical(kind=iwp) :: skip, exists real(kind=wp) :: del, dRp, O0, Oeq, R0p, R1p, Redm1, Redmx, Reqx, Rmax, Rmin, Rn1, Rout0, U, Umaxx, Uminx, xMass1, xMass2, xxx, & Rin(npin), Ein(npin) -character(len=2) :: At1x, At2x +character(len=2) :: At1x, At2x, DistUnit, EnerUnit character(len=4) :: word, Diatom, Diatomx ! For storing character data using gather/scatter DAFILE operations, it ! is imperative that the strings are aligned on integers. @@ -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 @@ -404,6 +404,16 @@ input: do iallrot = 1 case (tabinp(19)) + ! Distance units + Line = Get_Ln(LuIn) + call Get_S(1,DistUnit(1:2),1) + + case (tabinp(20)) + ! Energy units + Line = Get_Ln(LuIn) + call Get_S(1,EnerUnit(1:2),1) + + case (tabinp(21)) exit input case default -- GitLab From 27b1b47ec5ea70613ab902cec93ceb997ed0dbae Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Mon, 14 Nov 2022 02:05:30 +0100 Subject: [PATCH 20/42] Distance and energy units conversion --- src/vibrot/vibinp.F90 | 137 +++++++++++++++++++++++++++++++--- src/vibrot/vibrot_globals.F90 | 4 +- 2 files changed, 127 insertions(+), 14 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 12d6420047..bd7a293846 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -19,7 +19,7 @@ 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 + Titobs, Vibwvs, Vibwvs1, Vibwvs2, DistUnit, EnerUnit use Constants, only: Zero, One, Five, UTOAU use Definitions, only: wp, iwp, u6 @@ -32,7 +32,7 @@ integer(kind=iwp) :: LuIn, LuIn1, ntit1, ntit2, i, j, k, ii, kk, iadvi1, iadvi2, logical(kind=iwp) :: skip, exists real(kind=wp) :: del, dRp, O0, Oeq, R0p, R1p, Redm1, Redmx, Reqx, Rmax, Rmin, Rn1, Rout0, U, Umaxx, Uminx, xMass1, xMass2, xxx, & Rin(npin), Ein(npin) -character(len=2) :: At1x, At2x, DistUnit, EnerUnit +character(len=2) :: At1x, At2x character(len=4) :: word, Diatom, Diatomx ! For storing character data using gather/scatter DAFILE operations, it ! is imperative that the strings are aligned on integers. @@ -55,6 +55,8 @@ call SpoolInp(LuIn) Atom1 = '' Atom2 = '' +DistUnit = 0 +EnerUnit = 0 ipot = 0 ! Indicator for potential input ngrid = 199 ! Maximum number of grid points Rmin = One @@ -84,8 +86,6 @@ call RdNLst(LuIn,'VibRot') ! Read input data from input file -write(*,*)'Deyan molcas' - ntit1 = 0 skip = .false. input: do @@ -406,12 +406,12 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) - call Get_S(1,DistUnit(1:2),1) + call Get_I1(1,DistUnit) case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) - call Get_S(1,EnerUnit(1:2),1) + call Get_I1(1,EnerUnit) case (tabinp(21)) exit input @@ -589,12 +589,125 @@ if ((ipot == 0) .and. (ncase == 1)) then call Quit_OnUserError() end if -if (ipot /= 0) then - do i=1,nop - Rin(i) = Rin(i) * 1.8897259886 - Ein(i) = Ein(i) * 0.0000045563352812122295 - end do -end if +! Distance units + +select case (DistUnit) + + case (0) + ! 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 (1) + ! 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) * 1.8897259886 + end do + end if + + case (2) + ! 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) * 0.0188973 + 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 + +! Energy units + +select case (EnerUnit) + +case (0) + ! Energy units of hartrees, no need for conversion + write(u6,*) + write(u6,*) 'Energy provided in units of hartrees.' + write(u6,*) 'No conversion.' + +case (1) + ! Energy units of eV, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of electron volts.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.0367493 + end do + end if + +case (2) + ! Energy units of kcal/mol, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of kcal/mol.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.00159360264 + end do + end if + +case (3) + ! Energy units of kJ/mol, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of kJ/mol.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.00038087983 + end do + end if + +case (4) + ! Energy units of cm^(-1), convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of cm^(-1).' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 0.00000455633 + end do + end if + +case (5) + ! Energy units of MHz, convert to hartrees + write(u6,*) + write(u6,*) 'Distance provided in units of MHz.' + write(u6,*) 'Converting to hartrees.' + + if (ipot /= 0) then + do i=1,nop + Ein(i) = Ein(i) * 1.519829E-10 + 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 ! Print input potential diff --git a/src/vibrot/vibrot_globals.F90 b/src/vibrot/vibrot_globals.F90 index 97bc081864..dbe0f5d993 100644 --- a/src/vibrot/vibrot_globals.F90 +++ b/src/vibrot/vibrot_globals.F90 @@ -18,7 +18,7 @@ private integer(kind=iwp), parameter :: nRot_Max = 200, nobs = 10, npin = 500, npoint = 5000 integer(kind=iwp) :: J1A, J2A, lambda, ispc, iobs, nop, Vibwvs, iadvib, Vibwvs1, Vibwvs2, n0, nvib1, n02, nvib21, J1B, J2B, & - IfPrWf, iscale, iallrot + IfPrWf, iscale, iallrot, DistUnit, EnerUnit 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 @@ -26,6 +26,6 @@ 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 -- GitLab From 5ef2fd0de3c9919727b0d433ff0c0d48de7688ff Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 12 Dec 2022 17:03:18 -0500 Subject: [PATCH 21/42] Updated the example input in the vibrot docs because: - The previous input file did not have a complete potential, which made it impossible for anyone to use, without generating their own FeNi potential (a super inconvenience). This means people can't do a quick test to see what the output of the program is supposed to look like, without doing substantial work. - The new input file contains the new feature which allows people to use units that are more common in spectroscopy (which is where vibrot is most likely to be useful). - The new input also matches test/standard/016.input which is convenient. --- doc/source/users.guide/programs/vibrot.rst | 53 ++++++++++++++-------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 6c75ba8ca3..8d98d4031c 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -392,25 +392,40 @@ Input example :: &VIBROT - RoVibrational spectrum - Title = Vib-Rot spectrum for FeNi - Atoms = 0 Fe 0 Ni - 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 - Observable - Dipole Moment - 1.0 0.102354 - 1.1 0.112898 - [...] - Plot = 1.0 10.0 0.1 + RoVibrational spectrum + Title = H2 1Sg+ + 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 **Comments**: The vibrational-rotation spectrum for :math:`\ce{FeNi}` will be computed using the potential curve given in input. The 10 -- GitLab From 58d492eddef46e5f21c78de7e6648e10270b6f1e Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Sun, 18 Dec 2022 23:55:04 -0500 Subject: [PATCH 22/42] These two things: - added descriptions for the DistUnit and EnerUnit keywords. - Replaced comment about the example input file, with a new one which reflects the new example input file being used. --- doc/source/users.guide/programs/vibrot.rst | 34 +++++++++++++++++----- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 8d98d4031c..23c027e952 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -386,6 +386,28 @@ input: +:kword:`DISTunit` + Unit used for distances in the input potential. The default is `BOHR`. Other + options include `ANGSTOM` and `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 `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. + + .. xmldoc:: + %%Keyword: ENERunit + + Specifies the unit used for energies in the input potential. + + + Input example ............. @@ -427,12 +449,10 @@ Input example Orbital = 0 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 -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. +**Comments**: The vibrational-rotation spectrum for :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`\=0 to 3. The mass for +the most abundant isotope of :math:`\ce{H}` will be used. .. xmldoc:: -- GitLab From 76ed1d18a411fae7f7cbfb7370af8fbb191d7118 Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Mon, 19 Dec 2022 00:02:46 +0100 Subject: [PATCH 23/42] Use constants from constants.F90 --- src/vibrot/vibinp.F90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index bd7a293846..e9b7692279 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -20,8 +20,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, DistUnit, EnerUnit -use Constants, only: Zero, One, Five, UTOAU -use Definitions, only: wp, iwp, u6 +use Constants, only: Zero, One, Five, UTOAU, Angstrom, auToeV, auTokcalmol, auToHz, auTocm, cal_to_J +use Definitions, only: wp, iwp, r8, u6 implicit none integer(kind=iwp), intent(out) :: ncase, ngrid, nvib @@ -607,7 +607,7 @@ select case (DistUnit) if (ipot /= 0) then do i=1,nop - Rin(i) = Rin(i) * 1.8897259886 + Rin(i) = Rin(i) / Angstrom end do end if @@ -619,7 +619,7 @@ select case (DistUnit) if (ipot /= 0) then do i=1,nop - Rin(i) = Rin(i) * 0.0188973 + Rin(i) = Rin(i) * 1.0d-2 / Angstrom end do end if @@ -644,60 +644,60 @@ case (0) case (1) ! Energy units of eV, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of electron volts.' + write(u6,*) 'Energy provided in units of electron volts.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.0367493 + Ein(i) = Ein(i) / auToeV end do end if case (2) ! Energy units of kcal/mol, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of kcal/mol.' + write(u6,*) 'Energy provided in units of kcal/mol.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.00159360264 + Ein(i) = Ein(i) / auTokcalmol end do end if case (3) ! Energy units of kJ/mol, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of kJ/mol.' + write(u6,*) 'Energy provided in units of kJ/mol.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.00038087983 + Ein(i) = Ein(i) / (auTokcalmol * cal_to_J) end do end if case (4) ! Energy units of cm^(-1), convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of cm^(-1).' + write(u6,*) 'Energy provided in units of cm^(-1).' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 0.00000455633 + Ein(i) = Ein(i) / auTocm end do end if case (5) ! Energy units of MHz, convert to hartrees write(u6,*) - write(u6,*) 'Distance provided in units of MHz.' + write(u6,*) 'Energy provided in units of MHz.' write(u6,*) 'Converting to hartrees.' if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 1.519829E-10 + Ein(i) = Ein(i) * 1.0d6 / auToHz end do end if -- GitLab From 623a6f27b04f672fa208c34f3d5314a82e59a18c Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Mon, 19 Dec 2022 00:36:35 +0100 Subject: [PATCH 24/42] Add support for units provided by strings --- src/vibrot/vibinp.F90 | 26 +++++++++++++------------- src/vibrot/vibrot_globals.F90 | 3 ++- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index e9b7692279..327a4467fc 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -55,8 +55,8 @@ call SpoolInp(LuIn) Atom1 = '' Atom2 = '' -DistUnit = 0 -EnerUnit = 0 +DistUnit = '' +EnerUnit = '' ipot = 0 ! Indicator for potential input ngrid = 199 ! Maximum number of grid points Rmin = One @@ -406,12 +406,12 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) - call Get_I1(1,DistUnit) + DistUnit = Line case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) - call Get_I1(1,EnerUnit) + EnerUnit = Line case (tabinp(21)) exit input @@ -593,13 +593,13 @@ end if select case (DistUnit) - case (0) + 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 (1) + case ('ANGSTROM') ! Distance units of Angstroms, convert to Bohr radii write(u6,*) write(u6,*) 'Distance provided in units of Angstroms.' @@ -611,7 +611,7 @@ select case (DistUnit) end do end if - case (2) + case ('PICOMETER') ! Distance units of picometers, convert to Bohr radii write(u6,*) write(u6,*) 'Distance provided in units of picometers.' @@ -635,13 +635,13 @@ end select select case (EnerUnit) -case (0) +case ('', 'HARTREE') ! Energy units of hartrees, no need for conversion write(u6,*) write(u6,*) 'Energy provided in units of hartrees.' write(u6,*) 'No conversion.' -case (1) +case ('EV') ! Energy units of eV, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of electron volts.' @@ -653,7 +653,7 @@ case (1) end do end if -case (2) +case ('KCAL/MOL') ! Energy units of kcal/mol, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of kcal/mol.' @@ -665,7 +665,7 @@ case (2) end do end if -case (3) +case ('KJ/MOL') ! Energy units of kJ/mol, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of kJ/mol.' @@ -677,7 +677,7 @@ case (3) end do end if -case (4) +case ('CM-1') ! Energy units of cm^(-1), convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of cm^(-1).' @@ -689,7 +689,7 @@ case (4) end do end if -case (5) +case ('MHZ') ! Energy units of MHz, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of MHz.' diff --git a/src/vibrot/vibrot_globals.F90 b/src/vibrot/vibrot_globals.F90 index dbe0f5d993..e6009805e4 100644 --- a/src/vibrot/vibrot_globals.F90 +++ b/src/vibrot/vibrot_globals.F90 @@ -18,10 +18,11 @@ private integer(kind=iwp), parameter :: nRot_Max = 200, nobs = 10, npin = 500, npoint = 5000 integer(kind=iwp) :: J1A, J2A, lambda, ispc, iobs, nop, Vibwvs, iadvib, Vibwvs1, Vibwvs2, n0, nvib1, n02, nvib21, J1B, J2B, & - IfPrWf, iscale, iallrot, DistUnit, EnerUnit + IfPrWf, iscale, iallrot 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, & -- GitLab From 8abae7b3497baaea773f6bdef9476ee2d3040044 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 19 Dec 2022 00:09:34 -0500 Subject: [PATCH 25/42] These things: - Corrected typos in test/standard/016.input (e.g. ground state -> excited state, and aomic term symbols -> molecular term symbols) - For the ground state case alraedy in 016.input, I added the default unit specifications - For the ground state case alraedy in 016.input, I made a copy of it and used spectroscopic units instead of atomic units. --- test/standard/016.input | 50 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/test/standard/016.input b/test/standard/016.input index 2049b9e797..c7704c1b0f 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 -- GitLab From fc766deeebaedc499647e8433718640de9a3ab50 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 19 Dec 2022 01:09:30 -0500 Subject: [PATCH 26/42] Updated test/standard: - 016.input now has a check file appended to the end of it, after running "pymolcas verify 016 --generate" on the updated 016.input --- test/standard/016.input | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/standard/016.input b/test/standard/016.input index c7704c1b0f..36fbc30cd8 100644 --- a/test/standard/016.input +++ b/test/standard/016.input @@ -158,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 -- GitLab From 479e415448b22a006dd4d594e4d2ebb5229a16f7 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 09:56:29 -0500 Subject: [PATCH 27/42] Resolved the following vibrot documentaiton issues: - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890430 (typo in the word ANGSTROM in documentation) - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890430 (use of "SINGLE" instead of "CHOICE" in documentation) - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890433 (Added back the Observable and Plot keywords to the example input) --- doc/source/users.guide/programs/vibrot.rst | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 23c027e952..f7ce7d5fe8 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -388,9 +388,10 @@ input: :kword:`DISTunit` Unit used for distances in the input potential. The default is `BOHR`. Other - options include `ANGSTOM` and `PICOMETER`. + options include `ANGSTROM` and `PICOMETER`. The short form `PM` can also be used, + instead of `PICOMETER`. - .. xmldoc:: + .. xmldoc:: %%Keyword: DISTunit Specifies the unit used for distances in the input potential. @@ -401,7 +402,7 @@ input: Unit used for energies in the input potential. The default is `HARTREE`. Other options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. - .. xmldoc:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. @@ -447,12 +448,20 @@ Input example Vibrations = 3 Rotations = 0 3 Orbital = 0 + Observable + Dipole Moment + 1.0 0.102354 + 1.1 0.112898 + [...] + Plot = 1.0 10.0 0.1 Scale **Comments**: The vibrational-rotation spectrum for :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`\=0 to 3. The mass for -the most abundant isotope of :math:`\ce{H}` will be used. +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. .. xmldoc:: -- GitLab From 4fdecbcfe81ace9df379bbee011b76da845216e8 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 10:25:43 -0500 Subject: [PATCH 28/42] Fixed the following issues in vibinp (withiin vibrot): - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890442 (defaults are now set to BOHR and HARTREE instead of '' and '') - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890451 call UpCase(Line) was added for DistUnit and EnerUnit - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890454 (no longer explicitly accepting '' as input for DistUnit and EnerUnit) - https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1215890458 (added ability to allow 'PM" for picometers and 'MEGAHERZ' for MHz) --- src/vibrot/vibinp.F90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 327a4467fc..15b9181db2 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -55,8 +55,8 @@ call SpoolInp(LuIn) Atom1 = '' Atom2 = '' -DistUnit = '' -EnerUnit = '' +DistUnit = 'BOHR' +EnerUnit = 'HARTREE' ipot = 0 ! Indicator for potential input ngrid = 199 ! Maximum number of grid points Rmin = One @@ -406,11 +406,13 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) + call Upcase(line) DistUnit = Line case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) + call Upcase(line) EnerUnit = Line case (tabinp(21)) @@ -593,7 +595,7 @@ end if select case (DistUnit) - case ('', 'BOHR') + case ('BOHR') ! Distance units of Bohr radii, no need for conversion write(u6,*) write(u6,*) 'Distance provided in units of Bohr radii.' @@ -611,7 +613,7 @@ select case (DistUnit) end do end if - case ('PICOMETER') + case ('PICOMETER','PM') ! Distance units of picometers, convert to Bohr radii write(u6,*) write(u6,*) 'Distance provided in units of picometers.' @@ -635,7 +637,7 @@ end select select case (EnerUnit) -case ('', 'HARTREE') +case ('HARTREE') ! Energy units of hartrees, no need for conversion write(u6,*) write(u6,*) 'Energy provided in units of hartrees.' @@ -689,7 +691,7 @@ case ('CM-1') end do end if -case ('MHZ') +case ('MHZ','MEGAHERZ') ! Energy units of MHz, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of MHz.' -- GitLab From 94ed4b3b08968b28daef877492c5e9482ae8d43e Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 10:29:28 -0500 Subject: [PATCH 29/42] Fixed typo for Megahertz in vibinp.F90 within vibrot --- src/vibrot/vibinp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 15b9181db2..ba62d739d6 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -691,7 +691,7 @@ case ('CM-1') end do end if -case ('MHZ','MEGAHERZ') +case ('MHZ','MEGAHERTZ') ! Energy units of MHz, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in units of MHz.' -- GitLab From e49e36ee88859d85143afa5d834f68d982745ab7 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 10:31:52 -0500 Subject: [PATCH 30/42] Just made capitalization of L in Line consistent for: DistUnit and EnerUnit parts of vibinp.F90 within vibrot --- src/vibrot/vibinp.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index ba62d739d6..2338bfe35a 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -406,13 +406,13 @@ input: do case (tabinp(19)) ! Distance units Line = Get_Ln(LuIn) - call Upcase(line) + call Upcase(Line) DistUnit = Line case (tabinp(20)) ! Energy units Line = Get_Ln(LuIn) - call Upcase(line) + call Upcase(Line) EnerUnit = Line case (tabinp(21)) -- GitLab From 0600cfebdee3aeb22fecea38fa5d46049ab43617 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 20 Dec 2022 11:04:05 -0500 Subject: [PATCH 31/42] Added to vibrot documentation: - full dipole moment function for H2 --- doc/source/users.guide/programs/vibrot.rst | 27 ++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index f7ce7d5fe8..ea6834f312 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -449,10 +449,29 @@ Input example Rotations = 0 3 Orbital = 0 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 -- GitLab From 1762b8e93e6283f189152f971fe9705fd7f7ec4a Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Fri, 23 Dec 2022 13:19:13 +0100 Subject: [PATCH 32/42] First pass at converting distance units for observables --- src/vibrot/vibinp.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 2338bfe35a..151e6e2e5e 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -611,6 +611,9 @@ select case (DistUnit) do i=1,nop Rin(i) = Rin(i) / Angstrom end do + do i=1,iobs + RinO(1,i) = RinO(1,i) / Angstrom + end do end if case ('PICOMETER','PM') @@ -623,6 +626,9 @@ select case (DistUnit) do i=1,nop Rin(i) = Rin(i) * 1.0d-2 / Angstrom end do + do i=1,iobs + RinO(1,i) = RinO(1,i) * 1.0d-2 / Angstrom + end do end if case default -- GitLab From cea52301f6a7ac43b480475ad0bbc11f65d4d66b Mon Sep 17 00:00:00 2001 From: Deyan Mihaylov Date: Fri, 23 Dec 2022 14:36:19 +0100 Subject: [PATCH 33/42] Fix decimal notation --- src/vibrot/vibinp.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 151e6e2e5e..a7dd37e85d 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -624,10 +624,10 @@ select case (DistUnit) if (ipot /= 0) then do i=1,nop - Rin(i) = Rin(i) * 1.0d-2 / Angstrom + Rin(i) = Rin(i) * 1.0e-2_wp / Angstrom end do do i=1,iobs - RinO(1,i) = RinO(1,i) * 1.0d-2 / Angstrom + RinO(1,i) = RinO(1,i) * 1.0e-2_wp / Angstrom end do end if @@ -705,7 +705,7 @@ case ('MHZ','MEGAHERTZ') if (ipot /= 0) then do i=1,nop - Ein(i) = Ein(i) * 1.0d6 / auToHz + Ein(i) = Ein(i) * 1.0e6_wp / auToHz end do end if -- GitLab From 9b149c8b06fe978205f401a97259e4c123dce8f7 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Fri, 23 Dec 2022 15:03:36 -0500 Subject: [PATCH 34/42] Added LIST to vibrot.rst for: DistUnit and EnerUnit keywords as explained here: https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1220348740 --- doc/source/users.guide/programs/vibrot.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index ea6834f312..cd34a604fb 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -391,7 +391,7 @@ input: options include `ANGSTROM` and `PICOMETER`. The short form `PM` can also be used, instead of `PICOMETER`. - .. xmldoc:: + .. xmldoc:: %%Keyword: DISTunit Specifies the unit used for distances in the input potential. @@ -402,7 +402,7 @@ input: Unit used for energies in the input potential. The default is `HARTREE`. Other options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. - .. xmldoc:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. -- GitLab From a124b57025bb548fbd572e76b510fdf8f7cbb75f Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 27 Dec 2022 12:19:55 -0500 Subject: [PATCH 35/42] Switched vibrot.rst example potential: from ground electronic state to ^1Pi_u --- doc/source/users.guide/programs/vibrot.rst | 128 ++++++++++----------- 1 file changed, 64 insertions(+), 64 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index cd34a604fb..0f15fcfad1 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -414,71 +414,71 @@ Input example :: - &VIBROT - RoVibrational spectrum - Title = H2 1Sg+ - 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 - Observable - 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{H2}` -will be computed using the potential curve given in the input. The 3 +&VIBROT + RoVibrational spectrum + Title = H2 (^1 Pi_u) + Atoms = 0 H 0 H + Potential + 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 + 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 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`\=0 to 3. The mass for +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. -- GitLab From 2f324a7cbab9145970681705a24f22eba2e0ec00 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Tue, 27 Dec 2022 13:16:29 -0500 Subject: [PATCH 36/42] For unit conversion of observables in Vibrot: - The number of observables is iobs, and for each of those there's npobs radial points. Deyan had been looping over iobs and fixing the radial-point index at 1, instead of looping over both parameters. --- src/vibrot/vibinp.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index a7dd37e85d..b430f3400e 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -611,8 +611,10 @@ select case (DistUnit) do i=1,nop Rin(i) = Rin(i) / Angstrom end do - do i=1,iobs - RinO(1,i) = RinO(1,i) / Angstrom + do i = 1, iobs + do j = 1, npobs(i) + RinO(j,i) = RinO(j,i) / Angstrom + end do end do end if @@ -626,8 +628,10 @@ select case (DistUnit) do i=1,nop Rin(i) = Rin(i) * 1.0e-2_wp / Angstrom end do - do i=1,iobs - RinO(1,i) = RinO(1,i) * 1.0e-2_wp / Angstrom + 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 -- GitLab From b774632896957e6b964e9cce24b130d3326aac4d Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Fri, 30 Dec 2022 21:51:39 -0500 Subject: [PATCH 37/42] These things: As suggested and approved here: https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1222944358, the following lines were implemented: - KIND="CHOICE" LIST="BOHR,ANGSTROM,PICOMETER" - KIND="CHOICE" LIST="HARTREE,EV:ELECTRONVOLTS,KCAL/MOL,KJ/MOL,CM-1,MHZ" Also, ELECTRONVOLTS was added to vibinp as a possible input (so not just EV) Also, I fixed some other aspects of the language used in write statements related to the unit conversion, e.g. electronvolts without a space --- doc/source/users.guide/programs/vibrot.rst | 4 +-- src/vibrot/vibinp.F90 | 34 +++++++++++++++------- 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 0f15fcfad1..dfdc49a404 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -391,7 +391,7 @@ input: options include `ANGSTROM` and `PICOMETER`. The short form `PM` can also be used, instead of `PICOMETER`. - .. xmldoc:: + .. xmldoc:: %%Keyword: DISTunit Specifies the unit used for distances in the input potential. @@ -402,7 +402,7 @@ input: Unit used for energies in the input potential. The default is `HARTREE`. Other options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. - .. xmldoc:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index b430f3400e..9afa928cf4 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -650,14 +650,26 @@ select case (EnerUnit) case ('HARTREE') ! Energy units of hartrees, no need for conversion write(u6,*) - write(u6,*) 'Energy provided in units of hartrees.' + write(u6,*) 'Energy provided in units of hartree.' write(u6,*) 'No conversion.' +case ('ELECTRONVOLTS') + ! 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 ('EV') ! Energy units of eV, convert to hartrees write(u6,*) - write(u6,*) 'Energy provided in units of electron volts.' - write(u6,*) 'Converting to hartrees.' + write(u6,*) 'Energy provided in electronvolts.' + write(u6,*) 'Converting to hartree.' if (ipot /= 0) then do i=1,nop @@ -668,8 +680,8 @@ case ('EV') case ('KCAL/MOL') ! Energy units of kcal/mol, convert to hartrees write(u6,*) - write(u6,*) 'Energy provided in units of kcal/mol.' - write(u6,*) 'Converting to hartrees.' + write(u6,*) 'Energy provided in kcal/mol.' + write(u6,*) 'Converting to hartree.' if (ipot /= 0) then do i=1,nop @@ -680,8 +692,8 @@ case ('KCAL/MOL') case ('KJ/MOL') ! Energy units of kJ/mol, convert to hartrees write(u6,*) - write(u6,*) 'Energy provided in units of kJ/mol.' - write(u6,*) 'Converting to hartrees.' + write(u6,*) 'Energy provided in kJ/mol.' + write(u6,*) 'Converting to hartree.' if (ipot /= 0) then do i=1,nop @@ -692,8 +704,8 @@ case ('KJ/MOL') case ('CM-1') ! Energy units of cm^(-1), convert to hartrees write(u6,*) - write(u6,*) 'Energy provided in units of cm^(-1).' - write(u6,*) 'Converting to hartrees.' + write(u6,*) 'Energy provided in cm^(-1).' + write(u6,*) 'Converting to hartree.' if (ipot /= 0) then do i=1,nop @@ -704,8 +716,8 @@ case ('CM-1') case ('MHZ','MEGAHERTZ') ! Energy units of MHz, convert to hartrees write(u6,*) - write(u6,*) 'Energy provided in units of MHz.' - write(u6,*) 'Converting to hartrees.' + write(u6,*) 'Energy provided in MHz.' + write(u6,*) 'Converting to hartree.' if (ipot /= 0) then do i=1,nop -- GitLab From 3e159d095287291125046eefbe465219f0666e79 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Fri, 30 Dec 2022 21:55:48 -0500 Subject: [PATCH 38/42] Changed electronvolts to electronvolt (singular), also: previous commit wasn't precisely as I decsribed in the commit message, I'd done ELECTRONVOLTS instead of EV:ELECTRONVOLTS (per request) in GitLab discussion. --- doc/source/users.guide/programs/vibrot.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index dfdc49a404..3912cb2cab 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -402,7 +402,7 @@ input: Unit used for energies in the input potential. The default is `HARTREE`. Other options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. - .. xmldoc:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. -- GitLab From c68f288e8a91b128719912fc6442a4d22458d5a5 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Fri, 30 Dec 2022 22:00:02 -0500 Subject: [PATCH 39/42] ELECTRONVOLT and EV lines are: now combined in the vibinp code (and both singular rather than plural) --- src/vibrot/vibinp.F90 | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 9afa928cf4..5d0fbe1ee9 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -653,19 +653,7 @@ case ('HARTREE') write(u6,*) 'Energy provided in units of hartree.' write(u6,*) 'No conversion.' -case ('ELECTRONVOLTS') - ! 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 ('EV') +case ('EV','ELECTRONVOLT') ! Energy units of eV, convert to hartrees write(u6,*) write(u6,*) 'Energy provided in electronvolts.' -- GitLab From 8f22250d7e1d64b8ec0f3813b02af9748b8591be Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Sat, 31 Dec 2022 14:46:19 -0500 Subject: [PATCH 40/42] In vibinp.F90: moved tings around for the unit conversion code, in order to suppress warnings of the form: "CHARACTER expression will be truncated in assignment (9/180) at (1) [-Werror=character-truncation]" due to "DistUnit = Line" and "EnerUnit = Line" --- src/vibrot/vibinp.F90 | 292 ++++++++++++++++++++++++++---------------- 1 file changed, 182 insertions(+), 110 deletions(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 5d0fbe1ee9..18aace690a 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -407,13 +407,139 @@ input: do ! Distance units Line = Get_Ln(LuIn) call Upcase(Line) - DistUnit = 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 +! 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 @@ -591,135 +717,81 @@ if ((ipot == 0) .and. (ncase == 1)) then call Quit_OnUserError() end if -! Distance units -select case (DistUnit) - 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 -! Energy units -select case (EnerUnit) -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 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ! Print input potential -- GitLab From 219f36f21ea3566557ab2c8af7bfb992bca0ab80 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Sat, 31 Dec 2022 15:04:53 -0500 Subject: [PATCH 41/42] These things: - cleaned up vibinp.F90 a bit (whitespace leftover from previous edit) - responding to this: https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/585#note_1224710561 I've made Megahertz show up by default in the documentation, and MHz is now the "hidden" variable (although I've made it clear in the documentation that MHZ and EV can be used instead of MEGAHERTZ and ELECTRONVOLT --- doc/source/users.guide/programs/vibrot.rst | 6 +- src/vibrot/vibinp.F90 | 76 ---------------------- 2 files changed, 4 insertions(+), 78 deletions(-) diff --git a/doc/source/users.guide/programs/vibrot.rst b/doc/source/users.guide/programs/vibrot.rst index 3912cb2cab..206c255c40 100644 --- a/doc/source/users.guide/programs/vibrot.rst +++ b/doc/source/users.guide/programs/vibrot.rst @@ -400,9 +400,11 @@ input: :kword:`ENERunit` Unit used for energies in the input potential. The default is `HARTREE`. Other - options include `EV` (electron Volts), `KCAL/MOL`, `KJ/MOL`, `CM-1`, and `MHZ`. + 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:: + .. xmldoc:: %%Keyword: ENERunit Specifies the unit used for energies in the input potential. diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 18aace690a..8c7aaefc25 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -717,82 +717,6 @@ if ((ipot == 0) .and. (ncase == 1)) then call Quit_OnUserError() end if - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Print input potential if (ipot /= 0) then -- GitLab From 32a665e9144bafa9a24fe581017ebe0c5a95af87 Mon Sep 17 00:00:00 2001 From: Nike Dattani Date: Mon, 9 Jan 2023 00:57:06 -0500 Subject: [PATCH 42/42] No longer explicitly import r8: This should fix the NAG compiler issue: https://gitlab.com/Molcas/OpenMolcas/-/jobs/3542374039#L5308 I had pointed out the incorrect r8 here: https://gitlab.com/Molcas/OpenMolcas/-/merge_requests/580#note_1161174037 but it had actually been fixed 5 days earlier here: https://gitlab.com/Molcas/OpenMolcas/-/commit/2e4f878f566842978efc4185fadea9687be6db01 Somehow Deyan still had the r8 explicitly imported, but the place where it was accidentally used had been changed to wp. Strange. --- src/vibrot/vibinp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/vibrot/vibinp.F90 b/src/vibrot/vibinp.F90 index 8c7aaefc25..477ed9c210 100644 --- a/src/vibrot/vibinp.F90 +++ b/src/vibrot/vibinp.F90 @@ -21,7 +21,7 @@ use Vibrot_globals, only: Atom1, Atom2, dRo, EoutO, iad12, iad13, iadvib, iallro J2A, J2B, lambda, n0, n02, nop, npin, npobs, npoint, nRot_Max, nvib1, nvib21, Obsin, R0o, R1o, RinO, & 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, r8, u6 +use Definitions, only: wp, iwp, u6 implicit none integer(kind=iwp), intent(out) :: ncase, ngrid, nvib -- GitLab