Memory corrections
This commit is contained in:
julient31
109
examples/SPIN/in.spin.bfo
Normal file
109
examples/SPIN/in.spin.bfo
Normal file
@ -0,0 +1,109 @@
|
||||
###################
|
||||
#######Init########
|
||||
###################
|
||||
|
||||
clear
|
||||
#setting units to metal (Ang, picosecs, eV, ...):
|
||||
units metal
|
||||
|
||||
#setting dimension of the system (N=2 or 3):
|
||||
dimension 3
|
||||
|
||||
#setting boundary conditions. (p for periodic, f for fixed, ...)
|
||||
boundary p p f
|
||||
|
||||
#setting atom_style to spin:
|
||||
atom_style spin
|
||||
|
||||
#Define sort for paramagnetic simulations (if no pair interaction)
|
||||
#atom_modify sort 1000 4.0
|
||||
|
||||
#why?
|
||||
atom_modify map array
|
||||
|
||||
#newton off for pair spin in SEQNEI
|
||||
#newton off off
|
||||
|
||||
###########################
|
||||
#######Create atoms########
|
||||
###########################
|
||||
|
||||
#Lattice constant of sc Iron atoms of BFO
|
||||
lattice sc 3.96
|
||||
|
||||
#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 17.0 0.0 17.0 0.0 5.0
|
||||
|
||||
#Creating a simulation box based on the specified region. Entries: number of atom types and box ref.
|
||||
create_box 1 box
|
||||
|
||||
#Creating atoms (or molecules) on a lattice, or a single atom (or molecule), ...
|
||||
#Entries: atom type,
|
||||
create_atoms 1 box
|
||||
|
||||
#Replicating NxNxN the entire set of atoms
|
||||
#replicate 1 1 1
|
||||
|
||||
|
||||
#######################
|
||||
#######Settings########
|
||||
#######################
|
||||
|
||||
#Setting one or more properties of one or more atoms.
|
||||
#Setting mass
|
||||
mass 1 1.0
|
||||
#set group all mass 1.0
|
||||
#Setting spins orientation and moment
|
||||
set group all spin/random 11 2.50
|
||||
#set group all spin 2.50 1.0 0.0 0.0
|
||||
|
||||
#Magnetic exchange interaction coefficient for bulk fcc Cobalt
|
||||
pair_style pair/spin 5.7
|
||||
#type i and j | interaction type | cutoff | J1 (eV) | J2 (adim) | J3 (Ang) (for Exchange)
|
||||
pair_coeff * * exchange 5.7 -0.01575 0.0 1.965
|
||||
#type i and j | interaction type | cutoff | Int (eV) | [dx,dy,dz] (for DMI and ME)
|
||||
#pair_coeff * * dmi 2.6 0.001 0.0 0.0 1.0
|
||||
#pair_coeff * * me 4.0 0.000109 1.0 1.0 1.0
|
||||
#Mex10
|
||||
pair_coeff * * me 4.0 0.00109 1.0 1.0 1.0
|
||||
|
||||
#Define a skin distance, update neigh list every
|
||||
#neighbor 1.0 bin
|
||||
#neigh_modify every 10 check yes delay 20
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 check no delay 0
|
||||
|
||||
#Magnetic field fix
|
||||
#Type | Intensity (T or eV) | Direction
|
||||
#fix 1 all force/spin zeeman 1.0 0.0 0.0 1.0
|
||||
fix 1 all force/spin anisotropy 0.0000035 0.0 0.0 1.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.1 0.0 21
|
||||
#fix 2 all langevin/spin 0.0 0.0 0.0 21
|
||||
|
||||
#Magnetic integration fix
|
||||
fix 3 all nve/spin mpi
|
||||
|
||||
#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 500 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag_BFO.dat format %20.16g
|
||||
|
||||
#Setting the timestep for the simulation (in ps)
|
||||
timestep 0.0001
|
||||
|
||||
##################
|
||||
#######run########
|
||||
##################
|
||||
|
||||
#Dump the positions and spin directions of magnetic particles (vmd format)
|
||||
dump 1 all custom 5000 dump_spin_BFO.lammpstrj type x y z spx spy spz
|
||||
|
||||
#Running the simulations for N timesteps
|
||||
run 30000
|
||||
#run 1
|
||||
|
||||
120
examples/SPIN/in.spin.cobalt
Normal file
120
examples/SPIN/in.spin.cobalt
Normal file
@ -0,0 +1,120 @@
|
||||
###################
|
||||
#######Init########
|
||||
###################
|
||||
|
||||
clear
|
||||
units metal
|
||||
dimension 3
|
||||
boundary p p p
|
||||
#boundary f f f
|
||||
|
||||
#newton off
|
||||
|
||||
#setting atom_style to spin:
|
||||
atom_style spin
|
||||
|
||||
#Define sort for paramagnetic simulations (if no pair interaction)
|
||||
#atom_modify sort 1000 4.0
|
||||
|
||||
atom_modify map array
|
||||
|
||||
###########################
|
||||
#######Create atoms########
|
||||
###########################
|
||||
|
||||
lattice fcc 3.54
|
||||
region box block 0.0 5.0 0.0 5.0 0.0 5.0
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
|
||||
#######################
|
||||
#######Settings########
|
||||
#######################
|
||||
|
||||
#isolating 1 atom into a group
|
||||
group single_spin id 10
|
||||
|
||||
#Setting one or more properties of one or more atoms.
|
||||
mass 1 58.93
|
||||
|
||||
#Setting spins orientation and moment
|
||||
set group all spin/random 31 1.72
|
||||
#set group all spin 1.72 0.0 0.0 1.0
|
||||
#set group single_spin spin/random 11 1.72
|
||||
|
||||
velocity all create 200 4928459 rot yes dist gaussian
|
||||
|
||||
#Magneto-mechanic interactions for bulk fcc Cobalt
|
||||
#pair_style pair/spin/exchange 4.0
|
||||
#pair_style eam/alloy
|
||||
#pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0
|
||||
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
|
||||
|
||||
pair_coeff * * eam/alloy ../examples/SPIN/Co_PurjaPun_2012.eam.alloy Co
|
||||
#pair_coeff * * ../Co_PurjaPun_2012.eam.alloy Co
|
||||
|
||||
#pair_style pair/spin 4.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
|
||||
#pair_coeff * * exchange 4.0 0.0 0.003496 1.4885
|
||||
pair_coeff * * pair/spin/exchange 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.001 0.0 0.0 1.0
|
||||
#pair_coeff * * me 2.6 0.01 1.0 1.0 1.0
|
||||
|
||||
#type i and j | interaction type | cutoff | K1 (eV) | K2 (adim) | K3 (Ang) (for SOC)
|
||||
#pair_coeff * * pair/spin/soc/neel neel 4.0 -0.003330282 0.864159 2.13731
|
||||
|
||||
#Define a skin distance, update neigh list every
|
||||
#neighbor 1.0 bin
|
||||
#neigh_modify every 10 check yes delay 20
|
||||
neighbor 1.0 bin
|
||||
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 anisotropy 0.01 0.0 0.0 1.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
|
||||
|
||||
#Magnetic integration fix
|
||||
fix 3 all integration/spin mpi
|
||||
|
||||
#compute real time, total magnetization, magnetic energy, and spin temperature
|
||||
#Iteration | Time | Mx | My | Mz | |M| | Em | Tm
|
||||
|
||||
#Setting the timestep for the simulation (in ps)
|
||||
timestep 0.0001
|
||||
|
||||
##################
|
||||
#######run########
|
||||
##################
|
||||
|
||||
compute out_mag all compute/spin
|
||||
compute out_pe all pe
|
||||
compute out_ke all ke
|
||||
compute out_temp all temp
|
||||
|
||||
variable magz equal c_out_mag[4]
|
||||
variable magnorm equal c_out_mag[5]
|
||||
variable emag equal c_out_mag[6]
|
||||
variable tmag equal c_out_mag[7]
|
||||
variable mag_force equal f_1
|
||||
|
||||
thermo 10
|
||||
#thermo_style custom step time v_emag c_out_pe c_out_ke c_out_temp v_mag_force v_magnorm v_tmag etotal
|
||||
thermo_style custom step time v_tmag v_magz v_magnorm etotal
|
||||
thermo_modify format float %20.15g
|
||||
|
||||
#Dump the positions and spin directions of magnetic particles (vmd format)
|
||||
dump 1 all custom 500 dump_VSRSV.lammpstrj type x y z spx spy spz
|
||||
|
||||
#Running the simulations for N timesteps
|
||||
#run 100
|
||||
run 10000
|
||||
|
||||
126
examples/SPIN/in.spin.kagome
Normal file
126
examples/SPIN/in.spin.kagome
Normal file
@ -0,0 +1,126 @@
|
||||
###################
|
||||
#######Init########
|
||||
###################
|
||||
|
||||
clear
|
||||
#setting units to metal (Ang, picosecs, eV, ...):
|
||||
units metal
|
||||
|
||||
#setting dimension of the system (N=2 or 3):
|
||||
dimension 3
|
||||
|
||||
#setting boundary conditions. (p for periodic, f for fixed, ...)
|
||||
boundary p p p
|
||||
#boundary f f f
|
||||
|
||||
#setting atom_style to spin:
|
||||
atom_style spin
|
||||
|
||||
#Define sort for paramagnetic simulations (if no pair interaction)
|
||||
#atom_modify sort 1000 4.0
|
||||
|
||||
#why?
|
||||
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
|
||||
|
||||
#Test Kagome
|
||||
#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 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 5.0
|
||||
|
||||
#Creating a simulation box based on the specified region. Entries: number of atom types and box ref.
|
||||
create_box 1 box
|
||||
|
||||
#Creating atoms (or molecules) on a lattice, or a single atom (or molecule), ...
|
||||
#Entries: atom type,
|
||||
create_atoms 1 box
|
||||
|
||||
#Replicating NxNxN the entire set of atoms
|
||||
#replicate 1 1 1
|
||||
|
||||
|
||||
#######################
|
||||
#######Settings########
|
||||
#######################
|
||||
|
||||
#Setting one or more properties of one or more atoms.
|
||||
#Setting mass
|
||||
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
|
||||
|
||||
#Magnetic exchange interaction coefficient for bulk fcc Cobalt
|
||||
pair_style pair/spin 4.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.001 0.0 0.0 1.0
|
||||
#pair_coeff * * me 2.6 0.01 1.0 1.0 1.0
|
||||
|
||||
#Define a skin distance, update neigh list every
|
||||
#neighbor 1.0 bin
|
||||
#neigh_modify every 10 check yes delay 20
|
||||
neighbor 0.0 bin
|
||||
neigh_modify every 1 check no delay 0
|
||||
|
||||
#Magnetic field fix
|
||||
#Type | Intensity (T or eV) | Direction
|
||||
fix 1 all force/spin zeeman 1.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 1.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 mpi
|
||||
|
||||
#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 10 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag_Co_nodamp.dat format %20.16g
|
||||
|
||||
#Setting the timestep for the simulation (in ps)
|
||||
timestep 0.0001
|
||||
|
||||
##################
|
||||
#######run########
|
||||
##################
|
||||
|
||||
#Dump the positions and spin directions of magnetic particles (vmd format)
|
||||
dump 1 all custom 100 dump_spin_T100.lammpstrj type x y z spx spy spz
|
||||
|
||||
#Running the simulations for N timesteps
|
||||
run 2000
|
||||
#run 1
|
||||
|
||||
@ -74,11 +74,15 @@ void AtomVecSpin::grow(int n)
|
||||
type = memory->grow(atom->type,nmax,"atom:type");
|
||||
mask = memory->grow(atom->mask,nmax,"atom:mask");
|
||||
image = memory->grow(atom->image,nmax,"atom:image");
|
||||
|
||||
// allocating mech. quantities
|
||||
|
||||
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");
|
||||
|
||||
@ -58,7 +58,8 @@ class AtomVecSpin : public AtomVec {
|
||||
int write_data_hybrid(FILE *, double *);
|
||||
bigint memory_usage();
|
||||
|
||||
// clear mag. and mech. forces
|
||||
// clear magnetic and mechanic forces
|
||||
|
||||
void force_clear(int, size_t);
|
||||
|
||||
|
||||
@ -67,7 +68,7 @@ class AtomVecSpin : public AtomVec {
|
||||
int *type,*mask;
|
||||
imageint *image;
|
||||
double **x,**v,**f; // lattice quantities
|
||||
double *mumag,**sp, **fm; // spin quantities
|
||||
double *mumag,**sp,**fm; // spin quantities
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -36,7 +36,7 @@ using namespace MathConst;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg), mag(NULL)
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command");
|
||||
|
||||
@ -54,7 +54,6 @@ ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
ComputeSpin::~ComputeSpin()
|
||||
{
|
||||
memory->destroy(mag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -70,6 +69,12 @@ void ComputeSpin::init()
|
||||
void ComputeSpin::compute_vector()
|
||||
{
|
||||
int i, index;
|
||||
int countsp, countsptot;
|
||||
double mag[3], magtot[3];
|
||||
double magenergy, magenergytot;
|
||||
double tempnum, tempnumtot;
|
||||
double tempdenom, tempdenomtot;
|
||||
double spintemperature;
|
||||
|
||||
invoked_vector = update->ntimestep;
|
||||
|
||||
@ -94,6 +99,7 @@ void ComputeSpin::compute_vector()
|
||||
|
||||
// compute total magnetization and magnetic energy
|
||||
// compute spin temperature (Nurdin et al., Phys. Rev. E 61, 2000)
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (atom->mumag_flag && atom->sp_flag) {
|
||||
@ -144,9 +150,6 @@ void ComputeSpin::compute_vector()
|
||||
|
||||
void ComputeSpin::allocate()
|
||||
{
|
||||
memory->destroy(mag);
|
||||
memory->create(mag,4,"compute/spin:mag");
|
||||
memory->create(magtot,5,"compute/spin:mag");
|
||||
memory->create(vector,7,"compute/spin:vector");
|
||||
}
|
||||
|
||||
|
||||
@ -32,16 +32,7 @@ class ComputeSpin : public Compute {
|
||||
void compute_vector();
|
||||
|
||||
private:
|
||||
double *mag;
|
||||
double *magtot;
|
||||
double magenergy;
|
||||
double magenergytot;
|
||||
double tempnum,tempnumtot;
|
||||
double tempdenom,tempdenomtot;
|
||||
double spintemperature;
|
||||
double kb,hbar;
|
||||
int countsp;
|
||||
int countsptot;
|
||||
int usecenter;
|
||||
|
||||
void allocate();
|
||||
|
||||
@ -48,7 +48,7 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
|
||||
|
||||
if (narg < 7) error->all(FLERR,"Illegal force/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)
|
||||
// fix ID group force/spin magnitude (T or eV) style (zeeman or anisotropy) direction (3 cartesian coordinates)
|
||||
|
||||
// magnetic interactions coded for cartesian coordinates
|
||||
|
||||
@ -103,8 +103,6 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
|
||||
|
||||
FixForceSpin::~FixForceSpin()
|
||||
{
|
||||
memory->destroy(spi);
|
||||
memory->destroy(fmi);
|
||||
delete [] magstr;
|
||||
}
|
||||
|
||||
@ -136,7 +134,6 @@ void FixForceSpin::init()
|
||||
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
|
||||
}
|
||||
|
||||
// check variables
|
||||
if (magstr) {
|
||||
magvar = input->variable->find(magstr);
|
||||
if (magvar < 0)
|
||||
@ -149,10 +146,8 @@ void FixForceSpin::init()
|
||||
if (magfieldstyle != CONSTANT) varflag = EQUAL;
|
||||
|
||||
// set magnetic field components
|
||||
if (varflag == CONSTANT) set_magneticforce();
|
||||
|
||||
memory->create(spi,3,"forcespin:spi");
|
||||
memory->create(fmi,3,"forcespin:fmi");
|
||||
if (varflag == CONSTANT) set_magneticforce();
|
||||
|
||||
}
|
||||
|
||||
@ -174,6 +169,7 @@ void FixForceSpin::setup(int vflag)
|
||||
void FixForceSpin::post_force(int vflag)
|
||||
{
|
||||
// update gravity due to variables
|
||||
|
||||
if (varflag != CONSTANT) {
|
||||
modify->clearstep_compute();
|
||||
modify->addstep_compute(update->ntimestep + 1);
|
||||
@ -183,6 +179,7 @@ void FixForceSpin::post_force(int vflag)
|
||||
double **sp = atom->sp;
|
||||
double *mumag = atom->mumag;
|
||||
double **fm = atom->fm;
|
||||
double spi[3], fmi[3];
|
||||
const int nlocal = atom->nlocal;
|
||||
double scalar;
|
||||
|
||||
@ -214,17 +211,18 @@ void FixForceSpin::post_force(int vflag)
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixForceSpin::compute_zeeman(int i, double *fmi)
|
||||
|
||||
void FixForceSpin::compute_zeeman(int i, double fmi[3])
|
||||
{
|
||||
double *mumag = atom->mumag;
|
||||
fmi[0] += mumag[i]*hx;
|
||||
fmi[1] += mumag[i]*hy;
|
||||
fmi[2] += mumag[i]*hz;
|
||||
fmi[0] -= mumag[i]*hx;
|
||||
fmi[1] -= mumag[i]*hy;
|
||||
fmi[2] -= mumag[i]*hz;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixForceSpin::compute_anisotropy(int i, double * spi, double *fmi)
|
||||
void FixForceSpin::compute_anisotropy(int i, double spi[3], double fmi[3])
|
||||
{
|
||||
double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
|
||||
fmi[0] += scalar*Kax;
|
||||
@ -240,6 +238,7 @@ void FixForceSpin::post_force_respa(int vflag, int ilevel, int iloop)
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixForceSpin::set_magneticforce()
|
||||
{
|
||||
if (style == ZEEMAN) {
|
||||
@ -260,8 +259,10 @@ void FixForceSpin::set_magneticforce()
|
||||
|
||||
double FixForceSpin::compute_scalar()
|
||||
{
|
||||
double emag_all=0.0;
|
||||
double emag_all = 0.0;
|
||||
|
||||
// only sum across procs one time
|
||||
|
||||
if (eflag == 0) {
|
||||
MPI_Allreduce(&emag,&emag_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
eflag = 1;
|
||||
|
||||
@ -38,8 +38,8 @@ class FixForceSpin : public Fix {
|
||||
double compute_scalar();
|
||||
|
||||
int zeeman_flag, aniso_flag;
|
||||
void compute_zeeman(int, double *);
|
||||
void compute_anisotropy(int, double *, double *);
|
||||
void compute_zeeman(int, double fmi[3]);
|
||||
void compute_anisotropy(int, double spi[3], double fmi[3]);
|
||||
|
||||
protected:
|
||||
int style; // style of the magnetic force
|
||||
@ -57,18 +57,17 @@ class FixForceSpin : public Fix {
|
||||
char *magstr;
|
||||
|
||||
// zeeman field intensity and direction
|
||||
|
||||
double H_field;
|
||||
double nhx, nhy, nhz;
|
||||
double hx, hy, hz; // temp. force variables
|
||||
|
||||
// magnetic anisotropy intensity and direction
|
||||
|
||||
double Ka;
|
||||
double nax, nay, naz;
|
||||
double Kax, Kay, Kaz; // temp. force variables
|
||||
|
||||
// temp. spin variables
|
||||
double *spi, *fmi;
|
||||
|
||||
void set_magneticforce();
|
||||
|
||||
};
|
||||
|
||||
@ -55,7 +55,10 @@ enum{NONE,SPIN};
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixIntegrationSpin::FixIntegrationSpin(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
Fix(lmp, narg, arg),
|
||||
rsec(NULL), k(NULL), adv_list(NULL),
|
||||
lockpairspinexchange(NULL), lockpairspinsocneel(NULL), lockforcespin(NULL),
|
||||
locklangevinspin(NULL)
|
||||
{
|
||||
|
||||
if (narg != 4) error->all(FLERR,"Illegal fix integration/spin command");
|
||||
@ -78,37 +81,23 @@ FixIntegrationSpin::FixIntegrationSpin(LAMMPS *lmp, int narg, char **arg) :
|
||||
// error checks
|
||||
if (extra == SPIN && !atom->mumag_flag)
|
||||
error->all(FLERR,"Fix integration/spin requires spin attribute mumag");
|
||||
|
||||
magpair_flag = 0;
|
||||
soc_flag = 0;
|
||||
exch_flag = 0;
|
||||
exch_flag = soc_flag = 0;
|
||||
magforce_flag = 0;
|
||||
zeeman_flag = aniso_flag = 0;
|
||||
maglangevin_flag = 0;
|
||||
tdamp_flag = temp_flag = 0;
|
||||
|
||||
lockpairspin = NULL;
|
||||
lockpairspinexchange = NULL;
|
||||
lockpairspinsocneel = NULL;
|
||||
lockforcespin = NULL;
|
||||
locklangevinspin = NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
FixIntegrationSpin::~FixIntegrationSpin()
|
||||
{
|
||||
//delete lockpairspin;
|
||||
//delete lockforcespin;
|
||||
memory->destroy(xi);
|
||||
memory->destroy(k);
|
||||
memory->destroy(adv_list);
|
||||
|
||||
memory->destroy(sec);
|
||||
memory->destroy(rsec);
|
||||
memory->destroy(seci);
|
||||
|
||||
memory->destroy(rij);
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
|
||||
}
|
||||
|
||||
@ -127,21 +116,12 @@ int FixIntegrationSpin::setmask()
|
||||
void FixIntegrationSpin::init()
|
||||
{
|
||||
// set timesteps
|
||||
|
||||
dtv = update->dt;
|
||||
dtf = 0.5 * update->dt * force->ftm2v;
|
||||
dts = 0.25 * update->dt;
|
||||
|
||||
memory->create(xi,3,"integrations:xi");
|
||||
memory->create(sec,3,"integrations:sec");
|
||||
memory->create(rsec,3,"integrations:rsec");
|
||||
memory->create(seci,3,"integrations:seci");
|
||||
|
||||
memory->create(rij,3,"integrations:xi");
|
||||
memory->create(spi,3,"integrations:spi");
|
||||
memory->create(spj,3,"integrations:spj");
|
||||
memory->create(fmi,3,"integrations:fmi");
|
||||
memory->create(fmj,3,"integrations:fmj");
|
||||
|
||||
// set necessary pointers to interaction classes
|
||||
|
||||
if (strstr(force->pair_style,"pair/spin/exchange")) {
|
||||
exch_flag = 1;
|
||||
@ -192,9 +172,17 @@ void FixIntegrationSpin::init()
|
||||
}
|
||||
|
||||
|
||||
memory->create(rsec,3,"integration/spin:rsec");
|
||||
|
||||
// perform the sectoring if mpi integration
|
||||
|
||||
if (mpi_flag) sectoring();
|
||||
|
||||
// grow tables for k and adv_list
|
||||
|
||||
k = memory->grow(k,nsectors,"integration/spin:k");
|
||||
adv_list = memory->grow(adv_list,nsectors,atom->nmax,"integration/spin:adv_list");
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -202,6 +190,8 @@ void FixIntegrationSpin::init()
|
||||
void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
{
|
||||
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy;
|
||||
double xi[3];
|
||||
double spi[3], fmi[3];
|
||||
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
@ -215,8 +205,8 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
|
||||
// advance spin-lattice system, vsrsv
|
||||
// update half v for all particles
|
||||
// update half v for all atoms
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
@ -227,24 +217,22 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
}
|
||||
}
|
||||
|
||||
//#define SEC
|
||||
#define LIST
|
||||
// grow/shrink adv_list table to nlocal
|
||||
// init. tables for mpi-sector storage
|
||||
|
||||
#if defined LIST
|
||||
//printf("sectors = %d \n",nsectors);
|
||||
int adv_list[nsectors][nlocal];
|
||||
int k[nsectors];
|
||||
int p;
|
||||
adv_list = memory->grow(adv_list,nsectors,nlocal,"integration/spin:adv_list");
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
k[j] = 0;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
adv_list[j][i] = 0;
|
||||
}
|
||||
}
|
||||
int s, p;
|
||||
|
||||
// update half s for all particles
|
||||
// update half s for all atoms
|
||||
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) { // mpi seq. update
|
||||
if (mpi_flag) { // mpi seq. update
|
||||
int nseci;
|
||||
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
@ -256,89 +244,37 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts);
|
||||
adv_list[j][k[j]] = i;
|
||||
k[j]++;
|
||||
}
|
||||
}
|
||||
int ntest = 0;
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
ntest += k[j];
|
||||
}
|
||||
if (ntest != nlocal) error->all(FLERR,"error, S(k[j]) != nlocal");
|
||||
|
||||
for (int j = nsectors-1; j >= 0; j--) {
|
||||
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
for (int i = k[j]-1; i >= 0; i--) {
|
||||
p = adv_list[j][i];
|
||||
ComputeInteractionsSpin(p);
|
||||
AdvanceSingleSpin(p,dts,sp,fm);
|
||||
AdvanceSingleSpin(p,dts);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) { // serial seq. update
|
||||
comm->forward_comm(); // comm. positions of ghost atoms
|
||||
for (int i = 0; i < nlocal-1; i++){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts);
|
||||
}
|
||||
ComputeInteractionsSpin(nlocal-1);
|
||||
AdvanceSingleSpin(nlocal-1,2.0*dts,sp,fm); // advance half s for 1
|
||||
AdvanceSingleSpin(nlocal-1,2.0*dts); // advance half s for 1
|
||||
for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#if defined SEC
|
||||
// update half s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) { // mpi seq. update
|
||||
int nseci;
|
||||
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) { // serial seq. update
|
||||
comm->forward_comm(); // comm. positions of ghost atoms
|
||||
for (int i = 0; i < nlocal-1; i++){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
ComputeInteractionsSpin(nlocal-1);
|
||||
AdvanceSingleSpin(nlocal-1,2.0*dts,sp,fm); // advance half s for 1
|
||||
for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
// update x for all particles
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
x[i][0] += dtv * v[i][0];
|
||||
@ -348,17 +284,17 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
}
|
||||
|
||||
|
||||
#if defined LIST
|
||||
// update half s for all particles
|
||||
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) { // mpi seq. update
|
||||
if (mpi_flag) { // mpi seq. update
|
||||
int nseci;
|
||||
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < k[j]; i++) {
|
||||
p = adv_list[j][i];
|
||||
ComputeInteractionsSpin(p);
|
||||
AdvanceSingleSpin(p,dts,sp,fm);
|
||||
AdvanceSingleSpin(p,dts);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
|
||||
@ -366,69 +302,24 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
for (int i = k[j]-1; i >= 0; i--) {
|
||||
p = adv_list[j][i];
|
||||
ComputeInteractionsSpin(p);
|
||||
AdvanceSingleSpin(p,dts,sp,fm);
|
||||
AdvanceSingleSpin(p,dts);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) { // serial seq. update
|
||||
comm->forward_comm(); // comm. positions of ghost atoms
|
||||
for (int i = 0; i < nlocal-1; i++){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts);
|
||||
}
|
||||
ComputeInteractionsSpin(nlocal-1);
|
||||
AdvanceSingleSpin(nlocal-1,2.0*dts,sp,fm); // advance half s for 1
|
||||
AdvanceSingleSpin(nlocal-1,2.0*dts); // advance half s for 1
|
||||
for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
#endif
|
||||
|
||||
#if defined SEC
|
||||
// update half s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) { // mpi seq. update
|
||||
int nseci;
|
||||
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
|
||||
comm->forward_comm();
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) { // serial seq. update
|
||||
comm->forward_comm(); // comm. positions of ghost atoms
|
||||
for (int i = 0; i < nlocal-1; i++){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
ComputeInteractionsSpin(nlocal-1);
|
||||
AdvanceSingleSpin(nlocal-1,2.0*dts,sp,fm); // advance half s for 1
|
||||
for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
@ -436,8 +327,10 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
||||
void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
{
|
||||
const int nlocal = atom->nlocal;
|
||||
double xi[3], rij[3];
|
||||
double spi[3], spj[3];
|
||||
double fmi[3], fmj[3];
|
||||
|
||||
// force compute quantities
|
||||
int i,j,jj,inum,jnum,itype,jtype;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double **x = atom->x;
|
||||
@ -455,7 +348,6 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
}
|
||||
|
||||
double rsq, rd;
|
||||
double delx, dely, delz;
|
||||
double temp_cut, cut_2, inorm;
|
||||
temp_cut = cut_2 = inorm = 0.0;
|
||||
|
||||
@ -464,6 +356,7 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
int pair_compute_flag = 1;
|
||||
|
||||
// force computation for spin i
|
||||
|
||||
i = ilist[ii];
|
||||
|
||||
spi[0] = sp[i][0];
|
||||
@ -478,7 +371,8 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
// pair interaction
|
||||
// loop on neighbors for pair interactions
|
||||
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
|
||||
j = jlist[jj];
|
||||
@ -490,15 +384,15 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
delx = x[j][0] - xi[0];
|
||||
dely = x[j][1] - xi[1];
|
||||
delz = x[j][2] - xi[2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
|
||||
rij[0] = delx*inorm;
|
||||
rij[1] = dely*inorm;
|
||||
rij[2] = delz*inorm;
|
||||
rij[0] = x[j][0] - xi[0];
|
||||
rij[1] = x[j][1] - xi[1];
|
||||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
|
||||
temp_cut = 0.0;
|
||||
|
||||
@ -514,7 +408,7 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
temp_cut = lockpairspinsocneel->cut_soc_neel[itype][jtype];
|
||||
cut_2 = temp_cut*temp_cut;
|
||||
if (rsq <= cut_2) {
|
||||
lockpairspinsocneel->compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
//lockpairspinsocneel->compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
@ -525,23 +419,21 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
lockforcespin->compute_zeeman(i,fmi);
|
||||
}
|
||||
if (aniso_flag) { // aniso
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
lockforcespin->compute_anisotropy(i,spi,fmi);
|
||||
}
|
||||
}
|
||||
|
||||
if (maglangevin_flag) { // mag. langevin
|
||||
if (tdamp_flag) { // trans. damping
|
||||
if (tdamp_flag) { // transverse damping
|
||||
locklangevinspin->add_tdamping(spi,fmi);
|
||||
}
|
||||
if (temp_flag) { // temp.
|
||||
if (temp_flag) { // spin temperature
|
||||
locklangevinspin->add_temperature(fmi);
|
||||
}
|
||||
}
|
||||
|
||||
// replace the mag. force i by its new value
|
||||
// replace the magnethic force i by its new value
|
||||
|
||||
fm[i][0] = fmi[0];
|
||||
fm[i][1] = fmi[1];
|
||||
fm[i][2] = fmi[2];
|
||||
@ -551,6 +443,7 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixIntegrationSpin::sectoring()
|
||||
{
|
||||
int sec[3];
|
||||
double sublo[3],subhi[3];
|
||||
double* sublotmp = domain->sublo;
|
||||
double* subhitmp = domain->subhi;
|
||||
@ -591,18 +484,20 @@ void FixIntegrationSpin::sectoring()
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int FixIntegrationSpin::coords2sector(double *xi)
|
||||
|
||||
int FixIntegrationSpin::coords2sector(double x[3])
|
||||
{
|
||||
int nseci;
|
||||
int seci[3];
|
||||
double sublo[3];
|
||||
double* sublotmp = domain->sublo;
|
||||
for (int dim = 0 ; dim<3 ; dim++) {
|
||||
sublo[dim]=sublotmp[dim];
|
||||
}
|
||||
|
||||
double rix = (xi[0] - sublo[0])/rsec[0];
|
||||
double riy = (xi[1] - sublo[1])/rsec[1];
|
||||
double riz = (xi[2] - sublo[2])/rsec[2];
|
||||
double rix = (x[0] - sublo[0])/rsec[0];
|
||||
double riy = (x[1] - sublo[1])/rsec[1];
|
||||
double riz = (x[2] - sublo[2])/rsec[2];
|
||||
|
||||
seci[0] = static_cast<int>(rix);
|
||||
seci[1] = static_cast<int>(riy);
|
||||
@ -616,10 +511,12 @@ int FixIntegrationSpin::coords2sector(double *xi)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixIntegrationSpin::AdvanceSingleSpin(int i, double dtl, double **sp, double **fm)
|
||||
void FixIntegrationSpin::AdvanceSingleSpin(int i, double dtl)
|
||||
{
|
||||
int j=0;
|
||||
int *sametag = atom->sametag;
|
||||
double **sp = atom->sp;
|
||||
double **fm = atom->fm;
|
||||
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy,dts2;
|
||||
double cp[3],g[3];
|
||||
|
||||
@ -651,13 +548,15 @@ void FixIntegrationSpin::AdvanceSingleSpin(int i, double dtl, double **sp, doubl
|
||||
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;
|
||||
|
||||
// comm. sp[i] to atoms with same tag (serial algo)
|
||||
// comm. sp[i] to atoms with same tag (for serial algo)
|
||||
|
||||
if (mpi_flag == 0) {
|
||||
if (sametag[i] >= 0) {
|
||||
j = sametag[i];
|
||||
@ -692,6 +591,7 @@ void FixIntegrationSpin::final_integrate()
|
||||
int *mask = atom->mask;
|
||||
|
||||
// update half v for all particles
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
|
||||
@ -34,46 +34,35 @@ class FixIntegrationSpin : public Fix {
|
||||
virtual void initial_integrate(int);
|
||||
virtual void final_integrate();
|
||||
|
||||
// compute and advance single spin
|
||||
void ComputeInteractionsSpin(int);
|
||||
void AdvanceSingleSpin(int, double, double **, double **);
|
||||
void ComputeInteractionsSpin(int); // compute and advance single spin functions
|
||||
void AdvanceSingleSpin(int, double);
|
||||
|
||||
// sectoring operations
|
||||
void sectoring();
|
||||
int coords2sector(double *);
|
||||
void sectoring(); // sectoring operation functions
|
||||
int coords2sector(double x[3]);
|
||||
|
||||
protected:
|
||||
int extra, mpi_flag;
|
||||
|
||||
// velocity, force, and spin timesteps
|
||||
double dtv,dtf,dts;
|
||||
double dtv,dtf,dts; // velocity, force, and spin timesteps
|
||||
|
||||
// mag. interaction flags
|
||||
int magpair_flag;
|
||||
int soc_flag;
|
||||
int exch_flag;
|
||||
int magforce_flag;
|
||||
int magpair_flag; // magnetic pair flags
|
||||
int soc_flag, exch_flag;
|
||||
int magforce_flag; // magnetic force flags
|
||||
int zeeman_flag, aniso_flag;
|
||||
int maglangevin_flag;
|
||||
int maglangevin_flag; // magnetic langevin flags
|
||||
int tdamp_flag, temp_flag;
|
||||
|
||||
// pointers to interaction classes
|
||||
// pointers to magnetic interaction classes
|
||||
|
||||
class PairHybrid *lockhybrid;
|
||||
class PairSpin *lockpairspin;
|
||||
class PairSpinExchange *lockpairspinexchange;
|
||||
class PairSpinSocNeel *lockpairspinsocneel;
|
||||
class FixForceSpin *lockforcespin;
|
||||
class FixLangevinSpin *locklangevinspin;
|
||||
|
||||
// temporary variables
|
||||
double *xi, *rij;
|
||||
double *spi, *spj;
|
||||
double *fmi, *fmj;
|
||||
|
||||
// sectoring variables
|
||||
int nsectors;
|
||||
int *sec, *seci;
|
||||
int nsectors; // sectoring variables
|
||||
double *rsec;
|
||||
int *k, **adv_list;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -89,7 +89,7 @@ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) :
|
||||
}
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
//random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
random = new RanPark(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
@ -156,76 +156,37 @@ void FixLangevinSpin::setup(int vflag)
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
/*
|
||||
void FixLangevinSpin::post_force(int vflag)
|
||||
{
|
||||
double **sp = atom->sp;
|
||||
double **fm = atom->fm;
|
||||
int *mask = atom->mask;
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
// add the damping to the effective field of each spin
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
fmi[0] = fm[i][0];
|
||||
fmi[1] = fm[i][1];
|
||||
fmi[2] = fm[i][2];
|
||||
|
||||
if (tdamp_flag) {
|
||||
add_tdamping(spi,fmi);
|
||||
}
|
||||
|
||||
if (temp_flag) {
|
||||
add_temperature(fmi);
|
||||
}
|
||||
|
||||
fm[i][0] = fmi[0];
|
||||
fm[i][1] = fmi[1];
|
||||
fm[i][2] = fmi[2];
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
*/
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixLangevinSpin::add_tdamping(double *spi, double *fmi)
|
||||
void FixLangevinSpin::add_tdamping(double spi[3], double fmi[3])
|
||||
{
|
||||
double cpx = fmi[1]*spi[2] - fmi[2]*spi[1];
|
||||
double cpy = fmi[2]*spi[0] - fmi[0]*spi[2];
|
||||
double cpz = fmi[0]*spi[1] - fmi[1]*spi[0];
|
||||
|
||||
// taking away the transverse damping
|
||||
// adding the transverse damping
|
||||
|
||||
fmi[0] -= alpha_t*cpx;
|
||||
fmi[1] -= alpha_t*cpy;
|
||||
fmi[2] -= alpha_t*cpz;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixLangevinSpin::add_temperature(double *fmi)
|
||||
|
||||
void FixLangevinSpin::add_temperature(double fmi[3])
|
||||
{
|
||||
//#define GAUSSIAN_R
|
||||
#if defined GAUSSIAN_R
|
||||
// drawing gaussian random dist
|
||||
double rx = sigma*random->gaussian();
|
||||
double ry = sigma*random->gaussian();
|
||||
double rz = sigma*random->gaussian();
|
||||
#else
|
||||
|
||||
double rx = sigma*(-1.0+2.0*random->uniform());
|
||||
double ry = sigma*(-1.0+2.0*random->uniform());
|
||||
double rz = sigma*(-1.0+2.0*random->uniform());
|
||||
#endif
|
||||
|
||||
// adding the random field
|
||||
|
||||
fmi[0] += rx;
|
||||
fmi[1] += ry;
|
||||
fmi[2] += rz;
|
||||
|
||||
// adding Gilbert's prefactor
|
||||
// adding gilbert's prefactor
|
||||
|
||||
fmi[0] *= gil_factor;
|
||||
fmi[1] *= gil_factor;
|
||||
fmi[2] *= gil_factor;
|
||||
|
||||
@ -33,25 +33,22 @@ class FixLangevinSpin : public Fix {
|
||||
void setup(int);
|
||||
// virtual void post_force(int);
|
||||
void post_force_respa(int, int, int);
|
||||
void add_tdamping(double *, double *); // add transverse damping
|
||||
void add_temperature(double *); // add temperature
|
||||
void add_tdamping(double spi[3], double fmi[3]); // add transverse damping
|
||||
void add_temperature(double fmi[3]); // add temperature
|
||||
int tdamp_flag, ldamp_flag, temp_flag; // damping and temperature flags
|
||||
|
||||
protected:
|
||||
double *spi, *fmi;
|
||||
double alpha_t; // mag. transverse damping coeff.
|
||||
double alpha_l; // mag. longitudinal damping coeff.
|
||||
double dts; // timestep
|
||||
double alpha_t, alpha_l; // transverse and longitudunal damping coeff.
|
||||
double dts; // magnetic timestep
|
||||
double temp; // spin bath temperature
|
||||
double D,sigma; // bath intensity var.
|
||||
|
||||
double gil_factor; // Gilbert's prefactor
|
||||
double gil_factor; // gilbert's prefactor
|
||||
|
||||
char *id_temp;
|
||||
class Compute *temperature;
|
||||
|
||||
int nlevels_respa;
|
||||
//class RanMars *random;
|
||||
class RanPark *random;
|
||||
int seed;
|
||||
|
||||
|
||||
@ -1,637 +0,0 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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: Julien Tranchida (SNL)
|
||||
Aidan Thompson (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "pair_spin.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
single_enable = 0;
|
||||
exch_flag = 0;
|
||||
dmi_flag = 0;
|
||||
me_flag = 0;
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpin::~PairSpin()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_spin_exchange);
|
||||
memory->destroy(J1_mag);
|
||||
memory->destroy(J1_mech);
|
||||
memory->destroy(J2);
|
||||
memory->destroy(J3);
|
||||
|
||||
memory->destroy(cut_spin_dmi);
|
||||
memory->destroy(DM);
|
||||
memory->destroy(v_dmx);
|
||||
memory->destroy(v_dmy);
|
||||
memory->destroy(v_dmz);
|
||||
|
||||
memory->destroy(cut_spin_me);
|
||||
memory->destroy(ME);
|
||||
memory->destroy(v_mex);
|
||||
memory->destroy(v_mey);
|
||||
memory->destroy(v_mez);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(del);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xtmp,ytmp,ztmp;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_ex_2,cut_dmi_2,cut_me_2,cut_spin_pair_global2;
|
||||
double rsq,rd;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
cut_spin_pair_global2 = cut_spin_pair_global*cut_spin_pair_global;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
double *mumag = atom->mumag;
|
||||
double **sp = atom->sp;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// pair spin computations
|
||||
// loop over neighbors of my 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];
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
fi[0] = fi[1] = fi[2] = 0.0;
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
del[0] = del[1] = del[2] = 0.0;
|
||||
|
||||
del[0] = x[j][0] - xtmp;
|
||||
del[1] = x[j][1] - ytmp;
|
||||
del[2] = x[j][2] - ztmp;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rsq = del[0]*del[0] + del[1]*del[1] + del[2]*del[2];
|
||||
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// exchange interaction
|
||||
if (exch_flag) {
|
||||
cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
|
||||
if (rsq <= cut_ex_2) {
|
||||
compute_exchange(i,j,rsq,fmi,fmj,spi,spj);
|
||||
compute_exchange_mech(i,j,rsq,del,fi,fj,spi,spj);
|
||||
}
|
||||
}
|
||||
// dm interaction
|
||||
if (dmi_flag){
|
||||
cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
|
||||
if (rsq <= cut_dmi_2){
|
||||
compute_dmi(i,j,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
// me interaction
|
||||
if (me_flag){
|
||||
cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
|
||||
if (rsq <= cut_me_2){
|
||||
compute_me(i,j,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
/* if (newton_pair || j < nlocal) {*/
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_spin_pair_global2) {
|
||||
evdwl -= mumag[i]*spi[0]*fmi[0];
|
||||
evdwl -= mumag[i]*spi[1]*fmi[1];
|
||||
evdwl -= mumag[i]*spi[2]*fmi[2];
|
||||
}
|
||||
}
|
||||
|
||||
evdwl *= hbar;
|
||||
|
||||
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
|
||||
evdwl,ecoul,fi[0],fi[1],fi[2],del[0],del[1],del[2]);
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpin::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double Jex, ra;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
ra = rsq/J3[itype][jtype]/J3[itype][jtype];
|
||||
Jex = 4.0*J1_mag[itype][jtype]*ra;
|
||||
Jex *= (1.0-J2[itype][jtype]*ra);
|
||||
Jex *= exp(-ra);
|
||||
|
||||
fmi[0] += Jex*spj[0];
|
||||
fmi[1] += Jex*spj[1];
|
||||
fmi[2] += Jex*spj[2];
|
||||
|
||||
fmj[0] += Jex*spi[0];
|
||||
fmj[1] += Jex*spi[1];
|
||||
fmj[2] += Jex*spi[2];
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpin::compute_exchange_mech(int i, int j, double rsq, double *del, double *fi, double *fj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double Jex, Jex_mech, ra, rr;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
Jex = J1_mech[itype][jtype];
|
||||
|
||||
ra = rsq/J3[itype][jtype]/J3[itype][jtype];
|
||||
rr = sqrt(rsq)/J3[itype][jtype]/J3[itype][jtype];
|
||||
|
||||
Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra);
|
||||
Jex_mech *= 8.0*Jex*rr*exp(-ra);
|
||||
Jex_mech *= (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
|
||||
|
||||
fi[0] += Jex_mech*del[0];
|
||||
fi[1] += Jex_mech*del[1];
|
||||
fi[2] += Jex_mech*del[2];
|
||||
|
||||
fj[0] -= Jex_mech*del[0];
|
||||
fj[1] -= Jex_mech*del[1];
|
||||
fj[2] -= Jex_mech*del[2];
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpin::compute_dmi(int i, int j, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double dmix,dmiy,dmiz;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
dmix = DM[itype][jtype]*v_dmx[itype][jtype];
|
||||
dmiy = DM[itype][jtype]*v_dmy[itype][jtype];
|
||||
dmiz = DM[itype][jtype]*v_dmz[itype][jtype];
|
||||
|
||||
fmi[0] += spj[1]*dmiz-spj[2]*dmiy;
|
||||
fmi[1] += spj[2]*dmix-spj[0]*dmiz;
|
||||
fmi[2] += spj[0]*dmiy-spj[1]*dmix;
|
||||
|
||||
fmj[0] -= spi[1]*dmiz-spi[2]*dmiy;
|
||||
fmj[1] -= spi[2]*dmix-spi[0]*dmiz;
|
||||
fmj[2] -= spi[0]*dmiy-spi[1]*dmix;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpin::compute_me(int i, int j, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
double **sp = atom->sp;
|
||||
double **x = atom->x;
|
||||
double meix,meiy,meiz;
|
||||
double rx, ry, rz, inorm;
|
||||
|
||||
rx = x[j][0] - x[i][0];
|
||||
ry = x[j][1] - x[i][1];
|
||||
rz = x[j][2] - x[i][2];
|
||||
inorm = 1.0/sqrt(rx*rx+ry*ry+rz*rz);
|
||||
rx *= inorm;
|
||||
ry *= inorm;
|
||||
rz *= inorm;
|
||||
|
||||
meix = v_mey[itype][jtype]*rz - v_mez[itype][jtype]*ry;
|
||||
meiy = v_mez[itype][jtype]*rx - v_mex[itype][jtype]*rz;
|
||||
meiz = v_mex[itype][jtype]*ry - v_mey[itype][jtype]*rx;
|
||||
|
||||
meix *= ME[itype][jtype];
|
||||
meiy *= ME[itype][jtype];
|
||||
meiz *= ME[itype][jtype];
|
||||
|
||||
fmi[0] += spj[1]*meiz - spj[2]*meiy;
|
||||
fmi[1] += spj[2]*meix - spj[0]*meiz;
|
||||
fmi[2] += spj[0]*meiy - spj[1]*meix;
|
||||
|
||||
fmj[0] -= spi[1]*meiz - spi[2]*meiy;
|
||||
fmj[1] -= spi[2]*meix - spi[0]*meiz;
|
||||
fmj[2] -= spi[0]*meiy - spi[1]*meix;
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
memory->create(setflag,n+1,n+1,"pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
|
||||
memory->create(cut_spin_exchange,n+1,n+1,"pair:cut_spin_exchange");
|
||||
memory->create(J1_mag,n+1,n+1,"pair:J1_mag");
|
||||
memory->create(J1_mech,n+1,n+1,"pair:J1_mech");
|
||||
memory->create(J2,n+1,n+1,"pair:J2");
|
||||
memory->create(J3,n+1,n+1,"pair:J3");
|
||||
|
||||
memory->create(cut_spin_dmi,n+1,n+1,"pair:cut_spin_dmi");
|
||||
memory->create(DM,n+1,n+1,"pair:DM");
|
||||
memory->create(v_dmx,n+1,n+1,"pair:DM_vector_x");
|
||||
memory->create(v_dmy,n+1,n+1,"pair:DM_vector_y");
|
||||
memory->create(v_dmz,n+1,n+1,"pair:DM_vector_z");
|
||||
|
||||
memory->create(cut_spin_me,n+1,n+1,"pair:cut_spin_me");
|
||||
memory->create(ME,n+1,n+1,"pair:ME");
|
||||
memory->create(v_mex,n+1,n+1,"pair:ME_vector_x");
|
||||
memory->create(v_mey,n+1,n+1,"pair:ME_vector_y");
|
||||
memory->create(v_mez,n+1,n+1,"pair:ME_vector_z");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(del,3,"pair:del");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_spin_pair_global = force->numeric(FLERR,arg[0]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i+1; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j]) {
|
||||
cut_spin_exchange[i][j] = cut_spin_pair_global;
|
||||
cut_spin_dmi[i][j] = cut_spin_pair_global;
|
||||
cut_spin_me[i][j] = cut_spin_pair_global;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type spin pairs (only one for now)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::coeff(int narg, char **arg)
|
||||
{
|
||||
const double hbar = force->hplanck/MY_2PI;
|
||||
|
||||
if (!allocated) allocate();
|
||||
|
||||
if (strcmp(arg[2],"exchange")==0){
|
||||
if (narg != 7) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
exch_flag = 1;
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double j1 = (force->numeric(FLERR,arg[4]));
|
||||
const double j2 = force->numeric(FLERR,arg[5]);
|
||||
const double j3 = force->numeric(FLERR,arg[6]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_exchange[i][j] = rij;
|
||||
J1_mag[i][j] = j1/hbar;
|
||||
J1_mech[i][j] = j1;
|
||||
J2[i][j] = j2;
|
||||
J3[i][j] = j3;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else if (strcmp(arg[2],"dmi")==0) {
|
||||
if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
dmi_flag = 1;
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double dm = (force->numeric(FLERR,arg[4]))/hbar;
|
||||
double dmx = force->numeric(FLERR,arg[5]);
|
||||
double dmy = force->numeric(FLERR,arg[6]);
|
||||
double dmz = force->numeric(FLERR,arg[7]);
|
||||
|
||||
double inorm = 1.0/(dmx*dmx+dmy*dmy+dmz*dmz);
|
||||
dmx *= inorm;
|
||||
dmy *= inorm;
|
||||
dmz *= inorm;
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_dmi[i][j] = rij;
|
||||
DM[i][j] = dm;
|
||||
v_dmx[i][j] = dmx;
|
||||
v_dmy[i][j] = dmy;
|
||||
v_dmz[i][j] = dmz;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else if (strcmp(arg[2],"me")==0) {
|
||||
if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
me_flag = 1;
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double me = (force->numeric(FLERR,arg[4]))/hbar;
|
||||
double mex = force->numeric(FLERR,arg[5]);
|
||||
double mey = force->numeric(FLERR,arg[6]);
|
||||
double mez = force->numeric(FLERR,arg[7]);
|
||||
|
||||
double inorm = 1.0/(mex*mex+mey*mey+mez*mez);
|
||||
mex *= inorm;
|
||||
mey *= inorm;
|
||||
mez *= inorm;
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_me[i][j] = rij;
|
||||
DM[i][j] = me;
|
||||
v_mex[i][j] = mex;
|
||||
v_mey[i][j] = mey;
|
||||
v_mez[i][j] = mez;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
// check if Jex [][] still works for Ferrimagnetic exchange
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::init_style()
|
||||
{
|
||||
if (!atom->sp_flag || !atom->mumag_flag)
|
||||
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
#endif
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpin::init_one(int i, int j)
|
||||
{
|
||||
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
return cut_spin_pair_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j]) {
|
||||
if (exch_flag){
|
||||
fwrite(&J1_mag[i][j],sizeof(double),1,fp);
|
||||
fwrite(&J1_mech[i][j],sizeof(double),1,fp);
|
||||
fwrite(&J2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&J3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
if (dmi_flag) {
|
||||
fwrite(&DM[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_dmx[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_dmy[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_dmz[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_spin_dmi[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
if (me_flag) {
|
||||
fwrite(&ME[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_mex[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_mey[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_mez[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_spin_me[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
|
||||
allocate();
|
||||
|
||||
int i,j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++) {
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) {
|
||||
fread(&J1_mag[i][j],sizeof(double),1,fp);
|
||||
fread(&J1_mech[i][j],sizeof(double),1,fp);
|
||||
fread(&J2[i][j],sizeof(double),1,fp);
|
||||
fread(&J2[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_spin_exchange[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&J1_mag[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&J1_mech[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&J2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&J3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_spin_exchange[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_spin_pair_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpin::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_spin_pair_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_spin_pair_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
@ -1,99 +0,0 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(pair/spin,PairSpin)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SPIN_H
|
||||
#define LMP_PAIR_SPIN_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSpin : public Pair {
|
||||
public:
|
||||
PairSpin(class LAMMPS *);
|
||||
virtual ~PairSpin();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_exchange(int, int, double, double *, double *,double *, double *);
|
||||
void compute_exchange_mech(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_dmi(int, int, double *, double *, double *, double *);
|
||||
void compute_me(int, int, double *, double *, double *, double *);
|
||||
|
||||
int exch_flag, dmi_flag, me_flag;
|
||||
double cut_spin_pair_global;
|
||||
double cut_spin_dipolar_global;
|
||||
|
||||
double **cut_spin_exchange; // cutoff distance exchange
|
||||
double **cut_spin_dmi; // cutoff distance dmi
|
||||
double **cut_spin_me; // cutoff distance me
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
double hbar;
|
||||
|
||||
double **J1_mag, **J1_mech; // exchange coeffs Jij
|
||||
double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang
|
||||
|
||||
double **DM; // dmi coeff in eV
|
||||
double **v_dmx, **v_dmy, **v_dmz;// dmi direction
|
||||
|
||||
double **ME; // me coeff in eV
|
||||
double **v_mex, **v_mey, **v_mez;// me direction
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *del; // inter atomic distances
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args in pair_spin command
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Spin simulations require metal unit style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair spin requires atom attributes sp, mumag
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
@ -42,7 +42,7 @@ PairSpinExchange::PairSpinExchange(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
newton_pair_spin = 0; // no newton pair for now => to be corrected
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
@ -65,14 +65,6 @@ PairSpinExchange::~PairSpinExchange()
|
||||
memory->destroy(J2);
|
||||
memory->destroy(J3);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
@ -83,12 +75,13 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double xi[3], rij[3];
|
||||
double fi[3], fmi[3];
|
||||
double fj[3], fmj[3];
|
||||
double cut_ex_2,cut_spin_exchange_global2;
|
||||
double rsq,rd;
|
||||
double rsq,rd,inorm;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double spi[3],spj[3];
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
@ -109,14 +102,14 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// pair spin computations
|
||||
// loop over neighbors of my atoms
|
||||
// computation of the exchange interaction
|
||||
// loop over atoms and their neighbors
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
@ -124,6 +117,7 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
@ -138,15 +132,12 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rij[0] = x[j][0] - xi[0];
|
||||
rij[1] = x[j][1] - xi[1];
|
||||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
@ -154,7 +145,8 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// exchange interaction
|
||||
// compute exchange interaction
|
||||
|
||||
if (exch_flag) {
|
||||
cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
|
||||
if (rsq <= cut_ex_2) {
|
||||
@ -170,14 +162,15 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// check newton pair => needs correction
|
||||
// if (newton_pair || j < nlocal) {
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
//fm[j][0] += fmj[0];
|
||||
//fm[j][1] += fmj[1];
|
||||
//fm[j][2] += fmj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
@ -199,7 +192,8 @@ void PairSpinExchange::compute(int eflag, int vflag)
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinExchange::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj, double *spi, double *spj)
|
||||
|
||||
void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3], double fmj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -223,7 +217,8 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double *fmi,
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double *rij, double *fi, double *fj, double *spi, double *spj)
|
||||
|
||||
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double rij[3], double fi[3], double fj[3], double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -270,14 +265,6 @@ void PairSpinExchange::allocate()
|
||||
memory->create(J2,n+1,n+1,"pair:J2");
|
||||
memory->create(J3,n+1,n+1,"pair:J3");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
@ -320,6 +307,7 @@ void PairSpinExchange::coeff(int narg, char **arg)
|
||||
if (!allocated) allocate();
|
||||
|
||||
// set exch_mech_flag to 1 if magneto-mech simulation
|
||||
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
exch_mech_flag = 0;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
@ -335,15 +323,15 @@ void PairSpinExchange::coeff(int narg, char **arg)
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double j1 = (force->numeric(FLERR,arg[4]));
|
||||
const double rc = force->numeric(FLERR,arg[3]);
|
||||
const double j1 = force->numeric(FLERR,arg[4]);
|
||||
const double j2 = force->numeric(FLERR,arg[5]);
|
||||
const double j3 = force->numeric(FLERR,arg[6]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_exchange[i][j] = rij;
|
||||
cut_spin_exchange[i][j] = rc;
|
||||
J1_mag[i][j] = j1/hbar;
|
||||
if (exch_mech_flag) {
|
||||
J1_mech[i][j] = j1;
|
||||
@ -373,7 +361,7 @@ void PairSpinExchange::init_style()
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
// check this half/full request => needs correction/review
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
|
||||
@ -39,14 +39,14 @@ class PairSpinExchange : public Pair {
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_exchange(int, int, double, double *, double *,double *, double *);
|
||||
void compute_exchange_mech(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_exchange(int, int, double, double fmi[3], double fmj[3],double spi[3], double spj[3]);
|
||||
void compute_exchange_mech(int, int, double, double rij[3], double fmi[3], double fmj[3], double spi[3], double spj[3]);
|
||||
|
||||
int exch_flag; // mag. exchange flag
|
||||
int exch_mech_flag; // mech. exchange flags
|
||||
int exch_flag; // magnetic exchange flag
|
||||
int exch_mech_flag; // mechanic exchange flags
|
||||
|
||||
double cut_spin_exchange_global; // global exchange cutoff
|
||||
double **cut_spin_exchange; // cutoff distance exchange
|
||||
double **cut_spin_exchange; // cutoff distance per exchange
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
@ -55,11 +55,6 @@ class PairSpinExchange : public Pair {
|
||||
double **J1_mag, **J1_mech; // exchange coeffs Jij
|
||||
double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
|
||||
@ -64,14 +64,6 @@ PairSpinMe::~PairSpinMe()
|
||||
memory->destroy(v_mey);
|
||||
memory->destroy(v_mez);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
@ -81,12 +73,13 @@ PairSpinMe::~PairSpinMe()
|
||||
void PairSpinMe::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_me_2,cut_spin_me_global2;
|
||||
double rsq,rd;
|
||||
double evdwl, ecoul;
|
||||
double xi[3], rij[3];
|
||||
double spi[3], spj[3];
|
||||
double fi[3], fj[3];
|
||||
double fmi[3], fmj[3];
|
||||
double cut_me_2, cut_spin_me_global2;
|
||||
double rsq, rd, inorm;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -109,13 +102,13 @@ void PairSpinMe::compute(int eflag, int vflag)
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// magneto-electric computation
|
||||
// loop over neighbors of my atoms
|
||||
// loop over atoms and their neighbors
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
@ -123,6 +116,7 @@ void PairSpinMe::compute(int eflag, int vflag)
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
@ -139,13 +133,11 @@ void PairSpinMe::compute(int eflag, int vflag)
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rij[0] = x[j][0] - xi[0];
|
||||
rij[1] = x[j][1] - xi[1];
|
||||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
@ -154,6 +146,7 @@ void PairSpinMe::compute(int eflag, int vflag)
|
||||
jtype = type[j];
|
||||
|
||||
// me interaction
|
||||
|
||||
if (me_flag){
|
||||
cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
|
||||
if (rsq <= cut_me_2){
|
||||
@ -168,7 +161,7 @@ void PairSpinMe::compute(int eflag, int vflag)
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
// if (newton_pair || j < nlocal) { => to be corrected
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
@ -198,13 +191,13 @@ void PairSpinMe::compute(int eflag, int vflag)
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinMe::compute_me(int i, int j, double *fmi, double *fmj, double *spi, double *spj)
|
||||
|
||||
void PairSpinMe::compute_me(int i, int j, double fmi[3], double fmj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
double **sp = atom->sp;
|
||||
double **x = atom->x;
|
||||
double meix,meiy,meiz;
|
||||
double rx, ry, rz, inorm;
|
||||
@ -255,14 +248,6 @@ void PairSpinMe::allocate()
|
||||
memory->create(v_mey,n+1,n+1,"pair:ME_vector_y");
|
||||
memory->create(v_mez,n+1,n+1,"pair:ME_vector_z");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
@ -351,7 +336,7 @@ void PairSpinMe::init_style()
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
// check this half/full request => to be corrected/reviewed
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
|
||||
@ -39,7 +39,7 @@ class PairSpinMe : public Pair {
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_me(int, int, double *, double *, double *, double *);
|
||||
void compute_me(int, int, double fmi[3], double fmj[3], double spi[3], double spj[3]);
|
||||
|
||||
int me_flag; // me flag
|
||||
|
||||
@ -51,12 +51,7 @@ class PairSpinMe : public Pair {
|
||||
double hbar;
|
||||
|
||||
double **ME; // me coeff in eV
|
||||
double **v_mex, **v_mey, **v_mez;// me direction
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
double **v_mex, **v_mey, **v_mez; // me direction
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
@ -42,7 +42,7 @@ PairSpinSocDmi::PairSpinSocDmi(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
newton_pair_spin = 0; // no newton pair for now => to be corrected
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
@ -64,14 +64,6 @@ PairSpinSocDmi::~PairSpinSocDmi()
|
||||
memory->destroy(v_dmy);
|
||||
memory->destroy(v_dmz);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
@ -81,12 +73,13 @@ PairSpinSocDmi::~PairSpinSocDmi()
|
||||
void PairSpinSocDmi::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_dmi_2;
|
||||
double rsq,rd;
|
||||
double evdwl, ecoul;
|
||||
double xi[3], rij[3];
|
||||
double spi[3], spj[3];
|
||||
double fi[3], fj[3];
|
||||
double fmi[3], fmj[3];
|
||||
double cut_dmi_2, cut_spin_me_global2;
|
||||
double rsq, rd, inorm;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -113,9 +106,9 @@ void PairSpinSocDmi::compute(int eflag, int vflag)
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
@ -139,13 +132,11 @@ void PairSpinSocDmi::compute(int eflag, int vflag)
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rij[0] = x[j][0] - xi[0];
|
||||
rij[1] = x[j][1] - xi[1];
|
||||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
@ -197,7 +188,8 @@ void PairSpinSocDmi::compute(int eflag, int vflag)
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocDmi::compute_dmi(int i, int j, double *fmi, double *fmj, double *spi, double *spj)
|
||||
|
||||
void PairSpinSocDmi::compute_dmi(int i, int j, double fmi[3], double fmj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -238,14 +230,6 @@ void PairSpinSocDmi::allocate()
|
||||
memory->create(v_dmy,n+1,n+1,"pair:DM_vector_y");
|
||||
memory->create(v_dmz,n+1,n+1,"pair:DM_vector_z");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
@ -39,7 +39,7 @@ class PairSpinSocDmi : public Pair {
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_dmi(int, int, double *, double *, double *, double *);
|
||||
void compute_dmi(int, int, double fmi[3], double fmj[3], double spi[3], double spj[3]);
|
||||
|
||||
int dmi_flag; // dmi flag
|
||||
|
||||
@ -51,12 +51,7 @@ class PairSpinSocDmi : public Pair {
|
||||
double hbar;
|
||||
|
||||
double **DM; // dmi coeff in eV
|
||||
double **v_dmx, **v_dmy, **v_dmz;// dmi direction
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
double **v_dmx, **v_dmy, **v_dmz; // dmi direction
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
@ -42,7 +42,7 @@ PairSpinSocLandau::PairSpinSocLandau(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
newton_pair_spin = 0; // no newton pair for now => to be corrected
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
@ -58,20 +58,12 @@ PairSpinSocLandau::~PairSpinSocLandau()
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_soc_neel);
|
||||
memory->destroy(cut_soc_landau);
|
||||
memory->destroy(K1);
|
||||
memory->destroy(K1_mech);
|
||||
memory->destroy(K2);
|
||||
memory->destroy(K3);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
@ -81,12 +73,13 @@ PairSpinSocLandau::~PairSpinSocLandau()
|
||||
void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_soc_neel_2,cut_soc_global2;
|
||||
double rsq,rd;
|
||||
double evdwl, ecoul;
|
||||
double xi[3], rij[3];
|
||||
double spi[3], spj[3];
|
||||
double fi[3], fj[3];
|
||||
double fmi[3], fmj[3];
|
||||
double cut_soc_landau_2, cut_soc_global2;
|
||||
double rsq, rd, inorm;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -113,9 +106,9 @@ void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
@ -139,13 +132,11 @@ void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rij[0] = x[j][0] - xi[0];
|
||||
rij[1] = x[j][1] - xi[1];
|
||||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
@ -154,10 +145,10 @@ void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
jtype = type[j];
|
||||
|
||||
// compute mag. and mech. components of soc
|
||||
cut_soc_neel_2 = cut_soc_neel[itype][jtype]*cut_soc_neel[itype][jtype];
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
compute_soc_mech_neel(i,j,rsq,rij,fi,fj,spi,spj);
|
||||
cut_soc_landau_2 = cut_soc_landau[itype][jtype]*cut_soc_landau[itype][jtype];
|
||||
if (rsq <= cut_soc_landau_2) {
|
||||
compute_soc_landau(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
compute_soc_mech_landau(i,j,rsq,rij,fi,fj,spi,spj);
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
@ -167,7 +158,7 @@ void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
// if (newton_pair || j < nlocal) { => to be corrected
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
@ -178,7 +169,7 @@ void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
if (rsq <= cut_soc_landau_2) {
|
||||
evdwl -= spi[0]*fmi[0];
|
||||
evdwl -= spi[1]*fmi[1];
|
||||
evdwl -= spi[2]*fmi[2];
|
||||
@ -196,7 +187,8 @@ void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocLandau::compute_soc_neel(int i, int j, double rsq, double *rij, double *fmi, double *fmj, double *spi, double *spj)
|
||||
|
||||
void PairSpinSocLandau::compute_soc_landau(int i, int j, double rsq, double rij[3], double fmi[3], double fmj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -223,7 +215,8 @@ void PairSpinSocLandau::compute_soc_neel(int i, int j, double rsq, double *rij,
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocLandau::compute_soc_mech_neel(int i, int j, double rsq, double *rij, double *fi, double *fj, double *spi, double *spj)
|
||||
|
||||
void PairSpinSocLandau::compute_soc_mech_landau(int i, int j, double rsq, double rij[3], double fi[3], double fj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -279,21 +272,13 @@ void PairSpinSocLandau::allocate()
|
||||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
|
||||
memory->create(cut_soc_neel,n+1,n+1,"pair:cut_soc_neel");
|
||||
memory->create(K1,n+1,n+1,"pair:K1");
|
||||
memory->create(K1_mech,n+1,n+1,"pair:K1_mech");
|
||||
memory->create(K2,n+1,n+1,"pair:K2");
|
||||
memory->create(K3,n+1,n+1,"pair:K3");
|
||||
memory->create(cut_soc_landau,n+1,n+1,"pair/spin/soc/landau:cut_soc_landau");
|
||||
memory->create(K1,n+1,n+1,"pair/spin/soc/landau:K1");
|
||||
memory->create(K1_mech,n+1,n+1,"pair/spin/soc/landau:K1_mech");
|
||||
memory->create(K2,n+1,n+1,"pair/spin/soc/landau:K2");
|
||||
memory->create(K3,n+1,n+1,"pair/spin/soc/landau:K3");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
memory->create(cutsq,n+1,n+1,"pair/spin/soc/landau:cutsq");
|
||||
|
||||
}
|
||||
|
||||
@ -318,7 +303,7 @@ void PairSpinSocLandau::settings(int narg, char **arg)
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i+1; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j]) {
|
||||
cut_soc_neel[i][j] = cut_soc_global;
|
||||
cut_soc_landau[i][j] = cut_soc_global;
|
||||
}
|
||||
}
|
||||
|
||||
@ -335,7 +320,7 @@ void PairSpinSocLandau::coeff(int narg, char **arg)
|
||||
if (!allocated) allocate();
|
||||
|
||||
// set mech_flag to 1 if magneto-mech simulation
|
||||
//no longer correct: can be hybrid without magneto-mech
|
||||
//no longer correct: can be hybrid without magneto-mech
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
mech_flag = 0;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
@ -359,7 +344,7 @@ void PairSpinSocLandau::coeff(int narg, char **arg)
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_soc_neel[i][j] = rij;
|
||||
cut_soc_landau[i][j] = rij;
|
||||
K1[i][j] = k1/hbar;
|
||||
if (mech_flag) {
|
||||
K1_mech[i][j] = k1;
|
||||
@ -389,7 +374,7 @@ void PairSpinSocLandau::init_style()
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
// check this half/full request => to be corrected
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
@ -429,7 +414,7 @@ void PairSpinSocLandau::write_restart(FILE *fp)
|
||||
fwrite(&K1_mech[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_soc_neel[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_soc_landau[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -457,13 +442,13 @@ void PairSpinSocLandau::read_restart(FILE *fp)
|
||||
fread(&K1_mech[i][j],sizeof(double),1,fp);
|
||||
fread(&K2[i][j],sizeof(double),1,fp);
|
||||
fread(&K2[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_soc_neel[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_soc_landau[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&K1[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K1_mech[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_soc_neel[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_soc_landau[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -39,14 +39,14 @@ class PairSpinSocLandau : public Pair {
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_soc_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_soc_mech_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_soc_landau(int, int, double, double rij[3], double fmi[3], double fmj[3],double spi[3], double spj[3]);
|
||||
void compute_soc_mech_landau(int, int, double, double rij[3], double fi[3], double fj[3], double spi[3], double spj[3]);
|
||||
|
||||
int soc_neel_flag; // soc neel flag
|
||||
int mech_flag; // mech calc. flag
|
||||
|
||||
double cut_soc_global;
|
||||
double **cut_soc_neel; // cutoff distance exchange
|
||||
double **cut_soc_landau; // cutoff distance exchange
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
@ -55,11 +55,6 @@ class PairSpinSocLandau : public Pair {
|
||||
double **K1, **K1_mech; // exchange coeffs Kij
|
||||
double **K2, **K3; // K1 in eV, K2 adim, K3 in Ang
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
|
||||
@ -42,7 +42,7 @@ PairSpinSocNeel::PairSpinSocNeel(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
newton_pair_spin = 0; // no newton pair for now => to be corrected
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
@ -64,14 +64,6 @@ PairSpinSocNeel::~PairSpinSocNeel()
|
||||
memory->destroy(K2);
|
||||
memory->destroy(K3);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
@ -82,11 +74,12 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_soc_neel_2,cut_soc_global2;
|
||||
double rsq,rd;
|
||||
double xi[3], rij[3];
|
||||
double spi[3], spj[3];
|
||||
double fi[3], fj[3];
|
||||
double fmi[3], fmj[3];
|
||||
double cut_soc_neel_2, cut_soc_global2;
|
||||
double rsq, rd, inorm;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
@ -109,13 +102,13 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// pair spin computations
|
||||
// loop over neighbors of my atoms
|
||||
// loop over atoms and their neighbors
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
@ -123,6 +116,7 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
@ -139,13 +133,11 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rij[0] = x[j][0] - xi[0];
|
||||
rij[1] = x[j][1] - xi[1];
|
||||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
@ -153,7 +145,8 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// compute mag. and mech. components of soc
|
||||
// compute magnetic and mechanical components of soc_neel
|
||||
|
||||
cut_soc_neel_2 = cut_soc_neel[itype][jtype]*cut_soc_neel[itype][jtype];
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
@ -167,7 +160,7 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
// if (newton_pair || j < nlocal) { // => to be corrected
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
@ -196,7 +189,8 @@ void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocNeel::compute_soc_neel(int i, int j, double rsq, double *rij, double *fmi, double *fmj, double *spi, double *spj)
|
||||
|
||||
void PairSpinSocNeel::compute_soc_neel(int i, int j, double rsq, double rij[3], double fmi[3], double fmj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -212,18 +206,19 @@ void PairSpinSocNeel::compute_soc_neel(int i, int j, double rsq, double *rij, do
|
||||
scalar = rij[0]*spj[0]+rij[1]*spj[1]+rij[2]*spj[2];
|
||||
Kij_3 = Kij/3.0;
|
||||
|
||||
fmi[0] -= Kij*scalar*rij[0]-Kij_3*spj[0];
|
||||
fmi[1] -= Kij*scalar*rij[1]-Kij_3*spj[1];
|
||||
fmi[2] -= Kij*scalar*rij[2]-Kij_3*spj[2];
|
||||
fmi[0] -= Kij*scalar*rij[0] - Kij_3*spj[0];
|
||||
fmi[1] -= Kij*scalar*rij[1] - Kij_3*spj[1];
|
||||
fmi[2] -= Kij*scalar*rij[2] - Kij_3*spj[2];
|
||||
|
||||
fmj[0] += Kij*scalar*rij[0]+Kij_3*spi[0];
|
||||
fmj[1] += Kij*scalar*rij[1]+Kij_3*spi[1];
|
||||
fmj[2] += Kij*scalar*rij[2]+Kij_3*spi[2];
|
||||
fmj[0] += Kij*scalar*rij[0] + Kij_3*spi[0];
|
||||
fmj[1] += Kij*scalar*rij[1] + Kij_3*spi[1];
|
||||
fmj[2] += Kij*scalar*rij[2] + Kij_3*spi[2];
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocNeel::compute_soc_mech_neel(int i, int j, double rsq, double *rij, double *fi, double *fj, double *spi, double *spj)
|
||||
|
||||
void PairSpinSocNeel::compute_soc_mech_neel(int i, int j, double rsq, double rij[3], double fi[3], double fj[3], double spi[3], double spj[3])
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
@ -285,14 +280,6 @@ void PairSpinSocNeel::allocate()
|
||||
memory->create(K2,n+1,n+1,"pair:K2");
|
||||
memory->create(K3,n+1,n+1,"pair:K3");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
@ -335,7 +322,7 @@ void PairSpinSocNeel::coeff(int narg, char **arg)
|
||||
if (!allocated) allocate();
|
||||
|
||||
// set mech_flag to 1 if magneto-mech simulation
|
||||
//no longer correct: can be hybrid without magneto-mech
|
||||
//no longer correct: can be hybrid without magneto-mech => needs review/correction
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
mech_flag = 0;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
@ -377,7 +364,6 @@ void PairSpinSocNeel::coeff(int narg, char **arg)
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
@ -389,7 +375,7 @@ void PairSpinSocNeel::init_style()
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
// check this half/full request => to be verified
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
@ -469,7 +455,6 @@ void PairSpinSocNeel::read_restart(FILE *fp)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
@ -39,8 +39,8 @@ class PairSpinSocNeel : public Pair {
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_soc_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_soc_mech_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_soc_neel(int, int, double, double rij[3], double fmi[3], double fmj[3],double spi[3], double spj[3]);
|
||||
void compute_soc_mech_neel(int, int, double, double rij[3], double fi[3], double fj[3],double spi[3], double spj[3]);
|
||||
|
||||
int soc_neel_flag; // soc neel flag
|
||||
int mech_flag; // mech calc. flag
|
||||
@ -55,11 +55,6 @@ class PairSpinSocNeel : public Pair {
|
||||
double **K1, **K1_mech; // exchange coeffs Kij
|
||||
double **K2, **K3; // K1 in eV, K2 adim, K3 in Ang
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
|
||||
@ -98,7 +98,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
|
||||
rho = drho = e = de = cv = NULL;
|
||||
vest = NULL;
|
||||
|
||||
// USER-SPIN
|
||||
// SPIN package
|
||||
|
||||
mumag = NULL;
|
||||
sp = fm = NULL;
|
||||
|
||||
@ -172,7 +173,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
|
||||
radius_flag = rmass_flag = 0;
|
||||
ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
|
||||
|
||||
//Magnetic flags
|
||||
// magnetic flags
|
||||
|
||||
sp_flag = mumag_flag = 0;
|
||||
|
||||
vfrac_flag = 0;
|
||||
@ -428,7 +430,8 @@ 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
|
||||
// magnetic flags
|
||||
|
||||
sp_flag = mumag_flag = 0;
|
||||
|
||||
vfrac_flag = 0;
|
||||
|
||||
@ -62,6 +62,7 @@ class Atom : protected Pointers {
|
||||
int *ellipsoid,*line,*tri,*body;
|
||||
|
||||
// SPIN package
|
||||
|
||||
double *mumag, **sp;
|
||||
double **fm;
|
||||
|
||||
|
||||
Reference in New Issue
Block a user