two_BOBS

two_BOBS

diatom.f90

@@ -283,6 +283,7 @@ module diatom_module
  logical :: intensity = .false.
  logical :: frequency = .false.
  logical :: matelem = .false.
+ logical :: raman = .false.
  !
  end type actionT
  !
@@ -3058,7 +3059,7 @@ subroutine ReadInput
  !
  select case(w)
  !
- case('NONE','ABSORPTION','EMISSION','TM','DIPOLE-TM','RAMAN','POLARIZABILITY-TM','PARTFUNC')
+ case('NONE','ABSORPTION','EMISSION','TM','DIPOLE-TM','PARTFUNC')
  !
  intensity%action = trim(w)
  !
@@ -3084,6 +3085,10 @@ subroutine ReadInput
  !
  intensity%matelem = .true.
  !
+ case('RAMAN','POLARIZABILITY')
+ !
+ action%raman = .true.
+ !
  case('OVERLAP')
  !
  intensity%overlap = .true.
@@ -5763,6 +5768,8 @@ subroutine duo_j0(iverbose_,J_list_,enerout,quantaout,nenerout)
  stop
  endif
  !
+ vibTmat = 0 
+ !
  do jgrid = igrid,ngrid
  do kgrid=1,ngrid
  vibTmat(igrid,jgrid) = vibTmat(igrid,jgrid) + (12._rk)*(hstep**2)*(LobWeights(kgrid))*& 

dipole.f90

@@ -445,7 +445,7 @@ subroutine dm_intensity(Jval,iverbose)
  ! New RICHMOL format - one file for x,y,z
  !
  filename = &
- "matelem_MU_"//"_j"//trim(adjustl(char_jI))//"_j"//trim(adjustl(char_jF))//"_"//trim(intensity%linelist_file)//".rchm"
+ "matelem_MU"//"_j"//trim(adjustl(char_jI))//"_j"//trim(adjustl(char_jF))//"_"//trim(intensity%linelist_file)//".rchm"
  !
  call IOstart(trim(filename),richunit(indI,indF))
  open(unit = richunit(indI,indF), action = 'write',status='replace' , file = filename)

duo.f90

@@ -4,6 +4,7 @@ program duo
  use accuracy
  use refinement
  use dipole
+ use polarizability
 ! use compilation_details, only: write_compilation_details
  use header_info, only: write_logo
  use diatom_module,only : duo_j0,verbose,job,readinput,map_fields_onto_grid,action !, check_and_set_atomic_data
@@ -62,6 +63,14 @@ function isatty(fd) bind(c)
  !
  call map_fields_onto_grid(verbose)
  !
+ if (action%raman) then 
+ call pol_tranint
+ !
+ write(out, '(a)') '--End--'
+ stop
+ !
+ endif 
+ !
  !call define_jlist
  call dm_tranint
  !

functions.f90

@@ -157,6 +157,10 @@ subroutine define_analytical_field(ftype,fanalytical_field)
  !
  fanalytical_field => poten_two_coupled_EMOs
  !
+ case("TWO_COUPLED_BOBS")
+ !
+ fanalytical_field => poten_two_coupled_BOBs
+ !
  case("COUPLED_EMO_REPULSIVE")
  !
  fanalytical_field => poten_two_coupled_EMO_repulsive
@@ -1372,7 +1376,7 @@ function poten_two_coupled_EMOs(r,parameters) result(f)
  f2 = poten_EMO(r,parameters(nparams1+1:nparams1+nparams2))
  a = poten_cosh_polynom(r,parameters(nparams1+nparams2+1:nparams1+nparams2+nparams3))
  !
- discr = f1**2-2.0_rk*f1*f2+f2**2+4*a**2
+ discr = f1**2-2.0_rk*f1*f2+f2**2+4.0_rk*a**2
  !
  if (discr<-small_) then
  write(out,"('poten_two_coupled_EMOs: discriminant is negative')")
@@ -1386,6 +1390,38 @@ function poten_two_coupled_EMOs(r,parameters) result(f)
  !
  end function poten_two_coupled_EMOs
 
+
+ function poten_two_coupled_BOBs(r,parameters) result(f)
+ !
+ real(rk),intent(in) :: r ! geometry (Ang)
+ real(rk),intent(in) :: parameters(:) ! potential parameters
+ real(rk) :: f1,f2,f,f_switch,t_0,r_s,a_s
+ integer(ik) :: nparams1,nparams2,nparams3,icomponent,n
+ !
+ nparams1 = parameters(4)+6
+ nparams2 = parameters(nparams1+4)+6
+ nparams3 = 2
+ icomponent = parameters(nparams1+nparams2+nparams3+1)
+ !
+ N = nparams1+nparams2
+ !
+ f1 = poten_BOBLeRoy(r,parameters(1:nparams1))
+ f2 = poten_BOBLeRoy(r,parameters(nparams1+1:nparams1+nparams2))
+ !
+ r_s = parameters(N+1)
+ a_s = parameters(N+2)
+ !
+ f_switch = ( 1.0_rk+tanh(a_s*(r-r_s)) )*0.5_rk
+ !
+ f = 0 
+ !
+ if (icomponent==1) then 
+ f = f_switch*f2+f1*(1.0_rk-f_switch)
+ else
+ f = f_switch*f1+f2*(1.0_rk-f_switch)
+ endif
+ !
+ end function poten_two_coupled_BOBs
  !
  ! Morse/Long-Range, see Le Roy manuals
  !