Line data Source code
1 : !!****m* ABINIT/m_ingeo
2 : !! NAME
3 : !! m_ingeo
4 : !!
5 : !! FUNCTION
6 : !! Initialize geometry variables for the ABINIT code.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (XG, RC)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_ingeo
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_atomdata
28 : use m_sort
29 : use m_dtset
30 :
31 : use m_matrix, only : mati3inv, mati3det
32 : use m_symtk, only : chkorthsy, symrelrot, chkprimit, symatm, &
33 : & symmetrize_rprimd, symmetrize_tnons,symmetrize_xred
34 : use m_spgbuilder, only : gensymspgr, gensymshub, gensymshub4
35 : use m_symfind, only : symfind, symfind_expert, symanal, symlatt
36 : use m_geometry, only : mkradim, mkrdim, xcart2xred, xred2xcart, &
37 : & randomcellpos, metric, reduce2primitive, cart2spinaxis
38 : use m_parser, only : intagm, intagm_img, geo_t, geo_from_abivar_string, get_acell_rprim
39 : use m_numeric_tools, only : geteuler
40 :
41 : implicit none
42 :
43 : private
44 : !!***
45 :
46 : public :: ingeo ! Initialize geometry variables for the ABINIT code.
47 : public :: invacuum ! Determine whether there is vacuum along some of the primitive directions
48 : public :: checkspvec ! Check the consistency of spin-related input vectors with the spin quantization axis
49 : !!***
50 :
51 : contains
52 : !!***
53 :
54 : !!****f* m_ingeo/ingeo
55 : !!
56 : !! NAME
57 : !! ingeo
58 : !!
59 : !! FUNCTION
60 : !! Initialize geometry variables for the ABINIT code.
61 : !! 1) set up unit cell: acell, rprim and rprimd ; deduce Bravais lattice
62 : !! 2) (removed)
63 : !! 3) Set up the number of atoms (natrd) in the primitive set, to be read.
64 : !! 4) Read the type of each atom in the primitive set
65 : !! 5) Read coordinates for each atom in the primitive set
66 : !! 6) Eventually read the symmetries
67 : !! 7) Checks whether the geometry builder must be used,
68 : !! and call it if needed. Call eventually the symmetry builder and analyser
69 : !! Make the adequate transfers if the geometry
70 : !! builder is not needed.
71 : !! 8) Initialize the fixing of atoms, the initial velocities, and the initial atomic spin
72 : !!
73 : !! INPUTS
74 : !! berryopt == 4/14: electric field is on; berryopt = 6/7/16/17: electric displacement field is on
75 : !! iimage= index of the current image
76 : !! iout=unit number of output file
77 : !! jdtset=number of the dataset looked for
78 : !! lenstr=actual length of the string
79 : !! msym=default maximal number of symmetries
80 : !! natom=number of atoms
81 : !! nimage=number of images
82 : !! npsp=number of pseudopotentials (needed for the dimension of znucl)
83 : !! nspden=number of spin-density components
84 : !! nsppol=number of independent spin polarizations
85 : !! ntypat=number of type of atoms
86 : !! nzchempot=defines the use of a spatially-varying chemical potential along z
87 : !! pawspnorb=1 when spin-orbit is activated within PAW
88 : !! ratsph(1:ntypat)=radius of the atomic sphere
89 : !! string*(*)=character string containing all the input data. Initialized previously in instrng.
90 : !! supercell_latt(3)=supercell lattice
91 : !! comm: MPI communicator
92 : !!
93 : !! OUTPUT
94 : !! acell(3)=length of primitive vectors
95 : !! amu(ntypat)=mass of each atomic type
96 : !! bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
97 : !! chrgat(natom)=target charge for each atom. Not always used, it depends on the value of constraint_kind
98 : !! field_red(3)=applied field direction in reduced coordinates
99 : !! field_red_axial(3) = direction of applied axial (pseudo-)vector field
100 : !! (e.g., magnetic field) expressed in reduced coordinates
101 : !! genafm(3)=magnetic translation generator (in case of Shubnikov group type IV)
102 : !! iatfix(3,natom)=indices for atoms fixed along some (or all) directions
103 : !! jellslab=not zero if jellslab keyword is activated
104 : !! slabzbeg, slabzend= the z coordinates of beginning / end of the jellium slab
105 : !! mixalch(npspalch,ntypalch)=alchemical mixing factors
106 : !! nsym=actual number of symmetries
107 : !! nucdipmom(3,natom)=nuclear magnetic dipole moment of each atom in atomic units
108 : !! ptgroupma = magnetic point group number
109 : !! rprim(3,3)=dimensionless real space primitive translations
110 : !! spgroup=symmetry space group
111 : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
112 : !! symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
113 : !! symmorphi=if 0, only allows symmorphic symmetry operations
114 : !! symrel(3,3,1:msym)=symmetry operations in real space in terms
115 : !! of primitive translations
116 : !! tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
117 : !! tolsym=tolerance for the symmetry operations
118 : !! typat(natom)=type integer for each atom in cell
119 : !! vel(3,natom)=initial velocity of atoms in bohr/atomic time units
120 : !! vel_cell(3,3)=initial velocity of cell parameters in bohr/atomic time units
121 : !! xred(3,natom)=reduced dimensionless atomic coordinates
122 : !! znucl(1:npsp)=nuclear number of atom as specified in psp file
123 : !!
124 : !! SIDE EFFECTS
125 : !!
126 : !! NOTES
127 : !! the parameters ntypat and natom have already been read in indims,
128 : !! and were used to dimension the arrays needed here.
129 : !!
130 : !! TODO
131 : !! The dtset datastructure should NOT be an argument of this routine ... !
132 : !!
133 : !! MG: I completely agree. Abinit developers must learn that Fortran does not allow for aliasing!
134 : !!
135 : !! SOURCE
136 :
137 7233 : subroutine ingeo (acell,amu,atndlist,bravais,chrgat,dtset,field_red,field_red_axial,&
138 7233 : genafm,iatfix,iatnd,icoulomb,iimage,iout,jdtset,jellslab,lenstr,mixalch,&
139 : msym,natnd,natom,nimage,npsp,npspalch,nspden,nsym,ntypalch,ntypat,&
140 7233 : nucdipmom,nzchempot,pawspnorb,&
141 7233 : ptgroupma,ratsph,rprim,slabzbeg,slabzend,spgroup,spinat,spinat_cart,string,supercell_lattice,symafm,&
142 7233 : symmorphi,symrel,tnons,tolsym,typat,vel,vel_cell,xred,znucl,comm)
143 :
144 : !Arguments ------------------------------------
145 : !scalars
146 : integer,intent(in) :: iimage,iout,jdtset,lenstr,msym
147 : integer,intent(in) :: natnd,nimage,npsp,npspalch,nspden
148 : integer,intent(in) :: ntypalch,ntypat,nzchempot,pawspnorb,comm
149 : integer,intent(inout) :: natom,symmorphi
150 : integer,intent(out) :: icoulomb,jellslab,ptgroupma,spgroup !vz_i
151 : integer,intent(inout) :: nsym !vz_i
152 : real(dp),intent(out) :: slabzbeg,slabzend,tolsym
153 : character(len=*),intent(in) :: string
154 : !arrays
155 : integer,intent(in) :: supercell_lattice(3)
156 : integer,intent(out) :: bravais(11),iatfix(3,natom) !vz_i
157 : integer,intent(inout) :: iatnd(natnd),symafm(msym) !vz_i
158 : integer,intent(inout) :: symrel(3,3,msym) !vz_i
159 : integer,intent(out) :: typat(natom)
160 : real(dp),intent(inout) :: atndlist(3,natnd),chrgat(natom)
161 : real(dp),intent(inout) :: nucdipmom(3,natom),ratsph(ntypat)
162 : real(dp),intent(inout) :: spinat(3,natom), spinat_cart(3,natom)
163 : real(dp),intent(out) :: acell(3),amu(ntypat),field_red(3),field_red_axial(3)
164 : real(dp),intent(out) :: genafm(3),mixalch(npspalch,ntypalch)
165 : real(dp),intent(inout) :: rprim(3,3),tnons(3,msym) !vz_i
166 : real(dp),intent(out) :: vel(3,natom),vel_cell(3,3),xred(3,natom)
167 : real(dp),intent(in) :: znucl(npsp)
168 : type(dataset_type),intent(inout) :: dtset
169 :
170 : !Local variables-------------------------------
171 : character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )"
172 : character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) "
173 : !scalars
174 : integer :: bckbrvltt,brvltt,chkprim,chkprim_fake,expert_user
175 : integer :: fixed_mismatch,i1,i2,i3,iatom,iatom_supercell,idir,ierr,iexit,ii
176 : integer :: invar_z,ipsp,irreducible,isym,itranslat,itypat,jsym,marr,mismatch_fft_tnons,multi,multiplicity,natom_uc,natfix,natrd
177 : integer :: nobj,noncoll,nptsym,nsym_now,ntranslat,ntyppure,random_atpos,shubnikov,spgaxor,spgorig
178 : integer :: spgroupma,tgenafm,tnatrd,tread,try_primitive,tscalecart,tspgroupma,tread_geo,tread_cart
179 : integer :: txcart,txred,txrandom,use_inversion
180 : real(dp) :: amu_default,ucvol,sumalch
181 : character(len=1000) :: msg
182 7233 : character(len=lenstr) :: geo_string
183 : type(atomdata_t) :: atom
184 7233 : type(geo_t) :: geo
185 : !arrays
186 : integer :: bravais_reduced(11)
187 7233 : integer,allocatable :: intarr(:)
188 7233 : integer,allocatable :: is_translation(:)
189 7233 : integer,allocatable :: ptsymrel(:,:,:),typat_read(:)
190 : real(dp) :: angdeg(3), field_xred(3),gmet(3,3),gprimd(3,3),rmet(3,3),rcm(3)
191 : real(dp) :: rprimd(3,3),rprimd_read(3,3),rprimd_new(3,3),rprimd_primitive(3,3),scalecart(3)
192 7233 : real(dp),allocatable :: mass_psp(:),tnons_cart(:,:),tnons_new(:,:),translations(:,:)
193 7233 : real(dp),allocatable :: xcart(:,:),xcart_read(:,:),xred_read(:,:),dprarr(:)
194 :
195 : ! *************************************************************************
196 :
197 : !DEBUG
198 : !write(std_out,'(a)')' m_ingeo%ingeo : enter '
199 : !call flush(std_out)
200 : !ENDDEBUG
201 :
202 7233 : marr=max(12,3*natom,9*msym)
203 21699 : ABI_MALLOC(intarr,(marr))
204 21699 : ABI_MALLOC(dprarr,(marr))
205 :
206 : ! Try from geo_string
207 7233 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'structure', tread_geo, 'KEY', key_value=geo_string)
208 :
209 7233 : if (tread_geo /= 0) then
210 : ! Set up unit cell from external file.
211 55 : geo = geo_from_abivar_string(geo_string, comm)
212 220 : acell = one
213 715 : rprim = geo%rprimd
214 : !call exclude(lenstr, string, jdtset, iimage, "acell, rprim, angdeg, scalecar")
215 :
216 : else
217 : ! Set up unit cell from acell, rprim, angdeg
218 7178 : call get_acell_rprim(lenstr, string, jdtset, iimage, nimage, marr, acell, rprim)
219 : end if ! geo% or (acell, rprim, angdeg)
220 :
221 : ! Rescale rprim using scalecart (and set scalecart to one)
222 28932 : scalecart(1:3)=one
223 7233 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'scalecart',tscalecart,'LEN')
224 7263 : if(tscalecart==1) scalecart(1:3)=dprarr(1:3)
225 7233 : call intagm_img(scalecart,iimage,jdtset,lenstr,nimage,3,string,"scalecart",tscalecart,'LEN')
226 :
227 28932 : rprim(:,1)=scalecart(:)*rprim(:,1)
228 28932 : rprim(:,2)=scalecart(:)*rprim(:,2)
229 28932 : rprim(:,3)=scalecart(:)*rprim(:,3)
230 28932 : scalecart(:)=one
231 :
232 : ! Compute the multiplicity of the supercell
233 7233 : multiplicity=supercell_lattice(1)*supercell_lattice(2)*supercell_lattice(3)
234 :
235 7233 : if (tread_geo == 0) then
236 : ! Get the number of atom in the unit cell. Read natom from string
237 7178 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natom',tread,'INT')
238 :
239 : ! Might initialize natom from XYZ file
240 7178 : if (tread==0) call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'_natom',tread,'INT')
241 :
242 7178 : if(tread==1) natom_uc=intarr(1)
243 : else
244 55 : natom_uc = geo%natom
245 : end if
246 :
247 : ! Store the rprimd of the unit cell
248 7233 : call mkrdim(acell,rprim,rprimd_read)
249 :
250 : ! Multiply the rprim to get the rprim of the supercell
251 7233 : if(multiplicity > 1)then
252 4 : rprim(:,1) = rprim(:,1) * supercell_lattice(1)
253 4 : rprim(:,2) = rprim(:,2) * supercell_lattice(2)
254 4 : rprim(:,3) = rprim(:,3) * supercell_lattice(3)
255 : end if
256 :
257 : ! Compute different matrices in real and reciprocal space, also checks whether ucvol is positive.
258 7233 : call mkrdim(acell, rprim, rprimd)
259 7233 : call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
260 :
261 7233 : if (dtset%berryopt ==4) then
262 104 : do ii=1,3
263 338 : field_red(ii)=dot_product(dtset%efield(:),gprimd(:,ii))
264 : end do
265 : else if (dtset%berryopt == 6 ) then
266 0 : do ii=1,3
267 0 : field_red(ii)=dot_product(dtset%dfield(:),gprimd(:,ii))
268 0 : field_red(ii)=field_red(ii)+ dot_product(dtset%efield(:),gprimd(:,ii)) ! note: symmetry broken by D and E
269 : end do
270 : else if (dtset%berryopt == 14) then
271 0 : do ii=1,3
272 0 : field_red(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii))
273 : end do
274 : else if (dtset%berryopt == 16) then
275 0 : do ii=1,3
276 0 : field_red(ii)=dtset%red_dfield(ii)+dtset%red_efield(ii) ! symmetry broken by reduced d and e
277 : end do
278 : else if (dtset%berryopt == 17) then
279 0 : do ii=1,3
280 0 : field_red(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii))
281 0 : if(dtset%jfielddir(ii)==2) field_red(ii)=dtset%red_dfield(ii)
282 : end do
283 : end if
284 :
285 28932 : do ii = 1, 3
286 94029 : if (norm2(dtset%hspinfield_cart) > tol8) then
287 144 : field_red_axial(ii) = dot_product(dtset%hspinfield_cart(:), gprimd(:, ii))
288 : end if
289 : end do
290 :
291 : !tolsym = tol8
292 : !XG20200801 New default value for tolsym. This default value is also defined in m_invars1.F90
293 7233 : tolsym = tol5
294 : !if (tread_geo /= 0 .and. geo%filetype == "poscar") then
295 : ! tolsym = tol4
296 : ! ABI_COMMENT("Reading structure from POSCAR --> default value of tolsym is set to 1e-4")
297 : !end if
298 :
299 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'tolsym',tread,'DPR')
300 7233 : if(tread==1) tolsym=dprarr(1)
301 :
302 : ! Find a tentative Bravais lattice and its point symmetries (might not use them)
303 : ! Note that the Bravais lattice might not be the correct one yet (because the
304 : ! actual atomic locations might lower the symmetry obtained from the lattice parameters only)
305 21699 : ABI_MALLOC(ptsymrel,(3,3,msym))
306 7233 : call symlatt(bravais,dev_null,msym,nptsym,ptsymrel,rprimd,tolsym)
307 :
308 : ! 3) Possibly, initialize a jellium slab
309 7233 : jellslab=0
310 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'jellslab',tread,'INT')
311 7233 : if(tread==1) jellslab=intarr(1)
312 :
313 7233 : slabzbeg=zero
314 7233 : slabzend=zero
315 7233 : if(jellslab/=0)then
316 6 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slabzbeg',tread,'DPR')
317 6 : if(tread==1) slabzbeg=dprarr(1)
318 :
319 6 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'slabzend',tread,'DPR')
320 6 : if(tread==1) slabzend=dprarr(1)
321 : end if
322 :
323 : ! 4) Set up the number of atoms in the primitive set, to be read.
324 : ! This is the default
325 7233 : natrd=natom
326 7233 : if(multiplicity > 1) natrd = natom_uc
327 :
328 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natrd',tnatrd,'INT')
329 7233 : if(tnatrd==1) natrd=intarr(1)
330 :
331 7233 : if(natrd<1 .or. natrd>natom)then
332 0 : if(natrd>1 .and. multiplicity > 1) then
333 0 : if(natrd < natom)then
334 : write(msg, '(3a)' )&
335 0 : 'The number of atoms to be read (natrd) can not be used with supercell_latt.',ch10,&
336 0 : 'Action: Remove natrd or supercell_latt in your input file.'
337 0 : ABI_ERROR(msg)
338 : else
339 : write(msg,'(3a,I0,a,I0,a,I0,2a)')&
340 0 : 'The input variable supercell_latt is present',ch10,&
341 0 : 'thus a supercell ',supercell_lattice(1),' ',supercell_lattice(2),&
342 0 : ' ',supercell_lattice(3),' is generated',ch10
343 0 : ABI_WARNING(msg)
344 : end if
345 : else
346 : write(msg, '(3a,i0,a,i0,2a,a)' )&
347 0 : 'The number of atoms to be read (natrd) must be positive and not bigger than natom.',ch10,&
348 0 : 'This is not the case: natrd=',natrd,', natom=',natom,ch10,&
349 0 : 'Action: correct natrd or natom in your input file.'
350 0 : ABI_ERROR(msg)
351 : end if
352 : end if
353 :
354 : ! 5) Read the type and initial spin of each atom in the primitive set--------
355 21699 : ABI_MALLOC(typat_read,(natrd))
356 7233 : typat_read(1)=1
357 :
358 7233 : if (tread_geo == 0) then
359 7178 : call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'typat',tread,'INT')
360 :
361 : ! If not read, try the XYZ data
362 7178 : if(tread==0) call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'_typat',tread,'INT')
363 23224 : if(tread==1) typat_read(1:natrd)=intarr(1:natrd)
364 :
365 : else
366 385 : typat_read = geo%typat
367 : end if
368 :
369 23578 : do iatom=1,natrd
370 23578 : if(typat_read(iatom)<1 .or. typat_read(iatom)>ntypat )then
371 : write(msg,'(a,i0,a,i0,a,a,a,i0,a,a,a)')&
372 0 : 'The input type of atom number ',iatom,' is equal to ',typat_read(iatom),',',ch10,&
373 0 : 'while it should be between 1 and ntypat= ',ntypat,'.',ch10,&
374 0 : 'Action: change either the variable typat or the variable ntypat.'
375 0 : ABI_ERROR(msg)
376 : end if
377 : end do
378 :
379 : ! 6) Read coordinates for each atom in the primitive set--------
380 :
381 21699 : ABI_MALLOC(xcart_read,(3,natrd))
382 14466 : ABI_MALLOC(xred_read,(3,natrd))
383 :
384 7233 : random_atpos=0
385 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'random_atpos',txrandom,'INT')
386 7233 : if(txrandom==1) random_atpos=intarr(1)
387 7233 : if (random_atpos < 0 .or. random_atpos > 5) then
388 : write(msg,'(3a)')&
389 0 : 'Random positions is a variable defined between 0 and 5. Error in the input file. ',ch10,&
390 0 : 'Action: define one of these in your input file.'
391 0 : ABI_ERROR(msg)
392 : end if
393 :
394 : !if(nimage/=1 .and. iimage/=1)then
395 : !FIXME: should this be called outside the above end if?
396 : call randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd_read,typat_read,&
397 7233 : xred_read(:,1:natrd),znucl,acell)
398 : !This should not be printed if randomcellpos did nothing - it contains garbage. Spurious output anyway
399 : !end if
400 :
401 7233 : if (tread_geo /= 0) then
402 55 : txcart = 0; txrandom = 0; txred = 1
403 1210 : xred_read = geo%xred
404 :
405 : else
406 7178 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xred',txred,'DPR')
407 19822 : if (txred==1 .and. txrandom == 0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd) , [3, natrd])
408 7178 : call intagm_img(xred_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xred",txred,'DPR')
409 :
410 7178 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xcart',txcart,'LEN')
411 8868 : if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
412 7178 : call intagm_img(xcart_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xcart",txcart,'LEN')
413 :
414 : ! Might initialize xred from XYZ file
415 7178 : if (txred+txcart+txrandom==0) then
416 1 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'_xred',txred,'DPR')
417 1 : if (txred==1 .and. txrandom==0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
418 :
419 1 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'_xcart',txcart,'DPR')
420 3 : if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
421 : end if
422 :
423 : !TODO: Might initialize xred from getxred/xcart: NOT POSSIBLE YET. NEEDS INTER DTSET COMMUNICATION AT INVARS1 TIME
424 : ! if (txred+txcart+txrandom==0) then
425 : ! call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'getxred',txred,'INT')
426 : ! if (txred==1 .and. txrandom==0) xred_read(:,1:natrd) =
427 : !
428 : ! call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'getxcart',txcart,'INT')
429 : ! if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) =
430 : ! end if
431 : end if
432 :
433 7233 : if (txred + txcart + txrandom == 0) then
434 : write(msg, '(3a)' )&
435 0 : 'Neither xred nor xcart are present in input file. ',ch10,&
436 0 : 'Action: define one of these in your input file.'
437 0 : ABI_ERROR(msg)
438 : end if
439 :
440 7233 : if (txred==1) write(msg, '(a)' ) ' xred is defined in input file'
441 7233 : if (txcart ==1) write(msg, '(a)' ) ' xcart is defined in input file (possibly in Angstrom)'
442 7233 : if (txrandom ==1) write(msg, '(a)' ) ' xred as random positions in the unit cell'
443 7233 : if (txrandom ==1) write(msg, '(a)' ) ' xcart are defined from a random distribution '
444 7233 : call wrtout(std_out, msg)
445 :
446 7233 : if (txred + txcart + txrandom > 1)then
447 : write(msg, '(3a)' )&
448 0 : 'Too many input channels for atomic positions are defined.',ch10,&
449 0 : 'Action: choose to define only one of these.'
450 0 : ABI_ERROR(msg)
451 : end if
452 :
453 7233 : if (txred==1 .or. txrandom /=0 ) then
454 6387 : call wrtout(std_out,' ingeo: takes atomic coordinates from input array xred ')
455 6387 : call xred2xcart(natrd,rprimd_read,xcart_read,xred_read)
456 : else
457 846 : call wrtout(std_out,' ingeo: takes atomic coordinates from input array xcart')
458 846 : txcart=1
459 : end if
460 :
461 : !At this stage, the cartesian coordinates are known, for the atoms whose coordinates where read.
462 :
463 : ! Here, allocate the variable that will contain the completed
464 : ! sets of xcart, after the use of the geometry builder or the symmetry builder
465 21699 : ABI_MALLOC(xcart,(3,natom))
466 :
467 : !7) Eventually read the symmetries
468 : !Take care of the symmetries
469 :
470 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nsym',tread,'INT')
471 7233 : if(tread==1) nsym=intarr(1)
472 :
473 : ! Check that nsym is not negative
474 7233 : if (nsym<0) then
475 : write(msg, '(a,i0,4a)' )&
476 0 : 'Input nsym must be positive or 0, but was ',nsym,ch10,&
477 0 : 'This is not allowed.',ch10,'Action: correct nsym in your input file.'
478 0 : ABI_ERROR(msg)
479 : end if
480 : ! Check that nsym is not bigger than msym
481 7233 : if (nsym>msym) then
482 : write(msg, '(2(a,i0),5a)')&
483 0 : 'Input nsym = ',nsym,' exceeds msym = ',msym,'.',ch10,&
484 0 : 'This is not allowed.',ch10,'Action: correct nsym in your input file.'
485 0 : ABI_ERROR(msg)
486 : end if
487 7233 : if (multiplicity>1) then
488 1 : nsym = 1
489 1 : ABI_WARNING('Input nsym is now set to one due to the supercell_latt input')
490 : end if
491 :
492 : ! Read symmorphi
493 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'symmorphi',tread,'INT')
494 7233 : if(tread==1) symmorphi=intarr(1)
495 :
496 : ! Now, read the symmetry operations
497 7233 : if(nsym>0)then
498 1270 : call intagm(dprarr,intarr,jdtset,marr,9*nsym,string(1:lenstr),'symrel',tread,'INT')
499 1270 : if(nsym>1 .and. tread==0)then
500 : write(msg,'(3a)')&
501 0 : 'When nsym>1, symrel must be defined in the input file.',ch10,&
502 0 : 'Action: either change nsym, or define symrel in your input file.'
503 0 : ABI_ERROR(msg)
504 : end if
505 3106 : if(tread==1) symrel(:,:,1:nsym)=reshape( intarr(1:9*nsym) , [3, 3, nsym])
506 :
507 : ! Take care of tnons
508 41726 : tnons(:,1:nsym)=zero
509 1270 : call intagm(dprarr,intarr,jdtset,marr,3*nsym,string(1:lenstr),'tnons',tread,'DPR')
510 2348 : if(tread==1) tnons(:,1:nsym)=reshape( dprarr(1:3*nsym), [3, nsym])
511 :
512 1270 : if(symmorphi==0)then
513 58 : do isym=1,nsym
514 145 : if(sum(tnons(:,isym)**2)>tol6)then
515 : write(msg, '(5a,i0,a,3f8.4,3a)' )&
516 0 : 'When symmorphi /= 1, the vectors of translation (tnons)',ch10,&
517 0 : 'a symmetry operation must vanish.',ch10,&
518 0 : 'However, for the symmetry operation number ',isym,', tnons =',tnons(:,isym),'.',ch10,&
519 0 : 'Action: either change your list of allowed symmetry operations, or use the symmetry finder (nsym=0).'
520 0 : ABI_ERROR(msg)
521 : end if
522 : end do
523 : end if
524 :
525 : ! Take care of symafm
526 1270 : call intagm(dprarr,intarr,jdtset,marr,nsym,string(1:lenstr),'symafm',tread,'INT')
527 1318 : if(tread==1) symafm(1:nsym)=intarr(1:nsym)
528 : end if
529 :
530 :
531 : !8) Checks whether the geometry builder must be used, and call it if needed.
532 : !Call the symmetry builder and analyzer if needed.
533 :
534 : ! At this stage, nsym might still contain the default 0, msym contains the default dtset%maxnsym.
535 : ! The cartesian coordinates of the atoms of the primitive set are contained in xcart_read.
536 :
537 7233 : nobj=0
538 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nobj',tread,'INT')
539 7233 : if(tread==1) nobj=intarr(1)
540 7233 : if(nobj /= 0 .and. multiplicity > 1)then
541 : write(msg, '(3a)' )&
542 0 : 'nobj can not be used with supercell_latt.',ch10,&
543 0 : 'Action: Remove nobj or supercell_latt in your input file.'
544 0 : ABI_ERROR(msg)
545 : end if
546 :
547 : !If there are objects, chkprim will not be used immediately
548 : !But, if there are no objects, but a space group, it will be used directly.
549 : !Need first to check the value of expert_user
550 7233 : expert_user=0
551 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'expert_user',tread,'INT')
552 7233 : if(tread==1) expert_user=intarr(1)
553 25 : if(expert_user==0)then
554 7208 : chkprim=1
555 7208 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'chkprim',tread,'INT')
556 7208 : if(tread==1) chkprim=intarr(1)
557 : else
558 25 : chkprim=0
559 : endif
560 :
561 7233 : if(nobj/=0)then
562 :
563 : ! chrgat is read for each atom, from 1 to natom
564 4 : call intagm(dprarr,intarr,jdtset,marr,natom,string(1:lenstr),'chrgat',tread,'DPR')
565 4 : if(tread==1) then
566 0 : chrgat(1:natom) = dprarr(1:natom)
567 : end if
568 :
569 : ! Spinat is read for each atom, from 1 to natom
570 4 : call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'spinat',tread,'DPR')
571 4 : if(tread==1) then
572 0 : spinat(1:3,1:natom) = reshape( dprarr(1:3*natom) , [3, natom])
573 0 : dtset%spinat_in(1:3,1:natom) = spinat(1:3,1:natom)
574 : end if
575 :
576 4 : call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'spinat_cart',tread_cart,'DPR')
577 4 : if(tread_cart==1) spinat_cart(1:3,1:natom) = reshape( dprarr(1:3*natom) , [3, natom])
578 :
579 4 : call checkspvec('spinat',natom,dtset%spinaxis,tread,tread_cart,spinat,spinat_cart)
580 4 : if (tread == 0 .and. tread_cart == 1) dtset%spinat_in(1:3,1:natom) = spinat(1:3,1:natom)
581 :
582 : ! nucdipmom is read for each irreducible atom, from 1 to natom
583 1360 : nucdipmom=zero
584 4 : if(natnd > 0) then
585 0 : call intagm(dprarr,intarr,jdtset,marr,natnd,string(1:lenstr),'iatnd',tread,'INT')
586 0 : if(tread==1) iatnd(1:natnd)=intarr(1:natnd)
587 0 : call intagm(dprarr,intarr,jdtset,marr,3*natnd,string(1:lenstr),'atndlist',tread,'DPR')
588 0 : if(tread==1) atndlist(1:3,1:natnd)=reshape(dprarr(1:3*natnd),[3,natnd])
589 0 : do ii=1,natnd
590 0 : nucdipmom(1:3,iatnd(ii))=atndlist(1:3,ii)
591 : end do
592 : else
593 4 : call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'nucdipmom',tread,'DPR')
594 4 : if(tread==1)nucdipmom(1:3,1:natom) = reshape( dprarr(1:3*natom) , [3, natom])
595 : end if
596 :
597 : ! Will use the geometry builder
598 4 : if(tnatrd/=1 .and. nobj/=0)then
599 : write(msg, '(3a,i0,5a)' )&
600 0 : 'The number of atoms to be read (natrd) must be initialized',ch10,&
601 0 : 'in the input file, when nobj= ',nobj,'.',ch10,&
602 0 : 'This is not the case.',ch10,&
603 0 : 'Action: initialize natrd in your input file.'
604 0 : ABI_ERROR(msg)
605 : end if
606 :
607 4 : if(jellslab/=0)then
608 : write(msg, '(a,i0,3a)' )&
609 0 : 'A jellium slab cannot be used when nobj= ',nobj,'.',ch10,&
610 0 : 'Action: change one of the input variables jellslab or nobj in your input file.'
611 0 : ABI_ERROR(msg)
612 : end if
613 :
614 4 : call ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read)
615 :
616 : ! Finalize the computation of coordinates: produce xred.
617 4 : call xcart2xred(natom,rprimd,xcart,xred)
618 :
619 : else
620 : ! nobj==0
621 :
622 : ! chrgat is read for each irreducible atom, from 1 to natrd
623 7229 : call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'chrgat',tread,'DPR')
624 7291 : if(tread==1)chrgat(1:natrd) = dprarr(1:natrd)
625 :
626 : ! Spinat is read for each irreducible atom, from 1 to natrd
627 7229 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'spinat',tread,'DPR')
628 7229 : if(tread==1) then
629 5082 : spinat(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd])
630 10190 : dtset%spinat_in(1:3,1:natrd) = spinat(1:3,1:natrd)
631 : end if
632 :
633 7229 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'spinat_cart',tread_cart,'DPR')
634 7231 : if(tread_cart==1) spinat_cart(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natom])
635 7229 : call checkspvec('spinat',natrd,dtset%spinaxis,tread,tread_cart,spinat,spinat_cart)
636 7233 : if (tread == 0 .and. tread_cart == 1) dtset%spinat_in(1:3,1:natrd) = spinat(1:3,1:natrd)
637 :
638 : ! nucdipmom is read for each irreducible atom, from 1 to natrd
639 160877 : nucdipmom=zero
640 7229 : if(natnd > 0) then
641 14 : call intagm(dprarr,intarr,jdtset,marr,natnd,string(1:lenstr),'iatnd',tread,'INT')
642 32 : if(tread==1) iatnd(1:natnd)=intarr(1:natnd)
643 14 : call intagm(dprarr,intarr,jdtset,marr,3*natnd,string(1:lenstr),'atndlist',tread,'DPR')
644 42 : if(tread==1) atndlist(1:3,1:natnd)=reshape(dprarr(1:3*natnd),[3,natnd])
645 32 : do ii=1,natnd
646 86 : nucdipmom(1:3,iatnd(ii))=atndlist(1:3,ii)
647 : end do
648 : else
649 7215 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'nucdipmom',tread,'DPR')
650 7269 : if(tread==1)nucdipmom(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd])
651 : end if
652 :
653 : ! Compute xred/typat and spinat for the supercell
654 7229 : if(multiplicity > 1)then
655 1 : iatom_supercell = 0
656 3 : do i1 = 1, supercell_lattice(1)
657 7 : do i2 = 1, supercell_lattice(2)
658 14 : do i3 = 1, supercell_lattice(3)
659 52 : do iatom = 1, natom_uc
660 40 : iatom_supercell = iatom_supercell + 1
661 880 : xcart(:,iatom_supercell) = xcart_read(:,iatom) + matmul(rprimd_read,(/i1-1,i2-1,i3-1/))
662 40 : chrgat(iatom_supercell) = chrgat(iatom)
663 160 : spinat(1:3,iatom_supercell) = spinat(1:3,iatom)
664 160 : dtset%spinat_in(1:3,iatom_supercell) = dtset%spinat_in(1:3,iatom)
665 48 : typat(iatom_supercell) = typat_read(iatom)
666 : end do
667 : end do
668 : end do
669 : end do
670 1 : call xcart2xred(natom,rprimd,xcart,xred)
671 : else
672 : ! No supercell
673 7228 : call xcart2xred(natrd,rprimd,xcart_read,xred)
674 : end if
675 :
676 7229 : spgroup=0
677 7229 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroup',tread,'INT')
678 7229 : if(tread==1) spgroup=intarr(1)
679 :
680 7229 : if(spgroup/=0 .or. nsym/=0)then
681 :
682 2889 : if(jellslab/=0 .and. nsym/=1 .and. spgroup/=1)then
683 : write(msg, '(5a)' )&
684 0 : 'For the time being, a jellium slab can only be used',ch10,&
685 0 : 'either with the symmetry finder (nsym=0) or with the space group 1 (nsym=1)',ch10,&
686 0 : 'Action: change one of the input variables jellslab or nsym or spgroup in your input file.'
687 0 : ABI_ERROR(msg)
688 : end if
689 :
690 2889 : if(nzchempot/=0 .and. nsym/=1 .and. spgroup/=1)then
691 : write(msg, '(5a)' )&
692 0 : 'For the time being, a spatially-varying chemical potential can only be used',ch10,&
693 0 : 'either with the symmetry finder (nsym=0) or with the space group 1 (nsym=1)',ch10,&
694 0 : 'Action: change one of the input variables nzchempot or nsym or spgroup in your input file.'
695 0 : ABI_ERROR(msg)
696 : end if
697 :
698 7586 : typat(1:natrd)=typat_read(1:natrd)
699 :
700 2889 : if(spgroup/=0 .and. nsym/=0)then
701 : write(msg, '(a,i0,a,a,i0,a,a,a,a,a,a,a,a)' )&
702 0 : 'The spatial group number spgroup= ',spgroup,ch10,&
703 0 : 'is specified, as well as the number of symmetries nsym= ',nsym,ch10,&
704 0 : 'This is not allowed, as you can define the symmetries',ch10,&
705 0 : 'either using spgroup OR using nsym, but not both.',ch10,&
706 0 : 'Action: modify your input file',ch10,&
707 0 : '(either set spgroup to 0, or nsym to 0)'
708 0 : ABI_ERROR(msg)
709 : end if
710 :
711 2889 : brvltt=0
712 :
713 2889 : if(spgroup/=0)then
714 :
715 : ! Will generate the spatial group using spgroup
716 : ! Assign default values
717 1623 : spgaxor=1
718 1623 : spgorig=1
719 1623 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'brvltt',tread,'INT')
720 1623 : if(tread==1) brvltt=intarr(1)
721 1623 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgaxor',tread,'INT')
722 1623 : if(tread==1) spgaxor=intarr(1)
723 1623 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgorig',tread,'INT')
724 1623 : if(tread==1) spgorig=intarr(1)
725 :
726 : ! Treat the case of magnetic groups
727 1623 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroupma',tspgroupma,'INT')
728 1623 : if(tspgroupma==1) spgroupma=intarr(1)
729 1623 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'genafm',tgenafm,'DPR')
730 1626 : if(tgenafm==1) genafm(1:3)=dprarr(1:3)
731 1623 : if(tspgroupma/=0 .and. tgenafm/=0)then
732 : write(msg, '(a,i0,a,a,3es9.2,a,a,a,a,a,a,a,a)' )&
733 0 : 'The spatial group number spgroupma= ',spgroupma,ch10,&
734 0 : 'is specified, as well as the antiferromagnetic generator genafm=',genafm(1:3),ch10,&
735 0 : 'This is not allowed, as you can define the magnetic space group',ch10,&
736 0 : 'either using spgroupma OR using genafm, but not both.',ch10,&
737 0 : 'Action: modify your input file',ch10,&
738 0 : '(either define spgroupma or genafm)'
739 0 : ABI_ERROR(msg)
740 : end if
741 :
742 : ! TODO: all the symmetry generation operations should be in one big routine
743 :
744 : ! If spgroupma is defined, check whether it is consistent
745 : ! with spgroup, determine the Shubnikov type,
746 : ! and, for type IV, find the corresponding genafm
747 1623 : shubnikov=1
748 1623 : if(tspgroupma==1)then
749 1192 : call gensymshub(genafm,spgroup,spgroupma,shubnikov)
750 431 : else if(tgenafm==1)then
751 1 : shubnikov=4
752 : end if
753 :
754 : ! Generate the spatial group of symmetries in a conventional cell
755 : ! In case of Shubnikov space group type IV, only generate the
756 : ! Fedorov (non-magnetic) group. For Shubnikov type III space group,
757 : ! the magnetic part is generated here.
758 1623 : bckbrvltt=brvltt
759 1623 : if(brvltt==-1)brvltt=0
760 1623 : call gensymspgr(brvltt,msym,nsym,shubnikov,spgaxor,spgorig,spgroup,spgroupma,symafm,symrel,tnons)
761 :
762 : ! For shubnikov type IV groups,
763 : ! double the space group, using the antiferromagnetic translation generator
764 1623 : if(shubnikov==4)then
765 519 : call gensymshub4(genafm,msym,nsym,symafm,symrel,tnons)
766 : end if
767 :
768 : !write(std_out,*)' after gensymshub4, nsym =',nsym
769 : !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
770 : !do ii=1,nsym
771 : !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
772 : !end do
773 :
774 : ! If brvltt was -1 at input, one should now change the conventional cell
775 : ! to a primitive one, if brvltt/=1
776 1623 : if(bckbrvltt==-1 .and. brvltt/=1)then
777 : ! Will work with rprim only
778 384 : rprim(:,:)=rprimd(:,:)
779 384 : rprimd_new(:,:)=rprimd(:,:)
780 1536 : acell(:)=1.0_dp
781 : select case(brvltt)
782 : case(5)
783 112 : rprimd_new(:,2)=(rprim(:,2)+rprim(:,3))*0.5_dp
784 112 : rprimd_new(:,3)=(rprim(:,3)-rprim(:,2))*0.5_dp
785 : case(6)
786 0 : rprimd_new(:,1)=(rprim(:,1)+rprim(:,3))*0.5_dp
787 0 : rprimd_new(:,3)=(rprim(:,3)-rprim(:,1))*0.5_dp
788 : case(4)
789 456 : rprimd_new(:,1)=(rprim(:,1)+rprim(:,2))*0.5_dp
790 456 : rprimd_new(:,2)=(rprim(:,2)-rprim(:,1))*0.5_dp
791 : case(3)
792 248 : rprimd_new(:,1)=(rprim(:,2)+rprim(:,3))*0.5_dp
793 248 : rprimd_new(:,2)=(rprim(:,1)+rprim(:,3))*0.5_dp
794 248 : rprimd_new(:,3)=(rprim(:,1)+rprim(:,2))*0.5_dp
795 : case(2)
796 720 : rprimd_new(:,1)=(-rprim(:,1)+rprim(:,2)+rprim(:,3))*0.5_dp
797 720 : rprimd_new(:,2)=( rprim(:,1)-rprim(:,2)+rprim(:,3))*0.5_dp
798 720 : rprimd_new(:,3)=( rprim(:,1)+rprim(:,2)-rprim(:,3))*0.5_dp
799 : case(7)
800 0 : rprimd_new(:,1)=( rprim(:,1)*2.0_dp+rprim(:,2)+rprim(:,3))/3.0_dp
801 0 : rprimd_new(:,2)=(-rprim(:,1) +rprim(:,2)+rprim(:,3))/3.0_dp
802 384 : rprimd_new(:,3)=(-rprim(:,1)-rprim(:,2)*2.0_dp+rprim(:,3))/3.0_dp
803 : end select
804 384 : call symrelrot(nsym,rprimd,rprimd_new,symrel,tolsym)
805 : ! Produce xred in the new system of coordinates
806 384 : call xred2xcart(natrd,rprimd,xcart,xred)
807 384 : call xcart2xred(natrd,rprimd_new,xcart,xred)
808 : ! Produce tnons in the new system of coordinates
809 1152 : ABI_MALLOC(tnons_cart,(3,nsym))
810 384 : call xred2xcart(nsym,rprimd,tnons_cart,tnons)
811 384 : call xcart2xred(nsym,rprimd_new,tnons_cart,tnons)
812 384 : ABI_FREE(tnons_cart)
813 :
814 : ! write(std_out,*)' after change of coordinates, nsym =',nsym
815 : ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
816 : ! do ii=1,nsym
817 : ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
818 : ! end do
819 :
820 : ! Prune the symmetry operations: suppress those with
821 : ! exactly the same point and magnetic part
822 384 : nsym_now=1
823 14784 : do isym=2,nsym
824 263592 : irreducible=1
825 263592 : do jsym=1,nsym_now
826 3366168 : if(sum(abs(symrel(:,:,isym)-symrel(:,:,jsym)))==0 .and. symafm(isym)==symafm(jsym)) then
827 : irreducible=0
828 : exit
829 : end if
830 : end do
831 14784 : if(irreducible==1)then
832 5044 : nsym_now=nsym_now+1
833 65572 : symrel(:,:,nsym_now)=symrel(:,:,isym)
834 20176 : tnons(:,nsym_now)=tnons(:,isym)
835 5044 : symafm(nsym_now)=symafm(isym)
836 : end if
837 : end do
838 384 : nsym=nsym_now
839 : ! Translate tnons in the ]-0.5,0.5] interval
840 22096 : tnons(:,1:nsym)=tnons(:,1:nsym)-nint(tnons(:,1:nsym)-1.0d-8)
841 :
842 : ! DEBUG
843 : ! write(std_out,*)' after reduction, nsym =',nsym
844 : ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
845 : ! do ii=1,nsym
846 : ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
847 : ! end do
848 : ! ENDDEBUG
849 :
850 : ! Now that symrel, tnons and xred are expressed in the primitive
851 : ! axis system, update the geometric quantities
852 384 : rprimd(:,:)=rprimd_new(:,:)
853 384 : rprim(:,:)=rprimd_new(:,:)
854 384 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
855 384 : call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym)
856 : end if
857 :
858 : end if
859 :
860 2889 : if(natom/=natrd.and.multiplicity == 1)then
861 : ! Generate the full set of atoms from its knowledge in the irreducible part.
862 1588 : call fillcell(chrgat,natom,natrd,nsym,nucdipmom,spinat,spinat_cart,symafm,symrel,tnons,tolsym,typat,xred)
863 : ! Keep spinat_in consistent with symmetry-expanded spinat for output/printing.
864 98048 : dtset%spinat_in(1:3,1:natom) = spinat(1:3,1:natom)
865 : end if
866 :
867 : ! Check whether the symmetry operations are consistent with the lattice vectors
868 2889 : iexit=0
869 :
870 2889 : call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym)
871 :
872 : else
873 : ! spgroup==0 and nsym==0
874 :
875 : ! Here, spgroup==0 as well as nsym==0, so must generate
876 : ! the spatial group of symmetry. However, all the atom
877 : ! positions must be known, so the number
878 : ! of atoms to be read must equal the total number of atoms.
879 4340 : if(natrd/=natom .and. multiplicity== 1)then
880 : write(msg, '(a,i0,a,a,i0,a,a,a,a,a,a,a,a,a)' )&
881 0 : 'The number of atoms to be read (natrd)= ',natrd,ch10,&
882 0 : 'differs from the total number of atoms (natom)= ',natom,ch10,&
883 0 : 'while spgroup=0 and nsym=0.',&
884 0 : 'This is not allowed, since the information needed to',ch10,&
885 0 : 'generate the missing atomic coordinates is not available.',ch10,&
886 0 : 'Action: modify your input file',ch10,&
887 0 : '(either natrd, or natom, or spgroup, or nsym)'
888 0 : ABI_ERROR(msg)
889 : endif
890 15899 : if (multiplicity==1) typat(:)=typat_read(:)
891 :
892 : ! Find the symmetry operations: nsym, symafm, symrel and tnons.
893 : ! Use nptsym and ptsymrel, as determined by symlatt
894 : ! Will possibly correct xred and tnons.
895 :
896 4340 : noncoll=0; if (nspden == 4) noncoll=1
897 :
898 4340 : use_inversion=1
899 4340 : if (dtset%usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then
900 135 : ABI_COMMENT("Removing inversion and improper rotations from initial space group because of PAW + SOC")
901 : ! MMignolet: PAW can be used with inversion, however it results in seg faults in the dmft code. To enable when this is fixed...
902 135 : use_inversion=0
903 : end if
904 :
905 : ! Get field in reduced coordinates (reduced e/d field)
906 :
907 : field_xred(:)=zero
908 : if (dtset%berryopt ==4) then
909 : do ii=1,3
910 : field_xred(ii)=dot_product(dtset%efield(:),gprimd(:,ii))
911 : end do
912 : else if (dtset%berryopt == 6 ) then
913 : do ii=1,3
914 : field_xred(ii)=dot_product(dtset%dfield(:),gprimd(:,ii))
915 : field_xred(ii)=field_xred(ii)+ dot_product(dtset%efield(:),gprimd(:,ii)) ! note: symmetry broken by D and E
916 : end do
917 : else if (dtset%berryopt == 14) then
918 : do ii=1,3
919 : field_xred(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii))
920 : end do
921 : else if (dtset%berryopt == 16) then
922 : do ii=1,3
923 : field_xred(ii)=dtset%red_dfield(ii)+dtset%red_efield(ii) ! symmetry broken by reduced d and e
924 : end do
925 : else if (dtset%berryopt == 17) then
926 : do ii=1,3
927 : field_xred(ii)=dot_product(dtset%red_efieldbar(:),gmet(:,ii))
928 : if(dtset%jfielddir(ii)==2) field_xred(ii)=dtset%red_dfield(ii)
929 : end do
930 : end if
931 :
932 : ! Loop on trials to generate better point symmetries by relying on a primitive cell instead (possibly) of a non-primitive one,
933 : ! This loop has been disactivated, because it is not clear that one can generate a more complete set of point symmetries
934 : ! WITH INTEGER components of symrel from a primitive cell. One should allow non-integer components, but this would
935 : ! be a large departure from the current implementation. Still, the detection of the existence of the primitive cell
936 : ! and the corresponding Bravais lattice is activated.
937 4382 : do try_primitive=1,1
938 :
939 4340 : invar_z=0 ; if(jellslab/=0 .or. nzchempot/=0)invar_z=2
940 : call symfind_expert(gprimd,msym,natom,nptsym,nspden,nsym,&
941 : pawspnorb,dtset%prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,dtset%usepaw,xred,&
942 4340 : chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=dtset%field_red,invaraxial_red=dtset%field_red_axial,invar_z=invar_z)
943 :
944 4340 : chkprim_fake=-1
945 13020 : ABI_MALLOC(is_translation,(nsym))
946 4340 : call chkprimit(chkprim_fake, multi, nsym, symafm, symrel, is_translation)
947 :
948 4382 : if(multi/=1)then ! The cell is not primitive, get the point symmetries from a primitive cell.
949 117 : ntranslat=multi
950 351 : ABI_MALLOC(translations,(3,ntranslat))
951 117 : itranslat=0
952 17301 : do isym=1,nsym
953 17301 : if(is_translation(isym)==1)then
954 414 : itranslat=itranslat+1
955 1656 : translations(:,itranslat)=tnons(:,isym)
956 : endif
957 : enddo
958 117 : ABI_FREE(is_translation)
959 117 : call reduce2primitive(ntranslat, rprimd, rprimd_primitive, tolsym, translations)
960 117 : ABI_FREE(translations)
961 : !Find the Bravais lattice of the primitive cell, and the point symmetries (however, in the primitive basis)
962 117 : call symlatt(bravais_reduced,dev_null,msym,nptsym,ptsymrel,rprimd_primitive,tolsym)
963 : write(msg,'(2a,3(3es16.8,a),2(a,i4,a),3(a,3i4,a),a,i4)')&
964 117 : & ' The cell is not primitive. One could obtain a primitive cell using the following primitive vectors (rprimd) :',ch10,&
965 117 : & rprimd_primitive(1:3,1),ch10,&
966 117 : & rprimd_primitive(1:3,2),ch10,&
967 117 : & rprimd_primitive(1:3,3),ch10,&
968 117 : & ' This Bravais lattice has iholohedry =',bravais(1),ch10,&
969 117 : & ' center =',bravais(2),ch10,&
970 117 : & ' bravais(3:5) =',bravais(3:5),ch10,&
971 117 : & ' bravais(6:8) =',bravais(6:8),ch10,&
972 117 : & ' bravais(9:11)=',bravais(9:11),ch10,&
973 234 : & ' The number of point symmetries would be nptsym=',nptsym
974 117 : ABI_COMMENT(msg)
975 :
976 : !Convert the point symmetries to the non-primitive reduced coordinates
977 117 : call symrelrot(nsym, rprimd_primitive, rprimd, ptsymrel, tolsym, ierr)
978 : !Perhaps not all components of symrel are integers. This generates a return code, and precludes upgrading ptsymrel.
979 117 : if(ierr/=0)then
980 : write(msg,'(a)')&
981 75 : & ' Not all components of symrel are integers in the primitive cell coordinate system.'
982 75 : ABI_COMMENT(msg)
983 75 : exit
984 : endif
985 : else ! The cell is primitive
986 4223 : ABI_FREE(is_translation)
987 4223 : exit
988 : endif
989 :
990 : enddo ! try_primitive
991 :
992 : end if ! spgroup==0 and nsym==0
993 :
994 : ! Finalize the computation of coordinates: produce xcart
995 7229 : call xred2xcart(natom,rprimd,xcart,xred)
996 :
997 : end if ! check of existence of an object
998 :
999 7233 : ABI_FREE(ptsymrel)
1000 7233 : ABI_FREE(xcart_read)
1001 7233 : ABI_FREE(xcart)
1002 7233 : ABI_FREE(xred_read)
1003 7233 : ABI_FREE(typat_read)
1004 :
1005 7233 : call geo%free()
1006 :
1007 : ! Correct the default nsym value, if a symmetry group has not been generated.
1008 7233 : if (nsym==0) nsym=1
1009 :
1010 : !--------------------------------------------------------------------------------------------------------
1011 :
1012 7233 : icoulomb=0
1013 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'icoulomb',tread,'INT')
1014 7233 : if(tread==1)icoulomb=intarr(1)
1015 :
1016 : ! calculate the center of the atomic system such as to put the
1017 : ! atoms in the middle of the simulation box for the free BC case.
1018 7233 : if (icoulomb == 1) then
1019 0 : rcm(:)=zero
1020 0 : do iatom=1,natom
1021 0 : rcm(:)=rcm(:)+xred(:,iatom)
1022 : end do
1023 0 : rcm(:)=rcm(:)/real(natom,dp)-half
1024 0 : do iatom=1,natom
1025 0 : xred(:,iatom)=xred(:,iatom)-rcm(:)
1026 : end do
1027 : ! Also modify the tnons
1028 0 : do isym=1,nsym
1029 0 : tnons(:,isym)=matmul(symrel(:,:,isym),rcm(:))-rcm(:)+tnons(:,isym)
1030 : end do
1031 :
1032 0 : ABI_WARNING('icoulomb is 1 --> the average center of coordinates has been translated to (0.5,0.5,0.5)')
1033 : end if
1034 :
1035 : !========================================================================================================
1036 : !
1037 : ! At this stage, the cell parameters and atomic coordinates are known, as well as the symmetry operations
1038 : ! There has been a preliminary analysis of the holohedry (not definitive, though ...)
1039 : !
1040 : !========================================================================================================
1041 :
1042 : !DEBUG
1043 : !write(std_out,'(a)')' m_ingeo%ingeo : before symanal '
1044 : !call flush(std_out)
1045 : !ENDDEBUG
1046 :
1047 : ! Here, determine correctly the Bravais lattice and other space group or shubnikov group characteristics
1048 7233 : call symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym)
1049 :
1050 : !DEBUG
1051 : !write(std_out,'(a)')' m_ingeo%ingeo : after symanal '
1052 : !call flush(std_out)
1053 : !ENDDEBUG
1054 :
1055 : ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize the rprimd. Keep xred fixed.
1056 7233 : if(tolsym>1.00001e-8)then
1057 : ! Check whether the symmetry operations are consistent with the lattice vectors
1058 7205 : iexit=1
1059 :
1060 7205 : call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tol8)
1061 :
1062 7205 : if(iexit==-1)then
1063 : write(msg,'(5a,es11.3,15a)')&
1064 29 : 'It is observed that the input primitive vectors are not accurate:',ch10,&
1065 29 : 'the lattice is not left invariant within 1.0e-8 when applying symmetry operations.',ch10,&
1066 29 : 'However, they are only slightly inaccurate, as inaccuracies are within the input tolsym=', tolsym,ch10,&
1067 29 : 'In order to avoid spurious effects, the primitive vectors have been',ch10,&
1068 29 : 'symmetrized before storing them in the dataset internal variable.',ch10,&
1069 29 : 'So, do not be surprised by the fact that your input variables (acell, rprim, xcart, xred, ...)',ch10,&
1070 29 : 'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,&
1071 29 : & 'This is not a problem per se.',ch10,&
1072 29 : & 'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',&
1073 29 : & ' decrease tolsym to 1.0e-8 or lower.',ch10,&
1074 58 : 'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.'
1075 29 : ABI_WARNING(msg)
1076 :
1077 29 : call symmetrize_rprimd(bravais,nsym,rprimd,symrel,tol8)
1078 29 : call mkradim(acell,rprim,rprimd)
1079 :
1080 : !Needs one more resymmetrization, for the tnons
1081 87 : ABI_MALLOC(tnons_new,(3,nsym))
1082 : call symmetrize_xred(natom,nsym,symrel,tnons,xred,&
1083 29 : & fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym)
1084 2141 : tnons(:,1:nsym)=tnons_new(:,:)
1085 58 : ABI_FREE(tnons_new)
1086 :
1087 : end if
1088 :
1089 : end if
1090 :
1091 7233 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1092 : angdeg(1)=180.0_dp/pi * acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))
1093 : angdeg(2)=180.0_dp/pi * acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))
1094 : angdeg(3)=180.0_dp/pi * acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))
1095 : !write(std_out,'(a,3f14.8)') ' ingeo: angdeg(1:3)=',angdeg(1:3)
1096 :
1097 : !--------------------------------------------------------------------------------------
1098 :
1099 : !Finally prune the set of symmetry in case non-symmorphic operations must be excluded
1100 7233 : if(symmorphi==0)then
1101 93 : jsym=0
1102 2050 : do isym=1,nsym
1103 7921 : if(sum(tnons(:,isym)**2)<tol6)then
1104 1177 : jsym=jsym+1
1105 : ! This symmetry operation is non-symmorphic, and can be kept
1106 1177 : if(isym/=jsym)then
1107 9568 : symrel(:,:,jsym)=symrel(:,:,isym)
1108 2944 : tnons(:,jsym)=tnons(:,isym)
1109 736 : symafm(jsym)=symafm(isym)
1110 : end if
1111 : end if
1112 : end do
1113 93 : nsym=jsym
1114 : end if
1115 :
1116 : !call symmultsg(nsym,symafm,symrel,tnons)
1117 :
1118 : ! 9) initialize the list of fixed atoms, and initial velocities -----------------
1119 : ! Note: these inputs do not influence the previous generation of
1120 : ! symmetry operations. This might be changed in the future
1121 :
1122 : ! idir=0 is for iatfix , idir=1 is for iatfixx,
1123 : ! idir=2 is for iatfixy, idir=3 is for iatfixz
1124 162237 : iatfix(:,:)=0
1125 :
1126 36165 : do idir=0,3
1127 :
1128 28932 : if(idir==0)then
1129 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfix',tread,'INT')
1130 21699 : else if(idir==1)then
1131 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixx',tread,'INT')
1132 14466 : else if(idir==2)then
1133 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixy',tread,'INT')
1134 7233 : else if(idir==3)then
1135 7233 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixz',tread,'INT')
1136 : end if
1137 :
1138 : ! Use natfix also for natfixx,natfixy,natfixz
1139 28932 : natfix=0
1140 28932 : if(tread==1) natfix=intarr(1)
1141 :
1142 : ! Check the validity of natfix
1143 28932 : if (natfix<0 .or. natfix>natom) then
1144 : write(msg, '(a,a,a,i0,a,i4,a,a,a)' )&
1145 0 : 'The input variables natfix, natfixx, natfixy and natfixz must be',ch10,&
1146 0 : 'between 0 and natom (= ',natom,'), while one of them is ',natfix,'.',ch10,&
1147 0 : 'Action: correct that occurrence in your input file.'
1148 0 : ABI_ERROR(msg)
1149 : end if
1150 :
1151 : !Read iatfix
1152 28932 : if(idir==0)then
1153 7233 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfix',tread,'INT')
1154 21699 : else if(idir==1)then
1155 7233 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixx',tread,'INT')
1156 14466 : else if(idir==2)then
1157 7233 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixy',tread,'INT')
1158 7233 : else if(idir==3)then
1159 7233 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixz',tread,'INT')
1160 : end if
1161 :
1162 : ! If some iatfix was read, natfix must vanish
1163 28932 : if (natfix==0 .and. tread==1)then
1164 : write(msg, '(a,i1,5a)' )&
1165 0 : 'For direction ',idir,' the corresponding natfix is zero,',ch10,&
1166 0 : 'while iatfix specifies some atoms to be fixed.',ch10,&
1167 0 : 'Action: either specify a non-zero natfix(x,y,z) or suppress iatfix(x,y,z).'
1168 0 : ABI_ERROR(msg)
1169 : end if
1170 :
1171 : ! If natfix is non-zero, iatfix must be defined
1172 28932 : if (natfix>0 .and. tread==0)then
1173 : write(msg, '(a,i1,3a,i0,3a)' )&
1174 0 : 'For direction ',idir,' no iatfix has been specified,',ch10,&
1175 0 : 'while natfix specifies that some atoms to be fixed, natfix= ',natfix,'.',ch10,&
1176 0 : 'Action: either set natfix(x,y,z) to zero or define iatfix(x,y,z).'
1177 0 : ABI_ERROR(msg)
1178 : end if
1179 :
1180 36165 : if(tread==1)then
1181 593 : do ii=1,natfix
1182 : ! Checks the validity of the input iatfix
1183 326 : if (intarr(ii)<1 .or. intarr(ii)>natom) then
1184 : write(msg, '(a,a,a,i0,a,a,a)' )&
1185 0 : 'The input variables iatfix, iatfixx, iatfixy and iatfixz must be',ch10,&
1186 0 : 'between 1 and natom, while one of them is ',intarr(ii),'.',ch10,&
1187 0 : 'Action: correct that occurrence in your input file.'
1188 0 : ABI_ERROR(msg)
1189 : end if
1190 : ! Finally set the value of the internal iatfix array
1191 1695 : do iatom=1,natom
1192 1428 : if(intarr(ii)==iatom)then
1193 1217 : if(idir==0)iatfix(1:3,iatom)=1
1194 326 : if(idir/=0)iatfix(idir,iatom)=1
1195 : end if
1196 : end do
1197 : end do
1198 : end if
1199 :
1200 : end do
1201 :
1202 162237 : vel(:,:)=zero
1203 7233 : call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'vel',tread,'DPR')
1204 7239 : if(tread==1)vel(:,:)=reshape( dprarr(1:3*natom), [3, natom])
1205 7233 : call intagm_img(vel,iimage,jdtset,lenstr,nimage,3,natom,string,"vel",tread,'DPR')
1206 :
1207 7233 : vel_cell(:,:)=zero
1208 7233 : call intagm(dprarr,intarr,jdtset,marr,3*3,string(1:lenstr),'vel_cell',tread,'DPR')
1209 7233 : if(tread==1)vel_cell(:,:)=reshape( dprarr(1:9), [3,3])
1210 7233 : call intagm_img(vel_cell,iimage,jdtset,lenstr,nimage,3,3,string,"vel_cell",tread,'DPR')
1211 :
1212 : ! mixalch
1213 7233 : if(ntypalch>0)then
1214 26 : call intagm(dprarr,intarr,jdtset,marr,npspalch*ntypalch,string(1:lenstr),'mixalch',tread,'DPR')
1215 78 : if(tread==1) mixalch(1:npspalch,1:ntypalch)= reshape(dprarr(1:npspalch*ntypalch), [npspalch, ntypalch])
1216 52 : do itypat=1,ntypalch
1217 78 : sumalch=sum(mixalch(1:npspalch,itypat))
1218 52 : if(abs(sumalch-one)>tol10)then
1219 : write(msg, '(a,i0,2a,f8.2,4a)' )&
1220 0 : 'For the alchemical atom number ',itypat,ch10,&
1221 0 : 'the sum of the pseudopotential coefficients is',sumalch,ch10,&
1222 0 : 'while it should be one.',ch10,&
1223 0 : 'Action: check the content of the input variable mixalch.'
1224 0 : ABI_ERROR(msg)
1225 : end if
1226 : end do
1227 26 : call intagm_img(mixalch,iimage,jdtset,lenstr,nimage,npspalch,ntypalch,string,"mixalch",tread,'DPR')
1228 : end if
1229 :
1230 : ! amu (needs mixalch to be initialized ...)
1231 : ! Find the default mass
1232 21699 : ABI_MALLOC(mass_psp,(npsp))
1233 16331 : do ipsp=1,npsp
1234 9098 : call atomdata_from_znucl(atom,znucl(ipsp))
1235 9098 : amu_default = atom%amu
1236 16331 : mass_psp(ipsp)=amu_default
1237 : end do
1238 : ! When the pseudo-atom is pure, simple copy
1239 7233 : ntyppure=ntypat-ntypalch
1240 7233 : if(ntyppure>0)then
1241 16269 : amu(1:ntyppure)=mass_psp(1:ntyppure)
1242 : end if
1243 : ! When the pseudo-atom is alchemical, must make mixing
1244 7233 : if(ntypalch>0)then
1245 52 : do itypat=ntyppure+1,ntypat
1246 26 : amu(itypat)=zero
1247 104 : do ipsp=ntyppure+1,npsp
1248 78 : amu(itypat)=amu(itypat)+mixalch(ipsp-ntyppure,itypat-ntyppure)*mass_psp(ipsp)
1249 : end do
1250 : end do
1251 : end if
1252 7233 : ABI_FREE(mass_psp)
1253 :
1254 7233 : call intagm(dprarr,intarr,jdtset,marr,ntypat,string(1:lenstr),'amu',tread,'DPR')
1255 7765 : if(tread==1)amu(:)=dprarr(1:ntypat)
1256 7233 : call intagm_img(amu,iimage,jdtset,lenstr,nimage,ntypat,string,"amu",tread,'DPR')
1257 :
1258 7233 : ABI_FREE(intarr)
1259 7233 : ABI_FREE(dprarr)
1260 :
1261 : !DEBUG
1262 : !write(std_out,'(a)')' m_ingeo%ingeo : exit '
1263 : !call flush(std_out)
1264 : !ENDDEBUG
1265 :
1266 14466 : end subroutine ingeo
1267 : !!***
1268 :
1269 : !!****f* m_ingeo/ingeobld
1270 : !! NAME
1271 : !! ingeobld
1272 : !!
1273 : !! FUNCTION
1274 : !! The geometry builder.
1275 : !! Start from the types and coordinates of the primitive atoms
1276 : !! and produce the completed set of atoms, by using the definition
1277 : !! of objects, then application of rotation, translation and repetition.
1278 : !!
1279 : !! INPUTS
1280 : !! iout=unit number of output file
1281 : !! jdtset=number of the dataset looked for
1282 : !! lenstr=actual length of the string
1283 : !! natrd=number of atoms that have been read in the calling routine
1284 : !! natom=number of atoms
1285 : !! nobj=the number of objects
1286 : !! string*(*)=character string containing all the input data. Initialized previously in instrng.
1287 : !! typat_read(natrd)=type integer for each atom in the primitive set
1288 : !! xcart_read(3,natrd)=cartesian coordinates of atoms (bohr), in the primitive set
1289 : !!
1290 : !! OUTPUT
1291 : !! typat(natom)=type integer for each atom in cell
1292 : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
1293 : !!
1294 : !! SOURCE
1295 :
1296 4 : subroutine ingeobld (iout,jdtset,lenstr,natrd,natom,nobj,string,typat,typat_read,xcart,xcart_read)
1297 :
1298 : !Arguments ------------------------------------
1299 : !scalars
1300 : integer,intent(in) :: iout,jdtset,lenstr,natom,natrd,nobj
1301 : character(len=*),intent(in) :: string
1302 : !arrays
1303 : integer,intent(in) :: typat_read(natrd)
1304 : integer,intent(out) :: typat(natom)
1305 : real(dp),intent(in) :: xcart_read(3,natrd)
1306 : real(dp),intent(out) :: xcart(3,natom)
1307 :
1308 : !Local variables-------------------------------
1309 : character(len=*), parameter :: format01110 ="(1x,a6,1x,(t9,8i8) )"
1310 : character(len=*), parameter :: format01160 ="(1x,a6,1x,1p,(t9,3g18.10)) "
1311 : !scalars
1312 : integer :: belonga,belongb,iatom,iatrd,ii,irep,irep1,irep2,irep3,ivac,marr
1313 : integer :: natom_toberead,nread,objan,objbn,rotate,shift,tread,vacnum
1314 : real(dp) :: angle,cosine,norm2per,norma,normb,normper,project,sine
1315 : character(len=500) :: msg
1316 : !arrays
1317 : integer :: objarf(3),objbrf(3)
1318 4 : integer,allocatable :: objaat(:),objbat(:),typat_full(:),vaclst(:)
1319 : real(dp) :: axis2(3),axis3(3),axisa(3),axisb(3),objaax(6),objaro(4),objatr(12)
1320 : real(dp) :: objbax(6),objbro(4),objbtr(12),parall(3),perpen(3),rotated(3)
1321 : real(dp) :: vectora(3),vectorb(3)
1322 4 : real(dp),allocatable :: xcart_full(:,:)
1323 4 : integer,allocatable :: intarr(:)
1324 4 : real(dp),allocatable :: dprarr(:)
1325 :
1326 : ! *************************************************************************
1327 :
1328 4 : marr=max(12,3*natom)
1329 12 : ABI_MALLOC(intarr,(marr))
1330 12 : ABI_MALLOC(dprarr,(marr))
1331 :
1332 : !1) Set up the number of vacancies.
1333 :
1334 : !This is the default
1335 4 : vacnum=0
1336 4 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacnum',tread,'INT')
1337 4 : if(tread==1) vacnum=intarr(1)
1338 :
1339 4 : if (vacnum>0)then
1340 3 : ABI_MALLOC(vaclst,(vacnum))
1341 : ! Read list of atoms to be suppressed to create vacancies
1342 1 : call intagm(dprarr,intarr,jdtset,marr,vacnum,string(1:lenstr),'vaclst',tread,'INT')
1343 2 : if(tread==1) vaclst(:)=intarr(1:vacnum)
1344 1 : if(tread/=1)then
1345 : write(msg, '(a,a,a,a,a)' )&
1346 0 : & 'The array vaclst MUST be initialized in the input file',ch10,&
1347 0 : & 'when vacnum is non-zero.',ch10,&
1348 0 : & 'Action: initialize vaclst in your input file.'
1349 0 : ABI_ERROR(msg)
1350 : end if
1351 : end if
1352 :
1353 4 : natom_toberead=natom+vacnum
1354 :
1355 : !2) Set up list and number of atoms in objects, and the --------------
1356 : !operations to be performed on objects.
1357 :
1358 324 : write(msg,'(80a,a)')('=',ii=1,80),ch10
1359 4 : call wrtout(std_out,msg)
1360 4 : call wrtout(iout,msg)
1361 :
1362 4 : write(msg, '(a,a)' )'--ingeobld: echo values of variables connected to objects --------',ch10
1363 4 : call wrtout(std_out,msg)
1364 4 : call wrtout(iout,msg)
1365 :
1366 4 : if(vacnum>0)then
1367 1 : write(iout,format01110) 'vacnum',vacnum
1368 1 : write(std_out,format01110) 'vacnum',vacnum
1369 1 : write(iout,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:)
1370 1 : write(std_out,'(1x,a6,1x,(t9,20i3))') 'vaclst',vaclst(:)
1371 1 : write(iout, '(a)' ) ' '
1372 1 : write(std_out,'(a)' ) ' '
1373 : end if
1374 :
1375 4 : write(iout,format01110) 'nobj',nobj
1376 4 : write(std_out,format01110) 'nobj',nobj
1377 :
1378 4 : if(nobj/=1 .and. nobj/=2)then
1379 : write(msg, '(a,a,a,i8,a,a,a)' )&
1380 0 : & 'The number of object (nobj) must be either 1 or 2,',ch10,&
1381 0 : & 'while the input file has nobj=',nobj,'.',ch10,&
1382 0 : & 'Action: correct nobj in your input file.'
1383 0 : ABI_ERROR(msg)
1384 : end if
1385 :
1386 4 : if(nobj==1 .or. nobj==2)then
1387 :
1388 : ! Read the number of atoms of the object a
1389 4 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objan',tread,'INT')
1390 4 : if(tread==1) objan=intarr(1)
1391 :
1392 4 : if(tread/=1)then
1393 : write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1394 0 : & 'The number of atoms in object a (objan) must be initialized',ch10,&
1395 0 : & 'in the input file, when nobj=',nobj,'.',ch10,&
1396 0 : & 'This is not the case.',ch10,&
1397 0 : & 'Action: correct objan in your input file.'
1398 0 : ABI_ERROR(msg)
1399 : end if
1400 :
1401 4 : write(iout, '(a)' ) ' '
1402 4 : write(std_out,'(a)' ) ' '
1403 4 : write(iout,format01110) 'objan',objan
1404 4 : write(std_out,format01110) 'objan',objan
1405 :
1406 4 : if(objan<=1 .or. objan>natom)then
1407 : write(msg, '(a,a,a,a,a,i8,a,a,a)' )&
1408 0 : & 'The number of atoms in object a (objan) must be larger than 0',ch10,&
1409 0 : & 'and smaller than natom.',ch10,&
1410 0 : & 'It is equal to ',objan,', an unacceptable value.',ch10,&
1411 0 : & 'Action: correct objan in your input file.'
1412 0 : ABI_ERROR(msg)
1413 : end if
1414 :
1415 : ! Read list of atoms in object a
1416 4 : call intagm(dprarr,intarr,jdtset,marr,objan,string(1:lenstr),'objaat',tread,'INT')
1417 12 : ABI_MALLOC(objaat,(objan))
1418 26 : if(tread==1) objaat(1:objan)=intarr(1:objan)
1419 :
1420 4 : if(tread/=1)then
1421 : write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1422 0 : & 'The list of atoms in object a (objaat) must be initialized',ch10,&
1423 0 : & 'in the input file, when nobj=',nobj,'.',ch10,&
1424 0 : & 'This is not the case.',ch10,&
1425 0 : & 'Action: initialize objaat in your input file.'
1426 0 : ABI_ERROR(msg)
1427 : end if
1428 :
1429 4 : write(iout,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:)
1430 4 : write(std_out,'(1x,a6,1x,(t9,20i3))') 'objaat',objaat(:)
1431 :
1432 26 : do iatom=1,objan
1433 26 : if(objaat(iatom)<1 .or. objaat(iatom)>natom)then
1434 : write(msg, '(a,i8,a,a,i8,4a)' )&
1435 0 : & 'The input value of objaat for atom number ',iatom,ch10,&
1436 0 : & 'is equal to ',objaat(iatom),', an unacceptable value :',ch10,&
1437 0 : & 'it should be between 1 and natom. ',&
1438 0 : & 'Action: correct the array objaat in your input file.'
1439 0 : ABI_ERROR(msg)
1440 : end if
1441 : end do
1442 :
1443 4 : if(objan>1)then
1444 22 : do iatom=1,objan-1
1445 22 : if( objaat(iatom)>=objaat(iatom+1) )then
1446 : write(msg, '(a,i8,a,a,a,a,a,a)' )&
1447 0 : & 'The input value of objaat for atom number ',iatom,ch10,&
1448 0 : & 'is larger or equal to the one of the next atom,',ch10,&
1449 0 : & 'while this list should be ordered, and an atom cannot be repeated.',ch10,&
1450 0 : & 'Action: correct the array objaat in your input file.'
1451 0 : ABI_ERROR(msg)
1452 : end if
1453 : end do
1454 : end if
1455 :
1456 : ! Read repetition factors
1457 16 : objarf(1:3)=1
1458 4 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objarf',tread,'INT')
1459 16 : if(tread==1) objarf(1:3)=intarr(1:3)
1460 4 : write(iout,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:)
1461 4 : write(std_out,'(1x,a6,1x,(t9,20i3))') 'objarf',objarf(:)
1462 :
1463 4 : if(tread==1)then
1464 16 : do irep=1,3
1465 16 : if(objarf(irep)<1)then
1466 : write(msg, '(a,a,a,3i8,a,a,a)' )&
1467 0 : & 'The input values of objarf(1:3) must be positive,',ch10,&
1468 0 : & 'while it is ',objarf(1:3),'.',ch10,&
1469 0 : & 'Action: correct objarf in your input file.'
1470 0 : ABI_ERROR(msg)
1471 : end if
1472 : end do
1473 : end if
1474 :
1475 : ! Modify the number of atoms to be read
1476 4 : natom_toberead=natom_toberead-objan*(objarf(1)*objarf(2)*objarf(3)-1)
1477 :
1478 : ! Read rotations angles and translations
1479 4 : objaro(1:4)=0.0_dp
1480 4 : objatr(1:12)=0.0_dp
1481 4 : if (objarf(1)*objarf(2)*objarf(3) ==1) then
1482 0 : nread=1
1483 4 : else if (objarf(2)*objarf(3) ==1) then
1484 2 : nread=2
1485 2 : else if (objarf(3) ==1) then
1486 0 : nread=3
1487 : else
1488 2 : nread=4
1489 : end if
1490 4 : call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objaro',tread,'DPR')
1491 8 : if(tread==1) objaro(1:nread)=dprarr(1:nread)
1492 :
1493 4 : call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objatr',tread,'LEN')
1494 :
1495 40 : if(tread==1) objatr(1:3*nread)=dprarr(1:3*nread)
1496 4 : write(iout,format01160) 'objaro',objaro(1:4)
1497 4 : write(std_out,format01160) 'objaro',objaro(1:4)
1498 4 : write(iout,format01160) 'objatr',objatr(1:12)
1499 4 : write(std_out,format01160) 'objatr',objatr(1:12)
1500 : ! If needed, read axes, but default to the x-axis to avoid errors later
1501 4 : objaax(1:6)=0.0_dp ; objaax(4)=1.0_dp
1502 :
1503 4 : if(abs(objaro(1))+abs(objaro(2))+abs(objaro(3))+abs(objaro(4)) > 1.0d-10) then
1504 2 : call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objaax',tread,'LEN')
1505 14 : if(tread==1) objaax(1:6)=dprarr(1:6)
1506 2 : if(tread/=1)then
1507 : write(msg, '(a,a,a,a,a,a,a)' )&
1508 0 : & 'The axis of object a (objaax) must be initialized',ch10,&
1509 0 : & 'in the input file, when rotations (objaro) are present.',ch10,&
1510 0 : & 'This is not the case.',ch10,&
1511 0 : & 'Action: initialize objaax in your input file.'
1512 0 : ABI_ERROR(msg)
1513 : end if
1514 2 : write(iout,format01160) 'objaax',objaax(1:6)
1515 2 : write(std_out,format01160) 'objaax',objaax(1:6)
1516 : end if
1517 :
1518 16 : axisa(1:3)=objaax(4:6)-objaax(1:3)
1519 4 : norma=axisa(1)**2+axisa(2)**2+axisa(3)**2
1520 :
1521 4 : if(norma<1.0d-10)then
1522 : write(msg, '(5a)' )&
1523 0 : & 'The two points defined by the input array objaax are too',ch10,&
1524 0 : & 'close to each other, and will not be used to define an axis.',ch10,&
1525 0 : & 'Action: correct objaax in your input file.'
1526 0 : ABI_ERROR(msg)
1527 : end if
1528 16 : axisa(1:3)=axisa(1:3)/sqrt(norma)
1529 : end if ! End condition of existence of a first object
1530 :
1531 4 : if(nobj==2)then
1532 :
1533 : ! Read the number of atoms of the object b
1534 3 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'objbn',tread,'INT')
1535 3 : if(tread==1) objbn=intarr(1)
1536 :
1537 3 : if(tread/=1)then
1538 : write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1539 0 : & 'The number of atoms in object b (objbn) must be initialized',ch10,&
1540 0 : & 'in the input file, when nobj=',nobj,'.',ch10,&
1541 0 : & 'This is not the case.',ch10,&
1542 0 : & 'Action: initialize objbn in your input file.'
1543 0 : ABI_ERROR(msg)
1544 : end if
1545 :
1546 3 : write(iout, '(a)' ) ' '
1547 3 : write(std_out,'(a)' ) ' '
1548 3 : write(iout,format01110) 'objbn',objbn
1549 3 : write(std_out,format01110) 'objbn',objbn
1550 :
1551 3 : if(objbn<=1 .or. objbn>natom)then
1552 : write(msg, '(a,a,a,a,a,i8,a,a,a)' )&
1553 0 : & 'The number of atoms in object b (objbn) must be larger than 0',ch10,&
1554 0 : & 'and smaller than natom.',ch10,&
1555 0 : & 'It is equal to ',objbn,', an unacceptable value.',ch10,&
1556 0 : & 'Action: correct objbn in your input file.'
1557 0 : ABI_ERROR(msg)
1558 : end if
1559 :
1560 : ! Read list of atoms in object b
1561 3 : call intagm(dprarr,intarr,jdtset,marr,objbn,string(1:lenstr),'objbat',tread,'INT')
1562 9 : ABI_MALLOC(objbat,(objbn))
1563 :
1564 22 : if(tread==1) objbat(1:objbn)=intarr(1:objbn)
1565 3 : if(tread/=1)then
1566 : write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
1567 0 : & 'The list of atoms in object b (objbat) must be initialized',ch10,&
1568 0 : & 'in the input file, when nobj=',nobj,'.',ch10,&
1569 0 : & 'This is not the case.',ch10,&
1570 0 : & 'Action: initialize objbat in your input file.'
1571 0 : ABI_ERROR(msg)
1572 : end if
1573 :
1574 3 : write(iout,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:)
1575 3 : write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbat',objbat(:)
1576 :
1577 22 : do iatom=1,objbn
1578 22 : if(objbat(iatom)<1 .or. objbat(iatom)>natom)then
1579 : write(msg, '(a,i8,a,a,i8,a,a,a,a,a)' )&
1580 0 : & 'The input value of objbat for atom number ',iatom,ch10,&
1581 0 : & 'is equal to ',objbat(iatom),', an unacceptable value :',ch10,&
1582 0 : & 'it should be between 1 and natom. ',ch10,&
1583 0 : & 'Action: correct objbat in your input file.'
1584 0 : ABI_ERROR(msg)
1585 : end if
1586 : end do
1587 :
1588 3 : if(objbn>1)then
1589 19 : do iatom=1,objbn-1
1590 19 : if( objbat(iatom)>=objbat(iatom+1) )then
1591 : write(msg, '(a,i8,a,a,a,a,a,a)' )&
1592 0 : & 'The input value of objbat for atom number ',iatom,ch10,&
1593 0 : & 'is larger or equal to the one of the next atom,',ch10,&
1594 0 : & 'while this list should be ordered, and an atom cannot be repeated.',ch10,&
1595 0 : & 'Action: correct the array objbat in the input file.'
1596 0 : ABI_ERROR(msg)
1597 : end if
1598 : end do
1599 : end if
1600 :
1601 : ! Read repetition factors
1602 12 : objbrf(1:3)=1
1603 3 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'objbrf',tread,'INT')
1604 12 : if(tread==1) objbrf(1:3)=intarr(1:3)
1605 3 : write(iout,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:)
1606 3 : write(std_out,'(1x,a6,1x,(t9,20i3))') 'objbrf',objbrf(:)
1607 :
1608 3 : if(tread==1)then
1609 12 : do irep=1,3
1610 12 : if(objbrf(irep)<1)then
1611 : write(msg, '(a,a,a,3i8,a,a,a)' )&
1612 0 : & 'The input values of objbrf(1:3) must be positive,',ch10,&
1613 0 : & 'while it is ',objbrf(1:3),'.',ch10,&
1614 0 : & 'Action: correct objbrf in your input file.'
1615 0 : ABI_ERROR(msg)
1616 : end if
1617 : end do
1618 : end if
1619 :
1620 : ! Modify the number of atoms to be read
1621 3 : natom_toberead=natom_toberead-objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1)
1622 : ! Read rotations angles and translations
1623 3 : objbro(1:4)=0.0_dp
1624 3 : objbtr(1:12)=0.0_dp
1625 3 : if (objbrf(1)*objbrf(2)*objbrf(3) ==1) then
1626 0 : nread=1
1627 3 : else if (objbrf(2)*objbrf(3) ==1) then
1628 1 : nread=2
1629 2 : else if (objbrf(3) ==1) then
1630 0 : nread=3
1631 : else
1632 2 : nread=4
1633 : end if
1634 3 : call intagm(dprarr,intarr,jdtset,marr,nread,string(1:lenstr),'objbro',tread,'DPR')
1635 5 : if(tread==1) objbro(1:nread)=dprarr(1:nread)
1636 :
1637 3 : call intagm(dprarr,intarr,jdtset,marr,3*nread,string(1:lenstr),'objbtr',tread,'LEN')
1638 33 : if(tread==1) objbtr(1:3*nread)=dprarr(1:3*nread)
1639 :
1640 3 : write(iout,format01160) 'objbro',objbro(1:4)
1641 3 : write(std_out,format01160) 'objbro',objbro(1:4)
1642 3 : write(iout,format01160) 'objbtr',objbtr(1:12)
1643 3 : write(std_out,format01160) 'objbtr',objbtr(1:12)
1644 :
1645 : ! If needed, read axes, but default to the x-axis to avoid errors later
1646 3 : objbax(1:6)=0.0_dp ; objbax(4)=1.0_dp
1647 3 : if(abs(objbro(1))+abs(objbro(2))+abs(objbro(3))+abs(objbro(4)) > 1.0d-10) then
1648 1 : call intagm(dprarr,intarr,jdtset,marr,6,string(1:lenstr),'objbax',tread,'LEN')
1649 7 : if(tread==1) objbax(1:6)=dprarr(1:6)
1650 1 : if(tread/=1)then
1651 : write(msg, '(a,a,a,a,a,a,a)' )&
1652 0 : & 'The axis of object b (objbax) must be initialized',ch10,&
1653 0 : & 'in the input file, when rotations (objbro) are present.',ch10,&
1654 0 : & 'This is not the case.',ch10,&
1655 0 : & 'Action: initialize objbax in your input file.'
1656 0 : ABI_ERROR(msg)
1657 : end if
1658 1 : write(iout,format01160) 'objbax',objbax(1:6)
1659 1 : write(std_out,format01160) 'objbax',objbax(1:6)
1660 : end if
1661 12 : axisb(1:3)=objbax(4:6)-objbax(1:3)
1662 3 : normb=axisb(1)**2+axisb(2)**2+axisb(3)**2
1663 3 : if(normb<1.0d-10)then
1664 : write(msg, '(5a)' )&
1665 0 : & 'The two points defined by the input array objbax are too',ch10,&
1666 0 : & 'close to each other, and will not be used to define an axis.',ch10,&
1667 0 : & 'Action: correct objbax in your input file.'
1668 0 : ABI_ERROR(msg)
1669 : end if
1670 12 : axisb(1:3)=axisb(1:3)/sqrt(normb)
1671 :
1672 : ! Check whether both lists are disjoints. Use a very primitive algorithm.
1673 17 : do iatom=1,objan
1674 89 : do ii=1,objbn
1675 86 : if(objaat(iatom)==objbat(ii))then
1676 : write(msg, '(6a,i8,a,i8,3a)' )&
1677 0 : & 'The objects a and b cannot have a common atom, but it is',ch10,&
1678 0 : & 'found that the values of objaat and objbat ',&
1679 0 : & ' are identical, for their',ch10,&
1680 0 : & 'atoms number ',iatom,' and ',ii,'.',ch10,&
1681 0 : & 'Action: change objaat and/or objbat so that they have no common atom anymore.'
1682 0 : ABI_ERROR(msg)
1683 : end if
1684 : end do
1685 : end do
1686 : end if ! End condition of existence of a second object
1687 :
1688 : !Check whether the number of atoms to be read obtained by relying
1689 : !on natom, vacnum and the object definitions, or from natrd coincide
1690 4 : if(natrd/=natom_toberead)then
1691 : write(msg,'(11a,i0,a,i0,2a,i0,a)' )&
1692 0 : & ' ingeobld : ERROR -- ',ch10,&
1693 0 : & ' The number of atoms to be read (natrd) must be equal',ch10,&
1694 0 : & ' to the total number of atoms (natom), plus',ch10,&
1695 0 : & ' the number of vacancies (vacnum), minus',ch10,&
1696 0 : & ' the number of atoms added by the repetition of objects.',ch10,&
1697 0 : & ' This is not the case : natrd= ',natrd,', natom= ',natom,ch10,&
1698 0 : & ', vacnum= ',vacnum,';'
1699 0 : call wrtout(std_out,msg)
1700 :
1701 0 : if(nobj==1 .or. nobj==2) then
1702 : write(msg,'(a,i3,a,3i3,a,i5,a)' )&
1703 0 : & ' object a : objan=',objan,', objarf(1:3)=',objarf(1:3),&
1704 0 : & ' => adds ',objan*(objarf(1)*objarf(2)*objarf(3)-1),' atoms.'
1705 0 : call wrtout(std_out,msg)
1706 : end if
1707 :
1708 0 : if(nobj==2) then
1709 : write(msg,'(a,i3,a,3i3,a,i5,a)' )&
1710 0 : & ' object b : objbn=',objbn,', objbrf(1:3)=',objbrf(1:3),&
1711 0 : & ' => adds ',objbn*(objbrf(1)*objbrf(2)*objbrf(3)-1),' atoms.'
1712 0 : call wrtout(std_out,msg)
1713 : end if
1714 :
1715 : write(msg,'(3a)' )&
1716 0 : & ' Action : check the correspondence between natom+vacnum on one side,',ch10,&
1717 0 : & ' and natrd, objan, objbn, objarf and objbrf on the other side.'
1718 0 : ABI_ERROR(msg)
1719 : end if
1720 :
1721 : !6) Produce full set of atoms
1722 :
1723 : !Print the initial atom coordinates if the geometry builder is used
1724 4 : write(iout, '(/,a)' ) ' Cartesian coordinates of the primitive atoms '
1725 4 : write(std_out,'(/,a)' )' Cartesian coordinates of the primitive atoms '
1726 4 : write(iout,format01160) ' ',xcart_read(:,:)
1727 4 : write(std_out,format01160) ' ',xcart_read(:,:)
1728 :
1729 12 : ABI_MALLOC(typat_full,(natom+vacnum))
1730 12 : ABI_MALLOC(xcart_full,(3,natom+vacnum))
1731 :
1732 : !Use the work array xcart_full to produce full set of atoms,
1733 : !including those coming from repeated objects.
1734 4 : iatom=1
1735 93 : do iatrd=1,natrd
1736 :
1737 89 : belonga=0 ; belongb=0
1738 89 : if(nobj==1 .or. nobj==2)then
1739 : ! Determine whether the atom belongs to object a
1740 531 : do ii=1,objan
1741 531 : if(iatrd==objaat(ii))belonga=ii
1742 : end do
1743 : end if
1744 89 : if(nobj==2)then
1745 : ! Determine whether the atom belong to object b
1746 594 : do ii=1,objbn
1747 594 : if(iatrd==objbat(ii))belongb=ii
1748 : end do
1749 : end if
1750 :
1751 : !write(std_out,'(a,i5,a,i2,i2,a)' )' ingeobld : treating iatrd=',iatrd,', belong(a,b)=',belonga,belongb,'.'
1752 :
1753 : ! In case it does not belong to an object
1754 93 : if(belonga==0 .and. belongb==0)then
1755 192 : xcart_full(1:3,iatom)=xcart_read(1:3,iatrd)
1756 48 : typat_full(iatom)=typat_read(iatrd)
1757 48 : iatom=iatom+1
1758 : else
1759 :
1760 : ! Repeat, rotate and translate this atom
1761 41 : if(belonga/=0)then
1762 :
1763 : ! Treat object a
1764 : ! Compute the relative coordinate of atom with respect to first point of axis
1765 88 : vectora(1:3)=xcart_read(1:3,iatrd)-objaax(1:3)
1766 : ! Project on axis
1767 22 : project=vectora(1)*axisa(1)+vectora(2)*axisa(2)+vectora(3)*axisa(3)
1768 : ! Get the parallel part
1769 88 : parall(1:3)=project*axisa(1:3)
1770 : ! Get the perpendicular part, to be rotated
1771 88 : perpen(1:3)=vectora(1:3)-parall(1:3)
1772 : ! Compute the norm of the perpendicular part
1773 22 : norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2
1774 : ! Initialisation to avoid warnings even if used behind if rotate == 1.
1775 22 : normper = 0
1776 : ! It the norm is too small, there is not need to rotate
1777 22 : rotate=0
1778 22 : if(norm2per>=1.0d-18)then
1779 18 : rotate=1
1780 18 : normper=sqrt(norm2per)
1781 72 : axis2(1:3)=perpen(1:3)/normper
1782 : ! Get the vector perpendicular to axisa and axisa2
1783 18 : axis3(1)=axisa(2)*axis2(3)-axisa(3)*axis2(2)
1784 18 : axis3(2)=axisa(3)*axis2(1)-axisa(1)*axis2(3)
1785 18 : axis3(3)=axisa(1)*axis2(2)-axisa(2)*axis2(1)
1786 : end if
1787 :
1788 : ! Here the repetition loop
1789 60 : do irep3=1,objarf(3)
1790 130 : do irep2=1,objarf(2)
1791 260 : do irep1=1,objarf(1)
1792 : ! Here the rotation
1793 152 : if(rotate==1)then
1794 : ! Compute the angle of rotation
1795 : angle=objaro(1)+(irep1-1)*objaro(2) + &
1796 128 : & (irep2-1)*objaro(3)+(irep3-1)*objaro(4)
1797 128 : cosine=cos(angle/180.0*pi)
1798 128 : sine=sin(angle/180.0*pi)
1799 : rotated(1:3)=objaax(1:3)+parall(1:3)+&
1800 512 : & normper*(cosine*axis2(1:3)+sine*axis3(1:3))
1801 : else
1802 24 : rotated(1:3)=vectora(1:3)
1803 : end if
1804 : ! Here the translation
1805 : xcart_full(1:3,iatom)=rotated(1:3)+objatr(1:3)+&
1806 608 : & (irep1-1)*objatr(4:6)+(irep2-1)*objatr(7:9)+(irep3-1)*objatr(10:12)
1807 152 : typat_full(iatom)=typat_read(iatrd)
1808 222 : iatom=iatom+1
1809 : end do
1810 : end do
1811 : end do ! End the repetition loop
1812 :
1813 : else
1814 : ! If the atom belong to object b
1815 : ! Compute the relative coordinate of atom with respect to first point of axis
1816 76 : vectorb(1:3)=xcart_read(1:3,iatrd)-objbax(1:3)
1817 : ! Project on axis
1818 19 : project=vectorb(1)*axisb(1)+vectorb(2)*axisb(2)+vectorb(3)*axisb(3)
1819 : ! Get the parallel part
1820 76 : parall(1:3)=project*axisb(1:3)
1821 : ! Get the perpendicular part, to be rotated
1822 76 : perpen(1:3)=vectorb(1:3)-parall(1:3)
1823 : ! Compute the norm of the perpendicular part
1824 19 : norm2per=perpen(1)**2+perpen(2)**2+perpen(3)**2
1825 : ! Initialisation to avoid warnings even if used behind if rotate == 1.
1826 19 : normper = 0
1827 : ! It the norm is too small, there is not need to rotate
1828 19 : rotate=0
1829 19 : if(norm2per>=1.0d-18)then
1830 16 : rotate=1
1831 16 : normper=sqrt(norm2per)
1832 64 : axis2(1:3)=perpen(1:3)/normper
1833 : ! Get the vector perpendicular to axisb and axis2
1834 16 : axis3(1)=axisb(2)*axis2(3)-axisb(3)*axis2(2)
1835 16 : axis3(2)=axisb(3)*axis2(1)-axisb(1)*axis2(3)
1836 16 : axis3(3)=axisb(1)*axis2(2)-axisb(2)*axis2(1)
1837 : end if
1838 : ! Here the repetition loop
1839 54 : do irep3=1,objbrf(3)
1840 121 : do irep2=1,objbrf(2)
1841 242 : do irep1=1,objbrf(1)
1842 : ! Here the rotation
1843 140 : if(rotate==1)then
1844 : ! Compute the angle of rotation
1845 : angle=objbro(1)+(irep1-1)*objbro(2) + &
1846 120 : & (irep2-1)*objbro(3)+ (irep3-1)*objbro(4)
1847 120 : cosine=cos(angle/180.0*pi)
1848 120 : sine=sin(angle/180.0*pi)
1849 : rotated(1:3)=objbax(1:3)+parall(1:3)+&
1850 480 : & normper*(cosine*axis2(1:3)+sine*axis3(1:3))
1851 : else
1852 20 : rotated(1:3)=vectorb(1:3)
1853 : end if
1854 : ! Here the translation
1855 : xcart_full(1:3,iatom)=rotated(1:3)+objbtr(1:3)+&
1856 560 : & (irep1-1)*objbtr(4:6)+(irep2-1)*objbtr(7:9)+(irep3-1)*objbtr(10:12)
1857 140 : typat_full(iatom)=typat_read(iatrd)
1858 207 : iatom=iatom+1
1859 : end do
1860 : end do
1861 : end do ! End the repetition loop
1862 : end if ! Condition of belonging to object b
1863 : end if ! Condition of belonging to an object
1864 : end do ! Loop on atoms
1865 :
1866 : !Create the vacancies here
1867 4 : if(vacnum/=0)then
1868 : ! First label the vacant atoms as belonging to typat 0
1869 2 : do ivac=1,vacnum
1870 2 : typat_full(vaclst(ivac))=0
1871 : end do
1872 : ! Then compact the arrays
1873 1 : shift=0
1874 64 : do iatom=1,natom
1875 63 : if(typat_full(iatom+shift)==0) shift=shift+1
1876 64 : if(shift/=0)then
1877 252 : xcart_full(1:3,iatom)=xcart_full(1:3,iatom+shift)
1878 63 : typat_full(iatom)=typat_full(iatom+shift)
1879 : end if
1880 : end do
1881 : end if
1882 :
1883 : !Transfer the content of xcart_full and typat_full to the proper location
1884 1360 : xcart(:,1:natom)=xcart_full(:,1:natom)
1885 343 : typat(1:natom)=typat_full(1:natom)
1886 :
1887 4 : ABI_FREE(typat_full)
1888 4 : ABI_FREE(xcart_full)
1889 4 : if(allocated(objaat)) then
1890 4 : ABI_FREE(objaat)
1891 : end if
1892 4 : if(allocated(objbat)) then
1893 3 : ABI_FREE(objbat)
1894 : end if
1895 :
1896 4 : ABI_FREE(intarr)
1897 4 : ABI_FREE(dprarr)
1898 4 : if (vacnum>0) then
1899 1 : ABI_FREE(vaclst)
1900 : end if
1901 :
1902 4 : end subroutine ingeobld
1903 : !!***
1904 :
1905 : !!****f* m_ingeo/fillcell
1906 : !! NAME
1907 : !! fillcell
1908 : !!
1909 : !! FUNCTION
1910 : !! Computes the atomic position of all the atoms in the unit cell starting
1911 : !! with the symmetry operations and the atoms from the asymmetric unit cell.
1912 : !!
1913 : !! INPUTS
1914 : !! chrgat(natom)=target charge for each atom. Not always used, it depends on the value of constraint_kind
1915 : !! natrd = number of atoms in the asymmetric unit cell
1916 : !! natom = total number of atoms (to be checked)
1917 : !! nsym = number of symmetry operations
1918 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
1919 : !! symrel(3,3,nsym)=symmetry operations in real space in terms
1920 : !! of primitive translations
1921 : !! tnons(3,nsym)=nonsymmorphic translations for symmetry operations
1922 : !! tolsym=tolerance on symmetries
1923 : !! typat(1:natrd)=type integer for each atom in cell
1924 : !! xred(3,1:natrd)=reduced dimensionless atomic coordinates
1925 : !!
1926 : !! OUTPUT
1927 : !!
1928 : !! SIDE EFFECTS
1929 : !! At input, for the asymmetric unit cell
1930 : !! nucdipmom(3,1:natrd)=nuclear magnetic dipole moments of the atoms
1931 : !! spinat(3,1:natrd)=spin-magnetization of the atoms
1932 : !! spinat_cart(3,1:natrd)=spin-magnetization of the atoms (Cartesian)
1933 : !! typat(1:natrd)=type integer for each atom in cell
1934 : !! xred(3,1:natrd)=reduced dimensionless atomic coordinates
1935 : !!
1936 : !! At output, for the complete unit cell
1937 : !! nucdipmom(3,1:natom)=nuclear magnetic dipole moments of the atoms
1938 : !! spinat(3,1:natom)=spin-magnetization of the atoms
1939 : !! spinat_cart(3,1:natom)=spin-magnetization of the atoms (Cartesian)
1940 : !! typat(1:natom)=type integer for each atom in cell
1941 : !! xred(3,1:natom)=reduced dimensionless atomic coordinates
1942 : !!
1943 : !! SOURCE
1944 :
1945 1588 : subroutine fillcell(chrgat,natom,natrd,nsym,nucdipmom,spinat,spinat_cart,symafm,symrel,tnons,tolsym,typat,xred)
1946 :
1947 : !Arguments ------------------------------------
1948 : !scalars
1949 : integer,intent(in) :: natom,natrd,nsym
1950 : !arrays
1951 : integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
1952 : integer,intent(inout) :: typat(natom)
1953 : real(dp),intent(in) :: tolsym
1954 : real(dp),intent(in) :: tnons(3,nsym)
1955 : real(dp),intent(inout) :: chrgat(natom),nucdipmom(3,natom),spinat(3,natom),spinat_cart(3,natom),xred(3,natom)
1956 :
1957 : !Local variables ------------------------------
1958 : !scalars
1959 : integer :: curat,flagch,flageq,ii,iij,jj,kk
1960 : character(len=500) :: msg
1961 : !arrays
1962 3176 : integer :: bcktypat(nsym*natrd)
1963 3176 : real(dp) :: bckat(3),bcknucdipmom(3,nsym*natrd)
1964 3176 : real(dp) :: bckchrgat(nsym*natrd),bckspinat(3,nsym*natrd),bckspinat_cart(3,nsym*natrd),bckxred(3,nsym*natrd)
1965 :
1966 : ! *************************************************************************
1967 :
1968 : !DEBUG
1969 : !write(std_out,*)' fillcell : enter with nsym, natrd= ',nsym,natrd
1970 : !write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
1971 : !do ii=1,nsym
1972 : !write(std_out,'(i3,2x,9i3,3es12.2,i3)')ii,symrel(:,:,ii),tnons(:,ii),symafm(ii)
1973 : !end do
1974 : !write(std_out,*)' Describe the input atoms (index,typat,xred,spinat)'
1975 : !do jj=1,natrd
1976 : !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj)
1977 : !end do
1978 : !ENDDEBUG
1979 :
1980 1588 : curat=0
1981 :
1982 : !Cycle over all the symmetry operations
1983 23822 : do ii=1,nsym
1984 :
1985 : ! Cycle over all the atoms in the asymmetric unit cell
1986 49344 : do jj=1,natrd
1987 :
1988 : ! Symmetry operation application
1989 714616 : bckat(:)=matmul(symrel(:,:,ii),xred(:,jj))+tnons(:,ii)
1990 :
1991 : ! Normalization of the coordinates in [0,1)
1992 102088 : do iij=1,3
1993 101094 : do while (bckat(iij)<-tolsym)
1994 101094 : bckat(iij)=bckat(iij)+1.0d0
1995 : end do
1996 104263 : do while (bckat(iij)>=1.0d0-tolsym)
1997 2175 : bckat(iij)=bckat(iij)-1.0d0
1998 : end do
1999 : end do
2000 :
2001 : ! Check for duplicate atoms
2002 25522 : flagch=0
2003 324532 : do kk=1,curat
2004 300417 : flageq=0
2005 : if ( abs(bckxred(1,kk)-bckat(1))<tolsym .and. &
2006 300417 : & abs(bckxred(2,kk)-bckat(2))<tolsym .and. &
2007 : & abs(bckxred(3,kk)-bckat(3))<tolsym ) exit
2008 324532 : flagch=flagch+1
2009 : end do
2010 :
2011 47756 : if (flagch==curat) then
2012 : ! Add the obtained atom to the bckxred list
2013 24115 : curat=curat+1
2014 96460 : bckxred(:,curat)=bckat
2015 24115 : bcktypat(curat)=typat(jj)
2016 24115 : bckchrgat(curat)=chrgat(jj)
2017 96460 : bcknucdipmom(:,curat)=nucdipmom(:,jj)
2018 96460 : bckspinat(:,curat)=spinat(:,jj)*symafm(ii)
2019 96460 : bckspinat_cart(:,curat)=spinat_cart(:,jj)*symafm(ii)
2020 : end if
2021 :
2022 : end do
2023 : end do
2024 :
2025 : !DEBUG
2026 : !write(std_out,*)' fillcell : Proposed coordinates ='
2027 : !do ii=1,curat
2028 : !write(std_out,'(i4,3es16.6)' )ii,bckxred(:,ii)
2029 : !end do
2030 : !ENDDEBUG
2031 :
2032 1588 : if (curat>natom) then
2033 : write(msg, '(a,i3,a,a,i7,a,a,a,a)' )&
2034 0 : & 'The number of atoms obtained from symmetries, ',curat,ch10,&
2035 0 : & 'is greater than the input number of atoms, natom=',natom,ch10,&
2036 0 : & 'This is not allowed.',ch10,&
2037 0 : & 'Action: modify natom or the symmetry data in the input file.'
2038 0 : ABI_ERROR(msg)
2039 : end if
2040 :
2041 1588 : if (curat<natom) then
2042 : write(msg, '(a,i3,a,a,i7,a,a,a,a)' )&
2043 0 : & 'The number of atoms obtained from symmetries, ',curat,ch10,&
2044 0 : & 'is lower than the input number of atoms, natom=',natom,ch10,&
2045 0 : & 'This is not allowed.',ch10,&
2046 0 : & 'Action: modify natom or the symmetry data in the input file.'
2047 0 : ABI_ERROR(msg)
2048 : end if
2049 :
2050 : !Assignment of symmetry to xred
2051 98048 : xred(:,1:natom)=bckxred(:,1:natom)
2052 25703 : typat(1:natom)=bcktypat(1:natom)
2053 25703 : chrgat(1:natom)=bckchrgat(1:natom)
2054 98048 : nucdipmom(1:3,1:natom)=bcknucdipmom(1:3,1:natom)
2055 98048 : spinat(1:3,1:natom)=bckspinat(1:3,1:natom)
2056 98048 : spinat_cart(1:3,1:natom)=bckspinat_cart(1:3,1:natom)
2057 :
2058 : !DEBUG
2059 : !write(std_out,*)' fillcell : exit with natom=',natom
2060 : !write(std_out,*)' Describe the output atoms (index,typat,xred,spinat)'
2061 : !do jj=1,natom
2062 : !write(std_out,'(i3,2x,i3,6es12.2)')jj,typat(jj),xred(:,jj),spinat(:,jj)
2063 : !end do
2064 : !ENDDEBUG
2065 :
2066 1588 : end subroutine fillcell
2067 : !!***
2068 :
2069 : !!****f* m_ingeo/invacuum
2070 : !!
2071 : !! NAME
2072 : !! invacuum
2073 : !!
2074 : !! FUNCTION
2075 : !! Determine whether there is vacuum along some of the primitive directions in real space.
2076 : !!
2077 : !! INPUTS
2078 : !! jdtset=number of the dataset looked for
2079 : !! lenstr=actual length of the string
2080 : !! natom=number of atoms
2081 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
2082 : !! string*(*)=character string containing all the input data.
2083 : !! Initialized previously in instrng.
2084 : !! xred(3,natom)=reduced dimensionless atomic coordinates
2085 : !!
2086 : !! OUTPUT
2087 : !! vacuum(3)= for each direction, 0 if no vacuum, 1 if vacuum
2088 : !!
2089 : !! SOURCE
2090 :
2091 13932 : subroutine invacuum(jdtset,lenstr,natom,rprimd,string,vacuum,xred)
2092 :
2093 : !Arguments ------------------------------------
2094 : !scalars
2095 : integer,intent(in) :: jdtset,lenstr,natom
2096 : character(len=*),intent(in) :: string
2097 : !arrays
2098 : integer,intent(out) :: vacuum(3)
2099 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2100 :
2101 : !Local variables-------------------------------
2102 : !scalars
2103 : integer :: ia,ii,marr,tread
2104 : real(dp) :: max_diff_xred,ucvol,vacwidth,vacxred
2105 : !arrays
2106 13932 : integer,allocatable :: list(:)
2107 13932 : integer,allocatable :: intarr(:)
2108 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2109 13932 : real(dp),allocatable :: xred_sorted(:)
2110 13932 : real(dp),allocatable :: dprarr(:)
2111 :
2112 : ! *************************************************************************
2113 :
2114 : !Compute the maximum size of arrays intarr and dprarr
2115 13932 : marr=3
2116 13932 : ABI_MALLOC(intarr,(marr))
2117 13932 : ABI_MALLOC(dprarr,(marr))
2118 :
2119 : !Get metric quantities
2120 13932 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2121 :
2122 : !Read vacwidth, or set the default
2123 13932 : vacwidth=10.0_dp
2124 13932 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacwidth',tread,'LEN')
2125 13932 : if(tread==1) vacwidth=dprarr(1)
2126 :
2127 : !Read vacuum, or compute it using the atomic coordinates and vacwidth.
2128 13932 : vacuum(1:3)=0
2129 13932 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'vacuum',tread,'INT')
2130 :
2131 13932 : if(tread==1)then
2132 8 : vacuum(1:3)=intarr(1:3)
2133 : else
2134 : ! For each direction, determine whether a vacuum space exists
2135 41790 : ABI_MALLOC(list,(natom))
2136 41790 : ABI_MALLOC(xred_sorted,(natom))
2137 55720 : do ii=1,3
2138 : ! This is the minimum xred difference needed to have vacwidth
2139 167160 : vacxred=vacwidth*sqrt(sum(gprimd(:,ii)**2))
2140 : ! Project the reduced coordinate in the [0.0_dp,1.0_dp[ interval
2141 270222 : xred_sorted(:)=mod(xred(ii,:),1.0_dp)
2142 : ! list is dummy
2143 270222 : list(:)=0
2144 : ! Sort xred_sorted
2145 41790 : call sort_dp(natom,xred_sorted,list,tol14)
2146 41790 : if(natom==1)then
2147 : max_diff_xred=1.0_dp
2148 : else
2149 : ! Compute the difference between each pair of atom in the sorted order
2150 : max_diff_xred=0.0_dp
2151 220452 : do ia=1,natom-1
2152 220452 : max_diff_xred=max(max_diff_xred,xred_sorted(ia+1)-xred_sorted(ia))
2153 : end do
2154 : ! Do not forget the image of the first atom in the next cell
2155 33810 : max_diff_xred=max(max_diff_xred,1.0_dp+xred_sorted(1)-xred_sorted(ia))
2156 : end if
2157 55720 : if(vacxred<max_diff_xred+tol10)vacuum(ii)=1
2158 : end do
2159 13930 : ABI_FREE(list)
2160 13930 : ABI_FREE(xred_sorted)
2161 : end if
2162 :
2163 : !DEBUG
2164 : !write(std_out,*)' invacuum : vacuum=',vacuum(1:3)
2165 : !ENDDEBUG
2166 :
2167 13932 : ABI_FREE(intarr)
2168 13932 : ABI_FREE(dprarr)
2169 :
2170 13932 : end subroutine invacuum
2171 : !!***
2172 :
2173 : !!****f* m_ingeo/checkspvec
2174 : !!
2175 : !! NAME
2176 : !! checkspvec
2177 : !!
2178 : !! FUNCTION
2179 : !! Check the consistency of spin-related input vectors (such as spinat or
2180 : !! hspinfield) with the spin quantization axis (spinaxis), and convert
2181 : !! Cartesian coordnate to the local spinaxis coordinate when needed
2182 : !!
2183 : !! INPUTS
2184 : !! name=character string identifying the spin-related quantity: 'spinat' or 'hspinfield'
2185 : !! nitem=number of spin vectors to be treated: natom for spinat, 1 for hspinfield
2186 : !! spinaxis_in(3)=spin quantization axis
2187 : !! tread=integer flag (0 or 1), set to 1 if the local-frame quantity
2188 : !! (spinat or hspinfield) is provided in the input
2189 : !! tread_cart=integer flag (0 or 1), set to 1 if the Cartesian quantity
2190 : !! (spinat_cart or hspinfield_cart) is provided in the input
2191 : !! vec_local(nvec,nitem)=spin vectors defined in the local spin reference coordinates
2192 : !! vec_cart(nvec,nitem)=spin vectors defined in Cartesian coordinates
2193 : !!
2194 : !! OUTPUT
2195 : !! vec_local=updated consistent spin vectors in local spin reference coordinate
2196 : !! vec_cart=updated consistent spin vectors in cartesian coordinate
2197 : !!
2198 : !! SOURCE
2199 :
2200 14206 : subroutine checkspvec(name,nitem,spinaxis_in,tread,tread_cart,vec_local,vec_cart)
2201 :
2202 : !Arguments ------------------------------------
2203 : !scalars
2204 : integer,intent(in) :: nitem, tread, tread_cart
2205 : character(len=*),intent(in) :: name
2206 : !arrays
2207 : real(dp),intent(in) :: spinaxis_in(3)
2208 : real(dp),intent(inout) :: vec_local(3,nitem), vec_cart(3,nitem)
2209 :
2210 : !Local variables-------------------------------
2211 : !scalars
2212 : integer :: i, icase
2213 : real(dp) :: alpha, beta, norm
2214 : !arrays
2215 : real(dp) :: spinaxis(3), R(3,3), RT(3,3), v_cart(3), v_local(3)
2216 : logical :: trivial_axis
2217 : character(len=1000) :: msg
2218 :
2219 : ! *************************************************************************
2220 :
2221 : ! case 0: nothing provided
2222 15913 : if (tread == 0 .and. tread_cart == 0) return
2223 :
2224 1707 : spinaxis(:) = spinaxis_in(:)
2225 6828 : norm = sqrt(dot_product(spinaxis, spinaxis))
2226 1707 : if (norm <= tol8) then
2227 0 : spinaxis(:) = [zero, zero, one]
2228 0 : trivial_axis = .true.
2229 : else
2230 6828 : spinaxis(:) = spinaxis(:) / norm
2231 6825 : trivial_axis = all(abs(spinaxis(:) - [zero, zero, one]) < tol8)
2232 : end if
2233 :
2234 1707 : R(:,:) = zero
2235 1707 : R(1,1) = one; R(2,2) = one; R(3,3) = one
2236 :
2237 1707 : if (.not. trivial_axis) then
2238 1 : call geteuler(spinaxis, alpha, beta)
2239 1 : call cart2spinaxis(alpha, beta, R)
2240 : end if
2241 22191 : RT(:,:) = transpose(R)
2242 :
2243 : ! 3 cases: (tread, tread_cart)
2244 1707 : icase = 2*tread + tread_cart
2245 :
2246 : select case (icase)
2247 :
2248 : case (1) ! case 1: only _cart provided (tread = 0, tread_cart = 1)
2249 2 : do i = 1, nitem
2250 17 : vec_local(:,i) = matmul(R, vec_cart(:,i))
2251 : end do
2252 1706 : return
2253 :
2254 : case (2) ! case 2: only local provided (tread = 1, tread_cart = 0)
2255 1706 : if (.not. trivial_axis) then
2256 : write(msg,'(a)') 'Spinaxis is defined, but only ' // trim(name) // ' is present.' // ch10 // &
2257 : 'Action: please use ' // trim(name) // '_cart (Cartesian) instead of ' // trim(name) // &
2258 0 : 'when spinaxis is set.'
2259 0 : ABI_ERROR(msg)
2260 : end if
2261 3842 : do i = 1, nitem
2262 10250 : vec_cart(:,i) = vec_local(:,i)
2263 : end do
2264 : return
2265 :
2266 : case (3) ! case 3: both provided and check consistency: (tread = tread_cart = 1)
2267 0 : do i = 1, nitem
2268 0 : v_cart(:) = matmul(RT, vec_local(:,i))
2269 0 : v_local(:) = matmul(R, vec_cart(:,i))
2270 :
2271 0 : if (maxval(abs(v_local(:) - vec_local(:, i))) > tol8 .or. &
2272 0 : maxval(abs(v_cart(:) - vec_cart(:, i))) > tol8) then
2273 : write(msg,'(a,a,a,a,a,a,a,3(1x,es16.8),a,a,3(1x,es16.8),a,a,3(1x,es16.8),a,a)') &
2274 0 : 'Both ', trim(name), ' and ', trim(name), '_cart are set but inconsistent.', ch10, &
2275 0 : 'spinaxis =', spinaxis_in(1), spinaxis_in(2), spinaxis_in(3), ch10, &
2276 0 : trim(name)//' (local) =', vec_local(1,i), vec_local(2,i), vec_local(3,i), ch10, &
2277 0 : trim(name)//'_cart (cart) =', vec_cart(1,i), vec_cart(2,i), vec_cart(3,i), ch10, &
2278 0 : 'Action: make them consistent or provide only one of them.'
2279 0 : ABI_ERROR(msg)
2280 : end if
2281 : end do
2282 1707 : return
2283 :
2284 : end select
2285 :
2286 : end subroutine checkspvec
2287 : !!***
2288 :
2289 : end module m_ingeo
2290 : !!***
|
