jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Implementation of SeqNei V1 (Real) 

Still both Seq and SeqNei versions
In SeqNei, loop on Neighb not working yet
This commit is contained in:
julient31 2017-06-19 10:43:54 -06:00
parent bf5b3f96e9
commit b934621651
12 changed files with 3216 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

951
src/atom_vec_spin.cpp Normal file
View File

@ -0,0 +1,951 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "atom_vec_spin.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = 0;
mass_type = 1; //check why
//comm_x_only = 0;
comm_x_only = 1;
//comm_f_only = 1;
comm_f_only = 0;
size_forward = 7;
size_reverse = 6;
size_border = 11;
size_velocity = 3;
size_data_atom = 9; //to check later
size_data_vel = 4;
xcol_data = 4;
forceclearflag = 1;
atom->mumag_flag = atom->sp_flag = 1;
}
/* ----------------------------------------------------------------------
grow atom arrays
n = 0 grows arrays by a chunk
n > 0 allocates arrays to size n
------------------------------------------------------------------------- */
void AtomVecSpin::grow(int n)
{
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");
type = memory->grow(atom->type,nmax,"atom:type");
mask = memory->grow(atom->mask,nmax,"atom:mask");
image = memory->grow(atom->image,nmax,"atom:image");
x = memory->grow(atom->x,nmax,3,"atom:x");
v = memory->grow(atom->v,nmax,3,"atom:v");
f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
//Allocating mag. quantities
mumag = memory->grow(atom->mumag,nmax,"atom:mumag");
sp = memory->grow(atom->sp,nmax,4,"atom:sp");
fm = memory->grow(atom->fm,nmax*comm->nthreads,3,"atom:fm");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
}
/* ----------------------------------------------------------------------
reset local array ptrs
------------------------------------------------------------------------- */
void AtomVecSpin::grow_reset()
{
tag = atom->tag; type = atom->type;
mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f;
mumag = atom->mumag; sp = atom->sp; fm = atom->fm;
}
/* ----------------------------------------------------------------------
copy atom I info to atom J
------------------------------------------------------------------------- */
void AtomVecSpin::copy(int i, int j, int delflag)
{
tag[j] = tag[i];
type[j] = type[i];
mask[j] = mask[i];
image[j] = image[i];
x[j][0] = x[i][0];
x[j][1] = x[i][1];
x[j][2] = x[i][2];
v[j][0] = v[i][0];
v[j][1] = v[i][1];
v[j][2] = v[i][2];
mumag[j] = mumag[i];
sp[j][0] = sp[i][0];
sp[j][1] = sp[i][1];
sp[j][2] = sp[i][2];
sp[j][3] = sp[i][3];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag);
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
}
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_comm_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
dz = pbc[2]*domain->zprd;
}
if (!deform_vremap) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
}
}
}
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_comm_hybrid(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecSpin::unpack_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
void AtomVecSpin::unpack_comm_vel(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::unpack_comm_hybrid(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
}
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_reverse(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = f[i][0];
buf[m++] = f[i][1];
buf[m++] = f[i][2];
buf[m++] = fm[i][0];
buf[m++] = fm[i][1];
buf[m++] = fm[i][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecSpin::unpack_reverse(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
f[j][0] += buf[m++];
f[j][1] += buf[m++];
f[j][2] += buf[m++];
fm[j][0] += buf[m++];
fm[j][1] += buf[m++];
fm[j][2] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_border(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = mumag[j];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0];
dy = pbc[1];
dz = pbc[2];
}
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = mumag[j];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_border_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
m = 0;
if (pbc_flag == 0) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0];
buf[m++] = x[j][1];
buf[m++] = x[j][2];
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = mumag[j];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
if (domain->triclinic == 0) {
dx = pbc[0]*domain->xprd;
dy = pbc[1]*domain->yprd;
dz = pbc[2]*domain->zprd;
} else {
dx = pbc[0];
dy = pbc[1];
dz = pbc[2];
}
if (!deform_vremap) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = mumag[j];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = mumag[j];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
}
}
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::pack_border_hybrid(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = mumag[j];
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
return m;
}
/* ---------------------------------------------------------------------- */
void AtomVecSpin::unpack_border(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (i == nmax) grow(0);
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
tag[i] = (tagint) ubuf(buf[m++]).i;
type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i;
mumag[i] = buf[m++];
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->
unpack_border(n,first,&buf[m]);
}
/* ---------------------------------------------------------------------- */
void AtomVecSpin::unpack_border_vel(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
if (i == nmax) grow(0);
x[i][0] = buf[m++];
x[i][1] = buf[m++];
x[i][2] = buf[m++];
tag[i] = (tagint) ubuf(buf[m++]).i;
type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i;
mumag[i] = buf[m++];
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
}
if (atom->nextra_border)
for (int iextra = 0; iextra < atom->nextra_border; iextra++)
m += modify->fix[atom->extra_border[iextra]]->
unpack_border(n,first,&buf[m]);
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::unpack_border_hybrid(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
mumag[i] = buf[m++];
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
}
return m;
}
/* ----------------------------------------------------------------------
pack all atom quantities for shipping to another proc
xyz must be 1st 3 values, so that comm::exchange can test on them
------------------------------------------------------------------------- */
int AtomVecSpin::pack_exchange(int i, double *buf)
{
int m = 1;
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
buf[m++] = v[i][0];
buf[m++] = v[i][1];
buf[m++] = v[i][2];
buf[m++] = ubuf(tag[i]).d;
buf[m++] = ubuf(type[i]).d;
buf[m++] = ubuf(mask[i]).d;
buf[m++] = ubuf(image[i]).d;
buf[m++] = mumag[i];
buf[m++] = sp[i][0];
buf[m++] = sp[i][1];
buf[m++] = sp[i][2];
buf[m++] = sp[i][3];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
buf[0] = m;
return m;
}
/* ---------------------------------------------------------------------- */
int AtomVecSpin::unpack_exchange(double *buf)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
int m = 1;
x[nlocal][0] = buf[m++];
x[nlocal][1] = buf[m++];
x[nlocal][2] = buf[m++];
v[nlocal][0] = buf[m++];
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
tag[nlocal] = (tagint) ubuf(buf[m++]).i;
type[nlocal] = (int) ubuf(buf[m++]).i;
mask[nlocal] = (int) ubuf(buf[m++]).i;
image[nlocal] = (imageint) ubuf(buf[m++]).i;
mumag[nlocal] = buf[m++];
sp[nlocal][0] = buf[m++];
sp[nlocal][1] = buf[m++];
sp[nlocal][2] = buf[m++];
sp[nlocal][3] = buf[m++];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->
unpack_exchange(nlocal,&buf[m]);
atom->nlocal++;
return m;
}
/* ----------------------------------------------------------------------
size of restart data for all atoms owned by this proc
include extra data stored by fixes
------------------------------------------------------------------------- */
int AtomVecSpin::size_restart()
{
int i;
int nlocal = atom->nlocal;
int n = 16 * nlocal;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
for (i = 0; i < nlocal; i++)
n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
return n;
}
/* ----------------------------------------------------------------------
pack atom I's data for restart file including extra quantities
xyz must be 1st 3 values, so that read_restart can test on them
molecular types may be negative, but write as positive
------------------------------------------------------------------------- */
int AtomVecSpin::pack_restart(int i, double *buf)
{
int m = 1;
buf[m++] = x[i][0];
buf[m++] = x[i][1];
buf[m++] = x[i][2];
buf[m++] = ubuf(tag[i]).d;
buf[m++] = ubuf(type[i]).d;
buf[m++] = ubuf(mask[i]).d;
buf[m++] = ubuf(image[i]).d;
buf[m++] = v[i][0];
buf[m++] = v[i][1];
buf[m++] = v[i][2];
buf[m++] = mumag[i];
buf[m++] = sp[i][0];
buf[m++] = sp[i][1];
buf[m++] = sp[i][2];
buf[m++] = sp[i][3];
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
buf[0] = m;
return m;
}
/* ----------------------------------------------------------------------
unpack data for one atom from restart file including extra quantities
------------------------------------------------------------------------- */
int AtomVecSpin::unpack_restart(double *buf)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) {
grow(0);
if (atom->nextra_store)
memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
}
int m = 1;
x[nlocal][0] = buf[m++];
x[nlocal][1] = buf[m++];
x[nlocal][2] = buf[m++];
tag[nlocal] = (tagint) ubuf(buf[m++]).i;
type[nlocal] = (int) ubuf(buf[m++]).i;
mask[nlocal] = (int) ubuf(buf[m++]).i;
image[nlocal] = (imageint) ubuf(buf[m++]).i;
v[nlocal][0] = buf[m++];
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
mumag[nlocal] = buf[m++];
sp[nlocal][0] = buf[m++];
sp[nlocal][1] = buf[m++];
sp[nlocal][2] = buf[m++];
sp[nlocal][3] = buf[m++];
double **extra = atom->extra;
if (atom->nextra_store) {
int size = static_cast<int> (buf[0]) - m;
for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
}
atom->nlocal++;
return m;
}
/* ----------------------------------------------------------------------
create one atom of itype at coord
set other values to defaults
------------------------------------------------------------------------- */
void AtomVecSpin::create_atom(int itype, double *coord)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = 0;
type[nlocal] = itype;
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
mumag[nlocal] = 0.0;
sp[nlocal][0] = 0.0;
sp[nlocal][1] = 0.0;
sp[nlocal][2] = 0.0;
sp[nlocal][3] = 0.0;
atom->nlocal++;
}
/* ----------------------------------------------------------------------
unpack one line from Atoms section of data file
initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecSpin::data_atom(double *coord, imageint imagetmp, char **values)
{
int nlocal = atom->nlocal;
if (nlocal == nmax) grow(0);
tag[nlocal] = ATOTAGINT(values[0]);
type[nlocal] = atoi(values[1]);
if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
error->one(FLERR,"Invalid atom type in Atoms section of data file");
mumag[nlocal] = atof(values[2]);
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
sp[nlocal][0] = atof(values[6]);
sp[nlocal][1] = atof(values[7]);
sp[nlocal][2] = atof(values[8]);
sp[nlocal][3] = sqrt(sp[nlocal][0]*sp[nlocal][0] +
sp[nlocal][1]*sp[nlocal][1] +
sp[nlocal][2]*sp[nlocal][2]);
image[nlocal] = imagetmp;
mask[nlocal] = 1;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
atom->nlocal++;
}
/* ----------------------------------------------------------------------
unpack hybrid quantities from one line in Atoms section of data file
initialize other atom quantities for this sub-style
------------------------------------------------------------------------- */
int AtomVecSpin::data_atom_hybrid(int nlocal, char **values)
{
mumag[nlocal] = atof(values[0]);
sp[nlocal][0] = atof(values[1]);
sp[nlocal][1] = atof(values[2]);
sp[nlocal][2] = atof(values[3]);
sp[nlocal][3] = sqrt(sp[nlocal][0]*sp[nlocal][0] +
sp[nlocal][1]*sp[nlocal][1] +
sp[nlocal][2]*sp[nlocal][2]);
return 4;
}
/* ----------------------------------------------------------------------
pack atom info for data file including 3 image flags
------------------------------------------------------------------------- */
void AtomVecSpin::pack_data(double **buf)
{
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
buf[i][0] = ubuf(tag[i]).d;
buf[i][1] = ubuf(type[i]).d;
buf[i][2] = mumag[i];
buf[i][3] = x[i][0];
buf[i][4] = x[i][1];
buf[i][5] = x[i][2];
buf[i][6] = sp[i][0];
buf[i][7] = sp[i][1];
buf[i][8] = sp[i][2];
buf[i][9] = sp[i][3];
buf[i][10] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
buf[i][11] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
buf[i][12] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
}
}
/* ----------------------------------------------------------------------
pack hybrid atom info for data file
------------------------------------------------------------------------- */
int AtomVecSpin::pack_data_hybrid(int i, double *buf)
{
buf[0] = mumag[i];
buf[1] = sp[i][0];
buf[2] = sp[i][1];
buf[3] = sp[i][2];
buf[4] = sp[i][3];
return 5;
}
/* ----------------------------------------------------------------------
write atom info to data file including 3 image flags
------------------------------------------------------------------------- */
void AtomVecSpin::write_data(FILE *fp, int n, double **buf)
{
for (int i = 0; i < n; i++)
fprintf(fp,TAGINT_FORMAT \
" %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e "
"%-1.16e %d %d %d\n",
(tagint) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
buf[i][2],buf[i][3],buf[i][4],
buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],
(int) ubuf(buf[i][10]).i,(int) ubuf(buf[i][11]).i,
(int) ubuf(buf[i][12]).i);
}
/* ----------------------------------------------------------------------
write hybrid atom info to data file
------------------------------------------------------------------------- */
int AtomVecSpin::write_data_hybrid(FILE *fp, double *buf)
{
fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2],buf[3],buf[4]);
return 4;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
bigint AtomVecSpin::memory_usage()
{
bigint bytes = 0;
if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
if (atom->memcheck("fm")) bytes += memory->usage(fm,nmax*comm->nthreads,3);
if (atom->memcheck("mumag")) bytes += memory->usage(mumag,nmax);
if (atom->memcheck("sp")) bytes += memory->usage(sp,nmax,4);
return bytes;
}
void AtomVecSpin::force_clear(int n, size_t nbytes)
{
memset(&atom->f[0][0],0,3*nbytes);
memset(&atom->fm[0][0],0,3*nbytes);
}

