plot-mcdspectrum.f90

program plotspec

  ! copyright: Jochen Autschbach, SUNY Buffalo
  
  ! output graph data with MCD data generated by Molcas and a
  ! utility program, for plotting broadened spectra
  
  ! each excitation is broadened with a normalized Gaussian
  ! function

  ! input options with defaults are as follows:

  !  sigma = 0 ! use calculated (empirical) linewidth, otherwise
  !            ! sigma is given explicitly if nonzero 
  !  sharpen = 1.0d0 ! "sharpen" the linewidth by a factor 1/sharpen
  !                  ! (only for sigma=0)
  !  npoints = 200 ! no. of points to generate for the plot
  !  nexcit = 0 ! default no. of excitations to plot (determined by input
  !             ! or by symmetry info if zero
  !  border = 3d0 ! add this many sigmas at the plot boundaries
  !  lorentz = .f. ! use lorentz curve for lineshape instead of gauss


  implicit none

  double precision, dimension (:,:), allocatable :: results
  integer, dimension(100) :: nsyme, ndegen
  data nsyme/ 100*0/
  data ndegen/ 100*1/
  double precision :: tmp(9), prefac
  double precision, parameter :: one=1d0, zero=0d0, two=2d0, four=4d0, &
    explim = 500, &
    pi=3.14159265359d0, &
    au2cm=2.194746313702d5, &
    boltz=3.1668105d-6, &
    au2ev = 27.21138602d0, &
    xmcd = 0.001480279d0, &
    t2gauss = 10000d0
  
  logical invert, magnify, waveno, theta
  
  double precision :: window(2), sigma, sharpen, eshift, border, temp
  data window /0d0,0d0/
  integer :: npoints, nexcit, magnification
  character(1) :: term

  ! the namelist plot contains the input options. see above for defaults
  
  namelist /plot/ npoints, sigma, nsyme, nexcit, sharpen, eshift, &
     &            border, ndegen, invert, window,   &
     &            magnification, waveno, term, temp, theta

  ! Here is an example for some typical settings:
  ! &plot nsyme(1)=20,35, ndegen(1)=1,1, sigma=0.13, sharpen=1, npoints=300, 
  !  nexcit=50, invert=.t., waveno=.t. /

  double precision :: deltae, emax, emin, eexcit, fval, rval, t2de, &
    small, multiplier

  integer :: i, j, k, np, ll, isym, nsym, rinvert

  ! each excitation gets assigned a normalized Gauss curve, centered at x0.
  ! sigma is the full width parameter in the formula

  double precision:: gauss, dgauss, x, x0, sigma1, expo

  expo(x,x0,sigma) = (x-x0)**2 /(sigma**2)
  gauss(x,x0,sigma) = 1/(sqrt(pi)*sigma)*exp(-(x-x0)**2 /(sigma**2))
  dgauss(x,x0,sigma) = 2*(x0-x)/(sqrt(pi)*sigma**3)*exp(-(x-x0)**2 /(sigma**2))
 

  ! empirical linewidth from Brown et al., J. Chem. Soc (A) 1971, 751
  ! in wavenumbers

  sigma1(x0) = 1/au2cm * 7.5d0 * sqrt(x0 * au2cm)
  logical lcalsig, deriv

  ! =====================================================================

  write (6,'(/a/)') ' P L O T - M C D S P E C T R U M'

  open (7,file='impulses.dat',status='unknown')
  open (8,file='graph.dat',status='unknown')

  ! plot defaults. some can be changed by namelist input
  eshift = 0.d0 ! shift energies in the spectrum (not used)
  sigma = 0.0d0 ! use calculated (empirical) linewidth
  lcalsig = .true.
  sharpen = 1.0d0 ! "sharpen" the linewidth by a factor 1/sharpen
  npoints = 200
  nexcit = 0 ! default no. of excitations to plot (determined by input
             ! or by symmetry info (ALL)
  border = 3d0 ! add that many sigmas at the plot boundaries
  invert = .false. ! invert the CD spectrum
  magnify = .false. ! magnify a certain energy window
  magnification = 10 ! default magnification
  waveno = .true. ! energies and broadening in wavenumbers
  term = ''  ! which MCD term, MUST be supplied
  temp = 298d0 ! temperature
  theta = .false. ! if true, give [theta] in units per Gauss

  t2de = pi / (4500d0*log(10d0) ) ! [theta] to Delta epsilon conv

  ! how many data points do we have? no. of plot points, sigma?
  ! read the namelist "plot" from input file:
  read (5,plot)

  if (temp.le.zero) stop 'Temp set to zero or less. Aborting.'

  if (term.eq.'A') then
    prefac = -xmcd
    deriv = .true.
  else if (term.eq.'B') then
    prefac = xmcd
    deriv = .false.
  else if (term.eq.'C') then
    prefac = xmcd / (boltz * temp)
    deriv = .false.
  else
    stop 'MCD term spec is not A, B, or C. Aborting.'
  end if

  write (6,*) term//'-term spectrum'

  if (theta) then
    multiplier = one
  else
    multiplier = t2de * t2gauss
  end if
    
  write (6,*) multiplier
  
  if (sigma.ne.0d0) then
    lcalsig = .false.
    if (sigma.le.0d0) then
      write (0,*) 'linewidth < 0 requested. aborting ...'
      close (7,status='delete')
      close (8,status='delete')
      stop 'error termination'
    end if
    if (.not.waveno) then
      write (6,'(1X,a,f10.3,a)') 'using fixed width of ',sigma,' au'
    else
      write (6,'(1X,a,f10.3,a)') 'using fixed width of ',sigma,' cm^-1'
      sigma = sigma / au2cm
    end if
  else
    write (6,'(1X,a)') 'using empirical formula for linewidth' 
  end if

  if (invert) then
    rinvert = -one
  else
    rinvert = one
  end if

  ! magnify certain portion of the plot
  if (window(2).gt.zero) then
    if (window(1).gt.window(2)) stop 'window in error'
    magnify = .true.
  end if

  ! how many symmetries ?
  nsym = 0
  i = 0
  10 continue
  i = i+1
  if (nsyme(i).gt.0) then
    nsym = nsym + 1
    goto 10
  endif
  if (nsym.le.0) then
    write (0,*) 'error nsyme: no entries found. aborting ...'
    close (7,status='delete')
    close (8,status='delete')
    stop 'error termination'
  endif

  np = 0
  do i=1,nsym
    np = np + nsyme(i)
  enddo
  if (np.le.0) then
    write (0,*) 'error np: nothing to plot. aborting ...'
    close (7,status='delete')
    close (8,status='delete')
    stop 'error termination'
  endif

  if (nexcit.eq.0) then ! no input value for nexcit, set it to
    nexcit = np         ! its default = plot all excitations
  else ! check that nexcit is not too large
    if (nexcit.gt.np) then
      write (0,*) 'error: nexcit is greater than no. of energies. aborting'
      close (7,status='delete')
      close (8,status='delete')
      stop 'error termination'
    endif
  endif

  write (6,'(1X,a,i4,a,i4,a)') 'using ',nexcit,' out of ',np, &
     &' excitations for the graph data'

  allocate (results(np,9))
  results = 0

  write (7,*) '# E/cm^-1, MCD-term (au), MCD-term (D**2)'

  ll = 1
  do isym=1,nsym
 !   write (0,*) 'symmetry',isym,nsyme(isym)
    if (isym.eq.1) then
       read (5,'(//)') 
       
    else
       read (5,'(/)') 
       
    end if
    
    ! read the energies and C terms
    do i=ll,ll+nsyme(isym)-1
      read (5,*) results(i,1), results(i,6), results(i,7)
      !          write (6,*) ndum, results(i,1), results(i,2),rdum,results(i,6)
     
      ! write(6,*) results(i,1)*au2cm,results(i,7)
    enddo
        
    ! write "impulses" data
     
    write (7,*) '# Symmetry ',isym
    do i=ll,ll+nsyme(isym)-1

      ! multiply intensities with degeneracy of irrep
      results(i,2) = results(i,2) * float(ndegen(isym))
      results(i,6) = results(i,6) * float(ndegen(isym)) * rinvert
      results(i,7) = results(i,7) * float(ndegen(isym)) * rinvert

      !                    energy       MCD terms in au and D**2
      write (7,'(3e25.15)') results(i,1),results(i,6),results(i,7)
    enddo
    write (7,*) ; write (7,*)

    ll = ll+ nsyme(isym)
  enddo ! isym

  ! sort data in ascending energy (sum up all symmetries)
  ! (probably not a fast algorithm, but it seems to work...)

  do k=2,np
    ll = k-1
    small = results(ll,1)
    do i=k,np
      if (results(i,1).lt. small) then
!        write (6,*) k,i,small, results(i,1)
        tmp(:) = results(ll,:)
        results(ll,:) = results(i,:)
        results(i,:) = tmp(:)
        small = results(ll,1)
      endif
    enddo
  enddo

  do i=2,np
    if (results(i,1).lt.results(i-1,1)) then 
      write (0,*) results(1:np,1)
      write (0,*) 'sorting error, aborting...'
      close (7,status='delete')
      close (8,status='delete')
      stop 'error termination'
    endif
  enddo

  ! convert energies to au if we are in wavenumbers
  
  if (waveno) results(:,1) =  results(:,1) / au2cm

  ! determine plot range:
  emin = results(1,1)
  emax = results(nexcit,1) ! it's nexcit here, not np 

  !  write (6,*) emin,emax

  ! add  border*sigma of approx. 0.15 eV at the energy boundaries for plot:
  emin = emin - border * 0.15d0 / au2ev
  emax = emax + border * 0.15d0 / au2ev

  ! number of plot points defines step size deltae:
  deltae = abs(emax-emin)/float(npoints-1)
  
!  write (6,*) emin,emax,deltae

  write (8,*) '# E/cm^-1, zero , MCD. theta is ',theta
  do i=0,npoints
    fval = 0.d0
    x = emin + float(i) * deltae
    do j=1,np

      eexcit = results(j,1)
      if (lcalsig) sigma = sigma1(eexcit)/sharpen ! compute linewidth
      
      x0 = results(j,1) ! excitation energy = center of the curve

      rval = expo(x,x0,sigma)

      ! if (abs(rval).gt.explim) write (6,*) i,j, rval
      
      if (abs(rval).gt.explim) cycle
      
      if(abs(x-x0).lt.3*sigma) then
         if (deriv) then
            fval = fval + prefac *  x * results(j,6)*dgauss(x,x0,sigma) 
         else
            fval = fval + prefac *  x * results(j,6)*gauss(x,x0,sigma)        
         end if
      end if
   enddo

    if (waveno) x = x * au2cm

    fval = fval * multiplier
    write (8,'(3e25.15)') x,fval
 
  enddo

  ! we are done
  deallocate (results)
  close (7,status='keep')
  close (8,status='keep')


  write (6,*) 'Note: use the script "graph-mcd.sh" to plot the data'
  if (theta) write (6,*) 'and USE THE -t OPTION !'
  stop 'normal termination'
  end