diff --git a/examples/SPIN/gneb/README b/examples/SPIN/gneb/README index 7d388d898e..6f5db6d388 100644 --- a/examples/SPIN/gneb/README +++ b/examples/SPIN/gneb/README @@ -1,5 +1,10 @@ +Perform geodesic NEB calculations for spin configurations. +The two examples are: +- the magnetic switching of an iron nanoisland +- the collapse of a magnetic skyrmion + Run those examples as: -mpirun -np 4 lmp_mpi -in in.gneb.iron -partition 4x1 +mpirun -np 3 lmp_mpi -in in.gneb.iron -partition 3x1 You should be able to use any number of replicas >= 3. diff --git a/examples/SPIN/gneb/iron/in.gneb.iron b/examples/SPIN/gneb/iron/in.gneb.iron index 4cfbd723b7..c794292cfb 100644 --- a/examples/SPIN/gneb/iron/in.gneb.iron +++ b/examples/SPIN/gneb/iron/in.gneb.iron @@ -1,4 +1,3 @@ -# bcc iron in a 3d periodic box units metal dimension 3 @@ -11,6 +10,7 @@ atom_modify map array read_data initial.iron_spin # setting mass, mag. moments, and interactions for bcc iron +# (mass not necessary for fixed lattice calculation) mass 1 55.845 diff --git a/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion b/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion index aab6376e98..cbab56631b 100644 --- a/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion +++ b/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion @@ -1,4 +1,3 @@ -# bcc iron in a 3d periodic box units metal dimension 3 @@ -8,17 +7,12 @@ atom_style spin # necessary for the serial algorithm (sametag) atom_modify map array -#lattice sc 3.0 -#region box block 0.0 20.0 0.0 20.0 0.0 1.0 -#create_box 1 box -#create_atoms 1 box - read_data initial.skyrmion # setting mass, mag. moments, and interactions for bcc iron +# (mass not necessary for fixed lattice calculation) mass 1 55.845 -#set group all spin 2.2 -1.0 0.0 0.0 pair_style hybrid/overlay spin/exchange 3.1 spin/dmi 3.1 pair_coeff * * spin/exchange exchange 3.1 0.01593 0.06626915552 1.211 @@ -31,10 +25,8 @@ fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 anisotropy 5e-05 0.0 0.0 1.0 fix_modify 1 energy yes fix 2 all langevin/spin 0.0 0.0 21 fix 3 all neb/spin 1.0 -#fix 4 all nve/spin lattice no timestep 0.0001 -#run 0 compute out_mag all spin variable magx equal c_out_mag[1] @@ -49,10 +41,8 @@ thermo_modify format float %20.15g compute outsp all property/atom spx spy spz sp fmx fmy fmz variable u universe 1 2 3 4 -#dump 1 all custom 100 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7] dump 1 all custom 1 dump.$u type x y z c_outsp[1] c_outsp[2] c_outsp[3] min_style spin -min_modify alpha_damp 1.0 discret_factor 10.0 -neb/spin 1.0e-16 1.0e-16 1000 1000 10 final final.skyrmion -#neb/spin 1.0e-16 1.0e-16 10 10 10 final final.skyrmion +min_modify alpha_damp 1.0 discrete_factor 10.0 +neb/spin 1.0e-12 1.0e-12 10000 10000 10 final final.skyrmion diff --git a/src/SPIN/neb_spin.cpp b/src/SPIN/neb_spin.cpp index d981962a70..69c59e0484 100644 --- a/src/SPIN/neb_spin.cpp +++ b/src/SPIN/neb_spin.cpp @@ -637,52 +637,9 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) // initial, final and inter ang. values - double itheta,iphi,ftheta,fphi,ktheta,kphi; - double spix,spiy,spiz,spfx,spfy,spfz; - double spkx,spky,spkz,iknorm; - - spix = spi[0]; - spiy = spi[1]; - spiz = spi[2]; - - spfx = sploc[0]; - spfy = sploc[1]; - spfz = sploc[2]; - - iphi = itheta = fphi = ftheta = 0.0; - - iphi = acos(spiz); - if (sin(iphi) != 0.0) - itheta = acos(spix/sin(iphi)); - - fphi = acos(spfz); - if (sin(fphi) != 0.0) - ftheta = acos(spfx/sin(fphi)); - - kphi = iphi + fraction*(fphi-iphi); - ktheta = itheta + fraction*(ftheta-itheta); - - spkx = cos(ktheta)*sin(kphi); - spky = sin(ktheta)*sin(kphi); - spkz = cos(kphi); - - double knormsq = spkx*spkx + spky*spky + spkz*spkz; - if (knormsq != 0.0) - iknorm = 1.0/sqrt(knormsq); - - spkx *= iknorm; - spky *= iknorm; - spkz *= iknorm; - - //sploc[0] = spkx; - //sploc[1] = spky; - //sploc[2] = spkz; - - //double kx,ky,kz; + //double itheta,iphi,ftheta,fphi,ktheta,kphi; //double spix,spiy,spiz,spfx,spfy,spfz; - //double kcrossx,kcrossy,kcrossz,knormsq; - //double spkx,spky,spkz; - //double sdot,omega,iknorm; + //double spkx,spky,spkz,iknorm; //spix = spi[0]; //spiy = spi[1]; @@ -691,40 +648,89 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) //spfx = sploc[0]; //spfy = sploc[1]; //spfz = sploc[2]; - // - //kx = spiy*spfz - spiz*spfy; - //ky = spiz*spfx - spix*spfz; - //kz = spix*spfy - spiy*spfx; - //knormsq = kx*kx+ky*ky+kz*kz; - // - //if (knormsq != 0.0) { + //iphi = itheta = fphi = ftheta = 0.0; + + //iphi = acos(spiz); + //if (sin(iphi) != 0.0) + // itheta = acos(spix/sin(iphi)); + + //fphi = acos(spfz); + //if (sin(fphi) != 0.0) + // ftheta = acos(spfx/sin(fphi)); + + //kphi = iphi + fraction*(fphi-iphi); + //ktheta = itheta + fraction*(ftheta-itheta); + + //spkx = cos(ktheta)*sin(kphi); + //spky = sin(ktheta)*sin(kphi); + //spkz = cos(kphi); + + //double knormsq = spkx*spkx + spky*spky + spkz*spkz; + //if (knormsq != 0.0) // iknorm = 1.0/sqrt(knormsq); - // kx *= iknorm; - // ky *= iknorm; - // kz *= iknorm; - //} - // - //kcrossx = ky*spiz - kz*spiy; - //kcrossy = kz*spix - kx*spiz; - //kcrossz = kx*spiy - ky*spix; - - //sdot = spix*spfx + spiy*spfy + spiz*spfz; - - //omega = acos(sdot); - //omega *= fraction; - - //spkx = spix*cos(omega) + kcrossx*sin(omega); - //spky = spiy*cos(omega) + kcrossy*sin(omega); - //spkz = spiz*cos(omega) + kcrossz*sin(omega); - // - //iknorm = 1.0/sqrt(spkx*spkx+spky*spky+spkz*spkz); - //if (iknorm == 0.0) - // error->all(FLERR,"Incorrect rotation operation"); //spkx *= iknorm; //spky *= iknorm; //spkz *= iknorm; + + //sploc[0] = spkx; + //sploc[1] = spky; + //sploc[2] = spkz; + + double kx,ky,kz; + double spix,spiy,spiz,spfx,spfy,spfz; + double kcrossx,kcrossy,kcrossz,knormsq; + double kdots; + double spkx,spky,spkz; + double sdot,omega,iknorm,isnorm; + + spix = spi[0]; + spiy = spi[1]; + spiz = spi[2]; + + spfx = sploc[0]; + spfy = sploc[1]; + spfz = sploc[2]; + + kx = spiy*spfz - spiz*spfy; + ky = spiz*spfx - spix*spfz; + kz = spix*spfy - spiy*spfx; + + knormsq = kx*kx+ky*ky+kz*kz; + + if (knormsq != 0.0) { + iknorm = 1.0/sqrt(knormsq); + kx *= iknorm; + ky *= iknorm; + kz *= iknorm; + } + + kcrossx = ky*spiz - kz*spiy; + kcrossy = kz*spix - kx*spiz; + kcrossz = kx*spiy - ky*spix; + + kdots = kx*spix + ky*spiz + kz*spiz; + sdot = spix*spfx + spiy*spfy + spiz*spfz; + + omega = acos(sdot); + omega *= fraction; + + spkx = spix*cos(omega) + kcrossx*sin(omega); + spky = spiy*cos(omega) + kcrossy*sin(omega); + spkz = spiz*cos(omega) + kcrossz*sin(omega); + + spkx += kx*kdots*(1.0-cos(omega)); + spky += ky*kdots*(1.0-cos(omega)); + spkz += kz*kdots*(1.0-cos(omega)); + + isnorm = 1.0/sqrt(spkx*spkx+spky*spky+spkz*spkz); + if (isnorm == 0.0) + error->all(FLERR,"Incorrect rotation operation"); + + spkx *= isnorm; + spky *= isnorm; + spkz *= isnorm; sploc[0] = spkx; sploc[1] = spky;