Line data Source code
1 : !!****m* ABINIT/m_effective_potential
2 : !!
3 : !! NAME
4 : !! m_effective_potential
5 : !!
6 : !! FUNCTION
7 : !! Module for the effective potential
8 : !! Container type is defined, and destruction, print subroutines
9 : !! Contain also routine to evaluate the energy,forces and stresses
10 : !!
11 : !! COPYRIGHT
12 : !! Copyright (C) 2010-2026 ABINIT group (AM)
13 : !! This file is distributed under the terms of the
14 : !! GNU General Public Licence, see ~abinit/COPYING
15 : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
17 : !!
18 : !! SOURCE
19 :
20 : #if defined HAVE_CONFIG_H
21 : #include "config.h"
22 : #endif
23 :
24 : #include "abi_common.h"
25 :
26 : module m_effective_potential
27 :
28 : use defs_basis
29 : use defs_abitypes
30 : use m_errors
31 : use m_abicore
32 : use m_strain
33 : use m_ifc
34 : use m_supercell
35 : use m_phonons
36 : use m_ddb
37 : use m_polynomial_conf
38 : use m_polynomial_coeff
39 : use m_anharmonics_terms
40 : use m_harmonics_terms
41 : use m_xmpi
42 : use m_ewald
43 : use m_nctk
44 : use netcdf
45 : #if defined DEV_MS_SCALEUP
46 : use scup_global, only : global_calculate_energy, global_calculate_forces
47 : #endif
48 :
49 : use m_fstrings, only : replace, ftoa, itoa
50 : use m_io_tools, only : open_file, get_unit
51 : use m_dtfil, only : isfile
52 : use m_matrix, only : matr3inv
53 : use m_effpot_mpi, only : effpot_mpi_init,effpot_mpi_type,effpot_mpi_free
54 : use m_abihist, only : abihist
55 : use m_geometry, only : gred2fcart,fcart2gred, xcart2xred, xred2xcart, metric
56 : use m_crystal, only : crystal_t
57 : !use m_anaddb_dataset, only : anaddb_dataset_type, anaddb_dtset_free, outvars_anaddb, invars9
58 :
59 : implicit none
60 :
61 : public :: effective_potential_distributeResidualForces
62 : public :: effective_potential_evaluate
63 : public :: effective_potential_free
64 : public :: effective_potential_freeCoeffs
65 : public :: effective_potential_freempi
66 : public :: effective_potential_generateDipDip
67 : public :: effective_potential_getDisp
68 : public :: effective_potential_init
69 : public :: effective_potential_initmpi
70 : public :: effective_potential_copy
71 : public :: effective_potential_print
72 : public :: effective_potential_printSupercell
73 : public :: effective_potential_setCoeffs
74 : public :: effective_potential_setConfinement
75 : public :: effective_potential_setElastic3rd
76 : public :: effective_potential_setElastic4th
77 : public :: effective_potential_setElasticDispCoupling
78 : public :: effective_potential_setStrainPhononCoupling
79 : public :: effective_potential_setSupercell
80 : public :: effective_potential_writeAbiInput
81 : public :: effective_potential_writeXML
82 : public :: effective_potential_writeAnhHead
83 : !AM_EXPERIMENTAL
84 : public :: effective_potential_computeGradient
85 : ! public :: effective_potential_effpot2ddb
86 : ! public :: effective_potential_printPDOS
87 : public :: effective_potential_checkDEV
88 : public :: effective_potential_writeNETCDF
89 : public :: OPERATOR(==)
90 : !AM_EXPERIMENTAL
91 : !!***
92 :
93 : !!****t* m_effective_potential/effective_potential_type
94 : !! NAME
95 : !! effective_potential_type
96 : !!
97 : !! FUNCTION
98 : !! datatype for a effective potential constructed.
99 : !!
100 : !! SOURCE
101 :
102 : type, public :: effective_potential_type
103 :
104 : character(len=fnlen) :: name
105 : ! Name of the molecule (CaTiO3,...)
106 :
107 : type(crystal_t) :: crystal
108 : ! crystal type
109 : ! contains all information of the crystal
110 :
111 : real(dp):: energy
112 : ! Energy of the system (Hatree)
113 :
114 : real(dp), allocatable :: fcart(:,:)
115 : ! forces(3,natom)
116 : ! initial cartesian forces of the system
117 :
118 : real(dp) :: strten(6)
119 : ! strten(6)
120 : ! initial stresses (Ha/bohr^3) of the system
121 :
122 : type(harmonics_terms_type) :: harmonics_terms
123 : ! type with all information for harmonics terms
124 :
125 : type(anharmonics_terms_type) :: anharmonics_terms
126 : ! type with all information for anharmonics terms
127 :
128 : type(polynomial_conf_type) :: confinement
129 : ! type with all the information for the confinement
130 :
131 : type(supercell_type) :: supercell
132 : ! super cell type
133 : ! Store all the information of the suppercell
134 :
135 : logical :: has_anharmonicsTerms
136 : ! True : the aharmonic part is present
137 :
138 : ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
139 : ! This is for the parallelisation over the supercell
140 : type(effpot_mpi_type) :: mpi_ifc
141 : ! effpot_mpi_type with all the information for the IFC paralellisation
142 :
143 : type(effpot_mpi_type) :: mpi_coeff
144 : ! effpot_mpi_type with all the information for the polynomial coefficients paralellisation
145 : contains
146 : final :: effective_potential_finalizer
147 : end type effective_potential_type
148 : !!***
149 :
150 : interface operator (==)
151 : module procedure effective_potential_compare
152 : end interface
153 :
154 : CONTAINS !===========================================================================================
155 :
156 : !!****f* m_effective_potential/effective_potential_init
157 : !!
158 : !! NAME
159 : !! effective_potential_init
160 : !!
161 : !! FUNCTION
162 : !! Initialize effective_potential datatype
163 : !!
164 : !! INPUTS
165 : !! crytal<type(crystal_t)> = datatype with all the information for the crystal
166 : !! energy = energy of the reference structure
167 : !! ifcs <type(ifc_type)> = ifc type with cell,ewald short and total range of the ifcs
168 : !! ncoeff = number of coefficients in the polynomial
169 : !! nqpt = number of qpoints
170 : !! comm = mpi comunicator
171 : !! coeffs(ncoeff)<type(polynomial_coeff_type)> = optional,list of coefficients for the anharmonic part
172 : !! dynmat(2,3,natom,3,natom,nqpt) = optional,dynamical matrix for each qpoints
173 : !! epsilon_inf(3,3) = optional,dielectric tensor
174 : !! elastic_constants(6,6) = optional,elastic constants tensor
175 : !! elastic3rd(6,6,6) = optional,3 order derivatives with respect to to 3 strain
176 : !! elastic_displacement(6,6,3,natom)=optional, 3 order derivatives with respect to 2 strain and
177 : !! 1 Atom disp
178 : !! fcart(3,natom) = optional,optional,initial fcart in the structure
179 : !! strain_coupling(6,natom,3) = optional, internal strain coupling parameters
180 : !! strten(6) = optional,optional,initial strain in the structure
181 : !! name = optional, name of the structure
182 : !! phonon_strain(6) = optional,ifc type for the phonon-strain coupling (should be in anharmonics_terms)
183 : !! phfrq(3*natom,nqpt) = optional,phonons frequencies for each q points in Hartree/cm
184 : !! qpoints(3,nqpt) = optional,list of qpoints wavevectors
185 : !! has_anharmonicsTerms = optional, (default false) flag to set the anharmonics terms
186 : !! supercell<type(supercell_type)> = optional, supercell type to define
187 : !! zeff(3,natom) = optional,effective charges
188 : !!
189 : !! OUTPUT
190 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype to be initialized
191 : !!
192 : !! SOURCE
193 :
194 76 : subroutine effective_potential_init(crystal,eff_pot,energy,ifcs,ncoeff,nqpt,comm,&
195 152 : & coeffs,dynmat,elastic_constants,elastic3rd,&
196 : & elastic_displacement,epsilon_inf,&
197 114 : & fcart,strain_coupling,strten,name,phonon_strain,&
198 76 : & polynomial_conf,phfrq,qpoints,has_anharmonicsTerms,&
199 76 : & supercell,zeff)
200 :
201 : !Arguments ------------------------------------
202 : !scalars
203 : integer,intent(in) :: comm,ncoeff,nqpt
204 : real(dp),intent(in):: energy
205 : logical,optional,intent(in) :: has_anharmonicsTerms
206 : !arrays
207 : type(crystal_t),intent(in) :: crystal
208 : type(effective_potential_type), intent(out) :: eff_pot
209 : type(ifc_type),intent(in) :: ifcs
210 : character(len=fnlen), optional,intent(in) :: name
211 : real(dp),optional,intent(in) :: epsilon_inf(3,3),elastic_constants(6,6)
212 : real(dp),optional,intent(in) :: dynmat(:,:,:,:,:,:),qpoints(:,:),phfrq(:,:)
213 : real(dp),optional,intent(in) :: strain_coupling(:,:,:),zeff(:,:,:)
214 : type(supercell_type),optional,intent(in) :: supercell
215 : type(ifc_type),optional,intent(in) :: phonon_strain(6)
216 : real(dp),optional,intent(in) :: elastic3rd(6,6,6)
217 : real(dp),optional,intent(in) :: elastic_displacement(6,6,3,crystal%natom)
218 : type(polynomial_coeff_type),optional,intent(in) :: coeffs(:)
219 : real(dp),optional,intent(in) :: strten(6)
220 : real(dp),optional,intent(in) :: fcart(3,crystal%natom)
221 : type(polynomial_conf_type),optional,intent(in) :: polynomial_conf
222 :
223 : !Local variables-------------------------------
224 : !scalar
225 : character(len=500) :: msg
226 : !arrays
227 :
228 : ! *************************************************************************
229 :
230 : !1-Free the effective potential before filling it
231 76 : call effective_potential_free(eff_pot)
232 :
233 76 : if (present(name)) then
234 38 : eff_pot%name = name
235 : else
236 38 : eff_pot%name = ''
237 : end if
238 :
239 : !2-Perform some checks
240 76 : if (crystal%natom < 1) then
241 : write(msg, '(a,a,a,i10,a)' )&
242 0 : & 'The cell must have at least one atom.',ch10,&
243 0 : & 'The number of atom is ',crystal%natom,'.'
244 0 : ABI_BUG(msg)
245 : end if
246 :
247 76 : if (crystal%ntypat < 1) then
248 : write(msg, '(a,a,a,i10,a)' )&
249 0 : & 'The cell must have at least one type of atom.',ch10,&
250 0 : & 'The number of type of atom is ',crystal%ntypat,'.'
251 0 : ABI_BUG(msg)
252 : end if
253 :
254 : !3-Fill energy of the crystal (hartree)
255 76 : eff_pot%energy = energy
256 :
257 : !1-Fill the crystal
258 : !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
259 : call eff_pot%crystal%init(crystal%amu,crystal%space_group,crystal%natom,&
260 : & crystal%npsp,crystal%ntypat,crystal%nsym,crystal%rprimd,&
261 : & crystal%typat,crystal%xred,crystal%zion,crystal%znucl,&
262 : & crystal%timrev,.FALSE.,.FALSE.,crystal%title,&
263 76 : & symrel=crystal%symrel,tnons=crystal%tnons,symafm=crystal%symafm)
264 :
265 : !4-Fill harmonic part
266 76 : call harmonics_terms_init(eff_pot%harmonics_terms,ifcs,crystal%natom,ifcs%nrpt)
267 :
268 : !5-Init the anharmonics_terms to set the flag to false
269 76 : call anharmonics_terms_init(eff_pot%anharmonics_terms,crystal%natom,ncoeff)
270 :
271 : !5-Fill optional inputs
272 228 : ABI_MALLOC(eff_pot%fcart,(3,eff_pot%crystal%natom))
273 1652 : eff_pot%fcart = zero
274 76 : if(present(fcart))then
275 836 : eff_pot%fcart = fcart
276 : end if
277 :
278 76 : if(present(elastic_constants))then
279 3268 : eff_pot%harmonics_terms%elastic_constants(:,:) = elastic_constants
280 : end if
281 :
282 76 : if(present(epsilon_inf))then
283 988 : eff_pot%harmonics_terms%epsilon_inf(:,:) = epsilon_inf(:,:)
284 : end if
285 :
286 76 : if(present(dynmat).and.present(qpoints).and.present(phfrq))then
287 76 : call harmonics_terms_setDynmat(dynmat,eff_pot%harmonics_terms,crystal%natom,nqpt,phfrq,qpoints)
288 : end if
289 :
290 76 : if(present(strain_coupling))then
291 76 : call harmonics_terms_setInternalStrain(eff_pot%harmonics_terms,crystal%natom,strain_coupling)
292 : end if
293 :
294 76 : if(present(zeff))then
295 76 : call harmonics_terms_setEffectiveCharges(eff_pot%harmonics_terms,crystal%natom,zeff)
296 : end if
297 :
298 532 : eff_pot%strten = zero
299 76 : if(present(strten))then
300 266 : eff_pot%strten(:) = strten(:)
301 : end if
302 :
303 : !Set the flag for the strain coupling
304 76 : if(present(has_anharmonicsTerms)) then
305 76 : eff_pot%has_anharmonicsTerms = has_anharmonicsTerms
306 : else
307 0 : eff_pot%has_anharmonicsTerms = .false.
308 : end if
309 :
310 : !Allocation of phonon strain coupling array (3rd order)
311 76 : if(present(phonon_strain).and.has_anharmonicsTerms) then
312 : call anharmonics_terms_setStrainPhononCoupling(eff_pot%anharmonics_terms,crystal%natom,&
313 45 : & phonon_strain)
314 : end if
315 :
316 : !Set the 3rd order elastic tensor
317 76 : if(present(elastic3rd).and.has_anharmonicsTerms)then
318 45 : call anharmonics_terms_setElastic3rd(eff_pot%anharmonics_terms,elastic3rd)
319 : end if
320 :
321 : !TODO Comment Marcus: Uncomment and implement below if you want to use the elastic4th implementation, which is there, but
322 : ! not used anywere
323 : ! Below I just ensured that the has_elastic4th variable is set to FALSE as it was causing problems on some builders
324 : ! in the test farm
325 : ! Set the 4th order elastic tensor
326 : ! if(present(elastic4th).and.has_anharmonicsTerms)then
327 : ! call anharmonics_terms_setElastic3rd(eff_pot%anharmonics_terms,elastic4th)
328 : ! end if
329 :
330 : !MS Ensure has_elastic4th is .FALSE.
331 : ! eff_pot%anharmonics_terms%has_elastic4th=.FALSE.
332 :
333 :
334 : !Allocation of 3rd order with respecto to 2 strain and 1 atomic displacement
335 76 : if(present(elastic_displacement).and.has_anharmonicsTerms)then
336 : call anharmonics_terms_setElasticDispCoupling(eff_pot%anharmonics_terms,crystal%natom,&
337 38 : & elastic_displacement)
338 : end if
339 :
340 : !Allocation of the coefficients
341 76 : if(present(coeffs))then
342 27 : if(ncoeff /= size(coeffs))then
343 0 : ABI_BUG('ncoeff has not the same size than coeffs array')
344 : end if
345 27 : call effective_potential_setCoeffs(coeffs,eff_pot,ncoeff)
346 : end if
347 :
348 76 : if(present(supercell))then
349 38 : call effective_potential_setSupercell(eff_pot,comm,supercell=supercell)
350 : else
351 : ! call init_supercell(eff_pot%crystal%natom, (/1,0,0, 0,1,0, 0,0,1/), eff_pot%crystal%rprimd,&
352 : !& eff_pot%crystal%typat, eff_pot%crystal%xcart, eff_pot%crystal%znucl, eff_pot%supercell)
353 38 : call effective_potential_setSupercell(eff_pot,comm,ncell=(/1,1,1/))
354 : end if
355 :
356 : !Set the confinement potential
357 76 : if(present(polynomial_conf)) then
358 : call effective_potential_setConfinement(polynomial_conf%cutoff_disp,polynomial_conf%cutoff_strain,&
359 : & eff_pot,polynomial_conf%factor_disp,&
360 : & polynomial_conf%factor_strain,polynomial_conf%ndisp,&
361 : & polynomial_conf%power_disp,polynomial_conf%power_strain,&
362 0 : & polynomial_conf%need_confinement)
363 : end if
364 :
365 76 : end subroutine effective_potential_init
366 : !!***
367 :
368 :
369 : !!****f* m_effective_potential/effective_potential_initmpi
370 : !! NAME
371 : !! effective_potential_initmpi
372 : !!
373 : !! FUNCTION
374 : !! Initializes the mpi information for parallelism over supercell.
375 : !! Only the parallelisation over cell is done here.
376 : !! The parallelisation over cell and coeff is disable for now (experimental)
377 : !!
378 : !! INPUTS
379 : !! eff_pot<type(effective_potential_type)> = datatype for the effective potential
380 : !! comm = MPI communicator
381 : !!
382 : !! OUTPUT
383 : !! This is for the parallelisation over the supercell
384 : !! eff_pot%mpi_ifc%me_supercell = Index of my processor in the comm. over one cell
385 : !! eff_pot%mpi_ifc%my_ncell = Number of cell treated by current proc
386 : !! eff_pot%mpi_ifc%my_cells(:) = Number of the cells in the supercell treat by this CPU
387 : !! eff_pot%mpi_ifc%my_index_cells(:,:) = indexes of the cells in the supercell treat by this CPU
388 : !!
389 : !! SOURCE
390 :
391 221 : subroutine effective_potential_initmpi(eff_pot,comm)
392 :
393 : !Arguments ------------------------------------
394 : !scalars
395 : type(effective_potential_type),intent(inout) :: eff_pot
396 : integer,intent(in) :: comm
397 : !arrays
398 :
399 : !Local variables-------------------------------
400 : !scalars
401 : integer :: ndiv
402 : integer :: ncell
403 : !array
404 : integer :: cell_number(3)
405 : !character(len=500) :: msg
406 : ! ***********************************************************************
407 :
408 : !Set the number of cell in the supercell
409 221 : cell_number(1) = eff_pot%supercell%rlatt(1,1)
410 221 : cell_number(2) = eff_pot%supercell%rlatt(2,2)
411 221 : cell_number(3) = eff_pot%supercell%rlatt(3,3)
412 884 : ncell = product(cell_number(:))
413 :
414 : !Do some checks
415 884 : if (any(cell_number <= 0).or.ncell<=0) then
416 0 : ABI_ERROR('No supercell found for setting')
417 : end if
418 :
419 : !First mpi_ifc
420 221 : ndiv = 1
421 221 : call effpot_mpi_free(eff_pot%mpi_ifc)
422 : call effpot_mpi_init(eff_pot%harmonics_terms%ifcs%cell,cell_number,eff_pot%mpi_ifc,&
423 221 : & eff_pot%crystal%natom,ndiv,eff_pot%harmonics_terms%ifcs%nrpt,comm)
424 :
425 : !Second mpi_coeff
426 : ndiv = 1
427 221 : call effpot_mpi_free(eff_pot%mpi_coeff)
428 : call effpot_mpi_init(eff_pot%harmonics_terms%ifcs%cell,cell_number,eff_pot%mpi_coeff,&
429 221 : & eff_pot%crystal%natom,ndiv,eff_pot%harmonics_terms%ifcs%nrpt,comm)
430 :
431 221 : end subroutine effective_potential_initmpi
432 : !!***
433 :
434 : !!****f* m_effective_potential/effective_potential_copy
435 : !!
436 : !! NAME
437 : !! effective_potential_copy
438 : !!
439 : !! FUNCTION
440 : !! Copy one effective potential to another i.e. initialize the output effective potentiale
441 : !! eff_pot_out with the values the input effective potential eff_pot_in
442 : !!
443 : !! INPUTS
444 : !!
445 : !! eff_pot_in<type(effective_potential_type)> = effective_potential datatype to be initialized
446 : !!
447 : !! OUTPUT
448 : !! eff_pot_out<type(effective_potential_type)> = effective_potential datatype to be initialized
449 : !!
450 : !! SOURCE
451 :
452 38 : subroutine effective_potential_copy(eff_pot_out,eff_pot_in,comm)
453 :
454 : !Arguments ------------------------------------
455 : !scalars
456 : type(effective_potential_type),intent(out) :: eff_pot_out
457 : type(effective_potential_type),intent(in) :: eff_pot_in
458 : integer,intent(in) :: comm
459 : ! ***********************************************************************
460 :
461 : call effective_potential_init(eff_pot_in%crystal,eff_pot_out,eff_pot_in%energy,eff_pot_in%harmonics_terms%ifcs,&
462 : & eff_pot_in%anharmonics_terms%ncoeff,eff_pot_in%harmonics_terms%nqpt,comm,&
463 : & coeffs=eff_pot_in%anharmonics_terms%coefficients,dynmat=eff_pot_in%harmonics_terms%dynmat,&
464 : & elastic_constants=eff_pot_in%harmonics_terms%elastic_constants,&
465 : & elastic3rd=eff_pot_in%anharmonics_terms%elastic3rd,&
466 : & elastic_displacement=eff_pot_in%anharmonics_terms%elastic_displacement,&
467 : & epsilon_inf=eff_pot_in%harmonics_terms%epsilon_inf,&
468 : & fcart=eff_pot_in%fcart,strain_coupling=eff_pot_in%harmonics_terms%strain_coupling,&
469 : & strten=eff_pot_in%strten,name=eff_pot_in%name,&
470 : & phonon_strain=eff_pot_in%anharmonics_terms%phonon_strain,phfrq=eff_pot_in%harmonics_terms%phfrq,&
471 : & qpoints=eff_pot_in%harmonics_terms%qpoints,has_anharmonicsTerms=eff_pot_in%has_anharmonicsTerms,&
472 : & supercell=eff_pot_in%supercell,&
473 38 : & zeff=eff_pot_in%harmonics_terms%zeff)
474 :
475 :
476 38 : end subroutine effective_potential_copy
477 : !!***
478 :
479 : !****f* m_effective_potential/effective_potential_free
480 : !!
481 : !! NAME
482 : !! effective_potential_free
483 : !!
484 : !! FUNCTION
485 : !! deallocate all dynamic memory for this effective potential datatype
486 : !!
487 : !! INPUTS
488 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
489 : !!
490 : !! OUTPUT
491 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
492 : !!
493 : !! SOURCE
494 :
495 534 : subroutine effective_potential_free(eff_pot)
496 :
497 : !Arguments ------------------------------------
498 : !scalars
499 : !array
500 : type(effective_potential_type), intent(inout) :: eff_pot
501 : !Local variables-------------------------------
502 : !scalars
503 : !array
504 :
505 : ! *************************************************************************
506 :
507 534 : eff_pot%name = ''
508 534 : eff_pot%energy = zero
509 3738 : eff_pot%strten = zero
510 534 : eff_pot%has_anharmonicsTerms = .FALSE.
511 534 : eff_pot%anharmonics_terms%bounded = .FALSE.
512 :
513 534 : if(allocated(eff_pot%fcart)) then
514 2289 : eff_pot%fcart=zero
515 105 : ABI_FREE(eff_pot%fcart)
516 : end if
517 :
518 : ! Free others datatypes
519 534 : call anharmonics_terms_free(eff_pot%anharmonics_terms)
520 534 : call harmonics_terms_free(eff_pot%harmonics_terms)
521 534 : call eff_pot%supercell%free()
522 534 : call eff_pot%crystal%free()
523 534 : call effective_potential_freempi(eff_pot)
524 534 : call polynomial_conf_free(eff_pot%confinement)
525 :
526 534 : end subroutine effective_potential_free
527 : !!***
528 :
529 : !!****f* m_effective_potential/effective_potential_finalizer
530 : !!
531 : !! NAME
532 : !! effective_potential_finalizer
533 : !!
534 : !! FUNCTION
535 : !! Finalizer procedure for effective_potential_type to automatically free memory
536 : !!
537 : !! SOURCE
538 :
539 286 : subroutine effective_potential_finalizer(this)
540 : type(effective_potential_type), intent(inout) :: this
541 0 : call effective_potential_free(this)
542 286 : end subroutine effective_potential_finalizer
543 : !!***
544 :
545 : !****f* m_effective_potential/effective_potential_freeCoeffs
546 : !!
547 : !! NAME
548 : !! effective_potential_freeCoeffs
549 : !!
550 : !! FUNCTION
551 : !! deallocate all dynamic memory for the coefficients of
552 : !! the anharmonic part of the effective potential datatype
553 : !!
554 : !! INPUTS
555 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
556 : !!
557 : !! OUTPUT
558 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
559 : !!
560 : !! SOURCE
561 :
562 76 : subroutine effective_potential_freeCoeffs(eff_pot)
563 :
564 : !Arguments ------------------------------------
565 : !scalars
566 : !array
567 : type(effective_potential_type), intent(inout) :: eff_pot
568 : !Local variables-------------------------------
569 : !scalars
570 : !array
571 : ! *************************************************************************
572 :
573 76 : call anharmonics_terms_freeCoeffs(eff_pot%anharmonics_terms)
574 :
575 76 : end subroutine effective_potential_freeCoeffs
576 : !!***
577 :
578 : !****f* m_effective_potential/effective_potential_freempi
579 : !!
580 : !! NAME
581 : !! effective_potential_freempi
582 : !!
583 : !! FUNCTION
584 : !! deallocate all dynamic memory for mpi of supercell
585 : !!
586 : !! INPUTS
587 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
588 : !!
589 : !! OUTPUT
590 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
591 : !!
592 : !! SOURCE
593 :
594 534 : subroutine effective_potential_freempi(eff_pot)
595 :
596 : !Arguments ------------------------------------
597 : !scalars
598 : !array
599 : type(effective_potential_type), intent(inout) :: eff_pot
600 :
601 : !Local variables-------------------------------
602 : !scalars
603 : !array
604 :
605 : ! *************************************************************************
606 534 : call effpot_mpi_free(eff_pot%mpi_ifc)
607 534 : call effpot_mpi_free(eff_pot%mpi_coeff)
608 :
609 0 : end subroutine effective_potential_freempi
610 : !!***
611 :
612 : !****f* m_effective_potential/effective_potential_generateDipDip
613 : !!
614 : !! NAME
615 : !! effective_potential_generateDipDip
616 : !!
617 : !! FUNCTION
618 : !! Generate the supercell of the structure inside of the effective
619 : !! potential and fill the supercell datatype. Also adapt the harmonic
620 : !! part for the supercell (compute dipole-dipole interation)
621 : !!
622 : !! INPUTS
623 : !! eff_pot = effective potential datatype
624 : !! option = 0 Just generate supercell and fill effective potential
625 : !! 1 Regenerate the dipole dipo;le interaction
626 : !! ncell(3) = number of cell in the direction x, y and z
627 : !! comm=MPI communicator
628 : !!
629 : !! OUTPUT
630 : !! eff_pot<type(effective_potential_type)> = effective potential datatype
631 : !!
632 : !! SOURCE
633 :
634 67 : subroutine effective_potential_generateDipDip(eff_pot,ncell,option,asr,comm,file_option)
635 :
636 : !Arguments ------------------------------------
637 : !scalars
638 : integer,intent(in) :: option,asr
639 : integer,intent(in) :: comm
640 : !array
641 : integer,intent(in) :: ncell(3)
642 : type(effective_potential_type),intent(inout) :: eff_pot
643 : integer,optional,intent(in) :: file_option
644 : !Local variables-------------------------------
645 : !scalar
646 : integer,parameter :: master=0
647 : integer :: first_coordinate
648 : integer :: ia,i1,i2,i3,ii,ierr,irpt,irpt2,irpt_ref,min1,min2,min3
649 : integer :: min1f, min2f, min3f, max1f, max2f, max3f
650 : integer :: min1_cell,min2_cell,min3_cell,max1_cell,max2_cell,max3_cell
651 : integer :: max1,max2,max3,my_rank,natom_uc
652 : integer :: nproc,second_coordinate,size_tmp,sumg0
653 : integer :: my_nrpt,nrpt_alone
654 : real(dp) :: ucvol, eta
655 : character(len=500) :: msg
656 : logical :: iam_master , has_totFC
657 : !array
658 67 : integer,allocatable :: my_index_rpt(:,:)
659 67 : integer,allocatable :: bufsize(:),bufdisp(:)
660 67 : integer,allocatable :: my_irpt(:)
661 : real(dp) :: acell(3)
662 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
663 67 : real(dp),allocatable :: buff_ewald(:,:,:,:,:,:),dyew(:,:,:,:,:), dyewq0(:,:,:,:,:)
664 67 : real(dp),allocatable :: xred(:,:),xred_tmp(:,:),zeff_tmp(:,:,:),qdrp_cart(:,:,:,:)
665 67 : type(supercell_type) :: supercell
666 469 : type(ifc_type) :: ifc_tmp
667 67 : integer, allocatable :: full_cell(:,:)
668 : integer :: full_nrpt
669 67 : real(dp), allocatable :: full_cell_atmfrc(:,:,:,:,:), full_cell_short_atmfrc(:,:,:,:,:), full_cell_ewald_atmfrc(:,:,:,:,:)
670 :
671 : integer :: in_file_option
672 : ! *************************************************************************
673 :
674 : !0 MPI variables
675 67 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
676 : iam_master=.FALSE.
677 67 : iam_master = (my_rank == master)
678 67 : ierr=0
679 :
680 67 : in_file_option = 0
681 67 : if (present(file_option)) then
682 38 : in_file_option = file_option
683 : end if
684 : !0 Check the size of the cell
685 268 : do ia=1,3
686 268 : if(ncell(ia)<0.or.ncell(ia)>150)then
687 : write(msg, '(a,i0,a,i0,a,a,a,i0,a)' )&
688 0 : & 'ncell(',ia,') is ',ncell(ia),', which is lower than 0 of superior than 150.',&
689 0 : & ch10,'Action: correct ncell(',ia,').'
690 0 : ABI_ERROR(msg)
691 : end if
692 : end do
693 :
694 : call supercell%init(eff_pot%crystal%natom, &
695 : & (/ncell(1),0,0, 0,ncell(2),0, 0,0,ncell(3)/),&
696 : & eff_pot%crystal%rprimd,&
697 : & eff_pot%crystal%typat,&
698 : & eff_pot%crystal%xcart,&
699 670 : & eff_pot%crystal%znucl)
700 :
701 : !set variables
702 67 : natom_uc = eff_pot%crystal%natom
703 :
704 : !1 Store the information of the supercell of the reference structure into effective potential
705 67 : call supercell%copy(eff_pot%supercell)
706 : !2 Initialisation of new mpi over supercell
707 67 : call effective_potential_initmpi(eff_pot,comm)
708 :
709 : !3 Check if the bound of new cell correspond to the effective potential
710 : !only for option=zero
711 : !set min and max
712 : ! NGQPT
713 6342 : min1 = minval(eff_pot%harmonics_terms%ifcs%cell(1,:)) !negative maximum cells in direc 1 MARCUS
714 6342 : min2 = minval(eff_pot%harmonics_terms%ifcs%cell(2,:))
715 6342 : min3 = minval(eff_pot%harmonics_terms%ifcs%cell(3,:))
716 6342 : max1 = maxval(eff_pot%harmonics_terms%ifcs%cell(1,:)) !maximum cells in direc 1 MARCUS
717 6342 : max2 = maxval(eff_pot%harmonics_terms%ifcs%cell(2,:))
718 6342 : max3 = maxval(eff_pot%harmonics_terms%ifcs%cell(3,:))
719 67 : write(msg,'(5a,2I3,a,2I3,a,2I3,a)') ch10,' Bound for ifc SR:',&
720 134 : & ch10,ch10, " x=[",min1,max1,"], y=[",min2,max2,"] and z=[",min3,max3,"]"
721 67 : call wrtout(ab_out,msg,'COLL')
722 67 : call wrtout(std_out,msg,'COLL')
723 :
724 67 : if(option==0) then
725 28 : if(((max1-min1+1)/=ncell(1).and.&
726 : & (max2-min2+1)/=ncell(2).and.(max3-min3+1)/=ncell(3))) then
727 1640 : write(msg, '(90a,3I3,5a,3I3,3a)' )ch10,('-',i1=1,80),ch10,ch10,&
728 20 : & ' --- !WARNING:',ch10,&
729 20 : & ' dipdip is set to zero, the longe range interation might be wrong',ch10,&
730 20 : & ' because it is not recompute.',ch10,&
731 20 : & ' The previous harmonic part is build for ',(max1-min1+1),(max2-min2+1),(max3-min3+1)&
732 20 : &, ' cell.',ch10,&
733 20 : & ' Be sure than the dipole-dipole interation is correct ',ch10,&
734 20 : & ' for the supercell: ',ncell(:),' or set dipdip to 1',ch10,&
735 1660 : & ' ---'
736 20 : call wrtout(std_out,msg,"COLL")
737 : else
738 656 : write(msg,'(89a)')ch10,('-',i1=1,80),ch10,ch10,&
739 8 : & ' --- !WARNING:',ch10,&
740 8 : & ' dipdip is set to zero, the longe range interation is not recompute.',ch10,&
741 664 : & ' ---',ch10
742 8 : call wrtout(std_out,msg,"COLL")
743 : end if
744 :
745 :
746 2268 : write(msg,'(a,(80a))') ch10,('=',i1=1,80)
747 28 : call wrtout(ab_out,msg,'COLL')
748 28 : call wrtout(std_out,msg,'COLL')
749 :
750 : !4-Adapt harmonic part
751 156 : else if (option>=1.and.all(ncell(:)>0)) then
752 :
753 3159 : write(msg,'(a,(80a),3a)') ch10,('=',i1=1,80),ch10,' Generation of new ifc',ch10
754 39 : call wrtout(ab_out,msg,'COLL')
755 39 : call wrtout(std_out,msg,'COLL')
756 :
757 39 : irpt_ref = 0
758 39 : irpt = 0
759 39 : min1_cell = 0; max1_cell = 0
760 39 : min2_cell = 0; max2_cell = 0
761 39 : min3_cell = 0; max3_cell = 0
762 :
763 : ! ncell
764 39 : call findBound_supercell(min1_cell,max1_cell,ncell(1))
765 39 : call findBound_supercell(min2_cell,max2_cell,ncell(2))
766 39 : call findBound_supercell(min3_cell,max3_cell,ncell(3))
767 :
768 : write(msg, '(2a)' )&
769 39 : & ' dipdip is set to one, the dipole-dipole interation is recompute.'
770 39 : call wrtout(ab_out,msg,'COLL')
771 39 : call wrtout(std_out,msg,'COLL')
772 :
773 : !! Generate new bound
774 : ! !if(option==1)then
775 : ! ! If NGQPT > NCELL "=" dipdip range
776 : ! !if ((abs(min1) > abs(min1_cell)).or.(abs(max1) > abs(max1_cell)).or.&
777 : !& ! (abs(min2) > abs(min2_cell)).or.(abs(max2) > abs(max2_cell)).or.&
778 : !& ! (abs(min3) > abs(min3_cell)).or.(abs(max3) > abs(max3_cell))) then
779 : ! ! write(msg, '(6a,3I4,3a)' )ch10,&
780 : !& ! ' --- !WARNING',ch10,&
781 : !& ! ' The range of the dipole-dipole interaction is the same than the short-range.',ch10,&
782 : !& ! ' So the range of the total ifc is ',int((/(max1-min1+1),(max2-min2+1),(max3-min3+1)/),dp),' cell',ch10,&
783 : !& ! ' ---'
784 : ! ! call wrtout(std_out,msg,"COLL")
785 : ! ! if (abs(min1) < abs(min1_cell)) min1 = min1_cell
786 : ! ! if (abs(min2) < abs(min2_cell)) min2 = min2_cell
787 : ! ! if (abs(min3) < abs(min3_cell)) min3 = min3_cell
788 : ! ! if (abs(max1) < abs(max1_cell)) max1 = max1_cell
789 : ! ! if (abs(max2) < abs(max2_cell)) max2 = max2_cell
790 : ! ! if (abs(max3) < abs(max3_cell)) max3 = max3_cell
791 : !! ! If the cell is smaller, we redifine new cell to take into acount all atoms
792 : ! ! call supercell%free()
793 : ! ! call init_supercell(natom_uc,(/(max1-min1+1),0,0, 0,(max2-min2+1),0, 0,0,(max3-min3+1)/),&
794 : !& ! eff_pot%crystal%rprimd,eff_pot%crystal%typat,&
795 : !& ! eff_pot%crystal%xcart,eff_pot%crystal%znucl,supercell)
796 : !
797 : !! ! Store the information of the supercell of the reference structure into effective potential
798 : ! ! call effective_potential_setSupercell(eff_pot,comm,supercell=supercell)
799 : ! !else
800 : ! ! min1 = min1_cell ; min2 = min2_cell ; min3 = min3_cell
801 : ! ! max1 = max1_cell ; max2 = max2_cell ; max3 = max3_cell
802 : ! !end if
803 : ! !end if
804 :
805 : ! Print the new boundary
806 39 : write(msg,'(5a,2I3,a,2I3,a,2I3,4a)') ch10,' Bound for ifc (LR):',&
807 39 : & ch10,ch10, " x=[",min1_cell,max1_cell,"], y=[",min2_cell,max2_cell,"] and z=[",min3_cell,max3_cell,"]",ch10,ch10,&
808 78 : & " Computation of new dipole-dipole interaction."
809 39 : call wrtout(ab_out,msg,'COLL')
810 39 : call wrtout(std_out,msg,'COLL')
811 :
812 : ! Count the new number of ifc
813 128 : do i1=min1_cell,max1_cell
814 379 : do i2=min2_cell,max2_cell
815 907 : do i3=min3_cell,max3_cell
816 567 : irpt = irpt +1
817 818 : if(i1==0.and.i2==0.and.i3==0) irpt_ref = irpt
818 : end do
819 : end do
820 : end do
821 :
822 39 : ifc_tmp%nrpt = irpt
823 117 : ABI_MALLOC(ifc_tmp%cell,(3,ifc_tmp%nrpt))
824 2307 : ifc_tmp%cell(:,:) = 0
825 :
826 : ! Set MPI here and not at the begining because the number of cell is adjust just before
827 : ! Here we store in my_irpt a list with the number of each cell to be treat by this CPU
828 : ! Determine the number of cell for each CPU
829 117 : ABI_MALLOC(bufsize,(nproc))
830 78 : ABI_MALLOC(bufdisp,(nproc))
831 :
832 39 : nrpt_alone = mod(ifc_tmp%nrpt,nproc)
833 39 : my_nrpt = int(real(ifc_tmp%nrpt,sp)/nproc)
834 39 : if(my_rank >= (nproc-nrpt_alone)) then
835 2 : my_nrpt = my_nrpt + 1
836 : end if
837 :
838 : ! Initialisation of ifc temporary
839 234 : ABI_MALLOC(buff_ewald,(2,3,natom_uc,3,natom_uc,my_nrpt))
840 234 : ABI_MALLOC(ifc_tmp%short_atmfrc,(3,natom_uc,3,natom_uc,ifc_tmp%nrpt))
841 156 : ABI_MALLOC(ifc_tmp%ewald_atmfrc,(3,natom_uc,3,natom_uc,ifc_tmp%nrpt))
842 156 : ABI_MALLOC(ifc_tmp%atmfrc,(3,natom_uc,3,natom_uc,ifc_tmp%nrpt))
843 117 : ABI_MALLOC(my_irpt,(my_nrpt))
844 117 : ABI_MALLOC(my_index_rpt,(3,my_nrpt))
845 :
846 538 : my_irpt = 0
847 2035 : my_index_rpt(:,:) = 0
848 459550 : ifc_tmp%atmfrc(:,:,:,:,:) = zero
849 459550 : ifc_tmp%short_atmfrc(:,:,:,:,:) = zero
850 459550 : ifc_tmp%ewald_atmfrc(:,:,:,:,:) = zero
851 1091072 : buff_ewald(:,:,:,:,:,:) = zero
852 :
853 : ! Allocation of array
854 538 : do irpt = 1,my_nrpt
855 538 : if(my_rank >= (nproc-nrpt_alone))then
856 : my_irpt(irpt)=(int(real(ifc_tmp%nrpt,sp)/nproc))*(my_rank)+&
857 2 : & (my_rank - (nproc-nrpt_alone)) + irpt
858 : else
859 497 : my_irpt(irpt)=(my_nrpt)*(my_rank) + irpt
860 : end if
861 : end do
862 :
863 : irpt = 0
864 : ii = 0
865 128 : do i1=min1_cell,max1_cell
866 379 : do i2=min2_cell,max2_cell
867 907 : do i3=min3_cell,max3_cell
868 567 : ii = ii +1
869 567 : ifc_tmp%cell(1,ii) = i1; ifc_tmp%cell(2,ii) = i2; ifc_tmp%cell(3,ii) = i3;
870 27258 : if(any(my_irpt==ii))then
871 499 : irpt=irpt+1
872 499 : my_index_rpt(1,irpt) = i1;
873 499 : my_index_rpt(2,irpt) = i2;
874 499 : my_index_rpt(3,irpt) = i3;
875 : end if
876 : end do
877 : end do
878 : end do
879 :
880 : ! Allocate and initialize some array
881 117 : ABI_MALLOC(xred_tmp,(3,2*natom_uc))
882 117 : ABI_MALLOC(xred,(3,supercell%natom))
883 117 : ABI_MALLOC(zeff_tmp,(3,3,2*natom_uc))
884 117 : ABI_MALLOC(qdrp_cart,(3,3,3,2*natom_uc))
885 195 : ABI_MALLOC(dyew,(2,3,2*natom_uc,3,2*natom_uc))
886 156 : ABI_MALLOC(dyewq0,(2,3,natom_uc,3,natom_uc))
887 :
888 146311 : dyew = zero
889 36995 : dyewq0 = zero
890 14707 : xred(:,:) = zero
891 1591 : xred_tmp(:,:) = zero
892 5083 : zeff_tmp(:,:,:) = zero
893 15559 : qdrp_cart = zero
894 39 : sumg0 = 0
895 156 : acell = one
896 :
897 39 : call matr3inv(supercell%rprimd,gprimd)
898 : call xcart2xred(supercell%natom,supercell%rprimd,&
899 39 : & supercell%xcart,xred)
900 39 : call metric(gmet,gprimd,-1,rmet,supercell%rprimd,ucvol)
901 :
902 : ! Fill the atom position of the first cell (reference cell)
903 39 : first_coordinate = ((irpt_ref-1)*natom_uc) + 1
904 39 : second_coordinate = first_coordinate + natom_uc-1
905 815 : xred_tmp(:,1:natom_uc) = xred(:,first_coordinate:second_coordinate)
906 : ! Fill fake zeff array for ewald9
907 2561 : zeff_tmp(:,:,1:natom_uc) = eff_pot%harmonics_terms%zeff
908 2561 : zeff_tmp(:,:,natom_uc+1:2*natom_uc) = eff_pot%harmonics_terms%zeff
909 :
910 538 : do irpt=1,my_nrpt
911 499 : i1=my_index_rpt(1,irpt); i2=my_index_rpt(2,irpt); i3=my_index_rpt(3,irpt)
912 : ! Compute new dipole-dipole interaction
913 4335507 : dyew = zero
914 538 : if (i1==0.and.i2==0.and.i3==0) then
915 : call ewald9(acell,eff_pot%harmonics_terms%epsilon_inf,dyewq0,&
916 : & gmet,gprimd,natom_uc,real((/0,0,0/),dp),rmet,&
917 : & supercell%rprimd,sumg0,ucvol,xred_tmp(:,1:natom_uc),&
918 27 : & eff_pot%harmonics_terms%zeff,qdrp_cart,eta)
919 29935 : buff_ewald(:,:,:,:,:,irpt) = dyewq0
920 : else
921 472 : first_coordinate = ((my_irpt(irpt)-1)*natom_uc) + 1
922 472 : second_coordinate = first_coordinate + natom_uc - 1
923 : xred_tmp(:,natom_uc+1:2*natom_uc)=&
924 13468 : & xred(:,first_coordinate:second_coordinate)
925 : call ewald9(acell,eff_pot%harmonics_terms%epsilon_inf,dyew,gmet,gprimd,&
926 : & int(2*natom_uc),real((/0,0,0/),dp),&
927 : & rmet,supercell%rprimd,&
928 472 : & sumg0,ucvol,xred_tmp,zeff_tmp,qdrp_cart,eta)
929 : buff_ewald(:,:,:,:,:,irpt) = &
930 1061098 : & dyew(:,:,1:natom_uc,:,natom_uc+1:2*natom_uc)
931 : end if
932 : end do
933 :
934 : ! DEALLOCATION OF ARRAYS
935 39 : ABI_FREE(my_index_rpt)
936 39 : ABI_FREE(my_irpt)
937 39 : ABI_FREE(xred_tmp)
938 39 : ABI_FREE(xred)
939 39 : ABI_FREE(zeff_tmp)
940 39 : ABI_FREE(qdrp_cart)
941 39 : ABI_FREE(dyew)
942 39 : ABI_FREE(dyewq0)
943 :
944 : ! Set the bufsize for mpi allgather
945 120 : do ii = 1,nproc
946 81 : bufsize(ii) = int(real(ifc_tmp%nrpt,sp)/nproc)*3*natom_uc*3*natom_uc
947 120 : if(ii > (nproc-nrpt_alone)) then
948 6 : bufsize(ii) = bufsize(ii) + 3*natom_uc*3*natom_uc
949 : end if
950 : end do
951 :
952 39 : bufdisp(1) = 0
953 81 : do ii = 2,nproc
954 81 : bufdisp(ii) = bufdisp(ii-1) + bufsize(ii-1)
955 : end do
956 :
957 39 : size_tmp = 3*natom_uc*3*natom_uc*my_nrpt
958 39 : call xmpi_allgatherv(buff_ewald(1,:,:,:,:,:),size_tmp,ifc_tmp%ewald_atmfrc,bufsize,bufdisp, comm, ierr)
959 :
960 39 : ABI_FREE(bufsize)
961 39 : ABI_FREE(bufdisp)
962 39 : ABI_FREE(buff_ewald)
963 :
964 :
965 :
966 : ! Fill the short range part (calculated previously) only master
967 : ! if(iam_master)then
968 : ! do irpt=1,ifc_tmp%nrpt
969 : ! do irpt2=1,eff_pot%harmonics_terms%ifcs%nrpt
970 : ! if(eff_pot%harmonics_terms%ifcs%cell(1,irpt2)==ifc_tmp%cell(1,irpt).and.&
971 : !& eff_pot%harmonics_terms%ifcs%cell(2,irpt2)==ifc_tmp%cell(2,irpt).and.&
972 : !& eff_pot%harmonics_terms%ifcs%cell(3,irpt2)==ifc_tmp%cell(3,irpt).and.&
973 : !& any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)) > tol20)) then
974 : ! ifc_tmp%short_atmfrc(:,:,:,:,irpt) = &
975 : !& eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)
976 : ! end if
977 : ! end do
978 : ! end do
979 : ! end if
980 :
981 :
982 : !call xmpi_bcast(ifc_tmp%short_atmfrc, master, comm, ierr)
983 : ! Maybe useless
984 39 : call xmpi_bcast(eff_pot%harmonics_terms%ifcs%short_atmfrc, master, comm, ierr)
985 :
986 : ! Compute total ifc
987 : !ifc_tmp%atmfrc = ifc_tmp%short_atmfrc + ifc_tmp%ewald_atmfrc
988 : ! Set the full cell according to the largest box
989 :
990 39 : min1f=min(min1, min1_cell)
991 39 : min2f=min(min2, min2_cell)
992 39 : min3f=min(min3, min3_cell)
993 :
994 39 : max1f=max(max1, max1_cell)
995 39 : max2f=max(max2, max2_cell)
996 39 : max3f=max(max3, max3_cell)
997 :
998 39 : full_nrpt= (max1f-min1f+1) * (max2f-min2f+1) * (max3f-min3f+1)
999 :
1000 :
1001 : !full_nrpt = max(ifc_tmp%nrpt,eff_pot%harmonics_terms%ifcs%nrpt)
1002 :
1003 117 : ABI_MALLOC(full_cell,(3,full_nrpt))
1004 2012679 : ABI_CALLOC(full_cell_atmfrc,(3,natom_uc,3,natom_uc,full_nrpt)) ! Allocate and set to 0
1005 2012640 : ABI_CALLOC(full_cell_short_atmfrc,(3,natom_uc,3,natom_uc,full_nrpt)) ! Allocate and set to 0
1006 2012640 : ABI_CALLOC(full_cell_ewald_atmfrc,(3,natom_uc,3,natom_uc,full_nrpt)) ! Allocate and set to 0
1007 : ! if ( full_nrpt == ifc_tmp%nrpt ) then
1008 : ! full_cell = ifc_tmp%cell
1009 : ! else
1010 : ! full_cell = eff_pot%harmonics_terms%ifcs%cell
1011 : ! end if
1012 : irpt=1
1013 190 : do i1=min1f, max1f
1014 983 : do i2= min2f, max2f
1015 6413 : do i3= min3f, max3f
1016 21876 : full_cell(:, irpt)= [i1, i2, i3]
1017 6262 : irpt=irpt+1
1018 : end do
1019 : end do
1020 : end do
1021 :
1022 39 : has_totFC = .False.
1023 :
1024 5210 : do irpt = 1,eff_pot%harmonics_terms%ifcs%nrpt
1025 194310 : if(any(abs(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt)) > tol8))then
1026 5009 : has_totFC = .True.
1027 5009 : cycle
1028 : end if
1029 : end do
1030 :
1031 39 : if (has_totFC .and. in_file_option==1) then !
1032 : ! print *, ' >>>>> HAS TOTAL FC'
1033 : ! ! Copy LR into total_atmfrc
1034 86 : do irpt=1,ifc_tmp%nrpt ! LR IRPT
1035 1043 : do irpt2=1, full_nrpt
1036 : if( ifc_tmp%cell(1,irpt)==full_cell(1,irpt2).and.&
1037 957 : & ifc_tmp%cell(2,irpt)==full_cell(2,irpt2).and.&
1038 71 : & ifc_tmp%cell(3,irpt)==full_cell(3,irpt2) ) then
1039 : ! full_cell_atmfrc(:,:,:,:,irpt2) = ifc_tmp%ewald_atmfrc(:,:,:,:,irpt)
1040 15623 : full_cell_ewald_atmfrc(:,:,:,:,irpt2) = ifc_tmp%ewald_atmfrc(:,:,:,:,irpt)
1041 : end if
1042 : end do
1043 : end do
1044 :
1045 : ! Copy total FC into total_atmfrc
1046 258 : do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt ! SR IRPT
1047 5379 : do irpt2=1, full_nrpt
1048 : if( eff_pot%harmonics_terms%ifcs%cell(1,irpt)==full_cell(1,irpt2).and.&
1049 5121 : & eff_pot%harmonics_terms%ifcs%cell(2,irpt)==full_cell(2,irpt2).and.&
1050 243 : & eff_pot%harmonics_terms%ifcs%cell(3,irpt)==full_cell(3,irpt2) ) then
1051 67251 : full_cell_atmfrc(:,:,:,:,irpt2) = eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt)
1052 12431 : if(any(abs(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt) - full_cell_ewald_atmfrc(:,:,:,:,irpt2)) > tol8))then
1053 67251 : full_cell_short_atmfrc(:,:,:,:,irpt2) = eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt) - full_cell_ewald_atmfrc(:,:,:,:,irpt2)
1054 : end if
1055 : end if
1056 : end do
1057 : end do
1058 : else
1059 : ! Copy LR into total_atmfrc
1060 520 : do irpt=1,ifc_tmp%nrpt ! LR IRPT
1061 125280 : do irpt2=1, full_nrpt
1062 : if( ifc_tmp%cell(1,irpt)==full_cell(1,irpt2).and.&
1063 124760 : & ifc_tmp%cell(2,irpt)==full_cell(2,irpt2).and.&
1064 496 : & ifc_tmp%cell(3,irpt)==full_cell(3,irpt2) ) then
1065 443888 : full_cell_atmfrc(:,:,:,:,irpt2) = ifc_tmp%ewald_atmfrc(:,:,:,:,irpt)
1066 443888 : full_cell_ewald_atmfrc(:,:,:,:,irpt2) = ifc_tmp%ewald_atmfrc(:,:,:,:,irpt)
1067 : end if
1068 : end do
1069 : end do
1070 :
1071 : ! Copy SR into total_atmfrc
1072 4952 : do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt ! SR IRPT
1073 3245430 : do irpt2=1, full_nrpt
1074 : if( eff_pot%harmonics_terms%ifcs%cell(1,irpt)==full_cell(1,irpt2).and.&
1075 3240478 : & eff_pot%harmonics_terms%ifcs%cell(2,irpt)==full_cell(2,irpt2).and.&
1076 4928 : & eff_pot%harmonics_terms%ifcs%cell(3,irpt)==full_cell(3,irpt2) ) then
1077 1695152 : full_cell_atmfrc(:,:,:,:,irpt2) = full_cell_atmfrc(:,:,:,:,irpt2) + eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt)
1078 1695152 : full_cell_short_atmfrc(:,:,:,:,irpt2) = eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt)
1079 : end if
1080 : end do
1081 : end do
1082 : end if
1083 :
1084 : ! Count the rpt inferior to the tolerance
1085 39 : irpt2 = 0
1086 5508 : do irpt=1,full_nrpt
1087 54208 : if(any(abs(full_cell_atmfrc(:,:,:,:,irpt)) > tol8))then
1088 5323 : irpt2 = irpt2 + 1
1089 : end if
1090 : end do
1091 : !write(std_out,*) "irpt2: how many dipdip cells: ", irpt2
1092 :
1093 : ! Copy ifc into effective potential
1094 : ! !!!Warning eff_pot%harmonics_terms%ifcs only contains atmfrc,short_atmfrc,ewald_atmfrc,nrpt
1095 : ! and cell!! rcan,ifc%rpt,wghatm and other quantities
1096 : ! are not needed for effective potential!!!
1097 : ! Free ifc before copy
1098 39 : call eff_pot%harmonics_terms%ifcs%free()
1099 :
1100 : ! Fill the effective potential with new atmfr
1101 39 : eff_pot%harmonics_terms%ifcs%nrpt = irpt2
1102 195 : ABI_MALLOC(eff_pot%harmonics_terms%ifcs%atmfrc,(3,natom_uc,3,natom_uc,irpt2))
1103 156 : ABI_MALLOC(eff_pot%harmonics_terms%ifcs%short_atmfrc,(3,natom_uc,3,natom_uc,irpt2))
1104 156 : ABI_MALLOC(eff_pot%harmonics_terms%ifcs%ewald_atmfrc,(3,natom_uc,3,natom_uc,irpt2))
1105 117 : ABI_MALLOC(eff_pot%harmonics_terms%ifcs%cell,(3,irpt2))
1106 :
1107 39 : irpt2 = 0
1108 5508 : do irpt = 1,full_nrpt
1109 54208 : if(any(abs(full_cell_atmfrc(:,:,:,:,irpt)) > tol8))then
1110 5323 : irpt2 = irpt2 + 1
1111 1977227 : eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt2) = full_cell_atmfrc(:,:,:,:,irpt)
1112 1977227 : eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2) = full_cell_short_atmfrc(:,:,:,:,irpt)
1113 1977227 : eff_pot%harmonics_terms%ifcs%ewald_atmfrc(:,:,:,:,irpt2) = full_cell_ewald_atmfrc(:,:,:,:,irpt)
1114 21292 : eff_pot%harmonics_terms%ifcs%cell(:,irpt2) = full_cell(:,irpt)
1115 : end if
1116 : end do
1117 :
1118 : ! Free temporary ifc
1119 39 : call ifc_tmp%free()
1120 : !Deallocate temporary arrays
1121 39 : ABI_FREE(full_cell)
1122 39 : ABI_FREE(full_cell_atmfrc) ! Allocate and set to 0
1123 39 : ABI_FREE(full_cell_short_atmfrc) ! Allocate and set to 0
1124 39 : ABI_FREE(full_cell_ewald_atmfrc) ! Allocate and set to 0
1125 :
1126 : end if
1127 :
1128 67 : if(asr >= 0) then
1129 : ! Impose sum rule
1130 67 : call harmonics_terms_applySumRule(asr,eff_pot%harmonics_terms%ifcs,natom_uc)
1131 : end if
1132 :
1133 5427 : write(msg, '(a,(80a),a)' ) ch10,('=',ii=1,80)
1134 67 : call wrtout(ab_out,msg,'COLL')
1135 67 : call wrtout(std_out,msg,'COLL')
1136 :
1137 : ! Free suppercell
1138 67 : call supercell%free()
1139 :
1140 67 : end subroutine effective_potential_generateDipDip
1141 : !!***
1142 :
1143 : !****f* m_effective_potential/effective_potential_setCoeffs
1144 : !!
1145 : !! NAME
1146 : !! effective_potential_setCoeffs
1147 : !!
1148 : !! FUNCTION
1149 : !! Set the polynomial coefficients of the effective_potential
1150 : !!
1151 : !! INPUTS
1152 : !! coeffs(ncoeff)<type(polynomial_coeff)> = array of polynomial_coeff datatype
1153 : !! ncoeff = number of coefficient
1154 : !!
1155 : !! OUTPUT
1156 : !! eff_pot<type(effective_potential_type)> = datatype for effective potential
1157 : !!
1158 : !! SOURCE
1159 :
1160 96 : subroutine effective_potential_setCoeffs(coeffs,eff_pot,ncoeff)
1161 :
1162 : !Arguments ------------------------------------
1163 : !scalars
1164 : integer,intent(in) :: ncoeff
1165 : !array
1166 : type(effective_potential_type),intent(inout) :: eff_pot
1167 : type(polynomial_coeff_type),intent(in) :: coeffs(ncoeff)
1168 : !Local variables-------------------------------
1169 : !scalar
1170 : integer :: ii,jj
1171 : logical :: has_straincoupling
1172 : character(len=500) :: msg
1173 : !array
1174 : ! *************************************************************************
1175 :
1176 96 : if(ncoeff /= size(coeffs))then
1177 0 : ABI_BUG('ncoeff has not the same size than coeffs array')
1178 : end if
1179 :
1180 : ! Check if the strain coupling is present
1181 96 : has_straincoupling=.FALSE.
1182 2155 : do ii=1,ncoeff
1183 57970 : do jj=1,coeffs(ii)%nterm
1184 57874 : if (coeffs(ii)%terms(jj)%nstrain > 0) then
1185 14831 : has_straincoupling = .TRUE.
1186 : end if
1187 : end do
1188 : end do
1189 :
1190 : ! Set to false the strain coupling from the finite differences
1191 96 : if(has_straincoupling)then
1192 33 : if(eff_pot%has_anharmonicsTerms) then
1193 12 : write(msg, '(8a)' )ch10,&
1194 12 : & ' --- !WARNING',ch10,&
1195 12 : & ' There is strain coupling with the fitted coefficients,',ch10,&
1196 12 : & ' The previous contribution will be set to zero',ch10,&
1197 24 : & ' ---'
1198 12 : call wrtout(std_out,msg,"COLL")
1199 : end if
1200 33 : eff_pot%has_anharmonicsTerms = .FALSE.
1201 : end if
1202 :
1203 96 : call anharmonics_terms_setCoeffs(coeffs,eff_pot%anharmonics_terms,ncoeff)
1204 :
1205 96 : end subroutine effective_potential_setCoeffs
1206 : !!***
1207 :
1208 : !****f* m_effective_potential/effective_potential_setElastic3rd
1209 : !!
1210 : !! NAME
1211 : !! effective_potential_setElastic3rd
1212 : !!
1213 : !! FUNCTION
1214 : !! Set the 3rd order derivative of with respect to 3 strain
1215 : !!
1216 : !! INPUTS
1217 : !! elastics(6,6,6) = 3d order of elastics constant
1218 : !!
1219 : !! OUTPUT
1220 : !! eff_pot<type(effective_potential_type)> = datatype for effective potential
1221 : !!
1222 : !!
1223 : !! SOURCE
1224 :
1225 0 : subroutine effective_potential_setElastic3rd(eff_pot,elastics)
1226 :
1227 : !Arguments ------------------------------------
1228 : !scalars
1229 : !array
1230 : real(dp),intent(in) :: elastics(6,6,6)
1231 : type(effective_potential_type),intent(inout) :: eff_pot
1232 : !Local variables-------------------------------
1233 : !scalar
1234 : !array
1235 : ! *************************************************************************
1236 0 : call anharmonics_terms_setElastic3rd(eff_pot%anharmonics_terms,elastics)
1237 :
1238 0 : if(any(abs(eff_pot%anharmonics_terms%elastic3rd)> tol15)) then
1239 0 : eff_pot%has_anharmonicsTerms = .TRUE.
1240 : end if
1241 :
1242 0 : end subroutine effective_potential_setElastic3rd
1243 : !!***
1244 :
1245 : !****f* m_effective_potential/effective_potential_setElastic4th
1246 : !!
1247 : !! NAME
1248 : !! effective_potential_setElastic4th
1249 : !!
1250 : !! FUNCTION
1251 : !! Set the 4th order derivative of with respect to 4 strain
1252 : !!
1253 : !! INPUTS
1254 : !! elastics(6,6,6,6) = 4th order of elastics constant
1255 : !!
1256 : !! OUTPUT
1257 : !! eff_pot<type(effective_potential_type)> = datatype for effective potential
1258 : !!
1259 : !!
1260 : !! SOURCE
1261 :
1262 0 : subroutine effective_potential_setElastic4th(eff_pot,elastics)
1263 :
1264 : !Arguments ------------------------------------
1265 : !scalars
1266 : !array
1267 : real(dp),intent(in) :: elastics(6,6,6,6)
1268 : type(effective_potential_type),intent(inout) :: eff_pot
1269 : !Local variables-------------------------------
1270 : !scalar
1271 : !array
1272 : ! *************************************************************************
1273 0 : call anharmonics_terms_setElastic4th(eff_pot%anharmonics_terms,elastics)
1274 :
1275 0 : if(any(abs(eff_pot%anharmonics_terms%elastic4th)> tol15)) then
1276 0 : eff_pot%has_anharmonicsTerms = .TRUE.
1277 : end if
1278 :
1279 :
1280 0 : end subroutine effective_potential_setElastic4th
1281 : !!***
1282 :
1283 : !****f* m_effective_potential/effective_potential_setStrainPhononCoupling
1284 : !!
1285 : !! NAME
1286 : !! effective_potential_setStrainPhononCoupling
1287 : !!
1288 : !! FUNCTION
1289 : !! Set the strain phonon coupling of the effective_potential
1290 : !!
1291 : !! INPUTS
1292 : !! natom = number of atoms
1293 : !! strain_phonon(6)<type(ifc_type) = strain-phonon coupling
1294 : !!
1295 : !! OUTPUT
1296 : !! eff_pot<type(effective_potential_type)> = datatype for effective potential
1297 : !!
1298 : !! SOURCE
1299 :
1300 0 : subroutine effective_potential_setStrainPhononCoupling(eff_pot,natom,phonon_strain)
1301 :
1302 : !Arguments ------------------------------------
1303 : !scalars
1304 : integer, intent(in) :: natom
1305 : !array
1306 : type(ifc_type),intent(in) :: phonon_strain(6)
1307 : type(effective_potential_type),intent(inout) :: eff_pot
1308 : !Local variables-------------------------------
1309 : !scalar
1310 : integer :: ii
1311 : !array
1312 : ! *************************************************************************
1313 :
1314 0 : call anharmonics_terms_setStrainPhononCoupling(eff_pot%anharmonics_terms,natom,phonon_strain)
1315 0 : do ii=1,6
1316 0 : if(any(abs(eff_pot%anharmonics_terms%phonon_strain(ii)%atmfrc)> tol15)) then
1317 0 : eff_pot%has_anharmonicsTerms = .True.
1318 : end if
1319 : end do
1320 :
1321 0 : end subroutine effective_potential_setStrainPhononCoupling
1322 : !!***
1323 :
1324 : !****f* m_effective_potential/effective_potential_setElasticDispCoupling
1325 : !!
1326 : !! NAME
1327 : !! effective_potential_setElasticDispCoupling
1328 : !!
1329 : !! FUNCTION
1330 : !! Set the elastic constant displacement coupling of the effective_potential
1331 : !!
1332 : !! INPUTS
1333 : !! natom = number of atoms
1334 : !! elastic_displacement = (6,6,3,natom) array with elastic constant displacement coupling
1335 : !!
1336 : !! OUTPUT
1337 : !! eff_pot<type(effective_potential_type)> = datatype for effective potential
1338 : !!
1339 : !! SOURCE
1340 :
1341 0 : subroutine effective_potential_setElasticDispCoupling(eff_pot,natom,elastic_displacement)
1342 :
1343 : !Arguments ------------------------------------
1344 : !scalars
1345 : integer, intent(in) :: natom
1346 : !array
1347 : real(dp) :: elastic_displacement(6,6,3,natom)
1348 : type(effective_potential_type),intent(inout) :: eff_pot
1349 : !Local variables-------------------------------
1350 : !scalar
1351 : !array
1352 : ! *************************************************************************
1353 :
1354 0 : call anharmonics_terms_setElasticDispCoupling(eff_pot%anharmonics_terms,natom,elastic_displacement)
1355 0 : if(any(abs(eff_pot%anharmonics_terms%elastic_displacement)> tol15)) then
1356 0 : eff_pot%has_anharmonicsTerms = .True.
1357 : end if
1358 :
1359 0 : end subroutine effective_potential_setElasticDispCoupling
1360 : !!***
1361 :
1362 : !!****f* m_effective_potential/effective_potential_setConfinement
1363 : !!
1364 : !! NAME
1365 : !! effective_potential_setConfinement
1366 : !!
1367 : !! FUNCTION
1368 : !! Set the confinement in the effective_potential datatype
1369 : !!
1370 : !! INPUTS
1371 : !! cutoff_disp(6) = Cutoff array for the strain
1372 : !! cutoff_strain(ndisp) = Cutoff array for the atomic displacement
1373 : !! factor_disp = Factor to appy to the polynomial term of the confinement (displacement)
1374 : !! factor_strain = Factor to appy to the polynomial term of the confinement (strain)
1375 : !! ndisp = Number of displacement (atoms) for the cut off
1376 : !! power_disp = Power of the polynome related to the displacement
1377 : !! power_strain = Power of the polynome related to the strain
1378 : !! need_confinement = optional,Logical related to the necessity of the confinement
1379 : !!
1380 : !! OUTPUT
1381 : !! eff_pot<type(effective_potential_type)> = datatype for effective potential
1382 : !!
1383 : !! SOURCE
1384 :
1385 0 : subroutine effective_potential_setConfinement(cutoff_disp,cutoff_strain,eff_pot,factor_disp,&
1386 : & factor_strain,ndisp,power_disp,power_strain,&
1387 : & need_confinement)
1388 :
1389 : !Arguments ------------------------------------
1390 : !scalars
1391 : integer, intent(in) :: power_disp,power_strain,ndisp
1392 : real(dp),intent(in) :: factor_disp,factor_strain
1393 : logical,optional,intent(in) :: need_confinement
1394 : !arrays
1395 : real(dp),intent(in) :: cutoff_disp(ndisp),cutoff_strain(6)
1396 : type(effective_potential_type),intent(inout) :: eff_pot
1397 : !Local variables-------------------------------
1398 : !scalar
1399 : logical :: need_confinement_tmp
1400 : !arrays
1401 : !character(len=500) :: msg
1402 :
1403 : ! *************************************************************************
1404 :
1405 : !Checks
1406 0 : if (ndisp <= 0) then
1407 0 : ABI_ERROR('ndisp can not be inferior or equal to zero')
1408 : end if
1409 :
1410 : !First free the type
1411 0 : call polynomial_conf_free(eff_pot%confinement)
1412 :
1413 0 : need_confinement_tmp = .FALSE.
1414 0 : if (present(need_confinement)) need_confinement_tmp = need_confinement
1415 :
1416 : call polynomial_conf_init(cutoff_disp,cutoff_strain,factor_disp,factor_strain,ndisp,&
1417 : & eff_pot%confinement,power_disp,power_strain,&
1418 0 : & need_confinement=need_confinement_tmp)
1419 :
1420 :
1421 0 : end subroutine effective_potential_setConfinement
1422 : !!***
1423 :
1424 : !!****f* m_effective_potential/effective_potential_setSupercell
1425 : !!
1426 : !! NAME
1427 : !! effective_potential_setSupercell
1428 : !!
1429 : !! FUNCTION
1430 : !! Set the supercell type in the effective_potential type
1431 : !!
1432 : !! INPUTS
1433 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
1434 : !! comm = MPI communicator
1435 : !! ncell(3) = optional, size of the supercell
1436 : !! supercell = optional, supercell type to set to eff_pot
1437 : !!
1438 : !! OUTPUT
1439 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
1440 : !!
1441 : !! SOURCE
1442 :
1443 154 : subroutine effective_potential_setSupercell(eff_pot,comm,ncell,supercell)
1444 :
1445 : !Arguments ------------------------------------
1446 : !scalars
1447 : integer,intent(in) :: comm
1448 : !arrays
1449 : type(effective_potential_type),intent(inout) :: eff_pot
1450 : integer,optional,intent(in) :: ncell(3)
1451 : type(supercell_type),optional,intent(in) :: supercell
1452 : !Local variables-------------------------------
1453 : !scalar
1454 : !arrays
1455 : !character(len=500) :: msg
1456 :
1457 : ! *************************************************************************
1458 :
1459 : !Checks
1460 154 : if (.not.present(supercell).and..not.present(ncell)) then
1461 0 : ABI_ERROR(' You should at least set ncell of supercell type')
1462 : end if
1463 :
1464 154 : call eff_pot%supercell%free()
1465 :
1466 154 : if(present(supercell))then
1467 38 : call supercell%copy(eff_pot%supercell)
1468 : else
1469 : call eff_pot%supercell%init(eff_pot%crystal%natom, (/ncell(1),0,0, 0,ncell(2),0, 0,0,ncell(3)/), &
1470 : eff_pot%crystal%rprimd,eff_pot%crystal%typat,eff_pot%crystal%xcart,&
1471 1160 : eff_pot%crystal%znucl)
1472 : end if
1473 :
1474 : !Initialisation of new mpi over supercell
1475 154 : call effective_potential_initmpi(eff_pot,comm)
1476 :
1477 154 : end subroutine effective_potential_setSupercell
1478 : !!***
1479 :
1480 : !****f* m_effective_potential/effective_potential_print
1481 : !!
1482 : !! NAME
1483 : !! effective_potential_print
1484 : !!
1485 : !! FUNCTION
1486 : !! Print the effective_potential in the ouput
1487 : !!
1488 : !! INPUTS
1489 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
1490 : !! option = 0 no output
1491 : !! option = 1 only useful information
1492 : !!
1493 : !! OUTPUT
1494 : !!
1495 : !!
1496 : !! SOURCE
1497 :
1498 46 : subroutine effective_potential_print(eff_pot,option,filename)
1499 :
1500 : !Arguments ------------------------------------
1501 : !scalars
1502 : integer,intent(in) :: option
1503 : character(len=*),optional,intent(in) :: filename
1504 : !array
1505 : type(effective_potential_type),intent(inout) :: eff_pot
1506 : !Local variables-------------------------------
1507 : !scalar
1508 : integer :: ia,ii
1509 : real(dp):: fact
1510 : character(len=500) :: msg
1511 : !array
1512 : ! *************************************************************************
1513 :
1514 46 : if(option >= 1) then
1515 33 : if(present(filename)) then
1516 0 : write(msg, '(a,a,a,a,a,a)' )ch10,' The file ',trim(filename),&
1517 0 : & ' contains this effective potential for ',trim(eff_pot%name),':'
1518 : else
1519 33 : write(msg, '(a,a,a,a)' )ch10,' This effective potential contains ',&
1520 66 : & trim(eff_pot%name),':'
1521 : end if
1522 :
1523 33 : call wrtout(std_out,msg,'COLL')
1524 33 : call wrtout(ab_out,msg,'COLL')
1525 :
1526 : !**********************************************************************
1527 : ! Write basics values
1528 : !**********************************************************************
1529 :
1530 : write(msg,'(a,F20.10,2a,I3,2a,I4,2a,I4,2a,I3,2a)') &
1531 33 : & ' - Reference energy: ',eff_pot%energy ,ch10,&
1532 33 : & ' - Number of types of atoms: ',eff_pot%crystal%ntypat ,ch10,&
1533 33 : & ' - Number of atoms: ',eff_pot%crystal%natom ,ch10,&
1534 33 : & ' - Number of cells: ',eff_pot%harmonics_terms%ifcs%nrpt ,ch10,&
1535 33 : & ' - Number of qpoints: ',eff_pot%harmonics_terms%nqpt ,ch10,&
1536 66 : & ' - Primitive vectors (unit:Bohr): '
1537 33 : call wrtout(ab_out,msg,'COLL')
1538 33 : call wrtout(std_out,msg,'COLL')
1539 132 : do ii = 1,3
1540 99 : write(msg,'(3F12.6)') eff_pot%crystal%rprimd(1,ii),&
1541 99 : & eff_pot%crystal%rprimd(2,ii),&
1542 198 : & eff_pot%crystal%rprimd(3,ii)
1543 99 : call wrtout(ab_out,msg,'COLL')
1544 132 : call wrtout(std_out,msg,'COLL')
1545 : end do
1546 33 : write(msg,'(2a,3F12.6)') ' - acell (unit:Bohr):',ch10,one,one,one
1547 :
1548 33 : call wrtout(ab_out,msg,'COLL')
1549 33 : call wrtout(std_out,msg,'COLL')
1550 33 : write(msg,'(a)') ' - Dielectric tensor: '
1551 33 : call wrtout(ab_out,msg,'COLL')
1552 33 : call wrtout(std_out,msg,'COLL')
1553 132 : do ii=1,3
1554 99 : write(msg,'(3F12.6)')eff_pot%harmonics_terms%epsilon_inf(1,ii),&
1555 99 : & eff_pot%harmonics_terms%epsilon_inf(2,ii),&
1556 198 : & eff_pot%harmonics_terms%epsilon_inf(3,ii)
1557 99 : call wrtout(ab_out,msg,'COLL')
1558 132 : call wrtout(std_out,msg,'COLL')
1559 : end do
1560 33 : write(msg,'(a)') ' - Elastic tensor (unit:10^2GPa): '
1561 33 : call wrtout(ab_out,msg,'COLL')
1562 33 : call wrtout(std_out,msg,'COLL')
1563 33 : fact = HaBohr3_GPa / eff_pot%crystal%ucvol
1564 231 : do ii=1,6
1565 : write(msg,'(6F12.6)')&
1566 198 : & eff_pot%harmonics_terms%elastic_constants(1,ii)*fact/100,&
1567 198 : & eff_pot%harmonics_terms%elastic_constants(2,ii)*fact/100,&
1568 198 : & eff_pot%harmonics_terms%elastic_constants(3,ii)*fact/100,&
1569 198 : & eff_pot%harmonics_terms%elastic_constants(4,ii)*fact/100,&
1570 198 : & eff_pot%harmonics_terms%elastic_constants(5,ii)*fact/100,&
1571 396 : & eff_pot%harmonics_terms%elastic_constants(6,ii)*fact/100
1572 198 : call wrtout(ab_out,msg,'COLL')
1573 231 : call wrtout(std_out,msg,'COLL')
1574 : end do
1575 197 : do ia=1,eff_pot%crystal%natom
1576 164 : write(msg,'(a,I4,2a,F10.4,2a,F10.4,2a,3F12.6,2a)')' - Atoms',ia,ch10,&
1577 164 : & " - atomic number:",eff_pot%crystal%znucl(eff_pot%crystal%typat(ia)),ch10,&
1578 164 : & " - atomic mass:",eff_pot%crystal%amu(eff_pot%crystal%typat(ia)),ch10,&
1579 164 : & " - cartesian position:",eff_pot%crystal%xcart(:,ia),ch10,&
1580 328 : & " - Effective charges:"
1581 164 : call wrtout(ab_out,msg,'COLL')
1582 164 : call wrtout(std_out,msg,'COLL')
1583 689 : do ii = 1,3
1584 492 : write(msg,'(a,3(F12.6))') " ",eff_pot%harmonics_terms%zeff(:,ii,ia)
1585 492 : call wrtout(ab_out,msg,'COLL')
1586 656 : call wrtout(std_out,msg,'COLL')
1587 : end do
1588 : end do
1589 : end if
1590 :
1591 46 : end subroutine effective_potential_print
1592 : !!***
1593 :
1594 : !****f* m_effective_potential/effective_potential_printSupercell
1595 : !!
1596 : !! NAME
1597 : !! effective_potential_printSupercell
1598 : !!
1599 : !! FUNCTION
1600 : !! Print the supercell of the effective_potential
1601 : !! or if present the supercell as input
1602 : !! WARNING: need to be consistent with eff_pot
1603 : !!
1604 : !! INPUTS
1605 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
1606 : !! supercell<type(supercell_type)> = optional, supercell type to define
1607 : !!
1608 : !! OUTPUT
1609 : !!
1610 : !! SOURCE
1611 :
1612 0 : subroutine effective_potential_printSupercell(eff_pot,supercell)
1613 :
1614 : !Arguments ------------------------------------
1615 : !scalars
1616 : !array
1617 : type(effective_potential_type),target,intent(inout) :: eff_pot
1618 : type(supercell_type),optional,target,intent(in) :: supercell
1619 : !Local variables-------------------------------
1620 : !scalar
1621 : integer :: iatom,ii
1622 : character(len=500) :: msg
1623 : !array
1624 0 : real(dp), allocatable :: xred(:,:)
1625 : type(supercell_type),pointer :: supercell_tmp
1626 :
1627 : ! *************************************************************************
1628 :
1629 0 : if(present(supercell)) then
1630 : supercell_tmp => supercell
1631 : else
1632 0 : supercell_tmp => eff_pot%supercell
1633 : end if
1634 :
1635 0 : if(supercell_tmp%natom/= eff_pot%supercell%natom) then
1636 : write(msg, '(3a)' )&
1637 0 : & ' There is not the same numbers of atoms in the two supercell',ch10,&
1638 0 : & 'Action: modify the code'
1639 0 : ABI_BUG(msg)
1640 : end if
1641 :
1642 0 : ABI_MALLOC(xred,(3,supercell_tmp%natom))
1643 :
1644 : !**********************************************************************
1645 : ! Write basics values
1646 : !**********************************************************************
1647 :
1648 0 : write (msg, '(4a,I8,a)') ' Structure parameters of the supercell :',ch10,ch10,&
1649 0 : ' natom ', supercell_tmp%natom,ch10
1650 0 : call wrtout(ab_out,msg,'COLL')
1651 0 : call wrtout(std_out,msg,'COLL')
1652 :
1653 0 : write (msg, '(a)') ' znucl '
1654 0 : call wrtout(ab_out,msg,'COLL')
1655 0 : call wrtout(std_out,msg,'COLL')
1656 0 : write(msg,*) ''
1657 0 : do iatom = 1, size(eff_pot%crystal%znucl)
1658 0 : write (msg, '(a,I5)') trim(msg),int(eff_pot%crystal%znucl(iatom))
1659 0 : if (mod(iatom,6) == 0) then
1660 0 : call wrtout(ab_out,msg,'COLL')
1661 0 : call wrtout(std_out,msg,'COLL')
1662 0 : write(msg,*) ''
1663 : end if
1664 : end do
1665 0 : write (msg, '(2a)') trim(msg),ch10
1666 0 : call wrtout(ab_out,msg,'COLL')
1667 0 : call wrtout(std_out,msg,'COLL')
1668 0 : write (msg, '(a,I7,3a)') ' ntypat', size(eff_pot%crystal%znucl),ch10,ch10, ' typat '
1669 0 : call wrtout(ab_out,msg,'COLL')
1670 0 : call wrtout(std_out,msg,'COLL')
1671 :
1672 0 : write(msg,*) ''
1673 0 : do iatom = 1, supercell_tmp%natom
1674 0 : write (msg, '(a,I5)') trim(msg),&
1675 0 : & supercell_tmp%typat(supercell_tmp%atom_indexing(iatom))
1676 0 : if (mod(iatom,12) == 0)then
1677 0 : call wrtout(ab_out,msg,'COLL')
1678 0 : call wrtout(std_out,msg,'COLL')
1679 0 : write(msg,*) ''
1680 : end if
1681 : end do
1682 0 : write (msg, '(2a)') trim(msg),ch10
1683 0 : call wrtout(ab_out,msg,'COLL')
1684 0 : call wrtout(std_out,msg,'COLL')
1685 0 : write (msg, '(3a)') ' acell 1.0 1.0 1.0',ch10,ch10
1686 0 : write (msg, '(2a)') trim(msg),' rprim'
1687 0 : call wrtout(ab_out,msg,'COLL')
1688 0 : call wrtout(std_out,msg,'COLL')
1689 :
1690 0 : do ii = 1,3
1691 0 : write(msg,'(3E23.14,3E23.14,3E23.14)') supercell_tmp%rprimd(1,ii),&
1692 0 : & supercell_tmp%rprimd(2,ii),&
1693 0 : & supercell_tmp%rprimd(3,ii)
1694 0 : call wrtout(ab_out,msg,'COLL')
1695 0 : call wrtout(std_out,msg,'COLL')
1696 : end do
1697 :
1698 0 : write (msg, '(2a)') ch10,' xcart'
1699 0 : call wrtout(ab_out,msg,'COLL')
1700 0 : call wrtout(std_out,msg,'COLL')
1701 0 : do iatom = 1, supercell_tmp%natom
1702 0 : write (msg, '(3E23.14)') supercell_tmp%xcart(1,iatom),&
1703 0 : & supercell_tmp%xcart(2,iatom),&
1704 0 : & supercell_tmp%xcart(3,iatom)
1705 0 : call wrtout(ab_out,msg,'COLL')
1706 0 : call wrtout(std_out,msg,'COLL')
1707 : end do
1708 : call xcart2xred(supercell_tmp%natom,supercell_tmp%rprimd,&
1709 0 : & supercell_tmp%xcart,xred)
1710 0 : write (msg, '(2a)') ch10,' xred'
1711 0 : call wrtout(ab_out,msg,'COLL')
1712 0 : call wrtout(std_out,msg,'COLL')
1713 0 : do iatom = 1, supercell_tmp%natom
1714 0 : write (msg, '(3E23.14)') xred(1,iatom),&
1715 0 : & xred(2,iatom),&
1716 0 : & xred(3,iatom)
1717 0 : call wrtout(ab_out,msg,'COLL')
1718 0 : call wrtout(std_out,msg,'COLL')
1719 : end do
1720 :
1721 0 : write (msg, '(a)') ''
1722 0 : call wrtout(ab_out,msg,'COLL')
1723 0 : call wrtout(std_out,msg,'COLL')
1724 :
1725 : ! Deallocation array
1726 0 : ABI_FREE(xred)
1727 :
1728 0 : end subroutine effective_potential_printSupercell
1729 : !!***
1730 :
1731 : !!****f* m_effective_potential/effective_potential_writeXML
1732 : !! NAME
1733 : !! effective_potential_writeXML
1734 : !!
1735 : !! FUNCTION
1736 : !! This routine print the effective potential into xml format
1737 : !! Several options are available
1738 : !!
1739 : !! INPUTS
1740 : !! filename = the name of output file
1741 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
1742 : !! option = 0 Do nothing
1743 : !! = 1 Generate the XML file with:
1744 : !! - The system definition and the model (Harmonic + Anharmonic)
1745 : !! = 2 Generate two XML files with:
1746 : !! - The system definition and the model (Harmonic)
1747 : !! - The model (Anharmonic)
1748 : !! = 3 Generate one XML files with:
1749 : !! - The system definition and the model (Harmonic)
1750 : !! = 4 Generate one XML files with:
1751 : !! - The model (Anharmonic)
1752 : !!
1753 : !! OUTPUT
1754 : !!
1755 : !! SOURCE
1756 :
1757 29 : subroutine effective_potential_writeXML(eff_pot,option,filename,prt_dipdip)
1758 :
1759 : !Arguments ------------------------------------
1760 : !scalars
1761 : integer, intent(in) :: option
1762 : character(len=fnlen),optional,intent(in) :: filename
1763 : logical,optional,intent(in) :: prt_dipdip
1764 : !arrays
1765 : type(effective_potential_type), intent(in) :: eff_pot
1766 :
1767 : !Local variables-------------------------------
1768 : !scalar
1769 : integer :: ii,ia,ib,jj
1770 : integer :: iqpt,irpt,mu,nu
1771 : integer :: unit_xml
1772 : character(len=500) :: msg
1773 : character(len=fnlen) :: namefile
1774 : character(len=10) :: natom
1775 : logical :: new_file,need_prtdipdip
1776 : !arrays
1777 : real(dp) :: strain(9,6)
1778 :
1779 : ! *************************************************************************
1780 :
1781 290 : strain(:,1) = (/1,0,0,0,0,0,0,0,0/)
1782 290 : strain(:,2) = (/0,0,0,0,1,0,0,0,0/)
1783 290 : strain(:,3) = (/0,0,0,0,0,0,0,0,1/)
1784 290 : strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/)
1785 290 : strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/)
1786 290 : strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/)
1787 :
1788 29 : unit_xml = get_unit()
1789 29 : need_prtdipdip = .TRUE.
1790 29 : if(present(prt_dipdip)) need_prtdipdip = prt_dipdip
1791 : !Print only the reference system in xml format
1792 29 : if (option == 1 .or. option == 2 .or. option ==3) then
1793 :
1794 : ! convert natom in character
1795 27 : write (natom,'(I9)') eff_pot%crystal%natom
1796 :
1797 : ! Compute the name of the XML file
1798 27 : if(present(filename)) then
1799 27 : namefile=replace(trim(filename),".out","")
1800 10 : select case(option)
1801 : case(1)
1802 10 : namefile=trim(namefile)//"_model.xml"
1803 : case(2)
1804 10 : namefile=trim(namefile)//"_sys.xml"
1805 : case(3)
1806 27 : namefile=trim(namefile)//"_sys.xml"
1807 : end select
1808 : else
1809 0 : namefile='system.xml'
1810 : end if
1811 :
1812 27 : call isfile(namefile,'new')
1813 :
1814 27 : if (open_file(namefile,msg,unit=unit_xml,form="formatted",&
1815 : & status="new",action="write") /= 0) then
1816 0 : ABI_ERROR(msg)
1817 : end if
1818 :
1819 27 : write(msg,'(a,a,a)')ch10,&
1820 54 : & ' Generation of the xml file for the model in ',trim(namefile)
1821 :
1822 27 : call wrtout(ab_out,msg,'COLL')
1823 27 : call wrtout(std_out,msg,'COLL')
1824 :
1825 : ! Write header
1826 27 : WRITE(unit_xml,'("<?xml version=""1.0"" ?>")')
1827 27 : WRITE(unit_xml,'("<System_definition>")')
1828 :
1829 27 : WRITE(unit_xml,'(" <energy>")')
1830 27 : WRITE(unit_xml,'(E23.14)') (eff_pot%energy)
1831 27 : WRITE(unit_xml,'(" </energy>")')
1832 :
1833 27 : WRITE(unit_xml,'(" <unit_cell units=""bohrradius"">")')
1834 108 : do mu=1,3
1835 351 : WRITE(unit_xml,'(3(E23.14))') (eff_pot%crystal%rprimd(mu,nu),nu=1,3)
1836 : end do
1837 27 : WRITE(unit_xml,'(" </unit_cell>")')
1838 :
1839 27 : WRITE(unit_xml,'(" <epsilon_inf units=""epsilon0"">")')
1840 351 : WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%epsilon_inf)
1841 27 : WRITE(unit_xml,'(" </epsilon_inf>")')
1842 :
1843 27 : WRITE(unit_xml,'(" <elastic units=""hartree"">")')
1844 1161 : WRITE(unit_xml,'(6(E23.14))') (eff_pot%harmonics_terms%elastic_constants)
1845 27 : WRITE(unit_xml,'(" </elastic>")')
1846 :
1847 180 : do ia=1,eff_pot%crystal%natom
1848 : WRITE(unit_xml,'(" <atom mass=""",1F10.5,""" massunits=""atomicmassunit"">")') &
1849 153 : & eff_pot%crystal%amu(eff_pot%crystal%typat(ia))
1850 153 : WRITE(unit_xml,'(" <position units=""bohrradius"">")')
1851 612 : WRITE(unit_xml,'(3(E23.14))') (eff_pot%crystal%xcart(:,ia))
1852 153 : WRITE(unit_xml,'(" </position>")')
1853 153 : WRITE(unit_xml,'(" <borncharge units=""abs(e)"">")')
1854 1989 : WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%zeff(:,:,ia))
1855 153 : WRITE(unit_xml,'(" </borncharge>")')
1856 180 : WRITE(unit_xml,'(" </atom>")')
1857 : end do
1858 : ! Print the Ifc short range for each cell data is array 3*natom*3*natom
1859 : ! [ [x1 x2 ....]
1860 : ! [y1 y2 ....] for atom 1
1861 : ! [z1 z2 ....]
1862 : ! [x1 x2 ....]
1863 : ! [y1 y2 ....] for atom 2
1864 : ! [z1 z2 ....]
1865 : ! .... ]
1866 : ! Warning : The IFC are print in other order 1,mu,ia,nu,ib which is not fortran way
1867 : ! We do like that to because the previous script was in python
1868 : ! When you read ifc from XML file with fotran, you have to tranpose the matrix
1869 : !
1870 794 : do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt
1871 40723 : if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol20)) then
1872 719 : WRITE(unit_xml,'(" <local_force_constant units=""hartree/bohrradius**2"">")')
1873 719 : WRITE(unit_xml,'(" <data>")')
1874 5020 : do ia=1,eff_pot%crystal%natom
1875 17923 : do mu=1,3
1876 124770 : do ib=1,eff_pot%crystal%natom
1877 460371 : do nu=1,3
1878 : WRITE(unit_xml,'(e22.14)', advance="no")&
1879 447468 : & (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,ia,nu,ib,irpt))
1880 : end do
1881 : end do
1882 17204 : WRITE(unit_xml,'(a)')''
1883 : end do
1884 : end do
1885 719 : WRITE(unit_xml,'(" </data>")')
1886 719 : WRITE(unit_xml,'(" <cell>")')
1887 2876 : WRITE(unit_xml,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
1888 719 : WRITE(unit_xml,'(" </cell>")')
1889 719 : WRITE(unit_xml,'(" </local_force_constant>")')
1890 : end if
1891 : ! Print the IFC total for each cell, data is array 3*natom*3*natom
1892 : ! [ [x1 x2 ....]
1893 : ! [y1 y2 ....] for atom 1
1894 : ! [z1 z2 ....]
1895 : ! [x1 x2 ....]
1896 : ! [y1 y2 ....] for atom 2
1897 : ! [z1 z2 ....]
1898 : ! .... ]
1899 794 : if(need_prtdipdip)then
1900 30157 : if(all(abs(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt))<tol20)) then
1901 0 : if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol20)) then
1902 : write(msg, '(a,a,a,a)' )&
1903 0 : & ' There is no total range but short range in your effective potential',ch10,&
1904 0 : & 'Action: contact abinit group'
1905 0 : ABI_BUG(msg)
1906 : end if
1907 : else
1908 665 : WRITE(unit_xml,'(" <total_force_constant units=""hartree/bohrradius**2"">")')
1909 665 : WRITE(unit_xml,'(" <data>")')
1910 3886 : do ia=1,eff_pot%crystal%natom
1911 13549 : do mu=1,3
1912 56730 : do ib=1,eff_pot%crystal%natom
1913 197931 : do nu=1,3
1914 : WRITE(unit_xml,'(e22.14)', advance="no")&
1915 188268 : & (eff_pot%harmonics_terms%ifcs%atmfrc(mu,ia,nu,ib,irpt))
1916 : end do
1917 : end do
1918 12884 : WRITE(unit_xml,'(a)')''
1919 : end do
1920 : end do
1921 665 : WRITE(unit_xml,'(" </data>")')
1922 665 : WRITE(unit_xml,'(" <cell>")')
1923 2660 : WRITE(unit_xml,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
1924 665 : WRITE(unit_xml,'(" </cell>")')
1925 665 : WRITE(unit_xml,'(" </total_force_constant>")')
1926 : end if
1927 : end if
1928 : end do
1929 :
1930 134 : do iqpt=1,eff_pot%harmonics_terms%nqpt
1931 107 : WRITE(unit_xml,'(" <phonon>")')
1932 107 : WRITE(unit_xml,'(" <qpoint units=""2pi*G0"">")')
1933 428 : WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%qpoints(:,iqpt))
1934 107 : WRITE(unit_xml,'(" </qpoint>")')
1935 107 : WRITE(unit_xml,'(" <frequencies units=""reciprocal cm"">")')
1936 1514 : WRITE(unit_xml,'(3(e22.14))') (eff_pot%harmonics_terms%phfrq(:,iqpt))
1937 107 : WRITE(unit_xml,'(" </frequencies>")')
1938 107 : WRITE(unit_xml,'(" <dynamical_matrix units=""hartree/bohrradius**2"">")')
1939 576 : do ia=1,eff_pot%crystal%natom
1940 1983 : do mu=1,3
1941 9090 : do ib = 1,eff_pot%crystal%natom
1942 32139 : do nu=1,3
1943 30732 : WRITE(unit_xml,'(e22.14)',advance='no')(eff_pot%harmonics_terms%dynmat(1,nu,ib,mu,ia,iqpt))
1944 : end do
1945 : end do
1946 1876 : WRITE(unit_xml,'(a)')''
1947 : end do
1948 : end do
1949 107 : WRITE(unit_xml,'(" </dynamical_matrix>")')
1950 134 : WRITE(unit_xml,'(" </phonon>")')
1951 : end do
1952 : ! if phonon/forces strain is computed
1953 : jj = 1
1954 189 : do ii = 1,6
1955 162 : WRITE(unit_xml,'(" <strain_coupling voigt=""",I2,""">")') ii-1
1956 162 : WRITE(unit_xml,'(" <strain>")')
1957 1620 : WRITE(unit_xml,'(6(e12.4))') (strain(:,jj))
1958 162 : WRITE(unit_xml,'(" </strain>")')
1959 162 : WRITE(unit_xml,'(" <correction_force units=""hartree/bohrradius"">")')
1960 1080 : do ia=1,eff_pot%crystal%natom
1961 3672 : do mu=1,3
1962 : WRITE(unit_xml,'(e22.14)', advance="no")&
1963 3672 : & (eff_pot%harmonics_terms%strain_coupling(ii,mu,ia))
1964 : end do
1965 1080 : WRITE(unit_xml,'(a)')''
1966 : end do
1967 162 : WRITE(unit_xml,'(" </correction_force>")')
1968 162 : if (eff_pot%has_anharmonicsTerms)then
1969 72 : if (eff_pot%anharmonics_terms%has_elastic3rd) then
1970 0 : WRITE(unit_xml,'(" <elastic3rd units=""hartree"">")')
1971 0 : WRITE(unit_xml,'(6(E23.14))') (eff_pot%anharmonics_terms%elastic3rd(ii,:,:))
1972 0 : WRITE(unit_xml,'(" </elastic3rd>")')
1973 : end if
1974 72 : if (eff_pot%anharmonics_terms%has_elastic_displ) then
1975 0 : WRITE(unit_xml,'(" <correction_strain_force units=""hartree/bohrradius"">")')
1976 0 : do ia=1,eff_pot%crystal%natom
1977 0 : do mu=1,3
1978 0 : do nu=1,6
1979 : WRITE(unit_xml,'(e22.14)', advance="no")&
1980 0 : & (eff_pot%anharmonics_terms%elastic_displacement(ii,nu,mu,ia))
1981 : end do
1982 : end do
1983 0 : WRITE(unit_xml,'(a)')''
1984 : end do
1985 0 : WRITE(unit_xml,'(" </correction_strain_force>")')
1986 : end if
1987 72 : if (eff_pot%anharmonics_terms%has_strain_coupling) then
1988 0 : do irpt=1,eff_pot%anharmonics_terms%phonon_strain(ii)%nrpt
1989 0 : WRITE(unit_xml,'(" <correction_force_constant units=""hartree/bohrradius**2"">")')
1990 0 : WRITE(unit_xml,'(" <data>")')
1991 0 : do ia=1,eff_pot%crystal%natom
1992 0 : do mu=1,3
1993 0 : do ib=1,eff_pot%crystal%natom
1994 0 : do nu=1,3
1995 : WRITE(unit_xml,'(e22.14)', advance="no")&
1996 0 : & (eff_pot%anharmonics_terms%phonon_strain(ii)%atmfrc(mu,ia,nu,ib,irpt))
1997 : end do
1998 : end do
1999 0 : WRITE(unit_xml,'(a)')''
2000 : end do
2001 : end do
2002 0 : WRITE(unit_xml,'(" </data>")')
2003 0 : WRITE(unit_xml,'(" <cell>")')
2004 0 : WRITE(unit_xml,'(3(I4))') (eff_pot%anharmonics_terms%phonon_strain(ii)%cell(:,irpt))
2005 0 : WRITE(unit_xml,'(" </cell>")')
2006 0 : WRITE(unit_xml,'(" </correction_force_constant>")')
2007 : end do
2008 : end if!End if has_straincouplitn
2009 : end if!end Hasstrain_coupling
2010 162 : WRITE(unit_xml,'(" </strain_coupling>")')
2011 189 : jj = jj + 1
2012 : end do!end mu
2013 :
2014 27 : if(option /=1) WRITE(unit_xml,'("</System_definition>")')
2015 : ! Close file
2016 27 : CLOSE(unit_xml)
2017 :
2018 : end if!end option
2019 :
2020 : !Print the coefficients into XML file
2021 29 : new_file = .FALSE.
2022 29 : if (option==1 .or. option == 2 .or. option==4) then
2023 : ! Compute the name of the XML file
2024 22 : if(present(filename)) then
2025 22 : namefile=replace(trim(filename),".out","")
2026 10 : select case(option)
2027 : case(1)
2028 10 : new_file = .FALSE.
2029 10 : namefile=trim(namefile)//"_model.xml"
2030 : case(2)
2031 10 : new_file = .TRUE.
2032 10 : namefile=trim(namefile)//"_coeffs.xml"
2033 : case(4)
2034 2 : new_file = .TRUE.
2035 22 : namefile=trim(namefile)//"_coeffs.xml"
2036 : end select
2037 : else
2038 0 : namefile='coeffs.xml'
2039 : end if
2040 :
2041 22 : if(eff_pot%anharmonics_terms%ncoeff > 0) then
2042 : call polynomial_coeff_writeXML(eff_pot%anharmonics_terms%coefficients,&
2043 : & eff_pot%anharmonics_terms%ncoeff,namefile,unit=unit_xml,&
2044 15 : & newfile=new_file)
2045 : end if
2046 : end if!end option
2047 :
2048 29 : if(option==1)then
2049 : ! add the end of the file in the case option==1
2050 10 : open(unit=unit_xml,file=namefile,position="append")
2051 10 : WRITE(unit_xml,'("</System_definition>")')
2052 : ! Close file
2053 10 : CLOSE(unit_xml)
2054 : end if
2055 :
2056 29 : end subroutine effective_potential_writeXML
2057 : !!***
2058 :
2059 : !!****f* m_effective_potential/effective_potential_writeAbiInput
2060 : !! NAME
2061 : !! effective_potential_writeAbiInput
2062 : !!
2063 : !! FUNCTION
2064 : !! This routine print the effective potential into input of abinit
2065 : !! We can also apply a strain to the structure
2066 : !!
2067 : !! INPUTS
2068 : !! eff_pot<type(effective_potential_type)> = effective_potential datatype
2069 : !! filename = the name of input file
2070 : !! strain<(strain_type)> = optional,strain datatype if need to apply strain into rprim
2071 : !!
2072 : !! OUTPUT
2073 : !!
2074 : !! SOURCE
2075 :
2076 0 : subroutine effective_potential_writeAbiInput(eff_pot,filename,strain)
2077 :
2078 : !Arguments ------------------------------------
2079 : !scalars
2080 : type(strain_type),optional,intent(in) :: strain
2081 : character(len=fnlen),optional,intent(in) :: filename
2082 : !arrays
2083 : type(effective_potential_type), intent(in) :: eff_pot
2084 : !Local variables-------------------------------
2085 : !scalar
2086 : integer :: unit = 20
2087 : character(len=500) :: msg
2088 : character(len=fnlen) :: namefile
2089 : !arrays
2090 0 : real(dp) :: xred(3,eff_pot%crystal%natom)
2091 : type(strain_type) :: strain_tmp
2092 :
2093 : ! ************************************************************************
2094 :
2095 0 : if(present(strain)) then
2096 0 : strain_tmp = strain
2097 : else
2098 0 : call strain_init(strain_tmp)
2099 : end if
2100 :
2101 : ! try to open the file
2102 0 : if(present(filename)) then
2103 0 : namefile=filename
2104 : else
2105 0 : if(eff_pot%name /='') then
2106 0 : write(namefile,'(a)') 'structure_'//trim(eff_pot%name)
2107 : else
2108 0 : write(namefile,'(a)') 'structure'
2109 : end if
2110 0 : if (strain_tmp%name/='') then
2111 : write(namefile,'(a,a,a,a,a,a,a)') trim(namefile)//"_"//trim(strain_tmp%name)//"_"//&
2112 0 : & trim(itoa(strain_tmp%direction)),"_"//trim(ftoa(strain_tmp%delta))
2113 : end if
2114 :
2115 0 : namefile=trim(namefile)//".in"
2116 :
2117 : end if
2118 :
2119 0 : call isfile(namefile,'new')
2120 :
2121 0 : if (open_file(namefile,msg,unit=unit,form="formatted",status="new",action="write") /= 0) then
2122 0 : ABI_ERROR(msg)
2123 : end if
2124 :
2125 0 : write(msg,'(a,a,a,a)')ch10,&
2126 0 : & ' Generation of the input file in ',trim(namefile),ch10
2127 0 : call wrtout(ab_out,msg,'COLL')
2128 0 : call wrtout(std_out,msg,'COLL')
2129 :
2130 0 : write(unit,'("#Abinit Input for DFPT, this file contrains the keyword")')
2131 0 : write(unit,'("#To run DFPT calculation of")')
2132 0 : write(unit,'(a)') trim(eff_pot%name)
2133 0 : if (strain_tmp%direction /= 0) then
2134 0 : write(unit,'("# With a perturbation ")',advance="no")
2135 0 : write(unit,'(a)',advance="no") trim(strain_tmp%name)
2136 0 : write(unit,'(" in the direction : ")',advance="no")
2137 0 : write(unit,'(a)',advance='no') trim(itoa(strain_tmp%direction))
2138 0 : write(unit,'(" with the deformation : ")',advance="no")
2139 0 : write(unit,'(a)') trim(ftoa(strain_tmp%delta))
2140 : end if
2141 :
2142 0 : write(unit,'("")')
2143 0 : write(unit,'("ndtset 1 jdtset 1 2 3")')
2144 0 : write(unit,'("")')
2145 0 : write(unit,'("#DATASET1 GROUND STATE")')
2146 0 : write(unit,'("tolwfr1 = 1d-15")')
2147 0 : write(unit,'(" prtwf1 = 1")')
2148 0 : write(unit,'(" nline1 = 5")')
2149 :
2150 0 : write(unit,'("")')
2151 0 : write(unit,'("#DATASET2 DDK PERTURBATION")')
2152 0 : write(unit,'("getwfk2 = 1")')
2153 0 : write(unit,'(" iscf2 = -3")')
2154 0 : write(unit,'(" nline2 = 15")')
2155 0 : write(unit,'("nnsclo2 = 5")')
2156 0 : write(unit,'("kptopt2 = 2")')
2157 0 : write(unit,'(" nqpt2 = 1")')
2158 0 : write(unit,'(" qpt2 = 0 0 0 ")')
2159 0 : write(unit,'("rfelfd2 = 2")')
2160 0 : write(unit,'(" rfdir2 = 1 1 1 ")')
2161 0 : write(unit,'("tolwfr2 = 1.0d-20 ")')
2162 0 : write(unit,'(" prtwf2 = 1 ")')
2163 0 : write(unit,'(" ")')
2164 :
2165 0 : write(unit,'("#DATASET3 RF")')
2166 0 : write(unit,'(" getddk3 = 2")')
2167 0 : write(unit,'(" getwfk3 = 1")')
2168 0 : write(unit,'(" iscf3 = 7")')
2169 0 : write(unit,'(" kptopt3 = 2")')
2170 0 : write(unit,'(" nqpt3 = 1")')
2171 0 : write(unit,'(" qpt3 = 0 0 0")')
2172 0 : write(unit,'(" rfphon3 = 1")')
2173 0 : write(unit,'("rfatpol3 = 1 ")',advance='no')
2174 0 : write(unit,'(a)') itoa(eff_pot%crystal%natom)
2175 0 : write(unit,'(" rfelfd3 = 3")')
2176 0 : write(unit,'(" rfstrs3 = 3")')
2177 0 : write(unit,'(" rfdir3 = 1 1 1")')
2178 0 : write(unit,'(" tolvrs3 = 1.0d-8")')
2179 0 : write(unit,'("")')
2180 :
2181 0 : write(unit,'("#STRUCTURE")')
2182 0 : write(unit,'(" natom = ")',advance='no')
2183 0 : write(unit,'(a)') itoa(eff_pot%crystal%natom)
2184 0 : write(unit,'(" znucl =")',advance='no')
2185 0 : write(unit,'(10(F4.0))') (eff_pot%crystal%znucl)
2186 0 : write(unit,'("ntypat = ")',advance='no')
2187 0 : write(unit,'(a)') itoa(eff_pot%crystal%ntypat)
2188 0 : write(unit,'(" typat = ")',advance='no')
2189 0 : write(unit,'(10(I2))') (eff_pot%crystal%typat)
2190 0 : write(unit,'(" acell = 1 1 1")')
2191 0 : write(unit,'(" rprim ")')
2192 0 : write(unit,'(3(F20.10))') (matmul(eff_pot%crystal%rprimd,strain%strain))
2193 0 : write(unit,'(" xred ")')
2194 0 : call xcart2xred(eff_pot%crystal%natom,eff_pot%crystal%rprimd,eff_pot%crystal%xcart,xred)
2195 0 : write(unit,'(3(F15.10))') (xred)
2196 0 : write(unit,'(" ")')
2197 :
2198 0 : write(unit,'("#SCF")')
2199 0 : write(unit,'(" ecut = ")')
2200 0 : write(unit,'("pawecutdg = ")')
2201 0 : write(unit,'(" ecutsm = ")')
2202 0 : write(unit,'(" tolvrs = ")')
2203 0 : write(unit,'(" nband = ")')
2204 0 : write(unit,'(" ixc = ")')
2205 0 : write(unit,'(" occopt = ")')
2206 0 : write(unit,'(" nstep = ")')
2207 0 : write(unit,'(" kptopt = ")')
2208 0 : write(unit,'(" ngkpt = ")')
2209 0 : write(unit,'("")')
2210 0 : write(unit,'(" prtwf 0 prtden 0 prtdos 0")')
2211 :
2212 0 : close(unit)
2213 :
2214 0 : end subroutine effective_potential_writeAbiInput
2215 : !!***
2216 :
2217 :
2218 : !****f* m_effective_potential/effective_potential_evaluate
2219 : !!
2220 : !! NAME
2221 : !! effective_potential_evaluate
2222 : !!
2223 : !! FUNCTION
2224 : !! evaluate the energy, forces and stresses
2225 : !! of a structure with the effective potential.
2226 : !! For the structure you can only provide xred and rprimd,
2227 : !! this routine will compute the displacement and the strain.
2228 : !! Otherwise, you can provide the strain and the displacement
2229 : !!
2230 : !! INPUTS
2231 : !! eff_pot <type(effective_potential_type)> = effective potential datatype
2232 : !! natom = number of atoms in the supercell
2233 : !! rprimd(3,3) = dimensional primitive translations for real space
2234 : !! of the perturbed structure
2235 : !! xred(3,natom) = optional, reduced coordinates of the atoms of the structure
2236 : !! strain(6) = optional, strain apply wrt to the reference in eff_pot
2237 : !! displacement(3,natom) = optional, atomic displacement wrt to the reference in eff_pot
2238 : !! du_delta(6,3,natom) = optional, variation of displacement wrt to the strain
2239 : !! compute_anharmonic = optional, TRUE if the anhamornic part is needed
2240 : !! verbose = optional, FALSE for no verbose
2241 : !!
2242 : !! OUTPUTS
2243 : !! energy = energy of the structure (Ha)
2244 : !! fcart(3,natom) = forces in cartesian coordinates (Ha/Bohr)
2245 : !! gred(3,natom) = gradient wrt nuclear positions in reduced coordinates
2246 : !! strten(6) = stress tensor (Ha/Bohr^3)
2247 : !!
2248 : !! SOURCE
2249 :
2250 27038 : subroutine effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
2251 1922 : & displacement,du_delta,strain,xred,&
2252 : & compute_anharmonic,verbose,filename,elec_eval,efield_type,efield, &
2253 : & efield_lambda, efield_period,efield_phase,nefield,&
2254 : & efield_gmean,efield_gvel,efield_sigma,efield_background,time)
2255 :
2256 : !Arguments ------------------------------------
2257 : !scalars
2258 : integer, intent(in) :: natom
2259 : character(len=fnlen),optional,intent(in) :: filename
2260 : !array
2261 : type(effective_potential_type),intent(in) :: eff_pot
2262 : real(dp),intent(out) :: energy
2263 : real(dp),intent(out) :: fcart(3,natom)
2264 : real(dp),intent(out) :: gred(3,natom)
2265 : real(dp),intent(out) :: strten(6)
2266 : real(dp),intent(in) :: rprimd(3,3)
2267 : integer, intent(in),optional :: efield_type, nefield
2268 : real(dp),intent(in),optional :: xred(3,natom)
2269 : real(dp),intent(in),optional :: efield(3,*)
2270 : real(dp),intent(in),optional :: efield_background(3)
2271 : real(dp),intent(inout),optional :: efield_lambda(3,*)
2272 : real(dp),intent(in),optional :: efield_gmean(3)
2273 : real(dp),intent(in),optional :: efield_gvel(3)
2274 : real(dp),intent(in),optional :: efield_period
2275 : real(dp),intent(in),optional :: efield_phase(*)
2276 : real(dp),intent(in),optional :: efield_sigma
2277 : real(dp),intent(in),optional :: time
2278 : real(dp),intent(in),optional :: strain(6)
2279 : real(dp),intent(in),optional :: displacement(3,natom)
2280 : real(dp),intent(in),optional :: du_delta(6,3,natom)
2281 : logical,intent(in),optional :: verbose,compute_anharmonic,elec_eval
2282 : !Local variables-------------------------------
2283 : !scalar
2284 : integer :: alpha,ii,ia,mu,ncell,comm,natom_for_scup
2285 : real(dp):: energy_part,ucvol
2286 : logical :: has_strain,need_verbose
2287 : logical :: iam_master,need_anharmonic
2288 : logical :: err_eng, err_for
2289 : logical*1 :: update_dens
2290 : logical :: need_elec_eval
2291 : !integer :: icell,ierr,kk
2292 : integer, parameter:: master = 0
2293 : !array
2294 : type(strain_type) :: strain_t
2295 : integer :: sc_size(3)
2296 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2297 : real(dp) :: strain_tmp(6),strten_part(6)
2298 13519 : real(dp),allocatable :: disp_tmp(:,:)
2299 13519 : real(dp),allocatable :: du_delta_tmp(:,:,:)
2300 13519 : real(dp),allocatable :: fcart_part(:,:)
2301 13519 : real(dp),allocatable :: energy_coeff_part(:)
2302 13519 : real(dp),allocatable :: xcart(:,:)
2303 : !real(dp) :: ext_field(3),ext_field2(3),temp_pol(3)
2304 : character(len=500) :: msg
2305 : ! *************************************************************************
2306 :
2307 : !Hide SCALE-UP Variables
2308 : ABI_UNUSED(update_dens)
2309 :
2310 : !MPI variables
2311 13519 : comm = eff_pot%mpi_ifc%comm
2312 13519 : iam_master = (eff_pot%mpi_ifc%my_rank == master)
2313 :
2314 : ! Set variables
2315 13519 : ncell = eff_pot%supercell%ncells
2316 54076 : do ii = 1, 3
2317 54076 : sc_size(ii) = eff_pot%supercell%rlatt(ii,ii)
2318 : end do
2319 :
2320 13519 : natom_for_scup = natom
2321 :
2322 13519 : need_verbose = .TRUE.
2323 13519 : if(present(verbose)) then
2324 13519 : need_verbose = verbose
2325 : end if
2326 :
2327 13519 : need_anharmonic = .TRUE.
2328 13519 : if(present(compute_anharmonic))then
2329 3655 : need_anharmonic = compute_anharmonic
2330 : end if
2331 :
2332 13519 : need_elec_eval = .FALSE.
2333 13519 : if(present(elec_eval)) need_elec_eval = elec_eval
2334 :
2335 : ! Check some variables
2336 13519 : if (natom /= eff_pot%supercell%natom) then
2337 0 : write(msg,'(a,I7,a,I7,a)')' The number of atoms is not correct :',natom,&
2338 0 : & ' in argument istead of ',eff_pot%supercell%natom, ' in supercell'
2339 0 : ABI_ERROR(msg)
2340 : end if
2341 :
2342 13519 : if (present(displacement))then
2343 961 : if(size(displacement(1,:)) /= eff_pot%supercell%natom) then
2344 0 : write(msg,'(a,I7,a,I7,a)')' The number of atoms is not correct :',size(displacement(1,:)),&
2345 0 : & ' in displacement array instead of ',eff_pot%supercell%natom, ' in supercell'
2346 0 : ABI_ERROR(msg)
2347 : end if
2348 : end if
2349 13519 : if (present(du_delta))then
2350 961 : if(size(du_delta,3) /= eff_pot%supercell%natom) then
2351 0 : write(msg,'(a,I7,a,I7,a)')' The number of atoms is not correct :',size(du_delta,3),&
2352 0 : & ' in du_delta array instead of ',eff_pot%supercell%natom, ' in supercell'
2353 0 : ABI_ERROR(msg)
2354 : end if
2355 : end if
2356 54076 : do ii=1,3
2357 54076 : if(eff_pot%supercell%rlatt(ii,ii)<0.or.eff_pot%supercell%rlatt(ii,ii)>150)then
2358 : write(msg, '(a,i0,a,i2,a,a,a,i0,a)' )&
2359 0 : & 'eff_pot%supercell%rlatt(',ii,') is ',int(eff_pot%supercell%rlatt(ii,ii)),&
2360 0 : & ', which is lower than 0 or superior than 150.',ch10,'Action: correct ncell(',ii,').'
2361 0 : ABI_ERROR(msg)
2362 : end if
2363 : end do
2364 :
2365 13519 : if(need_verbose)then
2366 2056 : write(msg, '(a,a,a)' ) ch10,' Calculation of the energy with effective potential'
2367 2056 : call wrtout(ab_out,msg,'COLL')
2368 : end if
2369 :
2370 13519 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2371 : ! write(*,*) "---- UCVOL ----"
2372 : ! write(*,*) ucvol
2373 :
2374 : !--------------------------------------------
2375 : ! 1 - Set the perturbations and intialisation
2376 : !--------------------------------------------
2377 : ! Get strain
2378 13519 : has_strain = .FALSE.
2379 13519 : strain_tmp(:) = zero
2380 13519 : if (present(strain)) then
2381 961 : strain_tmp(:) = strain(:)
2382 : call strain_get(strain_t,mat_delta=reshape((/&
2383 : & strain_tmp(1),strain_tmp(6)/two,strain_tmp(5)/two,&
2384 : & strain_tmp(6)/two,strain_tmp(2),strain_tmp(4)/two,&
2385 9610 : & strain_tmp(5)/two,strain_tmp(4)/two,strain_tmp(3)/),(/3,3/)))
2386 961 : has_strain = .TRUE.
2387 : else
2388 : ! Compute the strain
2389 : call strain_get(strain_t,rprim=eff_pot%supercell%rprimd,&
2390 12558 : & rprim_def=rprimd,symmetrized=.FALSE.)
2391 12558 : if(need_verbose)then
2392 166536 : write(msg,'(80a)')('-',mu=1,80)
2393 2056 : call wrtout(std_out,msg,'COLL')
2394 2056 : call strain_print(strain_t)
2395 : end if
2396 12558 : if (strain_t%name /= "reference") then
2397 3676 : has_strain = .TRUE.
2398 14704 : do ii=1,3
2399 14704 : strain_tmp(ii) = strain_t%strain(ii,ii)
2400 : end do
2401 3676 : strain_tmp(4) = (strain_t%strain(2,3) + strain_t%strain(3,2))
2402 3676 : strain_tmp(5) = (strain_t%strain(3,1) + strain_t%strain(1,3))
2403 3676 : strain_tmp(6) = (strain_t%strain(2,1) + strain_t%strain(1,2))
2404 : else
2405 8882 : strain_tmp(:) = zero
2406 : end if
2407 : end if
2408 : ! Get displacement and the variation of the displacmeent wr to strain
2409 40557 : ABI_MALLOC(xcart,(3,natom))
2410 27038 : ABI_MALLOC(disp_tmp,(3,natom))
2411 40557 : ABI_MALLOC(du_delta_tmp,(6,3,natom))
2412 27038 : ABI_MALLOC(fcart_part,(3,natom))
2413 40557 : ABI_MALLOC(energy_coeff_part,(eff_pot%anharmonics_terms%ncoeff))
2414 36076719 : disp_tmp(:,:) = zero
2415 198361119 : du_delta_tmp(:,:,:) = zero
2416 13519 : if((.not.present(displacement).or..not.present(du_delta)).and.present(xred))then
2417 : ! Compute the displacement
2418 12558 : call xred2xcart(natom, rprimd, xcart, xred)
2419 : call effective_potential_getDisp(disp_tmp,du_delta_tmp,natom,rprimd,&
2420 : & eff_pot%supercell%rprimd,comm,xcart_hist=xcart,&
2421 : & xcart_ref=eff_pot%supercell%xcart,&
2422 : & compute_displacement = .not.present(displacement),&
2423 12558 : & compute_duDelta = .not.present(du_delta))
2424 : end if
2425 :
2426 : ! or set ftrom the arguments
2427 167279 : if(present(displacement)) disp_tmp(:,:) = displacement(:,:)
2428 859199 : if(present(du_delta))du_delta_tmp(:,:,:) = du_delta(:,:,:)
2429 : !Set to zero the outputs
2430 13519 : energy = zero
2431 36076719 : fcart(:,:) = zero
2432 13519 : strten(:) = zero
2433 :
2434 13519 : if(need_verbose)then
2435 166536 : write(msg, '(80a,2a)' ) ('-',mu=1,80),&
2436 168592 : & ch10,' Components of total energy (in Hartree) :'
2437 2056 : call wrtout(ab_out,msg,'COLL')
2438 2056 : call wrtout(std_out,msg,'COLL')
2439 : end if
2440 :
2441 : !------------------------------------
2442 : ! 2 - Transfert the reference values
2443 : !------------------------------------
2444 :
2445 : ! Set the value of the energy
2446 13519 : energy = eff_pot%energy * ncell
2447 :
2448 13519 : if(need_verbose)then
2449 2056 : write(msg, '(a,a,1ES24.16,a)' ) ch10,' Energy of the reference structure :',&
2450 4112 : & energy,' Hartree'
2451 2056 : call wrtout(ab_out,msg,'COLL')
2452 2056 : call wrtout(std_out,msg,'COLL')
2453 : end if
2454 :
2455 : !Set the value of the initial strees (should be zero if the system is relaxed)
2456 : !According to the original definition of the model, the reference should always be
2457 : ! a critical point of the PES so the first order derivative should not be take into account...
2458 : ! do ii=1,6
2459 : ! energy = energy + ncell * eff_pot%strten(ii)*strain_tmp(ii)
2460 : ! end do
2461 : ! strten(:) = ncell * eff_pot%strten(:)
2462 : !Set the value of the initial forces (should be zero if the system is relaxed)
2463 : ! ii = 1
2464 : ! This next part is not working, need to apply strain to the original fcart
2465 : ! do ia=1,eff_pot%supercell%natom
2466 : ! fcart(:,ia) = eff_pot%fcart(:,ii)
2467 : ! do mu=1,3
2468 : ! energy = energy + eff_pot%fcart(mu,ii)*disp_tmp(mu,ii)
2469 : ! end do
2470 : ! ii = ii + 1
2471 : ! if(ii > eff_pot%crystal%natom) ii = 1
2472 : ! end do
2473 :
2474 : !------------------------------------
2475 : ! 3 - Computation of the IFC part :
2476 : !------------------------------------
2477 :
2478 : energy_part = zero
2479 36076719 : fcart_part(:,:)= zero
2480 :
2481 : call harmonics_terms_evaluateIFC(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,:),disp_tmp,&
2482 : & energy_part,fcart_part,eff_pot%supercell%natom,&
2483 : & eff_pot%crystal%natom,eff_pot%mpi_ifc%my_ncell,&
2484 : & eff_pot%mpi_ifc%my_nrpt,eff_pot%mpi_ifc%my_atmrpt_index,&
2485 : & eff_pot%mpi_ifc%my_index_cells,sc_size,eff_pot%mpi_ifc%my_rpt,&
2486 13519 : & eff_pot%mpi_ifc%comm)
2487 :
2488 13519 : if(need_verbose)then
2489 2056 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the ifc part :',&
2490 4112 : & energy_part,' Hartree'
2491 2056 : call wrtout(ab_out,msg,'COLL')
2492 2056 : call wrtout(std_out,msg,'COLL')
2493 :
2494 2056 : if(abs(energy_part) < tol16 .and. eff_pot%anharmonics_terms%ncoeff == 0)then
2495 25 : write(msg, '(8a)' )ch10,&
2496 25 : & ' --- !WARNING!',ch10,&
2497 25 : & ' The harmonic part is negative, the simulation will diverge',ch10,&
2498 25 : & ' if the anharmonic part is not used',ch10,&
2499 50 : & ' ---'
2500 25 : call wrtout(std_out,msg,"COLL")
2501 : end if
2502 : end if
2503 :
2504 13519 : energy = energy + energy_part
2505 36076719 : fcart(:,:)= fcart(:,:) + fcart_part(:,:)
2506 :
2507 :
2508 : !----------------------------------------------------
2509 : ! 4 - Computation of the elastic part of the energy :
2510 : !----------------------------------------------------
2511 :
2512 13519 : energy_part = zero
2513 36076719 : fcart_part(:,:)= zero
2514 13519 : strten_part(:) = zero
2515 :
2516 : ! write(*,*) "--- STRTEN before harmonics evaluate --- "
2517 : ! write(*,*) strten(:)
2518 : call harmonics_terms_evaluateElastic(eff_pot%harmonics_terms%elastic_constants,disp_tmp,energy_part,&
2519 : & fcart_part,eff_pot%supercell%natom,eff_pot%crystal%natom,&
2520 13519 : & ncell,eff_pot%harmonics_terms%strain_coupling,strten_part,strain_tmp)
2521 :
2522 :
2523 13519 : if(has_strain.and.need_verbose)then
2524 975 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the elastic part :',&
2525 1950 : & energy_part,' Hartree'
2526 975 : call wrtout(ab_out,msg,'COLL')
2527 975 : call wrtout(std_out,msg,'COLL')
2528 : end if
2529 :
2530 13519 : energy = energy + energy_part
2531 36076719 : fcart(:,:) = fcart(:,:) + fcart_part(:,:)
2532 94633 : strten(:) = strten(:) + strten_part(:)
2533 : ! write(*,*) "--- STRTEN after harmonics evaluate --- "
2534 : ! write(*,*) strten(:)
2535 :
2536 : !------------------------------------
2537 : ! 5 - Treat 3rd order strain-coupling:
2538 : !------------------------------------
2539 13519 : if (need_anharmonic.and.eff_pot%has_anharmonicsTerms) then
2540 1104 : energy_part = zero
2541 1104 : strten_part(:) = zero
2542 648144 : fcart_part(:,:)= zero
2543 :
2544 : ! 1-Part due to the anharmonic strain
2545 : if(eff_pot%anharmonics_terms%has_elastic3rd.or.&
2546 1104 : & eff_pot%anharmonics_terms%has_elastic4th.or.&
2547 : & eff_pot%anharmonics_terms%has_elastic_displ)then
2548 : call anharmonics_terms_evaluateElastic(disp_tmp,energy_part,fcart_part,&
2549 : & eff_pot%supercell%natom,&
2550 : & eff_pot%crystal%natom,ncell,strten_part,strain_tmp,&
2551 : & elastic3rd=eff_pot%anharmonics_terms%elastic3rd,&
2552 : & elastic4th=eff_pot%anharmonics_terms%elastic4th,&
2553 0 : & elastic_displacement=eff_pot%anharmonics_terms%elastic_displacement)
2554 :
2555 0 : if(need_verbose)then
2556 0 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the anharmonic elastic part :',&
2557 0 : & energy_part,' Hartree'
2558 0 : call wrtout(ab_out,msg,'COLL')
2559 0 : call wrtout(std_out,msg,'COLL')
2560 : end if
2561 0 : energy = energy + energy_part
2562 1104 : strten(:) = strten(:) + strten_part(:)
2563 : end if
2564 :
2565 : ! 2-Part due to the strain-phonon coupling
2566 1104 : if (eff_pot%anharmonics_terms%has_strain_coupling) then
2567 : energy_part = zero
2568 0 : strten_part(:) = zero
2569 0 : fcart_part(:,:)= zero
2570 :
2571 : call anharmonics_terms_evaluateIFCStrainCoupling(eff_pot%anharmonics_terms%phonon_strain,&
2572 : & disp_tmp,energy_part,fcart_part,&
2573 : & eff_pot%supercell%natom,&
2574 : & eff_pot%crystal%natom,&
2575 : & sc_size,strain_tmp,&
2576 : & strten_part,eff_pot%mpi_ifc%my_cells,&
2577 : & eff_pot%mpi_ifc%my_ncell,&
2578 : & eff_pot%mpi_ifc%my_index_cells,&
2579 0 : & eff_pot%mpi_ifc%comm)
2580 :
2581 0 : if(need_verbose)then
2582 0 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the 3rd (strain-phonon coupling):',&
2583 0 : & energy_part,' Hartree'
2584 0 : call wrtout(ab_out,msg,'COLL')
2585 0 : call wrtout(std_out,msg,'COLL')
2586 : end if
2587 0 : energy = energy + energy_part
2588 0 : fcart = fcart + fcart_part
2589 0 : strten = strten + strten_part
2590 : end if
2591 :
2592 : end if
2593 :
2594 : !----------------------------------
2595 : ! 6 - Treat polynomial coefficient:
2596 : !----------------------------------
2597 12172 : if(need_anharmonic.and.eff_pot%anharmonics_terms%ncoeff > zero)then
2598 10265 : energy_part = zero
2599 35095225 : fcart_part(:,:) = zero
2600 10265 : strten_part(:) = zero
2601 : call polynomial_coeff_evaluate(eff_pot%anharmonics_terms%coefficients,disp_tmp,&
2602 : & energy_part,energy_coeff_part,fcart_part,eff_pot%supercell%natom,&
2603 : & eff_pot%crystal%natom,eff_pot%anharmonics_terms%ncoeff,&
2604 : & sc_size,strain_tmp,strten_part,eff_pot%mpi_coeff%my_ncell,&
2605 : & eff_pot%mpi_coeff%my_index_cells,eff_pot%mpi_coeff%comm,&
2606 10671 : & filename=filename)
2607 :
2608 10265 : if(need_verbose)then
2609 706 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the fitted coefficient :',&
2610 1412 : & energy_part,' Hartree'
2611 706 : call wrtout(ab_out,msg,'COLL')
2612 706 : call wrtout(std_out,msg,'COLL')
2613 : end if
2614 10265 : energy = energy + energy_part
2615 35095225 : fcart(:,:) = fcart(:,:) + fcart_part(:,:)
2616 71855 : strten(:) = strten(:) + strten_part(:)
2617 : end if
2618 :
2619 : !---------------------------------
2620 : ! 7 - Compute confinement
2621 : !---------------------------------
2622 13519 : if(eff_pot%confinement%need_confinement) then
2623 :
2624 : energy_part = zero
2625 :
2626 : call polynomial_conf_evaluate(disp_tmp,eff_pot%confinement%cutoff_disp,energy_part,&
2627 : & eff_pot%confinement%factor_disp,&
2628 : & eff_pot%confinement%factor_strain,fcart_part,strain_tmp,&
2629 : & eff_pot%confinement%cutoff_strain,strten_part,&
2630 : & eff_pot%confinement%power_disp,eff_pot%confinement%power_strain,&
2631 : & eff_pot%mpi_coeff%my_cells,&
2632 : & eff_pot%supercell%natom,eff_pot%crystal%natom,&
2633 : & eff_pot%mpi_coeff%my_ncell,eff_pot%mpi_coeff%my_index_cells,&
2634 0 : & eff_pot%mpi_coeff%comm)
2635 :
2636 0 : energy = energy + energy_part
2637 :
2638 0 : if(abs(energy_part) > tol10 .and. need_verbose )then
2639 0 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the confinement part :',&
2640 0 : & energy_part,' Hartree'
2641 0 : call wrtout(ab_out,msg,'COLL')
2642 0 : call wrtout(std_out,msg,'COLL')
2643 : end if
2644 : end if
2645 :
2646 : !-------------------------------------------
2647 : !!! External_ELECTRIC_FILED
2648 : ! 7.5 - Compute Forces and energies from external electric field
2649 : !------------------------------------------
2650 13519 : if(present(efield_type))then
2651 9864 : if (efield_type /= 0) then
2652 : call calculate_forces_efield(eff_pot,energy_part,eff_pot%mpi_coeff%comm,fcart_part,natom,&
2653 : & disp_tmp,xcart,efield_type,efield,efield_lambda,efield_period,&
2654 60 : & efield_phase,efield_gmean,efield_gvel,efield_sigma,efield_background,time,nefield)
2655 60 : if(need_verbose)then
2656 60 : write(msg, '(a,1ES24.16,a)' ) ' Energy from electric field is :',&
2657 120 : & energy_part,' Hartree'
2658 60 : call wrtout(ab_out,msg,'COLL')
2659 60 : call wrtout(std_out,msg,'COLL')
2660 :
2661 : end if
2662 60 : energy = energy + energy_part
2663 19260 : fcart = fcart + fcart_part
2664 : endif
2665 : end if
2666 : !-------------------------------------------
2667 : ! 8 - Compute electronic Part with SCALE-UP
2668 : !------------------------------------------
2669 :
2670 13519 : err_eng = .FALSE.
2671 13519 : err_for = .FALSE.
2672 :
2673 13519 : energy_part = 0
2674 36076719 : fcart_part = 0
2675 :
2676 13519 : if(need_elec_eval)then
2677 : #if defined DEV_MS_SCALEUP
2678 : update_dens = .TRUE.
2679 : write(msg,'(a)') ' wohoo i was here and call scale-up now!---STILL WOHOOO!---'
2680 : call wrtout(ab_out,msg,'COLL')
2681 : call wrtout(std_out,msg,'COLL')
2682 : err_eng = global_calculate_energy(energy_part,strain_tmp,disp_tmp,natom_for_scup)
2683 : err_for = global_calculate_forces(fcart_part,strain_tmp,disp_tmp,natom_for_scup,update_dens)
2684 : write(msg,'(a,1ES24.16,a)') 'The Energy of the electronic model is: ', energy_part,' unit?'
2685 : call wrtout(ab_out,msg,'COLL')
2686 : call wrtout(std_out,msg,'COLL')
2687 : #endif
2688 :
2689 0 : energy = energy + energy_part
2690 0 : do ia = 1, natom
2691 0 : fcart(:,ia) = fcart(:,ia) + fcart_part(:,ia)
2692 : end do
2693 : endif
2694 :
2695 :
2696 :
2697 : !-----------------------------------
2698 : ! 9 - Add stress due to forces on
2699 : ! atoms and variation of disp
2700 : ! due to strain
2701 : !-----------------------------------
2702 :
2703 : ! multiply forces by -1
2704 36076719 : fcart = -1 * fcart
2705 :
2706 :
2707 :
2708 :
2709 : !Calculatei
2710 13519 : strten_part(:) = zero
2711 9029319 : do ia = 1,eff_pot%supercell%natom
2712 36076719 : do mu = 1,3
2713 198347600 : do alpha=1,6
2714 189331800 : strten_part(alpha) = strten_part(alpha) - fcart(mu,ia) * du_delta_tmp(alpha,mu,ia)
2715 : end do
2716 : end do
2717 : end do
2718 94633 : strten(:) = strten(:) + strten_part(:)
2719 :
2720 : ! write(*,*) "--- STRTEN after du/deta application --- "
2721 : ! write(*,*) strten(:)
2722 :
2723 : !---------------------------------
2724 : ! 10 - Renormalize stress
2725 : !---------------------------------
2726 :
2727 : ! divide stess tensor by ucvol
2728 94633 : do alpha=1,6
2729 94633 : if (alpha < 4) then
2730 40557 : strten(alpha) = strten(alpha) * (1 + strain_tmp(alpha))/ ucvol
2731 : else
2732 40557 : strten(alpha) = strten(alpha) * (1 - strain_tmp(alpha)**2)/ ucvol
2733 : end if
2734 : end do
2735 :
2736 : ! write(*,*) "--- STRTEN after /ucvol --- "
2737 : ! write(*,*) strten(:)
2738 : ! Redistribute the residuale of the forces
2739 13519 : call effective_potential_distributeResidualForces(eff_pot,fcart,eff_pot%supercell%natom)
2740 :
2741 13519 : call fcart2gred(fcart,gred,rprimd,natom)
2742 :
2743 : !------------------------------------
2744 : ! 11 - Final Print:
2745 : !------------------------------------
2746 :
2747 13519 : if(need_verbose)then
2748 2056 : write(msg, '(2a,es21.14)' ) ch10,&
2749 4112 : & ' >>>>>>>>> Etotal= ',energy
2750 2056 : call wrtout(ab_out,msg,'COLL')
2751 2056 : call wrtout(std_out,msg,'COLL')
2752 :
2753 2056 : write(msg,'(2a,1p,e15.7,a)') ch10,' Unit cell volume ucvol=',ucvol+tol10,' bohr^3'
2754 2056 : call wrtout(std_out, msg,'COLL')
2755 :
2756 168592 : write(msg, '(a,80a,3a)' ) ch10,('-',mu=1,80),ch10,&
2757 170648 : & ' Cartesian components of stress tensor (hartree/bohr^3)'
2758 2056 : call wrtout(ab_out,msg,'COLL')
2759 2056 : call wrtout(std_out, msg,'COLL')
2760 :
2761 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2762 2056 : & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4)
2763 2056 : call wrtout(ab_out,msg,'COLL')
2764 2056 : call wrtout(std_out, msg,'COLL')
2765 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2766 2056 : & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5)
2767 2056 : call wrtout(ab_out,msg,'COLL')
2768 2056 : call wrtout(std_out, msg,'COLL')
2769 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2770 2056 : & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6)
2771 2056 : call wrtout(ab_out,msg,'COLL')
2772 2056 : call wrtout(std_out, msg,'COLL')
2773 :
2774 : ! Also output the pressure (minus one third the trace of the stress
2775 : ! tensor.
2776 2056 : write(msg, '(a,a,es12.4,a)' ) ch10,&
2777 2056 : & '-Cartesian components of stress tensor (GPa) [Pressure=',&
2778 4112 : & -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]'
2779 :
2780 2056 : call wrtout(ab_out,msg,'COLL')
2781 2056 : call wrtout(std_out, msg,'COLL')
2782 :
2783 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2784 2056 : & '- sigma(1 1)=',strten(1)*HaBohr3_GPa,&
2785 4112 : & ' sigma(3 2)=',strten(4)*HaBohr3_GPa
2786 2056 : call wrtout(ab_out,msg,'COLL')
2787 2056 : call wrtout(std_out, msg,'COLL')
2788 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2789 2056 : & '- sigma(2 2)=',strten(2)*HaBohr3_GPa,&
2790 4112 : & ' sigma(3 1)=',strten(5)*HaBohr3_GPa
2791 2056 : call wrtout(ab_out,msg,'COLL')
2792 2056 : call wrtout(std_out, msg,'COLL')
2793 : write(msg, '(a,1p,e16.8,a,1p,e16.8,a)' ) &
2794 2056 : & '- sigma(3 3)=',strten(3)*HaBohr3_GPa,&
2795 4112 : & ' sigma(2 1)=',strten(6)*HaBohr3_GPa
2796 2056 : call wrtout(ab_out,msg,'COLL')
2797 2056 : call wrtout(std_out, msg,'COLL')
2798 166536 : write(msg, '(80a,a)' ) ('-',mu=1,80),ch10
2799 2056 : call wrtout(ab_out,msg,'COLL')
2800 2056 : call wrtout(std_out,msg,'COLL')
2801 : end if
2802 :
2803 :
2804 :
2805 13519 : ABI_FREE(xcart)
2806 13519 : ABI_FREE(disp_tmp)
2807 13519 : ABI_FREE(du_delta_tmp)
2808 13519 : ABI_FREE(fcart_part)
2809 13519 : ABI_FREE(energy_coeff_part)
2810 :
2811 27999 : end subroutine effective_potential_evaluate
2812 :
2813 :
2814 : !****f* m_effective_potential/calculate_forces_efield
2815 : !!
2816 : !! NAME
2817 : !! calculate_forces_efield
2818 : !!
2819 : !! FUNCTION
2820 : !! Compute the action of the electric field
2821 : !! different values of efield_type allow different functional forms
2822 : !! including spatial inhomogeneous fields and time dependent fields
2823 : !! INPUTS
2824 : !!
2825 : !! effpot for the structural information about supercells etc
2826 : !! comm parallel communication
2827 : !! natom number of atoms on the supercell
2828 : !! disp_tmp instantaneous atomic displacements
2829 : !! xcart cartesian coordinates of the atoms
2830 : !! time instantaneous time
2831 : !! efield_type different electric field types
2832 : !! 1- homogeneous 2- AC Fields
2833 : !! 3- cos-like spatial modulation 4-cos-like spatial and time modulated
2834 : !! 5- Gaussian field 6-Double cos-like modulation
2835 : !! efield_lambda spatial modulation on cos-like fields
2836 : !! efield_period Time modulation on AC fields
2837 : !! efield_phase Phase shift for cos-like
2838 : !! efield_gmean Center of the Gaussian field
2839 : !! efield_gvel Velocity of the Gaussian field
2840 : !! efield_sigma Std deviation of the efield
2841 : !! efield_background Homogeneous background
2842 : !!
2843 : !! OUTPUT
2844 : !!
2845 : !! energy_part to update the energy contribution of the electric field
2846 : !! fcart_part to update the contribution to the forces
2847 : !!
2848 : !!
2849 60 : subroutine calculate_forces_efield(eff_pot,energy_part,comm,fcart_part,natom,disp_tmp,xcart,&
2850 : & efield_type,efield,efield_lambda,efield_period,efield_phase,&
2851 : & efield_gmean,efield_gvel,efield_sigma,efield_background,time,nefield)
2852 : !Arguments ------------------------------------
2853 : integer, intent(in) :: natom,comm
2854 : type(effective_potential_type),intent(in) :: eff_pot
2855 : real(dp),intent(inout) :: energy_part
2856 : real(dp),intent(inout) :: fcart_part(3,natom)
2857 : integer, intent(in) :: efield_type
2858 : integer, intent(in), optional :: nefield
2859 : real(dp),intent(in),optional :: efield(3,*)
2860 : real(dp),intent(in),optional :: efield_background(3)
2861 : real(dp),intent(inout),optional :: efield_lambda(3,*)
2862 : real(dp),intent(in),optional :: efield_gmean(3)
2863 : real(dp),intent(in),optional :: efield_gvel(3)
2864 : real(dp),intent(in),optional :: efield_period
2865 : real(dp),intent(in),optional :: efield_phase(*)
2866 : real(dp),intent(in),optional :: efield_sigma
2867 : real(dp),intent(in),optional :: time
2868 : real(dp),intent(in),optional :: disp_tmp(3,natom),xcart(3,natom)
2869 : !Arguments ------------------------------------
2870 : !Internal variables----------------------------
2871 : integer :: ii,ia
2872 : real(dp):: argumento,argumento2
2873 : integer :: icell,ierr,kk,ll
2874 : !array
2875 : real(dp) :: ext_field(3),ext_field2(3), temp_pol(3),efield_kk(3),pos_gcenter(3),ext_field_bck(3)
2876 : !Internal variables----------------------------
2877 70 : select case (efield_type)
2878 : case(1)
2879 10 : ext_field = zero
2880 10 : temp_pol = zero
2881 10 : energy_part = zero
2882 3210 : fcart_part(:,:) = zero
2883 10 : if (present(efield)) then
2884 40 : ext_field = -1 * efield(:,1) ! Use first field by default
2885 : endif
2886 10 : if (present(efield_background)) then
2887 40 : ext_field = ext_field - efield_background
2888 : endif
2889 :
2890 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2891 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2892 970 : do ia = 1,eff_pot%crystal%natom
2893 800 : kk = ii + ia
2894 : !We compute the polarization from the atomic displacements
2895 2400 : temp_pol = temp_pol + matmul(disp_tmp(:,kk),&
2896 23200 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2897 : !We compute the forces as the Born effective charges times the electric field
2898 20160 : fcart_part(:,kk) = matmul(ext_field, eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2899 : end do
2900 : end do
2901 : !We compute the energy as the Polarization times the electric field
2902 40 : energy_part = DOT_PRODUCT(ext_field(:),temp_pol(:))
2903 10 : call xmpi_sum(energy_part, comm, ierr)
2904 10 : call xmpi_sum(fcart_part, comm, ierr)
2905 : case (2)
2906 10 : if (present(efield) .and. present(nefield)) then
2907 40 : ext_field = -1 * efield(:,1) ! Use first field by default
2908 : endif
2909 10 : if (present(efield_background)) then
2910 40 : ext_field_bck = -1 * efield_background
2911 : endif
2912 10 : temp_pol = zero
2913 10 : energy_part = zero
2914 3210 : fcart_part(:,:) = zero
2915 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2916 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2917 970 : do ia = 1, eff_pot%crystal%natom
2918 800 : kk = ii + ia
2919 : !We compute the electric field value at the given time
2920 3200 : efield_kk(:)=cos(-two_pi/efield_period*time+efield_phase(1))*ext_field(:)+ext_field_bck(:)
2921 800 : temp_pol = matmul(disp_tmp(:,kk),&
2922 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2923 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2924 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
2925 : end do
2926 : end do
2927 10 : call xmpi_sum(energy_part, comm, ierr)
2928 10 : call xmpi_sum(fcart_part , comm, ierr)
2929 : case (3)
2930 10 : if (present(efield)) then
2931 40 : ext_field =-1* efield(:, 1)
2932 : endif
2933 10 : if (present(efield_background)) then
2934 40 : ext_field_bck =-1* efield_background
2935 : endif
2936 10 : temp_pol = zero
2937 10 : energy_part = zero
2938 3210 : fcart_part(:,:) = zero
2939 40 : do ll=1,3
2940 : !We update the value of the spatial modulation to convey with PBC in case
2941 : !the lattice has been updated
2942 40 : if(efield_lambda(ll,1)/=0.0 .and. abs(efield_lambda(ll,1))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
2943 10 : efield_lambda(ll,1)=sign(eff_pot%supercell%rprimd(ll,ll)/nint(eff_pot%supercell%rprimd(ll,ll)/efield_lambda(ll,1)),efield_lambda(ll,1))
2944 : endif
2945 : enddo
2946 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2947 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2948 970 : do ia = 1, eff_pot%crystal%natom
2949 800 : kk = ii + ia
2950 800 : argumento=0.0
2951 3200 : do ll=1,3
2952 3200 : if (efield_lambda(ll,1)/=0.0) then
2953 800 : argumento=argumento+two_pi*xcart(ll,kk)/efield_lambda(ll,1)
2954 : endif
2955 : enddo
2956 3200 : efield_kk(:)=cos(argumento+efield_phase(1))*ext_field(:)+ext_field_bck(:)
2957 800 : temp_pol = matmul(disp_tmp(:,kk),&
2958 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2959 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2960 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
2961 : end do
2962 : end do
2963 10 : call xmpi_sum(energy_part, comm, ierr)
2964 10 : call xmpi_sum(fcart_part , comm, ierr)
2965 : case (4)
2966 10 : if (present(efield)) then
2967 40 : ext_field =-1*efield(:, 1)
2968 : endif
2969 10 : if (present(efield_background)) then
2970 40 : ext_field_bck =-1* efield_background
2971 : endif
2972 10 : temp_pol = zero
2973 10 : energy_part = zero
2974 3210 : fcart_part(:,:) = zero
2975 40 : do ll=1,3
2976 40 : if(efield_lambda(ll,1)/=0.0 .and. abs(efield_lambda(ll,1))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
2977 10 : efield_lambda(ll,1)=sign(eff_pot%supercell%rprimd(ll,ll)/nint(eff_pot%supercell%rprimd(ll,ll)/efield_lambda(ll,1)),efield_lambda(ll,1))
2978 : endif
2979 : enddo
2980 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2981 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2982 970 : do ia = 1, eff_pot%crystal%natom
2983 800 : kk = ii + ia
2984 800 : argumento=0.0
2985 3200 : do ll=1,3
2986 3200 : if (efield_lambda(ll,1)/=0.0) then
2987 800 : argumento=argumento+two_pi*xcart(ll,kk)/efield_lambda(ll,1)
2988 : endif
2989 : enddo
2990 3200 : efield_kk(:)=cos(argumento-two_pi/efield_period*time+efield_phase(1))*ext_field(:)+ext_field_bck(:)
2991 800 : temp_pol = matmul(disp_tmp(:,kk),&
2992 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2993 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2994 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
2995 : end do
2996 : end do
2997 10 : call xmpi_sum(energy_part, comm, ierr)
2998 10 : call xmpi_sum(fcart_part , comm, ierr)
2999 : case (5)
3000 10 : if (present(efield) .and. present(nefield)) then
3001 40 : ext_field = -1 * efield(:,1) ! Use first field by default
3002 : endif
3003 10 : if (present(efield_background)) then
3004 40 : ext_field_bck =-1* efield_background
3005 : endif
3006 10 : temp_pol = zero
3007 10 : energy_part = zero
3008 3210 : fcart_part(:,:) = zero
3009 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
3010 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
3011 970 : do ia = 1, eff_pot%crystal%natom
3012 800 : kk = ii + ia
3013 3200 : pos_gcenter(:)=xcart(:,kk)-efield_gmean(:)-efield_gvel(:)*time
3014 3200 : do ll=1,3
3015 : !PBC for the Gaussian field
3016 : pos_gcenter(ll)=pos_gcenter(ll)-eff_pot%supercell%rprimd(ll,ll)*&
3017 3200 : &nint(pos_gcenter(ll)/eff_pot%supercell%rprimd(ll,ll))
3018 : enddo
3019 5600 : efield_kk(:)=EXP(-DOT_PRODUCT(pos_gcenter(:),pos_gcenter(:))/efield_sigma/efield_sigma)*ext_field(:)+ext_field_bck(:)
3020 800 : temp_pol = matmul(disp_tmp(:,kk),&
3021 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3022 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3023 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
3024 : end do
3025 : end do
3026 10 : call xmpi_sum(energy_part, comm, ierr)
3027 10 : call xmpi_sum(fcart_part , comm, ierr)
3028 : case (6)
3029 10 : if (present(efield) .and. present(nefield)) then
3030 10 : if (nefield >= 1) then
3031 40 : ext_field = -1 * efield(:,1)
3032 40 : ext_field2 = -1 * efield(:,2)
3033 : endif
3034 : endif
3035 10 : if (present(efield_background)) then
3036 40 : ext_field_bck = -1 * efield_background
3037 : endif
3038 10 : temp_pol = zero
3039 10 : energy_part = zero
3040 3210 : fcart_part(:,:) = zero
3041 40 : do ll=1,3
3042 40 : if(present(efield_lambda)) then
3043 30 : if(efield_lambda(ll,1)/=0.0 .and. abs(efield_lambda(ll,1))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
3044 10 : efield_lambda(ll,1)=sign(eff_pot%supercell%rprimd(ll,ll)/nint(eff_pot%supercell%rprimd(ll,ll)/efield_lambda(ll,1)),efield_lambda(ll,1))
3045 : endif
3046 30 : if(present(nefield) .and. nefield >= 2) then
3047 30 : if(efield_lambda(ll,2)/=0.0 .and. abs(efield_lambda(ll,2))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
3048 10 : efield_lambda(ll,2)=sign(eff_pot%supercell%rprimd(ll,ll)/nint(eff_pot%supercell%rprimd(ll,ll)/efield_lambda(ll,2)),efield_lambda(ll,2))
3049 : endif
3050 : endif
3051 : endif
3052 : enddo
3053 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
3054 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
3055 970 : do ia = 1, eff_pot%crystal%natom
3056 800 : kk = ii + ia
3057 800 : argumento=0.0
3058 800 : argumento2=0.0
3059 3200 : do ll=1,3
3060 2400 : if (present(efield_lambda) .and. efield_lambda(ll,1)/=0.0) then
3061 800 : argumento=argumento+two_pi*xcart(ll,kk)/efield_lambda(ll,1)
3062 : endif
3063 : if (present(efield_lambda) .and. present(nefield) &
3064 3200 : & .and. nefield >= 2 .and. efield_lambda(ll,2)/=0.0) then
3065 800 : argumento2=argumento2+two_pi*xcart(ll,kk)/efield_lambda(ll,2)
3066 : endif
3067 : enddo
3068 3200 : efield_kk(:)=cos(argumento+efield_phase(1))*ext_field(:)
3069 800 : if (present(nefield) .and. nefield >= 2) then
3070 3200 : efield_kk(:)=efield_kk(:)+cos(argumento2+efield_phase(2))*ext_field2
3071 : endif
3072 800 : if (present(efield_background)) then
3073 3200 : efield_kk(:)=efield_kk(:)+ext_field_bck(:)
3074 : endif
3075 800 : temp_pol = matmul(disp_tmp(:,kk),&
3076 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3077 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3078 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
3079 : end do
3080 : end do
3081 10 : call xmpi_sum(energy_part, comm, ierr)
3082 120 : call xmpi_sum(fcart_part , comm, ierr)
3083 : end select
3084 : !!***
3085 : !!***
3086 : !!Subroutine Inh Elec Field Ends-----
3087 : !!***
3088 : !!***
3089 60 : end subroutine calculate_forces_efield
3090 :
3091 :
3092 : !!***
3093 :
3094 : !****f* m_effective_potential/effective_potential_getDisp
3095 : !!
3096 : !! NAME
3097 : !! effective_potential_getDisp
3098 : !!
3099 : !! FUNCTION
3100 : !! Compute cartesian coordinates of the displacment
3101 : !! between two configurations.
3102 : !! You can provide reduce or cartesian coordinates for both perturbed structure
3103 : !! and reference structure
3104 : !!
3105 : !! INPUTS
3106 : !! natom = number of atoms
3107 : !! rprimd_hist = dimensional primitive translations for real space of the perturbed structure
3108 : !! rprimd_ref = dimensional primitive translations for real space of the reference structure
3109 : !! comm = mpi communicator
3110 : !! xcart_hist(3,natom) = optional, cartesian coordinates of the atoms in the perturbed structure
3111 : !! xred_hist(3,natom) = optional, reduced coordinates of the atoms in the perturbed structure
3112 : !! xred_ref(3,natom) = optional, reduced coordinates of the atoms in the reference structure
3113 : !! xcart_ref(3,natom) = optional, cartesian coordinates of the atoms in the reference structure
3114 : !! compute_displacement= optional, flag to compute the displacement array (default is true)
3115 : !! compute_duDelta = optional, flag to compute the du_delta array (default is true)
3116 : !! Be careful, if you specify compute_displacement=.false.,
3117 : !! you need to provide to this routine the good displacement array!
3118 : !!
3119 : !!
3120 : !! OUTPUT
3121 : !! displacement(3,natom) = cartesian atomic displacement between two configurations (Bohr)
3122 : !! du_delta(6,3,natom) = variation of the displacmeent wr to strain
3123 : !!
3124 : !! SOURCE
3125 :
3126 27038 : subroutine effective_potential_getDisp(displacement,du_delta,natom,rprimd_hist,rprimd_ref,comm,&
3127 26077 : & xcart_hist,xred_hist,xred_ref,xcart_ref,compute_displacement,&
3128 : & compute_duDelta)
3129 :
3130 : !Arguments ------------------------------------
3131 : !scalars
3132 : integer, intent(in) :: natom,comm
3133 : logical,optional,intent(in) :: compute_displacement,compute_duDelta
3134 : !array
3135 : real(dp),intent(in) :: rprimd_ref(3,3),rprimd_hist(3,3)
3136 : real(dp),intent(out) :: displacement(3,natom),du_delta(6,3,natom)
3137 : real(dp),intent(in),optional :: xred_hist(3,natom),xcart_hist(3,natom)
3138 : real(dp),intent(in),optional :: xred_ref(3,natom),xcart_ref(3,natom)
3139 : !Local variables-------------------------------
3140 : !scalar
3141 : integer :: ii,ia,ib,mu
3142 : integer :: ierr,nproc,my_rank,natom_alone,my_natom
3143 : character(len=500) :: msg
3144 : logical :: has_strain
3145 : logical :: need_displacement,need_duDelta
3146 : !array
3147 : integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
3148 13519 : integer,allocatable :: my_atoms(:)
3149 : type(strain_type) :: strain
3150 13519 : real(dp),allocatable :: xcart_hist_tmp(:,:),xcart_ref_tmp(:,:),xred_ref_tmp(:,:)
3151 : real(dp) :: strain_inv(3,3),strain_inv_u(3),strain_tmp(3,3)
3152 : ! *************************************************************************
3153 :
3154 13519 : if (.not.(present(xred_ref).or.present(xcart_ref))) then
3155 : write(msg, '(3a)' )&
3156 0 : & 'You need at least give xcart_ref or xred_ref '
3157 0 : ABI_ERROR(msg)
3158 : end if
3159 :
3160 13519 : if (.not.(present(xred_hist).or.present(xcart_hist))) then
3161 : write(msg, '(3a)' )&
3162 0 : & 'You need at least give xcart_hist or xred_hist '
3163 0 : ABI_ERROR(msg)
3164 : end if
3165 :
3166 13519 : need_duDelta = .TRUE.
3167 13519 : need_displacement = .TRUE.
3168 13519 : if(present(compute_duDelta)) need_duDelta = compute_duDelta
3169 13519 : if(present(compute_displacement)) need_displacement = compute_displacement
3170 13519 : if(.not.need_duDelta .and. .not.need_displacement) return
3171 :
3172 40557 : ABI_MALLOC(xcart_hist_tmp,(3,natom))
3173 27038 : ABI_MALLOC(xcart_ref_tmp,(3,natom))
3174 27038 : ABI_MALLOC(xred_ref_tmp,(3,natom))
3175 :
3176 : !--------------------------------------------
3177 : ! 0 - Set the MPI
3178 : !--------------------------------------------
3179 13519 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3180 13519 : natom_alone = mod(natom,nproc)
3181 13519 : my_natom = int(aint(real(natom,sp)/(nproc)))
3182 13519 : if(my_rank >= (nproc-natom_alone)) then
3183 200 : my_natom = my_natom + 1
3184 : end if
3185 40557 : ABI_MALLOC(my_atoms,(my_natom))
3186 8969039 : my_atoms = 0
3187 8969039 : do ii=1,my_natom
3188 8969039 : if(my_rank >= (nproc-natom_alone))then
3189 : my_atoms(ii)=(int(aint(real(natom,sp)/nproc)))*(my_rank)+&
3190 8500 : & (my_rank - (nproc-natom_alone)) + ii
3191 : else
3192 8947020 : my_atoms(ii)=(my_natom)*(my_rank) + ii
3193 : end if
3194 : end do
3195 :
3196 : !--------------------------------------------
3197 : ! 1 - Get the strain for this step
3198 : !--------------------------------------------
3199 13519 : has_strain = .FALSE.
3200 13519 : call strain_get(strain,rprim=rprimd_ref,rprim_def=rprimd_hist)
3201 13519 : if (strain%name /= "reference") then
3202 18548 : has_strain = .TRUE.
3203 60281 : strain_tmp = strain%strain
3204 4637 : strain_tmp(1,1) = strain_tmp(1,1) + 1.0
3205 4637 : strain_tmp(2,2) = strain_tmp(2,2) + 1.0
3206 4637 : strain_tmp(3,3) = strain_tmp(3,3) + 1.0
3207 : ! get (1+eta)^-1
3208 4637 : call matr3inv(strain_tmp,strain_inv)
3209 : else if (strain%name == "reference") then
3210 8882 : strain_inv(:,:) = zero
3211 8882 : strain_inv(1,1) = 1.0
3212 8882 : strain_inv(2,2) = 1.0
3213 8882 : strain_inv(3,3) = 1.0
3214 : end if
3215 : ! write(*,*) "---- STRAIN ----"
3216 : ! do ii = 1,3
3217 : ! write(*,*) strain%strain(ii,:)
3218 : ! enddo
3219 : ! write(*,*) "---- 1+STRAIN inv ----"
3220 : ! do ii = 1,3
3221 : ! write(*,*) strain_inv(ii,:)
3222 : ! enddo
3223 :
3224 : ! fill the history position
3225 13519 : if(present(xcart_hist)) then
3226 35921998 : xcart_hist_tmp(:,:) = xcart_hist(:,:)
3227 : else
3228 961 : call xred2xcart(natom, rprimd_hist, xcart_hist_tmp, xred_hist)
3229 : end if
3230 :
3231 : ! Fill the reference position and change the cartesian coordinates
3232 : ! if the rprimd is different
3233 13519 : if(has_strain) then
3234 4637 : if(present(xcart_ref)) then
3235 4637 : call xcart2xred(natom, rprimd_ref, xcart_ref, xred_ref_tmp)
3236 4637 : call xred2xcart(natom, rprimd_hist, xcart_ref_tmp, xred_ref_tmp)
3237 : else
3238 0 : call xred2xcart(natom, rprimd_hist, xcart_ref_tmp, xred_ref)
3239 : end if
3240 : else
3241 8882 : if(present(xcart_ref)) then
3242 34705746 : xcart_ref_tmp(:,:) = xcart_ref(:,:)
3243 : else
3244 0 : call xred2xcart(natom, rprimd_ref, xcart_ref_tmp, xred_ref)
3245 : end if
3246 : end if
3247 :
3248 : ! Compute displacement
3249 13519 : if(need_displacement)then
3250 36076719 : displacement(:,:) = zero
3251 : ! write(*,*) "----- DISP -----"
3252 9029319 : do ii = 1, natom
3253 36076719 : displacement(:,ii) = xcart_hist_tmp(:,ii) - xcart_ref_tmp(:,ii)
3254 : ! write(*,*) displacement(:,ii)
3255 : end do
3256 : end if
3257 :
3258 : ! Du_Delta after Equation A4 in 2017 Paper Carlos
3259 198361119 : du_delta=zero
3260 8969039 : do ia=1,my_natom
3261 8955520 : ib = my_atoms(ia)
3262 : !Calc (1+eta)^-1 * disp(ib)
3263 116421760 : strain_inv_u = MATMUL(strain_inv,displacement(:,ib))
3264 : ! write(*,*) "--- strain_inv_u ---"
3265 : ! write(*,*) strain_inv_u
3266 : ! fill du_delta_e
3267 62702159 : do ii = 1,6
3268 223888000 : do mu = 1,3
3269 161199360 : if(alpha(ii)==mu)then
3270 53733120 : du_delta(ii,mu,ib) = du_delta(ii,mu,ib) + half * strain_inv_u(beta(ii))
3271 : end if
3272 214932480 : if(beta(ii)==mu)then
3273 53733120 : du_delta(ii,mu,ib) = du_delta(ii,mu,ib) + half * strain_inv_u(alpha(ii))
3274 : end if
3275 : enddo
3276 : enddo
3277 : enddo
3278 13519 : call xmpi_sum(du_delta , comm, ierr)
3279 :
3280 13519 : ABI_FREE(my_atoms)
3281 13519 : ABI_FREE(xcart_hist_tmp)
3282 13519 : ABI_FREE(xcart_ref_tmp)
3283 13519 : ABI_FREE(xred_ref_tmp)
3284 :
3285 27038 : end subroutine effective_potential_getDisp
3286 : !!***
3287 :
3288 : !****f* m_effective_potential/effective_potential_distributeResidualForces
3289 : !!
3290 : !! NAME
3291 : !! effective_potential_distributeResidualForces
3292 : !!
3293 : !! FUNCTION
3294 : !! Distribute the residual forces in a weighted manner
3295 : !!
3296 : !! INPUTS
3297 : !! natom = number of atoms
3298 : !! eff_pot = effective potential structure
3299 : !!
3300 : !! OUTPUT
3301 : !! fcart = forces in cartesian coordinates
3302 : !!
3303 : !! SOURCE
3304 :
3305 13519 : subroutine effective_potential_distributeResidualForces(eff_pot,fcart,natom)
3306 :
3307 : !Arguments ------------------------------------
3308 : !scalars
3309 : integer, intent(in) :: natom
3310 : !array
3311 : type(effective_potential_type),intent(in) :: eff_pot
3312 : real(dp),intent(inout) :: fcart(3,natom)
3313 : !Local variables-------------------------------
3314 : !scalar
3315 : real(dp):: mass_ia,sum_mass
3316 : integer :: ia
3317 : !array
3318 : real(dp):: sum_f(3)
3319 :
3320 : ! *************************************************************************
3321 :
3322 9029319 : sum_f(1) = sum(fcart(1,:))
3323 9029319 : sum_f(2) = sum(fcart(2,:))
3324 9029319 : sum_f(3) = sum(fcart(3,:))
3325 13519 : sum_mass = zero
3326 :
3327 9029319 : do ia=1,natom
3328 9029319 : sum_mass = sum_mass + eff_pot%crystal%amu(eff_pot%supercell%typat(ia))
3329 : end do
3330 :
3331 9029319 : do ia=1,natom
3332 9015800 : mass_ia = eff_pot%crystal%amu(eff_pot%supercell%typat(ia))
3333 36076719 : fcart(:,ia) = fcart(:,ia) - (mass_ia/sum_mass) * sum_f(:)
3334 : end do
3335 :
3336 :
3337 13519 : end subroutine effective_potential_distributeResidualForces
3338 : !!***
3339 :
3340 : !****f* m_effective_potential/effective_potential_writeAnhHead
3341 : !!
3342 : !! NAME
3343 : !! effective_potential_writeAnhHead
3344 : !!
3345 : !! FUNCTION
3346 : !! Write Header of anharmonic_energy_terms.out file
3347 : !!
3348 : !! INPUTS
3349 : !! natom = number of atoms
3350 : !! eff_pot = effective potential structure
3351 : !!
3352 : !! OUTPUT
3353 : !!
3354 : !! SOURCE
3355 :
3356 1 : subroutine effective_potential_writeAnhHead(ncoeff,filename,anh_terms)
3357 :
3358 : !Arguments ------------------------------------
3359 : !scalars
3360 : integer, intent(in) :: ncoeff
3361 : !Strings/Characters
3362 : character(len=fnlen) :: filename
3363 : !array
3364 : type(anharmonics_terms_type ),intent(in) :: anh_terms
3365 : !Local variables-------------------------------
3366 : !scalar
3367 : integer :: icoeff,unit_out
3368 : !Strings/Characters
3369 : character(len=fnlen) :: name_file
3370 : character(len=200):: term_name
3371 : !array
3372 :
3373 : ! *************************************************************************
3374 :
3375 : ! Marcus: if wanted: analyze anharmonic terms of effective potential &&
3376 : ! and print anharmonic contribution to file anharmonic_energy_terms.out
3377 : ! Open File and write header
3378 1 : name_file=trim(filename)//'_anharmonic_terms_energy.dat'
3379 1 : unit_out = get_unit()
3380 1 : open(unit=unit_out,file=name_file,status='replace',form='formatted')
3381 1 : write(unit_out,*) '#---------------------------------------------#'
3382 1 : write(unit_out,*) '# Anharmonic Terms Energy Contribution #'
3383 1 : write(unit_out,*) '#---------------------------------------------#'
3384 1 : write(unit_out,*) ''
3385 1 : write(unit_out,'(A,I5)') 'Number of Terms: ', ncoeff
3386 1 : write(unit_out,*) ''
3387 1 : write(unit_out,'(A)') 'Terms Names'
3388 61 : do icoeff=1,ncoeff
3389 60 : term_name = anh_terms%coefficients(icoeff)%name
3390 61 : write(unit_out,'(I5,A,A)') icoeff,' ',trim(term_name)
3391 : enddo
3392 1 : write(unit_out,*) ''
3393 1 : write(unit_out,'(A)',advance='no') 'Cycle/Terms'
3394 61 : do icoeff=1,ncoeff
3395 61 : if(icoeff<ncoeff)then
3396 59 : write(unit_out,'(I5)',advance='no') icoeff
3397 : else
3398 1 : write(unit_out,'(I5)',advance='yes') icoeff
3399 : endif
3400 : enddo
3401 : !close(unit_out)
3402 :
3403 1 : end subroutine effective_potential_writeAnhHead
3404 :
3405 :
3406 : !AM_EXPERIMENTAL SECTION
3407 : !!****f* m_effective_potential/equal
3408 : !! NAME
3409 : !! equal
3410 : !!
3411 : !! FUNCTION
3412 : !! compare two effective potential
3413 : !!
3414 : !! INPUTS
3415 : !! e1<type(effective_potential_type)> = effective_potential datatype
3416 : !! e2<type(effective_potential_type)> = effective_potential datatype
3417 : !!
3418 : !! OUTPUT
3419 : !!
3420 : !! SOURCE
3421 :
3422 0 : pure function effective_potential_compare(e1,e2) result (res)
3423 :
3424 : !Arguments ------------------------------------
3425 : type(effective_potential_type), intent(in) :: e1,e2
3426 : logical :: res
3427 : ! *************************************************************************
3428 0 : res = .false.
3429 : if(e1%crystal%natom==e2%crystal%natom.and.&
3430 : & e1%harmonics_terms%ifcs%nrpt==e2%harmonics_terms%ifcs%nrpt.and.&
3431 : & e1%crystal%ntypat==e2%crystal%ntypat.and.&
3432 : & e1%harmonics_terms%nqpt==e2%harmonics_terms%nqpt.and.&
3433 0 : & abs(e1%energy-e2%energy)<tol16.and.&
3434 : & abs(e1%crystal%ucvol-e2%crystal%ucvol)<tol16) then
3435 0 : res = .true.
3436 : end if
3437 :
3438 0 : end function effective_potential_compare
3439 : !!***
3440 :
3441 : ! !****f* m_effective_potential/effective_potential_effpot2ddb
3442 : ! !!
3443 : ! !! NAME
3444 : ! !! effective_potential_effpot2ddb
3445 : ! !!
3446 : ! !! FUNCTION
3447 : ! !! Convert eff_pot into ddb datatype
3448 : ! !!
3449 : ! !! INPUTS
3450 : ! !! eff_pot = effective potential datatype
3451 : ! !!
3452 : ! !! OUTPUT
3453 : ! !! ddb = ddb with all information
3454 : ! !!
3455 : ! !! SOURCE
3456 :
3457 : ! subroutine effective_potential_effpot2ddb(ddb,crystal,eff_pot,ncell,nph1l,option,qph1l)
3458 :
3459 :
3460 : ! !Arguments ------------------------------------
3461 : ! !scalars
3462 : ! integer,intent(in) :: nph1l,option
3463 : ! !array
3464 : ! integer,intent(in) :: ncell(3)
3465 : ! real(dp),intent(in):: qph1l(3,nph1l)
3466 : ! type(effective_potential_type),intent(inout) :: eff_pot
3467 : ! type(ddb_type),intent(out) :: ddb
3468 : ! type(crystal_t),intent(out) :: crystal
3469 : ! !Local variables-------------------------------
3470 : ! !scalar
3471 : ! integer :: ii,jj,msym
3472 : ! real(dp):: ucvol
3473 :
3474 : ! ! type(anaddb_dataset_type) :: inp
3475 : ! !array
3476 : ! real(dp) :: gmet(3,3),rmet(3,3)
3477 : ! real(dp) :: gprimd(3,3),rprimd(3,3)
3478 : ! real(dp),allocatable :: xred(:,:)
3479 : ! ! character :: title(eff_pot%crystal%ntypat)
3480 : ! integer,allocatable :: symrel(:,:,:),symafm(:)
3481 : ! real(dp),allocatable :: tnons(:,:)
3482 :
3483 : ! ! *************************************************************************
3484 :
3485 : ! ! Number of 2dte blocks in present object
3486 : ! ! integer,allocatable :: flg(:,:)
3487 : ! ! flg(msize,nblok)
3488 : ! ! flag to indicate presence of a given block
3489 : ! ! integer,allocatable :: typ(:)
3490 : ! ! typ(nblok)
3491 : ! ! type of each block - ddk, dde, phonon etc...
3492 : ! ! real(dp),allocatable :: amu(:)
3493 : ! ! amu(ntypat)
3494 : ! ! mass of the atoms (atomic mass unit)
3495 : ! ! real(dp),allocatable :: nrm(:,:)
3496 : ! ! nrm(3,nblok)
3497 : ! ! norm of the q-points for each block - can be 0 to indicate a direction of approach to gamma
3498 : ! ! real(dp),allocatable :: qpt(:,:)
3499 : ! ! qpt(9,nblok)
3500 : ! ! q-point vector in reciprocal space (reduced lattice coordinates) for each block
3501 : ! ! real(dp),allocatable :: val(:,:,:)
3502 : ! ! val(2,msize,nblok)
3503 : ! ! values of the second energy derivatives in each block
3504 :
3505 : ! ! Useless value
3506 : ! ddb%nblok = -1
3507 :
3508 : ! !option = 1 just print ddb for 1 1 1 cell
3509 : ! if(option==1) then
3510 : ! ! Compute different matrices in real and reciprocal space, also
3511 : ! ! checks whether ucvol is positive.
3512 : ! call metric(gmet,gprimd,-1,rmet,eff_pot%crystal%rprimd,ucvol)
3513 :
3514 : ! ! Convert to rprimd
3515 : ! do ii=1,3
3516 : ! do jj=1,3
3517 : ! rprimd(ii,jj)=eff_pot%crystal%rprimd(ii,jj)
3518 : ! end do
3519 : ! end do
3520 :
3521 : ! ! Obtain reciprocal space primitive transl g from inverse trans of r
3522 : ! ! (Unlike in abinit, gprim is used throughout ifc; should be changed, later)
3523 : ! call matr3inv(rprimd,gprimd)
3524 :
3525 : ! ! transfert basic values
3526 : ! ddb%natom = eff_pot%crystal%natom
3527 : ! ddb%mpert = ddb%natom+6
3528 : ! ddb%msize = 3*ddb%mpert*3*ddb%mpert;
3529 : ! ddb%ntypat = eff_pot%crystal%ntypat
3530 : ! ddb%occopt = 3 ! default value
3531 : ! ddb%prtvol = 0 ! default value
3532 : ! ddb%rprim = rprimd ! dimensioless real space primitive vectors
3533 : ! ddb%gprim = gprimd ! dimensioless reciprocal space primitive vectors
3534 : ! ddb%acell = one
3535 : ! msym = 1
3536 : ! ! Setup crystal type
3537 : ! ABI_MALLOC(xred,(3,ddb%natom))
3538 : ! ! call xcar2xred(ddb%natom,eff_pot%crystal%rprimd,eff_pot%crystal%xcart,xred)
3539 : ! !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
3540 : ! ABI_MALLOC(symrel,(3,3,msym))
3541 : ! ABI_MALLOC(symafm,(msym))
3542 : ! ABI_MALLOC(tnons,(3,msym))
3543 :
3544 : ! ! call crystal_init(ddb%amu,crystal,1,ddb%natom,size(eff_pot%crystal%znucl),eff_pot%crystal%ntypat,1,&
3545 : ! !& eff_pot%crystal%rprimd,eff_pot%crystal%typat,xred,eff_pot%crystal%znucl,&
3546 : ! !& eff_pot%crystal%znucl,0,.FALSE.,.FALSE.,title)!,&
3547 : ! !& symrel=symrel,tnons=tnons,symafm=symafm)
3548 : ! ! call crystal_print(crystal)
3549 : ! ! stop
3550 : ! !TEST_AM
3551 : ! ABI_FREE(symrel)
3552 : ! ABI_FREE(symafm)
3553 : ! ABI_FREE(tnons)
3554 :
3555 : ! ABI_FREE(xred)
3556 :
3557 : ! else if (option==2) then
3558 : ! ! Compute different matrices in real and reciprocal space, also
3559 : ! ! checks whether ucvol is positive.
3560 : ! call metric(gmet,gprimd,-1,rmet,eff_pot%supercell%rprimd,ucvol)
3561 :
3562 : ! ! Convert to rprim (dimensionless)
3563 : ! do ii=1,3
3564 : ! do jj=1,3
3565 : ! rprimd(ii,jj)=eff_pot%supercell%rprimd(ii,jj)
3566 : ! end do
3567 : ! end do
3568 :
3569 : ! ! Obtain reciprocal space primitive transl g from inverse trans of r
3570 : ! ! (Unlike in abinit, gprim is used throughout ifc; should be changed, later)
3571 : ! call matr3inv(rprimd,gprimd)
3572 :
3573 : ! ! transfert basic values
3574 : ! ddb%natom = eff_pot%supercell%natom
3575 : ! ddb%ntypat = eff_pot%crystal%ntypat
3576 : ! ddb%mpert = ddb%natom+6
3577 : ! ddb%msize = 3*ddb%mpert*3*ddb%mpert;
3578 : ! ddb%occopt = 3 ! default value
3579 : ! ddb%prtvol = 0 ! default value
3580 : ! ddb%rprim = rprimd ! dimensioless real space primitive vectors
3581 : ! ddb%gprim = gprimd ! dimensioless reciprocal space primitive vectors
3582 : ! ddb%acell = one
3583 :
3584 : ! end if
3585 : ! !TEST_AM
3586 : ! !write(std_out,*)"natom ",ddb%natom
3587 : ! !write(std_out,*)"ntypat",ddb%ntypat
3588 : ! !write(std_out,*)"mpert",ddb%mpert
3589 : ! !write(std_out,*)"msize",ddb%msize
3590 : ! !write(std_out,*)"occopt",ddb%occopt
3591 : ! !write(std_out,*)"prtvol",ddb%prtvol
3592 : ! !write(std_out,*)"rprim",ddb%rprim
3593 : ! !write(std_out,*)"gprim",ddb%gprim
3594 : ! !write(std_out,*)"acell",ddb%acell
3595 : ! !TEST_AM
3596 :
3597 : ! end subroutine effective_potential_effpot2ddb
3598 : ! !!***
3599 :
3600 :
3601 : ! !****f* m_effective_potential/effective_potential_printPDOS
3602 : ! !!
3603 : ! !! NAME
3604 : ! !! effective_potential_printPDOS
3605 : ! !!
3606 : ! !! FUNCTION
3607 : ! !! Apply the acoustic sum rule on the effective potential
3608 : ! !!
3609 : ! !! INPUTS
3610 : ! !! eff_pot = effective potential datatype
3611 : ! !! option = 0 (default) do nothing
3612 : ! !! = 1 print PHFRQ for specific qgrid (need nph1l and qph1l)
3613 : ! !! = 2 print PHFRQ for supercell (q=gamma) (need nph1l and qph1l and ncell)
3614 : ! !! OUTPUT
3615 : ! !! eff_pot
3616 : ! !!
3617 : ! !! SOURCE
3618 :
3619 : ! subroutine effective_potential_printPDOS(eff_pot,filename,ncell,nph1l,option,qph1l)
3620 :
3621 : ! !Arguments ------------------------------------
3622 : ! !scalars
3623 : ! integer,intent(in) :: nph1l,option
3624 : ! !array
3625 : ! integer,intent(in) :: ncell(3)
3626 : ! real(dp),intent(in):: qph1l(3,nph1l)
3627 : ! type(effective_potential_type),intent(inout) :: eff_pot
3628 : ! character(len=fnlen),intent(in) :: filename
3629 : ! !Local variables-------------------------------
3630 : ! !scalar
3631 : ! ! integer :: lenstr
3632 : ! ! character(len=strlen) :: string
3633 : ! !array
3634 : ! type(crystal_t) :: Crystal
3635 : ! ! type(anaddb_dataset_type) :: inp
3636 : ! type(ddb_type) :: ddb
3637 : ! ! type(asrq0_t) :: asrq0
3638 :
3639 : ! ! *************************************************************************
3640 :
3641 : ! if (option > 0) then
3642 :
3643 : ! ! First: transfer into ddb datatype:
3644 : ! call effective_potential_effpot2ddb(ddb,Crystal,eff_pot,ncell,nph1l,option,qph1l)
3645 :
3646 : ! ! Setup fake anaddb_dataset
3647 : ! ! string = ''
3648 : ! ! lenstr = 0
3649 : ! ! call invars9(inp,lenstr,ddb%natom,string)
3650 : ! ! fill it with multibinit_dataset values
3651 : ! ! inp%prt_ifc = 1
3652 : ! ! inp%ifcflag = 1
3653 : ! ! inp%qph1l = qph1l
3654 : ! ! inp%nph1l = nph1l
3655 :
3656 : ! ! ! In case the interatomic forces are not calculated, the
3657 : ! ! ! ASR-correction (asrq0%d2asr) has to be determined here from the Dynamical matrix at Gamma.
3658 : ! ! if (inp%ifcflag == 0) then
3659 : ! ! asrq0 = ddb_get_asrq0(ddb, inp%asr, inp%rfmeth, crystal%xcart)
3660 : ! ! end if
3661 :
3662 : ! ! !MG: Note that I'm passing xmpi_comm_self here.
3663 : ! ! call mkphbs(eff_pot%harmonics_terms%ifcs,Crystal,inp,ddb,asrq0,filename,xmpi_comm_self)
3664 :
3665 : ! ! call asrq0_free(asrq0)
3666 :
3667 : ! end if
3668 :
3669 : ! end subroutine effective_potential_printPDOS
3670 : ! !!***
3671 :
3672 : !****f* m_effective_potential/effective_potential_computeGradient
3673 : !!
3674 : !! NAME
3675 : !! effective_potential_computeGradient
3676 : !!
3677 : !! FUNCTION
3678 : !! Compute finate differences on forces to compute dynmical matrix
3679 : !! at gamma for supercell
3680 : !!
3681 : !! INPUTS
3682 : !! eff_pot = effective potential datatype
3683 : !!
3684 : !! OUTPUT
3685 : !! dynmat = ddb with all information
3686 : !!
3687 : !! SOURCE
3688 :
3689 0 : subroutine effective_potential_computeGradient(delta,fcart_out,eff_pot,natom,ncell,option,comm)
3690 :
3691 : !Arguments ------------------------------------
3692 : !scalars
3693 : integer,intent(in) :: natom,option,comm
3694 : real(dp),intent(in) :: delta
3695 : !array
3696 : integer,intent(in) :: ncell(3)
3697 : type(effective_potential_type),intent(inout) :: eff_pot
3698 : real(dp),intent(out) :: fcart_out(3,natom)
3699 : !Local variables-------------------------------
3700 : !scalar
3701 : character(len=500) :: msg
3702 : integer :: ia,ib,ii,mu,nu,npt
3703 : real(dp):: delt,energy
3704 : !array
3705 : real(dp):: strten(6)
3706 0 : real(dp),allocatable :: disp(:,:),diff(:)
3707 0 : real(dp),allocatable :: gred(:,:),fcart(:,:),xred(:,:)
3708 :
3709 : ! *************************************************************************
3710 :
3711 : !Do Some check
3712 0 : do ii=1,3
3713 0 : if(abs(eff_pot%supercell%qphon(ii)-ncell(ii)) > tol16)then
3714 0 : call effective_potential_setSupercell(eff_pot,comm,ncell)
3715 : end if
3716 : end do
3717 :
3718 0 : write(msg,'(a,(80a),3a)') ch10,('-',ii=1,80),ch10,' Generation of the dynmical matrix by ',&
3719 0 : & 'finite differences'
3720 : ! call wrtout(ab_out,msg,'COLL')
3721 : ! call wrtout(std_out,msg,'COLL')
3722 :
3723 0 : npt = 3 !default
3724 0 : select case (option)
3725 : case (1)
3726 : ! write(msg,'(2a)') ch10,' Finite differences on 1 points '
3727 0 : npt = 2
3728 : case (2)
3729 : ! write(msg,'(2a)') ch10,' Finite differences on 3 points '
3730 0 : npt = 3
3731 : case (3)
3732 : ! write(msg,'(2a)') ch10,' Finite differences on 5 points '
3733 0 : npt = 5
3734 : end select
3735 :
3736 : ! call wrtout(ab_out,msg,'COLL')
3737 : ! call wrtout(std_out,msg,'COLL')
3738 :
3739 : ! Allocation of forces arrays
3740 :
3741 0 : ABI_MALLOC(disp,(3,natom))
3742 0 : ABI_MALLOC(diff,(npt))
3743 0 : ABI_MALLOC(gred,(3,natom))
3744 0 : ABI_MALLOC(fcart,(3,natom))
3745 0 : ABI_MALLOC(xred,(3,natom))
3746 :
3747 0 : fcart_out = zero
3748 :
3749 : call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3750 0 : & eff_pot%supercell%xcart,xred)
3751 :
3752 0 : do ia=1,eff_pot%supercell%natom
3753 0 : do mu=1,3
3754 0 : diff = zero
3755 0 : do ii=1,npt
3756 0 : delt = (-(npt/2+1)+ii) * delta
3757 0 : disp = zero
3758 0 : disp(mu,ia) = delt * eff_pot%supercell%rprimd(mu,mu)
3759 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,&
3760 : & strten,natom,eff_pot%supercell%rprimd,&
3761 : & displacement=disp,&
3762 0 : & compute_anharmonic=.FALSE.,verbose=.false.)
3763 :
3764 : ! diff(ii,:,:) = gred(:,:)
3765 0 : diff(ii) = energy
3766 : end do
3767 :
3768 0 : select case (option)
3769 : case (1)
3770 0 : fcart_out(mu,ia) = (diff(1)-diff(2)) / (delta)
3771 : case (2)
3772 0 : fcart_out(mu,ia) = (diff(3)-diff(1)) / (2*delta)
3773 : case (3)
3774 0 : fcart_out(mu,ia) = (-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta)
3775 : end select
3776 : end do
3777 : end do
3778 :
3779 : !TEST_AM
3780 : !Write the phonon into ddb format wavevector
3781 : !write(999, '(a,3es16.8,f6.1)' )' qpt',real((/0,0,0/),dp),1.0
3782 : !Write the matrix elements
3783 0 : do ib=1,eff_pot%supercell%natom
3784 : do nu=1,3
3785 : do ia=1,eff_pot%supercell%natom
3786 : do mu=1,3
3787 : ! write(999,'(4i4,2d22.14)')nu,ib,mu,ia,&
3788 : ! & dynmat(1,mu,ia,nu,ib),dynmat(2,mu,ia,nu,ib)
3789 : end do
3790 : end do
3791 : end do
3792 : end do
3793 : !TEST_AM
3794 :
3795 : ! Deallocation of arrays
3796 0 : ABI_FREE(disp)
3797 0 : ABI_FREE(diff)
3798 0 : ABI_FREE(gred)
3799 0 : ABI_FREE(fcart)
3800 0 : ABI_FREE(xred)
3801 :
3802 :
3803 0 : end subroutine effective_potential_computeGradient
3804 : !!***
3805 :
3806 : !!****f* m_effective_potential/effective_potential_checkDEV
3807 : !!
3808 : !! NAME
3809 : !! effective_potential_checkDEV
3810 : !!
3811 : !! FUNCTION
3812 : !! Routine for develloper Check by finite differences the equations in
3813 : !! effective_potential_evaluate need to provide HIST file, so you need to
3814 : !! activate the fit_process or bound_process to activate the reading of the HIST
3815 : !!
3816 : !! INPUTS
3817 : !! eff_pot<type(effective_potential)> = effective potential
3818 : !! hist<type(abihist)> = The history of the MD
3819 : !! natom = number of atom
3820 : !! ntime = number of time in the hist
3821 : !!
3822 : !! OUTPUT
3823 : !!
3824 : !! SOURCE
3825 :
3826 0 : subroutine effective_potential_checkDEV(eff_pot,hist,natom,ntime)
3827 :
3828 : !Arguments ------------------------------------
3829 : !scalars
3830 : integer, intent(in) :: natom,ntime
3831 : !arrays
3832 : type(effective_potential_type),intent(in) :: eff_pot
3833 : type(abihist),intent(in) :: hist
3834 : !Local variables-------------------------------
3835 : !scalar
3836 : integer :: ii,jj,ia,mu,npt,istep
3837 : ! integer :: ifirst
3838 : real(dp):: energy,delt,delta,ucvol
3839 : !arrays
3840 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),mat_def(3,3),identity(3,3)
3841 : real(dp):: strten(6),rprimd(3,3)
3842 : real(dp):: rprimd_def(3,3),rprimd_ref(3,3),deltalist(5)
3843 : real(dp):: strain(6),diff(5)
3844 0 : real(dp),allocatable :: xred(:,:),fcart(:,:),gred(:,:),disp(:,:),disp_red(:,:),du_delta(:,:,:)
3845 : integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
3846 : character(len=500) :: msg
3847 :
3848 : ! *************************************************************************
3849 :
3850 : !Do some checks
3851 0 : if(ntime /= hist%mxhist)then
3852 0 : write(msg,'(a)')'ntime is not correct'
3853 0 : ABI_BUG(msg)
3854 : end if
3855 :
3856 0 : if(natom /= size(hist%xred,2)) then
3857 0 : write(msg,'(a)')'natom is not correct'
3858 0 : ABI_BUG(msg)
3859 : end if
3860 :
3861 :
3862 0 : ABI_MALLOC(xred,(3,natom))
3863 0 : ABI_MALLOC(fcart,(3,natom))
3864 0 : ABI_MALLOC(gred,(3,natom))
3865 0 : ABI_MALLOC(disp,(3,natom))
3866 0 : ABI_MALLOC(disp_red,(3,natom))
3867 0 : ABI_MALLOC(du_delta,(6,3,natom))
3868 0 : xred = zero
3869 :
3870 : !option 1 => set the reference for the test
3871 : ! call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3872 : !& eff_pot%supercell%xcart,xred)
3873 : ! rprimd = eff_pot%supercell%rprimd
3874 :
3875 : !option 2 => set a specific step for the test
3876 0 : istep = 4
3877 0 : xred = hist%xred(:,:,istep)
3878 0 : rprimd = hist%rprimd(:,:,istep)
3879 :
3880 : rprimd_ref = eff_pot%supercell%rprimd
3881 0 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
3882 :
3883 0 : npt=5
3884 0 : delta = 0.001
3885 0 : deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/)
3886 : strain = zero
3887 :
3888 0 : do ia=1,natom
3889 0 : do mu=1,3
3890 0 : write(std_out,*) "atm: ",ia," dir: ",mu
3891 0 : do ii=1,npt
3892 0 : delt = deltalist(ii)
3893 :
3894 : ! Get the initial displacement
3895 : call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3896 : & eff_pot%supercell%rprimd,1,xred_hist=xred,&
3897 : & xcart_ref=eff_pot%supercell%xcart,&
3898 0 : & compute_displacement = .true.,compute_duDelta = .true.)
3899 :
3900 : ! Add the delta
3901 0 : call xcart2xred(natom, rprimd, disp, disp_red)
3902 0 : disp_red(mu,ia) = disp_red(mu,ia) + delt
3903 0 : call xred2xcart(natom, rprimd, disp, disp_red)
3904 :
3905 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3906 : & xred=xred,du_delta=du_delta,&
3907 0 : & displacement=disp,compute_anharmonic=.true.,verbose=.false.)
3908 0 : diff(ii) = energy
3909 :
3910 : end do
3911 :
3912 : ! Get the initial displacement
3913 : call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3914 : & eff_pot%supercell%rprimd,1,xred_hist=xred,&
3915 : & xcart_ref=eff_pot%supercell%xcart,&
3916 0 : & compute_displacement = .true.,compute_duDelta = .true.)
3917 :
3918 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3919 : & xred=xred,du_delta=du_delta,&
3920 0 : & displacement=disp,compute_anharmonic=.true.,verbose=.false.)
3921 :
3922 0 : write(std_out,*) "Analyti:",gred(mu,ia)
3923 0 : write(std_out,*) "FD :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta)
3924 0 : write(std_out,*) "Diff(%):",abs(100*(gred(mu,ia)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))&
3925 0 : & / (12*delta) )) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) ))
3926 :
3927 : end do
3928 : end do
3929 :
3930 :
3931 : ! Fill the identity matrix
3932 0 : identity = zero
3933 0 : forall(ii=1:3)identity(ii,ii)=1
3934 :
3935 0 : npt=5
3936 0 : delta = 0.0005
3937 0 : deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/)
3938 :
3939 0 : do jj=1,6
3940 0 : write(std_out,*) "strain ",jj
3941 0 : do ii=1,npt
3942 : strain = zero
3943 0 : delt = deltalist(ii)
3944 0 : mat_def = zero
3945 0 : strain(jj) = strain(jj) + delt
3946 :
3947 0 : mat_def(alpha(jj),beta(jj)) = mat_def(alpha(jj),beta(jj)) + half * strain(jj)
3948 0 : mat_def(beta(jj),alpha(jj)) = mat_def(beta(jj),alpha(jj)) + half * strain(jj)
3949 :
3950 0 : mat_def = mat_def + identity
3951 0 : rprimd_def = matmul(mat_def,rprimd)
3952 :
3953 : ! The two options should give the same result
3954 : ! Option 1 => compute the disps and provide them to evaluate
3955 : ! call effective_potential_getDisp(disp,du_delta,natom,rprimd_def,&
3956 : ! & rprimd_ref,1,xred_hist=xred,&
3957 : ! & xcart_ref=eff_pot%supercell%xcart,&
3958 : ! & compute_displacement = .true.,compute_duDelta = .true.)
3959 :
3960 : ! call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd_def,&
3961 : ! & xred=xred,du_delta=du_delta,&
3962 : ! & displacement=disp,strain=strain,&
3963 : ! & compute_anharmonic=.true.,verbose=.false.)
3964 :
3965 : ! Option 2 => compute the disps within evaluate
3966 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd_def,&
3967 0 : & xred=xred,compute_anharmonic=.true.,verbose=.false.)
3968 :
3969 :
3970 :
3971 0 : diff(ii) = energy
3972 :
3973 : end do
3974 :
3975 : ! The two options should give the same result
3976 : ! Option 1 => compute the disps and provide them to evaluate
3977 : ! call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3978 : ! & rprimd_ref,1,xred_hist=xred,&
3979 : ! & xcart_ref=eff_pot%supercell%xcart,&
3980 : ! & compute_displacement = .true.,compute_duDelta = .true.)
3981 :
3982 : ! call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3983 : ! & xred=xred,du_delta=du_delta,&
3984 : ! & displacement=disp,&
3985 : ! & compute_anharmonic=.true.,verbose=.false.)
3986 :
3987 : ! Option 2 => compute the disps within evaluate
3988 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3989 0 : & xred=xred,compute_anharmonic=.true.,verbose=.false.)
3990 :
3991 0 : write(std_out,*) "Analyti:",strten(jj)
3992 0 : write(std_out,*) "FD :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol
3993 0 : write(std_out,*) "Diff(%):",abs(100*(strten(jj)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))&
3994 0 : & / (12*delta) / ucvol)) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol))
3995 :
3996 : end do
3997 :
3998 0 : ABI_FREE(xred)
3999 0 : ABI_FREE(fcart)
4000 0 : ABI_FREE(gred)
4001 0 : ABI_FREE(disp)
4002 0 : ABI_FREE(disp_red)
4003 0 : ABI_FREE(du_delta)
4004 :
4005 0 : end subroutine effective_potential_checkDEV
4006 : !!***
4007 :
4008 : !!****f* m_effective_potential/effective_potential_writeNETCDF
4009 : !! NAME
4010 : !! effective_potential_writeNETCDF
4011 : !!
4012 : !! FUNCTION
4013 : !! This routine print the effective potential into netcdf format
4014 : !! Several options are available
4015 : !!
4016 : !! INPUTS
4017 : !! filename = the name of output file
4018 : !! eff_pot = datatype contains the effective potential
4019 : !! option = option for the format of the xml file
4020 : !! 1 print the xml for a system
4021 : !!
4022 : !! OUTPUT
4023 : !!
4024 : !! SOURCE
4025 :
4026 0 : subroutine effective_potential_writeNETCDF(eff_pot,option,filename)
4027 :
4028 : !Arguments ------------------------------------
4029 : !scalars
4030 : integer, intent(in) :: option
4031 : character(len=fnlen),optional,intent(in) :: filename
4032 : !arrays
4033 : type(effective_potential_type), intent(in) :: eff_pot
4034 :
4035 : !Local variables-------------------------------
4036 : !scalar
4037 : integer :: amu_id,bec_id,ifccell_id,epsinf_id,elastic_id
4038 : integer :: ifc_id,ifcs_id,natom_id,ntypat_id,nrpt_id,npsp_id,typat_id
4039 : integer :: six_id,two_id,xyz_id,znucl_id
4040 : integer :: ncerr,ncid,npsp
4041 : integer :: dimCids(2),dimEids(2),dimIids(6),dimPids(1),dimRids(2),dimXids(2)
4042 : integer :: etotal_id,rprimd_id,xcart_id
4043 : character(len=500) :: msg
4044 : character(len=fnlen) :: namefile
4045 : !arrays
4046 : real(dp) :: strain(9,6)
4047 : ! *************************************************************************
4048 :
4049 : strain(:,1) = (/1,0,0,0,0,0,0,0,0/)
4050 : strain(:,2) = (/0,0,0,0,1,0,0,0,0/)
4051 : strain(:,3) = (/0,0,0,0,0,0,0,0,1/)
4052 : strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/)
4053 : strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/)
4054 : strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/)
4055 :
4056 : !Print only the reference system in xml format
4057 0 : if (option == 1) then
4058 :
4059 0 : if(present(filename)) then
4060 0 : namefile=filename
4061 : else
4062 0 : namefile='ref.nc'
4063 : end if
4064 :
4065 0 : call isfile(namefile,'new')
4066 :
4067 0 : write(msg,'(a,a,a)')ch10,&
4068 0 : & ' Generation of the xml file for the reference structure in ',trim(namefile)
4069 :
4070 0 : call wrtout(ab_out,msg,'COLL')
4071 0 : call wrtout(std_out,msg,'COLL')
4072 :
4073 : ! 1. Create netCDF file
4074 0 : ncerr = nf90_create(path=trim(namefile),cmode=NF90_CLOBBER, ncid=ncid)
4075 0 : NCF_CHECK_MSG(ncerr,"create netcdf history file")
4076 :
4077 : ! 2. Define dimensions
4078 0 : ncerr = nf90_def_dim(ncid,"natom",eff_pot%crystal%natom,natom_id)
4079 0 : NCF_CHECK_MSG(ncerr," define dimension natom")
4080 :
4081 0 : ncerr = nf90_def_dim(ncid,"ntypat",eff_pot%crystal%ntypat,ntypat_id)
4082 0 : NCF_CHECK_MSG(ncerr," define dimension ntypat")
4083 :
4084 0 : ncerr = nf90_def_dim(ncid,"nrpt",eff_pot%harmonics_terms%ifcs%nrpt,nrpt_id)
4085 0 : NCF_CHECK_MSG(ncerr," define dimension ntypat")
4086 :
4087 0 : ncerr = nf90_def_var(ncid, "typat", NF90_INT, natom_id, typat_id)
4088 0 : NCF_CHECK_MSG(ncerr," define variable typat")
4089 :
4090 0 : npsp = size(eff_pot%crystal%znucl)
4091 0 : if (npsp /= eff_pot%crystal%ntypat) then
4092 0 : ABI_WARNING("HIST file does not support alchemical mixing!")
4093 : end if
4094 0 : ncerr = nf90_def_dim(ncid,"npsp",npsp,npsp_id)
4095 0 : NCF_CHECK_MSG(ncerr," define dimension npsp")
4096 :
4097 0 : ncerr = nf90_def_var(ncid, "znucl", NF90_DOUBLE, npsp_id, znucl_id)
4098 0 : NCF_CHECK_MSG(ncerr," define variable znucl")
4099 :
4100 0 : ncerr = nf90_def_dim(ncid,"xyz",3,xyz_id)
4101 0 : NCF_CHECK_MSG(ncerr," define dimension xyz")
4102 :
4103 0 : ncerr = nf90_def_dim(ncid,"six",6,six_id)
4104 0 : NCF_CHECK_MSG(ncerr," define dimension six")
4105 :
4106 0 : ncerr = nf90_def_dim(ncid,"two",2,two_id)
4107 0 : NCF_CHECK_MSG(ncerr," define dimension two")
4108 :
4109 : ! Dimensions for xcart,xred,fcart,gred and vel
4110 0 : dimXids = (/ xyz_id, natom_id /)
4111 : ! Dimensions for rprimd
4112 0 : dimRids = (/ xyz_id, xyz_id /)
4113 : ! Dimensions for ifc
4114 0 : dimIids = (/ 2, xyz_id, natom_id, xyz_id, natom_id, nrpt_id /)
4115 : ! Dimensions for position
4116 0 : dimPids = (/ xyz_id /)
4117 : ! Dimension for elastic constant
4118 0 : dimEids = (/six_id,six_id/)
4119 : ! Dimension for cell
4120 0 : dimCids = (/nrpt_id,3/)
4121 :
4122 : ! 3. Define variables and their attributes (units and mnemonics)
4123 : call ab_define_var(ncid, (/1/), etotal_id, NF90_DOUBLE,&
4124 0 : & "energy","Energy of the reference structure","Ha" )
4125 :
4126 : call ab_define_var(ncid, dimRids, rprimd_id, NF90_DOUBLE,&
4127 0 : & "rprimd","Real space PRIMitive translations, Dimensional","bohr" )
4128 :
4129 : call ab_define_var(ncid, dimRids, epsinf_id, NF90_DOUBLE,&
4130 0 : & "epsilon_inf","Dielectric tensor, Dimensional","epsilon_inf" )
4131 :
4132 : call ab_define_var(ncid, dimEids, elastic_id, NF90_DOUBLE,&
4133 0 : & "elastic","Elastic Constants, Dimensional","Ha" )
4134 :
4135 : call ab_define_var(ncid, dimRids, bec_id, NF90_DOUBLE,&
4136 0 : & "bec","Born Effective Charges, Dimensional","abs(e)" )
4137 :
4138 : call ab_define_var(ncid, dimXids, xcart_id, NF90_DOUBLE,&
4139 0 : & "xcart","vectors (X) of atom positions in CARTesian coordinates","bohr" )
4140 :
4141 : call ab_define_var(ncid, dimIids, ifcs_id, NF90_DOUBLE,&
4142 0 : & "IFCs","Interatomic Forces Constantes in real spaces (short range), Dimensional","Hatree/bohr**2" )
4143 :
4144 : call ab_define_var(ncid, dimIids, ifc_id, NF90_DOUBLE,&
4145 0 : & "IFC","Interatomic Forces Constantes in real spaces (total range), Dimensional","Hatree/bohr**2" )
4146 :
4147 : call ab_define_var(ncid, dimCids, ifccell_id, NF90_DOUBLE,&
4148 0 : & "cell","cell for the ifc, Dimensional","Dimensionless" )
4149 :
4150 : call ab_define_var(ncid, [ntypat_id], amu_id, NF90_DOUBLE,&
4151 0 : & "amu","Masses of each type of atom in atomic mass units", "" )
4152 :
4153 : ! 4. End define mode
4154 0 : ncerr = nf90_enddef(ncid)
4155 0 : NCF_CHECK_MSG(ncerr," end define mode")
4156 :
4157 : ! 5. Write variables
4158 0 : ncerr = nf90_put_var(ncid,etotal_id, eff_pot%energy)
4159 0 : NCF_CHECK_MSG(ncerr," write variable energy")
4160 :
4161 0 : ncerr = nf90_put_var(ncid,rprimd_id, eff_pot%crystal%rprimd)
4162 0 : NCF_CHECK_MSG(ncerr," write variable rprimd")
4163 :
4164 0 : ncerr = nf90_put_var(ncid,epsinf_id, eff_pot%harmonics_terms%epsilon_inf)
4165 0 : NCF_CHECK_MSG(ncerr," write variable epsilon_inf")
4166 :
4167 0 : ncerr = nf90_put_var(ncid,elastic_id , eff_pot%harmonics_terms%elastic_constants)
4168 0 : NCF_CHECK_MSG(ncerr," write variable elastic_constant")
4169 :
4170 0 : ncerr = nf90_put_var(ncid, bec_id, eff_pot%harmonics_terms%zeff)
4171 0 : NCF_CHECK_MSG(ncerr," write variable bec")
4172 :
4173 0 : ncerr = nf90_put_var(ncid,xcart_id, eff_pot%crystal%xcart)
4174 0 : NCF_CHECK_MSG(ncerr," write variable xcart")
4175 :
4176 0 : ncerr = nf90_put_var(ncid,ifccell_id, eff_pot%harmonics_terms%ifcs%cell)
4177 0 : NCF_CHECK_MSG(ncerr," write variable cell")
4178 :
4179 0 : ncerr = nf90_put_var(ncid,ifcs_id, eff_pot%harmonics_terms%ifcs%short_atmfrc)
4180 0 : NCF_CHECK_MSG(ncerr," write variable short ifc")
4181 :
4182 0 : ncerr = nf90_put_var(ncid,ifc_id, eff_pot%harmonics_terms%ifcs%atmfrc)
4183 0 : NCF_CHECK_MSG(ncerr," write variable total ifc")
4184 :
4185 :
4186 : ! 6. Close NetCDF file
4187 0 : ncerr = nf90_close(ncid)
4188 0 : NCF_CHECK_MSG(ncerr," close netcdf history file")
4189 : end if
4190 :
4191 0 : end subroutine effective_potential_writeNETCDF
4192 : !!***
4193 :
4194 : !AM_END_EXPERIMENTAL SECTION
4195 :
4196 1658 : end module m_effective_potential
4197 : !!***
|
