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

Implemetation of SeqNei Algo 1 

Still Seq and SeqNei versions
Loop on Neigh in SeqNei not working yet
This commit is contained in:
julient31 2017-06-19 10:40:36 -06:00
parent 4d375e72f0
commit bf5b3f96e9
10 changed files with 513 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
in.cobalt
View File

@ -10,8 +10,8 @@ units metal
dimension 3
#setting boundary conditions. (p for periodic, f for fixed, ...)
#boundary p p p
boundary f f f
boundary p p p
#boundary f f f
#setting atom_style, defines what can of atoms to use in simulation (atomic, molecule, angle, dipole, ...)
atom_style spin
@ -20,30 +20,35 @@ atom_style spin
atom_modify map array
#newton off for pair spin in SEQNEI
#newton off off
###########################
#######Create atoms########
###########################
#Lattice constant of fcc Cobalt
#lattice fcc 3.54
lattice sc 2.50
lattice fcc 3.54
#lattice sc 2.50
#Test Kagome
#variable a equal sqrt(3.0)
#variable d equal 1.0/(2.0*sqrt(3.0)+1)
#variable p_1 equal $d*sqrt(3.0)
#variable p_2 equal $d*2.2*sqrt(3.0)
#variable a equal sqrt(3.0)/8.0
#variable b equal 3.0*sqrt(3.0)/8.0
#variable c equal sqrt(3.0)/4.0
#lattice custom 1.0 a1 2.0 0.0 0.0 &
# a2 1.0 $a 0.0 &
# a3 0.0 0.0 1.0 &
# basis ${p_1} ${p_2} 0.0 &
# basis ${p_2} ${p_2} 0.0 &
# basis 0.5 0.0 0.0
#lattice custom 2.5 a1 1.0 0.0 0.0 &
# a2 0.0 1.0 0.0 &
# a3 0.0 0.0 1.0 &
# basis 0.0 $a 0.0 &
# basis 0.25 $a 0.0 &
# basis 0.375 0.0 0.0 &
# basis 0.25 $b 0.0 &
# basis 0.5 $b 0.0 &
# basis 0.625 $c 0.0
#Defining a geometric region of space. Sets ID(user's choice), style(block, sphere, ...), then, args depends on the style chosen
#(for block, one has x0, xf, y0, yf, z0, zf, in distance units)
region box block 0.0 5.0 0.0 5.0 0.0 1.0
region box block 0.0 2.0 0.0 2.0 0.0 2.0
#Creating a simulation box based on the specified region. Entries: number of atom types and box ref.
create_box 1 box
@ -65,27 +70,18 @@ create_atoms 1 box
mass 1 1.0
#set group all mass 1.0
#Setting spins orientation and moment
set group all spin/random 11 1.72
#set group all spin 1.72 1.0 0.0 0.0
#set group all spin/random 11 1.72
set group all spin 1.72 1.0 0.0 0.0
#Magnetic exchange interaction coefficient for bulk fcc Cobalt
#pair_style pair/spin/exchange 4.0
#type i and j | J1 | J2 | J3
#pair_coeff * * 0.0446928 0.003496 1.4885
#pair_coeff * * 0.0 0.003496 1.4885
#pair_style pair/spin/dmi 4.0
# type i and j | DM (in eV) | Directions
#pair_coeff * * 0.001 0.0 0.0 1.0
#pair_style hybrid/overlay pair/spin/exchange 4.0 pair/spin/dmi 4.0
pair_style pair/spin 4.0
#type i and j | interaction type | cutoff | J1 | J2 | J3
pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885
pair_coeff * * dmi 2.6 0.01 1.0 0.0 0.0
#type i and j | interaction type | cutoff | J1 (eV) | J2 (adim) | J3 (Ang) (for Exchange)
pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885
#type i and j | interaction type | cutoff | Int (eV) | [dx,dy,dz] (for DMI and ME)
#pair_coeff * * dmi 2.6 0.01 1.0 0.0 0.0
#pair_coeff * * me 2.6 0.01 1.0 1.0 1.0
#Fix Langevin spins (merging damping and temperature)
#Defines a cutof distance for the neighbors.
#Define a skin distance, update neigh list every
#neighbor 1.0 bin
#neigh_modify every 10 check yes delay 20
neighbor 0.0 bin
@ -93,13 +89,15 @@ neigh_modify every 1 check no delay 0
#Magnetic field fix
#Type | Intensity (T or eV) | Direction
fix 1 all force/spin zeeman 0.0 0.0 0.0 1.0
fix 1 all force/spin zeeman 10.0 0.0 0.0 1.0
#fix 1 all force/spin anisotropy 0.001 0.0 0.0 1.0
#fix 1 all force/spin anisotropy -0.1 0.0 0.0 1.0
#fix 1 all force/spin anisotropy 0.1 0.0 1.0 0.0
#Fix Langevin spins (merging damping and temperature)
#Temp | Alpha_trans | Alpha_long | Seed
fix 2 all langevin/spin 0.0 0.1 0.0 21
#fix 2 all langevin/spin 0.0 0.0 0.0 21
#fix 2 all langevin/spin 0.0 0.1 0.0 21
fix 2 all langevin/spin 0.0 0.0 0.0 21
#Magnetic integration fix
fix 3 all nve/spin
@ -107,13 +105,13 @@ fix 3 all nve/spin
#compute real time, total magnetization, magnetic energy, and spin temperature
#Iteration | Time | Mx | My | Mz | |M| | Em | Tm
compute mag all compute/spin
fix outmag all ave/time 1 1 50 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag.dat
fix outmag all ave/time 1 1 50 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag_seqnei_test.dat
#Defining a computation that will be performed on a group of atoms.
#Entries: ID(user assigned), group-ID(group of atoms to peform the sim on), style(temp, pe, ...), args
#Setting the timestep for the simulation
timestep 0.00005
timestep 0.00001
##################
#######run########
@ -123,6 +121,6 @@ timestep 0.00005
dump 1 all custom 50 dump_spin.lammpstrj type x y z spx spy spz
#Running the simulations for N timesteps
run 60000
#run 10
#run 200000
run 1

