Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_sampling
2 : : !! NAME
3 : : !! m_tdep_sampling
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the TDEP Sampling data type
7 : : !! which holds the set of configurations from which the IFC will be fit.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.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 : : module m_tdep_sampling
24 : :
25 : : use defs_basis
26 : : use m_errors
27 : : use m_abicore
28 : : use m_xmpi
29 : : use m_io_tools
30 : : use m_abihist, only : abihist
31 : : use m_matrix, only : matr3inv
32 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
33 : : use m_tdep_latt, only : Lattice_type, tdep_make_inbox
34 : : use m_tdep_supercell, only : tdep_Supercell_type
35 : :
36 : : implicit none
37 : : !!***
38 : :
39 : : !----------------------------------------------------------------------
40 : :
41 : : !!****t* m_tdep_sampling/tdep_Sampling_type
42 : : !! NAME
43 : : !! tdep_Sampling_type
44 : : !!
45 : : !! FUNCTION
46 : : !! Holds the set of configurations (positions, forces, energies and weights)
47 : : !! sampled along the MD trajectory, from which the IFC are fitted,
48 : : !! together with the derived displacements.
49 : : !!
50 : : !! SOURCE
51 : :
52 : : type tdep_Sampling_type
53 : :
54 : : integer :: natom
55 : : ! Number of atoms in the supercell
56 : :
57 : : integer :: natom_unitcell
58 : : ! Number of atoms in the unitcell
59 : :
60 : : integer :: my_nstep
61 : : ! Number of MD steps held locally
62 : :
63 : : integer :: nstep_tot
64 : : ! Total of MD steps
65 : :
66 : : integer, allocatable :: typat_unitcell(:)
67 : : ! typat_unitcell(natom_unitcell)
68 : : ! Atom type in the unitcell.
69 : :
70 : : integer, allocatable :: typat(:)
71 : : ! typat(natom)
72 : : ! Atom type in the supercell.
73 : :
74 : : !integer, allocatable :: iat_unitcell(:)
75 : : ! iat_unitcell(natom)
76 : : ! Atom index in the unitcell
77 : :
78 : : !integer, allocatable :: icell_supercell(:)
79 : : ! icell_supercell(natom)
80 : : ! Index of the supercell lattice vector.
81 : :
82 : : !integer, allocatable :: iat_supercell(:) ! GA: Not necessary
83 : : ! iat_supercell(natom)
84 : : ! Ordering of the atoms in the supercell (different from input ordering)
85 : :
86 : : double precision :: rprimd(3,3)
87 : : ! Primitve vectors of the supercell, in which xred is expressed.
88 : :
89 : : double precision, allocatable :: xred(:,:,:)
90 : : ! xred(3, natom, my_nstep)
91 : : ! Reduced position of the atom at each step,
92 : :
93 : : double precision, allocatable :: xred_average(:,:)
94 : : ! xred_average(3, natom)
95 : : ! Average position of the atoms
96 : :
97 : : double precision, allocatable :: fcart(:,:,:)
98 : : ! fcart(3, natom, my_nstep)
99 : : ! Cartesian forces at each step.
100 : :
101 : : double precision, allocatable :: etot(:)
102 : : ! etot(my_nstep)
103 : : ! Total energy at each step.
104 : :
105 : : double precision, allocatable :: weights(:)
106 : : ! weights(my_nstep)
107 : : ! The weight of each configuration for the fitting.
108 : : ! By default, these will be 1 / nstep_tot.
109 : :
110 : : double precision, allocatable :: Rlatt_cart(:,:,:)
111 : : ! Rlatt_cart(3, natom_unitcell, natom)
112 : : ! Cartesian lattice vector between a unitcell reference atom
113 : : ! and the ideal site of each atom.
114 : :
115 : : double precision, allocatable :: Rlatt_scaled(:,:,:)
116 : : ! Rlatt_scaled(3, natom_unitcell, natom)
117 : : ! Rlatt_cart divided by acell_unitcell.
118 : : ! May be non-integer for non-orthogonal lattice vectors.
119 : :
120 : : double precision, allocatable :: ucart(:,:,:)
121 : : ! ucart(3, natom, my_nstep)
122 : : ! Cartesian displacements of the atoms with respect to their equilibrium
123 : : ! positions (or mean displacement) at each time step.
124 : :
125 : : double precision, allocatable :: ucart0(:,:)
126 : : ! ucart0(3, natom)
127 : : ! Mean Cartesian displacements of the atoms with respect to their equilibrium position.
128 : :
129 : : double precision, allocatable :: Forces(:)
130 : : ! Forces(3*natom*my_nstep)
131 : : ! The cartesian forces for all configurations, as a flat array.
132 : :
133 : : contains
134 : :
135 : : procedure :: init => tdep_sampling_init_read
136 : : procedure :: free => tdep_sampling_free
137 : : procedure :: get_flat_index => tdep_sampling_get_flat_index
138 : : procedure :: shift_xred => tdep_sampling_shift_xred
139 : : procedure :: reorder => tdep_sampling_reorder
140 : : procedure :: rotate => tdep_sampling_rotate
141 : : procedure :: gather_inbox => tdep_sampling_gather_inbox
142 : :
143 : : end type tdep_Sampling_type
144 : : !!***
145 : :
146 : : public :: tdep_match_supercell
147 : : public :: tdep_match_unitcell
148 : : public :: tdep_write_xred_average_ideal
149 : : public :: tdep_write_ucart_fcart
150 : : public :: tdep_write_distance_rlatt
151 : :
152 : : contains
153 : :
154 : : !=====================================================================================================
155 : :
156 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_init_read
157 : : !! NAME
158 : : !! tdep_sampling_init_read
159 : : !!
160 : : !! FUNCTION
161 : : !! Initialize the TDEP Sampling object: allocate its arrays, then read the
162 : : !! reduced positions, cartesian forces and total energies of the MD
163 : : !! trajectory, either from the HIST.nc netcdf file or from the _xred.dat,
164 : : !! _fcart.dat and _etot.dat ASCII files. Only the steps between
165 : : !! Invar%nstep_min and Invar%nstep_max (every Invar%slice steps) that
166 : : !! belong to the present MPI process are stored.
167 : : !! The weights of the configurations are either read from the _weights.dat
168 : : !! file (Invar%use_weights=1) or set to 1/nstep_tot.
169 : : !! Finally, a global shift is applied to the reduced positions so that the
170 : : !! first atom of the unitcell (Invar%xred_unitcell) is anchored at the origin.
171 : : !!
172 : : !! INPUTS
173 : : !! Invar = Input object containing the input variables.
174 : : !! MPIdata = Info on MPI parallelism.
175 : : !! Hist = abihist object containing the MD trajectory
176 : : !! (only used if Invar%netcdf is true).
177 : : !!
178 : : !! OUTPUT
179 : : !!
180 : : !! SIDE EFFECTS
181 : : !! MD = TDEP Sampling object. On exit, its arrays are allocated and
182 : : !! MD%typat, MD%typat_unitcell, MD%xred, MD%fcart, MD%etot and
183 : : !! MD%weights are filled.
184 : : !!
185 : : !! SOURCE
186 : :
187 : 46 : subroutine tdep_sampling_init_read(MD,Invar,MPIdata,Hist)
188 : :
189 : : class(tdep_Sampling_type), intent(inout) :: MD
190 : : type(atdep_dataset_type), intent(in) :: Invar
191 : : type(MPI_enreg_type), intent(in) :: MPIdata
192 : : type(abihist), intent(in) :: Hist
193 : :
194 : : integer :: this_istep,istep,iatom,jstep
195 : : double precision :: tmp1,tmp2,tmp3
196 : : real(dp) :: multiplicity(3,3), multiplicity_invT(3,3)
197 : : real(dp) :: shift_unit(3), shift_scaled(3)
198 : :
199 : : ! FIXME Check that the number of step requested is not larger than available steps.
200 : :
201 : 46 : MD%natom = Invar%natom
202 : 46 : MD%natom_unitcell = Invar%natom_unitcell
203 : 46 : MD%my_nstep = Invar%my_nstep
204 : 46 : MD%nstep_tot = Invar%nstep_tot
205 : :
206 [ - + + - : 248 : ABI_CALLOC(MD%typat_unitcell, (MD%natom_unitcell))
+ - + - -
+ - + - +
+ + ]
207 [ - + + - : 7416 : ABI_CALLOC(MD%typat, (MD%natom))
+ - + - -
+ - + - +
+ + ]
208 [ - + + - : 589464 : ABI_CALLOC(MD%xred, (3,MD%natom,MD%my_nstep))
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
209 [ - + + - : 138 : ABI_MALLOC(MD%xred_average, (3,MD%natom))
+ - + - -
+ - + -
+ ]
210 [ - + + - : 589418 : ABI_CALLOC(MD%fcart, (3,MD%natom,MD%my_nstep))
- + + - -
+ - + - +
+ + + + +
+ ]
211 [ - + + - : 1130 : ABI_CALLOC(MD%etot, (MD%my_nstep))
+ - + - -
+ - + - +
+ + ]
212 [ - + + - : 1084 : ABI_CALLOC(MD%weights, (MD%my_nstep))
- + - + -
+ + + ]
213 [ - + + - : 70174 : ABI_CALLOC(MD%Rlatt_scaled, (3,MD%natom_unitcell,MD%natom))
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
214 [ - + + - : 70128 : ABI_CALLOC(MD%Rlatt_cart, (3,MD%natom_unitcell,MD%natom))
- + + - -
+ - + - +
+ + + + +
+ ]
215 [ - + + - : 589418 : ABI_CALLOC(MD%ucart, (3,MD%natom,MD%my_nstep))
- + + - -
+ - + - +
+ + + + +
+ ]
216 [ - + + - : 29204 : ABI_CALLOC(MD%ucart0, (3,MD%natom))
- + - + -
+ + + +
+ ]
217 [ + - - + : 441354 : ABI_CALLOC(MD%Forces, (3*MD%natom*MD%my_nstep))
+ - + - +
- - + - +
- + + + ]
218 : :
219 : : !ABI_CALLOC(MD%iat_unitcell, (MD%natom))
220 : : !ABI_CALLOC(MD%icell_supercell, (MD%natom))
221 : :
222 [ + + ]: 156 : MD%typat_unitcell(:) = Invar%typat_unitcell(:)
223 [ + + ]: 7324 : MD%typat(:) = Invar%typat(:)
224 : :
225 : : ! Read xred.dat, fcart.dat and etot.dat ASCII files or extract them from the HIST.nc netcdf file.
226 : 46 : write(ab_out,'(a)') ' '
227 : 46 : this_istep=0
228 : 46 : jstep=0
229 [ + + ]: 46 : if (Invar%use_weights.eq.1) then
230 : 1 : open(unit=30,file=trim(Invar%input_prefix)//'_weights.dat')
231 [ + - ]: 45 : else if (Invar%use_weights.eq.0) then
232 [ + + ]: 1017 : MD%weights=1.0d0/real(MD%nstep_tot)
233 : : endif
234 [ + + ]: 46 : if (Invar%netcdf) then
235 [ + + ]: 238 : do istep=Invar%nstep_min,Invar%nstep_max
236 [ + + ]: 238 : if (mod(istep-Invar%nstep_min,Invar%slice).ne.0) then
237 : : cycle
238 : : else
239 : 221 : jstep=jstep+1
240 [ - + ]: 221 : if (.not.MPIdata%my_step(jstep)) cycle
241 : 221 : this_istep=this_istep+1
242 [ + + + + ]: 79741 : MD%xred(:,:,this_istep) =Hist%xred (:,:,istep)
243 [ + + + + ]: 79741 : MD%fcart(:,:,this_istep)=Hist%fcart(:,:,istep)
244 : 231 : MD%etot(this_istep) =Hist%etot (istep)
245 : : end if
246 : : end do !istep
247 : 7 : write(ab_out,'(a)') ' The positions, forces and energies are extracted from the NetCDF file: HIST.nc'
248 : : else
249 : 39 : open(unit=60,file=trim(Invar%input_prefix)//'_fcart.dat')
250 : 39 : open(unit=50,file=trim(Invar%input_prefix)//'_xred.dat')
251 : 39 : open(unit=40,file=trim(Invar%input_prefix)//'_etot.dat')
252 [ - + ]: 39 : do istep=1,Invar%nstep_min-1
253 [ # # ]: 0 : if (Invar%use_weights.eq.1) then
254 : 0 : read(30,*) tmp1
255 : : endif
256 : 0 : read(40,*) tmp1
257 [ - - ]: 39 : do iatom=1,MD%natom
258 : 0 : read(50,*) tmp1,tmp2,tmp3
259 : 0 : read(60,*) tmp1,tmp2,tmp3
260 : : end do
261 : : end do
262 [ + + ]: 820 : do istep=Invar%nstep_min,Invar%nstep_max
263 [ + + ]: 820 : if (mod(istep-Invar%nstep_min,Invar%slice).ne.0) then
264 [ - + ]: 10 : if (Invar%use_weights.eq.1) then
265 : 0 : read(30,*) tmp1
266 : : endif
267 : 10 : read(40,*) tmp1
268 [ + + ]: 1290 : do iatom=1,Invar%natom
269 : 1280 : read(50,*) tmp1,tmp2,tmp3
270 : 1290 : read(60,*) tmp1,tmp2,tmp3
271 : : end do
272 : : else
273 : 771 : jstep=jstep+1
274 [ - + ]: 771 : if (.not.MPIdata%my_step(jstep)) then
275 [ # # ]: 0 : if (Invar%use_weights.eq.1) then
276 : 0 : read(30,*) tmp1
277 : : endif
278 : 0 : read(40,*) tmp1
279 [ # # ]: 0 : do iatom=1,Invar%natom
280 : 0 : read(50,*) tmp1,tmp2,tmp3
281 : 0 : read(60,*) tmp1,tmp2,tmp3
282 : : end do
283 : : else
284 : 771 : this_istep=this_istep+1
285 [ + + ]: 771 : if (Invar%use_weights.eq.1) then
286 : 20 : read(30,*) MD%weights(this_istep)
287 : : endif
288 : 771 : read(40,*) MD%etot(this_istep)
289 [ + + ]: 127963 : do iatom=1,MD%natom
290 : 127192 : read(50,*) MD%xred (1,iatom,this_istep),MD%xred (2,iatom,this_istep),MD%xred (3,iatom,this_istep)
291 : 127963 : read(60,*) MD%fcart(1,iatom,this_istep),MD%fcart(2,iatom,this_istep),MD%fcart(3,iatom,this_istep)
292 : : end do
293 : : end if !my_step
294 : : end if !slice
295 : : end do !istep
296 : 39 : close(40)
297 : 39 : close(50)
298 : 39 : close(60)
299 : 39 : write(ab_out,'(2a)') ' The positions, forces and energies are extracted from the ASCII files:',&
300 : 78 : & ' xred.dat, fcart.dat & etot.dat'
301 : : end if !netcdf
302 [ + + ]: 46 : if (Invar%use_weights.eq.1) then
303 : 1 : close(30)
304 : : end if
305 : :
306 : : ! Apply xred_unitcell shift to anchor atom 1 at the origin.
307 [ + + + + ]: 598 : multiplicity(:,:) = real(NINT(Invar%multiplicity(:,:)), kind=dp)
308 : 46 : call matr3inv(multiplicity, multiplicity_invT)
309 [ + + ]: 184 : shift_unit(:) = Invar%xred_unitcell(:,1) - Invar%xred_unitcell_original(:,1)
310 [ + + + + ]: 598 : shift_scaled = matmul(multiplicity_invT, shift_unit)
311 [ + + ]: 1038 : do istep=1,MD%my_nstep
312 [ + + ]: 148110 : do iatom=1,MD%natom
313 [ + + ]: 589280 : MD%xred(:,iatom,istep) = MD%xred(:,iatom,istep) + shift_scaled(:)
314 : : end do
315 : : end do
316 : :
317 : 46 : end subroutine tdep_sampling_init_read
318 : : !!***
319 : :
320 : : !=====================================================================================================
321 : :
322 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_free
323 : : !! NAME
324 : : !! tdep_sampling_free
325 : : !!
326 : : !! FUNCTION
327 : : !! Free the memory allocated in the TDEP Sampling object.
328 : : !!
329 : : !! SIDE EFFECTS
330 : : !! MD = TDEP Sampling object. All its allocated arrays are deallocated.
331 : : !!
332 : : !! SOURCE
333 : :
334 : 46 : subroutine tdep_sampling_free(MD)
335 : :
336 : : class(tdep_Sampling_type), intent(inout) :: MD
337 : :
338 [ - + ]: 46 : ABI_FREE(MD%typat)
339 [ - + ]: 46 : ABI_FREE(MD%typat_unitcell)
340 [ - + ]: 46 : ABI_FREE(MD%fcart)
341 [ - + ]: 46 : ABI_FREE(MD%etot)
342 [ - + ]: 46 : ABI_FREE(MD%weights)
343 [ - + ]: 46 : ABI_FREE(MD%Rlatt_scaled)
344 [ - + ]: 46 : ABI_FREE(MD%Rlatt_cart)
345 [ - + ]: 46 : ABI_FREE(MD%ucart)
346 [ - + ]: 46 : ABI_FREE(MD%ucart0)
347 [ - + ]: 46 : ABI_FREE(MD%Forces)
348 : :
349 [ - + ]: 46 : ABI_SFREE(MD%xred)
350 [ + - ]: 46 : ABI_SFREE(MD%xred_average)
351 : :
352 : 46 : end subroutine tdep_sampling_free
353 : : !!***
354 : :
355 : : !=====================================================================================================
356 : :
357 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_get_flat_index
358 : : !! NAME
359 : : !! tdep_sampling_get_flat_index
360 : : !!
361 : : !! FUNCTION
362 : : !! Return the index, in the flat arrays of dimension 3*natom*my_nstep
363 : : !! (such as MD%Forces), corresponding to a given cartesian direction,
364 : : !! atom and (local) MD step.
365 : : !!
366 : : !! INPUTS
367 : : !! MD = TDEP Sampling object.
368 : : !! icart = Cartesian direction (1 to 3).
369 : : !! iatom = Atom index (1 to natom).
370 : : !! istep = Local MD step index (1 to my_nstep).
371 : : !!
372 : : !! OUTPUT
373 : : !! ii = Flat index: icart + 3*(iatom-1) + 3*natom*(istep-1).
374 : : !!
375 : : !! SOURCE
376 : :
377 : 916992 : function tdep_sampling_get_flat_index(MD,icart,iatom,istep) result(ii)
378 : : class(tdep_Sampling_type), intent(in) :: MD
379 : : integer, intent(in) :: icart, iatom, istep
380 : : integer :: ii
381 : 916992 : ii = icart + 3*(iatom-1) + 3*MD%natom*(istep-1)
382 : 475776 : end function tdep_sampling_get_flat_index
383 : : !!***
384 : :
385 : : !=====================================================================================================
386 : :
387 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_shift_xred
388 : : !! NAME
389 : : !! tdep_sampling_shift_xred
390 : : !!
391 : : !! FUNCTION
392 : : !! Apply a constant shift to all reduced positions (MD%xred, and
393 : : !! MD%xred_average if allocated), at every MD step.
394 : : !!
395 : : !! INPUTS
396 : : !! shift(3) = Shift to apply, in reduced coordinates of the supercell.
397 : : !!
398 : : !! OUTPUT
399 : : !!
400 : : !! SIDE EFFECTS
401 : : !! MD = TDEP Sampling object. MD%xred and MD%xred_average are shifted.
402 : : !!
403 : : !! SOURCE
404 : :
405 : 46 : subroutine tdep_sampling_shift_xred(MD,shift)
406 : :
407 : : class(tdep_Sampling_type), intent(inout) :: MD
408 : : real(dp), intent(in) :: shift(3)
409 : : integer :: iatom,istep
410 : :
411 [ + + ]: 1038 : do istep=1, MD%my_nstep
412 [ + + ]: 148110 : do iatom=1,MD%natom
413 [ + + ]: 589280 : MD%xred(:,iatom,istep) = MD%xred(:,iatom,istep) + shift(:)
414 : : end do
415 : : end do
416 : :
417 [ + - ]: 46 : if (allocated(MD%xred_average)) then
418 [ + + ]: 7324 : do iatom=1,MD%natom
419 [ + + ]: 29158 : MD%xred_average(:,iatom) = MD%xred_average(:,iatom) + shift(:)
420 : : end do
421 : : end if
422 : :
423 : : ! Ideal positions are a property of SC (computed on the fly from it), not of MD, so nothing else to update here.
424 : :
425 : 46 : end subroutine tdep_sampling_shift_xred
426 : : !!***
427 : :
428 : : !=====================================================================================================
429 : :
430 : : !!****f* ABINIT/m_tdep_sampling/tdep_match_unitcell
431 : : !! NAME
432 : : !! tdep_match_unitcell
433 : : !!
434 : : !! FUNCTION
435 : : !! Find and apply a global shift to all xred so that the unitcell is
436 : : !! found within the set of xred_average.
437 : : !!
438 : : !! Every raw atom is tested against every unitcell site that shares its type.
439 : : !! Testing every atom against every type-matching site guarantees the correct
440 : : !! cluster wins outright: every atom contributes exactly one vote for the true
441 : : !! global shift (via its own true site), so that cluster totals natom
442 : : !! votes, strictly more than any single wrong (site,tested-site) mismatch
443 : : !! can gather (each such mismatch only ever collects votes from the one
444 : : !! subset of atoms that are truly the mismatched site).
445 : : !!
446 : : !! The clustering itself is done in UNITCELL-fractional coordinates, not
447 : : !! MD%xred_average's own (supercell-fractional) ones: a unitcell site
448 : : !! repeats once per PRIMITIVE cell, i.e. many times per supercell, and two
449 : : !! raw atoms that are both (a periodic copy of) the same unitcell site
450 : : !! generally differ by a fraction of a supercell period, not by a whole
451 : : !! one -- only in unitcell-fractional coordinates does "belongs to the
452 : : !! same unitcell site" reduce to a plain equality modulo 1.
453 : : !!
454 : : !! INPUTS
455 : : !! SC = Supercell object describing the ideal structure
456 : : !! (SC%xred_unitcell is used).
457 : : !!
458 : : !! OUTPUT
459 : : !!
460 : : !! SIDE EFFECTS
461 : : !! MD = TDEP Sampling object. On entry, MD%xred_average must already be
462 : : !! computed. On exit, MD%xred and MD%xred_average are shifted
463 : : !! (see tdep_sampling_shift_xred).
464 : : !!
465 : : !! SOURCE
466 : :
467 : 46 : subroutine tdep_match_unitcell(MD, SC)
468 : :
469 : : type(tdep_Sampling_type), intent(inout) :: MD
470 : : type(tdep_Supercell_type), intent(in) :: SC
471 : :
472 : : integer :: iatcell, jatcell, islot, iatom, icluster, ibest, ncluster, iatom_ref, islot_ref, navg
473 : : integer :: nclustermax
474 : : real(dp) :: tol, mindist, dnorm
475 : : real(dp) :: shift(3), tmp(3), diff(3), xred_unit(3), shift_ref(3), shift_sum(3)
476 : 46 : real(dp), allocatable :: cluster_shift(:,:)
477 : 46 : integer, allocatable :: cluster_count(:)
478 : 46 : integer, allocatable :: cluster_atom(:), cluster_slot(:)
479 : :
480 : : ! Clustering tolerance: half the smallest distance (unitcell-fractional,
481 : : ! minimum-image) separating any two distinct unitcell sites that share the
482 : : ! SAME type. This guarantees two genuinely distinct (and confusable, i.e.
483 : : ! same-type) unitcell sites are never merged into the same shift cluster,
484 : : ! while remaining generous relative to thermal noise on xred_average.
485 : : ! Sites of different types are never compared against each other in the
486 : : ! clustering loop below (the type check there rules that out), so they
487 : : ! must not constrain the tolerance either -- e.g. two unitcell sites of
488 : : ! different species placed at the very same reduced position (a legitimate
489 : : ! input) would otherwise force mindist, and hence tol, to zero.
490 : 46 : mindist = greatest_real
491 [ + + ]: 156 : do iatcell=1,MD%natom_unitcell
492 [ + + ]: 362 : do jatcell=iatcell+1,MD%natom_unitcell
493 [ + + ]: 206 : if (MD%typat_unitcell(jatcell) /= MD%typat_unitcell(iatcell)) cycle
494 [ + + ]: 744 : tmp = SC%xred_unitcell(:,jatcell) - SC%xred_unitcell(:,iatcell)
495 : 186 : call tdep_make_inbox(tmp,1,tol8)
496 [ + + ]: 874 : mindist = min(mindist, sqrt(sum(tmp**2)))
497 : : end do
498 : : end do
499 : : ! No two unitcell sites share a type (most common case, e.g. natom_unitcell==1):
500 : : ! there is nothing a candidate shift could be confused with, so a generous
501 : : ! tolerance is safe -- it only needs to comfortably exceed thermal noise.
502 [ + + ]: 46 : if (mindist == greatest_real) mindist = 0.8
503 : 46 : tol = half * mindist
504 : :
505 : : ! Cluster the candidate global shifts (in unitcell-fractional coordinates)
506 : : ! implied by every (raw atom, type-matching unitcell site) pairing:
507 : : ! shift_candidate(iatom,islot) = xred_unitcell(islot) - scaled_to_unit(xred_average(iatom)),
508 : : ! wrapped to the minimum image.
509 : 46 : nclustermax = MD%natom_unitcell*MD%natom_unitcell
510 [ - + + - : 138 : ABI_MALLOC(cluster_shift,(3,nclustermax))
+ - + - -
+ - + ]
511 [ - + + - : 138 : ABI_MALLOC(cluster_count,(nclustermax))
+ - + - -
+ - + ]
512 [ - + + - : 92 : ABI_MALLOC(cluster_atom,(nclustermax))
- + - + ]
513 [ - + + - : 92 : ABI_MALLOC(cluster_slot,(nclustermax))
- + - + ]
514 : 46 : ncluster = 0
515 : :
516 [ + + ]: 7324 : do iatom=1,MD%natom
517 : 7278 : xred_unit = SC%scaled_to_unit(MD%xred_average(:,iatom))
518 : :
519 [ + + ]: 23002 : do islot=1,MD%natom_unitcell
520 [ + + ]: 15678 : if (MD%typat(iatom) /= MD%typat_unitcell(islot)) cycle
521 : :
522 [ + + ]: 54072 : tmp = SC%xred_unitcell(:,islot) - xred_unit
523 : 13518 : call tdep_make_inbox(tmp,1,tol8)
524 : :
525 : 13518 : ibest = 0
526 [ + + ]: 49563 : do icluster=1,ncluster
527 : : ! Two independently-wrapped values sitting on
528 : : ! opposite sides of the [-0.5,0.5) boundary (e.g. +0.499 and -0.499)
529 : : ! represent the same cluster, so the DIFFERENCE must be re-wrapped
530 : : ! before measuring its norm.
531 [ + + ]: 197584 : diff = tmp-cluster_shift(:,icluster)
532 : 49396 : call tdep_make_inbox(diff,1,tol8)
533 [ + + ]: 197584 : dnorm = sqrt(sum(diff**2))
534 [ + + ]: 49563 : if (dnorm < tol) then
535 : : ibest = icluster
536 : : exit
537 : : end if
538 : : end do
539 : :
540 [ + + ]: 20796 : if (ibest > 0) then
541 : 13351 : cluster_count(ibest) = cluster_count(ibest) + 1
542 : : else
543 : 167 : ncluster = ncluster + 1
544 [ - + ]: 167 : if (ncluster > nclustermax) then
545 : 0 : ABI_ERROR('Found more distinct (unitcell site, global shift) combinations than expected.')
546 : : end if
547 [ + + ]: 668 : cluster_shift(:,ncluster) = tmp
548 : 167 : cluster_atom(ncluster) = iatom
549 : 167 : cluster_slot(ncluster) = islot
550 : 167 : cluster_count(ncluster) = 1
551 : : end if
552 : : end do
553 : : end do
554 : :
555 [ - + ]: 46 : if (ncluster == 0) then
556 : 0 : ABI_ERROR('Could not find any atom matching any of the unitcell atom types.')
557 : : end if
558 : :
559 [ - + + + : 259 : ibest = maxloc(cluster_count(1:ncluster), dim=1)
+ + ]
560 : 46 : iatom_ref = cluster_atom(ibest)
561 : 46 : islot_ref = cluster_slot(ibest)
562 : :
563 : : ! Compute a first estimate of the shift directly in supercell-fractional
564 : : ! coordinates (MD%xred's own units), from iatom_ref's own raw xred_average
565 : : ! -- NOT by converting cluster_shift(:,ibest) back through unit_to_scaled.
566 : : ! cluster_shift was wrapped to the minimum image in UNITCELL-fractional
567 : : ! coordinates, and a unitcell-fractional wrap boundary (e.g. a basis atom
568 : : ! sitting at exactly 0.5) does not correspond to a negligible ambiguity
569 : : ! once converted to supercell-fractional coordinates: it can differ by as
570 : : ! much as 1/N of a supercell period, N being the relevant HNF diagonal
571 : : ! entry. Recomputing directly in supercell-fractional coordinates and
572 : : ! wrapping only once, there, avoids that.
573 [ + + ]: 184 : shift_ref = SC%unit_to_scaled(SC%xred_unitcell(:,islot_ref)) - MD%xred_average(:,iatom_ref)
574 : 46 : call tdep_make_inbox(shift_ref,1,tol8)
575 : :
576 : : ! Refine shift_ref into the actual average over every (atom,slot) pairing
577 : : ! belonging to the winning cluster, not just iatom_ref alone: a single
578 : : ! atom's own thermal displacement at that instant is otherwise baked
579 : : ! wholesale into the global shift, which is precise enough to place every
580 : : ! atom on the right periodic image (find_nearest_atom's job) but noisier
581 : : ! than necessary as an origin convention -- averaging over the whole
582 : : ! cluster (typically a sizeable fraction of natom) suppresses that noise
583 : : ! by a further sqrt(cluster size).
584 : : !
585 : : ! The average is accumulated ENTIRELY in unitcell-fractional coordinates,
586 : : ! as deviations from cluster_shift(:,ibest) (already guaranteed small by
587 : : ! the cluster membership test itself). Averaging directly in
588 : : ! supercell-fractional coordinates does NOT work: different members of
589 : : ! the same cluster generally sit in different primitive cells, so their
590 : : ! raw supercell-fractional positions differ from each other by
591 : : ! cell-dependent fractions of a supercell period, not by small thermal
592 : : ! noise alone. The (small) averaged unitcell-fractional correction is
593 : : ! converted to supercell-fractional coordinates and added to shift_ref
594 : : ! only at the very end, which is safe precisely because it is small: it
595 : : ! does not, unlike a full shift value, sit near a wrap boundary.
596 : 46 : shift_sum = zero
597 : 46 : navg = 0
598 [ + + ]: 7324 : do iatom=1,MD%natom
599 [ + + ]: 23002 : do islot=1,MD%natom_unitcell
600 [ + + ]: 15678 : if (MD%typat(iatom) /= MD%typat_unitcell(islot)) cycle
601 : :
602 [ + + ]: 54072 : tmp = SC%xred_unitcell(:,islot) - SC%scaled_to_unit(MD%xred_average(:,iatom))
603 : 13518 : call tdep_make_inbox(tmp,1,tol8)
604 [ + + ]: 54072 : diff = tmp-cluster_shift(:,ibest)
605 : 13518 : call tdep_make_inbox(diff,1,tol8)
606 [ + + + + ]: 54072 : if (sqrt(sum(diff**2)) >= tol) cycle ! not a member of the winning cluster
607 : :
608 [ + + ]: 29112 : shift_sum = shift_sum + diff
609 : 22956 : navg = navg + 1
610 : : end do
611 : : end do
612 : :
613 [ + + + + ]: 368 : shift = shift_ref + SC%unit_to_scaled(shift_sum/real(navg,kind=dp))
614 : 46 : call tdep_make_inbox(shift,1,tol8)
615 : :
616 : 46 : call MD%shift_xred(shift)
617 : :
618 [ - + ]: 46 : ABI_FREE(cluster_shift)
619 [ - + ]: 46 : ABI_FREE(cluster_count)
620 : 46 : ABI_FREE(cluster_atom)
621 : 46 : ABI_FREE(cluster_slot)
622 : :
623 : 46 : end subroutine tdep_match_unitcell
624 : : !!***
625 : :
626 : : !=====================================================================================================
627 : :
628 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_reorder
629 : : !! NAME
630 : : !! tdep_sampling_reorder
631 : : !!
632 : : !! FUNCTION
633 : : !! Reorder all populated per-atom arrays held by MD, according to a gather-style index map:
634 : : !! new_array(:,iatom) = old_array(:,order(iatom)).
635 : : !!
636 : : !! Only the atom-indexed quantities that are populated by the time atoms get
637 : : !! reordered in tdep_match_supercell are handled here:
638 : : !! typat, xred_average, xred and fcart.
639 : : !! Quantities computed later (ucart, Forces, ucart0, Rlatt_cart, Rlatt_scaled)
640 : : !! are untouched. (xred_ideal and dij live in SC, which is never reordered.)
641 : : !!
642 : : !! INPUTS
643 : : !! order(natom) = Gather-style index map (see above).
644 : : !!
645 : : !! OUTPUT
646 : : !!
647 : : !! SIDE EFFECTS
648 : : !! MD = TDEP Sampling object. MD%typat, MD%xred_average, MD%xred and
649 : : !! MD%fcart are reordered.
650 : : !!
651 : : !! SOURCE
652 : :
653 : 46 : subroutine tdep_sampling_reorder(MD, order)
654 : :
655 : : class(tdep_Sampling_type),intent(inout) :: MD
656 : : integer,intent(in) :: order(MD%natom)
657 : :
658 : : integer :: istep,iatom
659 : 46 : integer, allocatable :: typat_tmp(:)
660 : 46 : real(dp), allocatable :: vec_tmp(:,:)
661 : :
662 [ - + + - : 138 : ABI_MALLOC(typat_tmp,(MD%natom))
+ - + - -
+ + - ]
663 [ + + ]: 7324 : do iatom=1,MD%natom
664 : 7324 : typat_tmp(iatom) = MD%typat(order(iatom))
665 : : end do
666 [ + + ]: 7324 : MD%typat(:) = typat_tmp(:)
667 [ - + ]: 46 : ABI_FREE(typat_tmp)
668 : :
669 [ - + + - : 138 : ABI_MALLOC(vec_tmp,(3,MD%natom))
+ - + - -
+ - + ]
670 : :
671 [ + + ]: 7324 : do iatom=1,MD%natom
672 [ + + ]: 29158 : vec_tmp(:,iatom) = MD%xred_average(:,order(iatom))
673 : : end do
674 [ + + + + ]: 29158 : MD%xred_average(:,:) = vec_tmp(:,:)
675 : :
676 [ + + ]: 1038 : do istep=1,MD%my_nstep
677 : :
678 [ + + ]: 148064 : do iatom=1,MD%natom
679 [ + + ]: 589280 : vec_tmp(:,iatom) = MD%xred(:,order(iatom),istep)
680 : : end do
681 [ + + + + ]: 589280 : MD%xred(:,:,istep) = vec_tmp(:,:)
682 : :
683 [ + + ]: 148064 : do iatom=1,MD%natom
684 [ + + ]: 589280 : vec_tmp(:,iatom) = MD%fcart(:,order(iatom),istep)
685 : : end do
686 [ + + + + ]: 589326 : MD%fcart(:,:,istep) = vec_tmp(:,:)
687 : :
688 : : end do
689 : :
690 [ - + ]: 46 : ABI_FREE(vec_tmp)
691 : :
692 : 46 : end subroutine tdep_sampling_reorder
693 : : !!***
694 : :
695 : : !=====================================================================================================
696 : :
697 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_gather_inbox
698 : : !! NAME
699 : : !! tdep_sampling_gather_inbox
700 : : !!
701 : : !! FUNCTION
702 : : !! Shift xred to keep atoms in the same unit cell at each step.
703 : : !! For each atom and each step, every reduced coordinate is shifted by
704 : : !! an integer (-1, 0 or 1) so as to be as close as possible to the
705 : : !! corresponding coordinate at the first MD step (which is communicated
706 : : !! to all MPI processes).
707 : : !!
708 : : !! INPUTS
709 : : !! MPIdata = Info on MPI parallelism.
710 : : !!
711 : : !! OUTPUT
712 : : !!
713 : : !! SIDE EFFECTS
714 : : !! MD = TDEP Sampling object. MD%xred is shifted.
715 : : !!
716 : : !! SOURCE
717 : :
718 : 46 : subroutine tdep_sampling_gather_inbox(MD,MPIdata)
719 : :
720 : : class(tdep_Sampling_type), intent(inout) :: MD
721 : : type(MPI_enreg_type), intent(in) :: MPIdata
722 : : integer :: natom,ii,iatom,istep,ierr
723 : : integer :: shift,shift_max,shift_best
724 : : double precision :: xi, dist, best_dist
725 : 46 : double precision, allocatable :: x0(:,:)
726 : :
727 : 46 : natom = MD%natom
728 [ - + + - : 138 : ABI_MALLOC(x0,(3,natom))
+ - + - -
+ - + ]
729 : :
730 : : ! Communicate xred at the first step
731 [ + + + + ]: 29158 : x0(:,:) = zero
732 [ + - ]: 46 : if (MPIdata%my_step(1)) then
733 [ + + + + ]: 29158 : x0(:,:) = MD%xred(:,:,1)
734 : : end if
735 : 46 : call xmpi_sum(x0,MPIdata%comm_step,ierr)
736 : :
737 : : ! Shift xred from all steps in the same unitcell as the first step
738 : 46 : shift_max = 1
739 [ + + ]: 1038 : do istep=1, MD%my_nstep
740 [ + + ]: 148110 : do iatom=1,natom
741 [ + + ]: 589280 : do ii=1,3
742 : 441216 : best_dist = abs(MD%xred(ii,iatom,istep) - x0(ii,iatom))
743 : 441216 : shift_best = 0
744 [ + + ]: 1764864 : do shift=-shift_max,shift_max
745 : 1323648 : xi = MD%xred(ii,iatom,istep) + shift
746 : 1323648 : dist = abs(xi - x0(ii,iatom))
747 [ + + ]: 1764864 : if (dist < best_dist) then
748 : 366 : best_dist = dist
749 : 366 : shift_best = shift
750 : : end if
751 : : end do
752 : 588288 : MD%xred(ii,iatom,istep) = MD%xred(ii,iatom,istep) + shift_best
753 : : end do
754 : : end do
755 : : end do
756 : :
757 [ - + ]: 46 : ABI_FREE(x0)
758 : :
759 : 46 : end subroutine tdep_sampling_gather_inbox
760 : : !!***
761 : :
762 : : !=====================================================================================================
763 : :
764 : : !!****f* ABINIT/m_tdep_sampling/tdep_sampling_rotate
765 : : !! NAME
766 : : !! tdep_sampling_rotate
767 : : !!
768 : : !! FUNCTION
769 : : !! Apply a rotation to the sampling.
770 : : !! Rotation_cart is a rotation matrix that has been applied to each primitive
771 : : !! vectors. Hence, the xred are unchanged, but the forces stored in Cartesian
772 : : !! coordinates must be rotated.
773 : : !!
774 : : !! INPUTS
775 : : !! rotation_cart(3,3) = Rotation matrix, in cartesian coordinates.
776 : : !!
777 : : !! OUTPUT
778 : : !!
779 : : !! SIDE EFFECTS
780 : : !! MD = TDEP Sampling object. MD%fcart is rotated.
781 : : !!
782 : : !! SOURCE
783 : :
784 : 46 : subroutine tdep_sampling_rotate(MD,rotation_cart)
785 : :
786 : : class(tdep_Sampling_type), intent(inout) :: MD
787 : : real(dp), intent(in) :: rotation_cart(3,3)
788 : : integer :: iatom,istep
789 : :
790 : : ! Apply rotation to fcart
791 [ + + ]: 1038 : do istep=1,MD%my_nstep
792 [ + + ]: 148110 : do iatom=1,MD%natom
793 [ + + + + : 2942432 : MD%fcart(:,iatom,istep) = MATMUL(rotation_cart, MD%fcart(:,iatom,istep))
+ + + + +
- + - ]
794 : : end do
795 : : end do
796 : :
797 : 46 : end subroutine tdep_sampling_rotate
798 : : !!***
799 : :
800 : : !====================================================================================================
801 : :
802 : : !!****f* ABINIT/m_tdep_sampling/tdep_match_supercell
803 : : !! NAME
804 : : !! tdep_match_supercell
805 : : !!
806 : : !! FUNCTION
807 : : !! Find the mapping between the atoms of the ideal (equilibrium) supercell SC
808 : : !! and the atoms of the input molecular dynamics (MD) trajectory, using their
809 : : !! average positions. Then compute the atom displacements (ucart) with
810 : : !! respect to the equilibrium positions at each time step of the MD, and a
811 : : !! few geometric quantities of SC (dij, Rlatt_cart/Rlatt_scaled) that are
812 : : !! needed downstream to build the shells and the IFC.
813 : : !!
814 : : !! Concretely, this routine, in order:
815 : : !! 1) Computes the average MD position of each atom (MD%xred_average).
816 : : !! 2) Applies a global shift to all MD positions so that they share the
817 : : !! same unitcell definition as SC (tdep_match_unitcell).
818 : : !! 3) Matches every raw MD atom to its nearest SC site (SC%find_nearest_atom)
819 : : !! and shifts its MD trajectory (xred and xred_average) by a supercell
820 : : !! lattice vector so that it sits on SC's fixed, centered ideal site.
821 : : !! 4) Reorders every per-atom quantity held by MD (typat, xred,
822 : : !! xred_average, fcart) to match the natural atom ordering of SC.
823 : : !! 5) Writes the (raw, unwrapped) ideal and average positions to
824 : : !! xred_average.xyz, and the supercell definition to a CIF file, for
825 : : !! diagnostic purposes.
826 : : !! 6) Computes the cartesian displacements MD%ucart at each MD step, the
827 : : !! flat cartesian force array MD%Forces, and the mean displacement
828 : : !! MD%ucart0. MD%xred is freed at the end of this step, since it is no
829 : : !! longer needed.
830 : : !! 7) Computes SC%dij, the distance (and cartesian vector) between every
831 : : !! pair of atoms in the supercell.
832 : : !! 8) Computes MD%Rlatt_cart and MD%Rlatt_scaled, the cartesian (resp.
833 : : !! scaled) lattice vector between each atom and its parent reference
834 : : !! unitcell atom. Their construction must match the convention expected
835 : : !! by module m_ifc (see the comment further below).
836 : : !!
837 : : !! INPUTS
838 : : !! Invar = Input object containing the input variables.
839 : : !! MPIdata = Info on MPI parallelism.
840 : : !!
841 : : !! OUTPUT
842 : : !!
843 : : !! SIDE EFFECTS
844 : : !! MD = TDEP Sampling object containing the input positions and forces.
845 : : !! Its per-atom quantities are reordered to match the natural
846 : : !! ordering of SC (see step 4 above), and xred_average, ucart,
847 : : !! ucart0, Forces, Rlatt_cart and Rlatt_scaled are computed
848 : : !! (see FUNCTION above for details). xred is freed.
849 : : !! SC = Supercell object describing the ideal structure. SC%dij is
850 : : !! computed (see FUNCTION above). SC itself is never reordered.
851 : : !!
852 : : !! NOTES
853 : : !!
854 : : !! SOURCE
855 : :
856 : 46 : subroutine tdep_match_supercell(Invar,MPIdata,MD,SC)
857 : :
858 : : type(atdep_dataset_type),intent(in) :: Invar
859 : : type(MPI_enreg_type),intent(in) :: MPIdata
860 : : type(tdep_Sampling_type),intent(inout) :: MD
861 : : type(tdep_Supercell_type),intent(inout) :: SC
862 : :
863 : : ! Scalars
864 : : logical :: debug
865 : : integer :: natom, natom_unitcell, nstep
866 : : integer :: ii,jj,iatcell,iatom,jatom,istep
867 : : integer :: ierr
868 : : integer :: icell,iat_unit,iat_super
869 : :
870 : : ! Arrays
871 : : integer :: R_at(3)
872 : : integer :: R_at_centered(3)
873 : : integer :: R_shift(3)
874 : : real(dp) :: R_shift_scaled(3)
875 : : real(dp) :: xred(3),xcart(3)
876 : : real(dp) :: xred_unit(3),xred_scaled(3)
877 : : real(dp) :: ucart(3)
878 : : real(dp) :: rprimd(3,3)
879 : : real(dp) :: xcart_average(3)
880 : : real(dp) :: Rlatt(3)
881 : : real(dp) :: R_cell(3)
882 : : real(dp) :: R_cell_ref(3)
883 : : real(dp) :: xred_ideal_iatom(3)
884 : : real(dp) :: origin(3)
885 : : logical :: freeze_dir(3)
886 : :
887 : 46 : integer, allocatable :: raw_of_natural(:)
888 : 46 : real(dp), allocatable :: xcart_ideal(:,:)
889 : 46 : real(dp), allocatable :: xred_ideal_raw(:,:)
890 : : character(len=500) :: msg
891 : :
892 : 46 : debug = (Invar%prtvol > 1)
893 : :
894 : 46 : write(ab_out,*)' '
895 : 46 : write(ab_out,*) '#############################################################################'
896 : 46 : write(ab_out,*) '###### Find the matching between ideal and average positions ###############'
897 : 46 : write(ab_out,*) '#############################################################################'
898 : :
899 : 46 : nstep = MD%my_nstep
900 : 46 : natom = MD%natom
901 : 46 : natom_unitcell = MD%natom_unitcell
902 : :
903 : : ! ===============================================================================
904 : : ! 1) Compute average positions from MD
905 : : !
906 : 46 : write(ab_out,*)' Computing average positions.'
907 [ + + + + ]: 29158 : MD%xred_average = zero
908 [ + + ]: 1038 : do istep=1,nstep
909 [ + + ]: 148110 : do iatom=1,natom
910 [ + + ]: 589280 : MD%xred_average(:,iatom) = MD%xred_average(:,iatom) + MD%xred(:,iatom,istep)
911 : : end do
912 : : end do
913 : 46 : call xmpi_sum(MD%xred_average,MPIdata%comm_step,ierr)
914 [ + + + + ]: 29158 : MD%xred_average(:,:) = MD%xred_average(:,:) / real(MD%nstep_tot)
915 : :
916 : :
917 : : ! ===============================================================================
918 : : ! 2) Find and apply a global shift to all xred to match the unitcell definition.
919 : : !
920 : : ! GA: This would not be necessary if the unitcell definition in the input file
921 : : ! matched the supercell constructed for the MD trajectory.
922 : : ! We could add this requirement explicitly and avoid this step.
923 : : !
924 : : ! TODO add an input variable to activate this call when needed.
925 : 46 : write(ab_out,*)' Computing global shift to match unitcell definition.'
926 : 46 : call tdep_match_unitcell(MD, SC)
927 : :
928 : :
929 : : ! ===============================================================================
930 : : ! 3) Match raw MD atoms to SC's natural sites and shift MD onto SC's fixed,
931 : : ! centered ideal positions (still in raw order, before any reordering).
932 : :
933 : : ! -----------------------------------------------------------------------------
934 : : ! Prepare second reordering to match legacy tdep_MatchIdeal2Average.
935 : : ! Find, in the raw MD order, the anchor atom that tdep_MatchIdeal2Average
936 : : ! uses to define its own reference order: the first atom whose type
937 : : ! matches the first unitcell atom. This only depends on the (untouched)
938 : : ! input ordering, so it can be found before any matching below.
939 : : !
940 : : !iatom_raw_ref = 0
941 : : !do iatom=1,natom
942 : : ! if (Invar%typat(iatom) == SC%typat_unitcell(1)) then
943 : : ! iatom_raw_ref = iatom
944 : : ! exit
945 : : ! end if
946 : : !end do
947 : : !if (iatom_raw_ref == 0) then
948 : : ! ABI_ERROR('Could not find an atom matching the unitcell basis type.')
949 : : !end if
950 : : ! ---------------------------------------------------------------------
951 : :
952 : 46 : write(ab_out,*)' Matching average atom positions to supercell sites.'
953 : :
954 [ - + + - : 7416 : ABI_CALLOC(raw_of_natural, (natom))
+ - + - -
+ - + +
+ ]
955 : :
956 [ + + ]: 7324 : do iatom=1,natom
957 : :
958 : 7278 : xcart = SC%red_scaled_to_cart(MD%xred_average(:,iatom))
959 : :
960 : 7278 : call SC%find_nearest_atom(xcart, MD%typat(iatom), iat_super, R_at, ucart)
961 : :
962 : : ! Shift the atom by a supercell lattice vector to match the SC definition.
963 : 7278 : R_at_centered = SC%get_R_at(iat_super)
964 [ + + ]: 29112 : R_shift = R_at_centered - R_at
965 [ + + - + ]: 29112 : R_shift_scaled = SC%unit_to_scaled(real(R_shift, kind=dp))
966 : :
967 [ + + ]: 29112 : MD%xred_average(:,iatom) = MD%xred_average(:,iatom) + R_shift_scaled
968 [ + + ]: 154350 : do istep=1,nstep
969 [ + + ]: 595566 : MD%xred(:,iatom,istep) = MD%xred(:,iatom,istep) + R_shift_scaled
970 : : end do
971 : :
972 : : ! --------------------------------------------------------------------
973 : : ! Check that no two atoms of the MD sampling have been assigned
974 : : ! to the same atom in the supercell.
975 [ - + ]: 7278 : if (raw_of_natural(iat_super) .ne. 0) then
976 : 0 : jatom = raw_of_natural(iat_super)
977 : 0 : call SC%get_iat_unitcell(jatom, icell, iat_unit)
978 : 0 : write(msg,*) 'More than one atom assigned to supercell site:',ch10, &
979 : 0 : & 'iatom = ', iatom, ch10,&
980 : 0 : & 'icell, iat_unit = ', icell, iat_unit, ch10,&
981 : 0 : & 'jatom = ', jatom, ch10,&
982 : 0 : & 'xcart (1) = ', xcart,ch10,&
983 [ # # ]: 0 : & 'xcart (2) = ', SC%get_xcart(iat_super)
984 : 0 : ABI_WARNING(msg)
985 : : end if
986 : 7324 : raw_of_natural(iat_super) = iatom
987 : : ! ----------------------------------------------------------------------
988 : :
989 : : end do
990 : :
991 : :
992 : : ! ===============================================================================
993 : : ! 4) Reorder atom positions, types, and forces to match the natural ordering of SC.
994 : : !
995 : 46 : write(ab_out,*)' Reordering atoms in MD trajectory.'
996 : 46 : call MD%reorder(raw_of_natural)
997 : :
998 : : ! ---------------------------------------------------------------------
999 : : !! Second reordering: natural -> reference
1000 : : !! Bring MD into the same atom enumeration as the legacy
1001 : : !! tdep_MatchIdeal2Average routine, anchored on the same physical atom
1002 : : !! (iat_super_ref) that routine would pick as its own reference.
1003 : :
1004 : : !ABI_CALLOC(natural_of_reference, (natom))
1005 : : !natural_of_reference = get_ifc_order(SC, iat_super_ref)
1006 : :
1007 : : !! Apply the reordering to the per-atom quantities currently held by MD
1008 : : !! (so far expressed in the natural SC order), bringing them into the
1009 : : !! reference order: vec_new(:,iatom) = vec_old(:,natural_of_reference(iatom)).
1010 : : !call MD%reorder(natural_of_reference)
1011 : :
1012 : : !! The "original MD order" bookkeeping array is not part of MD, so it is
1013 : : !! reordered separately, using the same gather convention.
1014 : : !ABI_MALLOC(raw_of_natural_tmp, (natom))
1015 : : !do iatom=1,natom
1016 : : ! raw_of_natural_tmp(iatom) = raw_of_natural(natural_of_reference(iatom))
1017 : : !end do
1018 : : !raw_of_natural(:) = raw_of_natural_tmp(:)
1019 : : !ABI_FREE(raw_of_natural_tmp)
1020 : :
1021 : : !! From this point on, all per-atom quantities in MD are in reference
1022 : : !! order. SC itself is never reordered, so any subsequent call into SC
1023 : : !! must translate via natural_of_reference(iatom).
1024 : :
1025 : : ! =============================================================================
1026 : : ! 5) Save ideal positions and write them to file.
1027 : :
1028 : : ! Save ideal positions
1029 [ - + + - : 138 : ABI_MALLOC(xred_ideal_raw,(3,natom))
+ - + - -
+ - + ]
1030 [ + + ]: 7324 : do iatom=1,natom
1031 : 7324 : xred_ideal_raw(:,iatom) = SC%get_xred_scaled(iatom)
1032 : : end do
1033 : :
1034 : : ! Note that xred_ideal_raw, used here, is not wrapped, unlike the ideal
1035 : : ! positions computed on the fly in step 8 (via
1036 : : ! SC%get_xred_scaled(iatom,wrap=.true.)). In previous implementation,
1037 : : ! this wrap happened before writing this file, which created an
1038 : : ! apparent mismatch.
1039 : 46 : write(ab_out,*)' Writing ideal and average positions in file xred_average.xyz'
1040 [ + - ]: 46 : if (MPIdata%iam_master) then
1041 : 46 : call tdep_write_xred_average_ideal(Invar,MD%xred_average,xred_ideal_raw,raw_of_natural)
1042 : : end if
1043 : :
1044 [ - + ]: 46 : ABI_FREE(xred_ideal_raw)
1045 [ - + ]: 46 : ABI_FREE(raw_of_natural)
1046 : :
1047 : : ! Write the supercell in file.
1048 [ + - ]: 46 : if (MPIdata%iam_master) then
1049 : 46 : call SC%write_cif(Invar%output_prefix)
1050 : : end if
1051 : :
1052 : : ! =============================================================================
1053 : : ! 6) Compute cartesian displacements at each MD step.
1054 : :
1055 : 46 : write(ab_out,*)' Computing cartesian displacements.'
1056 : :
1057 [ - + + - : 29204 : ABI_CALLOC(xcart_ideal,(3,natom))
- + - + +
+ + + ]
1058 [ + + ]: 46 : if (Invar%use_ideal_positions.eq.0) then
1059 [ + + ]: 1252 : do iatom=1,natom
1060 : : ! This is in scaled reduced coordinates
1061 [ + + ]: 4968 : xred = MD%xred_average(:,iatom)
1062 : 1242 : xcart = SC%red_scaled_to_cart(xred)
1063 [ + + ]: 4978 : xcart_ideal(:,iatom) = xcart
1064 : : end do
1065 : : else
1066 [ + + ]: 6072 : do iatom=1,natom
1067 : 6072 : xcart_ideal(:,iatom) = SC%get_xcart(iatom)
1068 : : end do
1069 : : end if
1070 : :
1071 : : ! Compute ucart for all MD steps
1072 [ + + ]: 1038 : do istep=1,nstep
1073 [ + + ]: 148110 : do iatom=1,natom
1074 [ + + ]: 588288 : xred = MD%xred(:,iatom,istep)
1075 : 147072 : xcart = SC%red_scaled_to_cart(xred)
1076 [ + + ]: 589280 : MD%ucart(:,iatom,istep) = xcart - xcart_ideal(:,iatom)
1077 : : end do
1078 : : end do
1079 : :
1080 : : ! Rearrange fcart into the Forces vector
1081 [ + + ]: 1038 : do istep=1,nstep
1082 [ + + ]: 148110 : do jatom=1,natom
1083 [ + + ]: 589280 : do ii=1,3
1084 : 441216 : jj = MD%get_flat_index(ii,jatom,istep)
1085 : 588288 : MD%Forces(jj) = MD%fcart(ii,jatom,istep)
1086 : : enddo
1087 : : enddo
1088 : : enddo
1089 : :
1090 : : ! Write ucart/fcart (for debugging)
1091 [ - + - - ]: 46 : if (debug.and.MPIdata%iam_master) then
1092 : 0 : write(ab_out,*)' Writing displacements and forces in file ucart_fcart.dat'
1093 : 0 : call tdep_write_ucart_fcart(Invar,MPIdata,natom,nstep,MD%typat,MD%ucart,MD%fcart)
1094 : : end if
1095 : :
1096 : : ! Save the average displacement, in case we want to reconstruct xcart from ucart
1097 [ + + ]: 46 : if (Invar%use_ideal_positions.eq.0) then
1098 [ + + + + ]: 4978 : MD%ucart0 = zero
1099 : : else
1100 [ + + ]: 6072 : do iatom=1,natom
1101 : 6036 : xcart_average = SC%red_scaled_to_cart(MD%xred_average(:,iatom))
1102 [ + + ]: 24180 : MD%ucart0(:,iatom) = xcart_average(:) - xcart_ideal(:,iatom)
1103 : : end do
1104 : : end if
1105 : :
1106 [ - + ]: 46 : ABI_FREE(xcart_ideal)
1107 : :
1108 : : ! We no longer need xred after this routine
1109 [ - + ]: 46 : ABI_FREE(MD%xred)
1110 : :
1111 : : ! ==========================================================================
1112 : : ! Here we should perform a check on ucart to make sure that it is not larger
1113 : : ! than rprimd, and warn the user when this is the case.
1114 : : ! In that case, unitcell attribution has failed, and the IFC will be wrong.
1115 : : ! One possible cause of this situation is that two atoms have swapped place
1116 : : ! during the MD simulation.
1117 : : ! In that case, we should be able to REPAIR this situation by swapping atom indices
1118 : : ! for certain steps.
1119 : : ! Otherwise, the user should curate their data.
1120 : :
1121 : :
1122 : : ! ==========================================================================
1123 : : ! 7) Compute all interatomic distances in the supercell.
1124 : : ! They are used to initialize the shell objects.
1125 : : !
1126 : : ! GA: Those quantities could be computed on-the-fly.
1127 : :
1128 : 46 : write(ab_out,*)' Computing all interatomic distances.'
1129 : 46 : call SC%calc_dij()
1130 : :
1131 : :
1132 : : ! ==========================================================================
1133 : : ! 8) Compute Rlatt_cart and Rlatt_scaled.
1134 : : ! These are the Cartesian lattice vector between a unitcell reference atom
1135 : : ! and the ideal site of each atom, and will be used in tdep_ifc2phi2
1136 : : ! and tdep_calc_gruneisen.
1137 : : !
1138 : : ! The construction of these arrays is quite convoluted, because they must
1139 : : ! match exactly the convention expected by module m_ifc, which compares
1140 : : ! lattice vectors against its own real-space box
1141 : : ! (Ifc%rcan/Ifc%rpt, built independently by canat9/make_bigbox in m_dynmat).
1142 : : !
1143 : : ! GA: TODO:
1144 : : ! Eventually, m_ifc should be refactored to consume the Supercell object
1145 : : ! convention directly. Then, this step could be removed.
1146 : :
1147 : : ! Compute supercell rprimd scaled by acell
1148 [ + + + + ]: 598 : rprimd = SC%rprimd_scaled
1149 [ + + ]: 184 : do ii=1,3
1150 [ + + ]: 598 : rprimd(:,ii) = rprimd(:,ii) / SC%acell(ii)
1151 : : end do
1152 : :
1153 : : ! Ideal positions are computed on the fly (SC%get_xred_scaled(iatom,...)).
1154 : : ! They are anchored on atom 1's ideal position and wrapped to [-0.5,0.5).
1155 : : ! This anchoring matters: wrapping is not linear, so wrap(x-origin) differs
1156 : : ! from wrap(x)-origin whenever origin is not itself a lattice vector (i.e.
1157 : : ! whenever the first unitcell site is not at the reduced-coordinate origin).
1158 : 46 : origin = SC%get_xred_scaled(1)
1159 : :
1160 : : ! Scaled position of atom 1 unitcell site, placed in the origin cell
1161 : 46 : call SC%get_iat_unitcell(1,icell,iat_unit)
1162 : 46 : xred_unit = SC%get_xred(SC%get_iat_supercell(1,iat_unit))
1163 : 46 : R_cell_ref(:) = SC%unit_to_scaled(xred_unit)
1164 : :
1165 : : ! When the supercell does not repeat the unitcell along one reduced direction
1166 : : ! (multiplicity=1 along the diagonal, 0 elsewhere in that row),
1167 : : ! force the per-cell reduced-coordinate component of Rlatt_red(:,1,:)
1168 : : ! to zero along that direction for every atom.
1169 [ + + ]: 184 : do ii=1,3
1170 : : freeze_dir(ii) = (abs(Invar%multiplicity(ii,ii)-one)<tol8) .and. &
1171 : : & (abs(Invar%multiplicity(ii,mod(ii,3)+1))<tol8) .and. &
1172 [ + + + + : 321 : & (abs(Invar%multiplicity(ii,mod(ii+1,3)+1))<tol8)
- + ]
1173 : : end do
1174 : :
1175 : : ! Compute Rlatt
1176 [ + + ]: 7324 : do iatom=1,natom
1177 : :
1178 : : ! Ideal position of this atom, anchored on atom 1 and wrapped;
1179 : : ! reused below for every iatcell.
1180 [ + + ]: 29112 : xred_ideal_iatom = SC%wrap_inbox(SC%get_xred_scaled(iatom) - origin)
1181 : :
1182 : 7278 : call SC%get_iat_unitcell(iatom,icell,iat_unit)
1183 : 7278 : xred_unit = SC%get_xred(SC%get_iat_supercell(1,iat_unit))
1184 [ + + ]: 36390 : R_cell(:) = xred_ideal_iatom(:) + R_cell_ref(:) - SC%unit_to_scaled(xred_unit)
1185 [ + + + + ]: 29112 : where (freeze_dir) R_cell = zero
1186 : :
1187 [ + + ]: 23002 : do iatcell=1,natom_unitcell
1188 : :
1189 : 15678 : xred_unit = SC%get_xred(SC%get_iat_supercell(1,iatcell))
1190 [ + + ]: 78390 : xred_scaled = SC%unit_to_scaled(xred_unit) - origin
1191 : 15678 : call tdep_make_inbox(xred_scaled,1,tol8)
1192 : :
1193 [ + + ]: 62712 : xred_scaled = xred_ideal_iatom(:) - xred_scaled
1194 : 15678 : Rlatt = R_cell
1195 : 15678 : call tdep_make_inbox(xred_scaled,1,tol8,Rlatt)
1196 : :
1197 [ + + ]: 62712 : MD%Rlatt_cart(:,iatcell,iatom) = SC%red_scaled_to_cart(Rlatt)
1198 [ + + + + : 258126 : MD%Rlatt_scaled(:,iatcell,iatom) = matmul(rprimd, Rlatt)
+ + + - ]
1199 : :
1200 : : end do
1201 : :
1202 : : end do
1203 : :
1204 : : ! Write distance/Rlatt for debugging.
1205 : 46 : call tdep_write_distance_rlatt(Invar,MPIdata,SC,natom,natom_unitcell,MD%Rlatt_cart,MD%Rlatt_scaled)
1206 : :
1207 : 46 : end subroutine tdep_match_supercell
1208 : : !!***
1209 : :
1210 : : !====================================================================================================
1211 : :
1212 : : !!****f* ABINIT/m_tdep_sampling/get_ifc_order
1213 : : !! NAME
1214 : : !! get_ifc_order
1215 : : !!
1216 : : !! FUNCTION
1217 : : !! Enumerate the atoms of the supercell in the reference order expected by
1218 : : !! the legacy IFC routines (i.e. the order in which tdep_MatchIdeal2Average,
1219 : : !! itself inherited from m_dynmat/make_bigbox, enumerates atoms when building
1220 : : !! its "bigbox").
1221 : : !!
1222 : : !! This can be combined with SC%reorder to bring any per-atom (3,natom) array,
1223 : : !! currently ordered according to the natural SC order, into the legacy
1224 : : !! (IFC-reference) order.
1225 : : !!
1226 : : !! INPUTS
1227 : : !! SC = Supercell object describing the ideal structure.
1228 : : !! iat_ref (optional) = natural SC atom (see get_iat_supercell/get_iat_unitcell)
1229 : : !! that should be placed at reference position 1. The legacy routine
1230 : : !! (tdep_MatchIdeal2Average) does not enumerate atoms in SC's own natural
1231 : : !! frame: it anchors its bigbox scan on whichever atom of the unitcell#1
1232 : : !! pattern it finds first in the MD trajectory (its "iatom_ref"), shifting
1233 : : !! the whole enumeration so that atom sits at reference position 1. Passing
1234 : : !! the natural SC atom corresponding to that same physical atom here
1235 : : !! reproduces the legacy ordering exactly. If absent, the enumeration is
1236 : : !! anchored on whatever site the bigbox scan naturally finds first (SC's own
1237 : : !! frame), which will generally NOT match the legacy order.
1238 : : !!
1239 : : !! OUTPUT
1240 : : !! order(natom) = order(iatom_ref) = iat_super, where iat_super is the natural
1241 : : !! (icell,iat_unit) index of the atom (see get_iat_supercell/get_iat_unitcell)
1242 : : !! that occupies reference position iatom_ref in the legacy enumeration.
1243 : : !!
1244 : : !! SOURCE
1245 : :
1246 [ # # ]: 0 : function get_ifc_order(SC, iat_ref) result(order)
1247 : :
1248 : : type(tdep_Supercell_type),intent(in) :: SC
1249 : : integer, intent(in), optional :: iat_ref
1250 : : integer :: order(SC%natom)
1251 : :
1252 : : integer, parameter :: max_ijk=20
1253 : : integer :: ii,jj,kk,iatcell,iatom
1254 : : integer :: jcell,jat_supercell
1255 : : integer :: R_at(3)
1256 : : integer :: R_shift(3)
1257 : : logical :: discard_R
1258 : : logical :: first
1259 : : real(dp) :: Rlatt(3),xred_try(3),xred_box(3)
1260 : :
1261 [ # # ]: 0 : order(:) = 0
1262 : 0 : R_shift = 0
1263 : 0 : first = .true.
1264 : :
1265 : 0 : iatom = 1
1266 [ # # ]: 0 : do ii=-max_ijk,max_ijk
1267 [ # # ]: 0 : do jj=-max_ijk,max_ijk
1268 [ # # ]: 0 : do kk=-max_ijk,max_ijk
1269 : :
1270 : 0 : Rlatt(1) = real(ii-1,kind=dp)
1271 : 0 : Rlatt(2) = real(jj-1,kind=dp)
1272 : 0 : Rlatt(3) = real(kk-1,kind=dp)
1273 : :
1274 : 0 : discard_R = .false.
1275 [ # # ]: 0 : do iatcell=1,SC%natom_unitcell
1276 : :
1277 [ # # ]: 0 : if (discard_R) cycle
1278 : :
1279 : : ! Compute the reduced position in the supercell.
1280 [ # # ]: 0 : xred_try(:) = Rlatt(:) + SC%xred_unitcell(:,iatcell)
1281 : 0 : xred_box = SC%unit_to_scaled(xred_try)
1282 : :
1283 : : ! If the first atom of the pattern is in the [0;1[ range then keep
1284 : : ! all the atoms of the pattern (even if the others are outside the
1285 : : ! box). Else, none are taken.
1286 [ # # ]: 0 : if (iatcell==1) then
1287 [ # # # # : 0 : if (minval(xred_box(:)).lt.0.d0.or.maxval(xred_box(:)).ge.(1.d0-tol12)) then
# # # # #
# # # # #
# # ]
1288 : : discard_R = .true.
1289 : : cycle
1290 : : end if
1291 : : end if
1292 : :
1293 [ # # ]: 0 : if (iatom.gt.(SC%natom+1)) then
1294 : 0 : ABI_ERROR('The number of atoms found in the bigbox exceeds natom')
1295 : : end if
1296 : :
1297 : : ! The very first atom kept by the scan defines reference position 1.
1298 : : ! Compute the shift that brings it into coincidence with iat_ref
1299 : : ! (if provided), and apply it to every subsequent site so the whole
1300 : : ! enumeration is anchored on iat_ref, exactly as the legacy routine
1301 : : ! anchors its own bigbox scan on iatom_ref.
1302 [ # # ]: 0 : if (first) then
1303 [ # # ]: 0 : if (present(iat_ref)) then
1304 [ # # ]: 0 : R_shift = NINT(Rlatt) - SC%get_R_at(iat_ref)
1305 : : end if
1306 : : first = .false.
1307 : : end if
1308 : :
1309 [ # # ]: 0 : R_at = NINT(Rlatt) - R_shift
1310 : 0 : jcell = SC%get_icell(R_at)
1311 : 0 : jat_supercell = SC%get_iat_supercell(jcell, iatcell)
1312 : :
1313 : 0 : order(iatom) = jat_supercell
1314 : 0 : iatom = iatom + 1
1315 : :
1316 : : end do
1317 : : end do
1318 : : end do
1319 : : end do
1320 : :
1321 [ # # ]: 0 : if (iatom.lt.SC%natom+1) then
1322 : 0 : ABI_ERROR('The number of atoms found in the big box is smaller than natom')
1323 : : end if
1324 : :
1325 : 0 : end function get_ifc_order
1326 : : !!***
1327 : :
1328 : : !====================================================================================================
1329 : :
1330 : : !!****f* ABINIT/m_tdep_sampling/tdep_write_xred_average_ideal
1331 : : !! NAME
1332 : : !! tdep_write_xred_average_ideal
1333 : : !!
1334 : : !! FUNCTION
1335 : : !! Write, for each atom, the ideal (equilibrium) position and the average
1336 : : !! MD position, both in reduced and cartesian coordinates, in the
1337 : : !! _xred_average.xyz file, for diagnostic purposes.
1338 : : !!
1339 : : !! INPUTS
1340 : : !! Invar = Input object containing the input variables.
1341 : : !! xred_average(3,natom) = Average reduced positions of the atoms.
1342 : : !! xred_ideal(3,natom) = Ideal reduced positions of the atoms.
1343 : : !! original_order(natom) = Index of each atom in the original ordering
1344 : : !! of the MD trajectory (written next to the
1345 : : !! average position).
1346 : : !!
1347 : : !! OUTPUT
1348 : : !! Only writing.
1349 : : !!
1350 : : !! SOURCE
1351 : :
1352 : 46 : subroutine tdep_write_xred_average_ideal(Invar,xred_average,xred_ideal,original_order)
1353 : :
1354 : : type(atdep_dataset_type), intent(in) :: Invar
1355 : : real(dp),intent(in) :: xred_average(3,Invar%natom)
1356 : : real(dp),intent(in) :: xred_ideal(3,Invar%natom)
1357 : : integer,intent(in) :: original_order(Invar%natom)
1358 : :
1359 : : integer :: unt
1360 : : integer :: iatom,jatom,ii,jj
1361 : : real(dp) :: rprimdt(3,3)
1362 : : real(dp) :: xred_A(3),xred_I(3),xcart_A(3),xcart_I(3)
1363 : :
1364 [ + + + + ]: 598 : rprimdt(:,:) = transpose(Invar%rprimd_md(:,:))
1365 : :
1366 : 46 : unt = get_unit()
1367 : 46 : open(unit=unt,file=trim(Invar%output_prefix)//'_xred_average.xyz')
1368 : 46 : write(unt,'(a,i4)') '# natom = ',Invar%natom
1369 : 46 : write(unt,'(a,i4)') '# natom_unitcell = ',Invar%natom_unitcell
1370 [ + - + + : 598 : write(unt,'(a,9(f4.1,1x))') '# multiplicity = ',((Invar%multiplicity(ii,jj),jj=1,3),ii=1,3 )
+ - + + ]
1371 : 46 : write(unt,'(a)') '#'
1372 : :
1373 : 46 : write(unt,'(a1,1x,a8,4x,a6,4x,3(a5,30x))') '#', 'position', 'iatom', 'xred ', 'xcart'
1374 : 46 : write(unt,'(a)')''
1375 : :
1376 : 46 : xred_I = zero
1377 : 46 : xred_A = zero
1378 [ + + ]: 7324 : do iatom=1,Invar%natom
1379 [ + + ]: 29112 : xred_A(:) = xred_average(:,iatom)
1380 [ + + ]: 29112 : xred_I(:) = xred_ideal(:,iatom)
1381 : :
1382 [ + + + + ]: 94614 : xcart_A = matmul(rprimdt, xred_A)
1383 [ + + + + ]: 94614 : xcart_I = matmul(rprimdt, xred_I)
1384 : :
1385 : 7278 : jatom = original_order(iatom)
1386 : :
1387 : : ! write(unt,'(2x,a8,4x,i6,2x,3(f10.6,1x),2x,3(f10.6,1x))')'Average',jatom,xred_A,xcart_A
1388 : : ! write(unt,'(2x,a8,4x,i6,2x,3(f10.6,1x),2x,3(f10.6,1x),2x,3(f10.6,1x))')'Ideal ',iatom,xred_I,xcart_I,R_red
1389 : 7278 : write(unt,'(2x,a8,4x,i6,2x,3(f10.6,1x),2x,3(f10.6,1x),2x,3(f10.6,1x))')'Ideal ',iatom,xred_I,xcart_I
1390 : 7278 : write(unt,'(2x,a8,4x,i6,2x,3(f10.6,1x),2x,3(f10.6,1x))')'Average',jatom,xred_A,xcart_A
1391 : 7324 : write(unt,'(a)')''
1392 : :
1393 : : end do
1394 : :
1395 : 46 : close(unt)
1396 : :
1397 : 46 : end subroutine tdep_write_xred_average_ideal
1398 : : !!***
1399 : :
1400 : : !====================================================================================================
1401 : :
1402 : : !!****f* ABINIT/m_tdep_sampling/tdep_write_ucart_fcart
1403 : : !! NAME
1404 : : !! tdep_write_ucart_fcart
1405 : : !!
1406 : : !! FUNCTION
1407 : : !! Write, for each MD step and each atom, the Cartesian displacement
1408 : : !! (ucart) and the Cartesian force (fcart), together with the relevant
1409 : : !! indices. This is meant for debugging/comparison purposes, e.g. to
1410 : : !! directly diff the output of tdep_match_supercell against that of
1411 : : !! tdep_MatchIdeal2Average for the same input trajectory.
1412 : : !!
1413 : : !! INPUTS
1414 : : !! Invar = Input object containing the input variables.
1415 : : !! MPIdata = Info on MPI parallelism.
1416 : : !! natom = Number of atoms in the supercell.
1417 : : !! nstep = Number of MD steps (locally held).
1418 : : !! typat(natom) = Atom type, in the same order as ucart/fcart.
1419 : : !! ucart(3,natom,nstep) = Cartesian displacements.
1420 : : !! fcart(3,natom,nstep) = Cartesian forces.
1421 : : !!
1422 : : !! NOTES
1423 : : !! ucart and fcart must already be expressed in the same atom ordering
1424 : : !! (this is not necessarily true of MD%fcart itself in every routine;
1425 : : !! callers should pass whichever arrays are mutually consistent).
1426 : : !!
1427 : : !! SOURCE
1428 : :
1429 : 0 : subroutine tdep_write_ucart_fcart(Invar,MPIdata,natom,nstep,typat,ucart,fcart)
1430 : :
1431 : : type(atdep_dataset_type),intent(in) :: Invar
1432 : : type(MPI_enreg_type),intent(in) :: MPIdata
1433 : : integer,intent(in) :: natom,nstep
1434 : : integer,intent(in) :: typat(natom)
1435 : : real(dp),intent(in) :: ucart(3,natom,nstep)
1436 : : real(dp),intent(in) :: fcart(3,natom,nstep)
1437 : :
1438 : : integer :: unt,istep,iatom
1439 : :
1440 [ # # ]: 0 : if (.not.MPIdata%iam_master) return
1441 : :
1442 : 0 : unt = get_unit()
1443 : 0 : open(unit=unt,file=trim(Invar%output_prefix)//'_ucart_fcart.dat')
1444 : :
1445 : 0 : write(unt,'(a)') '# istep iatom typat ucart(1:3) fcart(1:3)'
1446 : :
1447 [ # # ]: 0 : do istep=1,nstep
1448 [ # # ]: 0 : do iatom=1,natom
1449 : : write(unt,'(2(i6,1x),i4,4x,3(f16.10,1x),4x,3(f16.10,1x))') &
1450 : 0 : & istep,iatom,typat(iatom),ucart(:,iatom,istep),fcart(:,iatom,istep)
1451 : : end do
1452 : : end do
1453 : :
1454 : 0 : close(unt)
1455 : :
1456 : : end subroutine tdep_write_ucart_fcart
1457 : : !!***
1458 : :
1459 : : !====================================================================================================
1460 : :
1461 : : !!****f* ABINIT/m_tdep_sampling/tdep_write_distance_rlatt
1462 : : !! NAME
1463 : : !! tdep_write_distance_rlatt
1464 : : !!
1465 : : !! FUNCTION
1466 : : !! Write distance (interatomic distances between ideal positions) and
1467 : : !! MD%Rlatt_cart/MD%Rlatt_scaled (lattice vectors between the ideal site of
1468 : : !! each atom and each unitcell reference atom), for debugging/comparison
1469 : : !! purposes, e.g. to directly diff the output of tdep_match_supercell
1470 : : !! against that of tdep_MatchIdeal2Average for the same input trajectory.
1471 : : !! Both quantities are purely geometric (computed from ideal positions
1472 : : !! only), so unlike ucart/fcart they do not depend on which raw MD atom
1473 : : !! was matched to which site, and should already share a common atom
1474 : : !! ordering (the IFC-reference/box-scan order) in both routines.
1475 : : !!
1476 : : !! INPUTS
1477 : : !! Invar = Input object containing the input variables.
1478 : : !! MPIdata = Info on MPI parallelism.
1479 : : !! SC = Supercell object describing the ideal structure. SC%dij (distance
1480 : : !! and cartesian components between the ideal positions of pairs of
1481 : : !! atoms) is used.
1482 : : !! natom = Number of atoms in the supercell.
1483 : : !! natom_unitcell = Number of atoms in the unitcell.
1484 : : !! Rlatt_cart(3,natom_unitcell,natom) = Cartesian lattice vector between
1485 : : !! a unitcell reference atom and the ideal
1486 : : !! site of each atom.
1487 : : !! Rlatt_scaled(3,natom_unitcell,natom) = Rlatt_cart divided by acell_unitcell.
1488 : : !!
1489 : : !! SOURCE
1490 : :
1491 : 46 : subroutine tdep_write_distance_rlatt(Invar,MPIdata,SC,natom,natom_unitcell,Rlatt_cart,Rlatt_scaled)
1492 : :
1493 : : type(atdep_dataset_type),intent(in) :: Invar
1494 : : type(MPI_enreg_type),intent(in) :: MPIdata
1495 : : type(tdep_Supercell_type),intent(in) :: SC
1496 : : integer,intent(in) :: natom,natom_unitcell
1497 : : real(dp),intent(in) :: Rlatt_cart(3,natom_unitcell,natom)
1498 : : real(dp),intent(in) :: Rlatt_scaled(3,natom_unitcell,natom)
1499 : :
1500 : : integer :: unt,iatom,jatom,iatcell
1501 : :
1502 [ - + ]: 46 : if (.not.MPIdata%iam_master) return
1503 : :
1504 : 46 : unt = get_unit()
1505 : 46 : open(unit=unt,file=trim(Invar%output_prefix)//'_distance.dat')
1506 : 46 : write(unt,'(a)') '# iatom jatom dist rij(1:3)'
1507 [ + + ]: 7324 : do iatom=1,natom
1508 [ + + ]: 1558688 : do jatom=1,natom
1509 : 1558642 : write(unt,'(2(i6,1x),4(f16.10,1x))') iatom,jatom,SC%dij(iatom,jatom,1),SC%dij(iatom,jatom,2:4)
1510 : : end do
1511 : : end do
1512 : 46 : close(unt)
1513 : :
1514 : 46 : unt = get_unit()
1515 : 46 : open(unit=unt,file=trim(Invar%output_prefix)//'_rlatt.dat')
1516 : 46 : write(unt,'(a)') '# iatcell iatom Rlatt_cart(1:3) Rlatt_scaled(1:3)'
1517 [ + + ]: 7324 : do iatom=1,natom
1518 [ + + ]: 23002 : do iatcell=1,natom_unitcell
1519 : : write(unt,'(2(i6,1x),3(f16.10,1x),4x,3(f16.10,1x))') &
1520 : 22956 : & iatcell,iatom,Rlatt_cart(:,iatcell,iatom),Rlatt_scaled(:,iatcell,iatom)
1521 : : end do
1522 : : end do
1523 : 46 : close(unt)
1524 : :
1525 : : end subroutine tdep_write_distance_rlatt
1526 : : !!***
1527 : :
1528 : : !====================================================================================================
1529 : :
1530 [ # # # # : 147072 : end module m_tdep_sampling
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1531 : : !!***
|
