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
|
