forked from lijiext/lammps
Commit1 JT 082118
- created pppm_dipole_spin.h/cpp (child-class of pppm_dipole) - improved pair_spin_long.h/cpp - created documentation for pair_spin_long - new 3xN fm_long vector in atom_vec_spin (with associated comm)
This commit is contained in:
julient31
parent
5e287033f7
commit
8d79db03d3
Binary file not shown.
|
After Width: | Height: | Size: 12 KiB |
|
|
@ -0,0 +1,20 @@
|
|||
\documentclass[preview]{standalone}
|
||||
\usepackage{varwidth}
|
||||
\usepackage[utf8x]{inputenc}
|
||||
\usepackage{amsmath,amssymb,graphics,bm,setspace}
|
||||
|
||||
\begin{document}
|
||||
\begin{varwidth}{50in}
|
||||
\begin{equation}
|
||||
\mathcal{H}_{\rm long}=
|
||||
-\frac{\mu_{0} \left( \mu_B\right)^2}{4\pi}
|
||||
\sum_{i,j,i\neq j}^{N}
|
||||
\frac{g_i g_j}{r_{ij}^3}
|
||||
\Big(3
|
||||
\left(\bm{e}_{ij}\cdot \bm{s}_{i}\right)
|
||||
\left(\bm{e}_{ij}\cdot \bm{s}_{j}\right)
|
||||
-\bm{s}_i\cdot\bm{s}_j \Big)
|
||||
\nonumber
|
||||
\end{equation}
|
||||
\end{varwidth}
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 16 KiB |
|
|
@ -0,0 +1,23 @@
|
|||
\documentclass[preview]{standalone}
|
||||
\usepackage{varwidth}
|
||||
\usepackage[utf8x]{inputenc}
|
||||
\usepackage{amsmath,amssymb,graphics,bm,setspace}
|
||||
|
||||
\begin{document}
|
||||
\begin{varwidth}{50in}
|
||||
\begin{equation}
|
||||
\bm{F}_i =
|
||||
\frac{\mu_0 (\mu_B)^2}{4\pi} \sum_j
|
||||
\frac{g_i g_j}{r_{ij}^4}
|
||||
\Big[\big( (\bm{s}_i\cdot\bm{s}_j)
|
||||
-5(\bm{e}_{ij}\cdot\bm{s}_i)
|
||||
(\bm{e}_{ij}\cdot\bm{s}_j)\big) \bm{e}_{ij}+
|
||||
\big(
|
||||
(\bm{e}_{ij}\cdot\bm{s}_i)\bm{s}_j+
|
||||
(\bm{e}_{ij}\cdot\bm{s}_j)\bm{s}_i
|
||||
\big)
|
||||
\Big]
|
||||
\nonumber
|
||||
\end{equation}
|
||||
\end{varwidth}
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 9.2 KiB |
|
|
@ -0,0 +1,17 @@
|
|||
\documentclass[preview]{standalone}
|
||||
\usepackage{varwidth}
|
||||
\usepackage[utf8x]{inputenc}
|
||||
\usepackage{amsmath,amssymb,graphics,bm,setspace}
|
||||
|
||||
\begin{document}
|
||||
\begin{varwidth}{50in}
|
||||
\begin{equation}
|
||||
\bm{\omega}_i =
|
||||
\frac{\mu_0 (\mu_B)^2}{4\pi\hbar}\sum_{j}
|
||||
\frac{g_i g_j}{r_{ij}^3}
|
||||
\, \Big(
|
||||
3\,(\bm{e}_{ij}\cdot\bm{s}_{j})\bm{e}_{ij}
|
||||
-\bm{s}_{j} \Big) \nonumber
|
||||
\end{equation}
|
||||
\end{varwidth}
|
||||
\end{document}
|
||||
|
|
@ -0,0 +1,84 @@
|
|||
"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
|
||||
|
||||
pair_style spin/long command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
pair_style spin/long cutoff (cutoff)
|
||||
|
||||
cutoff = global cutoff pair (distance in metal units) :ulb,l
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
pair_style spin/long 10.0
|
||||
pair_coeff * * long 10.0
|
||||
pair_coeff 2 3 long 8.0 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
Style {pair/spin/long} computes interactions between pairs of particles
|
||||
that each have a magnetic spin.
|
||||
|
||||
:c,image(Eqs/pair_spin_long_range.jpg)
|
||||
|
||||
where si and sj are two magnetic spins of two particles with Lande factors
|
||||
gi and gj respectively, eij = (ri - rj)/|ri-rj| is the unit vector between
|
||||
sites i and j, mu0 the vacuum permeability, muB the Bohr magneton (muB =
|
||||
5.788 eV/T in metal units).
|
||||
|
||||
Style {pair/spin/long} computes magnetic precession vectors:
|
||||
|
||||
:c,image(Eqs/pair_spin_long_range_magforce.jpg)
|
||||
|
||||
with h the Planck constant (in metal units), and a mechanical force:
|
||||
|
||||
:c,image(Eqs/pair_spin_long_range_force.jpg)
|
||||
|
||||
|
||||
The following coefficient must be defined for each pair of atoms
|
||||
types via the "pair_coeff"_pair_coeff.html command as in the examples
|
||||
above, or in the data file or restart files read by the
|
||||
"read_data"_read_data.html or "read_restart"_read_restart.html
|
||||
commands, or by mixing as described below:
|
||||
|
||||
rc (distance units) :ul
|
||||
|
||||
with rc is the radius cutoff of the short-range component of the
|
||||
long-range interaction (see "(Cerda)"_#Cerda1 for more
|
||||
explanation).
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
The {pair/spin/long} style is part of the SPIN package. It is only
|
||||
enabled if LAMMPS was built with that package. See the
|
||||
"Making LAMMPS"_Section_start.html#start_3 section for more info.
|
||||
|
||||
The {pair/spin/long} style computes the short-range component of
|
||||
the dipole-dipole interaction. The functions evaluating the
|
||||
long-range component are part of the KSPACE package.
|
||||
They can be enabled only if LAMMPS was built with that package.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"atom_style spin"_atom_style.html, "pair_coeff"_pair_coeff.html,
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(Tranchida6)
|
||||
[(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
|
||||
Journal of Computational Physics, (2018).
|
||||
|
||||
:link(Cerda1)
|
||||
[(Cerda)] Cerda, Ballenegger, Lenz, and Holm, J Chem Phys, 129(23),
|
||||
234104 (2008).
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,213 @@
|
|||
/* -*- 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 KSPACE_CLASS
|
||||
|
||||
KSpaceStyle(pppm/dipole,PPPMDipole)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PPPM_DIPOLE_H
|
||||
#define LMP_PPPM_DIPOLE_H
|
||||
|
||||
#include "pppm.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PPPMDipole : public PPPM {
|
||||
public:
|
||||
PPPMDipole(class LAMMPS *, int, char **);
|
||||
virtual ~PPPMDipole();
|
||||
void init();
|
||||
void setup();
|
||||
void setup_grid();
|
||||
void compute(int, int);
|
||||
int timing_1d(int, double &);
|
||||
int timing_3d(int, double &);
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
void set_grid_global(double);
|
||||
//double newton_raphson_f();
|
||||
|
||||
void allocate();
|
||||
void allocate_peratom();
|
||||
void deallocate();
|
||||
void deallocate_peratom();
|
||||
void compute_gf_denom();
|
||||
|
||||
void slabcorr();
|
||||
|
||||
// grid communication
|
||||
|
||||
void pack_forward(int, FFT_SCALAR *, int, int *);
|
||||
void unpack_forward(int, FFT_SCALAR *, int, int *);
|
||||
void pack_reverse(int, FFT_SCALAR *, int, int *);
|
||||
void unpack_reverse(int, FFT_SCALAR *, int, int *);
|
||||
|
||||
// dipole
|
||||
|
||||
FFT_SCALAR ***densityx_brick_dipole,***densityy_brick_dipole,***densityz_brick_dipole;
|
||||
FFT_SCALAR ***vdxx_brick_dipole,***vdyy_brick_dipole,***vdzz_brick_dipole;
|
||||
FFT_SCALAR ***vdxy_brick_dipole,***vdxz_brick_dipole,***vdyz_brick_dipole;
|
||||
FFT_SCALAR ***ux_brick_dipole,***uy_brick_dipole,***uz_brick_dipole;
|
||||
FFT_SCALAR ***v0x_brick_dipole,***v1x_brick_dipole,***v2x_brick_dipole;
|
||||
FFT_SCALAR ***v3x_brick_dipole,***v4x_brick_dipole,***v5x_brick_dipole;
|
||||
FFT_SCALAR ***v0y_brick_dipole,***v1y_brick_dipole,***v2y_brick_dipole;
|
||||
FFT_SCALAR ***v3y_brick_dipole,***v4y_brick_dipole,***v5y_brick_dipole;
|
||||
FFT_SCALAR ***v0z_brick_dipole,***v1z_brick_dipole,***v2z_brick_dipole;
|
||||
FFT_SCALAR ***v3z_brick_dipole,***v4z_brick_dipole,***v5z_brick_dipole;
|
||||
FFT_SCALAR *work3,*work4;
|
||||
FFT_SCALAR *densityx_fft_dipole,*densityy_fft_dipole,*densityz_fft_dipole;
|
||||
class GridComm *cg_dipole;
|
||||
class GridComm *cg_peratom_dipole;
|
||||
int only_dipole_flag;
|
||||
double musum,musqsum,mu2;
|
||||
double find_gewald_dipole(double, double, bigint, double, double);
|
||||
double newton_raphson_f_dipole(double, double, bigint, double, double);
|
||||
double derivf_dipole(double, double, bigint, double, double);
|
||||
double compute_df_kspace_dipole(double);
|
||||
double compute_qopt_dipole();
|
||||
void compute_gf_dipole();
|
||||
void make_rho_dipole();
|
||||
void brick2fft_dipole();
|
||||
void poisson_ik_dipole();
|
||||
void poisson_peratom_dipole();
|
||||
void fieldforce_ik_dipole();
|
||||
void fieldforce_peratom_dipole();
|
||||
double final_accuracy_dipole(double dipole2);
|
||||
void musum_musq();
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Cannot (yet) use charges with Kspace style PPPMDipole
|
||||
|
||||
Charge-dipole interactions are not yet implemented in PPPMDipole so this
|
||||
feature is not yet supported.
|
||||
|
||||
E: Must redefine kspace_style after changing to triclinic box
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Kspace style requires atom attribute mu
|
||||
|
||||
The atom style defined does not have this attribute.
|
||||
|
||||
E: Cannot (yet) use kspace_modify diff ad with dipoles
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot (yet) use 'electron' units with dipoles
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot yet use triclinic cells with PPPMDipole
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot yet use TIP4P with PPPMDipole
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot use nonperiodic boundaries with PPPM
|
||||
|
||||
For kspace style pppm, all 3 dimensions must have periodic boundaries
|
||||
unless you use the kspace_modify command to define a 2d slab with a
|
||||
non-periodic z dimension.
|
||||
|
||||
E: Incorrect boundaries with slab PPPM
|
||||
|
||||
Must have periodic x,y dimensions and non-periodic z dimension to use
|
||||
2d slab option with PPPM.
|
||||
|
||||
E: PPPM order cannot be < 2 or > than %d
|
||||
|
||||
This is a limitation of the PPPM implementation in LAMMPS.
|
||||
|
||||
E: KSpace style is incompatible with Pair style
|
||||
|
||||
Setting a kspace style requires that a pair style with matching
|
||||
long-range dipole components be used.
|
||||
|
||||
W: Reducing PPPM order b/c stencil extends beyond nearest neighbor processor
|
||||
|
||||
This may lead to a larger grid than desired. See the kspace_modify overlap
|
||||
command to prevent changing of the PPPM order.
|
||||
|
||||
E: PPPM order < minimum allowed order
|
||||
|
||||
The default minimum order is 2. This can be reset by the
|
||||
kspace_modify minorder command.
|
||||
|
||||
E: PPPM grid stencil extends beyond nearest neighbor processor
|
||||
|
||||
This is not allowed if the kspace_modify overlap setting is no.
|
||||
|
||||
E: KSpace accuracy must be > 0
|
||||
|
||||
The kspace accuracy designated in the input must be greater than zero.
|
||||
|
||||
E: Could not compute grid size
|
||||
|
||||
The code is unable to compute a grid size consistent with the desired
|
||||
accuracy. This error should not occur for typical problems. Please
|
||||
send an email to the developers.
|
||||
|
||||
E: PPPM grid is too large
|
||||
|
||||
The global PPPM grid is larger than OFFSET in one or more dimensions.
|
||||
OFFSET is currently set to 4096. You likely need to decrease the
|
||||
requested accuracy.
|
||||
|
||||
E: Could not compute g_ewald
|
||||
|
||||
The Newton-Raphson solver failed to converge to a good value for
|
||||
g_ewald. This error should not occur for typical problems. Please
|
||||
send an email to the developers.
|
||||
|
||||
E: Non-numeric box dimensions - simulation unstable
|
||||
|
||||
The box size has apparently blown up.
|
||||
|
||||
E: Out of range atoms - cannot compute PPPM
|
||||
|
||||
One or more atoms are attempting to map their charge to a PPPM grid
|
||||
point that is not owned by a processor. This is likely for one of two
|
||||
reasons, both of them bad. First, it may mean that an atom near the
|
||||
boundary of a processor's sub-domain has moved more than 1/2 the
|
||||
"neighbor skin distance"_neighbor.html without neighbor lists being
|
||||
rebuilt and atoms being migrated to new processors. This also means
|
||||
you may be missing pairwise interactions that need to be computed.
|
||||
The solution is to change the re-neighboring criteria via the
|
||||
"neigh_modify"_neigh_modify command. The safest settings are "delay 0
|
||||
every 1 check yes". Second, it may mean that an atom has moved far
|
||||
outside a processor's sub-domain or even the entire simulation box.
|
||||
This indicates bad physics, e.g. due to highly overlapping atoms, too
|
||||
large a timestep, etc.
|
||||
|
||||
E: Using kspace solver PPPMDipole on system with no dipoles
|
||||
|
||||
Must have non-zero dipoles with PPPMDipole.
|
||||
|
||||
E: Must use kspace_modify gewald for system with no dipoles
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,750 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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 author: Stan Moore (SNL)
|
||||
Julien Tranchida (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <mpi.h>
|
||||
#include <cstring>
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <cmath>
|
||||
#include "pppm_dipole_spin.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "gridcomm.h"
|
||||
#include "neighbor.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "bond.h"
|
||||
#include "angle.h"
|
||||
#include "domain.h"
|
||||
#include "fft3d_wrap.h"
|
||||
#include "remap_wrap.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "update.h"
|
||||
|
||||
#include "math_const.h"
|
||||
#include "math_special.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using namespace MathSpecial;
|
||||
|
||||
#define MAXORDER 7
|
||||
#define OFFSET 16384
|
||||
#define LARGE 10000.0
|
||||
#define SMALL 0.00001
|
||||
#define EPS_HOC 1.0e-7
|
||||
|
||||
enum{REVERSE_SP};
|
||||
enum{FORWARD_SP,FORWARD_SP_PERATOM};
|
||||
|
||||
#ifdef FFT_SINGLE
|
||||
#define ZEROF 0.0f
|
||||
#define ONEF 1.0f
|
||||
#else
|
||||
#define ZEROF 0.0
|
||||
#define ONEF 1.0
|
||||
#endif
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp, int narg, char **arg) :
|
||||
PPPMDipole(lmp, narg, arg)
|
||||
{
|
||||
dipoleflag = 0;
|
||||
spinflag = 1;
|
||||
group_group_enable = 0;
|
||||
|
||||
cg_dipole = NULL;
|
||||
cg_peratom_dipole = NULL;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
free all memory
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
PPPMDipoleSpin::~PPPMDipoleSpin()
|
||||
{
|
||||
if (copymode) return;
|
||||
|
||||
deallocate();
|
||||
if (peratom_allocate_flag) deallocate_peratom();
|
||||
fft1 = NULL;
|
||||
fft2 = NULL;
|
||||
remap = NULL;
|
||||
cg_dipole = NULL;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
called once before run
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::init()
|
||||
{
|
||||
if (me == 0) {
|
||||
if (screen) fprintf(screen,"PPPMDipoleSpin initialization ...\n");
|
||||
if (logfile) fprintf(logfile,"PPPMDipoleSpin initialization ...\n");
|
||||
}
|
||||
|
||||
// error check
|
||||
|
||||
spinflag = atom->sp?1:0;
|
||||
|
||||
triclinic_check();
|
||||
|
||||
if (triclinic != domain->triclinic)
|
||||
error->all(FLERR,"Must redefine kspace_style after changing to triclinic box");
|
||||
|
||||
if (domain->dimension == 2) error->all(FLERR,
|
||||
"Cannot use PPPMDipoleSpin with 2d simulation");
|
||||
if (comm->style != 0)
|
||||
error->universe_all(FLERR,"PPPMDipoleSpin can only currently be used with "
|
||||
"comm_style brick");
|
||||
|
||||
if (!atom->sp) error->all(FLERR,"Kspace style requires atom attribute sp");
|
||||
|
||||
if (atom->sp && differentiation_flag == 1) error->all(FLERR,"Cannot (yet) use kspace_modify diff"
|
||||
" ad with spins");
|
||||
|
||||
if (spinflag && strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"'metal' units have to be used with spins");
|
||||
|
||||
if (slabflag == 0 && domain->nonperiodic > 0)
|
||||
error->all(FLERR,"Cannot use nonperiodic boundaries with PPPMDipoleSpin");
|
||||
if (slabflag) {
|
||||
if (domain->xperiodic != 1 || domain->yperiodic != 1 ||
|
||||
domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1)
|
||||
error->all(FLERR,"Incorrect boundaries with slab PPPMDipoleSpin");
|
||||
}
|
||||
|
||||
if (order < 2 || order > MAXORDER) {
|
||||
char str[128];
|
||||
sprintf(str,"PPPMDipoleSpin order cannot be < 2 or > than %d",MAXORDER);
|
||||
error->all(FLERR,str);
|
||||
}
|
||||
|
||||
// extract short-range Coulombic cutoff from pair style
|
||||
|
||||
triclinic = domain->triclinic;
|
||||
if (triclinic)
|
||||
error->all(FLERR,"Cannot yet use triclinic cells with PPPMDipoleSpin");
|
||||
|
||||
pair_check();
|
||||
|
||||
int itmp = 0;
|
||||
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
|
||||
if (p_cutoff == NULL)
|
||||
error->all(FLERR,"KSpace style is incompatible with Pair style");
|
||||
cutoff = *p_cutoff;
|
||||
|
||||
// kspace TIP4P not yet supported
|
||||
|
||||
if (tip4pflag)
|
||||
error->all(FLERR,"Cannot yet use TIP4P with PPPMDipoleSpin");
|
||||
|
||||
scale = 1.0;
|
||||
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
|
||||
mub = 5.78901e-5; // in eV/T
|
||||
mu_0 = 1.2566370614e-6; // in T.m/A
|
||||
mub2mu0 = mub * mub * mu_0; // in eV
|
||||
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
|
||||
spsum_spsq();
|
||||
natoms_original = atom->natoms;
|
||||
|
||||
// set accuracy (force units) from accuracy_relative or accuracy_absolute
|
||||
|
||||
// is two_charge_force still relevant for spin systems?
|
||||
|
||||
if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute;
|
||||
else accuracy = accuracy_relative * two_charge_force;
|
||||
|
||||
// free all arrays previously allocated
|
||||
|
||||
deallocate();
|
||||
if (peratom_allocate_flag) deallocate_peratom();
|
||||
|
||||
// setup FFT grid resolution and g_ewald
|
||||
// normally one iteration thru while loop is all that is required
|
||||
// if grid stencil does not extend beyond neighbor proc
|
||||
// or overlap is allowed, then done
|
||||
// else reduce order and try again
|
||||
|
||||
int (*procneigh)[2] = comm->procneigh;
|
||||
|
||||
GridComm *cgtmp = NULL;
|
||||
int iteration = 0;
|
||||
|
||||
while (order >= minorder) {
|
||||
if (iteration && me == 0)
|
||||
error->warning(FLERR,"Reducing PPPMDipoleSpin order b/c stencil extends "
|
||||
"beyond nearest neighbor processor");
|
||||
|
||||
compute_gf_denom();
|
||||
set_grid_global(sp2);
|
||||
set_grid_local();
|
||||
if (overlap_allowed) break;
|
||||
|
||||
cgtmp = new GridComm(lmp,world,1,1,
|
||||
nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in,
|
||||
nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out,
|
||||
procneigh[0][0],procneigh[0][1],procneigh[1][0],
|
||||
procneigh[1][1],procneigh[2][0],procneigh[2][1]);
|
||||
cgtmp->ghost_notify();
|
||||
if (!cgtmp->ghost_overlap()) break;
|
||||
delete cgtmp;
|
||||
|
||||
order--;
|
||||
iteration++;
|
||||
}
|
||||
|
||||
if (order < minorder) error->all(FLERR,"PPPMDipoleSpin order < minimum allowed order");
|
||||
if (!overlap_allowed && cgtmp->ghost_overlap())
|
||||
error->all(FLERR,"PPPMDipoleSpin grid stencil extends "
|
||||
"beyond nearest neighbor processor");
|
||||
if (cgtmp) delete cgtmp;
|
||||
|
||||
// adjust g_ewald
|
||||
|
||||
if (!gewaldflag) adjust_gewald();
|
||||
|
||||
// calculate the final accuracy
|
||||
|
||||
double estimated_accuracy = final_accuracy_dipole(sp2);
|
||||
|
||||
// print stats
|
||||
|
||||
int ngrid_max,nfft_both_max;
|
||||
MPI_Allreduce(&ngrid,&ngrid_max,1,MPI_INT,MPI_MAX,world);
|
||||
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
|
||||
|
||||
if (me == 0) {
|
||||
|
||||
#ifdef FFT_SINGLE
|
||||
const char fft_prec[] = "single";
|
||||
#else
|
||||
const char fft_prec[] = "double";
|
||||
#endif
|
||||
|
||||
if (screen) {
|
||||
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
|
||||
fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
|
||||
fprintf(screen," stencil order = %d\n",order);
|
||||
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
|
||||
estimated_accuracy);
|
||||
fprintf(screen," estimated relative force accuracy = %g\n",
|
||||
estimated_accuracy/two_charge_force);
|
||||
fprintf(screen," using %s precision FFTs\n",fft_prec);
|
||||
fprintf(screen," 3d grid and FFT values/proc = %d %d\n",
|
||||
ngrid_max,nfft_both_max);
|
||||
}
|
||||
if (logfile) {
|
||||
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
|
||||
fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
|
||||
fprintf(logfile," stencil order = %d\n",order);
|
||||
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
|
||||
estimated_accuracy);
|
||||
fprintf(logfile," estimated relative force accuracy = %g\n",
|
||||
estimated_accuracy/two_charge_force);
|
||||
fprintf(logfile," using %s precision FFTs\n",fft_prec);
|
||||
fprintf(logfile," 3d grid and FFT values/proc = %d %d\n",
|
||||
ngrid_max,nfft_both_max);
|
||||
}
|
||||
}
|
||||
|
||||
// allocate K-space dependent memory
|
||||
// don't invoke allocate peratom(), will be allocated when needed
|
||||
|
||||
allocate();
|
||||
cg_dipole->ghost_notify();
|
||||
cg_dipole->setup();
|
||||
|
||||
// pre-compute Green's function denominator expansion
|
||||
// pre-compute 1d charge distribution coefficients
|
||||
|
||||
compute_gf_denom();
|
||||
compute_rho_coeff();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute the PPPMDipoleSpin long-range force, energy, virial
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j;
|
||||
|
||||
// set energy/virial flags
|
||||
// invoke allocate_peratom() if needed for first time
|
||||
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = evflag_atom = eflag_global = vflag_global =
|
||||
eflag_atom = vflag_atom = 0;
|
||||
|
||||
if (evflag_atom && !peratom_allocate_flag) {
|
||||
allocate_peratom();
|
||||
cg_peratom_dipole->ghost_notify();
|
||||
cg_peratom_dipole->setup();
|
||||
}
|
||||
|
||||
// if atom count has changed, update qsum and qsqsum
|
||||
|
||||
if (atom->natoms != natoms_original) {
|
||||
spsum_spsq();
|
||||
natoms_original = atom->natoms;
|
||||
}
|
||||
|
||||
// return if there are no spins
|
||||
|
||||
if (spsqsum == 0.0) return;
|
||||
|
||||
// convert atoms from box to lamda coords
|
||||
|
||||
boxlo = domain->boxlo;
|
||||
|
||||
// extend size of per-atom arrays if necessary
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(part2grid);
|
||||
nmax = atom->nmax;
|
||||
memory->create(part2grid,nmax,3,"pppm_spin:part2grid");
|
||||
}
|
||||
|
||||
// find grid points for all my particles
|
||||
// map my particle charge onto my local 3d on-grid density
|
||||
|
||||
particle_map();
|
||||
make_rho_spin();
|
||||
|
||||
// all procs communicate density values from their ghost cells
|
||||
// to fully sum contribution in their 3d bricks
|
||||
// remap from 3d decomposition to FFT decomposition
|
||||
|
||||
cg_dipole->reverse_comm(this,REVERSE_SP);
|
||||
brick2fft_dipole();
|
||||
|
||||
// compute potential gradient on my FFT grid and
|
||||
// portion of e_long on this proc's FFT grid
|
||||
// return gradients (electric fields) in 3d brick decomposition
|
||||
// also performs per-atom calculations via poisson_peratom()
|
||||
|
||||
poisson_ik_dipole();
|
||||
|
||||
// all procs communicate E-field values
|
||||
// to fill ghost cells surrounding their 3d bricks
|
||||
|
||||
cg_dipole->forward_comm(this,FORWARD_SP);
|
||||
|
||||
// extra per-atom energy/virial communication
|
||||
|
||||
if (evflag_atom) {
|
||||
cg_peratom_dipole->forward_comm(this,FORWARD_SP_PERATOM);
|
||||
}
|
||||
|
||||
// calculate the force on my particles
|
||||
|
||||
fieldforce_ik_spin();
|
||||
|
||||
// extra per-atom energy/virial communication
|
||||
|
||||
if (evflag_atom) fieldforce_peratom_spin();
|
||||
|
||||
// sum global energy across procs and add in volume-dependent term
|
||||
|
||||
const double spscale = mub2mu0 * scale;
|
||||
const double g3 = g_ewald*g_ewald*g_ewald;
|
||||
|
||||
if (eflag_global) {
|
||||
double energy_all;
|
||||
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
energy = energy_all;
|
||||
|
||||
energy *= 0.5*volume;
|
||||
energy -= spsqsum*2.0*g3/3.0/MY_PIS;
|
||||
energy *= spscale;
|
||||
}
|
||||
|
||||
// sum global virial across procs
|
||||
|
||||
if (vflag_global) {
|
||||
double virial_all[6];
|
||||
MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
for (i = 0; i < 6; i++) virial[i] = 0.5*spscale*volume*virial_all[i];
|
||||
}
|
||||
|
||||
// per-atom energy/virial
|
||||
// energy includes self-energy correction
|
||||
|
||||
if (evflag_atom) {
|
||||
double **sp = atom->sp;
|
||||
double spx,spy,spz;
|
||||
int nlocal = atom->nlocal;
|
||||
int ntotal = nlocal;
|
||||
|
||||
if (eflag_atom) {
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
spx = sp[i][0]*sp[i][3];
|
||||
spy = sp[i][1]*sp[i][3];
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
eatom[i] *= 0.5;
|
||||
eatom[i] -= (spx*spx + spy*spy + spz*spz)*2.0*g3/3.0/MY_PIS;
|
||||
eatom[i] *= spscale;
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_atom) {
|
||||
for (i = 0; i < ntotal; i++)
|
||||
for (j = 0; j < 6; j++) vatom[i][j] *= 0.5*spscale;
|
||||
}
|
||||
}
|
||||
|
||||
// 2d slab correction
|
||||
|
||||
if (slabflag == 1) slabcorr();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
create discretized "density" on section of global grid due to my particles
|
||||
density(x,y,z) = charge "density" at grid points of my 3d brick
|
||||
(nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
|
||||
in global grid
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::make_rho_spin()
|
||||
{
|
||||
int l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz;
|
||||
FFT_SCALAR x0,y0,z0;
|
||||
FFT_SCALAR x1,y1,z1;
|
||||
FFT_SCALAR x2,y2,z2;
|
||||
|
||||
// clear 3d density array
|
||||
|
||||
memset(&(densityx_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0,
|
||||
ngrid*sizeof(FFT_SCALAR));
|
||||
memset(&(densityy_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0,
|
||||
ngrid*sizeof(FFT_SCALAR));
|
||||
memset(&(densityz_brick_dipole[nzlo_out][nylo_out][nxlo_out]),0,
|
||||
ngrid*sizeof(FFT_SCALAR));
|
||||
|
||||
// loop over my charges, add their contribution to nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
|
||||
double **sp = atom->sp;
|
||||
double spx,spy,spz;
|
||||
double **x = atom->x;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz);
|
||||
|
||||
spx = sp[i][0]*sp[i][3];
|
||||
spy = sp[i][1]*sp[i][3];
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
z0 = delvolinv * spx;
|
||||
z1 = delvolinv * spy;
|
||||
z2 = delvolinv * spz;
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
y0 = z0*rho1d[2][n];
|
||||
y1 = z1*rho1d[2][n];
|
||||
y2 = z2*rho1d[2][n];
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
x0 = y0*rho1d[1][m];
|
||||
x1 = y1*rho1d[1][m];
|
||||
x2 = y2*rho1d[1][m];
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
densityx_brick_dipole[mz][my][mx] += x0*rho1d[0][l];
|
||||
densityy_brick_dipole[mz][my][mx] += x1*rho1d[0][l];
|
||||
densityz_brick_dipole[mz][my][mx] += x2*rho1d[0][l];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
interpolate from grid to get magnetic field & force on my particles for ik
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::fieldforce_ik_spin()
|
||||
{
|
||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz;
|
||||
FFT_SCALAR x0,y0,z0;
|
||||
FFT_SCALAR ex,ey,ez;
|
||||
FFT_SCALAR vxx,vyy,vzz,vxy,vxz,vyz;
|
||||
|
||||
// loop over my charges, interpolate electric field from nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
|
||||
double **sp = atom->sp;
|
||||
double spx,spy,spz;
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz);
|
||||
|
||||
ex = ey = ez = ZEROF;
|
||||
vxx = vyy = vzz = vxy = vxz = vyz = ZEROF;
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
z0 = rho1d[2][n];
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
y0 = z0*rho1d[1][m];
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
x0 = y0*rho1d[0][l];
|
||||
ex -= x0*ux_brick_dipole[mz][my][mx];
|
||||
ey -= x0*uy_brick_dipole[mz][my][mx];
|
||||
ez -= x0*uz_brick_dipole[mz][my][mx];
|
||||
vxx -= x0*vdxx_brick_dipole[mz][my][mx];
|
||||
vyy -= x0*vdyy_brick_dipole[mz][my][mx];
|
||||
vzz -= x0*vdzz_brick_dipole[mz][my][mx];
|
||||
vxy -= x0*vdxy_brick_dipole[mz][my][mx];
|
||||
vxz -= x0*vdxz_brick_dipole[mz][my][mx];
|
||||
vyz -= x0*vdyz_brick_dipole[mz][my][mx];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// convert M-field to mech. and mag. forces
|
||||
|
||||
const double spfactor = mub2mu0 * scale;
|
||||
spx = sp[i][0]*sp[i][3];
|
||||
spy = sp[i][1]*sp[i][3];
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
f[i][0] += spfactor*(vxx*spx + vxy*spy + vxz*spz);
|
||||
f[i][1] += spfactor*(vxy*spx + vyy*spy + vyz*spz);
|
||||
f[i][2] += spfactor*(vxz*spx + vyz*spy + vzz*spz);
|
||||
|
||||
const double spfactorh = mub2mu0hbinv * scale;
|
||||
fm[i][0] += spfactorh*ex;
|
||||
fm[i][1] += spfactorh*ey;
|
||||
fm[i][2] += spfactorh*ez;
|
||||
|
||||
// create a new vector (in atom_spin style ?) to store long-range fm tables
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
interpolate from grid to get per-atom energy/virial
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::fieldforce_peratom_spin()
|
||||
{
|
||||
int i,l,m,n,nx,ny,nz,mx,my,mz;
|
||||
FFT_SCALAR dx,dy,dz,x0,y0,z0;
|
||||
FFT_SCALAR ux,uy,uz;
|
||||
FFT_SCALAR v0x,v1x,v2x,v3x,v4x,v5x;
|
||||
FFT_SCALAR v0y,v1y,v2y,v3y,v4y,v5y;
|
||||
FFT_SCALAR v0z,v1z,v2z,v3z,v4z,v5z;
|
||||
|
||||
// loop over my charges, interpolate from nearby grid points
|
||||
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
|
||||
// (dx,dy,dz) = distance to "lower left" grid pt
|
||||
// (mx,my,mz) = global coords of moving stencil pt
|
||||
|
||||
double **sp = atom->sp;
|
||||
double spx,spy,spz;
|
||||
double **x = atom->x;
|
||||
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
nx = part2grid[i][0];
|
||||
ny = part2grid[i][1];
|
||||
nz = part2grid[i][2];
|
||||
dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
|
||||
dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
|
||||
dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
|
||||
|
||||
compute_rho1d(dx,dy,dz);
|
||||
|
||||
ux = uy = uz = ZEROF;
|
||||
v0x = v1x = v2x = v3x = v4x = v5x = ZEROF;
|
||||
v0y = v1y = v2y = v3y = v4y = v5y = ZEROF;
|
||||
v0z = v1z = v2z = v3z = v4z = v5z = ZEROF;
|
||||
for (n = nlower; n <= nupper; n++) {
|
||||
mz = n+nz;
|
||||
z0 = rho1d[2][n];
|
||||
for (m = nlower; m <= nupper; m++) {
|
||||
my = m+ny;
|
||||
y0 = z0*rho1d[1][m];
|
||||
for (l = nlower; l <= nupper; l++) {
|
||||
mx = l+nx;
|
||||
x0 = y0*rho1d[0][l];
|
||||
if (eflag_atom) {
|
||||
ux += x0*ux_brick_dipole[mz][my][mx];
|
||||
uy += x0*uy_brick_dipole[mz][my][mx];
|
||||
uz += x0*uz_brick_dipole[mz][my][mx];
|
||||
}
|
||||
if (vflag_atom) {
|
||||
v0x += x0*v0x_brick_dipole[mz][my][mx];
|
||||
v1x += x0*v1x_brick_dipole[mz][my][mx];
|
||||
v2x += x0*v2x_brick_dipole[mz][my][mx];
|
||||
v3x += x0*v3x_brick_dipole[mz][my][mx];
|
||||
v4x += x0*v4x_brick_dipole[mz][my][mx];
|
||||
v5x += x0*v5x_brick_dipole[mz][my][mx];
|
||||
v0y += x0*v0y_brick_dipole[mz][my][mx];
|
||||
v1y += x0*v1y_brick_dipole[mz][my][mx];
|
||||
v2y += x0*v2y_brick_dipole[mz][my][mx];
|
||||
v3y += x0*v3y_brick_dipole[mz][my][mx];
|
||||
v4y += x0*v4y_brick_dipole[mz][my][mx];
|
||||
v5y += x0*v5y_brick_dipole[mz][my][mx];
|
||||
v0z += x0*v0z_brick_dipole[mz][my][mx];
|
||||
v1z += x0*v1z_brick_dipole[mz][my][mx];
|
||||
v2z += x0*v2z_brick_dipole[mz][my][mx];
|
||||
v3z += x0*v3z_brick_dipole[mz][my][mx];
|
||||
v4z += x0*v4z_brick_dipole[mz][my][mx];
|
||||
v5z += x0*v5z_brick_dipole[mz][my][mx];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
spx = sp[i][0]*sp[i][3];
|
||||
spy = sp[i][1]*sp[i][3];
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
if (eflag_atom) eatom[i] += spx*ux + spy*uy + spz*uz;
|
||||
if (vflag_atom) {
|
||||
vatom[i][0] += spx*v0x + spy*v0y + spz*v0z;
|
||||
vatom[i][1] += spx*v1x + spy*v1y + spz*v1z;
|
||||
vatom[i][2] += spx*v2x + spy*v2y + spz*v2z;
|
||||
vatom[i][3] += spx*v3x + spy*v3y + spz*v3z;
|
||||
vatom[i][4] += spx*v4x + spy*v4y + spz*v4z;
|
||||
vatom[i][5] += spx*v5x + spy*v5y + spz*v5z;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Slab-geometry correction term to dampen inter-slab interactions between
|
||||
periodically repeating slabs. Yields good approximation to 2D Ewald if
|
||||
adequate empty space is left between repeating slabs (J. Chem. Phys.
|
||||
111, 3155). Slabs defined here to be parallel to the xy plane. Also
|
||||
extended to non-neutral systems (J. Chem. Phys. 131, 094107).
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::slabcorr()
|
||||
{
|
||||
// compute local contribution to global spin moment
|
||||
|
||||
double **x = atom->x;
|
||||
double zprd = domain->zprd;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double spin = 0.0;
|
||||
double **sp = atom->sp;
|
||||
double spx,spy,spz;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
spin += spz;
|
||||
}
|
||||
|
||||
// sum local contributions to get global spin moment
|
||||
|
||||
double spin_all;
|
||||
MPI_Allreduce(&spin,&spin_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// compute corrections
|
||||
|
||||
const double e_slabcorr = MY_2PI*(spin_all*spin_all/12.0)/volume;
|
||||
const double spscale = mub2mu0 * scale;
|
||||
|
||||
if (eflag_global) energy += spscale * e_slabcorr;
|
||||
|
||||
// per-atom energy
|
||||
|
||||
if (eflag_atom) {
|
||||
double efact = spscale * MY_2PI/volume/12.0;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
eatom[i] += efact * spz * spin_all;
|
||||
}
|
||||
}
|
||||
|
||||
// add on mag. force corrections
|
||||
|
||||
double ffact = spscale * (-4.0*MY_PI/volume);
|
||||
double **fm = atom->fm;
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
fm[i][2] += ffact * spin_all;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute spsum,spsqsum,sp2
|
||||
called initially, when particle count changes, when spins are changed
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PPPMDipoleSpin::spsum_spsq()
|
||||
{
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
spsum = spsqsum = sp2 = 0.0;
|
||||
if (atom->sp_flag) {
|
||||
double **sp = atom->sp;
|
||||
double spx, spy, spz;
|
||||
double spsum_local(0.0), spsqsum_local(0.0);
|
||||
|
||||
// not exactly the good loop: need to add norm of spins
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
spx = sp[i][0]*sp[i][3];
|
||||
spy = sp[i][1]*sp[i][3];
|
||||
spz = sp[i][2]*sp[i][3];
|
||||
spsum_local += spx + spy + spz;
|
||||
spsqsum_local += spx*spx + spy*spy + spz*spz;
|
||||
}
|
||||
|
||||
MPI_Allreduce(&spsum_local,&spsum,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
MPI_Allreduce(&spsqsum_local,&spsqsum,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
sp2 = spsqsum * mub2mu0;
|
||||
}
|
||||
|
||||
if (sp2 == 0 && comm->me == 0)
|
||||
error->all(FLERR,"Using kspace solver PPPMDipoleSpin on system with no spins");
|
||||
}
|
||||
|
|
@ -13,28 +13,23 @@
|
|||
|
||||
#ifdef KSPACE_CLASS
|
||||
|
||||
KSpaceStyle(pppm/spin,PPPMSpin)
|
||||
KSpaceStyle(pppm/dipole/spin,PPPMDipoleSpin)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PPPM_DIPOLE_H
|
||||
#define LMP_PPPM_DIPOLE_H
|
||||
#ifndef LMP_PPPM_DIPOLE_SPIN_H
|
||||
#define LMP_PPPM_DIPOLE_SPIN_H
|
||||
|
||||
#include "pppm.h"
|
||||
#include "pppm_dipole.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PPPMSpin : public PPPM {
|
||||
class PPPMDipoleSpin : public PPPMDipole {
|
||||
public:
|
||||
PPPMSpin(class LAMMPS *, int, char **);
|
||||
virtual ~PPPMSpin();
|
||||
PPPMDipoleSpin(class LAMMPS *, int, char **);
|
||||
virtual ~PPPMDipoleSpin();
|
||||
void init();
|
||||
void setup();
|
||||
void setup_grid();
|
||||
void compute(int, int);
|
||||
int timing_1d(int, double &);
|
||||
int timing_3d(int, double &);
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
double hbar; // reduced Planck's constant
|
||||
|
|
@ -42,55 +37,15 @@ class PPPMSpin : public PPPM {
|
|||
double mu_0; // vacuum permeability
|
||||
double mub2mu0; // prefactor for mech force
|
||||
double mub2mu0hbinv; // prefactor for mag force
|
||||
void set_grid_global();
|
||||
double newton_raphson_f();
|
||||
|
||||
void allocate();
|
||||
void allocate_peratom();
|
||||
void deallocate();
|
||||
void deallocate_peratom();
|
||||
void compute_gf_denom();
|
||||
|
||||
void slabcorr();
|
||||
|
||||
// grid communication
|
||||
|
||||
void pack_forward(int, FFT_SCALAR *, int, int *);
|
||||
void unpack_forward(int, FFT_SCALAR *, int, int *);
|
||||
void pack_reverse(int, FFT_SCALAR *, int, int *);
|
||||
void unpack_reverse(int, FFT_SCALAR *, int, int *);
|
||||
|
||||
// spin
|
||||
|
||||
FFT_SCALAR ***densityx_brick_spin,***densityy_brick_spin,***densityz_brick_spin;
|
||||
FFT_SCALAR ***vdxx_brick_spin,***vdyy_brick_spin,***vdzz_brick_spin;
|
||||
FFT_SCALAR ***vdxy_brick_spin,***vdxz_brick_spin,***vdyz_brick_spin;
|
||||
FFT_SCALAR ***ux_brick_spin,***uy_brick_spin,***uz_brick_spin;
|
||||
FFT_SCALAR ***v0x_brick_spin,***v1x_brick_spin,***v2x_brick_spin;
|
||||
FFT_SCALAR ***v3x_brick_spin,***v4x_brick_spin,***v5x_brick_spin;
|
||||
FFT_SCALAR ***v0y_brick_spin,***v1y_brick_spin,***v2y_brick_spin;
|
||||
FFT_SCALAR ***v3y_brick_spin,***v4y_brick_spin,***v5y_brick_spin;
|
||||
FFT_SCALAR ***v0z_brick_spin,***v1z_brick_spin,***v2z_brick_spin;
|
||||
FFT_SCALAR ***v3z_brick_spin,***v4z_brick_spin,***v5z_brick_spin;
|
||||
FFT_SCALAR *work3,*work4;
|
||||
FFT_SCALAR *densityx_fft_spin,*densityy_fft_spin,*densityz_fft_spin;
|
||||
class GridComm *cg_spin;
|
||||
class GridComm *cg_peratom_spin;
|
||||
int only_spin_flag;
|
||||
double spsum,spsqsum,sp2;
|
||||
double find_gewald_spin(double, double, bigint, double, double);
|
||||
double newton_raphson_f_spin(double, double, bigint, double, double);
|
||||
double derivf_spin(double, double, bigint, double, double);
|
||||
double compute_df_kspace_spin();
|
||||
double compute_qopt_spin();
|
||||
void compute_gf_spin();
|
||||
void make_rho_spin();
|
||||
void brick2fft_spin();
|
||||
void poisson_ik_spin();
|
||||
void poisson_peratom_spin();
|
||||
void fieldforce_ik_spin();
|
||||
void fieldforce_peratom_spin();
|
||||
double final_accuracy_spin();
|
||||
void spsum_spsq();
|
||||
|
||||
};
|
||||
|
|
@ -102,9 +57,9 @@ class PPPMSpin : public PPPM {
|
|||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Cannot (yet) use charges with Kspace style PPPMSpin
|
||||
E: Cannot (yet) use charges with Kspace style PPPMDipoleSpin
|
||||
|
||||
Charge-spin interactions are not yet implemented in PPPMSpin so this
|
||||
Charge-spin interactions are not yet implemented in PPPMDipoleSpin so this
|
||||
feature is not yet supported.
|
||||
|
||||
E: Must redefine kspace_style after changing to triclinic box
|
||||
|
|
@ -123,11 +78,11 @@ E: Cannot (yet) use 'electron' units with spins
|
|||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot yet use triclinic cells with PPPMSpin
|
||||
E: Cannot yet use triclinic cells with PPPMDipoleSpin
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
E: Cannot yet use TIP4P with PPPMSpin
|
||||
E: Cannot yet use TIP4P with PPPMDipoleSpin
|
||||
|
||||
This feature is not yet supported.
|
||||
|
||||
|
|
@ -207,9 +162,9 @@ outside a processor's sub-domain or even the entire simulation box.
|
|||
This indicates bad physics, e.g. due to highly overlapping atoms, too
|
||||
large a timestep, etc.
|
||||
|
||||
E: Using kspace solver PPPMSpin on system with no spins
|
||||
E: Using kspace solver PPPMDipoleSpin on system with no spins
|
||||
|
||||
Must have non-zero spins with PPPMSpin.
|
||||
Must have non-zero spins with PPPMDipoleSpin.
|
||||
|
||||
E: Must use kspace_modify gewald for system with no spins
|
||||
|
||||
|
|
@ -48,7 +48,7 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
|
|||
comm_x_only = 0;
|
||||
comm_f_only = 0;
|
||||
size_forward = 7;
|
||||
size_reverse = 6;
|
||||
size_reverse = 9;
|
||||
size_border = 10;
|
||||
size_velocity = 3;
|
||||
size_data_atom = 9;
|
||||
|
|
@ -58,7 +58,6 @@ AtomVecSpin::AtomVecSpin(LAMMPS *lmp) : AtomVec(lmp)
|
|||
atom->sp_flag = 1;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
grow atom arrays
|
||||
n = 0 grows arrays by a chunk
|
||||
|
|
@ -88,6 +87,7 @@ void AtomVecSpin::grow(int n)
|
|||
|
||||
sp = memory->grow(atom->sp,nmax,4,"atom:sp");
|
||||
fm = memory->grow(atom->fm,nmax*comm->nthreads,3,"atom:fm");
|
||||
fm_long = memory->grow(atom->fm_long,nmax*comm->nthreads,3,"atom:fm_long");
|
||||
|
||||
if (atom->nextra_grow)
|
||||
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
|
||||
|
|
@ -103,7 +103,7 @@ void AtomVecSpin::grow_reset()
|
|||
tag = atom->tag; type = atom->type;
|
||||
mask = atom->mask; image = atom->image;
|
||||
x = atom->x; v = atom->v; f = atom->f;
|
||||
sp = atom->sp; fm = atom->fm;
|
||||
sp = atom->sp; fm = atom->fm; fm_long = atom->fm_long;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -342,6 +342,9 @@ int AtomVecSpin::pack_reverse(int n, int first, double *buf)
|
|||
buf[m++] = fm[i][0];
|
||||
buf[m++] = fm[i][1];
|
||||
buf[m++] = fm[i][2];
|
||||
buf[m++] = fm_long[i][0];
|
||||
buf[m++] = fm_long[i][1];
|
||||
buf[m++] = fm_long[i][2];
|
||||
}
|
||||
|
||||
return m;
|
||||
|
|
@ -361,6 +364,9 @@ void AtomVecSpin::unpack_reverse(int n, int *list, double *buf)
|
|||
fm[j][0] += buf[m++];
|
||||
fm[j][1] += buf[m++];
|
||||
fm[j][2] += buf[m++];
|
||||
fm_long[j][0] += buf[m++];
|
||||
fm_long[j][1] += buf[m++];
|
||||
fm_long[j][2] += buf[m++];
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -939,6 +945,7 @@ bigint AtomVecSpin::memory_usage()
|
|||
|
||||
if (atom->memcheck("sp")) bytes += memory->usage(sp,nmax,4);
|
||||
if (atom->memcheck("fm")) bytes += memory->usage(fm,nmax*comm->nthreads,3);
|
||||
if (atom->memcheck("fm_long")) bytes += memory->usage(fm_long,nmax*comm->nthreads,3);
|
||||
|
||||
return bytes;
|
||||
}
|
||||
|
|
@ -947,6 +954,7 @@ 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);
|
||||
memset(&atom->fm_long[0][0],0,3*nbytes);
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -68,10 +68,12 @@ class AtomVecSpin : public AtomVec {
|
|||
int *type,*mask;
|
||||
imageint *image;
|
||||
double **x,**v,**f; // lattice quantities
|
||||
double **sp,**fm; // spin quantities
|
||||
// sp[i][0-2] direction of the spin i
|
||||
|
||||
// spin quantities
|
||||
double **sp; // sp[i][0-2] direction of the spin i
|
||||
// sp[i][3] atomic magnetic moment of the spin i
|
||||
|
||||
double **fm; // fm[i][0-2] direction of magnetic precession
|
||||
double **fm_long; // storage of long-range spin prec. components
|
||||
};
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -441,8 +441,6 @@ void PairSpinExchange::allocate()
|
|||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
@ -527,4 +525,3 @@ void PairSpinExchange::read_restart_settings(FILE *fp)
|
|||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -13,19 +13,19 @@
|
|||
|
||||
/* ------------------------------------------------------------------------
|
||||
Contributing authors: Julien Tranchida (SNL)
|
||||
Stan Moore (SNL)
|
||||
|
||||
Please cite the related publication:
|
||||
Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018).
|
||||
Massively parallel symplectic algorithm for coupled magnetic spin dynamics
|
||||
and molecular dynamics. arXiv preprint arXiv:1801.10233.
|
||||
------------------------------------------------------------------------- */
|
||||
Aidan Thompson (SNL)
|
||||
|
||||
Please cite the related publication:
|
||||
Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018).
|
||||
Massively parallel symplectic algorithm for coupled magnetic spin dynamics
|
||||
and molecular dynamics. Journal of Computational Physics.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <cstring>
|
||||
|
||||
#include "pair_spin_long.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
|
|
@ -80,10 +80,154 @@ PairSpinLong::~PairSpinLong()
|
|||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(cut_spin_long);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_spin_long_global = force->numeric(FLERR,arg[0]);
|
||||
//cut_spin = force->numeric(FLERR,arg[0]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++) {
|
||||
for (j = i+1; j <= atom->ntypes; j++) {
|
||||
if (setflag[i][j]) {
|
||||
cut_spin_long[i][j] = cut_spin_long_global;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::coeff(int narg, char **arg)
|
||||
{
|
||||
if (!allocated) allocate();
|
||||
|
||||
// check if args correct
|
||||
|
||||
if (strcmp(arg[2],"long") != 0)
|
||||
error->all(FLERR,"Incorrect args in pair_style command");
|
||||
if (narg != 3)
|
||||
error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
double spin_long_cut_one = force->numeric(FLERR,arg[3]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
setflag[i][j] = 1;
|
||||
cut_spin_long[i][j] = spin_long_cut_one;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::init_style()
|
||||
{
|
||||
if (!atom->sp_flag)
|
||||
error->all(FLERR,"Pair spin requires atom/spin style");
|
||||
|
||||
// need a full neighbor list
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
|
||||
// checking if nve/spin is a listed fix
|
||||
|
||||
int ifix = 0;
|
||||
while (ifix < modify->nfix) {
|
||||
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
|
||||
ifix++;
|
||||
}
|
||||
if (ifix == modify->nfix)
|
||||
error->all(FLERR,"pair/spin style requires nve/spin");
|
||||
|
||||
// get the lattice_flag from nve/spin
|
||||
|
||||
for (int i = 0; i < modify->nfix; i++) {
|
||||
if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
|
||||
lockfixnvespin = (FixNVESpin *) modify->fix[i];
|
||||
lattice_flag = lockfixnvespin->lattice_flag;
|
||||
}
|
||||
}
|
||||
|
||||
// insure use of KSpace long-range solver, set g_ewald
|
||||
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Pair style requires a KSpace style");
|
||||
|
||||
g_ewald = force->kspace->g_ewald;
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinLong::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
cut_spin_long[j][i] = cut_spin_long[i][j];
|
||||
|
||||
return cut_spin_long_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
extract the larger cutoff if "cut" or "cut_coul"
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void *PairSpinLong::extract(const char *str, int &dim)
|
||||
{
|
||||
if (strcmp(str,"cut") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_spin_long_global;
|
||||
} else if (strcmp(str,"cut_coul") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_spin_long_global;
|
||||
} else if (strcmp(str,"ewald_order") == 0) {
|
||||
ewald_order = 0;
|
||||
ewald_order |= 1<<1;
|
||||
ewald_order |= 1<<3;
|
||||
dim = 0;
|
||||
return (void *) &ewald_order;
|
||||
} else if (strcmp(str,"ewald_mix") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &mix_flag;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::compute(int eflag, int vflag)
|
||||
|
|
@ -95,6 +239,7 @@ void PairSpinLong::compute(int eflag, int vflag)
|
|||
double evdwl,ecoul;
|
||||
double xi[3],rij[3];
|
||||
double spi[4],spj[4],fi[3],fmi[3];
|
||||
double local_cut2;
|
||||
double pre1,pre2,pre3;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
|
|
@ -156,26 +301,25 @@ void PairSpinLong::compute(int eflag, int vflag)
|
|||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
|
||||
|
||||
if (rsq < local_cut2) {
|
||||
r2inv = 1.0/rsq;
|
||||
rinv = sqrt(r2inv);
|
||||
|
||||
if (rsq < cut_spinsq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
|
||||
bij[0] = erfc * rinv;
|
||||
bij[1] = (bij[0] + pre1*expm2) * r2inv;
|
||||
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
|
||||
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
|
||||
bij[0] = erfc * rinv;
|
||||
bij[1] = (bij[0] + pre1*expm2) * r2inv;
|
||||
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
|
||||
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
|
||||
|
||||
compute_long(i,j,rij,bij,fmi,spi,spj);
|
||||
compute_long_mech(i,j,rij,bij,fmi,spi,spj);
|
||||
|
||||
}
|
||||
compute_long(i,j,rij,bij,fmi,spi,spj);
|
||||
compute_long_mech(i,j,rij,bij,fmi,spi,spj);
|
||||
}
|
||||
|
||||
// force accumulation
|
||||
|
|
@ -194,7 +338,7 @@ void PairSpinLong::compute(int eflag, int vflag)
|
|||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_spinsq) {
|
||||
if (rsq <= local_cut2) {
|
||||
evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] +
|
||||
spi[2]*fmi[2];
|
||||
evdwl *= hbar;
|
||||
|
|
@ -219,6 +363,7 @@ void PairSpinLong::compute_single_pair(int ii, double fmi[3])
|
|||
double r,rinv,r2inv,rsq;
|
||||
double grij,expm2,t,erfc;
|
||||
double bij[4],xi[3],rij[3],spi[4],spj[4];
|
||||
double local_cut2;
|
||||
double pre1,pre2,pre3;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
|
|
@ -267,25 +412,24 @@ void PairSpinLong::compute_single_pair(int ii, double fmi[3])
|
|||
rij[2] = x[j][2] - xi[2];
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
|
||||
|
||||
if (rsq < local_cut2) {
|
||||
r2inv = 1.0/rsq;
|
||||
rinv = sqrt(r2inv);
|
||||
|
||||
if (rsq < cut_spinsq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
|
||||
bij[0] = erfc * rinv;
|
||||
bij[1] = (bij[0] + pre1*expm2) * r2inv;
|
||||
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
|
||||
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
|
||||
bij[0] = erfc * rinv;
|
||||
bij[1] = (bij[0] + pre1*expm2) * r2inv;
|
||||
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
|
||||
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
|
||||
|
||||
compute_long(i,j,rij,bij,fmi,spi,spj);
|
||||
|
||||
}
|
||||
compute_long(i,j,rij,bij,fmi,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -361,111 +505,7 @@ void PairSpinLong::allocate()
|
|||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_spin = force->numeric(FLERR,arg[0]);
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 4 || narg > 5)
|
||||
error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
// check if args correct
|
||||
|
||||
if (strcmp(arg[2],"long") != 0)
|
||||
error->all(FLERR,"Incorrect args in pair_style command");
|
||||
if (narg != 3)
|
||||
error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinLong::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
double cut = cut_spin;
|
||||
return cut;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinLong::init_style()
|
||||
{
|
||||
if (!atom->sp_flag)
|
||||
error->all(FLERR,"Pair spin requires atom/spin style");
|
||||
|
||||
// need a full neighbor list
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
|
||||
// checking if nve/spin is a listed fix
|
||||
|
||||
int ifix = 0;
|
||||
while (ifix < modify->nfix) {
|
||||
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
|
||||
ifix++;
|
||||
}
|
||||
if (ifix == modify->nfix)
|
||||
error->all(FLERR,"pair/spin style requires nve/spin");
|
||||
|
||||
// get the lattice_flag from nve/spin
|
||||
|
||||
for (int i = 0; i < modify->nfix; i++) {
|
||||
if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
|
||||
lockfixnvespin = (FixNVESpin *) modify->fix[i];
|
||||
lattice_flag = lockfixnvespin->lattice_flag;
|
||||
}
|
||||
}
|
||||
|
||||
// insure use of KSpace long-range solver, set g_ewald
|
||||
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Pair style requires a KSpace style");
|
||||
|
||||
g_ewald = force->kspace->g_ewald;
|
||||
|
||||
cut_spinsq = cut_spin * cut_spin;
|
||||
memory->create(cut_spin_long,n+1,n+1,"pair:cut_spin_long");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -477,10 +517,14 @@ void PairSpinLong::write_restart(FILE *fp)
|
|||
write_restart_settings(fp);
|
||||
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (i = 1; i <= atom->ntypes; i++) {
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j]) {
|
||||
fwrite(&cut_spin_long[i][j],sizeof(int),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -495,11 +539,18 @@ void PairSpinLong::read_restart(FILE *fp)
|
|||
|
||||
int i,j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (i = 1; i <= atom->ntypes; i++) {
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) {
|
||||
fread(&cut_spin_long[i][j],sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_spin_long[i][j],1,MPI_INT,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -508,7 +559,7 @@ void PairSpinLong::read_restart(FILE *fp)
|
|||
|
||||
void PairSpinLong::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_spin,sizeof(double),1,fp);
|
||||
fwrite(&cut_spin_long_global,sizeof(double),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
|
|
@ -519,32 +570,9 @@ void PairSpinLong::write_restart_settings(FILE *fp)
|
|||
void PairSpinLong::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_spin,sizeof(double),1,fp);
|
||||
fread(&cut_spin_long_global,sizeof(double),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_spin,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_spin_long_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairSpinLong::extract(const char *str, int &dim)
|
||||
{
|
||||
if (strcmp(str,"cut") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_spin;
|
||||
} else if (strcmp(str,"cut_coul") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_spin;
|
||||
} else if (strcmp(str,"ewald_order") == 0) {
|
||||
ewald_order = 0;
|
||||
ewald_order |= 1<<1;
|
||||
ewald_order |= 1<<3;
|
||||
dim = 0;
|
||||
return (void *) &ewald_order;
|
||||
} else if (strcmp(str,"ewald_mix") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &mix_flag;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -49,6 +49,8 @@ class PairSpinLong : public PairSpin {
|
|||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
double cut_spin_long_global; // global long cutoff distance
|
||||
|
||||
protected:
|
||||
double hbar; // reduced Planck's constant
|
||||
|
|
@ -56,7 +58,8 @@ class PairSpinLong : public PairSpin {
|
|||
double mu_0; // vacuum permeability
|
||||
double mub2mu0; // prefactor for mech force
|
||||
double mub2mu0hbinv; // prefactor for mag force
|
||||
double cut_spin, cut_spinsq;
|
||||
|
||||
double **cut_spin_long; // cutoff distance long
|
||||
|
||||
double g_ewald;
|
||||
int ewald_order;
|
||||
|
|
|
|||
Loading…
Reference in New Issue