lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1380 from julient31/minimizer_spin 

Add a spin minimizer in the SPIN package
This commit is contained in:
Axel Kohlmeyer 2019-03-26 11:52:21 -04:00 committed by GitHub
parent 89bc79f657 19d25203ad
commit e2e4fe2cf7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
27 changed files with 769 additions and 131 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_all.txt
View File

@ -79,6 +79,7 @@ An alphabetic list of all general LAMMPS commands.
"minimize"_minimize.html,
"min_modify"_min_modify.html,
"min_style"_min_style.html,
"min_style spin"_min_spin.html,
"molecule"_molecule.html,
"ndx2group"_group2ndx.html,
"neb"_neb.html,

BIN
doc/src/Eqs/min_spin_damping.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 6.9 KiB

13
doc/src/Eqs/min_spin_damping.tex Normal file
View File

@ -0,0 +1,13 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath, amssymb, graphics, setspace}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
\frac{d \vec{s}_{i}}{dt} = \lambda\, \vec{s}_{i} \times\left( \vec{\omega}_{i} \times\vec{s}_{i} \right)
\nonumber
\end{equation}
\end{varwidth}
\end{document}

BIN
doc/src/Eqs/min_spin_timestep.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 5.8 KiB

14
doc/src/Eqs/min_spin_timestep.tex Normal file
View File

@ -0,0 +1,14 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath, amssymb, graphics, setspace}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
{\Delta t}_{\rm max} = \frac{2\pi}{\kappa
\left|\vec{\omega}_{\rm max} \right|}
\nonumber
\end{equation}
\end{varwidth}
\end{document}

1
doc/src/commands_list.txt
View File

@ -62,6 +62,7 @@ Commands :h1
mass
message
min_modify
min_spin
min_style
minimize
molecule

1
doc/src/lammps.book
View File

@ -175,6 +175,7 @@ mass.html
message.html
min_modify.html
min_style.html
min_spin.html
minimize.html
molecule.html
neb.html

19
doc/src/min_modify.txt
View File

@ -13,11 +13,15 @@ min_modify command :h3
min_modify keyword values ... :pre
one or more keyword/value pairs may be listed :ulb,l
keyword = {dmax} or {line}
keyword = {dmax} or {line} or {alpha_damp} or {discrete_factor}
{dmax} value = max
max = maximum distance for line search to move (distance units)
{line} value = {backtrack} or {quadratic} or {forcezero}
backtrack,quadratic,forcezero = style of linesearch to use :pre
backtrack,quadratic,forcezero = style of linesearch to use
{alpha_damp} value = damping
damping = fictitious Gilbert damping for spin minimization (adim)
{discrete_factor} value = factor
factor = discretization factor for adaptive spin timestep (adim) :pre
:ule
[Examples:]
@ -65,6 +69,17 @@ difference of two large values (energy before and energy after) and
that difference may be smaller than machine epsilon even if atoms
could move in the gradient direction to reduce forces further.
Keywords {alpha_damp} and {discrete_factor} only make sense when
a "min_spin"_min_spin.html command is declared.
Keyword {alpha_damp} defines an analog of a magnetic Gilbert
damping. It defines a relaxation rate toward an equilibrium for
a given magnetic system.
Keyword {discrete_factor} defines a discretization factor for the
adaptive timestep used in the {spin} minimization.
See "min_spin"_min_spin.html for more information about those
quantities.
Default values are {alpha_damp} = 1.0 and {discrete_factor} = 10.0.
[Restrictions:] none
[Related commands:]

65
doc/src/min_spin.txt Normal file
View File

