Skip to content

Commit

Permalink
Merge 72c991c into 9ea5456
Browse files Browse the repository at this point in the history
  • Loading branch information
@MSallermann
MSallermann authored Dec 18, 2018
2 parents 9ea5456 + 72c991c commit e0254ca
Show file tree
Hide file tree
Showing 6 changed files with 46 additions and 42 deletions.
26 changes: 12 additions & 14 deletions core/src/engine/Hamiltonian_Heisenberg.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -373,9 +373,7 @@ namespace Engine
if( jspin >= 0 )
{
Energy[ispin] -= 0.5 * mu_s[ispin] * mu_s[jspin] * mult / std::pow(ddi_magnitudes[i_pair], 3.0) *
(3 * spins[ispin].dot(ddi_normals[i_pair]) * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin].dot(spins[ispin]));
Energy[jspin] -= 0.5 * mu_s[ispin] * mu_s[jspin] * mult / std::pow(ddi_magnitudes[i_pair], 3.0) *
(3 * spins[ispin].dot(ddi_normals[i_pair]) * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin].dot(spins[ispin]));
(3 * spins[ispin].dot(ddi_normals[i_pair]) * spins[jspin].dot(ddi_normals[i_pair]) - spins[ispin].dot(spins[jspin]));
}
}
}
Expand Down Expand Up @@ -705,7 +703,6 @@ namespace Engine
if( jspin >= 0 )
{
gradient[ispin] -= mu_s[jspin] * skalar_contrib * (3 * ddi_normals[i_pair] * spins[jspin].dot(ddi_normals[i_pair]) - spins[jspin]);
gradient[jspin] -= mu_s[ispin] * skalar_contrib * (3 * ddi_normals[i_pair] * spins[ispin].dot(ddi_normals[i_pair]) - spins[ispin]);
}
}
}
Expand Down Expand Up @@ -1072,14 +1069,14 @@ namespace Engine

