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 :
|
