atomes 1.1.14
atomes: an atomic scale modeling tool box
Loading...
Searching...
No Matches
lattice.F90
Go to the documentation of this file.
1! This file is part of the 'atomes' software.
2!
3! 'atomes' is free software: you can redistribute it and/or modify it under the terms
4! of the GNU Affero General Public License as published by the Free Software Foundation,
5! either version 3 of the License, or (at your option) any later version.
6!
7! 'atomes' is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
8! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
9! See the GNU General Public License for more details.
10!
11! You should have received a copy of the GNU Affero General Public License along with 'atomes'.
12! If not, see <https://www.gnu.org/licenses/>
13!
14! Copyright (C) 2022-2024 by CNRS and University of Strasbourg
15!
20
21INTEGER (KIND=c_int) FUNCTION add_cells (NP, NPS, sizec) bind (C,NAME='add_cells_')
22
23USE parameters
24
25IMPLICIT NONE
26
27INTEGER (KIND=c_int), INTENT(IN) :: np, nps
28INTEGER (KIND=c_int), INTENT(IN), DIMENSION(3) :: sizec
29INTEGER :: pia, pib, pic, pid, pie, pif
30DOUBLE PRECISION, DIMENSION(3) :: lshift
31INTEGER, DIMENSION(:), ALLOCATABLE :: newlot
32DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: newpos
33
34INTERFACE
35 INTEGER FUNCTION send_pos(NPA, NPS, NLOT, POSTAB)
36 INTEGER, INTENT(IN) :: npa, nps
37 INTEGER, DIMENSION(:), INTENT(IN) :: nlot
38 DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: postab
39 END FUNCTION
40END INTERFACE
41
42pia = (sizec(1)+1)*(sizec(2)+1)*(sizec(3)+1)
43pib = np * pia
44
45if (allocated(newpos)) deallocate(newpos)
46allocate(newpos(pib,3,nps), stat=err)
47if (err .ne. 0) then
48 call show_error ("Impossible to allocate memory"//char(0), &
49 "Function: add_cells"//char(0), "Table: NEWPOS"//char(0))
50 add_cells=0
51 goto 001
52endif
53if (allocated(newlot)) deallocate(newlot)
54allocate(newlot(pib), stat=err)
55if (err .ne. 0) then
56 call show_error ("Impossible to allocate memory"//char(0), &
57 "Function: add_cells"//char(0), "Table: NEWLOT"//char(0))
58 add_cells=0
59 goto 001
60endif
61
62do pia=1, nps
63 pib=0
64 do pid=1, sizec(1)+1
65 do pie=1, sizec(2)+1
66 do pif=1, sizec(3)+1
67 lshift(1)=(pid-1)*the_box(1)%lvect(1,1) + (pie-1)*the_box(1)%lvect(2,1) + (pif-1)*the_box(1)%lvect(3,1)
68 lshift(2)=(pid-1)*the_box(1)%lvect(1,2) + (pie-1)*the_box(1)%lvect(2,2) + (pif-1)*the_box(1)%lvect(3,2)
69 lshift(3)=(pid-1)*the_box(1)%lvect(1,3) + (pie-1)*the_box(1)%lvect(2,3) + (pif-1)*the_box(1)%lvect(3,3)
70 do pic=1, np
71 pib=pib+1
72 newpos(pib,:,pia) = fullpos(pic,:,pia) + lshift(:)
73 newlot(pib) = lot(pic)
74 enddo
75 enddo
76 enddo
77 enddo
78enddo
79
80add_cells=0
81call init_data (pib, nsp, nps, 0)
82if (send_pos(pib, nps, newlot, newpos) .eq. 1) add_cells=1
83
84001 continue
85if (allocated(newpos)) deallocate(newpos)
86if (allocated(newlot)) deallocate(newlot)
87
88END FUNCTION
89
90INTEGER (KIND=c_int) FUNCTION shift_box_center (NP, NPS, cshift, REF) bind (C,NAME='shift_box_center_')
91
92USE parameters
93
94IMPLICIT NONE
95
96INTEGER (KIND=c_int), INTENT(IN) :: np, nps, ref
97real(kind=c_double), INTENT(IN), DIMENSION(3) :: cshift
98INTEGER :: pib, pic, pid
99DOUBLE PRECISION, DIMENSION(3,3) :: h_mat
100DOUBLE PRECISION, DIMENSION(3) :: tpo
101
102INTERFACE
103 INTEGER FUNCTION send_pos(NPA, NPS, NLOT, POSTAB)
104 INTEGER, INTENT(IN) :: npa, nps
105 INTEGER, DIMENSION(:), INTENT(IN) :: nlot
106 DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(IN) :: postab
107 END FUNCTION
108END INTERFACE
109
110h_mat(:,1) = the_box(1)%lvect(1,:)
111h_mat(:,2) = the_box(1)%lvect(2,:)
112h_mat(:,3) = the_box(1)%lvect(3,:)
113
114do pic=1, nps
115 do pib=1, np
116 do pid=1, 3
117 fullpos(pib,pid,pic) = fullpos(pib,pid,pic) + cshift(pid)
118 enddo
119 tpo=matmul(the_box(1)%lrecp, fullpos(pib,:,pic))
120 tpo=tpo-nint(tpo/0.5)
121 fullpos(pib,:,pic) = matmul(h_mat,tpo)
122
123 !FULLPOS(PIB,:,PIC) = FULLPOS(PIB,:,PIC) + cshift(2)
124 !FULLPOS(PIB,:,PIC) = FULLPOS(PIB,:,PIC) + cshift(3)
125 enddo
126enddo
127
128if (ref .eq. 1) then
129 shift_box_center = send_pos(np, nps, lot, fullpos)
130else
131 shift_box_center = 1
132endif
133
134END FUNCTION
135
136DOUBLE PRECISION FUNCTION f_dot_product (a, b)
137
138DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a, b
139INTEGER :: dim
140f_dot_product = 0.0d0
141
142do dim=1, 3
143 f_dot_product = f_dot_product + a(dim)*b(dim)
144enddo
145
146END FUNCTION
147
148SUBROUTINE f_cross_product (a, b, c)
149
150DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a, b
151DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: c
152
153c(:)=0.0d0
154
155c(1) = a(2)*b(3) - a(3)*b(2)
156c(2) = a(3)*b(1) - a(1)*b(3)
157c(3) = a(1)*b(2) - a(2)*b(1)
158
159END SUBROUTINE
160
161INTEGER (KIND=c_int) FUNCTION lattice (totl, lid, vectors, vmod, angles, lat, cfrac, apbc) bind (C,NAME='lattice_')
162
163!
164! Compute lattice angles from lattice vectors
165! Lattice vector modules
166! Lattice volume
167! Reciprocal lattice parameters
168!
169
170USE parameters
171
172IMPLICIT NONE
173
174real(kind=c_double), INTENT(IN), DIMENSION(3,3) :: vectors
175real(kind=c_double), INTENT(IN), DIMENSION(3) :: vmod
176real(kind=c_double), INTENT(INOUT), DIMENSION(3) :: angles
177INTEGER (KIND=c_int), INTENT(IN) :: totl, lid, lat, cfrac, apbc
178
179DOUBLE PRECISION :: alpha, beta, gama ! Lattice Angles
180DOUBLE PRECISION :: calpha, salpha, cbeta, sbeta, cgama, sgama ! Cosinus and Sinus
181DOUBLE PRECISION, DIMENSION(3) :: tmpla
182
183INTERFACE
184 DOUBLE PRECISION FUNCTION f_dot_product (a, b)
185 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a, b
186 END FUNCTION
187 SUBROUTINE f_cross_product (a, b, c)
188 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: a, b
189 DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT) :: c
190 END SUBROUTINE
191END INTERFACE
192
193! Transition from C/Gtk to Fortran90 !
194lattice = 0
195
196if (lid .eq. 0) then
197 if (allocated(the_box)) deallocate(the_box)
198 allocate(the_box(totl), stat=err)
199 if (err .ne. 0) then
200 call show_error ("Impossible to allocate memory"//char(0), &
201 "Function: lattice_"//char(0), "Type: THE_BOX"//char(0))
202 goto 001
203 endif
204 ncells = totl
205endif
206
207nbox => the_box(lid+1)
208
209nbox%GLASS=.false.
210nbox%CUBIC=.false.
211
212if (lat .gt. 0) then
213
214 do i=1, 3
215 do j=1, 3
216 nbox%lvect(j,i) = vectors(i,j)
217 enddo
218 enddo
219
220 nbox%modv(1) = vmod(1)
221 nbox%modv(2) = vmod(2)
222 nbox%modv(3) = vmod(3)
223 alpha = angles(1)
224 beta = angles(2)
225 gama = angles(3)
226
227 if (lat .eq. 1) then
228
229 if (alpha.eq.90.0 .and. beta.eq.90.0 .and. gama.eq.90.0) then
230 nbox%GLASS=.true.
231 if (nbox%modv(1).eq.nbox%modv(2) .and. nbox%modv(2).eq.nbox%modv(3)) nbox%CUBIC=.true.
232 endif
233
234 if (alpha.eq.90.0) then
235 alpha = pi/2.0d0
236 salpha = 1.0d0
237 calpha = 0.0d0
238 else
239 alpha = alpha*pi/180.0d0
240 salpha = sin(alpha)
241 calpha = cos(alpha)
242 endif
243 if (beta.eq.90.0) then
244 beta = pi/2.0d0
245 sbeta = 1.0d0
246 cbeta = 0.0d0
247 else
248 beta = beta*pi/180.0d0
249 sbeta = sin(beta)
250 cbeta = cos(beta)
251 endif
252 if (gama.eq.90.0) then
253 gama = pi/2.0d0
254 sgama = 1.0d0
255 cgama = 0.0d0
256 else
257 gama = gama*pi/180.0d0
258 sgama = sin(gama)
259 cgama = cos(gama)
260 endif
261
262 nbox%lvect(1,1) = nbox%modv(1)
263 nbox%lvect(1,2) = 0.0d0
264 nbox%lvect(1,3) = 0.0d0
265 nbox%lvect(2,1) = nbox%modv(2)*cgama
266 nbox%lvect(2,2) = nbox%modv(2)*sgama
267 nbox%lvect(2,3) = 0.d0
268 nbox%lvect(3,1) = nbox%modv(3)*cbeta
269 ltemp = (calpha - cbeta*cgama)/sgama
270 nbox%lvect(3,2) = nbox%modv(3)*ltemp
271 nbox%lvect(3,3) = nbox%modv(3)*sqrt(sbeta*sbeta - ltemp*ltemp)
272
273 else
274
275 nbox%modv(1)=sqrt(nbox%lvect(1,1)**2+nbox%lvect(1,2)**2+nbox%lvect(1,3)**2)
276 nbox%modv(2)=sqrt(nbox%lvect(2,1)**2+nbox%lvect(2,2)**2+nbox%lvect(2,3)**2)
277 nbox%modv(3)=sqrt(nbox%lvect(3,1)**2+nbox%lvect(3,2)**2+nbox%lvect(3,3)**2)
278
279 alpha= (nbox%lvect(3,1)*nbox%lvect(2,1)+ &
280 nbox%lvect(3,2)*nbox%lvect(2,2)+ &
281 nbox%lvect(3,3)*nbox%lvect(2,3))/(nbox%modv(2)*nbox%modv(3))
282 beta = (nbox%lvect(1,1)*nbox%lvect(3,1)+ &
283 nbox%lvect(1,2)*nbox%lvect(3,2)+ &
284 nbox%lvect(1,3)*nbox%lvect(3,3))/(nbox%modv(1)*nbox%modv(3))
285 gama = (nbox%lvect(1,1)*nbox%lvect(2,1)+ &
286 nbox%lvect(1,2)*nbox%lvect(2,2)+ &
287 nbox%lvect(1,3)*nbox%lvect(2,3))/(nbox%modv(1)*nbox%modv(2))
288
289 if (alpha.eq.0.0d0 .and. beta.eq.0.0d0 .and. gama.eq.0.0d0) then
290 nbox%GLASS=.true.
291 if (nbox%modv(1).eq.nbox%modv(2) .and. nbox%modv(2).eq.nbox%modv(3)) nbox%CUBIC=.true.
292 endif
293
294 if (alpha.eq.0.0d0) then
295 angles(1) = 90.0d0
296 alpha = pi/2.0d0
297 salpha = 1.0d0
298 calpha = 0.0d0
299 else
300 alpha = acos(alpha)
301 angles(1) = alpha*180.0d0/pi
302 salpha = sin(alpha)
303 calpha = cos(alpha)
304 endif
305 if (beta.eq.0.0d0) then
306 angles(2) = 90.0d0
307 beta = pi/2.0d0
308 sbeta = 1.0d0
309 cbeta = 0.0d0
310 else
311 beta = acos(beta)
312 angles(2) = beta*180.0d0/pi
313 sbeta = sin(beta)
314 cbeta = cos(beta)
315 endif
316 if (gama.eq.0.0d0) then
317 angles(3) = 90.0d0
318 gama = pi/2.0d0
319 sgama = 1.0d0
320 cgama = 0.0d0
321 else
322 gama = acos(gama)
323 angles(3) = gama*180.0d0/pi
324 sgama = sin(gama)
325 cgama = cos(gama)
326 endif
327
328 endif
329
330 if (alpha.eq.0.0d0 .and. beta.eq.0.0d0 .and. gama.eq.0.0d0) then
331 ! If some problems display a GTK error message.
332 call show_error ("Problem with the simulation box parameters"//char(0), &
333 "Computed angles are equal to 0.0d0"//char(0), "Function: lattice"//char(0))
334 goto 001
335 endif
336
337! write (*,*) NBOX%lvect(1,1), NBOX%lvect(1,2), NBOX%lvect(1,3)
338! write (*,*) NBOX%lvect(2,1), NBOX%lvect(2,2), NBOX%lvect(2,3)
339! write (*,*) NBOX%lvect(3,1), NBOX%lvect(3,2), NBOX%lvect(3,3)
340! write (*,*) NBOX%modv(1), NBOX%modv(2), NBOX%modv(3)
341! write (*,*) ALPHA*180/PI, BETA*180/PI, GAMA*180/PI
342
343 nbox%minv=min(nbox%modv(1),nbox%modv(2))
344 nbox%minv=min(nbox%minv,nbox%modv(3))
345 nbox%maxv=max(nbox%modv(1),nbox%modv(2))
346 nbox%maxv=max(nbox%maxv,nbox%modv(3))
347
348 nbox%VOLUME=(nbox%lvect(1,2)*nbox%lvect(2,3)-nbox%lvect(1,3)*nbox%lvect(2,2))*nbox%lvect(3,1) &
349 +(nbox%lvect(1,3)*nbox%lvect(2,1)-nbox%lvect(1,1)*nbox%lvect(2,3))*nbox%lvect(3,2) &
350 +(nbox%lvect(1,1)*nbox%lvect(2,2)-nbox%lvect(1,2)*nbox%lvect(2,1))*nbox%lvect(3,3)
351 nbox%VOLUME = abs(nbox%VOLUME)
352
353! Reciprocal lattice parameters !
354
355 call f_cross_product (nbox%lvect(2,:), nbox%lvect(3,:), tmpla)
356 nbox%lrecp(1,:) = tmpla / f_dot_product(nbox%lvect(1,:), tmpla)
357 call f_cross_product (nbox%lvect(3,:), nbox%lvect(1,:), tmpla)
358 nbox%lrecp(2,:) = tmpla / f_dot_product(nbox%lvect(2,:), tmpla)
359 call f_cross_product (nbox%lvect(1,:), nbox%lvect(2,:), tmpla)
360 nbox%lrecp(3,:) = tmpla / f_dot_product(nbox%lvect(3,:), tmpla)
361
362! modules of the reciprocal lattice vectors
363
364 nbox%modr(1)=sqrt(nbox%lrecp(1,1)**2+nbox%lrecp(1,2)**2+nbox%lrecp(1,3)**2)
365 nbox%modr(2)=sqrt(nbox%lrecp(2,1)**2+nbox%lrecp(2,2)**2+nbox%lrecp(2,3)**2)
366 nbox%modr(3)=sqrt(nbox%lrecp(3,1)**2+nbox%lrecp(3,2)**2+nbox%lrecp(3,3)**2)
367
368 nbox%modr(:) = 2.0d0*pi*nbox%modr(:)
369
370 nbox%minr=min(nbox%modr(1),nbox%modr(2))
371 nbox%minr=min(nbox%minr,nbox%modr(3))
372 nbox%maxr=max(nbox%modr(1),nbox%modr(2))
373 nbox%maxr=max(nbox%maxr,nbox%modr(3))
374
375! Creation of the matrix to convert fractional to cartesian coordinates
376
377 z=sqrt(abs(1 - calpha*calpha &
378 - cbeta*cbeta &
379 - cgama*cgama &
380 + 2*calpha*cbeta*cgama))
381 z=z/sgama
382 i = 0
383 if (cfrac > 0) i = cfrac - 1
384
385 nbox%fractocart(1,1)=nbox%modv(1)/(2.0**(i))
386 nbox%fractocart(1,2)=0.0d0
387 nbox%fractocart(1,3)=0.0d0
388 nbox%fractocart(2,1)=nbox%modv(2)*cgama/(2.0**(i))
389 nbox%fractocart(2,2)=nbox%modv(2)*sgama/(2.0**(i))
390 nbox%fractocart(2,3)=0.0d0
391 nbox%fractocart(3,1)=nbox%modv(3)*cbeta/(2.0**(i))
392 nbox%fractocart(3,2)=nbox%modv(3)*((calpha-cbeta*cgama)/sgama)/(2.0**(i))
393 nbox%fractocart(3,3)=nbox%modv(3)*z/(2.0**(i))
394 !write (6, *)
395 !write (6, '("Frac to cart matrix in lattice:")')
396 !write (6, '(f15.10,4x,f15.10,4x,f15.10)') NBOX%fractocart(1,1), NBOX%fractocart(1,2), NBOX%fractocart(1,3)
397 !write (6, '(f15.10,4x,f15.10,4x,f15.10)') NBOX%fractocart(2,1), NBOX%fractocart(2,2), NBOX%fractocart(2,3)
398 !write (6, '(f15.10,4x,f15.10,4x,f15.10)') NBOX%fractocart(3,1), NBOX%fractocart(3,2), NBOX%fractocart(3,3)
399 !write (6, *)
400
401! Creation of the matrix to convert cartesian to fractional coordinates
402
403 nbox%carttofrac(1,1)=1.0d0/nbox%fractocart(1,1)
404 nbox%carttofrac(1,2)=0.0d0
405 nbox%carttofrac(1,3)=0.0d0
406 nbox%carttofrac(2,1)=-cgama/(sgama*(nbox%modv(1)/(2.0**(i))))
407 nbox%carttofrac(2,2)=1.0d0/nbox%fractocart(2,2)
408 nbox%carttofrac(2,3)=0.0d0
409 nbox%carttofrac(3,1)=((nbox%modv(2)/(2.0**(i)))*(nbox%modv(3)/(2.0**(i))))/nbox%VOLUME
410 nbox%carttofrac(3,1)=nbox%carttofrac(3,1) * (calpha*cgama - cbeta)/salpha
411 nbox%carttofrac(3,2)=((nbox%modv(1)/(2.0**(i)))*(nbox%modv(3)/(2.0**(i))))/nbox%VOLUME
412 nbox%carttofrac(3,2)=nbox%carttofrac(3,2) * (cbeta*cgama - calpha)/sgama
413 nbox%carttofrac(3,3)=1.0d0/nbox%fractocart(3,3)
414
415 if (apbc .eq. 1) then
416 pbc=.true.
417 else
418 pbc=.false.
419 nbox%GLASS=.false.
420 nbox%CUBIC=.false.
421 endif
422
423else
424
425 nbox%VOLUME=0.0d0
426 nbox%minr=0.0d0
427 nbox%modv(:)=0.0d0
428 nbox%modr(:)=0.0d0
429 nbox%lvect(:,:)=0.0d0
430 nbox%lrecp(:,:)=0.0d0
431 nbox%fractocart(:,:)=0.0d0
432 nbox%GLASS=.false.
433 nbox%CUBIC=.false.
434 pbc=.false.
435
436endif
437
438real_density=0.0d0
439total_density=0.0d0
440if (nbox%VOLUME.ne.0.0) then
441 do i=1, nsp
442 real_density=real_density+nbspbs(i)*mass(i)
443 enddo
444 total_density = dble(na)/nbox%VOLUME
445 real_density=real_density/avogadro
446 real_density=real_density/nbox%VOLUME
447 real_density=real_density*10.0d0
448endif
449
450! To lattice_info_
451call lattice_info (lid, nbox%VOLUME, real_density, &
452 nbox%lvect, nbox%lrecp, nbox%modv, angles, &
453 nbox%fractocart, nbox%carttofrac)
454
455if (lid .eq. totl-1) then
456 mbox = 0.0d0
457 meanvol = 0.0d0
458 overall_cubic = .true.
459 do i=1, ncells
460 mbox = mbox + the_box(i)%minv
461 meanvol = meanvol + the_box(i)%VOLUME
462 if (overall_cubic .and. the_box(i)%CUBIC) then
463 overall_cubic = .true.
464 else
465 overall_cubic = .false.
466 endif
467 enddo
468 meanvol = meanvol / ncells
469 mbox = mbox / (2.0d0*ncells)
470endif
471
472lattice = 1
473
474001 continue
475
476END FUNCTION
min
#define min(a, b)
Definition global.h:75
max
#define max(a, b)
Definition global.h:74
show_error
void show_error(char *error, int val, GtkWidget *win)
show error message
Definition interface.c:293
lattice
integer(kind=c_int) function lattice(totl, lid, vectors, vmod, angles, lat, cfrac, apbc)
Definition lattice.F90:162
shift_box_center
integer(kind=c_int) function shift_box_center(np, nps, cshift, ref)
Definition lattice.F90:91
f_dot_product
double precision function f_dot_product(a, b)
Definition lattice.F90:137
add_cells
integer(kind=c_int) function add_cells(np, nps, sizec)
Definition lattice.F90:22
f_cross_product
subroutine f_cross_product(a, b, c)
Definition lattice.F90:149
parameters::fullpos
double precision, dimension(:,:,:), allocatable fullpos
Definition parameters.F90:546
parameters::ncells
integer ncells
Definition parameters.F90:157
parameters::total_density
double precision total_density
Definition parameters.F90:408
parameters::avogadro
double precision, parameter avogadro
Definition parameters.F90:397
parameters::mass
double precision, dimension(:), allocatable mass
Definition parameters.F90:436
parameters::ltemp
double precision ltemp
Definition parameters.F90:385
parameters::real_density
double precision real_density
Definition parameters.F90:408
parameters::nbspbs
integer, dimension(:), allocatable nbspbs
Definition parameters.F90:187
parameters::nbox
type(lattice), pointer nbox
Definition parameters.F90:668
parameters::z
double precision z
Definition parameters.F90:395
parameters::err
integer err
Definition parameters.F90:123
parameters::the_box
type(lattice), dimension(:), allocatable, target the_box
Definition parameters.F90:667
parameters::meanvol
double precision meanvol
Definition parameters.F90:407
parameters::mbox
double precision mbox
Definition parameters.F90:401
parameters::overall_cubic
logical overall_cubic
Definition parameters.F90:88
parameters::lot
integer, dimension(:), allocatable lot
Definition parameters.F90:186
parameters::pbc
logical pbc
Definition parameters.F90:67
parameters::nsp
integer nsp
Definition parameters.F90:147
parameters::pi
double precision, parameter pi
Definition parameters.F90:396
send_pos
integer function send_pos(npa, nps, nlot, postab)
Definition prepdata.F90:157