Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_dataset
2 : : !! NAME
3 : : !! m_tdep_dataset
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the data type holding the input variables of atdep
7 : : !! (atdep_dataset_type) and the routines reading and echoing them, as well
8 : : !! as the data type describing the MPI parallelism (MPI_enreg_type) and
9 : : !! the distribution of the MD steps over the processors.
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_tdep_dataset
26 : :
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_xmpi
31 : : use m_abihist
32 : : use m_parser
33 : : use netcdf
34 : : use m_nctk
35 : : use m_fstrings, only : inupper,ljust,next_token
36 : : use m_io_tools, only : file_exists
37 : : use m_ddb, only : ddb_type
38 : : use m_ddb_hdr, only : ddb_hdr_type
39 : : use m_crystal, only : crystal_t
40 : :
41 : : implicit none
42 : : !!***
43 : :
44 : : !----------------------------------------------------------------------
45 : :
46 : : !!****t* m_tdep_dataset/atdep_dataset_type
47 : : !! NAME
48 : : !! atdep_dataset_type
49 : : !!
50 : : !! FUNCTION
51 : : !! Holds the input variables of atdep (read from the input file, the
52 : : !! HIST.nc file and the DDB file), as well as a few derived quantities.
53 : : !!
54 : : !! SOURCE
55 : :
56 : : type atdep_dataset_type
57 : :
58 : : integer :: natom
59 : : ! Number of atoms in the supercell
60 : :
61 : : integer :: natom_unitcell
62 : : ! Number of atoms in the unit cell
63 : :
64 : : integer :: nshell_max
65 : : ! Maximum number of shells allowed (hard-coded)
66 : :
67 : : integer :: nstep_max
68 : : ! Last MD step used
69 : :
70 : : integer :: nstep_min
71 : : ! First MD step used
72 : :
73 : : integer :: nstep_tot
74 : : ! Total number of MD steps used: (nstep_max-nstep_min)/slice+1
75 : :
76 : : integer :: my_nstep
77 : : ! Number of MD steps held by the present processor (see tdep_init_MPIdata)
78 : :
79 : : integer :: ntypat
80 : : ! Number of atom types
81 : :
82 : : integer :: use_ideal_positions
83 : : ! Compute the displacements with respect to the average positions (0)
84 : : ! or to the ideal positions (1)
85 : :
86 : : integer :: ab_out
87 : : ! Unit number of the main output file
88 : :
89 : : integer :: std_out
90 : : ! Unit number of the standard output
91 : :
92 : : integer :: bzpath
93 : : ! Definition of the q-point path for the phonon dispersion:
94 : : ! > 0: number of special q-points given by their names (special_qpt)
95 : : ! < 0: number of q-points given by their reduced coordinates (qpt)
96 : :
97 : : integer :: ndivsm
98 : : ! Number of divisions of the smallest segment of the q-point path
99 : :
100 : : integer :: order
101 : : ! Maximum IFC order (2, 3, or 4)
102 : :
103 : : integer :: slice
104 : : ! Only one MD step every slice steps is used
105 : :
106 : : integer :: enunit
107 : : ! Unit of the phonon frequencies: 0=meV, 1=cm-1, 2=mHa, 3=THz
108 : :
109 : : integer :: readifc
110 : : ! Read the IFC from file: 0=no, 1=from _ifc_in.dat, 2=from _ifc_out.dat
111 : :
112 : : integer :: together
113 : : ! Compute the IFC of all orders simultaneously (1) or successively (0)
114 : :
115 : : integer :: alloy
116 : : ! Use the virtual crystal approximation (1) to merge the atom types
117 : : ! ityp_alloy1 and ityp_alloy2
118 : :
119 : : integer :: ityp_alloy1
120 : : ! First atom type merged when alloy=1
121 : :
122 : : integer :: ityp_alloy2
123 : : ! Second atom type merged when alloy=1
124 : :
125 : : integer :: nproc(2)
126 : : ! Number of processors for the distribution over shells (1) and over MD steps (2)
127 : :
128 : : integer :: bzlength
129 : : ! Number of segments of the q-point path whose lengths are given
130 : : ! explicitly (lgth_segments)
131 : :
132 : : integer :: ngqpt1(3)
133 : : ! q-point grid used to build the IFC and DDB objects
134 : :
135 : : integer :: ngqpt2(3)
136 : : ! q-point grid used to compute the phonon density of states
137 : :
138 : : integer :: bravais(11)
139 : : ! Bravais lattice (bravais(1)) and centering (bravais(2)), see symlatt
140 : :
141 : : integer :: use_weights
142 : : ! Read the weights of the MD configurations from the _weights.dat file (1),
143 : : ! or use the same weight for each of them (0)
144 : :
145 : : integer :: prtphdos
146 : : ! Print the phonon density of states (1)
147 : :
148 : : integer :: prtvol ! Verbosity level (see the ABINIT input variable of the same name)
149 : : integer :: symbec ! Symmetrize born effective charges
150 : :
151 : : integer, allocatable :: typat_unitcell(:)
152 : : ! typat_unitcell(natom_unitcell)
153 : : ! Atom type in the unitcell
154 : :
155 : : integer, allocatable :: typat(:)
156 : : ! typat(natom)
157 : : ! Atom type in the supercell
158 : :
159 : : integer, allocatable :: lgth_segments(:)
160 : : ! lgth_segments(bzlength)
161 : : ! Number of q-points along each segment of the q-point path
162 : :
163 : : logical :: dipdip
164 : : ! Whether the dipole-dipole (long-range) interaction is subtracted from
165 : : ! the dynamical matrix, a la [[anaddb:dipdip]] (same physical meaning,
166 : : ! but atdep only ever supports the on/off choice: unlike anaddb, there is
167 : : ! no dipdip=-1 fast-algorithm variant, and it is auto-enabled whenever
168 : : ! usable Born effective charges/dielectric tensor are found, see
169 : : ! tdep_read_input). Explicitly setting the input variable dipdip=0
170 : : ! forces it off even when such data is available.
171 : :
172 : : logical :: netcdf
173 : : ! Whether the MD trajectory is read from the HIST.nc netcdf file
174 : :
175 : : double precision :: angle_alpha
176 : : ! Angle (in degrees) defining the lattice, for some Bravais lattices
177 : :
178 : : double precision :: dosdeltae
179 : : ! Energy step of the phonon density of states (Ha)
180 : :
181 : : double precision :: rcut
182 : : ! Cutoff radius of the second-order IFC (bohr)
183 : :
184 : : double precision :: rcut3
185 : : ! Cutoff radius of the third-order IFC (bohr)
186 : :
187 : : double precision :: rcut4
188 : : ! Cutoff radius of the fourth-order IFC (bohr)
189 : :
190 : : double precision :: temperature
191 : : ! Temperature of the MD simulation (K)
192 : :
193 : : double precision :: tolread
194 : : ! Tolerance given together with readifc=1
195 : :
196 : : double precision :: rprimd_md(3,3)
197 : : ! Dimensioned primitive vectors
198 : :
199 : : double precision :: multiplicity(3,3)
200 : : ! Supercell scaling matrix (integers)
201 : :
202 : : double precision :: dielt(3,3)
203 : : ! Dielectric tensor
204 : :
205 : : double precision, allocatable :: amu(:)
206 : : ! amu(ntypat)
207 : : ! Mass of the atoms (atomic mass unit)
208 : :
209 : : double precision, allocatable :: zeff(:,:,:)
210 : : ! zeff(3,3,natom_unitcell)
211 : : ! Born effective charge tensor
212 : :
213 : : double precision, allocatable :: znucl(:)
214 : : ! znucl(ntypat)
215 : : ! Nuclear charge of each atom type
216 : :
217 : : double precision, allocatable :: qpt(:,:)
218 : : ! qpt(3,abs(bzpath))
219 : : ! Reduced coordinates of the q-points defining the path (bzpath < 0)
220 : :
221 : : double precision, allocatable :: xred_unitcell(:,:)
222 : : ! xred_unitcell(3, natom_unitcell)
223 : : ! Reduced equilibrium position of the atoms in the unitcell, anchored so
224 : : ! that atom 1 sits at the origin (0,0,0) and wrapped into [0,1).
225 : :
226 : : double precision, allocatable :: xred_unitcell_original(:,:)
227 : : ! xred_unitcell_original(3, natom_unitcell)
228 : : ! Reduced equilibrium position of the atoms in the unitcell,
229 : : ! as given in the input file, wrapped into (-0.5,0.5].
230 : :
231 : : character (len=2), allocatable :: special_qpt(:)
232 : : ! special_qpt(bzpath)
233 : : ! Names of the special q-points defining the path (bzpath > 0)
234 : :
235 : : character (len=fnlen) :: output_prefix
236 : : ! Root name of the output files
237 : :
238 : : character (len=fnlen) :: input_prefix
239 : : ! Root name of the input files (HIST.nc, DDB, xred.dat, ...)
240 : :
241 : : character (len=fnlen) :: output_file
242 : : ! Name of the main output file
243 : :
244 : : contains
245 : :
246 : : procedure :: print => tdep_write_input
247 : : procedure :: free => tdep_destroy_invar
248 : :
249 : : end type atdep_dataset_type
250 : : !!***
251 : :
252 : : !----------------------------------------------------------------------
253 : :
254 : : !!****t* m_tdep_dataset/MPI_enreg_type
255 : : !! NAME
256 : : !! MPI_enreg_type
257 : : !!
258 : : !! FUNCTION
259 : : !! Holds the information on the MPI parallelism: the 2D grid of processors
260 : : !! (shells x MD steps), the corresponding communicators, and the
261 : : !! distribution of the MD steps over the processors.
262 : : !!
263 : : !! SOURCE
264 : :
265 : : type MPI_enreg_type
266 : :
267 : : integer :: comm_shell
268 : : ! Communicator for the distribution over shells
269 : :
270 : : integer :: comm_step
271 : : ! Communicator for the distribution over MD steps
272 : :
273 : : integer :: comm_shellstep
274 : : ! Communicator over the whole 2D grid of processors (shells x steps)
275 : :
276 : : integer :: nproc
277 : : ! Total number of processors
278 : :
279 : : integer :: nproc_shell
280 : : ! Number of processors for the distribution over shells
281 : :
282 : : integer :: nproc_step
283 : : ! Number of processors for the distribution over MD steps
284 : :
285 : : integer :: master
286 : : ! Rank of the master processor
287 : :
288 : : integer :: me_shell
289 : : ! Rank of the present processor in comm_shell
290 : :
291 : : integer :: me_step
292 : : ! Rank of the present processor in comm_step
293 : :
294 : : integer, allocatable :: my_nshell(:)
295 : : ! Currently unused
296 : :
297 : : integer :: my_nstep
298 : : ! Number of MD steps held by the present processor
299 : :
300 : : logical :: iam_master
301 : : ! Whether the present processor is the master
302 : :
303 : : integer, allocatable :: nstep_acc(:)
304 : : ! nstep_acc(nproc_step+1)
305 : : ! Accumulated number of MD steps held by the previous processors
306 : :
307 : : integer, allocatable :: nstep_all(:)
308 : : ! nstep_all(nproc_step)
309 : : ! Number of MD steps held by each processor
310 : :
311 : : integer, allocatable :: shft_step(:)
312 : : ! shft_step(nproc_step)
313 : : ! Offset of the MD steps held by each processor (used in xmpi_gatherv)
314 : :
315 : : logical, allocatable :: my_shell(:)
316 : : ! Currently unused
317 : :
318 : : logical, allocatable :: my_step(:)
319 : : ! my_step(nstep_tot)
320 : : ! Whether each MD step is held by the present processor
321 : :
322 : : contains
323 : :
324 : : procedure :: init => tdep_init_MPIdata
325 : : procedure :: free => tdep_destroy_mpidata
326 : :
327 : : end type MPI_enreg_type
328 : : !!***
329 : :
330 : : public :: tdep_read_input
331 : : public :: tdep_write_input
332 : : public :: tdep_destroy_invar
333 : : public :: tdep_init_MPIdata
334 : : public :: tdep_destroy_mpidata
335 : :
336 : : contains
337 : :
338 : : !====================================================================================================
339 : :
340 : : !!****f* ABINIT/m_tdep_dataset/tdep_read_input
341 : : !! NAME
342 : : !! tdep_read_input
343 : : !!
344 : : !! FUNCTION
345 : : !! Read the input variables of atdep and store them in the Invar object.
346 : : !! The names of the input file, of the input/output prefixes and of the
347 : : !! main output file are either given by input_path (and the input
348 : : !! variables indata_prefix, outdata_prefix and output_file), or read from
349 : : !! the standard input (deprecated "files" file). The main output file is
350 : : !! then opened.
351 : : !! If the HIST.nc netcdf file is available, the MD trajectory is read
352 : : !! into Hist, together with natom, ntypat, typat, znucl and amu.
353 : : !! If the DDB file is available, the Born effective charges and the
354 : : !! dielectric tensor are read from it (and dipdip is enabled).
355 : : !! Then, the mandatory and optional input variables are read and checked,
356 : : !! the virtual crystal approximation is applied (if alloy=1), and the
357 : : !! unitcell is anchored so that atom 1 sits at the origin. Finally,
358 : : !! the total number of MD steps (nstep_tot) is computed.
359 : : !!
360 : : !! INPUTS
361 : : !! input_path = Path of the input file. If empty, the names of the files
362 : : !! are read from the standard input.
363 : : !!
364 : : !! OUTPUT
365 : : !! Hist = abihist object containing the MD trajectory
366 : : !! (only if the HIST.nc file is available).
367 : : !! Invar = Input object containing the input variables.
368 : : !!
369 : : !! SOURCE
370 : :
371 [ - + - + : 46 : subroutine tdep_read_input(input_path,Hist,Invar)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
372 : :
373 : : ! Arguments-------------------------------
374 : : character(len=*), intent(in):: input_path
375 : : type(atdep_dataset_type),intent(out) :: Invar
376 : : type(abihist), intent(out) :: Hist
377 : :
378 : : ! Local variables-------------------------
379 : : ! scalars
380 : : integer :: ncid, ncerr, me,ierr,master
381 : : integer :: nimage, mdtime, natom_id,nimage_id,time_id,xyz_id,six_id
382 : : integer :: ntypat_id,iatcell
383 : : integer :: ii,jj,shift,iatom,itypat,sum_alloy1,sum_alloy2
384 : : integer:: lenstr, marr, jdtset, tread
385 : : logical :: has_nimage
386 : : double precision :: dtion,amu_average,born_average
387 : : character(len=500) :: msg
388 : : character(len=fnlen) :: input_filename,hist_filename,ddb_filename
389 : : character(len=strlen):: string, raw_string
390 [ + + + + : 2392 : type(crystal_t):: crystal
+ + + + +
+ + + + +
+ + + + ]
391 : 46 : type(ddb_type):: ddb
392 [ + + + + : 2392 : type(ddb_hdr_type):: ddb_hdr
+ + + + +
+ + + + +
+ + + + ]
393 : : ! arrays
394 : 46 : integer, allocatable:: intarr(:)
395 : 46 : integer, allocatable :: typat_unitcell_tmp(:)
396 : : real(dp) :: zeff_average(3,3)
397 : 46 : real(dp), allocatable:: zeff_tmp(:,:,:)
398 : 92 : real(dp), allocatable :: xred_unitcell_tmp(:,:),amu_tmp(:),born_charge_tmp(:),znucl_tmp(:)
399 [ - + ]: 46 : real(dp), allocatable:: dprarr(:)
400 : :
401 : : ! *********************************************************************
402 : :
403 : : ! Define output files
404 : 46 : Invar%ab_out=ab_out
405 : 46 : Invar%std_out=std_out
406 : :
407 : : ! Define default values
408 : 46 : Invar%nshell_max=500
409 : 46 : Invar%angle_alpha=90.d0
410 : 46 : Invar%bzpath=0
411 : 46 : Invar%ndivsm=20
412 : 46 : Invar%order=2
413 : 46 : Invar%slice=1
414 : 46 : Invar%enunit=0
415 : 46 : Invar%together=1
416 : 46 : Invar%alloy=0
417 : 46 : Invar%ityp_alloy1=0
418 : 46 : Invar%ityp_alloy2=0
419 [ + + ]: 138 : Invar%nproc(:)=1
420 : 46 : Invar%bzlength=0
421 : 46 : Invar%readifc=0
422 : 46 : Invar%tolread=1.d-8
423 : 46 : Invar%dosdeltae=0.2_dp/Ha_cmm1
424 : 46 : Invar%prtvol=0
425 : 46 : Invar%dipdip=.false.
426 : 46 : Invar%netcdf=.false.
427 : 46 : Invar%use_ideal_positions=0
428 : 46 : Invar%use_weights=0
429 : 46 : Invar%prtphdos=1
430 : 46 : Invar%symbec=0
431 : : ! In order to have an accuracy better than 1meV
432 [ + + ]: 184 : Invar%ngqpt1(:)=8
433 [ + + ]: 184 : Invar%ngqpt2(:)=32
434 [ + + + + ]: 598 : Invar%dielt(:,:)=zero; Invar%dielt(1,1)=one; Invar%dielt(2,2)=one; Invar%dielt(3,3)=one
435 : :
436 : 46 : master = 0
437 : 46 : me = xmpi_comm_rank(xmpi_world)
438 [ + - ]: 46 : if (me==master) then
439 : :
440 [ - + ]: 46 : if (len_trim(input_path) == 0) then
441 : :
442 : 0 : write(std_out, "(2a)")" DeprecationWarning: ",ch10
443 : 0 : write(std_out, "(a)") " The files file has been deprecated in Abinit10 and will be removed in Abinit11."
444 : 0 : write(std_out, "(2a)")" Use the syntax `atdep t01.abi` where t01.abi is an atdep input,",ch10
445 : 0 : write(std_out, "(2a)")" and use input variables output_file, indata_prefix, outdata_prefix.",ch10
446 : :
447 [ # # ]: 0 : write(std_out,'(a)',err=10) ' Give name for input file '
448 [ # # ]: 0 : read(*, '(a)',err=10) input_filename
449 [ # # ]: 0 : if ( input_filename == "" ) input_filename='input.in'
450 [ # # ]: 0 : write(std_out, '(a)',err=10) '.'//trim(input_filename)
451 : : 10 continue
452 : : ! Check if a NetCDF file is available
453 [ # # ]: 0 : write(std_out,'(a)',err=11) ' Give root name for generic input files (NetCDF or ASCII)'
454 [ # # ]: 0 : read(*, '(a)',err=11) Invar%input_prefix
455 [ # # ]: 0 : write(std_out, '(a)',err=11) '.'//trim(Invar%input_prefix)
456 : : 11 continue
457 [ # # ]: 0 : write(std_out,'(a)', err=12)' Give root name for generic output files:'
458 [ # # ]: 0 : read (*, '(a)', err=12) Invar%output_prefix
459 [ # # ]: 0 : if ( Invar%output_prefix == "" ) Invar%output_prefix = 'atdep'
460 [ # # ]: 0 : write (std_out, '(a)', err=12 ) '.'//trim(Invar%output_prefix)
461 : : 12 continue
462 [ # # ]: 0 : Invar%output_file = trim(Invar%output_prefix)//'.abo'
463 : : else
464 [ - + ]: 46 : input_filename = input_path
465 : :
466 : : ! Read input
467 : 46 : string = repeat(" ", strlen)
468 : 46 : raw_string = repeat(" ", strlen)
469 : 46 : call instrng(input_filename, lenstr, 1, strlen, string, raw_string)
470 : : ! To make case-insensitive, map characters to upper case.
471 : 46 : call inupper(string(1:lenstr))
472 : :
473 : 46 : marr = 3
474 [ - + ]: 46 : ABI_MALLOC(intarr, (marr))
475 [ - + ]: 46 : ABI_MALLOC(dprarr, (marr))
476 : 46 : jdtset = 0
477 : :
478 : : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), &
479 : 46 : "indata_prefix", tread, 'KEY', key_value=Invar%input_prefix)
480 [ - + ]: 46 : if (tread == 0) then
481 : 0 : Invar%input_prefix = ''
482 : : end if
483 : 46 : write(std_out, "(2a)")"Root name for input files: ", trim(Invar%input_prefix)
484 : :
485 : : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), &
486 : 46 : "outdata_prefix", tread, 'KEY', key_value=Invar%output_prefix)
487 [ - + ]: 46 : if (tread == 0) then
488 : 0 : Invar%output_prefix = 'atdep'
489 : : end if
490 : 46 : write(std_out, "(2a)")"Root name for output files: ", trim(Invar%output_prefix)
491 : :
492 : : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), &
493 : 46 : "output_file", tread, 'KEY', key_value=Invar%output_file)
494 [ - + ]: 46 : if (tread == 0) then
495 [ # # ]: 0 : Invar%output_file = trim(Invar%output_prefix) // '.abo'
496 : : end if
497 : 46 : write(std_out, "(2a)")"Main output file: ", trim(Invar%output_file)
498 : :
499 [ - + ]: 46 : ABI_FREE(intarr)
500 [ - + ]: 46 : ABI_FREE(dprarr)
501 : :
502 : : end if
503 : :
504 : : ! Setup other filenames
505 [ - + ]: 46 : if ( Invar%input_prefix == "" ) then
506 : 0 : hist_filename='HIST.nc'
507 : 0 : ddb_filename='DDB'
508 : : else
509 [ + - ]: 46 : hist_filename=trim(Invar%input_prefix)//'_HIST.nc'
510 [ + - ]: 46 : ddb_filename=trim(Invar%input_prefix)//'_DDB'
511 : : end if
512 : :
513 : 46 : open(unit=ab_out,file=trim(Invar%output_file))
514 : :
515 : : end if !me
516 : :
517 : 46 : call xmpi_bcast(input_filename,master,xmpi_world,ierr)
518 : 46 : call xmpi_bcast(hist_filename,master,xmpi_world,ierr)
519 : 46 : call xmpi_bcast(ddb_filename,master,xmpi_world,ierr)
520 : 46 : call xmpi_bcast(Invar%output_prefix,master,xmpi_world,ierr)
521 : 46 : call xmpi_bcast(Invar%input_prefix,master,xmpi_world,ierr)
522 : :
523 : : #if defined HAVE_NETCDF
524 : : !Open netCDF file
525 : 46 : ncerr=nf90_open(path=trim(hist_filename),mode=NF90_NOWRITE,ncid=ncid)
526 [ + + ]: 46 : if(ncerr /= NF90_NOERR) then
527 : 39 : write(std_out,'(3a)') '-'//'Could not open ',trim(hist_filename),', starting from scratch'
528 : 39 : Invar%netcdf=.false.
529 : : else
530 : 7 : write(std_out,'(3a)') '-'//'Succesfully open ',trim(hist_filename),' for reading'
531 : 7 : write(std_out,'(a)') ' Extracting information from NetCDF file...'
532 : 7 : Invar%netcdf=.true.
533 : : end if
534 : :
535 [ + + ]: 46 : if ( Invar%netcdf) then
536 : : call get_dims_hist(ncid,Invar%natom,Invar%ntypat,nimage,mdtime,&
537 : 7 : & natom_id,ntypat_id,nimage_id,time_id,xyz_id,six_id,has_nimage)
538 [ - + + - : 30 : ABI_CALLOC(Invar%amu,(Invar%ntypat))
+ - + - -
+ - + - +
+ + ]
539 [ - + + - : 901 : ABI_CALLOC(Invar%typat,(Invar%natom))
+ - + - -
+ - + - +
+ + ]
540 [ - + + - : 30 : ABI_CALLOC(Invar%znucl,(Invar%ntypat))
+ - + - -
+ - + - +
+ + ]
541 : 7 : call read_csts_hist(ncid,dtion,Invar%typat,Invar%znucl,Invar%amu)
542 : :
543 : : ! Need to close NetCDF file because it is going to be reopened by read_md_hist
544 : 7 : ncerr = nf90_close(ncid)
545 : : ! .false. -> Velocities note used
546 : : ! .true. -> acell and rprimd may change (2017_04 only NVT/isoK used but maybe
547 : : ! .false. -> read all times
548 : : ! NPT one day ?)
549 : 21 : call read_md_hist(hist_filename,Hist,.false.,.true.,.false.)
550 : : end if
551 : : #endif
552 : :
553 : : ! =========================================================================== !
554 : : ! Read input file
555 : :
556 : 46 : string = repeat(" ", strlen)
557 : 46 : raw_string = repeat(" ", strlen)
558 : 46 : call instrng(input_filename, lenstr, 1, strlen, string, raw_string)
559 : : ! To make case-insensitive, map characters to upper case.
560 : 46 : call inupper(string(1:lenstr))
561 : :
562 : : ! marr is the aximum array size. It is thus a hard-coded maximum value
563 : : ! for 3 * (number of atoms)
564 : 46 : marr = 9000
565 [ - + ]: 46 : ABI_MALLOC(intarr, (marr))
566 [ - + ]: 46 : ABI_MALLOC(dprarr, (marr))
567 : 46 : jdtset = 0
568 : :
569 : : ! Mandatory input variables
570 : : ! -------------------------
571 : :
572 : : ! Bravais lattice
573 : 46 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'brav', tread, 'INT')
574 [ - + ]: 46 : if (tread == 0) then
575 : : write(msg,*)&
576 : 0 : 'Variable "brav" is mandatory, but was not found in input file.'
577 : 0 : ABI_ERROR(msg)
578 : : end if
579 [ + + ]: 138 : Invar%bravais(1:2) = intarr(1:2)
580 : :
581 : : ! Angle alpha
582 [ + + ]: 46 : if ((Invar%bravais(1).eq.2).or.(Invar%bravais(1).eq.5)) then
583 : 4 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'angle', tread, 'DPR')
584 [ - + ]: 4 : if (tread == 0) then
585 : : write(msg,*)&
586 : 0 : 'Variable "angle" is mandatory for this bravais lattice,',ch10,&
587 : 0 : 'but was not found in input file.'
588 : 0 : ABI_ERROR(msg)
589 : : end if
590 : 4 : Invar%angle_alpha = dprarr(1)
591 : : else
592 : 42 : Invar%angle_alpha=90.d0
593 : : end if
594 : :
595 : : ! natom_unitcell
596 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'natom_unitcell', tread, 'INT')
597 [ - + ]: 46 : if (tread == 0) then
598 : : write(msg,*)&
599 : 0 : 'Variable "natom_unitcell" is mandatory, but was not found in input file.'
600 : 0 : ABI_ERROR(msg)
601 : : end if
602 : 46 : Invar%natom_unitcell = intarr(1)
603 [ - + ]: 46 : if (3*Invar%natom_unitcell.gt.marr) then
604 : : write(msg,*)&
605 : 0 : 'Maximum number of atoms exceeded. Modify source code to circumvent problem.'
606 : 0 : ABI_ERROR(msg)
607 : : end if
608 : :
609 : : ! xred_unitcell
610 [ - + + - : 582 : ABI_CALLOC(Invar%xred_unitcell,(3,Invar%natom_unitcell))
+ - + - -
+ - + - +
+ + + + ]
611 : 46 : call intagm(dprarr, intarr, jdtset, marr, 3*Invar%natom_unitcell, string(1:lenstr), 'xred_unitcell', tread, 'DPR')
612 [ - + ]: 46 : if (tread == 0) then
613 : : write(msg,*)&
614 : 0 : 'Variable "xred_unitcell" is mandatory, but was not found in input file.'
615 : 0 : ABI_ERROR(msg)
616 : : end if
617 [ + + + + : 582 : Invar%xred_unitcell(:,:) = reshape(dprarr(1:3*Invar%natom_unitcell),(/3,Invar%natom_unitcell/))
+ + ]
618 [ + + ]: 184 : do ii=1,3
619 [ + + ]: 517 : do iatcell=1,Invar%natom_unitcell
620 [ + - + + ]: 471 : if ((Invar%xred_unitcell(ii,iatcell).le.(-0.5)).or.(Invar%xred_unitcell(ii,iatcell).gt.(0.5))) then
621 [ - + ]: 3 : do while (Invar%xred_unitcell(ii,iatcell).le.(-0.5))
622 : 0 : Invar%xred_unitcell(ii,iatcell)=Invar%xred_unitcell(ii,iatcell)+1.d0
623 : : end do
624 [ + + ]: 6 : do while (Invar%xred_unitcell(ii,iatcell).gt.(0.5))
625 : 6 : Invar%xred_unitcell(ii,iatcell)=Invar%xred_unitcell(ii,iatcell)-1.d0
626 : : end do
627 : : end if
628 : : end do
629 : : end do
630 : :
631 : : ! typat_unitcell
632 [ - + + - : 249 : ABI_CALLOC(Invar%typat_unitcell,(Invar%natom_unitcell))
+ - + - -
+ - + - +
+ + ]
633 : 46 : call intagm(dprarr, intarr, jdtset, marr, Invar%natom_unitcell, string(1:lenstr), 'typat_unitcell', tread, 'INT')
634 [ - + ]: 46 : if (tread == 0) then
635 : : write(msg,*)&
636 : 0 : 'Variable "typat_unitcell" is mandatory, but was not found in input file.'
637 : 0 : ABI_ERROR(msg)
638 : : end if
639 [ + + ]: 157 : Invar%typat_unitcell(:) = intarr(1:Invar%natom_unitcell)
640 : :
641 [ + + ]: 46 : if (Invar%netcdf) then
642 [ + + + + ]: 91 : Invar%rprimd_md(:,:)=TRANSPOSE(Hist%rprimd(:,:,Hist%ihist))
643 : : else
644 : : ! ntypat
645 : 39 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'ntypat', tread, 'INT')
646 [ + - ]: 39 : if (tread == 1) then
647 : 39 : Invar%ntypat = intarr(1)
648 : : else
649 : : write(msg,*)&
650 : 0 : 'The NetCDF file .nc is not used.',ch10,&
651 : 0 : 'The variable "ntypat" is thus mandatory.',ch10,&
652 : 0 : 'ACTION : Please modify your input file'
653 : 0 : ABI_ERROR(msg)
654 : : end if
655 : :
656 : : ! amu
657 [ - + + - : 166 : ABI_MALLOC(Invar%amu,(Invar%ntypat)); Invar%amu(:)=zero
+ - + - -
+ - + - +
+ + ]
658 : 39 : call intagm(dprarr, intarr, jdtset, marr, Invar%ntypat, string(1:lenstr), 'amu', tread, 'DPR')
659 [ + - ]: 39 : if (tread == 1) then
660 [ + - - + : 127 : Invar%amu = dprarr(1:Invar%ntypat)
- - + + ]
661 : : else
662 : : write(msg,*)&
663 : 0 : 'The NetCDF file .nc is not used.',ch10,&
664 : 0 : 'The variable "amu" is thus mandatory.',ch10,&
665 : 0 : 'ACTION : Please modify your input file'
666 : 0 : ABI_ERROR(msg)
667 : : end if
668 : :
669 : : ! rprimd_md
670 : 39 : call intagm(dprarr, intarr, jdtset, marr, 9, string(1:lenstr), 'rprimd', tread, 'LEN')
671 [ + - ]: 39 : if (tread == 1) then
672 [ + + + + ]: 507 : Invar%rprimd_md(:,:) = TRANSPOSE(reshape(dprarr(1:9),(/3, 3/)))
673 : : else
674 : : write(msg,*)&
675 : 0 : 'The NetCDF file .nc is not used.',ch10,&
676 : 0 : 'The variable "rprimd" is thus mandatory.',ch10,&
677 : 0 : 'ACTION : Please modify your input file'
678 : 0 : ABI_ERROR(msg)
679 : : end if
680 : :
681 : : ! natom
682 : 39 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'natom', tread, 'INT')
683 [ + - ]: 39 : if (tread == 1) then
684 : 39 : Invar%natom = intarr(1)
685 : : else
686 : : write(msg,*)&
687 : 0 : 'The NetCDF file .nc is not used.',ch10,&
688 : 0 : 'The variable "natom" is thus mandatory.',ch10,&
689 : 0 : 'ACTION : Please modify your input file'
690 : 0 : ABI_ERROR(msg)
691 : : end if
692 : :
693 : : ! typat
694 [ - + + - : 6515 : ABI_MALLOC(Invar%typat,(Invar%natom)); Invar%typat(:)=0
+ - + - -
+ - + - +
+ + ]
695 : 39 : call intagm(dprarr, intarr, jdtset, marr, Invar%natom, string(1:lenstr), 'typat', tread, 'INT')
696 [ + - ]: 39 : if (tread == 1) then
697 [ + - - + : 6476 : Invar%typat = intarr(1:Invar%natom)
- - + + ]
698 : : else
699 : : write(msg,*)&
700 : 0 : 'The NetCDF file .nc is not used.',ch10,&
701 : 0 : 'The variable "typat" is thus mandatory.',ch10,&
702 : 0 : 'ACTION : Please modify your input file'
703 : 0 : ABI_ERROR(msg)
704 : : end if
705 : :
706 : : end if
707 : :
708 : : ! multiplicity
709 : 46 : call intagm(dprarr, intarr, jdtset, marr, 9, string(1:lenstr), 'multiplicity', tread, 'DPR')
710 [ + - ]: 46 : if (tread == 1) then
711 [ + + + + ]: 598 : Invar%multiplicity(:,:) = TRANSPOSE(reshape(dprarr(1:9),(/3, 3/)))
712 : : else
713 : : write(msg,*)&
714 : 0 : 'Variable "multiplicity" is mandatory, but was not found in input file.'
715 : 0 : ABI_ERROR(msg)
716 : : end if
717 : :
718 : : ! temperature
719 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'temperature', tread, 'DPR')
720 [ + - ]: 46 : if (tread == 1) then
721 : 46 : Invar%temperature = dprarr(1)
722 : : else
723 : : write(msg,*)&
724 : 0 : 'Variable "temperature" is mandatory, but was not found in input file.'
725 : 0 : ABI_ERROR(msg)
726 : : end if
727 : :
728 : : ! nstep_max
729 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'nstep_max', tread, 'INT')
730 [ + - ]: 46 : if (tread == 1) then
731 : 46 : Invar%nstep_max = intarr(1)
732 : : else
733 : : write(msg,*)&
734 : 0 : 'Variable "nstep_max" is mandatory, but was not found in input file.'
735 : 0 : ABI_ERROR(msg)
736 : : end if
737 : :
738 : : ! nstep_min
739 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'nstep_min', tread, 'INT')
740 [ + - ]: 46 : if (tread == 1) then
741 : 46 : Invar%nstep_min = intarr(1)
742 : : else
743 : : write(msg,*)&
744 : 0 : 'Variable "nstep_min" is mandatory, but was not found in input file.'
745 : 0 : ABI_ERROR(msg)
746 : : end if
747 : :
748 : : ! rcut
749 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'rcut', tread, 'LEN')
750 [ + - ]: 46 : if (tread == 1) then
751 : 46 : Invar%rcut = dprarr(1)
752 : : else
753 : : write(msg,*)&
754 : 0 : 'Variable "rcut" is mandatory, but was not found in input file.'
755 : 0 : ABI_ERROR(msg)
756 : : end if
757 : :
758 : : ! znucl (optional)
759 [ + + ]: 46 : if (.not.Invar%netcdf) then
760 : 39 : call intagm(dprarr, intarr, jdtset, marr, Invar%ntypat, string(1:lenstr), 'znucl', tread, 'DPR')
761 [ + - ]: 39 : if (tread == 1) then
762 [ - + + - : 166 : ABI_MALLOC(Invar%znucl,(Invar%ntypat)); Invar%znucl(:)=zero
+ - + - -
+ - + - +
+ + ]
763 [ - + + + ]: 127 : Invar%znucl = dprarr(1:Invar%ntypat)
764 : : end if
765 : : end if
766 : :
767 : : ! =========================================================================== !
768 : : ! Read DDB file if available to retrieve the Born effective charges
769 : : ! and dielectric tensor. Enforce charge neutrality by equal redistribution of excess charge
770 : : ! =========================================================================== !
771 : :
772 : : ! symbec
773 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'symbec', tread, 'INT')
774 [ - + ]: 46 : if (tread == 1) then
775 : 0 : Invar%symbec = intarr(1)
776 : : end if
777 : :
778 : :
779 [ - + + - : 138 : ABI_MALLOC(Invar%zeff, (3,3,Invar%natom_unitcell))
+ - + - -
+ - + -
+ ]
780 [ + + + + : 1489 : Invar%zeff = zero
+ + ]
781 : :
782 [ + - - + ]: 46 : if (file_exists(ddb_filename).or.file_exists(nctk_ncify(ddb_filename))) then
783 : :
784 : 0 : call ddb%from_file(ddb_filename, ddb_hdr, crystal, xmpi_world, -1)
785 : : ! TODO Check that the crystal in the DDB matches the unit cell.
786 : :
787 : 0 : ii = ddb%get_dielt_zeff(crystal,1,1,0,Invar%dielt,Invar%zeff)
788 [ # # ]: 0 : if (ii/=0) Invar%dipdip = .true.
789 : :
790 : 0 : call crystal%free()
791 : 0 : call ddb_hdr%free()
792 : 0 : call ddb%free()
793 : :
794 : : !! Take the average born effective charge for each atom.
795 [ - - ]: 46 : if (Invar%symbec==1) then
796 [ # # ]: 0 : do iatom=1,Invar%natom_unitcell
797 : : born_average = zero
798 [ # # ]: 0 : do ii=1,3
799 : 0 : born_average = born_average + Invar%zeff(ii,ii,iatom)
800 : : end do
801 [ # # # # ]: 0 : Invar%zeff(:,:,iatom) = zero
802 [ # # ]: 0 : do ii=1,3
803 : 0 : Invar%zeff(ii,ii,iatom) = born_average / 3
804 : : end do
805 : : end do
806 : : end if
807 : : end if
808 : :
809 : :
810 : : ! =========================================================================== !
811 : : ! Optional input variables
812 : : ! =========================================================================== !
813 : :
814 : : ! prtvol
815 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'prtvol', tread, 'INT')
816 [ - + ]: 46 : if (tread == 1) then
817 : 0 : Invar%prtvol = intarr(1)
818 : : end if
819 : :
820 : : ! prtphdos
821 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'prtphdos', tread, 'INT')
822 [ - + ]: 46 : if (tread == 1) then
823 : 0 : Invar%prtphdos = intarr(1)
824 : : end if
825 : :
826 : : ! dosdeltae
827 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'dosdeltae', tread, 'ENE')
828 [ - + ]: 46 : if (tread == 1) then
829 : 0 : Invar%dosdeltae = dprarr(1)
830 : : end if
831 : :
832 : : ! use_ideal_positions
833 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'use_ideal_positions', tread, 'INT')
834 [ + + ]: 46 : if (tread == 1) then
835 : 42 : Invar%use_ideal_positions = intarr(1)
836 : : end if
837 : :
838 : : ! born_charge
839 : 46 : call intagm(dprarr, intarr, jdtset, marr, Invar%ntypat, string(1:lenstr), 'born_charge', tread, 'DPR')
840 [ + + ]: 46 : if (tread == 1) then
841 : 5 : Invar%dipdip=.true.
842 [ - + + - : 25 : ABI_MALLOC(born_charge_tmp,(Invar%ntypat)) ; born_charge_tmp(:) = zero
+ - + - -
+ - + +
+ ]
843 [ + + ]: 15 : born_charge_tmp(:) = dprarr(1:Invar%ntypat)
844 [ + + + + : 148 : Invar%zeff(:,:,:) = zero
+ + ]
845 [ + + ]: 16 : do iatom=1,Invar%natom_unitcell
846 : 11 : itypat = Invar%typat_unitcell(iatom)
847 [ + + ]: 49 : do ii=1,3
848 : 44 : Invar%zeff(ii,ii,iatom) = born_charge_tmp(itypat)
849 : : end do
850 : : end do
851 [ - + ]: 5 : ABI_FREE(born_charge_tmp)
852 : : end if
853 : :
854 : : ! dielec_constant
855 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'dielec_constant', tread, 'DPR')
856 [ + + ]: 46 : if (tread == 1) then
857 [ + + ]: 20 : do ii=1,3
858 : 20 : Invar%dielt(ii,ii) = dprarr(1)
859 : : end do
860 : : end if
861 : :
862 : : ! dipdip
863 : : ! Invar%dipdip is auto-enabled above whenever usable Born effective
864 : : ! charges/dielectric tensor are found (from the DDB file or born_charge).
865 : : ! Explicitly setting dipdip=0 overrides that and forces it off.
866 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'dipdip', tread, 'INT')
867 [ - + ]: 46 : if (tread == 1) then
868 [ # # ]: 0 : if (intarr(1) == 0) Invar%dipdip = .false.
869 : : end if
870 : :
871 : : ! bzpath
872 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'bzpath', tread, 'INT')
873 [ + + ]: 46 : if (tread == 1) then
874 : 10 : Invar%bzpath = intarr(1)
875 : :
876 [ - + ]: 10 : if (Invar%bzpath.lt.0) then
877 [ # # # # : 0 : ABI_MALLOC(Invar%qpt,(3,abs(Invar%bzpath))); Invar%qpt(:,:)=zero
# # # # #
# # # # #
# # ]
878 : :
879 : 0 : call intagm(dprarr, intarr, jdtset, marr, 1+3*abs(Invar%bzpath), string(1:lenstr), 'bzpath', tread, 'DPR')
880 [ # # # # : 0 : Invar%qpt(:,:) = reshape(dprarr(2:1+3*abs(Invar%bzpath)), (/3,abs(Invar%bzpath)/))
# # ]
881 : :
882 [ + - ]: 10 : else if (Invar%bzpath.gt.0) then
883 [ - + + - : 20 : ABI_MALLOC(Invar%special_qpt,(Invar%bzpath))
- + - + -
+ ]
884 : 10 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'special_qpt', tread, 'KEY', key_value=msg)
885 [ + - ]: 10 : if (tread == 1) then
886 : 10 : jj = 1
887 [ + + ]: 64 : do ii=1,Invar%bzpath
888 : 54 : Invar%special_qpt(ii) = " "
889 : 64 : ierr = next_token(msg,jj,Invar%special_qpt(ii))
890 : : !jj = 2*ii-1
891 : : !Invar%special_qpt(ii) = msg(jj:jj)
892 : : end do
893 : : else
894 : : write(msg,*)&
895 : 0 : 'Variable bzpath > 0, but variable "special_qpt" was not found in input file.'
896 : 0 : ABI_ERROR(msg)
897 : : end if
898 : : end if
899 : : end if
900 : :
901 : : ! ndivsm
902 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'ndivsm', tread, 'INT')
903 [ - + ]: 46 : if (tread == 1) then
904 : 0 : Invar%ndivsm = intarr(1)
905 : : end if
906 : :
907 : : ! order
908 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'order', tread, 'INT')
909 [ + + ]: 46 : if (tread == 1) then
910 : 16 : Invar%order = intarr(1)
911 [ + + ]: 16 : if (Invar%order.eq.3) then
912 : 7 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'order', tread, 'DPR')
913 : 7 : Invar%rcut3 = dprarr(2)
914 [ - + ]: 7 : if (Invar%rcut3.gt.Invar%rcut) then
915 : 0 : ABI_ERROR('The cutoff radius of the third order cannot be greater than the second order one.')
916 : : end if
917 [ + - ]: 9 : else if (Invar%order.eq.4) then
918 : 9 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'order', tread, 'DPR')
919 : 9 : Invar%rcut3 = dprarr(2)
920 : 9 : Invar%rcut4 = dprarr(3)
921 [ - + ]: 9 : if (Invar%rcut4.gt.Invar%rcut) then
922 : 0 : ABI_ERROR('The cutoff radius of the fourth order cannot be greater than the second order one.')
923 : : end if
924 : : else
925 : 0 : ABI_ERROR('Only the 3rd and 4th orders are allowed. Change your input file.')
926 : : end if
927 : : end if
928 : :
929 : : ! slice
930 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'slice', tread, 'INT')
931 [ + + ]: 46 : if (tread == 1) then
932 : 2 : Invar%slice = intarr(1)
933 : : end if
934 : :
935 : : ! enunit
936 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'enunit', tread, 'INT')
937 [ + + ]: 46 : if (tread == 1) then
938 : 38 : Invar%enunit = intarr(1)
939 : : end if
940 : :
941 : : ! nproc
942 : 46 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'nproc', tread, 'INT')
943 [ - + ]: 46 : if (tread == 1) then
944 [ # # ]: 0 : Invar%nproc(1:2) = intarr(1:2)
945 : : end if
946 : :
947 : : ! readifc
948 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'readifc', tread, 'INT')
949 [ + + ]: 46 : if (tread == 1) then
950 : 3 : Invar%readifc = intarr(1)
951 [ + + ]: 3 : if (Invar%readifc.eq.1) then
952 : 1 : call intagm(dprarr, intarr, jdtset, marr, 2, string(1:lenstr), 'readifc', tread, 'DPR')
953 : 1 : Invar%tolread = dprarr(1)
954 : : end if
955 : : end if
956 : :
957 : : ! alloy
958 : 46 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'alloy', tread, 'INT')
959 [ + + ]: 46 : if (tread == 1) then
960 : 1 : Invar%alloy = intarr(1)
961 : 1 : Invar%ityp_alloy1 = intarr(2)
962 : 1 : Invar%ityp_alloy2 = intarr(3)
963 : : end if
964 : :
965 : : ! together
966 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'together', tread, 'INT')
967 [ + + ]: 46 : if (tread == 1) then
968 : 2 : Invar%together = intarr(1)
969 : : end if
970 : :
971 : : ! bzlength
972 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'bzlength', tread, 'INT')
973 [ + + ]: 46 : if (tread == 1) then
974 : 2 : Invar%bzlength = intarr(1)
975 [ - + + - : 6 : ABI_MALLOC(Invar%lgth_segments,(Invar%bzlength))
+ - + - -
+ - + -
+ ]
976 : 2 : call intagm(dprarr, intarr, jdtset, marr, 1+Invar%bzlength, string(1:lenstr), 'bzlength', tread, 'DPR')
977 [ + - - + : 14 : Invar%lgth_segments = dprarr(2:1+Invar%bzlength)
- - + + ]
978 : : end if
979 : :
980 : : ! ngqpt1
981 : 46 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'ngqpt1', tread, 'INT')
982 [ + + ]: 46 : if (tread == 1) then
983 [ + + ]: 52 : Invar%ngqpt1 = intarr(1:3)
984 : : end if
985 : :
986 : : ! ngqpt2
987 : 46 : call intagm(dprarr, intarr, jdtset, marr, 3, string(1:lenstr), 'ngqpt2', tread, 'INT')
988 [ + + ]: 46 : if (tread == 1) then
989 [ + + ]: 92 : Invar%ngqpt2 = intarr(1:3)
990 : : end if
991 : :
992 : : ! use_weights
993 : 46 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'use_weights', tread, 'INT')
994 [ + + ]: 46 : if (tread == 1) then
995 : 1 : Invar%use_weights = intarr(1)
996 : : end if
997 : :
998 : 46 : ABI_FREE(intarr)
999 [ - + ]: 46 : ABI_FREE(dprarr)
1000 : :
1001 : : ! =========================================================================== !
1002 : : ! Perform some checks
1003 : :
1004 [ - + ]: 46 : if ((Invar%use_ideal_positions.ne.0).and.(Invar%use_ideal_positions.ne.1)) then
1005 : 0 : ABI_ERROR('STOP: THIS VALUE IS NOT ALLOWED FOR use_ideal_positions')
1006 : : end if
1007 : :
1008 : : ! Allowed values
1009 [ - + ]: 46 : if ((Invar%together.ne.1).and.(Invar%together.ne.0)) then
1010 : 0 : ABI_ERROR('STOP: The value of input variable TOGETHER is not allowed')
1011 : : end if
1012 [ - + ]: 46 : if ((Invar%alloy.ne.1).and.(Invar%alloy.ne.0)) then
1013 : 0 : ABI_ERROR('STOP: The value of input variable ALLOY is not allowed')
1014 : : end if
1015 [ + + ]: 46 : if (Invar%alloy.ge.1) then
1016 : : if ((Invar%ityp_alloy1.lt.1).or.(Invar%ityp_alloy2.lt.1).or.&
1017 [ + - + - : 1 : & (Invar%ityp_alloy1.gt.Invar%natom_unitcell).or.(Invar%ityp_alloy2.gt.Invar%natom_unitcell)) then
+ - - + ]
1018 : 0 : ABI_ERROR('STOP: The value of input variables IALLOY are not allowed')
1019 : : end if
1020 : : end if
1021 : :
1022 : : ! Incompatible variables :
1023 [ + + - + ]: 46 : if ((Invar%readifc.eq.1).and.(Invar%together.eq.1).and.(Invar%order.gt.2)) then
1024 : 0 : ABI_ERROR('STOP: readifc=1, together=1 and order=3 or 4 are incompatible')
1025 : : end if
1026 : :
1027 : : ! =========================================================================== !
1028 : : ! GA: TODO REMOVE ALLOYING FEATURE.
1029 : : ! We agreed this logic can be handled externally.
1030 : : !
1031 : : ! Treat virtual crystal approximation (VCA) when alloy=1.
1032 : : ! Redefine all the data depending on (n)typat(_unitcell) and natom_unitcell
1033 : :
1034 : : ! The mixing of zeff will only work for the isotropic case,
1035 : : ! that is, if it was specified with input variables born_charge
1036 : :
1037 [ + + ]: 46 : if (Invar%alloy.eq.1) then
1038 : 1 : sum_alloy1=0
1039 : 1 : sum_alloy2=0
1040 [ + + ]: 129 : do iatom=1,Invar%natom
1041 [ + + ]: 128 : if (Invar%typat(iatom).eq.Invar%ityp_alloy1) then
1042 : 96 : sum_alloy1=sum_alloy1+1
1043 : : end if
1044 [ + + ]: 129 : if (Invar%typat(iatom).eq.Invar%ityp_alloy2) then
1045 : 32 : sum_alloy2=sum_alloy2+1
1046 : : end if
1047 : : end do
1048 : : amu_average =(Invar%amu(Invar%ityp_alloy1)*sum_alloy1+&
1049 : 1 : & Invar%amu(Invar%ityp_alloy2)*sum_alloy2)/(sum_alloy1+sum_alloy2)
1050 [ - + ]: 1 : if (Invar%dipdip) then
1051 : 0 : zeff_average(:,:) = zero
1052 [ # # ]: 0 : do iatom=1,Invar%natom_unitcell
1053 [ # # ]: 0 : if (Invar%typat_unitcell(iatom).eq.Invar%ityp_alloy1) then
1054 [ # # # # ]: 0 : zeff_average(:,:) = zeff_average(:,:) + Invar%zeff(:,:,iatom) * sum_alloy1
1055 : : exit
1056 : : end if
1057 : : end do
1058 [ # # ]: 0 : do iatom=1,Invar%natom_unitcell
1059 [ # # ]: 0 : if (Invar%typat_unitcell(iatom).eq.Invar%ityp_alloy2) then
1060 [ # # # # ]: 0 : zeff_average(:,:) = zeff_average(:,:) + Invar%zeff(:,:,iatom) * sum_alloy2
1061 : : exit
1062 : : end if
1063 : : end do
1064 [ # # # # ]: 0 : zeff_average(:,:) = zeff_average(:,:) / (sum_alloy1 + sum_alloy2)
1065 [ # # ]: 0 : do iatom=1,Invar%natom_unitcell
1066 [ # # ]: 0 : if (Invar%typat_unitcell(iatom).eq.min(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1067 [ # # # # ]: 0 : Invar%zeff(:,:,iatom) = zeff_average(:,:)
1068 : : end if
1069 : : end do
1070 : : end if
1071 : 1 : shift=0
1072 [ + + ]: 3 : do iatom=1,Invar%natom_unitcell
1073 [ + + ]: 3 : if (Invar%typat_unitcell(iatom).lt.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1074 : 1 : Invar%typat_unitcell (iatom-shift)=Invar%typat_unitcell (iatom)
1075 [ + + ]: 4 : Invar%xred_unitcell(:,iatom-shift)=Invar%xred_unitcell(:,iatom)
1076 [ - + ]: 1 : if (Invar%dipdip) then
1077 [ # # # # ]: 0 : Invar%zeff(:,:,iatom-shift) = Invar%zeff(:,:,iatom)
1078 : : end if
1079 [ + - ]: 1 : else if (Invar%typat_unitcell(iatom).eq.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1080 : 1 : shift=shift+1
1081 : : else if (Invar%typat_unitcell(iatom).gt.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1082 : 0 : Invar%typat_unitcell (iatom-shift)=Invar%typat_unitcell (iatom) - 1
1083 [ # # ]: 0 : Invar%xred_unitcell(:,iatom-shift)=Invar%xred_unitcell(:,iatom)
1084 [ # # ]: 0 : if (Invar%dipdip) then
1085 [ # # # # ]: 0 : Invar%zeff(:,:,iatom-shift) = Invar%zeff(:,:,iatom)
1086 : : end if
1087 : : end if
1088 : : end do
1089 : 1 : Invar%natom_unitcell=Invar%natom_unitcell-shift
1090 [ + + ]: 129 : do iatom=1,Invar%natom
1091 [ + + ]: 129 : if (Invar%typat(iatom).ge.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1092 : 32 : Invar%typat(iatom)=Invar%typat(iatom) - 1
1093 : : end if
1094 : : end do
1095 [ + + ]: 3 : do itypat=1,Invar%ntypat
1096 [ + + ]: 3 : if (itypat.eq.min(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1097 : 1 : Invar%amu (itypat)=amu_average
1098 [ - + ]: 1 : else if (itypat.gt.max(Invar%ityp_alloy1,Invar%ityp_alloy2)) then
1099 : 0 : Invar%amu (itypat-1)=Invar%amu (itypat)
1100 : : end if
1101 : : end do
1102 : 1 : Invar%ntypat=Invar%ntypat-1
1103 [ - + + - : 3 : ABI_MALLOC(typat_unitcell_tmp,( Invar%natom_unitcell))
+ - + - -
+ - + ]
1104 [ - + + - : 3 : ABI_MALLOC(xred_unitcell_tmp ,(3,Invar%natom_unitcell))
+ - + - -
+ - + ]
1105 [ - + + - : 3 : ABI_MALLOC(amu_tmp ,( Invar%ntypat))
+ - + - -
+ - + ]
1106 [ - + ]: 1 : if (Invar%dipdip) then
1107 [ # # # # : 0 : ABI_MALLOC(zeff_tmp ,(3,3,Invar%natom_unitcell))
# # # # #
# # # ]
1108 : : end if
1109 [ + + ]: 2 : typat_unitcell_tmp (:)=Invar%typat_unitcell(1:Invar%natom_unitcell)
1110 [ + + + + ]: 5 : xred_unitcell_tmp(:,:)=Invar%xred_unitcell(:,1:Invar%natom_unitcell)
1111 [ + + ]: 2 : amu_tmp (:)=Invar%amu(1:Invar%ntypat)
1112 [ - + ]: 1 : if (Invar%dipdip) then
1113 [ # # # # : 0 : zeff_tmp(:,:,:) = Invar%zeff(:,:,1:Invar%natom_unitcell)
# # ]
1114 : : end if
1115 [ + - - + : 3 : ABI_REMALLOC(Invar%typat_unitcell,( Invar%natom_unitcell))
+ - + - +
- - + - +
- + ]
1116 [ + - - + : 3 : ABI_REMALLOC(Invar%xred_unitcell ,(3,Invar%natom_unitcell))
+ - + - +
- - + - +
- + ]
1117 [ + - - + : 3 : ABI_REMALLOC(Invar%amu ,( Invar%ntypat))
+ - + - +
- - + - +
- + ]
1118 [ + - - + : 3 : ABI_REMALLOC(Invar%zeff ,(3,3,Invar%natom_unitcell))
+ - + - +
- - + - +
- + ]
1119 [ + + ]: 2 : Invar%typat_unitcell (:)=typat_unitcell_tmp (:)
1120 [ + + + + ]: 5 : Invar%xred_unitcell(:,:)=xred_unitcell_tmp(:,:)
1121 [ + + ]: 2 : Invar%amu (:)=amu_tmp (:)
1122 [ - + ]: 1 : if (Invar%dipdip) then
1123 [ # # # # : 0 : Invar%zeff(:,:,:) =zeff_tmp(:,:,:)
# # ]
1124 : : end if
1125 [ + - ]: 1 : if (allocated(Invar%znucl)) then
1126 [ - + + - : 3 : ABI_MALLOC(znucl_tmp,(Invar%ntypat))
+ - + - -
+ - + ]
1127 [ + + ]: 2 : znucl_tmp(:) = Invar%znucl(1:Invar%ntypat)
1128 [ - + + - : 3 : ABI_REMALLOC(Invar%znucl,(Invar%ntypat))
+ - + - -
+ - + ]
1129 [ + + ]: 2 : Invar%znucl(:) = znucl_tmp(:)
1130 [ - + ]: 1 : ABI_FREE(znucl_tmp)
1131 : : end if
1132 [ - + ]: 1 : ABI_FREE(typat_unitcell_tmp)
1133 [ - + ]: 1 : ABI_FREE(xred_unitcell_tmp)
1134 [ - + ]: 1 : ABI_FREE(amu_tmp)
1135 [ - + ]: 1 : if (Invar%dipdip) then
1136 [ # # ]: 0 : ABI_FREE(zeff_tmp)
1137 : : end if
1138 : : end if
1139 : :
1140 : : ! Anchor the unitcell on atom 1: shift every atom by atom 1's own reduced
1141 : : ! position and wrap into [0,1). This is not just a local convenience: the
1142 : : ! symmetry finder (tdep_make_sym) needs atom 1 at the origin to find the
1143 : : ! space group, and this same convention is reused, unmodified, as the
1144 : : ! reference unitcell (Sym%xred_zero, Crystal%xred, SC%xred_unitcell) for the
1145 : : ! rest of the code. The original (as-given-in-input) positions are kept
1146 : : ! separately in xred_unitcell_original, so that the rigid shift between the
1147 : : ! two conventions can be applied explicitly wherever needed (see
1148 : : ! tdep_sampling_init_read).
1149 [ - + + - : 138 : ABI_MALLOC(Invar%xred_unitcell_original,(3,Invar%natom_unitcell))
+ - + - -
+ - + -
+ ]
1150 [ + + + + ]: 486 : Invar%xred_unitcell_original(:,:) = Invar%xred_unitcell(:,:)
1151 [ + + ]: 156 : do iatcell=1,Invar%natom_unitcell
1152 [ + + ]: 440 : Invar%xred_unitcell(:,iatcell) = Invar%xred_unitcell(:,iatcell) - Invar%xred_unitcell_original(:,1)
1153 [ + + ]: 486 : do ii=1,3
1154 [ + + ]: 330 : if (Invar%xred_unitcell(ii,iatcell).lt.0.d0) then
1155 : 66 : jj=int(Invar%xred_unitcell(ii,iatcell)+1.0)
1156 : : else
1157 : 264 : jj=int(Invar%xred_unitcell(ii,iatcell)+0.0)
1158 : : end if
1159 : 440 : Invar%xred_unitcell(ii,iatcell) = Invar%xred_unitcell(ii,iatcell) - real(jj)
1160 : : end do
1161 : : end do
1162 : :
1163 : : ! Compute Nstep as a function of the slice
1164 : 46 : Invar%nstep_tot=int(float(Invar%nstep_max-Invar%nstep_min)/float(Invar%slice)+1)
1165 : :
1166 [ - + - + : 92 : end subroutine tdep_read_input
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- - - - -
- - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
1167 : : !!***
1168 : :
1169 : : !====================================================================================================
1170 : :
1171 : : !!****f* ABINIT/m_tdep_dataset/tdep_write_input
1172 : : !! NAME
1173 : : !! tdep_write_input
1174 : : !!
1175 : : !! FUNCTION
1176 : : !! Echo the set of input variables to main output file.
1177 : : !!
1178 : : !! INPUTS
1179 : : !! Invar = Input object containing the input variables.
1180 : : !! unt = Unit number of the file to write in.
1181 : : !!
1182 : : !! OUTPUT
1183 : : !! Only writing.
1184 : : !!
1185 : : !! SOURCE
1186 : :
1187 : 46 : subroutine tdep_write_input(Invar, unt)
1188 : :
1189 : : class(atdep_dataset_type),intent(in) :: Invar
1190 : : integer, intent(in) :: unt
1191 : :
1192 : : integer :: ii,jj,kk
1193 : :
1194 : : ! *********************************************************************
1195 : :
1196 : 46 : write(unt,*) '#############################################################################'
1197 : 46 : write(unt,*) '######################### ECHO OF INPUT FILE ################################'
1198 : 46 : write(unt,*) '#############################################################################'
1199 : 46 : write(unt,*) ' '
1200 : 46 : write(unt,'(a)') ' ======================= Files ================================================'
1201 : 46 : write(unt,'(a1,a20,1x,a)')'-', ljust('indata_prefix',20),trim(Invar%input_prefix)
1202 : 46 : write(unt,'(a1,a20,1x,a)')'-', ljust('outdata_prefix',20),trim(Invar%output_prefix)
1203 : 46 : write(unt,'(a1,a20,1x,a)')'-', ljust('output_file',20),trim(Invar%output_file)
1204 : :
1205 : 46 : write(unt,'(a)') ' ======================= Define the unitcell ================================='
1206 : 46 : write(unt,'(1x,a20,1x,i4,1x,i4)') ljust('brav',20),Invar%bravais(1),Invar%bravais(2)
1207 [ + + ]: 46 : if ((Invar%bravais(1).eq.2).or.(Invar%bravais(1).eq.5)) then
1208 : 4 : write(unt,'(1x,a20,1x,f15.10)') 'angle',Invar%angle_alpha
1209 : : end if
1210 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('natom_unitcell',20),Invar%natom_unitcell
1211 : 46 : write(unt,'(1x,a20)') ljust('xred_unitcell',20)
1212 [ + + ]: 156 : do ii=1,Invar%natom_unitcell
1213 [ + - + + ]: 486 : write(unt,'(22x,3(f15.10,1x))') (Invar%xred_unitcell(jj,ii), jj=1,3)
1214 : : end do
1215 : : ! write(unt,'(1x,a20)') ljust('xred_unitcell_orig',20)
1216 : : ! do ii=1,Invar%natom_unitcell
1217 : : ! write(unt,'(22x,3(f15.10,1x))') (Invar%xred_unitcell_original(jj,ii), jj=1,3)
1218 : : ! end do
1219 : : ! write(unt,'(a)') &
1220 : : !& ' NOTE: The unitcell reduced coordinates (xred_unitcell) are shifted so that atom 1 sits at'
1221 : : ! write(unt,'(a)') &
1222 : : !& ' the origin (0,0,0); this convention is used throughout the rest of the code.'
1223 : : ! write(unt,'(a)') &
1224 : : !& ' xred_unitcell_original holds the positions exactly as given in the input file.'
1225 [ + + ]: 46 : if (Invar%alloy.eq.1) then
1226 : : write(unt,'(a)') &
1227 : 1 : & ' NOTE: Virtual crystal approximation applied (alloy=1); unitcell variables above'
1228 : : write(unt,'(a)') &
1229 : 1 : & ' reflect the merged (averaged) atom type.'
1230 : : end if
1231 [ + - + + ]: 156 : write(unt,'(1x,a20,20(1x,i4))') ljust('typat_unitcell',20),(Invar%typat_unitcell(jj),jj=1,Invar%natom_unitcell)
1232 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('ntypat',20),Invar%ntypat
1233 [ + - + + ]: 103 : write(unt,'(1x,a20,20(1x,f15.10))') ljust('amu',20),(Invar%amu(jj),jj=1,Invar%ntypat)
1234 [ + - ]: 46 : if (allocated(Invar%znucl)) then
1235 [ + - + + ]: 103 : write(unt,'(1x,a20,20(1x,f15.10))') ljust('znucl',20),(Invar%znucl(jj),jj=1,Invar%ntypat)
1236 : : end if
1237 : :
1238 : 46 : write(unt,'(a)') ' ======================= Define the supercell ================================'
1239 : 46 : write(unt,'(1x,a20)') ljust('rprimd',20)
1240 [ + + ]: 184 : do ii=1,3
1241 [ + - + + ]: 598 : write(unt,'(22x,3(f15.10,1x))') (Invar%rprimd_md(ii,jj),jj=1,3)
1242 : : end do
1243 : 46 : write(unt,'(1x,a20)') ljust('multiplicity',20)
1244 [ + + ]: 184 : do ii=1,3
1245 [ + - + + ]: 598 : write(unt,'(22x,3(f15.10,1x))') (Invar%multiplicity(ii,jj),jj=1,3)
1246 : : end do
1247 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('natom',20),Invar%natom
1248 : 46 : write(unt,'(1x,a20)') ljust('typat',20)
1249 [ - + + + ]: 792 : do ii=1,Invar%natom,10
1250 [ + + ]: 792 : if (ii+9.lt.Invar%natom) then
1251 [ + - + + ]: 7700 : write(unt,'(22x,10(i4,1x))') (Invar%typat(ii+jj-1),jj=1,10)
1252 : : else
1253 [ + - + + ]: 324 : write(unt,'(22x,10(i4,1x))') (Invar%typat(jj),jj=ii,Invar%natom)
1254 : : end if
1255 : : end do
1256 : :
1257 : 46 : write(unt,'(a)') ' ======================= Define computational details ========================'
1258 : 46 : write(unt,'(1x,a20,1x,i5)') ljust('nstep_max',20),Invar%nstep_max
1259 : 46 : write(unt,'(1x,a20,1x,i5)') ljust('nstep_min',20),Invar%nstep_min
1260 : 46 : write(unt,'(1x,a20,1x,i5)') ljust('nstep_tot',20),Invar%nstep_tot
1261 : 46 : write(unt,'(1x,a20,1x,f15.10)') ljust('rcut',20),Invar%rcut
1262 : 46 : write(unt,'(1x,a20,1x,f15.10)') ljust('temperature',20),Invar%temperature
1263 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('order',20),Invar%order
1264 [ + + ]: 46 : if (Invar%order.eq.3) then
1265 : 7 : write(unt,'(1x,a20,1x,f15.10)') ljust('rcut3',20),Invar%rcut3
1266 [ + + ]: 39 : else if (Invar%order.eq.4) then
1267 : 9 : write(unt,'(1x,a20,2(1x,f15.10))') ljust('rcut3/rcut4',20),Invar%rcut3,Invar%rcut4
1268 : : end if
1269 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('prtvol',20),Invar%prtvol
1270 : :
1271 : 46 : write(unt,'(a)') ' ======================= Optional input variables ============================'
1272 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('prtphdos',20),Invar%prtphdos
1273 : 46 : write(unt,'(1x,a20,1x,f15.10)') ljust('dosdeltae',20),Invar%dosdeltae
1274 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('use_ideal_positions',20),Invar%use_ideal_positions
1275 [ + + ]: 87 : write(unt,'(1x,a20,1x,i4)') ljust('dipdip',20),merge(1,0,Invar%dipdip)
1276 [ + + ]: 46 : if (Invar%dipdip) then
1277 : 5 : write(unt,'(1x,a20)') ljust('born_charge',20)
1278 [ + + ]: 16 : do kk=1,Invar%natom_unitcell
1279 : 11 : write(unt,'(24x,a,1x,i4)') 'Atom',kk
1280 [ + + ]: 49 : do ii=1,3
1281 [ + - + + ]: 143 : write(unt,'(22x,3(f15.10,1x))') (Invar%zeff(ii,jj,kk),jj=1,3)
1282 : : end do
1283 : : end do
1284 : 5 : write(unt,'(1x,a20)') ljust('dielec_constant',20)
1285 [ + + ]: 20 : do ii=1,3
1286 [ + - + + ]: 65 : write(unt,'(22x,3(f15.10,1x))') (Invar%dielt(ii,jj),jj=1,3)
1287 : : end do
1288 : 5 : write(unt,'(1x,a20,1x,i4)') ljust('symbec',20),Invar%symbec
1289 : : end if
1290 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('bzpath',20),Invar%bzpath
1291 [ - + ]: 46 : if (Invar%bzpath.lt.0) then
1292 : 0 : write(unt,'(a)') ' Q points as given in the input file:'
1293 [ # # ]: 0 : do jj=1,abs(Invar%bzpath)
1294 : 0 : write(unt,'(22x,3(f15.10,1x))') Invar%qpt(:,jj)
1295 : : end do
1296 [ + + ]: 46 : else if (Invar%bzpath.gt.0) then
1297 : 10 : write(unt,'(a,1x,a2,10("-",a2))') ' Special q-points: ',Invar%special_qpt(:)
1298 : : end if
1299 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('ndivsm',20),Invar%ndivsm
1300 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('slice',20),Invar%slice
1301 [ + + ]: 46 : if (Invar%enunit.eq.0) write(unt,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in meV)'
1302 [ + + ]: 46 : if (Invar%enunit.eq.1) write(unt,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in cm-1)'
1303 [ + + ]: 46 : if (Invar%enunit.eq.2) write(unt,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in mHa)'
1304 [ + + ]: 46 : if (Invar%enunit.eq.3) write(unt,'(1x,a20,1x,i4,1x,a)') ljust('enunit',20),Invar%enunit,'(Phonon frequencies in THz)'
1305 : 46 : write(unt,'(1x,a20,1x,i4,1x,i4)') ljust('nproc',20),Invar%nproc(1),Invar%nproc(2)
1306 [ + + ]: 46 : if (Invar%readifc.eq.1) then
1307 : 1 : write(unt,'(1x,a20,1x,i4,1x,f15.10)') ljust('readifc',20),Invar%readifc,Invar%tolread
1308 : : else
1309 : 45 : write(unt,'(1x,a20,1x,i4)') ljust('readifc',20),Invar%readifc
1310 : : end if
1311 : 46 : write(unt,'(1x,a20,1x,3(i4,1x))') ljust('alloy',20),Invar%alloy,Invar%ityp_alloy1,Invar%ityp_alloy2
1312 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('together',20),Invar%together
1313 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('bzlength',20),Invar%bzlength
1314 [ + + ]: 46 : if (Invar%bzlength.gt.0) then
1315 : 2 : write(unt,'(a,1x,i3,10("-",i3))') ' Length of BZ : ',Invar%lgth_segments(:)
1316 : : end if
1317 : 46 : write(unt,'(1x,a20,1x,3(i4,1x))') ljust('ngqpt1',20),Invar%ngqpt1(:)
1318 : 46 : write(unt,'(1x,a20,1x,3(i4,1x))') ljust('ngqpt2',20),Invar%ngqpt2(:)
1319 : 46 : write(unt,'(1x,a20,1x,i4)') ljust('use_weights',20),Invar%use_weights
1320 : :
1321 : 46 : end subroutine tdep_write_input
1322 : : !!***
1323 : :
1324 : : !====================================================================================================
1325 : :
1326 : : !!****f* ABINIT/m_tdep_dataset/tdep_destroy_invar
1327 : : !! NAME
1328 : : !! tdep_destroy_invar
1329 : : !!
1330 : : !! FUNCTION
1331 : : !! Free the memory allocated in the Invar object.
1332 : : !!
1333 : : !! SIDE EFFECTS
1334 : : !! Invar = Input object. All its allocated arrays are deallocated.
1335 : : !!
1336 : : !! SOURCE
1337 : :
1338 : 46 : subroutine tdep_destroy_invar(Invar)
1339 : :
1340 : : class(atdep_dataset_type), intent(inout) :: Invar
1341 : :
1342 [ - + ]: 46 : ABI_FREE(Invar%amu)
1343 [ - + ]: 46 : ABI_FREE(Invar%typat)
1344 [ - + ]: 46 : ABI_FREE(Invar%xred_unitcell)
1345 [ - + ]: 46 : ABI_FREE(Invar%xred_unitcell_original)
1346 [ - + ]: 46 : ABI_FREE(Invar%typat_unitcell)
1347 [ - + ]: 46 : if (Invar%bzpath.lt.0) then
1348 [ # # ]: 0 : ABI_FREE(Invar%qpt)
1349 [ + + ]: 46 : else if (Invar%bzpath.gt.0) then
1350 [ - + ]: 10 : ABI_FREE(Invar%special_qpt)
1351 : : end if
1352 [ + + ]: 46 : if (Invar%bzlength.gt.0) then
1353 [ - + ]: 2 : ABI_FREE(Invar%lgth_segments)
1354 : : end if
1355 [ + - ]: 46 : ABI_SFREE(Invar%zeff)
1356 [ + - ]: 46 : ABI_SFREE(Invar%znucl)
1357 : :
1358 : 46 : end subroutine tdep_destroy_invar
1359 : : !!***
1360 : :
1361 : : !====================================================================================================
1362 : :
1363 : : !!****f* ABINIT/m_tdep_dataset/tdep_init_MPIdata
1364 : : !! NAME
1365 : : !! tdep_init_MPIdata
1366 : : !!
1367 : : !! FUNCTION
1368 : : !! Initialize the MPIdata object: create the 2D cartesian grid of processors
1369 : : !! (nproc(1) processors for the shells times nproc(2) processors for the
1370 : : !! MD steps) and the corresponding communicators. If the product of
1371 : : !! nproc(1) and nproc(2) does not match the total number of processors,
1372 : : !! the parallelization is performed over the MD steps only.
1373 : : !! Then, distribute the nstep_tot MD steps over the processors, and write
1374 : : !! this distribution in the output file.
1375 : : !!
1376 : : !! OUTPUT
1377 : : !! MPIdata = Info on MPI parallelism, initialized.
1378 : : !!
1379 : : !! SIDE EFFECTS
1380 : : !! Invar = Input object containing the input variables.
1381 : : !! On exit, Invar%my_nstep is set.
1382 : : !!
1383 : : !! SOURCE
1384 : :
1385 : 46 : subroutine tdep_init_MPIdata(MPIdata,Invar)
1386 : :
1387 : : class(MPI_enreg_type), intent(out) :: MPIdata
1388 : : type(atdep_dataset_type), intent(inout) :: Invar
1389 : : integer :: ii,remain,ierr,iproc,istep
1390 : 46 : integer, allocatable :: tab_step(:)
1391 : : character(len=500) :: message
1392 : :
1393 : : #if defined HAVE_MPI
1394 : : integer :: dimcart,commcart_2d,me_cart_2d
1395 : : logical :: reorder
1396 : 46 : integer,allocatable :: coords(:),sizecart(:)
1397 [ + - ]: 46 : logical,allocatable :: periode(:), keepdim(:)
1398 : : #endif
1399 : :
1400 : : ! Check the number of processors
1401 : 46 : MPIdata%nproc_shell=Invar%nproc(1)
1402 : 46 : MPIdata%nproc_step =Invar%nproc(2)
1403 : 46 : MPIdata%nproc = xmpi_comm_size(xmpi_world)
1404 [ - + ]: 46 : if (MPIdata%nproc_step*MPIdata%nproc_shell.ne.MPIdata%nproc) then
1405 : 0 : ABI_WARNING('The parallelization is performed over steps')
1406 : 0 : MPIdata%nproc_step = xmpi_comm_size(xmpi_world)
1407 : : end if
1408 : :
1409 : 46 : MPIdata%master = 0
1410 : 46 : MPIdata%iam_master =.false.
1411 : : ! Initialize the MPIdata datastructure for sequential calculation
1412 [ + - ]: 46 : if (MPIdata%nproc.eq.1) then
1413 : 46 : MPIdata%comm_shell = xmpi_comm_null
1414 : 46 : MPIdata%comm_step = xmpi_comm_null
1415 : 46 : MPIdata%comm_shellstep = xmpi_comm_null
1416 : 46 : MPIdata%me_shell = 0
1417 : 46 : MPIdata%me_step = 0
1418 : 46 : MPIdata%my_nstep = Invar%nstep_tot
1419 : 46 : MPIdata%iam_master = (MPIdata%me_step == MPIdata%master)
1420 : : end if
1421 : :
1422 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1423 : : !!!!!!!!!!!!!!!!!!!!!! Parallel calculation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1424 : : !!!!!!!!!!!!!! Definition of the processor grid !!!!!!!!!!!!!!!!!!!!!!!!!!!!
1425 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1426 : : #if defined HAVE_MPI
1427 : : !FB if (MPIdata%nproc.eq.1) return
1428 : :
1429 : : ! Create the global cartesian 2D-communicator
1430 : 46 : dimcart=2
1431 [ - + ]: 46 : ABI_MALLOC(sizecart,(dimcart))
1432 [ - + ]: 46 : ABI_MALLOC(periode,(dimcart))
1433 : 46 : sizecart(1)=MPIdata%nproc_shell ! MPIdata%nproc_shell
1434 : 46 : sizecart(2)=MPIdata%nproc_step ! MPIdata%nproc_step
1435 [ + + ]: 138 : periode(:)=.false.;reorder=.false.
1436 : 46 : call MPI_CART_CREATE(xmpi_world,dimcart,sizecart,periode,reorder,commcart_2d,ierr)
1437 [ - + ]: 46 : ABI_FREE(periode)
1438 : 46 : ABI_FREE(sizecart)
1439 : :
1440 : : ! Find the index and coordinates of the current processor
1441 : 46 : call MPI_COMM_RANK(commcart_2d,me_cart_2d,ierr)
1442 [ - + + - : 138 : ABI_MALLOC(coords,(dimcart))
+ - + - -
+ - + ]
1443 : 46 : call MPI_CART_COORDS(commcart_2d,me_cart_2d,dimcart,coords,ierr)
1444 : 46 : MPIdata%me_shell=coords(1)
1445 : 46 : MPIdata%me_step =coords(2)
1446 : 46 : ABI_FREE(coords)
1447 [ + - + - ]: 46 : if ((MPIdata%me_shell == MPIdata%master).and.(MPIdata%me_step == MPIdata%master)) then
1448 : 46 : MPIdata%iam_master = .true.
1449 : : end if
1450 : :
1451 [ - + + - : 138 : ABI_MALLOC(keepdim,(dimcart))
+ - + - -
+ - + ]
1452 : : ! Create the communicator for shell distribution
1453 : 46 : keepdim(1)=.true.
1454 : 46 : keepdim(2)=.false.
1455 : 46 : call MPI_CART_SUB(commcart_2d,keepdim,MPIdata%comm_shell,ierr)
1456 : : ! Create the communicator for step distribution
1457 : 46 : keepdim(1)=.false.
1458 : 46 : keepdim(2)=.true.
1459 : 46 : call MPI_CART_SUB(commcart_2d,keepdim,MPIdata%comm_step,ierr)
1460 : : ! Create the communicator for shellstep distribution
1461 : 46 : keepdim(1)=.true.
1462 : 46 : keepdim(2)=.true.
1463 : 46 : call MPI_CART_SUB(commcart_2d,keepdim,MPIdata%comm_shellstep,ierr)
1464 : 46 : ABI_FREE(keepdim)
1465 : 46 : call xmpi_comm_free(commcart_2d)
1466 : :
1467 : : ! Write some data
1468 : 46 : write(message,'(a1,a20,2(1x,i4))')'-', 'Number of processors',MPIdata%nproc_shell,MPIdata%nproc_step
1469 : 46 : call wrtout(ab_out,message,'COLL')
1470 : : !FB write(message,'(a,2i5)') 'me_shell and me_step : ',MPIdata%me_shell,MPIdata%me_step
1471 : : !FB call wrtout(ab_out,message,'COLL')
1472 : : #endif
1473 : :
1474 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1475 : : !!!!!!!!!!!!!! Distribution over STEP processors !!!!!!!!!!!!!!!!!!!!!!!!!!!
1476 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1477 : 46 : MPIdata%my_nstep =int(Invar%nstep_tot/MPIdata%nproc_step)
1478 : 46 : remain=Invar%nstep_tot-MPIdata%nproc_step*MPIdata%my_nstep
1479 [ - + ]: 46 : do ii=1,remain
1480 [ - - ]: 46 : if ((ii-1).eq.MPIdata%me_step) MPIdata%my_nstep=MPIdata%my_nstep+1
1481 : : end do
1482 [ - + + - : 184 : ABI_MALLOC(MPIdata%nstep_all,(MPIdata%nproc_step)); MPIdata%nstep_all(:)=zero
+ - + - -
+ - + - +
+ + ]
1483 : 46 : call xmpi_allgather(MPIdata%my_nstep,MPIdata%nstep_all,MPIdata%comm_step,ierr)
1484 : 46 : write(ab_out,'(a)') ' '
1485 : 46 : write(ab_out,'(a,1x,i4)') ' All quantities are computed from nstep_min=',Invar%nstep_min
1486 : 46 : write(ab_out,'(a,1x,i4)') ' to nstep_max=',Invar%nstep_max
1487 [ + + ]: 46 : if (Invar%slice.ne.1) then
1488 : 2 : write(ab_out,'(a,1x,i4)') ' by using a slice=',Invar%slice
1489 : : end if
1490 : 46 : write(ab_out,'(a,1x,i4)') ' So, the real number of time steps is nstep=',Invar%nstep_tot
1491 [ - + ]: 46 : if (MPIdata%nproc_step.gt.1) then
1492 : 0 : write(ab_out,'(a,1000(1x,i5))') '-Distribution of number of steps wrt the number of processors=',MPIdata%nstep_all(:)
1493 : : end if
1494 : :
1495 [ + - - + : 230 : ABI_MALLOC(MPIdata%nstep_acc,(MPIdata%nproc_step+1)); MPIdata%nstep_acc(:)=zero
+ - + - +
- - + - +
- + + + ]
1496 : 46 : MPIdata%nstep_acc(1)=0
1497 [ + + ]: 92 : do ii=2,MPIdata%nproc_step+1
1498 : 92 : MPIdata%nstep_acc(ii)=MPIdata%nstep_acc(ii-1)+MPIdata%nstep_all(ii-1)
1499 : : end do
1500 [ - + ]: 46 : if (MPIdata%nstep_acc(MPIdata%nproc_step+1).ne.Invar%nstep_tot) then
1501 : 0 : write(std_out,*) 'STOP : pb in nstep_acc'
1502 : 0 : stop
1503 : : end if
1504 : :
1505 [ - + + - : 1130 : ABI_MALLOC(tab_step,(Invar%nstep_tot)); tab_step(:)=zero
+ - + - -
+ - + +
+ ]
1506 [ - + + - : 1084 : ABI_MALLOC(MPIdata%my_step ,(Invar%nstep_tot)); MPIdata%my_step (:)=.false.
- + - + -
+ + + ]
1507 [ + + ]: 92 : do iproc=1,MPIdata%nproc_step
1508 [ + + ]: 1084 : do istep=1,Invar%nstep_tot
1509 [ + - + - ]: 1038 : if ((istep.gt.MPIdata%nstep_acc(iproc)).and.(istep.le.MPIdata%nstep_acc(iproc+1))) then
1510 : 992 : tab_step(istep)=iproc-1
1511 : : end if
1512 : : end do
1513 : : end do
1514 [ + + ]: 1038 : do istep=1,Invar%nstep_tot
1515 : 1038 : MPIdata%my_step(istep) = (tab_step(istep) == MPIdata%me_step)
1516 : : end do
1517 : :
1518 [ - + + - : 184 : ABI_MALLOC(MPIdata%shft_step,(MPIdata%nproc_step)); MPIdata%shft_step(:)=zero
+ - + - -
+ - + - +
+ + ]
1519 : 46 : MPIdata%shft_step(1)=0
1520 [ - + ]: 46 : do ii=2,MPIdata%nproc_step
1521 : 46 : MPIdata%shft_step(ii)=MPIdata%shft_step(ii-1)+MPIdata%nstep_all(ii-1)
1522 : : end do
1523 [ - + ]: 46 : ABI_FREE(tab_step)
1524 : :
1525 : 46 : Invar%my_nstep=MPIdata%my_nstep
1526 : :
1527 : 46 : end subroutine tdep_init_MPIdata
1528 : : !!***
1529 : :
1530 : : !====================================================================================================
1531 : :
1532 : : !!****f* ABINIT/m_tdep_dataset/tdep_destroy_mpidata
1533 : : !! NAME
1534 : : !! tdep_destroy_mpidata
1535 : : !!
1536 : : !! FUNCTION
1537 : : !! Free the memory allocated in the MPIdata object.
1538 : : !!
1539 : : !! SIDE EFFECTS
1540 : : !! MPIdata = Info on MPI parallelism. Its allocated arrays are deallocated.
1541 : : !!
1542 : : !! SOURCE
1543 : :
1544 : 46 : subroutine tdep_destroy_mpidata(MPIdata)
1545 : :
1546 : : class(MPI_enreg_type), intent(inout) :: MPIdata
1547 : :
1548 [ - + ]: 46 : ABI_FREE(MPIdata%shft_step)
1549 [ - + ]: 46 : ABI_FREE(MPIdata%nstep_acc)
1550 [ - + ]: 46 : ABI_FREE(MPIdata%nstep_all)
1551 [ - + ]: 46 : ABI_FREE(MPIdata%my_step)
1552 : :
1553 : 46 : end subroutine tdep_destroy_mpidata
1554 : : !!***
1555 : :
1556 : : !====================================================================================================
1557 : :
1558 [ - + - - : 138 : end module m_tdep_dataset
+ + - + +
- - + - +
+ - - + -
+ - + - +
# # # # #
# # # ]
1559 : : !!***
|
