Branch data Line data Source code
1 : : !!****p*ABINIT/anaddb
2 : : !! NAME
3 : : !! anaddb
4 : : !!
5 : : !! FUNCTION
6 : : !! Main routine for analysis of the interatomic force constants and associated properties.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (XG,DCA,JCC,CL,XW,GA,MR)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : : !!
15 : : !! INPUTS
16 : : !! (main routine)
17 : : !!
18 : : !! OUTPUT
19 : : !! (main routine)
20 : : !!
21 : : !! SOURCE
22 : :
23 : : #if defined HAVE_CONFIG_H
24 : : #include "config.h"
25 : : #endif
26 : :
27 : : #include "abi_common.h"
28 : :
29 : 162 : program anaddb
30 : :
31 : : use, intrinsic :: iso_c_binding
32 : 162 : use defs_basis
33 : : use m_xmpi
34 : : use m_xomp
35 : : use m_abicore
36 : : use m_errors
37 : : use m_argparse
38 : : use m_nctk
39 : : use netcdf
40 : :
41 : : use m_build_info, only : abinit_version
42 : : use m_io_tools, only : open_file, flush_unit
43 : : use m_fstrings, only : int2char4, itoa, sjoin, strcat, inupper
44 : : use m_specialmsg, only : specialmsg_getcount, herald
45 : : use m_time, only : asctime, timein, timab, cwtime, cwtime_report
46 : : use m_dtfil, only : isfile
47 : : use m_crystal, only : crystal_t
48 : : use m_ddb, only : ddb_type, asrq0_t, ddb_lw_copy
49 : : use m_ddb_hdr, only : ddb_hdr_type
50 : : use m_ifc, only : ifc_type
51 : : use m_anaddb_dataset, only : anaddb_dataset_type
52 : : use m_anaddb_driver, only : anaddb_driver_type
53 : : use m_ddb_interpolate, only : ddb_interpolate
54 : : use m_elphon, only : elphon
55 : : use m_thmeig, only : thmeig
56 : : use m_symfind, only : symanal
57 : : use m_raman, only : ramansus, electrooptic
58 : : use m_ddb_diel, only : ddb_diel
59 : : use m_relaxpol, only : relaxpol
60 : : use m_ddb_elast, only : ddb_elast
61 : : use m_ddb_piezo, only : ddb_piezo
62 : : use m_ddb_internalstr, only : ddb_internalstr
63 : : use m_ddb_flexo, only : ddb_flexo
64 : : use m_ddb_magpen, only : ddb_magpen
65 : : use m_ddb_omega_interpol, only : ddb_omega_interpol
66 : : use m_lwf, only : run_lattice_wannier
67 : : use m_phonons, only : mkphbs
68 : : use m_gruneisen, only : gruns_anaddb
69 : :
70 : : implicit none
71 : :
72 : : !Local variables-------------------------------
73 : : integer, parameter:: master = 0
74 : : integer:: comm, ii, ierr
75 : : integer:: nproc, my_rank, ana_ncid
76 : : logical:: iam_master
77 : : real(dp):: tcpu, tcpui, twall, twalli !,cpu, wall, gflops
78 : : real(dp):: tsec(2)
79 : : integer:: units(2)
80 : : character(len=10):: procstr
81 : : character(len=24):: codename, start_datetime
82 : : ! character(len = strlen):: string, raw_string
83 : : character(len = fnlen):: worker_logfile
84 : : character(len=500):: msg
85 : : type(args_t):: args
86 [ + + ]: 486 : type(anaddb_dataset_type):: dtset
87 : 162 : type(anaddb_driver_type):: driver
88 [ + + + + : 8424 : type(crystal_t):: crystal
+ + + + +
+ + + + +
+ + + + ]
89 [ + + + + ]: 972 : type(ifc_type):: Ifc
90 : 162 : type(ddb_type):: ddb
91 : 162 : type(ddb_type):: ddb_lw
92 [ + + + + : 8424 : type(ddb_hdr_type):: ddb_hdr
+ + + + +
+ + + + +
+ + + + ]
93 : 162 : type(asrq0_t):: asrq0
94 : :
95 : : ! ========================================================================== !
96 : :
97 : : ! Change communicator for I/O (mandatory!)
98 : 162 : call abi_io_redirect(new_io_comm = xmpi_world)
99 : :
100 : : ! These units are defined in defs_basis
101 [ + + ]: 486 : units = [std_out, ab_out]
102 : :
103 : : ! Initialize MPI
104 : 162 : call xmpi_init()
105 : :
106 : : ! MPI variables
107 : 162 : comm = xmpi_world; nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
108 : 162 : iam_master = (my_rank == master)
109 : :
110 : : ! Parse command line arguments.
111 [ - + ]: 162 : args = args_parser(); if (args%exit /= 0) goto 100
112 : :
113 : : ! Initialize memory profiling if activated at configure time.
114 : : ! if a full report is desired, set the argument of abimem_init to "2" instead of "0" via the command line.
115 : : ! note that the file can easily be multiple GB in size so don't use this option normally
116 : : #ifdef HAVE_MEM_PROFILING
117 : : call abimem_init(args%abimem_level, limit_mb = args%abimem_limit_mb)
118 : : #endif
119 : :
120 : : ! Initialisation of the timing
121 : 162 : call timein(tcpui, twalli)
122 : :
123 [ + - ]: 162 : if (iam_master) then
124 : 162 : codename='ANADDB'//repeat(' ',18)
125 : 162 : call herald(codename, abinit_version, std_out)
126 : : end if
127 : :
128 : 162 : start_datetime = asctime()
129 : :
130 : : ! Zero out all accumulators of time and init timers
131 : 162 : call timab(1, 0, tsec)
132 : :
133 : : ! Initialise the code: write heading, and read names of files.
134 [ + - ]: 162 : if (iam_master) then
135 : 162 : call dtset%init(args%input_path)
136 : : end if
137 : :
138 : : ! Broadcast file names
139 : 162 : call dtset%bcast_files(comm)
140 : :
141 : : ! make log file for non-master procs
142 [ - + ]: 162 : if (.not. iam_master) then
143 : 0 : call int2char4(my_rank, procstr)
144 [ # # ]: 0 : ABI_CHECK((procstr(1:1)/='#'), 'Bug: string length too short!')
145 [ # # # # : 0 : worker_logfile = trim(dtset%filename_output) // "_LOG_P" // trim(procstr)
# # ]
146 [ # # ]: 0 : if (open_file(worker_logfile, msg, unit = std_out, form="formatted", action="write") /= 0) then
147 : 0 : ABI_ERROR(msg)
148 : : end if
149 : : end if
150 : :
151 : : ! ========================================================================== !
152 : : ! Read input variables
153 : 162 : call dtset%read_input(comm)
154 : :
155 [ - + ]: 162 : if (args%dry_run /= 0) then
156 : 0 : call wrtout(std_out, "Dry run mode. Exiting after have read the input")
157 : 0 : call dtset%free()
158 : 0 : goto 100
159 : : end if
160 : :
161 : : ! ========================================================================== !
162 : : ! Open output file
163 [ + - ]: 162 : if (iam_master) then
164 : 162 : call isfile(dtset%filename_output, 'new')
165 [ - + ]: 162 : if (open_file(dtset%filename_output, msg, unit=ab_out, form='formatted', status='new') /= 0) then
166 : 0 : ABI_ERROR(msg)
167 : : end if
168 : 162 : rewind (unit = ab_out)
169 : 162 : call herald(codename, abinit_version, ab_out)
170 : :
171 : : ! Echo the inputs to console and main output file
172 : 162 : call dtset%outvars(std_out)
173 : 162 : call dtset%outvars(ab_out)
174 : : else
175 : 0 : ab_out = dev_null
176 : : end if
177 : :
178 : : ! =========================================================================== !
179 : :
180 : : ! Initialize driver
181 : 162 : call driver%init(dtset)
182 : :
183 : : ! Read the DDB information and symmetrize partially the DDB
184 : 162 : write(msg, '(a, a)' )' read the DDB information and perform some checks',ch10
185 : 162 : call wrtout(units, msg)
186 : :
187 : 162 : call ddb%from_file(dtset%filename_ddb, ddb_hdr, crystal, comm, prtvol=dtset%prtvol)
188 : :
189 : : ! Change the bravais lattice if needed
190 : 162 : call ddb%set_brav(dtset%brav)
191 : : ! MR: a new ddb is necessary for the longwave quantities due to incompability of it with automatic reshapes
192 : : ! that ddb%val and ddb%flg experience when passed as arguments of some routines
193 : : ! Copy the long-wave ddb
194 [ + + ]: 162 : if (ddb_hdr%has_d3E_lw) then
195 : 16 : call ddb_lw_copy(ddb, ddb_lw, ddb_hdr)
196 : : end if
197 : :
198 : : ! Acoustic Sum Rule call
199 [ + + ]: 162 : if (dtset%flexoflag == 1 .and. dtset%asr == 6) then
200 : 2 : call driver%get_dcdq(dtset, ddb, ddb_lw, ddb_hdr)
201 : 2 : write(msg, '(a, a)' )' IFCs derivatives read',ch10
202 : 2 : call wrtout(units, msg)
203 : : end if
204 : : ! MR: Second- and third-order total energy derivatives calculated with the
205 : : ! magnetic penalty (constrained DFPT) are converted to physically relevant ones here.
206 [ + + ]: 162 : if (abs(dtset%magpen) > tol8) then
207 : : call ddb_magpen(ddb, ddb_lw, dtset%magpen, dtset%mpatpol, &
208 : : & dtset%mpdir, dtset%mpert, dtset%mpopt, Crystal%natom, dtset%prtvol, 1, Crystal%ucvol, dtset%timdisp, &
209 : 6 : & Crystal%xred)
210 : :
211 [ + - ]: 6 : if (dtset%freqflag/=0) then
212 : : call ddb_omega_interpol(Crystal%amu, ddb, ddb_lw, dtset%eta, dtset%prefix_outdata, &
213 : : & dtset%magpen, dtset%mpatpol, dtset%mpdir, dtset%mpert, dtset%mpopt, Crystal%natom, dtset%nfreq, Crystal%ntypat, &
214 : 6 : & dtset%freqflag, dtset%frmax, dtset%frmin, dtset%prtvol, Crystal%typat, Crystal%ucvol, Crystal%xred)
215 : : end if
216 : :
217 : : !Proceed with a normal anaddb run with relaxed- or fixed-spin quantities
218 [ - + ]: 6 : if (dtset%mpopt==1) then
219 [ # # # # : 0 : ddb%val= ddb%val_fs
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
220 [ + - ]: 6 : else if (dtset%mpopt==2) then
221 [ + - + - : 114284 : ddb%val= ddb%val_rs
+ - + - -
- - - - -
- - - - -
- - - + +
+ + + + ]
222 : : end if
223 : : end if
224 : :
225 : :
226 : : ! ! TODO: This is to maintain the previous behaviour in which all the arrays were initialized to zero.
227 : : ! ! In the new version asrq0%d2asr is always computed if the Gamma block is present
228 : : ! ! and this causes changes in [v5][t28]
229 : : ! if (.not. (dtset%ifcflag == 0 .or. dtset%instrflag /= 0 .or. dtset%elaflag /= 0)) then
230 : : ! asrq0%d2asr = zero
231 : : ! if (asrq0%asr == 3 .or. asrq0%asr == 4) then
232 : : ! asrq0%singular = zero; asrq0%uinvers = zero; asrq0%vtinvers = zero
233 : : ! end if
234 : : ! end if
235 : :
236 : : ! Acoustic Sum Rule
237 : : ! In case the interatomic forces are not calculated, the
238 : : ! ASR-correction (asrq0%d2asr) has to be determined here from the Dynamical matrix at Gamma.
239 : 162 : call asrq0%init(ddb, dtset%asr, dtset%rfmeth, crystal, dtset%sys_dim, driver%dcdq, driver%dcdqdq)
240 : :
241 : : ! Open netcdf output and write basic quantities
242 : 162 : call driver%open_write_nc(ana_ncid, dtset, crystal, comm)
243 : :
244 : : ! =========================================================================== !
245 : :
246 : : ! Compute dielectric tensor, Born effective charges, and quadrupoles.
247 [ + + ]: 162 : if (driver%do_electric_tensors) then
248 : 152 : call driver%electric_tensors(dtset, crystal, ddb, ddb_lw, ddb_hdr, ana_ncid, comm)
249 : : end if
250 : :
251 : : ! If low-dimensional systems, convert dielectric tensors if present
252 [ + + ]: 162 : if (dtset%sys_dim>1 .and. dtset%dipdip>0) then
253 : 4 : call driver%convertdim_dielt(crystal%rprimd, dtset%sys_dim,dtset%dielt_thick)
254 : : end if
255 : :
256 : : ! Structural response at fixed polarization
257 [ + + ]: 162 : if (dtset%polflag == 1) then
258 : 6 : call driver%structural_response(dtset, crystal, ddb)
259 : : end if
260 : :
261 : : ! Compute non-linear optical susceptibilities
262 : : ! and first-order change in the linear dielectric susceptibility
263 [ + + ]: 162 : if (dtset%nlflag > 0) then
264 : 8 : call driver%susceptibilities(dtset, ddb, ana_ncid, comm)
265 : : end if
266 : :
267 : : ! Interatomic force constants
268 [ + + ]: 162 : if (driver%do_ifc) then
269 : 108 : call driver%interatomic_force_constants(Ifc, dtset, crystal, ddb, ana_ncid, comm)
270 [ - + ]: 108 : if (dtset%flexoflag /= 1 .and. dtset%asr == 6) then
271 : 0 : write(msg, '(a, a)' )' IFCs derivatives computed from real-space moment',ch10
272 : 0 : call wrtout(units, msg)
273 : 0 : call asrq0%init(ddb, dtset%asr, dtset%rfmeth, crystal, dtset%sys_dim, driver%dcdq, driver%dcdqdq)
274 : : end if
275 : : end if
276 : :
277 : : ! Phonon density of states
278 [ + + ]: 162 : if (driver%do_phonon_dos) then
279 : 22 : call driver%phdos(dtset, crystal, Ifc, comm)
280 : : end if
281 : :
282 : : ! Phonon density of states and thermodynamical properties calculation
283 [ + + + + : 470 : if (dtset%ifcflag == 1 .and. any(dtset%thmflag==[1, 2])) then
+ + ]
284 : 8 : call driver%harmonic_thermo(dtset, crystal, Ifc, comm)
285 : : end if
286 : :
287 : : ! Phonon band structure
288 [ + + ]: 162 : if (driver%do_phonon_bs) then
289 : 130 : call mkphbs(Ifc, crystal, dtset, ddb, asrq0, dtset%prefix_outdata, comm)
290 : : end if
291 : :
292 : : ! DDB interpolation
293 [ + + ]: 162 : if (dtset%prtddb == 1 .and. dtset%ifcflag == 1) then
294 : 2 : call ddb_interpolate(Ifc, crystal, dtset, ddb, ddb_hdr, asrq0, comm)
295 : : end if
296 : :
297 : : ! =========================================================================== !
298 : : ! Electron-phonon section
299 : :
300 [ + + ]: 162 : if (dtset%elphflag == 1) then
301 : 30 : call elphon(dtset, crystal, Ifc, comm)
302 : : end if
303 : :
304 : : ! Thermal supercell calculation
305 [ + + + + : 2106 : if (sum(abs(dtset%thermal_supercell))>0 .and. dtset%ifcflag == 1) then
+ + ]
306 : 2 : call driver%thermal_supercell(dtset, crystal, ifc)
307 : : end if
308 : :
309 : : ! Thermal corrections to eigenvalues (old)
310 [ + + ]: 162 : if (dtset%thmflag >= 3 .and. dtset%thmflag <= 8) then
311 : 14 : call thmeig(dtset, ddb, crystal, ab_out, crystal%natom, dtset%mpert, dtset%msize, asrq0%d2asr, comm)
312 : : end if
313 : :
314 : : ! =========================================================================== !
315 : :
316 : : ! Compute the dielectric function and oscillator strength.
317 [ + + ]: 162 : if (driver%do_dielectric_q0) then
318 : 54 : call driver%dielectric_q0(dtset, crystal, ifc, ddb, asrq0, ana_ncid, comm)
319 : : end if
320 : :
321 : : ! Non-linear response: electrooptic and Raman (q = Gamma, TO modes only)
322 [ + + ]: 162 : if (dtset%nlflag == 1) then
323 : 4 : call driver%nonlinear_response(dtset, crystal, ana_ncid, comm)
324 : : end if
325 : :
326 : : ! Non-analyticity in the dynamical matrix
327 [ + + ]: 162 : if (driver%do_dielectric_nonana) then
328 : 48 : call driver%dielectric_nonana(dtset, crystal, ddb, ana_ncid, comm)
329 : : end if
330 : :
331 : : ! =========================================================================== !
332 : : ! Linear response with strain
333 : :
334 : : ! Internal strain (needed for the other linear response functions)
335 [ + + ]: 162 : if (dtset%instrflag /= 0) then
336 : 12 : call driver%internal_strain(dtset, ddb, asrq0)
337 : : end if
338 : :
339 : : ! Elastic tensor
340 [ + + ]: 162 : if (dtset%elaflag /= 0) then
341 : 12 : call driver%elastic_tensor(dtset, crystal, ddb, asrq0, ana_ncid)
342 : : end if
343 : :
344 : : ! Piezoelectric tensor
345 [ + + ]: 162 : if (dtset%piezoflag /= 0 .or. dtset%dieflag == 4 .or. dtset%elaflag == 4) then
346 : 10 : call driver%piezoelectric_tensor(dtset, crystal, ddb, ana_ncid)
347 : : end if
348 : :
349 : : ! Flexoelectric tensor
350 [ + + ]: 162 : if (dtset%flexoflag /= 0) then
351 : 4 : call driver%flexoelectric_tensor(dtset, crystal, ddb, ddb_lw, ddb_hdr, asrq0)
352 : : end if
353 : :
354 : : ! =========================================================================== !
355 : :
356 : : ! Gruneisen parameters
357 [ + + ]: 162 : if (dtset%gruns_nddbs /= 0) then
358 : 2 : call gruns_anaddb(dtset, comm)
359 : : end if
360 : :
361 : : ! Lattice Wannier functions
362 [ + + ]: 162 : if (dtset%ifcflag == 1 .and. dtset%lwfflag > 0 ) then
363 : 6 : call driver%lattice_wannier(dtset, crystal, Ifc, comm)
364 : : endif
365 : :
366 : : ! Output phonon frequencies for BoltzTrap
367 [ + + ]: 162 : if (iam_master .and. dtset%ifcflag == 1 .and. dtset%outboltztrap == 1) then
368 : 2 : call ifc%outphbtrap(crystal, dtset%ng2qpt, 1, dtset%q2shft, dtset%prefix_outdata)
369 : : end if
370 : :
371 : : ! =========================================================================== !
372 : : ! Close netcdf file
373 [ + - ]: 162 : if (iam_master) then
374 [ - + ]: 162 : NCF_CHECK(nf90_close(ana_ncid))
375 : : end if
376 : :
377 : : ! =========================================================================== !
378 : : ! Free memory
379 : 162 : call asrq0%free()
380 : 162 : call ifc%free()
381 : 162 : call crystal%free()
382 : 162 : call ddb%free()
383 : 162 : call ddb_hdr%free()
384 : 162 : call ddb_lw%free()
385 : 162 : call driver%free()
386 : 162 : call dtset%free()
387 : :
388 : : ! =========================================================================== !
389 : : ! Output timing and memory reports, then close output files
390 : :
391 : 162 : call timein(tcpu, twall)
392 : 162 : tsec(1)=tcpu-tcpui; tsec(2)=twall-twalli
393 : 162 : write(msg, '(a, i4, a, f13.1, a, f13.1)' )' Proc.',my_rank, ' individual time (sec): cpu=',tsec(1), ' wall=',tsec(2)
394 : 162 : call wrtout(std_out, msg)
395 : :
396 [ + - ]: 162 : if (iam_master) then
397 : 162 : write(ab_out, '(a, a, a, i4, a, f13.1, a, f13.1)' )'-',ch10, &
398 : 324 : '- Proc.',my_rank, ' individual time (sec): cpu=',tsec(1), ' wall=',tsec(2)
399 : : end if
400 : :
401 : 162 : call xmpi_sum(tsec, comm, ierr)
402 : :
403 : 162 : write(msg, '(a, (80a), a, a, a, f11.3, a, f11.3, a, a, a, a)' ) ch10, &
404 [ + + ]: 13284 : ('=',ii = 1, 80), ch10, ch10, &
405 : 162 : '+Total cpu time',tsec(1), ' and wall time',tsec(2), ' sec',ch10, ch10, &
406 [ - + ]: 13446 : ' anaddb : the run completed successfully.'
407 : 162 : call wrtout(units, msg)
408 : :
409 [ + - ]: 162 : if (iam_master) then
410 : : ! Write YAML document with the final summary.
411 : : ! we use this doc to test whether the calculation is completed.
412 : 162 : write(std_out, "(a)")"--- !FinalSummary"
413 : 162 : write(std_out, "(a)")"program: anaddb"
414 : 162 : write(std_out, "(2a)")"version: ",trim(abinit_version)
415 : 162 : write(std_out, "(2a)")"start_datetime: ",start_datetime
416 : 162 : write(std_out, "(2a)")"end_datetime: ",asctime()
417 : 162 : write(std_out, "(a, f13.1)")"overall_cpu_time: ",tsec(1)
418 : 162 : write(std_out, "(a, f13.1)")"overall_wall_time: ",tsec(2)
419 : 162 : write(std_out, "(a, i0)")"mpi_procs: ",xmpi_comm_size(xmpi_world)
420 : 162 : write(std_out, "(a, i0)")"omp_threads: ",xomp_get_num_threads(open_parallel=.True.)
421 : : !write(std_out, "(a, i0)")"num_warnings: ",nwarning
422 : : !write(std_out, "(a, i0)")"num_comments: ",ncomment
423 : 162 : write(std_out, "(a)")"..."
424 : 162 : call flush_unit(std_out)
425 : : end if
426 : :
427 : : ! Write information on file about the memory before ending mpi module, if memory profiling is enabled
428 : 162 : call abinit_doctor(dtset%filename_output)
429 : :
430 : 162 : call flush_unit(ab_out)
431 : 162 : call flush_unit(std_out)
432 : :
433 [ + - ]: 162 : if (iam_master) close(ab_out)
434 : :
435 : 162 : 100 call xmpi_end()
436 : :
437 : 0 : end program anaddb
438 : : !!***
439 : :
|