1
src/SPIN/atom_vec_spin.cpp
View File

@ -46,6 +46,7 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
forceclearflag = 1;
atom->mumag_flag = atom->sp_flag = 1;
}

68
src/SPIN/fix_force_spin.cpp
View File

@ -26,6 +26,7 @@
#include "math_const.h"
#include "error.h"
#include "force.h"
#include "memory.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -59,10 +60,13 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
Hx = Hy = Hz = 0.0;
Ka = 0.0;
Kax = Kay = Kaz = 0.0;
zeeman_flag = aniso_flag = 0;
if (strcmp(arg[3],"zeeman") == 0) {
if (narg != 8) error->all(FLERR,"Illegal fix zeeman command");
style = ZEEMAN;
zeeman_flag = 1;
H_field = force->numeric(FLERR,arg[4]);
Hx = force->numeric(FLERR,arg[5]);
Hy = force->numeric(FLERR,arg[6]);
@ -71,6 +75,7 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
} else if (strcmp(arg[3],"anisotropy") == 0) {
if (narg != 8) error->all(FLERR,"Illegal fix anisotropy command");
style = ANISOTROPY;
aniso_flag = 1;
Ka = force->numeric(FLERR,arg[4]);
Kax = force->numeric(FLERR,arg[5]);
Kay = force->numeric(FLERR,arg[6]);
@ -88,6 +93,8 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
FixForceSpin::~FixForceSpin()
{
memory->destroy(spi);
memory->destroy(fmi);
delete [] magstr;
}
@ -133,7 +140,10 @@ void FixForceSpin::init()
// set magnetic field components once and for all
if (varflag == CONSTANT) set_magneticforce();
memory->create(spi,3,"forcespin:spi");
memory->create(fmi,3,"forcespin:fmi");
}
/* ---------------------------------------------------------------------- */
@ -160,7 +170,7 @@ void FixForceSpin::post_force(int vflag)
set_magneticforce(); //Update value of the mag. field if time-dependent
}
double **x = atom->x;
// double **x = atom->x;
double **sp = atom->sp;
double *mumag = atom->mumag;
double **fm = atom->fm;
@ -169,22 +179,42 @@ void FixForceSpin::post_force(int vflag)
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;
}
}
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] += scalar*xmag;
fm[i][1] += scalar*ymag;
fm[i][2] += scalar*zmag;
}
}
for (int i = 0; i < nlocal; i++) {
fmi[0] = fmi[1] = fmi[2] = 0.0;
//if (style == ZEEMAN) {
if (zeeman_flag) {
compute_zeeman(i,fmi);
}
//if (style == ANISOTROPY) {
if (aniso_flag) {
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
compute_anisotropy(i,spi,fmi);
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
}
}
/* ---------------------------------------------------------------------- */
void FixForceSpin::compute_zeeman(int i, double *fmi)
{
double *mumag = atom->mumag;
fmi[0] += mumag[i]*xmag;
fmi[1] += mumag[i]*ymag;
fmi[2] += mumag[i]*zmag;
}
void FixForceSpin::compute_anisotropy(int i, double * spi, double *fmi)
{
double scalar = Kax*spi[0] + Kay*spi[1] + Kaz*spi[2];
fmi[0] += scalar*xmag;
fmi[1] += scalar*ymag;
fmi[2] += scalar*zmag;
}
/* ---------------------------------------------------------------------- */