90
src/atom_vec_spin.h Normal file
View File

@ -0,0 +1,90 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef ATOM_CLASS
AtomStyle(spin,AtomVecSpin)
#else
#ifndef LMP_ATOM_VEC_SPIN_H
#define LMP_ATOM_VEC_SPIN_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecSpin : public AtomVec {
public:
AtomVecSpin(class LAMMPS *);
void grow(int);
void grow_reset();
void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
void unpack_comm(int, int, double *);
void unpack_comm_vel(int, int, double *);
int unpack_comm_hybrid(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
int pack_border(int, int *, double *, int, int *);
int pack_border_vel(int, int *, double *, int, int *);
int pack_border_hybrid(int, int *, double *);
void unpack_border(int, int, double *);
void unpack_border_vel(int, int, double *);
int unpack_border_hybrid(int, int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, imageint, char **);
int data_atom_hybrid(int, char **);
void pack_data(double **);
int pack_data_hybrid(int, double *);
void write_data(FILE *, int, double **);
int write_data_hybrid(FILE *, double *);
bigint memory_usage();
//Test force clear
void force_clear(int, size_t);
private:
tagint *tag;
int *type,*mask;
imageint *image;
double **x,**v,**f; //MD quantities: position, velocity and force
double *mumag,**sp, **fm; //Magnetic quantities: mu, spin direction, magnetic force
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Per-processor system is too big
The number of owned atoms plus ghost atoms on a single
processor must fit in 32-bit integer.
E: Invalid atom type in Atoms section of data file
Atom types must range from 1 to specified # of types.
*/

145
src/compute_spin.cpp Normal file
View File

@ -0,0 +1,145 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <string.h>
#include "compute_spin.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "math_special.h"
#include "math_const.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace MathSpecial;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), mag(NULL)
{
if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command");
vector_flag = 1;
size_vector = 7;
extvector = 0;
init();
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeSpin::~ComputeSpin()
{
memory->destroy(mag);
}
/* ---------------------------------------------------------------------- */
void ComputeSpin::init()
{
hbar = force->hplanck/MY_2PI;
kb = force->boltz;
}
/* ---------------------------------------------------------------------- */
void ComputeSpin::compute_vector()
{
int i, index;
invoked_vector = update->ntimestep;
countsp = countsptot = 0.0;
mag[0] = mag[1] = mag[2] = mag[3] = 0.0;
magtot[0] = magtot[1] = magtot[2] = magtot[3] = 0.0;
magenergy = magenergytot = 0.0;
tempnum = tempnumtot = 0.0;
tempdenom = tempdenomtot = 0.0;
spintemperature = 0.0;
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mumag = atom->mumag;
double **sp = atom->sp;
double **fm = atom->fm;
double tx,ty,tz;
int nlocal = atom->nlocal;
// compute total magnetization and magnetic energy
// compute spin temperature; See Nurdin et al., Phys. Rev. E 61, 2000
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (atom->mumag_flag && atom->sp_flag) {
mag[0] += sp[i][0];
mag[1] += sp[i][1];
mag[2] += sp[i][2];
magenergy += mumag[i]*sp[i][0]*fm[i][0];
magenergy += mumag[i]*sp[i][1]*fm[i][1];
magenergy += mumag[i]*sp[i][2]*fm[i][2];
tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0];
tempnum += tx*tx+ty*ty+tz*tz;
tempdenom += sp[i][0]*sp[i][0]+sp[i][1]*sp[i][1]+sp[i][2]*sp[i][2];
countsp++;
}
}
else error->all(FLERR,"Compute spin/compute declared magnetic quantities (sp and mumag flags)");
}
MPI_Allreduce(mag,magtot,4,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&magenergy,&magenergytot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&tempnum,&tempnumtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&tempdenom,&tempdenomtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&countsp,&countsptot,1,MPI_DOUBLE,MPI_SUM,world);
double scale = 1.0/countsptot;
magtot[0] *= scale;
magtot[1] *= scale;
magtot[2] *= scale;
magtot[3] = sqrt(square(magtot[0])+square(magtot[1])+square(magtot[2]));
spintemperature = hbar*tempnumtot/2.0/kb/tempdenomtot;
vector[0] = invoked_vector*update->dt;
vector[1] = magtot[0];
vector[2] = magtot[1];
vector[3] = magtot[2];
vector[4] = magtot[3];
vector[5] = -0.5*magenergytot*hbar;
vector[6] = spintemperature;
}
/* ----------------------------------------------------------------------
free and reallocate arrays
------------------------------------------------------------------------- */
void ComputeSpin::allocate()
{
memory->destroy(mag);
memory->create(mag,4,"compute/spin:mag");
memory->create(magtot,5,"compute/spin:mag");
memory->create(vector,7,"compute/spin:vector");
}

