diff --git a/doc/src/Eqs/neb_spin_k.jpg b/doc/src/Eqs/neb_spin_k.jpg new file mode 100644 index 0000000000..add309694f Binary files /dev/null and b/doc/src/Eqs/neb_spin_k.jpg differ diff --git a/doc/src/Eqs/neb_spin_k.tex b/doc/src/Eqs/neb_spin_k.tex new file mode 100644 index 0000000000..f0ce8e180e --- /dev/null +++ b/doc/src/Eqs/neb_spin_k.tex @@ -0,0 +1,16 @@ +\documentclass[preview]{standalone} +\usepackage{varwidth} +\usepackage[utf8x]{inputenc} +\usepackage{amsmath, amssymb, graphics, setspace} + +\begin{document} +\begin{varwidth}{50in} + \begin{equation} + \vec{k}_i = + \frac{\vec{m}_i^I \times \vec{m}_i^F}{\left|\vec{m}_i^I + \times \vec{m}_i^F\right|} + %&{\rm ~if~}& \vec{m}_i^I \times \vec{m}_i^F + , \nonumber + \end{equation} +\end{varwidth} +\end{document} diff --git a/doc/src/Eqs/neb_spin_rodrigues_formula.jpg b/doc/src/Eqs/neb_spin_rodrigues_formula.jpg new file mode 100644 index 0000000000..66070f7bc5 Binary files /dev/null and b/doc/src/Eqs/neb_spin_rodrigues_formula.jpg differ diff --git a/doc/src/Eqs/neb_spin_rodrigues_formula.tex b/doc/src/Eqs/neb_spin_rodrigues_formula.tex new file mode 100644 index 0000000000..4a8347cd79 --- /dev/null +++ b/doc/src/Eqs/neb_spin_rodrigues_formula.tex @@ -0,0 +1,16 @@ +\documentclass[preview]{standalone} +\usepackage{varwidth} +\usepackage[utf8x]{inputenc} +\usepackage{amsmath, amssymb, graphics, setspace} + +\begin{document} +\begin{varwidth}{50in} + \begin{equation} + \vec{m}_i^{\nu} = \vec{m}_i^{I} \cos(\omega_i^{\nu}) + + (\vec{k}_i \times \vec{m}_i^{I}) \sin(\omega_i^{\nu}) + + (1.0-\cos(\omega_i^{\nu})) \vec{k}_i (\vec{k}_i\cdot + \vec{m}_i^{I}) + , \nonumber + \end{equation} +\end{varwidth} +\end{document} diff --git a/doc/src/neb_spin.txt b/doc/src/neb_spin.txt new file mode 100644 index 0000000000..33cb4cc2ed --- /dev/null +++ b/doc/src/neb_spin.txt @@ -0,0 +1,430 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +neb command :h3 + +[Syntax:] + +neb/spin etol ttol N1 N2 Nevery file-style arg keyword :pre + +etol = stopping tolerance for energy (energy units) :ulb,l +ttol = stopping tolerance for torque ( units) :l +N1 = max # of iterations (timesteps) to run initial NEB :l +N2 = max # of iterations (timesteps) to run barrier-climbing NEB :l +Nevery = print replica energies and reaction coordinates every this many timesteps :l +file-style = {final} or {each} or {none} :l + {final} arg = filename + filename = file with initial coords for final replica + coords for intermediate replicas are linearly interpolated + between first and last replica + {each} arg = filename + filename = unique filename for each replica (except first) + with its initial coords + {none} arg = no argument all replicas assumed to already have + their initial coords :pre +keyword = {verbose} +:ule + +[Examples:] + +neb/spin 0.1 0.0 1000 500 50 final coords.final +neb/spin 0.0 0.001 1000 500 50 each coords.initial.$i +neb/spin 0.0 0.001 1000 500 50 none verbose :pre + +[Description:] + +Perform a geodesic nudged elastic band (GNEB) calculation using multiple +replicas of a system. Two or more replicas must be used; the first +and last are the end points of the transition path. + +GNEB is a method for finding both the spin configurations and height +of the energy barrier associated with a transition state, e.g. +spins to perform a collective rotation from one energy basin to +another. +The implementation in LAMMPS follows the discussion in the +following paper: "(Bessarab)"_#Bessarab. + +Each replica runs on a partition of one or more processors. Processor +partitions are defined at run-time using the "-partition command-line +switch"_Run_options.html. Note that if you have MPI installed, you +can run a multi-replica simulation with more replicas (partitions) +than you have physical processors, e.g you can run a 10-replica +simulation on just one or two processors. You will simply not get the +performance speed-up you would see with one or more physical +processors per replica. See the "Howto replica"_Howto_replica.html +doc page for further discussion. + +NOTE: As explained below, a GNEB calculation performs a damped dynamics +minimization across all the replicas. The "spin"_min_spin.html +style minimizer has to be defined in your input script. + +When a GNEB calculation is performed, it is assumed that each replica +is running the same system, though LAMMPS does not check for this. +I.e. the simulation domain, the number of magnetic atoms, the +interaction potentials, and the starting configuration when the neb +command is issued should be the same for every replica. + +In a GNEB calculation each replica is connected to other replicas by +inter-replica nudging forces. These forces are imposed by the "fix +neb/spin"_fix_neb_spin.html command, which must be used in conjunction +with the neb command. +The group used to define the fix neb/spin command defines the +GNEB magnetic atoms which are the only ones that inter-replica springs +are applied to. +If the group does not include all magnetic atoms, then non-GNEB +magnetic atoms have no inter-replica springs and the torques they feel +and their precessional motion is computed in the usual way due only +to other magnetic atoms within their replica. +Conceptually, the non-GNEB atoms provide a background force field for +the GNEB atoms. +Their magnetic spins can be allowed to precess during the GNEB +minimization procedure. + +The initial spin configuration for each of the replicas can be +specified in different manners via the {file-style} setting, as +discussed below. Only atomic spins whose initial coordinates should +differ from the current configuration need to be specified. + +Conceptually, the initial and final configurations for the first +replica should be states on either side of an energy barrier. + +As explained below, the initial configurations of intermediate +replicas can be spin coordinates interpolated in a linear fashion +between the first and last replicas. This is often adequate for +simple transitions. For more complex transitions, it may lead to slow +convergence or even bad results if the minimum energy path (MEP, see +below) of states over the barrier cannot be correctly converged to +from such an initial path. In this case, you will want to generate +initial states for the intermediate replicas that are geometrically +closer to the MEP and read them in. + +################################################################### + +:line + +For a {file-style} setting of {final}, a filename is specified which +contains atomic and spin coordinates for zero or more atoms, in the +format described below. +For each atom that appears in the file, the new coordinates are +assigned to that atom in the final replica. Each intermediate replica +also assigns a new spin to that atom in an interpolated manner. +This is done by using the current direction of the spin at the starting +point and the read-in direction as the final point. +The angular distance between them is calculated, and the new direction +is assigned to be a fraction of the angular distance. + +NOTE: The "angular distance" between the starting and final point is +evaluated in the geodesic sense, as described in "(Bessarab)"_#Bessarab. + +NOTE: The angular interpolation between the starting and final point +is achieved using Rodrigues formula: + +:c,image(Eqs/neb_spin_rodrigues_formula.jpg) + +with m_i^I is the initial spin configuration for the spin i, +where the rotation and k_i is defined as: + +:c,image(Eqs/neb_spin_k.jpg) + +The distance between them is calculated, and the new position +is assigned to be a fraction of the distance. E.g. if there are 10 +replicas, the 2nd replica will assign a position that is 10% of the +distance along a line between the starting and final point, and the +9th replica will assign a position that is 90% of the distance along +the line. Note that for this procedure to produce consistent +coordinates across all the replicas, the current coordinates need to +be the same in all replicas. LAMMPS does not check for this, but +invalid initial configurations will likely result if it is not the +case. + +NOTE: The "distance" between the starting and final point is +calculated in a minimum-image sense for a periodic simulation box. +This means that if the two positions are on opposite sides of a box +(periodic in that dimension), the distance between them will be small, +because the periodic image of one of the atoms is close to the other. +Similarly, even if the assigned position resulting from the +interpolation is outside the periodic box, the atom will be wrapped +back into the box when the NEB calculation begins. + +For a {file-style} setting of {each}, a filename is specified which is +assumed to be unique to each replica. This can be done by using a +variable in the filename, e.g. + +variable i equal part +neb 0.0 0.001 1000 500 50 each coords.initial.$i :pre + +which in this case will substitute the partition ID (0 to N-1) for the +variable I, which is also effectively the replica ID. See the +"variable"_variable.html command for other options, such as using +world-, universe-, or uloop-style variables. + +Each replica (except the first replica) will read its file, formatted +as described below, and for any atom that appears in the file, assign +the specified coordinates to its atom. The various files do not need +to contain the same set of atoms. + +For a {file-style} setting of {none}, no filename is specified. Each +replica is assumed to already be in its initial configuration at the +time the neb command is issued. This allows each replica to define +its own configuration by reading a replica-specific data or restart or +dump file, via the "read_data"_read_data.html, +"read_restart"_read_restart.html, or "read_dump"_read_dump.html +commands. The replica-specific names of these files can be specified +as in the discussion above for the {each} file-style. Also see the +section below for how a NEB calculation can produce restart files, so +that a long calculation can be restarted if needed. + +NOTE: None of the {file-style} settings change the initial +configuration of any atom in the first replica. The first replica +must thus be in the correct initial configuration at the time the neb +command is issued. + +:line + +A NEB calculation proceeds in two stages, each of which is a +minimization procedure, performed via damped dynamics. To enable +this, you must first define a damped dynamics +"min_style"_min_style.html, such as {quickmin} or {fire}. The {cg}, +{sd}, and {hftn} styles cannot be used, since they perform iterative +line searches in their inner loop, which cannot be easily synchronized +across multiple replicas. + +The minimizer tolerances for energy and force are set by {etol} and +{ftol}, the same as for the "minimize"_minimize.html command. + +A non-zero {etol} means that the NEB calculation will terminate if the +energy criterion is met by every replica. The energies being compared +to {etol} do not include any contribution from the inter-replica +nudging forces, since these are non-conservative. A non-zero {ftol} +means that the NEB calculation will terminate if the force criterion +is met by every replica. The forces being compared to {ftol} include +the inter-replica nudging forces. + +The maximum number of iterations in each stage is set by {N1} and +{N2}. These are effectively timestep counts since each iteration of +damped dynamics is like a single timestep in a dynamics +"run"_run.html. During both stages, the potential energy of each +replica and its normalized distance along the reaction path (reaction +coordinate RD) will be printed to the screen and log file every +{Nevery} timesteps. The RD is 0 and 1 for the first and last replica. +For intermediate replicas, it is the cumulative distance (normalized +by the total cumulative distance) between adjacent replicas, where +"distance" is defined as the length of the 3N-vector of differences in +atomic coordinates, where N is the number of NEB atoms involved in the +transition. These outputs allow you to monitor NEB's progress in +finding a good energy barrier. {N1} and {N2} must both be multiples +of {Nevery}. + +In the first stage of NEB, the set of replicas should converge toward +a minimum energy path (MEP) of conformational states that transition +over a barrier. The MEP for a transition is defined as a sequence of +3N-dimensional states, each of which has a potential energy gradient +parallel to the MEP itself. The configuration of highest energy along +a MEP corresponds to a saddle point. The replica states will also be +roughly equally spaced along the MEP due to the inter-replica nudging +force added by the "fix neb"_fix_neb.html command. + +In the second stage of NEB, the replica with the highest energy is +selected and the inter-replica forces on it are converted to a force +that drives its atom coordinates to the top or saddle point of the +barrier, via the barrier-climbing calculation described in +"(HenkelmanB)"_#HenkelmanB. As before, the other replicas rearrange +themselves along the MEP so as to be roughly equally spaced. + +When both stages are complete, if the NEB calculation was successful, +the configurations of the replicas should be along (close to) the MEP +and the replica with the highest energy should be an atomic +configuration at (close to) the saddle point of the transition. The +potential energies for the set of replicas represents the energy +profile of the transition along the MEP. + +:line + +A few other settings in your input script are required or advised to +perform a NEB calculation. See the NOTE about the choice of timestep +at the beginning of this doc page. + +An atom map must be defined which it is not by default for "atom_style +atomic"_atom_style.html problems. The "atom_modify +map"_atom_modify.html command can be used to do this. + +The minimizers in LAMMPS operate on all atoms in your system, even +non-NEB atoms, as defined above. To prevent non-NEB atoms from moving +during the minimization, you should use the "fix +setforce"_fix_setforce.html command to set the force on each of those +atoms to 0.0. This is not required, and may not even be desired in +some cases, but if those atoms move too far (e.g. because the initial +state of your system was not well-minimized), it can cause problems +for the NEB procedure. + +The damped dynamics "minimizers"_min_style.html, such as {quickmin} +and {fire}), adjust the position and velocity of the atoms via an +Euler integration step. Thus you must define an appropriate +"timestep"_timestep.html to use with NEB. As mentioned above, NEB +will often converge more quickly if you use a timestep about 10x +larger than you would normally use for dynamics simulations. + +:line + +Each file read by the neb command containing atomic coordinates used +to initialize one or more replicas must be formatted as follows. + +The file can be ASCII text or a gzipped text file (detected by a .gz +suffix). The file can contain initial blank lines or comment lines +starting with "#" which are ignored. The first non-blank, non-comment +line should list N = the number of lines to follow. The N successive +lines contain the following information: + +ID1 x1 y1 z1 +ID2 x2 y2 z2 +... +IDN xN yN zN :pre + +The fields are the atom ID, followed by the x,y,z coordinates. The +lines can be listed in any order. Additional trailing information on +the line is OK, such as a comment. + +Note that for a typical NEB calculation you do not need to specify +initial coordinates for very many atoms to produce differing starting +and final replicas whose intermediate replicas will converge to the +energy barrier. Typically only new coordinates for atoms +geometrically near the barrier need be specified. + +Also note there is no requirement that the atoms in the file +correspond to the NEB atoms in the group defined by the "fix +neb"_fix_neb.html command. Not every NEB atom need be in the file, +and non-NEB atoms can be listed in the file. + +:line + +Four kinds of output can be generated during a NEB calculation: energy +barrier statistics, thermodynamic output by each replica, dump files, +and restart files. + +When running with multiple partitions (each of which is a replica in +this case), the print-out to the screen and master log.lammps file +contains a line of output, printed once every {Nevery} timesteps. It +contains the timestep, the maximum force per replica, the maximum +force per atom (in any replica), potential gradients in the initial, +final, and climbing replicas, the forward and backward energy +barriers, the total reaction coordinate (RDT), and the normalized +reaction coordinate and potential energy of each replica. + +The "maximum force per replica" is the two-norm of the 3N-length force +vector for the atoms in each replica, maximized across replicas, which +is what the {ftol} setting is checking against. In this case, N is +all the atoms in each replica. The "maximum force per atom" is the +maximum force component of any atom in any replica. The potential +gradients are the two-norm of the 3N-length force vector solely due to +the interaction potential i.e. without adding in inter-replica +forces. + +The "reaction coordinate" (RD) for each replica is the two-norm of the +3N-length vector of distances between its atoms and the preceding +replica's atoms, added to the RD of the preceding replica. The RD of +the first replica RD1 = 0.0; the RD of the final replica RDN = RDT, +the total reaction coordinate. The normalized RDs are divided by RDT, +so that they form a monotonically increasing sequence from zero to +one. When computing RD, N only includes the atoms being operated on by +the fix neb command. + +The forward (reverse) energy barrier is the potential energy of the +highest replica minus the energy of the first (last) replica. + +Supplementary information for all replicas can be printed out to the +screen and master log.lammps file by adding the verbose keyword. This +information include the following. The "path angle" (pathangle) for +the replica i which is the angle between the 3N-length vectors (Ri-1 - +Ri) and (Ri+1 - Ri) (where Ri is the atomic coordinates of replica +i). A "path angle" of 180 indicates that replicas i-1, i and i+1 are +aligned. "angletangrad" is the angle between the 3N-length tangent +vector and the 3N-length force vector at image i. The tangent vector +is calculated as in "(HenkelmanA)"_#HenkelmanA for all intermediate +replicas and at R2 - R1 and RM - RM-1 for the first and last replica, +respectively. "anglegrad" is the angle between the 3N-length energy +gradient vector of replica i and that of replica i+1. It is not +defined for the final replica and reads nan. gradV is the norm of the +energy gradient of image i. ReplicaForce is the two-norm of the +3N-length force vector (including nudging forces) for replica i. +MaxAtomForce is the maximum force component of any atom in replica i. + +When a NEB calculation does not converge properly, the supplementary +information can help understanding what is going wrong. For instance +when the path angle becomes acute, the definition of tangent used in +the NEB calculation is questionable and the NEB cannot may diverge +"(Maras)"_#Maras2. + + +When running on multiple partitions, LAMMPS produces additional log +files for each partition, e.g. log.lammps.0, log.lammps.1, etc. For a +NEB calculation, these contain the thermodynamic output for each +replica. + +If "dump"_dump.html commands in the input script define a filename +that includes a {universe} or {uloop} style "variable"_variable.html, +then one dump file (per dump command) will be created for each +replica. At the end of the NEB calculation, the final snapshot in +each file will contain the sequence of snapshots that transition the +system over the energy barrier. Earlier snapshots will show the +convergence of the replicas to the MEP. + +Likewise, "restart"_restart.html filenames can be specified with a +{universe} or {uloop} style "variable"_variable.html, to generate +restart files for each replica. These may be useful if the NEB +calculation fails to converge properly to the MEP, and you wish to +restart the calculation from an intermediate point with altered +parameters. + +There are 2 Python scripts provided in the tools/python directory, +neb_combine.py and neb_final.py, which are useful in analyzing output +from a NEB calculation. Assume a NEB simulation with M replicas, and +the NEB atoms labeled with a specific atom type. + +The neb_combine.py script extracts atom coords for the NEB atoms from +all M dump files and creates a single dump file where each snapshot +contains the NEB atoms from all the replicas and one copy of non-NEB +atoms from the first replica (presumed to be identical in other +replicas). This can be visualized/animated to see how the NEB atoms +relax as the NEB calculation proceeds. + +The neb_final.py script extracts the final snapshot from each of the M +dump files to create a single dump file with M snapshots. This can be +visualized to watch the system make its transition over the energy +barrier. + +To illustrate, here are images from the final snapshot produced by the +neb_combine.py script run on the dump files produced by the two +example input scripts in examples/neb. Click on them to see a larger +image. + +:image(JPG/hop1_small.jpg,JPG/hop1.jpg) +:image(JPG/hop2_small.jpg,JPG/hop2.jpg) + +:line + +[Restrictions:] + +This command can only be used if LAMMPS was built with the SPIN +package. See the "Build package"_Build_package.html doc +page for more info. + +:line + +[Related commands:] + +"min/spin"_min_spin.html, "fix neb/spin"_fix_neb_spin.html + +[Default:] + +none + +:line + +:link(Bessarab) +[(Bessarab)] Bessarab, Uzdin, Jonsson, Comp Phys Comm, 196, +335-347 (2015). diff --git a/examples/SPIN/gneb/iron/in.gneb.iron b/examples/SPIN/gneb/iron/in.gneb.iron index c794292cfb..95e7071cb0 100644 --- a/examples/SPIN/gneb/iron/in.gneb.iron +++ b/examples/SPIN/gneb/iron/in.gneb.iron @@ -7,11 +7,10 @@ atom_style spin # necessary for the serial algorithm (sametag) 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) +read_data initial.iron_spin mass 1 55.845 pair_style spin/exchange 3.5 diff --git a/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion b/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion index cbab56631b..cf55c9d1d4 100644 --- a/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion +++ b/examples/SPIN/gneb/skyrmion/in.gneb.skyrmion @@ -7,11 +7,10 @@ atom_style spin # necessary for the serial algorithm (sametag) atom_modify map array -read_data initial.skyrmion - # setting mass, mag. moments, and interactions for bcc iron # (mass not necessary for fixed lattice calculation) +read_data initial.skyrmion mass 1 55.845 pair_style hybrid/overlay spin/exchange 3.1 spin/dmi 3.1 diff --git a/src/SPIN/fix_neb_spin.h b/src/SPIN/fix_neb_spin.h index 9bbacc8bf0..8e016b2e23 100644 --- a/src/SPIN/fix_neb_spin.h +++ b/src/SPIN/fix_neb_spin.h @@ -60,20 +60,20 @@ class FixNEB_spin : public Fix { double **spprev,**spnext,**fmnext; double **springF; double **tangent; - double **xsend,**xrecv; // coords to send/recv to/from other replica - double **fsend,**frecv; // coords to send/recv to/from other replica - double **spsend,**sprecv; // sp to send/recv to/from other replica - double **fmsend,**fmrecv; // fm to send/recv to/from other replica - tagint *tagsend,*tagrecv; // ditto for atom IDs + double **xsend,**xrecv; // coords to send/recv to/from other replica + double **fsend,**frecv; // coords to send/recv to/from other replica + double **spsend,**sprecv; // sp to send/recv to/from other replica + double **fmsend,**fmrecv; // fm to send/recv to/from other replica + tagint *tagsend,*tagrecv; // ditto for atom IDs - // info gathered from all procs in my replica - double **xsendall,**xrecvall; // coords to send/recv to/from other replica - double **fsendall,**frecvall; // force to send/recv to/from other replica - double **spsendall,**sprecvall; // sp to send/recv to/from other replica - double **fmsendall,**fmrecvall; // fm to send/recv to/from other replica - tagint *tagsendall,*tagrecvall; // ditto for atom IDs + // info gathered from all procs in my replica + double **xsendall,**xrecvall; // coords to send/recv to/from other replica + double **fsendall,**frecvall; // force to send/recv to/from other replica + double **spsendall,**sprecvall; // sp to send/recv to/from other replica + double **fmsendall,**fmrecvall; // fm to send/recv to/from other replica + tagint *tagsendall,*tagrecvall; // ditto for atom IDs - int *counts,*displacements; // used for MPI_Gather + int *counts,*displacements; // used for MPI_Gather double geodesic_distance(double *, double *); void inter_replica_comm(); @@ -97,16 +97,16 @@ E: Potential energy ID for fix neb does not exist Self-explanatory. -E: Too many active NEB atoms +E: Too many active GNEB atoms UNDOCUMENTED -E: Too many atoms for NEB +E: Too many atoms for GNEB UNDOCUMENTED U: Atom count changed in fix neb -This is not allowed in a NEB calculation. +This is not allowed in a GNEB calculation. */ diff --git a/src/SPIN/neb_spin.cpp b/src/SPIN/neb_spin.cpp index 69c59e0484..77a94c5e84 100644 --- a/src/SPIN/neb_spin.cpp +++ b/src/SPIN/neb_spin.cpp @@ -111,6 +111,8 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, double **sp = atom->sp; int nlocal = atom->nlocal; + int temp_flag,rot_flag; + temp_flag = rot_flag = 0; int ii = 0; double spinit[3],spfinal[3]; for (int i = 0; i < nlocal; i++) { @@ -123,7 +125,7 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, spfinal[2] = buf_final[ii+2]; // interpolate intermediate spin states - + if (fraction == 0.0) { sp[i][0] = spinit[0]; sp[i][1] = spinit[1]; @@ -133,7 +135,8 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, sp[i][1] = spfinal[1]; sp[i][2] = spfinal[2]; } else { - initial_rotation(spinit,spfinal,fraction); + temp_flag = initial_rotation(spinit,spfinal,fraction); + rot_flag = MAX(temp_flag,rot_flag); sp[i][0] = spfinal[0]; sp[i][1] = spfinal[1]; sp[i][2] = spfinal[2]; @@ -141,6 +144,14 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, ii += 3; } + + // warning message if one or more couples (spi,spf) were aligned + // this breaks Rodrigues' formula, and an arbitrary rotation + // vector has to be chosen + + if ((rot_flag > 0) && (comm->me == 0)) + error->warning(FLERR,"arbitrary initial rotation of one or more spin(s)"); + } /* ---------------------------------------------------------------------- */ @@ -494,6 +505,8 @@ void NEB_spin::readfile(char *file, int flag) int ncount = 0; + int temp_flag,rot_flag; + temp_flag = rot_flag = 0; int nread = 0; while (nread < nlines) { nchunk = MIN(nlines-nread,CHUNK); @@ -566,7 +579,8 @@ void NEB_spin::readfile(char *file, int flag) sp[m][1] = spfinal[1]; sp[m][2] = spfinal[2]; } else { - initial_rotation(spinit,spfinal,fraction); + temp_flag = initial_rotation(spinit,spfinal,fraction); + rot_flag = MAX(temp_flag,rot_flag); sp[m][0] = spfinal[0]; sp[m][1] = spfinal[1]; sp[m][2] = spfinal[2]; @@ -588,6 +602,13 @@ void NEB_spin::readfile(char *file, int flag) nread += nchunk; } + // warning message if one or more couples (spi,spf) were aligned + // this breaks Rodrigues' formula, and an arbitrary rotation + // vector has to be chosen + + if ((rot_flag > 0) && (comm->me == 0)) + error->warning(FLERR,"arbitrary initial rotation of one or more spin(s)"); + // check that all atom IDs in file were found by a proc if (flag == 0) { @@ -621,69 +642,24 @@ void NEB_spin::readfile(char *file, int flag) } /* ---------------------------------------------------------------------- - initial configuration of spin sploc using Rodrigues' formula + initial configuration of intermediate spins using Rodrigues' formula interpolates between initial (spi) and final (stored in sploc) ------------------------------------------------------------------------- */ -void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) +int NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) { - // implementing initial rotation using atan2 - // this may not be a sufficient routine, need more accurate verifications + // no interpolation for initial and final replica - // interpolation only for intermediate replica + if (fraction == 0.0 || fraction == 1.0) return 0; - if (fraction == 0.0 || fraction == 1.0) return; - - // 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; - + int rot_flag = 0; 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; + double sidotsf,omega,iknorm,isnorm; spix = spi[0]; spiy = spi[1]; @@ -698,43 +674,73 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) 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; + sidotsf = spix*spfx + spiy*spfy + spiz*spfz; + + // if knormsq == 0.0, init and final spins are aligned + // Rodrigues' formula breaks, needs to define another axis k + + if (knormsq == 0.0) { + if (sidotsf > 0.0) { // spins aligned and in same direction + return 0; + } else if (sidotsf < 0.0) { // spins aligned and in opposite directions + + // defining a rotation axis + // first guess, k = spi x [100] + // second guess, k = spi x [010] + + if (spiy*spiy + spiz*spiz != 0.0) { // spin not along [100] + kx = 0.0; + ky = spiz; + kz = -spiy; + } else if (spix*spix + spiz*spiz != 0.0) { // spin not along [010] + kx = -spiz; + ky = 0.0; + kz = spix; + } else error->all(FLERR,"Incorrect initial rotation operation"); + rot_flag = 1; + } } - + 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 = acos(sidotsf); omega *= fraction; - spkx = spix*cos(omega) + kcrossx*sin(omega); - spky = spiy*cos(omega) + kcrossy*sin(omega); - spkz = spiz*cos(omega) + kcrossz*sin(omega); + // applying Rodrigues' formula + + spkx = spix*cos(omega); + spky = spiy*cos(omega); + spkz = spiz*cos(omega); + + spkx += kcrossx*sin(omega); + spky += kcrossy*sin(omega); + spkz += kcrossz*sin(omega); spkx += kx*kdots*(1.0-cos(omega)); spky += ky*kdots*(1.0-cos(omega)); spkz += kz*kdots*(1.0-cos(omega)); + // normalizing resulting spin vector + isnorm = 1.0/sqrt(spkx*spkx+spky*spky+spkz*spkz); if (isnorm == 0.0) - error->all(FLERR,"Incorrect rotation operation"); + error->all(FLERR,"Incorrect initial rotation operation"); spkx *= isnorm; spky *= isnorm; spkz *= isnorm; + // returns rotated spin + sploc[0] = spkx; sploc[1] = spky; sploc[2] = spkz; + + return rot_flag; } /* ---------------------------------------------------------------------- diff --git a/src/SPIN/neb_spin.h b/src/SPIN/neb_spin.h index 5988c04a3a..b7c20bc3a9 100644 --- a/src/SPIN/neb_spin.h +++ b/src/SPIN/neb_spin.h @@ -57,7 +57,7 @@ class NEB_spin : protected Pointers { double *fmaxatomInRepl; // force on an image void readfile(char *, int); - void initial_rotation(double *, double *, double); + int initial_rotation(double *, double *, double); void open(char *); void print_status(); }; diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp index 5e9ff7a39e..41430d230f 100644 --- a/src/SPIN/pair_spin_dmi.cpp +++ b/src/SPIN/pair_spin_dmi.cpp @@ -171,10 +171,11 @@ void PairSpinDmi::init_style() int ifix = 0; while (ifix < modify->nfix) { if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break; + if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break; ifix++; } if ((ifix == modify->nfix) && (comm->me == 0)) - error->warning(FLERR,"Using pair/spin style without nve/spin"); + error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin"); // get the lattice_flag from nve/spin diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index 71b4c2ebf6..0260a611cf 100644 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -162,7 +162,7 @@ void PairSpinExchange::init_style() ifix++; } if ((ifix == modify->nfix) && (comm->me == 0)) - error->warning(FLERR,"Using pair/spin style without nve/spin"); + error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin"); // get the lattice_flag from nve/spin diff --git a/src/SPIN/pair_spin_magelec.cpp b/src/SPIN/pair_spin_magelec.cpp index 6ff003521d..1f1488b93c 100644 --- a/src/SPIN/pair_spin_magelec.cpp +++ b/src/SPIN/pair_spin_magelec.cpp @@ -164,10 +164,11 @@ void PairSpinMagelec::init_style() int ifix = 0; while (ifix < modify->nfix) { if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break; + if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break; ifix++; } if ((ifix == modify->nfix) && (comm->me == 0)) - error->warning(FLERR,"Using pair/spin style without nve/spin"); + error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin"); // get the lattice_flag from nve/spin diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp index a39d6f3461..03041da17f 100644 --- a/src/SPIN/pair_spin_neel.cpp +++ b/src/SPIN/pair_spin_neel.cpp @@ -171,10 +171,11 @@ void PairSpinNeel::init_style() int ifix = 0; while (ifix < modify->nfix) { if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break; + if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break; ifix++; } if ((ifix == modify->nfix) && (comm->me == 0)) - error->warning(FLERR,"Using pair/spin style without nve/spin"); + error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin"); // get the lattice_flag from nve/spin