void Hamiltonian_Heisenberg::FFT_Dipole_Matrices(FFT::FFT_Plan & fft_plan_dipole, int img_a, int img_b, int img_c)
{
//prefactor of ddi interaction
// Prefactor of ddi interaction
scalar mult = C::mu_0 * C::mu_B * C::mu_B / ( 4*C::Pi * 1e-30 );

//size of original geometry
// Size of original geometry
int Na = geometry->n_cells[0];
int Nb = geometry->n_cells[1];
int Nc = geometry->n_cells[2];
//bravais vectors
// Bravais vectors
Vector3 ta = geometry->bravais_vectors[0];
Vector3 tb = geometry->bravais_vectors[1];
Vector3 tc = geometry->bravais_vectors[2];
Expand All @@ -1099,7 +1096,7 @@ namespace Engine
b_inter++;
inter_sublattice_lookup[i_b1 + i_b2 * geometry->n_cell_atoms] = b_inter;

//iterate over the padded system
// Iterate over the padded system
#pragma omp parallel for collapse(3)
for( int c = 0; c < n_cells_padded[2]; ++c )
{
Expand All @@ -1112,17 +1109,18 @@ namespace Engine
int c_idx = c < Nc ? c : c - n_cells_padded[2];
scalar Dxx = 0, Dxy = 0, Dxz = 0, Dyy = 0, Dyz = 0, Dzz = 0;
Vector3 diff;
//iterate over periodic images

// Iterate over periodic images
for( int a_pb = - img_a; a_pb <= img_a; a_pb++ )
{
for( int b_pb = - img_b; b_pb <= img_b; b_pb++ )
{
for( int c_pb = -img_c; c_pb <= img_c; c_pb++ )
{
diff = (a_idx + a_pb * Na + geometry->cell_atoms[i_b1][0] - geometry->cell_atoms[i_b2][0]) * ta
+ (b_idx + b_pb * Nb + geometry->cell_atoms[i_b1][1] - geometry->cell_atoms[i_b2][1]) * tb
+ (c_idx + c_pb * Nc + geometry->cell_atoms[i_b1][2] - geometry->cell_atoms[i_b2][2]) * tc;
diff = (a_idx + a_pb * Na) * ta
+ (b_idx + b_pb * Nb) * tb
+ (c_idx + c_pb * Nc) * tc
+ geometry->cell_atoms[i_b1]
- geometry->cell_atoms[i_b2];
if( diff.norm() > 1e-10 )
{
auto d = diff.norm();
Expand All @@ -1148,7 +1146,7 @@ namespace Engine
fft_dipole_inputs[idx + 4 * dipole_stride.comp] = Dyz;
fft_dipole_inputs[idx + 5 * dipole_stride.comp] = Dzz;

//We explicitly ignore the different strides etc. here
// We explicitly ignore the different strides etc. here
if( save_dipole_matrices && a < Na && b < Nb && c < Nc )
{
dipole_matrices[b_inter + n_inter_sublattice * (a + Na * (b + Nb * c))] << Dxx, Dxy, Dxz,
Expand Down
11 changes: 6 additions & 5 deletions core/src/engine/Hamiltonian_Heisenberg.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1073,11 +1073,12 @@ namespace Engine
{
for(int c_pb = -img[2]; c_pb <= img[2]; c_pb++)
{

diff = (a_idx + a_pb * n_cells[0] + cell_atoms[i_b1][0] - cell_atoms[i_b2][0]) * bravais_vectors[0]
+ (b_idx + b_pb * n_cells[1] + cell_atoms[i_b1][1] - cell_atoms[i_b2][1]) * bravais_vectors[1]
+ (c_idx + c_pb * n_cells[2] + cell_atoms[i_b1][2] - cell_atoms[i_b2][2]) * bravais_vectors[2];

diff = (a_idx + a_pb * n_cells[0]) * bravais_vectors[0]
+ (b_idx + b_pb * n_cells[1]) * bravais_vectors[1]
+ (c_idx + c_pb * n_cells[2]) * bravais_vectors[2]
+ cell_atoms[i_b1]
- cell_atoms[i_b2];

if(diff.norm() > 1e-10)
{
auto d = diff.norm();
Expand Down
35 changes: 23 additions & 12 deletions core/src/engine/Neighbours.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,8 @@ namespace Engine
{
auto pairs = pairfield(0);

if (radius > 1e-6)
// Check for a meaningful radius
if (std::abs(radius) > 1e-6)
{
Vector3 a = geometry.bravais_vectors[0];
Vector3 b = geometry.bravais_vectors[1];
Expand All @@ -157,36 +158,46 @@ namespace Engine

// This should give enough translations to contain all DDI pairs
int imax = 0, jmax = 0, kmax = 0;
if ( bounds_diff[0] > 0 )
imax = std::min(geometry.n_cells[0], (int)(1.1 * radius * geometry.n_cells[0] / bounds_diff[0]));
if ( bounds_diff[1] > 0 )
jmax = std::min(geometry.n_cells[1], (int)(1.1 * radius * geometry.n_cells[1] / bounds_diff[1]));
if ( bounds_diff[2] > 0 )
kmax = std::min(geometry.n_cells[2], (int)(1.1 * radius * geometry.n_cells[2] / bounds_diff[2]));

// If radius < 0 we take all pairs
if(radius > 0)
{
if ( bounds_diff[0] > 0 )
imax = std::min(geometry.n_cells[0] - 1, (int)(1.1 * radius * geometry.n_cells[0] / bounds_diff[0]));
if ( bounds_diff[1] > 0 )
jmax = std::min(geometry.n_cells[1] - 1, (int)(1.1 * radius * geometry.n_cells[1] / bounds_diff[1]));
if ( bounds_diff[2] > 0 )
kmax = std::min(geometry.n_cells[2] - 1, (int)(1.1 * radius * geometry.n_cells[2] / bounds_diff[2]));
} else {
imax = geometry.n_cells[0] - 1;
jmax = geometry.n_cells[1] - 1;
kmax = geometry.n_cells[2] - 1;
}

int i,j,k;
scalar dx;
Vector3 x0={0,0,0}, x1={0,0,0};

// Abort condidions for all 3 vectors
// Abort conditions for all 3 vectors
if (a.norm() == 0.0) imax = 0;
if (b.norm() == 0.0) jmax = 0;
if (c.norm() == 0.0) kmax = 0;

for (int iatom = 0; iatom < geometry.n_cell_atoms; ++iatom)
{
x0 = geometry.cell_atoms[iatom];
for (i = imax; i >= -imax; --i)

for (i = -imax; i <= imax; ++i)
{
for (j = jmax; j >= -jmax; --j)
for (j = -jmax; j <= jmax; ++j)
{
for (k = kmax; k >= -kmax; --k)
for (k = -kmax; k <= kmax; ++k)
{
for (int jatom = 0; jatom < geometry.n_cell_atoms; ++jatom)
{
x1 = geometry.cell_atoms[jatom] + i*a + j*b + k*c;
dx = (x0-x1).norm();
if (dx < radius)
if (dx < radius && dx > 1e-8) // Exclude self-interactions
{
pairs.push_back( {iatom, jatom, {i, j, k} } );
}
Expand Down
7 changes: 2 additions & 5 deletions core/src/engine/Vectormath.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,11 +99,8 @@ namespace Engine
std::vector<Data::vector2_t> basis_cell_points(geometry.n_cell_atoms + 3);
for(int i = 0; i < geometry.n_cell_atoms; i++)
{
for(int j=0; j<2; j++)
{
basis_cell_points[i].x = double(geometry.cell_atoms[i][j] * geometry.bravais_vectors[j][0]);
basis_cell_points[i].y = double(geometry.cell_atoms[i][j] * geometry.bravais_vectors[j][1]);
}
basis_cell_points[i].x = double(geometry.cell_atoms[i][0]);
basis_cell_points[i].y = double(geometry.cell_atoms[i][1]);
}

Vector3 basis_offset = geometry.cell_atoms[0][0] * geometry.bravais_vectors[0] + geometry.cell_atoms[0][1] * geometry.bravais_vectors[1];
Expand Down
7 changes: 2 additions & 5 deletions core/src/engine/Vectormath.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -297,11 +297,8 @@ namespace Engine
std::vector<Data::vector2_t> basis_cell_points(geometry.n_cell_atoms + 3);
for(int i = 0; i < geometry.n_cell_atoms; i++)
{
for(int j=0; j<2; j++)
{
basis_cell_points[i].x = double(geometry.cell_atoms[i][j] * geometry.bravais_vectors[j][0]);
basis_cell_points[i].y = double(geometry.cell_atoms[i][j] * geometry.bravais_vectors[j][1]);
}
basis_cell_points[i].x = double(geometry.cell_atoms[i][0]);
basis_cell_points[i].y = double(geometry.cell_atoms[i][1]);
}

Vector3 basis_offset = geometry.cell_atoms[0][0] * geometry.bravais_vectors[0] + geometry.cell_atoms[0][1] * geometry.bravais_vectors[1];
Expand Down
2 changes: 1 addition & 1 deletion core/test/input/physics_ddi.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ n_interaction_pairs 0
lattice_constant 1.0

### The bravais lattice type
bravais_lattice sc
bravais_lattice hex2d

### Number of basis cells along principal
### directions (a b c)
Expand Down

0 comments on commit e0254ca

Please sign in to comment.