Commit JT 040319

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.

# bcc iron in a 3d periodic box

units 		metal

dimension 	3

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

# bcc iron in a 3d periodic box

units 		metal

dimension 	3

# 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

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]

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

  // initial, final and inter ang. values 

  double itheta,iphi,ftheta,fphi,ktheta,kphi;

  //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 spix,spiy,spiz,spfx,spfy,spfz;

  double spkx,spky,spkz,iknorm;

  double kcrossx,kcrossy,kcrossz,knormsq;

  double kdots;

  double spkx,spky,spkz;

  double sdot,omega,iknorm,isnorm;

  spix = spi[0];

  spiy = spi[1];

  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);

  kx = spiy*spfz - spiz*spfy;

  ky = spiz*spfx - spix*spfz;

  kz = spix*spfy - spiy*spfx;

  spkx = cos(ktheta)*sin(kphi);

  spky = sin(ktheta)*sin(kphi);

  spkz = cos(kphi);

  knormsq = kx*kx+ky*ky+kz*kz;

  double knormsq = spkx*spkx + spky*spky + spkz*spkz;

  if (knormsq != 0.0)

  if (knormsq != 0.0) {

    iknorm = 1.0/sqrt(knormsq);

    kx *= iknorm;

    ky *= iknorm;

    kz *= iknorm;

  }

  spkx *= iknorm;

  spky *= iknorm;

  spkz *= iknorm;

  kcrossx = ky*spiz - kz*spiy;

  kcrossy = kz*spix - kx*spiz;

  kcrossz = kx*spiy - ky*spix;

  //sploc[0] = spkx;

  //sploc[1] = spky;

  //sploc[2] = spkz;

  kdots = kx*spix + ky*spiz + kz*spiz;

  sdot = spix*spfx + spiy*spfy + spiz*spfz;

  //double kx,ky,kz;

  //double spix,spiy,spiz,spfx,spfy,spfz;

  //double kcrossx,kcrossy,kcrossz,knormsq;

  //double spkx,spky,spkz;

  //double sdot,omega,iknorm;

  omega = acos(sdot);

  omega *= fraction;

  //spix = spi[0];

  //spiy = spi[1];

  //spiz = spi[2];

  spkx = spix*cos(omega) + kcrossx*sin(omega); 

  spky = spiy*cos(omega) + kcrossy*sin(omega); 

  spkz = spiz*cos(omega) + kcrossz*sin(omega); 

  //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;

  //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 += kx*kdots*(1.0-cos(omega));

  spky += ky*kdots*(1.0-cos(omega));

  spkz += kz*kdots*(1.0-cos(omega));

  //spkx *= iknorm;

  //spky *= iknorm;

  //spkz *= iknorm;

  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;