ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added:

Browse Source 
- For Paramag. simulations, the option "atom_modify" has to be set
ex: atom_modify sort 1000 4.0 (Freq,Dist).
- Actual time is now printed (c_mag[0] in compute_spin)
- Value of Gilbert's damping corrected
- Now even results for SD/Lammps comp. in purely paramg. or aniso. situations
- Pack and unpack reverse needed corrections (f only was set, not fm)
- Spin temperature is now computed (data c_mag[7] in spin_compute)

To do:
- Fcc with p p p bc is still not working
- If Zeeman/Aniso force not defined, error => to be removed
- Add DMI and ME (see if new file or add in the exchange file)
This commit is contained in:
julient31
2017-05-19 13:24:40 -06:00
parent 3168704858
commit af45d55b3f
9 changed files with 154 additions and 82 deletions
Download Patch File Download Diff File Expand all files Collapse all files

32
in.cobalt
View File

@ -10,14 +10,13 @@ units metal
dimension 3
#setting boundary conditions. (p for periodic, f for fixed, ...)
#boundary p p p
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
#atom_style hybrid sphere spin
atom_modify sort 1000 4.0
#Define sort for paramagnetic simulations (if no pair interaction)
#atom_modify sort 1000 4.0
###########################
#######Create atoms########
@ -25,10 +24,11 @@ atom_modify sort 1000 4.0
#Lattice constant of fcc Cobalt
lattice fcc 3.54
#lattice sc 3.54
#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 1.0 0.0 1.0 0.0 1.0
region box block 0.0 3.0 0.0 3.0 0.0 3.0
#Creating a simulation box based on the specified region. Entries: number of atom types and box ref.
create_box 1 box
@ -50,21 +50,21 @@ 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.7
set group all spin 1.7 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 4.0
#pair_coeff * * 0.0446928 0.003496 1.4885
pair_style pair/spin 4.0
pair_coeff * * 0.0446928 0.003496 1.4885
#pair_coeff * * 0.0 0.003496 1.4885
#Fix Langevin spins (merging damping and temperature)
#Defines a cutof distance for the neighbors.
#neighbor 0.5 bin
#Magnetic field fix
fix 1 all force/spin zeeman 10 0.0 0.0 1.0
fix 1 all force/spin zeeman 1.0 0.0 0.0 1.0
#fix 1 all force/spin zeeman 0.0 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.001 0.0 0.0 1.0
#Fix Langevin spins (merging damping and temperature)
fix 2 all langevin/spin 0.0 0.1 0.0 21
@ -75,22 +75,22 @@ fix 3 all nve/spin
#compute total magnetization and magnetic energy
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] 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.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.000005
##################
#######run########
##################
#Dump the positions and spin directions of magnetic particles (vmd format)
dump 1 all custom 100 dump_spin.lammpstrj type x y z spx spy spz
dump 1 all custom 50 dump_spin.lammpstrj type x y z spx spy spz
#Running the simulations for N timesteps
run 1000000
#run 100
run 200000
#run 2

34
src/SPIN/atom_vec_spin.cpp
View File

