merge trunk

src/canonical_configuration/main.f90

@@ -210,6 +210,12 @@ program canonical_configuration
         case (5) ! Siesta output
             fname = 'siesta_conf'//tochar(iconf, 4)
             call p%writetofile(trim(fname), opts%output_format, write_velocities=.true.)
+        case (6) ! QE output
+            fname = 'qe_conf'//tochar(iconf, 4)
+            call p%writetofile(trim(fname), opts%output_format, write_velocities=.false.)
+        case (7) ! QE output
+            fname = 'parsec_conf'//tochar(iconf, 4)
+            call p%writetofile(trim(fname), opts%output_format, write_velocities=.false.)
         end select
 
         ! just measure some stuff, for no good reason.

src/canonical_configuration/options.f90

@@ -66,9 +66,9 @@ subroutine parse(opts)
                  help_markdown="That is, use \( \sqrt{\frac{\hbar (2n+1) }{2 m \omega}} \) as the mean normal mode amplitudes instead of the classical \( \frac{1}{\omega}\sqrt{\frac{k_BT}{m}} \)", &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--output_format', switch_ab='-of', &
-                 help='Selects output format. 1 is VASP, 2 is Abinit, 4 is FHI-Aims, 5 is Siesta. Default 1.', &
-                 required=.false., act='store', def='1', choices='1,2,4,5', error=lo_status)
+    call cli%add(switch='--output_format', switch_ab='-of', hidden=.false., &
+                 help='Output format. 1 is VASP, 2 Abinit, 4 FHI-Aims, 5 Siesta, 6 QE and 7 Parsec', &
+                 required=.false., act='store', def='1', choices='1,2,4,5,6,7', error=lo_status)
     if (lo_status .ne. 0) stop
     call cli%add(switch='--mindist', &
                  help='What is the smallest distance between two atoms allowed, in units of the nearest neighbour distance.', &

src/dump_dynamical_matrices/Makefile

@@ -23,7 +23,7 @@ clean:
 	rm -f $(PROG) $(OBJS) *.mod $(OBJECT_PATH)*.mod
 
 
-$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o
+$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o ../../lib/libolle.a
 	$(F90) $(F90FLAGS) -c main.f90 -o $@
 $(OBJECT_PATH)options.o:
 	$(F90) $(F90FLAGS) -c options.f90 -o $@

src/dump_dynamical_matrices/main.f90

@@ -78,6 +78,8 @@ program dump_dynamical_matrices
 
 call fc%write_to_anaddb(uc, opts%qgrid, mw, mem)
 
+call fc%write_dynmat_to_qe(uc,opts%qgrid,mw,mem)
+
 !else
 !    ! automagically generate a grid. Probably a bad idea, since my version of gridgeneration and
 !    ! Abinits might be different.

src/extract_forceconstants/Makefile

@@ -26,7 +26,7 @@ clean:
 .f90.o:
 	$(F90) $(F90FLAGS) -c $< $(LIBS)
 
-$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o
+$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o
 	$(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@
-$(OBJECT_PATH)options.o:
+$(OBJECT_PATH)options.o: options.f90
 	$(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@

src/extract_forceconstants/main.f90

@@ -348,6 +348,11 @@ program extract_forceconstants
             do i = 1, 6
                 write (*, "(6(3X,F15.5))") lo_chop((fc2%elastic_constants_voigt(:, i) + fc2%elastic_constants_voigt_longrange(:, i))*lo_pressure_HartreeBohr_to_GPa, lo_tol)
             end do
+
+            if ( opts%dumpepw ) then
+                ! If desired, write IFCs in EPW format
+                call fc2%write_to_qe(uc,mw,mem)
+            endif
         end if
     end if
     if (map%have_fc_triplet) then

src/extract_forceconstants/options.f90

@@ -68,7 +68,8 @@ module options
     logical :: fake_dielectric = .false.
     ! developer mode
     logical :: devmode = .false.
-
+    ! Dump IFCs in EPW format
+    logical :: dumpepw = .false.
 contains
     procedure :: parse
 end type
@@ -248,6 +249,9 @@ subroutine parse(opts)
     call cli%add(switch='--developermode', switch_ab='-dev', hidden=.true., help='dev. mode', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
+    call cli%add(switch='--output_epw', switch_ab='-epw', hidden=.false., help='Write the (second order) interatomic force constants in XML format useable by EPW.', &
+                 required=.false., act='store_true', def='.false.', error=lo_status)
+    if (lo_status .ne. 0) stop
 
     ! actually parse it
     call cli%parse(error=lo_status)
@@ -306,6 +310,7 @@ subroutine parse(opts)
     call cli%get(switch='--temperature', val=opts%temperature)
     call cli%get(switch='--fakediel', val=opts%fake_dielectric)
     call cli%get(switch='--developermode', val=opts%devmode)
+    call cli%get(switch='--output_epw', val=opts%dumpepw)
 
     if (lo_status .ne. 0) stop
 

src/generate_structure/main.f90

@@ -113,6 +113,13 @@ program generate_structure
             ss%latticevectors = ss%latticevectors
             ss%inv_latticevectors = ss%inv_latticevectors
             call ss%writetofile('outfile.supercell_lammps_ipi', 3, transformationmatrix=tm)
+        case (6) ! QE
+            call ss%writetofile('outfile.supercell_qe', opts%outputformat)
+            write (*, *) '... wrote supercell in Quantum ESPRESSO format'
+        case (7) ! Parsec
+            call uc%writetofile('outfile.unitcell_parsec', opts%outputformat)
+            call ss%writetofile('outfile.supercell_parsec', opts%outputformat)
+            write (*, *) '... wrote supercell in PARSEC format'
         end select
     end if
 end block getsupercell

src/generate_structure/options.f90

@@ -64,8 +64,8 @@ subroutine parse(opts)
                  required=.false., act='store', def='5.0', error=lo_status)
     if (lo_status .ne. 0) stop
     call cli%add(switch='--output_format', switch_ab='-of', hidden=.false., &
-                 help='Output format. 1 is VASP, 2 Abinit, 4 FHI-Aims and 5 Siesta', &
-                 required=.false., act='store', def='1', choices='1,2,4,5', error=lo_status)
+                 help='Output format. 1 is VASP, 2 Abinit, 4 FHI-Aims, 5 Siesta, 6 QE and 7 Parsec', &
+                 required=.false., act='store', def='1', choices='1,2,4,5,6,7', error=lo_status)
     if (lo_status .ne. 0) stop
     call cli%add(switch='--thirdorder_cutoff', switch_ab='-rc3', hidden=.true., &
                  help='Cutoff for the third order force constants', &

src/libolle/type_crystalstructure_io.f90

@@ -376,8 +376,12 @@ module subroutine writetofile(p,filename,output_format,write_velocities,transfor
             end if
         case(4) ! FHI Aims
             call writetofile_aims(p,filename,write_velocities)
-        case(5) ! Siesta 
+        case(5) ! Siesta
             call writetofile_siesta(p,filename,write_velocities)
+        case(6) ! QE
+            call writetofile_qe(p,filename,write_velocities)
+        case(7) ! Parsec
+            call writetofile_parsec(p,filename,write_velocities)
         case default
             call lo_stop_gracefully(['Unknown output format: '//tochar(output_format)],lo_exitcode_io,__FILE__,__LINE__)
     end select
@@ -726,7 +730,7 @@ subroutine writetofile_siesta(p,filename,write_velocities)
         ! I don't know how to write velocities in Siesta. Someone should tell me that.
         ! MJV: fixed at least partly, by putting them in XV files (pos, veloc).
         ! not sure you can add them to the main fdf input files as well...
-        ! NB: there is a flag to use the XV file MD.UseSaveXV .true. but you have to fill 
+        ! NB: there is a flag to use the XV file MD.UseSaveXV .true. but you have to fill
         !   the coordinates block anyway, so might as well fill both here.
         !   Main advantage is that you can add the velocities in XV, but for configuration
         !   generation it's not important, you just get forces out.
@@ -785,7 +789,7 @@ subroutine writetofile_siesta(p,filename,write_velocities)
         write(u,*) "XC.functional GGA"
         write(u,*) "XC.authors    PBE"
         write(u,*)
- 
+
         write(u,*) "%block kgrid_Monkhorst_Pack"
         write(u,*) "   1   0   0   0."
         write(u,*) "   0   1   0   0."
@@ -803,25 +807,114 @@ subroutine writetofile_siesta(p,filename,write_velocities)
 
         if ( filename .ne. 'stdout' ) close(u)
 
-! add a XV format file as well: Uses bohr for the cartesian positions.
-!    what are the units for the velocity?
-!
-! latice as 3x3 block, with additional 3x3 block of 0s after
-! natom
-! itype zatom xcoord(3) veloc(3)
-!
-! NB: I _think_ they use Angstr for the positions, needs to be checked. 
-!    for the velocities, no idea what they expect...
+        ! add a XV format file as well: Uses bohr for the cartesian positions.
+        !    what are the units for the velocity?
+        !
+        ! latice as 3x3 block, with additional 3x3 block of 0s after
+        ! natom
+        ! itype zatom xcoord(3) veloc(3)
+        !
+        ! NB: I _think_ they use Angstr for the positions, needs to be checked.
+        !    for the velocities, no idea what they expect...
         u=open_file('out',filename//".XV")
         write(u,'(3E20.10,2x,3E20.10)') lo_choplarge(p%latticevectors(:,1),lo_sqtol), 0.000000000, 0.000000000, 0.000000000
         write(u,'(3E20.10,2x,3E20.10)') lo_choplarge(p%latticevectors(:,2),lo_sqtol), 0.000000000, 0.000000000, 0.000000000
         write(u,'(3E20.10,2x,3E20.10)') lo_choplarge(p%latticevectors(:,3),lo_sqtol), 0.000000000, 0.000000000, 0.000000000
-        write(u,*) tochar(p%na) 
+        write(u,*) tochar(p%na)
         opf="(1X,2(I5),2X,2(3F18.9,'  '))"
         do i=1,p%na
             write(u,opf) p%species(i), symbol_to_z(p%atomic_symbol(p%species(i))), lo_chop(p%rcart(:,i),lo_sqtol), lo_chop(p%v(:,i),lo_sqtol)
         enddo
     end subroutine
+
+    !> Writes a structure to QE output file or stdout
+    subroutine writetofile_qe(p,filename,write_velocities)
+        !> crystal structure
+        class(lo_crystalstructure), intent(in) :: p
+        !> filename
+        character(len=*), intent(in) :: filename
+        !> should be written. Default false.
+        logical, intent(in), optional :: write_velocities
+        ! local
+        integer :: i
+
+        real(flyt) :: latpar
+        character(len=1000) :: opf
+        integer :: u
+
+        ! file or stdout
+        if ( filename .eq. 'stdout' ) then
+            u=0
+        else
+            u=open_file('out',filename)
+        endif
+        write(u,*) "    nat             = ", tochar(p%na)
+        write(u,*) "    ntyp            = ", tochar(p%nelements)
+        write(u,*) "/"
+        write(u,*) "CELL_PARAMETERS bohr"
+        write(u,*) lo_choplarge(p%latticevectors(:,1),lo_sqtol)
+        write(u,*) lo_choplarge(p%latticevectors(:,2),lo_sqtol)
+        write(u,*) lo_choplarge(p%latticevectors(:,3),lo_sqtol)
+        write(u,*) " "
+        write(u,*) "ATOMIC_POSITIONS bohr"
+        opf="(1X,a,2X,3F18.9))"
+        do i=1, p%na
+          write(u,opf) trim(p%atomic_symbol(p%species(i))), lo_chop(p%rcart(:,i),lo_sqtol)
+        enddo
+        if ( filename .ne. 'stdout' ) close(u)
+    end subroutine
+
+    !> Writes a structure in PARSEC format
+    subroutine writetofile_parsec(p,filename,write_velocities)
+        !> crystal structure
+        class(lo_crystalstructure), intent(in) :: p
+        !> filename
+        character(len=*), intent(in) :: filename
+        !> should velocities be written. Default false.
+        logical, intent(in), optional :: write_velocities
+        ! local
+        integer :: i
+        !
+        real(flyt) :: latpar
+        logical :: vel
+        integer :: u
+
+        ! I don't have a velocity input format in parsec yet, so later problem
+        if ( present(write_velocities) ) then
+            vel=write_velocities
+        else
+            vel=.false.
+        endif
+
+        ! file or stdout
+        if ( filename .eq. 'stdout' ) then
+            u=0
+        else
+            u=open_file('out',filename)
+        endif
+
+        write(u,*) "# comment line"
+        write(u,*) "Bulk"
+        if ( p%info%unitcell_lattice_parameter .gt. 0.0_flyt ) then
+            latpar=p%info%unitcell_lattice_parameter
+            write(u,"(2X,F21.13)") latpar*lo_bohr_to_A
+            write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,1)/latpar,1E-12_flyt)
+            write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,2)/latpar,1E-12_flyt)
+            write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,3)/latpar,1E-12_flyt)
+        else
+            latpar=1.0_flyt
+            write(u,"(2X,F20.12)") latpar
+            write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,1)/latpar,1E-12_flyt)
+            write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,2)/latpar,1E-12_flyt)
+            write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,3)/latpar,1E-12_flyt)
+        endif
+        do i=1,p%na
+            write(u,"(1X,3(1X,E19.12),1X,A)") p%r(:,i),trim(p%atomic_symbol( p%species(i) ))
+        enddo
+
+        if ( filename .ne. 'stdout' ) close(u)
+    end subroutine
+
 end subroutine
 
 !> Read isotope distribution from file