71
src/compute_spin.h Normal file
View File

@ -0,0 +1,71 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(compute/spin,ComputeSpin)
#else
#ifndef LMP_COMPUTE_SPIN_H
#define LMP_COMPUTE_SPIN_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSpin : public Compute {
public:
ComputeSpin(class LAMMPS *, int, char **);
~ComputeSpin();
void init();
void compute_vector();
private:
double *mag;
double *magtot;
double magenergy;
double magenergytot;
double tempnum,tempnumtot;
double tempdenom,tempdenomtot;
double spintemperature;
double kb,hbar;
int countsp;
int countsptot;
int usecenter;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Chunk/atom compute does not exist for compute compute/spin
Self-explanatory.
E: Compute compute/spin does not use chunk/atom compute
The style of the specified compute is not chunk/atom.
*/

255
src/fix_force_spin.cpp Normal file
View File

@ -0,0 +1,255 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "fix_force_spin.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "respa.h"
#include "modify.h"
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "error.h"
#include "force.h"
#include "memory.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{ZEEMAN,ANISOTROPY};
enum{CONSTANT,EQUAL};
/* ---------------------------------------------------------------------- */
FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal fix spin command");
// 7 arguments for a force/spin fix command:
//(fix ID group force/spin magnitude (T or eV) style (zeeman or anisotropy) direction (3 cartesian coordinates)
//Magnetic interactions only coded for cartesian coordinates
dynamic_group_allow = 1;
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
respa_level_support = 1;
ilevel_respa = 0;
magstr = NULL;
magfieldstyle = CONSTANT;
H_field = 0.0;
Hx = Hy = Hz = 0.0;
Ka = 0.0;
Kax = Kay = Kaz = 0.0;
zeeman_flag = aniso_flag = 0;
if (strcmp(arg[3],"zeeman") == 0) {
if (narg != 8) error->all(FLERR,"Illegal fix zeeman command");
style = ZEEMAN;
zeeman_flag = 1;
H_field = force->numeric(FLERR,arg[4]);
Hx = force->numeric(FLERR,arg[5]);
Hy = force->numeric(FLERR,arg[6]);
Hz = force->numeric(FLERR,arg[7]);
magfieldstyle = CONSTANT;
} else if (strcmp(arg[3],"anisotropy") == 0) {
if (narg != 8) error->all(FLERR,"Illegal fix anisotropy command");
style = ANISOTROPY;
aniso_flag = 1;
Ka = force->numeric(FLERR,arg[4]);
Kax = force->numeric(FLERR,arg[5]);
Kay = force->numeric(FLERR,arg[6]);
Kaz = force->numeric(FLERR,arg[7]);
} else error->all(FLERR,"Illegal fix force/spin command");
degree2rad = MY_PI/180.0;
time_origin = update->ntimestep;
eflag = 0;
emag = 0.0;
}
/* ---------------------------------------------------------------------- */
FixForceSpin::~FixForceSpin()
{
memory->destroy(spi);
memory->destroy(fmi);
delete [] magstr;
}
/* ---------------------------------------------------------------------- */
int FixForceSpin::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= THERMO_ENERGY;
mask |= POST_FORCE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixForceSpin::init()
{
double hbar = force->hplanck/MY_2PI; //eV/(rad.THz)
double mub = 5.78901e-5; //in eV/T
double gyro = mub/hbar; //in rad.THz/T
H_field *= gyro; //in rad.THz
Ka /= hbar; //in rad.THz
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
}
// check variables
if (magstr) {
magvar = input->variable->find(magstr);
if (magvar < 0)
error->all(FLERR,"Variable name for fix magnetic field does not exist");
if (!input->variable->equalstyle(magvar))
error->all(FLERR,"Variable for fix magnetic field is invalid style");
}
varflag = CONSTANT;
if (magfieldstyle != CONSTANT) varflag = EQUAL;
// set magnetic field components once and for all
if (varflag == CONSTANT) set_magneticforce();
memory->create(spi,3,"forcespin:spi");
memory->create(fmi,3,"forcespin:fmi");
}
/* ---------------------------------------------------------------------- */
void FixForceSpin::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(ilevel_respa);
post_force_respa(vflag,ilevel_respa,0);
((Respa *) update->integrate)->copy_f_flevel(ilevel_respa);
}
}
/* ---------------------------------------------------------------------- */
void FixForceSpin::post_force(int vflag)
{
// update gravity due to variables
if (varflag != CONSTANT) {
modify->clearstep_compute();
modify->addstep_compute(update->ntimestep + 1);
set_magneticforce(); //Update value of the mag. field if time-dependent
}
// double **x = atom->x;
double **sp = atom->sp;
double *mumag = atom->mumag;
double **fm = atom->fm;
int nlocal = atom->nlocal;
double scalar;
eflag = 0;
emag = 0.0;
for (int i = 0; i < nlocal; i++) {
fmi[0] = fmi[1] = fmi[2] = 0.0;
//if (style == ZEEMAN) {
if (zeeman_flag) {
compute_zeeman(i,fmi);
}
//if (style == ANISOTROPY) {
if (aniso_flag) {
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
compute_anisotropy(i,spi,fmi);
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
}
}
/* ---------------------------------------------------------------------- */
void FixForceSpin::compute_zeeman(int i, double *fmi)
{
double *mumag = atom->mumag;
fmi[0] += mumag[i]*xmag;
fmi[1] += mumag[i]*ymag;
fmi[2] += mumag[i]*zmag;
}
void FixForceSpin::compute_anisotropy(int i, double * spi, double *fmi)
{
double scalar = Kax*spi[0] + Kay*spi[1] + Kaz*spi[2];
fmi[0] += scalar*xmag;
fmi[1] += scalar*ymag;
fmi[2] += scalar*zmag;
}
/* ---------------------------------------------------------------------- */
void FixForceSpin::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == ilevel_respa) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
//No acceleration for magnetic EOM, only a "magnetic force"
void FixForceSpin::set_magneticforce()
{
if (style == ZEEMAN) {
xmag = H_field*Hx;
ymag = H_field*Hy;
zmag = H_field*Hz;
}
if (style == ANISOTROPY) {
xmag = 2.0*Ka*Kax;
ymag = 2.0*Ka*Kay;
zmag = 2.0*Ka*Kaz;
}
}
/* ----------------------------------------------------------------------
potential energy in magnetic field
------------------------------------------------------------------------- */
double FixForceSpin::compute_scalar()
{
// only sum across procs one time
if (eflag == 0) {
MPI_Allreduce(&emag,&emag_all,1,MPI_DOUBLE,MPI_SUM,world);
eflag = 1;
}
return emag_all;
}

