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 : 88 : program atdep
24 : :
25 : 88 : 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_time, only : asctime, timein, timab
36 : : use m_ifc, only : ifc_type
37 : : use m_crystal, only : crystal_t
38 : : use m_ddb, only : ddb_type
39 : : use m_tdep_abitypes, only : Qbz_type, tdep_init_crystal, tdep_init_ifc, tdep_init_ddb, tdep_write_ddb, &
40 : : & tdep_destroy_qbz, tdep_ifc2phi2, tdep_read_ifc, tdep_write_ifc
41 : : use m_tdep_latt, only : Lattice_type, tdep_make_latt
42 : : use m_tdep_sym, only : tdep_make_sym, Symmetries_type, tdep_destroy_sym
43 : : use m_tdep_dataset, only : tdep_read_input, tdep_init_MPIdata, &
44 : : & tdep_destroy_mpidata, atdep_dataset_type, MPI_enreg_type, tdep_destroy_invar, version_string
45 : : use m_tdep_qpt, only : tdep_make_qptpath, Qpoints_type, tdep_destroy_qpt
46 : : use m_tdep_sampling, only : tdep_Sampling_type, tdep_sampling_init_read, tdep_sampling_free,&
47 : : & tdep_sampling_shift_xred, tdep_sampling_rotate, tdep_MatchIdeal2Average
48 : : use m_tdep_shell, only : Shell_type, tdep_init_shell2at, tdep_init_shell3at, tdep_init_shell4at, &
49 : : & tdep_init_shell1at, tdep_destroy_shell
50 : : use m_tdep_solver, only : tdep_Solver_type, tdep_solver_init, tdep_solver_free, tdep_solver_set_residual_forces, &
51 : : & tdep_calc_phi1fcoeff, tdep_calc_phi2fcoeff, tdep_calc_phi3fcoeff, tdep_calc_phi4fcoeff, &
52 : : & tdep_calc_MoorePenrose, tdep_calc_constraints
53 : : use m_tdep_model, only : tdep_Model_type, tdep_model_init, tdep_model_free, tdep_calc_model,&
54 : : & Phi2_type, tdep_init_phi2, tdep_destroy_phi2
55 : : use m_tdep_phi2, only : tdep_calc_phi2, tdep_write_phi2, tdep_calc_ftot2, &
56 : : & Eigen_type, tdep_init_eigen2nd, tdep_destroy_eigen2nd, tdep_calc_phi1, tdep_write_phi1
57 : : use m_tdep_phi3, only : tdep_calc_phi3ref, tdep_write_phi3, tdep_calc_ftot3, &
58 : : & tdep_calc_alpha_gamma, tdep_write_gruneisen
59 : : use m_tdep_phi4, only : tdep_calc_phi4ref, tdep_write_phi4, tdep_calc_ftot4
60 : : use m_tdep_phdos, only : tdep_calc_phdos,tdep_calc_elastic,tdep_calc_thermo
61 : : use m_tdep_utils, only : tdep_check_constraints, tdep_print_Aknowledgments
62 : :
63 : : implicit none
64 : :
65 : : integer :: print_mem_report
66 : : integer :: stdout,stdlog
67 : : integer :: iorder
68 : : real(dp) :: rotation(3,3)
69 : : real(dp) :: tcpu, tcpui, twall, twalli
70 : : real(dp) :: tsec(2)
71 : : character(len = 24):: start_datetime
72 : : type(args_t) :: args
73 : 88 : type(atdep_dataset_type) :: Invar
74 : 88 : type(MPI_enreg_type) :: MPIdata
75 : 88 : type(abihist) :: Hist
76 [ + + + + ]: 528 : type(ifc_type) :: Ifc
77 : 88 : type(ddb_type) :: DDB
78 [ + + + + : 4576 : type(crystal_t) :: Crystal
+ + + + +
+ + + + +
+ + + + ]
79 : : type(Lattice_type) :: Lattice
80 : 88 : type(Symmetries_type) :: Sym
81 : 88 : type(Qpoints_type) :: Qpt
82 : 88 : type(Qbz_type) :: Qbz
83 : 88 : type(Shell_type) :: Shell1at, Shell2at, Shell3at, Shell4at
84 : 88 : type(tdep_Sampling_type) :: MD
85 : 88 : type(tdep_Solver_type) :: Solver
86 : 88 : type(tdep_Model_type) :: Model
87 : 88 : type(Eigen_type) :: Eigen2nd_MP, Eigen2nd_path
88 [ + + + + ]: 1144 : type(phdos_t) :: PHdos
89 : :
90 : : !******************************************************************
91 : :
92 : : !==========================================================================================
93 : : !===================== Initialization & Reading ==========================================
94 : : !==========================================================================================
95 : : ! Change communicator for I/O (mandatory!)
96 : 88 : call abi_io_redirect(new_io_comm=xmpi_world)
97 : : ! Initialize MPI
98 : 88 : call xmpi_init()
99 : :
100 : : ! Initialisation of the timing
101 : 88 : call timein(tcpui, twalli)
102 : 88 : start_datetime = asctime()
103 : 88 : call timab(1, 0, tsec)
104 : :
105 : : ! Parse command line arguments.
106 [ - + ]: 88 : args = args_parser(); if (args%exit /= 0) goto 100
107 : :
108 : : ! Initialize memory profiling if activated at configure time.
109 : : ! if a full report is desired, set the argument of abimem_init to "2" instead of "0" via the command line.
110 : : ! note that the file can easily be multiple GB in size so don't use this option normally
111 : : #ifdef HAVE_MEM_PROFILING
112 : : call abimem_init(args%abimem_level, limit_mb=args%abimem_limit_mb)
113 : : #endif
114 : :
115 : : ! Read input values from the input.in input file
116 : 88 : call tdep_read_input(args%input_path,Hist,Invar)
117 : 88 : call tdep_init_MPIdata(Invar,MPIdata)
118 : 88 : call tdep_sampling_init_read(MD,Invar,MPIdata,Hist)
119 : 88 : call abihist_free(Hist)
120 : :
121 [ - + ]: 88 : if (args%dry_run /= 0) then
122 : 0 : call wrtout(std_out, "Dry run mode. Exiting after have read the input")
123 : 0 : call tdep_sampling_free(MD)
124 : 0 : call tdep_destroy_invar(Invar)
125 : 0 : call tdep_destroy_mpidata(MPIdata)
126 : 0 : goto 100
127 : : end if
128 : :
129 : : ! Initialize basic quantities
130 : 88 : print_mem_report = 1
131 : 88 : stdout = Invar%stdout
132 : 88 : stdlog = Invar%stdlog
133 : :
134 : : !==========================================================================================
135 : : !============== Define the ideal lattice, symmetries and Brillouin zone ===================
136 : : !==========================================================================================
137 : :
138 : : !Define all the quantities needed to buid the lattice (rprim*, acell*, brav*...)
139 : 88 : call tdep_make_latt(Invar,Lattice,rotation)
140 : :
141 : : !Compute all the symmetries coming from the bravais lattice
142 : 88 : call tdep_make_sym(Invar,Lattice,MPIdata,Sym)
143 : :
144 : : !Initialize the Brillouin zone and compute the q-points path
145 : 88 : call tdep_make_qptpath(Invar,Lattice,MPIdata,Qpt)
146 : :
147 : : !==========================================================================================
148 : : !============== Complete the initialization of the Sampling ===============================
149 : : !==========================================================================================
150 : :
151 : : !Shift xred to keep atoms in the same unit cell at each step.
152 : 88 : call tdep_sampling_shift_xred(MD, MPIdata)
153 : :
154 : : !Apply rotation to cartesian forces
155 : 88 : call tdep_sampling_rotate(MD, rotation)
156 : :
157 : : !Map the atoms of the input moledular dynamics to the ideal supercell
158 : 88 : call tdep_MatchIdeal2Average(MD,Invar,Lattice,Sym,MPIdata)
159 : :
160 : : !==========================================================================================
161 : : !============== Initialize the Shell1at datatype ==========================================
162 : : !==========================================================================================
163 : 88 : call tdep_init_shell1at(Shell1at,Invar,MD,Sym,MPIdata)
164 : :
165 : : !==========================================================================================
166 : : !============== Initialize the Shell2at datatype ==========================================
167 : : !==========================================================================================
168 : 88 : call tdep_init_shell2at(Shell2at,Invar,MD,Sym,MPIdata)
169 : :
170 : : !==========================================================================================
171 : : !============== Initialize the Shell3at datatype ==========================================
172 : : !==========================================================================================
173 [ + + ]: 88 : if (Invar%order.ge.3) then
174 : 28 : call tdep_init_shell3at(Shell3at,Invar,MD,Sym,MPIdata)
175 : : end if
176 : :
177 : : !==========================================================================================
178 : : !============== Initialize the Shell4at datatype ==========================================
179 : : !==========================================================================================
180 [ + + ]: 88 : if (Invar%order==4) then
181 : 16 : call tdep_init_shell4at(Shell4at,Invar,MD,Sym,MPIdata)
182 : : end if
183 : :
184 : : !==========================================================================================
185 : : !============== Initialize the TDEP Model datatype ========================================
186 : : !==========================================================================================
187 : 88 : call tdep_model_init(Model, Invar, Shell3at, Shell4at)
188 : :
189 : : !==========================================================================================
190 : : !============== Initialize the TDEP Solver datatype =======================================
191 : : !==========================================================================================
192 : 88 : call tdep_solver_init(Solver, Invar, Shell1at, Shell2at, Shell3at, Shell4at)
193 : :
194 : : !==========================================================================================
195 : : !============== Initialize Crystal, DDB, and IFC ABINIT Datatypes =========================
196 : : !==========================================================================================
197 : 88 : call tdep_init_crystal(Crystal,Invar,Lattice,Sym)
198 : 88 : call tdep_init_ddb(Crystal,DDB,Invar,Lattice,MPIdata,Qbz)
199 : 88 : call tdep_init_ifc(Crystal,DDB,Ifc,Invar,Lattice,MPIdata,Model%Phi2,MD%Rlatt_scaled,Shell2at,Sym)
200 : :
201 : : !==========================================================================================
202 : : !================= Copy the forces into the solver ========================================
203 : : !==========================================================================================
204 : : !Remove the supercell contribution (the "LR part") included in total forces
205 : : !before computing the "SR part". The full "LR part" will be added later.
206 : : ! GA: FIXME
207 : : ! Since Ifc%ewald_atfrc has not been initialized, Model%Forces are zero!
208 [ + + ]: 88 : if (Invar%loto) then
209 : 10 : call tdep_ifc2phi2(Ifc%dipdip,Ifc,Invar,Lattice,Invar%natom_unitcell,1,Model%Phi2,MD%Rlatt_scaled,Shell2at,Sym)
210 : 10 : call tdep_calc_ftot2(Model,Invar,Model%Phi2%Tot,MD%ucart)
211 : : end if
212 : :
213 : 88 : call tdep_solver_set_residual_forces(Solver, MD, Model)
214 : :
215 : : ! GA: I dont think this is even necessary
216 : : !if (Invar%loto) then
217 : : ! Model%Phi1(:) = zero
218 : : ! Model%Phi1Ui(:) = zero
219 : : ! Model%Phi2%SR(:,:) = zero
220 : : ! Model%Phi2%LR(:,:) = zero
221 : : ! Model%Phi2%Tot(:,:) = zero
222 : : ! Model%Phi2UiUj(:) = zero
223 : : ! Model%Forces(:) = zero
224 : : !end if
225 : :
226 : : !==========================================================================================
227 : : !================= Build fcoeff and compute constraints ===================================
228 : : !==========================================================================================
229 : :
230 [ + + ]: 88 : if (Invar%readifc.ne.1) then
231 : 86 : call tdep_calc_phi1fcoeff(Solver,Invar,Shell1at,Sym)
232 : 86 : call tdep_calc_phi2fcoeff(Solver,Invar,Shell2at,Sym,MD)
233 : : end if
234 : :
235 [ + + ]: 88 : if (Invar%order.ge.3) then
236 : 28 : call tdep_calc_phi3fcoeff(Solver,Invar,Shell3at,Sym,MD)
237 : : end if
238 : :
239 [ + + ]: 88 : if (Invar%order.eq.4) then
240 : 16 : call tdep_calc_phi4fcoeff(Solver,Invar,Shell4at,Sym,MD)
241 : : end if
242 : :
243 : : call tdep_calc_constraints(Solver,MD%distance,Invar,MPIdata,Sym,&
244 : 88 : & Shell1at,Shell2at,Shell3at,Shell4at)
245 : :
246 : :
247 : : !==========================================================================================
248 : : !============= Compute the pseudo inverse using the Moore-Penrose method ==================
249 : : !==========================================================================================
250 : :
251 : : !=================== If all the Orders are solved simultaneously ==========================
252 : :
253 [ + + ]: 88 : if (Invar%together.eq.1) then
254 : 86 : write(stdout,*) '############### (Solve simultaneously all the orders) #######################'
255 : :
256 [ + + ]: 86 : if (Invar%readifc.ne.1) then
257 : 84 : call tdep_calc_MoorePenrose(Solver,0,Invar,MPIdata)
258 : 84 : call tdep_calc_phi1(Solver,Shell1at,Sym,Model%Phi1)
259 : 84 : call tdep_calc_phi2(Solver,Shell2at,Sym,Model%Phi2%SR)
260 : : end if
261 : 86 : call tdep_calc_ftot2(Model,Invar,Model%Phi2%SR,MD%ucart)
262 : :
263 [ + + ]: 86 : if (Invar%order.ge.3) then
264 : 26 : call tdep_calc_phi3ref(Solver,Shell3at,Model%Phi3)
265 : 26 : call tdep_calc_ftot3(Model,Invar,Shell3at,MD%ucart,Sym)
266 : : end if
267 : :
268 [ + + ]: 86 : if (Invar%order.ge.4) then
269 : 14 : call tdep_calc_phi4ref(Solver,Shell4at,Model%Phi4)
270 : 14 : call tdep_calc_ftot4(Model,Invar,Shell4at,MD%ucart,Sym)
271 : : end if
272 : :
273 : : !=================== If all the Orders are solved successively ============================
274 : :
275 : : !ATTENTION : Le LR semble enleve a l'ordre 2 mais pas a l'ordre 3 et 4.
276 : : ! On repart de MD%Forces tout en bas et pas de Solver%Forces
277 : : ! (car les ordres 2 et 3 sont supprimes)
278 [ + - ]: 2 : else if (Invar%together.eq.0) then
279 : 2 : write(stdout,*) '################## (Solve successively each order) ##########################'
280 [ + + ]: 8 : do iorder=1,Invar%order-1
281 : 6 : write(stdout,*) ' For order=',iorder+1
282 [ - + ]: 6 : if (Invar%readifc.eq.1) cycle
283 : :
284 : 6 : call tdep_calc_MoorePenrose(Solver,iorder,Invar,MPIdata)
285 : :
286 [ + + ]: 6 : if (iorder.eq.1) then
287 [ + - ]: 2 : if (Invar%readifc.ne.1) then
288 : 2 : call tdep_calc_phi1(Solver,Shell1at,Sym,Model%Phi1)
289 : 2 : call tdep_calc_phi2(Solver,Shell2at,Sym,Model%Phi2%SR)
290 : : end if
291 : 2 : call tdep_calc_ftot2(Model,Invar,Model%Phi2%SR,MD%ucart)
292 : :
293 [ + + ]: 4 : else if (iorder.eq.2) then
294 : 2 : call tdep_calc_phi3ref(Solver,Shell3at,Model%Phi3)
295 : 2 : call tdep_calc_ftot3(Model,Invar,Shell3at,MD%ucart,Sym)
296 : :
297 [ + - ]: 2 : else if (iorder.eq.3) then
298 : 2 : call tdep_calc_phi4ref(Solver,Shell4at,Model%Phi4)
299 : 2 : call tdep_calc_ftot4(Model,Invar,Shell4at,MD%ucart,Sym)
300 : : end if
301 : :
302 : : ! Remove forces computed with the previous order
303 : 8 : call tdep_solver_set_residual_forces(Solver, MD, Model)
304 : :
305 : : end do ! iorder
306 : : end if
307 : :
308 : : ! Add the long-range part of the IFC
309 [ + + ]: 88 : if (Invar%loto) then
310 [ + - + - : 1847444 : Model%Phi2%Tot = Model%Phi2%LR + Model%Phi2%SR
+ - - - +
+ + + ]
311 : : end if
312 : :
313 : : ! Free some memory
314 : 88 : call tdep_solver_free(Solver)
315 : 88 : call tdep_destroy_shell(Shell1at)
316 : :
317 : : !==========================================================================================
318 : : !=================== Write the IFC and check the constraints ==============================
319 : : !==========================================================================================
320 : 88 : call tdep_write_phi1(Invar,Model%Phi1)
321 : 88 : call tdep_write_phi2(MD%distance,Invar,MPIdata,Model%Phi2%SR,Shell2at)
322 [ + + ]: 88 : if (Invar%order.ge.3) then
323 : 28 : call tdep_write_phi3(MD%distance,Invar,Model%Phi3,Shell3at,Sym)
324 : : end if
325 [ + + ]: 88 : if (Invar%order.ge.4) then
326 : 16 : call tdep_write_phi4(MD%distance,Invar,Model%Phi4,Shell4at,Sym)
327 : : end if
328 : :
329 : 88 : call tdep_check_constraints(Model,MD%distance,Invar,Sym,Shell3at,Shell4at)
330 : :
331 : : !==========================================================================================
332 : : !===================== Convert Phi2 into IFC object =======================================
333 : : !==========================================================================================
334 : :
335 : : call tdep_ifc2phi2(Ifc%dipdip,Ifc,Invar,Lattice,Invar%natom_unitcell,0,&
336 : 88 : & Model%Phi2,MD%Rlatt_scaled,Shell2at,Sym)
337 : :
338 : : !==========================================================================================
339 : : !===================== Compute the phonon spectrum, the DOS, ==============================
340 : : !===================== the dynamical matrix and write them ===============================
341 : : !==========================================================================================
342 : 88 : call tdep_init_eigen2nd(Eigen2nd_MP,Invar%natom_unitcell,Qbz%nqbz)
343 : 88 : call tdep_init_eigen2nd(Eigen2nd_path,Invar%natom_unitcell,Qpt%nqpt)
344 : :
345 : : call tdep_calc_phdos(Crystal,DDB,Eigen2nd_MP,Eigen2nd_path,Ifc,Invar,Lattice,MPIdata,Invar%natom,&
346 : 88 : & Invar%natom_unitcell,Model%Phi2,PHdos,Qbz,Qpt,MD%Rlatt_scaled,Shell2at,Sym)
347 : 88 : call tdep_destroy_shell(Shell2at)
348 : :
349 : :
350 : : ! Create a new DDB with the coarse q-point grid in the IBZ.
351 : 88 : call DDB%free()
352 : 88 : call Ifc%to_ddb(DDB,Crystal)
353 : 88 : call tdep_write_ddb(DDB,Crystal,Invar)
354 : :
355 : 88 : write(stdout,'(a)') ' See the dij.dat, omega.dat and eigenvectors files'
356 : 88 : write(stdout,'(a)') ' See also the DDB file'
357 : :
358 : : !==========================================================================================
359 : : !===================== Compute the elastic constants ======================================
360 : : !==========================================================================================
361 : 88 : call tdep_calc_elastic(Model%Phi2%SR,MD%distance,Invar,Lattice)
362 : :
363 : : !==========================================================================================
364 : : !=========== Compute U_0, the "free energy" from the model ===============
365 : : !==========================================================================================
366 : 88 : call tdep_calc_model(Model,MD,Invar,MPIdata)
367 : :
368 : : !==========================================================================================
369 : : !===================== Compute the thermodynamical quantities =============================
370 : : !==========================================================================================
371 : 88 : call tdep_calc_thermo(Invar,Lattice,MPIdata,PHdos,Model%U0)
372 : 88 : call PHdos%free()
373 : :
374 : : !==========================================================================================
375 : : !===================== CALCULATION OF THE 3rd ORDER =======================================
376 : : !==========================================================================================
377 [ + + ]: 88 : if (Invar%order>2) then
378 [ + - ]: 28 : if (MPIdata%iam_master) then
379 : 28 : call tdep_write_gruneisen(MD%distance,Eigen2nd_path,Invar,Model%Phi3,Qpt,MD%Rlatt_cart,Shell3at,Sym)
380 : : end if
381 : 28 : call tdep_calc_alpha_gamma(MD%distance,Eigen2nd_MP,Invar,Lattice,MPIdata,Model%Phi3,Qbz,MD%Rlatt_cart,Shell3at,Sym)
382 : :
383 : : !FB call tdep_calc_lifetime1(Crystal,MD%distance,Eigen2nd_MP,Ifc,Invar,Lattice,Model%Phi3,Qbz,MD%Rlatt_cart,Shell3at,Sym)
384 : : end if
385 : :
386 : : !==========================================================================================
387 : : !===================== Free memore ========================================================
388 : : !==========================================================================================
389 : :
390 : 88 : call tdep_destroy_eigen2nd(Eigen2nd_path)
391 : 88 : call tdep_destroy_eigen2nd(Eigen2nd_MP)
392 : :
393 [ + + ]: 88 : if (Invar%order>2) then
394 : 28 : call tdep_destroy_shell(Shell3at)
395 [ + + ]: 28 : if (Invar%order.eq.4) then
396 : 16 : call tdep_destroy_shell(Shell4at)
397 : : end if
398 : : end if
399 : :
400 : 88 : call Ifc%free()
401 : 88 : call DDB%free()
402 : 88 : call Crystal%free()
403 : 88 : call tdep_destroy_sym(Sym)
404 : 88 : call tdep_destroy_qbz(Qbz)
405 : 88 : call tdep_destroy_qpt(Qpt)
406 : 88 : call tdep_model_free(Model)
407 : 88 : call tdep_sampling_free(MD)
408 : 88 : call tdep_destroy_invar(Invar)
409 : 88 : call tdep_destroy_mpidata(MPIdata)
410 : :
411 : : !==========================================================================================
412 : : !===================== End the calculation ================================================
413 : : !==========================================================================================
414 : :
415 : 88 : call tdep_print_Aknowledgments(stdout)
416 : 88 : call flush_unit(stdout)
417 : :
418 : 88 : call timein(tcpu, twall)
419 : 88 : tsec(1)=tcpu-tcpui; tsec(2)=twall-twalli
420 : :
421 : :
422 : : ! Write YAML document with the final summary.
423 [ + - ]: 88 : if (MPIdata%iam_master) then
424 : 88 : write(stdlog, "(a)")""
425 : 88 : write(stdlog, "(a)")"--- !FinalSummary"
426 : 88 : write(stdlog, "(a)")"program: atdep"
427 : 88 : write(stdlog, "(2a)")"version: ", trim(version_string)
428 : 88 : write(stdlog, "(2a)")"start_datetime: ", start_datetime
429 : 88 : write(stdlog, "(2a)")"end_datetime: ", asctime()
430 : 88 : write(stdlog, "(a, f13.1)")"overall_cpu_time: ", tsec(1)
431 : 88 : write(stdlog, "(a, f13.1)")"overall_wall_time: ", tsec(2)
432 : 88 : write(stdlog, "(a, i0)")"mpi_procs: ", MPIdata%nproc
433 : 88 : write(stdlog, "(a)")"..."
434 : 88 : call flush_unit(stdlog)
435 : : end if
436 : :
437 : : !Memory analysis
438 [ + - ]: 88 : call abinit_doctor(trim(Invar%output_prefix), print_mem_report=print_mem_report)
439 : 88 : call flush_unit(stdlog)
440 : 88 : close(unit=stdout)
441 : 88 : 100 call xmpi_end()
442 : :
443 : 0 : end program atdep
444 : : !!***
|