@@ -935,7 +1028,6 @@ module subroutine write_structure_to_hdf5(p,filename,input_id)
     !> in case I want to write it as a part of another file
     integer(HID_T), intent(in), optional :: input_id
 
-    integer :: i
     character(len=2000) :: atomnames
     type(lo_hdf5_helper) :: h5
 

src/libolle/type_forceconstant_secondorder.f90

@@ -157,6 +157,8 @@ module type_forceconstant_secondorder
     procedure :: fake_forceconstant
     !> dump the forceconstant in Abinit anaddb format
     procedure :: write_to_anaddb
+    !> dump forceconstants into QE format
+    procedure :: write_to_qe
     !> Set parameters for Ewald summation and enforce hermiticity of Born effective charges
     procedure :: set_ewald_and_enforce_borncharge_hermiticity
     !> Return a dynamical matrix
@@ -171,6 +173,8 @@ module type_forceconstant_secondorder
     procedure :: supercell_longrange_dynamical_matrix_at_gamma
     !> Get the commensurate modes for a supercell
     procedure :: commensurate_modes
+    !>
+    procedure :: write_dynmat_to_qe
     !> Size in memory, in bytes
     procedure :: size_in_mem => fc2_size_in_mem
 end type
@@ -286,6 +290,23 @@ module subroutine write_to_anaddb(fc, uc, qgrid, mw, mem)
         type(lo_mpi_helper), intent(inout) :: mw
         type(lo_mem_helper), intent(inout) :: mem
     end subroutine
+    module subroutine write_to_qe(fc,uc,mw,mem)
+        !> forceconstant
+        class(lo_forceconstant_secondorder), intent(inout) :: fc
+        !> unitcell
+        type(lo_crystalstructure), intent(inout) :: uc
+        type(lo_mpi_helper), intent(inout) :: mw
+        type(lo_mem_helper), intent(inout) :: mem
+    end subroutine
+    module subroutine write_dynmat_to_qe(fc,uc,qgrid,mw,mem)
+        !> forceconstant
+        class(lo_forceconstant_secondorder), intent(inout) :: fc
+        !> unitcell
+        type(lo_crystalstructure), intent(inout) :: uc
+        integer, dimension(3),intent(in) :: qgrid
+        type(lo_mpi_helper), intent(inout) :: mw
+        type(lo_mem_helper), intent(inout) :: mem
+    end subroutine
 end interface
 
 contains