#include "precompilerdefinitions" program phonon_dispersion_relations !!{!src/phonon_dispersion_relations/manual.md!} use konstanter, only: flyt,lo_tol,lo_hartree_to_eV use helpers, only: walltime,lo_mpi_helper,tochar,lo_chop,lo_mean,MPI_IN_PLACE,MPI_DOUBLE_PRECISION,MPI_SUM use dump_data, only: lo_dump_gnuplot_2d_real use options, only: lo_opts use type_crystalstructure, only: lo_crystalstructure 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,lo_read_qmesh_from_file use type_phonon_dispersions, only: lo_phonon_dispersions use type_phonon_dos, only: lo_phonon_dos use type_phonon_bandstructure, only: lo_phonon_bandstructure ! local use densityplots use velocitydos ! implicit none ! for normal phonon things type(lo_opts) :: opts type(lo_mpi_helper) :: mw type(lo_forceconstant_secondorder) :: fc type(lo_forceconstant_thirdorder) :: fct type(lo_crystalstructure) :: uc type(lo_phonon_bandstructure) :: bs class(lo_qpoint_mesh), allocatable :: qp type(lo_phonon_dispersions) :: dr type(lo_phonon_dos) :: pd ! for phonon dos of alloys real(flyt) :: timer #ifdef AGRESSIVE_SANITY real(flyt), parameter :: toMB=1.0_flyt/1024.0_flyt**2 real(flyt), dimension(:), allocatable :: mem_per_rank #endif ! Read things from file init: block ! init MPI call mw%init() ! Get the command line arguments timer=walltime() call opts%parse() if ( mw%talk .eqv. .false. ) opts%verbosity=-10 ! Read the crystal structure and make sure I have all the symmetry information call uc%readfromfile('infile.ucposcar',verbosity=opts%verbosity) call uc%classify('wedge',timereversal=opts%timereversal) ! Maybe non-default isotope distribution if ( opts%readiso ) then if ( mw%talk ) write(*,*) '... reading isotope distribution from file' call uc%readisotopefromfile() endif ! Read the force constant call fc%readfromfile(uc,'infile.forceconstant',opts%verbosity) ! Perhaps some Gruneisen parameters, in that case we need third order force constants. if ( opts%gruneisen ) then call fct%readfromfile(uc,'infile.forceconstant_thirdorder') endif #ifdef AGRESSIVE_SANITY ! make space for the thing that holds memory per rank lo_allocate(mem_per_rank(mw%n)) mem_per_rank=0.0_flyt ! make a 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)=fc%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))") ' forceconstant:',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)=fct%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))") ' thirdorder forceconstant:',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)=bs%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))") ' bandstructure:',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 end block init ! Dispersions along a path? path: block real(flyt), dimension(3) :: v1 real(flyt) :: t0 integer :: i,verb character(len=1000) :: opf ! Get the dispersions on the path t0=walltime() if ( uc%na .ge. 4 .and. mw%talk ) then verb=opts%verbosity+2 else verb=opts%verbosity endif if ( opts%gruneisen ) then call bs%generate(uc,fc,timereversal=opts%timereversal,fct=fct,verbosity=verb,npts=opts%nq,readpathfromfile=opts%readpathfromfile,mpi_communicator=mw%comm) else call bs%generate(uc,fc,timereversal=opts%timereversal,verbosity=verb,npts=opts%nq,readpathfromfile=opts%readpathfromfile,correctionlevel=2,mpi_communicator=mw%comm) endif ! Write a small summary if ( mw%talk ) then write(*,*) ' ' write(*,*) ' Settings: ' write(*,*) ' mode Gruneisen parameters: ',opts%gruneisen write(*,*) ' number of bands: ',tochar(bs%nb) write(*,*) ' number of paths: ',tochar(bs%npath) write(*,*) ' number of points on each path: ',tochar(bs%npts) write(*,*) ' ' do i=1,bs%npath write(*,*) ' path '//tochar(i)//': from '//trim(bs%symb_q_start(i))//' to '//trim(bs%symb_q_end(i)) opf="(1X,' starting point: cartesian',3(F9.6,' '),'fractional',3(F9.6,' '))" v1=lo_chop( uc%cartesian_to_fractional( bs%segment(i)%r1, reciprocal=.true.,pbc=.false.) , lo_sqtol) write(*,opf) bs%segment(i)%r1,v1 opf="(1X,' ending point: cartesian',3(F9.6,' '),'fractional',3(F9.6,' '))" v1=lo_chop( uc%cartesian_to_fractional( bs%segment(i)%r2, reciprocal=.true.,pbc=.false.) , lo_sqtol) write(*,opf) bs%segment(i)%r2,v1 enddo write(*,*) ' ' write(*,"(' ... got the phonon bandstructure in ',F12.9,'s' )") walltime()-t0 ! Dump things to files for plotting call bs%write_dispersive_property(opts%enhet,'frequency','outfile.dispersion_relations',.false.) call bs%write_dispersive_property(opts%enhet,'groupvelocity','outfile.group_velocities',.false.) if ( opts%gruneisen ) call bs%write_dispersive_property(opts%enhet,'gruneisen','outfile.mode_gruneisen_parameters',.false.) endif end block path ! Get the full mesh and the density of states? dos: block real(flyt), dimension(:), allocatable :: temperatures real(flyt) :: f0,f1,f2,f3,t0 integer :: i,u,verb ! The dispersions along the path are done. Perhaps we want the full dispersion relations for something? if ( opts%fullmesh ) then if ( mw%talk ) then write(*,*) ' ' write(*,*) ' Calculating phonon dispersion across the Brillouin zone' endif ! Get a q-point mesh if ( opts%readqmesh ) then call lo_read_qmesh_from_file(qp,uc,'infile.qgrid.hdf5',verbosity=opts%verbosity) 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,mw=mw) end select endif ! Get the full dispersion relations t0=walltime() if ( uc%na .ge. 4 .and. mw%talk ) then verb=opts%verbosity+2 else verb=opts%verbosity endif call dr%generate(qp,fc,uc,verbosity=verb,timereversal=opts%timereversal,mpi_communicator=mw%comm) endif ! Perhaps we want the phonon density of states and related things if ( opts%dos ) then t0=walltime() if ( mw%talk .eqv. .false. ) then verb=0 else verb=max(1,opts%verbosity) endif call pd%generate(dr,qp,uc,verbosity=verb,integrationtype=opts%dosintegrationtype,adaptiveparameter=opts%sigma,& ndos=opts%dospoints,mpi_communicator=mw%comm) ! Dump it to file if ( mw%talk ) then call pd%write_to_file(uc,opts%enhet,'outfile.phonon_dos') call pd%write_to_hdf5(opts%enhet,'outfile.phonon_dos.hdf5') endif ! Get the density-plot thing !if ( mw%talk ) call omega_vs_norm_q_density(qp,dr,uc) endif ! A positive number of temperatures means I should dump energy. I do this not in parallel. if ( opts%trangenpts .gt. 0 .and. mw%talk ) then if ( .not. opts%fullmesh ) then write(*,*) 'I screwed up the heuristics' stop endif write(*,*) '... calculating free energy' lo_allocate(temperatures(opts%trangenpts)) call lo_linspace(opts%trangemin,opts%trangemax,temperatures) ! dump some free energies and entropies write(*,*) '' write(*,*) ' T(K) F(eV/atom) S(eV/K/atom) Cv(eV/K/atom)' u=open_file('out','outfile.free_energy') do i=1,opts%trangenpts f0=temperatures(i) f1=dr%phonon_free_energy(f0)*lo_hartree_to_eV f2=dr%phonon_entropy(f0)*lo_hartree_to_eV f3=dr%phonon_cv(f0)*lo_hartree_to_eV !f4=dr%phonon_free_energy_classical(f0) if ( abs(f1-123456789.0_flyt) .lt. lo_tol ) then ! return NaN if I have imaginary frequencies write(u,*) f0,' NaN NaN NaN' !' NaN' write(*,*) f0,' NaN NaN NaN' !' NaN' else write(u,"(1X,F12.5,4(1X,E18.11))") temperatures(i),f1,f2,f3 !,f4 write(*,"(1X,F12.5,4(1X,E18.11))") temperatures(i),f1,f2,f3 !,f4 endif enddo close(u) endif end block dos ! Print some file and wrap it up wrapup: block ! Write everything on the grid to file, for some inexplicable reason. if ( opts%dumpgrid ) then ! small sanity check if ( .not. opts%fullmesh ) then write(*,*) 'I screwed up the heuristics' stop endif if ( opts%gruneisen ) then call dr%gruneisen(qp,fct,uc,opts%verbosity) endif ! actually write it if ( mw%talk ) then call dr%write_to_hdf5(qp,uc,'outfile.grid_dispersions.hdf5') write(*,*) '... wrote data for the full grid' endif endif ! Dump everything to hdf5 if ( mw%talk ) then if ( opts%sortbands ) call bs%sort_bands() call bs%write_to_hdf5(uc,opts%enhet,'outfile.dispersion_relations.hdf5') endif #ifdef AGRESSIVE_SANITY ! make a 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)=fc%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))") ' forceconstant:',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)=fct%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))") ' thirdorder forceconstant:',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)=bs%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))") ' bandstructure:',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 if ( mw%talk ) then write(*,*) ' ' write(*,*) 'All done in ',tochar(walltime()-timer),'s' endif call mw%destroy() end block wrapup end program