Branch data 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 : 7284 : subroutine ingeo (acell,amu,atndlist,bravais,chrgat,dtset,field_red,field_red_axial,&
138 : 7284 : genafm,iatfix,iatnd,icoulomb,iimage,iout,jdtset,jellslab,lenstr,mixalch,&
139 : : msym,natnd,natom,nimage,npsp,npspalch,nspden,nsym,ntypalch,ntypat,&
140 : 7284 : nucdipmom,nzchempot,pawspnorb,&
141 : 7284 : ptgroupma,ratsph,rprim,slabzbeg,slabzend,spgroup,spinat,spinat_cart,string,supercell_lattice,symafm,&
142 : 7284 : 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 : 7284 : character(len=lenstr) :: geo_string
183 : : type(atomdata_t) :: atom
184 : 7284 : type(geo_t) :: geo
185 : : !arrays
186 : : integer :: bravais_reduced(11)
187 : 7284 : integer,allocatable :: intarr(:)
188 : 7284 : integer,allocatable :: is_translation(:)
189 : 7284 : 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 : 7284 : real(dp),allocatable :: mass_psp(:),tnons_cart(:,:),tnons_new(:,:),translations(:,:)
193 : 7284 : 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 : 7284 : marr=max(12,3*natom,9*msym)
203 [ - + + - : 21852 : ABI_MALLOC(intarr,(marr))
+ - - + -
+ ]
204 [ - + + - : 21852 : ABI_MALLOC(dprarr,(marr))
+ - - + -
+ ]
205 : :
206 : : ! Try from geo_string
207 : 7284 : call intagm(dprarr, intarr, jdtset, marr, 1, string(1:lenstr), 'structure', tread_geo, 'KEY', key_value=geo_string)
208 : :
209 [ + + ]: 7284 : 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 : 7229 : 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 [ + + ]: 29136 : scalecart(1:3)=one
223 : 7284 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'scalecart',tscalecart,'LEN')
224 [ + + + + ]: 7314 : if(tscalecart==1) scalecart(1:3)=dprarr(1:3)
225 : 7284 : call intagm_img(scalecart,iimage,jdtset,lenstr,nimage,3,string,"scalecart",tscalecart,'LEN')
226 : :
227 [ + + ]: 29136 : rprim(:,1)=scalecart(:)*rprim(:,1)
228 [ + + ]: 29136 : rprim(:,2)=scalecart(:)*rprim(:,2)
229 [ + + ]: 29136 : rprim(:,3)=scalecart(:)*rprim(:,3)
230 [ + + ]: 29136 : scalecart(:)=one
231 : :
232 : : ! Compute the multiplicity of the supercell
233 : 7284 : multiplicity=supercell_lattice(1)*supercell_lattice(2)*supercell_lattice(3)
234 : :
235 [ + + ]: 7284 : if (tread_geo == 0) then
236 : : ! Get the number of atom in the unit cell. Read natom from string
237 : 7229 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natom',tread,'INT')
238 : :
239 : : ! Might initialize natom from XYZ file
240 [ + + ]: 7229 : if (tread==0) call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'_natom',tread,'INT')
241 : :
242 [ + - ]: 7229 : 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 : 7284 : call mkrdim(acell,rprim,rprimd_read)
249 : :
250 : : ! Multiply the rprim to get the rprim of the supercell
251 [ + + ]: 7284 : 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 : 7284 : call mkrdim(acell, rprim, rprimd)
259 : 7284 : call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
260 : :
261 [ + - + - : 7284 : 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 [ + + ]: 29136 : do ii = 1, 3
286 [ + + + + : 94692 : 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 : 7284 : 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 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'tolsym',tread,'DPR')
300 [ + + ]: 7284 : 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 [ - + + - : 21852 : ABI_MALLOC(ptsymrel,(3,3,msym))
+ - + - -
+ - + ]
306 : 7284 : call symlatt(bravais,dev_null,msym,nptsym,ptsymrel,rprimd,tolsym)
307 : :
308 : : ! 3) Possibly, initialize a jellium slab
309 : 7284 : jellslab=0
310 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'jellslab',tread,'INT')
311 [ + + ]: 7284 : if(tread==1) jellslab=intarr(1)
312 : :
313 : 7284 : slabzbeg=zero
314 : 7284 : slabzend=zero
315 [ + + ]: 7284 : 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 : 7284 : natrd=natom
326 [ + + ]: 7284 : if(multiplicity > 1) natrd = natom_uc
327 : :
328 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natrd',tnatrd,'INT')
329 [ + + ]: 7284 : if(tnatrd==1) natrd=intarr(1)
330 : :
331 [ + - - + ]: 7284 : 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 [ - + + - : 21852 : ABI_MALLOC(typat_read,(natrd))
+ - + - -
+ - + ]
356 : 7284 : typat_read(1)=1
357 : :
358 [ + + ]: 7284 : if (tread_geo == 0) then
359 : 7229 : call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'typat',tread,'INT')
360 : :
361 : : ! If not read, try the XYZ data
362 [ + + ]: 7229 : if(tread==0) call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'_typat',tread,'INT')
363 [ + + + + ]: 23367 : if(tread==1) typat_read(1:natrd)=intarr(1:natrd)
364 : :
365 : : else
366 [ - + - - : 385 : typat_read = geo%typat
+ + ]
367 : : end if
368 : :
369 [ + + ]: 23721 : do iatom=1,natrd
370 [ + - - + ]: 23721 : 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 [ - + + - : 21852 : ABI_MALLOC(xcart_read,(3,natrd))
+ - + - -
+ - + ]
382 [ - + + - : 14568 : ABI_MALLOC(xred_read,(3,natrd))
- + - + ]
383 : :
384 : 7284 : random_atpos=0
385 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'random_atpos',txrandom,'INT')
386 [ + + ]: 7284 : if(txrandom==1) random_atpos=intarr(1)
387 [ - + ]: 7284 : 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 : 7284 : 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 [ + + ]: 7284 : if (tread_geo /= 0) then
402 : 55 : txcart = 0; txrandom = 0; txred = 1
403 [ + - + - : 1210 : xred_read = geo%xred
- + - - -
- - - + +
+ + ]
404 : :
405 : : else
406 : 7229 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xred',txred,'DPR')
407 [ + + + - : 19975 : if (txred==1 .and. txrandom == 0) xred_read(:,1:natrd) = reshape(dprarr(1:3*natrd) , [3, natrd])
+ + ]
408 : 7229 : call intagm_img(xred_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xred",txred,'DPR')
409 : :
410 : 7229 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'xcart',txcart,'LEN')
411 [ + + + - : 8919 : if (txcart==1 .and. txrandom==0) xcart_read(:,1:natrd) = reshape(dprarr(1:3*natrd), [3, natrd])
+ + ]
412 : 7229 : call intagm_img(xcart_read,iimage,jdtset,lenstr,nimage,3,natrd,string,"xcart",txcart,'LEN')
413 : :
414 : : ! Might initialize xred from XYZ file
415 [ + + ]: 7229 : 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 [ - + ]: 7284 : 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 [ + + ]: 7284 : if (txred==1) write(msg, '(a)' ) ' xred is defined in input file'
441 [ + + ]: 7284 : if (txcart ==1) write(msg, '(a)' ) ' xcart is defined in input file (possibly in Angstrom)'
442 [ + + ]: 7284 : if (txrandom ==1) write(msg, '(a)' ) ' xred as random positions in the unit cell'
443 [ + + ]: 7284 : if (txrandom ==1) write(msg, '(a)' ) ' xcart are defined from a random distribution '
444 : 7284 : call wrtout(std_out, msg)
445 : :
446 [ - + ]: 7284 : 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 [ + + + + ]: 7284 : if (txred==1 .or. txrandom /=0 ) then
454 : 6438 : call wrtout(std_out,' ingeo: takes atomic coordinates from input array xred ')
455 : 6438 : 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 [ - + + - : 21852 : ABI_MALLOC(xcart,(3,natom))
+ - + - -
+ - + ]
466 : :
467 : : !7) Eventually read the symmetries
468 : : !Take care of the symmetries
469 : :
470 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nsym',tread,'INT')
471 [ + + ]: 7284 : if(tread==1) nsym=intarr(1)
472 : :
473 : : ! Check that nsym is not negative
474 [ - + ]: 7284 : 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 [ - + ]: 7284 : 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 [ + + ]: 7284 : 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 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'symmorphi',tread,'INT')
494 [ + + ]: 7284 : if(tread==1) symmorphi=intarr(1)
495 : :
496 : : ! Now, read the symmetry operations
497 [ + + ]: 7284 : if(nsym>0)then
498 : 1284 : call intagm(dprarr,intarr,jdtset,marr,9*nsym,string(1:lenstr),'symrel',tread,'INT')
499 [ + + - + ]: 1284 : 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 [ + + + + ]: 3144 : if(tread==1) symrel(:,:,1:nsym)=reshape( intarr(1:9*nsym) , [3, 3, nsym])
506 : :
507 : : ! Take care of tnons
508 [ + + + + ]: 41796 : tnons(:,1:nsym)=zero
509 : 1284 : call intagm(dprarr,intarr,jdtset,marr,3*nsym,string(1:lenstr),'tnons',tread,'DPR')
510 [ + + + + ]: 2362 : if(tread==1) tnons(:,1:nsym)=reshape( dprarr(1:3*nsym), [3, nsym])
511 : :
512 [ + + ]: 1284 : if(symmorphi==0)then
513 [ + + ]: 74 : do isym=1,nsym
514 [ + + - + ]: 185 : 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 : 1284 : call intagm(dprarr,intarr,jdtset,marr,nsym,string(1:lenstr),'symafm',tread,'INT')
527 [ + + + + ]: 1332 : 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 : 7284 : nobj=0
538 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nobj',tread,'INT')
539 [ + + ]: 7284 : if(tread==1) nobj=intarr(1)
540 [ - + ]: 7284 : 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 : 7284 : expert_user=0
551 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'expert_user',tread,'INT')
552 [ + + ]: 7284 : if(tread==1) expert_user=intarr(1)
553 [ - + ]: 25 : if(expert_user==0)then
554 : 7259 : chkprim=1
555 : 7259 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'chkprim',tread,'INT')
556 [ + + ]: 7259 : if(tread==1) chkprim=intarr(1)
557 : : else
558 : 25 : chkprim=0
559 : : endif
560 : :
561 [ + + ]: 7284 : 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 : 7280 : call intagm(dprarr,intarr,jdtset,marr,natrd,string(1:lenstr),'chrgat',tread,'DPR')
624 [ + + + + ]: 7342 : if(tread==1)chrgat(1:natrd) = dprarr(1:natrd)
625 : :
626 : : ! Spinat is read for each irreducible atom, from 1 to natrd
627 : 7280 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'spinat',tread,'DPR')
628 [ + + ]: 7280 : if(tread==1) then
629 [ + + ]: 5106 : spinat(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natrd])
630 [ + + + + ]: 10254 : dtset%spinat_in(1:3,1:natrd) = spinat(1:3,1:natrd)
631 : : end if
632 : :
633 : 7280 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'spinat_cart',tread_cart,'DPR')
634 [ + + + + ]: 7282 : if(tread_cart==1) spinat_cart(1:3,1:natrd) = reshape( dprarr(1:3*natrd) , [3, natom])
635 : 7280 : call checkspvec('spinat',natrd,dtset%spinaxis,tread,tread_cart,spinat,spinat_cart)
636 [ + + + + : 7284 : 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 [ + + + + ]: 161296 : nucdipmom=zero
640 [ + + ]: 7280 : 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 : 7266 : call intagm(dprarr,intarr,jdtset,marr,3*natrd,string(1:lenstr),'nucdipmom',tread,'DPR')
650 [ + + + + ]: 7336 : 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 [ + + ]: 7280 : 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 : 7279 : call xcart2xred(natrd,rprimd,xcart_read,xred)
674 : : end if
675 : :
676 : 7280 : spgroup=0
677 : 7280 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'spgroup',tread,'INT')
678 [ + + ]: 7280 : if(tread==1) spgroup=intarr(1)
679 : :
680 [ + + + + ]: 7280 : if(spgroup/=0 .or. nsym/=0)then
681 : :
682 [ - + - - : 2903 : 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 [ - + - - : 2903 : 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 [ + + ]: 7614 : typat(1:natrd)=typat_read(1:natrd)
699 : :
700 [ + + - + ]: 2903 : 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 : 2903 : brvltt=0
712 : :
713 [ + + ]: 2903 : 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 [ + + + + ]: 2903 : 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 : 2903 : iexit=0
869 : :
870 : 2903 : 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 [ - + - - ]: 4377 : 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 [ + - + + ]: 16014 : 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 : 4377 : noncoll=0; if (nspden == 4) noncoll=1
897 : :
898 : 4377 : use_inversion=1
899 [ + + + + : 4377 : if (dtset%usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then
+ + ]
900 : 136 : 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 : 136 : 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 [ + + ]: 4417 : do try_primitive=1,1
938 : :
939 [ + + + + ]: 4377 : 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 : 4377 : chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=dtset%field_red,invaraxial_red=dtset%field_red_axial,invar_z=invar_z)
943 : :
944 : 4377 : chkprim_fake=-1
945 [ - + + - : 13131 : ABI_MALLOC(is_translation,(nsym))
+ - + - -
+ - + ]
946 : 4377 : call chkprimit(chkprim_fake, multi, nsym, symafm, symrel, is_translation)
947 : :
948 [ + + ]: 4417 : 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 : 77 : & ' Not all components of symrel are integers in the primitive cell coordinate system.'
982 : 77 : ABI_COMMENT(msg)
983 : 77 : exit
984 : : endif
985 : : else ! The cell is primitive
986 [ - + ]: 4260 : ABI_FREE(is_translation)
987 : 4260 : 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 : 7280 : call xred2xcart(natom,rprimd,xcart,xred)
996 : :
997 : : end if ! check of existence of an object
998 : :
999 [ - + ]: 7284 : ABI_FREE(ptsymrel)
1000 [ - + ]: 7284 : ABI_FREE(xcart_read)
1001 [ - + ]: 7284 : ABI_FREE(xcart)
1002 [ - + ]: 7284 : ABI_FREE(xred_read)
1003 [ - + ]: 7284 : ABI_FREE(typat_read)
1004 : :
1005 : 7284 : call geo%free()
1006 : :
1007 : : ! Correct the default nsym value, if a symmetry group has not been generated.
1008 [ - + ]: 7284 : if (nsym==0) nsym=1
1009 : :
1010 : : !--------------------------------------------------------------------------------------------------------
1011 : :
1012 : 7284 : icoulomb=0
1013 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'icoulomb',tread,'INT')
1014 [ - + ]: 7284 : 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 [ - + ]: 7284 : 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 : 7284 : 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 [ + + ]: 7284 : if(tolsym>1.00001e-8)then
1057 : : ! Check whether the symmetry operations are consistent with the lattice vectors
1058 : 7256 : iexit=1
1059 : :
1060 : 7256 : call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tol8)
1061 : :
1062 [ + + ]: 7256 : 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 : 7284 : 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 [ + + ]: 7284 : if(symmorphi==0)then
1101 : 101 : jsym=0
1102 [ + + ]: 2066 : do isym=1,nsym
1103 [ + + + + ]: 7961 : if(sum(tnons(:,isym)**2)<tol6)then
1104 : 1185 : jsym=jsym+1
1105 : : ! This symmetry operation is non-symmorphic, and can be kept
1106 [ + + ]: 1185 : 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 : 101 : 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 [ + + + + ]: 162656 : iatfix(:,:)=0
1125 : :
1126 [ + + ]: 36420 : do idir=0,3
1127 : :
1128 [ + + ]: 29136 : if(idir==0)then
1129 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfix',tread,'INT')
1130 [ + + ]: 21852 : else if(idir==1)then
1131 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixx',tread,'INT')
1132 [ + + ]: 14568 : else if(idir==2)then
1133 : 7284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natfixy',tread,'INT')
1134 [ + - ]: 7284 : else if(idir==3)then
1135 : 7284 : 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 : 29136 : natfix=0
1140 [ + + ]: 29136 : if(tread==1) natfix=intarr(1)
1141 : :
1142 : : ! Check the validity of natfix
1143 [ + - - + ]: 29136 : 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 [ + + ]: 29136 : if(idir==0)then
1153 : 7284 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfix',tread,'INT')
1154 [ + + ]: 21852 : else if(idir==1)then
1155 : 7284 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixx',tread,'INT')
1156 [ + + ]: 14568 : else if(idir==2)then
1157 : 7284 : call intagm(dprarr,intarr,jdtset,marr,natfix,string(1:lenstr),'iatfixy',tread,'INT')
1158 [ + - ]: 7284 : else if(idir==3)then
1159 : 7284 : 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 [ + + - + ]: 29136 : 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 [ + + - + ]: 29136 : 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 [ + + ]: 36420 : 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 [ + + + + ]: 162656 : vel(:,:)=zero
1203 : 7284 : call intagm(dprarr,intarr,jdtset,marr,3*natom,string(1:lenstr),'vel',tread,'DPR')
1204 [ + + + + ]: 7290 : if(tread==1)vel(:,:)=reshape( dprarr(1:3*natom), [3, natom])
1205 : 7284 : call intagm_img(vel,iimage,jdtset,lenstr,nimage,3,natom,string,"vel",tread,'DPR')
1206 : :
1207 : 7284 : vel_cell(:,:)=zero
1208 : 7284 : call intagm(dprarr,intarr,jdtset,marr,3*3,string(1:lenstr),'vel_cell',tread,'DPR')
1209 [ - + ]: 7284 : if(tread==1)vel_cell(:,:)=reshape( dprarr(1:9), [3,3])
1210 : 7284 : call intagm_img(vel_cell,iimage,jdtset,lenstr,nimage,3,3,string,"vel_cell",tread,'DPR')
1211 : :
1212 : : ! mixalch
1213 [ + + ]: 7284 : 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 [ - + + - : 21852 : ABI_MALLOC(mass_psp,(npsp))
+ - + - -
+ - + ]
1233 [ + + ]: 16464 : do ipsp=1,npsp
1234 : 9180 : call atomdata_from_znucl(atom,znucl(ipsp))
1235 : 9180 : amu_default = atom%amu
1236 : 16464 : mass_psp(ipsp)=amu_default
1237 : : end do
1238 : : ! When the pseudo-atom is pure, simple copy
1239 : 7284 : ntyppure=ntypat-ntypalch
1240 [ + + ]: 7284 : if(ntyppure>0)then
1241 [ + + ]: 16402 : amu(1:ntyppure)=mass_psp(1:ntyppure)
1242 : : end if
1243 : : ! When the pseudo-atom is alchemical, must make mixing
1244 [ + + ]: 7284 : 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 [ - + ]: 7284 : ABI_FREE(mass_psp)
1253 : :
1254 : 7284 : call intagm(dprarr,intarr,jdtset,marr,ntypat,string(1:lenstr),'amu',tread,'DPR')
1255 [ + + + + ]: 7816 : if(tread==1)amu(:)=dprarr(1:ntypat)
1256 : 7284 : call intagm_img(amu,iimage,jdtset,lenstr,nimage,ntypat,string,"amu",tread,'DPR')
1257 : :
1258 [ - + ]: 7284 : ABI_FREE(intarr)
1259 [ - + ]: 7284 : ABI_FREE(dprarr)
1260 : :
1261 : : !DEBUG
1262 : : !write(std_out,'(a)')' m_ingeo%ingeo : exit '
1263 : : !call flush(std_out)
1264 : : !ENDDEBUG
1265 : :
1266 [ - + - + : 14568 : 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 : 14034 : 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 : 14034 : integer,allocatable :: list(:)
2107 : 14034 : integer,allocatable :: intarr(:)
2108 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2109 : 14034 : real(dp),allocatable :: xred_sorted(:)
2110 : 14034 : real(dp),allocatable :: dprarr(:)
2111 : :
2112 : : ! *************************************************************************
2113 : :
2114 : : !Compute the maximum size of arrays intarr and dprarr
2115 : 14034 : marr=3
2116 [ - + ]: 14034 : ABI_MALLOC(intarr,(marr))
2117 [ - + ]: 14034 : ABI_MALLOC(dprarr,(marr))
2118 : :
2119 : : !Get metric quantities
2120 : 14034 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2121 : :
2122 : : !Read vacwidth, or set the default
2123 : 14034 : vacwidth=10.0_dp
2124 : 14034 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'vacwidth',tread,'LEN')
2125 [ + + ]: 14034 : if(tread==1) vacwidth=dprarr(1)
2126 : :
2127 : : !Read vacuum, or compute it using the atomic coordinates and vacwidth.
2128 : 14034 : vacuum(1:3)=0
2129 : 14034 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'vacuum',tread,'INT')
2130 : :
2131 [ + + ]: 14034 : 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 [ - + + - : 42096 : ABI_MALLOC(list,(natom))
+ - + - -
+ - + ]
2136 [ - + + - : 42096 : ABI_MALLOC(xred_sorted,(natom))
+ - + - -
+ - + ]
2137 [ + + ]: 56128 : do ii=1,3
2138 : : ! This is the minimum xred difference needed to have vacwidth
2139 [ + + ]: 168384 : vacxred=vacwidth*sqrt(sum(gprimd(:,ii)**2))
2140 : : ! Project the reduced coordinate in the [0.0_dp,1.0_dp[ interval
2141 [ + + ]: 271080 : xred_sorted(:)=mod(xred(ii,:),1.0_dp)
2142 : : ! list is dummy
2143 [ + + ]: 271080 : list(:)=0
2144 : : ! Sort xred_sorted
2145 : 42096 : call sort_dp(natom,xred_sorted,list,tol14)
2146 [ + + ]: 42096 : 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 [ + + ]: 220920 : do ia=1,natom-1
2152 : 220920 : 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 : 34032 : max_diff_xred=max(max_diff_xred,1.0_dp+xred_sorted(1)-xred_sorted(ia))
2156 : : end if
2157 [ + + ]: 56128 : if(vacxred<max_diff_xred+tol10)vacuum(ii)=1
2158 : : end do
2159 [ - + ]: 14032 : ABI_FREE(list)
2160 [ - + ]: 14032 : ABI_FREE(xred_sorted)
2161 : : end if
2162 : :
2163 : : !DEBUG
2164 : : !write(std_out,*)' invacuum : vacuum=',vacuum(1:3)
2165 : : !ENDDEBUG
2166 : :
2167 [ - + ]: 14034 : ABI_FREE(intarr)
2168 [ - + ]: 14034 : ABI_FREE(dprarr)
2169 : :
2170 : 14034 : 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 : 14308 : 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 [ + + + + ]: 16023 : if (tread == 0 .and. tread_cart == 0) return
2223 : :
2224 : 1715 : spinaxis(:) = spinaxis_in(:)
2225 [ + + ]: 6860 : norm = sqrt(dot_product(spinaxis, spinaxis))
2226 [ - + ]: 1715 : if (norm <= tol8) then
2227 : 0 : spinaxis(:) = [zero, zero, one]
2228 : 0 : trivial_axis = .true.
2229 : : else
2230 [ + + ]: 6860 : spinaxis(:) = spinaxis(:) / norm
2231 [ + + + + ]: 6857 : trivial_axis = all(abs(spinaxis(:) - [zero, zero, one]) < tol8)
2232 : : end if
2233 : :
2234 : 1715 : R(:,:) = zero
2235 : 1715 : R(1,1) = one; R(2,2) = one; R(3,3) = one
2236 : :
2237 [ + + ]: 1715 : if (.not. trivial_axis) then
2238 : 1 : call geteuler(spinaxis, alpha, beta)
2239 : 1 : call cart2spinaxis(alpha, beta, R)
2240 : : end if
2241 [ + + + + ]: 22295 : RT(:,:) = transpose(R)
2242 : :
2243 : : ! 3 cases: (tread, tread_cart)
2244 : 1715 : 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 : 1714 : return
2253 : :
2254 : : case (2) ! case 2: only local provided (tread = 1, tread_cart = 0)
2255 [ - + ]: 1714 : 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 [ + + ]: 3864 : do i = 1, nitem
2262 [ + + ]: 10314 : 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 [ + + - - ]: 1715 : return
2283 : :
2284 : : end select
2285 : :
2286 : : end subroutine checkspvec
2287 : : !!***
2288 : :
2289 : : end module m_ingeo
2290 : : !!***
|
