Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : :
6 : : #include "abi_common.h"
7 : :
8 : : module m_tdep_phi4
9 : :
10 : : use defs_basis
11 : : use m_errors
12 : : use m_abicore
13 : : use m_numeric_tools
14 : : use m_linalg_interfaces
15 : : use m_io_tools
16 : : use m_crystal, only : crystal_t
17 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
18 : : use m_tdep_latt, only : Lattice_type
19 : : use m_tdep_shell, only : Shell_type
20 : : use m_tdep_sym, only : Symmetries_type
21 : : use m_tdep_sampling, only : tdep_Sampling_type
22 : : use m_tdep_solver, only : tdep_Solver_type
23 : : use m_tdep_model, only : tdep_Model_type
24 : : use m_tdep_constraints, only : Constraints_type
25 : :
26 : : implicit none
27 : :
28 : : public :: tdep_calc_ftot4
29 : : public :: tdep_calc_phi4ref
30 : : public :: tdep_write_phi4
31 : : public :: tdep_build_phi4_3333
32 : :
33 : : contains
34 : :
35 : : !====================================================================================================
36 : :
37 : 8 : subroutine tdep_calc_ftot4(Model,Invar,Shell4at,ucart,Sym)
38 : :
39 : : type(tdep_Model_type),intent(inout) :: Model
40 : : type(atdep_dataset_type),intent(in) :: Invar
41 : : type(Shell_type),intent(in) :: Shell4at
42 : : type(Symmetries_type),intent(in) :: Sym
43 : : double precision, intent(in) :: ucart(3,Invar%natom,Invar%my_nstep)
44 : :
45 : : integer :: iatom,jatom,katom,latom,isym,itrans,ishell,iatshell
46 : : integer :: ii,jj,kk,ll,istep
47 : 8 : double precision, allocatable :: Phi4_3333(:,:,:,:)
48 : 8 : double precision, allocatable :: ucart_blas(:)
49 : 8 : double precision, allocatable :: ftot4(:,:)
50 : :
51 [ + - + + : 976 : ABI_MALLOC(Phi4_3333,(3,3,3,3)) ; Phi4_3333(:,:,:,:)=0.d0
+ + + + +
+ ]
52 [ + - - + : 88760 : ABI_MALLOC(ftot4,(3*Invar%natom,Invar%my_nstep)); ftot4(:,:)=0.d0
+ - - + +
- + - + -
- + - + +
+ + + ]
53 [ + + ]: 1484 : do iatom=1,Invar%natom
54 [ + + ]: 9316 : do ishell=1,Shell4at%nshell
55 : : ! Build the 3x3x3x3 IFC of an atom in this shell
56 [ + + ]: 7832 : if (Shell4at%neighbours(iatom,ishell)%n_interactions.eq.0) cycle
57 [ + + ]: 252016 : do iatshell=1,Shell4at%neighbours(iatom,ishell)%n_interactions
58 : 243732 : jatom=Shell4at%neighbours(iatom,ishell)%atomj_in_shell(iatshell)
59 : 243732 : katom=Shell4at%neighbours(iatom,ishell)%atomk_in_shell(iatshell)
60 : 243732 : latom=Shell4at%neighbours(iatom,ishell)%atoml_in_shell(iatshell)
61 : 243732 : isym =Shell4at%neighbours(iatom,ishell)%sym_in_shell(iatshell)
62 : 243732 : itrans=Shell4at%neighbours(iatom,ishell)%transpose_in_shell(iatshell)
63 : 243732 : call tdep_build_phi4_3333(isym,Model%Phi4(:,:,:,:,ishell),Phi4_3333,Sym,itrans)
64 : : ! Calculation of the force components (third order)
65 [ + + ]: 5126204 : do istep=1,Invar%my_nstep
66 [ + + ]: 19742292 : do ii=1,3
67 [ + + ]: 63370320 : do jj=1,3
68 [ + + ]: 190110960 : do kk=1,3
69 [ + + ]: 570332880 : do ll=1,3
70 : : ftot4(3*(iatom-1)+ii,istep)=ftot4(3*(iatom-1)+ii,istep)+&
71 : 526461120 : & Phi4_3333(ii,jj,kk,ll)*ucart(jj,jatom,istep)*ucart(kk,katom,istep)*ucart(ll,latom,istep)
72 : : end do !ll
73 : : end do !kk
74 : : end do !jj
75 : : end do !ii
76 : : end do !istep
77 : : end do !iatshell
78 : : end do !ishell
79 : : end do !iatom
80 [ - + ]: 8 : ABI_FREE(Phi4_3333)
81 [ + + + + ]: 88728 : ftot4(:,:)=ftot4(:,:)/6.d0
82 : :
83 [ - + + - : 4460 : ABI_MALLOC(ucart_blas ,(3*Invar%natom)) ; ucart_blas (:)=0.d0
+ - + - -
+ + - +
+ ]
84 [ + + ]: 168 : do istep=1,Invar%my_nstep
85 [ + + ]: 88720 : ucart_blas(:)=0.d0
86 [ + + ]: 29680 : do jatom=1,Invar%natom
87 [ + + ]: 118240 : do jj=1,3
88 : 118080 : ucart_blas(3*(jatom-1)+jj)=ucart(jj,jatom,istep)
89 : : end do
90 : : end do
91 : : call DGEMM('T','N',1,1,3*Invar%natom,1./4.d0,ftot4(:,istep),3*Invar%natom,ucart_blas,&
92 : 160 : & 3*Invar%natom,0.d0,Model%Phi4UiUjUkUl(istep),3*Invar%natom)
93 : : Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)=&
94 [ + + ]: 88728 : & Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)-ftot4(:,istep)
95 : : end do
96 [ - + ]: 8 : ABI_FREE(ucart_blas)
97 [ - + ]: 8 : ABI_FREE(ftot4)
98 : :
99 : 8 : end subroutine tdep_calc_ftot4
100 : :
101 : : !=====================================================================================================
102 : 8 : subroutine tdep_calc_phi4ref(Solver,Shell4at,Phi4_ref)
103 : :
104 : : type(tdep_Solver_type),intent(in) :: Solver
105 : : type(Shell_type),intent(in) :: Shell4at
106 : : double precision, intent(inout) :: Phi4_ref(3,3,3,3,Shell4at%nshell)
107 : :
108 : : integer :: ishell,ncoeff,ncoeff_prev
109 : : integer :: ii,jj,kk,ll,kappa
110 : 8 : double precision, allocatable :: Phi4_coeff(:)
111 : :
112 [ - + + - : 566 : ABI_CALLOC(Phi4_coeff, (Solver%ncoeff4th))
+ - + - -
+ - + +
+ ]
113 [ + + ]: 550 : Phi4_coeff(:) = Solver%theta(Solver%ncoeff1st+Solver%ncoeff2nd+Solver%ncoeff3rd+1:Solver%ntotcoeff)
114 : :
115 [ + + ]: 54 : do ishell=1,Shell4at%nshell
116 : : ! Build the 3x3x3x3 IFC per shell
117 : 46 : ncoeff =Shell4at%ncoeff(ishell)
118 : 46 : ncoeff_prev=Shell4at%ncoeff_prev(ishell)
119 : 46 : kappa=0
120 [ + + ]: 184 : do ii=1,3
121 [ + + ]: 598 : do jj=1,3
122 [ + + ]: 1794 : do kk=1,3
123 [ + + ]: 5382 : do ll=1,3
124 : 3726 : kappa=kappa+1
125 [ + + ]: 48870 : Phi4_ref(ii,jj,kk,ll,ishell)=sum(Shell4at%proj(kappa,1:ncoeff,ishell)*Phi4_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
126 : : end do
127 : : end do
128 : : end do
129 : : end do
130 : : ! Remove the rounding errors before writing (for non regression testing purposes)
131 [ + + ]: 192 : do ii=1,3
132 [ + + ]: 598 : do jj=1,3
133 [ + + ]: 1794 : do kk=1,3
134 [ + + ]: 5382 : do ll=1,3
135 [ + + ]: 4968 : if (abs(Phi4_ref(ii,jj,kk,ll,ishell)).lt.tol8) Phi4_ref(ii,jj,kk,ll,ishell)=zero
136 : : end do
137 : : end do
138 : : end do
139 : : end do
140 : : end do
141 : :
142 [ - + ]: 8 : ABI_FREE(Phi4_coeff)
143 : :
144 : 8 : end subroutine tdep_calc_phi4ref
145 : :
146 : : !=====================================================================================================
147 : 8 : subroutine tdep_write_phi4(distance,Invar,Phi4_ref,Shell4at,Sym)
148 : :
149 : : type(atdep_dataset_type),intent(in) :: Invar
150 : : type(Symmetries_type),intent(in) :: Sym
151 : : type(Shell_type),intent(in) :: Shell4at
152 : : double precision, intent(in) :: distance(Invar%natom,Invar%natom,4)
153 : : double precision, intent(in) :: Phi4_ref(3,3,3,3,Shell4at%nshell)
154 : :
155 : : integer :: ishell,isym,jatom,katom,latom
156 : : integer :: iatref,jatref,katref,latref,iatshell,itrans
157 : : integer :: ii,jj,kk
158 : : double precision :: tmp1,tmp2,tmp3
159 : 8 : double precision, allocatable :: Phi4_3333(:,:,:,:)
160 : :
161 : 8 : write(Invar%stdout,*) ' '
162 : 8 : write(Invar%stdout,*) '#############################################################################'
163 : 8 : write(Invar%stdout,*) '#### For each shell, list of coefficients (IFC), number of neighbours... ####'
164 : 8 : write(Invar%stdout,*) '#############################################################################'
165 : :
166 : : ! Write the IFCs in the data.out file (with others specifications:
167 : : ! number of atoms in a shell, Trace...)
168 [ + - + + : 976 : ABI_MALLOC(Phi4_3333,(3,3,3,3)) ; Phi4_3333(:,:,:,:)=0.d0
+ + + + +
+ ]
169 [ + + ]: 54 : do ishell=1,Shell4at%nshell
170 : 46 : iatref=Shell4at%iatref(ishell)
171 [ + - ]: 54 : if (Shell4at%neighbours(iatref,ishell)%n_interactions.ne.0) then
172 : 46 : jatref=Shell4at%jatref(ishell)
173 : 46 : katref=Shell4at%katref(ishell)
174 : 46 : latref=Shell4at%latref(ishell)
175 : 46 : write(Invar%stdout,'(a,i4,a,i4,a)') ' ======== NEW SHELL (ishell=',ishell,&
176 : 92 : & '): There are',Shell4at%neighbours(iatref,ishell)%n_interactions,' atoms on this shell'
177 [ + + ]: 1677 : do iatshell=1,Shell4at%neighbours(iatref,ishell)%n_interactions
178 : 1631 : jatom =Shell4at%neighbours(iatref,ishell)%atomj_in_shell(iatshell)
179 : 1631 : katom =Shell4at%neighbours(iatref,ishell)%atomk_in_shell(iatshell)
180 : 1631 : latom =Shell4at%neighbours(iatref,ishell)%atoml_in_shell(iatshell)
181 : 1631 : isym =Shell4at%neighbours(iatref,ishell)%sym_in_shell(iatshell)
182 : 1631 : itrans=Shell4at%neighbours(iatref,ishell)%transpose_in_shell(iatshell)
183 : 1631 : call tdep_build_phi4_3333(isym,Phi4_ref(:,:,:,:,ishell),Phi4_3333,Sym,itrans)
184 : 1631 : write(Invar%stdout,'(a,i4,a,i4)') ' For iatcell=',iatref,' ,with type=',mod(iatref-1,Invar%natom_unitcell)+1
185 : 1631 : write(Invar%stdout,'(a,i4,a,i4)') ' For jatom =',jatom ,' ,with type=',mod(jatom -1,Invar%natom_unitcell)+1
186 : 1631 : write(Invar%stdout,'(a,i4,a,i4)') ' For katom =',katom ,' ,with type=',mod(katom -1,Invar%natom_unitcell)+1
187 : 1631 : write(Invar%stdout,'(a,i4,a,i4)') ' For latom =',latom ,' ,with type=',mod(latom -1,Invar%natom_unitcell)+1
188 [ + + ]: 6524 : do ii=1,3
189 [ + + ]: 21203 : do jj=1,3
190 : : #if defined FC_NVHPC
191 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
192 : : #endif
193 : :
194 : 14679 : write(Invar%stdout,'(a,i2,i2,a)') ' Phi4^{',ii,jj,'kl}='
195 [ + + ]: 63609 : do kk=1,3
196 [ + + ]: 44037 : if (abs(Phi4_3333(ii,jj,1,kk)).lt.5.d-7) then
197 : 7416 : tmp1=0.d0
198 : : else
199 : 36621 : tmp1=Phi4_3333(ii,jj,1,kk)
200 : : end if
201 [ + + ]: 44037 : if (abs(Phi4_3333(ii,jj,2,kk)).lt.5.d-7) then
202 : 7306 : tmp2=0.d0
203 : : else
204 : 36731 : tmp2=Phi4_3333(ii,jj,2,kk)
205 : : end if
206 [ + + ]: 44037 : if (abs(Phi4_3333(ii,jj,3,kk)).lt.5.d-7) then
207 : 7410 : tmp3=0.d0
208 : : else
209 : 36627 : tmp3=Phi4_3333(ii,jj,3,kk)
210 : : end if
211 : 58716 : write(Invar%stdout,'(2x,3(f9.6,1x))') tmp1,tmp2,tmp3
212 : : end do
213 : : end do
214 : : end do
215 : 1631 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (i,j) vector components:', (distance(iatref,jatom,jj+1),jj=1,3)
216 : 1631 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (j,k) vector components:', (distance(jatom ,katom,jj+1),jj=1,3)
217 : 1631 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (j,k) vector components:', (distance(katom ,latom,jj+1),jj=1,3)
218 : 1631 : write(Invar%stdout,'(a,3(f9.6,1x))') ' (k,i) vector components:', (distance(latom,iatref,jj+1),jj=1,3)
219 : 1677 : write(Invar%stdout,*) ' '
220 : : end do !iatshell
221 : : end if !n_interactions
222 : : end do !ishell
223 [ - + ]: 8 : ABI_FREE(Phi4_3333)
224 : :
225 : 8 : end subroutine tdep_write_phi4
226 : :
227 : : !=====================================================================================================
228 : 245363 : subroutine tdep_build_phi4_3333(isym,Phi4_ref,Phi4_3333,Sym,itrans)
229 : :
230 : : type(Symmetries_type),intent(in) :: Sym
231 : : double precision, intent(in) :: Phi4_ref(3,3,3,3)
232 : : double precision, intent(out) :: Phi4_3333(3,3,3,3)
233 : : integer,intent(in) :: isym,itrans
234 : :
235 : : integer :: alpha,beta,gama,delta
236 : : integer :: ii,jj,kk,ll,ee,ff,gg,hh,mu,nu,xi,zeta
237 : : double precision :: Phi4_tmp(3,3,3,3)
238 : :
239 : :
240 : : ! Transform in the new basis wrt S_ref
241 : 245363 : Phi4_3333(:,:,:,:)=zero
242 [ + + ]: 981452 : do mu=1,3
243 [ + + ]: 3189719 : do alpha=1,3
244 [ + + ]: 9569157 : do nu=1,3
245 [ + + ]: 28707471 : do beta=1,3
246 [ + + ]: 86122413 : do xi=1,3
247 [ + + ]: 258367239 : do gama=1,3
248 [ + + ]: 775101717 : do zeta=1,3
249 [ + + ]: 2325305151 : do delta=1,3
250 : : #if defined FC_NVHPC
251 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
252 : : #endif
253 : : Phi4_3333(mu,nu,xi,zeta)=Phi4_3333(mu,nu,xi,zeta)+&
254 : : & Sym%S_ref(mu,alpha,isym,1)*Sym%S_ref(nu ,beta ,isym,1)*&
255 : 2146435524 : & Sym%S_ref(xi,gama ,isym,1)*Sym%S_ref(zeta,delta,isym,1)*Phi4_ref(alpha,beta,gama,delta)
256 : : end do
257 : : end do
258 : : end do
259 : : end do
260 : : end do
261 : : end do
262 : : end do
263 : : end do
264 : :
265 : : ! Take into account the 6 allowed permutations
266 : 245363 : Phi4_tmp(:,:,:,:)=Phi4_3333(:,:,:,:)
267 [ - + ]: 245363 : if ((itrans.lt.1).or.(itrans.gt.24)) then
268 : 0 : ABI_BUG('This value of the symmetry index is not permitted')
269 : : end if
270 [ + + ]: 981452 : do ii=1,3
271 [ + + ]: 3189719 : do jj=1,3
272 [ + + ]: 9569157 : do kk=1,3
273 [ + + ]: 28707471 : do ll=1,3
274 : : #if defined FC_NVHPC
275 : : if (itrans == -1) write(std_out, *)"NVHPC freezes here that is fixed by this print statement."
276 : : #endif
277 : :
278 [ + + ]: 19874403 : if (itrans==1) then ; ee=ii ; ff=jj ; gg=kk ; hh=ll ; endif !\Phi4_ijkl
279 [ + + ]: 19874403 : if (itrans==2) then ; ee=ii ; ff=kk ; gg=jj ; hh=ll ; endif !\Phi4_ikjl
280 [ + + ]: 19874403 : if (itrans==3) then ; ee=jj ; ff=ii ; gg=kk ; hh=ll ; endif !\Phi4_jikl
281 [ + + ]: 19874403 : if (itrans==4) then ; ee=jj ; ff=kk ; gg=ii ; hh=ll ; endif !\Phi4_jkil
282 [ + + ]: 19874403 : if (itrans==5) then ; ee=kk ; ff=ii ; gg=jj ; hh=ll ; endif !\Phi4_kijl
283 [ + + ]: 19874403 : if (itrans==6) then ; ee=kk ; ff=jj ; gg=ii ; hh=ll ; endif !\Phi4_kjil
284 : :
285 [ + + ]: 19874403 : if (itrans==7 ) then ; ee=ii ; ff=jj ; gg=ll ; hh=kk ; endif !\Phi4_ijlk
286 [ + + ]: 19874403 : if (itrans==8 ) then ; ee=ii ; ff=kk ; gg=ll ; hh=jj ; endif !\Phi4_iklj
287 [ + + ]: 19874403 : if (itrans==9 ) then ; ee=jj ; ff=ii ; gg=ll ; hh=kk ; endif !\Phi4_jilk
288 [ + + ]: 19874403 : if (itrans==10) then ; ee=jj ; ff=kk ; gg=ll ; hh=ii ; endif !\Phi4_jkli
289 [ + + ]: 19874403 : if (itrans==11) then ; ee=kk ; ff=ii ; gg=ll ; hh=jj ; endif !\Phi4_kilj
290 [ + + ]: 19874403 : if (itrans==12) then ; ee=kk ; ff=jj ; gg=ll ; hh=ii ; endif !\Phi4_kjli
291 : :
292 [ + + ]: 19874403 : if (itrans==13) then ; ee=ii ; ff=ll ; gg=jj ; hh=kk ; endif !\Phi4_iljk
293 [ + + ]: 19874403 : if (itrans==14) then ; ee=ii ; ff=ll ; gg=kk ; hh=jj ; endif !\Phi4_ilkj
294 [ + + ]: 19874403 : if (itrans==15) then ; ee=jj ; ff=ll ; gg=ii ; hh=kk ; endif !\Phi4_jlik
295 [ + + ]: 19874403 : if (itrans==16) then ; ee=jj ; ff=ll ; gg=kk ; hh=ii ; endif !\Phi4_jlki
296 [ + + ]: 19874403 : if (itrans==17) then ; ee=kk ; ff=ll ; gg=ii ; hh=jj ; endif !\Phi4_klij
297 [ + + ]: 19874403 : if (itrans==18) then ; ee=kk ; ff=ll ; gg=jj ; hh=ii ; endif !\Phi4_klji
298 : :
299 [ + + ]: 19874403 : if (itrans==19) then ; ee=ll ; ff=ii ; gg=jj ; hh=kk ; endif !\Phi4_lijk
300 [ + + ]: 19874403 : if (itrans==20) then ; ee=ll ; ff=ii ; gg=kk ; hh=jj ; endif !\Phi4_likj
301 [ + + ]: 19874403 : if (itrans==21) then ; ee=ll ; ff=jj ; gg=ii ; hh=kk ; endif !\Phi4_ljik
302 [ + + ]: 19874403 : if (itrans==22) then ; ee=ll ; ff=jj ; gg=kk ; hh=ii ; endif !\Phi4_ljki
303 [ + + ]: 19874403 : if (itrans==23) then ; ee=ll ; ff=kk ; gg=ii ; hh=jj ; endif !\Phi4_lkij
304 [ + + ]: 19874403 : if (itrans==24) then ; ee=ll ; ff=kk ; gg=jj ; hh=ii ; endif !\Phi4_lkji
305 : :
306 : 26499204 : Phi4_3333(ee,ff,gg,hh)=Phi4_tmp(ii,jj,kk,ll)
307 : : end do
308 : : end do
309 : : end do
310 : : end do
311 : :
312 : 245363 : end subroutine tdep_build_phi4_3333
313 : :
314 : : !=====================================================================================================
315 : : end module m_tdep_phi4
|