92
src/fix_force_spin.h Normal file
View File

@ -0,0 +1,92 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(force/spin,FixForceSpin)
#else
#ifndef LMP_FIX_FORCE_SPIN_H
#define LMP_FIX_FORCE_SPIN_H
#include "fix.h"
namespace LAMMPS_NS {
class FixForceSpin : public Fix {
friend class FixPour;
public:
FixForceSpin(class LAMMPS *, int, char **);
~FixForceSpin();
int setmask();
void init();
void setup(int);
virtual void post_force(int);
virtual void post_force_respa(int, int, int);
double compute_scalar();
int zeeman_flag, aniso_flag;
void compute_zeeman(int, double *);
void compute_anisotropy(int, double *, double *);
protected:
int style;
double xmag, ymag, zmag; //Magnetic force
double degree2rad;
int ilevel_respa;
int time_origin;
int eflag;
double emag, emag_all;
int varflag;
int magfieldstyle;
int magvar;
char *magstr;
double H_field; //Zeeman field intensity and direction
double Hx, Hy, Hz;
double Ka; //Magnetic anisotropy intensity and direction
double Kax, Kay, Kaz;
double *spi, *fmi; //Temp var. in compute
void set_magneticforce();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Variable name for fix force/spin does not exist
Self-explanatory.
E: Variable for fix force/spin is invalid style
Only equal-style variables can be used.
*/

200
src/fix_langevin_spin.cpp Normal file
View File

@ -0,0 +1,200 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Carolyn Phillips (U Mich), reservoir energy tally
Aidan Thompson (SNL) GJF formulation
------------------------------------------------------------------------- */
#include <mpi.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "fix_langevin_spin.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "comm.h"
#include "input.h"
#include "variable.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "group.h"
#include "math_const.h"
#include "random_park.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), id_temp(NULL), random(NULL)
{
if (narg != 7) error->all(FLERR,"Illegal fix langevin/spin command");
dynamic_group_allow = 1;
scalar_flag = 1;
global_freq = 1;
extscalar = 1;
nevery = 1;
temp = force->numeric(FLERR,arg[3]);
alpha_t = force->numeric(FLERR,arg[4]);
alpha_l = force->numeric(FLERR,arg[5]);
seed = force->inumeric(FLERR,arg[6]);
if (alpha_t < 0.0) error->all(FLERR,"Fix langevin/spin transverse damping must be >= 0.0");
if (alpha_l < 0.0) error->all(FLERR,"Fix langevin/spin transverse damping must be >= 0.0");
if (seed <= 0) error->all(FLERR,"Illegal fix langevin/spin seed must be > 0");
// initialize Marsaglia RNG with processor-unique seed
//random = new RanMars(lmp,seed + comm->me);
random = new RanPark(lmp,seed + comm->me);
}
/* ---------------------------------------------------------------------- */
FixLangevinSpin::~FixLangevinSpin()
{
delete random;
}
/* ---------------------------------------------------------------------- */
int FixLangevinSpin::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
mask |= THERMO_ENERGY;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixLangevinSpin::init()
{
// warn if any fix comes after this one
int after = 0;
int flag_force = 0;
int flag_lang = 0;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp("force/spin",modify->fix[i]->style)==0) flag_force = MAX(flag_force,i);
if (strcmp("langevin/spin",modify->fix[i]->style)==0) flag_lang = i;
}
if (flag_force >= flag_lang) error->all(FLERR,"Fix langevin/spin should come after all other spin fixes");
dts = update->dt;
Gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t));
double hbar = force->hplanck/MY_2PI; //eV/(rad.THz)
double kb = force->boltz;
D = (MY_2PI*Gil_factor*kb*temp)/hbar/dts;
sigma = sqrt(D);
}
/* ---------------------------------------------------------------------- */
void FixLangevinSpin::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixLangevinSpin::post_force(int vflag)
{
double **sp = atom->sp;
double **fm = atom->fm;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double sx, sy, sz;
double fmx, fmy, fmz;
double cpx, cpy, cpz;
double rx, ry, rz;
// add the damping to the effective field of each spin
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
sx = sp[i][0];
sy = sp[i][1];
sz = sp[i][2];
fmx = fm[i][0];
fmy = fm[i][1];
fmz = fm[i][2];
cpx = fmy*sz - fmz*sy;//Computing cross product
cpy = fmz*sx - fmx*sz;
cpz = fmx*sy - fmy*sx;
fmx -= alpha_t*cpx;//Taking the damping value away
fmy -= alpha_t*cpy;
fmz -= alpha_t*cpz;
fm[i][0] = fmx;
fm[i][1] = fmy;
fm[i][2] = fmz;
}
//apply thermal effects adding random fields to fm
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
#define GAUSSIAN_R
#if defined GAUSSIAN_R
rx = sigma*random->gaussian();//Drawing random distributions
ry = sigma*random->gaussian();
rz = sigma*random->gaussian();
#else
rx = sigma*(random->uniform() - 0.5);
ry = sigma*(random->uniform() - 0.5);
rz = sigma*(random->uniform() - 0.5);
#endif
fm[i][0] += rx;//Adding random field
fm[i][1] += ry;
fm[i][2] += rz;
fm[i][0] *= Gil_factor;//Multiplying by Gilbert's prefactor
fm[i][1] *= Gil_factor;
fm[i][2] *= Gil_factor;
}
}
/* ---------------------------------------------------------------------- */
void FixLangevinSpin::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}

