newly added ordered dipol field calculations, still preliminary
This commit is contained in:
parent
78a674886f
commit
383c7d0abd
40
dipol_fields/ordered_dipol_fields/LSCO.inp
Normal file
40
dipol_fields/ordered_dipol_fields/LSCO.inp
Normal file
@ -0,0 +1,40 @@
|
|||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
%
|
||||||
|
% LSCO.inp
|
||||||
|
%
|
||||||
|
% input file to calculate magnetic dipol field distribution for an ordered state
|
||||||
|
%
|
||||||
|
% data from Japanese Journal of Applied Physics, Part 1 (1988) 27, 1132-1137
|
||||||
|
%
|
||||||
|
% muon site from H.-U. Suter, Physica B 326 (2003) 329
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
% magnetic lattice related input
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
4 % number of magnetic atoms in the magnetic unit cell
|
||||||
|
2 % number of magnetic types in the magnetic unit cell
|
||||||
|
|
||||||
|
0.0000 0.0000 0.0000 1 % x, y, z, magnetic type
|
||||||
|
0.5000 0.5000 0.0000 2 % x, y, z, magnetic type
|
||||||
|
0.0000 0.5000 0.5000 1 % x, y, z, magnetic type
|
||||||
|
0.5000 0.0000 0.5000 2 % x, y, z, magnetic type
|
||||||
|
|
||||||
|
5.3535 5.4003 13.1441 % lattice parameters in angstroms (magnetic cell)
|
||||||
|
90.0000 90.0000 90.0000 % lattice angles in degrees
|
||||||
|
|
||||||
|
200.0000 % Lorentz sphere radius in angstroms
|
||||||
|
|
||||||
|
0.6800 0.0000 0.0000 % magnetic moment for type 1
|
||||||
|
-0.6800 0.0000 0.0000 % magnetic moment for type 2
|
||||||
|
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
% what to be calculated specifically follows now
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
|
||||||
|
point % tag showing that a single muon position to be calculated
|
||||||
|
0.1478 0.0000 0.1846 % muon site in lattice coordinates
|
||||||
|
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
|
% end
|
||||||
|
%--------------------------------------------------------------------------------
|
||||||
28
dipol_fields/ordered_dipol_fields/Makefile
Normal file
28
dipol_fields/ordered_dipol_fields/Makefile
Normal file
@ -0,0 +1,28 @@
|
|||||||
|
#---------------------------------------------------------------------------
|
||||||
|
# Makefile for the dipol field calculations
|
||||||
|
#
|
||||||
|
# Author: Andreas Suter, 2007-07-11
|
||||||
|
#
|
||||||
|
# $Id$
|
||||||
|
#---------------------------------------------------------------------------
|
||||||
|
|
||||||
|
|
||||||
|
# compiler
|
||||||
|
CC = g++
|
||||||
|
INCPATH = -I/usr/lib/qt-3.3/mkspecs/default -I$(QTDIR)/include
|
||||||
|
CFLAGS = -Wall -O2 -g
|
||||||
|
LDFLAGS = -Wl,-rpath,/usr/lib
|
||||||
|
|
||||||
|
# libraries
|
||||||
|
LIBS = -L$(QTDIR)/lib -lqt-mt -lbsd -lm -lutil
|
||||||
|
|
||||||
|
all: dipolar_ordered
|
||||||
|
|
||||||
|
dipolar_ordered: dipolar_ordered.o
|
||||||
|
$(CC) $(LDFLAGS) -o dipolar_ordered dipolar_ordered.o $(LIBS)
|
||||||
|
|
||||||
|
dipolar_ordered.o: dipolar_ordered.cpp
|
||||||
|
$(CC) $(CFLAGS) $(INCPATH) -c -o dipolar_ordered.o dipolar_ordered.cpp
|
||||||
|
|
||||||
|
clean:
|
||||||
|
rm -rf *.o *~
|
||||||
5
dipol_fields/ordered_dipol_fields/README
Normal file
5
dipol_fields/ordered_dipol_fields/README
Normal file
@ -0,0 +1,5 @@
|
|||||||
|
Andreas Suter, 2007/09/20
|
||||||
|
|
||||||
|
This is still extremly experimental and not documented well, still it works.
|
||||||
|
|
||||||
|
More to come ...
|
||||||
633
dipol_fields/ordered_dipol_fields/dipolar_ordered.cpp
Normal file
633
dipol_fields/ordered_dipol_fields/dipolar_ordered.cpp
Normal file
@ -0,0 +1,633 @@
|
|||||||
|
#include <iostream>
|
||||||
|
using namespace std;
|
||||||
|
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <math.h>
|
||||||
|
|
||||||
|
#include <qfile.h>
|
||||||
|
#include <qtextstream.h>
|
||||||
|
#include <qstring.h>
|
||||||
|
#include <qvaluevector.h>
|
||||||
|
|
||||||
|
#define DIPOLAR_SUCCESS 1
|
||||||
|
#define DIPOLAR_FILE_NOT_FOUND -1
|
||||||
|
#define DIPOLAR_FILE_OPEN_FAILED -2
|
||||||
|
#define DIPOLAR_WRONG_INPUT_FILE_FORMAT -3
|
||||||
|
|
||||||
|
#define CALC_POINT 0
|
||||||
|
#define CALC_PLANE 1
|
||||||
|
#define CALC_SPHERE 2
|
||||||
|
|
||||||
|
// The half the number of unit cells used to calculate the dipolar field.
|
||||||
|
// This should be absorbed into the input file in the future.
|
||||||
|
#define DIPOLAR_NO_UNIT_CELLS 50
|
||||||
|
|
||||||
|
#define PI 3.141592653589793238462643383279502884197
|
||||||
|
|
||||||
|
// where is this number from? Check units and comment!
|
||||||
|
#define UNITS 9.2741
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
typedef struct {
|
||||||
|
int line_no;
|
||||||
|
QString data;
|
||||||
|
} input_data;
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
typedef struct {
|
||||||
|
double pos[3];
|
||||||
|
int type;
|
||||||
|
} atom_pos;
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
typedef struct {
|
||||||
|
double component[3];
|
||||||
|
} mag_moment;
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
typedef struct {
|
||||||
|
int type;
|
||||||
|
double dipoleMatrix[3][3];
|
||||||
|
double field[3];
|
||||||
|
} dipole_type_info;
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*
|
||||||
|
* <p>calcTag: 0=single muon position, 1=dipole field in a plain, 2=dipole field on a sphere
|
||||||
|
*/
|
||||||
|
typedef struct {
|
||||||
|
int calcTag; /// tag showing what to be calculated
|
||||||
|
int noAtoms;
|
||||||
|
int noMagTypes;
|
||||||
|
QValueVector<atom_pos> atom;
|
||||||
|
double latticeParam[3];
|
||||||
|
double angle[3];
|
||||||
|
double lorentzSphereRadius;
|
||||||
|
QValueVector<mag_moment> magMoment;
|
||||||
|
double muPos[3];
|
||||||
|
double planeFix[3]; /// vector fixing a plane in the crystal coordinate system
|
||||||
|
double planeNormal[3]; /// normal plane vector
|
||||||
|
QValueVector<dipole_type_info> dipole;
|
||||||
|
} dipole_info;
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
void dipolar_ordered_syntax()
|
||||||
|
{
|
||||||
|
cout << endl << "syntax: dipolar_ordered input_file output_file [debug]";
|
||||||
|
cout << endl << " input_file: input file name";
|
||||||
|
cout << endl << " The structure of the input file is:";
|
||||||
|
cout << endl << " details to be written ...";
|
||||||
|
cout << endl;
|
||||||
|
cout << endl << " WARNING: all the coordinates (including the magnetic moments) needed";
|
||||||
|
cout << endl << " to be given in crystal structure coordinates. This is trivial";
|
||||||
|
cout << endl << " for structures with lattice angles all 90°, but might lead to";
|
||||||
|
cout << endl << " errors/confusions for lattice angles not equal to 90°";
|
||||||
|
cout << endl << " a typical example would be:";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << " %";
|
||||||
|
cout << endl << " % LSCO.inp";
|
||||||
|
cout << endl << " %";
|
||||||
|
cout << endl << " % input file to calculate magnetic dipol field distribution for an ordered state";
|
||||||
|
cout << endl << " %";
|
||||||
|
cout << endl << " % data from Japanese Journal of Applied Physics, Part 1 (1988) 27, 1132-1137 ";
|
||||||
|
cout << endl << " % ";
|
||||||
|
cout << endl << " % muon site from H.-U. Suter, Physica B 326 (2003) 329";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << " % magnetic lattice related input";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << " 4 % number of magnetic atoms in the magnetic unit cell";
|
||||||
|
cout << endl << " 2 % number of magnetic types in the magnetic unit cell";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " 0.0000 0.0000 0.0000 1 % x, y, z, magnetic type";
|
||||||
|
cout << endl << " 0.5000 0.5000 0.0000 2 % x, y, z, magnetic type";
|
||||||
|
cout << endl << " 0.0000 0.5000 0.5000 1 % x, y, z, magnetic type";
|
||||||
|
cout << endl << " 0.5000 0.0000 0.5000 2 % x, y, z, magnetic type";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " 5.3535 5.4003 13.1441 % lattice parameters in angstroms (magnetic cell)";
|
||||||
|
cout << endl << " 90.0000 90.0000 90.0000 % lattice angles in degrees";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " 200.0000 % Lorentz sphere radius in angstroms";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " 0.6800 0.0000 0.0000 % magnetic moment for type 1";
|
||||||
|
cout << endl << " -0.6800 0.0000 0.0000 % magnetic moment for type 2";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << " % what to be calculated specifically follows now ";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " point % tag showing that a single muon position to be calculated";
|
||||||
|
cout << endl << " 0.1478 0.0000 0.1846 % muon site in lattice coordinates";
|
||||||
|
cout << endl << "";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << " % end";
|
||||||
|
cout << endl << " %--------------------------------------------------------------------------------";
|
||||||
|
cout << endl << " output_file: output file name";
|
||||||
|
cout << endl << endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
int dipolar_ordered_read_input_file(char *file_name, dipole_info *info)
|
||||||
|
{
|
||||||
|
QFile file(file_name);
|
||||||
|
input_data line_data;
|
||||||
|
QValueVector<input_data> data;
|
||||||
|
|
||||||
|
bool error = false;
|
||||||
|
char line[256];
|
||||||
|
QString qline;
|
||||||
|
int status;
|
||||||
|
int param_switch = 0;
|
||||||
|
int line_no = 0;
|
||||||
|
atom_pos atom;
|
||||||
|
mag_moment magMoment;
|
||||||
|
|
||||||
|
// check if the file exists
|
||||||
|
if (!file.exists())
|
||||||
|
return DIPOLAR_FILE_NOT_FOUND;
|
||||||
|
|
||||||
|
// open file
|
||||||
|
if (!file.open(IO_ReadOnly))
|
||||||
|
return DIPOLAR_FILE_OPEN_FAILED;
|
||||||
|
|
||||||
|
// read all the data
|
||||||
|
int comment_pos;
|
||||||
|
while (file.readLine(&line[0], sizeof(line)) && !file.atEnd()) {
|
||||||
|
qline = line;
|
||||||
|
line_no++;
|
||||||
|
// ignore empty and comment lines
|
||||||
|
if (qline.stripWhiteSpace().isEmpty() || qline.startsWith("%"))
|
||||||
|
continue;
|
||||||
|
// strip comments within data line
|
||||||
|
comment_pos = qline.find("%");
|
||||||
|
if (comment_pos>0)
|
||||||
|
qline = qline.remove(comment_pos-1, qline.length()-comment_pos);
|
||||||
|
qline = qline.lower();
|
||||||
|
line_data.line_no = line_no;
|
||||||
|
line_data.data = qline;
|
||||||
|
data.append(line_data);
|
||||||
|
}
|
||||||
|
|
||||||
|
file.close();
|
||||||
|
|
||||||
|
// analyze data
|
||||||
|
QValueVector<input_data>::Iterator it;
|
||||||
|
for (it = data.begin(); it != data.end(); ++it) {
|
||||||
|
switch (param_switch) {
|
||||||
|
case 0: // no of magnetic atoms in the magnetic unit cell
|
||||||
|
status = sscanf(it->data.latin1(), "%d", &info->noAtoms);
|
||||||
|
if (status != 1) // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 1: // no of magnetic types in the magnetic unit cell
|
||||||
|
status = sscanf(it->data.latin1(), "%d", &info->noMagTypes);
|
||||||
|
if (status != 1) // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 2: // get all the atomic positions and the associated magnetic type
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf %d",
|
||||||
|
&atom.pos[0], &atom.pos[1], &atom.pos[2], &atom.type);
|
||||||
|
if (status != 4) { // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
} else {
|
||||||
|
info->atom.append(atom);
|
||||||
|
}
|
||||||
|
// check if all atom position are read
|
||||||
|
if (info->atom.count() == (unsigned int)info->noAtoms) {
|
||||||
|
param_switch++;
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
case 3: // lattice parameters
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf",
|
||||||
|
&info->latticeParam[0], &info->latticeParam[1], &info->latticeParam[2]);
|
||||||
|
if (status != 3) // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 4: // lattice angles
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf",
|
||||||
|
&info->angle[0], &info->angle[1], &info->angle[2]);
|
||||||
|
if (status != 3) // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 5: // Lorentz sphere radius in (A)
|
||||||
|
status = sscanf(it->data.latin1(), "%lf", &info->lorentzSphereRadius);
|
||||||
|
if (status != 1) // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 6: // get all components of the magnetic types
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf",
|
||||||
|
&magMoment.component[0], &magMoment.component[1], &magMoment.component[2]);
|
||||||
|
if (status != 3) { // likely to have wrong input file format
|
||||||
|
error = true;
|
||||||
|
} else {
|
||||||
|
info->magMoment.append(magMoment);
|
||||||
|
}
|
||||||
|
// check if all magnetic moments are read
|
||||||
|
if (info->magMoment.count() == (unsigned int)info->noMagTypes) {
|
||||||
|
param_switch++;
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
case 7: // tag what to be calculated
|
||||||
|
if (it->data.contains("point") > 0) {
|
||||||
|
info->calcTag = CALC_POINT;
|
||||||
|
} else if (it->data.contains("plane") > 0) {
|
||||||
|
info->calcTag = CALC_PLANE;
|
||||||
|
} else if (it->data.contains("sphere") > 0) {
|
||||||
|
info->calcTag = CALC_SPHERE;
|
||||||
|
} else {
|
||||||
|
error = true;
|
||||||
|
}
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 8: // depends on the calculation tag
|
||||||
|
switch (info->calcTag) {
|
||||||
|
case CALC_POINT: // get the muon position
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf",
|
||||||
|
&info->muPos[0], &info->muPos[1], &info->muPos[2]);
|
||||||
|
if (status != 3)
|
||||||
|
error = true;
|
||||||
|
break;
|
||||||
|
case CALC_PLANE:
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf",
|
||||||
|
&info->planeFix[0], &info->planeFix[1], &info->planeFix[2]);
|
||||||
|
if (status != 3)
|
||||||
|
error = true;
|
||||||
|
break;
|
||||||
|
case CALC_SPHERE:
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
case 9: // depends on the calculation tag
|
||||||
|
switch (info->calcTag) {
|
||||||
|
case CALC_PLANE:
|
||||||
|
status = sscanf(it->data.latin1(), "%lf %lf %lf",
|
||||||
|
&info->planeNormal[0], &info->planeNormal[1], &info->planeNormal[2]);
|
||||||
|
if (status != 3)
|
||||||
|
error = true;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
param_switch++;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (error) {
|
||||||
|
cout << endl << "ERROR in line " << it->line_no;
|
||||||
|
return DIPOLAR_WRONG_INPUT_FILE_FORMAT;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return DIPOLAR_SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
void dipolar_ordered_point_calc(dipole_info *info)
|
||||||
|
{
|
||||||
|
double distance_vector[3];
|
||||||
|
double direction_vector[3];
|
||||||
|
double distance; ///< distance from an atom to the muon
|
||||||
|
double distance2; ///< squared distance form an atom to the muon
|
||||||
|
double idist3; ///< 1/distance^3
|
||||||
|
double dcos[3]; ///< directional cosines
|
||||||
|
|
||||||
|
// init dipole vector structure
|
||||||
|
dipole_type_info dipole;
|
||||||
|
for (int ii=0; ii<3; ii++) {
|
||||||
|
for (int jj=0; jj<3; jj++) {
|
||||||
|
dipole.dipoleMatrix[ii][jj] = 0.0; // init dipole tensor
|
||||||
|
}
|
||||||
|
dipole.field[ii] = 0.0; // init dipole field
|
||||||
|
}
|
||||||
|
for (int i=0; i<info->noMagTypes; i++) {
|
||||||
|
dipole.type = i+1; // init magnetic type
|
||||||
|
info->dipole.append(dipole); // append to the vector
|
||||||
|
}
|
||||||
|
|
||||||
|
// calculate the directional cosines
|
||||||
|
for (int i=0; i<3; i++) {
|
||||||
|
dcos[i] = cos(info->angle[i]*PI/180.0);
|
||||||
|
}
|
||||||
|
|
||||||
|
// triple loop over the x, y, z unit cells
|
||||||
|
for (int uc_x = -DIPOLAR_NO_UNIT_CELLS; uc_x <= DIPOLAR_NO_UNIT_CELLS; uc_x++) {
|
||||||
|
for (int uc_y = -DIPOLAR_NO_UNIT_CELLS; uc_y <= DIPOLAR_NO_UNIT_CELLS; uc_y++) {
|
||||||
|
for (int uc_z = -DIPOLAR_NO_UNIT_CELLS; uc_z <= DIPOLAR_NO_UNIT_CELLS; uc_z++) {
|
||||||
|
// iterate all atoms in the unit cell
|
||||||
|
QValueVector<atom_pos>::iterator it;
|
||||||
|
for (it = info->atom.begin(); it != info->atom.end(); ++it) {
|
||||||
|
|
||||||
|
// calculate the distance vector of the magnetic atom to the muon
|
||||||
|
distance_vector[0] = info->latticeParam[0]*((it->pos[0]+uc_x) - info->muPos[0]); // "x" component
|
||||||
|
distance_vector[1] = info->latticeParam[1]*((it->pos[1]+uc_y) - info->muPos[1]); // "y" component
|
||||||
|
distance_vector[2] = info->latticeParam[2]*((it->pos[2]+uc_z) - info->muPos[2]); // "z" component
|
||||||
|
|
||||||
|
// calculate the real distance between the atom and the muon
|
||||||
|
distance2 = distance_vector[0]*distance_vector[0]+
|
||||||
|
distance_vector[1]*distance_vector[1]+
|
||||||
|
distance_vector[2]*distance_vector[2]+
|
||||||
|
distance_vector[0]*distance_vector[1]*dcos[0]+
|
||||||
|
distance_vector[0]*distance_vector[2]*dcos[1]+
|
||||||
|
distance_vector[1]*distance_vector[2]*dcos[2];
|
||||||
|
distance = sqrt(distance2);
|
||||||
|
|
||||||
|
// calculate the direction vector
|
||||||
|
direction_vector[0] = distance_vector[0]/distance;
|
||||||
|
direction_vector[1] = distance_vector[1]/distance;
|
||||||
|
direction_vector[2] = distance_vector[2]/distance;
|
||||||
|
|
||||||
|
// check if it is within the Lorentz sphere, if yes calculate the dipole field tensor
|
||||||
|
if (info->lorentzSphereRadius - distance > 0.0) {
|
||||||
|
idist3 = 1.0/(distance*distance2);
|
||||||
|
for (int ii=0; ii<3; ii++) {
|
||||||
|
for (int jj=0; jj<3; jj++) {
|
||||||
|
if (ii==jj) { // diagonal elements
|
||||||
|
info->dipole.at(it->type-1).dipoleMatrix[ii][jj] += idist3*(3*direction_vector[ii]*direction_vector[ii]-1);
|
||||||
|
} else { // off diagonal elements
|
||||||
|
info->dipole.at(it->type-1).dipoleMatrix[ii][jj] += idist3*(3*direction_vector[ii]*direction_vector[jj]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// calculate the field
|
||||||
|
for (int i=0; i<info->noMagTypes; i++) { // calculation for each magnetic type in the lattice
|
||||||
|
for (int k=0; k<3; k++) {
|
||||||
|
for (int l=0; l<3; l++) {
|
||||||
|
info->dipole.at(i).field[k] += info->dipole.at(i).dipoleMatrix[k][l]*info->magMoment.at(i).component[l];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
void dipolar_ordered_plane_calc(dipole_info *info)
|
||||||
|
{
|
||||||
|
cout << endl << "SORRY: not yet implemented ..." << endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
void dipolar_ordered_sphere_calc(dipole_info *info)
|
||||||
|
{
|
||||||
|
cout << endl << "SORRY: not yet implemented ..." << endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
int dipolar_ordered_write_output_file(char *input_file_name, char *output_file_name, dipole_info *info)
|
||||||
|
{
|
||||||
|
char str[128];
|
||||||
|
// copy the input file to the output file
|
||||||
|
sprintf(str, "cp %s %s", input_file_name, output_file_name);
|
||||||
|
system(str);
|
||||||
|
|
||||||
|
// open the new output file and append the results
|
||||||
|
QFile file(output_file_name);
|
||||||
|
|
||||||
|
if (!file.open(IO_WriteOnly | IO_Append))
|
||||||
|
return DIPOLAR_FILE_OPEN_FAILED;
|
||||||
|
|
||||||
|
QTextStream fout(&file);
|
||||||
|
|
||||||
|
double field_tot[3];
|
||||||
|
for (int i=0; i<3; i++)
|
||||||
|
field_tot[i] = 0.0;
|
||||||
|
|
||||||
|
fout << endl << "% --------- begin output data -----------------";
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " ALL THE COMPONENTS ARE GIVEN AS COMPONENTS OF THE CRYSTAL, i.e.";
|
||||||
|
fout << endl << " vec_B = B_a vec_a + B_b vec_b + B_c vec_c, etc., where e.g. vec_a is the";
|
||||||
|
fout << endl << " normalized vector along the a-axis of the crystal etc.";
|
||||||
|
fout << endl << " This needs to be kept in mind, especially for triklin, monoklin, ..., structures";
|
||||||
|
fout << endl;
|
||||||
|
QValueVector<dipole_type_info>::iterator it;
|
||||||
|
int i=1;
|
||||||
|
for (it = info->dipole.begin(); it != info->dipole.end(); ++it) {
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " THE DIPOLAR MATRIX FOR THE TYPE " << i++ << " in (kG/mu_B)";
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " D | 1 2 3";
|
||||||
|
fout << endl << " --|-----------------------------------";
|
||||||
|
for (int j=0; j<3; j++) {
|
||||||
|
sprintf(str, " %d | %+2.4f %+2.4f %+2.4f ",
|
||||||
|
j+1, UNITS*it->dipoleMatrix[j][0], UNITS*it->dipoleMatrix[j][1], UNITS*it->dipoleMatrix[j][2]);
|
||||||
|
fout << endl << str;
|
||||||
|
field_tot[j] += UNITS*it->field[j];
|
||||||
|
}
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " trace = " << UNITS*(it->dipoleMatrix[0][0]+it->dipoleMatrix[1][1]+it->dipoleMatrix[2][2]);
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " the dipolar magnetic field generated by the type " << it->type << " ions is: ";
|
||||||
|
fout << UNITS*sqrt(it->field[0]*it->field[0]+it->field[1]*it->field[1]+it->field[2]*it->field[2]) << " (kG)";
|
||||||
|
fout << endl << " with the components: " << UNITS*it->field[0] << ", " << UNITS*it->field[1] << ", " << UNITS*it->field[2] << " (kG)";
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " -----------------------";
|
||||||
|
}
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << " THE TOTAL DIPOLAR MAGNETIC FIELD IS ";
|
||||||
|
fout << sqrt(field_tot[0]*field_tot[0]+field_tot[1]*field_tot[1]+field_tot[2]*field_tot[2]) << " (kG)";
|
||||||
|
fout << endl << " with the components: " << field_tot[0] << ", " << field_tot[1] << ", " << field_tot[2] << " (kG)";
|
||||||
|
fout << endl;
|
||||||
|
fout << endl << "% --------- end output data -----------------";
|
||||||
|
|
||||||
|
file.close();
|
||||||
|
|
||||||
|
return DIPOLAR_SUCCESS;
|
||||||
|
}
|
||||||
|
|
||||||
|
//-----------------------------------------------------------------------
|
||||||
|
/**
|
||||||
|
* <p>
|
||||||
|
*/
|
||||||
|
int main(int argc, char *argv[])
|
||||||
|
{
|
||||||
|
int debug = 0;
|
||||||
|
dipole_info *info;
|
||||||
|
int status;
|
||||||
|
int i;
|
||||||
|
|
||||||
|
// check syntax
|
||||||
|
if (argc < 3) {
|
||||||
|
dipolar_ordered_syntax();
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (argc == 4)
|
||||||
|
sscanf(argv[3], "%d", &debug);
|
||||||
|
|
||||||
|
// allocate memory for the info structure
|
||||||
|
info = new dipole_info;
|
||||||
|
|
||||||
|
// init info structure
|
||||||
|
info->noAtoms = 0;
|
||||||
|
info->noMagTypes = 0;
|
||||||
|
info->latticeParam[0] = 0.0;
|
||||||
|
info->latticeParam[1] = 0.0;
|
||||||
|
info->latticeParam[2] = 0.0;
|
||||||
|
info->angle[0] = 0.0;
|
||||||
|
info->angle[1] = 0.0;
|
||||||
|
info->angle[2] = 0.0;
|
||||||
|
info->muPos[0] = 0.0;
|
||||||
|
info->muPos[1] = 0.0;
|
||||||
|
info->muPos[2] = 0.0;
|
||||||
|
info->lorentzSphereRadius = 0.0;
|
||||||
|
|
||||||
|
// read input file and fill info structure
|
||||||
|
status = dipolar_ordered_read_input_file(argv[1], info);
|
||||||
|
if (status != DIPOLAR_SUCCESS) {
|
||||||
|
switch (status) {
|
||||||
|
case DIPOLAR_FILE_NOT_FOUND:
|
||||||
|
cout << endl << "ERROR: couldn't find input file " << argv[1] << endl;
|
||||||
|
break;
|
||||||
|
case DIPOLAR_FILE_OPEN_FAILED:
|
||||||
|
cout << endl << "ERROR: couldn't open input file " << argv[1] << endl;
|
||||||
|
break;
|
||||||
|
case DIPOLAR_WRONG_INPUT_FILE_FORMAT:
|
||||||
|
cout << endl << "ERROR: wrong input file format ..." << endl;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
cout << endl << "unknown error ..." << endl;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
dipolar_ordered_syntax();
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (debug) {
|
||||||
|
cout << endl << "-----------------------";
|
||||||
|
cout << endl << "[debug] no of magnetic atoms in the unit cell = " << info->noAtoms;
|
||||||
|
cout << endl << "[debug] no of magnetic types in the unit cell = " << info->noMagTypes;
|
||||||
|
cout << endl << "[debug] list of the atom position (A) and its type";
|
||||||
|
cout << endl << "-----------------------";
|
||||||
|
i = 1;
|
||||||
|
for (QValueVector<atom_pos>::iterator it = info->atom.begin(); it != info->atom.end(); ++it) {
|
||||||
|
cout << endl << "[debug] " << i++ << ": ";
|
||||||
|
cout << it->pos[0] << ", ";
|
||||||
|
cout << it->pos[1] << ", ";
|
||||||
|
cout << it->pos[2] << ", ";
|
||||||
|
cout << it->type;
|
||||||
|
}
|
||||||
|
cout << endl << "-----------------------";
|
||||||
|
cout << endl << "[debug] lattice parameters (A) : ";
|
||||||
|
cout << info->latticeParam[0] << ", ";
|
||||||
|
cout << info->latticeParam[1] << ", ";
|
||||||
|
cout << info->latticeParam[2];
|
||||||
|
cout << endl << "[debug] angles (°): ";
|
||||||
|
cout << info->angle[0] << ", ";
|
||||||
|
cout << info->angle[1] << ", ";
|
||||||
|
cout << info->angle[2];
|
||||||
|
cout << endl << "[debug] Lorentz sphere radius (A): ";
|
||||||
|
cout << info->lorentzSphereRadius;
|
||||||
|
cout << endl << "[debug] magnetic moments of the various types";
|
||||||
|
cout << endl << "-----------------------";
|
||||||
|
i = 1;
|
||||||
|
for (QValueVector<mag_moment>::iterator it = info->magMoment.begin(); it != info->magMoment.end(); ++it) {
|
||||||
|
cout << endl << "[debug] " << i++ << ": ";
|
||||||
|
cout << it->component[0] << ", ";
|
||||||
|
cout << it->component[1] << ", ";
|
||||||
|
cout << it->component[2] << ", ";
|
||||||
|
}
|
||||||
|
cout << endl << "-----------------------";
|
||||||
|
if (info->calcTag == CALC_POINT) {
|
||||||
|
cout << endl << "[debug] calculate the dipolar field in a single point of the lattice.";
|
||||||
|
cout << endl << "[debug] muon position in lattice coordinates: ";
|
||||||
|
cout << info->muPos[0] << ", ";
|
||||||
|
cout << info->muPos[1] << ", ";
|
||||||
|
cout << info->muPos[2];
|
||||||
|
} else if (info->calcTag == CALC_PLANE) {
|
||||||
|
cout << endl << "[debug] calculate the dipolar field in a plane of the lattice.";
|
||||||
|
cout << endl << "[debug] vector to the plane: ";
|
||||||
|
cout << info->planeFix[0] << ", ";
|
||||||
|
cout << info->planeFix[1] << ", ";
|
||||||
|
cout << info->planeFix[2];
|
||||||
|
cout << endl << "[debug] plane normal vector: ";
|
||||||
|
cout << info->planeNormal[0] << ", ";
|
||||||
|
cout << info->planeNormal[1] << ", ";
|
||||||
|
cout << info->planeNormal[2];
|
||||||
|
} else if (info->calcTag == CALC_SPHERE) {
|
||||||
|
cout << endl << "[debug] calculate the dipolar field on a sphere.";
|
||||||
|
}
|
||||||
|
cout << endl << "-----------------------";
|
||||||
|
cout << endl << endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
// do the calculus
|
||||||
|
switch (info->calcTag) {
|
||||||
|
case CALC_POINT:
|
||||||
|
dipolar_ordered_point_calc(info);
|
||||||
|
break;
|
||||||
|
case CALC_PLANE:
|
||||||
|
dipolar_ordered_plane_calc(info);
|
||||||
|
break;
|
||||||
|
case CALC_SPHERE:
|
||||||
|
dipolar_ordered_sphere_calc(info);
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
cout << endl << "SEVERE ERROR: you never should get to here!?!?!?! Sorry ;-(" << endl;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
// write the output file
|
||||||
|
status = dipolar_ordered_write_output_file(argv[1], argv[2], info);
|
||||||
|
if (status != DIPOLAR_SUCCESS) {
|
||||||
|
cout << endl << "ERROR while trying to write to the output file " << argv[2];
|
||||||
|
cout << endl << endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
// clean up
|
||||||
|
if (info) {
|
||||||
|
info->atom.clear();
|
||||||
|
delete info;
|
||||||
|
info = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
return 1;
|
||||||
|
}
|
||||||
Loading…
x
Reference in New Issue
Block a user