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_constraints, only : Constraints_type, tdep_calc_orthonorm
38 : :
39 : : implicit none
40 : :
41 : : type tdep_Solver_type
42 : :
43 : : integer :: order
44 : : ! Maximum IFC order (2, 3, or 4)
45 : :
46 : : integer :: natom
47 : : ! Number of atoms in the supercell
48 : :
49 : : integer :: natom_unitcell
50 : : ! Number of atoms in the unitcell
51 : :
52 : : integer :: my_nstep
53 : : ! Number of MD steps
54 : :
55 : : integer :: ntotcoeff
56 : : ! Total number of coefficients
57 : :
58 : : integer :: ntotconst
59 : : ! Total number of constraints
60 : :
61 : : integer :: ncoeff1st
62 : : ! Number of 1st order coefficients
63 : :
64 : : integer :: ncoeff2nd
65 : : ! Number of 2nd order coefficients
66 : :
67 : : integer :: ncoeff3rd
68 : : ! Number of 3rd order coefficients
69 : :
70 : : integer :: ncoeff4th
71 : : ! Number of 4th order coefficients
72 : :
73 : : integer :: nconst_1st
74 : : ! Number of 1st order constraints
75 : :
76 : : integer :: nconst_2nd
77 : : ! Number of 2nd order constraints
78 : :
79 : : integer :: nconst_3rd
80 : : ! Number of 3rd order constraints
81 : :
82 : : integer :: nconst_4th
83 : : ! Number of 4th order constraints
84 : :
85 : : integer :: nconst_rot2nd
86 : : integer :: nconst_huang
87 : : integer :: nconst_dynmat
88 : : integer :: nconst_rot3rd
89 : : integer :: nconst_asr3rd
90 : : integer :: nconst_rot4th
91 : : integer :: nconst_asr4th
92 : :
93 : : double precision, allocatable :: fcoeff(:,:)
94 : : ! fcoeff(3*natom*my_nstep,ntotcoeff)
95 : : ! All the cartesian displacement matrices
96 : : ! at every order (u, u*u, u*u*u, u*u*u*u).
97 : :
98 : : double precision, allocatable :: const(:,:)
99 : : ! const(ntotconst,ntotcoeff)
100 : : ! The constraint matrices.
101 : :
102 : : double precision, allocatable :: theta(:)
103 : : ! theta(ntotcoeff)
104 : : ! The IFC coefficients at all orders, as a flat array.
105 : :
106 : : double precision, allocatable :: Forces(:)
107 : : ! Forces(3*natom*my_nstep)
108 : : ! The cartesian forces for all configurations, as a flat array.
109 : : ! These are weighted by the number of configurations.
110 : : ! In case of dipole-dipole interaction, the long-range part of the
111 : : ! forces should be removed.
112 : :
113 : : end type tdep_Solver_type
114 : :
115 : : public :: tdep_solver_init
116 : : public :: tdep_solver_free
117 : : public :: tdep_solver_set_residual_forces
118 : : public :: tdep_calc_phi1fcoeff
119 : : public :: tdep_calc_phi2fcoeff
120 : : public :: tdep_calc_phi3fcoeff
121 : : public :: tdep_calc_phi4fcoeff
122 : : public :: tdep_calc_MoorePenrose
123 : : public :: tdep_calc_constraints
124 : :
125 : : contains
126 : :
127 : : !=====================================================================================================
128 : :
129 : 44 : subroutine tdep_solver_init(Solver, Invar, Shell1at, Shell2at, Shell3at, Shell4at)
130 : :
131 : : type(tdep_Solver_type), intent(inout) :: Solver
132 : : type(atdep_dataset_type),intent(in) :: Invar
133 : : type(Shell_type),intent(in) :: Shell1at, Shell2at, Shell3at, Shell4at
134 : :
135 : 44 : Solver%order = Invar%order
136 : 44 : Solver%natom = Invar%natom
137 : 44 : Solver%natom_unitcell = Invar%natom_unitcell
138 : 44 : Solver%my_nstep = Invar%my_nstep
139 : :
140 : : !Rotational invariances (1st order)
141 : : ! constraints = 3
142 : 44 : Solver%nconst_1st = 3**2
143 : :
144 : : !Rotational invariances (2nd order) + Symetry of the Dynamical Matrix + Huang invariances
145 : : ! constraints = natom*3**2 + (3*natom_unitcell)**2 + 3**4
146 : 44 : Solver%nconst_rot2nd = 3**3*Solver%natom_unitcell
147 : 44 : Solver%nconst_dynmat = (3*Solver%natom_unitcell)**2
148 : 44 : Solver%nconst_huang = 3**4
149 : 44 : Solver%nconst_2nd = Solver%nconst_rot2nd + Solver%nconst_dynmat + Solver%nconst_huang
150 : :
151 : : !Rotational invariances (3rd order) + acoustic sum rules (3rd order)
152 : 44 : Solver%nconst_3rd=0
153 [ + + ]: 44 : if (Solver%order.ge.3) then
154 : : ! constraints = natom_unitcell*natom*3**4 + natom_unitcell*natom*3**3
155 : 14 : Solver%nconst_rot3rd = 3**4 * Solver%natom_unitcell * Solver%natom
156 : 14 : Solver%nconst_asr3rd = 3**3 * Solver%natom_unitcell * Solver%natom
157 : 14 : Solver%nconst_3rd = Solver%nconst_rot3rd + Solver%nconst_asr3rd
158 : : end if
159 : :
160 : : !Rotational invariances (4th order) + acoustic sum rules (4th order)
161 : 44 : Solver%nconst_4th=0
162 [ + + ]: 44 : if (Solver%order.ge.4) then
163 : : ! constraints = natom_unitcell*natom**2*3**5 + natom_unitcell*natom**2*3**4
164 : : !FB Solver%nconst_rot4th = 3**5*natom_unitcell*natom**2
165 : : !FB4TH Solver%nconst_asr4th = 3**4*natom_unitcell*natom**2
166 : 8 : Solver%nconst_rot4th = 0
167 : 8 : Solver%nconst_asr4th = 0
168 : : Solver%nconst_4th = Solver%nconst_rot4th + Solver%nconst_asr4th
169 : : end if
170 : :
171 : 44 : Solver%ncoeff1st = Shell1at%ntotcoeff
172 : 44 : Solver%ncoeff2nd = Shell2at%ntotcoeff
173 : 44 : Solver%ncoeff3rd = 0
174 : 44 : Solver%ncoeff4th = 0
175 [ + + ]: 44 : if (Solver%order.ge.3) Solver%ncoeff3rd = Shell3at%ntotcoeff
176 [ + + ]: 44 : if (Solver%order.ge.4) Solver%ncoeff4th = Shell4at%ntotcoeff
177 : :
178 : 44 : Solver%ntotcoeff = Solver%ncoeff1st + Solver%ncoeff2nd + Solver%ncoeff3rd + Solver%ncoeff4th
179 : 44 : Solver%ntotconst = Solver%nconst_1st + Solver%nconst_2nd + Solver%nconst_3rd + Solver%nconst_4th
180 : :
181 [ + - - + : 25945165 : ABI_CALLOC(Solver%fcoeff, (3*Solver%natom*Solver%my_nstep,Solver%ntotcoeff))
+ - - + +
- + - + -
- + - + -
+ + + +
+ ]
182 [ - + + - : 50744071 : ABI_CALLOC(Solver%const, (Solver%ntotconst,Solver%ntotcoeff))
- + + - +
- + - - +
- + - + +
+ + + ]
183 [ - + + - : 417588 : ABI_CALLOC(Solver%Forces, (3*Solver%natom*Solver%my_nstep))
+ - + - -
+ - + - +
+ + ]
184 [ - + + - : 2873 : ABI_CALLOC(Solver%theta, (Solver%ntotcoeff))
+ - + - -
+ - + - +
+ + ]
185 : :
186 : 44 : end subroutine tdep_solver_init
187 : :
188 : : !=====================================================================================================
189 : :
190 : 44 : subroutine tdep_solver_free(Solver)
191 : :
192 : : type(tdep_Solver_type), intent(inout) :: Solver
193 : :
194 [ - + ]: 44 : ABI_FREE(Solver%fcoeff)
195 [ - + ]: 44 : ABI_FREE(Solver%const)
196 [ - + ]: 44 : ABI_FREE(Solver%Forces)
197 [ - + ]: 44 : ABI_FREE(Solver%theta)
198 : :
199 : 44 : end subroutine tdep_solver_free
200 : :
201 : : !=====================================================================================================
202 : :
203 : : !!***
204 : : !!****f* ABINIT/m_tdep_solver/tdep_solver_set_residual_forces
205 : : !! NAME
206 : : !! tdep_solver_set_residual_forces
207 : : !!
208 : : !! FUNCTION
209 : : !! Compute the residual forces into the solver, that is, the differences
210 : : !! between the MD forces and the ones from the TDEP model.
211 : : !! The model forces may be zero, but they are non-zero when we remove
212 : : !! the long-range part of the IFC forces, or when we solve each IFC order
213 : : !! successively, and we remove the forces from the previously computed orders.
214 : : !!
215 : : !! INPUTS
216 : : !! Solver = TDEP Solver object that will compute the coefficients.
217 : : !! MD = TDEP Sampling object containing the input positions and forces.
218 : : !! Model = TDEP Model object containing the IFC and corresponding forces.
219 : : !!
220 : : !! OUTPUT
221 : : !!
222 : : !! SIDE EFFECTS
223 : : !! The following quantities in MD are computed:
224 : : !!
225 : : !! NOTES
226 : : !!
227 : : !! SOURCE
228 : :
229 : 47 : subroutine tdep_solver_set_residual_forces(Solver, MD, Model)
230 : :
231 : : type(tdep_Solver_type), intent(inout) :: Solver
232 : : type(tdep_Sampling_type), intent(in) :: MD
233 : : type(tdep_Model_type), intent(in) :: Model
234 : :
235 : : integer :: ii,jj,istep,iatom
236 : :
237 [ + + ]: 1059 : do istep=1,MD%my_nstep
238 [ + + ]: 151731 : do iatom=1,MD%natom
239 [ + + ]: 603700 : do ii=1,3
240 : 452016 : jj = ii + 3*(iatom-1) + 3*MD%natom*(istep-1)
241 : 602688 : Solver%Forces(jj) = (MD%Forces(jj) - Model%Forces(jj)) * MD%weights(istep)
242 : : end do
243 : : end do
244 : : end do
245 : :
246 : 47 : end subroutine tdep_solver_set_residual_forces
247 : :
248 : : !====================================================================================================
249 : :
250 : 43 : subroutine tdep_calc_phi1fcoeff(Solver,Invar,Shell1at,Sym)
251 : :
252 : : type(tdep_Solver_type), intent(inout) :: Solver
253 : : type(atdep_dataset_type),intent(in) :: Invar
254 : : type(Shell_type),intent(in) :: Shell1at
255 : : type(Symmetries_type),intent(in) :: Sym
256 : :
257 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,iatshell,iat_mod
258 : : integer :: icoeff,isym,mu,iatref
259 : : double precision :: terme
260 : :
261 : 43 : write(Invar%stdout,*) ' '
262 : 43 : write(Invar%stdout,*) '#############################################################################'
263 : 43 : write(Invar%stdout,*) '############## Fill the matrices used in the pseudo-inverse #################'
264 : 43 : write(Invar%stdout,*) '#############################################################################'
265 : :
266 : 43 : write(Invar%stdout,*) ' Compute the coefficients (at the 1st order) used in the Moore-Penrose...'
267 [ + + ]: 105 : do ishell=1,Shell1at%nshell
268 [ - + ]: 62 : if (Shell1at%neighbours(1,ishell)%n_interactions.eq.0) cycle
269 [ - - ]: 43 : do iatshell=1,Shell1at%neighbours(1,ishell)%n_interactions
270 : 0 : iatom=Shell1at%neighbours(1,ishell)%atomj_in_shell(iatshell)
271 : 0 : iat_mod=mod(iatom+Invar%natom_unitcell-1,Invar%natom_unitcell)+1
272 [ # # ]: 0 : if (iat_mod==1) cycle
273 : 0 : iatref=Shell1at%iatref(ishell)
274 : 0 : isym=Shell1at%neighbours(1,ishell)%sym_in_shell(iatshell)
275 : 0 : ncoeff =Shell1at%ncoeff(ishell)
276 : 0 : ncoeff_prev=Shell1at%ncoeff_prev(ishell)
277 [ - - ]: 62 : do mu=1,3
278 [ # # ]: 0 : do icoeff=1,ncoeff
279 [ # # ]: 0 : terme=sum(Sym%S_ref(mu,:,isym,1)*Shell1at%proj(:,icoeff,ishell))
280 [ # # ]: 0 : do istep=1,Invar%my_nstep
281 : : Solver%fcoeff(mu+3*(iatom-1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev)= &
282 : 0 : & Solver%fcoeff(mu+3*(iatom-1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev)+terme
283 : : ! Add all the other contributions, when iat_mod==1 (due to ASR)
284 : : Solver%fcoeff(mu+3*(iatom-iat_mod+1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev)= &
285 : 0 : & Solver%fcoeff(mu+3*(iatom-iat_mod+1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev)-terme
286 : : end do !istep
287 : : end do
288 : : end do
289 : : end do !iatshell
290 : : end do !ishell
291 : 43 : write(Invar%stdout,*) ' ------- achieved'
292 : :
293 : 43 : end subroutine tdep_calc_phi1fcoeff
294 : :
295 : : !====================================================================================================
296 : :
297 : 43 : subroutine tdep_calc_phi2fcoeff(Solver,Invar,Shell2at,Sym,MD)
298 : :
299 : : type(tdep_Solver_type), intent(inout) :: Solver
300 : : type(atdep_dataset_type),intent(in) :: Invar
301 : : type(Shell_type),intent(in) :: Shell2at
302 : : type(Symmetries_type),intent(in) :: Sym
303 : : type(tdep_Sampling_type), intent(in) :: MD
304 : :
305 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,jatom,iatshell
306 : : integer :: icoeff,isym
307 : : integer :: mu,nu,alpha,beta,itrans
308 : : double precision :: terme,temp
309 : : double precision :: udiff(3),SSu(3,9)
310 : 43 : double precision, allocatable :: SS_ref(:,:,:,:,:)
311 : :
312 : : ! For each couple of atoms, transform the Phi2 (3x3) ifc matrix using the symetry operation (S)
313 : : ! Note: iatom=1 is excluded in order to take into account the atomic sum rule (see below)
314 [ - + + - : 358210 : ABI_MALLOC(SS_ref,(3,9,3,Sym%nsym,2)); SS_ref(:,:,:,:,:)=zero
+ - + - +
- - + - +
+ + + + +
+ + + +
+ ]
315 [ + + ]: 1641 : do isym=1,Sym%nsym
316 [ + + ]: 6435 : do mu=1,3
317 [ + + ]: 20774 : do alpha=1,3
318 [ + + ]: 62322 : do nu=1,3
319 [ + + ]: 186966 : do beta=1,3
320 : 129438 : temp=Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu,beta,isym,1)
321 : 129438 : SS_ref(mu,beta+(alpha-1)*3,nu,isym,1)=temp
322 : 172584 : SS_ref(mu,alpha+(beta-1)*3,nu,isym,2)=temp
323 : : end do
324 : : end do
325 : : end do
326 : : end do
327 : : end do
328 : :
329 : 43 : write(Invar%stdout,*) ' Compute the coefficients (at the 2nd order) used in the Moore-Penrose...'
330 [ + + ]: 522 : do ishell=1,Shell2at%nshell
331 [ + + ]: 68344 : do iatom=1,Invar%natom
332 [ + + ]: 67822 : if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
333 [ + + ]: 462411 : do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
334 : 405630 : jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
335 [ + + ]: 405630 : if (iatom==jatom) cycle
336 : 399260 : isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
337 : 399260 : itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
338 : 399260 : ncoeff =Shell2at%ncoeff(ishell)
339 : 399260 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
340 : :
341 [ + + ]: 8424474 : do istep=1,Invar%my_nstep
342 : : ! In order to impose the acoustic sum rule we use (u(j)-u(i))==u_j^\nu
343 : 7957392 : udiff(1)=(MD%ucart(1,jatom,istep)-MD%ucart(1,iatom,istep))*MD%weights(istep)
344 : 7957392 : udiff(2)=(MD%ucart(2,jatom,istep)-MD%ucart(2,iatom,istep))*MD%weights(istep)
345 : 7957392 : udiff(3)=(MD%ucart(3,jatom,istep)-MD%ucart(3,iatom,istep))*MD%weights(istep)
346 : :
347 : : ! F_i^\mu(t)=\sum_{\alpha\beta,j,\nu}S^{\mu\alpha}.S^{\nu\beta}.\Phi_{ij}^{\alpha\beta}.u_j^\nu(t)
348 : 7957392 : SSu(:,:)=zero
349 [ + + ]: 31829568 : do nu=1,3
350 [ + + + + ]: 891227904 : SSu(:,:)=SSu(:,:)+SS_ref(:,:,nu,isym,itrans)*udiff(nu)
351 : : end do
352 [ + + ]: 32235198 : do mu=1,3
353 [ + + ]: 134399328 : do icoeff=1,ncoeff
354 [ + + ]: 1025697600 : terme=sum(SSu(mu,:)*Shell2at%proj(:,icoeff,ishell))
355 : : !FB write(Invar%stdlog,*) 'indices=', mu+3*(iatom-1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev
356 : : Solver%fcoeff(mu+3*(iatom-1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev)= &
357 : 126441936 : & Solver%fcoeff(mu+3*(iatom-1)+3*Invar%natom*(istep-1),icoeff+ncoeff_prev)+terme
358 : : end do
359 : : end do
360 : :
361 : : end do !istep
362 : : end do !iatshell
363 : : end do !iatom
364 : : end do !ishell
365 : 43 : write(Invar%stdout,*) ' ------- achieved'
366 [ - + ]: 43 : ABI_FREE(SS_ref)
367 : :
368 : 43 : end subroutine tdep_calc_phi2fcoeff
369 : :
370 : : !====================================================================================================
371 : :
372 : 14 : subroutine tdep_calc_phi3fcoeff(Solver,Invar,Shell3at,Sym,MD)
373 : :
374 : : type(tdep_Solver_type), intent(inout) :: Solver
375 : : type(atdep_dataset_type),intent(in) :: Invar
376 : : type(Shell_type),intent(in) :: Shell3at
377 : : type(Symmetries_type),intent(in) :: Sym
378 : : type(tdep_Sampling_type), intent(in) :: MD
379 : :
380 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,jatom,katom
381 : : integer :: icoeff,isym,itrans,iatshell
382 : : integer :: mu,nu,xi,alpha,beta,gama,iindex
383 : : double precision :: temp
384 : : double precision :: udiff_ki(3),udiff_ji(3)
385 : 14 : double precision, allocatable :: SSS_proj(:,:,:,:)
386 : : double precision :: SSS_tmp(27), proj_tmp(27)
387 : 14 : type(Constraints_type) :: Const
388 : :
389 [ - + + - : 3102 : ABI_MALLOC(Const%Sprod,(Sym%nsym,6))
+ - + - -
+ - + + +
+ + ]
390 [ + + ]: 510 : do isym=1,Sym%nsym
391 [ + + ]: 3486 : do itrans=1,6
392 [ - + - + : 2931856 : ABI_MALLOC(Const%Sprod(isym,itrans)%SSS,(3,27,3,3)); Const%Sprod(isym,itrans)%SSS(:,:,:,:)=zero
+ + + + +
+ + + ]
393 : : end do
394 : : end do
395 : :
396 : : ! For each couple of atoms, transform the Phi3 (3x3x3) ifc matrix using the symetry operation (S)
397 : : ! Note: iatom=1 is excluded in order to take into account the atomic sum rule (see below)
398 [ + + ]: 510 : do isym=1,Sym%nsym
399 [ + + ]: 1998 : do mu=1,3
400 [ + + ]: 6448 : do alpha=1,3
401 [ + + ]: 19344 : do nu=1,3
402 [ + + ]: 58032 : do beta=1,3
403 [ + + ]: 174096 : do xi=1,3
404 [ + + ]: 522288 : do gama=1,3
405 : 361584 : temp=Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu,beta,isym,1)*Sym%S_ref(xi,gama,isym,1)
406 : 361584 : Const%Sprod(isym,1)%SSS(mu,gama+(beta-1)*3+(alpha-1)*9,nu,xi)=temp !\Phi3_efg
407 : 361584 : Const%Sprod(isym,2)%SSS(mu,gama+(beta-1)*3+(alpha-1)*9,xi,nu)=temp !\Phi3_egf
408 : 361584 : Const%Sprod(isym,3)%SSS(nu,gama+(beta-1)*3+(alpha-1)*9,mu,xi)=temp !\Phi3_feg
409 : 361584 : Const%Sprod(isym,4)%SSS(nu,gama+(beta-1)*3+(alpha-1)*9,xi,mu)=temp !\Phi3_fge
410 : 361584 : Const%Sprod(isym,5)%SSS(xi,gama+(beta-1)*3+(alpha-1)*9,mu,nu)=temp !\Phi3_gef
411 : 482112 : Const%Sprod(isym,6)%SSS(xi,gama+(beta-1)*3+(alpha-1)*9,nu,mu)=temp !\Phi3_gfe
412 : : end do
413 : : end do
414 : : end do
415 : : end do
416 : : end do
417 : : end do
418 : : end do
419 : :
420 : 14 : write(Invar%stdout,*) ' Compute the coefficients (at the 3rd order) used in the Moore-Penrose...'
421 [ + + ]: 63 : do ishell=1,Shell3at%nshell
422 [ + + ]: 8391 : do iatom=1,Invar%natom
423 [ + + ]: 8328 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
424 [ + + ]: 86485 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
425 : 79388 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
426 : 79388 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
427 : : !FB if (iatom==jatom.or.iatom==katom) cycle
428 : 79388 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
429 : 79388 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
430 : 79388 : ncoeff =Shell3at%ncoeff(ishell)
431 : 79388 : ncoeff_prev=Shell3at%ncoeff_prev(ishell)+Solver%ncoeff2nd+Solver%ncoeff1st
432 : :
433 [ + + + - : 20913204 : ABI_MALLOC(SSS_proj,(3,3,3,ncoeff)) ; SSS_proj(:,:,:,:)=zero
+ - + + -
+ - + + +
+ + + + +
+ ]
434 [ + + ]: 317552 : do mu=1,3
435 [ + + ]: 1032044 : do nu=1,3
436 [ + + ]: 3096132 : do xi=1,3
437 [ + + ]: 60017328 : SSS_tmp(:)=Const%Sprod(isym,itrans)%SSS(mu,:,nu,xi)
438 [ + + ]: 16813620 : do icoeff=1,ncoeff
439 [ + + ]: 390758256 : proj_tmp(:)=Shell3at%proj(:,icoeff,ishell)
440 : 16099128 : SSS_proj(mu,nu,xi,icoeff)=DDOT(27,SSS_tmp,1,proj_tmp,1)
441 : : end do
442 : : end do
443 : : end do
444 : : end do
445 [ + + ]: 1667148 : do istep=1,Invar%my_nstep
446 : 1587760 : iindex=3*(iatom-1)+3*Invar%natom*(istep-1)
447 : : ! In order to impose the acoustic sum rule we use :
448 : : !FB udiff_ji(:)=MD%ucart(:,jatom,istep)-MD%ucart(:,iatom,istep)
449 : : !FB udiff_ki(:)=MD%ucart(:,katom,istep)-MD%ucart(:,iatom,istep)
450 [ + + ]: 6351040 : udiff_ji(:)=MD%ucart(:,jatom,istep)
451 [ + + ]: 6351040 : udiff_ki(:)=MD%ucart(:,katom,istep)
452 : : ! F_i^{\mu}(t)=\sum_{\alpha\beta\gamma,jk,\nu\xi} S^{\mu\alpha}.S^{\nu\beta}.S^{\xi\gamma}.
453 : : ! \Phi3_{ijk}^{\alpha\beta\gamma}.udiff_k^\xi(t).udiff_j^\nu(t)
454 [ + + ]: 6430428 : do nu=1,3
455 [ + + ]: 20640880 : do xi=1,3
456 : : Solver%fcoeff(iindex+1:iindex+3,ncoeff_prev+1:ncoeff_prev+ncoeff)= &
457 : : & Solver%fcoeff(iindex+1:iindex+3,ncoeff_prev+1:ncoeff_prev+ncoeff)+&
458 [ + + + + ]: 391203840 : & SSS_proj(1:3,nu,xi,1:ncoeff)*udiff_ji(nu)*udiff_ki(xi)/2.d0*MD%weights(istep)
459 : : end do
460 : : end do
461 : : end do !istep
462 [ - + ]: 87716 : ABI_FREE(SSS_proj)
463 : : end do !iatshell
464 : : end do !iatom
465 : : end do !ishell
466 : 14 : write(Invar%stdout,*) ' ------- achieved'
467 [ + + ]: 510 : do isym=1,Sym%nsym
468 [ + + ]: 3486 : do itrans=1,6
469 [ - + ]: 3472 : ABI_FREE(Const%Sprod(isym,itrans)%SSS)
470 : : end do
471 : : end do
472 [ + - + + : 3004 : ABI_FREE(Const%Sprod)
- + - + -
+ - + ]
473 : :
474 : 14 : end subroutine tdep_calc_phi3fcoeff
475 : :
476 : : !====================================================================================================
477 : :
478 : 8 : subroutine tdep_calc_phi4fcoeff(Solver,Invar,Shell4at,Sym,MD)
479 : :
480 : : type(tdep_Solver_type), intent(inout) :: Solver
481 : : type(atdep_dataset_type),intent(in) :: Invar
482 : : type(Shell_type),intent(in) :: Shell4at
483 : : type(Symmetries_type),intent(in) :: Sym
484 : : type(tdep_Sampling_type), intent(inout) :: MD
485 : :
486 : : integer :: ishell,ncoeff,ncoeff_prev,istep,iatom,jatom,katom,latom
487 : : integer :: icoeff,isym,iatshell,itrans,counter
488 : : integer :: mu,nu,xi,zeta,alpha,beta,gama,delta,iindex_l,iindex_h
489 : : integer :: ncoeff_prev_l,ncoeff_prev_h
490 : : double precision :: temp,SSSS_tmp(81),proj_tmp(81)
491 : 8 : double precision, allocatable :: SSSS_proj(:,:,:,:,:)
492 : 8 : type(Constraints_type) :: Const
493 : :
494 [ - + + - : 7032 : ABI_MALLOC(Const%Sprod,(Sym%nsym,24))
+ - + - -
+ - + + +
+ + ]
495 [ + + ]: 292 : do isym=1,Sym%nsym
496 [ + + ]: 7108 : do itrans=1,24
497 [ - + - + : 59899292 : ABI_MALLOC(Const%Sprod(isym,itrans)%SSSS,(3,81,3,3,3)); Const%Sprod(isym,itrans)%SSSS(:,:,:,:,:)=zero
+ + + + +
+ + + +
+ ]
498 : : end do
499 : : end do
500 : :
501 : : ! For each couple of atoms, transform the Phi4 (3x3x3) ifc matrix using the symetry operation (S)
502 : : ! Note: iatom=1 is excluded in order to take into account the atomic sum rule (see below)
503 [ + + ]: 292 : do isym=1,Sym%nsym
504 [ + + ]: 1144 : do mu=1,3
505 [ + + ]: 3692 : do alpha=1,3
506 [ + + ]: 11076 : do nu=1,3
507 [ + + ]: 33228 : do beta=1,3
508 [ + + ]: 99684 : do xi=1,3
509 [ + + ]: 299052 : do gama=1,3
510 [ + + ]: 897156 : do zeta=1,3
511 [ + + ]: 2691468 : do delta=1,3
512 : : #if defined FC_NVHPC
513 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
514 : : #endif
515 : :
516 : 1863324 : counter=delta+(gama-1)*3+(beta-1)*9+(alpha-1)*27
517 : : temp=Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu ,beta ,isym,1)*&
518 : 1863324 : & Sym%S_ref(xi,gama ,isym,1)*Sym%S_ref(zeta,delta,isym,1)
519 : 1863324 : Const%Sprod(isym,1 )%SSSS(mu,counter,nu,xi,zeta)=temp !\Phi4_efgh
520 : 1863324 : Const%Sprod(isym,2 )%SSSS(mu,counter,xi,nu,zeta)=temp !\Phi4_egfh
521 : 1863324 : Const%Sprod(isym,3 )%SSSS(nu,counter,mu,xi,zeta)=temp !\Phi4_fegh
522 : 1863324 : Const%Sprod(isym,4 )%SSSS(nu,counter,xi,mu,zeta)=temp !\Phi4_fgeh
523 : 1863324 : Const%Sprod(isym,5 )%SSSS(xi,counter,mu,nu,zeta)=temp !\Phi4_gefh
524 : 1863324 : Const%Sprod(isym,6 )%SSSS(xi,counter,nu,mu,zeta)=temp !\Phi4_gfeh
525 : :
526 : 1863324 : Const%Sprod(isym,7 )%SSSS(mu,counter,nu,zeta,xi)=temp !\Phi4_efhg
527 : 1863324 : Const%Sprod(isym,8 )%SSSS(mu,counter,xi,zeta,nu)=temp !\Phi4_eghf
528 : 1863324 : Const%Sprod(isym,9 )%SSSS(nu,counter,mu,zeta,xi)=temp !\Phi4_fehg
529 : 1863324 : Const%Sprod(isym,10)%SSSS(nu,counter,xi,zeta,mu)=temp !\Phi4_fghe
530 : 1863324 : Const%Sprod(isym,11)%SSSS(xi,counter,mu,zeta,nu)=temp !\Phi4_gehf
531 : 1863324 : Const%Sprod(isym,12)%SSSS(xi,counter,nu,zeta,mu)=temp !\Phi4_gfhe
532 : :
533 : 1863324 : Const%Sprod(isym,13)%SSSS(mu,counter,zeta,nu,xi)=temp !\Phi4_ehfg
534 : 1863324 : Const%Sprod(isym,14)%SSSS(mu,counter,zeta,xi,nu)=temp !\Phi4_ehgf
535 : 1863324 : Const%Sprod(isym,15)%SSSS(nu,counter,zeta,mu,xi)=temp !\Phi4_fheg
536 : 1863324 : Const%Sprod(isym,16)%SSSS(nu,counter,zeta,xi,mu)=temp !\Phi4_fhge
537 : 1863324 : Const%Sprod(isym,17)%SSSS(xi,counter,zeta,mu,nu)=temp !\Phi4_ghef
538 : 1863324 : Const%Sprod(isym,18)%SSSS(xi,counter,zeta,nu,mu)=temp !\Phi4_ghfe
539 : :
540 : 1863324 : Const%Sprod(isym,19)%SSSS(zeta,counter,mu,nu,xi)=temp !\Phi4_hefg
541 : 1863324 : Const%Sprod(isym,20)%SSSS(zeta,counter,mu,xi,nu)=temp !\Phi4_hegf
542 : 1863324 : Const%Sprod(isym,21)%SSSS(zeta,counter,nu,mu,xi)=temp !\Phi4_hfeg
543 : 1863324 : Const%Sprod(isym,22)%SSSS(zeta,counter,nu,xi,mu)=temp !\Phi4_hfge
544 : 1863324 : Const%Sprod(isym,23)%SSSS(zeta,counter,xi,mu,nu)=temp !\Phi4_hgef
545 : 2484432 : Const%Sprod(isym,24)%SSSS(zeta,counter,xi,nu,mu)=temp !\Phi4_hgfe
546 : :
547 : : end do
548 : : end do
549 : : end do
550 : : end do
551 : : end do
552 : : end do
553 : : end do
554 : : end do
555 : : end do
556 : :
557 : 8 : write(Invar%stdout,*) ' Compute the coefficients (at the 4th order) used in the Moore-Penrose...'
558 [ + + ]: 54 : do ishell=1,Shell4at%nshell
559 [ + + ]: 7886 : do iatom=1,Invar%natom
560 [ + + ]: 7832 : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
561 [ + + ]: 250586 : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
562 : 243732 : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
563 : 243732 : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
564 : 243732 : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
565 : 243732 : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
566 : 243732 : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
567 : 243732 : ncoeff =Shell4at%ncoeff(ishell)
568 : 243732 : ncoeff_prev=Shell4at%ncoeff_prev(ishell)+Solver%ncoeff3rd+Solver%ncoeff2nd+Solver%ncoeff1st
569 : 243732 : ncoeff_prev_l=ncoeff_prev+1
570 : 243732 : ncoeff_prev_h=ncoeff_prev+ncoeff
571 : : #if defined FC_NVHPC
572 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
573 : : #endif
574 [ - + + - : 699105928 : ABI_MALLOC(SSSS_proj,(3,3,3,3,ncoeff)) ; SSSS_proj(:,:,:,:,:)=zero
+ - + - -
+ - + + +
+ + + + +
+ + + ]
575 [ + + ]: 974928 : do mu=1,3
576 [ + + ]: 3168516 : do nu=1,3
577 [ + + ]: 9505548 : do xi=1,3
578 [ + + ]: 28516644 : do zeta=1,3
579 [ + + ]: 1618867944 : SSSS_tmp(:)=Const%Sprod(isym,itrans)%SSSS(mu,:,nu,xi,zeta)
580 [ + + ]: 493830108 : do icoeff=1,ncoeff
581 [ + + ]: 38335578264 : proj_tmp(:)=Shell4at%proj(:,icoeff,ishell)
582 : : ! SSSS_proj(mu,nu,xi,zeta,icoeff)=DDOT(81,Const%Sprod(isym,itrans)%SSSS(mu,:,nu,xi,zeta),1,proj(:,icoeff,ishell),1)
583 : 487249344 : SSSS_proj(mu,nu,xi,zeta,icoeff)=DDOT(81,SSSS_tmp,1,proj_tmp,1)
584 : : end do
585 : : end do
586 : : end do
587 : : end do
588 : : end do
589 [ + + ]: 5118372 : do istep=1,Invar%my_nstep
590 : 4874640 : iindex_l=3*(iatom-1)+3*Invar%natom*(istep-1)+1
591 : 4874640 : iindex_h=3*(iatom-1)+3*Invar%natom*(istep-1)+3
592 : : ! 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}.
593 : : ! \Phi4_{ijkl}^{\alpha\beta\gamma\delta}.u_l^\zeta(t).u_k^\xi(t).u_j^\nu(t)
594 [ + + ]: 19742292 : do nu=1,3
595 [ + + ]: 63370320 : do xi=1,3
596 [ + + ]: 190110960 : do zeta=1,3
597 : : Solver%fcoeff(iindex_l:iindex_h,ncoeff_prev_l:ncoeff_prev_h)= &
598 : : & Solver%fcoeff(iindex_l:iindex_h,ncoeff_prev_l:ncoeff_prev_h)+&
599 : : & 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 *&
600 [ + + + + ]: 12642341760 : & MD%weights(istep)
601 : : end do
602 : : end do
603 : : end do
604 : : end do !istep
605 [ - + ]: 251564 : ABI_FREE(SSSS_proj)
606 : : end do !iatshell
607 : : end do !iatom
608 : : end do !ishell
609 : 8 : write(Invar%stdout,*) ' ------- achieved'
610 [ + + ]: 292 : do isym=1,Sym%nsym
611 [ + + ]: 7108 : do itrans=1,24
612 [ - + ]: 7100 : ABI_FREE(Const%Sprod(isym,itrans)%SSSS)
613 : : end do
614 : : end do
615 [ + - + + : 6832 : ABI_FREE(Const%Sprod)
- + - + -
+ - + ]
616 : :
617 : 8 : end subroutine tdep_calc_phi4fcoeff
618 : :
619 : : !=====================================================================================================
620 : :
621 : 45 : subroutine tdep_calc_MoorePenrose(Solver,simult,Invar,MPIdata)
622 : :
623 : : type(tdep_Solver_type), intent(inout) :: Solver
624 : : type(atdep_dataset_type),intent(in) :: Invar
625 : : type(MPI_enreg_type), intent(in) :: MPIdata
626 : : integer, intent(in) :: simult
627 : :
628 : : integer :: INFO,ntotcoeff,ntotconst
629 : : integer :: natnstep,nconcoef,ierr,ncoeff_prev,nconst_prev,iconst,icoeff
630 : 45 : integer, allocatable :: IPIV(:)
631 : 45 : double precision, allocatable :: WORK(:)
632 : 45 : double precision, allocatable :: ffcoeff_tmp(:,:),fforces_tmp(:),b_const(:)
633 : 45 : double precision, allocatable :: A_tot(:,:),A_inv(:,:),b_tot(:),x_tot(:)
634 : :
635 : 45 : write(Invar%stdout,*) '################### And compute the pseudo-inverse ##########################'
636 : 45 : write(Invar%stdout,*) '#############################################################################'
637 : :
638 : 45 : natnstep = 3 * Solver%natom * Solver%my_nstep
639 : :
640 [ + + ]: 45 : if (simult.eq.0) then
641 : : ! Simultaneously (Invar%together=1)
642 : 42 : ncoeff_prev=0
643 : 42 : nconst_prev=0
644 : 42 : ntotcoeff=Solver%ntotcoeff
645 : 42 : ntotconst=Solver%ntotconst
646 [ + + ]: 3 : else if (simult.eq.1) then
647 : : ! Successively (Invar%together=0 and Invar%order=2)
648 : 1 : ncoeff_prev=0
649 : 1 : nconst_prev=0
650 : 1 : ntotcoeff=Solver%ncoeff1st +Solver%ncoeff2nd
651 : 1 : ntotconst=Solver%nconst_1st+Solver%nconst_2nd
652 [ + + ]: 2 : else if (simult.eq.2) then
653 : : ! Successively (Invar%together=0 and Invar%order=3)
654 : 1 : ncoeff_prev=Solver%ncoeff1st +Solver%ncoeff2nd
655 : 1 : nconst_prev=Solver%nconst_1st+Solver%nconst_2nd
656 : 1 : ntotcoeff=Solver%ncoeff3rd
657 : 1 : ntotconst=Solver%nconst_3rd
658 [ + - ]: 1 : else if (simult.eq.3) then
659 : : ! Successively (Invar%together=0 and Invar%order=4)
660 : 1 : ncoeff_prev=Solver%ncoeff1st +Solver%ncoeff2nd +Solver%ncoeff3rd
661 : 1 : nconst_prev=Solver%nconst_1st+Solver%nconst_2nd+Solver%nconst_3rd
662 : 1 : ntotcoeff=Solver%ncoeff4th
663 : 1 : ntotconst=Solver%nconst_4th
664 : : end if
665 : 45 : nconcoef=ntotcoeff+ntotconst
666 [ + + + + ]: 45 : if ((ntotconst.gt.0).and.(simult.ge.2)) then
667 [ - + + - : 23 : ABI_MALLOC(b_const,(ntotconst)) ; b_const(:)=0.d0
+ - - + -
+ + + ]
668 [ + + ]: 21 : do iconst=1,ntotconst
669 [ + + ]: 521 : do icoeff=1,ncoeff_prev
670 : : b_const(iconst)=b_const(iconst)+&
671 : 520 : & Solver%const(nconst_prev+iconst,icoeff)*Solver%theta(icoeff)
672 : : end do
673 : : end do
674 : : end if
675 : :
676 [ - + + - : 461308 : ABI_CALLOC(ffcoeff_tmp,(ntotcoeff,ntotcoeff))
- + + - +
- + - - +
- + + + +
+ ]
677 [ - + + - : 2780 : ABI_CALLOC(fforces_tmp,(ntotcoeff))
+ - + - -
+ - + +
+ ]
678 [ - + + - : 653574 : ABI_CALLOC(A_tot,(nconcoef,nconcoef))
- + + - +
- + - - +
- + + + +
+ ]
679 [ - + + - : 653574 : ABI_CALLOC(A_inv,(nconcoef,nconcoef))
- + + - -
+ + - + +
+ + ]
680 [ - + + - : 3194 : ABI_CALLOC(b_tot,(nconcoef))
+ - + - -
+ + - +
+ ]
681 [ - + + - : 3149 : ABI_CALLOC(x_tot,(nconcoef))
- + + - +
+ ]
682 : : call DGEMM('T','N',ntotcoeff,ntotcoeff,natnstep,2.d0,&
683 : : & Solver%fcoeff(:,ncoeff_prev+1:ncoeff_prev+ntotcoeff),natnstep,&
684 : : & Solver%fcoeff(:,ncoeff_prev+1:ncoeff_prev+ntotcoeff),natnstep,&
685 : 45 : & 0.d0,ffcoeff_tmp,ntotcoeff)
686 : : ! NOTE, we have to solve F_ij = -\sum_j \Phi_ij u_j, so we add a minus sign
687 : : call DGEMV('T',natnstep,ntotcoeff,-2.d0,&
688 : : & Solver%fcoeff(:,ncoeff_prev+1:ncoeff_prev+ntotcoeff),natnstep,&
689 : 45 : & Solver%Forces,1,0.d0,fforces_tmp,1)
690 : 45 : call xmpi_sum(ffcoeff_tmp,MPIdata%comm_step,ierr)
691 : 45 : call xmpi_sum(fforces_tmp,MPIdata%comm_step,ierr)
692 : :
693 [ + + + + ]: 461173 : A_tot(1:ntotcoeff,1:ntotcoeff)=ffcoeff_tmp(1:ntotcoeff,1:ntotcoeff)
694 [ - + ]: 45 : ABI_FREE(ffcoeff_tmp)
695 [ + + ]: 45 : if (ntotconst.gt.0) then
696 : : A_tot(ntotcoeff+1:nconcoef,1:ntotcoeff)=&
697 [ + + + + ]: 88222 : & Solver%const(nconst_prev+1:nconst_prev+ntotconst,ncoeff_prev+1:ncoeff_prev+ntotcoeff)
698 : : A_tot(1:ntotcoeff,ntotcoeff+1:nconcoef)=&
699 [ + + + + ]: 86382 : & transpose(Solver%const(nconst_prev+1:nconst_prev+ntotconst,ncoeff_prev+1:ncoeff_prev+ntotcoeff))
700 : : !FB ABI_FREE(Solver%const)
701 : : end if
702 [ + + ]: 2690 : b_tot(1:ntotcoeff)=fforces_tmp(:)
703 [ + + + + ]: 45 : if ((ntotconst.gt.0).and.(simult.ge.2)) then
704 [ + + ]: 21 : b_tot(ntotcoeff+1:nconcoef)=-b_const(1:ntotconst)
705 [ - + ]: 1 : ABI_FREE(b_const)
706 : : end if
707 [ - + ]: 45 : ABI_FREE(fforces_tmp)
708 : :
709 [ + - - + : 15205 : ABI_MALLOC(WORK, (5 * nconcoef)); WORK(:) = 0.d0
+ - + - +
- - + - +
+ + ]
710 [ - + + - : 3194 : ABI_MALLOC(IPIV, (nconcoef)); IPIV(:) = 0
+ - + - -
+ + - +
+ ]
711 [ + + + + ]: 653439 : A_inv(:,:) = A_tot(:,:)
712 : : !BEGIN DEBUG
713 : : !write(Invar%stdout,*) ' '
714 : : !write(Invar%stdout,*) ' The matrix A_inv is (before DGETRF):'
715 : : !do icoeff=1,nconcoef
716 : : ! write(Invar%stdout,*) (A_inv(icoeff,iconst), iconst=1, nconcoef)
717 : : !end do
718 : : !END DEBUG
719 : :
720 : : ! Check for small pivot elements
721 [ + + ]: 3059 : do icoeff=1,nconcoef
722 [ + + ]: 3059 : if (abs(A_inv(icoeff, icoeff)) < tol12) then
723 : 369 : write(Invar%stdlog,*) ' WARNING: Small pivot value at index ', icoeff, ' : ', A_inv(icoeff, icoeff)
724 : : ! A_inv(icoeff, icoeff) = tol14 ! Regularization to avoid numerical issues
725 [ + + + + : 50146203 : A_inv(icoeff, icoeff) = max(EPSILON(1.0_dp) * maxval(abs(A_inv)), tol12)
+ + + - -
+ - - ]
726 : : end if
727 : : end do
728 : :
729 : : ! Perform LU factorization
730 : 45 : call DGETRF(nconcoef, nconcoef, A_inv, nconcoef, IPIV, INFO)
731 [ - + ]: 45 : if (INFO.ne.0) then
732 : 0 : write(Invar%stdout,*) 'ERROR: Singular matrix detected in DGETRF. INFO=', INFO
733 : 0 : stop
734 : : end if
735 : :
736 : : ! Check for small pivot elements
737 [ + + ]: 3059 : do icoeff=1,nconcoef
738 [ - + ]: 3059 : if (abs(A_inv(icoeff, icoeff)) < tol12) then
739 : 0 : write(Invar%stdlog,*) ' WARNING: Small pivot value at index ', icoeff, ' : ', A_inv(icoeff, icoeff)
740 : : ! A_inv(icoeff, icoeff) = tol14 ! Regularization to avoid numerical issues
741 [ # # # # : 0 : A_inv(icoeff, icoeff) = max(EPSILON(1.0_dp) * maxval(abs(A_inv)), tol12)
# # # # #
# # # ]
742 : : end if
743 : : end do
744 : :
745 : : ! Compute matrix inverse using LU decomposition
746 : 45 : call DGETRI(nconcoef, A_inv, nconcoef, IPIV, WORK, 5 * nconcoef, INFO)
747 [ - + ]: 45 : if (INFO.ne.0) then
748 : 0 : write(Invar%stdout,*) 'ERROR: Matrix inversion failed in DGETRI. INFO=', INFO
749 : 0 : stop
750 : : end if
751 : : ! BEGIN DEBUG
752 : : !write(Invar%stdout,*) ' '
753 : : !write(Invar%stdout,*) ' The inverse matrix is (after DGETRI):'
754 : : !do icoeff=1,nconcoef
755 : : ! write(Invar%stdout,*) (A_inv(icoeff,iconst), iconst=1, nconcoef)
756 : : !end do
757 : : ! END DEBUG
758 : :
759 [ - + ]: 45 : ABI_FREE(WORK)
760 [ - + ]: 45 : ABI_FREE(IPIV)
761 : :
762 : 45 : call DGEMV('N',nconcoef,nconcoef,1.d0,A_inv,nconcoef,b_tot,1,0.d0,x_tot,1)
763 : 45 : write(Invar%stdout,*) ' The problem is solved'
764 : 45 : write(Invar%stdout,*) ' '
765 : : !BEGIN DEBUG
766 : : !write(Invar%stdout,*) ' The solutions are:'
767 : : !do icoeff=1,nconcoef
768 : : ! write(Invar%stdout,'(1x,i4,1x,f15.10)') icoeff,x_tot(icoeff)
769 : : !end do
770 : : !write(Invar%stdout,'(a,1x,f15.10)')' condition number=',maxval(x_tot(:))/minval(x_tot(:))
771 : : !END DEBUG
772 : :
773 [ + + ]: 2690 : Solver%theta(ncoeff_prev+1:ncoeff_prev+ntotcoeff)=x_tot(1:ntotcoeff)
774 [ - + ]: 45 : ABI_FREE(A_tot)
775 [ - + ]: 45 : ABI_FREE(A_inv)
776 [ - + ]: 45 : ABI_FREE(b_tot)
777 [ - + ]: 45 : ABI_FREE(x_tot)
778 : :
779 [ - + ]: 90 : end subroutine tdep_calc_MoorePenrose
780 : :
781 : :
782 : : !====================================================================================================
783 : :
784 : 44 : subroutine tdep_calc_constraints(Solver,distance,Invar,MPIdata,Sym,&
785 : : & Shell1at,Shell2at,Shell3at,Shell4at)
786 : :
787 : : type(tdep_Solver_type), intent(inout) :: Solver
788 : : type(atdep_dataset_type),intent(in) :: Invar
789 : : type(Symmetries_type),intent(in) :: Sym
790 : : type(MPI_enreg_type), intent(in) :: MPIdata
791 : : type(Shell_type),intent(in) :: Shell1at
792 : : type(Shell_type),intent(in) :: Shell2at
793 : : type(Shell_type),optional,intent(in) :: Shell3at
794 : : type(Shell_type),optional,intent(in) :: Shell4at
795 : : double precision, intent(in) :: distance(Invar%natom,Invar%natom,4)
796 : :
797 : : integer :: ishell,ncoeff,ncoeff_prev,iatom,jatom,katom,latom,iatshell,counter
798 : : integer :: icoeff,iconst,nconst_loc,iconst_loc,iconst_new,isym,itrans,ntotcoeff,iat_mod
799 : : integer :: mu,nu,xi,zeta,alpha,beta,gama,delta,lambda,natom_unitcell,natom,ii
800 : : double precision :: terme,temp,terme1,terme2,terme3,terme4
801 : 44 : double precision, allocatable :: SS_ref(:,:,:,:,:)
802 : 44 : double precision, allocatable :: vect(:,:)
803 : 44 : double precision, allocatable :: const_rot1st(:,:,:)
804 : 44 : double precision, allocatable :: const_rot2nd(:,:,:,:,:)
805 : 44 : double precision, allocatable :: const_dynmat(:,:,:,:,:)
806 : 44 : double precision, allocatable :: const_huang(:,:,:,:,:)
807 : : !FB double precision, allocatable :: const_asr4th(:,:,:,:,:,:)
808 : : !FB double precision, allocatable :: const_rot4th(:,:,:,:,:,:,:)
809 : 44 : type(Constraints_type) :: Const3,Const4
810 : : logical :: order2,order3,order4
811 : :
812 : : !TODO Move parts of this routine into m_tdep_constraints
813 : :
814 : 44 : natom_unitcell=Invar%natom_unitcell
815 : 44 : natom =Invar%natom
816 : :
817 : 44 : order2 = .false.
818 : 44 : order3 = .false.
819 : 44 : order4 = .false.
820 [ + - ]: 44 : if (Invar%order.ge.2) order2=.true.
821 [ + + ]: 44 : if (Invar%order.ge.3) order3=.true.
822 : : !FB4th if (Invar%order.ge.4) order4=.true.
823 : : if (Invar%order.ge.4) order4=.false.
824 : :
825 : 44 : write(Invar%stdout,*) ' '
826 : 44 : write(Invar%stdout,*) '#############################################################################'
827 : 44 : write(Invar%stdout,*) '###################### Compute the constraints ##############################'
828 : :
829 : : ! For each couple of atoms, transform the Phi2 (3x3) ifc matrix using the symetry operation (S)
830 [ + - ]: 44 : if (order2.or.order3) then
831 [ - + + - : 363592 : ABI_CALLOC(SS_ref,(3,9,3,Sym%nsym,2))
+ - + - +
- - + - +
+ + + + +
+ + + +
+ ]
832 [ + + ]: 1666 : do isym=1,Sym%nsym
833 [ + + ]: 6532 : do alpha=1,3
834 [ + + ]: 21086 : do mu=1,3
835 [ + + ]: 63258 : do beta=1,3
836 [ + + ]: 189774 : do nu=1,3
837 : 131382 : temp=Sym%S_ref(alpha,mu,isym,1)*Sym%S_ref(beta,nu,isym,1)
838 : 131382 : SS_ref(alpha,nu+(mu-1)*3,beta,isym,1)=temp
839 : 175176 : SS_ref(alpha,mu+(nu-1)*3,beta,isym,2)=temp
840 : : end do
841 : : end do
842 : : end do
843 : : end do
844 : : end do
845 : : end if
846 : : ! For each couple of atoms, transform the Phi3 (3x3x3) ifc matrix using the symetry operation (S)
847 [ + + ]: 44 : if (order3.or.order4) then
848 [ - + + - : 3102 : ABI_MALLOC(Const3%Sprod,(Sym%nsym,6))
+ - + - -
+ - + + +
+ + ]
849 [ + + ]: 510 : do isym=1,Sym%nsym
850 [ + + ]: 3486 : do itrans=1,6
851 [ - + - + : 2931856 : ABI_CALLOC(Const3%Sprod(isym,itrans)%SSS,(3,27,3,3))
+ + + + +
+ + + ]
852 : : end do
853 : : end do
854 [ + + ]: 510 : do isym=1,Sym%nsym
855 [ + + ]: 1998 : do alpha=1,3
856 [ + + ]: 6448 : do mu=1,3
857 [ + + ]: 19344 : do beta=1,3
858 [ + + ]: 58032 : do nu=1,3
859 [ + + ]: 174096 : do gama=1,3
860 [ + + ]: 522288 : do xi=1,3
861 : 361584 : temp=Sym%S_ref(alpha,mu,isym,1)*Sym%S_ref(beta,nu,isym,1)*Sym%S_ref(gama,xi,isym,1)
862 : 361584 : Const3%Sprod(isym,1)%SSS(alpha,xi+(nu-1)*3+(mu-1)*9,beta ,gama) =temp !\Phi3_efg
863 : 361584 : Const3%Sprod(isym,2)%SSS(alpha,xi+(nu-1)*3+(mu-1)*9,gama ,beta) =temp !\Phi3_egf
864 : 361584 : Const3%Sprod(isym,3)%SSS(beta ,xi+(nu-1)*3+(mu-1)*9,alpha,gama) =temp !\Phi3_feg
865 : 361584 : Const3%Sprod(isym,4)%SSS(beta ,xi+(nu-1)*3+(mu-1)*9,gama ,alpha)=temp !\Phi3_fge
866 : 361584 : Const3%Sprod(isym,5)%SSS(gama ,xi+(nu-1)*3+(mu-1)*9,alpha,beta) =temp !\Phi3_gef
867 : 482112 : Const3%Sprod(isym,6)%SSS(gama ,xi+(nu-1)*3+(mu-1)*9,beta ,alpha)=temp !\Phi3_gfe
868 : : end do
869 : : end do
870 : : end do
871 : : end do
872 : : end do
873 : : end do
874 : : end do
875 : : end if
876 : :
877 : : ! For each couple of atoms, transform the Phi4 (3x3x3x3) ifc matrix using the symetry operation (S)
878 : : if (order4) then
879 : : ABI_MALLOC(Const4%Sprod,(Sym%nsym,24))
880 : : do isym=1,Sym%nsym
881 : : do itrans=1,24
882 : : ABI_CALLOC(Const4%Sprod(isym,itrans)%SSSS,(3,81,3,3,3))
883 : : end do
884 : : end do
885 : : do isym=1,Sym%nsym
886 : : do alpha=1,3
887 : : do mu=1,3
888 : : do beta=1,3
889 : : do nu=1,3
890 : : do gama=1,3
891 : : do xi=1,3
892 : : do delta=1,3
893 : : do zeta=1,3
894 : : counter=zeta+(xi-1)*3+(nu-1)*9+(mu-1)*27
895 : : temp=Sym%S_ref(alpha,mu,isym,1)*Sym%S_ref(beta ,nu ,isym,1)*&
896 : : & Sym%S_ref(gama,xi ,isym,1)*Sym%S_ref(delta,zeta,isym,1)
897 : : Const4%Sprod(isym,1 )%SSSS(alpha,counter,beta,gama,delta)=temp !\Phi4_efgh
898 : : Const4%Sprod(isym,2 )%SSSS(alpha,counter,gama,beta,delta)=temp !\Phi4_egfh
899 : : Const4%Sprod(isym,3 )%SSSS(beta,counter,alpha,gama,delta)=temp !\Phi4_fegh
900 : : Const4%Sprod(isym,4 )%SSSS(beta,counter,gama,alpha,delta)=temp !\Phi4_fgeh
901 : : Const4%Sprod(isym,5 )%SSSS(gama,counter,alpha,beta,delta)=temp !\Phi4_gefh
902 : : Const4%Sprod(isym,6 )%SSSS(gama,counter,beta,alpha,delta)=temp !\Phi4_gfeh
903 : :
904 : : Const4%Sprod(isym,7 )%SSSS(alpha,counter,beta,delta,gama)=temp !\Phi4_efhg
905 : : Const4%Sprod(isym,8 )%SSSS(alpha,counter,gama,delta,beta)=temp !\Phi4_eghf
906 : : Const4%Sprod(isym,9 )%SSSS(beta,counter,alpha,delta,gama)=temp !\Phi4_fehg
907 : : Const4%Sprod(isym,10)%SSSS(beta,counter,gama,delta,alpha)=temp !\Phi4_fghe
908 : : Const4%Sprod(isym,11)%SSSS(gama,counter,alpha,delta,beta)=temp !\Phi4_gehf
909 : : Const4%Sprod(isym,12)%SSSS(gama,counter,beta,delta,alpha)=temp !\Phi4_gfhe
910 : :
911 : : Const4%Sprod(isym,13)%SSSS(alpha,counter,delta,beta,gama)=temp !\Phi4_ehfg
912 : : Const4%Sprod(isym,14)%SSSS(alpha,counter,delta,gama,beta)=temp !\Phi4_ehgf
913 : : Const4%Sprod(isym,15)%SSSS(beta,counter,delta,alpha,gama)=temp !\Phi4_fheg
914 : : Const4%Sprod(isym,16)%SSSS(beta,counter,delta,gama,alpha)=temp !\Phi4_fhge
915 : : Const4%Sprod(isym,17)%SSSS(gama,counter,delta,alpha,beta)=temp !\Phi4_ghef
916 : : Const4%Sprod(isym,18)%SSSS(gama,counter,delta,beta,alpha)=temp !\Phi4_ghfe
917 : :
918 : : Const4%Sprod(isym,19)%SSSS(delta,counter,alpha,beta,gama)=temp !\Phi4_hefg
919 : : Const4%Sprod(isym,20)%SSSS(delta,counter,alpha,gama,beta)=temp !\Phi4_hegf
920 : : Const4%Sprod(isym,21)%SSSS(delta,counter,beta,alpha,gama)=temp !\Phi4_hfeg
921 : : Const4%Sprod(isym,22)%SSSS(delta,counter,beta,gama,alpha)=temp !\Phi4_hfge
922 : : Const4%Sprod(isym,23)%SSSS(delta,counter,gama,alpha,beta)=temp !\Phi4_hgef
923 : : Const4%Sprod(isym,24)%SSSS(delta,counter,gama,beta,alpha)=temp !\Phi4_hgfe
924 : :
925 : : end do
926 : : end do
927 : : end do
928 : : end do
929 : : end do
930 : : end do
931 : : end do
932 : : end do
933 : : end do
934 : : end if
935 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
936 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Compute the constraints !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
937 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
938 : 44 : ntotcoeff=Solver%ntotcoeff
939 : : ! First order only
940 : 44 : write(Invar%stdout,*) '########################## At the 1st order #################################'
941 [ + - ]: 44 : if (order2) then
942 [ - + + - : 35765 : ABI_CALLOC(const_rot1st, (3,3,ntotcoeff))
+ - + - -
+ - + + +
+ + + + ]
943 [ - + + - : 541517 : ABI_CALLOC(const_rot2nd, (3,3,3,natom_unitcell,ntotcoeff))
- + + - +
- + - - +
- + + + +
+ + + + +
+ + ]
944 [ - + + - : 1855835 : ABI_CALLOC(const_dynmat, (3,3,natom_unitcell,natom_unitcell,ntotcoeff))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + +
+ ]
945 [ - + + - : 331793 : ABI_CALLOC(const_huang, (3,3,3,3,ntotcoeff))
+ - + - -
+ - + + +
+ + + + +
+ + + ]
946 [ + + ]: 107 : do ishell=1,Shell1at%nshell
947 [ - + ]: 63 : if (Shell1at%neighbours(1,ishell)%n_interactions.eq.0) cycle
948 [ - - ]: 44 : do iatshell=1,Shell1at%neighbours(1,ishell)%n_interactions
949 : 0 : iatom=Shell1at%neighbours(1,ishell)%atomj_in_shell(iatshell)
950 [ # # ]: 0 : if (iatom.ge.natom_unitcell) cycle
951 [ # # ]: 0 : if (iatom.eq.1) cycle
952 : 0 : isym=Shell1at%neighbours(1,ishell)%sym_in_shell(iatshell)
953 : 0 : ncoeff =Shell1at%ncoeff(ishell)
954 : 0 : ncoeff_prev=Shell1at%ncoeff_prev(ishell)
955 [ - - ]: 63 : do alpha=1,3
956 [ # # ]: 0 : do beta=1,3
957 [ # # ]: 0 : do icoeff=1,ncoeff
958 : : ! 1/ Rotational invariances (1st order)
959 [ # # ]: 0 : terme1=sum(Sym%S_ref(alpha,:,isym,1)*Shell1at%proj(:,icoeff,ishell))*distance(1,iatom,beta +1)
960 [ # # ]: 0 : terme2=sum(Sym%S_ref(beta ,:,isym,1)*Shell1at%proj(:,icoeff,ishell))*distance(1,iatom,alpha+1)
961 : : const_rot1st(alpha,beta,icoeff+ncoeff_prev)= &
962 : 0 : & const_rot1st(alpha,beta,icoeff+ncoeff_prev)+terme1-terme2
963 : :
964 : : ! 2/ Rotational invariances (for the 2nd order)
965 [ # # ]: 0 : do gama=1,3
966 : 0 : terme1=zero ; terme2=zero
967 [ # # # # ]: 0 : if (alpha.eq.gama) terme1=sum(Sym%S_ref(beta,:,isym,1)*Shell1at%proj(:,icoeff,ishell))
968 [ # # # # ]: 0 : if (alpha.eq.beta) terme2=sum(Sym%S_ref(gama,:,isym,1)*Shell1at%proj(:,icoeff,ishell))
969 : : const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)=&
970 : 0 : & const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)+terme1-terme2
971 : : const_rot2nd(alpha,beta,gama,1,icoeff+ncoeff_prev)=&
972 : 0 : & const_rot2nd(alpha,beta,gama,1,icoeff+ncoeff_prev)-terme1+terme2
973 : : end do
974 : : end do
975 : : end do
976 : : end do
977 : : end do !iatshell
978 : : end do !ishell
979 : :
980 : : ! First + second order
981 : 44 : write(Invar%stdout,*) '########################## At the 2nd order #################################'
982 [ + + ]: 544 : do ishell=1,Shell2at%nshell
983 [ + + ]: 2830 : do iatom=1,natom_unitcell
984 [ + + ]: 2286 : if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
985 [ + + ]: 6867 : do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
986 : 5268 : jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
987 [ + + ]: 5268 : if (iatom==jatom) cycle
988 : 5162 : isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
989 : 5162 : itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
990 : 5162 : ncoeff =Shell2at%ncoeff(ishell)
991 : 5162 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
992 : 5162 : iat_mod=mod(jatom+natom_unitcell-1,natom_unitcell)+1
993 : : ! 1/ Rotational invariances (2nd order). Number of constraints = natom_unitcell*3**2
994 [ + + ]: 20648 : do alpha=1,3
995 [ + + ]: 67106 : do beta=1,3
996 [ + + ]: 201318 : do gama=1,3
997 [ + + ]: 834156 : do icoeff=1,ncoeff
998 [ + + ]: 6483240 : terme1=sum(SS_ref(alpha,:,beta,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*distance(iatom,jatom,gama+1)
999 [ + + ]: 6483240 : terme2=sum(SS_ref(alpha,:,gama,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*distance(iatom,jatom,beta+1)
1000 : : const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)=&
1001 : 787698 : & const_rot2nd(alpha,beta,gama,iatom,icoeff+ncoeff_prev)+terme1-terme2
1002 : : end do
1003 : : end do
1004 : : end do
1005 : : end do
1006 : : ! 2/ Enforce the symetry of the dynamical matrix. Number of constraints = (3*natom_unitcell)**2
1007 : : ! Note that we are unable to enforce the symetry when iatom=jatom (We have to write the equations)
1008 [ + + ]: 20648 : do alpha=1,3
1009 [ + + ]: 67106 : do beta=1,3
1010 [ + + ]: 278052 : do icoeff=1,ncoeff
1011 : : terme=sum(SS_ref(alpha,:,beta,isym,itrans)*Shell2at%proj(:,icoeff,ishell))-&
1012 [ + + + + ]: 4106052 : & sum(SS_ref(beta,:,alpha,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1013 : : const_dynmat(alpha,beta,iatom,iat_mod,icoeff+ncoeff_prev)=&
1014 : 262566 : & const_dynmat(alpha,beta,iatom,iat_mod,icoeff+ncoeff_prev)+terme
1015 : : end do
1016 : : end do
1017 : : end do
1018 : : ! 3/ Huang invariances. Number of constraints = 3**4
1019 [ + + ]: 22934 : do alpha=1,3
1020 [ + + ]: 67212 : do beta=1,3
1021 [ + + ]: 201318 : do gama=1,3
1022 [ + + ]: 603954 : do lambda=1,3
1023 [ + + ]: 2502468 : do icoeff=1,ncoeff
1024 : : terme=sum(SS_ref(alpha,:,beta,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*&
1025 : : & distance(iatom,jatom,gama+1)*&
1026 : : & distance(iatom,jatom,lambda+1)-&
1027 : : & sum(SS_ref(gama,:,lambda,isym,itrans)*Shell2at%proj(:,icoeff,ishell))*&
1028 : : & distance(iatom,jatom,alpha+1)*&
1029 [ + + + + ]: 36954468 : & distance(iatom,jatom,beta+1)
1030 : : const_huang(alpha,beta,gama,lambda,icoeff+ncoeff_prev)=&
1031 : 2363094 : & const_huang(alpha,beta,gama,lambda,icoeff+ncoeff_prev)+terme
1032 : : end do
1033 : : end do
1034 : : end do
1035 : : end do
1036 : : end do
1037 : : end do !iatshell
1038 : : end do !iatom
1039 : : end do !ishell
1040 : : end if !order=1,2
1041 : :
1042 : : ! Third order
1043 [ + + ]: 44 : if (order3) then
1044 : 14 : write(Invar%stdout,*) '########################## At the 3rd order #################################'
1045 [ - + + - : 659640 : ABI_MALLOC(Const3%AsrRot3,(natom_unitcell,natom,ntotcoeff))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + ]
1046 [ + + ]: 42 : do iatom=1,natom_unitcell
1047 [ + + ]: 4930 : do jatom=1,natom
1048 [ + + ]: 457300 : do icoeff=1,ntotcoeff
1049 [ - + - + : 18095360 : ABI_MALLOC(Const3%AsrRot3(iatom,jatom,icoeff)%ABG, (3,3,3)); Const3%AsrRot3(iatom,jatom,icoeff)%ABG(:,:,:) =zero
+ + + + +
+ ]
1050 [ - + - + : 54743352 : ABI_MALLOC(Const3%AsrRot3(iatom,jatom,icoeff)%ABGD,(3,3,3,3)); Const3%AsrRot3(iatom,jatom,icoeff)%ABGD(:,:,:,:)=zero
+ + + + +
+ + + ]
1051 : : end do
1052 : : end do
1053 : : end do
1054 [ + + ]: 165 : do ishell=1,Shell2at%nshell
1055 [ + + ]: 476 : do iatom=1,natom_unitcell
1056 [ + + ]: 311 : if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1057 [ + + ]: 2266 : do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
1058 : 1856 : jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1059 [ + + ]: 1856 : if (iatom==jatom) cycle
1060 : 1828 : isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1061 : 1828 : itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1062 : 1828 : ncoeff =Shell2at%ncoeff(ishell)
1063 : 1828 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
1064 : 1828 : iat_mod=mod(jatom+natom_unitcell-1,natom_unitcell)+1
1065 : : ! 1/ Rotational invariances (coming from the 2nd order). Number of constraints = natom_unitcell*natom*3**3
1066 [ + - ]: 2139 : if (Invar%order.ge.3) then
1067 [ + + ]: 7312 : do alpha=1,3
1068 [ + + ]: 23764 : do beta=1,3
1069 [ + + ]: 71292 : do gama=1,3
1070 [ + + ]: 213876 : do lambda=1,3
1071 [ + + ]: 824364 : do icoeff=1,ncoeff
1072 : 626940 : terme1=zero ; terme2=zero ; terme3=zero ; terme4=zero ;
1073 [ + + + + ]: 2507760 : if (alpha.eq.lambda) terme1=sum(SS_ref(gama ,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1074 [ + + + + ]: 2507760 : if (beta.eq.lambda) terme2=sum(SS_ref(alpha ,:,gama ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1075 [ + + + + ]: 2507760 : if (alpha.eq.gama) terme3=sum(SS_ref(lambda,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1076 [ + + + + ]: 2507760 : if (beta.eq.gama) terme4=sum(SS_ref(alpha ,:,lambda,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1077 [ + + ]: 626940 : if (distance(iatom,jatom,1).lt.Invar%rcut3) then
1078 : : Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1079 : 34668 : & Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1+terme2-terme3-terme4
1080 : : end if
1081 : : Const3%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1082 : 775008 : & Const3%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)-terme1-terme2+terme3+terme4
1083 : : end do
1084 : : end do
1085 : : end do
1086 : : end do
1087 : : end do
1088 : : end if !proj3rd
1089 : : end do !iatshell
1090 : : end do !iatom
1091 : : end do !ishell
1092 [ + + ]: 63 : do ishell=1,Shell3at%nshell
1093 [ + + ]: 176 : do iatom=1,natom_unitcell
1094 [ + + ]: 113 : if (Shell3at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1095 [ + + ]: 1286 : do iatshell=1,Shell3at%neighbours(iatom,ishell)%n_interactions
1096 : 1144 : jatom=Shell3at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1097 : 1144 : katom=Shell3at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1098 : 1144 : isym =Shell3at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1099 : 1144 : itrans=Shell3at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1100 : 1144 : ncoeff =Shell3at%ncoeff(ishell)
1101 : 1144 : ncoeff_prev=Shell3at%ncoeff_prev(ishell)+Solver%ncoeff2nd+Solver%ncoeff1st
1102 : : ! 2/ Acoustic sum rules (3rd order). Number of constraints = natom_unitcell*natom*3**3
1103 [ + + ]: 4576 : do alpha=1,3
1104 [ + + ]: 14872 : do beta=1,3
1105 [ + + ]: 44616 : do gama=1,3
1106 [ + + ]: 236826 : do icoeff=1,ncoeff
1107 [ + + ]: 5477976 : terme =sum(Const3%Sprod(isym,itrans)%SSS(alpha,:,beta,gama)*Shell3at%proj(:,icoeff,ishell))
1108 : : Const3%AsrRot3(iatom,katom,icoeff+ncoeff_prev)%ABG(alpha,beta,gama)=&
1109 : 226530 : & Const3%AsrRot3(iatom,katom,icoeff+ncoeff_prev)%ABG(alpha,beta,gama)+terme
1110 : : end do
1111 : : end do
1112 : : end do
1113 : : end do
1114 : : ! 2/ Rotational invariances (coming from the 3rd order). Number of constraints = natom_unitcell*natom*3**4
1115 [ + + ]: 4689 : do alpha=1,3
1116 [ + + ]: 14872 : do beta=1,3
1117 [ + + ]: 44616 : do gama=1,3
1118 [ + + ]: 133848 : do lambda=1,3
1119 [ + + ]: 710478 : do icoeff=1,ncoeff
1120 : : terme1=sum(Const3%Sprod(isym,itrans)%SSS(alpha,:,beta,gama )&
1121 [ + + ]: 16433928 : & *Shell3at%proj(:,icoeff,ishell))*distance(iatom,katom,lambda+1)
1122 : : terme2=sum(Const3%Sprod(isym,itrans)%SSS(alpha,:,beta,lambda)&
1123 [ + + ]: 16433928 : & *Shell3at%proj(:,icoeff,ishell))*distance(iatom,katom,gama+1)
1124 : : Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1125 : 679590 : & Const3%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1-terme2
1126 : : end do
1127 : : end do
1128 : : end do
1129 : : end do
1130 : : end do
1131 : : end do !iatshell
1132 : : end do !iatom
1133 : : end do !ishell
1134 : : end if !order=3
1135 : :
1136 : : ! Fourth order
1137 : : if (order4) then
1138 : : write(Invar%stdout,*) '########################## At the 4th order #################################'
1139 : : ABI_MALLOC(Const4%AsrRot4,(natom_unitcell,natom,natom,ntotcoeff))
1140 : : do iatom=1,natom_unitcell
1141 : : do jatom=1,natom
1142 : : do katom=1,natom
1143 : : do icoeff=1,ntotcoeff
1144 : : ABI_MALLOC(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGD, (3,3,3,3))
1145 : : Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGD(:,:,:,:) =zero
1146 : : !FB ABI_MALLOC(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGDE,(3,3,3,3,3))
1147 : : !FB Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGDE(:,:,:,:,:)=zero
1148 : : end do
1149 : : end do
1150 : : end do
1151 : : end do
1152 : : !FB do ishell=1,Shell2at%nshell
1153 : : !FB do iatom=1,natom_unitcell
1154 : : !FB if (Shell2at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1155 : : !FB do iatshell=1,Shell2at%neighbours(iatom,ishell)%n_interactions
1156 : : !FB jatom=Shell2at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1157 : : !FB if (iatom==jatom) cycle
1158 : : !FB isym=Shell2at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1159 : : !FB itrans=Shell2at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1160 : : !FB ncoeff =Shell2at%ncoeff(ishell)
1161 : : !FB ncoeff_prev=Shell2at%ncoeff_prev(ishell)+Solver%ncoeff1st
1162 : : !FB iat_mod=mod(jatom+natom_unitcell-1,natom_unitcell)+1
1163 : : !FB! 1/ Rotational invariances (coming from the 2nd order). Number of constraints = natom_unitcell*natom*3**3
1164 : : !FB if (Invar%order.ge.3) then
1165 : : !FB do alpha=1,3
1166 : : !FB do beta=1,3
1167 : : !FB do gama=1,3
1168 : : !FB do lambda=1,3
1169 : : !FB do icoeff=1,ncoeff
1170 : : !FB terme1=zero ; terme2=zero ; terme3=zero ; terme4=zero ;
1171 : : !FB if (alpha.eq.lambda) terme1=sum(SS_ref(gama ,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1172 : : !FB if (beta.eq.lambda) terme2=sum(SS_ref(alpha ,:,gama ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1173 : : !FB if (alpha.eq.gama) terme3=sum(SS_ref(lambda,:,beta ,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1174 : : !FB if (beta.eq.gama) terme4=sum(SS_ref(alpha ,:,lambda,isym,itrans)*Shell2at%proj(:,icoeff,ishell))
1175 : : !FB if (distance(iatom,jatom,1).lt.Invar%rcut3) then
1176 : : !FB Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1177 : : !FB& Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1+terme2-terme3-terme4
1178 : : !FB end if
1179 : : !FB Const4%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1180 : : !FB& Const4%AsrRot3(iatom,iatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)-terme1-terme2+terme3+terme4
1181 : : !FB end do
1182 : : !FB end do
1183 : : !FB end do
1184 : : !FB end do
1185 : : !FB end do
1186 : : !FB end if !proj3rd
1187 : : !FB end do !iatshell
1188 : : !FB end do !iatom
1189 : : !FB end do !ishell
1190 : : do ishell=1,Shell4at%nshell
1191 : : do iatom=1,natom_unitcell
1192 : : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
1193 : : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
1194 : : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
1195 : : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
1196 : : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
1197 : : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
1198 : : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
1199 : : ncoeff =Shell4at%ncoeff(ishell)
1200 : : ncoeff_prev=Shell4at%ncoeff_prev(ishell)+Solver%ncoeff3rd+Solver%ncoeff2nd+Solver%ncoeff1st
1201 : : ! 2/ Acoustic sum rules (4th order). Number of constraints = natom_unitcell*natom**2*3**4
1202 : : do alpha=1,3
1203 : : do beta=1,3
1204 : : do gama=1,3
1205 : : do delta=1,3
1206 : : do icoeff=1,ncoeff
1207 : : terme =sum(Const4%Sprod(isym,itrans)%SSSS(alpha,:,beta,gama,delta)*Shell4at%proj(:,icoeff,ishell))
1208 : : Const4%AsrRot4(iatom,katom,latom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,delta)=&
1209 : : & Const4%AsrRot4(iatom,katom,latom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,delta)+terme
1210 : : end do
1211 : : end do
1212 : : end do
1213 : : end do
1214 : : end do
1215 : : !FB! 2/ Rotational invariances (coming from the 3rd order). Number of constraints = natom_unitcell*natom*3**4
1216 : : !FB do alpha=1,3
1217 : : !FB do beta=1,3
1218 : : !FB do gama=1,3
1219 : : !FB do lambda=1,3
1220 : : !FB do icoeff=1,ncoeff
1221 : : !FB terme1=sum(Const4%Sprod(isym,itrans)%SSS(alpha,:,beta,gama )&
1222 : : !FB &*Shell3at%proj(:,icoeff,ishell))*distance(iatom,katom,lambda+1)
1223 : : !FB terme2=sum(Const4%Sprod(isym,itrans)%SSS(alpha,:,beta,lambda)&
1224 : : !FB &*Shell3at%proj(:,icoeff,ishell))*distance(iatom,katom,gama+1)
1225 : : !FB Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)=&
1226 : : !FB& Const4%AsrRot3(iatom,jatom,icoeff+ncoeff_prev)%ABGD(alpha,beta,gama,lambda)+terme1-terme2
1227 : : !FB end do
1228 : : !FB end do
1229 : : !FB end do
1230 : : !FB end do
1231 : : !FB end do
1232 : : end do !iatshell
1233 : : end do !iatom
1234 : : end do !ishell
1235 : : end if !order=4
1236 : :
1237 [ + - ]: 44 : if (order2.or.order3) then
1238 [ - + ]: 44 : ABI_FREE(SS_ref)
1239 : : end if
1240 [ + + ]: 44 : if (order3.or.order4) then
1241 [ + + ]: 510 : do isym=1,Sym%nsym
1242 [ + + ]: 3486 : do itrans=1,6
1243 [ - + ]: 3472 : ABI_FREE(Const3%Sprod(isym,itrans)%SSS)
1244 : : end do
1245 : : end do
1246 [ + - + + : 2990 : ABI_FREE(Const3%Sprod)
- + - + -
+ - + ]
1247 : : end if
1248 : : if (order4) then
1249 : : do isym=1,Sym%nsym
1250 : : do itrans=1,24
1251 : : ABI_FREE(Const4%Sprod(isym,itrans)%SSSS)
1252 : : end do
1253 : : end do
1254 : : ABI_FREE(Const4%Sprod)
1255 : : end if
1256 : :
1257 : : ! Reduce the number of constraints by selecting the non-zero equations
1258 : 44 : write(Invar%stdout,*) '################## Reduce the number of constraints #########################'
1259 : 44 : iconst_new=0
1260 [ + - ]: 44 : if (order2) then
1261 : : ! 1/ For Rotational invariances (1st order)
1262 : 44 : iconst=0
1263 [ - + + - : 25241 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_1st)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1264 [ + + ]: 176 : do alpha=1,3
1265 [ + + ]: 572 : do beta=1,3
1266 : 396 : iconst=iconst+1
1267 [ + + ]: 25197 : vect(:,iconst)=const_rot1st(alpha,beta,:)
1268 : : end do
1269 : : end do
1270 : 44 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_1st,nconst_loc,vect)
1271 [ - + ]: 44 : if (nconst_loc.ne.0) then
1272 [ # # ]: 0 : do iconst_loc=1,nconst_loc
1273 : 0 : iconst_new=iconst_new+1
1274 [ # # ]: 0 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1275 : : end do
1276 : : end if
1277 [ - + ]: 44 : ABI_FREE(vect)
1278 [ - + ]: 44 : ABI_FREE(const_rot1st)
1279 : 44 : Solver%nconst_1st=nconst_loc
1280 : :
1281 : : ! 2/ For Rotational invariances (2nd order)
1282 : 44 : iconst=0
1283 [ - + + - : 366593 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_rot2nd)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1284 [ + + ]: 150 : do iatom=1,natom_unitcell
1285 [ + + ]: 468 : do alpha=1,3
1286 [ + + ]: 1378 : do beta=1,3
1287 [ + + ]: 4134 : do gama=1,3
1288 : 2862 : iconst=iconst+1
1289 [ + + ]: 367371 : vect(:,iconst)=const_rot2nd(alpha,beta,gama,iatom,:)
1290 : : end do
1291 : : end do
1292 : : end do
1293 : : end do
1294 : 44 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_rot2nd,nconst_loc,vect)
1295 [ + + ]: 44 : if (nconst_loc.ne.0) then
1296 [ + + ]: 113 : do iconst_loc=1,nconst_loc
1297 : 100 : iconst_new=iconst_new+1
1298 [ + + ]: 33151 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1299 : : end do
1300 : : end if
1301 [ - + ]: 44 : ABI_FREE(vect)
1302 [ - + ]: 44 : ABI_FREE(const_rot2nd)
1303 : 44 : Solver%nconst_rot2nd=nconst_loc
1304 : :
1305 : : ! 3/ For symetry of the dynamical matrix
1306 : 44 : iconst=0
1307 [ - + + - : 1278239 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_dynmat)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1308 [ + + ]: 150 : do iatom=1,natom_unitcell
1309 [ + + ]: 664 : do jatom=1,natom_unitcell
1310 [ + + ]: 2162 : do alpha=1,3
1311 [ + + ]: 6682 : do beta=1,3
1312 : 4626 : iconst=iconst+1
1313 [ + + ]: 1279605 : vect(:,iconst)=const_dynmat(alpha,beta,iatom,jatom,:)
1314 : : end do
1315 : : end do
1316 : : end do
1317 : : end do
1318 : 44 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_dynmat,nconst_loc,vect)
1319 [ + + ]: 44 : if (nconst_loc.ne.0) then
1320 [ + + ]: 69 : do iconst_loc=1,nconst_loc
1321 : 65 : iconst_new=iconst_new+1
1322 [ + + ]: 26885 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1323 : : end do
1324 : : end if
1325 [ - + ]: 44 : ABI_FREE(vect)
1326 [ - + ]: 44 : ABI_FREE(const_dynmat)
1327 : 44 : Solver%nconst_dynmat=nconst_loc
1328 : :
1329 : : ! 4/ For Huang invariances
1330 : 44 : iconst=0
1331 [ - + + - : 225761 : ABI_MALLOC(vect,(ntotcoeff,Solver%nconst_huang)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1332 [ + + ]: 176 : do alpha=1,3
1333 [ + + ]: 572 : do beta=1,3
1334 [ + + ]: 1716 : do gama=1,3
1335 [ + + ]: 5148 : do lambda=1,3
1336 : 3564 : iconst=iconst+1
1337 [ + + ]: 226773 : vect(:,iconst)=const_huang(alpha,beta,gama,lambda,:)
1338 : : end do
1339 : : end do
1340 : : end do
1341 : : end do
1342 : 44 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_huang,nconst_loc,vect)
1343 [ + + ]: 44 : if (nconst_loc.ne.0) then
1344 [ + + ]: 70 : do iconst_loc=1,nconst_loc
1345 : 51 : iconst_new=iconst_new+1
1346 [ + + ]: 6515 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1347 : : end do
1348 : : end if
1349 [ - + ]: 44 : ABI_FREE(vect)
1350 [ - + ]: 44 : ABI_FREE(const_huang)
1351 : 44 : Solver%nconst_huang=nconst_loc
1352 : 220 : Solver%nconst_2nd=Solver%nconst_rot2nd+Solver%nconst_dynmat+Solver%nconst_huang
1353 : : end if
1354 : :
1355 [ + + ]: 44 : if (order3) then
1356 : : ! 1/ For acoustic sum rules (3rd order)
1357 : 14 : iconst=0
1358 [ - + + - : 12346400 : ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_asr3rd)) ; vect (:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1359 [ + + ]: 42 : do iatom=1,natom_unitcell
1360 [ + + ]: 4930 : do jatom=1,natom
1361 [ + + ]: 19580 : do alpha=1,3
1362 [ + + ]: 63544 : do beta=1,3
1363 [ + + ]: 190632 : do gama=1,3
1364 : 131976 : iconst=iconst+1
1365 [ + + ]: 12390336 : do ii=1,ntotcoeff
1366 : 12346344 : vect(ii,iconst)=Const3%AsrRot3(iatom,jatom,ii)%ABG(alpha,beta,gama)
1367 : : end do
1368 : : end do
1369 : : end do
1370 : : end do
1371 : : end do
1372 : : end do
1373 : 14 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_asr3rd,nconst_loc,vect)
1374 [ + + ]: 14 : if (nconst_loc.ne.0) then
1375 [ + + ]: 75 : do iconst_loc=1,nconst_loc
1376 : 64 : iconst_new=iconst_new+1
1377 [ + + ]: 8571 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1378 : : end do
1379 : : end if
1380 [ - + ]: 14 : ABI_FREE(vect)
1381 : 14 : Solver%nconst_asr3rd=nconst_loc
1382 : :
1383 : : ! 2/ For Rotational invariances (3rd order)
1384 : 14 : iconst=0
1385 [ - + + - : 37039088 : ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_rot3rd)) ; vect(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
1386 [ + + ]: 42 : do iatom=1,natom_unitcell
1387 [ + + ]: 4930 : do jatom=1,natom
1388 [ + + ]: 19580 : do alpha=1,3
1389 [ + + ]: 63544 : do beta=1,3
1390 [ + + ]: 190632 : do gama=1,3
1391 [ + + ]: 571896 : do lambda=1,3
1392 : 395928 : iconst=iconst+1
1393 [ + + ]: 37171008 : do ii=1,ntotcoeff
1394 : 37039032 : vect(ii,iconst)=Const3%AsrRot3(iatom,jatom,ii)%ABGD(alpha,beta,gama,lambda)
1395 : : end do
1396 : : end do
1397 : : end do
1398 : : end do
1399 : : end do
1400 : : end do
1401 : : end do
1402 : 14 : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_rot3rd,nconst_loc,vect)
1403 [ + + ]: 14 : if (nconst_loc.ne.0) then
1404 [ + + ]: 101 : do iconst_loc=1,nconst_loc
1405 : 90 : iconst_new=iconst_new+1
1406 [ + + ]: 13631 : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1407 : : end do
1408 : : end if
1409 [ - + ]: 14 : ABI_FREE(vect)
1410 [ + + ]: 42 : do iatom=1,natom_unitcell
1411 [ + + ]: 4930 : do jatom=1,natom
1412 [ + + ]: 457300 : do icoeff=1,ntotcoeff
1413 [ - + ]: 452384 : ABI_FREE(Const3%AsrRot3(iatom,jatom,icoeff)%ABG)
1414 [ - + ]: 457272 : ABI_FREE(Const3%AsrRot3(iatom,jatom,icoeff)%ABGD)
1415 : : end do
1416 : : end do
1417 : : end do
1418 [ + - + + : 452398 : ABI_FREE(Const3%AsrRot3)
- + - + -
+ - + ]
1419 : 14 : Solver%nconst_rot3rd=nconst_loc
1420 : 42 : Solver%nconst_3rd=Solver%nconst_asr3rd+Solver%nconst_rot3rd
1421 : : end if
1422 : :
1423 : : if (order4) then
1424 : : ! 1/ For acoustic sum rules (4th order)
1425 : : iconst=0
1426 : : ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_asr4th)) ; vect (:,:)=zero
1427 : : do iatom=1,natom_unitcell
1428 : : do jatom=1,natom
1429 : : do katom=1,natom
1430 : : do alpha=1,3
1431 : : do beta=1,3
1432 : : do gama=1,3
1433 : : do delta=1,3
1434 : : iconst=iconst+1
1435 : : do ii=1,ntotcoeff
1436 : : vect(ii,iconst)=Const4%AsrRot4(iatom,jatom,katom,ii)%ABGD(alpha,beta,gama,delta)
1437 : : end do
1438 : : end do
1439 : : end do
1440 : : end do
1441 : : end do
1442 : : end do
1443 : : end do
1444 : : end do
1445 : : call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_asr4th,nconst_loc,vect)
1446 : : if (nconst_loc.ne.0) then
1447 : : do iconst_loc=1,nconst_loc
1448 : : iconst_new=iconst_new+1
1449 : : Solver%const(iconst_new,:)=vect(:,iconst_loc)
1450 : : end do
1451 : : end if
1452 : : ABI_FREE(vect)
1453 : : Solver%nconst_asr4th=nconst_loc
1454 : :
1455 : : !FB! 2/ For Rotational invariances (3rd order)
1456 : : !FB iconst=0
1457 : : !FB ABI_MALLOC(vect ,(ntotcoeff,Solver%nconst_rot3rd)) ; vect(:,:)=zero
1458 : : !FB do iatom=1,natom_unitcell
1459 : : !FB do jatom=1,natom
1460 : : !FB do alpha=1,3
1461 : : !FB do beta=1,3
1462 : : !FB do gama=1,3
1463 : : !FB do lambda=1,3
1464 : : !FB iconst=iconst+1
1465 : : !FB vect(:,iconst)=Const4%AsrRot3(iatom,jatom,:)%ABGD(alpha,beta,gama,lambda)
1466 : : !FB end do
1467 : : !FB end do
1468 : : !FB end do
1469 : : !FB end do
1470 : : !FB end do
1471 : : !FB end do
1472 : : !FB call tdep_calc_orthonorm(ntotcoeff,Solver%nconst_rot3rd,nconst_loc,vect)
1473 : : !FB if (nconst_loc.ne.0) then
1474 : : !FB do iconst_loc=1,nconst_loc
1475 : : !FB iconst_new=iconst_new+1
1476 : : !FB Solver%const(iconst_new,:)=vect(:,iconst_loc)
1477 : : !FB end do
1478 : : !FB end if
1479 : : !FB ABI_FREE(vect)
1480 : : !FB Solver%nconst_rot3rd=nconst_loc
1481 : : do iatom=1,natom_unitcell
1482 : : do jatom=1,natom
1483 : : do katom=1,natom
1484 : : do icoeff=1,ntotcoeff
1485 : : ABI_FREE(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGD)
1486 : : !FB ABI_FREE(Const4%AsrRot4(iatom,jatom,katom,icoeff)%ABGDE)
1487 : : end do
1488 : : end do
1489 : : end do
1490 : : end do
1491 : : ABI_FREE(Const4%AsrRot4)
1492 : : Solver%nconst_rot4th=0
1493 : : Solver%nconst_4th=Solver%nconst_asr4th+Solver%nconst_rot4th
1494 : : end if
1495 : :
1496 : : ! Finalize the orthonormalization
1497 : : !FB nconst=iconst_new
1498 : : !FB ABI_MALLOC(vect,(ntotcoeff,nconst)) ; vect(:,:)=zero
1499 : : !FB do iconst=1,nconst
1500 : : !FB vect(:,iconst)=Solver%const(iconst,:)
1501 : : !FB end do
1502 : : !FB ABI_FREE(Solver%const)
1503 : : !FB call tdep_calc_orthonorm(ntotcoeff,nconst,nconst_loc,vect)
1504 : : !FB Solver%ntotconst=nconst_loc
1505 : : !FB if (nconst_loc.ne.0) then
1506 : : !FB ABI_MALLOC(Solver%const ,(Solver%ntotconst,ntotcoeff)); Solver%const (:,:)=0.d0
1507 : : !FB do iconst_loc=1,nconst_loc
1508 : : !FB Solver%const(iconst_loc,:)=vect(:,iconst_loc)
1509 : : !FB end do
1510 : : !FB end if
1511 : : !FB ABI_FREE(vect)
1512 : 44 : Solver%ntotconst=iconst_new
1513 [ + - ]: 44 : if (MPIdata%iam_master) then
1514 : 44 : open(unit=16,file=trim(Invar%output_prefix)//'_constraints.dat')
1515 : 44 : write(16,*) ' ======== Constraints at the 1st order (Rotational Invariances) ========'
1516 : 44 : write(16,*) ' Number of constraints =',Solver%nconst_1st
1517 : 44 : write(16,*) ' ======== Constraints at the 2nd order (Rotational Invariances) ========'
1518 : 44 : write(16,*) ' Number of constraints =',Solver%nconst_rot2nd
1519 : 44 : write(16,*) ' ======== Constraints at the 2nd order (Dynamical Matrix) =============='
1520 : 44 : write(16,*) ' Number of constraints =',Solver%nconst_dynmat
1521 : 44 : write(16,*) ' ======== Constraints at the 2nd order (Huang) ========================='
1522 : 44 : write(16,*) ' Number of constraints =',Solver%nconst_huang
1523 [ + + ]: 44 : if (Invar%order.ge.3) then
1524 : 14 : write(16,*) ' ======== Constraints at the 3rd order (Acoustic sum rules) ============'
1525 : 14 : write(16,*) ' Number of constraints =',Solver%nconst_asr3rd
1526 : 14 : write(16,*) ' ======== Constraints at the 3rd order (Rotational Invariances) ========'
1527 : 14 : write(16,*) ' Number of constraints =',Solver%nconst_rot3rd
1528 : : end if
1529 [ + + ]: 44 : if (Invar%order.ge.4) then
1530 : 8 : write(16,*) ' ======== Constraints at the 4th order (Acoustic sum rules) ============'
1531 : 8 : write(16,*) ' Number of constraints =',Solver%nconst_asr4th
1532 : 8 : write(16,*) ' ======== Constraints at the 4th order (Rotational Invariances) ========'
1533 : 8 : write(16,*) ' Number of constraints =',Solver%nconst_rot4th
1534 : : end if
1535 : 44 : write(16,*) ' ======================================================================='
1536 : 44 : write(16,*) ' Total number of constraints =',Solver%ntotconst
1537 : 44 : close(16)
1538 : : end if
1539 : :
1540 [ - + - + : 44 : end subroutine tdep_calc_constraints
- + - + -
+ - + - -
- - - - -
- - + - -
- - - - -
- ]
1541 : :
1542 : : !====================================================================================================
1543 : :
1544 [ # # # # : 0 : end module m_tdep_solver
# # # # #
# # # # #
# # # # #
# ]
1545 : : !!***
|
