Branch data Line data Source code
1 : : !!****m* ABINIT/m_effective_potential_file
2 : : !! NAME
3 : : !! m_effective_potential_file
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains all routine to read the effective potential from files
7 : : !! Can also read coefficients from XML
8 : : !! (XML or DDB)
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2000-2026 ABINIT group (AM)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : :
26 : : module m_effective_potential_file
27 : :
28 : : use defs_basis
29 : : use m_errors
30 : : use m_abicore
31 : : use m_xmpi
32 : : use m_harmonics_terms
33 : : use m_anharmonics_terms
34 : : use m_effective_potential
35 : : use m_ifc
36 : : use m_ddb
37 : : use m_ddb_hdr
38 : : use netcdf
39 : :
40 : : use m_io_tools, only : open_file, get_unit
41 : : use m_geometry, only : xcart2xred, xred2xcart, metric
42 : : use m_symfind, only : symfind, symlatt
43 : : use m_crystal, only : crystal_t
44 : : use m_dynmat, only : dfpt_prtph
45 : : use m_abihist, only : abihist,abihist_init,abihist_free,abihist_copy,read_md_hist
46 : : use m_ddb_internalstr, only : ddb_internalstr
47 : :
48 : : implicit none
49 : :
50 : : public :: effective_potential_file_getDimCoeff
51 : : public :: effective_potential_file_getDimMD
52 : : public :: effective_potential_file_getDimSystem
53 : : public :: effective_potential_file_getDimStrainCoupling
54 : : public :: effective_potential_file_getType
55 : : public :: effective_potential_file_mapHistToRef
56 : : public :: effective_potential_file_read
57 : : public :: effective_potential_file_readDisplacement
58 : : public :: effective_potential_file_readMDfile
59 : : private :: coeffs_xml2effpot
60 : : private :: system_getDimFromXML
61 : : private :: system_xml2effpot
62 : : private :: system_ddb2effpot
63 : :
64 : : #ifndef HAVE_XML
65 : : private :: rdfromline
66 : : private :: rmtabfromline
67 : : private :: rdfromline_value
68 : : private :: elementfromline
69 : : #endif
70 : :
71 : : #if defined HAVE_XML
72 : : public :: effpot_xml_checkXML
73 : : public :: effpot_xml_getDimCoeff
74 : : public :: effpot_xml_readSystem
75 : : public :: effpot_xml_getValue
76 : : public :: effpot_xml_getAttribute
77 : : public :: effpot_xml_getDimSystem
78 : :
79 : : interface
80 : : subroutine effpot_xml_readSystem(filename,natom,&
81 : : & ntypat,nrpt,nqpt,amu,atmfrc,cell,dynmat,elastic_constants,&
82 : : & energy,epsilon_inf,ewald_atmfrc,&
83 : : & phfrq,rprimd,qph1l,short_atmfrc,typat,xcart,zeff)&
84 : : & bind(C,name="effpot_xml_readSystem")
85 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT
86 : : integer(C_INT) :: natom,ntypat,nrpt,nqpt
87 : : integer(C_INT) :: typat(natom)
88 : : integer(C_INT) :: cell(3,nrpt)
89 : : real(C_DOUBLE) :: energy
90 : : real(C_DOUBLE) :: dynmat(2,3,natom,3,natom,nqpt)
91 : : real(C_DOUBLE) :: phfrq(3*natom,nqpt),qph1l(3,nqpt)
92 : : real(C_DOUBLE) :: atmfrc(3,natom,3,natom,nrpt)
93 : : real(C_DOUBLE) :: short_atmfrc(3,natom,3,natom,nrpt)
94 : : real(C_DOUBLE) :: ewald_atmfrc(3,natom,3,natom,nrpt)
95 : : real(C_DOUBLE) :: amu(ntypat),rprimd(3,3),epsilon_inf(3,3)
96 : : real(C_DOUBLE) :: zeff(3,3,natom)
97 : : real(C_DOUBLE) :: elastic_constants(6,6),xcart(3,natom)
98 : : character(kind=C_CHAR) :: filename(*)
99 : : end subroutine effpot_xml_readSystem
100 : : end interface
101 : :
102 : : interface
103 : : subroutine effpot_xml_readStrainCoupling(filename,natom,&
104 : : & nrpt,voigt,elastic3rd,elastic_displacement,&
105 : : & strain_coupling,phonon_strain_atmfrc,phonon_straincell)&
106 : : & bind(C,name="effpot_xml_readStrainCoupling")
107 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT
108 : : integer(C_INT) :: natom
109 : : integer(C_INT) :: nrpt,voigt
110 : : integer(c_INT) :: phonon_straincell(3,nrpt)
111 : : real(C_DOUBLE) :: elastic3rd(6,6),elastic_displacement(6,3,natom)
112 : : real(C_DOUBLE) :: strain_coupling(3,natom)
113 : : real(C_DOUBLE) :: phonon_strain_atmfrc(3,natom,3,natom,nrpt)
114 : : character(kind=C_CHAR) :: filename(*)
115 : : end subroutine effpot_xml_readStrainCoupling
116 : : end interface
117 : :
118 : : interface
119 : : subroutine effpot_xml_readCoeff(filename,ncoeff,ndisp,nterm,&
120 : : & coefficient,atindx,cell,direction,power_disp,&
121 : : & power_strain,strain,weight)&
122 : : & bind(C,name="effpot_xml_readCoeff")
123 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT
124 : : character(kind=C_CHAR) :: filename(*)
125 : : integer(C_INT) :: ncoeff,ndisp,nterm
126 : : integer(C_INT) :: atindx(ncoeff,nterm,2,ndisp)
127 : : integer(C_INT) :: cell(ncoeff,nterm,3,2,ndisp)
128 : : integer(C_INT) :: direction(ncoeff,nterm,ndisp)
129 : : integer(C_INT) :: strain(ncoeff,nterm,ndisp)
130 : : integer(C_INT) :: power_disp(ncoeff,nterm,ndisp)
131 : : integer(C_INT) :: power_strain(ncoeff,nterm,ndisp)
132 : : real(C_DOUBLE) :: coefficient(ncoeff)
133 : : real(C_DOUBLE) :: weight(ncoeff,nterm)
134 : : end subroutine effpot_xml_readCoeff
135 : : end interface
136 : :
137 : : interface
138 : : subroutine effpot_xml_getDimSystem(filename,natom,ntypat,nqpt,nrpt1,nrpt2)&
139 : : & bind(C,name="effpot_xml_getDimSystem")
140 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_INT
141 : : integer(C_INT) :: natom,ntypat,nqpt,nrpt1,nrpt2
142 : : character(kind=C_CHAR) :: filename(*)
143 : : end subroutine effpot_xml_getDimSystem
144 : : end interface
145 : :
146 : : interface
147 : : subroutine effpot_xml_getDimStrainCoupling(filename,nrpt,voigt)&
148 : : & bind(C,name="effpot_xml_getDimStrainCoupling")
149 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_INT
150 : : integer(C_INT) :: voigt
151 : : integer(C_INT) :: nrpt
152 : : character(kind=C_CHAR) :: filename(*)
153 : : end subroutine effpot_xml_getDimStrainCoupling
154 : : end interface
155 : :
156 : : interface
157 : : subroutine effpot_xml_getDimCoeff(filename,ncoeff,nterm_max,ndisp_max)&
158 : : & bind(C,name="effpot_xml_getDimCoeff")
159 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT,C_PTR
160 : : character(kind=C_CHAR) :: filename(*)
161 : : ! character(kind=C_CHAR) :: name(*)
162 : : type(C_PTR) :: name
163 : : integer(C_INT) :: ncoeff,ndisp_max,nterm_max
164 : : end subroutine effpot_xml_getDimCoeff
165 : : end interface
166 : :
167 : :
168 : : interface
169 : : subroutine effpot_xml_checkXML(filename,name_root) &
170 : : & bind(C,name="effpot_xml_checkXML")
171 : : use, intrinsic :: iso_c_binding, only : C_CHAR
172 : : character(kind=C_CHAR) :: filename(*),name_root(*)
173 : : end subroutine effpot_xml_checkXML
174 : : end interface
175 : :
176 : : interface
177 : : subroutine effpot_xml_getValue(filename,name_value,value_result) &
178 : : & bind(C,name="effpot_xml_getValue")
179 : : use, intrinsic :: iso_c_binding, only : C_CHAR
180 : : character(kind=C_CHAR) :: filename(*),name_value(*)
181 : : character(kind=C_CHAR) :: value_result
182 : : end subroutine effpot_xml_getValue
183 : : end interface
184 : :
185 : : interface
186 : : subroutine effpot_xml_getAttribute(filename,name_value,name_attribute) &
187 : : & bind(C,name="effpot_xml_getAttribute")
188 : : use, intrinsic :: iso_c_binding, only : C_CHAR
189 : : character(kind=C_CHAR) :: filename(*),name_value(*),name_attribute(*)
190 : : end subroutine effpot_xml_getAttribute
191 : : end interface
192 : :
193 : : interface
194 : : subroutine effpot_xml_getNumberKey(filename,name_value,number) &
195 : : & bind(C,name="effpot_xml_getNumberKey")
196 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_INT
197 : : character(kind=C_CHAR) :: filename(*),name_value(*)
198 : : integer(C_INT) :: number
199 : : end subroutine effpot_xml_getNumberKey
200 : : end interface
201 : :
202 : : #endif
203 : :
204 : : integer,parameter :: XML_RECL = 50000
205 : : !!***
206 : :
207 : : CONTAINS !===========================================================================================
208 : :
209 : :
210 : : !****f* m_effective_potential_file/effective_potential_file_read
211 : : !!
212 : : !! NAME
213 : : !! effective_potential_file_read
214 : : !!
215 : : !! FUNCTION
216 : : !! tranfert file (XML or DDB) in effective potential type
217 : : !! Also transfert coefficient from xml file for ahnarmonic part
218 : : !!
219 : : !! INPUTS
220 : : !! filename = path of the file
221 : : !! hist<type(abihist)> = optional,The history of the MD (or snapshot of DFT)
222 : : !! inp<type(multibinit_dtset_type)> = optional,datatype with all the input variables (mantadory to
223 : : !! read DDB file)
224 : : !! comm=MPI communicator
225 : : !!
226 : : !! OUTPUT
227 : : !! eff_pot<type(effective_potential_type)> = datatype with all the information for effective potential
228 : : !!
229 : : !! SOURCE
230 : :
231 : 91 : subroutine effective_potential_file_read(filename,eff_pot,inp,comm,hist)
232 : :
233 : : use m_multibinit_dataset
234 : : use m_strain
235 : : use m_crystal, only : crystal_t
236 : : use m_dynmat, only : bigbx9
237 : :
238 : : !Arguments ------------------------------------
239 : : !scalars
240 : : integer,intent(in) :: comm
241 : : character(len=fnlen),intent(in) :: filename
242 : : !array
243 : : type(effective_potential_type), intent(inout) :: eff_pot
244 : : type(multibinit_dtset_type),optional,intent(in) :: inp
245 : 91 : type(ddb_type) :: ddb
246 [ + + + + : 4732 : type(crystal_t) :: Crystal
+ + + + +
+ + + + +
+ + + + ]
247 : : type(abihist),optional :: hist
248 : : !Local variables------------------------------
249 : : !scalars
250 : : integer :: ii,filetype,natom,ntypat,nqpt,nrpt
251 : : character(500) :: message
252 [ + + + + : 4732 : type(ddb_hdr_type) :: ddb_hdr
+ + + + +
+ + + + +
+ + + + ]
253 : : ! *************************************************************************
254 : :
255 : 91 : call effective_potential_file_getType(filename,filetype)
256 : :
257 [ + - ]: 91 : if (filetype/=0) then
258 : :
259 [ - + ]: 91 : if (.not.(present(inp))) then
260 : : write(message, '(4a)' )&
261 : 0 : & ' effective_potential_file_read: you need to give input file to compute ',&
262 : 0 : & 'the response fonction from DDB file ',ch10
263 : 0 : ABI_ERROR(message)
264 : : end if
265 : :
266 [ + + ]: 91 : if(filetype ==1) then
267 : : ! Read the DDB information, also perform some checks, and symmetrize partially the DDB
268 : 29 : write(message, '(3a)' )' Read the DDB information of the reference',&
269 : 58 : & ' system and perform some checks',ch10
270 : 29 : call wrtout(std_out,message,'COLL')
271 : 29 : call wrtout(ab_out,message,'COLL')
272 : :
273 : 29 : call effective_potential_file_getDimSystem(filename,comm,natom,ntypat,nqpt,nrpt)!
274 : :
275 : 29 : call ddb%from_file(filename,ddb_hdr,Crystal,comm)
276 : 29 : call ddb_hdr%free()
277 : :
278 : 29 : call ddb%set_brav(inp%brav)
279 : :
280 : : ! Transfert the ddb to the effective potential
281 : 29 : call system_ddb2effpot(Crystal,ddb,eff_pot,inp,comm)
282 : :
283 : : ! Free memory
284 : 29 : call ddb%free()
285 : 29 : call Crystal%free()
286 : :
287 : : ! Generate long rage interation for the effective potential for both type and generate supercell
288 [ + - ]: 29 : call effective_potential_generateDipDip(eff_pot,inp%dipdip_range,inp%dipdip,inp%asr,comm)
289 : :
290 : : ! If needed, print the effective potential into the output
291 [ + + ]: 29 : if (inp%prt_model>=3.or.inp%prt_model==-1) then
292 : 8 : call effective_potential_print(eff_pot,-1)
293 : : end if
294 : : end if
295 [ + + ]: 91 : if (filetype==2 .or.filetype==23) then ! xml file
296 : :
297 : : ! Free the effective potential before
298 : 38 : call effective_potential_free(eff_pot)
299 : :
300 : 38 : call system_xml2effpot(eff_pot,filename,comm,strcpling=inp%strcpling)
301 : :
302 : :
303 : : ! Assign the energy of the reference from input
304 [ - + ]: 38 : if(abs(inp%energy_reference)>tol16)then
305 : 0 : write(message,'(11a)') ch10,&
306 : 0 : & ' --- !WARNING',ch10,&
307 : 0 : & ' Energy of the reference structure is specify in ',ch10,&
308 : 0 : & ' the input file. The energy is set with',ch10,&
309 : 0 : & ' this value.',ch10,&
310 : 0 : & ' ---',ch10
311 : 0 : call wrtout(std_out,message,'COLL')
312 : 0 : eff_pot%energy = inp%energy_reference
313 : : end if
314 : :
315 : :
316 : : ! Generate long rage interation for the effective potential for both type and generate supercell
317 : 38 : call effective_potential_generateDipDip(eff_pot,inp%dipdip_range,inp%dipdip,inp%asr,comm,1)
318 : :
319 : : ! If needed, print the effective potential
320 : 38 : call effective_potential_print(eff_pot,inp%prt_model)
321 : : end if
322 [ + + ]: 91 : if (filetype==3 .or. filetype==23) then
323 : : ! Read the coefficient of the fit for the anharmonic part
324 : 28 : write(message, '(4a)' )ch10,' Read the coefficients of the polynomial fit from XML',&
325 : 56 : & ' and perform some checks',ch10
326 : 28 : call wrtout(std_out,message,'COLL')
327 : 28 : call wrtout(ab_out,message,'COLL')
328 : :
329 [ - + ]: 28 : if(eff_pot%anharmonics_terms%ncoeff/=0)then
330 : 0 : write(message,'(9a)') ch10,&
331 : 0 : & ' --- !WARNING',ch10,&
332 : 0 : & ' There is already fitted polynome set in the model',ch10,&
333 : 0 : & ' The previous coefficients will be remove',ch10,&
334 : 0 : & ' ---',ch10
335 : 0 : call wrtout(std_out,message,'COLL')
336 : : end if
337 : :
338 : 28 : call coeffs_xml2effpot(eff_pot,filename,comm)
339 : :
340 : : ! Assign the coeff number from input
341 [ + + ]: 28 : if(inp%ncoeff==0)then
342 : 20 : write(message,'(12a)') ch10,&
343 : 20 : & ' --- !WARNING',ch10,&
344 : 20 : & ' The values of the coefficients are set to 0',&
345 : 20 : & ' in the input file.',ch10,&
346 : 20 : & ' The values of the coefficients will be read in the XML',ch10,&
347 : 20 : & ' or might be fitted',ch10,&
348 : 40 : & ' ---',ch10
349 : 20 : call wrtout(std_out,message,'COLL')
350 [ + - - + : 20 : if(inp%fit_coeff <= 0 .and. &
+ + - + ]
351 : : & all(abs(eff_pot%anharmonics_terms%coefficients(:)%coefficient) <tol16)) then
352 : :
353 : 0 : write(message,'(12a)') ch10,&
354 : 0 : & ' --- !WARNING',ch10,&
355 : 0 : & ' The input for the fit process is set to 0 or -1',&
356 : 0 : & ' in the input file.',ch10,&
357 : 0 : & ' However, the values of the coefficients in the XMF files are zero,',ch10,&
358 : 0 : & ' So the coefficients can not be used',ch10,&
359 : 0 : & ' ---',ch10
360 : 0 : call wrtout(std_out,message,'COLL')
361 : :
362 : : end if
363 : : else
364 [ - + ]: 8 : if (eff_pot%anharmonics_terms%ncoeff /= inp%ncoeff)then
365 : : write(message, '(5a)' )&
366 : 0 : & ' The number of coefficients in the XML file is superior to the ',ch10,&
367 : 0 : & 'number of coefficients in the input ',ch10,&
368 : 0 : & 'Action: correct your input file or change the file'
369 : 0 : ABI_ERROR(message)
370 : : end if
371 [ + + ]: 40 : do ii = 1,eff_pot%anharmonics_terms%ncoeff
372 : : call polynomial_coeff_setCoefficient(inp%coefficients(ii),&
373 : 40 : & eff_pot%anharmonics_terms%coefficients(ii))
374 : : end do
375 : : end if
376 : :
377 [ - + ]: 63 : else if(filetype==4) then
378 [ # # ]: 0 : if(present(hist))then
379 : 0 : write(message,'(5a)')ch10,&
380 : 0 : & '-Reading the file ',trim(filename),ch10,&
381 : 0 : & ' with NetCDF in order to fit the polynomial coefficients'
382 : 0 : call wrtout(std_out,message,'COLL')
383 : 0 : call wrtout(ab_out,message,'COLL')
384 : 0 : call effective_potential_file_readMDfile(filename,hist,option=inp%ts_option)
385 : : else
386 : : write(message, '(3a)' )&
387 : 0 : & 'There is no hist argument ',ch10,&
388 : 0 : & 'Action: add hist argument'
389 : 0 : ABI_ERROR(message)
390 : : end if
391 : : end if
392 : : else
393 : : write(message, '(5a)' )&
394 : 0 : & ' The file ',trim(filename),' is not readable with Multibinit',ch10,&
395 : 0 : & ' Action: Change the file.'
396 : 0 : ABI_BUG(message)
397 : : end if
398 : :
399 [ - + - + : 91 : end subroutine effective_potential_file_read
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
400 : : !!***
401 : :
402 : : !!****f* m_effective_potential_file/effective_potential_file_getType
403 : : !!
404 : : !! NAME
405 : : !! effective_potential_file_getType
406 : : !!
407 : : !! FUNCTION
408 : : !! This routine test the xml or ddb file
409 : : !!
410 : : !! INPUTS
411 : : !! filename = names of the files
412 : : !!
413 : : !! OUTPUT
414 : : !! type_file = 0 no type found
415 : : !! 1 DDB file
416 : : !! 2 XML file the system definition and harmonic part
417 : : !! 3 XML file with polynomial coefficients
418 : : !! 23 XML file with both system definition and polynomial coefficients
419 : : !! 40 NetCDF file with history of MD or snapshot
420 : : !! 41 ASCII file with history of MD or snapshot
421 : : !!
422 : : !! SOURCE
423 : :
424 : 290 : subroutine effective_potential_file_getType(filename,filetype)
425 : :
426 : : !Arguments ------------------------------------
427 : : !scalars
428 : : character(len=fnlen),intent(in) :: filename
429 : : integer, intent(out) :: filetype
430 : : !arrays
431 : : !Local variables-------------------------------
432 : : !scalar
433 : : integer :: natom,nstep
434 : : integer :: ddbun = 666,ios=0
435 : : character(len=500) :: message
436 : : character (len=1000) :: line,readline
437 : : integer :: natom_id,time_id,xyz_id,six_id,ddb_version
438 : : integer :: ncid,ncerr
439 : : logical :: md_file
440 : : ! *************************************************************************
441 : :
442 : 290 : filetype = 0
443 : :
444 : 290 : ddbun = get_unit()
445 : :
446 [ - + ]: 290 : if (open_file(filename,message,unit=ddbun,form="formatted",status="old",action="read") /= 0) then
447 : 0 : ABI_ERROR(message)
448 : : end if
449 : :
450 : : !Check if the file is a XML file or a DDB and in this case, store the DDB code.
451 : 290 : ios = 0
452 [ + + ]: 26199 : do while ((ios==0))
453 : 25909 : read(ddbun,'(a)',iostat=ios) readline
454 : 25909 : call rmtabfromline(readline)
455 : 25909 : line=adjustl(readline)
456 [ + + ]: 26199 : if(line(3:13)=="xml version") then
457 [ + + ]: 380 : do while ((ios==0))
458 : 190 : read(ddbun,'(a)',iostat=ios) readline
459 : 190 : call rmtabfromline(readline)
460 : 190 : line=adjustl(readline)
461 [ + + - + ]: 190 : if(line(1:16)==char(60)//"Heff_definition".or.&
462 : : & line(1:17)==char(60)//"Terms_definition")then
463 : 72 : filetype = 3
464 : 72 : ios = -1
465 : : end if
466 [ + + ]: 380 : if(line(1:18)==char(60)//"System_definition") then
467 : 118 : filetype = 2
468 [ + + ]: 183318 : do while ((ios==0))
469 : 183200 : read(ddbun,'(a)',iostat=ios) readline
470 : 183200 : call rmtabfromline(readline)
471 : 183200 : line=adjustl(readline)
472 [ + + + - ]: 183200 : if(line(1:16)==char(60)//"Heff_definition".or.&
473 : 118 : & line(1:17)==char(60)//"Terms_definition")then
474 : 16 : filetype = 23
475 : 16 : ios = -1
476 : : end if
477 : : end do
478 : : end if
479 : : end do
480 [ + + ]: 25719 : else if(line(6:24)=="DERIVATIVE DATABASE") then
481 : 87 : filetype = 1
482 : 87 : ios = -1
483 : : end if
484 : : end do
485 : 290 : close(ddbun)
486 : :
487 [ + + ]: 303 : if(filetype/=0) return
488 : :
489 : : !try to read netcdf HIST file
490 : 13 : ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid)
491 [ + - ]: 13 : if(ncerr == NF90_NOERR) then
492 : 13 : md_file = .TRUE.
493 : 13 : ncerr = nf90_inq_dimid(ncid,"natom",natom_id)
494 [ - + ]: 13 : if(ncerr /= NF90_NOERR) md_file = .FALSE.
495 : 13 : ncerr = nf90_inq_dimid(ncid,"xyz",xyz_id)
496 [ - + ]: 13 : if(ncerr /= NF90_NOERR) md_file = .FALSE.
497 : 13 : ncerr = nf90_inq_dimid(ncid,"time",time_id)
498 [ - + ]: 13 : if(ncerr /= NF90_NOERR) md_file = .FALSE.
499 : 13 : ncerr = nf90_inq_dimid(ncid,"six",six_id)
500 [ + - ]: 13 : if(ncerr /= NF90_NOERR) md_file = .FALSE.
501 [ + - ]: 13 : if (md_file) then
502 : 13 : filetype = 40
503 : 13 : return
504 : : end if
505 : : end if
506 : 0 : ncerr = nf90_close(ncid)
507 : :
508 [ # # ]: 0 : if(filetype/=0) return
509 : :
510 : : ! Try to read netcdf DDB file
511 : 0 : ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid)
512 [ # # ]: 0 : if(ncerr==NF90_NOERR) then
513 : 0 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, 'ddb_version'), ddb_version)
514 [ # # ]: 0 : if (ncerr==NF90_NOERR) then
515 : 0 : filetype = 1
516 : 0 : return
517 : : end if
518 : : end if
519 : :
520 : : !Try to get the dim of MD ASCII file
521 : 0 : call effective_potential_file_getDimMD(filename,natom,nstep)
522 [ # # # # ]: 0 : if(natom /= 0 .and. nstep/=0) filetype = 41
523 : :
524 : 91 : end subroutine effective_potential_file_getType
525 : : !!***
526 : :
527 : : !!****f* m_effective_potential_file/effective_potential_file_getDimSystem
528 : : !!
529 : : !! NAME
530 : : !! effective_potential_file_getDimSystem
531 : : !!
532 : : !! FUNCTION
533 : : !! This routine test the xml or ddb file
534 : : !! Return the number of atoms/ntypat in the unit cell from ddb and xml
535 : : !! Return natom/ntypat/nqpt and nrpt if the file is XML file
536 : : !! In case of DDB file, you have to run bigbx9 to get nrpt
537 : : !!
538 : : !! INPUTS
539 : : !! filename = names of the files
540 : : !! comm = MPI communicator
541 : : !!
542 : : !! OUTPUT
543 : : !! natom = number of atoms
544 : : !! ntypat= number of type of atoms
545 : : !! nqpt = number of q points
546 : : !! nrpt = number of rpt points
547 : : !!
548 : : !! SOURCE
549 : :
550 : 134 : subroutine effective_potential_file_getDimSystem(filename,comm,natom,ntypat,nqpt,nrpt)
551 : :
552 : : !Arguments ------------------------------------
553 : : !scalars
554 : : character(len=fnlen),intent(in) :: filename
555 : : integer,intent(out) :: natom,ntypat,nqpt,nrpt
556 : : integer,intent(in) :: comm
557 : : !arrays
558 : :
559 : : !Local variables-------------------------------
560 : : !scalar
561 : : integer :: filetype
562 : : ! integer :: dimekb,lmnmax,mband,mtyp,msym,nblok,nkpt,usepaw
563 : : character(len=500) :: message
564 [ + + + + : 6968 : type(ddb_hdr_type) :: ddb_hdr
+ + + + +
+ + + + +
+ + + + ]
565 : : ! *************************************************************************
566 : :
567 : 134 : natom = 0
568 : 134 : ntypat= 0
569 : 134 : nqpt = 0
570 : 134 : nrpt = 0
571 : :
572 : 134 : call effective_potential_file_getType(filename,filetype)
573 : :
574 [ + + ]: 134 : if(filetype==1)then
575 : 58 : write(message, '(6a)' )ch10,' The file ',trim(filename),ch10,&
576 : 116 : & ' is DDB file (extraction of the number of atoms)',ch10
577 : 58 : call wrtout(std_out,message,'COLL')
578 : :
579 : : write(message, '(8a)' )&
580 : 58 : & ' The file ',trim(filename),ch10,' is ddb file only the number of atoms is read,',&
581 : 116 : & 'if you want to predic the number of cell (nrpt)',ch10,' use bigbx9 routines',ch10
582 : 58 : call wrtout(std_out,message,'COLL')
583 : :
584 : 58 : call ddb_hdr%open_read(filename,comm,dimonly=1)
585 : 58 : natom = ddb_hdr%natom
586 : 58 : ntypat = ddb_hdr%ntypat
587 : 58 : call ddb_hdr%free()
588 : :
589 : : ! Must read some value to initialze array (nprt for ifc)
590 : : ! call bigbx9(inp%brav,dummy_cell,0,1,inp%ngqpt,inp%nqshft,nrpt,ddb%rprim,dummy_rpt)
591 : :
592 [ + - ]: 76 : else if (filetype==2 .or. filetype==23) then
593 : 76 : write(message, '(5a)' )ch10,' The file ',trim(filename),&
594 : 152 : & ' is XML file (extraction of all information)'
595 : 76 : call wrtout(std_out,message,'COLL')
596 : :
597 : 76 : call system_getDimFromXML(filename,natom,ntypat,nqpt,nrpt)
598 : :
599 : : else
600 : : write(message, '(a,a,a,a)' )&
601 : 0 : & ' The file ',trim(filename),' is not compatible with multibinit',ch10
602 : 0 : ABI_ERROR(message)
603 : : end if
604 : :
605 : : ! TODO hexu: temporarily disabled. Discuss with alex how to do this properly.
606 : : ! Do some checks
607 : : ! if (natom < 1) then
608 : : ! write(message, '(a,a,a,a,a)' )&
609 : : !& ' Unable to read the number of atom from ',trim(filename),ch10,&
610 : : !& 'This file is not compatible with multibinit',ch10
611 : : ! ABI_ERROR(message)
612 : : ! end if
613 : : !
614 : : ! if (filetype==2 .or. filetype==23) then
615 : : !
616 : : ! if (natom < 1) then
617 : : ! write(message, '(a,a,a)' )&
618 : : !& ' Unable to read the number of atom from ',trim(filename),ch10
619 : : ! ABI_ERROR(message)
620 : : ! end if
621 : : !
622 : : ! if (nrpt < 1) then
623 : : ! write(message, '(a,a,a)' )&
624 : : !& ' Unable to read the number of rpt points ',trim(filename),ch10
625 : : ! ABI_ERROR(message)
626 : : ! end if
627 : : !
628 : : ! if (ntypat < 1) then
629 : : ! write(message, '(a,a,a)' )&
630 : : !& ' Unable to read the number of type of atoms ',trim(filename),ch10
631 : : ! ABI_ERROR(message)
632 : : ! end if
633 : : !
634 : : ! end if
635 : :
636 [ - + - + : 134 : end subroutine effective_potential_file_getDimSystem
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
637 : : !!***
638 : :
639 : : !!****f* m_effective_potential_file/effective_potential_file_getDimCoeff
640 : : !!
641 : : !! NAME
642 : : !! effective_potential_file_getDimCoeff
643 : : !!
644 : : !! FUNCTION
645 : : !! This routine test the xml with polynomial coefficients
646 : : !! Return the number of coefficients and the maximum number of displacement/strain
647 : : !!
648 : : !! INPUTS
649 : : !! filename = names of the files
650 : : !!
651 : : !! OUTPUT
652 : : !! ncoeff = number of coefficient for the polynome
653 : : !! nterm(ncoeff) = number terms per coefficient
654 : : !! ndisp(nterm,ncoeff) = number displacement per term
655 : : !!
656 : : !! SOURCE
657 : :
658 : 28 : subroutine effective_potential_file_getDimCoeff(filename,ncoeff,ndisp_max,nterm_max)
659 : :
660 : : !Arguments ------------------------------------
661 : : !scalars
662 : : character(len=fnlen),intent(in) :: filename
663 : : integer,intent(out) :: ncoeff,ndisp_max,nterm_max
664 : : !Local variables-------------------------------
665 : : !scalar
666 : : integer :: filetype
667 : : #ifndef HAVE_XML
668 : : integer :: count,count2
669 : : integer :: funit = 1,ios=0
670 : : logical :: found,found2
671 : : #endif
672 : : !arrays
673 : : #ifndef HAVE_XML
674 : : character (len=XML_RECL) :: line,readline
675 : : #endif
676 : : character(len=500) :: message
677 : : ! *************************************************************************
678 : :
679 : 28 : call effective_potential_file_getType(filename,filetype)
680 : :
681 [ + - ]: 28 : if (filetype==3 .or. filetype==23) then
682 : 28 : write(message, '(2a)' )' Extraction of the number of coefficient in the XML ',&
683 : 56 : & trim(filename)
684 : 28 : call wrtout(std_out,message,'COLL')
685 : :
686 : 28 : ncoeff = 0
687 : 28 : nterm_max = 0
688 : 28 : ndisp_max = 0
689 : :
690 : : #if defined HAVE_XML
691 : : ! Read with libxml the number of coefficient
692 : : call effpot_xml_getDimCoeff(char_f2c(trim(filename)),ncoeff,nterm_max,ndisp_max)
693 : : #else
694 : : ! Read by hand
695 : : ! Start a reading loop
696 : 28 : found=.false.
697 : : ncoeff = 0
698 : :
699 [ - + ]: 28 : if (open_file(filename,message,unit=funit,form="formatted",status="old",&
700 : : & action="read") /= 0) then
701 : 0 : ABI_ERROR(message)
702 : : end if
703 : :
704 : : ! First parse to know the number of coefficients
705 : 28 : ios = 0
706 [ + + ]: 5038 : do while (ios == 0)
707 : 5010 : read(funit,'(a)',iostat=ios) readline
708 [ + + ]: 5038 : if(ios == 0)then
709 : 4982 : call rmtabfromline(readline)
710 : 4982 : line=adjustl(readline)
711 : : ! Need test with char(9) because the old version of XML file
712 : : ! from old script includes tarbulation at the begining of each line
713 [ + + ]: 4982 : if (line(1:12)==char(60)//'coefficient') then
714 : 1802 : ncoeff=ncoeff+1
715 : 1802 : count = 0
716 : 1802 : found = .false.
717 : : do while(.not.found)
718 : 49004 : read(funit,'(a)',iostat=ios) readline
719 : 49004 : call rmtabfromline(readline)
720 : 49004 : line=adjustl(readline)
721 [ + + ]: 49004 : if (line(1:5)==char(60)//'term') then
722 : 47202 : count = count +1
723 : 47202 : found2 = .false.
724 : 47202 : count2 = 0
725 : : do while(.not.found2)
726 : 501636 : read(funit,'(a)',iostat=ios) readline
727 : 501636 : call rmtabfromline(readline)
728 : 501636 : line=adjustl(readline)
729 [ + + ]: 501636 : if (line(1:13)==char(60)//'displacement') then
730 : 109854 : count2 = count2 + 1
731 [ + + ]: 391782 : else if (line(1:7)==char(60)//'strain') then
732 : 15018 : count2 = count2 + 1
733 [ + + ]: 376764 : else if (line(1:6)==char(60)//'/term') then
734 [ + + ]: 47202 : if (count2 > ndisp_max) ndisp_max = count2
735 : : found2 = .true.
736 : : else
737 : 329562 : cycle
738 : : end if
739 : : end do
740 [ - + ]: 1802 : else if (line(1:13)==char(60)//'/coefficient') then
741 [ + + ]: 1802 : if (count > nterm_max) nterm_max = count
742 : : found = .true.
743 : : else
744 : 0 : cycle
745 : : end if
746 : : end do
747 : : cycle
748 : : end if
749 : : end if
750 : : end do
751 : :
752 : 28 : close(funit)
753 : : #endif
754 : :
755 : : else
756 : : ! Maybe one day add an other type of file...
757 : : write(message, '(a,a,a,a)' )&
758 : 0 : & ' The file ',trim(filename),' is not compatible with multibinit',ch10
759 : 0 : ABI_ERROR(message)
760 : : end if
761 : :
762 : : ! Do some checks
763 [ - + ]: 28 : if (ncoeff < 1) then
764 : : write(message, '(5a)' )&
765 : 0 : & ' Unable to read the number of coeff from ',trim(filename),ch10,&
766 : 0 : & ' This file is not compatible with multibinit',ch10
767 : 0 : ABI_ERROR(message)
768 : : end if
769 : :
770 : 28 : end subroutine effective_potential_file_getDimCoeff
771 : : !!***
772 : :
773 : :
774 : : !!****f* m_effective_potential_file/effective_potential_file_getDimStrainCoupling
775 : : !!
776 : : !! NAME
777 : : !! effective_potential_file_getDimStrainCoupling
778 : : !!
779 : : !! FUNCTION
780 : : !! Return the number of nrpt for specific strain coupling from xml system file
781 : : !!
782 : : !! INPUTS
783 : : !! filename = names of the files
784 : : !! voigt = voigt notation of the strain
785 : : !!
786 : : !! OUTPUT
787 : : !! nrpt = number of rpt points
788 : : !!
789 : : !! SOURCE
790 : :
791 : 228 : subroutine effective_potential_file_getDimStrainCoupling(filename,nrpt,voigt)
792 : :
793 : : !Arguments ------------------------------------
794 : : !scalars
795 : : character(len=fnlen),intent(in) :: filename
796 : : integer,intent(in) :: voigt
797 : : integer,intent(out) :: nrpt
798 : : !Local variables-------------------------------
799 : : !scalar
800 : : #ifndef HAVE_XML
801 : : integer :: irpt,ivoigt
802 : : integer :: funit = 1,ios=0
803 : : logical :: found
804 : : #endif
805 : : !arrays
806 : : #ifndef HAVE_XML
807 : : character (len=XML_RECL) :: line,readline,strg,strg1
808 : : character(len=500) :: message
809 : : #endif
810 : : ! *************************************************************************
811 : :
812 : 228 : nrpt = 0
813 : :
814 : : #if defined HAVE_XML
815 : : ! Read with libxml the number of coefficient
816 : : call effpot_xml_getDimStrainCoupling(char_f2c(trim(filename)),nrpt,voigt)
817 : : #else
818 : : ! Read by hand
819 : : ! Start a reading loop
820 : 228 : found=.false.
821 : :
822 [ - + ]: 228 : if (open_file(filename,message,unit=funit,form="formatted",status="old",&
823 : : & action="read") /= 0) then
824 : 0 : ABI_ERROR(message)
825 : : end if
826 : :
827 : : ! First parse to know the number of atoms
828 [ + - + + ]: 349458 : do while (ios == 0.and.(.not.found))
829 : 349230 : read(funit,'(a)',iostat=ios) readline
830 [ - + ]: 349458 : if(ios == 0)then
831 : 349230 : call rmtabfromline(readline)
832 : 349230 : line=adjustl(readline)
833 [ + + ]: 349230 : if ((line(1:16)=='<strain_coupling')) then
834 : 798 : read(funit,'(a)',iostat=ios) readline
835 : 798 : call rdfromline("voigt",line,strg)
836 : 798 : strg1=trim(strg)
837 : 798 : read(strg1,*) ivoigt
838 [ + + ]: 798 : if (ivoigt == voigt)then
839 : : irpt = 0
840 : : do while (.not.found)
841 : 2592 : read(funit,'(a)',iostat=ios) readline
842 : 2592 : call rmtabfromline(readline)
843 : 2592 : line=adjustl(readline)
844 [ - + ]: 2592 : if ((line(1:26)=='<correction_force_constant')) then
845 : 0 : irpt = irpt + 1
846 : 0 : cycle
847 : : end if
848 [ + + ]: 2592 : if ((line(1:17)=='</strain_coupling')) then
849 : 228 : found = .TRUE.
850 : 228 : nrpt = irpt
851 : : cycle
852 : : end if
853 : : end do
854 : : else
855 : : cycle
856 : : end if
857 : : end if
858 : : end if
859 : : end do
860 : :
861 : 228 : close(funit)
862 : : #endif
863 : :
864 : 228 : end subroutine effective_potential_file_getDimStrainCoupling
865 : : !!***
866 : :
867 : : !!****f* m_effective_potential_file/effective_potential_file_getDimMD
868 : : !!
869 : : !! NAME
870 : : !! effective_potential_file_getDimMD
871 : : !!
872 : : !! FUNCTION
873 : : !! Read MD FILE (HIST or ASCII) and return the dimensions
874 : : !! (natom and nstep)
875 : : !!
876 : : !! INPUTS
877 : : !! filename = path of the file
878 : : !!
879 : : !! OUTPUT
880 : : !! natom = number of atoms
881 : : !! nstep = number of MD steps
882 : : !!
883 : : !! SOURCE
884 : :
885 : 0 : subroutine effective_potential_file_getDimMD(filename,natom,nstep)
886 : :
887 : : !Arguments ------------------------------------
888 : : !scalars
889 : : integer,intent(out) :: natom,nstep
890 : : !arrays
891 : : character(len=fnlen),intent(in) :: filename
892 : : !Local variables-------------------------------
893 : : !scalar
894 : : integer :: ia,natm_old,natm_new
895 : : integer :: nenergy,nrprimd
896 : : integer :: ios=0,ios2=0,ios3=0
897 : : integer :: unit_md=24
898 : : logical :: compatible,netcdf
899 : : integer :: natom_id,time_id,xyz_id,six_id
900 : : integer :: ncid,ncerr
901 : : character(len=5) :: char_tmp
902 : : !arrays
903 : : character (len=10000) :: readline,line
904 : : character(len=500) :: msg
905 : : ! *************************************************************************
906 : :
907 : 0 : natom = 0
908 : 0 : nstep = 0
909 : : !try to read netcdf
910 : 0 : netcdf = .false.
911 : 0 : ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid)
912 [ # # ]: 0 : if(ncerr == NF90_NOERR) then
913 : 0 : netcdf = .TRUE.
914 : 0 : ncerr = nf90_inq_dimid(ncid,"natom",natom_id)
915 [ # # ]: 0 : if(ncerr /= NF90_NOERR) netcdf = .FALSE.
916 : 0 : ncerr = nf90_inq_dimid(ncid,"xyz",xyz_id)
917 [ # # ]: 0 : if(ncerr /= NF90_NOERR) netcdf = .FALSE.
918 : 0 : ncerr = nf90_inq_dimid(ncid,"time",time_id)
919 [ # # ]: 0 : if(ncerr /= NF90_NOERR) netcdf = .FALSE.
920 : 0 : ncerr = nf90_inq_dimid(ncid,"six",six_id)
921 [ # # ]: 0 : if(ncerr /= NF90_NOERR) netcdf = .FALSE.
922 [ # # ]: 0 : if(netcdf)then
923 : 0 : ncerr = nf90_inquire_dimension(ncid,natom_id,char_tmp,natom)
924 [ # # ]: 0 : NCF_CHECK_MSG(ncerr," inquire dimension ID for natom")
925 : 0 : ncerr = nf90_inquire_dimension(ncid,time_id,char_tmp,nstep)
926 [ # # ]: 0 : NCF_CHECK_MSG(ncerr," inquire dimension ID for time")
927 : : end if
928 : : end if
929 : :
930 [ # # ]: 0 : if(.not.netcdf) then
931 : : ! try to read ASCII file...
932 [ # # ]: 0 : if (open_file(filename,msg,unit=unit_md,form="formatted",&
933 : : & status="old",action="read") /= 0) then
934 : 0 : ABI_ERROR(msg)
935 : : end if
936 : :
937 : : ! Start a reading loop in fortran to get the dimension of the file
938 : 0 : rewind(unit=unit_md)
939 : 0 : ios = 0
940 : 0 : nstep = 0
941 : 0 : nrprimd = 0
942 : 0 : natm_old= 0
943 : 0 : natm_new= 0
944 : 0 : nenergy = 0
945 : 0 : compatible = .TRUE.
946 : :
947 [ # # ]: 0 : do while ((ios==0))
948 : : ! special treatment of the first step
949 [ # # ]: 0 : if(nstep==0)then
950 : 0 : ios2 = 0
951 [ # # ]: 0 : do while ((ios2==0))
952 : 0 : read(unit_md,'(a)',iostat=ios) readline
953 : 0 : line=adjustl(readline)
954 : 0 : call elementfromline(line,ia)
955 [ # # ]: 0 : if (ia==1)then
956 : 0 : nstep = nstep + 1
957 : 0 : ios2 = 1
958 : : end if
959 : : end do
960 : : end if
961 : 0 : read(unit_md,'(a)',iostat=ios) readline
962 [ # # ]: 0 : if(ios == 0)then
963 : 0 : line=adjustl(readline)
964 : 0 : call elementfromline(line,ia)
965 [ # # ]: 0 : if (ia==1)then
966 : 0 : nenergy = nenergy + 1
967 : 0 : nrprimd = 0
968 [ # # ]: 0 : else if(ia==3)then
969 : 0 : nrprimd = nrprimd + 1
970 : : end if
971 [ # # ]: 0 : if(nrprimd == 3)then
972 : 0 : ios3 = 0
973 : 0 : natm_new = 0
974 [ # # ]: 0 : do while ((ios3==0))
975 : 0 : read(unit_md,'(a)',iostat=ios3) readline
976 [ # # ]: 0 : if(ios3==0)then
977 : 0 : line=adjustl(readline)
978 : 0 : call elementfromline(line,ia)
979 [ # # ]: 0 : if(ia==1)then
980 [ # # ]: 0 : if(nstep==1) then
981 : : natm_old = natm_new
982 : : else
983 [ # # ]: 0 : if(natm_old /= natm_new) compatible = .FALSE.
984 : : end if
985 : 0 : ios3 = 1
986 : 0 : ios2 = 1
987 : 0 : nstep = nstep + 1
988 : : end if
989 [ # # ]: 0 : if(ia==6)then
990 : 0 : natm_new = natm_new + 1
991 : : end if
992 : : end if!end if ios3
993 : : end do
994 : : end if ! end if nrprimd
995 : : end if! end if os1
996 : : end do
997 : :
998 : 0 : natom = natm_new - 1
999 [ # # ]: 0 : if(nstep /= nenergy) compatible = .FALSE.
1000 [ # # ]: 0 : if(natom <= 0) compatible = .FALSE.
1001 [ # # ]: 0 : if(nstep <= 0) compatible = .FALSE.
1002 : :
1003 [ # # ]: 0 : if(.not.compatible)then
1004 : 0 : natom = 0
1005 : 0 : nstep = 0
1006 : : end if
1007 : : end if! end if not netcdf
1008 : :
1009 : 0 : end subroutine effective_potential_file_getDimMD
1010 : : !!***
1011 : :
1012 : : !!****f* m_effective_potential_file/system_getDimFromXML
1013 : : !! NAME
1014 : : !! system_getDimFromXML
1015 : : !!
1016 : : !! FUNCTION
1017 : : !! Open xml file of effective potentiel, then reads the variables that
1018 : : !! must be known in order to dimension the arrays before complete reading
1019 : : !!
1020 : : !! INPUTS
1021 : : !! character(len=*) filnam: name of input or output file
1022 : : !!
1023 : : !! OUTPUT
1024 : : !! natom=number of atoms
1025 : : !! ntypat=number of atom types
1026 : : !! nrpt =number of real space points used to integrate IFC
1027 : : ! nph1l =number of wavevectors for phonon
1028 : : !!
1029 : : !! SOURCE
1030 : :
1031 : 76 : subroutine system_getDimFromXML(filename,natom,ntypat,nph1l,nrpt)
1032 : :
1033 : : !Arguments ------------------------------------
1034 : : !scalars
1035 : : character(len=fnlen),intent(in) :: filename
1036 : : integer, intent(out) :: natom,ntypat,nph1l,nrpt
1037 : : !arrays
1038 : : !Local variables-------------------------------
1039 : : !scalar
1040 : : integer :: nrpt1,nrpt2
1041 : : real :: itypat
1042 : : character(len=500) :: message
1043 : : #ifndef HAVE_XML
1044 : : integer :: funit = 1,ios = 0
1045 : : integer :: iatom
1046 : : logical :: found
1047 : : character (len=XML_RECL) :: line,readline
1048 : : character (len=XML_RECL) :: strg,strg1
1049 : : #endif
1050 : : !arrays
1051 : : #ifndef HAVE_XML
1052 : 76 : integer,allocatable :: typat(:)
1053 : : #endif
1054 : : ! *************************************************************************
1055 : :
1056 : : !Open the atomicdata XML file for reading
1057 : 76 : write(message,'(5a)') ' system_getDimFromXML :',&
1058 : 76 : & '-Opening the file ',trim(filename),' to read dimensions',&
1059 : 152 : & ' (before initialisation)'
1060 : :
1061 : 76 : call wrtout(std_out,message,'COLL')
1062 : :
1063 : 76 : natom = 0
1064 : 76 : ntypat= 0
1065 : 76 : nph1l = 0
1066 : 76 : nrpt = 0
1067 : 76 : nrpt1 = 0
1068 : 76 : nrpt2 = 0
1069 : 76 : itypat= 0
1070 : :
1071 : : !Open the atomicdata XML file for reading
1072 : :
1073 : : #if defined HAVE_XML
1074 : : !Read with libxml
1075 : : call effpot_xml_getDimSystem(char_f2c(trim(filename)),natom,ntypat,nph1l,nrpt1,nrpt2)
1076 : : #else
1077 : : !Read by hand
1078 : :
1079 : : !Start a reading loop
1080 : 76 : found=.false.
1081 : :
1082 [ - + ]: 76 : if (open_file(filename,message,unit=funit,form="formatted",status="old",&
1083 : : & action="read") /= 0) then
1084 : 0 : ABI_ERROR(message)
1085 : : end if
1086 : :
1087 : : !First parse to know the number of atoms
1088 [ + + ]: 297112 : do while ((ios==0).and.(.not.found))
1089 : 297036 : read(funit,'(a)',iostat=ios) readline
1090 [ + + ]: 297112 : if(ios ==0)then
1091 : 296960 : call rmtabfromline(readline)
1092 : 296960 : line=adjustl(readline)
1093 : :
1094 : : !Need test with char(9) because the old version of XML file
1095 : : !from old script includes tarbulation at the begining of each line
1096 [ + + ]: 296960 : if (line(1:5)==char(60)//'atom') then
1097 : 408 : natom=natom+1
1098 : 408 : cycle
1099 : : end if
1100 : :
1101 [ + + ]: 296552 : if (line(1:21)==char(60)//'local_force_constant') then
1102 : 2256 : nrpt1 = nrpt1+1
1103 : 2256 : cycle
1104 : : end if
1105 : :
1106 [ + + ]: 294296 : if (line(1:21)==char(60)//'total_force_constant') then
1107 : 2166 : nrpt2 = nrpt2+1
1108 : 2166 : cycle
1109 : : end if
1110 : :
1111 [ + + ]: 292130 : if (line(1:7)==char(60)//'qpoint') then
1112 : 284 : nph1l = nph1l+1
1113 : 284 : cycle
1114 : : end if
1115 : : end if
1116 : : end do
1117 : :
1118 : : !second parse to get the number of typat
1119 [ - + + - : 228 : ABI_MALLOC(typat,(natom))
+ - + - -
+ - + ]
1120 [ + + ]: 484 : typat = 0
1121 : 76 : iatom = 0
1122 : :
1123 : 76 : rewind(funit)
1124 : : !Start a reading loop
1125 : 76 : ios = 0
1126 : 76 : found = .false.
1127 : :
1128 [ + - + + ]: 5980 : do while ((ios==0).and.(.not.found))
1129 : 5904 : read(funit,'(a)',iostat=ios) readline
1130 [ - + ]: 5980 : if(ios == 0)then
1131 : 5904 : call rmtabfromline(readline)
1132 : 5904 : line=adjustl(readline)
1133 : :
1134 [ + + ]: 5904 : if (line(1:5)==char(60)//'atom') then
1135 : 408 : iatom = iatom + 1
1136 : 408 : call rdfromline("mass",line,strg)
1137 : 408 : strg1=trim(strg)
1138 : 408 : read(unit=strg1,fmt=*) itypat
1139 [ + + + + : 2128 : if (.not.any(typat==int(itypat))) then
+ + ]
1140 : 196 : ntypat= ntypat+1
1141 : : end if
1142 : 408 : typat(iatom) = int(itypat)
1143 : : end if
1144 : :
1145 [ + + ]: 5904 : if (line(1:6)==char(60)//'local') then
1146 : 76 : found=.true.
1147 : : end if
1148 : : end if
1149 : : end do
1150 : :
1151 : 76 : close(funit)
1152 [ - + ]: 76 : ABI_FREE(typat)
1153 : :
1154 : : #endif
1155 : :
1156 : : !Check the RPT
1157 [ + + ]: 76 : if (nrpt2/=nrpt1) then
1158 [ + + ]: 46 : if(nrpt1> 0 .and. nrpt2== 0) then
1159 : : continue;
1160 [ - + ]: 32 : else if (nrpt1==0.and.nrpt2>=0) then
1161 : 0 : write(message, '(5a)' )ch10,&
1162 : 0 : & ' WARNING: the number of local IFC is set to 0 ',ch10,&
1163 : 0 : & ' Dipdip must be set to zero',ch10
1164 : 0 : call wrtout(std_out,message,'COLL')
1165 [ + - ]: 32 : else if (nrpt2 > nrpt1) then
1166 : 32 : write(message, '(2a,I0,3a,I0,5a)' )ch10,&
1167 : 32 : & ' WARNING: the number of total IFC (',nrpt2,') is not equal to the ',ch10,&
1168 : 32 : & ' the number of short range IFC (',nrpt1,') in ',trim(filename),ch10,&
1169 : 64 : & ' the missing ifc will be set to zero',ch10
1170 : 32 : call wrtout(std_out,message,'COLL')
1171 : : else if(nrpt1>nrpt2)then
1172 : 0 : write(message, '(2a,I0,3a,I0,5a)' )ch10,&
1173 : 0 : & ' The number of total IFC (',nrpt2,') is inferior to ',ch10,&
1174 : 0 : & ' the number of short range IFC (',nrpt1,') in ',trim(filename),ch10,&
1175 : 0 : & ' This is not possible',ch10
1176 : 0 : ABI_BUG(message)
1177 : : end if
1178 : : end if
1179 : :
1180 : : !nrpt is the max between local and total:
1181 : 76 : nrpt = max(nrpt1,nrpt2)
1182 : :
1183 : :
1184 : 76 : end subroutine system_getDimFromXML
1185 : : !!***
1186 : :
1187 : : !!****f* m_effective_potential_file/system_xml2effpot
1188 : : !! NAME
1189 : : !! system_xml2effpot
1190 : : !!
1191 : : !! FUNCTION
1192 : : !! Open xml file of effective potentiel, then reads the variables
1193 : : !! and store them in effective potentential type
1194 : : !!
1195 : : !! INPUTS
1196 : : !! eff_pot<type(effective_potential_type)> = datatype with all the information for effective potential
1197 : : !! comm=MPI communicator
1198 : : !! character(len=*) filnam: name of input or output file
1199 : : !! strcpling = optional,logical to disable the strcpling
1200 : : !!
1201 : : !! OUTPUT
1202 : : !! eff_pot<type(effective_potential_type)> = datatype with all the information for effective potential
1203 : : !!
1204 : : !! SOURCE
1205 : :
1206 : 38 : subroutine system_xml2effpot(eff_pot,filename,comm,strcpling)
1207 : :
1208 : : use m_atomdata
1209 : : use m_multibinit_dataset, only : multibinit_dtset_type
1210 : : use m_ab7_symmetry
1211 : :
1212 : : !Arguments ------------------------------------
1213 : : !scalars
1214 : : character(len=*),intent(in) :: filename
1215 : : integer, intent(in) :: comm
1216 : : integer, optional,intent(in) :: strcpling
1217 : : !arrays
1218 : : type(effective_potential_type), intent(inout) :: eff_pot
1219 : :
1220 : : !Local variables-------------------------------
1221 : : !scalar
1222 : : integer :: ierr,ii,itypat,my_rank,msym,natom,ncoeff,nrpt,nrpt_scoupling
1223 : : integer :: ntypat,nph1l,nptsym,npsp,nproc,nsym,space_group,timrev,use_inversion,voigt
1224 : : real(dp):: energy,tolsym,ucvol
1225 : : character(len=500) :: message
1226 : : integer,parameter :: master=0
1227 : : logical :: has_anharmonics
1228 : : logical :: iam_master
1229 : : #ifndef HAVE_XML
1230 : : integer :: funit = 1,ios=0
1231 : : integer :: iatom,iamu,iph1l,irpt,irpt1,irpt2,irpt3,jj,mu,nu
1232 : : real(dp):: amu
1233 : : logical :: found,found2,short_range,total_range
1234 : : character (len=XML_RECL) :: line,readline
1235 : : character (len=XML_RECL) :: strg,strg1
1236 : : logical :: has_straincoupling
1237 : : #endif
1238 : : !arrays
1239 : : integer :: bravais(11)
1240 : 38 : integer,allocatable :: typat(:)
1241 : 38 : integer,allocatable :: symrel(:,:,:),symafm(:),ptsymrel(:,:,:)
1242 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1243 : : real(dp) :: elastic_constants(6,6),elastic3rd(6,6,6),epsilon_inf(3,3)
1244 : 38 : real(dp),allocatable :: all_amu(:),cell_local(:,:),cell_total(:,:)
1245 : 38 : real(dp),allocatable :: elastic_displacement(:,:,:,:),dynmat(:,:,:,:,:,:)
1246 : 38 : real(dp),allocatable :: local_atmfrc(:,:,:,:,:),total_atmfrc(:,:,:,:,:)
1247 : 38 : real(dp),allocatable :: spinat(:,:),strain_coupling(:,:,:),phfrq(:,:),qph1l(:,:),tnons(:,:)
1248 : 38 : real(dp),allocatable :: xcart(:,:),xred(:,:),zeff(:,:,:),znucl(:),zion(:)
1249 : 38 : character(len=132),allocatable :: title(:)
1250 [ + + + + ]: 228 : type(ifc_type) :: ifcs
1251 : 38 : type(ifc_type),dimension(:),allocatable :: phonon_strain
1252 [ + + + + : 1976 : type(crystal_t) :: crystal
+ + + + +
+ + + + +
+ + + + ]
1253 : : type(atomdata_t) :: atom
1254 : : #ifdef HAVE_XML
1255 : : real(dp),allocatable :: phonon_strain_atmfrc(:,:,:,:,:)
1256 : : integer,allocatable :: phonon_straincell(:,:)
1257 : : #endif
1258 : : #ifndef HAVE_XML
1259 : 38 : real(dp),allocatable :: work2(:,:)
1260 : : #endif
1261 : : ! *************************************************************************
1262 : :
1263 : : !Open the atomicdata XML file for reading
1264 : 38 : write(message,'(a,a)')'-Opening the file ',trim(filename)
1265 : :
1266 : 38 : call wrtout(ab_out,message,'COLL')
1267 : 38 : call wrtout(std_out,message,'COLL')
1268 : :
1269 : 38 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1270 : 38 : iam_master = (my_rank == master)
1271 : :
1272 : : !Get Dimention of system and allocation/initialisation of array
1273 : 38 : call effective_potential_file_getDimSystem(filename,comm,natom,ntypat,nph1l,nrpt)
1274 : 38 : gmet= zero; gprimd = zero; rmet = zero; rprimd = zero
1275 : 38 : elastic_constants = zero; epsilon_inf = zero; ncoeff = 0
1276 [ - + + - : 114 : ABI_MALLOC(all_amu,(ntypat))
+ - + - -
+ - + ]
1277 [ - + + - : 114 : ABI_MALLOC(cell_local,(3,nrpt))
+ - + - -
+ - + ]
1278 [ - + + - : 76 : ABI_MALLOC(cell_total,(3,nrpt))
- + - + ]
1279 [ - + + - : 114 : ABI_MALLOC(elastic_displacement,(6,6,3,natom))
+ - + - -
+ - + ]
1280 [ - + + - : 228 : ABI_MALLOC(ifcs%atmfrc,(3,natom,3,natom,nrpt))
+ - - + +
- - + + -
+ - + - -
+ - + ]
1281 [ - + + - : 114 : ABI_MALLOC(ifcs%cell,(3,nrpt))
+ - + - -
+ - + ]
1282 [ - + + - : 152 : ABI_MALLOC(ifcs%short_atmfrc,(3,natom,3,natom,nrpt))
- + + - -
+ + - - +
- + ]
1283 [ - + + - : 152 : ABI_MALLOC(ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt))
- + + - -
+ + - - +
- + ]
1284 [ - + + - : 114 : ABI_MALLOC(strain_coupling,(6,3,natom))
+ - + - -
+ - + ]
1285 [ - + + - : 152 : ABI_MALLOC(total_atmfrc,(3,natom,3,natom,nrpt))
- + + - -
+ + - - +
- + ]
1286 [ - + + - : 152 : ABI_MALLOC(local_atmfrc,(3,natom,3,natom,nrpt))
- + + - -
+ + - - +
- + ]
1287 [ - + + - : 228 : ABI_MALLOC(dynmat,(2,3,natom,3,natom,nph1l))
+ - - + +
- + + + -
+ - + + -
+ - + ]
1288 [ - + + - : 114 : ABI_MALLOC(typat,(natom))
+ - + - -
+ - + ]
1289 [ + - - + : 152 : ABI_MALLOC(phfrq,(3*natom,nph1l))
+ - + + +
- + - + +
- + - + ]
1290 [ + + + - : 114 : ABI_MALLOC(qph1l,(3,nph1l))
+ - + + -
+ - + ]
1291 [ - + + - : 114 : ABI_MALLOC(xcart,(3,natom))
+ - + - -
+ - + ]
1292 [ - + + - : 76 : ABI_MALLOC(xred,(3,natom))
- + - + ]
1293 [ - + + - : 114 : ABI_MALLOC(zeff,(3,3,natom))
+ - + - -
+ - + ]
1294 [ - + + - : 76 : ABI_MALLOC(zion,(ntypat))
- + - + ]
1295 [ - + + - : 76 : ABI_MALLOC(znucl,(ntypat))
- + - + ]
1296 : :
1297 [ - + + + : 456 : ABI_MALLOC(phonon_strain,(6))
+ + + + ]
1298 : 38 : nrpt_scoupling = 0
1299 [ + + ]: 266 : do ii = 1,6
1300 : : ! Get The size of the strainPhonon-coupling
1301 : 228 : call effective_potential_file_getDimStrainCoupling(filename,nrpt_scoupling,ii-1)
1302 [ - + + - : 1140 : ABI_MALLOC(phonon_strain(ii)%atmfrc,(3,natom,3,natom,nrpt_scoupling))
- + + - +
- - - + -
- + - + -
+ - + ]
1303 [ + - - - : 684 : ABI_MALLOC(phonon_strain(ii)%cell,(3,nrpt_scoupling))
+ - - + -
+ - + -
+ ]
1304 : 228 : phonon_strain(ii)%nrpt = nrpt_scoupling
1305 [ - + - - : 228 : phonon_strain(ii)%atmfrc = zero
- - - - -
- ]
1306 [ - + - - ]: 266 : phonon_strain(ii)%cell = 0
1307 : : end do
1308 : :
1309 [ + + ]: 136 : all_amu(:) = zero
1310 [ + + + + : 97068 : dynmat(:,:,:,:,:,:) = zero
+ + + + +
+ + + ]
1311 [ + + + + ]: 5126 : cell_local(:,:) = 99D99
1312 [ + + + + ]: 5126 : cell_total(:,:) = 99D99
1313 : 38 : elastic3rd(:,:,:) = zero
1314 [ + + + + : 26558 : elastic_displacement(:,:,:,:) = zero
+ + + + ]
1315 : 38 : ifcs%nrpt = nrpt
1316 [ + + + + : 633086 : ifcs%atmfrc(:,:,:,:,:) = zero
+ + + + +
+ ]
1317 [ + + + + ]: 5126 : ifcs%cell(:,:) = 0
1318 [ + + + + : 633086 : ifcs%ewald_atmfrc(:,:,:,:,:) = zero
+ + + + +
+ ]
1319 [ + + + + : 633086 : ifcs%short_atmfrc(:,:,:,:,:) = zero
+ + + + +
+ ]
1320 [ + + + + : 4526 : strain_coupling(:,:,:) = zero
+ + ]
1321 [ + + + + ]: 2016 : phfrq = zero
1322 [ + + + + ]: 606 : qph1l = 0
1323 [ + + + + ]: 854 : xcart = zero
1324 [ + + + + : 2690 : zeff = zero
+ + ]
1325 [ + + ]: 136 : znucl = zero
1326 : :
1327 [ + + ]: 38 : if(iam_master)then
1328 : : !Open the atomicdata XML file for reading
1329 : : #if defined HAVE_XML
1330 : :
1331 : : write(message,'(a,a,a,a)')'-Reading the file ',trim(filename),&
1332 : : & ' with LibXML library'
1333 : :
1334 : : call wrtout(ab_out,message,'COLL')
1335 : : call wrtout(std_out,message,'COLL')
1336 : :
1337 : : ! Read with libxml library
1338 : : call effpot_xml_readSystem(char_f2c(trim(filename)),natom,ntypat,nrpt,nph1l,all_amu,&
1339 : : & ifcs%atmfrc,ifcs%cell,dynmat,elastic_constants,energy,&
1340 : : & epsilon_inf,ifcs%ewald_atmfrc,phfrq,rprimd,qph1l,&
1341 : : & ifcs%short_atmfrc,typat,xcart,zeff)
1342 : :
1343 : : ! convert atomic mass unit to znucl
1344 : : do itypat=1,ntypat
1345 : : do ii=1,103
1346 : : call atomdata_from_znucl(atom,real(ii,dp))
1347 : : if (abs((real(atom%amu,sp)-real(all_amu(itypat),sp))&
1348 : : & /real(all_amu(itypat),sp)*100)<0.1) then
1349 : : znucl(itypat) = atom%znucl
1350 : : exit
1351 : : end if
1352 : : end do
1353 : : end do
1354 : :
1355 : : ! Get the Phonon Strain coupling
1356 : : do voigt = 1,6
1357 : : nrpt_scoupling = phonon_strain(voigt)%nrpt
1358 : : ABI_MALLOC(phonon_straincell,(3,nrpt_scoupling))
1359 : : ABI_MALLOC(phonon_strain_atmfrc,(3,natom,3,natom,nrpt_scoupling))
1360 : :
1361 : : ! Get The value
1362 : : call effpot_xml_readStrainCoupling(char_f2c(trim(filename)),natom,nrpt_scoupling,(voigt-1),&
1363 : : & elastic3rd(voigt,:,:),elastic_displacement(voigt,:,:,:),&
1364 : : & strain_coupling(voigt,:,:),&
1365 : : & phonon_strain_atmfrc,phonon_straincell)
1366 : :
1367 : : ! Check if the 3rd order strain_coupling is present
1368 : : has_anharmonics = .FALSE.
1369 : : if(any(elastic3rd>tol10).or.any(elastic_displacement>tol10)) has_anharmonics = .TRUE.
1370 : : phonon_strain(voigt)%atmfrc(:,:,:,:,:) = phonon_strain_atmfrc(:,:,:,:,:)
1371 : : phonon_strain(voigt)%cell(:,:) = phonon_straincell(:,:)
1372 : : if(any(phonon_strain(voigt)%atmfrc > tol10)) has_anharmonics = .TRUE.
1373 : :
1374 : : ABI_FREE(phonon_straincell)
1375 : : ABI_FREE(phonon_strain_atmfrc)
1376 : : end do
1377 : : #else
1378 : :
1379 : : ! Read by hand
1380 : 22 : write(message,'(a,a,a,a)')'-Reading the file ',trim(filename),&
1381 : 44 : & ' with Fortran'
1382 : :
1383 : 22 : call wrtout(ab_out,message,'COLL')
1384 : 22 : call wrtout(std_out,message,'COLL')
1385 : :
1386 [ - + ]: 22 : if (open_file(filename,message,unit=funit,form="formatted",&
1387 : : & status="old",action="read") /= 0) then
1388 : 0 : ABI_ERROR(message)
1389 : : end if
1390 : :
1391 : : !Start a reading loop in fortran
1392 : 22 : rewind(unit=funit)
1393 : 22 : found=.false.
1394 : :
1395 : 22 : iatom = 1
1396 : 22 : iamu = 1
1397 : 22 : itypat = 1
1398 : 22 : irpt = 1
1399 : 22 : irpt1 = 0
1400 : 22 : irpt2 = 0
1401 : 22 : iph1l = 1
1402 : 22 : amu = zero
1403 : 22 : short_range = .false.
1404 : 22 : total_range = .false.
1405 : 22 : has_straincoupling = .FALSE.
1406 [ + + ]: 94533 : do while ((ios==0).and.(.not.found))
1407 : 94511 : read(funit,'(a)',iostat=ios) readline
1408 [ + + ]: 94533 : if(ios == 0)then
1409 : 94489 : call rmtabfromline(readline)
1410 : 94489 : line=adjustl(readline)
1411 [ + + ]: 94489 : if (.not.has_straincoupling) then
1412 : :
1413 [ + + ]: 5163 : if ((line(1:7)=='<energy')) then
1414 : 22 : call rdfromline_value('energy',line,strg)
1415 [ - + ]: 22 : if (strg/="") then
1416 : 0 : strg1=trim(strg)
1417 : 0 : read(strg1,*) energy
1418 : : else
1419 : 22 : read(funit,'(a)',iostat=ios) readline
1420 : 22 : call rmtabfromline(readline)
1421 : 22 : line=adjustl(readline)
1422 : 22 : call rdfromline_value('energy',line,strg)
1423 [ - + ]: 22 : if (strg/="") then
1424 : 0 : strg1=trim(strg)
1425 : : else
1426 : 22 : strg1=trim(line)
1427 : : end if
1428 : 22 : read(strg1,*) energy
1429 : : end if
1430 : : cycle
1431 : : end if
1432 : :
1433 [ + + ]: 5141 : if ((line(1:10)=='<unit_cell')) then
1434 : 22 : call rdfromline_value('unit_cell',line,strg)
1435 [ - + ]: 22 : if (strg/="") then
1436 : 0 : strg1=trim(strg)
1437 : 0 : read(strg1,*) (rprimd(1,mu),mu=1,3)
1438 : 0 : read(funit,*) (rprimd(2,mu),mu=1,3)
1439 : : else
1440 [ + + ]: 66 : do nu=1,2
1441 : 66 : read(funit,*) (rprimd(nu,mu),mu=1,3)
1442 : : end do
1443 : : end if
1444 : 22 : read(funit,'(a)',iostat=ios) readline
1445 : 22 : call rmtabfromline(readline)
1446 : 22 : line=adjustl(readline)
1447 : 22 : call rdfromline_value('unit_cell',line,strg)
1448 [ - + ]: 22 : if (strg/="") then
1449 : 0 : strg1=trim(strg)
1450 : : else
1451 : 22 : strg1=trim(line)
1452 : : end if
1453 : 22 : read(strg1,*) (rprimd(3,mu),mu=1,3)
1454 : 22 : cycle
1455 : : end if
1456 : :
1457 [ + + ]: 5119 : if ((line(1:12)=='<epsilon_inf')) then
1458 : 22 : call rdfromline_value('epsilon_inf',line,strg)
1459 [ - + ]: 22 : if (strg/="") then
1460 : 0 : strg1=trim(strg)
1461 : 0 : read(strg1,*) (epsilon_inf(mu,1),mu=1,3)
1462 : 0 : read(funit,*) (epsilon_inf(mu,2),mu=1,3)
1463 : : else
1464 [ + + ]: 66 : do nu=1,2
1465 : 66 : read(funit,*) (epsilon_inf(mu,nu),mu=1,3)
1466 : : end do
1467 : : end if
1468 : 22 : read(funit,'(a)',iostat=ios) readline
1469 : 22 : call rmtabfromline(readline)
1470 : 22 : line=adjustl(readline)
1471 : 22 : call rdfromline_value('epsilon_inf',line,strg)
1472 [ - + ]: 22 : if (strg/="") then
1473 : 0 : strg1=trim(strg)
1474 : : else
1475 : 22 : strg1=trim(line)
1476 : : end if
1477 : 22 : read(strg1,*) (epsilon_inf(mu,3),mu=1,3)
1478 : 22 : cycle
1479 : : end if
1480 : :
1481 [ + + ]: 5097 : if ((line(1:8)=='<elastic')) then
1482 : 22 : call rdfromline_value('elastic',line,strg)
1483 [ - + ]: 22 : if (strg/="") then
1484 : 0 : strg1=trim(strg)
1485 : 0 : read(strg1,*) (elastic_constants(mu,1),mu=1,6)
1486 [ # # ]: 0 : do nu=2,5
1487 : 0 : read(funit,*) (elastic_constants(mu,nu),mu=1,6)
1488 : : end do
1489 : : else
1490 [ + + ]: 132 : do nu=1,5
1491 : 132 : read(funit,*) (elastic_constants(mu,nu),mu=1,6)
1492 : : end do
1493 : : end if
1494 : 22 : read(funit,'(a)',iostat=ios) readline
1495 : 22 : call rmtabfromline(readline)
1496 : 22 : line=adjustl(readline)
1497 : 22 : call rdfromline_value('elastic',line,strg)
1498 [ - + ]: 22 : if (strg/="") then
1499 : 0 : strg1=trim(strg)
1500 : : else
1501 : 22 : strg1=trim(line)
1502 : : end if
1503 : 22 : read(strg1,*) (elastic_constants(mu,6),mu=1,6)
1504 : 22 : cycle
1505 : : end if
1506 : :
1507 [ + + ]: 5075 : if ((line(1:5)=='<atom')) then
1508 : 132 : call rdfromline("mass",line,strg)
1509 : 132 : strg1=trim(strg)
1510 : 132 : read(strg1,*) amu
1511 [ + + + + : 396 : if (.not.any(abs(all_amu-amu)<tol16)) then
+ + ]
1512 : 58 : all_amu(iamu) = amu
1513 : 58 : typat(iatom) = int(amu)
1514 : : !convert atomic mass unit to znucl
1515 [ + - ]: 1260 : do ii=1,103
1516 : 1260 : call atomdata_from_znucl(atom,real(ii,dp))
1517 : 1260 : if (abs((real(atom%amu,sp)-real(amu,sp))&
1518 [ + + ]: 1260 : & /real(amu,sp)*100)<0.1) then
1519 : 58 : znucl(iamu) = atom%znucl
1520 : 58 : exit
1521 : : end if
1522 : : end do
1523 : 58 : iamu = iamu +1
1524 : : end if
1525 [ + + ]: 496 : do itypat=1,ntypat
1526 [ + + ]: 496 : if(abs(amu-all_amu(itypat))<tol16) then
1527 : 132 : typat(iatom) = itypat
1528 : : end if
1529 : : end do
1530 : : cycle
1531 : : end if
1532 : :
1533 [ + + ]: 4943 : if ((line(1:9)=='<position')) then
1534 : 132 : call rdfromline_value('position',line,strg)
1535 [ - + ]: 132 : if (strg/="") then
1536 : 0 : strg1=trim(strg)
1537 : 0 : read(strg1,*)(xcart(mu,iatom),mu=1,3)
1538 : : else
1539 : 132 : read(funit,'(a)',iostat=ios) readline
1540 : 132 : call rmtabfromline(readline)
1541 : 132 : line=adjustl(readline)
1542 : 132 : call rdfromline_value('position',line,strg)
1543 [ - + ]: 132 : if (strg/="") then
1544 : 0 : strg1=trim(strg)
1545 : : else
1546 : 132 : strg1=trim(line)
1547 : : end if
1548 : 132 : read(strg1,*)(xcart(mu,iatom),mu=1,3)
1549 : : end if
1550 : : cycle
1551 : : end if
1552 : :
1553 [ + + ]: 4811 : if ((line(1:11)=='<borncharge')) then
1554 : 132 : call rdfromline_value('borncharge',line,strg)
1555 [ - + ]: 132 : if (strg/="") then
1556 : 0 : strg1=trim(strg)
1557 : 0 : read(strg1,*) (zeff(mu,1,iatom),mu=1,3)
1558 : 0 : read(funit,*) (zeff(mu,2,iatom),mu=1,3)
1559 : : else
1560 [ + + ]: 396 : do nu=1,2
1561 : 396 : read(funit,*) (zeff(mu,nu,iatom),mu=1,3)
1562 : : end do
1563 : : end if
1564 : 132 : read(funit,'(a)',iostat=ios) readline
1565 : 132 : line=adjustl(readline)
1566 : 132 : call rdfromline_value('borncharge',line,strg)
1567 [ - + ]: 132 : if (strg/="") then
1568 : 0 : strg1=trim(strg)
1569 : : else
1570 : 132 : strg1=trim(line)
1571 : : end if
1572 : 132 : read(strg1,*) (zeff(mu,3,iatom),mu=1,3)
1573 : 132 : cycle
1574 : : end if
1575 : :
1576 [ + + ]: 4679 : if ((line(1:7)==char(60)//char(47)//'atom'//char(62))) then
1577 : 132 : iatom=iatom+1
1578 : 132 : cycle
1579 : : end if
1580 : :
1581 [ + + ]: 4547 : if ((line(1:12)=='<local_force')) then
1582 : 648 : found2 = .False.
1583 : 648 : irpt1 = irpt1 + 1
1584 : : do while (.not.found2)
1585 : 1944 : read(funit,'(a)',iostat=ios) readline
1586 : 1944 : call rmtabfromline(readline)
1587 : 1944 : line=adjustl(readline)
1588 [ + + ]: 1944 : if ((line(1:5)=='<data')) then
1589 : 648 : call rdfromline_value('data',line,strg)
1590 [ - + ]: 648 : if (strg/="") then
1591 [ # # # # : 0 : ABI_MALLOC(work2,(3*natom,3*natom))
# # # # #
# # # # #
# # # # #
# ]
1592 : 0 : strg1=trim(strg)
1593 : 0 : read(strg1,*) (work2(1,nu),nu=1,3*natom)
1594 [ # # ]: 0 : do mu=2,3*natom-1
1595 : 0 : read(funit,*)(work2(mu,nu),nu=1,3*natom)
1596 : : end do
1597 : 0 : read(funit,'(a)',iostat=ios) readline
1598 : 0 : call rmtabfromline(readline)
1599 : 0 : line=adjustl(readline)
1600 : 0 : call rdfromline_value('data',line,strg)
1601 [ # # ]: 0 : if (strg/="") then
1602 : 0 : strg1=trim(strg)
1603 : : else
1604 : 0 : strg1=trim(line)
1605 : : end if
1606 [ # # # # ]: 0 : read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
1607 [ # # ]: 0 : local_atmfrc(:,:,:,:,irpt1) = reshape(work2,(/3,natom,3,natom/))
1608 [ # # ]: 0 : ABI_FREE(work2)
1609 : : else
1610 [ + - - + : 2592 : ABI_MALLOC(work2,(3*natom,3*natom))
+ - + - -
+ + - + -
+ - - + -
+ - + ]
1611 [ + + ]: 12726 : do mu=1,3*natom
1612 : 12726 : read(funit,*)(work2(mu,nu),nu=1,3*natom)
1613 : : end do
1614 [ + + ]: 3240 : local_atmfrc(:,:,:,:,irpt1) = reshape(work2,(/3,natom,3,natom/))
1615 [ - + ]: 648 : ABI_FREE(work2)
1616 : : end if
1617 : : end if
1618 [ + + ]: 1944 : if ((line(1:5)=='<cell')) then
1619 : 648 : call rdfromline_value('cell',line,strg)
1620 [ - + ]: 648 : if (strg/="") then
1621 : 0 : strg1=trim(strg)
1622 : 0 : read(strg1,*)(cell_local(mu,irpt1),mu=1,3)
1623 : : else
1624 : 648 : read(funit,*)(cell_local(mu,irpt1),mu=1,3)
1625 : : end if
1626 : : found2 = .TRUE.
1627 : : cycle
1628 : : end if
1629 : : end do
1630 : : end if
1631 : :
1632 [ + + ]: 4547 : if ((line(1:12)=='<total_force')) then
1633 : 531 : irpt2 = irpt2 + 1
1634 : 531 : found2 = .False.
1635 : : do while (.not.found2)
1636 : 1593 : read(funit,'(a)',iostat=ios) readline
1637 : 1593 : call rmtabfromline(readline)
1638 : 1593 : line=adjustl(readline)
1639 [ + + ]: 1593 : if ((line(1:5)=='<data')) then
1640 : 531 : call rdfromline_value('data',line,strg)
1641 [ - + ]: 531 : if (strg/="") then
1642 [ # # # # : 0 : ABI_MALLOC(work2,(3*natom,3*natom))
# # # # #
# # # # #
# # # # #
# ]
1643 : 0 : strg1=trim(strg)
1644 : 0 : read(strg1,*) (work2(1,nu),nu=1,3*natom)
1645 [ # # ]: 0 : do mu=2,3*natom-1
1646 : 0 : read(funit,*)(work2(mu,nu),nu=1,3*natom)
1647 : : end do
1648 : 0 : read(funit,'(a)',iostat=ios) readline
1649 : 0 : call rmtabfromline(readline)
1650 : 0 : line=adjustl(readline)
1651 : 0 : call rdfromline_value('data',line,strg)
1652 [ # # ]: 0 : if (strg/="") then
1653 : 0 : strg1=trim(strg)
1654 : : else
1655 : 0 : strg1=trim(line)
1656 : : end if
1657 [ # # # # ]: 0 : read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
1658 [ # # ]: 0 : total_atmfrc(:,:,:,:,irpt2) = reshape(work2,(/3,natom,3,natom/))
1659 [ # # ]: 0 : ABI_FREE(work2)
1660 : : else
1661 [ + - - + : 2124 : ABI_MALLOC(work2,(3*natom,3*natom))
+ - - + +
- + - + -
- + - + -
+ ]
1662 [ + + ]: 8208 : do mu=1,3*natom
1663 : 8208 : read(funit,*)(work2(mu,nu),nu=1,3*natom)
1664 : : end do
1665 [ + + ]: 2655 : total_atmfrc(:,:,:,:,irpt2) = reshape(work2,(/3,natom,3,natom/))
1666 [ - + ]: 531 : ABI_FREE(work2)
1667 : : end if
1668 : : end if
1669 [ + + ]: 1593 : if ((line(1:5)=='<cell')) then
1670 : 531 : call rdfromline_value('cell',line,strg)
1671 [ - + ]: 531 : if (strg/="") then
1672 : 0 : strg1=trim(strg)
1673 : 0 : read(strg1,*)(cell_total(mu,irpt2),mu=1,3)
1674 : : else
1675 : 531 : read(funit,*)(cell_total(mu,irpt2),mu=1,3)
1676 : : end if
1677 : : found2 = .TRUE.
1678 : : cycle
1679 : : end if
1680 : : end do
1681 : : end if
1682 : :
1683 [ + + ]: 4547 : if ((line(1:7)=='<qpoint')) then
1684 : 74 : call rdfromline_value('qpoint',line,strg)
1685 [ - + ]: 74 : if (strg/="") then
1686 : 0 : strg1=trim(strg)
1687 : 0 : read(strg1,*)(qph1l(mu,iph1l),mu=1,3)
1688 : : else
1689 : 74 : read(funit,*) (qph1l(mu,iph1l),mu=1,3)
1690 : : end if
1691 : : end if
1692 : :
1693 [ + + ]: 4547 : if ((line(1:12)=='<frequencies')) then
1694 : 74 : call rdfromline_value('frequencies',line,strg)
1695 [ - + ]: 74 : if (strg/="") then
1696 : 0 : strg1=trim(strg)
1697 : 0 : read(strg1,*)(phfrq(mu,iph1l),mu=1,3*natom)
1698 : : else
1699 [ + + ]: 410 : do nu=1,natom
1700 [ + - + + ]: 1418 : read(funit,*) (phfrq(((nu-1)*3)+mu,iph1l),mu=1,3)
1701 : : end do
1702 : : end if
1703 : : end if
1704 : :
1705 [ + + ]: 4547 : if ((line(1:17)=='<dynamical_matrix')) then
1706 : 74 : call rdfromline_value('dynamical_matrix',line,strg)
1707 [ - + ]: 74 : if (strg/="") then
1708 [ # # # # : 0 : ABI_MALLOC(work2,(3*natom,3*natom))
# # # # #
# # # # #
# # # # #
# ]
1709 : 0 : strg1=trim(strg)
1710 : 0 : read(strg1,*) (work2(nu,1),nu=1,3*natom)
1711 [ # # ]: 0 : do mu=2,3*natom-1
1712 : 0 : read(funit,*)(work2(nu,mu),nu=1,3*natom)
1713 : : end do
1714 : 0 : read(funit,'(a)',iostat=ios) readline
1715 : 0 : call rmtabfromline(readline)
1716 : 0 : line=adjustl(readline)
1717 : 0 : call rdfromline_value('dynamical_matrix',line,strg)
1718 [ # # ]: 0 : if (strg/="") then
1719 : 0 : strg1=trim(strg)
1720 : : else
1721 : 0 : strg1=trim(line)
1722 : : end if
1723 [ # # # # ]: 0 : read(strg1,*) (work2(nu,3*natom),nu=1,3*natom)
1724 [ # # # # : 0 : dynmat(1,:,:,:,:,iph1l) = reshape(work2,(/3,natom,3,natom/))
# # # # #
# ]
1725 [ # # ]: 0 : ABI_FREE(work2)
1726 : : else
1727 [ + - - + : 296 : ABI_MALLOC(work2,(3*natom,3*natom))
+ - - + +
- + - + -
- + - + -
+ ]
1728 [ + + ]: 1082 : do mu=1,3*natom
1729 : 1082 : read(funit,*)(work2(nu,mu),nu=1,3*natom)
1730 : : end do
1731 [ + + + + : 26002 : dynmat(1,:,:,:,:,iph1l) = reshape(work2,(/3,natom,3,natom/))
+ + + + +
+ ]
1732 [ - + ]: 74 : ABI_FREE(work2)
1733 : : end if
1734 : : end if
1735 : :
1736 [ + + ]: 4547 : if ((line(1:8)==char(60)//char(47)//'phonon')) then
1737 : 74 : iph1l = iph1l +1
1738 : : end if
1739 : :
1740 [ + + ]: 4547 : if ((line(1:16)=='<strain_coupling')) then
1741 : 22 : read(funit,'(a)',iostat=ios) readline
1742 : 22 : call rdfromline("voigt",line,strg)
1743 : 22 : strg1=trim(strg)
1744 : 22 : read(strg1,*) voigt
1745 : 22 : voigt = voigt + 1 ! 0 to 5 in the xml
1746 : 22 : has_straincoupling = .true.
1747 : 22 : irpt = 1
1748 : : end if
1749 : :
1750 : : else
1751 : : ! Now treat the strain phonon coupling part
1752 [ + + ]: 89326 : if ((line(1:16)=='<strain_coupling')) then
1753 : 110 : read(funit,'(a)',iostat=ios) readline
1754 : 110 : call rdfromline("voigt",line,strg)
1755 : 110 : strg1=trim(strg)
1756 : 110 : read(strg1,*) voigt
1757 : 110 : voigt = voigt + 1 ! 0 to 5 in the xml
1758 : 110 : irpt = 1
1759 : 110 : cycle
1760 : : end if
1761 : :
1762 [ - + ]: 89216 : if(voigt>6)then
1763 : 0 : write(message, '(4a)' )ch10,&
1764 : 0 : & ' WARNING: the number of strain phonon coupling is superior to 6 in ',trim(filename),ch10
1765 : 0 : call wrtout(std_out,message,'COLL')
1766 : 0 : exit
1767 : : end if
1768 : :
1769 [ + + ]: 89216 : if ((line(1:22)=='<correction_force unit')) then
1770 : 132 : call rdfromline_value('correction_force',line,strg)
1771 [ - + ]: 132 : if (strg/="") then
1772 [ # # # # : 0 : ABI_MALLOC(work2,(3,natom))
# # # # #
# # # #
# ]
1773 : 0 : strg1=trim(strg)
1774 : 0 : read(strg1,*) (work2(nu,1),nu=1,3)
1775 [ # # ]: 0 : do mu=2,natom-1
1776 : 0 : read(funit,*)(work2(nu,mu),nu=1,3)
1777 : : end do
1778 : 0 : read(funit,'(a)',iostat=ios) readline
1779 : 0 : call rmtabfromline(readline)
1780 : 0 : line=adjustl(readline)
1781 : 0 : call rdfromline_value('correction_force',line,strg)
1782 [ # # ]: 0 : if (strg/="") then
1783 : 0 : strg1=trim(strg)
1784 : 0 : read(strg1,*) (work2(nu,natom),nu=1,3)
1785 : : else
1786 : 0 : strg1=trim(line)
1787 : 0 : read(strg1,*) (work2(nu,natom),nu=1,3)
1788 : : end if
1789 [ # # # # ]: 0 : strain_coupling(voigt,:,:) = work2(:,:)
1790 [ # # ]: 0 : ABI_FREE(work2)
1791 : : else
1792 [ - + + - : 396 : ABI_MALLOC(work2,(3,natom))
+ - + - -
+ - + -
+ ]
1793 [ + + ]: 924 : do mu=1,natom
1794 : 924 : read(funit,*)(work2(nu,mu),nu=1,3)
1795 : : end do
1796 [ + + + + ]: 3300 : strain_coupling(voigt,:,:) = work2(:,:)
1797 [ - + ]: 132 : ABI_FREE(work2)
1798 : : end if
1799 : : end if
1800 : :
1801 [ - + ]: 89216 : if ((line(1:11)=='<elastic3rd')) then
1802 : 0 : call rdfromline_value('elastic3rd',line,strg)
1803 [ # # ]: 0 : if (strg/="") then
1804 : 0 : strg1=trim(strg)
1805 : 0 : read(strg1,*) (elastic3rd(voigt,mu,1),mu=1,6)
1806 [ # # ]: 0 : do nu=2,5
1807 : 0 : read(funit,*) (elastic3rd(voigt,mu,nu),mu=1,6)
1808 : : end do
1809 : : else
1810 [ # # ]: 0 : do nu=1,5
1811 : 0 : read(funit,*) (elastic3rd(voigt,mu,nu),mu=1,6)
1812 : : end do
1813 : : end if
1814 : 0 : read(funit,'(a)',iostat=ios) readline
1815 : 0 : call rmtabfromline(readline)
1816 : 0 : line=adjustl(readline)
1817 : 0 : call rdfromline_value('elastic3rd',line,strg)
1818 [ # # ]: 0 : if (strg/="") then
1819 : 0 : strg1=trim(strg)
1820 : 0 : read(strg1,*) (elastic3rd(voigt,mu,6),mu=1,6)
1821 : : else
1822 : 0 : strg1=trim(line)
1823 : 0 : read(strg1,*) (elastic3rd(voigt,mu,6),mu=1,6)
1824 : : end if
1825 : 0 : has_anharmonics = .true.
1826 : 0 : cycle
1827 : : end if
1828 : :
1829 [ - + ]: 89216 : if ((line(1:29)=='<correction_strain_force unit')) then
1830 : 0 : call rdfromline_value('correction_strain_force',line,strg)
1831 [ # # ]: 0 : if (strg/="") then
1832 [ # # # # : 0 : ABI_MALLOC(work2,(3*6,natom))
# # # # #
# # # #
# ]
1833 : 0 : strg1=trim(strg)
1834 : 0 : read(strg1,*) (work2(nu,1),nu=1,3*6)
1835 [ # # ]: 0 : do mu=2,natom-1
1836 : 0 : read(funit,*)(work2(nu,mu),nu=1,3*6)
1837 : : end do
1838 : 0 : read(funit,'(a)',iostat=ios) readline
1839 : 0 : call rmtabfromline(readline)
1840 : 0 : line=adjustl(readline)
1841 : 0 : call rdfromline_value('correction_strain_force',line,strg)
1842 [ # # ]: 0 : if (strg/="") then
1843 : 0 : strg1=trim(strg)
1844 : 0 : read(strg1,*) (work2(nu,natom),nu=1,3*6)
1845 : : else
1846 : 0 : strg1=trim(line)
1847 : 0 : read(strg1,*) (work2(nu,natom),nu=1,3*6)
1848 : : end if
1849 [ # # # # : 0 : elastic_displacement(voigt,:,:,:) = reshape(work2(:,:),(/6,3,natom/))
# # # # ]
1850 [ # # ]: 0 : ABI_FREE(work2)
1851 : : else
1852 [ # # # # : 0 : ABI_MALLOC(work2,(3*6,natom))
# # # # #
# # # #
# ]
1853 [ # # ]: 0 : do mu=1,natom
1854 : 0 : read(funit,*)(work2(nu,mu),nu=1,3*6)
1855 : : end do
1856 [ # # # # : 0 : elastic_displacement(voigt,:,:,:) = reshape(work2(:,:),(/6,3,natom/))
# # # # ]
1857 [ # # ]: 0 : ABI_FREE(work2)
1858 : : end if
1859 : : end if
1860 : :
1861 [ - + ]: 89216 : if ((line(1:26)=='<correction_force_constant')) then
1862 : : found2=.false.
1863 : : do while (.not.found2)
1864 : 0 : read(funit,'(a)',iostat=ios) readline
1865 : 0 : call rmtabfromline(readline)
1866 : 0 : line=adjustl(readline)
1867 [ # # ]: 0 : if ((line(1:5)=='<data')) then
1868 : 0 : call rdfromline_value('data',line,strg)
1869 [ # # ]: 0 : if (strg/="") then
1870 [ # # # # : 0 : ABI_MALLOC(work2,(3*natom,3*natom))
# # # # #
# # # # #
# # # # #
# ]
1871 : 0 : strg1=trim(strg)
1872 : 0 : read(strg1,*) (work2(1,nu),nu=1,3*natom)
1873 [ # # ]: 0 : do mu=2,3*natom-1
1874 : 0 : read(funit,*)(work2(mu,nu),nu=1,3*natom)
1875 : : end do
1876 : 0 : read(funit,'(a)',iostat=ios) readline
1877 : 0 : call rmtabfromline(readline)
1878 : 0 : line=adjustl(readline)
1879 : 0 : call rdfromline_value('data',line,strg)
1880 [ # # ]: 0 : if (strg/="") then
1881 : 0 : strg1=trim(strg)
1882 [ # # # # ]: 0 : read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
1883 : : else
1884 : 0 : strg1=trim(line)
1885 [ # # # # ]: 0 : read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
1886 : : end if
1887 : : phonon_strain(voigt)%atmfrc(:,:,:,:,irpt) = &
1888 [ # # # # : 0 : & reshape(work2,(/3,natom,3,natom/))
# # # # #
# ]
1889 [ # # ]: 0 : ABI_FREE(work2)
1890 : : else
1891 [ # # # # : 0 : ABI_MALLOC(work2,(3*natom,3*natom))
# # # # #
# # # # #
# # # # #
# ]
1892 [ # # ]: 0 : do mu=1,3*natom
1893 : 0 : read(funit,*)(work2(mu,nu),nu=1,3*natom)
1894 : : end do
1895 : : phonon_strain(voigt)%atmfrc(:,:,:,:,irpt) =&
1896 [ # # # # : 0 : & reshape(work2,(/3,natom,3,natom/))
# # # # #
# ]
1897 [ # # ]: 0 : ABI_FREE(work2)
1898 : : end if
1899 : 0 : has_anharmonics = .true.
1900 : : end if
1901 [ # # ]: 0 : if ((line(1:5)=='<cell')) then
1902 : 0 : call rdfromline_value('cell',line,strg)
1903 [ # # ]: 0 : if (strg/="") then
1904 : 0 : strg1=trim(strg)
1905 [ # # # # ]: 0 : read(strg1,*)(phonon_strain(voigt)%cell(mu,irpt),mu=1,3)
1906 : : else
1907 [ # # # # ]: 0 : read(funit,*)(phonon_strain(voigt)%cell(mu,irpt),mu=1,3)
1908 : : end if
1909 : 0 : irpt = irpt + 1
1910 : : found2=.true.
1911 : : cycle
1912 : : end if
1913 : : end do
1914 : : end if
1915 : :
1916 [ + + ]: 89216 : if ((line(1:17)==char(60)//char(47)//'strain_coupling')) then
1917 : : ! set nrpt for the previous value of strain
1918 : 132 : phonon_strain(voigt)%nrpt = irpt - 1
1919 : : ! restart the calculation of nrpt
1920 : : end if
1921 : : end if
1922 : : end if
1923 : : end do
1924 : :
1925 : :
1926 : : ! Reorder the ATMFRC
1927 : : ! Case 1: only local in the xml
1928 [ + + ]: 22 : if (irpt1>0 .and. irpt2==0) then
1929 [ + + + + ]: 763 : ifcs%cell(:,:) = int(cell_local(:,:))
1930 [ + + + + : 306916 : ifcs%atmfrc(:,:,:,:,:) = zero !local_atmfrc(:,:,:,:,:)
+ + + + +
+ ]
1931 [ + + + + : 306916 : ifcs%short_atmfrc(:,:,:,:,:) = local_atmfrc(:,:,:,:,:)
+ + + + +
+ ]
1932 [ + + + + : 306916 : ifcs%ewald_atmfrc(:,:,:,:,:) = zero
+ + + + +
+ ]
1933 : :
1934 : : ! Case 2: only total in the xml
1935 [ - + ]: 15 : else if(irpt1==0 .and. irpt2>0)then
1936 [ # # # # ]: 0 : ifcs%cell(:,:) = int(cell_total(:,:))
1937 [ # # # # : 0 : ifcs%atmfrc(:,:,:,:,:) = total_atmfrc(:,:,:,:,:)
# # # # #
# ]
1938 [ # # # # : 0 : ifcs%short_atmfrc(:,:,:,:,:) = zero
# # # # #
# ]
1939 [ # # # # : 0 : ifcs%ewald_atmfrc(:,:,:,:,:) = zero !total_atmfrc(:,:,:,:,:)
# # # # #
# ]
1940 : :
1941 : : ! Case 3: local + total in the xml
1942 [ + - ]: 15 : else if (irpt1>0 .and. irpt2>0)then
1943 [ + - ]: 15 : if(irpt1 <= irpt2)then
1944 : 15 : irpt3 = 0
1945 [ + + ]: 546 : do ii=1,irpt2
1946 [ + + ]: 2124 : ifcs%cell(:,ii) = int(cell_total(:,ii))
1947 [ + + + + : 159699 : ifcs%atmfrc(:,:,:,:,ii) = total_atmfrc(:,:,:,:,ii)
+ + + + ]
1948 [ + + ]: 36753 : do jj=1,irpt1
1949 [ + + + + : 46965 : if (all(abs(int(cell_local(:,jj))-ifcs%cell(:,ii))<tol16)) then
+ + ]
1950 [ + + + + : 144651 : ifcs%short_atmfrc(:,:,:,:,ii) = local_atmfrc(:,:,:,:,jj)
+ + + + ]
1951 : 459 : irpt3 = irpt3 + 1
1952 : : end if
1953 : : end do
1954 : : end do
1955 [ - + ]: 15 : if(irpt3 /= irpt1)then
1956 : 0 : write(message, '(4a)' )ch10,&
1957 : 0 : & ' There is several similar short IFC in ',trim(filename),ch10
1958 : 0 : ABI_BUG(message)
1959 : : end if
1960 : : else
1961 : 0 : write(message, '(2a,I5,3a,I5,5a)' )ch10,&
1962 : 0 : & ' The number of total IFC (',irpt2,') is inferior to ',ch10,&
1963 : 0 : & ' the number of short range IFC (',irpt1,') in ',trim(filename),ch10,&
1964 : 0 : & ' This is not possible',ch10
1965 : :
1966 : 0 : ABI_BUG(message)
1967 : : end if
1968 : : end if
1969 : :
1970 : : ! Do some checks
1971 [ + + + - : 154 : if (any(typat==0)) then
- + ]
1972 : : write(message, '(a,a,a)' )&
1973 : 0 : & ' Unable to read the type of atoms ',trim(filename),ch10
1974 : 0 : ABI_ERROR(message)
1975 : : end if
1976 : :
1977 [ + + + - : 80 : if (any(abs(znucl)<tol16)) then
- + ]
1978 : : write(message, '(a,a,a)' )&
1979 : 0 : & ' Unable to read the atomic number ',trim(filename),ch10
1980 : 0 : ABI_ERROR(message)
1981 : : end if
1982 : :
1983 [ + + + - : 80 : if (any(abs(all_amu)<tol16)) then
- + ]
1984 : : write(message, '(a,a,a)' )&
1985 : 0 : & ' Unable to read the atomic mass ',trim(filename),ch10
1986 : 0 : ABI_ERROR(message)
1987 : : end if
1988 : :
1989 : 22 : close(unit=funit)
1990 : :
1991 : : #endif
1992 : :
1993 : : end if !End if master
1994 : :
1995 : : !MPI BROADCAST
1996 : 38 : call xmpi_bcast(energy,master, comm, ierr)
1997 : 38 : call xmpi_bcast(all_amu,master, comm, ierr)
1998 : 38 : call xmpi_bcast(dynmat,master, comm, ierr)
1999 : 38 : call xmpi_bcast(elastic_constants,master, comm, ierr)
2000 : 38 : call xmpi_bcast(epsilon_inf,master, comm, ierr)
2001 : 38 : call xmpi_bcast(ifcs%nrpt,master, comm, ierr)
2002 : 38 : call xmpi_bcast(ifcs%atmfrc,master, comm, ierr)
2003 : 38 : call xmpi_bcast(ifcs%cell,master, comm, ierr)
2004 : 38 : call xmpi_bcast(ifcs%ewald_atmfrc,master, comm, ierr)
2005 : 38 : call xmpi_bcast(ifcs%short_atmfrc,master, comm, ierr)
2006 : 38 : call xmpi_bcast(strain_coupling,master, comm, ierr)
2007 : 38 : call xmpi_bcast(phfrq,master, comm, ierr)
2008 : 38 : call xmpi_bcast(qph1l,master, comm, ierr)
2009 : 38 : call xmpi_bcast(typat,master, comm, ierr)
2010 : 38 : call xmpi_bcast(rprimd,master, comm, ierr)
2011 : 38 : call xmpi_bcast(xcart,master, comm, ierr)
2012 : 38 : call xmpi_bcast(zeff,master, comm, ierr)
2013 : 38 : call xmpi_bcast(znucl,master, comm, ierr)
2014 [ + + ]: 266 : do ii = 1,6
2015 : 228 : call xmpi_bcast(phonon_strain(ii)%nrpt ,master, comm, ierr)
2016 : 228 : call xmpi_bcast(phonon_strain(ii)%atmfrc ,master, comm, ierr)
2017 : 266 : call xmpi_bcast(phonon_strain(ii)%cell ,master, comm, ierr)
2018 : : end do
2019 : 38 : call xmpi_bcast(elastic3rd ,master, comm, ierr)
2020 : 38 : call xmpi_bcast(elastic_displacement ,master, comm, ierr)
2021 : 38 : call xmpi_bcast(has_anharmonics ,master, comm, ierr)
2022 : :
2023 : : !Fill somes others variables
2024 : 38 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2025 : 38 : call xcart2xred(natom,rprimd,xcart,xred)
2026 : :
2027 : : !Re-generate symmetry operations from the lattice and atomic coordinates
2028 : 38 : tolsym=tol8
2029 : 38 : msym = 384
2030 [ - + + - : 114 : ABI_MALLOC(spinat,(3,natom))
+ - + - -
+ - + ]
2031 [ - + ]: 38 : ABI_MALLOC(ptsymrel,(3,3,msym))
2032 [ - + ]: 38 : ABI_MALLOC(symafm,(msym))
2033 [ - + ]: 38 : ABI_MALLOC(symrel,(3,3,msym))
2034 [ - + ]: 38 : ABI_MALLOC(tnons,(3,msym))
2035 : 38 : use_inversion=1
2036 [ + + + + ]: 854 : spinat = 0;
2037 [ + + + + : 189734 : symrel = 0;
+ + ]
2038 [ + + ]: 14630 : symafm = 0;
2039 [ + + + + ]: 58406 : tnons = 0 ;
2040 : 38 : space_group = 0;
2041 : 38 : call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym)
2042 : 38 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2043 : : call symfind(gprimd,msym,natom,nptsym,0,nsym,&
2044 : 38 : & 0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred)
2045 : :
2046 : : !Initialisation of crystal
2047 : 38 : npsp = ntypat; timrev = 1
2048 [ - + + - : 114 : ABI_MALLOC(title, (ntypat))
+ - + - -
+ - + ]
2049 [ + + ]: 136 : do ii=1,ntypat
2050 : 136 : write(title(ii),'(a,i0)')"No title for typat ",ii
2051 : : end do
2052 : :
2053 : : !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
2054 : : call crystal%init(all_amu,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,&
2055 : : & zion,znucl,timrev,.FALSE.,.FALSE.,title,&
2056 : 38 : & symrel=symrel(:,:,1:nsym),tnons=tnons(:,1:nsym),symafm=symafm(1:nsym))
2057 : :
2058 : : !amu is not fill in crystal_init...
2059 [ + + ]: 136 : Crystal%amu(:) = all_amu(:)
2060 : :
2061 [ - + ]: 38 : ABI_FREE(symrel)
2062 [ - + ]: 38 : ABI_FREE(symafm)
2063 [ - + ]: 38 : ABI_FREE(tnons)
2064 [ - + ]: 38 : ABI_FREE(spinat)
2065 [ - + ]: 38 : ABI_FREE(ptsymrel)
2066 : :
2067 : : !if strcpling is set to 0 by the user, need to set the flag to false for
2068 : : !the initialisation of the effective potential
2069 [ + - ]: 38 : if (present(strcpling))then
2070 [ - + ]: 38 : if(strcpling == 0 )then
2071 : 0 : has_anharmonics = .FALSE.
2072 : : end if
2073 : : end if
2074 : :
2075 : : !Initialisation of eff_pot
2076 : : call effective_potential_init(crystal,eff_pot,energy,ifcs,ncoeff,nph1l,comm,&
2077 : : & dynmat=dynmat,elastic_constants=elastic_constants,&
2078 : : & elastic3rd=elastic3rd,elastic_displacement=elastic_displacement,&
2079 : : & epsilon_inf=epsilon_inf,strain_coupling=strain_coupling,&
2080 : : & phonon_strain=phonon_strain,phfrq=phfrq,qpoints=qph1l,&
2081 : 38 : & has_anharmonicsTerms=has_anharmonics,zeff=zeff)
2082 : :
2083 : : !DEALLOCATION OF ARRAYS
2084 [ - + ]: 38 : ABI_FREE(all_amu)
2085 [ - + ]: 38 : ABI_FREE(cell_local)
2086 [ - + ]: 38 : ABI_FREE(cell_total)
2087 [ - + ]: 38 : ABI_FREE(total_atmfrc)
2088 [ - + ]: 38 : ABI_FREE(local_atmfrc)
2089 [ - + ]: 38 : ABI_FREE(ifcs%atmfrc)
2090 [ - + ]: 38 : ABI_FREE(ifcs%cell)
2091 [ - + ]: 38 : ABI_FREE(ifcs%short_atmfrc)
2092 [ - + ]: 38 : ABI_FREE(ifcs%ewald_atmfrc)
2093 [ - + ]: 38 : ABI_FREE(dynmat)
2094 [ - + ]: 38 : ABI_FREE(strain_coupling)
2095 [ - + ]: 38 : ABI_FREE(phfrq)
2096 [ - + ]: 38 : ABI_FREE(qph1l)
2097 [ - + ]: 38 : ABI_FREE(title)
2098 [ - + ]: 38 : ABI_FREE(typat)
2099 [ - + ]: 38 : ABI_FREE(xcart)
2100 [ - + ]: 38 : ABI_FREE(xred)
2101 [ - + ]: 38 : ABI_FREE(zeff)
2102 [ - + ]: 38 : ABI_FREE(zion)
2103 [ - + ]: 38 : ABI_FREE(znucl)
2104 [ + + ]: 266 : do ii = 1,6
2105 : 228 : phonon_strain(ii)%nrpt = nrpt
2106 [ - + - - : 228 : phonon_strain(ii)%atmfrc = zero
- - - - -
- ]
2107 [ - + - - ]: 228 : phonon_strain(ii)%cell = 0
2108 [ - + ]: 228 : ABI_FREE(phonon_strain(ii)%atmfrc)
2109 [ - + ]: 266 : ABI_FREE(phonon_strain(ii)%cell)
2110 : : end do
2111 [ + - + + : 266 : ABI_FREE(phonon_strain)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
2112 [ - + ]: 38 : ABI_FREE(elastic_displacement)
2113 : :
2114 : : !DEALLOCATION OF TYPES
2115 : 38 : call ifcs%free()
2116 : 38 : call crystal%free()
2117 : :
2118 [ - + - + : 114 : end subroutine system_xml2effpot
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
2119 : : !!***
2120 : :
2121 : : !!****f* m_effective_potential_file/system_ddb2effpot
2122 : : !!
2123 : : !! NAME
2124 : : !! system_ddb2effpot
2125 : : !!
2126 : : !! FUNCTION
2127 : : !! Transfert ddb into effective potential structure.
2128 : : !! Also calculate the IFC
2129 : : !!
2130 : : !! INPUTS
2131 : : !! crytal<type(crystal_t)> = datatype with all the information for the crystal
2132 : : !! ddb<type(ddb_type)> = datatype with the ddb
2133 : : !! inp<type(multibinit_dtset_type)> = datatype with the input variables of multibinit
2134 : : !! comm = MPI communicator
2135 : : !!
2136 : : !! OUTPUT
2137 : : !! effective_potantial<type(effective_potential_type)> = effective_potential datatype to be initialized
2138 : : !!
2139 : : !! SOURCE
2140 : :
2141 : 29 : subroutine system_ddb2effpot(crystal,ddb, effective_potential,inp,comm)
2142 : :
2143 : 38 : use m_dynmat
2144 : :
2145 : : use m_copy, only : alloc_copy
2146 : : use m_crystal, only : crystal_t
2147 : : use m_multibinit_dataset, only : multibinit_dtset_type
2148 : :
2149 : : !Arguments ------------------------------------
2150 : : !scalars
2151 : : integer,intent(in) :: comm
2152 : : !arrays
2153 : : type(ddb_type),intent(inout) :: ddb
2154 : : type(effective_potential_type), intent(inout) :: effective_potential
2155 : : type(crystal_t),intent(in) :: crystal
2156 : : type(multibinit_dtset_type),intent(in) :: inp
2157 : :
2158 : : !Local variables-------------------------------
2159 : : !scalar
2160 : : real(dp):: wcount1,wcount2
2161 : : integer :: chneut,i1,i2,i3,ia,ib,iblok,idir1,idir2,ierr,ii,ipert1,iphl1
2162 : : integer :: ipert2,irpt,irpt2,ivarA,ivarB,max1,max2,max3,min1,min2,min3
2163 : : integer :: msize,mpert,natom,nblok,nrpt_new,nrpt_new2,rftyp,selectz
2164 : : integer :: my_rank,nproc,prt_internalstr
2165 : : logical :: iam_master
2166 : : integer,parameter :: master=0
2167 : : integer :: nptsym,nsym
2168 : : integer :: msym = 384, use_inversion = 1, space_group
2169 : : real(dp):: max_phfq,eta, tolsym = tol8
2170 : : !arrays
2171 : : integer :: bravais(11),cell_number(3),cell2(3)
2172 : : integer :: shift(3),rfelfd(4),rfphon(4),rfstrs(4)
2173 : 29 : integer,allocatable :: cell_red(:,:)
2174 : : real(dp):: dielt(3,3),elast_clamped(6,6),fact
2175 : : real(dp):: red(3,3),qphnrm(3),qphon(3,3)
2176 : 29 : real(dp),allocatable :: blkval(:,:,:,:,:,:),d2asr(:,:,:,:,:)
2177 : 29 : real(dp),allocatable :: instrain(:,:),zeff(:,:,:),qdrp_cart(:,:,:,:)
2178 : 58 : real(dp),pointer :: atmfrc_red(:,:,:,:,:),wghatm_red(:,:,:)
2179 : : character(len=500) :: message
2180 : 29 : type(asrq0_t) :: asrq0
2181 [ + + + + ]: 174 : type(ifc_type) :: ifc
2182 : 29 : real(dp),allocatable :: d2cart(:,:,:,:,:),displ(:)
2183 : 58 : real(dp),allocatable :: eigval(:,:),eigvec(:,:,:,:,:),phfrq(:)
2184 : 29 : real(dp),allocatable :: spinat(:,:),tnons(:,:)
2185 : 29 : integer,allocatable :: symrel(:,:,:),symafm(:),ptsymrel(:,:,:)
2186 : : ! *************************************************************************
2187 : :
2188 : : !0 MPI variables
2189 : 29 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2190 : : iam_master=.FALSE.
2191 : 29 : iam_master = (my_rank == master)
2192 : :
2193 : : !Free the eff_pot before filling
2194 : 29 : call effective_potential_free(effective_potential)
2195 : :
2196 : : !Initialisation of usefull values
2197 : 29 : natom = ddb%natom
2198 : 29 : nblok = ddb%nblok
2199 : 29 : mpert= ddb%mpert
2200 : 29 : msize=3*mpert*3*mpert;
2201 : :
2202 : : !Tranfert the ddb into usable array (ipert and idir format like in abinit)
2203 [ - + + - : 174 : ABI_MALLOC(blkval,(2,3,mpert,3,mpert,nblok))
+ - - + +
- - + + -
+ - + - -
+ - + ]
2204 : :
2205 [ + + + + : 4992469 : blkval = 0
+ + + + +
+ + + ]
2206 [ + + - + ]: 116 : if(size(ddb%val) /= 2*3*mpert*3*mpert*nblok ) then
2207 : 0 : ABI_BUG("Size of ddb%val is not consistent.")
2208 : : endif
2209 [ + + ]: 203 : blkval = reshape(ddb%val,(/2,3,mpert,3,mpert,nblok/))
2210 : :
2211 : :
2212 : : !**********************************************************************
2213 : : ! Transfert crystal values
2214 : : !**********************************************************************
2215 : : ! Re-generate symmetry operations from the lattice and atomic coordinates
2216 [ - + + - : 87 : ABI_MALLOC(spinat,(3,natom))
+ - + - -
+ - + ]
2217 [ - + + - : 87 : ABI_MALLOC(ptsymrel,(3,3,msym))
+ - + - -
+ - + ]
2218 [ - + + - : 87 : ABI_MALLOC(symafm,(msym))
+ - + - -
+ - + ]
2219 [ - + + - : 58 : ABI_MALLOC(symrel,(3,3,msym))
- + - + ]
2220 [ - + + - : 87 : ABI_MALLOC(tnons,(3,msym))
+ - + - -
+ - + ]
2221 [ + + + + : 201085 : spinat = zero; symrel = 0; symafm = 0; tnons = zero ; space_group = 0;
+ + + + +
+ + + + +
+ + ]
2222 : 29 : call symlatt(bravais,std_out,msym,nptsym,ptsymrel,crystal%rprimd,tolsym)
2223 : : call symfind(crystal%gprimd,msym,crystal%natom,nptsym,0,nsym,&
2224 : : & 0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,&
2225 : 29 : & crystal%typat,use_inversion,crystal%xred)
2226 [ - + ]: 29 : if(crystal%nsym/=nsym)then
2227 : 0 : write(message,'(4a,I0,3a,I0,3a)') ch10,&
2228 : 0 : & ' --- !WARNING:',ch10,&
2229 : 0 : & ' There is ',nsym,' found for the crystal',ch10,&
2230 : 0 : & ' but ',crystal%nsym,' found in the DDB',ch10,&
2231 : 0 : & ' ---'
2232 : 0 : call wrtout(std_out,message,'COLL')
2233 : : end if
2234 : : call effective_potential%crystal%init(ddb%amu,&
2235 : : & space_group,crystal%natom,crystal%npsp,&
2236 : : & crystal%ntypat,nsym,crystal%rprimd,&
2237 : : & crystal%typat,crystal%xred,crystal%zion,&
2238 : : & crystal%znucl,crystal%timrev,crystal%use_antiferro,&
2239 : : & .FALSE.,crystal%title,&
2240 : : & symrel=symrel,tnons=tnons,&
2241 : 29 : & symafm=symafm)
2242 : :
2243 [ - + ]: 29 : ABI_FREE(spinat)
2244 [ - + ]: 29 : ABI_FREE(ptsymrel)
2245 [ - + ]: 29 : ABI_FREE(symafm)
2246 [ - + ]: 29 : ABI_FREE(symrel)
2247 [ - + ]: 29 : ABI_FREE(tnons)
2248 : :
2249 : : !**********************************************************************
2250 : : ! Transfert energy from input file
2251 : : !**********************************************************************
2252 [ + + ]: 2378 : write(message, '(2a,(80a),6a)') ch10,('=',ii=1,80),ch10,ch10,&
2253 [ + - ]: 2407 : & ' Extraction of the energy of the structure (unit: Hartree)',ch10
2254 : 29 : call wrtout(std_out,message,'COLL')
2255 : 29 : call wrtout(ab_out,message,'COLL')
2256 [ + + ]: 29 : if (ddb%get_etotal(effective_potential%energy) == 0) then
2257 [ - + ]: 17 : if(abs(inp%energy_reference) < tol16)then
2258 : : write(message,'(5a)')&
2259 : 0 : & ' Warning : Energy of the reference structure is not specify in',&
2260 : 0 : & ' the input file.',ch10,' Energy will set to zero',ch10
2261 : 0 : call wrtout(std_out,message,'COLL')
2262 : 0 : effective_potential%energy = zero
2263 : : else
2264 : 17 : effective_potential%energy = inp%energy_reference
2265 : : end if
2266 : : else
2267 [ + + ]: 12 : if(abs(inp%energy_reference) > tol16)then
2268 : : write(message,'(6a)')&
2269 : 1 : & ' Warning : Energy of the reference structure is specify in',&
2270 : 1 : & ' the input file.',ch10,' and in the DDB.',&
2271 : 2 : & ' The value of the energy is set with the value from the input file',ch10
2272 : 1 : call wrtout(std_out,message,'COLL')
2273 : 1 : effective_potential%energy = inp%energy_reference
2274 : : end if
2275 : : end if
2276 : 29 : write(message,'(a,es25.12)') ' Energy = ',&
2277 : 58 : & effective_potential%energy
2278 : 29 : call wrtout(std_out,message,'COLL')
2279 : 29 : call wrtout(ab_out,message,'COLL')
2280 : :
2281 : : !**********************************************************************
2282 : : ! Dielectric Tensor and Effective Charges
2283 : : !**********************************************************************
2284 [ - + + - : 87 : ABI_MALLOC(zeff,(3,3,natom))
+ - + - -
+ - + ]
2285 [ - + + - : 87 : ABI_MALLOC(qdrp_cart,(3,3,3,natom))
+ - + - -
+ - + ]
2286 [ - + + - : 58 : ABI_MALLOC(effective_potential%harmonics_terms%zeff,(3,3,natom))
- + - + -
+ ]
2287 : :
2288 : 29 : rftyp = 1 ! Blocks obtained by a non-stationary formulation.
2289 : 29 : chneut = 1 ! The ASR for effective charges is imposed
2290 : 29 : selectz = 0 ! No selection of some parts of the effective charge tensor
2291 : 29 : iblok = ddb%get_dielt_zeff(crystal,rftyp,chneut,selectz,dielt,zeff)
2292 [ + + + + : 6109 : qdrp_cart = zero
+ + + + ]
2293 [ + + + + : 377 : if (iblok /=0 .and. maxval(abs(dielt)) < 10000) then
+ + + - -
+ + + +
- ]
2294 [ + + + + ]: 247 : effective_potential%harmonics_terms%epsilon_inf = dielt
2295 [ + - + - : 1364 : effective_potential%harmonics_terms%zeff = zeff
+ - - + -
- + + + +
+ + ]
2296 : : else
2297 : 10 : effective_potential%harmonics_terms%epsilon_inf(1,1) = one
2298 : 10 : effective_potential%harmonics_terms%epsilon_inf(2,2) = one
2299 : 10 : effective_potential%harmonics_terms%epsilon_inf(3,3) = one
2300 [ + + + + : 660 : effective_potential%harmonics_terms%zeff = zero
+ + ]
2301 : : end if
2302 : :
2303 : : !**********************************************************************
2304 : : ! Look after the blok no. that contains the stress tensor
2305 : : !**********************************************************************
2306 [ + + ]: 2378 : write(message, '(a,a,(80a),a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2307 [ + - ]: 2407 : & ' Extraction of the stress tensor (unit: GPa) and forces (unit: Ha/bohr)'
2308 : 29 : call wrtout(std_out,message,'COLL')
2309 : 29 : call wrtout(ab_out,message,'COLL')
2310 : :
2311 [ - + + - : 58 : ABI_MALLOC(effective_potential%fcart,(3,natom))
- + - + -
+ ]
2312 [ + + + + ]: 637 : effective_potential%fcart = zero
2313 [ + + ]: 203 : effective_potential%strten = zero
2314 : :
2315 [ + + ]: 116 : qphon(:,1)=zero
2316 : 29 : qphnrm(1)=zero
2317 [ + + ]: 87 : rfphon(1:2)=0
2318 [ + + ]: 87 : rfelfd(1:2)=0
2319 [ + + ]: 87 : rfstrs(1:2)=0
2320 : 29 : rftyp=4
2321 : :
2322 : 29 : call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2323 : :
2324 [ + + ]: 29 : if (iblok /=0) then
2325 [ + + + - : 84 : if(any(abs(inp%strten_reference)>tol16))then
- + ]
2326 : 0 : write(message,'(10a)') ch10,&
2327 : 0 : & ' --- !WARNING:',ch10,&
2328 : 0 : & ' The stress tensor of the reference structure is specify in the',ch10,&
2329 : 0 : & ' input file and in the DDB. The value of the stress tensor is set',ch10,&
2330 : 0 : & ' with the value from the input file',ch10,&
2331 : 0 : & ' ---'
2332 : 0 : call wrtout(std_out,message,'COLL')
2333 : 0 : call wrtout(ab_out,message,'COLL')
2334 [ # # ]: 0 : effective_potential%strten = inp%strten_reference
2335 : : else
2336 : : ! firts give the corect stress values store in hartree
2337 : : ! diagonal parts
2338 : 12 : effective_potential%strten(1)=blkval(1,1,natom+3,1,1,iblok) * crystal%ucvol
2339 : 12 : effective_potential%strten(2)=blkval(1,2,natom+3,1,1,iblok) * crystal%ucvol
2340 : 12 : effective_potential%strten(3)=blkval(1,3,natom+3,1,1,iblok) * crystal%ucvol
2341 : : ! the shear parts
2342 : 12 : effective_potential%strten(4)=blkval(1,1,natom+4,1,1,iblok) * crystal%ucvol
2343 : 12 : effective_potential%strten(5)=blkval(1,2,natom+4,1,1,iblok) * crystal%ucvol
2344 : 12 : effective_potential%strten(6)=blkval(1,3,natom+4,1,1,iblok) * crystal%ucvol
2345 : : end if
2346 : : ! Get forces
2347 [ + + + + ]: 308 : effective_potential%fcart(:,1:natom) = blkval(1,:,1:natom,1,1,iblok)
2348 : : else
2349 [ + + + - : 119 : if(all(abs(inp%strten_reference(:))<tol16))then
+ - ]
2350 : 17 : write(message,'(8a)') ch10,&
2351 : 17 : & ' --- !WARNING:',ch10,&
2352 : 17 : & ' The stress tensor of the reference structure is not specify',ch10,&
2353 : 17 : & ' The stress tensor will be set to zero',ch10,&
2354 : 34 : & ' ---'
2355 : 17 : call wrtout(std_out,message,'COLL')
2356 : 17 : call wrtout(ab_out,message,'COLL')
2357 [ + + ]: 119 : effective_potential%strten = zero
2358 : : else
2359 [ # # ]: 0 : effective_potential%strten = inp%strten_reference
2360 : : end if
2361 : : end if
2362 : :
2363 [ + + + + : 131 : if(any(abs(effective_potential%strten(:)) >tol16))then
+ + ]
2364 : 12 : write(message, '(3a)' )ch10,&
2365 : 24 : & ' Cartesian components of forces (hartree/bohr)',ch10
2366 : 12 : call wrtout(ab_out,message,'COLL')
2367 : 12 : call wrtout(std_out, message,'COLL')
2368 [ + + ]: 86 : do ii = 1, natom
2369 : : write(message, '(I4,a,3(e16.8))' ) &
2370 : 74 : & ii,' ',effective_potential%fcart(:,ii)
2371 : :
2372 : 74 : call wrtout(ab_out,message,'COLL')
2373 : 86 : call wrtout(std_out, message,'COLL')
2374 : : end do
2375 : :
2376 : 12 : write(message, '(a,a)' )ch10,&
2377 : 24 : & ' Cartesian components of stress tensor (hartree/bohr^3)'
2378 : 12 : call wrtout(ab_out,message,'COLL')
2379 : 12 : call wrtout(std_out, message,'COLL')
2380 : : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2381 : 12 : & ' sigma(1 1)=',effective_potential%strten(1) / crystal%ucvol,&
2382 : 24 : & ' sigma(3 2)=',effective_potential%strten(4) / crystal%ucvol
2383 : 12 : call wrtout(ab_out,message,'COLL')
2384 : 12 : call wrtout(std_out, message,'COLL')
2385 : : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2386 : 12 : & ' sigma(2 2)=',effective_potential%strten(2) / crystal%ucvol,&
2387 : 24 : & ' sigma(3 1)=',effective_potential%strten(5) / crystal%ucvol
2388 : 12 : call wrtout(ab_out,message,'COLL')
2389 : 12 : call wrtout(std_out, message,'COLL')
2390 : : write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2391 : 12 : & ' sigma(3 3)=',effective_potential%strten(3) / crystal%ucvol,&
2392 : 24 : & ' sigma(2 1)=',effective_potential%strten(6) / crystal%ucvol
2393 : 12 : call wrtout(ab_out,message,'COLL')
2394 : 12 : call wrtout(std_out, message,'COLL')
2395 : 12 : write(message, '(a)' ) ' '
2396 : 12 : call wrtout(ab_out,message,'COLL')
2397 : 12 : call wrtout(std_out, message,'COLL')
2398 : : end if
2399 : :
2400 : : !**********************************************************************
2401 : : ! Elastic tensors at Gamma Point
2402 : : !**********************************************************************
2403 [ + + ]: 2378 : write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2404 [ + - ]: 2407 : & ' Extraction of the clamped elastic tensor (unit:10^2GPa)',ch10
2405 : 29 : call wrtout(std_out,message,'COLL')
2406 : 29 : call wrtout(ab_out,message,'COLL')
2407 : :
2408 : : ! look after the blok no.iblok that contains the elastic tensor
2409 [ + + ]: 116 : qphon(:,1)=zero
2410 : 29 : qphnrm(1)=zero
2411 [ + + ]: 87 : rfphon(1:2)=0
2412 [ + + ]: 87 : rfelfd(1:2)=0
2413 [ + + ]: 87 : rfstrs(1:2)=3 ! Need uniaxial both stresses and shear stresses
2414 : 29 : rftyp=1 ! Blocks obtained by a non-stationary formulation.
2415 : : ! for both diagonal and shear parts
2416 : 29 : call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2417 : :
2418 [ + - ]: 29 : if (iblok /=0) then
2419 : : ! extraction of the elastic constants from the blkvals (GPa)
2420 [ + + ]: 203 : do ivarA=1,6
2421 [ + + ]: 1247 : do ivarB=1,6
2422 : : ! because the elastic constant is 6*6,
2423 : : ! so we should judge if the idir is larger than 3
2424 : : ! or not
2425 [ + + ]: 1044 : if(ivarA>3) then
2426 : 522 : idir1=ivarA-3
2427 : 522 : ipert1=natom+4 !for the shear modulus
2428 : : else if(ivarA<=3) then
2429 : 522 : idir1=ivarA
2430 : 522 : ipert1=natom+3 !for the diagonal part
2431 : : end if
2432 [ + + ]: 1044 : if(ivarB>3) then
2433 : 522 : idir2=ivarB-3
2434 : 522 : ipert2=natom+4 !for the shear modulus
2435 : : else if(ivarB<=3) then
2436 : 522 : idir2=ivarB
2437 : 522 : ipert2=natom+3 !for the diagonal part
2438 : : end if
2439 : 1218 : elast_clamped(ivarA,ivarB) = blkval(1,idir1,ipert1,idir2,ipert2,iblok)
2440 : : end do
2441 : : end do
2442 : 29 : fact=HaBohr3_GPa / crystal%ucvol
2443 [ + + ]: 203 : do ivarA=1,6
2444 : 174 : write(message,'(6f12.7)')elast_clamped(ivarA,1)*fact/100.00_dp,&
2445 : 174 : & elast_clamped(ivarA,2)*fact/100.00_dp,&
2446 : 174 : & elast_clamped(ivarA,3)*fact/100.00_dp,&
2447 : 174 : & elast_clamped(ivarA,4)*fact/100.00_dp,&
2448 : 174 : & elast_clamped(ivarA,5)*fact/100.00_dp,&
2449 : 348 : & elast_clamped(ivarA,6)*fact/100.00_dp
2450 : 174 : call wrtout(std_out,message,'COLL')
2451 : 203 : call wrtout(ab_out,message,'COLL')
2452 : : end do
2453 : :
2454 : : ! Set the clamped tensor into the effective potentiel
2455 [ + + + + ]: 1247 : effective_potential%harmonics_terms%elastic_constants = elast_clamped
2456 : :
2457 : : else
2458 : :
2459 : 0 : write(message,'(3a)')ch10,&
2460 : 0 : & ' Warning : Elastic Tensor is set to zero (not available in the DDB)'
2461 : 0 : call wrtout(std_out,message,'COLL')
2462 : 0 : call wrtout(ab_out,message,'COLL')
2463 : :
2464 : : ! Set the clamped tensor to zero into the effective potentiel (not available in the DDB)
2465 [ # # # # ]: 0 : effective_potential%harmonics_terms%elastic_constants = zero
2466 : : end if
2467 : :
2468 : : !**********************************************************************
2469 : : ! Acoustic Sum Rule
2470 : : !***************************************************************************
2471 : : ! ASR-correction (d2asr) has to be determined here from the Dynamical matrix at Gamma.
2472 [ - + + - : 145 : ABI_MALLOC(d2asr,(2,3,natom,3,natom))
+ - - + +
- + - + -
- + - + ]
2473 : :
2474 [ + + ]: 2378 : write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2475 [ + - ]: 2407 : & ' Calculation of acoustic sum rule',ch10
2476 : 29 : call wrtout(std_out,message,'COLL')
2477 : 29 : call wrtout(ab_out,message,'COLL')
2478 : :
2479 : : ! Find the Gamma block in the DDB (no need for E-field entries)
2480 [ + + ]: 116 : qphon(:,1)=zero
2481 : 29 : qphnrm(1)=zero
2482 [ + + ]: 87 : rfphon(1:2)=1
2483 : 29 : rfelfd(:)=0
2484 : 29 : rfstrs(:)=0
2485 : 29 : rftyp=inp%rfmeth
2486 : :
2487 : 29 : call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2488 : :
2489 [ + + + + : 31477 : d2asr = zero
+ + + + +
+ ]
2490 [ + - ]: 29 : if (iblok /=0) then
2491 : 29 : call asria_calc(inp%asr,d2asr,ddb%val(:,:,iblok),ddb%mpert,ddb%natom)
2492 : : end if
2493 : :
2494 : : ! Acoustic sum rule imposition (not yet applied)
2495 : 29 : call asrq0%init(ddb, inp%asr, inp%rfmeth, crystal, inp%sys_dim)
2496 : :
2497 : : !**********************************************************************
2498 : : ! Interatomic Forces Calculation
2499 : : !**********************************************************************
2500 : : ! ifc to be calculated for interpolation
2501 [ + + ]: 2378 : write(message, '(a,a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,&
2502 [ + - ]: 2407 : & ' Calculation of the interatomic forces from DDB',ch10
2503 : 29 : call wrtout(std_out,message,'COLL')
2504 : 29 : call wrtout(ab_out,message,'COLL')
2505 : :
2506 : : call ifc%init(crystal,ddb,inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,&
2507 : : & inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,effective_potential%harmonics_terms%zeff,qdrp_cart,&
2508 : 29 : & inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,inp%sys_dim,comm)
2509 : : !***************************************************************************
2510 : : ! Interpolation of the dynamical matrix for each qpoint from ifc
2511 : : !***************************************************************************
2512 : :
2513 [ - + + - : 116 : ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
- + + - +
- + - - +
- + ]
2514 [ + - - + : 87 : ABI_MALLOC(displ,(2*3*natom*3*natom))
+ - + - +
- - + -
+ ]
2515 [ - + + - : 58 : ABI_MALLOC(eigval,(3,natom))
- + - + ]
2516 [ - + + - : 87 : ABI_MALLOC(eigvec,(2,3,natom,3,natom))
- + + - -
+ - + ]
2517 [ + - - + : 87 : ABI_MALLOC(phfrq,(3*natom))
+ - + - +
- - + -
+ ]
2518 : :
2519 [ - + + - : 145 : ABI_MALLOC(effective_potential%harmonics_terms%dynmat,(2,3,natom,3,natom,inp%nph1l))
- + + - -
+ + - + -
+ - - + -
+ - + ]
2520 [ - + + - : 116 : ABI_MALLOC(effective_potential%harmonics_terms%phfrq,(3*natom,inp%nph1l))
- + + - +
- + - - +
- + - + ]
2521 [ - + + - : 87 : ABI_MALLOC(effective_potential%harmonics_terms%qpoints,(3,inp%nph1l))
+ - + - -
+ - + -
+ ]
2522 : :
2523 [ + + ]: 2378 : write(message,'(a,(80a),3a)')ch10,('=',ii=1,80),ch10,ch10,&
2524 [ + - ]: 2407 : & ' Calculation of dynamical matrix for each ph1l points '
2525 : 29 : call wrtout(ab_out,message,'COLL')
2526 : 29 : call wrtout(std_out,message,'COLL')
2527 : :
2528 : : !Transfer value in effective_potential structure
2529 : 29 : effective_potential%harmonics_terms%nqpt = inp%nph1l
2530 [ + + + + ]: 369 : effective_potential%harmonics_terms%qpoints(:,:) = inp%qph1l(:,:)
2531 : :
2532 : : ! Store the highest frequency
2533 : 29 : max_phfq = zero
2534 : :
2535 [ + + ]: 114 : do iphl1=1,inp%nph1l
2536 : :
2537 : : ! Initialisation of the phonon wavevector
2538 [ + + ]: 340 : qphon(:,1)=inp%qph1l(:,iphl1)
2539 : 85 : if (inp%nph1l /= 0) qphnrm(1) = inp%qnrml1(iphl1)
2540 : :
2541 : : ! Get d2cart using the interatomic forces and the
2542 : : ! long-range coulomb interaction through Ewald summation
2543 : : call gtdyn9(ddb%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip,ifc%dyewq0,d2cart,crystal%gmet,&
2544 : : & ddb%gprim,mpert,natom,ifc%nrpt,qphnrm(1),qphon(:,1),crystal%rmet,ddb%rprim,ifc%rpt,&
2545 : : & ifc%trans,crystal%ucvol,ifc%wghatm,crystal%xred,zeff,qdrp_cart,ifc%ewald_option,eta,&
2546 : 85 : & xmpi_comm_self,ifc%sys_dim)
2547 : :
2548 : : ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
2549 : : call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,crystal%indsym,&
2550 : : & mpert,crystal%nsym,natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),qphon,&
2551 : 85 : & crystal%rprimd,inp%symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol)
2552 : :
2553 : : ! Write the phonon frequencies
2554 : 85 : call dfpt_prtph(displ,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
2555 : :
2556 : : ! Store the highest frequency in cmm-1
2557 [ + - - + : 1298 : max_phfq = max(maxval(phfrq*Ha_cmm1),max_phfq)
- - + + ]
2558 : :
2559 [ + + + + : 59309 : effective_potential%harmonics_terms%dynmat(:,:,:,:,:,iphl1) = d2cart(:,:,:natom,:,:natom)
+ + + + +
+ ]
2560 [ + + ]: 1242 : effective_potential%harmonics_terms%phfrq(:,iphl1) = phfrq(:) * Ha_cmm1
2561 : :
2562 : : end do
2563 : :
2564 : 29 : write(message, '(2a,f15.7,a)' ) ch10,&
2565 : 58 : & ' The highest frequency found is ',max_phfq,' cm-1'
2566 : 29 : call wrtout(std_out,message,'COLL')
2567 : :
2568 [ - + ]: 29 : ABI_FREE(d2cart)
2569 [ - + ]: 29 : ABI_FREE(displ)
2570 [ - + ]: 29 : ABI_FREE(eigval)
2571 [ - + ]: 29 : ABI_FREE(eigvec)
2572 [ - + ]: 29 : ABI_FREE(phfrq)
2573 : :
2574 : : !**********************************************************************
2575 : : ! Transfert inter-atomic forces constants in reduced coordinates
2576 : : !**********************************************************************
2577 : :
2578 : : !Reorder cell from canonical coordinates to reduced coordinates (for multibinit)
2579 : : !store the number of ifc before rearrangement
2580 : :
2581 : : ! Store the sum of the weight of IFC for the final check
2582 : 29 : wcount1 = 0
2583 [ + + ]: 5112 : do irpt=1,ifc%nrpt
2584 [ + + + + ]: 166452 : wcount1 = wcount1 + sum(ifc%wghatm(:,:,irpt))
2585 : : end do
2586 : :
2587 : : !Set the maximum and the miminum for the bound of the cell
2588 [ + + + + ]: 10195 : max1 = maxval(ifc%cell(1,:)); min1 = minval(ifc%cell(1,:))
2589 [ + + + + ]: 10195 : max2 = maxval(ifc%cell(2,:)); min2 = minval(ifc%cell(2,:))
2590 [ + + + + ]: 10195 : max3 = maxval(ifc%cell(3,:)); min3 = minval(ifc%cell(3,:))
2591 : 29 : cell_number(1) = max1 - min1 + 1
2592 : 29 : cell_number(2) = max2 - min2 + 1
2593 : 29 : cell_number(3) = max3 - min3 + 1
2594 : :
2595 : : ! set the new number of cell, sometimes, in canonical coordinates,
2596 : : ! some cell are delete but they exist in reduced coordinates.
2597 [ + + ]: 116 : nrpt_new = product(cell_number(:))
2598 : :
2599 : : ! Allocate temporary array
2600 [ - + + - : 174 : ABI_MALLOC(atmfrc_red,(3,natom,3,natom,nrpt_new))
+ - - + +
- - + + -
+ - + - -
+ - + ]
2601 [ - + + - : 145 : ABI_MALLOC(wghatm_red,(natom,natom,nrpt_new))
- + + - -
+ + - + -
+ - - + -
+ ]
2602 [ - + + - : 87 : ABI_MALLOC(cell_red,(3,nrpt_new))
+ - + - -
+ - + -
+ ]
2603 : :
2604 [ + + + + : 169320 : wghatm_red(:,:,:) = zero
+ + ]
2605 : :
2606 [ + + ]: 29 : if(iam_master)then
2607 [ + + ]: 137 : do ia=1,natom
2608 [ + + ]: 1001 : do ib=1,natom
2609 : :
2610 : : ! Simple Lattice
2611 [ + - ]: 864 : if (inp%brav==1) then
2612 : : ! In this case, it is better to work in reduced coordinates
2613 : : ! As rcan is in canonical coordinates, => multiplication by gprim
2614 [ + + ]: 3456 : do ii=1,3
2615 : : red(1,ii)= ifc%rcan(1,ia)*ddb%gprim(1,ii) + &
2616 : : & ifc%rcan(2,ia)*ddb%gprim(2,ii) + &
2617 : 2592 : & ifc%rcan(3,ia)*ddb%gprim(3,ii)
2618 : : red(2,ii)= ifc%rcan(1,ib)*ddb%gprim(1,ii) + &
2619 : : ifc%rcan(2,ib)*ddb%gprim(2,ii) + &
2620 : 3456 : & ifc%rcan(3,ib)*ddb%gprim(3,ii)
2621 : : end do
2622 : : end if
2623 : :
2624 : : ! Get the shift of cell
2625 [ + + ]: 3456 : shift(:) = int(anint(red(2,:) - crystal%xred(:,ib)) - anint(red(1,:) - crystal%xred(:,ia)))
2626 : :
2627 [ + + ]: 133224 : do irpt=1,ifc%nrpt
2628 : :
2629 [ + + ]: 528976 : cell2(:)= int(ifc%cell(:,irpt) + shift(:))
2630 : :
2631 : : ! Use boundary condition to get the right cell
2632 [ + + + - ]: 132244 : if (cell2(1) < min1 .and. cell2(1) < max1) then
2633 : 915 : cell2(1) = cell2(1) + cell_number(1)
2634 [ + + + + ]: 131329 : else if (cell2(1) > min1 .and. cell2(1) > max1) then
2635 : 915 : cell2(1) = cell2(1) - cell_number(1)
2636 : : end if
2637 : :
2638 [ + + + - ]: 132244 : if (cell2(2) < min2 .and. cell2(2) < max2) then
2639 : 960 : cell2(2) = cell2(2) + cell_number(2)
2640 [ + + + + ]: 131284 : else if (cell2(2) > min2 .and. cell2(2) > max2) then
2641 : 960 : cell2(2) = cell2(2) - cell_number(2)
2642 : : end if
2643 : :
2644 [ + + + - ]: 132244 : if (cell2(3) < min3 .and. cell2(3) < max3) then
2645 : 948 : cell2(3) = cell2(3) + cell_number(3)
2646 [ + + + + ]: 131296 : else if (cell2(3) > min3 .and. cell2(3) > max3) then
2647 : 948 : cell2(3) = cell2(3) - cell_number(3)
2648 : : end if
2649 : :
2650 : : irpt2=1
2651 [ + + ]: 1198440 : do i1=min1,max1
2652 [ + + ]: 10360972 : do i2=min2,max2
2653 [ + + ]: 91175516 : do i3=min3,max3
2654 : : if (i1 == cell2(1) .and.&
2655 [ + + + + : 80946788 : i2 == cell2(2) .and.&
+ + ]
2656 : : i3 == cell2(3)) then
2657 : 132244 : wghatm_red(ia,ib,irpt2) = ifc%wghatm(ia,ib,irpt)
2658 [ + + + + ]: 1719172 : atmfrc_red(:,ia,:,ib,irpt2) = ifc%atmfrc(:,ia,:,ib,irpt)
2659 : 132244 : cell_red(1,irpt2) = i1
2660 : 132244 : cell_red(2,irpt2) = i2
2661 : 132244 : cell_red(3,irpt2) = i3
2662 : : end if
2663 : 90110184 : irpt2 = irpt2 + 1
2664 : : end do
2665 : : end do
2666 : : end do
2667 : : end do
2668 : : end do
2669 : : end do
2670 : : end if
2671 : :
2672 : 29 : call xmpi_bcast(atmfrc_red,master, comm, ierr)
2673 : 29 : call xmpi_bcast(wghatm_red,master, comm, ierr)
2674 : 29 : call xmpi_bcast(cell_red,master, comm, ierr)
2675 : :
2676 : : ! Copy ifc into effective potential
2677 : : ! !!Warning eff_pot%ifcs only contains atmfrc,short_atmfrc,ewald_atmfrc,,nrpt and cell!!
2678 : : ! rcan,ifc%rpt,wghatm and other quantities
2679 : : ! are not needed for effective potential!!!
2680 : 29 : call ifc%free()
2681 : 29 : call effective_potential%harmonics_terms%ifcs%free()
2682 : :
2683 : : ! Only conserve the necessary points in rpt
2684 : 29 : nrpt_new2 = 0
2685 [ + + ]: 5220 : do irpt = 1, nrpt_new
2686 [ + + + + : 169320 : if (abs(sum(wghatm_red(:,:,irpt))) >tol16) then
+ + ]
2687 : 5003 : nrpt_new2 = nrpt_new2 + 1
2688 : : end if
2689 : : end do
2690 : :
2691 : : ! Set the new number of rpt
2692 : 29 : effective_potential%harmonics_terms%ifcs%nrpt = nrpt_new2
2693 : :
2694 : : ! Allocation of the final arrays
2695 [ - + + - : 174 : ABI_MALLOC(effective_potential%harmonics_terms%ifcs%atmfrc,(3,natom,3,natom,nrpt_new2))
+ - - + +
- - + + -
+ - + - -
+ - + -
+ ]
2696 [ - + + - : 116 : ABI_MALLOC(effective_potential%harmonics_terms%ifcs%short_atmfrc,(3,natom,3,natom,nrpt_new2))
- + + - -
+ + - - +
- + - + ]
2697 [ - + + - : 116 : ABI_MALLOC(effective_potential%harmonics_terms%ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt_new2))
- + + - -
+ + - - +
- + - + ]
2698 [ - + + - : 87 : ABI_MALLOC(effective_potential%harmonics_terms%ifcs%cell,(3,nrpt_new2))
+ - + - -
+ - + -
+ ]
2699 [ - + + - : 145 : ABI_MALLOC(effective_potential%harmonics_terms%ifcs%wghatm,(natom,natom,nrpt_new2))
- + + - -
+ + - + -
+ - - + -
+ - + ]
2700 : :
2701 : 29 : irpt2 = 0
2702 [ + + ]: 5220 : do irpt = 1,nrpt_new
2703 [ + + + + : 169320 : if (abs(sum(wghatm_red(:,:,irpt))) > tol16) then
+ + ]
2704 : 5003 : irpt2 = irpt2 + 1
2705 : : ! Apply weight on each R point
2706 [ + + ]: 30335 : do ia=1,effective_potential%crystal%natom
2707 [ + + ]: 164743 : do ib=1,effective_potential%crystal%natom
2708 [ + + + + ]: 1772636 : atmfrc_red(:,ia,:,ib,irpt) = atmfrc_red(:,ia,:,ib,irpt)*wghatm_red(ia,ib,irpt)
2709 : : end do
2710 : : end do
2711 [ + + ]: 20012 : effective_potential%harmonics_terms%ifcs%cell(:,irpt2) = cell_red(:,irpt)
2712 [ + + + + : 1719227 : effective_potential%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt2) = atmfrc_red(:,:,:,:,irpt)
+ + + + ]
2713 [ + + ]: 5003 : if (inp%dipdip == 1) then
2714 : : effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)=&
2715 [ + + + + : 1520030 : & atmfrc_red(:,:,:,:,irpt)
+ + + + ]
2716 : : else
2717 [ + + + + : 199197 : effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2) = zero
+ + + + ]
2718 : : end if
2719 [ + + + + : 1719227 : effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)=atmfrc_red(:,:,:,:,irpt)
+ + + + ]
2720 [ + + + + : 1719227 : effective_potential%harmonics_terms%ifcs%ewald_atmfrc(:,:,:,:,irpt2) = zero
+ + + + ]
2721 [ + + + + ]: 164743 : effective_potential%harmonics_terms%ifcs%wghatm(:,:,irpt2) = wghatm_red(:,:,irpt)
2722 : : end if
2723 : : end do
2724 : :
2725 : :
2726 [ - + ]: 29 : ABI_FREE(atmfrc_red)
2727 [ - + ]: 29 : ABI_FREE(wghatm_red)
2728 [ - + ]: 29 : ABI_FREE(cell_red)
2729 : :
2730 : : ! Final check
2731 : 29 : wcount2 = 0
2732 [ + + ]: 5032 : do irpt = 1, effective_potential%harmonics_terms%ifcs%nrpt
2733 [ + + + + ]: 164772 : wcount2 = wcount2 + sum(effective_potential%harmonics_terms%ifcs%wghatm(:,:,irpt))
2734 : : end do
2735 : :
2736 [ - + ]: 29 : if (abs(wcount1-wcount2)/(wcount1+wcount2)>tol8) then
2737 : 0 : write(message,'(2a,es15.4,a,es15.4,a,es15.4)')'The total wghatm has changed',ch10,&
2738 : 0 : & wcount1,' before and ', wcount2, ' now, difference being ',wcount1-wcount2
2739 : 0 : ABI_BUG(message)
2740 : : end if
2741 : :
2742 : :
2743 : : !**********************************************************************
2744 : : ! Internal strain tensors at Gamma point
2745 : : !**********************************************************************
2746 [ + + ]: 2378 : write(message, '(a,a,(80a),a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2747 [ + - ]: 2407 : & ' Calculation of the internal-strain tensor'
2748 : 29 : call wrtout(std_out,message,'COLL')
2749 : 29 : call wrtout(ab_out,message,'COLL')
2750 [ + - - + : 116 : ABI_MALLOC(instrain,(3*natom,6))
+ - + - +
- + - - +
- + ]
2751 : : ! looking after the no. of blok that contains the internal strain tensor
2752 [ + + ]: 116 : qphon(:,1)=zero
2753 : 29 : qphnrm(1)=zero
2754 [ + + ]: 87 : rfphon(1:2)=0
2755 : : rfelfd(1:2)=0
2756 [ + + ]: 87 : rfstrs(1:2)=3
2757 : 29 : rftyp=1
2758 : 29 : call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2759 : :
2760 [ - + + - : 87 : ABI_MALLOC(effective_potential%harmonics_terms%strain_coupling,(6,3,natom))
+ - + - -
+ - + -
+ ]
2761 [ + + + + : 3373 : effective_potential%harmonics_terms%strain_coupling = zero
+ + ]
2762 : :
2763 [ + - ]: 29 : if (iblok /=0) then
2764 : :
2765 : : ! then print the internal strain tensor (only the force one)
2766 : 29 : prt_internalstr=1
2767 : : call ddb_internalstr(inp%asr,ddb%val,d2asr,iblok,instrain,&
2768 : 29 : & ab_out,mpert,natom,nblok,prt_internalstr)
2769 : :
2770 [ + + ]: 203 : do ipert1=1,6
2771 [ + + ]: 1115 : do ipert2=1,natom
2772 [ + + ]: 3822 : do idir2=1,3
2773 : 2736 : ii=3*(ipert2-1)+idir2
2774 : : effective_potential%harmonics_terms%strain_coupling(ipert1,idir2,ipert2)=&
2775 : 3648 : & (-1.0_dp)*instrain(ii,ipert1)
2776 : : end do
2777 : : end do
2778 : : end do
2779 : : else
2780 : 0 : write(message,'(3a)')ch10,&
2781 : 0 : & ' Warning : Internal strain is set to zero (not available in the DDB)'
2782 : 0 : call wrtout(std_out,message,'COLL')
2783 : 0 : call wrtout(ab_out,message,'COLL')
2784 : : end if
2785 : : !-------------------------------------------------------------------------------------
2786 : : ! DEALLOCATION OF ARRAYS
2787 [ - + ]: 29 : ABI_FREE(blkval)
2788 [ - + ]: 29 : ABI_FREE(zeff)
2789 [ - + ]: 29 : ABI_FREE(qdrp_cart)
2790 [ - + ]: 29 : ABI_FREE(instrain)
2791 [ - + ]: 29 : ABI_FREE(d2asr)
2792 : 29 : call asrq0%free()
2793 : :
2794 : 29 : write(message,'(a)')ch10
2795 : 29 : call wrtout(std_out,message,'COLL')
2796 : 29 : call wrtout(ab_out,message,'COLL')
2797 : :
2798 [ - + - + : 87 : end subroutine system_ddb2effpot
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
2799 : : !!***
2800 : :
2801 : : !!****f* m_effective_potential_file/coeffs_xml2effpot
2802 : : !! NAME
2803 : : !! coeffs_xml2effpot
2804 : : !!
2805 : : !! FUNCTION
2806 : : !! Open xml file of effective potentiel, then reads the variables
2807 : : !! and store them in effective potentential type
2808 : : !!
2809 : : !! INPUTS
2810 : : !! filename = path of input or output file
2811 : : !! comm=MPI communicator
2812 : : !!
2813 : : !! OUTPUT
2814 : : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
2815 : : !!
2816 : : !! SOURCE
2817 : :
2818 : 28 : subroutine coeffs_xml2effpot(eff_pot,filename,comm)
2819 : :
2820 : 29 : use m_atomdata
2821 : : use m_polynomial_coeff
2822 : : use m_polynomial_term
2823 : : use m_crystal, only : symbols_crystal
2824 : : #if defined HAVE_XML
2825 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_PTR,c_f_pointer
2826 : : #endif
2827 : :
2828 : : !Arguments ------------------------------------
2829 : : !scalars
2830 : : character(len=*),intent(in) :: filename
2831 : : integer, intent(in) :: comm
2832 : : !arrays
2833 : : type(effective_potential_type), intent(inout) :: eff_pot
2834 : :
2835 : : !Local variables-------------------------------
2836 : : !scalar
2837 : : integer :: ii,jj,my_rank,ndisp,ncoeff,nterm_max,nstrain,ndisp_max,nproc,nterm
2838 : : ! character(len=200),allocatable :: name(:)
2839 : : character(len=200) :: name
2840 : : #ifdef HAVE_XML
2841 : : integer :: icoeff,iterm
2842 : : #endif
2843 : :
2844 : : #ifndef HAVE_XML
2845 : : integer :: funit = 1,ios = 0
2846 : : integer :: icoeff,idisp,istrain,iterm,mu
2847 : : logical :: found,found2,displacement
2848 : : character (len=XML_RECL) :: line,readline
2849 : : character (len=XML_RECL) :: strg,strg1
2850 : : #endif
2851 : : character(len=500) :: message
2852 : : character(len=264) :: filename_tmp
2853 : 28 : character(len=5),allocatable :: symbols(:)
2854 : : integer,parameter :: master=0
2855 : : logical :: iam_master
2856 : : logical :: debug
2857 : : !arrays
2858 : 28 : real(dp),allocatable :: coefficient(:),weight(:,:)
2859 : 28 : integer,allocatable :: atindx(:,:,:,:), cell(:,:,:,:,:),direction(:,:,:),power_disp(:,:,:)
2860 : 28 : integer,allocatable :: strain(:,:,:),power_strain(:,:,:)
2861 : 28 : type(polynomial_coeff_type),dimension(:),allocatable :: coeffs
2862 : 28 : type(polynomial_term_type),dimension(:,:),allocatable :: terms
2863 : : ! *************************************************************************
2864 : :
2865 : 28 : filename_tmp = trim(filename)
2866 : : !Open the atomicdata XML file for reading
2867 : 28 : write(message,'(a,a)')'-Opening the file ',trim(filename_tmp)
2868 : :
2869 : 28 : call wrtout(ab_out,message,'COLL')
2870 : 28 : call wrtout(std_out,message,'COLL')
2871 : :
2872 : 28 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2873 : 28 : iam_master = (my_rank == master)
2874 : :
2875 : : !Get Dimention of system and allocation/initialisation of array
2876 : : ncoeff = 0
2877 : 28 : nterm = 0
2878 : 28 : ndisp = 0
2879 : 28 : nstrain = 0
2880 : 28 : call effective_potential_file_getDimCoeff(filename,ncoeff,ndisp_max,nterm_max)
2881 : :
2882 : : ! Do some checks
2883 [ - + ]: 28 : if (nterm_max<=0) then
2884 : : write(message, '(a,a,a)' )&
2885 : 0 : & ' Unable to read the number of terms in ',trim(filename),ch10
2886 : 0 : ABI_ERROR(message)
2887 : : end if
2888 : :
2889 [ - + ]: 28 : if (ndisp_max<=0) then
2890 : : write(message, '(a,a,a)' )&
2891 : 0 : & ' Unable to read the number of displacement in ',trim(filename),ch10
2892 : 0 : ABI_ERROR(message)
2893 : : end if
2894 : :
2895 : : !Allocation ov the polynomial coeff type
2896 [ - + + - : 1886 : ABI_MALLOC(coeffs,(ncoeff))
+ - + - -
+ - + + +
+ - - + -
- - - - -
- - - - -
- - - -
- ]
2897 : :
2898 [ + + ]: 28 : if(iam_master)then
2899 : :
2900 : : #if defined HAVE_XML
2901 : : write(message,'(3a)')'-Reading the file ',trim(filename),&
2902 : : & ' with LibXML library'
2903 : : #else
2904 : 20 : write(message,'(3a)')'-Reading the file ',trim(filename),&
2905 : 40 : & ' with Fortran'
2906 : : #endif
2907 : 20 : call wrtout(ab_out,message,'COLL')
2908 : 20 : call wrtout(std_out,message,'COLL')
2909 : :
2910 : :
2911 [ - + + - : 60 : ABI_MALLOC(symbols,(eff_pot%crystal%natom))
+ - + - -
+ - + -
+ ]
2912 : : ! Get the symbols arrays
2913 : : call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,&
2914 : 20 : & eff_pot%crystal%npsp,symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl)
2915 : :
2916 : :
2917 : : !Read with libxml librarie
2918 : : #if defined HAVE_XML
2919 : :
2920 : : ABI_MALLOC(terms,(ncoeff,nterm_max))
2921 : : ABI_MALLOC(atindx,(ncoeff,nterm_max,2,ndisp_max))
2922 : : ABI_MALLOC(coefficient,(ncoeff))
2923 : : ABI_MALLOC(cell,(ncoeff,nterm_max,3,2,ndisp_max))
2924 : : ABI_MALLOC(direction,(ncoeff,nterm_max,ndisp_max))
2925 : : ABI_MALLOC(strain,(ncoeff,nterm_max,ndisp_max))
2926 : : ABI_MALLOC(power_disp,(ncoeff,nterm_max,ndisp_max))
2927 : : ABI_MALLOC(power_strain,(ncoeff,nterm_max,ndisp_max))
2928 : : ABI_MALLOC(weight,(ncoeff,nterm_max))
2929 : :
2930 : : ! Read the values of this term with libxml
2931 : : call effpot_xml_readCoeff(char_f2c(trim(filename)),ncoeff,ndisp_max,nterm_max,&
2932 : : & coefficient,atindx,cell,direction,power_disp,power_strain,&
2933 : : & strain,weight)
2934 : : ! In the XML the atom index begin to zero
2935 : : ! Need to shift for fortran array
2936 : : atindx(:,:,:,:) = atindx(:,:,:,:) + 1
2937 : :
2938 : : do icoeff=1,ncoeff
2939 : : do iterm=1,nterm_max
2940 : : ! Initialisation of the polynomial_term structure with the values from the
2941 : : call polynomial_term_init(atindx(icoeff,iterm,:,:),cell(icoeff,iterm,:,:,:),&
2942 : : & direction(icoeff,iterm,:),ndisp_max,ndisp_max,terms(icoeff,iterm),&
2943 : : & power_disp(icoeff,iterm,:),power_strain(icoeff,iterm,:),&
2944 : : & strain(icoeff,iterm,:),weight(icoeff,iterm),check=.true.)
2945 : : end do
2946 : : ! Initialisation of the polynomial_coefficent structure with the values
2947 : : call polynomial_coeff_init(coefficient(icoeff),nterm_max,coeffs(icoeff),&
2948 : : & terms(icoeff,:), check=.true., debug_str="init from xml")
2949 : : ! Set the name of the coefficient
2950 : :
2951 : : ! Get the name of this coefficient and set it
2952 : : ! Try to find the index of the term corresponding to the interation in the
2953 : : ! reference cell (000) in order to compute the name correctly...
2954 : : ! If this coeff is not in the ref cell, take by default the first term
2955 : : if(coeffs(icoeff)%nterm > 0)then
2956 : : call polynomial_coeff_getName(name,coeffs(icoeff),symbols,recompute=.true.)
2957 : : call polynomial_coeff_setName(name,coeffs(icoeff))
2958 : : end if
2959 : :
2960 : : ! Free them all
2961 : : do iterm=1,nterm_max
2962 : : call polynomial_term_free(terms(icoeff,iterm))
2963 : : end do
2964 : : end do
2965 : :
2966 : : #else
2967 [ - + + - : 1884 : ABI_MALLOC(terms,(1,nterm_max))
+ - + - -
+ - + + +
+ + - + -
+ - + - +
- + - + -
+ - + ]
2968 [ - + + - : 60 : ABI_MALLOC(atindx,(1,1,2,ndisp_max))
+ - + - -
+ - + ]
2969 [ - + ]: 20 : ABI_MALLOC(coefficient,(1))
2970 [ - + + - : 60 : ABI_MALLOC(cell,(1,1,3,2,ndisp_max))
+ - + - -
+ - + ]
2971 [ - + + - : 60 : ABI_MALLOC(direction,(1,1,ndisp_max))
+ - + - -
+ - + ]
2972 [ - + + - : 40 : ABI_MALLOC(strain,(1,1,ndisp_max))
- + - + ]
2973 [ - + + - : 40 : ABI_MALLOC(power_disp,(1,1,ndisp_max))
- + - + ]
2974 [ - + + - : 40 : ABI_MALLOC(power_strain,(1,1,ndisp_max))
- + - + ]
2975 [ - + - + ]: 20 : ABI_MALLOC(weight,(1,1))
2976 : : ! Loop over the file
2977 : : ! Read the values of all the terms with fortran
2978 [ - + ]: 20 : if (open_file(filename,message,unit=funit,form="formatted",&
2979 : : & status="old",action="read") /= 0) then
2980 : 0 : ABI_ERROR(message)
2981 : : end if
2982 : :
2983 : : ! Start a reading loop in fortran
2984 : 20 : rewind(unit=funit)
2985 : 20 : ios = 0
2986 : 20 : found=.false.
2987 : :
2988 : : ! Initialisation of counter
2989 : 20 : icoeff = 0
2990 : :
2991 : : ! Parser
2992 [ + + ]: 4966 : do while (ios==0)
2993 : 4946 : read(funit,'(a)',iostat=ios) readline
2994 [ + + ]: 4966 : if (ios == 0) then
2995 : 4926 : call rmtabfromline(readline)
2996 : 4926 : line=adjustl(readline)
2997 [ + + ]: 4926 : if ((line(1:12)==char(60)//'coefficient')) then
2998 : : ! Read headers of coefficient
2999 : 1770 : call rdfromline('text',line,strg)
3000 [ + - ]: 1770 : if (strg/="") then
3001 : 1770 : name=trim(strg)
3002 : : end if
3003 : 1770 : call rdfromline('value',line,strg)
3004 [ + + ]: 1770 : if (strg/="") then
3005 : 1754 : strg1=trim(strg)
3006 : 1754 : read(strg1,*) coefficient(1)
3007 : : else
3008 : 16 : coefficient(1) = zero
3009 : : end if
3010 : : ! End read headers of coefficient
3011 : : ! Reset counter
3012 : 8833 : found = .false.
3013 [ + + + + : 220723 : atindx = 0; cell = 0 ; direction = 0
+ + + + +
+ + + + +
+ + + + +
+ + + +
+ ]
3014 [ + + + + : 65337 : strain = 0; power_strain = 0; power_disp = 0
+ + + + +
+ + + + +
+ + + + ]
3015 : 1770 : iterm = 0
3016 : 1770 : idisp = 0
3017 : 1770 : istrain = 0
3018 : 1770 : nterm = 0
3019 : : do while (.not.found)
3020 : 48012 : read(funit,'(a)',iostat=ios) readline
3021 : 48012 : call rmtabfromline(readline)
3022 : 48012 : line=adjustl(readline)
3023 [ + + ]: 48012 : if ((line(1:13)==char(60)//'/coefficient')) then
3024 : : found= .true.
3025 : : cycle
3026 : : end if
3027 [ - + ]: 46242 : if ((line(1:5)==char(60)//'term')) then
3028 : 46242 : nterm = nterm + 1
3029 : 46242 : ndisp = 0
3030 : 46242 : nstrain = 0
3031 : 46242 : idisp = 0
3032 : 46242 : istrain = 0
3033 : 46242 : displacement = .true.
3034 : 46242 : call rdfromline('weight',line,strg)
3035 [ + - ]: 46242 : if (strg/="") then
3036 : 46242 : strg1=trim(strg)
3037 : 46242 : read(strg1,*) weight
3038 : : end if
3039 [ + + ]: 214908 : do while(displacement)
3040 : 168666 : read(funit,'(a)',iostat=ios) readline
3041 : 168666 : call rmtabfromline(readline)
3042 : 168666 : line=adjustl(readline)
3043 [ + + ]: 168666 : if ((line(1:6)==char(60)//'/term')) then
3044 : 46242 : displacement = .false.
3045 : : end if
3046 [ + + ]: 168666 : if ((line(1:7)==char(60)//'strain')) then
3047 : 14682 : nstrain = nstrain + 1
3048 : 14682 : istrain = istrain + 1
3049 : 14682 : call rdfromline('power',line,strg)
3050 [ + - ]: 14682 : if (strg/="") then
3051 : 14682 : strg1=trim(strg)
3052 : 14682 : read(strg1,*) power_strain(1,1,istrain)
3053 : : end if
3054 : 14682 : call rdfromline('voigt',line,strg)
3055 [ + - ]: 14682 : if (strg/="") then
3056 : 14682 : strg1=trim(strg)
3057 : 14682 : read(strg1,*) strain(1,1,istrain)
3058 : : end if
3059 : : end if
3060 [ + + ]: 214908 : if ((line(1:18)==char(60)//'displacement_diff')) then
3061 : 107742 : ndisp = ndisp + 1
3062 : 107742 : idisp = idisp + 1
3063 : 107742 : found2=.true.
3064 : 107742 : call rdfromline('atom_a',line,strg)
3065 [ + - ]: 107742 : if (strg/="") then
3066 : 107742 : strg1=trim(strg)
3067 : 107742 : read(strg1,*) atindx(1,1,1,idisp)
3068 : : end if
3069 : 107742 : call rdfromline('atom_b',line,strg)
3070 [ + - ]: 107742 : if (strg/="") then
3071 : 107742 : strg1=trim(strg)
3072 : 107742 : read(strg1,*) atindx(1,1,2,idisp)
3073 : : end if
3074 : 107742 : call rdfromline('direction',line,strg)
3075 [ + - ]: 107742 : if (strg/="") then
3076 : 107742 : strg1=trim(strg)
3077 [ + + ]: 107742 : if (trim(strg1).eq."x") direction(1,1,idisp) = 1
3078 [ + + ]: 107742 : if (trim(strg1).eq."y") direction(1,1,idisp) = 2
3079 [ + + ]: 107742 : if (trim(strg1).eq."z") direction(1,1,idisp) = 3
3080 : : end if
3081 : 107742 : call rdfromline('power',line,strg)
3082 [ + - ]: 107742 : if (strg/="") then
3083 : 107742 : strg1=trim(strg)
3084 : 107742 : read(strg1,*) power_disp(1,1,idisp)
3085 : : end if
3086 : : do while(found2)
3087 : 323226 : read(funit,'(a)',iostat=ios) readline
3088 : 323226 : call rmtabfromline(readline)
3089 : 323226 : line=adjustl(readline)
3090 [ + + ]: 323226 : if ((line(1:7)==char(60)//'cell_a')) then
3091 : 107742 : call rdfromline_value('cell_a',line,strg)
3092 [ + - ]: 107742 : if (strg/="") then
3093 : 107742 : strg1=trim(strg)
3094 : 107742 : read(strg1,*) (cell(1,1,mu,1,idisp),mu=1,3)
3095 : : else
3096 : 0 : read(funit,'(a)',iostat=ios) readline
3097 : 0 : call rmtabfromline(readline)
3098 : 0 : line=adjustl(readline)
3099 : 0 : call rdfromline_value('cell_a',line,strg)
3100 [ # # ]: 0 : if (strg/="") then
3101 : 0 : strg1=trim(strg)
3102 : 0 : read(strg1,*)(cell(1,1,mu,1,idisp),mu=1,3)
3103 : : else
3104 : 0 : strg1=trim(line)
3105 : 0 : read(strg1,*)(cell(1,1,mu,1,idisp),mu=1,3)
3106 : : end if
3107 : : end if
3108 : : end if
3109 [ + + ]: 323226 : if ((line(1:7)==char(60)//'cell_b')) then
3110 : 107742 : call rdfromline_value('cell_b',line,strg)
3111 [ + - ]: 107742 : if (strg/="") then
3112 : 107742 : strg1=trim(strg)
3113 : 107742 : read(strg1,*) (cell(1,1,mu,2,idisp),mu=1,3)
3114 : : else
3115 : 0 : read(funit,'(a)',iostat=ios) readline
3116 : 0 : call rmtabfromline(readline)
3117 : 0 : line=adjustl(readline)
3118 : 0 : call rdfromline_value('cell_b',line,strg)
3119 [ # # ]: 0 : if (strg/="") then
3120 : 0 : strg1=trim(strg)
3121 : 0 : read(strg1,*)(cell(1,1,mu,2,idisp),mu=1,3)
3122 : : else
3123 : 0 : strg1=trim(line)
3124 : 0 : read(strg1,*)(cell(1,1,mu,2,idisp),mu=1,3)
3125 : : end if
3126 : : end if
3127 : : end if
3128 [ + + ]: 323226 : if ((line(1:19)==char(60)//'/displacement_diff')) then
3129 : : found2=.false.
3130 : : end if
3131 : : end do
3132 : : end if
3133 : : end do!end do while displacement
3134 : : ! In the XML the atom index begin to zero
3135 : : ! Need to shift for fortran array
3136 [ + + + + ]: 599490 : atindx(1,1,:,:) = atindx(1,1,:,:) + 1
3137 : : ! Initialisation of the polynomial_term structure with the values from the
3138 : : ! previous step
3139 : 46242 : iterm = iterm + 1
3140 : : call polynomial_term_init(atindx(1,1,:,:),cell(1,1,:,:,:),&
3141 : : & direction(1,1,:),ndisp,nstrain,terms(1,iterm),&
3142 : : & power_disp(1,1,:),power_strain(1,1,:),&
3143 : 46242 : & strain(1,1,:),weight(1,1),check=.true.)
3144 : : end if!end if term
3145 : : end do!end do while found (coeff)
3146 : :
3147 : : ! Initialisation of the polynomial_coefficent structure with the values from the
3148 : : ! previous step
3149 : 1770 : icoeff = icoeff + 1
3150 : 1770 : call polynomial_coeff_init(coefficient(1),nterm,coeffs(icoeff),terms(1,:), debug_str="init from xml fortran", check=.true.)
3151 : 1770 : call polynomial_coeff_getName(name,coeffs(icoeff),symbols,recompute=.true.)
3152 : 1770 : call polynomial_coeff_setName(name,coeffs(icoeff))
3153 : : ! Deallocation of the terms array for this coefficient
3154 [ + + ]: 84330 : do jj=1,nterm_max
3155 : 84330 : call polynomial_term_free(terms(1,jj))
3156 : : end do
3157 : : end if!end if line = coefficient
3158 : : end if!end if ios==0
3159 : : end do!end do while on file
3160 : :
3161 : 20 : close(unit=funit)
3162 : :
3163 : : #endif
3164 [ + - + + : 932 : ABI_FREE(terms)
- + - + -
+ - + - +
- + - + -
+ ]
3165 [ - + ]: 20 : ABI_FREE(atindx)
3166 [ - + ]: 20 : ABI_FREE(coefficient)
3167 [ - + ]: 20 : ABI_FREE(cell)
3168 [ - + ]: 20 : ABI_FREE(direction)
3169 [ - + ]: 20 : ABI_FREE(strain)
3170 [ - + ]: 20 : ABI_FREE(power_disp)
3171 [ - + ]: 20 : ABI_FREE(power_strain)
3172 [ - + ]: 20 : ABI_FREE(weight)
3173 [ - + ]: 20 : ABI_FREE(symbols)
3174 : : end if !End if master
3175 : :
3176 : : !9-MPI BROADCAST
3177 [ + + ]: 1830 : do ii=1,ncoeff
3178 : 1830 : call polynomial_coeff_broadcast(coeffs(ii),master, comm)
3179 : : end do
3180 : :
3181 : : !10-checks
3182 : :
3183 : : !11-debug print
3184 : 28 : debug = .FALSE.
3185 : : if(debug)then
3186 : : do ii=1,ncoeff
3187 : : do jj=1,coeffs(ii)%nterm
3188 : : #if defined HAVE_XML
3189 : : write(200+my_rank,*)"ii,jj,ndisp,nterm",ii,jj,coeffs(ii)%nterm,coeffs(ii)%terms(jj)%ndisp
3190 : : write(200+my_rank,*)"atindx",coeffs(ii)%terms(jj)%atindx
3191 : : write(200+my_rank,*)"cell1",coeffs(ii)%terms(jj)%cell(:,1,:)
3192 : : write(200+my_rank,*)"cell2",coeffs(ii)%terms(jj)%cell(:,2,:)
3193 : : write(200+my_rank,*)"direction",coeffs(ii)%terms(jj)%direction
3194 : : write(200+my_rank,*)"power_disp",coeffs(ii)%terms(jj)%power_disp
3195 : : write(200+my_rank,*)"weight",coeffs(ii)%terms(jj)%weight
3196 : : #else
3197 : : write(300+my_rank,*)"ii,jj,ndisp,nterm",ii,jj,coeffs(ii)%nterm,coeffs(ii)%terms(jj)%ndisp
3198 : : write(300+my_rank,*)"atindx",coeffs(ii)%terms(jj)%atindx
3199 : : write(300+my_rank,*)"cell1",coeffs(ii)%terms(jj)%cell(:,1,:)
3200 : : write(300+my_rank,*)"cell2",coeffs(ii)%terms(jj)%cell(:,2,:)
3201 : : write(300+my_rank,*)"direction",coeffs(ii)%terms(jj)%direction
3202 : : write(300+my_rank,*)"power_disp",coeffs(ii)%terms(jj)%power_disp
3203 : : write(300+my_rank,*)"weight",coeffs(ii)%terms(jj)%weight
3204 : : #endif
3205 : : end do
3206 : : end do
3207 : : #if defined HAVE_XML
3208 : : close(200+my_rank)
3209 : : #else
3210 : : close(300+my_rank)
3211 : : #endif
3212 : : end if
3213 : :
3214 : : !12-Initialisation of eff_pot
3215 : 28 : call effective_potential_setCoeffs(coeffs,eff_pot,ncoeff)
3216 : :
3217 : : !13-Deallocation of type
3218 [ + + ]: 1830 : do ii=1,ncoeff
3219 : 1830 : call polynomial_coeff_free(coeffs(ii))
3220 : : end do
3221 [ + - + + : 1830 : ABI_FREE(coeffs)
- + - - -
- - - - -
- - - - -
- - - -
+ ]
3222 : :
3223 : :
3224 [ - + - + : 84 : end subroutine coeffs_xml2effpot
- + ]
3225 : : !!***
3226 : :
3227 : : !!****f* m_effective_potential_file/effective_potential_file_readMDfile
3228 : : !!
3229 : : !! NAME
3230 : : !! effective_potential_file_readMDfile
3231 : : !!
3232 : : !! FUNCTION
3233 : : !! Read MD FILE (HIST or ASCII)
3234 : : !!
3235 : : !! INPUTS
3236 : : !! filename = path of the file
3237 : : !! option,optional = 0 (default), the stress is printed in the MD File
3238 : : !! 1, the force on the cell is printed in the MD File (-1 * stress),
3239 : : !! in this case, we multiply the stress by -1 in order to get the stresse
3240 : : !!
3241 : : !! OUTPUT
3242 : : !! hist<type(abihist)> = datatype with the history of the MD
3243 : : !!
3244 : : !! SOURCE
3245 : :
3246 : 13 : subroutine effective_potential_file_readMDfile(filename,hist,option)
3247 : :
3248 : : !Arguments ------------------------------------
3249 : : !scalars
3250 : : integer,optional :: option
3251 : : !arrays
3252 : : type(abihist),intent(inout) :: hist
3253 : : character(len=fnlen),intent(in) :: filename
3254 : : !Local variables-------------------------------
3255 : : !scalar
3256 : : integer :: ia,ii,mu,nu,natom,nstep,type,option_in
3257 : : integer :: ios=0, unit_md=24
3258 : : !arrays
3259 : : character (len=10000) :: readline,line
3260 : : real(dp) :: tmp(6)
3261 : 13 : real(dp),allocatable :: xcart(:,:)
3262 : : ! *************************************************************************
3263 : :
3264 : 13 : call effective_potential_file_getType(filename,type)
3265 : :
3266 : 13 : option_in = 0
3267 : : if(present(option))then
3268 : : option_in = option
3269 : : end if
3270 [ + - ]: 13 : if(type==40)then
3271 : : ! Netcdf type
3272 : 13 : call read_md_hist(filename,hist,.FALSE.,.FALSE.,.FALSE.)
3273 : :
3274 [ # # ]: 0 : else if(type==41)then
3275 : :
3276 : : ! ASCII file
3277 : 0 : call effective_potential_file_getDimMD(filename,natom,nstep)
3278 : :
3279 : 0 : ii = 1
3280 : 0 : ios = 0
3281 : :
3282 [ # # # # : 0 : ABI_MALLOC(xcart,(3,natom))
# # # # #
# # # ]
3283 : 0 : call abihist_free(hist)
3284 : 0 : call abihist_init(hist,natom,nstep,.FALSE.,.FALSE.)
3285 : :
3286 : : ! Start a reading loop in fortran
3287 : 0 : rewind(unit=unit_md)
3288 [ # # # # ]: 0 : do while ((ios==0).and.ii<=nstep)
3289 : 0 : read(unit_md,'(a)',iostat=ios) readline
3290 : 0 : read(unit_md,'(a)',iostat=ios) readline
3291 : 0 : line=adjustl(readline)
3292 : 0 : read(line,*) hist%etot(ii)
3293 : : hist%etot(ii) = hist%etot(ii)
3294 [ # # ]: 0 : do mu=1,3
3295 : 0 : read(unit_md,'(a)',iostat=ios) readline
3296 : 0 : line=adjustl(readline)
3297 [ # # # # ]: 0 : read(line,*) (hist%rprimd(nu,mu,ii),nu=1,3)
3298 : : end do
3299 [ # # ]: 0 : do ia=1,natom
3300 : 0 : read(unit_md,'(a)',iostat=ios) readline
3301 : 0 : line=adjustl(readline)
3302 : 0 : read(line,*) (tmp(mu),mu=1,6)
3303 [ # # ]: 0 : xcart(:,ia) = tmp(1:3)
3304 [ # # ]: 0 : hist%fcart(:,ia,ii) = tmp(4:6)
3305 : : end do
3306 : 0 : call xcart2xred(natom,hist%rprimd(:,:,ii),xcart(:,:),hist%xred(:,:,ii))
3307 : 0 : read(unit_md,'(a)',iostat=ios) readline
3308 : 0 : line=adjustl(readline)
3309 [ # # # # ]: 0 : read(line,*) (hist%strten(mu,ii),mu=1,6)
3310 : 0 : ii = ii + 1
3311 : : end do
3312 [ # # ]: 0 : do ii=1,nstep
3313 [ # # ]: 0 : do mu=1,3
3314 [ # # ]: 0 : hist%acell(mu,:) = hist%rprimd(mu,mu,ii)
3315 : : end do
3316 : : end do
3317 : 0 : close(unit_md)
3318 [ # # ]: 0 : ABI_FREE(xcart)
3319 : :
3320 : : end if!end if type
3321 : :
3322 [ + - - + ]: 13 : if((type==40 .or. type==41).and.option == 1)then
3323 : : ! multiply by -1 if the current strten -1*stress, we need only stress...
3324 [ # # # # ]: 0 : hist%strten(:,:) = -1 * hist%strten(:,:)
3325 : : end if
3326 : :
3327 : :
3328 : :
3329 : 28 : end subroutine effective_potential_file_readMDfile
3330 : : !!***
3331 : :
3332 : : !!****f* m_effective_potential_file/effective_potential_file_mapHistToRef
3333 : : !!
3334 : : !! NAME
3335 : : !! effective_potential_file_mapHistToRef
3336 : : !!
3337 : : !! FUNCTION
3338 : : !! Generate the supercell in the effective potential according to the size of the
3339 : : !! supercell in the hist file
3340 : : !! Check if the hist file match to reference supercell in the effective potential
3341 : : !! If not, the hist file is reordering
3342 : : !!
3343 : : !! INPUTS
3344 : : !! eff_pot<type(effective_potential)> = effective potential
3345 : : !! hist<type(abihist)> = The history of the MD
3346 : : !! comm = MPI communicator
3347 : : !!
3348 : : !! OUTPUT
3349 : : !! hist<type(abihist)> = The history of the MD
3350 : : !!
3351 : : !! SOURCE
3352 : :
3353 : 34 : subroutine effective_potential_file_mapHistToRef(eff_pot,hist,comm,iatfix,verbose,sc_size)
3354 : :
3355 : : !Arguments ------------------------------------
3356 : : !scalars
3357 : : integer,intent(in) :: comm
3358 : : logical,optional,intent(in) :: verbose
3359 : : !arrays
3360 : : type(effective_potential_type),intent(inout) :: eff_pot
3361 : : type(abihist),intent(inout) :: hist
3362 : : integer,optional,allocatable,intent(inout) :: iatfix(:,:)
3363 : : integer,optional,intent(in) :: sc_size(3)
3364 : : !Local variables-------------------------------
3365 : : !scalar
3366 : : integer :: factE_hist,ia,ib,ii,jj,natom_hist,ncells,nstep_hist
3367 : : real(dp):: factor,ratio
3368 : : logical :: revelant_factor,need_map,need_verbose,need_fixmap
3369 : : !arrays
3370 : : real(dp) :: rprimd_hist(3,3),rprimd_ref(3,3)
3371 : : integer :: ncell(3),scale_cell(3)
3372 : 34 : integer,allocatable :: shift(:,:),iatfix_tmp(:,:)
3373 : 34 : integer,allocatable :: list_map(:) !blkval(:),
3374 : 34 : real(dp),allocatable :: xred_ref(:,:) ! xred_hist(:,:),
3375 : 34 : real(dp),allocatable :: list_dist(:),list_reddist(:,:),list_absdist(:,:)
3376 : : character(len=500) :: msg
3377 : 34 : type(abihist) :: hist_tmp
3378 : : ! *************************************************************************
3379 : :
3380 : : !Set optional values
3381 : 34 : need_verbose = .false.
3382 : 34 : need_fixmap = .FALSE.
3383 [ + + ]: 34 : if (present(verbose)) need_verbose = verbose
3384 [ + + ]: 34 : if (present(iatfix)) need_fixmap = .TRUE.
3385 : :
3386 : :
3387 : 34 : natom_hist = size(hist%xred,2)
3388 : 34 : nstep_hist = size(hist%xred,3)
3389 : :
3390 : : ! Try to set the supercell according to the hist file
3391 [ + + + + ]: 442 : rprimd_ref(:,:) = eff_pot%crystal%rprimd
3392 [ + + + + ]: 442 : rprimd_hist(:,:) = hist%rprimd(:,:,1)
3393 : :
3394 : :
3395 [ + + ]: 34 : if(present(sc_size))then
3396 : 12 : ncell(:) = sc_size
3397 : : else
3398 [ + + ]: 88 : do ia=1,3
3399 : 66 : scale_cell(:) = 0
3400 [ + + ]: 264 : do ii=1,3
3401 [ + + ]: 264 : if(abs(rprimd_ref(ii,ia)) > tol10)then
3402 : 66 : scale_cell(ii) = nint(rprimd_hist(ii,ia) / rprimd_ref(ii,ia))
3403 : : end if
3404 : : end do
3405 : : ! Check if the factor for the supercell is revelant
3406 : 66 : revelant_factor = .TRUE.
3407 [ + + ]: 264 : do ii=1,3
3408 [ + + ]: 198 : if(abs(scale_cell(ii)) < tol10) cycle
3409 : 66 : factor = abs(scale_cell(ii))
3410 [ + + ]: 264 : do jj=ii,3
3411 [ + + ]: 132 : if(abs(scale_cell(jj)) < tol10) cycle
3412 [ + - ]: 396 : if(abs(abs(scale_cell(ii))-abs(scale_cell(jj))) > tol10) revelant_factor = .FALSE.
3413 : : end do
3414 : : end do
3415 [ - + ]: 88 : if(.not.revelant_factor)then
3416 : : write(msg, '(3a)' )&
3417 : 0 : & 'unable to map the hist file ',ch10,&
3418 : 0 : & 'Action: check/change your MD file'
3419 : 0 : ABI_ERROR(msg)
3420 : : else
3421 : 66 : ncell(ia) = int(factor)
3422 : : end if
3423 : : end do
3424 : : end if
3425 : :
3426 [ + + ]: 136 : ncells = product(ncell)
3427 : :
3428 : : !Check if the energy stored in the hist is revelant, sometimes some MD files gives
3429 : : !the energy of the unit cell... This is not suppose to happen... But just in case...
3430 [ + + ]: 705 : do ii=1,nstep_hist
3431 [ + - ]: 705 : if(abs(eff_pot%energy)>tol12)then
3432 : 671 : ratio=hist%etot(ii) / eff_pot%energy
3433 [ + - ]: 671 : if(abs(ratio)<real(huge(factE_hist))*half)then
3434 : 671 : factE_hist = nint(ratio)
3435 [ - + ]: 671 : if(factE_hist == 1) then
3436 : : ! In this case we mutiply the energy of the hist by the number of cell
3437 : 0 : hist%etot(ii) = hist%etot(ii) * ncells
3438 : : end if
3439 [ + + ]: 671 : if(factE_hist /=1 .and. factE_hist /= ncells)then
3440 : 31 : write(msg, '(4a,I0,a,I0,2a,I0,3a,I0,3a)' )ch10,&
3441 : 31 : & ' --- !WARNING',ch10,&
3442 : 31 : & ' The energy of the history step ',ii,' seems to be with multiplicity of ',factE_hist,ch10,&
3443 : 31 : & ' However, the multiplicity of the cell is ',ncells,'.',ch10,&
3444 : 31 : & ' Please check the energy of the step ',ii,ch10,&
3445 : 62 : & ' ---',ch10
3446 [ + - ]: 31 : if(need_verbose) call wrtout(std_out,msg,'COLL')
3447 : : endif
3448 : : else
3449 : 0 : write(msg, '(4a,i0,3a,es16.6,5a)' )ch10,&
3450 : 0 : & ' --- !WARNING',ch10,&
3451 : 0 : & ' The energy of the history step ',ii,' is apparently not initialized.',ch10,&
3452 : 0 : & ' Its current value is',hist%etot(ii),ch10,&
3453 : 0 : & ' This does not allow to perform checking on the multiplicity of the cell ',ch10,&
3454 : 0 : & ' ---',ch10
3455 [ # # ]: 0 : if(need_verbose) call wrtout(std_out,msg,'COLL')
3456 : : end if
3457 : : end if
3458 : : end do
3459 : :
3460 : :
3461 : : !Set the new supercell datatype into the effective potential reference
3462 : 34 : call effective_potential_setSupercell(eff_pot,comm,ncell)
3463 : :
3464 : : !allocation
3465 [ - + + - : 102 : ABI_MALLOC(shift,(3,natom_hist))
+ - + - -
+ - + ]
3466 [ - + + - : 102 : ABI_MALLOC(list_map,(natom_hist))
+ - + - -
+ - + ]
3467 [ - + + - : 102 : ABI_MALLOC(list_reddist,(3,natom_hist))
+ - + - -
+ - + ]
3468 [ - + + - : 68 : ABI_MALLOC(list_absdist,(3,natom_hist))
- + - + ]
3469 [ - + + - : 102 : ABI_MALLOC(list_dist,(natom_hist))
+ - + - -
+ - + ]
3470 [ - + + - : 68 : ABI_MALLOC(xred_ref,(3,natom_hist))
- + - + ]
3471 : :
3472 : : !Putting maping list to zero
3473 [ + + ]: 1534 : list_map = 0
3474 : :
3475 : : !Fill xcart_ref/hist and xred_ref/hist
3476 : :
3477 : : call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3478 : 34 : & eff_pot%supercell%xcart,xred_ref) ! Get xred_ref
3479 : :
3480 : :
3481 [ + + ]: 1534 : do ia=1,natom_hist !Loop over all reference atoms
3482 : : ! Put temporary lists to zero
3483 [ + + + + ]: 310300 : list_reddist = 0
3484 [ + + + + ]: 310300 : list_absdist = 0
3485 [ + + ]: 78700 : list_dist = 0
3486 [ + + + + ]: 310300 : shift = 0
3487 [ + + ]: 78700 : do ib=1,natom_hist !Loop over all atoms of distorted structure
3488 : : !Calculate list of reduced distance between reference atom ia and all others
3489 [ + + ]: 308800 : list_reddist(:,ib) = hist%xred(:,ib,1) - xred_ref(:,ia)
3490 : : !If the distorted atom is further away than half the unit cell shift it.
3491 [ + + ]: 77200 : if(list_reddist(1,ib) > 0.5)then
3492 : 6439 : list_reddist(1,ib) = 1 - list_reddist(1,ib)
3493 : 6439 : shift(1,ib) = -1
3494 : : end if
3495 [ + + ]: 77200 : if(list_reddist(2,ib) > 0.5)then
3496 : 6370 : list_reddist(2,ib) = 1 - list_reddist(2,ib)
3497 : 6370 : shift(2,ib) = -1
3498 : : end if
3499 [ + + ]: 77200 : if(list_reddist(3,ib) > 0.5)then
3500 : 2541 : list_reddist(3,ib) = 1 - list_reddist(3,ib)
3501 : 2541 : shift(3,ib) = -1
3502 : : end if
3503 [ + + ]: 77200 : if(list_reddist(1,ib) < -0.5)then
3504 : 9714 : list_reddist(1,ib) = -1 - list_reddist(1,ib)
3505 : 9714 : shift(1,ib) = 1
3506 : : end if
3507 [ + + ]: 77200 : if(list_reddist(2,ib) < -0.5)then
3508 : 9410 : list_reddist(2,ib) = -1 - list_reddist(2,ib)
3509 : 9410 : shift(2,ib) = 1
3510 : : end if
3511 [ + + ]: 77200 : if(list_reddist(3,ib) < -0.5)then
3512 : 4543 : list_reddist(3,ib) = -1 - list_reddist(3,ib)
3513 : 4543 : shift(3,ib) = 1
3514 : : end if
3515 : 77200 : list_absdist(1,ib) = (rprimd_hist(1,1)+rprimd_hist(2,1)+rprimd_hist(3,1))*list_reddist(1,ib)
3516 : 77200 : list_absdist(2,ib) = (rprimd_hist(1,2)+rprimd_hist(2,2)+rprimd_hist(3,2))*list_reddist(2,ib)
3517 : 77200 : list_absdist(3,ib) = (rprimd_hist(1,3)+rprimd_hist(2,3)+rprimd_hist(3,3))*list_reddist(3,ib)
3518 : 78700 : list_dist(ib) = sqrt(abs(list_absdist(1,ib))**2 + abs(list_absdist(2,ib))**2 + abs(list_absdist(3,ib))**2 )
3519 : : end do !ib
3520 : : !find the closest atom ib
3521 [ + - - + : 78700 : list_map(ia) = minloc(list_dist,DIM=1)
+ + + + ]
3522 : : !If the closest atom ib was shifted, apply and store the shift
3523 [ + + + + : 6018 : if(any(shift(:,list_map(ia)) /= 0))then
+ + ]
3524 [ + + ]: 16 : hist%xred(1,list_map(ia),:)= hist%xred(1,list_map(ia),:) + 1*shift(1,list_map(ia))
3525 [ + + ]: 16 : hist%xred(2,list_map(ia),:)= hist%xred(2,list_map(ia),:) + 1*shift(2,list_map(ia))
3526 [ + + ]: 16 : hist%xred(3,list_map(ia),:)= hist%xred(3,list_map(ia),:) + 1*shift(3,list_map(ia))
3527 : : end if
3528 : : !TEST MS
3529 : : !write(*,*) 'Atom', ia,' of reference is matche with', list_map(ia)
3530 : : !write(*,*) 'xred_ref(',xred_ref(:,ia),'), xred_hist(',hist%(:,list_map(ia)),')'
3531 : : end do ! ia
3532 : :
3533 [ + + ]: 34 : if(need_verbose) then
3534 : 13 : write(msg,'(2a,I3,a,I3,a,I3)') ch10,&
3535 : 26 : & ' The size of the supercell for the fit is ',ncell(1),' ',ncell(2),' ',ncell(3)
3536 : 13 : call wrtout(std_out,msg,'COLL')
3537 : 13 : call wrtout(ab_out,msg,'COLL')
3538 : : end if
3539 : :
3540 : :
3541 [ + + + - : 1534 : if(any(list_map(:)==0))then
- + ]
3542 : : write(msg, '(5a)' )&
3543 : 0 : & 'Unable to map the molecular dynamic file ',ch10,&
3544 : 0 : & 'on the reference supercell structure',ch10,&
3545 : 0 : & 'Action: change the MD file'
3546 : 0 : ABI_ERROR(msg)
3547 : : end if
3548 : :
3549 : 34 : need_map = .FALSE.
3550 [ + + ]: 1534 : do ia=1,natom_hist
3551 [ + + ]: 1534 : if(list_map(ia) /= ia) need_map = .TRUE.
3552 : : end do
3553 [ + + ]: 34 : if(need_map)then
3554 [ + + ]: 21 : if(need_verbose) then
3555 : 1 : write(msg, '(11a)' )ch10,&
3556 : 1 : & ' --- !WARNING',ch10,&
3557 : 1 : & ' The ordering of the atoms in the _HIST.nc file is different,',ch10,&
3558 : 1 : & ' of the one built by multibinit. The _HIST.nc file will be mapped,',ch10,&
3559 : 1 : & ' to the ordering of multibinit.',ch10,&
3560 : 2 : & ' ---',ch10
3561 : 1 : call wrtout(ab_out,msg,'COLL')
3562 : 1 : call wrtout(std_out,msg,'COLL')
3563 : : end if
3564 : :
3565 : : ! Allocate hist datatype
3566 : 21 : call abihist_init(hist_tmp,natom_hist,nstep_hist,.false.,.false.)
3567 : : ! copy all the information
3568 [ + + ]: 512 : do ia=1,nstep_hist
3569 : 491 : hist%ihist = ia
3570 : 491 : hist_tmp%ihist = ia
3571 : 512 : call abihist_copy(hist,hist_tmp)
3572 : : end do
3573 : 21 : hist_tmp%mxhist = nstep_hist
3574 : :
3575 : : ! reoder array
3576 [ + + ]: 861 : do ia=1,natom_hist
3577 [ + + + + ]: 79400 : hist_tmp%xred(:,ia,:) = hist%xred(: ,list_map(ia),:)
3578 [ + + + + ]: 79400 : hist_tmp%fcart(:,ia,:) = hist%fcart(:,list_map(ia),:)
3579 [ + + + + ]: 79421 : hist_tmp%vel(:,ia,:) = hist%vel(:,list_map(ia),:)
3580 : : end do
3581 : :
3582 : : ! free the old hist and reinit
3583 : 21 : call abihist_free(hist)
3584 : 21 : call abihist_init(hist,natom_hist,nstep_hist,.false.,.false.)
3585 : : ! copy the temporary hist into output
3586 [ + + ]: 512 : do ia=1,nstep_hist
3587 : 491 : hist%ihist = ia
3588 : 491 : hist_tmp%ihist = ia
3589 : 512 : call abihist_copy(hist_tmp,hist)
3590 : : end do
3591 : : hist_tmp%mxhist = nstep_hist
3592 : 21 : call abihist_free(hist_tmp)
3593 : :
3594 : : !map also fixes if present
3595 [ + + ]: 21 : if(need_fixmap)then
3596 [ - + + - : 2 : ABI_MALLOC(iatfix_tmp,(3,natom_hist))
- + - + ]
3597 [ + + ]: 41 : do ia=1,natom_hist
3598 [ + + ]: 161 : iatfix_tmp(:,ia) = iatfix(:,list_map(ia))
3599 : : end do
3600 [ + - + - : 162 : iatfix = iatfix_tmp
- + - - +
+ + + ]
3601 [ - + ]: 1 : ABI_FREE(iatfix_tmp)
3602 : : end if
3603 : : end if !need map
3604 : :
3605 : : !deallocation
3606 [ - + ]: 34 : ABI_FREE(shift)
3607 [ - + ]: 34 : ABI_FREE(list_map)
3608 [ - + ]: 34 : ABI_FREE(list_dist)
3609 [ - + ]: 34 : ABI_FREE(list_reddist)
3610 [ - + ]: 34 : ABI_FREE(list_absdist)
3611 [ - + ]: 34 : ABI_FREE(xred_ref)
3612 [ - + - + : 34 : end subroutine effective_potential_file_mapHistToRef
- + - + -
+ - + - +
- + - + -
+ - + ]
3613 : : !!***
3614 : :
3615 : :
3616 : : !****f* m_effective_potential_file/effective_potential_file_readDisplacement
3617 : : !!
3618 : : !! NAME
3619 : : !! effective_potential_file_readDisplacement
3620 : : !!
3621 : : !! FUNCTION
3622 : : !! Read a displacement ASCII file
3623 : : !!
3624 : : !! INPUTS
3625 : : !! filename = path of the file
3626 : : !! natom = number of atoms in the cell
3627 : : !! nstep = number of time step
3628 : : !!
3629 : : !! OUTPUT
3630 : : !! disp(3,natom_sc) = atomics displacement between configuration and the reference
3631 : : !!
3632 : : !! SOURCE
3633 : :
3634 : 0 : subroutine effective_potential_file_readDisplacement(filename,disp,nstep,natom)
3635 : :
3636 : : !Arguments ------------------------------------
3637 : : !scalars
3638 : : integer, intent(in) :: natom,nstep
3639 : : character(len=fnlen),intent(in) :: filename
3640 : : !array
3641 : : real(dp),intent(out) :: disp(nstep,3,natom)
3642 : : !Local variables------------------------------
3643 : : !scalars
3644 : : integer :: ios = 0
3645 : : integer :: ia,istep,mu
3646 : : character(500) :: message
3647 : : character (len=500000) :: line,readline
3648 : : integer :: funit = 666
3649 : : ! *************************************************************************
3650 : :
3651 [ # # ]: 0 : if (open_file(filename,message,unit=funit,form="formatted",&
3652 : : status="old",action="read") /= 0) then
3653 : 0 : ABI_ERROR(message)
3654 : : end if
3655 : :
3656 : 0 : write(message, '(2a)' ) " Read displacements from ", trim(filename)
3657 : :
3658 : 0 : call wrtout(ab_out,message,'COLL')
3659 : 0 : call wrtout(std_out,message,'COLL')
3660 [ # # ]: 0 : do istep=1,nstep
3661 [ # # ]: 0 : do ia=1,natom
3662 : 0 : read(funit,'(a)',iostat=ios) readline
3663 : 0 : line=adjustl(readline)
3664 : 0 : read(unit=line,fmt=*) (disp(istep,mu,ia),mu=1,3)
3665 : 0 : write(111,'(3(es21.12))') disp(istep,:,ia)
3666 : : end do
3667 : : end do
3668 : :
3669 : 0 : close(funit)
3670 : :
3671 : 0 : end subroutine effective_potential_file_readDisplacement
3672 : : !!***
3673 : :
3674 : : !!****f* m_effective_potential_file/elementfromline
3675 : : !! NAME
3676 : : !! elementfromline
3677 : : !!
3678 : : !! FUNCTION
3679 : : !! Read the number of element of a line
3680 : : !!
3681 : : !! INPUTS
3682 : : !! line= string from which the data are read
3683 : : !!
3684 : : !! OUTPUT
3685 : : !! nelement = number of element in the line
3686 : : !!
3687 : : !! SOURCE
3688 : :
3689 : 0 : subroutine elementfromline(line,nelement)
3690 : :
3691 : : !Arguments ---------------------------------------------
3692 : : character(len=*), intent(in) :: line
3693 : : integer, intent(out) :: nelement
3694 : : !Local variables ---------------------------------------
3695 : : integer :: ii,n
3696 : : logical :: element
3697 : : ! *********************************************************************
3698 : :
3699 : : !Set the output
3700 : 0 : nelement = 0
3701 : 0 : n = len_trim(line)
3702 : 0 : element = .false.
3703 [ # # ]: 0 : do ii=1,n
3704 [ # # # # ]: 0 : if(.not.element.and.line(ii:ii) /="") then
3705 : : element=.true.
3706 : : else
3707 [ # # ]: 0 : if((element.and.line(ii:ii) =="")) then
3708 : 0 : element=.false.
3709 : 0 : nelement = nelement + 1
3710 : : end if
3711 : : end if
3712 [ # # ]: 0 : if((element.and.ii==n)) nelement = nelement + 1
3713 : : end do
3714 : :
3715 : 0 : end subroutine elementfromline
3716 : : !!***
3717 : :
3718 : : !!****f* m_effective_potential_file/rdfromline
3719 : : !! NAME
3720 : : !! rdfromline
3721 : : !!
3722 : : !! FUNCTION
3723 : : !! Read the value of a keyword from a XML line
3724 : : !! Same function than m_pawxmlps/paw_rdfromline.F90
3725 : : !!
3726 : : !! INPUTS
3727 : : !! keyword= keyword which value has to be read
3728 : : !! line= string from which the data are read (line from a XML)
3729 : : !!
3730 : : !! OUTPUT
3731 : : !! output= (string) value of the keyword
3732 : : !!
3733 : : !! SOURCE
3734 : :
3735 : 511584 : subroutine rdfromline(keyword,line,output)
3736 : :
3737 : : !Arguments ---------------------------------------------
3738 : : character(len=*), intent(in) :: keyword,line
3739 : : character(len=*), intent(out) :: output
3740 : : !Local variables ---------------------------------------
3741 : 511584 : character(len=len(line)) :: temp
3742 : : integer :: pos,pos2
3743 : :
3744 : : ! *********************************************************************
3745 : :
3746 : 511584 : output=""
3747 : 511584 : pos=index(line,trim(keyword))
3748 [ + + ]: 511584 : if (pos>0) then
3749 [ + - + - ]: 511568 : temp=line(pos+len_trim(keyword):len_trim(line))
3750 : 511568 : pos=index(temp,char(34))
3751 [ + - ]: 511568 : if (pos>0) then
3752 : 511568 : pos2=index(temp(pos+1:len_trim(temp)),char(34))
3753 [ + - ]: 511568 : if (pos2>0) then
3754 [ + - + - ]: 511568 : output=temp(pos+1:pos+pos2-1)
3755 : : end if
3756 : : end if
3757 : : end if
3758 : :
3759 [ + - ]: 511584 : end subroutine rdfromline
3760 : : !!***
3761 : :
3762 : :
3763 : : !!****f* m_effective_potential_file/rmtabfromline
3764 : : !! NAME
3765 : : !! rmtabfromline
3766 : : !!
3767 : : !! FUNCTION
3768 : : !! Read remove tab from the begining of line
3769 : : !!
3770 : : !! INPUTS
3771 : : !! line= string from which the data are read (line from a XML)
3772 : : !!
3773 : : !! OUTPUT
3774 : : !! output= line without tab
3775 : : !!
3776 : : !! SOURCE
3777 : :
3778 : 2062683 : recursive subroutine rmtabfromline(line)
3779 : :
3780 : : !Arguments ---------------------------------------------
3781 : : character(len=*), intent(inout) :: line
3782 : : !Local variables ---------------------------------------
3783 : : integer :: pos
3784 : :
3785 : : ! *********************************************************************
3786 : :
3787 : 2062764 : pos=index(line,char(9))
3788 : 2062764 : if (pos==1) then
3789 [ + - + - ]: 81 : line = line(2:len_trim(line))//" "
3790 : : call rmtabfromline(line)
3791 : : end if
3792 : :
3793 [ + + ]: 2062764 : end subroutine rmtabfromline
3794 : : !!***
3795 : :
3796 : :
3797 : :
3798 : : !!****f* m_effective_potential_file/rdfromline_value
3799 : : !! NAME
3800 : : !! rdfromline
3801 : : !!
3802 : : !! FUNCTION
3803 : : !! Read the value of a keyword from a XML line
3804 : : !!
3805 : : !! INPUTS
3806 : : !! keyword= keyword which value has to be read
3807 : : !! line= string from which the data are read (line from a XML)
3808 : : !!
3809 : : !! OUTPUT
3810 : : !! output= (string) value of the keyword
3811 : : !!
3812 : : !! SOURCE
3813 : :
3814 : 218900 : subroutine rdfromline_value(keyword,line,output)
3815 : :
3816 : : !Arguments ---------------------------------------------
3817 : : character(len=*), intent(in) :: keyword,line
3818 : : character(len=*), intent(out) :: output
3819 : : !Local variables ---------------------------------------
3820 : 218900 : character(len=len(line)) :: temp
3821 : : integer :: pos,pos2
3822 : :
3823 : : ! *********************************************************************
3824 : :
3825 : 218900 : output=""
3826 : 218900 : pos=index(line,trim(keyword))
3827 [ + + ]: 218900 : if (pos==2) then
3828 : 218548 : pos=pos+len_trim(keyword)
3829 : 218548 : pos=pos+index(line(pos:len_trim(line)),char(62))
3830 [ + - + - ]: 218548 : temp=line(pos:len_trim(line))
3831 : 218548 : pos2=index(temp,char(60))
3832 [ + + ]: 218548 : if (pos2>0) then
3833 [ + - + - ]: 215484 : output=line(pos:pos+pos2-2)
3834 : : else
3835 [ + - + - ]: 3064 : output=line(pos:len_trim(line))
3836 : : end if
3837 : : else
3838 [ - + ]: 352 : if(pos>2)then
3839 [ # # # # ]: 0 : output=line(1:pos-3)
3840 : : end if
3841 : : end if
3842 [ + - ]: 218900 : end subroutine rdfromline_value
3843 : : !!***
3844 : :
3845 : :
3846 : : !!****f* m_effpot_xml/char_f2c
3847 : : !! NAME
3848 : : !! char_f_to_c
3849 : : !!
3850 : : !! FUNCTION
3851 : : !! Helper function to convert a Fortran string to a C string
3852 : : !! Based on a routine by Joseph M. Krahn
3853 : : !!
3854 : : !! INPUTS
3855 : : !! f_string=Fortran string
3856 : : !!
3857 : : !! OUTPUT
3858 : : !! c_string=C string
3859 : : !!
3860 : : !! SOURCE
3861 : :
3862 : : #if defined HAVE_XML
3863 : :
3864 : : function char_f2c(f_string) result(c_string)
3865 : :
3866 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_NULL_CHAR
3867 : : !Arguments ------------------------------------
3868 : : character(len=*),intent(in) :: f_string
3869 : : character(kind=C_CHAR,len=1) :: c_string(len_trim(f_string)+1)
3870 : : !Local variables -------------------------------
3871 : : integer :: ii,strlen
3872 : : !! *************************************************************************
3873 : : strlen=len_trim(f_string)
3874 : : forall(ii=1:strlen)
3875 : : c_string(ii)=f_string(ii:ii)
3876 : : end forall
3877 : : c_string(strlen+1)=C_NULL_CHAR
3878 : : end function char_f2c
3879 : : !!***
3880 : :
3881 : : !----------------------------------------------------------------------
3882 : :
3883 : : !!****f* m_effpot_xml/char_c2f
3884 : : !! NAME
3885 : : !! char_c_to_f
3886 : : !!
3887 : : !! FUNCTION
3888 : : !! Helper function to convert a C string to a Fortran string
3889 : : !! Based on a routine by Joseph M. Krahn
3890 : : !!
3891 : : !! INPUTS
3892 : : !! c_string=C string
3893 : : !!
3894 : : !! OUTPUT
3895 : : !! f_string=Fortran string
3896 : : !!
3897 : : !! SOURCE
3898 : :
3899 : : subroutine char_c2f(c_string,f_string)
3900 : :
3901 : : use, intrinsic :: iso_c_binding, only : C_CHAR,C_NULL_CHAR
3902 : : !Arguments ------------------------------------
3903 : : character(kind=C_CHAR,len=1),intent(in) :: c_string(*)
3904 : : character(len=*),intent(out) :: f_string
3905 : : !Local variables -------------------------------
3906 : : integer :: ii
3907 : : !! *************************************************************************
3908 : : ii=1
3909 : : do while(c_string(ii)/=C_NULL_CHAR.and.ii<=len(f_string))
3910 : : f_string(ii:ii)=c_string(ii) ; ii=ii+1
3911 : : end do
3912 : : if (ii<len(f_string)) f_string(ii:)=' '
3913 : : end subroutine char_c2f
3914 : : !!***
3915 : : #endif
3916 : :
3917 : : end module m_effective_potential_file
3918 : : !!***
|