106
src/fix_langevin_spin.h Normal file
View File

@ -0,0 +1,106 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(langevin/spin,FixLangevinSpin)
#else
#ifndef LMP_FIX_LANGEVIN_SPIN_H
#define LMP_FIX_LANGEVIN_SPIN_H
#include "fix.h"
namespace LAMMPS_NS {
class FixLangevinSpin : public Fix {
public:
FixLangevinSpin(class LAMMPS *, int, char **);
virtual ~FixLangevinSpin();
int setmask();
void init();
void setup(int);
virtual void post_force(int);
void post_force_respa(int, int, int);
protected:
//First mag. quantities
int transv_damp_flag, long_damp_flag; //Flags for transverse or longitudinal mag. dampings
double alpha_t, alpha_l; //Transverse and long. damping value
double dts,temp,D,sigma;//timestep, temp, noise
double Gil_factor;//Gilbert's prefactor
char *id_temp;
class Compute *temperature;
int nlevels_respa;
//class RanMars *random;
class RanPark *random;
int seed;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix langevin period must be > 0.0
The time window for temperature relaxation must be > 0
W: Energy tally does not account for 'zero yes'
The energy removed by using the 'zero yes' flag is not accounted
for in the energy tally and thus energy conservation cannot be
monitored in this case.
E: Variable for fix langevin is invalid style
It must be an equal-style variable.
E: Cannot zero Langevin force of 0 atoms
The group has zero atoms, so you cannot request its force
be zeroed.
E: Fix langevin variable returned negative temperature
Self-explanatory.
E: Could not find fix_modify temperature ID
The compute ID for computing temperature does not exist.
E: Fix_modify temperature ID does not compute temperature
The compute ID assigned to the fix must compute temperature.
W: Group for fix_modify temp != fix group
The fix_modify command is specifying a temperature computation that
computes a temperature on a different group of atoms than the fix
itself operates on. This is probably not what you want to do.
*/

562
src/fix_nve_spin.cpp Normal file
View File

@ -0,0 +1,562 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "fix_nve_spin.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
#include "respa.h"
#include "force.h"
#include "error.h"
#include "math_vector.h"
#include "math_extra.h"
#include "math_const.h"
#include "modify.h"
//Add headers (see delete later)
#include "pair.h"
#include "timer.h"
#include "integrate.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "pair_spin.h"
#include "memory.h"
#include "fix_force_spin.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
using namespace MathExtra;
enum{NONE,SPIN};
/* ---------------------------------------------------------------------- */
FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
FixNVE(lmp, narg, arg)
{
if (narg < 3) error->all(FLERR,"Illegal fix nve/spin command");
time_integrate = 1;
extra = NONE;
int iarg = 2;
if (strcmp(arg[iarg],"nve/spin") == 0) {
if (iarg+1 > narg) error->all(FLERR,"Illegal fix nve/spin command");
extra = SPIN;
}
// error checks
if (extra == SPIN && !atom->mumag_flag)
error->all(FLERR,"Fix nve/spin requires spin attribute mumag");
#if defined SEQNEI
lockpairspin = NULL;
lockforcespin = NULL;
exch_flag = dmi_flag = me_flag = 0;
zeeman_flag = aniso_flag = 0;
#endif
}
/* ---------------------------------------------------------------------- */
FixNVESpin::~FixNVESpin(){
#if defined SEQNEI
//delete lockpairspin;
//delete lockforcespin;
memory->destroy(spi);
memory->destroy(fmi);
memory->destroy(fmj);
#endif
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::init()
{
FixNVE::init();
dts = update->dt;
#if defined SEQNEI
lockpairspin = (PairSpin *) force->pair;
memory->create(spi,3,"nves:spi");
memory->create(fmi,3,"nves:fmi");
memory->create(fmj,3,"nves:fmj");
int iforce;
for (iforce = 0; iforce < modify->nfix; iforce++)
if (strstr(modify->fix[iforce]->style,"force/spin")) break;
lockforcespin = (FixForceSpin *) modify->fix[iforce];
exch_flag = lockpairspin->exch_flag;
dmi_flag = lockpairspin->dmi_flag;
me_flag = lockpairspin->me_flag;
zeeman_flag = lockforcespin->zeeman_flag;
aniso_flag = lockforcespin->aniso_flag;
#endif
/*int idamp;
for (idamp = 0; idamp < modify->nfix; idamp++)
if (strstr(modify->fix[idamp]->style,"damping/spin")) break;
if (idamp == modify->nfix)
error->all(FLERR,"Integration of spin systems requires use of fix damping (set damping to 0.0 for NVE)");
lockspindamping = (FixSpinDamping *) modify->fix[idamp];
alpha_t = lockspindamping->get_damping(0);
*/
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::initial_integrate(int vflag)
{
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy;
double cp[3],g[3];
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **sp = atom->sp;
double **fm = atom->fm;
double *rmass = atom->rmass;
double *mass = atom->mass;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
int *type = atom->type;
int *mask = atom->mask;
// Advance half spins all particles
//See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011
if (extra == SPIN) {
#if defined SEQNEI
for (int i = 0; i < nlocal; i++){
ComputeSpinInteractionsNei(i);
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#endif
#if defined SEQ
for (int i = 0; i < nlocal; i++){
AdvanceSingleSpin(i,0.5*dts,sp,fm);
ComputeSpinInteractions();
}
#endif
}
// update half v all particles
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (rmass) dtfm = dtf / rmass[i];
else dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
// update x for all particles
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += 0.5 * dtv * v[i][0];
x[i][1] += 0.5 * dtv * v[i][1];
x[i][2] += 0.5 * dtv * v[i][2];
}
}
}
#if defined SEQNEI
/* ---------------------------------------------------------------------- */
void FixNVESpin::ComputeSpinInteractionsNei(int ii)
{
int nflag,sortflag;
int nlocal = atom->nlocal;
int n_post_integrate = modify->n_post_integrate;
int n_pre_exchange = modify->n_pre_exchange;
int n_pre_neighbor = modify->n_pre_neighbor;
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
int n_end_of_step = modify->n_end_of_step;
bigint ntimestep;
ntimestep = update->ntimestep;
//Force compute quantities
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double **x = atom->x;
double **sp = atom->sp;
double **fm = atom->fm;
int *type = atom->type;
int newton_pair = force->newton_pair;
inum = lockpairspin->list->inum;
ilist = lockpairspin->list->ilist;
numneigh = lockpairspin->list->numneigh;
firstneigh = lockpairspin->list->firstneigh;
double xtmp,ytmp,ztmp;
double rsq,rd,delx,dely,delz;
double cut_ex_2, cut_dmi_2, cut_me_2;
cut_ex_2 = cut_dmi_2 = cut_me_2 = 0.0;
int eflag = 1;
int vflag = 0;
int pair_compute_flag = 1;
if (atom->sortfreq > 0) sortflag = 1;
else sortflag = 0;
if (n_post_integrate) modify->post_integrate();
timer->stamp(Timer::MODIFY);
// regular communication vs neighbor list rebuild
nflag = neighbor->decide();
if (nflag == 0) {
timer->stamp();
comm->forward_comm();
timer->stamp(Timer::COMM);
} else {
if (n_pre_exchange) {
timer->stamp();
modify->pre_exchange();
timer->stamp(Timer::MODIFY);
}
//if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
if (sortflag && ntimestep >= atom->nextsort) atom->sort();
comm->borders();
//if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
if (n_pre_neighbor) {
modify->pre_neighbor();
timer->stamp(Timer::MODIFY);
}
neighbor->build();
timer->stamp(Timer::NEIGH);
}
///////Force computation for spin i/////////////
i = ilist[ii];
//Clear atom i
fm[i][0] = fm[i][1] = fm[i][2] = 0.0;
timer->stamp();
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
fmi[0] = fmi[1] = fmi[2] = 0.0;
fmj[0] = fmj[1] = fmj[2] = 0.0;
jlist = firstneigh[i];
jnum = numneigh[i];
// printf("Test inum: %g \n",inum);
/*
//Pair interaction
for (int jj = 0; jj < inum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
itype = type[ii];
jtype = type[j];
if (exch_flag) {
cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]);
if (rsq <= cut_ex_2) {
lockpairspin->compute_exchange(i,j,rsq,fmi,fmj);
}
}
if (dmi_flag) {
cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]);
if (rsq <= cut_dmi_2) {
lockpairspin->compute_dmi(i,j,fmi,fmj);
}
}
if (me_flag) {
cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]);
if (rsq <= cut_me_2) {
lockpairspin->compute_me(i,j,fmi,fmj);
}
}
}
*/
//Pair interaction
int natom = nlocal + atom->nghost;
for (int k = 0; k < natom; k++) {
delx = xtmp - x[k][0];
dely = ytmp - x[k][1];
delz = ztmp - x[k][2];
rsq = delx*delx + dely*dely + delz*delz;
itype = type[ii];
jtype = type[k];
if (exch_flag) {
cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]);
if (rsq <= cut_ex_2) {
lockpairspin->compute_exchange(i,k,rsq,fmi,fmj);
}
}
if (dmi_flag) {
cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]);
if (rsq <= cut_dmi_2) {
lockpairspin->compute_dmi(i,k,fmi,fmj);
}
}
if (me_flag) {
cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]);
if (rsq <= cut_me_2) {
lockpairspin->compute_me(i,k,fmi,fmj);
}
}
}
//post force
if (zeeman_flag) {
lockforcespin->compute_zeeman(i,fmi);
}
if (aniso_flag) {
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
lockforcespin->compute_anisotropy(i,spi,fmi);
}
//Replace the force by its new value
fm[i][0] = fmi[0];
fm[i][1] = fmi[1];
fm[i][2] = fmi[2];
}
#endif
/* ---------------------------------------------------------------------- */
void FixNVESpin::ComputeSpinInteractions()
{
int nflag,sortflag;
int nlocal = atom->nlocal;
int n_post_integrate = modify->n_post_integrate;
int n_pre_exchange = modify->n_pre_exchange;
int n_pre_neighbor = modify->n_pre_neighbor;
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
int n_end_of_step = modify->n_end_of_step;
bigint ntimestep;
ntimestep = update->ntimestep;
//int eflag = update->integrate->eflag;
//int vflag = update->integrate->vflag;
int eflag = 1;
int vflag = 0;
int pair_compute_flag = 1;
if (atom->sortfreq > 0) sortflag = 1;
else sortflag = 0;
if (n_post_integrate) modify->post_integrate();
timer->stamp(Timer::MODIFY);
// regular communication vs neighbor list rebuild
nflag = neighbor->decide();
if (nflag == 0) {
timer->stamp();
comm->forward_comm();
timer->stamp(Timer::COMM);
} else {
if (n_pre_exchange) {
timer->stamp();
modify->pre_exchange();
timer->stamp(Timer::MODIFY);
}
//if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
if (sortflag && ntimestep >= atom->nextsort) atom->sort();
comm->borders();
//if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
if (n_pre_neighbor) {
modify->pre_neighbor();
timer->stamp(Timer::MODIFY);
}
neighbor->build();
timer->stamp(Timer::NEIGH);
}
// force computations
// important for pair to come before bonded contributions
// since some bonded potentials tally pairwise energy/virial
// and Pair:ev_tally() needs to be called before any tallying
size_t nbytes;
nbytes = sizeof(double) * nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;
atom->avec->force_clear(0,nbytes);
timer->stamp();
if (n_pre_force) {
modify->pre_force(vflag);
timer->stamp(Timer::MODIFY);
}
if (pair_compute_flag) {
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
}
/*if (kspace_compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(Timer::KSPACE);
}*/
if (n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
// reverse communication of forces
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm)
{
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
int *type = atom->type;
int *mask = atom->mask;
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy;
double cp[3],g[3];
cp[0] = cp[1] = cp[2] = 0.0;
g[0] = g[1] = g[2] = 0.0;
fm2 = (fm[i][0]*fm[i][0])+(fm[i][1]*fm[i][1])+(fm[i][2]*fm[i][2]);
fmsq = sqrt(fm2);
energy = (sp[i][0]*fm[i][0])+(sp[i][1]*fm[i][1])+(sp[i][2]*fm[i][2]);
cp[0] = fm[i][1]*sp[i][2]-fm[i][2]*sp[i][1];
cp[1] = fm[i][2]*sp[i][0]-fm[i][0]*sp[i][2];
cp[2] = fm[i][0]*sp[i][1]-fm[i][1]*sp[i][0];
g[0] = sp[i][0]+cp[0]*dts;
g[1] = sp[i][1]+cp[1]*dts;
g[2] = sp[i][2]+cp[2]*dts;
g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts*dts;
g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts*dts;
g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts*dts;
g[0] /= (1+0.25*fm2*dts*dts);
g[1] /= (1+0.25*fm2*dts*dts);
g[2] /= (1+0.25*fm2*dts*dts);
sp[i][0] = g[0];
sp[i][1] = g[1];
sp[i][2] = g[2];
//Renormalization (may not be necessary)
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = 1.0/sqrt(msq);
sp[i][0] *= scale;
sp[i][1] *= scale;
sp[i][2] *= scale;
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::final_integrate()
{
double dtfm,msq,scale,fm2,fmsq,energy;
double cp[3],g[3];
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **sp = atom->sp;
double **fm = atom->fm;
double *rmass = atom->rmass;
double *mass = atom->mass;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
int *type = atom->type;
int *mask = atom->mask;
// update half v for all particles
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (rmass) dtfm = dtf / rmass[i];
else dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
}
}
// Advance half spins all particles
//See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011
if (extra == SPIN) {
#if defined SEQNEI
for (int i = nlocal-1; i >= 0; i--){
ComputeSpinInteractionsNei(i);
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#endif
#if defined SEQ
for (int i = nlocal-1; i >= 0; i--){
AdvanceSingleSpin(i,0.5*dts,sp,fm);
ComputeSpinInteractions();
}
#endif
}
}

90
src/fix_nve_spin.h Normal file
View File

@ -0,0 +1,90 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(nve/spin,FixNVESpin)
#else
#ifndef LMP_FIX_NVE_SPIN_H
#define LMP_FIX_NVE_SPIN_H
#include "fix_nve.h"
namespace LAMMPS_NS {
class FixNVESpin : public FixNVE {
public:
FixNVESpin(class LAMMPS *, int, char **);
virtual ~FixNVESpin();
void init();
virtual void initial_integrate(int);
void AdvanceSingleSpin(int, double, double **, double **);
virtual void final_integrate();
//#define SEQ
#define SEQNEI
void ComputeSpinInteractions();
void ComputeSpinInteractionsNei(int);
protected:
int extra;
double dts;
//double alpha_t;
#if defined SEQNEI
private:
int exch_flag, dmi_flag, me_flag;
int zeeman_flag, aniso_flag;
class PairSpin *lockpairspin;
double *spi, *fmi, *fmj; //Temp var. for compute
class FixForceSpin *lockforcespin;
#endif
//private:
//class FixSpinDamping *lockspindamping;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/sphere requires atom style sphere
Self-explanatory.
E: Fix nve/sphere update dipole requires atom attribute mu
An atom style with this attribute is needed.
E: Fix nve/sphere requires extended particles
This fix can only be used for particles of a finite size.
E: Fix nve/sphere dlm must be used with update dipole
The DLM algorithm can only be used in conjunction with update dipole.
*/

558
src/pair_spin.cpp Executable file
View File

@ -0,0 +1,558 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include "pair_spin.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "math_const.h"
#include "error.h"
#include "update.h"
#include <string.h>
//Add. lib. for full neighb. list
#include "neigh_request.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
exch_flag = 0;
dmi_flag = 0;
me_flag = 0;
}
/* ---------------------------------------------------------------------- */
PairSpin::~PairSpin()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cut_spin_exchange);
memory->destroy(J_1);
memory->destroy(J_2);
memory->destroy(J_2);
memory->destroy(cut_spin_dmi);
memory->destroy(DM);
memory->destroy(v_dmx);
memory->destroy(v_dmy);
memory->destroy(v_dmz);
memory->destroy(cut_spin_me);
memory->destroy(ME);
memory->destroy(v_mex);
memory->destroy(v_mey);
memory->destroy(v_mez);
memory->destroy(fmi);
memory->destroy(fmj);
memory->destroy(cutsq);
}
}
/* ---------------------------------------------------------------------- */
void PairSpin::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl,ecoul;
double xtmp,ytmp,ztmp;
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
double cut_ex_2,cut_dmi_2,cut_me_2;
double rsq,rd,delx,dely,delz;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **fm = atom->fm;
double *mumag = atom->mumag;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Pair spin computations
// Loop over neighbors of my itoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
//Loop on Neighbors
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
fmi[0] = fmi[1] = fmi[2] = 0.0;
fmj[0] = fmj[1] = fmj[2] = 0.0;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; //square or inter-atomic distance
itype = type[i];
jtype = type[j];
//Exchange interaction
if (exch_flag) {
cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
if (rsq <= cut_ex_2) {
compute_exchange(i,j,rsq,fmi,fmj);
}
}
//DM interaction
if (dmi_flag){
cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
if (rsq <= cut_dmi_2){
compute_dmi(i,j,fmi,fmj);
}
}
//ME interaction
if (me_flag){
cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
if (rsq <= cut_me_2){
compute_me(i,j,fmi,fmj);
}
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (newton_pair || j < nlocal) {
fm[j][0] += fmj[0];
fm[j][1] += fmj[1];
fm[j][2] += fmj[2];
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj)
{
int *type = atom->type;
int itype, jtype;
double **sp = atom->sp;
double dmix,dmiy,dmiz;
double Jex, ra;
itype = type[i];
jtype = type[j];
ra = rsq/J_3[itype][jtype]/J_3[itype][jtype];
Jex = 4.0*J_1[itype][jtype]*ra;
Jex *= (1.0-J_2[itype][jtype]*ra);
Jex *= exp(-ra);
fmi[0] += Jex*sp[j][0];
fmi[1] += Jex*sp[j][1];
fmi[2] += Jex*sp[j][2];
fmj[0] += Jex*sp[i][0];
fmj[1] += Jex*sp[i][1];
fmj[2] += Jex*sp[i][2];
}
/* ---------------------------------------------------------------------- */
void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj)
{
int *type = atom->type;
int itype, jtype;
double **sp = atom->sp;
double dmix,dmiy,dmiz;
itype = type[i];
jtype = type[j];
dmix = DM[itype][jtype]*v_dmx[itype][jtype];
dmiy = DM[itype][jtype]*v_dmy[itype][jtype];
dmiz = DM[itype][jtype]*v_dmz[itype][jtype];
fmi[0] += sp[j][1]*dmiz-sp[j][2]*dmiy;
fmi[1] += sp[j][2]*dmix-sp[j][0]*dmiz;
fmi[2] += sp[j][0]*dmiy-sp[j][1]*dmix;
fmj[0] -= sp[i][1]*dmiz-sp[i][2]*dmiy;
fmj[1] -= sp[i][2]*dmix-sp[i][0]*dmiz;
fmj[2] -= sp[i][0]*dmiy-sp[i][1]*dmix;
}
/* ---------------------------------------------------------------------- */
void PairSpin::compute_me(int i, int j, double *fmi, double *fmj)
{
int *type = atom->type;
int itype, jtype;
itype = type[i];
jtype = type[j];
double **sp = atom->sp;
double **x = atom->x;
double meix,meiy,meiz;
double rx, ry, rz, inorm;
rx = x[j][0] - x[i][0];
ry = x[j][1] - x[i][1];
rz = x[j][2] - x[i][2];
inorm = 1.0/sqrt(rx*rx+ry*ry+rz*rz);
rx *= inorm;
ry *= inorm;
rz *= inorm;
meix = v_mey[itype][jtype]*rz - v_mez[itype][jtype]*ry;
meiy = v_mez[itype][jtype]*rx - v_mex[itype][jtype]*rz;
meiz = v_mex[itype][jtype]*ry - v_mey[itype][jtype]*rx;
meix *= ME[itype][jtype];
meiy *= ME[itype][jtype];
meiz *= ME[itype][jtype];
fmi[0] += sp[j][1]*meiz - sp[j][2]*meiy;
fmi[1] += sp[j][2]*meix - sp[j][0]*meiz;
fmi[2] += sp[j][0]*meiy - sp[j][1]*meix;
fmj[0] -= sp[i][1]*meiz - sp[i][2]*meiy;
fmj[1] -= sp[i][2]*meix - sp[i][0]*meiz;
fmj[2] -= sp[i][0]*meiy - sp[i][1]*meix;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSpin::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cut_spin_exchange,n+1,n+1,"pair:cut_spin_exchange");
memory->create(J_1,n+1,n+1,"pair:J_1");
memory->create(J_2,n+1,n+1,"pair:J_2");
memory->create(J_3,n+1,n+1,"pair:J_3");
memory->create(cut_spin_dmi,n+1,n+1,"pair:cut_spin_dmi");
memory->create(DM,n+1,n+1,"pair:DM");
memory->create(v_dmx,n+1,n+1,"pair:DM_vector_x");
memory->create(v_dmy,n+1,n+1,"pair:DM_vector_y");
memory->create(v_dmz,n+1,n+1,"pair:DM_vector_z");
memory->create(cut_spin_me,n+1,n+1,"pair:cut_spin_me");
memory->create(ME,n+1,n+1,"pair:ME");
memory->create(v_mex,n+1,n+1,"pair:ME_vector_x");
memory->create(v_mey,n+1,n+1,"pair:ME_vector_y");
memory->create(v_mez,n+1,n+1,"pair:ME_vector_z");
memory->create(fmi,3,"pair:fmi");
memory->create(fmj,3,"pair:fmj");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSpin::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2)
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
if (strcmp(update->unit_style,"metal") != 0)
error->all(FLERR,"Spin simulations require metal unit style");
cut_spin_pair_global = force->numeric(FLERR,arg[0]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) {
cut_spin_exchange[i][j] = cut_spin_pair_global;
cut_spin_dmi[i][j] = cut_spin_pair_global;
cut_spin_me[i][j] = cut_spin_pair_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type spin pairs (only one for now)
------------------------------------------------------------------------- */
void PairSpin::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if (strcmp(arg[2],"exchange")==0){
if (narg != 7) error->all(FLERR,"Incorrect args in pair_style command");
exch_flag = 1;
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double rij = force->numeric(FLERR,arg[3]);
double J1 = force->numeric(FLERR,arg[4]);
double J2 = force->numeric(FLERR,arg[5]);
double J3 = force->numeric(FLERR,arg[6]);
double hbar = force->hplanck/MY_2PI;
J1 /= hbar;
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut_spin_exchange[i][j] = rij;
J_1[i][j] = J1;
J_2[i][j] = J2;
J_3[i][j] = J3;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
} else if (strcmp(arg[2],"dmi")==0) {
if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command");
dmi_flag = 1;
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double rij = force->numeric(FLERR,arg[3]);
double dm = force->numeric(FLERR,arg[4]);
double dmx = force->numeric(FLERR,arg[5]);
double dmy = force->numeric(FLERR,arg[6]);
double dmz = force->numeric(FLERR,arg[7]);
double inorm = 1.0/(dmx*dmx+dmy*dmy+dmz*dmz);
dmx *= inorm;
dmy *= inorm;
dmz *= inorm;
double hbar = force->hplanck/MY_2PI;
dm /= hbar;
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut_spin_dmi[i][j] = rij;
DM[i][j] = dm;
v_dmx[i][j] = dmx;
v_dmy[i][j] = dmy;
v_dmz[i][j] = dmz;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
} else if (strcmp(arg[2],"me")==0) {
if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command");
me_flag = 1;
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double rij = force->numeric(FLERR,arg[3]);
double me = force->numeric(FLERR,arg[4]);
double mex = force->numeric(FLERR,arg[5]);
double mey = force->numeric(FLERR,arg[6]);
double mez = force->numeric(FLERR,arg[7]);
double inorm = 1.0/(mex*mex+mey*mey+mez*mez);
mex *= inorm;
mey *= inorm;
mez *= inorm;
double hbar = force->hplanck/MY_2PI;
me /= hbar;
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut_spin_me[i][j] = rij;
DM[i][j] = me;
v_mex[i][j] = mex;
v_mey[i][j] = mey;
v_mez[i][j] = mez;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
} else error->all(FLERR,"Incorrect args in pair_style command");
//Check if Jex [][] still works for Ferrimagnetic exchange
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairSpin::init_style()
{
if (!atom->sp_flag || !atom->mumag_flag)
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
neighbor->request(this,instance_me);
#define FULLNEI
#if defined FULLNEI
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
#endif
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSpin::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return cut_spin_pair_global;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSpin::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
if (exch_flag){
fwrite(&J_1[i][j],sizeof(double),1,fp);
fwrite(&J_2[i][j],sizeof(double),1,fp);
fwrite(&J_3[i][j],sizeof(double),1,fp);
fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp);
}
if (dmi_flag) {
fwrite(&DM[i][j],sizeof(double),1,fp);
fwrite(&v_dmx[i][j],sizeof(double),1,fp);
fwrite(&v_dmy[i][j],sizeof(double),1,fp);
fwrite(&v_dmz[i][j],sizeof(double),1,fp);
fwrite(&cut_spin_dmi[i][j],sizeof(double),1,fp);
}
if (me_flag) {
fwrite(&ME[i][j],sizeof(double),1,fp);
fwrite(&v_mex[i][j],sizeof(double),1,fp);
fwrite(&v_mey[i][j],sizeof(double),1,fp);
fwrite(&v_mez[i][j],sizeof(double),1,fp);
fwrite(&cut_spin_me[i][j],sizeof(double),1,fp);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSpin::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&J_1[i][j],sizeof(double),1,fp);
fread(&J_2[i][j],sizeof(double),1,fp);
fread(&J_2[i][j],sizeof(double),1,fp);
fread(&cut_spin_exchange[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&J_1[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&J_2[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&J_3[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_spin_exchange[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSpin::write_restart_settings(FILE *fp)
{
fwrite(&cut_spin_pair_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSpin::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_spin_pair_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_spin_pair_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}

96
src/pair_spin.h Executable file
View File

@ -0,0 +1,96 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(pair/spin,PairSpin)
#else
#ifndef LMP_PAIR_SPIN_H
#define LMP_PAIR_SPIN_H
#include "pair.h"
namespace LAMMPS_NS {
class PairSpin : public Pair {
public:
PairSpin(class LAMMPS *);
virtual ~PairSpin();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void compute_exchange(int, int, double, double *, double *);
void compute_dmi(int, int, double *, double *);
void compute_me(int, int, double *, double *);
//Test for seq. integ.
//protected:
int exch_flag,dmi_flag,me_flag;
double cut_spin_pair_global;
double cut_spin_dipolar_global;
double **cut_spin_exchange; //cutting distance exchange
double **cut_spin_dmi; //cutting distance dmi
double **cut_spin_me; //cutting distance me
//Test for seq. integ.
protected:
double **J_1, **J_2, **J_3; //exchange coeffs Jij
//J1 in eV, J2 adim and J3 in Ang
double **DM;
double **v_dmx, **v_dmy, **v_dmz;//DMI coeffs
//DM int. in eV, v direction
double **ME;
double **v_mex, **v_mey, **v_mez;//ME coeffs
//ME in eV, v direction
double *fmi, *fmj; //Temp var. in compute
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_style command
Self-explanatory.
E: Spin simulations require metal unit style
Self-explanatory.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair spin requires atom attributes sp, mumag
The atom style defined does not have these attributes.
*/