Branch data Line data Source code
1 : : !!****p* ABINIT/atdep
2 : : !! NAME
3 : : !! atdep
4 : : !!
5 : : !! FUNCTION
6 : : !! Calculations of phonons using molecular dynamic simulations.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (FB,JB,GA)
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 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : 92 : program atdep
24 : :
25 : 92 : use defs_basis
26 : : use m_abicore
27 : : use m_phonons
28 : : use m_errors
29 : : use m_abi_linalg
30 : : use m_xmpi
31 : : use m_abihist
32 : : use m_io_tools
33 : : use m_argparse
34 : :
35 : : use m_build_info, only : abinit_version
36 : : use m_time, only : asctime, timein, timab, cwtime, cwtime_report
37 : : use m_ifc, only : ifc_type
38 : : use m_crystal, only : crystal_t
39 : : use m_ddb, only : ddb_type
40 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type, tdep_read_input
41 : : use m_tdep_latt, only : Lattice_type, tdep_make_latt
42 : : use m_tdep_sym, only : Symmetries_type, tdep_SearchS_1at
43 : : use m_tdep_supercell, only : tdep_Supercell_type
44 : : use m_tdep_shell, only : Shell_type
45 : : use m_tdep_qpt, only : Qbz_type, Qpoints_type
46 : : use m_tdep_phi2, only : Eigen_type
47 : : use m_tdep_sampling, only : tdep_Sampling_type, tdep_match_supercell
48 : : use m_tdep_solver, only : tdep_Solver_type
49 : : use m_tdep_model, only : tdep_Model_type
50 : : use m_tdep_phdos, only : tdep_calc_phdos
51 : : use m_tdep_utils, only : tdep_print_header, tdep_print_Aknowledgments
52 : : use m_tdep_abitypes, only : tdep_init_crystal, tdep_init_ddb, tdep_write_ddb, &
53 : : & tdep_init_ifc, tdep_read_ifc, tdep_write_ifc, &
54 : : & tdep_ifc2phi2, tdep_test_ifc2phi2
55 : :
56 : : implicit none
57 : :
58 : : integer :: print_mem_report
59 : : integer :: iorder
60 : : real(dp) :: rotation(3,3)
61 : : real(dp) :: tcpu, tcpui, twall, twalli
62 : : real(dp) :: tsec(2)
63 : : character(len = 24):: start_datetime
64 : : integer, parameter :: ntimers = 11
65 : : integer :: itimer
66 : : real(dp) :: cpu, wall, gflops
67 : : character(len=64) :: timer_label(ntimers)
68 : : real(dp) :: timer_wall(ntimers)
69 : : type(args_t) :: args
70 : 92 : type(atdep_dataset_type) :: Invar
71 : 92 : type(MPI_enreg_type) :: MPIdata
72 : 92 : type(abihist) :: Hist
73 [ + + + + ]: 552 : type(ifc_type) :: Ifc
74 : 92 : type(ddb_type) :: DDB
75 [ + + + + : 4784 : type(crystal_t) :: Crystal
+ + + + +
+ + + + +
+ + + + ]
76 : : type(Lattice_type) :: Lattice
77 : 92 : type(tdep_Supercell_type) :: SC
78 : 92 : type(Symmetries_type) :: Sym
79 : 92 : type(Qpoints_type) :: Qpt
80 : 92 : type(Qbz_type) :: Qbz
81 : 92 : type(Shell_type) :: Shell1at, Shell2at, Shell3at, Shell4at
82 : 92 : type(tdep_Sampling_type) :: MD
83 : 92 : type(tdep_Solver_type) :: Solver
84 : 92 : type(tdep_Model_type) :: Model
85 : 92 : type(Eigen_type) :: Eigen2nd_MP, Eigen2nd_path
86 [ + + + + ]: 1196 : type(phdos_t) :: PHdos
87 : :
88 : : !******************************************************************
89 : :
90 : : !==========================================================================================
91 : : !===================== Initialization & Reading ==========================================
92 : : !==========================================================================================
93 : : ! Change communicator for I/O (mandatory!)
94 : 92 : call abi_io_redirect(new_io_comm=xmpi_world)
95 : : ! Initialize MPI
96 : 92 : call xmpi_init()
97 : :
98 : : ! Initialisation of the timing
99 : 92 : call timein(tcpui, twalli)
100 : 92 : start_datetime = asctime()
101 : 92 : call timab(1, 0, tsec)
102 : :
103 : : ! Initialisation of the timing breakdown (reported to std_out, see record_timer below).
104 : 92 : itimer = 0
105 : 92 : call cwtime(cpu, wall, gflops, "start")
106 : :
107 : : ! Parse command line arguments.
108 [ - + ]: 92 : args = args_parser(); if (args%exit /= 0) goto 100
109 : :
110 : : ! Initialize memory profiling if activated at configure time.
111 : : ! if a full report is desired, set the argument of abimem_init to "2" instead of "0" via the command line.
112 : : ! note that the file can easily be multiple GB in size so don't use this option normally
113 : : #ifdef HAVE_MEM_PROFILING
114 : : call abimem_init(args%abimem_level, limit_mb=args%abimem_limit_mb)
115 : : #endif
116 : :
117 : : ! Read main input file and hist file and initalize Invar
118 : 92 : call tdep_read_input(args%input_path,Hist,Invar)
119 : :
120 : : ! Print header
121 : 92 : call tdep_print_header(ab_out)
122 : :
123 : : ! Echo the input variables to the main output file.
124 : 92 : call Invar%print(ab_out)
125 : :
126 : 92 : call MPIdata%init(Invar)
127 : 92 : call MD%init(Invar,MPIdata,Hist)
128 : 92 : call abihist_free(Hist)
129 : :
130 [ - + ]: 92 : if (args%dry_run /= 0) then
131 : 0 : call wrtout(std_out, "Dry run mode. Exiting after have read the input")
132 : 0 : call MD%free()
133 : 0 : call Invar%free()
134 : 0 : call MPIdata%free()
135 : 0 : goto 100
136 : : end if
137 : :
138 : : ! Initialize basic quantities
139 : 92 : print_mem_report = 1
140 : :
141 : 92 : call record_timer("Initialization and reading input")
142 : :
143 : : !==========================================================================================
144 : : !============== Define the ideal lattice, symmetries and Brillouin zone ===================
145 : : !==========================================================================================
146 : :
147 : : !Define all the quantities needed to buid the lattice (rprim*, acell*, brav*...)
148 : 92 : call tdep_make_latt(Invar,Lattice,rotation)
149 : :
150 : : !Compute all the symmetries coming from the bravais lattice
151 : 92 : call Sym%init(Invar,MPIdata,Lattice)
152 : :
153 : : !Build the Crystal object describing the unit cell.
154 : 92 : call tdep_init_crystal(Crystal,Invar,Lattice,Sym)
155 : :
156 : 92 : call SC%init(Crystal,Invar%multiplicity,Lattice%acell_unitcell)
157 : 92 : call SC%print(std_out)
158 : :
159 : : !Initialize the Brillouin zone and compute the q-points path
160 : 92 : call Qpt%init(Invar,MPIdata,Lattice)
161 : :
162 : 92 : call record_timer("Lattice, symmetry and Brillouin zone setup")
163 : :
164 : : !==========================================================================================
165 : : !======== Complete the initialization of the Sampling and Symmetry objects ================
166 : : !==========================================================================================
167 : :
168 : : !Shift xred to keep atoms in the same unit cell at each step.
169 : 92 : call MD%gather_inbox(MPIdata)
170 : :
171 : : !Apply rotation to cartesian forces
172 : 92 : call MD%rotate(rotation)
173 : :
174 : : !Map the atoms of the input moledular dynamics to the ideal supercell
175 : 92 : call tdep_match_supercell(Invar,MPIdata,MD,SC)
176 : :
177 : : !Find the symmetry operations between atoms of the supercell
178 : 92 : call tdep_SearchS_1at(Invar,MPIdata,Sym,SC)
179 : :
180 : 92 : call record_timer("Supercell matching and symmetry search")
181 : :
182 : : !==========================================================================================
183 : : !============== Initialize all Shell datatypes ============================================
184 : : !==========================================================================================
185 : :
186 : : !Shells for 1st order coefficients
187 : 92 : call Shell1at%init(Invar,MPIdata,SC,Sym,1)
188 : :
189 : : !Shells for 2nd order coefficients
190 : 92 : call Shell2at%init(Invar,MPIdata,SC,Sym,2)
191 : :
192 : : !Shells for 3rd order coefficients
193 [ + + ]: 92 : if (Invar%order.ge.3) call Shell3at%init(Invar,MPIdata,SC,Sym,3)
194 : :
195 : : !Shells for 4th order coefficients
196 [ + + ]: 92 : if (Invar%order.eq.4) call Shell4at%init(Invar,MPIdata,SC,Sym,4)
197 : :
198 : 92 : call record_timer("Shell initialization")
199 : :
200 : : !==========================================================================================
201 : : !======================== Initialize main TDEP datatypes ==================================
202 : : !==========================================================================================
203 : :
204 : 92 : call Model%init(Invar, Lattice, Shell3at, Shell4at)
205 : 92 : call Solver%init(Invar, Shell1at, Shell2at, Shell3at, Shell4at)
206 : :
207 : 92 : call tdep_init_ddb(Invar,MPIdata,Crystal,DDB,Lattice,Qbz)
208 : 92 : call tdep_init_ifc(Invar,MPIdata,Crystal,DDB,Ifc,Lattice,Model%Phi2,MD%Rlatt_scaled,Shell2at,Sym)
209 : :
210 : 92 : call record_timer("Model, Solver, Crystal, DDB and IFC initialization")
211 : :
212 : : !==========================================================================================
213 : : !================= Copy the forces into the solver ========================================
214 : : !==========================================================================================
215 : : !Remove the supercell contribution (the "LR part") included in total forces
216 : : !before computing the "SR part". The full "LR part" will be added later.
217 [ + + ]: 92 : if (Invar%dipdip) then
218 : 10 : call tdep_ifc2phi2(Ifc%dipdip,Ifc,Invar,Lattice,Invar%natom_unitcell,1,Model%Phi2,MD%Rlatt_scaled,Shell2at,Sym)
219 : 10 : call Model%calc_ftot2(Invar,Model%Phi2%Tot,MD%ucart)
220 : : end if
221 : :
222 : 92 : call Solver%set_residual_forces(MD, Model)
223 : :
224 : : !==========================================================================================
225 : : !================= Build fcoeff and compute constraints ===================================
226 : : !==========================================================================================
227 : :
228 [ + + ]: 92 : if (Invar%readifc.ne.1) then
229 : 90 : call Solver%calc_phi1fcoeff(Shell1at,Sym)
230 : 90 : call Solver%calc_phi2fcoeff(Shell2at,Sym,MD)
231 : : end if
232 : :
233 [ + + ]: 92 : if (Invar%order.ge.3) then
234 : 32 : call Solver%calc_phi3fcoeff(Shell3at,Sym,MD)
235 : : end if
236 : :
237 [ + + ]: 92 : if (Invar%order.eq.4) then
238 : 18 : call Solver%calc_phi4fcoeff(Shell4at,Sym,MD)
239 : : end if
240 : :
241 : : call Solver%calc_constraints(Invar,MPIdata,SC,Sym,&
242 : 92 : & Shell1at,Shell2at,Shell3at,Shell4at)
243 : :
244 : :
245 : : !==========================================================================================
246 : : !============= Compute the pseudo inverse using the Moore-Penrose method ==================
247 : : !==========================================================================================
248 : :
249 : : !----------------- If all the Orders are solved simultaneously
250 [ + + ]: 92 : if (Invar%together.eq.1) then
251 : 90 : write(ab_out,*) '############### (Solve simultaneously all the orders) #######################'
252 : :
253 [ + + ]: 90 : if (Invar%readifc.ne.1) then
254 : 88 : call Solver%calc_MoorePenrose(MPIdata,0)
255 : 88 : call Solver%calc_phi1(Shell1at,Sym,Model%Phi1)
256 : 88 : call Solver%calc_phi2(Shell2at,Sym,Model%Phi2%SR)
257 : : end if
258 : 90 : call Model%calc_ftot2(Invar,Model%Phi2%SR,MD%ucart)
259 : :
260 [ + + ]: 90 : if (Invar%order.ge.3) then
261 : 30 : call Solver%calc_phi3ref(Shell3at,Model%Phi3)
262 : 30 : call Model%calc_ftot3(Invar,Shell3at,MD%ucart,Sym)
263 : : end if
264 : :
265 [ + + ]: 90 : if (Invar%order.ge.4) then
266 : 16 : call Solver%calc_phi4ref(Shell4at,Model%Phi4)
267 : 16 : call Model%calc_ftot4(Invar,Shell4at,MD%ucart,Sym)
268 : : end if
269 : :
270 : : !----------------- If all the Orders are solved successively
271 : : !
272 : : !ATTENTION : Le LR semble enleve a l'ordre 2 mais pas a l'ordre 3 et 4.
273 : : ! On repart de MD%Forces tout en bas et pas de Solver%Forces
274 : : ! (car les ordres 2 et 3 sont supprimes)
275 : :
276 [ + - ]: 2 : else if (Invar%together.eq.0) then
277 : 2 : write(ab_out,*) '################## (Solve successively each order) ##########################'
278 : :
279 [ + + ]: 8 : do iorder=1,Invar%order-1
280 : 6 : write(ab_out,*) ' For order=',iorder+1
281 [ - + ]: 6 : if (Invar%readifc.eq.1) cycle
282 : :
283 : 6 : call Solver%calc_MoorePenrose(MPIdata,iorder)
284 : :
285 [ + + ]: 6 : if (iorder.eq.1) then
286 [ + - ]: 2 : if (Invar%readifc.ne.1) then
287 : 2 : call Solver%calc_phi1(Shell1at,Sym,Model%Phi1)
288 : 2 : call Solver%calc_phi2(Shell2at,Sym,Model%Phi2%SR)
289 : : end if
290 : 2 : call Model%calc_ftot2(Invar,Model%Phi2%SR,MD%ucart)
291 : :
292 [ + + ]: 4 : else if (iorder.eq.2) then
293 : 2 : call Solver%calc_phi3ref(Shell3at,Model%Phi3)
294 : 2 : call Model%calc_ftot3(Invar,Shell3at,MD%ucart,Sym)
295 : :
296 [ + - ]: 2 : else if (iorder.eq.3) then
297 : 2 : call Solver%calc_phi4ref(Shell4at,Model%Phi4)
298 : 2 : call Model%calc_ftot4(Invar,Shell4at,MD%ucart,Sym)
299 : : end if
300 : :
301 : : ! Remove forces computed with the previous order
302 : 8 : call Solver%set_residual_forces(MD, Model)
303 : :
304 : : end do ! iorder
305 : : end if
306 : :
307 : : ! Add the long-range part of the IFC
308 [ + + ]: 92 : if (Invar%dipdip) then
309 [ + - + - : 1847444 : Model%Phi2%Tot = Model%Phi2%LR + Model%Phi2%SR
+ - - - +
+ + + ]
310 : : end if
311 : :
312 : : ! Free some memory
313 : 92 : call Solver%free()
314 : 92 : call Shell1at%free()
315 : :
316 : 92 : call record_timer("Coefficient fitting and constraints")
317 : :
318 : : !==========================================================================================
319 : : !=================== Write the IFC and check the constraints ==============================
320 : : !==========================================================================================
321 : 92 : call Model%write_phi1()
322 : 92 : call Model%write_phi2(Invar,MPIdata,SC,Shell2at)
323 : :
324 [ + + ]: 92 : if (Invar%order.ge.3) call Model%write_phi3(SC,Shell3at,Sym)
325 [ + + ]: 92 : if (Invar%order.ge.4) call Model%write_phi4(SC,Shell4at,Sym)
326 : :
327 : 92 : call Model%check_constraints(Invar,SC,Sym,Shell3at,Shell4at)
328 : :
329 : : !==========================================================================================
330 : : !================= Convert Phi2 into IFC and DDB objects ==================================
331 : : !==========================================================================================
332 : :
333 : : call tdep_ifc2phi2(Ifc%dipdip,Ifc,Invar,Lattice,Invar%natom_unitcell,0,&
334 : 92 : & Model%Phi2,MD%Rlatt_scaled,Shell2at,Sym)
335 : :
336 : : ! Write Ifc%atmfrc in the ifc_out.dat file
337 [ + - ]: 92 : if (MPIdata%iam_master) then
338 : 92 : call tdep_write_ifc(Crystal,Ifc,Invar,Invar%natom_unitcell,0)
339 : : end if
340 : :
341 : : ! For test purpose: read the previous IFC from ifc_out.dat and write it in ifc_check.dat
342 [ + + ]: 92 : if (Invar%readifc.eq.2) then
343 : : call tdep_test_ifc2phi2(Invar,MPIdata,Lattice,Crystal,Sym,Shell2at,Model%Phi2,DDB,&
344 : 4 : & MD%Rlatt_scaled,Invar%natom,Invar%natom_unitcell)
345 : : end if
346 : :
347 : : ! Create a new DDB with the coarse q-point grid in the IBZ.
348 : 92 : call DDB%free()
349 : 92 : call Ifc%to_ddb(DDB,Crystal)
350 : 92 : call tdep_write_ddb(DDB,Crystal,Invar)
351 : :
352 : 92 : call record_timer("Writing force constants to IFC and DDB")
353 : :
354 : : !==========================================================================================
355 : : !===================== Compute the phonon spectrum, the DOS, ==============================
356 : : !===================== the dynamical matrix and write them ===============================
357 : : !==========================================================================================
358 : 92 : call Eigen2nd_MP%init(Invar%natom_unitcell,Qbz%nqbz)
359 : 92 : call Eigen2nd_path%init(Invar%natom_unitcell,Qpt%nqpt)
360 : :
361 : 92 : call tdep_calc_phdos(Invar,MPIdata,Crystal,Ifc,Qbz,Qpt,Eigen2nd_MP,Eigen2nd_path,PHdos)
362 : 92 : call Shell2at%free()
363 : :
364 : 92 : call record_timer("Phonon spectrum and DOS")
365 : :
366 : : !==========================================================================================
367 : : !================= Elastic constants, free energy and thermodynamics ======================
368 : : !==============================================e===========================================
369 : :
370 : : !Compute the elastic constantsD======================================
371 : 92 : call Model%calc_elastic(SC,Invar)
372 : :
373 : : !Compute U_0, the "free energy" from the model
374 : 92 : call Model%compute_energy(Invar,MPIdata,MD)
375 : :
376 : : !Compute thermodynamical quantities
377 : 92 : call Model%calc_thermo(Invar,MPIdata,PHdos)
378 : :
379 : : !Free some memory
380 : 92 : call PHdos%free()
381 : :
382 : 92 : call record_timer("Elastic constants, free energy and thermodynamics")
383 : :
384 : : !==========================================================================================
385 : : !===================== Calculation of 3rd order quantities ================================
386 : : !==========================================================================================
387 [ + + ]: 92 : if (Invar%order>2) then
388 : 32 : call Model%calc_gruneisen(Invar,MPIdata,SC,Eigen2nd_path,Qpt,MD%Rlatt_cart,Shell3at,Sym)
389 : 32 : call Model%calc_alpha_gamma(Invar,MPIdata,SC,Eigen2nd_MP,Qbz,MD%Rlatt_cart,Shell3at,Sym)
390 : : !FB call tdep_calc_lifetime1(Crystal,SC%dij,Eigen2nd_MP,Ifc,Invar,Lattice,Model%Phi3,Qbz,MD%Rlatt_cart,Shell3at,Sym)
391 : : end if
392 : :
393 : 92 : call record_timer("Third-order (Gruneisen) properties")
394 : :
395 : : !==========================================================================================
396 : : !===================== Free memory ========================================================
397 : : !==========================================================================================
398 : :
399 [ + + ]: 92 : if (Invar%order.ge.3) call Shell3at%free()
400 [ + + ]: 92 : if (Invar%order.eq.4) call Shell4at%free()
401 : :
402 : 92 : call Invar%free()
403 : 92 : call MPIdata%free()
404 : 92 : call SC%free()
405 : 92 : call MD%free()
406 : 92 : call Model%free()
407 : 92 : call Ifc%free()
408 : 92 : call DDB%free()
409 : 92 : call Crystal%free()
410 : 92 : call Sym%free()
411 : :
412 : 92 : call Qbz%free()
413 : 92 : call Qpt%free()
414 : 92 : call Eigen2nd_MP%free()
415 : 92 : call Eigen2nd_path%free()
416 : :
417 : 92 : call tdep_print_Aknowledgments(ab_out)
418 : 92 : call flush_unit(ab_out)
419 : 92 : call record_timer("Cleanup and finalization")
420 : :
421 : : !==========================================================================================
422 : : !===================== End the calculation ================================================
423 : : !==========================================================================================
424 : :
425 : 92 : call timein(tcpu, twall)
426 : 92 : tsec(1)=tcpu-tcpui; tsec(2)=twall-twalli
427 : :
428 : : ! Write YAML document with the final summary.
429 [ + - ]: 92 : if (MPIdata%iam_master) then
430 : :
431 : : ! Write YAML document with the timing breakdown.
432 : 92 : write(std_out, "(a)")""
433 : 92 : write(std_out, "(a)")"--- !Timing"
434 : 92 : write(std_out, "(a, f13.1)")"total_wall_time: ", tsec(2)
435 : 92 : write(std_out, "(a)")"breakdown:"
436 [ + + ]: 1104 : do itimer=1,ntimers
437 : : write(std_out, "(a, a, a, f13.1, a)") &
438 [ + - ]: 1104 : & " - {label: ", '"'//trim(timer_label(itimer))//'"', ", wall: ", timer_wall(itimer), "}"
439 : : end do
440 : 92 : write(std_out, "(a)")"..."
441 : :
442 : :
443 : 92 : write(std_out, "(a)")""
444 : 92 : write(std_out, "(a)")"--- !FinalSummary"
445 : 92 : write(std_out, "(a)")"program: atdep"
446 : 92 : write(std_out, "(2a)")"version: ", trim(abinit_version)
447 : 92 : write(std_out, "(2a)")"start_datetime: ", start_datetime
448 : 92 : write(std_out, "(2a)")"end_datetime: ", asctime()
449 : 92 : write(std_out, "(a, f13.1)")"overall_cpu_time: ", tsec(1)
450 : 92 : write(std_out, "(a, f13.1)")"overall_wall_time: ", tsec(2)
451 : 92 : write(std_out, "(a, i0)")"mpi_procs: ", MPIdata%nproc
452 : 92 : write(std_out, "(a)")"..."
453 : 92 : call flush_unit(std_out)
454 : : end if
455 : :
456 : : !Memory analysis
457 [ + - ]: 92 : call abinit_doctor(trim(Invar%output_prefix), print_mem_report=print_mem_report)
458 : 92 : call flush_unit(std_out)
459 : 92 : close(unit=ab_out)
460 : 92 : 100 call xmpi_end()
461 : :
462 : : contains
463 : :
464 : : ! Close the timer opened by the previous record_timer call (or by the initial cwtime "start" call).
465 : 2024 : subroutine record_timer(label)
466 : :
467 : : character(len=*),intent(in) :: label
468 : :
469 : 1012 : itimer = itimer + 1
470 : 1012 : timer_label(itimer) = label
471 : 1012 : call cwtime_report(label, cpu, wall, gflops, out_wall=timer_wall(itimer))
472 : :
473 [ + - ]: 1104 : end subroutine record_timer
474 : :
475 : : end program atdep
476 : : !!***
|
