anchor.diff:                 ener_vib = energies_vib(ivib)
                  ener_rot = energy-ener_vib
                  !
  -               if (ener_rot<small_) then
  +               if (ener_rot<-small_) then
                    !
                    ! reassign vibrational energy to the lowest  
                    !
  
positive.code:                   !
                    ! let's assume that negative rotational energies can be made positive without breaking the physics
                    ener_rot = abs(energyi-ener_vib)
                  else
                    ivib = ivib_state(ilevelf)
                    ener_vib = energies_vib(ivib)
                    ener_rot = abs(energyf-ener_vib)
                  endif
               endif
               !
  
negative.code: 
                 elseif(vibpopulation_do) then
                    ! apply vibrational population
                    !read(quantum_numbers(QN%dens_col-4,ilevelf),*) ndensity
                    !
                    abscoef=emcoef*acoef*gtot(ilevelf)*exp(-c2/temp*ener_rot)*ndensity*&
                            tranfreq/(partfunc)
  