Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_solver
2 : : !! NAME
3 : : !! m_tdep_solver
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the TDEP Solver data type
7 : : !! which fits the IFC on the forces and displacements of the sampling.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_tdep_solver
24 : :
25 : : use defs_basis
26 : : use m_errors
27 : : use m_abicore
28 : : use m_xmpi
29 : : use m_io_tools
30 : : use m_numeric_tools
31 : : use m_linalg_interfaces
32 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
33 : : use m_tdep_shell, only : Shell_type
34 : : use m_tdep_sym, only : Symmetries_type
35 : : use m_tdep_sampling, only : tdep_Sampling_type
36 : : use m_tdep_model, only : tdep_Model_type
37 : : use m_tdep_phi2, only : tdep_build_phi2_33
38 : : use m_tdep_constraints, only : Constraints_type, tdep_calc_orthonorm
39 : : use m_tdep_supercell, only : tdep_Supercell_type
40 : :
41 : : implicit none
42 : : !!***
43 : :
44 : : !!****t* m_tdep_solver/tdep_Solver_type
45 : : !! NAME
46 : : !! tdep_Solver_type
47 : : !!
48 : : !! FUNCTION
49 : : !! Holds the linear system to be solved in order to fit the IFC: the matrix
50 : : !! of the displacement products (fcoeff), the forces, the constraints on the
51 : : !! IFC coefficients, and the solution (theta).
52 : : !!
53 : : !! SOURCE
54 : :
55 : : type tdep_Solver_type
56 : :
57 : : integer :: order
58 : : ! Maximum IFC order (2, 3, or 4)
59 : :
60 : : integer :: natom
61 : : ! Number of atoms in the supercell
62 : :
63 : : integer :: natom_unitcell
64 : : ! Number of atoms in the unitcell
65 : :
66 : : integer :: my_nstep
67 : : ! Number of MD steps
68 : :
69 : : integer :: ntotcoeff
70 : : ! Total number of coefficients
71 : :
72 : : integer :: ntotconst
73 : : ! Total number of constraints
74 : :
75 : : integer :: ncoeff1st
76 : : ! Number of 1st order coefficients
77 : :
78 : : integer :: ncoeff2nd
79 : : ! Number of 2nd order coefficients
80 : :
81 : : integer :: ncoeff3rd
82 : : ! Number of 3rd order coefficients
83 : :
84 : : integer :: ncoeff4th
85 : : ! Number of 4th order coefficients
86 : :
87 : : integer :: nconst_1st
88 : : ! Number of 1st order constraints
89 : :
90 : : integer :: nconst_2nd
91 : : ! Number of 2nd order constraints
92 : :
93 : : integer :: nconst_3rd
94 : : ! Number of 3rd order constraints
95 : :
96 : : integer :: nconst_4th
97 : : ! Number of 4th order constraints
98 : :
99 : : integer :: nconst_rot2nd
100 : : ! Number of rotational invariance constraints at the 2nd order
101 : :
102 : : integer :: nconst_huang
103 : : ! Number of Huang invariance constraints
104 : :
105 : : integer :: nconst_dynmat
106 : : ! Number of constraints enforcing the symmetry of the dynamical matrix
107 : :
108 : : integer :: nconst_rot3rd
109 : : ! Number of rotational invariance constraints at the 3rd order
110 : :
111 : : integer :: nconst_asr3rd
112 : : ! Number of acoustic sum rule constraints at the 3rd order
113 : :
114 : : integer :: nconst_rot4th
115 : : ! Number of rotational invariance constraints at the 4th order
116 : :
117 : : integer :: nconst_asr4th
118 : : ! Number of acoustic sum rule constraints at the 4th order
119 : :
120 : : double precision, allocatable :: fcoeff(:,:)
121 : : ! fcoeff(3*natom*my_nstep,ntotcoeff)
122 : : ! All the cartesian displacement matrices
123 : : ! at every order (u, u*u, u*u*u, u*u*u*u).
124 : :
125 : : double precision, allocatable :: const(:,:)
126 : : ! const(ntotconst,ntotcoeff)
127 : : ! The constraint matrices.
128 : :
129 : : double precision, allocatable :: theta(:)
130 : : ! theta(ntotcoeff)
131 : : ! The IFC coefficients at all orders, as a flat array.
132 : :
133 : : double precision, allocatable :: Forces(:)
134 : : ! Forces(3*natom*my_nstep)
135 : : ! The cartesian forces for all configurations, as a flat array.
136 : : ! These are weighted by the number of configurations.
137 : : ! In case of dipole-dipole interaction, the long-range part of the
138 : : ! forces should be removed.
139 : :
140 : : contains
141 : :
142 : : procedure :: init => tdep_solver_init
143 : : procedure :: free => tdep_solver_free
144 : : procedure :: set_residual_forces => tdep_solver_set_residual_forces
145 : : procedure :: calc_phi1fcoeff => tdep_calc_phi1fcoeff
146 : : procedure :: calc_phi2fcoeff => tdep_calc_phi2fcoeff
147 : : procedure :: calc_phi3fcoeff => tdep_calc_phi3fcoeff
148 : : procedure :: calc_phi4fcoeff => tdep_calc_phi4fcoeff
149 : : procedure :: calc_MoorePenrose => tdep_calc_MoorePenrose
150 : : procedure :: calc_constraints => tdep_calc_constraints
151 : : procedure :: calc_phi1 => tdep_calc_phi1
152 : : procedure :: calc_phi2 => tdep_calc_phi2
153 : : procedure :: calc_phi3ref => tdep_calc_phi3ref
154 : : procedure :: calc_phi4ref => tdep_calc_phi4ref
155 : :
156 : : end type tdep_Solver_type
157 : : !!***
158 : :
159 : : contains
160 : :
161 : : !=====================================================================================================
162 : :
163 : : !!****f* ABINIT/m_tdep_solver/tdep_solver_init
164 : : !! NAME
165 : : !! tdep_solver_init
166 : : !!
167 : : !! FUNCTION
168 : : !! Initialize the Solver object: set the dimensions, compute the number of
169 : : !! IFC coefficients at each order (from the shells) and the maximum number
170 : : !! of constraints at each order, and allocate (and set to zero) the arrays.
171 : : !! Note that the constraints at the 4th order are currently disabled.
172 : : !!
173 : : !! INPUTS
174 : : !! Invar = Input object containing the input variables.
175 : : !! Shell1at = Shell object for first-order interactions.
176 : : !! Shell2at = Shell object for second-order interactions.
177 : : !! Shell3at = Shell object for third-order interactions (only used if order >= 3).
178 : : !! Shell4at = Shell object for fourth-order interactions (only used if order == 4).
179 : : !!
180 : : !! OUTPUT
181 : : !!
182 : : !! SIDE EFFECTS
183 : : !! Solver = TDEP Solver object, initialized.
184 : : !!
185 : : !! SOURCE
186 : :
187 : 46 : subroutine tdep_solver_init(Solver, Invar, Shell1at, Shell2at, Shell3at, Shell4at)
188 : :
189 : : class(tdep_Solver_type), intent(inout) :: Solver
190 : : type(atdep_dataset_type),intent(in) :: Invar
191 : : type(Shell_type),intent(in) :: Shell1at, Shell2at, Shell3at, Shell4at
192 : :
193 : 46 : Solver%order = Invar%order
194 : 46 : Solver%natom = Invar%natom
195 : 46 : Solver%natom_unitcell = Invar%natom_unitcell
196 : 46 : Solver%my_nstep = Invar%my_nstep
197 : :
198 : : !Rotational invariances (1st order)
199 : : ! constraints = 3
200 : 46 : Solver%nconst_1st = 3**2
201 : :
202 : : !Rotational invariances (2nd order) + Symetry of the Dynamical Matrix + Huang invariances
203 : : ! constraints = natom*3**2 + (3*natom_unitcell)**2 + 3**4
204 : 46 : Solver%nconst_rot2nd = 3**3*Solver%natom_unitcell
205 : 46 : Solver%nconst_dynmat = (3*Solver%natom_unitcell)**2
206 : 46 : Solver%nconst_huang = 3**4
207 : 46 : Solver%nconst_2nd = Solver%nconst_rot2nd + Solver%nconst_dynmat + Solver%nconst_huang
208 : :
209 : : !Rotational invariances (3rd order) + acoustic sum rules (3rd order)
210 : 46 : Solver%nconst_3rd=0
211 [ + + ]: 46 : if (Solver%order.ge.3) then
212 : : ! constraints = natom_unitcell*natom*3**4 + natom_unitcell*natom*3**3
213 : 16 : Solver%nconst_rot3rd = 3**4 * Solver%natom_unitcell * Solver%natom
214 : 16 : Solver%nconst_asr3rd = 3**3 * Solver%natom_unitcell * Solver%natom
215 : 16 : Solver%nconst_3rd = Solver%nconst_rot3rd + Solver%nconst_asr3rd
216 : : end if
217 : :
218 : : !Rotational invariances (4th order) + acoustic sum rules (4th order)
219 : 46 : Solver%nconst_4th=0
220 [ + + ]: 46 : if (Solver%order.ge.4) then
221 : : ! constraints = natom_unitcell*natom**2*3**5 + natom_unitcell*natom**2*3**4
222 : : !FB Solver%nconst_rot4th = 3**5*natom_unitcell*natom**2
223 : : !FB4TH Solver%nconst_asr4th = 3**4*natom_unitcell*natom**2
224 : 9 : Solver%nconst_rot4th = 0
225 : 9 : Solver%nconst_asr4th = 0
226 : : Solver%nconst_4th = Solver%nconst_rot4th + Solver%nconst_asr4th
227 : : end if
228 : :
229 : 46 : Solver%ncoeff1st = Shell1at%ntotcoeff
230 : 46 : Solver%ncoeff2nd = Shell2at%ntotcoeff
231 : 46 : Solver%ncoeff3rd = 0
232 : 46 : Solver%ncoeff4th = 0
233 [ + + ]: 46 : if (Solver%order.ge.3) Solver%ncoeff3rd = Shell3at%ntotcoeff
234 [ + + ]: 46 : if (Solver%order.ge.4) Solver%ncoeff4th = Shell4at%ntotcoeff
235 : :
236 : 46 : Solver%ntotcoeff = Solver%ncoeff1st + Solver%ncoeff2nd + Solver%ncoeff3rd + Solver%ncoeff4th
237 : 46 : Solver%ntotconst = Solver%nconst_1st + Solver%nconst_2nd + Solver%nconst_3rd + Solver%nconst_4th
238 : :
239 [ + - - + : 27857931 : ABI_CALLOC(Solver%fcoeff, (3*Solver%natom*Solver%my_nstep,Solver%ntotcoeff))
+ - - + +
- + - + -
- + - + -
+ + + +
+ ]
240 [ - + + - : 56727645 : ABI_CALLOC(Solver%const, (Solver%ntotconst,Solver%ntotcoeff))
- + + - +
- + - - +
- + - + +
+ + + ]
241 [ - + + - : 441354 : ABI_CALLOC(Solver%Forces, (3*Solver%natom*Solver%my_nstep))
+ - + - -
+ - + - +
+ + ]
242 [ - + + - : 3029 : ABI_CALLOC(Solver%theta, (Solver%ntotcoeff))
+ - + - -
+ - + - +
+ + ]
243 : :
244 : 46 : end subroutine tdep_solver_init
245 : : !!***
246 : :
247 : : !=====================================================================================================
248 : :
249 : : !!****f* ABINIT/m_tdep_solver/tdep_solver_free
250 : : !! NAME
251 : : !! tdep_solver_free
252 : : !!
253 : : !! FUNCTION
254 : : !! Free the memory allocated in the Solver object.
255 : : !!
256 : : !! SIDE EFFECTS
257 : : !! Solver = TDEP Solver object. All its allocated arrays are deallocated.
258 : : !!
259 : : !! SOURCE
260 : :
261 : 46 : subroutine tdep_solver_free(Solver)
262 : :
263 : : class(tdep_Solver_type), intent(inout) :: Solver
264 : :
265 [ - + ]: 46 : ABI_FREE(Solver%fcoeff)
266 [ - + ]: 46 : ABI_FREE(Solver%const)
267 [ - + ]: 46 : ABI_FREE(Solver%Forces)
268 [ - + ]: 46 : ABI_FREE(Solver%theta)
269 : :
270 : 46 : end subroutine tdep_solver_free
271 : : !!***
272 : :
273 : : !=====================================================================================================
274 : :
275 : : !!****f* ABINIT/m_tdep_solver/tdep_solver_set_residual_forces
276 : : !! NAME
277 : : !! tdep_solver_set_residual_forces
278 : : !!
279 : : !! FUNCTION
280 : : !! Compute the residual forces into the solver, that is, the differences
281 : : !! between the MD forces and the ones from the TDEP model.
282 : : !! The model forces may be zero, but they are non-zero when we remove
283 : : !! the long-range part of the IFC forces, or when we solve each IFC order
284 : : !! successively, and we remove the forces from the previously computed orders.
285 : : !!
286 : : !! INPUTS
287 : : !! MD = TDEP Sampling object containing the input positions and forces.
288 : : !! Model = TDEP Model object containing the IFC and corresponding forces.
289 : : !!
290 : : !! OUTPUT
291 : : !!
292 : : !! SIDE EFFECTS
293 : : !! Solver = TDEP Solver object that will compute the coefficients.
294 : : !! On exit, Solver%Forces contains the residual forces,
295 : : !! weighted by the weight of each configuration (MD%weights).
296 : : !!
297 : : !! SOURCE
298 : :
299 : 49 : subroutine tdep_solver_set_residual_forces(Solver, MD, Model)
300 : :
301 : : class(tdep_Solver_type), intent(inout) :: Solver
302 : : type(tdep_Sampling_type), intent(in) :: MD
303 : : type(tdep_Model_type), intent(in) :: Model
304 : :
305 : : integer :: ii,jj,istep,iatom
306 : :
307 [ + + ]: 1101 : do istep=1,MD%my_nstep
308 [ + + ]: 159693 : do iatom=1,MD%natom
309 [ + + ]: 635420 : do ii=1,3
310 : 475776 : jj = MD%get_flat_index(ii,iatom,istep)
311 : 634368 : Solver%Forces(jj) = (MD%Forces(jj) - Model%Forces(jj)) * MD%weights(istep)
312 : : end do
313 : : end do
314 : : end do
315 : :
316 : 49 : end subroutine tdep_solver_set_residual_forces
317 : : !!***
318 : :
319 : : !====================================================================================================
320 : :
321 : : !!****f* ABINIT/m_tdep_solver/tdep_calc_phi1fcoeff
322 : : !! NAME
323 : : !! tdep_calc_phi1fcoeff
324 : : !!
325 : : !! FUNCTION
326 : : !! Fill the columns of the matrix Solver%fcoeff corresponding to the
327 : : !! first-order IFC coefficients, i.e. the derivatives of the forces with
328 : : !! respect to these coefficients, for each MD step. The acoustic sum rule
329 : : !! is imposed by subtracting the contribution of each atom from the first
330 : : !! atom of its unitcell.
331 : : !! INPUTS
332 : : !! Shell1at = Shell object for first-order interactions.
333 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
334 : : !!
335 : : !! OUTPUT
336 : : !!
337 : : !! SIDE EFFECTS
338 : : !! Solver = TDEP Solver object. The columns of Solver%fcoeff corresponding
339 : : !! to the first-order coefficients are filled.
340 : : !!
341 : : !! SOURCE
342 : :
343 : 45 : subroutine tdep_calc_phi1fcoeff(Solver,Shell1at,Sym)
344 : :
345 : : class(tdep_Solver_type), intent(inout) :: Solver
346 : : type(Shell_type),intent(in) :: Shell1at
347 : : type(Symmetries_type),intent(in) :: Sym
348 : :
349 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,iatshell,iat_mod
350 : : integer :: icoeff,isym,mu,iatref
351 : : double precision :: terme
352 : :
353 : 45 : write(ab_out,*) ' '
354 : 45 : write(ab_out,*) '#############################################################################'
355 : 45 : write(ab_out,*) '############## Fill the matrices used in the pseudo-inverse #################'
356 : 45 : write(ab_out,*) '#############################################################################'
357 : :
358 : 45 : write(ab_out,*) ' Compute the coefficients (at the 1st order) for Moore-Penrose inversion.'
359 [ + + ]: 109 : do ishell=1,Shell1at%nshell
360 [ - + ]: 64 : if (Shell1at%neighbours(1,ishell)%n_interactions.eq.0) cycle
361 [ - - ]: 45 : do iatshell=1,Shell1at%neighbours(1,ishell)%n_interactions
362 : 0 : iatom=Shell1at%neighbours(1,ishell)%atomj_in_shell(iatshell)
363 : 0 : iat_mod=mod(iatom+Solver%natom_unitcell-1,Solver%natom_unitcell)+1
364 [ # # ]: 0 : if (iat_mod==1) cycle
365 : 0 : iatref=Shell1at%iatref(ishell)
366 : 0 : isym=Shell1at%neighbours(1,ishell)%sym_in_shell(iatshell)
367 : 0 : ncoeff =Shell1at%ncoeff(ishell)
368 : 0 : ncoeff_prev=Shell1at%ncoeff_prev(ishell)
369 [ - - ]: 64 : do mu=1,3
370 [ # # ]: 0 : do icoeff=1,ncoeff
371 [ # # ]: 0 : terme=sum(Sym%S_ref(mu,:,isym,1)*Shell1at%proj(:,icoeff,ishell))
372 [ # # ]: 0 : do istep=1,Solver%my_nstep
373 : : Solver%fcoeff(mu+3*(iatom-1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev)= &
374 : 0 : & Solver%fcoeff(mu+3*(iatom-1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev)+terme
375 : : ! Add all the other contributions, when iat_mod==1 (due to ASR)
376 : : Solver%fcoeff(mu+3*(iatom-iat_mod+1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev)= &
377 : 0 : & Solver%fcoeff(mu+3*(iatom-iat_mod+1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev)-terme
378 : : end do !istep
379 : : end do
380 : : end do
381 : : end do !iatshell
382 : : end do !ishell
383 : :
384 : 45 : end subroutine tdep_calc_phi1fcoeff
385 : : !!***
386 : :
387 : : !====================================================================================================
388 : :
389 : : !!****f* ABINIT/m_tdep_solver/tdep_calc_phi2fcoeff
390 : : !! NAME
391 : : !! tdep_calc_phi2fcoeff
392 : : !!
393 : : !! FUNCTION
394 : : !! Fill the columns of the matrix Solver%fcoeff corresponding to the
395 : : !! second-order IFC coefficients, i.e. the derivatives of the forces with
396 : : !! respect to these coefficients, for each MD step:
397 : : !! F_i^mu(t) = - sum_{j,nu} sum_{alpha,beta} S^{mu alpha} S^{nu beta}
398 : : !! Phi2_{ij}^{alpha beta} (u_j^nu(t) - u_i^nu(t))
399 : : !! where the second-order IFC are expanded on the independent coefficients
400 : : !! of each shell (Shell2at%proj). The acoustic sum rule is imposed by using
401 : : !! the difference of the displacements u_j - u_i.
402 : : !! INPUTS
403 : : !! Shell2at = Shell object for second-order interactions.
404 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
405 : : !! MD = TDEP Sampling object containing the displacements (MD%ucart) and
406 : : !! the weights of the configurations.
407 : : !!
408 : : !! OUTPUT
409 : : !!
410 : : !! SIDE EFFECTS
411 : : !! Solver = TDEP Solver object. The columns of Solver%fcoeff corresponding
412 : : !! to the second-order coefficients are filled.
413 : : !!
414 : : !! SOURCE
415 : :
416 : 45 : subroutine tdep_calc_phi2fcoeff(Solver,Shell2at,Sym,MD)
417 : :
418 : : class(tdep_Solver_type), intent(inout) :: Solver
419 : : type(Shell_type),intent(in) :: Shell2at
420 : : type(Symmetries_type),intent(in) :: Sym
421 : : type(tdep_Sampling_type), intent(in) :: MD
422 : :
423 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,jatom,iatshell
424 : : integer :: icoeff,isym
425 : : integer :: mu,nu,alpha,beta,itrans
426 : : double precision :: terme,temp
427 : : double precision :: udiff(3),SSu(3,9)
428 : 45 : double precision, allocatable :: SS_ref(:,:,:,:,:)
429 : :
430 : : ! For each couple of atoms, transform the Phi2 (3x3) ifc matrix using the symetry operation (S)
431 : : ! Note: iatom=1 is excluded in order to take into account the atomic sum rule (see below)
432 [ - + + - : 374350 : ABI_MALLOC(SS_ref,(3,9,3,Sym%nsym,2)); SS_ref(:,:,:,:,:)=zero
+ - + - +
- - + - +
+ + + + +
+ + + +
+ ]
433 [ + + ]: 1715 : do isym=1,Sym%nsym
434 [ + + ]: 6725 : do mu=1,3
435 [ + + ]: 21710 : do alpha=1,3
436 [ + + ]: 65130 : do nu=1,3
437 [ + + ]: 195390 : do beta=1,3
438 : 135270 : temp=Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu,beta,isym,1)
439 : 135270 : SS_ref(mu,beta+(alpha-1)*3,nu,isym,1)=temp
440 : 180360 : SS_ref(mu,alpha+(beta-1)*3,nu,isym,2)=temp
441 : : end do
442 : : end do
443 : : end do
444 : : end do
445 : : end do
446 : :
447 : 45 : write(ab_out,*) ' Compute the coefficients (at the 2nd order) for Moore-Penrose inversion.'
448 : :
449 [ + + ]: 546 : do ishell=1,Shell2at%nshell
450 [ + + ]: 72724 : do iatom=1,Solver%natom
451 [ + + ]: 72178 : if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
452 [ + + ]: 507937 : do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
453 : 446778 : jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
454 [ + + ]: 446778 : if (iatom==jatom) cycle
455 : 440012 : isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
456 : 440012 : itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
457 : 440012 : ncoeff =Shell2at%ncoeff(ishell)
458 : 440012 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
459 : :
460 [ + + ]: 9284622 : do istep=1,Solver%my_nstep
461 : : ! In order to impose the acoustic sum rule we use (u(j)-u(i))==u_j^\nu
462 : 8772432 : udiff(1)=(MD%ucart(1,jatom,istep)-MD%ucart(1,iatom,istep))*MD%weights(istep)
463 : 8772432 : udiff(2)=(MD%ucart(2,jatom,istep)-MD%ucart(2,iatom,istep))*MD%weights(istep)
464 : 8772432 : udiff(3)=(MD%ucart(3,jatom,istep)-MD%ucart(3,iatom,istep))*MD%weights(istep)
465 : :
466 : : ! F_i^\mu(t)=\sum_{\alpha\beta,j,\nu}S^{\mu\alpha}.S^{\nu\beta}.\Phi_{ij}^{\alpha\beta}.u_j^\nu(t)
467 : 8772432 : SSu(:,:)=zero
468 [ + + ]: 35089728 : do nu=1,3
469 [ + + + + ]: 982512384 : SSu(:,:)=SSu(:,:)+SS_ref(:,:,nu,isym,itrans)*udiff(nu)
470 : : end do
471 [ + + ]: 35536506 : do mu=1,3
472 [ + + ]: 149067456 : do icoeff=1,ncoeff
473 [ + + ]: 1139777280 : terme=sum(SSu(mu,:)*Shell2at%proj(:,icoeff,ishell))
474 : : !FB write(std_out,*) 'indices=', mu+3*(iatom-1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev
475 : : Solver%fcoeff(mu+3*(iatom-1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev)= &
476 : 140295024 : & Solver%fcoeff(mu+3*(iatom-1)+3*Solver%natom*(istep-1),icoeff+ncoeff_prev)+terme
477 : : end do
478 : : end do
479 : :
480 : : end do !istep
481 : : end do !iatshell
482 : : end do !iatom
483 : : end do !ishell
484 [ - + ]: 45 : ABI_FREE(SS_ref)
485 : :
486 : 45 : end subroutine tdep_calc_phi2fcoeff
487 : : !!***
488 : :
489 : : !====================================================================================================
490 : :
491 : : !!****f* ABINIT/m_tdep_solver/tdep_calc_phi3fcoeff
492 : : !! NAME
493 : : !! tdep_calc_phi3fcoeff
494 : : !!
495 : : !! FUNCTION
496 : : !! Fill the columns of the matrix Solver%fcoeff corresponding to the
497 : : !! third-order IFC coefficients, i.e. the derivatives of the forces with
498 : : !! respect to these coefficients, for each MD step:
499 : : !! F_i^mu(t) = - 1/2 sum_{j,k,nu,xi} sum_{alpha,beta,gamma}
500 : : !! S^{mu alpha} S^{nu beta} S^{xi gamma}
501 : : !! Phi3_{ijk}^{alpha beta gamma} u_j^nu(t) u_k^xi(t)
502 : : !! where the third-order IFC are expanded on the independent coefficients
503 : : !! of each shell (Shell3at%proj).
504 : : !! INPUTS
505 : : !! Shell3at = Shell object for third-order interactions.
506 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
507 : : !! MD = TDEP Sampling object containing the displacements (MD%ucart) and
508 : : !! the weights of the configurations.
509 : : !!
510 : : !! OUTPUT
511 : : !!
512 : : !! SIDE EFFECTS
513 : : !! Solver = TDEP Solver object. The columns of Solver%fcoeff corresponding
514 : : !! to the third-order coefficients are filled.
515 : : !!
516 : : !! SOURCE
517 : :
518 : 16 : subroutine tdep_calc_phi3fcoeff(Solver,Shell3at,Sym,MD)
519 : :
520 : : class(tdep_Solver_type), intent(inout) :: Solver
521 : : type(Shell_type),intent(in) :: Shell3at
522 : : type(Symmetries_type),intent(in) :: Sym
523 : : type(tdep_Sampling_type), intent(in) :: MD
524 : :
525 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,jatom,katom
526 : : integer :: icoeff,isym,itrans,iatshell
527 : : integer :: mu,nu,xi,alpha,beta,gama,iindex
528 : : double precision :: temp
529 : : double precision :: udiff_ki(3),udiff_ji(3)
530 : 16 : double precision, allocatable :: SSS_proj(:,:,:,:)
531 : : double precision :: SSS_tmp(27), proj_tmp(27)
532 : 16 : type(Constraints_type) :: Const
533 : :
534 [ - + + - : 3552 : ABI_MALLOC(Const%Sprod,(Sym%nsym,6))
+ - + - -
+ - + + +
+ + ]
535 [ + + ]: 584 : do isym=1,Sym%nsym
536 [ + + ]: 3992 : do itrans=1,6
537 [ - + - + : 3357448 : ABI_MALLOC(Const%Sprod(isym,itrans)%SSS,(3,27,3,3)); Const%Sprod(isym,itrans)%SSS(:,:,:,:)=zero
+ + + + +
+ + + ]
538 : : end do
539 : : end do
540 : :
541 : : ! For each couple of atoms, transform the Phi3 (3x3x3) ifc matrix using the symetry operation (S)
542 : : ! Note: iatom=1 is excluded in order to take into account the atomic sum rule (see below)
543 [ + + ]: 584 : do isym=1,Sym%nsym
544 [ + + ]: 2288 : do mu=1,3
545 [ + + ]: 7384 : do alpha=1,3
546 [ + + ]: 22152 : do nu=1,3
547 [ + + ]: 66456 : do beta=1,3
548 [ + + ]: 199368 : do xi=1,3
549 [ + + ]: 598104 : do gama=1,3
550 : 414072 : temp=Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu,beta,isym,1)*Sym%S_ref(xi,gama,isym,1)
551 : 414072 : Const%Sprod(isym,1)%SSS(mu,gama+(beta-1)*3+(alpha-1)*9,nu,xi)=temp !\Phi3_efg
552 : 414072 : Const%Sprod(isym,2)%SSS(mu,gama+(beta-1)*3+(alpha-1)*9,xi,nu)=temp !\Phi3_egf
553 : 414072 : Const%Sprod(isym,3)%SSS(nu,gama+(beta-1)*3+(alpha-1)*9,mu,xi)=temp !\Phi3_feg
554 : 414072 : Const%Sprod(isym,4)%SSS(nu,gama+(beta-1)*3+(alpha-1)*9,xi,mu)=temp !\Phi3_fge
555 : 414072 : Const%Sprod(isym,5)%SSS(xi,gama+(beta-1)*3+(alpha-1)*9,mu,nu)=temp !\Phi3_gef
556 : 552096 : Const%Sprod(isym,6)%SSS(xi,gama+(beta-1)*3+(alpha-1)*9,nu,mu)=temp !\Phi3_gfe
557 : : end do
558 : : end do
559 : : end do
560 : : end do
561 : : end do
562 : : end do
563 : : end do
564 : :
565 : 16 : write(ab_out,*) ' Compute the coefficients (at the 3rd order) used in the Moore-Penrose...'
566 [ + + ]: 73 : do ishell=1,Shell3at%nshell
567 [ + + ]: 9913 : do iatom=1,Solver%natom
568 [ + + ]: 9840 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
569 [ + + ]: 106113 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
570 : 97496 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
571 : 97496 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
572 : : !FB if (iatom==jatom.or.iatom==katom) cycle
573 : 97496 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
574 : 97496 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
575 : 97496 : ncoeff =Shell3at%ncoeff(ishell)
576 : 97496 : ncoeff_prev=Shell3at%ncoeff_prev(ishell)+Solver%ncoeff2nd+Solver%ncoeff1st
577 : :
578 [ + + + - : 28137288 : ABI_MALLOC(SSS_proj,(3,3,3,ncoeff)) ; SSS_proj(:,:,:,:)=zero
+ - + + -
+ - + + +
+ + + + +
+ ]
579 [ + + ]: 389984 : do mu=1,3
580 [ + + ]: 1267448 : do nu=1,3
581 [ + + ]: 3802344 : do xi=1,3
582 [ + + ]: 73706976 : SSS_tmp(:)=Const%Sprod(isym,itrans)%SSS(mu,:,nu,xi)
583 [ + + ]: 22305096 : do icoeff=1,ncoeff
584 [ + + ]: 526266720 : proj_tmp(:)=Shell3at%proj(:,icoeff,ishell)
585 : 21427632 : SSS_proj(mu,nu,xi,icoeff)=DDOT(27,SSS_tmp,1,proj_tmp,1)
586 : : end do
587 : : end do
588 : : end do
589 : : end do
590 [ + + ]: 2047416 : do istep=1,Solver%my_nstep
591 : 1949920 : iindex=3*(iatom-1)+3*Solver%natom*(istep-1)
592 : : ! In order to impose the acoustic sum rule we use :
593 : : !FB udiff_ji(:)=MD%ucart(:,jatom,istep)-MD%ucart(:,iatom,istep)
594 : : !FB udiff_ki(:)=MD%ucart(:,katom,istep)-MD%ucart(:,iatom,istep)
595 [ + + ]: 7799680 : udiff_ji(:)=MD%ucart(:,jatom,istep)
596 [ + + ]: 7799680 : udiff_ki(:)=MD%ucart(:,katom,istep)
597 : : ! F_i^{\mu}(t)=\sum_{\alpha\beta\gamma,jk,\nu\xi} S^{\mu\alpha}.S^{\nu\beta}.S^{\xi\gamma}.
598 : : ! \Phi3_{ijk}^{\alpha\beta\gamma}.udiff_k^\xi(t).udiff_j^\nu(t)
599 [ + + ]: 7897176 : do nu=1,3
600 [ + + ]: 25348960 : do xi=1,3
601 : : Solver%fcoeff(iindex+1:iindex+3,ncoeff_prev+1:ncoeff_prev+ncoeff)= &
602 : : & Solver%fcoeff(iindex+1:iindex+3,ncoeff_prev+1:ncoeff_prev+ncoeff)+&
603 [ + + + + ]: 524605440 : & SSS_proj(1:3,nu,xi,1:ncoeff)*udiff_ji(nu)*udiff_ki(xi)/2.d0*MD%weights(istep)
604 : : end do
605 : : end do
606 : : end do !istep
607 [ - + ]: 107336 : ABI_FREE(SSS_proj)
608 : : end do !iatshell
609 : : end do !iatom
610 : : end do !ishell
611 [ + + ]: 584 : do isym=1,Sym%nsym
612 [ + + ]: 3992 : do itrans=1,6
613 [ - + ]: 3976 : ABI_FREE(Const%Sprod(isym,itrans)%SSS)
614 : : end do
615 : : end do
616 [ + - + + : 3440 : ABI_FREE(Const%Sprod)
- + - + -
+ - + ]
617 : :
618 : 16 : end subroutine tdep_calc_phi3fcoeff
619 : : !!***
620 : :
621 : : !====================================================================================================
622 : :
623 : : !!****f* ABINIT/m_tdep_solver/tdep_calc_phi4fcoeff
624 : : !! NAME
625 : : !! tdep_calc_phi4fcoeff
626 : : !!
627 : : !! FUNCTION
628 : : !! Fill the columns of the matrix Solver%fcoeff corresponding to the
629 : : !! fourth-order IFC coefficients, i.e. the derivatives of the forces with
630 : : !! respect to these coefficients, for each MD step:
631 : : !! F_i^mu(t) = - 1/6 sum_{j,k,l,nu,xi,zeta} sum_{alpha,beta,gamma,delta}
632 : : !! S^{mu alpha} S^{nu beta} S^{xi gamma} S^{zeta delta}
633 : : !! Phi4_{ijkl}^{alpha beta gamma delta} u_j^nu(t) u_k^xi(t) u_l^zeta(t)
634 : : !! where the fourth-order IFC are expanded on the independent coefficients
635 : : !! of each shell (Shell4at%proj).
636 : : !! INPUTS
637 : : !! Shell4at = Shell object for fourth-order interactions.
638 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
639 : : !! MD = TDEP Sampling object containing the displacements (MD%ucart) and
640 : : !! the weights of the configurations.
641 : : !!
642 : : !! OUTPUT
643 : : !!
644 : : !! SIDE EFFECTS
645 : : !! Solver = TDEP Solver object. The columns of Solver%fcoeff corresponding
646 : : !! to the fourth-order coefficients are filled.
647 : : !!
648 : : !! SOURCE
649 : :
650 : 9 : subroutine tdep_calc_phi4fcoeff(Solver,Shell4at,Sym,MD)
651 : :
652 : : class(tdep_Solver_type), intent(inout) :: Solver
653 : : type(Shell_type),intent(in) :: Shell4at
654 : : type(Symmetries_type),intent(in) :: Sym
655 : : type(tdep_Sampling_type), intent(in) :: MD
656 : :
657 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,jatom,katom,latom
658 : : integer :: icoeff,isym,iatshell,itrans,counter
659 : : integer :: mu,nu,xi,zeta,alpha,beta,gama,delta,iindex_l,iindex_h
660 : : integer :: ncoeff_prev_l,ncoeff_prev_h
661 : : double precision :: temp,SSSS_tmp(81),proj_tmp(81)
662 : 9 : double precision, allocatable :: SSSS_proj(:,:,:,:,:)
663 : 9 : type(Constraints_type) :: Const
664 : :
665 [ - + + - : 8211 : ABI_MALLOC(Const%Sprod,(Sym%nsym,24))
+ - + - -
+ - + + +
+ + ]
666 [ + + ]: 341 : do isym=1,Sym%nsym
667 [ + + ]: 8309 : do itrans=1,24
668 [ - + - + : 70023116 : ABI_MALLOC(Const%Sprod(isym,itrans)%SSSS,(3,81,3,3,3)); Const%Sprod(isym,itrans)%SSSS(:,:,:,:,:)=zero
+ + + + +
+ + + +
+ ]
669 : : end do
670 : : end do
671 : :
672 : : ! For each couple of atoms, transform the Phi4 (3x3x3) ifc matrix using the symetry operation (S)
673 : : ! Note: iatom=1 is excluded in order to take into account the atomic sum rule (see below)
674 [ + + ]: 341 : do isym=1,Sym%nsym
675 [ + + ]: 1337 : do mu=1,3
676 [ + + ]: 4316 : do alpha=1,3
677 [ + + ]: 12948 : do nu=1,3
678 [ + + ]: 38844 : do beta=1,3
679 [ + + ]: 116532 : do xi=1,3
680 [ + + ]: 349596 : do gama=1,3
681 [ + + ]: 1048788 : do zeta=1,3
682 [ + + ]: 3146364 : do delta=1,3
683 : : #if defined FC_NVHPC
684 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
685 : : #endif
686 : :
687 : 2178252 : counter=delta+(gama-1)*3+(beta-1)*9+(alpha-1)*27
688 : : temp=Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu ,beta ,isym,1)*&
689 : 2178252 : & Sym%S_ref(xi,gama ,isym,1)*Sym%S_ref(zeta,delta,isym,1)
690 : 2178252 : Const%Sprod(isym,1 )%SSSS(mu,counter,nu,xi,zeta)=temp !\Phi4_efgh
691 : 2178252 : Const%Sprod(isym,2 )%SSSS(mu,counter,xi,nu,zeta)=temp !\Phi4_egfh
692 : 2178252 : Const%Sprod(isym,3 )%SSSS(nu,counter,mu,xi,zeta)=temp !\Phi4_fegh
693 : 2178252 : Const%Sprod(isym,4 )%SSSS(nu,counter,xi,mu,zeta)=temp !\Phi4_fgeh
694 : 2178252 : Const%Sprod(isym,5 )%SSSS(xi,counter,mu,nu,zeta)=temp !\Phi4_gefh
695 : 2178252 : Const%Sprod(isym,6 )%SSSS(xi,counter,nu,mu,zeta)=temp !\Phi4_gfeh
696 : :
697 : 2178252 : Const%Sprod(isym,7 )%SSSS(mu,counter,nu,zeta,xi)=temp !\Phi4_efhg
698 : 2178252 : Const%Sprod(isym,8 )%SSSS(mu,counter,xi,zeta,nu)=temp !\Phi4_eghf
699 : 2178252 : Const%Sprod(isym,9 )%SSSS(nu,counter,mu,zeta,xi)=temp !\Phi4_fehg
700 : 2178252 : Const%Sprod(isym,10)%SSSS(nu,counter,xi,zeta,mu)=temp !\Phi4_fghe
701 : 2178252 : Const%Sprod(isym,11)%SSSS(xi,counter,mu,zeta,nu)=temp !\Phi4_gehf
702 : 2178252 : Const%Sprod(isym,12)%SSSS(xi,counter,nu,zeta,mu)=temp !\Phi4_gfhe
703 : :
704 : 2178252 : Const%Sprod(isym,13)%SSSS(mu,counter,zeta,nu,xi)=temp !\Phi4_ehfg
705 : 2178252 : Const%Sprod(isym,14)%SSSS(mu,counter,zeta,xi,nu)=temp !\Phi4_ehgf
706 : 2178252 : Const%Sprod(isym,15)%SSSS(nu,counter,zeta,mu,xi)=temp !\Phi4_fheg
707 : 2178252 : Const%Sprod(isym,16)%SSSS(nu,counter,zeta,xi,mu)=temp !\Phi4_fhge
708 : 2178252 : Const%Sprod(isym,17)%SSSS(xi,counter,zeta,mu,nu)=temp !\Phi4_ghef
709 : 2178252 : Const%Sprod(isym,18)%SSSS(xi,counter,zeta,nu,mu)=temp !\Phi4_ghfe
710 : :
711 : 2178252 : Const%Sprod(isym,19)%SSSS(zeta,counter,mu,nu,xi)=temp !\Phi4_hefg
712 : 2178252 : Const%Sprod(isym,20)%SSSS(zeta,counter,mu,xi,nu)=temp !\Phi4_hegf
713 : 2178252 : Const%Sprod(isym,21)%SSSS(zeta,counter,nu,mu,xi)=temp !\Phi4_hfeg
714 : 2178252 : Const%Sprod(isym,22)%SSSS(zeta,counter,nu,xi,mu)=temp !\Phi4_hfge
715 : 2178252 : Const%Sprod(isym,23)%SSSS(zeta,counter,xi,mu,nu)=temp !\Phi4_hgef
716 : 2904336 : Const%Sprod(isym,24)%SSSS(zeta,counter,xi,nu,mu)=temp !\Phi4_hgfe
717 : :
718 : : end do
719 : : end do
720 : : end do
721 : : end do
722 : : end do
723 : : end do
724 : : end do
725 : : end do
726 : : end do
727 : :
728 : 9 : write(ab_out,*) ' Compute the coefficients (at the 4th order) used in the Moore-Penrose...'
729 [ + + ]: 58 : do ishell=1,Shell4at%nshell
730 [ + + ]: 8538 : do iatom=1,Solver%natom
731 [ + + ]: 8480 : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
732 [ + + ]: 257501 : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
733 : 249996 : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
734 : 249996 : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
735 : 249996 : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
736 : 249996 : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
737 : 249996 : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
738 : 249996 : ncoeff =Shell4at%ncoeff(ishell)
739 : 249996 : ncoeff_prev=Shell4at%ncoeff_prev(ishell)+Solver%ncoeff3rd+Solver%ncoeff2nd+Solver%ncoeff1st
740 : 249996 : ncoeff_prev_l=ncoeff_prev+1
741 : 249996 : ncoeff_prev_h=ncoeff_prev+ncoeff
742 : : #if defined FC_NVHPC
743 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
744 : : #endif
745 [ - + + - : 703567840 : ABI_MALLOC(SSSS_proj,(3,3,3,3,ncoeff)) ; SSSS_proj(:,:,:,:,:)=zero
+ - + - -
+ - + + +
+ + + + +
+ + + ]
746 [ + + ]: 999984 : do mu=1,3
747 [ + + ]: 3249948 : do nu=1,3
748 [ + + ]: 9749844 : do xi=1,3
749 [ + + ]: 29249532 : do zeta=1,3
750 [ + + ]: 1660473432 : SSSS_tmp(:)=Const%Sprod(isym,itrans)%SSSS(mu,:,nu,xi,zeta)
751 [ + + ]: 497480940 : do icoeff=1,ncoeff
752 [ + + ]: 38579472504 : proj_tmp(:)=Shell4at%proj(:,icoeff,ishell)
753 : : ! SSSS_proj(mu,nu,xi,zeta,icoeff)=DDOT(81,Const%Sprod(isym,itrans)%SSSS(mu,:,nu,xi,zeta),1,proj(:,icoeff,ishell),1)
754 : 490731048 : SSSS_proj(mu,nu,xi,zeta,icoeff)=DDOT(81,SSSS_tmp,1,proj_tmp,1)
755 : : end do
756 : : end do
757 : : end do
758 : : end do
759 : : end do
760 [ + + ]: 5249916 : do istep=1,Solver%my_nstep
761 : 4999920 : iindex_l=3*(iatom-1)+3*Solver%natom*(istep-1)+1
762 : 4999920 : iindex_h=3*(iatom-1)+3*Solver%natom*(istep-1)+3
763 : : ! F_i^{\mu}(t)=\sum_{\alpha\beta\gamma\delta,jkl,\nu\xi\zeta} S^{\mu\alpha}.S^{\nu\beta}.S^{\xi\gamma}.S^{\zeta\delta}.
764 : : ! \Phi4_{ijkl}^{\alpha\beta\gamma\delta}.u_l^\zeta(t).u_k^\xi(t).u_j^\nu(t)
765 [ + + ]: 20249676 : do nu=1,3
766 [ + + ]: 64998960 : do xi=1,3
767 [ + + ]: 194996880 : do zeta=1,3
768 : : Solver%fcoeff(iindex_l:iindex_h,ncoeff_prev_l:ncoeff_prev_h)= &
769 : : & Solver%fcoeff(iindex_l:iindex_h,ncoeff_prev_l:ncoeff_prev_h)+&
770 : : & SSSS_proj(1:3,nu,xi,zeta,1:ncoeff)*MD%ucart(nu,jatom,istep)*MD%ucart(xi,katom,istep)*MD%ucart(zeta,latom,istep)/6.d0 *&
771 [ + + + + ]: 12726167040 : & MD%weights(istep)
772 : : end do
773 : : end do
774 : : end do
775 : : end do !istep
776 [ - + ]: 258476 : ABI_FREE(SSSS_proj)
777 : : end do !iatshell
778 : : end do !iatom
779 : : end do !ishell
780 [ + + ]: 341 : do isym=1,Sym%nsym
781 [ + + ]: 8309 : do itrans=1,24
782 [ - + ]: 8300 : ABI_FREE(Const%Sprod(isym,itrans)%SSSS)
783 : : end do
784 : : end do
785 [ + - + + : 7986 : ABI_FREE(Const%Sprod)
- + - + -
+ - + ]
786 : :
787 : 9 : end subroutine tdep_calc_phi4fcoeff
788 : : !!***
789 : :
790 : : !=====================================================================================================
791 : :
792 : : !!****f* ABINIT/m_tdep_solver/tdep_calc_MoorePenrose
793 : : !! NAME
794 : : !! tdep_calc_MoorePenrose
795 : : !!
796 : : !! FUNCTION
797 : : !! Solve the constrained least-squares problem giving the IFC coefficients:
798 : : !! minimize |F + fcoeff.theta|^2 under the constraints const.theta = 0,
799 : : !! using Lagrange multipliers. The corresponding linear system
800 : : !! | 2 fcoeff^T fcoeff const^T | | theta | | -2 fcoeff^T F |
801 : : !! | const 0 | | lambda | = | b |
802 : : !! is built (summing over the MD steps held by all the processors), and
803 : : !! solved by inverting its matrix (LU decomposition). Small diagonal
804 : : !! elements are regularized.
805 : : !! The coefficients of all orders can be solved simultaneously, or order by
806 : : !! order. In the latter case, only the coefficients of the current order are
807 : : !! solved for, and b = -const_lower.theta_lower contains the contribution of
808 : : !! the coefficients already computed at the lower orders (b=0 otherwise).
809 : : !!
810 : : !! INPUTS
811 : : !! MPIdata = Info on MPI parallelism.
812 : : !! simult = Selects the coefficients to be solved:
813 : : !! 0 = all orders simultaneously (Invar%together=1)
814 : : !! 1 = 1st and 2nd orders (Invar%together=0)
815 : : !! 2 = 3rd order (Invar%together=0)
816 : : !! 3 = 4th order (Invar%together=0)
817 : : !!
818 : : !! OUTPUT
819 : : !!
820 : : !! SIDE EFFECTS
821 : : !! Solver = TDEP Solver object. On entry, Solver%fcoeff, Solver%Forces and
822 : : !! Solver%const must be computed. On exit, the corresponding
823 : : !! coefficients of Solver%theta are computed.
824 : : !!
825 : : !! SOURCE
826 : :
827 : 47 : subroutine tdep_calc_MoorePenrose(Solver,MPIdata,simult)
828 : :
829 : : class(tdep_Solver_type), intent(inout) :: Solver
830 : : type(MPI_enreg_type), intent(in) :: MPIdata
831 : : integer, intent(in) :: simult
832 : :
833 : : integer :: INFO,ntotcoeff,ntotconst
834 : : integer :: natnstep,nconcoef,ierr,ncoeff_prev,nconst_prev,iconst,icoeff
835 : 47 : integer, allocatable :: IPIV(:)
836 : 47 : double precision, allocatable :: WORK(:)
837 : 47 : double precision, allocatable :: ffcoeff_tmp(:,:),fforces_tmp(:),b_const(:)
838 : 47 : double precision, allocatable :: A_tot(:,:),A_inv(:,:),b_tot(:),x_tot(:)
839 : :
840 : 47 : write(ab_out,*) '################### And compute the pseudo-inverse ##########################'
841 : 47 : write(ab_out,*) '#############################################################################'
842 : :
843 : 47 : natnstep = 3 * Solver%natom * Solver%my_nstep
844 : :
845 [ + + ]: 47 : if (simult.eq.0) then
846 : : ! Simultaneously (Invar%together=1)
847 : 44 : ncoeff_prev=0
848 : 44 : nconst_prev=0
849 : 44 : ntotcoeff=Solver%ntotcoeff
850 : 44 : ntotconst=Solver%ntotconst
851 [ + + ]: 3 : else if (simult.eq.1) then
852 : : ! Successively (Invar%together=0 and Invar%order=2)
853 : 1 : ncoeff_prev=0
854 : 1 : nconst_prev=0
855 : 1 : ntotcoeff=Solver%ncoeff1st +Solver%ncoeff2nd
856 : 1 : ntotconst=Solver%nconst_1st+Solver%nconst_2nd
857 [ + + ]: 2 : else if (simult.eq.2) then
858 : : ! Successively (Invar%together=0 and Invar%order=3)
859 : 1 : ncoeff_prev=Solver%ncoeff1st +Solver%ncoeff2nd
860 : 1 : nconst_prev=Solver%nconst_1st+Solver%nconst_2nd
861 : 1 : ntotcoeff=Solver%ncoeff3rd
862 : 1 : ntotconst=Solver%nconst_3rd
863 [ + - ]: 1 : else if (simult.eq.3) then
864 : : ! Successively (Invar%together=0 and Invar%order=4)
865 : 1 : ncoeff_prev=Solver%ncoeff1st +Solver%ncoeff2nd +Solver%ncoeff3rd
866 : 1 : nconst_prev=Solver%nconst_1st+Solver%nconst_2nd+Solver%nconst_3rd
867 : 1 : ntotcoeff=Solver%ncoeff4th
868 : 1 : ntotconst=Solver%nconst_4th
869 : : end if
870 : 47 : nconcoef=ntotcoeff+ntotconst
871 [ + + + + ]: 47 : if ((ntotconst.gt.0).and.(simult.ge.2)) then
872 [ - + + - : 23 : ABI_MALLOC(b_const,(ntotconst)) ; b_const(:)=0.d0
+ - - + -
+ + + ]
873 [ + + ]: 21 : do iconst=1,ntotconst
874 [ + + ]: 521 : do icoeff=1,ncoeff_prev
875 : : b_const(iconst)=b_const(iconst)+&
876 : 520 : & Solver%const(nconst_prev+iconst,icoeff)*Solver%theta(icoeff)
877 : : end do
878 : : end do
879 : : end if
880 : :
881 [ - + + - : 472990 : ABI_CALLOC(ffcoeff_tmp,(ntotcoeff,ntotcoeff))
- + + - +
- + - - +
- + + + +
+ ]
882 [ - + + - : 2936 : ABI_CALLOC(fforces_tmp,(ntotcoeff))
+ - + - -
+ - + +
+ ]
883 [ - + + - : 676132 : ABI_CALLOC(A_tot,(nconcoef,nconcoef))
- + + - +
- + - - +
- + + + +
+ ]
884 [ - + + - : 676132 : ABI_CALLOC(A_inv,(nconcoef,nconcoef))
- + + - -
+ + - + +
+ + ]
885 [ - + + - : 3407 : ABI_CALLOC(b_tot,(nconcoef))
+ - + - -
+ + - +
+ ]
886 [ - + + - : 3360 : ABI_CALLOC(x_tot,(nconcoef))
- + + - +
+ ]
887 : : call DGEMM('T','N',ntotcoeff,ntotcoeff,natnstep,2.d0,&
888 : : & Solver%fcoeff(:,ncoeff_prev+1:ncoeff_prev+ntotcoeff),natnstep,&
889 : : & Solver%fcoeff(:,ncoeff_prev+1:ncoeff_prev+ntotcoeff),natnstep,&
890 : 47 : & 0.d0,ffcoeff_tmp,ntotcoeff)
891 : : ! NOTE, we have to solve F_ij = -\sum_j \Phi_ij u_j, so we add a minus sign
892 : : call DGEMV('T',natnstep,ntotcoeff,-2.d0,&
893 : : & Solver%fcoeff(:,ncoeff_prev+1:ncoeff_prev+ntotcoeff),natnstep,&
894 : 47 : & Solver%Forces,1,0.d0,fforces_tmp,1)
895 : 47 : call xmpi_sum(ffcoeff_tmp,MPIdata%comm_step,ierr)
896 : 47 : call xmpi_sum(fforces_tmp,MPIdata%comm_step,ierr)
897 : :
898 [ + + + + ]: 472849 : A_tot(1:ntotcoeff,1:ntotcoeff)=ffcoeff_tmp(1:ntotcoeff,1:ntotcoeff)
899 [ - + ]: 47 : ABI_FREE(ffcoeff_tmp)
900 [ + + ]: 47 : if (ntotconst.gt.0) then
901 : : A_tot(ntotcoeff+1:nconcoef,1:ntotcoeff)=&
902 [ + + + + ]: 92802 : & Solver%const(nconst_prev+1:nconst_prev+ntotconst,ncoeff_prev+1:ncoeff_prev+ntotcoeff)
903 : : A_tot(1:ntotcoeff,ntotcoeff+1:nconcoef)=&
904 [ + + + + ]: 90867 : & transpose(Solver%const(nconst_prev+1:nconst_prev+ntotconst,ncoeff_prev+1:ncoeff_prev+ntotcoeff))
905 : : !FB ABI_FREE(Solver%const)
906 : : end if
907 [ + + ]: 2842 : b_tot(1:ntotcoeff)=fforces_tmp(:)
908 [ + + + + ]: 47 : if ((ntotconst.gt.0).and.(simult.ge.2)) then
909 [ + + ]: 21 : b_tot(ntotcoeff+1:nconcoef)=-b_const(1:ntotconst)
910 [ - + ]: 1 : ABI_FREE(b_const)
911 : : end if
912 [ - + ]: 47 : ABI_FREE(fforces_tmp)
913 : :
914 [ + - - + : 16236 : ABI_MALLOC(WORK, (5 * nconcoef)); WORK(:) = 0.d0
+ - + - +
- - + - +
+ + ]
915 [ - + + - : 3407 : ABI_MALLOC(IPIV, (nconcoef)); IPIV(:) = 0
+ - + - -
+ + - +
+ ]
916 [ + + + + ]: 675991 : A_inv(:,:) = A_tot(:,:)
917 : : !BEGIN DEBUG
918 : : !write(ab_out,*) ' '
919 : : !write(ab_out,*) ' The matrix A_inv is (before DGETRF):'
920 : : !do icoeff=1,nconcoef
921 : : ! write(ab_out,*) (A_inv(icoeff,iconst), iconst=1, nconcoef)
922 : : !end do
923 : : !END DEBUG
924 : :
925 : : ! Check for small pivot elements
926 [ + + ]: 3266 : do icoeff=1,nconcoef
927 [ + + ]: 3266 : if (abs(A_inv(icoeff, icoeff)) < tol12) then
928 : 424 : write(std_out,*) ' WARNING: Small pivot value at index ', icoeff, ' : ', A_inv(icoeff, icoeff)
929 : : ! A_inv(icoeff, icoeff) = tol14 ! Regularization to avoid numerical issues
930 [ + + + + : 50895944 : A_inv(icoeff, icoeff) = max(EPSILON(1.0_dp) * maxval(abs(A_inv)), tol12)
+ + + - -
+ - - ]
931 : : end if
932 : : end do
933 : :
934 : : ! Perform LU factorization
935 : 47 : call DGETRF(nconcoef, nconcoef, A_inv, nconcoef, IPIV, INFO)
936 [ - + ]: 47 : if (INFO.ne.0) then
937 : 0 : write(ab_out,*) 'ERROR: Singular matrix detected in DGETRF. INFO=', INFO
938 : 0 : stop
939 : : end if
940 : :
941 : : ! Check for small pivot elements
942 [ + + ]: 3266 : do icoeff=1,nconcoef
943 [ - + ]: 3266 : if (abs(A_inv(icoeff, icoeff)) < tol12) then
944 : 0 : write(std_out,*) ' WARNING: Small pivot value at index ', icoeff, ' : ', A_inv(icoeff, icoeff)
945 : : ! A_inv(icoeff, icoeff) = tol14 ! Regularization to avoid numerical issues
946 [ # # # # : 0 : A_inv(icoeff, icoeff) = max(EPSILON(1.0_dp) * maxval(abs(A_inv)), tol12)
# # # # #
# # # ]
947 : : end if
948 : : end do
949 : :
950 : : ! Compute matrix inverse using LU decomposition
951 : 47 : call DGETRI(nconcoef, A_inv, nconcoef, IPIV, WORK, 5 * nconcoef, INFO)
952 [ - + ]: 47 : if (INFO.ne.0) then
953 : 0 : write(ab_out,*) 'ERROR: Matrix inversion failed in DGETRI. INFO=', INFO
954 : 0 : stop
955 : : end if
956 : : ! BEGIN DEBUG
957 : : !write(ab_out,*) ' '
958 : : !write(ab_out,*) ' The inverse matrix is (after DGETRI):'
959 : : !do icoeff=1,nconcoef
960 : : ! write(ab_out,*) (A_inv(icoeff,iconst), iconst=1, nconcoef)
961 : : !end do
962 : : ! END DEBUG
963 : :
964 [ - + ]: 47 : ABI_FREE(WORK)
965 [ - + ]: 47 : ABI_FREE(IPIV)
966 : :
967 : 47 : call DGEMV('N',nconcoef,nconcoef,1.d0,A_inv,nconcoef,b_tot,1,0.d0,x_tot,1)
968 : 47 : write(ab_out,*) ' Moore-Penrose inversion complete.'
969 : 47 : write(ab_out,*) ' The effective IFC have been computed.'
970 : 47 : write(ab_out,*) ' '
971 : : !BEGIN DEBUG
972 : : !write(ab_out,*) ' The solutions are:'
973 : : !do icoeff=1,nconcoef
974 : : ! write(ab_out,'(1x,i4,1x,f15.10)') icoeff,x_tot(icoeff)
975 : : !end do
976 : : !write(ab_out,'(a,1x,f15.10)')' condition number=',maxval(x_tot(:))/minval(x_tot(:))
977 : : !END DEBUG
978 : :
979 [ + + ]: 2842 : Solver%theta(ncoeff_prev+1:ncoeff_prev+ntotcoeff)=x_tot(1:ntotcoeff)
980 [ - + ]: 47 : ABI_FREE(A_tot)
981 [ - + ]: 47 : ABI_FREE(A_inv)
982 [ - + ]: 47 : ABI_FREE(b_tot)
983 [ - + ]: 47 : ABI_FREE(x_tot)
984 : :
985 [ - + ]: 94 : end subroutine tdep_calc_MoorePenrose
986 : : !!***
987 : :
988 : : !====================================================================================================
989 : :
990 : : !!****f* ABINIT/m_tdep_solver/tdep_calc_constraints
991 : : !! NAME
992 : : !! tdep_calc_constraints
993 : : !!
994 : : !! FUNCTION
995 : : !! Build the constraints on the IFC coefficients:
996 : : !! - 1st order: rotational invariances,
997 : : !! - 2nd order: rotational invariances, symmetry of the dynamical matrix
998 : : !! and Huang invariances,
999 : : !! - 3rd order: acoustic sum rules and rotational invariances,
1000 : : !! - 4th order: acoustic sum rules (currently disabled).
1001 : : !! For each kind of constraint, the set of constraint vectors is
1002 : : !! orthonormalized (see tdep_calc_orthonorm) to keep only the independent
1003 : : !! ones, which are stored in Solver%const. The number of constraints of each
1004 : : !! kind is written in the _constraints.dat file.
1005 : : !!
1006 : : !! INPUTS
1007 : : !! Invar = Input object containing the input variables.
1008 : : !! MPIdata = Info on MPI parallelism.
1009 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
1010 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1011 : : !! Shell1at = Shell object for first-order interactions.
1012 : : !! Shell2at = Shell object for second-order interactions.
1013 : : !! Shell3at = (optional) Shell object for third-order interactions.
1014 : : !! Shell4at = (optional) Shell object for fourth-order interactions.
1015 : : !!
1016 : : !! OUTPUT
1017 : : !!
1018 : : !! SIDE EFFECTS
1019 : : !! Solver = TDEP Solver object. On exit, Solver%const is filled, and the
1020 : : !! numbers of constraints (Solver%nconst_* and Solver%ntotconst)
1021 : : !! are updated to the numbers of independent constraints.
1022 : : !!
1023 : : !! SOURCE
1024 : :
1025 : 46 : subroutine tdep_calc_constraints(Solver,Invar,MPIdata,SC,Sym,&
1026 : : & Shell1at,Shell2at,Shell3at,Shell4at)
1027 : :
1028 : : class(tdep_Solver_type), intent(inout) :: Solver
1029 : : type(atdep_dataset_type),intent(in) :: Invar
1030 : : type(MPI_enreg_type), intent(in) :: MPIdata
1031 : : type(tdep_Supercell_type),intent(in) :: SC
1032 : : type(Symmetries_type),intent(in) :: Sym
1033 : : type(Shell_type),intent(in) :: Shell1at
1034 : : type(Shell_type),intent(in) :: Shell2at
1035 : : type(Shell_type),optional,intent(in) :: Shell3at
1036 : : type(Shell_type),optional,intent(in) :: Shell4at
1037 : :
1038 : : integer :: ishell,ncoeff,ncoeff_prev,iatom,jatom,katom,latom,iatshell,counter
1039 : : integer :: icoeff,iconst,nconst_loc,iconst_loc,iconst_new,isym,itrans,ntotcoeff,iat_mod
1040 : : integer :: mu,nu,xi,zeta,alpha,beta,gama,delta,lambda,natom_unitcell,natom,ii
1041 : : double precision :: terme,temp,terme1,terme2,terme3,terme4
1042 : 46 : double precision, allocatable :: SS_ref(:,:,:,:,:)
1043 : 46 : double precision, allocatable :: vect(:,:)
1044 : 46 : double precision, allocatable :: const_rot1st(:,:,:)
1045 : 46 : double precision, allocatable :: const_rot2nd(:,:,:,:,:)
1046 : 46 : double precision, allocatable :: const_dynmat(:,:,:,:,:)
1047 : 46 : double precision, allocatable :: const_huang(:,:,:,:,:)
1048 : : !FB double precision, allocatable :: const_asr4th(:,:,:,:,:,:)
1049 : : !FB double precision, allocatable :: const_rot4th(:,:,:,:,:,:,:)
1050 : 46 : type(Constraints_type) :: Const3,Const4
1051 : : logical :: order2,order3,order4
1052 : :
1053 : : !TODO Move parts of this routine into m_tdep_constraints
1054 : :
1055 : 46 : natom_unitcell=Solver%natom_unitcell
1056 : 46 : natom =Solver%natom
1057 : :
1058 : 46 : order2 = .false.
1059 : 46 : order3 = .false.
1060 : 46 : order4 = .false.
1061 [ + - ]: 46 : if (Solver%order.ge.2) order2=.true.
1062 [ + + ]: 46 : if (Solver%order.ge.3) order3=.true.
1063 : : !FB4th if (Solver%order.ge.4) order4=.true.
1064 : : if (Solver%order.ge.4) order4=.false.
1065 : :
1066 : 46 : write(ab_out,*) ' '
1067 : 46 : write(ab_out,*) '#############################################################################'
1068 : 46 : write(ab_out,*) '###################### Compute the constraints ##############################'
1069 : :
1070 : : ! For each couple of atoms, transform the Phi2 (3x3) ifc matrix using the symetry operation (S)
1071 [ + - ]: 46 : if (order2.or.order3) then
1072 [ - + + - : 379732 : ABI_CALLOC(SS_ref,(3,9,3,Sym%nsym,2))
+ - + - +
- - + - +
+ + + + +
+ + + +
+ ]
1073 [ + + ]: 1740 : do isym=1,Sym%nsym
1074 [ + + ]: 6822 : do alpha=1,3
1075 [ + + ]: 22022 : do mu=1,3
1076 [ + + ]: 66066 : do beta=1,3
1077 [ + + ]: 198198 : do nu=1,3
1078 : 137214 : temp=Sym%S_ref(alpha,mu,isym,1)*Sym%S_ref(beta,nu,isym,1)
1079 : 137214 : SS_ref(alpha,nu+(mu-1)*3,beta,isym,1)=temp
1080 : 182952 : SS_ref(alpha,mu+(nu-1)*3,beta,isym,2)=temp
1081 : : end do
1082 : : end do
1083 : : end do
1084 : : end do
1085 : : end do
1086 : : end if
1087 : : ! For each couple of atoms, transform the Phi3 (3x3x3) ifc matrix using the symetry operation (S)
1088 [ + + ]: 46 : if (order3.or.order4) then
1089 [ - + + - : 3552 : ABI_MALLOC(Const3%Sprod,(Sym%nsym,6))
+ - + - -
+ - + + +
+ + ]
1090 [ + + ]: 584 : do isym=1,Sym%nsym
1091 [ + + ]: 3992 : do itrans=1,6
1092 [ - + - + : 3357448 : ABI_CALLOC(Const3%Sprod(isym,itrans)%SSS,(3,27,3,3))
+ + + + +
+ + + ]
1093 : : end do
1094 : : end do
1095 [ + + ]: 584 : do isym=1,Sym%nsym
1096 [ + + ]: 2288 : do alpha=1,3
1097 [ + + ]: 7384 : do mu=1,3
1098 [ + + ]: 22152 : do beta=1,3
1099 [ + + ]: 66456 : do nu=1,3
1100 [ + + ]: 199368 : do gama=1,3
1101 [ + + ]: 598104 : do xi=1,3
1102 : 414072 : temp=Sym%S_ref(alpha,mu,isym,1)*Sym%S_ref(beta,nu,isym,1)*Sym%S_ref(gama,xi,isym,1)
1103 : 414072 : Const3%Sprod(isym,1)%SSS(alpha,xi+(nu-1)*3+(mu-1)*9,beta ,gama) =temp !\Phi3_efg
1104 : 414072 : Const3%Sprod(isym,2)%SSS(alpha,xi+(nu-1)*3+(mu-1)*9,gama ,beta) =temp !\Phi3_egf
1105 : 414072 : Const3%Sprod(isym,3)%SSS(beta ,xi+(nu-1)*3+(mu-1)*9,alpha,gama) =temp !\Phi3_feg
1106 : 414072 : Const3%Sprod(isym,4)%SSS(beta ,xi+(nu-1)*3+(mu-1)*9,gama ,alpha)=temp !\Phi3_fge
1107 : 414072 : Const3%Sprod(isym,5)%SSS(gama ,xi+(nu-1)*3+(mu-1)*9,alpha,beta) =temp !\Phi3_gef
1108 : 552096 : Const3%Sprod(isym,6)%SSS(gama ,xi+(nu-1)*3+(mu-1)*9,beta ,alpha)=temp !\Phi3_gfe
1109 : : end do
1110 : : end do
1111 : : end do
1112 : : end do
1113 : : end do
1114 : : end do
1115 : : end do
1116 : : end if
1117 : :
1118 : : ! For each couple of atoms, transform the Phi4 (3x3x3x3) ifc matrix using the symetry operation (S)
1119 : : if (order4) then
1120 : : ABI_MALLOC(Const4%Sprod,(Sym%nsym,24))
1121 : : do isym=1,Sym%nsym
1122 : : do itrans=1,24
1123 : : ABI_CALLOC(Const4%Sprod(isym,itrans)%SSSS,(3,81,3,3,3))
1124 : : end do
1125 : : end do
1126 : : do isym=1,Sym%nsym
1127 : : do alpha=1,3
1128 : : do mu=1,3
1129 : : do beta=1,3
1130 : : do nu=1,3
1131 : : do gama=1,3
1132 : : do xi=1,3
1133 : : do delta=1,3
1134 : : do zeta=1,3
1135 : : counter=zeta+(xi-1)*3+(nu-1)*9+(mu-1)*27
1136 : : temp=Sym%S_ref(alpha,mu,isym,1)*Sym%S_ref(beta ,nu ,isym,1)*&
1137 : : & Sym%S_ref(gama,xi ,isym,1)*Sym%S_ref(delta,zeta,isym,1)
1138 : : Const4%Sprod(isym,1 )%SSSS(alpha,counter,beta,gama,delta)=temp !\Phi4_efgh
1139 : : Const4%Sprod(isym,2 )%SSSS(alpha,counter,gama,beta,delta)=temp !\Phi4_egfh
1140 : : Const4%Sprod(isym,3 )%SSSS(beta,counter,alpha,gama,delta)=temp !\Phi4_fegh
1141 : : Const4%Sprod(isym,4 )%SSSS(beta,counter,gama,alpha,delta)=temp !\Phi4_fgeh
1142 : : Const4%Sprod(isym,5 )%SSSS(gama,counter,alpha,beta,delta)=temp !\Phi4_gefh
1143 : : Const4%Sprod(isym,6 )%SSSS(gama,counter,beta,alpha,delta)=temp !\Phi4_gfeh
1144 : :
1145 : : Const4%Sprod(isym,7 )%SSSS(alpha,counter,beta,delta,gama)=temp !\Phi4_efhg
1146 : : Const4%Sprod(isym,8 )%SSSS(alpha,counter,gama,delta,beta)=temp !\Phi4_eghf
1147 : : Const4%Sprod(isym,9 )%SSSS(beta,counter,alpha,delta,gama)=temp !\Phi4_fehg
1148 : : Const4%Sprod(isym,10)%SSSS(beta,counter,gama,delta,alpha)=temp !\Phi4_fghe
1149 : : Const4%Sprod(isym,11)%SSSS(gama,counter,alpha,delta,beta)=temp !\Phi4_gehf
1150 : : Const4%Sprod(isym,12)%SSSS(gama,counter,beta,delta,alpha)=temp !\Phi4_gfhe
1151 : :
1152 : : Const4%Sprod(isym,13)%SSSS(alpha,counter,delta,beta,gama)=temp !\Phi4_ehfg
1153 : : Const4%Sprod(isym,14)%SSSS(alpha,counter,delta,gama,beta)=temp !\Phi4_ehgf
1154 : : Const4%Sprod(isym,15)%SSSS(beta,counter,delta,alpha,gama)=temp !\Phi4_fheg
1155 : : Const4%Sprod(isym,16)%SSSS(beta,counter,delta,gama,alpha)=temp !\Phi4_fhge
1156 : : Const4%Sprod(isym,17)%SSSS(gama,counter,delta,alpha,beta)=temp !\Phi4_ghef
1157 : : Const4%Sprod(isym,18)%SSSS(gama,counter,delta,beta,alpha)=temp !\Phi4_ghfe
1158 : :
1159 : : Const4%Sprod(isym,19)%SSSS(delta,counter,alpha,beta,gama)=temp !\Phi4_hefg
1160 : : Const4%Sprod(isym,20)%SSSS(delta,counter,alpha,gama,beta)=temp !\Phi4_hegf
1161 : : Const4%Sprod(isym,21)%SSSS(delta,counter,beta,alpha,gama)=temp !\Phi4_hfeg
1162 : : Const4%Sprod(isym,22)%SSSS(delta,counter,beta,gama,alpha)=temp !\Phi4_hfge
1163 : : Const4%Sprod(isym,23)%SSSS(delta,counter,gama,alpha,beta)=temp !\Phi4_hgef
1164 : : Const4%Sprod(isym,24)%SSSS(delta,counter,gama,beta,alpha)=temp !\Phi4_hgfe
1165 : :
1166 : : end do
1167 : : end do
1168 : : end do
1169 : : end do
1170 : : end do
1171 : : end do
1172 : : end do
1173 : : end do
1174 : : end do
1175 : : end if
1176 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1177 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Compute the constraints !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1178 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179 : 46 : ntotcoeff=Solver%ntotcoeff
1180 : : ! First order only
1181 : 46 : write(ab_out,*) '########################## At the 1st order #################################'
1182 [ + - ]: 46 : if (order2) then
1183 [ - + + - : 37721 : ABI_CALLOC(const_rot1st, (3,3,ntotcoeff))
+ - + - -
+ - + + +
+ + + + ]
1184 [ - + + - : 553675 : ABI_CALLOC(const_rot2nd, (3,3,3,natom_unitcell,ntotcoeff))
- + + - +
- + - - +
- + + + +
+ + + + +
+ + ]
1185 [ - + + - : 1864095 : ABI_CALLOC(const_dynmat, (3,3,natom_unitcell,natom_unitcell,ntotcoeff))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + +
+ ]
1186 [ - + + - : 349949 : ABI_CALLOC(const_huang, (3,3,3,3,ntotcoeff))
+ - + - -
+ - + + +
+ + + + +
+ + + ]
1187 [ + + ]: 111 : do ishell=1,Shell1at%nshell
1188 [ - + ]: 65 : if (Shell1at%neighbours(1,ishell)%n_interactions.eq.0) cycle
1189 [ - - ]: 46 : do iatshell=1,Shell1at%neighbours(1,ishell)%n_interactions
1190 : 0 : iatom=Shell1at%neighbours(1,ishell)%atomj_in_shell(iatshell)
1191 [ # # ]: 0 : if (iatom.ge.natom_unitcell) cycle
1192 [ # # ]: 0 : if (iatom.eq.1) cycle
1193 : 0 : isym=Shell1at%neighbours(1,ishell)%sym_in_shell(iatshell)
1194 : 0 : ncoeff =Shell1at%ncoeff(ishell)
1195 : 0 : ncoeff_prev=Shell1at%ncoeff_prev(ishell)
1196 [ - - ]: 65 : do alpha=1,3
1197 [ # # ]: 0 : do beta=1,3
1198 [ # # ]: 0 : do icoeff=1,ncoeff
1199 : : ! 1/ Rotational invariances (1st order)
1200 [ # # ]: 0 : terme1=sum(Sym%S_ref(alpha,:,isym,1)*Shell1at%proj(:,icoeff,ishell))*SC%dij(1,iatom,beta +1)
1201 [ # # ]: 0 : terme2=sum(Sym%S_ref(beta ,:,isym,1)*Shell1at%proj(:,icoeff,ishell))*SC%dij(1,iatom,alpha+1)
1202 : : const_rot1st(alpha,beta,icoeff+ncoeff_prev)= &
1203 : 0 : & const_rot1st(alpha,beta,icoeff+ncoeff_prev)+terme1-terme2
1204 : :
1205 : : ! 2/ Rotational invariances (for the 2nd order)
1206 [ # # ]: 0 : do gama=1,3
1207 : 0 : terme1=zero ; terme2=zero
1208 [ # # # # ]: 0 : if (alpha.eq.gama) terme1=sum(Sym%S_ref(beta,:,isym,1)*Shell1at%proj(:,icoeff,ishell))
1209 [ # # # # ]: 0 : if (alpha.eq.beta) terme2=sum(Sym%S_ref(gama,:,isym,1)*Shell1at%proj(:,icoeff,ishell))
1210 : : const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)=&
1211 : 0 : & const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)+terme1-terme2
1212 : : const_rot2nd(alpha,beta,gama,1,icoeff+ncoeff_prev)=&
1213 : 0 : & const_rot2nd(alpha,beta,gama,1,icoeff+ncoeff_prev)-terme1+terme2
1214 : : end do
1215 : : end do
1216 : : end do
1217 : : end do
1218 : : end do !iatshell
1219 : : end do !ishell
1220 : :
1221 : : ! First + second order
1222 : 46 : write(ab_out,*) '########################## At the 2nd order #################################'
1223 [ + + ]: 568 : do ishell=1,Shell2at%nshell
1224 [ + + ]: 2898 : do iatom=1,natom_unitcell
1225 [ + + ]: 2330 : if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1226 [ + + ]: 7341 : do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
1227 : 5676 : jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1228 [ + + ]: 5676 : if (iatom==jatom) cycle
1229 : 5566 : isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1230 : 5566 : itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1231 : 5566 : ncoeff =Shell2at%ncoeff(ishell)
1232 : 5566 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
1233 : 5566 : iat_mod=mod(jatom+natom_unitcell-1,natom_unitcell)+1
1234 : : ! 1/ Rotational invariances (2nd order). Number of constraints = natom_unitcell*3**2
1235 [ + + ]: 22264 : do alpha=1,3
1236 [ + + ]: 72358 : do beta=1,3
1237 [ + + ]: 217074 : do gama=1,3
1238 [ + + ]: 899136 : do icoeff=1,ncoeff
1239 [ + + ]: 6987600 : terme1=sum(SS_ref(alpha,:,beta,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*SC%dij(iatom,jatom,gama+1)
1240 [ + + ]: 6987600 : terme2=sum(SS_ref(alpha,:,gama,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*SC%dij(iatom,jatom,beta+1)
1241 : : const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)=&
1242 : 849042 : & const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)+terme1-terme2
1243 : : end do
1244 : : end do
1245 : : end do
1246 : : end do
1247 : : ! 2/ Enforce the symetry of the dynamical matrix. Number of constraints = (3*natom_unitcell)**2
1248 : : ! Note that we are unable to enforce the symetry when iatom=jatom (We have to write the equations)
1249 [ + + ]: 22264 : do alpha=1,3
1250 [ + + ]: 72358 : do beta=1,3
1251 [ + + ]: 299712 : do icoeff=1,ncoeff
1252 : : terme=sum(SS_ref(alpha,:,beta,isym,itrans)*Shell2at%proj(:,icoeff,ishell))-&
1253 [ + + + + ]: 4425480 : & sum(SS_ref(beta,:,alpha,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1254 : : const_dynmat(alpha,beta,iatom,iat_mod,icoeff+ncoeff_prev)=&
1255 : 283014 : & const_dynmat(alpha,beta,iatom,iat_mod,icoeff+ncoeff_prev)+terme
1256 : : end do
1257 : : end do
1258 : : end do
1259 : : ! 3/ Huang invariances. Number of constraints = 3**4
1260 [ + + ]: 24594 : do alpha=1,3
1261 [ + + ]: 72468 : do beta=1,3
1262 [ + + ]: 217074 : do gama=1,3
1263 [ + + ]: 651222 : do lambda=1,3
1264 [ + + ]: 2697408 : do icoeff=1,ncoeff
1265 : : terme=sum(SS_ref(alpha,:,beta,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*&
1266 : : & SC%dij(iatom,jatom,gama+1)*&
1267 : : & SC%dij(iatom,jatom,lambda+1)-&
1268 : : & sum(SS_ref(gama,:,lambda,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*&
1269 : : & SC%dij(iatom,jatom,alpha+1)*&
1270 [ + + + + ]: 39829320 : & SC%dij(iatom,jatom,beta+1)
1271 : : const_huang(alpha,beta,gama,lambda,icoeff+ncoeff_prev)=&
1272 : 2547126 : & const_huang(alpha,beta,gama,lambda,icoeff+ncoeff_prev)+terme
1273 : : end do
1274 : : end do
1275 : : end do
1276 : : end do
1277 : : end do
1278 : : end do !iatshell
1279 : : end do !iatom
1280 : : end do !ishell
1281 : : end if !order=1,2
1282 : :
1283 : : ! Third order
1284 [ + + ]: 46 : if (order3) then
1285 : 16 : write(ab_out,*) '########################## At the 3rd order #################################'
1286 [ - + + - : 742520 : ABI_MALLOC(Const3%AsrRot3,(natom_unitcell,natom,ntotcoeff))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + ]
1287 [ + + ]: 48 : do iatom=1,natom_unitcell
1288 [ + + ]: 5728 : do jatom=1,natom
1289 [ + + ]: 513248 : do icoeff=1,ntotcoeff
1290 [ - + - + : 20301440 : ABI_MALLOC(Const3%AsrRot3(iatom,jatom,icoeff)%ABG, (3,3,3)); Const3%AsrRot3(iatom,jatom,icoeff)%ABG(:,:,:) =zero
+ + + + +
+ ]
1291 [ - + - + : 61417536 : ABI_MALLOC(Const3%AsrRot3(iatom,jatom,icoeff)%ABGD,(3,3,3,3)); Const3%AsrRot3(iatom,jatom,icoeff)%ABGD(:,:,:,:)=zero
+ + + + +
+ + + ]
1292 : : end do
1293 : : end do
1294 : : end do
1295 [ + + ]: 189 : do ishell=1,Shell2at%nshell
1296 [ + + ]: 544 : do iatom=1,natom_unitcell
1297 [ + + ]: 355 : if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1298 [ + + ]: 2740 : do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
1299 : 2264 : jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1300 [ + + ]: 2264 : if (iatom==jatom) cycle
1301 : 2232 : isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1302 : 2232 : itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1303 : 2232 : ncoeff =Shell2at%ncoeff(ishell)
1304 : 2232 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
1305 : 2232 : iat_mod=mod(jatom+natom_unitcell-1,natom_unitcell)+1
1306 : : ! 1/ Rotational invariances (coming from the 2nd order). Number of constraints = natom_unitcell*natom*3**3
1307 [ + - ]: 2587 : if (Solver%order.ge.3) then
1308 [ + + ]: 8928 : do alpha=1,3
1309 [ + + ]: 29016 : do beta=1,3
1310 [ + + ]: 87048 : do gama=1,3
1311 [ + + ]: 261144 : do lambda=1,3
1312 [ + + ]: 1018008 : do icoeff=1,ncoeff
1313 : 776952 : terme1=zero ; terme2=zero ; terme3=zero ; terme4=zero ;
1314 [ + + + + ]: 3107808 : if (alpha.eq.lambda) terme1=sum(SS_ref(gama ,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1315 [ + + + + ]: 3107808 : if (beta.eq.lambda) terme2=sum(SS_ref(alpha ,:,gama ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1316 [ + + + + ]: 3107808 : if (alpha.eq.gama) terme3=sum(SS_ref(lambda,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1317 [ + + + + ]: 3107808 : if (beta.eq.gama) terme4=sum(SS_ref(alpha ,:,lambda,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1318 [ + + ]: 776952 : if (SC%dij(iatom,jatom,1).lt.Invar%rcut3) then
1319 : : Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1320 : 43740 : & Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1+terme2-terme3-terme4
1321 : : end if
1322 : : Const3%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1323 : 957744 : & Const3%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)-terme1-terme2+terme3+terme4
1324 : : end do
1325 : : end do
1326 : : end do
1327 : : end do
1328 : : end do
1329 : : end if !proj3rd
1330 : : end do !iatshell
1331 : : end do !iatom
1332 : : end do !ishell
1333 [ + + ]: 73 : do ishell=1,Shell3at%nshell
1334 [ + + ]: 202 : do iatom=1,natom_unitcell
1335 [ + + ]: 129 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1336 [ + + ]: 1506 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
1337 : 1340 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1338 : 1340 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1339 : 1340 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1340 : 1340 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1341 : 1340 : ncoeff =Shell3at%ncoeff(ishell)
1342 : 1340 : ncoeff_prev=Shell3at%ncoeff_prev(ishell)+Solver%ncoeff2nd+Solver%ncoeff1st
1343 : : ! 2/ Acoustic sum rules (3rd order). Number of constraints = natom_unitcell*natom*3**3
1344 [ + + ]: 5360 : do alpha=1,3
1345 [ + + ]: 17420 : do beta=1,3
1346 [ + + ]: 52260 : do gama=1,3
1347 [ + + ]: 297126 : do icoeff=1,ncoeff
1348 [ + + ]: 6968808 : terme =sum(Const3%Sprod(isym,itrans)%SSS(alpha,:,beta,gama)*Shell3at%proj(:,icoeff,ishell))
1349 : : Const3%AsrRot3(iatom,katom,icoeff+ncoeff_prev)%ABG(alpha,beta,gama)=&
1350 : 285066 : & Const3%AsrRot3(iatom,katom,icoeff+ncoeff_prev)%ABG(alpha,beta,gama)+terme
1351 : : end do
1352 : : end do
1353 : : end do
1354 : : end do
1355 : : ! 2/ Rotational invariances (coming from the 3rd order). Number of constraints = natom_unitcell*natom*3**4
1356 [ + + ]: 5489 : do alpha=1,3
1357 [ + + ]: 17420 : do beta=1,3
1358 [ + + ]: 52260 : do gama=1,3
1359 [ + + ]: 156780 : do lambda=1,3
1360 [ + + ]: 891378 : do icoeff=1,ncoeff
1361 : : terme1=sum(Const3%Sprod(isym,itrans)%SSS(alpha,:,beta,gama )&
1362 [ + + ]: 20906424 : & *Shell3at%proj(:,icoeff,ishell))*SC%dij(iatom,katom,lambda+1)
1363 : : terme2=sum(Const3%Sprod(isym,itrans)%SSS(alpha,:,beta,lambda)&
1364 [ + + ]: 20906424 : & *Shell3at%proj(:,icoeff,ishell))*SC%dij(iatom,katom,gama+1)
1365 : : Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1366 : 855198 : & Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1-terme2
1367 : : end do
1368 : : end do
1369 : : end do
1370 : : end do
1371 : : end do
1372 : : end do !iatshell
1373 : : end do !iatom
1374 : : end do !ishell
1375 : : end if !order=3
1376 : :
1377 : : ! Fourth order
1378 : : if (order4) then
1379 : : write(ab_out,*) '########################## At the 4th order #################################'
1380 : : ABI_MALLOC(Const4%AsrRot4,(natom_unitcell,natom,natom,ntotcoeff))
1381 : : do iatom=1,natom_unitcell
1382 : : do jatom=1,natom
1383 : : do katom=1,natom
1384 : : do icoeff=1,ntotcoeff
1385 : : ABI_MALLOC(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGD, (3,3,3,3))
1386 : : Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGD(:,:,:,:) =zero
1387 : : !FB ABI_MALLOC(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGDE,(3,3,3,3,3))
1388 : : !FB Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGDE(:,:,:,:,:)=zero
1389 : : end do
1390 : : end do
1391 : : end do
1392 : : end do
1393 : : !FB do ishell=1,Shell2at%nshell
1394 : : !FB do iatom=1,natom_unitcell
1395 : : !FB if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1396 : : !FB do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
1397 : : !FB jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1398 : : !FB if (iatom==jatom) cycle
1399 : : !FB isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1400 : : !FB itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1401 : : !FB ncoeff =Shell2at%ncoeff(ishell)
1402 : : !FB ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
1403 : : !FB iat_mod=mod(jatom+natom_unitcell-1,natom_unitcell)+1
1404 : : !FB! 1/ Rotational invariances (coming from the 2nd order). Number of constraints = natom_unitcell*natom*3**3
1405 : : !FB if (Solver%order.ge.3) then
1406 : : !FB do alpha=1,3
1407 : : !FB do beta=1,3
1408 : : !FB do gama=1,3
1409 : : !FB do lambda=1,3
1410 : : !FB do icoeff=1,ncoeff
1411 : : !FB terme1=zero ; terme2=zero ; terme3=zero ; terme4=zero ;
1412 : : !FB if (alpha.eq.lambda) terme1=sum(SS_ref(gama ,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1413 : : !FB if (beta.eq.lambda) terme2=sum(SS_ref(alpha ,:,gama ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1414 : : !FB if (alpha.eq.gama) terme3=sum(SS_ref(lambda,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1415 : : !FB if (beta.eq.gama) terme4=sum(SS_ref(alpha ,:,lambda,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1416 : : !FB if (SC%dij(iatom,jatom,1).lt.Invar%rcut3) then
1417 : : !FB Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1418 : : !FB& Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1+terme2-terme3-terme4
1419 : : !FB end if
1420 : : !FB Const4%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1421 : : !FB& Const4%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)-terme1-terme2+terme3+terme4
1422 : : !FB end do
1423 : : !FB end do
1424 : : !FB end do
1425 : : !FB end do
1426 : : !FB end do
1427 : : !FB end if !proj3rd
1428 : : !FB end do !iatshell
1429 : : !FB end do !iatom
1430 : : !FB end do !ishell
1431 : : do ishell=1,Shell4at%nshell
1432 : : do iatom=1,natom_unitcell
1433 : : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1434 : : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1435 : : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1436 : : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1437 : : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
1438 : : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1439 : : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1440 : : ncoeff =Shell4at%ncoeff(ishell)
1441 : : ncoeff_prev=Shell4at%ncoeff_prev(ishell)+Solver%ncoeff3rd+Solver%ncoeff2nd+Solver%ncoeff1st
1442 : : ! 2/ Acoustic sum rules (4th order). Number of constraints = natom_unitcell*natom**2*3**4
1443 : : do alpha=1,3
1444 : : do beta=1,3
1445 : : do gama=1,3
1446 : : do delta=1,3
1447 : : do icoeff=1,ncoeff
1448 : : terme =sum(Const4%Sprod(isym,itrans)%SSSS(alpha,:,beta,gama,delta)*Shell4at%proj(:,icoeff,ishell))
1449 : : Const4%AsrRot4(iatom,katom,latom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,delta)=&
1450 : : & Const4%AsrRot4(iatom,katom,latom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,delta)+terme
1451 : : end do
1452 : : end do
1453 : : end do
1454 : : end do
1455 : : end do
1456 : : !FB! 2/ Rotational invariances (coming from the 3rd order). Number of constraints = natom_unitcell*natom*3**4
1457 : : !FB do alpha=1,3
1458 : : !FB do beta=1,3
1459 : : !FB do gama=1,3
1460 : : !FB do lambda=1,3
1461 : : !FB do icoeff=1,ncoeff
1462 : : !FB terme1=sum(Const4%Sprod(isym,itrans)%SSS(alpha,:,beta,gama )&
1463 : : !FB &*Shell3at%proj(:,icoeff,ishell))*SC%dij(iatom,katom,lambda+1)
1464 : : !FB terme2=sum(Const4%Sprod(isym,itrans)%SSS(alpha,:,beta,lambda)&
1465 : : !FB &*Shell3at%proj(:,icoeff,ishell))*SC%dij(iatom,katom,gama+1)
1466 : : !FB Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1467 : : !FB& Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1-terme2
1468 : : !FB end do
1469 : : !FB end do
1470 : : !FB end do
1471 : : !FB end do
1472 : : !FB end do
1473 : : end do !iatshell
1474 : : end do !iatom
1475 : : end do !ishell
1476 : : end if !order=4
1477 : :
1478 [ + - ]: 46 : if (order2.or.order3) then
1479 [ - + ]: 46 : ABI_FREE(SS_ref)
1480 : : end if
1481 [ + + ]: 46 : if (order3.or.order4) then
1482 [ + + ]: 584 : do isym=1,Sym%nsym
1483 [ + + ]: 3992 : do itrans=1,6
1484 [ - + ]: 3976 : ABI_FREE(Const3%Sprod(isym,itrans)%SSS)
1485 : : end do
1486 : : end do
1487 [ + - + + : 3424 : ABI_FREE(Const3%Sprod)
- + - + -
+ - + ]
1488 : : end if
1489 : : if (order4) then
1490 : : do isym=1,Sym%nsym
1491 : : do itrans=1,24
1492 : : ABI_FREE(Const4%Sprod(isym,itrans)%SSSS)
1493 : : end do
1494 : : end do
1495 : : ABI_FREE(Const4%Sprod)
1496 : : end if
1497 : :
1498 : : ! Reduce the number of constraints by selecting the non-zero equations
1499 : 46 : write(ab_out,*) '################## Reduce the number of constraints #########################'
1500 : 46 : iconst_new=0
1501 [ + - ]: 46 : if (order2) then
1502 : : ! 1/ For Rotational invariances (1st order)
1503 : 46 : iconst=0
1504 [ - + + - : 26617 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_1st)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1505 [ + + ]: 184 : do alpha=1,3
1506 [ + + ]: 598 : do beta=1,3
1507 : 414 : iconst=iconst+1
1508 [ + + ]: 26571 : vect(:,iconst)=const_rot1st(alpha,beta,:)
1509 : : end do
1510 : : end do
1511 : 46 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_1st,nconst_loc,vect)
1512 [ - + ]: 46 : if (nconst_loc.ne.0) then
1513 [ # # ]: 0 : do iconst_loc=1,nconst_loc
1514 : 0 : iconst_new=iconst_new+1
1515 [ # # ]: 0 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1516 : : end do
1517 : : end if
1518 [ - + ]: 46 : ABI_FREE(vect)
1519 [ - + ]: 46 : ABI_FREE(const_rot1st)
1520 : 46 : Solver%nconst_1st=nconst_loc
1521 : :
1522 : : ! 2/ For Rotational invariances (2nd order)
1523 : 46 : iconst=0
1524 [ - + + - : 374809 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_rot2nd)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1525 [ + + ]: 156 : do iatom=1,natom_unitcell
1526 [ + + ]: 486 : do alpha=1,3
1527 [ + + ]: 1430 : do beta=1,3
1528 [ + + ]: 4290 : do gama=1,3
1529 : 2970 : iconst=iconst+1
1530 [ + + ]: 375615 : vect(:,iconst)=const_rot2nd(alpha,beta,gama,iatom,:)
1531 : : end do
1532 : : end do
1533 : : end do
1534 : : end do
1535 : 46 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_rot2nd,nconst_loc,vect)
1536 [ + + ]: 46 : if (nconst_loc.ne.0) then
1537 [ + + ]: 121 : do iconst_loc=1,nconst_loc
1538 : 108 : iconst_new=iconst_new+1
1539 [ + + ]: 33751 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1540 : : end do
1541 : : end if
1542 [ - + ]: 46 : ABI_FREE(vect)
1543 [ - + ]: 46 : ABI_FREE(const_rot2nd)
1544 : 46 : Solver%nconst_rot2nd=nconst_loc
1545 : :
1546 : : ! 3/ For symetry of the dynamical matrix
1547 : 46 : iconst=0
1548 [ - + + - : 1283719 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_dynmat)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1549 [ + + ]: 156 : do iatom=1,natom_unitcell
1550 [ + + ]: 678 : do jatom=1,natom_unitcell
1551 [ + + ]: 2198 : do alpha=1,3
1552 [ + + ]: 6786 : do beta=1,3
1553 : 4698 : iconst=iconst+1
1554 [ + + ]: 1285101 : vect(:,iconst)=const_dynmat(alpha,beta,iatom,jatom,:)
1555 : : end do
1556 : : end do
1557 : : end do
1558 : : end do
1559 : 46 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_dynmat,nconst_loc,vect)
1560 [ + + ]: 46 : if (nconst_loc.ne.0) then
1561 [ + + ]: 71 : do iconst_loc=1,nconst_loc
1562 : 67 : iconst_new=iconst_new+1
1563 [ + + ]: 27027 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1564 : : end do
1565 : : end if
1566 [ - + ]: 46 : ABI_FREE(vect)
1567 [ - + ]: 46 : ABI_FREE(const_dynmat)
1568 : 46 : Solver%nconst_dynmat=nconst_loc
1569 : :
1570 : : ! 4/ For Huang invariances
1571 : 46 : iconst=0
1572 [ - + + - : 238081 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_huang)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1573 [ + + ]: 184 : do alpha=1,3
1574 [ + + ]: 598 : do beta=1,3
1575 [ + + ]: 1794 : do gama=1,3
1576 [ + + ]: 5382 : do lambda=1,3
1577 : 3726 : iconst=iconst+1
1578 [ + + ]: 239139 : vect(:,iconst)=const_huang(alpha,beta,gama,lambda,:)
1579 : : end do
1580 : : end do
1581 : : end do
1582 : : end do
1583 : 46 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_huang,nconst_loc,vect)
1584 [ + + ]: 46 : if (nconst_loc.ne.0) then
1585 [ + + ]: 76 : do iconst_loc=1,nconst_loc
1586 : 56 : iconst_new=iconst_new+1
1587 [ + + ]: 6927 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1588 : : end do
1589 : : end if
1590 [ - + ]: 46 : ABI_FREE(vect)
1591 [ - + ]: 46 : ABI_FREE(const_huang)
1592 : 46 : Solver%nconst_huang=nconst_loc
1593 : 230 : Solver%nconst_2nd=Solver%nconst_rot2nd+Solver%nconst_dynmat+Solver%nconst_huang
1594 : : end if
1595 : :
1596 [ + + ]: 46 : if (order3) then
1597 : : ! 1/ For acoustic sum rules (3rd order)
1598 : 16 : iconst=0
1599 [ - + + - : 13856896 : ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_asr3rd)) ; vect (:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1600 [ + + ]: 48 : do iatom=1,natom_unitcell
1601 [ + + ]: 5728 : do jatom=1,natom
1602 [ + + ]: 22752 : do alpha=1,3
1603 [ + + ]: 73840 : do beta=1,3
1604 [ + + ]: 221520 : do gama=1,3
1605 : 153360 : iconst=iconst+1
1606 [ + + ]: 13907952 : do ii=1,ntotcoeff
1607 : 13856832 : vect(ii,iconst)=Const3%AsrRot3(iatom,jatom,ii)%ABG(alpha,beta,gama)
1608 : : end do
1609 : : end do
1610 : : end do
1611 : : end do
1612 : : end do
1613 : : end do
1614 : 16 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_asr3rd,nconst_loc,vect)
1615 [ + + ]: 16 : if (nconst_loc.ne.0) then
1616 [ + + ]: 91 : do iconst_loc=1,nconst_loc
1617 : 78 : iconst_new=iconst_new+1
1618 [ + + ]: 9731 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1619 : : end do
1620 : : end if
1621 [ - + ]: 16 : ABI_FREE(vect)
1622 : 16 : Solver%nconst_asr3rd=nconst_loc
1623 : :
1624 : : ! 2/ For Rotational invariances (3rd order)
1625 : 16 : iconst=0
1626 [ - + + - : 41570560 : ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_rot3rd)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1627 [ + + ]: 48 : do iatom=1,natom_unitcell
1628 [ + + ]: 5728 : do jatom=1,natom
1629 [ + + ]: 22752 : do alpha=1,3
1630 [ + + ]: 73840 : do beta=1,3
1631 [ + + ]: 221520 : do gama=1,3
1632 [ + + ]: 664560 : do lambda=1,3
1633 : 460080 : iconst=iconst+1
1634 [ + + ]: 41723856 : do ii=1,ntotcoeff
1635 : 41570496 : vect(ii,iconst)=Const3%AsrRot3(iatom,jatom,ii)%ABGD(alpha,beta,gama,lambda)
1636 : : end do
1637 : : end do
1638 : : end do
1639 : : end do
1640 : : end do
1641 : : end do
1642 : : end do
1643 : 16 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_rot3rd,nconst_loc,vect)
1644 [ + + ]: 16 : if (nconst_loc.ne.0) then
1645 [ + + ]: 129 : do iconst_loc=1,nconst_loc
1646 : 116 : iconst_new=iconst_new+1
1647 [ + + ]: 15805 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1648 : : end do
1649 : : end if
1650 [ - + ]: 16 : ABI_FREE(vect)
1651 [ + + ]: 48 : do iatom=1,natom_unitcell
1652 [ + + ]: 5728 : do jatom=1,natom
1653 [ + + ]: 513248 : do icoeff=1,ntotcoeff
1654 [ - + ]: 507536 : ABI_FREE(Const3%AsrRot3(iatom,jatom,icoeff)%ABG)
1655 [ - + ]: 513216 : ABI_FREE(Const3%AsrRot3(iatom,jatom,icoeff)%ABGD)
1656 : : end do
1657 : : end do
1658 : : end do
1659 [ + - + + : 507552 : ABI_FREE(Const3%AsrRot3)
- + - + -
+ - + ]
1660 : 16 : Solver%nconst_rot3rd=nconst_loc
1661 : 48 : Solver%nconst_3rd=Solver%nconst_asr3rd+Solver%nconst_rot3rd
1662 : : end if
1663 : :
1664 : : if (order4) then
1665 : : ! 1/ For acoustic sum rules (4th order)
1666 : : iconst=0
1667 : : ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_asr4th)) ; vect (:,:)=zero
1668 : : do iatom=1,natom_unitcell
1669 : : do jatom=1,natom
1670 : : do katom=1,natom
1671 : : do alpha=1,3
1672 : : do beta=1,3
1673 : : do gama=1,3
1674 : : do delta=1,3
1675 : : iconst=iconst+1
1676 : : do ii=1,ntotcoeff
1677 : : vect(ii,iconst)=Const4%AsrRot4(iatom,jatom,katom,ii)%ABGD(alpha,beta,gama,delta)
1678 : : end do
1679 : : end do
1680 : : end do
1681 : : end do
1682 : : end do
1683 : : end do
1684 : : end do
1685 : : end do
1686 : : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_asr4th,nconst_loc,vect)
1687 : : if (nconst_loc.ne.0) then
1688 : : do iconst_loc=1,nconst_loc
1689 : : iconst_new=iconst_new+1
1690 : : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1691 : : end do
1692 : : end if
1693 : : ABI_FREE(vect)
1694 : : Solver%nconst_asr4th=nconst_loc
1695 : :
1696 : : !FB! 2/ For Rotational invariances (3rd order)
1697 : : !FB iconst=0
1698 : : !FB ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_rot3rd)) ; vect(:,:)=zero
1699 : : !FB do iatom=1,natom_unitcell
1700 : : !FB do jatom=1,natom
1701 : : !FB do alpha=1,3
1702 : : !FB do beta=1,3
1703 : : !FB do gama=1,3
1704 : : !FB do lambda=1,3
1705 : : !FB iconst=iconst+1
1706 : : !FB vect(:,iconst)=Const4%AsrRot3(iatom,jatom,:)%ABGD(alpha,beta,gama,lambda)
1707 : : !FB end do
1708 : : !FB end do
1709 : : !FB end do
1710 : : !FB end do
1711 : : !FB end do
1712 : : !FB end do
1713 : : !FB call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_rot3rd,nconst_loc,vect)
1714 : : !FB if (nconst_loc.ne.0) then
1715 : : !FB do iconst_loc=1,nconst_loc
1716 : : !FB iconst_new=iconst_new+1
1717 : : !FB Solver%const(iconst_new,:)=vect(:,iconst_loc)
1718 : : !FB end do
1719 : : !FB end if
1720 : : !FB ABI_FREE(vect)
1721 : : !FB Solver%nconst_rot3rd=nconst_loc
1722 : : do iatom=1,natom_unitcell
1723 : : do jatom=1,natom
1724 : : do katom=1,natom
1725 : : do icoeff=1,ntotcoeff
1726 : : ABI_FREE(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGD)
1727 : : !FB ABI_FREE(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGDE)
1728 : : end do
1729 : : end do
1730 : : end do
1731 : : end do
1732 : : ABI_FREE(Const4%AsrRot4)
1733 : : Solver%nconst_rot4th=0
1734 : : Solver%nconst_4th=Solver%nconst_asr4th+Solver%nconst_rot4th
1735 : : end if
1736 : :
1737 : : ! Finalize the orthonormalization
1738 : : !FB nconst=iconst_new
1739 : : !FB ABI_MALLOC(vect,(ntotcoeff,nconst)) ; vect(:,:)=zero
1740 : : !FB do iconst=1,nconst
1741 : : !FB vect(:,iconst)=Solver%const(iconst,:)
1742 : : !FB end do
1743 : : !FB ABI_FREE(Solver%const)
1744 : : !FB call tdep_calc_orthonorm(ntotcoeff,nconst,nconst_loc,vect)
1745 : : !FB Solver%ntotconst=nconst_loc
1746 : : !FB if (nconst_loc.ne.0) then
1747 : : !FB ABI_MALLOC(Solver%const ,(Solver%ntotconst,ntotcoeff)); Solver%const (:,:)=0.d0
1748 : : !FB do iconst_loc=1,nconst_loc
1749 : : !FB Solver%const(iconst_loc,:)=vect(:,iconst_loc)
1750 : : !FB end do
1751 : : !FB end if
1752 : : !FB ABI_FREE(vect)
1753 : 46 : Solver%ntotconst=iconst_new
1754 [ + - ]: 46 : if (MPIdata%iam_master) then
1755 : 46 : open(unit=16,file=trim(Invar%output_prefix)//'_constraints.dat')
1756 : 46 : write(16,*) ' ======== Constraints at the 1st order (Rotational Invariances) ========'
1757 : 46 : write(16,*) ' Number of constraints =',Solver%nconst_1st
1758 : 46 : write(16,*) ' ======== Constraints at the 2nd order (Rotational Invariances) ========'
1759 : 46 : write(16,*) ' Number of constraints =',Solver%nconst_rot2nd
1760 : 46 : write(16,*) ' ======== Constraints at the 2nd order (Dynamical Matrix) =============='
1761 : 46 : write(16,*) ' Number of constraints =',Solver%nconst_dynmat
1762 : 46 : write(16,*) ' ======== Constraints at the 2nd order (Huang) ========================='
1763 : 46 : write(16,*) ' Number of constraints =',Solver%nconst_huang
1764 [ + + ]: 46 : if (Solver%order.ge.3) then
1765 : 16 : write(16,*) ' ======== Constraints at the 3rd order (Acoustic sum rules) ============'
1766 : 16 : write(16,*) ' Number of constraints =',Solver%nconst_asr3rd
1767 : 16 : write(16,*) ' ======== Constraints at the 3rd order (Rotational Invariances) ========'
1768 : 16 : write(16,*) ' Number of constraints =',Solver%nconst_rot3rd
1769 : : end if
1770 [ + + ]: 46 : if (Solver%order.ge.4) then
1771 : 9 : write(16,*) ' ======== Constraints at the 4th order (Acoustic sum rules) ============'
1772 : 9 : write(16,*) ' Number of constraints =',Solver%nconst_asr4th
1773 : 9 : write(16,*) ' ======== Constraints at the 4th order (Rotational Invariances) ========'
1774 : 9 : write(16,*) ' Number of constraints =',Solver%nconst_rot4th
1775 : : end if
1776 : 46 : write(16,*) ' ======================================================================='
1777 : 46 : write(16,*) ' Total number of constraints =',Solver%ntotconst
1778 : 46 : close(16)
1779 : : end if
1780 : :
1781 [ - + - + : 46 : end subroutine tdep_calc_constraints
- + - + -
+ - + - -
- - - - -
- - + - -
- - - - -
- ]
1782 : : !!***
1783 : :
1784 : : !====================================================================================================
1785 : :
1786 : : !!****f* m_tdep_solver/tdep_calc_phi1
1787 : : !! NAME
1788 : : !! tdep_calc_phi1
1789 : : !!
1790 : : !! FUNCTION
1791 : : !! Build the first-order IFC vector from the solved coefficients (Solver%theta).
1792 : : !! The IFC of the reference atom of each shell is obtained from the
1793 : : !! independent coefficients, transformed to every atom of the shell using
1794 : : !! the symmetry operations, and the acoustic sum rule is imposed.
1795 : : !!
1796 : : !! INPUTS
1797 : : !! Solver = TDEP Solver object, containing the solved coefficients.
1798 : : !! Shell1at = Shell object for first-order interactions.
1799 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1800 : : !!
1801 : : !! OUTPUT
1802 : : !! Phi1(3*natom) = First-order IFC.
1803 : : !!
1804 : : !! SOURCE
1805 : :
1806 : 45 : subroutine tdep_calc_phi1(Solver,Shell1at,Sym,Phi1)
1807 : :
1808 : : class(tdep_Solver_type),intent(in) :: Solver
1809 : : type(Shell_type),intent(in) :: Shell1at
1810 : : type(Symmetries_type),intent(in) :: Sym
1811 : : double precision,intent(out) :: Phi1(3*Solver%natom)
1812 : :
1813 : : integer :: ishell,isym,iatom,ncoeff,ncoeff_prev
1814 : : integer :: nshell,ii,iatshell,iat_mod
1815 : 45 : double precision,allocatable :: Phi1_coeff(:)
1816 : 45 : double precision,allocatable :: Phi1_3(:),Phi1_ref(:,:)
1817 : :
1818 : 45 : nshell=Shell1at%nshell
1819 [ + - - - : 135 : ABI_MALLOC(Phi1_coeff,(Solver%ncoeff1st))
+ - - + -
+ - + ]
1820 [ - + ]: 45 : Phi1_coeff(:) = Solver%theta(1:Solver%ncoeff1st)
1821 [ - + + - : 391 : ABI_MALLOC(Phi1_ref,(3,nshell)); Phi1_ref(:,:)=zero
+ - + - -
+ - + + +
+ + ]
1822 [ + - + + ]: 225 : ABI_MALLOC(Phi1_3,(3)) ; Phi1_3(:)=0.d0
1823 [ + + ]: 109 : do ishell=1,nshell
1824 : : ! Build the 3x3 IFC per shell
1825 : 64 : ncoeff =Shell1at%ncoeff(ishell)
1826 : 64 : ncoeff_prev=Shell1at%ncoeff_prev(ishell)
1827 [ + + ]: 256 : do ii=1,3
1828 [ - + ]: 256 : Phi1_ref(ii,ishell)=sum(Shell1at%proj(ii,1:ncoeff,ishell)*Phi1_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
1829 : : end do
1830 : : ! Build the vector-IFC of an atom in this shell
1831 [ - + ]: 64 : if (Shell1at%neighbours(1,ishell)%n_interactions.eq.0) cycle
1832 [ - - ]: 45 : do iatshell=1,Shell1at%neighbours(1,ishell)%n_interactions
1833 : 0 : iatom=Shell1at%neighbours(1,ishell)%atomj_in_shell(iatshell)
1834 : 0 : isym =Shell1at%neighbours(1,ishell)%sym_in_shell(iatshell)
1835 [ # # ]: 0 : do ii=1,3
1836 [ # # ]: 0 : Phi1_3(ii)=sum(Sym%S_ref(ii,:,isym,1)*Phi1_ref(:,ishell))
1837 : : end do
1838 [ - - ]: 64 : Phi1((iatom-1)*3+1:(iatom-1)*3+3)=Phi1_3(:)
1839 : : end do !iatshell
1840 : : end do !ishell
1841 : : ! Acoustic sum rule
1842 [ + + ]: 180 : do ii=1,3
1843 [ + + ]: 20478 : do iatom=1,Solver%natom
1844 : 20298 : iat_mod=mod(iatom+Solver%natom_unitcell-1,Solver%natom_unitcell)+1
1845 [ + + ]: 20298 : if (iat_mod==1) cycle
1846 : 20433 : Phi1((iatom-iat_mod+1)*3+ii)=Phi1((iatom-iat_mod+1)*3+ii)-Phi1((iatom-1)*3+ii)
1847 : : end do
1848 : : end do
1849 [ - + ]: 45 : ABI_FREE(Phi1_coeff)
1850 [ - + ]: 45 : ABI_FREE(Phi1_3)
1851 [ - + ]: 45 : ABI_FREE(Phi1_ref)
1852 : :
1853 : : ! Remove the rounding errors before writing (for non regression testing purposes)
1854 [ + + ]: 20343 : do ii=1,3*Solver%natom
1855 [ + - ]: 20343 : if (abs(Phi1(ii)).lt.tol8) Phi1(ii)=zero
1856 : : end do
1857 : :
1858 : 45 : end subroutine tdep_calc_phi1
1859 : : !!***
1860 : :
1861 : : !====================================================================================================
1862 : :
1863 : : !!****f* m_tdep_solver/tdep_calc_phi2
1864 : : !! NAME
1865 : : !! tdep_calc_phi2
1866 : : !!
1867 : : !! FUNCTION
1868 : : !! Build the second-order IFC matrix from the solved coefficients (Solver%theta).
1869 : : !! The 3x3 IFC of the reference pair of each shell is obtained from the
1870 : : !! independent coefficients, transformed to every pair of the shell using
1871 : : !! the symmetry operations (see tdep_build_phi2_33), the matrix is
1872 : : !! symmetrized, and the acoustic sum rule is imposed.
1873 : : !!
1874 : : !! INPUTS
1875 : : !! Solver = TDEP Solver object, containing the solved coefficients.
1876 : : !! Shell2at = Shell object for second-order interactions.
1877 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
1878 : : !!
1879 : : !! OUTPUT
1880 : : !! Phi2(3*natom,3*natom) = Second-order IFC, in the supercell.
1881 : : !!
1882 : : !! SOURCE
1883 : :
1884 : 45 : subroutine tdep_calc_phi2(Solver,Shell2at,Sym,Phi2)
1885 : :
1886 : : class(tdep_Solver_type),intent(in) :: Solver
1887 : : type(Shell_type),intent(in) :: Shell2at
1888 : : type(Symmetries_type),intent(in) :: Sym
1889 : : double precision,intent(out) :: Phi2(3*Solver%natom,3*Solver%natom)
1890 : :
1891 : : integer :: ishell,isym,eatom,fatom,ncoeff,ncoeff_prev
1892 : : integer :: nshell,ii,jj,kk,kappa,iatshell,itrans
1893 : 45 : double precision,allocatable :: Phi2_33(:,:),Phi2_ref(:,:,:)
1894 : 45 : double precision,allocatable :: Phi2_coeff(:)
1895 : :
1896 : 45 : nshell=Shell2at%nshell
1897 [ - + + - : 135 : ABI_MALLOC(Phi2_coeff,(Solver%ncoeff2nd))
+ - + - -
+ - + ]
1898 [ + + ]: 2059 : Phi2_coeff(:) = Solver%theta(Solver%ncoeff1st+1:Solver%ncoeff1st+Solver%ncoeff2nd)
1899 [ - + + - : 6648 : ABI_MALLOC(Phi2_ref,(3,3,nshell)); Phi2_ref(:,:,:)=zero
+ - + - -
+ - + + +
+ + + + ]
1900 [ + - + + : 630 : ABI_MALLOC(Phi2_33,(3,3)) ; Phi2_33(:,:)=0.d0
+ + ]
1901 [ + + ]: 546 : do ishell=1,nshell
1902 : : ! Build the 3x3 IFC per shell
1903 : 501 : ncoeff =Shell2at%ncoeff(ishell)
1904 : 501 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)
1905 : 501 : kappa=0
1906 [ + + ]: 2004 : do ii=1,3
1907 [ + + ]: 6513 : do jj=1,3
1908 : 4509 : kappa=kappa+1
1909 [ + + ]: 24138 : Phi2_ref(ii,jj,ishell)=sum(Shell2at%proj(kappa,1:ncoeff,ishell)*Phi2_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
1910 : : end do
1911 : : end do
1912 [ + + ]: 72724 : do eatom=1,Solver%natom
1913 : : ! Build the 3x3 IFC of an atom in this shell
1914 [ + + ]: 72178 : if (Shell2at%neighbours(eatom,ishell)%n_interactions.eq.0) cycle
1915 [ + + ]: 507937 : do iatshell=1,Shell2at%neighbours(eatom,ishell)%n_interactions
1916 : 446778 : fatom=Shell2at%neighbours(eatom,ishell)%atomj_in_shell(iatshell)
1917 : 446778 : isym =Shell2at%neighbours(eatom,ishell)%sym_in_shell(iatshell)
1918 : 446778 : itrans=Shell2at%neighbours(eatom,ishell)%transpose_in_shell(iatshell)
1919 [ + + ]: 446778 : if (fatom.lt.eatom) cycle
1920 : 226772 : call tdep_build_phi2_33(isym,Phi2_ref(:,:,ishell),Phi2_33,Sym,itrans)
1921 : : ! Symetrization of the Phi2 matrix
1922 [ + + + + ]: 2948036 : Phi2((eatom-1)*3+1:(eatom-1)*3+3,3*(fatom-1)+1:3*(fatom-1)+3)=Phi2_33(:,:)
1923 [ + + ]: 979266 : do ii=1,3
1924 [ + + ]: 3168042 : do jj=1,3
1925 : 2721264 : Phi2((fatom -1)*3+ii,3*(eatom-1)+jj)=Phi2_33(jj,ii)
1926 : : end do
1927 : : end do
1928 : : end do !iatshell
1929 : : end do !eatom
1930 : : end do !ishell
1931 : : ! Acoustic sum rule
1932 [ + + ]: 6811 : do eatom=1,Solver%natom
1933 [ + + ]: 27109 : do jj=1,3
1934 [ + + ]: 87958 : do kk=1,3
1935 [ + + ]: 11684172 : do fatom=1,Solver%natom
1936 [ + + ]: 11602980 : if (fatom==eatom) cycle
1937 : : Phi2((eatom-1)*3+jj,(eatom-1)*3+kk)=Phi2((eatom-1)*3+jj,3*(eatom-1)+kk)&
1938 : 11663874 : & -Phi2((eatom-1)*3+jj,3*(fatom-1)+kk)
1939 : : enddo
1940 : : enddo
1941 : : enddo
1942 : : enddo
1943 [ - + ]: 45 : ABI_FREE(Phi2_coeff)
1944 [ - + ]: 45 : ABI_FREE(Phi2_33)
1945 [ - + ]: 45 : ABI_FREE(Phi2_ref)
1946 : :
1947 : : ! Remove the rounding errors before writing (for non regression testing purposes)
1948 [ + + ]: 20343 : do ii=1,3*Solver%natom
1949 [ + + ]: 11623323 : do jj=1,3*Solver%natom
1950 [ + + ]: 11623278 : if (abs(Phi2(ii,jj)).lt.tol8) Phi2(ii,jj)=zero
1951 : : end do
1952 : : end do
1953 : :
1954 : 45 : end subroutine tdep_calc_phi2
1955 : : !!***
1956 : :
1957 : : !====================================================================================================
1958 : :
1959 : : !!****f* m_tdep_solver/tdep_calc_phi3ref
1960 : : !! NAME
1961 : : !! tdep_calc_phi3ref
1962 : : !!
1963 : : !! FUNCTION
1964 : : !! Build the third-order IFC reference tensor (per shell) from the solved
1965 : : !! coefficients (Solver%theta).
1966 : : !!
1967 : : !! INPUTS
1968 : : !! Solver = TDEP Solver object, containing the solved coefficients.
1969 : : !! Shell3at = Shell object for third-order interactions.
1970 : : !!
1971 : : !! OUTPUT
1972 : : !!
1973 : : !! SIDE EFFECTS
1974 : : !! Phi3_ref(3,3,3,nshell) = Third-order IFC of the reference triplet of each shell.
1975 : : !!
1976 : : !! SOURCE
1977 : :
1978 : 16 : subroutine tdep_calc_phi3ref(Solver,Shell3at,Phi3_ref)
1979 : :
1980 : : class(tdep_Solver_type),intent(in) :: Solver
1981 : : type(Shell_type),intent(in) :: Shell3at
1982 : : double precision, intent(inout) :: Phi3_ref(3,3,3,Shell3at%nshell)
1983 : :
1984 : : integer :: ishell,ncoeff,ncoeff_prev
1985 : : integer :: ii,jj,kk,kappa
1986 : 16 : double precision, allocatable :: Phi3_coeff(:)
1987 : :
1988 [ + + + - : 273 : ABI_CALLOC(Phi3_coeff,(Solver%ncoeff3rd))
+ - + + -
+ - + +
+ ]
1989 [ + + ]: 241 : Phi3_coeff(:) = Solver%theta(Solver%ncoeff1st+Solver%ncoeff2nd+1:Solver%ncoeff1st+Solver%ncoeff2nd+Solver%ncoeff3rd)
1990 : :
1991 [ + + ]: 73 : do ishell=1,Shell3at%nshell
1992 : : ! Build the 3x3x3 IFC per shell
1993 : 57 : ncoeff =Shell3at%ncoeff(ishell)
1994 : 57 : ncoeff_prev=Shell3at%ncoeff_prev(ishell)
1995 : 57 : kappa=0
1996 [ + + ]: 228 : do ii=1,3
1997 [ + + ]: 741 : do jj=1,3
1998 [ + + ]: 2223 : do kk=1,3
1999 : 1539 : kappa=kappa+1
2000 [ + + ]: 8127 : Phi3_ref(ii,jj,kk,ishell)=sum(Shell3at%proj(kappa,1:ncoeff,ishell)*Phi3_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
2001 : : end do
2002 : : end do
2003 : : end do
2004 : : ! Remove the rounding errors before writing (for non regression testing purposes)
2005 [ + + ]: 244 : do ii=1,3
2006 [ + + ]: 741 : do jj=1,3
2007 [ + + ]: 2223 : do kk=1,3
2008 [ + + ]: 2052 : if (abs(Phi3_ref(ii,jj,kk,ishell)).lt.tol8) Phi3_ref(ii,jj,kk,ishell)=zero
2009 : : end do
2010 : : end do
2011 : : end do
2012 : : end do
2013 [ - + ]: 16 : ABI_FREE(Phi3_coeff)
2014 : :
2015 : 16 : end subroutine tdep_calc_phi3ref
2016 : : !!***
2017 : :
2018 : : !====================================================================================================
2019 : :
2020 : : !!****f* m_tdep_solver/tdep_calc_phi4ref
2021 : : !! NAME
2022 : : !! tdep_calc_phi4ref
2023 : : !!
2024 : : !! FUNCTION
2025 : : !! Build the fourth-order IFC reference tensor (per shell) from the solved
2026 : : !! coefficients (Solver%theta).
2027 : : !!
2028 : : !! INPUTS
2029 : : !! Solver = TDEP Solver object, containing the solved coefficients.
2030 : : !! Shell4at = Shell object for fourth-order interactions.
2031 : : !!
2032 : : !! OUTPUT
2033 : : !!
2034 : : !! SIDE EFFECTS
2035 : : !! Phi4_ref(3,3,3,3,nshell) = Fourth-order IFC of the reference quadruplet of each shell.
2036 : : !!
2037 : : !! SOURCE
2038 : :
2039 : 9 : subroutine tdep_calc_phi4ref(Solver,Shell4at,Phi4_ref)
2040 : :
2041 : : class(tdep_Solver_type),intent(in) :: Solver
2042 : : type(Shell_type),intent(in) :: Shell4at
2043 : : double precision, intent(inout) :: Phi4_ref(3,3,3,3,Shell4at%nshell)
2044 : :
2045 : : integer :: ishell,ncoeff,ncoeff_prev
2046 : : integer :: ii,jj,kk,ll,kappa
2047 : 9 : double precision, allocatable :: Phi4_coeff(:)
2048 : :
2049 [ - + + - : 583 : ABI_CALLOC(Phi4_coeff, (Solver%ncoeff4th))
+ - + - -
+ - + +
+ ]
2050 [ + + ]: 565 : Phi4_coeff(:) = Solver%theta(Solver%ncoeff1st+Solver%ncoeff2nd+Solver%ncoeff3rd+1:Solver%ntotcoeff)
2051 : :
2052 [ + + ]: 58 : do ishell=1,Shell4at%nshell
2053 : : ! Build the 3x3x3x3 IFC per shell
2054 : 49 : ncoeff =Shell4at%ncoeff(ishell)
2055 : 49 : ncoeff_prev=Shell4at%ncoeff_prev(ishell)
2056 : 49 : kappa=0
2057 [ + + ]: 196 : do ii=1,3
2058 [ + + ]: 637 : do jj=1,3
2059 [ + + ]: 1911 : do kk=1,3
2060 [ + + ]: 5733 : do ll=1,3
2061 : 3969 : kappa=kappa+1
2062 [ + + ]: 50328 : Phi4_ref(ii,jj,kk,ll,ishell)=sum(Shell4at%proj(kappa,1:ncoeff,ishell)*Phi4_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
2063 : : end do
2064 : : end do
2065 : : end do
2066 : : end do
2067 : : ! Remove the rounding errors before writing (for non regression testing purposes)
2068 [ + + ]: 205 : do ii=1,3
2069 [ + + ]: 637 : do jj=1,3
2070 [ + + ]: 1911 : do kk=1,3
2071 [ + + ]: 5733 : do ll=1,3
2072 [ + + ]: 5292 : if (abs(Phi4_ref(ii,jj,kk,ll,ishell)).lt.tol8) Phi4_ref(ii,jj,kk,ll,ishell)=zero
2073 : : end do
2074 : : end do
2075 : : end do
2076 : : end do
2077 : : end do
2078 : :
2079 [ - + ]: 9 : ABI_FREE(Phi4_coeff)
2080 : :
2081 : 9 : end subroutine tdep_calc_phi4ref
2082 : : !!***
2083 : :
2084 : : !====================================================================================================
2085 : :
2086 [ # # # # : 0 : end module m_tdep_solver
# # # # #
# # # # #
# # # # #
# ]
2087 : : !!***
|