@ -44,7 +44,6 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
forceclearflag = 1;
atom->mumag_flag = atom->sp_flag = 1;
//atom->rmass_flag = 1;
}
@ -138,6 +137,7 @@ int AtomVecSpin::pack_comm(int n, int *list, double *buf,
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) {
@ -157,6 +157,7 @@ int AtomVecSpin::pack_comm(int n, int *list, double *buf,
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
}
return m;
@ -180,6 +181,7 @@ int AtomVecSpin::pack_comm_vel(int n, int *list, double *buf,
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];
@ -203,6 +205,7 @@ int AtomVecSpin::pack_comm_vel(int n, int *list, double *buf,
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];
@ -219,6 +222,7 @@ int AtomVecSpin::pack_comm_vel(int n, int *list, double *buf,
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;
@ -246,6 +250,7 @@ int AtomVecSpin::pack_comm_hybrid(int n, int *list, double *buf)
buf[m++] = sp[j][0];
buf[m++] = sp[j][1];
buf[m++] = sp[j][2];
buf[m++] = sp[j][3];
}
return m;
}
@ -265,6 +270,7 @@ void AtomVecSpin::unpack_comm(int n, int first, double *buf)
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
}
}
@ -283,6 +289,7 @@ void AtomVecSpin::unpack_comm_vel(int n, int first, double *buf)
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++];
@ -301,6 +308,7 @@ int AtomVecSpin::unpack_comm_hybrid(int n, int first, double *buf)
sp[i][0] = buf[m++];
sp[i][1] = buf[m++];
sp[i][2] = buf[m++];
sp[i][3] = buf[m++];
}
return m;
}
@ -317,6 +325,9 @@ int AtomVecSpin::pack_reverse(int n, int first, double *buf)
buf[m++] = f[i][0];
buf[m++] = f[i][1];
buf[m++] = f[i][2];
buf[m++] = fm[i][0];
buf[m++] = fm[i][1];
buf[m++] = fm[i][2];
}
return m;
}
@ -333,6 +344,9 @@ void AtomVecSpin::unpack_reverse(int n, int *list, double *buf)
f[j][0] += buf[m++];
f[j][1] += buf[m++];
f[j][2] += buf[m++];
fm[j][0] += buf[m++];
fm[j][1] += buf[m++];
fm[j][2] += buf[m++];
}
}
@ -855,9 +869,10 @@ void AtomVecSpin::pack_data(double **buf)
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;
buf[i][9] = sp[i][3];
buf[i][10] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
buf[i][11] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
buf[i][12] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
}
}
@ -871,8 +886,9 @@ int AtomVecSpin::pack_data_hybrid(int i, double *buf)
buf[1] = sp[i][0];
buf[2] = sp[i][1];
buf[3] = sp[i][2];
buf[4] = sp[i][3];
return 4;
return 5;
}
/* ----------------------------------------------------------------------
@ -887,9 +903,9 @@ void AtomVecSpin::write_data(FILE *fp, int n, double **buf)
"%-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);
buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9],
(int) ubuf(buf[i][10]).i,(int) ubuf(buf[i][11]).i,
(int) ubuf(buf[i][12]).i);
}
/* ----------------------------------------------------------------------
@ -898,7 +914,7 @@ void AtomVecSpin::write_data(FILE *fp, int n, double **buf)
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]);
fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2],buf[3],buf[4]);
return 4;
}

23
src/SPIN/compute_spin.cpp
View File

@ -20,9 +20,12 @@
#include "memory.h"
#include "error.h"
#include "math_special.h"
#include "math_const.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace MathSpecial;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
@ -32,7 +35,7 @@ ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command");
vector_flag = 1;
size_vector = 6;
size_vector = 7;
extvector = 0;
init();
@ -52,6 +55,8 @@ ComputeSpin::~ComputeSpin()
void ComputeSpin::init()
{
hbar = force->hplanck/MY_2PI;
kb = force->boltz;
}
/* ---------------------------------------------------------------------- */
@ -66,6 +71,9 @@ void ComputeSpin::compute_vector()
mag[0] = mag[1] = mag[2] = mag[3] = 0.0;
magtot[0] = magtot[1] = magtot[2] = magtot[3] = 0.0;
magenergy = magenergytot = 0.0;
tempnum = tempnumtot = 0.0;
tempdenom = tempdenomtot = 0.0;
spintemperature = 0.0;
double **x = atom->x;
int *mask = atom->mask;
@ -74,6 +82,7 @@ void ComputeSpin::compute_vector()
double *mumag = atom->mumag;
double **sp = atom->sp;
double **fm = atom->fm;
double tx,ty,tz;
int nlocal = atom->nlocal;
@ -87,6 +96,11 @@ void ComputeSpin::compute_vector()
magenergy += mumag[i]*sp[i][0]*fm[i][0];
magenergy += mumag[i]*sp[i][1]*fm[i][1];
magenergy += mumag[i]*sp[i][2]*fm[i][2];
tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0];
tempnum += tx*tx+ty*ty+tz*tz;
tempdenom += sp[i][0]*sp[i][0]+sp[i][1]*sp[i][1]+sp[i][2]*sp[i][2];
countsp++;
}
}
@ -95,6 +109,8 @@ void ComputeSpin::compute_vector()
MPI_Allreduce(mag,magtot,4,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&magenergy,&magenergytot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&tempnum,&tempnumtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&tempdenom,&tempdenomtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&countsp,&countsptot,1,MPI_DOUBLE,MPI_SUM,world);
double scale = 1.0/countsptot;
@ -102,6 +118,7 @@ void ComputeSpin::compute_vector()
magtot[1] *= scale;
magtot[2] *= scale;
magtot[3] = sqrt(square(magtot[0])+square(magtot[1])+square(magtot[2]));
spintemperature = hbar*tempnumtot/2.0/kb/tempdenomtot;
vector[0] = invoked_vector*update->dt;
vector[1] = magtot[0];
@ -109,6 +126,8 @@ void ComputeSpin::compute_vector()
vector[3] = magtot[2];
vector[4] = magtot[3];
vector[5] = magenergytot;
vector[6] = spintemperature;
}
/* ----------------------------------------------------------------------
@ -120,6 +139,6 @@ void ComputeSpin::allocate()
memory->destroy(mag);
memory->create(mag,4,"compute/spin:mag");
memory->create(magtot,5,"compute/spin:mag");
memory->create(vector,6,"compute/spin:vector");
memory->create(vector,7,"compute/spin:vector");
}

4
src/SPIN/compute_spin.h
View File

@ -36,6 +36,10 @@ class ComputeSpin : public Compute {
double *magtot;
double magenergy;
double magenergytot;
double tempnum,tempnumtot;
double tempdenom,tempdenomtot;
double spintemperature;
double kb,hbar;
int countsp;
int countsptot;
int usecenter;

18
src/SPIN/fix_force_spin.cpp
View File

@ -77,8 +77,6 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
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;
@ -177,19 +175,17 @@ void FixForceSpin::post_force(int vflag)
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]);
fm[i][0] += scalar*xmag;
fm[i][1] += scalar*ymag;
fm[i][2] += scalar*zmag;
}
}
printf("test force. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//printf("test force. 1;i=0, fx=%g, fy=%g, fz=%g, mumag=%g \n",fm[0][0],fm[0][1],fm[0][2],mumag[0]);
//printf("Field force compute, fm[0][2]=%g \n",fm[0][2]);
}
@ -211,6 +207,12 @@ void FixForceSpin::set_magneticforce()
ymag = H_field*Hy;
zmag = H_field*Hz;
}
if (style == ANISOTROPY) {
xmag = 2.0*Ka*Kax;
ymag = 2.0*Ka*Kay;
zmag = 2.0*Ka*Kaz;
}
}

5
src/SPIN/fix_langevin_spin.cpp
View File

@ -107,7 +107,7 @@ void FixLangevinSpin::init()
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));
Gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t));
double hbar = force->hplanck/MY_2PI; //eV/(rad.THz)
double kb = force->boltz;
@ -189,6 +189,9 @@ void FixLangevinSpin::post_force(int vflag)
fm[i][2] *= Gil_factor;
}
//printf("test langevin 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//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]);

17
src/SPIN/fix_nve_spin.cpp
View File

@ -15,7 +15,6 @@
#include <stdio.h>
#include <string.h>
#include "fix_nve_spin.h"
//#include "fix_damping_spin.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
@ -121,9 +120,7 @@ void FixNVESpin::initial_integrate(int vflag)
// 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) {
if (mask[i] & groupbit) {
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]);
@ -134,10 +131,6 @@ void FixNVESpin::initial_integrate(int vflag)
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;
@ -146,9 +139,9 @@ void FixNVESpin::initial_integrate(int vflag)
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));
g[0] /= (1+0.25*fm2*dts*dts);
g[1] /= (1+0.25*fm2*dts*dts);
g[2] /= (1+0.25*fm2*dts*dts);
sp[i][0] = g[0];
sp[i][1] = g[1];
@ -161,7 +154,7 @@ void FixNVESpin::initial_integrate(int vflag)
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. 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);
}

77
src/SPIN/pair_spin.cpp
View File

@ -41,14 +41,13 @@ 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);
}
}
@ -56,23 +55,24 @@ PairSpin::~PairSpin()
void PairSpin::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl,ecoul;
double xtmp,ytmp,ztmp,fmix,fmiy,fmiz,fmjx,fmjy,fmjz,omx,omy,omz;
double cut, Jex, ra;
double rsq,rd,delx,dely,delz;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **fm = atom->fm;
double *mumag = atom->mumag;
double **sp = atom->sp;
int *type = atom->type;
int 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;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
@ -81,6 +81,7 @@ void PairSpin::compute(int eflag, int vflag)
// Pair spin computations
// Loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
@ -89,9 +90,20 @@ void PairSpin::compute(int eflag, int vflag)
jlist = firstneigh[i];
jnum = numneigh[i];
//printf(":::::::Loop for atom numb. %d, jnum=%d ::::::: \n",i,jnum);
//printf("Test print real atom: i=%d, sx=%g, sy=%g, sz=%g \n",i,sp[i][0],sp[i][1],sp[i][2]);
//printf("Test print real atom: i=%d, rx=%g, ry=%g, rz=%g \n",i,x[i][0],x[i][1],x[i][2]);
int testcount=0;
//Exchange interaction
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
fmix = fmiy = fmiz = 0.0;
fmjx = fmjy = fmjz = 0.0;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
@ -103,8 +115,6 @@ void PairSpin::compute(int eflag, int vflag)
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;
@ -116,22 +126,47 @@ void PairSpin::compute(int eflag, int vflag)
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];
fmjx = Jex*sp[i][0];
fmjy = Jex*sp[i][1];
fmjz = Jex*sp[i][2];
//printf("Neighb pair: i=%d, j=%d \n",i,j);
//printf("Test print ghost/real neib atom: i=%d, j=%d, sx=%g, sy=%g, sz=%g \n",i,j,sp[j][0],sp[j][1],sp[j][2]);
//printf("Test g/r neib pair: i=%d, j=%d, rx=%g, ry=%g, rz=%g \n",i,j,x[j][0],x[j][1],x[j][2]);
//printf("Atom i: %d of type %d, Atom j: %d of type %d, \n",i,itype,j,jtype);
//printf("Exchange pair (%d,%d), Jij=%g, rij=%g \n",i,j,Jex,rd);
testcount++;
}
//printf("Test print ghost/real atom: j=%d, sx=%g, sy=%g, sz=%g \n",j,sp[j][0],sp[j][1],sp[j][2]);
//printf("Test print ghost/real atom: j=%d, rx=%g, ry=%g, rz=%g \n",j,x[j][0],x[j][1],x[j][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;
if (newton_pair || j < nlocal) {
fm[j][0] += fmjx;
fm[j][1] += fmjy;
fm[j][2] += fmjz;
}
//printf("Val fm %d: [%g,%g,%g] \n",i,fm[i][0],fm[i][1],fm[i][2]);
//printf("Val fm %d: [%g,%g,%g] \n",j,fm[j][0],fm[j][1],fm[j][2]);
}
// printf("Test count %d \n",testcount);
}
//printf("New pair val: %d \n",newton_pair);
//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]);
//printf("::::::::::::::::::::::::: End loop ::::::::::::::::::::::::\n");
if (vflag_fdotr) virial_fdotr_compute();
}

2
src/set.cpp
View File

@ -860,7 +860,7 @@ void Set::set(int keyword)
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;
mumag[i] = dvalue;
}
// set quaternion orientation of ellipsoid or tri or body particle