From d53def58538a1835c14428b13ea4a4064dd12456 Mon Sep 17 00:00:00 2001 From: julient31 Date: Fri, 26 May 2017 08:40:57 -0600 Subject: [PATCH] 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 --- in.cobalt | 23 +++-- src/SPIN/atom_vec_spin.cpp | 18 ++-- src/SPIN/compute_spin.cpp | 2 +- src/SPIN/fix_nve_spin.cpp | 201 +++++++++++++++++++++++++++++-------- src/SPIN/fix_nve_spin.h | 7 +- src/SPIN/pair_spin.cpp | 99 +++++++++++------- src/SPIN/pair_spin.h | 6 ++ 7 files changed, 256 insertions(+), 100 deletions(-) diff --git a/in.cobalt b/in.cobalt index 06b5508fde..1e58798c16 100644 --- a/in.cobalt +++ b/in.cobalt @@ -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 diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp index bda03f9b27..3807ac0e56 100644 --- a/src/SPIN/atom_vec_spin.cpp +++ b/src/SPIN/atom_vec_spin.cpp @@ -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); } diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index 9a0ba20496..bd1b09fb44 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -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; } diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index f156f76c8a..20aaf4afef 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -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;inlocal; + 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]; } } - } diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 55fffa7fb0..82e9822bf1 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -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; }; } diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index a2c75cbc90..d684a7800d 100755 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -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 diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index abdc1978c1..fd5822c163 100755 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -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