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

First commit for the SPIN package. 

Changes to come:
-Exchange interaction computation to check (loop on neighbors),
-Temperature/random fluctuations to correct (effects too strong),
-Physical results to check,
-Add final interactions (DMI, ME, Dipolar),
-Compute spin temperature (Nurdin and Ma formslisms),
-Work on MPI parallelization,
-Ewald sums to implement (see with Stan's pakage),
-See for prefered magnetic axis (Mitchell's idea),
This commit is contained in:
julient31 2017-05-12 11:43:30 -06:00
parent d5ec76290b
commit 9b0f8a0c55
18 changed files with 2801 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

934
src/SPIN/atom_vec_spin.cpp Normal file
View File

@ -0,0 +1,934 @@
/* ----------------------------------------------------------------------
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;
comm_x_only = 0;
comm_f_only = 1;
size_forward = 6;
size_reverse = 3;
size_border = 11;
size_velocity = 3;
size_data_atom = 9;
size_data_vel = 4;
xcol_data = 4;
forceclearflag = 1;
atom->mumag_flag = atom->sp_flag = 1;
//atom->rmass_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];
}
} 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];
}
}
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++] = 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++] = 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];
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];
}
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++];
}
}
/* ---------------------------------------------------------------------- */
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++];
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++];
}
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];
}
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++];
}
}
/* ---------------------------------------------------------------------- */
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] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
buf[i][10] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
buf[i][11] = 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];
return 4;
}
/* ----------------------------------------------------------------------
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],
(int) ubuf(buf[i][9]).i,(int) ubuf(buf[i][10]).i,
(int) ubuf(buf[i][11]).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",buf[0],buf[1],buf[2],buf[3]);
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;
}
//Test force clear in spin
void AtomVecSpin::force_clear(int n, size_t nbytes)
{
memset(&atom->fm[0][0],0,3*nbytes);
}

90
src/SPIN/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.
*/

134
src/SPIN/compute_spin.cpp Normal file
View File

@ -0,0 +1,134 @@
/* ----------------------------------------------------------------------
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"
using namespace LAMMPS_NS;
using namespace MathSpecial;
/* ---------------------------------------------------------------------- */
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 = 5;
extvector = 0;
init();
allocate();
}
/* ---------------------------------------------------------------------- */
ComputeSpin::~ComputeSpin()
{
memory->destroy(mag);
}
/* ---------------------------------------------------------------------- */
void ComputeSpin::init()
{
}
/* ---------------------------------------------------------------------- */
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;
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;
int nlocal = atom->nlocal;
// compute total magnetization
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];
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(&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]));
// compute total magnetic energy
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (atom->mumag_flag && atom->sp_flag) {
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];
}
else error->all(FLERR,"Compute spin/compute declared magnetic quantities (sp and mumag flags)");
}
}
MPI_Allreduce(&magenergy,&magenergytot,1,MPI_DOUBLE,MPI_SUM,world);
vector[0] = magtot[0];
vector[1] = magtot[1];
vector[2] = magtot[2];
vector[3] = magtot[3];
vector[4] = magenergytot;
}
/* ----------------------------------------------------------------------
free and reallocate arrays
------------------------------------------------------------------------- */
void ComputeSpin::allocate()
{
memory->destroy(mag);
memory->create(mag,4,"compute/spin:mag");
memory->create(magtot,5,"compute/spin:mag");
vector = magtot;
}

67
src/SPIN/compute_spin.h Normal file
View File

@ -0,0 +1,67 @@
/* -*- 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;
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.
*/

229
src/SPIN/fix_force_spin.cpp Normal file
View File

