From bc4a93e049433631864a43d8cbf69ecad42f8d05 Mon Sep 17 00:00:00 2001 From: Ravindra Shinde Date: Mon, 8 Feb 2021 09:39:18 +0100 Subject: [PATCH] initial template files for reading the input --- src/Makefile | 42 + src/parser.f90 | 285 ++++++ src/read_cards.f90 | 1938 +++++++++++++++++++++++++++++++++++++ src/read_input.f90 | 77 ++ src/read_namelists.f90 | 2055 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 4397 insertions(+) create mode 100644 src/Makefile create mode 100644 src/parser.f90 create mode 100644 src/read_cards.f90 create mode 100644 src/read_input.f90 create mode 100644 src/read_namelists.f90 diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..0603f0f --- /dev/null +++ b/src/Makefile @@ -0,0 +1,42 @@ +#/a makefile for modules + +include ../make.inc + +# location of needed modules +modflags=$(basemod_flags) \ + $(mod_flag)../elpa/src + +# list of modules + +modules = \ +io_files.o \ +io_global.o \ +parameters.o \ +parser.o \ +read_cards.o \ +read_input.o \ +read_namelists.o \ +io_base.o +#hdf5_qe.o\ +#h5_module.o\ + + +# list of subroutines and functions (not modules) + +objs = \ +inpfile.o \ +plot_io.o \ +test_input_file.o + +all : libio08.a + +## the following is needed only for lapack compiled from sources + +libio08.a: $(modules) $(objs) + $(ar) $(arflags) $@ $? + $(ranlib) $@ + +clean : + - /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.l + +include make.depend diff --git a/src/parser.f90 b/src/parser.f90 new file mode 100644 index 0000000..5ce07b4 --- /dev/null +++ b/src/parser.f90 @@ -0,0 +1,285 @@ +! +! +! ... subroutine field_count: accepts two string (one of them is optional) +! and one integer and count the number of fields +! in the string separated by a blank or a tab +! character. if the optional string is specified +! (it has anyway len=1) it is assumed as the +! separator character. +! ignores any character following the exclamation +! mark (fortran comment) +! +! ... subroutine con_cam: counts the number of fields in a string +! separated by the optional character +! +! ... subroutine field_compare: accepts two strings and one integer. counts the +! fields contained in the first string and +! compares it with the integer. +! if they are less than the integer calls the +! routine error and show by the second string the +! name of the field where read-error occurred. +! +! +!---------------------------------------------------------------------------- +module parser + !---------------------------------------------------------------------------- + ! + use io_global, only : stdout + use kinds, only : dp + ! + private + ! + public :: parse_unit, field_count, read_line, get_field + public :: version_parse, version_compare + ! + integer :: parse_unit = 5 ! normally 5, but can be set otherwise + ! + contains + ! + ! + !-------------------------------------------------------------------------- + subroutine field_count( num, line, car ) + !-------------------------------------------------------------------------- + ! + implicit none + ! + integer, intent(out) :: num + character(len=*), intent(in) :: line + character(len=1), optional, intent(in) :: car + character(len=1) :: sep1, sep2 + integer :: j + ! + ! + num = 0 + ! + if ( .not. present(car) ) then + ! + sep1 = char(32) ! ... blank character + sep2 = char(9) ! ... tab character + ! + do j = 2, max( len( line ), 256 ) + ! + if ( line(j:j) == '!' .or. line(j:j) == char(0) ) then + ! + if ( line(j-1:j-1) /= sep1 .and. line(j-1:j-1) /= sep2 ) then + ! + num = num + 1 + ! + end if + ! + exit + ! + end if + ! + if ( ( line(j:j) == sep1 .or. line(j:j) == sep2 ) .and. & + ( line(j-1:j-1) /= sep1 .and. line(j-1:j-1) /= sep2 ) ) then + ! + num = num + 1 + ! + end if + ! + end do + ! + else + ! + sep1 = car + ! + do j = 2, max( len( line ), 256 ) + ! + if ( line(j:j) == '!' .or. & + line(j:j) == char(0) .or. line(j:j) == char(32) ) then + ! + if ( line(j-1:j-1) /= sep1 ) num = num + 1 + ! + exit + ! + end if + ! + if ( line(j:j) == sep1 .and. line(j-1:j-1) /= sep1 ) num = num + 1 + ! + end do + ! + end if + ! + return + ! + end subroutine field_count + ! + ! + !-------------------------------------------------------------------------- + subroutine read_line( line, nfield, field, end_of_file, error ) + !-------------------------------------------------------------------------- + ! + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + use io_global, only : ionode, ionode_id + ! + implicit none + ! + character(len=*), intent(out) :: line + character(len=*), optional, intent(in) :: field + integer, optional, intent(in) :: nfield + logical, optional, intent(out) :: end_of_file, error + logical :: tend, terr + ! + ! + if( len( line ) < 256 ) then + call errore(' read_line ', ' input line too short ', max(len(line),1) ) + end if + ! + tend = .false. + terr = .false. + if ( ionode ) then +30 read (parse_unit, fmt='(a256)', err=15, end=10) line + if( line == ' ' .or. line(1:1) == '#' ) go to 30 + go to 20 +10 tend = .true. + go to 20 +15 terr = .true. +20 continue + end if + ! + call mp_bcast( tend, ionode_id, intra_image_comm ) + call mp_bcast( terr, ionode_id, intra_image_comm ) + call mp_bcast( line, ionode_id, intra_image_comm ) + ! + if( present(end_of_file) ) then + end_of_file = tend + else if( tend ) then + call infomsg(' read_line ', ' end of file ' ) + end if + if( present(error) ) then + error = terr + else if( terr ) then + call infomsg(' read_line ', ' read error ' ) + end if + if( present(field) .and. .not.(tend.or.terr) ) & + &call field_compare( line, nfield, field ) + ! + end subroutine read_line + ! + ! + !-------------------------------------------------------------------------- + subroutine field_compare( str, nf, var ) + !-------------------------------------------------------------------------- + ! + implicit none + ! + character(len=*), intent(in) :: var + integer, intent(in) :: nf + character(len=*), intent(in) :: str + integer :: nc + ! + call field_count( nc, str ) + ! + if( nc < nf ) & + call errore( ' field_compare ', & + & ' wrong number of fields: ' // trim( var ), 1 ) + ! + return + ! + end subroutine field_compare + ! + ! + !-------------------------------------------------------------------------- + subroutine con_cam(num, line, car) + !-------------------------------------------------------------------------- + character(len=*) :: line + character(len=1) :: sep + character(len=1), optional :: car + integer :: num, j + + num = 0 + if (len(line) .gt. 256 ) then + write( stdout,*) 'riga ', line + write( stdout,*) 'lunga ', len(line) + num = -1 + return + end if + + write( stdout,*) '1riga ', line + write( stdout,*) '1lunga ', len(line) + if ( .not. present(car) ) then + sep=char(32) !char(32) is the blank character + else + sep=car + end if + + do j=2, max(len(line),256) + if ( line(j:j) == '!' .or. line(j:j) == char(0)) then + return + end if + if ( (line(j:j) .eq. sep) .and. & + (line(j-1:j-1) .ne. sep) ) then + num = num + 1 + end if + end do + return + end subroutine con_cam + ! + !-------------------------------------------------------------------------- + subroutine get_field(n, field, str, sep) + !-------------------------------------------------------------------------- + ! extract whitespace-separated nth block from string + implicit none + integer,intent(in) :: n + character(len=*),intent(out) :: field + character(len=*),intent(in) :: str + character(len=1),optional,intent(in) :: sep + integer :: i,j,z ! block start and end + integer :: k ! block counter + character(len=1) :: sep1, sep2 + !print*, "------------- parser start -------------" + !print '(3a)', "string: -->", str,"<--" + if(present(sep)) then + sep1 = sep + sep2 = sep ! redundant, but easy + else + sep1 = char(32) ! ... blank character + sep2 = char(9) ! ... tab char + endif + ! + k = 1 ! counter for the required block + ! + do i = 1,len(str) + ! look for the beginning of the required block + z = max(i-1,1) + !print '(2a1,3i4,2l)', str(i:i), str(z:z), i,z,k,n,& + ! (str(i:i) == sep1 .or. str(i:i) == sep2), (str(z:z) /= sep1 .and. str(z:z) /= sep2) + if( k == n) exit + if( (str(i:i) == sep1 .or. str(i:i) == sep2) & + .and. & + (str(z:z) /= sep1 .and. str(z:z) /= sep2) & + ) & + k = k+1 + enddo + ! + !print*, "i found: ",i + do j = i,len(str) + ! look for the beginning of the next block + z = max(j-1,1) + if( (str(j:j) == sep1 .or. str(j:j) == sep2) & + .and. & + (str(z:z) /= sep1 .and. str(z:z) /= sep2) & + ) & + k = k+1 + if( k >n) exit + enddo + !print*, "j found: ",j + ! + if (j <= len(str)) then + ! if we are here, the reqired block was followed by a separator + ! and another field, we have to trash one char (a separator) + field = trim(adjustl(str(i:j-1))) + !print*, "taking: ",i,j-2 + else + ! if we are here, it was the last block in str, we have to take + ! all the remaining chars + field = trim(adjustl(str(i:len(str)))) + !print*, "taking from ",i + endif + !print*, "------------- parser end -------------" + + end subroutine get_field + +end module parser diff --git a/src/read_cards.f90 b/src/read_cards.f90 new file mode 100644 index 0000000..26505ad --- /dev/null +++ b/src/read_cards.f90 @@ -0,0 +1,1938 @@ +! +! copyright (c) 2002-2014 quantum espresso group +! this file is distributed under the terms of the +! gnu general public license. see the file `license' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +!--------------------------------------------------------------------------- +module read_cards_module + !--------------------------------------------------------------------------- + ! + ! ... this module handles the reading of cards from standard input + ! ... original version written by carlo cavazzoni + ! + use kinds, only : dp + use io_global, only : stdout + use wy_pos, only : wypos + use parser, only : field_count, read_line, get_field, parse_unit + use io_global, only : ionode, ionode_id + ! + use input_parameters + ! + ! + implicit none + ! + save + ! + private + ! + public :: read_cards + ! + ! ... end of module-scope declarations + ! + ! ---------------------------------------------- + ! +contains + ! + ! ... read cards .... + ! + ! ... subroutines + ! + !---------------------------------------------------------------------- + subroutine card_default_values( ) + !---------------------------------------------------------------------- + ! + use autopilot, only : init_autopilot + ! + implicit none + ! + ! + ! ... mask that control the printing of selected kohn-sham occupied + ! ... orbitals, default allocation + ! + call allocate_input_iprnks( 0, nspin ) + nprnks = 0 + ! + ! ... simulation cell from standard input + ! + trd_ht = .false. + rd_ht = 0.0_dp + ! + ! ... reference simulation cell from standard input + ! + ref_cell = .false. + rd_ref_ht = 0.0_dp + ! + ! ... constraints + ! + nconstr_inp = 0 + constr_tol_inp = 1.e-6_dp + ! + ! ... ionic mass initialization + ! + atom_mass = 0.0_dp + ! + ! ... k-points + ! + k_points = 'gamma' + tk_inp = .false. + nkstot = 1 + nk1 = 0 + nk2 = 0 + nk3 = 0 + k1 = 0 + k2 = 0 + k3 = 0 + ! + ! ... electronic states + ! + tf_inp = .false. + ! + ! ... ion_velocities + ! + tavel = .false. + ! + call init_autopilot() + ! + return + ! + end subroutine card_default_values + ! + ! + !---------------------------------------------------------------------- + subroutine read_cards ( prog, unit ) + !---------------------------------------------------------------------- + ! + use autopilot, only : card_autopilot + ! + implicit none + ! + integer, intent(in), optional :: unit + ! + character(len=2) :: prog ! calling program ( pw, cp, wa ) + character(len=256) :: input_line + character(len=80) :: card + character(len=1), external :: capital + logical :: tend + integer :: i + ! + ! read_line reads from unit parse_unit + ! + if (present(unit)) then + parse_unit = unit + else + parse_unit = 5 + end if + ! + call card_default_values( ) + ! +100 call read_line( input_line, end_of_file=tend ) + ! + if( tend ) goto 120 + if( input_line == ' ' .or. input_line(1:1) == '#' .or. & + input_line(1:1) == '!' ) goto 100 + ! + read (input_line, *) card + ! + do i = 1, len_trim( input_line ) + input_line( i : i ) = capital( input_line( i : i ) ) + enddo + ! + if ( trim(card) == 'autopilot' ) then + ! + call card_autopilot( input_line ) + if ( prog == 'pw' .and. ionode ) & + write( stdout,'(a)') 'warning: card '//trim(input_line)//' ignored' + ! + elseif ( trim(card) == 'atomic_species' ) then + ! + call card_atomic_species( input_line ) + ! + elseif ( trim(card) == 'atomic_positions' ) then + ! + call card_atomic_positions( input_line, prog ) + ! + elseif ( trim(card) == 'atomic_forces' ) then + ! + call card_atomic_forces( input_line ) + ! + elseif ( trim(card) == 'constraints' ) then + ! + call card_constraints( input_line ) + ! + elseif ( trim(card) == 'dipole' ) then + ! + call errore('read_cards','card dipole no longer existing',1) + ! + elseif ( trim(card) == 'esr' ) then + ! + call errore('read_cards','card esr no longer existing',1) + ! + elseif ( trim(card) == 'k_points' ) then + ! + if ( ( prog == 'cp' ) ) then + if( ionode ) & + write( stdout,'(a)') 'warning: card '//trim(input_line)//' ignored' + else + call card_kpoints( input_line ) + endif + ! + elseif ( trim(card) == 'additional_k_points' ) then + ! + call card_add_kpoints( input_line ) + + elseif ( trim(card) == 'occupations' ) then + ! + call card_occupations( input_line ) + ! + elseif ( trim(card) == 'cell_parameters' ) then + ! + call card_cell_parameters( input_line ) + ! + elseif ( trim(card) == 'ref_cell_parameters' ) then + ! + call card_ref_cell_parameters( input_line ) + ! + elseif ( trim(card) == 'atomic_velocities' ) then + ! + call card_ion_velocities( input_line ) + ! + elseif ( trim(card) == 'ksout' ) then + ! + call card_ksout( input_line ) + if ( ( prog == 'pw' ) .and. ionode ) & + write( stdout,'(a)') 'warning: card '//trim(input_line)//' ignored' + ! + elseif ( trim(card) == 'plot_wannier' ) then + ! + call card_plot_wannier( input_line ) + + elseif ( trim(card) == 'wannier_ac' .and. ( prog == 'wa' )) then + ! + call card_wannier_ac( input_line ) + + else + ! + if ( ionode ) & + write( stdout,'(a)') 'warning: card '//trim(input_line)//' ignored' + ! + endif + ! + ! ... end of loop ... ! + ! + goto 100 + ! +120 continue + ! + return + ! + end subroutine read_cards + + ! + ! ... description of the allowed input cards + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! atomic_species + ! + ! set the atomic species been read and their pseudopotential file + ! + ! syntax: + ! + ! atomic_specie + ! label(1) mass(1) psfile(1) + ! ... ... ... + ! label(n) mass(n) psfile(n) + ! + ! example: + ! + ! atomic_species + ! o 16.0 o.blyp.upf + ! h 1.00 h.fpmd.upf + ! + ! where: + ! + ! label(i) ( character(len=4) ) label of the atomic species + ! mass(i) ( real ) atomic mass + ! ( in u.m.a, carbon mass is 12.0 ) + ! psfile(i) ( character(len=80) ) file name of the pseudopotential + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_atomic_species( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: is, ip, ierr + character(len=4) :: lb_pos + character(len=256) :: psfile + ! + ! + if ( taspc ) then + call errore( ' card_atomic_species ', ' two occurrences', 2 ) + endif + if ( ntyp > nsx ) then + call errore( ' card_atomic_species ', ' nsp out of range ', ntyp ) + endif + ! + do is = 1, ntyp + ! + call read_line( input_line ) + read( input_line, *, iostat=ierr ) lb_pos, atom_mass(is), psfile + call errore( ' card_atomic_species ', & + 'cannot read atomic specie from: '//trim(input_line), abs(ierr)) + atom_pfile(is) = trim( psfile ) + lb_pos = adjustl( lb_pos ) + atom_label(is) = trim( lb_pos ) + ! + do ip = 1, is - 1 + if ( atom_label(ip) == atom_label(is) ) then + call errore( ' card_atomic_species ', & + & ' two occurrences of the same atomic label ', is ) + endif + enddo + ! + enddo + taspc = .true. + ! + return + ! + end subroutine card_atomic_species + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! atomic_positions + ! + ! set the atomic positions in the cell + ! + ! syntax: + ! + ! atomic_positions (units_option) + ! label(1) tau(1,1) tau(2,1) tau(3,1) mbl(1,1) mbl(2,1) mbl(3,1) + ! label(2) tau(1,2) tau(2,2) tau(3,2) mbl(1,2) mbl(2,2) mbl(3,2) + ! ... ... ... ... ... + ! label(n) tau(1,n) tau(2,n) tau(3,n) mbl(1,3) mbl(2,3) mbl(3,3) + ! + ! example: + ! + ! atomic_positions (bohr) + ! o 0.0099 0.0099 0.0000 0 0 0 + ! h 1.8325 -0.2243 -0.0001 1 1 1 + ! h -0.2243 1.8325 0.0002 1 1 1 + ! + ! where: + ! + ! units_option == crystal position are given in scaled units + ! units_option == bohr position are given in bohr + ! units_option == angstrom position are given in angstrom + ! units_option == alat position are given in units of alat + ! + ! label(k) ( character(len=4) ) atomic type + ! tau(:,k) ( real ) coordinates of the k-th atom + ! mbl(:,k) ( integer ) mbl(i,k) > 0 the i-th coord. of the + ! k-th atom is allowed to be moved + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_atomic_positions( input_line, prog ) + ! + use wrappers, only: feval_infix + ! + implicit none + ! + character(len=256) :: input_line + character(len=2) :: prog + character(len=4) :: lb_pos + integer :: ia, k, is, nfield, idx, rep_i + logical, external :: matches + logical :: tend + real(dp) :: inp(3) + integer :: fieldused + ! + integer :: ifield, ierr + real(dp) :: field_value + character(len=256) :: field_str, error_msg, wp + ! + ! + if ( tapos ) then + call errore( 'card_atomic_positions', 'two occurrences', 2 ) + endif + if ( .not. taspc ) then + call errore( 'card_atomic_positions', & + & 'atomic_species must be present before', 2 ) + endif + if ( ntyp > nsx ) then + call errore( 'card_atomic_positions', 'nsp out of range', ntyp ) + endif + if ( nat < 1 ) then + call errore( 'card_atomic_positions', 'nat out of range', nat ) + endif + ! + call allocate_input_ions(ntyp,nat) + ! + rd_if_pos = 1 + ! + sp_pos = 0 + rd_pos = 0.0_dp + na_inp = 0 + lsg=.false. + ! + if ( matches( "crystal_sg", input_line ) ) then + atomic_positions = 'crystal' + lsg=.true. + elseif ( matches( "crystal", input_line ) ) then + atomic_positions = 'crystal' + elseif ( matches( "bohr", input_line ) ) then + atomic_positions = 'bohr' + elseif ( matches( "angstrom", input_line ) ) then + atomic_positions = 'angstrom' + elseif ( matches( "alat", input_line ) ) then + atomic_positions = 'alat' + else + if ( trim( adjustl( input_line ) ) /= 'atomic_positions' ) then + call errore( 'read_cards ', & + & 'unknown option for atomic_position: '& + & // input_line, 1 ) + endif + call infomsg( 'read_cards ', & + & 'deprecated: no units specified in atomic_positions card' ) + if ( prog == 'cp' ) atomic_positions = 'bohr' + if ( prog == 'pw' ) atomic_positions = 'alat' + call infomsg( 'read_cards ', & + & 'atomic_positions: units set to '//trim(atomic_positions) ) + endif + ! + reader_loop : do ia = 1,nat + ! + call read_line( input_line, end_of_file = tend ) + if ( tend ) call errore( 'read_cards', & + 'end of file reading atomic positions', ia ) + ! + call field_count( nfield, input_line ) + ! + ! read atom symbol (column 1) + ! + call get_field(1, lb_pos, input_line) + lb_pos = trim(lb_pos) + ! + error_msg = 'error while parsing atomic position card.' + ! + ! read field 2 (atom x coordinate or wyckoff position symbol) + ! + call get_field(2, field_str, input_line) + ! + ! check if position ia is expressed in wyckoff positions + ! + idx = len_trim(field_str) + if ( lsg .and. (idx < 4) .and. & + ( iachar(field_str(idx:idx)) > 64 .and. & + iachar(field_str(idx:idx)) < 123 ) ) then + ! + ! wyckoff positions + ! + if ( nfield < 3 .and. nfield > 8 ) & + call errore( 'read_cards', 'wrong number of columns ' // & + & 'in atomic_positions', ia ) + wp=field_str + inp(:)=1.d5 + ! + do k = 3,min(nfield,5) + ! read k-th field (coordinate k-2) + call get_field(k, field_str, input_line) + inp(k-2) = feval_infix(ierr, field_str ) + call errore('card_atomic_positions', error_msg, ierr) + enddo + ! + call wypos(rd_pos(1,ia),wp,inp,space_group, & + uniqueb,rhombohedral,origin_choice) + ! + ! count how many fields were used to find wyckoff positions + ! + fieldused=2 + if ( any (rd_pos(1:3,ia)==inp(1)) ) fieldused=fieldused+1 + if ( any (rd_pos(2:3,ia)==inp(2)) ) fieldused=fieldused+1 + if ( rd_pos( 3,ia)==inp(3) ) fieldused=fieldused+1 + ! + else + ! + ! no wyckoff positions + ! + if ( nfield /= 4 .and. nfield /= 7 ) & + call errore( 'read_cards', 'wrong number of columns ' // & + & 'in atomic_positions', ia ) + ! + ! field just read is coordinate x + ! + rd_pos(1,ia) = feval_infix(ierr, field_str ) + call errore('card_atomic_positions', error_msg, ierr) + do k = 3,4 + ! read fields 3 and 4 (atom y and z coordinate) + call get_field(k, field_str, input_line) + rd_pos(k-1,ia) = feval_infix(ierr, field_str ) + call errore('card_atomic_positions', error_msg, ierr) + end do + ! + fieldused=4 + ! + endif + ! read constraints if present (last 3 fields) + if ( nfield-fieldused > 0 .and. nfield-fieldused /= 3 ) & + call errore( 'read_cards', 'unexpected number of columns ' // & + & 'in atomic_positions', ia ) + do k = fieldused+1, nfield + call get_field(k, field_str, input_line) + read(field_str, *) rd_if_pos(k-fieldused,ia) + enddo + ! + match_label: do is = 1, ntyp + ! + if ( trim(lb_pos) == trim( atom_label(is) ) ) then + ! + sp_pos(ia) = is + exit match_label + ! + endif + ! + enddo match_label + ! + if( ( sp_pos(ia) < 1 ) .or. ( sp_pos(ia) > ntyp ) ) then + ! + call errore( 'read_cards', 'species '//trim(lb_pos)// & + & ' in atomic_positions is nonexistent', ia ) + ! + endif + ! + is = sp_pos(ia) + ! + na_inp(is) = na_inp(is) + 1 + ! + enddo reader_loop + ! + tapos = .true. + ! + + return + ! + end subroutine card_atomic_positions + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! atomic_forces + ! + ! read external forces (in atomic units) from standard input + ! + ! syntax: + ! + ! atomic_forces + ! label fx(1) fy(1) fz(1) + ! ..... + ! label fx(n) fy(n) fz(n) + ! + ! example: + ! + ! ??? + ! + ! where: + ! + ! label (character(len=4)) atomic label + ! fx(:), fy(:) and fz(:) (real) x, y and z component of the external force + ! acting on the ions whose coordinate are given + ! in the same line in card atomic_position + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_atomic_forces( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: ia, k, nfield + character(len=4) :: lb + ! + ! + if( tforces ) then + call errore( ' card_atomic_forces ', ' two occurrences ', 2 ) + endif + ! + if( .not. tapos ) then + call errore( ' card_atomic_forces ', & + & ' atomic_species must be present before ', 2 ) + endif + ! + rd_for = 0.0_dp + ! + do ia = 1, nat + ! + call read_line( input_line ) + call field_count( nfield, input_line ) + if ( nfield == 4 ) then + read(input_line,*) lb, ( rd_for(k,ia), k = 1, 3 ) + elseif( nfield == 3 ) then + read(input_line,*) ( rd_for(k,ia), k = 1, 3 ) + else + call errore( ' iosys ', ' wrong entries in atomic_forces ', ia ) + endif + ! + enddo + ! + tforces = .true. + ! + return + ! + end subroutine card_atomic_forces + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! k_points + ! + ! use the specified set of k points + ! + ! syntax: + ! + ! k_points (mesh_option) + ! n + ! xk(1,1) xk(2,1) xk(3,1) wk(1) + ! ... ... ... ... + ! xk(1,n) xk(2,n) xk(3,n) wk(n) + ! + ! example: + ! + ! k_points + ! 10 + ! 0.1250000 0.1250000 0.1250000 1.00 + ! 0.1250000 0.1250000 0.3750000 3.00 + ! 0.1250000 0.1250000 0.6250000 3.00 + ! 0.1250000 0.1250000 0.8750000 3.00 + ! 0.1250000 0.3750000 0.3750000 3.00 + ! 0.1250000 0.3750000 0.6250000 6.00 + ! 0.1250000 0.3750000 0.8750000 6.00 + ! 0.1250000 0.6250000 0.6250000 3.00 + ! 0.3750000 0.3750000 0.3750000 1.00 + ! 0.3750000 0.3750000 0.6250000 3.00 + ! + ! where: + ! + ! mesh_option == automatic k points mesh is generated automatically + ! with monkhorst-pack algorithm + ! mesh_option == crystal k points mesh is given in stdin in scaled + ! units + ! mesh_option == tpiba k points mesh is given in stdin in units + ! of ( 2 pi / alat ) + ! mesh_option == gamma only gamma point is used ( default in + ! cpmd simulation ) + ! mesh_option == tpiba_b as tpiba but the weights gives the + ! number of points between this point + ! and the next + ! mesh_option == crystal_b as crystal but the weights gives the + ! number of points between this point and + ! the next + ! mesh_option == tpiba_c the code expects three k points + ! k_0, k_1, k_2 in tpiba units. + ! these points define a rectangle + ! in reciprocal space with vertices k_0, k_1, + ! k_2, k_1+k_2-k_0: k_0 + \alpha (k_1-k_0)+ + ! \beta (k_2-k_0) with 0<\alpha,\beta < 1. + ! the code produces a uniform mesh n1 x n2 + ! k points in this rectangle. n1 and n2 are + ! the weights of k_1 and k_2. the weight of k_0 + ! is not used. useful for contour plots of the + ! bands. + ! mesh_option == crystal_c as tpiba_c but the k points are given + ! in crystal coordinates. + ! + ! + ! n ( integer ) number of k points + ! xk(:,i) ( real ) coordinates of i-th k point + ! wk(i) ( real ) weights of i-th k point + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_kpoints( input_line ) + ! + use bz_form, only : transform_label_coord + use cell_base, only : cell_base_init, celldm_cb => celldm + implicit none + ! + character(len=256) :: input_line, buffer + integer :: i, j + integer :: nkaux + integer, allocatable :: wkaux(:) + real(dp), allocatable :: xkaux(:,:) + integer :: npk_label, nch + character(len=3), allocatable :: letter(:) + integer, allocatable :: label_list(:) + real(dp) :: delta, wk0 + real(dp) :: dkx(3), dky(3) + logical, external :: matches + logical :: tend,terr + logical :: kband = .false. + logical :: kband_plane = .false. + ! + ! + if ( tkpoints ) then + call errore( ' card_kpoints ', ' two occurrences', 2 ) + endif + ! + if ( matches( "automatic", input_line ) ) then + ! automatic generation of k-points + k_points = 'automatic' + elseif ( matches( "crystal", input_line ) ) then + ! input k-points are in crystal (reciprocal lattice) axis + k_points = 'crystal' + if ( matches( "_b", input_line ) ) kband=.true. + if ( matches( "_c", input_line ) ) kband_plane=.true. + elseif ( matches( "tpiba", input_line ) ) then + ! input k-points are in 2pi/a units + k_points = 'tpiba' + if ( matches( "_b", input_line ) ) kband=.true. + if ( matches( "_c", input_line ) ) kband_plane=.true. + elseif ( matches( "gamma", input_line ) ) then + ! only gamma (k=0) is used + k_points = 'gamma' + else + ! by default, input k-points are in 2pi/a units + k_points = 'tpiba' + endif + ! + if ( k_points == 'automatic' ) then + ! + ! ... automatic generation of k-points + ! + nkstot = 0 + call read_line( input_line, end_of_file = tend, error = terr ) + if (tend) goto 10 + if (terr) goto 20 + read(input_line, *, end=10, err=20) nk1, nk2, nk3, k1, k2 ,k3 + if ( k1 < 0 .or. k1 > 1 .or. & + k2 < 0 .or. k2 > 1 .or. & + k3 < 0 .or. k3 > 1 ) call errore & + ('card_kpoints', 'invalid offsets: must be 0 or 1', 1) + if ( nk1 <= 0 .or. nk2 <= 0 .or. nk3 <= 0 ) call errore & + ('card_kpoints', 'invalid values for nk1, nk2, nk3', 1) + allocate ( xk(3,1), wk(1) ) ! prevents problems with debug flags + ! ! when init_startk is called in iosys + elseif ( ( k_points == 'tpiba' ) .or. ( k_points == 'crystal' ) ) then + ! + ! ... input k-points + ! + call read_line( input_line, end_of_file = tend, error = terr ) + if (tend) goto 10 + if (terr) goto 20 + read(input_line, *, end=10, err=20) nkstot + ! + if (kband) then +! +! only the initial and final k points of the lines are given +! + nkaux=nkstot + allocate(xkaux(3,nkstot), wkaux(nkstot)) + allocate ( letter(nkstot) ) + allocate ( label_list(nkstot) ) + npk_label=0 + do i = 1, nkstot + call read_line( input_line, end_of_file = tend, error = terr ) + if (tend) goto 10 + if (terr) goto 20 + do j=1,256 ! loop over all characters of input_line + if ((ichar(input_line(j:j)) < 58 .and. & ! a digit + ichar(input_line(j:j)) > 47) & + .or. ichar(input_line(j:j)) == 43 .or. & ! the + sign + ichar(input_line(j:j))== 45 .or. & ! the - sign + ichar(input_line(j:j))== 46 ) then ! a dot . +! +! this is a digit, therefore this line contains the coordinates of the +! k point. we read it and exit from the loop on the characters +! + read(input_line,*, end=10, err=20) xkaux(1,i), & + xkaux(2,i), xkaux(3,i), wk0 + wkaux(i) = nint ( wk0 ) ! beware: wkaux is integer + exit + elseif ((ichar(input_line(j:j)) < 123 .and. & + ichar(input_line(j:j)) > 64)) then +! +! this is a letter, not a space character. we read the next three +! characters and save them in the letter array, save also which k point +! it is +! + npk_label=npk_label+1 + read(input_line(j:),'(a3)') letter(npk_label) + label_list(npk_label)=i +! +! now we remove the letters from input_line and read the number of points +! of the line. the next two line should account for the case in which +! there is only one space between the letter and the number of points. +! + nch=3 + if ( ichar(input_line(j+1:j+1))==32 .or. & + ichar(input_line(j+2:j+2))==32 ) nch=2 + buffer=input_line(j+nch:) + read(buffer,*,err=20) wkaux(i) + exit + endif + enddo + enddo + if ( npk_label > 0 ) then + call cell_base_init ( ibrav, celldm, a, b, c, cosab, & + cosac, cosbc, trd_ht, rd_ht, cell_units ) + call transform_label_coord(ibrav, celldm_cb, xkaux, letter, & + label_list, npk_label, nkstot, k_points, point_label_type ) + end if + + deallocate(letter) + deallocate(label_list) + ! count k-points first + nkstot=sum(wkaux(1:nkaux-1))+1 + do i=1,nkaux-1 + if (wkaux(i)==0) nkstot=nkstot+1 + enddo + allocate ( xk(3,nkstot), wk(nkstot) ) + ! + ! generate the points along the lines + ! + call generate_k_along_lines(nkaux, xkaux, wkaux, xk, wk, nkstot) + ! + ! workaround: discard current wk (contains the length of k-path, + ! never used), replace with wk=1 so that band occupations (wg) + ! are correctly written to file - needed by berkeleygw interface + ! + wk(:) = 1.0_dp + deallocate(xkaux) + deallocate(wkaux) + ! + elseif (kband_plane) then +! +! generate a uniform mesh of k points on the plane defined by +! the origin k_0, and two vectors applied in k_0, k_1 and k_2. +! + if (nkstot /= 3) call errore ('card_kpoints', & + 'option _c requires 3 k points',i) + nkaux=nkstot + allocate(xkaux(3,nkstot), wkaux(nkstot)) + do i = 1, nkstot + call read_line( input_line, end_of_file = tend, error = terr ) + if (tend) goto 10 + if (terr) goto 20 + read(input_line,*, end=10, err=20) xkaux(1,i), xkaux(2,i), & + xkaux(3,i), wk0 + wkaux(i) = nint ( wk0 ) ! beware: wkaux is integer + enddo + ! count k-points first + nkstot = wkaux(2) * wkaux(3) + allocate ( xk(3,nkstot), wk(nkstot) ) + call generate_k_in_plane(nkaux, xkaux, wkaux, xk, wk, nkstot) + deallocate(xkaux) + deallocate(wkaux) + else +! +! reads on input the k points +! + allocate ( xk(3, nkstot), wk(nkstot) ) + do i = 1, nkstot + call read_line( input_line, end_of_file = tend, error = terr ) + if (tend) goto 10 + if (terr) goto 20 + read(input_line,*, end=10, err=20) xk(1,i),xk(2,i),xk(3,i),wk(i) + enddo + endif + ! + elseif ( k_points == 'gamma' ) then + ! + nkstot = 1 + allocate ( xk(3,1), wk(1) ) + xk(:,1) = 0.0_dp + wk(1) = 1.0_dp + ! + endif + ! + tkpoints = .true. + tk_inp = .true. + ! + return +10 call errore ('card_kpoints', ' end of file while reading ' & + & // trim(k_points) // ' k points', 1) +20 call errore ('card_kpoints', ' error while reading ' & + & // trim(k_points) // ' k points', 1) + ! + end subroutine card_kpoints + + subroutine card_add_kpoints( input_line ) + use additional_kpoints, only : nkstot_add, xk_add + implicit none + character(len=*),intent(in) :: input_line + character(len=256) :: input_line_aux + real(dp),allocatable :: xk_old(:,:), wk_old(:) + integer :: nk1_old, nk2_old, nk3_old, nkstot_old + integer :: k1_old, k2_old, k3_old + logical, external :: matches + ! + if(.not.allocated(xk) .or. .not.allocated(wk))& + call errore("add_kpoints", "additional_k_points must appear after k_points",1) + if(.not.tkpoints) & + call errore("add_kpoints", "additional_k_points must appear after k_points",2) + if(matches( "automatic", input_line )) & + call errore("add_kpoints", "additional_k_points cannot be 'automatic'", 3) + + ! back-up existing points + nkstot_old = nkstot + allocate(xk_old(3,nkstot_old)) + allocate(wk_old(nkstot_old)) + xk_old = xk + wk_old = wk + nk1_old = nk1 + nk2_old = nk2 + nk3_old = nk3 + k1_old = k1 + k2_old = k2 + k3_old = k3 + deallocate(xk,wk) + + ! prepare to read k-points again + nkstot = 0 + input_line_aux = trim(adjustl(input_line)) + input_line_aux = input_line_aux(12:) + tkpoints = .false. + call card_kpoints(input_line_aux) + ! + ! backup new points to module + nkstot_add = nkstot + if(nkstot_add==0) call errore("add_kpoints", "no new k_points?",1) + allocate(xk_add(3,nkstot_add)) + xk_add = xk + + ! put back previous stuff + deallocate(xk, wk) + nkstot = nkstot_old + allocate(xk(3,nkstot)) + allocate(wk(nkstot)) + xk = xk_old + wk = wk_old + nk1 = nk1_old + nk2 = nk2_old + nk3 = nk3_old + k1 = k1_old + k2 = k2_old + k3 = k3_old + deallocate(xk_old,wk_old) + + return + end subroutine card_add_kpoints + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! occupations + ! + ! use the specified occupation numbers for electronic states. + ! note that you should specify 10 values per line maximum! + ! + ! syntax (nspin == 1): + ! + ! occupations + ! f(1) .... .... f(10) + ! f(11) .... f(nbnd) + ! + ! syntax (nspin == 2): + ! + ! occupations + ! u(1) .... .... u(10) + ! u(11) .... u(nbnd) + ! d(1) .... .... d(10) + ! d(11) .... d(nbnd) + ! + ! example: + ! + ! occupations + ! 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 + ! 2.0 2.0 2.0 2.0 2.0 1.0 1.0 + ! + ! where: + ! + ! f(:) (real) these are the occupation numbers + ! for lda electronic states. + ! + ! u(:) (real) these are the occupation numbers + ! for lsd spin == 1 electronic states + ! d(:) (real) these are the occupation numbers + ! for lsd spin == 2 electronic states + ! + ! note, maximum 10 values per line! + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_occupations( input_line ) + ! + use wrappers, only: feval_infix + ! + implicit none + ! + character(len=256) :: input_line, field_str + integer :: is, nx10, i, j, nspin0 + integer :: nfield, nbnd_read, nf, ierr + logical :: tef + ! + ! + if ( tocc ) then + call errore( ' card_occupations ', ' two occurrences', 2 ) + endif + nspin0=nspin + if (nspin == 4) nspin0=1 + ! + allocate ( f_inp ( nbnd, nspin0 ) ) + do is = 1, nspin0 + ! + nbnd_read = 0 + do while ( nbnd_read < nbnd) + call read_line( input_line, end_of_file=tef ) + if (tef) call errore('card_occupations',& + 'missing occupations, end of file reached',1) + call field_count( nfield, input_line ) + ! + do nf = 1,nfield + nbnd_read = nbnd_read+1 + if (nbnd_read > nbnd ) exit + call get_field(nf, field_str, input_line) + ! + f_inp(nbnd_read,is) = feval_infix(ierr, field_str ) + call errore('card_occupations',& + 'error parsing occupation: '//trim(field_str), nbnd_read*ierr) + enddo + enddo + ! + enddo + ! + tf_inp = .true. + tocc = .true. + ! + return + ! + end subroutine card_occupations + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! cell_parameters + ! + ! use the specified cell dimensions + ! + ! syntax: + ! + ! cell_parameters (cell_option) + ! ht(1,1) ht(1,2) ht(1,3) + ! ht(2,1) ht(2,2) ht(2,3) + ! ht(3,1) ht(3,2) ht(3,3) + ! + ! cell_option == alat lattice vectors in units of alat + ! cell_option == bohr lattice vectors in bohr + ! cell_option == angstrom lattice vectors in angstrom + ! + ! example: + ! + ! cell_parameters + ! 24.50644311 0.00004215 -0.14717844 + ! -0.00211522 8.12850030 1.70624903 + ! 0.16447787 0.74511792 23.07395418 + ! + ! where: + ! + ! ht(i,j) (real) cell dimensions ( in a.u. ), + ! note the relation with lattice vectors: + ! ht(1,:) = a1, ht(2,:) = a2, ht(3,:) = a3 + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_cell_parameters( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: i, j + logical, external :: matches + ! + ! + if ( tcell ) then + call errore( ' card_cell_parameters ', ' two occurrences', 2 ) + endif + ! + if ( matches( "bohr", input_line ) ) then + cell_units = 'bohr' + elseif ( matches( "angstrom", input_line ) ) then + cell_units = 'angstrom' + elseif ( matches( "alat", input_line ) ) then + cell_units = 'alat' + else + cell_units = 'none' + call infomsg( 'read_cards ', & + & 'deprecated: no units specified in cell_parameters card' ) + ! cell parameters are set in cell_base_init + endif + ! + do i = 1, 3 + call read_line( input_line ) + read(input_line,*) ( rd_ht( i, j ), j = 1, 3 ) + enddo + ! + trd_ht = .true. + tcell = .true. + ! + return + ! + end subroutine card_cell_parameters + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! ref_cell_parameters + ! + ! use the specified cell dimensions + ! + ! syntax: + ! + ! ref_cell_parameters (cell_option) + ! rd_ref_ht(1,1) rd_ref_ht(1,2) rd_ref_ht(1,3) + ! rd_ref_ht(2,1) rd_ref_ht(2,2) rd_ref_ht(2,3) + ! rd_ref_ht(3,1) rd_ref_ht(3,2) rd_ref_ht(3,3) + ! + ! cell_option == alat lattice vectors in units of alat set by ref_alat keyword (default) + ! cell_option == bohr lattice vectors in bohr + ! cell_option == angstrom lattice vectors in angstrom + ! + ! example: + ! + ! ref_cell_parameters + ! 24.50644311 0.00004215 -0.14717844 + ! -0.00211522 8.12850030 1.70624903 + ! 0.16447787 0.74511792 23.07395418 + ! + ! where: + ! + ! rd_ref_ht(i,j) (real) cell dimensions ( in a.u. ), + ! note the relation with reference lattice vectors: + ! rd_ref_ht(1,:) = ref_a1, rd_ref_ht(2,:) = ref_a2, rd_ref_ht(3,:) = re_a3 + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_ref_cell_parameters( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: i, j + logical, external :: matches + ! + ! + if ( ref_cell ) then + call errore( ' card_reference_cell_parameters ', ' two occurrences', 2 ) + endif + ! + if ( matches( "bohr", input_line ) ) then + ref_cell_units = 'bohr' + elseif ( matches( "angstrom", input_line ) ) then + ref_cell_units = 'angstrom' + else + ref_cell_units = 'alat' + endif + ! + do i = 1, 3 + call read_line( input_line ) + read(input_line,*) ( rd_ref_ht( i, j ), j = 1, 3 ) + enddo + ! + ref_cell = .true. + ! + return + ! + end subroutine card_ref_cell_parameters + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! atomic_velocities + ! + ! read velocities (in atomic units) from standard input + ! + ! syntax: + ! + ! atomic_velocities + ! label(1) vx(1) vy(1) vz(1) + ! .... + ! label(n) vx(n) vy(n) vz(n) + ! + ! example: + ! + ! ??? + ! + ! where: + ! + ! label (character(len=4)) atomic label + ! vx(:), vy(:) and vz(:) (real) x, y and z velocity components of + ! the ions + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_ion_velocities( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: ia, k, is, nfield + character(len=4) :: lb_vel + ! + ! + if( tionvel ) then + call errore( ' card_ion_velocities ', ' two occurrences', 2 ) + endif + ! + if( .not. tapos ) then + call errore( ' card_ion_velocities ', & + & ' atomic_species must be present before ', 2 ) + endif + ! + rd_vel = 0.0_dp + sp_vel = 0 + ! + if ( ion_velocities == 'from_input' ) then + ! + tavel = .true. + ! + do ia = 1, nat + ! + call read_line( input_line ) + call field_count( nfield, input_line ) + if ( nfield == 4 ) then + read(input_line,*) lb_vel, ( rd_vel(k,ia), k = 1, 3 ) + else + call errore( ' iosys ', & + & ' wrong entries in atomic_velocities ', ia ) + endif + ! + match_label: do is = 1, ntyp + if ( trim( lb_vel ) == atom_label(is) ) then + sp_vel(ia) = is + exit match_label + endif + enddo match_label + ! + if ( sp_vel(ia) < 1 .or. sp_vel(ia) > ntyp ) then + call errore( ' iosys ', ' wrong label in ion_velocities ', ia ) + endif + ! + enddo + ! + endif + ! + tionvel = .true. + ! + return + ! + end subroutine + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! constraints + ! + ! ionic constraints + ! + ! syntax: + ! + ! constraints + ! nconstr constr_tol + ! constr_type(.) constr(1,.) constr(2,.) ... { constr_target(.) } + ! + ! where: + ! + ! nconstr(integer) number of constraints + ! + ! constr_tol tolerance for keeping the constraints + ! satisfied + ! + ! constr_type(.) type of constrain: + ! 1: for fixed distances ( two atom indexes must + ! be specified ) + ! 2: for fixed planar angles ( three atom indexes + ! must be specified ) + ! + ! constr(1,.) constr(2,.) ... + ! + ! indices object of the constraint, as + ! they appear in the 'position' card + ! + ! constr_target target for the constrain ( in the case of + ! planar angles it is the cos of the angle ). + ! this variable is optional. + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_constraints( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: i, nfield + ! + ! + if ( tconstr ) call errore( 'card_constraints', 'two occurrences', 2 ) + ! + call read_line( input_line ) + ! + call field_count( nfield, input_line ) + ! + if ( nfield == 1 ) then + ! + read( input_line, * ) nconstr_inp + ! + elseif ( nfield == 2 ) then + ! + read( input_line, * ) nconstr_inp, constr_tol_inp + ! + else + ! + call errore( 'card_constraints', 'too many fields', nfield ) + ! + endif + write(stdout,'(5x,a,i4,a,f12.6)') & + 'reading',nconstr_inp,' constraints; tolerance:', constr_tol_inp + ! + call allocate_input_constr() + ! + do i = 1, nconstr_inp + ! + call read_line( input_line ) + ! + read( input_line, * ) constr_type_inp(i) + ! + call field_count( nfield, input_line ) + ! + if ( nfield > nc_fields + 2 ) & + call errore( 'card_constraints', & + 'too many fields for this constraint', i ) + ! + select case( constr_type_inp(i) ) + case( 'type_coord', 'atom_coord' ) + ! + if ( nfield == 5 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_inp(4,i) + ! + write(stdout,'(7x,i3,a,i3,a,i2,a,2f12.6)') i, & + ') '//constr_type_inp(i)(1:4),int(constr_inp(1,i)) ,& + ' coordination wrt type:', int(constr_inp(2,i)), & + ' cutoff distance and smoothing:', constr_inp(3:4,i) + elseif ( nfield == 6 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_inp(4,i), & + constr_target_inp(i) + ! + constr_target_set(i) = .true. + ! + write(stdout,'(7x,i3,a,i3,a,i2,a,2f12.6,a,f12.6)') i, & + ') '//constr_type_inp(i)(1:4),int(constr_inp(1,i)) , & + ' coordination wrt type:', int(constr_inp(2,i)), & + ' cutoff distance and smoothing:', constr_inp(3:4,i), & + '; target:', constr_target_inp(i) + else + ! + call errore( 'card_constraints', 'type_coord, ' // & + & 'atom_coord: wrong number of fields', nfield ) + ! + endif + ! + case( 'distance' ) + ! + if ( nfield == 3 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i) + ! + write(stdout,'(7x,i3,a,2i3)') & + i,') distance between atoms: ', int(constr_inp(1:2,i)) + elseif ( nfield == 4 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_target_inp(i) + ! + constr_target_set(i) = .true. + ! + write(stdout,'(7x,i3,a,2i3,a,f12.6)') i, & + ') distance between atoms: ', int(constr_inp(1:2,i)), & + '; target:', constr_target_inp(i) + else + ! + call errore( 'card_constraints', & + & 'distance: wrong number of fields', nfield ) + ! + endif + ! + case( 'planar_angle' ) + ! + if ( nfield == 4 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i) + ! + write(stdout, '(7x,i3,a,3i3)') & + i,') planar angle between atoms: ', int(constr_inp(1:3,i)) + elseif ( nfield == 5 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_target_inp(i) + ! + constr_target_set(i) = .true. + ! + write(stdout, '(7x,i3,a,3i3,a,f12.6)') i, & + ') planar angle between atoms: ', int(constr_inp(1:3,i)), & + '; target:', constr_target_inp(i) + else + ! + call errore( 'card_constraints', & + & 'planar_angle: wrong number of fields', nfield ) + ! + endif + ! + case( 'torsional_angle' ) + ! + if ( nfield == 5 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_inp(4,i) + ! + write(stdout, '(7x,i3,a,4i3)') & + i,') torsional angle between atoms: ', int(constr_inp(1:4,i)) + elseif ( nfield == 6 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_inp(4,i), & + constr_target_inp(i) + ! + constr_target_set(i) = .true. + ! + write(stdout, '(7x,i3,a,4i3,a,f12.6)') i, & + ') torsional angle between atoms: ', int(constr_inp(1:4,i)),& + '; target:', constr_target_inp(i) + else + ! + call errore( 'card_constraints', & + & 'torsional_angle: wrong number of fields', nfield ) + ! + endif + ! + case( 'bennett_proj' ) + ! + if ( nfield == 5 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_inp(4,i) + ! + write(stdout, '(7x,i3,a,i3,a,3f12.6)') i, & + ') bennet projection of atom ', int(constr_inp(1,i)), & + ' along vector:', constr_inp(2:4,i) + elseif ( nfield == 6 ) then + ! + read( input_line, * ) constr_type_inp(i), & + constr_inp(1,i), & + constr_inp(2,i), & + constr_inp(3,i), & + constr_inp(4,i), & + constr_target_inp(i) + ! + constr_target_set(i) = .true. + ! + write(stdout, '(7x,i3,a,i3,a,3f12.6,a,f12.6)') i, & + ') bennet projection of atom ', int(constr_inp(1,i)), & + ' along vector:', constr_inp(2:4,i), & + '; target:', constr_target_inp(i) + else + ! + call errore( 'card_constraints', & + & 'bennett_proj: wrong number of fields', nfield ) + ! + endif + ! + case default + ! + call errore( 'card_constraints', 'unknown constraint ' // & + & 'type: ' // trim( constr_type_inp(i) ), 1 ) + ! + end select + ! + enddo + ! + tconstr = .true. + ! + return + ! + end subroutine card_constraints + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! ksout + ! + ! enable the printing of kohn sham states + ! + ! syntax ( nspin == 2 ): + ! + ! ksout + ! nu + ! iu(1) iu(2) iu(3) .. iu(nu) + ! nd + ! id(1) id(2) id(3) .. id(nd) + ! + ! syntax ( nspin == 1 ): + ! + ! ksout + ! ns + ! is(1) is(2) is(3) .. is(ns) + ! + ! example: + ! + ! ??? + ! + ! where: + ! + ! nu (integer) number of spin=1 states to be printed + ! iu(:) (integer) indexes of spin=1 states, the state iu(k) + ! is saved to file ks_up.iu(k) + ! + ! nd (integer) number of spin=2 states to be printed + ! id(:) (integer) indexes of spin=2 states, the state id(k) + ! is saved to file ks_dw.id(k) + ! + ! ns (integer) number of lda states to be printed + ! is(:) (integer) indexes of lda states, the state is(k) + ! is saved to file ks.is(k) + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_ksout( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + integer :: i, s, nksx + type occupancy_type + integer, pointer :: occs(:) + end type occupancy_type + type(occupancy_type), allocatable :: is(:) + ! + if ( tksout ) then + call errore( ' card_ksout ', ' two occurrences', 2 ) + endif + ! + nprnks = 0 + nksx = 0 + ! + allocate ( is (nspin) ) + ! + do s = 1, nspin + ! + call read_line( input_line ) + read(input_line, *) nprnks( s ) + ! + if ( nprnks( s ) < 1 ) then + call errore( ' card_ksout ', ' wrong number of states ', 2 ) + endif + ! + allocate( is(s)%occs( 1:nprnks(s) ) ) + ! + call read_line( input_line ) + read(input_line, *) ( is(s)%occs(i), i = 1, nprnks( s ) ) + ! + nksx = max( nksx, nprnks( s ) ) + ! + enddo + ! + call allocate_input_iprnks( nksx, nspin ) + ! + do s = 1, nspin + ! + do i = 1, nprnks( s ) + ! + iprnks( i, s ) = is(s)%occs(i) + ! + enddo + ! + deallocate( is(s)%occs ) + ! + enddo + ! + deallocate( is ) + ! + tksout = .true. + ! + return + ! + end subroutine + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! plot wannier + ! + ! needed to specify the indices of the wannier functions that + ! have to be plotted + ! + ! syntax: + ! + ! plot_wannier + ! index1, ..., indexn + ! + ! where: + ! + ! index1, ..., indexn are indices of the wannier functions + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_plot_wannier( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + logical, external :: matches + ! + integer :: i, ib + character(len=6) :: i_char + character(len=6), external :: int_to_char + ! + ! + if ( twannier ) & + call errore( 'card_plot_wannier', 'two occurrences', 2 ) + ! + if ( nwf > 0 ) then + ! + if ( nwf > nwf_max ) & + call errore( 'card_plot_wannier', 'too many wannier functions', 1 ) + ! + call read_line( input_line ) + ! + ib = 0 + ! + do i = 1, nwf_max + ! + i_char = int_to_char( i ) + ! + if ( matches( ' ' // trim( i_char ) // ',', & + ' ' // trim( input_line ) // ',' ) ) then + ! + ib = ib + 1 + ! + if ( ib > nwf ) & + call errore( 'card_plot_wannier', 'too many indices', 1 ) + ! + wannier_index(ib) = i + ! + endif + ! + enddo + ! + endif + ! + twannier = .true. + ! + return + ! + end subroutine card_plot_wannier + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + ! + ! + ! template + ! + ! this is a template card info section + ! + ! syntax: + ! + ! template + ! rvalue ivalue + ! + ! example: + ! + ! ??? + ! + ! where: + ! + ! rvalue (real) this is a real value + ! ivalue (integer) this is an integer value + ! + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_template( input_line ) + ! + implicit none + ! + character(len=256) :: input_line + ! + ! + if ( ttemplate ) then + call errore( ' card_template ', ' two occurrences', 2 ) + endif + ! + ! .... code here + ! + ttemplate = .true. + ! + return + ! + end subroutine + ! + ! + !------------------------------------------------------------------------ + ! begin manual + !---------------------------------------------------------------------- + !wannier_ac + !wannier# 1 10.5 15.7 2 + !atom 1 + !d 1 0.45 + !p 3 0.55 + !wannier# 2 10.5 15.7 1 + !atom 3 + !p 1 0.8 + !spin#2: + !wannier# 1 10.5 15.7 2 + !atom 1 + !d 1 0.45 + !p 3 0.55 + !wannier# 2 10.5 15.7 1 + !atom 3 + !p 1 0.8 + !---------------------------------------------------------------------- + ! end manual + !------------------------------------------------------------------------ + ! + subroutine card_wannier_ac( input_line ) + ! + use wannier_new, only: nwan + + implicit none + ! + character(len=256) :: input_line + integer :: i,j,k, nfield, iwan, ning, iatom,il,im,ispin + logical :: tend + real :: c, b_from, b_to + character(len=10) :: text, lo + + ispin = 1 + ! + do i = 1, nwan + ! + call read_line( input_line, end_of_file = tend ) + ! + if ( tend ) & + call errore( 'read_cards', & + 'end of file reading trial wfc composition', i ) + ! + call field_count( nfield, input_line ) + ! + if ( nfield == 4 ) then + read(input_line,*) text, iwan, b_from, b_to + ning = 1 + elseif ( nfield == 5 ) then + read(input_line,*) text, iwan, b_from, b_to, ning + else + call errore( 'read_cards', & + 'wrong format', nfield ) + endif + if(iwan/=i) call errore( 'read_cards', 'wrong wannier order', iwan) + + ! read atom number + call read_line( input_line, end_of_file = tend ) + read(input_line,*) text, iatom + ! + wan_data(iwan,ispin)%iatom = iatom + wan_data(iwan,ispin)%ning = ning + wan_data(iwan,ispin)%bands_from = b_from + wan_data(iwan,ispin)%bands_to = b_to + ! + do j=1, ning + call read_line( input_line, end_of_file = tend ) + ! + if ( tend ) & + call errore( 'read_cards', & + 'not enough wavefunctions', j ) + if (ning==1) then + read(input_line,*) lo,im + c = 1.d0 + else + read(input_line,*) lo,im,c + endif + + select case(trim(lo)) + case('s') + il = 0 + case('p') + il = 1 + case('d') + il = 2 + case('f') + il = 3 + case default + call errore( 'read_cards', & + 'wrong l-label', 1 ) + end select + + wan_data(iwan,ispin)%ing(j)%l = il + wan_data(iwan,ispin)%ing(j)%m = im + wan_data(iwan,ispin)%ing(j)%c = c + enddo + enddo + + !is there spin 2 information? + call read_line( input_line, end_of_file = tend ) + ! + if ( .not. tend ) then + read(input_line,*) text + if ( trim(text) == 'spin#2:') then ! ok, there is spin 2 data + ispin = 2 + ! + do i = 1, nwan + ! + call read_line( input_line, end_of_file = tend ) + ! + if ( tend ) & + call errore( 'read_cards', & + 'end of file reading trial wfc composition', i ) + ! + call field_count( nfield, input_line ) + ! + if ( nfield == 4 ) then + read(input_line,*) text, iwan, b_from, b_to + ning = 1 + elseif ( nfield == 5 ) then + read(input_line,*) text, iwan, b_from, b_to, ning + else + call errore( 'read_cards', & + 'wrong format', nfield ) + endif + if(iwan/=i) call errore( 'read_cards', 'wrong wannier order', iwan) + + ! read atom number + call read_line( input_line, end_of_file = tend ) + read(input_line,*) text, iatom + ! + wan_data(iwan,ispin)%iatom = iatom + wan_data(iwan,ispin)%ning = ning + wan_data(iwan,ispin)%bands_from = b_from + wan_data(iwan,ispin)%bands_to = b_to + ! + do j=1, ning + call read_line( input_line, end_of_file = tend ) + ! + if ( tend ) & + call errore( 'read_cards', & + 'not enough wavefunctions', j ) + if (ning==1) then + read(input_line,*) lo,im + c = 1.d0 + else + read(input_line,*) lo,im,c + endif + + select case(trim(lo)) + case('s') + il = 0 + case('p') + il = 1 + case('d') + il = 2 + case('f') + il = 3 + case default + call errore( 'read_cards', & + 'wrong l-label', 1 ) + end select + + wan_data(iwan,ispin)%ing(j)%l = il + wan_data(iwan,ispin)%ing(j)%m = im + wan_data(iwan,ispin)%ing(j)%c = c + enddo + enddo + else + ! oups - that is not our data - lets's move one line up in input file + ! not sure that a direct access to the parce_unit is safe enougth + backspace(parse_unit) + endif + else + ! ok, that's the end of file. but i will move one line up + ! for a correct handling of eof in the parent read_cards subroutine + ! otherwise (at least with gfortran on mac) there will be the read error + backspace(parse_unit) + endif + ! + return + ! + end subroutine card_wannier_ac +end module read_cards_module diff --git a/src/read_input.f90 b/src/read_input.f90 new file mode 100644 index 0000000..0d8461d --- /dev/null +++ b/src/read_input.f90 @@ -0,0 +1,77 @@ +! +! copyright (c) 2011 quantum espresso group +! this file is distributed under the terms of the +! gnu general public license. see the file `license' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! read input data in pw and cp from text file (xml file to be implemented) +! note aug 2018 (pg): reading of old xml input file using iotk deleted +! +!---------------------------------------------------------------------------- +module read_input + !--------------------------------------------------------------------------- + ! + use kinds, only: dp + ! + implicit none + save + ! + private + public :: read_input_file, has_been_read + ! + logical :: has_been_read = .false. + ! + contains + ! + !------------------------------------------------------------------------- + subroutine read_input_file ( prog, input_file_ ) + !------------------------------------------------------------------------- + ! + use input_parameters, only : reset_input_checks + use read_namelists_module, only : read_namelists + use read_cards_module, only : read_cards + use io_global, only : ionode, ionode_id, qestdin + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + use open_close_input_file, only : open_input_file, close_input_file + ! + implicit none + ! + character(len=*), intent (in) :: prog + character(len=*), intent (in) :: input_file_ + ! + logical :: xmlinput + integer :: ierr + ! + if ( ionode ) ierr = open_input_file( input_file_, xmlinput ) + call mp_bcast( ierr, ionode_id, intra_image_comm ) + if ( ierr > 0 ) call errore('read_input', 'opening input file',ierr) + call mp_bcast( xmlinput, ionode_id, intra_image_comm ) + ! + call reset_input_checks () + ! + if ( xmlinput ) then + ! + call errore('read_input', 'xml input disabled',1) + ! + else + ! + ! ... read namelists + ! + call read_namelists( prog, qestdin ) + ! + ! ... read cards (requires in input only first two letters of prog) + ! + call read_cards ( prog(1:2), qestdin ) + ! + end if + if ( ionode) ierr = close_input_file( ) + ! + has_been_read = .true. + ! + return + ! + end subroutine read_input_file + ! +end module read_input diff --git a/src/read_namelists.f90 b/src/read_namelists.f90 new file mode 100644 index 0000000..6c63e1c --- /dev/null +++ b/src/read_namelists.f90 @@ -0,0 +1,2055 @@ +! +! copyright (c) 2002-2020 quantum espresso group +! this file is distributed under the terms of the +! gnu general public license. see the file `license' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +!--------------------------------------------- +! tb +! included gate related stuff, search 'tb' +!--------------------------------------------- +! +!---------------------------------------------------------------------------- +module read_namelists_module + !---------------------------------------------------------------------------- + ! + ! ... this module handles the reading of input namelists + ! ... written by carlo cavazzoni, with many additions + ! -------------------------------------------------- + ! + use kinds, only : dp + use input_parameters + ! + implicit none + ! + save + ! + private + ! + real(dp), parameter :: sm_not_set = -20.0_dp + ! + public :: read_namelists, sm_not_set + public :: check_namelist_read ! made public upon request of a.jay + ! fixme: should the following ones be public? + public :: control_defaults, system_defaults, & + electrons_defaults, wannier_ac_defaults, ions_defaults, & + cell_defaults, press_ai_defaults, wannier_defaults, control_bcast,& + system_bcast, electrons_bcast, ions_bcast, cell_bcast, & + press_ai_bcast, wannier_bcast, wannier_ac_bcast, control_checkin, & + system_checkin, electrons_checkin, ions_checkin, cell_checkin, & + wannier_checkin, wannier_ac_checkin, fixval + ! + ! ... end of module-scope declarations + ! + ! ---------------------------------------------- + ! + contains + ! + !=-----------------------------------------------------------------------=! + ! + ! variables initialization for namelist control + ! + !=-----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine control_defaults( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: temp_string + ! + ! + if ( prog == 'pw' ) then + title = ' ' + calculation = 'scf' + else + title = 'md simulation' + calculation = 'cp' + end if + + verbosity = 'default' + if( prog == 'pw' ) restart_mode = 'from_scratch' + if( prog == 'cp' ) restart_mode = 'restart' + nstep = 50 + if( prog == 'pw' ) iprint = 100000 + if( prog == 'cp' ) iprint = 10 + if( prog == 'pw' ) isave = 0 + if( prog == 'cp' ) isave = 100 + ! + tstress = .false. + tprnfor = .false. + tabps = .false. + ! + if( prog == 'pw' ) dt = 20.0_dp + if( prog == 'cp' ) dt = 1.0_dp + ! + ndr = 50 + ndw = 50 + ! + ! ... use the path specified as outdir and the filename prefix + ! ... to store output data + ! + call get_environment_variable( 'espresso_tmpdir', outdir ) + if ( trim( outdir ) == ' ' ) outdir = './' + if( prog == 'pw' ) prefix = 'pwscf' + if( prog == 'cp' ) prefix = 'cp' + ! + ! ... directory containing the pseudopotentials + ! + call get_environment_variable( 'espresso_pseudo', pseudo_dir ) + if ( trim( pseudo_dir ) == ' ') then + call get_environment_variable( 'home', pseudo_dir ) + pseudo_dir = trim( pseudo_dir ) // '/espresso/pseudo/' + end if + ! + ! ... max number of md steps added to the xml file. needs to be limited for very long + ! md simulations + call get_environment_variable('max_xml_steps', temp_string) + if ( trim(temp_string) .ne. ' ') read(temp_string, *) max_xml_steps + refg = 0.05_dp + max_seconds = 1.e+7_dp + ekin_conv_thr = 1.e-6_dp + etot_conv_thr = 1.e-4_dp + forc_conv_thr = 1.e-3_dp + disk_io = 'default' + dipfield = .false. + gate = .false. !tb + lberry = .false. + gdir = 0 + nppstr = 0 + wf_collect = .true. + lelfield = .false. + lorbm = .false. + nberrycyc = 1 + lecrpa = .false. + tqmmm = .false. + ! + saverho = .true. + memory = 'default' + ! + lfcpopt = .false. + lfcpdyn = .false. + ! + call get_environment_variable( 'qexml', input_xml_schema_file ) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist system + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine system_defaults( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + ! + ibrav = -1 + celldm = (/ 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp /) + a = 0.0_dp + b = 0.0_dp + c = 0.0_dp + cosab = 0.0_dp + cosac = 0.0_dp + cosbc = 0.0_dp + nat = 0 + ntyp = 0 + nbnd = 0 + tot_charge = 0.0_dp + tot_magnetization = -1 + ecutwfc = 0.0_dp + ecutrho = 0.0_dp + nr1 = 0 + nr2 = 0 + nr3 = 0 + nr1s = 0 + nr2s = 0 + nr3s = 0 + nr1b = 0 + nr2b = 0 + nr3b = 0 + occupations = 'fixed' + smearing = 'gaussian' + degauss = 0.0_dp + nspin = 1 + nosym = .false. + nosym_evc = .false. + force_symmorphic = .false. + use_all_frac = .false. + noinv = .false. + ecfixed = 0.0_dp + qcutz = 0.0_dp + q2sigma = 0.01_dp + input_dft = 'none' + ecutfock = -1.0_dp + starting_charge = 0.0_dp +! +! ... set starting_magnetization to an invalid value: +! ... in pw starting_magnetization must be set for at least one atomic type +! ... (unless the magnetization is set in other ways) +! ... in cp starting_magnetization must remain unset +! + starting_magnetization = sm_not_set + + if ( prog == 'pw' ) then + starting_ns_eigenvalue = -1.0_dp + u_projection_type = 'atomic' + end if + ! + ! .. dft + u and its extensions + ! + lda_plus_u = .false. + lda_plus_u_kind = 0 + hubbard_u = 0.0_dp + hubbard_u_back = 0.0_dp + hubbard_v = 0.0_dp + hubbard_j0 = 0.0_dp + hubbard_j = 0.0_dp + hubbard_alpha = 0.0_dp + hubbard_alpha_back = 0.0_dp + hubbard_beta = 0.0_dp + hubbard_parameters = 'input' + reserv = .false. + reserv_back = .false. + backall = .false. + lback = -1 + l1back = -1 + hub_pot_fix = .false. + step_pen=.false. + a_pen=0.0_dp + sigma_pen=0.01_dp + alpha_pen=0.0_dp + ! + ! ... exx + ! + ace=.true. + n_proj = 0 + localization_thr = 0.0_dp + scdm=.false. + scdmden=1.0d0 + scdmgrd=1.0d0 + nscdm=1 + ! + ! ... electric fields + ! + edir = 1 + emaxpos = 0.5_dp + eopreg = 0.1_dp + eamp = 0.0_dp + ! tb gate related variables + zgate = 0.5 + relaxz = .false. + block = .false. + block_1 = 0.45 + block_2 = 0.55 + block_height = 0.0 + ! + ! ... postprocessing of dos & phonons & el-ph + ! + la2f = .false. + ! + ! ... non collinear program variables + ! + lspinorb = .false. + lforcet = .false. + starting_spin_angle=.false. + noncolin = .false. + lambda = 1.0_dp + constrained_magnetization= 'none' + fixed_magnetization = 0.0_dp + b_field = 0.0_dp + angle1 = 0.0_dp + angle2 = 0.0_dp + report =-1 + ! + no_t_rev = .false. + ! + assume_isolated = 'none' + ! + one_atom_occupations=.false. + ! + spline_ps = .false. + ! + real_space = .false. + ! + ! ... dft-d, tkatchenko-scheffler, xdm + ! + vdw_corr = 'none' + london = .false. + london_s6 = 0.75_dp + london_rcut = 200.00_dp + london_c6 = -1.0_dp + london_rvdw = -1.0_dp + ts_vdw = .false. + ts_vdw_isolated = .false. + ts_vdw_econv_thr = 1.e-6_dp + xdm = .false. + xdm_a1 = 0.0_dp + xdm_a2 = 0.0_dp + dftd3_version = 3 + dftd3_threebody = .true. + ! + ! ... esm + ! + esm_bc='pbc' + esm_efield=0.0_dp + esm_w=0.0_dp + esm_a=0.0_dp + esm_zb=-2.0_dp + esm_nfit=4 + esm_debug=.false. + esm_debug_gpmax=0 + ! + ! ... fcp + ! + fcp_mu = 0.0_dp + fcp_mass = 10000.0_dp + fcp_tempw = 0.0_dp + fcp_relax = 'lm' + fcp_relax_step = 0.5_dp + fcp_relax_crit = 0.001_dp + fcp_mdiis_size = 4 + fcp_mdiis_step = 0.2_dp + ! + ! ... wyckoff + ! + space_group=0 + uniqueb = .false. + origin_choice = 1 + rhombohedral = .true. + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist electrons + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine electrons_defaults( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + ! + emass = 400.0_dp + emass_cutoff = 2.5_dp + orthogonalization = 'ortho' + ortho_eps = 1.e-9_dp + ortho_max = 300 + electron_maxstep = 100 + scf_must_converge = .true. + ! + ! ... ( 'sd' | 'cg' | 'damp' | 'verlet' | 'none' | 'diis' | 'cp-bo' ) + ! + electron_dynamics = 'none' + electron_damping = 0.1_dp + ! + ! ... ( 'zero' | 'default' ) + ! + electron_velocities = 'default' + ! + ! ... ( 'nose' | 'not_controlled' | 'rescaling') + ! + electron_temperature = 'not_controlled' + ekincw = 0.001_dp + fnosee = 1.0_dp + ampre = 0.0_dp + grease = 1.0_dp + conv_thr = 1.e-6_dp + diis_size = 4 + diis_nreset = 3 + diis_hcut = 1.0_dp + diis_wthr = 0.0_dp + diis_delt = 0.0_dp + diis_maxstep = 100 + diis_rot = .false. + diis_fthr = 0.0_dp + diis_temp = 0.0_dp + diis_achmix = 0.0_dp + diis_g0chmix = 0.0_dp + diis_g1chmix = 0.0_dp + diis_nchmix = 3 + diis_nrot = 3 + diis_rothr = 0.0_dp + diis_ethr = 0.0_dp + diis_chguess = .false. + mixing_mode = 'plain' + mixing_fixed_ns = 0 + mixing_beta = 0.7_dp + mixing_ndim = 8 + diagonalization = 'david' + diago_thr_init = 0.0_dp + diago_cg_maxiter = 20 + diago_ppcg_maxiter = 20 + diago_david_ndim = 2 + diago_full_acc = .false. + ! + sic = 'none' + sic_epsilon = 0.0_dp + sic_alpha = 0.0_dp + force_pairing = .false. + ! + fermi_energy = 0.0_dp + n_inner = 2 + niter_cold_restart=1 + lambda_cold=0.03_dp + rotation_dynamics = "line-minimization" + occupation_dynamics = "line-minimization" + rotmass = 0.0_dp + occmass = 0.0_dp + rotation_damping = 0.0_dp + occupation_damping = 0.0_dp + ! + tcg = .false. + maxiter = 100 + passop = 0.3_dp + niter_cg_restart = 20 + etresh = 1.e-6_dp + ! + epol = 3 + efield = 0.0_dp + epol2 = 3 + efield2 = 0.0_dp + efield_cart(1)=0.d0 + efield_cart(2)=0.d0 + efield_cart(3)=0.d0 + efield_phase='none' + ! + occupation_constraints = .false. + ! + adaptive_thr = .false. + conv_thr_init = 0.1e-2_dp + conv_thr_multi = 0.1_dp + ! + ! ... cp-bo ... + tcpbo = .false. + emass_emin = 200.0_dp + emass_cutoff_emin = 6.0_dp + electron_damping_emin = 0.35_dp + dt_emin = 4.0_dp + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist wannier_ac + ! + !---------------------------------------------------------------------- + subroutine wannier_ac_defaults( prog ) + !---------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + ! + plot_wannier = .false. + use_energy_int = .false. + print_wannier_coeff = .false. + nwan = 0 + constrain_pot = 0.d0 + plot_wan_num = 0 + plot_wan_spin = 1 + ! + return + ! + end subroutine + + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist ions + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine ions_defaults( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + ! ... ( 'sd' | 'cg' | 'damp' | 'verlet' | 'none' | 'bfgs' | 'beeman' ) + ! + ion_dynamics = 'none' + ion_radius = 0.5_dp + ion_damping = 0.1_dp + ! + ! ... ( 'default' | 'from_input' ) + ! + ion_positions = 'default' + ! + ! ... ( 'zero' | 'default' | 'from_input' ) + ! + ion_velocities = 'default' + ! + ! ... ( 'nose' | 'not_controlled' | 'rescaling' | 'berendsen' | + ! 'andersen' | 'initial' ) + ! + ion_temperature = 'not_controlled' + ! + tempw = 300.0_dp + fnosep = -1.0_dp + fnosep(1) = 1.0_dp + nhpcl = 0 + nhptyp = 0 + ndega = 0 + tranp = .false. + amprp = 0.0_dp + greasp = 1.0_dp + tolp = 100.0_dp + ion_nstepe = 1 + ion_maxstep = 100 + delta_t = 1.0_dp + nraise = 1 + ! + refold_pos = .false. + remove_rigid_rot = .false. + ! + upscale = 100.0_dp + pot_extrapolation = 'atomic' + wfc_extrapolation = 'none' + ! + ! ... bfgs defaults + ! + bfgs_ndim = 1 + trust_radius_max = 0.8_dp ! bohr + trust_radius_min = 1.e-4_dp ! bohr + trust_radius_ini = 0.5_dp ! bohr + w_1 = 0.01_dp + w_2 = 0.50_dp + ! + l_mplathe=.false. + n_muller=0 + np_muller=1 + l_exit_muller=.false. + + + return + ! + end subroutine + ! + ! + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist cell + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine cell_defaults( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + ! + cell_parameters = 'default' + ! + ! ... ( 'sd' | 'pr' | 'none' | 'w' | 'damp-pr' | 'damp-w' | 'bfgs' ) + ! + cell_dynamics = 'none' + ! + ! ... ( 'zero' | 'default' ) + ! + cell_velocities = 'default' + press = 0.0_dp + wmass = 0.0_dp + ! + ! ... ( 'nose' | 'not_controlled' | 'rescaling' ) + ! + cell_temperature = 'not_controlled' + temph = 0.0_dp + fnoseh = 1.0_dp + greash = 1.0_dp + ! + ! ... ('all'* | 'volume' | 'x' | 'y' | 'z' | 'xy' | 'xz' | 'yz' | 'xyz' ) + ! + cell_dofree = 'all' + cell_factor = 0.0_dp + cell_nstepe = 1 + cell_damping = 0.1_dp + press_conv_thr = 0.5_dp + treinit_gvecs = .false. + ! + return + ! + end subroutine + ! + ! + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist press_ai + ! + !=----------------------------------------------------------------------=! + ! + !---------------------------------------------------------------------- + subroutine press_ai_defaults( prog ) + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + abivol = .false. + abisur = .false. + pvar = .false. + fill_vac = .false. + cntr = .false. + scale_at = .false. + t_gauss = .false. + jellium = .false. + + p_ext = 0.0_dp + p_in = 0.0_dp + p_fin = 0.0_dp + surf_t = 0.0_dp + rho_thr = 0.0_dp + dthr = 0.0_dp + step_rad = 0.0_dp + delta_eps = 0.0_dp + delta_sigma = 0.0_dp + r_j = 0.0_dp + h_j = 0.0_dp + + n_cntr = 0 + axis = 3 + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! variables initialization for namelist wannier + ! + !----------------------------------------------------------------------- + subroutine wannier_defaults( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + ! + ! + wf_efield = .false. + wf_switch = .false. + ! + sw_len = 1 + ! + efx0 = 0.0_dp + efy0 = 0.0_dp + efz0 = 0.0_dp + efx1 = 0.0_dp + efy1 = 0.0_dp + efz1 = 0.0_dp + ! + wfsd = 1 + ! + wfdt = 5.0_dp + maxwfdt = 0.30_dp + wf_q = 1500.0_dp + wf_friction = 0.3_dp +!======================================================================= +!exx_wf related + exx_neigh = 60 + vnbsp = 0 + exx_poisson_eps = 1.e-6_dp + exx_dis_cutoff = 8.0_dp + exx_ps_rcut_self = 6.0_dp + exx_ps_rcut_pair = 5.0_dp + exx_me_rcut_self = 10.0_dp + exx_me_rcut_pair = 7.0_dp + exx_use_cube_domain = .false. +!======================================================================= + ! + nit = 10 + nsd = 10 + nsteps = 20 + ! + tolw = 1.e-8_dp + ! + adapt = .true. + ! + calwf = 3 + nwf = 0 + wffort = 40 + ! + writev = .false. + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist control + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine control_bcast() + !----------------------------------------------------------------------- + ! + use io_global, only : ionode_id + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + call mp_bcast( title, ionode_id, intra_image_comm ) + call mp_bcast( calculation, ionode_id, intra_image_comm ) + call mp_bcast( verbosity, ionode_id, intra_image_comm ) + call mp_bcast( restart_mode, ionode_id, intra_image_comm ) + call mp_bcast( nstep, ionode_id, intra_image_comm ) + call mp_bcast( iprint, ionode_id, intra_image_comm ) + call mp_bcast( isave, ionode_id, intra_image_comm ) + call mp_bcast( tstress, ionode_id, intra_image_comm ) + call mp_bcast( tprnfor, ionode_id, intra_image_comm ) + call mp_bcast( tabps, ionode_id, intra_image_comm ) + call mp_bcast( dt, ionode_id, intra_image_comm ) + call mp_bcast( ndr, ionode_id, intra_image_comm ) + call mp_bcast( ndw, ionode_id, intra_image_comm ) + call mp_bcast( outdir, ionode_id, intra_image_comm ) + call mp_bcast( wfcdir, ionode_id, intra_image_comm ) + call mp_bcast( prefix, ionode_id, intra_image_comm ) + call mp_bcast( max_seconds, ionode_id, intra_image_comm ) + call mp_bcast( ekin_conv_thr, ionode_id, intra_image_comm ) + call mp_bcast( etot_conv_thr, ionode_id, intra_image_comm ) + call mp_bcast( forc_conv_thr, ionode_id, intra_image_comm ) + call mp_bcast( pseudo_dir, ionode_id, intra_image_comm ) + call mp_bcast( refg, ionode_id, intra_image_comm ) + call mp_bcast( disk_io, ionode_id, intra_image_comm ) + call mp_bcast( tefield, ionode_id, intra_image_comm ) + call mp_bcast( tefield2, ionode_id, intra_image_comm ) + call mp_bcast( dipfield, ionode_id, intra_image_comm ) + call mp_bcast( lberry, ionode_id, intra_image_comm ) + call mp_bcast( gdir, ionode_id, intra_image_comm ) + call mp_bcast( nppstr, ionode_id, intra_image_comm ) + call mp_bcast( point_label_type, ionode_id, intra_image_comm ) + call mp_bcast( wf_collect, ionode_id, intra_image_comm ) + call mp_bcast( lelfield, ionode_id, intra_image_comm ) + call mp_bcast( lorbm, ionode_id, intra_image_comm ) + call mp_bcast( nberrycyc, ionode_id, intra_image_comm ) + call mp_bcast( saverho, ionode_id, intra_image_comm ) + call mp_bcast( lecrpa, ionode_id, intra_image_comm ) + call mp_bcast( tqmmm, ionode_id, intra_image_comm ) + call mp_bcast( vdw_table_name,ionode_id, intra_image_comm ) + call mp_bcast( memory, ionode_id, intra_image_comm ) + call mp_bcast( lfcpopt, ionode_id, intra_image_comm ) + call mp_bcast( lfcpdyn, ionode_id, intra_image_comm ) + call mp_bcast( input_xml_schema_file, ionode_id, intra_image_comm ) + call mp_bcast( gate, ionode_id, intra_image_comm ) !tb + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist system + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine system_bcast() + !----------------------------------------------------------------------- + ! + use io_global, only : ionode_id + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + call mp_bcast( ibrav, ionode_id, intra_image_comm ) + call mp_bcast( celldm, ionode_id, intra_image_comm ) + call mp_bcast( a, ionode_id, intra_image_comm ) + call mp_bcast( b, ionode_id, intra_image_comm ) + call mp_bcast( c, ionode_id, intra_image_comm ) + call mp_bcast( cosab, ionode_id, intra_image_comm ) + call mp_bcast( cosac, ionode_id, intra_image_comm ) + call mp_bcast( cosbc, ionode_id, intra_image_comm ) + call mp_bcast( nat, ionode_id, intra_image_comm ) + call mp_bcast( ntyp, ionode_id, intra_image_comm ) + call mp_bcast( nbnd, ionode_id, intra_image_comm ) + call mp_bcast( tot_charge, ionode_id, intra_image_comm ) + call mp_bcast( tot_magnetization, ionode_id, intra_image_comm ) + call mp_bcast( ecutwfc, ionode_id, intra_image_comm ) + call mp_bcast( ecutrho, ionode_id, intra_image_comm ) + call mp_bcast( nr1, ionode_id, intra_image_comm ) + call mp_bcast( nr2, ionode_id, intra_image_comm ) + call mp_bcast( nr3, ionode_id, intra_image_comm ) + call mp_bcast( nr1s, ionode_id, intra_image_comm ) + call mp_bcast( nr2s, ionode_id, intra_image_comm ) + call mp_bcast( nr3s, ionode_id, intra_image_comm ) + call mp_bcast( nr1b, ionode_id, intra_image_comm ) + call mp_bcast( nr2b, ionode_id, intra_image_comm ) + call mp_bcast( nr3b, ionode_id, intra_image_comm ) + call mp_bcast( occupations, ionode_id, intra_image_comm ) + call mp_bcast( smearing, ionode_id, intra_image_comm ) + call mp_bcast( degauss, ionode_id, intra_image_comm ) + call mp_bcast( nspin, ionode_id, intra_image_comm ) + call mp_bcast( nosym, ionode_id, intra_image_comm ) + call mp_bcast( nosym_evc, ionode_id, intra_image_comm ) + call mp_bcast( noinv, ionode_id, intra_image_comm ) + call mp_bcast( force_symmorphic, ionode_id, intra_image_comm ) + call mp_bcast( use_all_frac, ionode_id, intra_image_comm ) + call mp_bcast( ecfixed, ionode_id, intra_image_comm ) + call mp_bcast( qcutz, ionode_id, intra_image_comm ) + call mp_bcast( q2sigma, ionode_id, intra_image_comm ) + call mp_bcast( input_dft, ionode_id, intra_image_comm ) + + ! ... exx + + call mp_bcast( ace, ionode_id, intra_image_comm ) + call mp_bcast( localization_thr, ionode_id, intra_image_comm ) + call mp_bcast( scdm, ionode_id, intra_image_comm ) + call mp_bcast( scdmden, ionode_id, intra_image_comm ) + call mp_bcast( scdmgrd, ionode_id, intra_image_comm ) + call mp_bcast( nscdm, ionode_id, intra_image_comm ) + call mp_bcast( n_proj, ionode_id, intra_image_comm ) + call mp_bcast( nqx1, ionode_id, intra_image_comm ) + call mp_bcast( nqx2, ionode_id, intra_image_comm ) + call mp_bcast( nqx3, ionode_id, intra_image_comm ) + call mp_bcast( exx_fraction, ionode_id, intra_image_comm ) + call mp_bcast( screening_parameter, ionode_id, intra_image_comm ) + call mp_bcast( gau_parameter, ionode_id, intra_image_comm ) + call mp_bcast( exxdiv_treatment, ionode_id, intra_image_comm ) + call mp_bcast( x_gamma_extrapolation, ionode_id, intra_image_comm ) + call mp_bcast( yukawa, ionode_id, intra_image_comm ) + call mp_bcast( ecutvcut, ionode_id, intra_image_comm ) + call mp_bcast( ecutfock, ionode_id, intra_image_comm ) + ! + call mp_bcast( starting_charge, ionode_id, intra_image_comm ) + call mp_bcast( starting_magnetization, ionode_id, intra_image_comm ) + call mp_bcast( starting_ns_eigenvalue, ionode_id, intra_image_comm ) + call mp_bcast( u_projection_type, ionode_id, intra_image_comm ) + call mp_bcast( lda_plus_u, ionode_id, intra_image_comm ) + call mp_bcast( lda_plus_u_kind, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_u, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_u_back, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_j0, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_j, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_v, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_alpha, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_alpha_back, ionode_id, intra_image_comm ) + call mp_bcast( hubbard_beta, ionode_id, intra_image_comm ) + call mp_bcast( hub_pot_fix, ionode_id,intra_image_comm ) + call mp_bcast( hubbard_parameters, ionode_id,intra_image_comm ) + call mp_bcast( reserv, ionode_id,intra_image_comm ) + call mp_bcast( reserv_back, ionode_id,intra_image_comm ) + call mp_bcast( backall, ionode_id,intra_image_comm ) + call mp_bcast( lback, ionode_id,intra_image_comm ) + call mp_bcast( l1back, ionode_id,intra_image_comm ) + call mp_bcast( step_pen, ionode_id, intra_image_comm ) + call mp_bcast( a_pen, ionode_id, intra_image_comm ) + call mp_bcast( sigma_pen, ionode_id, intra_image_comm ) + call mp_bcast( alpha_pen, ionode_id, intra_image_comm ) + call mp_bcast( edir, ionode_id, intra_image_comm ) + call mp_bcast( emaxpos, ionode_id, intra_image_comm ) + call mp_bcast( eopreg, ionode_id, intra_image_comm ) + call mp_bcast( eamp, ionode_id, intra_image_comm ) + call mp_bcast( la2f, ionode_id, intra_image_comm ) + ! + ! ... non collinear broadcast + ! + call mp_bcast( lspinorb, ionode_id, intra_image_comm ) + call mp_bcast( lforcet, ionode_id, intra_image_comm ) + call mp_bcast( starting_spin_angle, ionode_id, intra_image_comm ) + call mp_bcast( noncolin, ionode_id, intra_image_comm ) + call mp_bcast( angle1, ionode_id, intra_image_comm ) + call mp_bcast( angle2, ionode_id, intra_image_comm ) + call mp_bcast( report, ionode_id, intra_image_comm ) + call mp_bcast( constrained_magnetization, ionode_id, intra_image_comm ) + call mp_bcast( b_field, ionode_id, intra_image_comm ) + call mp_bcast( fixed_magnetization, ionode_id, intra_image_comm ) + call mp_bcast( lambda, ionode_id, intra_image_comm ) + ! + call mp_bcast( assume_isolated, ionode_id, intra_image_comm ) + call mp_bcast( one_atom_occupations, ionode_id, intra_image_comm ) + call mp_bcast( spline_ps, ionode_id, intra_image_comm ) + ! + call mp_bcast( vdw_corr, ionode_id, intra_image_comm ) + call mp_bcast( ts_vdw, ionode_id, intra_image_comm ) + call mp_bcast( ts_vdw_isolated, ionode_id, intra_image_comm ) + call mp_bcast( ts_vdw_econv_thr, ionode_id, intra_image_comm ) + call mp_bcast( london, ionode_id, intra_image_comm ) + call mp_bcast( london_s6, ionode_id, intra_image_comm ) + call mp_bcast( london_rcut, ionode_id, intra_image_comm ) + call mp_bcast( london_c6, ionode_id, intra_image_comm ) + call mp_bcast( london_rvdw, ionode_id, intra_image_comm ) + call mp_bcast( xdm, ionode_id, intra_image_comm ) + call mp_bcast( xdm_a1, ionode_id, intra_image_comm ) + call mp_bcast( xdm_a2, ionode_id, intra_image_comm ) + ! + call mp_bcast( no_t_rev, ionode_id, intra_image_comm ) + ! + ! ... esm method broadcast + ! + call mp_bcast( esm_bc, ionode_id, intra_image_comm ) + call mp_bcast( esm_efield, ionode_id, intra_image_comm ) + call mp_bcast( esm_w, ionode_id, intra_image_comm ) + call mp_bcast( esm_a, ionode_id, intra_image_comm ) + call mp_bcast( esm_zb, ionode_id, intra_image_comm ) + call mp_bcast( esm_nfit, ionode_id, intra_image_comm ) + call mp_bcast( esm_debug, ionode_id, intra_image_comm ) + call mp_bcast( esm_debug_gpmax, ionode_id, intra_image_comm ) + ! + ! ... fcp + ! + call mp_bcast( fcp_mu, ionode_id, intra_image_comm ) + call mp_bcast( fcp_mass, ionode_id, intra_image_comm ) + call mp_bcast( fcp_tempw, ionode_id, intra_image_comm ) + call mp_bcast( fcp_relax, ionode_id, intra_image_comm ) + call mp_bcast( fcp_relax_step, ionode_id, intra_image_comm ) + call mp_bcast( fcp_relax_crit, ionode_id, intra_image_comm ) + call mp_bcast( fcp_mdiis_size, ionode_id, intra_image_comm ) + call mp_bcast( fcp_mdiis_step, ionode_id, intra_image_comm ) + ! + ! + ! ... space group information + ! + call mp_bcast( space_group, ionode_id, intra_image_comm ) + call mp_bcast( uniqueb, ionode_id, intra_image_comm ) + call mp_bcast( origin_choice, ionode_id, intra_image_comm ) + call mp_bcast( rhombohedral, ionode_id, intra_image_comm ) + ! + ! tb - gate broadcast + ! + call mp_bcast( zgate, ionode_id, intra_image_comm ) + call mp_bcast( relaxz, ionode_id, intra_image_comm ) + call mp_bcast( block, ionode_id, intra_image_comm ) + call mp_bcast( block_1, ionode_id, intra_image_comm ) + call mp_bcast( block_2, ionode_id, intra_image_comm ) + call mp_bcast( block_height, ionode_id, intra_image_comm ) + + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist electrons + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine electrons_bcast() + !----------------------------------------------------------------------- + ! + use io_global, only : ionode_id + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + call mp_bcast( emass, ionode_id, intra_image_comm ) + call mp_bcast( emass_cutoff, ionode_id, intra_image_comm ) + call mp_bcast( orthogonalization, ionode_id, intra_image_comm ) + call mp_bcast( electron_maxstep, ionode_id, intra_image_comm ) + call mp_bcast( scf_must_converge, ionode_id, intra_image_comm ) + call mp_bcast( ortho_eps, ionode_id, intra_image_comm ) + call mp_bcast( ortho_max, ionode_id, intra_image_comm ) + call mp_bcast( electron_dynamics, ionode_id, intra_image_comm ) + call mp_bcast( electron_damping, ionode_id, intra_image_comm ) + call mp_bcast( electron_velocities, ionode_id, intra_image_comm ) + call mp_bcast( electron_temperature, ionode_id, intra_image_comm ) + call mp_bcast( conv_thr, ionode_id, intra_image_comm ) + call mp_bcast( ekincw, ionode_id, intra_image_comm ) + call mp_bcast( fnosee, ionode_id, intra_image_comm ) + call mp_bcast( startingwfc, ionode_id, intra_image_comm ) + call mp_bcast( ampre, ionode_id, intra_image_comm ) + call mp_bcast( grease, ionode_id, intra_image_comm ) + call mp_bcast( startingpot, ionode_id, intra_image_comm ) + call mp_bcast( diis_size, ionode_id, intra_image_comm ) + call mp_bcast( diis_nreset, ionode_id, intra_image_comm ) + call mp_bcast( diis_hcut, ionode_id, intra_image_comm ) + call mp_bcast( diis_wthr, ionode_id, intra_image_comm ) + call mp_bcast( diis_delt, ionode_id, intra_image_comm ) + call mp_bcast( diis_maxstep, ionode_id, intra_image_comm ) + call mp_bcast( diis_rot, ionode_id, intra_image_comm ) + call mp_bcast( diis_fthr, ionode_id, intra_image_comm ) + call mp_bcast( diis_temp, ionode_id, intra_image_comm ) + call mp_bcast( diis_achmix, ionode_id, intra_image_comm ) + call mp_bcast( diis_g0chmix, ionode_id, intra_image_comm ) + call mp_bcast( diis_g1chmix, ionode_id, intra_image_comm ) + call mp_bcast( diis_nchmix, ionode_id, intra_image_comm ) + call mp_bcast( diis_nrot, ionode_id, intra_image_comm ) + call mp_bcast( diis_rothr, ionode_id, intra_image_comm ) + call mp_bcast( diis_ethr, ionode_id, intra_image_comm ) + call mp_bcast( diis_chguess, ionode_id, intra_image_comm ) + call mp_bcast( mixing_fixed_ns, ionode_id, intra_image_comm ) + call mp_bcast( mixing_mode, ionode_id, intra_image_comm ) + call mp_bcast( mixing_beta, ionode_id, intra_image_comm ) + call mp_bcast( mixing_ndim, ionode_id, intra_image_comm ) + call mp_bcast( tqr, ionode_id, intra_image_comm ) + call mp_bcast( tq_smoothing, ionode_id, intra_image_comm ) + call mp_bcast( tbeta_smoothing, ionode_id, intra_image_comm ) + call mp_bcast( diagonalization, ionode_id, intra_image_comm ) + call mp_bcast( diago_thr_init, ionode_id, intra_image_comm ) + call mp_bcast( diago_cg_maxiter, ionode_id, intra_image_comm ) + call mp_bcast( diago_ppcg_maxiter, ionode_id, intra_image_comm ) + call mp_bcast( diago_david_ndim, ionode_id, intra_image_comm ) + call mp_bcast( diago_full_acc, ionode_id, intra_image_comm ) + call mp_bcast( sic, ionode_id, intra_image_comm ) + call mp_bcast( sic_epsilon , ionode_id, intra_image_comm ) + call mp_bcast( sic_alpha , ionode_id, intra_image_comm ) + call mp_bcast( force_pairing , ionode_id, intra_image_comm ) + ! + ! ... ensemble-dft + ! + call mp_bcast( fermi_energy, ionode_id, intra_image_comm ) + call mp_bcast( n_inner, ionode_id, intra_image_comm ) + call mp_bcast( niter_cold_restart, ionode_id, intra_image_comm ) + call mp_bcast( lambda_cold, ionode_id, intra_image_comm ) + call mp_bcast( rotation_dynamics, ionode_id, intra_image_comm ) + call mp_bcast( occupation_dynamics,ionode_id, intra_image_comm ) + call mp_bcast( rotmass, ionode_id, intra_image_comm ) + call mp_bcast( occmass, ionode_id, intra_image_comm ) + call mp_bcast( rotation_damping, ionode_id, intra_image_comm ) + call mp_bcast( occupation_damping, ionode_id, intra_image_comm ) + ! + ! ... conjugate gradient + ! + call mp_bcast( tcg, ionode_id, intra_image_comm ) + call mp_bcast( maxiter, ionode_id, intra_image_comm ) + call mp_bcast( etresh, ionode_id, intra_image_comm ) + call mp_bcast( passop, ionode_id, intra_image_comm ) + call mp_bcast( niter_cg_restart, ionode_id, intra_image_comm ) + ! + ! ... electric field + ! + call mp_bcast( epol, ionode_id, intra_image_comm ) + call mp_bcast( efield, ionode_id, intra_image_comm ) + ! + call mp_bcast( epol2, ionode_id, intra_image_comm ) + call mp_bcast( efield2, ionode_id, intra_image_comm ) + call mp_bcast( efield_cart, ionode_id, intra_image_comm ) + call mp_bcast( efield_phase, ionode_id, intra_image_comm ) + ! + ! ... occupation constraints ... + ! + call mp_bcast( occupation_constraints, ionode_id, intra_image_comm ) + ! + ! ... real space ... + call mp_bcast( real_space, ionode_id, intra_image_comm ) + call mp_bcast( adaptive_thr, ionode_id, intra_image_comm ) + call mp_bcast( conv_thr_init, ionode_id, intra_image_comm ) + call mp_bcast( conv_thr_multi, ionode_id, intra_image_comm ) + ! + ! ... cp-bo ... + call mp_bcast( tcpbo, ionode_id, intra_image_comm ) + call mp_bcast( emass_emin, ionode_id, intra_image_comm ) + call mp_bcast( emass_cutoff_emin, ionode_id, intra_image_comm ) + call mp_bcast( electron_damping_emin, ionode_id, intra_image_comm ) + call mp_bcast( dt_emin, ionode_id, intra_image_comm ) + ! + return + ! + end subroutine + ! + ! + !=----------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist ions + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine ions_bcast() + !----------------------------------------------------------------------- + ! + use io_global, only: ionode_id + use mp, only: mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + call mp_bcast( ion_dynamics, ionode_id, intra_image_comm ) + call mp_bcast( ion_radius, ionode_id, intra_image_comm ) + call mp_bcast( ion_damping, ionode_id, intra_image_comm ) + call mp_bcast( ion_positions, ionode_id, intra_image_comm ) + call mp_bcast( ion_velocities, ionode_id, intra_image_comm ) + call mp_bcast( ion_temperature, ionode_id, intra_image_comm ) + call mp_bcast( tempw, ionode_id, intra_image_comm ) + call mp_bcast( fnosep, ionode_id, intra_image_comm ) + call mp_bcast( nhgrp, ionode_id, intra_image_comm ) + call mp_bcast( fnhscl, ionode_id, intra_image_comm ) + call mp_bcast( nhpcl, ionode_id, intra_image_comm ) + call mp_bcast( nhptyp, ionode_id, intra_image_comm ) + call mp_bcast( ndega, ionode_id, intra_image_comm ) + call mp_bcast( tranp, ionode_id, intra_image_comm ) + call mp_bcast( amprp, ionode_id, intra_image_comm ) + call mp_bcast( greasp, ionode_id, intra_image_comm ) + call mp_bcast( tolp, ionode_id, intra_image_comm ) + call mp_bcast( ion_nstepe, ionode_id, intra_image_comm ) + call mp_bcast( ion_maxstep, ionode_id, intra_image_comm ) + call mp_bcast( delta_t, ionode_id, intra_image_comm ) + call mp_bcast( nraise, ionode_id, intra_image_comm ) + call mp_bcast( refold_pos, ionode_id, intra_image_comm ) + call mp_bcast( remove_rigid_rot, ionode_id, intra_image_comm ) + call mp_bcast( upscale, ionode_id, intra_image_comm ) + call mp_bcast( pot_extrapolation, ionode_id, intra_image_comm ) + call mp_bcast( wfc_extrapolation, ionode_id, intra_image_comm ) + ! + ! ... bfgs + ! + call mp_bcast( bfgs_ndim, ionode_id, intra_image_comm ) + call mp_bcast( trust_radius_max, ionode_id, intra_image_comm ) + call mp_bcast( trust_radius_min, ionode_id, intra_image_comm ) + call mp_bcast( trust_radius_ini, ionode_id, intra_image_comm ) + call mp_bcast( w_1, ionode_id, intra_image_comm ) + call mp_bcast( w_2, ionode_id, intra_image_comm ) + ! + call mp_bcast(l_mplathe, ionode_id, intra_image_comm ) + call mp_bcast(n_muller, ionode_id, intra_image_comm ) + call mp_bcast(np_muller, ionode_id, intra_image_comm ) + call mp_bcast(l_exit_muller, ionode_id, intra_image_comm ) + + + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist cell + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine cell_bcast() + !----------------------------------------------------------------------- + ! + use io_global, only: ionode_id + use mp, only: mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + call mp_bcast( cell_parameters, ionode_id, intra_image_comm ) + call mp_bcast( cell_dynamics, ionode_id, intra_image_comm ) + call mp_bcast( cell_velocities, ionode_id, intra_image_comm ) + call mp_bcast( cell_dofree, ionode_id, intra_image_comm ) + call mp_bcast( press, ionode_id, intra_image_comm ) + call mp_bcast( wmass, ionode_id, intra_image_comm ) + call mp_bcast( cell_temperature, ionode_id, intra_image_comm ) + call mp_bcast( temph, ionode_id, intra_image_comm ) + call mp_bcast( fnoseh, ionode_id, intra_image_comm ) + call mp_bcast( greash, ionode_id, intra_image_comm ) + call mp_bcast( cell_factor, ionode_id, intra_image_comm ) + call mp_bcast( cell_nstepe, ionode_id, intra_image_comm ) + call mp_bcast( cell_damping, ionode_id, intra_image_comm ) + call mp_bcast( press_conv_thr, ionode_id, intra_image_comm ) + call mp_bcast( treinit_gvecs, ionode_id, intra_image_comm ) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist press_ai + ! + !=----------------------------------------------------------------------=! + ! + !---------------------------------------------------------------------- + subroutine press_ai_bcast() + !---------------------------------------------------------------------- + ! + use io_global, only: ionode_id + use mp, only: mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + ! + call mp_bcast( abivol, ionode_id, intra_image_comm ) + call mp_bcast( abisur, ionode_id, intra_image_comm ) + call mp_bcast( t_gauss, ionode_id, intra_image_comm ) + call mp_bcast( cntr, ionode_id, intra_image_comm ) + call mp_bcast( p_ext, ionode_id, intra_image_comm ) + call mp_bcast( surf_t, ionode_id, intra_image_comm ) + call mp_bcast( pvar, ionode_id, intra_image_comm ) + call mp_bcast( p_in, ionode_id, intra_image_comm ) + call mp_bcast( p_fin, ionode_id, intra_image_comm ) + call mp_bcast( delta_eps, ionode_id, intra_image_comm ) + call mp_bcast( delta_sigma, ionode_id, intra_image_comm ) + call mp_bcast( fill_vac, ionode_id, intra_image_comm ) + call mp_bcast( scale_at, ionode_id, intra_image_comm ) + call mp_bcast( n_cntr, ionode_id, intra_image_comm ) + call mp_bcast( axis, ionode_id, intra_image_comm ) + call mp_bcast( rho_thr, ionode_id, intra_image_comm ) + call mp_bcast( dthr, ionode_id, intra_image_comm ) + call mp_bcast( step_rad, ionode_id, intra_image_comm ) + call mp_bcast( jellium, ionode_id, intra_image_comm ) + call mp_bcast( r_j, ionode_id, intra_image_comm ) + call mp_bcast( h_j, ionode_id, intra_image_comm ) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist wannier + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine wannier_bcast() + !----------------------------------------------------------------------- + ! + use io_global, only: ionode_id + use mp, only: mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + call mp_bcast( wf_efield, ionode_id, intra_image_comm ) + call mp_bcast( wf_switch, ionode_id, intra_image_comm ) + call mp_bcast( sw_len, ionode_id, intra_image_comm ) + call mp_bcast( efx0, ionode_id, intra_image_comm ) + call mp_bcast( efy0, ionode_id, intra_image_comm ) + call mp_bcast( efz0, ionode_id, intra_image_comm ) + call mp_bcast( efx1, ionode_id, intra_image_comm ) + call mp_bcast( efy1, ionode_id, intra_image_comm ) + call mp_bcast( efz1, ionode_id, intra_image_comm ) + call mp_bcast( wfsd, ionode_id, intra_image_comm ) + call mp_bcast( wfdt, ionode_id, intra_image_comm ) + call mp_bcast( maxwfdt, ionode_id, intra_image_comm ) + call mp_bcast( wf_q, ionode_id, intra_image_comm ) + call mp_bcast( wf_friction, ionode_id, intra_image_comm ) + call mp_bcast( nit, ionode_id, intra_image_comm ) + call mp_bcast( nsd, ionode_id, intra_image_comm ) + call mp_bcast( nsteps, ionode_id, intra_image_comm ) + call mp_bcast( tolw, ionode_id, intra_image_comm ) + call mp_bcast( adapt, ionode_id, intra_image_comm ) + call mp_bcast( calwf, ionode_id, intra_image_comm ) + call mp_bcast( nwf, ionode_id, intra_image_comm ) + call mp_bcast( wffort, ionode_id, intra_image_comm ) + call mp_bcast( writev, ionode_id, intra_image_comm ) +!================================================================= +!exx_wf related + call mp_bcast( exx_neigh, ionode_id, intra_image_comm ) + call mp_bcast( exx_poisson_eps, ionode_id, intra_image_comm ) + call mp_bcast( exx_dis_cutoff, ionode_id, intra_image_comm ) + call mp_bcast( exx_ps_rcut_self, ionode_id, intra_image_comm ) + call mp_bcast( exx_ps_rcut_pair, ionode_id, intra_image_comm ) + call mp_bcast( exx_me_rcut_self, ionode_id, intra_image_comm ) + call mp_bcast( exx_me_rcut_pair, ionode_id, intra_image_comm ) + call mp_bcast( exx_use_cube_domain, ionode_id, intra_image_comm ) + call mp_bcast( vnbsp, ionode_id, intra_image_comm ) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------------=! + ! + ! broadcast variables values for namelist wannier_new + ! + !=----------------------------------------------------------------------------=! + ! + !---------------------------------------------------------------------- + subroutine wannier_ac_bcast() + !---------------------------------------------------------------------- + ! + use io_global, only: ionode_id + use mp, only: mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + ! + call mp_bcast( plot_wannier,ionode_id, intra_image_comm ) + call mp_bcast( use_energy_int,ionode_id, intra_image_comm ) + call mp_bcast( print_wannier_coeff,ionode_id, intra_image_comm ) + call mp_bcast( nwan, ionode_id, intra_image_comm ) + call mp_bcast( plot_wan_num,ionode_id, intra_image_comm ) + call mp_bcast( plot_wan_spin,ionode_id, intra_image_comm ) +! call mp_bcast( wan_data,ionode_id, intra_image_comm ) + call mp_bcast( constrain_pot, ionode_id, intra_image_comm ) + return + ! + end subroutine + + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist control + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine control_checkin( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = ' control_checkin ' + integer :: i + logical :: allowed = .false. + ! + ! + do i = 1, size( calculation_allowed ) + if( trim(calculation) == calculation_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) & + call errore( sub_name, ' calculation "'// & + & trim(calculation)//'" not allowed ',1) + if( ndr < 50 ) & + call errore( sub_name,' ndr out of range ', 1 ) + if( ndw > 0 .and. ndw < 50 ) & + call errore( sub_name,' ndw out of range ', 1 ) + if( nstep < 0 ) & + call errore( sub_name,' nstep out of range ', 1 ) + if( iprint < 1 ) & + call errore( sub_name,' iprint out of range ', 1 ) + + if( prog == 'pw' ) then + if( isave > 0 ) & + call infomsg( sub_name,' isave not used in pw ' ) + else + if( isave < 1 ) & + call errore( sub_name,' isave out of range ', 1 ) + end if + + if( dt < 0.0_dp ) & + call errore( sub_name,' dt out of range ', 1 ) + if( max_seconds < 0.0_dp ) & + call errore( sub_name,' max_seconds out of range ', 1 ) + + if( ekin_conv_thr < 0.0_dp ) then + if( prog == 'pw' ) then + call infomsg( sub_name,' ekin_conv_thr not used in pw ') + else + call errore( sub_name,' ekin_conv_thr out of range ', 1 ) + end if + end if + + if( etot_conv_thr < 0.0_dp ) & + call errore( sub_name,' etot_conv_thr out of range ', 1 ) + if( forc_conv_thr < 0.0_dp ) & + call errore( sub_name,' forc_conv_thr out of range ', 1 ) + if( prog == 'cp' ) then + if( dipfield ) & + call infomsg( sub_name,' dipfield not yet implemented ') + if( lberry ) & + call infomsg( sub_name,' lberry not implemented yet ') + if( gdir /= 0 ) & + call infomsg( sub_name,' gdir not used ') + if( nppstr /= 0 ) & + call infomsg( sub_name,' nppstr not used ') + end if + ! + if( prog == 'pw' .and. trim( restart_mode ) == 'reset_counters' ) then + call infomsg ( sub_name, ' restart_mode == reset_counters' // & + & ' not implemented in pw ' ) + end if + ! + if( refg < 0 ) & + call errore( sub_name, ' wrong table interval refg ', 1 ) + ! + if( ( prog == 'cp' ) .and. ( trim(memory) == 'small' ) .and. wf_collect ) & + call errore( sub_name, ' wf_collect = .true. is not allowed with memory = small ', 1 ) + + allowed = .false. + do i = 1, size( memory_allowed ) + if( trim(memory) == memory_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) & + call errore(sub_name, ' memory "' // trim(memory)//'" not allowed',1) + ! tb + if ( gate .and. tefield .and. (.not. dipfield) ) & + call errore(sub_name, ' gate cannot be used with tefield if dipole correction is not active', 1) + if ( gate .and. dipfield .and. (.not. tefield) ) & + call errore(sub_name, ' dipole correction is not active if tefield = .false.', 1) + + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist system + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine system_checkin( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = ' system_checkin ' + integer :: i + logical :: allowed + ! + ! + if( ( ibrav /= 0 ) .and. (celldm(1) == 0.0_dp) .and. ( a == 0.0_dp ) ) & + call errore( ' iosys ', & + & ' invalid lattice parameters ( celldm or a )', 1 ) + ! + if( nat < 0 ) & + call errore( sub_name ,' nat less than zero ', max( nat, 1) ) + ! + if( ntyp < 0 ) & + call errore( sub_name ,' ntyp less than zero ', max( ntyp, 1) ) + if( ntyp < 0 .or. ntyp > nsx ) & + call errore( sub_name , & + & ' ntyp too large, increase nsx ', max( ntyp, 1) ) + ! + if( nspin < 1 .or. nspin > 4 .or. nspin == 3 ) & + call errore( sub_name ,' nspin out of range ', max(nspin, 1 ) ) + ! + if( ecutwfc < 0.0_dp ) & + call errore( sub_name ,' ecutwfc out of range ',1) + if( ecutrho < 0.0_dp ) & + call errore( sub_name ,' ecutrho out of range ',1) + ! + if( prog == 'cp' ) then + if( degauss /= 0.0_dp ) & + call infomsg( sub_name ,' degauss is not used in cp ') + end if + ! + if( ecfixed < 0.0_dp ) & + call errore( sub_name ,' ecfixed out of range ',1) + if( qcutz < 0.0_dp ) & + call errore( sub_name ,' qcutz out of range ',1) + if( q2sigma < 0.0_dp ) & + call errore( sub_name ,' q2sigma out of range ',1) + if( prog == 'cp' ) then + if( any(starting_magnetization /= sm_not_set ) ) & + call infomsg( sub_name ,& + & ' starting_magnetization is not used in cp ') + if( la2f ) & + call infomsg( sub_name ,' la2f is not used in cp ') + if( any(hubbard_alpha /= 0.0_dp) ) & + call infomsg( sub_name ,' hubbard_alpha is not used in cp ') + if( nosym ) & + call infomsg( sub_name ,' nosym not implemented in cp ') + if( nosym_evc ) & + call infomsg( sub_name ,' nosym_evc not implemented in cp ') + if( noinv ) & + call infomsg( sub_name ,' noinv not implemented in cp ') + end if + ! + ! ... control on sic variables + ! + if ( sic /= 'none' ) then + ! + if (sic_epsilon > 1.0_dp ) & + call errore( sub_name, & + & ' invalid sic_epsilon, greater than 1.',1 ) + if (sic_epsilon < 0.0_dp ) & + call errore( sub_name, & + & ' invalid sic_epsilon, less than 0 ',1 ) + if (sic_alpha > 1.0_dp ) & + call errore( sub_name, & + & ' invalid sic_alpha, greater than 1.',1 ) + if (sic_alpha < 0.0_dp ) & + call errore( sub_name, & + & ' invalid sic_alpha, less than 0 ',1 ) + ! + if ( .not. force_pairing ) & + call errore( sub_name, & + & ' invalid force_pairing with sic activated', 1 ) + if ( nspin /= 2 ) & + call errore( sub_name, & + & ' invalid nspin with sic activated', 1 ) + if ( tot_magnetization /= 1._dp ) & + call errore( sub_name, & + & ' invalid tot_magnetization_ with sic activated', 1 ) + ! + endif + ! + ! ... control on exx variables + ! + do i = 1, size( exxdiv_treatment_allowed ) + if( trim(exxdiv_treatment) == exxdiv_treatment_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) call errore(sub_name, & + ' invalid exxdiv_treatment: '//trim(exxdiv_treatment), 1 ) + ! + if ( trim(exxdiv_treatment) == "yukawa" .and. yukawa <= 0.0 ) & + call errore(sub_name, ' invalid value for yukawa', 1 ) + ! + if ( trim(exxdiv_treatment) == "vcut_ws" .and. ecutvcut <= 0.0 ) & + call errore(sub_name, ' invalid value for ecutvcut', 1 ) + ! + if ( x_gamma_extrapolation .and. ( trim(exxdiv_treatment) == "vcut_ws" .or. & + trim(exxdiv_treatment) == "vcut_spherical" ) ) & + call errore(sub_name, ' x_gamma_extrapolation cannot be used with vcut', 1 ) + ! + ! tb - gate check + ! + if ( gate .and. tot_charge == 0 ) & + call errore(sub_name, ' charged plane (gate) to compensate tot_charge of 0', 1) + return + ! + ! ... control on fcp variables + ! + allowed = .false. + do i = 1, size(fcp_relax_allowed) + if( trim(fcp_relax) == fcp_relax_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) & + call errore(sub_name, ' fcp_relax '''//trim(fcp_relax)//''' not allowed ', 1) + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist electrons + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine electrons_checkin( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = ' electrons_checkin ' + integer :: i + logical :: allowed = .false. + ! + ! + do i = 1, size(electron_dynamics_allowed) + if( trim(electron_dynamics) == & + electron_dynamics_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) & + call errore( sub_name, ' electron_dynamics "'//& + & trim(electron_dynamics)//'" not allowed ',1) + if( emass <= 0.0_dp ) & + call errore( sub_name, ' emass less or equal 0 ',1) + if( emass_cutoff <= 0.0_dp ) & + call errore( sub_name, ' emass_cutoff less or equal 0 ',1) + if( ortho_eps <= 0.0_dp ) & + call errore( sub_name, ' ortho_eps less or equal 0 ',1) + if( ortho_max < 1 ) & + call errore( sub_name, ' ortho_max less than 1 ',1) + if( fnosee <= 0.0_dp ) & + call errore( sub_name, ' fnosee less or equal 0 ',1) + if( ekincw <= 0.0_dp ) & + call errore( sub_name, ' ekincw less or equal 0 ',1) + if( occupation_constraints ) & + call errore( sub_name, ' occupation_constraints not yet implemented ',1) + +! + return + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist ions + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine ions_checkin( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = ' ions_checkin ' + integer :: i + logical :: allowed = .false. + ! + ! + allowed = .false. + do i = 1, size(ion_dynamics_allowed) + if( trim(ion_dynamics) == ion_dynamics_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) & + call errore( sub_name, ' ion_dynamics "'// & + & trim(ion_dynamics)//'" not allowed ',1) + if( tempw <= 0.0_dp ) & + call errore( sub_name,' tempw out of range ',1) + if( fnosep( 1 ) <= 0.0_dp ) & + call errore( sub_name,' fnosep out of range ',1) + if( nhpcl > nhclm ) & + call infomsg ( sub_name,' nhpcl should be less than nhclm') + if( nhpcl < 0 ) & + call infomsg ( sub_name,' nhpcl out of range ') + if( ion_nstepe <= 0 ) & + call errore( sub_name,' ion_nstepe out of range ',1) + if( ion_maxstep < 0 ) & + call errore( sub_name,' ion_maxstep out of range ',1) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist cell + ! + !=----------------------------------------------------------------------=! + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine cell_checkin( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = ' cell_checkin ' + integer :: i + logical :: allowed = .false. + ! + ! + do i = 1, size(cell_dynamics_allowed) + if( trim(cell_dynamics) == & + cell_dynamics_allowed(i) ) allowed = .true. + end do + if( .not. allowed ) & + call errore( sub_name, ' cell_dynamics "'// & + trim(cell_dynamics)//'" not allowed ',1) + if( wmass < 0.0_dp ) & + call errore( sub_name,' wmass out of range ',1) + if( prog == 'cp' ) then + if( cell_factor /= 0.0_dp ) & + call infomsg( sub_name,' cell_factor not used in cp ') + end if + if( cell_nstepe <= 0 ) & + call errore( sub_name,' cell_nstepe out of range ',1) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist wannier + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine wannier_checkin( prog ) + !----------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = 'wannier_checkin' + ! + if ( calwf < 1 .or. calwf > 5 ) & + call errore( sub_name, ' calwf out of range ', 1 ) + ! + if ( wfsd < 1 .or. wfsd > 3 ) & + call errore( sub_name, ' wfsd out of range ', 1 ) ! + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! check input values for namelist wannier_new + ! + !=----------------------------------------------------------------------=! + ! + !---------------------------------------------------------------------- + subroutine wannier_ac_checkin( prog ) + !-------------------------------------------------------------------- + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = 'wannier_new_checkin' + ! + ! + if ( nwan > nwanx ) & + call errore( sub_name, ' nwan out of range ', 1 ) + + if ( plot_wan_num < 0 .or. plot_wan_num > nwan ) & + call errore( sub_name, ' plot_wan_num out of range ', 1 ) + + if ( plot_wan_spin < 0 .or. plot_wan_spin > 2 ) & + call errore( sub_name, ' plot_wan_spin out of range ', 1 ) + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! set values according to the "calculation" variable + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine fixval( prog ) + !----------------------------------------------------------------------- + ! + use constants, only : e2 + ! + implicit none + ! + character(len=2) :: prog ! ... specify the calling program + character(len=20) :: sub_name = ' fixval ' + ! + ! + select case( trim( calculation ) ) + case ('scf', 'ensemble') + if( prog == 'cp' ) then + electron_dynamics = 'damp' + ion_dynamics = 'none' + cell_dynamics = 'none' + end if + case ('nscf', 'bands') + if( prog == 'cp' ) occupations = 'bogus' + if( prog == 'cp' ) electron_dynamics = 'damp' + case ( 'cp-wf' ) + if( prog == 'cp' ) then + electron_dynamics = 'damp' + ion_dynamics = 'damp' + end if + if ( prog == 'pw' ) & + call errore( sub_name, ' calculation ' // & + & trim( calculation ) // ' not implemented ', 1 ) + case ( 'vc-cp-wf' ) + if( prog == 'cp' ) then + electron_dynamics = 'verlet' + ion_dynamics = 'verlet' + cell_dynamics = 'pr' + else if( prog == 'pw' ) then + call errore( sub_name, ' calculation ' // & + & trim( calculation ) // ' not implemented ', 1 ) + end if + ! +!========================================================================= +!lingzhu kong + case ( 'cp-wf-nscf' ) + if( prog == 'cp' ) then + occupations = 'fixed' + electron_dynamics = 'damp' + ion_dynamics = 'damp' + end if + if ( prog == 'pw' ) & + call errore( sub_name, ' calculation ' // & + & trim( calculation ) // ' not implemented ', 1 ) +!========================================================================= + case ('relax') + if( prog == 'cp' ) then + electron_dynamics = 'damp' + ion_dynamics = 'damp' + else if( prog == 'pw' ) then + ion_dynamics = 'bfgs' + end if + case ( 'md', 'cp' ) + if( prog == 'cp' ) then + electron_dynamics = 'verlet' + ion_dynamics = 'verlet' + else if( prog == 'pw' ) then + ion_dynamics = 'verlet' + end if + case ('vc-relax') + if( prog == 'cp' ) then + electron_dynamics = 'damp' + ion_dynamics = 'damp' + cell_dynamics = 'damp-pr' + else if( prog == 'pw' ) then + ion_dynamics = 'bfgs' + cell_dynamics= 'bfgs' + end if + case ( 'vc-md', 'vc-cp' ) + if( prog == 'cp' ) then + electron_dynamics = 'verlet' + ion_dynamics = 'verlet' + cell_dynamics = 'pr' + else if( prog == 'pw' ) then + ion_dynamics = 'beeman' + end if + ! + case default + ! + call errore( sub_name,' calculation '// & + & trim(calculation)//' not implemented ', 1 ) + ! + end select + ! + if ( prog == 'pw' ) then + ! + if ( calculation == 'nscf' .or. calculation == 'bands' ) then + ! + startingpot = 'file' + startingwfc = 'atomic+random' + ! + else if ( restart_mode == "from_scratch" ) then + ! + startingwfc = 'atomic+random' + startingpot = 'atomic' + ! + else + ! + startingwfc = 'file' + startingpot = 'file' + ! + end if + ! + else if ( prog == 'cp' ) then + ! + startingwfc = 'random' + startingpot = ' ' + ! + end if + ! + if ( trim( sic ) /= 'none' ) then + force_pairing = ( nspin == 2 .and. ( tot_magnetization==0._dp .or. & + tot_magnetization==1._dp ) ) + end if + ! + return + ! + end subroutine + ! + !=----------------------------------------------------------------------=! + ! + ! namelist parsing main routine + ! + !=----------------------------------------------------------------------=! + ! + !----------------------------------------------------------------------- + subroutine read_namelists( prog_, unit ) + !----------------------------------------------------------------------- + ! + ! this routine reads data from standard input and puts them into + ! module-scope variables (accessible from other routines by including + ! this module, or the one that contains them) + ! ---------------------------------------------- + ! + ! ... declare modules + ! + use io_global, only : ionode, ionode_id + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + ! + ! ... declare variables + ! + character(len=*) :: prog_ ! specifies the calling program, allowed: + ! prog = 'pw' pwscf + ! prog = 'cp' cp + ! prog = 'pw+ipi' pwscf + i-pi + ! + integer, intent(in), optional :: unit + ! + ! ... declare other variables + ! + character(len=2) :: prog + integer :: ios + integer :: unit_loc=5 + ! + ! ... end of declarations + ! + ! ---------------------------------------------- + ! + if(present(unit)) unit_loc = unit + ! + prog = prog_(1:2) ! allowed: 'pw' or 'cp' + if( prog /= 'pw' .and. prog /= 'cp' ) & + call errore( ' read_namelists ', ' unknown calling program ', 1 ) + ! + ! ... default settings for all namelists + ! + call control_defaults( prog ) + call system_defaults( prog ) + call electrons_defaults( prog ) + call ions_defaults( prog ) + call cell_defaults( prog ) + ! + ! ... here start reading standard input file + ! + ! + ! ... control namelist + ! + ios = 0 + if( ionode ) then + read( unit_loc, control, iostat = ios ) + end if + call check_namelist_read(ios, unit_loc, "control") + ! + call control_bcast( ) + call control_checkin( prog ) + ! + ! ... fixval changes some default values according to the value + ! ... of "calculation" read in control namelist + ! + call fixval( prog ) + ! + ! ... system namelist + ! + ios = 0 + if( ionode ) then + read( unit_loc, system, iostat = ios ) + end if + call check_namelist_read(ios, unit_loc, "system") + ! + call system_bcast( ) + ! + call system_checkin( prog ) + ! + ! ... electrons namelist + ! + ios = 0 + if( ionode ) then + read( unit_loc, electrons, iostat = ios ) + end if + call check_namelist_read(ios, unit_loc, "electrons") + ! + call electrons_bcast( ) + call electrons_checkin( prog ) + ! + ! ... ions namelist - must be read only if ionic motion is expected, + ! ... or if code called by i-pi via run_driver + ! + ios = 0 + if ( ionode ) then + if ( ( trim( calculation ) /= 'nscf' .and. & + trim( calculation ) /= 'bands' ) .or. & + ( trim( prog_ ) == 'pw+ipi' ) ) then + read( unit_loc, ions, iostat = ios ) + end if + ! + ! scf might (optionally) have &ions :: ion_positions = 'from_file' + ! + if ( (ios /= 0) .and. trim( calculation ) == 'scf' ) then + ! presumably, not found: rewind the file pointer to the location + ! of the previous present section, in this case electrons + rewind( unit_loc ) + read( unit_loc, electrons, iostat = ios ) + end if + ! + end if + ! + call check_namelist_read(ios, unit_loc, "ions") + ! + call ions_bcast( ) + call ions_checkin( prog ) + ! + ! ... cell namelist + ! + ios = 0 + if( ionode ) then + if( trim( calculation ) == 'vc-relax' .or. & + trim( calculation ) == 'vc-cp' .or. & + trim( calculation ) == 'vc-md' .or. & + trim( calculation ) == 'vc-md' .or. & + trim( calculation ) == 'vc-cp-wf') then + read( unit_loc, cell, iostat = ios ) + end if + end if + call check_namelist_read(ios, unit_loc, "cell") + ! + call cell_bcast() + call cell_checkin( prog ) + ! + ios = 0 + if( ionode ) then + if (tabps) then + read( unit_loc, press_ai, iostat = ios ) + end if + end if + call check_namelist_read(ios, unit_loc, "press_ai") + ! + call press_ai_bcast() + ! + ! ... wannier namelist + ! + call wannier_defaults( prog ) + ios = 0 + if( ionode ) then + if( trim( calculation ) == 'cp-wf' .or. & + trim( calculation ) == 'vc-cp-wf' .or. & + trim( calculation ) == 'cp-wf-nscf') then + read( unit_loc, wannier, iostat = ios ) + end if + end if + call check_namelist_read(ios, unit_loc, "wannier") + ! + call wannier_bcast() + call wannier_checkin( prog ) + ! + ! ... wannier_new namelist + ! + call wannier_ac_defaults( prog ) + ios = 0 + if( ionode ) then + if( use_wannier ) then + read( unit_loc, wannier_ac, iostat = ios ) + end if + end if + call check_namelist_read(ios, unit_loc, "wannier_ac") + ! + call wannier_ac_bcast() + call wannier_ac_checkin( prog ) + ! + return + ! + end subroutine read_namelists + ! + subroutine check_namelist_read(ios, unit_loc, nl_name) + use io_global, only : ionode, ionode_id + use mp, only : mp_bcast + use mp_images, only : intra_image_comm + ! + implicit none + integer,intent(in) :: ios, unit_loc + character(len=*) :: nl_name + character(len=512) :: line + integer :: ios2 + ! + if( ionode ) then + ios2=0 + if (ios /=0) then + backspace(unit_loc) + read(unit_loc,'(a512)', iostat=ios2) line + end if + end if + + call mp_bcast( ios2, ionode_id, intra_image_comm ) + if( ios2 /= 0 ) then + call errore( ' read_namelists ', ' could not find namelist &'//trim(nl_name), 2) + endif + ! + call mp_bcast( ios, ionode_id, intra_image_comm ) + call mp_bcast( line, ionode_id, intra_image_comm ) + if( ios /= 0 ) then + call errore( ' read_namelists ', & + ' bad line in namelist &'//trim(nl_name)//& + ': "'//trim(line)//'" (error could be in the previous line)',& + 1 ) + end if + ! + end subroutine check_namelist_read + ! +end module read_namelists_module