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