@ -0,0 +1,229 @@
/* ----------------------------------------------------------------------
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"
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;
if (strcmp(arg[3],"zeeman") == 0) {
if (narg != 8) error->all(FLERR,"Illegal fix zeeman command");
style = ZEEMAN;
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;
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");
//printf("test field in creator: H=%g, Hx=%g, Hy=%g, Hz=%g \n",H_field,Hx,Hy,Hz);
degree2rad = MY_PI/180.0;
time_origin = update->ntimestep;
eflag = 0;
emag = 0.0;
}
/* ---------------------------------------------------------------------- */
FixForceSpin::~FixForceSpin()
{
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();
}
/* ---------------------------------------------------------------------- */
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;
if (style == ZEEMAN) {
for (int i = 0; i < nlocal; i++) {
fm[i][0] += mumag[i]*xmag;
fm[i][1] += mumag[i]*ymag;
fm[i][2] += mumag[i]*zmag;
// emag -= (sp[i][0]*xmag + sp[i][1]*ymag + sp[i][2]*zmag);
}
}
if (style == ANISOTROPY) {
for (int i = 0; i < nlocal; i++) {
scalar = Kax*sp[i][0] + Kay*sp[i][1] + Kaz*sp[i][2];
fm[i][0] -= Ka*scalar*Kax;
fm[i][1] -= Ka*scalar*Kay;
fm[i][2] -= Ka*scalar*Kaz;
//emag -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
}
}
//printf("test force. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//printf("Field force compute, fm[0][2]=%g \n",fm[0][2]);
}
/* ---------------------------------------------------------------------- */
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"
//(keeping set_magneticforce in case of time--dependent mag. field implementation)
void FixForceSpin::set_magneticforce()
{
if (style == ZEEMAN) {
xmag = H_field*Hx;
ymag = H_field*Hy;
zmag = H_field*Hz;
}
}
/* ----------------------------------------------------------------------
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;
}

86
src/SPIN/fix_force_spin.h Normal file
View File

@ -0,0 +1,86 @@
/* -*- 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();
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;
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.
*/

204
src/SPIN/fix_langevin_spin.cpp Normal file
View File

@ -0,0 +1,204 @@
/* ----------------------------------------------------------------------
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 = alpha_t/(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;
// apply transverse magnetic damping to spins
// add the damping to the effective field
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
sx = sp[i][0];//Getting Mag. and Mag. force components
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;
}
//printf("test damping. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//apply thermal effects
//add random field to fm
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
rx = sigma*random->gaussian();//Drawing random distributions
ry = sigma*random->gaussian();
rz = sigma*random->gaussian();
//rx = sigma*(random->uniform() - 0.5);
//ry = sigma*(random->uniform() - 0.5);
//rz = sigma*(random->uniform() - 0.5);
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;
}
//printf("test rand var: %g, sigma=%g \n",(random->uniform()-0.5),sigma);
//printf("test rand var: %g, sigma=%g \n",random->gaussian(),sigma);
//printf("test random 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//printf("test dt: %g, sigma=%g \n",dts,random->gaussian(),sigma);
}
/* ---------------------------------------------------------------------- */
void FixLangevinSpin::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}

106
src/SPIN/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.
*/

211
src/SPIN/fix_nve_spin.cpp Normal file
View File

@ -0,0 +1,211 @@
/* ----------------------------------------------------------------------
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 "fix_damping_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"
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");
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::init()
{
FixNVE::init();
dts = update->dt;
/*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 *rmass = atom->rmass;
double *mass = atom->mass;
double **sp = atom->sp;
double **fm = atom->fm;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
int *type = atom->type;
int *mask = atom->mask;
// 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 half 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];
}
}
// update sp for all particles
if (extra == SPIN) {
// Advance spins
//See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (sp[i][3] > 0.0) {
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];
//cp[0] = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
//cp[1] = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
//cp[2] = sp[i][0]*fm[i][1]-sp[i][1]*fm[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+(fmsq*dts*0.5)*(fmsq*dts*0.5));
g[1] /= (1+(fmsq*dts*0.5)*(fmsq*dts*0.5));
g[2] /= (1+(fmsq*dts*0.5)*(fmsq*dts*0.5));
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;
// printf("test fix integ. 1;i=%d, fx=%g, fy=%g, fz=%g \n",i,fm[i][0],fm[i][1],fm[i][2]);
//printf("test fix integ.; i=%d, sx=%g, sy=%g, sz=%g, norm=%g \n",i,sp[i][0],sp[i][1],sp[i][2],scale);
}
}
//printf("test fix integ. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
}
/* ---------------------------------------------------------------------- */
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 *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 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];
}
}
// 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];
}
}
}

76
src/SPIN/fix_nve_spin.h Normal file
View File

