diff --git a/examples/SPIN/in.spin.bfo b/examples/SPIN/in.spin.bfo new file mode 100644 index 0000000000..aee3454765 --- /dev/null +++ b/examples/SPIN/in.spin.bfo @@ -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 + diff --git a/examples/SPIN/in.spin.cobalt b/examples/SPIN/in.spin.cobalt new file mode 100644 index 0000000000..8147aad5fc --- /dev/null +++ b/examples/SPIN/in.spin.cobalt @@ -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 + diff --git a/examples/SPIN/in.spin.kagome b/examples/SPIN/in.spin.kagome new file mode 100644 index 0000000000..c51c35ff73 --- /dev/null +++ b/examples/SPIN/in.spin.kagome @@ -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 + diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index 214d7752b6..8edf8211ca 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -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"); diff --git a/src/SPIN/atom_vec_spin.h b/src/SPIN/atom_vec_spin.h index aedc3efb3e..c0f245ba27 100644 --- a/src/SPIN/atom_vec_spin.h +++ b/src/SPIN/atom_vec_spin.h @@ -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); @@ -66,8 +67,8 @@ class AtomVecSpin : public AtomVec { tagint *tag; int *type,*mask; imageint *image; - double **x,**v,**f; // lattice quantities - double *mumag,**sp, **fm; // spin quantities + double **x,**v,**f; // lattice quantities + double *mumag,**sp,**fm; // spin quantities }; diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index ce60cf373f..a7069fa808 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -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"); } diff --git a/src/SPIN/compute_spin.h b/src/SPIN/compute_spin.h index a20b956b01..59f0ce2876 100644 --- a/src/SPIN/compute_spin.h +++ b/src/SPIN/compute_spin.h @@ -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(); diff --git a/src/SPIN/fix_force_spin.cpp b/src/SPIN/fix_force_spin.cpp index b0e38989c9..921a9964bb 100644 --- a/src/SPIN/fix_force_spin.cpp +++ b/src/SPIN/fix_force_spin.cpp @@ -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; } @@ -124,19 +122,18 @@ int FixForceSpin::setmask() void FixForceSpin::init() { - const double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) - const double mub = 5.78901e-5; // in eV/T - const double gyro = mub/hbar; // in rad.THz/T + const double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) + const double mub = 5.78901e-5; // in eV/T + const double gyro = mub/hbar; // in rad.THz/T - H_field *= gyro; // in rad.THz - Ka /= hbar; // in rad.THz + H_field *= gyro; // in rad.THz + Ka /= hbar; // in rad.THz if (strstr(update->integrate_style,"respa")) { ilevel_respa = ((Respa *) update->integrate)->nlevels-1; if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); } - // check variables if (magstr) { magvar = input->variable->find(magstr); if (magvar < 0) @@ -149,11 +146,9 @@ 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"); - } /* ---------------------------------------------------------------------- */ @@ -174,15 +169,17 @@ 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); - set_magneticforce(); // update mag. field if time-dep. + set_magneticforce(); // update mag. field if time-dep. } double **sp = atom->sp; double *mumag = atom->mumag; - double **fm = atom->fm; + 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; diff --git a/src/SPIN/fix_force_spin.h b/src/SPIN/fix_force_spin.h index 7da3ece415..1f2b36d8e7 100644 --- a/src/SPIN/fix_force_spin.h +++ b/src/SPIN/fix_force_spin.h @@ -38,11 +38,11 @@ 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 + int style; // style of the magnetic force double degree2rad; double hbar; @@ -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 + 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 + double Kax, Kay, Kaz; // temp. force variables - // temp. spin variables - double *spi, *fmi; - void set_magneticforce(); }; diff --git a/src/SPIN/fix_integration_spin.cpp b/src/SPIN/fix_integration_spin.cpp index a790aa7ea6..4e7bc961c4 100644 --- a/src/SPIN/fix_integration_spin.cpp +++ b/src/SPIN/fix_integration_spin.cpp @@ -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(sec); + memory->destroy(k); + memory->destroy(adv_list); + 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; @@ -190,11 +170,19 @@ void FixIntegrationSpin::init() if (locklangevinspin->tdamp_flag == 1) tdamp_flag = 1; if (locklangevinspin->temp_flag == 1) temp_flag = 1; } + - + 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,26 +217,24 @@ 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 + for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal comm->forward_comm(); k[j] = 0; for (int i = 0; i < nlocal; i++) { @@ -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 + } 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 - for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal + 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,87 +284,42 @@ void FixIntegrationSpin::initial_integrate(int vflag) } -#if defined LIST - // update half s for all particles + // 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 + 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 + 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 + } 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 - for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal + 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,19 +384,19 @@ 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; - if (exch_flag) { // exchange + if (exch_flag) { // exchange temp_cut = lockpairspinexchange->cut_spin_exchange[itype][jtype]; cut_2 = temp_cut*temp_cut; if (rsq <= cut_2) { @@ -510,38 +404,36 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii) } } - if (soc_flag) { // soc + if (soc_flag) { // soc 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); } } } - if (magforce_flag) { // mag. forces - if (zeeman_flag) { // zeeman + if (magforce_flag) { // mag. forces + if (zeeman_flag) { // zeeman lockforcespin->compute_zeeman(i,fmi); } - if (aniso_flag) { // aniso - spi[0] = sp[i][0]; - spi[1] = sp[i][1]; - spi[2] = sp[i][2]; + if (aniso_flag) { // aniso lockforcespin->compute_anisotropy(i,spi,fmi); } } - if (maglangevin_flag) { // mag. langevin - if (tdamp_flag) { // trans. damping + if (maglangevin_flag) { // mag. langevin + 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(rix); seci[1] = static_cast(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]; diff --git a/src/SPIN/fix_integration_spin.h b/src/SPIN/fix_integration_spin.h index 1156cc576f..dfeef599a3 100644 --- a/src/SPIN/fix_integration_spin.h +++ b/src/SPIN/fix_integration_spin.h @@ -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 - class PairHybrid *lockhybrid; - class PairSpin *lockpairspin; + // pointers to magnetic interaction classes + + class PairHybrid *lockhybrid; 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; }; diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index 1af32ff5cb..b017e20ae3 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -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); } @@ -135,8 +135,8 @@ void FixLangevinSpin::init() gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t)); dts = update->dt; - double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) - double kb = force->boltz; // eV/K + double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) + double kb = force->boltz; // eV/K D = (MY_2PI*alpha_t*gil_factor*kb*temp); D /= (hbar*dts); sigma = sqrt(2.0*D); @@ -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; diff --git a/src/SPIN/fix_langevin_spin.h b/src/SPIN/fix_langevin_spin.h index 3690037225..d95f5c08c7 100644 --- a/src/SPIN/fix_langevin_spin.h +++ b/src/SPIN/fix_langevin_spin.h @@ -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 - int tdamp_flag, ldamp_flag, temp_flag; // damping and temperature flags + 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 temp; // spin bath temperature - double D,sigma; // bath intensity var. - - double gil_factor; // Gilbert's prefactor + 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 char *id_temp; class Compute *temperature; int nlevels_respa; - //class RanMars *random; class RanPark *random; int seed; diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp deleted file mode 100755 index e4e8e2e15a..0000000000 --- a/src/SPIN/pair_spin.cpp +++ /dev/null @@ -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 -#include -#include - -#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); -} diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h deleted file mode 100755 index f73be7d466..0000000000 --- a/src/SPIN/pair_spin.h +++ /dev/null @@ -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. - -*/ diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index 208a82ea30..9080fce51d 100755 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -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; @@ -64,14 +64,6 @@ PairSpinExchange::~PairSpinExchange() memory->destroy(J1_mech); 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); diff --git a/src/SPIN/pair_spin_exchange.h b/src/SPIN/pair_spin_exchange.h index 14cdf7a959..38cd03570f 100755 --- a/src/SPIN/pair_spin_exchange.h +++ b/src/SPIN/pair_spin_exchange.h @@ -39,26 +39,21 @@ 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_global; // global exchange cutoff + double **cut_spin_exchange; // cutoff distance per exchange 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 *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 **J1_mag, **J1_mech; // exchange coeffs Jij + double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang void allocate(); }; diff --git a/src/SPIN/pair_spin_me.cpp b/src/SPIN/pair_spin_me.cpp index b48abe93e6..3c6f4fa7ed 100755 --- a/src/SPIN/pair_spin_me.cpp +++ b/src/SPIN/pair_spin_me.cpp @@ -63,14 +63,6 @@ PairSpinMe::~PairSpinMe() 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(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); diff --git a/src/SPIN/pair_spin_me.h b/src/SPIN/pair_spin_me.h index fb216a61ea..31e0d8b560 100755 --- a/src/SPIN/pair_spin_me.h +++ b/src/SPIN/pair_spin_me.h @@ -39,24 +39,19 @@ 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 + int me_flag; // me flag - double cut_spin_me_global; // global me cutoff - double **cut_spin_me; // me cutoff distance + double cut_spin_me_global; // global me cutoff + double **cut_spin_me; // me cutoff distance protected: int newton_pair_spin; 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 **ME; // me coeff in eV + double **v_mex, **v_mey, **v_mez; // me direction void allocate(); }; diff --git a/src/SPIN/pair_spin_soc_dmi.cpp b/src/SPIN/pair_spin_soc_dmi.cpp index eb4f957c8e..341f230fd3 100755 --- a/src/SPIN/pair_spin_soc_dmi.cpp +++ b/src/SPIN/pair_spin_soc_dmi.cpp @@ -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; @@ -63,14 +63,6 @@ PairSpinSocDmi::~PairSpinSocDmi() memory->destroy(v_dmx); 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"); } diff --git a/src/SPIN/pair_spin_soc_dmi.h b/src/SPIN/pair_spin_soc_dmi.h index a28d1772d7..627001135e 100755 --- a/src/SPIN/pair_spin_soc_dmi.h +++ b/src/SPIN/pair_spin_soc_dmi.h @@ -39,24 +39,19 @@ 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 + int dmi_flag; // dmi flag - double cut_spin_dmi_global; // short range pair cutoff - double **cut_spin_dmi; // cutoff distance dmi + double cut_spin_dmi_global; // short range pair cutoff + double **cut_spin_dmi; // cutoff distance dmi protected: int newton_pair_spin; 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 **DM; // dmi coeff in eV + double **v_dmx, **v_dmy, **v_dmz; // dmi direction void allocate(); }; diff --git a/src/SPIN/pair_spin_soc_landau.cpp b/src/SPIN/pair_spin_soc_landau.cpp index 13cd83868b..2beebb47f5 100755 --- a/src/SPIN/pair_spin_soc_landau.cpp +++ b/src/SPIN/pair_spin_soc_landau.cpp @@ -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,19 +58,11 @@ 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); } } } diff --git a/src/SPIN/pair_spin_soc_landau.h b/src/SPIN/pair_spin_soc_landau.h index c23739c071..491cb9a581 100755 --- a/src/SPIN/pair_spin_soc_landau.h +++ b/src/SPIN/pair_spin_soc_landau.h @@ -39,26 +39,21 @@ 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 + 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; double hbar; - 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 + double **K1, **K1_mech; // exchange coeffs Kij + double **K2, **K3; // K1 in eV, K2 adim, K3 in Ang void allocate(); }; diff --git a/src/SPIN/pair_spin_soc_neel.cpp b/src/SPIN/pair_spin_soc_neel.cpp index e8cd96fece..94ac35848b 100755 --- a/src/SPIN/pair_spin_soc_neel.cpp +++ b/src/SPIN/pair_spin_soc_neel.cpp @@ -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; @@ -63,14 +63,6 @@ PairSpinSocNeel::~PairSpinSocNeel() 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); } @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/SPIN/pair_spin_soc_neel.h b/src/SPIN/pair_spin_soc_neel.h index ce9c76eb49..4584ddb9eb 100755 --- a/src/SPIN/pair_spin_soc_neel.h +++ b/src/SPIN/pair_spin_soc_neel.h @@ -39,26 +39,21 @@ 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 + 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_neel; // cutoff distance exchange protected: int newton_pair_spin; double hbar; - 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 + double **K1, **K1_mech; // exchange coeffs Kij + double **K2, **K3; // K1 in eV, K2 adim, K3 in Ang void allocate(); }; diff --git a/src/atom.cpp b/src/atom.cpp index a847719e7e..86ef309f09 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -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; diff --git a/src/atom.h b/src/atom.h index 8e7f9811ac..285f067e24 100644 --- a/src/atom.h +++ b/src/atom.h @@ -62,6 +62,7 @@ class Atom : protected Pointers { int *ellipsoid,*line,*tri,*body; // SPIN package + double *mumag, **sp; double **fm;