1+ module extended_geom
2+ use constants, only: real12
3+ use misc_linalg, only: modu, LUinv, cross
4+ use rw_geom, only: basis_type, species_type
5+ implicit none
6+
7+
8+ private
9+
10+ public :: extended_basis_type
11+
12+
13+ type, extends(basis_type) :: extended_basis_type
14+ real (real12) :: max_extension
15+ integer :: nspec_images
16+ integer :: num_images
17+ type (species_type), dimension (:), allocatable :: image_species
18+ contains
19+ procedure , pass(this) :: create_images
20+ end type extended_basis_type
21+
22+ contains
23+
24+ subroutine create_images (this , max_bondlength )
25+ class(extended_basis_type), intent (inout ) :: this
26+ real (real12), intent (in ) :: max_bondlength
27+
28+ type (species_type), dimension (this% nspec) :: image_species
29+
30+ integer :: is, ia, js, i, j, k
31+ integer :: amax, bmax, cmax
32+ real (real12), dimension (3 ) :: vtmp1
33+
34+
35+ !- --------------------------------------------------------------------------
36+ ! get the maximum number of lattice vectors to consider
37+ ! NOTE: this is not perfect
38+ ! won't work for extremely acute/obtuse angle cells
39+ ! (due to diagonal path being shorter than individual lattice vectors)
40+ !- --------------------------------------------------------------------------
41+ amax = ceiling (max_bondlength/ modu(this% lat(1 ,:)))
42+ bmax = ceiling (max_bondlength/ modu(this% lat(2 ,:)))
43+ cmax = ceiling (max_bondlength/ modu(this% lat(3 ,:)))
44+
45+
46+ ! !! WARNING, NEED IGNORE LIST IN HERE TO ONLY APPLY TO ATOMS WE WANT TO EXTEND !!!
47+ ! !! needs and update_images subroutine !!!
48+ spec_loop1: do is= 1 ,this% nspec
49+ allocate ( image_species(is)% atom( &
50+ this% spec(is)% num* (2 * amax+1 )* (2 * bmax+1 )* (2 * cmax+1 ), &
51+ size (this% spec(is)% atom,2 ) ) &
52+ )
53+ image_species(is)% num = 0
54+ image_species(is)% mass = this% spec(is)% mass
55+ image_species(is)% charge = this% spec(is)% charge
56+ image_species(is)% radius = this% spec(is)% radius
57+ image_species(is)% name = this% spec(is)% name
58+ atom_loop1: do ia= 1 ,this% spec(is)% num
59+ do i=- amax,amax+1 ,1
60+ vtmp1(1 ) = this% spec(is)% atom(ia,1 ) + real (i, real12)
61+ do j=- bmax,bmax+1 ,1
62+ vtmp1(2 ) = this% spec(is)% atom(ia,2 ) + real (j, real12)
63+ do k=- cmax,cmax+1 ,1
64+ vtmp1(3 ) = this% spec(is)% atom(ia,3 ) + real (k, real12)
65+ if ( get_distance_from_unit_cell(vtmp1, this% lat) .le. max_bondlength ) then
66+ ! add the image to the list
67+ image_species(is)% num = this% spec(is)% num + 1
68+ this% image_species(is)% atom(image_species(is)% num,:3 ) = vtmp1
69+ end if
70+ end do
71+ end do
72+ end do
73+ end do atom_loop1
74+ end do spec_loop1
75+
76+
77+ this% nspec_images = count (image_species% num.gt. 0 )
78+ allocate (this% image_species(this% nspec_images))
79+ js = 0
80+ do is = 1 , this% nspec
81+ js = js + 1
82+ this% image_species(js)% num = image_species(is)% num
83+ this% image_species(js)% mass = image_species(is)% mass
84+ this% image_species(js)% charge = image_species(is)% charge
85+ this% image_species(js)% radius = image_species(is)% radius
86+ this% image_species(js)% name = image_species(is)% name
87+ allocate (this% image_species(js)% atom(image_species(is)% num,size (image_species(is)% atom,2 )))
88+ this% image_species(js)% atom(:,:) = image_species(is)% atom(:this% nspec_images,:)
89+ end do
90+
91+ end subroutine create_images
92+
93+
94+
95+
96+
97+ function get_distance_from_unit_cell (point , lattice , closest_point , is_cartesian ) result(distance)
98+ implicit none
99+ ! Input
100+ real (real12), intent (in ) :: point(3 ) ! Point in 3D space
101+ real (real12), intent (in ) :: lattice(3 ,3 ) ! 3x3 array representing lattice vectors as columns
102+ ! Output
103+ real (real12), intent (out ), optional :: closest_point(3 ) ! Closest point on the unit cell surface
104+ logical , optional , intent (in ) :: is_cartesian ! Flag indicating whether the point is in Cartesian coordinates
105+ real (real12) :: distance ! Shortest distance to the unit cell surface
106+ ! Local variables
107+ real (real12), dimension (3 ) :: point_
108+ real (real12), dimension (3 ,3 ) :: inverse_lattice
109+ real (real12), dimension (3 ) :: normal
110+ real (real12), dimension (3 ) :: plane_point
111+ real (real12), dimension (3 ) :: projection, closest_point_
112+ real (real12), dimension (3 ) :: inverse_projection
113+ real (real12) :: min_distance
114+ logical :: is_outside = .false.
115+ integer :: i, j, k
116+ integer , dimension (3 ) :: index_list = [1 , 2 , 3 ]
117+ logical :: is_cartesian_ = .false.
118+
119+
120+
121+ if (present (is_cartesian)) is_cartesian_ = is_cartesian
122+
123+ inverse_lattice = LUinv(lattice)
124+ if (is_cartesian_) then
125+ ! Convert point to fractional coordinates
126+ ! point_ = matmul(LUinv(lattice), point)
127+ point_ = point
128+ else
129+ point_ = matmul (lattice, point)
130+ end if
131+
132+ min_distance = huge (1._real12 )
133+
134+ ! get projection of point onto each face of the lattice
135+ ! get the length of the projection vector
136+ ! if negative, then the point is inside the unit cell
137+ ! if positive, then the point is outside the unit cell
138+ ! if the projection falls outside of the cell edges, use edge or corner distances
139+
140+ do i = 1 , 3
141+ index_list = cshift (index_list, 1 )
142+ plane_point = 0._real12
143+ do j = 1 , 2
144+ normal = (- 1._real12 )** j * cross(lattice(index_list(2 ),:), lattice(index_list(3 ),:))
145+ projection = project_point_onto_plane(point_, plane_point, normal)
146+
147+ ! check if point minus projection is negative
148+ ! if so, it is on the wrong side of the plane and should be ignored
149+ if ( dot_product (point_ - projection, normal) .lt. 0._real12 ) cycle
150+ is_outside = .true.
151+
152+ ! check if projection is outside the surface
153+
154+ inverse_projection = matmul (inverse_lattice, projection)
155+ if ( any ( inverse_projection .lt. 0._real12 ) .or. &
156+ any ( inverse_projection .gt. 1._real12 ) ) then
157+ ! projection is outside the surface
158+ ! check if the projection is outside the edges
159+ ! if it is, then the closest point is the edge or corner
160+ ! if it is not, then the closest point is the projection
161+ do k = 1 , 3
162+ if ( inverse_projection(k) .lt. 0._real12 ) then
163+ inverse_projection(k) = 0._real12
164+ else if ( inverse_projection(k) .gt. 1._real12 ) then
165+ inverse_projection(k) = 1._real12
166+ end if
167+ end do
168+ projection = matmul (lattice, inverse_projection)
169+ end if
170+ distance = norm2(point_ - projection)
171+ if ( distance .lt. min_distance ) then
172+ min_distance = distance
173+ closest_point_ = projection
174+ end if
175+
176+ ! ! makes it apply to the next iteration
177+ plane_point = plane_point + lattice(index_list(1 ),:)
178+ end do
179+ end do
180+
181+ if ( is_outside ) then
182+ distance = min_distance
183+ else
184+ distance = 0._real12
185+ end if
186+
187+ if ( present (closest_point) ) then
188+ if (is_cartesian_) then
189+ closest_point = closest_point_
190+ else
191+ closest_point = matmul (inverse_lattice, closest_point_)
192+ end if
193+ end if
194+
195+ end function get_distance_from_unit_cell
196+
197+
198+
199+
200+ function project_point_onto_plane (point , plane_point , normal ) result(output)
201+ implicit none
202+ real (real12), dimension (3 ), intent (in ) :: point
203+ real (real12), dimension (3 ), intent (in ) :: plane_point
204+ real (real12), dimension (3 ), intent (in ) :: normal
205+ real (real12), dimension (3 ) :: output
206+
207+ real (real12) :: distance
208+ real (real12), dimension (3 ) :: vector_to_plane
209+
210+ vector_to_plane = point - plane_point
211+
212+ distance = dot_product (vector_to_plane, normal) / dot_product (normal, normal)
213+
214+ output = point - distance * normal
215+
216+ end function project_point_onto_plane
217+
218+
219+ end module extended_geom
0 commit comments