#include "precompilerdefinitions" program lineshape !!{!src/lineshape/manual.md!} use konstanter, only: flyt use type_crystalstructure, only: lo_crystalstructure use type_forceconstant_secondorder, only: lo_forceconstant_secondorder use type_forceconstant_thirdorder, only: lo_forceconstant_thirdorder use type_forceconstant_fourthorder, only: lo_forceconstant_fourthorder use type_qpointmesh, only: lo_qpoint_mesh,lo_generate_qmesh,lo_read_qmesh_from_file,lo_get_small_group_of_qpoint use type_phonon_dispersions, only: lo_phonon_dispersions use type_phonon_dos, only: lo_phonon_dos use type_phonon_bandstructure, only: lo_phonon_bandstructure use helpers, only: open_file,lo_mpi_helper,walltime,lo_chop !use scatteringrates, only: lo_listofscatteringrates use phonondamping use options use io implicit none type(lo_opts) :: opts type(lo_mpi_helper) :: mw type(lo_crystalstructure) :: uc type(lo_forceconstant_secondorder) :: fc type(lo_forceconstant_thirdorder) :: fct type(lo_forceconstant_fourthorder) :: fcf type(lo_phonon_dispersions) :: dr,drd type(lo_phonon_dos) :: pd class(lo_qpoint_mesh), allocatable :: qp,qpd real(flyt) :: timer_init,timer_total ! Init MPI! call mw%init() ! Start timers timer_total=walltime() timer_init=walltime() ! some options call opts%parse() ! set up all possible meshes init: block ! only be verbose on the first rank if ( .not. mw%talk ) opts%verbosity=0 ! Read structure call uc%readfromfile('infile.ucposcar') call uc%classify('wedge',timereversal=.true.) if ( mw%talk ) write(*,*) '... using ',tochar(mw%n),' MPI ranks' if ( mw%talk ) write(*,*) '... read unitcell poscar' if ( opts%readiso ) then if ( mw%talk ) write(*,*) '... reading isotope distribution from file' call uc%readisotopefromfile() endif ! Read forceconstants call fc%readfromfile(uc,'infile.forceconstant') if ( mw%talk ) write(*,*) '... read second order forceconstant' if ( opts%thirdorder ) then call fct%readfromfile(uc,'infile.forceconstant_thirdorder') if ( mw%talk ) write(*,*) '... read third order forceconstant' endif if ( opts%fourthorder ) then call fcf%readfromfile(uc,'infile.forceconstant_fourthorder') if ( mw%talk ) write(*,*) '... read fourth order forceconstant' endif ! Get a q-mesh for the integrations if ( opts%readqmesh ) then if ( mw%talk ) write(*,*) '... getting q-mesh from file' call lo_read_qmesh_from_file(qp,uc,'infile.qgrid.hdf5') else if ( mw%talk ) write(*,*) '... generating q-mesh' select case(opts%meshtype) case(1) call lo_generate_qmesh(qp,uc,opts%qgrid,'monkhorst',verbosity=opts%verbosity,timereversal=opts%timereversal) case(2) call lo_generate_qmesh(qp,uc,opts%qgrid,'fft',verbosity=opts%verbosity,timereversal=opts%timereversal) case(3) call lo_generate_qmesh(qp,uc,opts%qgrid,'wedge',verbosity=opts%verbosity,timereversal=opts%timereversal) end select endif ! Dispersions for everyone! call dr%generate(qp,fc,uc,timereversal=opts%timereversal,verbosity=opts%verbosity,mpi_communicator=mw%comm) if ( mw%talk ) write(*,*) '... got the full dispersion relations' timer_init=walltime()-timer_init end block init ! Now everything is read and set, time to calculate something. Three main modes: single point, ! bandstructure or dos. Only do one of them, gets annoying and confusing otherwise, I think. if ( opts%oneqpoint ) then onepoint: block type(lo_phonon_dispersions_qpoint) :: ompoint type(lo_qpoint) :: qpoint type(lo_phonon_selfenergy) :: se real(flyt), dimension(3) :: v0 ! Do the rest on a single rank. No point in using MPI at all. if ( mw%talk ) then write(*,*) '' write(*,*) 'Caltulating lineshapes for a single q-point' opts%verbosity=opts%verbosity+1 endif ! coordinates for this single point if ( trim(opts%highsymmetrypoint) .ne. 'none' ) then v0=uc%coordinate_from_high_symmetry_point_label(opts%highsymmetrypoint) else v0=uc%fractional_to_cartesian(opts%qpoint,reciprocal=.true.) endif qpoint%v=lo_chop(v0,lo_sqtol) qpoint%w=lo_chop(v0-uc%bz%gshift(v0+lo_degenvector),lo_sqtol) call lo_get_small_group_of_qpoint(qpoint,uc) ! harmonic properties at this point call ompoint%generate(fc,uc,qpoint) ! get the self-energy call se%generate(qpoint,ompoint,uc,fc,fct,fcf,qp,dr,opts,mw) ! dump it to file if ( mw%talk ) then call write_lineshape_to_hdf5(se,opts%enhet,opts%temperature,'outfile.lineshape.hdf5') endif end block onepoint endif if ( opts%qpointpath ) then qppath: block type(lo_phonon_bandstructure) :: bs type(lo_phonon_selfenergy) :: se integer :: npts,i,j if ( mw%talk ) then write(*,*) '' write(*,*) 'Calculating lineshapes on a path in the BZ.' endif ! Generate the path, and the harmonic things on said path npts=0 do i=1,opts%nq_on_path npts=npts+opts%stride if ( npts .ge. opts%nq_on_path-opts%stride ) exit enddo npts=npts+1 call bs%generate(uc,fc,timereversal=opts%timereversal,npts=npts,readpathfromfile=opts%readpathfromfile,verbosity=opts%verbosity,mpi_communicator=mw%comm) ! Calculate the actual thingy. Quite long, so moved to it's own routine. call spectralfunction_along_path(bs,uc,fc,fct,fcf,qp,dr,opts,mw) ! Dump this if ( mw%talk ) then call write_lineshapes_along_path_to_hdf5(bs,opts%enhet,'outfile.sqe.hdf5') call bs%write_to_hdf5(uc,opts%enhet,'outfile.dispersion_relations.hdf5') endif end block qppath endif if ( opts%phonondos ) then if ( mw%talk ) write(*,*) '' if ( mw%talk ) write(*,*) 'Calculating the phonon dos' select case(opts%meshtype) case(1) call lo_generate_qmesh(qpd,uc,opts%qgrid_dos,'monkhorst',verbosity=opts%verbosity,timereversal=opts%timereversal) case(2) call lo_generate_qmesh(qpd,uc,opts%qgrid_dos,'fft',verbosity=opts%verbosity,timereversal=opts%timereversal) case(3) call lo_generate_qmesh(qpd,uc,opts%qgrid_dos,'wedge',verbosity=opts%verbosity,timereversal=opts%timereversal) end select ! and dispersions on this grid call drd%generate(qpd,fc,uc,timereversal=opts%timereversal,verbosity=opts%verbosity,mpi_communicator=mw%comm) ! Get the actual dos call get_intensity_as_dos(pd,qpd,drd,uc,fc,fct,fcf,qp,dr,opts,mw) ! Dump the dos to file if ( mw%talk ) then call pd%write_to_file(uc,opts%enhet,'outfile.lineshape_phonon_dos') endif endif ! All done, print timings if ( mw%talk ) then timer_total=walltime()-timer_total write(*,*) ' ' write(*,*) ' Timings:' write(*,*) ' initialization:',timer_init write(*,*) ' total:',timer_total endif ! Kill MPI call mw%destroy() end program