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