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

Field compute error apparently corrected. The issue was related to the reverse communication. 

To do:
- Remove all checks/prints used to debug
- Check all the flag set in the atom_vec_spin creator (very important for the reverse comm)
- Code DMI/ME interactions
- Start to work on parallel implementation of the integration
This commit is contained in:
julient31 2017-05-26 08:40:57 -06:00
parent af45d55b3f
commit d53def5853
7 changed files with 256 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
in.cobalt
View File

@ -18,17 +18,20 @@ 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 constant of fcc Cobalt
lattice fcc 3.54
#lattice sc 3.54
#lattice sc 2.50
#lattice bcc 3.54
#Defining a geometric region of space. Sets ID(user's choice), style(block, sphere, ...), then, args depends on the style chosen
#(for block, one has x0, xf, y0, yf, z0, zf, in distance units)
region box block 0.0 3.0 0.0 3.0 0.0 3.0
region box block 0.0 2.0 0.0 2.0 0.0 2.0
#Creating a simulation box based on the specified region. Entries: number of atom types and box ref.
create_box 1 box
@ -60,20 +63,24 @@ pair_coeff * * 0.0446928 0.003496 1.4885
#Fix Langevin spins (merging damping and temperature)
#Defines a cutof distance for the neighbors.
#neighbor 0.5 bin
#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
fix 1 all force/spin zeeman 1.0 0.0 0.0 1.0
fix 1 all force/spin zeeman 10.0 0.0 0.0 1.0
#fix 1 all force/spin zeeman 0.0 0.0 0.0 1.0
#fix 1 all force/spin anisotropy 0.001 0.0 0.0 1.0
#Fix Langevin spins (merging damping and temperature)
fix 2 all langevin/spin 0.0 0.1 0.0 21
#fix 2 all langevin/spin 0.0 0.0 0.1 21
#fix 2 all langevin/spin 0.0 0.0 0.0 21
#Magnetic integration fix
fix 3 all nve/spin
#compute total magnetization and magnetic energy
#Ite Time Mx My Mz |M| Em Tm
compute mag all compute/spin
fix outmag all ave/time 1 1 50 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag.dat
@ -81,7 +88,7 @@ fix outmag all ave/time 1 1 50 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_ma
#Entries: ID(user assigned), group-ID(group of atoms to peform the sim on), style(temp, pe, ...), args
#Setting the timestep for the simulation
timestep 0.000005
timestep 0.00005
##################
#######run########
@ -91,6 +98,6 @@ timestep 0.000005
dump 1 all custom 50 dump_spin.lammpstrj type x y z spx spy spz
#Running the simulations for N timesteps
run 200000
#run 2
run 30000
#run 1

18
src/SPIN/atom_vec_spin.cpp
View File