@ -0,0 +1,65 @@
"LAMMPS WWW Page"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Commands_all.html)
:line
min_style spin command :h3
[Syntax:]
min_style spin :pre
[Examples:]
min_style spin :pre
[Description:]
Apply a minimization algorithm to use when a "minimize"_minimize.html
command is performed.
Style {spin} defines a damped spin dynamics with an adaptive
timestep, according to:
:c,image(Eqs/min_spin_damping.jpg)
with lambda a damping coefficient (similar to a Gilbert
damping).
Lambda can be defined by setting the {alpha_damp} keyword with the
"min_modify"_min_modify.html command.
The minimization procedure solves this equation using an
adaptive timestep. The value of this timestep is defined
by the largest precession frequency that has to be solved in the
system:
:c,image(Eqs/min_spin_timestep.jpg)
with {|omega|_{max}} the norm of the largest precession frequency
in the system (across all processes, and across all replicas if a
spin/neb calculation is performed).
Kappa defines a discretization factor {discrete_factor} for the
definition of this timestep.
{discrete_factor} can be defined with the "min_modify"_min_modify.html
command.
NOTE: The {spin} style replaces the force tolerance by a torque
tolerance. See "minimize"_minimize.html for more explanation.
[Restrictions:]
This minimization procedure is only applied to spin degrees of
freedom for a frozen lattice configuration.
[Related commands:]
"min_style"_min_style.html, "minimize"_minimize.html,
"min_modify"_min_modify.html
[Default:]
The option defaults are {alpha_damp} = 1.0 and {discrete_factor} =
10.0.

7
doc/src/min_style.txt
View File

@ -11,11 +11,12 @@ min_style command :h3
min_style style :pre
style = {cg} or {hftn} or {sd} or {quickmin} or {fire} :ul
style = {cg} or {hftn} or {sd} or {quickmin} or {fire} or {spin} :ul
[Examples:]
min_style cg
min_style spin
min_style fire :pre
[Description:]
@ -61,6 +62,10 @@ the velocity non-parallel to the current force vector. The velocity
of each atom is initialized to 0.0 by this style, at the beginning of
a minimization.
Style {spin} is a damped spin dynamics with an adaptive
timestep.
See the "min/spin"_min_spin.html doc page for more information.
Either the {quickmin} and {fire} styles are useful in the context of
nudged elastic band (NEB) calculations via the "neb"_neb.html command.

7
doc/src/minimize.txt
View File

@ -103,6 +103,13 @@ the line search fails because the step distance backtracks to 0.0
the number of outer iterations or timesteps exceeds {maxiter}
the number of total force evaluations exceeds {maxeval} :ul
NOTE: the "minimization style"_min_style.html {spin} replaces
the force tolerance {ftol} by a torque tolerance.
The minimization procedure stops if the 2-norm (length) of the
global torque vector (defined as the cross product between the
spins and their precession vectors omega) is less than {ftol},
or if any of the other criteria are met.
NOTE: You can also use the "fix halt"_fix_halt.html command to specify
a general criterion for exiting a minimization, that is a calculation
performed on the state of the current system, as defined by an

55
examples/SPIN/spinmin/in.spinmin.bfo Normal file
View File

@ -0,0 +1,55 @@
# bfo in a 3d periodic box
units metal
dimension 3
boundary p p f
atom_style spin
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice sc 3.96
region box block 0.0 34.0 0.0 34.0 0.0 1.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for bcc iron
mass 1 1.0
set group all spin/random 11 2.50
pair_style hybrid/overlay spin/exchange 6.0 spin/magelec 4.5 spin/dmi 4.5
pair_coeff * * spin/exchange exchange 6.0 -0.01575 0.0 1.965
#pair_coeff * * spin/magelec magelec 4.5 0.000109 1.0 1.0 1.0
pair_coeff * * spin/magelec magelec 4.5 0.00109 1.0 1.0 1.0
pair_coeff * * spin/dmi dmi 4.5 0.00005 1.0 1.0 1.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
#fix 1 all precession/spin zeeman 0.001 0.0 0.0 1.0 anisotropy 0.01 1.0 0.0 0.0
fix 1 all precession/spin anisotropy 0.0000033 0.0 0.0 1.0
fix_modify 1 energy yes
timestep 0.0001
compute out_mag all spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo 50
thermo_style custom step time v_magnorm v_emag v_tmag temp etotal
thermo_modify format float %20.15g
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 1 all custom 50 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]
min_style spin
min_modify alpha_damp 1.0 discrete_factor 10.0
minimize 1.0e-10 0.0 1000 100

55
examples/SPIN/spinmin/in.spinmin.iron Normal file
View File

