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 27038 : real(dp) :: disp_tmp(3,natom)
2297 27038 : real(dp) :: du_delta_tmp(6,3,natom)
2298 27038 : real(dp) :: fcart_part(3,natom)
2299 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2300 : real(dp) :: strain_tmp(6),strten_part(6)
2301 27038 : real(dp) :: energy_coeff_part(eff_pot%anharmonics_terms%ncoeff)
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 36076719 : disp_tmp(:,:) = zero
2411 198361119 : du_delta_tmp(:,:,:) = zero
2412 13519 : if((.not.present(displacement).or..not.present(du_delta)).and.present(xred))then
2413 : ! Compute the displacement
2414 12558 : call xred2xcart(natom, rprimd, xcart, xred)
2415 : call effective_potential_getDisp(disp_tmp,du_delta_tmp,natom,rprimd,&
2416 : & eff_pot%supercell%rprimd,comm,xcart_hist=xcart,&
2417 : & xcart_ref=eff_pot%supercell%xcart,&
2418 : & compute_displacement = .not.present(displacement),&
2419 12558 : & compute_duDelta = .not.present(du_delta))
2420 : end if
2421 :
2422 : ! or set ftrom the arguments
2423 167279 : if(present(displacement)) disp_tmp(:,:) = displacement(:,:)
2424 859199 : if(present(du_delta))du_delta_tmp(:,:,:) = du_delta(:,:,:)
2425 : !Set to zero the outputs
2426 13519 : energy = zero
2427 36076719 : fcart(:,:) = zero
2428 13519 : strten(:) = zero
2429 :
2430 13519 : if(need_verbose)then
2431 166536 : write(msg, '(80a,2a)' ) ('-',mu=1,80),&
2432 168592 : & ch10,' Components of total energy (in Hartree) :'
2433 2056 : call wrtout(ab_out,msg,'COLL')
2434 2056 : call wrtout(std_out,msg,'COLL')
2435 : end if
2436 :
2437 : !------------------------------------
2438 : ! 2 - Transfert the reference values
2439 : !------------------------------------
2440 :
2441 : ! Set the value of the energy
2442 13519 : energy = eff_pot%energy * ncell
2443 :
2444 13519 : if(need_verbose)then
2445 2056 : write(msg, '(a,a,1ES24.16,a)' ) ch10,' Energy of the reference structure :',&
2446 4112 : & energy,' Hartree'
2447 2056 : call wrtout(ab_out,msg,'COLL')
2448 2056 : call wrtout(std_out,msg,'COLL')
2449 : end if
2450 :
2451 : !Set the value of the initial strees (should be zero if the system is relaxed)
2452 : !According to the original definition of the model, the reference should always be
2453 : ! a critical point of the PES so the first order derivative should not be take into account...
2454 : ! do ii=1,6
2455 : ! energy = energy + ncell * eff_pot%strten(ii)*strain_tmp(ii)
2456 : ! end do
2457 : ! strten(:) = ncell * eff_pot%strten(:)
2458 : !Set the value of the initial forces (should be zero if the system is relaxed)
2459 : ! ii = 1
2460 : ! This next part is not working, need to apply strain to the original fcart
2461 : ! do ia=1,eff_pot%supercell%natom
2462 : ! fcart(:,ia) = eff_pot%fcart(:,ii)
2463 : ! do mu=1,3
2464 : ! energy = energy + eff_pot%fcart(mu,ii)*disp_tmp(mu,ii)
2465 : ! end do
2466 : ! ii = ii + 1
2467 : ! if(ii > eff_pot%crystal%natom) ii = 1
2468 : ! end do
2469 :
2470 : !------------------------------------
2471 : ! 3 - Computation of the IFC part :
2472 : !------------------------------------
2473 :
2474 : energy_part = zero
2475 36076719 : fcart_part(:,:)= zero
2476 :
2477 : call harmonics_terms_evaluateIFC(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,:),disp_tmp,&
2478 : & energy_part,fcart_part,eff_pot%supercell%natom,&
2479 : & eff_pot%crystal%natom,eff_pot%mpi_ifc%my_ncell,&
2480 : & eff_pot%mpi_ifc%my_nrpt,eff_pot%mpi_ifc%my_atmrpt_index,&
2481 : & eff_pot%mpi_ifc%my_index_cells,sc_size,eff_pot%mpi_ifc%my_rpt,&
2482 13519 : & eff_pot%mpi_ifc%comm)
2483 :
2484 13519 : if(need_verbose)then
2485 2056 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the ifc part :',&
2486 4112 : & energy_part,' Hartree'
2487 2056 : call wrtout(ab_out,msg,'COLL')
2488 2056 : call wrtout(std_out,msg,'COLL')
2489 :
2490 2056 : if(abs(energy_part) < tol16 .and. eff_pot%anharmonics_terms%ncoeff == 0)then
2491 25 : write(msg, '(8a)' )ch10,&
2492 25 : & ' --- !WARNING!',ch10,&
2493 25 : & ' The harmonic part is negative, the simulation will diverge',ch10,&
2494 25 : & ' if the anharmonic part is not used',ch10,&
2495 50 : & ' ---'
2496 25 : call wrtout(std_out,msg,"COLL")
2497 : end if
2498 : end if
2499 :
2500 13519 : energy = energy + energy_part
2501 36076719 : fcart(:,:)= fcart(:,:) + fcart_part(:,:)
2502 :
2503 :
2504 : !----------------------------------------------------
2505 : ! 4 - Computation of the elastic part of the energy :
2506 : !----------------------------------------------------
2507 :
2508 13519 : energy_part = zero
2509 36076719 : fcart_part(:,:)= zero
2510 13519 : strten_part(:) = zero
2511 :
2512 : ! write(*,*) "--- STRTEN before harmonics evaluate --- "
2513 : ! write(*,*) strten(:)
2514 : call harmonics_terms_evaluateElastic(eff_pot%harmonics_terms%elastic_constants,disp_tmp,energy_part,&
2515 : & fcart_part,eff_pot%supercell%natom,eff_pot%crystal%natom,&
2516 13519 : & ncell,eff_pot%harmonics_terms%strain_coupling,strten_part,strain_tmp)
2517 :
2518 :
2519 13519 : if(has_strain.and.need_verbose)then
2520 975 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the elastic part :',&
2521 1950 : & energy_part,' Hartree'
2522 975 : call wrtout(ab_out,msg,'COLL')
2523 975 : call wrtout(std_out,msg,'COLL')
2524 : end if
2525 :
2526 13519 : energy = energy + energy_part
2527 36076719 : fcart(:,:) = fcart(:,:) + fcart_part(:,:)
2528 94633 : strten(:) = strten(:) + strten_part(:)
2529 : ! write(*,*) "--- STRTEN after harmonics evaluate --- "
2530 : ! write(*,*) strten(:)
2531 :
2532 : !------------------------------------
2533 : ! 5 - Treat 3rd order strain-coupling:
2534 : !------------------------------------
2535 13519 : if (need_anharmonic.and.eff_pot%has_anharmonicsTerms) then
2536 1104 : energy_part = zero
2537 1104 : strten_part(:) = zero
2538 648144 : fcart_part(:,:)= zero
2539 :
2540 : ! 1-Part due to the anharmonic strain
2541 : if(eff_pot%anharmonics_terms%has_elastic3rd.or.&
2542 1104 : & eff_pot%anharmonics_terms%has_elastic4th.or.&
2543 : & eff_pot%anharmonics_terms%has_elastic_displ)then
2544 : call anharmonics_terms_evaluateElastic(disp_tmp,energy_part,fcart_part,&
2545 : & eff_pot%supercell%natom,&
2546 : & eff_pot%crystal%natom,ncell,strten_part,strain_tmp,&
2547 : & elastic3rd=eff_pot%anharmonics_terms%elastic3rd,&
2548 : & elastic4th=eff_pot%anharmonics_terms%elastic4th,&
2549 0 : & elastic_displacement=eff_pot%anharmonics_terms%elastic_displacement)
2550 :
2551 0 : if(need_verbose)then
2552 0 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the anharmonic elastic part :',&
2553 0 : & energy_part,' Hartree'
2554 0 : call wrtout(ab_out,msg,'COLL')
2555 0 : call wrtout(std_out,msg,'COLL')
2556 : end if
2557 0 : energy = energy + energy_part
2558 1104 : strten(:) = strten(:) + strten_part(:)
2559 : end if
2560 :
2561 : ! 2-Part due to the strain-phonon coupling
2562 1104 : if (eff_pot%anharmonics_terms%has_strain_coupling) then
2563 : energy_part = zero
2564 0 : strten_part(:) = zero
2565 0 : fcart_part(:,:)= zero
2566 :
2567 : call anharmonics_terms_evaluateIFCStrainCoupling(eff_pot%anharmonics_terms%phonon_strain,&
2568 : & disp_tmp,energy_part,fcart_part,&
2569 : & eff_pot%supercell%natom,&
2570 : & eff_pot%crystal%natom,&
2571 : & sc_size,strain_tmp,&
2572 : & strten_part,eff_pot%mpi_ifc%my_cells,&
2573 : & eff_pot%mpi_ifc%my_ncell,&
2574 : & eff_pot%mpi_ifc%my_index_cells,&
2575 0 : & eff_pot%mpi_ifc%comm)
2576 :
2577 0 : if(need_verbose)then
2578 0 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the 3rd (strain-phonon coupling):',&
2579 0 : & energy_part,' Hartree'
2580 0 : call wrtout(ab_out,msg,'COLL')
2581 0 : call wrtout(std_out,msg,'COLL')
2582 : end if
2583 0 : energy = energy + energy_part
2584 0 : fcart = fcart + fcart_part
2585 0 : strten = strten + strten_part
2586 : end if
2587 :
2588 : end if
2589 :
2590 : !----------------------------------
2591 : ! 6 - Treat polynomial coefficient:
2592 : !----------------------------------
2593 12172 : if(need_anharmonic.and.eff_pot%anharmonics_terms%ncoeff > zero)then
2594 10265 : energy_part = zero
2595 35095225 : fcart_part(:,:) = zero
2596 10265 : strten_part(:) = zero
2597 : call polynomial_coeff_evaluate(eff_pot%anharmonics_terms%coefficients,disp_tmp,&
2598 : & energy_part,energy_coeff_part,fcart_part,eff_pot%supercell%natom,&
2599 : & eff_pot%crystal%natom,eff_pot%anharmonics_terms%ncoeff,&
2600 : & sc_size,strain_tmp,strten_part,eff_pot%mpi_coeff%my_ncell,&
2601 : & eff_pot%mpi_coeff%my_index_cells,eff_pot%mpi_coeff%comm,&
2602 10671 : & filename=filename)
2603 :
2604 10265 : if(need_verbose)then
2605 706 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the fitted coefficient :',&
2606 1412 : & energy_part,' Hartree'
2607 706 : call wrtout(ab_out,msg,'COLL')
2608 706 : call wrtout(std_out,msg,'COLL')
2609 : end if
2610 10265 : energy = energy + energy_part
2611 35095225 : fcart(:,:) = fcart(:,:) + fcart_part(:,:)
2612 71855 : strten(:) = strten(:) + strten_part(:)
2613 : end if
2614 :
2615 : !---------------------------------
2616 : ! 7 - Compute confinement
2617 : !---------------------------------
2618 13519 : if(eff_pot%confinement%need_confinement) then
2619 :
2620 : energy_part = zero
2621 :
2622 : call polynomial_conf_evaluate(disp_tmp,eff_pot%confinement%cutoff_disp,energy_part,&
2623 : & eff_pot%confinement%factor_disp,&
2624 : & eff_pot%confinement%factor_strain,fcart_part,strain_tmp,&
2625 : & eff_pot%confinement%cutoff_strain,strten_part,&
2626 : & eff_pot%confinement%power_disp,eff_pot%confinement%power_strain,&
2627 : & eff_pot%mpi_coeff%my_cells,&
2628 : & eff_pot%supercell%natom,eff_pot%crystal%natom,&
2629 : & eff_pot%mpi_coeff%my_ncell,eff_pot%mpi_coeff%my_index_cells,&
2630 0 : & eff_pot%mpi_coeff%comm)
2631 :
2632 0 : energy = energy + energy_part
2633 :
2634 0 : if(abs(energy_part) > tol10 .and. need_verbose )then
2635 0 : write(msg, '(a,1ES24.16,a)' ) ' Energy of the confinement part :',&
2636 0 : & energy_part,' Hartree'
2637 0 : call wrtout(ab_out,msg,'COLL')
2638 0 : call wrtout(std_out,msg,'COLL')
2639 : end if
2640 : end if
2641 :
2642 : !-------------------------------------------
2643 : !!! External_ELECTRIC_FILED
2644 : ! 7.5 - Compute Forces and energies from external electric field
2645 : !------------------------------------------
2646 13519 : if(present(efield_type))then
2647 9864 : if (efield_type /= 0) then
2648 : call calculate_forces_efield(eff_pot,energy_part,eff_pot%mpi_coeff%comm,fcart_part,natom,&
2649 : & disp_tmp,xcart,efield_type,efield,efield_lambda,efield_period,&
2650 60 : & efield_phase,efield_gmean,efield_gvel,efield_sigma,efield_background,time,nefield)
2651 60 : if(need_verbose)then
2652 60 : write(msg, '(a,1ES24.16,a)' ) ' Energy from electric field is :',&
2653 120 : & energy_part,' Hartree'
2654 60 : call wrtout(ab_out,msg,'COLL')
2655 60 : call wrtout(std_out,msg,'COLL')
2656 :
2657 : end if
2658 60 : energy = energy + energy_part
2659 19260 : fcart = fcart + fcart_part
2660 : endif
2661 : end if
2662 : !-------------------------------------------
2663 : ! 8 - Compute electronic Part with SCALE-UP
2664 : !------------------------------------------
2665 :
2666 13519 : err_eng = .FALSE.
2667 13519 : err_for = .FALSE.
2668 :
2669 13519 : energy_part = 0
2670 36076719 : fcart_part = 0
2671 :
2672 13519 : if(need_elec_eval)then
2673 : #if defined DEV_MS_SCALEUP
2674 : update_dens = .TRUE.
2675 : write(msg,'(a)') ' wohoo i was here and call scale-up now!---STILL WOHOOO!---'
2676 : call wrtout(ab_out,msg,'COLL')
2677 : call wrtout(std_out,msg,'COLL')
2678 : err_eng = global_calculate_energy(energy_part,strain_tmp,disp_tmp,natom_for_scup)
2679 : err_for = global_calculate_forces(fcart_part,strain_tmp,disp_tmp,natom_for_scup,update_dens)
2680 : write(msg,'(a,1ES24.16,a)') 'The Energy of the electronic model is: ', energy_part,' unit?'
2681 : call wrtout(ab_out,msg,'COLL')
2682 : call wrtout(std_out,msg,'COLL')
2683 : #endif
2684 :
2685 0 : energy = energy + energy_part
2686 0 : do ia = 1, natom
2687 0 : fcart(:,ia) = fcart(:,ia) + fcart_part(:,ia)
2688 : end do
2689 : endif
2690 :
2691 :
2692 :
2693 : !-----------------------------------
2694 : ! 9 - Add stress due to forces on
2695 : ! atoms and variation of disp
2696 : ! due to strain
2697 : !-----------------------------------
2698 :
2699 : ! multiply forces by -1
2700 36076719 : fcart = -1 * fcart
2701 :
2702 :
2703 :
2704 :
2705 : !Calculatei
2706 13519 : strten_part(:) = zero
2707 9029319 : do ia = 1,eff_pot%supercell%natom
2708 36076719 : do mu = 1,3
2709 198347600 : do alpha=1,6
2710 189331800 : strten_part(alpha) = strten_part(alpha) - fcart(mu,ia) * du_delta_tmp(alpha,mu,ia)
2711 : end do
2712 : end do
2713 : end do
2714 94633 : strten(:) = strten(:) + strten_part(:)
2715 :
2716 : ! write(*,*) "--- STRTEN after du/deta application --- "
2717 : ! write(*,*) strten(:)
2718 :
2719 : !---------------------------------
2720 : ! 10 - Renormalize stress
2721 : !---------------------------------
2722 :
2723 : ! divide stess tensor by ucvol
2724 94633 : do alpha=1,6
2725 94633 : if (alpha < 4) then
2726 40557 : strten(alpha) = strten(alpha) * (1 + strain_tmp(alpha))/ ucvol
2727 : else
2728 40557 : strten(alpha) = strten(alpha) * (1 - strain_tmp(alpha)**2)/ ucvol
2729 : end if
2730 : end do
2731 :
2732 : ! write(*,*) "--- STRTEN after /ucvol --- "
2733 : ! write(*,*) strten(:)
2734 : ! Redistribute the residuale of the forces
2735 13519 : call effective_potential_distributeResidualForces(eff_pot,fcart,eff_pot%supercell%natom)
2736 :
2737 13519 : call fcart2gred(fcart,gred,rprimd,natom)
2738 :
2739 : !------------------------------------
2740 : ! 11 - Final Print:
2741 : !------------------------------------
2742 :
2743 13519 : if(need_verbose)then
2744 2056 : write(msg, '(2a,es21.14)' ) ch10,&
2745 4112 : & ' >>>>>>>>> Etotal= ',energy
2746 2056 : call wrtout(ab_out,msg,'COLL')
2747 2056 : call wrtout(std_out,msg,'COLL')
2748 :
2749 2056 : write(msg,'(2a,1p,e15.7,a)') ch10,' Unit cell volume ucvol=',ucvol+tol10,' bohr^3'
2750 2056 : call wrtout(std_out, msg,'COLL')
2751 :
2752 168592 : write(msg, '(a,80a,3a)' ) ch10,('-',mu=1,80),ch10,&
2753 170648 : & ' Cartesian components of stress tensor (hartree/bohr^3)'
2754 2056 : call wrtout(ab_out,msg,'COLL')
2755 2056 : call wrtout(std_out, msg,'COLL')
2756 :
2757 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2758 2056 : & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4)
2759 2056 : call wrtout(ab_out,msg,'COLL')
2760 2056 : call wrtout(std_out, msg,'COLL')
2761 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2762 2056 : & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5)
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(3 3)=',strten(3),' sigma(2 1)=',strten(6)
2767 2056 : call wrtout(ab_out,msg,'COLL')
2768 2056 : call wrtout(std_out, msg,'COLL')
2769 :
2770 : ! Also output the pressure (minus one third the trace of the stress
2771 : ! tensor.
2772 2056 : write(msg, '(a,a,es12.4,a)' ) ch10,&
2773 2056 : & '-Cartesian components of stress tensor (GPa) [Pressure=',&
2774 4112 : & -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]'
2775 :
2776 2056 : call wrtout(ab_out,msg,'COLL')
2777 2056 : call wrtout(std_out, msg,'COLL')
2778 :
2779 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2780 2056 : & '- sigma(1 1)=',strten(1)*HaBohr3_GPa,&
2781 4112 : & ' sigma(3 2)=',strten(4)*HaBohr3_GPa
2782 2056 : call wrtout(ab_out,msg,'COLL')
2783 2056 : call wrtout(std_out, msg,'COLL')
2784 : write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) &
2785 2056 : & '- sigma(2 2)=',strten(2)*HaBohr3_GPa,&
2786 4112 : & ' sigma(3 1)=',strten(5)*HaBohr3_GPa
2787 2056 : call wrtout(ab_out,msg,'COLL')
2788 2056 : call wrtout(std_out, msg,'COLL')
2789 : write(msg, '(a,1p,e16.8,a,1p,e16.8,a)' ) &
2790 2056 : & '- sigma(3 3)=',strten(3)*HaBohr3_GPa,&
2791 4112 : & ' sigma(2 1)=',strten(6)*HaBohr3_GPa
2792 2056 : call wrtout(ab_out,msg,'COLL')
2793 2056 : call wrtout(std_out, msg,'COLL')
2794 166536 : write(msg, '(80a,a)' ) ('-',mu=1,80),ch10
2795 2056 : call wrtout(ab_out,msg,'COLL')
2796 2056 : call wrtout(std_out,msg,'COLL')
2797 : end if
2798 :
2799 :
2800 :
2801 13519 : ABI_FREE(xcart)
2802 :
2803 27999 : end subroutine effective_potential_evaluate
2804 :
2805 :
2806 : !****f* m_effective_potential/calculate_forces_efield
2807 : !!
2808 : !! NAME
2809 : !! calculate_forces_efield
2810 : !!
2811 : !! FUNCTION
2812 : !! Compute the action of the electric field
2813 : !! different values of efield_type allow different functional forms
2814 : !! including spatial inhomogeneous fields and time dependent fields
2815 : !! INPUTS
2816 : !!
2817 : !! effpot for the structural information about supercells etc
2818 : !! comm parallel communication
2819 : !! natom number of atoms on the supercell
2820 : !! disp_tmp instantaneous atomic displacements
2821 : !! xcart cartesian coordinates of the atoms
2822 : !! time instantaneous time
2823 : !! efield_type different electric field types
2824 : !! 1- homogeneous 2- AC Fields
2825 : !! 3- cos-like spatial modulation 4-cos-like spatial and time modulated
2826 : !! 5- Gaussian field 6-Double cos-like modulation
2827 : !! efield_lambda spatial modulation on cos-like fields
2828 : !! efield_period Time modulation on AC fields
2829 : !! efield_phase Phase shift for cos-like
2830 : !! efield_gmean Center of the Gaussian field
2831 : !! efield_gvel Velocity of the Gaussian field
2832 : !! efield_sigma Std deviation of the efield
2833 : !! efield_background Homogeneous background
2834 : !!
2835 : !! OUTPUT
2836 : !!
2837 : !! energy_part to update the energy contribution of the electric field
2838 : !! fcart_part to update the contribution to the forces
2839 : !!
2840 : !!
2841 60 : subroutine calculate_forces_efield(eff_pot,energy_part,comm,fcart_part,natom,disp_tmp,xcart,&
2842 : & efield_type,efield,efield_lambda,efield_period,efield_phase,&
2843 : & efield_gmean,efield_gvel,efield_sigma,efield_background,time,nefield)
2844 : !Arguments ------------------------------------
2845 : integer, intent(in) :: natom,comm
2846 : type(effective_potential_type),intent(in) :: eff_pot
2847 : real(dp),intent(inout) :: energy_part
2848 : real(dp),intent(inout) :: fcart_part(3,natom)
2849 : integer, intent(in) :: efield_type
2850 : integer, intent(in), optional :: nefield
2851 : real(dp),intent(in),optional :: efield(3,*)
2852 : real(dp),intent(in),optional :: efield_background(3)
2853 : real(dp),intent(inout),optional :: efield_lambda(3,*)
2854 : real(dp),intent(in),optional :: efield_gmean(3)
2855 : real(dp),intent(in),optional :: efield_gvel(3)
2856 : real(dp),intent(in),optional :: efield_period
2857 : real(dp),intent(in),optional :: efield_phase(*)
2858 : real(dp),intent(in),optional :: efield_sigma
2859 : real(dp),intent(in),optional :: time
2860 : real(dp),intent(in),optional :: disp_tmp(3,natom),xcart(3,natom)
2861 : !Arguments ------------------------------------
2862 : !Internal variables----------------------------
2863 : integer :: ii,ia
2864 : real(dp):: argumento,argumento2
2865 : integer :: icell,ierr,kk,ll
2866 : !array
2867 : real(dp) :: ext_field(3),ext_field2(3), temp_pol(3),efield_kk(3),pos_gcenter(3),ext_field_bck(3)
2868 : !Internal variables----------------------------
2869 70 : select case (efield_type)
2870 : case(1)
2871 10 : ext_field = zero
2872 10 : temp_pol = zero
2873 10 : energy_part = zero
2874 3210 : fcart_part(:,:) = zero
2875 10 : if (present(efield)) then
2876 40 : ext_field = -1 * efield(:,1) ! Use first field by default
2877 : endif
2878 10 : if (present(efield_background)) then
2879 40 : ext_field = ext_field - efield_background
2880 : endif
2881 :
2882 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2883 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2884 970 : do ia = 1,eff_pot%crystal%natom
2885 800 : kk = ii + ia
2886 : !We compute the polarization from the atomic displacements
2887 2400 : temp_pol = temp_pol + matmul(disp_tmp(:,kk),&
2888 23200 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2889 : !We compute the forces as the Born effective charges times the electric field
2890 20160 : fcart_part(:,kk) = matmul(ext_field, eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2891 : end do
2892 : end do
2893 : !We compute the energy as the Polarization times the electric field
2894 40 : energy_part = DOT_PRODUCT(ext_field(:),temp_pol(:))
2895 10 : call xmpi_sum(energy_part, comm, ierr)
2896 10 : call xmpi_sum(fcart_part, comm, ierr)
2897 : case (2)
2898 10 : if (present(efield) .and. present(nefield)) then
2899 40 : ext_field = -1 * efield(:,1) ! Use first field by default
2900 : endif
2901 10 : if (present(efield_background)) then
2902 40 : ext_field_bck = -1 * efield_background
2903 : endif
2904 10 : temp_pol = zero
2905 10 : energy_part = zero
2906 3210 : fcart_part(:,:) = zero
2907 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2908 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2909 970 : do ia = 1, eff_pot%crystal%natom
2910 800 : kk = ii + ia
2911 : !We compute the electric field value at the given time
2912 3200 : efield_kk(:)=cos(-two_pi/efield_period*time+efield_phase(1))*ext_field(:)+ext_field_bck(:)
2913 800 : temp_pol = matmul(disp_tmp(:,kk),&
2914 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2915 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2916 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
2917 : end do
2918 : end do
2919 10 : call xmpi_sum(energy_part, comm, ierr)
2920 10 : call xmpi_sum(fcart_part , comm, ierr)
2921 : case (3)
2922 10 : if (present(efield)) then
2923 40 : ext_field =-1* efield(:, 1)
2924 : endif
2925 10 : if (present(efield_background)) then
2926 40 : ext_field_bck =-1* efield_background
2927 : endif
2928 10 : temp_pol = zero
2929 10 : energy_part = zero
2930 3210 : fcart_part(:,:) = zero
2931 40 : do ll=1,3
2932 : !We update the value of the spatial modulation to convey with PBC in case
2933 : !the lattice has been updated
2934 40 : if(efield_lambda(ll,1)/=0.0 .and. abs(efield_lambda(ll,1))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
2935 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))
2936 : endif
2937 : enddo
2938 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2939 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2940 970 : do ia = 1, eff_pot%crystal%natom
2941 800 : kk = ii + ia
2942 800 : argumento=0.0
2943 3200 : do ll=1,3
2944 3200 : if (efield_lambda(ll,1)/=0.0) then
2945 800 : argumento=argumento+two_pi*xcart(ll,kk)/efield_lambda(ll,1)
2946 : endif
2947 : enddo
2948 3200 : efield_kk(:)=cos(argumento+efield_phase(1))*ext_field(:)+ext_field_bck(:)
2949 800 : temp_pol = matmul(disp_tmp(:,kk),&
2950 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2951 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2952 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
2953 : end do
2954 : end do
2955 10 : call xmpi_sum(energy_part, comm, ierr)
2956 10 : call xmpi_sum(fcart_part , comm, ierr)
2957 : case (4)
2958 10 : if (present(efield)) then
2959 40 : ext_field =-1*efield(:, 1)
2960 : endif
2961 10 : if (present(efield_background)) then
2962 40 : ext_field_bck =-1* efield_background
2963 : endif
2964 10 : temp_pol = zero
2965 10 : energy_part = zero
2966 3210 : fcart_part(:,:) = zero
2967 40 : do ll=1,3
2968 40 : if(efield_lambda(ll,1)/=0.0 .and. abs(efield_lambda(ll,1))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
2969 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))
2970 : endif
2971 : enddo
2972 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
2973 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
2974 970 : do ia = 1, eff_pot%crystal%natom
2975 800 : kk = ii + ia
2976 800 : argumento=0.0
2977 3200 : do ll=1,3
2978 3200 : if (efield_lambda(ll,1)/=0.0) then
2979 800 : argumento=argumento+two_pi*xcart(ll,kk)/efield_lambda(ll,1)
2980 : endif
2981 : enddo
2982 3200 : efield_kk(:)=cos(argumento-two_pi/efield_period*time+efield_phase(1))*ext_field(:)+ext_field_bck(:)
2983 800 : temp_pol = matmul(disp_tmp(:,kk),&
2984 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2985 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
2986 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
2987 : end do
2988 : end do
2989 10 : call xmpi_sum(energy_part, comm, ierr)
2990 10 : call xmpi_sum(fcart_part , comm, ierr)
2991 : case (5)
2992 10 : if (present(efield) .and. present(nefield)) then
2993 40 : ext_field = -1 * efield(:,1) ! Use first field by default
2994 : endif
2995 10 : if (present(efield_background)) then
2996 40 : ext_field_bck =-1* efield_background
2997 : endif
2998 10 : temp_pol = zero
2999 10 : energy_part = zero
3000 3210 : fcart_part(:,:) = zero
3001 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
3002 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
3003 970 : do ia = 1, eff_pot%crystal%natom
3004 800 : kk = ii + ia
3005 3200 : pos_gcenter(:)=xcart(:,kk)-efield_gmean(:)-efield_gvel(:)*time
3006 3200 : do ll=1,3
3007 : !PBC for the Gaussian field
3008 : pos_gcenter(ll)=pos_gcenter(ll)-eff_pot%supercell%rprimd(ll,ll)*&
3009 3200 : &nint(pos_gcenter(ll)/eff_pot%supercell%rprimd(ll,ll))
3010 : enddo
3011 5600 : efield_kk(:)=EXP(-DOT_PRODUCT(pos_gcenter(:),pos_gcenter(:))/efield_sigma/efield_sigma)*ext_field(:)+ext_field_bck(:)
3012 800 : temp_pol = matmul(disp_tmp(:,kk),&
3013 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3014 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3015 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
3016 : end do
3017 : end do
3018 10 : call xmpi_sum(energy_part, comm, ierr)
3019 10 : call xmpi_sum(fcart_part , comm, ierr)
3020 : case (6)
3021 10 : if (present(efield) .and. present(nefield)) then
3022 10 : if (nefield >= 1) then
3023 40 : ext_field = -1 * efield(:,1)
3024 40 : ext_field2 = -1 * efield(:,2)
3025 : endif
3026 : endif
3027 10 : if (present(efield_background)) then
3028 40 : ext_field_bck = -1 * efield_background
3029 : endif
3030 10 : temp_pol = zero
3031 10 : energy_part = zero
3032 3210 : fcart_part(:,:) = zero
3033 40 : do ll=1,3
3034 40 : if(present(efield_lambda)) then
3035 30 : if(efield_lambda(ll,1)/=0.0 .and. abs(efield_lambda(ll,1))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
3036 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))
3037 : endif
3038 30 : if(present(nefield) .and. nefield >= 2) then
3039 30 : if(efield_lambda(ll,2)/=0.0 .and. abs(efield_lambda(ll,2))<2.0*eff_pot%supercell%rprimd(ll,ll)) then
3040 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))
3041 : endif
3042 : endif
3043 : endif
3044 : enddo
3045 170 : do icell = 1,eff_pot%mpi_coeff%my_ncell
3046 160 : ii = eff_pot%mpi_coeff%my_index_cells(4,icell)
3047 970 : do ia = 1, eff_pot%crystal%natom
3048 800 : kk = ii + ia
3049 800 : argumento=0.0
3050 800 : argumento2=0.0
3051 3200 : do ll=1,3
3052 2400 : if (present(efield_lambda) .and. efield_lambda(ll,1)/=0.0) then
3053 800 : argumento=argumento+two_pi*xcart(ll,kk)/efield_lambda(ll,1)
3054 : endif
3055 : if (present(efield_lambda) .and. present(nefield) &
3056 3200 : & .and. nefield >= 2 .and. efield_lambda(ll,2)/=0.0) then
3057 800 : argumento2=argumento2+two_pi*xcart(ll,kk)/efield_lambda(ll,2)
3058 : endif
3059 : enddo
3060 3200 : efield_kk(:)=cos(argumento+efield_phase(1))*ext_field(:)
3061 800 : if (present(nefield) .and. nefield >= 2) then
3062 3200 : efield_kk(:)=efield_kk(:)+cos(argumento2+efield_phase(2))*ext_field2
3063 : endif
3064 800 : if (present(efield_background)) then
3065 3200 : efield_kk(:)=efield_kk(:)+ext_field_bck(:)
3066 : endif
3067 800 : temp_pol = matmul(disp_tmp(:,kk),&
3068 17600 : & eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3069 20000 : fcart_part(:,kk) = matmul(efield_kk , eff_pot%harmonics_terms%zeff(:,:,eff_pot%supercell%atom_indexing(kk)))
3070 3360 : energy_part=energy_part+DOT_PRODUCT(efield_kk(:),temp_pol(:))
3071 : end do
3072 : end do
3073 10 : call xmpi_sum(energy_part, comm, ierr)
3074 120 : call xmpi_sum(fcart_part , comm, ierr)
3075 : end select
3076 : !!***
3077 : !!***
3078 : !!Subroutine Inh Elec Field Ends-----
3079 : !!***
3080 : !!***
3081 60 : end subroutine calculate_forces_efield
3082 :
3083 :
3084 : !!***
3085 :
3086 : !****f* m_effective_potential/effective_potential_getDisp
3087 : !!
3088 : !! NAME
3089 : !! effective_potential_getDisp
3090 : !!
3091 : !! FUNCTION
3092 : !! Compute cartesian coordinates of the displacment
3093 : !! between two configurations.
3094 : !! You can provide reduce or cartesian coordinates for both perturbed structure
3095 : !! and reference structure
3096 : !!
3097 : !! INPUTS
3098 : !! natom = number of atoms
3099 : !! rprimd_hist = dimensional primitive translations for real space of the perturbed structure
3100 : !! rprimd_ref = dimensional primitive translations for real space of the reference structure
3101 : !! comm = mpi communicator
3102 : !! xcart_hist(3,natom) = optional, cartesian coordinates of the atoms in the perturbed structure
3103 : !! xred_hist(3,natom) = optional, reduced coordinates of the atoms in the perturbed structure
3104 : !! xred_ref(3,natom) = optional, reduced coordinates of the atoms in the reference structure
3105 : !! xcart_ref(3,natom) = optional, cartesian coordinates of the atoms in the reference structure
3106 : !! compute_displacement= optional, flag to compute the displacement array (default is true)
3107 : !! compute_duDelta = optional, flag to compute the du_delta array (default is true)
3108 : !! Be careful, if you specify compute_displacement=.false.,
3109 : !! you need to provide to this routine the good displacement array!
3110 : !!
3111 : !!
3112 : !! OUTPUT
3113 : !! displacement(3,natom) = cartesian atomic displacement between two configurations (Bohr)
3114 : !! du_delta(6,3,natom) = variation of the displacmeent wr to strain
3115 : !!
3116 : !! SOURCE
3117 :
3118 27038 : subroutine effective_potential_getDisp(displacement,du_delta,natom,rprimd_hist,rprimd_ref,comm,&
3119 26077 : & xcart_hist,xred_hist,xred_ref,xcart_ref,compute_displacement,&
3120 : & compute_duDelta)
3121 :
3122 : !Arguments ------------------------------------
3123 : !scalars
3124 : integer, intent(in) :: natom,comm
3125 : logical,optional,intent(in) :: compute_displacement,compute_duDelta
3126 : !array
3127 : real(dp),intent(in) :: rprimd_ref(3,3),rprimd_hist(3,3)
3128 : real(dp),intent(out) :: displacement(3,natom),du_delta(6,3,natom)
3129 : real(dp),intent(in),optional :: xred_hist(3,natom),xcart_hist(3,natom)
3130 : real(dp),intent(in),optional :: xred_ref(3,natom),xcart_ref(3,natom)
3131 : !Local variables-------------------------------
3132 : !scalar
3133 : integer :: ii,ia,ib,mu
3134 : integer :: ierr,nproc,my_rank,natom_alone,my_natom
3135 : character(len=500) :: msg
3136 : logical :: has_strain
3137 : logical :: need_displacement,need_duDelta
3138 : !array
3139 : integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
3140 13519 : integer,allocatable :: my_atoms(:)
3141 : type(strain_type) :: strain
3142 27038 : real(dp) :: xcart_hist_tmp(3,natom),xcart_ref_tmp(3,natom)
3143 0 : real(dp) :: xred_ref_tmp(3,natom),strain_inv(3,3),strain_inv_u(3),strain_tmp(3,3)
3144 : ! *************************************************************************
3145 :
3146 13519 : if (.not.(present(xred_ref).or.present(xcart_ref))) then
3147 : write(msg, '(3a)' )&
3148 0 : & 'You need at least give xcart_ref or xred_ref '
3149 0 : ABI_ERROR(msg)
3150 : end if
3151 :
3152 13519 : if (.not.(present(xred_hist).or.present(xcart_hist))) then
3153 : write(msg, '(3a)' )&
3154 0 : & 'You need at least give xcart_hist or xred_hist '
3155 0 : ABI_ERROR(msg)
3156 : end if
3157 :
3158 13519 : need_duDelta = .TRUE.
3159 13519 : need_displacement = .TRUE.
3160 13519 : if(present(compute_duDelta)) need_duDelta = compute_duDelta
3161 13519 : if(present(compute_displacement)) need_displacement = compute_displacement
3162 13519 : if(.not.need_duDelta .and. .not.need_displacement) return
3163 :
3164 : !--------------------------------------------
3165 : ! 0 - Set the MPI
3166 : !--------------------------------------------
3167 13519 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3168 13519 : natom_alone = mod(natom,nproc)
3169 13519 : my_natom = int(aint(real(natom,sp)/(nproc)))
3170 13519 : if(my_rank >= (nproc-natom_alone)) then
3171 200 : my_natom = my_natom + 1
3172 : end if
3173 40557 : ABI_MALLOC(my_atoms,(my_natom))
3174 8969039 : my_atoms = 0
3175 8969039 : do ii=1,my_natom
3176 8969039 : if(my_rank >= (nproc-natom_alone))then
3177 : my_atoms(ii)=(int(aint(real(natom,sp)/nproc)))*(my_rank)+&
3178 8500 : & (my_rank - (nproc-natom_alone)) + ii
3179 : else
3180 8947020 : my_atoms(ii)=(my_natom)*(my_rank) + ii
3181 : end if
3182 : end do
3183 :
3184 : !--------------------------------------------
3185 : ! 1 - Get the strain for this step
3186 : !--------------------------------------------
3187 13519 : has_strain = .FALSE.
3188 13519 : call strain_get(strain,rprim=rprimd_ref,rprim_def=rprimd_hist)
3189 13519 : if (strain%name /= "reference") then
3190 18548 : has_strain = .TRUE.
3191 60281 : strain_tmp = strain%strain
3192 4637 : strain_tmp(1,1) = strain_tmp(1,1) + 1.0
3193 4637 : strain_tmp(2,2) = strain_tmp(2,2) + 1.0
3194 4637 : strain_tmp(3,3) = strain_tmp(3,3) + 1.0
3195 : ! get (1+eta)^-1
3196 4637 : call matr3inv(strain_tmp,strain_inv)
3197 : else if (strain%name == "reference") then
3198 8882 : strain_inv(:,:) = zero
3199 8882 : strain_inv(1,1) = 1.0
3200 8882 : strain_inv(2,2) = 1.0
3201 8882 : strain_inv(3,3) = 1.0
3202 : end if
3203 : ! write(*,*) "---- STRAIN ----"
3204 : ! do ii = 1,3
3205 : ! write(*,*) strain%strain(ii,:)
3206 : ! enddo
3207 : ! write(*,*) "---- 1+STRAIN inv ----"
3208 : ! do ii = 1,3
3209 : ! write(*,*) strain_inv(ii,:)
3210 : ! enddo
3211 :
3212 : ! fill the history position
3213 13519 : if(present(xcart_hist)) then
3214 35921998 : xcart_hist_tmp(:,:) = xcart_hist(:,:)
3215 : else
3216 961 : call xred2xcart(natom, rprimd_hist, xcart_hist_tmp, xred_hist)
3217 : end if
3218 :
3219 : ! Fill the reference position and change the cartesian coordinates
3220 : ! if the rprimd is different
3221 13519 : if(has_strain) then
3222 4637 : if(present(xcart_ref)) then
3223 4637 : call xcart2xred(natom, rprimd_ref, xcart_ref, xred_ref_tmp)
3224 4637 : call xred2xcart(natom, rprimd_hist, xcart_ref_tmp, xred_ref_tmp)
3225 : else
3226 0 : call xred2xcart(natom, rprimd_hist, xcart_ref_tmp, xred_ref)
3227 : end if
3228 : else
3229 8882 : if(present(xcart_ref)) then
3230 34705746 : xcart_ref_tmp(:,:) = xcart_ref(:,:)
3231 : else
3232 0 : call xred2xcart(natom, rprimd_ref, xcart_ref_tmp, xred_ref)
3233 : end if
3234 : end if
3235 :
3236 : ! Compute displacement
3237 13519 : if(need_displacement)then
3238 36076719 : displacement(:,:) = zero
3239 : ! write(*,*) "----- DISP -----"
3240 9029319 : do ii = 1, natom
3241 36076719 : displacement(:,ii) = xcart_hist_tmp(:,ii) - xcart_ref_tmp(:,ii)
3242 : ! write(*,*) displacement(:,ii)
3243 : end do
3244 : end if
3245 :
3246 : ! Du_Delta after Equation A4 in 2017 Paper Carlos
3247 198361119 : du_delta=zero
3248 8969039 : do ia=1,my_natom
3249 8955520 : ib = my_atoms(ia)
3250 : !Calc (1+eta)^-1 * disp(ib)
3251 116421760 : strain_inv_u = MATMUL(strain_inv,displacement(:,ib))
3252 : ! write(*,*) "--- strain_inv_u ---"
3253 : ! write(*,*) strain_inv_u
3254 : ! fill du_delta_e
3255 62702159 : do ii = 1,6
3256 223888000 : do mu = 1,3
3257 161199360 : if(alpha(ii)==mu)then
3258 53733120 : du_delta(ii,mu,ib) = du_delta(ii,mu,ib) + half * strain_inv_u(beta(ii))
3259 : end if
3260 214932480 : if(beta(ii)==mu)then
3261 53733120 : du_delta(ii,mu,ib) = du_delta(ii,mu,ib) + half * strain_inv_u(alpha(ii))
3262 : end if
3263 : enddo
3264 : enddo
3265 : enddo
3266 13519 : call xmpi_sum(du_delta , comm, ierr)
3267 :
3268 13519 : ABI_FREE(my_atoms)
3269 :
3270 27038 : end subroutine effective_potential_getDisp
3271 : !!***
3272 :
3273 : !****f* m_effective_potential/effective_potential_distributeResidualForces
3274 : !!
3275 : !! NAME
3276 : !! effective_potential_distributeResidualForces
3277 : !!
3278 : !! FUNCTION
3279 : !! Distribute the residual forces in a weighted manner
3280 : !!
3281 : !! INPUTS
3282 : !! natom = number of atoms
3283 : !! eff_pot = effective potential structure
3284 : !!
3285 : !! OUTPUT
3286 : !! fcart = forces in cartesian coordinates
3287 : !!
3288 : !! SOURCE
3289 :
3290 13519 : subroutine effective_potential_distributeResidualForces(eff_pot,fcart,natom)
3291 :
3292 : !Arguments ------------------------------------
3293 : !scalars
3294 : integer, intent(in) :: natom
3295 : !array
3296 : type(effective_potential_type),intent(in) :: eff_pot
3297 : real(dp),intent(inout) :: fcart(3,natom)
3298 : !Local variables-------------------------------
3299 : !scalar
3300 : real(dp):: mass_ia,sum_mass
3301 : integer :: ia
3302 : !array
3303 : real(dp):: sum_f(3)
3304 :
3305 : ! *************************************************************************
3306 :
3307 9029319 : sum_f(1) = sum(fcart(1,:))
3308 9029319 : sum_f(2) = sum(fcart(2,:))
3309 9029319 : sum_f(3) = sum(fcart(3,:))
3310 13519 : sum_mass = zero
3311 :
3312 9029319 : do ia=1,natom
3313 9029319 : sum_mass = sum_mass + eff_pot%crystal%amu(eff_pot%supercell%typat(ia))
3314 : end do
3315 :
3316 9029319 : do ia=1,natom
3317 9015800 : mass_ia = eff_pot%crystal%amu(eff_pot%supercell%typat(ia))
3318 36076719 : fcart(:,ia) = fcart(:,ia) - (mass_ia/sum_mass) * sum_f(:)
3319 : end do
3320 :
3321 :
3322 13519 : end subroutine effective_potential_distributeResidualForces
3323 : !!***
3324 :
3325 : !****f* m_effective_potential/effective_potential_writeAnhHead
3326 : !!
3327 : !! NAME
3328 : !! effective_potential_writeAnhHead
3329 : !!
3330 : !! FUNCTION
3331 : !! Write Header of anharmonic_energy_terms.out file
3332 : !!
3333 : !! INPUTS
3334 : !! natom = number of atoms
3335 : !! eff_pot = effective potential structure
3336 : !!
3337 : !! OUTPUT
3338 : !!
3339 : !! SOURCE
3340 :
3341 1 : subroutine effective_potential_writeAnhHead(ncoeff,filename,anh_terms)
3342 :
3343 : !Arguments ------------------------------------
3344 : !scalars
3345 : integer, intent(in) :: ncoeff
3346 : !Strings/Characters
3347 : character(len=fnlen) :: filename
3348 : !array
3349 : type(anharmonics_terms_type ),intent(in) :: anh_terms
3350 : !Local variables-------------------------------
3351 : !scalar
3352 : integer :: icoeff,unit_out
3353 : !Strings/Characters
3354 : character(len=fnlen) :: name_file
3355 : character(len=200):: term_name
3356 : !array
3357 :
3358 : ! *************************************************************************
3359 :
3360 : ! Marcus: if wanted: analyze anharmonic terms of effective potential &&
3361 : ! and print anharmonic contribution to file anharmonic_energy_terms.out
3362 : ! Open File and write header
3363 1 : name_file=trim(filename)//'_anharmonic_terms_energy.dat'
3364 1 : unit_out = get_unit()
3365 1 : open(unit=unit_out,file=name_file,status='replace',form='formatted')
3366 1 : write(unit_out,*) '#---------------------------------------------#'
3367 1 : write(unit_out,*) '# Anharmonic Terms Energy Contribution #'
3368 1 : write(unit_out,*) '#---------------------------------------------#'
3369 1 : write(unit_out,*) ''
3370 1 : write(unit_out,'(A,I5)') 'Number of Terms: ', ncoeff
3371 1 : write(unit_out,*) ''
3372 1 : write(unit_out,'(A)') 'Terms Names'
3373 61 : do icoeff=1,ncoeff
3374 60 : term_name = anh_terms%coefficients(icoeff)%name
3375 61 : write(unit_out,'(I5,A,A)') icoeff,' ',trim(term_name)
3376 : enddo
3377 1 : write(unit_out,*) ''
3378 1 : write(unit_out,'(A)',advance='no') 'Cycle/Terms'
3379 61 : do icoeff=1,ncoeff
3380 61 : if(icoeff<ncoeff)then
3381 59 : write(unit_out,'(I5)',advance='no') icoeff
3382 : else
3383 1 : write(unit_out,'(I5)',advance='yes') icoeff
3384 : endif
3385 : enddo
3386 : !close(unit_out)
3387 :
3388 1 : end subroutine effective_potential_writeAnhHead
3389 :
3390 :
3391 : !AM_EXPERIMENTAL SECTION
3392 : !!****f* m_effective_potential/equal
3393 : !! NAME
3394 : !! equal
3395 : !!
3396 : !! FUNCTION
3397 : !! compare two effective potential
3398 : !!
3399 : !! INPUTS
3400 : !! e1<type(effective_potential_type)> = effective_potential datatype
3401 : !! e2<type(effective_potential_type)> = effective_potential datatype
3402 : !!
3403 : !! OUTPUT
3404 : !!
3405 : !! SOURCE
3406 :
3407 0 : pure function effective_potential_compare(e1,e2) result (res)
3408 :
3409 : !Arguments ------------------------------------
3410 : type(effective_potential_type), intent(in) :: e1,e2
3411 : logical :: res
3412 : ! *************************************************************************
3413 0 : res = .false.
3414 : if(e1%crystal%natom==e2%crystal%natom.and.&
3415 : & e1%harmonics_terms%ifcs%nrpt==e2%harmonics_terms%ifcs%nrpt.and.&
3416 : & e1%crystal%ntypat==e2%crystal%ntypat.and.&
3417 : & e1%harmonics_terms%nqpt==e2%harmonics_terms%nqpt.and.&
3418 0 : & abs(e1%energy-e2%energy)<tol16.and.&
3419 : & abs(e1%crystal%ucvol-e2%crystal%ucvol)<tol16) then
3420 0 : res = .true.
3421 : end if
3422 :
3423 0 : end function effective_potential_compare
3424 : !!***
3425 :
3426 : ! !****f* m_effective_potential/effective_potential_effpot2ddb
3427 : ! !!
3428 : ! !! NAME
3429 : ! !! effective_potential_effpot2ddb
3430 : ! !!
3431 : ! !! FUNCTION
3432 : ! !! Convert eff_pot into ddb datatype
3433 : ! !!
3434 : ! !! INPUTS
3435 : ! !! eff_pot = effective potential datatype
3436 : ! !!
3437 : ! !! OUTPUT
3438 : ! !! ddb = ddb with all information
3439 : ! !!
3440 : ! !! SOURCE
3441 :
3442 : ! subroutine effective_potential_effpot2ddb(ddb,crystal,eff_pot,ncell,nph1l,option,qph1l)
3443 :
3444 :
3445 : ! !Arguments ------------------------------------
3446 : ! !scalars
3447 : ! integer,intent(in) :: nph1l,option
3448 : ! !array
3449 : ! integer,intent(in) :: ncell(3)
3450 : ! real(dp),intent(in):: qph1l(3,nph1l)
3451 : ! type(effective_potential_type),intent(inout) :: eff_pot
3452 : ! type(ddb_type),intent(out) :: ddb
3453 : ! type(crystal_t),intent(out) :: crystal
3454 : ! !Local variables-------------------------------
3455 : ! !scalar
3456 : ! integer :: ii,jj,msym
3457 : ! real(dp):: ucvol
3458 :
3459 : ! ! type(anaddb_dataset_type) :: inp
3460 : ! !array
3461 : ! real(dp) :: gmet(3,3),rmet(3,3)
3462 : ! real(dp) :: gprimd(3,3),rprimd(3,3)
3463 : ! real(dp),allocatable :: xred(:,:)
3464 : ! ! character :: title(eff_pot%crystal%ntypat)
3465 : ! integer,allocatable :: symrel(:,:,:),symafm(:)
3466 : ! real(dp),allocatable :: tnons(:,:)
3467 :
3468 : ! ! *************************************************************************
3469 :
3470 : ! ! Number of 2dte blocks in present object
3471 : ! ! integer,allocatable :: flg(:,:)
3472 : ! ! flg(msize,nblok)
3473 : ! ! flag to indicate presence of a given block
3474 : ! ! integer,allocatable :: typ(:)
3475 : ! ! typ(nblok)
3476 : ! ! type of each block - ddk, dde, phonon etc...
3477 : ! ! real(dp),allocatable :: amu(:)
3478 : ! ! amu(ntypat)
3479 : ! ! mass of the atoms (atomic mass unit)
3480 : ! ! real(dp),allocatable :: nrm(:,:)
3481 : ! ! nrm(3,nblok)
3482 : ! ! norm of the q-points for each block - can be 0 to indicate a direction of approach to gamma
3483 : ! ! real(dp),allocatable :: qpt(:,:)
3484 : ! ! qpt(9,nblok)
3485 : ! ! q-point vector in reciprocal space (reduced lattice coordinates) for each block
3486 : ! ! real(dp),allocatable :: val(:,:,:)
3487 : ! ! val(2,msize,nblok)
3488 : ! ! values of the second energy derivatives in each block
3489 :
3490 : ! ! Useless value
3491 : ! ddb%nblok = -1
3492 :
3493 : ! !option = 1 just print ddb for 1 1 1 cell
3494 : ! if(option==1) then
3495 : ! ! Compute different matrices in real and reciprocal space, also
3496 : ! ! checks whether ucvol is positive.
3497 : ! call metric(gmet,gprimd,-1,rmet,eff_pot%crystal%rprimd,ucvol)
3498 :
3499 : ! ! Convert to rprimd
3500 : ! do ii=1,3
3501 : ! do jj=1,3
3502 : ! rprimd(ii,jj)=eff_pot%crystal%rprimd(ii,jj)
3503 : ! end do
3504 : ! end do
3505 :
3506 : ! ! Obtain reciprocal space primitive transl g from inverse trans of r
3507 : ! ! (Unlike in abinit, gprim is used throughout ifc; should be changed, later)
3508 : ! call matr3inv(rprimd,gprimd)
3509 :
3510 : ! ! transfert basic values
3511 : ! ddb%natom = eff_pot%crystal%natom
3512 : ! ddb%mpert = ddb%natom+6
3513 : ! ddb%msize = 3*ddb%mpert*3*ddb%mpert;
3514 : ! ddb%ntypat = eff_pot%crystal%ntypat
3515 : ! ddb%occopt = 3 ! default value
3516 : ! ddb%prtvol = 0 ! default value
3517 : ! ddb%rprim = rprimd ! dimensioless real space primitive vectors
3518 : ! ddb%gprim = gprimd ! dimensioless reciprocal space primitive vectors
3519 : ! ddb%acell = one
3520 : ! msym = 1
3521 : ! ! Setup crystal type
3522 : ! ABI_MALLOC(xred,(3,ddb%natom))
3523 : ! ! call xcar2xred(ddb%natom,eff_pot%crystal%rprimd,eff_pot%crystal%xcart,xred)
3524 : ! !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
3525 : ! ABI_MALLOC(symrel,(3,3,msym))
3526 : ! ABI_MALLOC(symafm,(msym))
3527 : ! ABI_MALLOC(tnons,(3,msym))
3528 :
3529 : ! ! call crystal_init(ddb%amu,crystal,1,ddb%natom,size(eff_pot%crystal%znucl),eff_pot%crystal%ntypat,1,&
3530 : ! !& eff_pot%crystal%rprimd,eff_pot%crystal%typat,xred,eff_pot%crystal%znucl,&
3531 : ! !& eff_pot%crystal%znucl,0,.FALSE.,.FALSE.,title)!,&
3532 : ! !& symrel=symrel,tnons=tnons,symafm=symafm)
3533 : ! ! call crystal_print(crystal)
3534 : ! ! stop
3535 : ! !TEST_AM
3536 : ! ABI_FREE(symrel)
3537 : ! ABI_FREE(symafm)
3538 : ! ABI_FREE(tnons)
3539 :
3540 : ! ABI_FREE(xred)
3541 :
3542 : ! else if (option==2) then
3543 : ! ! Compute different matrices in real and reciprocal space, also
3544 : ! ! checks whether ucvol is positive.
3545 : ! call metric(gmet,gprimd,-1,rmet,eff_pot%supercell%rprimd,ucvol)
3546 :
3547 : ! ! Convert to rprim (dimensionless)
3548 : ! do ii=1,3
3549 : ! do jj=1,3
3550 : ! rprimd(ii,jj)=eff_pot%supercell%rprimd(ii,jj)
3551 : ! end do
3552 : ! end do
3553 :
3554 : ! ! Obtain reciprocal space primitive transl g from inverse trans of r
3555 : ! ! (Unlike in abinit, gprim is used throughout ifc; should be changed, later)
3556 : ! call matr3inv(rprimd,gprimd)
3557 :
3558 : ! ! transfert basic values
3559 : ! ddb%natom = eff_pot%supercell%natom
3560 : ! ddb%ntypat = eff_pot%crystal%ntypat
3561 : ! ddb%mpert = ddb%natom+6
3562 : ! ddb%msize = 3*ddb%mpert*3*ddb%mpert;
3563 : ! ddb%occopt = 3 ! default value
3564 : ! ddb%prtvol = 0 ! default value
3565 : ! ddb%rprim = rprimd ! dimensioless real space primitive vectors
3566 : ! ddb%gprim = gprimd ! dimensioless reciprocal space primitive vectors
3567 : ! ddb%acell = one
3568 :
3569 : ! end if
3570 : ! !TEST_AM
3571 : ! !write(std_out,*)"natom ",ddb%natom
3572 : ! !write(std_out,*)"ntypat",ddb%ntypat
3573 : ! !write(std_out,*)"mpert",ddb%mpert
3574 : ! !write(std_out,*)"msize",ddb%msize
3575 : ! !write(std_out,*)"occopt",ddb%occopt
3576 : ! !write(std_out,*)"prtvol",ddb%prtvol
3577 : ! !write(std_out,*)"rprim",ddb%rprim
3578 : ! !write(std_out,*)"gprim",ddb%gprim
3579 : ! !write(std_out,*)"acell",ddb%acell
3580 : ! !TEST_AM
3581 :
3582 : ! end subroutine effective_potential_effpot2ddb
3583 : ! !!***
3584 :
3585 :
3586 : ! !****f* m_effective_potential/effective_potential_printPDOS
3587 : ! !!
3588 : ! !! NAME
3589 : ! !! effective_potential_printPDOS
3590 : ! !!
3591 : ! !! FUNCTION
3592 : ! !! Apply the acoustic sum rule on the effective potential
3593 : ! !!
3594 : ! !! INPUTS
3595 : ! !! eff_pot = effective potential datatype
3596 : ! !! option = 0 (default) do nothing
3597 : ! !! = 1 print PHFRQ for specific qgrid (need nph1l and qph1l)
3598 : ! !! = 2 print PHFRQ for supercell (q=gamma) (need nph1l and qph1l and ncell)
3599 : ! !! OUTPUT
3600 : ! !! eff_pot
3601 : ! !!
3602 : ! !! SOURCE
3603 :
3604 : ! subroutine effective_potential_printPDOS(eff_pot,filename,ncell,nph1l,option,qph1l)
3605 :
3606 : ! !Arguments ------------------------------------
3607 : ! !scalars
3608 : ! integer,intent(in) :: nph1l,option
3609 : ! !array
3610 : ! integer,intent(in) :: ncell(3)
3611 : ! real(dp),intent(in):: qph1l(3,nph1l)
3612 : ! type(effective_potential_type),intent(inout) :: eff_pot
3613 : ! character(len=fnlen),intent(in) :: filename
3614 : ! !Local variables-------------------------------
3615 : ! !scalar
3616 : ! ! integer :: lenstr
3617 : ! ! character(len=strlen) :: string
3618 : ! !array
3619 : ! type(crystal_t) :: Crystal
3620 : ! ! type(anaddb_dataset_type) :: inp
3621 : ! type(ddb_type) :: ddb
3622 : ! ! type(asrq0_t) :: asrq0
3623 :
3624 : ! ! *************************************************************************
3625 :
3626 : ! if (option > 0) then
3627 :
3628 : ! ! First: transfer into ddb datatype:
3629 : ! call effective_potential_effpot2ddb(ddb,Crystal,eff_pot,ncell,nph1l,option,qph1l)
3630 :
3631 : ! ! Setup fake anaddb_dataset
3632 : ! ! string = ''
3633 : ! ! lenstr = 0
3634 : ! ! call invars9(inp,lenstr,ddb%natom,string)
3635 : ! ! fill it with multibinit_dataset values
3636 : ! ! inp%prt_ifc = 1
3637 : ! ! inp%ifcflag = 1
3638 : ! ! inp%qph1l = qph1l
3639 : ! ! inp%nph1l = nph1l
3640 :
3641 : ! ! ! In case the interatomic forces are not calculated, the
3642 : ! ! ! ASR-correction (asrq0%d2asr) has to be determined here from the Dynamical matrix at Gamma.
3643 : ! ! if (inp%ifcflag == 0) then
3644 : ! ! asrq0 = ddb_get_asrq0(ddb, inp%asr, inp%rfmeth, crystal%xcart)
3645 : ! ! end if
3646 :
3647 : ! ! !MG: Note that I'm passing xmpi_comm_self here.
3648 : ! ! call mkphbs(eff_pot%harmonics_terms%ifcs,Crystal,inp,ddb,asrq0,filename,xmpi_comm_self)
3649 :
3650 : ! ! call asrq0_free(asrq0)
3651 :
3652 : ! end if
3653 :
3654 : ! end subroutine effective_potential_printPDOS
3655 : ! !!***
3656 :
3657 : !****f* m_effective_potential/effective_potential_computeGradient
3658 : !!
3659 : !! NAME
3660 : !! effective_potential_computeGradient
3661 : !!
3662 : !! FUNCTION
3663 : !! Compute finate differences on forces to compute dynmical matrix
3664 : !! at gamma for supercell
3665 : !!
3666 : !! INPUTS
3667 : !! eff_pot = effective potential datatype
3668 : !!
3669 : !! OUTPUT
3670 : !! dynmat = ddb with all information
3671 : !!
3672 : !! SOURCE
3673 :
3674 0 : subroutine effective_potential_computeGradient(delta,fcart_out,eff_pot,natom,ncell,option,comm)
3675 :
3676 : !Arguments ------------------------------------
3677 : !scalars
3678 : integer,intent(in) :: natom,option,comm
3679 : real(dp),intent(in) :: delta
3680 : !array
3681 : integer,intent(in) :: ncell(3)
3682 : type(effective_potential_type),intent(inout) :: eff_pot
3683 : real(dp),intent(out) :: fcart_out(3,natom)
3684 : !Local variables-------------------------------
3685 : !scalar
3686 : character(len=500) :: msg
3687 : integer :: ia,ib,ii,mu,nu,npt
3688 : real(dp):: delt,energy
3689 : !array
3690 : real(dp):: strten(6)
3691 0 : real(dp),allocatable :: disp(:,:),diff(:)
3692 0 : real(dp),allocatable :: gred(:,:),fcart(:,:),xred(:,:)
3693 :
3694 : ! *************************************************************************
3695 :
3696 : !Do Some check
3697 0 : do ii=1,3
3698 0 : if(abs(eff_pot%supercell%qphon(ii)-ncell(ii)) > tol16)then
3699 0 : call effective_potential_setSupercell(eff_pot,comm,ncell)
3700 : end if
3701 : end do
3702 :
3703 0 : write(msg,'(a,(80a),3a)') ch10,('-',ii=1,80),ch10,' Generation of the dynmical matrix by ',&
3704 0 : & 'finite differences'
3705 : ! call wrtout(ab_out,msg,'COLL')
3706 : ! call wrtout(std_out,msg,'COLL')
3707 :
3708 0 : npt = 3 !default
3709 0 : select case (option)
3710 : case (1)
3711 : ! write(msg,'(2a)') ch10,' Finite differences on 1 points '
3712 0 : npt = 2
3713 : case (2)
3714 : ! write(msg,'(2a)') ch10,' Finite differences on 3 points '
3715 0 : npt = 3
3716 : case (3)
3717 : ! write(msg,'(2a)') ch10,' Finite differences on 5 points '
3718 0 : npt = 5
3719 : end select
3720 :
3721 : ! call wrtout(ab_out,msg,'COLL')
3722 : ! call wrtout(std_out,msg,'COLL')
3723 :
3724 : ! Allocation of forces arrays
3725 :
3726 0 : ABI_MALLOC(disp,(3,natom))
3727 0 : ABI_MALLOC(diff,(npt))
3728 0 : ABI_MALLOC(gred,(3,natom))
3729 0 : ABI_MALLOC(fcart,(3,natom))
3730 0 : ABI_MALLOC(xred,(3,natom))
3731 :
3732 0 : fcart_out = zero
3733 :
3734 : call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3735 0 : & eff_pot%supercell%xcart,xred)
3736 :
3737 0 : do ia=1,eff_pot%supercell%natom
3738 0 : do mu=1,3
3739 0 : diff = zero
3740 0 : do ii=1,npt
3741 0 : delt = (-(npt/2+1)+ii) * delta
3742 0 : disp = zero
3743 0 : disp(mu,ia) = delt * eff_pot%supercell%rprimd(mu,mu)
3744 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,&
3745 : & strten,natom,eff_pot%supercell%rprimd,&
3746 : & displacement=disp,&
3747 0 : & compute_anharmonic=.FALSE.,verbose=.false.)
3748 :
3749 : ! diff(ii,:,:) = gred(:,:)
3750 0 : diff(ii) = energy
3751 : end do
3752 :
3753 0 : select case (option)
3754 : case (1)
3755 0 : fcart_out(mu,ia) = (diff(1)-diff(2)) / (delta)
3756 : case (2)
3757 0 : fcart_out(mu,ia) = (diff(3)-diff(1)) / (2*delta)
3758 : case (3)
3759 0 : fcart_out(mu,ia) = (-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta)
3760 : end select
3761 : end do
3762 : end do
3763 :
3764 : !TEST_AM
3765 : !Write the phonon into ddb format wavevector
3766 : !write(999, '(a,3es16.8,f6.1)' )' qpt',real((/0,0,0/),dp),1.0
3767 : !Write the matrix elements
3768 0 : do ib=1,eff_pot%supercell%natom
3769 : do nu=1,3
3770 : do ia=1,eff_pot%supercell%natom
3771 : do mu=1,3
3772 : ! write(999,'(4i4,2d22.14)')nu,ib,mu,ia,&
3773 : ! & dynmat(1,mu,ia,nu,ib),dynmat(2,mu,ia,nu,ib)
3774 : end do
3775 : end do
3776 : end do
3777 : end do
3778 : !TEST_AM
3779 :
3780 : ! Deallocation of arrays
3781 0 : ABI_FREE(disp)
3782 0 : ABI_FREE(diff)
3783 0 : ABI_FREE(gred)
3784 0 : ABI_FREE(fcart)
3785 0 : ABI_FREE(xred)
3786 :
3787 :
3788 0 : end subroutine effective_potential_computeGradient
3789 : !!***
3790 :
3791 : !!****f* m_effective_potential/effective_potential_checkDEV
3792 : !!
3793 : !! NAME
3794 : !! effective_potential_checkDEV
3795 : !!
3796 : !! FUNCTION
3797 : !! Routine for develloper Check by finite differences the equations in
3798 : !! effective_potential_evaluate need to provide HIST file, so you need to
3799 : !! activate the fit_process or bound_process to activate the reading of the HIST
3800 : !!
3801 : !! INPUTS
3802 : !! eff_pot<type(effective_potential)> = effective potential
3803 : !! hist<type(abihist)> = The history of the MD
3804 : !! natom = number of atom
3805 : !! ntime = number of time in the hist
3806 : !!
3807 : !! OUTPUT
3808 : !!
3809 : !! SOURCE
3810 :
3811 0 : subroutine effective_potential_checkDEV(eff_pot,hist,natom,ntime)
3812 :
3813 : !Arguments ------------------------------------
3814 : !scalars
3815 : integer, intent(in) :: natom,ntime
3816 : !arrays
3817 : type(effective_potential_type),intent(in) :: eff_pot
3818 : type(abihist),intent(in) :: hist
3819 : !Local variables-------------------------------
3820 : !scalar
3821 : integer :: ii,jj,ia,mu,npt,istep
3822 : ! integer :: ifirst
3823 : real(dp):: energy,delt,delta,ucvol
3824 : !arrays
3825 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),mat_def(3,3),identity(3,3)
3826 0 : real(dp):: fcart(3,natom),gred(3,natom),strten(6),rprimd(3,3)
3827 : real(dp):: rprimd_def(3,3),rprimd_ref(3,3),deltalist(5)
3828 0 : real(dp):: disp(3,natom),disp_red(3,natom),strain(6),du_delta(6,3,natom),diff(5)
3829 0 : real(dp),allocatable :: xred(:,:)
3830 : integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
3831 : character(len=500) :: msg
3832 :
3833 : ! *************************************************************************
3834 :
3835 : !Do some checks
3836 0 : if(ntime /= hist%mxhist)then
3837 0 : write(msg,'(a)')'ntime is not correct'
3838 0 : ABI_BUG(msg)
3839 : end if
3840 :
3841 0 : if(natom /= size(hist%xred,2)) then
3842 0 : write(msg,'(a)')'natom is not correct'
3843 0 : ABI_BUG(msg)
3844 : end if
3845 :
3846 :
3847 0 : ABI_MALLOC(xred,(3,natom))
3848 0 : xred = zero
3849 :
3850 : !option 1 => set the reference for the test
3851 : ! call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3852 : !& eff_pot%supercell%xcart,xred)
3853 : ! rprimd = eff_pot%supercell%rprimd
3854 :
3855 : !option 2 => set a specific step for the test
3856 0 : istep = 4
3857 0 : xred = hist%xred(:,:,istep)
3858 0 : rprimd = hist%rprimd(:,:,istep)
3859 :
3860 : rprimd_ref = eff_pot%supercell%rprimd
3861 0 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
3862 :
3863 0 : npt=5
3864 0 : delta = 0.001
3865 0 : deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/)
3866 : strain = zero
3867 :
3868 0 : do ia=1,natom
3869 0 : do mu=1,3
3870 0 : write(std_out,*) "atm: ",ia," dir: ",mu
3871 0 : do ii=1,npt
3872 0 : delt = deltalist(ii)
3873 :
3874 : ! Get the initial displacement
3875 : call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3876 : & eff_pot%supercell%rprimd,1,xred_hist=xred,&
3877 : & xcart_ref=eff_pot%supercell%xcart,&
3878 0 : & compute_displacement = .true.,compute_duDelta = .true.)
3879 :
3880 : ! Add the delta
3881 0 : call xcart2xred(natom, rprimd, disp, disp_red)
3882 0 : disp_red(mu,ia) = disp_red(mu,ia) + delt
3883 0 : call xred2xcart(natom, rprimd, disp, disp_red)
3884 :
3885 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3886 : & xred=xred,du_delta=du_delta,&
3887 0 : & displacement=disp,compute_anharmonic=.true.,verbose=.false.)
3888 0 : diff(ii) = energy
3889 :
3890 : end do
3891 :
3892 : ! Get the initial displacement
3893 : call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3894 : & eff_pot%supercell%rprimd,1,xred_hist=xred,&
3895 : & xcart_ref=eff_pot%supercell%xcart,&
3896 0 : & compute_displacement = .true.,compute_duDelta = .true.)
3897 :
3898 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3899 : & xred=xred,du_delta=du_delta,&
3900 0 : & displacement=disp,compute_anharmonic=.true.,verbose=.false.)
3901 :
3902 0 : write(std_out,*) "Analyti:",gred(mu,ia)
3903 0 : write(std_out,*) "FD :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta)
3904 0 : write(std_out,*) "Diff(%):",abs(100*(gred(mu,ia)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))&
3905 0 : & / (12*delta) )) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) ))
3906 :
3907 : end do
3908 : end do
3909 :
3910 :
3911 : ! Fill the identity matrix
3912 0 : identity = zero
3913 0 : forall(ii=1:3)identity(ii,ii)=1
3914 :
3915 0 : npt=5
3916 0 : delta = 0.0005
3917 0 : deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/)
3918 :
3919 0 : do jj=1,6
3920 0 : write(std_out,*) "strain ",jj
3921 0 : do ii=1,npt
3922 : strain = zero
3923 0 : delt = deltalist(ii)
3924 0 : mat_def = zero
3925 0 : strain(jj) = strain(jj) + delt
3926 :
3927 0 : mat_def(alpha(jj),beta(jj)) = mat_def(alpha(jj),beta(jj)) + half * strain(jj)
3928 0 : mat_def(beta(jj),alpha(jj)) = mat_def(beta(jj),alpha(jj)) + half * strain(jj)
3929 :
3930 0 : mat_def = mat_def + identity
3931 0 : rprimd_def = matmul(mat_def,rprimd)
3932 :
3933 : ! The two options should give the same result
3934 : ! Option 1 => compute the disps and provide them to evaluate
3935 : ! call effective_potential_getDisp(disp,du_delta,natom,rprimd_def,&
3936 : ! & rprimd_ref,1,xred_hist=xred,&
3937 : ! & xcart_ref=eff_pot%supercell%xcart,&
3938 : ! & compute_displacement = .true.,compute_duDelta = .true.)
3939 :
3940 : ! call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd_def,&
3941 : ! & xred=xred,du_delta=du_delta,&
3942 : ! & displacement=disp,strain=strain,&
3943 : ! & compute_anharmonic=.true.,verbose=.false.)
3944 :
3945 : ! Option 2 => compute the disps within evaluate
3946 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd_def,&
3947 0 : & xred=xred,compute_anharmonic=.true.,verbose=.false.)
3948 :
3949 :
3950 :
3951 0 : diff(ii) = energy
3952 :
3953 : end do
3954 :
3955 : ! The two options should give the same result
3956 : ! Option 1 => compute the disps and provide them to evaluate
3957 : ! call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3958 : ! & rprimd_ref,1,xred_hist=xred,&
3959 : ! & xcart_ref=eff_pot%supercell%xcart,&
3960 : ! & compute_displacement = .true.,compute_duDelta = .true.)
3961 :
3962 : ! call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3963 : ! & xred=xred,du_delta=du_delta,&
3964 : ! & displacement=disp,&
3965 : ! & compute_anharmonic=.true.,verbose=.false.)
3966 :
3967 : ! Option 2 => compute the disps within evaluate
3968 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3969 0 : & xred=xred,compute_anharmonic=.true.,verbose=.false.)
3970 :
3971 0 : write(std_out,*) "Analyti:",strten(jj)
3972 0 : write(std_out,*) "FD :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol
3973 0 : write(std_out,*) "Diff(%):",abs(100*(strten(jj)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))&
3974 0 : & / (12*delta) / ucvol)) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol))
3975 :
3976 : end do
3977 :
3978 0 : end subroutine effective_potential_checkDEV
3979 : !!***
3980 :
3981 : !!****f* m_effective_potential/effective_potential_writeNETCDF
3982 : !! NAME
3983 : !! effective_potential_writeNETCDF
3984 : !!
3985 : !! FUNCTION
3986 : !! This routine print the effective potential into netcdf format
3987 : !! Several options are available
3988 : !!
3989 : !! INPUTS
3990 : !! filename = the name of output file
3991 : !! eff_pot = datatype contains the effective potential
3992 : !! option = option for the format of the xml file
3993 : !! 1 print the xml for a system
3994 : !!
3995 : !! OUTPUT
3996 : !!
3997 : !! SOURCE
3998 :
3999 0 : subroutine effective_potential_writeNETCDF(eff_pot,option,filename)
4000 :
4001 : !Arguments ------------------------------------
4002 : !scalars
4003 : integer, intent(in) :: option
4004 : character(len=fnlen),optional,intent(in) :: filename
4005 : !arrays
4006 : type(effective_potential_type), intent(in) :: eff_pot
4007 :
4008 : !Local variables-------------------------------
4009 : !scalar
4010 : integer :: amu_id,bec_id,ifccell_id,epsinf_id,elastic_id
4011 : integer :: ifc_id,ifcs_id,natom_id,ntypat_id,nrpt_id,npsp_id,typat_id
4012 : integer :: six_id,two_id,xyz_id,znucl_id
4013 : integer :: ncerr,ncid,npsp
4014 : integer :: dimCids(2),dimEids(2),dimIids(6),dimPids(1),dimRids(2),dimXids(2)
4015 : integer :: etotal_id,rprimd_id,xcart_id
4016 : character(len=500) :: msg
4017 : character(len=fnlen) :: namefile
4018 : !arrays
4019 : real(dp) :: strain(9,6)
4020 : ! *************************************************************************
4021 :
4022 : strain(:,1) = (/1,0,0,0,0,0,0,0,0/)
4023 : strain(:,2) = (/0,0,0,0,1,0,0,0,0/)
4024 : strain(:,3) = (/0,0,0,0,0,0,0,0,1/)
4025 : strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/)
4026 : strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/)
4027 : strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/)
4028 :
4029 : !Print only the reference system in xml format
4030 0 : if (option == 1) then
4031 :
4032 0 : if(present(filename)) then
4033 0 : namefile=filename
4034 : else
4035 0 : namefile='ref.nc'
4036 : end if
4037 :
4038 0 : call isfile(namefile,'new')
4039 :
4040 0 : write(msg,'(a,a,a)')ch10,&
4041 0 : & ' Generation of the xml file for the reference structure in ',trim(namefile)
4042 :
4043 0 : call wrtout(ab_out,msg,'COLL')
4044 0 : call wrtout(std_out,msg,'COLL')
4045 :
4046 : ! 1. Create netCDF file
4047 0 : ncerr = nf90_create(path=trim(namefile),cmode=NF90_CLOBBER, ncid=ncid)
4048 0 : NCF_CHECK_MSG(ncerr,"create netcdf history file")
4049 :
4050 : ! 2. Define dimensions
4051 0 : ncerr = nf90_def_dim(ncid,"natom",eff_pot%crystal%natom,natom_id)
4052 0 : NCF_CHECK_MSG(ncerr," define dimension natom")
4053 :
4054 0 : ncerr = nf90_def_dim(ncid,"ntypat",eff_pot%crystal%ntypat,ntypat_id)
4055 0 : NCF_CHECK_MSG(ncerr," define dimension ntypat")
4056 :
4057 0 : ncerr = nf90_def_dim(ncid,"nrpt",eff_pot%harmonics_terms%ifcs%nrpt,nrpt_id)
4058 0 : NCF_CHECK_MSG(ncerr," define dimension ntypat")
4059 :
4060 0 : ncerr = nf90_def_var(ncid, "typat", NF90_INT, natom_id, typat_id)
4061 0 : NCF_CHECK_MSG(ncerr," define variable typat")
4062 :
4063 0 : npsp = size(eff_pot%crystal%znucl)
4064 0 : if (npsp /= eff_pot%crystal%ntypat) then
4065 0 : ABI_WARNING("HIST file does not support alchemical mixing!")
4066 : end if
4067 0 : ncerr = nf90_def_dim(ncid,"npsp",npsp,npsp_id)
4068 0 : NCF_CHECK_MSG(ncerr," define dimension npsp")
4069 :
4070 0 : ncerr = nf90_def_var(ncid, "znucl", NF90_DOUBLE, npsp_id, znucl_id)
4071 0 : NCF_CHECK_MSG(ncerr," define variable znucl")
4072 :
4073 0 : ncerr = nf90_def_dim(ncid,"xyz",3,xyz_id)
4074 0 : NCF_CHECK_MSG(ncerr," define dimension xyz")
4075 :
4076 0 : ncerr = nf90_def_dim(ncid,"six",6,six_id)
4077 0 : NCF_CHECK_MSG(ncerr," define dimension six")
4078 :
4079 0 : ncerr = nf90_def_dim(ncid,"two",2,two_id)
4080 0 : NCF_CHECK_MSG(ncerr," define dimension two")
4081 :
4082 : ! Dimensions for xcart,xred,fcart,gred and vel
4083 0 : dimXids = (/ xyz_id, natom_id /)
4084 : ! Dimensions for rprimd
4085 0 : dimRids = (/ xyz_id, xyz_id /)
4086 : ! Dimensions for ifc
4087 0 : dimIids = (/ 2, xyz_id, natom_id, xyz_id, natom_id, nrpt_id /)
4088 : ! Dimensions for position
4089 0 : dimPids = (/ xyz_id /)
4090 : ! Dimension for elastic constant
4091 0 : dimEids = (/six_id,six_id/)
4092 : ! Dimension for cell
4093 0 : dimCids = (/nrpt_id,3/)
4094 :
4095 : ! 3. Define variables and their attributes (units and mnemonics)
4096 : call ab_define_var(ncid, (/1/), etotal_id, NF90_DOUBLE,&
4097 0 : & "energy","Energy of the reference structure","Ha" )
4098 :
4099 : call ab_define_var(ncid, dimRids, rprimd_id, NF90_DOUBLE,&
4100 0 : & "rprimd","Real space PRIMitive translations, Dimensional","bohr" )
4101 :
4102 : call ab_define_var(ncid, dimRids, epsinf_id, NF90_DOUBLE,&
4103 0 : & "epsilon_inf","Dielectric tensor, Dimensional","epsilon_inf" )
4104 :
4105 : call ab_define_var(ncid, dimEids, elastic_id, NF90_DOUBLE,&
4106 0 : & "elastic","Elastic Constants, Dimensional","Ha" )
4107 :
4108 : call ab_define_var(ncid, dimRids, bec_id, NF90_DOUBLE,&
4109 0 : & "bec","Born Effective Charges, Dimensional","abs(e)" )
4110 :
4111 : call ab_define_var(ncid, dimXids, xcart_id, NF90_DOUBLE,&
4112 0 : & "xcart","vectors (X) of atom positions in CARTesian coordinates","bohr" )
4113 :
4114 : call ab_define_var(ncid, dimIids, ifcs_id, NF90_DOUBLE,&
4115 0 : & "IFCs","Interatomic Forces Constantes in real spaces (short range), Dimensional","Hatree/bohr**2" )
4116 :
4117 : call ab_define_var(ncid, dimIids, ifc_id, NF90_DOUBLE,&
4118 0 : & "IFC","Interatomic Forces Constantes in real spaces (total range), Dimensional","Hatree/bohr**2" )
4119 :
4120 : call ab_define_var(ncid, dimCids, ifccell_id, NF90_DOUBLE,&
4121 0 : & "cell","cell for the ifc, Dimensional","Dimensionless" )
4122 :
4123 : call ab_define_var(ncid, [ntypat_id], amu_id, NF90_DOUBLE,&
4124 0 : & "amu","Masses of each type of atom in atomic mass units", "" )
4125 :
4126 : ! 4. End define mode
4127 0 : ncerr = nf90_enddef(ncid)
4128 0 : NCF_CHECK_MSG(ncerr," end define mode")
4129 :
4130 : ! 5. Write variables
4131 0 : ncerr = nf90_put_var(ncid,etotal_id, eff_pot%energy)
4132 0 : NCF_CHECK_MSG(ncerr," write variable energy")
4133 :
4134 0 : ncerr = nf90_put_var(ncid,rprimd_id, eff_pot%crystal%rprimd)
4135 0 : NCF_CHECK_MSG(ncerr," write variable rprimd")
4136 :
4137 0 : ncerr = nf90_put_var(ncid,epsinf_id, eff_pot%harmonics_terms%epsilon_inf)
4138 0 : NCF_CHECK_MSG(ncerr," write variable epsilon_inf")
4139 :
4140 0 : ncerr = nf90_put_var(ncid,elastic_id , eff_pot%harmonics_terms%elastic_constants)
4141 0 : NCF_CHECK_MSG(ncerr," write variable elastic_constant")
4142 :
4143 0 : ncerr = nf90_put_var(ncid, bec_id, eff_pot%harmonics_terms%zeff)
4144 0 : NCF_CHECK_MSG(ncerr," write variable bec")
4145 :
4146 0 : ncerr = nf90_put_var(ncid,xcart_id, eff_pot%crystal%xcart)
4147 0 : NCF_CHECK_MSG(ncerr," write variable xcart")
4148 :
4149 0 : ncerr = nf90_put_var(ncid,ifccell_id, eff_pot%harmonics_terms%ifcs%cell)
4150 0 : NCF_CHECK_MSG(ncerr," write variable cell")
4151 :
4152 0 : ncerr = nf90_put_var(ncid,ifcs_id, eff_pot%harmonics_terms%ifcs%short_atmfrc)
4153 0 : NCF_CHECK_MSG(ncerr," write variable short ifc")
4154 :
4155 0 : ncerr = nf90_put_var(ncid,ifc_id, eff_pot%harmonics_terms%ifcs%atmfrc)
4156 0 : NCF_CHECK_MSG(ncerr," write variable total ifc")
4157 :
4158 :
4159 : ! 6. Close NetCDF file
4160 0 : ncerr = nf90_close(ncid)
4161 0 : NCF_CHECK_MSG(ncerr," close netcdf history file")
4162 : end if
4163 :
4164 0 : end subroutine effective_potential_writeNETCDF
4165 : !!***
4166 :
4167 : !AM_END_EXPERIMENTAL SECTION
4168 :
4169 1658 : end module m_effective_potential
4170 : !!***
|