8
src/SPIN/fix_force_spin.h
View File

@ -37,6 +37,10 @@ class FixForceSpin : public Fix {
virtual void post_force_respa(int, int, int);
double compute_scalar();
int zeeman_flag, aniso_flag;
void compute_zeeman(int, double *);
void compute_anisotropy(int, double *, double *);
protected:
int style;
@ -58,8 +62,10 @@ class FixForceSpin : public Fix {
double Ka; //Magnetic anisotropy intensity and direction
double Kax, Kay, Kaz;
double *spi, *fmi; //Temp var. in compute
void set_magneticforce();
};
}

32
src/SPIN/fix_langevin_spin.cpp
View File

@ -145,25 +145,25 @@ void FixLangevinSpin::post_force(int vflag)
// add the damping to the effective field of each spin
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
sx = sp[i][0];
sy = sp[i][1];
sz = sp[i][2];
fmx = fm[i][0];
fmy = fm[i][1];
fmz = fm[i][2];
sx = sp[i][0];
sy = sp[i][1];
sz = sp[i][2];
fmx = fm[i][0];
fmy = fm[i][1];
fmz = fm[i][2];
cpx = fmy*sz - fmz*sy;//Computing cross product
cpy = fmz*sx - fmx*sz;
cpz = fmx*sy - fmy*sx;
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;
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;
fm[i][0] = fmx;
fm[i][1] = fmy;
fm[i][2] = fmz;
}
//apply thermal effects adding random fields to fm

369
src/SPIN/fix_nve_spin.cpp
View File

