-
Notifications
You must be signed in to change notification settings - Fork 7
/
atom.f90
203 lines (173 loc) · 7.41 KB
/
atom.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
module Class_atom
use Class_helper
use Constants
use mpi
implicit none
enum, bind(c) !> A or B site in graphene
enumerator :: A_site, B_site, no_site
end enum
type atom
real(8) :: m_phi !> azimuthal spin angle \f$\phi\f$
!> see german wikipedia, not english
real(8) :: m_theta !> polar spin angle \f$\theta\f$
!> see german wikipedia, not english
real(8), dimension(3) :: pos !> Position in RS in atomic units
integer(4) :: site_type !> A or B site
integer , allocatable :: neigh_idx(:) !> index of neighbour atom
real(8), allocatable :: neigh_conn(:,:) !> real space connection to neighbour.
integer(4), allocatable :: conn_type(:) !> type of connection
!> First index connection, second element of connection.
integer :: me, nProcs
contains
procedure :: get_m_cart => get_m_cart
procedure :: set_sphere => set_sphere
procedure :: set_m_cart => set_m_cart
procedure :: free_atm => free_atm
procedure :: compare_to_root => compare_to_root
end type atom
contains
subroutine free_atm(self)
implicit none
class(atom) :: self
if(allocated(self%neigh_idx)) deallocate(self%neigh_idx)
if(allocated(self%neigh_conn)) deallocate(self%neigh_conn)
end subroutine free_atm
function get_m_cart(self) result(coord)
implicit none
Class(atom), intent(in) :: self
real(8), dimension(3) :: coord
! assume r = 1
coord(1) = sin(self%m_theta) * cos(self%m_phi)
coord(2) = sin(self%m_theta) * sin(self%m_phi)
coord(3) = cos(self%m_theta)
end function get_m_cart
function init_ferro_z(p_pos, site) result(self)
implicit none
type(atom) :: self
real(8), intent(in) :: p_pos(3)
integer, optional :: site
integer :: ierr(2)
call MPI_Comm_size(MPI_COMM_WORLD, self%nProcs, ierr(1))
call MPI_Comm_rank(MPI_COMM_WORLD, self%me, ierr(2))
call check_ierr(ierr, self%me, "init_ferro_z call")
if(present(site)) then
self%site_type = site
else
self%site_type = no_site
endif
self%m_phi = 0d0
self%m_theta = 0d0
self%pos = p_pos
end function init_ferro_z
subroutine set_m_cart(self,x,y,z)
implicit none
class(atom) :: self
real(8), intent(in) :: x,y,z
if(abs(my_norm2([x,y,z]) - 1d0) > 1d-4) then
call error_msg("Spin not normed", abort=.True.)
endif
self%m_theta = acos(z / sqrt(x*x + y*y + z*z))
self%m_phi = atan2(y,x)
end subroutine set_m_cart
subroutine set_sphere(self, phi, theta)
implicit none
class(atom) :: self
real(8), intent(in) :: phi, theta
self%m_phi = phi
self%m_theta = theta
!write(*,*) "self%m_theta",self%m_theta
end subroutine set_sphere
function compare_to_root(self) result(success)
implicit none
class(atom) :: self
real(8) :: tmp, tmp_p(3)
integer :: ierr(10), tmp_i
integer, allocatable :: tmp_ivec(:)
integer(4) :: tmp_i4
integer(4), allocatable :: tmp_i4vec(:)
real(8), allocatable :: tmp_rmtx(:,:)
logical :: success
success = .True.
! compare angles
if(self%me == root) tmp = self%m_phi
call MPI_Bcast(tmp, 1, MPI_REAL8, root, MPI_COMM_WORLD, ierr(1))
if(abs(tmp - self%m_phi) > 1d-12) then
call error_msg("m_phi doesn't match", abort=.True.)
success = .False.
endif
if(self%me == root) tmp = self%m_theta
call MPI_Bcast(tmp, 1, MPI_REAL8, root, MPI_COMM_WORLD, ierr(2))
if(abs(tmp - self%m_theta) > 1d-12) then
call error_msg("m_theta doesn't match", abort=.True.)
success = .False.
endif
! compare positions
if(self%me == root) tmp_p = self%pos
call MPI_Bcast(tmp_p, 3, MPI_REAL8, root, MPI_COMM_WORLD, ierr(3))
if(my_norm2(tmp_p - self%pos) > 1d-11) then
call error_msg("pos doesn't match", abort=.True.)
success = .False.
endif
! compare site_types
if(self%me == root) tmp_i4 = self%site_type
call MPI_Bcast(tmp_i4, 1, MPI_INTEGER4, root, MPI_COMM_WORLD, ierr(4))
if(tmp_i4 /= self%site_type) then
call error_msg("site_type doesn't match", abort=.True.)
success = .False.
endif
! compare neighbours
if(self%me == root) tmp_i = size(self%neigh_idx)
call MPI_Bcast(tmp_i, 1, MPI_INTEGER, root, MPI_COMM_WORLD, ierr(5))
if(tmp_i /= size(self%neigh_idx)) then
call error_msg("size(neigh_idx) doesn't match", abort=.True.)
success = .False.
endif
allocate(tmp_ivec(size(self%neigh_idx)))
if(self%me == root) tmp_ivec = self%neigh_idx
call MPI_Bcast(tmp_ivec, size(self%neigh_idx), MPI_INTEGER, &
root, MPI_COMM_WORLD, ierr(6))
if(any(tmp_ivec /= self%neigh_idx)) then
write (*,*) self%me, "neigh_idx", self%neigh_idx
write (*,*) self%me, "tmp_ivec", tmp_ivec
call error_msg("neigh_idx doesn't match", abort=.True.)
success = .False.
endif
deallocate(tmp_ivec)
if(self%me == root) tmp_i = size(self%neigh_conn)
call MPI_Bcast(tmp_i, 1, MPI_INTEGER, root, MPI_COMM_WORLD, ierr(7))
if(tmp_i /= size(self%neigh_conn)) then
call error_msg("size(neigh_conn) doesn't match", abort=.True.)
success = .False.
endif
allocate(tmp_rmtx(size(self%neigh_conn, dim=1), &
size(self%neigh_conn, dim=2)))
if(self%me == root) tmp_rmtx = self%neigh_conn
call MPI_Bcast(tmp_rmtx, size(tmp_rmtx), MPI_REAL8, &
root, MPI_COMM_WORLD, ierr(8))
if(mtx_norm(tmp_rmtx - self%neigh_conn) > 1d-11) then
call error_msg("neigh_conn doesn't match", abort=.True.)
success = .False.
endif
deallocate(tmp_rmtx)
! compare conn types
if(self%me == root) tmp_i = size(self%conn_type)
call MPI_Bcast(tmp_i, 1, MPI_INTEGER, root, MPI_COMM_WORLD, ierr(9))
if(tmp_i /= size(self%conn_type)) then
call error_msg("size(conn_type) doesn't match", abort=.True.)
success = .False.
endif
allocate(tmp_i4vec(size(self%conn_type)))
if(self%me == root) tmp_i4vec = self%conn_type
call MPI_Bcast(tmp_i4vec, size(tmp_i4vec), MPI_INTEGER4, &
root, MPI_COMM_WORLD, ierr(10))
if(any(tmp_i4vec /= self%conn_type)) then
call error_msg("conn_type doesn't match", abort=.True.)
success = .False.
endif
deallocate(tmp_i4vec)
if(.not. success) then
call error_msg("Atom test failed", abort=.True.)
endif
call check_ierr(ierr, self%me, "compare atoms")
end function compare_to_root
end module