@ -0,0 +1,55 @@
# bcc iron in a 3d periodic box
units metal
dimension 3
boundary p p f
atom_style spin
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice bcc 2.8665
region box block 0.0 4.0 0.0 4.0 0.0 4.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for bcc iron
mass 1 55.845
set group all spin/random 31 2.2
#set group all spin 2.2 1.0 1.0 -1.0
pair_style spin/exchange 3.5
pair_coeff * * exchange 3.4 0.02726 0.2171 1.841
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
#fix 1 all precession/spin zeeman 0.001 0.0 0.0 1.0 anisotropy 0.01 1.0 0.0 0.0
fix 1 all precession/spin anisotropy 0.0001 0.0 0.0 1.0
fix_modify 1 energy yes
timestep 0.0001
compute out_mag all spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magx equal c_out_mag[1]
variable magy equal c_out_mag[2]
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo 100
thermo_style custom step time v_magx v_magz v_magnorm v_tmag etotal
thermo_modify format float %20.15g
compute outsp all property/atom spx spy spz sp fmx fmy fmz
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]
min_style spin
min_modify alpha_damp 1.0 discrete_factor 10.0
minimize 1.0e-10 1.0e-10 100000 1000

1
src/SPIN/README
View File

@ -9,6 +9,7 @@ atom in the system
* implementing magnetic pair interactions and magnetic forces
* thermostating and applying a transverse damping to the magnetic spins
* computing and outputing magnetic quantities
* minimizing the energy or total torque of a magnetic system
The different options provided by this package are explained in the
LAMMPS documentation.

16
src/SPIN/compute_spin.cpp
View File

@ -84,7 +84,7 @@ void ComputeSpin::compute_vector()
invoked_vector = update->ntimestep;
countsp = countsptot = 0.0;
countsp = countsptot = 0.0;
mag[0] = mag[1] = mag[2] = mag[3] = 0.0;
magtot[0] = magtot[1] = magtot[2] = magtot[3] = 0.0;
magenergy = magenergytot = 0.0;
@ -96,7 +96,7 @@ void ComputeSpin::compute_vector()
double **sp = atom->sp;
double **fm = atom->fm;
double tx,ty,tz;
int nlocal = atom->nlocal;
// compute total magnetization and magnetic energy
@ -105,16 +105,16 @@ void ComputeSpin::compute_vector()
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (atom->sp_flag) {
mag[0] += sp[i][0];
mag[1] += sp[i][1];
mag[2] += sp[i][2];
magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
mag[0] += sp[i][0];
mag[1] += sp[i][1];
mag[2] += sp[i][2];
magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0];
tempnum += tx*tx+ty*ty+tz*tz;
tempdenom += sp[i][0]*fm[i][0]+fm[i][1]*sp[i][1]+sp[i][2]*fm[i][2];
countsp++;
tempdenom += sp[i][0]*fm[i][0]+fm[i][1]*sp[i][1]+sp[i][2]*fm[i][2];
countsp++;
}
}
else error->all(FLERR,"Compute compute/spin requires atom/spin style");

6
src/SPIN/fix_langevin_spin.cpp
View File

@ -131,8 +131,8 @@ void FixLangevinSpin::init()
gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t));
dts = update->dt;
double hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
double kb = force->boltz; // eV/K
double hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
double kb = force->boltz; // eV/K
D = (MY_2PI*alpha_t*gil_factor*kb*temp);
D /= (hbar*dts);
sigma = sqrt(2.0*D);
@ -158,7 +158,7 @@ void FixLangevinSpin::add_tdamping(double spi[3], double fmi[3])
double cpx = fmi[1]*spi[2] - fmi[2]*spi[1];
double cpy = fmi[2]*spi[0] - fmi[0]*spi[2];
double cpz = fmi[0]*spi[1] - fmi[1]*spi[0];
// adding the transverse damping
fmi[0] -= alpha_t*cpx;

82
src/SPIN/fix_nve_spin.cpp
View File