@ -0,0 +1,76 @@
/* -*- 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 {
friend class FixSpinDamping;
public:
FixNVESpin(class LAMMPS *, int, char **);
virtual ~FixNVESpin() {}
void init();
virtual void initial_integrate(int);
virtual void final_integrate();
protected:
int extra;
double dts;
double alpha_t;
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.
*/

339
src/SPIN/pair_spin.cpp Executable file
View File

@ -0,0 +1,339 @@
/* ----------------------------------------------------------------------
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>
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
}
/* ---------------------------------------------------------------------- */
PairSpin::~PairSpin()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cut_spin_exchange);
memory->destroy(cut_spin_dipolar);
memory->destroy(J_1);
memory->destroy(J_2);
memory->destroy(J_2);
memory->destroy(cutsq);
}
}
/* ---------------------------------------------------------------------- */
void PairSpin::compute(int eflag, int vflag)
{
double **x = atom->x;
double **fm = atom->fm;
double *mumag = atom->mumag;
double **sp = atom->sp;
int *type = atom->type;
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,fmix,fmiy,fmiz,fmjx,fmjy,fmjz,omx,omy,omz;
double rsq,rd,delx,dely,delz;
int *ilist,*jlist,*numneigh,**firstneigh;
double cut, Jex, ra;
double evdwl,ecoul;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// Pair spin computations
// Loop over neighbors of my atoms
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];
//Exchange interaction
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
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
rd = sqrt(rsq); //Inter-atomic distance
cut = cut_spin_exchange_global;
if (rd <= cut) {
itype = type[i];
jtype = type[j];
fmix = fmiy = fmiz = 0.0;
fmjx = fmjy = fmjz = 0.0;
ra = (rd/J_3[itype][jtype])*(rd/J_3[itype][jtype]);
Jex = 4.0*J_1[itype][jtype]*ra;
Jex *= (1.0-J_2[itype][jtype]*ra);
Jex *= exp(-ra);
Jex *= mumag[ii]*mumag[jj];
fmix = Jex*sp[j][0];
fmiy = Jex*sp[j][1];
fmiz = Jex*sp[j][2];
//fmjx = Jex*sp[i][0];
//fmjy = Jex*sp[i][1];
//fmjz = Jex*sp[i][2];
}
fm[i][0] += fmix;
fm[i][1] += fmiy;
fm[i][2] += fmiz;
//fm[j][0] += fmjx;
//fm[j][1] += fmjy;
//fm[j][2] += fmjz;
}
}
//printf("vals exchange: Jx=%g, Jy=%g, Jz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//printf("test exchange. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
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(cut_spin_dipolar,n+1,n+1,"pair:cut_spin_dipolar");
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(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,"electron") == 0)
error->all(FLERR,"Cannot (yet) use 'electron' units with spins");
cut_spin_exchange_global = force->numeric(FLERR,arg[0]);
if (narg == 1) cut_spin_dipolar_global = cut_spin_exchange_global;
else cut_spin_dipolar_global = force->numeric(FLERR,arg[1]);
// 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_exchange_global;
cut_spin_dipolar[i][j] = cut_spin_dipolar_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type spin pairs (only one for now)
------------------------------------------------------------------------- */
void PairSpin::coeff(int narg, char **arg)
{
if (narg != 5)
error->all(FLERR,"Incorrect number of args for pair spin coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double J1 = force->numeric(FLERR,arg[2]);
double J2 = force->numeric(FLERR,arg[3]);
double J3 = force->numeric(FLERR,arg[4]);
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++) {
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 for spinpair coefficients");
//Simple (Anti)Ferromagnetic exchange for now.
//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);
}
/* ----------------------------------------------------------------------
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_exchange_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]) {
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);
}
}
}
/* ----------------------------------------------------------------------
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_exchange_global,sizeof(double),1,fp);
fwrite(&cut_spin_dipolar_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_exchange_global,sizeof(double),1,fp);
fread(&cut_spin_dipolar_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_spin_dipolar_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}

75
src/SPIN/pair_spin.h Executable file
View File

@ -0,0 +1,75 @@
/* -*- 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 *);
protected:
double cut_spin_exchange_global, cut_spin_dipolar_global; //Global cutting distance
double **cut_spin_exchange; //cutting distance for each exchange interaction
double **cut_spin_dipolar; //cutting distance for the dipolar interaction
double **J_1, **J_2, **J_3; //coefficients for computing the exchange interaction Jij
//J1 is an energy (in eV), J2 is adim and J3 is a distance (in Ang)
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_style command
Self-explanatory.
E: Cannot (yet) use 'electron' units with spins
This feature is not yet supported.
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.
*/

12
src/atom.cpp
View File

@ -98,6 +98,10 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
rho = drho = e = de = cv = NULL;
vest = NULL;
// USER-SPIN
mumag = NULL;
sp = fm = NULL;
// USER-DPD
uCond = uMech = uChem = uCG = uCGnew = NULL;
@ -167,6 +171,9 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
omega_flag = torque_flag = angmom_flag = 0;
radius_flag = rmass_flag = 0;
ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
//Magnetic flags
sp_flag = mumag_flag = 0;
vfrac_flag = 0;
spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0;
@ -421,6 +428,9 @@ void Atom::create_avec(const char *style, int narg, char **arg, int trysuffix)
radius_flag = rmass_flag = 0;
ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
//Magnetic flags
sp_flag = mumag_flag = 0;
vfrac_flag = 0;
spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0;
cs_flag = csforce_flag = vforce_flag = etag_flag = 0;
@ -494,7 +504,7 @@ AtomVec *Atom::new_avec(const char *style, int trysuffix, int &sflag)
AtomVecCreator avec_creator = (*avec_map)[style];
return avec_creator(lmp);
}
//printf("test entries function: %s, %d, %d \n ",style, trysuffix, &sflag);
error->all(FLERR,"Unknown atom style");
return NULL;
}

