Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : #include "abi_common.h"
7 :
8 : module m_tdep_dataset
9 :
10 : use defs_basis
11 : use m_errors
12 : use m_abicore
13 : use m_xmpi
14 : use m_abihist
15 : use m_parser
16 : use netcdf
17 : use m_nctk
18 : use m_fstrings, only : inupper,ljust,next_token
19 : use m_io_tools, only : file_exists
20 : use m_ddb, only : ddb_type
21 : use m_ddb_hdr, only : ddb_hdr_type
22 : use m_crystal, only : crystal_t
23 :
24 : implicit none
25 :
26 : character(len=6),public,parameter :: version_string = ' 4.0'
27 :
28 : type atdep_dataset_type
29 :
30 : integer :: natom
31 : ! Number of atoms in the supercell
32 :
33 : integer :: natom_unitcell
34 : ! Number of atoms in the unit cell
35 :
36 : integer :: nshell_max
37 : integer :: nstep_max
38 : integer :: nstep_min
39 : integer :: nstep_tot
40 : integer :: my_nstep
41 : integer :: ntypat
42 : integer :: use_ideal_positions
43 : integer :: stdout
44 : integer :: stdlog
45 : integer :: bzpath
46 : integer :: order
47 : integer :: slice
48 : integer :: enunit
49 : integer :: readifc
50 : integer :: together
51 : integer :: alloy
52 : integer :: ityp_alloy1
53 : integer :: ityp_alloy2
54 : integer :: nproc(2)
55 : integer :: bzlength
56 : integer :: ngqpt1(3)
57 : integer :: ngqpt2(3)
58 : integer :: bravais(11)
59 : integer :: use_weights
60 : integer :: prtphdos
61 : integer :: symbec ! Symmetrize born effective charges
62 : integer, allocatable :: typat_unitcell(:)
63 : integer, allocatable :: typat(:)
64 : integer, allocatable :: lgth_segments(:)
65 : logical :: debug
66 : logical :: loto ! TODO rename dipdip
67 : logical :: netcdf
68 : double precision :: angle_alpha
69 : double precision :: dosdeltae
70 : double precision :: rcut
71 : double precision :: rcut3
72 : double precision :: rcut4
73 : double precision :: temperature
74 : double precision :: tolread
75 : double precision :: tolinbox
76 : double precision :: tolmatch
77 : double precision :: tolmotif
78 :
79 : double precision :: rprimd_md(3,3)
80 : ! Dimensioned primitive vectors
81 :
82 : double precision :: multiplicity(3,3)
83 : ! Supercell scaling matrix (integers)
84 :
85 : double precision :: dielt(3,3)
86 : ! Dielectric tensor
87 :
88 : double precision, allocatable :: amu(:)
89 : ! amu(natom_unitcell)
90 : ! Mass of the atoms (atomic mass unit)
91 :
92 : double precision, allocatable :: zeff(:,:,:)
93 : ! zeff(3,3,natom_unitcell)
94 : ! Born effective charge tensor
95 :
96 : double precision, allocatable :: znucl(:)
97 : double precision, allocatable :: qpt(:,:)
98 : double precision, allocatable :: xred_unitcell(:,:)
99 :
100 : character (len=2), allocatable :: special_qpt(:)
101 : character (len=fnlen) :: output_prefix
102 : character (len=fnlen) :: input_prefix
103 : character (len=fnlen) :: output_file
104 :
105 : end type atdep_dataset_type
106 :
107 : type MPI_enreg_type
108 :
109 : integer :: comm_shell
110 : integer :: comm_step
111 : integer :: comm_shellstep
112 : integer :: nproc
113 : integer :: nproc_shell
114 : integer :: nproc_step
115 : integer :: master
116 : integer :: me_shell
117 : integer :: me_step
118 : integer, allocatable :: my_nshell(:)
119 : integer :: my_nstep
120 : logical :: iam_master
121 : integer, allocatable :: nstep_acc(:)
122 : integer, allocatable :: nstep_all(:)
123 : integer, allocatable :: shft_step(:)
124 : logical, allocatable :: my_shell(:)
125 : logical, allocatable :: my_step(:)
126 :
127 : end type MPI_enreg_type
128 :
129 : public :: tdep_read_input
130 : public :: tdep_destroy_invar
131 : public :: tdep_init_MPIdata
132 : public :: tdep_destroy_mpidata
133 :
134 : contains
135 :
136 : !====================================================================================================
137 :
138 44 : subroutine tdep_read_input(input_path,Hist,Invar)
139 :
140 : ! Arguments-------------------------------
141 : character(len=*), intent(in):: input_path
142 : type(atdep_dataset_type),intent(out) :: Invar
143 : type(abihist), intent(out) :: Hist
144 :
145 : ! Local variables-------------------------
146 : ! scalars
147 : integer :: values(8)
148 : integer :: ncid, ncerr, me,ierr,master
149 : integer :: nimage, mdtime, natom_id,nimage_id,time_id,xyz_id,six_id
150 : integer :: ntypat_id,iatcell
151 : integer :: ii,jj,shift,iatom,itypat,sum_alloy1,sum_alloy2
152 : integer:: lenstr, marr, jdtset, tread
153 : logical :: has_nimage
154 : double precision :: dtion,amu_average,born_average
155 : character (len=8) :: date
156 : character (len=10) :: time
157 : character (len=5) :: zone
158 : character(len=500) :: msg
159 : character(len=fnlen) :: input_filename,hist_filename,ddb_filename
160 : character(len=strlen):: string, raw_string
161 2288 : type(crystal_t):: crystal
162 44 : type(ddb_type):: ddb
163 2288 : type(ddb_hdr_type):: ddb_hdr
164 : ! arrays
165 : character(len=3),parameter :: month_names(12)=(/'Jan','Feb','Mar','Apr','May','Jun',&
166 : & 'Jul','Aug','Sep','Oct','Nov','Dec'/)
167 44 : integer, allocatable:: intarr(:)
168 44 : integer, allocatable :: typat_unitcell_tmp(:)
169 : real(dp) :: zeff_average(3,3)
170 44 : real(dp), allocatable:: zeff_tmp(:,:,:)
171 88 : real(dp), allocatable :: xred_unitcell_tmp(:,:),amu_tmp(:),born_charge_tmp(:),znucl_tmp(:)
172 44 : real(dp), allocatable:: dprarr(:)
173 :
174 : ! *********************************************************************
175 :
176 : ! Define output files
177 44 : Invar%stdout=ab_out
178 44 : Invar%stdlog=std_out
179 :
180 : ! Define default values
181 44 : Invar%nshell_max=500
182 44 : Invar%angle_alpha=90.d0
183 44 : Invar%bzpath=0
184 44 : Invar%order=2
185 44 : Invar%slice=1
186 44 : Invar%enunit=0
187 44 : Invar%together=1
188 44 : Invar%alloy=0
189 44 : Invar%ityp_alloy1=0
190 44 : Invar%ityp_alloy2=0
191 132 : Invar%nproc(:)=1
192 44 : Invar%bzlength=0
193 44 : Invar%tolread=1.d-8
194 44 : Invar%tolmotif=5.d-2
195 44 : Invar%tolinbox=5.d-2
196 44 : Invar%tolmatch=5.d-2
197 44 : Invar%dosdeltae=0.2_dp/Ha_cmm1
198 44 : Invar%debug=.false.
199 44 : Invar%loto=.false.
200 44 : Invar%netcdf=.false.
201 44 : Invar%use_ideal_positions=0
202 44 : Invar%use_weights=0
203 44 : Invar%prtphdos=1
204 44 : Invar%symbec=0
205 : ! In order to have an accuracy better than 1meV
206 176 : Invar%ngqpt1(:)=8
207 176 : Invar%ngqpt2(:)=32
208 572 : Invar%dielt(:,:)=zero; Invar%dielt(1,1)=one; Invar%dielt(2,2)=one; Invar%dielt(3,3)=one
209 :
210 44 : master = 0
211 44 : me = xmpi_comm_rank(xmpi_world)
212 44 : if (me==master) then
213 :
214 44 : if (len_trim(input_path) == 0) then
215 :
216 0 : write(std_out, "(2a)")" DeprecationWarning: ",ch10
217 0 : write(std_out, "(a)") " The files file has been deprecated in Abinit10 and will be removed in Abinit11."
218 0 : write(std_out, "(2a)")" Use the syntax `atdep t01.abi` where t01.abi is an atdep input,",ch10
219 0 : write(std_out, "(2a)")" and use input variables output_file, indata_prefix, outdata_prefix.",ch10
220 :
221 0 : write(Invar%stdlog,'(a)',err=10) ' Give name for input file '
222 0 : read(*, '(a)',err=10) input_filename
223 0 : if ( input_filename == "" ) input_filename='input.in'
224 0 : write(Invar%stdlog, '(a)',err=10) '.'//trim(input_filename)
225 : 10 continue
226 : ! Check if a NetCDF file is available
227 0 : write(Invar%stdlog,'(a)',err=11) ' Give root name for generic input files (NetCDF or ASCII)'
228 0 : read(*, '(a)',err=11) Invar%input_prefix
229 0 : write(Invar%stdlog, '(a)',err=11) '.'//trim(Invar%input_prefix)
230 : 11 continue
231 0 : write(Invar%stdlog,'(a)', err=12)' Give root name for generic output files:'
232 0 : read (*, '(a)', err=12) Invar%output_prefix
233 0 : if ( Invar%output_prefix == "" ) Invar%output_prefix = 'atdep'
234 0 : write (Invar%stdlog, '(a)', err=12 ) '.'//trim(Invar%output_prefix)
235 : 12 continue
236 0 : Invar%output_file = trim(Invar%output_prefix)//'.abo'
237 : else
238 44 : input_filename = input_path
239 :
240 : ! Read input
241 44 : string = repeat(" ", strlen)
242 44 : raw_string = repeat(" ", strlen)
243 44 : call instrng(input_filename, lenstr, 1, strlen, string, raw_string)
244 : ! To make case-insensitive, map characters to upper case.
245 44 : call inupper(string(1:lenstr))
246 :
247 44 : marr = 3
248 44 : ABI_MALLOC(intarr, (marr))
249 44 : ABI_MALLOC(dprarr, (marr))
250 44 : jdtset = 0
251 :
252 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), &
253 44 : "indata_prefix", tread, 'KEY', key_value=Invar%input_prefix)
254 44 : if (tread == 0) then
255 0 : Invar%input_prefix = ''
256 : end if
257 44 : write(Invar%stdlog, "(2a)")"- Root name for input files: ", trim(Invar%input_prefix)
258 :
259 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), &
260 44 : "outdata_prefix", tread, 'KEY', key_value=Invar%output_prefix)
261 44 : if (tread == 0) then
262 0 : Invar%output_prefix = 'atdep'
263 : end if
264 44 : write(Invar%stdlog, "(2a)")"- Root name for output files: ", trim(Invar%output_prefix)
265 :
266 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), &
267 44 : "output_file", tread, 'KEY', key_value=Invar%output_file)
268 44 : if (tread == 0) then
269 0 : Invar%output_file = trim(Invar%output_prefix) // '.abo'
270 : end if
271 44 : write(Invar%stdlog, "(2a)")"- Main output file: ", trim(Invar%output_file)
272 :
273 44 : ABI_FREE(intarr)
274 44 : ABI_FREE(dprarr)
275 :
276 : end if
277 :
278 : ! Setup other filenames
279 44 : if ( Invar%input_prefix == "" ) then
280 0 : hist_filename='HIST.nc'
281 0 : ddb_filename='DDB'
282 : else
283 44 : hist_filename=trim(Invar%input_prefix)//'_HIST.nc'
284 44 : ddb_filename=trim(Invar%input_prefix)//'_DDB'
285 : end if
286 :
287 44 : open(unit=Invar%stdout,file=trim(Invar%output_file))
288 :
289 : end if !me
290 :
291 44 : call xmpi_bcast(input_filename,master,xmpi_world,ierr)
292 44 : call xmpi_bcast(hist_filename,master,xmpi_world,ierr)
293 44 : call xmpi_bcast(ddb_filename,master,xmpi_world,ierr)
294 44 : call xmpi_bcast(Invar%output_prefix,master,xmpi_world,ierr)
295 44 : call xmpi_bcast(Invar%input_prefix,master,xmpi_world,ierr)
296 :
297 : #if defined HAVE_NETCDF
298 : !Open netCDF file
299 44 : ncerr=nf90_open(path=trim(hist_filename),mode=NF90_NOWRITE,ncid=ncid)
300 44 : if(ncerr /= NF90_NOERR) then
301 38 : write(Invar%stdlog,'(3a)') '-'//'Could not open ',trim(hist_filename),', starting from scratch'
302 38 : Invar%netcdf=.false.
303 : else
304 6 : write(Invar%stdlog,'(3a)') '-'//'Succesfully open ',trim(hist_filename),' for reading'
305 6 : write(Invar%stdlog,'(a)') ' Extracting information from NetCDF file...'
306 6 : Invar%netcdf=.true.
307 : end if
308 :
309 44 : if ( Invar%netcdf) then
310 : call get_dims_hist(ncid,Invar%natom,Invar%ntypat,nimage,mdtime,&
311 6 : & natom_id,ntypat_id,nimage_id,time_id,xyz_id,six_id,has_nimage)
312 26 : ABI_CALLOC(Invar%amu,(Invar%ntypat))
313 682 : ABI_CALLOC(Invar%typat,(Invar%natom))
314 26 : ABI_CALLOC(Invar%znucl,(Invar%ntypat))
315 6 : call read_csts_hist(ncid,dtion,Invar%typat,Invar%znucl,Invar%amu)
316 :
317 : ! Need to close NetCDF file because it is going to be reopened by read_md_hist
318 6 : ncerr = nf90_close(ncid)
319 : ! .false. -> Velocities note used
320 : ! .true. -> acell and rprimd may change (2017_04 only NVT/isoK used but maybe
321 : ! .false. -> read all times
322 : ! NPT one day ?)
323 18 : call read_md_hist(hist_filename,Hist,.false.,.true.,.false.)
324 : end if
325 : #endif
326 :
327 : ! =========================================================================== !
328 : ! Read input file
329 :
330 44 : string = repeat(" ", strlen)
331 44 : raw_string = repeat(" ", strlen)
332 44 : call instrng(input_filename, lenstr, 1, strlen, string, raw_string)
333 : ! To make case-insensitive, map characters to upper case.
334 44 : call inupper(string(1:lenstr))
335 :
336 : ! marr is the aximum array size. It is thus a hard-coded maximum value
337 : ! for 3 * (number of atoms)
338 44 : marr = 9000
339 44 : ABI_MALLOC(intarr, (marr))
340 44 : ABI_MALLOC(dprarr, (marr))
341 44 : jdtset = 0
342 :
343 : ! Mandatory input variables
344 : ! -------------------------
345 :
346 : ! Bravais lattice
347 44 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'brav', tread, 'INT')
348 44 : if (tread == 0) then
349 : write(msg,*)&
350 0 : 'Variable "brav" is mandatory, but was not found in input file.'
351 0 : ABI_ERROR(msg)
352 : end if
353 132 : Invar%bravais(1:2) = intarr(1:2)
354 :
355 : ! Angle alpha
356 44 : if ((Invar%bravais(1).eq.2).or.(Invar%bravais(1).eq.5)) then
357 4 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'angle', tread, 'DPR')
358 4 : if (tread == 0) then
359 : write(msg,*)&
360 0 : 'Variable "angle" is mandatory for this bravais lattice,',ch10,&
361 0 : 'but was not found in input file.'
362 0 : ABI_ERROR(msg)
363 : end if
364 4 : Invar%angle_alpha = dprarr(1)
365 : else
366 40 : Invar%angle_alpha=90.d0
367 : end if
368 :
369 : ! natom_unitcell
370 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'natom_unitcell', tread, 'INT')
371 44 : if (tread == 0) then
372 : write(msg,*)&
373 0 : 'Variable "natom_unitcell" is mandatory, but was not found in input file.'
374 0 : ABI_ERROR(msg)
375 : end if
376 44 : Invar%natom_unitcell = intarr(1)
377 44 : if (3*Invar%natom_unitcell.gt.marr) then
378 : write(msg,*)&
379 0 : 'Maximum number of atoms exceeded. Modify source code to circumvent problem.'
380 0 : ABI_ERROR(msg)
381 : end if
382 :
383 : ! xred_unitcell
384 560 : ABI_CALLOC(Invar%xred_unitcell,(3,Invar%natom_unitcell))
385 44 : call intagm(dprarr, intarr, jdtset, marr, 3*Invar%natom_unitcell, string(1:lenstr), 'xred_unitcell', tread, 'DPR')
386 44 : if (tread == 0) then
387 : write(msg,*)&
388 0 : 'Variable "xred_unitcell" is mandatory, but was not found in input file.'
389 0 : ABI_ERROR(msg)
390 : end if
391 560 : Invar%xred_unitcell(:,:) = reshape(dprarr(1:3*Invar%natom_unitcell),(/3,Invar%natom_unitcell/))
392 176 : do ii=1,3
393 497 : do iatcell=1,Invar%natom_unitcell
394 453 : if ((Invar%xred_unitcell(ii,iatcell).le.(-0.5)).or.(Invar%xred_unitcell(ii,iatcell).gt.(0.5))) then
395 3 : do while (Invar%xred_unitcell(ii,iatcell).le.(-0.5))
396 0 : Invar%xred_unitcell(ii,iatcell)=Invar%xred_unitcell(ii,iatcell)+1.d0
397 : end do
398 6 : do while (Invar%xred_unitcell(ii,iatcell).gt.(0.5))
399 6 : Invar%xred_unitcell(ii,iatcell)=Invar%xred_unitcell(ii,iatcell)-1.d0
400 : end do
401 : end if
402 : end do
403 : end do
404 :
405 : ! typat_unitcell
406 239 : ABI_CALLOC(Invar%typat_unitcell,(Invar%natom_unitcell))
407 44 : call intagm(dprarr, intarr, jdtset, marr, Invar%natom_unitcell, string(1:lenstr), 'typat_unitcell', tread, 'INT')
408 44 : if (tread == 0) then
409 : write(msg,*)&
410 0 : 'Variable "typat_unitcell" is mandatory, but was not found in input file.'
411 0 : ABI_ERROR(msg)
412 : end if
413 151 : Invar%typat_unitcell(:) = intarr(1:Invar%natom_unitcell)
414 :
415 44 : if (Invar%netcdf) then
416 78 : Invar%rprimd_md(:,:)=TRANSPOSE(Hist%rprimd(:,:,Hist%ihist))
417 : else
418 : ! ntypat
419 38 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'ntypat', tread, 'INT')
420 38 : if (tread == 1) then
421 38 : Invar%ntypat = intarr(1)
422 : else
423 : write(msg,*)&
424 0 : 'The NetCDF file .nc is not used.',ch10,&
425 0 : 'The variable "ntypat" is thus mandatory.',ch10,&
426 0 : 'ACTION : Please modify your input file'
427 0 : ABI_ERROR(msg)
428 : end if
429 :
430 : ! amu
431 162 : ABI_MALLOC(Invar%amu,(Invar%ntypat)); Invar%amu(:)=zero
432 38 : call intagm(dprarr, intarr, jdtset, marr, Invar%ntypat, string(1:lenstr), 'amu', tread, 'DPR')
433 38 : if (tread == 1) then
434 124 : Invar%amu = dprarr(1:Invar%ntypat)
435 : else
436 : write(msg,*)&
437 0 : 'The NetCDF file .nc is not used.',ch10,&
438 0 : 'The variable "amu" is thus mandatory.',ch10,&
439 0 : 'ACTION : Please modify your input file'
440 0 : ABI_ERROR(msg)
441 : end if
442 :
443 : ! rprimd_md
444 38 : call intagm(dprarr, intarr, jdtset, marr, 9, string(1:lenstr), 'rprimd', tread, 'LEN')
445 38 : if (tread == 1) then
446 494 : Invar%rprimd_md(:,:) = TRANSPOSE(reshape(dprarr(1:9),(/3, 3/)))
447 : else
448 : write(msg,*)&
449 0 : 'The NetCDF file .nc is not used.',ch10,&
450 0 : 'The variable "rprimd" is thus mandatory.',ch10,&
451 0 : 'ACTION : Please modify your input file'
452 0 : ABI_ERROR(msg)
453 : end if
454 :
455 : ! natom
456 38 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'natom', tread, 'INT')
457 38 : if (tread == 1) then
458 38 : Invar%natom = intarr(1)
459 : else
460 : write(msg,*)&
461 0 : 'The NetCDF file .nc is not used.',ch10,&
462 0 : 'The variable "natom" is thus mandatory.',ch10,&
463 0 : 'ACTION : Please modify your input file'
464 0 : ABI_ERROR(msg)
465 : end if
466 :
467 : ! typat
468 6332 : ABI_MALLOC(Invar%typat,(Invar%natom)); Invar%typat(:)=0
469 38 : call intagm(dprarr, intarr, jdtset, marr, Invar%natom, string(1:lenstr), 'typat', tread, 'INT')
470 38 : if (tread == 1) then
471 6294 : Invar%typat = intarr(1:Invar%natom)
472 : else
473 : write(msg,*)&
474 0 : 'The NetCDF file .nc is not used.',ch10,&
475 0 : 'The variable "typat" is thus mandatory.',ch10,&
476 0 : 'ACTION : Please modify your input file'
477 0 : ABI_ERROR(msg)
478 : end if
479 :
480 : end if
481 :
482 : ! multiplicity
483 44 : call intagm(dprarr, intarr, jdtset, marr, 9, string(1:lenstr), 'multiplicity', tread, 'DPR')
484 44 : if (tread == 1) then
485 572 : Invar%multiplicity(:,:) = TRANSPOSE(reshape(dprarr(1:9),(/3, 3/)))
486 : else
487 : write(msg,*)&
488 0 : 'Variable "multiplicity" is mandatory, but was not found in input file.'
489 0 : ABI_ERROR(msg)
490 : end if
491 :
492 : ! temperature
493 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'temperature', tread, 'DPR')
494 44 : if (tread == 1) then
495 44 : Invar%temperature = dprarr(1)
496 : else
497 : write(msg,*)&
498 0 : 'Variable "temperature" is mandatory, but was not found in input file.'
499 0 : ABI_ERROR(msg)
500 : end if
501 :
502 : ! nstep_max
503 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'nstep_max', tread, 'INT')
504 44 : if (tread == 1) then
505 44 : Invar%nstep_max = intarr(1)
506 : else
507 : write(msg,*)&
508 0 : 'Variable "nstep_max" is mandatory, but was not found in input file.'
509 0 : ABI_ERROR(msg)
510 : end if
511 :
512 : ! nstep_min
513 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'nstep_min', tread, 'INT')
514 44 : if (tread == 1) then
515 44 : Invar%nstep_min = intarr(1)
516 : else
517 : write(msg,*)&
518 0 : 'Variable "nstep_min" is mandatory, but was not found in input file.'
519 0 : ABI_ERROR(msg)
520 : end if
521 :
522 : ! rcut
523 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'rcut', tread, 'LEN')
524 44 : if (tread == 1) then
525 44 : Invar%rcut = dprarr(1)
526 : else
527 : write(msg,*)&
528 0 : 'Variable "rcut" is mandatory, but was not found in input file.'
529 0 : ABI_ERROR(msg)
530 : end if
531 :
532 : ! debug_mode (optional)
533 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'debug_mode', tread, 'INT')
534 44 : if (tread == 1) then
535 2 : if (intarr(1).ne.0) Invar%debug = .true.
536 : end if
537 :
538 : ! =========================================================================== !
539 : ! Read DDB file if available to retrieve the Born effective charges
540 : ! and dielectric tensor. Enforce charge neutrality by equal redistribution of excess charge
541 : ! =========================================================================== !
542 :
543 : ! symbec
544 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'symbec', tread, 'INT')
545 44 : if (tread == 1) then
546 0 : Invar%symbec = intarr(1)
547 : !write(Invar%stdout,'(1x,a20,1x,i4)') ljust('symbec',20),Invar%symbec
548 : end if
549 :
550 :
551 132 : ABI_MALLOC(Invar%zeff, (3,3,Invar%natom_unitcell))
552 1435 : Invar%zeff = zero
553 :
554 44 : if (file_exists(ddb_filename).or.file_exists(nctk_ncify(ddb_filename))) then
555 :
556 0 : call ddb%from_file(ddb_filename, ddb_hdr, crystal, xmpi_world, -1)
557 : ! TODO Check that the crystal in the DDB matches the unit cell.
558 :
559 0 : ii = ddb%get_dielt_zeff(crystal,1,1,0,Invar%dielt,Invar%zeff)
560 0 : if (ii/=0) Invar%loto = .true.
561 :
562 0 : call crystal%free()
563 0 : call ddb_hdr%free()
564 0 : call ddb%free()
565 :
566 : !! Take the average born effective charge for each atom.
567 44 : if (Invar%symbec==1) then
568 0 : do iatom=1,Invar%natom_unitcell
569 : born_average = zero
570 0 : do ii=1,3
571 0 : born_average = born_average + Invar%zeff(ii,ii,iatom)
572 : end do
573 0 : Invar%zeff(:,:,iatom) = zero
574 0 : do ii=1,3
575 0 : Invar%zeff(ii,ii,iatom) = born_average / 3
576 : end do
577 : end do
578 : end if
579 : end if
580 :
581 :
582 : ! =========================================================================== !
583 : ! Output header and mandatory input variables
584 : ! =========================================================================== !
585 :
586 : ! Write version, copyright, date...
587 44 : write(Invar%stdout,*) ' '
588 :
589 44 : if (Invar%debug) then
590 2 : write(Invar%stdout,'(a,a,a)') '.Version ', version_string,' of ATDEP (Debug)'
591 : else
592 42 : write(Invar%stdout,'(a,a,a)') '.Version ', version_string,' of ATDEP'
593 : end if
594 :
595 44 : write(Invar%stdout,'(a)') '.Copyright (C) 1998-2026 ABINIT group (FB,JB,GA).'
596 44 : write(Invar%stdout,'(a)') ' ABINIT comes with ABSOLUTELY NO WARRANTY.'
597 44 : write(Invar%stdout,'(a)') ' It is free software, and you are welcome to redistribute it'
598 44 : write(Invar%stdout,'(a)') ' under certain conditions (GNU General Public License,'
599 44 : write(Invar%stdout,'(a)') ' see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt).'
600 44 : write(Invar%stdout,*) ' '
601 44 : write(Invar%stdout,'(a)') ' ABINIT is a project of the Universite Catholique de Louvain,'
602 44 : write(Invar%stdout,'(a)') ' Corning Inc. and other collaborators, see'
603 44 : write(Invar%stdout,'(a)') ' ~abinit/doc/developers/contributors.txt .'
604 44 : write(Invar%stdout,'(a)') ' Please read https://docs.abinit.org/theory/acknowledgments for suggested'
605 44 : write(Invar%stdout,'(a)') ' acknowledgments of the ABINIT effort.'
606 44 : write(Invar%stdout,'(a)') ' For more information, see http://www.abinit.org .'
607 :
608 44 : call date_and_time(date,time,zone,values)
609 44 : write(Invar%stdout,'(/,a,i2,1x,a,1x,i4,a)') '.Starting date : ',values(3),month_names(values(2)),values(1),'.'
610 :
611 44 : write(Invar%stdout,*) ' '
612 44 : write(Invar%stdout,*) '#############################################################################'
613 44 : write(Invar%stdout,*) '######################### ECHO OF INPUT FILE ################################'
614 44 : write(Invar%stdout,*) '#############################################################################'
615 :
616 44 : write(Invar%stdout,'(a)') ' ======================= Define the unitcell ================================='
617 44 : write(Invar%stdout,'(1x,a20,1x,i4,1x,i4)') ljust('brav',20),Invar%bravais(1),Invar%bravais(2)
618 44 : if ((Invar%bravais(1).eq.2).or.(Invar%bravais(1).eq.5)) then
619 4 : write(Invar%stdout,'(1x,a20,1x,f15.10)') 'angle',Invar%angle_alpha
620 : end if
621 44 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('natom_unitcell',20),Invar%natom_unitcell
622 44 : write(Invar%stdout,'(1x,a20)') ljust('xred_unitcell',20)
623 151 : do ii=1,Invar%natom_unitcell
624 472 : write(Invar%stdout,'(22x,3(f15.10,1x))') (Invar%xred_unitcell(jj,ii), jj=1,3)
625 : end do
626 151 : write(Invar%stdout,'(1x,a20,20(1x,i4))') ljust('typat_unitcell',20),(Invar%typat_unitcell(jj),jj=1,Invar%natom_unitcell)
627 44 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('ntypat',20),Invar%ntypat
628 100 : write(Invar%stdout,'(1x,a20,20(1x,f15.10))') ljust('amu',20),(Invar%amu(jj),jj=1,Invar%ntypat)
629 :
630 : ! znucl (optional)
631 44 : if (.not.Invar%netcdf) then
632 38 : call intagm(dprarr, intarr, jdtset, marr, Invar%ntypat, string(1:lenstr), 'znucl', tread, 'DPR')
633 38 : if (tread == 1) then
634 162 : ABI_MALLOC(Invar%znucl,(Invar%ntypat)); Invar%znucl(:)=zero
635 124 : Invar%znucl = dprarr(1:Invar%ntypat)
636 86 : write(Invar%stdout,'(1x,a20,20(1x,f15.10))') ljust('znucl',20),(Invar%znucl(jj),jj=1,Invar%ntypat)
637 : end if
638 : end if
639 :
640 :
641 44 : write(Invar%stdout,'(a)') ' ======================= Define the supercell ================================'
642 44 : write(Invar%stdout,'(1x,a20)') ljust('rprimd',20)
643 176 : do ii=1,3
644 572 : write(Invar%stdout,'(22x,3(f15.10,1x))') (Invar%rprimd_md(ii,jj),jj=1,3)
645 : end do
646 44 : write(Invar%stdout,'(1x,a20)') ljust('multiplicity',20)
647 176 : do ii=1,3
648 572 : write(Invar%stdout,'(22x,3(f15.10,1x))') (Invar%multiplicity(ii,jj),jj=1,3)
649 : end do
650 44 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('natom',20),Invar%natom
651 44 : write(Invar%stdout,'(1x,a20)') ljust('typat',20)
652 750 : do ii=1,Invar%natom,10
653 750 : if (ii+9.lt.Invar%natom) then
654 7282 : write(Invar%stdout,'(22x,10(i4,1x))') (Invar%typat(ii+jj-1),jj=1,10)
655 : else
656 306 : write(Invar%stdout,'(22x,10(i4,1x))') (Invar%typat(jj),jj=ii,Invar%natom)
657 : end if
658 : end do
659 :
660 44 : write(Invar%stdout,'(a)') ' ======================= Define computational details ========================'
661 44 : write(Invar%stdout,'(1x,a20,1x,i5)') ljust('nstep_max',20),Invar%nstep_max
662 44 : write(Invar%stdout,'(1x,a20,1x,i5)') ljust('nstep_min',20),Invar%nstep_min
663 44 : write(Invar%stdout,'(1x,a20,1x,f15.10)') ljust('rcut',20),Invar%rcut
664 44 : write(Invar%stdout,'(1x,a20,1x,f15.10)') ljust('temperature',20),Invar%temperature
665 :
666 :
667 : ! =========================================================================== !
668 : ! Optional input variables
669 : ! TODO separate reading optional input variables from writing in the output.
670 :
671 44 : write(Invar%stdout,'(a)') ' ======================= Optional input variables ============================'
672 :
673 : ! prtphdos
674 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'prtphdos', tread, 'INT')
675 44 : if (tread == 1) then
676 0 : Invar%prtphdos = intarr(1)
677 0 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('prtphdos',20),Invar%use_ideal_positions
678 : end if
679 :
680 : ! dosdeltae
681 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'dosdeltae', tread, 'ENE')
682 44 : if (tread == 1) then
683 0 : Invar%dosdeltae = dprarr(1)
684 0 : write(Invar%stdout,'(1x,a20,1x,f15.10)') ljust('dosdeltae',20),Invar%dosdeltae
685 : end if
686 :
687 : ! use_ideal_positions
688 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'use_ideal_positions', tread, 'INT')
689 44 : if (tread == 1) then
690 40 : Invar%use_ideal_positions = intarr(1)
691 40 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('use_ideal_positions',20),Invar%use_ideal_positions
692 : end if
693 :
694 : ! born_charge
695 44 : call intagm(dprarr, intarr, jdtset, marr, Invar%ntypat, string(1:lenstr), 'born_charge', tread, 'DPR')
696 44 : if (tread == 1) then
697 5 : Invar%loto=.true.
698 25 : ABI_MALLOC(born_charge_tmp,(Invar%ntypat)) ; born_charge_tmp(:) = zero
699 15 : born_charge_tmp(:) = dprarr(1:Invar%ntypat)
700 148 : Invar%zeff(:,:,:) = zero
701 16 : do iatom=1,Invar%natom_unitcell
702 11 : itypat = Invar%typat_unitcell(iatom)
703 49 : do ii=1,3
704 44 : Invar%zeff(ii,ii,iatom) = born_charge_tmp(itypat)
705 : end do
706 : end do
707 : ! GA: TODO eventually change the way this is reported
708 5 : write(Invar%stdout,'(1x,a20,20(1x,f15.10))') ljust('born_charge',20),(born_charge_tmp(jj),jj=1,Invar%ntypat)
709 5 : ABI_FREE(born_charge_tmp)
710 : end if
711 :
712 : ! dielec_constant
713 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'dielec_constant', tread, 'DPR')
714 44 : if (tread == 1) then
715 20 : do ii=1,3
716 20 : Invar%dielt(ii,ii) = dprarr(1)
717 : end do
718 5 : write(Invar%stdout,'(1x,a20,1x,f15.10)') ljust('dielec_constant',20),Invar%dielt(1,1)
719 : end if
720 :
721 : ! bzpath
722 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'bzpath', tread, 'INT')
723 44 : if (tread == 1) then
724 10 : Invar%bzpath = intarr(1)
725 10 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('bzpath',20),Invar%bzpath
726 :
727 10 : if (Invar%bzpath.lt.0) then
728 0 : ABI_MALLOC(Invar%qpt,(3,abs(Invar%bzpath))); Invar%qpt(:,:)=zero
729 :
730 0 : call intagm(dprarr, intarr, jdtset, marr, 1+3*abs(Invar%bzpath), string(1:lenstr), 'bzpath', tread, 'DPR')
731 0 : Invar%qpt(:,:) = reshape(dprarr(2:1+3*abs(Invar%bzpath)), (/3,abs(Invar%bzpath)/))
732 :
733 0 : write(Invar%stdout,'(a)') ' Q points as given in the input file:'
734 0 : do jj=1,abs(Invar%bzpath)
735 0 : write(Invar%stdout,'(22x,3(f15.10,1x))') Invar%qpt(:,jj)
736 : end do
737 :
738 10 : else if (Invar%bzpath.gt.0) then
739 20 : ABI_MALLOC(Invar%special_qpt,(Invar%bzpath))
740 10 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'special_qpt', tread, 'KEY', key_value=msg)
741 10 : if (tread == 1) then
742 10 : jj = 1
743 64 : do ii=1,Invar%bzpath
744 54 : Invar%special_qpt(ii) = " "
745 64 : ierr = next_token(msg,jj,Invar%special_qpt(ii))
746 : !jj = 2*ii-1
747 : !Invar%special_qpt(ii) = msg(jj:jj)
748 : end do
749 10 : write(Invar%stdout,'(a,1x,a2,10("-",a2))') ' Special q-points: ',Invar%special_qpt(:)
750 : else
751 : write(msg,*)&
752 0 : 'Variable bzpath > 0, but variable "special_qpt" was not found in input file.'
753 0 : ABI_ERROR(msg)
754 : end if
755 : end if
756 : end if
757 :
758 : ! order
759 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'order', tread, 'INT')
760 44 : if (tread == 1) then
761 14 : Invar%order = intarr(1)
762 14 : if (Invar%order.eq.3) then
763 6 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'order', tread, 'DPR')
764 6 : Invar%rcut3 = dprarr(2)
765 6 : write(Invar%stdout,'(1x,a20,1x,i4,1x,f15.10)') ljust('order',20),Invar%order,Invar%rcut3
766 6 : if (Invar%rcut3.gt.Invar%rcut) then
767 0 : ABI_ERROR('The cutoff radius of the third order cannot be greater than the second order one.')
768 : end if
769 8 : else if (Invar%order.eq.4) then
770 8 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'order', tread, 'DPR')
771 8 : Invar%rcut3 = dprarr(2)
772 8 : Invar%rcut4 = dprarr(3)
773 8 : write(Invar%stdout,'(1x,a20,1x,i4,2(1x,f15.10))') ljust('order',20),Invar%order,Invar%rcut3,Invar%rcut4
774 8 : if (Invar%rcut4.gt.Invar%rcut) then
775 0 : ABI_ERROR('The cutoff radius of the fourth order cannot be greater than the second order one.')
776 : end if
777 : else
778 0 : ABI_ERROR('Only the 3rd and 4th orders are allowed. Change your input file.')
779 : end if
780 : end if
781 :
782 : ! slice
783 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'slice', tread, 'INT')
784 44 : if (tread == 1) then
785 2 : Invar%slice = intarr(1)
786 2 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('slice',20),Invar%slice
787 : end if
788 :
789 : ! enunit
790 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'enunit', tread, 'INT')
791 44 : if (tread == 1) then
792 36 : Invar%enunit = intarr(1)
793 36 : if (Invar%enunit.eq.0) write(Invar%stdout,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in meV)'
794 36 : if (Invar%enunit.eq.1) write(Invar%stdout,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in cm-1)'
795 36 : if (Invar%enunit.eq.2) write(Invar%stdout,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in mHa)'
796 36 : if (Invar%enunit.eq.3) write(Invar%stdout,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in THz)'
797 : end if
798 :
799 : ! nproc
800 44 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'nproc', tread, 'INT')
801 44 : if (tread == 1) then
802 0 : Invar%nproc(1:2) = intarr(1:2)
803 0 : write(Invar%stdout,'(1x,a20,1x,i4,1x,i4)') ljust('nproc',20),Invar%nproc(1),Invar%nproc(2)
804 : end if
805 :
806 : ! readifc
807 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'readifc', tread, 'INT')
808 44 : if (tread == 1) then
809 3 : Invar%readifc = intarr(1)
810 3 : if (Invar%readifc.eq.1) then
811 1 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'readifc', tread, 'DPR')
812 1 : Invar%tolread = dprarr(1)
813 1 : write(Invar%stdout,'(1x,a20,1x,i4,1x,f15.10)') ljust('readifc',20),Invar%readifc,Invar%tolread
814 : else
815 2 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('readifc',20),Invar%readifc
816 : end if
817 : end if
818 :
819 : ! alloy
820 44 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'alloy', tread, 'INT')
821 44 : if (tread == 1) then
822 1 : Invar%alloy = intarr(1)
823 1 : Invar%ityp_alloy1 = intarr(2)
824 1 : Invar%ityp_alloy2 = intarr(3)
825 1 : write(Invar%stdout,'(1x,a20,1x,3(i4,1x))') ljust('alloy',20),Invar%alloy,Invar%ityp_alloy1,Invar%ityp_alloy2
826 : end if
827 :
828 : ! together
829 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'together', tread, 'INT')
830 44 : if (tread == 1) then
831 1 : Invar%together = intarr(1)
832 1 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('together',20),Invar%together
833 : end if
834 :
835 : ! bzlength
836 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'bzlength', tread, 'INT')
837 44 : if (tread == 1) then
838 2 : Invar%bzlength = intarr(1)
839 2 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('bzlength',20),Invar%bzlength
840 6 : ABI_MALLOC(Invar%lgth_segments,(Invar%bzlength))
841 2 : call intagm(dprarr, intarr, jdtset, marr, 1+Invar%bzlength, string(1:lenstr), 'bzlength', tread, 'DPR')
842 14 : Invar%lgth_segments = dprarr(2:1+Invar%bzlength)
843 2 : write(Invar%stdout,'(a,1x,i3,10("-",i3))') ' Length of BZ : ',Invar%lgth_segments(:)
844 : end if
845 :
846 : ! ngqpt1
847 44 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'ngqpt1', tread, 'INT')
848 44 : if (tread == 1) then
849 44 : Invar%ngqpt1 = intarr(1:3)
850 11 : write(Invar%stdout,'(1x,a20,1x,3(i4,1x))') ljust('ngqpt1',20),Invar%ngqpt1(:)
851 : end if
852 :
853 : ! ngqpt2
854 44 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'ngqpt2', tread, 'INT')
855 44 : if (tread == 1) then
856 48 : Invar%ngqpt2 = intarr(1:3)
857 12 : write(Invar%stdout,'(1x,a20,1x,3(i4,1x))') ljust('ngqpt2',20),Invar%ngqpt2(:)
858 : end if
859 :
860 : ! tolmotifinboxmatch
861 : !call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'tolmotifinboxmatch', tread, 'DPR')
862 : !if (tread == 1) then
863 : ! Invar%tolmotif = dprarr(1)
864 : ! Invar%tolinbox = dprarr(2)
865 : ! Invar%tolmatch = dprarr(3)
866 : ! write(Invar%stdout,'(1x,a20,f10.5)') ljust('tolmotif',20),Invar%tolmotif
867 : ! write(Invar%stdout,'(1x,a20,f10.5)') ljust('tolinbox',20),Invar%tolinbox
868 : ! write(Invar%stdout,'(1x,a20,f10.5)') ljust('tolmatch',20),Invar%tolmatch
869 : !end if
870 :
871 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'tolmotif', tread, 'DPR')
872 44 : if (tread == 1) then
873 4 : Invar%tolmotif = dprarr(1)
874 4 : write(Invar%stdout,'(1x,a20,f10.5)') ljust('tolmotif',20),Invar%tolmotif
875 : end if
876 :
877 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'tolinbox', tread, 'DPR')
878 44 : if (tread == 1) then
879 4 : Invar%tolinbox = dprarr(1)
880 4 : write(Invar%stdout,'(1x,a20,f10.5)') ljust('tolinbox',20),Invar%tolinbox
881 : end if
882 :
883 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'tolmatch', tread, 'DPR')
884 44 : if (tread == 1) then
885 4 : Invar%tolmatch = dprarr(1)
886 4 : write(Invar%stdout,'(1x,a20,f10.5)') ljust('tolmatch',20),Invar%tolmatch
887 : end if
888 :
889 : ! use_weights
890 44 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'use_weights', tread, 'INT')
891 44 : if (tread == 1) then
892 1 : Invar%use_weights = intarr(1)
893 1 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('use_weights',20),Invar%use_weights
894 : end if
895 :
896 44 : ABI_FREE(intarr)
897 44 : ABI_FREE(dprarr)
898 :
899 : ! =========================================================================== !
900 : ! Output other information and perform some checks
901 :
902 44 : if (Invar%use_ideal_positions.eq.0) then
903 10 : write(Invar%stdout,'(a)') ' USE AVERAGE POSITIONS TO COMPUTE SPECTRUM'
904 34 : else if (Invar%use_ideal_positions.eq.1) then
905 34 : write(Invar%stdout,'(a)') ' USE IDEAL POSITIONS TO COMPUTE SPECTRUM'
906 : else
907 0 : write(Invar%stdout,'(a)') ' STOP: THIS VALUE IS NOT ALLOWED FOR use_ideal_positions'
908 : end if
909 44 : if (Invar%loto) write(Invar%stdout,'(a)') ' USE NON-ANALYTICAL CORRECTIONS (LO-TO)'
910 :
911 : ! Allowed values
912 44 : if ((Invar%together.ne.1).and.(Invar%together.ne.0)) then
913 0 : ABI_ERROR('STOP: The value of input variable TOGETHER is not allowed')
914 : end if
915 44 : if ((Invar%alloy.ne.1).and.(Invar%alloy.ne.0)) then
916 0 : ABI_ERROR('STOP: The value of input variable ALLOY is not allowed')
917 : end if
918 44 : if (Invar%alloy.ge.1) then
919 : if ((Invar%ityp_alloy1.lt.1).or.(Invar%ityp_alloy2.lt.1).or.&
920 1 : & (Invar%ityp_alloy1.gt.Invar%natom_unitcell).or.(Invar%ityp_alloy2.gt.Invar%natom_unitcell)) then
921 0 : ABI_ERROR('STOP: The value of input variables IALLOY are not allowed')
922 : end if
923 : end if
924 :
925 : ! Incompatible variables :
926 44 : if ((Invar%readifc.eq.1).and.(Invar%together.eq.1).and.(Invar%order.gt.2)) then
927 0 : ABI_ERROR('STOP: readifc=1, together=1 and order=3 or 4 are incompatible')
928 : end if
929 :
930 : ! =========================================================================== !
931 : ! Treat virtual crystal approximation (VCA) when alloy=1.
932 : ! Redefine all the data depending on (n)typat(_unitcell) and natom_unitcell
933 :
934 : ! The mixing of zeff will only work for the isotropic case,
935 : ! that is, if it was specified with input variables born_charge
936 :
937 44 : if (Invar%alloy.eq.1) then
938 1 : sum_alloy1=0
939 1 : sum_alloy2=0
940 129 : do iatom=1,Invar%natom
941 128 : if (Invar%typat(iatom).eq.Invar%ityp_alloy1) then
942 96 : sum_alloy1=sum_alloy1+1
943 : end if
944 129 : if (Invar%typat(iatom).eq.Invar%ityp_alloy2) then
945 32 : sum_alloy2=sum_alloy2+1
946 : end if
947 : end do
948 : amu_average =(Invar%amu(Invar%ityp_alloy1)*sum_alloy1+&
949 1 : & Invar%amu(Invar%ityp_alloy2)*sum_alloy2)/(sum_alloy1+sum_alloy2)
950 1 : if (Invar%loto) then
951 0 : zeff_average(:,:) = zero
952 0 : do iatom=1,Invar%natom_unitcell
953 0 : if (Invar%typat_unitcell(iatom).eq.Invar%ityp_alloy1) then
954 0 : zeff_average(:,:) = zeff_average(:,:) + Invar%zeff(:,:,iatom) * sum_alloy1
955 : exit
956 : end if
957 : end do
958 0 : do iatom=1,Invar%natom_unitcell
959 0 : if (Invar%typat_unitcell(iatom).eq.Invar%ityp_alloy2) then
960 0 : zeff_average(:,:) = zeff_average(:,:) + Invar%zeff(:,:,iatom) * sum_alloy2
961 : exit
962 : end if
963 : end do
964 0 : zeff_average(:,:) = zeff_average(:,:) / (sum_alloy1 + sum_alloy2)
965 0 : do iatom=1,Invar%natom_unitcell
966 0 : if (Invar%typat_unitcell(iatom).eq.min(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
967 0 : Invar%zeff(:,:,iatom) = zeff_average(:,:)
968 : end if
969 : end do
970 : end if
971 1 : shift=0
972 3 : do iatom=1,Invar%natom_unitcell
973 3 : if (Invar%typat_unitcell(iatom).lt.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
974 1 : Invar%typat_unitcell (iatom-shift)=Invar%typat_unitcell (iatom)
975 4 : Invar%xred_unitcell(:,iatom-shift)=Invar%xred_unitcell(:,iatom)
976 1 : if (Invar%loto) then
977 0 : Invar%zeff(:,:,iatom-shift) = Invar%zeff(:,:,iatom)
978 : end if
979 1 : else if (Invar%typat_unitcell(iatom).eq.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
980 1 : shift=shift+1
981 : else if (Invar%typat_unitcell(iatom).gt.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
982 0 : Invar%typat_unitcell (iatom-shift)=Invar%typat_unitcell (iatom) - 1
983 0 : Invar%xred_unitcell(:,iatom-shift)=Invar%xred_unitcell(:,iatom)
984 0 : if (Invar%loto) then
985 0 : Invar%zeff(:,:,iatom-shift) = Invar%zeff(:,:,iatom)
986 : end if
987 : end if
988 : end do
989 1 : Invar%natom_unitcell=Invar%natom_unitcell-shift
990 129 : do iatom=1,Invar%natom
991 129 : if (Invar%typat(iatom).ge.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
992 32 : Invar%typat(iatom)=Invar%typat(iatom) - 1
993 : end if
994 : end do
995 3 : do itypat=1,Invar%ntypat
996 3 : if (itypat.eq.min(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
997 1 : Invar%amu (itypat)=amu_average
998 1 : else if (itypat.gt.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
999 0 : Invar%amu (itypat-1)=Invar%amu (itypat)
1000 : end if
1001 : end do
1002 1 : Invar%ntypat=Invar%ntypat-1
1003 3 : ABI_MALLOC(typat_unitcell_tmp,( Invar%natom_unitcell))
1004 3 : ABI_MALLOC(xred_unitcell_tmp ,(3,Invar%natom_unitcell))
1005 3 : ABI_MALLOC(amu_tmp ,( Invar%ntypat))
1006 1 : if (Invar%loto) then
1007 0 : ABI_MALLOC(zeff_tmp ,(3,3,Invar%natom_unitcell))
1008 : end if
1009 2 : typat_unitcell_tmp (:)=Invar%typat_unitcell(1:Invar%natom_unitcell)
1010 5 : xred_unitcell_tmp(:,:)=Invar%xred_unitcell(:,1:Invar%natom_unitcell)
1011 2 : amu_tmp (:)=Invar%amu(1:Invar%ntypat)
1012 1 : if (Invar%loto) then
1013 0 : zeff_tmp(:,:,:) = Invar%zeff(:,:,1:Invar%natom_unitcell)
1014 : end if
1015 3 : ABI_REMALLOC(Invar%typat_unitcell,( Invar%natom_unitcell))
1016 3 : ABI_REMALLOC(Invar%xred_unitcell ,(3,Invar%natom_unitcell))
1017 3 : ABI_REMALLOC(Invar%amu ,( Invar%ntypat))
1018 3 : ABI_REMALLOC(Invar%zeff ,(3,3,Invar%natom_unitcell))
1019 2 : Invar%typat_unitcell (:)=typat_unitcell_tmp (:)
1020 5 : Invar%xred_unitcell(:,:)=xred_unitcell_tmp(:,:)
1021 2 : Invar%amu (:)=amu_tmp (:)
1022 1 : if (Invar%loto) then
1023 0 : Invar%zeff(:,:,:) =zeff_tmp(:,:,:)
1024 : end if
1025 1 : if (allocated(Invar%znucl)) then
1026 3 : ABI_MALLOC(znucl_tmp,(Invar%ntypat))
1027 2 : znucl_tmp(:) = Invar%znucl(1:Invar%ntypat)
1028 3 : ABI_REMALLOC(Invar%znucl,(Invar%ntypat))
1029 2 : Invar%znucl(:) = znucl_tmp(:)
1030 1 : ABI_FREE(znucl_tmp)
1031 : end if
1032 1 : ABI_FREE(typat_unitcell_tmp)
1033 1 : ABI_FREE(xred_unitcell_tmp)
1034 1 : ABI_FREE(amu_tmp)
1035 1 : if (Invar%loto) then
1036 0 : ABI_FREE(zeff_tmp)
1037 : end if
1038 :
1039 1 : write(Invar%stdout,'(a)') ' ==================== Virtual Crystal Approximation =========================='
1040 1 : write(Invar%stdout,'(a)') ' ================ Several input variables are modified ======================='
1041 1 : write(Invar%stdout,'(a)') ' --> Beginning of the modifications'
1042 1 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('ntypat',20),Invar%ntypat
1043 1 : write(Invar%stdout,'(1x,a20,1x,i4)') ljust('natom_unitcell',20),Invar%natom_unitcell
1044 2 : write(Invar%stdout,'(1x,a20,20(1x,f15.10))') ljust('amu',20),(Invar%amu(jj),jj=1,Invar%ntypat)
1045 2 : write(Invar%stdout,'(1x,a20,20(1x,i4))') ljust('typat_unitcell',20),(Invar%typat_unitcell(jj),jj=1,Invar%natom_unitcell)
1046 1 : write(Invar%stdout,'(1x,a20)') ljust('xred_unitcell',20)
1047 2 : do ii=1,Invar%natom_unitcell
1048 5 : write(Invar%stdout,'(22x,3(f15.10,1x))') (Invar%xred_unitcell(jj,ii), jj=1,3)
1049 : end do
1050 1 : if (Invar%loto) then
1051 : ! GA: TODO eventually change the way this is reported
1052 : ! At the moment it is inconsistent with previous report.
1053 0 : write(Invar%stdout,'(1x,a20,20(1x,f15.10))') ljust('born_charge',20),(Invar%zeff(1,1,jj),jj=1,Invar%natom_unitcell)
1054 : end if
1055 1 : write(Invar%stdout,'(1x,a20)') ljust('typat',20)
1056 14 : do ii=1,Invar%natom,10
1057 14 : if (ii+9.lt.Invar%natom) then
1058 132 : write(Invar%stdout,'(22x,10(i4,1x))') (Invar%typat(ii+jj-1),jj=1,10)
1059 : else
1060 9 : write(Invar%stdout,'(22x,10(i4,1x))') (Invar%typat(jj),jj=ii,Invar%natom)
1061 : end if
1062 : end do
1063 1 : write(Invar%stdout,'(a)') ' --> End of the modifications'
1064 1 : write(Invar%stdout,'(a)') ' '
1065 : end if
1066 :
1067 : ! Compute Nstep as a function of the slice
1068 44 : Invar%nstep_tot=int(float(Invar%nstep_max-Invar%nstep_min)/float(Invar%slice)+1)
1069 44 : write(Invar%stdlog,*) 'nstep_tot=',Invar%nstep_tot
1070 :
1071 :
1072 88 : end subroutine tdep_read_input
1073 :
1074 : !====================================================================================================
1075 44 : subroutine tdep_destroy_invar(Invar)
1076 :
1077 : type(atdep_dataset_type), intent(inout) :: Invar
1078 :
1079 44 : ABI_FREE(Invar%amu)
1080 44 : ABI_FREE(Invar%typat)
1081 44 : ABI_FREE(Invar%xred_unitcell)
1082 44 : ABI_FREE(Invar%typat_unitcell)
1083 44 : if (Invar%bzpath.lt.0) then
1084 0 : ABI_FREE(Invar%qpt)
1085 44 : else if (Invar%bzpath.gt.0) then
1086 10 : ABI_FREE(Invar%special_qpt)
1087 : end if
1088 44 : if (Invar%bzlength.gt.0) then
1089 2 : ABI_FREE(Invar%lgth_segments)
1090 : end if
1091 44 : ABI_SFREE(Invar%zeff)
1092 44 : ABI_SFREE(Invar%znucl)
1093 :
1094 44 : end subroutine tdep_destroy_invar
1095 :
1096 : !====================================================================================================
1097 :
1098 44 : subroutine tdep_init_MPIdata(Invar,MPIdata)
1099 :
1100 : type(atdep_dataset_type), intent(inout) :: Invar
1101 : type(MPI_enreg_type), intent(out) :: MPIdata
1102 : integer :: ii,remain,ierr,iproc,istep
1103 44 : integer, allocatable :: tab_step(:)
1104 : character(len=500) :: message
1105 :
1106 : #if defined HAVE_MPI
1107 : integer :: dimcart,commcart_2d,me_cart_2d
1108 : logical :: reorder
1109 44 : integer,allocatable :: coords(:),sizecart(:)
1110 44 : logical,allocatable :: periode(:), keepdim(:)
1111 : #endif
1112 :
1113 : ! Check the number of processors
1114 44 : MPIdata%nproc_shell=Invar%nproc(1)
1115 44 : MPIdata%nproc_step =Invar%nproc(2)
1116 44 : MPIdata%nproc = xmpi_comm_size(xmpi_world)
1117 44 : if (MPIdata%nproc_step*MPIdata%nproc_shell.ne.MPIdata%nproc) then
1118 0 : ABI_WARNING('The parallelization is performed over steps')
1119 0 : MPIdata%nproc_step = xmpi_comm_size(xmpi_world)
1120 : end if
1121 :
1122 44 : MPIdata%master = 0
1123 44 : MPIdata%iam_master =.false.
1124 : ! Initialize the MPIdata datastructure for sequential calculation
1125 44 : if (MPIdata%nproc.eq.1) then
1126 44 : MPIdata%comm_shell = xmpi_comm_null
1127 44 : MPIdata%comm_step = xmpi_comm_null
1128 44 : MPIdata%comm_shellstep = xmpi_comm_null
1129 44 : MPIdata%me_shell = 0
1130 44 : MPIdata%me_step = 0
1131 44 : MPIdata%my_nstep = Invar%nstep_tot
1132 44 : MPIdata%iam_master = (MPIdata%me_step == MPIdata%master)
1133 : end if
1134 :
1135 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1136 : !!!!!!!!!!!!!!!!!!!!!! Parallel calculation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1137 : !!!!!!!!!!!!!! Definition of the processor grid !!!!!!!!!!!!!!!!!!!!!!!!!!!!
1138 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1139 : #if defined HAVE_MPI
1140 : !FB if (MPIdata%nproc.eq.1) return
1141 :
1142 : ! Create the global cartesian 2D-communicator
1143 44 : dimcart=2
1144 44 : ABI_MALLOC(sizecart,(dimcart))
1145 44 : ABI_MALLOC(periode,(dimcart))
1146 44 : sizecart(1)=MPIdata%nproc_shell ! MPIdata%nproc_shell
1147 44 : sizecart(2)=MPIdata%nproc_step ! MPIdata%nproc_step
1148 132 : periode(:)=.false.;reorder=.false.
1149 44 : call MPI_CART_CREATE(xmpi_world,dimcart,sizecart,periode,reorder,commcart_2d,ierr)
1150 44 : ABI_FREE(periode)
1151 44 : ABI_FREE(sizecart)
1152 :
1153 : ! Find the index and coordinates of the current processor
1154 44 : call MPI_COMM_RANK(commcart_2d,me_cart_2d,ierr)
1155 132 : ABI_MALLOC(coords,(dimcart))
1156 44 : call MPI_CART_COORDS(commcart_2d,me_cart_2d,dimcart,coords,ierr)
1157 44 : MPIdata%me_shell=coords(1)
1158 44 : MPIdata%me_step =coords(2)
1159 44 : ABI_FREE(coords)
1160 44 : if ((MPIdata%me_shell == MPIdata%master).and.(MPIdata%me_step == MPIdata%master)) then
1161 44 : MPIdata%iam_master = .true.
1162 : end if
1163 :
1164 132 : ABI_MALLOC(keepdim,(dimcart))
1165 : ! Create the communicator for shell distribution
1166 44 : keepdim(1)=.true.
1167 44 : keepdim(2)=.false.
1168 44 : call MPI_CART_SUB(commcart_2d,keepdim,MPIdata%comm_shell,ierr)
1169 : ! Create the communicator for step distribution
1170 44 : keepdim(1)=.false.
1171 44 : keepdim(2)=.true.
1172 44 : call MPI_CART_SUB(commcart_2d,keepdim,MPIdata%comm_step,ierr)
1173 : ! Create the communicator for shellstep distribution
1174 44 : keepdim(1)=.true.
1175 44 : keepdim(2)=.true.
1176 44 : call MPI_CART_SUB(commcart_2d,keepdim,MPIdata%comm_shellstep,ierr)
1177 44 : ABI_FREE(keepdim)
1178 44 : call xmpi_comm_free(commcart_2d)
1179 :
1180 : ! Write some data
1181 44 : write(message,'(a21,2(1x,i4))') '-Number of processors',MPIdata%nproc_shell,MPIdata%nproc_step
1182 44 : call wrtout(Invar%stdout,message,'COLL')
1183 : !FB write(message,'(a,2i5)') 'me_shell and me_step : ',MPIdata%me_shell,MPIdata%me_step
1184 : !FB call wrtout(Invar%stdout,message,'COLL')
1185 : #endif
1186 :
1187 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1188 : !!!!!!!!!!!!!! Distribution over STEP processors !!!!!!!!!!!!!!!!!!!!!!!!!!!
1189 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1190 44 : MPIdata%my_nstep =int(Invar%nstep_tot/MPIdata%nproc_step)
1191 44 : remain=Invar%nstep_tot-MPIdata%nproc_step*MPIdata%my_nstep
1192 44 : do ii=1,remain
1193 44 : if ((ii-1).eq.MPIdata%me_step) MPIdata%my_nstep=MPIdata%my_nstep+1
1194 : end do
1195 176 : ABI_MALLOC(MPIdata%nstep_all,(MPIdata%nproc_step)); MPIdata%nstep_all(:)=zero
1196 44 : call xmpi_allgather(MPIdata%my_nstep,MPIdata%nstep_all,MPIdata%comm_step,ierr)
1197 44 : write(Invar%stdout,'(a)') ' '
1198 44 : write(Invar%stdout,'(a,1x,i4)') ' All quantities are computed from nstep_min=',Invar%nstep_min
1199 44 : write(Invar%stdout,'(a,1x,i4)') ' to nstep_max=',Invar%nstep_max
1200 44 : if (Invar%slice.ne.1) then
1201 2 : write(Invar%stdout,'(a,1x,i4)') ' by using a slice=',Invar%slice
1202 : end if
1203 44 : write(Invar%stdout,'(a,1x,i4)') ' So, the real number of time steps is nstep=',Invar%nstep_tot
1204 44 : if (MPIdata%nproc_step.gt.1) then
1205 0 : write(Invar%stdout,'(a,1000(1x,i5))') '-Distribution of number of steps wrt the number of processors=',MPIdata%nstep_all(:)
1206 : end if
1207 :
1208 220 : ABI_MALLOC(MPIdata%nstep_acc,(MPIdata%nproc_step+1)); MPIdata%nstep_acc(:)=zero
1209 44 : MPIdata%nstep_acc(1)=0
1210 88 : do ii=2,MPIdata%nproc_step+1
1211 88 : MPIdata%nstep_acc(ii)=MPIdata%nstep_acc(ii-1)+MPIdata%nstep_all(ii-1)
1212 : end do
1213 44 : if (MPIdata%nstep_acc(MPIdata%nproc_step+1).ne.Invar%nstep_tot) then
1214 0 : write(Invar%stdlog,*) 'STOP : pb in nstep_acc'
1215 0 : stop
1216 : end if
1217 :
1218 1084 : ABI_MALLOC(tab_step,(Invar%nstep_tot)); tab_step(:)=zero
1219 1040 : ABI_MALLOC(MPIdata%my_step ,(Invar%nstep_tot)); MPIdata%my_step (:)=.false.
1220 88 : do iproc=1,MPIdata%nproc_step
1221 1040 : do istep=1,Invar%nstep_tot
1222 996 : if ((istep.gt.MPIdata%nstep_acc(iproc)).and.(istep.le.MPIdata%nstep_acc(iproc+1))) then
1223 952 : tab_step(istep)=iproc-1
1224 : end if
1225 : end do
1226 : end do
1227 996 : do istep=1,Invar%nstep_tot
1228 996 : MPIdata%my_step(istep) = (tab_step(istep) == MPIdata%me_step)
1229 : end do
1230 :
1231 176 : ABI_MALLOC(MPIdata%shft_step,(MPIdata%nproc_step)); MPIdata%shft_step(:)=zero
1232 44 : MPIdata%shft_step(1)=0
1233 44 : do ii=2,MPIdata%nproc_step
1234 44 : MPIdata%shft_step(ii)=MPIdata%shft_step(ii-1)+MPIdata%nstep_all(ii-1)
1235 : end do
1236 44 : ABI_FREE(tab_step)
1237 :
1238 44 : Invar%my_nstep=MPIdata%my_nstep
1239 :
1240 44 : end subroutine tdep_init_MPIdata
1241 :
1242 : !====================================================================================================
1243 :
1244 44 : subroutine tdep_destroy_mpidata(MPIdata)
1245 :
1246 : type(MPI_enreg_type), intent(inout) :: MPIdata
1247 :
1248 44 : ABI_FREE(MPIdata%shft_step)
1249 44 : ABI_FREE(MPIdata%nstep_acc)
1250 44 : ABI_FREE(MPIdata%nstep_all)
1251 44 : ABI_FREE(MPIdata%my_step)
1252 :
1253 44 : end subroutine tdep_destroy_mpidata
1254 :
1255 : !====================================================================================================
1256 :
1257 0 : end module m_tdep_dataset
|
