diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp new file mode 100644 index 0000000000..6d12a7d4e3 --- /dev/null +++ b/src/SPIN/atom_vec_spin.cpp @@ -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 +#include +#include +#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 (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); +} + + diff --git a/src/SPIN/atom_vec_spin.h b/src/SPIN/atom_vec_spin.h new file mode 100644 index 0000000000..f4b13926f0 --- /dev/null +++ b/src/SPIN/atom_vec_spin.h @@ -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. + +*/ diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp new file mode 100644 index 0000000000..9ff12ac3dd --- /dev/null +++ b/src/SPIN/compute_spin.cpp @@ -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 +#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; +} + diff --git a/src/SPIN/compute_spin.h b/src/SPIN/compute_spin.h new file mode 100644 index 0000000000..c4011d5ded --- /dev/null +++ b/src/SPIN/compute_spin.h @@ -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. + +*/ diff --git a/src/SPIN/fix_force_spin.cpp b/src/SPIN/fix_force_spin.cpp new file mode 100644 index 0000000000..0163b32eb9 --- /dev/null +++ b/src/SPIN/fix_force_spin.cpp @@ -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 +#include +#include +#include +#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; +} diff --git a/src/SPIN/fix_force_spin.h b/src/SPIN/fix_force_spin.h new file mode 100644 index 0000000000..75805a7734 --- /dev/null +++ b/src/SPIN/fix_force_spin.h @@ -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. + +*/ diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp new file mode 100644 index 0000000000..ee8f5c1c49 --- /dev/null +++ b/src/SPIN/fix_langevin_spin.cpp @@ -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 +#include +#include +#include +#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); +} + diff --git a/src/SPIN/fix_langevin_spin.h b/src/SPIN/fix_langevin_spin.h new file mode 100644 index 0000000000..ffd94745e2 --- /dev/null +++ b/src/SPIN/fix_langevin_spin.h @@ -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. + +*/ diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp new file mode 100644 index 0000000000..4b9597e988 --- /dev/null +++ b/src/SPIN/fix_nve_spin.cpp @@ -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 +#include +#include +#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]; + } + } + +} diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h new file mode 100644 index 0000000000..55fffa7fb0 --- /dev/null +++ b/src/SPIN/fix_nve_spin.h @@ -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. + + +*/ diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp new file mode 100755 index 0000000000..975dcfc012 --- /dev/null +++ b/src/SPIN/pair_spin.cpp @@ -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 +#include +#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 + +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); +} diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h new file mode 100755 index 0000000000..abdc1978c1 --- /dev/null +++ b/src/SPIN/pair_spin.h @@ -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. + +*/ diff --git a/src/atom.cpp b/src/atom.cpp index ae36a8884b..a847719e7e 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -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; } diff --git a/src/atom.h b/src/atom.h index abfa6e5eb5..8e7f9811ac 100644 --- a/src/atom.h +++ b/src/atom.h @@ -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; diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 34ab218204..04ffc3b29b 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -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) diff --git a/src/dump_custom.h b/src/dump_custom.h index 1420d69b9b..c96148f275 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -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); diff --git a/src/set.cpp b/src/set.cpp index 77cae4cb09..6758b74445 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -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) { diff --git a/src/verlet.cpp b/src/verlet.cpp index 019f3f2f05..7bac6d7647 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -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); }