@ -183,7 +183,7 @@ void FixNVESpin::init()
npairs = pair->instance_total;
for (int i = 0; i<npairs; i++) {
if (force->pair_match("spin",0,i)) {
npairspin ++;
npairspin ++;
}
}
}
@ -203,8 +203,8 @@ void FixNVESpin::init()
} else if (npairspin > 1) {
for (int i = 0; i<npairs; i++) {
if (force->pair_match("spin",0,i)) {
spin_pairs[count] = (PairSpin *) force->pair_match("spin",0,i);
count++;
spin_pairs[count] = (PairSpin *) force->pair_match("spin",0,i);
count++;
}
}
}
@ -264,8 +264,8 @@ void FixNVESpin::init()
void FixNVESpin::initial_integrate(int /*vflag*/)
{
double dtfm;
double **x = atom->x;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
@ -291,32 +291,32 @@ void FixNVESpin::initial_integrate(int /*vflag*/)
// update half s for all atoms
if (sector_flag) { // sectoring seq. update
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
if (sector_flag) { // sectoring seq. update
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
comm->forward_comm();
int i = stack_foot[j];
while (i >= 0) {
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
i = forward_stacks[i];
AdvanceSingleSpin(i);
i = forward_stacks[i];
}
}
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
comm->forward_comm();
int i = stack_head[j];
while (i >= 0) {
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
i = backward_stacks[i];
AdvanceSingleSpin(i);
i = backward_stacks[i];
}
}
} else if (sector_flag == 0) { // serial seq. update
comm->forward_comm(); // comm. positions of ghost atoms
for (int i = 0; i < nlocal; i++){ // advance quarter s for nlocal
} else if (sector_flag == 0) { // serial seq. update
comm->forward_comm(); // comm. positions of ghost atoms
for (int i = 0; i < nlocal; i++){ // advance quarter s for nlocal
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
}
for (int i = nlocal-1; i >= 0; i--){ // advance quarter s for nlocal
for (int i = nlocal-1; i >= 0; i--){ // advance quarter s for nlocal
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
}
@ -336,32 +336,32 @@ void FixNVESpin::initial_integrate(int /*vflag*/)
// update half s for all particles
if (sector_flag) { // sectoring seq. update
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
if (sector_flag) { // sectoring seq. update
for (int j = 0; j < nsectors; j++) { // advance quarter s for nlocal
comm->forward_comm();
int i = stack_foot[j];
while (i >= 0) {
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
i = forward_stacks[i];
AdvanceSingleSpin(i);
i = forward_stacks[i];
}
}
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
for (int j = nsectors-1; j >= 0; j--) { // advance quarter s for nlocal
comm->forward_comm();
int i = stack_head[j];
while (i >= 0) {
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
i = backward_stacks[i];
AdvanceSingleSpin(i);
i = backward_stacks[i];
}
}
} else if (sector_flag == 0) { // serial seq. update
comm->forward_comm(); // comm. positions of ghost atoms
for (int i = 0; i < nlocal; i++){ // advance quarter s for nlocal-1
} else if (sector_flag == 0) { // serial seq. update
comm->forward_comm(); // comm. positions of ghost atoms
for (int i = 0; i < nlocal; i++){ // advance quarter s for nlocal-1
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
}
for (int i = nlocal-1; i >= 0; i--){ // advance quarter s for nlocal-1
for (int i = nlocal-1; i >= 0; i--){ // advance quarter s for nlocal-1
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
}
@ -384,10 +384,10 @@ void FixNVESpin::setup_pre_neighbor()
void FixNVESpin::pre_neighbor()
{
double **x = atom->x;
double **x = atom->x;
int nlocal = atom->nlocal;
if (nlocal_max < nlocal) { // grow linked lists if necessary
if (nlocal_max < nlocal) { // grow linked lists if necessary
nlocal_max = nlocal;
backward_stacks = memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");
forward_stacks = memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");
@ -399,7 +399,7 @@ void FixNVESpin::pre_neighbor()
}
int nseci;
for (int j = 0; j < nsectors; j++) { // stacking backward order
for (int j = 0; j < nsectors; j++) { // stacking backward order
for (int i = 0; i < nlocal; i++) {
nseci = coords2sector(x[i]);
if (j != nseci) continue;
@ -407,7 +407,7 @@ void FixNVESpin::pre_neighbor()
stack_head[j] = i;
}
}
for (int j = nsectors-1; j >= 0; j--) { // stacking forward order
for (int j = nsectors-1; j >= 0; j--) { // stacking forward order
for (int i = nlocal-1; i >= 0; i--) {
nseci = coords2sector(x[i]);
if (j != nseci) continue;
@ -453,11 +453,11 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
// update langevin damping and random force
if (maglangevin_flag) { // mag. langevin
if (tdamp_flag) { // transverse damping
if (maglangevin_flag) { // mag. langevin
if (tdamp_flag) { // transverse damping
locklangevinspin->add_tdamping(spi,fmi);
}
if (temp_flag) { // spin temperature
if (temp_flag) { // spin temperature
locklangevinspin->add_temperature(fmi);
}
}
@ -567,7 +567,7 @@ void FixNVESpin::AdvanceSingleSpin(int i)
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]);
energy = (sp[i][0]*fm[i][0])+(sp[i][1]*fm[i][1])+(sp[i][2]*fm[i][2]);
dts2 = dts*dts;
dts2 = dts*dts;
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];
@ -576,18 +576,18 @@ void FixNVESpin::AdvanceSingleSpin(int i)
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*dts2;
g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
g[0] /= (1+0.25*fm2*dts2);
g[1] /= (1+0.25*fm2*dts2);
g[2] /= (1+0.25*fm2*dts2);
sp[i][0] = g[0];
sp[i][1] = g[1];
sp[i][2] = g[2];
sp[i][2] = g[2];
// renormalization (check if necessary)
@ -616,9 +616,9 @@ void FixNVESpin::AdvanceSingleSpin(int i)
/* ---------------------------------------------------------------------- */
void FixNVESpin::final_integrate()
{
{
double dtfm;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;

9
src/SPIN/fix_precession_spin.cpp
View File

@ -170,7 +170,14 @@ void FixPrecessionSpin::setup(int vflag)
/* ---------------------------------------------------------------------- */
void FixPrecessionSpin::post_force(int /*vflag*/)
void FixPrecessionSpin::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixPrecessionSpin::post_force(int /* vflag */)
{
// update mag field with time (potential improvement)

1
src/SPIN/fix_precession_spin.h
View File

@ -33,6 +33,7 @@ class FixPrecessionSpin : public Fix {
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);

333
src/SPIN/min_spin.cpp Normal file
View File

@ -0,0 +1,333 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------
Contributing authors: Julien Tranchida (SNL)
Please cite the related publication:
------------------------------------------------------------------------- */
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "min_spin.h"
#include "universe.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "output.h"
#include "timer.h"
#include "error.h"
#include "modify.h"
#include "math_special.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
// EPS_ENERGY = minimum normalization for energy tolerance
#define EPS_ENERGY 1.0e-8
#define DELAYSTEP 5
/* ---------------------------------------------------------------------- */
MinSpin::MinSpin(LAMMPS *lmp) : Min(lmp) {}
/* ---------------------------------------------------------------------- */
void MinSpin::init()
{
alpha_damp = 1.0;
discrete_factor = 10.0;
Min::init();
dts = dt = update->dt;
last_negative = update->ntimestep;
}
/* ---------------------------------------------------------------------- */
void MinSpin::setup_style()
{
double **v = atom->v;
int nlocal = atom->nlocal;
// check if the atom/spin style is defined
if (!atom->sp_flag)
error->all(FLERR,"min/spin requires atom/spin style");
for (int i = 0; i < nlocal; i++)
v[i][0] = v[i][1] = v[i][2] = 0.0;
}
/* ---------------------------------------------------------------------- */
int MinSpin::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"alpha_damp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
alpha_damp = force->numeric(FLERR,arg[1]);
return 2;
}
if (strcmp(arg[0],"discrete_factor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
discrete_factor = force->numeric(FLERR,arg[1]);
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
set current vector lengths and pointers
called after atoms have migrated
------------------------------------------------------------------------- */
void MinSpin::reset_vectors()
{
// atomic dof
// size sp is 4N vector
nvec = 4 * atom->nlocal;
if (nvec) spvec = atom->sp[0];
nvec = 3 * atom->nlocal;
if (nvec) fmvec = atom->fm[0];
if (nvec) xvec = atom->x[0];
if (nvec) fvec = atom->f[0];
}
/* ----------------------------------------------------------------------
minimization via damped spin dynamics
------------------------------------------------------------------------- */
int MinSpin::iterate(int maxiter)
{
bigint ntimestep;
double fmdotfm;
int flag,flagall;
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
energy_force(0);
dts = evaluate_dt();
// apply damped precessional dynamics to the spins
advance_spins(dts);
eprevious = ecurrent;
ecurrent = energy_force(0);
neval++;
//// energy tolerance criterion
//// only check after DELAYSTEP elapsed since velocties reset to 0
//// sync across replicas if running multi-replica minimization
if (update->etol > 0.0 && ntimestep-last_negative > DELAYSTEP) {
if (update->multireplica == 0) {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
return ETOL;
} else {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return ETOL;
}
}
// magnetic torque tolerance criterion
// sync across replicas if running multi-replica minimization
if (update->ftol > 0.0) {
fmdotfm = fmnorm_sqr();
if (update->multireplica == 0) {
if (fmdotfm < update->ftol*update->ftol) return FTOL;
} else {
if (fmdotfm < update->ftol*update->ftol) flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return FTOL;
}
}
// output for thermo, dump, restart files
if (output->next == ntimestep) {
timer->stamp();
output->write(ntimestep);
timer->stamp(Timer::OUTPUT);
}
}
return MAXITER;
}
/* ----------------------------------------------------------------------
evaluate max timestep
---------------------------------------------------------------------- */
double MinSpin::evaluate_dt()
{
double dtmax;
double fmsq;
double fmaxsqone,fmaxsqloc,fmaxsqall;
int nlocal = atom->nlocal;
double **fm = atom->fm;
// finding max fm on this proc.
fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
for (int i = 0; i < nlocal; i++) {
fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];
fmaxsqone = MAX(fmaxsqone,fmsq);
}
// finding max fm on this replica
fmaxsqloc = fmaxsqone;
MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
// finding max fm over all replicas, if necessary
// this communicator would be invalid for multiprocess replicas
fmaxsqall = fmaxsqloc;
if (update->multireplica == 1) {
fmaxsqall = fmaxsqloc;
MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
}
if (fmaxsqall == 0.0)
error->all(FLERR,"Incorrect fmaxsqall calculation");
// define max timestep by dividing by the
// inverse of max frequency by discrete_factor
dtmax = MY_2PI/(discrete_factor*sqrt(fmaxsqall));
return dtmax;
}
/* ----------------------------------------------------------------------
geometric damped advance of spins
---------------------------------------------------------------------- */
void MinSpin::advance_spins(double dts)
{
int nlocal = atom->nlocal;
double **sp = atom->sp;
double **fm = atom->fm;
double tdampx,tdampy,tdampz;
double msq,scale,fm2,energy,dts2;
double cp[3],g[3];
dts2 = dts*dts;
// loop on all spins on proc.
for (int i = 0; i < nlocal; i++) {
// calc. damping torque
tdampx = -alpha_damp*(fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
tdampy = -alpha_damp*(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);
tdampz = -alpha_damp*(fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]);
// apply advance algorithm (geometric, norm preserving)
fm2 = (tdampx*tdampx+tdampy*tdampy+tdampz*tdampz);
energy = (sp[i][0]*tdampx)+(sp[i][1]*tdampy)+(sp[i][2]*tdampz);
cp[0] = tdampy*sp[i][2]-tdampz*sp[i][1];
cp[1] = tdampz*sp[i][0]-tdampx*sp[i][2];
cp[2] = tdampx*sp[i][1]-tdampy*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] += (tdampx*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
g[1] += (tdampy*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
g[2] += (tdampz*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
g[0] /= (1+0.25*fm2*dts2);
g[1] /= (1+0.25*fm2*dts2);
g[2] /= (1+0.25*fm2*dts2);
sp[i][0] = g[0];
sp[i][1] = g[1];
sp[i][2] = g[2];
// renormalization (check if 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;
// no comm. to atoms with same tag
// because no need for simplecticity
}
}
/* ----------------------------------------------------------------------
compute and return ||mag. torque||_2^2
------------------------------------------------------------------------- */
double MinSpin::fmnorm_sqr()
{
int nlocal = atom->nlocal;
double tx,ty,tz;
double **sp = atom->sp;
double **fm = atom->fm;
// calc. magnetic torques
double local_norm2_sqr = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
ty = (fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);
tz = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]);
local_norm2_sqr += tx*tx + ty*ty + tz*tz;
}
// no extra atom calc. for spins
if (nextra_atom)
error->all(FLERR,"extra atom option not available yet");
double norm2_sqr = 0.0;
MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world);
return norm2_sqr;
}

59
src/SPIN/min_spin.h Normal file
View File

@ -0,0 +1,59 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef MINIMIZE_CLASS
MinimizeStyle(spin,MinSpin)
#else
#ifndef LMP_MIN_SPIN_H
#define LMP_MIN_SPIN_H
#include "min.h"
namespace LAMMPS_NS {
class MinSpin : public Min {
public:
MinSpin(class LAMMPS *);
~MinSpin() {}
void init();
void setup_style();
int modify_param(int, char **);
void reset_vectors();
int iterate(int);
double evaluate_dt();
void advance_spins(double);
double fmnorm_sqr();
private:
// global and spin timesteps
double dt;
double dts;
double alpha_damp; // damping for spin minimization
double discrete_factor; // factor for spin timestep evaluation
double *spvec; // variables for atomic dof, as 1d vector
double *fmvec; // variables for atomic dof, as 1d vector
bigint last_negative;
};
}
#endif
#endif

52
src/SPIN/pair_spin_dmi.cpp
View File

@ -173,8 +173,8 @@ void PairSpinDmi::init_style()
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
ifix++;
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
if ((ifix == modify->nfix) && (comm->me == 0))
error->warning(FLERR,"Using pair/spin style without nve/spin");
// get the lattice_flag from nve/spin
@ -238,7 +238,7 @@ void PairSpinDmi::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **sp = atom->sp;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
@ -295,32 +295,32 @@ void PairSpinDmi::compute(int eflag, int vflag)
// compute dmi interaction
if (rsq <= local_cut2) {
compute_dmi(i,j,eij,fmi,spj);
if (lattice_flag) {
compute_dmi_mech(i,j,rsq,eij,fi,spi,spj);
}
compute_dmi(i,j,eij,fmi,spj);
if (lattice_flag) {
compute_dmi_mech(i,j,rsq,eij,fi,spi,spj);
}
}
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
@ -342,7 +342,7 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
double delx,dely,delz;
double spj[3];
int i,j,jnum,itype,jtype,ntypes;
int j,jnum,itype,jtype,ntypes;
int k,locflag;
int *jlist,*numneigh,**firstneigh;
@ -350,7 +350,7 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
itype = type[ii];
@ -360,20 +360,20 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
// if interaction applies to type ii,
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
//i = ilist[ii];
@ -422,7 +422,7 @@ void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double
{
int *type = atom->type;
int itype, jtype;
double dmix, dmiy, dmiz;
double dmix, dmiy, dmiz;
itype = type[i];
jtype = type[j];
@ -444,7 +444,7 @@ void PairSpinDmi::compute_dmi_mech(int i, int j, double rsq, double /*eij*/[3],
{
int *type = atom->type;
int itype, jtype;
double dmix,dmiy,dmiz;
double dmix,dmiy,dmiz;
itype = type[i];
jtype = type[j];
double csx,csy,csz,cdmx,cdmy,cdmz,irij;
@ -509,7 +509,7 @@ void PairSpinDmi::write_restart(FILE *fp)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&DM[i][j],sizeof(double),1,fp);
fwrite(&DM[i][j],sizeof(double),1,fp);
fwrite(&v_dmx[i][j],sizeof(double),1,fp);
fwrite(&v_dmy[i][j],sizeof(double),1,fp);
fwrite(&v_dmz[i][j],sizeof(double),1,fp);

46
src/SPIN/pair_spin_exchange.cpp
View File

@ -160,8 +160,8 @@ void PairSpinExchange::init_style()
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
ifix++;
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
if ((ifix == modify->nfix) && (comm->me == 0))
error->warning(FLERR,"Using pair/spin style without nve/spin");
// get the lattice_flag from nve/spin
@ -223,7 +223,7 @@ void PairSpinExchange::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **sp = atom->sp;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
@ -278,32 +278,32 @@ void PairSpinExchange::compute(int eflag, int vflag)
// compute exchange interaction
if (rsq <= local_cut2) {
compute_exchange(i,j,rsq,fmi,spj);
compute_exchange(i,j,rsq,fmi,spj);
if (lattice_flag) {
compute_exchange_mech(i,j,rsq,eij,fi,spi,spj);
}
compute_exchange_mech(i,j,rsq,eij,fi,spi,spj);
}
}
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
@ -325,7 +325,7 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3])
double delx,dely,delz;
double spj[3];
int i,j,jnum,itype,jtype,ntypes;
int j,jnum,itype,jtype,ntypes;
int k,locflag;
int *jlist,*numneigh,**firstneigh;
@ -343,24 +343,24 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3])
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
// if interaction applies to type ii,
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
if (locflag == 1) {
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
@ -388,7 +388,7 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3])
compute_exchange(ii,j,rsq,fmi,spj);
}
}
}
}
}
/* ----------------------------------------------------------------------

32
src/SPIN/pair_spin_magelec.cpp
View File

@ -166,8 +166,8 @@ void PairSpinMagelec::init_style()
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
ifix++;
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
if ((ifix == modify->nfix) && (comm->me == 0))
error->warning(FLERR,"Using pair/spin style without nve/spin");
// get the lattice_flag from nve/spin
@ -332,7 +332,7 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
double delx,dely,delz;
double spj[3];
int i,j,jnum,itype,jtype,ntypes;
int j,jnum,itype,jtype,ntypes;
int k,locflag;
int *jlist,*numneigh,**firstneigh;
@ -340,7 +340,7 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
itype = type[ii];
@ -350,42 +350,42 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
// if interaction applies to type ii,
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
if (locflag == 1) {
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
jlist = firstneigh[ii];
jnum = numneigh[ii];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
local_cut2 = cut_spin_magelec[itype][jtype]*cut_spin_magelec[itype][jtype];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
delx = xi[0] - x[j][0];
dely = xi[1] - x[j][1];
delz = xi[2] - x[j][2];
@ -394,7 +394,7 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
eij[0] = -inorm*delx;
eij[1] = -inorm*dely;
eij[2] = -inorm*delz;
if (rsq <= local_cut2) {
compute_magelec(ii,j,eij,fmi,spj);
}

18
src/SPIN/pair_spin_neel.cpp
View File

@ -173,8 +173,8 @@ void PairSpinNeel::init_style()
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
ifix++;
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
if ((ifix == modify->nfix) && (comm->me == 0))
error->warning(FLERR,"Using pair/spin style without nve/spin");
// get the lattice_flag from nve/spin
@ -343,7 +343,7 @@ void PairSpinNeel::compute_single_pair(int ii, double fmi[3])
double xi[3], rij[3], eij[3];
double spi[3], spj[3];
int i,j,jnum,itype,jtype,ntypes;
int j,jnum,itype,jtype,ntypes;
int k,locflag;
int *jlist,*numneigh,**firstneigh;
@ -351,7 +351,7 @@ void PairSpinNeel::compute_single_pair(int ii, double fmi[3])
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
itype = type[ii];
@ -361,20 +361,20 @@ void PairSpinNeel::compute_single_pair(int ii, double fmi[3])
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
// if interaction applies to type ii,
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
if (locflag == 1) {

6
src/min.cpp
View File

@ -655,7 +655,11 @@ void Min::modify_params(int narg, char **arg)
else if (strcmp(arg[iarg+1],"forcezero") == 0) linestyle = 2;
else error->all(FLERR,"Illegal min_modify command");
iarg += 2;
} else error->all(FLERR,"Illegal min_modify command");
} else {
int n = modify_param(narg-iarg,&arg[iarg]);
if (n == 0) error->all(FLERR,"Illegal fix_modify command");
iarg += n;
}
}
}

1
src/min.h
View File

@ -38,6 +38,7 @@ class Min : protected Pointers {
int request(class Pair *, int, double);
virtual bigint memory_usage() {return 0;}
void modify_params(int, char **);
virtual int modify_param(int, char **) {return 0;}
double fnorm_sqr();
double fnorm_inf();