Branch data Line data Source code
1 : : !!****m* ABINIT/m_supercell
2 : : !! NAME
3 : : !! m_supercell
4 : : !!
5 : : !! FUNCTION
6 : : !! Module for using a supercell, in particular for phonon displacement freezing.
7 : : !! Container type is defined, and destruction, print subroutines as well as the central supercell_init
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2010-2026 ABINIT group (MJV, DJA)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public Licence, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
15 : : !!
16 : : !! SOURCE
17 : :
18 : : #if defined HAVE_CONFIG_H
19 : : #include "config.h"
20 : : #endif
21 : :
22 : : #include "abi_common.h"
23 : :
24 : : module m_supercell
25 : :
26 : : use defs_basis
27 : : use m_errors
28 : : use m_abicore
29 : :
30 : : use m_matrix, only : matr3inv
31 : : use m_copy, only : alloc_copy
32 : : use m_io_tools, only : open_file
33 : : use m_fstrings, only : int2char4, write_num, itoa, sjoin
34 : : use m_numeric_tools, only : isdiagmat
35 : :
36 : : implicit none
37 : :
38 : : private
39 : :
40 : : !!***
41 : :
42 : : !!****t* m_supercell/supercell_type
43 : : !! NAME
44 : : !! supercell_type
45 : : !!
46 : : !! FUNCTION
47 : : !! structure for a supercell constructed from a basic rprimd and xcart, with indexing to original atoms
48 : : !! The supercell may not be oriented the same way as the original cell, if you can reduce it by symmetry
49 : : !!
50 : : !! SOURCE
51 : :
52 : : type, public :: supercell_type
53 : : integer :: natom_primcell
54 : : ! number of atoms in primitive cell
55 : : integer :: natom
56 : : ! number of atoms in supercell
57 : : integer :: ntypat
58 : : ! number of atom types
59 : : integer :: ncells
60 : : ! number of unit cells in supercell
61 : : integer :: rlatt(3,3)
62 : : ! matrix for multiplicity of supercell (for the time being must be diagonal)
63 : : real(dp) :: rprimd(3,3)
64 : : ! new lattice vectors for supercell
65 : : real(dp) :: qphon(3)
66 : : ! phonon q vector used to generate scell, if any
67 : : character(len=3) :: xyz_order
68 : : ! Order used to build the supercell.
69 : : real(dp), allocatable :: xcart(:,:)
70 : : ! (3, natom) Cartesian positions of atoms
71 : : real(dp), allocatable :: xcart_ref(:,:)
72 : : ! (3, natom) equilibrium Cartesian positions of atoms
73 : : integer, allocatable :: atom_indexing(:)
74 : : ! (natom) indexes original atom: 1..natom_primcell
75 : : integer, allocatable :: uc_indexing(:,:)
76 : : ! (3, natom) indexes unit cell atom is.
77 : : integer, allocatable :: typat(:)
78 : : ! (natom) type of each atom in the supercell.
79 : : real(dp), allocatable :: znucl(:)
80 : : ! (ntypat) nuclear charges of species
81 : : integer, allocatable :: rvecs(:,:)
82 : : ! supercell vectors
83 : :
84 : : contains
85 : : procedure :: init_for_qpt => supercell_init_for_qpt
86 : : procedure :: init => supercell_init
87 : : procedure :: freeze_displ => supercell_freeze_displ
88 : : procedure :: copy => supercell_copy
89 : : procedure :: free => supercell_free
90 : : procedure :: print_for_qpt => supercell_print_for_qpt
91 : : procedure :: print_abinit => supercell_print_abinit
92 : : procedure :: write_xsf => supercell_write_xsf
93 : : end type supercell_type
94 : :
95 : : public :: distance_supercell
96 : : public :: findBound_supercell
97 : : public :: getPBCIndexes_supercell
98 : : public :: mksupercell ! computes atomic positons, magnetic ordering of supercell
99 : : !!***
100 : :
101 : : CONTAINS !===========================================================================================
102 : :
103 : : !!****f* m_supercell/supercell_init_for_qpt
104 : : !!
105 : : !! NAME
106 : : !! supercell_init_for_qpt
107 : : !!
108 : : !! FUNCTION
109 : : !! Initialize scell structure, from unit cell vectors, and atoms, based on qpoint chosen
110 : : !!
111 : : !! INPUTS
112 : : !! natom_primcell = number of atoms in primitive cell
113 : : !! qphon(3) = phonon wavevector
114 : : !! find smallest supercell which will accomodate phonon qphon = (1/2,1/2,1/2)
115 : : !! rprimd_primcell(3,3) = real space lattice vectors (bohr)
116 : : !! typat_primcell = types of atoms
117 : : !! xcart_primcell(3,natom) = cartesian positions of atoms in primitive cell
118 : : !! znucl = nuclear charges for all species
119 : : !! ordering = if true, typat will be 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 ....
120 : : !! if false, typat will be 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 ....
121 : : !!
122 : : !! OUTPUT
123 : : !! scell = supercell structure to be initialized
124 : : !!
125 : : !! SOURCE
126 : :
127 : 6 : subroutine supercell_init_for_qpt(scell, natom_primcell, qphon, rprimd_primcell, &
128 [ - + + - ]: 6 : typat_primcell, xcart_primcell, znucl, ordering)
129 : :
130 : : !Arguments ------------------------------------
131 : : !scalars
132 : : class(supercell_type), intent(out) :: scell
133 : : integer, intent(in) :: natom_primcell
134 : : logical,optional,intent(in) :: ordering
135 : : !arrays
136 : : integer , intent(in) :: typat_primcell(natom_primcell)
137 : : real(dp), intent(in) :: qphon(3)
138 : : real(dp), intent(in) :: znucl(:)
139 : : real(dp), intent(in) :: rprimd_primcell(3,3)
140 : : real(dp), intent(in) :: xcart_primcell(3,natom_primcell)
141 : :
142 : : !Local variables-------------------------------
143 : : !scalar
144 : : integer :: ii, maxsc, iscmult
145 : : real(dp) :: qbymult
146 : : !arrays
147 : : integer :: rlatt(3,3) ! number of primitive cells in each direction for the supercell
148 : : character(len=500) :: msg
149 : : ! *************************************************************************
150 : :
151 : : ! maximum number of unit cells in a given direction
152 : 6 : maxsc = 10
153 : :
154 : : ! find smallest supercell which will accomodate phonon.
155 : : ! FIXME: for the moment, just get smallest multiple along each direction, with an upper bound
156 : 6 : rlatt = 0
157 : 6 : rlatt(1,1) = -1
158 : 6 : rlatt(2,2) = -1
159 : 6 : rlatt(3,3) = -1
160 [ + + ]: 24 : do ii=1,3
161 [ + - ]: 21 : do iscmult=1,maxsc
162 : 21 : qbymult = qphon(ii)*iscmult
163 [ + + ]: 21 : if (abs(qbymult - int(qbymult)) < tol10) then
164 : 18 : rlatt(ii,ii) = iscmult
165 : 18 : exit
166 : : end if
167 : : end do
168 [ - + ]: 24 : if (rlatt(ii,ii) == -1) then
169 : 0 : write(msg,'(a,I0,a,I0,2a,3E20.10)')' No supercell found with less than ', &
170 : 0 : maxsc,' unit cells in direction ', ii, ch10, ' qphon = ', qphon
171 : 0 : ABI_ERROR(msg)
172 : : end if
173 : : end do
174 : :
175 [ - + ]: 6 : if (present(ordering)) then
176 : 0 : call scell%init(natom_primcell, rlatt, rprimd_primcell, typat_primcell, xcart_primcell, znucl, ordering)
177 : : else
178 : 6 : call scell%init(natom_primcell, rlatt, rprimd_primcell, typat_primcell, xcart_primcell, znucl)
179 : : end if
180 : :
181 [ + + ]: 24 : scell%qphon = qphon
182 : :
183 : 6 : end subroutine supercell_init_for_qpt
184 : : !!***
185 : :
186 : : !!****f* m_supercell/supercell_init
187 : : !! NAME
188 : : !! supercell_init
189 : : !!
190 : : !! FUNCTION
191 : : !! Initialize scell structure, from unit cell vectors, and atoms, based on rlatt multiplicity matrix
192 : : !!
193 : : !! INPUTS
194 : : !! natom_primcell = number of atoms in primitive cell
195 : : !! rlatt(3,3) = multiplicity of primtive unit cells in supercell
196 : : !! rprimd_primcell(3,3) = real space lattice vectors (bohr)
197 : : !! typat_primcell(natom) = types of all atoms in primitive cell
198 : : !! xcart_primcell(3,natom) = cartesian positions of atoms in primitive cell
199 : : !! znucl = nuclear charges for all species
200 : : !! [ordering] = if true, typat will be 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 ....
201 : : !! if false, typat will be 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 ....
202 : : !! [xyz_order]= Order used to build the supercell.
203 : : !! "zyx" if one should move along z first.
204 : : !! "xyz" if one should move along x first. This value should be used for plotting purposes
205 : : !!
206 : : !! OUTPUT
207 : : !! scell = supercell structure to be initialized
208 : : !!
209 : : !! SOURCE
210 : :
211 [ - + + - ]: 212 : subroutine supercell_init(scell, natom_primcell, rlatt, rprimd_primcell, typat_primcell, xcart_primcell, znucl, &
212 : : ordering, xyz_order) ! optional
213 : :
214 : : !Arguments ------------------------------------
215 : : !scalars
216 : : class(supercell_type), intent(out) :: scell
217 : : integer, intent(in) :: natom_primcell
218 : : logical,optional,intent(in) :: ordering
219 : : character(len=*),optional,intent(in) :: xyz_order
220 : : !arrays
221 : : integer , intent(in) :: rlatt(3,3)
222 : : integer , intent(in) :: typat_primcell(natom_primcell)
223 : : real(dp), intent(in) :: znucl(:)
224 : : real(dp), intent(in) :: rprimd_primcell(3,3)
225 : : real(dp), intent(in) :: xcart_primcell(3,natom_primcell)
226 : :
227 : : !Local variables-------------------------------
228 : : !scalars
229 : : integer :: iatom_supercell, i1,i2,i3, iatom, icell
230 : : ! *************************************************************************
231 : :
232 [ - + ]: 212 : if (.not. isdiagmat(rlatt)) then
233 : 0 : ABI_ERROR('rlatt is not diagonal.')
234 : : end if
235 : :
236 [ + + ]: 212 : scell%xyz_order = "zyx"; if (present(xyz_order)) scell%xyz_order = xyz_order(1:3)
237 : :
238 : 212 : scell%natom_primcell = natom_primcell
239 [ + + + + ]: 2756 : scell%rlatt = rlatt
240 : 212 : scell%ncells = rlatt(1,1)*rlatt(2,2)*rlatt(3,3)
241 [ + + ]: 848 : scell%rprimd(:,1) = rprimd_primcell(:,1) * rlatt(1,1)
242 [ + + ]: 848 : scell%rprimd(:,2) = rprimd_primcell(:,2) * rlatt(2,2)
243 [ + + ]: 848 : scell%rprimd(:,3) = rprimd_primcell(:,3) * rlatt(3,3)
244 : :
245 : : !call metric(scell%gmet, scell%gprimd, -1, scell%rmet, scell%rprimd, scell%ucvol)
246 : :
247 : 212 : scell%ntypat = size(znucl)
248 [ - + + - : 636 : ABI_MALLOC(scell%znucl,(scell%ntypat))
+ - + - -
+ - + -
+ ]
249 [ + + ]: 751 : scell%znucl(:) = znucl(:)
250 : :
251 : : ! number of atoms in full supercell
252 : 212 : scell%natom= natom_primcell*scell%ncells
253 [ - + + - : 636 : ABI_MALLOC(scell%xcart,(3,scell%natom))
+ - + - -
+ - + -
+ ]
254 [ - + + - : 636 : ABI_MALLOC(scell%xcart_ref,(3,scell%natom))
+ - + - -
+ - + -
+ ]
255 [ - + + - : 636 : ABI_MALLOC(scell%typat,(scell%natom))
+ - + - -
+ - + -
+ ]
256 [ - + + - : 636 : ABI_MALLOC(scell%atom_indexing,(scell%natom))
+ - + - -
+ - + -
+ ]
257 [ - + + - : 636 : ABI_MALLOC(scell%uc_indexing,(3,scell%natom))
+ - + - -
+ - + -
+ ]
258 [ - + + - : 636 : ABI_MALLOC(scell%rvecs, (3, scell%ncells))
+ - + - -
+ - + -
+ ]
259 : :
260 : 212 : iatom_supercell = 0; icell =0
261 : :
262 : 199 : select case (scell%xyz_order)
263 : : case ("zyx")
264 : : ! legacy mode.
265 [ + + ]: 604 : do i1 = 1, rlatt(1,1)
266 [ + + ]: 1607 : do i2 = 1, rlatt(2,2)
267 [ + + ]: 3523 : do i3 = 1, rlatt(3,3)
268 : 3118 : call build_()
269 : : end do
270 : : end do
271 : : end do
272 : :
273 : : case ("xyz")
274 [ + + ]: 90 : do i3 = 1, rlatt(3,3)
275 [ + + ]: 749 : do i2 = 1, rlatt(2,2)
276 [ + + ]: 7965 : do i1 = 1, rlatt(1,1)
277 : 7888 : call build_()
278 : : end do
279 : : end do
280 : : end do
281 : :
282 : : case default
283 [ + + - ]: 212 : ABI_ERROR(sjoin("Invalid xyz_order", scell%xyz_order))
284 : : end select
285 : :
286 [ - + ]: 212 : ABI_CHECK_IEQ(iatom_supercell, scell%natom, "iatom_supercell /= scell%natom")
287 : :
288 [ + - + - : 107088 : scell%xcart = scell%xcart_ref
- + - - -
- - - + +
+ + ]
289 [ + + ]: 848 : scell%qphon = zero
290 : :
291 [ - + ]: 212 : if (present(ordering)) then
292 [ # # ]: 0 : if (ordering) call order_supercell_typat(scell)
293 : : end if
294 : :
295 : : contains
296 : 9344 : subroutine build_()
297 [ + + ]: 37376 : icell = icell+1; scell%rvecs(:,icell) = [i1-1, i2-1, i3-1]
298 [ + + ]: 36010 : do iatom = 1, natom_primcell
299 : 26666 : iatom_supercell = iatom_supercell + 1
300 [ + + ]: 106664 : scell%uc_indexing(:,iatom_supercell) = [i1-1, i2-1, i3-1]
301 : : scell%xcart_ref(:,iatom_supercell) = xcart_primcell(:,iatom) &
302 [ + + + + : 586652 : + matmul(rprimd_primcell,scell%uc_indexing(:,iatom_supercell))
+ + + + +
- + - +
+ ]
303 : 26666 : scell%atom_indexing(iatom_supercell) = iatom
304 : 36010 : scell%typat(iatom_supercell) = typat_primcell(iatom)
305 : : end do
306 : 9344 : end subroutine build_
307 : :
308 : : end subroutine supercell_init
309 : : !!***
310 : :
311 : : !!****f* m_supercell/order_supercell_typat
312 : : !!
313 : : !! NAME
314 : : !! order_supercell_typat
315 : : !!
316 : : !! FUNCTION
317 : : !! Re-order atoms in place for types
318 : : !!
319 : : !! INPUTS
320 : : !! scell = supercell structure with reference atomic positions etc...
321 : : !!
322 : : !! OUTPUT
323 : : !! scell = supercell structure: typat, xcart and so on will be updated
324 : : !!
325 : : !! SOURCE
326 : :
327 : 0 : subroutine order_supercell_typat(scell)
328 : :
329 : : !Arguments ------------------------------------
330 : : !scalars
331 : : class(supercell_type), intent(inout) :: scell
332 : :
333 : : !Local variables-------------------------------
334 : : integer :: itypat, iatom_supercell, iatom
335 : 0 : type(supercell_type) :: scell_tmp
336 : : ! *************************************************************************
337 : :
338 : 0 : call scell%copy(scell_tmp)
339 : :
340 : 0 : iatom_supercell = 0
341 [ # # ]: 0 : do itypat = 1, scell%ntypat
342 [ # # ]: 0 : do iatom = 1, scell%natom
343 [ # # ]: 0 : if (scell_tmp%typat(iatom) /= itypat) cycle
344 : 0 : iatom_supercell = iatom_supercell + 1
345 [ # # ]: 0 : scell%xcart(:,iatom_supercell) = scell_tmp%xcart(:,iatom)
346 [ # # ]: 0 : scell%xcart_ref(:,iatom_supercell) = scell_tmp%xcart_ref(:,iatom)
347 : 0 : scell%atom_indexing(iatom_supercell) = scell_tmp%atom_indexing(iatom)
348 [ # # ]: 0 : scell%uc_indexing(:,iatom_supercell) = scell_tmp%uc_indexing(:,iatom)
349 : 0 : scell%typat(iatom_supercell) = scell_tmp%typat(iatom)
350 : : end do
351 : : end do
352 : :
353 : 0 : call scell_tmp%free()
354 : :
355 [ # # # # : 0 : end subroutine order_supercell_typat
# # # # #
# # # #
# ]
356 : : !!***
357 : :
358 : :
359 : : !!****f* m_supercell/freeze_displ_supercell
360 : : !!
361 : : !! NAME
362 : : !! freeze_displ_supercell
363 : : !!
364 : : !! FUNCTION
365 : : !! Freeze a specific displacement phonon field into the supercell scell
366 : : !!
367 : : !! INPUTS
368 : : !! displ = phonon displacement vectors for this mode
369 : : !! freeze_displ = desired amplitude for phonon displacement along displ.
370 : : !! for thermal displacement use sqrt[ (1/2 + bose_einstein(freq,T)) / freq ]
371 : : !! scell = supercell structure with reference atomic positions etc...
372 : : !!
373 : : !! OUTPUT
374 : : !! scell = supercell structure: xcart will be updated with phonon displacement
375 : : !!
376 : : !! SOURCE
377 : :
378 : 1446 : subroutine supercell_freeze_displ(scell, displ, freeze_displ)
379 : :
380 : : !Arguments ------------------------------------
381 : : !scalars
382 : : class(supercell_type), intent(inout) :: scell
383 : : real(dp), intent(in) :: freeze_displ
384 : : !arrays
385 : : real(dp), intent(in) :: displ(2,3*scell%natom_primcell)
386 : :
387 : : !Local variables-------------------------------
388 : : integer :: iatom, ipratom
389 : : complex(dp) :: expqdotr, j=cmplx(zero,one)
390 : : complex(dp) :: phase
391 : 2892 : complex(dp) :: zdispl(3,scell%natom_primcell)
392 : : ! *************************************************************************
393 : :
394 : : zdispl = (cmplx(reshape(displ(1,:), (/3,scell%natom_primcell/)),&
395 [ + + + + : 41358 : reshape(displ(2,:), (/3,scell%natom_primcell/))))
+ - + + +
+ ]
396 : :
397 : : ! fix gauge by imposing real displacement for first atom in first direction
398 : : ! multiply by normalized complex conjugate of first element
399 : : ! NB 6 March 2018: this may be imposing a positive (not just real) displacement for 1st atom along x!!!
400 : : ! That might be problematic below, though for the thermal displacement method freeze_displ swaps sign for each new mode
401 : 1446 : phase = cmplx(one,zero)
402 [ + + ]: 1446 : if (abs(zdispl(1,1)) > tol10) then
403 : 697 : phase = conjg(zdispl(1,1)) / abs(zdispl(1,1))
404 : : end if
405 : :
406 [ + + ]: 69234 : do iatom = 1, scell%natom
407 : : expqdotr = exp(j*two_pi*(scell%qphon(1)*scell%uc_indexing(1,iatom) &
408 : : +scell%qphon(2)*scell%uc_indexing(2,iatom) &
409 : 67788 : +scell%qphon(3)*scell%uc_indexing(3,iatom)))
410 : :
411 : : ! this is offset in zdispl vector due to primitive cell atom position
412 : 67788 : ipratom = scell%atom_indexing(iatom)
413 : :
414 : : !add real part of displacement times Bloch phase
415 : : scell%xcart(:,iatom) = scell%xcart(:,iatom) &
416 [ + + ]: 272598 : & + freeze_displ * real(expqdotr * zdispl(:,ipratom) * phase)
417 : :
418 : : ! scell%xcart(:,iatom) = scell%xcart(:,iatom) &
419 : : !& + freeze_displ * cos(qdotr) * displ(1,ipratom+1:ipratom+3) &
420 : : !& - freeze_displ * sin(qdotr) * displ(2,ipratom+1:ipratom+3)
421 : : end do
422 : :
423 : 1446 : end subroutine supercell_freeze_displ
424 : : !!***
425 : :
426 : : !****f* m_supercell/supercell_print_for_qpt
427 : : !!
428 : : !! NAME
429 : : !! supercell_print_for_qpt
430 : : !!
431 : : !! FUNCTION
432 : : !! output atomic positions, supercell vectors, etc... to a file. single qpoint and mode.
433 : : !!
434 : : !! INPUTS
435 : : !! freq = phonon frequency for mode jmode
436 : : !! jmode = mode which has been frozen into xcart contained in scell
437 : : !! outfile_radix = radix of file name to be written to
438 : : !!
439 : : !! OUTPUT
440 : : !! printing to file
441 : : !!
442 : : !! SOURCE
443 : :
444 : 36 : subroutine supercell_print_for_qpt(scell, freq, jmode, outfile_radix)
445 : :
446 : : !Arguments ------------------------------------
447 : : !scalars
448 : : class(supercell_type), intent(in) :: scell
449 : : real(dp), intent(in) :: freq
450 : : integer, intent(in) :: jmode
451 : : character(len=*), intent(in) :: outfile_radix
452 : :
453 : : !Local variables-------------------------------
454 : : !scalar
455 : : character(len=fnlen) :: filename
456 : : character(len=10) :: jmodestring
457 : : character(len=80) :: title1, title2
458 : : character(len=5) :: qphonstring1, qphonstring2, qphonstring3
459 : : ! *************************************************************************
460 : :
461 : : ! add suffix with mode and qpoint
462 : 36 : call int2char4(jmode, jmodestring)
463 [ - + ]: 36 : ABI_CHECK((jmodestring(1:1)/='#'),'Bug: string length too short!')
464 : :
465 : : ! qphonstring should be like 0.000_0.000_0.000
466 : 36 : call write_num(scell%qphon(1),qphonstring1,'(F5.3)')
467 : 36 : call write_num(scell%qphon(2),qphonstring2,'(F5.3)')
468 : 36 : call write_num(scell%qphon(3),qphonstring3,'(F5.3)')
469 : : filename = trim(outfile_radix) // "_qpt_" // qphonstring1 // "_" // qphonstring2 // &
470 [ + - ]: 36 : "_" // qphonstring3 // "_mode_" // trim(jmodestring)
471 : :
472 : 36 : write (title1, '(a,3E20.10)') '# phonon q point : ', scell%qphon
473 : 36 : write (title2, '(a,I7,a,E20.10)') '# phonon mode number : ', jmode, ' frequency ', freq
474 : :
475 : 36 : call scell%print_abinit(filename, title1, title2)
476 : :
477 : 36 : end subroutine supercell_print_for_qpt
478 : : !!***
479 : :
480 : : !****f* m_supercell/supercell_print_abinit
481 : : !! NAME
482 : : !! supercell_print_abinit
483 : : !!
484 : : !! FUNCTION
485 : : !! output atomic positions, supercell vectors, etc... to a file
486 : : !! in Abinit input format.
487 : : !!
488 : : !! INPUTS
489 : : !! filename = filename
490 : : !! title1 = first line of description of contents
491 : : !! title2 = second line of description of contents
492 : : !! scell = supercell structure with data to be output
493 : : !!
494 : : !! OUTPUT
495 : : !! printing to file
496 : : !!
497 : : !! SOURCE
498 : :
499 : 46 : subroutine supercell_print_abinit(scell, filename, title1, title2)
500 : :
501 : : !Arguments ------------------------------------
502 : : !scalars
503 : : class(supercell_type), intent(in) :: scell
504 : : character(len=fnlen), intent(in) :: filename
505 : : character(len=80), intent(in) :: title1
506 : : character(len=80), intent(in) :: title2
507 : :
508 : : !Local variables-------------------------------
509 : : !scalar
510 : : integer :: scunit, iatom
511 : : character(len=500) :: msg
512 : : real(dp) :: xred(3), gprimd(3,3)
513 : : ! *************************************************************************
514 : :
515 [ - + ]: 46 : if (open_file(filename, msg, newunit=scunit, status="unknown", action="write") /= 0) then
516 : 0 : ABI_ERROR(msg)
517 : : end if
518 : :
519 : : ! print header
520 : 46 : write (scunit, '(a)') '#'
521 : 46 : write (scunit, '(a)') '# anaddb file with frozen phonon mode in supercell'
522 : 46 : write (scunit, '(a)') '# !!! Do not forget to adjust nband !!! '
523 : 46 : write (scunit, '(a)') '#'
524 : 46 : write (scunit, '(a)') title1
525 : 46 : write (scunit, '(a)') title2
526 : 46 : write (scunit, '(a,3(3I7,2x))') '# supercell rlatt is ', scell%rlatt
527 : 46 : write (scunit, '(a,I7,a)') '# and has ', scell%ncells, ' primitive unit cells '
528 : 46 : write (scunit, '(a)') '#'
529 : 46 : write (scunit, '(a)') '# lattice vectors for supercell :'
530 : 46 : write (scunit, '(a,I7)') 'natom ', scell%natom
531 : 46 : write (scunit, *)
532 : 46 : write (scunit, '(a)') 'znucl '
533 [ + + ]: 114 : do iatom = 1, size(scell%znucl)
534 : 68 : write (scunit, '(I5)', ADVANCE="NO") int(scell%znucl(iatom))
535 [ - + ]: 114 : if (mod(iatom,6) == 0) write (scunit, *)
536 : : end do
537 : 46 : write (scunit, *)
538 : 46 : write (scunit, *)
539 : 46 : write (scunit, '(a,I7)') 'ntypat', scell%ntypat
540 : 46 : write (scunit, '(a)') 'typat '
541 [ + + ]: 634 : do iatom = 1, scell%natom
542 : 588 : write (scunit, '(I5)', ADVANCE="NO") scell%typat(iatom)
543 [ + + ]: 634 : if (mod(iatom,6) == 0) write (scunit, *)
544 : : end do
545 : 46 : write (scunit, *)
546 : 46 : write (scunit, '(a)') 'acell 1.0 1.0 1.0'
547 : 46 : write (scunit, '(a)') 'rprim'
548 : 46 : write (scunit, '(3E20.10)') scell%rprimd(:,1)
549 : 46 : write (scunit, '(3E20.10)') scell%rprimd(:,2)
550 : 46 : write (scunit, '(3E20.10)') scell%rprimd(:,3)
551 : 46 : write (scunit, *)
552 : 46 : write (scunit, '(a)') 'xcart'
553 [ + + ]: 634 : do iatom = 1, scell%natom
554 : 634 : write (scunit, '(3E20.10)') scell%xcart(:,iatom)
555 : : end do
556 : : ! for information, also print xred for atoms inside full supercell
557 : 46 : call matr3inv(scell%rprimd, gprimd)
558 : : ! TODO: check this transpose is correct in some asymetric case
559 [ + + + + : 1196 : gprimd = transpose(gprimd)
+ + + + ]
560 : 46 : write (scunit, '(a)') '# for information, add xred as well'
561 : 46 : write (scunit, '(a)') '# xred'
562 [ + + ]: 634 : do iatom = 1, scell%natom
563 [ + + + + : 7644 : xred = matmul (gprimd, scell%xcart(:,iatom))
+ - ]
564 : 634 : write (scunit, '(a, 3E20.10)') '# ', xred
565 : : end do
566 : :
567 : : ! close file
568 : 46 : close(scunit)
569 : :
570 : 46 : end subroutine supercell_print_abinit
571 : : !!***
572 : :
573 : : !****f* m_supercell/supercell_copy
574 : : !!
575 : : !! NAME
576 : : !! supercell_copy
577 : : !!
578 : : !! FUNCTION
579 : : !! copy supercell structure
580 : : !!
581 : : !! INPUTS
582 : : !! scell_in = supercell structure with data to copy
583 : : !!
584 : : !! OUTPUT
585 : : !! scell = supercell structure with data to be output
586 : : !!
587 : : !! SOURCE
588 : :
589 : 105 : subroutine supercell_copy(scell_in, scell_copy)
590 : :
591 : : !Arguments ------------------------------------
592 : : class(supercell_type), intent(in) :: scell_in
593 : : class(supercell_type), intent(inout) :: scell_copy
594 : : ! *************************************************************************
595 : :
596 : 105 : call scell_copy%free()
597 : :
598 : 105 : scell_copy%natom_primcell = scell_in%natom_primcell
599 : 105 : scell_copy%natom = scell_in%natom
600 : 105 : scell_copy%ntypat = scell_in%ntypat
601 : 105 : scell_copy%ncells = scell_in%ncells
602 [ + + + + ]: 1365 : scell_copy%rlatt = scell_in%rlatt
603 [ + + + + ]: 1365 : scell_copy%rprimd = scell_in%rprimd
604 [ + + ]: 420 : scell_copy%qphon = scell_in%qphon
605 [ - + ]: 105 : call alloc_copy(scell_in%xcart , scell_copy%xcart)
606 [ - + ]: 105 : call alloc_copy(scell_in%xcart_ref , scell_copy%xcart_ref)
607 [ - + ]: 105 : call alloc_copy(scell_in%atom_indexing, scell_copy%atom_indexing)
608 [ - + ]: 105 : call alloc_copy(scell_in%uc_indexing , scell_copy%uc_indexing)
609 [ - + ]: 105 : call alloc_copy(scell_in%typat , scell_copy%typat)
610 [ - + ]: 105 : call alloc_copy(scell_in%znucl , scell_copy%znucl)
611 : :
612 : 105 : end subroutine supercell_copy
613 : : !!***
614 : :
615 : : !!****f* m_effective_potential/getPBCIndexes_supercell
616 : : !! NAME
617 : : !!
618 : : !! FUNCTION
619 : : !! Get the index of the cell by using PBC
620 : : !!
621 : : !! INPUTS
622 : : !! index = index of the cell into the supercell
623 : : !! ncell = number of total cell
624 : : !!
625 : : !! OUTPUT
626 : : !! index = index of the cell into the supercell with PBC
627 : : !!
628 : : !! SOURCE
629 : :
630 : 4411437984 : subroutine getPBCIndexes_supercell(index,ncell)
631 : :
632 : : !Arguments ---------------------------------------------
633 : : integer, intent(inout) :: index(3)
634 : : integer, intent(in) :: ncell(3)
635 : :
636 : : !Local variables ---------------------------------------
637 : : integer :: ii
638 : : ! *********************************************************************
639 : :
640 [ + + ]: 17645751936 : do ii=1,3
641 [ + + ]: 13848999534 : do while (index(ii) > ncell(ii))
642 : 13848999534 : index(ii) = index(ii) - ncell(ii)
643 : : end do
644 [ + + ]: 18203804122 : do while (index(ii) <= 0)
645 : 558052186 : index(ii) = index(ii) + ncell(ii)
646 : : end do
647 : : end do
648 : :
649 : 4411437984 : end subroutine getPBCIndexes_supercell
650 : : !!***
651 : :
652 : : !****f* m_supercell/findBound_supercell
653 : : !! NAME
654 : : !! findBound_supercell
655 : : !!
656 : : !! FUNCTION
657 : : !! compute the bound of the supercell by considering the 0 0 0 (reference)
658 : : !! in the center of the supercell.
659 : : !! for example: (4 4 4) => min = -1 and max = 2
660 : : !!
661 : : !! INPUTS
662 : : !! ncell(3) = size of the supercell (for example 3 3 3)
663 : : !!
664 : : !! OUTPUT
665 : : !! min = minimun of the range
666 : : !! max = maximum of the range
667 : : !!
668 : : !! SOURCE
669 : :
670 : 204 : subroutine findBound_supercell(min, max, ncell)
671 : :
672 : : !Arguments ---------------------------------------------
673 : : integer, intent(inout) :: min,max
674 : : integer, intent(in) :: ncell
675 : :
676 : : ! *********************************************************************
677 [ - + ]: 204 : if(abs(max)>abs(min)) then
678 [ # # ]: 0 : max=(ncell)/2; min=-max; if(mod(ncell,2)==0) max = max -1
679 : : else
680 [ + + ]: 204 : min=-(ncell)/2; max=-min; if(mod(ncell,2)==0) min= min +1
681 : : end if
682 : :
683 : 204 : end subroutine findBound_supercell
684 : : !!***
685 : :
686 : : !!****f* m_supercell/distance_supercell
687 : : !! NAME
688 : : !!
689 : : !! FUNCTION
690 : : !! compute the distance_supercell betwen 2 atoms in different cell
691 : : !!
692 : : !! INPUTS
693 : : !! xcart1(3) = cartesian coordinates of the first atom
694 : : !! xcart1(3) = cartesian coordinates of the second atom
695 : : !! rprimd(3,3) = primitive lattice vectors
696 : : !! cell1(3) = index of the cell of the first atom (for example -1 0 2)
697 : : !! cell2(3) = index of the cell of the second atom (for example 0 0 2)
698 : : !!
699 : : !! OUTPUT
700 : : !! distance_supercell = distance_supercell between the 2 atoms
701 : : !!
702 : : !! SOURCE
703 : : !!
704 : :
705 : 0 : pure real(dp) function distance_supercell(xcart1,xcart2,rprimd,cell1,cell2) result(dist)
706 : :
707 : : !Arguments ------------------------------------
708 : : real(dp),intent(in):: rprimd(3,3)
709 : : real(dp),intent(in):: xcart1(3),xcart2(3)
710 : : integer,intent(in) :: cell1(3),cell2(3)
711 : :
712 : : !Local variables -------------------------------
713 : : real(dp) :: rpt1(3),rpt2(3)
714 : : integer :: mu
715 : : ! *************************************************************************
716 : :
717 [ # # ]: 0 : do mu=1,3
718 : 0 : rpt1(mu) = cell1(1)*rprimd(mu,1)+cell1(2)*rprimd(mu,2)+cell1(3)*rprimd(mu,3)
719 : 0 : rpt2(mu) = cell2(1)*rprimd(mu,1)+cell2(2)*rprimd(mu,2)+cell2(3)*rprimd(mu,3)
720 : : end do
721 : :
722 : : dist = ((xcart2(1)+rpt2(1)-xcart1(1)-rpt1(1))**2+&
723 : : (xcart2(2)+rpt2(2)-xcart1(2)-rpt1(2))**2+&
724 : 0 : (xcart2(3)+rpt2(3)-xcart1(3)-rpt1(3))**2)**0.5
725 : :
726 : 0 : end function distance_supercell
727 : : !!***
728 : :
729 : : !****f* m_supercell/supercell_free
730 : : !!
731 : : !! NAME
732 : : !! supercell_free
733 : : !!
734 : : !! FUNCTION
735 : : !! deallocate all dynamic memory for this supercell structure
736 : : !!
737 : : !! SOURCE
738 : :
739 : 889 : subroutine supercell_free(scell)
740 : :
741 : : !Arguments ------------------------------------
742 : : class(supercell_type), intent(inout) :: scell
743 : : ! *************************************************************************
744 : :
745 [ + + ]: 889 : ABI_SFREE(scell%xcart)
746 [ + + ]: 889 : ABI_SFREE(scell%xcart_ref)
747 [ + + ]: 889 : ABI_SFREE(scell%typat)
748 [ + + ]: 889 : ABI_SFREE(scell%atom_indexing)
749 [ + + ]: 889 : ABI_SFREE(scell%uc_indexing)
750 [ + + ]: 889 : ABI_SFREE(scell%znucl)
751 [ + + ]: 889 : ABI_SFREE(scell%rvecs)
752 : :
753 : 889 : end subroutine supercell_free
754 : : !!***
755 : :
756 : : !!****f* m_supercell/mksupercell
757 : : !! NAME
758 : : !! mksupercell
759 : : !!
760 : : !! FUNCTION
761 : : !! computes atomic positons, magnetic ordering of supercell
762 : : !!
763 : : !! INPUTS
764 : : !! magv_org (optional) magnetic ordering of atoms in primitive cell,
765 : : !! ordering of atoms given als 1 and -1, if not given fm is assumed
766 : : !! xred_org relative position of atoms in primitive cell
767 : : !! rprimd_org unit cell dimensions of primitive cell
768 : : !! natom=number of atoms in unit cell
769 : : !! option= 1 output ion-ion distances / 2 output ordering of ion-ion distances / 3 output variables in varlist
770 : : !! according to ion-ion distances * magnetic ordering
771 : : !!
772 : : !! OUTPUT
773 : : !! magv_sc magnetic ordering of atoms in supercell
774 : : !! xred_sc relative position of atoms in supercell
775 : : !! rprimd_sc unit cell dimensions of supercell
776 : : !!
777 : : !! SOURCE
778 : :
779 : 45 : subroutine mksupercell(xred_org,magv_org,rprimd_org,nat_org,nat_sc,xred_sc,magv_sc,rprimd_sc,ext,prtvol)
780 : :
781 : : !Arguments ------------------------------------
782 : : !scalars
783 : : integer,intent(in) :: nat_org,nat_sc
784 : : integer,intent(in),optional :: prtvol
785 : : !arrays
786 : : real(dp),intent(in) :: rprimd_org(3,3)
787 : : integer,intent(in) :: ext(3)
788 : : real(dp),intent(in) :: xred_org(3,nat_org)
789 : : real(dp),intent(out) :: xred_sc(3,nat_sc)
790 : : real(dp),intent(out) :: magv_sc(nat_sc)
791 : : real(dp),intent(out) :: rprimd_sc(3,3)
792 : : integer,intent(in),optional :: magv_org(nat_org)
793 : :
794 : : !Local variables-------------------------------
795 : : !scalars
796 : : integer :: prtvoll,ix,iy,iz,nprcl,iprcl,jdim,iatom
797 : : !arrays
798 : 90 : real(dp) :: magvv_org(nat_org)
799 : 45 : real(dp),allocatable :: transv(:,:,:)
800 : : ! *************************************************************************
801 : :
802 [ + - ]: 45 : if (present(magv_org)) then
803 [ + + ]: 110 : magvv_org=magv_org
804 : : else
805 [ # # # # : 0 : magvv_org=(/ (1, iatom=1,nat_org) /)
# # ]
806 : : end if
807 : :
808 : : if (present(prtvol)) then
809 : : prtvoll=prtvol
810 : : else
811 : : prtvoll=1
812 : : end if
813 : :
814 [ + + + + : 990 : rprimd_sc=reshape((/ (rprimd_org(ix,:)*ext(ix) ,ix=1,3) /),(/3,3 /))
+ + ]
815 [ + + ]: 180 : nprcl=product(ext)
816 [ - + + - : 180 : ABI_MALLOC(transv,(3,nat_org,nprcl))
- + + - +
- + - - +
- + ]
817 : :
818 [ + + + + : 64355 : transv=reshape((/ (((((/ ix,iy,iz /),iatom=1,nat_org),ix=0,ext(1)-1),iy=0,ext(2)-1),iz=0,ext(3)-1) /), (/ 3, nat_org,nprcl/) )
+ + + - +
+ + - + +
+ + + - +
+ + + + +
+ - + - +
- - + - -
+ + + + +
+ ]
819 : :
820 : : !write(std_out,*)'mksupercell: xred_org ' ,xred_org
821 [ + + ]: 3365 : do iprcl=1,nprcl
822 [ + + + + ]: 21080 : xred_sc(:,1+(iprcl-1)*nat_org:iprcl*nat_org)=xred_org+transv(:,:,iprcl)
823 [ + + ]: 7805 : magv_sc(1+(iprcl-1)*nat_org:iprcl*nat_org)=magv_org
824 : : end do
825 : :
826 [ + + ]: 180 : do jdim=1,3
827 [ + + ]: 13500 : xred_sc(jdim,:)=xred_sc(jdim,:)/ext(jdim)
828 : : end do
829 : :
830 : : !write(std_out,*)'mksupercell: xred_sc ', xred_sc
831 : : !write(std_out,*)'mksupercell: magv_sc ', magv_sc
832 : :
833 [ - + ]: 45 : ABI_FREE(transv)
834 : :
835 [ + - ]: 45 : end subroutine mksupercell
836 : : !!***
837 : :
838 : : !****f* m_supercell/supercell_write_xsf
839 : : !! NAME
840 : : !! supercell_write_xsf
841 : : !!
842 : : !! FUNCTION
843 : : !! output atomic positions, supercell vectors, etc... to xsf_filenamt
844 : : !!
845 : : !! INPUTS
846 : : !! xsf_filename = filename
847 : : !!
848 : : !! OUTPUT
849 : : !! printing to file
850 : : !!
851 : : !! SOURCE
852 : :
853 : 4 : subroutine supercell_write_xsf(scell, xsf_filename)
854 : :
855 : : !Arguments ------------------------------------
856 : : !scalars
857 : : class(supercell_type), intent(in) :: scell
858 : : character(len=*), intent(in) :: xsf_filename
859 : :
860 : : !Local variables-------------------------------
861 : : integer :: ount, ix, iy, iatom
862 : : character(len=500) :: msg
863 : : ! *************************************************************************
864 : :
865 [ - + ]: 4 : if (open_file(xsf_filename, msg, newunit=ount, status="unknown", action="write") /= 0) then
866 : 0 : ABI_ERROR(msg)
867 : : end if
868 : :
869 : : ! Note: Don't put comments because Vesta on my Mac does not like them!
870 : : !write (ount, '(a)')"#", trim(title)
871 : : !write (ount, '(a,3(3I7,2x))') '# supercell rlatt is ', scell%rlatt
872 : : !write (ount, '(a,I0,a)') '# and has ', scell%ncells, ' primitive unit cells '
873 : : !write (ount, '(a)') '#'
874 : :
875 : 4 : write(ount,'(1X,A)') 'DIM-GROUP'
876 : 4 : write(ount,*) '3 1'
877 : 4 : write(ount,'(1X,A)') 'PRIMVEC'
878 : : !write(ount, "(a)")"# these are primitive lattice vectors (in Angstroms)"
879 [ + + ]: 16 : do iy = 1,3
880 [ + - + + ]: 52 : write(ount, '(3(ES17.10,2X))') (Bohr_Ang * scell%rprimd(ix,iy), ix=1,3)
881 : : end do
882 : 4 : write(ount, "(1X, a)")"PRIMCOORD"
883 : 4 : write(ount, "(i0,1x,i0)") scell%natom, 1 ! # The second number is always 1 for PRIMCOORD coordinates.
884 : :
885 [ + + ]: 4560 : do iatom=1,scell%natom
886 : : write(ount, '(i9, 6(3X,ES17.10))') &
887 : 4556 : NINT(scell%znucl(scell%typat(iatom))), & ! WARNING alchemy not supported by XCrysden
888 [ + + + + ]: 36452 : scell%xcart_ref(:,iatom) * Bohr_Ang, (scell%xcart(:,iatom) - scell%xcart_ref(:,iatom)) * Bohr_Ang
889 : : end do
890 : :
891 : 4 : close(ount)
892 : :
893 : 4 : end subroutine supercell_write_xsf
894 : : !!***
895 : :
896 [ - + - - : 27320 : end module m_supercell
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + ]
897 : : !!***
|
