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_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
32 : use m_tdep_latt, only : Lattice_type, tdep_make_inbox
33 : use m_tdep_sym, only : Symmetries_type, tdep_SearchS_1at
34 :
35 : implicit none
36 :
37 : type tdep_Sampling_type
38 :
39 : integer :: natom
40 : ! Number of atoms in the supercell
41 :
42 : integer :: natom_unitcell
43 : ! Number of atoms in the unitcell
44 :
45 : integer :: my_nstep
46 : ! Number of MD steps held locally
47 :
48 : integer :: nstep_tot
49 : ! Total of MD steps
50 :
51 : integer, allocatable :: typat_unitcell(:)
52 : ! typat_unitcell(natom_unitcell)
53 : ! Atom type in the unitcell.
54 :
55 : integer, allocatable :: typat(:)
56 : ! typat(natom)
57 : ! Atom type in the supercell.
58 :
59 : double precision, allocatable :: xred_unitcell(:,:)
60 : ! xred_unitcell(3, natom_unitcell)
61 : ! Reduced equilibrium position of the atoms in the unitcell.
62 :
63 : double precision, allocatable :: xred_ideal(:,:)
64 : ! xred_ideal(3, natom)
65 : ! Reduced equilibrium position of the atoms in the supercell.
66 :
67 : double precision, allocatable :: xred(:,:,:)
68 : ! xred(3, natom, my_nstep)
69 : ! Reduced position of the atom at each step.
70 :
71 : double precision, allocatable :: fcart(:,:,:)
72 : ! fcart(3, natom, my_nstep)
73 : ! Cartesian forces at each step.
74 :
75 : double precision, allocatable :: etot(:)
76 : ! etot(my_nstep)
77 : ! Total energy at each step.
78 :
79 : double precision, allocatable :: weights(:)
80 : ! weights(my_nstep)
81 : ! The weight of each configuration for the fitting.
82 : ! By default, these will be 1 / nstep_tot.
83 :
84 : double precision, allocatable :: Rlatt_cart(:,:,:)
85 : ! Rlatt_cart(3, natom_unitcell, natom)
86 : ! Cartesian coordinate of the lattice vectors of the unitcell
87 : ! within the supercell, for each atom.
88 :
89 : double precision, allocatable :: Rlatt_scaled(:,:,:)
90 : ! Rlatt_scaled(3, natom_unitcell, natom)
91 : ! Rlatt_cart divided by acell_unitcell.
92 :
93 : double precision, allocatable :: ucart(:,:,:)
94 : ! ucart(3, natom, my_nstep)
95 : ! Cartesian displacements of the atoms with respect to their equilibrium
96 : ! positions at each time step.
97 :
98 : double precision, allocatable :: distance(:,:,:)
99 : ! distance(natom, natom, 4)
100 : ! Distance between the ideal positions of the atoms in the supercell,
101 : ! (norm, and cartesian components).
102 :
103 : double precision, allocatable :: Forces(:)
104 : ! Forces(3*natom*my_nstep)
105 : ! The cartesian forces for all configurations, as a flat array.
106 :
107 : end type tdep_Sampling_type
108 :
109 : public :: tdep_sampling_init_read
110 : public :: tdep_sampling_free
111 : public :: tdep_sampling_shift_xred
112 : public :: tdep_sampling_rotate
113 : public :: tdep_MatchIdeal2Average
114 : public :: tdep_write_xred_average
115 :
116 : contains
117 :
118 : !=====================================================================================================
119 :
120 44 : subroutine tdep_sampling_init_read(MD,Invar,MPIdata,Hist)
121 :
122 : type(tdep_Sampling_type), intent(inout) :: MD
123 : type(atdep_dataset_type), intent(in) :: Invar
124 : type(MPI_enreg_type), intent(in) :: MPIdata
125 : type(abihist), intent(in) :: Hist
126 :
127 : integer :: this_istep,istep,iatom,jstep
128 : double precision :: tmp1,tmp2,tmp3
129 :
130 44 : MD%natom = Invar%natom
131 44 : MD%natom_unitcell = Invar%natom_unitcell
132 44 : MD%my_nstep = Invar%my_nstep
133 44 : MD%nstep_tot = Invar%nstep_tot
134 :
135 238 : ABI_CALLOC(MD%typat_unitcell, (MD%natom_unitcell))
136 7014 : ABI_CALLOC(MD%typat, (MD%natom))
137 556 : ABI_CALLOC(MD%xred_unitcell, (3,MD%natom_unitcell))
138 27660 : ABI_CALLOC(MD%xred_ideal, (3,MD%natom))
139 557736 : ABI_CALLOC(MD%xred, (3,MD%natom,MD%my_nstep))
140 557692 : ABI_CALLOC(MD%fcart, (3,MD%natom,MD%my_nstep))
141 1084 : ABI_CALLOC(MD%etot, (MD%my_nstep))
142 1040 : ABI_CALLOC(MD%weights, (MD%my_nstep))
143 5917156 : ABI_CALLOC(MD%distance, (MD%natom,MD%natom,4))
144 66602 : ABI_CALLOC(MD%Rlatt_scaled, (3,MD%natom_unitcell,MD%natom))
145 66558 : ABI_CALLOC(MD%Rlatt_cart, (3,MD%natom_unitcell,MD%natom))
146 557692 : ABI_CALLOC(MD%ucart, (3,MD%natom,MD%my_nstep))
147 417588 : ABI_CALLOC(MD%Forces, (3*MD%natom*MD%my_nstep))
148 :
149 150 : MD%typat_unitcell(:) = Invar%typat_unitcell(:)
150 468 : MD%xred_unitcell(:,:) = Invar%xred_unitcell(:,:)
151 6926 : MD%typat(:) = Invar%typat(:)
152 :
153 : ! Read xred.dat, fcart.dat and etot.dat ASCII files or extract them from the HIST.nc netcdf file.
154 44 : write(Invar%stdout,'(a)') ' '
155 44 : this_istep=0
156 44 : jstep=0
157 44 : if (Invar%use_weights.eq.1) then
158 1 : open(unit=30,file=trim(Invar%input_prefix)//'_weights.dat')
159 43 : else if (Invar%use_weights.eq.0) then
160 975 : MD%weights=1.0d0/real(MD%nstep_tot)
161 : endif
162 44 : if (Invar%netcdf) then
163 217 : do istep=Invar%nstep_min,Invar%nstep_max
164 217 : if (mod(istep-Invar%nstep_min,Invar%slice).ne.0) then
165 : cycle
166 : else
167 201 : jstep=jstep+1
168 201 : if (.not.MPIdata%my_step(jstep)) cycle
169 201 : this_istep=this_istep+1
170 62441 : MD%xred(:,:,this_istep) =Hist%xred (:,:,istep)
171 62441 : MD%fcart(:,:,this_istep)=Hist%fcart(:,:,istep)
172 211 : MD%etot(this_istep) =Hist%etot (istep)
173 : end if
174 : end do !istep
175 6 : write(Invar%stdout,'(a)') ' The positions, forces and energies are extracted from the NetCDF file: HIST.nc'
176 : else
177 38 : open(unit=60,file=trim(Invar%input_prefix)//'_fcart.dat')
178 38 : open(unit=50,file=trim(Invar%input_prefix)//'_xred.dat')
179 38 : open(unit=40,file=trim(Invar%input_prefix)//'_etot.dat')
180 38 : do istep=1,Invar%nstep_min-1
181 0 : if (Invar%use_weights.eq.1) then
182 0 : read(30,*) tmp1
183 : endif
184 0 : read(40,*) tmp1
185 38 : do iatom=1,MD%natom
186 0 : read(50,*) tmp1,tmp2,tmp3
187 0 : read(60,*) tmp1,tmp2,tmp3
188 : end do
189 : end do
190 799 : do istep=Invar%nstep_min,Invar%nstep_max
191 799 : if (mod(istep-Invar%nstep_min,Invar%slice).ne.0) then
192 10 : if (Invar%use_weights.eq.1) then
193 0 : read(30,*) tmp1
194 : endif
195 10 : read(40,*) tmp1
196 1290 : do iatom=1,Invar%natom
197 1280 : read(50,*) tmp1,tmp2,tmp3
198 1290 : read(60,*) tmp1,tmp2,tmp3
199 : end do
200 : else
201 751 : jstep=jstep+1
202 751 : if (.not.MPIdata%my_step(jstep)) then
203 0 : if (Invar%use_weights.eq.1) then
204 0 : read(30,*) tmp1
205 : endif
206 0 : read(40,*) tmp1
207 0 : do iatom=1,Invar%natom
208 0 : read(50,*) tmp1,tmp2,tmp3
209 0 : read(60,*) tmp1,tmp2,tmp3
210 : end do
211 : else
212 751 : this_istep=this_istep+1
213 751 : if (Invar%use_weights.eq.1) then
214 20 : read(30,*) MD%weights(this_istep)
215 : endif
216 751 : read(40,*) MD%etot(this_istep)
217 124343 : do iatom=1,MD%natom
218 123592 : read(50,*) MD%xred (1,iatom,this_istep),MD%xred (2,iatom,this_istep),MD%xred (3,iatom,this_istep)
219 124343 : read(60,*) MD%fcart(1,iatom,this_istep),MD%fcart(2,iatom,this_istep),MD%fcart(3,iatom,this_istep)
220 : end do
221 : end if !my_step
222 : end if !slice
223 : end do !istep
224 38 : close(40)
225 38 : close(50)
226 38 : close(60)
227 38 : write(Invar%stdout,'(2a)') ' The positions, forces and energies are extracted from the ASCII files:',&
228 76 : & ' xred.dat, fcart.dat & etot.dat'
229 : end if !netcdf
230 44 : if (Invar%use_weights.eq.1) then
231 1 : close(30)
232 : end if
233 :
234 44 : end subroutine tdep_sampling_init_read
235 :
236 : !=====================================================================================================
237 :
238 44 : subroutine tdep_sampling_free(MD)
239 :
240 : type(tdep_Sampling_type), intent(inout) :: MD
241 :
242 44 : ABI_FREE(MD%typat_unitcell)
243 44 : ABI_FREE(MD%typat)
244 44 : ABI_FREE(MD%xred_unitcell)
245 44 : ABI_FREE(MD%xred_ideal)
246 44 : ABI_FREE(MD%xred)
247 44 : ABI_FREE(MD%fcart)
248 44 : ABI_FREE(MD%etot)
249 44 : ABI_FREE(MD%weights)
250 44 : ABI_FREE(MD%distance)
251 44 : ABI_FREE(MD%Rlatt_scaled)
252 44 : ABI_FREE(MD%Rlatt_cart)
253 44 : ABI_FREE(MD%ucart)
254 44 : ABI_FREE(MD%Forces)
255 :
256 44 : end subroutine tdep_sampling_free
257 :
258 : !=====================================================================================================
259 :
260 : ! Shift xred to keep atoms in the same unit cell at each step.
261 44 : subroutine tdep_sampling_shift_xred(MD,MPIdata)
262 :
263 : type(tdep_Sampling_type), intent(inout) :: MD
264 : type(MPI_enreg_type), intent(in) :: MPIdata
265 : integer :: natom,ii,iatom,istep,ierr
266 : integer :: shift,shift_max,shift_best
267 : double precision :: xi, dist, best_dist
268 44 : double precision, allocatable :: x0(:,:)
269 :
270 44 : natom = MD%natom
271 132 : ABI_MALLOC(x0,(3,natom))
272 :
273 : ! Communicate xred at the first step
274 27572 : x0(:,:) = zero
275 44 : if (MPIdata%my_step(1)) then
276 27572 : x0(:,:) = MD%xred(:,:,1)
277 : end if
278 44 : call xmpi_sum(x0,MPIdata%comm_step,ierr)
279 :
280 : ! Shift xred from all steps in the same unitcell as the first step
281 44 : shift_max = 1
282 996 : do istep=1, MD%my_nstep
283 140148 : do iatom=1,natom
284 557560 : do ii=1,3
285 417456 : best_dist = abs(MD%xred(ii,iatom,istep) - x0(ii,iatom))
286 417456 : shift_best = 0
287 1669824 : do shift=-shift_max,shift_max
288 1252368 : xi = MD%xred(ii,iatom,istep) + shift
289 1252368 : dist = abs(xi - x0(ii,iatom))
290 1669824 : if (dist < best_dist) then
291 366 : best_dist = dist
292 366 : shift_best = shift
293 : end if
294 : end do
295 556608 : MD%xred(ii,iatom,istep) = MD%xred(ii,iatom,istep) + shift_best
296 : end do
297 : end do
298 : end do
299 :
300 44 : ABI_FREE(x0)
301 :
302 44 : end subroutine tdep_sampling_shift_xred
303 :
304 : !=====================================================================================================
305 :
306 44 : subroutine tdep_sampling_rotate(MD,rotation_cart)
307 :
308 : type(tdep_Sampling_type), intent(inout) :: MD
309 : integer :: iatom,istep
310 : double precision :: rotation_cart(3,3)
311 :
312 : ! Apply rotation to fcart
313 996 : do istep=1,MD%my_nstep
314 140148 : do iatom=1,MD%natom
315 2783992 : MD%fcart(:,iatom,istep) = MATMUL(rotation_cart, MD%fcart(:,iatom,istep))
316 : end do
317 : end do
318 :
319 44 : end subroutine tdep_sampling_rotate
320 :
321 : !=====================================================================================================
322 :
323 : !!***
324 : !!****f* ABINIT/m_tdep_sampling/tdep_MatchIdeal2Average
325 : !! NAME
326 : !! tdep_MatchIdeal2Average
327 : !!
328 : !! FUNCTION
329 : !! Find the mapping between the atoms in the ideal (equilibrium) supercell,
330 : !! and the atoms of the input moledular dynamics using their average positions.
331 : !! Then compute the atom displacements with respect to the equilibrium positions
332 : !! at each time step of the MD.
333 : !!
334 : !! INPUTS
335 : !! MD = TDEP Sampling object containing the input positions and forces.
336 : !! Invar = Input object containing the input variables.
337 : !! Lattice = Lattice object describing the ideal structure.
338 : !! Sym = Symetries object describing all the symmetry operations of the crystal.
339 : !! MPIdata = Info on MPI parallelism.
340 : !!
341 : !! OUTPUT
342 : !!
343 : !! SIDE EFFECTS
344 : !! The following quantities in MD are computed:
345 : !!
346 : !! distance = Distance between the ideal positions of the atoms in the supercell,
347 : !! (norm, and cartesian components).
348 : !! Forces = Cartesian forces on the atoms at each time steps, as a flat array.
349 : !! ucart = Cartesian displacements of the atoms with respect to their equilibrium
350 : !! positions at each time step.
351 : !! Rlatt_cart = Cartesian coordinate of the lattice vectors of the unitcell
352 : !! within the supercell, for each atom.
353 : !! This array seems to have an extra dimension, for algorithmic simplicity.
354 : !! Rlatt_scaled = Rlatt_cart divided by acell_unitcell.
355 : !! These are used when reading an IFC file, to compare with the R vectors
356 : !! that are stored in the file.
357 : !!
358 : !! Some of the reduced positions of the atoms MD%xred are shifted by a supercell
359 : !! lattice vector in order to re-center the crystal.
360 : !!
361 : !! NOTES
362 : !!
363 : !! SOURCE
364 :
365 44 : subroutine tdep_MatchIdeal2Average(MD,Invar,Lattice,Sym,MPIdata)
366 :
367 : type(tdep_Sampling_type),intent(inout) :: MD
368 : type(atdep_dataset_type),intent(inout) :: Invar
369 : type(Lattice_type),intent(in) :: Lattice
370 : type(Symmetries_type),intent(inout) :: Sym
371 : type(MPI_enreg_type),intent(in) :: MPIdata
372 :
373 : integer :: ii,jj,kk,max_ijk,iatcell,jatcell,iatom,jatom,eatom,fatom,istep
374 : integer :: iatom_ref,ierr
375 : integer :: ndir_match,natom_match
376 : double precision :: tmp(3),tmp1(3),tmp2(3),Rlatt(3),xred_tmp(3),rprimd_md_tmp(3,3),distance_tmp(3)
377 44 : double precision, allocatable :: dist_unitcell(:,:,:),xcart_average(:,:)
378 44 : double precision, allocatable :: fcart_tmp(:,:,:),ucart_tmp(:,:,:)
379 44 : double precision, allocatable :: xred_average(:,:)
380 44 : double precision, allocatable :: xred_center(:,:)
381 44 : double precision, allocatable :: Rlatt_red (:,:,:)
382 44 : double precision, allocatable :: xred_ideal(:,:)
383 : ! double precision, allocatable :: distance_average(:,:,:)
384 44 : integer, allocatable :: FromIdeal2Average(:)
385 44 : double precision, allocatable :: xcart(:,:,:)
386 44 : double precision, allocatable :: xcart_ideal(:,:)
387 : logical :: ok,must_shift,discard_R
388 : character(len=500) :: msg
389 :
390 44 : write(Invar%stdout,*)' '
391 44 : write(Invar%stdout,*) '#############################################################################'
392 44 : write(Invar%stdout,*) '###### Find the matching between ideal and average positions ###############'
393 44 : write(Invar%stdout,*) '#############################################################################'
394 :
395 : !==========================================================================================
396 : !======== 1/ Determine ideal positions and distances ======================================
397 : !==========================================================================================
398 44 : write(Invar%stdout,*)' Determine ideal positions and distances...'
399 : ! Define the bigbox with ideal positions
400 66602 : ABI_CALLOC(Rlatt_red ,(3,MD%natom_unitcell,MD%natom))
401 27660 : ABI_CALLOC(xred_ideal,(3,MD%natom))
402 44 : max_ijk=20
403 44 : iatom=1
404 1848 : do ii=-max_ijk,max_ijk
405 75812 : do jj=-max_ijk,max_ijk
406 3108292 : do kk=-max_ijk,max_ijk
407 :
408 3032524 : Rlatt(1)=real(ii-1)
409 3032524 : Rlatt(2)=real(jj-1)
410 3032524 : Rlatt(3)=real(kk-1)
411 :
412 3032524 : discard_R = .false.
413 10412114 : do iatcell=1,MD%natom_unitcell
414 :
415 7305626 : if (discard_R) cycle
416 :
417 : ! Compute the reduced positions
418 12141444 : tmp(:) = Rlatt(:) + MD%xred_unitcell(:,iatcell)
419 3035361 : call DGEMV('T',3,3,1.d0,Lattice%multiplicitym1(:,:),3,tmp(:),1,0.d0,xred_tmp(:),1)
420 :
421 : ! If the first atom of the pattern is in the [0;1[ range then keep all the
422 : ! atoms of the pattern (even if the others are outside the box). Else,
423 : ! none are taken.
424 3035361 : if (iatcell==1) then
425 27292716 : if (minval(xred_tmp(:)).lt.0.d0.or.maxval(xred_tmp(:)).ge.(1.d0-tol12)) then
426 : discard_R = .true.
427 : cycle
428 : end if
429 : end if
430 :
431 : !GA: Why natom+1 ?
432 6882 : if (iatom.gt.(MD%natom+1)) then
433 0 : ABI_ERROR('The number of atoms found in the bigbox exceeds natom' )
434 : end if
435 :
436 27528 : xred_ideal(:,iatom) = xred_tmp(:)
437 6882 : call DGEMV('T',3,3,1.d0,Lattice%multiplicitym1(:,:),3,Rlatt(:),1,0.d0,Rlatt_red(:,1,iatom),1)
438 7309671 : iatom = iatom + 1
439 : end do
440 : end do
441 : end do
442 : end do
443 :
444 44 : if (iatom.lt.MD%natom+1) then
445 0 : ABI_ERROR('The number of atoms found in the big box is smaller than natom')
446 : end if
447 :
448 : ! Compute the distances between ideal positions in the SUPERcell
449 6926 : do eatom=1,MD%natom
450 1479234 : do fatom=1,MD%natom
451 5889232 : tmp(:)=xred_ideal(:,fatom)-xred_ideal(:,eatom)
452 1472308 : call tdep_make_inbox(tmp,1,1d-4)
453 19140004 : rprimd_md_tmp(:,:)=Lattice%rprimd_md(:,:)
454 5889232 : distance_tmp(:)=MD%distance(eatom,fatom,2:4)
455 1472308 : call DGEMV('T',3,3,1.d0,rprimd_md_tmp,3,tmp,1,0.d0,distance_tmp,1)
456 5889232 : MD%distance(eatom,fatom,2:4)=distance_tmp(:)
457 5889232 : do ii=1,3
458 : ! Remove the rounding errors before writing (for non regression testing purposes)
459 4416924 : if (abs(MD%distance(eatom,fatom,ii+1)).lt.tol8) MD%distance(eatom,fatom,ii+1)=zero
460 5889232 : MD%distance(eatom,fatom,1)=MD%distance(eatom,fatom,1)+(MD%distance(eatom,fatom,ii+1))**2
461 : end do
462 1472308 : MD%distance(eatom,fatom,1)=MD%distance(eatom,fatom,1)**0.5
463 1479190 : MD%distance(eatom,fatom,1)=tol12 * dint(MD%distance(eatom,fatom,1) / tol12)
464 : end do
465 : end do
466 :
467 : ! Compute the distances between ideal positions in the UNITcell
468 2212 : ABI_MALLOC(dist_unitcell,(MD%natom_unitcell,MD%natom_unitcell,3)); dist_unitcell(:,:,:)=zero
469 150 : do iatcell=1,MD%natom_unitcell
470 664 : do jatcell=1,MD%natom_unitcell
471 2056 : tmp(:) = xred_ideal(:,jatcell)-xred_ideal(:,iatcell)
472 514 : call tdep_make_inbox(tmp,1,tol8)
473 2162 : dist_unitcell(iatcell,jatcell,:) = tmp(:)
474 : end do
475 : end do
476 :
477 : !==========================================================================================
478 : !======== 2/ Find the matching between the ideal and average ==============================
479 : !======== (from the MD simulations) positions. ==========================================
480 : !==========================================================================================
481 : ! NOTE: - xred_center is used to find the matching with the ideal positions
482 : ! - xred_average is used to compute the displacements (from MD trajectories)
483 : ! The difference between those two is that xred_center will be shifted to bring
484 : ! one of the average positions at the origin, for an easier comparison with
485 : ! xred_ideal. Some shifts by a supercell lattice vector will be computed
486 : ! from the difference between xred_center and xred_ideal, and those shifts
487 : ! will be applied to xred_average and xred at all steps.
488 :
489 44 : write(Invar%stdout,*)' Compute average positions...'
490 27660 : ABI_CALLOC(xred_average,(3,MD%natom))
491 27616 : ABI_CALLOC(xred_center,(3,MD%natom))
492 : ! Average positions from MD (on nstep steps)
493 996 : do istep=1,MD%my_nstep
494 140148 : do iatom=1,MD%natom
495 557560 : xred_average(:,iatom)=xred_average(:,iatom)+MD%xred(:,iatom,istep)
496 : end do
497 : end do
498 44 : call xmpi_sum(xred_average,MPIdata%comm_step,ierr)
499 27572 : xred_average(:,:) = xred_average(:,:) / real(MD%nstep_tot)
500 :
501 : ! Search the basis of atoms in the supercell
502 : ! in order to find iatom_ref
503 44 : write(Invar%stdout,*)' Search the unitcell basis of atoms in the MD trajectory...'
504 44 : ok=.false.
505 27572 : xred_center(:,:) = xred_average(:,:)
506 44 : iatcell=1
507 51 : do iatom=1,MD%natom
508 51 : if (MD%typat(iatom).ne.MD%typat_unitcell(iatcell)) cycle
509 51 : natom_match = 0
510 7797 : do jatom=1,MD%natom
511 :
512 30984 : tmp(:)=xred_center(:,jatom)-xred_center(:,iatom)
513 7746 : call tdep_make_inbox(tmp,1,Invar%tolinbox)
514 :
515 24087 : do jatcell=1,MD%natom_unitcell
516 16403 : if (MD%typat(jatom).ne.MD%typat_unitcell(jatcell)) cycle
517 14263 : ndir_match = 0
518 57052 : do ii=1,3
519 57052 : if (abs(tmp(ii)-dist_unitcell(iatcell,jatcell,ii)).le.Invar%tolmotif) then
520 5723 : ndir_match=ndir_match+1
521 : end if
522 : end do
523 21896 : if (ndir_match==3) then
524 113 : natom_match = natom_match + 1
525 113 : exit
526 : end if
527 : end do
528 : end do
529 51 : if (natom_match.eq.MD%natom_unitcell) then
530 : iatom_ref = iatom
531 : ok=.true.
532 : exit
533 7 : else if (natom_match.gt.MD%natom_unitcell) then
534 0 : write(msg,'(5a)') 'Too many atoms match the unit cell.',ch10,&
535 0 : 'Perhaps the value of tolmotif is too large,',ch10,&
536 0 : 'or the value of tolinbox is too small.'
537 0 : ABI_ERROR(msg)
538 : endif
539 : end do
540 : if (.not.ok) then
541 0 : call tdep_write_xred_average(Invar,MPIdata,Lattice,xred_ideal,xred_center)
542 0 : write(msg,'(3a)') 'The basis of atoms written in input.in file does not appear in the MD trajectory.',ch10,&
543 0 : 'Perhaps, you can adjust the tolerance (tolmotif).'
544 0 : ABI_ERROR(msg)
545 : end if
546 44 : ABI_FREE(dist_unitcell)
547 :
548 44 : write(Invar%stdout,*)' Compare ideal and average positions using PBC...'
549 : ! Modification of xred and Rlatt tabs
550 : ! for averaged quantities: xred_center, xred_average, xred
551 : ! 1/ The "iatom_ref" atom is put in (0.0;0.0;0.0)
552 176 : tmp(:) = xred_center(:,iatom_ref)
553 6926 : do jatom=1,MD%natom
554 27572 : xred_center(:,jatom) = xred_center(:,jatom) - tmp(:)
555 : end do
556 : ! 2/ All the atoms are put in the range [-0.5;0.5[ (use of PBC)
557 6926 : do jatom=1,MD%natom
558 27528 : tmp(:)=xred_center(:,jatom)
559 6882 : call tdep_make_inbox(tmp,1,Invar%tolinbox,xred_center(:,jatom))
560 6882 : call tdep_make_inbox(tmp,1,Invar%tolinbox,xred_average(:,jatom))
561 146078 : do istep=1,MD%my_nstep
562 146034 : call tdep_make_inbox(tmp,1,Invar%tolinbox,MD%xred(:,jatom,istep))
563 : end do
564 : end do
565 : ! Modification of xred and Rlatt tabs
566 : ! for ideal quantities: Rlatt_red et xred_ideal
567 : ! 1/ The atom 1 is put in (0.0;0.0;0.0)
568 176 : tmp1(:)=xred_ideal(:,1)
569 176 : tmp2(:)=Rlatt_red(:,1,1)
570 6926 : do jatom=1,MD%natom
571 27528 : xred_ideal(:,jatom)= xred_ideal(:,jatom) -tmp1(:)
572 27572 : Rlatt_red (:,1,jatom)=Rlatt_red (:,1,jatom)-tmp2(:)
573 : end do
574 : ! 2/ All the atoms are put in the range [-0.5;0.5[ (use of PBC)
575 6926 : do jatom=1,MD%natom
576 27528 : tmp(:)=xred_ideal(:,jatom)
577 6882 : call tdep_make_inbox(tmp,1,tol8,xred_ideal(:,jatom))
578 6926 : call tdep_make_inbox(tmp,1,tol8,Rlatt_red(:,1,jatom))
579 : !FB call tdep_make_inbox(Rlatt_red(:,1,jatom),1,tol8)
580 : end do
581 :
582 : ! When the multiplicity equals 1 along one direction, there is some trouble
583 : ! To clean!!!!!!!
584 176 : do ii=1,3
585 : if ((Invar%multiplicity(ii,ii).eq.1).and.(Invar%multiplicity(ii,mod(ii ,3)+1).eq.0)&
586 176 : & .and.(Invar%multiplicity(ii,mod(ii+1,3)+1).eq.0)) then
587 65 : Rlatt_red(ii,1,:)=0.d0
588 1 : write(Invar%stdout,*) 'WARNING: multiplicity=1 for ii=',ii
589 : end if
590 : end do
591 :
592 : ! Define Rlatt for all the atoms in the basis (Rlatt_red varies as a function of iatcell)
593 44 : if (MD%natom_unitcell.gt.1) then
594 93 : do iatcell=2,MD%natom_unitcell
595 32109 : Rlatt_red(:,iatcell,:)=Rlatt_red(:,1,:)
596 : end do
597 : end if
598 6926 : do iatom=1,MD%natom
599 21812 : do iatcell=1,MD%natom_unitcell
600 59544 : tmp(:)=xred_ideal(:,iatom)-xred_ideal(:,iatcell)
601 21768 : call tdep_make_inbox(tmp,1,tol8,Rlatt_red(:,iatcell,iatom))
602 : end do
603 : end do
604 44 : if (Invar%debug) then
605 12 : do iatcell=1,MD%natom_unitcell
606 10 : write(Invar%stdout,*) 'For iatcell=',iatcell
607 1036 : do jatom=1,MD%natom
608 1034 : write(Invar%stdout,'(a,i4,a,3(f16.10,1x))') 'For jatom=',jatom,', Rlatt=',Rlatt_red(1:3,iatcell,jatom)
609 : end do
610 : end do
611 : end if
612 :
613 : ! Matching between Ideal and Average positions: xred_ideal and xred_center
614 : ! Then, write them in the xred_average.xyz file.
615 44 : write(Invar%stdout,*)' Write the xred_average.xyz file with ideal and average positions...'
616 7014 : ABI_CALLOC(FromIdeal2Average,(MD%natom))
617 6926 : do iatom=1,MD%natom
618 6882 : ok =.false.
619 739595 : do jatom=1,MD%natom
620 739595 : if (MD%typat(iatom).ne.MD%typat_unitcell(mod(jatom-1,MD%natom_unitcell)+1)) cycle
621 624459 : must_shift=.false.
622 624459 : ndir_match=0
623 2497836 : do ii=1,3
624 2497836 : if (abs(xred_center(ii,iatom)-xred_ideal(ii,jatom)).le.Invar%tolmatch) then
625 246347 : ndir_match=ndir_match+1
626 : else if ((abs(xred_center(ii,iatom)-xred_ideal(ii,jatom)-1.d0).le.Invar%tolmatch) &
627 1627030 : & .or.(abs(xred_center(ii,iatom)-xred_ideal(ii,jatom)+1.d0).le.Invar%tolmatch)) then
628 1114 : ndir_match=ndir_match+1
629 1114 : must_shift=.true.
630 : endif
631 : end do
632 624459 : if (ndir_match==3.and..not.must_shift) then
633 6882 : FromIdeal2Average(jatom)=iatom
634 : ok=.true.
635 : exit
636 617577 : else if (ndir_match==3.and.must_shift) then
637 0 : do ii=1,3
638 0 : if (abs(xred_center(ii,iatom)-xred_ideal(ii,jatom)-1.d0).le.Invar%tolmatch) then
639 0 : xred_center(ii,iatom)=xred_center(ii,iatom)-1d0
640 0 : xred_average(ii,iatom)=xred_average(ii,iatom)-1d0
641 0 : do istep=1,MD%my_nstep
642 0 : MD%xred(ii,iatom,istep)=MD%xred(ii,iatom,istep)-1d0
643 : end do
644 0 : FromIdeal2Average(jatom)=iatom
645 0 : else if (abs(xred_center(ii,iatom)-xred_ideal(ii,jatom)+1.d0).le.Invar%tolmatch) then
646 0 : xred_center(ii,iatom)=xred_center(ii,iatom)+1d0
647 0 : xred_average(ii,iatom)=xred_average(ii,iatom)+1d0
648 0 : do istep=1,MD%my_nstep
649 0 : MD%xred(ii,iatom,istep)=MD%xred(ii,iatom,istep)+1d0
650 : end do
651 0 : FromIdeal2Average(jatom)=iatom
652 : end if
653 : end do
654 : ok=.true.
655 : exit
656 : end if
657 : end do
658 44 : if (.not.ok) then
659 0 : write(Invar%stdlog,*) 'Problem to find the average position for iatom=',iatom
660 0 : write(Invar%stdlog,*) ' Reasons:'
661 0 : write(Invar%stdlog,*) ' 1/ One atom jump to another equilibrium position'
662 0 : write(Invar%stdlog,*) ' 2/ The system is no more solid'
663 0 : write(Invar%stdlog,*) ' 3/ Perhaps, you can adjust the tolerance (tolmatch)'
664 0 : write(Invar%stdlog,*) ' xred_center=',(xred_center(ii,iatom),ii=1,3)
665 0 : do eatom=1,MD%natom
666 0 : write(Invar%stdlog,'(a,1x,3(f10.6,1x))') 'I',xred_ideal (:,eatom)
667 0 : write(Invar%stdlog,'(a,1x,3(f10.6,1x))') 'C',xred_center(:,eatom)
668 : end do
669 0 : ABI_ERROR('Problem to find the average position')
670 : end if
671 : end do
672 :
673 : ! WARNING: VERY IMPORTANT: The positions are displayed/sorted
674 : ! (and used in the following) according to ideal positions xred_ideal.
675 44 : call tdep_write_xred_average(Invar,MPIdata,Lattice,xred_ideal,xred_center,FromIdeal2Average)
676 44 : ABI_FREE(xred_center)
677 :
678 : !====================================================================================
679 : !====================== END OF REDUCED COORDINATES ==================================
680 : !====================================================================================
681 : ! a/ Get cartesian coordinates from reduced ones
682 : ! b/ Compute ucart and fcart tabs
683 : ! c/ The atoms are sorted according the IDEAL arrangement
684 : ! The correspondance function is contained in: FromIdeal2Average
685 : ! WARNING : Consequently the arrangement of the xcart* tabs is not modified.
686 44 : write(Invar%stdout,*)' Compute cartesian coordinates and forces...'
687 557736 : ABI_MALLOC(xcart ,(3,MD%natom,MD%my_nstep)); xcart(:,:,:)=0.d0
688 27660 : ABI_MALLOC(xcart_ideal ,(3,MD%natom)) ; xcart_ideal(:,:)=0.d0
689 27616 : ABI_MALLOC(xcart_average,(3,MD%natom)) ; xcart_average(:,:)=0.d0
690 557692 : ABI_MALLOC(ucart_tmp ,(3,MD%natom,MD%my_nstep)); ucart_tmp(:,:,:)=0.d0
691 6926 : do iatom=1,MD%natom
692 6882 : call DGEMV('T',3,3,1.d0,Lattice%rprimd_md(:,:),3,xred_ideal (:,iatom),1,0.d0,xcart_ideal (:,iatom),1)
693 6882 : call DGEMV('T',3,3,1.d0,Lattice%rprimd_md(:,:),3,xred_average(:,iatom),1,0.d0,xcart_average(:,iatom),1)
694 21812 : do iatcell=1,MD%natom_unitcell
695 21768 : call DGEMV('T',3,3,1.d0,Lattice%rprimd_md(:,:),3,Rlatt_red(:,iatcell,iatom),1,0.d0,MD%Rlatt_cart(:,iatcell,iatom),1)
696 : end do
697 : end do
698 996 : do istep=1,MD%my_nstep
699 140148 : do iatom=1,MD%natom
700 139152 : jatom = FromIdeal2Average(iatom)
701 : call DGEMV('T',3,3,1.d0,Lattice%rprimd_md(:,:),3,MD%xred(:,jatom,istep),&
702 139152 : & 1,0.d0,xcart(:,jatom,istep),1)
703 140104 : if (Invar%use_ideal_positions.eq.0) then
704 115648 : ucart_tmp(:,iatom,istep) = xcart(:,jatom,istep) - xcart_average(:,jatom)
705 : else
706 440960 : ucart_tmp(:,iatom,istep) = xcart(:,jatom,istep) - xcart_ideal(:,iatom)
707 : end if
708 : end do
709 : end do
710 44 : ABI_FREE(xred_average)
711 44 : ABI_FREE(xcart)
712 44 : ABI_FREE(xcart_ideal)
713 44 : ABI_FREE(xcart_average)
714 :
715 : ! Rearrangement of the fcart tabs in column --> MD%Forces
716 557736 : ABI_CALLOC(fcart_tmp,(3,MD%natom,MD%my_nstep))
717 996 : do istep=1,MD%my_nstep
718 140148 : do iatom=1,MD%natom
719 557560 : fcart_tmp(:,iatom,istep) = MD%fcart(:,FromIdeal2Average(iatom),istep)
720 : end do
721 : end do
722 996 : do istep=1,MD%my_nstep
723 140148 : do jatom=1,MD%natom
724 557560 : do ii=1,3
725 417456 : jj = ii + 3*(jatom-1) + 3*MD%natom*(istep-1)
726 417456 : MD%Forces(jj) = fcart_tmp(ii,jatom,istep)
727 556608 : MD%ucart(ii,jatom,istep) = ucart_tmp(ii,jatom,istep)
728 : enddo
729 : enddo
730 : enddo
731 44 : ABI_FREE(FromIdeal2Average)
732 44 : ABI_FREE(ucart_tmp)
733 44 : ABI_FREE(fcart_tmp)
734 :
735 : ! Define Rlatt_scaled, fulfilling the definition of mkphdos (ABINIT routine)
736 176 : do ii=1,3
737 572 : rprimd_md_tmp(ii,:) = Lattice%rprimd_md(ii,:) / Lattice%acell_unitcell(ii)
738 : end do
739 6926 : do iatom=1,MD%natom
740 21812 : do iatcell=1,MD%natom_unitcell
741 21768 : call DGEMV('T',3,3,1.d0,rprimd_md_tmp,3,Rlatt_red(:,iatcell,iatom),1,0.d0,MD%Rlatt_scaled(:,iatcell,iatom),1)
742 : end do
743 : end do
744 :
745 : ! Find the symetry operation between 2 atoms
746 44 : call tdep_SearchS_1at(Invar,MPIdata,Sym,xred_ideal)
747 27572 : MD%xred_ideal(:,:)=xred_ideal(:,:)
748 44 : ABI_FREE(xred_ideal)
749 44 : ABI_FREE(Rlatt_red)
750 :
751 44 : end subroutine tdep_MatchIdeal2Average
752 :
753 : !====================================================================================================
754 :
755 44 : subroutine tdep_write_xred_average(Invar,MPIdata,Lattice,&
756 44 : xred_ideal,xred_center,&
757 44 : FromIdeal2Average)
758 : type(atdep_dataset_type), intent(in) :: Invar
759 : type(MPI_enreg_type),intent(in) :: MPIdata
760 : type(Lattice_type),intent(in) :: Lattice
761 : double precision,intent(in) :: xred_ideal(3,Invar%natom)
762 : double precision,intent(in) :: xred_center(3,Invar%natom)
763 : integer,intent(in),optional :: FromIdeal2Average(Invar%natom)
764 :
765 : integer :: unt
766 : !integer :: natom,natom_unitcell
767 : integer :: iatom,jatom,ii,jj
768 : !logical :: with_xcart
769 44 : integer,allocatable :: ideal2average(:)
770 : double precision :: rprimd(3,3)
771 : double precision :: xred_C(3),xred_I(3),xcart_C(3),xcart_I(3)
772 :
773 44 : if (MPIdata%iam_master) then
774 :
775 572 : rprimd(:,:) = Lattice%rprimd_md(:,:)
776 :
777 132 : ABI_MALLOC(ideal2average,(Invar%natom))
778 6926 : ideal2average(:)=0
779 44 : if (present(FromIdeal2Average)) then
780 6926 : ideal2average(:) = FromIdeal2Average(:)
781 : else
782 0 : do iatom=1,Invar%natom
783 0 : ideal2average(iatom) = iatom
784 : end do
785 : end if
786 :
787 44 : unt=31
788 44 : open(unit=unt,file=trim(Invar%output_prefix)//'_xred_average.xyz')
789 44 : write(unt,'(a,i4)') '# natom = ',Invar%natom
790 44 : write(unt,'(a,i4)') '# natom_unitcell = ',Invar%natom_unitcell
791 572 : write(unt,'(a,9(f4.1,1x))') '# multiplicity = ',((Lattice%multiplicity(ii,jj),jj=1,3),ii=1,3 )
792 44 : write(unt,'(a)') '#'
793 :
794 44 : write(unt,'(a1,1x,a8,2x,a6,2x,2(a5,30x))') '#', 'position', 'iatom', 'xred ', 'xcart'
795 44 : write(unt,'(a)')''
796 :
797 44 : xred_I = zero
798 44 : xred_C = zero
799 6926 : do iatom=1,Invar%natom
800 6882 : jatom = ideal2average(iatom)
801 27528 : xred_I = xred_ideal (:,iatom)
802 27528 : xred_C = xred_center(:,jatom)
803 :
804 6882 : xcart_I(:)=zero
805 6882 : xcart_C(:)=zero
806 6882 : call DGEMV('T',3,3,1.d0,rprimd(:,:),3,xred_I,1,0.d0,xcart_I,1)
807 6882 : call DGEMV('T',3,3,1.d0,rprimd(:,:),3,xred_C,1,0.d0,xcart_C,1)
808 :
809 6882 : write(unt,'(2x,a6,4x,i6,2x,3(f10.6,1x),2x,3(f10.6,1x))')'Ideal ',iatom,xred_I,xcart_I
810 6882 : write(unt,'(2x,a6,4x,i6,2x,3(f10.6,1x),2x,3(f10.6,1x))')'Center',jatom,xred_C,xcart_C
811 6926 : write(unt,'(a)')''
812 :
813 : end do
814 :
815 44 : close(unt)
816 44 : ABI_FREE(ideal2average)
817 : end if
818 :
819 44 : end subroutine tdep_write_xred_average
820 :
821 : !====================================================================================================
822 :
823 139152 : end module m_tdep_sampling
824 : !!***
|