7
src/atom.h
View File

@ -61,6 +61,10 @@ class Atom : protected Pointers {
double *radius,*rmass;
int *ellipsoid,*line,*tri,*body;
// SPIN package
double *mumag, **sp;
double **fm;
// PERI package
double *vfrac,*s0;
@ -146,6 +150,9 @@ class Atom : protected Pointers {
int rho_flag,e_flag,cv_flag,vest_flag;
int dpd_flag,edpd_flag,tdpd_flag;
//USER-SPIN package
int mumag_flag,sp_flag;
// USER-SMD package
int smd_flag;

147
src/dump_custom.cpp
View File

@ -41,6 +41,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS,
IX,IY,IZ,
VX,VY,VZ,FX,FY,FZ,
Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER,
MUMAG,SPX,SPY,SPZ,SP,
OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
TQX,TQY,TQZ,
COMPUTE,FIX,VARIABLE,INAME,DNAME};
@ -83,8 +84,8 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
pack_choice = new FnPtrPack[nfield];
vtype = new int[nfield];
memory->create(field2index,nfield,"dump:field2index");
memory->create(argindex,nfield,"dump:argindex");
field2index = new int[nfield];
argindex = new int[nfield];
buffer_allow = 1;
buffer_flag = 1;
@ -201,8 +202,8 @@ DumpCustom::~DumpCustom()
delete [] pack_choice;
delete [] vtype;
memory->destroy(field2index);
memory->destroy(argindex);
delete [] field2index;
delete [] argindex;
delete [] idregion;
memory->destroy(thresh_array);
@ -245,15 +246,11 @@ DumpCustom::~DumpCustom()
for (int i = 1; i <= ntypes; i++) delete [] typenames[i];
delete [] typenames;
if(vformat) {
for (int i = 0; i < size_one; i++) delete [] vformat[i];
delete [] vformat;
}
for (int i = 0; i < size_one; i++) delete [] vformat[i];
delete [] vformat;
if(format_column_user) {
for (int i = 0; i < size_one; i++) delete [] format_column_user[i];
delete [] format_column_user;
}
for (int i = 0; i < size_one; i++) delete [] format_column_user[i];
delete [] format_column_user;
delete [] columns;
}
@ -857,6 +854,36 @@ int DumpCustom::count()
for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i];
ptr = dchoose;
nstride = 1;
} else if (thresh_array[ithresh] == MUMAG) {//Magnetic properties
if (!atom->mumag_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = atom->mumag;
nstride = 1;
} else if (thresh_array[ithresh] == SPX) {
if (!atom->sp_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = &atom->sp[0][0];
nstride = 4;
} else if (thresh_array[ithresh] == SPY) {
if (!atom->sp_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = &atom->sp[0][1];
nstride = 4;
} else if (thresh_array[ithresh] == SPZ) {
if (!atom->sp_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = &atom->sp[0][2];
nstride = 4;
} else if (thresh_array[ithresh] == SP) {
if (!atom->sp_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
ptr = &atom->sp[0][3];
nstride = 4;
} else if (thresh_array[ithresh] == OMEGAX) {
if (!atom->omega_flag)
error->all(FLERR,
@ -1285,6 +1312,35 @@ int DumpCustom::parse_fields(int narg, char **arg)
pack_choice[i] = &DumpCustom::pack_mu;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"mumag") == 0) {//Magnetic properties
if (!atom->mumag_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[i] = &DumpCustom::pack_mumag;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"spx") == 0) {
strcpy(arg[iarg],"vx");
if (!atom->sp_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[i] = &DumpCustom::pack_spx;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"spy") == 0) {
strcpy(arg[iarg],"vy");
if (!atom->sp_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[i] = &DumpCustom::pack_spy;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"spz") == 0) {
strcpy(arg[iarg],"vz");
if (!atom->sp_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[i] = &DumpCustom::pack_spz;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"sp") == 0) {
if (!atom->sp_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
pack_choice[i] = &DumpCustom::pack_sp;
vtype[i] = DOUBLE;
} else if (strcmp(arg[iarg],"radius") == 0) {
if (!atom->radius_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
@ -1787,6 +1843,13 @@ int DumpCustom::modify_param(int narg, char **arg)
else if (strcmp(arg[1],"muy") == 0) thresh_array[nthresh] = MUY;
else if (strcmp(arg[1],"muz") == 0) thresh_array[nthresh] = MUZ;
else if (strcmp(arg[1],"mu") == 0) thresh_array[nthresh] = MU;
//Magnetic quantities
else if (strcmp(arg[1],"mumag") == 0) thresh_array[nthresh] = MUMAG;
else if (strcmp(arg[1],"spx") == 0) thresh_array[nthresh] = SPX;
else if (strcmp(arg[1],"spy") == 0) thresh_array[nthresh] = SPY;
else if (strcmp(arg[1],"spz") == 0) thresh_array[nthresh] = SPZ;
else if (strcmp(arg[1],"sp") == 0) thresh_array[nthresh] = SP;
else if (strcmp(arg[1],"radius") == 0) thresh_array[nthresh] = RADIUS;
else if (strcmp(arg[1],"diameter") == 0) thresh_array[nthresh] = DIAMETER;
@ -2701,6 +2764,66 @@ void DumpCustom::pack_mu(int n)
}
}
/* ---------------------------------------------------------------------- */
//Magnetic quantities
void DumpCustom::pack_mumag(int n)
{
double *mumag = atom->mumag;
for (int i = 0; i < nchoose; i++) {
buf[n] = mumag[clist[i]];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_spx(int n)
{
double **sp = atom->sp;
for (int i = 0; i < nchoose; i++) {
buf[n] = sp[clist[i]][0];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_spy(int n)
{
double **sp = atom->sp;
for (int i = 0; i < nchoose; i++) {
buf[n] = sp[clist[i]][1];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_spz(int n)
{
double **sp = atom->sp;
for (int i = 0; i < nchoose; i++) {
buf[n] = sp[clist[i]][2];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_sp(int n)
{
double **sp = atom->sp;
for (int i = 0; i < nchoose; i++) {
buf[n] = sp[clist[i]][3];
n += size_one;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_radius(int n)

5
src/dump_custom.h
View File

@ -178,6 +178,11 @@ class DumpCustom : public Dump {
void pack_muy(int);
void pack_muz(int);
void pack_mu(int);
void pack_mumag(int); //Magnetic quantities
void pack_spx(int);
void pack_spy(int);
void pack_spz(int);
void pack_sp(int);
void pack_radius(int);
void pack_diameter(int);

90
src/set.cpp
View File

@ -44,7 +44,7 @@ using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,
DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
MESO_E,MESO_CV,MESO_RHO,EDPD_TEMP,EDPD_CV,CC,SMD_MASS_DENSITY,
SMD_CONTACT_RADIUS,DPDTHETA,INAME,DNAME};
@ -215,6 +215,34 @@ void Set::command(int narg, char **arg)
setrandom(DIPOLE_RANDOM);
iarg += 3;
} else if (strcmp(arg[iarg],"spin") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else xvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else yvalue = force->numeric(FLERR,arg[iarg+3]);
if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) varparse(arg[iarg+4],4);
else zvalue = force->numeric(FLERR,arg[iarg+4]);
if (!atom->sp_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(SPIN);
iarg += 5;
} else if (strcmp(arg[iarg],"spin/random") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
dvalue = force->numeric(FLERR,arg[iarg+2]);
if (!atom->sp_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
if (dvalue <= 0.0)
error->all(FLERR,"Invalid dipole length in set command");
setrandom(SPIN_RANDOM);
iarg += 3;
} else if (strcmp(arg[iarg],"quat") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
@ -821,6 +849,21 @@ void Set::set(int keyword)
mu[i][2]*mu[i][2]);
}
// set magnetic moments
else if (keyword == SPIN) {
double **sp = atom->sp;
double *mumag = atom->mumag;
double sp_norm = sqrt(xvalue*xvalue+yvalue*yvalue+zvalue*zvalue);
sp[i][0] = xvalue/sp_norm;
sp[i][1] = yvalue/sp_norm;
sp[i][2] = zvalue/sp_norm;
sp[i][3] = sqrt(sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1] +
sp[i][2]*sp[i][2]); //Should be 1 for atomic spins
mumag[i] = sp_norm;
}
// set quaternion orientation of ellipsoid or tri or body particle
// set quaternion orientation of ellipsoid or tri or body particle
// enforce quat rotation vector in z dir for 2d systems
@ -981,6 +1024,51 @@ void Set::setrandom(int keyword)
}
}
// set spin moments to random orientations in 3d or 2d
// spin length is fixed to unity
} else if (keyword == SPIN_RANDOM) {
double **sp = atom->sp;
double *mumag = atom->mumag;
int nlocal = atom->nlocal;
double sp_sq,scale;
if (domain->dimension == 3) {
for (i = 0; i < nlocal; i++)
if (select[i]) {
random->reset(seed,x[i]);
sp[i][0] = random->uniform() - 0.5;
sp[i][1] = random->uniform() - 0.5;
sp[i][2] = random->uniform() - 0.5;
sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1] + sp[i][2]*sp[i][2];
scale = 1.0/sqrt(sp_sq);
sp[i][0] *= scale;
sp[i][1] *= scale;
sp[i][2] *= scale;
sp[i][3] = sqrt(sp_sq);
mumag[i] = dvalue;
count++;
}
} else {
for (i = 0; i < nlocal; i++)
if (select[i]) {
random->reset(seed,x[i]);
sp[i][0] = random->uniform() - 0.5;
sp[i][1] = random->uniform() - 0.5;
sp[i][2] = 0.0;
sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1];
scale = 1.0/sqrt(sp_sq);
sp[i][0] *= scale;
sp[i][1] *= scale;
sp[i][3] = sqrt(sp_sq);
mumag[i] = dvalue;
count++;
}
}
// set quaternions to random orientations in 3d and 2d
} else if (keyword == QUAT_RANDOM) {

3
src/verlet.cpp
View File

@ -75,6 +75,7 @@ void Verlet::init()
torqueflag = extraflag = 0;
if (atom->torque_flag) torqueflag = 1;
if (atom->avec->forceclearflag) extraflag = 1;
if (atom->mumag_flag) extraflag = 1;
// orthogonal vs triclinic simulation box
@ -385,6 +386,8 @@ void Verlet::force_clear()
if (nbytes) {
memset(&atom->f[0][0],0,3*nbytes);
//test memset for fm
//memset(&atom->fm[0][0],0,3*nbytes);
if (torqueflag) memset(&atom->torque[0][0],0,3*nbytes);
if (extraflag) atom->avec->force_clear(0,nbytes);
}