@ -30,15 +30,17 @@ using namespace LAMMPS_NS;
AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
{
molecular = 0;
mass_type = 1;
mass_type = 1; //check why
comm_x_only = 0;
comm_f_only = 1;
size_forward = 6;
size_reverse = 3;
//comm_x_only = 0;
comm_x_only = 1;
//comm_f_only = 1;
comm_f_only = 0;
size_forward = 7;
size_reverse = 6;
size_border = 11;
size_velocity = 3;
size_data_atom = 9;
size_data_atom = 9; //to check later
size_data_vel = 4;
xcol_data = 4;
@ -318,7 +320,6 @@ int AtomVecSpin::unpack_comm_hybrid(int n, int first, double *buf)
int AtomVecSpin::pack_reverse(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
@ -337,7 +338,6 @@ int AtomVecSpin::pack_reverse(int n, int first, double *buf)
void AtomVecSpin::unpack_reverse(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
@ -941,9 +941,9 @@ bigint AtomVecSpin::memory_usage()
return bytes;
}
//Test force clear in spin
void AtomVecSpin::force_clear(int n, size_t nbytes)
{
memset(&atom->f[0][0],0,3*nbytes);
memset(&atom->fm[0][0],0,3*nbytes);
}

2
src/SPIN/compute_spin.cpp
View File

@ -125,7 +125,7 @@ void ComputeSpin::compute_vector()
vector[2] = magtot[1];
vector[3] = magtot[2];
vector[4] = magtot[3];
vector[5] = magenergytot;
vector[5] = -0.5*magenergytot*hbar;
vector[6] = spintemperature;
}

201
src/SPIN/fix_nve_spin.cpp
View File

@ -26,6 +26,9 @@
#include "math_const.h"
#include "modify.h"
#include "pair.h"
#include "timer.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
@ -105,7 +108,7 @@ void FixNVESpin::initial_integrate(int vflag)
v[i][2] += dtfm * f[i][2];
}
}
// update half x for all particles
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -115,55 +118,168 @@ void FixNVESpin::initial_integrate(int vflag)
}
}
size_t nbytes;
nbytes = sizeof(double) * nlocal;
int eflag = 3;
// update sp for all particles
if (extra == SPIN) {
// Advance spins
//See Omelyan et al., PRL 86, 2001 and P.W. Ma et al, PRB 83, 2011
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
cp[0] = cp[1] = cp[2] = 0.0;
g[0] = g[1] = g[2] = 0.0;
fm2 = (fm[i][0]*fm[i][0])+(fm[i][1]*fm[i][1])+(fm[i][2]*fm[i][2]);
fmsq = sqrt(fm2);
energy = (sp[i][0]*fm[i][0])+(sp[i][1]*fm[i][1])+(sp[i][2]*fm[i][2]);
cp[0] = fm[i][1]*sp[i][2]-fm[i][2]*sp[i][1];
cp[1] = fm[i][2]*sp[i][0]-fm[i][0]*sp[i][2];
cp[2] = fm[i][0]*sp[i][1]-fm[i][1]*sp[i][0];
g[0] = sp[i][0]+cp[0]*dts;
g[1] = sp[i][1]+cp[1]*dts;
g[2] = sp[i][2]+cp[2]*dts;
g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts*dts;
g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts*dts;
g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts*dts;
g[0] /= (1+0.25*fm2*dts*dts);
g[1] /= (1+0.25*fm2*dts*dts);
g[2] /= (1+0.25*fm2*dts*dts);
sp[i][0] = g[0];
sp[i][1] = g[1];
sp[i][2] = g[2];
//Renormalization (may not be necessary)
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = 1.0/sqrt(msq);
sp[i][0] *= scale;
sp[i][1] *= scale;
sp[i][2] *= scale;
//printf("test fix integ. 1;i=%d, fx=%g, fy=%g, fz=%g \n",i,fm[i][0],fm[i][1],fm[i][2]);
//printf("test fix integ.; i=%d, sx=%g, sy=%g, sz=%g, norm=%g \n",i,sp[i][0],sp[i][1],sp[i][2],scale);
}
}
#define CONC
#if defined CONC
for (int i = 0; i < nlocal; i++){
AdvanceSingleSpin(i,dts,sp,fm);
}
#endif
//#define SEQ
#if defined SEQ
//advance N-1 spins to half
for (int i = 0; i < nlocal-1; i++){
//Recomp field
atom->avec->force_clear(0,nbytes);
timer->stamp();
modify->pre_force(vflag);
timer->stamp(Timer::PAIR);
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
comm->reverse_comm();
timer->stamp(Timer::COMM);
modify->post_force(vflag);
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
//advance N spin
//Recomp field
atom->avec->force_clear(0,nbytes);
timer->stamp();
modify->pre_force(vflag);
timer->stamp(Timer::PAIR);
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
comm->reverse_comm();
timer->stamp(Timer::COMM);
modify->post_force(vflag);
AdvanceSingleSpin(nlocal-1,dts,sp,fm);
//advance N-1 spins to half
for (int i = nlocal-2; i >= 0; i--){
//Recomp field
atom->avec->force_clear(0,nbytes);
timer->stamp();
modify->pre_force(vflag);
timer->stamp(Timer::PAIR);
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
comm->reverse_comm();
timer->stamp(Timer::COMM);
modify->post_force(vflag);
AdvanceSingleSpin(i,0.5*dts,sp,fm);
}
#endif
}
#define FORCE_PRINT
#if defined FORCE_PRINT
FILE* file_force=NULL;
file_force=fopen("spin_force_Lammps.dat","a");
fprintf(file_force,"---------------------------------- \n");
for(int i=0;i<nlocal;i++){
fprintf(file_force,"%d %lf %lf %lf \n",i,fm[i][0],fm[i][1],fm[i][2]);
}
if (file_force!=NULL) fclose(file_force);
#endif
//printf("test fix integ. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm)
{
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
int *type = atom->type;
int *mask = atom->mask;
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy;
double cp[3],g[3];
cp[0] = cp[1] = cp[2] = 0.0;
g[0] = g[1] = g[2] = 0.0;
fm2 = (fm[i][0]*fm[i][0])+(fm[i][1]*fm[i][1])+(fm[i][2]*fm[i][2]);
fmsq = sqrt(fm2);
energy = (sp[i][0]*fm[i][0])+(sp[i][1]*fm[i][1])+(sp[i][2]*fm[i][2]);
cp[0] = fm[i][1]*sp[i][2]-fm[i][2]*sp[i][1];
cp[1] = fm[i][2]*sp[i][0]-fm[i][0]*sp[i][2];
cp[2] = fm[i][0]*sp[i][1]-fm[i][1]*sp[i][0];
//#define ALG_MA //Ma algo
#if defined ALG_MA
double A[3];
A[0]= A[1] = A[2] = 0.0;
double xi=fmsq*dts;
double zeta=0.0; // this is because omega contains already the transverse damping
double chi=energy/fmsq;
double expo=exp(2.0*zeta);
double K1=1.0+expo+chi*(1.0-expo);
double K=2.0*exp(zeta)/K1;
double Ktrigo=1.0+0.25*xi*xi;
double cosinus=(1.0-0.25*xi*xi)/Ktrigo; //cos(xi)
double sinus=xi/Ktrigo; //sin(xi)
A[0]=K*cosinus;
A[1]=K*sinus;
A[2]=(1.0-expo+chi*(1.0+expo-2.0*exp(zeta)*cosinus))/K1;
g[0] = A[0]*sp[i][0]+A[1]*cp[0]/fmsq+A[2]*fm[i][0]/fmsq;
g[1] = A[0]*sp[i][1]+A[1]*cp[1]/fmsq+A[2]*fm[i][0]/fmsq;
g[2] = A[0]*sp[i][2]+A[1]*cp[2]/fmsq+A[2]*fm[i][0]/fmsq;
#endif
#define ALG_OM //Omelyan algo
#if defined ALG_OM
g[0] = sp[i][0]+cp[0]*dts;
g[1] = sp[i][1]+cp[1]*dts;
g[2] = sp[i][2]+cp[2]*dts;
g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts*dts;
g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts*dts;
g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts*dts;
g[0] /= (1+0.25*fm2*dts*dts);
g[1] /= (1+0.25*fm2*dts*dts);
g[2] /= (1+0.25*fm2*dts*dts);
#endif
sp[i][0] = g[0];
sp[i][1] = g[1];
sp[i][2] = g[2];
//Renormalization (may not be necessary)
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = 1.0/sqrt(msq);
sp[i][0] *= scale;
sp[i][1] *= scale;
sp[i][2] *= scale;
}
/* ---------------------------------------------------------------------- */
void FixNVESpin::final_integrate()
@ -180,7 +296,7 @@ void FixNVESpin::final_integrate()
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
int *type = atom->type;
int *mask = atom->mask;
// update half x for all particles
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -200,5 +316,4 @@ void FixNVESpin::final_integrate()
v[i][2] += dtfm * f[i][2];
}
}
}