@ -26,8 +26,15 @@
#include "math_const.h"
#include "modify.h"
//Add headers (see delete later)
#include "pair.h"
#include "timer.h"
#include "integrate.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "pair_spin.h"
#include "memory.h"
#include "fix_force_spin.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -57,7 +64,24 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
// error checks
if (extra == SPIN && !atom->mumag_flag)
error->all(FLERR,"Fix nve/spin requires spin attribute mumag");
#if defined SEQNEI
lockpairspin = NULL;
lockforcespin = NULL;
exch_flag = dmi_flag = me_flag = 0;
zeeman_flag = aniso_flag = 0;
#endif
}
/* ---------------------------------------------------------------------- */
FixNVESpin::~FixNVESpin(){
#if defined SEQNEI
//delete lockpairspin;
//delete lockforcespin;
memory->destroy(spi);
memory->destroy(fmi);
memory->destroy(fmj);
#endif
}
/* ---------------------------------------------------------------------- */
@ -68,6 +92,26 @@ void FixNVESpin::init()
dts = update->dt;
#if defined SEQNEI
lockpairspin = (PairSpin *) force->pair;
memory->create(spi,3,"nves:spi");
memory->create(fmi,3,"nves:fmi");
memory->create(fmj,3,"nves:fmj");
int iforce;
for (iforce = 0; iforce < modify->nfix; iforce++)
if (strstr(modify->fix[iforce]->style,"force/spin")) break;
lockforcespin = (FixForceSpin *) modify->fix[iforce];
exch_flag = lockpairspin->exch_flag;
dmi_flag = lockpairspin->dmi_flag;
me_flag = lockpairspin->me_flag;
zeeman_flag = lockforcespin->zeeman_flag;
aniso_flag = lockforcespin->aniso_flag;
#endif
/*int idamp;
for (idamp = 0; idamp < modify->nfix; idamp++)
if (strstr(modify->fix[idamp]->style,"damping/spin")) break;
@ -101,9 +145,18 @@ void FixNVESpin::initial_integrate(int vflag)
// Advance half spins all particles
//See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011
if (extra == SPIN) {
for (int i = 0; i < nlocal; i++){
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#if defined SEQNEI
for (int i = 0; i < nlocal; i++){
ComputeSpinInteractionsNei(i);
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#endif
#if defined SEQ
for (int i = 0; i < nlocal; i++){
AdvanceSingleSpin(i,0.5*dts,sp,fm);
ComputeSpinInteractions();
}
#endif
}
// update half v all particles
@ -125,18 +178,293 @@ void FixNVESpin::initial_integrate(int vflag)
x[i][2] += 0.5 * dtv * v[i][2];
}
}
}
//#define FORCE_PRINT
#if defined FORCE_PRINT
FILE* file_force=NULL;
file_force=fopen("spin_force_Lammps.dat","a");
fprintf(file_force,"---------------------------------- \n");
for(int i=0;i<nlocal;i++){
fprintf(file_force,"%d %lf %lf %lf \n",i,fm[i][0],fm[i][1],fm[i][2]);
}
if (file_force!=NULL) fclose(file_force);
#if defined SEQNEI
/* ---------------------------------------------------------------------- */
void FixNVESpin::ComputeSpinInteractionsNei(int ii)
{
int nflag,sortflag;
int nlocal = atom->nlocal;
int n_post_integrate = modify->n_post_integrate;
int n_pre_exchange = modify->n_pre_exchange;
int n_pre_neighbor = modify->n_pre_neighbor;
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
int n_end_of_step = modify->n_end_of_step;
bigint ntimestep;
ntimestep = update->ntimestep;
//Force compute quantities
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double **x = atom->x;
double **sp = atom->sp;
double **fm = atom->fm;
int *type = atom->type;
int newton_pair = force->newton_pair;
inum = lockpairspin->list->inum;
ilist = lockpairspin->list->ilist;
numneigh = lockpairspin->list->numneigh;
firstneigh = lockpairspin->list->firstneigh;
double xtmp,ytmp,ztmp;
double rsq,rd,delx,dely,delz;
double cut_ex_2, cut_dmi_2, cut_me_2;
cut_ex_2 = cut_dmi_2 = cut_me_2 = 0.0;
int eflag = 1;
int vflag = 0;
int pair_compute_flag = 1;
if (atom->sortfreq > 0) sortflag = 1;
else sortflag = 0;
if (n_post_integrate) modify->post_integrate();
timer->stamp(Timer::MODIFY);
// regular communication vs neighbor list rebuild
nflag = neighbor->decide();
if (nflag == 0) {
timer->stamp();
comm->forward_comm();
timer->stamp(Timer::COMM);
} else {
if (n_pre_exchange) {
timer->stamp();
modify->pre_exchange();
timer->stamp(Timer::MODIFY);
}
//if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
if (sortflag && ntimestep >= atom->nextsort) atom->sort();
comm->borders();
//if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
if (n_pre_neighbor) {
modify->pre_neighbor();
timer->stamp(Timer::MODIFY);
}
neighbor->build();
timer->stamp(Timer::NEIGH);
}
///////Force computation for spin i/////////////
i = ilist[ii];
//Clear atom i
fm[i][0] = fm[i][1] = fm[i][2] = 0.0;
timer->stamp();
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
fmi[0] = fmi[1] = fmi[2] = 0.0;
fmj[0] = fmj[1] = fmj[2] = 0.0;
jlist = firstneigh[i];
jnum = numneigh[i];
// printf("Test inum: %g \n",inum);
/*
//Pair interaction
for (int jj = 0; jj < inum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
itype = type[ii];
jtype = type[j];
if (exch_flag) {
cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]);
if (rsq <= cut_ex_2) {
lockpairspin->compute_exchange(i,j,rsq,fmi,fmj);
}
}
if (dmi_flag) {
cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]);
if (rsq <= cut_dmi_2) {
lockpairspin->compute_dmi(i,j,fmi,fmj);
}
}
if (me_flag) {
cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]);
if (rsq <= cut_me_2) {
lockpairspin->compute_me(i,j,fmi,fmj);
}
}
}
*/
//Pair interaction
int natom = nlocal + atom->nghost;
for (int k = 0; k < natom; k++) {
delx = xtmp - x[k][0];
dely = ytmp - x[k][1];
delz = ztmp - x[k][2];
rsq = delx*delx + dely*dely + delz*delz;
itype = type[ii];
jtype = type[k];
if (exch_flag) {
cut_ex_2 = (lockpairspin->cut_spin_exchange[itype][jtype])*(lockpairspin->cut_spin_exchange[itype][jtype]);
if (rsq <= cut_ex_2) {
lockpairspin->compute_exchange(i,k,rsq,fmi,fmj);
}
}
if (dmi_flag) {
cut_dmi_2 = (lockpairspin->cut_spin_dmi[itype][jtype])*(lockpairspin->cut_spin_dmi[itype][jtype]);
if (rsq <= cut_dmi_2) {
lockpairspin->compute_dmi(i,k,fmi,fmj);
}
}
if (me_flag) {
cut_me_2 = (lockpairspin->cut_spin_me[itype][jtype])*(lockpairspin->cut_spin_me[itype][jtype]);
if (rsq <= cut_me_2) {
lockpairspin->compute_me(i,k,fmi,fmj);
}
}
}
//post force
if (zeeman_flag) {
lockforcespin->compute_zeeman(i,fmi);
}
if (aniso_flag) {
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
lockforcespin->compute_anisotropy(i,spi,fmi);
}
//Replace the force by its new value
fm[i][0] = fmi[0];
fm[i][1] = fmi[1];
fm[i][2] = fmi[2];
}
#endif
/* ---------------------------------------------------------------------- */
void FixNVESpin::ComputeSpinInteractions()
{
int nflag,sortflag;
int nlocal = atom->nlocal;
int n_post_integrate = modify->n_post_integrate;
int n_pre_exchange = modify->n_pre_exchange;
int n_pre_neighbor = modify->n_pre_neighbor;
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
int n_end_of_step = modify->n_end_of_step;
bigint ntimestep;
ntimestep = update->ntimestep;
//int eflag = update->integrate->eflag;
//int vflag = update->integrate->vflag;
int eflag = 1;
int vflag = 0;
int pair_compute_flag = 1;
if (atom->sortfreq > 0) sortflag = 1;
else sortflag = 0;
if (n_post_integrate) modify->post_integrate();
timer->stamp(Timer::MODIFY);
// regular communication vs neighbor list rebuild
nflag = neighbor->decide();
if (nflag == 0) {
timer->stamp();
comm->forward_comm();
timer->stamp(Timer::COMM);
} else {
if (n_pre_exchange) {
timer->stamp();
modify->pre_exchange();
timer->stamp(Timer::MODIFY);
}
//if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
if (sortflag && ntimestep >= atom->nextsort) atom->sort();
comm->borders();
//if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
if (n_pre_neighbor) {
modify->pre_neighbor();
timer->stamp(Timer::MODIFY);
}
neighbor->build();
timer->stamp(Timer::NEIGH);
}
// force computations
// important for pair to come before bonded contributions
// since some bonded potentials tally pairwise energy/virial
// and Pair:ev_tally() needs to be called before any tallying
size_t nbytes;
nbytes = sizeof(double) * nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;
atom->avec->force_clear(0,nbytes);
timer->stamp();
if (n_pre_force) {
modify->pre_force(vflag);
timer->stamp(Timer::MODIFY);
}
if (pair_compute_flag) {
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
}
/*if (kspace_compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(Timer::KSPACE);
}*/
if (n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
// reverse communication of forces
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
}
/* ---------------------------------------------------------------------- */
@ -218,8 +546,17 @@ void FixNVESpin::final_integrate()
// Advance half spins all particles
//See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011
if (extra == SPIN) {
for (int i = 0; i < nlocal; i++){
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#if defined SEQNEI
for (int i = nlocal-1; i >= 0; i--){
ComputeSpinInteractionsNei(i);
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#endif
#if defined SEQ
for (int i = nlocal-1; i >= 0; i--){
AdvanceSingleSpin(i,0.5*dts,sp,fm);
ComputeSpinInteractions();
}
#endif
}
}

20
src/SPIN/fix_nve_spin.h
View File

@ -25,24 +25,34 @@ FixStyle(nve/spin,FixNVESpin)
namespace LAMMPS_NS {
class FixNVESpin : public FixNVE {
friend class FixSpinDamping;
public:
FixNVESpin(class LAMMPS *, int, char **);
virtual ~FixNVESpin() {}
virtual ~FixNVESpin();
void init();
virtual void initial_integrate(int);
void AdvanceSingleSpin(int, double, double **, double **);
virtual void final_integrate();
//Sorting atoms/spins routine
void SortSpins();
//#define SEQ
#define SEQNEI
void ComputeSpinInteractions();
void ComputeSpinInteractionsNei(int);
protected:
int extra;
double dts;
//double alpha_t;
#if defined SEQNEI
private:
int exch_flag, dmi_flag, me_flag;
int zeeman_flag, aniso_flag;
class PairSpin *lockpairspin;
double *spi, *fmi, *fmj; //Temp var. for compute
class FixForceSpin *lockforcespin;
#endif
//private:
//class FixSpinDamping *lockspindamping;
};

41
src/SPIN/pair_spin.cpp
View File

@ -25,6 +25,9 @@
#include "update.h"
#include <string.h>
//Add. lib. for full neighb. list
#include "neigh_request.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -34,7 +37,9 @@ PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
exch_flag = 0;
dmi_flag = 0;
dmi_flag = 0;
me_flag = 0;
}
/* ---------------------------------------------------------------------- */
@ -76,7 +81,7 @@ void PairSpin::compute(int eflag, int vflag)
double evdwl,ecoul;
double xtmp,ytmp,ztmp;
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
double cut_ex,cut_dmi,cut_me;
double cut_ex_2,cut_dmi_2,cut_me_2;
double rsq,rd,delx,dely,delz;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -120,28 +125,27 @@ void PairSpin::compute(int eflag, int vflag)
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
itype = type[i];
jtype = type[j];
//Exchange interaction
if (exch_flag) {
cut_ex = cut_spin_exchange[itype][jtype];
if (rd <= cut_ex) {
cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
if (rsq <= cut_ex_2) {
compute_exchange(i,j,rsq,fmi,fmj);
}
}
//DM interaction
if (dmi_flag){
cut_dmi = cut_spin_dmi[itype][jtype];
if (rd <= cut_dmi){
cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
if (rsq <= cut_dmi_2){
compute_dmi(i,j,fmi,fmj);
}
}
//ME interaction
if (me_flag){
cut_me = cut_spin_me[itype][jtype];
if (rd <= cut_me){
cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
if (rsq <= cut_me_2){
compute_me(i,j,fmi,fmj);
}
}
@ -162,7 +166,7 @@ void PairSpin::compute(int eflag, int vflag)
}
/* ---------------------------------------------------------------------- */
inline void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj)
void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj)
{
int *type = atom->type;
int itype, jtype;
@ -184,13 +188,12 @@ inline void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, d
fmj[0] += Jex*sp[i][0];
fmj[1] += Jex*sp[i][1];
fmj[2] += Jex*sp[i][2];
}
/* ---------------------------------------------------------------------- */
inline void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj)
void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj)
{
int *type = atom->type;
int itype, jtype;
double **sp = atom->sp;
@ -212,7 +215,7 @@ inline void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj)
}
/* ---------------------------------------------------------------------- */
inline void PairSpin::compute_me(int i, int j, double *fmi, double *fmj)
void PairSpin::compute_me(int i, int j, double *fmi, double *fmj)
{
int *type = atom->type;
int itype, jtype;
@ -247,8 +250,6 @@ inline void PairSpin::compute_me(int i, int j, double *fmi, double *fmj)
fmj[1] -= sp[i][2]*meix - sp[i][0]*meiz;
fmj[2] -= sp[i][0]*meiy - sp[i][1]*meix;
// printf("test val fmi=%g, fmj=%g \n",fmi[2],fmj[2]);
}
/* ----------------------------------------------------------------------
@ -312,6 +313,7 @@ void PairSpin::settings(int narg, char **arg)
if (setflag[i][j]) {
cut_spin_exchange[i][j] = cut_spin_pair_global;
cut_spin_dmi[i][j] = cut_spin_pair_global;
cut_spin_me[i][j] = cut_spin_pair_global;
}
}
@ -438,6 +440,13 @@ void PairSpin::init_style()
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
neighbor->request(this,instance_me);
#define FULLNEI
#if defined FULLNEI
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
#endif
}
/* ----------------------------------------------------------------------

17
src/SPIN/pair_spin.h
View File

@ -39,19 +39,22 @@ class PairSpin : public Pair {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
inline void compute_exchange(int, int, double, double *, double *);
inline void compute_dmi(int, int, double *, double *);
inline void compute_me(int, int, double *, double *);
void compute_exchange(int, int, double, double *, double *);
void compute_dmi(int, int, double *, double *);
void compute_me(int, int, double *, double *);
protected:
//Test for seq. integ.
//protected:
int exch_flag,dmi_flag,me_flag;
double cut_spin_pair_global;
double cut_spin_dipolar_global;
double **cut_spin_exchange; //cutting distance exchange
double **cut_spin_dmi; //cutting distance dmi
double **cut_spin_me; //cutting distance me
double **cut_spin_dmi; //cutting distance dmi
double **cut_spin_me; //cutting distance me
//Test for seq. integ.
protected:
double **J_1, **J_2, **J_3; //exchange coeffs Jij
//J1 in eV, J2 adim and J3 in Ang
double **DM;

1
src/modify.cpp
View File

@ -1,4 +1,5 @@
/* ----------------------------------------------------------------------
eoundary p f f
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov