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