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 : !!***
|
