#include "precompilerdefinitions" program thermal_conductivity !!{!src/thermal_conductivity/manual.md!} use konstanter, only: flyt,lo_temperaturetol,lo_status,lo_kappa_au_to_SI use helpers, only: walltime,lo_mpi_helper,tochar,lo_linspace,lo_logspace,lo_mean,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,MPI_SUM use type_crystalstructure, only: lo_crystalstructure use type_mdsim, only: lo_mdsim use type_forceconstant_secondorder, only: lo_forceconstant_secondorder use type_forceconstant_thirdorder, only: lo_forceconstant_thirdorder use type_qpointmesh, only: lo_qpoint_mesh,lo_generate_qmesh use type_phonon_dispersions, only: lo_phonon_dispersions use type_phonon_dos, only: lo_phonon_dos use dump_data, only: lo_dump_gnuplot_2d_real ! unique use options, only: lo_opts use scatteringstrengths, only: calculate_scattering_amplitudes use pbe, only: get_kappa,calculate_qs,get_selfconsistent_solution use phononevents, only: lo_threephononevents,lo_find_all_scattering_events use mfp, only: lo_mfp,get_cumulative_plots,write_cumulative_plots implicit none ! Standard, from libolle type(lo_opts) :: opts type(lo_forceconstant_secondorder) :: fc type(lo_forceconstant_thirdorder) :: fct type(lo_phonon_dispersions) :: dr type(lo_phonon_dos) :: pd type(lo_crystalstructure) :: uc class(lo_qpoint_mesh), allocatable :: qp type(lo_mpi_helper) :: mw ! Unique type(lo_threephononevents) :: sc type(lo_mfp) :: mf ! Small stuff real(flyt), dimension(:,:), allocatable :: thermal_cond real(flyt), dimension(:), allocatable :: temperatures ! timers real(flyt) :: timer_init,timer_count,timer_matrixelements,timer_scf real(flyt) :: timer_kappa,timer_qs,timer_cumulative,tt0 #ifdef AGRESSIVE_SANITY real(flyt), parameter :: toMB=1.0_flyt/1024.0_flyt**2 real(flyt), dimension(:), allocatable :: mem_per_rank #endif ! Set up all harmonic properties. That involves reading all the input file, ! creating grids, getting the harmonic properties on those grids. initharmonic: block integer :: i,j ! Start MPI and timers call mw%init() #ifdef AGRESSIVE_SANITY ! make space for the thing that holds memory per rank lo_allocate(mem_per_rank(mw%n)) mem_per_rank=0 #endif ! Get options call opts%parse() if ( mw%r .ne. 0 ) opts%verbosity=0 tt0=walltime() timer_init=walltime() timer_qs=0.0_flyt timer_kappa=0.0_flyt timer_scf=0.0_flyt timer_cumulative=0.0_flyt if ( mw%talk ) write(*,*) '... using ',tochar(mw%n),' MPI ranks' ! There is a bunch of stuff that all ranks need, first the unit cell: call uc%readfromfile('infile.ucposcar',verbosity=opts%verbosity) call uc%classify('wedge',timereversal=opts%timereversal) if ( mw%talk ) write(*,*) '... read unitcell poscar' ! Perhaps non-natural isotope distribution if ( opts%readiso ) then if ( mw%talk ) write(*,*) '... reading isotope distribution from file' call uc%readisotopefromfile() if ( mw%talk ) then do i=1,uc%na do j=1,uc%isotope(i)%n write(*,"(' isotope: ',I2,' concentration: ',F8.5,' mass: ',F12.6)") & j,uc%isotope(i)%conc(j),uc%isotope(i)%mass(j) enddo write(*,"(' element: ',A2,' mean mass: ',F12.6,' mass disorder parameter',F12.9)") & trim(uc%atomic_symbol( uc%species(i) )),uc%isotope(i)%mean_mass,& uc%isotope(i)%disorderparameter enddo endif elseif ( opts%verbosity .gt. 0 ) then do i=1,uc%na do j=1,uc%isotope(i)%n write(*,"(' isotope: ',I2,' concentration: ',F8.5,' mass: ',F12.6)") & j,uc%isotope(i)%conc(j),uc%isotope(i)%mass(j) enddo write(*,"(' element: ',A2,' mean mass: ',F12.6,' mass disorder parameter',F12.9)") & trim(uc%atomic_symbol( uc%species(i) )),uc%isotope(i)%mean_mass,& uc%isotope(i)%disorderparameter enddo endif ! Read the force constants call fc%readfromfile(uc,'infile.forceconstant') if ( mw%talk ) write(*,*) '... read second order forceconstant' call fct%readfromfile(uc,'infile.forceconstant_thirdorder') if ( mw%talk ) write(*,*) '... read third order forceconstant' ! Get q-point mesh call lo_generate_qmesh(qp,uc,opts%qgrid,'fft',timereversal=opts%timereversal,nosym=.not.opts%qpsymmetry,verbosity=opts%verbosity) ! Distribute the q-points over MPI call qp%divide_across_mpi(mw%comm) ! Get frequencies in the whole BZ if ( mw%talk ) then write(*,*) '... getting the full dispersion relations' endif call dr%generate(qp,fc,uc,correctionlevel=2,timereversal=opts%timereversal,mpi_communicator=mw%comm,verbosity=opts%verbosity) ! Also the phonon DOS, for diagnostics call pd%generate(dr,qp,uc,opts%mfppts*2,opts%integrationtype,opts%sigma,opts%verbosity,mw%comm) ! Make sure it's stable if ( dr%is_unstable(qp) ) then write(*,*) '' write(*,*) 'FOUND UNSTABLE MODES. WILL STOP NOW.' call mpi_barrier(mw%comm,mw%error) call mpi_finalize(lo_status) stop endif #ifdef AGRESSIVE_SANITY ! make an initial memory report if ( mw%talk ) then write(*,*) '' write(*,*) 'MEMORY DIAGNOSTICS (mean,max,min in MB)' endif mem_per_rank=0.0_flyt mem_per_rank(mw%r+1)=uc%size_in_mem()*toMB call mpi_allreduce(MPI_IN_PLACE,mem_per_rank,mw%n,MPI_DOUBLE_PRECISION,MPI_SUM,mw%comm,mw%error) if ( mw%talk ) then write(*,"(1X,A,3(2X,F12.5))") ' structure:',lo_mean(mem_per_rank),maxval(mem_per_rank),minval(mem_per_rank) endif mem_per_rank=0.0_flyt mem_per_rank(mw%r+1)=qp%size_in_mem(qp)*toMB call mpi_allreduce(MPI_IN_PLACE,mem_per_rank,mw%n,MPI_DOUBLE_PRECISION,MPI_SUM,mw%comm,mw%error) if ( mw%talk ) then write(*,"(1X,A,3(2X,F12.5))") ' q-mesh:',lo_mean(mem_per_rank),maxval(mem_per_rank),minval(mem_per_rank) endif mem_per_rank=0.0_flyt mem_per_rank(mw%r+1)=dr%size_in_mem()*toMB call mpi_allreduce(MPI_IN_PLACE,mem_per_rank,mw%n,MPI_DOUBLE_PRECISION,MPI_SUM,mw%comm,mw%error) if ( mw%talk ) then write(*,"(1X,A,3(2X,F12.5))") ' dispersions:',lo_mean(mem_per_rank),maxval(mem_per_rank),minval(mem_per_rank) endif mem_per_rank=0.0_flyt mem_per_rank(mw%r+1)=pd%size_in_mem()*toMB call mpi_allreduce(MPI_IN_PLACE,mem_per_rank,mw%n,MPI_DOUBLE_PRECISION,MPI_SUM,mw%comm,mw%error) if ( mw%talk ) then write(*,"(1X,A,3(2X,F12.5))") ' phonon dos:',lo_mean(mem_per_rank),maxval(mem_per_rank),minval(mem_per_rank) endif #endif ! now I have all harmonic things, stop the init timer timer_init=walltime()-timer_init end block initharmonic ! Get the integration weights and matrix elements weights_elements: block real(flyt) :: t0 t0=walltime() timer_count=walltime() call lo_find_all_scattering_events(sc,qp,dr,uc,mw,opts%sigma,opts%thres,opts%integrationtype,& opts%correctionlevel,opts%mfp_max,opts%isotopescattering) call mpi_barrier(mw%comm,mw%error) ! stop counting timer, start matrixelement timer timer_count=walltime()-timer_count timer_matrixelements=walltime() ! Calculate scattering amplitudes t0=walltime() call calculate_scattering_amplitudes(uc,qp,sc,dr,fct,mw) call mpi_barrier(mw%comm,mw%error) ! stop matrix element timer, start some other timer timer_matrixelements=walltime()-timer_matrixelements if ( mw%talk ) write(*,*) 'Counted and got scattering amplitudes in ',tochar( walltime()-t0 ) end block weights_elements ! Make space and initialize everything to calculate thermal conductivity initkappa: block integer :: i ! Get thermal conductivity allocate(thermal_cond(10,opts%trangenpts)) thermal_cond=0.0_flyt ! The temperature axis allocate(temperatures(opts%trangenpts)) if ( opts%logtempaxis ) then call lo_logspace(opts%trangemin,opts%trangemax,temperatures) else call lo_linspace(opts%trangemin,opts%trangemax,temperatures) endif ! Setup the mean-free-path vs kappa plots, as well as the thin-film plots ! how many points on the x-axis? mf%np=opts%mfppts ! how many temperature? mf%nt=opts%trangenpts ! one plot for each temperature allocate(mf%temp(mf%nt)) ! thin film conductivity !if ( opts%thinfilm ) then ! allocate(mf%film(mf%nt)) !endif ! Make some space to keep intermediate values do i=1,qp%nq_irr allocate(dr%iq(i)%p_plus( dr%nb )) allocate(dr%iq(i)%p_minus( dr%nb )) allocate(dr%iq(i)%p_iso( dr%nb )) allocate(dr%iq(i)%qs( dr%nb )) allocate(dr%iq(i)%linewidth( dr%nb )) allocate(dr%iq(i)%F0( 3,dr%nb )) allocate(dr%iq(i)%Fn( 3,dr%nb )) allocate(dr%iq(i)%mfp( 3,dr%nb )) allocate(dr%iq(i)%scalar_mfp( dr%nb )) dr%iq(i)%linewidth=0.0_flyt dr%iq(i)%p_plus=0.0_flyt dr%iq(i)%p_minus=0.0_flyt dr%iq(i)%p_iso=0.0_flyt dr%iq(i)%qs=0.0_flyt dr%iq(i)%F0=0.0_flyt dr%iq(i)%Fn=0.0_flyt dr%iq(i)%mfp=0.0_flyt dr%iq(i)%scalar_mfp=0.0_flyt enddo do i=1,qp%nq_tot allocate(dr%aq(i)%kappa(3,3,dr%nb)) dr%aq(i)%kappa=0.0_flyt enddo end block initkappa ! Iteratively solve the BTE for each temperature. Additionally calculate mean free path plots ! and things like that. getkappa: block real(flyt), dimension(3,3) :: kappa,m0 real(flyt) :: t0 integer :: i if ( mw%talk ) then write(*,*) '' write(*,*) 'THERMAL CONDUCTIVITY' if ( opts%scfiterations .eq. 0 ) then write(*,"(3X,A11,6(1X,A14))") 'Temperature','kxx ','kyy ','kzz ','kxy ','kxz ','kyz ' endif endif ! Main loop over temperatures to solve the BTE do i=1,opts%trangenpts ! I might get a silly tiny temperature: if ( temperatures(i) .lt. lo_temperaturetol ) then kappa=0.0_flyt thermal_cond(1,i)=temperatures(i) thermal_cond(2:10,i)=0.0_flyt cycle endif ! Scattering rates t0=walltime() call calculate_qs(qp,sc,dr,temperatures(i),mw) timer_qs=timer_qs+walltime()-t0 ! Get the self-consistent solution call mpi_barrier(mw%comm,mw%error) if ( opts%scfiterations .gt. 0 ) then if ( mw%talk ) then write(*,*) '' write(*,*) 'Temperature: ',tochar(temperatures(i)) write(*,"(1X,A4,6(1X,A14),2X,A10)") 'iter','kxx ','kyy ','kzz ','kxy ','kxz ','kyz ','DeltaF/F' endif t0=walltime() call get_selfconsistent_solution(sc,dr,qp,uc,mw,temperatures(i),opts%scfiterations,opts%scftol) timer_scf=timer_scf+walltime()-t0 call get_kappa(dr,qp,uc,temperatures(i),kappa) m0=kappa*lo_kappa_au_to_SI if ( mw%talk ) write(*,"(5X,6(1X,F14.4),2X,E10.3)") m0(1,1),m0(2,2),m0(3,3),m0(1,2),m0(1,3),m0(2,3) else call get_kappa(dr,qp,uc,temperatures(i),kappa) m0=kappa*lo_kappa_au_to_SI if ( mw%talk ) write(*,"(1X,F12.3,6(1X,F14.4),2X,E10.3)") temperatures(i),m0(1,1),m0(2,2),m0(3,3),m0(1,2),m0(1,3),m0(2,3) endif ! Store thermal conductivity tensor thermal_cond(1,i)=temperatures(i) thermal_cond(2,i)=kappa(1,1)*lo_kappa_au_to_SI thermal_cond(3,i)=kappa(2,2)*lo_kappa_au_to_SI thermal_cond(4,i)=kappa(3,3)*lo_kappa_au_to_SI thermal_cond(5,i)=kappa(1,3)*lo_kappa_au_to_SI thermal_cond(6,i)=kappa(2,3)*lo_kappa_au_to_SI thermal_cond(7,i)=kappa(1,2)*lo_kappa_au_to_SI thermal_cond(8,i)=kappa(3,1)*lo_kappa_au_to_SI thermal_cond(9,i)=kappa(3,2)*lo_kappa_au_to_SI thermal_cond(10,i)=kappa(2,1)*lo_kappa_au_to_SI ! Calculate the cumulative plots t0=walltime() call mpi_barrier(mw%comm,mw%error) call get_cumulative_plots(mf%temp(i),qp,dr,pd,uc,opts%mfppts,temperatures(i),opts%sigma,kappa,mw) ! And the thin film thing !if ( mw%r .eq. 0 .and. opts%thinfilm ) then ! call calculate_thinfilm_kappa(mf%film(i),qp,dr,mf%np,uc,temperatures(i),& ! [0.0_flyt,0.0_flyt,1.0_flyt],[0.0_flyt,1.0_flyt,0.0_flyt]) !endif timer_cumulative=timer_cumulative+walltime()-t0 call mpi_barrier(mw%comm,mw%error) enddo end block getkappa ! dump things to file and print timings finalize_and_write: block real(flyt) :: t0 ! Write thermal conductivity to file if ( mw%talk ) call lo_dump_gnuplot_2d_real(thermal_cond,'outfile.thermal_conductivity',& ylabel='\kappa W/mK',xlabel='Temperature (K)') ! Write the cumulative kappa if ( mw%talk ) call write_cumulative_plots(mf,pd,'thz') ! Maybe dump data on a grid if ( mw%talk .and. opts%dumpgrid .and. opts%trangenpts .eq. 1 ) then write(*,*) '... dumping auxiliary data to files:' call dr%write_to_hdf5( qp,uc,'outfile.grid_thermal_conductivity.hdf5',temperatures(1) ) endif ! sum up the total time if ( mw%talk ) tt0=walltime()-tt0 ! Print timings if ( mw%talk ) then t0=timer_init+timer_count+timer_matrixelements+timer_qs+timer_kappa+timer_scf+timer_cumulative write(*,*) ' ' write(*,*) 'Timings:' write(*,"(A,F12.3,A,F7.3,A)") ' initialization:',timer_init, ' s, ', real(timer_init*100/tt0),'%' write(*,"(A,F12.3,A,F7.3,A)") ' integration weights:',timer_count, ' s, ', real(timer_count*100/tt0),'%' write(*,"(A,F12.3,A,F7.3,A)") ' matrix elements:',timer_matrixelements,' s, ', real(timer_matrixelements*100/tt0),'%' write(*,"(A,F12.3,A,F7.3,A)") ' QS calculation:',timer_qs, ' s, ', real(timer_qs*100/tt0),'%' write(*,"(A,F12.3,A,F7.3,A)") ' kappa:',timer_kappa, ' s, ', real(timer_kappa*100/tt0),'%' write(*,"(A,F12.3,A,F7.3,A)") ' self consistency:',timer_scf, ' s, ', real(timer_scf*100/tt0),'%' write(*,"(A,F12.3,A,F7.3,A)") ' cumulative plots:',timer_cumulative, ' s, ', real(timer_cumulative*100/tt0),'%' write(*,"(A,F12.3,A)") ' total:',tt0,' seconds' endif end block finalize_and_write ! And we are done! call mpi_barrier(mw%comm,mw%error) call mpi_finalize(lo_status) end program