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