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