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