7
src/SPIN/fix_nve_spin.h
View File

@ -32,15 +32,16 @@ class FixNVESpin : public FixNVE {
virtual ~FixNVESpin() {}
void init();
virtual void initial_integrate(int);
void AdvanceSingleSpin(int, double, double **, double **);
virtual void final_integrate();
protected:
int extra;
double dts;
double alpha_t;
//double alpha_t;
private:
class FixSpinDamping *lockspindamping;
//private:
//class FixSpinDamping *lockspindamping;
};
}

99
src/SPIN/pair_spin.cpp
View File

@ -78,9 +78,28 @@ void PairSpin::compute(int eflag, int vflag)
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
//#define GHOST_TAG
tagint *tag = atom->tag;
//#define GH_PR
#if defined GH_PR
FILE* file=NULL;
file=fopen("spin_ghosts_lammps.lammpstrj","w");
fprintf(file,"ITEM: TIMESTEP\n");
fprintf(file,"0.0 \n");
fprintf(file,"ITEM: NUMBER OF ATOMS\n");
fprintf(file,"n \n");
fprintf(file,"ITEM: BOX BOUNDS\n");
for(int d=0; d<3; d++) fprintf(file,"%lf %lf\n",-10.0,10.0);
fprintf(file,"ITEM: ATOMS type x y z vx vy vz\n");
#endif
// Pair spin computations
// Loop over neighbors of my atoms
// Loop over neighbors of my itoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
@ -90,13 +109,6 @@ void PairSpin::compute(int eflag, int vflag)
jlist = firstneigh[i];
jnum = numneigh[i];
//printf(":::::::Loop for atom numb. %d, jnum=%d ::::::: \n",i,jnum);
//printf("Test print real atom: i=%d, sx=%g, sy=%g, sz=%g \n",i,sp[i][0],sp[i][1],sp[i][2]);
//printf("Test print real atom: i=%d, rx=%g, ry=%g, rz=%g \n",i,x[i][0],x[i][1],x[i][2]);
int testcount=0;
//Exchange interaction
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
@ -120,7 +132,6 @@ void PairSpin::compute(int eflag, int vflag)
Jex = 4.0*J_1[itype][jtype]*ra;
Jex *= (1.0-J_2[itype][jtype]*ra);
Jex *= exp(-ra);
Jex *= mumag[ii]*mumag[jj];
fmix = Jex*sp[j][0];
fmiy = Jex*sp[j][1];
@ -129,46 +140,62 @@ void PairSpin::compute(int eflag, int vflag)
fmjx = Jex*sp[i][0];
fmjy = Jex*sp[i][1];
fmjz = Jex*sp[i][2];
#if defined GH_PR
fprintf(file,"%d %lf %lf %lf %lf %lf %lf\n",j,x[j][0],x[j][1],x[j][2],sp[j][0],sp[j][1],sp[j][2]);
#endif
//printf("Neighb pair: i=%d, j=%d \n",i,j);
//printf("Test print ghost/real neib atom: i=%d, j=%d, sx=%g, sy=%g, sz=%g \n",i,j,sp[j][0],sp[j][1],sp[j][2]);
//printf("Test g/r neib pair: i=%d, j=%d, rx=%g, ry=%g, rz=%g \n",i,j,x[j][0],x[j][1],x[j][2]);
//printf("Atom i: %d of type %d, Atom j: %d of type %d, \n",i,itype,j,jtype);
//printf("Exchange pair (%d,%d), Jij=%g, rij=%g \n",i,j,Jex,rd);
testcount++;
#if defined GHOST_TAG
printf("Test print ghost/real neib atom: i=%d, j=%d, tag=%d, fx=%g, fy=%g, fz=%g \n",i,j,tag[j],fmjx,fmjy,fmjz);
#endif
}
//printf("Test print ghost/real atom: j=%d, sx=%g, sy=%g, sz=%g \n",j,sp[j][0],sp[j][1],sp[j][2]);
//printf("Test print ghost/real atom: j=%d, rx=%g, ry=%g, rz=%g \n",j,x[j][0],x[j][1],x[j][2]);
fm[i][0] += fmix;
fm[i][1] += fmiy;
fm[i][2] += fmiz;
if (newton_pair || j < nlocal) {
//#define INDEX
#if defined INDEX
int index = atom->map(tag[j]);
fm[index][0] += fmjx;
fm[index][1] += fmjy;
fm[index][2] += fmjz;
#else
if (newton_pair || j < nlocal) {
fm[j][0] += fmjx;
fm[j][1] += fmjy;
fm[j][2] += fmjz;
}
#endif
//printf("Val fm %d: [%g,%g,%g] \n",i,fm[i][0],fm[i][1],fm[i][2]);
//printf("Val fm %d: [%g,%g,%g] \n",j,fm[j][0],fm[j][1],fm[j][2]);
}
// printf("Test count %d \n",testcount);
}
}
//printf("New pair val: %d \n",newton_pair);
//printf("vals exchange: Jx=%g, Jy=%g, Jz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//printf("test exchange. 1;i=0, fx=%g, fy=%g, fz=%g \n",fm[0][0],fm[0][1],fm[0][2]);
//printf("::::::::::::::::::::::::: End loop ::::::::::::::::::::::::\n");
//Test print atoms + ghosts
#if defined GH_PR
if (file!=NULL) fclose(file);
#endif
// printf("::::::::::::::::::::::::: End loop ::::::::::::::::::::::::\n");
if (vflag_fdotr) virial_fdotr_compute();
//#define TAG_NALL
#if defined TAG_NALL
int nall = nlocal + atom->nghost;
for (i = 0; i < nall; i++) {
int num_neig = numneigh[i];
int itag = tag[i];
printf("atom %d has %d neighbs and tag numb %d \n",i,num_neig,itag);
}
#endif
#if defined TRANS_FORCE
transferfm(fm);
#endif
}
/* ----------------------------------------------------------------------
allocate all arrays

6
src/SPIN/pair_spin.h
View File

@ -37,6 +37,12 @@ class PairSpin : public Pair {
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
//Test transf. force
//#define TRANS_FORCE
#if defined TRANS_FORCE
void transferfm(double **);
#endif
protected:
double cut_spin_exchange_global, cut_spin_dipolar_global; //Global cutting distance